title: Estimating variance when only the estimates are available

text: Suppose you have a large set of estimates (predictors) presumably created with a generalized linear model (GLM) with gamma error and an inverse link (possibly a log link). You know that there are 3 variables used in the model: state, type, and territory (zip code groups). You do not know the betas ( _i.e._ , estimated coefficients); you can only estimate them by comparing the estimates. For example, the prediction for territory 1 in a particular state is $x\%$ different from territory 2, so there is a beta that creates that difference in the output. Is there a way--via simulation, closed form, or whatever--to create a bound for the variance of each estimate? I know this is a very strange question, but this does have useful practical applications and your thoughts would be greatly appreciated. **The question again condensed:** You have estimates from a GLM with gamma error and inverse link. You know the estimates vary by state, type and territory. Can you estimate or find an outer bound for the variance of each estimate using any technique you can conjure up?


title: Hypothesis test on variance of normal sample

text: This question is a follow-up to this discussion: http://math.stackexchange.com/questions/634914/distribution-of-likelihood- ratio-in-a-test-on-the-unknown-variance-of-a- normal-s/635371?noredirect=1#comment1339933_635371 Could someone audit the reasoning below? I am trying to derive the distribution of the likelihood ratio statistic for the hypothesis test below. Let $X_1 ... X_{n}$ be a random sample from a $N(\mu,\sigma^2)$ distribution, where **$\mu$ is known and $\sigma^2$ is unknown**. I want to test the hypothesis $H_0 : \sigma^2 = \sigma_{0}^{2} $ vs. $H_1 : \sigma^2 \neq \sigma_{0}^{2}$ (and, trivially, $\sigma^2 >0$). The generic joint pdf for the n independent random variables (ie. the likelihood function for the random sample) is: $L=\prod_{i=1}^{n} \large\left(\frac{1}{\sqrt{2\pi\sigma^2}}\right)\cdot e^{-\frac{(X_i - \mu)^2}{2\sigma^2}}= \large\left(\frac{1}{\sqrt{2\pi\sigma^2}}\right)^{n}\cdot e^{-\frac{\sum_{i=1}^{n} (X_i - \mu)^2}{2\sigma^2}}$ Under the null hypothesis, the maximum value taken by $L$ is: $\large\left(\frac{1}{\sqrt{2\pi\sigma_{0}^{2}}}\right)^{n}\cdot e^{-\frac{\sum_{i=1}^{n} (X_i - \mu)^2}{2\sigma_{0}^{2}}}$ If we do not constrain $\sigma^{2}$ to be equal to $\sigma_{0}^{2}$, then $L$ is maximised by the maximum likelihood estimator $\hat{\sigma}^{2}=\frac{\sum_{i=1}^{n} (X_i - \mu)^2}{n}$ In this case, the maximum likelihood becomes: $\large\left(\frac{1}{\sqrt{2\pi\hat{\sigma_{0}}^{2}}}\right)^{n}\cdot e^{-\frac{\sum_{i=1}^{n} (X_i - \mu)^2}{2\hat{\sigma_{0}}^{2}}}$ Setting these as numerator and denominator, respectively, I get the following likelihood ratio statistic $\Lambda = \LARGE\frac{\large\left(\frac{1}{\sqrt{2\pi\sigma_{0}^{2}}}\right)^{n}\cdot e^{-\frac{\sum_{i=1}^{n} (X_i - \mu)^2}{2\sigma_{0}^{2}}}}{\large(\frac{1}{\sqrt{2\pi\hat{\sigma_{0}}^{2}}})^{n}\cdot e^{-\frac{\sum_{i=1}^{n} (X_i - \mu)^2}{2\hat{\sigma_{0}}^{2}}}} \\\ = \frac{ \left(\frac{1}{\sqrt{2\pi\sigma_0^2}}\right)^{n/2}\cdot e^{-\frac{\sum_{i=1}^n (X_i - \mu)^2}{2\sigma_0^2}}} {\left(\frac{1}{\sqrt{2\pi\hat{\sigma}^{2}}}\right)^{n}\cdot e^{-(n/2)}} = \left(\frac {\hat{\sigma}^2}{\sigma_0^2}\right)^{n/2} \cdot \exp\left \\{-\frac 12\left(\frac{\sum_{i=1}^n (X_i - \mu)^2}{\sigma_0^2}-n\right)\right\\}$ Using the expression for $\hat \sigma^2$ we get $$\Lambda = \left(\frac {\hat{\sigma}^2}{\sigma_0^2}\right)^{n/2} \cdot \exp\left \\{-\frac n2\left(\frac {\hat{\sigma}^2}{\sigma_0^2}-1\right)\right\\}$$ Now, under the null, the random variable denoted $$z_i^2 = \left(\frac {x_i - \mu}{\sigma_0}\right)^2 \sim \chi^2(1)$$ We have $$ \frac {\hat{\sigma}^2}{\sigma_0^2} = \frac 1n\sum_{i=1}^{n} \left(\frac {x_i - \mu}{\sigma_0}\right)^2 = \frac 1n \sum_{i=1}^{n}z_i^2$$ So we can write $$\Lambda = \left(\frac 1n \sum_{i=1}^{n}z_i^2 \right)^{n/2} \cdot \exp\left \\{-\frac n2\left( \frac 1n \sum_{i=1}^{n}z_i^2-1\right)\right\\}$$ Taking minus log we have $$-\ln \Lambda = -\frac n2 \ln \left(\frac 1n \sum_{i=1}^{n}z_i^2 \right) +\frac n2\left( \frac 1n \sum_{i=1}^{n}z_i^2-1\right)$$ Manipulating the second term in the RHS, $$= -\frac n2 \ln \left(\frac 1n \sum_{i=1}^{n}z_i^2 \right) + \sqrt {\frac 1 2} \left( \frac {\sum_{i=1}^{n}(z_i^2-1)}{\sqrt {2}}\right)$$ and multiplying throughout by $\sqrt {2}$ we obtain $$-\sqrt {2} \ln \Lambda = -\frac{n}{\sqrt{2}} \ln \left(\frac 1n \sum_{i=1}^{n}z_i^2 \right) + \left( \frac {\sum_{i=1}^{n}(z_i^2-1)}{\sqrt {2}}\right)$$ $$= -\frac{1}{\sqrt{2}} \ln \left[\left(\frac 1n \sum_{i=1}^{n}z_i^2\right)^n \right] + \left( \frac {\sum_{i=1}^{n}(z_i^2-1)}{\sqrt {2}}\right)$$ The second term in the RHS is a standardized sum of i.i.d $\chi^2(1)$ random variables, and this quantity will converge to a $N(0,1)$. Then (abusing notation a bit), $$\operatorname{plim}\left(-\sqrt {2} \ln \Lambda\right) = \operatorname{plim}\left(-\frac{1}{\sqrt{2}} \right)\cdot \operatorname{plim}\left\\{\ln\left[ \left(\frac 1n \sum_{i=1}^{n}z_i^2 \right)\right]^n \right\\} + \operatorname{plim}\left( \frac {\sum_{i=1}^{n}(z_i^2-1)}{\sqrt {2}}\right) $$ $$= -\frac{1}{\sqrt{2}} \cdot \left\\{\operatorname{plim}\ln\left[\left(1 \right)^n\right] \right\\}+ N(0,1)=-\frac{1}{\sqrt{2}} \cdot\infty + N(0,1) = -\infty$$ At this point, I am completely stuck as I can find neither the test statistic nor its distribution. Could someone flag any mistakes I have made and possibly help me find the test statistic and its distribution?


title: exact calculation of SD on R-squared

text: In R package "psychometrics" an estimate of SE of R squared of               sersq <- sqrt((4*rsq*(1-rsq)^2*(n-k-1)^2)/((n^2-1)*(n+3)))      with n sample size, and k number of parameters if sample size greater than 60 is found. Does anyone have a formula for smaller sample size or an exact formula? I have been through sos - but only found the above? - the real underlying problem is whether comparing models aiming at explaining here a series of concentrations based solely on the R-squared of the individual model predictions vis a vis the measured values is a healthy method


title: Thompson sampling with multivariate posterior distribution

text: I'm implementing Thompson sampling for a multi armed bandit problem (see http://en.wikipedia.org/wiki/Thompson_sampling). The underlying Bayesian model is a Bayesian Linear Regression, which has a conjugate Normal-inverse gamma (NIG) prior for the regression coefficients and the standard error. This NIG distribution is multivariate over the regression coefficients, but I'm not sure how this fits into the Thompson sampling paradigm. Traditionally, you are supposed to sample from the posterior distribution and choose the sample that maximizes the expected reward (i.e. if you are running a beta-binomial model describing the probability of success, then you chose the highest sample from the beta distributions). But, with a multivariate sample, I'm not really sure how to choose the sample that maximizes the expected reward. Wondering if anyone who's familiar with Thompson sampling or multi-armed bandits could help. Thanks!


title: Logistic regression, goodness of fit interpretation

text: I'm having some trouble interpretting whether the model is a good fit. Below is an extract from some output from R.               (Dispersion parameter for binomial family taken to be 1)              Null deviance: 3946.2  on 2999  degrees of freedom     Residual deviance: 3924.3  on 2991  degrees of freedom     AIC: 3942.3          Number of Fisher Scoring iterations: 4          > 1 - pchisq(mylogit$null.deviance - mylogit$deviance, mylogit$df.null - mylogit$df.residual)     [1] 0.004985122      So, as the p-value is low, this would indicate that model is a good fit ie. it is significant right ? Many thanks


title: How to conduct in-sample forecasting?

text: Whilst this may be a fundamental and basic question, I need to pin down a solid understanding of exactly _what_ in-sample forecasting is and entails. I understand out-of-sample as taking a period of training data and forecasting future values based on the regression of data from that period. However, in- sample seems to be seldom defined without reference to not being out-of- sample. Take the example of a time series from 1900-2000 for which I have a regression specification. I want to make in-sample forecasts for 1990-2000. What does this entail?


title: Kernel PCA vs. k-means

text: Could someone compare k-means and kernel PCA in the domain of vector quantization (memory, speed, effectivity, ...)?


title: Equivalence between least squares and MLE in Gaussian model

text: I am new to Machine Learning, and am trying to learn it on my own. Recently I was reading through some lecture notes and had a basic question. Slide 13 says that "Least Square Estimate is same as Maximum Likelihood Estimate under a Gaussian model". It seems like it is something simple, but I am unable to see this. Can someone please explain what is going on here? I am interested in seeing the Math. I will later try to see the probabilistic viewpoint of Ridge and Lasso regression also, so if there are any suggestions that will help me, that will be much appreciated also. Thanks in advance.


title: Examples of when confidence interval and credible interval coincide

