title: Probability Curve from logit model

text: Related my earlier question here. Having "mastered" linear regression, I'm trying to learn everything I can about logistic regression and am having issues turning largely "useless" coefficients into meaningful information. I asked in my previous question about graphing the probability curve for every permutation of a logit model. However, I was working on just plotting the main curve and was having some issues. If I was running a logit with just one predictor, I'd run the following:               mod1 = glm(factor(won) ~ as.numeric(bid), data=mydat, family=binomial(link="logit"))          all.x <- expand.grid(won=unique(won), bid=unique(bid))     y.hat.new <- predict(mod1, newdata=all.x, type="response")     plot(bid<-000:250,predict(mod1,newdata=data.frame(bid<-c(000:250)),type="response"), lwd=5, col="blue", type="l")      However, what about with multiple predictors? I tried to use the mtcars data set to mess around and couldn't get it. Any suggestion on how to plot the main probability curve for the logit model.               head(mtcars)          m1 = glm(vs ~ disp + wt, data=mtcars, family=binomial(link="logit"))     summary(m1)          all.x <- expand.grid(vs=unique(mtcars$vs), disp=unique(mtcars$disp), wt=unique(mtcars$wt))          y.hat.new <- predict(m1, newdata=all.x, type="response")     plot(disp<-000:250,predict(m1,newdata=data.frame(disp<-c(000:250), wt<-c(0,250)),type="response"), lwd=5, col="blue", type="l")      EDIT = OR is my previous question what I need to do. Since Y = B0 + B1X1 + B2X2, a one unit change in X1 is associated with a exp(B1) change in Y, regardless of the value of B2. Then would it be possible that my original question is really all that should be done. My appologies on my stupidity if that is indeed the case.


title: Confusion in parameter estimation of conditional random field

text: I have a certain confusion while taking the derivative of the log likelihood of the conditional random field. As given in this paper http://people.cs.umass.edu/~mccallum/papers/crf-tutorial.pdf I mean while calculating the gradient with respect to a parameter $$ \lambda_k $$ the derivate of the logZ term bring this term. $$ p(y,y'|x^{(i)}) $$ I am not sure how this term appeared while calculating the gradient. It is equation 1.22 of page 12 in the paper http://people.cs.umass.edu/~mccallum/papers/crf-tutorial.pdf Can anyone please provide some insights?


title: How to spot if model is zero inflated?

text: I'm running a linear regression. The response variable is a proportion, and has quite a lot of zeros. The predictor variable is normal distributed. My two variables look something like these:               set.seed(50)     y <- sample(c(rep(0, 20), seq(0, 1, by=0.01)), 100)     set.seed(50)     x <- rnorm(100)      And my model looks something like this:               summary(lm(y ~ x))          Call:     lm(formula = y ~ x)          Residuals:          Min       1Q   Median       3Q      Max      -0.42550 -0.31907 -0.02064  0.29710  0.59214           Coefficients:                  Estimate Std. Error t value Pr(>|t|)         (Intercept)  0.410990   0.033378  12.313   <2e-16 ***     x           -0.006252   0.033574  -0.186    0.853         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 0.3316 on 98 degrees of freedom     Multiple R-squared:  0.0003537, Adjusted R-squared:  -0.009847      F-statistic: 0.03468 on 1 and 98 DF,  p-value: 0.8527      Is there a statistic that will indicate if the model is zero inflated?


title: Time series clustering: Fourier transform and PCA

text: I have biological time series (9 years long) of the biomass of species which logically exhibit a seasonal pattern. I would like to cluster them into a few groups based on their typical seasonal evolution (e.g. spring vs. summer species). To do so, I was advised to use Fourier transform in order to decompose their signal into `N` harmonics (e.g. 3: annual, bi-annual and tri- annual seasonal cycles) and use the amplitudes and phases of these in a Principal Components Analysis (PCA; which would work as the harmonics are orthogonal/uncorrelated). I know there are already some similar subjects in this Forum, yet some aspects remain unclear to me. My questions are: (1) When I reconstruct the time evolution from the `N` first harmonics computed from the Discrete Fourier Transform (DFT), the explained variability of the original signal (the R² of the linear model between recomposed signal and the original data) is sometimes only 0.40 (`N=3`) or 0.60 (`N=5`). In your experience, does it mean the data are not suited for this approach, does that invalidate the approach? Is there more pre-processing I could do to fix that (e.g., smoothing the signals, …)? Some species exhibit sudden increases spaced by total absence, and I wonder if this doesn’t call for the need of higher frequency harmonics; should I expect difficulties there and how to tackle them? (2) Beside DFT which appears limited here, I considered using continuous Fourier Transform through a Fast Fourier Transform (FFT) algorithm and working on the power spectrum of each time series. I wonder if this could allow me to select `N'` so-called “harmonics” by selecting the `N'` highest peaks in the periodogram and then calculating the corresponding amplitude and phase to be used in a following PCA... Does that make sense? How to concretely use the info given by a FFT algorithm in R (such as `fft()` or `spec.pgram()`) in order to run a subsequent PCA (or any other clustering method)? [any R code snippet would be very welcome] (3) How to reconstruct the signal from selected harmonics in the continuous case (FFT)? I can easily do this in the DFT case, but I am stupidly blocked in the continuous case… Any R code snippet is of course very welcome. Any help regarding these questions would be very appreciated. Links toward concrete examples, especially with associated R code, would be very helpful too (as well as method name or keywords). Thank you. PS: in case it is useful: The time series are of equal length and pre- processed to have uniform sampling intervals; stationarity may be assumed; no long-term trend is in the way. I divided the time series in 52 equally-spaced observations per year (i.e., 468 observations over the 9 years).


title: Why and how does the inclusion of random effects in mixed models influence the fixed-effect intercept term?

text: The question is best illustrated by this example which uses a dataset (in library `faraway`) and `lme4` library (both in R). This intercept-only model is only for illustrative purposes.               library(lme4)     library(faraway)     data(epilepsy)     log(mean(epilepsy$seizures))#expected intercept in intercept only model = 2.5544     (Ep1a=glm(seizures~1,family=poisson,data=epilepsy))#intercept term =2.554 as expected.     (Ep1b=glmer(seizures~1+(1|id),family=poisson,data=epilepsy))#intercept term =2.214.          -------results----     1. simple glm     Call:  glm(formula = seizures ~ 1, family = poisson, data = epilepsy)     Coefficients:(Intercept)  2.554            Generalized linear mixed model fit by maximum likelihood ['glmerMod']     Family: poisson ( log )     Formula: seizures ~ 1 + (1 | id)      Random effects:      Groups Name        Std.Dev.      id     (Intercept) 0.7795       Number of obs: 295, groups: id, 59     Fixed Effects:     (Intercept)             2.214        My understanding is that the inclusion of the random term (id) tells the model that there is a repeated measures across subject (in this case). I can understand that this allows for the non-independence of the data: there are fewer than n=295 independent data points. But why does the fixed-effect intercept value decrease? Are there circumstances where it would increase? I note from the following website: http://www.danielezrajohnson.com/glasgow_workshop.R the following in relation to a model unrelated to the one i've specified above (suggest you search for "average speaker"): > "...this model has random effects for speaker and word. The fixed effects > reported are for a sort of average speaker and word. However, word, > especially, tends to be a very skewed variable. There will always be a few > very common words, that may favor or disfavor the response. The mixed model > largely counteracts this weighting." In my real example, all the coefficients (from a two categorical factor, crossed design with random-intercept-and-slope random effect - such that there are DistClaz dependencies within ID but Treat is independent) are considerably less than the corresponding mean values for those factor combinations as apparent in the raw data.               Results-----------------------------------     Generalized linear mixed model fit by maximum likelihood ['glmerMod']     Family: poisson ( log )     Formula: RCP ~ Treated * DistClaz + (1 + DistClaz | ID)            AIC       BIC    logLik  deviance       69246.20  69307.43 -34611.10  69222.20      Random effects:      Groups    Name        Variance Std.Dev. Corr             ID        (Intercept) 12.945   3.598                               DistClazAZE 11.418   3.379    -0.63                      DistClazRef  8.769   2.961    -0.80  0.71     Number of obs: 1215, groups: ID, 344     Fixed effects:                                Estimate Std. Error z value Pr(>|z|)         (Intercept)                  2.2528     0.4240   5.313 1.08e-07 ***     TreatedTreated              -2.1244     0.4916  -4.321 1.55e-05 ***     DistClazAZE                  1.5925     0.5506   2.892  0.00382 **      DistClazRef                  2.6302     0.3614   7.277 3.41e-13 ***     TreatedTreated:DistClazAZE   0.3613     0.6110   0.591  0.55428         TreatedTreated:DistClazRef   0.7068     0.4201   1.682  0.09248 .       ---     Correlation of Fixed Effects:                 (Intr) TrtdTr DsCAZE DstClR TT:DCA     TreatedTrtd -0.862                                 DistClazAZE -0.500  0.431                          DistClazRef -0.828  0.714  0.556                   TrtdT:DCAZE  0.450 -0.525 -0.901 -0.501            TrtdTrt:DCR  0.712 -0.830 -0.478 -0.860  0.583     ------------------------------------------END      These are the log(means) for the combinations:               Results-------------------------------------                  CE      AZE      Ref     No EMB  5.454940 5.808012 6.273650     Treated 3.626005 4.387088 5.253717     -------------------------------------------END      The fixed-effect intercept value (2.2528) relates to No EMB:CE (mean =5.4549)- but is this in relation to a 'sort-of-average' ID? Any pointers would be much appreciated. Thank you.


title: How many eligible bachelors in a city?

text: This is a very simple question posed to me by a friend of mine. I know it's a statistical analysis problem, but I suck at math. Given the total population of $x$ within a metropolitan area, what would be the formula to determine approximately how many single heterosexual men or women there are between age($\min$) and age($\max$)? Now as a secondary goal, how would you go about factoring in/out other relationship qualifications (such as education, employment status, etc.) from that subset of people? I am sure there are standard tools statisticians use to solve problems like this. Please explain the algorithm, bearing in mind that if you use anything more than basic algebra you'll probably lose me :-)


title: normalizing a dataset

