title: Significance of the component in a multivariate linear model ("within-design")

text: I'm trying to test the significance of the "component" effect in a multivariate regression model. I'm not sure what is the right way. Using R, I have tried a way with `lm()` and another way with `gls()`, and they don't yield compatible results. **Please note that this is not a question about which methodology is the right one to use to analyze my data. By the way these are simulated data. My question is about the understanding in mathematical terms of the R procedures I use.** The dataset:               > str(dat)     'data.frame':   31 obs. of  5 variables:      $ group: Factor w/ 5 levels "1","5","2","3",..: 1 1 1 1 1 1 1 1 3 3 ...          $ id   : Factor w/ 8 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 1 3 ...      $ x    : num  2.5 3 3 4 1.2 3.8 3.9 4 2.5 2.9 ...          $ y    : num  2.6 3.8 3.9 3.8 1.6 5.2 1.3 3.6 4 3.2 ...      $ z    : num  3.1 3.6 4.9 3.8 2.1 6 2.1 2.9 4.2 2.9 ...     > head(dat,10)        group id   x   y   z     1      1  1 2.5 2.6 3.1     2      1  2 3.0 3.8 3.6     3      1  3 3.0 3.9 4.9     4      1  4 4.0 3.8 3.8     5      1  5 1.2 1.6 2.1     6      1  6 3.8 5.2 6.0     7      1  7 3.9 1.3 2.1     8      1  8 4.0 3.6 2.9     9      2  1 2.5 4.0 4.2     10     2  3 2.9 3.2 2.9      I convert this dataset into "long format" for graphics (and later for `gls()`):               dat$subject <- dat$group : dat$id         dat.long <- reshape(dat, direction="long", varying=list(3:5),           idvar="subject", v.names="value", timevar="component", times=c("x","y","z"))         dat.long$component <- factor(dat.long$component)          xyplot(value ~ component | group, data=dat.long,          pch=16,          strip = strip.custom(strip.names=TRUE,var.name="group" ), layout=c(5,1))      ![enter image description here](http://i.stack.imgur.com/14KtA.png) Each individual of each group has $3$ repeated measures $x$,$y$,$z$ (I should join the points in the graphic to see the repeated measures). I want to fit a MANOVA model using group as factor and $(x,y,z)$ is the multivariate response: $$ \begin{pmatrix} x_{ij} \\ y_{ij} \\ z_{ij} \end{pmatrix} \sim {\cal N}_3\left( \begin{matrix} \mu_{i1} \\ \mu_{i2} \\ \mu_{i3} \end{matrix}, \Sigma \right), \quad i=1,\ldots,5 $$ (of course we could use the default R parameterization $\mu_{ik}=\mu_{1k} + \alpha_{ik}$ by considering `group1`as the "intercept" for each response but I prefer "my" parameterization). This model is fitted as follows using `lm()`:               ###  multivariate least-squares fitting  ###     mfit <- lm( cbind(x,y,z)~group, data=dat )      I think the model can also be fitted with `gls()` as follows (but with a different fitting procedure) :               ### generalized least-squares fitting  ###     library(nlme)     gfit <- gls(value ~ group*component, data=dat.long, correlation=corSymm(form= ~ 1 | subject))      Recall that `subject = group:id` is the identifier of the individuals. The `correlation=corSymm(form= ~ 1 | subject)` argument means that the responses $x$, $y$, $z$ for each individual are correlated. Here `corSymm` means a general, "unrestricted", covariance structure (termed as "unstructured" in SAS language). To check that `mfit` and `gfit` are equivalent, we can check for instance that we can deduce the estimated parameters of `mfit` from the estimated parameters of `gfit`and vice-versa (so the "mean" parameters have exactly the same fitted values):               > coef(mfit)                       x          y          z     (Intercept)  3.1750  3.2250000  3.5625000     group5      -0.9500 -0.4750000  0.1125000     group2      -1.0750 -0.5678571 -0.2339286     group3      -0.7875 -0.1000000  0.1875000     group4      -0.3750  0.4000000 -0.0125000     > coef(gfit)           (Intercept)            group5            group2            group3              3.1750000        -0.9500000        -1.0750000        -0.7875000                 group4        componenty        componentz group5:componenty             -0.3750000         0.0500000         0.3875000         0.4750000      group2:componenty group3:componenty group4:componenty group5:componentz              0.5071429         0.6875000         0.7750000         1.0625000      group2:componentz group3:componentz group4:componentz              0.8410714         0.9750000         0.3625000       Now I want to test the "component effect". Rigorously speaking, writing the model as $$ \begin{pmatrix} x_{ij} \\ y_{ij} \\ z_{ij} \end{pmatrix} \sim {\cal N}_3\left( \begin{matrix} \mu_{i1} \\ \mu_{i2} \\ \mu_{i3} \end{matrix}, \Sigma \right), $$ I want to test the hypothesis $\boxed{H_0\colon \{\mu_{1i}=\mu_{2i}=\mu_{3i} \quad \forall i=1,2,3,4,5 \}}$. Below are my attempts, one attempt with `gfit` and two attempts with `mfit()`:               ###########################################     ## testing significance of the component ##     ###########################################          > ### with gfit  ###     > anova(gfit)     Denom. DF: 78                      numDF  F-value p-value     (Intercept)         1 504.5226  <.0001     group               4   0.7797  0.5418     component           2  11.7073  <.0001     group:component     8   0.7978  0.6063     >      > ### with mfit ###     > library(car)     >      > # first attempt :      > idata <- data.frame(component=c("x","y","z"))     > ( av.ok <- Anova(mfit, idata=idata, idesign=~component, type="III") )          Type III Repeated Measures MANOVA Tests: Pillai test statistic                     Df test stat approx F num Df den Df    Pr(>F)         (Intercept)      1   0.84396  140.625      1     26 5.449e-12 ***     group            4   0.10369    0.752      4     26    0.5658         component        1   0.04913    0.646      2     25    0.5328         group:component  4   0.22360    0.818      8     52    0.5901         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      >      > # second attempt :     > linearHypothesis(mfit, "(Intercept) = 0", idata=idata, idesign=~component, iterms="component")           Response transformation matrix:       component1 component2     x          1          0     y          0          1     z         -1         -1          Sum of squares and products for the hypothesis:                component1 component2     component1    1.20125    1.04625     component2    1.04625    0.91125          Sum of squares and products for error:                component1 component2     component1   31.46179   14.67696     component2   14.67696   21.42304          Multivariate Tests:                       Df test stat  approx F num Df den Df  Pr(>F)     Pillai            1 0.0491253 0.6457903      2     25 0.53277     Wilks             1 0.9508747 0.6457903      2     25 0.53277     Hotelling-Lawley  1 0.0516632 0.6457903      2     25 0.53277     Roy               1 0.0516632 0.6457903      2     25 0.53277      With `anova(gfit)` the component is significant, but not with my two attempts using `mfit` and the `car` package. I know that `gls()` use a different fitting method than `lm()` but this is surely not the cause of the difference. So my questions are :   * did I do something wrong ?   * which method tests my $H_0$ hypothesis ?   * what is the $H_0$ hypothesis of the other methods ? And I have an auxiliary question: how to get $\hat\Sigma$ with `mfit` and `gfit` ? ## Update 1 Below is a reproducible example which simulates the dataset. Now I think I understand : both ANOVA methods are correct (the first one with `anova(gfit)` and the second one with `Anova(mfit, ...)`, **and they yield very close results when using the type II sum of squares in`Anova(mfit, ...)`**. For the above example:               > Anova(mfit, idata=idata, idesign=~component, type="II")          Type II Repeated Measures MANOVA Tests: Pillai test statistic                     Df test stat approx F num Df den Df    Pr(>F)         (Intercept)      1   0.94691   463.72      1     26 < 2.2e-16 ***     group            4   0.10369     0.75      4     26 0.5657691         component        1   0.50344    12.67      2     25 0.0001584 ***     group:component  4   0.22360     0.82      8     52 0.5900848         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1       is very close to               > anova(gfit)     Denom. DF: 78                      numDF  F-value p-value     (Intercept)         1 504.5226  <.0001     group               4   0.7797  0.5418     component           2  11.7073  <.0001     group:component     8   0.7978  0.6063      Below is the reproducible code with the data sampler (I simulate uncorrelated repeated measures but it suffices to include a covariance matrix in the `rmvnorm()` function to simulate correlated repeated measures) :               library(mvtnorm)     library(nlme)     library(car)          # set data parameters      I <- 5 # number of groups     J <- 16 # number of individuals per group     dat <- data.frame(         group = gl(I,J),         id = gl(J,1,I*J),         x=NA,          y=NA,          z=NA     )     Mu <- c(1:I) # group means of components (assuming E(x)=E(y)=E(z) in each group)          # simulates data:      for(i in 1:I){         which.group.i <- which(dat$group==i)         dat[which.group.i,c("x","y","z")] <- round(rmvnorm(n=J, mean=rep(Mu[i],3)), 1)     }          dat$subject <- droplevels( dat$group : dat$id )         dat.long <- reshape(dat, direction="long", varying=list(3:5),           idvar="subject", v.names="value", timevar="component", times=c("x","y","z"))         dat.long$component <- factor(dat.long$component)          # multivariate least-squares fitting      mfit <- lm( cbind(x,y,z)~group, data=dat )          # gls fitting     dat.long$order.xyz <- as.numeric(dat.long$component)     gfit <- gls(value ~ group*component , data=dat.long, correlation=corSymm(form=  ~ order.xyz | subject))           # compares ANOVA :      anova(gfit)     idata <- data.frame(component=c("x","y","z"))     Anova(mfit, idata=idata, idesign=~component, type="II")     Anova(mfit, idata=idata, idesign=~component, type="III")      So now I wonder which type of sum of squares is the more appropriate one for my real study... but this is another question ## Update 2 About my question _"how to get $\hat\Sigma$"_ , here is the answer for `gls()`:               > getVarCov(gfit)     Marginal variance covariance matrix             [,1]    [,2]    [,3]     [1,] 0.92909 0.47739 0.24628     [2,] 0.47739 0.92909 0.53369     [3,] 0.24628 0.53369 0.92909       Standard Deviations: 0.96389 0.96389 0.96389       That shows that **`mfit`and `gfit` were not equivalent models**: `gfit`assumes the same variance for the three components. In order to fit a fully unrestricted covariance matrix for the repeated measures, we have to type:               gfit2 <- gls(value ~ group*component , data=dat.long,          correlation=corSymm(form=  ~ 1 | subject),          weights=varIdent(form = ~1 | component))          > summary(gfit2)     Generalized least squares fit by REML       Model: value ~ group * component        Data: dat.long             AIC      BIC    logLik       264.0077 313.4986 -111.0038          Correlation Structure: General      Formula: ~1 | subject       Parameter estimate(s):      Correlation:        1     2         2 0.529           3 0.300 0.616     Variance function:      Structure: Different standard deviations per stratum      Formula: ~1 | component       Parameter estimates:            x        y        z      1.000000 1.253534 1.169335           ....          Residual standard error: 0.8523997       But yet I don't understand the extracted covariance matrix given by `getVarCov()` (but this is not important since we get this matrix with `summary(gfit2)`):                  > getVarCov(gfit2)         Error in t(S * sqrt(vars)) :            dims [product 9] do not match the length of object [0]         > getVarCov(gfit2, individual="1:1")         Marginal variance covariance matrix                 [,1]    [,2]    [,3]         [1,] 0.72659 0.48164 0.25500         [2,] 0.48164 1.14170 0.65562         [3,] 0.25500 0.65562 0.99349           Standard Deviations: 0.8524 1.0685 0.99674          > getVarCov(gfit2, individual="1:2")         Marginal variance covariance matrix                 [,1]    [,2]    [,3]         [1,] 1.14170 0.56319 0.27337         [2,] 0.56319 0.99349 0.52302         [3,] 0.27337 0.52302 0.72659           Standard Deviations: 1.0685 0.99674 0.8524       Unfortunately, the `anova(gfit2)` table is not as close to `Anova(mfit, ..., type="II")` as `anova(gfit)`:               > anova(gfit2)     Denom. DF: 78                      numDF  F-value p-value     (Intercept)         1 498.1744  <.0001     group               4   1.0514  0.3864     component           2  13.1801  <.0001     group:component     8   0.8310  0.5780          > Anova(mfit, idata=idata, idesign=~component, type="II")          Type II Repeated Measures MANOVA Tests: Pillai test statistic                     Df test stat approx F num Df den Df    Pr(>F)         (Intercept)      1   0.94691   463.72      1     26 < 2.2e-16 ***     group            4   0.10369     0.75      4     26 0.5657691         component        1   0.50344    12.67      2     25 0.0001584 ***     group:component  4   0.22360     0.82      8     52 0.5900848         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           > anova(gfit)     Denom. DF: 78                      numDF  F-value p-value     (Intercept)         1 504.5226  <.0001     group               4   0.7797  0.5418     component           2  11.7073  <.0001     group:component     8   0.7978  0.6063