text: In the wikipedia article on Credible Interval, it says: > For the case of a single parameter and data that can be summarised in a > single sufficient statistic, it can be shown that the credible interval and > the confidence interval will coincide if the unknown parameter is a location > parameter (i.e. the forward probability function has the form Pr(x | μ) = > f(x − μ) ), with a prior that is a uniform flat distribution;[5] and also if > the unknown parameter is a scale parameter (i.e. the forward probability > function has the form Pr(x | s) = f(x / s) ), with a Jeffreys' prior [5] — > the latter following because taking the logarithm of such a scale parameter > turns it into a location parameter with a uniform distribution. But these > are distinctly special (albeit important) cases; in general no such > equivalence can be made." Could people give specific examples of this? When does the 95% CI actually correspond to "95% chance", thus "violating" the general definition of CI?


title: Whats on the causal path?

text: I have an experiment that perturbs variable x and causes a change in variable z. There is a concurrent change in variable y. How can I determine whether variable y is on the causal path between x and z or if it is an irrelevant epiphenomenon?


title: What is Sigma in Gibbons Ross Shanken (GRS) test?

text: I am trying to figure out what $\Sigma$ means in this GRS formula $F_{GRS} = \frac{T - K - N}{N} \frac{\hat{\alpha}^T \hat{\Sigma}^{-1} \hat{\alpha}} {1 + \hat{\mu}_{K}^T \hat{\Sigma}_{K}^{-1} \hat{\mu}_{K}} $ The goal is to determine whether the alphas are statistically significant from 0. First, we start with a matrix of factors and a matrix of test assets. We regress each excess returns of each test asset on the factors: $r_t^e = \hat{\alpha} + \hat{B} r_{Kt}^{e} + \hat{e_t} $ Alphas from the above regression are the $\hat{\alpha}$. T is the number of observations in each factor (or test asset). K is the number of factors. N is the number of test assets. $ \hat{\Sigma}_{K}^{-1} $ is the covariance matrix or all the factors which is the KxK $ \hat{\mu}_{K} $ is means of each factor which is Kx1. I think I have the above right, but I am not sure what $ \hat{\Sigma} $ is. I think it is either a covariance matrix of test asset returns or covariance matrix of residuals associated with each test asset from the regression equation above. I tried to convince myself that those are the same thing, but I believe that this is the case only when coefficients on all factors are 0. I looked at this thread this and it seems that $ \hat{\Sigma} $ is a covariance matrix of test asset returns, but I am not sure whether I am reading that code properly. Any help is appreciated. Thank You.


title: Estimating nested copula parameters in R

text: Using the R copula package, is there a built-in way to estimate the theta parameters of a nested Archimedean copula (ideally together with the marginals) based on empirical data? In the non-nested case, parameters can be estimated conveniently using:               mydata = read.table("/home/xinli/Desktop/data.csv")         x <- mydata[,1]     y <- mydata[,2]     z <- mydata[,3]     M1 <- mvdc(claytonCopula(5,dim=3),c("norm","norm","norm"),     +      list(list(0,1), list(0,1), list(0,1)))     f1 <- fitMvdc(cbind(x,y,z),M1,       +      start=c(mean(x),sd(x),mean(y),sd(y),mean(z),sd(z),5),      +      method = "Nelder", optim.control=list(trace= TRUE, REPORT= 2))      * * * Given a nested copula, I could not replicate an equivalent outcome:               cop <- onacopula("C", C(5, 1, C(5, c(2,3))))     M2 <- mvdc(cop,c("norm","norm","norm"),list( list(0,1), list(0,1), list(0,1)))          Error in validObject(.Object) :      invalid class “mvdc” object: invalid object for slot "copula" in class      "mvdc": got class "outer_nacopula", should be or extend class "copula"      So, the mvdc constructor accepts copula objects, but my nested copula is of class outer_nacopula. This does not seem to be the right way. Next, let's try estimating only the copula parameters:                f2 <- fitCopula(copula=cop, method="itau", data=pobs(cbind(x,y,z)))           Error in fitCopula.itau(copula, data) :       no slot of name "parameters" for this object of class "outer_nacopula"      Ok, literally the same problem. While previously, the type safety rejected my outer_nacopula, here there does not seem to be a check, but the function clearly expects a copula object again. Moving on to dedicated nested copula estimators:                f3 <- emle(u=pobs(cbind(x,y,z)), cop=cop)           Error: max(cop@comp) == d is not TRUE      At this point, the function checks whether the data dimensionality equals that of the copula. Unfortunately, this is done via cop@comp, ignoring the nested sub-copula that combines the remaining dimensions 2 & 3\. The methods emde and estim.misc have the same problem. This is where my journey currently ends. If anyone can spot my mistake or point me to the correct way of doing this, I would be very grateful.


title: output of a non linear neuron

text: If a neuron uses a non linear activation function such as a sigmoid function, then the output of that neuron can be any value between 0 and 1. suppose if the activation function results in value like 0.6, what will be the output of that neuron? Should it be 0.6 or 1?


title: What is meant by the "level" of a time series?

text: In much of the literature I'm studying it's one of those terms that occurs frequently yet without a rigorous definition to be found. Specifically, I am told: > For time-indexed random variables (RVs) $\\{X_t\\}$, the **additive > decomposition model** is given as > > $$X_t = {ll}(X_{t-1}, X_{t-2}, \ldots) + {fc}(X_{t-1}, X_{t-2}, \ldots, > \varepsilon_t, \varepsilon_{t-1}, \ldots)$$ > > where > >   * $ll$ is the **long-term level** , which is a stochastic process and can > be visualised as a smoothed version of $\\{X_t\\}$, not to be confused with > **trends** which are deterministic patterns >   * $fc$ is the **fluctuation component** which represents changes in > **local level** , assumed **stationary** and with **zero mean level** >   * $\\{\varepsilon_t\\}$ are **innovations** , and are IID mean-zero RVs > But what is the difference in meaning between **trend** vs. **long-term level** vs. **local level** vs. **mean level**? Additionally, aren't the **fluctutation component** and **innovations** modelling the same thing, which is the noise associated with each observation? So why complicate things by including both?


title: How to test for and deal with regression toward the mean?

text: I am working with a large dataset of behavioral data that I am treating (post- hoc) as a time-series experimental design to look for reliable change in a single dependent variable as a result of a treatment. The data comes from user's interaction with a website over 10 years. There is an overall improvement from time 1 to time 2, p<.001. But there is a regression to the mean effect (also see here) such that those with low DV at time 1 increases at time 2, and those with high DV at time 1 decreases at time 2. (This can be seen clearly on graphs.) I don't know how to proceed with the analysis. Can I quantify the regression to the mean effect and from that determine how it is that the treatment effect exceeds the regression effect? Here are some additional details about the study:   * the data for the study come from a support group website where users write about life problems. volunteer counsellors read the entries and respond to the users with support and advice. there is 10 years worth of data; n=~200,000.   * my research is a natural experiment because it works with the website data that was not collected with research in mind. biggest issue with that: no control group.   * the volunteer counsellors also tag (privately) the user's written entry with topic, attribute, and severity labels. I conducted a survey of the counsellors asking them to rate the relative severity of these tags (e.g. 'depression-panic'~0, 'school-worry'~2, 'relationships-happy'~5). The survey results and the tags applied to each writing entry were used to derive a simple proxy for the user's state at time of writing. this proxy was normalized across the sample, has a quasi-normal distribution, and it is treated as the IV.   * the treatment is simply use of the site (writing about life problems & receiving social support), so DV from writing entry 1 to entry 2. the main effect is that IV does increase overall from entry 1 to 2, but as described there is a regression to the mean effect.    * after establishing a main effect, I am interested in looking into a variety of interacting variables: user's language choices, details of website interaction, timing of counsellor response, etc.


title: What is a seasonality of a data stream ? How to use it to estimate a stream attributes ?

text: I am trying to understand what is the rigorous definition of seasonality is, but couldn't find it. Is it a mathematical tool at all? I need to analyze a stream coming from one source and be able to estimate the data coming from the source. I need to take into account the seasonality of the data and try to prototype the source behavior. let's say the stream is `(time, data)` tuples.


title: How to make representative sample set from a large overall dataset?

text: What are the statistical techniques to create a sample set, which is representative of the entire population (with a known confidence level)? Also,   * How to validate, if the sample fits the overall dataset?   * Is it possible, without parsing the entire dataset (which could be billions of records)?


title: Sum of truncated normal distribution

text: I need to estimate what's the probability that the sum of the height of 3 women + 7 men is greater than 17m considering that men follow a normal distribution [mu=1.8 sd=0.10] and women [mu=1.6 sd=0.08], with a catch: I have PRE-selected only men with height lower than 1.9m and women with height lower than 1.75m. Any ideas on how to proceed? Thanks


title: Getting started with biclustering

text: I have been doing some casual internet research on biclusters. (I have read the Wiki article several times.) So far, it seems as if there are few definitions or standard terminology.   1. I was wondering if there were any standard papers or books that anybody who is interested in algorithms for finding biclusters should read.    2. Is it possible to say what is the state of the art in the field? I was intrigued by the notion of finding biclusters using genetic algorithms, so I would appreciate comments on that approach in particular in the context of other approaches.   3. Usually in clustering, the goal is to partition the data-set into groups where each element is in some group. Do bicluster algorithms also seek to put all elements in a particular group?


title: Using an SVM for feature selection

text: Lets say I have a highly dimensional classification problem with a lot of noise, and I want to improve my results by removing some of the noisy variables. I've read several papers on using SVMs for feature selection, but I'm at a loss as to how to implement this in R. Are there pre-existing packages that do this, or am I going to have to roll my own?


title: Incorporating interval-censored event times into standard survival models

text: I am trying to model the time until some event occurs for individuals observed over a 24 month period. For about 75% of people, no event occurs. For 15% of people, we know exact time of the event. For the other 10%, we only know a time window in which it occurs. So some sample "survival times" might be:               >24      2.5      >24      5.0      0 to 6     6 to 12      >24      18 to 24       and so on. Can this type of data be accommodated within "standard" survival models? If so, some hint at how to construct the `Surv()` object in `R` and how to pass it to `survreg` would be appreciated.


title: Can full conditionals determine the joint distribution?

text: I heard that all the full conditionals (as used in Gibbs sampling) can determine the joint distribution. But I don't understand why and how. Or did I mis-hear? Thanks!


title: Estimate Xmax for power law heads

text: poweRlaw is able to estimate the $x_{min}$ and $\alpha$ for power law tails, but I need to estimate the $x_{max}$ and $\alpha$ for a power law head. Can I take the reciprocal for each data point, and estimate $x_{min}$ and $\alpha$ for the reciprocal distribution? Is there any issue with this approach?


title: Clarification about normal distribution issues

