title: Combining and ranking standard deviations

text: 1. I am trying to take three independent indicators and combine them into a single composite measure.   2. For discussion, let's say that this is a composite measure for 'economic performance', and that the measure contains three independent indicators: worklessness, qualifications and unemployment. Let's also assume that there is a complete data set, with scores for each of the 50 local areas.   3. To compare and contrast those 50 local areas, we believe that the most robust way is to work out the standard deviation for each indicator (i.e. working out the standard deviation [s.d.] for each of the three indicators, for each of the 50 local areas).   4. The problem that we face, however, is when we come to combine the indicators together into a single composite measure. I understand from reading elsewhere (including on stats.stackexchange.com) that standard deviations cannot simply be added together; and that the correct method is to square each of the standard deviations individually to obtain the variance; add the variances together, and then divide this by the number of variances (i.e. divide by three in this case) and then taking the square root of this number.   5. If it was okay to add standard deviations together, there would be no problems, as we would simply be able to add standard deviations together, rank the resulting value, and quintile them.   6. However, if we follow (4.) above to combine standard deviations, the majority of the results now fall within one standard deviation of the mean. The smallest level of granularity is now 1 s.d., making it impossible to rank / quintile them, except assume that anything within 1 s.d. is in the mid quintile.   7. Is there any way I can justify simply adding together standard deviations, so that they can actually be ranked and quintiled, instead of following the 'proper' method?


title: Defining the "uniformity" of a dataset

text: I am working on a few algorithms where I have a list of $N$ samples. Currently I have plotted these into a histogram and have a view of how uniform the values are distributed within an interval, which is quite good as a visualization, although I need a comparable value of how uniform the dataset is, in order to measure how robust it is compared to my other algorithms. I have been looking at chi-squared test, but could not figure out how it would become helpful in my usecase? Sample from dataset:               8725     462     1492     972     9941     8235     8220     6949     1252      Code for importing data and applying chi-squared in `R`:               mydata = read.csv2("/opt/doc/stat/uniform_test_1.csv")     x <- sapply(mydata, as.numeric)     chisq.test(x)      Result: `X-squared = 1664769844, df = 999998, p-value < 2.2e-16`


title: Multiple imputation with the Amelia package

text: I have a general question about the Amelia package. I'm no mathematician or statistician, but I had to use R and impute and analyze some data, and Amelia showed results that fitted my expectations. I'll have to defend my choice soon, but I haven't totally grasped what Amelia does. I'm particularly interested in a simple as possible explanation in how Amelia imputation works. I've read that it uses a bootstrapping-based algorithm, but how does it chose the values? The data had mainly value >0 (chemical concentrations, water temperature and pH-value).


title: High Standard Errors for Coefficients imply Model is bad?

text: Suppose we have a regression model. If we get estimated of some of the coefficients and the standard errors are high, does this mean that the model is wrong/bad? How exactly do statistical packages choose regression models (in particular ordinal regression)?


title: Correlated time-series

text: First, I would like to apology if my vocabulary is not correct. I am not statistician (and not mother tongue English speaker either). So here is my problem : For n subjects I got 8 values (human pressure at 8 different times). I would like to know if my n subjects are presenting the "same pattern" ( i don't know which word to use ... are they all correlated?) Here is what I did (and would like to know if it is statistically correct) for each couple of subject (so sum(1:n-1) couples) I am computing their correlation. And with this sum(1:n-1) correlation values I am running a multiple comparison permutation t-test. If the p-value < 0.05 I consider that they are following "the same pattern". I know that there is a problem of dependence. But, I cant see how to do it in an other way. Any suggestion? Thanks a lot!


title: Segmentation using cluster analysis in SPSS

text: I am doing a segmentation project and am struggling with cluster analysis in spss right now. Could you please help me get this answered: How do I determine the _quality_ of the clustering in spss? In many articles/ tutorials I've read it's advisable to run a hierarchical clustering to determine the number of clusters based on agglomeration schedule and a dendogram - and then to do k-means clustering. Let's say the 1st step results are not clear and I am hesitant between 4 and 5 cluster-solutions. I can try both with k-means method - but how do I see which one's best? Same goes for any re-running of k-means clustering procedure, since every time the output is slightly different. Thanks a lot for any info! I would also be grateful for link to any good ready tutorials on cluster analysis in spss. What I've found so far is very random and limited... Articles on cluster analysis are not enough for me, because then I don't know how to run the tests mentioned in them in spss :-[ Cheers!


title: Large coefficients and std. errors

text: I run a fixed effect model with `Stata` and because my dependent variable is a large number (max of 12 million and min of - 4 million), I got large coefficients for all independent variables, although most variables are significant for example coefficients   * 122896.2   22478.689   * 21096.348   412.539 is it a problem thanks in advance


title: Updating probabilities given sample

text: Simplified version of our problem. Assume each person can be one of three categories.   * The true distribution of the entire population, as a whole, of being in the i'th category is (10%, 20%, 30%, 40%)   * Your prior belief of Person A of which category he/she belongs to is (25%, 35%, 15%, 25%)   * A sample is taken from the population. Person A belongs in the sample. The distribution of the sample, as a whole, is (5%, 25%, 20%, 50%) (Purpose chose non-uniform numbers to avoid any ambiguity) If we are given the proportion within each category (5%, 25%, 20%, 50%) we can update Person A's belief to (12.8%, 44.8%, 10.2%, 32%) Question now becomes: What if we weren't given the entire proportion within each category [C1,C2,C3,C4], but linear combinations instead? How would we update the belief if just given ([C1+C2],[C3+C4]) for (30%,70%) to continue? Or ([C1+C2],[C1+C3]), etc. (Example: Assume four categories are all combinations of Men/Women and Young/Old. Instead of giving the proportion within each four combination the output only gives total just Men or total just Young, or some combination thereof)


title: SEM: sample covariance can not be inverted

text: I am trying to fit a structural equations model, which looks like this: {model<- 'ambition=~NA*imprich+ipshabt+ipsuces+iprspot+0*ipeqopt+0*iphlppl altruism = ~0*imprich+0*ipshabt+0*ipsuces+0*iprspot+ipeqopt+iphlppl ambition ~~ 1*ambition altruism ~~ 1*altruism ipshabt ~~ iprspot imprich ~~ ipsuces ambition ~ altruism+agea+male+cntry_nl+cntry_be altruism ~ agea+male+cntry_nl+cntry_be' fit<-sem(model,data=q) summary(fit,fit.measures=TRUE,standardized=TRUE, modindices=FALSE)} But I keep getting the following error message: Warning message: In lav_samplestats_from_data(lavdata = lavdata, missing = lavoptions$missing, : lavaan WARNING: sample covariance can not be inverted I have already understood that is might be due to a too small sample size. My sample size is however 3532. Shouldn't this be sufficient? Can somebody tell me what to do please. Many thanks in advance, Pieter


title: How to use the Huber/White estimator of covariances in a generalized linear mixed model (glmmPQL) in R?

text: An analysis was implemented in SPSS 22 that uses the "Generalized Linear Mixed Models" feature of the program. Now I am looking for a way to port this to R. I use the glmmPQL() function of the MASS package and got the same results as in SPSS when using "model-based covariances". But, the data requires to use the "robust covariances" instead. In the SPSS manual I discovered what is the difference between both methods of estimating the covariances: > Covariance matrix. The model-based estimator is the negative of the > generalized inverse of the Hessian matrix. The robust (also called the > Huber/White/sandwich) estimator is a "corrected" model-based estimator that > provides a consistent estimate of the covariance, even when the > specification of the variance and link functions is incorrect. Source: SPSS > manual How can I use the Huber/White estimator in R?


title: Predicted probabilities from simulated betas and hypothetical data after conditional logit?

text: I'm working with conditional a conditional logit model to avoid bias that comes with FE logit models, when it comes to generating some hypothetical substantive effects, however, I run into trouble. My normal process (following methods employed by simcf, CLARIFY, etc packages):   1. Estimate model in R. `m<-glm(y~x, family="binomial")`   2. Use mvrnorm in R to get several thousand simulated coefficients: `betas <\- mvrnorm(10000,m$coefficients, vcov(m))`   3. Come up with some hypothetical values for the x's and calculate the predicted probability off the simulated betas: `x <\- .5 simmu <\- beats%*%x estimate <\- mean(1/(1+ exp(-(simmu)))` When it comes to conditional logit using a similar method will get you the unconditional prediction, but this is normally extremely high and not the quantity of most interest. What I want, I think, is the conditional prediction. Stata calculates it like this: ![enter image description here](http://i.stack.imgur.com/dc7Pw.gif) Where $i$ is the group we're conditioning on and T is the set of observations within Y. When calculating these probabilities from observed data, this seems simple enough. But how does one calculate these conditional probabilities for hypothetical scenarios?


title: Can I do two factor repeated measures ANOVA with differing numbers of measures per group?

text: I want to compare reaction time in four groups but I have three measures (pretest,post test and delayed post test) only in one of the groups. it was impossible to give post test and delayed post test in the other two groups and in the last group just possible to give pretest and post test, and only the first group received treatment. Now I want to see whether reaction times has decrees during time? Now since in three of the groups measurement has been down incompletely and no treatment has been given, can I use the previous measures for the lacking measures and use repeated measure ANOVA?


title: Repeated measure nested within individuals?

text: All, (I'm a beginner so my question could be really dumb...) So, I have a set of data that contains 30 schools' variable A and variable B. These two variables were measured repeatedly each year for 10 years. I'm interested in knowing whether these two variables are correlated, but was confused by the repeated measurement that is involved. I've only learned basic ANOVA and regression, so I searched online and found this "repeated measure nested within individuals..." Is this the correct method I should use? Any help is greatly appreciated! Thanks a lot!!


title: Expected coverage of a set of random samples

text: For a set A, I'm running 8 independent random samples, each with a probability of 1/8=12.5% and is without replacement. I know that the set formed by the union of these 8 samples will be of a size between 12.5% and 100%. Is there an expected value for the size of this union?


title: aov() error strata degrees of freedom

text: I have the following data frame:               df=structure(list(S = structure(c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,      3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L,      7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L, 10L,      10L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L,      14L, 14L, 14L, 14L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 1L,      1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 5L,      5L, 5L, 5L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 8L, 8L, 8L, 8L, 9L,      9L, 9L, 9L, 10L, 10L, 10L, 10L, 11L, 11L, 11L, 11L, 12L, 12L,      12L, 12L, 13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L, 15L, 15L, 15L,      15L, 16L, 16L, 16L, 16L), .Label = c("s1", "s2", "s3", "s4",      "s5", "s6", "s7", "s8", "s9", "s10", "s11", "s12", "s13", "s14",      "s15", "s16"), class = "factor"), A = structure(c(1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("a1",      "a2"), class = "factor"), B = structure(c(1L, 1L, 2L, 2L, 1L,      1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L,      1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L,      1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L,      1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L,      1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L,      1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L,      1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L,      1L, 2L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, 2L, 2L), .Label = c("b1",      "b2"), class = "factor"), C = structure(c(1L, 2L, 1L, 2L, 1L,      2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,      2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,      2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,      2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,      2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,      2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,      2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,      2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("c1",      "c2"), class = "factor"), y = c(1, 1, 1, 1, 1, 1, 0, 1, 1, 1,      0, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,      1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1,      1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1,      0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1,      0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,      1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0)), .Names = c("S", "A",      "B", "C", "y"), row.names = c(NA, 128L), class = "data.frame")      Using deviation coding ( `options(contrasts=c("contr.sum","contr.poly"))`) I have been attempting a simple mixed-design ANOVA in R with `aov()`. However, I am not getting the correct error strata, nor the correct degrees of freedom. In a very similar post, the issue was resolved by converting all independent variables into factors; in my data, all are already factors. There is also much written about how to specify error strata for `aov()`; none of the oft- conflicting guidelines have solved my problem. `A` is between-subjects, nests the id variable `S` and is crossed with crossed, within-subjects `B` and `C`. The data is ("are") completely balanced. My best guess for the R code that I need for my analysis is:               aov(y~A*B*C+Error(S),df)      which, as this post here and, e.g., this page, lead me, perhaps erroneously, to believe should be equivalent to:               aov(y~A*B*C+Error(S/(B*C)),df)      It is not:               > summary(aov(y~A*B*C+Error(S),df))          Error: S               Df Sum Sq Mean Sq F value Pr(>F)     Residuals 15  1.125   0.075                         Error: Within                Df Sum Sq Mean Sq F value   Pr(>F)         A           1 10.125  10.125  59.476 7.44e-12 ***     B           1  0.781   0.781   4.589   0.0345 *       C           1  0.031   0.031   0.184   0.6692         A:B         1  0.031   0.031   0.184   0.6692         A:C         1  0.281   0.281   1.652   0.2015         B:C         1  0.500   0.500   2.937   0.0895 .       A:B:C       1  1.125   1.125   6.608   0.0116 *       Residuals 105 17.875   0.170                          ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1               > summary(aov(y~A*B*C+Error(S/(B*C)),df))          Error: S               Df Sum Sq Mean Sq F value Pr(>F)     Residuals 15  1.125   0.075                         Error: S:B               Df Sum Sq Mean Sq F value Pr(>F)       B          1 0.7812  0.7812   5.952 0.0276 *     Residuals 15 1.9687  0.1312                      ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Error: S:C               Df Sum Sq Mean Sq F value Pr(>F)     C          1  0.031 0.03125   0.146  0.708     Residuals 15  3.219 0.21458                         Error: S:B:C               Df Sum Sq Mean Sq F value Pr(>F)     B:C        1   0.50  0.5000   2.308   0.15     Residuals 15   3.25  0.2167                         Error: Within               Df Sum Sq Mean Sq F value   Pr(>F)         A          1 10.125  10.125  64.371 4.47e-11 ***     A:B        1  0.031   0.031   0.199  0.65740         A:C        1  0.281   0.281   1.788  0.18621         A:B:C      1  1.125   1.125   7.152  0.00963 **      Residuals 60  9.437   0.157                          ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      To start, in neither case do I get `A` listed with its correct denominator $df$, $$ df_{S/A} = a(n-1) = 2(16-1) = 30 $$ Nor can the error strata for `B`, `C` and `B`x`C` ("`S:B`", "`S:C`" and "`S:B:C`", respectively) be correct, based on the $df$ output; both `B` and `C` are of two levels, so: $$ df_{BxS/A}, df_{CxS/C} = a(b-1)(n-1) = 2(2-1)(16-1) = 2(1)(15) = 30 $$ $$ df_{BxCxS/A} = (b-1)(c-1)(n-1)a = (2-1)(2-1)(16-1)2 = (1)(1)(15)2 = 30 $$ Assuredly either my code is wrong or my math is, or both. Please help.


title: What is the best source to start learning SAS?

text: I'm looking for some book recommendations, video lectures, tutorials or anything that can help me learn SAS.


title: Quantitative ,longitudinal, observational study:Population, sample and instrument

text: I am studying US public company financial performance for the last 2 decades. I am searching for evidence of "sustainable superior performance", thus to isolate those companies that display such characteristics. I will use the Compustat database. Questions: 1) Is the total population the universe of US public company with their past, present and future performance **or** the population is the subset of those displaying superior performance (subject of the study)? 2) Because all the companies will be observed to spot those with certain characteristics (i.e. superior performance), what is the sample? Can I consider the company present in the Compustat database (independently if they have the traits I am seeking or not) to be the sample? 3) What is the Instrument I need to ensure validity and reliability? -The Compustat database with its screening capabilities? -The criteria I will use to assess "sustained, superior performance"? (like sales growth above industry average and so on) - The statistical technique used to find such outlier? (like Cluster Analysis for instance) Thanks a lot for your help.


title: Log Likelihood and Score Vector of Conditional Grouped Continuous Model (CGCM) when the data is given

text: I am interested in the relationship between weight Y_i at mating and the number Z_i of lambs born in a flock of n=25 female sheep.Assume Z_i can be 0(no lambs born) 1 and 2, then I know that Y_i* will be linked to Z_i by the following threshold model:   1. **Z_i** = 0 if Y_i* < a    2. **Z_i** = 1 if a <= Y_i* < b   3. **Z_i** = 2 if Y_i* >= b With (Y_i* ,Y_i)**T distributed N_2 with mean=(0,mu)* _T and variance_ *T=(1 r ; r sigma**2) And I know that Y_i and Z_i are jointly distributed according to the **CGCM**. My question is if I know the data (weight,litter size) How can i find   * **log likelihood**   * **score vector** Then I need to use score vector in a built-in optimization function in R


title: Non Stationary series, VAR Model

text: Hi I am working with a multivariate time series data consisting of 1) Demand Data 2)Sales Data 3) Rainfall Data , all available from 2010-2013,at monthly level. Approach: I am trying to estimate the effect of rainfall on demand and sales of the product using VAR( Vector Autoregression) model. Questions : Is it correct to check individual series for stationarity ? What to do if there is no cointegration present among the series, i.e ( in My case: Demand series is non stationary, Rainfall and Sales series are stationary)   How to proceed with the model. Also if the eigen value test shows that VAR system is stable , does it imply that over all VAR is stationary and it does not matter whether individual time series are stationary or not, Should we go for differencing ? Thanks Nany