text: I have a dataset that has budget numbers for various organizations. The numbers range from less than a million to hundreds of millions. The dataset also has information about various departments in the organization. For example it has the number of employees, the number of assets they maintain, etc. I want to be able to compare these organizations against one another, using a final numerical score, which is based on the variables of budget, employees etc. This is what I am thinking $$\frac{1}{\Sigma_i w_i} \Sigma_i w_ix_i$$ where $x_i$ is the Normalized Variable (budget, number of employees, etc) and $w_i$ is the weight assigned to each variable. The exact function to calculate the final score can be changed. I need help coming up with a normalization function that will not only normalize each variable to some specified range (say 0-1 or 1-10, it doesn't matter as long as all variables are normalized to that range), but it will also be resilient to changes in the data or to the addition of removal of an organization.


title: Accepted threshold to use for accepting results after empirically estimating a FDR from repeated randomization of the original data set?

text: ## Question: Is there a "rule of thumb" for choosing a False Discovery Rate threshold for accepting results when empirically estimating an FDR from repeated randomization of the original data set to abolish the "real" relationships? ## Example: Lets say that I have three linked data sets: `A` = $\{A_0 \dots A_n\}$ `B` = $\{B_0 \dots B_n\}$ `C` = $\{C_0 \dots C_n\}$ The sets are linked in that comparisons are only valid between members of: `comparison_sets_real` = $\{A_i, B_i, C_i\}$ I have a distribution (`REAL`) of results generated using `comparison_sets_real`. I then shuffle the order of sets `A`, `B`, and `C` such that `comparison_sets_null` = $\{A_r, B_r, C_r\}$ represent the same data set but with the real relationships abolished. A large number of randomized distributions are generated (~1000) and a median distribution (`NULL`) is derived from the set of randomized distributions. From comparing the `REAL` and `NULL` distributions using a reversed cumulative histogram we can generate a progressively more stringent FDR threshold [1]. **Is there a standard guideline similar to the $\alpha$ of p-values?** ### Footnotes:   1. In the provided images, the randomization replication is set at 15 to reduce calculation times. Actual results would be generated with ~1000 replications.


title: Regression equation formatting

text: I have a simple question for you, which has to do with style. Since I am a novice in writing research papers, I have the small issue of not knowing how to represent an equation in an acceptable way. I am running an OLS regression, for which I use year fixed effects. Since I do not wish to write up every b*year-dummy combination, what is the appropriate way to write this in a paper? I am thinking of writing it as: yi = β0 +β1xi+ ... + SUM(b*year-dummy) + e Is there a standard way to write this, or is it common to just say the regression included time fixed effects in writing? Your answers are much apprech :)


title: Validation accuracy larger than training accuracy

