title: How to test whether the probability of a given colour of skittle is equal to that of another colour?

text: From a family pack skittles, some of the colours of different skittles are counted:   * 41 green,    * 48 purple   * 100 red    * 63 yellow.  Now I want to test if the probabilty of green is equal to the probability of purple; H0: P(green) = P(purple).   * How do I test this is R?    * Do I use a chi square test?


title: How can I check if my time series data is zero mean, stationary and independent identically distributed?

text: I have time series data which is international monthly tourist arrival to Malaysia (N=264). My objective is to forecast tourist arrival for 6 months ahead. After analyzing my data pattern, I found that using Box-Jenkins method gives me the best forecasting model. The question is that, at such preliminary stage, when using Box-Jenkins methodology time series data must follow certain criteria: mean is zero, variance is constant and the data is stationary. I used graphical methods to check if my time series data is stationary, has zero mean and constant variance. However, how to check whether my time series data is independent and identically distributed (iid)? What method should i used?


title: Do the null hypotheses used for benjamini-hochberg have to be identical?

text: There are some cases where I would like to correct for multiple comparisons where the null hypotheses vary. For example, if I am performing a chi-squared test or a fisher's exact test on multiple contingency tables of different sizes, the null hypotheses vary from table to table. Can I still apply the Benjamini-Hochberg procedure for multiple comparisons using the p-values from these tests? If not, what should be done in this situation?


title: What is the difference between linearly dependent and linearly correlated?

text: Please explain what is the difference between if two variables are _linearly dependent_ or _linearly correlated_. I looked up the wikipedia article but didn't get a proper example. Please explain it with example.


title: VAR model and Granger causality issues

text: Is it necessary to create a stable vector autoregressive (VAR) model satisfying all the necessary checks to create it and then and only can we say that Granger causality holds true?


title: Linear model with constraints

text: I'm quite new to R and I have a following problem: I have a simple 2-factor linear model:               # factor1 has 8 categorical values, factor2 has 6 categories     Rate ~ factor1 + factor2      model1 <- lm(Rate~factor1+factor2, data=myData)      And want to put constraints SUM of factor1 coefficients = 0, the same for factor2. None of the manuals gives any clue how to do this. I found a link to similar problem here but it is different and I couldnt figure out how to modify it...


title: How should I objectively test my program results?

text: I have developed two differing methods in MATLAB which aim to analyse a pop song and then automatically create a 30 second audio thumbnail (a preview clip) containing part of the chorus section. Both methods have varying results:   1. The first method can create a thumbnail for each track, managing to find a chorus section in 40 out of 50 tested songs   2. The second method only managed to work on 30 out of the 50 songs, and it found the chorus section 21 times out the 30. Obviously I know which method is superior, but I need to describe and explain the results in a report which requires the demonstration of proper statistical testing. Other academic papers have previously used an f-test to do this, but because their methods are vastly superior, their aims are usually involve the detection of chorus onset times with 100% accuracy. My aim is more relaxed as I am just looking for the generated thumbnails to contain any part of the chorus, regardless of onset. Can anyone suggest some objective tests that I could possibly explore with regards to my project? This is my first time conducting an investigation like this so my experience/knowledge is incredibly low. Thank you! **edit:** Can I possibly conduct an f-test? Recall = Number of thumbnails produced / Total number of tested songs Precision = Number of choruses detected / Total number of thumbnails Would this work or am I completely wrong?


title: Difference between kappa and delta in mixed-effects model