title: Testing significances in Likert type items

text: I have a questionnaire of 18 questions with Likert-type answers, going from 1-5. One of them was administered to parents, one to their children. I need to test 2 thing: a) within each group, are the answers between each combination of two questions statistically significant? I know this sounds a bit strange, but I need to see if any of the questions is more 'important' than the other. Do I need to use nonparametric tests and test each possible pair of questions, or could I create a correlation matrix (Spearman's, Kendall's) to do this? Or is PCA the way to go? b) I need to test if there is a significant difference on questions between groups? E.g. do parents and children generally give the same answers to Q1, Q2, etc. I'm doing this in R, and generally have an idea which functions to use for which test, but I really need advice on my methodology. Thanks!


title: Supervised learning : Hebb learning rule doubts

text: (A) In this book "Introduction to Neural Network using Matlab 6.0 - S. N. Sivanandam, S. N Deepa" http://books.google.ca/books?id=jJTN8RPgyXgC&printsec=frontcover#v=onepage&q=genetic&f=false, it is mentioned that Hebb is a supervised learning algorithm in Page21. Later, under Learning Techniques, Chapter 3, Section 3.3 Hebbian is categorized under Unsupervised learning. This is confusing !! (B) Again in this impelementation NN Hebb Learning, there is a target vector and the weight is updated by cross-multiplying the transpose of the input with the target (Did not quite understand). (C) In matlab               learnh calculates the weight change dW for a given neuron from the neuron's input P, output A, and learning rate LR according to the Hebb learning rule:          dw = lr*a*p'      This equation shows that Hebb is an unsupervised algorithm So, my questions are (Q1) Is Hebb supervised or unsupervised learning algorithm? (Q2) What is the correct equation? (Q3) Can somebody provide a simple code for Hebb learning? (Q4) What are the other supervised learning algorithm? (Q5) Does evolutionary algorithm for learning like Genetic algorithm, particle swarm fall under supervised? Thank you


title: Comparing Elo Ranking system with another ranking system

text: I'm a guy who is interested in computer science and algorithms. I've spent a bit of time recently developing my own ranking system based on PageRank. I'd like to compare it to the Elo system, in particular, the chess Elo ratings of the top 50 chess grandmasters. I've got 2 datasets. All games played between the top 50 in 2012 and 2013. I've produced my own set of rankings (using a modified PageRank algorithm) based on the 2012 dataset, and I would like to see how well this ranking predicts the results of the 2013 dataset. I would also like to know how well Elo predicts the 2013 dataset, to see which ranking is more accurate. I do not know how to do this in a meaningful way. One problem I see is that a chess game between two players, P and Q, can be either a win, loss or draw. How do I use rankings to create some sort of distribution for the expected result of a game between P and Q that also accounts for draws? And even if I have such a distribution (eg P wins = 50%, Q wins 20%, draw = 30%) and I have a game played between P and Q that I want to predict, how do I select from this distribution fairly? Any thoughts or insights would be helpful. If it helps, my ranking system will assign each player a value between 0-1, eg   * 0.05832816: Karjakin, Sergey   * 0.04895855: Grischuk, Alexander   * 0.04840682: Carlsen, Magnus   * 0.04751164: Ivanchuk, Vassily   * 0.04692380: Nakamura, Hikaru   * 0.04330322: Morozevich, Alexander   * 0.04281909: Aronian, Levon


title: Proof that CRF models and logistic models are convex functions

text: Where can I find a good proof that CRF based models and logistic regression based models are convex? Is there a general trick to test/prove that a model or objective function is convex?


title: Comparison of Slice sampling and Gibbs sampling

text: To me, the two are similar in the sense that slice sampling is just Gibbs sampling for the uniform distribution over the area under the plot of the density function. Is that right? I was wondering if someone can compare between slice sampling and Gibbs sampling. For example, in terms of rate of convergence, which one is better? If you can think about other aspects, please feel free to reply. Thanks and regards!


title: Can I say independent components of ICA are also linear independent?

text: In other words, does the independent components form a vector basis? EDIT:Let'me try to clarify the question If I calculate $n$ independent components for a given set of $n$ dimensional observations, are those components also linearly independent (property required to form a vector basis)?


