title: Building reliable glmnet model and constructing predictions

text: I am new to glmnet. I have over 15,000 predictors and a binary outcome for approximately 400 samples. I am having trouble finding sources online to describe explicitly if what I am planning to do is correct. My current analysis plan includes:   1. Obtain the training parameters. Use the caret package to repeatedly cross validate and find alpha and lambda with maximum ROC. I am only testing alphas between 0.5 and 1 because I want to lean towards LASSO analysis. Code adapted from : ##http://moderntoolmaking.blogspot.com/2011_05_01_archive.html              MyTrainControl=trainControl(       method = "repeatedCV",       number=10,       repeats=5,       returnResamp = "all",       classProbs = TRUE,       summaryFunction=twoClassSummary     )          model <- train(cf.x.df, cf.y.factor,method='glmnet',            metric = "ROC",            tuneGrid = expand.grid(.alpha=seq(0.5,1, by=0.1),.lambda=seq(0,0.3,by=0.01)),            trControl=MyTrainControl)        2. Use training parameters, bootstrap N times to obtain a frequency of non-zero coefficients for each variable. Assume these are the best variables to include in the model.   3. Reduce the data to only these parameters and use previous training parameters train in caret package to bootstrap for prediction accuracy. Although I know this is somewhat controversial, is there a simple approach to obtaining the coefficients for each of the bootstraps in (3) so I can calculate confidence intervals for the coefficients at the same time? I know I can create a single model using a "training" data set and the alpha and lambda obtained as training parameters, and test prediction in a "testing" data set. But it seems so arbitrary to only run the model once. At the same time, I would like to know more about the train function in caret. Can I specify the variables to be included in the model, or should I just create a new data set with only these variables?


title: Interpreting the lm output for football data and corresponding QQ Plot

text: I had the possibility to try some regression on Football data. I receive the following ouput: ![lm output](http://i.stack.imgur.com/zvY5v.png) I receive a quite good R² but I wonder if it´s important that the intercept has no significance and a rather high std. error. I also get the following qqPlot and I assume, my linear model is acceptable. ![QQPlot](http://i.stack.imgur.com/v9ar9.png) Maybe someone could help me to Interpret this result and clarify my doubts about the intercept.


title: Survey Data Analysis - summing questions on different Likert/ordinal scales

text: I am working on a survey analysis looking at the effect of an intervention with three components (each having multiple sub-parts). Questions ask respondents how frequently they used each component. Two of the questions have scales ranging from 1 to 5 (never/rarely/often/almost always/always), while the other question has a scale ranging from 1 to 4 (never/1-2 times/3-4 times/5 or more times). Here are my questions:   1. Is it okay to sum all three questions together despite being on different scales? If not, how should I remedy the scale differences given that they are all measuring the same underlying construct (how the overall intervention was adopted)?   2. Each question has a different number of sub-parts (Q1 has 4, Q2 has 5, and Q3 has 6) b/c the number of items in each part of the intervention varies. Therefore, since the answer scale of Q1 is 1 to 5 and it has 4 sub-parts, the theoretical sum for any respondent could range from 4 (lowest) to 20 (highest), but for Q2, the range is obviously different (5 to 25). Q3 is on a different range (6 to 24) because the answer scale range is 1 to 4 and it has 6 sub-parts. How can I standardize the ranges? Do I even need to? I hope these questions make sense and thanks in advance for any help! It is greatly appreciated!


title: Convergence in distribution of the maximum of a sequence of random variables

text: How to find the following: Let $X_1$, $X_2$, $X_3$,..., $X_n$, be i.i.d with chi-square distribution with one-degree of freedom. Find $a_n$ and $b_n$ such that $ a_n(\max_i X_i - b_n)$ converges in distribution to a nondegenerate random variable. I thought about Central limit theorem but I don't think it is the useful here.


title: Finding overlap in two similar populations

text: I run a few different projects throughout India and want to understand how many unique individuals I've provided services to total. I don't collect names or unique identifiers, so all I know is a) the amount of people that have accessed services per project and b) the total population of India, which can further be refined by states or cities. However, there will be some individuals who access services from more than one project, thus, there will be overlap. What sort of statistical test can I apply to understand the probability of overlap between projects?


title: Bivariate normal distribution and its distribution function as correlation coefficient $\rightarrow \pm 1$

text: I am not sure what happens to a bivariate normal distribution when $|\rho| \rightarrow 1$. Is the distribution well defined then? Moreover, when $$ \Phi \left(\frac{x_1}{\sigma_1}, \frac{x_2}{\sigma_2}, \rho \right) = 1 \text{ or } 0?$$ Namely, is it true that if $\Bigl \lvert \frac{x_1}{\sigma_1} \Bigr \lvert \leq M$ and $\Bigl \lvert \frac{x_2}{\sigma_2} \Bigr \lvert \leq L$, where $M, L < \infty$, then $$ \Phi \left(\frac{x_1}{\sigma_1}, \frac{x_2}{\sigma_2}, \rho \right) \neq 1 \text{ and } 0? $$ Can it be stated without any restrictions on $\rho$? If you know any reference, where it is explained well, I would be very grateful.


title: Imbalanced training dataset and Random Forest regression model

text: I have a large dataset (>300,000 observations) that represent the distance (RMSD) between proteins. I'm building a regression model (Random Forest) that is supposed to predict the distance between any two proteins. My problem is that I'm more interested in close matches (short distances), however my data distribution is highly biased such that the majority of the distances are large. I don't really care how good the model will be able to predict large distances, so I want to make sure that the model will be able to accurately predict the distance of close models. However, when I train the model on the full data the performance of the model isn't good, so I wonder what is the best sampling way I can do such that I can guarantee that the model will predict the close matches distance as much accurately as possible and at the same time now to stratify the data so much since unfortunately this biased data distribution represent the real world data distribution that I'm going to validate and test the model on. The following is my data distribution where the first column represents the distances and the second column represent the number of observations in this distance range:               Distance  Observations     0          330     1          1903     2          12210     3          35486     4          54640     5          62193     6          60728     7          47874     8          33666     9          21640     10         12535     11         6592     12         3159     13         1157     14         349     15         86     16         12


title: Computing p-values

text: I am a computer engineer, living in a digital world. I have no probabilistic background besides a probability course in the University and tried to grasp the Fisher's idea and have got the following questions. ~~Why is this room called Ten Fold and Corss validation?~~ Wikipedia says that "To compute a p-value from the test statistic, one must simply sum (or integrate over) the probabilities of more extreme events occurring". Thereby, we are warned against interpreting p-values as type I errors. But, AFAIK, $a$-values-based Type I errors just mean exactly that you cut the tail distribution whereas correct result resides in it, which is exactly Type 1 error. Secondly, I could not find the formulas or examples computing p-values. Instead WP tells that they were invented by Fisher, who `Rather than using a table of p-values, inverted the CDF, publishing a list of values of the test statistic for given fixed p-values.` Does it say that I must know the p-value in advance? How do I compute it if I must know it in advance? Third, Wikipedia gives a hint on a way to compute the p-value: "Hypothesis tests, such as Student's t-test, typically produce test statistics whose sampling distributions under the null hypothesis are known. To compute a p-value from the test statistic, one must simply sum (or integrate over) the probabilities of more extreme events occurring." This begs a question: do we need convert t- into p-value and apply Fisher test if we can use the t-value and Student test immediately? Why Fisher is better? Finally, Wikipedia says that the formula for p-value is $2/2^n$. But, it is stated in coin flipping context and it is not clear if it is true in every context. Thereby, it points out that rolling two dies, which yield 6+6=12 has probability of 1/36, which is sufficient to overcome the 5% threshold but it is not statistically significant. This leaves me with uncertainty: how can p-value determine the statistical significance after that?


title: R's MANOVA homoscedasticity assumption

text: Does R's `manova` function assume homoscedasticity among all outcomes or only among the subgroups (according to the predictor) of each outcome?


title: Choosing one variable from each of 3 buckets of variables

text: I have a regression model that looks like the following               glm.nb(formula = y ~ Gender + Age + x1 + x2 + x3, data = df)      In my problem, there are 20 possible choices of variables for `x1`, 20 possible choices for `x2`, and 20 possible choices for `x3`. `Gender` and `Age` must be in the model. This leaves me with 20*20*20 = 8,000 possible regressions. I was able to create a program that ran all of these regressions and deliver me the lowest AIC, but I was wondering if there was a library that already does this. I do not consider what I will find to be the "best" model in any statistical manner, but I do find this exercise useful for exploring my data. I have already attempted using bestglm and leaps. I do not believe these programs allow for specifying the choice of variable from multiple bucket of variables.


title: Influence of a time series on another one/clearing up trend or seasonality

text: I've got data as follows               | Date         |    Turnover    |    Day of week    |    Max Temp |     -------------------------------------------------------------------     | 2013-07-08  |    133          |    Monday         |    22       |     | 2013-07-09  |    150          |    Tuesday        |    23       |     | 2013-07-10  |    161          |    Wednesday      |    22       |     | 2013-07-11  |    127          |    Thursday       |    18       |      Now what I need to do is to find out how much influence the Max Temperature has on the turnover. The effects of the time series like trends and day of week (periodicity) should be cleared up before. So that in the end I want to be able to tell that the max temperature influences the turnover either positive or negative when it is between 10-15, 15-20, etc. What I already did with R is to do a time series decomposition. My next plan was to then take the irregular "noise" component from the time series (which as I understand is without any periodicity or trend) and see how the temperature influences the turnover. But here I don't know which statistical method is best suited for this. I need some method to cluster the data maybe. Any hints?


title: How to compare variability within and between groups?

text: I re-edited this question for clarity:   1. I have coordinates obtained by tracking birds with GPS (all points in plot).    2. I've used these points to perform **kernel density estimation** to reveal the areas with most dense bird activity (areas "1" and "2" in plot).   3. I also know coordinates of nesting trees where each nest was protected by female and male (big blue points with labels ("A" ..."J").   4. On each nest were done nest-defense recordings (10 min each recording, on three various intruders, twice per the same intruder)= 60 min total for each nest.   5. Each bird (female and male) exhibit 4 behaviours: "attack", "threat", "jump","check". I have all these values (duration).  ![enter image description here](http://i.stack.imgur.com/ZWRYP.png) Thanks to David Robinson for his time. I've done mistake in refer! I want to refer to **Table 2**. in this article: onlinelibrary.wiley.com/doi/10.1111/jbi.12048/full In article is used analysis of molecular variance (AMOVA). May I use something similar to my data in order to explain: **"source of variation"**   * among areas   * among nesting trees within clusters   * among individuals (female vs. male) within nesting tree Is it a problem that this AMOVA is (probably) "designed" on genetic data? Which procedure in appropriate for my data? Or should I use nested ANOVA? How? Thanks.


