title: How to test for mediation with a continuous mediator and DV, but a categorical predictor?

text: I want to explore whether attitudes mediate the relationship between employment status and intention to apply to graduate school. (See image below). I planned to use the Sobel Test to explore the mediation, but I'm not sure I can. How does one/how can one test for mediation with a categorical predictor and both a continuous mediator variable and dependent variable? Prefer answers oriented to SPSS. ![enter image description here](http://i.stack.imgur.com/TXQJ5.png)


title: Why doesn't Random Forest handle missing values in predictors?

text: What are theoretical reasons to not handle missing values? Gradient boosting machines, regression trees handle missing values. Why doesn't Random Forest do that?


title: Optimization parameter for classification

text: I do not have enough knowledge about optimization. My problem is simple. Lets say I have 100 classes. Each class contains some instances (images). The feature vector obtained from a particular image depends on a parameter P. I need to learn this parameter automatically by maximizing the classification accuracy. As I have no idea about optimization, could you please suggest me about how I can do this. Thank you.


title: Visualising a linear model with 6 predictors in R

text: I've recently been wondering how to visualise a linear model with >1 predictor. Today I came across a paper that managed to visualise a linear model with 6 predictors (and one response variable) in a scatterplot. The model was similar to this: weight ~ days since 1st Sept + site + age + year + PC1 + arrival date Weight was a continuous response variable. Days since 1st Sept was a integer predictor, site was a factor predictor (2 levels), age was a factor predictor (2 levels), year was a factor predictor (2 levels), PC1 was a continuous predictor and arrival date was a integer predictor. On the y axis of the scatterplot was values of weight predicted by the model. On the x axis of the scatterplot was days since 1st Sept. Eight lines of best fit showed predicted values at each level of the three factor predictors. PC1 and arrival date were 'held constant at their mean values'. The plot looked like this (made using Photoshop): ![enter image description here](http://i.stack.imgur.com/x9Sfe.jpg) I'm hoping someone can recreate a plot similar to this using a dataset installed in R? I'm particularly interested in how to keep two predictors 'constant at their mean values', to enable visualisation of the other predictors. The mtcars dataset might be appropriate.


title: NEGM Implementation in R

text: I am applying some chaos tests in financial time series in R and I am interested to run a NEGM (Nychka, Ellner, Gallant, McAffrey) Test for a largest Lyapunov exponent based on Jacobian methods. Do you know any package that contains this test, or at least some implementation?


title: Inflated p-value after adjusting for covariates, in GWAS

text: I am working on some GWAS (Genome-Wide Association Studies) now. A genome scan was done for all the SNPs, with first 3 principal components adjusted (PCs are used for adjusting ethnicity effect) and the QQ plot looks fine (most p-values lay on the diagonal line with a few signals off the line). However, when I adjusted for clinical covariates (treatment, clinical stage) and rerun the genome scan, the QQ plot looks weird. Usually, when you adjust for covariates, you expect to see the p-values lay under the diagonal line (or fewer hits than unadjusted model). In my model, after adjusting for clinical covariates, the QQ plot is way off the diagonal line and almost all the SNPs are false positives. I don't know how to explain this. One important thing is that some clinical covariates are very very significant in the clinical model (the model with no SNP information), and the p-value goes to 10E-11. I have 250 samples and 750K SNPs in my dataset.


title: Understanding the predictive distribution in gaussian linear regression

text: I'm reading through the Gaussian Process book http://www.gaussianprocess.org/gpml/chapters/RW2.pdf and there's one section here I don't quite understand (page 11). The author says: > "the predictive distribution is given by averaging the output of all > possible linear models wrt the Gaussian posterior" $$ \begin{aligned} p(f_*|x_*,X,y) &= \int p(f_*|x_*,w)~p(w|X,y)~dw \\ &=\mathcal N\left(\frac{1}{\sigma_n^2}x^T_*A^{-1}Xy,~x^T_*A^{-1}x_*\right) \end{aligned} $$ What does this mean, exactly? I understand that the purpose of using Gaussians is to be able to calculate uncertainty for a prediction, but I'm unclear how the "averaging of the output" doesn't end up with just a mean value of the weight. And how were the parameters for the mean and covariance derived?


title: The difference of kernels in SVM?

text: Can someone please tell me the difference between the kernels in SVM:   1. Linear    2. Polynomial    3. Gaussian (RBF)   4. Sigmoid  Because as we know that kernel is used to mapped our input space into high dimensionality feature space. And in that feature space, we find the linearly separable boundary.. When are they are used (under what condition) and why?


title: Which graphs are best suited to visualize multiple pre- and post-test scores between multiple groups?

text: I have a relatively small sample ($N = 64$) with the following variables:   * `ScoreA_I`: continuous pre-test response variable    * `ScoreA_II`: continuous post-test response variable    * `ScoreB_I`: continuous pre-test response variable    * `ScoreB_II`: continuous post-test response variable    * `sex`: 2 level categorical predictor variable    * `location`: 2 level categorical predictor variable  I'd like to visualize differences between the different groups (sex and location) for the scores. I think I probably want one graph for `sex * ScoreA_I-II * ScoreB_I-II`, another for `location * ScoreA_I-II * ScoreB_I- II`, and another one for the combination `sex * location * Scores`. An option would be to calculate the differences of the scores and plot them with error bars. But I'd like to also show the achieved scores for each group. What graphs are best suitable (and can be ideally done in SPSS)?


title: Do the mean and the variance always exist for exponential family distributions?

text: Assume a scalar random variable $X$ belongs to a vector-parameter exponential family with p.d.f. $$ f_X(x|\boldsymbol \theta) = h(x) \exp\left(\sum_{i=1}^s \eta_i({\boldsymbol \theta}) T_i(x) - A({\boldsymbol \theta}) \right) $$ where ${\boldsymbol \theta} = \left(\theta_1, \theta_2, \cdots, \theta_s \right )^T$ is the parameter vector and $\mathbf{T}(x)= \left(T_1(x), T_2(x), \cdots,T_s(x) \right)^T$ is the joint sufficient statistic. It can be show that the mean and the variance for each $T_i(x)$ exist. However, do the mean and the variance for $X$ (i.e. $E(X)$ and $Var(X)$) always exist as well? If not, is there an example of an exponential family distribution of this form whose mean and variable do not exist? Thank you.


title: Establishing a strategy to log data processing steps

text: I am building a database schema for a web application that will collect user feedback on restaurants. As the collection of user feedback grows for each restaurant, I will be analyzing it to assign properties to the restaurant record that will be exposed publicly in the web app. Some example processing algorithms I plan on implementing are:   * For each restaurant with X "up votes" from independent IP addresses, set property RECOMMENDED to TRUE.    * Merge restaurant records with equivalent NAME and ADDRESS fields. More important than those examples, however, is the fact that I don't know exactly how I will process the data. I'm sure that once I begin receiving user feedback, I'll want to tweak my processing routines, or establish new processing steps. It strikes me as a good idea to include in each restaurant record a list of all processing actions taken against it. **Thus, my question:** **Are there papers or other resources that you could direct me towards that discuss best practices for how to structure logs for data processing actions in a flexible, useful structure?** I've looked around on Google, but all the resources I could find discussed how to process log data, rather than logging the process. Thanks in advance.


title: How do I test for equal variances in a nested design?

text: I've got a dataset where stations (3) are nested within beaches (5). Every station is replicated 6 times per beach, so that I have a total of 90 observations. I now want to test for equal variances, so I chose the Bartlett test in R. I somehow got myself confused; what would be my reference variable for my variances? The replicates? Stations, beaches? Is it correct to assume that I have to test the variance of stations nested in beaches?    **Output of Bartlett test in R:** > bartlett.test(Items_m2 ~ Beach/Station, data=beachlitter)               Bartlett test of homogeneity of variances     data:  Items_m2 by Beach by Station      Bartlett's K-squared = 20.6793, df = 4, p-value = 0.0003666         **Screenshot of the data** (where transect is treated as replicates) ![enter image description here](http://i.stack.imgur.com/DWsf8.png)   **Dataset to download (.csv)** https://drive.google.com/file/d/0B5vU4bWlwahrcTZSLVNIaTlDQ0U/edit?usp=sharing


title: Efficient estimators

text: let $x_{1},..,x_{n}$ be iid with cauchy distribution. Let $y=median(x_{1},..,x_{n})$ , $z=MLE$. I get that $y$ is not an asymptotically efficient estimator to the location of a cauchy distribution because the asymptotic approximation to the distribution of the sample median has variance ${\pi}^2/4$ and minimum variance is $1/2$. Is $z$ asymptotically efficient estimator?


title: regression with kNN on dataset with categorical variables

text: I am trying to train a regression model for dataset with 500k observations and 3 features. The features are categorical and have 50, 50 and 100 levels. Is (generally) kNN appropriate for this kind of task? I am using R. I tried to turn my categorical variables into dummy variables but I end up with very large and sparse data set. I am using data.matrix for conversion and it sets the matrix to double by default. Is there a way to set it to logical instead?


title: Backshift operator property not clear

text: In my introductory book on time series analysis the backshift operator $\mathbf{B}$ is introduced using the following definition: $$ \mathbf{B}x_t=x_{t-1} $$ Then the author sets off to derive some properties of random walks, and there is this step, that is a complete mystery to me: $$ x_t=(1-\mathbf{B})^{-1}w_t\Rightarrow x_t=(1+\mathbf{B}+\mathbf{B}^2+\mathbf{B}^3+\dots)w_t $$ I do not understand how $(1-\mathbf{B})^{-1}$ is equivalent to $(1+\mathbf{B}+\mathbf{B}^2+\mathbf{B}^3+\dots)$. Maybe in the book the author did not provide some important properties of the backshift operator that were used deriving the above step.


title: Cochran's Q & McNemar tests together

text: I have binomial outcome data for three groups that I am analyzing with Cochran's Q test. As discussed at this site, Cochran. However, I have to repeat this seven different times because I am analyzing performance of three diagnostic tests to see if they can find resistance to seven different drugs. When there are significant differences between the three groups, I perform a series of three paired post-hoc McNemar tests to further analyze the data. However, my question is when there are not significant differences, do I still perform the three paired post-hoc tests ? Part of me says why not, these McNemar's are valid tests. However, another part of me figures that if this is the case, why do the Cochran's Q test at all ? Could some one point me in the correct direction ? Thanks !


title: Recursive feature extraction with random forests following by statistical inference?

text: I have a small dataset of n = 200 with about 80 predictors. The outcome I am working with is binary. I'm interested in doing a two-step process where: (1) I conducted recursive feature extraction using random forests to pull out the predictors that maximize predictive ability. (2) Examining that subset of predictors more closely. I'd like to find the direction of their association with the outcome and if possible the magnitude of their effect size. I don't think random forests can do this, right? Would it make sense to run standard logistic regression on the subset of variables that the RF pulled out in step one? So ultimately I'm not interested in prediction. I'd like to use machine learning to select a subset of features to examine more closely, but at the end of it I'm interested in inference and effect estimation (in what direction is predictor 31 associated with the outcome). Any ideas?


title: ABC: Why not use the distance measure as a pseudo-likelihood instead?

text: I've read about the ABC rejection algorithm when not being able to calculate the likelihood directly, and my question is: if we have to introduce a distance measure $\rho(D,D')$ anyways, why not use that measure as a pseudo- likelihood to weight the $\theta$ that generated $D'$ instead of thresholding on an arbitrary value of $\epsilon$? It seems like this would be much more efficient in high dimensional data spaces where you are unlikely to 'hit' close to your original dataset very often. I realize that this approach is close (identical?) to assuming some measurement error model ($\rho$) for the likelihood. But is that really worse than the approximation error from the thresholding approach?