title: 95% confidence interval for a given data set

text: I have a data set $\mathcal{S}$, which is univariate. I would like to discover the data points which can be thought of as outliers of $\mathcal{S}$. I use $\bar{\mathcal{S}}+2*\sigma(\mathcal{S})$ and $\bar{\mathcal{S}}-2*\sigma(\mathcal{S})$ to compute the upper and lower threshold in terms of 95% confidence interval. $\bar{\mathcal{S}}$ is the mean of $\mathcal{S}$ and $\sigma(\mathcal{S})$ is the standard deviation. I know it assumes the normal distribution. I would like to know are there any loopholes for this approach from the statistical analysis perspective. Moreover, are there any other better statistical approaches to give this kind of threshold.


title: Terminology - Comparing Z Scores

text: Let's say I have some sample of people with a measure of interest for each one, let's call it `v`. Let's say I measure `v` twice for the same set of people, one in January and one in December. Now I want to assess whether a person's `v` value improved from January to December. Let's assume that a high `v` value means they did well. I thought of the following analyses, and I think I have seen it somewhere, and I wanted to know what the name of this is. Consider taking each `v` value in January and giving each person a z-score. That is,               January_z_i = (January_v_i - mean(January_v))/sd(January_v))      where `January_z_i` is the z-score for a single person, `January_v_i` is the `v` value for a single person in January, and `January_v` represents the entire list of `v` values for January. Now for December, I go ahead and do the following               December_z_i = (December_v_i - mean(January_v))/sd(January_v)      Notice how I used the `mean` and `sd` of `January_v` for the December group. The logic is that, assuming the same sample space as measured in January, did that same person improve within that space if you consider their score in December. What is the technical name for this analysis. Sorry I didn't explain it more succinctly.


title: Relationship between $\text{Cov}(x_i^2, e_i^2)$, the asymptotic variance of b under homoscedasticity and heteroscedasticity?

text: I am trying to figure out the relationship between $\text{Cov}(x_i^2, e_i^2), V$ and $V_0$, where: $V=$ asymptotic variance of $\sqrt{n(\hat{\beta}-β)}$ under heteroskedasticity, and   $V_0=$ asymptotic variance of $\sqrt{n(\hat{\beta}-β)}$ under homoscedasticity Also, I am given that that the expected value of x squared, E(x_i^2), is equal to one. From this, I get that: $V = 1/N ∑〖(x_i^2 \cdot e_i^2)〗$   $V_0=σ^2$ I am certain that this is right. Now, using this I get that $\text{Cov}(x_i^2, e_i^2) = V - 1/N ∑〖(e_i^2)〗$ Now to get $V_0$ into this equation, by mathematical manipulation, I obtained $1/N ∑〖(e_i^2)〗= V_0$, and thus that $\text{Cov}(x_i^2, e_i^2) = V-V_0$ Is this correct? I am afraid that the substitution of $1/N ∑〖(e_i^2)〗$ for $V_0$ is only possible if I assume homoscedasticity. Is this worry justified? If so, how do I proceed? Grateful for a response. In case anyone has the energy to read this, here is how I derived the above (I apologize for the notation, but I do not know how to change it): Cov(x_i^2 〖,e〗_i^2 )=E(x_i^2 e_i^2 )-E(x_i^2 )E(e_i^2 )                                 =V-E(x_i^2 )E(e_i^2 )                            =V-E(e_i^2 ),  as E(x_i^2 )=1                            =V-1/N ∑〖(e_i^2)〗      Now, e_i=y_i-X_i'β ̂ = X_i'β+ε_i-X_i'β ̂=ε_i+X_i'β-X_i'β ̂, e_i^2=ε_i^2+(β-β ̂ )'X_iX_i'(β-β ̂ )+2ε_i X_i'(β-β ̂ ), σ ̂^2=1/N ∑(e_i^2 )=1/N ∑(ε_i^2)+(β-β ̂ )'1/N ∑(X_iX_i') (β-β ̂) +2 1/N ∑ε_iX_i'(β-β ̂ )                             (1)              (2)                               (3)      (2) and (3) converge in probability to zero, as it is assumed that beta-hat is a consistent estimator for beta. (1) converges in probability to σ^2. Therefore, 1/N ∑(e_i^2 ) converges in probability to σ^2 Thus, we have that Cov(x_i^2,e_i^2 )=E(x_i^2 e_i^2 )-E(x_i^2 )E(e_i^2 )=V-1/N ∑(e_i^2 ) =V-V0


title: Fixed effect Logit with R

text: I would like to perform a Fixed effect logit estimation in R. Can someone point out a package that can do the job? Note: For the time being I'm not really interested in the random effect. Update; Essentially I wonder if there is the plm package for a binary response model. Here is some documentation for the plm package: http://cran.r-project.org/web/packages/plm/index.html   http://dss.princeton.edu/training/Panel101.pdf


title: Wild Bootstrap F-test

text: Is anyone aware of code for a wild bootstrap F-test? Namely, testing joint significance of several coefficients where the standard errors of each coefficient are computed using a wild bootstrap. The following paper exposits the theory, but I wanted to see if anyone has implemented this in practice? Baltagi, B., Kao, C., Na, S. (2013) Testing for cross-sectional dependence in a panel factor model using the wild bootstrap F test _Statistical Papers_ , _54_ :1067–1094.


title: How to use one time series to predict another?

text: **Note:** There is another post with very similar title, but it is really on the author's specific problem with submarines. Mine is more general. I have asked a similar question before but turns out it wasn't a good question, so I am asking a modified one. (a practice recommended by the meta site I think) ## Question starts: I have time series $X_t$ and $Y_t$, $t=0,...,N$. I want to develop a model to use $X_t$ to predict $Y_t$. This is not a forecasting question, because in the remote future I won't observe $Y_t$. So I won't be using the history of $Y_t$ to predict its future. So I thought the easiest approach is to use linear regression `lm(Y~x)`, possibly with lagged covariates, except the residuals $\epsilon_t$ are very clearly auto-correlated. So i thought I could try Rob Hyndman's OLS with ARIMA residuals. But then this is where I got stuck. Say you have trained your model now and the residuals follow an AR(1). When you end up using it to predict at the remote future $t=N+M$, $M$ large, you can easily get the fixed part $X_{N+M}\hat{\beta}$, but   1. How do you get $\hat{\epsilon}_{N+M}$? Do you set it to be 0? Or do you do the crazy thing of forecasting it $M$ timestep from $t=N$?    2. Say you set it to 0, what about $\hat{\epsilon}_{N+M+1}$? Do you start using your AR(1) developed with a 0 for $\hat{\epsilon}_{N+M}$?   3. On that note, what if the residuals are ARIMA? Do you set the first difference to be 0?   4. Is there a better approach? How about `nlme:gls` with `corStruct=corARMA()` (although it seems to be doing the same thing)?  [edited: this is realy a general methodology question. But if you want something more concrete, say both $X_t$ and $Y_t$ are daily data over 4 years. So maybe seasonality is there.]


title: How many quarters of observations are suggested before and after an event to determine the events significance in a time series analysis?

text: I'm reviewing a multiple regression time series analysis, and I'm concerned about the (possibly insufficient) sample size. The goal of the analysis is to determine whether an event during the time period had an effect on the dependent variable. The event is dummy coded (0 before the event, 1 for every period thereafter). The data we have is on a quarterly basis. I'm concerned that the sample size is so small that a single observation could skew the results. The DV is a macroeconomic indicator. How many quarters of data are suggested before and after an event to determine statistical significance?


title: What is the best way to proof your time series has a deterministic linear trend

text: In case you suspect your time series has a linear trend, what is the best way to prove it? If you just regress it against time, you ignore the auto correlation of the time series so I assume that is a bit foolish. We learned to do unit root tests, and in those test you remove the auto correlation by including lags of the time series. So that problem is solved. But I am struggling a bit with the unit root test which includes a trend ( `ur.df(type="trend")` in R ) , because the intercept is also included. But this intercept is like a drift, doesn't it remove the significance of the trend? Why don't we test for unit root with only the presence of a trend and not the intercept?


title: Use PCA for pattern recognition

text: I have always used PCA as a preprocessing techniques. There are some people that use it for pattern recognition. Once I used pca to project my data in a 2 dimensional space and then I observed that the samples where divided in 2 groups. What other kind of analysis can be done using PCA? Or at least what are the most used analysis done with PCA? A good answer should have some small theoretical description and some examples. R code is highly appreciate.


title: Estimating parameters of Student's t-distribution

text: What are the maximum-likelihood estimators for the parameters of Student's t-distribution? Do they exist in closed form? A quick Google search didn't give me any results. Today I am interested in the univariate case, but probably I will have to extend the model to multiple dimensions. EDIT: I am actually mostly interested in the location and scale parameters. For now I can assume that the degrees of freedom parameter is fixed, and possibly use some numeric scheme to find the optimal value later.


title: Stratification vs. Adjusting in Survival Analysis

text: Suppose you want to compare survival of males and females within three groups. What is the difference between performing a log-rank test for each group (stratification) vs adjusting for it? Does the stratification allow for the survival curves in each group to differ?


title: Fitting an ANOVA with too many predictors in R. No p-values but how can I still interpret the model?

text: So I've fit an ANOVA model to 32 predictors and 32 data points. The summary output simply gives me the SS for each predictor and no F-stat or p-value. Is there any way to still make some sort of discussion as to the validity or effectiveness of the model?


title: Interpretation of R's output for binomial regression

text: I'm quite new on this with binomial data tests, but needed to do one and now I´m not sure how to interpret the outcome. The y-variable, the response variable, is binomial and the explanatory factors are continuous. This is what I got when summarizing the outcome:               glm(formula = leaves.presence ~ Area, family = binomial, data = n)          Deviance Residuals:      Min      1Q  Median      3Q     Max       -1.213  -1.044  -1.023   1.312   1.344            Coefficients:                             Estimate Std. Error z value Pr(>|z|)      (Intercept)           -0.3877697  0.0282178 -13.742  < 2e-16 ***     leaves.presence        0.0008166  0.0002472   3.303 0.000956 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      (Dispersion parameter for binomial family taken to be 1)          Null deviance: 16662  on 12237  degrees of freedom     Residual deviance: 16651  on 12236  degrees of freedom     (314 observations deleted due to missingness)     AIC: 16655     Number of Fisher Scoring iterations: 4      There's a number of things I don't get here, what does this really say:                                       Estimate Std. Error z value Pr(>|z|)      (Intercept)           -0.3877697  0.0282178 -13.742  < 2e-16 ***     leaves.presence        0.0008166  0.0002472   3.303 0.000956 ***      And what does AIC and Number of Fisher Scoring iterations mean?               > fit     Call:  glm(formula = Lövförekomst ~ Areal, family = binomial, data = n)          Coefficients:     (Intercept)        Areal       -0.3877697    0.0008166            Degrees of Freedom: 12237 Total (i.e. Null);  12236 Residual     (314 observations deleted due to missingness)     Null Deviance:      16660      Residual Deviance: 16650        AIC: 16650      And here what does this mean:               Coefficients:     (Intercept)        Areal       -0.3877697    0.0008166