title: Some help on genome analysis

text: I have an exam tomorrow about statistics for genome analysis and I'm having some troubles with some questions from an example exam. All help would be very welcome! These are the questions:   1. The single step maxT procedure as explained in the course notes is a nonparametric exact correction and takes the dependence between tests into account. Why do we introduce the false discovery rate if this single step maxT procedure has such nice properties?   2. After preprocessing the probe-level data you obtain an expression value for each probeset on each array. Suppose that you use these probeset intensities to conduct e.g. a t-test for differential gene expression, while ignoring the variability eduction caused by preprocessing. Is this problematic (from a statistical point of view)? Why / why not?


title: Factor models with small noises

text: The standard factor model formulation is $y=W x+\epsilon$ where $x \sim \mathcal{N}(0, I)$, $\epsilon \sim\mathcal{N}(0, \Sigma)$. $W$ and $\Sigma$ are typically estimated from MLE. The solution can be obtained numerically; in general there are no analytical solutions. Question: assume that $\Sigma$ belongs to some class of psd matrices, such as matrices of bounded norm, or with bounded trace. Is there an analytical solution to factor models in which the noise is vanishing, i.e. $y=W x+ n^{-1}\epsilon$ with $n\rightarrow \infty$? By this, I mean than the solutions $(W_n, \Sigma_n)$ converge and that possibly the solution can be expressed in closed form. I know the answer is in the affirmative when it is known a priori that $\Sigma = I$, but I would be happy to see that it holds more generally. Pointers to literature welcome.


title: Estimated marginal means significant, interaction not