title: Clustering method that is robust against inhomogeneous cluster sizes

text: I have data clustered in various sets whose amount of member variables and occupied area varies, the former among multiple magitudes. Here is an **example** : _Set A is very concentrated and has 100000 Elements. Set B is slightly more widespread and contains 1000 Elements. Set C is covers a comparably big area and has 10 Elements._ I am working on a one dimensional axis. Since I do not know the number of sets I want to obtain, I am looking for a clustering algorithm that returns a variable number of clusters. So far, hirarchical clustering has given good results when working with two different magnitudes of cluster sizes, but has failed me with three magnitudes, as in the example above.


title: What is the proper way to measure error in collected data for calibration purposes?

text: I got into a discussion with a colleague today about some of our data sensors. Our sensors weigh objects anywhere from 10 to 100 tonnes. currently we calculate the error as: $$ error = \frac{\left (x_{field} - x_{true} \right )}{x_{true}} $$ The problem is that if the error is consistent in terms of value, the percentages aren't. For example (fictitious data):               True_wt Field_wt    error       diff_wt     12000   13500       -0.125      -1500     34000   32500       0.0441      1500     65000   63000       0.0307      2000     24000   26000       -0.083      -2000     18000   19500       -0.083      -1500      My point was that estimating the error in such a way for our sensors is incorrect since the real issue is that they're off by `x` amount and not necessary off by a percentage. **Question:** In situations such as the example above, what would be a, or some of the, way(s) that the error in measurement could be measured in? **Edit 1: To reflect on @Whuber's comment** The context of the error is to solve the calibration of the sensor itself. The idea is to calibrate the field sensor to return a value that is very close to the true weight, so that over the year or so the sensor is in the ground, it stores as accurate and precise of weights as possible


title: 2x2x5 repeated measures ANOVA: significant 3-way interaction

text: I have 3 within-subject factors, namely `offset` (1px, ..., 5px), `side` (left, right) and `color` (red, green), which define the characteristics of the stimulus in a reaction time experiment. The DV is reaction time `RT`. The design is fully balanced. I ran a repeated measures ANOVA in R, like this:               options(contrasts = c("contr.helmert", "contr.poly"))          simon.aov <- aov(median.RT ~ color*side*offset + Error(VP / (color*side*offset)), data=dfa)      The results revealed a significant main effect of the `color`, as well as a significant interaction `color x side` and a significant 3-way interaction `color x side x offset`.   My primary focus lies on the interactions. **Specifically, I want to know on which of the 5 offsets (i.e. on which levels of the third factor) the 2-way interaction`color x side` reaches significance.** I am by no means familiar with post-hoc contrasts and multiple comparisons, but this question is the gist of the thesis that I'm working on. So my progress depends on an adequate test to examine this question. I highly appreciate any help on which test to run, and how to do this most efficiently in R. ## Edit: I'm sorry I didn't provide any plots earlier. @John: Here is the plot you requested. ![Plot 1 of 3-way anova results](http://i.stack.imgur.com/Z6alg.png) However, I believe, that this following plot rather clarifies my question: ![Plot 2 of 3-way anova results](http://i.stack.imgur.com/wCtKH.png) It seems like there is no `color x side` interaction at the first 3 levels of `offset`, but this interaction emerges at `offset` 4 and 5. This is what the plot seems to imply, however I don't know how to prove it statistically.


title: Why is it bad if number of dimensions / factors > sample size?

text: I've been told (read) this many times, but I never understood why it's bad for the number of dimensions in your data, or the number of explanatory variables in your model to be higher than your number of samples. Why is this the case?


title: Does the formula $X|(Y,Z) =(X|Y)|(Z|Y)$ hold?

text: Suppose we have three absolutely continuous random vectors $X, Y, Z$. Can anyone prove the following formula? $$X|(Y,Z) =(X|Y)|(Z|Y)$$ Put differently $p(X|Y,Z)=p(\hat{X}|\hat{Z})$, where $\hat{X}=X|Y$ and $\hat{Z}=Z|Y$. **Example of use** I want to determine distribution of $X|(Y,Z)$. Suppose, I know following \begin{align} \hat{X}=X|Y \sim N(\mu_x,\Sigma_x)\\ \hat{Z}=Z|Y \sim N(\mu_z,\Sigma_z) \end{align} So using the properies of normal distribution I get $$(\hat{X}|\hat{Z}=z) \sim N(\mu_x + \Sigma_{xz}\Sigma_z^{-1}(z-\mu_z),\Sigma_x-\Sigma_{xz}\Sigma_z^{-1}\Sigma_{zx})$$ with $Cov(\hat{X},\hat{Z})=\Sigma_{xz}$. Did I determine the distribution of $X|(Y,Z)$? **Update** I have found exactly the same type of derivation that I wanted to make on slides here: http://stanford.edu/class/ee363/lectures/kf.pdf - slides 14, 15 where this problem is stated as obvious.


title: Is there an analytical expression for the distribution of the max of a normal k sample?

text: For example:               k <- 100     R <- 10000     max.g <- numeric(R)          for(i in 1:R) max.g [i] <- max(rnorm(k))          hist(max.g)  # We can see it's right tailed...      I remember once encountering that there is a name for this type of distributions, but the name alludes me.


title: Fractal alternative to correlation

text: I am looking for a fractal-based statistical measure which could be used as alternative to correlation between two variables (I know that hurst exponent can be used for auto-correlation). Is anyone aware of such measures?


title: modeling changing items in market basket analysis

text: I'm performing a "dynamic" market basket analysis where a customer transaction is traced from the beginning (here it has an initial basket of product) to the finalized purchase (where there is the final basket) but he can change the items in the basket time by time and it is all traced. How can I model such dynamics?


title: When is Maximum Likelihood the same as Least Squares

text: In this paper on p315: http://www.ssc.upenn.edu/~fdiebold/papers/paper55/DRAfinal.pdf They explain that they use Levenberg Marquardt (LM) (along with BHHH) to maximize the likelihood. However as I understand it LM can only be used to solve Least Squares (LS) problems? Are the LS and MLE solutions the same for this type of problem? I know that when the errors are normal like in OLS then the solutions are the same. Here the processes being estimated are AR(1) so the errors are normal even though the overall process is not. Can I still treat the MLE and LS solution interchangeably in this situation? In which case can I just apply LM to solve the the LS solution safe in the knowledge that the optimal LS parameters are also the ones that will solve the MLE problem? Or does the LM algorithm have to be changed in some way so that it can be applied directly to the MLE estimation? If so how? Kind Regards Baz


title: Reverse scoring when question is stated in a negative fashion

text: In my survey, I had a question that asked respondents why they did NOT engage in a particular behavior (x). Some of the answers were stated in a positive way, such as "because of age" while some were stated in a negative way "I did NOT know how to do x". I am trying to reverse code because I am using these statements in a factor analysis and I would like the means that I compute to make sense (so that higher values indicate higher agreement with a positively stated statement). The problem is that when I run Cronbach's alpha test, the values I get make sense when I do NOT reverse code and they are negative when I DO reverse code. I am wondering if it is because my question was stated in a negative way to start with ("why did you not engage in x behavior?")? Myself and my supervisor feel that in this case it doesn't really make sense to reverse code, but I'm wondering how that will affect the means that I calculate? I will be using these means subsequently in a logistic regression.


title: Why does random-looking data give a really good model?

text: I ran `lm()` (I'm using R, but I don't think anything else in this question is R specific) then realized I'd used the wrong input data, basically meaningless. Yet, I got a really good model:               Residuals:           Min        1Q    Median        3Q       Max      -0.069654 -0.003899 -0.000381  0.004722  0.083622           Coefficients:                   Estimate Std. Error t value Pr(>|t|)         (Intercept)  1.795e-02  1.800e-03   9.970   <2e-16 ***     x           -3.676e-06  3.848e-07  -9.554   <2e-16 ***     ---          Residual standard error: 0.01416 on 11338 degrees of freedom     Multiple R-squared: 0.007986,   Adjusted R-squared: 0.007899      F-statistic: 91.28 on 1 and 11338 DF,  p-value: < 2.2e-16       I don't quite have permission to upload the data, and also it is large, but here is a plot of y and x (NB. x is on the vertical axis): ![My real data](http://i.stack.imgur.com/90O35.png) Okay, the break at roughly x==5200 is interesting, but basically it looks to me like for any value of x, y could be either side of y==0 with equal probability. So I don't see how it could possibly find a model with such a low p-value. To help understand, I tried to do the same with random-looking data:               d=data.frame(         y=(rnorm(10000,0,0.1))^2,         x=5000+runif(10000)*500         )     d$y=ifelse(runif(10000)<0.5,-d$y,d$y)      But `lm(y~x,data=d)` gives the terrible model I expected:               Residuals:           Min        1Q    Median        3Q       Max      -0.165179 -0.004537  0.000104  0.004742  0.134604           Coefficients:                   Estimate Std. Error t value Pr(>|t|)     (Intercept)  1.439e-03  6.324e-03   0.228    0.820     x           -2.967e-07  1.203e-06  -0.247    0.805          Residual standard error: 0.01737 on 9998 degrees of freedom     Multiple R-squared: 6.082e-06,  Adjusted R-squared: -9.394e-05      F-statistic: 0.06081 on 1 and 9998 DF,  p-value: 0.8052       That random data looks like this: ![random data](http://i.stack.imgur.com/Q4tK4.png) So, my real data and my random data behave completely differently when given to `lm()`, despite "x" looking equally useless as a predictor for "y" in both cases. It must be one of:   * I'm misreading the linear model output   * I'm asking the wrong question   * The "x" in my real data is genuinely a good predictor for "y"


title: Handling unbalanced data using SMOTE - NO BIG DIFFERENCE?

text: I have a classification problem with 2 classes. I have nearly 5000 samples, each of which is represented as vector with 570 features. The positive class samples are nearly 600. Meaning, I have a 1:8 ratio of positive and negative samples in the dataset. This imbalance in the dataset is mitigated using SMOTE. Subsequently, classification with 10 fold CV is performed. I get a f-measure of 0.91. To study the effect of imbalance in the dataset, I tried using the data with imbalance itself (i.e. without SMOTE). This time around, I observed a f-measure of 0.92. I understand the difference is using accuracy and f-measure to interpret the classifier predictions and since I have an unbalanced dataset, I chose to use f-measure. There seems to be no big difference in the end result whether or not I have an imbalance in the dataset in my case. In this context, I have the following questions:   * Why there seems to be no big difference in the f-measure in both the cases?   * It could be noted that, after I used SMOTE to mitigate the imbalance, the dataset becomes balanced and still I use f-measure to evaluate the classification results. Is it right to use f-measure in this case or should I use accuracy?   * SMOTE does oversampling of the minority class. Similarly, down sampling (or undersampling) the majority class could also rectify the imbalance. Why this methods is not preferred (If I may say so)? What effect does under sampling have on the classifier subjected to training and accuracy compared to oversampling.


title: Standardizing variables at item or index level: Does it make a difference?

text: I'm running some multiple group CFA models comparing covariance structure by race/ethnicity and have survey data from 6th, 8th, 10th and 12th graders. My supervisor has told me to combine 6th and 8th grade to run a middle school only model and to combine 10th and 12th grade for a high school only model. This makes good sense to me developmentally as a number of our variables mean different things at 6th grade as they do at 12th grade (attachment to parents, dating relationships etc). My concern is that I shouldn't just combine 6th and 8th grade without doing a within-grade standardization first. My primary question is at what level to standardize. A number of our variables of interest are composites of 2 to 6 items. Most of these indicies have been tested for reliability and validity by previous survey research. I could: a) standardize the items that make up the indicies within grade before combining them into indicies; OR b) standardize the indicies within grade and proceed from there. I just learning the statistical theory behind a lot of these things so I don't really have enough information to answer this question as of yet. Any thoughts would be much appreciated.