title: Will this have a closed form expression: Conjugate priors

text: This is continuing the discussion about the posterior probability densities and conjugate priors. So, my posterior distribution over the parameters is defined as: $$ P(w, \phi, \lambda | y) \propto P(\phi) P(w|\lambda) P(\lambda) \prod_{i}P(y|w, \phi, \lambda) $$ Now, $\phi$ is the parameter that models the noise precision. Noise is assumed to be global zero mean and IID distributed. So, the mean is known and we have an unknown precision and the conjugate prior is given by a Gamma distribution. So, $P(\phi) = Ga(\phi; a, b)$. Now, the distribution over $w$ where both the mean and precision are unknown, we can write $P(w, \lambda)$ as $P(w|\lambda)P(\lambda)$. As described in Bishop's book, to achieve conjugacy, $P(w|\lambda)$ is modelled as a zero mean MVN whose precision is a function of $\lambda$ given by $(\lambda \Lambda)^{-1}$. Also, $P(\lambda)$ is modelled using a gamma distribution i.e. $P(\lambda) = Ga(\lambda; c, d)$. So, now my questions are about what the posterior distribution. 1: Due to the use of such conjugate priors, is the posterior over all the parameters $P(w, \phi, \lambda | y)$ easier to $\textit{sample}$ from or is there a closed form expression for it which obviates the need for sampling all together and I can estimate the sufficient statistics analytically? 2: Can I say something that this joint posterior will have some particular form or can I only say something about the marginal distributions i.e. $P(w)$ will be normally distributed posteriorly, $P(\phi)$ will be gamma distributed etc. What about the sufficient statistics of these marginals? Can they be computed analytically? What about dependence/coupling of the parameters in the posterior distribution? I am guessing since we have coupling of $\lambda$ and $w$ in our model, it should be coupled in the posterior distribution as well. However, to perform any sort of approximate inference on such a system, do we $\textit{have}$ to using mean field approximations and assume independence between the parameters in the posterior distribution?


title: Sample size too large for t-test?

text: All Drugs   CASE    CONTROL     AVG         11.07   7.88     ST DEV      3.58    2.95     # SUBJECTS  7399    3537          NM DRUGS             AVG         5.49    3     ST DEV      4.33    2.51     # SUBJECTS  6831    2507      Ho - Cases and controls are prescribed the same numbers of medication alpha =.05 I've got some data like those above and I wish to compare if the cases took significantly more drugs than the controls. These tables are two separate, but related data. When I ran a t test the respective t-stats were like 49 and 39. Should I be using a different test here?


title: Separating different frequencies from time series

text: I have a time series, which is produced from a sensor. This sensor detect sound in the real World. My aim is to separate sounds with different frequencies from this time series. With which technique can I do this job?


title: Repeat normality check after modifying group cut-off?

text: I have a set of data which consists of one independent variable (2 groups) and one dependent variable. I successfully checked for normality (each group separately) and conducted a t-test. Now I want to check for robustness by changing the cut-off between groups. Do I have to also do another normality check first with the different group cut-off?


title: Omitted Variable Bias, verification in Gretl

text: I am trying to verify the expression for Omitted Variable Bias (OVB) as given e.g. in Wooldridge: $\tilde{\beta_1} = \hat{\beta_1} + \hat{\beta_2} \cdot \tilde{\delta_1}$, where $\tilde{\delta_1}$ is the estimated slope of the regression of $x_2$ on $x_1$. Choosing the housing price data (hprice1.gdt) from Wooldridge available for gretl, I obtain the following estimates for the relevant regression coefficients:               Model1 (price ~ sqft)       ------------------------       const       11.2041                sqrft       0.140211                          Model2 (price ~ sqft + bdrms)       ------------------------       const       -19.315              sqrft       0.128436                       bdrms       15.1982                           Model3 (bdrms ~ sqft)       ------------------------       const       2.00808              sqrft       0.000774748                   So $\tilde{\beta_1}=0.140211$, $\tilde{\delta_1} = 0.000774748$, $\hat{\beta_2}=15.1982$ and $\hat{\beta_1}=0.128436$ But these values do not quite match the expression from above: $0.124836+15.1982*0.000774748 = 0.1366108 \neq 0.140211$ Where is my misunderstanding ? Is the formula above valid only in the limit of infinite sample size?


title: Confirming multicollinearity with small number predictors (after initial model construction)

text: I am working on a dataset with with 100 or so X variables and the effects on a single Y. As with previous threads this contains a lot of fairly significant correlations. I considered using PCA but did not really want to loose the original variables by creating linear combinations (factors).   X variables also in this case have some nasty distributions and some are categorical. However, my Y is nicely normal. With this I decided to use bootstrapping with partition analysis. Great, a few predictors were produced that all seemed to add to the explanation of Y, however on further investigation it is looking as if only one or 2 predictors are responsible for the Y variation. Partial correlation appears to be the reason why others are included. Indeed, when run together in a linear model VIFs are below 10 but many are around 5 Now putting these variables into a PCA produces suspicious results: PCA appears to pull 3 factors of worth with loadings of these. Rotated factor loadings               Factor 1        Factor 2        Factor 3     0.5710189449    0.5594104391    0.0092385097     0.0573348956    0.128427451    0.8688140622     -0.029701697    0.8398533899    0.3350420004     -0.024820744    -0.931953121    0.1360273344     0.434265868      0.8133639053   0.022057529     0.8103764005    0.3739830113    -0.225363624     0.4622020154    0.2581755128    -0.208011559     0.8560481143    0.3198727429    -0.191875778     0.7308835871    -0.08016601     0.2661637317     0.8815834293    0.2432798016    -0.120213308     0.8507444962    -0.044814927    0.3006428831     -0.857600499    -0.205246822    -0.036006619     0.6900628882    0.1578752523    0.0377927821     0.7131252931    -0.289263321    0.3166515143      Is this effective saying that if I do not want to use recombined factors I only really have a couple of variables?   I am unsure if I can use PCA this way.


title: Comparison between fat-matrix

text: I have two data sets obtained from two different kind of chemical analysis: GC-MS and IR. Both the chemical analysis were performed using the same samples. My aim is to know if there are some relations between GC-MS data and IR data. The data sets are two 'fat' matrix where: number of columns (objects) > number of rows (different samples). R is the way I use to analyze them. Which is the better approach for this kind of work?


title: impact coding legitimity and rescaling

text: I have implemented a version of "impact coding" : http://www.win- vector.com/blog/2012/07/modeling-trick-impact-coding-of-categorical-variables- with-many-levels/. I have not see that name anywhere else even though I found a paper on the subject : http://dl.acm.org/citation.cfm?id=507538 that can be found with a ps extension on the web. I was wondering two things. Is this preprocessing legit ? I mean does anybody use it and have good results ? By testing it on some toy datasets I found It could improve my performances on some cases. But I think it may have a too strong effect for some variables, especially when using random forest. For example : the Ticket variable in Titanic Dataset because it poorly generalise (as I used the implementation with lambda = n_i / (n_i+m), lambda is a decreasing function of m so random forest will 'ignore' this parameter and it become hard not to overfit the training set). My second question is : does anyone know if it would be better to rescale the newly created variable if I were to use a lasso, svm, or any method that requires comparable inputs (and if I had other continuous features) ?


title: Doing EOF analysis in R

text: I am very new to R and statistics as a whole so this may be a very simple question. I am trying to carry out empirical orthogonal function (EOF) analysis of sea-level pressure (SLP) data to determine the amplitude of the North Atlantic Oscillation over time. I have a matrix (1000,756), with 1000 years of winter mean SLP (i.e. from December to March) at 756 grid points in the North Atlantic. How do I go about doing an EOF analysis on this data? I have spent the last week trying to do a crash course in R but am a little tight for time so I thought I would ask! Thanks in advance for any help.


title: Normalising higher moments of features for machine learning?

text: It's quite common to normalise different feature vectors so that they have the same mean and variance (usually `(X-mean(X))/sd(X)`), so that the changes in the smaller variables aren't swamped. Is the same ever done for higher moments (skewness, kurtosis, ...), for data where those moments are extreme? If so, are there any precautions that need to be taken before doing so? Edit: the question is fairly general, I should qualify that I'm doing clustering analysis, and would love answers relating to the field. However, more general answers are welcome!


title: How to compare models with different distributional assumptions for response variable in GLM?

text: Let's say I have measurements $Y$ which are all positive, and the distribution seems to be somewhat skewed. I'm modelling $Y$ in GLM framework. Now I could set my GLM using different distributional assumptions for $Y$, i.e. I could set family as normal or gamma, or I could log-transform $Y$ and use normal distribution for that. The problem is, how do I compare these models and figure out which is the best? My thinking is that I cannot use AIC or BIC because of different distributions (or in log-transformation case different response variable values). S what wold be the correct way? Only thing I have come up is checking the histograms of residuals and seeing what they look. Edit: Thinking more, if likelihood is the probability of the data given the model, and if all models have same number of parameters, I can directly compare the log-likehoods, which must contain all the constants also. Can somebody confirm this?


title: p-value of hypothesis "a home court advantage exists"

text: **The practical POV** : I sample 35 matches: 15 turn out as home win, 10 as tie and 10 as home lose. I want to conduct a test for the hypothesis "A home court advantage does exist." with a significance level of say 5%. What would be a reasonable test statistic / calculation of the p-value to accept or reject the hypothesis? **Clarification:** The matches I sample are independent of a specific team. I just take a number of matches and count for the home teams - wins, ties and defeats. Due to symmetry the guest teams will have deterministically corresponding results and hence can be ignored. * * * **The generalized POV** : $A_1$ be the random variable for the number of matches the local team won in a total sample of $N$ matches. Naturally this is a bernoulli experiment (given we only assume the possibility of win and lose) and the **p-value** of the hypothesis "a home court advantage exists - i.e. the probability of a home win $p_1$ is larger than 0.5." given $a_1$ observed wins is: $$\Pr[A_1 >= a_1 | p_1 = 0.5]=\sum_{k={a_1}}^N{N\choose{k}}(0.5)^k$$ So far so good (though please correct me if I am wrong or my wording is problematic) - now the aspect that puzzles me is how to factor in ties. My idea would be to create a test statistic $\tilde{A}_1 := A_1 + \lfloor{A_0/2}\rfloor$ with $A_1$ and $A_0$ being the random variables for number of home wins and ties in $N$ matches. $$\Pr[\tilde{A}_1 >= \tilde{a}_1 | \tilde{p}_1 = 0.5]=\sum_{k={a_1}}^N{N\choose{k}}(0.5)^k$$ Now in my mind this makes sense somehow because it reflects the indecisive and symmetric nature of a tie. On the other hand I am not sure how to bend the bernoulli experiment to reflect the notion of $\tilde{p}_1$ being the probability of a home win in the second model (which means home win or "50% of a tie") while keeping up the pretence of formal correctness. Does this make sense?


