title: Changing reference level for contrasts changes results in R 3.0.2/lme4 1.1-2 vs. R 2.15.0/lme4 .999999-0

text: I hope this is an appropriate forum to post this question. I recently upgraded my R software from 2.15.0 to 3.0.2. I also upgraded the lme4 package from .999999-0 to 1.1-2. After doing so, the results from one of my linear mixed models analyses have changed a bit unexpectedly. In some respects, I was expecting some change, as the lme4 developers very clearly stated that they had made some significant changes to some fundamental components in the package. However, the changes that I am seeing (described below) make me think that something else is awry. I will start by explaining the experimental design, which is quite simple and then the issue at hand. My experiment is a basic repeated measures design. I used 24 "Items" that each appeared in three different "Conditions" (SmallClause_Som, NoSmallClause, SmallClause_NoSom). Levels of Condition were rotated across three presentation lists such that each Subject (45 total, each assigned to a particular list) only saw one level of each item. I used lmer() for the analysis. Condition was entered in as a Fixed effect and "Subject" and "Item" were entered as Random effects. The problem: Using the current version of R 3.0.2 and lme4 1.1-2 with NoSmallClause as the reference level (and no weighting on any of the contrasts), the ConditionSmallClause_Som/NoSmallClause contrast produces a t value of 1.680. But, when I change reference level to SmallClause_Som (to observe the one remaining contrast) I get not only a change in the polarity of the effect (plus to minus, as expected), but the values change as well. When I use R 2.15.0 and lme4 .999999-0 (on another computer), I do not experience this issue. I get slightly different values, but they do not change (apart from the polarity) when I change reference level. My colleague also tried my analysis for me using R 3.0.2 and a version of lme4 (pre version 1.0) (I don't know exactly which version, but it was before the major changes) and he also does not experience the issue. R 2.15.0 lme4 1.1-2 (older) output:               > #set ref level     > test$Condition <- relevel(test$Condition, ref="NoSmallClause")          > #Model 4: Random slopes by Subject and Item     > test.lmer4=lmer(Rating~Condition+(1+Condition|Subject)+(1+Condition|Item),test, REML=FALSE)     > summary(test.lmer4)          Linear mixed model fit by maximum likelihood      Formula: Rating ~ Condition + (1 + Condition | Subject) + (1 + Condition |      Item)         Data: test       AIC  BIC logLik deviance REMLdev      3930 4010  -1949     3898    3902     Random effects:      Groups   Name                       Variance   Std.Dev. Corr                Subject  (Intercept)                0.98998765 0.994981                              ConditionSmallClause_Som   0.00203374 0.045097 -1.000                       ConditionSmallClause_NoSom 0.00019873 0.014097  1.000 -1.000       Item     (Intercept)                0.96231875 0.980978                              ConditionSmallClause_Som   0.89924400 0.948285 -0.020                       ConditionSmallClause_NoSom 0.62128577 0.788217 -0.256  0.361       Residual                            1.68810777 1.299272                    Number of obs: 1080, groups: Subject, 45; Item, 24          Fixed effects:                            Estimate Std. Error t value     (Intercept)                  2.9583     0.2584  11.447     ConditionSmallClause_Som     0.3639     0.2165   1.680     ConditionSmallClause_NoSom   0.1472     0.1878   0.784          Correlation of Fixed Effects:                 (Intr) CnSC_S     CndtnSmlC_S -0.116            CndtnSmC_NS -0.260  0.392          > #set ref level     > test$Condition <- relevel(test$Condition, ref="SmallClause_Som")     >      > #Model 4: Random slopes by Subject and Item     > test.lmer4=lmer(Rating~Condition+(1+Condition|Subject)+(1+Condition|Item),test, REML=FALSE)     > summary(test.lmer4)          Linear mixed model fit by maximum likelihood      Formula: Rating ~ Condition + (1 + Condition | Subject) + (1 + Condition |      Item)         Data: test        AIC  BIC logLik deviance REMLdev      3930 4010  -1949     3898    3902     Random effects:      Groups   Name                       Variance  Std.Dev. Corr                Subject  (Intercept)                0.9023239 0.949907                              ConditionNoSmallClause     0.0020340 0.045099 1.000                        ConditionSmallClause_NoSom 0.0035039 0.059194 1.000  1.000        Item     (Intercept)                1.8238288 1.350492                              ConditionNoSmallClause     0.8992237 0.948274 -0.687                       ConditionSmallClause_NoSom 0.9804329 0.990168 -0.604  0.670       Residual                            1.6881050 1.299271                    Number of obs: 1080, groups: Subject, 45; Item, 24          Fixed effects:                                Estimate Std. Error t value     (Intercept)                  3.3222     0.3174  10.468     ConditionNoSmallClause      -0.3639     0.2165  -1.680     ConditionSmallClause_NoSom  -0.2167     0.2243  -0.966          Correlation of Fixed Effects:                 (Intr) CndNSC     CndtnNSmllC -0.588            CndtnSmC_NS -0.521  0.638      R 3.0.2 and lme4 1.1-2 (newer) output:               > #set ref level     > test$Condition <- relevel(test$Condition, ref="NoSmallClause")          > #Model 4: Random slopes by Subject and Item     > summary(test.lmer4)          Linear mixed model fit by maximum likelihood ['lmerMod']     Formula: Rating ~ Condition + (1 + Condition | Subject) + (1 + Condition |      Item)         Data: test                 AIC       BIC    logLik  deviance       3942.557  4022.312 -1955.278  3910.557           Random effects:      Groups   Name                       Variance Std.Dev. Corr             Subject  (Intercept)                0.9522   0.9758                             ConditionSmallClause_NoSom 0.1767   0.4204    0.03                     ConditionSmallClause_Som   0.1760   0.4196   -0.15  0.92      Item     (Intercept)                1.2830   1.1327                             ConditionSmallClause_NoSom 0.7782   0.8822   -0.41                     ConditionSmallClause_Som   1.4901   1.2207    0.09  0.41      Residual                            1.6466   1.2832                   Number of obs: 1080, groups: Subject, 45; Item, 24          Fixed effects:                                Estimate Std. Error t value     (Intercept)                  2.9583     0.2814  10.512     ConditionSmallClause_NoSom   0.1472     0.2133   0.690     ConditionSmallClause_Som     0.3639     0.2741   1.327          Correlation of Fixed Effects:                 (Intr) CSC_NS     CndtnSmC_NS -0.357            CndtnSmlC_S -0.007  0.451     > #anova (test.lmer3, test.lmer4)     >      > #set ref level     > test$Condition <- relevel(test$Condition, ref="SmallClause_Som")     >      > #Model 4: Random slopes by Subject and Item     > test.lmer4=lmer(Rating~Condition+(1+Condition|Subject)+(1+Condition|Item),test, REML=FALSE)     > summary(test.lmer4)     Linear mixed model fit by maximum likelihood ['lmerMod']     Formula: Rating ~ Condition + (1 + Condition | Subject) + (1 + Condition |      Item)         Data: test                 AIC       BIC    logLik  deviance       3951.357  4031.113 -1959.679  3919.357           Random effects:      Groups   Name                       Variance Std.Dev. Corr             Subject  (Intercept)                0.88980  0.9433                             ConditionNoSmallClause     0.04299  0.2073   0.83                      ConditionSmallClause_NoSom 0.01562  0.1250   0.90  0.67       Item     (Intercept)                2.39736  1.5483                             ConditionNoSmallClause     0.72053  0.8488   -0.04                     ConditionSmallClause_NoSom 1.87804  1.3704   -0.16  0.53      Residual                            1.65166  1.2852                   Number of obs: 1080, groups: Subject, 45; Item, 24          Fixed effects:                                Estimate Std. Error t value     (Intercept)                  3.3222     0.3525   9.425     ConditionNoSmallClause      -0.3639     0.2004  -1.816     ConditionSmallClause_NoSom  -0.2167     0.2963  -0.731          Correlation of Fixed Effects:                 (Intr) CndNSC     CndtnNSmllC -0.045            CndtnSmC_NS -0.160  0.514      My question: What is going on here? Why is changing the reference level producing a shift from 1.327 to -1.816 in the t scores for the new version of lme4 whereas it produces the same (disregarding sign) value of 1.680/-1.680 in the old version's t scores? Only the older version seems to make sense to me. 1) Am I specifying my model incorrectly for the new version of lme4? 2) Am I missing some basic fundamental fact about how contrasts work? That is, is it possible to get different values just from changing the reference level? (the correlation values look a bit odd in the newer output). 3) Is this a bug in lme4? 4) Some other explanation...? I have had some other odd issues as well with this same analysis using lme4 1.1-2. For example, if I don't clear the workspace and re-run an analysis, the values also will change between analyses (and also within the analysis as I change the reference level). This never happened to me on the earlier version (and it still does not happen when I run it on the earlier version now). I hope someone can help with this. I found two other similar questions online (after much searching) but neither had any informative responses. Thanks DT


title: How to show that the variance of Sequential Importance Sampling estimates increase with the dimension?

text: I am trying to understand the Particle Filter and the motivation to use it over the regular Sequential Importance Sampling. As far as I understand until now: 1- We try to estimate the expectation of a function $f_n:\mathbb{R}^{n} \mapsto \mathbb{R}$ under the multivariate density $\pi_{n}(x_{1:n})$, namely $E_{\pi_{n}(x_{1:n})}[f_n(x_{1:n})]$ by using $N$ samples. Typically this is a dynamic system and $n$ in $\pi_{n}(x_{1:n})$ increases in each step, for example as in the filtering density of a Hidden Markov Model. 2- In each step, we sample from an incremental proposal density like $q_{n}(x_{n}|x_{1:n-1})$. In total, we have $N$ samples from an $n$ dimensional density $q_{n}(x_{1:n}) = q_{1}(x_{1})\prod_{k=2}^{n}q_{k}(x_{k}|x_{1:k-1})$ at time $n$. 3-We calculate the estimate $\frac{1}{N}\sum_{i=1}^{N}f_n(x_{1:n}^{i})w_n(x_{1:n}^{i})$ where $w_n(x_{1:n}^{i})$ is the weight function with $w_n(x_{1:n}^{i}) = \frac{\pi_{n}(x^i_{1:n})}{q_{n}(x^i_{1:n})}$. Now, I understand that we cannot calculate the exact optimal proposal density at each step which minimizes the estimate variance, since we have only the freedom of selecting $q_{n}(x_{n}|x_{1:n-1})$, the rest of the proposal is already fixed. What I don't understand is, why the estimate variance must increase as the dimension $n$ increases; this is explained as the primary motivation for applying Particle Filter solutions. In order to see that, I calculated the variance of the Monte Carlo estimate $Var[\frac{1}{N}\sum_{i=1}^{N}f_n(x_{1:n}^{i})w_n(x_{1:n}^{i})]$ which turns out to be $E_{q_{n}(x_{1:n})}[f_n^2(x_{1:n})w_n^2(x_{1:n})] - \mu_n$. ($\mu_n$ is the actual expectation we want to estimate). This means in order the variance to increase, $E_{q_{n}(x_{1:n})}[f_n^2(x_{1:n})w_n^2(x_{1:n})] - \mu_n$ must increase as a function of $n$. But as a mathematical expression I am not able to show how this is guaranteed to increase with $n$: This is a very general expression which is dependent on the form of the function $f_n$, the proposal distribution $q_n$ and original distribution $\pi_n$. So I cannot see how the variance increases with the dimension $n$ here. Am I doing something wrong in my approach and calculations or am I missing something important here? Please help me on this issue. Thanks in advance.


