title: Validity of Log-rank p-value when there are zero events in one group?

text: I am looking to do a Kaplan Meier analysis to look at survival in two groups. I have 21 people in group A and 17 people in group B. 5 people in group B died at various times, 0 people in group A died, ie all my data is censored for Group A. I am getting a log-rank p-value of 0.028 when I do the KM analysis. My question is can I still use the KM for survival analysis given the absence of events in one group and if so, is the p-value valid. Thanks in advance


title: Discretizing a continuous variable in SAS using a decision tree

text: I have a variable X15 that varies from .001 to 10,000 or so and I'd like to construct classes that are unbalanced to another variable, Y, that is either 1 or 0, so that each class will have more of (Y=0) or of (Y=1). All I want is the output from some PROC in 9.3 that will tell me what the bins should be. I've tried PROC SPLIT, PROC DMSPLIT, and PROC ARBORETUM, but PROC PRINTing all the output data sets doesn't give me a clue. How do I get the splitting values for, say, 5 splits into 6 bins of my X15 data using a decision tree? An arbitrary number of splitting values is preferred, I just can't figure out how to extract that data from a PROC statement.


title: Show $Y$ converges to $a$

text: Given: $f_{Y_{(1)}}(y) = nbe^{-nb(y-a)}$, where $b> 0$ and $y \geq a$. Show that as $n \rightarrow\infty$, $Y_{(1)}$ converges to $a$ in probability. I have calculated $E[Y_{(1)}] = \frac{1}{nb} + a$ Which theorem should I apply to show convergence. I was trying to use Chebyshev's inquality. [EDIT] If I want to find to what $Y_{(1)}$ converges in distribution, is this the right way to do it: $F_{Y_{(1)}} = 1-e^{-nb(y-a)}$ As $n \rightarrow \infty $ $F_{Y_{(1)}} = 1, y < a$ $F_{Y_{(1)}} = 0, y \geq a$ $P(|Y_{(1)}|< y) = P(|Y_{(1)}|< \epsilon)$ [replacing $y$ with $\epsilon$] = $1-e^{-nb(\epsilon-a)}$ As $n \rightarrow \infty $ $P(|Y_{(1)}|< y) \rightarrow 1$; So $ Y_{(1)} \rightarrow Y$ in distribution. So the limiting distribution is degenerate. Please let me know if this approach is correct.


title: Naive SE vs Time Series SE: which statistics should I report after Bayesian estimation?

text: I am new to Bayesian estimation. When I do some estimations with JAGS, I find there are statistics called Naive SE and Time Series SE. What exactly do they mean? Is it necessary that I report one or both of them as part of the estimation result?


title: Convergence in distribution

text: For a statistic $T_n = \frac{1}{n} \sum_{i=1}^nY_i - \frac{1}{a}$. Prove directly (without CLT) that scaled and appropriately shifted version of $T_n$ converges in distribution to $N(0,1)$. [EDIT] $f(y|a,b)=ae^{-a(y-b)}$ for $ y\geq b$ How should I approach the problem? [EDIT] I thought that if I could find the expected value and variance of $T_n$ and then represent it in $N(\mu,\sigma^2)$ for $f(y)$, $E[Y] = b+1/a$, $Var [Y] = 1/a^2$ For $Tn$, $Var[T_n] = Var[\frac{1}{n} \sum_{i=1}^nY_i - \frac{1}{a}]$ = $a^2/n$ $E[T_n] = E [\frac{1}{n} \sum_{i=1}^nY_i - \frac{1}{a}] = b$ [EDIT] $M_X(t)=\frac{ae^{bt}}{(a-t)}$ How should I go from here?


title: How can I test the difference of two Weibull distributions?

text: I have two Weibull distribution sets from two wind datasets in order to check whether they are same. I thought a 2 sample t-test would be applicable but I couldn't find any ways to do that on the Internet. Does anyone know what type of test is applicable to my purpose? and what R function can you recommend? Plus, if it turned out that there is a difference between the two datasets, can I just fit a linear line between the datasets?


title: Sane stepwise regression?

text: Suppose I want to build a binary classifier. I have several thousand features and only a few 10s of samples. From domain knowledge, I have a good reason to believe that the class label can be accurately predicted using only a few features, but I have no idea _which_ ones. I also want the final decision rule to be easy to interpret/explain, further necessitating a small number of features. Certain subsets of my features are highly correlated, so selecting the most predictive few independently wouldn't work. I also want to be able to meaningfully do hypothesis testing on my features. Is the following stepwise regression procedure reasonable under these conditions:   1. Given the features already in the model (or just the intercept on the first iteration), select the feature that produces the largest log likelihood ratio when added to the model. Use the likelihood ratio chi-square test to calculate a nominal P-value for each hypothesis test performed in this selection. The null here is that adding the extra variable to the model provides no additional predictive ability. The alternative is that it does increase predictive abilityl   2. Treat the hypotheses tested in Step 1 of each iteration as a family and calculate the false discovery rate for the smallest P-value (for the feature selected) using something like Benjamini-Hochberg.   3. Goto 1 unless some stopping criteria are met.   4. Report the false discovery rates for the individual features, but **_not_** the P-value for the model as a whole (since this will be massively inflated). Each of these multiple testing corrected P-values represents the statistical significance of that feature _given_ all of the features previously added to the model. Does doing something like this under these circumstances successfully avoid all of the typical criticisms of stepwise regression? Are the false discovery rates calculated in this way reasonable?


title: Need to convert duplicate column elements to a unique element in R

text: I need to convert the table A to table B. How can I do that using R? **TABLE A**               Y 10     Y 12     Y 18      X 22      X 12     Z 11     Z 15      **TABLE B**               X 22 12     Y 10 12 18     Z 11 15


title: How can one use the predict function on a lm object where the IVs have been dynamically scaled?

text: **Note: I've updated the example case code, there were some errors in the previous version** Cross posted to R-help, because I half suspect this is 'unexpected behaviour'. I want to predict values from an existing lm (linear model, e.g. lm.obj) result in R using a new set of predictor variables (e.g. newdata). Specifically, I am interested in the predicted y value at the mean, 1 SD above of the mean, and 1 SD below the mean for each predictor. However, it seems that because my linear models was made by calling scale() on the target predictor that predict exits with an error, "Error in scale(xxA, center = 9.7846094491829, scale = 0.959413568556403) : object 'xxA' not found". By debugging predict, I can see that the error occurs in a call to model.frame. By debugging model frame I can see the error occurs with this command: variables <\- eval(predvars, data, env); it seems likely that the error is because predvars looks like this: list(scale(xxA, center = 10.2058714830537, scale = 0.984627257169526), scale(xxB, center = 20.4491690881149, scale = 1.13765718273923)) An example case:                  dat <- data.frame(xxA = rnorm(20,10), xxB = rnorm(20,20))        dat$out <- with(dat,xxA+xxB+xxA*xxB+rnorm(20,20))        lm.res.scale <- lm(out ~ scale(xxA)*scale(xxB),data=dat)        my.data <- lm.res.scale$model #load the data from the lm object        newdata <- expand.grid(X1=c(-1,0,1),X2=c(-1,0,1))        names(newdata) <- c("scale(xxA)","scale(xxB)")        newdata$Y <- predict(lm.res.scale,newdata)      Is there something I could do before passing newdata or lm.obj to predict() that would prevent the error? I tried: From the help file it looks like I might be able to do something with the terms, argument but I haven't quite figured out what I would need to do. Alternatively, is there a fix for model.frame that would prevent the error? Should predict() behave this way? **Additional Details:** However, I really want a solution that, in one step will provide values like:               coef(lm.res.scale)[1]+     coef(lm.res.scale)[2]*newdata[,1]+     coef(lm.res.scale)[3]*newdata[,2]+     coef(lm.res.scale)[4]*newdata[,1]*newdata[,2]      I think that should be exactly what predict() should do. That is, I think my example code should be equivalent to:                  dat <- data.frame(xxA = rnorm(20,10), xxB = rnorm(20,20))        dat$out <- with(dat,xxA+xxB+xxA*xxB+rnorm(20,20))        #rescaling outside of lm        X1 <- with(dat,as.vector(scale(xxA)))        X2 <- with(dat,as.vector(scale(xxB)))        y  <- with(dat,out)        lm.res.correct <- lm(y~X1*X2)        my.data <- lm.res.correct$model #load the data from the lm object        newdata <- expand.grid(X1=c(-1,0,1),X2=c(-1,0,1))        #No need to rename newdata as it matches my lm object already        newdata$Y <- predict(lm.res.correct,newdata)      Notably, adjusting my formula to include as.vector() does not solve the problem with my attempt to use predict() directly with newdata.


title: Computing similarity in LSA

text: I have a question regarding Latent Semantic Analysis (LSA) introduced in http://lsa.colorado.edu/papers/JASIS.lsi.90.pdf . On page 14, schemas for calculating similarities between different types of objects are explained. My question is - what's the point of $S^2$ in these formulas, for example in calculating documet-document similarity by $D\times{S^2}\times{D^T}$ - why not simply use $D\times{D^T}$. There is some explanation given below that 'the positions of the points are the same except that each of the axes has been stretched or shrunk in proportion to the corresponding diagonal element of $S'$. I don't really get it, so I will be very grateful for some clarification. How, in fact, should I calculate the cosine similarity between two documents in this representation? Dot product is not enough, since rows of $T$ and $D$ are not of unit length, so I guess it must be scaled by something - but what exactly? By the product of both vector lengths, as the cosine similarity definition says? Or should I somehow use the $S^2$ in here as well? I will be very grateful for a clear explanation.


title: R package for symbolic data analysis

text: Most methods for symbolic data analyis are currently implemented in the SODAS software. Are there any R packages for symbolic data except clamix and clusterSim?


title: Help understand kNN for multi-dimensional data

text: I understand the premise of kNN algorithm for spatial data. And I know I can extend that algorithm to be used on any continuous data variable (or nominal data with Hamming Distance). However, what strategies are used when dealing with higher dimensional data? For example, say I have a table of data (x[1], x[2], x[3], ..., x[n]) and I want to build a set of classifiers to predict one of those columns (say x[n]). Using kNN algorithm I would pick any two columns from the remaining columns (x[1]-x[n-1]) to train against. So say I could pick x[1] and x[2] and build a classifier off those. Or I could pick x[1] and x[4], or I could pick x[5] and x[8], etc. I could even pick just a single column and build a classifiers off that, or 3 columns and build a classifiers off that. Is there an advantage to using higher dimensions (2D, 3D, etc) or should you just build x-1 single dimension classifiers and aggregate their predictions in some way? Since building all of these classifiers from all potential combinations of the variables would be computationally expensive. How could I optimize this search to find the the best kNN classifiers from that set? And, once I find a series of classifiers what's the best way to combine their output to a single prediction? Voting might be the simplest answer to this question. Or weighting each vote by error rates from the training data for each classifier. How do most implementations apply kNN to a more generalized learning?