title: Frontier analysis techniques course

text: I am currently examining efficiency measurement techniques in Kenyan hospitals. We don't have a lot of resources, tutorials or courses here and hence very challenging to undertake the analysis if not an expert. Anyone please advice on any intensive short courses that concentrate on stochastic frontier analysis and data envelopment analysis. The course(s) can be located in any country.


title: Simple linear regression with a random predictor

text: We understand a SLR model as $$y_i = \alpha + \beta x_i + \varepsilon_i$$ with $\varepsilon_i$ i.i.d with equal variance. Suppose I have two instruments measuring a common entity, say, density of samples of different liquids, and assume that we do not know the exact density of each liquid (nor is there any way to accurately measure it). My objective is to construct a model that will summarise the difference between the two instruments. Here is how I go about it: (for simplicity sake, all models are SLR models) Let readings from instrument 1 and instrument 2 when measuring the $i$-th liquid be $y_i$ and $w_i$ respectively, and the actual density of the $i$-th liquid be $x_i$ Then we have: $$y_i=\alpha_1+\beta_1 x_i + \varepsilon_{1i} \\ w_i=\alpha_2+\beta_2 x_i + \varepsilon_{2i}$$ where $\varepsilon_{1i}$ and $\varepsilon_{2i}$ are i.i.d and equal variance w.r.t $i$. Rearranging from the second equation we have: $$x_i= \frac{w_i- \alpha_2- \varepsilon_{2i}}{\beta_2 } $$ Plugging into the first equation, $$ \begin{align*} y_i &=\alpha_1+\beta_1 \frac{w_i- \alpha_2- \varepsilon_{2i}}{\beta_2 } + \varepsilon_{1i} \\ &= \alpha_1 - \frac{\beta_1}{\beta_2} \alpha_2 + \frac{\beta_1}{\beta_2} w_i + \varepsilon_{1i}-\frac{\beta_1}{\beta_2}\varepsilon_{2i}\\ &= \alpha + \beta w_i + \varepsilon_i \end{align*} $$ where $\beta = \frac{\beta_1}{\beta_2}$, $\alpha = \alpha_1 - \beta \alpha_2$, and $ \varepsilon_i = \varepsilon_{1i}-\beta \varepsilon_{2i}$, and $ \varepsilon_i $ are consequently i.i.d. with equal variance w.r.t. $i$. I would like to ask,   1. Is this a valid SLR model for my objective? Do I just treat readings from instrument 2 as an 'independent variable' and proceed like a normal SLR analysis?   2. If so, can I interpret the slope($\beta$) as the ratio of the individual slopes of the two instruments?   3. Can I then test the hypothesis that the two instruments are similar by testing $\alpha = 0$ and $ \beta =1 $?   4. Suppose there are N total liquids we are interested in. All N liquids have been measured by instrument 1. However, we only have time/money to measure n ($<$N) liquids using instrument 2. How should I select the n liquids from N such that I would minimise variance for $\beta$ and $\alpha$?


title: Helmert Transformations