title: How to choose the training, cross-validation, and test set sizes for small sample-size data?

text: Assume I have a small sample size, e.g. N=100, and two classes. How should I choose the training, cross-validation, and test set sizes for machine learning? I would intuitively pick   * Training set size as 50   * Cross validation set size 25, and   * Test size as 25. But probably this makes more or less sense. How should I really decide these values? May I try different options (though I guess it is not so preferable... increased possibility of over learning)? What if I had more than two classes?


title: Proof for Sufficient and complete statistic (Shao)

text: Please can you help me, with this question: Let $X$ be a random variable with a distribution $P_\theta$ in $\{P_\theta : \theta \in \Theta\}$, $f_\theta$ be the p.d.f of $P_\theta$ w.r.t a measure $\nu$, $A$ be an event, and $\mathcal{P}_A=\{f_\theta I_A/P_{\theta} (A): \theta \in \Theta\}$.   * Show that if $T(X)$ is sufficient for $P_\theta \in \mathcal{P}$, then it is sufficient for $P_\theta \in \mathcal {P}_A$.   * Show that if $T$ is sufficient and complete for $P_\theta \in \mathcal{P}$, then it is complete for $P_\theta \in \mathcal{P}_A$.


title: Is there a kind of random walk process for wich the reciprocal of the variance is a sum of reciprocals of variances of compound processes?

text: **Context:** The variance of a sum of independent random walks is a sum of their variances: $\sigma^2 = \sigma_1^2 + \sigma_2^2$. In case of a dependent random walks with bivariate normal distribution it will be $\sigma^2 = \sigma_1^2 + \sigma_2^2 + 2 \rho \sigma_1 \sigma_2$. **Question:** Is there any known class of dependent random walk processes for which the reciprocal of variance is a sum of reciprocals of variances of each random walk process: $1/\sigma^2 = 1/\sigma_1^2 + 1/\sigma_2^2$ ? Can one infer any useful statistical information about such processes, given the relation holds true? **Remark:** The relation $1/\sigma_{total} = 1/\sigma_{diffusion-controlled- process} + 1/\sigma_{advection-controlled-process}$ follows from the coupling theory of eddy dispersion [1] (formula 2.10-8): the total dispersion of a chromatographic peak is due to diffusional-based and advectional-base mechanisms. This relation describes a corpus of the experimental data quite well, however Giddings derived the theory with a lot of assumptions and to me it seems more like an empirical correlation, which lacks rigorous statistical reasoning. I am trying to reverse engineer it and find out what kind of a relationship between two processes may lead to this peculiar rule for the sum of the variances. [1] Giddings, J. C. (1965). Dynamics of Chromatography, Part 1: Principles and Theory. New York, NY, USA: Marcel Dekker.


title: Can collaborative filtering be cast as a classical regression problem?

text: Having the Netflix challenge in mind: collaborative filtering is typically presented as a matrix dimension reduction. My question is how does the problem relate to classical regression (supervised learning) problems? Can it be seen as an instance of these or does it differ in some fundamental way?


title: Principal component Regression Using R

text: Can anyone explain principal component regression with the help of an example and the code in R? How to interpret the result of a principal component regression? How to find the individual effect of independent variables on dependent variable in case of PCA since we deal with components?


title: Noncentrality parameter in two-sample-t-test (Power analysis)

text: I am working with R, doing a post-hoc power analysis. Can anyone explain to me why these two commands generate different outputs? I thought they should not!               > 1-pt(tcrit-nc, df)     [1] 0.4248602     > 1-pt(tcrit, df, ncp=nc)     [1] 0.430632      I've got ...               n = 18     df = (n-1) * 2 = 34     tcrit = qt(0.05, df, lower = F)     nc = sqrt(n/2)*0.5 = 1.5 # 0.5 is the effect size      Edit: I found out, that R uses different algorithms to estimate the probabilities, depending on if the ncp-paramter is set or not. So can anybody confirm that this should be the only reason for the probabilities to differ?


title: Difference between MANOVA and two-way ANOVA

text: I have read the articles of both of these statistical methods on wikipedia, and as far as I can tell they are used in the same cases. I could not tell in how far they work the same way, since the two-way ANOVA article is rather elliptic. Could you tell me:   * Whether there is a difference between these two methods   * What the difference is   * What use cases is each method best suited for ? Many Thaks,


title: ARMA parameters

text: I am currently doing a short term forecast using ARIMA model.I have been following Box and Jenkins method and to choose the best ARIMA parameters to do do my forecast I tested various(p,q) combinations and took the one with the lowest aic(akaike information criteria). I assigned a d value as the number of differenciation required. My problem is that my forecast is converging into a straight line after 2 hours. Did I choose a wrong method? Is it because I have to try higer AR(p) and MA(q) parameters. I put a threshold for AR(p) at 4 and for MA(q) at 2. I saw in some articles, they were going upt to 24 to do there forecast. Is there a limit for choosing parameters? Can you give me some references please?


title: Domain analysis in surveys

