title: Generating non-normal correlated data for mediation analysis?

text: For a simulation study I have to generate data according to a trivariate mediation model (Baron & Kennys causal step approach). Assuming $X$, $M$, and $Y$ are continuous variables of sample size $n$, the following regression equations constitute the data generating process: $M = \alpha + aX + \epsilon$ $Y = \alpha + cX + bM + \epsilon$ Additionally, $X$ is not normally distributed, but is sampled from a distribution with skewness. I chose log-normal distribution in my current R implementation (see below). Sample size is also highly variable, starting at $n = 15$ and going up to $n = 15.000$. Here is my admittedly really primitive and inefficient code: Sample from `rlnorm()` until a random sample with desired skewness is found (for clarification: `tol` refers to a tolerance parameter set to $.01$).               library(moments)     s <- 1 # desired skewness     tol <- .01     n <- 20     a <- b <- c <- .59     skewness <- sqrt(log(-1+2^(1/3)/(2+s^2+sqrt(4*s^2+s^4))^(1/3)+(2+s^2+sqrt(4*s^2+s^4))^(1/3)/2^(1/3)))          repeat{X <- rlnorm(n=n, meanlog=0, sdlog=skewness)         if (skewness(X) > (skewness - tol) && skewness(X) < (skewness + tol)) {break}     }       Then calculate $M$ and $Y$ from $X$ in a similar fashion and according to the regression equations above (where $a$, $b$ and $c$ are the correlations between the variables):               repeat{         M <- a * X + rnorm(n)         if (cor(X, M) > (a - tol) && cor(X, M) < (a + tol)) {break}     }          repeat{         Y <- c * X + b * M + rnorm(n)         if (cor(Y, X) > (c - tol) && cor(Y, X) < (c + tol) &&             cor(Y, M) > (b - tol) && cor(Y, M) < (b + tol)) {break}     }      Again, this code is woefully ineffective and is only applicable for very low $n$, as the computation time is pretty low for small samples. Once higher sample sizes are reached (about $n=200$), it stops working altogether. So here are my actual questions:   * Is there a way to generate intercorrelated data ($X$, $M$ and $Y$) when $X$ isn't normally distributed and therefore `mvrnorm` from the `MASS` package can't be used? Also the correlation between the variables has to be either an exact value, or within a certain range ($\pm$ tolerance).    * Assuming it can be done theoretically (which I am quite sure of), how can it be implemented in `R`?


title: Comparing transaction counts between samples

text: I have a treatment group and two control groups for customers and the count of their transactions:               Treatment Group     Customer | Transaction Count     -----------------------------     1        | 0     2        | 2     3        | 3     4        | 0 ...          Control Group 1     Customer | Transaction Count     -----------------------------     1        | 1     2        | 2     3        | 0     4        | 3 ...      Each group has approximately 50 customers, ranging from 0-20 transactions, with low numbers (0-5) being most common. I want to test the hypothesis that customers in the treatment group have the same number of transactions when compared to customers in either control group. I keep coming back to the Kruskal-Wallis test, since the counts are definitely not distributed normally. However, I'm a bit worried since there would be a high number of ties (for example, lots of zeros), and my understanding is that the rankings should be clearly ordered. Is there a more appropriate hypothesis test for this question, or should I stick with the Kruskal-Wallis test?


title: Why, when I scale my data set, glmnet gives error?

text: I'm using glmnet for building the regression models. My data are already log- transformed. when I scale my data set (zero mean, and SD=1), I get the following error:               Error in elnet(x, is.sparse, ix, jx, y, weights, offset, type.gaussian,  :        NA/NaN/Inf in foreign function call (arg 5)      But without scaling data, glmnet run perfectly. would someone help me to understand why this happens?


title: Confidence bands for QQ line

text: This question doesn't specifically pertain to `R`, but I chose to use `R` to illustrate it. Consider the code for producing confidence bands round a (normal) qq-line:               library(car)     library(MASS)     b0<-lm(deaths~.,data=road)     qqPlot(b0$resid,pch=16,line="robust")      I'm looking for an explanation of (or alternative a link to a paper/online document explaining) how these confidence bands are constructed (I have seen a reference to Fox 2002 in the help files of R but sadly I don't have this book handy). My question will be made more precise with an example. Here's how `R` computes these particular CI's (I've shortened/simplified the code used in `car::qqPlot`)               x<-b0$resid     good<-!is.na(x)     ord<-order(x[good])     ord.x<-x[good][ord]     n<-length(ord.x)     P<-ppoints(n)     z<-qnorm(P)     plot(z,ord.x,type="n")     coef<-coef(rlm(ord.x~z))     a<-coef[1]     b<-coef[2]     abline(a,b,col="red",lwd=2)     conf<-0.95     zz<-qnorm(1-(1-conf)/2)     SE<-(b/dnorm(z))*sqrt(P*(1-P)/n)     #[WHY?]     fit.value<-a+b*z     upper<-fit.value+zz*SE     lower<-fit.value-zz*SE     lines(z,upper,lty=2,lwd=2,col="red")     lines(z,lower,lty=2,lwd=2,col="red")      The question is: what is the justification for the formula used to compute these SE (e.g. the line `SE<-(b/dnorm(z))*sqrt(P*(1-P)/n)`). FWIW this formula is very different from the formula of the usual confidence bands used in linear regression


title: combining classifiers trained on different data

text: Could you please recommend methods for combining classifiers, which were trained on different patients? For example, there are 10 patients and for each of them I trained a binary classifier. Now I want to combine the classifiers I already have by taking weighted average of their predictions, so my question is how to learn these weights? It resembles stacking a bit (like here), but I'm not sure. I would be grateful for any other proposals.


title: Mixed modelling with lmer

text: I've read multiple posts on model fitting using `lmer()`, but I saw multiple versions of using the error term which made it confusing for me. I have a factorial design with 2×2×2 with one between and two within factors. I do an eyetracking task in which participants are in one of two sound conditions (`cond`). Each participant does the same task in which they see/don't see prime pictures (`cat`). I estimate whether people look more at matching vs. non-matching targets (`means_side`). My basic model is (`cat`= prime vs. control):               fm1 <- lmer(Attention ~ cond*cat*means_side + (1|pp) , data)      Since I saw a lot of variation in the different category pictures I used (four different prime pictures and a control picture, coded as `category`) I want to allow this variable to vary in my model, so I used:               fm2 <- lmer(Attention ~ cond*cat*means_side + (1|pp) + (1|category) , data)      OR               fm3 <- lmer(Attention ~ cond*cat*means_side + (1|pp) + (category|cat) , data)      I am not sure which one to choose, or whether I simply have to use category instead of `cat` in the model specifications. Anybody have an idea?


title: Adding mean and standard deviation to each other

text: I make a certain measurement and I would like to arrive at a single descriptive number that would give me an idea of both the average value of the measurement, as well as the deviation from the mean. In fact, the deviation (general extent - how far do the values go?) interests me more than the mean, but I am aware that interpreting the standard deviation value without the mean will lead me astray. I am considering adding mean and standard deviation together (mean + stdiv = new measure) - how horridly wrong will that be from the statistical point of view? As requested, here are the details: my goal is to measure the extent of saturation in the hidden layer outputs of a neural network. Hidden layer consists of hidden units, where every unit receives weighted real-valued inputs and produces a real-valued output within a certain range (typically [0,1]). An individual output is produced for every pattern in a dataset. In other words, I have a data set (outputs of the hidden units for all patterns in the data set, recorded n times, where n is the number of algorithm iterations), and I would like to see how the average activation value (hidden unit output) changes from iteration to iteration, as well as how much it fluctuates.


title: Issue on prediction with FinalModel of RandomForest in R using the CARET package

text: I use the caret package for training a randomForest object with 10x10CV.               library(caret)     tc <- trainControl("repeatedcv", number=10, repeats=10, classProbs=TRUE, savePred=T)      RFFit <- train(Defect ~., data=trainingSet, method="rf", trControl=tc, preProc=c("center", "scale"))      after that, i test the randomForest on a testSet (new data)               RF.testSet$Prediction <- predict(RFFit, newdata=testSet)      the confusion Matrix shows me, that the model is isnt that bad.               confusionMatrix(data=RF.testSet$Prediction, RF.testSet$Defect)                   Reference         Prediction   0   1                  0 886 179                  1  53 126                  Accuracy : 0.8135                            95% CI : (0.7907, 0.8348)     No Information Rate : 0.7548               P-Value [Acc > NIR] : 4.369e-07                               Kappa : 0.4145       i now want to test the $finalModel and i think it should give me the same result, but somehow i receive               > RF.testSet$Prediction <- predict(RFFit$finalModel, newdata=RF.testSet)     >  confusionMatrix(data=RF.testSet$Prediction, RF.testSet$Defect)     Confusion Matrix and Statistics                    Reference     Prediction   0   1              0 323  66              1 616 239                         Accuracy : 0.4518                                95% CI : (0.4239, 0.4799)         No Information Rate : 0.7548                   P-Value [Acc > NIR] : 1                                           Kappa : 0.0793       what am i missing? edit @topepo : i also learned a another randomForest without the preProcessed option and got another result:               RFFit2 <- train(Defect ~., data=trainingSet, method="rf", trControl=tc)     testSet$Prediction2 <- predict(RFFit2, newdata=testSet)     confusionMatrix(data=testSet$Prediction2, testSet$Defect)          Confusion Matrix and Statistics                    Reference     Prediction   0   1              0 878 174              1  61 131                         Accuracy : 0.8111                                95% CI : (0.7882, 0.8325)         No Information Rate : 0.7548                   P-Value [Acc > NIR] : 1.252e-06                                   Kappa : 0.4167


title: Use of autoregressive metric for ARIMA clustering and analysis

text: I wonder if anyone has put into use the autoregressive metric for ARIMA clustering proposed by Corduas and Piccolo (2008). The authors define the distance autoregressive metric between two processes $X, Y$ as: $$\text{d}(X,Y)=\sum_{j=1}^m \left| (\pi_j,x-\pi_j,y)^2\right|\qquad (1)$$ for $j=1$ to $m$. So in standard ARIMA notation, the squared distance of two $X,Y$ processes as: $$\text{d}(X,Y)^2=(\phi x-\theta x)^2/(1-\theta x^2)+(\phi y-\theta y)^2/(1-\theta y^2)-2(\phi y-\theta y)^2/(1-\theta y\theta x)\qquad (2)$$ So I have a few questions.   1. Does $j$ in (1) refer to the number $m$ ACF and PACF values taken for $N$ lags? If that's the case why do we have to fit in the time series a particular ARIMA model beforehands?   2. For an ARIMA(0,0,1) model the $\phi$ parameter in (2) is 1?    3. For an ARIMA(1,0,0) model the $\theta$ parameter in (2) is 0?   4. In case of a very large value of $\theta x$ or $\theta y$ (close to 1) we get a "ridiculously" large value for squared distance (as described by (2)). Does this have a qualitative meaning? I'm writing the necessary routines in R (for distance matrix calculation). I'm pretty confident there is nothing implemented yet but you never know for sure. However if anyone knows something more or is interested in this field I'll gladly share my experience.