text: Let $Y_f, Y_1,..., Y_n$ denote measurements taken from a father and his n sons. It is assumed that $Y_f, Y_1,...,Y_n$ have equal variance with father- son correlation $ρ_f$ and equal correlation $ρ_S$ between sons. Transform $Y_f, Y_1,..., Y_n$ linearly to $Z_0, Z_1,..., Z_n$ to remove the correlations. The hint is to let $Z_0=Y_f+c(Y_1+...+Y_n)$ where c is a suitably chosen constant. Here is what I have tried. I let $Z_0, Z_1,..., Z_n$ be the Helmert transformations of $Y_f, Y_1,..., Y_n$ where $Z_j=H_jY$ for j=1,2,...,n and thus, ![enter image description here](http://i.stack.imgur.com/pi4FQ.png) where $Y = Y_1,..., Y_n$ I am stuck until this step and i don't have a clue on how to find the constant c to solve the problem. Thanks for the help.


title: Use the improper prior $p(v) \propto 1/v$ into Jags

text: I know that one can approximate this density ($p(v) \propto 1/v$) using its truncated version and implement it this way:                  B~dunif(log(BInf),log(BSup))        v<-exp(B)      but I would like to use the exact form (I checked that under this prior my posterior is proper).


title: Can I simply subtract an index from another?

text: I am currently using a measure of the price level (CPI) from which I would like to remove the effects of food prices (also an index). Notice that food prices are included in the CPI (so that CPI is the total index). First, I made sure that both were based in the same way (that is, the base value of 100 is on the same time point). Now, if I simply subtract the food price index from the CPI, I end up with small values hovering around 2, so I am quite sure it is not the way. I tried the following: CPI+(CPI-food price) and CPI-(CPI-food price) in order to get rid of the food price index. However, I am not sure whether it is the first or the second method since, obviously, both give me a value of 100 on the base date. My hunch is, it is the CPI-(CPI-food price) but I would like your take on this.


title: Anova with repeated measures?

text: I have conducted a study of children with leukemia where I have looked at the change in BMI z score during treatment with glucocorticoids. I have done some different tests and have seen that BMI z score increases significantly during treatment, but I would like to know whether there is a correlation between the BMI z score change and the amount of doses they have recieved. Each day they recieve 1 dose. I don't have BMI z score for each day for each subject. Eg. BMI z scores for subject 1 are available at day 1, 30, 55,60, 100, 105 & 150, while the BMI z scores for subject 2 are available at day 1, 19, 25, 40, 69, 89, 99, 101, 120, 130, 135, 140 & 150. A friend of mine has suggest an anova with repeated measures where patients are a random factor, doses are a fixed variable and BMI z score are the response variable. Do you agree and if you do how do I preform such a test?


title: How to get correlation between two categorical variable and a categorical variable and continuous variable?

text: I am building a regression model and I need to calculate the below to check for correlations   1. Correlation between 2 Multi level categorical variables or 2 binary variables   2. Correlation between a Multi level categorical variable and continuous variable    3. VIF(variance inflation factor) for a Multi level categorical variables with other variables in the model I believe its wrong to use Pearson correlation coefficient for the above scenarios because Pearson only works for 2 continuous variables. Please answer the below questions   1. Which correlation coefficient works best for the above cases ?    2. VIF calculation only works for continuous data so what is the alternative?    3. What are the assumptions I need to check before I use the correlation coefficient you suggest?    4. How to implement them in SAS & R?


title: Information retrieval performance measures for unknown test collection?

text: I am evaluating a web search relevance feedback algorithm. The algorithm uses Bing API as source for its result sets. To evaluate the algorithm I will be conducting a user study. In the end I will have binary judgements, whether the documents are a good answer to the query. Since the web as a test collection is big, I can't calculate recall, because I don't know the amount relavant documents, which are not found. The only measure I can think of is precision, because it uses only the result set and the judgements. Question: **Are there other performance measures (besides precision), which would work on documents returned by a web search?**


title: Model selection and model performance in logistic regression

text: I have a question about model selection and model performance in logistic regression. I have three models that are based on three different hypotheses. The first two models (lets name them z and x) only have one explanatory variable in each model, and the third (lets name it w) is more complicated. I’m using AIC for variable selection for the w model and then AIC for comparing which of the three models that explain the dependent variable best. I’ve found that the w model has the lowest AIC and now want to do some performance statistics on that model to get some idea about the predictive power of the model. Since all I know is that this model is better than the other two but not how good it is. Since I’ve used all data to learn the model (to be able to compare all three models) how do I go about with model performance? From what I’ve gathered I can’t just do a k-fold cross validation on the final model I got from model selection using AIC but need to start from the beginning with all explanatory variables included, is this correct? I’d think that it is the final model I’ve chosen with AIC that I want to know how well it performs, but do realize that I’ve trained on all data so the model might be biased. So if I should start from the beginning with all explanatory variables in all folds I will get different final models for some folds, can I just choose the model from the fold which gave the best predictive power and apply that to the full data set to compare AIC with the two other models (z and x)? Or how does it work? Second part of my question is a basic question about over-parameterization. I have 156 data points, 52 is 1’s the rest are 0’s. I have 14 explanatory variables to choose from for the w model, I realize that I can’t include all due to over-parameterization, I’ve read that you should only use 10% of the group of the dependent variable with fewest observations which only would be 5 for me. I’m trying to answer a question in ecology, is it ok to select the starting variables which I think explains the dependent best simply based on ecology? Or how do I choose the starting explanatory variables? Doesn’t feel right to completely exclude some variables. So I really have three questions:   * Could it be ok to test performance on a model trained on the full data set with cross-validation?   * If not, how do I choose the final model when doing cross-validation?   * How do I choose the starting variables so I want over-parameterize? Sorry for my messy questions and my ignorance. I know that similar questions have been asked but still feel a little confused. Appreciate any thoughts and suggestions.


Let $Y_1, Y_2, Y_3$ be **i.i.d** continuous random variables. For $i = 1, 2$ define $U_i$ as $U_i=1$ if $Y_{i+1} > Y_i,$ $=0$ ,otherwise Find the mean and variance of $U_1 + U_2$ I can't find the distribution of $U_1 + U_2$.


title: Determing noise type and level of noise

text: **Context** I have a set of data that was collected from several inertial measurement units (orientation and acceleration data). I want to determine to what extent an inference method degrades when the data becomes noisy (or, should I say "noisier"). **Questions** How do I determine what type of noise to use? Is there a way in which I can find this out from the data itself? What range of noise levels would be appropriate for testing the inference method's robustness? My guess is that it would depend on the answer to the question: "how noisy would you expect the environment/system to get"? Again, given only the collected data, what is the best approach in determining these levels?


title: What distribution pluffy to buy for an aspiring econometrician?

text: My girlfriend is an Actuarial Analyst at a large insurance company in the Netherlands and because we'll soon have our two year anniversary, I thought of gifts for her. On Proof: Math is beautiful I discovered these Distribution pluffies. So here's my question: What distribution is of the most relevance in the field of an econometrician? The available pluffies are:   * Standard Normal Distribution   * t Distribution   * Chi-Square Distribution   * Log Normal Distribution   * Continuous Uniform Distribution   * Weibull Distribution   * Cauchy Distribution   * Poisson Distribution   * Gumbel Distribution   * Erlang Distribution Any help much appreciated. **EDIT:** Thanks a lot for all the suggestions despite this being just off- topic! I'll get her the t Distribution pluffy.


title: What does inverse-chi-square in Fisher method (classifying) exactly do?

text: In the book “Programming Collective Intelligence” Segaran explains the Fisher method for categorizing text as an alternative to Naive Bayes classifier. The Fisher method uses inverse-chi-square-distribution, which I do not really understand. I watched this video found on stats.stackexchange about chi-square- distribution to understand at least the “forward” function: http://www.youtube.com/watch?v=dXB3cUGnaxQ Segaran explains in his book that they use inverse chi-square to somehow get a probability “that a random set of probabilities would return such a high number”. With high number he means that an item fitting a specific category has many features with high probabilities in that category. Somehow he also seems to take into account that “if the probabilities were independent and random, the result of this calculation would fit a chi-squared distribution”. But as he mentioned before the words are not independant (which is also a false assumption at naive bayes). So how does this then work? And if I understand it right now, the inverse chi-square function somehow checks if many of my words have a high probability of being in the text and only if all words have such a high probability it returns a high over-all probability? I’m sort of confused. PS: The whole paragraph: “Fisher showed that if the probabilities were independent and random, the result of this calculation would fit a chi-squared distribution. You would expect an item that doesn’t belong in a particular category to contain words of varying feature probabilities for that category (which would appear somewhat random), and an item that does belong in that category to have many features with high probabilities. By feeding the result of the Fisher calculation to the inverse chi-square function”, ou get the probability that a random set of probabilities would return such a high number.”


title: What is the way to represent factor variables in scikit-learn while using Random Forests?

text: I am solving a classification problem using Random Forests. For that I have decided to use Python library scikit-learn. But I am new to both Random Forest algorithm and this tool. My data contains many factor variables. I googled for that and found out that it's not right to give numerical values to factor variables like we do in linear regression, as it will treat it as continuous variable and give wrong result. But I could not find anything about how to deal with factor variables in scikit-learn. Please tell me the options to use or point me to some document where I can get it.


title: Symbolic computer algebra for statistics

text: What functionality should exist in a CAS that was specifically geared toward Statistics? Symbolic algebra systems like Mathematica and Maple are often used for calculus, logic, and physics problems but are rarely used for statistics. Why is this? What statistical constructs could be added to a symbolic algebra system to improve its use in this field? What are some specific code samples that many people would like to be able to do. Please think about the following three users: research statistician, non- statistics researcher using statistics in another field (such as biology), statistics student. I'll be working on SymPy's statistics code over the next few months and would like to solicit input for desired functionality. The things I use are not necessarily what the broader community uses.


title: Using R's phyper to get the probability of list overlap

text: I know that I, and others, sometimes get confused by the hypergeometric distribution (HD) as it pertains to overlapping lists. This is because the HD is usually described with the "balls in an urn" metaphor and not using "overlapping lists." What is the proper way to calculate the p-value, according to the hypergeometric distribution, for the overlap of $B$ and $C$ in the lists below, ideally using the `phyper` function in R, where   * $A$ contains all of the genes in the genome   * $B$ is one subset of genes in the genome   * $C$ is another subset of genes in the genome?


title: Standard error on median for exponential distribution

text: I am trying to find the standard error on the median, $\sigma_\tilde{x}$, for a sample, $X_i$, of a population whose pdf could be modelled as $\lambda e^{-\lambda x} $ if normalized. To make sure we are all on the same page. $\mu = 1/ \lambda $, $\sigma = 1 / \lambda$ and $\tilde{x} = ln(2) / \lambda$ for the pdf. Clearly I can calculate $\mu$, $\sigma$, $\sigma_\mu = \sigma / \sqrt{N}$ and $\tilde{x}$ from my sample but the $\sigma_\tilde{x}$ is more complicated. As explained in Central limit theorem for sample medians, $\sigma_\tilde{x} = \frac{1}{2 \sqrt{N} f(\tilde{x})}$, where $ f(\tilde{x})$ is the value of the pdf at the median. Substituting $\tilde{x} = ln(2) / \lambda$ into the pdf $f(\tilde{x}) = \lambda e^{- \lambda \tilde{x} } = \lambda / 2$ so the formula above gives $\sigma_\tilde{x} = \frac{1}{\lambda\sqrt{N}}$ The problem at hand is that I do not know how to calculate $\lambda$. $\lambda=1/\mu$ and $\lambda=1/ \sigma$ so I could just take the calculation of $\mu$ or $\sigma$ from my sample. Would they give the same exact value on the same set? Furthermore, if $\lambda=1/ \sigma$ is taken then $\sigma_\tilde{x} = \sigma_\mu$ which I have a hard time believing is a coincidence. Alternatively I could fit the sample to the assumed form of the pdf to get $\lambda$. I would think that an appropriate choice of fitting algorithm would return the same value. So what is the "best" way to obtain $\lambda$ from $X_i$ and why?


title: Why does SAS Enterprise Miner keep all dummy variables for a coded categorical variable in stepwise logistic regression?

text: SAS Enterprise Miner nicely creates coded dummy variables for any categorical variables used in a logistic regression model. When it performs a variable selection using stepwise sequential selection in the Regression node, however, if one of the dummy variables is included in the regression model, all of the other dummy variables are then also automatically included, even if they are not found to be predictive of the target. Here's a snippet of the node Results output after the stepwise selection showing that the dummy variables for some of the levels of the Industry variable are significant in the model, but others are not.               Parameter      DF    Estimate       Error    Chi-Square    Pr > ChiSq                  Intercept       1     -9.2383      1.9222         23.10        <.0001     IMP_REP_Age     1      0.3594      0.0938         14.69        0.0001     IMP_UnionSubs           No         1      0.5114      0.1472         12.07        0.0005          Industry       Agriculture   1      1.4439      0.1871         59.54        <.0001       Construction  1      1.2982      0.2228         33.97        <.0001       Finance       1     -0.3826      0.2536          2.28        0.1313       IT            1     -0.1355      0.2641          0.26        0.6080       Professional  1      0.3569      0.3469          1.06        0.3037                        Public Sector 1     -2.3698      0.3522         45.28        <.0001                     Retail        1     -1.3766      0.5483          6.30        0.0120                                 Occupation Type        Casual       1      0.0260      0.2499          0.01        0.9171          Employed     1     -0.9068      0.1828         24.61        <.0001        So, for example, the Industry-Agriculture variable seems predictive of the target, but the Industry-IT variables does not. All seven dummy variables for the seven levels of the Industry variable are included in the final model, however. It seems to me that in the stepwise selection the dummy variables should be treated as individual variables rather than as a group. Does anyone know why SAS Enterprise Miner does it differently?


title: Acceptance rejection method for a binomial distribution

text: I have a R-code for acceptance rejection method for a beta-distribtuion but would like to use it for a binomial(30, 0.4) distribution instead. Here is the code:               sample.x = runif(100000,0,1)     accept = c()          for(i in 1:length(sample.x)){      U = runif(1, 0, 1)      if(dunif(sample.x[i], 0, 1)*3*U <= dbeta(sample.x[i], 6, 3)) {        accept[i] = 'Yes'      }      else if(dunif(sample.x[i],0,1)*3*U > dbeta(sample.x[i], 6, 3)) {       accept[i] = 'No'      }     }      Is it possible to modify it so that it works for my binomial distribution instead? Thanks in advance! Edit: So I'm trying to implement Glen's answer to my question, but I can't get it to work as I want. Does anyone have a solution? Here's my new r-code:               sample.x = runif(10000,0,1)     accept = c()          for(i in 1:length(sample.x)){      U = runif(1, 0, 1600)      if(dunif(1, 0, 1)*1600*U <= dbinom(1, 30, 0.4)) {        accept[i] = 'Yes'      }      else if(dunif(1,0,1)*1600*U > dbinom(1, 30, 0.4)) {       accept[i] = 'No'      }     }


title: Consequence of violating independence assumption of ANOVA

text: I would like some help with the following problem. I have 40 subjects. On each subject I take a measurement at 25 body sites. The measurement is a continuous variable that varies between 0 and 1 and appears to be normally distributed. I want do a test to see if body site statistically significantly affects the measurement. I also have demographic data on the subjects but haven't tried to analyze that yet. I think the correct thing to do would be to use repeated measures ANOVA or MANOVA to do the analysis because the data is paired. The problem is that I have missing data and my basic understanding is that I would have to exclude any subject that has data missing for even one body site. If I did this, I would only have 6 subjects left. If I excluded any body site with data missing then I would only have 8 body sites left. So my thought was that I could use regular ANOVA but just lose some power to see a difference. However, this seems like it violates the assumption of independence for ANOVA. So my main questions are as follows:   1. Does using Regular ANOVA on matched data violate the independence assumption, and if so is the only consequence a decreased likelihood of rejecting the null hypothesis or is there more to it than that? Also does anyone know of a reference I could site regarding this?   2. Does anyone have any suggestions for a better way to analyze my data? By the way I usually use JMP or SPSS for analysis. Thanks! Billy


title: How to present the gain in explained variance thanks to the correlation of Y and X?

text: I'm searching how to (visually) explain simple linear correlation to first year students. The classical way to visualize would be to give an Y~X scatter plot with a straight regression line. Recently, I came by the idea of extending this type of graphics by adding to the plot 3 more images, leaving me with: the scatter plot of y~1, then of y~x, resid(y~x)~x and lastly of residuals(y~x)~1 (centered to the mean) Here is an example of such a visualization: ![alt text](http://i.stack.imgur.com/Pe2ul.png) And the R code to produce it:               set.seed(345)     x <- runif(50) * 10     y <- x +rnorm(50)               layout(matrix(c(1,2,2,2,2,3 ,3,3,3,4), 1,10))     plot(y~rep(1, length(y)), axes = F, xlab = "", ylim = range(y))     points(1,mean(y), col = 2, pch = 19, cex = 2)     plot(y~x, ylab = "", )     abline(lm(y~x), col = 2, lwd = 2)          plot(c(residuals(lm(y~x)) + mean(y))~x, ylab = "", ylim = range(y))     abline(h =mean(y), col = 2, lwd = 2)          plot(c(residuals(lm(y~x)) + mean(y))~rep(1, length(y)), axes = F, xlab = "", ylab = "", ylim = range(y))     points(1,mean(y), col = 2, pch = 19, cex = 2)      Which leads me to my question: I would appreciate any **suggestions on how this graph can be enhanced** (either with text, marks, or any other type of relevant visualizations). Adding relevant R code will also be nice. One direction is to add some information of the R^2 (either by text, or by somehow adding lines presenting the magnitude of the variance before and after the introduction of x) Another option is to highlight one point and showing how it is "better explained" thanks to the regression line. Any input will be appreciated. Cheers, Tal


title: Fisher's exact test or chi-square test

text: I have a 2x4 table with nominal data (the columns are simply YES/NO, the rows are four categories)               Category A: 7, 13     Category B: 15, 5     Category C: 15, 5     Category D: 19, 1       I am hoping to test the significance of a couple of the categories to each other (2x2), but also to assess the significance of the whole table (2x4, although I am not really sure how to interpret the meaning of the significance that may result here). I understand that as my sample sizes are small (as you can see, Category D features "1" under NO, and each category has only 20 people) that I should be using Fisher's Exact Test. Is this correct? Can Fisher's also be used in a 2x4? And what does it mean if this result is significant?


title: determining drivers of loss rate change between two time periods (binomial regression)

text: I've read through the R Book and I think is the closest I've come to the problem. I have two time periods - say TEST1 and TEST2 where I have aggregated samples of data of how many 'animals' were lost in that time period compared to the starting number of animals. I would like to identify what's driving the difference in the outcomes in TEST1 and TEST2. Is it reasonable to compare the interaction coefficients on the covariates with the test? For example test:animal+test:state This is the mock data I'm using to get comfortable with the model.                   require(data.table)         df <- data.table(expand.grid(test=c("TEST1", "TEST2"), animal=c("cat", "dog", "mouse", "llama"), state = state.abb))         set.seed(5)              df[,Ct:= as.integer(runif(nrow(df), 1, 1000))]              df[,Lost:=as.integer(Ct*.1)]         df[animal=='cat', Lost:=as.integer(Ct*.5)]         df[test=='TEST2', Lost:=as.integer(Lost * 1.2)]         df[test == 'TEST1' &( substring(state,1,1)=='M' | substring(state,1,1) == 'N'), Lost:=as.integer(Lost * .8)]                   model <- glm(cbind(df[,Lost], df[,Ct])~df[,animal]+df[,test]+df[,test]:df[,animal]+df[,test]:df[,state], data=df, family=binomial)         summary(model)      The model is correctly identifying some of the differences in the tests:                                                       Estimate Std. Error z value Pr(>|z|)             (Intercept)                       -6.958e-01  7.106e-02  -9.792  < 2e-16 ***         df[, animal]dog                   -1.621e+00  2.691e-02 -60.248  < 2e-16 ***         df[, animal]mouse                 -1.620e+00  2.528e-02 -64.062  < 2e-16 ***         df[, animal]llama                 -1.619e+00  2.487e-02 -65.092  < 2e-16 ***         df[, test]TEST2                    1.851e-01  8.568e-02   2.160  0.03073 *           df[, animal]mouse:df[, test]TEST2 -4.016e-03  3.434e-02  -0.117  0.90691             df[, animal]llama:df[, test]TEST2 -4.552e-03  3.453e-02  -0.132  0.89513             ....         df[, test]TEST1:df[, state]ME     -2.262e-01  9.096e-02  -2.487  0.01287 *            ....                   Null deviance: 22533.373  on 399  degrees of                  Residual deviance:    14.249  on 294  degrees of freedom         AIC: 2483.2      I'm interpreting that as no incremental effect of animal on TEST1/TEST2 (since the interaction was not significant). I see an significant interaction on state for test1 for any state with M/N as the first letter (as expected). I'm also seeing TEST2 have a general higher level. Is this the right way to approach the problem? Ideally I would like to report back which variables contributed most to the difference in test1/test2 (or time period 1 and time period 2).


title: Combining critic scores into single rating

text: I have a collection of scores from critics for various restaurants. Some restaurants are rated by multiple critics and some by only a single critic. Each critic rates restaurants using different scales such as 1 - 5 (stars), 1 - 10, 1 - 100 etc. Most critics rate a restaurant only once although they may revisit it at a certain point (probably this complication can be ignored and I will use just the most recent rating). Also some critics only use a subset of the range (ie 60 - 100). This may be due to the fact that they only review reasonably good restaurants and hence don't use the bottom end of the range, or possibly other reasons. In addition the distribution of ratings given by an individual critic are often quite skewed (left or right). Can anyone suggest some good ways of combining individual critics scores into a single score ? I thought about standardizing each critics score and then taking the average. However this creates a couple of problems:   * There is the case where multiple critics give a high rating and one might give a very low rating. If so should I ignore the outlier or give them equal weighting? Or should I assign a given critic some sort of reliability rating and use that to create a weighted average. How might one calculate this reliability ?   * I'm not sure I should drop the information that a critic only uses a subrange of the scale such as 60 - 100 instead of 0 - 100 as this may reflect the fact that they are being selective about which restaurants they review. In other words an above average z-score for critic A (who picks only good restaurants), may mean something quite different than an above average z-score for critic B (who dines at a wider range of restaurants).  Sorry for the vague question, but does anyone have any suggestions about how I approach this ? Is there anyway to use the fact that I have multiple ratings for a large number of restaurants to somehow calibrate the critics' ratings to a common scale ? Or is this nonsensical. Should I be looking at some sort of random effects model ? Thanks for any pointers or suggestions.


title: Minima between modes at different time points

text: I have measurements for a number of different discrete groups (e.g., Group A, Group B, . . . Group Z). I believe that the measurements for Groups A and B are smaller than those for groups C through Z. I would like to know the best way to demonstrate this, taking sample size into account. So far I've been using ANOVA. The main reason I'm writing is that I would like to establish a legitimate threshold for "small." I'd like to find a point (or the points) that best represents minima between modes. Also, though my measurements are continuous, the data are all from discrete points. The minimum/cutoff that I'm looking for would be the one that best represents a minimum across all time points. I'm looking for something along the lines of this but that would optimize the minimum across timepoints. I'd vastly prefer something that works in R.


title: Parameter Estimates

text: I have copied my Parameter estimates table below and am not sure how to interpret it or what it means. Any help would be greatly appreciated! Warmly, Carolynn Dependent Variable Parameter B Std. Error t Sig. 95% Confidence Interval Lower Bound Upper Bound neutral pre stai Intercept MASCscore [Gender0F1M=0] [Gender0F1M=1] STAI anxiety pre Intercept MASCscore [Gender0F1M=0] [Gender0F1M=1] neutral POST SCORES Intercept MASCscore [Gender0F1M=0] [Gender0F1M=1] STAI anxiety post Intercept MASCscore [Gender0F1M=0] [Gender0F1M=1] 29.363 1.902 15.439 .000 25.575 33.151 .070 .045 1.566 .122 -.019 .160 2.086 1.381 1.511 .135 -.664 4.836 0a . . . . . 35.076 1.875 18.702 .000 31.341 38.812 -.039 .044 -.888 .377 -.128 .049 1.076 1.362 .790 .432 -1.636 3.788 0a . . . . . 24.827 1.653 15.020 .000 21.534 28.119 .174 .039 4.455 .000 .096 .252 -.934 1.200 -.778 .439 -3.324 1.457 0a . . . . . 27.052 2.128 12.715 .000 22.815 31.289 .122 .050 2.433 .017 .022 .222 .639 1.545 .413 .680 -2.438 3.715 0a . . . . . a. Page


title: How to use error term in AR (2) model for predicting future values?

text: We use turbidity to estimate suspended-sediment concentration (SSC)- our data was serially correlated. We ran an ARMA process and ended up with a AR (2) model. Our equation in log form is: estlogi(SSC) = 1.35149 + 0.785099logi(Turbidity) - 0.00095 + 0.7492ei-1+ 0.2052ei-2 How can we input the ei-1 & ei-2 values in real time, that is, how do we enter residuals for current or future values when our turbidity/SSC data pairs are in the past (2010 to 2011). Someone did suggest an R script that will estimate K steps ahead for a specified ARMA model and past data- does anyone know this script?


title: Time series dynamic poisson regression

text: I have a time series count data by customers that I would like to regress on past months items (count) sold and promotional effects (current and past). Below is an example, and the dataset has one record per customer. Y= number of items (count) sold **Customer 1:** Y/Count data (current month) = Y/count data(past month) + promo data (current month) + promo data (past month) **Customer 2:** Y/Count data (current month) = Y/count data(past month) + promo data (current month) + promo data (past month) I'm not sure what is the right model to use. I could use Dynamic regression/ARIMAX if my response was a continous variable. But the response variable in my hand is a count data, so I'm not sure if ARIMAX is the right model for this problem. In addition, I have lot of zeros in my response variable (similar to zero inflated poisson regression). I did search in standard textbooks and the internet, I did not come across any regrssion model with time series count data. I would really appreciate if anyone could recommend an appropriate model for this problem. Thanks


title: Data Sets suitable for k-means

text: My task is to determine which ones of these datasets given in the picture are suitable for k-means. My script says that k-means usually performs well on concave structures, however the sizes of the obvious clusters of case 1 differ very strongly, which would indicate a possible bad performance of k-means, right? So my pick here would be case 2, since the sizes of the 3 clusters are almost identical and the shapes are - however not concave - very dense and case 3 has a lot of noise, which k-means can't handle. Are my thoughts correct? ![enter image description here](http://i.stack.imgur.com/Uk23U.png) Note, that i don't actually know any values for the dataset, so it is supposed to be clear by looking at the picture if k-means is suitable


title: PCA matrix-decomposition in image recognition

text: I am implementing an image recognition program as an exercise but I am comfused in the following point.I have a data set M with dimensions $NxD$ where $N<<D$ so I have used dimension reduction to compute the eigen- decomposition of the covariance matrix $M'*M$ From this matrix I can't compute the eigen decomposition because it is very large but I can find the first N components from the matrix $M*M'$ and $Vi=Ei*M$ (Ei are the eigenvectors found from the matrix M*M') is the eigenvector for the original covariance matrix $M'M$.Is this procedure correct?I have found and other sources in the web using the similar procedure (Bishop:Machine Learning and Pattern Recognition,Chapter 12.1.4) and from this webpage:http://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=8&ved=0CEkQFjAH&url=http%3A%2F%2Fwww.face- rec.org%2Falgorithms%2FPCA%2Fjcn.pdf&ei=VjmZUKSwGurD0QW0woCoBA&usg=AFQjCNEr3xWDIJcmhi_ROgXQY61lLjnQvA&sig2=BfRoGZxt2n_MPWkQxeX9lg