text: I was performing an experiment but got a higher validation accuracy than training accuracy. I've got a 39 mice data and performed leave one out cross- validation. The validation accuracy was 100%. But when I repeat the experiment on the whole training data (with optimal parameters selected by cv), the highest training accuracy is only 87.18%. I am pretty sure there's no problem with my code but this really confused me :( Could anyone have ideas on potential problems or explanations? Thanks!


title: Dataset for multi class perceptron

text: I am developing a multi-class perceptron algorithm and was wondering if there are any datasets that could be used to test a multi-class perceptron? - A dataset where the classes are linearly separable and have at least 100 or more instances for training?


title: How can a quadratic relationship be written two ways

text: In this answer, under: "Why include the linear term?" in the answer, it mentions the relationship can be written two ways. Later it says that the point where $x=b$ is the vertex of the parabola. > It is illuminating to notice that a quadratic relationship can be written in > two ways: > > $y=a_0+a_1x+a_2x^2=a_2(x−b)^2+c$   >  (where, equating coefficients, we find $−2a_2b=a_1$ and $a_2b^2+c=a_0$). > The value $x=b$ corresponds to a global extremum of the relationship > (geometrically, it locates the vertex of a parabola). > > If you do not include the linear term $a_1x$, the possibilities are reduced > to   >  $y=a_0+a_2x^2=a_2(x−0)^2+c$   >  (where now, obviously, $c=a_0$ and it is assumed the model contains a > constant term $a_0$). That is, you force $b=0$. > > In light of this, question #1 comes down to whether you are certain that the > global extremum must occur at $x=0$. If you are, then you may safely omit > the linear term $a_1x$. Otherwise, you must include it. My question is: Where did the betas come from, or how did you get: $a_2(x−b)^2+c$ from $y=a_0+a_1x+a_2x^2$


title: Why can't a Self Organizing Map to recognize handwriting train to acceptable error rate when number of samples is greater than output neurons?

text: I have a simple Self organizing map that I have been working on. It has an output neuron for all 26 letters, and when I provide it with 26 samples (one for each letter) it quickly is trained to an error rate below 1%. It is fairly accurate, but when I want to add more samples, one for each output neuron, the lowest possible error rate shoots up to around 37%. Even when I use a decaying learning rate it does this. Even If I add one more sample so that I have 27 instead of 26 the training gets stuck on around 16% error rate. I don't understand why as soon as I have more samples than the number of output neurons I cannot train it properly. I am using the subjective learning method, yet the same happens when I switch over to the additive method. I am still a beginner with neural networks so maybe I am not understanding something properly. Any help would be greatly appreciated!


title: What level of measurement is a factor score?

text: After doing a factor analysis, I decide to use that factor to test a hypothesis. However, I am not sure what data type of the Factor score the is (i.e., nominal , scale, etc.). Can I do a t-test or ANOVA as usual?


title: Confidence intervals around a centroid with modified Gower similarity

text: I would like to obtain 95% confidence intervals for centroids based on Gower similarity between some mulivariate samples (community data from sediment cores). I have so far used the `vegan{}` package in R to obtain modified Gower similarity between cores (based on Anderson 2006; now included in R as part of `vegdist()`). Does anyone know how I can calculate 95% confidence intervals for the centroids of, for example, sampling sites, based on modified Gower similarity? Additionally, if possible, I would like to plot these 95% CIs on a PCO that shows the centroids, so it's evident if they're overlapping. To get modified Gower similarity, I used:               dat.mgower <- vegdist(decostand(dat, "log"), "altGower")      But as far as I know, you don't get centroids from `vegdist()`. I need to get centroids, then 95% CIs, then plot them... in R. Help! Anderson, M. J., K. E. Ellingsen, and B. H. McArdle. 2006. Multivariate dispersion as a measure of beta diversity. Ecology Letters 9:683–693.


title: Calculate one median for data from five experimental repetitions

text: I have an experiment that produces a few thousand data points which I then have repeated 5 times. Now I wonder how to calculate one single median to summarize all of them.   1. I could combine the data from the five repetitions and calculate a median for that.   2. I could calculate a median for each experiment and then average them.   3. I could calculate a median for each experiment and then take the median of those.   4. Some other procedure I haven't thought of. Which procedure make the most sense as a way to combine the data? How would I then go about producing a confidence interval for the median? I assume some sort of bootstrap procedure makes sense, since the distribution of the data is not normal.


title: Correlated alpha-stable processes with different alpha values

text: Can one define a pair of stochastic processes with independent alpha-stable increments and different alpha parameter (index of stability) values (e.g. a Brownian motion and a Cauchy or Levy process) which are correlated?


title: Is centering needed when bootstrapping the sample mean?

text: When reading about how to approximate the distribution of the sample mean I came across the nonparametric bootstrap method. Apparently one can approximate the distribution of $\bar{X}_n-\mu$ by the distribution of $\bar{X}_n^*-\bar{X}_n$, where $\bar{X}_n^*$ denotes the sample mean of the bootstrap sample. My question then is: Do I need the centering? What for? Couldn't I just approximate $\mathbb{P}\left(\bar{X}_n \leq x\right)$ by $\mathbb{P}\left(\bar{X}_n^* \leq x\right)$?


title: How to put data in time series regression

text: I tried to get a result in regression on `LIBSVM`, and I get the same problem on single feature data or multivariate data. Suppose I try to find the price of something, and I have the data for $N$ days. I give the day ($1,2,\ldots,N-1$) price for $X$ as independent variables, and I train it with the day ($2,3,\ldots,N$) price as the dependent variable. So I try to get the next day's price always. When I get the solution ($Y$) on test data, surely I expect the data price for the days ($2 \ldots N$ for the test data). But I realized that the $Y$ values are always seems like the X values. When comparing the solution $Y$ with $X$ values, ($1,2,\ldots,N-1$) days, results are more successful than ($2,3,\ldots,N$) days. In other words I expect the next days values but the result is more successful on the same day again. What is wrong with the approach?


title: A problem on a.s. convergence

text: I'm preparing for an exam and I came across this problem from old exams. I'm really clueless on how to solve it. Consider a sequence of random variables $\{X_n\}_{n=1} ^\infty$ defined on the probability space $([0,1],B[0,1],\lambda)$ where $\lambda$ is Lebesgue measure. Define $X_n(\omega) = 1_{[1/2-(2n)^{-1}, 1/2+(2n)^{-1} ]}(\omega)$ for all $\omega \in [0,1]$. Prove that $X_n \xrightarrow{a.s.}0$ as $n \rightarrow \infty$.


title: How to investigate a 3-way interaction?

text: There is a significant 3-way interaction in a data-set I'm working with. The interaction involves both categorical and quantitative variables. I have been pointed towards simple slopes and this website but I find the explanations lacking. I only have a basic background in statistics, and Googling for other examples has not been very helpful. Any insight as to how to begin to understand this interaction will be very welcome, especially using R rather than SPSS. Thank you very much in advance, from a lower-year college student who has suddenly found himself over his head!


title: fitting an exponential function using least squares vs. generalized linear model vs. nonlinear least squares

text: I have a data set that represents exponential decay. I would like to fit an exponential function $y = Be^{ax}$ to this data. I've tried log transforming the response variable and then using least squares to fit a line; using a generalized linear model with a log link function and a gamma distribution around the response variable; and using nonlinear least squares. I get a different answer for my two coefficients with each method, although they are all similar. Where I have confusion is I'm not sure which method is the best to use and why. Can someone please compare and contrast these methods? Thank you.


title: Should I be using "pooled variance" in this situation?

text: Two variables were known to me: `dark sample` and `light sample` To see if measured distance of sample components affected the shade of samples, I measured component distance of three dark samples and three light samples.   * I measured each sample three times - for repeatability -    * took the mean of the measurements for each sample.    * Then I took the Std Dev of the 6 means    * and then took the pooled variance of the 6 samples. soooo.....               light#1 - 4.56mm, 4.57mm, 4.56mm = mean 4.56       light#2 - 4.47mm, 4.45mm, 4.43mm = mean 4.45         light#3 - 4.37mm, 4.36mm, 4.37mm = mean 4.37        dark#1  - 4.47mm, 4.45mm, 4.42mm = mean 4.45        dark#2  - 4.40mm, 4.41mm, 4.42mm = mean 4.41        dark#3  - 4.41mm, 4.44mm, 4.41mm = mean 4.42        * Then I took Std Dev of Light means `= 0.04` and Std Dev of Dark means `= 0.02`   * Then I took the pooled variance of the samples: $$(3-1)(0.04) + (3-1)(0.02) / (3-1) + (3-1) = 0.03$$ Is that the correct use of pooled variance?


title: Following up a three-way interaction found with (mixed) model comparison

text: I hope this isn't a silly question, I'd like some advice on following up a threeway interaction in a mixed effects model. I've been building my models incrementally, like this:               lmer0<- lmer(predval ~ 1 + (1 + Quantile|Subj), data = d)     lmer1<- update(lmer0, .~. + epoch)     lmer2<- update(lmer1, .~. + Quantile)     lmer3<- update(lmer2., .~. + Congruence)     lmer4<- update(lmer3, .~. + FixDur)     lmer5<- update(lmer4, .~. + Quantile:Congruence)     lmer6<- update(lmer5, .~. + Quantile:FixDur)     lmer7<- update(lmer6, .~. + Congruence:FixDur)     lmer8<- update(lmer7, .~. + Quantile:Congruence:FixDur)      Here's what I get when I test these models against each other. All interaction terms improve the model (indicated by decreasing AIC/BIC values). This printout indicates the fully interactive model (lmer8) is the best model so far.                       Df     AIC     BIC   logLik     Chisq Chi Df Pr(>Chisq)         lmer0   5 2062522 2062573 -1031256                                     lmer1    6 2053156 2053217 -1026572 9368.2240      1  < 2.2e-16 ***     lmer2    7 2053144 2053215 -1026565   14.2730      1  0.0001581 ***     lmer3    8 2051857 2051938 -1025921 1288.7862      1  < 2.2e-16 ***     lmer4   11 2051604 2051716 -1025791  259.3176      3  < 2.2e-16 ***     lmer5   12 2051606 2051728 -1025791    0.0284      1  0.8662848         lmer6   15 2051168 2051320 -1025569  443.8425      3  < 2.2e-16 ***     lmer7   18 2051152 2051335 -1025558   21.9526      3  6.673e-05 ***     lmer8   21 2050697 2050911 -1025327  460.9543      3  < 2.2e-16 ***      It's the last term in the model that really interests me - the three way interaction. Unfortunately I can't embed the image here (newbie), but you can find a plot of the data shows **Congruence** (shape) by **Quantile** (Colour gradient/x-axis) by **Fixation Duration** (facet_grid) here: ![enter image description here](http://i.stack.imgur.com/H6ao2.png) What I think I should do now is as follows:   * Test the Congruence x Quantile 2-way at each level of FixDur (i.e. Begin new model comparison at each level of FixDur, that's each column in the plot above).    * For those two-way interactions that are significant (they all are), follow with t-tests at each Quantile (i.e. 20 paired t-tests per FixDur). I'm not too sure how to conduct these using lmer, I'm basically just using the t.test function. The obvious issue with this will be multiple comparisons. Essentially I want to do 20 t-tests at each level of Fixation Duration (so I can say something about the timecourse of my effects). That's a lot of t-tests, so using a Bonferroni or Scheffe correction is going to kill me! Does anyone have a suggestion on what method or approach would be suitable for dealing with this problem? It's a hunch, but I feel as though the co-occurence of significant p-vals should count for something. That is, if the quantiles that yield a t-vals significant at .05 occur all in a row (i.e. show the effect emerging, as in the first column of the graph), that should give us confidence that they're not type 1 errors. Am I going down a dangerous road? This is a terribly long post i'm sorry. I'd welcome all thoughts though on any part of it! My thanks to you in advance.


title: Variance estimates in k-fold cross-validation

text: K-fold cross-validation can be used to estimate the generalization capability of a given classifier. Can I (or should I) also compute a pooled variance from all validation runs in order to obtain a better estimate of its variance? If not, why? I have found papers which do use the pooled standard deviation across cross- validation runs. I have also found papers explicitly stating there is no universal estimator for the validation variance. However, I have also found papers showing some variance estimators for the generalization error (I am still reading and trying to comprehend this one). What do people really do (or report) in practice? **EDIT:** When CV is used to measure the crude classification error (i.e. either a sample has been labeled correctly or it hasn't; e.g. true or false) then it may not make sense to talk about a pooled variance. However, I am talking about the case in which the statistic we are estimating does have a variance defined. So, for a given fold, we can end up with both a value for the statistic and a variance estimate. It does not seems right to discard this information and consider only the average statistic. And while I am aware I could build a variance estimate using bootstrap methods, (if I am not very wrong) doing so would still ignore the fold variances and take only the statistic estimates into consideration (plus requiring much more computation power).


title: Can I use a measure when it has low reliability?

text: I would like to use a 5-item subscale to measure the effects of a behavioural intervention program on the use of behavioural strategies. I obtain significant results when I analyse it in models however the reliability of this measure is very low at .5 and the mean inter-item correlation is below .2. Can I use this measure or should I exclude it from all subsequent analyses? Thank you.


title: Am I using the wrong statistical test?

text: I am using $\chi^2$ to analyze a set of data gathered over 4 time periods. The two variables are independent--each 'event' is categorized as 'good' or 'bad'. We want to prove that the intervention caused an increase in the good events and decrease in the 'bad' events, or at least that the increase in good events outpaces the bad events. Here is a sample of what the data might look like:                      good    bad      m1     1035    2278     m2     1152    2643     m3     1189    2917     m4     1125    2974      A $\chi^2$ test returns 15 which would have me believe the result is significant, but the change doesn't look very significant at all. In fact, I don't think there is a case to be made that there wasn't any improvement at all. I think that this may be due to the fact that my $N$ varies in each measurement (though I thought $\chi^2$ controlled for that). It also seems strange to me that if I simply move the decimal over two places (thereby changing my N but not the proportion of good results to bad) then I get $\chi^2=.15$ which is closer to what I would expect. I know I need to change my approach drastically. How can I make this comparison? Should I be using two independent student's t-tests? Or maybe just use the correlation coefficient for each column and the coefficient of determination? Thanks in advance for your help and please let me know if anything is unclear. EDIT: Here is a thought: should I set up benchmarks for expected improvement, something like 3-5% increase (proportional to $N$) in each measure? This way my expected values in my $\chi^2$ might more accurately reflect what I am trying to show. EDIT for clarification: the % of good events compared to bad over the measures:   31%, 30%, 28%, 27% Clearly there has not been improvement with regard to the proportion of good to bad so what test should I be using to demonstrate this?


title: Subjective Bayesian's care for real world validation and classical statistician's worry about CI related paradoxes for a given data set?

text: I was thinking about CI and subjective Bayesian and I have following two questions:   1. If a subjective ( _not objective_ ) Bayesian would care if her predictions don't do well in the real world.   2. A classical statistician would not care if her confidence statement is (obviously) wrong for a given data set (as in Welch's Paradox, where conditioning is on ancillary statistics leads to the resolution of pathological behavior). I think my answer for 1. is YES and 2. is NO. But I don't know if I am thinking along the right lines. Can I have some more insights? * * * **UPDATE** **Welch's example** : This example works for any $n$, but we will take $n=2$ for simplicity. $X_1, X_2 \sim U(\theta - 1/2, \theta +1/2)$ (iid), $\theta \in R$. This implies $X_1 - \theta \sim U(-1/2, 1/2)$ (iid). $(X_1 + X_2) /2 - \theta$ (note that this is _NOT_ a statistic) has a distribution independent of $\theta$. We can choose c > 0 s.t. $Prob_{\theta} [-c \le (X_1 + X_2) /2 - \theta \le c] = 1- \alpha (~0.99)$, implying $((X_1 + X_2) /2 - c, (X_1 + X_2) /2 + c)$ is the 99% CI of $\theta$. The interpretation of this CI is: if we sample repeatedly, we will get different $(X_1 + X_2) /2$ and (at least) 99% times it will contain true \theta. **But for a particular set of $X_1, X_2$, we can't say if the CI contains $\theta$**. Now, consider the following data: $X_1 =0$ and $X_2=1$, as $|X_1 - X_2|=1$, we know FOR SURE (Prob =1) that the interval $(X_1, X_2)$ contains theta (one possible criticism, $P{|X_1 - X_2|=1} = 0$, but we can handle it mathematically and I won't discuss it). (Better details are in Pratt, 1961; Lehman, Chap 10, 2nd Edition, Prob 27, 28; Kiefer, 1977; Berger and Wolpert, 1988) Thanks, S.


title: Do negative probabilities/probability amplitudes have applications outside quantum mechanics?

text: Quantum Mechanics has generalized probability theory to negative/imaginary numbers, mostly to explain interference patterns, wave/particle duality and generally weird things like that. It can be seen more abstractly, however, as a noncommutative generalisation of Bayesian probability (quote from Terrence Tao). I'm curious about these things, though by no means an expert. Does this have any applications outside Quantum Mechanics? Just curious.


title: Analysing choices pattern

text: we have a process in which, at each step, a set of elements are presented to user, the user choses one, his choice is recorded and next round starts with a new set of elements.   For example:   1. {20,50,80} and the user chose 80   2. {80,110,140} and the user chose 80 again   3. ... We want to verify the hypothesis whether the user prefers large elements. You can assume we have many steps. How does the following suggestion fare ?   * at each step compute MEAN, and STD of the presented values.   * normalize elements and the choice to be in standard units, i.e.               standard_value = (value-MEAN) / STD      for example: in (1) above the used chose -1 and in (2) he chose +1   * count how many times each standard unit is chosen and draw a scatter plot. So for example we'll plot points like: (-3.5,0), (-1.3, 2), (0,4), (1,7), (2,10), Is this too simplistic ?


title: How to get percentiles from empirical density in R?

text: The `density()` function in `R` allows me to enter observations and get an empirical density that I can plot x and y values. I like it because it allows me to weight observations according to how important they are, and it allows me to specify the smoothing bandwidth I want. My question is once I run the `density()` function how do I obtain percentiles from this density? Note this isnt the same as just getting the sample percentiles from my data, because I'd like to use weights on the observations.


title: test validation versus k-fold cross validation

text: I am attempting to use a neural network, after using other machine learning algorithms. I am using the RSNNS package (I am willing to use / evaluate other packages) that's part of R. I would like to get a precision that's at least 66%. I split the data in a training and test sets, with 4/5 of the data in the training. I then trained models using different network layouts and learning rates, using the same training set each time. I selected the parameters that gave precision >66% and the largest F-measure on the test set. The parameters I selected gave a precision of 70% on the test set. I then took the data and did a 10-fold cross validation using the same network layout and learning rate. With this k-fold cross validation, I get a precision that is just above 50% (which is similar to the other learning algorithms I used). My question is, is the 70% precision accurate with the test set? Is my k-fold validation possibly finding local optima, and not giving an accurate precision? **EDIT** Since this seems like an important point I left out, there are 2 classes, positive and negative. It's 18% positive, and 81% negative. There are about 550 cases. Following Matteo's suggestion in the comments, I ran the network again multiple times. I just used the best selected parameters, because the neural network takes some time to run. I split into training (80%) and test (20%) sets again, except did 10 random splits using `sample`. Since it's random, some of the data appears more often in the training sets than the test sets. Using this, the precision ranged anywhere from 30% to 70%. When I averaged the 10 runs together, it came out to just above 50% precision. I am leaning towards saying that is the best precision I can get using this data set, since the earlier machine learning algorithms gave a similar precision (data not shown).


title: Confidence Interval and Signifcance test Question

text: I have done the first part of the question already I just need help on the second part. You read in a U.S Census Bureau that a 99% confidence interval for the mean income in 2005 of American households headed by a college educated person at least 25 years old was 100,272±1651. Based on this interval, can you reject the null hypothesis that the mean income in this group is $95,000? What is the alternate hypothesis of the test? What is its significance level? I answered that you can reject the null hypothesis at the 0.01 level as confidence interval is between 98621 and 101923, and 95,000 is not in that in that range. I am confused about the alternate hypothesis is it Ha>95,000 or is it Ha does not equal 95,000. Also How do I get the signicance level of the null hypothesis? Thank you very much everyone.


title: Combining two Odds ratios with the same control group

text: I am doing a meta-analysis of observational studies. As the disease I am working on has different stages, some of the studies reported ORs (95%CI) for each stage separately in comparison to a single (shared) control population. I wonder if we can combine these two ORs (and their 95%CI) given that they reflect two case groups but a single control group. Some other info:   * The sample size for each case group is different.   * We have the sample size for each of the groups but not the frequency of exposure (otherwise we could use just the raw data).   * **Please note that I have no problem combining the ORs, but simply combining two ORs without taking into account that they both have the same control group will cause narrowing of the 95%CI and overweighting of the study in the analysis**. Any ideas?


title: How can i show mathematically Partial Least Square Regression is better than other Ordinary Least Square Regression?

text: I want to develop techniques for attribute selection (important independent variable X) using Partial least square 2 regression(PLS2R) for a large data sets .Initially i tried using multivariate linear regression using many independent variable (X) and many dependent variable (Y).So i got different regression model for Y1,Y2,Y3 and so on . Here Y1,Y2 and Y3 are separately modeled and attributes are selected separately for each Y1. Instead of modelling separately , i want to store my responses in single matrix Y and model with independent variables matrix X.so i will get a Partial least square regression model or a two-block model Y~X. Does use of block model using PLS2R(simultaneously using X and Y matrix ) provides more information for attribute selection in comparison to other classical regression techniques ? And Can it be shown mathematically that PLS2R method is better than other multivariate regression method? If you have any good references or suggestions please let me know ,i would appreciate it . Thanks in Advance


title: Does one group respond differently from others

text: I'm looking for the best way to test the hypothesis that two groups are different considering the relationship between two variables. There is a population of individuals and, for each individual $x$, I know a value $m(x)$ from 0 to 100 and a binary value $s(x)$. I know there is a strong correlation between $m$ and $s$. (Higher $m$ means $s$ is more likely to be 1.) Each individual also is a member of one of a small finite set of groups, $g(x)$. By binning on $m$ and plotting bin number versus the percentage of $s$s that are 1, separately for each group, I can see that there appear to be differences: E.g. those in group $A$ with low $m$ scores are more likely to have $s=1$ than those with low $m$ scores who are not in group $A$. I can also see that members of group $B$ seem to have a lower correlation between $m$ and $s$ than members of other groups. But is this observed difference significant? So, for each group $A$ I'd like to test the hypothesis that being in group $A$ means that the relationship between $m$ and $s$ is significantly different than for those not in $A$. The null hypothesis is that members of group $A$ are just the same as those who are not in group $A$. My question then is: what sort of analysis should I look into applying in order to test this sort of hypothesis?


title: Is hold-out validation a better approximation of "getting new data" than k-fold CV?

text: I've been rethinking an answer I gave to a question a couple weeks ago Hold-out cross-validation produces a single test set that can be used repeatedly for demonstration. We all seem to agree that this is in many ways a negative feature, since the one held-out set could turn out to be non- representative through randomness. Moreover, you could end up overfitting to the _test_ data in the same way you can overfit to the training data. However, it seems to me that the static nature of a held-out sample is a better approximation of "getting more data" than k-fold CV, and avoids the issue of averaging across folds. I can't, however, come up with any statistical basis for this feeling I have. Is there any logic in my intuition? For instance, what I have in mind for an upcoming project is first using hold- out validation to build and test a model, then as a validation step re-drawing the hold-out set several times to show that my estimates of prediction error (on the test set) are robust to sampling error in the test set. Is this a bad idea for any reason? This question was asked before but never received an answer.


title: Identify category for combination of words from categories of component words

text: I have the following matrix:                                  category                       I    II   III                   1  70   300     5         word       2  20   400    45                   3  10   300     0      It shows the number of times that a word_i appears in a phrase of category j, for example, word_1 appears 300 times in category_II. I would like to know the probability of a given phrase to be of category I, II, or III. E.g. if the phrase consist of word_1 and word_2. Two approaches I tried to calculate the probability of the phrase to be of category I:   1. 70/(70+300+5) = 0.186   2. normalization for each category -> 0.7/(0.7+0.3+0.1) = 0.64 In the first case I am biased by the large number of case that fall into category II. In the second, I ignore totally that I am (should be) able to predict case of category II since I do have a lot of statistics. **Goal** : I would like to have a classifier that avoids being fooled by the large statistics of category II, while being precise when predicting category II. My gut feeling is that I need some measure likelihood here.


title: How to calculate prediction sensitivity/specificty when [True/False/Don't know] outcomes are possible?

text: I have the following dataset:               OXXO      The output of classification method #1 is:               OOXO      Classification method #2 has a prediction confidence estimator and it's output is:               O~X~      where ~ indicates that method #2 refuses to produce a prediction due to a low confidence level. Assuming that no answer is better than a wrong answer, how should I compare the predictive power of these two methods?


title: Training a convolutional neural network

text: Based on my research on convolution neural networks, every other layer in such a network has a subsampling operation, in which the resolution of the image is reduced so as to improve generalization of the network. So, a CNN could consist of an alternation of convolution and subsampling layers. However, when using backpropagation to train a convolutional neural network, I don't quite understand how one would train a convolutional layer. When you are training a convolution layer, don't you need the weights of the next layer to calculate delta? Based on my understanding of backprop, the equation to find delta of a hidden layer is as follows: `Neuron.Delta = Neuron.Output * (1 - Neuron.Output) * ErrorFactor` But, to find ErrorFactor, you need the weights of the connections between the current layer and the next one. And if the next layer is a subsampling layer, then there will not be any weights to use to calculate delta. My current solution to this problem is to simply look at the next layer's weights to calculate the delta of the current layer. So, if layer 1 is a convolution layer, layer 2 is a subsampling layer, and layer 3 is a convolution layer, I would look at the weights connecting layers 2 and 3 to calculate the delta at layer 1. Is this a correct understanding of how to train a convolution neural network?


title: Use correlation to compare two time series?

text: I have two time series data sets which contain hourly-intervalled, monthly, and yearly household electricity consumption in kWh. One data set is produced by a simulation, the other gathered from the real-world. My aim is to validate the simulated output by using the data gathered from the real-world. I want to measure the similarity between these data sets, and be able to say if these are statistically similar. My first intuition is to use a correlation coefficient such a Pearson product moment correlation. But from what I read in previous posts that in general the correlation coefficient between two time- series may be a very poor metric. I'm not very keen on statistics related to time series, but would something like a cross-correlation or maybe ARIMA do the trick? Could someone please point me in the direction for a technique which I can use in SPSS?


title: PCA and SVD matlab

text: I am having some problem in computing SVD and PCA in Matlab. I do not know if I am doing theoretical mistakes or programming mistakes. Starting with a data matrix $X$ PCA computes the eigenvalue $\lambda_i$ of the matrix $X^TX/(n-1)$. On the other side for SVD $X=U\Sigma V^*$ and so $X^T X=V \Sigma^T U^T U\Sigma V^T=V \Sigma^T\Sigma V^T$ where $\Sigma$ is diagonal with elements $\sigma_i$. So $\lambda_i=\sigma_i^2 /(n-1)$. Now I try a matlab example:                   load hald;         [u s v]=svd(ingredients);         sigma=cov(ingredients);         [a,b]=eig(sigma);         disp('sigma')         disp(diag(s))         disp('lambda')         disp(diag(b))      > >               sigma       211.3369        77.2356        28.4597        10.2667          lambda         0.2372        12.4054        67.4964       517.7969      The obtained value do not respect the original equation. Where is the mistake? Thanks :D


title: Given a set of sub-graphs, how to infer the underlying graph?

text: There are two questions I feel puzzled in the recent period.   1. Firstly, if given a set of sub-graphs, which are sampled from the underlying graph, how can I infer the underlying graph given the set of observed sub-graphs? If you can refer some good papers to me, I will appreciate it very much. This is just one case.    2. Another case is that a set of sub-graphs are sampled from several underlying graphs: How can I infer the underlying graphs given the set of observed sub-graphs? I welcome any comments.


title: Chi-Squared Formula

text: I have several nonlinear curves of the transmission spectrum of a semiconductor. I got a fit to each of the curve. To get the chi squared for each my teacher gave me the formula $$\chi^2=\sum_{i=1}^{N}\frac{(T_{exp}-T_{fit})^2}{N\sigma^2_{exp}} $$ where $$\sigma_{exp}=\sum_{i=1}^{N}\frac{\sigma_{i}}{N}.$$ However, reading some statisics book I haven't seen this formula anywhere and have only seen $$\chi^2=\sum_{i=1}^{N}\frac{(T_{exp}-T_{fit})^2}{\sigma^2_{i}}.$$ Can someone explain the top, or is the formula just wrong?


title: Large data analysis - learning resources

text: My question is very simple: which learning resources (books, courses, online courses, and so on) about "large data analysis" would you suggest to a graduate with a strong background in Machine Learning and Computer Science? **Edit:** I am looking for something similar to section 17 of this very well known online course: https://www.coursera.org/course/ml (you can click on preview lectures to check it out) but in very much detail and many other different methods/algorithms.


title: Measures of multidimensional spread or variance

text: What's a good measure of spread in a multidimensional space? In a single dimension _variance_ would be the measure I need, but in a multidimensional space I need more than just variances. Note that in a single dimension, I need something different than the _range_ covered by the data points. Consider a hypercube representing the space of all possible data points (all attributes are limited to values in the interval $(0,1)$). I need a measure of spread which is optimal when all corners of the cube are equally populated (and there are no other data points within the cube). When only two opposed corners would be populated, all variances would be maximal, but this is not what I need. Also a population which covers all corners, but also contains points more toward the middle of the cube should have a lower spread. My first hunch would be to sum up all values on the main diagonal of the covariance matrix, and substract all other values (all proper covariances). This idea is very ad hoc, however, and I don't know whether this is thinking in the right direction. Please help me find a good measure of spread/variance in a multidimensional space.


title: Weeding out multi-collinear categorical variables

text: I have a huge amount of features that are categorical variables and I'm trying to find a system for weeding out categorical variables that are close to being multicollinear. Is vif a reasonable measure in the context of categorical variables? Would pearson's chi-squared be a good measure of association between categorical variables? I've also heard mention of polychoric correlation. Edit 1: For vif on categorical variables, I'd need to use the r-squared from a linear probability model right? I'm assuming you can't use pseudo-r-squared since it's a completely different animal.


title: Standard latent Dirichlet allocation corpora?

text: I have an implementation of latent Dirichlet allocation and I was wondering if there are any standard corpora that people have used to test LDA. Also, I would like to know how to get them.


title: Is an I-GARCH process stationary?

text: So I know a GARCH process is stationary, however I'm having doubts for the I-GARCH process as there is need for the coefficients of the process to sum up to 1, which does not satisfy the stationary condition. What's the difference between an I-GARCH and a normal GARCH?Why would you use one and not the other?


title: Test incorrect functional form when residuals have non-normal distribution

text: J. B. Ramsey (in _"Tests for specification errors in classical linear least- squares regression analysis." Journal of the Royal Statistical Society. 1969_ ) says that the RESET test assumes that the residuals are normally distributed. If one wants to test the incorrect functional form of a model but the residuals have a **non-normal** distribution, how can it be done? Ramsey also says that " _the cases where mis-specification leads to a non- normal distribution of û [residuals] are to be discussed in a later paper_ ". Does any one know which paper is this?


title: Paired t-test for binary data

text: I have one sample with n=170 and two binary variables (A,B) that can take as a value 1 or 0, where 1 counts as a success and 0 counts as a failure. What I want to know is whether the means of these two variables are equal. To find this out I generate a new variable that takes the difference between these two variables called C, so C = B-A. I then compute the p-value for the hypothesis that C is normally distributed with the Shapiro-Wilk test and I find a p-value of .96, so I choose not to reject this hypothesis. Apart from that the difference is normally distributed, I am not worried about the other assumptions required for a paired t-test. **Question:** _Can I use the paired t-test in this circumstance or is it a mistake to use the Shapiro-Wilk test for binary data to check for normality and should I use the Wilcoxon sign rank test instead?_ I would much prefer to use the t-test, because I believe it has a higher power than the Wilcoxon sign rank test, but that higher power pretty much does not matter if the test used is the wrong one. Cheers, Martin


title: Comparing category distributions

text: I have categorical count data (for categories C1 to C3, but potentially several more categories) for two datasets:                       | --- Dataset 1 --- | --- Dataset 2 --- |             |  C1    C2    C3   |  C1    C2    C3   |     Item 1  |  0     200   300  |  0     2      3   |     Item 2  |  0     200   300  |  5     0      0   |      The total number of data-points in each dataset is different (500 and 5 in this example). What statistical test should I use to determine if the distribution of counts for each item across the categories is the same between the two datasets? For example, the distribution of Item 1 is the same across the two datasets, but the distribution of Item 2 is not. I will test each item separately.


title: Adding prediction on the test corpus to the training corpus

text: I have the following train and test sets:               train: 1000 obs. of  40 variables     test: 9000 obs. of  40 variables      Using 75% of `train` data, I could create a model that achieves 93% of success in classifying the `test` data. Is it possible, as I have seen suggested here, to improve the prediction results by adding prediction on the test corpus to the training corpus? (source wrote: "predict the labels of the test corpus and use them as part of the training corpus"). Intuitively I assume that since I do not have the class results for `test`, this means there is no new information to derive from the test data. But the source mentioned above seem to state that it does improve prediction results. Can it be true, and how would you explain that?


title: Does the sandwich estimator in GEE protect against both correlation misspecification and heteroscedasticity?

text: The relative merits of GEE with exchangeable correlation or GEE with independence and the sandwich estimate have been discussed, but I couldn't find a post specifically addressing my question. I have analyzed longitudinal data using GEE with AR-1 correlation structure and the sandwich/robust estimates of the standard errors. In the case of longitudinal data, the sandwich estimates are used to protect against miss-specifications of the correlation structure. The AR-1 was chosen based on the auto-correlation function of the data since it will produce a more efficient estimate than assuming independence. However, the univariate distrubitions of the outcomes are positively skewed and the residuals from the models exhibit some heteroscedasticity. For GEE with longitudinal (or more generally, clustered) data, does the sandwich estimate protect against heteroscedasticity (what it was developed to do) in addition to its intended use to protect against miss-specificaiton of the correlation structure? Does reducing the heteroscedastcity with a log transformation of the outcome offer more efficiency for the estimate or some other benefit? I would especially appreciate a reference to a journal article addressing this issue.


title: Exponential family where set of natural parameters has empty interior

text: In my math-stat class we have a theorem that goes: Let $\{P_\theta : \theta \in \Theta\}$ be a $k$ parameter exponential family (i.e. the density of a member of this family can be written as $f(\mathbf{x}) = h(\mathbf{x})e^{\sum_{i=1}^k c_i(\theta)T_i(\mathbf{x}) - B(\theta)}$), suppose $\{c_1(\theta),\ldots,c_k(\theta)\}$ (what I am calling the natural parameters) has a **non-empty** interior. Then $T(\mathbf{X}) = \{(T_1(\mathbf{X}),\ldots, T_k(\mathbf{X})\}$ is complete and sufficient. This theorem is often used to identify UMVUE (Uniform minimum variance unbiased estimators) using the Lehmann-Scheffe Theorem. I haven't run into any good examples of an exponential family where the set of natural parameters has an empty interior, can you guys provide a good example?


title: Using R - lme() with Log dataset

text: this is a question regarding using **lme()** in **R** to calculate between- group and within-group variation in a dataset containing _log values_. I have read that when using lme() to calculate sigma (variation) values in a dataset, if the values in the dataset (read in R from a .csv file and stored as a data frame) are log, they are to be logged before using lme(); i.e., if we have a data frame with x, y values like: x= (1, 1, 2, 2, 3, 3) and corresponding y values: y= (4.5e-8, 4.2e-7, 3.2e-7, 3.1e-7, 3e-5, 3.4e-5). For such a dataset, we log all the Y-values _before using lme()_ on this data frame. Now, my question is, if we wish to plot this x, y data and its mean on a chart and use the sigma values to calculate their max. and min. deviation, do we log the mean and original dataset as well or modify the sigma values that R returns ? Currently, I am plotting the original values ( _not_ log) as a chart and am calculating their mean as well **before** logging them. Now, this results in the sigma values being much higher than the actual dataset. Should I not log the values at all before using lme() ? Or, alternatively, should I modify the returned sigma values (take their Inverse or something) and use that in my calculation as that would give a more realistic value with respect to the dataset and its mean ? If anyone has any suggestions, kindly let me know soon. Thanks.


title: First-differences versus endline analysis with an RCT

text: Quick question that I hope someone can help with. Imagine you have data from a randomized controlled trial (RCT) and you are interested in examining the impact of the treatment on a test score. Say that you administer the same test to participants prior to the treatment and after the treatment. I have presented means and standard deviations for each group in the table below. The control group is scoring a little higher on the pre-treatment test (but the difference between the control and treatment is not statistically significant. i.e., we are confident that the treatment indicator is an exogenous variable). We then examine the test scores after the treatment and find that the treatment effect is not statistically significant. However, if we take first differences, we find that the treatment effect is statistically significant. Which results should one trust? The endline only analysis or the first- difference analysis? * * *                               | Treatment   |  Control   | H0: mean(Treatment) - mean(Control) = 0       ------------------------------------------------------------------------------------          Before        |   60.61     |   63.41    |  Fail to reject                     |  (16.13)    |  (17.86)   |       After         |   64.02     |   62.80    |  Fail to reject                     |  (11.42)    |  (14.76)   |     After - Before  |    3.41     |   -0.61    |  Reject                                 |  (13.99)    |  (15.10)   |       The three test that are being carried out are :   1) $Treatment_{pre}$ vs $Control_{pre}$ (pre-treatment scores). This is to check if the groups are equivalent prior to receiving treatment.   2) $Treatment_{post}$ vs $Control_{post}$ (post-treatment scores). This is to check if the treatment has had an effect.   3) $(Treatment_{post} – Treatment_{pre})$ vs $(Control_{post} – Control_{pre})$. This is an alternative way of checking if the treatment is having an effect. Given that we fail to reject the null in test (1), I'm wondering whether test (2) or test (3) is more appropriate for testing whether the treatment is having an effect?