title: Can I use multiple regression when I have mixed categorical and continuous predictors?

text: It looks like you can use coding for one categorical variable, but I have two categorical and one continuous predictor variable. Can i use multiple regression for this in SPSS and if so how? thanks!


title: ANOVA repeated measures when data are non normal and non homogenous

text: > **Possible Duplicate:**   >  ANOVA assumption normality/normal distribution of residuals I have a dataset storing measurements of seed area, colour and roundness over 5 time points broken down by 2 treatment groups. There are 12 replicates for each gruop over the time points. I wanted to do a normal ANOVA (Area ~ group*time + Error(id)) and a repeated measures anlysis and a Tukey test (time and group) but my Area in my data is not normal but it is homegenous, colour is the other way round and roundess is neither normal nor homogenous. I used boxcox to transform the data but I still get the same results. What do you recommend to analyse this data? Thanks in advance? Luana


title: Weight a rating system to favor items rated highly by more people over items rated highly by fewer people?

text: Thanks in advance for bearing with me, I am not a statistician of any kind and don't know how to describe what I'm imagining, so Google isn't helping me here... I'm including a rating system in a web application I'm working on. Each user can rate each item exactly once. I was imagining a scale with 4 values: "strongly dislike", "dislike", "like", and "strongly like", and I had planned on assigning these values of -5, -2, +2, and +5 respectively. Now, if every item was going to have the same number of ratings, then I would be quite comfortable with this scoring system as clearly differentiating the most liked and least liked items. However, the items will not have the same number of ratings, and the disparity between the number of votes on different photos may be quite dramatic. In that case, comparing the cumulative scores on two items means that an old item with a lot of mediocre ratings is going to have a much higher score than an exceptional new item with many fewer votes. So, the first obvious thing I thought of us to take an average... but now if an item has only one rating of "+5" it has a better average than an item that has a score of 99 "+5" ratings and 1 "+2" rating. Intuitively that isn't an accurate representation of the popularity of an item. I imagine this problem is common and you guys don't need me to belabor it with more examples, so I'll stop at this point and elaborate in comments if needed. My questions are:   1. What is this kind of problem called, and is there a term for the techniques used to solve it? I'd like to know this so I can read up on it.   2. If you happen to know of any lay-friendly resources on the subject, I'd very much appreciate a link.   3. Finally, I'd appreciate any other suggestions about how to effectively collect and analyze this kind of data. Thank you very much for your time!


title: Help relating survey answers to hypothesis

text: I am statistically challenged and trying to put my survey results into proper words. I am surveying if more games result in any type of aggression or not. I am using ANOVA and I am not even sure if I am required to use it ... Questions like 1) Do you feel aggressive when you lose a game a. Yes (value 1) b. No (value 2) c. Sometimes (value 3) 2) Do you get obsessed with a game if you lose a couple of games in a row a. Yes, I have to defeat him once (value 1) b. No (value 2) c. Sometimes (value 3) Now these two variables define the aggression combined for each specific age group. The question is, how would I relate these answers to my hypothesis? If I do ANOVA, I can tell if one group is significantly different from another. If they are, then these two variables don't define the level of aggression but if they aren't then it means these two questions relate to each other (correct me if I am wrong here). That being said, how do these two questions actually relate to my hypothesis? Do I have to do some correlation?


title: Choice of a less true model over the truer model if it predicts better and my purpose is prediction

text: Sometimes a less true model predicts better than a truer model (When will a less true model predict better than a truer model?). So should I choose a less true model over the truer model if it predicts better and my purpose is prediction? Similarly, should I choose a model which violates the assumptions (e.g normality or homoscedasticity of the residuals in linear regression) over the model which satisfies all of the assumptions if it predicts better and my purpose is prediction?


title: Paired t-test, Kruskal-Wallis rank sum test and chisquare test. Which test should I choose?

text: I observed moth oviposition in orchard for 4 moths (400 clusters per week). The moth oviposit on leaf and apple. I want to compare the means between oviposit on leaf and apple. I suppose to use paired t-test, however, the result of shapiro.test showed my data do not fit normality. I was wondering if I should use Kruskal-Wallis rank sum test for this case? Moreover, I did some treatment on apples to prevent the ovipostion. I observed eggs on treated and untreated clusters weekly (treated and untreated apple trees were in the same orchard). Then I got a table like:                              Treated  untreated     Egg Observed       40       1852     Non-observed     6000      12600      Can I use a chisquare test in this case?


title: Testing for Poisson process

text: I have some discrete times of events and I would like to do a test to see if they are likely to have come from a homogeneous Poisson process. From www.stat.wmich.edu/wang/667/classnotes/pp/pp.pdf I see " REMARK 6.3 ( TESTING POISSON ) The above theorem may also be used to test the hypothesis that a given counting process is a Poisson process. This may be done by observing the process for a fixed time t. If in this time period we observed n occurrences and if the process is Poisson, then the unordered occurrence times would be independently and uniformly distributed on (0, t]. Hence, we may test if the process is Poisson by testing the hypothesis that the n occurrence times come from a uniform (0, t] population. This may be done by standard statistical procedures such as the Kolmogorov-Smirov test." However I don't quite understand what to do in practice. Say my times are               [1, 7, 18, 22, 41, 43, 66, 73, 86, 92]      and the time interval I chose was from $1$ to $100$. How exactly do I do the Kolmogorov-Smirov test in this example?


title: $\chi^2$ test of homogeneity for three-way contingency table

text: I have a three-way contingency table in which the marginal totals for two sides are fixed and for the third are random. I'm wondering how to perform a chi-square test for homogeneity for such a three-way contingency table. **Example** ![enter image description here](http://i.stack.imgur.com/r3QGt.png) Assume both Trt and Gender marginal totals are fixed. Thanks


title: Finding Optimal threshold from ROC curve (SVM classifier in Matlab)

text: I plotted the ROC curve for the SVM classifier performance in Matlab. I would need to find the optimal threshold value from the ROC curve and need to set the SVM at that particular threshold for future classification. Is there any Matlab function that can be used to automatically calculate the optimal threshold value (where sensitivity & specificity is max or above 0.85) and set the SVM at with the same threshold? Any help on this would really be appreciated. Thanks in advance.


title: Bootstrapped parameter and fit estimates with non-normality for structural equation models

text: ### Context: Within the context of structural equation modelling, I have non-normality according to the Mardia test but univariate indices of skewness and kurtosis are less than 2.0. ### Questions:   * Should parameter estimates (coefficient estimates) be evaluated using bootstrapping (1000 replicates) with bias-corrected methods?    * In place of the traditional chi-square test, should the Bollen-Stine bootstrapped version be used?


title: Different score range when calculating area of under curve in ROC curves

text: I have two classifiers which try to classify the same data sets. In order to check the efficiency of the classifiers I intend to plot the curves and calculate the AUC value. The concern is that one of the classifiers always produce positive values as the scores (actually it is the p value) but the other classifier produces both positive and negative values as the scores(it has its own producing score method rather than producing pvalues). Should I take the absolute values for the second scores when plotting the roc curves and calculating the AUC or I can use scores as they are (both positive and negative). When I make the absolute values for the second classifier to plot the curves the AUC amount decreases. Thank you in advance.


title: Advanced regression in SPSS

text: I have data of several school classes (with several pupils). I want to calculate expected/predicted values of class A for some items using a regression (a multiple linear regression in that special case) for which I want to estimate the regression parameters/weights using all classes except A. Lets say for an item $I$, I want to estimate $f_{model}$ as $f_{model}\left( x_0, \dots, x_n \right) = I + \varepsilon$ by using data $x_0, \dots, x_n$ from all pupils $p\in \\{\mbox{pupils from all classes}/(p\in A)\\}$ where $A$ is the set of pupils in the school class $A$. Having $f$ I'd like to calculate $Y_{predicted} = f_{model}\left( y_0, \dots,y_n \right)$ for the data $y_0, \dots, y_n$ from all pupils $p' \in A$. Although I know how to do that for one case (write down the paramters and calculate the regression manually), I'd like to know whether there is an automated way to do that, e.g. a function or method or a checkbox to check. Even the name of an extra package would be sufficiant, probably my university already has a licence for it.


title: Meaning of a t-test comparing regression output from two samples