title: How to compare weighted multivariate linear models?

text: I've got a set of multivariate regression models, with weights, that I'm trying to compare in R. Looks like:               f0 <- lm(cbind(Y1,Y2,Y3,Y4) ~ co1 + co2 + co3 + co4 + x1, dat, weight=wt)     f1 <- lm(cbind(Y1,Y2,Y3,Y4) ~ co1 + co2 + co3 + co4 + x1 + x2, dat, weight=wt)     f2 <- lm(cbind(Y1,Y2,Y3,Y4) ~ co1 + co2 + co3 + co4 + x1 + poly(x2,2), dat, weight=wt)      I get perfectly reasonable and interpretable fits, but I'd like to be able to answer a question like, "overall, is there an effect of x2 on the dependent measures?" In other contexts, with nested models, I've done something like:               anova(f0, f1, f2)      and used a Chi-square test. But in this case, I get this error:               Error in SSD.mlm(object) : 'mlm' objects with weights are not supported      So, what alternatives do I have to compare these models? Thanks!


title: Issues when data is unavailable when sampling a finite population without replacement

text: I am specifying a sample of weekly work records to be entered into a spreadsheet for an estimation of means on two variables. There are about 250 employees, I am only interested in weeks between a specific date and the present, the total universe of possible work weeks is 28,538. Because we wanted to assure that the sample was representative by employee and calendar quarter, the population was stratified by calendar quarter and employee, and random sampling with replacement was used within the strata to define a sample of 2190 employee weeks. I recently learned that the data for the first year of the time period we're interested in is unavailable. The unavailable population data includes 50 weeks out of 341 weeks in the period of interest, and also includes 23 out of the 250 employees. Overall, of the 2190 lines specified in the sample as designed, about 175 lines are unavailable. A biased sample representing one employee has already been drawn from the population, the obtained sample variance on the statistic of interest was 115.37, n=276. I realize that the unavailability of data for a certain time period systematically biases the sample's representativeness across calendar quarters, and that bias is introduced by the fact that 23 employees are not represented at all. My question is, what other concerns are introduced by this limitation in the availability of data from which to sample? My statistics experience comes from laboratory experimental design; I now find myself charged with questions regarding sampling and power in a naturalistic environment, any insight would be most helpful. Thanks.


title: Can every retrospective cohort study be transformed into a case control study?

text: A retrospective cohort study basically is classifying people based on exposure status. The disease and exposure have already occurred. So we can also classify people based on cases and look at whether they were exposed or not. Thus can we change a retrospective cohort study into a case-control study?


title: Discrete event simulation with 2D data

text: I want to write a discrete event simulation. I downloaded bus traces (long. and lati. (x and y coordinates)) data from crowded group at Darmouth. When I run statistics tests, x and y data are correlated, so I cannot consider them as two separate distributions. Traditionally, we have a probability distribution, and we generate random numbers according to this distribution, this time I have 2D distribution. I am not clear about fitting how to fit this data to a distribution and generate random numbers to continue other steps in my DES. Anyone have ideas? I would really appreciate.


title: Extending the birthday paradox to more than 2 people

text: In the traditional birthday paradox the question is "what are the chances that two or more people in a group of n people share a birthday". I'm stuck on a problem which is an extension of this. Instead of knowing the probability that two people share a birthday I need to extend the question to know what is the probability that x or more people share a birthday. With x=2 you can do this by calculating the probability that no two people share a birthday and subtract that from 1, but I don't think I can extend this logic to larger numbers of x. To further complicate this I also need a solution which will work for very large numbers for n (millions) and x (thousands).


title: Regarding boosting, bagging and bootstrapping

text: How to understand the relationships, comparative advantages, and comparative disadvantages of boosting, bootstrapping and bagging in terms of their respective applications in data mining.


title: How to test for normality of growth disturbances in chemo treatment?

text: I'm a med student, conducting a retrospective analysis of weight/growth disturbances during chemo treatment.   **I wonder, if I should:**   * assume, that growth is a variable normally distributed across the population   * test for (or visually assess) the normality once (when?)   * test for (or visually assess) normality separately, for each sample/each test performed **Why is that a problem for me:**   * there are many samples and subsamples (gender, treatment protocol, time point of measurment)   * samples have different means (supposedly effect of the treatment or small sample size)   * sample sizes are different (ranging ~5-60)   * treatment can hypothetically change the distribution of the variable (however I doubt it, since I've already generated some histograms and sometimes (bigger sample size) values look normally distributed)   * each analysis will exclue different cases (database is not complete, I must perform pairwise deletion)   * tests for normality are useless with small samples? **Here are some examples.** : http://imgur.com/a/zb2hM     * X axis is an SDS for weight(it must compared indirectly, since patients are at different age; population's SD and mean are known for the specified age)   * Samples differ with sample size, case's gender, treatment etc. ### Questions   * Should I test for normality or assume normality?   * Should I test for normality separately in each sample or just once?   * How should I test for normality?


title: Why isn't the Dantzig selector popular in applied statistics?

text: Lasso-like methods have become pretty common in applied statistics but the Dantzig selector remains unpopular despite having great properties (minimax optimality). Why hasn't it become more popular?


title: How to standardise cases within groups to perform ranking across groups and over time?

text: I am trying to create a method to rank ships according to their fuel burn rates. There are two different classes of ships, with different burn rates, so I am trying to create a mean and standard deviation for each class, then rank the ships according to how many standard deviations they are below the mean of their class. I will need to do this each month and each quarter. There are three ships in one category and 7 in the other. In the category with three ships, there are also 9 ships in another area that I do not wish to rank, but have data for. In the category with 7 ships, there are 23 other ships I have data for, but do not wish to include in the ranking. Questions:   1. Is this the correct way to rank ships in two different classes? Any suggestions?    2. Should the mean and standard deviation remain constant, or change every month / quarter?    3. should the standard deviation and mean calculations include the other ships that I will not be ranking?    4. I have data that lists the ships' average burn rates over a three year period, but these averages are not all based on the same amount of data (some ships were underway all year while some were underway only for a month).   Can I use this data to create the standard deviation and mean, or do I need the monthly / quarterly data?


title: How to rank order performance of four systems using aggregates of pairwise assessments of system performance?

text: I am currently designing an experiment to evaluate 4 iterations of a system. As it is rather difficult for a human to judge the output of the system by assigning it a score or several scores, I plan on giving the following task: Give the human evaluator the outputs of two systems on the same input and ask him to decide which system's output is better in certain regards. My final desired result is to have a ranked list of system iterations and/or pairwise comparisons stating which system iteration is better. How should I proceed with the experimental design? What significance test should I use? **Update** : It seems that sign test could be appropriate to test whether one system is truly better than the other. **Update2** : Does obtaining the ranking of 4 system iterations necessitate performing _m_ trials for each pair or could transitivity be used, as in if A > B and B > C, then A > C, to reduce the number of required trials?


title: Selecting cases for analysis based on multiply imputed values