text: When using the survey procs in SAS on imputed dated, is it best to include the imputation variable on the table or domain statement, or should it be included on the by statement? I found one paper (http://support.sas.com/resources/papers/proceedings10/265-2010.pdf) that included it on the domain statement which seems kind of strange to me since imputations are not really subsets of the data in the same way something like a grouping based on age or sex or race is. Can someone explain? Thx.


title: Clustering of 10's of millions of high dimensional data

text: I have a set of 50 million text snippets and I would like to create some clusters out of them. The dimensionality might be somewhere between 60k-100k. The average text snippet length would be 16 words. As you can imagine, the frequency matrix would be pretty sparse. I am looking for a software package / libray / sdk that would allow me to find those clusters. I had tried CLUTO in the past but this seems a very heavy task for CLUTO. From my research online I found that BIRCH is an algorithm that can handle such problems, but, unfortunately, I couldn't find any BIRCH implementation software online (I only found a couple of ad-hoc implementations, like assignment projects, that lacked any sort of documentation whatsoever). Any suggestions? **EDIT** : Here is some explanation on what "dimensionality" and "similarity" is in the context of my problem: Each text-snippet is a sentence taken from a huge text corpus (news articles). A vector is created for each sentence by projecting each text snippet onto a N-dimensional space where each **dimension** corresponds to a word. So, if a sentence is:               "I like red apples more than green apples"       then the corresponding vector would have dimensions "I":1, "like":1, "apples":2 "more":1 etc... _I will decide which of the millions of distinct words will make up the space upon which the projection will be made. Typically the top-N most frequent words will be selected for this._ Then, the **similarity** will simply be the distance between the vectors of two sentences. (usually I normalize the vectors first and then take the dot product between two vectors - this is known as "cosine similarity")


title: P value for chi-squared goodness of fit of geometric distribution doesn't make sense from R fitdist

text: Trying to teach myself R and found the following conundrum. I wrote a script to simulate rolling a die and adding the value until a "1" is rolled. This worked well (eventually). The mean and variance matched with the theoretical (1-p)/p. I took one of the runs of 1000 samples and tried to fit a geometric distribution using the R function fitdist from package fitdistrplus. The fitdist seemed to fit reasonably with a "prob" estimate of .139 and std error of .004. Graphs look OK. I then ran gofstat on the output, and the Chi-squared p-value was 10^-11? By the definition of how the data were created, it is a geometric distribution, so why did Chi square fail? Or do I have a fundamental misunderstanding of the underlying maths? I'm not using scripts, but just running commands from the console. Thanks. edit: code added   Note: I started calculating the values, rather than rolls, then realized it came out to 4 times the rolls (4 being average roll on a d6). Doh. It was working so I left it in.               samples <- 1000     runs <- 10     vecAttacks <- 1:samples     vecRolls <- 1:samples     superAve <- 1:runs     superRoll <- 1:runs     for (j in 1:runs){     for (i in 1:samples){       attack <- 0       numrolls <- 0 ##numrolls is the number of successful rolls       repeat{         roll <- sample(6,1,TRUE)         if(roll == 1){           vecAttacks[i] <- attack           vecRolls[i] <- numrolls           break         } else {           attack <- attack + roll           numrolls <- numrolls + 1         }       }     }     AvgAttack <- mean(vecAttacks)     AvgRolls <- mean(vecRolls)     superAve[j] <- AvgAttack     superRoll[j] <- AvgRolls     }     print("sup=")     print (mean(superAve))     print (mean(superRoll))          > source('~/.active-rstudio-document')     [1] "sup="     [1] 20.3182     [1] 5.0747          > swap <- fitdist(vecRolls, distr = "geom", method = "mle")     > swap     Fitting of the distribution ' geom ' by maximum likelihood      Parameters:           estimate  Std. Error     prob 0.1714678 0.004935428     > gofstat(swap)     Chi-squared statistic:  13.06858      Degree of freedom of the Chi-squared distribution:  10      Chi-squared p-value:  0.2198627      Chi-squared table:           obscounts theocounts     <= 0  172.00000  171.46777     <= 1  141.00000  142.06657     <= 2  120.00000  117.70673     <= 3   86.00000   97.52382     <= 4   99.00000   80.80163     <= 5   63.00000   66.94676     <= 6   65.00000   55.46755     <= 8   64.00000   84.03321     <= 10  61.00000   57.68592     <= 12  42.00000   39.59940     <= 15  41.00000   37.38888     > 15   46.00000   49.31177          Goodness-of-fit criteria                                1-mle-geom     Aikake's Information Criterion   5346.514     Bayesian Information Criterion   5351.422


title: Python alternative to Factorie

text: Is there alternative / replacement in Python of Factorie (http://factorie.cs.umass.edu/index.html) - In particular I am looking for a tool in Python that allows me to create arbitrarily structured CRFs for a custom NER task. I have found https://pystruct.github.io/ but it looks that does not support arbitrarly factor graphs yet.


title: What theory can explain why after several rounds of bootstrapping, the result come to converge?

text: When I use bootstrapping (updating the negative samples generated by existed classifier model and then retraining classifier model for next round; And do this for several rounds) in machine learning, I find that after several round (maybe 10 rounds), the result come to converge and stable even if the classifier model still can generate some negative sample. What theory can explain this theory that bootstrapping will come to converge and the termination existed? Is there a margin of it? Thanks.


title: A statistics book that explains using more images than equations

text: I have become interested in statistics, but I must admit that it has been a long time since I have used math seriously. Sometimes I understand what equations mean, but sometimes I can't follow them. I like the answer given here that uses the image with an arrow: Conceptual understanding of root mean squared error and mean bias deviation. Would you have any advice regarding serious statistical books on modeling, PCA, simulation, etc., that use figures more than equations? I'm hoping for something like a 'stats for dummies' book, a book that would make stats accessible for everyone?.


title: Which machine learning blog do you recommend?

text: In the spirit of What statistical blogs would you recommend? I propose a collection of machine learning blogs worth reading.


title: Collinear variables in Multiclass LDA training

text: I'm training a Multi-class LDA classifier with **8 classes** of data. While performing training, I get a warning of: " **Variables are collinear** " I'm getting a **training** accuracy of over **90%**. I'm using **scikits-learn** library in **Python** do train and test the Multi- class data. I get decent **testing** accuracy too(about **85%-95%** ). I don't understand what the error/warning means. Please help me out.


title: Multivariate time series model evaluation with conditional moments

text: Consider multivariate time series models that estimate potentially time- varying conditional means, variances, and correlations (one type of model might be a VAR(p)+Garch(1,1)+DCC Gaussian Copula model). I will simulate from one "true" distribution and compare other potential models to it. I don't particularly care about RMSE because I care just as much about the conditional covariance matrix as the conditional means. For the sake of simplicity, assume that the conditional mean and conditional covariance are sufficient to explain "true" distribution. I can calculate the conditional mean vector and conditional covariance matrix of each model and the true model at any given point in time. Alternately, I can calculate the forecasted conditional means and covariances on a rolling basis.   1. Is there any reason to prefer the Kullback-Leibler divergence or Bhattacharyya distance as comparison methodologies or should I just do both? Alternately, is there a different method that is even more appropriate.    2. In the constant multivariate Gaussian case, is it possible to make statistical inferences based on these comparison methodologies, such as that the divergence is not statistically bigger than zero or that the divergence of one is bigger than another?    3. If I calculate the conditional means and variances at any given point in time, I could calculate the comparison statistic in each period. This would produce a time series of divergences/distances. Is it reasonable to make statistical inferences in this case by just assuming it is a typical time series?


title: What is the difference between a partial likelihood, profile likelihood and marginal likelihood?

text: I see these terms being used and I keep getting them mixed up. Is there a simple explanation of the differences between them?


title: Choosing the number of clusters in hierarchical agglomerative clustering

text: I have a set of points that I want to cluster into groups according to a number of features computed. I have distance matrix containing the distances between all different pairs of points. I have tried K-Means, and DBSCAN first but since I have no clue about that number of clusters underlying in the data, required for K-Means or the optimal Epsilon and MinPnts required for DBSCAN. I decided to turn to Hierarchical Agglomeration Clustering since it doesn't require setting any parameters for clustering. I used R to perform HAC. Now to turn the resultant dendrogram into a number of groups of points (flat clusters), I want to choose which level to cut the tree at (the same as choosing the number of clusters). I managed to do this for different values for the number of clusters and evaluate the Silhouette value for each clustering results and therefore I can choose the number of clusters that gives me the maximum silhouette value as the optimal number of clusters. What happens is that I find that the Maximum silhouette value (0.8) obtained is for a number of clusters = 5 but the cluster sizes are not very good (one cluster is > 900 points, second one is 5 points, and the other three are single point for each). I was using group average for the HAC cluster-cluster distance.. I tried also Ward's method (which is supposed to give clusters of equal sizes) but although this one give me somehow better distribution for the number of points in each cluster... The silhouette value for the clusters produced is really low (0.1). Do you suggest any good methods for doing clustering without requiring parameters? or to obtain the optimal values for clustering parameters (dendrogram cutting level for HAC, (epsilon, MinPnts) for DBSCAN)?


title: Intercept significant but not the variables in GLM

text: I got these results from a GLM from R, but I am not sure what it means that the intercept is significant but none of my variables actually are. So is nothing really significant? Here are my results:               Groups Name        Variance   Std.Dev.  Corr         Plot   (Intercept) 1.4842e+07 3852.4696                     Year        3.6783e+00    1.9179 -1.000      Number of obs: 121, groups: Plot, 50          Fixed effects:                        Estimate Std. Error z value Pr(>|z|)        (Intercept)        -3.90239    1.36721  -2.854  0.00431 **     Wald_ManagBUBAdi  -16.50851 9563.66009  -0.002  0.99862        Wald_ManagBUBAj     0.45015    2.31227   0.195  0.84564        Wald_ManagBUBAun    1.63915    1.46892   1.116  0.26447        Wald_ManagBUdi     -0.03527    2.25725  -0.016  0.98753        Wald_ManagBUplen    1.14402    1.52959   0.748  0.45450        Wald_ManagFIBA      2.92878    1.53228   1.911  0.05595 .


title: Solving for regression parameters in closed-form vs gradient descent

text: In Andrew Ng's machine learning course, he introduces linear regression and logistic regression, and shows how to fit the model parameters using gradient descent and Newton's method. I know gradient descent can be useful in some applications of machine learning (e.g., backpropogation), but in the more general case is there any reason why you wouldn't solve for the parameters in closed form-- i.e., by taking the derivative of the cost function and solving via Calculus? What is the advantage of using an iterative algorithm like gradient descent over a closed-form solution in general, when one is available?


title: Forecasting handbooks

text: In engineering, we usually have Handbooks that pretty much dictate the state of the practice. These books are usually devoid of theory and focus on the applied methodology. Is there a forecasting Handbook out there? that solely focuses on the technique and not the background?


title: Minimum variance unbiased estimator

text: **A**. Determine the restrictions needed on the $a$'s so that $z=\sum_{i=1}^{n}a_iX_i$ will be an unbiased estimator $E\left[X\right]$, where $X_1,X_2,\dotsc,X_n$ represents a random sample of $X$. **Answer:** The $a$'s must sum to one. **B**. What is the best set of $a$'s to choose from **A** if $z=\sum_{i=1}^{n}a_iX_i$ is to be an unbiased estimator of $E\left[X\right]$ with _minimum variance_. **Answer:** $a_i=1/n$. The first answer I've shown I understand. However I'm a bit stumped as to how to arrive at the answer to part **B**. **EDIT:** I'm not exactly sure what the question is asking. What I think it's asking is: find $$\left\{a_i:\mathrm{var}\left(\sum_{i=1}^{n}a_iX_i\right)=\frac{1}{n\sigma^2}\right\},$$ where $\sigma^2$ is the expected value of the squared derivative of some arbitrary log likelihood function of $X$ given some parameter $\theta$. I know that $$\mathrm{var}\left(X\right)=E\left[\left(X-\theta\right)^2\right].$$ How do I compute the variance of $z$? $\mathrm{var}\left(z\right) = E\left[z^2\right] - \theta^2$. How do I evaluate $E\left[z^2\right]$ without a density function? Thanks!


title: How to estimate a model with fixed and random effects for a long panel dataset?

text: NOTE: I am using Stata for doing this. I have a long panel dataset, meaning my N is much smaller than my T. I have N = 5, T = 61. I tried to estimate my model, but I get an error related to the fact that I do not have enough degrees of freedom for calculating as many variables as I have. This is what I get:               Wald chi2(4)       =         .     Prob > chi2        =         .      Stata's help tells me the following: > The VCE you have just estimated is not of sufficient rank to perform the > model test. As discussed in [R] test, the model test with clustered or > survey data is distributed as `F(k,d-k+1) or chi2(k)`, where k is the number > of constraints and d=number of clusters or d=number of PSUs minus the number > of strata. Because the rank of the VCE is at most d and the model test > reserves 1 degree of freedom for the constant, at most d-1 constraints can > be tested, so k must be less than d. The model that you just fit does not > meet this requirement. > > To simplify the remaining discussion, let's consider the case of clustered > data. This discussion applies to survey estimation in general by > substituting, "PSUs - strata" for "clusters". > > There is no mechanical problem with your model, but you need to consider > carefully whether any of the reported standard errors mean anything. The > theory that justifies the standard error calculation is asymptotic in the > number of clusters, and we have just established that you are estimating at > least as many parameters as you have clusters. > > That concern aside, the model test statistic issue is that you cannot > simultaneously test that all coefficients are zero because there is not > enough information. You could test a subset, but not all, and so Stata > refuses to report the overall model test statistic. > > Here note the degrees of freedom reported for the chi2 or F. You might see > chi2(6) or F(6, 5). If you were to count the number of coefficients that > would be constrained to 0 in a model test in this case, you would find that > number to be greater than 6. You could find out what that number is by re- > estimating the model parameters without the `vce(robust)` and `vce(cluster > clustvar)` options. In any case, the 6 reported is the maximum number of > coefficients that could be simultaneously tested. Reading the book "Microeconometrics using Stata", the author clearly indicates that the `xtreg` command, with the `vce(cluster id)` option for calculating robust errors, is mostly appropriate for short panels. An alternative is to use the command `xtregar` for estimating random and fixed effects, which is based on an AR(1) process for the errors. However, we have tested, and our errors do not show an AR(1) pattern (`xtserial` test). The `xtregar` command has the option `rhof(#)`, where # indicates the desired rho value (AR(rho)). So the questions are two: 1) What is the right way to calculate Fixed and Random Effects for a long panel? 2) Would specifying rho = 0 completely eliminate the AR(1) process for the errors, and leave us with an estimation that fits our data?


title: continuous response variable with lots of zero values

text: I'm trying to determine the effects of four treatments on flowers on the volume of nectar available to pollinators. Problem is, I have about 50% zero volumes across treatment categories. What is the best approach to comparing mean volumes between treatments?


title: Top R package for PLS regression?

text: I'm very new to **R** and PLS-regression. I would like to know, based on your experience, which **R** packages for PLS-regression are most highly recommended. My area of application is chemistry.


title: R: Package 'sphet' on spatial autoregressive models

text: I'm using the function `gstslshet` from the package `sphet` in R, but I'm wondering what the outcome is. Is it a model which can be used to make new predictions? I've put my summary output below (I am new to this kind of modelling).                Generalized stsls          Call:     gstslshet(formula = price_value ~ numberofbedrooms + energyconsumption,          data = regression, listw = listw)          Residuals:        Min. 1st Qu.  Median    Mean 3rd Qu.    Max.      -422000  -91100  -29500     -42   69800  701000           Coefficients:                          Estimate  Std. Error t-value  Pr(>|t|)         (Intercept)       -6.9897e+04  2.3622e+04 -2.9589  0.003087 **      numberofbedrooms   4.0633e+04  2.8040e+03 14.4911 < 2.2e-16 ***     energyconsumption -5.0070e+01  2.5359e+01 -1.9744  0.048332 *       lambda             8.2850e-01  4.8803e-02 16.9763 < 2.2e-16 ***     rho               -6.1301e-01  8.4064e-02 -7.2922 3.049e-13 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Wald test that rho and lambda are both zero:      Statistics: 9.3902 p-val: 0.0021815


title: Estimation of log-likelihood via importance sampling

text: I am looking at a model trained with stochastic gradient variational Bayes. In this paper an importance sampler is proposed to estimate the likelihood: $$p(x) \approx {1 \over S} \sum_{s=1}^S {p(x|z_s) p(z_s) \over q(z_s|x)}$$ with $h_s \sim q(z_s|x)$. It is not entirely necessary what the quantities mean in order to understand the question, but I will give some info. We are looking at a graphical model of the form $z \rightarrow x$, with $z$ being continuous. The map from $z$ to $x$ is a nonlinear regressor representing a conditional distribution $p(x|z)$. $p(z)$ on the other hand is just a prior. $q(z|x)$ is a "recognition model", which approximates the intractable $p(z|x)$ with another nonlinear regressor. Consider that these distributions are given and fixed, their estimation is a different matter--it is a variational method which also gives us an upper bound on the negative log-likelihood. Now, it seems I have reproduced the results on their models and one data set. Yet, I have two practical questions. 1) How to test an importance sampler? Since $p(x)$ is intractable, I guess it can only be done for small models. The route I went is to assert its correctness on factor analysis, where $p(x)$ is tractable and I made up a very simple $q(x|h)$. 2) The other question is that I am ultimately interested in the negative log- likelihood for comparative reasons. Thus I arrived at $$-\log p(x) = \log S - \log (\sum_{s=1}^S \exp(\log p(x|z) + \log p(z) - \log q(z|x)).$$ A nice side effect is that this quantity can be computed numerically stable with the log-sum-exp trick. However, depending on the sample size $S$, I get the following results for my negative log likelihood:               1     -> 29.6682416864     20    -> 22.144055709761528     50    -> 21.067476078084795     100   -> 20.458267754505115     200   -> 19.901382220921288     1000  -> 18.911037074948325     10000 -> 17.8730836533      ![Results of importance sample given sample size.](http://i.stack.imgur.com/Xxr1G.png) I do not understand this asymptotic behaviour. Does it make sense? Or is sth wrong in the deviation of the log-likelihood sampling. Furthermore, I have an upper bound on that value, which is ~14. Thus, if everything is correct this should converge to a value lower than that.


title: How can I prevent Overflow for the Gamma function?

text: I have a Dirichlet distribution from which I'm maximizing $\alpha$ by calculating the log likelihood with the following equation $p\left(s|\alpha\right) = \sum\limits_{i=1}^O\log\left(\frac{n_i!}{\prod\limits_{j=1}^kn_{ij}!}\frac{\Gamma\left(\alpha_0\right)}{\Gamma\left(\sum\limits_{j=1}^k\left(n_{ij}+a_j\right)\right)} \prod\limits_{j=1}^k \frac{\Gamma\left(n_{ij}+\alpha_j\right)}{\Gamma\left(\alpha_j\right)}\right)$ I have a python script to calculate the log likelihood. It works for small values, but when $n_{ij}$ gets high $\Gamma\left(\sum\limits_{j=1}^k\left(n_{ij}+a_j\right)\right)$ and $\Gamma\left(n_{ij}+\alpha_j\right)$ python runs into overflow error (the max for math.gamma is between 171 and 172). Is there a way to change the likelihood function so that $\Gamma\left(\sum\limits_{j=1}^k\left(n_{ij}+a_j\right)\right)$ and $\Gamma\left(n_{ij}+\alpha_j\right)$ don't run into overflow errors. Can I just divide the values in $\Gamma$ by for example 100, if I then compare the likelihood between $p(s|\alpha)$ with different $\alpha$ values do the best $\alpha$ values still have the highest likelihood?


title: Distribution of random vector divided by it's norm

text: Does anyone know anything about the distribution of a random vector divided by it's norm? More specifically, if we know the distribution of a vector valued random variable $x$, what can we say about the distributions of $\frac{\vec{x}}{||\vec{x}||_2}$ and $\frac{\vec{x}}{||\vec{x}||^2_2}$ 1) In general 2) If $\vec{x} \sim Normal_p(\vec{\mu}, \Sigma)$


title: How to understand missing data mechanisms using DAGs

text: I'm wondering when it is possible to distinguish between a data generating process that is missing at random (MAR) vs one that is not missing at random (NMAR) by analyzing a directed acyclic graph (DAG). For example, suppose we have the following: $$ \begin{matrix} T& \longrightarrow & Y \\ \downarrow & & \\ R_{Y} \end{matrix} $$ where $R$ is the indicator for missingness on the outcome variable $Y$ and $T$ is treatment and always observed. In this case it's clear that $R$ is conditionally independent of $Y$ given $T$, which is what's required for making a determination of MAR (at least according to most definitions I've seen). But suppose we had a slightly different graph: $$ \begin{matrix} T& \longrightarrow & Y \\ \downarrow & & \\ L & & \\ \downarrow \\ R_{Y} \end{matrix} $$ where $L$ is an unobserved/latent variable. Suppose that $R_Y=0$ ($Y$ is missing) if $L$ is greater than some threshold value. Like in the first graph, we still have a situation in which missingness on $Y$ is conditionally independent of $Y$ given $T$ (which is always observed). But now missingness on $Y$ does depend on the unobserved value of $L$.   * Does that make this NMAR?    * Or is it still MAR since even though we don't observe $L$ we know that it depends on the always observed variable $T$?


title: Optimal orthogonal polynomial chaos basis functions for log-normally distributed random variables

text: I hope this is the appropriate venue for this type of question. If not, please feel free to migrate! :) I'm trying to solve a stochastic partial differential equation of the form $$\alpha(\omega)\nabla^2u=f$$ where $\alpha(\omega)$ represents a random field that is log-normally distributed, i.e. it has a probability density function $$pdf(x)=\frac{1}{x\sqrt{2\pi\sigma}}e^{-\frac{(ln(x)-\mu)^2}{2\sigma^2}}.$$ I want to represent the solution of this problem as a polynomial chaos expansion $$u=\sum_{i=0}^p u_i(x)\Psi_i(\xi)$$ where $u_i(x)$ is a deterministic coefficient and $\Psi_i(\xi)$ are orthogonal polynomials in terms of a random variable $\xi$ with the same log-normal probability density function. According to Xiu & Karniadakis (2002), certain orthogonal polynomial bases give optimal (exponential) convergence of finite expansions to the true solution $u$. For instance, hermite polynomials are optimal for gaussian distributions, legendre polynomials for uniform distributions, laguerre for gamma distributions etc (see the above paper, bottom of page 8). What is the corresponding optimal polynomial basis for log-normal distributions?