title: Data Mining Methods with Nested Data

text: I have three general questions that I am really struggling to answer: a) When we have nested data (e.g. employees nested in departments which are nested in companies which are in specific geographical regions) how data mining methods adjust to nesting? b) Are data mining methods ""immune"" to nested data because of their goal (i.e. predictive analysis)? c)Can we use data mining methods (e.g. neural networks, random forests) as precursors of variable selection before we run multilevel analysis? Please accept my apologies if the general orientations of these questions are inconvenient to you. Thank you VERY much for your help. I am looking forward to receive your wise views.


title: What's the difference between correlation and simple linear regression?

text: In particular, I am referring to the Pearson product-moment correlation coefficient.


title: Dealing with multicollinearity by smoothing coefficients?

text: I would like to perform a linear regression. However, the predictor variables are families of variables indexed by time. Let's say the regression problem is: target ~ x(1)+x(2)+...+x(40)+y(1)+...+y(40) and imagine the x(t) and y(t) are measurements in 40 successive years. The variables inside the families x(t),y(t) are thus highly collinear. The usual way to deal with this problem would be variable selection. However, in the problem at hand it is much desired that the coefficients the x(t) resp. y(t) variables are chosen smoothly in time, so setting some of them to 0 is not a good choice. My idea was to choose the coefficients for the variable families as the values of penalised splines with knots every 5 years. Is there a way to do this in R?   Or is there another way of smoothing the coefficients without the use of splines?


title: How to get the required data?

text: The input is a bunch of numbers. How can I calculate a threshold that there are 10% of the totally amount of numbers are above this threshold. I think this is equally as to calculate the distribution function, but what I need is the threshold number, not the


title: "Philosophical" data contamination question

text: I'm going through Yaser Abu-Mostafa's Learning from Data course, and I'm having some trouble getting my head around data contaimination. So we know stuff like VC-analysis can give us some gaurentee that hypotheses in a set $\mathcal{H}$ will generalize, but the point was made that $\mathcal{H}$ must be chosen without examining the data. Intuitively this makes perfect sense. Eg, if I look at my data, notice some random 17th order polynomial $p$ fits it, and choose $\mathcal{H} = \\{p\\}$, we'd have very strong generalization gaurentees even without something more involved like VC- analysis. But clearly this isn't _really_ going to generalize. I'm just wondering where in the argument for generalization this is happening though. Ie, where in VC-analysis arguments do we assume that $\mathcal{H}$ was chosen without looking at our data?


title: Developing a statistical test to ascertain better "fit"

text: In a data set with thousands of data points, I am testing different short-term and longer term data outputs based on 5 rolling data points all the way to 100 rolling data points (which each value being a separate column in excel: 5, 6,..., 100). The test I developed (with rudimentary knowledge of the whole thing) is to check which are the better "fits" of these rolling data outputs (kinda like moving averages, but not quite), in other words, if column 5, 17 or 98. In theory, the better fits should produce lower standard deviations, in the conventional calculation. An analogy for better understanding is something like this: if I do a 5 point moving average on a five period sin wave, the output should be zero or flat, the same for a 6 point moving average on a six period sin wave, and so on all the way to 100. This is the same as saying that the standard deviation in each case should be zero. So the test really is choosing the lowest standard deviation in each column from 5 all the way to 100, which should then provide the best "fit" (to go back to the analogy, fit with the "sin wave"). However, as the statistical inclined will quickly note (am afraid am not one of them for sure), this method presents two major difficulties:   1. The longer data sets, ie closer to 100 rolling data points, have larger absolute deviations, thus producing higher standard deviations. So a lower standard deviation in the shorter data series may not be indicative of better "fitting", just lower absolute deviations;    2. Each column has a different number of rolling data points, so once again the numbers and statistical fit may not be comparable between the columns.  So my question is therefore: is there a statistical measure better than the standard deviation I used above, such that it is not sensitive to absolute deviations and number of data points? Thanks!


title: Calculating significance of difference between two binary values on one dataset

text: Suppose my sample comprises 350 instances. Each instance has two independent binary (`pos` and `neg`) evaluations, **A** and **B**. If my alternative hypothesis is that A is more often positive than B, how would I express whether the difference is statistically significant? * * * To be more precise on the data: This question came to mind when looking at patient blood serum data. These are 350 blood samples positive against **X**. Each sample may or may not also be positive against A or B, whereas A and B are both more often positive if X is positive, but biologically are not dependent on one another. So my initial thought was to look whether A is more often positive than B. In the example it was, A is positive in ~80% and B in ~20%. So I thought about assessing this difference for significance and this is where I got stuck. :)


title: How to calculate forecast error (confidence intervals) for ongoing periods?

text: I often need to forecast for future periods in monthly series of data. Formulas are available to calculate the confidence interval at alpha for the next period in the time series, but this never includes how to treat the second period, and third, etc. I'd visually imagine that if any forecast was graphed with upper and lower confidence intervals, generally those intervals should exponentially increase or decrease against the mean forecast, as uncertainty is a cumulative force. Let's say I had unit sale of Apr = 10 May = 8 June = 11 July = 13 and no other context such as seasonality or population data We need to forecast (albeit blindly) August, September, October. What method would you use? and more importantly here, how will you measure confidence for September and October? Sorry that this might be a simple question for some experts - I have been digging far for a clear answer, and I'm sure this is something all amateurs like me would love to understand.


title: What are the metrics to evaluate a machine learning algorithm

text: I would like to know what are the various techniques and metrics used to evaluate how accurate/good an algorithm is and how to use a given metric to derive a conclusion about a ML model.


title: Skills hard to find in machine learners?

text: It seems that data mining and machine learning became so popular that now almost every CS student knows about classifiers, clustering, statistical NLP ... etc. So it seems that finding data miners is not a hard thing nowadays. My question is: What are the skills that a data miner could learn that would make him different than the others? To make him a not-so-easy-to-find-someone- like-him kind of person.


title: Why is the logistic regression cost function scaled by the number of examples?

text: I sometimes see that the cost function, along with the regularizer is divided by 1/2m where m is the number of examples. When we are trying to find the minimum of the cost, why does scaling by this amount matter? It doesn't affect where the minimum is. Example: ![enter image description here](http://i.stack.imgur.com/V54zN.png)


title: What does that mean that two time series are colinear?

text: I am familiar with the concept of cointegration. But I hear sometimes people talking about colinearity (or collinearity) for time series. A set of points is collinear if they are on the same line. But what does that mean for time series? Is it exactly the same as cointegration of order 1? Or is there something stronger/different in the concept of collinearity?


title: Hypothesis Testing

text: > A Lab has been asked to evaluate the claim that drinking water in a local > restaurant has a lead concentration of 6 parts per billion (ppb). Repeated > measurements follow a normal distribution and the population standard > deviation is taken to be 0.25 ppb. α = 0.01. > > • A sample of three measurements is taken and finds: 6.79; 6.13; 7.17 > >   1. Is there evidence to suggest that the lead concentration is different > from 6 ppb? > This is a homework question that I'm really stuck on, so if someone could point me in the right direction that would be great. I realise that it is a hypothesis testing question and that the null hypothesis is that µ = 6 ppb while the alternative is that µ does not equal 6 ppb. However, I really don't know what to do with the three sample measurements. **Do I get the mean of the three figures or do I use each in a separate hypothesis test?**


title: Random Forest - how to know if variables affect positively or negatively

text: I'm running a `RandomForest` in R on a set of data with many variables. Using `varImpPlot()` I know how important is each variable to explain the target variable. However, I would also like to know whether a **_variable affects positively or negatively the target value_**. For example, in predicting personal income, I may get that "age" explains a lot, but I would like to know whether "higher age -> higher income", or whether "lower age -> higher income". Is there any way to know this with Random Forest?


title: Bayesian inference over an unknown variance

text: I am observing a random variable $X \in \mathbb{R}$ which can be assumed to be normally distributed with mean $\mu$ and variance $\sigma^2$. I am interested in fitting a posterior distribution over the **unknown variance** which, according to wikipedia, can be given in closed form by an inverse-gamma. My question is, given samples $\\{x_1,...,x_n\\}$ how can I fit this distribution, that is, calculate the parameters $\alpha$ and $\beta$? **edit:** I found this pdf, which gives the equations for the parameter updates, and much more relevant info.


title: Visualise data with different sample sizes