text: This is probably pretty basic question, but how come I can get totally non- significant 3-way interaction, while all the marginal means are significant? See the 3 way interaction in the table. ![enter image description here](http://i.stack.imgur.com/XcVUJ.png)


title: Fligner.test with R

text: Why does this:               fligner.test(a ~ b + c + d, data=df)      and this               fligner.test(a ~ b + d, data=df)      always equal for whatever data. In the meantime, this               fligner.test(a ~ b + c + d, data=df)      and this               fligner.test(a ~ c + b + d, data=df)      do not equal. Here is an example of data               set.seed(11)          df = data.frame(a=rnorm(12), b=factor(rep(c(1,2,3),4)), c=factor(rep(c(1,2),6)), d=rnorm(12))      * * * It seems that only the first variable of the formula is considered. Does it mean that if we want to test the homoscedasticity of such model > aov(a ~ c + b + d, data=df) we should run:               filgner.test(a ~ b, data=df)     filgner.test(a ~ c, data=df)     filgner.test(a ~ d, data=df)      and if all of them are not significantly heteroscedastic, we can accept that our model respect the assumption of homoscedasticity. Is it correct? Note that the same effect is found with the bartlett.test(..)


title: given mean and standard deviation, find value of given percentile

text: I'm running a DnD campaign and made a Random NPC generator in Excel. It's pretty slick, but the random height/weight calculator is... flat. RAND(X,Y) flat. So I'd like to add some Bell Curve calculations, because why not. I've entered the average min/max heights and weights for each race, and can easily calculate the mean. Now, if I assume those average min/max values are 1 standard deviation, I've got the ingredients of a formula. (the average height listed for Human was pretty close to real life 1SD range for a human male) It's been... years... since I took a stats course, and I remember almost nothing.   Bell Curve and SD is about it. I spent a few hours poking around, but didn't find any answers, or at least I don't think I did. X__X So, given a Mean, and ranges for 1SD, I'd like to enter a percentile (1d100) and arrive at a height.   For example, a 16% would be -1SD, 50% would be exactly the median, and 85% would be +1SD (following the 68–95–99.7 rule). Hopefully that makes sense. Given the criteria above, could someone please explain what Excel formula to use? Preferably using a built-in one.


title: Measuring homogeneity across different spatial aggregations of data

text: I'm working with dataset of individual households that I aggregate into 'areas' using several different spatial configurations, from smaller to bigger. These areas are then characterized by four variables (two categorical, two continuous). I'd like to see what effects these different aggregations have on the dataset. Particularly, I'd like to estimate what the differences in homogeneity are when I move from one spatial resolution to another. What would be the best way to approach this problem? Is there any measure I could use for this purpose?


title: Does the sum of several distributions become more central or approximated to Normal

text: As the classic CLT states Xs follow the same distribution, then the sum of Xs approximate to Normal distribution. But what about several Xs follows different distributions (maybe the same class but not the same value of parameters or totally different kinds)? It seems that the sum will become more central or even more Normal "likely". Maybe some kinds of "Generalized" CLT has interests about that. Is there any theory doing study about such kind of research? What characteristics if the Xs follow, we will have such results. What kinds of assumptions fails, then we have another results? We know that sum of the uniforms is triangle distribution, while sum of the Poissons is still Poisson. Sum of stable family (including Normal, Cauchy) is still stable distribution. And about student-t, the sum is not student-t but not knowing what kind it is. But it seem that the fat-tail effect of the sum is thinner. So could you please tell me some more general rules about such aspects. Thank you.


title: How to draw estimates based on variance covariance matrix?

text: Suppose I fitted a logistic model and get the estimates as well as their `vcov` matrix. I would realize this: draw length($\beta_s$) independent $\mathcal N(0,1)$ values to create a random vector $z$, then $\beta^*$=$\hat{\beta}+Az$ where $A$ is the upper triangular matrix of the Cholesky decomposition matrix ($\hat{V}=A'A$). How can I draw the $\beta^*$ using the `vcov` matrix? Here is an example:               set.seed(123)     df <- data.frame(y=rbinom(100,1,0.5),                      x1=rnorm(100,10,2),                      x2=rbinom(100,20,0.6))          fit <- glm(y~x1+x2, data=df, family="binomial")      $\hat{\beta}$:               coef(summary(fit))                 Estimate Std. Error  z value Pr(>|z|)     (Intercept)   0.1482    1.57451  0.09413   0.9250     x1           -0.1710    0.10962 -1.55984   0.1188     x2            0.1181    0.09047  1.30567   0.1917      vcov matrix:               vcov(fit)     (Intercept)      2.4791 -0.1225802 -0.1020612     x1              -0.1226  0.0120164  0.0003505     x2              -0.1021  0.0003505  0.0081844      Would somebody know how to draw new $\beta^*$ using the `coef` and `vcov`?


title: Compare regression coefficient for different dependent variables

text: I am currently reviewing a paper in which the authors perform the following analysis. They use the same set of independent variables to explain 3 different variables. The hypotheses are of the form: X has a stronger positive effect on Y1 than on Y2 and Y3. The authors run separate regression models for the three DV's (Y1, Y2, Y3), each time with the same IV's. The sample is also identical, they just use a different scale to define the DV for each regression model. so: Y1 ~ a1 + b1X1 + c1X2 + d1X3 + e1 Y2 ~ a2 + b2X1 + c2X2 + d2X3 + e2 Y3 ~ a3 + b3X1 + c3X2 + d3X3 + e3 with hypotheses such as: c1 > c2 and c3, and all > 0 The authors then simple compare c1 with c2 and c3 and conclude that their hypothesis is supported because the fitted c1 is indeed positive and larger than c2 and c3. To me, this feels incomplete. Seeing that c1 is larger than c2 doesn't mean the difference is statistically significant. How should I respond in my review? I'd like to comment on this if the analysis is incorrect, but would also like to be able to suggest how the analysis could be done correctly. I hope you can help.


title: Help finding method of estimating regression parameters for a given performance

text: Pardon my ignorance, for I have very basic training in statistics, but have never taken a regression course. I ran an experiment where participants are in one of two conditions (between subjects). They categorized an object that had values for rotation and spacing, both ranging from 1 to 100. Each participant saw 40 of those values twice. There were a lot of participants. What I want to know is for what values of rotation (or spacing) were the participants as a whole closest to chance (.5, is it or is it not in the category) in each condition, and whether those values differ significantly from one another. What is this process called, and what are some recommended packages or methods in SAS, R, SPSS, Mathematica, or Stata (or really anything) I could use to go about figuring it out. Any help rephrasing the question is also very welcome.


title: Whether lead content readings follow a Normal distribution or not

text: > Studies indicated that drinking water supplied by some old lead-lined city > piping systems contain harmful levels of lead. It appears that the > distribution of lead content readings for individual water specimens has > mean .033 mg/L and SD .10 mg/L. > > Explain why is it obvious that content readings are **not** normally > distributed. I cannot see the answer instantly. Why can't we perceive this distribution as normal with very high variance?


title: How to collect related words with specified one?

text: Assume I have a word such as `statistics`, and research articles or text books about statistics. My question is how to detect top 100 related, highly associated words, with `Statistics`. My intention to extract words from PDF documents and then determine top correlated words with `statistics`. How can I determine which words are correlated with `statistics`?


title: How to prune a decision tree properly in R

text: I have a sample of 12,500 observations and 12 explanatory variables. I want to build a pruning decision tree, to do that I am using the `rpart` function and then the `prune` function. My problem is that I am not sure how to build the full tree, I am using               rpart(*, *, *, control = (minbucket=x, minsplit=y, cp=z))     # the * are irrelevant to my question      to build the full tree, but I don't know what `x`, `y` and `z` should I use. All the examples I have found are with much fewer observations compared to my case.


title: How do I test differences in reaction times?

text: I have data from an experiment in which participants had to respond to stimuli in two different conditions in a within-subject design. Obviously, I can fairly easily test for differences in accuracy between the two conditions. However, my initial idea was to check for differences in reaction times. When a subject incorrectly answers, I feel the associated reaction time is meaningless. What is the 'standard' way to do this?


title: Does normal distribution assumption improve the quality of results of a meta-analysis of effect-sizes?

text: Different methods of meta-analysis have different assumptions and one of the popular assumption that has been invoked so often is that of normal distribution. can we work out meta-analytic results without this assumption especially when a meta-analysis utilizes a number of studies scattered over a range of populations.


title: Evaluation of IR approach without test collections

text: I am trying to evaluate an information retrieval approach. More specifically, it's a query expansion algorithm, based on topic distribution in retrieved documents. I want to evaluate my approach using the world wide web. Since most (all?) evaluation approaches use a test collection, where the whole corpus is known, I am having a hard time to find an evaluation approach which works only on retrieved documents. So, I am looking for information on this topic. Maybe a paper or a keyword?


title: Can you put more than one factor in ANCOVA analysis?

text: I am running ANCOVA on SPSS using post-operation score (continuous) as Dependent variable, pre-operation score (continuous) as Covariate, and 2 groups (categorial 1&0) as Fixed factors. This is the usual scenario. But can I add another 2 variables one is categorial (gender) and one continuous (age) to the Fixed Factors all at same time? My idea is to check which one of these factors (groups, gender, age) has significant effect and higher partial η2.


title: Statistical Tests for comparing groups

text: I am a biologist and I have developed an algorithm that has the ability to identify and quantify fluorescent objects within the nucleus of a cell when imaged in a fixed state under a microscope. The program will first identity nuclei, then analyze each nuclei individually and determine a threshold which will be used to decide what to consider an object (e.i. 10 nuclei, 10 unique thresholds). There are other programs that can do something similar but it makes a single threshold for an entire image and one threshold is used for all nuclei (e.i. 10 nuclei, one threshold for all nuclei). Ultimately I will compare the manual count (control) against my automated count. But first, I will compare the control against 3 thresholds which will be set at the value 2.5%, median and 97.5% of the measured median intensities of all nuclei within a sample (n = 100). Since I'm counting the number of objects within an area, and from previous knowledge I know the data follows a Poisson distribution and I need to uses a non parametric Friedman test followed by a post hoc Dunnets test to compare every group to the control. I also plan to correlate each each group against the control use an orthogonal comparison. Is there any other type of statistical test that would be interesting to use when comparing the manual count vs thresholds set at 2.5,median and 95.5%? What about when comparing my automated count against the other groups?


title: Permutation tests and multiple testing

text: I'm using permutation tests, based on **random forest** importance measures, to perform **variables selection** in my dataset. Using $99$ permutations, my permutation tests $p$-values are comprised between $0.01$ and $1$. As I am using this procedure with more than $300$ covariates, I considered the issue of multiple testing. However, with all the $p$-value adjustment procedures I considered (BH, q-value, etc.) none of my $p$-values remained significant. Is there an easy way to deal with **multiple testing** when using **permutation tests**?


title: Importance of multivariate normality assumption for BIC-like sparse model selection inference with PCA

text: I am reading a paper for robust, sparse PCA in which they propose a BIC-like criterion for selecting the appropriate value of the sparsity parameter $\lambda$. They define this as: $$BIC(\lambda)=\frac{\tilde{RV}}{RV}+df(\lambda)\frac{log(n)}{n}$$ where $(\tilde{RV}, RV)$ refer to the total robust variance of the residuals matrix obtained from a sparse PCA and an unconstrained PCA. An assumption of of PCA is that the data are jointly normally distributed. I often read that this assumption is largely ignorable when performing data exploration, and the comments usually make the point that PCA is not typically used for inference. But here it is. Focusing on the RV term in the equation, this seems to be the proxy for the likelihood function in a normal BIC. I understand the second term, which penalizes for model complexity. My question is that if my data is not jointly normally distributed, will this potentially invalidate inferences made from this information criterion (which is in model selection serves the role of a p-value: (Equivalence of AIC and p-values in model selection), and will these problems arise because the RV term is inappropriate? I notice that I get selected values of $\lambda$ which often seem too high to be intuitively reasonable, indicating that perhaps the RV term is smaller than it should be.


title: Does correlation imply a (relatively) constant difference or ratio between the two variables?

text: Most of you will have been entertained by the recent blog post showing how random things in the world correlate without, of course, necessarily being causally related. Looking at those graphs, I was wondering in my mind whether, for two vectors X and Y to be correlated, it is the _difference_ or the _ratio_ between their respective elements that needs to be relatively constant between elements. In the graphs, the two variables are always on very different scales, and so it is the difference that is constant rather than the ration. Did a quick check in Statistica and turns out that either is sufficient: ![enter image description here](http://i.stack.imgur.com/SUzm8.png) ![enter image description here](http://i.stack.imgur.com/IsK06.png) Or am I making a mistake here? Looked at the definition of Cov(X,Y) but didn't help me understand better.


title: Standardizing count variables in panel data with overdispersion - R or Stata

text: I'm running a regression where the dependent (response) variable is a highly dispersed (slightly zero-inflated) count and the explanatory (independent or predictor) variables are continuous, counts as well as binary. I read a couple of interesting contributions on this topic on SO (here, here, and here as well as the articles / blogs they refer to), yet they generally deal either explicitly or implicitly with response variables that are continuous. So I have a couple of questions. Any help and (partial) answers would be much appreciated.   1. Given overdispersion of the response variable (variance >> mean) and a mean that varies over time, what is the best way to 'change' the response variable?   2. @Gung argued in one of his/her contributions that centering or scaling predictors is especially useful when you plan to interact them with other variables or in case of quadratic terms. This would reduce collinearity. Now in my data, I find that the correlation between the mean-centred variables and their quadratic terms is indeed smaller (from .88 to .80), however, the correlation between the different mean-centred variables x1_centered and x2_centered is actually bigger. Hence I'm not sure about which effect is worse...   3. It is often suggested (e.g. Angrist & Pischke, 2009 - Mostly harmless econometrics) to mean-center the response variables in panel data as this would account for fixed effects. But this raises 2 questions:      1. Which mean do you use in panel data? Using the mean of the entire sample would imply knowing the future before it unfolds. Using a rolling mean that includes the lagged variables for each year of regression might work but I have never seen it done in a paper.     2. What is the interpretation of a mean-centered count variable. In my case I look at patent citations and mean-centered would mean I have negative values suddenly   4. Finally, by mean-centering the response count variable (i.e. if there is a good reason to do so, I have to move away from a negative binomial regression because the nature of the response variable changes from positive integers to positive and negative non-integers. As such, I think the points made by @MansT regarding least squares estimators (see here) are not valid (I don't imply @MansT argues they are!!). Clearly, IF I have to change the regression (from glm, familiy = quasibinomial()) to something else, this will affect all the betas. Ok, I know it's a long question so I can only hope some of you can contribute to parts of it. Thanks in advance! Simon


title: significance Compairing two independent groups with single categorical variable

text: I want to see the significance between two group (One with Albandazole and without Albandazole) inrespect to their CT scane Resolution with three catories(Compltele,partially resolved, Calcified)


title: ROC curve crossing the diagonal

text: I am running a binary classifier at the moment. When I plot the ROC curve I get a good lift at the beginning then it changes direction and crosses the diagonal then of course back up, making the curve a tilted S like shape. What can be an interpretation/explanation to this effect? Thanks


title: Bayes rule and conditional independence

text: I have two conditionally independent random variables $A$, $B$ such that $$ P(A,B\mid C) = P(A\mid C)P(B\mid C) . $$ I have to find posterior formula $P(C \mid A,B)$. My result with a straigthforward application of Bayes rule is $$ P(C \mid A,B) = \frac{P(B\mid C)P(A\mid C)P(A)}{P(A\cap B)} . $$ with few variants (e.g. get an intersection on numerator). But I can't get the lecturer's solution that is $$ \frac{P(B\mid C)P(C\mid A)}{P(B\mid A)} . $$


title: Transforming data with both negative and positive skew

text: I am attempting to run a MANCOVA of memory assessment data. My IV has three levels (no memory impairment, moderate memory impairment, and severe memory impairment). The covariate is education. The DV is 7 memory assessment scores, with the scores being recorded as either percentiles or scaled scores. The problem I am encountering is meeting the assumption of univariate and multivariate normality with my data. I have several outliers (mostly because I am comparing patients with varying degrees of memory impairments on measures assessing memory). The IV groups are rather small (22, 26, and 27). I have been using the DeCarlo SPSS Macro to assess multivariate normality. When I run the macro without transformed data, all of the subjects in my no memory impairment group have negatively skewed data ranging from (-.5 to -1.7) and all of the subjects in my 2 memory impaired groups have mostly positively skewed data ranging from (.6 to 1.03) for the moderately impaired and (.7-2.7) for the severely impaired. Further the Omnibus test of multivariate normality is significant for all IV levels, and the omnibus test of univariate normality is significant for all but one of the DVs in the severely impaired group. Any suggestions of transformations I can use to ensure my data is meeting the normality assumption?


title: Mediation analysis interpretation

text: I am carrying out a mediation analysis using PROCESS in spss. According to Field a mediation effect occurs if the confidence intervals for the indirect effect do not include 0 - this is the case but the predictor to mediator pathway is non significant, although there is a sign effect for predictor to outcome, and sign effect mediator to outcome. Does this still mean there is a mediator effect?


title: What is the difference between “factors” and “covariate” in terms of ANOVA?

text: I am a bit confused on the term "covariate". It seems like the term can mean two different things. In ANCOVA, the term is used for the third variable that is not directly related to the experiment. For ex: age or IQ on performance study (comparing) between Male and female in Standardized test i.e. IQ is used as a covariate. In ANOVA/regression design, "covariate" just refers to factors/Ind. var? I may have completely misunderstood this. Can anyone give a simple example of the term "covariate" used in different context?


title: How to understand SARIMAX intuitively?

text: I'm trying to understand a paper about electric load forecasting but I'm struggling with the concepts inside, specially the **SARIMAX** model. This model is used to the predict the load and uses many statistical concepts that I do not understand (I'm an undergrad computer science student -- you can consider me a layperson in statistics). It is not necessary for me to understand completely how it works, but I'd like to at least understand intuitively what is happening. I've been trying to split SARIMAX into smaller pieces and trying to understand each of these pieces separately and then putting them together. Can you guys help me? Here's what I have so far. I started with AR and MA. **AR** : _Autoregressive_. I have learned what a regression is, and from my understanding, it simply answers the question: given a set of values / points, how can I find a model that explains these values? So we have, for example, the linear regression, which tries to find a line that can explain all these points. An autoregression is a regression that tries to explain the values using their previous values. **MA** : _Moving Average_. I'm actually quite lost here. I know what a moving average is, but the moving average model doesn't seem to have anything to do with the "normal" moving average. The formula for the model seems awkwardly similar to AR and I can't seem to understand any of the concepts I find in the internet. What is the purpose of MA? What's the difference between MA and AR? So now we have ARMA. The **I** then comes from _Integrated_ , which as far as I have understood, simply serves the purpose of allowing the ARMA model to have a tendency, either increasing or decreasing. (Is this equivalent to saying ARIMA allows it to be non-stationary?) Now comes the **S** from _seasonal_ , which adds periodicity to ARIMA, which basically says, for example in the case of load forecasting, that the load looks very similar everyday at 6 PM. Finally the **X** , from _exogenous_ variables, which basically allows external variables to be considered in the model, such as weather forecasts. So we finally have SARIMAX! Are my explanations OK? Recognize that these explanations aren't required to be rigorously correct. Can someone explain me what MA does intuitively?


title: Is there a statistical technique to perform this comparison? How could I do this on Stata?

text: I've been busy on a work where I have to compare two financial econometrics models on the determinants of financial leverage (panel data). These have only few control variables and the dependent variable in common, and some other independent variables differ. So, these are:   1) $\text{Leverage} = a + bX_1 + bX_2 + bX_3 + bZ_1 + bZ_2 $ 2) $\text{Leverage} = a + bX_1 + bX_2 + bX_3 + bZ_3 + bZ_4 + bZ_5 $ Can you suggest me any technique to compare these models, if possible? How could I do this on Stata?


title: Forecasting price based on historic time series

text: I work for a material stock holder. I am trying to quantify how the probability of not getting a specified minimum of gross profit, changes over time. This problem basically boils down to assigning probabilities to the set of future market prices. I was hoping that I could use time series analysis on the company's purchase and sales records (they span the last ten years, so I should hopefully have a decent dataset), but I am unsure about of how to go solve this. I have done some googling and from what I can see the an ARIMA modeling is what I should be looking into, but again I am unsure. If anyone could give me some pointers that would be great. The price of the material we stock fluctuate over time, (of course, as a result of the changes in the factors that influence the market price.) Plotting the prices into a histogram, it is clear that it is not a normal distribution. (this is the case for both deflated and nominal prices) Which, from my standpoint is expected, prices are influenced by business cycles and the evolution of the production costs. If anyone could give me a pointer in the right direction of where to start looking / reading, I would be very grateful.


title: Confidence size and coverage probability in a confidence set?

text: Let $\theta \in \Theta \subseteq \mathbb{R^d}$ be the parameter of interest and let $\theta_0$ be the true population parameter value. Let $n$ be the sample size. Let $CS_n$ be the confidence set for $\theta_0$ of level $(1-\alpha)$ obtaining by inverting a test, i.e. $$ CS_n=\\{\theta \in \Theta | T_n(\theta)\leq c_{1-\alpha}(\theta) \\} $$ where $T_n(\theta)$ is an appropriate test statistic and $c_{1-\alpha}(\theta)$ is the associated critical value. Some definitions: (1) Exact coverage probability of $CS_n$: $\mathbb{P}(\theta_0 \in CS_n)$. It should be $\geq (1-\alpha)$ (2) Asymptotic coverage probability of $CS_n$: $\liminf_{n\rightarrow \infty}\mathbb{P}(\theta_0 \in CS_n)$. It should be $\geq (1-\alpha)$ (3) Exact confidence size of $CS_n$: $inf_{\theta \in \Theta} \mathbb{P}(T_n(\theta)\leq c_{1-\alpha}(\theta))$ (4) Asymptotic confidence size of $CS_n$: $\liminf_{n\rightarrow \infty}(inf_{\theta \in \Theta} \mathbb{P}(T_n(\theta)\leq c_{1-\alpha}(\theta)))$ Questions: (a) Why in (2) we need $\liminf$ instead of $\lim $? (b) Relation between (1) and (3)? (c) Relation between (2) and (4)? (d) Why in (3) and (4) we need to take $\inf_{\theta \in \Theta}$?


title: Finding the integral of a fitted function

text: I have a function obtained by fitting some data, and I do not have access to the data itself. The fitting parameters of the function have confidence bounds. I need to obtain an expression for the integral of this function with corresponding errors. I have thought of solving the problem with a Monte Carlo approach. By integrating perturbations of the fitting function by varying the fitting parameters (assuming some probability distribution, e.g. Gaussian), I could then fit this mock data by a new function, and obtain new fitting parameters with new confidence intervals. Is there a more general procedure, or is this a decent approach? This approach assumes no covariance among the fitting parameters which may or may not be a bad assumption. Any good sources on how to approach this problem are greatly appreciated. **EXAMPLE** Suppose we have the fitted function $f(x)=x^{\alpha} \exp\left(- x^{\beta}\right),\alpha=1\pm0.02,\beta=3\pm0.3$. We then wish to obtain an expression for $F(t)=\int^{\infty}_tf(x) dx, t>0$. How do I find the best fit to this function, and how do the errors of the $f$ fit propagate to $F$?


title: Discrete variables in regression model?

text: I know that in theory for regression both the Y and factors should be continuous variables. However, I have some factors that are discrete but show both correlation and would fit a regression model. I am looking at energy consumption and my factors are the number of calls, the data transferred, temperature, customers, number of buildings. The only continuous variable is the temperature. Any advice if my multiple regression model is still applicable even though I have factors temperature and subscribers? How about correlation? Would I be able to say that there is in fact say correlation between the number of customers and energy consumption? Thanks for your help!


title: Smoothing with fda package when time grids differ

text: In fda package, function smooth.basis is used to construct a functional data object by smoothing data using a roughness penalty. The usage can be found here.               smooth.basis(argvals=1:n, y, fdParobj, wtvec=NULL, fdnames=NULL,                  covariates=NULL, method="chol", dfscale=1, returnMatrix=FALSE)      In the documentation, it said > _**argvals:_** Argument values can even vary from one variable to another, > in which case they are input as an array with dimensions corresponding to > observation points, curves and variables, respectively. Note, however, that > the number of observation points per curve and per variable may NOT vary. If > it does, then curves and variables must be smoothed individually rather than > by a single call to this function. **My question is how to set it up to use the function when we have subjects with different number of observations observed at different time series.** For example, subject 1 has 10 observations at time week 1 to week 10, but subject 2 has only 5 observations at week 1.5, week 2.5, week 3.5, week 4.5, week 5.5. Do I need to use the same basis system like b-spline to smooth two curves to make them comparable or consistent? Can anyone give me some examples to illustrate it? Thanks.


title: A special case of GMM estimation in R

text: I want to estimate the forward looking version of the Taylor rule equation using the iterative nonlinear GMM: ![Model](http://latex.codecogs.com/gif.latex?i_t%20%3D%20%5Calpha%20+%20%281-%5Crho%29%5Cbeta%5Cpi_%7Bt+1%7D%20+%20%281-%5Crho%29%5Cgamma%20x_t%20+%20%5Crho%20i_%7Bt-1%7D%20+%20%5Cvarepsilon_t) I have the data for all the variables in the model, namely ![pi](http://latex.codecogs.com/gif.latex?%5Cpi) (inflation rate), ![x](http://latex.codecogs.com/gif.latex?x) (unemployment gap) and ![i](http://latex.codecogs.com/gif.latex?i) (effective federal funds rate) and what I am trying to estimate is the set of parameters ![i](http://latex.codecogs.com/gif.latex?%5Cbeta), ![i](http://latex.codecogs.com/gif.latex?%5Cgamma) and ![i](http://latex.codecogs.com/gif.latex?%5Crho). Where I need help is in the usage of the `gmm()` function in the `{gmm}` R package. I 'think' that the parameters of the function that I need are the parameters:               gmm(g, x, type = "iterative",...)      where `g` is the formula (so, the model stated above), `x` is the data vector (or matrix) and `type` is the type of GMM to use. My problem is with the data matrix parameter. I do not know the way in which to construct it (not that I don't know of matrices in R and all the examples I have seen on the internet are not similar to what I am attempting to do here. Also, this is my first time using the `gmm()` function in R. Is there anything else I need to know?


title: How to perform a basic forecasting model from pooled cross-sectional timeseries data in SPSS?

text: For 86 companies and for 103 days, I have collected (i) tweets (variable `hbVol0`) about each company and (ii) pageviews for the corporate wikipedia page (`wikiVol0`). The dependent variable is each company's stock trading volume (`stockVol0`). My data is structured as follows:               company  date  hbVol0   wikiVol0  stockVol0  comp1  comp2 ... comp86  marketRet     -------------------------------------------------------------------------------     1        1     200        150     2423325      1      0   ...   0     -2.50     1        2     194        152     2455343      1      0   ...   0     -1.45     .        .      .          .         .         .      .   ...   .     1       103    205        103     2563463      1      0   ...   0      1.90     2        1     752        932     7434124      0      1   ...   0     -2.50     2        2     932        823     7464354      0      1   ...   0     -1.45     .        .      .          .         .         .      .   ...   .     .        .      .          .         .         .      .   ...   .     86      103     3          55      32324       0      0   ...   1      1.90      As you can see, the dataset is a pooled cross-sectional time series data. Although other models would be better suited, my thesis coach was OK with me running an OLSon it. After taking the log for `stockVol0` and `wikiVol0` and taking lag 1 and 2 for all variables, this is one of the best-fitting regressions: ![enter image description here](http://i.stack.imgur.com/4KSkQ.png) Now, my coach tells me "to do some predictive modelling". Since we have had very little statistics classes, I am unfamiliar with this. My coach told me "to do an easy random forest model", but SPSS (which is the only software I have and am somewhat used to) doesn't have a "random forest" model in it's options. It does support ARIMA, an I have been told it can do forecasting as well. So I tried, but SPSS interprets all (103 x 86) 8,858 rows in the dataset as being unique dates. But there are only 103 days (see `date` variable in dataset). **Question 1:** can anybody tell me how to tell SPSS I have 86 observations for each of the 103 days and not 8,858 dates? **Question 2:** am I right that ARIMA is a forecasting model? If not, what is a basic forecasting model that I could apply and/or is random forest available in SPSS under a different name? P.S. I have asked around among my friends, and 90% of my fellow student graduate without having any predictive modeling in their thesis. It's just that my coach really wants it in there, so a basic model will suffice. I hope some people here can help me on my way since my thesis is due in 2 weeks and I am totally lost here...


title: Inverse Incidence Rate Ratio

text: I wondered if it is possible to calculate the inverse of an IRR like we can for odds ratio (inverse = 1/x). For example: IRR = 1.40, 95% CI = 1.18 -1.65 Inverse IRR: = 1/1.40 = 0.71 Thank you (Using Stata 11)


title: Basic Probability Question

text: Given 3 phillips head screws, and assuming there are 5 different orientations that are possible when the screw is finally placed into the wood, what is the probability of all 3 screws having the exact same orientation/alignment ? In other words the slots of the screws are all lined up the same.


title: Multicollinearity: is a matrix scatterplot enough to rule it out?

text: I am analyzing an environmental data set containing one response variable (R) and many explanatory variables. The explanatory variables are either factors (F1, F2, F3) or continuous variables (V1, V2, V3…). I wanted to conduct a linear regression in the form lm (R~F1+F2+F3+V1+V2+V3+F1*F2+F2*F3…) The matrix scatterplot did not show any collinearity issue between my explanatory variable, but from previous analysis (ANOVA) I know that the factors are significantly explaining variation in the continuous variables. (i.e. aov (V1~(F1+F2+F3)^3), all factors significant; (V2~(F1+F2+F3)^3), all factors significant….) The vif in the linear regression model for each explanatory variable is <3. So my questions are: 1) If an ANOVA analysis points out a relationship between subsets of explanatory variables are these to be considered collinear? 2) Should then I remove all the continuous variables from the analysis?


title: 2SLS with two instruments for one endogenous variable in MATLAB

text: I have one endogenous variable and two instruments for it, and I want to calculate my beta with the direct (one step) matrix formula $\beta_2sls = X' Z(Z'X)^{-1}Z'X^{-1}X'Z(Z'Z)^{-1}Z'Y$ But if I have two instruments for one endogenous variable X and Z are not the same length. Any ideas? Thanks!


title: What type of statistic should I use to compare ratios in my experiment?

text: I have performed an experiment to analyse the ratio of intronic versus exonic levels for 1000 transcripts after treated the cells with a drug and compare my results with untreated cells . Each transcript is equivalent to a gene, while intronic regions represent a part of the transcript that is rapidly remove and exons the part of the transcript that encode for proteins. The drug is thought to potentially increase the levels of intronic regions compare to exonic regions. I have compared the ratio of intron/exon signals in treated versus untreated cells in two independent experiments (two biological replicates). What kind of statistic test should I use to see if there is any significant increase in intronic/exonic regions after drug treatment? I was performing a student´s t-test comparing the ratio of the intron/exon at biological replicates 1 and 2 treated with the drug with untreated cells, is that correct? Thanks a lot for your help


title: does rstandard standardize in z?

text: I'm new to R, so please be gentle. I was under the impression that `rstandard(model)` returns the z-scores of the residuals in `model`. However, when I standardized the residuals myself in z, the result was different. In fact, `rstandard(model)` had a mean different from 0 and had a standard deviation different from 1. The differences seem to be nonnegligible. What exactly does `rstandard(model)` do? Or am I doing something wrong here?               > x = rnorm(30, mean = 100, sd = 15)     > y = rnorm(30, mean = 100, sd = 15)     > model = lm(y ~ x)     > z.res = residuals(model)/sd(residuals(model)) #standardizing it myself     > rstandard(model) - z.res #difference between rstandard and what i did                 1             2             3             4             5      -4.422354e-04  1.556269e-04 -4.576832e-03 -1.274350e-03  1.048068e-01                  6             7             8             9            10      -2.333922e-02  1.820134e-02 -3.307542e-03  3.368978e-02 -1.804108e-04                 11            12            13            14            15      -1.100621e-01 -1.343715e-03 -1.300427e-03  1.509862e-03  3.246602e-03                 16            17            18            19            20       3.734255e-03 -1.821539e-06 -1.153190e-02 -1.713254e-06 -2.185101e-02                 21            22            23            24            25      -2.681935e-02  2.562472e-03 -4.721114e-02 -1.084481e-04 -3.430827e-03                 26            27            28            29            30       4.149684e-04  7.705807e-04  2.166815e-03  2.537837e-02  4.182761e-04      > mean(z.res)      [1] -9.428041e-18      #as expected of z-scores, mean is 0     > sd(z.res)     [1] 1      #as expected of z-scores standard deviation is 1     > mean(rstandard(model))      [1] -0.001990908     #not really 0     > sd(rstandard(model))      [1] 1.019699     #not really 1      Also, the way I understood Standardized residuals in R's lm output, `rstandard` is actually studentized residuals. But isn't there already `rstudent`? I'm using R version 2.14.1 in Xubuntu 12.04. Thank you.


title: Visualizing large file-based or Redis-in-memory stored large datasets (millions of data points)

text: I am very active at StackExchange's QuantFinance forum but thought this question is more suitable to be asked here. I am generating large time series data and store them in-memory in Redis (alternatively could also save to disk in any format) and like to visualize the data efficiently as part of my data mining task. The data is of the following structure (but I could adjust the structure to fit the format of a capable visualization package): column 1 key/index : millisecond precision time stamp (.Net long type) column 2 - n : data containing time series which could be of type string, bool, double, int,... I look for a visualization package with the following features:   * Handle time series data that contain millions of data points   * Able to either load the data set from file or access memory-mapped data (preferably Redis)   * Freely select which variable to chart on x- or y-axes.   * Chart several series on the same chart   * Jump to specific windows within the complete data set to only chart subsets (by start and end time (column1)   * Different chart types (line, financial,...)   * Quickly scroll through the complete data set by either loading the complete data set into memory or loading on demand (visible data only).   * Flexible zooming in/out to highlight specific windows   * Should be GUI based I have looked at a couple solutions but am so far not satisfied:   * R/RStudio : Capable analytical packages but very poor charting capabilities. Need to specify in scrips each time I want to look at a different subset of data, poor zooming capabilities. But great to get the data in through memory map packages.   * RapidMiner : Looks promising in terms of visualization capabilities but so far not sold on data imports (I am not aware it works with reading data from Redis in-memory server, or memory mapped files). Also, I am a bit skeptical in regards to performance. I need to process data sets that potentially contain tens of millions of data-points.    * C# Charting libraries: I generate the data within C# and store in Redis In-Memory. So far only Sci-Charts (commercial charting library) looks like it could handle the size of data. Downside, I would need to write a lot of code to achieve the above mentioned charting features.   * Teafiles -> TeaHouse: Proprietary data format for time series. Seems to have excellent charting capabilities but I am a bit worried that it is slightly too close-source (in terms of not being able to extend charting capabilities). Also, I would need to save the data in the TeaFile specific format and save to disk, no Redis support. While this is not a requirement, I am a bit put off because this is a pure charting application which does not offer the slightest time series analysis tools.  I am running on Windows (store data from within a C# application) and therefore am not able to consider Linux or Unix based solutions. Any pointers or ideas what I should look or focus on? Again I like R a lot, but just have not come across GUI based charting packages that could handle what I am looking to do. But I am sure I have not come across all packages that offer charting capabilities. Thanks a lot for your recommendations.


title: How to test whether correlation measures differ when controlling or not for a third variable?

text: If the correlation between demographic dissimilarity and satisfaction is $r=.-14$ and the partial correlation, with career development partialled out, between demographic dissimilarity and satisfaction is $r=-.06$ across a very large sample of size $n$, what is the appropriate test to determine if these correlations are significantly different? Steiger (1980) appears to be the authoritative article on this but the article and most other sources assume one is comparing the correlation between x/y and v/y with dependent groups. This one has stumped our school's statistic resource center.


title: Multilevel model vs. separate models for each level

text: What are the advantages and disadvantages of running separate models vs. multilevel modeling? More particularly, suppose a study examined patients nested within doctors' practices nested within countries. What are the advantages/disadvantages of running separate models for each country vs. a three level nested model?


title: R's coxph won't converge when I include factor (categorical) variables

text: I have a dataset of 371 observations. When I run `coxph` with numeric variables it works fine. However, when I try to add factor (categorical) variables it returns “Ran out of iterations and the model did not converge”. Of note, when I restructure all factors to binary variables with dummy and use glmnet-lasso the model converges. Here are examples of the code and output (including summary description of the variables):               > maxSTree.cox <- coxph (Surv(time,status)~Chemo_Simple, data=dataset)          Warning message:     In fitter(X, Y, strats, offset, init, control, weights = weights,  :       Ran out of iterations and did not converge          > summary (dataset$Chemo_Simple)              Anthra-HDAC       Anthra-Plus       ArsenicAtra              ATRA           ATRA-GO                      0               163                 2                12                 0                 2               ATRA_IDA Demeth-HistoneDAC          Flu-HDAC     Flu-HDAC-plus         HDAC-Clof         HDAC-only                      0                34                37                 4                24                 1              HDAC-Plus        LowArac+/-       LowDAC-Clof         MYLO+IL11    No Rx in MDACC            Phase1                      4                 8                30                 5                 1                 5                    SCT    StdARAC-Anthra      StdAraC-Plus          Targeted         VNP40101M                      0                 0                 0                13                23


title: How to add standard error to plots in ggplot2 with R?

text: I'm in the process of learning R, in the hope of replacing everything I do in SPSS/Sigmplot with R. It's going well so far :) I've got to the point of running a repeated-measures ANOVA, but have come unstuck when trying to plot the results I've worked out how to plot a set of means using ggplot2, but now I'm unsure of how to plot the standard error as error bars. I've seen a number of guides with different implementations, and none of them seem to be appropriate (or even agree with each other). Many people use standard deviations, which is not what I am after. Others have different methods of computing the standard error, so I'm unsure of the best way to proceed. What I have so far is this:               qplot(CATEGORIES, means, shape=factor(ANOTHER_CATEGORY), facets=MORE_CATEGORIES ~ ., data=alldata)      I was wondering if someone could point me in the right direction in terms of how to get the standard errors from a repeated-measures ANOVA in R, and then how to translate this into error bars in ggplot? Thanks!


title: How to set up and estimate a multinomial logit model in R?

text: I ran a multinomial logit model in JMP and got back results which included the AIC as well chi-squared p-values for each parameter estimate. The model has one categorical outcome and 7 categorical explanatory vars. I then fit what I thought would build the same model in R, using the `multinom` function in the nnet package. The code was basically:               fit1 <- multinom(y ~ x1+x2+...xn,data=mydata);     summary(fit1);      However, the two give different results. With JMP the AIC is 2923.21, and with `nnet::multinom` the AIC is 3116.588. So my **first question** is: Is one of the models wrong? The second thing is, JMP gives chi-squared p-values for each parameter estimate, which I need. Running summary on the multinom `fit1` does not - it just gives the estimates, AIC and Deviance. My **second question** is thus: Is there a way to get the p-values for the model and estimates when using `nnet::multinom`? I know mlogit is another R package for this and it looks like its output includes the p-values; however, I have not been able to run `mlogit` using my data. I think I had the data formatted right, but it said I had an invalid formula. I used the same formula that I used for `multinom`, but it seems like it requires a different format using a pipe and I don't understand how that works. Thanks.


title: estimating period and dealing with Non negative values in forecasting

text: When I read time series in a ts object and put a period: 1) tr <\- ts(data[,4],frequency=). This works for two different periods and decomposes perfectly to show (downward) trend, seasonality and error. How do I know which is the correct period. 2) When I use a ETS or STLF function in the forecast package on the above ts object, the summary shows: Model Information: ETS(A,N,N) Why is that? We do have a seasonal+trend component here? 3) What can we do to deal with negative values when the training set is non-negative and negative values make no sense. Thanks, Ratheen


title: What statistical test to use to examine effect of condition on states over time?

text: I need to know the types of tests I need to be thinking about for this stats project. I will then go and look up the specific tests you recommend. The data is from all 50 states for the years 1960-2010. I think the best way to illustrate the data is to show an example (I just made up numbers):               State           Condition       Score     Alabama1995     Condition 1     .534     Alabama1996     Condition 1     .343     Alabama1997     Condition 1     .545     Alabama1998     Condition 2      NA     Alabama1999     Condition 2     .344     ....     Alaska1995      Condition 2     .434     Alaska1996      Condition 2     .444     Alaska1997      Condition 2      NA     Alaska1998      Condition 3     .456     Alaska1999      Condition 3     .678     Alaska2000      Condition 3     .437     ....     Maine1995       Condition 3     .343     Maine1996       Condition 3     .342     Maine1997       Condition 3     .543     Maine1998       Condition 1     .439      So what I want to know is if the conditions and scores are statistically different. Do states when under condition 3 have a significant difference from states when under 2 or 1? Let's say a state in a given year was either condition 1, 2, 3. Not all states could have a score computed for them, but most did. I want to evaluate the following types of things:   * What is the relationship between conditions and scores?   * What is the relationship between time and conditions and scores?   * When a state changed from condition 1 to condition 2, was the change in score significant?   * Does condition 1 overall have a higher/lower score than other conditions? Are there other things besides simple significance tests I should be thinking about? Ideally something that R (Rcmdr!) can do nicely. Stata, SPSS are also available to me.


title: Public databases of learned HMM models for NLP

text: I understand that HMM models model language with Parts of Speech (POS) as hidden states and words as observations. These HMM models are usually learned from large text corpora, and many of these corpora are publicly available. Where can I find such models with their parameters, i.e., the list of POS or other hidden states, observed words, the transition probabilities, and emission probabilities. I don't care whether the models is learned statistically from a large corpus or from expert knowledge or some combination. I need trained models that I can use. Are there any public sources for such knowledge?


title: Correlated variables in kmeans clustering

text: I have a dataset with 3 variables: A, B and C. Now, A and B are ordinal variables (i.e.; the result of two questions measured using a 5-point Likert), whereas B is continuous. A and B are also correlated, Spearman rho = .50, p-value = 0.0046 I want to partition my dataset in 3 cluster using kmeans (the default R implementation). Does the fact that some of the variables in my dataset are correlated violates any assumptions for running the algorithm?


title: neural network output layer for binary classification

text: I'm using a neural network for a binary classification problem. Is it better to have one neuron in the output layer or to use two, i.e. one for each class?


title: Mixed models with R - convert from SAS code

text: I have this SAS code running a mixed model:               proc mixed data = Data;       class Treatment SubjectID Technique;       model Time = Treatment;       random Technique;       random SubjectID(Technique);       lsmeans Treatment / cl pdiff=control('New') adjust=dunnett;     run;      I want to run this analysis in R, what is the equivalent R code ? Thank you ! (I find the R code for mixed models to be extremely confusing).


title: Is there a branch of statistics that tries to explain "why" the dataset has certain statistical properties?

text: Suppose I have a big dataset and I compute some statistical summary of it - e.g., the correlation of one dimension with another. I think a reasonable question to ask would be "what data points explain this result" - e.g., perhaps, is it because there were two huge outliers that explain the whole correlation? Or are all the points approximately equally important? Or, say, I'm computing the mean of a response variable, measured at two different values of a parameter. I find that the means are approximately the same. Why is it the same: is it because the distributions of the variable are the same, or is it because the mean in both cases is determined by a big cluster that determines the whole mean and obscures the differences in the rest of the distribution? I guess, in general, I'm interested in computing sensitivity of a statistical summary (not necessarily a single number) w.r.t. the data points and parameters involved in the computation. This notion of sensitivity would help me, on one hand, avoid meaningless results (which are entirely explained by abnormal measurements), and on the other hand, it would guide my further exploration. So: is there a branch of statistics that studies this kind of sensitivity? If yes, what are some useful methods from that branch / what further reading would you recommend? E.g., I imagine, it might be useful for data visualization methods - not just draw the data, but colorize each point according to the sensitivity of some metric w.r.t. that point? I tried googling for things like "explanatory statistics", "sensitivity of statistical summaries" etc., but did not find a lot. I found Uncertainty and sensitivity analysis but this is not quite what I'm looking for - I'm not interested in sensitivity of one variable of the dataset w.r.t. the other, I'm interested in making deeper answers to questions like "do these two variables correlate". P.S. Some more googling yielded the keyword "input importance" and some data visualization methods, e.g. http://vis.cs.ucdavis.edu/papers/TVCG_Chan_GSS.pdf . But I crave for advice from experts :)


title: Panel unit root test (IPS test)

text: My question is very straight forward. For a long time, I have been looking for an R package that contains commands on a very popular panel unit root test which assumes cross-section independence, but I have found none. The test is the Im-Pesaran-Shin (IPS) test. Any help will be highly welcome. Thank you.


title: Calculate Average Rate of Change for a Set of Points

text: I have two lists containing X/Y values which are used to plot points on a graph. The X-values are always sequential (1, 2, 3, 4...) but the Y values can vary. For example:               var x = new List<int>(new[] {0, 1, 2, 3, 4, 5, 6, 7});     var y = new List<int>(new[] {0, 3, 7, 12, 14, 15, 17, 19});     var avgY = // What goes here?      With these values (the Y-values being more important in this case), how can I calculate the ~~average rate of change~~ linear regression line so that I can plot another "average" line (or "prediction" line if we were to extend the graph)?


title: K-Means Clustering variable results depending on sort order of data

text: I'm using PASWStatistics18 and K-Means clustering data that I'd previously clustered using an earlier version of SPSS... My problem is that I'm getting different cluster results depending on how the data file is sorted. I can't remember this being an issue in the earlier versions of SPSS, but I can't repeat the same results, using the same data file, again! Only by using the saved cluster centres can I get close to the original result. Have others had similar isues? Am I doing something wrong?


title: Test selection for group comparison

text: We are running a study which investigates platelet aggregation changes during three trimesters of pregnancy and postnatal period. Then, we compare each of these phases to control non-pregnant women. The sample size is 46 including 10 per each trimester, 10 postnatal and 6 as control. The cases we recruited during these phases are not same women (i.e first trimester women are different from those of second trimester and so on). We used nine aggregation agonist for each group. We used Kruskal-Wallis Test to compare these five groups to each other as they are completely different then if we find p-value is significant, we performed PostHoc Mann-Whitney Test to find which two groups are exactly different in aggregation point of view. The questions here are: 1\. Is it the right way to analyse the data in this study using the above tests and are we justified? 2\. What are the limitations of using these tests? 3\. Is there any better way to use in this research? I will be very grateful for any kind of explanation. Thank you very much in advance.


title: Probability that at least one person at a party will accidentally choose their own gift

text: Thirty people are invited to a holiday party, where a gift exchange will take place. Each person brings a gift to put into the pot, and at the end of the night each person selects--at random--one gift from the pot to take home. The probability of one (1) person accidentally choosing their own gift from the pot is equal to $\frac{1}{30} \approx 3.33\%$, obviously. What is the probability that **at least** one person will accidentally choose their own gift? In other words, how unlikely is it that all 30 guests will take home a new present?


title: SEM: Collinearity between two latent variables that are used to predict a third latent variable

text: In some structural equation models that I use in my bachelor thesis, there is a substantial correlation between two latent variables that are used to predict a third latent variable. Now I know there are several ways of quantifying multicollinearity when it comes to observable variables, but what about exogenous latent variables? Can I simply use the same indices that I use for observable variables in regressions analysis (such as the VIF)? And are there any guidelines as to when collinearity between latent variables becomes a problem? I would really appreciate your advice!


title: Robust cluster method for mixed data in R

text: I'm looking to cluster a small data set (64 observations of 4 interval variables and a single three-factor categorical variable). Now, I'm quite new to cluster analysis, but I am aware that there has been considerable progress since the days when hierarchical clustering or k-means were the only available options. In particular, it seems that new methods of model based clustering are available that, as pointed out by chl, enable the use of "goodness-of-fit indices to decide about the number of clusters or classes". However, the standard R package for model based clustering `mclust` apparently will not fit models with mixed data types. The `fpc` model will, but has trouble fitting a model, I suspect because of the non-gaussian nature of the continuous variables. Should I continue with the model-based approach? I'd like to continue to use R if possible. As I see it, I have a few options:   1. Convert the three-level categorical variable into two dummy variables and use `mclust`. I'm unsure if this will bias the results, but if not this is my preferred option.   2. Transform the continuous variables somehow and use the `fpc` package.   3. Use some other R package I haven't yet encountered.   4. Create a dissimilarity matrix using Gower's measure and use traditional hierarchical or relocation cluster techniques. Does the stats.se hivemind have any suggestions here?


title: Two Step estimation of treatment effects

text: I want to estimate the treatment effect of a selection Model without a preimplemented package in R. I allready have the probit estimates for the inverse mills ratio. How can I estimate the switching regression Modell by ols to estimate the treatment effect? Any advice would be great. Thanks.


title: If two time series $X$ and $Z$ follow $0 \leq Z \leq X$, can we say that $\text{var}(Z) \leq \text{var}(X)$?

text: **Now I see it can't hold. Thank you for the counter examples... You guys rule!** **Thank you very much for your comments!** I added, however, some observations that were missing. Most importantly is the fact that we can assume that there exist a positive covariance between X and Y. At first, it seemed to me that it would be easy to demonstrate... but I still did not manage to solve this problem. Can you guys give me a hand? Suppose we make use of $\mathbf{i)}$ a time series $X = [x_1,...,x_N]$ containing only positive entries (i.e. $0 \leq x_i$ for all $i$), $\mathbf{ii)}$ a vector of weights of the same of length given by $Y = [y_1,...,y_N]$ where $0 \leq y_i \leq 1$ for all $i$ to build $\mathbf{iii)}$ a time series $Z = [z_1,...,z_N]$, where the $i$th term is given by $z_i = x_i y_i$, i.e., $Z = [x_1 y_1,...,x_N y_N]$. Clearly, as $Y \in [0,1]$, we have that $0 \leq Z \leq X$ for all $i$. $\mathbf{Question)}$ Can we demonstrate that $\text{var}(Z) \leq \text{var}(X)$? For example, if $X = [2, 6, 99, 12, 3, 1]$ and $Y = [0.34, 0.01, 0.2, 1, 0.3, 0.17]$, we have $ Z = [x_1 y_1,...,x_N y_N] = [0.68, 0.06, 19.8, 12, 0.9, 0.17]$ $ \widehat{\sigma}^{2}_{X} = 1494.70$ $ \widehat{\sigma}^{2}_{Z} = 69.81$ $\mathbf{Important} \text{ } \mathbf{observations}$: 1) $X$ and $Y$ are stationary, ergodic random processes 2) $X$ is not a constant time series, in a sense that $\text{var}(X) \geq 0$ 3) It can be assumed that $\text{var}(X) \geq \text{var}(Y) \geq 0$ 4) There exist a positive covariance $X$ between and $Y$   * Possible implication of 4)? As $0 \leq Z \leq X$, we could define a given time series $W \geq 0$ such as $Z + W = X$. Thus, $\text{var}(X) = \text{var}(Z + W) = \text{var}(Z) + \text{var}(W) + 2\text{cov}(Z,W)$. Note that if $\text{cov}(Z,W) \geq 0$ then $\text{var}(X) \geq \text{var}(Z)$ because $\text{var}(W)$ is also greater than zero. Does the fact that $\text{cov}(X,Y) \geq 0$ infer that $\text{cov}(Z,W) \geq 0$? There is any condition that guarantees $\text{cov}(Z,W) > 0$ Why I was so convinced about $\text{var}(Z) \leq \text{var}(X)$? In the application I am interested in, I have observed that the relation $\text{var}(Z) \leq \text{var}(X)$ is attended at every time I run my algorithm. If I cannot demonstrate that $\text{var}(Z) \leq \text{var}(X)$ holds given the observations 1) to 4), I would like to know what is forcing that relation, like, for example, $\text{cov}(Z,W) \geq 0$ as mentioned above. Thanks again for the replies! Cheers


title: How to find potential classification clues given two sets of classified data

text: I am new to data science and R, and using the following Kaggle data to learn. The input is a matrix with `40` columns of data, for 9.000 lines. The output is a `class`: `0` or `1`. I have split the `train` data into two matrix: one for Class 0, one for Class 1. My intuition tells me to look at the 40 columns and find one (or several) that help best to differentiate Class 0 elements from Class 0 elements. What is/are the statistical technique(s) for this? Should I compare the two distributions, for each column of data. Or how do you recommend to approach such a problem?


title: Multidimensional Scaling - external scale regression using standardized or unstandardized weights

text: I am having trouble determining whether to use standardized or raw regression weights when using multiple regression to interpret Multidimensional Scaling plots. Different textbooks seem to advocate different approaches and I am unable to determine why. In detail, using SPSS I have obtained a MDS ordination for the dataset (249 objects) in 2 dimensions. The stimuli were Likert scale survey responses and I used the mean of the responses along each hypothesised dimension as my external scales. Then I ran linear regression using the external scale as the criterion and the MDS coordinate vectors as predictors. As output I have both raw regression weights (non-zero intercept) and the standardized coefficients with a R value of 0.98 showing good fit. I would like to plot a unit vector of the linear regression line on the MDS plot, so I can determine the angle for rotation. The raw coefficients have a non-zero intercept while the standardized coefficients do not. My question is: Should I use unstandardized (ignoring the constant weight) or standardized coefficients on the MDS plot and why is the other approach wrong? Thanks for any help! I'm sure its a simple answer but explained in contradictory ways in the various textbooks.


title: What is the difference between descriptive and inferential statistics?

text: My understanding was that descriptive statistics quantitatively described features of a data sample, while inferential statistics made inferences about the populations from which samples were drawn. However, the wikipedia page for statistical inference states: > For the most part, statistical inference makes propositions about > populations, using data drawn from the population of interest via some form > of random sampling. The "for the most part" has made me think I perhaps don't properly understand these concepts. Are there examples of inferential statistics that don't make propositions about populations?


title: Similarity algorithms

text: Let's say I have 300 restaurants that I want to compare to each other on the basis of a "similarity score". To try and determine similarity scores, I pick a reference restaurant and pick 3 other restaurants at random and ask a user "Which of these 3 restaurants is most similar to the reference restaurant?". The user selects one of the 3 restaurants and the data is recorded. How many times does this question have to be answered before I can reasonably start to use the data to determine similarity between restaurants?


title: Understanding over-dispersion as it relates to the Poisson and the Neg. Binomial

text: I am developing a Poisson-family glm model in R for a dataset that I have. This dataset has 650 entries with two measures of exposure. The model, though not that relevant to the question, is: $$\ln(E(y_i)) = \ln(\beta_1) + \beta_2 \ln(\text{exp}_1) + \beta_3 \ln(\text{exp}_2)$$ After developing the model in R using a Poisson Distribution, the `summary()` function yields:               Call:     glm(formula = Y ~ log(Vmaj) + log(Vmin), family = poisson(link = "log"),         data = data.raw)          Deviance Residuals:          Min       1Q   Median       3Q      Max           -8.3003  -2.1364  -0.2094   1.8477   5.9722            Coefficients:                      Estimate Std. Error z value Pr(>|z|)             (Intercept) -6.66172    0.24970  -26.68   <2e-16 ***         log(Vmaj)    0.55473    0.02132   26.02   <2e-16 ***         log(Vmin)    0.46638    0.01297   35.96   <2e-16 ***         ---         Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for poisson family taken to be 1)          Null deviance: 7063.6  on 649  degrees of freedom     Residual deviance: 4757.6  on 647  degrees of freedom      Running the same model, but using a negative binomial distribution yields:               Call:         glm.nb(formula = Y ~ log(Vmaj) + log(Vmin), data = data.raw,             init.theta = 5.642678818, link = log)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.3076  -0.8128  -0.0806   0.6368   2.1555            Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept) -6.63440    0.64881 -10.225   <2e-16 ***     log(Vmaj)    0.55507    0.05599   9.914   <2e-16 ***     log(Vmin)    0.46321    0.03258  14.217   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for Negative Binomial(5.6427) family taken to be 1)              Null deviance: 983.84  on 649  degrees of freedom     Residual deviance: 678.02  on 647  degrees of freedom     AIC: 5473.7          Number of Fisher Scoring iterations: 1                    Theta:  5.643            Std. Err.:  0.362            2 x log-likelihood:  -5465.748       One problem that I keep being told about with using a Poisson based distribution is over-dispersion. I understand that a measure of over- dispersion can be found by dividing the Residual (Scaled) deviance by the degrees of freedom. This measure would make the negative binomial look substantial better than the Poisson based distribution for dealing with over- dispersion. I understand that over-dispersion, from the POV of the Poisson, is fundamentally related to the way in which it assumes that the mean = variance. Consequently, I understand that the Negative binomial should have a larger variance, and so the standard error terms should be higher and the z-statistic lower. What I do not understand is why the models have essentially identical coefficients, and are, for all intents and purposes _identical_. My understanding of the process is that the Negative Binomial is a better model, which to me means more accurate. However, when I look at these results I do not understand how using the Negative Binomial has improved the situation, despite the fact that cursory examination based on this shows the Negative Binomial has better addressed the problem over-dispersion.