text: I'm doing the following self-help revision question on Normal distribution. I don't quite seem to understand:   * the answer in Part II and how it is derived   * for Part III, why is is that p = 0.1075    * for Part III, why does P( Y = 5 ) not have a multiplication with (1 - 0.1075) like P( Y = 4 ) Appreciate some advice and directions please. **Question** ![enter image description here](http://i.stack.imgur.com/tOHW6.jpg) **Answer Key** ![enter image description here](http://i.stack.imgur.com/OJR0q.jpg)


title: Power to detect significant correlation given measurement uncertainty

text: First of all, apologies if some of the technical terminology in here is not accurate - my experience with statistics is fairly limited. To illustrate my problem, let's say I wanted to design a study to calculate a linear correlation between IQ and head circumference (HC) as measured with a tape. The human makes a measurement error of 0.1 cm on average, and errors are randomly distributed. How do I determine how many subjects need to be enrolled in the study in order for me to detect a statistically significant correlation between IQ and HC? Is it possible to determine this a priori, or is this dependent on how strong the correlation turns out to be? Thanks for clearing this up!


title: Lévy stable vs. extreme value distributions

text: I'm trying to understand the advantages (if any) of employing the Generalized Extreme Value distribution (GEV) vs. a stable distribution in the context of understanding the probability of crossing a large threshold. The stable is limiting for averages, whereas the GEV is limiting for extremes. Both distributions have heavy-tailed forms (relevant for crossing high threshold). But it seems to me that the one should employ the GEV because 1) it focuses on the largest values, which are of interest; 2) it can describe data w/ more types of tailedness (can have a bounded upper tail, while the stable can only be heavy or thin). **Are these differences correct?** Assume large values and thresholds of are interest, and someone has analyzed a time series and found that it follows a Lévy-stable distribution (i.e., is heavy-tailed). **Is it sufficient (or satisfactory) to to draw conclusions about the probability of observing a certain extreme value using the stable, or is there more information to be gained by using the GEV?**


title: Is it necessary to perform a transformation on proportion data if it's reasonably well behaved?

text: If you're modelling a proportion response against numerous predictors that are also proportions, is it necessary to transform the response if the standard OLS model is seemingly well behaved? By well behaved I mean:   * None of the fitted values are outside the range [0,1] (In fact they are fairly accurate)    * Residuals look good I believe arcsine transform is typically used in this scenario to make the data look normal, but what if this is not needed? Also, say the data wasn't normal, would a transform still be necessary if one were modelling the proportions with the Random Forest technique? Cheers


title: Expected value of a transformation

text: A simple question. If $Y=\frac{1}{X}$ and I know $f_X(x)$, is it true that $E(Y) = E(1/X) = \int_{-\infty}^\infty \frac{1}{x}f_X(x) dx$?


title: What kind of distribution is $f_X(x) = 2 \lambda \pi x e^{-\lambda \pi x ^2}$?

text: What kind of function is: $f_X(x) = 2 \lambda \pi x e^{-\lambda \pi x ^2}$ Is this a common distribution? I am trying to find a confidence interval of $\lambda$ using the estimator $\hat{\lambda}=\frac{n}{\pi \sum^n_{i=1} X^2_i}$ and I am struggling to prove if this estimator has Asymptotic Normality. Thanks


title: Laplace rule of succession for a Poisson process

text: Suppose that I measure for an interval of time $t$ a Poisson process that has unknown average rate $\mu>0$ and I observe no event. What could be a reasonable estimate of $\mu$ from a Bayesian point of view?


title: What are the assumptions (and H0) for Wilcoxon signed-rank test?

text: I am working on the assumption page of the Wilcoxon signed-rank test (the wilcox "paired test") in Wikipedia. I was able to locate a reference for the assumption that I wrote there, which are: Let $Z_i=X_i – Y_i$ for $i = 1, \ldots , n$.   1. The differences $Z_i$ are assumed to be independent.    2. Each $Z_i$ comes from the same continuous population, and is symmetric about a common median θ .   3. The values which $X_i$ and $Y_i$ represent are ordered (at least the ordinal level of measurement), so the comparisons "greater than", "less than", and "equal to" are useful. However, I get this nagging feeling that I am missing something here. Can anyone correct/expand on this? Thanks.


title: Simpson's paradox or confounding?

text: Consider a scenario where a two-way contingency table is analyzed by a chi- squared test of independence and a significant result is found. Now, it turns out that this table is an aggregation of data from two subgroups which are heterogenous and one much bigger than the other. When analysed at the sub- group level, both groups give a non-significant result. How is this best explained ? Edit: I managed to find the paper where I read about this. They explain it as Simpson's Paradox. http://www.amstat.org/publications/jse/secure/v7n3/datasets.morrell.cfm


title: Should I run measurement invariance if confirmatory factor analysis is not affirmative?

text: I have a repeated measures design for a Likert scale with 10 items. This scale is assumed to be unidimensional and this is one of the research questions along with the scales measurement invariance (MI) over time. I conduct the analysis as follows: First assessment is randomly split into two pieces, I apply PCA in the one half, CFA in the second half and then I wish to make MI comparisons between the second half of the first assessment and the second assessment. PCA is supporting a one factor construct but not very strongly (λ = 2.7, Prop: 73%) while CFA is problematic to provide one factor structure (c2/df > 3, RMSEA = 0.12, GFI = 0.87). Under that circumstances does it make sense to run and provide also the necessary MI comparisons between the two assessments?


title: Probability that items chosen randomly are defective?

text: 6 light bulbs are chosen at random from 17 bulbs of which 6 are defective. (a) What is the probability that exactly 2 are defective? (b) What is the probability that at most 1 is defective?


title: Appropriate sample statistics to provide a reliable estimator for an Erlang-like parent population

text: We would like to measure (in near real-time) the completion latency of requests to our queue-based system. The volume of inbound requests is sufficiently high that we would like to sample/aggregate and report sample statistics back on regular (3 second basis). Our problem is that we don't know what statistics to capture from the sample and how to combine them to provide a reliable set of estimators for the parent population. We know that the distribution of the parent population is not normal, we suspect the closest is a gamma function, e.g. the Erlang distribution.


title: How the Pearson's Chi Squared Test works

text: Following a recent down vote I have been trying to check my understanding of the Pearson Chi Squared test. I usually use the chi squared statistic (or reduced chi squared statistic) for fitting or checking the resulting fit. In this case the variance is not usually the expected number of counts in a table or histogram but some experimentally determined variance. Either way, I was always under the impression that the test still used the asymptotic normality of the multinomial PDF (i.e. my test statistic is $$Q = (n-Nm)^\top V^{-1}(n-Nm)$$ and $(n-Nm)$ is asymptotically multinormal where $V$ is is covariance matrix). Therefore $Q$ has a chi-squared distribution given large $n$ so using the expected number of counts as the denominator in the statistic becomes valid for large $n$. Its possible that this is only true for histograms, I haven't analysed a small table of data in years. Is there a more subtle argument that I am missing? I would be interested in a reference, or even better a short explanation. (Although its possible I just got down voted for omitting the word asymptotic, which I concede is rather important.)


title: Choosing regression model

text: How can I choose between these two regression models? `R` outputs: **Regression 1**               Call:     lm(formula = log.h ~ year)          Residuals:          Min       1Q   Median       3Q      Max      -1.24004 -0.45221 -0.05301  0.44744  1.42060           Coefficients:                   Estimate Std. Error t value Pr(>|t|)         (Intercept) -1.588e+02  1.565e+01  -10.15 9.47e-13 ***     year         8.809e-02  7.859e-03   11.21 4.65e-14 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.6396 on 41 degrees of freedom     Multiple R-squared: 0.7539, Adjusted R-squared: 0.7479      F-statistic: 125.6 on 1 and 41 DF,  p-value: 4.654e-14       ![Residuals vs predicted values](http://i.stack.imgur.com/aVVst.jpg) ![Observed values vs predicted values](http://i.stack.imgur.com/IFWgm.jpg) * * * **Regression 2:**               Call:     lm(formula = log.h ~ year + I(year^2) + I(year^3))          Residuals:         Min      1Q  Median      3Q     Max      -1.4044 -0.3966  0.1141  0.3447  1.2321           Coefficients:                   Estimate Std. Error t value Pr(>|t|)       (Intercept) -1.075e+06  4.663e+05  -2.304   0.0266 *     year         1.622e+03  7.027e+02   2.308   0.0264 *     I(year^2)   -8.162e-01  3.529e-01  -2.313   0.0261 *     I(year^3)    1.369e-04  5.909e-05   2.317   0.0259 *     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           Residual standard error: 0.58 on 39 degrees of freedom     Multiple R-squared: 0.8075, Adjusted R-squared: 0.7927      F-statistic: 54.54 on 3 and 39 DF,  p-value: 5.093e-14       ![enter image description here](http://i.stack.imgur.com/Svo88.jpg) ![enter image description here](http://i.stack.imgur.com/JvUaE.jpg) **Thoughts:** The point is the first model is good, $R^{2}$ is good, p-values are great, but there are hints of heteroskedasticity of residuals vs predicted values. The second one has a better adjusted $R^{2}$, a great $p$-value of regression, and the scatterplots have got better (not drastically, but better). On the other hand, the $p$-values of the coefficients are a lot higher than before. Considering all this, how can I choose between them?


title: Confusion related to conditional distribution from latent variables and observations

text: I am a bit confused about how this expression was derived where Z are the latent variables and X the observations. ![enter image description here](http://i.stack.imgur.com/RU4pl.png) I think it is something like this but just to be sure $p(z_{nk}|X,Z_{-nk}) = \frac{p(z_{nk}, X, Z_{-nk})}{p(X,Z_{-nk})}$ $= \frac{p(X|z_{nk},Z_{-nk})*p(z_{nk}|p(Z_{-nk})) * p(Z_{-nk})}{p(X|Z_{-nk})*p(Z_{-nk})}$ $\propto p(X|Z)*p(z_{nk}|Z_{-nk})$ Is it correct? My question is don't I need to set $z_{nk} = 1$ inside $p(X|Z)$ since $z_{nk}$ is a part of it? ![enter image description here](http://i.stack.imgur.com/1wDgB.png)


title: Algorithm for determining performance speedup/slowdown in a code change vs. historical data?

text: I work on a large programming team, and I run a suite of performance tests on every change that is made in our program, which basically measure time it takes to run the test. For every code change, we run these tests, and we calculate whether the change caused the test to run slower by doing a two- sample t test (against the results from the previous code change). This works decently, but the problem is that we only have a small number of sample data points, generally 5 per test, per code change. There are about 400 individual measurements that we track performance on, so we see some noise in our results (i.e. the t-test will yield a small p value for tests which are not actually any faster/slower due to the code change). Even though we have a small number of sample points on each code change, we have a very large history of results. I want to use this historical data to help us, but I'm not sure how I can. A problem I'm worried about is that any code change may cause the tests to run faster or slower, so just blindly aggregating historical data will yield a poor result. Are there any statistical tests that will help me with this? For a little more info: Most of the time code changes do not have any impact on performance, and those of them that do cause performance tests to run faster or slower only do so on a handful of the 400 tests. Which means that for any given test, it could be hundreds of code changes before a change actually makes the test run faster or slower. To clarify, I want to figure out when a code change actually causes a test to run faster or slower. What options do I have?


title: Categorization scheme of regression techniques

text: There exist many techniques that can be used for regression task. Some are from statistics (e.g. generalized least squares) and some from machine learning (e.g. MLP), some are linear and some nonlinear, some are parametric and some nonparametric. Are there available any categorization schemes of regression techniques?


title: Distribution check using R

text: I want to check the probability distribution for my data using R.In some site I found `fitdistr` function in package `MASS` can do that. To check that I have generated 105 random Poisson numbers & run the `fitdistr` function to check whether it is able to identify or not. I used following code               library(MASS)     zpois=rpois(105,0.1)          fitdistr(zpois, 't')$loglik     fitdistr(zpois, 'normal')$loglik     fitdistr(zpois, 'logistic')$loglik     fitdistr(zpois, 'weibull')$loglik     fitdistr(zpois, 'gamma')$loglik     fitdistr(zpois, 'lognormal')$loglik     fitdistr(zpois, 'exponential')$loglik     fitdistr(zpois, 'Poisson')$loglik     fitdistr(zpois, 'negative binomial')$loglik      I found that it is giving lowest value for normal distribution. I know that the large sample approximation of Poisson distribution is normal distribution but I don't want the large sample approximation. I want to see the exact distribution. Can you help me to guide the suitable function in R using which I can get the exact distribution fit?


title: Strict exogeneity and lagged variables

text: I am confused why strict exogeneity must be violated when we have lagged time series variables. My understanding of strict exogeneity is that a variable must be uncorrelated with error terms in all periods. But isn't exogeneity always a necessary assumption for estimation? If $x_t$ and $u_t$ are uncorrelated, and $x_{t-1}$ and $u_{t-1}$ are uncorrelated, how would it violate strict exogeneity if we have a specification with $x_t$ and $x_{t-1}$?


title: Is Fisher overestimating the variance here?

text: [NB: the title of this post used to be "Please help me decipher Fisher on ANOVA". I've edited the question to make it more focused, and modified the title accordingly.] Since anything written by R. A. Fisher seems to be fertile ground for misunderstanding, I will quote him at length. At the beginning of ch VII ("Intraclass correlations and the analysis of variance"), on pp. 213-214 of his _Statistical methods for research workers_ (1973, 14th ed. rev. & enl.), Fisher writes: >     If we have measurements of $n^\prime$ pairs of brothers, we may > ascertain the correlation between brothers in two slightly different ways. > In the first place we may divide the brothers into two classes, as for > instance elder brother and younger brother, and find the correlation between > these two classes exactly as we do with parent and child. If we proceed in > this manner we shall find the mean of the measurements of the elder > brothers, and separately that of the younger brothers. Equally the standard > deviations about the mean are found separately for the two classes. The > correlation so obtained, being that between two classes of measurements, is > termed for distinctness an **interclass** correlation. Such a procedure > would be imperative if the quantities to be correlated were, for example, > the _ages_ , or some characteristic sensibly dependent upon age, at a fixed > date. On the other hand, we may not know, in each case, which measurement > belongs to the elder and which to the younger brother, or, such a > distinction may be quite irrelevant to our purpose ; in these cases it is > usual to use a common mean derived from all the measurements, and a common > standard deviation about that mean. If $x_1,\; {x^\prime}_{\\!1}; \;\; > x_2,\;{x^\prime}_{\\!2}; \;\; \cdots > ;\;\;x_{n^\prime},\;{x^\prime}_{\\!n^\prime}$ are the pairs of measurements > given, we calculate $$\begin{array} & & \\\ \overline{x} & = & > \frac{1}{2n^\prime} \mathrm{S}(x + x^\prime), \\\ s^2 & = & > \frac{1}{2n}\\{\mathrm{S}(x - \overline{x})^2 + \mathrm{S}(x^\prime - > \overline{x})^2 \\}, \\\ r & = & \frac{1}{ns^2} > \mathrm{S}\\{(x-\overline{x})(x^\prime - \overline{x})\\}. \end{array}$$ > When this is done, $r$ is distinguished as an **intraclass** correlation, > since we have treated all the brothers as belonging to the same class, and > having the same mean and standard deviation. (…) (I've made the point to reproduce Fisher's words, punctuation, typography, and notation _exactly_. You may take this to mean that a "sic" applies to the entire passage quoted above, including all the mathematical expressions.) If I were to _attempt_ to interpret these expressions using more modern notation (and displaying a modicum of elementary courtesy towards the reader), I'd write $$\begin{array} & & \\\ \overline{x} & = & \frac{1}{2n^\prime} \sum_{i=1}^{n^\prime}(x_i + {x^\prime}_{\\!i}), \\\ s^2 & = & \frac{1}{2n} \left\\{ \sum_{i=1}^{n^\prime}(x_i - \overline{x})^2 + \sum_{i=1}^{n^\prime}({x^\prime}_{\\!i} - \overline{x})^2 \right\\}, \\\ r & = & \frac{1}{ns^2} \sum_{i=1}^{n^\prime}(x_i-\overline{x})({x^\prime}_{\\!i} - \overline{x}). \end{array}$$ Now, if we _assume_ that $n = n^\prime - 1$, then the last two definitions become $$\begin{array} & & \\\ s^2 & = & \frac{1}{2n^\prime - 2} \left\\{ \sum_{i=1}^{n^\prime}(x_i - \overline{x})^2 + \sum_{i=1}^{n^\prime}({x^\prime}_{\\!i} - \overline{x})^2 \right\\}, \\\ r & = & \frac{1}{(n^\prime - 1)s^2} \sum_{i=1}^{n^\prime}(x_i-\overline{x})({x^\prime}_{\\!i} - \overline{x}). \end{array}$$ In this form, the definition of $r$ matches that of a sample correlation coefficient for the case in which both variables have the same estimated mean $\overline{x}$ and estimated variance $s^2$. But I'm still thrown off by the expression for this estimated variance $s^2$. According to Fisher's description $s^2$ should be the square of "a _common_ standard deviation about [a _common_ ] mean". Therefore, I expected that the denominator in its pre-factor would be $2n^\prime - 1$. IOW, I expected the whole expression would be $$ s^2 = \frac{1}{2n^\prime - 1} \left\\{ \sum_{i=1}^{n^\prime}(x_i - \overline{x})^2 + \sum_{i=1}^{n^\prime}({x^\prime}_{\\!i} - \overline{x})^2 \right\\}\,. $$ But if $n = n^\prime - 1$, then Fisher uses a denominator of $2n = 2n^\prime - 2$. Does anyone know why? Since the quote above comes from a 14th revised edition, it seems to me unlikely that this is merely a typo. Furthermore, in the following page (p. 215) he extends the discussion to the the case of "trios" of brothers, and there he gives the expression: $$ s^2 = \frac{1}{3n}\\{\mathrm{S}(x - \overline{x})^2 + \mathrm{S}(x^\prime - \overline{x})^2 + \mathrm{S}(x^{\prime\prime} - \overline{x})^2\\}\,, $$ or, in my notation, and still assuming that $n = n^\prime - 1$, $$ s^2 = \frac{1}{3n^\prime - 3} \left\\{ \sum_{i=1}^{n^\prime}(x_i - \overline{x})^2 + \sum_{i=1}^{n^\prime}({x^\prime}_{\\!i} - \overline{x})^2 + \sum_{i=1}^{n^\prime}({x^{\prime\prime}}_{\\!i} - \overline{x})^2 \right\\}\,. $$ Therefore, if we're dealing with typos here, they are _matched typos_. Of course, for a sufficiently large $n^\prime$, the difference between a denominator of $2n^\prime - 1$ and one of $2n^\prime - 2$ is negligible, but my understanding is that Fisher is explicitly assuming that we are working with values of $n^\prime$ that are in general too small to justify this approximation. In fact, if the approximation $n^\prime \approx n^\prime - 1$ were valid, Fisher's notational distinction between $n^\prime$ and $n = n^\prime - 1$ would be hard to understand.


title: ANOVA / Kruskal-Wallis: one of four groups has different distribution

text: I'm using Kruskal-Wallis to as a non-parametric ANOVA of a non-normal distribution. My understanding is that it assumes each group of the independent variable has the same shape. As it turns out, three of my four groups have approximately the same shape, however the fourth does not. The fourth is roughly the inverse of the first two. Does this invalidate significance? Are there any steps I should take to mitigate its impact? Update based on comment and further steps taken: @Miroslav Sabo pointed me to a similar question he asked about using Mann- Whitney when certain assumptions are not met. That answer suggested welch's correction.


title: Best way to present a random forest in a publication?

text: I am using the random forest algorithm as a robust classifier of two groups in a microarray study with 1000s of features.   * What is the best way to present the random forest so that there is enough information to make it reproducible in a paper?   * Is there a plot method in R to actually plot the tree, if there are a small number of features?   * Is the OOB estimate of error rate the best statistic to quote? Many thanks


How outliers can effect the correlation results (Spearman and Pearson)? I have two series of data. The source of the data are different but both are temperature that are measured by different instruments with different accuracy. Please possibly give answers with reference.


A little background about my data used. I have 15 predictors, 14 are categorical and one is continuous. One with 4 levels, another one with 3 levels and rest are with two levels each, 1 or 2. The response variable in my model is a four-level categorical variable. I used the last level as the reference category. I tried combining the levels of the response variable and predictors and also excluded some of the predictors but still get a warning: `"Unexpected singularities in the Hessian matrix are encountered. This indicates that either some predictor variables should be excluded or some categories should be merged"`, with Standard Error >2. Could anyone please tell me what I should do to solve this problem?


title: What's the bayesian equivalent of a general goodness of fit test?

text: I have two data sets, one from a set of physical observations (temperatures), and one from an ensemble of numerical models. I'm doing a perfect-model analysis, assuming the model ensemble represents a true, independent sample, and checking to see if the observations are drawn from that distribution. The statistic I've calculated is normalised, and should theoretically be a standard normal distribution. Of course it's not perfect, so I want to test for goodness of fit. Using frequentist reasoning, I could calculate a Cramér-von Mises statistic (or Kolmogorov-Smirnov, etc.), or similar, and look up the value in a table to get a p-value, to help me decide how unlikely the value I see is, given the observations are the same as the model. What would the Bayesian equivalent of this process be? That is, how do I quantify the strength of my belief that these two distributions (my calculated statistic and the standard normal) are different?


title: What is the analysis techniques for comparing two sets of ratings

text: Say, I have collected 50 comments. For each comment, I am planning to get 10 ratings from non-expert raters for the following criteria: "interesting vs. not interesting". Based on the plurality of those ratings, I will determine if a comment is "interesting" or not. For example, 6 out of 10 raters rated comment #1 as "interesting", I assume comment #1 is interesting The thing that complicates the matter a bit is that I also have experts (the ones who are trained/versed in giving such ratings) rate those comments in a similar fashion as mentioned above. My question is: I would also like to know if I can use non-experts (less costly) to rate the comments instead of paying more for experts. What kind of statistical/mathematical calculation can I do to find that out? Is there a calculation that tells us how much of two sets of ratings overlap/agree? Thank you very much.


title: Trees generated by multi-type branching processes in n steps

text: I am trying to develop some algorithm to compute probabilities in multi-type branching trees, and I doubt I am doing this right... Let us consider a multi-type branching process with two types, denoted by 0 and 1. The process starts in state 0 with probability 1, so that the root vertex of any tree generated by this process has state 0. A vertex in state 0 generates two vertices in state 1 with probability 1. A vertex in state 1 generates either two vertices in state 0 (with probability 0.5) or one vertex in state 0 and one vertex in state 1 (with probability 0.5). We denote those three possible transitions by:   * A: 0 -> 1 1 (1.)   * B: 1 -> 0 0 (0.5)   * C: 1 -> 0 1 (0.5) (their probabilities are between brackets). In the sequel, a "tree" will refer to an unordered tree, _i.e._ two isomorphic ordered trees such that the tree isomorphism preserves the labels will be considered as a same mathematical object (following the principle in Chi (2004), p. 1993, paragraph 3 - see link at the end of the post). Now let us consider the particular labeled tree of height three, composed by one root vertex in state 0 with   * one first child in state 1 that itself has one child in state 0 and one child in state 1   * one second child in state 1 that itself has two children in state 0 This tree should be depicted in some file I am not allowed to post, so I should draw it in ascii mode with limited guarantee on the result:                          0                  |            -------------            |           |            1           1          |           |        ---------   ---------        |       |   |       |        0       1   0       0        I would like to compute the probability that this tree was generated by the above multi-type branching process after 2 generation steps (the generation of the root vertex does not count as a step). Using equation 2 p. 1994 in Chi (2004), this probability should be 0.25, since each transition A, B and C is applied once. However, each possible tree of height 3 among the three trees generated by this process has probability 0.25, and the sum of the probabilities is 0.75 instead of 1. Another possibility is to consider that in this given tree, the set of children (0, 0) of vertex 1 may have been generated by any vertex in state 1 (and same principle for the set of children (0, 1)), so that the tree probability is in fact 0.25+0.25 = 0.5. Finally, how to compute the probability that a given tree _t_ of height _n_ was generated by a multi-type branching process after _n-1_ generation steps ? Can equation 2 p. 1994 in Chi (2004) be used ? Or do we have to compute the number of trees that are isomorphic to _t_ in some sense ? Or do we have to give up the idea that isomorphic trees are equivalent representations of same object ? Thanks for your help ! JB Ref.   Z. Chi. Limit laws of estimators for critical multi-type Galton–Watson processes. Ann. Appl. Probab. Volume 14, Number 4 (2004), 1992-2015. * * * As a beginning of answer to the suggestion of Aniko, I think the following equivalence relation could lead to equivalence classes of trees such that formula (2) could be applied. It consists in adding constraints to the usual notion of isomorphism (by "isomorphism", I mean "an isomorphism preserving labels"). I denote by $I(v)$ the label of vertex $v$, and $t_v$ the complete subtree rooted in $v$. A complete subtree rooted in $v$ is the subgraph induced by the descendants of $v$. The property $P(t)$ is defined inductively by: $P(t)$ is satisfied iif $t$ has a single vertex, or if for every pair ${v, v'}$ of children of the root of $t$ such that $I(v)$ = $I(v')$, $t_v$ and $t_{v'}$ are isomorphic and they satisfy $P(t_{v'})$ and $P(t_{v})$. The two trees $t$ and $t'$ are in the same equivalence class iif [they are isomorphic, and property $P$ is satisfied for $t$ and $t'$] or [$t=t'$].


title: Determining regression tree quality

text: I'm attempting to use the `rpart` R package, and I'm having difficulty figuring out how to determine the quality of a given tree output. For most linear models I would just examine p-values and $r^2$ values to determine whether the model performs satisfactorally. Is there a similar number for decision trees, or is the only performance metric available how well it can fit the data? (Note: I'm trying to use a decision tree to fit to a continuous dataset. I guess a very relevant separate question is, whether a decision tree is appropriate for that type of data?)


title: calculating recall and precision confusion

text: I got 100 text, 9 labeled as + and 91 labeled as -. My classifier predicted, 45 as + and 55 as - So I have been asking myself this question. Out of 45+ predicted text, how many are actually +: answer = 9. Out of 55- predicted text, how many are actually -: answer = 3. Precision = 9/12 Recall = 9/45 Are these corrects? I ask because I have looked at the wiki and other websites for examples but their predicted values tend to be less than the labeled values and in my case it is the opposite so I'm tearing my hairs out as how to calculate them. Thanks. Edit: After rechecking my figures, this is my confusion matrix:                                 Labeled                      | +  | -  |                      -----------        Predicted | + | 9  | 36 |                  | - | 0  | 55 |      Edit: Realised my mistake, I have flipped the confusion matrix wrong way around... So now recall= 9/9 = 1 and precision= 9/45 = 0.2


title: Metropolis-Hastings algorithms used in practice

text: I was reading Christian Robert's Blog today and quite liked the new Metropolis-Hastings algorithm he was discussing. It seemed simple and easy to implement. Whenever I code up MCMC, I tend to stick with very basic MH algorithms, such as independent moves or random walks on the log scale. What MH algorithms do people routinely use? In particular:   * Why do you use them?    * In some sense you must think that they are optimal - after all you use them routinely! So how do you judge optimality: ease-of-coding, convergence, ...  I'm particularly interested in what is used in practice, i.e. when you code up your own schemes.


title: Doubts regarding neural network training using k-fold cross validation?

text: Facing the following problem: Using k-fold cross validation to select network parameters. Early stopping on the held-out fold to obtain the 'best' weights in terms of classification error. Repeat for k-folds. Use lowest average error (from these k runs) to select best network. Now I have selected the best network but how do I train it on the whole training set for deployment? In the k runs, best 'weights' were obtained at different iterations. Do I use a separate test set as stopping criteria ? If so would this indicate fairly the generalization performance?


title: Does caret train function for glmnet cross-validate for both alpha and lambda?

text: Does the R `caret` package cross-validate over both `alpha` and `lambda` for the `glmnet` model? Running this code,               eGrid <- expand.grid(.alpha = (1:10) * 0.1,                           .lambda = (1:10) * 0.1)          Control <- trainControl(method = "repeatedcv",repeats = 3,verboseIter =TRUE)          netFit <- train(x =train_features, y = y_train,               method = "glmnet",               tuneGrid = eGrid,               trControl = Control)      The training log looks like this.               Fold10.Rep3: alpha=1.0, lambda=NA       What does `lambda=NA` mean?


title: Criteria to set STL s.window width

text: Using `R` to perform STL decomposition, `s.window` controls how rapidly the seasonal component can change. Small values allow more rapid change. Setting the seasonal window to be infinite is equivalent to forcing the seasonal component to be periodic (i.e., identical across years). **My questions:**   1. If I have a monthly time series (that is frequency equal to $12$), what criteria should be used to set `s.window`?   2. Is there any link between that and the time series frequency?


title: Mahout Scability

text: Do you know any real world examples of how much Mahout can scale? I wonder how much it can scale in collaborative filtering, clustering, and classification ?


title: Fitting a distribution with ties in the data

text: I'd like to fit some standard distributions, say lognormal, to a set of data which unfortunately has ties. The precision of the measurement simply isn't high enough. Of course you can still fit the data but all goodness-of-fit test like KS and the like throw warnings - which is natural, since this case cannot happen for continuous distributions. I tried to add some arbitrary small random $\epsilon$ to every data point to remove the ties without altering the data too much. As expected this results in ridiculously low p-values. On the other hand: The data consists of some measurements per year. Taking only the average value/year one can fit a distribution very well. So I think it's really because of the ties. Is there any magic solution to this problem?


title: Simulation of multivariate correlated variables

text: I am trying to recreate the simulation method used in the paper _Relaxing the Rule of Ten Events per Variable in Logistic and Cox Regression_ by Vittinghoff and McCulloch. The setup is the following: We are dealing with one (primary) binary predictor/covariate $X_1$ and several other continuous adjustment variables/covariates $X_2,\ldots,X_n$. These continuous variables follow a multivariate normal distribution with mean zero, unit variance and pairwise correlation of $0.25$, i.e. $$ (X_2,\ldots,X_n)\sim \mathcal{N}_{n-1}(\mathbf{0},\Sigma), $$ where $\Sigma$ is such that $\mathrm{Corr}(X_i,X_j)=0.25$ for $i,j=2,\ldots,n$ and $i\neq j$ and such that $\mathrm{Var}(X_i)=1$ for $i=2,\ldots,n$. Now they're saying the following: > The binary primary predictor was generated with expected prevalence of > $0.1$, $0.25$, or $0.5$ and multiple correlation with the covariates of $0$, > $0.25$, or $0.75$. My interpretation of this is that $X_1\sim b(1,p)$, where $p$ is varied through $0.1$, $0.25$ and $0.5$. But I do not understand the part with the multiple correlation. Can anyone elaborate on this? What does it mean mathematically and how do I simulate such $X_1,X_2,\ldots,X_n$? Thanks in advance. * * * After thinking about this a little further, I came up with the following: Suppose we want the binary predictor to have expected prevalence $p$ and the multiple correlation with the other continuous covariates to be $0.75$. Let $$ (\tilde{X}_1,X_2,\ldots,X_n)\sim \mathcal{N}_n(\mathbf{0},\mathbf{S}), $$ where $\mathbf{S}$ is the correlation matrix given by $$ \mathbf{S}= \begin{pmatrix} 1 & 0.75 & 0.75 & \cdots & 0.75 \\\ 0.75 & 1 & 0.25 & \cdots & 0.25 \\\ \vdots & 0.25 & 1 & \cdots & \vdots \\\ 0.75 & 0.25 &\cdots & \cdots & 1 \end{pmatrix}. $$ If we let $a\in\mathbb{R}$ be given such that $P(\tilde{X}_1>a)=p$ and $X_1=1_{\\{\tilde{X}_1>a\\}}$, then $X_1\sim b(1,p)$. But has $X_1$ the same correlation properties as $\tilde{X}_1$? Furthermore, the $\mathbf{S}$ given above only says that $\tilde{X}_1$ has a pairwise correlation with $(X_2,\ldots,X_n)$ of $0.75$ and not a multiple correlation of $0.75$.


title: Comparing 2 classifiers with unlimited training data

text: I would like to compare 2 text classifiers C1 and C2, which can be trained with "unlimited" noisy training datasets, meaning that you can use as much data as you want for training, such data being very noisy. I am thinking of the 2 options below. **Which one would be better?** 1) Set up a datasize size interval 0-z, z being a very large number. Identify dataset size x at which C1 scores maximum accuracy (eg using a learning curve). Then train C2 with such dataset size x. 2) Identify dataset size y at which C1 stops learning (eg very small slope of the learning curve). Then train C2 with such dataset size y. Please not that I am not asking about method for comparing classifiers, but how to decide the training dataset size at which 2 classifiers should be compared. **Edit:** Below is the learning curve for one classifier. X axis is training dataset size, Y axis is accuracy. For the purpose of comparing two classifiers, should I pick the size at which highest accuracy has been found (red), or the size at which the accuracy remains (more or less) stable (green)? ![enter image description here](http://i.stack.imgur.com/wDykB.png) **Edit2** : @cbeleites @Dikran Marsupial I am sorry for the lack of information:   * I have 3 classes in my problem.   * The above learning curve was generated by increasing training dataset size, starting from 1,000 examples/class(3,000 in total) until 140,000 examples/class(420,000 in total); for each increase step 1,000 new examples/class(3,000 in total) are being added, and the trained model is tested for every iteration with the same test dataset, composed of 350 examples/class (1,050 examples in total)   * Training dataset samples were automatically labeled, test samples were manually labeled


title: Enrichment analysis by gene duplication level

text: **Biological Background** Over time, some plant species tend to duplicate their entire genomes, gaining an additional copy of each gene. Due to the instability of this setup, many of these genes are then deleted, and the genome rearranges itself and stabilises, ready to duplicate again. These duplication events are associated with speciation and invasion events, and the theory is that the duplication helps plants to adapt faster to their new environments. Lupinus, a genus of flowering plant, invaded the Andes in one of the most rapid speciation events ever detected, and what's more, it seems to have more duplicate copies in its genome than the most closely related genus, Baptisia. **And now the mathematical problem:** The genomes of a member of Lupinus and a member of Baptisia have been sequenced, providing raw data about 25,000 genes in each species. By querying against a database of genes of known function, I now have a "best guess" for what functions that gene might perform - so for example, Gene1298 might be associated with "fructose metabolism, salt stress response, cold stress response". I want to know, if there was a duplication event between Baptisia and Lupinus, whether gene loss took place at random, or whether genes performing particular functions were more likely to be kept or deleted. I have a script which will output a table like the one shown below. L * is a count of all Lupinus genes associated with the function. L 1+ is a count of lupinus genes associated with the function where at least one duplicate copy exists. I can get it to produce L 2+, L 3+ etc., although L 1+ is a much more reliable group than L 2+ due to the sequencing process.               Function            | L *  | L 1+ | B *  | B 1+ |     fructose metabolism | 1000 | 994  | 1290 | 876  |     salt stress         | 56   | 45   | 90   | 54   |     etc.      What I would like to do is to test, for each gene function, whether there are more or fewer genes with duplicates than might be expected purely by chance in Lupinus and Baptisia, and whether Lupinus differs from Baptisia in the ratio of observed to expected. **The best thing I have so far** Previous studies on different species have used Enrichment Analysis, with Fisher's Exact Test and FDR correction for multiple sampling, to do a contingency test on each row. It would be nice to improve on this; I'm not sure this sounds like the best way to do it. Glen_b has suggested using a GLM to analyse the data; I have played around with GLMs in JMP8, which has been interesting, but I will admit to not really understanding them. That said, I'm trying to use R instead now. **What am I using this for?** This was originally supposed to be as part of a short research project I'm doing at university, but has now spanned off into an enormous genome annotation project. Why? Because bioinformatics is cool. Being able to take a string of A,T,C and G and use it to infer information about events which happened millions of years ago is amazing. Needless to say, I am not going to try and submit any kindly provided answer as my own work. I would be happy to include an acknowledgement in the paper if I use a method suggested here in the submitted work.


title: Multinomial log. regression in SPSS with a dummy independent variable

text: This is my test design which is simply a grouping variable: ![enter image description here](http://i.stack.imgur.com/8bTwB.jpg) When I choose to regress a gender variable (w/m) onto this grouping variable with a multinomial logisitic regression model I'll observe ![enter image description here](http://i.stack.imgur.com/W2z0K.jpg) Why does SPSS list both categories of gender with one being redundant rather then choosing one by default (the last category) as the reference? If there's no reference for the gender variable the interpretation of the intercept gets rather confusing. Is there simply a mistake involved or do I miss something?


title: How to calculate goodness of fit in glm (R)

text: I have the following result from running glm function. How can I interpret the following values:   * Null deviance   * Residual deviance   * AIC Do they have something to do with the goodness of fit? Can I calculate some goodness of fit measure from these result such as R-square or any other measure?               Call:     glm(formula = tmpData$Y ~ tmpData$X1 + tmpData$X2+              tmpData$X3 + as.numeric(tmpData$X4) + tmpData$X5 +          tmpData$X6 + tmpData$X7)          Deviance Residuals:           Min        1Q    Median        3Q       Max       -0.52628  -0.24781  -0.02916   0.25581   0.48509            Coefficients:                                   Estimate Std. Error  t value Pr(>|t|)         (Intercept         -1.305e-01  1.391e-01   -0.938   0.3482         tmpData$X1         -9.999e-01  1.059e-03 -944.580   <2e-16 ***         tmpData$X2         -1.001e+00  1.104e-03 -906.787   <2e-16 ***     tmpData$X3         -5.500e-03  3.220e-03   -1.708   0.0877 .           tmpData$X4         -1.825e-05  2.716e-05   -0.672   0.5017         tmpData$X5          1.000e+00  5.904e-03  169.423   <2e-16 ***         tmpData$X6          1.002e+00  1.452e-03  690.211   <2e-16 ***     tmpData$X7          6.128e-04  3.035e-04    2.019   0.0436 *       ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for gaussian family taken to be 0.08496843)              Null deviance: 109217.71  on 3006  degrees of freedom     Residual deviance:    254.82  on 2999  degrees of freedom       (4970 observations deleted due to missingness)     AIC: 1129.8          Number of Fisher Scoring iterations: 2


title: PCA with "zero-spike" variables

text: I have a large number of predictors that are hypothesized to be important in determining my binary outcome variable (here's a bit more about my goal, the predictors and the outcome). The problem is that considering the number of "positive cases" that I have, I have too many predictors (I have about 450 cases in total, about 100 positive cases, and more than 20 predictors.) I'd like to use PCA to reduce my data, but many of my predictors are non- negative with a spike at zero and a long tail. For instance, the value for "the number of cigarettes per day smoked in the past 30 days" is zero for about 75% of the respondents and 99+ for a few. I understand that there's nothing conceptually wrong with using PCA with such variables. But PCA finds the directions of maximum variance in data and I'm not sure that variance is the right measure of dispersion for this kind of data. (These variables also cause all sorts of errors and warnings to be generated when I use other data reduction software.) I'd like to find a non-response-variable-driven way to deal with these variables. What I've tried already is to use $x \rightarrow log(x + 1)$ transformation followed by restricted cubic spline transforms (as implemented by $transcan$ function in Hmisc package in R), but they don't improve the situation much. I've seen this paper by Royston et al. that explicitly deals with this problem, but their approach is response-variable driven, which is something I'd like to avoid. I'd love to hear everyone's thoughts, experiences, and suggestions. Thanks!


title: Why does MICE fail to impute multilevel data with 2l.norm and 2l.pan?

text: Why does MICE fail to impute multilevel data with 2l.norm and 2l.pan in this situation ? Here is a reproducible example:               require(foreign)     require(mice)     require(pan)          dt.fail <- read.csv("http://goo.gl/pg8um")     dt.fail$X <- NULL          dt.fail$out <- as.factor(dt.fail$out )     dt.fail$grp<- as.factor(dt.fail$grp)     dt.fail$v1<- as.factor(dt.fail$v1)     dt.fail$v2<- as.factor(dt.fail$v2)     dt.fail$v3 <- as.factor(dt.fail$v3)     dt.fail$v7<- as.factor(dt.fail$v7)     dt.fail$v8 <- as.factor(dt.fail$v8)     dt.fail$v9 <- as.factor(dt.fail$v9)     dt.fail$v11 <- as.factor(dt.fail$v11)     dt.fail$v12 <- as.factor(dt.fail$v12)          dt.fail <- dt.fail[!is.na(dt.fail$grp),]          PredMatrix <- quickpred(dt.fail)     PredMatrix['CTP',] <- c(1,-2,0,0,0,0,0,0,0,0,1,0,1,1,0,2)               impute = mice(     data=dt.fail,      m = 1,      maxit = 1,         imputationMethod = c(             "logreg",   # out         "",     # grp   ----> cluster grouping factor             "pmm",  # v1             "polyreg",  # v2             "logreg",   # v3             "pmm",  # v4             "logreg",   # v5             "logreg",   # v6         "polyreg",  # v7              "polyreg",  # v8              "polyreg",  # v9              "polyreg",  # v10             "",     # v11 ----> complete             "",     # v12 ----> complete             "2l.pan",   # CTP ----> multilevel imputation             ""),        # const ----> needed for multilevel impuitation     predictorMatrix = PredMatrix, seed = 101     )      This produces the following error:               Error in order(dfr$group) : argument 1 is not a vector      Using the `2l.norm` method, it produces the following error:               Error in factor(x[, type == (-2)], labels = 1:n.class) :  invalid labels; length 20592 should be 1 or 2      Using `pmm` there is no error


title: Regression for poisson process in R

text: I have a series of samples of varying length, and the number of bugs created in those time samples. Reading the literature, this is often modeled as a Poisson process. If you write it like: $$P(k,t)=\frac{e^{-\lambda t}(\lambda t)^k}{k!}$$ I know $k$ and $t$ and would like to regress to find $\lambda$. Is there an R package to do this? The standard `glm` works for Poisson regression, but I can't find anything which allows for samples having non-unit times.


title: Correlation of lognormal processes

text: Proposition: If X(t) is a lognormal process then corr(X(t),1/X(t))=-1.   1. What is a lognormal process? What are the differences between a lognormal process and a lognormal distribution?   2. Can anybody prove the above proposition?


title: Analysing a multiple choice question

text: I asked my questionnaire's respondents whether they are aware of a particular subject and provided them with a list of four statements to choose from. I need to analyse whether their awareness differs among these respondents, say by comparing this question with the demographics (variables having an ordinal or nominal scale). I would like to know which statistical test I should use. Can I use the Chi squared test?


title: Revolutions in statistics for the last 50 years?

text: Let's say chat is a better place for such discussions. --mbq What areas of statistics have been substantially revolutionised in the last 50 years? For example, about 40 years ago, Akaike with colleagues revolutionised the area of statistical model discrimination. About 10 years ago, Hyndman with colleagues revolutionised the area of exponential smoothing. About XX years ago, ... How do I possibly continue the list, with years and names please? By statistics I mean its all four types from Bartholomew's 1995 presidential address, Chambers's greater and lesser statistics together, as featuring in Hand's recent presidential address on 'Modern statistics' and so on - anything professionally relevant.


title: Analysis of cross correlation between point-processes

text: I would like an advice on a analysis method I am using, to know if it it statistically sound. I have measured two point processes $T^1 = t^1_1, t^1_2, ..., t^1_n$ and $T^2 = t^2_1, t^2_2, ..., t^2_m$ and I want to determine if the events in $T^1$ are somehow correlated to the events in $T^2$. One of the methods that I have found in the literature is that of constructing a cross-correlation histogram: for each $t^1_n$ we find the delay to all the events of $T^2$ that fall in a given window of time (before and after $t^1_n$), and then we construct an histogram of all these delays. If the two processes are not correlated I would expect a flat histogram, as the probability of having an event in $T^2$ after (or before) an event in $T^1$ is equal at all delays. On the other hand if there is a peak in the histogram, this suggests that the two point process are somehow influencing each other (or, at least, have some common input). Now, this is nice and good, but how do I determine whether the histograms do have a peak (I have to say that for my particular set of data they're clearly flat, but still it would be nice to have a statistical way of confirming that)? So, here what I've done: I've repeated the process of generating the histogram for several (1000) times keeping $T^1$ as it is and using a "shuffled" version of $T^2$. To shuffle $T^2$ I calculate the intervals between all the events, shuffle them and sum them to reconstitute a new point process. In R I simply do this with:               times2.swp <- cumsum(sample(diff(times2)))      So, I end up with 1000 new histogram, that show me the density of events in $T^{2*}$ compared to $T^1$. For each bin of these histogram (they're all binned in the same way) I calculate the density of 95% of the histogram. In other words I'm saying, for instance: at time delay 5 ms, in 95% of the shuffled point processes there is a probability x of finding an event in $T^{2*}$ after an event in $T^1$. I would then take this 95% value for all of the time delays and use it as some "confidence limit" (probably this is not the correct term) so that anything that goes over this limit in the original histogram can be considered a "true peak". **Question 1** : is this method statistically correct? If not how would you tackle this problem? **Question 2** : another thing that I want to see is whether there is a "longer" type of correlation of my data. For instance there may be similar changes in the rate of events in the two point processes (note that they may have quite different rates), but I'm not sure how to do that. I thought of creating an "envelope" of each point process using some sort of smoothing kernel and then performing a cross-correlation analysis of the two envelopes. Could you suggest any other possible type of analysis? Thank you and sorry for this very long question.


title: Factor analysis on a "rotating" subset of survey questions

text: To help develop a survey instrument, I'd like to perform an exploratory factor analysis on a pool of 50-100 survey questions. (I'm open to discussing whether this is the best strategy for the project, but to keep the scope of this question reasonable, let's assume this is in fact a reasonable thing to do.) I expect to recover some four to eight factors. Because of logistics of administering the surveys, I'd rather administer a randomly-selected subset of 20-30 questions instead of the whole thing. Using standard statistical software, is it still reasonable to do EFA on the kind of data this would produce? Or do I need to administer the entire survey to every respondent for it to work?


title: What is the relation or difference between reliability and validity?

text: A test should be designed in a manner to get information to be believed in by investigator and academia world .


title: R: Computation of quantization error of Self-Organizing Maps

text: In the Kohonen package the quantization error is implemented as:               data(yeast)     ## take only complete cases     X <- yeast[[3]][apply(yeast[[3]], 1, function(x) sum(is.na(x))) == 0,]     yeast.som <- som(X, somgrid(5, 8, "hexagonal"))     ## quantization error:     mean(yeast.som$distances)      By definition die quantization error is the average distance to its closest node. An original implementation could look like this:               part.quant<- sapply(split(as.data.frame(yeast.som$data), yeast.som$unit.classif), ss)     mean.quant<- sum(part.quant) / nrow(X)      with               ss <- function(x) sum(scale(x, scale = FALSE)^2)      I thought I would get the same results, however they where slightly different. What did I miss out and how to correct it?


title: Is it possible to do a test of significance for a string occurrence in two datasets

text: I have two sets of data compiled from different sources. Both sets of data contain words with frequencies of occurrence. I would like to check if a certain word exists in both datasets and if it does, whether it is possible to perform some kind of test of significance statistically to prove the word is significant. ### For example:               word = 'apple'     dict1 = {'oranges': 45, 'apple': 34, ..., 'x': y}     dict2 = {'apple': 165, 'orange': 12, ..., 'x': y}      If the word "apple" appears in both datasets (dict1 and dict2), then calculate a test of significance for the word apple. ### Edit:   1. First, I want to check if both words exist in both datasets.   2. For example, if dict has 1000 words, and I arrange the words according to frequencies, I will get some kind of graph. If the top word has a frequency of 13,000 and "apple" has a frequency of 34, I would like to test whether the gap (13,000 - 34) is too huge and the word "apple" didn't appear enough times compared to the top word to be considered significant. However, if 80% of the words fall within the frequency of 20-50, then it's not a good idea to say "apple" is not significant.   3. I have 2 datasets where the word "apple" may appear. So, I need to make sure "apple" does not fall too low in the frequency value in either or both datasets.


title: Sample size for binary text classification

text: In my application as a Binary text classification, one class has around 36,000 sample and another one has around 300 samples.   * I under-sample the first class. So, each class will have ~ 300 samples.   * TfidfVectorizer with max_features = 5000 (in sklearn) is used for text vectorization(in sklearn).   * SGDClassifier (in sklearn) is used for classification (loss='hinge', penalty="elasticnet", shuffle=True, class_weight='auto')   * and here is the result      * score on training data: 1     * score on test data: 0.83     * confusion_matrix: * * *                 [[112  21]       [ 17  86]]                            precision    recall  f1-score   support              class-0          0.87      0.84      0.85       133         class-1          0.80      0.83      0.82       103             avg / total       0.84      0.84      0.84       236      Things that make me worry are that how can I get this good accuracy.Number of samples are too low compare to number of features (5000) and classifier should over-fitted but nothing is happening (as far as I detect). If I am doing something wrong? what is the minimum number of samples for training a text classifier?


title: Whether to use EFA or PCA to assess dimensionality of a set of Likert items

text: This follows on from my previous question on assessing reliability. I designed a questionnaire (six 5-points Likert items) to evaluate the attitude of a group of users toward a product. I would like to estimate the reliability of the questionnaire for example computing Cronbach's alpha or lambda6. So, I need to check the dimensionality of my scale. I have seen that some people use PCA to find out the number of dimensions (e.g. the principal components), other people prefer to use EFA.   * Which is the most suitable approach?   * Further, if I find more than one principal component or more than one latent factor does this mean that I am measuring more than one construct or several aspects of the same construct?


title: Selecting optimal number of input features and optimal number of hidden layers for a MLP?

text: What is the best way to select parameters for a binary neural network classifier? More specifically I have 265 features ranked according to Mutual Information Criterion. I have to determine the optimal number of inputs and the optimal number of hidden layer nodes for use in a multi-layer perceptron (MLP). Note: The selecting of the optimal number of hidden nodes and optimal number of input features are not independent.


title: Are studentized deleted residuals a form of k-fold cross validation when K=N?

text: Are Studentized deleted residuals a form of k-fold cross validation when K=N? (this question is asked in the context of the discussion here)


title: Testing the importance of an item among a finite set of items

text: Consider the following survey question: > >     Q: How would you classify the importance for you of the following 5 > items: >   >     A >     B >     C >     D >     E >   > > Assign to each item a number in the set {1,2,3,4,5}, with 1 meaning the > highest importance and 5 meaning the lowest importance; the number used to > an item cannot be used in any other item.   * How can one test whether item A is the most important item for a sample with 100 individuals?


title: Comparing two datasets (of the same physical quantity) - what do I learn from this graph?

text: I'm having a hard time understanding what the authors of this paper (pdf) want to tell me with this graph (Fig. 2,3 (shown below) and 4, right): ![relative error vs st. dev. of difference](http://i.stack.imgur.com/0Wti8.png) [Caption: A comparison between the standard deviation of the differences (green line) and the standard deviation of all GOMOS (red line) and LIDAR (blue line) ozone profiles.] _Note:_ contrary to convention, the measured quantity is plotted on the x-axis, not y-axis. The problem at hand is as following: Two instruments measure the same physical quantity (in this case vertical distribution of ozone in the atmosphere). Their corresponding (instrument characterized) standard deviations (averaged over if I understand correctly) are displayed in red and blue. The difference between the two datasets is computed and the standard deviation of the _difference data set_ is plotted in green. What does this graph mean? Is this a way to say that the two datasets are in agreement? If yes, on what basis? What would other extremes mean?   * stdev of the difference is much **larger** than the measurement stdevs   I think: measurement stdevs are too optimistic.   * stdev of the difference is much **smaller** than the measurement stdevs   I think: measurement stdevs are too conservative. Maybe an example would help me to understand this better. Thanks


title: Calculate mean of one variable while controlling for another using regression

text: This is probably a very naive question... I'd like to estimate "adjusted" or "conditional" means for a variable (i'm unsure of the correct terminology). My data are on cortisol levels (dependent variable) in rabbits (n=56). I have many measurements at different times of the day, over many months. I'd like to calculate mean weekly values of cortisol for each individual rabbit so these can be used as a predictor in another model for which I only have weekly data. Rather than calculate the means from the raw data, i'd like to control for the time of day the samples were taken (this can influence the measurement). I thought i'd regress time of day (in minutes from 00:00 each day) on cortisol level and then extract the fitted values and calculate the weekly mean for each rabbit from these. Would this give me the estimated mean for cortisol, while controlling for time of day? I can't share my data, but i've created a similar mock up using the iris data set. Here I fit a model, extract the fitted values and then calculate "adjusted" means for each species while controlling for the predictor. Am I right in thinking the difference between these means and the ones for the raw data (below) reflect the adjustment made when controlling for the independent variable?               data(iris)          fit <- lm(Sepal.Length ~ Petal.Length, data = iris)     summary(fit)          with(iris, plot(Sepal.Length ~ Petal.Length, col = as.numeric(Species), asp = 1))     abline(coef(fit))          iris$fitted <- fitted(fit)          with(iris, aggregate(fitted, list(Species), mean))     #      Group.1       x     # 1     setosa  4.9044     # 2 versicolor  6.0486     # 3  virginica  6.5769          with(iris, aggregate(Sepal.Length, list(Species), mean))     #      Group.1      x     # 1     setosa  5.006     # 2 versicolor  5.936     # 3  virginica  6.588


title: What would be the most appropriate significance test for my scenario?

text: I have three groups of users $-$ $G1$, $G2$ and $Control$. The users in each of these three groups are different but are carefully selected and have similar properties. I treat each group with a unique treatment. $G1$ is treated by $T1$, $G2$ is treated by $T2$ and the $Control$ group is not treated at all. For each of these groups, I measure a variable $X$ before the treatment and after the treatment. I have the following questions - [1] I would like to measure if my treatment (both $T1$ and $T2$) made a difference. In this case, I need to compare the variable $X$ for users in groups $G1 + G2$ to $Control$. Since I do not have any knowledge of the underlying distribution of data, I plan to use a Wilcoxon Signed Rank Test. Is this correct? If not, which other test would be appropriate? [2] I want to measure which of the two treatments ($T1$ or $T2$) performed better. How do I do this? Should I perform a statistical test between users in $G1$ and $G2$ OR perform a statistical test between $G1$ and $Control$ AND $G2$ and $Control$ and interpret from the corresponding result?


title: Comparing means across sub-groups?

text: Bit of a stats newbie question here. I've got about 10,000 rows of data about pupil performance, annotated with school and teacher details, that look as follows:               Pupil      School     Grade(1-5)    Teacher     123        Morley     4             Adams     124        Westlake   3             Morris     125        Morley     5             Adams     126        Westlake   2             Philips     127        St Xavier  4             Smith      I've calculated that the overall mean grade is 3.4 with standard deviation 1.0. Now, say I want to work out whether one school (or one teacher) has overall a significantly lower or higher performance than average. If I calculate that the overall mean grade at Morley is 3.9, with s.d. 0.8, the overall mean grade at Westlake is 3.2, with s.d. 0.7, and the overall mean at St Xavier is 3.8, with s.d. 0.5, what test can I use to say whether any of these schools have a significantly different performance from the average? From my knowledge so far, I'm guessing I want to use some form of t-test, probably Student's t-test - although I'm not sure how to adapt it for multiple groups. (Of course, there are lots of factors influencing why grades might be lower at one school than another and yet that school might still be performing well - I'm not getting that deep here. I just want to know how to compare means across sub-groups.) Many thanks for your help.


title: Regression analysis with 1 dependent and 1 independent variable

text: Is it possible to have a regression analysis with 1 dependent and only 1 independent variable?


title: In SVM, what are the labels and how do you get them from the data?

text: I'm working on a school project and have decided to use SVM for stock market prediction. I have a 1000x5 matrix of stock quotes containg data for open, close, high, low, volume data. From what I know, the labels should consist of +1/-1, for either up or down. How do I get this label data? I'm using python and I am somewhat lost. I'm ultimately trying to predict the stock price a few days out in the future, or at least one day out.


title: Probability the sum of numbered balls drawn from a box is odd

text: A box contains 100 balls, numbered from 1 to 100. If 3 balls are selected at random and with replacement from the box, what is the probability that the sum of the 3 numbers on the balls selected from the box will be odd? The answer to this question is 1/2. This is a a high level GMAT sample question and the working out in the answer shown is pretty long. When I picked 1/2 so quickly, the system assumed I guessed. My question is was it a lucky guess? Or can we say   * if we were to pick one ball instead of 3 the answer would have been 1/2?   * same answer if we were to pick 2 balls? If so, can we apply the same principle to 3 balls? Or is the rationale incorrect? What about if we were to pick 4 or 5 balls? Thank you.


title: What is a tertile average of a group?

text: I came across a paper , where they have taken tertile average of a hormone level in a group. I am wondering what difference does calculating tertile average make? Thank You!


title: Combining several survival experiments

text: In four independent experiments, two strains (A and B) were compared in their respective survival curves. Two experiments show a clear-cut result (A survives at 100%, almost half B dies). Two further experiments show no significance (log-rank / Mantel-Cox test and Breslow/ Gehan-Wilcoxon), but even in these insignificant experiments strain A always had 100% survival and all individuals that die belong to strain B. Unfortunately, I have no experience with survival tests. How could one combine the results of the four experiments other than just pool the data? Pooling the data seems to me like a wrong choice, because there clearly is an effect of the experiment. To make the matter more complex, there is a parameter D (dosage) that varied slightly between the four experiments. It _looks_ (without testing) that there is an effect of dosage on survival. How should one test it? (the current solution proposed by the experimentalists is to repeat the experiments with more individuals) EDIT: What I did until now was to combine the p-values from the experiments using Stouffer's method (Fisher's method gives an even lower p-value). The combined p-value in Stouffer's test is $ < 10^{-4}$.


title: How to fit a TAR-GARCH model in R

text: I'm trying to understand whether it is possible to fit a threshold AR model with a GARCH effect using R. So far, I only found the functions `garchFit` (which also allows to "merge" ARIMA and GARCH models) and `tar`. In my book (by Tsay) I found the code to fit such models on RATS, but I'm asked to use R. If that is impossible, I'd ask to suggest a solution to bypass this problem. I'm studying the continuously compounded log-returns of S&P500 index, and I successfully (at least from the tests!) fitted a GARCH model already. This because the series showed heteroskedasticity in the variance. But the plot also shows some leverage effect, so I thought that a TAR (better than a TGARCH) could explain this phenomenon.   Is it reasonable to try merging those two models? Should I leave the GARCH? Thanks.


title: Interval censoring

text: I ran an interval censor survival curve with R, JMP and SAS. They both gave me identical graphs, but the tables differed a bit. This is the table JMP gave me.               Start Time  End Time    Survival    Failure SurvStdErr     .            14.0000      1.0000     0.0000     0.0000     16.0000      21.0000      0.5000     0.5000     0.2485     28.0000      36.0000      0.5000     0.5000     0.2188     40.0000      59.0000      0.2000     0.8000     0.2828     59.0000      91.0000      0.2000     0.8000     0.1340     94.0000     .             0.0000     1.0000     0.0000      This is the table SAS gave me:               Obs Lower Upper Probability Cum Probability Survival Prob Std.Error     1    14    16      0.5          0.5             0.5        0.1581     2    21    28      0.0          0.5             0.5        0.1581     3    36    40      0.3          0.8             0.2        0.1265     4    91    94      0.2          1.0             0.0        0.0      R had a smaller output. The graph was identical, and the output was:               Interval (14,16] -> probability 0.5     Interval (36,40] -> probability 0.3     Interval (91,94] -> probability 0.2      My problems are:   1. I don't understand the differences   2. I don't know how to interpret the results...   3. I don't understand the logic behind the method. If you could assist me, especially with the interpretation, it would be a great help. I need to summarize the results in a couple of lines and not sure how to read the tables. I should add that the sample had 10 observations only, unfortunately, of intervals in which events happened. I didn't want to use the midpoint imputation method which is biased. But I have two intervals of (2,16], and the first person not to survive is failed at 14 in the analysis, so I don't know how it does what it does. Graph: ![enter image description here](http://i.stack.imgur.com/g59l3.png)


title: Precision and recall are equal when the size is same

text: I have a set of users. A classification algorithm is applied on all users, and I take (call `analyzedExperts`) a set of users which are binary classified (expert & non-expert). And I use another method to evaluate this algorithm. That does the same thing, and I take another set of users which are also binary classified (call `realExperts`). But if I want to measure the precision and recall, I take the same result for both. The sets `analyzedExperts` and `realExperts` have both the same size of data. I don't understand why they are same, and don't know whether it is normal. P.S. I'm not sure whether the precision and recall is a good way to measure the evaluating the results. **EDIT:** Thus, the question is: if they are equalsized, precision and recall have to be same? Suppose they have 3 users in common (True positive). FN and FP will be always same because they both have the same size. What implicits that the precision and recall will be same. Second question might be then: does `realExperts` has to have greater size? Or is it not the good place to use _precision_ or _recall_? ![enter image description here](http://i.stack.imgur.com/1bwxX.png)


title: Maximum likelihood estimation and density estimation

text: Let's consider a general signal processing estimation problem where the measurements are modeled as $${\bf x}[n]={\bf s(\theta)}+{\bf w}[n],$$ where ${\bf w}$ is a non-Gaussian r.v. (noise term) and ${\bf s(\theta)}$ is a deterministic nonlinear signal model. Since the noise is non-Gaussian I am wondering how one could design a Maximum Likelihood estimator for the parameter vector ${\bf \theta}$ if an estimate of the pdf of ${\bf w}$ was obtained by density estimation or another nonparametric technique. I have looked into empirical likelihood but I am not sure if it applies to this problem. Is this a standard problem in the context on non-parametric statistics? References are welcome.


title: Evaluating matrix factorization algorithms for Netflix

text: I've been trying to implement Simon Funk's movie recommendation algorithm explained here. I understand how the user and item factors are computed. However the evaluation method is not clearly explained. I checked this post also but there was no evaluation method there. My Question: Let's assume we divide the the user-rating matrix R into training(95%) and testing data-sets(5%): R_train and R_test and we use Simon Funk's algorithm to find the user factor matrix (U) and movie factor matrix (M) using the training data matrix: R_train=U*M. How do we evaluate algorithm using R_test?


title: Model selection of multinomial logistic regression for 2x3 contingency table

text: I want to analyze a 2x3 contingency table using multinomial logistic regression and I hope to be able to do it in Matlab or in R. I have looked around in old threads, but haven't been able to find a solution. My predictor can only take two values (male,female) and the response is death by one of three different diseases (stroke,cancer,other). Example data could look like this:                          Stroke    Cancer    Other     Male:        22         2       19     Female:      30        20       35      Using the matlab functions mnrfit and mnrval, I can fit a nomimal multinomial logistic model to the data and get the probability of dying of either of three diseases, for males and females - 6 values all in all: $p_{stroke}^{male},p_{cancer}^{male}$, $p_{other}^{male}$ and $p_{stroke}^{female}$, $p_{cancer}^{female}$, $p_{other}^{female}$. I can also calculate the upper and lower boundary of the 95% confidence interval (CI) of the probabilities from my fit. Since the 95% CIs of $p_{cancer}^{male}$ and $p_{cancer}^{female}$ do not overlap, I conclude that females are more significantly more likely to die from cancer than males are. Is this a correct interpretation of the model? Futhermore: 51% of males die from stroke, but only 35% of females die from stroke, so I thought that there would perhaps be a significant difference here also. But since the 95% CIs of $p_{stroke}^{male}$ and $p_{stroke}^{female}$ overlap, do I have to conclude that there is no gender difference in the probability of dying from stroke? I had the following idea for an improved analysis: I would like to fit several multinomial logistic models, and select the best one using BIC or some similar meausure. For example, I would like to fit these three models (and also other permutations): Model 1: $p_{stroke}^{male} = p_{stroke}^{female}$, $p_{cancer}^{male} \neq p_{cancer}^{female}$ and $p_{other}^{male} = p_{other}^{female}$ Model 2: $p_{stroke}^{male} \neq p_{stroke}^{female}$, $p_{cancer}^{male} \neq p_{cancer}^{female}$ and $p_{other}^{male} = p_{other}^{female}$ Model 3: $p_{stroke}^{male} \neq p_{stroke}^{female}$, $p_{cancer}^{male} \neq p_{cancer}^{female}$ and $p_{other}^{male} \neq p_{other}^{female}$ If "Model 2" is the best model, then I would conclude, that there is a gender difference both in the probability of dying from cancer and the probability of dying from stroke. Is this a sound way to analyze the data? And can someone suggest some literature or a worked example how to program this in matlab or R? Is there a package where you can fit multinomial models under the constraint that some probabilites must be the same? I'm new to this kind of modeling, but I wish to learn. Thanks so much in advance!


title: A Poisson EM clustering algorithm in R?

text: I've been looking for a Poisson expectation-maximization clustering implemented in R. Does anyone know if there is one available or can point me in the right direction? thanks! and please let me know if there's a better way/place to ask questions.


title: Effects of different changes on test statistic

text: Will each of the following increase, decrease, or have no effect on the value of the test statistic in a one-independent sample t-test?   1. The sample size is doubled    2. The level of significance is reduced   3. The sample variance increases    4. The difference between the sample and population mean increased This was given as part of an assignment and I can't find a set of rules that helps with this question.