text: I have a large public dataset describing traffic crashes. I am interested in events occurring in a certain type of street intersection only. The dataset has a lot of missing values, and I plan to use multiple imputation (I'm using SAS). My plan has been to use the full dataset to do the MI, then select the cases I want to analyze, then do the analysis. Three variables provide information for selecting the cases in which I am interested. But here is the problem: two of these variables have missing values, so there will likely be uncertainty as to whether these cases belong in my analysis. I have looked and looked for a description or example of this online but found nothing. I'm not sure at this point if I'm making some kind of naive student error (my first time doing MI), or whether people don't really talk about this. It's easy to think of ad hoc ways of proceeding (e.g., impute only from cases that are complete on the selection variables--which seems almost as bad as just analyzing only complete cases--or selecting cases in which some threshold proportion of imputed values for a case meet the selection criterion). However, I'd like to do this "by the book," to the degree there is a book. Does anyone have thoughts about how I should select cases for analysis? (A few details: Most of the variables are categorical, and the missing pattern is decidedly nonmonotonic. I have been planning to convert the categorical variables into dummy binaries and use MCMC for the imputation. There is a monotonic-ish missing pattern among the case selection variables, so I could do MI in two steps, using the discriminant analysis method, which at least gives discrete values within a dataset, for the selection variables, but this still does not solve the likely problem of variability of values across MI datasets in the selection variables. As you can see, without a case selection method in mind ahead of time, I'm just grasping for things that make the problem seem simpler.)


title: Mahalanobis distance and percentage of the distribution represented

text: In a one dimensional normal distribution, it is really handy to know that 68% of the data are within one standard deviation, 95% lies within two standard deviations, etc. My question is about the higher dimensional version of this. Can anyone give me rules of thumb similar to the ones we use so often with standard devs, i.e., in dimension `N`, what percentage of a normal distribution is at Mahalanobis distance 1, 2, 3.... I kind of assume someone has come up with handy rules for this, even if they are only good approximations. The point is to know at what Mahalanobis distance I will have excluded x% of the data. **Notes** I worked out the distribution of the Mahalanobis distance for a normal distribution in N dimensional space to be $$ c \cdot r^{N-1} \cdot \exp\Big(\frac{-r^2}{2}\Big) $$ where $c$ is a constant. It is easy to check that the most common distance for a point will be $\sqrt{N-1}$. This means $N > 1$ is quite different than the one dimensional case because the the most frequent distance is not distance zero. (Yes, I know one can integrate this by hand using successive integration by parts until one gets an answer in terms of $erf$ ... but (1) that sounds painful for unspecified $N$, and (2) it has got to already be known)


title: Exponentiated logistic regression coefficient different than odds ratio

text: As I understand it, the exponentiated beta value from a logistic regression is the odds ratio of that variable for the dependent variable of interest. However, the value does not match the manually calculated odds ratio. My model is predicting stunting (a measure of malnutrition) using, amongst other indicators, insurance.               // Odds ratio from LR, being done in stata     logit stunting insurance age ... etc.      or_insurance = exp(beta_value_insurance)          // Odds ratio, manually calculated     odds_stunted_insured = num_stunted_ins/num_not_stunted_ins     odds_stunted_unins = num_stunted_unins/num_not_stunted_unins     odds_ratio = odds_stunted_ins/odds_stunted_unins      What is the conceptual reason for these values being different? Controlling for other factors in the regression? Just want to be able to explain the discrepancy.


title: AIC & BIC number interpretation

text: I am looking for examples of how to interpret AIC (Akaike information criterion) and BIC (Bayesian information criterion) estimates. Can negative difference between BICs be interpreted as the posterior odds of one model over the other? How can I put this into words? For example the BIC = -2 may imply that the odds of the better model over the other model are approximately $e^2= 7.4$? Any basic advice is appreciated by this neophyte.


title: Is a different CV arrangement the same as a validation set?

text: I have a smallish dataset ~ 1500 rows X 500 columns. I've been using a standard 5 fold CV setup where row 1 = CV set1, row2 = CV set2, ... row 6 = CV set1,etc. I'm at the point where I'm trying to do some feature reduction work/parameter optimization and concerned about over-fitting by using the same CV setup to find parameters and then get error measures. My question is this: if I create a separate CV indexing setup where row 1 = CV set 5, row 2 = CV set 3,..., row6 = CV set 1, etc. does this function similarly to having a separate validation set or do I need to explicitly have a 50/25/25 train/validation/test setup?


title: Fitting models to empirical data: do fits for random samples converge on the true fit?

text: I'm trying to fit a mathematical/theoretical model to empirical data, but the dataset is impractically big to fit all at once. (Specifically: I'm fitting a power law model using the methods of Clauset, Shalizi, & Newman (2007) - but I think the specific model doesn't matter for this question.) Is it true that, if I fit multiple 1% random subsets of the data to the best- fit models for those subsets, that the mean of the model parameters will converge on the best-fit models for the full dataset? This seems to make sense as a consequence of the Central Limit Theorem, but lots of things seem to make sense that aren't true. Thanks!


title: Plotting selected explanatory variables in multiple linear regression with interaction terms

text: I have a question about plotting multiple linear regression results. I would like two graphs, one that shows the estimated relationship between d15N and avomc, and another that shows the relationship between d15N and N. My questions are: What are the equations giving the fitted values for these two relationships? Does it make sense to plot such a thing?               Call:     lm(formula = d15N ~ avomc + N + fPlant + fGrazing + fblock +      fPlant:fGrazing, data = D_data)          Residuals:      Min       1Q   Median       3Q      Max      -2.50874 -0.47156  0.02093  0.56252  2.04284           Coefficients:                                Estimate Std. Error t value Pr(>|t|)         (Intercept)                   -0.428674   0.848159  -0.505 0.614178         avomc                        -17.400385   6.620888  -2.628 0.009689 **      N                              0.541445   0.192159   2.818 0.005645 **      fPlantPcrantzii               -1.260810   0.293716  -4.293 3.56e-05 ***     fPlantSvirgaurea              -2.170915   0.277706  -7.817 2.15e-12 ***     fGrazingout                   -0.002856   0.263470  -0.011 0.991368         fblock2                        1.203011   0.313434   3.838 0.000198 ***     fblock3                        1.053795   0.319100   3.302 0.001258 **      fblock4                       -0.153207   0.294686  -0.520 0.604076         fblock5                        0.379873   0.292868   1.297 0.197051         fblock6                       -1.015164   0.292780  -3.467 0.000727 ***     fPlantPcrantzii:fGrazingout   -0.875229   0.398690  -2.195 0.030037 *       fPlantSvirgaurea:fGrazingout  -0.040601   0.361892  -0.112 0.910857       Thanks for your time and help!


title: Two sample tests for time-series segments

text: I have a data set consisting of a number of time-series segments, i.e each segment contains n number of contiguous samples from a time-series. Are there any good statistical tests to determine how likely two segments come from the same probability distribution?


title: Testing for useful variables in a "net lift model"

text: I am often involved in modeling the Net lift, aka Uplift, aka incremental response of direct marketing campaigns. In a nutshell, this approach looks to model and thus select for marketing those individuals who require promotion in order to take the desired action (e.g. order a product). Given prospect A and B, if A will buy if we send a direct mail letter and B will buy anyway, we target A. Think of "swing voters" in political elections - you put resources against the persuadable undecided. I am looking for affirmation or suggestions on the following approach to help understand which variables are important. Here is an example of a binary predictor variable. Lets say that 20,000 direct mail letters were sent (10,000 to those with predictor1 =0 and 10,000 with predictor1 =1) and 2,000 control customers were selected and held out to not receive the letter. In this case, the difference in response rate (treated group - control) is 1.5% for predictor1 =0 and -0.3% for predictor1=1. We would conclude that those with predictor1=0 are the better candidates to send a letter to - as the letter actually decreased the point estimate of response for those with a value of 1 for predictor1. The difference in the odds ratio is 1.413 (treated - control). ![enter image description here](http://i.stack.imgur.com/roCGQ.jpg) Here is an example where predictor 1 is not useful. ![enter image description here](http://i.stack.imgur.com/xAi8l.jpg) Is it accurate to test the value of a given predictor on this difference in response rates (for any type of predictor, nit just a simple binary one) using a standard logistic regression? For example,the first scenario shown above would result in a significant interaction variable of treated*predictor1:               dat<-data.frame(treated=c(1,1,0,0),predictor=c(0,1,0,1), responded=c(250,122,10,15), notresponded = c(9750,9879,990,985))     mod<-glm(cbind(responded,notresponded)~predictor*treated, data=dat,family=binomial)     summary(mod)           Coefficients:                       Estimate Std. Error z value Pr(>|z|)         (Intercept)        -4.5951     0.3178 -14.458  < 2e-16 ***     predictor           0.4105     0.4107   1.000  0.31754         treated             0.9316     0.3242   2.873  0.00406 **      predictor:treated  -1.1411     0.4255  -2.682  0.00733 **       While the second scenario does not:               dat<-data.frame(treated=c(1,1,0,0),predictor=c(0,1,0,1), responded=c(130,122,10,15), notresponded = c(9870,9879,990,985))     mod<-glm(cbind(responded,notresponded)~predictor*treated, data=dat,family=binomial)     summary(mod)               Coefficients:                       Estimate Std. Error z value Pr(>|z|)         (Intercept)        -4.5951     0.3178 -14.458   <2e-16 ***     predictor           0.4105     0.4107   1.000    0.318         treated             0.2654     0.3299   0.805    0.421         predictor:treated  -0.4750     0.4299  -1.105    0.269        **Question:** Issues with multiple testing aside, does this not represent a viable method for testing variables for net lift modeling? **Should I exclude the main effects and only have the interaction term?** This is a principle question and speaks to the interpretation of an interaction - given that I am looking for "significant" variables which could be used to select customers to market to in order to maximize the difference between treated and control. Can one also use this framework to test for deeper interactions (treated*predictor1*predictor2) to see if we need to really look at predictor1 and predictor2 combinations?


title: Fewer variables have higher R-squared value in logistic regression

text: I am testing out 3 modeling approaches for malnutrition in children. Theoretically, distal determinants (education,poverty) operate through proximal determinants (water, sanitation) in determining malnutrition rates. The three logistic models, where stunting is a binary indicator for malnutrition, are:               // Proximal determinants only: both binary indicators     stunting ~ water + sanitation          // Distal determinants only: both categorical indicators     stunting ~ i.education + i.poverty          // Both proximal and distal determinants     stunting ~ water + sanitation + i.education + i.poverty      I am surprised to find that the r-squared value of the second model is higher than the third model, as calculated by the correlation between the predicted and actual values (stata):               predict predicted, xb     corr predicted stunting     local rsq = r(rho)      While I expected the strength of the relationship and statistical significance of the more proximal causes to decrease (as they were soaked up by the distal causes), I expected the combined model to have higher explanatory power (as measured by r-squared). Does anyone have any explanation as to why the second model has the most explanatory power? Let me know if I can provide additional information for answering this question.


title: Johanson test conditions and Breusch-Godfrey LM test

text: I am a student from Belgium and I am making a thesis about the relationship between credit aggregates and property prices. I examine the Granger causality between the two variables and I also do some conintegration tests. I have a question about the latter. What are the conditions for doing a Johansen cointegration test? Do the residuals need to be tested for serial correlation before you can de the test? Because sometimes I just don't find a model where this is the case. Or do the residuals don't need to be uncorrelated? Are there other conditions regarding the Johansen test? When I test for residual autocorrelation I use the Breusch-Godfrey LM test. What is the lag-order that I need to choose for this test? I have 158 observations in my time series. If I do the varselection in levels I always used a maximum lag of 12 and for the Breusch-Godfrey LM test I always used 6. Is this correct? If so: where do I find references to support these conclusions?


title: t-test for more than 30 samples?

text: I'm conducting a study that requires to compare the means of dependent groups and independent groups. I thought that t-test would be the best test. But I'm having doubts since I know that t-test is used when the samples are less than 30. So can I still use t-test even if my I used 150 samples for both groups? Thanks


title: What's the name of this quantile curve estimation technique?

text: I'm still working with the same dataset as in a previous question here. Quantile regression has taken me some way, and now I'm going further. The underlying relationship between these x and y variables has no reason to be linear. Monotonic increasing (plus noise), but not linear, and not necessarily _strictly_ increasing - it can plateau in sections. So, I've applied a new algorithm, to try to pick out the underlying curve. The thing is, I've just developed this by playing with the data. I'm pretty sure I've reinvented someone's wheel, as it's not a particularly deep algorithm, and there's little new under the sun. (and I'm not yet sure if the wheel is square or round). Here's the curve I've derived from some of the data from the previous question. ![](http://i.stack.imgur.com/btVdw.png) Here's how I did it. I've got a set of several thousand paired observations $(x,y)$. First I create a set of points $X$ by dividing the x-axis up across N evenly- spaced points, with $x_{k+1}=x_k+I$, i.e. interval width $I$ between them. For each value $x \in X$, I take the subset of observations which have an x value in the range $[x-mI,x+mI]$, for some value of smoothing parameter $m$: in the case illustrated, $m=5$. Then, I take a particular quantile $q$ of the y-values of those observations - e.g. the median ($q=0.5$) - and that's the y-value of the estimated curve at point $x$. For any $m>0.5$, the subsets of observations for two neighbouring points $x_k, x_{k+1}$, overlap. So, there are three parameters to choose in any curve construction: $I,m,q$. N is calculated just by dividing the range of observed $x$ by $I$ and rounding up. **What is the name of the technique that I've reinvented? And is it a known Good Thing, or a Bad Thing?** It seems useful, and I'd like to write up my work to date; so I'd like to get the name of the algorithm right, and check the literature to see if this approach is sane and (ideally) robust. This isn't spline-smoothing (no knots), and it doesn't seem to fit into the class of general additive models either. It does seem to handle the vicious heterskedasticity in my data pretty well. It's not quantile regression, because it's definitely not fitting a straight line (and does seem somewhat superior to the quantile regression I did previously). I don't yet have a rigorous basis for optimising the smoothing parameters; I do have some faith that a meaningful metric could be defined.


title: Limiting distributions identified as functions of Brownian motion or stochastic integrals

text: I am teaching a stochastic processes course to MA stat students, and to stay on topic I would like some examples of limiting distributions in stat that are identified as functions of brownian motion or stochastic integrals. I have thought of kolmogorov-smirnov & dickey-fuller, neither of which uses stochastic integration. Examples have to be short, especially, not require lots of stats backround.


title: Standard measurements in excel

text: I have a set of data in MS Excel 2013 and I want to do very common calculations on it![screenshot](http://i.stack.imgur.com/9wITY.png) I got a bit confused with the different formulae available and am not quite sure if I chose the correct ones. Could you take a look at it? Thank you very much!


title: Two-sample test of proportions

text: In my research, I aim to compare websites and need to identify if the websites differ significantly from each other. I have divided the website into two categories. I have sampled 45 websites in the first category and 30 websites in the the other category. I need to identify if the first set of websites differ significantly from the second set. My data is largely nominal, as I just count the presence or absence of certain elements on the webpages. So, during my data collection stage, I counted the Presence or Absence of certain elements on these website, and calculated the the percentage of occurrence of these elements, for the two sets of websites. When I had the percentage of each elements for both the sets, I applied the Two Sample Test of Proportions, and tested my hypotheses (Is one set different from the other in terms of element A, B, C...etc). I had 60 elements, and I applied the test separately 60 times for each element. My results do make sense, but I am just not sure if the test is a valid measure for my research. Can someone please help.


title: Bootstrapping and comparing multiple proportions

text: **EDIT 1: Thanks to @gung for pointing out that if the bag initially had 324 sweets, then as Mother hands them out, there are consecutively fewer amounts she can hand out. To simplify things (because I'm interested in understanding the basics instead of something too advanced), he offered a rewording.** **EDIT 2: This is not homework. I'm just trying to learn some statistics in a fun way** **COMMENT: From all of the kind comments/answers thus far, it occurs to me that I really need to work on defining my questions in the future much better because I didn't realise how much the answer could vary depending on how it's posed (very useful lesson!). I initially thought this was well posed but it's good to learn these things and hopefully improve, thank you.** Let's say a a mother has a bag containing an infinite number of sweets. We keep a track of how many sweets she gives to each of her children named "A", "B", ..., "G" from this bag, and how many she gives out to an _unknown_ number of adults. When the mother comes up to each person, she opens a new, identical bag of sweets, reaches in & grabs some for that individual. At the end of the day, she has given out a total of 324 sweets!               > (DF <- data.frame(A=15, B=4, C=1, D=4, E=44, F=4, G=1, Adults=251, Total=324))             A B C D  E F G Adults Total       15 4 1 4 44 4 1    251   324      **Question 1:** Her children were given a total of 73 sweets out of 324. The adults were give 251 sweets out of 324. They children want to know if the adults were given a statistically significant proportion of the sweets even though their mother claims that everyone had an equal chance of getting the same number of sweets and it was just luck that some got more than others. So basically we want to compare the _entire_ group of children against the _entire_ group of adults, not individuals. **Question 2:** Did one or more of the children get a statistically significant number of more sweets when compared to their siblings (i.e. we ignore the Adults for this question and only consider the children "A", ..., "G")? And if so, which one(s) got a statistically significant proportion? Mother again claims everyone had an equal chance of getting the same number of sweets, and it was just luck that some got more than others. **Additional Information:** I want to use a bootstrapping approach because I seem to get a better feel for probability when using simulation (NB: I am not a statistician, just doing this for fun). Below I will give my approach to the first question which I think is OK but would appreciate some advice on if I'm doing it correctly. I have no idea how to approach the second question because there are multiple proportions involved. * * * My approach to Question One: **Hypotheses:** My null hypothesis is that children and adults had an equal chance of getting the same number of sweets. My alternative hypothesis is that Mother is more generous with giving out sweets to adults than to her children. **Research Question:** How often might we observe a proportion of sweets as small as the children got, if the mother really was handing out sweets fairly? **Methodology:** We can test this directly by creating a hypothetical population that embodies the null hypothesis, with 50% "children" and 50% "adults" and repeatedly drawing random samples of 324 from it, with replacement, recording the results. We then look at how many of these samples have 73 or fewer "children" nodes in them (this is our p-value).               replicates <- 999     size <- 324     runs <- list(replicates)     reference <- 73/size          for(i in 1:replicates){       runs[[i]] <- table(sample(x=c("children", "adult"), size=size, replace=TRUE, prob=c(0.5, 0.5)))      }          bootstraps <- as.data.frame(do.call(rbind, runs))     #   adult children     # 1   166      158     # 2   171      153     # 3   158      166     # 4   151      173     # 5   156      168     # etc.          sum(bootstraps$children <= reference) / replicates     #[1] 0      **Conclusion:** My p-value is effectively zero, which means there is strong evidence that Mother has been handing out sweets unfairly because getting 73 or fewer sweets is very unlikely to have happened by luck alone. As for question 2, I'm drawing a complete blank :(


I need to prove that the limiting distribution of $(n-1)S^2/\sigma^2$ is a Normal, where $S^2$ is the sample variance. However, I have no clue on how to do it. We know that $(n-1)S^2/\sigma^2 \sim \chi^2(n-1)$, however, I'm not even sure this fact can help me... Any help would be appreciated EDIT:I needed to add some info. Sorry. $X_1,...,X_n$ are iid Normal random variables, with mean $\mu$, and variance $\sigma^2$


title: PRESS statistic for ridge regression

text: In ordinary least squares, regressing a target vector $y$ against a set of predictors $X$, the hat matrix is computed as $$H = X (X^tX)^{-1} X^t$$ and the PRESS (predicted residual sum of squares) is calculated by $$SS_P = \sum_i \left( \frac{e_i}{1-h_{ii}}\right)^2$$ where $e_i$ is the $i$th residual and the $h_{ii}$ are the diagonal elements of the hat matrix. In ridge regression with penalty coefficient $\lambda$, the hat matrix is modified to be $$H = X (X^t X + \lambda I)^{-1} X^t$$ Can the PRESS statistic be calculated in the same way, using the modified hat matrix?


title: Determining Sample Size to Show Increase in Disease Prevalence

text: A prof and I are trying to figure out how many samples a study would need to show a certain % increase in disease prevalence, but I'm not sure I'm doing this correctly. Basically we want to know how many samples would be needed to show, with 90% confidence, a 5% increase in brain tumors due to cell phone usage assuming the prevalence of brain tumors is 1/10,000. He doesn't think an adequate number of samples have been taken in previous studies to give this confidence, but I'm not sure how to calculate the numbers myself. My thinking right now is that I should calculate two confidence intervals using the equation CI = p +- z*sqrt(p(1-p)/n) where z=1.96, p is 1/10,000 in one case and 1.05/10,000 in the other case (5% higher), then find the n necessary such that the two CI don't overlap. Again, I'm not sure if this is the way I should be approaching the problem and any guidance would be appreciated.


title: How to determine which variable or combination of the variables are affecting to the predictor variable?

text: I have one dependent variable name as "win ration" of the deal contested and more than 30 independent variables, all are categorical variable name as role of the customer, geo, region, and 27 competencies marks between 1 to 12 for each competency (this is like the performance parameter "1" is high and "12" is low rating) name as comp1,comp2,....,comp27 My question is that which is the best model or predictor model to filter out which all competencies and other variables are really affection the win rate. for this i used beta regression but non of the competencies are coming out significant and when i perform the step wise for dim reduction, this method is not working on beta reg. In this case all competencies are treated as quantitative variable, the reason is that all are bound between 1 to 12 and if i will categorize this all 27 competencies will have 12 categories Please help me to do this analysis


title: Multiple response question with duplicated answers

text: I have a multiple response (categories) question (Q11, ..., Q15). My problem is that the data-entry clerck repeats the answers in the same case. Example:                      Q11 Q12 Q13 Q14 Q15     Case 1 001 003 015 001 022 (001 was type twice)     Case 2 089 032 089 089 014 (089 was typed three times)      I need to find a way to identified the cases with these errors. Thanks,


title: Reproduce linear discriminant analysis projection plot

text: I'm struggling with projection points in linear discriminant analysis (LDA). Many books on multivariate statistical methods illustrate the idea of the LDA with the figure below. ![figure-1](http://i.stack.imgur.com/9k7iT.png) The problem description is as follows. First we need to draw decision boundary, add perpendicular line and than plot projections of data points on it. I wonder how to add projection points to the perpendicular line. Any suggestions / pointers?


title: Chow test or not?

text: I am trying to set up an automatic screen to detect structural breaks in large numbers of time series. The time series are weekly and represent behaviour of customers. I have set up a Chow test. I use the most recent 4 weeks and compare it to the preceding 22. I want to know if their recent behaviour has significantly diverged from the preceding behaviour. My question is this: > Is the Chow Test the most appropriate test for this question? If this is not the most appropriate test how can I determine what test is the most appropriate one?


title: How to do Chow test for constancy of parameters across 2 groups

text: I am a beginner in statistics and poor in mathematics. I am trying to to assess effect of intervention in one state versus another using annual data.               My data are     State 1         State 2      Cases   Deaths  Cases   Deaths     2004    1125    5   2024    254     2005    1213    5   1978    209     2006    1003    4   2294    217     2007    1425    6   2312    249     2008    1172    4   1528    197     2009    1092    3   1683    204     2010    1316    4   2024    218      When I was stuck for the correct statistical procedure, one senior member "irishstat" advised me the following very convincing analysis for which I am ever grateful: X1 is the number of cases and Y is the number of deaths. X2 is the empirically identified point of anomaly; (2009 .. period 6 for State1 and 2004 .. period 1 for State2. Outlier detection led to identifying one anomalous data point for each state reflecting some unknown background variable thus yielding a more robust estimate of the mortality rates. Analysis of State1               State1 Y(T) = -.65649     +[X1(T)][.0046)] CASES +[X2(T)][-1.3608)] PULSE6 I~P00006STATE1 + [A(T)]     Suggesting an unusually low mortality rate for 2009     Analysis of State2     State2 Y(T) = 123.55     +[X1(T)][(+ .0468)] CASES +[X2(T)][(+ 35.7590)] PULSE1 I~P00001STATE2 + [A(T)]     Suggesting an unusually high mortality rate for 2004      This leads to estimating two cleansed data points               STATE    YEAR   Y OBSERVED   Y ADJUSTED     STATE1   2009   3            4.36     STATE2   2004   254          218.24      Replacing these two observed possibly errant values possibly due to some unspecified concomitant factor (“lurking Variable”) one computes a rate of.0046 for STATE1 and .0468 for STATE2. My problem now is how to do The Chow Test for constancy of parameters across groups to check for rejection of the null hypothesis of equal coefficients. I have SPSS v19. Kindly advise me step by step.


title: Estimating variability of unseen factor

text: I'm looking at binomial data where I believe that the probability of the outcome is the product of two independent factors. If you think of it as a two step decision,   1. At the first step, there is a certain probability $p$ that you will have a success, and $1-p$ that you will have a failure. If you have a failure at this step, there is no further opportunity for success.    2. Then, at the second step, there is a certain probability $q$ of success. So, the overall probability of success is $p\times q$. I'm assuming, for now, that $p$ and $q$ vary independently across individuals. Both $p$ and $q$ have some variance across individuals, but not necessarily the same one. Due to the nature of the data, I can estimate the variance of $q$, and the variance of $p\times q$. From this, I would like to estimate the variance of $p$. I just can't wrap my mind around how.


title: allowing for different variances between groups versus allowing for random slopes in lme

text: I have a data set that looks like this: Genotype Condition Trait A 1 0.0007 B 1 0.005 A 2 0.0003 B 2 0.0097 I'm using lme (from nlme package in R) to fit a linear mixed model to the data allowing for random intercepts and random slopes (or not):                   model1 <- lme(trait ~ Condition, random = ~ 1 + Condition|Genotype, data=data, method = "REML")      or                   model2 <- lme(trait ~ Condition, random = ~ 1|Genotype, data=data, method = "REML")      Then I use an anova [anova(model1,model2); note that I always use REML and keep the same fixed effects] to decide if the random slopes are important or not for the fit of the model. When looking at the residuals plot I noticed that the fact that my data has a lot of zeros, and exhibits different variances between conditions, is probably affecting the modelling. So I decided to try fitting a model that allows for different variances by Condition to see if it would improve the fit.                   model3 <- lme(trait ~ Condition, random = ~ 1 + Condition|Genotype, weight=varIdent(form=~1|Condition), data=data, method = "REML")      and                   model4 <- lme(trait ~ Condition, random = ~ 1|Genotype, weight=varIdent(form=~1|Condition), data=data, method = "REML")      When I looked into the new models I verified that model3 is no different from model1, and that model4 is also identical to model1 and model3. The difference is that model4 exhibits better dispersal of residuals for Condition 2.                   > anova(model1,model3)                Model df       AIC       BIC  logLik   Test L.Ratio p-value         model1     1  6 -166.0884 -149.8855 89.0442                                model3     2  7 -164.0884 -145.1850 89.0442 1 vs 2       0       1              > anova(model2,model4)                Model df       AIC       BIC   logLik   Test  L.Ratio p-value         model2     1  4 -131.5512 -120.7493 69.77559                                 model4     2  5 -168.0884 -154.5860 89.04420 1 vs 2 38.53722  <.0001              > anova(model3,model4)               Model df       AIC       BIC  logLik   Test      L.Ratio p-value         model3     1  5 -168.0884 -154.586 89.0442                                     model4     2  7 -164.0884 -145.185 89.0442 1 vs 2 5.684342e-13       1      My question is: is there a conceptual reason why allowing for random slopes per genotype is similar to allowing for different variances between conditions? Is allowing for different variances between conditions somehow allowing for different slopes per genotypes? Or is this just a coincidence for my data? Thank you.


title: Conditional probability and instantiation

text: This questions consists of two parts that are quite similar and concern conditional probability. Firstly, I would like to confirm the calculation of the conditional probability when we know the value of a random variable. Let's assume that we have a conditional matrix which gives us the values for $P(D|A\&B)$ and that B can take the values $(b0, b1)$. If we know that $B = b0$, then if we have to calculate $P(D|A)$, it would just be $P(D|A\&b0)$, right? Let's now assume that we have a conditional probability matrix giving us $P(F|G\&H)$. If we assume that $F$ is independent from $H$, how can we calculate $P(F|G)$ only using $P(F|G\&H)$?


title: Median of medians as robust mean of means?

text: The location and scale of a normally distributed data can be estimated by sampling the data then taking the mean of the sample means and standard deviations, respectively. For non-normal (heavy-tailed) data, is it correct to take the median of the sample medians and IQR/MAD, instead? That is, is it correct to use the median of sample medians as a robust estimator of location similar to the mean of sample means for normal data?


title: Can someone help me with the interpretation of this cox proportional hazards model and how to cross-validate it?

text: The question I'm asking is related to this question here and here. And I apologize for asking so many questions here, as I am thoroughly a stats novice and probably my MD thinking is clouding my ability to interpret my findings. So, because of the input I've received and following some contemplation, I've decided to run a Cox proportional hazards model on this dataset regarding the treatment of patients and how this affects their time of wound healing. As the duration "wound healing" is affected by censoring, following @Glen_b and @Frank Harrell's input, I believe running a normal linear and/or Poisson regression model (as the days are counted and the Poisson's residual qqplots actually have a more linear/normal distribution) isn't quite feasible, because neither of these models actually account for the fact, that when patients are healed, they drop out of the "study" and thus are "censored". Therefore I've decided to explore Cox proportional hazards modelling a bit further. Although not originally planned from a study endpoint perspective (long story), I've decided to chose the need for operative intervention (because the wound fails to heal) as the time to event here (also because I need an event variable in the `coxph` and `Surv` function in R). The study sample size is 4918 patients, of which 575 required operative treatment (11.7%). Now I'm not even sure if this is a high enough event rate as the independent variables are n = 14 and if I go by the "rule of thumb" of 10 events per predictor variable, I'd actually require 1400 events (although I read here, that this rule can sometimes be relaxed), but that's also one of the reasons why I'm posting my question here. I've run a multivariable Cox proportional hazards model, adjusting for those variables that may influence the decision to "treat" a patient and or/wound healing dynamics as well as the treatment the patients received. The output looks as follows:               library(survival)     multicoxph<-coxph(Surv(Timetoheal,Operation=="Yes")~FA     +Age+Woundsurface+Mechanism+Sup+Mid+Deep)     summary(multicoxph)     Call:     coxph(formula = Surv(Timetoheal, Operation == "Yes") ~ FA + Age +          Woundsurface + Mechanism + Sup + Mid + Deep)            n= 4877, number of events= 575         (41 observations deleted due to missingness)                                   coef exp(coef)  se(coef)      z Pr(>|z|)         FAInadequate        -0.225084  0.798449  0.089056 -2.527   0.0115 *       Age                 -0.012107  0.987966  0.002401 -5.042 4.60e-07 ***     Woundsurface         0.042953  1.043889  0.017338  2.477   0.0132 *       Mechanism1          -0.341706  0.710557  0.168422 -2.029   0.0425 *       Mechanism2          -0.169782  0.843848  0.312529 -0.543   0.5870         Mechanism3          -0.421703  0.655929  0.523798 -0.805   0.4208         Mechanism4          -0.131178  0.877062  0.175603 -0.747   0.4551         Mechanism5           0.071292  1.073895  0.313708  0.227   0.8202         Mechanism6          -1.132178  0.322331  0.473741 -2.390   0.0169 *       Mechanism7          -1.216408  0.296292  0.308236 -3.946 7.94e-05 ***     Mechanism8          -0.011187  0.988876  0.162881 -0.069   0.9452         SupYes              -1.476173  0.228511  1.121847 -1.316   0.1882         MidYes              -0.056176  0.945373  1.010415 -0.056   0.9557         DeepYes              1.568644  4.800135  1.002629  1.565   0.1177         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1                              exp(coef) exp(-coef) lower .95 upper .95     FAInadequate           0.7984     1.2524   0.67057    0.9507     Age                    0.9880     1.0122   0.98333    0.9926     Woundsurface           1.0439     0.9580   1.00901    1.0800     Mechanism1             0.7106     1.4073   0.51079    0.9885     Mechanism2             0.8438     1.1850   0.45734    1.5570     Mechanism3             0.6559     1.5246   0.23496    1.8311     Mechanism4             0.8771     1.1402   0.62167    1.2374     Mechanism5             1.0739     0.9312   0.58067    1.9861     Mechanism6             0.3223     3.1024   0.12737    0.8157     Mechanism7             0.2963     3.3750   0.16194    0.5421     Mechanism8             0.9889     1.0112   0.71862    1.3608     SupYes                 0.2285     4.3762   0.02535    2.0598     MidYes                 0.9454     1.0578   0.13048    6.8497     DeepYes                4.8001     0.2083   0.67269   34.2526          Concordance= 0.746  (se = 0.017 )     Rsquare= 0.066   (max possible= 0.769 )     Likelihood ratio test= 331.9  on 14 df,   p=0     Wald test            = 240.3  on 14 df,   p=0     Score (logrank) test = 288.5  on 14 df,   p=0      Now if I understand my results correctly, this is telling me that the variable `FAInadequate` reduces the hazards for the event which in this case is "requirement for an operation", right? Now this is a bit bizarre, because if I perform simple comparative statistics, I find that FAInadequate patients have longer healing times and more operations as shown by these results:               t.test(log(Timetoheal+0.5)~FA, var.equal=T)              Two Sample t-test          data:  log(Timetoheal + 0.5) by FA     t = -7.9285, df = 4916, p-value = 2.724e-15     alternative hypothesis: true difference in means is not equal to 0     95 percent confidence interval:      -0.2773227 -0.1673666     sample estimates:       mean in group Adequate mean in group Inadequate                      2.338534                 2.560879       and this               require(gmodels)     CrossTable(Operation, FA, format="SPSS", chisq=T)             Cell Contents     |-------------------------|     |                   Count |     | Chi-square contribution |     |             Row Percent |     |          Column Percent |     |           Total Percent |     |-------------------------|          Total Observations in Table:  4877                        | FA          Operation|   Adequate | Inadequate |  Row Total |      -------------|------------|------------|------------|               No |      2583  |      1719  |      4302  |                   |     2.764  |     3.835  |            |                   |    60.042% |    39.958% |    88.210% |                   |    91.143% |    84.141% |            |                   |    52.963% |    35.247% |            |      -------------|------------|------------|------------|              Yes |       251  |       324  |       575  |                   |    20.682  |    28.690  |            |                   |    43.652% |    56.348% |    11.790% |                   |     8.857% |    15.859% |            |                   |     5.147% |     6.643% |            |      -------------|------------|------------|------------|     Column Total |      2834  |      2043  |      4877  |                   |    58.109% |    41.891% |            |      -------------|------------|------------|------------|               Statistics for All Table Factors               Pearson's Chi-squared test      ------------------------------------------------------------     Chi^2 =  55.971     d.f. =  1     p =  7.354796e-14           Pearson's Chi-squared test with Yates' continuity correction      ------------------------------------------------------------     Chi^2 =  55.29973     d.f. =  1     p =  1.03483e-13                       Minimum expected frequency: 240.8704       The fact of the matter is, that in plain "medical" terms, the longer it takes a wound to heal, the worse the outcome for the patient, so I would've originally expected the inadequate treatment to increase the hazards, especially in light of the simple t-test cross-tabulation results. Or could it be, that in light of this, running a proportional hazards function is incorrect and I should be looking at cumulative hazards? If that is the case how can I run a cumulative hazards function in R? Further, I'm a bit concerned, that the likelihood ratio, logrank and Wald test all resulted in p-values = 0. Finally, how would you ideally chose to cross-validate this model and adjust for the slightly "too low" event per variable rate? I'm sorry for all the questions, but this has been thoroughly been driving me nuts over the last couple of days, as I can't really wrap my head around the output for some reason. I appreciate any help greatly. Thanks.


title: IV First and Second Stage Graphs and Interpretations

text: First of all, many thanks to those who have kindly answered my previous questions, I am learning a lot! Suppose I have a model y = x1 + x2 + x3 x1 is endogenous and has instrument Z1 that is relevant and exogenous conditional upon X2 and X3. So First Stage: X1 = Z1 + X2 + X3 Reduced Form: Y = Z1 + X2 + X3 Structural Equation Y = (X1 = Z1) X1 X2 My questions are: Is it convention to plot the coefficient on Z1 against X1 to illustrate relevance? Also to plot the coefficient on Z1 in the reduced form on y. Additionally, I am confused about the relationship in the structural equation. I understand the relationship of the coefficient on Z1 in the first stage on X1, but in the second stage we substitute x1 not only with the coefficient on Z1 in the first stage but the fitted value of X1, which entails the coefficients on X2 and X3 in addition to Z1 - I don't see why and this clouds my understanding of who the instrument impacts and to what magnitude.


title: Lag order for Granger Causality Test

text: Suppose I'm considering several independent variables for possible inclusion in an ARIMAX model I'm developing. Before fitting different variables, I'd like to screen out variables that exhibit reverse causality by using a Granger test (I'm using the `granger.test` function from the `MSBVAR` package in R, although, I believe other implimentations work similarly). How do I determine how many lags should be tested? The R function is: granger.test(y, p), where y is a df or matrix, and p is the lags. > The null hypothesis is that the past p values of X do not help in predicting > the value of Y. Is there any reason not to select a very high lag here (other than the loss of observations)? Note that I have already differenced _every_ time series in my df, based on the order of integration of my dependent time series. (e.g., differencing my _dependent_ time series once made it stationary. Therefore, I also differenced all "independept" time series once.)


title: What test should I use for hypothesis with many variables?

text: I have a survey with a list of questions in which respondents can pick their perceived personality traits. Then I have a question that asks respondents about their favorite in-game character's perceived personality and a question that asks respondents about their least favorite in-game character's perceived personality trait. My hypothesis is: People will prefer characters that they perceive to be similar to themselves, and dislike ones that they perceive to be different. What is the best way that I can conduct this statistic test? Example: One respondent may answer like this:   * Extrovert (1-7): 2   * Dependable (1-7): 5   * Reserved (1-7): 4 Favorite Character:   * Extrovert: no   * Dependable: yes   * Reserved: no Least Favorite Character:   * Extrovert: no   * Dependable: no   * reserved: no I know how to test simple hypotheses about two simple numeric variables, but this looks a lot more complex.


title: What are the units of the second derivative?

text: Generally, what are the units of the second derivative? In this (very specific) example I have plotted height against gap probability, what should I put in place of the "?". ![enter image description here](http://i.stack.imgur.com/fITpL.png)


title: Sparse PCA adjusted variance explained

text: In the original Sparse PCA paper Sparse Principal Component Analysis ZOU, HASTIE, TIBSHIRANI they describe a way to compute the adjusted variance explained by computing QR decomposition of the Z matrix (size n times k) where n is the number of observations and k is the low dimensional space. Then the adjusted variance is just trace(R^2). I tried doing that for probabilistic sparse PCA that I implemented but it turns that the adj. variance explained is a very high number ~700. We know that it should be between 0 and 1, so probably I am missing some normalization constant. Does anyone have any idea on what could be going wrong or how to normalize. My original data dimensions are 100*50 i.e. n=100, p=50 and k=10 Thanks!


title: Is it valid to compare p-values from test statistics with different DF?

text: I need to compare results from data with different residual DF as my $x$ variable has different levels. The following is just an example (in R, for demonstration purpose, but this is not a R question):               # first case      set.seed (123)     data1 <- data.frame (y = rnorm (100, 5, 2),       x = sample (c("A", "B"), 100, replace = T))     anova(lm(y~ x, data = data1))     Analysis of Variance Table          Response: y               Df Sum Sq Mean Sq F value Pr(>F)     x          1   2.07  2.0669  0.6177 0.4338     Residuals 98 327.89  3.3459                         # second case:       set.seed (123)     data2 <- data.frame (y = rnorm (100, 5, 2),       x = sample (c("A", "B", "C", "D", "E"), 100, replace = T))     anova(lm(y~ x, data = data2))     Analysis of Variance Table          Response: y               Df Sum Sq Mean Sq F value Pr(>F)     x          4   4.89  1.2224  0.3572 0.8384     Residuals 95 325.07  3.4218             Here I have two different DF for the residuals (95 vs. 96) and $x$ (1 vs. 4): Is it valid to compare p-values as such? I know that the F-test considers $x$ and residual while calculating p-value. Is there any extra-caution needed?


title: How to draw mean, median and mode lines in R that end at density?

text: I have drawn a skewed distribution, using these commands:               x <- seq(-2.5, 10, length=1000000)     hx5 <- rnorm(x,0,1) + rexp(x,1/5) # tau=5 (rate = 1/tau)     plot(density(hx5), xlim=c(-2.5,10), type="l", col="green",          xlab="x", main="ExGaussian  curve",lwd=2)      Now I want to draw three lines, for the mean, the mode and the median of the distribution. If I simply write, for example:               abline(v=median(hx5))      the line go out of the curve, but I want to end the line with the density point of the parameter. So, my problem is: **how can I find the values of the density at the mean, mode and median of my observations in order to set the correct coordinates for drawing?**


title: Combine MultiClassClassifier and MetaCost in WEKA and make cost matrices for each binary classifier

text: Ok, I have a classification problem in which i need to classify instances into one of 5 classes. The class distribution of the classes is imbalanced however. These are the frequencies of instances for each class:               class 1: 129      class 2: 2      class 3: 187      class 4: 18      class 5: 285      So this is a multi-class classification problem. With imbalanced class distributions like this, I read papers that recommend to use cost sensitive classifiers. Most literature that I read however has pretty much set out how to do this with binary classification problems but not for multi-class problems. One important aspect, making a cost matrix, is relatively easy when you have a binary class problem. For example if you have 10 A and 90 B instances, you can make a cost matrix that penalizes classifying A as B with 9, while penalizing classifying B as A only with 1:               0 1     9 0      So I would like to be able to make a binary classifier problem of my multiclass problem. For this I use the following things in weka: First of all I use the metaclassifier: MultiClassClassifier. In this MultiClassClassifier I put a MetaCost classifier, which on its turn uses for example a J48 classifier. My problem is, that I need to define a cost matrix in that MetaCost classifier. I want that cost matrix to be build differently for each round that the MultiClassClassifier makes. Is there any way to do this?


title: Sample size determination for polling with multiple choice type questions

text: > **Possible Duplicate:**   >  Question regarding sampling, estimation and accuray Let's say there are N sets of diamonds, containing 5 in each set. There can be anywhere between 1 to 5 defective diamonds in each set. The goal is to get an estimate of total defective diamonds in all sets combined by randomly sampling M sets. What formula/theory should I use to calculate the sample size M, the estimate and the accuracy of the estimation?


title: How should I interpret the interaction term in a Cox proportional hazards model?

text: How should I interpret the following interaction term of 2 continuous predictors in the output of a Cox proportional hazards model? The Hazard ratio for the interaction of X and Y is >1, which means its log (the original coefficient) is 0-1 (~0.16). The individual items have a HR less than one, and coefficients of X= -0.18 and Y=-0.11.                   |   Variable                   | HR (s.e.)     | p value       -----------------------------------------------------------           1 A (5 points)                 |0.756 (0.088)  |    0.001 |               2 B (5 points)                 |1.379 (0.11)   |    0.001 |               3 X  (10 points)               |0.837 (0.033)  |    0.0011|               4 Y  (1 point)                 |0.895 (0.03)   |     0.001|               5 X (10 points)x Y (1 point)   |1.016 (0.006)  |    0.011 |       The effect of a 10 point increase in X, with Y=0 is to decrease the "death" rate by 16%. The effect of a 1 point increase in Y, with X=0, is to decrease the death rate by 10.5%. What is the effect of a one-point increase in Y on the effect of a 10-point increase in X on the death rate? X has a range of 0 to 90. Y has a range of 0 to 10. With a one-point increase in Y, does the effect of a 10-point increase in X increase from 16% to (16% + 1.6%) = 17.6%, or does it decrease by 1.6% to 14.4%? Thought I had it straight but now very stuck here.


title: Linear model overfitting  due to too many covariates

text: My study design involves a control and 2 test groups plus some covariates. Each group consists of around 20 observations. In total I look at around 1,000 variables. I created a linear model using the `lm()` function in R including 2 covariates. After that I thought I would include another covariate because doing a PCA plot earlier showed a slight effect on that covariate. However, after adding this covariate to the model 50% of the significant hits are now different. I was actually assuming that it would pretty much identical as the effect was hardly seen in the PCA. Could it be that I have overfitted the model? Or is the effect simple just not shown in the PCA plot but is there? I just compared the two models using `anova(lm1, lm2)` and the p-value is significant which I think means that the third covariate adds significant information to the model?               lm1 <- lm(var ~ factor_of_interest + cov1 + cov2)     lm2 <- lm(var ~ factor_of_interest + cov1 + cov2 + cov3)     anova(lm1, lm2)


title: What kind of cut-offs are used for publication bias tests in meta-analysis?

text: I have performed meta-analysis and included four different tests for publication bias in the analysis - classic fail-safe N, Orwin's fail-safe N, Duval and Tweedie's trim and fill as well as funnel plots. The tests did not seem to indicate publication bias in my material - according to classic fail- safe N 23 papers would be needed to bring the result to p > 0.05, according to Orwin's fail-safe N 8 papers would be needed to reduce the OR into a trivial one, trim and fill test did not show adjustments and funnel plot was symmetrical. However, one of my reviewers wants me to include cut-offs for publication bias tests in the Methods. I have tried to look for any cut-offs that would be commonly used with these tests but in many cases even the interpretation of results is not very exact. How should I proceed with this problem, what should I answer to the reviewer? Thanks!


title: How to interpret OOB and confusion matrix for random forest?

text: I got a an R script from someone to run a random forest model. I modified and run it with some employee data. We are trying to predict voluntary separations. Here is some additional info: this is a classification model were 0 = employee stayed, 1= employee terminated, we are currently only looking at a dozen predictor variables, the data is "unbalanced" in that the term'd records make up about 7% of the total record set. I run the model with various mtry and ntree selections but settled on the below. The OOB is 6.8% which I think is good but the confusion matrix seems to tell a different story for predicting terms since the error rate is quite high at 92.79% **Am I right in assuming that I can't rely on and use this model because the high error rate for predicting terms? or is there something also I can do to use RF and get a smaller error rate for predicting terms?**                FOREST_model <- randomForest(theFormula, data=trainset, mtry=3, ntree=500, importance=TRUE, do.trace=100)     ntree      OOB      1      2       100:   6.97%  0.47% 92.79%       200:   6.87%  0.36% 92.79%       300:   6.82%  0.33% 92.55%       400:   6.80%  0.29% 92.79%       500:   6.80%  0.29% 92.79%     > print(FOREST_model)          Call:      randomForest(formula = theFormula, data = trainset, mtry = 3,      ntree = 500, importance = TRUE, do.trace = 100)                     Type of random forest: classification                          Number of trees: 500     No. of variables tried at each split: 3                  OOB estimate of  error rate: 6.8%     Confusion matrix:          0  1 class.error     0 5476 16 0.002913328     1  386 30 0.927884615     > nrow(trainset)     [1] 5908


title: Calculating E[log(g(x))]

text: Let $g(x)=a+bx$ and let $x$ have known distribution $f(x)$, say. Is it possible to calculate: $$\mathbb{E}[\log(g(x))] = \int f(x) \log{g(x)}\,\mathrm{d}x= \int \log(a + b x) f(x) \,\mathrm{d}x?$$ Many thanks for all of the help in advance.


title: Predictive Modeling - Should we care about mixed modeling?

text: For predictive modeling, do we need to concern ourselves with statistical concepts such as random effects and non independence of observations (repeated measures)? For example.... I have data from 5 direct mail campaigns (occured over the course of a year) with various attributes and a flag for purchase. Ideally, I would use all this data combined to build a model for purchase given customer attributes at the time of the campaign. The reason is that the event of purchase is rare and I would like to use as much information as possible. There is a chance that a given customer could be on anywhere from 1 to 5 of the campaigns - meaning there is not independence between the records. Does this matter when using: 1) A machine learning approach (e.g. tree, MLP, SVM) 2) A statistical approach(logistic regression)? **ADD:** My thought about predictive modeling is if the model works, use it. So that I have never really considered the importance of assumptions. Thinking about the case I describe above got me wondering. Take machine learning algorithms such as a `MLP and SVM`. These are used successfully to model a binary event such as my example above but also time series data that are clearly correlated. However, many use loss functions that are likelihoods and derived assuming the errors are iid. For example, gradient boosted trees in R `gbm` uses deviance loss functions that are derived from the binomial (Page 10).


title: Power analysis for ordinal logistic regression

text: I am looking for a program (in R or SAS or standalone, if free or low cost) that will do power analysis for ordinal logistic regression.


title: $E(\min_{1 \leq i \leq m} |x - y_i|)$

text: I have faced a statistical problem in a computational biology question. I will explain it in statistical language. Let's say we have the results of rolling one red die, $x$, and $m$ black dice $y_1,\ldots,y_m$. These two types of dice are physical identical and both have $N$ faces. The question is finding the average expected minimum absolute differences of values between the red and black ones. $E(\min_{1 \leq i \leq m}|m-y_i|)$ where $x$ and $y$s are discrete and uniformly distributed. My main difficultly is mainly on $E(\min)$ of $M$ random variables. Is there any R function or package regarding such a computation? **Update:** Here is my simulation in R:               nx = 1 # number of red dice , x     ny = 4 # number of black dices, y     N = 7 # number of faces in dices          itr = 14000 # iteration      dis.vec = rep(0,itr)     for (cnt in 1:itr)     {     pool = rep(0, N)     xz = sample(c(1:N),nx)     yz = sample(c(1:N),ny)     dis.vec[cnt] = min(abs(yz - xz)) #calculate the minimum     }     print(mean(dis.vec)) #average!


title: Need help to understand the meaning of this formula

text: I need help to understand the meaning of the following formula: $$\sum_{i=1}^n\sum_{j\neq i,j=1}^n(x_i-\mu_i)(x_j-\mu_j)c_{ij}$$ where C is positive semidefinite matrix, $\mu_i$ is the mean of $x_i$. I was thinking, since $c_{ij}$ is the covariance between $x_i$ and $x_j$, this formula kind of sums up the correlation between all the pairs of $x_i$ and $x_j$. Does that make sense? Do you have any other interpretation of the formula? Many thanks in advance, **Edited v1** Matrix $C$ is simply a positive semidefinite matrix, not necessarily covariance matrix. $\mu $ was given very general as well, not mentioned to be mean. That was my own interpretation of it. If you know of another formula using $\mu$ with different meaning, please do tell. Many thanks for Emre and Nesp for the response. So far, it really makes sense. **Edited v2** Thank you all very much for your responses. The Mahalanobis distance is a very good explanation of the formula given. However, I wonder if we take off the diagonal (all the i=j elements) from the formula, does it still make any sense? Without the condition $i\neq j$, the summation can even be negative, right? In that case, the interpretation of distance seems to be not appropriate anymore. (I made the modification to the formula above)


title: Frequentist ability to use observations

text: I'm not an expert in statistics, but I gather there is disagreement whether a "frequentist" or "Bayesian" interpretation of probability is the "right" one. From Wagenmakers et. al p. 183: > Consider a uniform distribution with mean μ and width 1. Draw two values > randomly from this distribution, label the smallest one $s$ and the largest > one $l$, and check whether the mean μ lies in between $s$ and $l$. If this > procedure is repeated very many times, the mean μ will lie in between s and > l in half of the cases. Thus, $(s, l)$ gives a 50% frequentist confidence > interval for μ. But suppose that for a particular draw, $s = 9.8$ and $l = > 10.7$. The difference between these values is 0.9, and this covers 9/10th of > the range of the distribution. Hence, for these particular values of $s$ and > $l$ we can be 100% confident that $s < μ < l$, even though the frequentist > confidence interval would have you believe you should only be 50% confident. Are there really people who believe that there is only 50% confidence in this case or is it a straw man? I guess more generally, the book seems to be saying that frequentists can't express a conditional claim like "Given s = 9.8 and l = 10.7, $s<\mu<l$ with probability 1". Is it true that conditioning implies Bayesian reasoning?


title: Use of propensity score in a case-control study

text: The theory of propensity score (PS) suggests that it should only be used for the cohort study because PS matches the "treated/exposed" to the "un- treated/un-exposed" groups. However, cases and controls are the outcomes (not the exposures) in a case-control study. What are the pitfalls if PS is used in a case-control study?


title: Clustering the elliptic Fourier series coefficients

text: What do you think about the best way to perform a search for two discrete group membership in data that consist of four Fourier series coefficients like this:               $an      [1] -1.099759e+02 -3.277411e+01 -2.069982e+01  1.294945e+00 -4.557681e+00      [6]  3.659611e+00 -2.266350e+00  2.290051e+00 -1.165922e+00  4.290447e-01     [11] -5.157878e-01 -4.373163e-02 -9.627400e-02 -2.160347e-01 -8.664855e-02     [16] -7.125669e-02 -1.661091e-01  5.034174e-02 -5.773743e-02  9.070195e-02     [21] -1.103436e-01  4.968477e-02 -3.843040e-02 -3.539704e-02 -2.001452e-03     [26] -6.103696e-02  5.080494e-03 -8.438730e-02  9.728283e-02 -1.201637e-01     [31]  5.214383e-03 -3.097534e-02          $bn      [1] -9.950280480 21.060718135  1.749972609  3.125969744 -1.221718994      [6] -0.366008050  0.134939868 -1.433288279  1.286034494 -0.793100554     [11]  1.024338370 -0.645426109  0.752002518 -0.354023866  0.371995040     [16] -0.079819647  0.093705335 -0.085753179  0.062387444 -0.007351952     [21]  0.112815714 -0.071472761  0.180375050 -0.155871666  0.165280059     [26] -0.084238927  0.157971637 -0.115017123  0.058533444  0.021303459     [31] -0.009807206 -0.037191448          $cn      [1] 186.42631186 -48.82963771 -13.05530543 -17.63800716  -8.78793762      [6]  -5.13169720  -4.05998782  -1.53832276  -1.92963513  -1.04525653     [11]  -0.80438591  -0.98701575  -0.24698734  -1.00053803  -0.21660353     [16]  -0.85765069  -0.32822392  -0.54232592  -0.39971080  -0.22425141     [21]  -0.46861898  -0.05909403  -0.45459051  -0.03218630  -0.39596045     [26]  -0.06020157  -0.25675128  -0.08580023  -0.16002496  -0.14539536     [31]  -0.13190927  -0.12939389          $dn      [1] -60.193804743  17.089146986  -5.359459446   9.834089998  -1.841459861      [6]   2.894124403  -1.229890034  -0.060334681  -0.383461709  -0.414984724     [11]   0.161235406  -0.209664090   0.133169983   0.114066031  -0.113437218     [16]   0.326989075  -0.353357830   0.366442458  -0.379583013   0.327422310     [21]  -0.253858236   0.177601217  -0.120553856   0.010128884   0.009787827     [26]  -0.048527754   0.075992199  -0.090581996   0.043151283  -0.050768870     [31]   0.022838609   0.004311564      when I have 26 of such sets (and choosing the first 8 of each coefficient excluding the first one). These coefficients describe shape in mountain goat horns, and I need the best way to separate them according to sex which is accounting for horn curvature and other shape features. Now, PCA or LDA separate the extreme shapes well but what about the intermediate ones? I tried numerous partitional clustering approaches but found no consistent grouping pattern with them. Also the majority of parametric statistic methods do not work because of high degree of correlation and non-normality of such data. Any advice?


title: Estimating synergy of treatments in time series data

text: I am studying the growth (in volume) of tumors in mice when they are given 1) saline solution (control), 2) treatment A, 3) treatment B, and 4) treatment A and B. The treatments are designed to inhibit the growth of the tumors. I am hoping to estimate whether there is a "synergistic" relationship between A and B. By synergy, I mean that the combination treatment of A and B leads to what we would call, above some uncertainty threshold, a greater than additive effect. My data will be in time series format, as the measurements will be made longitudinally over many days, with replicates in each group. Here's a quick literature review. Many studies cite and use Chou's method, but that will not work for us because it requires a separate dose-response curve for each drug, which we will not have. I also found one study that uses a regression-based model (Tan et al), which employs the model $ y_i = \alpha + y_{i0}1_{12}\beta_0 + x^1_i\beta_1 + x_i^2\beta_2 + + x_i^3\beta_3 + \epsilon_i $ where $y_i$ is the log tumor size, the $x_i$'s refer to the quantity of the treatments used, and the $\beta$'s refer to the parameters. They then do expectation-maximization to get MLE's of the parameters. In this setting, $\beta_3$ is the crucial interaction parameter which we want to estimate as either significantly less than zero or not. However, the authors do not mention access to their implementation of the method and it seems difficult to re-implement. Moreover, our treatments will be binary rather than two-leveled as in their data set. I also found a study that uses an ANOVA (Soundararajan et al) and presumably post-hoc tests, but I didn't exactly understand their methods. I am not looking to reinvent the wheel and I would prefer to use existing R packages and/or software. If the latter, I'd prefer open-source, to boost reproducibility. Any suggestions?