title: Calculating the probability of being correct due to chance

text: I'm looking for a third opinion on the probability of making a particular decision due to chance alone. I have a task where people read a paragraph of text, then they are shown it a second time, except some of those times one word has changed. People are then asked whether they detected a change. If they say no then that's the end of the task. If they say yes, then they are given a set of nine words and are asked to respond which word was changed. If people correctly guess that a change has occurred AND which word was changed, then a single "correct" is registered. If they correctly detect a change and wrongly guess the word that has changed then an "incorrect" is registered (similarly, it is incorrect if a word has changed but the person did not detect it). What I'm trying to work out is what is the probability of correctly answering that a word has change, AND which word it was, due to chance alone. My first thought was that it was 1/18. That's a 1/2 chance of firstly correctly guessing that a change had occurred, followed by a 1/9 chance of correctly guessing the word. Someone I've spoken to has suggested that the chance rate is actually 1/10. In effect, "no, a change has not occurred" becomes the tenth option. So, is the chance of being correct 1/18, 1/10, or are we both wrong?


title: How to resample samples according to a known distribution

text: I am encountering the following problem: I have been given a number of samples, and my task is to select those from the distribution that we are interested. Suppose the distribution that I am interested in is a Gaussian with a known mean and variance. Then how to only select samples from this Gaussian but not others? I know that I can construct a hypothesis test based on the prior knowledge about this Gaussian, but it does not work that well, especially I have many samples that are actually not from the interested Gaussian distribution, but fall into the confidence interval of the Gaussian distribution. I guess there should be some systematic way to solve this problem. Does anyone has any clue? Actually, when I plot the sample histogram, I can clearly see multiple bumps. Since I know the Gaussian that I am interested, there should be only one bump located at the Gaussian mean, and all other bumps are caused by the samples not from the Gaussian. Is there any method allows me to resample over the given samples and generates a new set of samples following the Gaussian distribution that I am interested in.