text: I have a following stochastic model describing evolution of a process ($Y$) in space and time. Ds and Dt are domain in space (2D with $x$ and $y$ axes) and time (1D with $t$ axis). This model is usually known as mixed-effects model or components-of-variation models ![enter image description here](http://i.stack.imgur.com/xq9XD.jpg) I am currently developing Y as follow:               %# Time parameters     T=1:1:20; % input     nT=numel(T);          %# Grid and model parameters     nRow=100;     nCol=100;               [Grid.Nx,Grid.Ny,Grid.Nt] = meshgrid(1:1:nCol,1:1:nRow,T);          xPower=0.1;     tPower=1;     noisePower=1;     detConstant=1;          deterministic_mu = detConstant.*(((Grid.Nt).^tPower)./((Grid.Nx).^xPower));          beta_s = randn(nRow,nCol); % mean-zero random effect representing location specific variability common to all times          gammaTemp = randn(nT,1);          for t = 1:nT         gamma_t(:,:,t) = repmat(gammaTemp(t),nRow,nCol); % mean-zero random effect representing time specific variability common to all locations     end          var=0.1;% noise has variance = 0.1     for t=1:nT         kappa_st(:,:,t) = sqrt(var)*randn(nRow,nCol);     end          for t=1:nT         Y(:,:,t) = deterministic_mu(:,:,t) + beta_s + gamma_t(:,:,t) + kappa_st(:,:,t);     end      Can someone help explain, through some illustration using Matlab, if I am correctly producing $Y$? Also, how to produce delta in the expression for $Y$?


title: Regarding Transition Probability Matrix

text: I am calculating transition probability matrix (TPM) using `hmmestimate` command of MATALB. The simple vector `a=[ 1 1 1 1 3 3 3 2 1 3]` and for this the TPM is               hmmestimate(a,a)          ans =     0.6000         0    0.4000     1.0000         0         0          0    0.3333    0.6667      I believe it is calculated as : (1,3) element:               (# of 1->3 transitions)/(total # of 1) = 2/5 = 0.4      Same rule follows for all the elements. Then why not for element (3,2) and (3,3). For example: (3,2) element:               (# of 3->2 transitions)/(total # of 3) = 1/4 = 0.25       but the transition matrix calculated by MATLAB entry 0.33. For `(3,3)`, according to my calculation the entry should have been `(2/4)=0.5` but it is `0.6667` What am I missing?


title: Probability of Largest Sample Observation

text: > Let $X_1, X_2, X_3$ be three random variable following a normal distribution > $N(6,4)$. What is the probability that the largest observation exceeds 8? > Hint: $Y = \max(X_1, X_2, X_3)$. Here is what I tried: $P(Y\leq y) = P(X_i\leq y)$ for $i=1,2,3$. Then $$ P(Y\leq y) = P(X_1 \leq y)*P(X_2 \leq y) * P(X_3 \leq y). $$ To get $P(X_ i \leq y)$, I use the $P(X > 8) = 1 - P(X < 8) = 1 - 0.841 = 0.159$. Then I plug in $0.159$ to the previous equation: $P(X_1 \leq y) * P(X_2 \leq y) *P(X_3 \leq y)$, so that $(0.159)(0.159)(0.159) = 0.0040$. However, my teacher said this is wrong.


title: Are the marginal distributions of a multivariate distribution the corresponding univariate distributions?

text: Are the marginal distributions of a multivariate distribution necessarily the corresponding univariate distributions? For example:   * Every marginal distribution of a multivariate normal distribution must be a univariate normal distribution.   * So it is for the multinomial distribution (which is the multivariate generalization of binomial distribution).   * So it is for Dirichlet distribution (which is the multivariate generalization of beta distribution)   * Not sure for multivariate t distribution. Those are all the multivariate distributions I have seen so far. I am not sure for other multivariate distributions I haven't seen either. So I wonder if it is true for all the multivariate distributions? This may involve a related question:   How is a univariate distribution generalized to multivariate generally?


title: Can Dirichlet and multinomial distributions be defined from their univariate distributions this way?

text: For Dirichlet distribution:   * $X := [x_1, \cdots, x_K] \in [0,1]^K$, $x_i \sim {\rm beta}(\alpha_i, \beta_i)$ and $\sum x_i = 1$. Can we say the distribution of $X$ is a Dirichlet distribution?    * If no, will adding $\sum_{i=1}^K \beta_i = (K-1) \sum_{i=1}^K \alpha_i$ work?    * Or will replace $\sum x_i = 1$ with $\sum_{i=1}^K \beta_i = (K-1) \sum_{i=1}^K \alpha_i$ work? For multinomial distribution:   * $X := [x_1, \cdots, x_K] \in \{0,1, \dots, n\}^K$, $x_i \sim {\rm binomial}(p_i)$ and $\sum x_i = n$ for some constant natural number $n$. Can we say the distribution of $X$ is a multinomial distribution?    * If no, will adding $\sum_{i=1}^K p_i = 1$ work?   * Or will replace $\sum x_i = n$ with $\sum_{i=1}^K p_i = 1$ work?


title: Using GLM with transformed data

text: I have a dataset with outcomes from a Vasicek distribution (see this pdf) and some covariates. Re-expressing the Vasicek's pdf into the exponential family form requires me to transform my data, i.e. instead of $y.\theta + c(\theta) + d(y)$, I have $a(y).\theta + c(\theta) + d(y)$, where $a(y) = \Phi^{-1}(y)$ is the inverse normal cdf function. Will transforming my data by $\Phi^{-1}(y)$ and then fitting it to a GLM violate any assumptions?


title: Am I getting something wrong? A question about a paper on the F-test

text: First, I'm sorry for the long post, but I needed a second opinion from you, the experts, about this problem! I was reading a paper by Gallagher (2006), where he puts an example on "how one may correctly perform an F-test and how some softwares perform poorly on this task" (which could have been a pretty interesting papers for students, but I think I found an error, which maybe it'll make it even more interesting!). He poses the following example: > _Consider the following two artificial data sets:_ > >> **Sample 1** : 1,2,3,3.3 $\Rightarrow$ $s_1^2=1.08917$ >> >> **Sample 2** : 1,2,3 $\Rightarrow$ $s_2^2=1$ > > (where $s_1^2$ and $s_2^2$ are the estimated sample variances). This is ok (thanks to @Procastinator!). However, here's what bothers me: > _(...) it may be erroneously assumed that $F=1.08917$ lies in the upper tail > of the F distribution on 3 and 2 degrees of freedom, as $F >1$. Clearly, > from (1)_ (where the equation shows that $P(F_{3,2}>1.08917)=0.5114$) _this > is not the case (...)"_ Correct me if I'm wrong, but the mode for the F-distribution is $$\hat{f}=\frac{\nu_2(\nu_1-2)}{\nu_1(\nu_2+2)}$$ which for $\nu_1=3$ and $\nu_2=2$, gives $\hat{f}=1/6=0.17$. Because $F>\hat{f}$, then **actually the value resides in the upper tail of the distribution**. There's nothing strange with having a probability greather than $0.5$ because the PDF of the F-distribution for $3$ and $2$ degrees of freedom is actually pretty skewed. **Did I miss something here?** [1]: The F test for Comparing Two Normal Variances: Correct and Incorrect Calculation of the Two-Sided p-value?


title: Improving transformation of dependent variable and robust regression

text: In a multiple regression with 16k cases 2 IV (non-normally distributed) and one dependent variable that is also not normally distributed. DV see below: ![enter image description here](http://i.stack.imgur.com/zQ524.png)![enter image description here](http://i.stack.imgur.com/tOYbl.png) I've tried three ways of transforming the DV (sqrt/ln/log) to normal but the K-S statistic reports it clearly as non-normal . ![enter image description here](http://i.stack.imgur.com/TffAZ.png) The best I could come up with was the ln or log transformation, which clearly shows a trend in the detrended qq plot. Any ideas how to "fine-tune" the ln transformation to remove that trend? ![enter image description here](http://i.stack.imgur.com/UIPzZ.png) ![enter image description here](http://i.stack.imgur.com/yi5YS.png) Now I've performed the same ln-transformation which on both the IV which were also not normally distributed. The result of the transformation is shown below: ![enter image description here](http://i.stack.imgur.com/5vnvU.png) ![enter image description here](http://i.stack.imgur.com/rCgw5.png) ![enter image description here](http://i.stack.imgur.com/ZdSA2.png) ![enter image description here](http://i.stack.imgur.com/MoraJ.png) My approach was instead of trying harder to better transform the IV and DV, to bootstrap the regression in spss and and use the ln_RT_vol_in as the DV and ln_AT_vol_in and ln_FF_vol_in as the predictors. The residuals look fine at least the histogram. The regression without bootstrapping:               REGRESSION      /MISSING LISTWISE      /STATISTICS COEFF OUTS CI(95) R ANOVA COLLIN TOL      /CRITERIA=PIN(.05) POUT(.10)      /NOORIGIN      /DEPENDENT LN_RT_vol_in      /METHOD=ENTER LN_AT_vol_in LN_AT_bin_in_deg      /PARTIALPLOT ALL      /SCATTERPLOT=(*ZRESID ,*ZPRED)      /RESIDUALS DURBIN HISTOGRAM(ZRESID) NORMPROB(ZRESID).      ![enter image description here](http://i.stack.imgur.com/uAu05.png) ![enter image description here](http://i.stack.imgur.com/9d28b.png) And with bootstrapping. I have the same B values but slightly bigger confidence intervals for the bootstrapped version and a slightly less see below.               BOOTSTRAP      /SAMPLING METHOD=SIMPLE      /VARIABLES TARGET=LN_RT_vol_in INPUT=  LN_AT_vol_in LN_AT_bin_in_deg      /CRITERIA CILEVEL=95 CITYPE=PERCENTILE  NSAMPLES=1000      /MISSING USERMISSING=EXCLUDE.      ![enter image description here](http://i.stack.imgur.com/mIQz1.png) What do you think of this approach? Do you have any comments on how to improve the transformation in order to maybe get rid of the trend in the de-trended versions of the qqplots? Another problem is that the bootstrapping doesnt seem to fit the 1000 samples in my memory instead only 49 see below. What can I do about it apart from increasing my mem?


title: Significance of interactions

text: It was suggested to me recently that the significance of an interaction term in a glm has to be higher than a main effect. For example, p<0.05 is commonly thought of as significant for a main effect but a two-way interaction has to be higher (p<0.025 or something) and a three-way even higher. However, I can't find any literature on this so I am unsure. Are they getting confused with multiple comparisons I wonder? Or am I confused and the two things are linked?


title: How to interpret result for MultiModelByRegression in RapidMiner?

text: I couldn't find any information in the documentation of rapidminer. I have a data set with the following attributes: `a,b,c,d,e`. The types are: `numerical`, `binomial`, `binomial`, `binomial`, `binomial`. Binomial values are given as `{true, false}`. The last one is the label I want to be able to predict. So the value to predict is a true/false decision. I understand that I have to use logistic regression for that. My module chain looks like this.               Read CSV -> Set Role -> Nominal To Binary -> Classification by Regression      `Set Role`: I can only set one attribute as the predictor attribute. How can I set all other input attributes as predictor attribute too? `Nominal to Binary`: The binomial values are given as Strings "true" and "false". That's why I do the conversion here. The output for MultiModelByRegression is:               MultiModelByRegression (prediction model for label is_helpful)                16.889 * b ^ 1.000         + 9.553 * c ^ 2.000         + 0.102 * a ^ 1.000         - 71.438           76.078 * b ^ 4.000         + 38.618 * c ^ 1.000         + 0.082 * a ^ 1.000         - 88.701      The `Performance Vector Output` is:                   true false  true true   class precision     pred. false 2706    129 95.45%     pred. true  636 40  5.92%     class recall    80.97%  23.67%        I know how to interpret the above. All of this is done with the training set, which is already labelled with the correct classes. How do I apply the test set at this point? I suppose I need the test set to somehow evaluate the results of my classifier, right? Anyway, I am really confused here and I would appreciate any kind of help.


title: Standardized mean decomposition (Oaxaca-Blinder decomposition)

text: How does standardized mean decomposition differ from the simpler dummy variable regression? How does the average in mean outcomes or interpretation of results attributable to a particular treatment or dummy differ when standard decomposition is used, as opposed to dummy variable regression?


title: Distinguishing missing at random (MAR) from missing completely at random (MCAR)

text: I've had these two explained multiple times. They continue to cook my brain. Missing Not at Random makes sense to be, and Missing Completely at Random makes sense...it's the Missing at Random that doesn't as much. What gives rise to data that would be MAR but not MCAR?


title: Regression analysis or another test for analyzing binomial data over 9 months?

text: I have 3 binomial data questions asked over 9 months to over 1000 diabetic patients. I would like to test the hypothesis that each of these 3 (quality of care) questions asked to the patients increased over 9 months. These patients who went into a hospital for a medical reason other than their diabetes were treated not in full the way diabetic patients should be checked on and maintained. To test the quality of care at the hospital these patients were asked 3 yes,no questions once they left the hospital to see if the hospital workers did what they were supposed to for diabetic patients. What test is best for analyzing this data and finding out if the quality of care increased (these 3 questions were answered yes) over the reported nine months?


title: Expected value and variance of arithmetic mean of random variables

text: The question is > Let $X_1,...,X_n$ be drawn iid from $Beta(0.1,0.5)$. Let $\bar{X} = > \frac{1}{n}\sum^n_{i=1} X_i$. > > a) Derive $\mathbb{E}(\bar{X})$ and $\mathbb{V}(\bar{X})$ I know how to get the expectation and variance of a random variable from a given distribution. But I'm not how to get the expectation and variance of $\bar{X}%$. Not sure where to being actually.


title: What is $\bar\mu_n$?

text: Am I right in thinking that it is the average of the sum of $n$ different populations means? Here it is used in the context that confused me. It's the Chebychev WLLN, apparently. "If $x_i, i = 1, . . ., n$ is a sample of $n$ observations such that $E[x_i] = \mu_i < \infty$ and Var[$x_i] = \sigma_i^2$ such that $\bar\sigma_i^2/n = (1/n^2)\Sigma_i \sigma_i^2 \rightarrow 0$ as $n \rightarrow \infty$ then $plim(\bar x_n - \bar\mu_n$) = 0." Is this saying that each sample of $i$, corresponds to it's own population of $i$ (from above, $E[x_i] = \mu_i$) and as the sample get bigger we have to average over populations? So if I were to draw the random variable 1 from a population of {1,2,3} and the random variable 4 from the population {4,5,6} then the population of my sample 1,4 is {1,2,3,4,5,6}?


title: Understanding ANOVA formula in R and choosing the most appropriate model

text: I am new to R (and by extension statistics), and am having a tough time in making sense of this repeated measures ANOVA formula, please bear with me. I think I understand what the formula implies, but am not too sure on its application.               aov.out = aov(Factor1~Factor2*Factor3 + Error(Factor4), data=Data)      So I looked up this link What is the meaning of operators in regression or anova formulas in R and assumed that the forumla is basically determining the dependence of `Factor1` with all possible interactions of `Factors2,3` and having `Factor4` as an error? Now one of my colleagues gave me a formula to use for my experiment set. Now suppose in an experimental situation, We have measured `Factor1` (lets call it mass gain), `n` times (Factor3) with a certain number of people (Factor4) divided into treated and control groups (Factor2). So does this ANOVA make sense?                  aov.out = aov(MassGain~ExperimentalGroups*Trial + Error(Individual), data=Data)      It is assuming that the Mass Gain depends on which experimental groups he belongs to, and the trials undergone by that person and the error accounts for the inter individual variance. Now considering that the trials took place in the same time frame (not separate days) having that as an independent variable does not make sense, as during each trial we have ensured that the conditions remain the same. I was thinking of modifying the formula to reflect this, so that it either looks something like this: `aov.out = aov(MassGain~ExperimentalGroups + Error(Individual + Trial), data=Data)` which takes into account the error due to individual variances and the variances that may have occurred during the trial. Is additive the right thing to use here? I mean in the sense that the individual's are going to be different from each other, but does the individual variance also interact with the repeated measurements of the Mass Gain? Or should I completely ignore the trial variable, since we are keeping the conditions same when we take our observations , and have something like this? `aov.out = aov(MassGain~ExperimentalGroups + Error(Individual), data=Data)` Do note that I am not a statistician, merely a biologist trying to apply the right technique to look at my data, and hence please do correct me where I am wrong.


title: Randomization testing in Bayesian spatial scan statistics

text: I was reading about Bayesian Spatial Scan statistics paper. I have this confusion about why randomization testing is not necessary in this approach. The paper says that it is by construction. But, I didn't get how it is so?


title: Multivariate normality plots

text: To check multivariate normality, which is the best plot - a chi square qq plot for the entire set of variables or as many qqplots as the number of variables?


title: Conceptual questions on Entropy and estimation

text: Learning Informative Statistics: A Nonparametric Approach  paper presents an approach to parameter estimation by entropy minimization. There are other related works "Minimum-entropy estimation in semi-parametric models" download link ( http://dl.acm.org/citation.cfm?id=1195853). The rationale provided is that minimization of error entropy is equivalent to maximization of likelihood. I am new to this area and find it hard to understand the intuition behind why entropy of the error minimization will yield the parameters. What happens when entropy is minimized?   1. What happens when Shannon Entropy is maximized? Entropy (Shannon's) is the uncertainty = average information or uncertainty (unsure).    2. And what happens when entropy is minimized and    3. What is the meaning of minimizing entropy of error?


title: Confidence interval scoring

text: Question migrated form http://stackoverflow.com/questions/10019576/confidence- interval-scoring-with-programming-languages#comment12812602_10019576 I used a CI scoring algorithm from http://evanmiller.org/how-not-to-sort-by- average-rating.html, The data is most liked programmings and most disliked programming languages of each person. According to the scores it should seem that `Python` should have a relative higher score, but according the the score data it is in the middle. Am I doing something wrong? `Joel Cornett` mentioned `You need to use Student's T curve, not the normal distribution for smaller sample sizes.` Is there a way to fix this? Here is my gist https://gist.github.com/2305834, which contains the `data`, `scores` and the `script` that generates the scores.


title: about test for multivariate normality with kolmogorov-smirnov test and $\chi^2$ test in R program

text: how can I test for multivariate normality with kolmogorov-smirnov test and $\chi^2$ test in R program?


title: Incorporating utility function with error function in a neural network model

text: I want to weight errors differently based on domain knowledge by introducing a cost function. I am trying to solve a regression problem with neural networks and the delta error (`predicted output - actual output`) is multiplied by the cost weight depending on custom defined cost function. For a certain error range, the cost weight is 100 and this is leading to really high parameter values being learned. Should I do some sort of normalization over the actual output values, delta error before propagating them back for gradient computation.


title: Interpretation of $R^2$ in simple regression

text: I have run a linear regression with five independent variables and two dependent variable I am using EPS and ROE as my dependent variable and Board meeting, Board size, leadership, Independend Non-Executive, Audit committee as my Independent variable and my resulting Model is shown as:               Model Summary 1    (Dv: ROE)          R. Square 0.39     Adjusted R-Square = 0.014     Error of the Estimate = 19.4807565     a. Predictors: (Constant), Board meeting, Board size, leadership,           Independent Non-Executive, Audit committee.           Model Summary 2    (Dv: EPS)     R   .308a                R Square = .095     Adjusted R Square = .072     Std. Error of the Estimate = 35.308797     a. Predictors: (Constant), Board meeting, Board size, leadership,           Independent Non-Executive, Audit committee       I got confused regarding how to explain the correlation of my independent and dependent variables. I'm also wondering if my results are any good. Thanks.


title: How well the experiment would do in excluding alternative hypothesis?

text: The experiment has $N$ observations $N_i$ which are Poisson Distributed random variables (i.e. we have histogram with $N$ bins. The width of bins is not fixed so I can adjust it ) . I have two hypothesis.The null hypothesis predicts $\lambda_i^{null}$ values of corresponding parameters of Poisson Distribution, alternative predicts $\lambda_i^{alt}$. The more data we collect, the better confidence level of excluding of alternative hypothesis we get. But how to calculate this confidence level as a function of size of data we collect? In other words I would like to estimate how many data should I collect to exclude alternative hypothesis at specified confidence level. Firstly, I thought about using chi-squared test and calculate the expected value of chi-square i.e. $\left<\chi^2\right>$. But when I tried different width of bins to obtain the optimal binning, I discovered that the less bins we use the better result we get. So I think it's not the optimal test for this task because we don't fully use information from predictions when we have just a few bins. I would like to know which tests are more appropriate for such problems. I would also be grateful for any kind of advice regarding the literature on this subject. P.S. I've asked the same question at physics.stackexchange, but with some physical details.


title: Regarding the size of training data for building classifier

text: When we build a classifier, like SVM or Naive Bayesian, are there any generic rules or theoretical derivations on the size of training data set? For example, to train a SVM-based classifier, what should be the minimum size of training data in terms of feature space and some target performance metrics, such as precision and recall?


title: Repeated-measures ANCOVA in R?

text: experts, I would like to ask for advice regarding the analysis of a dataset I am currently working on. In the experiment subjects were tested in a reaction time task (RT as dependent variable). In random order every subject was tested in two treatment conditions (factor 'A': placebo vs. drug). Also for subject in a third session the level of a blood value was measured (a numeric covariate 'B'). This covariate 'B' is stable within a subject and believed to be a predictor for the effect of factor 'A'. Moreover subjects were grouped according to their genotype (factor 'C'). Also I coded the first or second testing in factor 'D' in order to model potential learning effects. Please find below some sample code.               # some example data     subject <- c(1,1,2,2,3,3,4,4)     A <- rep(c('placebo', 'drug'), 4) # factor A     B <- c(1,1,4,4,8,8,12,12) # covariate B     C <- c('x','x','y','y','x','x','y','y')     D <- c('1st','2nd','1st','2nd','2nd','1st','2nd','1st')     RT <- c(2,12,1,16,2,26,3,39)          data <- as.data.frame(cbind(subject, A, B, C, D, RT))     data$B <- as.numeric(as.character(data$B))     data$RT <- as.numeric(as.character(data$RT))      I would know how to estimate the overall effect of the covariate 'B' using lm() and how to estimate the overall effect of factor 'A' using aov().               # Model effect of covariate B using lm()     data.lm <- lm(RT~B, data=data)     summary(data.lm)          # Model effect of factor A using aov()     data.aov <- aov(RT~A+Error(subject/A),data=data)     summary(data.aov)      From the comments I understood that the correct ANCOVA model to test effects of A, B and their interaction would be as followed:               # ANCOVA for effects of A, B and their interaction     data.aov2 <- aov(RT~ B * A + Error(subject/A), data = data)     summary(data.aov2)      However if I would like to look at B, C and their interaction the output from summary() does not provide significances anymore which makes me assume the model is not correct:               # ANOVA for effects B, C and their interaction     data.aov3 <- aov(RT ~ B * C + Error(subject), data=data)     summary(data.aov3)      Also of course I would like to know whether it is possible to estimate the most comprehensive model including all factors (A, B, C, D) and their interaction (?). Probably the dataset is too small to estimate this. But how can I know whether my data is sufficient to allow for the estimation of such a complicated model?               # ANOVA for effects B, C and their interaction     data.aov4 <- aov(RT ~ A * B * C * D + Error(subject/A), data=data)     summary(data.aov4)      Also I would like to know whether the same models can be implemented in the same way in case there are repeated-measures of the dependent variable 'RT'? E.g.:               data2 <- rbind(data, data, data, data, data)     data2$RT <- rnorm(nrow(data2))          # ANOVA for effects B, C and their interaction     data.aov5 <- aov(RT ~ A * B * C * D + Error(subject/A), data=data2)     summary(data.aov5)      Any help regarding this (including hints for an implementation in R as well as literature) is highly appreciated. Many thanks! Jokel


title: Coefficient of variation of the product of a random variable and a constant

text: I have a random variable $X$ and I need to calculate the coefficient of variation ($\text{CV}$) of the product of a constant, $c$, and $X$. In other words $\text{CV}(X \cdot c)$. To do this I can apply the following: \begin{align} \text{CV}(X \cdot c) &= \frac{\text{sd}(X \cdot c)}{E(X \cdot c)}\\[1em] \text{CV}^{2}(X \cdot c) &= \frac{\text{Var}(X \cdot c)}{E^2(X \cdot c)}\\[1em] \text{CV}^{2}(X \cdot c) &= \frac{c^2\text{Var}(X)}{c^2E^2(X)}\\[1em] \text{CV}^{2}(X \cdot c) &= \frac{\text{Var}(X)}{E^2(X)}\\[1em] \text{CV}(X \cdot c) &= \frac{\text{sd}(X)}{E(X)}\\ \end{align} So basically multiplying a random variable by a constant doesn't change the CV. Is my math correct?


title: Confidence with probability values

text: I am constructing a predictive model for a binary outcome, which gives me results in a "probability-prediction" fashion such as "0.3-YES" or "0.4-NO". The model is working perfectly. My question here is what would be the probability that I should trust? is any value below 0.35 or above 0.65 confident enough? I know this varies according to my data, but I am asking about the typically used value to consider my prediction as "confident". For example in hypothesis testing, a 0.05 p value is the most commonly used value for significance of the existence of a difference between populations. What is the typically used values for a probability that ranges between 0-1? Cheers,


title: Urn probability function

text: Suppose I have an urn with an infinite number of balls which can be either red or white. I do not know what the proportion of each colour is, but I _do_ know it's a fixed proportion. After drawing $N$ balls, I have observed $r$ red ones and $w$ white ones. I believe the probability that I will observe a red ball on the next draw is given by Laplace's Law of Succession, $\frac{r+1}{N+2}$. However, how sure should I be of that? That is, before I drew any balls, I believed any proportion other than $0$ or $1$ was the true one. After I drew those $N$ balls, what should be my estimate pdf over the possible values for the proportion of red balls in the urn?


title: Do interactions in mixed designs inflate the main effect of fixed factor?

text: Explain in clear way why in a mixed design the presence of an interaction between a random and a fixed factor inflates the estimate of the main effect of the fixed factor.


title: How to perform a Wilcoxon signed rank test for survival data in R?

text: Say you have survival data like this:               obs <- data.frame(       time = c(floor(runif(100) * 30), floor((runif(100)^2) * 30)),       status = c(rbinom(100, 1, 0.2), rbinom(100, 1, 0.7)),       group = gl(2,100)       )      To perform a standard log rank test, one can use               survdiff(Surv(time, status) ~ group, data = obs, rho = 0)      right? But what about other test? How could you perform a Wilcoxon signed rank test, a Peto test or a Fleming-Harrington test? R provides a possibility to perform a Wilcoxon test, however I didn't find how to let it take censoring into account. Furthermore the doc states that setting `rho = 1` would make the test a "Peto & Peto modification of the Gehan-Wilcoxon test". But is this the same as the Peto test?


title: Robustly standardize residuals in MM regression

text: Does anyone know how we can robustly standardize the residuals in MM regression? First we perform MM regression and then obtain the residuals: how can we robustly standardize the residuals obtained from MM regression? I have found the method for least median squares (LMS) and least trimmed mean squares (LTS) in which the scale of the errors is estimated using a formula and then the residuals will be divided by that estimated scale. But for MM regression I could not find a formula for estimating the scale of the errors in order to standardize the residuals.


title: Very different scale parameter estimates in Poisson regression

text: **The background** : I'm analysing survival data using a Poisson model. I've splitted the data on 2 time-scales (attained age and calendar year). Attained age is modelled using flexible parametric functions, calendar year is a categorical variable and I've also included the interaction terms between attained age and calendar year. **The problem** : looking at the goodness-of-fit criteria, I see that the deviance divided by the degrees of freedom ( _df_ ) is 0.03, while the Pearson's $\chi^2$ divided by the _df_ is 5.84. The fact that the former is <1 and the latter is >1 implies that the standard errors associated to my parameter estimates are deflated, if I estimate the scale parameter using _Pearson's $\chi^2$ / df_ , or inflated, if I estimate the scale parameter using _Deviance / df_. This obviously leads to completely different results. What could be wrong here? **Edit (04 Jul 2012)** : I think this question can be reformulated also as: with clustered data (i.e. dependent data within each cluster), are the goodness-of-fit criteria for overdispersion (namely Pearson's $\chi^2$ divided by the residual degrees of freedom and Deviance divided by the residual degrees of freedom) still interpretable? (With iid data, >1 indicates _overdispersion_ , while <1 indicates _underdispersion_ ). **Edit (19 Aug 2013)** : After delving into this topic a little bit more, I no longer think that in this particular scenario the use of those goodness-of-fit criteria makes sense. In fact, we're only exploiting the GLM framework to fit a survival model characterized by a non-constant baseline incidence rate (see also: piecewise exponential model). We're not assuming that the outcome has a Poisson distribution. For this reason, the question as reformulated in my last edit (04 Jul 2012) is also wrong, as we don't have clustered data here. See this answer for more information.


title: What is the dimension (or units) of a CDF and PDF?

text: Given a continuous random variable $X$, what are the units of the PDF and CDF of $X$?


title: What are the most commonly used predictive models when dealing with binary data?

text: I know everybody uses logistic regression as the starting point, but I'm curious to know: What are the other commonly used predictive models when data is primarily binary?


title: Determining suitable time lag in Granger causality test

text: In order to determine the right time lag for the Granger causality test, I applied a cross correlation function (CCF) among EU and National legislation and the CCF figure is attached. I would be grateful if you could kindly help me in resolving the following concerns:   * I am wondering what time-lag you recommend based on the attached figure?   * To find the greatest value in CCF figure, should I consider absolute values or just the raw values?   * If we consider the negative values in CCF Figure, then I assume the result of Granger causality test would show the negative causality, correct?   * I am interested in one side of the causality (Effect of EU on National and Not the vice versa), so, should I only consider the bound [-5...0]?


title: Score Variance of supplementary individuals in PCA

text: I have posted an almost similar question last week but I want to get some more of the theoretical background Difference in Variance of factor scores for supplementary and active observations PCA. I am doing PCA on genetic data (number of variables >> number of rows) and all data points take values in {0,1,2}. The aim of doing the PCA is to isolate the effect of different origins (People from Europe have different genetic structure than people from Asia). I have run the PCA with several different packages on two kinds of different data sets. One with little structure ( all the samples had the same origin) and one including samples from two different world regions. I use a subset of the data to find PCs and then apply the loading matrix to supplementary samples. What I observe is the following. In the case of the data set with little structure, as expected the eigenvalues are almost similar and so the proportion of total population variance explained by the first principal component is small. Now if I apply the loadings matrix to supplementary individuals that have been standardized using the mean of the active sample the score variance of these supplementary individuals is much smaller (3 to 5 times smaller). I then did the PCA only on the supplementary individuals and it turns out that the PC1 of the two subsets are poorly correlated (~ $0.07$). When I do the exact same thing on data including samples from two different world regions than the difference in Score Variance between the active and supplementary individuals is much smaller ( $\frac{Var(Score_{passive})}{Var(Score_{active})} = 0.9$) and the PC1 of the two subset are strongly correlated (~ $0.95$). However the the Score Variance of the passive sample is still smaller than the one of the active sample. Now my question: Is this a global result? Is the Score Variance of the supplementary individuals always smaller than the one of the active? I'd really like to understand but I'm having trouble getting the mathematics and Intuition behind that straight.


title: Bayesian Belief Network - directions of arcs between nodes

text: I generated a BBN below based on environmental variables and a response of some organism. ![BBN](http://i.stack.imgur.com/GpZh6.png) My aim here is to see how environmental variables (A-H in a graph above) interact with each other and influence organisms (phenotypic response). The sensitivity of "organism" is as follows: Sensitivity of 'organism' to a finding at another node: Node Variance Percent Variance of ---- Reduction Beliefs   Target 0.00393 100 0.2004898   A 0.003929 100 0.2004673   B 0.001801 45.8 0.0919060   C 0.001537 39.1 0.0783964   D 0.001203 30.6 0.0613732   E 0.001027 26.1 0.0523963   F 0.0009696 24.7 0.0494689   G 0.0009023 23 0.0460363   H 0.0007741 19.7 0.0394967 All good except a direction of arc between A and B. The direction is actually counter-intuitive. But when I reverse the direction, the sensitivity drops. Is it feasible to say B potentially influence "organism" through interaction with A? Based on this model, is it ok to say B influences A? But how can I explain the direction of arc between A and B? My guess is that a probability of A can be best predicted by a combination of E, F and G, but not with B, and information on A helps predict a probability of B. So the direction between A and B is not necessarily mean that A affects B (i.e. B would affect A and also indirectly affect organism). I would appreciate any comments. Thanks!


title: How to verify the reliability of a survey

text: My question is about the **reliability of the survey data** (survey answers). Do you know any statistical methods for that purpose ? Maybe there are methods which use the data of the demographic questions to test the reliability !? _Let's say that a survey is reliable, if I administer it many times and get roughly the same distribution of results each time_. Now, my idea is the following: Instead of repeating the survey many times, to check whether or not it's reliable. I could just split the the result set (according to the demographic data) _n_ times, and then compare the the _n_ results to each other. Can I do this ? What should I consider ? I'm pretty sure there are already methods/approaches that handles this. It would be nice to know this methods/approaches, so that my idea get's a mathematical ground.


title: Chi-squared distribution for dice not returning expected values?

text: I'm trying to verify that my dice are reasonably fair and balanced. I'm not sure I'm doing it right, though, since I'm basically teaching myself stats without references. To achieve this, I've settled on the chi-squared test. Here's what I'm doing: I have a sample size $N$ of rolls. On a dC, where $C$ is the number of sides, the expected frequency for each roll 1 thru $C$ is $\frac{N}{C}$. Then I compute the $\chi^{2}$ value for the distribution from $1$ to $C$. Then, I compute the chi squared distribution (in Excel using `CHIDIST`) with degrees of freedom equal to $C-1$. Here's an example case:   * Take a d4, and assume each face was rolled equally 1500 times.   * Compute the chi-squared value for this, with each face and an expected frequency of 1500   * Compute the chi squared distribution using `chidist` and 3 as the degrees of freedom The problem is, this is returning nonsensical values. For instance, on a d4, if I input 6000 evenly distributed rolls, `chidist` returns `.801252`. No matter how high I make $N$, even if there's an even distribution, I'm not getting below `.801252`. I know I'm doing something wrong here, but I'm not sure what it is. The chi squared distribution value should be falling to within any sort of significance. Any help? * * * Edit: Per request, I'm posting a screenshot of my Excel layout and listing the formulas used. I've hidden the irrelevant rows. ![enter image description here](http://i.stack.imgur.com/XPmpw.png) Formulas used:   * `C2:C5 = B22/4` (Frequency)   * `B22 = SUM(B2:B5)` (Sum)   * (StDev & StdErr are not relevant)   * `B25 = CHITEST(B2:B5, B91:B94)` (Chi squared)   * `B26 = CHIDIST(B25, B1-1)` (Chi distribution) Here's some raw data:                  4  -  6     1| 56    36     2| 45    33     3| 57    33     4| 42    40     5|       29     6|       29


title: Using quantile regression to predict probability of surpassing threshold

text: Consider a continuous response $Y$ and design matrix vector $\mathbf{X}$. These are related through some function $f(X) = Y$. Suppose that I am interested in estimating the probability that $Y \leq 0.1$ conditional on observing $\mathbf{X}$. I want to use quantile regression to do this - **can I confirm that this is a legitimate methodology**? We have quantiles $\tau \in [0,1]$ and after estimating our quantile regression for each $\tau$ we have our quantile estimates $\mathbf{q} := \{\hat{Q}(\tau) : \tau \in \{0.01,0.02,...,0.99\}\}$. I want to select the $\tau$ such that $\hat{Q}(\tau) \approx 0.1$. When I find such a $\hat{Q}(\tau)$ it seems to then follow naturally that $P(Y \leq 0.1) = \tau$. The reason is that my model has estimated the $\tau$-th quantile to be $0.1$, which is point on the x-axis in $Y$'s pdf that I need to find to be able to determine $P(Y \leq 0.1)$. In practice this may not work since an estimated quantile can be lower for higher $\tau$ under some $\mathbf{X}$. Not looking for logistic regression with a discretized response as a solution (since I already know about this).


title: Comparing correlations of two groups

text: I have two groups of correlations and I need to test if there's a significant difference between these two groups. Each group consists of correlations between two different variables. I transformed each Pearson r to z according to Fisher r-to-z transformation and calculated a mean of correlations per group. But I'm not sure what to do next. I read about Pearson-Filon test but its explanation doesn't go beyond 4 different correlations (comparing a correlation between a-b and between c-d) but in my case there are many more correlations in every group (21 correlations in one group (between 7 different variables)) and 65 correlations in another group (between 12 different variables)). The correlations in two groups were calculated on the same sample. I was suggested to do a t-test and compare the mean correlations of two groups but I'm not sure about that. Thank you!


title: Is likelihood ratio test the only way to build hypothesis tests?

text: Usually we can construct likelihood ratio for testing the Null hypothesis and alternative hypothesis: The likelihood ratio test $ P( l(\beta_{1}) / l(\beta_{2}) ) < \alpha $ is the rejecting region for the null hypothesis. Then the inequality would reduce to a formula with sufficient statistics as the variable. My question is : is likelihood ratio test the only way to build up hypothesis test? If not, what would be the alternatives? If there are multiple alternative methods, when would likelihood ratio test be preferred ?


title: When designing experiments why is it necessary to have multiple groups run in parallel?

text: Let's say I wish to test the effectiveness of a new technique of teaching as opposed to an old one. Why is it recommended to have a control group running the experiment in parallel with the group trying out the new technique (what is that group called :)? What is the loss if these experiments are done non- parallely? That is, we try method 1 in year 1 and method 2 in year 2 and so on. Now, the set of students would be different but that's true even if the experiment is run in parallel. I understand that certain confounding variables could come into play and we'd like to minimize/eliminate them as much as possible. But if the costs of running the experiment with two or groups is prohibitively expensive then one may adopt for a one group at a time experiment. So if we strive to keep as many variables of the environment "constant" (to minimize confounding), can we conduct any meaningful statistical analyses on the outcomes of the experiments? (Ex.: t-tests, ANOVA etc.,) What can we say about the validity of that tests of statistical significance in this case or does it not matter?


title: Interpreting CCF correlation in R

text: I am using `ccf` to find a correlation between 2 time series. I am getting a plot that looks like that: ![enter image description here](http://i.stack.imgur.com/S7f65.png) Note that I am mainly interested in correlation for the lag=0. Questions:   1. Do interpret it correctly that there is a cross-correlation for the lag=0, as for this lag the cross-correlation is above the dotted line?   2. How should I interpret the level of cross-correlation in this example, is this significant (as I interpret it right now, there is a small cross-correlation)?   3. How can I extract only `acf` value for lag=0?


title: Retrieving minimum width that contains specified fraction of all values

text: I'd like to find the min/max boundaries of a sliding window of minimum size that contains a certain fraction of the total number of elements in an array or collection of numbers. Example: taking integers to make it easier to explain, say the proportion of elements that we're looking for is 50% on this array: [1,1,2,3,3,3,3,3,3,4,5,6,7,8,9,10] The mystery function would return something like ( 2, 4 ), meaning that the minimum value is 2 and the max value is 4 (let's assume it's inclusive and note that these are VALUES, not indices). That little window contains 8 of the 16 values and it's only 2 units wide, the narrowest such window that contains half the number of values. Note: the quartiles are [1,1,2,3], [3,3,3,3], [3,4,5,6], [7,8,9,10] so they are NOT what I'm looking for. I can code it...may take some time to get it nice. Hopefully somebody has run across it before. Does it have a name? Is there an existing method that somebody knows about? I am interested in Python implementations. UPDATE: Thanks to Glen_b who gave the names 'short-half' and 'shortest interval' I was able to come up with this: Find probability density intervals


title: Maximizing true positives

text: I'm working on a classification model which is to be built on slightly imbalanced data (nothing extreme, about 60-to-40). Now, I thought to use random forest for this and get some discouraging results (which is not related to this question). Something I just realized is that all I really need for my particular case is a good positive predictive value, namely **true positives/(true positives + false positives)**. That is, if it says True, it should be True with very high probability, and I don't care if it identifies positively only a small fraction of the events (ideally, there would be some threshold variable determining this portion). I think this problem set-up could have some name, but I am not aware of it. (I am also aware of methods used in, i.e. spam/fraud detection, but those are usually built to minimize false negatives, AFAIK) So, is there a family of classification methods designed in particular for such task? What is considered the best out-of-the-box method for this?


title: Covariance matrix equality

text: The (unbiased) sample covariance matrix $$\mathbf{S}=\dfrac{1}{n-1}\sum_{j=1}^{n}(\mathbf{X}_{j}-\bar{\mathbf{X}})(\mathbf{X}_{j}-\bar{\mathbf{X}})^{T}$$ can be rewritten as $$\mathbf{S}=\dfrac{1}{n-1}\mathbf{X}^{T}\mathbf{X}-\dfrac{1}{n(n-1)}\mathbf{X}^{T}\mathbf{1}\mathbf{1}^{T}\mathbf{X}$$ where $$\mathbf{1}=\left(\begin{array}{c} 1\\ \vdots\\ 1 \end{array}\right)_{(n\times1)}.$$ One (tedious) way of proving this is to expand out the left hand side and the right hand side of the equality, and showing that the entries of the matrices match. I've done this successfully. My question: is there a neater / more concise way to prove such an equality?


title: Difference between significance and contribution

text: "What is the difference between "significance" and "contribution" of individual variable to a outcome variable? How to address this while creating a prediction model basically.I understand contribution made by a variable towards the outcome, which shoudl be considered while creating a model but not sure of how "significance" is different from this. Thanks a ton.


title: Computing baseline probability

text: What is the meaning of the term Baseline probability in an experiment? How is it computed, say for a binary classifier? How to measure the performance of a classifier according to a given Baseline probability?


title: Probability usage in a customer event

text: I have a probability question for our customer event. Users choose tickets that contains one word from alphabet {A,E,L,P}. The ticket is not given back and user can choose as many tickets as he wants. There is no maximum number of involved users (but in real there will max. 1000). In an event we need there must be min. 20 winners = 20 people have collected word "APPLE". User that wants to win he needs to collect word APPLE. Obviously for word "APPLE" they need to collect five tickets - 1xA, 2xP, 1xL, 1xE. **Question:** How many tickets (with one letter from mentioned alphabet) we should print to be satisfied that there will be in total minimum 20 winners? And how many particular letter we should print? (e.g. letter P must be selected twice for one word)


title: Data transformation and confidence intervals for mean difference

text: I have a between-group independent variable with two level (A and B) and a dependent variable _Y_ that I transformed in order to normalize the distribution of the residual. I used a Box Cox transformation with the following formula : (( _Y_ ^3,169833)-1)/3,169833 With this transformation I observed the following results ( _Y_ ' refers to the transformed variable) ![enter image description here](http://i.stack.imgur.com/knunF.png) I used this formula to back transform my data:( _Y_ '*3,169833+1)^(1/3,169833) and I obtained the means and CI in red in the figure. My question is how can I get the difference between my two means and the 95% CI for that difference? I tried to calculate the difference between the trasnformed means, to get its 95% CI and then to transform back. However I obtained an improbable result. I guess that maybe this is because I transformed _Y_ and not the difference between my to groups. Is there a formula to compute the CI of a difference between means from the two CI of the two means? Thank you very much in advance Best regards **Update after Alexy's answer** I have a question about the fact that "inferences about f(x) are not the same thing as inferences about x". I agree with that. However I thought that it was for that reason that we have to back transform the transformed data in order to interpret the result more easily. For example, Bland & Altman (1996)* have suggested that we can back-transformed the confidence limits of a log transform data using the anti-log in order interpret the results in the original unit. Do you agree with this possibility? Or does the fact that "inferences about f(x) are not the same thing as inferences about x" preclude this possibility? I you agree with Bland & Altman (1996), do you think that we can approximate the sd for x by following the following steps? ![enter image description here](http://i.stack.imgur.com/5fNS0.png) Does it make a little sense or do I am totally wrong in doing that? *Of course I have forgotten that Bland & Altman (1996) have underlined the impossibility to back transform the standard deviation. **2nd Update** So I tried to compute the sd for back-transformed _Y_ from the back- transformed CI and the critical _t_ value and the square root of the sample size (see the results above). I wanted to find a method in order to check the validity of this method. My assumption was that if the method is valid, the inference made with transformed data would have to lead to the same conclusion as the inference made with the back-transformed data. If so the t-test with the two sets of data should lead to the same results and a standardized version of the effect size estimate (point and CI) should also be similar with the two sets of data. So I performed a t-test for independent samples with the transformed data and another t-test for independent samples with the back transformed data (using the sd computed from the above-mentioned formula). I also computed the Cohen's _d_ and their 95% CI from the two sets of data. The following table presents the results of my analyses. ![enter image description here](http://i.stack.imgur.com/Zx4kU.png) From that results it seems that the inference made with the back-transformed data and the approximate sd is consistent with those made with the transformed data. Of course the inference made with the back-transformed data is very very slightly liberal compared with the inference made with the transformed data. A more central issue with the formula is that as Alexy said, "inferences about f(x) are not the same thing as inferences about x". However, while keeping this fact in mind, it is easier to interpret the back-transformed data than the transformed data and since the MOE is similar for the two sets of data why we could not use the back-transformed data? I aknowledge that the obtained back-transformed mean difference is not equivalent to the arithmetic mean difference. So I think that if we decide to use this method to compute the back-transformed mean difference and its CI, we have to remind that fact to the reader. What do you think about my solution to the problem? Do you see other limits of my solution? Do you think that I underestimate the limits that I mentioned previously? Thanks again for your time and for your answers. **Reference** Bland, J. M., & Altman, D. G. (1996). Statistics notes: Transformations, means, and confidence intervals. _BMJ_ , 312(7032), 1079–1079. doi:10.1136/bmj.312.7038.1079.


title: how to create labeled term to document matrix for classification for TMG

text: I already asked the same question but not got any reply. Please any one help me. I have created term-to-document matrix for document collection. I want to classify them with TMG (Term to matrix generator). For this I have to input labeled term-to-document matrix. I don't know how to created a labeled such matrix. Pleas anyone, if have idea, do reply.


title: Does residual autocorrelation beyond the first lag have any implication in a regular regression?

text: I have seen the review of a multiple regression analysis using time series with a quarterly frequency. The original modeler advanced that the model's residuals were not autocorrelated by disclosing a Durbin Watson score that was reasonably close to 2.0. But, the model reviewer using an autocorrelation analysis more often seen in ARIMA, tested both the autocorrelation and partial autocorrelation going back up to 15 quarter lags. Next, he calculated the Confidence Interval associated with such autocorrelation. Given that the sample was reasonably large with about 250 quarterly-data observations, the resulting CI of those correlations was pretty narrow. As a result, the model reviewer indicated that the regression model suffered from residual autocorrelation because some of the lags beyond t-1 were associated with residual autocorrelation that were deemed statistically significant (beyond the 95% CI, where an autocorrelation could be deemed not different than zero). I thought the above residual analysis was unconventional and more appropriate for an ARIMA modeling effort. I also wondered about the relevance of any residual autocorrelation beyond the 1st quarter lag. The latter is the only one associated with truly chained or linked datapoints. At the 1 quarter lag, if the residuals are reasonably highly correlated, they could definitely trend by having the same sign for numerous consecutive periods. This is much less the case or often not at all with any greater quarter lag. The residuals will not trend as indicated. I figured this is probably why the DB score focused on the first period lag appears to be the dominant metric of residual autocorrelation. Also, let's say you do have residual autocorrelation even upon the very first lag. Is'nt there a rule or something that states if you have a large enough sample (maybe > 100 observations). The residual autocorrelation will not matter anyway because due to the large sample the statistical significance of the regression coefficients will be pretty much unaffected. This may be an offshoot of the Central Limit Theorem. What do you think?


title: Does heteroskedasticity matter if you have a large enough sample?

text: Let's say you run a regression with over 200 observations. Would this reasonably large sample mitigate the impact of residuals heteroskedasticity as an offshoot of the Central Limit Theorem, or something similar. As a result, the statistical significance of the regression coefficients would not be in question.


title: significant difference in mineral concentration between chemical analysis and values declared on packaging?

text: I want to find out whether there was any significant difference in mineral concentration between chemical analysis and values declared on packaging. After carrying out Shapiro-Wilk tests to assess the normality of the data obtained in the lab, results indicted that data was both parametric and non- parametric. I have carried out a sample t test to assess the parametric data. However, what non-parametric test would be suitable to find whether there was any significant difference in mineral concentration between chemical analysis and reference values?


title: Visibly different curves but non significant likelihood ratio test

text: I have run a series of probit regressions for growth variables across 3 environments. I then graphed the resulting curves as probability density functions (3 curves per graph one for each environment). Then I ran likelihood ratio tests to determine whether grouping by environments was a better fit than not grouping. I have done this for a series of variables. With one exception, my LRT are non significant and these results accord well with the graphs (the curves heavily overlap). But for one, my curves are well defined and have clearly separated peaks. Two of the curves overlap a bit, but not as much as in other graphs, and the third curve is clearly separate and the distribution is visually less variable than the other two. These curves are visually more different than in any of the graphs with significant LRT results. My question is why would the LRT not be significant if the curves are so visibly different from one another (even more so than some of the significant models)?


title: The Intercept terms for Ridge regression, lasso , pcr and PLS (elements of statistical learning)

text: In table 3.3 (page 63) of the elements of statistical learning book, the intercept terms for Ridge regression, lasso , pcr and PLS differ. However, according to the theory in the book, these models should all have the same $\hat{\beta_0} = \bar{y}$. How are the intercepts estimated in the table ? Note that: all these models are applied onto the same dataset, where the inputs are centered. ![enter image description here](http://i.stack.imgur.com/DD4LY.png)


title: Linear regression slope confidence interval - how to choose a cutoff with few points?

text: I am performing a large number of linear regressions (around 2000) through the origin. Each regression is on a different number of points (between 2 and 1000), and I am using robust regression (with the rlm package in R) to estimate the slope of each regression line. Clearly I should not trust the results of the regressions on 2 points, but I am trying to decide which regressions have results that are "good enough" and which do not. I can calculate the 95% confidence interval of the slope $\beta$ as: $\beta \in [\hat{\beta} - s_\hat{\beta}*t_{n-2} , \hat{\beta} + s_\hat{\beta}*t_{n-2}] $ where $s_\hat{\beta}$ is the standard error of the slope, and $t_{n-2}$ is the 97.5th quantile of a t distribution with n-2 degrees of freedom. I have calculated the 95% confidence intervals for each of the 2000 regression slopes, but how can I choose on a reasonable cutoff for the 95% interval such that the regressions above the cutoff are "good enough"? Should I be comparing ratio of the width of the confidence interval over the slope for different regression rather than just the width of the confidence interval?


title: Jointly estimation of model or single step estimation?

text: I have financial data (return) and use the following model: ![returnmodel](http://i.stack.imgur.com/iUeit.png) where $R_t$ is the return at time t, $\mu$ is set to zero, $\sigma$ is the volatility and $\epsilon$ is an innovation process. The $\sigma^2$ is estimated with a volatility model. The distribution of the innovations is mostly Gaussian, more advanced models use e.g. generalized hyperbolic distribution. My question is now, how to estimate the parameters of the models? For example, consider the easiest combination: The volatility at timepoint t is estimated with a simple moving average, which uses exponentially weights. The innovation process is said to be Gaussian. I do the estimation as follows:   1. I calculate the volatility at each time point using the original return series.   2. I calculate the standardized residuals ($R_t/\sigma_t)$ and use these values to estimate a normal distribution using ML.   3. I look at the estimates of ML, if the model is correct specified, the mean should be zero and the standard deviation/variance equal to zero, because the standard model assumes the innovation process to be Gaussian with mean zero and Variance of one.   4. Likely this will not be the case, so I use another model, e.g. the generalized hyperbolic distribution and hopefully see, that the generalized hyperbolic distribution seems to match the first four empirical moments better than the Gaussian distribution So in this case, the procedure would be the same but If I assume the volatility process to follow a GARCH process and the innovations to follow a certain distribution, e.g. Gaussian, generalized hyperbolic etc. do I have to do jointly estimation? So that I use a method of R which does the estimation of the volatility process and the innovation process? Or Can I do the following, according to the steps above: 1. Calculate the volatility using a GARCH process. 2. Use the standardized residuals ($R_t/\sigma_t$) to do ML estimation of the distribution Will this "single step" estimation procedure lead to the same estimates as the "jointly" estimation procedure, which most GARCH procedures do?


title: Need help finding UMVUE for a Poisson Distribution

text: My fellow classmates and I are stuck on a homework problem that is a three part problem to find the UMVUE of a Poisson distribution. The problem goes like this: Let X ~ Pois$(\lambda$), and we want to estimate $\theta=e^{-\lambda}$. a) Find the Rao-Cramer lower bound for an unbiased estimator for $\theta$ (based on a sample of size n). Hint: The differentiation needs to be with respect to $\theta$, not $\lambda$. So the majority of us used the Fisher Information Index to get: $I_1(\theta) = \Sigma ((\frac{d}{d\theta}ln(e^{-\lambda}\frac{\lambda^x}{x!}))^2*e^{-\lambda}\frac{\lambda^x}{x!})$ We then substituted in our ln $e^{-\lambda}$ for $\theta$. $I_1(\theta) = \Sigma ((\frac{d}{d\theta}ln(\theta\frac{\lambda^x}{x!}))^2*e^{-\lambda}\frac{\lambda^x}{x!})$ After that, we broke out the ln function: $I_1(\theta) = \Sigma ((\frac{d}{d\theta}ln(\theta) + ln(\lambda^x) - ln(x!))^2*e^{-\lambda}\frac{\lambda^x}{x!})$ Then we took the derivative with respect to $\theta$. $I_1(\theta) = \Sigma ((\frac{1}{\theta})^2*e^{-\lambda}\frac{\lambda^x}{x!})$ This is where we have a little difference of opinion because one of us came up with the final answer as being: $\frac{ne^{\lambda}}{\theta}$ He took the final answer by cancelling one of the $\theta$s and then using the Taylor series to get $e^{\lambda}$ Where two of us got the answer as being: $\frac{n}{\theta^2}$ We don't know which is right because either one could make sense. If someone could please let us know on this, we would be grateful. Part b we are also stuck on as well. b) Show that $T=\Sigma X_i$ is a complete sufficient statistics for $\theta=e^{-\lambda}$ Now, from our notes, we have that a probability distribution belongs to an exponential family if it can be written in the form $e^{A(x)B(\theta)+C(X)+D(\theta)}$ So we start out and we get: $e^{-\lambda}\frac{\lambda^x}{x!}$ $=e^{-\lambda}e^{ln(\frac{\lambda^x}{x!})}$ $=e^{-\lambda}e^{xln(\lambda)-ln(x!)}$ This is where we get stuck because we are not sure 1) where to put the T into this equation and 2) how are we supposed to relate this back to the theta. So far, we have only come up with: $=e^{ln(\theta+\Sigma X_i ln(\lambda)-ln(\Sigma X_i!)}$ We know we are probably completely and totally wrong on this. When it comes to part c, of our equation, we are a bit stuck as to what to do. Part c goes like: If we use $T^*=(\frac{n-1}{n})^{\Sigma X_i}$ as an estimator for $\theta=e^(-\lambda)$, as the Rao-Blackwell Theorem indicated, combine your findings in Parts (a) and (b) to determine Var($T^*$). We know that this is a UVMUE that can probably be solved by the Lehmann- Scheffe theorem. The thing is that we just don't know how to really apply it or how we can used what we find from parts a and b to get Var($T^*$). If anyone can help out with any parts of this problem, we will be very grateful because we are very, very, very stuck on this. Thanks!


title: Derivative of the modified Bessel function of third kind with respect to its order in R

text: How can I find the derivative of the modified Bessel function of third kind ($K_{\nu}(x)$) with respect to its order ($\nu$)? Specifically, how can I evaluate \begin{equation} \frac{\partial K_{\nu}(x)}{\partial \nu} \end{equation} in R, where $\nu$ and $x$ are assumed real to be real. R package Bessel evaluates Bessel functions for a variety of arguments ($x$) and order ($\nu$). There is also R package fAsianOptions that evaluates the Bessel functions for integer order and has functions for evaluating \begin{equation} \frac{\partial K_{\nu}(x)}{\partial x} \; (\text{note } x \text{ not } \nu). \end{equation} A related, but not similar, question is here. The recurrence relation between the $\frac{\partial K_{\nu}(x)}{\partial \nu}$ and other Bessel function can be found here. Using asymptotic expansions \begin{equation} \frac{\partial K_{\nu}(x)}{\partial \nu} = K_{\nu}(x)\bigg[1 - \frac{1}{2\nu} - \log \big(\frac{ex}{2\nu} \big)\bigg]. \end{equation} I was wondering if there are better approximations.


title: How to model a biased coin with time varying bias?

text: Models of biased coins typically have one parameter $\theta = P(\text{Head} | \theta)$. One way to estimate $\theta$ from a series of draws is to use a beta prior and compute posterior distribution with binomial likelihood. In my settings, because of some weird physical process, my coin properties are slowly changing and $\theta$ becomes a function of time $t$. My data is a set of ordered draws i.e. $\{H,T,H,H,H,T,...\}$. I can consider that I have only one draw for each $t$ on a discrete and regular time grid. How would you model this? I'm thinking of something like a Kalman filter adapted to the fact that hidden variable is $\theta$ and keeping the binomial likelihood. What could I use to model $P(\theta(t+1)|\theta(t))$ to keep inference tractable? _Edit following answers (thanks!)_ : I would like to model $\theta(t)$ as a Markov Chain of order 1 like it is done in HMM or Kalman filters. The only assumption I can make is that $\theta(t)$ is smooth. I could write $P(\theta(t+1)|\theta(t)) = \theta(t) + \epsilon$ with $\epsilon$ a small Gaussian noise (Kalman filter idea), but this would break the requirement that $\theta$ must remain in $[0,1]$. Following idea from @J Dav, I could use a probit function to map the real line to $[0,1]$, but I have the intuition that this would give a non-analytical solution. A beta distribution with mean $\theta(t) $ and a wider variance could do the trick. I'm asking this question since I have the feeling that this problem is so simple that it must have been studied before.


title: Regression and simple SEM: identification and model comparison

text: I have two problems I am trying to solve. I am using MPLUS, but I have very general questions for which MPLUS knowledge will not be necessary. Terminology: y = DV, x = IV **(1) Simple regression (N=1300)** I want to compare two linear regressions. In the first model, a metric DV is predicted by 14 IVs as well as 2 covariates (age, sex). In the second model, regression coefficients of the 14 IVs are constrained to be equal, only age and sex are freely estimated. In MPLUS, a simplified example would look like this ("on" marks regressions, and the brackets set regression weights equal): Model I:               y on x1;     y on x2;     y on x3;     ...     y on sex age;      Model II:               y on x1 (1);     y on x2 (1);     y on x3 (1);     ...     y on sex age;      The first regression has 0 DF. Therefor, fit cannot be assessed, which in turns makes model comparison impossible. How do I get the first model identified so I can compare the two models? Does a linear regression always have 0 DF? I thought about perhaps constraining two of the fourteen regression weights that are very close to each other equal - is that a valid or common approach? For instance, x2 and x14 have a similar nonsignificant regression coefficient estimate of .050 which I could constraint equal to identify the model. If I do that, I get a chi-square of 0.000 for that model with 1 DF -- can that be correct? The constrained model has a chi-square of 398 with 13 DF. My next question would be how to compare these models. MPLUS outputs AIC, BIC, RMSEA, SMRM, CFI/TLI, and I understand how to interpret them generally. However, am I allowed to interpret them for a simple regression? Am I allowed to interpret them in a model that is not identified? And as for the model test, is it correct that I simply use the chi-square values in the two models, build the difference, multiply it by 2 and look up the value in a chi-square table for the difference between the degrees of freedom of the models? **(2) SEM (N=3500)** The second problem is somewhat similar, but involves more complex models. I am new to structural equation modeling, but assume these models to be SEMs (?). Model I consists of 9 regressions with metric DVs, each regression has the same 6 predictors. The model allows for all variables to be correlated (also y with y and x with x). The DV is a the second measurement point of y (yt2), and the regressions control for the first measurement point of (yt1). So, we're really predicting changes in y. Model I estimates all regressions freely, while model II constrains each x to affect all y in the same way. Again, a simplified example: Model I:               y1t2 on x1 x2 x3 y1t1;     y2t2 on x1 x2 x3 y2t2;     ...     y9t2 on x1 x2 x3 y9t2;      Model II:               y1t1 on x1 x2 x3 y1t2 (1 2 3 4);     y2t1 on x1 x2 x3 y2t2 (1 2 3 5);     ...     y9t1 on x1 x2 x3 y9t2 (1 2 3 6);      Is it valid to compare these two models in the way mentioned above (chi-square difference test)? Can fit indices be interpreted normally? Model I has a chi-square of 197.293 and 72 DF, model II 617.376 and 120 DF, so that would be 120-72 DF, 2*(617.376-197.293) --> p<.001 ?


title: Error in cor(y): allocMatrix: too many elements specified

text: x <- read.table('file_name',header=TRUE, row.names=1)     y <- t(x)     y <- data.frame(y)     row.names(y) <- names(x)     names(y) <- row.names(x)     library(corrplot)     corr <- cor(y)      par(ask = TRUE)     corrplot(corr, order = "hclust")      I'm trying use corrplot on my dataset. The original dataset has 25000 rows and 100 columns. I transposed to 25000 columns and 100 rows and finally applied corrplot. And failed! If any help would be appreciated. [note:1st column is and first row are names and others are values (like -0.9 or 2.5 or 14449.45 etc.)]


title: Is it possible to run Arellano-Bond estimator when i have three time periods?

text: I have a hard time figuring out if I can use Arellano-Bond GMM when I hahve 3 time periods. In the first period there is 37 observations, 47 in the second period and 60 in the thrid period. However, the way I understand it, to run the GMM Arrelano-Bond it requires 4 time periods? I tríed to run it with two lags in stata and get the following message:   J(): 3300 argument out of range _xtdpd_getpinfo(): - function returned error : - function returned error When I run it with 1 lag i get the following message: artests() cannot exceed the number of time periods


title: probability of the next number of events occurring when the true probability is known

text: > **Possible Duplicate:**   >  Probability of getting between What is the probability that at least 24 of the next 50 people like to swim when the true probability of people liking to swim is 35%? I thought it might be a conditional probability but I'm not really sure...


title: Accuracy and other measures of performance for Self-Organizing Map

text: I have a Self-Organizing Map and I would like to compute something like its accuracy, precision and specificity to evaluate the performance of my classifier in prediction. About accuracy, I used two indices: Adjusted Random Index and Adjusted Mutual Information. Are they good choices? How can I also compute specificity, precision and recall?


I know that some of the questions I am asking here have been answered in a general case in the two questions I am referring to in the problem section. Nonetheless, I am asking for a very specific case here. --> The representation of a high-dimensional data set by a low number of data points **Background:** I am currently working on a European electricity sector model, which should allow the user to model the development of the sector for the foreseeable future. As the model can grow very big with a lot of input data (and thus be almost infeasible), I am trying to reduce my input data to a minimum. Up until now this problem is solved by choosing snapshots of the data set in a look-on fashion and thus choosing by a sense of proportion and not a scientific method. **Problem:** In particular, I am trying to best represent a high-dimensional (numerical) data set, with 144 dimensions and 8760 data points, by a small number of representative data points via the use of a (or better some) data- mining method(s). My data set actually consists of 144 time series, each describing a feature of the European electricity sector. In other words, I have electrical load, Photovoltaics-infeed, Wind-On- as well as Wind-Offshore-infeed for 36 European countries (hourly --> 8760h). Ergo 4*36 = 144 dimensions. As I already described, I would love to find a (or better some) best-possible representation(s) of these 8760 hours in order to significantly reduce the number of points in time. The best-case solution would be, that I would be able to find multiple methods in order to solve the above-mentioned Problem. A next step would be to compare the found methods and their solutions, respectively... **Status-Quo:** Unfortunately I am rather new to the data-mining/unsupervised- learning topic. I read a lot over the last few weeks and tried to get acquainted with different methods. Nonetheless, I have serious doubts as to the feasibility of my undertaking. Having a look at this question 1 and this question 2 made me wonder how much use "classical" methods, meaning partitioning and hierarchical methods can be in my case. Especially regarding the curse of dimensionality. In order to get around this specific problem I ended up reading a lot about sub-space clustering. As I am pretty sure that all "my" dimensions are relevant and I need a snapshot of all dimensions at any given time, I don't see how this kind of clustering can have viable (wholesome) solutions to the problem. Somehow, due to my skepticism as to what extend my problem is really feasible with regard to the usage of unsupervised learning I am feeling a little bit like I am at a dead end here and don't really see how this problem can be solved in a meaningful manner. Because of that, I would be very delighted if some of you could help me with my questions and in doing that shed some light into the dark of Clustering high-dimensional data sets. **Questions:** I) Does it make "sense" to use a (plain) partitioning and/or hierarchical method in my specific case? (Sorry, if I cannot infer this from the answers to the other questions. But I really mean in my particular case!) II) Am I wrong to assume that sub-space clustering could give a lot of clusters in subspaces, without really computing clusters (in all dimensions) that allow me to wholesomely consider the European electricity sector? III) Apart from the already mentioned methods, which other methods could be feasible in this particular setting? What can I read into?


title: Help clarify the implication of normality in an Ordinary Least Square (OLS) Regression

text: > If the dependent variable is normally distributed for a fixed set of > predictor values, then the residual values should be normally distributed > with a mean of 0. I have two questions based on the above paragraph: First, I don't understand why the above paragraph must be true. Could someone please provide more details why the above must be true? Would be great if the answer clearly states what properties of random variables are being used to arrive at the above conclusion and any other critical assumptions not explicitly stated in the paragraph. For simplicity, please assume a simple linear regression. The second question is I am not really sure what the phrase "...for a fixed set of predictor values.." mean. Does this mean that: Suppose we have the following pairs of observations ($x$ is the predictor value, $y$ is the dependent value): $$( x_1 , y_1) , (x_2,y_2), (x_3,y_3)$$ Does this mean that for each $x_i$, there will be corresponding predicted value $y_i$ and all of these $y_i$ will follow a normal distribution with possibly different means for $y_i$ but same variance, where $i \in [1,3]$?


title: Detecting heterogeneity in groups

text: This is a mockup of a dataset I am currently working on:               > df <- data.frame(name      = LETTERS[1:3],                        trials    = c(40, 60, 60),                        successes = c(22, 30, 45))     > df$success.rate = df$successes / df$trials     > df       name trials successes success.rate     1    A     40        22         0.55     2    B     60        30         0.50     3    C     60        45         0.75      I have a few hundred groups of varied size like this one, and what I want to do is figure out if the success rates varies between individuals. My first though was to run a chi-squared test for goodness of fit to a uniform distribution, like so:               > with(df, chisq.test(success.rate))                  Chi-squared test for given probabilities          data:  success.rate     X-squared = 0.0583, df = 2, p-value = 0.9713          Warning message:     In chisq.test(success.rate) : Chi-squared approximation may be incorrect      Which I find problematic since it doesn't take into account the fact that each individual has a different number of trials (not to mention using proportions is problematic in itself here). So my next step was to run it like this:               > with(df, chisq.test(x = successes, p = trials, rescale.p = TRUE))                  Chi-squared test for given probabilities          data:  successes     X-squared = 3.3711, df = 2, p-value = 0.1853      Which is an indication that the weighted number of successes within this group is uniformly distributed. However, intuitively my boss and I see a relevant difference between individuals B and C. Both have the same number of trials and yet one has had success in 15 more cases than the other. We believe that difference is unlikely due to chance and should therefore show up on some statistic. I believe that the chi-squared test conducted above answers a different question. After all, **it's one thing to test if the heights of a group of people are normally distributed, and a completely different thing to say that the guy measuring 1,80 m is as tall as the guy with 1,60 m**. As alternatives to detecting the difference in performance of those individuals, we have already considered and tested many things from simple range checking (we would raise red flags on groups where the range of success rates was higher than 0.20) to crude confidence interval setting ($Q1 + 1.5 I\!Q\!R$, $\mu \pm 1.96 \sigma$ and other combinations thereof). I've also tried testing for other distributions (e.g., Normal through `shapiro.test(df$success.rate)`), but since my samples are usually so small it's no wonder I'll almost always get a high p-value. The amount of groups we're analyzing and the diversity of sizes (which can range from 2 to 15) has made us look towards a non-graphic and non-parametric solution, but I am unsure about which would be the most appropriate.


title: Detecting fishy data

text: Just a little thought I've been having. If we rolled a fair dice 60 times, and got 60 sixes in a row, we would (wrongly?) definitely assume that something fishy's going on. Is there any statistical measure which could raise our attention to this potential problem? Or is it simply not a problem, since this sequence is just as likely as any other sequence of 60 rolls (e.g. 10 occurrences of each of the six numbers, or 60 fives, etc.)? I'm guessing we could use the weak law of large numbers to calculate what the probability of the average observation being different from the expected value of the dice, by a value of 2.5 in this case, after 60 trials is?


title: Scoring items which are not easily compared

text: First of all, I apologize since this question has probably been asked many times and is easily answered. However, as a statistics amateur I simply couldn't figure out what keywords are relevant to my question. Suppose you have 100 merchants and 100 products. Each merchant sells a certain range of products, ranging from only one product to all 100 products. Also, products are sold in widely different proportions, which differ among merchants, and are subject to the merchant's individual (irrational) preferences. Whenever a merchant makes a "pitch" on the market, we observe whether or not he manages to sell the product he's pitching. We assume the probability of success depends (a) on the skill of the merchant and (b) the attractiveness of the product. The products' prices are fixed, so that's not a factor. The data we have consists of millions of pitches. For each pitch, we know whether or not it was successful, the merchant, and the product. Obviously, if we compare merchants by their average success rate, this information is useless because every merchant sells different products. Likewise, if we compare products, we gain no information since every product is sold by different merchants. **What we want is a skill score for each merchant, which is independent of the products the merchant is selling, and an attractiveness score for each product, which is independent of the merchants who are selling it.** I don't need a comprehensive explanation, just some keywords to point me in the right direction. I literally have no idea where to start. Edit: Note that our assumption is that the product attractiveness is merchant- independent and the merchant skill is product-independent, i.e. there are no merchants which are better at selling certain products but worse at selling others.


title: Survival analysis - Using parameter estimates of a fully-parameterized PH model into an AFT model

text: I want to predict an event based on transactional data. I'm only interested in predicting either future hazard or time-to-event. My understanding is that to predict in the future I have to use a fully-parameterized proportional hazard model with time-varying covariates. By "fully parameterized" I mean that I assume the underlying distribution of the event in time. So the PH model has two parts: the baseline hazard (whose exact formula I assumed, so this is not a Cox model) and the effet of the covariates (`e^(linear combination of the covariates)`). The value of the second part (the effect of the covariates) is varying in time until now; after that, I assume the value of the covariates stays constant. As for the first part, the baseline hazard, I know its function, so I can extrapolate the value of the baseline hazard in the future. My first question: does it make sense to make predictions in the future in that manner with a PH model? Also, I'm interested in predicting not only the hazard but the time-to-event. For that I need to create an AFT model. My understanding is that even though theoretically possible it is computationally too intensive to create an AFT model with time-varying covariates. So I want to build an AFT model with the last slice in time of my transactional data set. But I want to take into account my whole history to get my parameters. The idea would be to obtain the survival function $S(t) = e^{(-h(t) \cdot t)}$ from the hazard function $h(t)$ of a fitted PH model. Isolate t in the formula and express it in function of $S(t)$ and $h(t)$: $T(t)=- \frac{ln(S(t))}{h(t)}$. Then I fix t at the last point in time, so the function doesn't depend on time but is using parameter estimates that were derived by taking time into account. It is important to note that I'm not actually fitting the AFT model. The reason why I'm trying to use the estimates obtained with the fully-parameterized PH model with time- varying covariates into a non-fitted AFT model is that I believe these betas are "better" since they were derived based on a larger data set that takes the time dimension into account. My second question is: can I do that or is there a silly step in the above development? Thanks in advance for your help and sorry for the long question!


title: In Stata, why do the stcox CI differ when using margins?

text: I'm trying to understand how the `margins` and `marginsplot` work after `stcox` command in Stata. So let's start with a dumb example:               webuse stan3     stset     stcox i.posttran i.surg      That gives me result               ------------------------------------------------------------------------------               _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------       1.posttran |   1.173034   .3444712     0.54   0.587     .6597023    2.085801        1.surgery |   .3449043   .1482592    -2.48   0.013     .1485267    .8009264     ------------------------------------------------------------------------------      Now, when I use               margins, over( posttran surgery )      The output is               Predictive margins                                Number of obs   =        172     Model VCE    : OIM          Expression   : Relative hazard, predict()     over         : posttran surgery          ----------------------------------------------------------------------------------                      |            Delta-method                      |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]     -----------------+----------------------------------------------------------------     posttran#surgery |                 0 0  |          1          .        .       .            .           .                 0 1  |   .3449043   .1482592     2.33   0.020     .0543217    .6354869                 1 0  |   1.173034   .3444712     3.41   0.001     .4978826    1.848185                 1 1  |   .4045843   .2057569     1.97   0.049     .0013083    .8078604     ----------------------------------------------------------------------------------      What is puzzling me here is the fact that estimates and their Std. errors are the same as when using `stcox`. However - the CI, z and _p_ are now different. What is causing such behaviour? * * * Update: I also posted the question on Statalist.


title: Which Statistical Test

text: I have a variable (VAR1 )with two possibilities, let's call them V and N. I asked my experiment participants to detect the VAR1 value. My aim is to compare the correctness rate of answers by participants (e.g. when the value of VAR1 is V, the participants can answer V or N, but V is correct, so if 70 percent choose V, correctness rate is 70%). To explain more: These are all the possibilities: VAR1 value: V Correctness of Detection by participants: C1 VAR1 value: N Correctness of Detection by participants: C2 Question: Which statistical test should I use to see in which case Correctness rate is significantly more? What will happen if I add another independent variable VAR2 and compare the whole result? like this: VAR1 value: V VAR2 value:V Correctness Rate: C1 VAR1 value: V VAR2 value:N Correctness Rate: C2 VAR1 value: N VAR2 value:V Correctness Rate: C3 VAR1 value: N VAR2 value:N Correctness Rate: C4 In which category Correctness Rate is significantly more?


title: Categorical logit Predictor with too many different levels

text: One of the predictors I had in a logit model is "City". Problem is this categorical variable has too many factor levels. e.g. In a Sample of $\sim 3000$ there are already $\sim 200$ different cities. Is it fair to still retain City as a predictor or should I purge it entirely from the model? An alternative is to retain, say, the top five most common cities and then code all the rest as a new factor level "Others". The top city occurs $\sim 60$ times but by the fifth common city this occurance drops down to 30. Some cities occur only once or twice in the dataset. PS. One problem I face (if I retain all factor levels) is that certain levels occur in the validation set but not in the training set. Then the model complains at validation time. PPS. On more reading, I found sugesstions to use combine.levels() from the Hmisc package in `R`. Maybe that will work, though not sure how exactly yet. Is there an elegant way to deal with these issues?


title: Do we need Overlap/Common Support in case of a parametric regression?

text: If I want to make a causal statement based on selection on observables. One typically assumes "Common Support" (/"Overlap") - which means that for any value of the confounding variables X a unit i can be potentially observed with treatment (D=1) and without treatment (D=0). I.e.:               0<P(D=1|X=x)<1       where `P(D=1|X=x)` denotes the probability of receiving treatment conditional dependend on the confounding variables X. In case of non-parametric (semi- parametric) estimation (matching on X or on the propensity score) this assumption is crucial. However, I am wondering whether this assumption has to hold if I want to estimate the treatment effect in a parametric regression (e.g. a simple multivariate linear model fitted by OLS). Then the assumptions on the error term (e.g. normality) enable us to extrapolate the results/treatment effects to regions without common support anyway - hence, we might also consider the units i without common support in the estimation(?). Please let me know about the necessity of this assumption in case of the parametric estimation. Thanks in advance.


title: Heterogeneous Treatment Effects - How to test differences in the ATE?

text: I want to conduct a simple propensity score estimation where the treatment $D_i$ is a binary variable ($D_i=1$ individual $i$ participates in the labor market program, zero otherwise). I estimate the propensity score using a simple probit model using various explanatory variables (including gender). Then I want to calculate the average treatment effect (ATE) (and the average treatment effect on the treated, ATET) using e.g. caliper matching with replacement or a simple 1:1 matching with replacement. The question: how can I check whether the ATEs of e.g. male and female individuals differ, i.e. whether $$(\text{ATE|i is male} - \text{ATE|i is female}) = 0$$ How can I make inference (whether the difference between the ATEs is different from zero)? Would you recommend bootstrap? If this was the case how would I do it?


title: Getting the number of active employees per month

text: In large companies with thousands of employees, how to get the number of employees per month? everyday there could be new employees, resigned or terminated employees, etc. If I want to get the number of the employees in a given month it will be confusing since the number of employees differs from one day to another in the same month. Do I just calculate the average number (total employee in all days/total days) or can I get the highest and lowest and divide by two? which is better or is there a better and more accurate way?


title: Linear Discriminant Analysis and non-normal distributed data

text: If I understand correctly, a Linear Discriminant Analysis (LDA) assumes normal distributed data, independent features, and identical covariances for every class for the optimality criterion. Since the mean and variance is estimated from the training data, isn't it already a violation? I found a quotation in an article (Li, Tao, Shenghuo Zhu, and Mitsunori Ogihara. “Using Discriminant Analysis for Multi-Class Classification: An Experimental Investigation.” Knowledge and Information Systems 10, no. 4 (2006): 453–72. ) > "linear discriminant analysis frequently achieves good performances in the > tasks of face and object recognition, even though the assumptions of common > covariance matrix among groups and normality are often violated (Duda, et > al., 2001)" -- unfortunately, I couldn't find the corresponding section in Duda et. al. "Pattern Classification". Any experiences or thoughts about using LDA (vs. Regularized LDA or QDA) for non-normal data in context of dimensionality reduction?


title: VAR / VEC in levels or difference depending on Cointegration

text: Thanks in advance. I have four I(1) variables I'm trying to model by VAR/VEC. I know that it is only okay to model non-stationary variables in levels only if they are cointegrated. What I would clarification on is if they ALL have to be cointegrated with each other or is it sufficient that there only needs to be one cointegration? When I run a Johansen Test (eigenvalue statistics) I get 2 cointegrating vectors. Obviously it doesn't tell me which ones, I'll have to probably do the Engle-Granger test for that. In any case, is this enough to run a VAR or VEC in levels for these variables if out of four variables I get only two error- correcting vectors? Thanks.


title: How can I find a level $\alpha$ most powerful test?

text: Let $X_1,X_2,\dots ,X_n$ be a random sample from a distribution with pdf $f(x,\theta)$. Find a level $\alpha$ most powerful test of $H:\theta=\theta_0$ against $K:\theta = \theta_1$ when $$f(x,\theta)=\theta x^2I_{(\theta, \infty)}(x),\theta_0 \neq \theta_1$$ I consider $$\frac{p_1(x)}{p_0(x)}=\dfrac{\theta_1^n(\Pi x_i)^2 \Pi(I_{(\theta_1, \infty)}(x_i))}{\theta_0^n(\Pi x_i)^2 \Pi(I_{(\theta_0, \infty)}(x_i))}\\=(\theta_1/\theta_0)^n\dfrac{\Pi(I_{(\theta_1, \infty)}(x_i))}{ \Pi(I_{(\theta_0, \infty)}(x_i))}$$ Then I need to consider $\frac{p_1(x)}{p_0(x)}>k$ for a Neyman-Pearson test but I am stuck. Please help.