title: How to measure the strength of association between binary variables where the majority of values is zero?

text: I have two binary variables, $X_1$ and $X_2$, which take value in {0, 1}. The overwhelming majority of their values are 0. I wanted to know if there is a way to appropriately measure their association given that the majority of their values are zero. Note that I am also opened to dependence tests as well. I was looking at the _phi_ coefficient, but it seems at a high level, the _phi_ coefficient is proportional to the product of the agreements minus the product of the disagreements. In my data set, intuitively, I get the sense, thus, that the _phi_ coefficient will likely be very large. Are there more appropriate association measures or independence tests for this type of data?


title: Can I use a two proportion z-test?

text: I have a response, "Satisfaction" that goes on a scale from 0-10 (where 0=dissatisfied and 10= very satisfied), and I have a predictor variable that also uses the same scale. I recoded the response so that 0-5 is "not satisfied" and 6-10 is "satisfied". I did the same for the predictor variable. I used SPSS to get the frequencies for each (there are 4 different combinations). My question is, can I use a two proportion z-test to test whether there is a significant difference in the predictor at the two different levels of my response? Thank you


title: Why do random effect models require the effects to be uncorrelated with the input variables, while fixed effect models allow correlation?

text: From Wikipedia > There are two common assumptions made about the individual specific effect, > the random effects assumption and the fixed effects assumption. The random > effects assumption (made in a random effects model) is that the individual > specific effects are **uncorrelated** with the independent variables. The > fixed effect assumption is that the individual specific effect is > **correlated** with the independent variables. If the random effects > assumption holds, the random effects model is more efficient than the fixed > effects model. However, if this assumption does not hold (i.e., if the > Durbin–Watson test fails), the random effects model is not consistent. I was wondering why random effect models require the random effects to be uncorrelated with the input variables, while fixed effect models allow the effects to be correlated with the input variable? Thanks!


title: How to analyse a Likert-type scale with a very small sample size (n=15)?

text: **Data description:** I have data from questionnaire data on 15 participants (small sample). Parents rated the frequency that they rough and tumble with their child: (1) Never, (2) not very often, (3) once per month, (4) once or more per week, (5) every day. **Initial analyses:** So far I have created several bar charts and am explaining them using percentages (e.g., 78% of the boys rough and tumble with one or other or both of their parents once or more per week this figure is only 60% for the girls). ### Questions   * Should I use medians instead of means?   * Is the sample size of 15 too small for meaningful results?   * Can I perform a meaningful Chi-square or other statistical test?


title: Markov model for time series, going back n periods?

text: I realise there are a lot of questions about markov here, but as we say in Dutch, I couldn't see the threes through the forest. I have a sequence of intervals between subsequent notes (pitches). Like: 2 -3 -1 0 2 3 etc I suspect that the next interval is not only dependant on the previous interval, but let's say, dependant on the previous two intervals. (Later on, I might also look at intervals that occur at the same time, so other voices). Is Markov still appropriate for this? Or more specifically hidden markov right? Or should I use some totally different approach? How would I go about this? I can't just summarise occurances of intervals to create the TM. Should I create a fake attribute that would be two subsequent intervals? Is there some software that could help me do this easily? I have the book on R, so I could delve into that. But i wanted to be sure about my approach and possible easier solutions before I look into that in more detail. I am also reading this very interesting intro, but my brain wanted a plain and simple explanation before going too deep into mathematics.


title: AFT model with frailty for right censored life time data

text: Let say I have a 'kidney catheter' data set. Data are about the recurrence times to infection, at the point of insertion of the catheter, for kidney patients using portable dialysis equipment. Catheters may be removed for reasons other than infection, in which case the observation is censored. Each patient has exactly 2 observations. I want to fit a Accelerated Failure Time model with random effect (frailty). By using built in survival package R code `survreg()`, I can do it like this. _Log-normal AFT model with Gaussian frailty:_ `mopdel1 <- survreg(Surv(time, status)~ age+ sex+ frailty.gaussian(id), dist="lognormal", kidney)` _Log-normal AFT model with Gamma frailty:_               model2 <- survreg(Surv(time, status) ~ age+ sex+ frailty.gamma(id),  dist="lognormal", kidney)      Does variance of random effect of these two models correct? Which estimation procedure are they using to find out the variance of random effect?


title: Why does one report statistical power only when results are non significant?

text: It is well known that a study with low statistical power has a reduced chance of detecting a true effect. But if low power also reduces the likelihood that a statistically significant result reflects a true effect (as shown in Power failure: why small sample size undermines the reliability of neuroscience), why does one report statistical power only when results are non significant? As I can't comment yet, I edit my question in response to @Glen_b and @Penguin_Knight : I should have said "discuss" instead of "report". Of course, one should do the power calculation a priori. But my question is really: Why people who get small p-value with low statistical power don't moderate their findings?


title: Sample size calculation to demonstrate reliability of physical test method

text: I'm preparing a procedure for validating attribute physical test methods (e.g. vision system) by demonstrating the reliability of the test method to detect known good/bad samples.I have sourced the below sample size formula which proposes to be based off Binomial Distribution:                n = (0.5 x Xsq(1-C,2(r+1)))/(1-R)      where n is sample size, C is confidence, r is number of failures, R is Reliability Image of formula added below (may not appear): ![Image:Sample size calculation using Binomial](http://i.stack.imgur.com/Q74AZ.jpg) I can find no other reference for the derivation of this formula and am reluctant to use it without understanding where it comes from. The sample sizes calculated from this formula align with sample sizes I have used previously so I am confident it is correct but I'd really like to know how it is derived from the Binomial Distribution.


title: Should I include demographics in my regression?

text: I am doing a study to try to explain certain behaviour using 6 predictors. I have used a questionnaire to gather the data. The questionnaire included a demographics section eg age, gender, job position. I an not sure whether I should include the demographic variables in the regression or not. I didn't originally plan to since my literature review focused on the 6 predictors and not on differences in demographics. I have found that two of the demographics are in fact significant in the GLM. What should I do? Any feedback is appreciated.


title: Is it meaningful to identify interaction effect when a main effect in a model (main effects only, without interaction terms) is not significant?

text: Let's say there's one DV (`Y`) and three IVs (`X1`, `X2`, `X3`), and among IVs, `X1` is a dummy variable. In the regression model without interaction terms, the results can be represented like this:               Y ~ X1 + X2 + X3          X1 : non-significant     X2 : significant     X3 : significant      In this case, is it meaningful to check some interaction terms (e.g. `X1` $\cdot$ `X2` or `X1` $\cdot$ `X3`)? At first I thought I don't have to because the main effect of `X1` indicates non-significance. But I'm afraid I'm missing something important.


title: How can I interpret interaction effect in my response surface results?

text: in my thesis i had a significant interaction effect between two variables and their main effect is insignificant. A variable had a p-value of 0.636 B variable had a p-value of 0.219 A*B had a p-value of 0.022 A is employee training practices, B is Autonomy, and the response is explicit knowledge creation. my thesis is about evaluation of the impact of quality management practices on knowledge creation in a public corporation in Jordan. how can i explain this result ???


title: Sum of independent Wishart with same degrees of freedom but different scale matrices

text: Is there any result showing that a sum of independent Wishart with same degrees of freedom but different scale matrices is a Wishart? For example, if I have two random variables: $$ Y \sim W_p(n,\sigma_1)\ \\\ X \sim W_p(n, \sigma_2), $$ where $p$ is the dimension of matrices, will $Y+X$ be equal to $W_p(n, \sigma_1 + \sigma_2)$?


title: Gradient descent vs lm() function in R?

text: I'm going through the videos in Andrew Ng's free online machine learning course at Stanford. He discusses Gradient Descent as an algorithm to solve linear regression and writing functions in Octave to perform it. Presumably I could rewrite those functions in R, but my question is doesn't the lm() function already give me the output of linear regression? Why would I want to write my own gradient descent function? Is there some advantage or is it purely as a learning exercise? Does lm() do gradient descent?


title: How do you explain the difference between relative risk and absolute risk?

text: The other day I had a consultation with an epidemiologist. She is an MD with a public health degree in epidemiology and has a lot of statistical savvy. She mentors her research fellows and residents and helps them with statistical issues. She understands hypothesis testing pretty well. She had a typical problem of comparing two groups to see if there is a difference in there risk related to getting congestive heart failure (CHF). She tested the mean difference in the proportion of subjects getting CHF. The p-value was 0.08. Then she also decided to look at the relative risk and got a p-value of 0.027. So she asked why is one significant and the other not. Looking at 95% two- sided confidence intervals for the difference and for the ratio she saw that the mean difference interval contained 0 but the upper confidence limit for the ratio was less than 1. So why do we get inconsistent results. My answer while technically correct was not very satisfactory. I said "These are different statistics and can give different results. The p-values are both in the area of marginally significant. This can easily happen." I think there must be better ways to answer this in laymen's terms to physicians to help them understand the difference between testing relative risk vs absolute risk. In epi studies this problem comes up a lot because they often look at rare events where the incidence rates for both groups are very small and the sample sizes are not very large. I have been think about this a little and have some Ideas that I will share. But first I would like to hear how some of you would handle this. I know that many of you work or consult in the medical field and have probably faced this issue. What would you do?


title: Differences between PROC Mixed and lme / lmer in R - degrees of freedom

