title: Problem with informative censoring

text: I am reading "Monte Carlo Statistical Methods" by Robert and Cassella, and problem 1.3 asks > In example 1.1, the distribution of the random variable $Z=\min(X,Y)$ was of > interest. Derive the distribution of $Z$ in the following case of > _informative censoring,_ where $Y\sim N(\theta,\sigma^2)$ and $X\sim > N(\theta,\theta^2\sigma^2)$. Pay attention to the identifiability issues. Now I am almost certainly missing something, but my understanding is that informative censoring happens when $X$ and $Y$ are not independent. However, just knowing that they are not independent is not enough information to get the joint distribution, but if they are independent, I do not see any identifiability issues. _Added:_ If $X$ and $Y$ are independent it is straightforward but tedious to write down the distribution of $Z$, the tedium exacerbated by the fact that the distribution of $X$ is a delta function when $\theta$ is $0$. For a given distribution however, we can find $\theta$ as the (obviously unique) third quartile of the distribution, and given $\theta$, $\sigma^2$ is just a scaling parameter, so there are no identifiability issues that I can see. So, in summary, my questions are:   * What precisely is the definition of _informative censoring_ and why is the censoring in this exercise informative?   * If we are meant to take $X$ and $Y$ as independent, what are the identifiability issues that need to be paid attention to?  _Further addition_ With Ocram's explanation of _informative censoring_ it is now clear that the identifiability issues that neede to be paid attention to were that there were not any. If the parameters for the failure and censoring distributions were separate then there would be identifiability issues as we could swap the two distributions and get the same result. If someone more knowledgeable than I is feeling particularly Quixotic, please consider clarifying the Censoring wikipedia page.


