title: Two-day Music chart predictions and the serial correlation effect on the prediction

text: I have about 64000 music Charts ranked by their usage frequency. I want to have a future two-day prediction frequency and eventually its rank for each music chart using its past 21 days usage frequencies but the frequencies are obviously correlated but because of scale, I dont want to use ARIMA related models. I am not sure what would be the effect of correlated residuals on my predictions. Constraints:   1-Even though the frequencies are correlated (time series data), I don't want to use ARIMA related models as for a production scale, they aren't stable and break a lot (Matrix singularities,..) 2- Because of the scale, it is very much preferable to just use the 21 data points for each record, independent of other records (no Mixed effect model) 3- The model must outperform the current method ( Using today's record rank as the predicted rank for two days later (base model). The base model for the top rank records is fairly good and hard to beat as the daily ranks don't fluctuate much. This is what I have been doing. 1- Use this model: freq = b0 + b1 Day(-2) i.e. I use days 1:21 to predict the frequencies of days 3:23 (day1 for freq on day3,...). A weighted least square approach   lm(freq ~ poly(day,2),data=df,weights=wgh) is used to counter the non-constant variance. 2- Since this was not good enough, I used a weighted average of predicted frequency and the base frequency. For example if ( 1198,1234) are the predicted frequencies of day 23 and the actual frequency on day 21 (two days earlier(the base frequency)) respectively.   My final predicted freq for day 23 will be   ** w freq_pred[23] + (1-w) freq[21] = w * 1198 + (1-w) 1234** for some w. That is the only way I can beat the base prediction for top records.   How can I improve my model? How will the unaccounted serial correlation will affect my results?   Are there other suggestions? Is there a big problem with this approach?


title: Sample size and cross-validation methods for Cox regression predictive models

text: I have a question I would like to pose to the community. I have recently been asked to provide statistical analysis for a **tumor marker prognostic study**. I have primarily used these two references to guide my analysis:   1. McShane LM, et al. Reporting recommendations for tumor marker prognostic studies (REMARK). J Natl Cancer Inst. 2005 Aug 17; 97(16):1180-4.   2. Simon RM, et al. Using cross-validation to evaluate predictive accuracy of survival risk classifiers based on high-dimensional data. Brief Bioinform. 2011 May; 12(3):203-14. Epub 2011 Feb 15. I have summarized the study and my analyses below. I would appreciate any comments, suggestions, or criticisms. ## Study background: Some patients with cancer X experience early relapse after treatment. The clinical prognostic score currently used by doctors does not do a good job of predicting clinical outcome in these patients. It would therefore be useful to identify biological prognostic markers that add value above and beyond this standard score. The goal of this study is to discover such a biomarker. ## Study methods: ### Pre-selection of candidate biomarkers Twelve biomarkers associated with cancer X were identified in a previous study. We attempted to validate the association between these 12 candidates and cancer X in an independent sample of patients/tumors, described below. ### Univariate validation of pre-selected candidate biomarkers Levels of these biomarkers were measured in a set 220 patients/tumors. [Note: I have masked the data and made them available for public download as a *.csv file. The file has the following columns: “ID”, a unique identifier for each patient; “PS”, the prognostic score for each patient, with 1 indicating a good prognosis and 2 indicating a bad prognosis; “m1” to “m12”, levels of each tumor marker; “time”, in months; and “event”, where 0 indicates that the observation is censured and 1 indicates that treatment failure occurred.] Univariable Cox regression models with time to death as the dependent variable were built for each of the 12 biomarkers (n = 220 observations, number of events = 91).                   Risk  LCI  UCI pValue     1   0.93 0.86 1.02 0.1088     2   0.93 0.88 0.99 0.0215     3   0.99 0.92 1.05 0.6528     4   0.93 0.87 1.00 0.0468     5   0.93 0.88 0.98 0.0055     6   0.97 0.92 1.01 0.1202     7   0.91 0.83 0.99 0.0297     8   0.98 0.90 1.07 0.6972     9   0.99 0.92 1.06 0.7841     10  1.01 0.91 1.11 0.9149     11  0.96 0.87 1.05 0.3837     12  0.90 0.83 0.97 0.0047      Using a threshold p value of 0.05/12 = 0.004, none of the results were significant. ### Multivariable analyses It was decided to fit a model to the data by inputting all 12 biomarkers at once into a stepwise Cox regression algorithm using ten-fold cross-validation. After building ten models on the ten different training sets, time-dependent ROC curves were built to allow selection of optimal cutoff points to identify two groups of patients, “high” and “low” risk. Cut points that minimized “1 - TP + FP” were selected. These ten models were then asked to make predictions about the corresponding patients in the validation groups. These patients were then classified into “high” and “low” risk groups and plotted on a single, cross-validated Kaplan Meier curve. ## Conclusions The confidence intervals of the high and low risk curves significantly overlapped, suggesting that the identified biomarkers were not useful prognostic markers. Our study therefore has not identified any significant univariate or multivariate associations between these markers and patient prognosis. ## Questions for the community Have I gone about analyzing my data in the correct manner? If you had been the statistician on this study, would you have done something differently? Prior to performing the validation analyses, sample size and power calculations were not performed to determine the number of samples to include and the detectable effect size. I would like to perform these analyses now to guide future studies. Can someone tell me how to do this? What I am really interested in is whether these biomarkers provide predictive information above and beyond the clinical prognostic score. From what I understand, this would entail making three different models: (1) a model with clinical covariates only, (2) a biomarker model with biomarker covariates only, and (3) a biomarker/clinical model based on both types of covariates. So far I have made models 1 (not shown above; it was unable to differentiate between high and low risk patients in our sample either) and 2 (shown above). Because 1 and 2 were not significant, I didn’t make model 3. Should I do this any way? Any additional comments about analytical concerns would be greatly appreciated! Please feel free to download the masked data and have a look yourself.


title: Modeling what should be a logistic regression but has no negative responses

text: I have a data set of reported food-borne illnesses and we're trying to determine what environmental conditions during food cultivation led to high bacterial counts in the food, and thus caused the illnesses. Unfortunately, I only have data of foods that caused confirmed illnesses. I requested that we go back and "randomly" sample from food tags that did not cause a reported illness but am not allowed to do so for various reasons. Even that would have had problems (because just because an illness is not reported doesn't mean it didn't occur), but at least this would have given me some negative observations. I was originally planning to model these data using a logistic regression but I am stuck at what to do now. Without negative observations, I can only really provide univariate descriptive statistics, right? I'm hoping that someone else has had this problem and perhaps there's some model I haven't heard of before that can handle this. Thank you.


title: Comparing entries in correlation matrices

text: I am trying to compare entries in a 10x10 Spearman correlation matrix to determine whether one entry is closer to another entry than it is to a third entry. The caveat is that I actually have 1000 of these correlation matrices, not one, and the matrices are not statistically independent. To give some background: I have 1000 molecules and 10 ways of comparing them.   Each row and column in the matrices correspond to one of the ways of comparing the molecules and row-column pair is the correlation between the comparison methods. Each matrix corresponds to one molecule being compared to all the other molecules, so each cell in that matrix corresponds to comparing two methods using that molecule as a reference against the 1000 molecules. If I had one matrix, I would have converted the entries using Fisher's Z transformation the way it was done here: Comparing correlation coefficients With multiple matrices, should I a) Take a mean of each of the cells (spearman's coeficients) across the 1000 matrices and then perform the transformation on the matrix of the means? or b) First perform the Fisher's Z transformation on all cells in every matrix, then take a mean of the cells over 1000 matrices and use the new matrix of means of Fisher's Z transformations for comparisons? c) Something else? Edit1: So far I have this: ![Comparing two molecular similarity methods](http://i.stack.imgur.com/pVhTG.png) The bottom left histogram corresponds to a histogram of Spearman correlations between two comparison methods, where each molecule is picked as a reference. Hence, there are 1000 points in the histogram. If I picked just one molecule as a reference to produce one matrix, I wouldn't know which one it is best to pick and why I should ignore others. Edit2: Can I just treat them as populations of tau and rho coefficients and say that the mean of one population is significantly larger or smaller than the mean of another population? In other words, have 10 populations of 1000 members each and try to determine whether any are significantly larger or smaller than others?


title: What is percentage of variance in PCA?

text: I am trying to understand Principal Component Analysis (PCA). I found a webpage on PCA that introduces it and the concept of the percentage of variance. However, I am very confused about what "Percentage of Variance" (POV) means. Here are the questions I have:   1. Is there a formal definition or mathematics formula to define "percentage of variance"?    2. Do we only use POV for calculating PCA? Is POV used somewhere else?    3. Is POV the same as "Proportion of Variance Explained"? (There are many similar terms online which really confuse me.)


title: Likelihood for Poisson data

text: In my book, it says: > Independent random variables $X_1, X_2, \dots, X_n$ are modeled by a Poisson > distribution with mean $\lambda > 0$. The likelihood for $\lambda$ based on > data $\mathbf{x}=(x_1,x_2,\dots\,x_n)^T$ is $$L(\lambda)=\prod_{i=1}^n > \frac{\lambda^{x_i}e^{-\lambda}}{x_i!}$$ $$=k \lambda^{n > \bar{x}}e^{-n\lambda}$$ where $k$ is a constant. At first I was a bit confused by this, but I think I get it now. The $e^{-n\lambda}$ bit is easy. For the $\lambda^{n \bar{x}}$ part I think they had $\lambda^{x_1+x_2+\dots,x_n}=\lambda^{\sum{x_i}}$ so $\bar{x}=\frac{\sum{x_i}}{n}$ and $\sum{x_i}=n \bar{x}$ and then $k= \frac{1}{\prod_{i=1}^nx_i!}$ which doesn't depend on $\mathbf{\lambda}$. Have I missed anything ?


title: What's the forward stagewise regression algorithm?

text: Maybe it's just that I'm tired, but I'm having trouble trying to understand the Forward Stagewise Regression algorithm. From "Elements of Statistical Learning" page 60: > Forward-stagewise regression (FS) is even more constrained than forward- > stepwise regression. It starts like forward-stepwise regression, with an > intercept equal to [the mean of] y , and centered predictors with coeﬃcients > initially all 0. > > At each step the algorithm identiﬁes the variable most correlated with the > current residual. It then computes the simple linear regression coeﬃcient of > the residual on this chosen variable, and then adds it to the current co- > eﬃcient for that variable. This is continued till none of the variables have > correlation with the residuals—i.e. the least-squares ﬁt when N > p. So, is this the algorithm?:               b[1]=mean(y)     b[2..n]=0     r=(y-X*b)     index, maxCorr = max(transpose(r)*X)     while(abs(maxCorr) > someThreshold)       b[index]=b[index]+regress(r,X[1..n][index])       r=(y-X*b)       index, maxCorr = max(transpose(r)*X)      Where b is a column-vector of the coefficients, X is a matrix of inputs, and y is a column-vector of outputs. I.e. y=X*b+error. Asking because this algorithm gives me only a few non-zero coefficients on the dataset I'm testing it on (with threshold=.0001), and the prediction accuracy isn't very good at all.


title: Fuzzy K-means - Cluster Sizes

text: I'm trying to do fuzzy k-means clustering on a dataset using the cmeans function (R) . The problem Im facing is that the sizes of clusters are not as I would like them to be. This is done by calculating the cluster to which the observations are "closest".               cl$size      [1]   108    31   192    51   722 18460    67  1584   419 17270      Here we see that for 10 clusters we have two huge clusters and a lot of very small ones. Does this imply that two clusters are optimal in any way? If I do regular K-means 10 segments look very well, with good sizes and their intepretation makes a lot of sense but I would like to try fuzzy correctly. I just started exploring this fuzzy clustering so any help and pointers are overly welcome.


title: Choose among poststratification/calibration weights

text: I have data from a survey, stratified by country with a roughly equal number of participants in each country. Also, I have calibration totals, three dimensions (say, age, sex and work status) plus the country dimensions. A few of the category combinations of the calibration data have not been covered by the survey. I have merged categories using two ad-hoc strategies to allow for poststratification. For the first strategy, I combined "neighboring" categories (e.g., age 50 and age 55 of male unemployed); the second strategy consisted of collapsing entire categories (e.g., treating the unemployed and the "working at home" as one stratum). In addition, I calibrated against marginal totals using GREG estimation. My question is: Is there some general guidance which weights would be preferable, and why? The client provided us with indicators for which the weights should be checked. The point estimates vary slightly after the reweighting, but the variance by country _increases_ for some countries and decreases for others. Can I decide from the data which weights to prefer?


title: What is the effect of merging categories on logistic regression estimates?

text: I have a logistic regression with a number of predictors variables including a factor with, say, 5 categories. The estimates for that factor compare categories 2-5 in turn to the reference category zero. If I merge 2 of the non-reference categories (say merge 2 and 3) then the estimates for the other change slightly. I think this is because estimates will be correlated in a non-orthognal design but I'd like confirmation or a more thorough explanation please. ### Edit to be honest my colleague has such a problem with conditional logistic regression but the same issue occurs with ordinary logistic regression and I suspect in regressoin generally. for myself I'd like to know for general knowledge what the issue is (as well as to reassure my colleague I hope it is OK to post code (I hope it doesn't get too scrambled - I indented the code) but this illustrates the occurence with the generally available (e.g. in `R` package `MASS` from the Venables and Ripley book) dataset `birthwt` where binary outcome "low birthweight" is modelled by logistic regression - I have a factor `ftv` (which I make groups 0-4 by merging then I create s copy `fvt2` where groups 1 and 2 are merged to be "2". The estimates in both analyses compare group zero to the other groups and you can see below the comparison of group 4 to 0 changes (slightly - I've put them in bold). Why? the data for those 2 groups in unaltered but I can see there is an issue of correlated estimates here. Many thanks for any advice.               library(MASS)     data(birthwt)     birthwt$ftv<-ifelse(birthwt$ftv>3,4,birthwt$ftv)     # change 4 or more to 4     birthwt$ftv2<-birthwt$ftv     # create duplicate     birthwt$ftv2<-ifelse(birthwt$ftv==1,2,birthwt$ftv)     # then merge categories in duplicate     table(birthwt$ftv,birthwt$ftv2)      # 1 and 2 for ftv become 2 in ftv2                   0   2   3   4     0 100   0   0   0     1   0  47   0   0     2   0  30   0   0     3   0   0   7   0     4   0   0   0   5               low_log<-glm(low ~  ht + ui + lwt +as.factor(ftv), data= birthwt,family = binomial)     coef(summary(low_log))                         Estimate  Std. Error    z value    Pr(>|z|)     (Intercept)      1.15892582 0.852443200  1.3595344 0.173977321     ht               1.87181914 0.710991788  2.6326874 0.008471227      ui               0.90472684 0.442643299  2.0439185 0.040961597     lwt             -0.01627838 0.006638716 -2.4520372 0.014204999     as.factor(ftv)1 -0.55718506 0.419152552 -1.3293133 0.183744632     as.factor(ftv)2 -0.39179846 0.501082892 -0.7819035 0.434271296     as.factor(ftv)3  0.76139653 0.837278585  0.9093706 0.363154537     as.factor(ftv)4 **-0.17724595** 1.183684076 -0.1497409 0.880969020          low_log<-glm(low ~  ht + ui + lwt +as.factor(ftv2), data= birthwt,family = binomial)     coef(summary(low_log))                          Estimate  Std. Error    z value    Pr(>|z|)     (Intercept)       1.14052105 0.850496353  1.3410064 0.179918375     ht                1.84783639 0.704472332  2.6230078 0.008715729     ui                0.91241494 0.441875434  2.0648691 0.038935385     lwt              -0.01612754 0.006620628 -2.4359530 0.014852614     as.factor(ftv2)2 -0.49476082 0.354944684 -1.3939097 0.163344887     as.factor(ftv2)3  0.76274256 0.836785227  0.9115153 0.362023910     as.factor(ftv2)4 **-0.17931225** 1.183211033 -0.1515471 0.879544125      Actually if you remove the other variables - `ht`, `ui` and `lwt`, then the estimates/contrasts for groups 4 and 3 vs 0 are unaffected! The estimates are always the difference in log odds between the 2 groups I am uneasy this is well understood (at least by others).


title: Junction tree theorem on cycles

text: In page 31 of David Sontag's 2010 MIT PhD thesis 1 ( _Approximate Inference in Graphical Models using LP Relaxations_ ) we read that: > The junction tree theorem guarantees that every joint distribution on the > variables in the cycle can be equivalently factored into triplet > distributions along the triangulation that are consistent with one another, > and vice versa. I partly know the junction tree theorem, but I don't see why "triplets". Even if we group any two consecutive vertices, it should give us consistent answer. Any comments? Am I misunderstanding the JT theorem?


title: How is this sampling called? (getting a "representative" sample set for teaching)

text: I know of random sampling and stratified sampling, but what I am not sure what is the name of the type of sampling I need: I wish to work-out an example with my students of simple correlation/regression. I found the following nice dataset: "When do Babies Start to Crawl", but it has too many observations (12). And what I wish is something even smaller that we can do calculations on (let's say 6). I would like to sample these points, but I would rather do it in such a way so that the result of the analysis would remain as similar as possible to what I would get from the larger dataset. So of course the SE would change, but I would like the correlation to be as close as possible, and if possible to have the range stay similar (I'm ignoring the outlier issue). Does this type of situation have a known name which I can use when searching for good solutions? Thanks.


title: Community finding algorithms in large heterogeneous networks

text: Consider a network that consists of vertices with various meanings. For example: stack overflow users, keywords and user location when asking/answering a question. In this network, when a user asks a question its ID node is linked to several keyword vertices and to a vertex representing users' geographical location. I would like to ask a general and vague question: are there any meaningful communities in this network? One approach would be to link (lets say) users using other vertices and then analyzing the resulting graph of users. However, are there approaches that retain the heterogeneous types of information in the graph? Are there metric measurements that take into consideration various types of nodes?


title: Treating ordinal variables as continuous for regression problems

text: In the social sciences I have encountered that it is common to treat ordinal variables as continuous, for example variables originating from rating or Likert scales (strongly disagree, disagree, agree, strongly agree). This topic has been discussed for example in this post from 2010: Under what conditions should Likert scales be used as ordinal or interval data? I am looking for a more formal comparison/evaluation especially in the context of regression modeling. Rhemtulla et al. (2012) examine the performance of treating ordinal variables as continuous and make recommendations for structural equation models (SEM). I am not very familiar with SEMs, so I'm not sure if their results would also apply to regression problems. Does anyone know about similar studies/literature in the context of regression? **Edit:** Just to answer the question below: I'm mainly interested in the case where the outcome variable is ordinal (with possibly an ordinal covariate).


title: Topic models evaluation in Gensim

text: I've been experimenting with LDA topic modelling using Gensim. I couldn't seem to find any topic model evaluation facility in Gensim, which could report on the perplexity of a topic model on held-out evaluation texts thus facilitates subsequent fine tuning of LDA parameters (e.g. number of topics). It would be greatly appreciated if anyone could shed some light on how I can perform topic model evaluation in Gensim. This question had also been posted on Stackoverflow.


title: How to present a large dataset of polynomial multivariable regression in a useful format

text: I've produced a large dataset of results from 3 variable inputs. eg: x1, x2, x3 -> y y(x1,x2,x3) I've gone over the problem logically and it was clear the problem was non linear. To this end I've used a polynomial collection based on x1, x2, and x3 (eg. x1^2 and x1*x2) to perform linear regression. I've got a great set of results. my p-values for each variable (treating x1 and x1^2 as separate variables) are all around 0. My R^2 value is 0.995. I'm pretty happy with this outcome. **However, how can I present this data in a useful manner?** I can't work out how to graph the data in any useful way, and I have about 400 different combinations of x1, x2 and x3, so I can't really present the whole table. The equation and statistics are OK but I feel that I should have something which shows the data in a little more detail. Thanks! Edit: I left some ambiguity in '400 combinations' The 400 combinations are 400 different combinations used for x1, x2, x3 as inputs. Not 400 different combination of x1x2, x1x3 x1^2x3 etc.


title: Creating a bound for response variable using softImpute?

text: I'm working on the Netflix challenge in R and was I'm curious if there is a way to create a bound for the possible responses that the softImpute algorithm can predict. It doesn't look like there is a built in functionality for that but perhaps somebody here could give me a suggestion. The reason I need a boundary is because I'm trying to incorporate the extra data (movie type, age, etc.) and those variables are all 1's and 0's once I model.matrix the dataset. Once I get the dataset into 'tidy' form (observations on rows, variables on columns) it becomes different from the user x movie form it was originally in. Not sure if that's a smart way to look at the data but its massive (millions of rows) so I'm having some trouble effectively playing around with it. I've already played around with softImpute on purely the users x movies data frame and don't feel like I can improve much without incorporating the extra data.


title: Interpretation C value in linear SVM

text: My C value is very low (close to 0). **Does this mean that my feature (dimensions) have no real separative (and thus predictive) value?** (As the SVM basically chooses to ignore the training data because C is almost zero and thus doesn't penalize misclassifications)? PS: if relevant to your answer:   * Several Cs (2^-10 to 20^10) were tested on a validation dataset and the model with the best accuracy was chosen to execute on the testing dataset.   * the values of my feature (dimensions) are small (between -1 and +1) because of the nature of the data.   * Their weight ( _w values_ ) is small too.


title: Best way to classify a set through a single feature

text: I need to classify a single dataset through a numeric value. I added below a simple dataset to explain what I need. **Restriction:** Category has two values: `1` or `2` The question is **"What is the best`T` score to classify new records through `T` score".** Sample data:               -------------------------------     App     Category       T score     -------------------------------     X       1              180     Y       1              75     Z       2              220     A       2              120     B       1              180


title: Determining market share from multiple choice questions on a survey

text: Say, for example, that I want to determine the market share or relative popularity of coffee houses within a certain population through a survey. What is the best way to write a question that accurately measures this? Some issues that I am thinking of: I don't want to ask a single choice question ("Which coffee house do you go to?"), because the coffee houses are not mutually exclusive. I may enjoy more than one equally often. If I ask a multiple choice question, then I can't really get a real "market share," since the sum of the proportions of people who go to each coffee house sum to over 100%. I can only say that "x% of people in this population go to this coffee house." Is it possible to ask a series of questions ("Which coffee house do you prefer the most?" "Which coffee house do you prefer the second most?", or "Rank the following coffee houses")? Can I ask a multiple choice question then rebase the proportion to the total number of responses? For example, if I have 100 respondents, but they selected 200 coffee houses (because each respondent said they went to two coffee houses, maybe), can I calculate a frequency table based on 200, the number of selections, instead of the number of respondents?


title: How to check that a sample suits multi-dimensional uniform distribution?

text: I have a 3-dimensional sample $(X_k,Y_k,Z_k), k=1, \ldots, N$ which I suspect to be uniform on some parallelepiped in $R^3$ (i.e. a set of the form [a;b]X[c;d]X[e;f], where numbers a,b,c,d,e,f are unknown).   1. How should I estimate numbers a, b, c, d, e, f? Obviously I can try MLE, but then my estimates are biased. Does unbiased estimates exist?   2. How can I check that my sample is indeed uniform?


title: Updating the lasso fit with new observations

text: I am fitting an L1-regularized linear regression to a very large dataset (with n>>p.) The variables are known in advance, but the observations arrive in small chunks. I would like to maintain the lasso fit after each chunk. I can obviously re-fit the entire model after seeing each new set of observations. This, however, would be pretty inefficient given that there is a lot of data. The amount of _new_ data that arrives at each step is very small, and the fit is unlikely to change much between steps. Is there anything I can do to reduce the overall computational burden? I was looking at the LARS algorithm of Efron et al., but would be happy to consider any other fitting method if it can be made to "warm-start" in the way described above. Notes:   1. I am mainly looking for an algorithm, but pointers to existing software packages that can do this may also prove insightful.   2. In addition to the current lasso trajectories, the algorithm is of course welcome to keep other state. > Bradley Efron, Trevor Hastie, Iain Johnstone and Robert Tibshirani, Least > Angle Regression, _Annals of Statistics_ (with discussion) (2004) 32(2), 407 > --499.


title: How to identify the formula for calculating Cohen's d?

text: If I am just asked to calculate Cohen's d what do I do? Do I take the sample mean from the population mean and then divide that by the sample std deviation?


title: When reporting a CFA, is it desirable to assess univariate normality in addition to multivariate normality?

text: I am using SPSS AMOS to do a factor analysis and it produces statistics related to both univariate and multivariate normality. p35 of the AMOS Users Guide states that "[T]he observed variables must meet some distributional requirements. If the observed variables have a multivariate normal distribution, that will suffice". I have also been told in person by a friend that univariate normality should also be reported. In addition to that, this paper states on p19 that "[A]ssessment of univariate and multivariate normality should be provided." Is it the case that if I have a problem with univariate normality then that will also be reflected in the multivariate normality statistic? If so, would that be a good argument for me just reporting multivariate normality?


title: Q: plot glmm fixed and random effects (glmer in package lme4) using ggplot2

text: I am trying to visualize the results from a glmm that I ran with the lme4 package.               M8<-glmer(abundance~ Mom+ Mom*settlment2 + (1|Pop) + (1|obs), family=poisson, data=glm)      My response variable are count data, therefore I have chosen a poisson distribution. Also, my data are overdispersed, so I have added in observation as a random effect to help account for this. However, my data still seem to be violating independence, which is why I am using a generalised linear mixed model. However, the built-in R predict function breaks down for lmer type objects, so I using the predictSE.mer command from the package (AICcmodavg).               plot(abundance, settlment2, xlab="settlment", ylab="abundance", data=glm) #plot my raw data     MyData<-data.frame(settlment2=seq(from=0, to = 350, by=15.3), glm$Mom) #make a new data frame          G<-predictSE.mer(M8, newdata=MyData, type="link", se.fit=TRUE, level=0, print.matrix=FALSE) #now to calculate fitted means         F<-exp(G$fit)     lines(MyData$settlment2, F, lty=1)      ![plot with glmer model](http://i.stack.imgur.com/w928g.png) This seems to work somewhat well, but I am having trouble producing data.frames correctly as I increase the number of fixed factors and interactions. So, I guess what I am looking for is a way to do this in ggplot2. It works well with lme objects from nlme but not glmer objects. For example using lme:               m11<-lme(abundance~settlment2*Pop + Pop + settlment2, random=~1|Mom, data=glm,  method="ML", na.action=na.exclude)          g0 <-ggplot(glm,aes(x=settlment2,y=abundance,colour=Pop))+geom_point(aes(group=Pop),alpha=0.3)+     theme_bw()          g0 + geom_line(data=transform(glm,abundance=predict(m11)), lty=1)      ![plot from lme model](http://i.stack.imgur.com/f7hBK.png) Essentially, what I would like to do is use ggplot2 so that I can break up the results like the second graph, into my two populations (Pop) but using the code above for glmer and not with lme. I have been reading Mixed Effects Models and Extension in Ecology in R (Zuur et al.) trying to reproduce their deer data (pg. 326) or (pg. 329) but instead of probabilities on the Y-axis, I would like just predicted values. I have looked into the function plotLMER.fnc in the package LMERConvienienceFunctions as well, but I would really like to know if this is possible to do in ggplot2 as I am much more familiar with that package. Thanks in advance for the advice!


title: Obtaining the Standard Error of Weighted Data in SPSS

text: I'm trying to find confidence intervals for the means of various variables in a database using SPSS, and I've run into a spot of trouble. The data is weighted, because each of the people who was surveyed represents a different portion of the overall population. For example, one young man in our sample might represent 28000 young men in the general population. The problem is that SPSS seems to think that the young man's database entries each represent 28000 measurements when they actually just represent one, and this makes SPSS think we have much more data than we actually do. As a result SPSS is giving very very low standard error estimates and very very narrow confidence intervals. I've tried fixing this by dividing every weight value by the mean weight. This gives plausible figures and an average weight of 1, but I'm not sure the resulting numbers are actually correct. Is my approach sound? If not, what should I try? I've been using the Explore command to find mean and standard error (among other things), in case it matters.


title: Sampling from a long tailed distribution

text: I have about 500K records from a dataset that consists of counts (so every record is a count of something, like the number of attempts by an IP address to connect to a website). I know, a priori, that each record belongs to one of two groups (for example malicious IPs and legitimate IPs), but I don't know which one without manually inspecting the record and determining the grouping. I also know that the distribution of counts is very long tailed (both looking at my 500K sample and a priori knowledge.) In order to get a better idea of how the grouping changes as the counts value increases, I'd like to take a reasonably small sample of records (say, no more than 1000) to examine manually. What is the right way to sample this data, considering that a simple random sample would miss the tail?


title: How can I generate ensembles using Copula?

text: I just got buried in a mountain of math and papers, and I am so confused. So basically, I have data (Probabilistic QPE) from which I have extracted 99 quantiles (1%, 2%, ... 100%) to sample the CDF of QPE for each point. I would like to use Copula (any) to generate ensembles. And I am not sure how to proceed. I have read the joy of Copulas and still no headway. I have seen some examples, where they fit the copula to absolute values. Could you please give me some direction as to how I should proceed? TLDR: I would like to fit an n-dimensional multivariate Gaussian copula to the CDF.


title: Recommendation System with Graph Database

text: Hi I'm looking into creating a recommendation system based on a graph database for an ecommerce website. Basically something like amazon, "if my customer's have bought these products recommend some others." I was hoping someone might be able to give me some pointers on what types of recommendation systems to check out, hopefully some examples. Also it would be helpful to get some clues on how to deal with the cold start problem, and if there are multiple recommendations how to decide which to use. Any advice is greatly appreciated.


title: Probability of heads in a biased coin

text: Given $N$ flips of the same coin resulting in $k$ occurrences of 'heads', what is the probability density function of heads-probability of the coin?


title: Boosted Regression Trees in R

text: As a general statistics question: When using BRTs should you avoid using strongly correlated variables the way you would in say, GLMs or multiple linear regression? Any useful references would be greatly appreciated. Cheers, Israel


title: Simple multiple choice for test statistic and significance level

text: A study aims to determine if a new tutoring program (method 1) is better than standard preparatory program (method 2) to prepare students for the SAT. It measures the proportion of students who score in the top quarter of test takers. Results: 53% of sampled students using the standard program (method 2) scored in the top quarter of test takers, while 45% of sampled students using the new program (method 1) scored in the top quarter of test takers. A hypothesis test will be performed using a 5% significance level.   1. Will the test statistic be positive or negative?   2. True or false: At a 5% significance level, the results are statistically significant. I put negative for #1 and false for #2, but I got both of them wrong. Can someone explain how I solve this type of question in the future and, most important, avoid choosing incorrectly?


title: Identify if there is a relationship between two data sets

text: I am by no means a mathematician (I am a software developer by trade), but I am trying to find out if there is a relationship between two data sets. I have referral sources the provide us customers, and I am trying to find out if there is a relationship between the number of calls we make to the referral source and the number of customers they give us, and if so, what is the best number of calls to make to yield the most number of patients. I have a year's worth of call and customer data to use, but am not sure how to proceed. If someone could point me in the right direction, I would be very appreciative. Thanks.


title: How to interpret multiple factors in model output in R

text: I am well aware how to read the model summary in R for a regression model when a factor is included. The "first" level, in terms of ABC, is regarded as the base level to which all further levels of that factor are compared to. In an ANOVA-style model the baseline value is found in the intercept (equals the mean response value for base-level class). However, if one or several factors are mixed with continuous predictors, then how can I see what the base-level values are at all? to what would I compare? In the model output below, there are two factors:   1. LandUse (4 Levels)   2. Type_LU (4 Levels) I see now for example that `LandUseLow` is 0.35 units higher than the base- line `LandUseHigh`. Would one now simply look at the mean response for the class `LandUseHigh` and compare? Is it that simple?               My_model:                        Estimate Std. Error Adjusted SE z value Pr(>|z|)     (Intercept)      4.6086772  1.7754606   1.7773711   2.593  0.00951 **     DiTempRange     -0.1409464  0.0764872   0.0765969   1.840  0.06575 .     LandUseLow       0.3520743  0.5989777   0.5997903   0.587  0.55721     LandUseMedium    0.2741413  0.3668149   0.3675811   0.746  0.45579     LandUseNone     -1.0128945  0.5735342   0.5744652   1.763  0.07787 .     MAP              0.0048810  0.0009128   0.0009142   5.339    1e-07 ***     Rivier           0.3502782  0.2252743   0.2257546   1.552  0.12076     TempRange        0.0823410  0.0606546   0.0607942   1.354  0.17560     Tmean           -0.1762862  0.0994486   0.0996353   1.769  0.07684 .     TYPE_LUconserva -0.9312487  0.4770681   0.4781244   1.948  0.05145 .     TYPE_LUprivate  -0.4839229  0.3289011   0.3296201   1.468  0.14207     TYPE_LUstate     0.0004062  0.4079678   0.4089744   0.001  0.99921     logVRM           0.1370973  0.1140166   0.1142342   1.200  0.23008     logTWI          -0.0735540  0.4195267   0.4202589   0.175  0.86106     logDAH           1.7132823  3.5937000   3.6028996   0.476  0.63441


title: Error in fitting negative binomial regression model in R when replicating published results (works in Stata)

text: I'm trying to replicate the results of the first model of this article: > Hultman, Lisa, Jacob Kathman, and Megan Shannon. 2013. “United Nations > Peacekeeping and Civilian Protection in Civil War.” _American Journal of > Political Science_ 57(4): 875–91. Replication material can be found here: http://thedata.harvard.edu/dvn/dv/ajps/faces/study/StudyPage.xhtml?studyId=87987&tab=files The Stata code provided in the do file runs without problems. However, when I try to replicate the model in R. The error strongly resembles the one in this CV question. Here's the R code I've written to replicate the results:               library(MASS)     library(foreign)          pko <- read.dta("HKS_AJPS_2013.dta")          pko_model1 <- glm.nb(osvAll ~ troopLag +                             policeLag +                             militaryobserversLag +                             brv_AllLag +                             osvAllLagDum  +                             incomp +                             epduration +                             lntpop, data = pko, link = log)      This produces the following error:                 Error in glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,  :        NA/NaN/Inf in 'x'      If I include the `control=glm.control(trace=10,maxit=100)` option in `glm.nb` it produces the following output:               Deviance = 75787029 Iterations - 1     Deviance = 28247900 Iterations - 2     Deviance = 11043902 Iterations - 3     Deviance = 4952253 Iterations - 4     Deviance = 2896062 Iterations - 5     Deviance = 2286069 Iterations - 6     Deviance = 2152722 Iterations - 7     Deviance = 2135621 Iterations - 8     Deviance = 2134804 Iterations - 9     Deviance = 2134801 Iterations - 10     Deviance = 2134801 Iterations - 11     theta.ml: iter 0 'theta = 0.000609'     theta.ml: iter1 theta =0.00120256     theta.ml: iter2 theta =0.00234778     theta.ml: iter3 theta =0.00449211     theta.ml: iter4 theta =0.0082914     theta.ml: iter5 theta =0.0143798     theta.ml: iter6 theta =0.0225089     theta.ml: iter7 theta =0.0302781     theta.ml: iter8 theta =0.0342821     theta.ml: iter9 theta =0.0349533     theta.ml: iter10 theta =0.0349683     Initial value for 'theta': 0.034968     Deviance = 3634.951 Iterations - 1     Deviance = 1160161 Iterations - 2     Error in glm.fitter(x = X, y = Y, w = w, etastart = eta, offset = offset,  :        NA/NaN/Inf in 'x'      If I exclude the `epduration` variable or the `incomp` variable, the error disappears and I can roughly replicate the results from the article, but the parameter estimates of course vary because I don't include all variables in the model (and I don't use robust, clustered standard errors in R). Two questions:   1. Why does this run in Stata without any complaint, but not in R?   2. How can I make this work in R? The answers to this question suggest a possible solution by first estimating a Poisson regression and then feeding the results as starting parameter values into an MLE estimation. Yet, I haven't been able to make this work in R.  I realize that this might be a duplicate to the existing CrossValidated question but since I'm having a similar problem with different data, this might be a more general problem.


title: When are $(X, Y)$ jointly normal given that both $X$ and $Y$ are normal?

text: We know that $X$, $Y$ are normal does not guarantee $(X, Y)$ is jointly normal. A typical example is: $X=Z$, and $Y=ZU$, where $Z$ standard normal, $P(U=1)=P(U=-1)=1/2$, and $Z, U$ are independent. My question is: what would be a condition for $(X, Y)$ to be jointly normal given that both $X$, and $Y$ are normal?


title: Summation of conditional probability in Bayesian network

text: I have the joint distribution of the bayesian network defined as $P(CERFD)=P(C)\times P(R)\times P(E|CR) \times P(F|E) \times P(D|F)$ and I am trying to calculate $P(D=d|C=c)$, below are my workings and I am stuck at the last step: $P(D=d|C=c)= \frac{P(D=d,C=c)}{P(C=c)} $ $=\frac{\sum_{ERF}P(C=c)\times P(R)\times P(E|C=c,R) \times P(F|E) \times P(D=d|F)} {\sum_{DERF}P(C=c)\times P(R)\times P(E|C=c,R) \times P(F|E) \times P(D|F)}$ $=\frac{\sum_{ERF}P(C=c)\times P(R)\times P(E|C=c,R) \times P(F|E) \times P(D=d|F)} {\sum_{DERF}P(C=c)\times P(R)\times P(E|C=c,R) \times \sum_F P(F|E) \times \sum_D P(D|F)}$ $=\frac{\sum_{ERF}P(C=c)\times P(R)\times P(E|C=c,R) \times P(F|E) \times P(D=d|F)} {\sum_{ER}P(C=c)\times \sum_R P(R)\times \sum_E P(E|C=c,R) \times 1}$ $=\frac{\sum_{ERF}P(C=c)\times P(R)\times P(E|C=c,R) \times P(F|E) \times P(D=d|F)} {P(C=c)\times 1}$ I don't even know how to simplify the numerator at all... This is not homework, I am learning probabilistic graphical model from coursera and I'm testing my own understanding.


title: How to interpret non-significant effect of a covariate in ANCOVA?

text: I need to use ANCOVA to analyse effect of executive functions (EF-covariate) on scores collected during a test. I want to compare this effect between a control and schizophrenics group (variable diagnoza). I conducted an analysis of covariance and I have some questions. Here is SPSS output for my ANCOVA model: ![](http://i.imgur.com/xA1jg.jpg) Is it a problem that the effect of EF is not significant? Can I interpret this result as follows: There is a significant difference in MIE between groups when the effect of EF is controlled?


title: How should outliers be dealt with in linear regression analysis?

text: Often times a statistical analyst is handed a set dataset and asked to fit a model using a technique such as linear regression. Very frequently the dataset is accompanied with a disclaimer similar to "Oh yeah, we messed up collecting some of these data points -- do what you can". This situation leads to regression fits that are heavily impacted by the presence of outliers that may be erroneous data. Given the following:   * It is dangerous from both a scientific and moral standpoint to throw out data for no reason other than it "makes the fit look bad".   * In real life, the people who collected the data are frequently not available to answer questions such as "when generating this data set, which of the points did you mess up, exactly?" What statistical tests or rules of thumb can be used as a basis for excluding outliers in linear regression analysis? Are there any special considerations for multilinear regression?


title: Time series for count data, with counts < 20

text: I recently started working for a tuberculosis clinic. We meet periodically to discuss the number of TB cases we're currently treating, the number of tests administered, etc. I'd like to start modeling these counts so that we're not just guessing whether something is unusual or not. Unfortunately, I've had very little training in time series, and most of my exposure has been to models for very continuous data (stock prices) or very large numbers of counts (influenza). But we deal with 0-18 cases per month (mean 6.68, median 7, var 12.3), which are distributed like this: ![alt text](http://img827.imageshack.us/img827/1927/activetbcases.png) ![alt text](http://img827.imageshack.us/img827/4348/tbcasedistribution.png) I've found a few articles that address models like this, but I'd greatly appreciate hearing suggestions from you - both for approaches and for R packages that I could use to implement those approaches. **EDIT:** mbq's answer has forced me to think more carefully about what I'm asking here; I got too hung-up on the monthly counts and lost the actual focus of the question. What I'd like to know is: does the (fairly visible) decline from, say, 2008 onward reflect a downward trend in the overall number of cases? It looks to me like the number of cases monthly from 2001-2007 reflects a stable process; maybe some seasonality, but overall stable. From 2008 through the present, it looks like that process is changing: the overall number of cases is declining, even though the monthly counts might wobble up and down due to randomness and seasonality. How can I test if there's a real change in the process? And if I can identify a decline, how could I use that trend and whatever seasonality there might be to estimate the number of cases we might see in the upcoming months? Whew. Thanks for bearing with me.


title: Analysis of speed-accuracy tradeoff in lmer

text: I am analyzing data from an experiment in which participants had to respond to Yes/No-questions of two types as quickly as possible. Hence, there are two dependent variables: response time and correctness of the response. There are 80 responses for every participant in each of the conditions. My hypothesis is that one type of question leads to slower or less correct responses. **Can anybody think of a reason why the following analysis might be flawed?** I determined, for every participant, the probability of giving a correct response before a certain point in time, i.e., $P(RT < t_n \wedge response~correct)$. As far as I can see, this corresponds to the ECDF multiplied by the proportion of correct responses. I removed all data points that equal 0 and fit a linear mixed-effects model with the following specification (in lme4): $logit(P_{corrRT}) \sim time + condition + time:condition + (1|participant)$ As far as I can see, a main effect of condition or an interaction between condition and time would mean that there are differences in either RT or accuracy. Am I making any incorrect implicit assumptions, or does this analysis make sense?


title: Problem with PCA

text: I am trying to do some PC analysis on my data coming from lipids measurements in different samples. I only have one factor: if samples are diabetic or non- dibetic. Here is the PCA graph I get: ![enter image description here](http://i.stack.imgur.com/LbinR.png) As you can see, data are not clustered at all! What should I do for that? Is there any other clustration analysis I can use instead of PCA?


title: Free statistical textbooks

text: Are there any free statistical textbooks available?


title: Logistic Regression: wald chi-square

text: I am running a logistic regression with a binary dependent variable and 5 class independent variables. One information the software I use returns, is the following: Effect / DF / Wald / Pr > ChiSq Var1 / 1 / 150 / <.0001 Var2 / 3 / 119 / <.0001 Var3 / 8 / 157 / <.0001 Var4 / 6 / 1553 / <.0001 Var5 / 4 / 15975 / <.0001 Concerning the data above I have two questions: I) How can I interpret the fact that for var5 the value of the wald chi-square is so much higher than the values of the remaining variables? And is it a bad indicator of the regression quality? II) I saw the impact of a regression variable calculated as the division between the wald chi-square value of the variable and the total sum of all the variables wald chi-square values (shown in percentage….is this example var5 would represent 88.9% of the probability resulting from the logistic regression), but found no statistical grounds for this methodology. Anyone knows this methodology or others? Any guidance will be helpful


title: Measure of similarity/distance of data points in geographic space

text: Given two points $p_1=(x_1,y_1,t_1)$ and $p_2=(x_2,y_2,t_2)$, where $x$ and $y$ refer to the geographic coordinates in the plane, and $t$ to some measured value. Two distance measures to evaluate the similarity between these two points come to my mind: $$d_1(p_1,p_2) = \sqrt{ (x_1-x_2)^2+(y_1-y_2)^2+(t_1-t_2)^2 }$$ $$d_2(p_1,p_2) = \sqrt{ (x_1-x_2)^2+(y_1-y_2)^2 }+ \sqrt{ (t_1-t_2)^2 }$$ Measure $d_1$ is simply Euclidean distance in 3d-space, while $d_2$ is the sum between spatial distance and attribute distance. Which measure does make more sense and should be applied for e.g. clustering?


title: What languages are commonly used in medical statistics?

text: A doctor friend has asked if I would teach them some basic programming, because they think it will be useful when they start to do research in the future. To give some background - they don't know what language their lab-mates will use (because they don't know what lab they are joining) but they would like to learn something that will be generally useful. They are not interested in doing statistical research, but rather applying statistical techniques to medical data. There are some obvious options (I know enough about all of these to teach them)   * Python   * R   * MATLAB and a few less obvious ones (I don't know anything about these)   * Stata   * SAS   * SPSS I'd like to know what languages are most commonly used in medical statistics - do people use general purpose languages like Python, or more specialized languages like R and MATLAB, or statistical software like Stata/SAS/SPSS? I should clarify that I'm thinking of someone whose primary job is as a doctor/medical researcher who will be applying statistics, not someone who is a professional statistician or software developer who happens to be working in medicine (so e.g. I think that C/C++/Java would all be poor choices).


title: R auto.arima() with non-stationary covariates

text: I want to fit an ARMA model with covariates to a non-stationary time series. I have daily measurements for water flow for 4 stations (S1-S4) and the time series is not stationary, so I will have to take a first difference. I want to find an ARMA model that takes as covariates the other 3 stations. What I have done:   a) m1 = auto.arima(S1, d=1, max.p, max.q, xreg=data.frame(S2,S3,S4))   b) m2 = auto.arima(diff(S1), d=0, max.p, max.q, xreg=data.frame(diff(S2),diff(S3),diff(S4))) My question: which model is correct, considering the series are non- stationary? I will appreciate any suggestion. Thank you,   Patricia


title: Generalized linear mixed models: model selection

text: This question/topic came up in a discussion with a colleague and I was looking for some opinions on this: I am modeling some data using a random effects logistic regression, more precisely a random intercept logistic regression. For the fixed effects I have 9 variables that are of interest and come into consideration. I would like to do some sort of model selection to find the variables that are significant and give the “best” model (main effects only). My first idea was to use the AIC to compare different models but with 9 variables I was not too exciting to compare 2^9=512 different models (keyword: data dredging). I discussed this with a colleague and he told me that he remembered reading about using stepwise (or forward) model selection with GLMMs. But instead of using a p-value (e.g. based on a likelihood ratio test for GLMMs), one should use the AIC as entry/exit criterion. I found this idea very interesting, but I did not find any references that further discussed this and my colleague did not remember where he read it. Many books suggest using the AIC to compare models but I did not find any discussion about using this together with a stepwise or forward model selection procedure. So I have basically two questions:   1. Is there anything wrong with using the AIC in a stepwise model selection procedure as entry/exit criterion? If yes, what would be the alternative?   2. Do you have some references that discuss the above procedure that (also as reference for a final report? Best, Emilia


title: random forest classification in R - no separation in training set

text: Originally posted on Stack Overflow, but suggested to move here... I'm new to machine learning, but I've tried to perform a Random Forest classification (randomForest package in R) on some metabolomics data with bad results. My normal approach in this case would be to employ a PLS-DA strategy. However, I decided to try both RF and SVM as there are some publications highly recommending these machine learning approaches for Omics data. In my case, 'X' is a 16*100 data frame (16 individuals with 100 recorded features/predictors) read from a CSV file. 'Y' is a factor vector (length=16) with 8 'high' and 8 'low'. In both PLS-DA and SVM (both linear and radial kernel) I get excellent separation. However, I get 3 misclassifications out of 16 in the RF model. The RF model looks like: RFA1=randomForest(X,Y)               ## read file and fix data frame     in.data = read.csv2(file='Indata.csv', header = FALSE, skip=5)[,-4] # Col 1-3 are identifiers. Features/predictors from col 4     names(in.data)=names(read.csv2(file='Indata.csv',header=T)[,-4])     # str(in.data)      # $ ID                       : Factor w/ 27 levels "2","3","4","5",..: 2 3 4 6 8 10 20 23 5 11 ...          # $ period                   : Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...      # $ consumption          : Factor w/ 2 levels "high","low": 1 1 1 1 1 1 1 1 2 2 ...          # $ FEATURES...          ## Sort DF into X (features) and Y (classifier based on consumption)     y = in.data$consumption                   # Classifier based on high vs low consumption     x = in.data[,-1:-3]                       # 100 features/predictors into X NB Contains many NAs     nr=nrow(x)     nc=ncol(x)     x.na = as.data.frame(is.na(x))            # Find NAs in X     col.min=apply(x,2,min,na.rm=T)            # Find min value per feature (omitting NAs)     ## Deal with zero/missing data-situation     x2=x                                      # Compute new x2 matrix without NA     for (i in 1:nc) {         x2[x.na[,i],i]=col.min[i]             # Substitute missing data with col.min     }          ## Make classifiers according to period (A vs B)     a.ind = in.data$period=='A'         b.ind = in.data$period=='B'          ## Choose data from period A only & transform/scale X     x2a=x2[a.ind,]                 # Original data     x2a.scale=scale(x2a)           # Scaled     x2a.log=log(x2a)               # Log-transformed     x2a.logscale=scale(log(x2a))   # Log-transformed and scaled     ya=y[a.ind]          ## Perform analysis for period A     library(randomForest)     (rfa1=randomForest(x2a,ya))     (rfa2=randomForest(x2a.scale,ya))     (rfa3=randomForest(x2a.log,ya))     (rfa4=randomForest(x2a.logscale,ya))      This generates output like:               Call:      randomForest(x = x2a, y = ya)                     Type of random forest: classification                      Number of trees: 500     No. of variables tried at each split: 10                  OOB estimate of  error rate: 18.75%     Confusion matrix:          high low class.error     high    6   2       0.250     low     1   7       0.125      I have played around with both mtry (5-50) and ntree (500-2000) with no apparent success. I've also tried combinations of transforms and scaling of 'X'. But as I understand it, RF is a non-parametric method and as such, transformations and scaling won't do anything for the results. For comparison, using the exact same data, PLS-DA using SIMCA13 provides excellent separation already in the 1st component. SVM using the kernlab package in R provides 0 training error. At this stage I'm not looking at validation or using test sets. I want to first make sure I get good classification on my training set. I'm sure I'm missing something, but I don't know what. I hope to have supplied sufficient information to describe the problem. Thanks in advance for any help! Sincerely, Calle


title: Why is correlation obtained from nlme different from Pearson correlation?

text: Fixed effects: Ratio ~ ADF + CP + FCM + DMI + DIM                       Value  Std.Error  DF   t-value p-value     (Intercept)  3.1199808 0.16237303 158 19.214896  0.0000     ADF         -0.0265626 0.00406990 158 -6.526603  0.0000     CP          -0.0534021 0.00539108 158 -9.905636  0.0000     FCM         -0.0149314 0.00353524 158 -4.223598  0.0000     DMI          0.0072318 0.00498779 158  1.449894  0.1491     DIM         -0.0008994 0.00019408 158 -4.634076  0.0000          Correlation:          (Intr) ADF    CP     FCM    DMI        ADF -0.628                                 CP  -0.515  0.089                          FCM -0.299  0.269 -0.203                   DMI -0.229 -0.145  0.083 -0.624            DIM -0.113  0.127 -0.061  0.010 -0.047      I obtained these value from `lme` function for mixed model. However, I am wondering why the values for correlation of fixed effect provided here differ from what I obtain using `cor.test` function (Pearson correlation)


title: How exactly does one “control for other variables”?

text: Here is the article that motivated this question: Does impatience make us fat? I liked this article, and it nicely demonstrates the concept of “controlling for other variables” (IQ, career, income, age, etc) in order to best isolate the true relationship between just the 2 variables in question. Can you explain to me **_how_** you actually control for variables on a typical data set? E.g., if you have 2 people with the same impatience level and BMI, but different incomes, how do you treat these data? Do you categorize them into different subgroups that do have similar income, patience, and BMI? But, eventually there are dozens of variables to control for (IQ, career, income, age, etc) How do you then aggregate these (potentially) 100’s of subgroups? In fact, I have a feeling this approach is barking up the wrong tree, now that I’ve verbalized it. Thanks for shedding any light on something I've meant to get to the bottom of for a few years now...!


title: How does controlling for variables work in multiple regression?

text: My question is not about statistical programming. We know how to code the software in a regression with many independent variables. My question is about how the computer software controls for all variables at the same time mathematically. Simple correlation of two variables is easy to understand, as each variable has the other one's full attention. In multiple regression, however, multiple independent variables are "looking at" the dependent variable at the same time, and each one looks to "subtract out" or explain the variation, but they can't all be in operation with the dependent variable at the same time. If a single given independent variable is given priority to be the first to interact with the dependent variable, and it "takes away" its share of the dependent variable, then what it takes away won't be rightfully eligible to be explained by other independent variables, as it justly should. My intuitive sense of what happens must be mistaken, and I need an answer that does not involve formulas. The challenge for the teaching profession is how to explain complex ideas narratively or metaphorically for analysts that don't have a fast, intuitive sense for symbolic explanation. A related question is this: is the process iterative in that initial approximations (of betas) are made, and they incrementally adjusted as each of the other independent variables add more information to the initial beta estimates?


title: Derivation of conditional distribution from other two distributions

text: $$Y|X=x \sim N(x,1)\\\X\sim N(\mu,\sigma^2 )$$ What distribution does $X|Y=y$ follow? My initial startegy was to $f_{Y|X}f_X=f_{X,Y}$ and solve for $f_{X|Y}=f_{X,Y}/f_{Y}$ . Computing for $f_{X,Y}$, I get the following: $$f_{X,Y}=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\\{-\frac{(x-\mu)^2}{2\sigma^2}\right\\}\frac{1}{\sqrt{2\pi}}\exp\left\\{-\frac{(y-x)^2}{2}\right\\}$$ And in trying to compute for $f_{Y}$, I was trying to integrate above w.r.t $x$, but I was stuck. I am not sure if that's integratable, and if this is a right approach to solve this question. I am curious if there is some kind of tricks/insights I am missing.


title: Correlation using JMP & Excel

text: Real life situation. I have data from different cities saying on average how fast they build a home. I want to correlate this to the amount of construction activity they do & number of homes constructed per construction personnel. I also have data saying total # of construction personnel & total construction cost in respective cities. I have ran all this in JMP & Excel and am getting a correlation of up to .20, and if I take out outliers I have correlations up to .60. My question is: how can I get correlation up to .75 or above, if possible? I am also trying to work with multivariate regression and see if I can get better results.


title: Finding the maximum of unparameterized 2D data

text: I have a set "Z" of values, one for each combination of (X,Y), given on a grid. I would like to estimate the location of the maximum of Z (expected to lie near the center of the dark blue pixels). The maximum need not lie on one of the discrete (x,y) pairs defined by the vectors. The array was created such that there is not really a direct way to define the analytic function for Z(X,Y). Currently, the only thing that I know to do is marginalize over X and Y, respectively, fit a 1-D function (perhaps a Gaussian) to each in order to get an explicit analytic expression, minimize along each direction independently, and then just use the x^ and y^ of the minimum values for each of these as the coordinates of the 2D minimum. Is there a better way to do this (conceptually)? Is there a way to do this using the Z array directly, that does not rely on approximating X and Y as specific distributions? Even in this case, the minimum does not have to fall on the grid defined by X and Y. Thank you! * * * _**Additional Information regarding Z_** :   * Z is essentially an 2D array of log-likelihood values, and should not have zeros; all elements of Z are less than zero.


Why can a polynomial of degree $>2$ not be a cumulant generating function? I read somewhere that this is impossible but can't retrieve the source. The answer by StasK to Higher order generalization of the multivariate normal distribution mentions a related statement ''Once you depart from zero third cumulant, all higher order cumulants have to be non-zero, as well: there is no distribution for which $\kappa_4=0$ if $\kappa_3\ne 0$," but also gives no source.


title: R model.tables() incorrect means – possible bug?

text: I was doing ANOVA in SPSS and then in R and to my huge surprise, the results from the latest version of R were incorrect. When I use the function `model.tables(x,"means")` to get descriptive statistics, the independent variable means by the second dependent are slightly incorrect (e.g. 129 instead of 130.27). My question is what could cause the problem? I am a novice to R but using the same data, SPSS gets the result correctly, so something is obviously wrong. `head(data)`:                 skupina pohlavie zodpovedny     1       1        1        152     2       1        1        118     3       2        2         88     4       2        1        140      Code:               x <- aov(zodpovedny ~ pohlavie*skupina,data=data)     model.tables(x,"means")      Problem illustrated: ![This is unfortunate.](http://i.stack.imgur.com/g19cG.jpg)


title: Determining whether two samples come from the same population and which hypothesis test is appropriate

text: I have a study conducted by a team of doctors: two groups of people were treated for stress urinary incontinence in two ways: TVT and abbrevo. Each group consisted of 23 people. The hypothesis is as follows: **neither method is better than the other.** To begin, I want to check the comparability of groups. That is, to check whether these two samples are a random sample from the same population. **Which criteria is best to use?** Secondly, **how do I formulate a strictly statistical hypothesis test and what method?** Initially I inclined to the Mann-Whitney U test, but Wikipedia says that this criteria can not be used when a lot of duplicates. We have information about   1. Age of patient   2. Patients are females   3. Amount of births and pregnancies   4. Mass index   5. Menopause(yes/no)    6. **UDI-6 (before operation)**   7. **IIQ-7 (before operation)**   8. **Pain (yes/no, before operation)**   9. Diagnosis (stress urinary incontinence, recurrent incontinence, or mixed form)   10. Prolapse during op. (yes/no)   11. Complications during operation   12. Postoperative complications   13. Complication (after 6 month)   14. Cough test (after 6 month)   15. **Pain (yes/no)**   16. **UDI-6**   17. **IIQ-7**


title: Comparison of Cox models

text: I am new to statistics and SPSS v 20. I am doing survival analysis - KM and Cox regression with SPSS. I am trying to make Cox models with different variables - clinical, genetic etc. I am interested to know how well each of these models are able to predict the occurrence of metastases. My questions are:   1. Is it possible to do this with SPSS?   2. Is there a way to save the model predicted outcome while doing Cox regression? I saw that when I do multinomial logistic regression in SPSS, I can save the predicted category and compare it with the actual category. I don't have that option in Cox.   3. Is it possible to make ROC curves for different Cox models and then compare the area under the curve for different models? I am not sure I am going about this the right way.


title: Can AIC determine which data better fit the same model?

text: Many moons ago, I asked how to differentiate between two very similar non- linear fits and which was better. Finally got that all straightened out after many headaches and three different software packages -- calculating AIC, AICc, and BIC and all give a consistent result. Now I'm trying to apply a slight modification to the data and re-fit. The question is, since $\chi^2$ is fairly useless on very non-linear models, can I use AIC to tell me if the original data or new data are better fits to that model? As in, I use an $f(x)$ model, fit it with dataset1, and then I fit the SAME model to dataset2, and I calculate AIC/AICc/BIC of both fits. Is that a meaningful comparison, or can they only be used to rate between two different fits to the same data?


title: Code Problem in Lognormal Bayesian Model in JAGS/WinBUGS

text: I have a data set with 853 units. Below, a small sample of it: > y                [1]  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000     [10]  0.000000  6.751283  0.000000  0.000000 11.927481  0.000000  0.000000  0.000000  0.000000     [19]  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000     [28]  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000  0.000000     [37]  0.000000  0.000000  0.000000  0.000000  5.342471  2.286760  0.000000  0.000000  3.464403     [46]  3.952309 16.030779 17.942967  0.000000  0.000000      > head(x)                   gdp.cap     bf.cap     1  8.789771 0.03911118     2 10.204732 0.02341257     3  6.890112 0.04160150     4  9.178957 0.03892352     5  9.384171 0.03880186     6  9.562188 0.04165802      I'm calling JAGS from R to run the following code: 1) In JAGS               model {              # Lognormal likelihood         for(i in 1:N) {              y[i] ~ dlnorm(mu[i], tau)             mu[i] <- beta0 + inprod(x[i,], beta[])         }              # Prior for betas         beta0 ~ dnorm(mbeta0, precbeta0)              for (j in 1:K) { # b are the number betas to be estimated             beta[j] ~ dnorm(m[j], prec[j])         }              # Prior for precision of Y         tau ~ dgamma(tau.a, tau.b)         sigma <- 1/sqrt(tau)          }      2) In R:               dt <- read.csv("data.csv", header=T, sep=',')     dt <- na.omit(dt)          # Organize Data     y <- dt[,30]     x <- as.matrix(dt[,c(31,46)])          # Data:       data <- list(N=length(y), K=2, y=y, x=x, m=c(-2,2), prec=c(.20, .20), tau.a=1, tau.b=2, mbeta0=0, precbeta0=.01)     inits <- rep(list(list(beta0=0, beta=c(1,1), tau=1)),5)     param <- c("beta0", "beta", "sigma")          # Sample Posterior     library(rjags)     sim <- jags.model(file="Bayesian/lognormal.bug", data=data, inits=inits, n.chains=5, n.adapt=1000)      I don't know why I get the following error:               Compiling model graph     Resolving undeclared variables     Allocating nodes     Graph Size: 39212          Initializing model     Deleting model          Error in jags.model(file = "Bayesian/lognormal.bug", data = data, inits = inits,  :      Error in node y[1]     Observed node inconsistent with unobserved parents at initialization      Any help, please?


title: Adding a history flag variable to panel data with possible auto-correlation

text: I have a pretty basic question. A client has given me a data set with tens of thousands of subjects. Each subject has from one to about ten separate records. These records cover successive non-overlapping time periods of variable duration, up to one year. There is a response variable, Y >= 0. The client is interested in seeing, among other things, whether non-zero response within the previous three years has any predictive value / correlation with response in the current year. We have decided to create a binary variable, A, with the value of 1 if Y > 0 for any record in the previous three years (for each individual subject), 0 otherwise. Note that the oldest record for each subject will always have A = 0. If most of my subjects had long historical trails, then I'd be willing to drop the first three years for each subject, and just do the analysis of Y versus A for all the records which have a full three year history. However, for this data set, most individuals don't have three years of history, so many of the A variable values are partially truncated. Of course, this could bias the analysis. What is the best way to carry out this analysis? Thanks!


title: MLE estimation of spatial effects radius

text: I am trying to identify the maximum likelihood estimates in an SDM model (a hedonic home price model, with observations being 5,000 individual homes), $$y=\rho W y + X\beta + WX\lambda + \epsilon$$ The likelihood function of this model is $$\ln(\mathcal{L}) = \dfrac{n\ln(\pi\sigma^2)}{2} + \ln|I - \rho W| + \dfrac{e'e}{2\sigma^2} $$ with $e= y - \rho W y - X\beta - WX\lambda$ I have been able to successfully estimate parameters $\rho, \beta,\lambda$ using the `lagsarlm` function in the `spdep` package. In my case, $W$ is an inverse distance matrix limited to a radius $r$, or $$W_{ij} = \begin{cases} 1/d_{ij} & \text{for $d_{ij} \leq r$ miles}\\\ 0 & \text{for $d_{ij} > r$ miles} \end{cases}$$ I would like to identify $r$ simultaneously with $\rho, \beta,\lambda$. I have even written a likelihood function to do this in R (making use of the neighbors and lag functions in `spdep`)               sdm.lik <- function(params, y, X, GeoPoints){       n <- length(y)       rho <- params[1]       radius <- params[2]       sigma2 <- params[3]       delta <- params[-c(1,2,3)]       # make W subject to radius       points.dnn <- dnearneigh(GeoPoints, 0, radius*5280)       dists      <- nbdists(points.dnn, GeoPoints)       dists.inv  <- lapply(dists, function(x) 1/x)       Wlist <- nb2listw(points.dnn, glist=dists.inv, style="W", zero.policy=TRUE)          W <- listw2mat(Wlist)       # Lag X variables       Z <- cbind(X[,1], lag.listw(Wlist, X[,1], zero.policy=TRUE))        for(i in 2:ncol(X)){         Z <- cbind(Z, X[,i], lag.listw(Wlist, X[,1], zero.policy=TRUE))       }       # Calculate log likelihood       e <- y - ((rho * W) %*% y) - (Z %*% delta)       lagmatrix <- diag(1,nrow=length(y), ncol=length(y)) - rho*W       logL <- ( -0.5*n*log(pi * sigma2) + log( det(lagmatrix) )                                                -(t(e) %*% e) / (2*sigma2) )       return(-logL)     }      Optimizing this function is far too slow to be practical, however. Unsurprisingly, the highest cost operation is the log-determinant.   1. Do you have tips for computing this determinant more quickly?   2. Does `lagsarlm` already have functionality to do this?   3. Can I hack it to make it do this?


title: Factor analysis with severly skewed ordinal data and censored ordinal data

text: I am aiming to run initial exploratory factor analysis in one sample and then confirmatory in another sample. My indicators are ordinal and so I planned to generate a polychoric correlation matrix and use this when estimating the models in R and openMx. I am aware though that the polychoric correlations assume an underlying normal distribution. Most of my variables are highly skewed (i.e., I am measuring low frequency behaviours and so there is a large preponderance of zeros) and some variables have clear right censoring (i.e., the last category is often something like "more than x times"). I am wondering whether there is any recommendations for dealing with this? I have a massive sample (15000) so is it justifiable to use the central limit theorem and not worry about the observed distribution of the indicators (and just use polychoric correlations)? Will I have biased or invalid parameter estimates if I just use the polychoric correlations anyway (could I defend this?). I am reluctant to use IRT because I believe the data is multidimensional with factors being at least reasonably correlated. I am also predicting that there is a second order factor that will explain all the correlated subfactors. I also have no firm a priori predictions about how the items will load on unidimensional factors (this is my understanding of IRT, i.e., that it assumes unidimensionality but happy to be instructed otherwise). Any guidance or correction of me reasoning would be greatly appreciated. Is there a solution to this in R? Cheers!


title: How to test counter factuals

text: Let's say I have the following data on leads, monthly media spend, and clicks               Month  Leads   Media     Clicks     Jan     150    1000       500     Feb     200    1000       550      March   300    1200       800     ...      Let's say I run a linear regression where y is leads and the predictors are media and clicks. That's good, I know the relationships between these variable and can generate some lags to produce predictions. But what if I had spent 500 (or 2000 or 0) on media, what would have occurred. How do I perform this type of 'counter-factual' analysis where I attempt to find the results of a model if the actual value from one or two of the predictors was lower or higher? What is the standard approach (aka statistically proper approach)? Is it just a matter to "adjusting" the data to the 'new' number and rerunning the regression? or maybe simulating a regression 100+ times with 100+ different values for media?


title: Checking for seasonality in data set

text: I'm working on housing data set for a project for a undergraduate time series course. I'm trying to see if there is seasonality in the data. I used the following commands but do not know how to read the graphs to see if there is seasonality.               spec.pgram(data) #non-parametric estimate of the spectrum     spec.pgram(data, log="no")      ![enter image description here](http://i.stack.imgur.com/t1o8m.png) This is the data set I'm working with. http://www.quandl.com/FRED/HOUST- Housing-Starts-Total-New-Privately-Owned-Housing-Units-Started


title: Estimator of autocovariance in a wide-sense stationary process

text: From Wikipedia http://en.wikipedia.org/wiki/Ergodic_process > One can discuss the ergodicity of various properties of a stochastic > process. For example, a wide-sense stationary process $x(t)$ has mean > $m_x(t)= E[x(t)]$ and autocovariance $r_x(\tau) = E[(x(t)-m_x(t)) > (x(t+\tau)-m_x(t+\tau))]$ which do not change with time. One way to estimate > the mean is to perform a time average: $$ \hat{m}_x(t)_{T} = \frac{1}{2T} > \int_{-T}^{T} x(t) \, dt.$$ > > If $\hat{m}_x(t)_{T}$ converges in squared mean to $m_x(t)$ as $T > > \rightarrow \infty$, then the process x(t) is said to be mean-ergodic[1] or > mean-square ergodic in the first moment.[2] > > Likewise, one can estimate the autocovariance $r_x(\tau)$ by performing a > time average: $$ \hat{r}_x(\tau) = \frac{1}{2T} \int_{-T}^{T} > [x(t+\tau)-m_x(t+\tau)] [x(t)-m_x(t)] \, dt.$$ > > If this expression converges in squared mean to the true autocovariance > $r_x(\tau) = E[(x(t+\tau)-m_x(t+\tau)) (x(t)-m_x(t))]$, then the process is > said to be autocovariance-ergodic or mean-square ergodic in the second > moment.[2] > > A process which is ergodic in the first and second moments is sometimes > called ergodic in the wide sense. I wonder if $m_x$ should be $\hat{m}_x$ instead in the definition of $\hat{r}_x(\tau)$, i.e. $$ \hat{r}_x(\tau) = \frac{1}{2T} \int_{-T}^{T} [x(t+\tau)-m_x(t+\tau)] [x(t)-m_x(t)] \, dt$$ should be $$ \hat{r}_x(\tau) = \frac{1}{2T} \int_{-T}^{T} [x(t+\tau)-\hat{m}_x(t+\tau)] [x(t)-\hat{m}_x(t)] \, dt$$ Thanks!


title: Error in books formula concerning the log-rank test?

text: I am currently reading the following book:   * Clinical Trials: A Practical Guide to Design, Analysis, and Reporting by Duolao Wang, Ameet Bakhai (2006) On page 243-244 there is an example about the log-rank test. The log-rank test statistic is denoted by the chi square symbol and calculated in the following way: ![enter image description here](http://i.stack.imgur.com/cQLEk.png) and further on it is given: ![t2](http://i.stack.imgur.com/W3nWj.png) and ![t2](http://i.stack.imgur.com/RiI9z.png) I do not want to focus on the idea behind these formulas and what each symbol stands for. I give you the following table, so you do not need to calculate the $d_j$ and the $n_{1j}$ resp. $n_j$. ![t4](http://i.stack.imgur.com/9Hed2.png) It is maybe difficult to read, but the Treatment A colmuns are referred to index 1 and the Treatment B columns are reffered to index 2. **In my opinion the given formula to calculate the $e_{1j}$ is wrong.** E.g. let's consider the first 2 j's: $e_{11}=\frac{n_{11}-d_1}{n_1}=\frac{17-1}{34}=0.47$ and this is $!=0.5$ $e_{12}=\frac{n_{12}-d_2}{n_2}=\frac{16-1}{31}=0.484$ and this is $!=0.5161$ In my opinion, the correct formula is: $e_{1j}=\frac{d_j}{n_j}*n_{1j}$ I check this for Time (months) =23 (j=7) and 27 (j=8): $e_{17}=\frac{d_7}{n_7}*n_{17}=\frac{1}{20}*13=0.65$ $e_{18}=\frac{d_8}{n_8}*n_{18}=\frac{1}{19}*13=0.6842$ And indeed, these are the values given in the table. **What do you say?** EDIT: Of course, I should say that the table is from the book.


title: R2jags does not remove the burn in part sometimes?

text: I found out that the function `jags()` in the `R2jags` package sometimes does not remove the burn in part even with the option `n.burnin=##`. Here is a very simple example in `R` (a simple linear model):               library(R2jags)          N <- 1000     y <- rnorm(N)     x <- rnorm(N)     data <- list("N", "y", "x")          inits <- function(){list(beta0=rnorm(1), beta1=rnorm(1), tau=1)}     parameters <- c("beta0", "beta1", "tau")      the model `m.bug` is like this:               model{     for (i in 1:N){     y[i] ~ dnorm(mu[i], tau)     mu[i] <- beta0 + beta1*x[i]     }          beta0 ~ dnorm(0, 0.00001)     beta1 ~ dnorm(0, 0.00001)     tau ~ dgamma(0.001, 0.001)     sigma2 <- 1/tau     }      using "jags()" in `R2jags` package like this:               m <- jags(data, inits, parameters, "m.bug",        n.chains=3, n.iter=2000, n.burnin=1000, n.thin=1)      My question: In the output of `m`, the posterior estimates are based on the right number of interations (1000), but if we check the traceplot (using `traceplot(m)`), the burnin part seems still there(e.g. the first few values for "tau" are not converged). Why? and there is only one progress bar (uaually two, one for burn-in, one for the rest). also, if I change `n.iter=2000, n.burnin=1000` to `n.iter=2001000, n.burnin=2000000` the elapsed time does not change, which is "too fast" for so many iterations. ps. I used `R` version 2.15.2 and `R2jags` version 0.03-08


title: When does AIC lose its power to discriminate models?

text: There are two simple questions at the end, but I think it is also useful to share the background that motivated them. It comes from this question on an unexpected forecast from the fully automatic methodology behind the forecast::ets function in R. The code, plot and forecast are given below for convenience:               library(forecast);options(scipen=999)     usage <- ts(scan('http://cl.ly/102L0j3o1p2m0m3p0t2o/usage'), frequency = 24)          plot(f1<-forecast(m1<-HoltWinters(usage), h = 168))     plot(f2<-forecast(m2<-ets(usage), h = 168));AIC(m2) #replication OK      ![enter image description here](http://i.stack.imgur.com/CJ4N3.png)               plot(f3<-forecast(m3<-ets(usage,additive.only=TRUE), h = 168))     plot(f4<-forecast(m4<-ets(usage,additive.only=TRUE,damped=FALSE), h = 168))      Just by looking at the plot of the data it appears to me immediately that ARIMA(0,2,2) or ETS(AAN) will be among the best non-seasonal models (and their point forecasts will not differ much). Following the usual advice that the AIC of only a small set of potentially useful models should be compared, I can see no reason to consider multiplicative models here, nor can I see how a damped forecast will be useful for me. With this information in, the "best" ets model m4 and its associated forecast f4 is what was expected, but the process is not a fully automatic one. With many time series, which I would not have time or desire to look at, I would hardly have a better option than to blindly use ets(data). The ets documentation page assures me that _"The methodology is fully automatic. The only required argument for ets is the time series. The model is chosen automatically if not specified"_. The above example is not the first one to show that the methodology that works is at best semi-automatic and the fully automatic one may fail to work as expected even in simpler cases. One way to rectify the problem is to consider n $\ge$ 1 models with suitably low values of AIC as equally supported by the data. I agree with RJH ( _"Although the model may not give the best AIC, it may give forecasts that are more useful to you."_ ) that the most useful forecasts need not be those from the model with the smallest value of the AIC. But then we have the question of when the AIC is useful for selecting a model that forecasts best. It can be deduced from, for example, quotes in this question that if the AICs of all models considered differ by no more than some small number c, then the AIC loses its power to distinguish between these models and other factors/ideas must be considered. But how small is c? Those working with Akaike (Y Sakamoto and M Ishiguro and G. Kitagawa) in the book entitled "Akaike Information Criterion statistics" (in the section entitled "Some remarks on the use of the AIC") mentioned c=1. The number c=2 is often mentioned (e.g. a quote from Brian Ripley can be added to those two linked to above). The number c=4 was mentioned by Chris Chatfield in one of his books. I have not seen anything explicit on the values of c greater than four, but this probably depends on the variability of data, sampling error of the deviance and related factors. In the above example, AIC(m3,m4) suggests that the AIC of the model with a more useful forecast is greater than the AIC of the other model by "only" 78348.75-78337.22=11.52. Are there any formulae, guidelines or rules of thumb for useful values of c given data? Have values of c greater than four been mentioned in the literature?


title: income elasticity + regression + household survey data

text: I'll try to be concise and up to the point: **Problem** : estimate income elasticities for various products from household survey data. **Given** : I have to big complex household survey datasets that uses two stage stratified sampling. However, I have the design specified (strata variable) for only one, and replicate weights for the other data set. **The idea** : is to estimate some for of a regression `expenditure_for_product ~ household's income` and then use the coefficients to estimate it's elasticity **What is unclear** : To estimate simple statistics as mean expenditure for a product or its S.E. one has to properly account for the survey design (especially for variance estimation). But I am not sure how if at all I need to adjust it in my regression case because the literature survey in fact does not provide a clear answer ( whether the weights (and correspondingly the survey design) should be used in my regression case) **Software** : there is in fact R package `survey` which is capable of incorporating (not sure exactly how, as help files does not specify or I can't locate it) the survey design into regressions but this of course does not mean I actually need to do it in my case. **A small digression** : For example, say that `log(expenditure_for_product) ~ log(household's income)` captures the relationship nicely. In that case, the regression coefficient would be all that I need. However, as there are many zero purchases reported by households for goods, I cannot model it in logs. I had an idea of aggregating the data into percentiles so that there would be 100 points with `(average income of x percentile, average expenditure for some product of x percentile)`. But how can if at all should adjust this idea for survey design to get correct estimates? **In essence** : In essence I am looking for a theoretically correct* way to calculate income elasticities given complex household survey data and thus references on this matter or suggestions would be very helpful. *Note, I cannot incorporate price effects because I do not have price information available, so I realize that it will not be too precise. Thank you!


title: Gaussian Process , selecting the hyperparameters

text: I am using Gaussian Process regression toolbox from the site http://www.gaussianprocess.org/gpml/code/matlab/doc/ I was able to use implement the code in matlab easily, following the guide lines. However I am unable to get a coo-relation coefficient value. I have a training data set with n=64 and d=5. I use the following mean, covariance and likelihood functions.The testing data set is with n=30, d=6.               mean = {@meanPoly ,2}; hyp.mean = [1;1;1;1;1;1;1;1;1;1;1;1]     cov = {'covSEiso'};ell=1/4, sf=1 hyp.cov = log ([ell;sf])     lik = @likGauss; sn = 0.1; hyp.lik = log(sn)      Then using minimize function as in demoregression provided               hyp = minimize(hyp, @gp, -50, @infExact, mean, cov, lik, x, y);     [ymu ys2 fmu fs2 lp] = gp(hyp, @infExact, mean, cov, lik, x, y, xs, ys)      I have tried varying the hyper parameters manually but its time consuming and hasn't led to any fruitful results till now. I would really appreciate if you could show me error I am making in implementing the code.


title: Hausman-Newey test for serial correlation in Poisson with Fixed Effects

text: The article from Hausman, Hall, and Griliches (1984) "Econometric Models for Count Data with an Application to the Patents-R&D Relationship" has become the canonical example for conditional MLE of count panel model with fixed effects. In the Appendix B, Hausman and Newey develop a test for serial correlation for this fixed effects model. I need to implement or find an implementation of this test, and I'm having some trouble figuring out how to calculate $\hat{V}_m$, the asymptotic variance-covariance of this vector $m_i(\beta)$. The last sentence says that B.3 is analogous to a Lagrange multiplier test and can be computed via a regression. What regression? It's not at all clear to me. Barring that I'm trying to compute $\hat{V}_m$. The piece I'm missing now is how to compute $U_i{(\hat{\beta})}$. Is this just the score vector of the multinomial likelihood in 2.5? I think I must be mistaken, because it doesn't seem to me that this conformable with $m_i(\hat{\beta})$? Link to the paper


title: Why is the square root of the product of two regression coefficients, one from X predicting Y and vice versa equal to the correlation?

text: If there's a regression model where $Y = a + bX$ where $a = 1.6$ and $b=0.4$, it has a correlation coefficient of $r = 0.60302$. If $X$ and $Y$ are then switched around and the equation becomes $X = c + dY$ where $c=0.4545$ and $d=0.9091$, it also has an $r$ value of $0.60302$. I'm hoping someone can explain why $r=(d\times b)^{0.5}$ is also $0.60302$.


title: Trouble using pareto levy stable distribution software

text: I'm using an arcane free program off the internet called "stable.exe" trying to fit a stable distribution curve to a dataset, but I'm having trouble entering the dataset file into the program. When the program asks me to input the file, I enter "stable test 1", which is the name of the dataset (a notepad document consisting of a column of 20 numbers, in the same working directory as the program). It then spits out an "error reading file" message. I've tried specifying the directory as `C:\...etc...\stable test 1`, and so on, to no avail. If someone has had experience with this nifty yet somewhat enigmatic program, or has fortran knowledge in general, I would appreciate your advice. Also, the read me file for the program says that the program will not read the last number in the dataset unless I include a "carriage return" after it. Can someone explain what a "carriage return" is.


title: Paired multiarm bandit

text: I have a set of independent experiments with different distributions and I'm trying to determine which has the highest mean payoff. I would like to treat this as a multi-arm bandit problem, but the setting is a bit unusual. At each iteration I choose 2 of the experiments, each randomly produces a set of items, then a single item is chosen out of the combined item sets and the experiment that produced the chosen item is the winner. It is possible for both experiments to produce the same (or similar) set of items and if the winning item was in both experiments, then its a draw. I can keep track of the likelihood of overlap for all pairs of experiments if it makes sense that experiments producing identical item sets offer little information and experiments producing wildly different items contribute significantly. My initial idea was to choose the experiment pair that maximizes $$ \text{UCB}(E_i) \text{UCB}(E_j) \hat{diff}(E_i, E_j) $$ but I have trouble justifying this or finding any bounds. Has anyone seen something similar or have any idea how I should proceed?


title: Are confidence intervals open or closed intervals?

text: I have a question about confidence intervals. > In general, are confidence intervals open or closed?


title: Language of Stata's ADO files

text: What language Stata ADO files are written in? I am looking at svy.ado and curious on what language it is? Is it something built by the Stata company?


title: How to estimate age of a person using their given name?

text: How we could infer the date of birth of a given name? How should I model this problem? Should I infer a single value describing the age itself?


title: Obtaining final classifier in RTextTools

text: I have a question concerning the RTextTools package. When I was reading textbooks and wikipedia I imagined that text classification is all about defining a classificator. In RTextTools you define learning and testing documents in one command. For example               data <- data[sample(1:3100,size=100,replace=FALSE),]     matrix <- create_matrix(cbind(data$Title,data$Subject), language="english",     removeNumbers=TRUE, stemWords=FALSE, weighting=weightTfIdf)     corpus <- create_corpus(matrix,data$Topic.Code,trainSize=1:75, testSize=76:100,virgin=FALSE)     models <- train_models(corpus, algorithms=c("MAXENT","SVM"))      We create a corpus from one dataset, where like first 75 documents is for learning and the rest is for training. But when it comes to applying text classification theories for practical use, you want already predefined text classificator to be applied for different datasets. For now the only solution i managed to find is to train with the same documents and then use different test data. But it sounds very unconfortable when it comes to large texts. So my question is how to create predefined classificator and then apply it to different datasets with R? It would be preferable to use algorithms from RTextTools.


title: Adjusting for tilt of the earth

text: I have rewritten the old question (below) to hopefully make things a bit clearer. Basically I think that the temperature of the earth should be normally distributed but is not due to the ‘seasonal tilt’ and curvature of the earth’s surface. I have added some histograms for different ‘latitude bands’ of the temperatures and rainfall in mm here (I have included rainfall here but please ignore it for now). For example, [-10 ... 10] corresponds to latitude band (-)10 degrees south to 10 degrees north inclusively. Of course temperature depends on latitude but this is due to the fact that it actually depends on seasonal tilt and curvature IMHO. At the end of the day I am after some transformation which removes the effect of seasonal tilt/curvature to arrive at some global normal distribution for all four seasons. Does this make sense? Thanks. Christian Old question: I have some data of the earth's surface (temperature + rainfall in mm): ![enter image description here](http://i.stack.imgur.com/A9wlr.jpg) Clearly the earth's tilt affects the modality/normality of, for example, the temperature distributions. I am just wondering whether there is a way to adjust for this to make the (combined?) data more normal/less modal? I am not a statistician ... so not sure whether this is possible? Thanks in advance. Christian


title: Transformation for stabilizing variance in time series

text: In most forecasting packages there are two transformations easily available for stabilizing variance: the square root and the logarithnic. Is there any procedure (e.g. test) that can tell the user which one of the two transformations to choose? Finally, the variance stabilization transofrmation should be apllied to time series data (if they exhibit multiplicative variance) no matter what is the model to be used (e.g. it should be used in ARIMA, exponential smoohting, unobserverd components, classical decomposition etc if the time series exhibit multiplicative variance)? Thnaks in advance, Andreas


title: The logic behind MSE

text: Why is it customary to use mean square error rather than mean absolute error (absolute value of target - output)? Specifically in neural networks - does it have anything to do with its convenient derivative?


title: How to make a randomForest algorithm cost-sensitive?

text: Having used randomForest in R to produce a fairly successful classifier is there any way to emphasise sensitivity over specificity, for example, if the cost of missing a disease is much greater than diagnosing a false positive?