title: How do I estimate sample size for a single-sample equivalence test?

text: * From literature I have a theoretically proposed value for human performance in a particular task.    * I want to experimentally measure performance in this task across a sample of humans and test whether the result I get is equivalent to the theoretical value.    * It seems that a Single Sample Equivalence Test can do this. Given that I do not know population variance, or mean, only sample variance and mean, how can I go about a power analysis to calculate the sample size I need to use? I have so far looked into the R packages PowerTOST and Equivalence, but have been confused somewhat by the fact that "single-sample" TOST is quite rare, I don't understand what the "theta0" value is and the fact co-efficient of variance seems to be calculated from population rather than the sample. If possible I would like to use R to do this, though I am open to other free solutions. Many thanks.


title: Alternative to 2-way ANOVA when variances are unequal

text: I'm working on a question on 2-way ANOVA where there are 2 treatments and 1 blocking factor. The residual plots suggest unequal variance. What is another parametric or non-parametric alternative method to 2-way ANOVA?


title: non constant variable suggested power transformation question

text: > ncvTest(alm)     Non-constant Variance Score Test      Variance formula: ~ fitted.values      Chisquare = 121.2316    Df = 1     p = 3.400245e-28           > spreadLevelPlot(alm)     Suggested power transformation:  4.428269      I am having issues adjusting my regression formula based on what the results of the non-constant variable test shows. How do I implement the suggested power transformation here? Obviously with such a low p value this is heteroscedastic. I have run a robust standard error linear regression already and it does not change the BP test p value. I have also performed a variance inflation factor test to see if multicolinearity is an issue here, which it does not seem to be. Any and all help is appreciated!.