title: Poisson counting process parameter

text: Two quick questions:   1. What's the maximum likelihood estimator of the parameter of an homogeneous Poisson counting process?   2. To estimate $\lambda$ I'm currently using number of events/total time, $N(t)/\Delta t$. What distribution does this statistic follow?


title: Testing equality of proportions

text: I am comparing three proportions for age categories for categorical variables that have also categories. I used Chi-square test but I found that SPSS gave me a note which is for example: a. 2 cells (33.3%) have expected count less than 5. The minimum expected count is .23. What is the solution>


title: Cross-correlation in Matlab

text: What is the difference between:   * 1) xcov(x,y,10,'unbiased')/sqrt(xcov(x,x,0,'unbiased')*xcov(y,y,0,'unbiased'));   * 2) xcorr(x,y,10,'unbiased');   * 3) [ **A** , B] = crosscorr(x,y,10); ? I think (but I am not sure) that "crosscorr" removes the means and "xcorr" doesn't, but i don't get why the first one gives a completely different answer compared to the others. Am I using a wrong formula in 1)? Thank you very much


title: Geometric interpretation of generalized linear model

text: For linear model $y=x*\beta+e$, we can have a nice geometric interpretation of estimated model via OLS: $\hat{y}=x*\hat{\beta}+\hat{e}$. $\hat{y}$ is the projection of y onto the space spanned by x and residual $\hat{e}$ is perpendicular to this space spanned by x. Now, my question is: is there any geometric interpretation of generalized linear model (logistic regression, Poission, survival)? I am very curious about how to interpret the estimated binary logistic regression model $logit(\hat{p})=x*\hat{\beta}$ geometrically, in a similar way as linear model. It even does not have an error term. I found one talk about geometric Interpretation for Generalized Linear Models. http://statweb.stanford.edu/~lpekelis/talks/13_obs_studies.html#(7). Unfortunately, the figures are not available and it is quite hard to picture. Any help, referencing, and suggestion will be greatly appreciated!!!