title: Is it possible to test for a Gaussian distribution?

text: Suppose that I have $n$ independent samples from a probability distribution. I believe that the distribution is Gaussian ($G$) but I don't know; then my question is: does there exist a test such that I can interpret the result of the test to be something like "the probability that my samples are the output of a non-Gaussian distribution T with $d(G,T)>\epsilon$ is smaller than $x$%", where $d$ is the variational distance. I have Googled and found the Kolmogorov-Smirnov test and some others, but I'm not sure if I can use the results of those tests to make a strong statement like the one above.


title: What are the standard hypothesis tests in any type of regression?

text: In general, what the the standard hypothesis tests in all types of regression? We get some coefficients and then test the hypotheses that these coefficients are significant (i.e not equal to $0$)? For example, in logistic regression the p-value indicates that the regression coefficients are significant or not? Likewise with linear regression? > Question: Suppose coefficients $\beta_0$ and $\beta_1$ are significant in > linear regression. Can we deduce anything about the significance of > $\beta_0+\beta_1$? What about $\beta_0= \beta_1$?


title: How do I get "V-shaped" distributed random numbers from uniformly distributed numbers?

text: I have 1000 uniformly distributed random numbers. How do I manipulate them to get a V-shaped histogram?


title: Time series (stochastic process) estimating parameters using characteristic function

text: I have a time series of assets ${A_1, A_2, ..., A_n}$, which is described by a sophisticated distribution having the following characteristic function: $\phi(u; t;\theta)$, where $\theta$ is a vector of unknown parameters. I need to esimate vector of unknown parameters $\theta$ of characteristic function. I tried to find a PDF using the inverse Fourier transform to use the maximum likelihood method, but the characteristic function is too complicated for that. I also thought about building the empirical characteristic function using the time series of assets and to estimate parameters using the least square method, but I do not know how to build the empirical characteristic function, because time series is not just a sample of random variables, it is a random process that depends on time.


title: How to interpret a box plot?