title: Subgroup analyses

text: 1200 pregnant patients were asked by the doctor/nurse if they smoke or not. Once the patient admits to smoking she is referred to quit smoking clinic. In the quit smoking clinic they are offered nicotine replacement therapy. The questions are:   1. In the women who use nicotine replacement therapy do their babies do better than ones who refuse nicotine replacement therapy?    2. How do we test statistical significance?


title: Can I test a proportion (k success out of n) against an expected probability obtained from a logistic GLM?

text: I have a logistic regression model which is obtained from literature about the probability of developing a toxicity ($y_{i}=1$ or $0$) when patients are given dose $X_{i}$. $X_{i}$ is a continuous variable. The model can be written as: $Y_{i} \sim B(1,\pi_{i})$ $E(Y_{i})=\pi_{i}$ and $\operatorname{var}(Y_{i})=\pi_{i}(1-\pi_{i})$ $\operatorname{logit}(\pi_{i})=\alpha+\beta X_{i}$ and the values of $\alpha$ and $\beta$ (and their confidence intervals) are known. Now I have a sample of patients (from a different population than the model population) with known $X_{i}$ and $Y_{i}$. In general, I would like to test whether the new patient population behaves significantly different than model above. One method that I can think of is to compare the sample's proportion of toxicity against the expected probability of toxicity computed from the model. Suppose that the proportion of toxicity in the sample is $\hat{\pi}$. Given $X_{i}, i=1,...n$ in the sample, I can compute the expected probability of toxicity as $\frac{\sum\pi_{i}}{n}$. The null hypothesis is $H_0:\hat{\pi}=\frac{\sum\pi_{i}}{n}$ This equals a binomial test. My questions are: 1) Assuming the model (excluding the confidence interval) from the literature is the ground truth, is this a proper test? Do I need to assume that each $X_{i}$ in the population has equal probability? 2) Or I should test $\hat{\pi}$ against the expected probability of toxicity from the theoretical model population? Given equal probability of $X_{i}$, thus equal probability of $\pi_{i}$, $E(\pi_{i})=0.5$, which sounds like I am comparing the sample proportion against 0.5, no matter what the model is. 3) Are there any other methods of comparing the patient population to the given model, e.g. hypothesis based on $\beta$?