title: Advantage of central moment over moment?

text: Is there any advantage of using "central moments" over "moments" when approximating a distribution to a known distribution using moment matching? I have noticed that in lot of papers.


title: Repeatability and measurement error from and between observers

text: I have 3 observers that each take 2 measurements (length and weight) on 100 individuals; these procedures are repeated once (i.e., the same measurements are taken on the same 100 individuals by the same 3 observers), so that the data set is duplicated (i.e., early reading and late reading).   1. What is the best way to figure out how each individual observer's measurements vary between the late and early trial measurements?   2. How can I best compare how close or different the measurements of length (or weight) differ among the 3 observers?


title: Resources about probability proportional to size (PPS) sampling method

text: Following the question here, Probability Proportional to Size (PPS) Sampling Method seems appropriate in constructing a sample when the primary units are of different sizes in the target population. As I understand, this is a self-weighting technique, so does not require any special consideration in the data analysis stage. It is also a probability sampling technique so meets the assumption of correlation and regression models. I am looking for some good resources that can help me better understand this technique. Perhaps a flowchart that illustrates this method would be great to show to my 'non-statistical' target audience.


title: Reference for random forests

text: I would like to understand how do the Boruta package work. Could you suggest some references for the theoretical aspect of so-called random forests? Thanks. Below are two illustrative examples of why am I attracted by the Boruta algorithm. _First example_ :               > set.seed(666)     > # simulates data     > # y does not depend on x4     > # y depends on x1, x2, x3 only through x3     > x1 <- rnorm(50); x2 <- rnorm(50) ; x3 <- (x1+x2)^2; x4 <- rnorm(50)     > y <- x3+rnorm(50,0.1)     >     > # lm() only indicates x3 is "important"     > summary(lm(y~x1+x2+x3+x4))          Call:     lm(formula = y ~ x1 + x2 + x3 + x4)          Residuals:         Min      1Q  Median      3Q     Max      -3.0000 -0.7137  0.0352  0.7082  1.7918           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)   0.1875     0.1874   1.000    0.323         x1           -0.1541     0.1469  -1.049    0.300         x2            0.1153     0.1949   0.591    0.557         x3            0.9097     0.0501  18.160   <2e-16 ***     x4            0.1263     0.1518   0.832    0.410         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 1.123 on 45 degrees of freedom     Multiple R-squared: 0.9013,     Adjusted R-squared: 0.8925      F-statistic: 102.7 on 4 and 45 DF,  p-value: < 2.2e-16           > # Boruta indicates x1, x2, x3 are important     > Boruta(y~x1+x2+x3+x4, maxRuns=500)     Boruta performed 174 randomForest runs in 10.409 secs.             3 attributes confirmed important: x1 x2 x3             1 attributes confirmed unimportant: x4      _Second example_ :               > set.seed(421)     > # simulates data     > # y does not depend on u1     > # y does not depend on u2     > # but y depends on u1+u2     > nsims <- 100     > u1 <- runif(nsims)     > u2 <- runif(nsims)     > x <- (u1+u2)-floor(u1+u2)     > y <- rnorm(nsims, x,.05)     >      > # lm() does not detect some dependence     > summary(fit <- lm(y~u1))          Call:     lm(formula = y ~ u1)          Residuals:         Min      1Q  Median      3Q     Max      -0.5411 -0.2234  0.0005  0.2197  0.4901           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)  0.55721    0.05902   9.442 1.97e-15 ***     u1          -0.07858    0.09714  -0.809    0.421         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.2668 on 98 degrees of freedom     Multiple R-squared: 0.006633,   Adjusted R-squared: -0.003503      F-statistic: 0.6544 on 1 and 98 DF,  p-value: 0.4205           > summary(fit <- lm(y~u2))          Call:     lm(formula = y ~ u2)          Residuals:          Min       1Q   Median       3Q      Max      -0.53996 -0.21855  0.01298  0.22406  0.49940           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)  0.47022    0.05124   9.178 7.38e-15 ***     u2           0.09435    0.09298   1.015    0.313         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.2663 on 98 degrees of freedom     Multiple R-squared: 0.0104,     Adjusted R-squared: 0.000299      F-statistic:  1.03 on 1 and 98 DF,  p-value: 0.3127           > summary(fit <- lm(y~u1+u2))          Call:     lm(formula = y ~ u1 + u2)          Residuals:          Min       1Q   Median       3Q      Max      -0.53400 -0.22071  0.00699  0.21612  0.54375           Coefficients:                 Estimate Std. Error t value Pr(>|t|)         (Intercept)  0.51609    0.06840   7.545 2.45e-11 ***     u1          -0.09978    0.09859  -1.012    0.314         u2           0.11176    0.09455   1.182    0.240         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.2663 on 97 degrees of freedom     Multiple R-squared: 0.02074,    Adjusted R-squared: 0.0005478      F-statistic: 1.027 on 2 and 97 DF,  p-value: 0.3619           >      > # Boruta() does     > Boruta(y~u1)     Boruta performed 44 randomForest runs in 6.328125 secs.             No attributes has been deemed important             1 attributes confirmed unimportant: u1     > Boruta(y~u2)     Boruta performed 20 randomForest runs in 2.8125 secs.             No attributes has been deemed important             1 attributes confirmed unimportant: u2     > Boruta(y~u1+u2)     Boruta performed 48 randomForest runs in 6.796875 secs.             2 attributes confirmed important: u1 u2             No attributes has been deemed unimportant