title: New output fields of glm.summary (in R)?

text: When I do a (logistic) regression in R, I run something like this:               mydata <- read.csv("data.csv")     mylogit <- glm(a ~ c+d, data = mydata, family="binomial")     summary(mylogit)      As of a few months ago, the output for the coefficients might look like this:               Call:     glm(formula = a ~ c + d, family = "binomial", data = mydata)     ...     Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  -1.6476     0.1898  -8.680  < 2e-16 ***     c             2.4558     0.3414   7.194 6.29e-13 ***     d             2.3783     0.4466   5.326 1.01e-07 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      Trying it today (with a newer version of R), the output looks like the following:               Call:     glm(formula = a ~ c + d, family = "binomial", data = mydata)     ...     Coefficients: (1 not defined because of singularities)                 Estimate Std. Error z value Pr(>|z|)         (Intercept)  -1.6709     0.1924  -8.683  < 2e-16 ***     c1            2.4961     0.3476   7.181 6.94e-13 ***     cc           18.2370   979.6100   0.019    0.985         d1            2.4524     0.4630   5.296 1.18e-07 ***     dd                NA         NA      NA       NA         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      What do the "c1", "cc", etc fields mean? I can't seem to find this any documentation, but perhaps I am looking in the wrong places?


title: Understanding dummy (manual or automated) variable creation in GLM