title: What's the difference between observable fixed effect and control variable?

text: I am confused about the exact definitions here. Assuming I have a cross-sectional regression, let's say, Wage on Education and I additionally control for observable characteristics with a set of dummies or variables like intelligence level, age, parent's education level, urbanization area, gender, race, work experience etc.   1. Does this mean I used dummies/variables to " **control** for observable fixed effects" that I obtained through my data collection? (Is, for example, parent's education level thus an **observable** fixed effect?)   2. Do **unobservable** fixed effects like ability (as often quoted in the literature) are then said to be "controlled" by proxies through my dummies like intelligence level, experience etc.? Sorry for the basic question here. In a way I'd like to know the exact difference between controlling for a variable, observable and unobservable characteristics, observable and unobservable fixed effects. Thanks.


title: Explain Statistics: Matching formulas for chi square tests

text: The chi square formula given in my book is: $\chi^2 = \frac{(n-1)s^2}{\sigma^2} $ At first I admitted to feeling some kind of strangeness to this formula, but after a few hours I realized that $\chi^2$ was $\displaystyle \sum_{i=1}^n \left(\frac{X_i - \mu}{\sigma}\right)^2 - \left( \frac{\bar{x} - \mu}{\frac{\sigma}{\sqrt{n}}} \right)^2$ in disguise. However, there's a formula in my book that says that $\chi^2 = \displaystyle \frac{\sum_{i=1}^n (O_i - E_i)^2}{E_i}$ where $O_i$ is the observed value and $E_i$ is the expected value. I want to show that this expression is equivalent to the one above: $$\displaystyle \sum_{i=1}^n \left(\frac{X_i - \mu}{\sigma}\right)^2 - \left( \frac{\bar{x} - \mu}{\frac{\sigma}{\sqrt{n}}} \right)^2$$ This equation comes from the fact that $(n-1)s^2 = \displaystyle \sum_{i=1}^n ((X_i - \mu)+(\mu - \bar{x}))^2$


title: Difference between White's reality check and Monte Carlo permutation

text: In finance, two procedures are sometimes used to help evaluate the likelihood that positive results were obtained strictly due to data mining:   * White's Reality Check (WRC), described in White's 2000 Econometrica paper   * Monte Carlo Permutation (pair matching of historical data with signals without replacement) These procedures apply when a researcher compares, for example, numerous trading rules, and selects the one that performs best on the historical data. The question is whether the best performing rule is much better than what one would expect from the best rule selected from among useless, uninformative rules. I read the following passage in a textbook explaining the difference between WRC and MCP: > Both WRC and MC test the null hypothesis that all rules in a data-mined > universe are useless. In the context of WRC, a useless rule has an expected > return equal to zero. In the context of MCP, useless refers to a rule whose > output values are randomly paired with the market's future change. Aren't these two null hypotheses equivalent?


title: Proportional odds assumption in ordinal logistic regression in R with the packages VGAM and rms

text: An assumption of the ordinal logistic regression is the proportional odds assumption. Using R and the 2 packages mentioned I have 2 ways to check that but I have questions in each one. 1) Using the rms package Given the next commands               library(rms)     ddist <- datadist(Ki67,Cyclin_E)     options(datadist='ddist')     f <- lrm(grade ~Ki67+Cyclin_E);f     sf <- function(y)     c('Y>=1'=qlogis(mean(y >= 1)),'Y>=2'=qlogis(mean(y >= 2)),'Y>=3'=qlogis(mean(y >= 3)))     s <- summary(grade ~Ki67+Cyclin_E, fun=sf)     plot(s,which=1:3,pch=1:3,xlab='logit',main='',xlim=c(-2.5,2.5))      I have               lrm(formula = grade ~ Ki67 + Cyclin_E)          Frequencies of Missing Values Due to Each Variable        grade     Ki67 Cyclin_E             0        0        3                                     Model Likelihood     Discrimination    Rank Discrim.                                 Ratio Test            Indexes          Indexes                 Obs            42    LR chi2     11.38    R2       0.268    C       0.728          1             11    d.f.            2    g        1.279    Dxy     0.456          2             15    Pr(> chi2) 0.0034    gr       3.592    gamma   0.458          3             16                         gp       0.192    tau-a   0.308         max |deriv| 1e-07                         Brier    0.166                                             Coef    S.E.   Wald Z Pr(>|Z|)     y>=2     -0.1895 0.8427 -0.22  0.8221       y>=3     -2.0690 0.9109 -2.27  0.0231       Ki67      0.0971 0.0330  2.94  0.0033       Cyclin_E -0.0076 0.0227 -0.33  0.7387       The `s` table gives: (unfortunately I don't know how to upload a graph made in R)               grade    N=45          +--------+-------+--+----+---------+----------+     |        |       |N |Y>=1|Y>=2     |Y>=3      |     +--------+-------+--+----+---------+----------+     |Ki67    |[ 2, 9)|12|Inf |0.6931472|-1.0986123|     |        |[ 9,16)|12|Inf |0.3364722|-2.3978953|     |        |[16,24)|10|Inf |2.1972246| 0.0000000|     |        |[24,44]|11|Inf |2.3025851| 1.5040774|     +--------+-------+--+----+---------+----------+     |Cyclin_E|[ 3,16)|15|Inf |1.0116009|-0.1335314|     |        |[16,22)| 7|Inf |1.7917595|-0.9162907|     |        |[22,33)|10|Inf |1.3862944|-0.8472979|     |        |[33,80]|10|Inf |0.4054651|-0.4054651|     |        |Missing| 3|Inf |      Inf| 0.6931472|     +--------+-------+--+----+---------+----------+     |Overall |       |45|Inf |1.1284653|-0.4054651|     +--------+-------+--+----+---------+----------+      Where for the Ki67 I see that 3 out of the 4 differences `logit(P[Y> = 2])-logit(P[Y> = 3])` are close to 2. Only the last one is quite lower (around 0.8). But here Ki67 is continuous and not categorical so I don't know if the results of the table are correct and there isn't any p-value to decide. By the way I run the above in SPSS and I didn't reject the assumption. 2) Using the VGAM package Here using the next commands I have the model under the assumption of proportional odds               library(VGAM)     fit1 <- vglm(grade ~Ki67+Cyclin_E,family=cumulative(parallel=T))     summary(fit1)      And the results               Coefficients:                     Estimate Std. Error  z value     (Intercept):1  0.1894723   0.820442  0.23094     (Intercept):2  2.0690395   0.886732  2.33333     Ki67          -0.0970972   0.032423 -2.99467     Cyclin_E       0.0075887   0.021521  0.35261          Number of linear predictors:  2           Names of linear predictors: logit(P[Y< = 1]), logit(P[Y< = 2])          Dispersion Parameter for cumulative family:   1          Residual deviance: 79.86801 on 80 degrees of freedom          Log-likelihood: -39.93401 on 80 degrees of freedom          Number of iterations: 5       While using the next commands I have the model without the assumption of proportional odds               fit2 <- vglm(grade ~Ki67+Cyclin_E,family=cumulative(parallel=F))      where unfortunately i receice the next message > Warning message: In vglm.fitter(x = x, y = y, w = w, offset = offset, Xm2 = > Xm2, : convergence not obtained in 30 iterations However if I type `summary(fit2)` I get results but again I don't know if they are correct. My intention was to use the next commands and get the answer but know I doubt if this is correct (by the way if I do it I get `p-value=0.6`.               pchisq(deviance(fit1)-deviance(fit2),     df=df.residual(fit1)-df.residual(fit2),lower.tail=FALSE)      So, regarding the methods mentioned above does anyone knows whether the results I get are valid or in the case of the VGAM package is there any way to increase the number of itterations?Is there any other way to check it?


title: approximation to maximum and minimum function : soft-min and soft-max

text: The approximation to the function max(x) can be written as a NOISY-OR as given below: $ max_k(x) = 1-\prod_k(1-x) $ Are there any way to approximate min(x) ?


title: Evaluation of a GLS fit

text: My first post in this useful community. Ive been trying to use GLS methods in non-linear mixed effects models and have found this post, Prediction with GLS extremely informative and helpful. Some of my questions are probably more of a follow up to the above post which discusses the predictive capability of GLS. Below I will state my question and give a brief description of why I would like to know more.   * How does one evaluate a GLS fit for a **real dataset** as compared to a regular fit? Using simulations I could show that GLS is robust to error model mis- specification, especially when the data is sparse or there is deviation from normality (heavy-tailed distribution or time-varying distribution of residuals). For a simulated model, there is decreased bias and improved precision of the model parameters in a gls fit compared to the regular nlme fit. However, I don't understand how would I evaluate my gls fit for a real dataset as compared to a regular fit. Standard errors/Goodness of Fit/etc..?   * At which stage of model building does one test a GLS fit? Carrol and Ruppert state that when there is enough evidence of heteroschedasticity of your variance, then GLS estimators are very efficient. I obviously could see that in simulations studies, but how do I justify a better fit in a **real dataset** when the parameter estimates of the model from a gls fit are not very different from those obtained after a regular fit (even though the gls estimator is more efficient). This obviously comes back to the question, why (or actually why not) use GLS estimation at all.   * The estimation of residual variance parameters is independent of response prediction. T/F? If the variance changes with your predictions (e.g. a constant CV model), is my understanding of the above statement correct that the estimation of the residual variance parameters is independent of model prediction of responses (or vice-versa). I believe this is the basic principle behind GLS.   * Can anyone suggest good references whether gls is design dependent ( both structure and density of observations)   * Any comments on gls performance/robustness for variable selection. Given bullet 3 above, can we expect gls to be more sensitive in picking up predictors for a response variable. Thanks in advance, VJ


title: Logistic regression versus t-test, which is better for this comparative-descriptive study?

text: Which type of test is most appropriate for this comparative-descriptive study? I have two groups of subjects (250 in each). My goal is to determine if the two subject groups are different and if so, how. I'm using an online survey of 50 statements (embedded are 5 statements about 10 different possible factors). Subjects rate these statements on a Likert scale from "strongly agree" to "strongly disagree". Do I have to use a t-test or a logistic regression to determine categorical differences, and predictors?


title: Joint PMF for two Geometric distribution variables

text: I am interested to know how to calculate the joint probability mass function for two independent geometric random variables. Suppose two variables X1 and X2 are independent, such that `Xi∼Geometric(theta)`, how to find the joint pmf distribution of X1 and X2. I am not sure but I think it should be the product of pmf of both mass function. Also, how should I calculate the probability of the event where kth trial being the the first success/failure for both the variables or k1th trial for X1 and k2th trial for X2 ? Thanks.


title: Significant difference real or due to the internal variability?

text: In my data I have 9 different sets of data for 2 different groups. Each one of these datasets is the same measurement changing over the time. If I make a graph, I can see 9 lines for each group. I did an average of the 9 lines for each group, and I can get 2 lines with SD and N, one for each group. How do I know if one of the 9 lines I used is just an outlier? I want to show that the variability is not making appear a difference between these 2 groups but is a real difference. Maybe an ANOVA, but I'm not sure.


title: Logistic regression residual analysis

text: This question is sort of general and long-winded, but please bear with me. In my application, I have many datasets, each consisting of ~20,000 datapoints with ~50 features and a single dependent binary variable. I am attempting to model the datasets using regularized logistic regression (R package glmnet) As part of my analysis, I've created residual plots as follows. For each feature, I sort the datapoints according to the value of that feature, divide the datapoints into 100 buckets, and then compute the average output value and the average prediction value within each bucket. I plot these differences. Here is an example residual plot: ![](http://i.imgur.com/GXKYq.png) In the above plot, the feature has a range of [0,1] (with a heavy concentration at 1). As you can see, when the feature value is low, the model appears to be biased towards overestimating the likelihood of a 1-output. For example, in the leftmost bucket, the model overestimates the probability by about 9%. Armed with this information, I'd like to alter the feature definition in a straightforward manner to roughly correct for this bias. Alterations like replacing $x \rightarrow \sqrt{x}$ or $x \rightarrow f_a(x) = \cases{a & if $x<a$ \cr x & else}$ How can I do this? I'm looking for a general methodology so that a human could quickly scroll through all ~50 plots and make alterations, and do this for all datasets and repeat often to keep models up-to-date as the data evolves over time. As a general question, is this even the right approach? Google searches for "logistic regression residual analysis" don't return many results with good practical advice. They seem to be fixated on answering the question, "Is this model a good fit?" and offer various tests like Hosmer-Lemeshow to answer. But I don't care about whether my model is good, I want to know how to make it better!


title: How are Taguchi orthogonal arrays derived?

text: I am implementing statistical analysis in a C# program and am using the Taguchi Method to do it (the Taguchi arrays are about a 1/4th of the way down the page). The article in the above link mentions that you can look up smaller Taguchi arrays and derive others, but I can't seem to find any resource that explains how to derive these Taguchi arrays. I can hard-code the smaller arrays in my program, but I want to be able to derive the others. Is there a resource or link that explains, as simply as possible, how to derive Taguchi arrays? Maybe would someone here be able to explain how to derive them?


title: Solving simultaneous linear equations with R

text: Am trying to solve a system of linear equations involving three variables using `R`. However, I get an error as shown below:               > A = matrix(data=c(-0.4, 0.4, 0.15, 0.3, -0.55, 0.6, 0.1, 0.15, -0.75), nrow=3, ncol=3, byrow=TRUE)     > A          [,1]  [,2]  [,3]     [1,] -0.4  0.40  0.15     [2,]  0.3 -0.55  0.60     [3,]  0.1  0.15 -0.75     >      > b = matrix(data=c(0, 0, 0), nrow=3, ncol=1, byrow=TRUE)     > b          [,1]     [1,]    0     [2,]    0     [3,]    0     >      > solve(A, b)     Error in solve.default(A, b) :        system is computationally singular: reciprocal condition number = 1.04615e-17      When I work them out using substitution method, I get the variable values as: 0.456, 0.4027, 0.1413 respectively (not exactly precise values). The sum of all the variables should be 1. How could I solve such a system? Background: These equations come from solving a Discrete Time Markov Chain. In reality, there could be around 20 such simultaneous equations on 20 variables.


title: Identifying differences among calibration curves: ANCOVA?

text: I've got a large number of calibration curves (i.e. instrument response vs. analyte concentration), and I'd like to check whether they are different from each other. Exactly one curve is associated with each unique combination of the levels of two factors (day and treatment). So, in R, the data look like this:               set.seed(0)     day = gl(2, 6, labels=c("day1", "day2"))     treat = gl(2, 3, 12, labels=c("A", "B"))     conc = rep(0:2, 4)     response = conc + rnorm(12)/2     df <- data.frame(day=day, treat=treat, conc=conc, response=response)      I'd like to know whether one or both of the factors influences the slope of the calibration curve (I don't care so much about the intercept). If they're different, I will use a different calibration curve for each corresponding dataset; if they're not different, I will average them to increase the precision of the measurements that are based on the calibration curves. Is it correct to use ANCOVA to test for differences among the slopes, using the interaction of `day` and `treat`? If so, is the correct R code               model <- lm(response ~ conc + day:treat)      EDIT: To be clear: along with each calibration curve, I have measured data, which I will calibrate using the calibration curve. I am trying to understand whether I should use a separate calibration curve for each timepoint & treatment (which would improve accuracy, if there is an effect of timepoint and/or treatment) or whether it is acceptable to average the calibration curves (which would improve precision). Therefore I'd like to know whether treatment, timepoint, or treatment and timepoint influence the slope of the calibration curve. Actual data as a csv file is here (the number of days will increase in the future). EDIT2: Thanks to all involved in the conversation below. To clarify, I'm trying to understand what model will tell me whether there are significant differences in the slopes of the response-vs-conc curve, by either day OR treatment OR both. It might be instructive to look at my actual data, posted here, rather than the fake data I generated for simplicity.               require(ggplot2)     df <- read.csv("http://dl.dropbox.com/u/48901983/IntrxnDataSo.csv")     ggplot(df, aes(x=conc, y=response, colour=day)) +             geom_point() + geom_line() +            facet_wrap(~treatment)      ![calib_by_treatment](http://i.stack.imgur.com/GutB7.png) So, it looks like there is an effect of day on the slope of the response-vs- conc plot.               ggplot(df, aes(x=conc, y=response, colour=treatment)) +             geom_point() + geom_line() +            facet_wrap(~day)      ![calib_by_day](http://i.stack.imgur.com/VYFJ5.png) It is not as clear that there is an effect of treatment on the slope of the response-vs-conc plot.               int <- ddply(df, .(treatment, day), function(x) coefficients(lm(x$response~x$conc))[1])     names(int)[3] <- "Intercept"     ggplot(int, aes(x=treatment, y=Intercept, colour=day)) + geom_point()      ![enter image description here](http://i.stack.imgur.com/9Xsvj.png) It also looks like there is an effect of treatment on the intercept of the response-vs-conc plot.


title: Is it possible to estimate the odds of winning a multi-entry contest, when I don't know the breakdown of entries?

text: Suppose I am entered into a contest, with the following rules:   * Every person may get up to 6 entries   * All the entries will be pooled, and 25% of the entries will be selected to be winners, with a maximum of 25.   * Each person can only win once, regardless of the number of their entries. If someone's name gets drawn again, it is discarded and a new name drawn.   * I know how many entries I have (the maximum, 6)   * I know how many total entries there are, broken down by type of entry   * I do _not_ know how many of the entries are repeat entries by the same person. The count of entries by type is as follows: > Type 1: 42 Type 2: 72 Type 3: 119 Type 4: 217 Type 5: 156 Type 6: 178 Is it possible to estimate my odds of winning in this situation? I'm a bit confused by the fact that I can't predict how the early winners will effect my chances, since I don't know how many entries each winner will remove from the pool. I'm interested in the solution given the data set, but I'm also interested in the proper procedure/algorithm for calculating it.


title: How can I compute Pearson's $\chi^2$ test statistic for lack of fit on a logistic regression model in R?

text: The likelihood ratio (a.k.a. deviance) $G^2$ statistic and lack-of-fit (or goodness-of-fit) test is fairly straightforward to obtain for a logistic regression model (fit using the `glm(..., family = binomial)` function) in R. However, it can be easy to have some cell counts end up low enough that the test is unreliable. One way to verify the reliability of the likelihood ratio test for lack of fit is to compare its test statistic and _P_ -value to those of Pearson's chi square (or $\chi^2$) lack-of-fit test. Neither the `glm` object nor its `summary()` method report the test statistic for Pearson's chi square test for lack of fit. In my search, the only thing I came up with is the `chisq.test()` function (in the `stats` package): its documentation says "`chisq.test` performs chi-squared contingency table tests and goodness-of-fit tests." However, the documentation is sparse on how to perform such tests: > If `x` is a matrix with one row or column, or if `x` is a vector and `y` is > not given, then a _goodness-of-fit_ test is performed (`x` is treated as a > one-dimensional contingency table). The entries of `x` must be non-negative > integers. In this case, the hypothesis tested is whether the population > probabilities equal those in `p`, or are all equal if `p` is not given. I'd imagine that you could use the `y` component of the `glm` object for the `x` argument of `chisq.test`. However, you can't use the `fitted.values` component of the `glm` object for the `p` argument of `chisq.test`, because you'll get an error: "`probabilities must sum to 1.`" How can I (in R) at least calculate the Pearson $\chi^2$ test statistic for lack of fit without having to run through the steps manually?


title: How do I deal with a violation of sample independence in a Pearson's correlation?

text: There was a concern raised by reviewers of a manuscript of mine over the proper execution of a Pearson's correlation. In brief, this was undertaken in order to determine the relationship between the extent of wheel running (y axis) and ethanol intake (x axis) across three, separate 10 day periods in 7 animals. In the paper, the correlational plots for each 10 day-period had 70 data points: One point for each day and each animal across 10 days of experimentation. The reviewers, however, appropriately pointed out that this is a violation of the assumption of sample independence for Pearson's test, and I should have had only 7 points, which would reflect the means of my two variables for each individual animal across 10 days. Is this appropriate or is there a means of accounting for repeated sampling with a correlation test?


title: External validation of a regression in Stata

text: What I'm essentially trying to do is a temporal external validation of a Cox Proportional Hazards model and also a logistic regression model on the newest year of a dataset that was not included during the development of these models. What I'd like to do is apply these models to the validation dataset and then obtain the new C-statistic/AUROC and Pseudo $R^2$ values for comparison. The question I have is once I run my regressions, how do I actually apply that model to the training dataset in Stata (using version 12), so that I can then work on calculating these values?


title: Within scatter matrix linear discriminant analysis

text: I am trying to implement Linear Discriminant Analysis. I have 10 classes and each class has 3 observations at various instances:               class 1 = {{a1,a2,a3}               {b1,b2,b3}               {c1,c2,c3}}      `a,b,c` are 3 observations found at various instances `a1,a2,a3`. Class 1 is a 3*3 Matrix! Now I have to find the mean observation of each class. For example, I have to find the mean of class 1:               A = (a1+a2+a3)/3     B = (b1+b2+b3)/3     C = (c1+c2+c3)/3     mean of class1 = (A+B+C)/3      I am confused at this point, kindly help me to solve this?


title: How to put values over bars in barplot in R

text: This question is related to my previous question. I would like to put values over bars in barplot. I am beginner in plotting in R.


title: How to interpret the characteristic roots of moment equation of a AR(2) model?

text: I am learning the financial time series using the book 'Analysis of financial time series' by Ruey Tsay. In chapter 2, they introduced AR(2) models. The moment equation (which is the function between the autocorrelations), is a second-order polynomials. And by solving the equation, if there is a real solution with two roots, he said the AR(2) model can be regarded as an AR(1) model operates on top of another AR(1) model.And the ACF of r is then a mixture of two exponential decays. I kind of get the idea, because you can write the moment equation in this way: (1-w1B)(1-w2B)=0 Compare to AR(1), its 1-w1B = 0, (B is the back-shift operator, or you can use L, the lag operator, the same thing), and AR(1) model, ACF is an exponential decay. So for AR(2), this looks like two exponential decays. But my question is, I didn't learn mixture of exponentials before, how do I know the relation between the lag-3 autocorrelation and lag-2, or lag-1 from the roots w1, and w2? Or put it this way, what is the parameters for the mixtured exponentials from w1, w2 I calculated?


title: Why does statistical significance increase with data, BUT the effects may not be meaningful?

text: I have read that _"When you get more and more data, you can find statistically significant differences **wherever you look** "_ My question is:   1. Why is this the case? (any intuitive **examples** that show this behavior?)   2. Why do such increases in statistical difference do not necessarily imply that the observed effects are meaningful / important?


title: How to compare two groups with multiple measurements for each individual with R?

text: I have a problem like the following: 1) There are six measurements for each individual with large within-subject variance 2) There are two groups (Treatment and Control) 3) Each group consists of 5 individuals 4) I want to perform a significance test comparing the two groups to know if the group means are different from one another. The data looks like this: ![http://s10.postimg.org/p9krg6f3t/examp.png](http://i.stack.imgur.com/55V9J.png) And I have run some simulations using this code which does t tests to compare the group means. The group means were calculated by taking the means of the individual means. **This ignores within-subject variability** :                n.simulations<-10000         pvals=matrix(nrow=n.simulations,ncol=1)         for(k in 1:n.simulations){           subject=NULL           for(i in 1:10){             subject<-rbind(subject,as.matrix(rep(i,6)))           }           #set.seed(42)                #Sample Subject Means           subject.means<-rnorm(10,100,2)                #Sample Individual Measurements           values=NULL           for(sm in subject.means){             values<-rbind(values,as.matrix(rnorm(6,sm,20)))           }                out<-cbind(subject,values)                #Split into GroupA and GroupB           GroupA<-out[1:30,]           GroupB<-out[31:60,]                #Add effect size to GroupA           GroupA[,2]<-GroupA[,2]+0                colnames(GroupA)<-c("Subject", "Value")           colnames(GroupB)<-c("Subject", "Value")                #Calculate Individual Means and SDS           GroupA.summary=matrix(nrow=length(unique(GroupA[,1])), ncol=2)           for(i in 1:length(unique(GroupA[,1]))){             GroupA.summary[i,1]<-mean(GroupA[which(GroupA[,1]==unique(GroupA[,1])[i]),2])             GroupA.summary[i,2]<-sd(GroupA[which(GroupA[,1]==unique(GroupA[,1])[i]),2])           }           colnames(GroupA.summary)<-c("Mean","SD")                     GroupB.summary=matrix(nrow=length(unique(GroupB[,1])), ncol=2)           for(i in 1:length(unique(GroupB[,1]))){             GroupB.summary[i,1]<-mean(GroupB[which(GroupB[,1]==unique(GroupB[,1])[i]),2])             GroupB.summary[i,2]<-sd(GroupB[which(GroupB[,1]==unique(GroupB[,1])[i]),2])           }           colnames(GroupB.summary)<-c("Mean","SD")                Summary<-rbind(cbind(1,GroupA.summary),cbind(2,GroupB.summary))           colnames(Summary)[1]<-"Group"                pvals[k]<-t.test(GroupA.summary[,1],GroupB.summary[,1], var.equal=T)$p.value         }      And here is code for plots:               #Plots     par(mfrow=c(2,2))     boxplot(GroupA[,2]~GroupA[,1], col="Red", main="Group A",              ylim=c(.9*min(out[,2]),1.1*max(out[,2])),             xlab="Subject", ylab="Value")     stripchart(GroupA[,2]~GroupA[,1], vert=T, pch=16, add=T)     #abline(h=mean(GroupA[,2]), lty=2, lwd=3)          for(i in 1:length(unique(GroupA[,1]))){       m<-mean(GroupA[which(GroupA[,1]==unique(GroupA[,1])[i]),2])       ci<-t.test(GroupA[which(GroupA[,1]==unique(GroupA[,1])[i]),2])$conf.int[1:2]            points(i-.2,m, pch=15,cex=1.5, col="Grey")       segments(i-.2,                ci[1],i-.2,                ci[2], lwd=4, col="Grey"       )     }     legend("topleft", legend=c("Individual Means +/- 95% CI"), bty="n", pch=15, lwd=3, col="Grey")               boxplot(GroupB[,2]~GroupB[,1], col="Light Blue", main="Group B",              ylim=c(.9*min(out[,2]),1.1*max(out[,2])),             xlab="Subject", ylab="Value")     stripchart(GroupB[,2]~GroupB[,1], vert=T, pch=16, add=T)     #abline(h=mean(GroupB[,2]), lty=2, lwd=3)          for(i in 1:length(unique(GroupB[,1]))){       m<-mean(GroupB[which(GroupB[,1]==unique(GroupB[,1])[i]),2])       ci<-t.test(GroupB[which(GroupB[,1]==unique(GroupB[,1])[i]),2])$conf.int[1:2]            points(i-.2,m, pch=15,cex=1.5, col="Grey")       segments(i-.2,                ci[1],i-.2,                ci[2], lwd=4, col="Grey"       )     }     legend("topleft", legend=c("Individual Means +/- 95% CI"), bty="n", pch=15, lwd=3, col="Grey")               boxplot(Summary[,2]~Summary[,1], col=c("Red","Light Blue"), xlab="Group", ylab="Average Value",             ylim=c(.9*min(Summary[,2]),1.1*max(Summary[,2])),             main="Individual Averages")     stripchart(Summary[,2]~Summary[,1], vert=T, pch=16, add=T)          points(.9, mean(GroupA.summary[,1]), pch=15,cex=1.5, col="Grey")     segments(.9,              t.test(GroupA.summary[,1])$conf.int[1],.9,                  t.test(GroupA.summary[,1])$conf.int[2], lwd=4, col="Grey"     )          points(1.9, mean(GroupB.summary[,1]), pch=15,cex=1.5, col="Grey")     segments(1.9,              t.test(GroupB.summary[,1])$conf.int[1],1.9,                  t.test(GroupB.summary[,1])$conf.int[2], lwd=4, col="Grey"     )     legend("topleft", legend=c("Group Means +/- 95% CI"), bty="n", pch=15, lwd=3, col="Grey")               hist(pvals, breaks=seq(0,1,by=.05), col="Grey",          main=c(paste("# sims=", n.simulations),                 paste("% Sig p-values=",100*length(which(pvals<0.05))/length(pvals)))     )      Now, it seems to me that because each individual mean is an estimate itself, that we should be less certain about the group means than shown by the 95% confidence intervals indicated by the bottom-left panel in the figure above. Thus the p-values calculated are underestimating the true variability and should lead to increased false-positives if we wish to extrapolate to future data. So what is the correct way to analyze this data? **Bonus:** The example above is a simplification. For the actual data: 1) The within-subject variance is positively correlated with the mean. 2) Values can only be multiples of two. 3) The individual results are not roughly normally distributed. They suffer from zero floor effect, and have long tails at the positive end. 4) Number of Subjects in each group are not necessarily equal. Previous literature has used the t-test ignoring within-subject variability and other nuances as was done for the simulations above. Are these results reliable? If I can extract some means and standard errors from the figures how would I calculate the "correct" p-values. **EDIT:** Ok, here is what _actual_ data looks like. There is also three groups rather than two: ![enter image description here](http://i.stack.imgur.com/k1xWd.png) dput() of data:               structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,      1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,      2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,      2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,      3, 3, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3,      3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6,      6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10,      10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12,      12, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 15, 15, 15,      15, 15, 15, 16, 16, 16, 16, 16, 16, 17, 17, 17, 17, 17, 17, 18,      18, 18, 18, 18, 18, 2, 0, 16, 2, 16, 2, 8, 10, 8, 6, 4, 4, 8,      22, 12, 24, 16, 8, 24, 22, 6, 10, 10, 14, 8, 18, 8, 14, 8, 20,      6, 16, 6, 6, 16, 4, 2, 14, 12, 10, 4, 10, 10, 8, 4, 10, 16, 16,      2, 8, 4, 0, 0, 2, 16, 10, 16, 12, 14, 12, 8, 10, 12, 8, 14, 8,      12, 20, 8, 14, 2, 4, 8, 16, 10, 14, 8, 14, 12, 8, 14, 4, 8, 8,      10, 4, 8, 20, 8, 12, 12, 22, 14, 12, 26, 32, 22, 10, 16, 26,      20, 12, 16, 20, 18, 8, 10, 26), .Dim = c(108L, 3L), .Dimnames = list(         NULL, c("Group", "Subject", "Value")))      **EDIT 2:** In response to Henrik's answer: So if I instead perform anova followed by TukeyHSD procedure on the individual averages as shown below, I could interpret this as underestimating my p-value by about 3-4x? My goal with this part of the question is to understand how I, as a reader of a journal article, can better interpret previous results given their choice of analysis method. For example they have those "stars of authority" showing me 0.01>p>.001. So if i accept 0.05 as a reasonable cutoff I should accept their interpretation? The only additional information is mean and SEM.               #Get Invidual Means     summary=NULL     for(i in unique(dat[,2])){     sub<-which(dat[,2]==i)     summary<-rbind(summary,cbind(     dat[sub,1][3],     dat[sub,2][4],     mean(dat[sub,3]),     sd(dat[sub,3])     )     )     }     colnames(summary)<-c("Group","Subject","Mean","SD")          TukeyHSD(aov(summary[,3]~as.factor(summary[,1])+ (1|summary[,2])))          #      Tukey multiple comparisons of means     #        95% family-wise confidence level     #         #    Fit: aov(formula = summary[, 3] ~ as.factor(summary[, 1]) + (1 | summary[, 2]))     #         #    $`as.factor(summary[, 1])`     #             diff       lwr       upr     p adj     #    2-1 -0.672619 -4.943205  3.597967 0.9124024     #    3-1  7.507937  1.813822 13.202051 0.0098935     #    3-2  8.180556  2.594226 13.766885 0.0046312      **EDIT 3:** I think we are getting close to my understanding. Here is the simulation described in the comments to @Stephane:               #Get Subject Means     means<-aggregate(Value~Group+Subject, data=dat, FUN=mean)          #Initialize "dat2" dataframe     dat2<-dat          #Initialize within-Subject sd     s<-.001     pvals=matrix(nrow=10000,ncol=2)          for(j in 1:10000){     #Sample individual measurements for each subject     temp=NULL     for(i in 1:nrow(means)){     temp<-c(temp,rnorm(6,means[i,3], s))     }          #Set new values     dat2[,3]<-temp          #Take means of sampled values and fit to model     dd2 <- aggregate(Value~Group+Subject, data=dat2, FUN=mean)     fit2 <- lm(Value~Group, data=dd2)          #Save sd and pvalue     pvals[j,]<-cbind(s,anova(fit2)[[5]][5])          #Update sd     s<-s+.001     }          plot(pvals[,1],pvals[,2], xlab="Within-Subject SD", ylab="P-value")      ![enter image description here](http://i.stack.imgur.com/gMMDY.png)


title: Does ARIMA require normally distributed errors or normally distributed input data?

text: I have two questions related to time series forecasting with ARIMA:   1. Does ARIMA require normally distributed errors or normally distributed input data ?   2. Are there any assumptions on input time series data for ARIMA model and exogenous variables for ARIMAX model ?


title: What's the approximate distribution? Replace the true mean with sample mean

text: If say for a random variable $X$, I have observation of $x_1,x_2,x_3,\ldots,x_n$. Let $m$ be the sample mean, and $s$ be the sample standard deviation. Does the new random variable $(X-m)/s$ follow some distribution? It's not $t$-distribution I guess, since for it to be $t$ distributed, $m$ needs to be replaced by true mean. Can statistics expert shed some light on this?


title: Interaction wipes out my direct effects in regression (non zero variable)

text: I have the following regression $children = \beta_0 + \beta_1 \log(earnings) + \beta_2 grandparents + \epsilon$ and $\beta_1>0$ with $p$=0.01 and $\beta_2>0$ with $p$=0.01, and N is large (N>10.000) and grandparents takes values 0,1,2,3,4. Then I add the interaction term ($\log(earnings)*grandparents$) to equation 1, such that: $children = \beta_0 + \beta_1 \log( earnings) + \beta_2 grandparents+ \beta_3 \log( earnings)*grandparents + \epsilon$ and $\beta_1>0$ with $p$=0.01, $\beta_2$ is no longer statistically significant and also $\beta_3$ is not statistically significant. I do not understand how to interpret the results and if the interaction term wipes out the direct effect of grandparents since $\log(earnings)$ is always different from 0. > > > > There is a way to test the stat. sign. of the effect of Grandparents > in the interacted model? (Thanks Maarten for your previous answer)


title: Statistical test to compare two ratios from two independent models

text: What statistical test can I use to compare two ratios from two independent samples. The ratios are after to before results. I need to compare the after/before ratios for two independent models and show whether they are have significant difference or not. Please help!


title: Unbalanced mixed effect ANOVA for repeated measures

text: I have data from patients treated with 2 different kinds of treatments during surgery. I need to analyze its effect on heart rate. The heart rate measurement is taken every 15 minutes. Given that the surgery length can be different for each patient, each patient can have between 7 and 10 heart rate measurements. So an unbalanced design should be used. I'm doing my analysis using R. And have been using the ez package to do repeated measure mixed effect ANOVA. But I do not know how to analyse unbalanced data. Can anyone help? Suggestions on how to analyze the data are also welcomed. Update:   As suggested, I fitted the data using the `lmer` function and found that the best model is:               heart.rate~ time + treatment + (1|id) + (0+time|id) + (0+treatment|time)      with the following result:               Random effects:      Groups   Name        Variance   Std.Dev. Corr         id       time        0.00037139 0.019271              id       (Intercept) 9.77814104 3.127002              time     treat0      0.09981062 0.315928                       treat1      1.82667634 1.351546 -0.504       Residual             2.70163305 1.643665             Number of obs: 378, groups: subj, 60; time, 9          Fixed effects:                  Estimate Std. Error t value     (Intercept) 72.786396   0.649285  112.10     time         0.040714   0.005378    7.57     treat1       2.209312   1.040471    2.12          Correlation of Fixed Effects:            (Intr) time       time   -0.302            treat1 -0.575 -0.121      Now I'm lost at interpreting the result. Am I right in concluding that the two treatments differed in affecting heart rate? What does the correlation of -504 between treat0 and treat1 means?