text: _Note : this question is a repost, as my previous question had to be deleted for legal reasons._ * * * While comparing PROC MIXED from SAS with the function `lme` from the `nlme` package in R, I stumbled upon some rather confusing differences. More specifically, the degrees of freedom in the different tests differ between `PROC MIXED` and `lme`, and I wondered why. Start from the following dataset (R code given below) :   * ind : factor indicating the individual where the measurement is taken   * fac : organ where measurement is taken   * trt : factor indicating the treatment   * y : some continuous response variable The idea is to build the following simple models : `y ~ trt + (ind)` : `ind` as a random factor `y ~ trt + (fac(ind))` : `fac` nested in `ind` as a random factor Note that the last model should cause singularities, as there's only 1 value of `y` for every combination of `ind` and `fac`. **First Model** In SAS, I build the following model :               PROC MIXED data=Data;         CLASS ind fac trt;         MODEL y = trt /s;         RANDOM ind /s;     run;      According to tutorials, the same model in R using `nlme` should be :               > require(nlme)     > options(contrasts=c(factor="contr.SAS",ordered="contr.poly"))     > m2<-lme(y~trt,random=~1|ind,data=Data)      Both models give the same estimates for the coefficients and their SE, but when carrying out an F test for the effect of `trt`, they use a different amount of degrees of freedom :               SAS :      Type 3 Tests of Fixed Effects      Effect Num DF Den DF     F  Value Pr > F      trt         1      8  0.89        0.3724           R :      > anova(m2)                 numDF denDF  F-value p-value     (Intercept)     1     8 70.96836  <.0001     trt             1     6  0.89272  0.3812      _Question1:_ What is the difference between both tests? Both are fitted using REML, and use the same contrasts. NOTE: I tried different values for the DDFM= option (including BETWITHIN, which theoretically should give the same results as lme) **Second Model** In SAS :               PROC MIXED data=Data;         CLASS ind fac trt;         MODEL y = trt /s;         RANDOM fac(ind) /s;     run;      The equivalent model in R should be :               > m4<-lme(y~trt,random=~1|ind/fac,data=Data)      In this case, there are some very odd differences :   * R fits without complaining, whereas SAS notes that the final hessian is not positive definite (which doesn't surprise me a bit, see above)   * The SE on the coefficients differ (is smaller in SAS)   * Again, the F test used a different amount of DF (in fact, in SAS that amount = 0) SAS output :               Effect     trt Estimate Std Error  DF t Value Pr > |t|      Intercept        0.8863    0.1192  14    7.43 <.0001      trt       Cont  -0.1788    0.1686   0   -1.06 .       R Output :               > summary(m4)     ...     Fixed effects: y ~ trt                     Value Std.Error DF   t-value p-value     (Intercept)  0.88625 0.1337743  8  6.624963  0.0002     trtCont     -0.17875 0.1891855  6 -0.944840  0.3812     ...      (Note that in this case, the F and T test are equivalent and use the same DF.) Interestingly, when using `lme4` in R the model doesn't even fit :               > require(lme4)     > m4r <- lmer(y~trt+(1|ind/fac),data=Data)     Error in function (fr, FL, start, REML, verbose)  :        Number of levels of a grouping factor for the random effects     must be less than the number of observations      _Question 2_ : What is the difference between these models with nested factors? Are they specified correctly and if so, how comes the results are so different? * * * Simulated Data in R :               Data <- structure(list(y = c(1.05, 0.86, 1.02, 1.14, 0.68, 1.05, 0.22,      1.07, 0.46, 0.65, 0.41, 0.82, 0.6, 0.49, 0.68, 1.55), ind = structure(c(1L,      2L, 3L, 1L, 3L, 4L, 4L, 2L, 5L, 6L, 7L, 8L, 6L, 5L, 7L, 8L), .Label = c("1",      "2", "3", "4", "5", "6", "7", "8"), class = "factor"), fac = structure(c(1L,      1L, 1L, 2L, 2L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 1L, 2L, 2L), .Label = c("l",      "r"), class = "factor"), trt = structure(c(2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("Cont",      "Treat"), class = "factor")), .Names = c("y", "ind", "fac", "trt"     ), row.names = c(NA, -16L), class = "data.frame")      Simulated Data :                  y ind fac   trt     1.05   1   l Treat     0.86   2   l Treat     1.02   3   l Treat     1.14   1   r Treat     0.68   3   r Treat     1.05   4   l Treat     0.22   4   r Treat     1.07   2   r Treat     0.46   5   r  Cont     0.65   6   l  Cont     0.41   7   l  Cont     0.82   8   l  Cont     0.60   6   r  Cont     0.49   5   l  Cont     0.68   7   r  Cont     1.55   8   r  Cont


title: Is there an unbiased estimator of an Hellinger distance?

text: In a setting where one observes $X_1,\ldots,X_n$ distributed from a distribution with density $f$, I wonder if there is an unbiased estimator (based on the $X_i$'s) of the Hellinger distance to another distribution with density $f_0$, namely $$ \mathfrak{H}(f,f_0) = \left\\{ 1 - \int_\mathcal{X} \sqrt{f(x)f_0(x)} \text{d}x \right\\}^{1/2}\,. $$


title: How to interpret the output of survival analysis?

text: I am beginner of survival analysis, so I am confused about some basic ideas in this area. I think the basic assumption in survival analysis is that time of event for a subject can be estimated given all relevant time-invariable attributes of subjects. Therefore, ideally, if we can regress over the failure time correctly, we can then predict the expected failure time of subjects. Since regression is hard, researchers transformed it to ranking problem. That is, they want to **rank subjects according to their failure time**. Furthermore, the performance of ranking is measured by concordance index (c-index). My question is, "Even with a perfect ranking list (larger score for longer survival time), how can we determine the failure risk of each subject at time _T_?" @mbq, @ocram :)


title: Using expectation to detect bias

text: I was going through Penn State's online notes and noticed this expression: $ v^2 = \frac{1}{n} \sum_{i=1}^n (y_i - \bar{y})^2$ In the line below it they stated that the $E[v^2] = (1 - \frac{1}{n})\sigma^2$. I was wondering how would you get that? Would it be wrong for me to say that since $\sum_{i=1}^n (y_i - \bar{y})^2 = (n-1)s^2$ the $E[v^2] = E[\frac{(n-1)s^2}{n}] = (1 - \frac{1}{n})\sigma^2$ ? But this line of reasoning forces me to assume that $E[s^2] = \sigma^2$ and I don't even know why that's true.


title: Minimum number of observations for multiple linear regression

text: I am doing multiple linear regression. I have 21 observations and 5 variables. My aim is just finding the relation between variables   1. Is my data set enough to do multiple regression?   2. The t-test result revealed 3 of my variables are not significant. Do I need to do my regression again with the significant variables (or my first regression is enough to get conclusion)? My correlation matrix is as follow                     var 1   var 2    var 3   var 4   var 5     Y     var 1   1.0     0.0       0.0   -0.1    -0.3    -0.2     var 2   0.0     1.0       0.4    0.3    -0.4    -0.4     var 3   0.0     0.4       1.0    0.7    -0.7    -0.6     var 4  -0.1     0.3       0.7    1.0    -0.7    -0.9     var 5  -0.3    -0.4      -0.7   -0.7    1.0      0.8     Y      -0.2    -0.4      -0.6   -0.9    0.8      1.0      var 1 and var 2 are continues variables and var 3 to 5are categorical variables and y is my dependent variable . **It should be mentioned the important variable which has been considered in the literature as the most influential factor on my dependent variable is not also among my regression variables due to my data limitation. Does still make sense to do regression without this important variable?** here is my confidence interval                   Varibales   Regression Coefficient  Lower 95% C.L.  Upper 95% C.L.         Intercept   53.61                       38.46        68.76         var 1       -0.39                      -0.97         0.19         var 2       -0.01                      -0.03         0.01         var 3        5.28                      -2.28         12.84         var 4       -27.65                     -37.04       -18.26         **var 5      11.52                      0.90         22.15**


title: t-test for partially paired and partially unpaired data

text: An investigator wishes to produce a combined analysis of several datasets. In some datasets there are paired observations for treatment A and B. In others there are unpaired A and/or B data. I am looking for a reference for an adaptation of the t-test, or for a likelihood ratio test, for such partially paired data. I am willing to (for now) to assume normality with equal variance and that the population means for A are the same for each study (and likewise for B).


title: Tracking down the assumptions made by SciPy's ttest_ind() function

text: I'm trying to write my own Python code to compute t-statistics and p-values for one and two tailed independent t tests. I can use the normal approximation, but for the moment I am trying to just use the t-distribution. I've been unsuccessful in matching the results of SciPy's stats library on my test data. I could use a fresh pair of eyes to see if I'm just making a dumb mistake somewhere. Note, this isn't so much of a coding question as it is a "why isn't this computation yielding the right t-stat?" I give the code for completeness, but don't expect any software advice. Just help understanding why this isn't right. My code:               import numpy as np     import scipy.stats as st          def compute_t_stat(pop1,pop2):              num1 = pop1.shape[0]; num2 = pop2.shape[0];              # The formula for t-stat when population variances differ.         t_stat = (np.mean(pop1) - np.mean(pop2))/np.sqrt( np.var(pop1)/num1 + np.var(pop2)/num2 )              # ADDED: The Welch-Satterthwaite degrees of freedom.         df = ((np.var(pop1)/num1 + np.var(pop2)/num2)**(2.0))/(   (np.var(pop1)/num1)**(2.0)/(num1-1) +  (np.var(pop2)/num2)**(2.0)/(num2-1) )               # Am I computing this wrong?         # It should just come from the CDF like this, right?         # The extra parameter is the degrees of freedom.              one_tailed_p_value = 1.0 - st.t.cdf(t_stat,df)         two_tailed_p_value = 1.0 - ( st.t.cdf(np.abs(t_stat),df) - st.t.cdf(-np.abs(t_stat),df) )                       # Computing with SciPy's built-ins         # My results don't match theirs.         t_ind, p_ind = st.ttest_ind(pop1, pop2)              return t_stat, one_tailed_p_value, two_tailed_p_value, t_ind, p_ind      **Update:** After reading a bit more on the Welch's t-test, I saw that I should be using the Welch-Satterthwaite formula to calculate degrees of freedom. I updated the code above to reflect this. With the new degrees of freedom, I get a closer result. My two-sided p-value is off by about 0.008 from the SciPy version's... but this is still much too big an error so I must still be doing something incorrect (or SciPy distribution functions are very bad, but it's hard to believe they are only accurate to 2 decimal places). **Second update:** While continuing to try things, I thought maybe SciPy's version automatically computes the Normal approximation to the t-distribution when the degrees of freedom are high enough (roughly > 30). So I re-ran my code using the Normal distribution instead, and the computed results are actually further away from SciPy's than when I use the t-distribution.


title: Will AIC and R-square rank models similarly if the number of variables is equal?

text: This is a fairly basic question: If a number of regression models are being evaluated for fit to a dataset, and each model has 3 variables, will AIC and R-square rank the models similarly? I'm assuming the answer is yes, because the penalty aspect of AIC for added variables would be a non-issue. But, I'm unclear if this would always be the case. Thanks


title: R: How to "control" for another variable in Linear Mixed Effects Regression model?