text: If a factor variable (e.g. gender with levels M and F) is used in the glm formula, dummy variable(s) are created, and can be found in the glm model summary along with their associated coefficients (e.g. genderM) If, instead of relying on R to split up the factor in this way, the factor is encoded in a series of numeric 0/1 variables (e.g. genderM (1 for M, 0 for F), genderF (1 for F, 0 for M) and these variables are then used as numeric variables in the glm formula, would the coefficient result be any different? Basically the question is: does R use a different coefficient calculation when working with factor variables versus numeric variables? Follow-up question (possibly answered by the above): besides just the efficiency of letting R create dummy variables, is there any problem with re- coding factors as a series of numeric 0,1 variables and using those in the model instead?


title: If $X$ is lognormally distributed, what is the distribution of $1 / (1 + X)$?

text: Let $X$ be lognormal with parameters $\mu$ and $\sigma$ (such that $\log(X)$ is Gaussian with mean $\mu$ and variance $\sigma^2$). What is the distribution of $1 / (X + 1)$? I am wondering whether it is a "simple" parametric distribution. If the "+1" were replaced with zero, the problem becomes quite easy.


title: Robust regression and Sandwich estimators

text: Can you give me an example of the use of sandwich estimators in order to perform a robust regression? I can see the example in `?sandwich` but I don't quite understand how we can go from               lm(a ~ b, data) # R-coded      to an estimate and a p.value resulting of a robust regression using the variance-covariance matrix returned by the function `sandwich`.


title: What is the distribution of an outer product of random values mod P?

text: I haven't been able to answer this daunting question for days now and I am not able to find the right resource for it, so maybe someone here might help me out. Consider an outer product of two **uniform** random vectors **mod p** , where the vector consists of independent uniform random values modulo a prime number $p$. Vectors $A$ and $B$ are independent of each other. $$ \begin{pmatrix} a_0 & a_1 & a_2 &.. &a_n\end{pmatrix} \mod p$$ and $$ \begin{pmatrix} b_0 & b_1 & b_2 &.. &b_n\end{pmatrix} \mod p$$ such that the outer product is $$ \begin{pmatrix} a_0b_0 & a_1b_0 & a_2b_0 &...\\a_0b_1 & a_1b_1 & a_2b_1 &...\\a_0b_2 & a_1b_2 & a_2b_2 &...\\... &... &... &...\end{pmatrix} \mod p$$ What is the distribution of this outer product? Are the values uniform or not and are they dependent or not? Can we also get it to be uniform and independent (w.h.p) by using some other distribution for the random vectors or maybe by adding up more (k) random outer products?


title: Biological data, 3x3 independent assays

text: this is an example experiment design for which I would need some help (the values are the same between assays but should be different). ![enter image description here](http://i.stack.imgur.com/IC9Ey.png) I want to test if the use of drug 1 or drug 2 has a significant effect on the output being measured. The control group was not exposed to any drug. The assay is then repeated 2 more times for statistic significance. Each assay is done at different days with different cells and media. The problem comes when a statistical test is to be chosen. Here are my thoughts:   * For each assay individually, ANOVA should not be used because I don't want to compare 3 groups, I want to compare the treated groups with the one control group. So, I would do t-test for Control Vs. D1 and Control vs. D2. Is this correct?  But then, there's the 3 assays considered all together and the representation of the data. Doing a normalization, with Control being 100% and doing the mean of the 3 assays according to the % is a good way to show the data but statistics wise I don't know how to do this. Pooling all measurements for each group and calculate the mean (doing then the t-tests again) doesn't seem appropriate because even the control group usually shows variations. Normalizing, with Control being 100%, seems even wronger, as the Control group(s) won't have any deviation values. I thought about introducing the Mean-SD-N for each 3 assays but Graphpad doesn't seem to accept that for t-test and 1wayANOVA.


title: Prove that serial correlation in error terms makes $X_t e_t$ a non-MDS

text: > Given the simple regression $$Y_t = \beta X_t + e_t$$, with $$e_t = > \epsilon_t + \theta_1 \epsilon_{t-1} + \theta_2 \epsilon_{t-2}, \epsilon_t > \sim WN(0, \sigma^2)$$ and $$X_t iid, X_s \perp e_t \forall s,t$$ Is ${X_t > e_t}$ a Martingale Difference Sequence w.r.t $F_{t-1}=(X_{t-1}, e_{t-1}, > X_{t-2}, e_{t-2}, ...)$? This is my attempt so far: $E(X_t e_t | F_{t-1}) = E(X_t E(e_t|F_{t-1}, X_t)|F_{t-1})$ In the serially uncorrelated case, $E(e_t|F_{t-1}, X_t)=0$, leading to a MDS. In this serially correlated case, I suspect that it will be different from 0 but cannot calculate the expectation exactly. Ultimate it boils down to finding $E(\epsilon_{t-1}|F_{t-1})$ and $E(\epsilon_{t-2}|F_{t-1})$, which is where I am stuck.


title: What are the assumptions for repeated-measures MANCOVA?

text: Are these assumptions correct and the list complete for **repeated-measures MANCOVA**?   * sphericity   * multivariate normality   * homogeneity of variance   * homogeneity of covariance matrices   * homogeneity of regression slopes   * independence of independent variable and covariate


title: is it possible to use partitioned data(train&test) together with cross-fold validation?

text: I have used SPSS Clementine in order to train a classifier, for this I have used a partition node with 2 parts(train and test),then using a c5-tree and cross-fold validation. I did this because I think SPSS separates training data and does a 10-fold validation to build a best model and then test it by test data. would you please verify or criticize this? this is the explanation of SPSS about how cross-fold works: "Cross-validate. If this option is selected, C5.0 will use a set of models built on subsets of the training data to estimate the accuracy of a model built on the full dataset. This is useful if your dataset is too small to split into traditional training and testing sets. The cross-validation models are discarded after the accuracy estimate is calculated. You can specify the number of folds, or the number of models used for cross-validation." And if my idea is not true(to use both partitions and cross-validation) what has happened in SPSS? Is there any priority or not? thanks all


title: transformation to normality of the dependent variable in multiple regression

text: Is it really important to normalize dependent variables in multiple regression or are there any exceptions? My model is providing better results with more significant hypothesis when the DVs are not normalized (transformed). Will be grateful for the comments.


title: Convert Likert scale to pairwise?

text: I want to find a way to convert the Likert scale to Pairwise (AHP/ANP) scale, because my questionnaire contains 30 criteria to compare, and making a pairwise questionnaire will contain (30*(30-1))/2 questions, which are very long to answer for the respondents. Does anyone know how to deal with that? or any other method to weigh my criteria?


title: Probability in Normal Distribution

text: I've a question here. Basically I have the mean and standard deviation of a variable (IQ). I want to find the probability that 8 out of 10 people randomly selected have IQ scores between 85.0 and 122.5. I also have the following information worked out:   * mean = 100   * standard distribution = 15   * n = 10   * p = 0.7745 Should I be using Poisson to do this? Appreciate any directions and advice please.


title: How to model football (soccer) scores

text: I'm trying to model football scores using the poisson distribution. I'm calculating a lambda for each team and then multiplying the probabilities of two poisson calculations to get the probability of a specific score. e.g. poison(lambdaHome, 1)*poisson(lambdaAway, 2) to calculate the probability of a 1-2 result at the end of the match. The problem I'm having is that the poisson distribution assumes that the second goal is just as likely as the third goal is just as likely as a fourth goal. But in reality this shouldn't be true. If a match has a score of 2-0 then intuitively it is less likely for the team to score another goal is it would probably spend more time on defence to keep the score than it would do attacking to score more goals. Is there a probability distribution similar to poisson, but rather than having a constant probability for all goals, makes each goal less and less likely?


title: Normality assumption in linear regression

text: As an assumption of linear regression, the normality of the distribution of the error is sometimes wrongly "extended" or interpreted as the need for normality of the y or x. Is it possible to construct a scenario/dataset that where the X and Y are non- normal but the error term is and therefore the obtained linear regression estimates are valid?


title: Calculating a popularity score for news articles based on user comments

text: My research project is to come up with a tool for analyzing publication and subscription patterns of Sri Lankan news web sites. As a part of my project, my tool should be capable of finding out the followings, in addition to analyzing the publication patterns of news websites. • Mostly subscribed news topics (most popular news topics) based on user comments. Usage of user comments to find out the popularity was a requirement of the Supervisor of the project. • News topics which increase popularity with respect to time (for example: fewer comments were made on university strike, but all of a sudden the number of comments increases) According to the implementation, my tool extracts news articles and the associated user comments on a daily basis from several news websites. The tool continuously check whether each article will further receive comments by tracking the article for “X” number of days (“X” is a user defined value and can be changed by the user). I have collected following statistics for each article that I have extracted from Month April to present. • Published date of the article • Number of days that the article has received comments • Number of comments (excluding replies) made by the readers. • Number of replied comments made by the readers. • Number of unique readers (nick names) who have made comments and reply comments • The entire total number of comments and replies per article Is there a better statistical methodology which I could use for calculating a popularity score for each news article, using some of the factors listed above. Instead of calculating the total number of comments that each article has received and sorting them to select the “k” top most articles.


title: Compare dispersion between two distributions with thick tails

text: What tests would I use to see if two distributions are 'significantly' different in dispersion if they are thick tailed, like a t-distribution? thanks


title: Unsupervised Dimensional reduction for mixed data types

text: I have a data set with about 50K rows and 100 columns. You can consider every row to be representing one restaurant. My goal is to calculate dissimilarities between all the restaurants - Gower's coefficient. Of those 100 columns (features), a few of them are numeric data and nominal data. The problem is the other columns (about 90) are very sparse binary data (1/0). I do think that those 90 columns of binary data can be reduced to some smaller number of columns, so that the computational time can be reduced significantly. But I don't know what method I should use to reduce such a large amount of binary data. Can anyone give me some suggestions? It will be most helpful if you can provide me some references and R code.


title: SVM versus Bayesian regression example(s)?

text: I am trying to track down examples where some basic problems have been tackled via both classical Machine Learning algorithms and more formal statistical methods. In particular, I'm interested in examples where SVM has been contrasted with hierarchical Bayesian methods. (but I'll happily take non-SVM examples.) FWIW- I ran across a blog post that suggested that Bayesian stat methods provided additional insight beyond that provided by SVM, but foolishly didn't bookmark.


title: Is extra sum of squares $SSR(X_p|X_1,...X_{p-1})$ in multiple regression always positive?

text: An extra sum of squares $SSR(X_p|X_1,...X_{p-1})$, assuming that no pair of predictor variables are perfectly correlated, measures the marginal reduction in the error sum of squares. Eventually one can view an extra sum of squares as measuring the marginal increase in the regression sum of squares when one or several predictor variables are added to the regression model. Assuming the number of observational data $n$ is lager than the number of predictor variables $p$. Can I have a rigourous proof that $SSTO=\sum_{i=1}^n (Y_i-\bar{Y})^2\geq SSR(X_p|X_1,...X_{p-1})\geq0$ ? * * * Applying LSE $SSE(X_1,...,X_{p-1})=\sum(Y_i- b_0-b_1X_{1,i}-...-b_{p-1}X_{p-1,i})^2=\sum(Y_i- b_0-b_1X_{1,i}-...-b_{p-1}X_{p-1,i}-0\times X_{p,i})^2\geq\sum(Y_i- b_0'-b_1'X_{1,i}-...-b_{p-1}'X_{p-1,i}-b_p'X_{p,i})^2=SSE(X_1,...,X_{p-1},X_p)$ $SSR(X_p|X_1,...X_{p-1})=SSE(X_1,...,X_{p-1})-SSE(X_1,...,X_{p-1},X_p)\geq 0$ Proved.


title: When does a distribution not have a mean or a variance?

text: I believe I read today a phrase which went something like this: > If a distribution has a mean and a variance ... So I guess that means some distributions do not have means or variances? I fiend that a bit difficult to understand - for example: I could get data which is said to follow such a mean-less and variance-less distribution and as an estimation compute the mean and variance of such a data set. could someone help shed some light on what I am missing?


title: How do I test if two random variables are independent

text: I'm trying to understand the concept of independence experimentally using computer simulation. I imagine to understand this concept:   1. I should be able to generate three data sets from random Variables A, B and C where A and B are independent and A and C are not independent.   2. I should be able to run some test where I compare data from A and B (and A and C) and the test gives me a measure of independence. Could I get pointers of a resource (paper, textbook chapter, code, blog post, etc.) where there is sufficient information for me to experimentally understand the concept of independence.


title: Friedman's test and post-hoc analysis

text: I am carrying out a statistical analysis for my research. I am using a Friedman's test with post hoc analysis. At present I am using the function `friedman.test.with.post.hoc` available for R software. This function reports "maxT" value instead of the chi-square value. Someone can explain me what is maxT and its relationship with Chi-square? Thank you.


title: How to remove trend with no look ahead bias?

text: I would like to explore the different ways one can **detrend** a time series **without look ahead bias**. I wanted to use the Hodrick Prescott filter, which seems like a quite good frequency filter, but it is based on an optimization method, and I understand that it may give strange and volatile results at the border. Wavelet smoothing on a rolling window would be another option, but again border effects can be huge (the data is copied by symmetry which is horrible for the precision of the technic at the edge). Any idea or comments? PS: The subject has already been discussed here, I know. But I would like to dig a bit more on a more precise question.


title: Variance partitioning - why be cautious?

text: I'm about to use variance partitioning to interpret my results of a given model and across models and have come across various criticisms of it most notably by Pedhazur (1982, 1997). Also, the criticisms are of both the approaches to VP - commonality analysis and incremental partitioning of variance. From what I understand the former is not as "dangerous" as the latter. In short, the criticisms are related to the high dependence on R^2 and multicollinearity. Furthermore, the criticism is about the how meaningful it is to attribute explanatory power to individual variables in a multivariate situation. Since I will be using this model for the first time, and even after reading other posts on this and other forums/literature, I'm having trouble wrapping my head around the criticisms. Any clear explanation of the criticisms would be appreciated.


title: Likelihood ratio test vs. $\chi^2$/Z-test for comparing binomial datasets

text: When does/can one use the likelihood ratio significance test instead of Fisher's exact test or its Pearson $\chi^2$ approximation for comparing two binomial datasets? Given two binomial datasets (distributions), I'm seeing the LR test being used to compare one distribution against the global (combined) distribution. Usually I apply Fisher's test for comparing one dataset against the other. I realize that LR testing is of the Neyman-Pearson school, which assumes a fully specified alternative model as well as null model. E.g., in the LR test Wikipedia page example, it's being used to compare two binomial datasets (# heads/tails for two coins). Why not use the $\chi^2$ test to compare the two samples against each other? What are the conceptual differences in these two approaches? When do I use which? And when is it appropriate to compare one sample against not its complement but the global dataset?


title: Using an OLS coefficient to estimate a non-linear coefficient

text: Using OLS, I've estimated the following equation: $y_i = \alpha_0 + \alpha_1 X_i + \epsilon_i$ I know that theoretically, the following should be true: $y_i = a + (1-e^{-\lambda 60}) X_i$ Is there any way, having an estimate of $\alpha_1$ I can translate it to an estimate of $\lambda$? As a follow up, if this is not possible without some difficulty, If I knew the distribution of $\alpha_1$ was a normal distribution with some mean and variance, is there a way to describe what the form of the distribution of $\lambda$ would be? I feel like it would be a log-linear distribution, but I'm not sure what the mean/variance would be.


title: What is perplexity?

text: I came across term **perplexity** which refers to the log-averaged inverse probability on unseen data. Wikipedia article on **perplexity** doesnt give a intuitive meaning for the same. This perplexity measure was used in pLSA paper. Can anyone explain the need and intuitive meaning of **perplexity** measure?


title: Unsupervised Bayesian naive Bayes

text: I'm reading a paper _Gibbs sampling for the uninitiated_. In this paper, the authors try to use Gibbs sampling for a bayesian naive bayes model. They formalize the model as a graphical model in page 8. And in the example, they are trying to predict the emotion(sentiment) of a document. However, what I don't understand is that, they claim without label $L$, using Gibbs sampling could still sample all the parameters needed, including $L$. I'm not sure how should I interpret this. Without training label, it's essentially a clustering problem, but if not using labels, how should we interpret the learnt label $L$? Thanks in advance.


title: Are there any inverse distribution graph that looks like this?

text: I want to generate random numbers according to a distribution curve like the inverse cummulative normal curve. But the curve should allow parameters to adjust the parameters as shown in the image, namely, the slopes (A, B, C), min and max values, between 0.0 and 1.0 in a meaningful manner like in terms of variance, mean etc. the usual stuff. Is there a well known solution for this curve in terms of statistical distribution? ![Inverse Distribution Graph](http://i.stack.imgur.com/0KCbo.jpg)


title: Calculate log-likelihood "by hand" for generalized nonlinear least squares regression (nlme)

text: I'm trying to calculate the log-likelihood for a generalized nonlinear least squares regression for the function $f(x)=\frac{\beta_1}{(1+\frac x\beta_2)^{\beta_3}}$ optimized by the `gnls` function in the R package `nlme`, using the variance covariance matrix generated by distances on a a phylogenetic tree assuming Brownian motion (`corBrownian(phy=tree)` from the `ape` package). The following reproducible R code fits the gnls model using x,y data and a random tree with 9 taxa:               require(ape)     require(nlme)     require(expm)     tree <- rtree(9)     x <- c(0,14.51,32.9,44.41,86.18,136.28,178.21,262.3,521.94)     y <- c(100,93.69,82.09,62.24,32.71,48.4,35.98,15.73,9.71)     data <- data.frame(x,y,row.names=tree$tip.label)     model <- y~beta1/((1+(x/beta2))^beta3)     f=function(beta,x) beta[1]/((1+(x/beta[2]))^beta[3])     start <- c(beta1=103.651004,beta2=119.55067,beta3=1.370105)     correlation <- corBrownian(phy=tree)     fit <- gnls(model=model,data=data,start=start,correlation=correlation)     logLik(fit)       I would like to calculate the log-likelihood "by hand" (in R, but without use of the `logLik` function) based on the estimated parameters obtained from `gnls` so it matches the output from `logLik(fit)`. NOTE: I am not trying to estimate parameters; I just want to calculate log-likelihood of the parameters estimated by the `gnls` function (although if someone has a reproducible example of how to estimate parameters without `gnls`, I would be very interested in seeing it!). I'm not really sure how to go about doing this in R. The linear algebra notation described in Mixed-Effects Models in S and S-Plus (Pinheiro and Bates) is very much over my head and none of my attempts have matched `logLik(fit)`. Here are the details described by Pinheiro and Bates: The log-likelihood for the generalized nonlinear least squares model $y_i=f_i(\phi_i,v_i)+\epsilon_i$ where $\phi_i=A_i\beta$ is calculated as follows: $l(\beta,\sigma^2,\delta|y)=-\frac 12 \Bigl\{ N\log(2\pi\sigma^2)+\sum\limits_{i=1}^M{\Bigl[\frac{||y_i^*-f_i^*(\beta)||^2}{\sigma^2}+\log|\Lambda_i|\Bigl]\Bigl\}}$ where $N$ is the number of observations, and $f_i^*(\beta)=f_i^*(\phi_i,v_i)$. $\Lambda_i$ is positive-definite, $y_i^*=\Lambda_i^{-T/2}y_i$ and $f_i^*(\phi_i,v_i)=\Lambda_i^{-T/2}f_i(\phi_i,v_i)$ For fixed $\beta$ and $\lambda$, the ML estimator of $\sigma^2$ is $\hat\sigma(\beta,\lambda)=\sum\limits_{i=1}^M||y_i^*-f_i^*(\beta)||^2 / N$ and the profiled log-likelihood is $l(\beta,\lambda|y)=-\frac12\Bigl\{N[\log(2\pi/N)+1]+\log\Bigl(\sum\limits_{i=1}^M||y_i^*-f_i^*(\beta)||^2\Bigl)+\sum\limits_{i=1}^M\log|\Lambda_i|\Bigl\}$ which is used with a Gauss-Seidel algorithm to find the ML estimates of $\beta$ and $\lambda$. A less biased estimate of $\sigma^2$ is used: $\sigma^2=\sum\limits_{i=1}^M\Bigl|\Bigl|\hat\Lambda_i^{-T/2}[y_i- f_i(\hat\beta)]\Bigl|\Bigl|^2/(N-p)$ where $p$ represents the length of $\beta$. I have compiled a list of specific questions that I am facing:   1. What is $\Lambda_i$? Is it the distance matrix produced by `big_lambda <\- vcv.phylo(tree)` in `ape`, or does it need to be somehow transformed or parameterized by $\lambda$, or something else entirely?   2. Would $\sigma^2$ be `fit$sigma^2`, or the equation for the less biased estimate (the last equation in this post)?   3. Is it necessary to use $\lambda$ to calculate log-likelihood, or is that just an intermediate step for parameter estimation? Also, how is $\lambda$ used? Is it a single value or a vector, and is it multiplied by all of $\Lambda_i$ or just off-diagonal elements, etc.?   4. What is $||y-f(\beta)||$? Would that be `norm(y-f(fit$coefficients,x),"F")` in the package `Matrix`? If so, I'm confused about how to calculate the sum $\sum\limits_{i=1}^M||y_i^*-f_i^*(\beta)||^2$, because `norm()` returns a single value, not a vector.   5. How does one calculate $\log|\Lambda_i|$? Is it `log(diag(abs(big_lambda)))` where `big_lambda` is $\Lambda_i$, or is it `logm(abs(big_lambda))` from the package `expm`? If it is `logm()`, how does one take the sum of a matrix (or is it implied that it is just the diagonal elements)?   6. Just to confirm, is $\Lambda_i^{-T/2}$ calculated like this: `t(solve(sqrtm(big_lambda)))`?   7. How are $y_i^*$ and $f_i^*(\beta)$ calculated? Is it either of the following: `y_star <\- t(solve(sqrtm(big_lambda))) %*% y` and `f_star <\- t(solve(sqrtm(big_lambda))) %*% f(fit$coefficients,x)` or would it be `y_star <\- t(solve(sqrtm(big_lambda))) * y` and `f_star <\- t(solve(sqrtm(big_lambda))) * f(fit$coefficients,x)` ? If all of these questions are answered, in theory, I think the log-likelihood should be calculable to match the output from `logLik(fit)`. Any help on any of these questions would be greatly appreciated. If anything needs clarification, please let me know. Thanks! **UPDATE** : I have been experimenting with various possibilities for the calculation of the log-likelihood, and here is the best I have come up with so far. `logLik_calc` is consistently about 1 to 3 off from the value returned by `logLik(fit)`. Either I'm close to the actual solution, or this is purely by coincidence. Any thoughts?                 C <- vcv.phylo(tree) # variance-covariance matrix       tC <- t(solve(sqrtm(C))) # C^(-T/2)       log_C <- log(diag(abs(C))) # log|C|       N <- length(y)       y_star <- tC%*%y        f_star <- tC%*%f(fit$coefficients,x)       dif <- y_star-f_star         sigma_squared <-  sum(abs(y_star-f_star)^2)/N       # using fit$sigma^2 also produces a slightly different answer than logLik(fit)       logLik_calc <- -((N*log(2*pi*(sigma_squared)))+            sum(((abs(dif)^2)/(sigma_squared))+log_C))/2


title: Threshold models and flu epidemic recognition

text: I'm fooling around with threshold time series models. While I was digging through what others have done, I ran across the CDC's site for flu data. http://www.cdc.gov/flu/weekly/ About 1/3 of the way down the page is a graph titled "Pneumonia and Influenza Mortality....". It shows the actuals in red, and two black seasonal series. The top seasonal series is labeled "Epidemic Threshold" and appears to be some constant percent/amount above the "Seasonal Baseline" series. My first question is: Is that really how they determine when to publicly say we're in an epidemic (some percent above baseline)? It looks to me like they're in the noise range, not to mention the "other factors" influence that is obviously not accounted for in that baseline series. To me, there are way too many false positives. My second question is: Can you point me to any real world examples/publications of threshold models (hopefully in R)?


title: Test for relative proportions

text: This is my first question here. I have experimental results for two datasets where I want to check if the difference is significant or not. Here is a description of my experiment: Experiment: I have a set of 70 genes where I showed that 12 of them are evolving in a particular manner. Now I wanted to check if these 70 genes are biased to carry that evolutionary characteristics compared to all the genes in that species. Unfortunately, due to experimental restriction I could not replicate this experiment on all the genes. So I selected a random set of 100 genes (without replacement) and got the results that 8 out of 100 genes were carrying the same characteristics. Now to check if this difference is statistically significant or not, I initially used binomial test on the data (using R prop.test). But someone told me that this is not the right test for your data. Can someone please help me with the correct test Thanks, RT


title: How to use Confidence Intervals to find the true mean within a percentage

text: I'm working through a practice problem for my Stats homework. We're using Confidence Intervals to find a range that the true mean lies within. I'm having trouble understanding how to find the required sample size to estimate the true mean within something like +- 0.5%. I understand how to work the problem when the range is given as a number, such as +- 0.5 mm. How do I handle percentages?


title: Decision to add all predictive variables in Kaplan-Meier and Cox PH

text: In R's `survival` and `rms` packages, Kaplan-Meier (`Surv` and `survift`) and Cox Proportional Hazards (`cph`) are implemented. I have been experimenting with these two by adding one-by-one the predictors I have in my dataset. These are not time-dependent meaning for one record, there is only one observation for these predictors which do not depend on a time basis. The format for KM requires a `Surv` object to be made first before given to the `survfit` function for further processing.               fit      <- Surv(time,status)~X1+X2+X3+ ... +Xn     model.km <- survfit(fit)      `cph` also uses a similar way to create a Cox model.               model.cox = cph(fit, surv=TRUE, x=TRUE, y=TRUE)      * * * Do I have to add all of my predictors in the dataset into the `Surv` in order to get a legitimate or wholesome representation of my data? I always try to add each one at a time into the functions and look out for ones that will cause an error which I remove to go further the creation. If one of these predictors are not welcomed into `Surv`, `survfit` and `cph`, I instantly remove them from the creation and not the dataset. Is this method that I do correct?


title: How can I compute regression for several longitudinal data sets (thus, with auto-correlated error)?

text: My actual project is a bit complicated, but I'll explain by analogy (which I hope facilitates response): I have 3 substances, say water, motor oil, and ethanol. For each substance, I have 5 samples in a beaker (total 15 beakers). I heat all the beakers on a hot-plate up to 70 degrees Celsius, and over the next hour, I measure the temperature of the fluid in each beaker at 5 minute intervals. Newtonian cooling provides me with a good prediction about these temperature data, namely that the temperature of the fluid in each cup should follow an exponential distribution: y = a + e^(-kt) where a is room temperature. I want to estimate the value k for each substance and test the hypothesis that k1 > k2 > k3 (1, 2, 3 corresponding to my three substances). The natural method of estimating k seems to be computing a non-linear regression on each substance's data, or possibly log-transforming all the data and then just computing a simple linear regression. However, there are problems. Some questions:   1. Given the obvious autocorrelation in the longitudinal data (confirmed by my (P)ACF plots of course), must I compute an AR term and filter my data prior to computing the regression?   2. Assuming I compute this autoregression term, how do I compute it for five independent sets of data (the five beakers of a give substance)? I could average the five beakers together and then compute the regression, but this screws up my AR term (assuming I need one) and also throws off my estimate of the actual within-beaker variance from the model.   3. What completely wrong-headed assumption(s) have I worked in here...?


title: Why is k-fold cross validation a better idea than k-times resampling true validation?

text: I'm currently working through a machine learning textbook and just read a bit about k-fold cross validation, and I am wondering the following. I want to estimate a parameter, e.g. a penalty parameter for a penalized likelihood method. In order to do this, I can do two different things:   1. I sample the training data so that I get $k$ equally large folds, and for each fold I use the other folds as training data to get estimates for $y$ and I compare these estimates with the actual $y$ from the fold in question. This, I do for every interesting choice of my parameter, and choose the parameter which has the least error, averaged over all folds and all members of each fold.   2. I sample the training data so I get 2 equally large sets, one of which I use as training data to predict the error of the other set. For every interesting lambda, I note the average error. Then, I re-sample the data so I get 2 (different) equally large sets, where I repeat the above procedure. I sample $k$ times in total, and average over these to get an estimate to the best parameter. The second approach looks rather naive, and I am wondering if there is something wrong with it. Are there reasons, generally speaking, why one would prefer method 1 over method 2? Are there computational reasons, or even statistical ones?


title: Using IQR to find Mean without outliers

text: I am trying to figure out a away to calculate the average a set of data in the most robust way that I can, while staying accurate (this has to do with billing software). For example, let's say I have the monthly minutes available to me. I need to find a good monthly average in order to figure out how many minutes are used in a day on average. I decided to use IQR in order to eliminate any outliers due to not properly getting a value one month or a month with very high usage. **The following is a list of numbers that I used:** > > 0,0,0,0,100,200,200,200,300,300,300,300,300,300,300,300,300,300,400,400,400,400,400,400,500,500,500,500,500,500,500,500,600,600,600,600,800,1000,1200,1200,1400,1400,1900 These are the values I was able to get: > Median: 400   > q1: 300   >  q3: 600   >  IQR: 300   >  rangeLow: -50   >  rangeHigh: 850   >  Average after IQR: 359.459459459459   >  Regular Average: 497.674418604651 > > rangeLow = median - IQR*1.5   > rangeHigh = median + IQR*1.5 Does this seem like a decently accurate way of calculating this? I was also thinking about excluding 0's since a lot of those occur due to error. Is there a better way for me to go about this?