text: I did a literature review and ended up with a lot of data in the kind of "10 out of 20 studies found x while 4 out of 8 studies found y". To be more precise; I looked up the contribution of seafood to environmental impact categories in different studies. Here is a short extraction of my data: ![enter image description here](http://i.stack.imgur.com/orXjw.png) I would like to represent the data somehow now. For example I would like to have a graph that shows how many cases of "Fishing" contributes "33%" compared to "Processing", "Packaging" and so on. How can I deal with the different sample sizes? Should I select a minimum number of studies to include and if yes, is there any criteria to set that limit? I don't need to do statistics but would like to have some simple graphs to allow readers a very quick first overview.


title: Plotting a heatmap based on clustering in R

text: id       sam1  sam2  sam3  sam4      gene1       33    23    88    98      gene2        0     0    99    95      gene3       77   100    44    65      gene4        0     0     0     0      gene5      100   100   100   100       :       :      gene20000   58    33    78    56      I have 20K genes (rows) and four samples (columns) which is given in percentage. Each gene is common to four samples. The hypothesis is, if any of gene has 100% in any of sample it is said to be present in that sample and if 0% it is absent. (For eg. gene3 is present in sample2 (sam2) and gene4 is absent in all samples.) Since I have 20K genes, I would like to do clustering technique. Can I do biclustering and produce a heatmap? If so, please let me know how to do in R?


title: Geometric interpretation of penalized linear regression II

text: An older question gives an intuitive explanation of how penalized linear regression works, using two separate contours: one for the least square objective, one for the penalty term (i.e. regularization). It also says that the sum of the least squares and regularization objectives is minimized when the contours meet for the first time. It seems, however, that to "calculate" the point where they meet "for the first time", we need to know the rate of growth of both contours. Unfortunately, I can't figure out how to arrive at it. **UPDATE:** At the bottom of page 5 of these notes on variable selection and regularization, the authors say that $$ argmin \sum{(y_i - \alpha - \beta'x_i)^2}$$ subject to $$\sum{|\beta_j|} < t$$ is equivalent to minimizing $$\frac{1}{2n}\sum{(y_i - \alpha - \beta'x_i)^2} + \lambda \sum{|\beta_j|}$$. If we look at the contour growing from this point of view, it means that:   1. We fix the contour plot for the penalty term, somehow (?) calculating $t$ from $\lambda$.   2. We only grow the contour for the least square objective.  Is this reasoning correct? If yes, could you point me to some literature that discusses this in more detail, especially:   1. The equivalence between the minimization of the function with inequality constraints and the minimization of the objective function in the more familiar form with the $\lambda$ parameter.   2. The calculation of $t$ given $\lambda$.


title: How could a Tukey HSD test be more signif then the uncorrected P value of t.test?

text: I came by the post "Post-hoc Pairwise Comparisons of Two-way ANOVA" (responding to this post), which shows the following:               dataTwoWayComparisons <- read.csv("http://www.dailyi.org/blogFiles/RTutorialSeries/dataset_ANOVA_TwoWayComparisons.csv")          model1 <- aov(StressReduction~Treatment+Age, data =dataTwoWayComparisons)     summary(model1) # Treatment is signif          pairwise.t.test(dataTwoWayComparisons$StressReduction, dataTwoWayComparisons$Treatment, p.adj = "none")     # no signif pair          TukeyHSD(model1, "Treatment")     # mental-medical   is the signif pair.      (Output is attached bellow) Could someone please explain why the Tukey HSD is able to find a significant pairing while the paired (unadjusted pvalue) t-test fails in doing so? Thanks. * * * Here is the code output               > model1 <- aov(StressReduction~Treatment+Age, data =dataTwoWayComparisons)     > summary(model1) # Treatment is signif                 Df Sum Sq Mean Sq F value    Pr(>F)         Treatment    2     18   9.000      11 0.0004883 ***     Age          2    162  81.000      99     1e-11 ***     Residuals   22     18   0.818                           ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      >      > pairwise.t.test(dataTwoWayComparisons$StressReduction, dataTwoWayComparisons$Treatment, p.adj = "none")                  Pairwise comparisons using t tests with pooled SD           data:  dataTwoWayComparisons$StressReduction and dataTwoWayComparisons$Treatment                    medical mental     mental   0.13    -          physical 0.45    0.45            P value adjustment method: none      > # no signif pair     >      > TukeyHSD(model1, "Treatment")       Tukey multiple comparisons of means         95% family-wise confidence level          Fit: aov(formula = StressReduction ~ Treatment + Age, data = dataTwoWayComparisons)          $Treatment                      diff         lwr        upr     p adj     mental-medical      2  0.92885267 3.07114733 0.0003172     physical-medical    1 -0.07114733 2.07114733 0.0702309     physical-mental    -1 -2.07114733 0.07114733 0.0702309          > # mental-medical   is the signif pair.


title: Tukey for GLM can't find data in model

text: I fit a GLM to a dataset. Now I want to see where the difference between my groups is, so I tried to run a Tukey HSD as a post-hoc test. Because of it is a GLM, I can't use `TukeyHSD()`. So I tried to run a Tukey test with:               summary(glht(my.mod, mcp(treatment="Tukey")))        When I tried this, I get the following error:               Variable(s) ‘treatment’ have been specified in ‘linfct’ but cannot be found in ‘model’!      This is what it looks like in R:               > recol        Sample At..Ind. At.Sp.   X treatment     1      C1        1      1   0         C     2      C2        2      2   0         C     3      C3        1      1   0         C     4     FS1        1      1  50        FS     5     FS2        4      3  80        FS     6     FS3        0      0  80        FS     7     FL1        3      3 280        FL     8     FL2        4      3 310        FL     9     FL3        2      2 220        FL     10    VS1        2      2  40        VS     11    VS2        0      0  50        VS     12    VS3        1      1   5        VS     13    VL1        3      3 400        VL     14    VL2        7      6 600        VL     15    VL3        3      3 500        VL     16    MS1        1      1  30        MS     17    MS2        3      3  30        MS     18    MS3        0      0  30        MS     19    ML1        2      1 400        ML     20    ML2        7      7 300        ML     21    ML3        7      7 300        ML     >  my.mod=glm(Arachnida~as.factor(treatment),data=soort)     > summary(my.mod)          Call:     glm(formula = Arachnida ~ as.factor(treatment), data = soort)          Deviance Residuals:          Min       1Q   Median       3Q      Max       -2.3333  -0.3333   0.0000   0.0000   2.6667            Coefficients:                              Estimate Std. Error t value Pr(>|t|)        (Intercept)             3.107e-16  6.901e-01   0.000  1.00000        as.factor(treatment)FL -1.277e-16  9.759e-01   0.000  1.00000        as.factor(treatment)FS -3.511e-16  9.759e-01   0.000  1.00000        as.factor(treatment)ML  3.000e+00  9.759e-01   3.074  0.00825 **     as.factor(treatment)MS  3.333e-01  9.759e-01   0.342  0.73775        as.factor(treatment)VL  2.333e+00  9.759e-01   2.391  0.03141 *      as.factor(treatment)VS  3.333e-01  9.759e-01   0.342  0.73775        ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for gaussian family taken to be 1.428571)              Null deviance: 48.571  on 20  degrees of freedom     Residual deviance: 20.000  on 14  degrees of freedom     AIC: 74.571          Number of Fisher Scoring iterations: 2          > glht(my.mod, mcp(treatment="Tukey"))     Error in mcp2matrix(model, linfct = linfct) :      Variable(s) ‘treatment’ have been specified in ‘linfct’ but cannot be found in ‘model’      > summary(glht(my.mod, mcp(treatment="Tukey")))     Error in mcp2matrix(model, linfct = linfct) :      Variable(s) ‘treatment’ have been specified in ‘linfct’ but cannot be found in ‘model’       Could anyone help me with my problem?


title: Power calculation; determine required difference in means

text: I want to perform a power analysis for a 2 sample t-test for the difference in means. However, instead of computing the required sample size, I would like to compute the required difference in the means given the sample size. The sample size is fixed at $n=34$ (i.e. 17 per group). Due to different reasons, it is not possible to increase the sample size. Therefore, I was asked to perform a power analysis to determine the required difference in means. The investigator wants to determine if it would be worthwhile to even do the study or if the required difference is just too big and unrealistic to achieve. For a given type I error and power I can compute the required effect size. If I further assume the standard deviations in both groups (from prior studies) I can easily compute the required difference. However, I'm wondering if this is reasonable or whether there are problems that I'm missing with this approach.


title: Logistic regression for questionnaires on 7 points Likert scale

text: I want to analyze some data of a questionnaire on presence (I deal with experiments on virtual reality). The questionnaire was provided twice to particpants performing the experiment, after each of the two provided conditions. The 6 questions of the questionnaire were evaluated by participants on a 7-points Likert scale. Following what I found in bibliography, I have to analyze the data in the following way in order to find if the differences between the two conditions are significant: I have to count the number of answers that have a score of 6 or 7, then I will have for each condition a variable equal to the mean of the count of 6 and 7 scores among all the 6 questions. Afterwards, I have to treat those two variables as binomially distributed for a logistic regression on group. My problem is that I don´t know how to perform the logistic regression. I tried to study but I do not understood how to apply it to my case. In addition I have also problems in undestanding how to perform it in R. Do you also have an example in R? Let's say that the two variables are (meand +-std):               Count_condition1 = 1 +- 1.7 and Count_condition1 = 2 +- 2.0


title: How to rigorously test if variables are drawn from a certain distribution

text: Say I have a list of numbers $X = \\{x_1, x_2, \dots, x_n\\}$, and I expect them to be drawn from a certain distribution. For my case it is the Binomial distribution $P(x) = \binom{n}{x}p^x(1-p)^{n-x}$, but I think I general answer would be the most helpful. What is the standard and most rigorous way I could determine if the samples $X$ are from that distribution?


title: to test the difference between 3 vaules

text: I measured half maximal effective concentration (EC50) of three different RNA construct and I want to test whether they are statistically significance or not. Which test should I use?


title: Trend detection quantitative models

text: What are the best quantitative models for trend detection? I.e. market trend.


title: Regarding the quantile regression via optimization approach

text: The quantile regression is defined through the optimization approach. But I am not clear how does the function of $\rho_{\tau}(u)$ related to the $\tau$-th quantile. Or in other words, how to derive the equation of the $\tau$-th sample quantile $=argmin \sum\rho_{\tau}$ ![enter image description here](http://i.stack.imgur.com/DuRiO.jpg)


title: Interaction in logistic models with R. Using * operator or create an interaction variable?

text: I want to test the presence of an interaction term in a logistic regression with glm(). The formula is:               Gest.hypertension ~ ART.conc + Parity + Age + Smoke + Nulliparity + Syst.disease.type      With Age the only continuous variable and the others categorical. The interaction term is ART.conc (2 leves) x Parity (3 levels), which are also the variables of interest. I found two ways of representing the interaction in glm();   * By explicitly defining ART.conc * Parity, which create both simple and interaction effect.   * By creating an artificial new variable with interaction(ART.conc, Parity) and adding the simple terms. This are the results, omitting the other controls variables (OR and p value presented):               model1: Gest.hypertension ~ ART.conc * Parity + ...          ART.conc->yes                       3.35     0.073 .     Parity->2                           7.15     0.001 **     Parity->3                             38     <0.001 ***     ART.conc->yes:Parity->2            0.262     0.054 .     ART.conc->yes:Parity->3           0.0532     0.009 **      * * *               model2: Gest.hypertension ~ interaction(ART.conc, Parity) + ART.conc + Parity + ...          interaction(ART.conc, Parity)->yes.1  3.35    0.072590 .       interaction(ART.conc, Parity)->no.2   7.15    0.001206 **      interaction(ART.conc, Parity)->yes.2  6.26    0.002698 **      interaction(ART.conc, Parity)->no.3   38.05    0.000202 ***     interaction(ART.conc, Parity)->yes.3  6.78    0.011071 *       ART.concyes                              NA         NA         Parity2                                  NA         NA         Parity3                                  NA         NA      so the two models give more or less the same results one one of the two terms is zero and that's ok. I have troubles understanding what happened when the two terms are different from zero. For example in the second model that keeping parity == 2 and passing ART.conc from no to yes, we pass from an OR of 7.15 to 6.26. Instead, applying what I understood of logistic regressions and interactions, the coeff for the case (ART.conc, Parity)->yes.2 should be $$OR_{parity.2}*OR_{art.yes:parity.2} = 1.87$$ very different from 6.26... Therefore, I'd like to know what I'm missing and which of the model should I use. Thanks! UPDATE!!! Thanks to @gavin simpson answer I thought I solved the problem, and ditched the interaction() trick. But then I noticed a pretty silly error by me. We said that to interpretate the interaction you do $$OR_{simple.eff}*OR_{interaction.eff}$$ Taking in example the case ART.yes/Parity.2. I said that to compute the effect of interaction of ART on Parity with the classical method you should do $$OR_{art.yes}*OR_{art.yes:parity.2} ≈ 3.35*0.262 ≈ 0.87$$ But if you perform $$OR_{parity.2}*OR_{art.yes:parity.2} ≈ 7.15*0.262 ≈ 1.87$$ I was troubled by this disparity. I need a number that could describe the overall risk of the pathology in the above mentioned case, and instead i got two numbers that i didn't know how to describe. Trying to solve this I noticed my error. To compute the overall effect, is not enough to multiply one simple effect for the interaction effect, but you should add the other simple effect too! $$OR_{simple.eff.1}*OR_{simple.eff.2}*OR_{interaction.eff}$$ or $$OR_{art.yes}*OR_{parity.2}*OR_{art.yes:parity.2} ≈ 3.35*7.15*0.262 ≈ 6.26$$ 6.26 is the same OR for the `interaction(ART.conc, Parity)->yes.2` that is produced using the interaction trick. Therefore, using interaction(), with `drop=T` option and not putting at all the simple effects by themselves (to avoid the NAs), you have that model 1 and model 2 give the same results, with model two including the already computed effect of the interaction.


title: Literature for Cross Validation on Sparse Data?

text: I've read a lot about Cross Validation to estimate prediction error, specifically for selecting the number of components in a PCA model (I'm not doing SVD/PCA, but it's very similar), but I can't find much literature on using CV for really, really large sparse data. Could anyone help me out with paper links, personal experience, or other techniques for estimating prediction error? Also, are there any special concerns for performing CV on Poisson type data?


title: how do you do regression analysis on advertising impact

text: i have data that includes clicks, spend, signups and date. for 1 week, i turn off advertising spend to see what clicks and signups are. the next week, i turn advertising back on to see what the new clicks or signups are. Given this 2 sets of data, how can i run regression analysis to see how impactful advertising is? Should i run a regression analysis on y=(signups_week2-signups_week1) and x=(spend_week2-spend_week1) ? Thanks, J


title: Are explanatory variables considered random in PCA?

text: One of properties of PCA states that the sum of the variances of the principal components is equal to the sum of the variances of the explanatory variables. I wonder how to interpret this as I've always thought that we do not consider $X$'s as random variables. I'm quite new to probability theory and I need to get it straight: **are explanatory variables random variables (or do we consider them fixed)?** And if we don't consider them random, how is it possible to apply the variance operator to a variable that is not stochastic?


title: How to graph quantitative and qualitative improvement in human health related measures?

text: I'm creating an app that logs daily data from my users. Data includes:   * waist circumference, weight (Quantitative data)    * certain conditions like psoriasis, migraines, fatigue (Qualitative data rated from 1 to 5).  Showing charts for each entitiy is fairly easy.   * How can I take both data and show an overall improvement graph/chart based on each day?    * Which type of graph/chart will work in this situation? I'll be using google charts to visualize my data (i.e., see this this gallery of Google Charts).


title: Conditional expectation of $X$ given $Z = X + Y$

text: Suppose I have two independent normal variables $X$ and $Y$ with known mean and variance. Defining $Z = X+Y$, what is the most straightforward way to compute $\mathbb{E}\left[X|Z\right]$? I am writing \begin{align} \mathbb{E}\left[X|Z\right] & = \int_x x f_{X}(x|Z=z) \mathrm{d}x \\\ & = \int_x x \frac{f_{XZ}(x,z)}{f_Z(z)} \mathrm{d}x \\\ & = \int_x x \frac{f_{Y}(z-x)f_X(x)}{f_Z(z)} \mathrm{d}x, \\\ \end{align} but I am not sure if the best (and only?) way is to compute this expression using our knowledge of the pdf's. Thanks.


title: How to calculate probability for n dependent events

text: I have a dataset containing some approx. 11 million instances, all of which can be labelled as either class A or class B. I know a priori that approx. 1,000 of these instances belong to class A, and the rest are class B. However none of the instances are actually labelled. So, I would like to be able to sample these 11 million instances, and know the probability of my sample containing an instance from class A. For instance: If I take 30,000 random samples without replacement from the 11 million instances (which are i.i.d), what is the probability that my sample will contain an instance from class A (size 1,000 instances)? Is this even possible to calculate? Also is it possible for me to approximately bound how many instances from class A would be in this sample (for varying sample sizes)? Is there a formula I can simply plug my values into? Apologies if this problem seems simple, but stats has never been a strength of mine. This isn't a homework task, but a work problem. I'm trying to build a training set from a domain outside my own, whereby I need to minimise the prevalence of class A in my training set. Class A are positive examples which need to be removed - unfortunately when the data was first processed the data labels were discarded :-/ . It would take years to label them manually and even that process is prone to human error! Thanks. Update: For those with a similar problem I found methods in the Apache Commons Mathematics Library for calculating binomials and more importantly the Hypergeometic Distribution which be of use. The hyper() command available in R (suggested by @whuber) also proved very useful. Would certainly recommend. For instance to produce a plot you can use: plot(dhyper(0:10,1000,11000000-1000,30000),type="b",xlab="x label",ylab="y label", main="Title") Where 0:10 is the range. 1000 is the number of positive marbles in the urn. 11000000-1000 is the number of negative marbles in the urn. 30000 is the number of draws. This will produce a plot with lines and points which is type "b". Other types are described online (I would post more links, but I can't with such low reputation). Thanks again for the help everyone.


title: Ensemble model performs better with worse performing consitutent models?

text: I have a forecast model I am developing that uses some very unreliable input data, missing data (due to sensors or comms failures) is the rule, not an exception. The quantity being forecast is a daily average over a region, so the average of all hourly measurements from all sensors. This forecast is made just before each day boundary by a human forecaster, assisted by a model (I am working on an improvement to our existing model). Now, the main part of the model takes recent history of daily averages and uses some fairly standard time-series forecasting techniques to predict the next day. This works okay, however it is well known the that last few hours of each day are particularly strong predictors of the next days value and the human forecasters have been using this knowledge for some time. In order to incorporate that technique into the model, I take the daily value time-series prediction and combine it with the last few hours of data for the day just ending, then train weights for each of these for a weighted average ensemble prediction (I optimise over RMSE). With the weight for the daily model fixed at unity, the sum of the hourly model weights typically come to around 0.5-0.8 or so. I retrain the weights each day using only the most recent 100 days of data because there are biases in each sensor (due to the underlying physics of the system as well as sensor issues) that change over time. The reason I used a weighted average and not some other much more powerful technique is due to the missing data. Most days have about 5-10% of sensor data missing when the forecast is being made (but not the same 5-10% each day). I figured a weighted average should be robust in this case, with the weights (hopefully) corresponding to the relative predictive power of each model. Now, this worked pretty well however I thought I could do better. Knowing that the hourly data has unique (but unknown a priori) biases, I thought I could fit a simple linear model for each hour for each station, using the last 100 days of data for that station and hour as the independent variable and the daily average for the following day as the dependant variable. Then I'd take the predictions from each of these very simple models, include the daily model based predictions and take a weighted average of those. In other words instead of just using the raw hourly data, I'm using a model trained using that data. This would surely produce a better result. But no, the weighted average forecast performs more poorly in this case. I simply cannot understand why. I checked and for an un-weighted average, using the predictions from the linear models produces a better result (as you'd expect) but once I train the weights, the model using just the raw data performs better. Interestingly I find that in the linear model predictors case, more of the weights (around 60-70% of them typically) actually get fixed at zero in the optimisation, whereas for the raw data approach less than 50% get fixed to zero (I am using L-BFGS-B minimisation [in scipy in anyone is interested]). It was actually that fact that so many got fixed to zero in the raw data case that lead me to try the linear fit approach. I figured that the badly biased sensors were being ignored and if I fixed the bias then perhaps that would be able to contribute more meaningfully. I'm not sure if I'm just doing something a bit silly, or if there is a fundamental concept that might explain this. When might you expect an ensemble forecast to perform better with worse input models?


title: How to obtain the variance of my dependent variable in a linear regression with R?

text: How can I obtained the estimated variance of a linear model when using R, i.e. \begin{equation} \widehat{var(y)}. \end{equation}


title: Understanding a measure of convergence of MCMC simulations

text: I am trying to better understand better the Gelman/Rubin measure of convergence of MCMCs. The method starts off by defining two quantities: $B$ and $W$. $B$ is said to be the between chain variance (as typically $m$ chains are run in parallel), and $W$ is the average within chain variance. I understand $W$ - it is merely the average of all the variances in each of the individual chains. $B$, at first glance, appears to measure the average deviation across all chains of their average from the overall mean. However, $B$ is multiplied by $n$ - the sample size - which I don't quite understand. Texts say that this is because there are $n$ points in each chain, however, I was hoping someone could provide a little intuition here? Further, $W$ is combined with $B$ in a weighted average, with $W$'s weight being $(n-1)/n$, and $B$'s weight simply being $(1/n)$. This is then typically divided by $W$, and the square root taken. If $B$ equals $W$ I can see that this ratio is then $1$. However, I don't see how this ratio ever gets far away from $1$ for moderately large sample sizes. It must be the case that $B>>W$, part of which I can understand if the chains are initially overdispersed. Still, the fact that $B$ is multiplied by $(1/n)$, makes it look like this part of the weighted average should disappear at a fast rate, even if $B$ is large. Basically, I think my question could be condensed into the following 'Why is this measure of convergence reasonable? Given that as $n$ gets large, ceteris paribus, the role played by the between chain variance appears to become smaller.'


title: Are time series variances additive?

text: I am trying to measure and quantify risk, variance, and standard deviation over a time period $T$. It is broken into two sub-periods $t_1$ and $t_2$. $X_1$ is the time series for $t_1$, and $X_2$ is the time series for $t_2$, where $t_1+t_2=T$. We can assume the $X_1$ and $X_2$ are sub-time series of a longer time series. All I know is that the variances over the two sub time periods are, say $S_1$ for $t_1$ and $S_2$ for $t_2$; we do not know what $X_1$ and $X_2$ are. I want to find the variance over the entire time period $T$ if $S_1$ is the variance for time period $t_1$, and $S_2$ is the variance for time period $t_2$, and $t_1 + t_2 =T$. My rationale so far: If $S_1=S_2=S$, then the variances over $t_1$ and $t_2$ are $S_1t_1$ and $S_2t_2$ respectively. Since $S_1=S_2$, it simply becomes $S_1T$ as $t_1+t_2=T$. The standard deviation over the time period is $S\sqrt{T}$(squaring both terms and adding, then take the square root again). If $S_1$ and $S_2$ are not the same, then is the formula for standard deviation:   $$ s_T=\sqrt{S_1t_1 + S_2t_2 + 2\sqrt{t_1t_2}{\rm Cov}(X_1, X_2)} $$ or do I add like it is a portfolio, with result being:   $$ s_T=\sqrt{\frac{S_1t_1^2}{T^2} + \frac{S_2t_2^2}{T^2} + \frac{2t_1t_2}{T^2}{\rm Cov}(X_1, X_2)} $$ where $s$ is the standard deviation? I have been told to use square root of time to change time periods for standard deviation. However, if we think of the time series as a portfolio with weightings $t_1/T$ and $t_2/T$, then the weights are not the square root of time with respect to variance but just the straight time weights as in the portfolio risk formula: $$ s_p=\sqrt{\frac{S_1h_1^2}{H^2} + \frac{S_2h_2^2}{H^2} + \frac{2h_1h_2}{H^2}{\rm Cov}(X_1, X_2)} $$ where $\frac{h_1}{H}+\frac{h_2}{H}=1$ I hope this clarifies what I am asking. Can someone please explain to me the proper way to think about this? Do I add variances with time squared or just time?


title: Sum to Zero Constraint GAM Factor Interaction

text: I was under the impression that the smooths fit with `mgcv` were made identifiable through a sum to zero constraint - i.e. if one was to sum the smooth over the values of its covariates, it would equal zero. This question gives some background, with an excellent answer from @Gavin Simpson. With single smooths, this seems to be the case:               set.seed(1234);      dat <- gamSim(4);           mod1<-gam(y ~ s(x1)+s(x2), data=dat)      p<-predict(mod1,dat,type="terms")     sum(p[,1]) #[1] -8.729996e-16     sum(p[,2]) #[1] -4.956868e-12      But when a smooth is interacted with a factor, this is no longer the case. Is this an error in my understanding, or is something else going on?               mod2<-gam(y ~ s(x2, by=fac), data=dat)      p<-predict(mod2,dat,type="terms")     head(p)     sum(p[,1]) #[1] 1.923649     sum(p[,2]) #[1] 6.496321     sum(p[,3]) #[1] 45.36179


title: Metric for cross-comparing vectors of multiple values for a few class labels

text: There are 20 individuals who get tested for a particular task on 3 consecutive days. They are considered to be iid. These results are stored in a 60 element vector. There is a similar test on 3 consecutive days for a different task involving the same 20 people producing a 60 element vector. I would like to select which vectors are the most interesting. Any interesting features would be useful. What I have done is produce and aggregate (total sum) of each vector and rank them. Effectively there are 3 data points which there are multiple tests of. Each vector has comparable values and do not need normalization etc. Are there any recommendations of metrics on each vector that I can do to cross compare the values? I would not wish to do a very extensive statistical analysis, something relatively basic to extract the information at present.


title: When using higher order terms in multivariable regression what is the effect on the P-value

text: I have a large dataset, lets say 500 lines of data. As a quick workaround to higher order regression I've made educated guesses about what higher order terms I would expect in my regression equation. I've then included these as separate variables in my linear regression. For example:               x(1) = x(1)     x(2) = x(2)     x(3) = x(3)     x(4) = x(1)^2     x(5) = x(1)*x(2)      etc. My question is since I can get a t stat and hence a P-value for each of these x(1) - x(5) how do I judge the significance of x(1) since x(1) is a component of x(4) and x(5). Or can I treat them as separate?


title: Shall I trust AIC (non-full model) or slope (full model)?

text: The purpose to run regressions for butterfly richness again 5 environmental variables is to show the importance rank of the independent variables mainly by AIC. In non-full models, they reveal that variable A tends to be more influential than the others by delta AIC. However, in the full model, the regression coefficient of variable A is slightly second to that of variable B. (R-square of the full model is 1.43) The conflicting outcomes (non-full model by AIC and full model by slope) seems to make it difficult to ascertain that variable A is the variable mostly weighted. Please kindly suggest which criterion should be relied on for the specified purpose or any further test should be carried out. Thank you.


title: "Paired data" / "repeated measures" in 2x2 experiment

text: I have test subjects answering questions, the data will be quantitative. The questions can be divided to four categories by two factors that have two levels. For each question a factor A, B, A and B or none are "active". So each subject answers one of these four questions. I am interested in how the factors A and B affect the answer. So, this is resembles 2x2 factorial experiment, BUT the data are "paired", in the sense that each subject has answered all of the four questions. Do you have suggestions for statistical tests I can use? Do you know how situations like this has been dealt with before? EDIT: I found out that "repeated measures design" might be relevant.


title: Stacking up panel data to use OLS over it

text: I'm posting this question for someone. He has the model: (1) $Y_{it} = a + b X_{1,it} + c X_{2,it} + d X_{3,it} + e_{it}$ $i = 1,2,...,N;\phantom{...} t = 1,2,3,4,5,6,7$ Variables are approximately normal and we suspect no violation in OLS assumptions. The coefficient estimates are not subscripted with $i$ because the intention is to get panel-like estimates. $t$ represents a year and $i$ represents a company. Instead of using pooled OLS he wants to stack up the data. So, if we consider the matrix $\bf{Y}$ with number of rows $N$ and number of columns $7$ where $Y_{it}$ is some element in $\bf{Y}$ $\forall i \wedge \forall t$, he wants to stack this such that it's a $7*N$ row, $1$ column vector as the response variable to be estimated in usual cross- sectional OLS (after doing the same for the design matrix, of course). Is this bad? (Looking for general answers on its weaknesses, or whether it's fine).


title: What is the correct way to calculate the explained variance of each EOF as calculated from a gappy data set?

text: I am trying to determine the correct amount of variance explained by each mode of an Empirical Orthogonal Function (EOF) analysis (similar to "PCA") as applied to a gappy data set. (i.e., containing NaNs). The following question builds on an earlier one that I had regarding the differing results obtained from the decomposition of the data set's covariance matrix using either `eigen` or `svd`. In essence, the problem is that I have read that both decompositions can be used interchangeably for obtaining the EOFs from a square covariance matrix. This does seem to be the case when the data set is not gappy (as illustrated below):               ###Make complete and gappy data set     set.seed(1)     x <- 1:100     y <- 1:100     grd <- expand.grid(x=x, y=y)          #complete data     z <- matrix(rnorm(dim(grd)[1]), length(x), length(y))     image(x,y,z, col=rainbow(100))          #gappy data     zg <- replace(z, sample(seq(z), length(z)*0.5), NaN)     image(x,y,zg, col=rainbow(100))               ###Covariance matrix decomposition     #complete data     C <- cov(scale(z), use="pair")     E <- eigen(C)     S <- svd(C)          #sum of lambda     sum(E$values)         sum(S$d)     sum(diag(C))      The sum of lambda in both `eigen` and `svd` equals the sum of the diagonal of the covariance matrix. So far, so good - Both methods explain the correct amount of variance. The next example does the same routine for a gappy version of the data set (50% `NaN`s):               #gappy data (50%)     Cg <- cov(scale(zg), use="pair")     Eg <- eigen(Cg)     Sg <- svd(Cg)          #sum of lambda     sum(Eg$values)         sum(Sg$d)     sum(diag(Cg))      And here we see that the lambda values calculated by `svd` are greater than the sum of the diagonal of the covariance matrix. Those calculated by `eigen` are equal. However, because the covariance matrix is no longer positive definite, there are some negative trailing lambda values. In my previous question I showed that this tendency becomes greater with increasing gappiness.   So, I can live with this if need be, but now I'm concerned about how to correctly assign how much of the data set's variance is explained by each EOF. This should be lambda/sum(lambda). When I plot the cumulative explained variance of the EOFs, you will see the problem - because the `eigen` decomposition contains some negative eigenvalues, the slope of cumulative explained variance is steeper and bell-shaped:               #cumulative explained variance of the EOFs     E.cumexplvar <- cumsum(E$values/sum(E$values))     S.cumexplvar <- cumsum(S$d/sum(S$d))     Eg.cumexplvar <- cumsum(Eg$values/sum(Eg$values))     Sg.cumexplvar <- cumsum(Sg$d/sum(Sg$d))               ###plot the cumulative explained variance     png("cumexplvar.png", width=8, height=4, units="in", res=200)     par(mfcol=c(1,2))     YLIM <- range(c(E.cumexplvar, S.cumexplvar, Eg.cumexplvar, Sg.cumexplvar))          plot(E.cumexplvar, t="o", col=1, ylim=YLIM, xlab="EOF", ylab="cum. expl. var.", main="non-gappy")     points(S.cumexplvar, t="o", pch=2, col=2)     abline(h=1, col=8, lty=2)     legend("bottomright", legend=c("Eigen", "SVD"), col=c(1,2), pch=c(1,2), lty=1)          plot(Eg.cumexplvar, t="o", col=1, ylim=YLIM, xlab="EOF", ylab="cum. expl. var.", main="gappy")     points(Sg.cumexplvar, t="o", pch=2, col=2)     abline(h=1, col=8, lty=2)     legend("bottomright", legend=c("Eigen", "SVD"), col=c(1,2), pch=c(1,2), lty=1)          dev.off()      ![cumulative explained variance](http://i.stack.imgur.com/SXNL1.png) The problem may be that I should be using the sum of the absolute eigenvalues to assign their explained variance, but this also leaves me to wonder how to interpret the explained variance of the negative eigenvalues. I would be very grateful for any insight, as this is not an issue that I have come across in any reference regarding EOF as applied to gappy data.


title: Do you perform cross validation on bootstrap samples?

text: Let's say that I have to chose the best model to predict a variable but my sample size is small. I would like to resample my data using the bootstrap, run each model and evaluate its prediction error with cross validation on each bootstrap sample. In the end I will chose the model with the lowest average prediction error. Is this a acceptable procedure? If yes is there something I should be aware of while using this procedure?


title: Mixed repeated measures Anova

text: I'm trying to assess the influence of fertilisation and irrigation on biodiversity in a complete randomized block design. I have two treatments, one control and two blocks. Treatments went on for three years so I have three repeated measures. I want to find out if biodiversity is significantly different for the treatments; #1 in the single years and #2 over all three years. I think that repeated measures anova should do the job, however, I'm struggling on how to include the blocks in the calculation. Would it be appropriate to calculate a mean over the blocks before doing an anova? Is there another way to do this?


title: Inverse function for a non-decreasing CDF

text: For a CDF that is not strictly increasing, i.e. its inverse is not defined, define the quantile function $$F^{-1} (u) =\inf \\{x: F(x) \geq u \\},\quad 0<u<1\. $$ Where U has a uniform $(0,1)$ distribution. Prove that the random variable $F^{-1} (u)$ has cdf $F(x)$. In case of a strictly increasing CDF the proof is quite easy because the inverse is defined. Define $X=F^{-1} (u)$ $$ P\left[X<x \right]= P \left [F^{-1} (U) \leq x \right]= P \left[U \leq F(x) \right] =F(x) $$ But how do I accomondate the nondecreasing CDF whose inverse is given by the quantile function? I am a begginer so any help is welcome. Thank you.


title: Binomial testing for mutation patterns

text: I am looking for a certain type of mutation (DNA base A going to DNA base B when it is followed by C) in a dataset of 24 samples. I've noticed that this has happened 3 times out of the 4 times a gene has been mutated in the dataset (1 mutation per sample). I've also got counts of how many times that kind of mutation occurs across all the genes in the genome every time they've been mutated. Is the use of a binomial test comparing 3/4 to the background probability as defined by the number of A>B changes when followed by C in the entirety of the dataset/ the four samples in which that change has been observed a valid method to infer if the frequency of A>B mutations I've seen is greater than expected by chance? In R the code I'm looking to use would be               binom.test(c(3,1),p=(background))      where background = number of A>B changes in CX/ number of CX sequences that have undergone mutations for the four samples. Best wishes, Ankur.


title: Multiple imputation on new data in R

text: I am looking for a R package that can do multiple imputation on 2 sets of data in the same fashion. That is, I would like to multiply impute the training set and then impute the test set in the same way the training set was done. I have looked at `Amelia` and `mice` and others but can't figure out how to perform this process. The `caret` package has imputation where both the training and test set are imputed in the same way. I am looking to do this in the same way with multiple imputation. Thanks for any advice.


title: Multiple Imputation and two-level data

text: I have a question on multiple imputation where one variable is the sum of several sub-groups. I have about 5 variables with a significant level of missingness. However I have a sixth variable which is the sum of the 5 sub- groups and it is missing considerably fewer observations. To give you an idea of data structure:               sg1     sg2     sg3     sg4     sg5   total     NA      NA      NA      NA      NA    30000     5000    20000   10000   4000    6000  45000      Is there a specific method for imputing strata like this? I'm using R, so if you have any R specific advice I would greatly appreciate that as well. Thanks.


title: Validation of clustering results through correlation maps

text: How can I compare correlation maps independently to the number of clusters in terms of measuring the 'quality' of well separating (uncorrelated) clusters, i.e. a criterion to maximize the intra-cluster agreement and the inter-cluster disagreement?


title: What is the purpose of working on a logit scale in partial dependence plots?

text: What is the purpose of working on a logit scale in partial dependence plots (in binary classification)? One could simply go about as follows:   1. Grow a forest   2. Suppose `x` has `v` distinct values in the training data set. Construct `v` data sets as follows. For each of the `v` values of `x` make up a new data `se` where `x` only takes on that value, leaving all the other values untouched.   3. For each of the `v` data sets thus obtained, predict the response using random forest.   4. For each of the `v` data sets, average these predictions   5. Plot `v` and the corresponding averaged predictions Instead of just doing that, in literature they go on and replace step 4 and 5 by (see also partialPlot in R package randomForest): Step 4.Compute `average(log(predictions)-0.5(log(predicions) + log(1-predictions)))=average(0.5 logit(predictions))` Step 5.Plot `v` and the corresponding computed `1/2 logit(predictions)` Why?


title: Oversampling correction for multinomial logistic regression

text: When modeling rare events with logistic regression, oversampling is a common method to reduce computation complexity (i.e., keep all the rare positive cases but just a subsample of negative cases). After model fitting, adding a offset to the intercept term is a common method to correct the event probability to reflect the original sample proportion. The offset is equal to log( r1*(1-p1) / (1-r1)*p1 ), where r1 is the proportion of rare events in the oversampled data and p1 is the proportion in the original data. What is the equivalent formula with multinomial logistic regression, where 1 or more classes is oversampled?


title: cforest and randomForest classification prediction error

text: I used cforest and randomForest for a 300 rows and 9 columns dataset and received good (almost overfitted - error equal to zero) results for randomForest and big prediction errors for cforest classifiers. What is the main difference between these two procedures? I admit that for cforest I used any possible input parameters combination e.g. the best one, but still with big classification errors, was `cforest_control(savesplitstats = TRUE, ntree=100, mtry=8, mincriterion=0, maxdepth=400, maxsurrogate = 1)`. For very big datasets (about 10000 rows and 192 columns) randomForest and cforest have almost the same errors (the former slightly better on the same level as radial kernel svms), but for the mentioned small one for my surprise there is no way to improve cforest prediction accuracy...


title: Metric for probability based classification

text: I am doing a system for classifiying documents. The project demands the use of probability based output. So a sample will have a probability for belonging to each class. For now I use logistic regression, but this could be subject to change. So I don't want to do an R^2 approximation. I also don't want to use standard metrics for classification like F-measure because it doesn't work with probabiliies. I have no idea what is the custom metric in this situation. Any ideas?


title: Finding k-representatives from a 2 dimensional data set

text: I have a data set for second hand vehicle sales with 2 variables price and kms covered. If we plot these data in a 2-dimensional plane, the points will form scattered chart. I would like to select k samples (k - no of samples given by the user at run time) from the total data, maybe some 5,000 records. The samples selected must ensure that they are real representatives of the 5,000 records and the samples selected must be linked with the other underlying data. So when we select a particular sample, it should be possible to get the data for which the sample acts as representative. This may be done in more than one level also as done in subspace clustering. Random sampling, density biased sampling or sorting by attributes will not help. Because a vehicle with more kms run may be priced high because it is in good condition. And a vehicle with less kms covered may be priced lower because there may some problem with the vehicle. So, if there is any statistic method available to solve this problem please help me.


title: Essential data checking tests

text: In my job role I often work with other people's datasets, non-experts bring me clinical data and I help them to summarise it and perform statistical tests. The problem I am having is that the datasets I am brought are almost always riddled with typos, inconsistencies, and all sorts of other problems. I am interested to know if other people have standard tests which they do to try to check any datasets that come in. I used to draw histograms of each variable just to have a look but I now realise there are lots of horrible errors that can survive this test. For example, I had a repeated measures dataset the other day where, for some individuals, the repeated measure was identical at Time 2 as it was at Time 1. This was subsequently proved to be incorrect, as you would expect. Another dataset had an individual who went from being very severely disordered (represented by a high score) to being problem-free, represented by 0's across the board. This is just impossible, although I couldn't prove it definitively. So what basic tests can I run on each dataset to make sure that they don't have typos and they don't contain impossible values? Thanks in advance!


title: Matlab's quadprog equivalent in python?

text: For my SVM algorithm I need to do an optimization in standard QP form. In Matlab I would use 'quadprog' with the 'interior-point-convex' algorithm from the Optimization toolbox. What is an equivalent (preferably free for acadamic use) library for python?


title: Correlation of log-normal random variables

text: Given $X_1$ and $X_2$ normal random variables with correlation coefficient $\rho$, how do I find the correlation between following lognormal random variables $Y_1$ and $Y_2$? $Y_1 = a_1 \exp(\mu_1 T + \sqrt{T}X_1)$ $Y_2 = a_2 \exp(\mu_2 T + \sqrt{T}X_2)$ Now, if $X_1 = \sigma_1 Z_1$ and $X_2 = \sigma_1 Z_2$, where $Z_1$ and $Z_2$ are standard normals, from the linear transformation property, we get: $Y_1 = a_1 \exp(\mu_1 T + \sqrt{T}\sigma_1 Z_1)$ $Y_2 = a_2 \exp(\mu_2 T + \sqrt{T}\sigma_2 (\rho Z_1 + \sqrt{1-\rho^2}Z_2)$ Now, how to go from here to compute correlation between $Y_1$ and $Y_2$?


title: Strategy for analysing stock history data from multiple files with SPSS

text: I have a folder with a hundred comma separated files (CSVs) of which the filename equals a company stock exchange followed by it's symbol, delimited with a "_", eg: NASDAQ_MSFT.csv Each file contains historical daily stock information, eg. a csv file looks like:               Date,Open,High,Low,Close,Volume     29-Dec-00,21.97,22.91,21.31,21.69,93999000     28-Dec-00,22.56,23.12,21.94,22.28,75565800     27-Dec-00,23.06,23.41,22.50,23.22,66881000     ...     5-Jan-00,55.56,58.19,54.69,56.91,62712600     4-Jan-00,56.78,58.56,55.94,56.31,52866600     3-Jan-00,58.69,59.31,56.00,58.34,51680600      Now, I want to analyse this information (analyse the companies as if the company name is a column field rather than a new list of fields). But there are a couple of issues:   1. Each company has it's own file with the symbol and exchange as filename   2. Some companies start at different dates than others. Eg. some might have a range of 30 days while others have a range of 30 months (every step is still 1 day difference though). I am using SPSS as my analysis tool. My question is, how can I import these files to perform senseful analystical operations on them? For example, I wish to see the average slope of the open price of all companies together, etc.


title: R highfrequency package

text: I'm newbie to R, I'm using it for my thesis on trade classification algorithms. I've some problems with the package `highfrequency`. I've converted my txt data in `xts` format using the code:               library(highfrequency)           library(timeDate)          from="1990-11-01";          to="1990-11-02";          datasource = "~/raw_data";          datadestination = "~/xts_data";          convert(from, to, datasource, datadestination, trades=TRUE, quotes=FALSE,          ticker=c("ALO"), dir=TRUE,  extension="txt",header=FALSE,tradecolnames=NULL,     quotecolnames=NULL, format="%Y%m%d %H:%M:%S");          xts_data = TAQLoad(tickers="ALO", from="1990-11-01", to="1990-11-02",trades=T, quotes=F,datasource=datadestination)     head(xts_data)      but then when I load them on the console they appear like this without giving me                                   SYMBOL EX PRICE SIZE COND CORR G127     <NA>                ""     "" ""    ""   ""   ""   ""       1970-01-01 01:00:01 ""     "" ""    ""   ""   ""   ""       1970-01-01 01:00:00 ""     "" ""    ""   ""   ""   ""       1970-01-01 01:00:00 ""     "" ""    ""   ""   ""   ""       1970-01-01 01:00:01 ""     "" ""    ""   ""   ""   ""       <NA>                ""     "" ""    ""   ""   ""   ""        nothing! That's a problem I think because the next step I'd like to do is clean the data, I've tried with this code:               data("sample_tdataraw");          dim(sample_tdataraw);          tdata_afterfirstcleaning = tradesCleanup(tdataraw=sample_tdataraw,exchanges="N");          tdata_afterfirstcleaning$report;      but what I've to put `data("xts_data")`,`data("AC_trades")`?Because `sample_tdatarw` is generic,right? I'm using quite old data (144 stocks NYSE TORQ from Nov90 to Jan91), they are the only free TAQ data that I could find, they're taken from a zipped file found in internet, they originally were in ASC format, that I could not read so I've opened it with txt and I've saved in this format, maybe the problem is due to the format of the data? They're like this (date,time,exchange,etc) but I've tried even with (ticker,date,time,exc,etc)               19901102  9:31:59  N  18.87500   2000         0   6139      Additionally, when I run the package highfrequency I get this messagge:               The following objects are masked from ‘package:base’:              as.Date, as.Date.numeric      Thank you for the help!!


title: How to report effect size measures r and r-squared and what is a non-technical explanation of them?

text: If $r$ is the effect size for the correlation between $A$ and $B$, then $r^2$ is the amount of variance in $B$ that can be attributed to variable $A$.   1. Is it important to report both indexes in a report, or just one or the other?    2. How do you explain them in plain English (for a non-statistical audience)?


title: Restoring original distribution from noisy observations

text: There is a known set of pairs $(y_i, \sigma_i)$ such that $y_i = x_i + \sigma_i N_i$ $N_i \sim \mathbf{N}(0,1) $ for all $i$ $x_i \sim \rho$ for all $i$ where $y$ is observed value, $x$ is true value. How to restore the distribution $\rho$ ? It is not neccessarily normal.


title: Find matching samples via Canonical Correlation Analysis (CCA)

text: I want to use Canonical Correlation Analysis (CCA) to identify relationships between two sets of variables X and Y. The CCA should give a score (highest correlation) between two samples of X and Y. I tried to implement it via scikit-learn like that:               X = [[1, 0, 0], [1, 1, 0], [1, 1, 1], [2, 0, 0]]     Y = [[0, 0, 1], [1, 0, 0], [2, 2, 2], [3, 5, 4]]          cca = CCA(n_components=2)     cca.fit(X, Y)          for x in X:         print "----------"         print x         for y in Y:             print str(y) + " : " + str(cca.score(x,y))      The output for the first element in X:               [1, 0, 0]     [0, 0, 1] : 0.35461498401     [1, 0, 0] : -0.0502507710089     [2, 2, 2] : 0.0     [3, 5, 4] : -22.2417510911      But the result is not as expected, there is no highest correlation between [1, 0, 0] in X and Y. But the score returns something else back: "Returns the coefficient of determination R^2 of the prediction". How is it possible via CCA to find the highest correlation matching pairs (X#n, Y#m)? Is it possible via scikit-learn, or I have to use another library? Thanks in advance.


title: If $X, Y$ are jointly standard normal with correlation $r$, and $a, b$ are constants, what is $p(Y < b | X < a)$?

text: The application here is interpreting the correlation coefficient $r$ in terms of $X$'s ability to predict $Y$ for extreme values of $X$. For example, if $r = .8$, then what is $p(Y < 0 | X < -2)$?


title: how to perform factor analysis for clustered longitudinal binary data?

text: I measure a longitudinal binary outcome (correctness of detection, 0: incorrect, 1: correct) with respect to 5 different experimental conditions (1 baseline and 4 treatments). The outcome is always measured at the same 10 time points. Each of the 9 subjects participated in all 5 conditions. Additionally, for each subject and condition, the experiment was replicated 12 times. I therefore end up with 9*5*12=540 binary longitudinal series (= trials of 10 points each). My aim is to assess the influence of the experimental condition on my binary outcome. I need to build a model that would take into consideration the correlation along time for a given trial and the correlation among trials for a given subject. I am considering a 3 levels logistic models where 10 consecutive binary measurements (level 1) are obtained on replicates (level 2) which are clustered into subjects (level 3). My only level 1 covariate would be the time of measurement and as level 2 covariate, I consider the experimental condition. I don't have any level 3 covariate per se, but still want the model to consider between-subject variability. As there are many ways to look at those data, I'd be grateful to receive some input to know whether I'm on the right track and some insight as to which R package is most suited for the analysis I want to perform.


title: Mediation analysis using binary IV and difference scores

text: I have one IV ($A$) and two DVs ($B$,$C$). $A$ is a binary, experimentally manipulated variable. Each subject has two scores on $B$ and two scores on $C$, corresponding to the experimental manipulations $A$ and ~$A$. I would like to say that $A$ causes an increase in $B$, and increases in $B$ lead to an increase in $C$. $A$ should lead to increases in $C$, if and only if it leads to increases in $B$ (full mediation). Paired t-tests confirm $B$ is higher for $A$ vs ~$A$, and $C$ is higher for $A$ vs ~$A$. $\Delta B$ is highly correlated with $\Delta C$. I'm relying on difference scores under the assumption that $B$ and $C$ are not necessarily correlated because of individual differences in the baseline of $B$. How can I test for mediation here? Since $A$ is binary, there is no "$\Delta A$" with which to run a regression analysis. I'm not sure where to go next; any help is appreciated!


title: Computing chi-square for large tables with some expected cell counts equal to zero

text: What should I do if the expected value in a Chi-square goodness-of-fit test is zero? I know there's Fisher's test but I have a very large table!


title: What is the most appropriate way to test significance of effect difference across two groups?

text: Here is the question: After running two separate OLS regressions, one for each group (male vs. female), we get coefficient b1 for a specific IV in "male" regression, and coefficient b2 for the same IV in "female" regression. What would be the most appropriate way to test the hypothesis: Ho: b1>b2? As far as I know, Wald is a widely used test, but it only check whether the effect difference is different from zero (Ho: b1-b2=0) Thanks!


title: What method to test if the mean is constant?

text: I have a time series with prices (around 800 prices). I need to test if the mean is constant during all the series. I think I should subdivide the series in groups, but what method do I have to use to check if the mean is constant if the distribution is _not_ normal?


title: Does the product of two p.s.d kernel matrices result in a kernel matrix?

text: In a ML setting, where $a_1,..., a_n$ are a set of training points. A kernel function is a function $κ$ that gives the inner product between two vectors in the feature space: $κ(a_i, a_j ) = ψ(a_i) · ψ(a_j )$ where $ψ$ is a non-linear function. We can form an $n × n$ matrix $K$, called the kernel matrix, whose $(i,j)$ entry corresponds to $κ(a_i , a_j )$. Here, $K$ is symmetric and positive semi-definite and can be factorized as $K=XX^T$ where $X=U\lambda^{1/2}$ given that $K=U \lambda U^T$ is the spectral decomposition of the kernel matrix. Question:- Now given two kernel matrices $K_1$ and $K_2$, is $K_p=K_1K_2$ a kernel matrix?


title: R: adehabitatHS compana test - getting NaN for Wilk's lambda

text: I'm getting an "NaN" test result using the R package adehabitatHS, and the function "compana". The compana test compares habitat availability with habitat usage by animals. The test performs a compositional analysis (Aebischer et al. 1993), yielding a value for Wilk's lambda and a p-value. I'm running the test on empirical data for rattlesnakes (16 animals, 8 habitat categories). For some reason the compana test returns: lambda = NaN, p-value = 1.0. I've examined my data for the standard causes of NaNs, and I don't see any problems, such as missing values (NAs), negative proportions, or row sums > 100 percent. Moreover, the error arises when I use the full data set (16 rows, 8 columns), but not when I use a subset (missing one column). It does not seem to matter which column(s) is excluded, so long as only 7 or less are included. Details of my data are below. Any help or suggestions would be greatly appreciated. Platform: OSX 10.7.4 (Lion) R Version: R 2.15.1 GUI 1.52 Leopard build 64-bit (6188) R CODE: proportions for habitat used (matrix):               used.m         Atr   Cre    DW  Dist  Gra   Mar   Mes   Opn     1   2.76 19.31  4.14 17.24 0.00 28.28  0.00  0.00     2  15.12  1.16  3.49  5.81 2.33 32.56 25.58  3.49     3  57.35 13.24  1.47  0.00 6.62  0.00  8.09 12.50     4  41.18  0.98  2.94  0.00 3.92  0.00 29.41 18.63     5   9.52 17.01 23.13  0.00 0.00 37.41  3.40  0.00     6   2.54 67.01  0.00  0.51 0.51  6.09  3.05  0.00     7  10.45 15.67  5.22  5.97 0.00  0.00 33.58 10.45     8   6.34  0.00  0.00  4.23 0.00  0.00  0.00 21.83     9   0.66 18.54  1.32 39.74 0.00  8.61  2.65  0.66     10 20.75  0.00  0.00  0.00 0.00  0.00  0.00 49.06     11 20.14 18.06 12.50  8.33 0.00  5.56 22.22  9.72     12  0.00  0.00  0.00 58.04 0.00  0.00 13.99  0.00     13 44.23  4.81  0.00  0.00 0.00  0.00 28.85 19.23     14 13.04  2.61  1.74  2.61 0.00 19.13  6.09 25.22     15 35.96  3.51  0.00  2.63 0.00  1.75 25.44  0.88     16  0.00  0.00  0.00 35.14 0.00  7.21  0.00  8.11      proportions for habitat available (matrix):               avail.m         Atr   Cre    DW  Dist   Gra   Mar   Mes   Opn     1   1.11 39.18 15.50  0.29  0.00 30.18  4.15  3.43     2  18.96  1.02  0.87  7.40  4.36 37.80  6.17 10.56     3  65.75  3.85  0.57  0.00 12.36  0.00  9.97  3.43     4  20.52  5.44  2.01  0.00  0.08  0.00 10.99 44.22     5  11.16 43.55 10.34  0.00  0.00 18.11  9.60  3.17     6  25.24  9.88  2.77  2.44  4.95 12.90 11.44  2.19     7   5.47 21.13  6.67  0.57  0.00  0.10 30.73  9.59     8   6.23  0.00  0.00  4.85  0.00  0.00  0.00  5.72     9   0.85  1.97  0.09 33.52  0.00  5.58  1.73  0.84     10 18.97  0.00  0.00  0.00  0.00  0.00  0.00  5.59     11  0.16 44.41 19.16  3.53  0.00  1.28 21.65  4.13     12  6.67  3.89  0.31 13.63  0.00  0.00 20.56  0.00     13 21.40  2.38  2.66  0.00  0.00  0.00 21.23 46.96     14 42.62  9.36  1.49  3.89  0.43  3.61 20.79  5.80     15 33.32 12.90  4.17  4.09  0.15  3.96 12.32  4.94     16  5.65  0.00  0.00  1.36  0.00  0.59  0.00  9.71      load the adehabitatHS library               library(adehabitatHS)      run the compana test on full data set (all 8 columns)               fulldata.comp <- compana(used.m, avail.m)      check result               fulldata.comp$test     Lambda      P         NaN      1       now run compana on partial data (columns 1-7)               partdata.comp <- compana(used.m[, 1:7], avail.m[, 1:7])      check result               partdata.comp$test        Lambda         P       0.3319731 0.6280000


title: Assessing logistic regression models

text: This question arises from my actual confusion about how to decide if a logistic model is good enough. I have models that use the state of pairs individual-project two years after they are formed as a dependent variable. The outcome is successful (1) or not (0). I have independent variables measured at the time of formation of the pairs. My aim is to test whether a variable, which I hypothesized would influence the success of the pairs has an effect on that success, controlling for other potential influences. In the models, the variable of interest is significant. The models were estimated using the `glm()` function in `R`. To assess the quality of the models, I have done a few things: `glm()` gives you the `residual deviance`, the `AIC` and the `BIC` by default. In addition, I have calculated the error rate of the model and plotted the binned residuals.   * The complete model has a smaller residual deviance, AIC and BIC than the other models that I have estimated (and that are nested in the complete model), which leads me to think that this model is "better" than the others.    * The error-rate of the model is fairly low, IMHO (as in Gelman and Hill, 2007, pp.99):   `error.rate <- mean((predicted>0.5 & y==0) | (predicted<0.5 & y==1)`, at around 20%. So far so good. But when I plot the binned residual (again following Gelman and Hill's advice), a large portion of the bins fall outside of the 95% CI: ![Binned Residuals plot](https://lh5.googleusercontent.com/-DhQ3a9hTVoE/Tr1H-Csj_JI/AAAAAAAAAC0/eYXlUlkc6ic/s550/binned.res.jpeg) That plot leads me to think there is something utterly wrong about the model. Should that lead me to throw the model away? Should I acknowledge that the model is imperfect but keep it and interpret the effect of the variable of interest? I have toyed around with excluding variables in turn, and also some transformation, without really improving the binned residuals plot. **Edit:**   * At the moment, the model has a dozen predictors and 5 interaction effects.    * The pairs are "relatively" independent of each other in the sense that they are all formed during a short period of time (but not stricly speaking, all simultaneously) and there are a lot of projects (13k) and a lot of individuals (19k), so a fair proportion of projects are only joined by one individual (there are about 20000 pairs).


title: Rules of thumb for creating composite scales

text: Normally when I create an index I basically just take the or averages of the individual variables that comprise the index. But I'm wondering if one should / could not, take the product? E.g., a scale "friendship" is comprised of the variables 1) "fun to be with" and 2) "good to talk to when I am sad". (just a hypothetical example). Lets say that variable 1 is answered on a 1-5 scale, and and variable 2 is answered on a 1-10 scale. Lets say that the two variables correlate nicely with a Pearson coefficient of 0.6. Now, with both variables scaled from eg. 0 to 1, what would be the best friendship measure? The average of the two variables or the product? If I use the product measure, friendship has a much larger variance than if I use the standard averages measure. Is there a simple mathematical reason why I should not take the product? Or does it depend on theoretical assumptions about "friendship"?


title: How can I evaluate the performance of a system that generates word clusters?

text: The word2vec tool uses deep learning to compute vector representations of words. They've mentioned that - "The word vectors can be also used for deriving word classes from huge data sets. This is achieved by performing K-means clustering on top of the word vectors." When I execute the code, it seems that the tool computes word vectors from a dataset and clusters it into 200 classes. What if I wanted to use another method of clustering? How would I measure how 'good' these classes are, or how well the clustering algorithm is working? Also, it seems to me that the word vector representation would be a very sparse vector space. Couldn't the clustering be improved by PCA or some method of dimension reduction? Again, how can I test how well the clustering is performing? If I can figure out a way to evaluate, I can change the code and add PCA and see if it does better. EDIT - I found some papers that mention that word clustering can be used in document classification. The idea is to use the word clusters instead of a bag-of-words as features for classification, and this reduces the feature space while preserving 'redundant' features. I'm looking more into this, but it seems like I could use the word clusters on a document classification task and see if it performs better if I change the clustering algorithm. Would this work? A follow-up question - how would I implement this approach of using word clusters as features?


title: When does counterbalancing task participation order take care of practice effects?

text: I'm thinking of a study in which participants are randomly assigned to do task A then task B, or vice-versa. I figured there will be a problem if doing Task A is really good practice for Task B, but doing task B isn't good practice for Task A. However, let's say I have good grounds for thinking Task A's practice effect on Task B will be the same as Task B's practice effect on Task A. If I can assume that much am I clear to proceed? Or is there something else I should be worried about?


title: Evaluating Time Series Prediction Performance

text: I have a Dynamic Naive Bayes Model trained on a couple of temporal variables. The output of the model is the prediction of `P(Event) @ t+1`, estimated at each `t`. The plot of `P(Event)` versus `time` is as given in the figure below. In this figure, the _black line_ represents `P(Event)` as predicted by my model; the _horizontal red line_ represents the prior probability of the event happening; and the _dotted vertical lines_ represent the (five) event occurrences on the time series. Ideally, I wish to see the predicted `P(Event)` peak prior to observing any events and remain close to zero when there is no prospect of an event. ![P\(Event\) versus Time Graph](http://i.stack.imgur.com/ldhhT.png) I want to be able to report how well my model (the black line) performs in predicting the event occurrences. An obvious candidate to compare my model with is the prior probability of event (the red line), which -if used as a predictor- would predict the same probability value for all `t`. What is **the best formal method** to achieve this comparison? **P.S:** I am currently using the (intuitive) scoring as coded below, where an overall lower score indicates better prediction performance. I found that it is actually quite difficult to beat the prior with this scoring:               # Get prediction performance     model_score = 0; prior_score=0;           for t in range(len(timeSeries)):             if(timeSeries[t]== event):  # event has happened           cur_model_score = 1- prob_prediction[t];            cur_prior_score = 1 - prior        else: # no event           cur_model_score = prob_prediction[t] - 0;           cur_prior_score = prior - 0;             model_score = model_score + abs(cur_model_score);        prior_score = prior_score + abs(cur_prior_score);


title: Confusion on a Statistical Statement

text: I am confused about the statistical statement I am making based on the following scenario: **Scenario:** The membership of a professional organisation is only open to the top 10% of students in a district. There are 400 members, coming from (say) 45 schools. **Statement** : The organisation represents the top 10% of students in the district. My confusion is that School A may have lower quality students who make the top 10% than School B. e.g.   * School A has students who scored 41, 42, 43, 44 and 45 out of 50 (top 10%)   * School B has students who scored 46, 47, 48, 49 and 50 out of 50 (top 10%) In this case, can they be the top 10% of students (especially when School B has students who scored 41, 42, 43, 44 and 45)? I am using just 2 schools for simplicity.


title: Why do we need deviance

text: Deviance > is a generalization of the idea of using the sum of squares of residuals in > ordinary least squares to cases where model-fitting is achieved by maximum > likelihood. As I understand it, deviance is a statistic by which we could judge if a model is good or better. However, IMHO, we've already have RSS (for regression) and classification error rate (for classification), and by these 2 statistics, I think we've been able to judge a model quality. So why do I have to computed deviance?


title: Normalized root mean squared error (NRMSE) vs root mean squared error (RMSE)

text: The response values in my data set (100 data points) are all positive integers (should not be either negative or zero values). I have developed two statistical models: Linear Regression (LR) and K Nearest Neighbor (KNN, 2 neighbours) using the data set in R. The R methods I have used are lm() and knn.reg(). To select between these two models, I have conducted 10 fold cross- validation test and first computed root mean squared error (RMSE). Although the LR model is giving negative prediction values for several test data points, its RMSE is low compared to KNN. When I see the prediction values of KNN, they are positive and for me it makes sense to use KNN over LR although its RMSE is higher. Moreover, when I used Normalized RMSE (http://en.wikipedia.org/wiki/Root-mean-square_deviation), KNN has low NRMSE compared to LR. Furthermore, I would like to define "prediction accuracy" of the models as (100 - NRMSE) as it looks like we can consider NRMSE as percentage error. Please let me know the above methodology I am following is fine or not. Thank you.


title: detect line in scatter

text: I have repeated samples from the following process: Most samples will only contain points that are randomly distributed on a 2-dimensional plane. Sometimes, however, the sample will contain not only randomly distributed points but also some points that will be arranged roughly on a line. I say roughly because the points I observe may not all be on the line exactly but be randomly displaced somewhat. Here's the R code for a little toy example to demonstrate what I mean:               line_active <- TRUE          line_n <- 50     noise_n <- 100     perturbation.sd <- 0.01          # background noise     noise_x <- runif(noise_n)     noise_y <- runif(noise_n)          # line coordinates     start <- list(x=0.2,y=0.1)     end <- list(x=0.8,y=0.9)          if(line_active) {       line_pos <- runif(line_n)       line_x <- start$x+line_pos*(end$x - start$x)           line_y <- start$y+line_pos*(end$y - start$y)            line_x_perturbed <- line_x+rnorm(length(line_x),0,perturbation.sd)       line_y_perturbed <- line_y+rnorm(length(line_y),0,perturbation.sd)     } else {       line_x_perturbed <- c()       line_y_perturbed <- c()     }          points <- data.frame(x=c(noise_x,line_x_perturbed),y=c(noise_y,line_y_perturbed))          # plot!     plot(points$x,points$y)     points(line_x_perturbed,line_y_perturbed,col='red')      My questions are:   1. How would I go about detecting whether or not there is a line in any sample I'm looking at?   2. How would I go about detecting the slope and intercept of the line if it is there?   3. How much more tricky does this problem get if the noise does not come from some known distribution but has to be estimated?


title: Relation between omnibus test and multiple comparison?

text: Wikipedia says > Methods which rely on **an omnibus test** before proceeding to **multiple > comparisons**. Typically these methods require a significant ANOVA/Tukey's > range test before proceeding to multiple comparisons. These methods have > "weak" control of Type I error. Also > The F-test in ANOVA is an example of an omnibus test, which tests the > overall significance of the model. Significant F test means that among the > tested means, at least two of the means are significantly different, but > this result doesn't specify exactly what means are different one from the > other. Actually, testing means' differences has been made by the quadratic > rational F statistic ( F=MSB/MSW). In order to determine which mean differ > from another mean or which contrast of means are significantly different, > Post Hoc tests (Multiple Comparison tests) or planned tests should be > conducted after obtaining a significant omnibus F test. It may be consider > using the simple Bonferroni correction or other suitable correction. So an omnibus test is used to test the overall significance, while multiple comparison is to find which differences are significant. But if I understand correctly, the main purpose of multiple comparison is to test the overall significance, and it can also find which differences are significant. In other words, multiple comparison can do what an omnibus can do. Then why do we need an omnibus test? Thanks and regards!


title: need derivation of join-count variance (spatial autocorrelation stat), know where it is?

text: I am using interlibrary loan to get Cliff and Ord's book Spatial Processes, but the semester just ended and it is slow now. On page 18ish of this book, Cliff and Ord show how the variance for the join-count statistic is calculated and also how this variance was derived. I have the variance equation, but need to see the derivation because I need to calculate variance of join-counts WHILE taking edge effects into account. So, I must modify the variance but can only do this if I understand the steps in its derivation very well (which is hard to do when looking at a reduced form equation). Time is of the essence, has anyone seen an online reference that shows this derivation and can point me to it? It might take quite a while to have this book in hand, given the situation. Also, buying a book like this for a few pages is a bit extravagant given the cost of books like this. Thanks. I doubt anyone can help, it is too odd a question likely.