text: Essentially, I have two collinear variables which could be seen as either random or as fixed effects, a dependent variable I'm fitting the model to, and a variable that's assuredly a random effect. **Dependent var:** Number of neuron spikes (`FiringRate`) in a specific region of mousebrain **Fixed effects:** **1)** `Time` at which data sample was taken (on a linear scale in days -- so day two would be 2, day 5 would be 5, and so on) **2)** The `Age` of the mouse in days (so there's definitely collinearity between this and the `Time` variable, but there are enough mice of different ages to make this worthwhile as a separate variable) **Random effect:** `Subject` \-- "Name" (ID number) of the mouse Essentially, I'm wondering if it would be appropriate to run two LMEs. In the first, I'd treat `Age` and `Subject` as random variables in order to control for the effects of `Age` (and thus the collinearity between `Age` and `Time`) and see if Time is a significant predictor of the # of spikes (dependent variable). In the second, I'd enter `Time` and `Subject` as random variables to see if `Age` was a significant predictor.               library(lme4)     a = lmer(FiringRate ~ Time + (1|Age) + (1|Subject))     b = lmer(FiringRate ~ Age + (1|Time) + (1|Subject))


title: Goodness of fit for 2D histograms

text: I have two sets of data representing stars parameters: an observed one and a modeled one. With these sets I create what is called a two-color-diagram (TCD). A sample can be seen here: ![histogramas](http://i.stack.imgur.com/qAOdF.png) **A** being the observed data and **B** the data extracted from the model (never mind the black lines, the dots represent the data) I have only one **A** diagram, but can produce as many different **B** diagrams as I want, and what I need is to keep the one that best fits **A**. So what I need is a reliable way to check the goodness of fit of diagram **B** (model) to diagram **A** (observed). Right now what I do is I create a 2D histogram or _grid_ (that's what I call it, maybe it has a more proper name) for each diagram by binning both axis (100 bins for each) Then I go through each cell of the grid and I find the **absolute difference** in counts between **A** and **B** for that particular cell. After having gone through all the cells, I sum the values for each cell and so I end up with a single positive parameter representing the goodness of fit ($gf$) between **A** and **B**. The closest to zero, the better the fit. Basically, this is what that parameter looks like: $gf = \sum_{ij} |a_{ij}-b_{ij}|$ ; where $a_{ij}$ is the number of stars in diagram **A** for that particular cell (determined by $ij$) and $b_{ij}$ is the number for **B**. This is what those $(a_{ij}-b{ij})$ differences in counts in each cell look like in the _grid_ I create (note that I'm not using absolute values of $(a_{ij}-b{ij})$ in this image but I **do** use them when calculating the $gf$ parameter): ![hess](http://i.stack.imgur.com/0ZAqy.png) The problem is that I've been advised that this might not be a good estimator, mainly because other than saying _this fit is better than this other one because the parameter is lower_ , I really can't say anything more. * * * ## **Important** : (thanks @PeterEllis for bringing this up) 1- Points in **B** are not related one-to-one with points in **A**. That's an important thing to keep in mind when searching for the best fit: the number of points in **A** and **B** is **not** necessarily the same and the _goodness of fit_ test should also account for this discrepancy and try to minimize it. 2- The number of points in every **B** data set (model output) I try to fit to **A** is **not** fixed. * * * I've seen the _**Chi-Squared**_ test used in some cases: $\sum_i (O_i-E_i)^2/E_i$; where $O_i$ is observed frequency (model) and $E_i$ is expected frequency (observation). but the problem is: what do I do if $E_i$ is zero? As you can see in the image above, if I create a grid of those diagrams in that range there will be **lots** of cells where $E_i$ is zero. Also, I've read some people recommend a _**log likelihood Poisson**_ test to be applied in cases like this where histograms are involved. If this is correct I'd **really** appreciate it if someone could instruct me on how to use that test to this particular case (remember, my knowledge of statistics is pretty poor, so please keep it as simple as you can :)


title: Estimating k in d=kv

text: This example was taken from Mathematical Statistics : A Unified Introduction (ISBN 9780387227696), page 58, under the section 'The Principle of Least Squares'. I think my problem has more to do with algebra but since this is a statistics book, I thought I should post here... Anyway, I need to estimate k base on this equation: ![alt text](http://i.stack.imgur.com/JQM6B.png) The goal is to reach here: ![alt text](http://i.stack.imgur.com/nqyzM.png) I managed to follow the example until I got stuck here: ![alt text](http://i.stack.imgur.com/1Xs42.png) The explanation to reach to the goal from above is described as follows: ![alt text](http://i.stack.imgur.com/j2UyX.png) I know how to get k but I totally don't understand how the 'middle term' was eliminated. Please help. I will provide more details if needed.


title: How to decide on the correct number of clusters?

text: We find the cluster centers and assign points into different bins in k-means clustering  which is a very simple algorithm and is found almost in every machine learning material on the net. But the lacking and more important part in my opinion is the choice of correct k. What is the best value for it? And, what is meant by saying _best_? I use MATLAB for scientific computing and looking at silhouette plots is given as a way to decide on the correct number of clusters as discussed here. But, I would be more interested in Bayesian approaches. Any suggestions are appreciated.


title: Does a sparse training set adversely affect an SVM?

text: I'm trying to classify messages into different categories using an SVM. I've compiled a list of desirable words/symbols from the training set. For each vector, which represents a message, I set the corresponding row to `1` if the word is present: > "corpus" is: [mary, little, lamb, star, twinkle] > > first message: "mary had a little lamb" -> [1 1 1 0 0] > > second message: "twinkle little star" -> [0 1 0 1 1] I think this is fairly common setup with SVM, but my question is, with thousands of words in the set, what if there's only 1-2 words per message that actually show up? Is the linear dependence of my set of training vectors going to adversely affect the ability of the algorithm to converge?


title: Why do we assume that the error is normally distributed?

text: I wonder why do we use the Gaussian assumption when modelling the error. In Stanford's ML course, Prof. Ng describes it basically in two manners:   1. It is mathematically convenient. (It's related to Least Squares fitting and easy to solve with pseudoinverse)   2. Due to the Central Limit Theorem, we may assume that there are lots of underlying facts affecting the process and the sum of these individual errors will tend to behave like in a zero mean normal distribution. In practice, it seems to be so. I'm interested in the second part actually. The Central Limit Theorem works for iid samples as far as I know, but we can not guarantee the underlying samples to be iid. Do you have any ideas about the Gaussian assumption of the error?


title: Conditional probability estimation in multinomial regression?

text: I have been analyzing data where the response consists in ordered categories (6) of damage to structures, with a number of damage predictors, thus I have been using multinomial regression. However, the model can be improved. I know that the probability of a given structure to be at damage level 6 will have an influence on the lower damage levels. This is because, damage level 6 corresponds to collapse of the structure which components become projectiles that most of the time affect the other structures. My problem is, I am not sure how to include this information in the model. I thought about using something like, p(ds=DSi ¦ DSn=1) for each previous damage level? DSn being the last damage state, "collapse". But how to implement this in a GLM? Or maybe there is another way to do this all together. Some pointers would be greatly appreciated!


title: Meta-analysis comparing difference between pre vs post treatment group vs pre vs post control/placebo group

text: I am hoping to obtain ES of 20 studies, each having pre vs post treatment compared to pre vs post placebo (i.e., is there a way to obtain ES through difference in means between the two groups and combined all 20 studies after that to get the total ES). I have data which include means for both groups pre and post, standard deviation of all 4 and number of each group. The program below is perfect to use: http://gemini.gmu.edu/cebcp/EffectSizeCalculator/d/mean-gain-scores-pre- post-r.html The only problem is that I would need to get correlation coefficient. Can I use this other program to assess "r" between pre and post treatment rather than "r" for treatment vs control? http://gemini.gmu.edu/cebcp/EffectSizeCalculator/r/means-and-standard- deviations.html Thanks for the feedback.


title: Probit regression in R using panel data - Do we need to normalize the predictors?

text: I am trying to run a probit regression using panel data by first computing the log likelihood and then using the `optim()` function to optimize. I have a couple of questions on this and would appreciate tips and suggestions from everyone.   1. Scale of predictor variables: The predictor variables that I have vary in scale significantly. Some of the variables are in tens while some are in '000s. When I am trying to run the optim function, it frequently gives the error message: vmmin is not finite. So I have to be extremely careful in choosing the initial values. And some of these initial values are as small as 0.0001. Is this standard in probit? Or do we need to normalize the predictor variables before using probit?   2. I have a few variables that are computed within the estimation code, i.e., some of the parameters that are being estimated are used to calculate these inbuilt covariates. When I run the full model with all such inbuilt parameters, the optim function gives me a singular hessian matrix. Again, I am not sure what does singularity of the hessian implies. When I add the inbuilt covariates one at a time, the hessian is not singular.  Did anyone encounter such problems before. Please advise.


title: Centering when using splines in R

text: I am having trouble understanding why centering seems to only work with simple linear models and not with splines for example. I am using centering to report the estimated group differences at different $x$, but also statistical values (ignoring multiple comparisons for the moment).               set.seed(1)          # simulate data     N <- 10     x <- rep(seq(0.2,1,0.2),N)     group <- factor(rep(c('I','II'),each=length(x)/N))     y <- -x^2 + 2*x*as.numeric(group) + rnorm(length(x),mean=0,sd=0.1)     d <- data.frame(group,x,y)          # fit a linear model with x-group interaction     l <- lm(y~x*group,data=d)     d$lmfit <- fitted(l)     coef(l)['groupII'] # group difference at x==0     #     groupII      #  -0.1097071           library(ggplot2)     ggplot(d,aes(x,y,colour=group)) + geom_point() + geom_line(aes(x,lmfit,colour=group))      The plot confirms the reported small group difference `groupII` of 0.05 at $x=0$ if we were to extrapolate back to 0. Now let us centre the data at $x=1$ and estimate the group difference there.               # center data at x==1 and refit     l <- lm(y~I(x-1)*group,data=d)     coef(l)['groupII'] # group difference at x==1     #   groupII      #  2.08525       In agreement with the plot the difference is about 2. Now let us fit a spline model.               # fit data with splines     library(splines)     l <- lm(y~ns(x,2)*group,data=d)     d$lmsplinefit <- fitted(l)         coef(l)['groupII'] # group difference at x==0.2         #     groupII          #  0.2987893          # compare to: d$lmsplinefit[6] - d$lmsplinefit[1]          ggplot(d,aes(x,y,colour=group)) + geom_point() + geom_line(aes(x,lmsplinefit,colour=group))      Interestingly, the spline fit reports the group difference at the first $x$, i.e. $x=0.2$. If we try to centre at $x=1$ we get the same result, i.e. the difference at $x=0.2$.               l <- lm(y~ns(I(x-1),2)*group,data=d)     coef(l)['groupII']     # same result as un-centered data, i.e. 0.2987893      Why is that? And is there a way to show the group difference at a different $x$? Btw, centering $x$ manually before the model fit does not make a difference.


title: Accuracy outcomes puzzle when using logit for prediction

text: I'm building a predictive logit model **for a rare event** using 5-fold cross- validation. I have three years of data. I build a model with a training set of years 1 and 2, then I test on the test set with year 1 predictors and year 2 outcome. Then I test the model "out-of-sample" using year 2 predictors and year 3 outcome. My process for building the model is as follows:   1. In each cv-run I build a model on the training data.  `fit = glm(y ~ ., family=binomial("logit"), data = train)`   2. Then I save a matrix of `n x k` 0-1 predictions, where `n` is the number of observations and `k` is the number of discrimination thresholds I am testing.    3. Calculate accuracy statistics with the test set for all `k` vectors of 0-1 predictions. There are two puzzles: **First, why does logit produce such high accuracy (with high specificity and low sensitivity) predicting on the test set in year 2?** **Second, why do the accuracy statistics reverse when predicting year 3? That is, when predicting year 3 from year 2 data, I get very high sensitivity, but low accuracy and specificity -- the opposite of predictions for the test set.** I have noticed that across _almost_ all the cv iterations, across all discrimination thresholds, accuracy is always greater than 50% when predicting the test set. **Is this an artifact of logit or something else (e.g., my rare- event data, my modeling choices)?** This is different from any other approach -- lasso, OLS, random forests -- for the same rare-event prediction, where very low discrimination thresholds produce high sensitivity, but low accuracy and low specificity. By _discrimination threshold_ , I mean the predicted probability of event occurrence (logit fitted value) above which I predict that the rare event will occur in the next time period. With all the other approaches, I do not observe such major switching of the relative magnitudes of the three accuracy statistics. That is, with any other method, the accuracy statistics for predicting year 3 are pretty close to accuracy statistics from predicting onto the test set. Only with logit are those two sets of accuracy statistics dramatically different. Unfortunately, I cannot upload the data. Complete R code is below. _But note that I do not think this is a code question, but rather a matter of the underlying statistical approach, which is why I posted here at CV rather than SO._               #### LOGIT          ## Section     #### Run the main model     ####    Logit     # - split sample into k folds     # -- * start loop     # - fit model k times,  each time on k-1 folds     # - for k models we have n predicted probabilities (n x k matrix of pred. prob.)     # - average pred. prob. of occurence across the k models,  save n long vector of pred. prob.     # - use the current discrimination threshold to classify pred. prob. as 0, 1 predications     # - compare these predictions to the observed y values     # -- * end loop,  run this loop for many discrimination thresholds               ## Define parameters of model     ROC_seq = seq(from = 0.01,  to = 0.3,  by = 0.01)     CV.runs = 200          # Define holder matrix for results          output = matrix(NA,  nrow = CV.runs,  ncol = 6)     colnames(output) = c("accuracy.rate",  "true.pos.rate",  "true.neg.rate",  "FPTP",  "FNTP",  "threshold")          # Save 2012 predictions to matrix below for each model and we'll average them later     # This is 200 rows,  where each row has a set of predictions for all nrow(data_years2_3) communities          predictions_2012 =  matrix(NA,  nrow = CV.runs,  ncol = nrow(data_years2_3))          for(j in 1:CV.runs){     set.seed(06510+j)          #       set.seed(06511)     data_randomization <- data_years1_2[sample(1:nrow(data_years1_2),  nrow(data_years1_2),  replace=FALSE), ]     dim = nrow(data_years1_2)     splitID =       dim - round(dim/5)     train <- data_randomization[1:splitID, ]     colnames(train) = colnames(data_years1_2)     test  <- data_randomization[(splitID+1):dim, ]     colnames(test) = colnames(data_years1_2)     # Remove to unit ID column     train = data.frame(train[, -2])     test = data.frame(test[, -2])               #- run k-fold CV using glmnet nfolds = 10 by default     fit = glm(ytest ~ .,  family=binomial("logit"),  data = train)          pred2010  =  fit[["fitted.values"]]          #define matrix     # each column is the predictions given one threshold,  across units     # each row is a vector of predictions across thresholds,  given one unit           ROC_results = matrix(NA,  nrow = length(ROC_seq),  ncol = 6)     colnames(ROC_results) = c("accuracy.rate",  "true.pos.rate",  "true.neg.rate",  "FPTP",  "FNTP",  "threshold")          for(m in 1:length(ROC_seq)){         pred2010_01 = ifelse(pred2010> ROC_seq[m], 1, 0)             errors = table(train[, 1],  pred2010_01)         ROC_results[m, 1] <- (errors[1, 1]+errors[2, 2])/sum(errors) #accuracy.rate         ROC_results[m, 2] <- (errors[2, 2]/(errors[2, 2]+errors[2, 1])) #true.pos.rate         ROC_results[m, 3] <- (errors[1, 1]/(errors[1, 1]+errors[1, 2])) #true.neg.rate         ROC_results[m, 4] <- (errors[1, 2]/errors[2, 2]) #FPTP         ROC_results[m, 5] <- (errors[2, 1]/errors[2, 2]) #FNTP         ROC_results[m, 6] <- ROC_seq[m] #threshold         }          ROC_results = subset(ROC_results,  ROC_results[, "accuracy.rate"]>0.5)          threshold = ROC_results[1, "threshold"]          ##Now apply the optimal threshold from this run to the test data          pred2010_test=predict(fit,    newdata = data.frame(test))          pred2010_test_01 = ifelse(pred2010_test> threshold, 1, 0)                 errors = table(test[, 1],  pred2010_test_01)         output[j, 1] <- (errors[1, 1]+errors[2, 2])/sum(errors) #accuracy.rate         output[j, 2] <- (errors[2, 2]/(errors[2, 2]+errors[2, 1])) #true.pos.rate         output[j, 3] <- (errors[1, 1]/(errors[1, 1]+errors[1, 2])) #true.neg.rate         output[j, 4] <- (errors[1, 2]/errors[2, 2]) #FPTP         output[j, 5] <- (errors[2, 1]/errors[2, 2]) #FNTP         output[j, 6] <- threshold #threshold          ### Save predictions_2012 for this model          predictions_2012[j, ] = predict.glm(fit,    newdata = data.frame(data_years2_3),  type = "response")          }          output_mean = apply(output, 2, mean)     output_sd = apply(output, 2, sd)          #Now use na.omit to calculate FPTP and FNTP     infID = is.infinite(output)      output[infID] <- NA      output_mean2 = apply(na.omit(output), 2, mean)     output_sd2 = apply(na.omit(output), 2, sd)          ### Calculate true out of sample accuracy     pred2012 = apply(predictions_2012, 2, mean)     pred2012_SD = apply(predictions_2012, 2, sd)     pred2012_01 = ifelse(pred2012> output_mean["threshold"], 1, 0)            performance2012 = c()         errors = table(data_years2_3[, 1],  pred2012_01)         performance2012[1] <- (errors[1, 1]+errors[2, 2])/sum(errors) #accuracy.rate         performance2012[2] <- (errors[2, 2]/(errors[2, 2]+errors[2, 1])) #true.pos.rate         performance2012[3] <- (errors[1, 1]/(errors[1, 1]+errors[1, 2])) #true.neg.rate         performance2012[4] <- (errors[1, 2]/errors[2, 2]) #FPTP         performance2012[5] <- (errors[2, 1]/errors[2, 2]) #FNTP          names(performance2012) = c("accuracy.rate",  "true.pos.rate",  "true.neg.rate",  "FPTP",  "FNTP")


title: 10% dip in February for metrics that count?

text: February usually has 28 days, unlike it's neighbors January and March which have 31. It seems that most countable things will exhibit, on average, a 10% dip in February for the missing 3 days. This isn't a matter of pure statistics, but it's something I just noticed. Is this known? Is it ever adjusted for charts?


title: Theoretical expected value and variance

text: Let $X$ be a random variable having expected value $\mu$ and variance $\sigma^2$. Find the Expected Value and Variance of $Y = \frac{X−\mu}{\sigma}$. I would like to show some progress I've made so far, but honestly I've been thinking about this problem for the past few days but just have no idea where to start. Any hint or insight on a starting point would be much appreciated. Thanks!


title: Tips when regressing growth rates

text: I have 20 months of Year over Year growth rates for a X independent variable and Y dependent variable. The correlation between these two variables is 0.72. I would like to predict Y using X for future periods. The challenge is, if for a given month Y grew a lot the same month last year, it faces a tougher comparable. My X variable is not aware that Y faces a tougher comparable (the correlations are not perfect, so they are operating on slightly different scales). Is there an elegant solution to minimize out of sample error? Should I try including the growth rate from last year's month as another independent variable? Or something else?


title: Correction for multiple testing with unequal sample sizes

text: I perform several thousand tests on a genomic dataset (Mann Whitney U test). However, there is great variability in the sample size for each test (per test the number of measurements (independent) varies between 4 and 20 with outliers up to 75). I understand that tests with more measurements (samples) have more power and are more likely to survive multiple testing (I use FDR at the moment). (i.e. the resulting list will be biased towards regions with larger sample size) However, a reviewer pointed out to me that, apart from this bias, the p-values are "not comparable" and that standard adjustment procedures like FDR can not be used at all. I understand the bias, but if I accept that (and describe it in the interpretation of the results as I do extensively), I do not see a problem.


title: Information theoretic alternative to hypergeometric test

text: Is there any information-theoretic alternative to the hypergeometric test (or Fisher's exact test)? In other word, is it possible to calculate AIC values for data which one would classically analyse with these tests?


title: How to check for normal distribution using Excel for performing a t-test?

text: I want to know **how to check a data set for normality in Excel, just to verify that the requirements for using a t-test are being met**. For the right tail, is it appropriate to just calculate a mean and standard deviation, add 1, 2 & 3 standard deviations from the mean to create a range then compare that to the normal 68/95/99.7 for the standard normal distribution after using the norm.dist function in excel to test each standard deviation value. Or is there a better way to test for normality?


title: How to approximate error on Chi-Squared Fit when bin counts are zero

text: I am using a galaxy image simulator that provides a 2D histogram that has the number of photons per pixel (bin) $N$. I am currently using a least-squares residual: $\sum_{bin}(f_{data}-f_{model})^2$ for the fitter (The Levenberg-Marquardt Algorithm) to approximate certain parameters of this "image". This least-squares fit does not divide by the uncertainty in each bin. Now when I decided to approximate the uncertainty in each bin by $\sqrt{N}$ such as: $\sum_{bin}\frac{(f_{data}-f_{model})^2}{\sqrt{N}}$ then I have to divide by zero in bins that contain zero entries. The way I remedy this is that I make the error equal to 1. Using this form, I receive negative biases for a subset of my parameters. How can I remedy this? And, as a comment pointed out, the second equation is not a chi-squared value?


title: Number of Predictors and Classification Algorithm

text: In general, is it better to include more predictors in algorithms such as SVMs and random forests compared to logistic regression? It seems that when we add more predictors to logistic regression, the improvement is not significant. However, when we add more predictors to SVMs and random forests, the improvement is significant. Why is this?


title: logistic regression for modelling

text: ![enter image description here](http://i.stack.imgur.com/OQpth.png) I have these data plotted above. The explanatory variable represents intensity levels of ground shaking at different locations in an earthquake, and the response variable represents amounts of compensation awarded to damaged buildings. I want to fit a model to these data so that how much compensation is likely to be awarded can be estimated given an intensity level. Can I use logistic regression here? How should I go about building the model for these data? Thanks in advance.


title: regression with non-independent data

text: I will be performing regression on subjects total scores from 2 player games (prisoners dilemma) that they will be playing. I am aware that including both players score from a game will cause problems due to non-independence. Is there a way that to deal with this apart from randomly picking one subject from each game for the analysis (and so losing half the data). Is there a way to introduce this into the model instead, perhaps as a random effect?


title: Probit analysis using 13 cases and checking normal distribution

text: Can I use probit if my dependent variable has data that looks something like this, 1 case is equal to 0, and 12 cases are equal to 1? Thanks


title: Technique to remove factor structure from panel data

text: Let's say I have panel data (i.e. time-series data and cross-sectional). Factors drive the bulk of the data as well as some noise. Naturally I can express this factor structure using a covariance matrix - or a set of eigenvectors and corresponding eigenvalues. I would like to create a new panel dataset that has some particular (or all) factors removed (i.e. de-factored). What's a good way to do this in R? My preference is to use factor analysis rather than PCA.


title: Citation for permutation test

text: I am using a permutation test to verify the statistical significance of an improvement in multi-class classification performance. What citation would be appropriate? Would it be fine to cite "Permutation tests: ..." by Philip Good, 2000?


title: Question about the frequency definition of probability

text: In the book, it is written that the relative frequency or the frequency ratio gradually tends to become more or less constant as N becomes larger and larger. It is also mentioned that this is an empirical fact. My question is: why should the frequency ratio tend to a constant in the long run ? I mean, there isn't any _convergence proof_ working after all (is there ?).


title: Why t distribution in confidence level for the mean?

text: This is a basic question, but when I was asked about it I only could give a weak answer. That's why I am asking it here. If we want to calculate the confidence level for the mean there is the formula: $(\bar{x}-t_{n-1}s_{\bar{x}},\bar{x}+t_{n-1}s_{\bar{x}})$ I know that the t distribution comes in, when the variance is unkown and has to be estimated. So the assumption that sample means are normally distributed is violated. The sample means follow the t distribution. But why? For me there is a point missing or at least it is not intuitive to me.


title: gain score vs. residualized change score when analysing correlates of change (NO group comparison)

text: Im interested in the relationship between change in so-called dynamic risk factors (e.g. low self-control) and recidivism, as well as the relationship between prison climate and change in dynamic risk factors. Thus, change is the IV in the first question, the DV in the second question. Dynamic risk is measured at two points of time. In both situations, im wondering if i should use a simple difference score or a residualized change score (or using the posttest score while controlling for pretest score, respectively). I've read a bunch of articles on this topic (e.g. this or this or this) as well as some other posts in this forum (e.g. here or here), but i'm still lacking a convincing solution to my specific problem (mainly, because most literature on this topic focusses on adjusting for baseline differences in group comparisons). Following most literature, i would generally prefer using the change score, but there is one important pitfall: dynamic risk at pretest is correlated with recidivism (otherwise, it wouldnt be a risk factor), and those with high scores at pretest show larger changes than those with low scores at pretest (to which extent this amounts to regression to the mean would be another nice topic). Therefore, those exhibiting more prosocial change are also recidivating more, something which is not really meaningful to interpret. On the other hand, controlling for the baseline seems somewhat wrong to me because the fact that people with high (negative) baseline scores are more amenable to change is not simply a statistical problem, but of theoretical and contentual interest. So, what to use and why? Thanks for your help!


title: Expected error of best possible linear fit?

text: In my textbook, there is a statement mentioned on the topic of linear regression/machine learning, without a proof or rigorious justification, which is simply quoted as, > Consider a noisy target, $ y = (w^{*})^T \textbf{x} + \epsilon $, for > generating the data, where $\epsilon$ is a noise term with zero mean and > $\sigma^2$ variance, independently generated for every example > $(\textbf{x},y)$. **The expected error of the best possible linear fit to > this target is thus** $\sigma^2$. Here my concerns are, How do we know what the best possible linear fit is ? Assuming in some way we estimated best linear fit, is it $y = (w^{*})^T \textbf{x}$ ? If so why can not other parameter $w_2^{*}$ be the best fitting parameter ?


title: Can I use the Stationary wavlet transforming directly with GARCH models

text: I am trying to estimate and forecast the garch models directly with wavelet transforming but I am having an estimation error in the estimation step "constraints: 'Boundary Constraints Active; Errors may be Inaccurate" The of C in mean equation is in minus and the always the optimization stopped after the 10 iterations. I guess that I am doing something wrong. Although I am following the steps as they are exactly mentioned in the wavelet toolbox in Matlab. I am using the SWT with Db5 and level 4. Daily stock index data which transformed into continues return series. The length of sample is long enough for garch, 4336 observations and extended via the extension tool in Matlab( with special and carful extension for SWT in wavelet toolbox). Also I checked the data after doing the extension just to make sure that the series shape and characteristics have not significantly changed.


title: Parametrization of Gamma and Negative Binomial in R

text: I have some Poisson data {${y_1,...,y_n}$} and a Gamma prior, and I wish to construct a predictive posterior distribution. As I understand, if my Gamma hyperparameters are $\alpha$ (the prior number of occurrences) and $\beta$ (the prior number of observations), the predictive posterior is: ${y\space|\space y_1...y_n,\alpha,\beta}\space\tilde\space NegBin({\sum\limits_{i=1}^n y_i+\alpha}, \frac{1}{1+\beta+n})$ However, because of the multiple parametrizations of both Gamma and negative binomial, my attempts to actually implement this in R were so far futile. Am I right that in the above formula $\alpha$ is `shape` and $\beta$ is `rate` in `dgamma()`? Am I right then that $\sum\limits_{i=1}^n y_i+\alpha$ is `size` in `dnbinom()`? And what is `prob` then? I thought it was $\frac{1}{1+\beta+n}$, but the resulting distribution doesn't seem to make sense.


title: Method to compare 2 models of different years (not by a dummy variable)

text: I'm a student in economics having trouble with a research in SPSS (BINARY LOGISTIC REGRESSION). I'm using 2 different datsasets with different samples, same variables, each dataset represents a different year (economic recession year VS non-economic recession year)(SPSS). I want to investigate the effect of my IV on my DepV in each dataset + I want to compare each individual IV for the economic recession year to see if there was a significance greater/ smaller effect as compared to the non-economic recession year. **PROBLEM** : The database where i got my data from (NSSBF) uses in the 2 years different definitions in some variables. I have to keep this in mind when interpreting the results. (example: 1993: > 50% in hands of family = family firm , 2003: 100% in hands of family = family firm). If I copy the 2nd dataset into the first and make a dummy variable for year (to do an interaction effect with my IV), I get the two different definitions mixed in the IV?. So then my interaction is correct but not my single effect of family firms (without comparison of the years). Can someone please help me with this? You have to excuse me for my limited knowledge in statistics, but I'm trying for days to figure this one out. Thank you!


title: Post-hoc analysis on Friedman's test

text: I am using a Friedman's test with post-hoc analysis. I am using the function `friedman.test.with.post.hoc` available for R software. Lets say we have 6 wines and we want an ordered list of the wines according to the judges' rates (best wine, ... ,worst wine). So, we perform the Friedman test. Suppose $p < \alpha$, we reject the null hypothesis. In other words, there is a significant difference between the wines given the judges' rates. In order to identify which wines are different, we run the post-hoc analysis. By running the post-hoc test, we have that there is only significant difference between the pairs: Wine A/Wine C and Wine A/Wine D. Given that, how do you conclude the order of the wines? Should I use any other tools? Thanks.


title: How to implement a fractional polynomial transformation in R for logistic regression

text: I'm working on a data set modeling road kills (0 = random point, 1 = road kill) as a function of a number of habitat variables. Following Hosmer and Lemeshow, I've examined each continuous predictor variable for linearity, and a couple appear nonlinear. I'd like to try a fractional polynomial transformation for each, also following Hosmer and Lemeshow, and have looked at the R package mfp, but I'm having trouble coming up with (and understanding) the R code that will correctly transform the variable. Can anyone suggest R code that would help me accomplish the concepts on p. 101 - 102 of Hosmer and Lemeshow's Applied Logistic Regression (2000). Thanks!


title: How can I calculate leaves and nodes of a C 4.5 decision tree?

text: I have a given amount of attributes, for example 1024, and now want to calculate the amount of leaves and nodes C 4.5 produces (approximately). Is there something like that for binary trees - a formula for example?


title: Is there an R package for constrained regression that has both a formula and predict interface similar to that of lm() and glm()?

text: Is there an R package for constrained regression, that has both a formula interface and a `predict` method? Two existing options are pcls from mgcv and nnls from nnls, but neither one appears to have a formula interface or a predict method.


title: Compare shapes with Kolmogorov-Smirnov test

text: I have two time series, both 100 points length: I would like to perform the goodness of fit to know if these two series have same or similar shape or if they differs significantly. I have read about the K/S test for 1 sample and it seems it compares the serie with a normal or gaussian distribution. I do not want to compare my two series with two standard gaussian shapes, I only want to know if they are similar or different. I have found this class for the K/S test http://www.jsc.nildram.co.uk/api/index.html Here are my questions: 1) is the K/S two sample test the best one for non parametric test with series.length()< 150 points? 2) are there other statistical test to compare shapes? 3) the jsc Java class has some different methods, which one should I have to use to return how much my series differs? Thanks all.


title: Some of my predictors are on very different scales - do I need to transform them before fitting a linear regression model?

text: I would like to run linear regression over a multi-dimensional data set. There exist differences among different dimensions in terms of their magnitude of order. For instance, dimension 1 generally has value range of [0, 1], and dimension 2 has value range of [0, 1000]. Do I need to do any transformation to ensure the data ranges for different dimensions are in the same scale? If it has to, are there any guidance for this kind of transformation?


title: How to calculate the rowMeans with some single rows in data?

text: I tried to group a set of elements and to calculate the row means If my list has more than one element, it works fine:               tapply(colnames(myMA), c(1,1,1,2,2,2,3,3,4,4), list)         myMAmean <- sapply(myList, function(x) rowMeans(myMA[,x]))      However, in my data some of the rows are unique:               myList <- tapply(colnames(myMA), c(1,1,1,2,2,2,3,3,4,5), list)     myMAmean <- sapply(myList, function(x) rowMeans(myMA[,x]))      Notice the 4 & 5 it makes them as unique so if I run this it says               "Error in rowMeans(myMA[, x]) :        'x' must be an array of at least two dimensions"      I don't want to run a loop. Is there any solutions for this?


title: Why do I get different heteroscedasticity robust standard errors in R when using the plm package?

text: I am now writing my bachelors thesis and I have come across some difficulties. I am about to do some panel regressions with time and entity fixed effects and I would therefore like to use the plm package. But when I do add fixed effects and want to have heteroscedasticity robust standard errors they seem to be incorrect. Does anyone know why the HC standard errors differ? Here is my code               # Load data     load(file="panel")     attach(panel)          # Load packages     library(lmtest)     library(plm)               # Create two models. The lm.model is a linear model and as the     # LAND variable is a factor variable representing countries     # (Land = Country in Swedish) this model will have entity fixed     # effects. In the plm.model the plm package is used and     # individual effects and within model is turned on (which is     # the same as entity fixed effects)     lm.model<-lm(NETTOSPARANDE ~ EURO + LAND, data=panel)     plm.model<-plm(NETTOSPARANDE ~ EURO, index=c("LAND","ÅR"), effect="individual", model="within", data=panel)          # When looking at the coefficents without heteroscadisity robust     # standard errors they are identical. They do also have the same     # value in stata.     coeftest(lm.model)[1:2,]     coeftest(plm.model)          # But when looking at the coefficents using heteroscadisity     # robust standard errors the lm.model and the plm.model produces     # different standard errors.     coeftest(lm.model, vcov.=vcovHC(lm.model, method="white2", type="HC1"))[1:2,]     coeftest(plm.model, vcov.=vcovHC(plm.model, method="white2", type="HC1"))      If you want to test the data it can be found here (the panel file) [1]: https://sourceforge.net/projects/emumoralhazard/files/ _R-data_ And here is my output               1> # Load data     1> load(file="panel")          1> attach(panel)          1> # Load packages     1> library(lmtest)     Loading required package: zoo          1> library(plm)     Loading required package: kinship     Loading required package: survival     Loading required package: splines     Loading required package: nlme     Loading required package: lattice     [1] "kinship is loaded"     Loading required package: Formula     Loading required package: MASS     Loading required package: sandwich          1> # Create two models. The lm.model is a linear model and as the     1> # LAND variabel is a factor variable representing countries     1> # (Land = Country in swedish) this model will have entity fixed     1> # effects. In the plm.model the plm package is used and     1> # individual effects and within model is turned on (which is     1> # the same as entity fixed effects)     1> lm.model<-lm(NETTOSPARANDE ~ EURO + LAND, data=panel)          1> plm.model<-plm(NETTOSPARANDE ~ EURO, index=c("LAND","ÅR"), effect="individual", model="within", data=panel)          1> # When looking at the coefficients without heteroscedasticity robust     1> # standard errors they are identical. They do also have the same     1> # value in Stata.     1> coeftest(lm.model)[1:2,]                  Estimate Std. Error   t value     Pr(>|t|)     (Intercept) -3.731024  0.7731778 -4.825570 1.726921e-06     EURO1        2.187170  0.4076720  5.365024 1.112984e-07          1> coeftest(plm.model)          t test of coefficients:                Estimate Std. Error t value  Pr(>|t|)         EURO1  2.18717    0.40767   5.365 1.113e-07 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1                1> # But when looking at the coefficients using heteroscedasticity      1> # robust standard errors the lm.model and the plm.model produces     1> # different standard errors.     1> coeftest(lm.model, vcov.=vcovHC(lm.model, method="white2", type="HC1"))[1:2,]                  Estimate Std. Error    t value     Pr(>|t|)     (Intercept) -3.731024  0.3551280 -10.506138 5.102122e-24     EURO1        2.187170  0.3386029   6.459395 2.009894e-10          1> coeftest(plm.model, vcov.=vcovHC(plm.model, method="white2", type="HC1"))          t test of coefficients:                Estimate Std. Error t value  Pr(>|t|)         EURO1  2.18717    0.33849  6.4615 1.983e-10 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


title: How to calculate probability for a multivariate/continuous event

text: I'm brushing up my stats for an AI course that is going to show me the Gaussian mixture model. My prof said I'm supposed to know how to "calculate the probability of a multivariate and continuous event X under the normal distribution". I have forgot this multivariate density function. Can someone please remind/show me how to do this? Thank you.


title: What is the conjugate of the noncentral hypergeometric distribution?

text: I want to find the conjugate prior distribution for the Fisher's noncentral hypergeometric distribution? Basically I want to integrate the parameters out of the distributions to arrive at the Bayesian Likelihood of observing some data with the noncentral hypergeometric distribution.


title: How to obtain aggregate estimate of a factor with various levels in General Linear Mixed Models

text: I run general linear mixed model with factor 'Substitutability' with seven levels in SPSS. It gave six level estimates. I want to calculate aggregate estimate for 'substitutability'. How to get? Please help.


title: Why some people test regression-like model assumptions on their raw data and other people test them on the residual?

text: I am a Phd Student in experimental psychology and I try hard to improve my skills and knowledge about how to analyze my data. Until my 5th year in Psychology, I thought that the regression-like models (e.g., ANOVA) assume the following things:   * normality of the data   * variance homogeneity for the data and so on My undergraduate courses lead me to believe that the assumptions were about the data. However in my 5th year, some of my instructors underlined the fact that the assumptions are about the error (estimated by the residuals) and not the raw data. Recently I was talking about the assumptions question with some of my colleagues who also admited that they discovered the importance of checking the assumptions on the residual only in their last years of university. If I understand well, the regression-like models make assumptions on the error. Thus it makes sense to check the assumptions on the residuals. If so, why some people check the assumptions on the raw data? Is it because such checking procedure approximate what we would obtain by checking the residual? I would be very intersting in a disccussion about this issue with some people who have more accurate knowledge than my colleagues and I. I thank you in advance for your answers.


title: Why does differencing time-series introduce negative autocorrelation

text: For example, this mentioned here: link. I also saw this in my data. I wonder - does anyone know a good reference where this is explained and justified more rigorously with some math and for some more-or-less wide class of processes? Here are some plots: Series with positive autocorrelation before differencing ![Series with positive autocorrelation before differencing](http://i.stack.imgur.com/TFOk3.gif) Autocorrelation before differencing ![Autocorrelation before differencing](http://i.stack.imgur.com/eRfd3.gif) Series after differencing ![Series after differencing](http://i.stack.imgur.com/lIctm.gif) Autocorrelation after differencing ![Autocorrelation after differencing](http://i.stack.imgur.com/EYj6j.gif)


title: Detecting outliers in circular data?

text: I've a day of the year data about number of event occurence in different sites: ![enter image description here](http://i.stack.imgur.com/MM62W.png) A day of the year is circular data. I know that a usual detecting of outliers, for example by boxplot is no use here: ![enter image description here](http://i.stack.imgur.com/S2qDe.png) How can I detect outliers in this kind of situations?


title: Correct way to compare two (very) different regression models?

text: I'm working with some piecewise linear regression models, and I'd like to compare their predictions with those produced by multiple (weighted) linear regression models. Both models describe the same physical system, but have very different parameterizations of the independent variable. The two different parameterizations are such that the averages of the independent variables are much different--that is, call $x_1$ the independent variables under the first parameterization, and $x_2$ the ind. vars. under the second parameterization. Generally, I have $\mathbb{E}x_1 \gg \mathbb{E}x_2$. This (in turn) means that the model coefficients may be very different. Additionally, it is sometimes the case that the piecewise linear regression may have a segment with slope = 0 and intercept = 0, which would seem to cause a problem for a statistic like CVRMSE. The _best_ way I can think to compare these two models is to just use a training and a test set, but then, I'm not sure what statistic I should compute to (algorithmically) say "this one is better". Is there a better way to discriminate between these two models in an a priori way?


title: Regression with complex causal structure

text: I have already run a whole bunch of OLS and found the following Regress P= B1*L+e_1, found b1<0 Regress X= B2*L+e_2, found b2>0 Regress X= B3*P+e_3, found b3<0 I want to build a case with complex causal structure, in which there are two channels that L positively affects X. 1.Reducing P increases X 2.L positively affects X on its own 3.In addition, L positively affects X through the effect that increasing L reduces P, and we know reducing P increases X, thus this is the second channel that L positively affects X. Two questions. A. What OLS or other type of regression equations[I guess it takes more than 1 equation] should I run to prove my case, ie, to prove (if not possible to prove, strongly suggest) this whole system of causal strucuture? B. What regression diagnostics should I use in SAS to show my logic behind the causal strcuture (L increases X through 2 channels )are ok? Thanks.


title: How to show unbiased estimator of combination of bernoulli and normal variables?

text: $X_1, X_2, \ldots, X_n$ is a random sample from $\mathrm{Bernoulli}(\theta)$, $\epsilon_1, \epsilon_2, \ldots, \epsilon_n$ are independent $\mathcal N(0, \sigma^2)$, independent of $X_i$. Define $Y_i = \theta X_i + \epsilon_i$, for $i = 1, 2, \ldots,n$. Define estimating function $$ \psi[\theta;(X,Y)] = \sum_{i =1}^n(Y_i - \theta X_i) \>. $$ a) Show that $\psi[\theta;(X,Y)]$ is unbiased estimator of $\theta$ if $\mathbb E(\psi[\theta;(X,Y)])=0$. b) Find the estimator $\hat\theta$ such that $\psi[\hat\theta;(X,Y)] = 0$. Is $\hat\theta$ unbiased?


title: Fitting an ordinal response variable

text: Suppose $Y$ is an ordinal response variable with 4 levels (i.e. $Y$ can take on the values $1,2,3$ or $4$). We want to model $Y$ with $X_1, X_2$ using a proportional odds model. If $Y$ takes on the values of $1,2,3$ and $4$ an equal number of times (e.g. $Y = 1$ 4 times, $Y = 2$ 4 times, $Y = 3$ 4 times, and $Y= 4$ 4 times), will the estimates of some of the coefficients be 0?


title: Best terminology to describe a time-series problem

text: I have a problem that I would like to investigate regarding time-series data, but as I'm inexperienced in the field of statistics, I am unsure of the best terminology to describe my problem (so I can try to investigate further myself). _Any suggestions on approaches are of course welcome!_ **The problem can be summarised as:** We receive several sets of deliveries on a periodic basis (although the timing can be imperfect and the periodicity varies between sets). For a particlar set of deliveries, sometimes a delivery has been missed. Say the deliveries are typically quarterly, then ignoring the slight inaccuracy in delivery dates, then you would get something like this:   * **X** = delivered item   * **0** = expected delivery missed   * N = no expected delivery **X** N N **X** N N **0** N N **X** N N **0** N N **X** N N etc I'm interested in how this can be analysed automatically, so given that one won't know the period in advance, I'd like to find a way to **a)** figure out the period, then **b)** see which expected deliveries are missing on that basis. There is an added complication that in some cases deliveries are made on a cycle that isn't continuous, e.g. in a repeating three month period, they are delivered roughly on the first and second month and then not on the third. Where I'm stumbling right away is that I don't know the formal terminology that one would use in describing this type of problem. Is this best described a survival analysis problem with repeating events? Or just a discrete time- series problem? Is there a special name for these sorts of binary event type time-series? (and also where the time measurements (i.e. the deliveries) do not occur at the same time - my understanding is that usually time-series data is gathered at the same interval. Anyway, I'm keen to help myself if at all possible but if there are suggestions on the best approach (e.g. spectral analysis approaches?) I'd be grateful for those too. Do ask for further clarification if I've missed something crucial or should give more detail. Many thanks, Neil


title: How to report the ratio of two sets of experimental results?

text: I have measured the time taken to solve a problem by algorithm $X$ and by algorithm $Y$. It takes a quite long time, so I have only 10 data for each algorithm: $$ X : ( x_1, x_2, \dots , x_{10}) \\\ Y : ( y_1, y_2, \dots , y_{10}) $$ **EDIT:** The problem they are solving is **randomized**. I generated 10 instances of the problem using 10 different random seeds. The 10 computational times correspond to these 10 problem instances. In this sense, **the data are paired**. The change of a seed does not change the difficulty of the problem very much. **END OF EDIT** I have computed the ratio of the averages: $$ avg = \frac{\sum_{k=1}^{10} x_k }{\sum_{k=1}^{10} y_k } $$ This however, does not convey any information on how **precise** the ratio is. One possible way is to estimate the **standard deviation**. According to this answer, the average of iid random variables is **asymptotically normal** and therefore the ratio has **asymptotical Cauchy distribution** , whose standard deviation is _infinite_. This does not satisfy me, especially since I have only 10 data. Then according to this answer I should approximate the standard deviation using **Taylor series**. This answer looks better, but it still doesn't feel right. The distribution of a ratio is intuitively highly **assymetrical around 1.** (you have only the interval $(0; 1)$ to capture the fact that algorithm $X$ is faster, but the entire $(1 ; \infty)$ to capture the fact that $Y$ is faster). So even a well estimated standard deviation can be of little use. It would be better to provide some sort of **confidence interval**. For example: the ratio is 1,5 with an assymetrical confidence interval of (1,3 ; 2,8). But I have no idea how to estimate this since I do not know the distribution of my data. **EDIT2:** Here are my data:               X       Y     111536  160134     111165  164850     112494  165844     115959  166409     121296  161755     119948  167781     119172  168666     117330  169766     116661  166518     129311  169884      **EDIT3:** To answer the question (in comment) of **D L Dahly** **why not just report that one algorithm is faster in all instances** For the brevity of the question I did not mention that I have actually 84 sets of data described in this question. 2 problems x 6 dimensions of the problem x 7 possible sizes of the problem. In some instances **X** is faster, in some **Y** is faster and in some instances the results are inconclusive. I do not necessarily need confidence intervals or standard deviations. I just want to provide the reader with **something richer than just averages**. The reader should have a sense of how well the average represents the experimental results.


title: Repeated measures one way ANOVA inconsistent with paired t- test result?

text: I have a repeated measures one way ANOVA with a significant interaction, but no significant main effects. I ran a paired samples t-test to further my analysis and found a significant difference between two conditions in my data. Is it okay to report both the paired samples t- test result and the ANOVA result together is the results section or can I only go with reporting the ANOVA interaction? It is my understanding that the paired samples t-test is more 'sensitive' to differences between conditions then the ANOVA, right?


title: Detecting outliers using correlogram

text: If there is an outlier in a time series, how does its correlogram behave? Is it possible to find outliers using a correlogram? **EDIT** I have such a Time series:               ts <- c(1,2,3,1,2,2,30,40,3,2,4,1,3,2,3,1,1,2)      its correlogram looks like this: ![enter image description here](http://i.stack.imgur.com/wGJ4S.png) I have an outlier in my time series, Is it possible to detect these outlier from correlogram?


title: Is there a method for testing follow-up individual factor and item invariance using the lava an or semTools packages in R?

text: I've found the "measurementInvariance" command from the semTools package EXTREMELY helpful recently, but now I am wondering if there is a way to conduct follow-up factor-level and item-level comparisons using a similar command? This would be particularly helpful as the "measurementInvariance" command automatically uses the robust estimator fit indices to compare groups (including calculating the chisq difference test), while I believe the alternative method of testing in lavaan (anova) defaults to using the "ML" estimates. For example:               HS.model <- 'visual  =~ x1 + x2 + x3     textual =~ x4 + x5 + x6     speed   =~ x7 + x8 + x9'     fit <- cfa(HS.model, data=HolzingerSwineford1939, estimator="MLR", std.lv=TRUE, group="sex")          measurementInvariance(HS.model, data=HolzingerSwineford1939, group="sex", estimator="MLR")      Once I determine there is variance (in factor loadings for example), I may want to compare each individual factor (and perhaps each item) to determine where the variability is found. I do know how to conduct this test (I think), but it uses the "ML" estimates to calculate the chisq difference test, instead of the robust fit measures that "measurementInvariance" uses. Here is how I currently would test for loading invariance on a single factor (e.g. "visual"), but would receive the "ML" fit measures in the anova test of chisq diff:               HS.model <- 'visual  =~ x1 + x2 + x3     textual =~ x4 + x5 + x6     speed   =~ x7 + x8 + x9'     fit <- cfa(HS.model, data=HolzingerSwineford1939, estimator="MLR", std.lv=TRUE, group="sex")          HS.model.b <- 'visual  =~ label("a", "a")*x1 + label("b", "b")*x2 + label("c", "c")*x3     textual =~ x4 + x5 + x6     speed   =~ x7 + x8 + x9'          fit.b <- cfa(HS.model.b, data=HolzingerSwineford1939, estimator="MLR", std.lv=TRUE, group="sex")     anova(fit.b, fit, SB.classic=TRUE)


title: Approximation of the variance of the first order statistic (min) of normal random variates

text: I'm looking for a closed form approximation of the variance of the minimum order statistic for normal random variates. Can anyone point me to a reference, or an approximation? I've seen the post about the approximate expected value of the min order statistic (Approximate order statistics for normal random variables), but I'm looking for the variance. Thanks!


title: How to extract "MLR" fit measures generated by the Lavaan package of R

text: I am estimating some Confirmatory Factor Analysis (CFA) models using the Lavaan package and I hoping to extract the fitMeasures to export out of the model to a spreadsheet. This is easily done using the `fitMeasures` command, but it automatically defaults to exporting the "ML" fit measures, even if the model was run using "MLR". Is there a way to extract these `fitMeasures` from the robust model? Here is my example script:               ## The famous Holzinger and Swineford (1939) example     HS.model <- 'visual  =~ x1 + x2 + x3     textual =~ x4 + x5 + x6     speed   =~ x7 + x8 + x9'     fit <- cfa(HS.model, data=HolzingerSwineford1939, estimator="MLR", std.lv=TRUE)     fitMeasures(fit, fit.measures = c("chisq", "cfi", "rmsea", "aic"))      Running the summary reveals that `fitMeasures` produces "ML" fit measures, not Robust fit measures:               summary(fit, standardized=TRUE, fit.measures=TRUE, rsquare=TRUE)      Is there a way I can obtain the Robust fit measures? Thanks in-advance for your help! * * * **I think I found the solution to this (below):** I used the following in my example above to extract the robust (scaled) fit measures:               fits<-inspect(fit, "fitMeasures")     fits[c(5:7,24,44,34)]      I hope it helps someone out there and please let me know if there is a better way to accomplish this!


title: Comparing population proportions vs hypothesis testing

text: The author of a famous online course takes a problem where 642 of 1000 polled men and 591 of 1000 polled women responded _yes_. The question is if men are more likely to agree or not, using 95% confidence interval. He starts by computing the difference of the means, $\bar p_1 - \bar p_2 = .642 -. 591 = 0.051$ and its standard deviation (i.e. the standard error of the .051), which is a sum of variances $\sigma = \sigma(\bar p_1) + \sigma(\bar p_2) = \sqrt{1/1000 (p_1(1-p_1) + p_2(1-p_2))} = 0.021715$. This difference of means is normally distributed and, according to z-table, 1.96 standard deviations contain 95% of such mean differences. We thus can say that the difference is 0.051 ± 1.96 * 0.021715 = 0.051 ± .042562 or there is 95% confidence that $p_1-p_2$ lies between .008 and .094. Since zero is excluded from the region, the 95% confidence interval is purely positive and we are sure that men are more likely to say _yes_ than women. The problem is that after this conclusion, author conducts something similar yet different: he tests the hypothesis that $\bar p_1 = \bar p_2$! That is, he takes the population average, $p = (\bar p_1 + \bar p_2)/2$, computes its error by summing variances for men and women. BTW, I wonder why claiming that there are 2000 individuals averaged he uses formula $\sqrt{2p(1-p) \over 1000} = 0.0217$ instead of $\sqrt{p(1-p) \over 2000}$? The difference is 2x. Anyway, he takes the z-score = $0.051/0.0217=2.35$ deviations afterwards, which is outside the 95% confidence interval (of 1.96 deviations) and, therefore, it is unlikely that means, $\bar p_1$ and $\bar p_2$, are equal. I do not understand: what is the difference between the methods and which of the tests should I use?


title: Alignment and comparison of two unimodal and one uniformly distributed datasets

text: This question is similar to the following question: Normalize 3 irregulary distributed datasets and make their datapoints statistically relevant to each other describes similar problem, but is more basic, without the lame solutions I had given there, therefore, I think more appropriate for discussion. I am sorry. I am really stuck here. New to stats. I will simplify as much as I can explaining what I want to do. Have N learn-cases each one with 3 parameters: A, B and C. Their values all fit within the (-1, 1) interval, but are differently distributed along this interval. A and B have unimodal distribution. A has a more sharper peak and is situated at the left of the interval, B has smoother peak, and situated at the right side of the interval, but are generally non-symmetrical to each other. The C parameter has an uniform distribution. Let's say I have ONLY ONE sample learn-case also having the A, B and C properties with some values in the interval (-1, 1), I wish to check all my dataset learncases for the best similarity with this sample learn-case. **The questions** :   1. How would you do that? What comparison method ? Excluding neural networks please.   2. Before comparison do I have to make some normalization / alignment of A related to B, or to C as they are not proportional to one another?


title: How to properly present McNemar's Test Statistics

text: I have read an article and I wonder if that is a proper way to present a McNemar's Test statistics. The study try to measure the change in knowledge before and after a lesson. They performed McNemar's test to calculate the p-value, and tabulate it. I don't feel comfortable when reading the table because: McNemar's test just shows if the difference in the count of the discordant pairs (Pre-T + Post-F) OR (Pre-F + Post-T) are statistically significant. It does not give any clue on whether the question improve or worsen attendees' knowledge. But the table seems to claim that the education program caused an increased number of attendees improving their knowledge, support by McNemar's test, which I don't feel comfortable about that statement. I wonder if the count of the discordant pairs are included in the table (number of Pre-T + Post-F and number of Pre-F + Post-T), will give a different picture. My question is: By aggregating the counts in the 2x2 table into Pre-correct count and Post- correct and then present it with the McNemar's test statistics, is it appropriate? I tried to find one similar article in Lancet or NEMJ but I can't. If it is not appropriate, how can it be presented in a more proper way? Thanks.


title: XY correlation, time-efficiently

text: We have 25m XY pairs to be correlated within R. Data is bank financial data and smooth over Time ; observations for x and y are 32 quarters each. Testing 25m rships exhaustively will take forever ; this task is easily over- engineered. We’ll use the best of XY relationships to predict. (fn1) We’re swamped by choice in R packages (forecast, fastVAR, tseries, nlts etc.). All seem hard to compare. Can someone help us with which package(s) is/are most apt ? That is, which can time-efficiently test for correlation given properties of the data ? We’ll save lots of time with fine advice. (And yes, to conserve time, we’ll employ foreach & doSNOW.) Steve fn1 : X and Y are series both 50 quarters in length. Y is offset 8 quarters forward, such that X (periods 9 to 50) is compared to Y (1 to 42). Best fit (x actual vs x fitted using y actual) is determined. Method identified is used with Y (43 to 50) to predict X for the final 8 quarters.


title: Qualitative variable coding in regression leads to "singularities"

text: I have an independent variable called "quality"; this variable has 3 modalities of response (bad quality; medium quality; high quality). I want to introduce this independent variable into my multiple linear regression. When I have a binary independent variable (dummy variable, I can code `0` / `1`) it is easy to introduce it into a multiple linear regression model. But with 3 modalities of response, I have tried to code this variable like this :               Bad quality      Medium quality      High quality               0                1                  0          1                0                  0          0                0                  1          0                1                  0      But there is a problem when I try to do my multiple linear regression: the modality `Medium quality` gives me `NA`:               Coefficients: (1 not defined because of singularities)       How can I code this variable "quality" with 3 modalities? Do I have to create a variable as a factor (`factor` in `R`) but then can I introduce this factor in a multiple linear regression?


title: Binary logistic regression - Categorical function leaving out a category?

text: I'm still learning logistic regression, so hopefully my question makes sense. I have 10 independent variables and one dependent. The dependent and 3 of the independent variables are dichotomous (with no=0 and yes=1). The other independent variables all have a comparison of 0 and 3 options. When I go to analyze - regression - binary logistic, I put in all my variables and then click the categorical function in the top right corner. I highlight the 7 variables that are not dichotomous, and make sure they are set to compare to the first category. (I'm using SPSS 21 if that helps) After the model runs, everything is fine (all the variables are comparing to the 0), but I noticed that regardless of how many variables I have in the model, or the order they appear in the model, the last categorical variable will always be missing a category (although it is comparing the included categories to 0). Any ideas for why this is happening? Should I just create dummy variables for these rather than using the categorical function?


title: Bayesians: slaves of the likelihood function?

text: In his book "All of Statistics", Prof. Larry Wasserman presents the following Example (11.10, page 188). Suppose that we have a density $f$ such that $f(x)=c\,g(x)$, where $g$ is a **known** (nonnegative, integrable) function, and the normalization constant $c>0$ is **unknown**. We are interested in those cases where we can't compute $c=1/\int g(x)\,dx$. For example, it may be the case that $f$ is a pdf over a very high-dimensional sample space. It is well known that there are simulation techniques that allow us to sample from $f$, even though $c$ is unknown. Hence, the puzzle is: How could we estimate $c$ from such a sample? Prof. Wasserman describes the following Bayesian solution: let $\pi$ be some prior for $c$. The likelihood is $$ L_x(c) = \prod_{i=1}^n f(x_i) = \prod_{i=1}^n \left(c\,g(x_i)\right) = c^n \prod_{i=1}^n g(x_i) \propto c^n \, . $$ Therefore, the posterior $$ \pi(c\mid x) \propto c^n \pi(c) $$ does not depend on the sample values $x_1,\dots,x_n$. Hence, a Bayesian can't use the information contained in the sample to make inferences about $c$. Prof. Wasserman points out that "Bayesians are slaves of the likelihood function. When the likelihood goes awry, so will Bayesian inference". My question for my fellow stackers is: Regarding this particular example, what went wrong (if anything) with Bayesian methodology? P.S. As Prof. Wasserman kindly explained in his answer, the example is due to Ed George.


title: Can something show bias but not be significant?

text: I was reading a research paper and read this statement:               "A slight bias (although not significant) for     50/30 microRNAs to target genes with related      functions using miRanda (Figure 2B)"      However, if there is a bias: let's say one mean is greater than the other, but that it is not significant, can you say there is a bias at all? Couldn't it be (equally?) likely that the mean be lower than the other given that the two distributions aren't significantly indistinguishable? I understand that you may argue that if there was more data then that mean may eventually become significant, but the dataset is in the thousands and is likely to not to get much larger? Does the term 'slight bias' have any statistical meaning?


title: Can I define an $R^2$-like measure in this way when predicting with exotic loss functions?

text: It is common in regression to see $R^2$ formulated as follows: $$R^2\equiv 1 - {SS_{\rm err}\over SS_{\rm tot}},$$ where $SS_\text{err}=\sum_i (y_i - f_i)^2$ and $SS_\text{tot}=\sum_i (y_i-\bar{y})^2$. This exact formulation seems to make the most sense if the goal of the regression is indeed to minimize the residual sum of squares. However, sometimes other, more exotic loss functions are more appropriate for a variety of reasons, and some predictive modeling algorithms might be tried to minimize this loss function. I would like to define my $R^2$-like measure in the following way: $$R^2= 1 - \frac{\sum_il(y_i, f_i)}{\sum_il(y_i, \bar{y})},$$ where $l(y_i, f_i)$ denotes the "loss" occurred in some sense if $f_i$ is predicted but the actual answer is $y_i$. In the case of OLS regression, $l(y_i, f_i) = (y_i - f_i)^2$, but this allows for other predictive models based on least absolute deviations or even more exotic functions based upon the problem domain. Is this valid to do? Is anything gained with this more general formulation of $R^2$? Is it still valid to call it $R^2$? Does this add any useful interpretations, or is it more likely to just be confusing? Should the traditional $R^2$ be reported as well?


title: How to calculate the combined standard error of two models?

text: Let's say you have a model that predicts annual GDP growth. It has a standard error of 2%. You have a second model that predicts annual stock returns. It has a standard error of 6%. What is the standard error when you estimate stock returns based on estimates of GDP growth?


title: Can one have two random variables, perfectly correlated, but with different variances (as percent of their mean)?

text: I'm thinking of the difference between two random variables, e.g. the spread between two stock prices.


title: Factoring simple Markov network

text: Let $X$ by the joint distribution of the random variables $A$, $B$, $C$, and $D$. Let $(A \perp B) \mid (C, D)$ and $(C \perp D) \mid (A, B)$. ![Four Node Markov Network](http://i.stack.imgur.com/yrGho.png) I understand that this distribution should factor over the four pairwise cliques. $\Pr(X) = \frac{1}{Z}\phi_1(A,D)\phi_2(A,C)\phi_3(C,B)\phi_4(D,B)$ However, I'd like to see how, starting with $\Pr(X)$, we could factor the distribution into appropriate functions. I've been knocking my head against this most of the day, so any hints appreciated.


title: Set of data and averaging/standard deviation

text: I have a set of data (positive real numbers < 1) in five categories. My aim is to show that the data in the last category is bigger than the other categories for a range of examples (the data set), but the last category isn't necessarily more important than the others. I thought the best thing is to average each category and then plot the graph of it which works. But I feel this is too simple and is it good if I look at the standard deviation of each category to show that there isn't too much difference and the average-taking is valid? Also, I was going to use a scatter plot, is this OK? What would you recommend?


title: Transformation methods to make skewed distributions close to each other

text: I have a group of samples and observations in each sample constitute a right skewed distribution. So this group of samples constitute a group of skewed distributions. However, these skewed distributions don't correspond to the same distribution curve, or the distribution curves of different samples deviate from each other. My goal is to make these distributions close to each other through transforming the original data. I plan to try two methods to do the transformation. The first method is mean-based method in which each original value is divided by the mean value of the distribution. The second method is z-score. I call the two methods normalization methods under the goal of making different distributions more close to each other. After the transformation I can compare the effect of the two methods.   Before doing the transformation I want to read some related literatures about the normalization effect of mean-based method and z-score for skewed distributions. But after my thorough search in google I can't find any related papers. My question is has anyone seen related papers about this topic and can you give me some hints on how to find them? Welcome any comments on this problem too.


title: How to make panel VAR analysis in Eviews?

text: I want to perform panel VAR analysis in Eviews but I am not sure which is the correct option as there isn't any built in option in the software. Could you please advice what are the exact steps for performing panel VAR in Eviews?


title: Calculation of incidence rate for epidemiological study in hospital

text: I have kinda puzzled when I heard from the other about the calculation of incidence rate From Kenneth Rothman's Modern Epidemiology, the incidence rate is calculated as `number of cases in a fix period of time` divided by `person-time at risk at that period of time`, which mean if the patient has the disease in the middle of the month, only the first half of patient days will be included in the denominator, but not for the second half. Am I correct? Problem arises when I saw someone use the following formula for calculation: `number of cases in a fix period of time`/`total person-time at the period of time` If this is used, then the incidence rate will be underestimated, because the denominator now is bigger than it should be. But this formula seems to be used often, I wonder what this formula does. Thanks!