title: Lizards probability

text: I am trying to learn basic probability and I would like to get the answer of the questions below. I honestly have no idea on how to even start solving it. Could someone show me or give me some advice on how to tackle it? Wild lizards are found in a particular mountain range in Australia at an abundance of 1.4 per 6 square miles. The standard deviation of lizard abundance is 0.4 per 5 squares miles. Assume that lizard abundance is a random variable that has a bell-shaped distribution. a. If you pick a random 12 square mile area in this mountain range, what is the likelihood that there are more than 7 lizards in this area? b. If you pick a random 5 square mile area in this mountain range, what is the likelihood that there are between 1 and 3 lizards in this particular 5 square mile area?


title: What statistical model is used to calculate the test results for GWO?

text: Google Website Optimizer (GWO) is a tool provided by Google to do A/B and MVT experiments on websites. This has been an unanswered question for a long time so I thought I'd ask it here and see if I can get any help. Here is some documentation (clues) that Google has published about the statistics used:   * All about statistics   * Fractional versus Full Factorial analysis Example Data (summary pdf):   * XML   * CSV I (as well as many others) are interested in discovering the math used by Google, because no one has been able to reproduce their calculations yet. Thank you in advance for any help that you can provide.


title: Regression output and Fisher-Snedecor distribution

text: I'm working on regression models in STATISTICA application and I need to know what is Fisher-Snedecor distribution for and how to analyze my regression model in this distribution. What the significance level means? What is v1 and v2? I need an explanation and little tutorial on real data. Thanks in advance!


title: How do I use Q-Learning to update values?

text: I am having trouble using the Q-learning algorithm to solve the problem in the link below. http://math.stackexchange.com/questions/357397/how-do-i-use-q-learning-to- update-values Any help would be appreciated (couldn't post image as I don't have enough badge points)


title: Comparing the variance of paired observations

text: I have $N$ paired observations ($X_i$, $Y_i$) drawn from a common unknown distribution, which has finite first and second moments, and is symmetric around the mean. Let $\sigma_X$ the standard deviation of $X$ (unconditional on $Y$), and $\sigma_Y$ the same for Y. I would like to test the hypothesis $H_0$: $\sigma_X = \sigma_Y$ $H_1$: $\sigma_X \neq \sigma_Y$ Does anyone know of such a test? I can assume in first analysis that the distribution is normal, although the general case is more interesting. I am looking for a closed-form solution. Bootstrap is always a last resort.


title: How to compare performance of multiple regression models?

text: I am taking a little data analysis course using SAS software and I need help with pretty much the basics. There is this full description in the American Statistical Association page. Basically I have the data obtained from caught fish (species, weight, length etc)                   Obs    Spec  Wt        Lt1     Lt2     Lt3     HtP    WhP    Sex     ----+----1----+----2----+----3----+----4----+----5----+----6----+----7----+-         1      1     242.0     23.2    25.4    30.0    38.4   13.4   NA         2      1     290.0     24.0    26.3    31.2    40.0   13.8   NA         3      1     340.0     23.9    26.5    31.1    39.8   15.1   NA         4      1     363.0     26.3    29.0    33.5    38.0   13.3   NA         5      1     430.0     26.5    29.0    34.0    36.6   15.1   NA         ...      I want to predict weight using other known characteristics. How can I prove of disprove that the same regression function would work for all fish species? Would checking hypothesis that the weight distributions (for different species) are statistically different be enough? Then I have a bunch of models apparently suggested by experts               Weight=a+b*(Length3*Height*Width)+epsilon     Log(Weight)=a+b1*Length3+epsilon     Weight^(1/3)=a+b1*Length3+epsilon     Log(Weight)=a+b1*Length3+b2*Height+b3*Width+epsilon     Weight^(1/3)=a+b1*Length3+b2*Height+b3*Width+epsilon     Weight=a*Length3^b1*Height^b2*Width^b3+epsilon      How should I compare them? Simply look for the best root mean square error?


title: web analytic: average session length based on visits numbers per minute

text: At hand, I only have the number of visits per minute for one web page (one day period), is it possible to generate the average session length of such page based on these numbers? Or and references I can check for related information? Thank you in advance.


title: Interpretation of mean in this example

text: I recently presented a national test and the company in charge of preparing the test then does a standardization to provide the final scores for each person. These are the values they gave at the end of the page where the grades were posted:               N = 259     Mean = 72     StdDev = 7      I don't remember anything about this, and my career doesn't deal with any of this. I want to know how can I interpret these results with those values.   * What does it mean that the mean was 72?    * What does it mean that the standard deviation value is 7?    * What does it mean if, for example, I obtained an 88 or a 78 on the test?


title: Mixed thoughts about the Durbin-Wu-Hausman Test. Is it really feasible?

text: I just studied (again) some statements of the Durbin-Wu-Hausman (DWH) test. **After I thought about the theory behind the DWH test, is it not basically kind of useless since I cannot tell if there's actually a problem regarding endogeneity or not?** For completeness, the hypothesis of the test: $H_0: \ \ y = X\boldsymbol{\beta} + u, \ \ \ \ u \sim IID\left(0,\sigma^2I_n\right), \ \ E(X^Tu)=0$ $H_1: \ \ y = X\boldsymbol{\beta} + u, \ \ \ \ u \sim IID\left(0,\sigma^2I_n\right), \ \ E(W^Tu)=0$ Where $W$ is the matrix of instruments. If I'm not mistaken, we need to define some instruments in order to perform the DWH test even if we are not sure about the potential endogeneity of some of our regressors $X$. Now, in case I observe a significant outcome of the DWH test, that is, I reject $H_0$, two interpretations are possible:   1. The first is that at least some columns of $X$ in $H_0$ are endogenous such the $H_1$ is correctly specified. This interpretation would the desirable one.   2. The second involves the regressions $y = X\beta + P_W Y\delta + u$ where $Y$ are the columns of $X$ which are assumed to be endogenous (i.e. a submatrix of $X$) and $P_W$ is the orthogonal projection matrix onto the column space of $W$. If now $\hat{\delta}$ (estimated via OLS) differs significantly from $0$ (tested via F-test) then we would know that the nature of the $H_0$ model is not that columns of $X$ are actually endogenous rather then the linear combinations of the instruments, i.e., $P_WY$ have explanatory power for the dependent variable $y$. This would actually tell us that the instruments which we used are not actually good instruments because there is not supposed to be any explanatory power within the instruments if the regressors $X$ are involved.


title: Latent Dirichlet Allocation (LDA): What exactly is inferred?

text: I am working my way through LDA and I think I got they main idea of it. Please correct me if I am wrong. Given the Plate notation: ![LDA model](http://i.stack.imgur.com/8nzTh.png) The variables $\alpha$ and $\beta$ are Dirichlet distribution parameters. The variable $Z_{d,n}$ assigns observed word $W_{d,n}$ to topic $\phi_k$, which is a distribution over words. Variable $\theta_d$ is the document-specific topic distribution. Both distributions $\theta_d$ and $\phi_k$ are drawn from Dirichlet distributions. Now, only $W_{d,n}$ is observed and can be "directly" calculated. My question: What exactly is inferred/calculated with e.g. Gibbs sampling, variational Inference and so on? * * * For instance: For a Gaussian Naive Bayes classifier one assumes that the likelihood of each feature is Gaussian. In other words each feature has a Gaussian distribution: $$ P(x) = \frac{1}{{\sigma \sqrt {2\pi } }} e^{ \frac{ - ( {x - \mu } )^2} {2\sigma ^2 } } $$ To find this distribution $\sigma$ and $\mu$ have to be determined which is pretty straight forward. However, plainly said: What Numbers do I determine for LDA?


title: How to calculate turbulence and similarity index with TraMineR?

text: I am new to sequence analysis and wondering whether TraMineR can calculate Elzinga's turbulence and similarity index or I should learn to use his CHESA software to do this? I appreciate any thoughts as I learn to do sequence analysis.


title: How do I get the p-value when using Java's OLSMultipleLinearRegression?

text: I'm doing my regression using Apache's OLSMultipleLinearRegression class. How do I get the p-value? Is there another name for p-value that I should be looking for in the docs? EDIT: I have the regression params and the standard error for each param. Is it reasonable to calculate my own p-values? Any reason this wouldn't be in the Apache library?


title: Asymptotic assumptions in OLS

text: Let us assume the following data generating process: (1) $ \ \ y = X\beta + u$ where $y$ and $u$ are $n\times 1$ vectors, $X$ is a $n\times k$-matrix with $rk(X)=k$ and $\beta$ is a $k\times 1$-vector of coefficients. For the OLS- estimator to be unbiased we need to assume that the conditional expectation of $u$ given $X$ i.e. $E(u|X)=0$ holds for (1). If the data are actually generated with (1) and $\beta_0$ is the true value for $\beta$ then we can write $E(\hat{\beta}|X) = \beta_0 + (X^TX)^{-1}X^TE(u|X) = \beta_0$. **My question is: does the assumption $E(u|X)$ already imply by the law of large numbers (LLN) that $\text{plim}_{n\to\infty} \ \frac{1}{n}X^Tu =0$?** My guess is yes because by LLN $\text{plim}_{n\to\infty} \ \frac{1}{n}X^Tu = \text{lim}_{n\to\infty}\frac{1}{n}E(X^Tu) = 0$ ? The reason for this question is that if the assumption $E(u|X)=0$ does imply that $\text{plim}_{n\to\infty} \ \frac{1}{n}X^Tu =0$ than the asymptotic consistency of $\hat{\beta}$ follows directly from this particular assumption.