title: Is using a questionnaire score (EuroQol's EQ-5D) with a bimodal distribution as outcome in linear regression a problem?

text: There is currently a debate whether the EQ-5D score that has a **ceiling** problem and a **bimodal** distribution can be used in a linear regression model or not. ### Background The score is **very** simple and frequently used to assess patient's _health related quality of life_ and consists out of five questions where each has 3 possible answers (there is a newer with 5 answers but it's less used). The score is commonly used in national registries as a patient reported outcome measure (PROM) and is very convenient because the questions are easy to answer and the completeness is therefore good. ### Continuous score The score is created by using a "tariff" where unique combinations of the 5 variables translate into a continuous-like variable but with the above mentioned limitations. I'm not sure how they decide on the tariff score but the different combinations of the answers combine into a unique value, for instance if you have answered _best health_ on all five categories you get a code of 11111 that gives the maximum of 1.000. If you've answered best on the first 2 question and worst on the last 3 you have a code of 11333 and get a score of -0.066. The score is country-adjusted and ranges between -0.594 to 1.000 in my Swedish tariff. ### The Paretian calculation In most orthopaedic studies we have a preoperative score and a postoperative score. By comparing the two models as the _Paretian Classification of Health Change_ suggests we get four possible outcomes; no change, worse, improved, or mixed change. Mixed meaning that one category became improved while another one deteriorated. As I understand the Paretian outcome is best analyzed using a multinomial logistic model. ## My questions   * When having large datasets of > 10 000 patients does it matter that the score is not normally distributed and is the Paretian way of analyzing the score better?    * Scores like this are very frequently used today - what are the limitations? # Update After taking all these wise arguments and discussing them closer with our statistician I got some interesting input:   1. In large sample the central limit theorem will kick in as long as the sample isn't heavily skewed   2. If the score itself has a flaw (as the EQ-5D score) it might not be right to expect a normal distribution because the bimodality is not due to a subgroup but due to a score feature (I think this is a different way of putting what @whuber wrote: "... The residuals will closely reflect that error distribution"   3. The normality of the sample helps in calculating the p-value/confidence interval and this could be circumvented by using bootstrapping   4. Using ordinal regression and leaving out the mixed group we can validate the results from the linear regression - i.e. show that the predictors behave similarly when used "non-parametric"


title: Log-likelihood of multivariate Poisson distribution

text: I'm having difficulty getting the gradient of the log-likelihood of a multivariate Poisson distribution. Here's how I have it setup:   1. A finite set of finite-dimensional vectors $T$ with elements $\mathbf{t}$   2. $d$ functions $\left\\{f_1,f_2,\dotsc,f_d\right\\}$ with compact support.   3. The parameter of the multivariate Poisson is given by $\lambda_{\mathbf{t}}\left(\boldsymbol{\theta}\right) = \sum_{k=1}^{d}\theta_k f_k\left(\mathbf{t}\right)$.   4. A sample from this distribution looks like this: $y_\mathbf{t}\sim\textrm{ Poisson}\left(\exp\left(\lambda_{\mathbf{t}}\left(\boldsymbol{\theta}\right)\right)\right)$   5. Multivariate Poisson likelihood function: $$L\left(\boldsymbol\theta\right)=\prod_{\mathbf{t}\in T}\frac{\exp\left(-\lambda_{\mathbf{t}}\left(\boldsymbol{\theta}\right)\right)\left(\lambda_{\mathbf{t}}\left(\boldsymbol{\theta}\right)\right)^{y_\mathbf{t}}}{y_\mathbf{t}!}$$ Here's where I am: \begin{align*} l\left(\boldsymbol\theta\right)&=\sum_{\mathbf{t}\in T}\log\frac{\exp\left(-\lambda_{\mathbf{t}}\left(\boldsymbol{\theta}\right)\right)\left(\lambda_{\mathbf{t}}\left(\boldsymbol\theta\right)\right)^{y_{\mathbf{t}}}}{y_{\mathbf{t}}!}\\\ &\ldots\textrm{ a little bit of algebra later }\\\\\ &=\sum_{\mathbf{t}\in T}\left(-\lambda_\mathbf{t}\left(\boldsymbol\theta\right) + y_\mathbf{t}\log\left(\lambda_\mathbf{t}\left(\boldsymbol\theta\right)\right)\right)-\log\left(y_\mathbf{t}!\right) \end{align*} What is the next step to take in terms of the derivatives? I'm not sure how to take derivatives with respect to $\boldsymbol\theta$ (i.e., what is the resulting type from $\frac{\mathrm{d}}{\mathrm{d}\,\boldsymbol\theta}\left(-\lambda_\mathbf{t}\left(\boldsymbol\theta\right)\right)$; is it a matrix, a vector, etc.). I would appreciate it if people's answers gave as little away about the problem as possible, I'd like to be able to finish deriving the equation myself; I just need a little push in the right direction. Much appreciated!


title: Correlation under uncertainty

text: I've got two sets of measurements whose uncertainty has been quantified: $x_1..x_n$ and $y_1..y_n$. For each $x_i$ and $y_i$, I have not only a point estimate but also an interval that expresses my uncertainty in the accuracy of the measurement. Now, I'd like to compute a correlation between $x$ and $y$. Is there a reasonable way of taking the intervals into account?


title: Instrumental variables with dichotomous instrument, variable, and outcome? Can I use TSLS?

text: Say I have a binary treatment, X, and I'm interested in its (local average treatment) effect on a binary outcome variable, Y. I find (binary) instrument T that directly affects X but not Y. How do I estimate this? I get that linear models like TSLS might not be ideal, but Angrist and Pischke's "Mostly Harmless Econometrics" seems to suggest that the effect estimates from something like a bivariate probit will be similar and maybe it isn't worth the trouble. So I'm wondering: do I need to run something other than TSLS here? And if so, is there an R package for that? Thank you very much.


title: Estimating changes of choice probabilities with choice removal in nested logit

text: Suppose that I have a nested logit model. I can estimate the probability of making a particular choice for each decision maker. Now, I want to remove one or more of the choices (though each nest keeps at least one option). How can I estimate/simulate the new choice probabilities? As a simple example, I might estimate a logit model (non-nested in this case, but you get the idea) over the choice of fruit. Consider an individual that actually chose a banana over an apple and an orange. I can estimate the probability that he would choose each kind of fruit. Suppose that a banana is no longer a choice. Now, what is the chance that he chooses the apple? And the orange? Thanks.


title: How is the minimum of a set of random variables distributed?

text: If $X_1, ..., X_n$ are independent identically-distributed random variables, what can be said about the distribution of $\min(X_1, ..., X_n)$ in general?


title: Should PCA or correlations be examined first in the context of correlated predictors in regression?

text: I'm looking for clarification: Do we have to run correlation analysis before carrying out principal component analysis (PCA), or is this implicitly subsumed under PCA framework? I am on my way to write a research plan:   * Isolate predictors with high collinear relationship to the process.   * Run PCA to convert correlated predictors to uncorrelated principal component variables. I am wondering if correlation analysis has to be done separately on the dataset in order to choose highly correlated predictors before applying PCA.


title: What is the difference between tapped delay line and sliding window in a neural network?

text: I'm trying to model some time series data, and I've been reading about tapped delay line and sliding window to transform the input data. In my understanding a sliding window with windows size 1 shifts the inputs by one at each time step and feeds them back to a feed forward neural network.               target| input          2| 1 0          3| 2 1          4| 3 2          5| 4 3          6| 5 4          7| 6 5      In my understanding, a tapped delay line with delay 2 takes the two most recent inputs and outputs and feeds them back to a recurrent neural network               y(t) = f(x(t-1),y(t-1))      Is there a difference? Why are there different names? Please correct me if I'm wrong.


title: Correct formula for standardized Student's t-distribution

text: I am wondering about the correct formula of the standardized Student's t distribution. In the rugarch package on page 15 it is given as: > Substituting $\frac{(\nu-2)}{\nu}$ into 49 we obtain the standardized > Student's distribution: $$ > f\left(\frac{x-\mu}{\sigma}\right)=\frac{1}{\sigma}f(z)=\frac{1}{\sigma}\frac{\Gamma\left(\frac{\nu+1}{2}\right)}{\sqrt{(\nu-2)\pi}\Gamma\left(\frac{\nu}{2}\right)}\left(1+\frac{z^2}{(\nu-2)}\right)^{-\left(\frac{\nu+1}{2}\right)} > $$ whereas from a book (Carol Alexander Market Risk Analysis, Quantiative Methods in Finance, Vol. 1 page 99) I have the following formula: > Setting $\mu=0$ and $\beta\nu=\nu-2$ in (I.3.55) given the density for the > _standardized Student t distribution_ , i.i. the Student _t_ distribution > with zero expercation and unit variance. This that the density function $$ > f_\nu(x)=((\nu-2)\pi)^{-1/2}\Gamma\left(\frac{\nu}{2}\right)^{-1}\Gamma\left(\frac{\nu+1}{2}\right)\left(1+\frac{x^2}{\nu-2}\right)^{-\left(\frac{\nu+1}{2}\right)} > $$ Now I don't know which of them is the correct one? Why is there an additional $\frac{1}{\sigma}$ in the first formula?


title: OLS regression - robust estimates for parameter's variance

text: I'm estimating a model for corporate social responsibility (not important). I have found my variable of interest significant at 5% confidence level. My sample is $N=84$, cross-section. For this I used OLS estimation. Then I doubted if I should have estimated the model with robust estimators (that is, HAC, HC1, ...), and did it. It turned out that the parameter is now only weakly significant (i.e. not at 5%, but at 10%). My questions are:   1. with such a small sample $N=84$, is it appropriate to use robust estimators?    2. can I test for heteroskedasticity some other way, so I can rely on my non-robust estimation?


title: Calculating influence for PLS with outer loadings and beta coefficient of latent variable

text: I'm calculating a penalized least squares regression (PLS). Two influence variables are connected to a latent variable. This latent variable has an influence (beta coefficient) of 0.6 on the response (or latent variable of the response). Both influence variables are correlated (outer loading) with 0.6 and 0.8 to the latent variable. Does anyone know how to (back)-calculate a beta coefficient for each influence variable and not just the latent variable connected to the influence variables? I'm using a mix of the smartPLS programming tool and R. EDIT: http://stackoverflow.com/questions/24446300/sempls-package- backcalculating-influence


title: Directional research hypotheses vs. directional hypothesis testing

text: When I took a methods course, I was told to avoid making a directional hypothesis. I understand the reason. Now however, I often see hypotheses like these: "We hypothesized that reaction time will decrease due to the treatment". It's rare to read research papers that have non-directional hypotheses (at least in my limited readings). If the researcher makes a directional hypothesis, I would assume that he will need to run a one-tailed test. From what I remember, if the results go the opposite direction--i.e., reaction time increases--our researcher failed to reject the null. What if an ANOVA is used instead? Does that negate the directionality of the hypothesis? From my understanding, the ANOVA is a non-directional test. On a related note, would rewording the hypothesis to "we hypothesized that reaction time **may** decrease due to the treatment" make it non-directional? Do researchers "change" their hypothesis after the fact to make it non- directional because they found an opposite effect? Sounds rather unethical to me.


title: Maximum likelihood estimation for mixed Poisson and Gaussian data

text: **Background** I've been doing a little bit of work lately on maximum likelihood estimation (MLE), for cases where the data is normally-distributed and also for cases where the data is Poisson distributed. For these two cases, the likelihood $L$ is given by: $$ L\left ( \mathbf{a} \right ) = \prod_{i} \mathrm{P}\left ( c_{i};m_{\mathbf{a}}\left ( x_{i} \right ) \right ) $$ Following the notation and steps of the first reference below, $\mathrm{P}\left ( c_{i};m_{\mathbf{a}}\left ( x_{i} \right ) \right )$ is the probability that a measurement gives $c_{i}$ if the true value is given by the model $m_{\mathbf{a}}\left ( x_{i} \right )$, where $\mathbf{a}$ is the set of parameters for the model. For the two different distributions, the reference gives: $$ L_{G}\left ( \mathbf{a} \right ) = \prod_{i} \frac{1}{\sqrt{2\pi \sigma_{i} ^{2}}} \mathrm{e}^{-\frac{\left ( c_{i}-m_{\mathbf{a}}\left ( x_{i} \right ) \right )^{2}}{2\sigma_{i} ^{2}}} $$ for Gaussian distributed data, and $$ L_{P}\left ( \mathbf{a} \right ) = \prod_{i} \frac{\left [m_{\mathbf{a}}\left ( x_{i} \right )\right ]^{c_{i}}}{c_{i}!}\mathrm{e}^{-m_{\mathbf{a}}\left ( x_{i} \right )} $$ for Poisson-distributed data. I'm happy following the subsequent steps to get a likelihood ratio test out. The steps involve calculating the negative-log of the likelihood, then finding a ratio based on maximising $L$ for a given model with respect to the global maximum of the likelihood $L$, i.e. $$ \frac{\mathrm{max}_{\mathbf{a}}\;L(\mathbf{c}|m_{\mathbf{a}})}{\mathrm{max}\;L(\mathbf{c}|\mathbf{m})} $$ **Edit** This question is related: What is distribution of $Z = X + Y$ where $X$ is Poisson distributed and $Y$ is normally distributed? **Question** I'm comfortable with the MLE for Gaussian and Poisson data. The scenario I'm interested in is for mixed Poisson-Gaussian data. An example might be photons arriving at a detector, which will exhibit Poisson noise, and the signal from the detector then subsequently being corrupted by Gaussian noise (e.g. thermal noise in the electronics). The noise model is then a mix of Poisson and Gaussian noise. What I want to do is test whether a Poisson, Gaussian or mixed model is most appropriate for various parameters (I don't necessarily know whether the Poisson data can be approximated by a Gaussian by the way, hence the question - it's what I'm trying to test!). What would an appropriate likelihood function $L$ be? Is it just a sum of the Poisson and Gaussian likelihoods? If so, the log-likelihood will be a bit tricky to simplify, I believe? I'm sure I'm missing something with this problem, I just haven't been able to make the leap from the two cases described above (I should add that it's been a bit of time since I did statistics in any serious depth). Any help much appreciated! **References** For this problem, I've been following the paper _"Comparison of maximum likelihood estimation and chi-square statistics applied to counting experiments"_ by T. Hauschild and M. Jentschel (2001). Link and DOI: 10.1016/S0168-9002(00)00756-7


title: STRING KERNELS FOR LIBSVM

text: I'm working on a protein classification problem and i'm using edit distance kernel defined in libSVM. Now, for instance, the implementation of spectrum kernel is very difficult, but i want try to test another kernel and preferably the combination between edit distance and another kernel too. Are there some c++ implementations of string kernels? I know that there is Shogun toolbox and maybe using precomputed kernel i can do it, but, if it is possible, i want copy a string kernel function into svm.cpp directly!


title: In Bayesian hypothesis testing, do the prior model probabilities have to be equal?

text: The posterior odds is the product of the Bayes factor and prior odds: $\frac{p(M_1|data)}{p(M_2|data)}=\frac{p(data|M_1)}{p(data|M_2)}\times\frac{p(M_1)}{p(M_2)}$. I was under the impression that the prior model odds could be specified as anything by the analyst, but then I read a passage in a psychology paper, Myung and Pitt (1997, p. 84) which states, > p(M1) represents the probability that the model M1 is a true description of > the events under study before data are collected. For real-world problems, > this probability may not exist, or not be known even if it does. We assume > equal model priors; that is, p(M1)/p(M2). This assumption is often made by > researchers in Bayesian model selection...and is an unavoidable > simplification of a problem for which there is not as yet a satisfactory > solution. An implication of the assumption is that model selection is based > solely on the Bayes factor. The decision to ignore model priors does not > undermine use of the Bayesian approach. It appears that they are stating that specifying equal prior probabilities is necessary (i.e., "unavoidable"), but this seems contrary to other things I have read and the structure of the equation above. So must the prior odds be set to 1? Myung, I. J., & Pitt, M. A. (1997). Applying Occam’s razor in modeling cognition: A Bayesian approach. _Psychonomic Bulletin & Review_, _4_ , p. 79-95.


title: Dirichlet sample by normalising Gamma RVs

text: I know that if you sample $K$ random variables $(X_1, X_2, \dots, X_K)$ from Gamma distributions using shape parameters $(\alpha_1, \alpha_2, \dots \alpha_K)$ and a scale parameter $\theta = 1$ such that $X_i \sim \Gamma(\alpha_i,\theta) = \Gamma(\alpha_i,1)$ then $\left(\frac{X_1}{\sum_{i = 1}^K X_i}, \frac{X_2}{\sum_{i = 1}^K X_i}, \dots, \frac{X_K}{\sum_{i = 1}^K X_i}\right) \sim \textrm{Dir}(\alpha_1, \alpha_2, \dots, \alpha_K)$ where $\textrm{Dir}(\alpha_1, \alpha_2, \dots, \alpha_K)$ is a Dirichlet distribution with a concentration parameter $(\alpha_1, \alpha_2, \dots, \alpha_K)$. My question is does this result hold for any scale parameter $\theta > 0$? The proofs I've read (e.g. http://mayagupta.org/publications/FrigyikKapilaGuptaIntroToDirichlet.pdf), which use the change-of-variables formula, seem to set $\theta = 1$.


title: (R) IF statement: stop & warn if FALSE else continue

text: this feels like it should be easy but I can't think my way around it - grateful for any tips. I'm making a function and before it does any of the hard stuff I need it to check that all the column names listed in the 'samples' dataset are also present in the 'grids' dataset (the function maps one onto the other). all(names(samples[expvar]) %in% names(grids)) This does that: the code within all() asks if all the names in the list ('expvar') of columns in 'samples' are also names in 'grids'. The output for a correct length=3 expvar would be TRUE TRUE TRUE. 'all' asks if all are TRUE, so the output here is TRUE. I want to make an IF statement along the lines of: if(all(names(samples[expvar]) %in% names(grids)) = FALSE) {stop("Not all expvar column names found as column names in grids")} No else needed, it'll just carry on. The problem is that the '= FALSE' is redundant because all() is a logically evaluable statement... is there a "carry on" function, e.g. if(all(etc)) CARRYON else {stop("warning")} ? Or can anyone think of a way I can restructure this to make it work? Thanks in advance!


title: The coefficient of one variable in Cox regression becames negative when doing multiple variable regression

text: I am doing Cox proportional regression. I first did Cox regression with only one variable `var`; its coefficient 0.8752721 was about what I expected and is easy to interpret in my situation. But when I added some other variables, the coefficient of `var` in the model became negative (-0.894304). I have performed the stepwise process. The negative coefficient seems wrong and is hard to interpret in my problem. Where am I wrong? Thanks.


title: Numerical estimation of MLE in Python -- normal distribution and gradient is close to zero away from the mean

text: I am exploring how to model a data set using normal distributions with both mean and variance defined as linear functions of independent variables. Something like $\mathcal{N} \sim \left (f(x), g(x)\right)$. I generate a random sample like this:               def draw(x):         return norm(5 * x + 2, 3 *x + 4).rvs(1)[0]      So I want to retrieve 5, 2 and 4 as the parameters for my distribution. I generate my sample:               smp = np.zeros((100,2))          for i in range(0, len(smp)):         smp[i][0] = i         smp[i][1] = draw(i)      The likelihood function is:               def lh(p):         p_loc_b0 = p[0]         p_loc_b1 = p[1]         p_scl_b0 = p[2]         p_scl_b1 = p[3]              l = 1         for i in range(0, len(smp)):             x = smp[i][0]             y = smp[i][1]             l = l * norm(p_loc_b0 + p_loc_b1 * x, p_scl_b0 + p_scl_b1 * x).pdf(y)              return -l      So the parameters for the linear functions used in the model are given in the p 4-variable vector. Using scipy.optimize, I can solve for the MLE parameters using an extremely low `xtol`, and already giving the solution as the starting point:               fmin(lh, x0=[2,5,3,4], xtol=1e-35)      Which does not work to well:               Warning: Maximum number of function evaluations has been exceeded.     array([ 3.27491346,  4.69237042,  5.70317719,  3.30395462])      Raising the `xtol` to higher values does no good. So i try using a starting solution far from the real solution:               >>> fmin(lh, x0=[1,1,1,1], xtol=1e-8)     Optimization terminated successfully.              Current function value: -0.000000              Iterations: 24              Function evaluations: 143     array([ 1.,  1.,  1.,  1.])      Which makes me think: PDF are largely clustered around the mean, and have very low gradients only a few standard deviations away from the mean, which must be not too good for numerical methods. So how does one go about doing these kind of numerical estimation in functions where gradient is very near to zero away from the solution?


title: Confidence intervals based on the CLT: ever useful?

text: Suppose, for concreteness, that I am trying to estimate the mean of a population using a random sample of size $N$. Many elementary books discuss forming a confidence interval for the population mean by using the central limit theorem to argue that the sample mean is approximately normally distributed. However, the central limit theorem is a theorem about the limit as $N \to \infty$. But if $N$ were really large then the width of the confidence interval would be really small and giving just a point estimate would be good enough. So it seems that we are implicitly assuming that there is a range of $N$ for which the CLT is not too bad an assumption but for which $N$ is not so large as to shrink the confidence interval to practically a point. My question is: is there any basis for making this implicit assumption? Is there any way to judge what this range of useful $N$ is for a particular application?


title: What is the meaning of "All models are wrong, but some are useful"

text: > "Essentially, all models are wrong, but some are useful." > > \--- Box, George E. P.; Norman R. Draper (1987). Empirical Model-Building > and Response Surfaces, p. 424, Wiley. ISBN 0471810339. What exactly is the meaning of the above phrase?


title: How to choose between ANOVA and regression?

text: I am doing a thesis project evaluating differences in attitudes towards trans students. The only manipulation was which vignette participants read. My variables are as follows: DV1: First impressions score (continuous) DV2: Social distance score (continuous) IV: description of the vignette character: transsexual, cross-dresser, or control. I also have several potential predictors. Two of which I do not think I have sufficient data to report on (Sexual orientation and religiosity). The others are gender, level of contact (4 categories), and casual attribution (2 categories). Because I am primarily interested in the differences between groups, I feel as though I should conduct an ANOVA. However, I can't see any straightforward way of incorporating my other predictors into such an analysis. Alternatively, I could run a hierarchical multiple regression to look first at the effect of group and then the other predictors (ordered by theoretical importance). However all my predictors are categorical variables and I'm not really sure how to dummy code to do this properly. Would appreciate any help!


title: Constructing a Bayesian Classifier in R

text: Does anybody know the way to construct a Bayesian classifier in R for two bi- variate gaussian distributions of which the means and variances are known? The two classes are equi-probable and the variables are independent from each other. Any help is deeply appreciated.


title: Difference between marginal and conditional models

text: A _marginal model_ accounts for the correlation within each cluster. A _conditional model_ also takes into account the correlation within each cluster. My questions are:   1. Does a marginal model models main effects across a population whereas a conditional model models main effects within a cluster and across a population?    2. The interpretation of the coefficients of a marginal model is basically the same as "regular model." But what about the coefficients of a conditional model?


title: significant lags in causality test?

text: significant lags in causality test? [Cross-posted to R-Sig-Finance] Hi all, I have a question regarding Granger causality test. I have searched on Internet intensively and extensively and have used both the "grangertest" function in lmtest package and the "causality" function in vars package. If we see that x causes y, is there a way to obtain the significant lags of x which help forecast y? Thanks!


title: Understanding published data: z-ratio for proportions

text: I am trying to understand a published analysis. This is the data of interest:                     D1>0    D1<0     D2>0   7        2     9     D2<0   9        15   24     total  16       17   33      The author notes that 17/33 is 51.5% and states: > "we expect about 50% of the D1's to be negative, and that is what we > actually observe here (z=.08, p=n.s)". Now I assume that the z comparison is 17/33 vs 16.5/33 (chance?) and that the statistic employed is the z-ratio for the significance of the difference between two proportions. However my own calculation for this comes out as z=.123. Can anyone help with the discrepancy? UPDATE: The author has responded and admits that the published z is an error. He gives the corrected answer as z=0.17. I am still not sure how he worked it out though: > We are not considering the difference between two proportions, since they > are not independent 16/33 = 17/33 -1. We are interested whether the > proportion of D1<0 significantly deviated from the expected value of 0.5 The nearest I can get is R prop.test function. Is this directly comparable?               prop.test(c(17, 16.5), c( 33, 33))          2-sample test for equality of proportions with continuity correction     data:  c(17, 16.5) out of c(33, 33)      X-squared = 0, df = 1, p-value = 1     alternative hypothesis: two.sided      95 percent confidence interval:      -0.2411995  0.2715025      sample estimates:        prop 1    prop 2      0.5151515 0.5000000


title: Types of moments used in the method of moments?

text: In Wikipedia, the method of moments uses only ordinary moments: > One starts with deriving equations that related the population moments > (i.e., the expected values of powers of the random variable under > consideration) to the parameters of interest What if the method of moments uses other types of moments, such as central moments or standardized moments? If they are not used or preferred, why? Thanks!


title: Is there justification for using cross validation scores as model averaging weights?

text: Bayesian model averaging uses approximate Bayes factors. Some researchers use AIC to weight models. Is there justification for using, say, the Brier score, median absolute deviation, or other such scores computed from out of sample predictions to construct weights for model averaging instead of the various information criteria? I played around with this a while back when averaging the results of popular election forecast models.


title: Intuitive explanation for dividing in (n-1) when calculating sd?

text: So I've been asked today in class why you divide the sum of square error with $(n-1)$ instead of with $n$, when calculating the sd. I said I am not going to answer it in class (since I didn't wanna go into unbiased estimators), but later I wondered - **is there** an intuitive explanation for this?!


title: How to parameterize the ratio of two normally distributed variables, or the inverse of one?

text: **Problem:** I am parameterizing distributions for use as a priors and data in a Bayesian meta-analysis. The data are provided in the literature as summary statistics, almost exclusively assumed to be normally distributed (although none of the variables can be < 0, some are ratios, some are mass, and etc.). I have come across two cases for which I have no solution. Sometimes the parameter of interest is the inverse of the data or the ratio of two variables. **Examples:**   1. the ratio of two normally distributed variables:      * data: mean and sd for percent nitrogen and percent carbon      * parameter: ratio of carbon to nitrogen.   2. the inverse of a normally distributed variable:      * data: mass/area     * parameter: area/mass My current approach is to use simulation: e.g. for a set of percent carbon and nitrogen data with means: xbar.n,c, variance: se.n,c, and sample size: n.n, n.c:               set.seed(1)     per.c <- rnorm(100000, xbar.c, se.c*n.c) # percent C     per.n <- rnorm(100000, xbar.n, se.n*n.n) # percent N      I want to parameterize ratio.cn = perc.c/perc.n               # parameter of interest     ratio.cn <- perc.c / perc.n      Then choose the best fit distributions with range $0 \rightarrow \infty$ for my prior               library(MASS)     dist.fig <- list()     for(dist.i in c('gamma', 'lognormal', 'weibull')) {         dist.fit[[dist.i]] <- fitdist(ratio.cn, dist.i)     }      **Question:** Is this a valid approach? Are there other / better approaches? Thanks in advance! **Update:** the Cauchy distribution, which is defined as the ratio of two normals with $\mu=0$, has limited utility since I would like to estimate variance. Perhaps I could calculate the variance of a simulation of n draws from a Cauchy? I did find the following closed-form approximations but I haven't tested to see if they give the same results... Hayya et al, 1975 $$\hat{\mu}_{y:x} = \mu_y/mu_x + \sigma^2_x * \mu_y / \mu_x^3 + cov(x,y) * \sigma^2_x * \sigma^2_y / \mu_x^2$$ $$\hat{\sigma}^2_{y:x} = \sigma^2_x\times\mu_y / mu_x^4 + \sigma^2_y / mu_x^2 - 2 * cov(x,y) * \sigma^2_x * \sigma^2_y / mu_x^3$$ Hayya, J. and Armstrong, D. and Gressis, N., 1975. A note on the ratio of two normally distributed variables. Management Science 21: 1338--1341


title: What is the Mean and Standard Deviation of the division of two random variables?

text: I have two normally-distributed independent random variables X and Y and I need to calculate its division Z. As far as I understand the mean of Z is $\mu_Z = \frac{\mu_X}{\mu_Y}$, but I don't know how to calculate the Standard Deviation $\sigma_Z$. Is $\sigma_Z = \frac{\sigma_X}{\sigma_Y}$?


title: Quick Problem with maximizing profit using different distributions

text: I need a little help with a problem I am working on. So here's the situation, a seller can produce a apple for $1.00$ and I need to find the optimal price of selling the apple and expected profit per buyer if the distribution of the value of the apple is EXP(1). So I just need some confirmation on what I am doing. What I have is this: > Expected profit = $(v-1)\cdot (1-F(v))$ > > Optimal price = derivative of expected profit w.r.t. to $v$ set equal to > zero then solve for $v$ I hope this makes sense guys. Like I said I just want make sure I am doing this right or else I am going to have to read the chapter for the 3rd time.


title: Show average instead of median in boxplot

text: When plotting a boxplot with python matplotblib, the lines halfway the plot is the median of the distribution. Is there a possibility to instead have the line at the average. Or to plot it next to it in a different style. Also, because it is common for the line to be the median, will it really confuse my readers if I make it the average (off course I will add a note what the middle line is)?


title: Change detection for beginners

text: I've seen Change detection algorithm - likelihood ratio but I am afraid my question is more basic. I have a sequence $(x_j)_{j=1..N}$ of random observations. I know, these observations are not all from the same random variable, but the random variable changes at several times. I.e. there is a $m$ with $1 < m < N$ and a monotone function $f:\\{1..N\\}\to\\{1..m\\}$ so that $x_j$ is sampled from the discrete distribution $(X_{f(j)})$. All $X_k, k=1..m$ have the same range and only vary in their distribution. My goal is to identify the jumps in $f$ as precise as possible and subsequently to obtain empirical estimates for the unknown distribution of the $X_k$. To make things practically more complicated, the distributions of the $X_k$ seem to be fairly similar. To make a sample, if the $X_k$ were rolls of loaded dice, the probabilites would be more like changing from $(0.17, 0.16, 0.18, 0.18, 0.16, 0.15)$ to (say) $(0.16, 0.16, 0.16, 0.16, 0.18, 0.18)$. The number of rolls might be as few as about 1000 between two jumps of $f$. I might be able to get rough estimates (a few 100 sample indexes wide) for the jumps of $f$ from a different information source, but I'd rather not rely on that. So far we are finding those $f$ jumps manually be educated guessing and then performing G-Tests between the different empirical $X'_k, X'_{k+1}$ to verify a significant change. I'd like to learn about methods which can attack this problem, if feasible.


title: Highly significant positive binary IV: Interaction always positive too

text: I'm working on a project where the main independent variable is binary. We then interact this variable with different continuous variables. The binary variable is by far the most significant variable in the model (t stats ~20) and the variable is always positive. Now, when I introduce an interaction term between this highly significant binary variable and a continuous variable, the interaction is virtually always significant and positive. I'm afraid that the highly significant binary var drives my results. Would this be a valid concern? Perhaps you have some suggestions in how to address this issue.


title: How can I calculate cut-off points from a normal distribution?

text: I'm trying to calculate the upper percentage points for the 0.99 percentile for samples drawn from a normal distribution, with a sample size of 500. How can I calculate the expected values for skewness, kurtosis, studentized range, which would be expected for a drawing of 500 values?


title: Continuous and categorical variables in SPSS GLM

text: I'm currently writing the results part of my thesis and I'm extremely bad at statistics. My variables are   * Dependent variable: Response time (with two levels)   * Factor: manipulation variable (exposing to either one screen or another screen).   * Covariate: 7 point scale. When checking assumptions I found an interaction between the covariate and the independent/factor, resulting in violating of the homogeneity of the slopes. My teacher said that I have to use the manipulation variable as a categorical variable and the covariate as a “continuous variable” in a multivariate regression analysis in the general linear model (GLM; full factorial). She gave my other fellow students a syntax but forgot me. Is there somebody who knows have I can get the covariate to be a continuous variable and how to conduct a GLM (so not an ANOVA)?


title: Multiple regression with categorical variables

text: I have a data set with 3 continuous variables and 3 categorical variables. I know that I have to create dummy variables for categorical variables but I am wondering if it is a must to do that for all 3 categorical variables in this case?


title: If linear regression is related to Pearson's correlation, are there any regression techniques related to Kendall's and Spearman's correlations?

text: Maybe this question is naive, but: If linear regression is closely related to Pearson's correlation coefficient, are there any regression techniques closely related to Kendall's and Spearman's correlation coefficients?


title: Has anyone used the Marascuilo procedure for comparing multiple proportions?

text: The Marascuilo procedure as described here seems to be a test that addresses the issue of multiple comparisons for proportions when you want to test which specific proportions are different from each other after rejecting the null in an overall chi-square test. However, I am not very familiar with this test. So, my questions:   1. What nuances (if any) should I worry about when using this test?   2. I know of at least two other approaches (see below) to address the same issue. Which test is the 'better' approach?:     * Performing "partitioned chi square" refereed in this answer by @Brett Magill      * Using a Holm–Bonferroni method to adjust p-values.


title: Can someone explain the concept of 'exchangeability'?

text: I see the concept of 'exchangeability' being used in different contexts (e.g., bayesian models) but I have never understood the term very well.   1. What does this concept mean?   2. Under what circumstances is this concept invoked and why?


title: Probability understanding 52 C 5

text: Hi I have a problem undertanding this: Example: In a deck of 52 cards, 5 cards are chosen.               What is the probability that all 5 cards have diﬀerent face values?     total number of outcomes = 52 C 5      total number of face value combinations =      13 C 5      total number of suit possibilities, with replacement = 4^5      Now I get the 52 C 5, the total number of possibilities, and I get the 13 C 5 the number of face cards in the deck, so:                P( 5 diff face values) = (13 C 5)(4^5)/(52 C 5)      Why? What exactly does the 4^5 is doing there? Help.


title: Is it possible to have xerror increased in a tree using rpart?

text: I am new to R and rpart package. When I plot the tree using rpart:               > temp_control <- rpart.control(xval=10, minbucket=2, minsplit=4, cp=0.0001)     > dfit <- rpart(Target~., data = temp_data, method = 'class', control=temp_control)     > printcp(dfit)      Then I get :                         CP nsplit rel error xerror     xstd     1 0.00189329      0   1.00000 1.0000 0.040140     2 0.00172117     28   0.92255 1.0861 0.041708     3 0.00114745     32   0.91566 1.0947 0.041861     4 0.00098353     41   0.90534 1.1102 0.042133     5 0.00086059     48   0.89845 1.1274 0.042433     6 0.00043029     62   0.88640 1.1515 0.042849     7 0.00034423     75   0.87952 1.1635 0.043055     8 0.00028686     80   0.87780 1.1687 0.043142     9 0.00010000     89   0.87263 1.1807 0.043346      Why does xerror increase with the growth of a tree? Do I need some more adjustment of the parameters? Also, I am wondering how the root node error is calculated. Is it only related to a certain dataset? Does it have any relationship with parameters setting? Actually, I tried "anova" method although my response variable is categorical (Y/N). I just changed them to 0/1 and run "anova", then I can get :                        CP    nsplit rel error  xerror     xstd     1   3.1473e-02      0   1.00000 1.00025 0.037408     2   1.1506e-02      1   0.96853 0.97164 0.035479     3   5.6396e-03      2   0.95702 0.96528 0.035172     4   4.6137e-03      3   0.95138 0.96970 0.035029     5   4.4412e-03      6   0.93754 0.97246 0.035019     6   4.3751e-03      7   0.93310 0.97006 0.034915     7   4.1352e-03     10   0.91997 0.97109 0.034912     8   3.5702e-03     11   0.91584 0.97316 0.034847     9   3.0148e-03     14   0.90513 0.96819 0.034671     10  2.5334e-03     15   0.90211 0.96872 0.034725     11  2.2789e-03     16   0.89958 0.96959 0.034753     12  2.2342e-03     17   0.89730 0.97437 0.034829     13  1.8732e-03     18   0.89507 0.98647 0.035104     14  1.8401e-03     19   0.89319 0.99511 0.035199      Anyone has any idea about this?


title: Does lm() use partial correlation - R Squared Change?

text: I come from an SPSS background and am attempting to move to `R` for it's superior flexibility and data manipulation abilities. I have some concerns however as to whether the `lm()` is really using partial correlations. I'm basically trying to run a linear regression, using something similar to the "enter" setting in SPSS, which essentially builds the model one variable at a time, reporting the change in $R^2$ with each additional variable. This allows you to determine how much predictive power each variable adds to the model. When I run the same analysis in `R` however, I don't get any information on the $R^2$ contributed by individual variables, and I'm not even sure that it's using partial corrrelations to calculate the p-values it's reporting! My code follows:               summary(m1 <- lm(totalprop ~ cos(Angle) + Alignment + colour +        Angle*Alignment, dataset))      My questions:   1. Does R use partial correlations to determine reported p-values from `lm()`?   2. How can I make `R` report change in $R^2$ with each additional variable.   3. How can I make `R` act like SPSS in calulating the model piece by piece? Is this possible without running multiple iterations of the lm() function? If not, how does one control for the effects of covariates in R?


title: Missing data and SVM

text: I know that in multiple imputation, one runs a regression for each imputed data set. The final regression coefficients are averages of the coefficients obtained from the imputed data sets. Can you do the same thing with support vector machines? That is, can you run a support vector machine model on each of the imputed data sets and then get a final "averaged" result?


title: How does the standard error work?

text: I have been looking into the inner-workings of the standard error recently, and I found myself unable to understand how it works. My understanding of the standard error is that it is the standard deviation of the distribution of sample means. My questions are: • how do we know the standard error is the standard deviation of the sample means when we usually take just a single sample? • why doesnt the equation to calculate the standard error mirror the standard deviation equation for a single sample?


title: Standard error for expected value

text: I have a computer simulation from which I get a sample of values of some microscopic quantity $X$, i.e. $\\{ x_1,\ldots,x_N \\}$. I'm interested in estimating the expected value of $X$, i.e. $E[ X]\approx\frac{1}{N}\sum x_i$, and the variance associated with this **estimate**. Now, do I understand correctly that:   1. taking the estimate of the variance of $X$, $Var[X]\approx\frac{1}{N}\sum (x_i-\bar{x})^2$, is NOT what I'm looking for (which in this notation would be $Var[E[X]]$)?   2. I need to make use of some statistical technique, such as bootstrapping?


title: Standard error, margin of error, polls and binomial distribution

text: So, when doing polling, we use binomial distribution to capture the data (in case there are only two choices.). My question is, let us say that estimated percentage according to the poll that supported Obama was 47%. Then we would be able to set $p = 0.47$ and calculate standard error when given the number of samples, $n$. But the true standard error will not be this. So, is it OK to use estimated standard error right away?


title: Very easy question : Sample Variance

text: > **Possible Duplicate:**   >  How does the standard error work? Let's say I have an unbiased estimator, which is the sample mean, to estimate $\theta$. When can you say that $$VAR[\overline{X}] = \frac{\sigma^2}{n} $$ ? If so yes, if not: when is this true?


title: What is the estimate of sample mean variance?

text: > **Possible Duplicate:**   >  How does the standard error work? What is the estimate of sample mean variance and how do you get it? Are my understandings of the following correct? $var(Y)=\frac{1}{N} \sum (X_i-\mu_Y)^2$ population variance, how much the individuals vary from the mean. $var(\hat Y)=\frac{1}{n-1}\sum(x_i-\bar y)^2$ sample variance, how much the individuals of the sample vary from the sample mean. $var(\bar y)=\frac{var(Y)}{n}$ sample mean variance, how much the sample mean can vary from the true mean. $var(\hat{\bar y})=$ estimate of sample mean variance $Y={X_1...X_N}$ is a population with N individuals $y={x_1...x_n}$ is a sample of the population with n observations $\mu$ is the mean $\bar y=\sum x_i/n$ is the mean of the sample.


title: two-sample Kolmogorov-Smirnov test in R p-value confusion

text: I'm confused about the appropriate interpretation of p-values returned by the two-sample Kolmogorov-Smirnov test (ks.test) in R. In slide 23 of this presentation about non-parametric two-sample tests, the author states that when analyzing the ks.test results:               ks.test(male, female)     Two-sample Kolmogorov-Smirnov test     data: male and female      D = 0.8333, p-value = 0.02597      the p-value > needs to be multiplied by 2 for a 2-tail test. Thus, P = 0.05194 Is that true? If we used the original p = 0.02597, we would reject the hypothesis that the distributions similar, because p < 0.05, correct? Whereas if we multiply it by 2, the p would suggest that there is no difference between distributions, since p > 0.05? What am I missing?


title: Smoothing dirty data?

text: I'm trying to do analysis on revenue data and due to circumstances beyond my pay grade, there are large spikes in the data around the time financial reports are due, followed by corresponding low values immediately after. Is there a preferred algorithm or methodology to compensate for this and get more meaningful data? It would seem that a simple moving average would be inappropriate for the task, given that I know there are discontinuities at time edges. Is there a better smoothing strategy?


title: Why use chi-square (or other) statistic over regression modeling for 2x2 contingency tables?

text: I have a question regarding the use of logistic/log-linear regression vs. contingency test statistics, such as chi-square. Can someone explain to me why it would ever be preferable to use test statistics over logistic (or log-linear) regression modeling?


title: Definitions of Prediction vs. Predictor

text: I am writing an article which includes discussion of the MMSE estimator of the posterior predictive distribution. Since I use this term quite frequently, I was considering referring to this estimator as the prediction, however, I am thinking that the prediction connotes an estimate rather than an estimator: For the difference between the two, wikipedia states: In statistics, an estimator is a rule for calculating an estimate of a given quantity based on observed data: thus the rule and its result (the estimate) are distinguished. On the other hand, I considered using the word predictor, however, this word seems to connote a variable that is used to help predict another variable, not an estimation function. Any suggestions on which term (prediction/predictor) is more reasonable to use to denote such response estimators?


title: What are some methods for generating simulated time series data for use in modeling?

text: I have a data set which consists of 50 observed years for which I have date and inflow values between a river and a reservoir. The data is formatted as follows:               Day     Inflows     123      4356.0      I have 50 years of observed data which includes daily values for each day. The plotted data looks as follows: ![Plot of inflow values for each day of 2005](http://i.stack.imgur.com/DYCQH.png) Ignoring the small outlying patches of green for now, I'm looking for a method of simulating the green portion of the curve (which represents the spring runoff period) so that I can do a series of model runs with simulated spring inflows. I'm hoping that someone might be able to point me towards a 'standard' method for approaching this kind of problem. There is extensive variation between the years in the data set (drought years show very little change between the blue and green portions of the year). What are some methods for approaching generating data like this? Let me know if there is any other information that I can provide. I'm quite new to this kind of problem so I'm unsure what would be most useful.


title: Is there a way to "regularize" the Johansen cointegration test?

text: I'm using a Johansen cointegration test to check for cointegration among a large number of time series. I've found that while the eigenvectors look great in-sample, the cointegrating relationships often don't hold up out of sample. However, if I manually "regularize" the eigenvectors by setting the smaller weights to zero, then the OOS performance is much better. Is there a more robust or principled way to perform this regularization?


title: How to select a subset of variables from my original long list in order to perform logistic regression analysis?

text: My situation:   * small sample size: 116   * binary outcome variable   * long list of explanatory variables: 44   * explanatory variables did not come from the top of my head; their choice was based on the literature. Statistical test chosen: logistic regression I need to find the variables that best explain variations in the outcome variable (I am not interested in making predictions). The problem: This question is a follow-up on the 2 questions listed below. From them, I got that performing automated stepwise regression has its downsides. Anyway, it seems that my sample size would be too small for that. It seems that my sample is also too small to enter all variables at once (using the SPSS 'Enter' method). This leaves me with my issue unresolved: how can I select a subset of variables from my original long list in order to perform multivariate logistic regression analysis? UPDATE1: I am not an statistician, so I would appreciate if jargons can be reduced to the minimum. I am working with SPSS and am not familiar with other packages, so options that could be run with that software would be highly preferable. UPDATE2: It seems that SPSS does not support LASSO for logistic regression. So following one of your suggestions, I am now struggling with R. I have passed through the basics, and managed to run a univariate logistic regression routine successfully using the glm code. But as I tried glmnet with the same dataset, I am receiving an error message. How could I fix it? Below is the code I used, followed by the error message:               data1 <- read.table("C:\\\data1.csv",header=TRUE,sep=";",na.string=99:9999)          y <- data1[,1]          x <- data1[,2:45]          glmnet(x,y,family="binomial",alpha=1)            **in lognet(x, is.sparse, ix, jx, y, weights, offset, alpha, nobs,  :      (list) object cannot be coerced to type 'double'**      UPDATE3: I got another error message, now related to missing values. My question concerning that matter is here.   * Can univariate linear regression be used to identify useful variables for a subsequent multiple logistic regression?   * Algorithms for automatic model selection


title: Given the moment generating function of $X$, find the distribution of $X$

text: $X$ has MGF $m_{x}(t)=e^{2t}$ for any $t \in \mathbb{R}$. What is the distribution of $X$? My instructor went about it like this: $m_{x}(t)=E(e^{xt})=Ee^{2t}+e^{2t}$ Here I am really confused, how did $E(e^{xt})$ become what's on the right hand side of the equation? I wonder if I copied it down wrong...but can someone please explain what he was trying to get at here? So $P(X=2)=1$ Also, this is assuming that $X$ is a continuous variable, correct? How can I tell from the problem that $X$ is continuous, not discrete?


title: Using random forest in MATLAB

text: I am having issues in using random forests in MATLAB. I have features of size 2000 and around 4000 data points. I am trying to learn how to compute random forests in MATLAB using the library Random Forest. However, I guess it's too slow to use this method. I have set the number of trees to 500 and `mtry` to 720 and it is taking ages. Is it natural to take this long? How long does it normally take?


title: Neural Nets: cost functions

text: Multi-Layer Neural Nets' implementation of one of the commercial packages contains the following cost functions: squared error, cross-entropy, maximum likelihood and perceptron convergence. While first 2 of them are very well discussed in the literature, the last two are not so clear. How max. likelihood is defined and what is "perceptron convergence" cost function?


title: Kalman filtering confusion

text: Say I have a very simple model \begin{align} &C_t = C_{t-1} + z_t\\\ &X_t = C_t + v_t \end{align} everything is 1-dimensional, $z_t$ and $v_t$ are are uncorrelated and white noises, and $X_t$ is the observed value. I generate a simulated data sequence with some predefined paramenter for variance, then I use Kalman filter with the same parameter to do the filtering. Then I check the output $C_t$, it doesn't look like anything like random walk, $dC_t = C_t - C_{t-1}$ are serially correlated. What do I miss here? Or rather, is the filtered sequence necessarily satisfy model assumption? The following matlab code implement the algorithm               verror = 2;     vobjs = 3;     N = 10000;     x = randn(1,N)*sqrt(verror);     n = randn(1,N)*sqrt(vobjs);     obs = cumsum(x) + n;     %plot(obs)     filtered = zeros(size(x));          filtered(1) = 0;     s0 = verror;     for i = 2 : N         Kgain = (s0+verror)/(s0+vobjs+verror);         filtered(i) = (1-Kgain)*filtered(i-1) + Kgain*(obs(i));         s0 = (1-Kgain)*(s0+verror);     end     var(diff(filtered(1000:end)))     var(filtered(1000:end) - obs(1000:end))     % The filtered sequence is 0 auto-correlated.     autocorr(filtered(1000:end) - obs(1000:end))      The code is now fixed. there is no auto-correlation. What about variance? The variances are 2.0226 and 1.3694 respectively.


title: Output of a system with a white Gaussian process as an input

text: I want to present a question which seems to predict a theorem I know. Hence I guess I'm missing something and would be happy to understand what I'm missing. Here's is the system: ![enter image description here](http://i.stack.imgur.com/lvYXg.png)   The signals $ {X}_{1}, {X}_{2} $ are Independent White Gaussian Noise Process.   The filters $ {H}_{1}(f), {H}_{2}(f) $ are identical, namely, $ {H}_{1}(f) = {H}_{2} $.   The phase, $ \theta \sim U[0, 2\pi] $ and independent of any other parameter or signal. The questions are   1. Are $ {Y}_{1}(t), \ {Y}_{2}(t) $ Gaussian Processes?   2. Is $ r(t) $ Gaussian Process? The answers, according to a book are $ {Y}_{1}, \ {Y}_{2} $ aren't Gaussian Processes and $ r(t) $ is a Gaussian Process.   Yet $ {Y}_{1}, \ {Y}_{2} $ are obviously independent (Edit: They are dependent, see answers and comments below) and according to Cramer's Theorem:   http://en.wikipedia.org/wiki/Normal_distribution#Infinite_divisibility_and_Cram.C3.A9r.27s_theorem   If the sum of two independent R.V. is Gaussian then each of them must be Gaussian.   Yet it's pretty clear to me that Y1 / Y2 aren't Gaussian.   I wrote the following MATLAB Simulation:               numRealizations = 1e4;     numSamples      = 1e3;          h1Coef = ones(10, 1) / 10;     h2Coef = h1Coef;          harmonicFreq    = 10; %<! [Hz]     timeVector      = [0:(numSamples - 1)];     randomPhaseVec  = 2 * pi * rand(numRealizations, 1);          % White Noise is Ergodic     % Each Row is a realization (Over the Time Vector)     x1SignalMatrix  = randn(numRealizations, numSamples);     x2SignalMatrix  = randn(numRealizations, numSamples);          % Filtering along the columns     z1SignalMatrix = filter(h1Coef, 1, x1SignalMatrix, [], 2);     z2SignalMatrix = filter(h2Coef, 1, x2SignalMatrix, [], 2);          y1SignalMatrix = z1SignalMatrix .* cos(bsxfun(@plus, (2 * pi * harmonicFreq * timeVector), randomPhaseVec));     y2SignalMatrix = z2SignalMatrix .* sin(bsxfun(@plus, (2 * pi * harmonicFreq * timeVector), randomPhaseVec));          rSignalMatrix = y1SignalMatrix + y2SignalMatrix;          timeIndex           = randi([1, numSamples], [1, 1]);     realizationIndex    = randi([1, numRealizations], [1, 1]);          figure;     hist(x1SignalMatrix(:, timeIndex), 500);     figure;     normplot(x1SignalMatrix(:, timeIndex));               figure;     hist(z1SignalMatrix(:, timeIndex), 500);     figure;     normplot(z1SignalMatrix(:, timeIndex));               figure;     hist(y1SignalMatrix(:, timeIndex), 500);     figure;     normplot(y1SignalMatrix(:, timeIndex));          figure;     hist(y1SignalMatrix(realizationIndex, :), 500);     figure;     normplot(y1SignalMatrix(realizationIndex, :));               figure;     hist(rSignalMatrix(:, timeIndex), 500);     figure;     normplot(rSignalMatrix(:, timeIndex));          figure;     hist(rSignalMatrix(realizationIndex, :), 500);     figure;     normplot(rSignalMatrix(realizationIndex, :));      This is the distribution of Y1 over the Ensemble: ![The Distribution of Y1 over the Ensemble](http://i.stack.imgur.com/HQx70.png) Which is clearly not Gaussian as I expected. This is the distribution of Y2 over the time (One Realization): ![The Distribution of Y2 over Time - 1 Realization](http://i.stack.imgur.com/OiWdr.png) Now, the results for $ r(t) $: The Distribution of the Ensemble: ![enter image description here](http://i.stack.imgur.com/FB4BP.png) The Distribution of 1 Realization: ![enter image description here](http://i.stack.imgur.com/HWZK8.png) Well, for Y1, I can see why the ensemble won't Gaussian as for each 't' the sample is multiplied by a different factor hence each of them, though is Gaussian, is coming from a different distribution.   I even wrote a proof why it doesn't comply with Gaussian Distribution properties (The fourth Moment) - http://imgur.com/l017RrR. Yet I can't get how $ r(t) $ is Gaussian.   What am I missing? Thank You.


title: Converting R data frame to mids object with as.mids returns - error ini$imp[[i]] : subscript out of bounds

text: I am currently working on a project addressing individuals' sexual partnering behavior. There are missing data - which I impute in R using mice. The imputations go off without issue. The trouble begins when I need to ensure that a few impossible values are not present after the imputations. Specifically, in same-sex sexual partnering encounters males cannot have a positive response on the variable 'occurrence of vaginal intercourse.' For a number of reasons the analyses for the project are divided by sex from the outset. So here I am only dealing with males (as in the starting data frame has already been subset by sex, for which there were no missing values). Following the imputations, this anatomical impossibility shows up in a few places. My thoughts were to simply edit the values since there are very few (a couple across 5 imputed data sets with several thousand cases). I convert the mids object to a data frame, edit the values manually, and then attempt to convert back to a mids object using as.mids. At which point I get the message:               Error in ini$imp[[i]] : subscript out of bounds      The original data frame I'm dealing with at this point ("data" consisting of only male respondents) contains the following variables (plus a lot of others that aren't relevant here) -               respondentID, age, race, samesexENC, mutualmast, oralsex, vaginalsex, analsex, enjoy      After using `mice()` to do the imputations, the `mids` object is `dataIMP`. At which point I do the conversion to a data frame using `complete()`. I then subset down to identify the problematic rows and correct `vaginalsex` to "No".               # convert to a data frame in order to do edits     EDITthis <- complete(dataIMP,action="long",include=TRUE)     #      # subset to identify samesexENC that included vaginalsex     q <- subset(EDITthis, vaginalsex=="Yes")     q <- subset(q, samesexENC=="Yes")     #      # edit the cells in the frame directly so      # that these equal "No" on vaginal (which is column 29)     EDITthis[5955, 29]="No"       EDITthis[19656, 29]="No"      We're all clear at this point. Pretty simple stuff even for my beginner level R fluency. Then I want to convert back to `mids` and do the analyses, so I use `as.mids()` \-               # convert back to mids object     dataFINAL <- as.mids(EDITthis)      This is where my full-steam-ahead day of data analysis came to a halt yesterday with the previously mentioned error message. I first went directly to the mice documentation. Looking at the arguments to `as.mids()` there doesn't seem to be any issues. The defaults for the `.id` and `.imp` arguments correspond to the correct columns in `EDITthis`. I went back to the documentation for `complete()` as well and cannot find any issues there. Finally I called to `as.mids` to eyeball the code behind the function. I believe I understand what's going on there, but nothing tipped me off to any potential issues in my code. Finally, as a safety measure, I re-installed the package all together, walked back through the entire process from the beginning, and the results were exactly the same - everything is fine up to the point where I want to turn the `EDITthis` data frame back into a `mids` object. Probably a comma or quotes off somewhere, but I can't ferret it out. Thanks in advance!


title: Are there any alternatives to simulation for determining the distribution of number of events from two dependent non-homogeneous Poisson processes?

text: A "state of the art" model for the distribution of goals scored in a soccer match is that of Dixon and Robinson (1998) "A Birth Process Model for Association Football Matches" which accounts for two key phenomenon: 1) More goals are scored at the end of matches than at the start (hypothesised to be due to the fatigue suffered by both teams) 2) Scoring rates depend upon the current score line for a myriad of reasons such as teams with a lead becoming complacent or teams preferring to play out for a draw rather than risk a loss by going for the win The model assumes that the goals scored by the home and away teams in a match follow non-homogeneous Poisson processes. Let $t$ denote the time elapsed in a match, normalised to fall between $0$ and $1$, the $x$-length vector $\vec{t_H}$ denote the times at which the home team scored goals and the $y$-length vector $\vec{t_A}$ denote the times at which the away team scored goals. The likelihood for the match is then $$ L(\vec{t_H},\vec{t_A}) = \exp\left(-\int_0^1 \lambda(t) dt\right)\frac{\prod_{i=1}^{x} \lambda({t_{H}}_i)}{x!}exp\left(-\int_0^1 \mu(t) dt\right) \frac{\prod_{j=1}^{y} \mu({t_{A}}_j)}{y!} $$ where $\lambda(t)$ is the scoring rate for the home team at time $t$ dependent on a combination of time homogeneous factors (e.g. home team attacking ability versus away team defending ability, home advantage) and time inhomogeneous factors (e.g. score line at time $t$). Similarly for $\mu(t)$. The two processes are dependent because when a team scores the score line changes and the scoring rates are themselves score line dependent. The likelihood can easily be evaluated by carrying out the integration in the exponent numerically. Hence it is straightforward to calculate the parameters of the model (team abilities, home advantage, time effect, score line parameters etc.) via maximum likelihood. In terms of prediction, obvious quantities of interest are:   * $P(x > y)$: home team wins   * $P(x < y)$: away team wins   * $P(x = y)$: draw   * Probability of particular score lines, e.g. $P(x=1,y=0)$   * Probability of total goals in the match, e.g. $P((x+y) < 2.5)$ To calculate these quantities (approximately) given a set of model parameters, we could use Monte Carlo methods to generate matches according to these processes and then calculate the frequencies of each final score. Simulating from the processes is relatively straightforward, by generating goals from a single enveloping homogeneous Poisson process in conjunction with rejection sampling and then distributing them to the home or away team accordingly. The drawback to this approach is, obviously, the computational burden of Monte Carlo simulation. Consider attempting to make predictions in real-time as matches are being played, of which there may be many happening simultaneously, and it quickly becomes a cause for concern. My question, therefore, is whether there are any alternative approaches we can consider which do not incur such as high computational cost (even if they rely on an approximation that sacrifices accuracy for ease of calculation)? * * * _For clarity, I am not looking for (basic) suggestions on how to efficiently implement the Monte Carlo simulation which I have already written in multi- threaded C, uses quasi-random numbers which have been pre-generated using unrolling and exploits piecewise thinning to achieve a very high acceptance rate. If you think there is still scope for a dramatic performance increase then of course I am all ears but really I am looking for a fundamentally different approach!_


title: How is the formula for the Standard error of the slope in linear regression derived?

text: As stated in many textbooks, the Standard error of the slope in linear regression with one variable is $\sqrt{\frac{s^2}{SSX}}$ or some rewrite, ${s^2}$ being the error variance and ${SSX}$ being the sum of the ${x}$-squares. Can anybody help with an explicit proof?


title: Meaning of 'independence' in regression model

text: I am familiar with the following condition for statistical independence where $x$ and $y$ are two different (random) variables: \begin{equation} E(xy)=E(x)\times E(y) \end{equation} Now, I came across a paper (Appendix B, p.1311) where it mentions that if $h$ and $^{\partial h}/_{\partial z}$ are independent then following holds: \begin{equation} E\left(h \cdot \frac{\partial h}{\partial z}\right)=E(h)\times E\left(\frac{\partial h}{\partial z}\right) \end{equation} I would like to know: What does it mean if $h$ is a dependent, continuous variable--such as smoking expenditure (in dollars)--and $z$ is a continuous independent variable--say income in dollars--in regression model. Note that $^{\partial h}/_{\partial z}$ is the partial derivative of $h$ with respect to $z$ which is equal to $dh/h$ (percentage change in $h$)


title: Identifying sequential patterns

text: I am working with sequence data which are long lists of malware win-api calls. I am trying to cast the problem of identifying 'malware behavior' into one of finding sequential patterns. I treat each api call as a single item Itemset. The number of different possible items (api calls) is quite large. Now, when I apply the SPADE algorithm (see also, Zaki, SPADE: An Efficient Algorithm for Mining Frequent Sequences, Machine Learning, 42, 31–60, 2001) I run into memory problems. Is there a better alternative way to find sequential patterns among large high vocabulary sequences?


title: How accurate sum of kernel function needs to be, so that we can use it in Mean shift algorithm (may be for image segmentation)?

text: Mean shift is a procedure for locating the maxima of a density function given discrete data sampled from that function. It is useful for detecting the modes of this density. This is an iterative method. How accurate sum of kernels needs to be so that this algorithm works? For example - Is it okay if it performs poorly at tails and better at mode?


title: How can I compare the relative effects between values of a categorical IV?

text: I am a social science phd student trying to figure out a statistical test for small conference presentation I am working on. However, I realized I don't know how to run the model I want. Some background on my question: **DV:** Change in a person's weight (continuous variable) **IV:** Which fruit someone ate for breakfast (categorical variable with four categories: Guava, raspberries, Pears, and strawberries) **Statistical Method:** Simple OLS Regression According to my theory, eating guavas should have the greatest effect on someone's weight compared to the other fruit, then raspberries, then pears, and with strawberries having the least effect on someone's weight. To test my theory, my original plan was to run a simple OLS regression with 3 binary variables (one for each fruit, with strawberries being the baseline). However, I just realized that that sort of model would only tell me if each fruit had a greater effect than strawberries, but it would not tell me that the effect of guava was greater than raspberries, or if raspberries had a greater effect than pears, etc... Is there a simple method in an OLS regression (in STATA) that would allow me to test if guava's had the greatest effect, raspberries had a second greatest effect, pears had the third greatest effect, and strawberries had the least effect? Thanks! FYI: No my paper is not really about fruit and weight loss. Those variables were made up for the purposes of this question.


title: Incorporating new data in classification

text: I am trying to predict class membership (classification) and have a database of 40k observations (people) with 94 attributes. I have also asked a thousand people from the database to fill in a more extensive survey and have 400 further attributes for them. I don't want to waste the time of everyone on the database with filling in the entire survey. I'd like to work out for each person what questions would give me a reasonable estimate of their class membership. Essentially I want to build a decision tree for every person where the nodes in the decision tree are the questions to ask. Because values for the 94 attributes will be different from person to person I would expect everyone's decision tree to be different, but because some questions are much more predictive of the class I would expect most trees to be similar. Any help or pointers appreciated!


title: Independence test

text: I'm working on Bayesian Network and I need to find a broad range of statistical test for testing independence and conditional independence between 2 variables with a potential conditioning set of important size. My data is a mix of normal (not the most common), non normal (the most common), continuous (the most common), discrete, with dependences being not linear. So far I was using Z-fisher test as I found a nice implementation in the `bnt` toolbox for MATLAB (I'm not so good in developing) but it assumes linearity and normality which are really heavy assumptions. I found two or three implementation of the Hilbert Schmidt Independence Criteria but, unfortunately, they perform quite poorly. Do you have some advice? Pointers? I would like to design a kind of super class of test which will have access to a specific test depending on the nature of the parameters (testing continuous independent of discrete conditional on a continuous variable is still a bit unclear for me).


title: Critical effect sizes and power for paired t test

text: Is there some sort of easy way to calculate the required sample size for a paired t-test with an effect size of 2 standard deviations? In my line of work, we often use an effect size of 2SD from the mean for power and sample size calculation with independent t-tests but I'm not sure how to adapt that to a paired t-test. I'm not even sure if that typed of effect size makes sense for a dependent t-test. Thanks for the help.


title: Difference between model / line?

text: I have a dataset. I'm asked to write out the linear regression model. I'm also asked to calculate the least squares regression line. What's the difference?? [[EDIT]]


title: One sided estimates

text: I am looking for a method to estimate an unobservable process in the following set up: at some moments of time (which I don't control) I am getting an observation, that the value of the process is above (or below) some value (which is also no in my control). My goal is assuming that the unobservable process is stationary, get an estimate of its current value. Thanks a lot!


title: RMSE vs. Coefficient of Determination

text: I am evaluating a physical model and would like to know which one of the methods I should be using here (between RMSE and Coefficient of Determination R2) The problem is as follows: I have a function that outputs predictions for the input value x, $\overline{y_x}= f(x)$. I also have the actual observation for that value that I call $y_x$. My question is what are the pros and cons of either RMSE or $R^2$. I have seen both of them being used in papers for the problem that I am working on.


title: Which kind of statistical exact test should I use?

text: I am testing a new kind of medicine against an older medicine in a very limited population. All the test subjects will first recieve the older, ineffective medicine. They will only recieve the new medicine if the old medicine was ineffective. My data so far: None of 4 test subjects were cured by the old medicine, so all were given the new medicine. All 4 test subjects were cured by the new medicine. There are only two outcomes: cure or no cure. I guess I cannot use a Fisher's exact test for this study, as the data are paired. Which test(s) could I use for this kind of experiment?


title: Ratio of Range to IQR vs. Coefficient of Variation -- which is the more useful robust measure?

text: For a given set of data, spread is often calculated either as the standard deviation or as the IQR (inter-quartile range). Whereas a `standard deviation` is normalised (z-scores, etc.) and so can be used to compare the spread from two different populations, this is not the case with the IQR since the samples from two different populations could have values at two quite different scales,                e.g.       Pop A:  100, 67, 89, 75, 120, ...      Pop B:  19, 22, 43, 8, 12, ...      What I'm after is a robust (non-parametric) measure that I can use to compare the variation within different populations. Choice 1: `IQR / Median` \-- this would be by analogy to the coefficient of variation, i.e. to $ \frac{\sigma}{\mu}$. Choice 2: `Range / IQR` Question: Which is the more meaningful measure for comparing variation between populations? And if it is Choice 1, is Choice 2 useful for anything / meaningful, or is it a fundamentally flawed measure?


title: Gaussian Processes: How to use GPML for multi-dimensional output

text: Is there a way to perform Gaussian Process Regression on multidimensional output (possibly correlated) using GPML? In the demo script I could only find a 1D example. A similar question on CV that tackles case of multidimensional input. * * * I went through their book to see if I could find anything. In the 9th chapter of this book (section 9.1), they have mentioned this case of multiple outputs. They have mentioned a couple of ways to deal with this, One - using a correlated noise process and Two - Cokriging (Correlated prior). I still don't know, how I can incorporate any of these ideas into the GPML framework. * * * Also, are there any other GP libraries/frameworks that support multi- dimensional output?


title: Reference for dummy variable regression for repeated measurements

text: I want to analyze some repeated measures data using a dummy variable for each measurement unit. This is one approach to deal with the fact that the data is not i.i.D. when repeated-measurements are present (in addition to classical repeated-measures ANOVA, multivariate ANOVA and mixed modeling, see Keselman & Algina, 2001), but I am missing a reference for that. Therefore my question is:   * **What is a reference for coding subjects as a fixed effect in repeated-measures analysis?** Keselman, H. J., Algina, J., & Kowalchuk, R. K. (2001). The analysis of repeated measures designs: A review. British Journal of Mathematical and Statistical Psychology, 54(1), 1–20. doi:10.1348/000711001159357


title: Prediction and explanation of user rating based on multiple criteria

text: I'm trying to figure out a way to both predict how a user would rate a certain document, as well as an explanation of why certain documents are rated a certain way. A user is represented by:   * favorite topics;    * average length of the documents the user read;    * the distribution of the times of day a user reads the documents;    * location;    * age. With this information I want to:   1. Predict how much a user will like a document (based on the document's topics, length, the time of day, etc)   2. Get explanations like: documents with this length are usually read around this time; or, documents about these topics are usually read at these locations. I hope this question is not too broad. I just need some starting tips.


title: How can I simulate a normal distribution from means and 95% confidence limits?

text: Well, I have these values: $\mathrm{Mean}=5.77$; $95\%\mathrm{CL}=5.30-6.11$ How can I simulate a normal distribution using these parameters? I do not have standard deviations. How can I obtain random values (for example, 1000 of them) that fit to these parameters? (like using the `rnorm` function when having Standard deviations except $95\%CLs$)


title: How to assess predictive power of set of categorical predictors of a binary outcome? Calculate probabilities or logistic regression?

text: I am trying to determine if simple probabilities will work for my problem or if it will be better to use (and learn about) more sophisticated methods like logistic regression. The response variable in this problem is a binary response (0, 1). I have a number of predictor variables that are all categorical and unordered. I am trying to determine which combinations of the predictor variables yield the highest proportion of 1's. Do I need logistic regression? How would it be an advantage to just calculating proportions in my sample set for each combination of the categorical predictors?


title: Multivariate ordered logit or probit

text: I have two ordinal dependent variables, each having three response levels. You can use an ordered logit or probit model for such data if you have one dependent variable. I've seen some papers about multivariate ordered regression, and wonder if there are prepackaged functions in any of the usual stats software environments to do this. I am most proficient in R and Stata. Thank you.


title: Accuracy of each path in a decision tree in Matlab

text: I would like to know the accuracy of each path in a decision tree in Matlab. I can build a decision tree in Matlab by:               ctree = ClassificationTree.fit(X,Y)      Where `X` is a matrix of instances (rows are observations, columns are variables), and `Y` is the response variable. Here's an example of a simple output tree: ![Simple example](https://lh5.googleusercontent.com/-Dij9WpLXnAc/T9jlc1rli7I/AAAAAAAALjg/Lb6Kq_90UVs/s800/simple_tree.png) For example, the above tree has three different paths:   1. (x1 < .5) => class 0   2. (x1 >= .5) And (x2 < .5) => class 0   3. (x1 >= .5) And (x2 >= .5) => class 1 Now I would like to get the accuracy/confidence of each path. Is there a way to do that? For example, I want to know how many instances were classified as `0` **through the first path** and they were actually class 0 (I can provide the same training set `X`,`Y` for that) If I use the function `resubLoss`:               resuberror = resubLoss(ctree)      I will get the overall accuracy by resubmitting the same training set. I'm confused why I could not find what I'm looking for, it seems that it's super useful/important.


title: What is a good AUC for a precision-recall curve?

text: Because I have a very imbalanced dataset (9% positive outcomes), I decided a precision-recall curve was more appropriate than an ROC curve. I obtained the analogous summary measure of area under the P-R curve (.49, if you're interested) but am unsure of how to interpret it. I've heard that .8 or above is what a good AUC for ROC is, but would the general cutoffs be the same for the AUC for a precision-recall curve?


title: How do you calculate the standard deviation for Roulette?

text: I've seen quoted in a number of places that the standard deviation for a 1-number bet on a 38 number Roulette wheel (0, 00, 1, 2, .., 36) is **5.76**. I can't seem to find a single source that shows how this calculation is made. I am a student, this is for further understanding of a homework problem, so please assume I'm clueless and be as detailed as possible.


title: How to correctly apply a linear trendline equation

text: I'm a programmer, and unfortunately not a good statistician. I am trying to apply a linear equation (`y=mx+b`) to display a trendline on a graph. Here's a sample of what I'm able to produce: ![enter image description here](http://i.stack.imgur.com/qKX8C.png) Where the Y-axis is a percent relative to the goal (not in currency, but in a percent), and the X-axis is the date (represented discreetly in this image). I've tweaked the graph to show the trendline as an extension from the last day (today) as to visually combine 2 lines into 1 (as opposed to 2 lines running together). From what I understand, we can calculate a trendline using two or more data points. And that's where my question comes in. I'm using the entire current dataset (July 2-24) to calculate my trendline. So in the image, that would be all data points on the blue line. For example, the first point (7/2/13) is represented as `{1,0.28}`; the second as `{2,0.01}`; and so forth. **Am I correctly applying the trendline equation by using that dataset?** As an additional part to the question: should I be using another form of trendline equation, like the polynomial? And, finally, could my linear equation be derived by using two data points from the same time last year? For example, July 26-27's forecast is solved for by using the same dates' historical data from the year prior.


title: Mixed effect models: Why are variance components incongruent with predictive power in presence of hidden interactions?

text: I have noticed that if there are interactions between hidden variables not in the model, then the variance estimates are inflated greatly compared to the predictive power of the model itself, and I'm trying to understand why this is. Here is an illustrative example:               library(lme4)     library(lmerTest)          x1 <- sample(c(0,1),1000,replace=T)     x2 <- sample(c(0,1),1000,replace=T)     y <- x1*x2 + rnorm(1000,sd=0.1)      Now I execute lmer               > lmer(y~0+(1|x1))     Linear mixed model fit by maximum likelihood  ['merModLmerTest']     Formula: y ~ 0 + (1 | x1)           AIC       BIC    logLik  deviance  df.resid       825.2395  835.0550 -410.6198  821.2395       998      Random effects:      Groups   Name        Std.Dev.      x1       (Intercept) 0.3731        Residual             0.3626       Number of obs: 1000, groups:  x1, 2     No fixed effect coefficients       So the model says that x1 explains roughly 50% of the variance components               > my_model <- lmer(y~0+(1|x1))     > var_expl <- as.data.frame(VarCorr(my_model))     > var_expl <- (var_expl[var_expl$grp == 'x1',4])/(var_expl[var_expl$grp == 'x1',4]+sum(var_expl[var_expl$grp != 'x1',4]))     > var_expl     [1] 0.5143418      However, the model built using x1 has only ~35% linear variance explained in reality               > cor(predict(my_model,as.data.frame(x1)),y)^2     [1] 0.3464968      So if I were to interpret this as an experiment, I would say that there is a hidden confound x2 that affects the results of x1 directly... but I know that x1 and x2 are uncorrelated. Again, trying to interpret this practically - there is a hidden factor that happens randomly if I do the experiment but I have not noticed it. (Note, adding a main effect to x1 and an interaction term gives very comparable results) Does this mean that interpreting the variance coefficients in a case with hidden confounds does not tell you anything at all about the predictive power of x1 in this case as a variable? I would like to know if there's a way to correct for such hidden factors and when it is okay to interpret the variance components as percent variance explained (assuming that I only have one variable to look at, but potential confounders that I may not be able to know explicitly)


title: Analysis method for count data with unequal sample sizes and categorical predictors?

text: I have a problem I've been going over and over for months to find the right statistical analysis method. I'm planning to execute the analysis in R, so any mention of appropriate packages is also appreciated. I appreciate any advice anyone can offer. Here's an analogy to my data. Consider that I have 5 colleges. My hypothesis is that one particular college has awarded more degrees than the other colleges. So I want to test whether the number of degrees awarded by the different colleges differs significantly. My data includes 3 data sets as follows: 1)The total number of degrees earned by students from each college, since the opening of each respective college. This is not a random sample, but the actual count. All colleges have different opening dates, so they differ in age. I do have the age of each college if it turns out to be relevant. This data includes second, third, and fourth degrees earned by individual students in the count. So the data includes two columns (college name, and number of degrees awarded), and 5 rows (one for each college). 2)The number of degrees per student per school. I have the list of about 5000 students in the total data set, and the number of degrees each student obtained at each school. Some students attended more than one of the schools, and many earned degrees from more than one school. Every degree ever awarded by these five colleges is represented in this data set as well. This data set does not include students who attended, but never earned a degree at these five colleges. This data set has six columns (one for student ID, and five for college name), and about 5000 rows (one for each student). The totals for each college column are equal to the values in the first data set. 3)The total number of students attending each school over the 30 year period, including those who did not earn a degree. The total number of students, of course, varies at each school. So five rows (one for each college) and two columns (college name, and total students attended in the history of the college). So I have count data with unequal sample sizes. My categorical predictor needs to be the college, as that's what I'm interested in addressing. I've considered different types of regression models, but I'm not sure if these would be most appropriate. For a regression model I guess my response would simply be the number of degrees for each school. Or maybe the proportion of degrees with respect to the total students per school. Would this be sufficient to correct for the unequal sample sizes? Any thoughts on which regression model would be most appropriate? Binomial doesn't work, because I need to take into account students with more than one degree, so it's not a binomial response. I'm not sure what type of distribution should be assumed, as picking one to test against seems very arbitrary to me. Would some type of contingency table be more appropriate? Thanks in advance for any suggestions you can offer.


title: Detecting An Image is Bright or Dark?

text: I use Java and I have some PNG images. I want to put a caption on my images. However I want that: if an image is bright I want to put a black text on image and if an image is dark I want to put a white text on it. How can I detect an image is bright or dark. This should be a fast algorithm and speed is a bit more important than accuracy.


title: Importance of optimizing the correct loss function

text: I want to understand the importance of optimizing the correct loss function. Say that I am building a linear regression model $p$ for predicting some values $y_1,\ldots,y_n$. I choose to fit my linear model such that it minimizes the mean squared errors. Now I send my model off to a statistical prediction competition, where they instead of using MSE as an error metric uses mean absolute error (so no squares). How will this affect my models predictive power? In general, what can be said about the importance of optimizing the "correct" loss function? **Edit:** If one should put the question in a more specific context, then that context would be predictive modelling competitions like those featured on Kaggle.com. I want to understand the importance of choosing models and loss functions which correspond to the evaluation metric for the competition. One reason for this is this comment by the winner of a Kaggle competition.


title: The between estimator in panel data

text: When you model the between estimator of you panel data, you regress the averages of the explanatory variables of the subjects against the averages of the outcome variables of the subjects. But in this regression model, do you have to include an intercept? Out of our textbook: > The between estimator exploits the cross-sectional dimension (differences > between units) of the data by regressing the individual averages of y on the > individual averages of x and a constant using OLS I have played around with the example data of the book of Gujarati: Basic in econometrics, chapter 16. You can find the data here: http://shazam.econ.ubc.ca/student/gujarati/table15.4 This is the plot I made (the colors represent 3 different companies): On the plot you can see the Pooled OLS regression data and the fixed effects (aka within estimator) regression data. How would the between estimator look like? Does it have one intercept or three? And what if all the data points of the companies lie exactly above each other, so they have the same average x? ![Fixed effects plus pooled OLS](http://i.stack.imgur.com/gFxyN.jpg)


title: Regression inconsistent results?

text: I have a question regarding the findings in an article which I don't fully grasp. The authors examined the relationships between variables measured at different time points. They found that a variable V1 predicted positively a variable V2; V2 predicted positively V3; However V1 predicted negatively V3. Is this possible, that although V1 positively predicts V2, it predicts negatively V3? V1 and V2 were measures of before/after an event, collected through the same scale. V1 is time 1 variable, V2 and V3, time 2. Any insights much regarded!


title: Let f(x) be some PDF, and F(x) be its CDF. Shouldn't F(x)=.5 give us the expected value of f(x)?

text: I was playing around in R and have gotten myself very confused about the relation between probability distributions, their expected values, and their cumulative distribution functions. Say we're interested in the expected value of the gamma distribution. Here's a hacky way to calculate it in R:               library("pracma")     m=2     s=1     grid = seq(0,50,.01)     y = dgamma(grid, m, s)     (expected = trapz(grid,grid*y))     [1] 2      Intuitively, the expected value of the distribution is something like its 'center of mass': a point where half the weight of the distribution is on one side, and half is on the other side. If this reasoning is right (I thought to myself), then this point should be identical to the point where the cumulative distribution function is equal to .50. Why? Well, because at the X value where the CDF is .5, a random variable has 50% change of being less than X and a 50% chance of being greater than X-- that is, X is the point where half the weight of the distribution lies on either side. _(EDIT: If I had to pinpoint where my misconception occurred, I would say it was here. Center of mass isn't just about proportion of weight, as the answers below have reminded me.)_ I decided to follow this hunch:               (halfway  = qgamma(.5, m, s))     [1] 1.678347      Nope, my hunch was wrong. The numerical difference seems to be too systematic to be the result of using numerical integration (seems to happen on any asymmetric distribution, to varying degrees, try for yourself). How am I thinking about this incorrectly?


title: Analysis of the Residuals vs Fitted

text: I have a model for which I gathered 10 observations from each person, a total of 25 people, then 250 observations. Well, this is part of my summary of the model,               > summary(m)          Call:     lm(formula = fmla, data = mydata)          Residuals:         Min      1Q  Median      3Q     Max      -9.3311 -3.8480 -0.3134  3.3273 13.4413           Residual standard error: 5.246 on 216 degrees of freedom     Multiple R-squared:  0.7702,    Adjusted R-squared:  0.7351      F-statistic: 21.94 on 33 and 216 DF,  p-value: < 2.2e-16      Looking at this results may seem that the model is pretty significant, but when I plot the Residuals vs Fitted I get this plot, ![enter image description here](http://i.stack.imgur.com/JyyIc.png) Which suggests me that there is a pattern in the plot. I figure this is about the 10 observations for each person. Can anyone help me analyse this plot?


title: Equivalent of paired $t$-test for dispersion

text: I'd like to compare, in a within subject design, the dispersion of two conditions. At first I thought about taking the standard deviation by individual and condition, and compute a paired $t$-test on these values comparing conditions. However I'm not sure about the validity of this method, so I'd like to know which measure you would take at the individual level, and which test you would run on it at the group level.


title: expected values determined with model parameters estimated from a nlme analysis

text: I'm kind new into nonlinear mixed model theory and I've seen that you cannot determine expected values of your response variable by simply inserting the estimated parameters into your model equation, because the model is nonlinear in the random components. Why is that? What is the theory behind that? And, what should I do when I want to insert the parameters back to the model? For example, let's say I have the Brody growth curve:               Wit = A - (A-Wo)e^-kti      where, W is weight of individual i at time t, A, Wo, and k are parameters. Once I run a nlme analysis (e.g. with nlme/lme4 in R), what should I consider to put the estimated parameters back in the model so I can get predicted values of Wit? Many thanks.


title: What's the typical range of possible values for the shrinkage parameter in penalized regression?

text: In lasso or ridge regression, one has to specify a shrinkage parameter, often called by $\lambda$ or $\alpha$. This value is often chosen via cross validation by checking a bunch of different values on training data and seeing which yields the best e.g. $R^2$ on test data. What is the range of values one should check? Is it $(0,1)$?


title: Considering non-i.i.d. covariates in random forests

text: Random forests are theoretically funded on the assumption that the data are i.i.d. realizations from a multivariate random vector $(X_1, \ldots, X_p, Y)$. Does it make sense to use random forests (for prediction and for feature selection based on variable importance) when the predictive variables are not i.i.d. from a multivariate vector $(X_1, \ldots, X_p)$ ? For instance when there is a "trend" in some covariates or when covariates are controlled during the experiments (a fixed design) ? Obviously it makes no sense to consider a variable importance as an estimate of the population variable importance in such a situation. But can we consider that nevertheless, the variable importances are relevant ?


title: How to determinate standard deviation above and below mean classes programatically

text: I have an array of float from wich I calculate the mean and the standard deviation. The user can specify a number of class varying between 2 and 8. The min and max values of each class is calculated by a number of standard deviation above and below the mean. Ex:               nbr_classes = 4;     mean = 10;     min_value = 1     max_value = 10     std_deviation = 2          class_1_min = min_value     class_1_max = mean - 2*std_deviation          class_2_min = class_1_max     class_2_max = mean - std_deviation          class_3_min = class_2_max     class_3_max = mean + std_deviation          class_4_min = class_3_max     class_4_max = max_value      I would like to create a function that would do that programatically, knowing that number of classes can be odd