text: Suppose two samples A and B of size N = 20 are taken from a population of pairs $(Y_i, X_i)$, and separate OLS regressions are calculated from each sample for the model: $$Y_i=\beta_1+\beta_2X_i+\varepsilon_i$$ Considering say the slope coefficient, this yields two estimates $\hat\beta_{2A}$ and $\hat\beta_{2B}$ and their associated standard errors $s_{2A}$ and $s_{2B}$. It is possible then to apply a t-test of a null hypothesis that I provisionally state as follows: > $H_0$: The difference between the mean of the sampling distribution of the > estimate of $\beta_2$ associated with sample A and the mean of the sampling > distribution of the estimate of $\beta_2$ associated with sample B is less > than $\Delta$. Here $\Delta$ is the size of difference that would be of practical concern given the purpose of the model (I've been guided here by the answers to this question). Since the standard errors are likely to differ, the appropriate test appears to be Welch’s t-test, the test statistic being (assuming $\hat\beta_{2A} > \hat\beta_{2B}$): $$t=\frac{\hat\beta_{2A}-\hat\beta_{2B}-\Delta}{\sqrt(s_{2A}^2+s_{2B}^2)}$$ I can see how to do the calculation, and see in general terms that if the null hypothesis were to be rejected, that might suggest departures from the assumptions of the classical linear regression model or non-randomness in sample selection. However, I am puzzled as to what _exactly_ is tested by such a t-test. **Question** : How could my formulation of $H_0$ be improved to indicate more precisely what this application of a t-test actually tests? As formulated, it sems to assume that it makes sense to refer to the sampling distribution _of a particular sample_. This seems wrong because a sampling distribution (of a statistic such as the mean) is a property of _repeated_ samples, not of just one sample. If the two samples had been drawn from distinct sub-populations, then it would make sense to refer to the sampling distributions of the means for each sub-population. But this is not the case here: both samples are drawn from the whole population.


title: How to apply Neural Network to time series forecasting?

text: I'm new to machine learning, and I have been trying to figure out how to apply neural network to time series forecasting. I have found resource related to my query, but I seem to still be a bit lost. I think a basic explanation without too much detail would help. Let's say I have some price values for each month over a few years, and I want to predict new price values. I could get a list of prices for the last few months, and then try to find similar trends in the past using K-Nearest- Neighbor. I could them use the rate of change or some other property of the past trends to try and predict new prices. How I can apply neural network to this same problem is what I am trying to find out.


title: Logistic regressioin - contradictory results?

text: I am using the glmer method in lme4 to build a mixed effect logistic regression model, the model is as follows:                textbook.usage.glm <- glmer(textbook.usageSession ~ session.week * condition.player  + (session.week|group.name),family="binomial",data=dfAN)      When I run the summary                 summary(textbook.usage.glm)      I got the following results:                     Generalized linear mixed model fit by the Laplace approximation            Formula: textbook.usageSession ~ session.week * condition.player + (session.week |      group.name)            Data: dfAN            AIC   BIC logLik deviance           104.9 123.1 -45.45     90.9           Random effects:           Groups     Name         Variance Std.Dev. Corr              group.name (Intercept)  0.45445  0.67413                      session.week 0.17793  0.42181  -0.484            Number of obs: 99, groups: group.name, 6               Fixed effects:                                       Estimate Std. Error z value Pr(>|z|)            (Intercept)                            4.6851     3.4607   1.354   0.1758            session.week                          -1.2304     0.8017  -1.535   0.1248            condition.playerDEFAULT                4.7845     4.2266   1.132   0.2576            session.week:condition.playerDEFAULT  -3.1617     1.4026  -2.254   0.0242 *          ---          Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1               Correlation of Fixed Effects:             (Intr) sssn.w c.DEFA         session.wek -0.976                       cnd.DEFAULT -0.787  0.765                s.:.DEFAULT  0.515 -0.506 -0.926      From the above statistics, we can see only one significance for the interaction between session.week and condition.player. But there is no significance on the main effect (condition.player) alone. However, if we further build three incremental models as follows, and run anova to test the significance of each variable                   textbook.usage.glm.2 <- glmer(textbook.usageSession ~ session.week +  condition.player + (session.week|group.name),family="binomial",data=dfAN)         textbook.usage.glm.1 <- glmer(textbook.usageSession ~ session.week + (session.week|group.name),family="binomial",data=dfAN)         textbook.usage.glm.0 <- glmer(textbook.usageSession ~ 1 + (session.week|group.name),family="binomial",data=dfAN)        anova(textbook.usage.glm.0,textbook.usage.glm.1,textbook.usage.glm.2,textbook.usage.glm)      We get                 Data: dfAN       Models:       textbook.usage.glm.0: textbook.usageSession ~ 1 + (session.week | group.name)       textbook.usage.glm.1: textbook.usageSession ~ session.week + (session.week | group.name)       textbook.usage.glm.2: textbook.usageSession ~ session.week + condition.player + (session.week |        textbook.usage.glm.2:     group.name)       textbook.usage.glm: textbook.usageSession ~ session.week * condition.player + (session.week |        textbook.usage.glm:     group.name)                      Df    AIC    BIC  logLik   Chisq Chi Df Pr(>Chisq)           textbook.usage.glm.0  4 126.57 136.95 -59.284                                     textbook.usage.glm.1  5 125.43 138.41 -57.716  3.1376      1    0.07651 .         textbook.usage.glm.2  6 109.48 125.05 -48.738 17.9548      1  2.262e-05 ***       textbook.usage.glm    7 104.91 123.07 -45.452  6.5716      1    0.01036 *         ---       Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       Here, we see, by adding condition.player in the model 2, it is significant as well, in addition to significance to the interaction effect. How can I report the result? Is the condition.player variable significant or not? I've seen many logistic regression examples with the first approach, what is the difference between that and the second approach where I compare models?


title: What are some ways to graphically display non-normal distributions in Excel?

text: The ultimate goal is to show users, at a glance, if their data is normally distributed. The first attempt is a kludge that plots the data in a frequency graph. Then, the observed mean and standard deviation are used to build a "normal curve" graph. The frequency chart is laid over the the normal-curve chart and put next to some key statistics. The frequency chart also colors positive bins green and negative bins red. It looks like this: ![enter image description here](http://i.stack.imgur.com/vPNKj.png) I understand the fallacy of this approach, but for now it's practical. What is a better way to approach this issue?


title: How to combine data from 5 surveys from the same population spanning 10 years

text: I have results from 5 surveys each 2 years apart and let us assume that no subjects are selected in more than one survey. The sampling method used in these surveys are biased and I have sampling weights calculated(with respect to the population) for each data point in each study. The question is, how would I be able to combine the 5 datasets and have the weights recalculated so as to obtain one giant dataset for analysis on this population? Also, what should I do if subjects appear in more than one survey? # Updates/Further Elaboration: thank you @user30523, here are some more infomation that might be useful: Suppose I wish to find out the estimated distribution of height across the population using these 5 datasets. In some data, younger people are oversampled because of the location where the survey are conducted. Let's assume the weights are calculated with respect to their age. Eg. assuming 2% of the population are 15 years old, and the location of the survey is at a mall where 15-year-olds made up 5% of all shoppers, then sampling weight for an subject aged 15 in that survey would be calculated as 0.02 / 0.05 = 0.4. For simplicity, each person in the mall has equal chance of being surveyed and all participants complied when asked. Given that 5 surveys are conducted in 5 different malls and each has their set of weights calculated in the same way, how would I then be able to combine all 5 datasets and recalculate the sampling weights? P.S: I'm new to the topic on sampling weights so do correct me if I have made errors in the way I have calculated the weights.


title: Non-normality and heterogeneity in ANCOVA

text: I am performing an ANCOVA. My dependent variable is the number of days a student is absent from school. As one can imagine, this variable is not normal, so I performed a log10 transformation. When I run the analysis with the transformed variable, I find no significant relationship with my variable of interest. However, when I run it with the original dependent variable the relationship is significant. Also, the homogeneity assumption is violated with the original dependent variable, but not with the transformed variable. Given the violations, I was under the impression that statistical power should be reduced. Why is it that I only find significant results when these assumptions are violated? Can someone help explain this to me? I should also mention that my sample sizes are quite unequal and in this case the smaller cell has a smaller variance on the dependent variable.


title: Suggest an exhaustive procedure that will find a separating vector for linearly separable pattern in a finite number of steps

text: I read in a pattern classification text, that if we consider weight vectors whose components are integer valued, the perceptron procedure would terminate in a finite number of steps. What is the intuition and theory behind this? **EDIT:** It is a homework problem from the text book Pattern classification by Duda, Hart, Stork, and I just wanted a hint, so posted the above query. Here is the actual homework problem: > Let ${y_{1},...,y_{n}}$ be a finite set of linearly separable samples in d > dimensions. Suggest an exhaustive procedure that will find a separating > vector in a finite number of steps. (You might wish to consider weight > vectors whose components are integer valued.)


title: Adding correlation coefficients of time series

text: I have computed correlation coefficients for 90 day increments of a time series (i.e., one coefficient for days 0-90, 91-180, etc.). It was computationally very expensive to compute these and I would like to analyze correlations over 6 month periods, 1 year periods etc. Is it possible to some how to combine these coefficients? I.e., can I get the correlation over 180 days as 0.5*c1+0.5*c2? I know in the general case this cannot be done, but because the coefficients are computed from the same number of samples, is it possible?


title: Zero-inflated count models in R: what is the real advantage?

text: For analysing zero-inflated bird counts I'd like to apply zero-inflated count models using the R package pscl. However, having a look at the example provided in the documentation for one of the main functions (?zeroinfl), I begin doubting what's the real advantage of these models. According to the sample code given there, I calculated standard poisson, quasi-poisson and negative bionomial models, simple zero-inflated poisson and negative binomial models and zero-inflated poisson and negative-binomial models with regressors for the zero component. Then I inspected the histograms of the observed and the fitted data. (Here's the code for replicating that.)               library(pscl)     data("bioChemists", package = "pscl")          ## standard count data models     fm_pois  <- glm(art ~ .,    data = bioChemists, family = poisson)     fm_qpois <- glm(art ~ .,    data = bioChemists, family = quasipoisson)     fm_nb    <- glm.nb(art ~ ., data = bioChemists)          ## with simple inflation (no regressors for zero component)     fm_zip  <- zeroinfl(art ~ . | 1, data = bioChemists)     fm_zinb <- zeroinfl(art ~ . | 1, data = bioChemists, dist = "negbin")          ## inflation with regressors     fm_zip2  <- zeroinfl(art ~ fem + mar + kid5 + phd + ment | fem + mar + kid5 + phd +                           ment, data = bioChemists)     fm_zinb2 <- zeroinfl(art ~ fem + mar + kid5 + phd + ment | fem + mar + kid5 + phd +                           ment, data = bioChemists, dist = "negbin")          ## histograms     breaks <- seq(-0.5,20.5,1)     par(mfrow=c(4,2))     hist(bioChemists$art,  breaks=breaks)     hist(fitted(fm_pois),  breaks=breaks)     hist(fitted(fm_qpois), breaks=breaks)     hist(fitted(fm_nb),    breaks=breaks)     hist(fitted(fm_zip),   breaks=breaks)     hist(fitted(fm_zinb),  breaks=breaks)     hist(fitted(fm_zip2),  breaks=breaks)     hist(fitted(fm_zinb2), breaks=breaks)!      ![Histogram of observed and fitted data](http://i.stack.imgur.com/BUYsT.png) I can't see any fundamental difference between the different models (apart from that the example data don't appear very "zero-inflated" to me...); actually none of the models yields a halfway reasonable estimation of the number of zeros. Can anyone explain what's the advantage of the zero-inflated models? I suppose there must have been a reason to choose this as the example for the function.


title: Logistic input with Gaussian noise

text: Both the logistic function and standard deviation are usually denoted $\sigma$. I'll use $\sigma(x) = 1/(1+\exp(-x))$ and $s$ for standard deviation. I have a logistic neuron with a random input whose mean $\mu$ and standard deviation $s$ I know. I hope the difference from the mean can be approximated well by some Gaussian noise. So, with a slight abuse of notation, assume it produces $\sigma(\mu + N(0,s^2))=\sigma(N(\mu,s^2))$. **What is the expected value of $\sigma(N(\mu,s^2))$?** The standard deviation $s$ might be large or small compared with $\mu$ or $1$. A good closed form approximation for the expected value would be almost as good as a closed form solution. I don't think a closed form solution exists. This can be viewed as a convolution, and the characteristic function for the logistic density is known ($\pi t ~\text{csch} ~\pi t$), but I'm not sure how much that helps. The inverse symbolic calculator was unable to recognize the density at $0$ of the convolution of the logistic distribution's density and a standard normal distribution, which suggests but does not prove that there is no simple elementary integral. More circumstantial evidence: In some papers on adding Gaussian input noise to neural networks with logistic neurons, the papers didn't give closed form expressions either. This question arose in trying to understand the error in the mean field approximation in Boltzman machines.


title: Why does a mixed design using R's aov() need the between subject factors specified more than once?

text: An R cookbook http://wiki.stdout.org/rcookbook/Statistical%20analysis/ANOVA/ has an example of using aov() for mixed design ANOVAs. I'll copy it here:               data <- read.table(header=T, con <- textConnection('      subject sex   age before after            1   F   old    9.5   7.1            2   M   old   10.3  11.0            3   M   old    7.5   5.8            4   F   old   12.4   8.8            5   M   old   10.2   8.6            6   M   old   11.0   8.0            7   M young    9.1   3.0            8   F young    7.9   5.2            9   F   old    6.6   3.4           10   M young    7.7   4.0           11   M young    9.4   5.3           12   M   old   11.6  11.3           13   M young    9.9   4.6           14   F young    8.6   6.4           15   F young   14.3  13.5           16   F   old    9.2   4.7           17   M young    9.8   5.1           18   F   old    9.9   7.3           19   F young   13.0   9.5           20   M young   10.2   5.4           21   M young    9.0   3.7           22   F young    7.9   6.2           23   M   old   10.1  10.0           24   M young    9.0   1.7           25   M young    8.6   2.9           26   M young    9.4   3.2           27   M young    9.7   4.7           28   M young    9.3   4.9           29   F young   10.7   9.8           30   M   old    9.3   9.4     '))     close(con)      Then reshape it:               library(reshape2)          # Make sure subject column is a factor     data$subject <- factor(data$subject)           # Convert it to long format     data.long <- melt(data, id = c("subject","sex","age"), # keep these columns the same                   measure = c("before","after"),       # Put these two columns into a new column                   variable.name="time")                # Name of the new column          # subject sex   age   time value     #       1   F   old before   9.5     #       2   M   old before  10.3     #...      Now analyze using a mixed anova:               aov.after.age.time <- aov(value ~ age*time + Error(subject/time), data=data.long)     summary(aov.after.age.time)      But when there are more than two predictor variables, the R examples show that the between subject factors are added again after the error term:               #e.g., from R cookbook     #aov.bww <- aov(y ~ b1*b2*w1 + Error(subject/(w1)) + b1*b2, data=data.long)          # which would translate in our case as:     aov.bww <- aov(value ~ sex*age*time + Error(subject/time) + sex*age, data=data.long)     summary(aov.bww)      But why is b1*b2, or in our case sex*age, specified twice? It doesn't seem to make a difference when we remove them after the Error() term:               aov.bww2 <- aov(value ~ sex*age*time + Error(subject/time), data=data.long)     summary(aov.bww2)      Can anyone explain why the examples have those extra terms? The R manual just has this example, where the between factors are not specified twice:               # fm <- aov(yield ~ v + n*p*k + Error(farms/blocks), data=farm.data)      * * * **Edit:** I have checked the references from the R Cookbook and found other web sites also specify the terms twice in their mixed design examples. See here: http://www.personality-project.org/R/r.anova.html where they have the example:               aov.ex5 = aov.ex5 = aov(Recall ~ (Task*Valence*Gender*Dosage) +     Error(Subject/(Task*Valence)) + (Gender*Dosage), data.example5 )      and see here http://www.statmethods.net/stats/anova.html with their example:               # Two Within Factors W1 W2, Two Between Factors B1 B2     fit <- aov(y ~ (W1*W2*B1*B2) + Error(Subject/(W1*W2)) + (B1*B2),     data=mydataframe)      Which is presumably where the cookbook got their info from.


title: A logistic problem about decision theory

text: Belows is the question and the solution for part c which is the part that i don't understand. Can someone explains to me? I don't quite get how it gets $3\over 7$ and why does it needs to?Any hint about that??? Transcribed from two images: Dora Bender and Joey Mark graduated together at Business School. Both hope to open a coffee shop around the corner in the financial center. Using a system approach, Dora and Joey have identified three strategies. Strategy $1$ is to invest in a fairly expensive coffee brewer with high-quality coffee bean imported from Italy. In a favorable market, they should be able to obtain a net profit of $\$12,000$ in a month. If the market is unfavorable, they can lose $\$8,000$. Strategy $2$ is to purchase a less expensive brewer. In a favorable market, they could earn a return of $\$8,000$, whereas in an unfavorable market, they would incur a loss of $\$5,000$. The third strategy is to do nothing. Basically, Joey is a risk taker, but Dora is a risk avoider. (a) What type of decision procedure should Joey use? What would Joey's decision be? (b) Sketch the utility curve vs. monetary value for Dora and Joey. (c) During the decision making period, Dora and Joey have been approached by a marketing firm. PriceAd Marketing, who claimed that they could provide technical analysis with certainty for the market condition for the coming year. Using the Expected Value of Perfect Information (EVPI) technique, determine the most Joey and Dora would be willing to pay for PriceAd? $\begin{array}{} \text{(c) } &0 \le p \le \frac37\quad&\text{EVPI} = 5000 - 1000 p\\\ &\frac37 \lt p \lt 1\quad&\text{EVPI} = 8000 - 8000 p\\\ &p=0.5\quad&\text{EVPI} = \$4000 \end{array}$


title: Calculating some kind of confidence or error rate for a set of binary data

text: Newbie question here. I'm calculating error rates from software test results. Basically for any particular test run it's going to pass or fail, in this case, because of inherent race conditions in the software in question. We're running about ~50k overall tests a week (some particular tests might be run 4000 times a week, some 200, with all sorts of numbers in between), so I'm starting to quantify the results of these into failure rates. However, clearly the answer from failing 3 out of 200 tests is a less certain value that failing 30 out of 2000. So given a series of binary data, with no known baseline, is there a standard method using the sample size to provide some error rate? The other challenge that this data has is we are doing rolling binning. So calculating the failure rates every day. Which means that at some point early in the day there are only a few events. It would be good to be able to quantify the difference in confidence between 4/4 success results vs. 100/100 success results. Part of my thinking was that you could quantify this as 1/n (n being population size) as you could never know the answer better than the "resolution" you have. However, that's just instincts kicking in, with no statistical backing in there.


title: Lasso on Negative Binomial Regression Model

text: Is there anyway that I can perform LASSO with Negative Binomial Regression on R? I am performing a negative binomial regression on my dataset because the data are too dispersed to impose poisson regression. Meanwhile, I am also facing some multicollinearity problem. I already tried using `glmnet` with `family = poisson`, but the data is not fitting very well (for both alpha = 0 and alpha = 1)...I honestly don't know what to do to analyze this big mess of data :/ thank you EDIT: here is variance-covariance table of the negative binomial fit                      8.392729e+18  1.239178e+06  -3.624090e+05  1.896258e+17  -3.702521e+17            1.239178e+06  1.119052e-04   5.201989e-06 -1.877590e+05  -2.558095e+05           -3.624090e+05  5.201989e-06   5.179343e-06 -8.021543e+04  -1.436381e+05            1.896258e+17 -1.877590e+05  -8.021543e+04  2.193290e+17   6.413947e+16           -3.702521e+17 -2.558095e+05  -1.436381e+05  6.413947e+16   2.142183e+17


title: Probability that a given Normal Distribution is Maximum among others

text: You are given the mean and standard deviations of N normal distributions x1,x2...xn What is the probability that x1 is maximum? ie. Find P(x1>x2,x3..xn) How do I go about solving this? x1,x2,x3 etc are independent. Any help will be appreaciated, Thanks!


title: Permutation test for factor analysis

text: We have a survey instrument and are interested in assessing dimensionality of it. Looking at plots of multidimensional scaling, it appears as though there are, perhaps, 3 distinct dimensions to the survey since there are 3 seemingly well defined clusters of responses. When I perform a Scree plot, I observe there are 7 dimensions for which the eigenvalues are greater than one, so the Kaiser-Guttman indicates this anticonservative high dimensionality. Doing a parallel analysis by generating independent random normal values of the same shape as the original matrix of responses gives a much more conservative 4 dimensions before eigenvalues become consistent with those of randomly generated data. However, the parallel analysis doesn't factor in the coding and format of the data, and differing distributions of possible covariances under the null hypothesis. It makes sense to do a permutation test of the original response data by permuting values in each column. This maintains the same mean and standard deviation for these values, but finds a sampling distribution for the correlation matrix under a null hypothesis of completely independent data. Has this been explored before? Is there a name for such a test? Are there caveats to this approach?


title: Question about computing Bayes Error - with or without loss function?

text: I am new to Bayesian Decision Theory and don't understand the following concept: So from what I understood, the Bayes error is used to report the performance of a Bayes classifier in terms of the probability of making and error. From the conditional error probabilities ![conditional error](http://sebastianraschka.com/_my_resources/images/equations/cond_error.png) we can obtain the total probability of making an error (the probability to mis-classify). ![error](http://sebastianraschka.com/_my_resources/images/equations/error.png) Now, if my Bayes classifier was designed to minimize the overall risk, I have a loss function that gives penalties to certain decisions. ![conditional risk](http://sebastianraschka.com/_my_resources/images/equations/cond_risk.png) ![overall risk](http://sebastianraschka.com/_my_resources/images/equations/overall_risk.png) So, if my classifier includes such a loss function when I optimize my classifier for minimum overall risk, shouldn't be the Bayes error also include the loss function term? Hope you can help me here, because I think I am missing something here ...       **EDIT:** I'll try to express my problem using an 2D-classification problem: Let's assume I have two pdfs (e.g., p(x|c1) and p(x|c2) ) with slight overlap. And mis-classifying a pattern as c2 where it truly belongs to c1 is more costly than vice-versa. In this case I would assign a higher loss to **"classify pattern x as c1 when it is truly c2"** than **"classify pattern x as c2 when it is truly c1"** in order to calculate and minimize the overall risk. I would therefore increase the probability to classify a pattern x as c2 over c1 due to the minimum risk optimization. Isn't this something I have to also include in p(error)?


title: Weibull distribution parameters $k$ and $c$ for wind speed data

text: I have 5 years worth of mean wind speed meteorological data. I want to compare the observed wind speed data against a Weibull distribution on a graph of **the probability density distribution.** I have read about the different methods to find parameters $k$ (shape) and $c$ (scale) of a Weibull distribution, and I found that the maximum-likelihood method is better but I am confused. What should my steps be?


title: Correlation among repeated measures - I need an explanation

text: I am asked in G*Power to enter the _Correlation among repeated measures_. I have repeated an experiment with the same subjects under 3 conditions (Set1, Set2 and Set3). I calculate the correlation this way:               lala=cbind(         Set1.Weber,       Set2.Weber,       Set3.Weber     )     cor(lala)      And get                          Set1.Weber Set2.Weber Set3.Weber     Set1.Weber  1.0000000  0.3683676  0.1283023     Set2.Weber  0.3683676  1.0000000 -0.0959547     Set3.Weber  0.1283023 -0.0959547  1.0000000      as a result. Is that correct? In G*Power there is only **one** parameter for the correlation. Should I use the average correlation. In the example:               (0.37+0.13+0.10)/3      ? Here is a screenshot: ![enter image description here](http://i.stack.imgur.com/mUcVm.png) **Edit:** As a response to the very good answers so far I want to specify the question: * * * Do I understand that correct that the ANOVA assumes that the correlation should be equal in the matrix and not like in my case (3 different values)? So it maybe not super correct but I can use an average as a poormans solution?! I could use (|a|+|b|+|c|)/3 or sqrt(a²+b³+c²)?! * * *


title: What are pitfalls of bootstrapping on random sample of master data?

text: Will I obtain seriously biased results if I use bootstrapping on a subsample of a larger dataset? Rather than drawing 100 bootstrap samples from a dataset of 50 million + records, which could hog server resources, I'm thinking of first drawing a 5% random sample, w/ replacement, of records from the master dataset. Then build bootstrapped 95% confidence intervals for several statistics of interest using the 5% sample. So in essence I'll be bootstrapping from a bootstrap sample. I'm working with health data which is clustered according to provider and episode of care, hence drawing a proper representative 5% sample will require some care. Some earlier advice has pointed me in the right direction ( Proper bootstrapping technique for clustered data? ).


title: Derivative of $H$ with respect to $W$ when performing generalized linear squares

text: I am trying to solve a generalized linear squares model with the following form: $\hat{Y}= X(X'\Omega^{-1}WX)^{-1}X'\Omega^{-1}WY $ $ H= X(X'\Omega^{-1}WX)^{-1}X'\Omega^{-1}W $ $ \Omega$ is the co-variance matrix. $ W$ is a diagonal matrix that weights a given observation $w_i$ and is normally $I$. Occasionally a single or multiple diagonal entries in $W$ are updated. How should a change in $W$ update $H$ without having to recompute $H$ entirely? This question is related to the following: Derivative of a linear model. Thanks for the input.


title: Constrained Maximum Lilkihood Estimation

text: Given the joint PDF $h_{X,Y}(x,y)$ of the correlated variables $X$ and $Y$, and given that the joint PDF is a function of a parameter $\lambda$ such that $|\lambda|\le 1$, I need to find the maximum likelihood estimator for the parameter $\lambda$. The Joint PDF is given by $$h_{X,Y}(x,y) = f_X(x)f_Y(y) \left[1+\lambda\left(2F_X(x)-1\right)\left(2F_Y(y)-1\right)\right],\quad |\lambda|\le 1 $$ The marginal PDFs $f_X(x)$ and $f_Y(y)$ are known as well as the marginal CDFs $F_X(x)$ and $F_Y(y)$. So how can I apply constrained maximum likelihood estimation to obtain the estimated value of $\lambda$?


title: Finding Jeffreys prior for the multivariate normal distribution. Conditional probability algebra

text: For the multivariate normal model, find Jeffrey's prior of the form $p_j(\theta, \Sigma|y1,...,yn)$, $p_j(\theta|\Sigma, y_1,...,y_n)$, and $p_j(\Sigma|y_1,...,y_n)$. Attempt: I know that $p_j(\theta, \Sigma|y_1,...,y_n) \propto p(\theta, \Sigma)xp(y_1,...,y_n|\theta, \Sigma)$, where $p_j(\theta, \Sigma) \propto |\Sigma|^{-(p+2)/2}$ and $p(y_1,...,y_n|\theta, \Sigma) = \frac{1}{(2 \pi)^{p/2}} \frac{1}{\sqrt |\Sigma|} exp(-\frac{1}{2}(X - \theta)^T \Sigma^{-1} (X - \theta))$ This leads me to think that $p_j(\theta, \Sigma|y_1,...,y_n) = |\Sigma|^{-(p+2)/2} \frac{1}{(2 \pi)^{p/2}} \frac{1}{\sqrt |\Sigma|} exp(-\frac{1}{2}(X - \theta)^T \Sigma^{-1} (X - \theta))$. I am not sure how to arrange it into $p_j(\theta, \Sigma|y_1,...,y_n)$ into $p_j(\theta|\Sigma, y_1,...,y_n)$ and $p_j(\Sigma|y_1,...,y_n).$ For my last question, why cannot $p_j(\theta, \Sigma)$ be a probability density for $(\theta, \Sigma).$


title: How can we measure between-experimenter reliability for experimental data?

text: **Description of my situation:** I will run an experiment in a few weeks with a total sample size of _N_ = 70. Ideally I would prefer to have only one experimenter to collect all the data but I have to do with seven experimenters. So I will randomly assign the 70 participants to the 7 experimenters in order to control a potential experimenter effect. **My question:** Is there any way to compute a between-experimenter reliability index for my data as in the case of inter-rater reliability index? Thanks in advance for your answers.


title: Interpretation of results of power analysis with regard to the null hypothesis

text: I have recently read a paper, in which the authors conducted a $2*3$ Between- within-ANOVA. The alternative hypothesis was formulated with regard to an interaction effect, which did not reach the significance level of $\alpha = .05$. The authors then went on to state that "the null hypothesis can be accepted for an effect size of $f^2 = .15$ because of sufficient power ($1 - \beta = .97$ for $\alpha = .05$$". As far as I know, _power analysis_ is usually conducted to compute a conditional probability, namely the probability of _getting a significant result ($H_1$)_ given _sample size ($N$)_ , _effect size (in this case $f^2$)_ and _significance level ($\alpha$)_. Mathematically, this should be expressed as $$Power = P(H_1 | N,f^2,\alpha)$$. The statement that the null hypothesis can be accepted sounds more like the _probablility that the null hypothesis is true_ to me. Is the interpretation of the authors too far-reaching? Or put in a different way: Can one conclude something about the null hypothesis from the results of a power analysis?


title: Difficulty with auto.arima function results

text: I am new to R and the ARIMA model and I am attemping to forecast 1440 values into the future using a base of roughly 5000 numbers. It is data extracted roughly every minute from a machine log(performance values). The intend to forecast 1 day into the future, which explains the 1440 values(as they are minutes). Here is my result using the following commands: datats<-c(data); arima<-auto.arima(datats); fcast<-forecast(arima, h=1440); ![forecast](http://i.imgur.com/N7aCqvx.png) The prediction begins at the flat line on the right hand side. Forecast method: ARIMA(0,1,1) Model Information: Series: datats ARIMA(0,1,1) Coefficients: ma1 -0.9373 s.e. 0.0071 sigma^2 estimated as 86737: log likelihood=-21221.46 AIC=42446.93 AICc=42446.93 BIC=42458.93 Error measures: ME RMSE MAE MPE MAPE MASE Training set 0.6506441 294.4619 196.7211 -59.85254 85.45473 0.7637028 ACF1 Training set 0.01519673 Dataset here: http://pastebin.com/92ssDExn Is the issue too little past values? To many values to be predicted? Any information or advice would be extremely welcomed, any other information required will be provided.


title: Mean and single prediction bands for custom nonlineal CDF in Mathematica

text: I have a CDF obtained with the code:                dist = ProbabilityDistribution[{"CDF",     Exp[-a Exp[-x b] - c Exp[-x  d]]}, {x, 0, Infinity},                  Assumptions -> {{0 < a < 10}, {0 < b < 1},   {0 < c < 10}, {0 < d <   1}}]        res = FindDistributionParameters[data,  dist, {{a, 3.8}, {b, 0.006}, {c, 0.08}, {d, 0.0002}}]      I can plot the confidence and prediction intervals. But for the custom CDF obtained with `ProbabilityDistribution`, I tried this:                Plot[{Func[x], bands90[x]}, {x, 0, 3000}, Filling -> {2 -> {1}}]      Without results. Does anyone know how to get the confidence and prediction bands in the case of CDF obtained with that function, as well as the ANOVA or the table that one can shows with `ParamterTable` (when fitting with `NonlinearModelFit`)? As I see in some software is possible obtain more quality information for Nonlinear fitting that for distribution or custom distribution fitting, so I trying to construct my own code. To start getting the variance-covariance matrix I followed the code found here: Standard errors for maximum likelihood estimates in FindDistributionParameters And I elimitated the last part (Sqrt[Diagonal[cov]/len]]) as I also want to abtain the covariance matrix, so my code is:                 covariance[data, dist, paramlist, mleRule] := Block[{len, infmat,       cov}, len = Length[data]; (compute negative of expected Fisher       information) infmat = -D[LogLikelihood[dist, data], {paramlist,       2}]/len /. mleRule; (invert to get asymptotic covariance) cov =       Inverse[infmat]]      In this way I got the central part of the Hat matrix, $(X^\top X)^{-1}$. But somebody can verify if I'm in the right way? If so, I have to obtain $X$ and $X^\top$, and then I have to add the variance $s^2$ to continue with the guide I found here


title: What other normalizing transformations are commonly used beyond the common ones like square root, log, etc.?

text: In the analysis of test scores (e.g., in Education or Psychology), common analysis techniques often assume that data are normally distributed. However, perhaps more often than not, scores tend to deviate sometimes wildly from normal. I am familiar with some basic normalizing transformations, like: square roots, logarithms, reciprocal transformations for reducing positive skew, reflected versions of the above for reducing negative skew, squaring for leptokurtic distributions. I have heard of arcsine transformations and power transformations, though I am not really knowledgeable about them. So, I am curious as to **what other transformations are commonly used by analysts?**


title: Constructing smoothing splines with cross-validation

text: Can someone provide me with a book or online reference on how to construct smoothing splines with cross-validation? I have a programming and undergraduate level mathematics background. I would also appreciate an overview of whether this is smoothing technique is a good one for smoothing data and whether there are any disadvantages of which a non-statistician needs to be aware.


title: Detecting strong currents in a sparse directed graph

text: I have a very large, sparse, weighted, directed graph. The structure is such that it mainly consists of strings of nodes connected with highly weighted edges. These strings can be connected by weak edges or fully unconnected.   * Is there an algorithm that will detect these highways / currents / linear clusters and group the nodes together?   * Otherwise, does this type of graph structure have a specific name that will help me research it more?


title: Intuition behind logistic regression

text: Recently I began studying machine learning, however I failed to grasp the intuition behind logistic regression. The following are the facts about logistic regression that I understand.   1. As the basis for hypothesis we use sigmoid function. I do understand why it's _a_ correct choice, however why it's the _only_ choice I don't understand. Hypothesis represents the probability that the appropriate output is $1$, therefore the domain of our function should be $[0,1]$, this is the only property of sigmoid function I found useful and appropriate here, however many functions satisfy this property. In addition, sigmoid function has a derivative in this form $f(x)(1-f(x))$, but I don't see the utility of this special form in logistic regression.  **Question** : what so special about sigmoid function, and why we cannot use any other function with domain $[0,1]$?   2. The cost function consists of two parameters ${\rm Cost}(h_{\theta}(x),y)=-\log(h_{\theta}(x))$ if $y=1, {\rm Cost}(h_{\theta}(x),y)=-\log(1-h_{\theta}(x))$ if $y=0$. In the same was as above, I do understand why it's correct, however why is it the only form? For example, why couldn't $|h_{\theta(x)}-y|$ be a good choice for the cost function? **Question** : what is so special about the above form of cost function; why cannot we use another form? I would appreciate if you could share your understanding of logistic regression.


title: How to use sum of square reduction test, SSRT or other methods for comparing two nonlinear regression curves in R

text: as the title How to use sum of square reduction test, SSRT or other methods for comparing two nonlinear regression curves (for example, 4-parameters logistic regression curves, S curves from drc package) in R? I am a biology student and beginner of R, and I wanna get something like p-value to describe two nonlinear regression curves. If you could tell me some R package can solve my problem, I still feel very happy. Thanks


title: How to deal with very s-shaped qq plot in a linear mixed model (lme4 1.0.4)?

text: I have a dataset in which I have measured a variable (`diff`) on 38 participants (`sub`) in two different conditions (`cond.lag`). Now I am interested in whether or not the slope of the `position` of the item differs between the two conditions. That is I am interested in whether or not there is an interaction between `cond.lag` and `position` (which I center on 0 prior).               require(lme4)     options(contrasts=c('contr.sum', 'contr.poly'))          # read data     dat <- read.table("http://pastebin.com/raw.php?i=MmNQigRv", colClasses = c(NA, rep("factor", 2), rep("numeric", 2)))     # center position     dat$pos.centered <- scale(dat$position, scale = FALSE)          # fit the model     m1 <- lmer(diff ~ cond.lag * pos.centered + (cond.lag * pos.centered|sub), dat)     print(summary(m1), corr = FALSE)     ## [...]     ## Fixed effects:     ##                          Estimate Std. Error t value     ## (Intercept)             0.0819639  0.0121378   6.753     ## cond.lag1              -0.0122033  0.0155427  -0.785     ## pos.centered            0.0031775  0.0006429   4.942     ## cond.lag1:pos.centered -0.0011495  0.0011351  -1.013      Although it appears as if I don't find an interaction (only a main effect of the slope of position), I am unsure of wether or not this model makes sense as I find a very unusual (I expect heavy tailed) qq plot:               qqnorm(resid(m1))     qqline(resid(m1))      ![enter image description here](http://i.stack.imgur.com/QEWoj.png) **Questions** :   * Does this model make any sense albeit the seemingly bad behaved residuals?   * What can I do to remove those heavy tails or obtain a better model?


title: Checking if a nominal variable is important in a GLM model

text: Is it possible to put a nominal explanatory variable as a continuous variable in a GLM in R (instead of telling R that its a factor with different levels) to see if this variable (treatment) is a significant driver of the response variable that we are interested in before looking at the levels? Thanks!


title: Statistical significance for correlated significance tests

text: I refer to this article http://www.scientificamerican.com/article.cfm?id=the- mind-reading-salmon. Is there a p-value equivalent for a set of significance tests, especially in the case where some of the tests are correlated? By correlated I mean that they are not completely independent, and considering them completely independent would underestimate the significance of a good test.


title: Gibbs sampling for a simple linear model -- need help with the likelihood function

text: So in order to better acquaint myself with Gibbs sampling, I've been working on a fairly simple linear model, written in Python/R. Basically, I have 2-dimensional input data (the _x i_) and a scalar output vector (the _y i_). I am looking to fit a beta vector, i.e. _β T * xi = yi \+ εi_ ( _ε i_ is the noise). So I decided I'd use a Gaussian prior for the betas, plus an inverse gamma prior for their variance, giving me a posterior distribution that is Gaussian (thanks conjugacy!). And so I used the formula from the Wikipedia page on conjugate distributions to get the likelihood function I need to generate random samples of my beta coefficients and their variances, given the data (by their notation, _x i_'s): So, it would seem in order to do Gibbs sampling I simply iterate through _β 1_, _γ 1_, _β 2_, and _γ 2_, generating samples from the distributions listed on that wikipedia page. My confusion is, what exactly are my data points (the "samples" I'm using to calculate the distribution parameters, which I need to sample my own betas, i.e. the _x i_'s on the Wikipedia page)??? The way I see it, it could be one of two things. For discussion purposes, let's discuss _β 1_ for now.   1. First thought is, I'll just iterate over each _i_ , generating my data points using _(y i \- β2*xi,2) / xi,1_ — basically, I subtract out the influence of the second factor and divide the "leftovers" by the first factor, to obtain the impact of the first factor on the response variable. Then I can just use that population to get my distribution parameters, and finally can sample my _β 1_ and _γ 1_.   2. Second thought is, I'll just go through each _i_ , and generate my data points as ( _y i / xi_), meaning I do NOT subtract out the influence of the other factors on the response variable; it treats everything independently. So far I've been using #1, but even with artificial data sets, I'm finding that it is implying enormous beta variances that don't seem to make any sense — BUT, they are darn close to the sample variances in the populations I'm generating (my "data points" I use to calc those distribution parameters). What am I doing wrong? Is #1 or #2 the right way to do a linear model? Or am I missing something entirely? Thanks in advance for your help! Please let me know if any info would help.


title: DLM results looking wonky

text: I am teaching myself DLM's using R's `dlm` package and have two strange results. I am modeling a time series using three combined elements: a trend (`dlmModPoly`), seasonality (`dlmModTrig`), and moving seasonality (`dlmModReg`). The first strange result is with the `$f` (one-step-ahead foreacast) result. Most of this forecast appears to be one month behind the actual data, which I believe I've seen in many examples of one-step-ahead forecasting online and in books. The strange thing is that the moving seasonality is NOT similarly lagged, but hits exactly where it should. Is this normal? If I use the result's `$m` to manually assemble the componenet, everything lines up perfectly, so it's weird, though it makes sense in a way: the moving seasonality has exogenous data to help it while the rest of the forecast does not. (Still, it'd be nice to simply `lag` the resulting `$f` and see a nice match.) More troubling is the difference I see if I change the degree of `dlmModPoly`'s polynomial (from 1 to 2) in an attempt to get a smoother level. This introduces a huge spike in all three components at month 9. The spikes all basically cancel out in the composite, but obviously make each piece, say the level or the seasonality, look rather ridiculous there. Is this just one of those things that happens and I should be prepared to throw away the result's first year of data as "break-in"? Or is it an indication that something is wrong? (Even in the degree 1 polynomial case, the first year's moving seasonality's level is a bit unsettled, but no huge spike as when I use a degree 2 polynomial.) Here is my R code:               lvl0 <- log (my.data[1])     slp0 <- mean (diff (log (my.data)))          buildPTR2 <- function (x)        {        pm <- dlmModPoly (order=1, dV=exp (x[1]), dW=exp (x[2]), m0=lvl0)        tm <- dlmModTrig (s=12, dV=exp (x[1]), q=2, dW=exp (x[3:4]))        rm <- dlmModReg (moving.season, dV=exp (x[1]))             ptrm <- pm + tm + rm        return (ptrm)        }          mlptr2 <- dlmMLE (log (my.data), rep (1, 6), buildPTR2)     dptr2 <- buildPTR2 (mlptr2$par)     dptrf2 <- dlmFilter (log (my.data), dptr2)          tsdiag (dptrf2)          buildPTR3 <- function (x)        {        pm <- dlmModPoly (order=2, dV=exp (x[1]), dW=c(0, exp (x[2])), m0=c(lvl0, slp0))        tm <- dlmModTrig (s=12, dV=exp (x[1]), q=2, dW=exp (x[3:4]))        rm <- dlmModReg (moving.season, dV=exp (x[1]))             ptrm <- pm + tm + rm        return (ptrm)        }          mlptr3 <- dlmMLE (log (my.data), rep (1, 8), buildPTR3)     dptr3 <- buildPTR3 (mlptr3$par)     dptrf3 <- dlmFilter (log (my.data), dptr3)       Per the follow-on question: the data itself is monthly data for 10 years, with each month being the weekly average attendance at a theatrical production. The data definitely has seasonal and moving seasonal effects. I want to model the trend and the seasonal effects to give the management some insight, and to prepare for forecasting. (Which is not directly possible with `dlm` when you include a `dlmModReg` component, though that's the next step.) (I am trying to use an order=2 polynomial component that I believe creates an IRW trend, which is supposed to be nicely smooth.) If it matters, my moving seasonality is a yearly Big Bash Gala event that can fall in two different months, and I indicate it with 0 for most months and 1 for months in which the Big Bash falls.


title: Optimize Weights used in a Weighted Average

text: I want to calculate a weighted average, mostly for illustration purposes. However, I have an outcome variable that is Y/N and would like to optimize the weights relative to this outcome. What are some basic ways of doing this? Basically, I want/hope to use past data to drive the weights instead of arbitrary assignment.


title: Goodness-of-fit of model derived from multivariate logistic regression

text: I have a question regarding the relationship between the value of Chi square and df in determining the goodness-of-fit in a model derived from multivariate logistic regression. If the N= 290, Chi square =26.57, p=0.003 in a model with 16 variables (df is not reported), can one tell if this is legit or would this model represent overfitting or multicollinearity? Hope that makes sense. I really appreciate any help in making sense of these results.


title: Individual-level parameter estimates from ordered choice regressions

text: I have got a question regarding ordered choice regressions in **R**. I have several demographic variables with which I want to explain the ordered choice of individuals within a survey in an **ordered choice** ( **probit** or **logit** , this is not important) framework. Standard ordered choice estimations of course just give me aggregate parameter estimates. For my task it would however be useful to estimate or extract "hypothetical" individual- level parameter estimates (betas) for a certain independent variable and each individual in the survey. I have experimented with hierarchical Bayes algorithms provided by the **bayesm** and **ChoiceModelR**. Correct me if I am wrong but I think these techniques also demand that individuals appear several times within a survey and are confronted with different choice situations, so that one can estimate the influence of certain attributes on the individuals choices. My data however don't have any panel structure. I was also experimenting with Bayesian inference in example by the **MCMCoprobit** function in the **MCMCpack** package, but this function just simulates betas. I can't however, as far as I know, attribute them to certain individuals in the survey, which would be good. I would be very glad if somebody could give me a hint, sometimes already a catchword is helpful to google the correct solution!


title: Statistical test for finding the effect magnitude of different causes

text: I am fairly new to statistics and would wish to find out if there is a statistical test that tells me the magnitude/contribution/importance of different causes over a perceived effect. For example, let's say that we have different flaws, and these different flaws all contribute to an error count. The flaws are the causes, and the error count is the effect of those causes. Each flaw contributes more or less to the total error count. Is there any statistic tool that computes, with a certain degree of confidence, what would be the magnitude of each flaw in the final error count? Which flaw contributes the most to the error? The data looks like this (header and one example row):               flaw1 flaw2 flaw3 error_count      30    20    60        8      Thanks


title: Optimum value of minlev in function combine.levels

text: I have a dataset of 306967 rows and 23 columns, and I am building a regression model for predicting a factor based upon 6 other iid variables. Doing this I encountered an issue of large no of levels present in one of the iid variables, so to sort out it I used `combine.levels` in `Hmisc` package. The syntax of the general formula is `combine.levels(var,minlev=0.05)`   _Note 0.05 is the default value_.   I used 0.001 as it reduces my levels from 19881 to just 121 and my model ran well. Is it advisable to use 0.001 as value for it? Is there any technique to judge it and select its value?


title: Forecast using multivariate multiple linear regression using R for multiple dependent variables

text: I'm trying to do predict using multiple linear regression in R. I have been able to do the multiple regression bit, by converting raw data to data table. However, when i'm trying to use predict function, i'm unable to do it at multiple variable level, for each SKU ... as given below.               SKU Sales   Week    AP  Discount     101451  1308    2010_19 3.48    0.1     101451  1678    2010_20 3.5 0     101451  2003    2010_21 3.5 0.12     495172  168 2012_18 5.36    0     495172  138 2012_19 5.15    0     495172  97  2012_20 4.99    0     495172  377 2012_21 3.5 0.26     586606  130 2012_7  7.24    0     586606  143 2012_8  6.46    0.07     586606  1098    2012_9  4.99    0.28      The code I have used is (mlr is the rawdata)               library(data.table)     mlrt<- data.table(SKU=mlr$SKU,Sales=mlr$Sales,Week=mlr$Week,AP=mlr$AP,Discount=mlr$Discount)          # Run MLR model, by SKU,Display multiple Intercepts and slopes, SKU wise     mlrt[,list(intercept=coef(lm(Sales~AP+Discount))   [1],SLOPE1=coef(lm(Sales~AP+Discount))[2],SLOPE2=coef(lm(Sales~AP+Discount))[3]),by=SKU]      I'm unable to figure out how to predict using `predict` function here, for each `SKU` separately. (`Sales` is DV and `AP` and `Discount` are IDV)


title: Confused about region of rejection vs P-value.

text: What exactly is the difference between the Region of Rejection and P-value, when doing Hypothesis Testing? Read this online: > To find a rejection region, work backwards from the level of significance > (α) to a value of the test statistic (call it $z^*$). Values of the test > statistic that are farther from the hypothesized value of $\mu_x$ (i.e., > values of $z$ where $z^*$ is between $\mu_x$ and $z$) will lead to rejection > of the null hypothesis. That sounds alot like what a p-value is? Maybe I am confused about what exactly a p-value is.


title: Transforming data for canonical correlation analysis

text: Okay, I'm a stats newbie so I'll try to be as specific and clear as possible. I have a set of predictor variables (2 predictor variables) and a set of response variables (7 response variables). I am thus looking to conduct a canonical correlation analysis to investigate the amount of variance in the response variables that is explained by the predictor variables (please let me know if this is not described accurately!). First off, I am assessing if any of my variables should be transformed/standardised in some way. Three of my response variables are percentages, the rest of the variables (response and predictors) are continuous. I am guessing I need to transform the percentage data. However, I am aware that arcsine transformations and such like are now not considered useful. If so, what is the best way to deal with these? Also, my data appears severely non-normal (assessed using qq plots and shapiro-wilk tests). However, since I have a sample of more than 300 independent observations, can I assume normality via the Central Limit Theorem? Please let me know if there's anything I can provide to make things clearer in answering my question.


title: How to interpret 2-way and 3-way interaction in lmer?

text: I have a problem with interpreting 2-way and 3-way interactions in lmer. My DV is height which is a continuous variable. All IVs are categorical variables. The first factor is animal, either rat or lion. The second factor is sex, either male or female. The third factor is color: red, white, or yellow. I get confused with interpreting the output:               Fixed effects:                                       Estimate Std. Error t value     (Intercept)                       164.6888     7.8180  21.065     rat                               -14.1342     8.2889  -1.705     sexmale                           -16.0883    10.0071  -1.608     colorred                            0.5776     6.2473   0.092     coloryellow                        -14.4048     6.1025  -2.360     rat:sexmale                         15.3645    11.8567   1.296     rat:colorred                        12.5258     4.4028   2.845     rat:coloryellow                     10.3136     4.3196   2.388     sexmale:colorred                     2.0272     5.2773   0.384     sexmale:coloryellow                  5.7643     5.1669   1.116     rat:sexmale:colorred                -5.5144     6.2838  -0.878     rat:sexmale:coloryellow              0.9735     6.1690   2.158      According to Vasishth _et al._ (2007), the significance of fixed effects can be judged from the absolute t value; if it is higher than 2, then that factor is significant. In interpreting this output, I choose only factors which are significant. Please check if my interpretations are correct:   1. `coloryellow` = The height of subjects are lower when they like yellow, and are higher if they like white.    2. `rat:colorred` = The effect of rat preference enhances the preference of red, and these two promote height of subjects.    3. `rat:sexmale:coloryellow` = The effect of rat preference, being male, enhances the preference of yellow, and subjects who like rat and yellow and are male have higher height.  From these interpretations, I would like to ask: if I would like to know the effect of `lion:sexfemale:colorred`, and `rat:sexmale:colorred` compared to `rat:sexfemale:coloorred`, do I have to run new statistics?


title: References for Bayesian change point magnitude estimation

text: This paper by Adams and MacKay provides a computationally efficient algorithm for the online estimation of change points in a data sequence. I am looking for similar material that estimates not only the location of changes (or, in the case of Adams and MacKay, the size of run lengths) but also the magnitude of these changes. For example, in a Bernoulli process, suppose I want to update my beliefs about the parameter $p$, the length $r_t$ of a run at time $t$, _and_ the distribution of possible change magnitudes. I've come across some related papers but they either focus on a specific stochastic process or they rely on retrospective segmentation of the data sequence (they are offline). I am looking for online estimation. I'm guessing it wouldn't be too hard to extend the Adams and MacKay algorithm to do this, but I'd like to make sure somebody else hasn't done it first. Thanks.


title: How does missing data affect analysis of clinical trials?

text: _I was helping a friend with some problems for a clinical trials administrations course they are taking, most of it was simple statistics the last question was a little more tricky. We got the answer from google but it made me think a bit the quesiton was to do with patients dropping out of a trial with time without achieveing the requisite follow-up period. I was wondering about a simper case._ Consider a simple trial where I compare two different interventions. Patients are assigned to one of two interventions randomly, 100 into each arm. My end point is a single continuous variable. If, for whatever reason, only 94 patients in one arm and 97 in the other are suitable for analysis does this affect my analysis beyond the obvious increase in random uncertainty?


title: How to describe multi-channel sequence objects in TraMineR?

text: I have created a multi-channel object here: https://github.com/aronlindberg/VOSS-Sequencing- Toolkit/blob/master/creating_multichannel_object.R However, when I try to use various functions to describe the multi-channel sequences, I get errors saying that what I'm trying to process is not a sequence object. I created my sequence object with seqdistmc(). Is it OK to then convert this object to a common sequence object using seqdef()? This would allow me to use other TraMineR functions to describe my sequences. Best, Aron P.S. Can someone with over 300 in reputation create the tag "traminer" please?


title: What is a good way to measure the "linearity" of a dataset?

text: I have an empirically gathered dataset which relates two variables. Over a small range the relationship appears linear, however over a larger range there is clearly some second order polynomial relationship as can be seen in the image at http://imgur.com/W7f9p. I'm trying to get a measure of linearity for different ranges considered. E.g. at 20 < x < 60 or 100 < x < 120 it is very linear, but at 20 < x < 180 it is not very linear. I have tried to fit a straight line to the data and calculate the R^2 data (goodness of fit) but this shows that the straight line over the larger range has a better fit than over the smaller range. While this may be true with MS Excel, from the image it is clear that the larger range is less linear...if you hold the side of a piece of paper against the points. Is there a better way to measure the "linearity" of a dataset?


title: What is the distribution of the variance of a sample from an unknown distribution?

text: I am sampling from a parameter with unknown distribution. I would like to calculate a 95% CI for the standard deviation of the sample. @cardinal provides a nice general solution for calculating a CI in his [answer] to my previous question, Calculating required sample size, precision of variance estimate?. And @erik-p provides an estimate of the standard deviation of the variance of the sample. However, in order to calculate the 95%CI for the sample variance, it seems that I need to know the distribution of the sample variance. Is it possible to calculate such an estimate without knowing the distribution from which the sample was taken? A related question is Reference for $Var[s^2]=\sigma^4 \left(\frac{2}{n-1} + \frac{\kappa}{n}\right)$?


title: Reference for $\mathrm{Var}[s^2]=\sigma^4 \left(\frac{2}{n-1} + \frac{\kappa}{n}\right)$?

text: In his answer to my previous question, @Erik P. gives the expression $$ \mathrm{Var}[s^2]=\sigma^4 \left(\frac{2}{n-1} + \frac{\kappa}{n}\right) \>, $$ where $\kappa$ is the _excess kurtosis_ of the distribution. A reference to the Wikipedia entry on the distribution of the sample variance is given, but the wikipedia page says "citation needed". My primary question is, is there a reference for this formula? Is it 'trivial' to derive, and if so, can it be found in a textbook? (@Erik P. couldn't find it in Mathematical statistics and data analysis nor I in Statistical Inference by Casella and Berger. Even though the topic is covered. It would be nice to have a textbook reference, but even more useful to have a (the) primary reference. (A related question is: What is the distribution of the variance of a sample from an unknown distribution?) **Update** : @cardinal pointed out another equation on math.SE: $$ \mathrm{Var}(S^2)={\mu_4\over n}-{\sigma^4\,(n-3)\over n\,(n-1)} $$ where $\mu_4$ is the fourth central moment. Is there some way that to rearranged the equations and resolve the two, or is the equation in the title wrong?


title: Statistic to find out association between ordinal data

text: I have a 4x3 contigency table with ordinal data for a research finding in tuberculosis, I analysed cytology patterns and afb patterns; cytological pattern $A,B,C,D$ (independent variable) and afb pattern 1,2,3 (dependent variable) and put a 4x3 table with number of cases in corresponding cell. Which statistical test can I use to find out if pattern A is associated with afb pattern 1 and pattern C/D is associated with afb pattern 3? I cannot use chi square as one cell has no cases and two cells have < 5 cases


title: boxplot : computed values combined with raw data

text: Given computed values for a box plot (min, lq, median, uq, max) and original raw data size is it possible to merge that with raw data. given:               box1 = [min:1,lq:2,median:5,uq:8,max:10]     size1 = 50 #size of dataset that gave me the box1 calculated values     data2 = [1,1,6,...]      How can I combine both into a new box plot. Is this possible, or do I need the original raw data and have to merge it with the new raw data before calculating the box plot values.


title: nonparametric confidence interval about the difference of means for unpaired data

text: Consider two samples ${(x_i)}_{i=1}^m$ and ${(y_i)}_{i=1}^n$. Assume that the $x_i$ are independent replicates from a distribution with expectation $\mu_X$ and similarly the $y_i$ are independent replicates from a distribution with expectation $\mu_Y$, and also assume that the two samples are independent. Is it possible to get a boostrap confidence interval about the difference of means $\mu_X-\mu_Y$ ? Or is there another nonparametric way to get such a confidence interval ? EDIT : Oops - I've just seen this topic Confidence interval for the difference of two means using boot package in R Nevertheless I am interested in understanding why the method is correct. This is not a classical bootstrap procedure, isn't it ? Here we sample separately in each data sample. This is not like the classical bootstrap which is the case when there is only one data sample.


title: How to sample transactions to estimate customer-level effects

text: Let's say I am trying to estimate the treatment effect of rude waiters on how often customers frequent a restaurant. To simplify, let's assume that rudeness is binary, observable, and random. I consider the population of all transactions in October, and I see that some customers have eaten multiple meals and others only once. My outcome will be number of meals in November-Dec period. I would like to construct a customer-level cross section out of my panel data since the estimation method I am using gets complicated with non-independent observations. I have three options available to me:   1. Grab a random transaction for each customer. Use the waiter from that transaction to determine if a diner was treated or not.   2. Aggregate the October data by calculating the # of rude waiters/# of meals for each customer   3. For each person, grab a random rude transactions if there was a rude one at any time in October. Otherwise, grab a random transaction.  Unfortunately (2) is off the table since there are other attributes of the meal that are hard to aggregate. Method (1) has the problem that 2 people who both dined twice and had rude transactions half of the time might wind up in separate groups, which seems strange since they have the exact same history. Method (3) seems the most sensible, but it is giving me an estimate that have the wrong sign, and seem fragile. I have three related questions:   * Is there anything wrong with (3) from a theoretical sampling perspective?    * Can I rationalize using (1), since on average the mis-assignment bias might cancel out?   * Can I get my cross section in some other way, perhaps by adding weights of some sort?


title: Backward stepwise regression with cross validation in R

text: I would like to do model selection using backward stepwise procedure and cross validation. https://www.otexts.org/fpp/5/3 I have used `stepAIC` in `MASS` package for predictors selection and I would like to see whether cross validation is a better criteria. In other words, I am trying to use cross-validation as a criterion in backward selection procedure _instead_ of AIC. The procedure should start with a full model, drop one predictor at a time, keep the variable if it improves CV and continue until no improvements can be made. I think that the `rfe` function in the `caret` package may be able to do the job. I would also like to do bootstrapping if possible. If so, how should I specify it? Or are there any functions that can do that? I have tried writing my own function, but it seems to be beyond my ability. Many thanks.


title: test for correlation of a multivariate non-normal distribution

text: Assume I want to generate a sample of n-dimensional vector $z$ with a given correlation matrix $Corr(z)=R$. Some of the margins of its distribution are not normal (if you're interested here's the motivation behind my question), i.e. $z$ doesn't follow the multivariate normal distribution.   Let's assume that I have some imperfect method and I want to check if the sample correlation matrix is close to the theoretical one, i.e. I'd like to test the null hypothesis $H_0:R_{sample}=R$ vs $H_1:R_{sample}\neq R$   How can I do it?


title: Finding the distribution of a statistic

text: Studying for a test. Couldn't answer this one. Let $X_{1,i},X_{2,i},X_{3,i}, i=1,\ldots,n$ be iid $\mathcal{N}(0,1)$ random variables. Define $W_i = (X_{1,i} + X_{2,i}X_{3,i})/\sqrt{1 + X_{3,i}^2}, i = 1, \ldots, n$, and $\overline{W}_n = n^{-1}\sum_{i=1}^nW_i$, $S_n^2 = (n-1)^{-1}\sum_{i=1}^n(W_i - \overline{W}_n)^2, n \ge 2.$ What is the distribution of $\overline{W}_n$, $S_n^2$? How do I get an idea of the best method to use when starting a problem like this?


title: residualize binary outcome variable

text: Does it make sense and what is the correct approach to residualize a binary variable? For a continuous variable `y`, I simply run a regression that predicts `y` as a function of `x` and subtract `y_hat=xb` (the predicted valable of y from the regression) from `y`. Can I use the same procedure for a binary variable `y` so that `y_hat=invlogit(xb)` and the residualized variable ranges from -1 to 1?


title: EM algorithm for fitting GMM in multivariate case in R

text: I am currently playing around with Gaussian Mixture Models in order to model stock returns. Part of all this is using the EM algorithm to obtain MLE of parameters. I have found a package in R (mixtools) that provides the functions normalmixEM and mvnormalmixEM. I have tried it but I don't understand the output for mvnormalmix (my input consisted of a 200x2 matrix): Why do I get for two components two 2x1 mu-vectors and why do I get two 2x2 covariance matrices? Isn't that one too many? If you could clear that up or give me a link where the output of the multivariate case is explained (because I find it in general a bit confusing) I would very much appreciate it.


title: Analyzing trends in distributions

text: I have around 9000 cases; per case and per year, I have a value ranging from 0 to 1 explaining the distribution of events (for exact context, see my other question) Example data would look like this > $$ \begin{array}{c|c|c|c} \text{Case}&2002 & 2003 & 2004 & Etc\\\\\hline > \text{#1}& 0.1411 & 0.5421 & 0.6321\\\\\hline \text{#2}& 0.1121 & 0.1522 & > 0.1135\\\\\hline \text{#3}& 0.9122 & 1 & 0.1231\\\\\hline \text{#4}& 0.3331 > & 0.1353 & 0.1243\\\ \end{array} $$ And I'd like to analyze the trends for the change of the distribution of these values through out the years. **Ultimately, what I am hoping to find out is, whether the distribution is changing from concentrated (high cell values) to divided (low cell values) and vice versa.** But have no idea where to start with my limited statistics/math knowledge. Any push in the right direction is appreciated!


title: Assessing fit with identity line in Q-Q plot

text: I'm using a QQ plot to asses the similarity between two observed distributions. I need a number that quantifies how much the quantile-quantile plot _separates_ itself (is there a better word for this concept?) from the identity line `y=x` (IL). The `R^2` coefficient is of no use since it compares the points with the best fit line, not the IL. I was thinking of simply summing the absolute distances of each quantile to the identity line so that the larger this number the more I can say they are separated from the IL. Is there a better/recommended way to go about this? I'm just looking for a number that can be used to tell at a glance how far away my quantiles are from the IL. Example of one of my QQ plots in graph below. In this example, I should get a somewhat extreme (whether it is high or low) value of whatever statistic I end up using (the data set for the plot can be downloaded here). ![enter image description here](http://i.stack.imgur.com/KmFsF.png)


title: Comparing importance of different sets of predictors

text: I was advising a research student with a particular problem, and I was keen to get the input of others on this site. ### Context: The researcher had three types of predictor variables. Each type contained a different number of predictor variables. Each predictor was a continuous variable:   * Social: S1, S2, S3, S4 (i.e., four predictors)   * Cognitive: C1, C2 (i.e., two predictors)   * Behavioural: B1, B2, B3 (i.e., three predictors) The outcome variable was also continuous. The sample included around 60 participants. The researcher wanted to comment about which type of predictors were more important in explaining the outcome variable. This was related to broader theoretical concerns about the relative importance of these types of predictors. ### Questions   * What is a good way to assess the relative importance of one set of predictors relative to another set?   * What is a good strategy for dealing with the fact that there are different numbers of predictors in each set?   * What caveats in interpretation might you suggest? Any references to examples or discussion of techniques would also be most welcome.


title: How can I interpret this scatterplot?

text: Please help me to interpret this graph in terms of correlation type. Which type involves two vectors in directions observed below (see graph screenshot)? Thanks in advance. ![enter image description here](http://i.stack.imgur.com/DoIm7.png)


title: Predicting value $y({\bf{s_{0}}})$ with multivariate normal distribution

text: Suppose I have a multivariate normal ${\bf{Y}}|{\bf{\theta}} \sim {\bf{MVN}}(X {\bf{\beta}}, \sigma^{2}H(\phi))$ where ${\bf{Y}}$ is a set of observations ${\bf{Y}} = \\{y({\bf{s}}_{1}),y({\bf{s}}_{2}),... ,y({\bf{s}}_{n})\\}$ and $H(\phi)$ is a covariance matrix $H(\phi) = \rho({\bf{s}}_{i} - {\bf{s}}_{j};\phi)$. Now I want to predict a value $y_{0}$, so as it is usual, the predictive posterior distribution is $$p(y_{0}|{\bf{y}}, X, x_{0}) = \int p(y_{0}|{\bf{y}},{\bf{\theta}}, x_0)p({\bf{\theta}}|{\bf{y}}, X)d\theta$$ Usually, this is described as a product of $\text{Likelihood} \times \text{Posterior}$. However, in this case, the likelihood was initially expressed as a multivariate normal but I just want to predict the value $y({\bf{s}})$ at some point ${\bf{s}}_{0}$, so obviously I can't express a single point as a multivariate normal. What should I do? As an alternative, I should be able to use a conditional distribution $y_{0} | \bf{Y}$ from a joint distribution $\left( y_{0}, {\bf{Y}} \right)$ in the following way: $$\left( \begin{array}{ccc} Y_{1} \\\ Y_{2} \end{array} \right) \sim N\left(\left( \begin{array}{ccc} \mu_{1} \\\ \mu_{2} \end{array} \right), \left( \begin{array}{ccc} \Omega_{11} & \Omega_{12} \\\ \Omega_{21} & \Omega_{22} \end{array} \right)\right)$$ The problem is that this in some way breaks with the idea of using $\text{Likelihood} \times \text{Posterior}$ for the predictive posterior distribution. Thanks!


title: Binomial model of true/false test taker

text: I have a very large number of true/false questions. I want to take a sample of these questions and use them to test a subject. I would like to infer the proportion of questions out of the larger set that the subject will answer correctly, based on the number the subject answered correctly in the sample. Would it be a valid approach to model the subject's answers on a sample of questions from the full set as a Binomial experiment, assuming that the sample questions are selected randomly? If so, what assumptions am I making in using this model? It seems to me that I would be required to assume that the subject's answers to each question is in some way fixed, so that, for a particular question, she would choose the same answer regardless of when she was asked.