title: Assessing Seasonal Covariates in a Seasonal ARIMA Time Series Model

text: I am using R to develop an ARIMA model to evaluate the influence of several seasonal covariates (e.g., meteorological data) upon the incidence of a seasonal disease. I have weekly data available and have set the period equal to 52 weeks. Using the `auto.arima` function the disease of interest has a form of `ARIMA(0,1,2)(0,0,1)`. Using the `TSA` package, I can then evaluate the influence of each of the covariates:               covariates <- data.frame(covariate1, covariate2, covariate3)     model <- arima(disease, order=c(0,1,2), seasonal=list(order=c(0,0,1), period=52), xreg=covariates)      However, I am worried that this approach is identifying spurious associations between the covariates, each of which has a seasonal component... Is it more appropriate to decompose the covariate data and subtract the seasonal component before fitting the ARIMA model?               covariate1_components <- decompose(covariate1)     covariate1_adjust <- covariate1 - covariate1_components$seasonal     [...]     covariates_adjust <- data.frame(covariate1_adjust, covariate2_adjust, covariate3_adjust)     model2 <- arima(disease, order=c(0,1,2), seasonal=list(order=c(0,0,1), period=52), xreg=covariates_adjust)      Any thoughts on which of the two approaches would be preferable for evaluating seasonal covariates?