title: Finding the effect of time from a multilevel model?

text: I would like to determine the effect of AGE or equivalently "time" on a continuous outcome variable while controlling for gender. I have repeated measures data across 6 time points. In order to take into account the correlation of the repeated measures, I am thinking of using a random slops/random intercepts model multilevel model where "measurements at different ages" are clustered within individuals. I am a bit confused about interpretation of model coefficients. What coefficients should I look at to the extract the effect of a one unit increase in time/age on my outcome? N26


title: Dropping unused levels in facets with ggplot2

text: Is it possible to drop levels that are not used in ggplot2s facets? This is my code:               tab = as.data.frame(cbind(groups = mtcars$cyl, names = row.names(mtcars), val = mtcars$mpg, N = mtcars$disp))     tab$N = as.numeric(tab$N)          ggplot(tab, aes(names,val)) +      geom_point() + coord_flip() +      theme_bw() +      facet_grid(groups ~ ., drop=TRUE)#, scales="free", as.table=F, space="free")      I tried the `drop=T` switch but it doesnt help. What am I doing wrong? Thanks for any help.


title: Regression vs classification within decision trees in OpenCV

text: I have been looking into the random trees and decision trees in OpenCv and something that still confuses me is the differences between regression and classification. Can anyone outline what the differences are, when would one use one over the other? The reason for me asking is that i am using classification for labeling a multilabel class of a dataset, and then i ran across the code for the article "Object Classification using Hough Forest" and its source code. In their code, the methods are called regression and its also doing a multiclass labeling, so now i am not sure if I know the difference? A little more about my case: I am extracting N descriptors from images, each descriptor are different computed features from an image patch, forming a 1x8096 descriptor. Each descriptor are labelled with a label 0:4. I then use OpenCV RandomTrees to train and classify. Looking at the code, its using the methods for ordered classification. It also have categorical classification/regression and categorical regression. So from what i read, would the right term not be that im using a regression, my values can be ordered continuous: > In brief: > >   1. Classification trees have dependent variables that are categorical and > unordered >   2. Regression trees have dependent variables that are continuous values or > ordered whole values. > http://www.differencebetween.com/difference-between-classification-and-vs- regression/