text: I have some data in which there are 5 categorical explanatory variables (`concern, breath, weath, sleep, act`) and 1 continuous response variable (`tto`). In addition, every categorical explanatory variable is divided into 5 levels which show how strong a person feels about it. level 1 and level 5 show the perfect and worst states respectively. I was advised to create a box plot to see the relationship between the explanatory variables and the response variable. The plot is given below. However, I do not know how to read a box plot. Can any one please help me interpret it. ![enter image description here](http://i.stack.imgur.com/HmxBX.png)


title: Advantages and disadvantages of parametric and non-parametric models

text: I have a large simulated loss data (from catastrophic models developed at my school) to calculate some extreme quantiles. Previously they used non- parametric methods to do this (find the point estimate for these extreme quantiles and CIs). I am using parametric models (extreme value theory, fat tail distributions, etc.) to do it. I have been thinking about the pros and cons for these two methods. I would appreciate if someone could provide some summaries of parametric and non-parametric models, their advantages and disadvantages.


title: Symmetrised Kullback - Leibler divergence

text: I have trouble understanding KL divergence, where P is probability mass function of true distribution of data and Q is the approximation of P. The definition of KL divergence is: $$D_{\mathrm{kl}}(P || Q) = \sum_i \ln\left( \frac{P_i}{Q_i}\right)P_i$$ If I want a symmetrised KL divergence, should it look like the following ? $$D_{\mathrm{kl}}(P || Q) + D_{\mathrm{kl}}(Q||P) = \sum_i \ln\left( \frac{P_i}{Q_i}\right)P_i + \sum_i \ln\left( \frac{Q_i}{P_i}\right)Q_i$$ Also, which package should I use in R to compute the KL divergence for **discrete distributions** ? Flexmix or FNN ? Or should I just write my own R function for this? And what if I don't really know $P$, the true distribution of data?


title: Why do the 95% confidence limits in ARIMA models widen at the forecasts?

text: Can someone please explain why when I do an ARIMA model the forecast's 95% confidence interval widen?


title: Test hypothesis "This profit is significantly different than the average"

text: I am writing a paper about Algo-trading and my supervisor asked me to test the hypothesis that "the return of the trading strategy is significantly different than the average". To reach this result, it has been implemented a software doing the following:   1. Compute the profit by trading on a certain stock for a certain period of time according to the pre-set trading rules;    2. Generate a big sample (1,000) of "random trading strategies". For reaching this, it has been implemented a Linear Congruential Generator which is returning random numbers between 0 and 1; these numbers are hence used, over the same stock and over the same time-period, to decide if trading or not and then compute the profit of each trading session.  The output looks like follows: ![enter image description here](http://i.stack.imgur.com/BuS8D.png) I've studied statistics long time ago so, even if I've realized myself this question is very basic, I'm not able to remind.... how can I say if the profit 12.96% got by applying the trading rules is significantly different from the sample mean 13.05% of 10,000 trading sessions randomly executed? **EDIT:** Any further suggestion about the sample-size or other factors might make the statistics more valid would be highly appreciated.


title: LRT yielding a test statistic with extreme behavior

text: I have what seems like a textbook binary hypothesis testing scenario: an observation $X$ drawn from one of two probability distributions $\{\Pi_0,\Pi_1\}$ that are known to me, and a need to classify $X$ as having been drawn either according to $\Pi_0$ or $\Pi_1$ (without loss of generality, let's assume equal prior probabilities). The standard texts on the detection theory (such as Van Trees) prescribe constructing a likelihood ratio test (LRT): $$\Lambda(X)\triangleq\frac{p_{x|H_1}(X|H_1)}{p_{x|H_0}(X|H_0)}$$ where $H_0$ and $H_1$ denote the hypotheses where $X$ is, respectively, either drawn from $\Pi_0$ or $\Pi_1$. Since, as I mentioned, both $\Pi_0$ and $\Pi_1$ are known, and the prior probabilities are equal, the likelihood functions $p_{x|H_1}(X|H_1)$ and $p_{x|H_0}(X|H_0)$ are well-defined. Thus, the test reduces to computing $\Lambda(X)$ and comparing to a threshold. Denote by $E_0[\Lambda(X)]$ and $E_1[\Lambda(X)]$ the expected value of $\Lambda(X)$ corresponding to, respectively, $X\sim\Pi_0$ and $X\sim\Pi_1$. Similarly, denote by $\operatorname{Var}_0[\Lambda(X)]$ and $\operatorname{Var}_1[\Lambda(X)]$ the variance of $\Lambda(X)$ corresponding to, respectively, $X\sim\Pi_0$ and $X\sim\Pi_1$. Now, in my particular scenario, under certain parameters for $\Pi_1$ (and the corresponding likelihood function $p_{x|H_1}(X|H_1)$), I am finding that both $\operatorname{Var}_0[\Lambda(X)]$ and $\operatorname{Var}_1[\Lambda(X)]$ are inifinite. Furthermore, I am finding that with these parameters the expected value of $\Lambda(X)$ when $X\sim\Pi_1$, $E_1[\Lambda(X)]$, is infinite. What does this mean? Per Glen_b's comment, this is not necessarily a trivial detection problem where $\Pi_0$ is always discriminated from $\Pi_1$ without any errors. So, what can be said (if anything) about the hypothesis test between $H_0$ and $H_1$ based on these facts about $\Lambda(X)$?


title: Parametric Anova Test to find differences?

text: H0: Male and Female members do not differ in average Awareness, Attitude and Loyalty Ha: Male and Female members differ in average Awareness, Attitude and Loyalty. STEP 1: Used a Levene's Test to find for equal variances According to the Levene's test - P value across all variables from the Levene’s test of ‘Awareness’ is 0.270, ‘Attitude’ is 0.533 and ‘Loyalty’ is 0.905 and is not significant so we assume that there is an equality of variances. A Levene’s test has verified the equality of variances in the samples. Hence, we cannot reject null hypothesis about equality of variances. The variances for the two populations are the same. STEP 2: Ran a Parametric test. P-value is 0.021 for ‘Awareness’, 0.064 for ‘Attitude’ and 0.372 for ‘Loyalty’. Should I reject the null hypothesis as a whole or only reject Awareness?


title: What factor loading is recommended to retaining items in a factor analysis?

text: _I received the following question via email and thought it would be suited to this site:_ > I have a debate with a friend about factor loading and squared multiple > correlation. ... In my debate with my colleague, I argue that regression > weights/factor loading above 0.6 is enough for a construct. I mean for a > construct which factor loading below 0.6 should be deleted because of the > convergent validity problem. However, the colleague of mine argued that we > also need to see squared multiple correlation. For squared multiple > correlation below 0.7 we have to delete the items. So, now I am getting > confused since in my research dataset, I found many items with squared > multiple correlation below 0.7 but factor loading above 0.6. Should I delete > these items? > > So please advise me in your knowledge whether we really need to see squared- > multiple-correlation or factor loading is enough?. In addition, actually > what is the function/utility of squared-multiple-correlation. Is there a cut > off value for squared-multiple-correlation? Thus, I distil from this email a core question: **What decision rules should be used regarding squared multiple correlations and factor loadings when deciding on whether to retain an item of a question in the context of factor analysis?**


title: Not sure how to interpret "positive" results when using more than one control group

text: So we're starting to do A/B testing in my company (it's a mobile app), and my CEO is pretty skeptical about the whole thing (I guess I can't blame him, especially after the results I'm seeing). One of our first tests is to delay showing ads for 1 day vs. 2 days. The "control" in this case is 2 days, which is what we currently do. Because the first time around in our testing, we made our group sizes too small the results were very noisy with lots of outliers. So this time around, we created big groups (~8K users in each group) and I put in 2 control groups. Here is the output from running aov in R for one of the metrics we are tracking:               Call:     aov(formula = connections ~ as.factor(group_id), data = userdata_exp104)          Residuals:         Min      1Q  Median      3Q     Max       -2.609  -2.262  -2.103  -0.262 113.391           Coefficients:                            Estimate Std. Error t value Pr(>|t|)         (Intercept)              2.6095     0.1275  20.469  < 2e-16 ***     as.factor(group_id)472  -0.5060     0.1799  -2.812  0.00493 **      as.factor(group_id)473  -0.3474     0.1798  -1.933  0.05333 .       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 6.117 on 6947 degrees of freedom     Multiple R-squared:  0.001189,  Adjusted R-squared:  0.0009014      F-statistic: 4.135 on 2 and 6947 DF,  p-value: 0.01604      In this case group_id 473 is actually the test group (with 1 day delay in ads) while the other two groups are the "controls". You can see that one of the controls (group_id 472) appears to be highly statistically significantly different from the other control. I really don't know how to interpret this. In theory, the two control groups should have been the same. I don't know what to tell my CEO, other than we found no difference for the test group. Any ideas?


title: Training and test data for Holt Winters Method

text: I am new to forecasting. I have weekly data. I want to forecast using the Holt-Winters method. Should I make training and test subsets of the data? How important are the training and test data?


title: Variance of a smoothed AR(1) process

text: The query I have relates to calculating the variance of AR(1) processes that are smoothed with a simple moving average. So: In an AR(1) process of the form: $$ X_t=c+\varphi X_{t-1}+\varepsilon_t, $$ the variance can be calculated as: $$ \text{var}(X_t)=E(X^2_t)-\mu^2=\frac{\sigma^2_{\varepsilon}}{1-\varphi^2}, $$ where $\sigma$ is the standard deviation of $\varepsilon$ (a white noise), and $\varphi$ is the variable defining the autocorrelation properties of the AR(1) process. I would like to be able to calculate the variance of the AR(1) process after smoothing by a simple unweighted moving average for various window sizes. I have thus far looked at the problem analytically with data, and as the window size of the moving average increases the variance obviously falls, and the way in which it falls (i.e. rate and shape) is dependent on $\varphi$. Clearly, when the moving average window is equal to $N$ (the size of the dataset analysed) then the variance is 0. Ultimately therefore, is there a way of determining the expected variance of the AR(1) in terms of $\varphi$, $N$, and the size of the moving average window? Any pointers gratefully received


title: Should I bin my continuous input variables for neural networks?

text: And what is the optimal method to do this? Is there a package for this in R?


title: Random effects in mixed models

text: I am not very familiar with mixed models --- I have been reviewing the various tutorials on the web but still am not sure how to specify my model, even as a starting point for comparison. I have a relatively simple design:   * two non-random groups (normal/obese) -- randomized within each group to two conditions (control/treatment). And the DV measured at pre and post-test. So I am interested in looking at differences between pre and post, for the treatment vs. control and also between the two groups. I am using SPSS -- this is the syntax I've come up with so far.               MIXED Score BY Group Treatment Time     /CRITERIA=CIN(95) MXITER(200) MXSTEP(10) SCORING(5) SINGULAR(0.000000000001) HCONVERGE(0,     ABSOLUTE) LCONVERGE(0, ABSOLUTE) PCONVERGE(0.000001, ABSOLUTE)     /FIXED=Group Treatment Group*Treatment Time Group*Time Treatment*Time Group*Treatment*Time | SSTYPE(3)     /METHOD=REML     /PRINT=G R SOLUTION TESTCOV     /RANDOM=INTERCEPT | SUBJECT(ID) COVTYPE(VC)     /REPEATED=Time | SUBJECT(ID) COVTYPE(UN)      Some of my questions --   * I have specified a random intercept for subject ID, but do I need to account somehow that these are random effects from TWO separate populations (normal weight and obese)? If so, how would I do this?   * It seems that I could estimate the time effect using either the `REPEATED` or `RANDOM` command (by adding time to the `RANDOM`). I can't figure out the difference of what that means conceptually and how to decide which fits my data or design better.   * Is it correct to list `TIME` as a fixed factor? I have also seen in one of the books I'm looking at that variables have been coded as 0/1 indicators and entered as covariates. I have no idea why this was done -- I can't seem to find an explanation. However, If I do this, it changes my fixed effects (even with all the same effects and interactions specified). And changes on whether the pre or post is specified as 0/1. So how does mixed models treat factors vs covariates differently? What could be the advantage of treating a variable as a covariate rather than a factor? I noted that although the tests of fixed effects changed, the estimates using `EMMEANS` or `TEST` remained the same, which I also thought was weird.   * Does one need to account for baseline group differences, and if so how would this be done? Since my data is in long format I no longer have a pre-test measurement that I can put in as a covariate in post-test scores. Sorry if these are dumb questions -- I am much more used to ANCOVA and GLM.


title: Does control get us closer to or farther from causation?

text: In logistic regression with an N of 40,000, purchase decision is unrelated to price. However, with certain demographic variables controlled, price can show a positive coefficient of meaningful strength. (Searches for 2-, 3-, and 4-way interactions involving price have yielded nothing.) The goal here is to maximize the leverage one can get out of discounting. There seems to be a real risk of confusing the "real" relationship with one that, in the words of Elazar Pedhazur, has an "air of fantasy about it" due to statistical manipulation (regression control). Any suggestions as to how to proceed?


title: Perfectly correlated (normal) random variables

text: I am not sure in the terminology, so I will simply try to explain the situation that I would like to model as I see it. Suppose there is a set of random variables. The variables are correlated in such a way that they deviate from their expected values into the same direction all together. By this I mean that they can be either all together larger then their expectation or all together lower. Is it possible to model such a dependency with a multivariate normal random variable $\mathbf{X} \sim \mathcal{N}(\mathbf{\mu}, \mathbf{\Sigma})$, assuming the knowledge about the marginal distributions of the components $\mathbf{X}_i \sim \mathcal{N}(\mu_i, \sigma^2_i)$? How to construct $\mathbf{\Sigma}$ is this situation? Thank you. Best wishes, Ivan


title: Interpretation and validation of a Cox proportional hazards regression model using R in plain English

text: Can someone explain my Cox model to me in plain English? I fitted the following Cox regression model to **all** of my data using the `cph` function. My data are saved in an object called `Data`. The variables `w`, `x`, and `y` are continuous; `z` is a factor of two levels. Time is measured in months. Some of my patients are missing data for variable `z` ( _NB_ : I have duly noted Dr. Harrell's suggestion, below, that I impute these values so as to avoid biasing my model, and will do so in the future).               > fit <- cph(formula = Surv(time, event) ~ w + x + y + z, data = Data, x = T, y = T, surv = T, time.inc = 12)          Cox Proportional Hazards Model     Frequencies of Missing Values Due to Each Variable     Surv(time, event)    w    x    y    z                      0    0    0    0   14                           Model Tests          Discrimination                                                  Indexes             Obs       152   LR chi2      8.33    R2       0.054         Events     64   d.f.            4    g        0.437         Center 0.7261   Pr(> chi2) 0.0803    gr       1.548                         Score chi2   8.07                                           Pr(> chi2) 0.0891                                                   Coef    S.E.   Wald Z   Pr(>|Z|)              w      -0.0133  0.0503    -0.26     0.7914                x      -0.0388  0.0351    -1.11     0.2679                y      -0.0363  0.0491    -0.74     0.4600                z=1     0.3208  0.2540     1.26     0.2067      I also tried to test the assumption of proportional hazards by using the `cox.zph` command, below, but do not know how to interpret its results. Putting `plot()` around the command gives an error message.                cox.zph(fit, transform="km", global=TRUE)                 rho chisq      p      w      -0.1125 1.312 0.2520      x       0.0402 0.179 0.6725      y       0.2349 4.527 0.0334      z=1     0.0906 0.512 0.4742      GLOBAL      NA 5.558 0.2347      * * * ### First Problem   * Can someone explain the results of the above output to me in plain English? I have a medical background and no formal training in statistics. ### Second Problem   * As suggested by Dr. Harrell, I would like to internally validate my model by performing 100 iterations of 10-fold cross-validation using the `rms` package (from what I understand, this would entail building `100 * 10 = 1000` different models and then asking them to predict the survival times of patients that they had never seen). I tried using the `validate` function, as shown.              > v1 <- validate(fit, method="crossvalidation", B = 10, dxy=T)     > v1           index.orig training    test optimism index.corrected  n     Dxy      -0.2542  -0.2578 -0.1356  -0.1223         -0.1320 10     R2        0.0543   0.0565  0.1372  -0.0806          0.1350 10     Slope     1.0000   1.0000  0.9107   0.0893          0.9107 10     D         0.0122   0.0128  0.0404  -0.0276          0.0397 10     U        -0.0033  -0.0038  0.0873  -0.0911          0.0878 10     Q         0.0155   0.0166 -0.0470   0.0636         -0.0481 10     g         0.4369   0.4424  0.6754  -0.2331          0.6700 10      How do you perform the 100x resampling? I think my above code only performs the cross-validation once.   * I then wanted to know how good my model was at prediction. I tried the following:              > c_index <- abs(v1[1,5])/2 + 0.5     > c_index     [1] 0.565984      Does this mean that my model is only very slightly better than flipping a coin? ### Third Problem Dr. Harrell points out that I have assumed linearity for the covariate effects, and that the number of events in my sample is just barely large enough to fit a reliable model if all covariate effects happen to be linear.   * Does this mean that I should include some sort of interaction term in my model? If so, any advice as to what to put?


title: Standard errors of regression coefficients based on sample size

text: For any particular nonlinear regression: $$Y_i = f(\mathbb{x_i},\theta) + \epsilon_i, i=1,...,n$$ I currently have standard errors for each of the $\theta_j$ obtained via the Gauss-Newton algorithm Is there a systematic relationship between these standard errors and the sample size such that I would be able to predict the standard error if I had a different sample size? Secondly, forgetting nonlinear regression, can this be done easily for linear regression? Is it just proportional to $\frac{1}{\sqrt{n}}$?


title: Which variables are driving correlations within groups

text: I'm running an analysis on a few data sets that each typically have 100-200 cases measured across 120-160 variables - something similar to looking at gene expressions. Each variable is a non-centered score for expression of a particular attribute frequency for each case. In many cases though, any given attribute is likely to be 0 (i.e. very sparse for most). The cases typically fall into 2-4 natural groups, and I'm trying to figure out how to find out which high-expression attributes are most representative for each group and/or which ones are driving the distinctions between group memberships. I've been experimenting with using correlation clustering and the resulting groups _do_ appear to match the natural groups assigned by other means, so now I'm just looking for a way to "unpack" those clusters in terms of their attribute expressions to find out which attributes/variables are the most influential ones in each of my samples. Given the number of variables involved, it seems like the usual approaches like PCA or discriminate or factor analysis would be very cumbersome. So far, I haven't really found much information on how to deal with variable influence that would fit situations when there are scores or hundreds of variables. Any suggestions?


title: Is the exponential distribution a good model for this data?

text: I'm trying to determine if the exponential distribution is a good model for a data set that I'm exploring. It doesn't have to be precise. I'm using the data for capacity planning (if it's a good fit) and for my own learning. I found a data set that I have ready access to and I wanted to test to see it fit a known distribution. To test if it is a good fit, I've generated 2 plots from a data set with about 1.5M entries. In the first plot, I take the data and generate an empirical CDF and plot the log of the CCDF and the linear model using `geom_smooth`. I'm doing this because if it is exponential (y ~= e^−λx), then the log of both sides gives log y ~= -λx which should be a straight line with slope -λ.               data <- sort(data)      valid <- data[1:(length(data)*0.96)]       ecdf <- (1:length(valid))/length(valid)     df <- data.frame(cbind(x=valid, y=ecdf))            # if this is exponential, get the slope and the x intercept     lm(log(1-y)~x, data=df, subset=seq(1,length(valid)-1))          ggplot(aes(x=x, y=log(1-y)), data=df[1:length(valid)-1,]) +         geom_line(color="blue") +         geom_smooth(method="lm", color="red", size=0.3) +         labs(x="Data", y="CCDF of Data Log(y)", title="Line of Best Fit on Log(CCDF) of Data")      Here is the result. It mostly fits except at the tail. Should I conclude that the exponential distribution is a (mostly) good fit? ![enter image description here](http://i.stack.imgur.com/g1w92.png) Edit: here is the summary data from the model generated by `lm` for the Log(CCDF) plot:               Residuals:         Min      1Q  Median      3Q     Max      -9.7550 -0.0750  0.0044  0.1171  0.1328           Coefficients:                   Estimate Std. Error t value Pr(>|t|)         (Intercept)  1.239e-01  2.388e-04   518.8   <2e-16 ***     x           -1.868e-03  2.989e-07 -6248.6   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.1907 on 1473365 degrees of freedom     Multiple R-squared: 0.9636,     Adjusted R-squared: 0.9636      F-statistic: 3.904e+07 on 1 and 1473365 DF,  p-value: < 2.2e-16       Next I attempted to generate a Q-Q plot by generating a data set that is the same size as the empirical data set (about 1.5M) using a function that generates random numbers from an exponential distribution. From the `lm` above I used the x output and switched the sign as my lambda. I sorted both of them and plotted the sampled data as x and the empirical data as y. Here is the result. Again, it's not quite x=y. Should I conclude that the exponential distribution is a good enough fit? ![enter image description here](http://i.stack.imgur.com/sQhk6.png) Edit: here is the summary data for the model from the Q-Q plot as generated by `lm`.               Residuals:         Min      1Q  Median      3Q     Max      -49.256  -9.214  -0.069  10.060 268.769           Coefficients:                   Estimate Std. Error t value Pr(>|t|)         (Intercept) -1.757e+01  4.086e-02  -429.9   <2e-16 ***     x            1.264e+00  6.078e-05 20803.7   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 32.82 on 1488715 degrees of freedom     Multiple R-squared: 0.9966,     Adjusted R-squared: 0.9966      F-statistic: 4.328e+08 on 1 and 1488715 DF,  p-value: < 2.2e-16       Thanks everyone for all of the helpful comments!


title: OK to use residual sum of squares for cross-validation of binary outcome?

text: For an OLS model the mean squared error can be used to assess the fit of the trained model on the validation data. What is the equivalent for a logistic regression model? Can I simply use the following residual sum of squares function? $RSS=\Sigma^N_{i=1}(y_i-{\hat p})^2$


title: Why do R function `eigen()' and Armadillo's `eig_sym()` give different results

text: I am trying to compute eigenvalues in C++ using the `Armadillo` function `eig_sym` via RcppArmadillo. The results are not entirely the same as the output of the R function `eigen()`: In R:               set.seed(1)     X=matrix(sample(1:25), 5)     X          #     [,1] [,2] [,3] [,4] [,5]     #[1,]    7   18    4   19    6     #[2,]    9   22    3   25   16     #[3,]   14   12   24    8    2     #[4,]   20   11   21   23   15     #[5,]    5    1   13   10   17          Xcov=cov(X)     eigen(Xcov)$values             #[1]  1.585160e+02  7.475128e+01  5.938207e+01  3.250609e+00 -4.293203e-15      In C++: //--------------------------C++ Code ------------------------               #include <RcppArmadillo.h>     // [[Rcpp::depends(RcppArmadillo)]]          using namespace Rcpp;     using namespace arma;          // [[Rcpp::export]]     arma::vec eigenval(arma::mat M) {         arma::vec values=arma::eig_sym(M);         return values;     }      I then compile the code in R using `sourceCpp()`:               sourceCpp("eigenval.cpp")     sort(as.vector(eigenval(Xcov)), TRUE)     #1.585160e+02 7.475128e+01 5.938207e+01 3.250609e+00 1.065789e-14      Compare the above results with the R function result (repeated below), we see that the last value is different.               #[1]  1.585160e+02  7.475128e+01  5.938207e+01  3.250609e+00 -4.293203e-15      Not sure if this is the right place to ask, but I wonder if anyone has any idea about the difference.


title: Fitting a simple JAGS model with RSTAN

text: I'm trying to fit a simple exponential model for left censored data using RSTAN to replicate something I did in JAGS. The JAGS model is:               model{       for(i in 1:N){         isAboveLOD[i] ~ dinterval(x[i], LOD[i])         x[i] ~ dexp(lambda)       }            lambda ~ dgamma(0.001, 0.001)     }      A full explanation of the simulated data can be found at http://jamescurran.co.nz/2014/06/bayesian-modelling-of-left-censored-data- using-jags/, but the essence of it is that the data had 9691 censored values (out of 10,000) and that the data was originally sampled from an Exp(1.05) distribution, with left censoring performed at log(29). When I run this model with RSTAN, I essentially get no movement in the chain, and hence the estimate of lambda is way off (0.3 in my first run). I'm sure the error is elementary, but I'd like your advice. My RSTAN code (including data simulation) is given below.               expCode = '       data {         int<lower=0> nObs;         int<lower=0> nCens;         real<lower=0> yObs[nObs];         real<lower=0> U;       }       parameters {         real<lower=0, upper=U> yCens[nCens];         real<lower=0> lambda;       }       model {         for (n in 1:nObs)           yObs[n] ~ exponential(lambda);         for (n in 1:nCens)           yCens[n] ~ exponential(lambda);         lambda ~ gamma(0.001, 0.001);       }       '               set.seed(35202)     x = rexp(10000, rate = 1.05)          ## set all the censored values to NA's     x[x < log(29)] = NA          yObs = x[!is.na(x)]     yCens = x[is.na(x)]          dataList = list(nObs = length(yObs), nCens = length(yCens), U = log(29), yObs = yObs)     library(rstan)     fit = stan(model_code = expCode, data = dataList, iter = 1000, chains = 1)


title: About cross-validation for machine learning

text: Assume I have 1000 samples of data. I split the data randomly into training and test sets of size 800 and 200, respectively. Now, I train a classifier using the training set, and then evaluate the performance of the classifier using the test set. Assume I am not satisfied with the results. So, I change some parameters, and do this again. Is the approach wrong? Should I have split the 800 samples into another two sets, fine-tune model, and only then try it on the test set? What if my sample is too small for two splits? The data is split randomly, so I assume this would not be completely wrong?


title: Is there a general definition of the effect size?

text: The `effect-size` tag has no wiki. The wikipedia page about the effect size does not provide a precise general definition. And I have never seen a general definition of the _effect size_. However when reading some discussions such as this one I am under the impression that people have in mind a general notion of effect size, _in the context of statistical tests_. I have already seen that the standardized mean $\theta=\mu/\sigma$ is termed as the effect size for a normal model ${\cal N}(\mu,\sigma^2)$ as well as the standardized mean difference $\theta=(\mu_1-\mu_2)/\sigma$ for a "two Gaussian means" model. But how about a general definition ? The interesting property shared by the two examples above is that, as far as I can see, _the power depends on the parameters only through $\theta$ and is an increasing function of $|\theta|$_ when we consider the usual tests for $H_0:\{\mu=0\}$ in the first case and $H_0:\{\mu_1=\mu_2\}$ in the second case. Is this property the underlying idea behind the notion of effect size ? That would mean that the effect size is defined up to a monotone one-to-one transformation ? Or is there a more precise general definition ?


title: Statistical Distance Within Multivariate Distribution

text: Suppose you have some multivariate normal distribution $X \sim N(\mu,\Sigma)$. Is there a good way to calculate a measure of distance between an arbitrary $X_{i}$ and $X_{j}$ of $X$ (I suppose the simplest approach would be to begin with a two dimensional distribution, which can be obtained through the conditional multivariate normal for any arbitrary large distribution)? Beyond just correlation, I want to be able to take into account the differences in the means and variances as well when thinking about how similar or different they are. Some approaches I'm familiar with don't seem to be what I'm looking for (Mahalanobis distance can provide a distance for a random sample from the distribution $X$ and the Kullback–Leibler divergence or Bhattacharyya distance would allow you to calculate how different $X$ would be from some other distribution $Y$).


title: Count data with continuous variables

text: Let's assume one has looked at the occurence of a plant species at different altitude, different temperature and environment. Here are its data:                              Temperature    environment   altitude         plant1        18.1           mud         812         plant2        15.3          field        754         plant3        17.4           mud         213         plant4        15.2          forest       678         plant5        16.6          field        1023         etc...      This guy wants to know if the abondance of plants (number of plants) that one can find is dependent on the three variables. He could run a GLM with poisson distributed errors for the `environment` variable with these data:               Number of plants     environment          311               forest          102                mud           71               field          etc...      How could he do in order to evaluate the effect of continuous variables on the abondance of plants without having to cut the variable into chunks (`Temperature`: [10.1:13], [13.1:16], [16.1:19])? What type of model can incorporates continuous and categorial variables in order to makes sense of count data? Would it make sense to run two Kolmogorov-Smirnov tests on `Temperature` and `altitude` to check if these variables are uniformly distributed and a test GLM, poisson distributed errors for `environment`? But then what happen if it tends to have more `mud` in high altitude? He would need one model with all variables.


title: Model for survival data with an additional response

text: I have a dataset looking like:               > dat       time censoring  x     1   20         0 NA     2   10         0 NA     3   60         0 NA     4   80         0 NA     5  100         1  5     6   10         0 NA     7  100         1  6      (Of course the real dataset is bigger). The "time" column is the time-to-death of the individuals. For those indivduals who are still alive, we also have another response variable "x". This variable "x" is not a time-to-event, this is a quantitative variable like an indicator of the severity of a disease. My question is: what statistical model could we use for such a dataset ? I want to jointly model all the columns of the dataset. Moreover I have in fact such a dataset for two groups of individuals, and the goal is to compare the two groups.


title: Using days and months as levels in hierarchy

text: I have two different sources for order / purchase related data, one is at individual (purchase) level, the other is aggregated daily. Both sources contain different _types_ of data points, albeit all related to purchases. I want to use this data to model changes in purchase patterns. My question is can I use individual vs daily records as seperate levels? Does it make sense to use time itself as a level in a hierarchy? Until now I've mostly seen schools, states, cities, etc. being used as levels, I never saw examples using time. The alternative, I guess, is to summarize points in the specific source so it becomes daily, but I am guessing I could be losing some valueable information that way. Any ideas would be appreciated,


title: Why do I get wildly different results for poly(raw=T) vs. poly()?

text: I want to model two different time variables, some of which are heavily collinear in my data (age + cohort = period). Doing this I ran into some trouble with `lmer` and and interactions of `poly()`, but it's probably not limited to `lmer`, I got the same results with `nlme` IIRC. Obviously, my understanding of what the poly() function does is lacking. I understand what `poly(x,d,raw=T)` does and I thought without `raw=T` it makes orthogonal polynomials (I can't say I really understand what that means), which makes fitting easier, but doesn't let you interpret the coefficients directly.   I read that because I'm using the predict function, the predictions should be the same. But they aren't, even when the models converge normally. I'm using centered variables and I first thought that maybe the orthogonal polynomial leads to higher fixed effect correlation with the collinear interaction term, but it seems comparable. I've pasted two model summaries over here. These plots hopefully illustrate the extent of the difference. I used the predict-function which is only available in the dev. version of lme4 (heard about it here), but the fixed effects are the same in the CRAN version (and they also seem off by themselves, e.g. ~ 5 for the interaction when my DV has a range of 0-4). The lmer call was               cohort2_age =lmer(churchattendance ~      poly(cohort_c,2,raw=T) * age_c +      ctd_c + dropoutalive + obs_c + (1+ age_c |PERSNR), data=long.kg)      The prediction was fixed effects only, on fake data (all other predictors=0) where I marked the range present in the original data as extrapolation=F.               predict(cohort2_age,REform=NA,newdata=cohort.moderates.age)      I can provide more context if need be (I didn't manage to produce a reproducible example easily, but can of course try harder), but I think this is a more basic plea: explain the `poly()` function to me, pretty please. ### Raw polynomials ![Raw polynomials](http://i.stack.imgur.com/2FIIO.jpg) ### Orthogonal polynomials (clipped, nonclipped at Imgur) ![Orthogonal polynomials](http://i.stack.imgur.com/CvVdC.jpg)


title: Mean$\pm$SD or Median$\pm$MAD to summarise a highly skewed variable?

text: I'm working on highly skewed data, so I'm using the median instead of the mean to summarise the central tendency. I'd like to have a measure of dispersion While I often see people reporting **mean $\pm$ standard deviation** or **median$\pm$quartiles** to summarise the central tendency, is it ok to report **median $\pm$median absolute dispersion (MAD)**? Are there potential issues with this approach? I would find this approach more compact and intuitive than reporting lower and upper quartiles, especially in large tables full of figures. Thanks for any hint! Mulone


title: Regression technique for data comprised of categorical explanatory variables & a continuous response variable

text: i suppose one way to characterize data is by a combination of the variable types that comprises it:                Continuous/Continuous  |  Continuous/Discrete      -----------------------|---------------------      Discrete/Discrete      |  Discrete/Continuous      Each of the four cells comprising the 2x2 table just above is comprised of two descriptors, one for the Explanatory Variables (EV), which we'll assume are all of the same type, and one for the Response Variable (RV):               Explanatory Variable Type / Response Variable Type      The columns represent the EV type; the rows, the CV type. Nearly all of the data i see can be placed in either of the two cells that comprise the _first row_. So for instance, _OLS Regression_ is a suitable model type for data in row1/col1; and for data in row1/col2, _Logistic Regression_ is an appropriate model choice. It's the second row, and in particular _row2/col2_ , that my Question is directed to. I'm aware of a few regression techniques like _ordinal regression_ which handle a particular type of discrete variables (ordinary or rank, 1st, 2nd, 3rd,....) but i am interested in techniques for handling discrete data more generally, given that most categorical variables do not have an implicit ordinal relationship among the values that comprise them. For instance:               Sex | City_of_Residence | Car_Make&Model | Married? | DUI? | Prior_3P_Claims?      F  |   Cleveland       |  Chevy Camaro  |   No     |  No  |   Yes      And the Response Variable is _continuous_ --e.g., building a model to predict the quotes offered by major auto insurers-- _a price_.


title: ARIMA modeling with non-stationary daily data, at what lag do I difference the series?

text: If I have daily data set and it's a non-stationary series, then what is the lag I have to consider for the first difference 1 or 7?


title: How can I compute a repeated measurement ANOVA with a mixed model?

text: It is my task to program a mixed model. The data is from 20 experimental subjects, who were measured under 2 factors. Factor A has 7 levels, factor B has 3 levels. So, that for every subject I got 21 data points. The data points for one subject are repeated. I analysed the data with repeated measurement ANOVA just to know the influence of factor A and B. Now I want to do my task. But I don't know what the model looks like. Because I have a random factor and two fixed effects. I wanted to compute the model about the variance components.   * Do I have to use the formula for a three-way ANOVA with mixed effects?   * How would I compute it in a regression model?   * Does my model look OK? $Y = Z * \beta + X * b$   $Y$ ... data points   $Z$ ... design matrix for random effects (20 x 20)   $X$ ... design matrix for fixed effects (20 x 21)   **EDIT:** Yes, I have 21 data points. I used ANOVA with complete repeated measurement (also known as ANOVA with one observation per cell).


title: Difference between arithmetic vs geometric random walk

text: I have read about arithmetic and geometric random walks. What is the difference between them?


title: Multivariate Interpolation Approaches

text: Is there a good, modern treatment covering the various methods of multivariate interpolation, including which methodologies are typically best for particular types of problems? I'm interested in a solid statistical treatment including error estimates under various model assumptions. An example: Shepard's method Say we're sampling from a multivariate normal distribution with unknown parameters. What can we say about the standard error of the interpolated estimates? I was hoping for a pointer to a general survey addressing similar questions for the various types of multivariate interpolations in common use.


title: Variance Estimation in case of nonrespondents

text: I saw in the book of Rubin (1987) that an increase in variance of estimation will occur due to nonresponses. But I wonder the reason behind this. Thanks for your shares!


title: Justifying unsupervised clustering using Random Forest?

text: I have been looking at ways to carry out unsupervised clustering of data with both numeric and nominal (but not ordinal) variables. I also suspect non- linearity in the data. A possible solution would appear to be as follows: (1) use Random Forest in unsupervised mode to get a matrix of pairwise similarities (or 'proximities') followed by (2) applying your favourite clustering algorithm to this matrix. This approach was advocated by Tao Shi, Steve Horvath and colleagues: http://labs.genetics.ucla.edu/horvath/RFclustering/RFclustering.htm (See their technical report or corresponding journal articles for a detailed description of the methodology.) However, what I find confusing about the proposed approach is that (as I understand it) it works via separating the original data from a permuted version of the data (with variable values randomly permuted between instances) such that the significant variables which contribute most to determining the proximities are those which are highly correlated with other variables in the original data. (I have had a cursory look at this using the randomForest implementation in R and this seems consistent with my understanding.) Have I misunderstood something? If not, how can this be justified? In the area I work in (modelling of chemical data using quantitative-structure activity relationships), correlation between predictor variables (descriptors) is usually taken as indicative of (some degree of) redundancy and the discarding of (one of a pair of) correlated descriptors is advocated. I can appreciate that discarding non-completely redundant descriptors may loose information, but I don't understand why focusing on correlated descriptors would make sense when trying to determine 'logical' clustering of the data? Does anyone have any justification for this?


title: Vocabulary: do we measure actual values or observations?

text: Consider that $\theta$ is an hidden parameter and one has an observation such that $O$: $$ O \sim N(\theta,\sigma^2). $$ My question concerning vocabulary: do we measure $\theta$ and it gives us $O$? (so we measure the true value) or do we measure $O$ ? (so we measure the observation) I am looking for unquestionable sources.


title: Karlin-Rubin theorem in hypothesis testing

text: In karlin-Rubin theorem, after showing that power function is an increasing function, we take hypothesis H0: Theta=Theta0 vs H1: Theta=Theta1. then how do they show that {sup(thet<= theta0) power function} is = alpha? I am not able to understand this from any of the research paper.


title: what is the intuition behind stationarity condition for AR(p) process?

text: i get that you have to find the roots of the characteristic polynomial but can someone explain the intuition behind the roots must be outside the unit circle? what is a unit circle? before anyone downvotes me saying theres another similar discussion, i have read it and i dont understand what it means for roots to lie outside a unit circle.


title: Inferential methods on large panel data with sparse clusters and rare outcomes

text: Using longitudinal survey data on children using psychotropic medications, we are interested in estimating associations with medication classes, their persistence and adherence (longitudinal exposures), and certain rare outcomes such as suicidal idealization. There are clusters defined at the child level with no cluster having more than three repeated measures, most clusters have only one observation. There is left truncation of child observations relative to their starting specific drugs and significant recall bias in their self- reported "start date". What statistical methods would be appropriate for estimating the relative rates of rare outcomes in this cohort? A couple designs come to mind:   1. Conditional logistic regression with propensity (baseline-risk) matched case-control design: Use time as a stratification factor and estimate a med-class time interaction   2. Nested case-control: similar design, match based on risk at baseline, but use survival methods for time-to-event analysis. The problem with 1 is that I am not sure if adjustment appropriately accounts for confounding due to time. Some kids may be asymptomatic because of their adherence (or, conversely, adhered to meds because they were asymptomatic) so potentially high risk kids are identified as low risk due to successful treatment over time. The problem with 2 is that we eliminate from risk-sets children who experience the event early and their symptoms reduce over time which is very likely. I don't know nor do I think some time based frailty could/would be appropriate for analysis of this data.


title: Proper splitting of data set for Ensemble methods

text: I have 10,000 documents. Each document has a label ($Y$) that is either $0$ or $1$ (the 0-1 split is pretty much 50/50 over my 10,000 documents). Each document has 10 fields. Each field can have any number of words in it. I create 10 feature-spaces by fitting a tf-idf over each field individually over all documents. It looks something like this:                      For f_1:                                      For f_2:                        34,974                                        113,351               ------------------------------                ------------------------------            |                        |   |                |                        |   |            |                        |   |                |                        |   |     10,000 |            X_1         | Y |         10,000 |            X_2         | Y |            |                        |   |                |                        |   |            |                        |   |                |                        |   |            ------------------------------                ------------------------------      On, and on for each field. Each matrix will have 10,000 rows, but a different number of columns. The $Y$ column is always the same. I'm interested in using each of these matrices as the input to a classifier, and using some ensemble of them to predict the labels $Y$. My initial approach was to choose a random 70% of the 10,000 documents and set that as the training set, and then use the other 30% as my predicting set. My plan was to train a logistic regression model on 70% of $X_1$ and then have that model predict the labels of the remaining 30% of $X_1$ to give me $Y_{lr_1}'$. I would use the same 70% and train a random forest, and then have the random forest predict the 30% to give $Y'_{rf_1}$. I would use the same 70%/30% of rows to train/predict a logistic regression and random forest on $X_2$ through $X_{10}$. In the end I would have some matrix:                                                   20                                      -------------------------------------------------------------           |                                                       |   |           |                                                       |   |     3,000 | Y'_lr1  Y'_rf1  Y'_lr2  Y'_rf2  ...  Y'_lr10  Y'_rf10 | Y |           |                                                       |   |           |                                                       |   |           -------------------------------------------------------------      I then trained a final logistic regression to predict the $Y$ from the 20 $Y'$. Is this a normal technique? Should I be training many models, and do a further voting stage to get each $Y'$? Any help is appreciated. Most of the sources I find talk about drawing N rows with replacement, and then merging those models, but in this case, I have many models training on different features. I don't know how much of a difference this makes. I'm using sklearn on Python if that makes any advice easier to relay! Sorry for the length, but I wanted to be detailed. Thanks!


title: Interpreting significant difference when additional data added

text: I have two questions on interpreting the statistical analysis of data collected from production runs to produce an experimental material. We made four production runs, with some changes in each run. Each run produced three pieces of material, and each piece of material was tested in multiple places to measure the property of interest. So, the measurement results look like this:               Run 1 Piece A -- 25 measurements;     Run 1 Piece B -- 25 measurements;     Run 1 Piece C -- 25 measurements;          Run 2 Pieces A, B, and C -- 25 measurements on each          Run 3 Pieces A, B, and C -- 25 measurements on each      In order to determine if there existed a significant variation among the production runs, I performed Duncan's test using three groups of three replicates, taking the average measurements from one piece as one “replicate” (i.e. the value for group "Run 1" and replicate "Piece A" was the average of the 25 measurements). So the **first question** is: is it valid to take the average of the measurements from one piece and call it one replicate? I could say I have 9 groups with 25 replicates each, but my goal is to compare batches and not pieces independently. The Duncan’s test showed Run 2 to be significantly different from the other two runs. To complicate things, we made a fourth run:               Run 4 Piece A -- 9 measurements;     Run 4 Piece B -- 9 measurements;     Run 4 Piece C -- 9 measurements.      The fourth group had much higher averages and standard deviations from the others. When I repeat Duncan’s test adding in this fourth group of data, the test shows the fourth group to be significantly different from the other three, but the other three are no longer different from each other. So, I am faced with interpreting the significant difference. With only three runs Run 2 was significantly different from Runs 1 and 3. After doing the fourth run, Run 2 was no longer different from Runs 1 and 3, just different from Run 4. I understand mathematically what occurred – that the results from Run 4 were so far different that they pushed the first three runs together. But what is the practical interpretation ( **second question** )? Is Run 2 really different from Runs 1 and 3? Any thoughts on this problem are greatly appreciated!


title: Can First Differencing Cause Negative Serial Correlation

text: Ex. series, say stock prices 103 101 102 **150** 101 102 100 First differenced 2 1 **48** **-49** 1 -2 Notice you could guess a very large negative number following the very large positive in the first differenced numbers. Is it the case that negative serial correlation is common in first differenced data and do you correct for this using MA models? (including the lagged error term?) I couldn't find any hints about this on a search so perhaps something is wrong in my logic and I need a vacation..


title: How to restore market making algorithm?

text: I need your advice to know how to restore certain market maker algorithm based on historical data. I have dataset with limit order book information at discrete time periods. For every time moment it contains best bid/ask prices and volume of fist 10 asks and 10 bids. Also it has bid and ask prices - the choice of the market maker. It's known that market maker's decision at time $t$ depends on order book condition at $t$ and maybe previous time moments. I tried to use machine learning regression methods like reg.tree ensembles, neural nets and SVR to educate a model based on given features but when I tried to check prediction quality on cross-validation it showed very bad quality. I think that I should properly extract features from my data but I don't know how to do it. Or maybe should I use different model and not regression algorithms?


title: How to reduce the dimension of $10^8$ vectors

text: I have $10^8$ vectors in $1000$ dimensions each. I would like to drastically reduce their dimensions. However PCA seems computationally infeasible. Are there near linear time methods to do dimensionality reduction?


title: Clustering binary categorical data

text: I have some data where I have certain classes (c1, c2, c3, c4 ...) and the data comprises of binary vectors where 1 and 0 denote that an entry belongs to a class or not. The number of classes will be > 200. The data will look like this:               c1    c2    c3    c4    ...      1     0     0     0    ...      0     1     1     0      Would this data come under "Categorical" type? Details:   * **Sample Size** : ~20000   * **No. of classes** : 300   * **Data Matrix Sparsity** : 99.52%   * **Problem Statement** : The classes that I am talking about are medical services provided by Hospitals. If a hospital provides the service we just put 1 or else 0 in the binary vector. I want to cluster similar hospitals on the basis of their services. I tried out PCA for dimension reduction on this dataset and I even got good clusters with DBSCAN but I read that for categorical sparse data PCA is not recommended and also Euclidian distance as the distance measure is not good. I am planning to use MCA (Multiple Correspondence Analysis) but I cannot figure out how am I supposed to represent the data for that.


title: How to determine whether a data set is continuously growing?

text: I sampled some data at certain sampling rate (e.g., 1 sample per minute). I want to know whether the values in this data set are continuously growing or they are just fluctuating around certain value. How can I do that? In addition, I would like to predict the future values to be sampled. Based on the prediction (whether it will grow along with certain confidence level), I can decide whether to continue the experiment to sample more data or simply stop the experiment. **Update** : I think this is exactly what I need: Trend Esitmation **Note:** I know very little (if not nothing) about statistics, so bear with me and do correct me if you find me using wrong terminology.


title: Analyzing repeated measures ANOVA with two groups

text: I have a data set that looks like this:               GROUP                    WEEK1 WEIGHT(kg)          WEEK2 WEIGHT(kg)         WEEK3 WEIGHT(kg)     Control       Subject C1             40                        42                       44       Subject C2             52                        57                       57        Subject C3             55                        55                       60       Subject C4             45                        47                       49          Experimental       Subject E1             38                        40                       49       Subject E2             42                        47                       50        Subject E3             55                        60                       68       Subject E4             65                        77                       88      My question is: What tool shall I use to test my hypothesis that the two groups (control and experimental) do not significantly differ with regard to weight? Please note that the weights are repeated measures. I understand that if I were to examine the data per group (`Control` first, then `Experimental`), I'd perform a repeated measures ANOVA for control and another for experimental.


title: P-value approximation from an underlying unknown distribution

text: I am calculating a certain statistic(relative entropy) of a certain sequence alignment over DNA sequences. The statistic has a certain distribution in the entire genome of which I only have 1000 sequences of length 1000. I want to test the significance of a certain value of the statistic but I obviously do not know the distribution of the statistic so I can't use t-test and such and I also have a sample of size 1000 from the population. I guess I could use non-parametric methods like bootstrapping and such to approximate the p-value of a certain statistic value but the research papers that I have been reading seem to use methods from Large Deviation Statistics and Fast Fourier Transforms. I have not have time to learn these concepts on my own and just need a nudge in the right direction. Any ideas ?


title: Measuring correlation between two matrices

text: I have two matrices A and B of the same dimensions. Each cell a(i,j) represents some parameter I'm interested in, and b(i,j) represents the corresponding value in the other matrix. The cells are not interrelated. What's the best way to represent how close these two matrices are? The entries in these matrices are the number of occurrences of some event out of a large sample of events. If this is the case, does the problem make more sense?


title: Sequential clustering with an unknown number of clusters

text: Is there any clustering algorithms that:   * does not assume the number of clusters to be known, and   * processes the data in only one pass by considering it as a continuously arriving data stream (and we do not know the data set size beforehand)?


title: Alternative to spherical K-Means for clustering large high dimensional dataset

text: What are some alternatives to Spherical K-Means for clustering very large datasets of high dimension? I'm looking for something that will be fast even on large datasets, and preferably will not require you to specify the number of clusters to find. Another property that would be nice (I'm not quite sure the name of it) is that the algorithm can calculate the clusters of a dataset A, and then later find the best cluster for a new point b which wasn't in A, without having to recalculate the clustering.


title: Measuring Regression to the Mean in Hitting Home Runs

text: Anyone that follows baseball has likely heard about the out-of-nowhere MVP- type performance of Toronto's Jose Bautista. In the four years previous, he hit roughly 15 home runs per season. Last year he hit 54, a number surpassed by only 12 players in baseball history. In 2010 he was paid 2.4 million and he's asking the team for 10.5 million for 2011. They're offering 7.6 million. If he can repeat that in 2011, he'll be easily worth either amount. But what are the odds of him repeating? How hard can we expect him to regress to the mean? How much of his performance can we expect was due to chance? What can we expect his regression-to-the-mean adjusted 2010 totals to be? How do I work it out? I've been playing around with the Lahman Baseball Database and squeezed out a query that returns home run totals for all players in the previous five seasons who've had at least 50 at-bats per season. The table looks like this (notice Jose Bautista in row 10)                    first     last hr_2006 hr_2007 hr_2008 hr_2009 hr_2010     1    Bobby    Abreu      15      16      20      15      20     2   Garret Anderson      17      16      15      13       2     3  Bronson   Arroyo       2       1       1       0       1     4  Garrett   Atkins      29      25      21       9       1     5     Brad   Ausmus       2       3       3       1       0     6     Jeff    Baker       5       4      12       4       4     7      Rod  Barajas      11       4      11      19      17     8     Josh     Bard       9       5       1       6       3     9    Jason Bartlett       2       5       1      14       4     10    Jose Bautista      16      15      15      13      54      and the full result (232 rows) is available here. I really don't know where to start. Can anyone point me in the right direction? Some relevant theory, and R commands would be especially helpful. Thanks kindly Tommy Note: The example is a little contrived. Home runs definitely aren't the best indicator of a player's worth, and home run totals don't consider the varying number of chances per season that a batter has the chance to hit home runs (plate appearances). Nor does it reflect that some players play in more favourable stadiums, and that league average home runs change year over year. Etc. Etc. If I can grasp the theory behind accounting for regression to the mean, I can use it on more suitable measures than HRs.


title: How to paste SPSS output using the LIST command into SPSS - the output all pastes into one cell

text: If I use select if and then list some variables, the output appears as a plain text table in the output viewer. When I paste this into Excel it all goes into one cell. I've tried the copy special options and none of them seem to make a difference. Also the columns of data are all separated by different numbers of spaces, which means I can't do a find and replace to switch the spaces for commas or tabs. How can I get this output so that I can paste it into excel with each variable in its own column?


title: PIT on a sample with m bins, and KS test used to estimate a good value for m

text: I know about PIT, but this works only when you know the distribution, or at least have a strong hint. What I am trying to achieve is to transform a given sample into an equivalent sample with continuous standard uniform distribution. I have a sample of size $n$. I choose arbitrary a value $m$, and estimate $m+1$ quantiles (for example, if $m=4$, I compute quantiles for $\{0, .25, .5, .75, 1\}$). The procedure is described in Wikipedia. Using quantiles I transform each $x_i$. If $x_i$ happens to be exactly a computed quantile, than I know precisely its equivalent value, otherwise I interpolate the equivalent value linearly. I've done a small simulation. I build a random sample from standard normal with $10^6$ values. I applied the described transformation and make KS test for some values of $m$. The results looks like:                 m         D   p-value     100  0.006090   0.000000000000 ***     200  0.003151   0.000000004733 ***     300  0.001991   0.000720875707 ***     400  0.001484   0.024403417075 *     500  0.001057   0.213437843144      It looks like I can do a PIT on sample with only 500 interpolation points. The question is: can I use Kolmogorov-Smirnov one sample goodness-of-fit to find a proper value for $m$ (the number computed quantiles)?


title: processing a sound file for analysis of different spoken languages

text: So i have sound files for 5 languages by 2 person, thus my input data has 10 sound files. Now i want to cluster them based on the languages (thus 5 clusters) and not based on the speaker/voice ( which in case would create only two clusters). but my problem is how do i separate the "language information" from the "speaker characteristics/speaker voice". Is there something in "pca" to solve this problem ? For example if theres a software/method which would convert/change the voice tone of both the speakers in same tone, then that way we can differentiate only on the basis of the language (and not the voice of speakers)


title: calculating distribution of $U=-\ln Y$

text: suppose $X_1,\ldots,X_n$ be random sample of continuous distribution with distribution function $F$ if $$Y=\prod_{i=1}^k F(X_i) \prod_{i=k+1}^n [1-F(X_i)]$$ how can I calculate distribution of $U=-\ln Y$


title: Linearly dependent features

text: I have a matrix A of 1000 observations (rows) and 100 features (cols). I would like to find:   1. Linearly dependent features so that I can remove them and simplify the problem. rank(A) gives me 88, which I assume means that 12 of the features are linearly dependent. Am I right?   2. After the above step, how do I determine which 12 out of the 100 columns are linearly dependent? I know there is no unique answer. But does that mean I can choose any 12 columns?   3. Let's say I choose to remove the last 12 columns. But before removing them, I what to find the 12 linear combinations that compute to the last 12 columns. How do I get these? So far I have tried using Matlab's PCA, QR and SVD, but each of them give different matrices and I don't know how to use these matrices to get what I want.


title: Estimating parameters for univariate skew t

text: How can I solve the MLE for the skew-t distribution via EM? I am comfortable with the EM methods for t, so could someone show it for the skew-t?


title: How to de-trend the data before doing MDS?

text: I am given a set of regions' macro-economical development indices, such as number of unemployed, gross value added, and so on. Each index naturally is correlated with a size of the region. I want to visualize differences between the regions preferably with multidimensional scaling. Before I start, I need to transform each index to remove the dependency on the size of the region, otherwise the results would be quite one-dimensional. Obvious way to do it is to divide each index by the "size of the region". The problem is that the there is no consensus about what the "size of the region" really is: it could be an area, population, total income, or something in between. * * * How I plan to do it: I am thinking about doing Principal Component Analysis on the data (possible log-transformed) to remove the first component (which should be a "size of the region"), and doing the metric MDS on the remaining components/factors. Does my attempt it makes sense?


title: How to calculate the percentage of a sample that is delayed

text: I'm working on a thesis and I'm using a bit statistics to process the data. I have a certain population and in one sample (N=250) I measure the time it takes to reach a goal. Then I change something in the population that delays that time. I take another sample (N=250) and compare those two. I use a T test to get the difference between the averages of the two samples. So I know how much the average delay is. Now I'd also like to know how much percentage of the second sample has a delay. Is it possible to calculate this? And how? EDIT We are simulating a swarm of robots. One of them, the searcher, has to find a target by using the information of the surrounding robots. We do this 250 times and record the time it takes for the searcher to get to the target. This is our reference data. Then we add an enemy who sabotages the searcher and delays the time. We again run this 250 times and record the data. We are sure that some times the searcher gets to the target before the enemy can sabotage it because of the algorithm.


title: Poker shared cards change odds

text: I wrote a poker simulator for seven card texas hold 'em, and I found a counter-intuitive result: If player A starts with a pair of aces against player B who has a random hand, I get a different %win if they share five cards on the board (as it works in Texas hold'em) than if they each have separately drawn hands of seven cards. In the standard Texas Hold'em procedure player A has an 85% chance of winning, but if I draw two seven card hands separately (instead of sharing five cards) then player A has a 77% chance of winning. In both scenarios the cards are drawn from a deck of 50 cards (ie a normal deck excluding the two Ace's I've given player A at the beginning). Between the two scenarios the probabilities I measure of player A and player B having a high card, a pair, two pairs etc are unchanged .. but somehow player A wins more often in one scenario than the other. I've tested lots of times in different ways, and I'm convinced it's not due to a bug in my program. Can someone help me with an intuitive picture of what is going on?


title: What are my choices for transforming data that are not normally distributed?

text: I have two sets of data comparing two variables: x and y. I'd like to transform and compare both sets either in `r` or `matlab`. **Scatter plot** ![enter image description here](http://i.stack.imgur.com/8sB1P.png) Neither are normal and come from (by the looks of it) different distributions) **QQPLOT of data 1** ![enter image description here](http://i.stack.imgur.com/1uGg2.png) **QQPLOT of data 2** ![enter image description here](http://i.stack.imgur.com/y9eCK.png) I would like to transform both data sets to compare them on a level playing field but I know they are not normally distributed, though I don't know from what distribution they come. What are my choices?


title: Output of hyperbfit?

text: I want to fit a hyperbolid distribution to my data, in my notation, I have the density \begin{align*} H(l;\alpha,\beta,\mu,\delta)&=\frac{\sqrt{\alpha^2-\beta^2}}{2\alpha \delta K_1 (\delta\sqrt{\alpha^2-\beta^2})} exp\left(-\alpha\sqrt{\delta^2+(l-\mu)^2}+\beta(l-\mu)\right) \end{align*} with R using HyperbolicDist package and the hyperbFit command I get the output parameter estimates:               pi    zeta    delta     mu     0.09    0.2   0.002   -0.002      can someone explain me, what this parameterization is in my case? So I do not have a pi and zeta?


title: Measuring the performance of a classifier that estimates bin counts for a Poisson process

text: I have some data that results from the convolution of a Poisson random process with a decaying exponential, plus some additive Gaussian noise. I'm working on designing a classifier that returns a score which should be proportional to the number of Poisson hits for each time bin from the noisy data. The data looks like this: ![enter image description here](http://i.stack.imgur.com/AOmfs.png) I'd like to be able to quantify how good my classifier is at predicting the counts for each time bin, given the noisy data. If I treat the true set of bin counts as a boolean vector (i.e. making no distinction between bins with multiple hits and bins with single hits), then it's relatively straightforward to compute an ROC curve for my classifier and quantify the area under the curve (AUC) as my measure of performance: ![enter image description here](http://i.stack.imgur.com/rBcNm.png) However, this gives a misleadingly optimistic measure of performance, since my goal is to estimate the actual bin counts rather than just the presence/absence of any hits. Is there a way to generalize ROC analysis to deal with my case, where there could, in principle, be infinitely many different integer bin counts?


title: Estimation and functional space

text: In the first chapter of the book _Algebraic Geometry and Statistical Learning Theory_ which talks about the convergence of estimations in different functional space, it mentions that the Bayesian estimation corresponds to the Schwartz distribution topology, whereas the maximum likelihood method corresponds to the sup-norm: ![excerpt](http://i.stack.imgur.com/zqeOQ.jpg) Can anyone explain this in a further detail? Thank you.


title: Reducing the number of variables in a multiple regression

text: I have a large data set consisting of the values of several hundred financial variables that could be used in a multiple regression to predict the behavior of an index fund over time. I would like to reduce the number of variables to ten or so while still retaining as much predictive power as possible. _Added: The reduced set of variables needs to be a subset of the original variable set in order to preserve the economic meaning of the original variables. Thus, for example, I shouldn't end up with linear combinations or aggregates of the original variables._ Some (probably naive) thoughts on how to do this:   1. Perform a simple linear regression with each variable and choose the ten with the largest $R^2$ values. Of course, there's no guarantee that the ten best individual variables combined would be the best group of ten.   2. Perform a principal components analysis and try to find the ten original variables with the largest associations with the first few principal axes.  I don't think I can perform a hierarchical regression because the variables aren't really nested. Trying all possible combinations of ten variables is computationally infeasible because there are too many combinations. > Is there a standard approach to tackle this problem of reducing the number > of variables in a multiple regression? It seems like this would be a sufficiently common problem that there would be a standard approach. ~~A very helpful answer would be one that not only mentions a standard method but also gives an overview of how and why it works. Alternatively, if there is no one standard approach but rather multiple ones with different strengths and weaknesses, a very helpful answer would be one that discusses their pros and cons.~~ whuber's comment below indicates that the request in the last paragraph is too broad. Instead, I would accept as a good answer a list of the major approaches, perhaps with a very brief description of each. Once I have the terms I can dig up the details on each myself.


title: Using Fractional Polynomials for Logistic Regression Modelling in R

text: I am learning logistic regression modeling using the book "Applied Logistic Regression" by Hosmer. In chpaters, he suggested using Fractional Polynomials for fitting continuous variable which does not seems to be related to logit in linear fashion. I tried the `mfp` package and can give exactly the same verbose as the book. But I don't know how to write the transformed variable based on the output of fractional polynomials. The book only shows example of the transformed variable when $J=2$ with $p_1=0$ and $p_2=-0.5$ (page 101) and when $J=2$ with $p_1=2$ and $p_2=2$ (page 101), But what about the others? Currently my case is $J=2$ with $p_1=-1$ and $p_2=-1$. I know little about fractional polynomials and the book seems not giving sufficient hits on this part. Can anyone refer me to some place which I can know how to write the polynomial? Thanks.


title: Log transformation of regression

text: I have a model needs to transform to log function to perform a linear model, how can I move from $\log(E(Y))= -1.146 + 0.747x$ to get $\hat y = 0.318(1.0776)^x$?


title: How to calculate t-values from percent ranks?

text: A norm-table from a psychological test is giving percent ranks (PR). How can I transform those into t-values by calculation? (Assuming normal distribution.) I found a lot of references on how to transform z --> t, but none on PR --> t. Way back in university, we had tables with the bell-curve, to look up PR, z-values, t-values and stanines - but how to calculate?


title: What is a "thin" svd?

text: I've seen lots of mention of this in the implementation of SVDs in various programming environments. What does it actually mean?


title: model selection with glmnet

text: I am trying to fit a multinomial logit model using glmnet. I have a few questions:   1. How is the baseline category specified?   2. Looking at the model coefficients using coef.glmnet, I'm thinking that many are given by a dot. I assume this means the coefficient was set to zero by LASSO, so the variable was dropped. However, I'm finding that variables are dropped this way even when lambda=0, so there is no regularizing term. Can someone explain what's going on here?   3. Sometimes I find that a variable is dropped in some but not all of the logits. How should this be interpreted?  Here is a working example:               library(glmnet)     data(iris)     X <- model.matrix(Species ~., data=iris)     y <- iris$Species     fit <- glmnet (X, y, family="multinomial")     coef(fit, s=0)


title: Effective Sample Size for posterior

text: I am trying to implement unsuccessfully a function in matlab, to compute the effective sample size after a MCMC chain, with a posterior with 3 coefficients. ![enter image description here](http://i.stack.imgur.com/2LaS6.png) Source: Sims MCMC $ VAR(1) / Y_t=\mu +AY_{t-1}+ \epsilon_t$ I think that I need to compute $ T\Gamma(0)/\sum_{i=1}^{n}(\Gamma(J)$ If a VAR(1) is stationary $\Gamma(J)=A^j\Gamma(0) $. So the first step is to estimate the model. Second with the matrixes found, I need to compute $\Gamma(0)=V(y_t)=E(y_t'y_t)=A\Gamma(0)A'+\Omega $. Third compute the effective sample size This is my try, a very basic code in Matlab to estimate the VAR(1) coefficients with the posterior load in betaProbit.                   Varprobit=betaProbit; % The posterior betas of x0 x1 x2         Varprobit_1=zeros(ndraws,3);% x0 x1 x2 ,t-1         for j=2:ndraws            Varprobit_1(j,:)=Varprobit(j-1,:);         end         epsilons=randn(ndraws,3);% generate Random normal numbers ( errors)         XVAR1=ones(ndraws,1),Varprobit_1(:,1),Varprobit_1(:,2),...         Varprobit_1(:,3),epsilons(:,1)];         XVAR2=[ones(ndraws,1) ,Varprobit_1(:,1),Varprobit_1(:,2),...         Varprobit_1(:,3),epsilons(:,2)];         XVAR3=[ones(ndraws,1) ,Varprobit_1(:,1),Varprobit_1(:,2),...         Varprobit_1(:,3),epsilons(:,3)];         % Estimate VAR coefficients;         BetasEc1= inv(XVAR1'*XVAR1)*XVAR1'*Varprobit(:,1);         BetasEc2= inv(XVAR2'*XVAR2)*XVAR2'*Varprobit(:,2);         BetasEc3= inv(XVAR3'*XVAR3)*XVAR3'*Varprobit(:,2)   ;          A=[BetasEc1(2:4,1)';BetasEc2(2:4,1)';BetasEc3(2:4,1)'];         OMEGA=eye(3,3);         Gam=Varprobit'*Varprobit;           acfj= VARacfLagj( Gam ,A, OMEGA,j);         %In this step I need to construct a function to solve Gamma(0)      The key issue is to find a function to solve $\Gamma(0)=V(y_t)=E(y_t'y_t)=A\Gamma(0)A'+\Omega $. Thanks for your help.


title: Should I use t-test on highly skewed and discrete data?

text: I have samples from a highly skewed dataset about users' participation (e.g.: number of posts), that have different sizes (but not less than 200) and I want to compare their mean. For that, I'm using two-sample unpaired t-tests(and t-tests with the Welch’s factor, when the samples had different variances). As I have heard that, for really large samples, it doesn't matter that the sample are not normal distributed. My metrics are discrete, they are counts of the number of each user's participation. Of course we have those users who participate much more than the others, but I'm not considering them as outliers. Here are the data description: https://docs.google.com/spreadsheets/d/1WhSKgYIuP35eRsukHVoUFUlITNwO_RRcYoOoR9EmXHg/edit?usp=sharing My problem: someone, reviewing what I've done, said that the tests I am using were not suitable for my data. They suggested to log-transform my samples before using the t-tests. I do know that I can't log-transform these, because all of them have zero- values on the samples. My guess is, if I can't use t-test, I should use the Mann Whitney U test. Are they wrong? Am I wrong? If they are wrong, is there a book or scientific paper which I could cite/show them? If I am wrong, which test should I use?


title: ML algorithm to find optimal control parameter

text: I have a training dataset $(X, y) \rightarrow z$. Where $X$ is an $n$th dimensional real vector, $y$ is an integer number in $\{1, 2, 3\}$, and $z$ is a real number. I am looking for machine learning algorithm which solves the next problem: $X$ is given and we need to choose $y$ for maximization $z$. I would prefer if it would be something like decision trees. I would really appreciate some advice on that.


title: What is the difference between entropy and deviance?

text: In terms of classification task using decision trees, the formula for these looks almost the same. So, how are they different/same? what is the purpose of each in terms of impurity measure? $\text{Entropy}~(p_1,p_2) = -\sum p_i \log (p_i); i= 1,2;$ $p_i $ are fractions. Say, if I have 2 Yes and 3 No in a node, $p_1=2/5$, $p_2=3/5$. $\text{Deviance}~D= - 2\sum n_k \log(p_k) ;~k $ is the class in each leaf. Both are used as impurity measures. But I am not able to understand the difference between these.


title: Does party package in R provide out-of-bag estimates of error for Random Forest models?

text: I'm a new R user, and also new to Random Forest modeling. I cannot seem to figure out how to obtain the out-of-bag (OOB) error estimates for cforest models built with the Party Package in R. In the randomForest package, the OOB error estimates are displayed if you simply "print" the model object, but the party package doesn't work the same way. ### Run random forest model using randomForest package:               > SBrf<- randomForest(formula = factor(SB_Pres) ~ SST + Chla + Dist2Shr + DaylightHours + Bathy + Slope + MoonPhase + factor(Region), data = SBrfImpute, ntree = 500, replace = FALSE, importance = TRUE)     > print(SBrf)          Call:      randomForest(formula = factor(SB_Pres) ~ SST + Chla + Dist2Shr + DaylightHours + Bathy + Slope + MoonPhase + factor(Region),      data = SBrfImpute, ntree = 500, replace = FALSE, importance = TRUE)                     Type of random forest: classification                          Number of trees: 500     No. of variables tried at each split: 2                  OOB estimate of  error rate: 23.67%     Confusion matrix:         0   1 class.error     0 823 127   0.1336842     1 211 267   0.4414226      ### Run random forest model using party package:               > SBcf<- cforest(formula = factor(SB_Pres) ~ SST + Chla + Dist2Shr+ DaylightHours + Bathy + Slope + MoonPhase + factor(Region), data = bll_SB_noNA, control = cforest_unbiased())     > print(SBcf)          Random Forest using Conditional Inference Trees     Number of trees:  500           Response:  factor(SB_Pres)      Inputs:  SST, Chla, Dist2Shr, DaylightHours, Bathy, Slope, MoonPhase, factor(Region)      Number of observations:  534       I've read through the manuals and vignettes but can't seem to find an answer. Does anyone know how to retrieve the OOB error estimates once you have run a random forest model using the party package? Or am I completely missing some very important difference between the two packages that results in no OOB error estimates for random forest models built with the party package?