title: What high level statically typed/compiled languages are good for statistical computing?

text: I have varying levels of experience with Matlab, Mathematica, R, Octave, Python (SciPy) and Julia for statistical/numerical programming. These are all dynamically typed languages and are not compiled (if we ignore JIT compilation). While it is quick and easy to write code in these languages, and some (notably Mathematica) are incredibly powerful, code written in such languages can be brittle to programming errors that can be detected at compile time in statically typed (compiled) languages, and to API changes. Having just been bitten by the latter in some R code I'd written only a year ago (minor changes in the ggplot2 API), I can really see the value in the static and compiled approach. But unless I'm overlooking something obvious, there are no modern high level languages for statistical computing that are also statically typed and compiled. I wouldn't count C/C++/FORTRAN among them. (Perhaps Java or Scala via a good set of libraries?) So, are there statically typed, compiled languages that are particularly good for statistical computing (either in and of themselves, or via a particularly good library)? Thanks.


title: Create variables with pre-specified correlation

text: I have a vector: `a<-rnorm(100)` I would like to create 5 variables that are correlated with `a` to a different extent. For example each of the 5 variables should have a correlation of `.5, .6, .7, .8, .9` with `a` respectively How can this be done?


title: Significance test for *differences of correlations*

text: I have data for two classes A and B for N participants (same group -> dependent) in an experiment. Both classes are described by a feature vector. For each participant i have a certain number of samples (=feature vectors) for each class. Now, in a first step I am interested in the difference of the average within- class sample-to-sample correlation of A (that is, the sample-to-sample reliability of A) and the average between-class sample-to-sample correlation (r_AB; that is, the similarity of A and B) for each participant separately. I could do the same for B, but for simplicity I focus on A. Thus I subtract z-transformed r_AB (=z_AB) from z-transformed r_A (=z_A) for each participant: diffz_A = z_A - z_AB   (let's suppose diffz_A has now N values, i.e. one z-transformed correlation difference for each participant) Now I want to look at two time-points 1 and 2, such that i have two instances of those correlation differences: diffz_A1 and diffz_A2 **How can I test whether the two correlation differences diffz_A1 and diffz_A2 are significantly different?**   Upon literature consultation I found that there is a specific test for testing the difference of two dependent non-overlapping correlations (Steiger et al. 1980). However, it's absolutely unclear how to apply this test to _differences of correlations_ , which in turn are partly based on within-class correlations on the participant level.. So my questions are:   * do i need a complicated z-test in first place for testing between two _differences of correlations_? Or can I use a simple t-test?   * If no, which test should i use and how?   * and somewhat unrelated: if i want to transform a correlation difference back from z-space to r-space (using tanh), can i apply tanh simply to the obtained correlation difference in z-space (i.e. meandiffr = tanh(diffz_A)) or do i have to apply tanh separately to the two terms in the subtraction formulae?


title: Consequences of the Neyman-Pearson lemma?

text: Suppose $X_{1},...X_{n}$ are independently, identically distributed Bernoulli random quantities with parameter $k$.      Consider the hypothesis test:      $H_{0}: k = k_{0}$ vs $H_{1} : k = k_{1} \ \ $where $k_{1} \gt k_{0}$.      Suppose, using the Neyman-Pearson lemma, we obtain a critical region of the form $C^{*} \ = \ \left\\{ (x_{1},...,x_{n} : \sum\limits_{i=1}^{n} x_{i} \geq c) \right\\}$      where c is the critical value, to perform this test such that $P(\sum\limits_{i=1}^{n} x_{i} \geq c | H_{0} \ true) = \alpha$ and the sampling distribution of $\sum\limits_{i=1}^n X_{i}$ is $Bin(n,k)$      My question is, what are the implications (eg advantages) of using the Neyman- Pearson lemma to derive this critical region? I'm not really sure how to go about it so any help would be greatly appreciated.


title: Regression tree algorithm with linear regression models in each leaf

text: **Short version:** I'm looking for an R package that can build decision trees whereas each leaf in the decision tree is a full Linear Regression model. AFAIK, the library `rpart` creates decision trees where the dependent variable is constant in each leaf. Is there another library (or a `rpart` setting I'm not aware of) that can build such trees? **Long version:** I'm looking for an algorithm that builds a decision tree based on a training data set. Each decision in the tree splits the training data set into two parts, according to a condition on one of the independent variables. The root of the tree contains the full data set, and each item in the data set is contained in exactly one leaf node. The algorithm goes like this:   1. Begin with the full dataset, which is the root node of the tree. Pick this node and call it $N$.   2. Create a Linear Regression model on the data in $N$.   3. If $R^2$ of $N$'s linear model is higher than some threshold $\theta_{R^2}$, then we're done with $N$, so mark $N$ as a leaf and jump to step 5.   4. Try $n$ random decisions, and pick the one that yields the best $R^2$ in the subnodes:      * Pick a random independent variable $v_i$, as well as a random threshold $\theta_i$.     * The decision $v_i \leq \theta_i$ splits the data set of $N$ into two new nodes, $\hat{N}$ and $\tilde{N}$.     * Create Linear Regression models on both $\hat{N}$ and $\tilde{N}$, and calculate their $R^2$ (call them $\hat{r}$ and $\tilde{r}$).     * From all those $n$ tuples $(v_i, \theta_i, \hat{r}, \tilde{r})$, select the one with the maximal $min(\hat{r}, \tilde{r})$. This yields a new decision in the tree, and $N$ has two new subnodes $\hat{N}$ and $\tilde{N}$.   5. We have finished processing $N$. Pick a new node $N$ which has not yet been processed and go back to step 2. If all nodes have been processed, the algorithm ends. This will recursively build a decision tree which splits the data into smaller parts, and calculates a Linear Model on each of those parts. Step 3 is the exit condition, which prevents the algorithm from overfitting. Of course, there are other possible exit conditions:   * Exit if $N$'s depth in the tree is above $\theta_{depth}$   * Exit if the data set in $N$ is smaller than $\theta_{data set}$ **Is there such an algorithm in an R package?**


title: YES/ NO ANSWERS

text: ![enter image description here](http://i.stack.imgur.com/GxSgU.jpg)I have a question. I set up a hypothesis: The greater the number of franchise units within a franchised company, the greater the use of the intranet. On the one hand we have the number of franchise units of each franchise company, on the other hand if the company using intranet or not (answers: YES / NO). Could I establish a correlation between date or do some other statistical methods? Can I code Yes as 1, and No as 0? Thank you in advance for your reply.


title: Understanding central limit theorem

text: What is wrong with the following sentence: > The Central Limit Theorem implies that, as the sample size grows, the error > distribution approaches normality. Am I correct by saying that it should in stead state the MEAN of the sample error approaches zero as sample size grows?