title: Seasonal vs non-seasonal coefficients in R ARIMA

text: Let's say I have the two following ARIMA models:   1. ARIMA(7,1,1) (no seasonality)    2. ARIMA(6,1,1)(1,0,0)7 (seasonality of period 7). Are they conceptually the same? If so, why is that when I model them using R, I get different models and different fits?


title: How to test the effect of condition on change scores for a possibly ordinal dependent variable?

text: I have a categorical independent variable (two levels: condition 1; condition 2) and ordinal(?) dependent variables (numerical magnitude shifts (e.g. -3 if there was a decrease in number magnitude from 4 to 1). I now want to assess whether the condition predicts a change in the dependent variable. It is a repeated measurements design. Now I have got the following questions:   1. Can I conduct a linear regression for DV = ordinal (categorical)?   2. Can I conduct a linear regression for repeated measurements design?


title: In practice, how is the penalisation matrix for splines, created by smoothCon (package mgcv), specified?

text: Evaluating a smooth object with `smoothCon` provides, besides several other things, the ordinary "untouched" spline bases               objx <- te(x,bs="ps",k=4)     Con.mgcv.x <- smoothCon(objx,dat,knots=NULL,scale.penalty=FALSE)     Mx.mgcv <- Con.mgcv.x[[1]]$margin[[1]]     Mx.mgcv$X       the reparameterized version $\boldsymbol{X}\prime$               Con.mgcv.x[[1]]$X      which corresponds to the smooth parameterized in terms of function values at a set of knots: $f_x(x) = \boldsymbol{X}\boldsymbol{\alpha}$, set $\boldsymbol{\alpha}\prime=\boldsymbol{\Gamma}\boldsymbol{\alpha}$ where $\Gamma_{ij} = \text{b}_i\left(x=x_j\right)$, $\boldsymbol{X}\prime = \boldsymbol{X}\boldsymbol{\Gamma}^{-1}$ the $\boldsymbol{\Gamma}^{-1}$ Matrix               Con.mgcv.x[[1]]$XP[[1]]      Now in textbooks such as Wood, Generalized Additive Models: An Introduction with R the penalisation matrix $\boldsymbol{S}$ should also be reparameterized as $\boldsymbol{S}\prime = \boldsymbol{\Gamma}^{-T}\boldsymbol{S}\boldsymbol{\Gamma}^{-1}$. If simply checking the output of `smoothCon` with what one would get if building the $\boldsymbol{\Gamma}$ matrix themself and calculating $\boldsymbol{S}\prime$ the results differ though $\boldsymbol{X}\prime$ is essentially equal. What is the reason for that? Is some further rescaling done since my $\boldsymbol{S}\prime$ matrix is rougly $10000$-times larger then what is given by `smoothCon`.


title: Interpret t-values when not assuming normal distribution of the error term

text: Assume that you have a regression with a whole set of variables and you know that the residuals are not normal distributed. So you just estimate a regression using OLS to find the best linear fit. For this you disclaim the assumption of normal distributed error terms. After the estimation you have 2 "significant" coefficients. But how can anyone interpret these coefficients? So there is no way to say: "These coefficients are significant", although the Hypothesis $\beta=0$ can be declined with a high t-statistic (because of disclaiming normal error assumption). But what to do in this case? How would you argue?


title: Why does noisy data result in better prediction performance?

text: I have tested a regression framework's robustness to noise and I have noticed in some cases that adding noise improves the prediction performance and in other cases the performance degrades. What could be the reasons for this? If there are multiple reasons, how to I determine which is the cause? Edit: Some more details about what I am doing. The framework uses ridge regression. The inputs are vectors of extracted image features. The outputs are vectors of angles (in degrees, -180 to 180). To test for robustness to noise I am applying 3 levels of noise (white additive Gaussian noise) to the angles (targets) proportional to the individual angle variances (2%, 5%, and 10% of the variance of each angle). I have noticed that in some observations, adding a small amount of noise (2-5%) leads to a small improvement in performance and in one case, all levels of noise give improvement. In my tests, the regularisation term is fixed across all noise levels, and I have ran each noise level test several times to take into account the fluctuations of the random noise. Also, I have two broad sets of observation data. The first set was observed relatively accurately, however the second set was more complex (significantly more heterogeneous, leading to notable performance degradation relative to the first set) and exhibited a number of minor errors due to the observation technique being more limited than that which was used with the first set. In the first set, the phenomena of better performance through adding noise did not occur. However, sometimes more noise was better than less noise. If more information is required to better answer the question, I'd be happy to provide it.


title: Multivariate multi-level analysis in nlme

text: **The question** I have a dataset which I think requires a multivariate multilevel analysis. I am unsure both of the appropriate model and of how to fit it with `R`. I have come up with a tentative model, but my understanding of the math is so superficial that I cannot tell whether my analysis is "right" or whether it includes blatant errors. I would appreciate any insight on the model design or the model specification in `R`. **The study design** The question is whether the architectural design of a clinic will influence the outcome of a pathology for permanent residents in this clinic. We have collected data on 13 symptoms for 8 patients per clinic in 21 clinics. There is a patient-level IV `medication` and a clinic-level IV `architecture`. All variables are continuous-ish. The 13 symptoms are correlated +.20 on average, which I think indicates a multivariate multilevel analysis is appropriate. **The data** To run the multivariate analysis with `nlme` I have standardized my DVs, stacked these 13 DVs in a single column, and added a categorical dummy variable to flag which row corresponds to which symptom. It looks like this:                Clinic Patient Symptom    Score    Medication  Architecture      1            1   EP1      0.12         1               3.2        1            1   EP2      0.11         1               3.2        1            1   EP3      0.13         1               3.2        1            2   EP1      0.56         4               3.2        1            2   EP2      0.67         4               3.2        1            2   EP3      0.23         4               3.2        2            3   EP1      0.22         3               5.1        2            3   EP2      0.25         3               5.1        2            3   EP3      0.14         3               5.1        2            4   EP1      0.78         6               5.1        2            4   EP2      0.89         6               5.1        2            4   EP3      0.11         6               5.1        **The model design**   * To run the analysis as multivariate, I use both `symptom` and `symptom:architecture` as IVs and I remove the intercept in both the fixed and random parts of the model. I do not include the main effect of `architecture` as an IV.    * The effect of `medication` should be the same within all clinics, so there is no random effect for this variable.    * I do not want to constrain equality between the effect of `architecture` on the different `symptoms`.    * Due to the multivariate nature of the analysis, I expect the residuals to be correlated, with different correlations between the 13 different symptoms; therefore I specify the covariance structure of residuals as `corSymm` (non-zero but unstructured, if I get this correctly).    * I also expect heteroscedasticity between the different symptoms (there should be more variance on certain symptoms), so I add the option `weights` as `(~ 1|symptoms)`.  **The end result** This is the model I come up with:               model1 = lme(fixed = Score ~ symptom + medication:symptom + architecture:symptom + medication:architecture:symptom - 1,       + random = ~ symptom - 1 | patient/clinic,       + correlation = corSymm,       + weights=varIdent(form= ~ 1|symptoms)       + method = "ML")        In order to test the effect of the architectural variables, I would then compare this model to the following constrained model, dropping all the terms related to architecture:               model2 = lme(fixed = Score ~ symptom + medication:symptom - 1,     + random = ~ symptom - 1 | patient/clinic,     + correlation = corSymm,     + weights=varIdent(form= ~ 1|symptoms)     + method = "ML")      I would then run this comparison with the command `anova(model1, model2)` and compare the log-likelihood of the two models. Overall, do these model design and r specification look correct to you? Thank you so much for your help!


title: Relation between categories

text: I'm trying to analyze an experiment, in which every subject is asked the same question in different conditions. I then describe each of their answers into categories describing their answers in detail. I'm trying to find relations between different categories. For example, several users whose answers can be described by category A can also often be described with category B but is not related at all with category C. I thought that Factor Analysis could find the relations between categories, but it seems that for non-numeric data, Correspondence Analysis is more appropriate, is this correct? Thank you.


title: Evaluate adequacy of hidden Markov models

text: I am trying to fit an HMM to a dataset resembling:               id0: A, A, A, C, B, A, C, B     id1: C, B, A, A, C, C     id2: B, A, A, C, B, B, A      In total there are three characters (A, B, C) and about 40k of such rows. I'm using R's `HMMFit` from the `RHmm` package. As input to `HMMFit` one has to select the desired number of hidden states, so I can generate different models for different numbers of hidden states. In addition to the model, `HMMFit` reports the log-likelihood, the BIC criterium and the AIC criterium. How does one evaluate the adequacy of such a model? Should one infer the model from training data as opposed to the entire dataset like I do now? If yes, how does one check if the trained model describes the test data? Given different potentially adequate models generated for different numbers of hidden states, which one should I prefer?


title: Analyzing residual plot vs independent variables plot

text: Why do we analyze residual plot in regression analysis and NOT between two individual variables? For example when checking for normality, heteroscedasticity etc. we don't analyze two individual variables but residual plot, why so?


title: Less significant findings using conditional logistic regression

text: I used regular logistic regression on my dataset and got a few significant hits. However, since the data is 1:1 case-control matched data I decided to try using conditional logistic regression (`clogit()` in the `survival` package of `R`). However, none of the results are significant any more. Is this normal? I was hoping to gain power by switching to conditional logistic regression.                m<-clogit(PHENO==2 ~ x + as.factor(COVAR[,1]) + strata(COVAR[,2]) )      This is a GWAS (genome wide association study). `PHENO` is either 1 (control) or 2 (case). `x` is number of copies of the minor allele (0, 1, or 2) in the SNP of interest, `COVAR[,1]` is a race covariate and `COVAR[,2]` is the matching covariate. This test is performed for each SNP (~300,000) and the resulting p-values are adjusted for FDR. Perhaps the matching of the cases and controls was not very good? Matching: case and control pairs had: 1) same sex. 2) BMI > threshold. 3) same self-reported ethnicity. 4) lived within 50 miles of procurement site. Perhaps I should use GLMM?


title: Overfitting in GLMM

text: I have fit a mixed model with `lmer()` and am left with 4 significant interaction terms. There were found by removing the interaction term and comparing with the full model using `anova(fm1, fm2)`. I have also left the single terms in the model. When I come to report the chi-square for the single terms by removing them and using `anova(fm1, fm2)` I am told there is 0 df and chi-square is 0. I don't see how this is possible: is my model over-fit?


title: (Why) should bootstrap sampling distribution for logistic regression slope be conditional on $S=\sum Y_j$?

text: I am working my way through Chapter 4 (Tests) in Davison & Hinkley's (D&H's) book "Bootstrap Methods and their Applications", and have a question about one of their examples. ## Davison & Hinkley's proposed test Example 4.1 proposes a permutation test for the (alternative) hypothesis of non-zero slope in a simple logistic regression. Quoting from page 141: > **Example 4.1 (Logistic regression)** Suppose that $y_1, \ldots , y_n$ are > independent binary outcomes, with corresponding scalar covariate values > $x_1, \ldots , x_n$ and that we wish to test whether or not $x$ influences > $y$. If our chosen model is the logistic regression model $$ \log > \frac{\Pr(Y_j = 1 | x_j)}{\Pr(Y_j = 0 | x_j} = \lambda + \psi x_j, \quad > j=1,\ldots,n, $$ then the null hypothesis is $H_0 : \psi=0.$ Under $H_0$ the > sufficient statistic for $\lambda$ is $S=\sum Y_j$ and $T=\sum x_j Y_j$ is > the natural test statistic > > [...snipped text...] > > The null distribution of $Y_1,\ldots,Y_n$ given $S = s$ is uniform over $n > \choose s$ permutations of $y_1, \ldots, y_n$. Rather than compute all of those permutations to get the exact distribution of $T$, they propose leaving the $x_i$ fixed and generating $R$ permutations of the $y_i$ to approximate $T$'s exact distribution. ## My naïve test That all seems reasonable, but my thought would have been to simulate a slightly different null distribution. Given that $\psi = 0$ under the null hypothesis, $$ \Pr(Y_j = 1| x_j) = \frac{e^{\lambda}}{e^{\lambda} + 1}, \quad j=1,\ldots,n, $$ Leaving $x_1, \ldots , x_n$ fixed, I would produce $R$ sets of $y^*_1, \ldots , y^*_n$ with each $y^*_i \sim \textrm{Bernoulli}(\frac{e^{\lambda}}{e^{\lambda} + 1})$. Instead of (but probably equivalently to) D&H's $T=\sum x_j Y_j$, I would take as my test statistic the estimate of $\psi$ from a fitted logistic regression of the bootstrap sample. ## Questions The two tests above depend on different distributions of bootstrap samples. Mine, for instance, may produce some samples in which all of the $y^*_i$ equal 1, whereas all of D&H's samples contain the same number of 1's. The salient difference seems to be that D&H's sampling distribution is conditional on $S=\sum Y_j$ whereas mine is not. Is one of these approaches preferable to the other in all situations? They seem to assume/model different data-generation processes in the originally sampled populations. Is that correct? If so, when should I choose one approach over the other?


title: Standard Deviation of percentages

text: Say I have a bunch of weighted coins that should show heads 90% of the time. I take a bunch of them, throw them in a box, and calculate the percentage of coins that landed on heads. I do this a bunch of times and calculate a standard deviation of 5% and a mean of 90%. Assuming a normal distribution, I then (ignorantly) decide in subsequent runs I can plan on:   * 34.1% (1 STD DEV) of the time I will see 90% to 95% heads   * 13.6% (2 STD DEV) of the time I will see 95% to 100% heads   * **2.1% (3 STD DEV) of the time I will see 100% to 105% heads** This last entry, of course, is impossible. I am assuming the entire analysis is fundamentally wrong. What is the correct way to interperet the standard deviation when dealing with situations like the above? (Finally, I'm a statistics noob so feel free to set me straight if I'm way out in left field!)


title: Multicollinearity between categorical and continuous variable

text: Simplified version of my model:               glm(cbind(young, adults) ~ as.factor(month) + effort, family = "binomial")      i.e., I study proportion of young as a dependent variable on month (or season), taking into account observer effort. However, the observer effort is dependent on the month: ![enter image description here](http://i.stack.imgur.com/9VZ2I.png) How to solve this problem? I want to take into account both variables. I was looking into literature but haven't found any solution. My naive idea is to compute mean effort for each month and instead of taking effort, take difference between effort and this mean. But this is just a naive idea, I would like to hear your advice. Thanks! **EDIT** - response to Scortchi question - no, not actually:               > m = glm(cbind(young, adults) ~ as.factor(month) + effort, family = "quasibinomial")     > summary(m)          Call:     glm(formula = cbind(young, adults) ~ as.factor(month) + effort, family = "quasibinomial")          Deviance Residuals:          Min       1Q   Median       3Q      Max       -2.6829  -1.1138   0.0000   0.9717   4.0090            Coefficients:                          Estimate Std. Error t value Pr(>|t|)         (Intercept)        -0.7868764  0.2170738  -3.625 0.000306 ***     as.factor(month)2   0.8857561  0.2606780   3.398 0.000710 ***     as.factor(month)3   0.7055741  0.2918895   2.417 0.015843 *       as.factor(month)4   0.3943665  0.3269973   1.206 0.228138         as.factor(month)5   0.4831113  0.3730987   1.295 0.195713         as.factor(month)6  -0.5217349  0.5027560  -1.038 0.299676         as.factor(month)7   0.1612901  0.4333682   0.372 0.709851         as.factor(month)8   0.5114890  0.3545159   1.443 0.149444         as.factor(month)9   0.7741060  0.3126087   2.476 0.013466 *       as.factor(month)10  0.6601093  0.2609937   2.529 0.011608 *       as.factor(month)11  0.4891778  0.2647303   1.848 0.064967 .       as.factor(month)12  0.4743091  0.2565709   1.849 0.064849 .       effort              0.0032506  0.0007976   4.075 5.02e-05 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for quasibinomial family taken to be 1.270962)              Null deviance: 1518.0  on 878  degrees of freedom     Residual deviance: 1447.8  on 866  degrees of freedom       (750 observations deleted due to missingness)     AIC: NA          Number of Fisher Scoring iterations: 4          Warning message:     In summary.glm(m) :       observations with zero weight not used for calculating dispersion     >      >      > m = glm(cbind(young, adults) ~ as.factor(month), family = "quasibinomial")     > summary(m)          Call:     glm(formula = cbind(young, adults) ~ as.factor(month), family = "quasibinomial")          Deviance Residuals:          Min       1Q   Median       3Q      Max       -3.1142  -1.1266   0.0000   0.9235   3.6484            Coefficients:                        Estimate Std. Error t value Pr(>|t|)         (Intercept)         -0.2296     0.1590  -1.444  0.14890         as.factor(month)2    0.8610     0.2059   4.181 3.09e-05 ***     as.factor(month)3    1.0887     0.2062   5.280 1.51e-07 ***     as.factor(month)4    0.4184     0.2374   1.762  0.07822 .       as.factor(month)5    0.1495     0.3086   0.485  0.62802         as.factor(month)6   -0.6177     0.3872  -1.595  0.11091         as.factor(month)7   -0.4636     0.3666  -1.265  0.20622         as.factor(month)8    0.1089     0.2976   0.366  0.71440         as.factor(month)9    0.4932     0.2490   1.980  0.04787 *       as.factor(month)10   0.6322     0.2096   3.016  0.00261 **      as.factor(month)11   0.4919     0.2152   2.286  0.02243 *       as.factor(month)12   0.2296     0.2127   1.079  0.28071         ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for quasibinomial family taken to be 1.309215)              Null deviance: 2400.1  on 1345  degrees of freedom     Residual deviance: 2310.4  on 1334  degrees of freedom     AIC: NA          Number of Fisher Scoring iterations: 4          Warning message:     In summary.glm(m) :       observations with zero weight not used for calculating dispersion     >      >      >      > m = glm(cbind(young, adults) ~ effort, family = "quasibinomial")     > summary(m)          Call:     glm(formula = cbind(young, adults) ~ effort, family = "quasibinomial")          Deviance Residuals:         Min      1Q  Median      3Q     Max       -2.718  -1.119   0.000   1.011   4.236            Coefficients:                  Estimate Std. Error t value Pr(>|t|)         (Intercept) -0.326899   0.102473  -3.190  0.00147 **      effort       0.003827   0.000688   5.563 3.52e-08 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1           (Dispersion parameter for quasibinomial family taken to be 1.268467)              Null deviance: 1518.0  on 878  degrees of freedom     Residual deviance: 1475.4  on 877  degrees of freedom       (750 observations deleted due to missingness)     AIC: NA          Number of Fisher Scoring iterations: 4          Warning message:     In summary.glm(m) :       observations with zero weight not used for calculating dispersion


title: Explanation for the thresholds in the sequential probability ratio test

text: Recently I have learned about sequential analysis especially sequential probability ratio tests (after I have struggled a lot with cumulation of alpha - errrors). See also this question Sequential hypothesis testing in basic science. My question is: What is the derviation or explanation of the tresholds of the stopping-rule (Again see Sequential probability ratio tests - Theory)? How do these thresholds prevent the errors from accumulation? I am talking about the scheme: $a < S_i < b$ where $S_i$ is the likelihood-ratio and $a:=\frac{\beta}{1-\alpha}$ and $b:=\frac{1-\beta}{\alpha}$. Why are a and b set this way ? I'd prefer an intuitive explanation, but a mathematical one using not too many rarely known concepts is fine, too.


title: What's the difference between a probability and a proportion?

text: Say I have eaten hamburgers every Tuesday for years. You could say that I eat hamburgers 14% of the time, or that the probability of me eating a hamburger in a given week is 14%. What are the main differences between probabilities and proportions? Is a probability an expected proportion? Are probabilities uncertain and proportions are guaranteed?


title: Shannon entropy for non-stationary and non-linear signal

text: I would like to quantify the amount of uncertainty in a given message, but the signal I work with is non-stationary and non-linear. Is it possible to apply Shannon entropy for such signal?


title: How to know if a list of prices are mean-reverting?

text: I have two historical price lists with the following columns > **data - price** Now i have to create the ratio between the prices of these lists: > **list A:** 01/01/2011 10.50 > > **list B:** 01/01/2011 23.89 I compare the date of the lists, if the day is the same I find the ratio, doing: > **ratio** = 10.50 / 23.89 Ok... I do this division for each price of the lists. The result should be:               0.4395   0.4400   0.4289   0.4361      Now the question is: **How could I check if this serie (ratios) is mean reverting?**


title: Is there a boxplot variant for Poisson distributed data?

text: I'd like to know if there is a boxplot variant adapted to Poisson distributed data (or possibly other distributions)? With a Gaussian distribution, whiskers placed at L = Q1 - 1.5 IQR and U = Q3 + 1.5 IQR, the boxplot has the property that there will be roughly as many low outliers (points below L) as there are high outliers (points above U). If the data is Poisson distributed however, this does not hold anymore because of the positive skewness we get _Pr(X <L) < Pr(X>U)_. Is there an alternate way to place the whiskers such that it would 'fit' a Poisson distribution? Thanks for your help !


title: Are Archimedean copulas useless for representing multivariate data?

text: Following Hofert et al.'s paper "Likelihood inference for Archimedean copulas in high dimensions under known margins," (http://dl.acm.org/citation.cfm?id=2263953) I wrote a script in Matlab to produce estimates of Archimidean copulas in high dimensions. The procedure in their paper allows computing the derivatives of the generator functions without resorting to iterative methods; just closed forms. So I thought it a would idea to employ it as an alternative to the use of elliptical copulas for some multivariate modeling. I used it to represent data with extreme dependence using Gumbel copulas. However, when I employed it to represent more than say four variates, the estimates weren't as good as with elliptical copulas! In a paper that I subsequently wrote a reviewer made me notice that this was because _Gumbel copula has only one parameter_ , which turns it out to be useless to represent high dimensional data. I should employ, he/she said a vine copula in those cases (with added complexity); not a one-parmeter copula. Thereafter my question is, Are Archimidean copulas really useful to represent high dimensional data? Are they useless except for maybe homogeneous dependence structures, where the dependence might be equal between pairs of variates? Any example of a _good_ representation of multivariate data?


title: How is EM used in the sense of data mining on images?

text: I understand how EM is used in the sense of estimating the Gaussian model that underlies a set of data, but its unclear how this is applicable. I am trying to understand how EM might be used to perform data mining on any sort of computer vision/image processing task. This is the domain I am most familiar with, so it is helpful when learning about new algorithms to frame them in a familiar framework. So one again, how and where is expectation maximization used in data mining images? Thank you


title: Maximum entropy classifier and sentiment analysis

text: I am doing a project work in sentiment analysis (on Twitter data) using machine learning approach. In order to find the 'best' way to this I have experimented with naive Bayesian and maximum entropy classifier by using unigrams, bigrams and unigram and bigrams together. I'm using the SharpEntropy library for ME, and an own implementation for the NB. My question is related to the ME classifier. In one of my tests I've trained it using both unigrams and bigrams because, these outperformed the unigram only and bigram only solutions. I used 100 iterations and no cutoff (= minimum appearance for using it in the training). After the initial training I've started the validation and the strange thing I encountered is the following: if the test data is not converted to the unigram-bigram representation my classifier was 4% more accurate (73%->77%) than the case when I did the transformation. So basically my classifier is trained using uni- and bigrams, but the test data remained in the unigram representation and I got better results. (unigram training - unigram test data resulted in 58% accuracy, bigram only training - bigram only test resulted in 68% accuray) What could be the reason for this behavior? I am looking for some help to get a deeper understanding for this.


title: theoretical profile correlation

text: My data (n > 4000) contains the results of two separate tests (instruments). I have a profile made from likert scales for one of the instruments (S1, S2, ... S20). I want to find the corresponding profile for the 2nd instrument that would correlate perfectly (or near perfectly) with the first instrument. It seems there should be a permutation that would find this unknown profile.


title: Clustering mixed variables in SAS

text: I have following variables in my dataset:   1. Working hours (numerical:ordinal)   2. Effectiveness (categorical:ordinal ; 4 values-> (poor,average,good,best))   3. Satisfaction (categorical:ordinal ; 4 values-> (poor,average,good,best)) I want to cluster the data on the basis of how good is my worker. I am expecting 4-5 clusters effectively. I ran fastclus in sas after normalising my data (mean=0 std=1) But i doubt the k means will work on categorical variables. I have getting an Rsqaure value of 62% when i have 5 clusters. Any help on how to go about this problem?


title: Is segmented (piecewise) regression possible with a dichotomous dependent variable?

text: I have created a model of American cities with a dichotomous Y variable using logistic regression. I have theoretical reasons to believe that the model will differ significantly between larger and smaller cities. I have shown this using arbitrary breakpoints, but wanted to know if a segmented (piecewise) regression strategy can identify a breakpoint for me in a model with a binary dependent variable. (I have used the "nl" command successfully in Stata for a different project when the Y was continuous.) Can I use this strategy to statistically identify a breakpoint? Or should I build the case for a number I pick otherwise? Many thanks.


title: Calculation of standard variance $s^2$

text: First picture is the question and its answer key: ![enter image description here](http://i.stack.imgur.com/3qWCe.jpg) The second picture is my solution: ![enter image description here](http://i.stack.imgur.com/cIlt5.jpg) I know that $S^2=\frac{\sum(x_i-\bar x)^2}{n-1}$ But in the answer key, $s^2$ cannot be found in that way. I guess, there is a mistake here. Or, ıs there a rule that I dont know? Please clarify this. Thank you.


title: Reading a crosstabs

text: I am confused and struggling to read this simple crosstabs. ![Crosstabs](http://i.imgur.com/NNdbMsP.jpg) I am trying to create a profile of downloaders. One of the demographics I am looking at is age. Should I be looking at the percentage within Age or percentage within Type of downloader? For example, Mixed downloaders are most likely to be 18-21 (50%)? Or should it be they are most likely to be 18-21 (68.4%)?? Could someone help explain please Thanks


title: What is the best single measure of brain networks based on functional MRI data?

text: I want to analyze functional MRI (fMRI) data as follows:   1. I am comparing the brain networks of two groups of subjects (patients and matched).   2. For each subject I have a correlation (Pearson) matrix with the correlations between multiple brain regions (above 50000 correlations).   3. I want to get a single measure to represent every subject (maybe the mean of correlations of each subject network). **What is the most suitable measure?** Should I use Fisher transformation to average the correlation matrix?    4. I assume that the variance of the two groups is different and that the variance in each group is similar: **How can check that without doing multiple comparisons?**


title: Quadratic weighted kappa versus linear weighted kappa

text: When should I use quadratic weighted kappa or linear weighted kappa? I have two observers evaluating the classes of a number of objects. The classes are fail, pass1, pass2, and excellent (ordinal scale). The errors in classification between "fail" or "excellent" and the different degrees of "pass" are more severe than errors between the classes of pass (pass1 and pass2). Could I define the values of the classes as "fail = 1", "pass1= 20", "pass2=25", "excellent=40" and use linear weighted kappa (higher penalties on the extremes in a intuitive/subjective way)? Or should I use quadratic weighted kappa? Why should I prefer one kind over the other?


title: Piecewise orthogonal regression

text: I'm trying to calculate a single breakpoint on $x$ from several thousand $(x,y)$ points. Due to the density distribution as per the 2D histogram below, the least squares method is proving unsuitable. An orthogonal regression is a much better fit for the lower $x$-axis, but is it possible to conduct a piecewise orthogonal regression? How could I do this, preferably in R? ![enter image description here](http://i.stack.imgur.com/xb1fd.png)


title: Clustering elements by access counts in sessions

text: I have a set of sessions and urls that have been accessed in each of these sessions and frequencies with which they have been accessed. I've put them in a matrix-like representation. Imagine I have the following "Pageview matrix":               COLUMN HEADINGS          books placement resources br aca      Each row represents a session. Here is an example of the records:               4 5 0 2 2     1 2 1 7 3     1 3 6 1 6      saved in a `.txt` file Can I give this as an input to a k-means program and obtain clusters based on the highest frequency of occurrence? How do I use it? If not k-means, what other cluster method can I use?


title: Counterexample for interaction and parallel curves?

text: It is said that if the plots of the hypothetical responses are not parallel, but crossed, there is interaction. Suppose we have two factors. Is it possible that the plots cross but we do not have interaction? That is more reasonable when the plots are close to each other.   I noticed the converse is true. We may have interaction even if neither the curves of factor A on factor B nor factor B on factor A intersect. In a book I read this happens when the interaction is removable, which means that there is another cognate independent variable.


title: What is the correct way of estimating the proportion of individuals in a population from a count of their individual parts?

text: My original question was going to be : How do you estimate the proportion of species in a population from count data of their individual body parts, provided that you can identify each part as belonging to a specific species. **Is it just enough to take the average number of body parts found per species?** (keeping in mind some species might have a higher frequency of a certain body part than others). This is related to: Estimating Diversity of operon types using HMMs across metagenomes: Mann-Whitney? Kruskal-Wallis? or other? Edit: More Explanation **EXAMPLE 1** - you take samples of animals that are coming out in pieces. A normal animal called Species 1 has 1 leg, 1 nose, and 1 eye, Each of these traits you can see and count. Species 2 has 2 traits you can identify: ears (of which it also has only 1), and wings (Also only 1). (assume this is true for K Species, with K_y traits) You take 1 sample from x locations each time you count traits. the following is an example from 1 location: 1 leg, 5 noses, 12 eyes, 6 ears, and 8 wings... Which means the following Number of Species 1's possible are : Minimum 12 individuals ; Maximum 18 individuals (1+5+12): number of Species 2's possible are : Minimum 8 individuals ; Maximum 14 individuals (8+6). So in a bag with 1 eye and 1ear I know I have 1 individual of each species(no variation), in a bag with 1 ear 1eye 2 wings and 5 legs I know I have a maximum of 6 individuals (5+1) of species 1 and a min of 5. Species 2, max of 4 min of 3. So here we have possible limits for Species 1 and Species 2 at sample location x1.   What is the correct way to estimate and compare the proportions of each animal. We identify each part as belonging to a specific species; the eyes, noses and legs (any eye nose or legs found) belong to species 1, the highest count of individual traits means I have a bigger possible number of individuals belonging to Species 1, the more traits mean I have a larger range to predict the number of individuals. **EXAMPLE 2** You could also think of them as real animals, say chickens, lambs, Piglets with parts in a single bag. You can easily tell the difference between the three animals limbs, you know that if you get 18 chicken wings and only 3 drumsticks, then the maximum number of chicken slaughtered (assuming the drumsticks, and wings all come from different animals) is 18+3. Where as the minimum possible number of chicken slaughtered is where the wings are paired 18/2 (9 individuals), and the drumsticks are paired 3/2 (so only a min of 2 whole individuals) per chicken, **and** the wings and the drumsticks come from the same chickens(and some limbs are lost obviously). So the min# of individuals is the minimum number of pairs [9 whole chicken], and the more likely possible #individuals is the sum of all pairs[9+2 chickens],(but not necessarily???) So as you can see possible limits of chicken found by calculating their individual parts is : min- 9 chickens, and maximum 21 chickens…… the true number is probably somewhere in the middle(THIS is what I would like to calculate with a high degree of certainty, It could be 11 but I'm not sure of this. I want to be sure of this). Now let's count the number of pork chops, and pig heads…. etc Except With the first example the question is simpler, **because each trait is exclusive and only occurs once for species. different species will also have different traits** (You cannot count pork-chops and lamb chops and consider them from the same species.) Should I : -- Calculate a midpoint/mean as (( Maximum - Minimum )/(number of traits)) + Minimum for each species. The calculate the Standard deviation (for each) and compare the Coefficient of variances across each sample x? -- Take the Midpoint/mean calculated above, rank them and do non parametric pairwise comparisons of my x's -- Take the Maximum once and minimum once, rank each proceed non parametrically once for the upper limit once for the lower limit. --- There is probably a way to more accurately find the "true" number (or something close to it) of Individuals per species... If you have any ideas, please let me know. (I realize 1 sample is not nearly enough, if you would like more information see the link, I can only get 1 sample per location, i'd like to treat each trait as a separate sampling unit, these traits aren't really from animals etc..)


title: Empirical Bayes and random effects?

text: I recently happend to read about empirical Bayes (some introduction by Cassella?) and it looked a lot like random effects model; in that both have estimates shrunken to global mean... But I havent read it throughly... Does any one has any insight about the similarity and differences between them?


title: How do you assess a three-way interaction in a generalized multi-level model?

text: I am unsure which statistical tests I should use in order to follow-up my three-way interaction. I used to split my file and run the analysis twice per group, but as I read in the paper by Hayes & Matthes, 2009, _Behavior Research Methods_, this is not a good approach. I have a centered continuous variable (questionnaire score), one dichotomous variable and one variable called facetype with three categories. I have 6 trials per facetype condition per participant, so each participant has multiple lines in my SPSS file. My dependent variable is binary. I use the generalized linear model to analyze this data. The fixed effects show a significant three-way interaction, and I would like to interpret the interaction by for example testing the significance of one slope. I have no idea how to do this with multilevel data.               GENLINMIXED       /DATA_STRUCTURE SUBJECTS=id*trial       /FIELDS TARGET=in_team_yes_no TRIALS=NONE OFFSET=NONE       /TARGET_OPTIONS REFERENCE=0 DISTRIBUTION=BINOMIAL LINK=LOGIT       /FIXED  EFFECTS=questionnairescore Treatment Facetype questionnairescore*Treatment questionnairescore*Facetype Treatment*Facetype **questionnairescore*Treatment*Face**type USE_INTERCEPT=TRUE       /RANDOM USE_INTERCEPT=TRUE SUBJECTS=id COVARIANCE_TYPE=VARIANCE_COMPONENTS       /RANDOM EFFECTS=trial USE_INTERCEPT=FALSE COVARIANCE_TYPE=VARIANCE_COMPONENTS       /BUILD_OPTIONS TARGET_CATEGORY_ORDER=ASCENDING INPUTS_CATEGORY_ORDER=ASCENDING MAX_ITERATIONS=100 CONFIDENCE_LEVEL=95 DF_METHOD=RESIDUAL COVB=MODEL.


title: Can sub-optimality of various hierarchical clustering methods be assessed or ranked?

text: Classic agglomerative hierarchical clustering methods are based on a greedy algorithm. This means that they (many of them) are prone to give sub-optimal solutions instead of the global optimum result, especially on later steps of agglomeration. To clarify: each of the agglomerative methods makes the best choice of which two clusters to merge on a given step $q$, the choice that minimizes (to value $\delta_q$) the coefficient of colligation $\delta$ on the step; however, it is not impossible that if the algo had chosen _not the best but a bit worse_ choice on _some previous_ step(s) it then would have been able to minimize the coefficient on step $q$ to value less than that $\delta_q$ while preserving monotonical growth of $\delta$. Globally optimal solution on step $q$ corresponds to the absolutely minimal such way achievable value $\delta_q$ with the restriction remaining that $\delta_q \ge \delta_{q-1} \ge \delta_{q-2} \cdots$. Risk of sub-optimality is the main reason why hierarchical clustering is commonly not recommended with large amount of objects (say, more than 300): the researcher typically wants few clusters and therefore looks at later steps, but if a later step is, say, a 1000th one, it is more suspicious for having gone astray from globally optimal partition possible on the 1000th step than, say, a 100th step gone from globally optimal partition possible on the 100th step. This looks a matter of dictum. My question is what do you think about **which among well-known agglomerative methods** - single linkage, complete linkage, average linkage (within groups), average linkage (between groups), centroid, median, Ward's sum-of-squares - I'm mentioning those in SPSS, but there are other similar variants as well - **are more and which are less prone to the sub-optimality flaw** as the step number grows. Intuition suggests to me that single or complete linkage aren't prone at all and always give their globally best solutions, because these methods don't engage in computation of any statistics derivative from raw distances (e.g. centroids). I may be right as well as not. What are about the rest methods? Can someone here attempt to analyze the (relative) degree of local-optimum risk of the above concrete greedy algorithms? **Illustration.** Fortunately, for _Ward_ method we _have_ clues to observe how sub-optimality accumulates as $q$ grows. Because there is a well-known iterative method that attempts to optimize exactly the same function $\delta$ as Ward does; this is _K-means_ clustering: both try to minimize pooled within-cluster sum-of-squares. We can do Ward clustering and on each step save cluster centres, and use those as initial centres in K-means procedure. Will K-means improve the Ward's solutions in terms of sum-of-squares? Note. This illustration tests Ward's observed $\delta_q$ values not against _its own_ global optimum values which we can't know because we are unable to try, for each step, all innumerable alternative reconstructions of earlier steps in order to find out the minimal possible $\delta_q$; it tests against _K-means_ optimality (which is almost its global here because initial centres provided by Ward can be regarded reasonably good start for iterations). Here is some data (5 not well separated clusters, 183 points). The data underwent both Ward clustering session and K-means clustering sessions, as described. ![enter image description here](http://i.stack.imgur.com/T2emJ.jpg) Results (value of $\delta$) for clusters 50 through 2 are shown below. Although the 2 curves are close to each other, K-means' values are somewhat better. That is, K-means is more optimal than Ward. ![enter image description here](http://i.stack.imgur.com/mVzaT.jpg) The last picture plots $\delta^{Ward}/\delta^{Kmeans}$. The trend is clearly seen: as the number of clusters diminishes (that is, step $q$ for Ward grows), Ward tends to become less optimal relative to K-means. ![enter image description here](http://i.stack.imgur.com/qe33Z.jpg) One tempting question would be whether the relative sub-optimality of Ward (the value of $\delta^{Ward}/\delta^{Kmeans}$) on these latest 49 steps of agglomeration would be higher if I analyzed in the same fashion 1830 data points instead of just 183.


title: Unequal variances but equal means

text: Assume the sample ${(x_i)}_{i=1}^{n_1} \sim_{\text{iid}} {\cal N}(\mu, \sigma_1^2)$ is independent of the sample ${(y_i)}_{i=1}^{n_2} \sim_{\text{iid}} {\cal N}(\mu, \sigma_2^2)$. What are the available methods to get a confidence interval about the common mean $\mu$ ? In my case I have $n_1=n_2$. I would be satisfied by an answer for this case but I'm also interested in the general case. **EDIT:** Ooops.. I finally found this paper with Google. Sorry.


title: Modeling success rate with gaussian distribution

text: In many papers I see data representing a rate of success (i.e a number between 0 and 1) modeled as a gaussian. This is clearly a sin (the range of variation of the gaussian is all of R), but how bad is that sin? Under what assumptions would you say it is tolerable?


title: Reasons besides prediction to build models?

text: Joshua Epstein wrote a paper titled "Why Model?" available at http://www.santafe.edu/media/workingpapers/08-09-040.pdf in which gives 16 reasons:   1. Explain (very distinct from predict)   2. Guide data collection   3. Illuminate core dynamics   4. Suggest dynamical analogies   5. Discover new questions   6. Promote a scientific habit of mind   7. Bound (bracket) outcomes to plausible ranges   8. Illuminate core uncertainties.   9. Offer crisis options in near-real time   10. Demonstrate tradeoffs / suggest efficiencies   11. Challenge the robustness of prevailing theory through perturbations   12. Expose prevailing wisdom as incompatible with available data   13. Train practitioners   14. Discipline the policy dialogue   15. Educate the general public   16. Reveal the apparently simple (complex) to be complex (simple) (Epstein elaborates on many of the reasons in more detail in his paper.) I would like to ask the community:   * are there are additional reasons that Epstein did not list?   * is there a more elegant way to conceptualize (a different grouping perhaps) these reasons?   * are any of Epstein's reasons flawed or incomplete?   * are their clearer elaborations of these reasons?


title: R fit restricted AR(p) model

text: I have a question about using R to fit an AR model. If we want fit a AR(p) model, the equation will be $Y_t = φ_1Y_{t-1} + φ_2Y_{t-2} + ... + φ_pY_{t-p} + Z_t$.   What about I only want to fit the model like $Y_t = φ_1Y_{t-1} + φ_{11}Y_{t-11} + Z_t$? ($Z_t$ is white noise).   Does anyone know how to fit this model in R?


title: assigning students to different specialities

text: We are building a system that helps a university and its potential students choose the right field of study.(excuse my study-related-vocabulary! I am not a native English speaker). Applicants have many fields to choose from(not the same number of fields for every student). Access to some fields is based on some conditions(which also differs depending on the field and the applicants previous studies). We need help to determine the fairest way to give each student the right field. A statistical solution, guidance, keywords to search for... Etc. An example could be helpful: Say we have option field1, field2, and field3. to access field1 the applicant should have studied previous1, previous2, or previous3 which is further conditioned by having an "A grade". field2 is accessible to every applicant. field3 is accessible to previous1 with grade B or more. Applicants should submit their choices in order of preferences depending on their previous studies achievements and the conditions of the fields to be studied. Any help is appreciated. And we are open to any further explanations.


title: When/why does the central tendency of a resampling simulation markedly differ from the observed value?

text: Should one always expect the central tendency (i.e., mean and/or median) of a bootstrapped sample to be similar to the observed value? In this particular case I have responses that are distributed exponentially for subjects across two conditions (I didn't run the experiment, I only have the data). I have been tasked with boot strapping the effect size (in terms of Cohen's d, the one-sample formula, i.e. $\bar{M_D}\over{s_D}$ where is the sample estimate of the population standard deviation. The forumla for this is provided in Rosenthal & Rosnow (2008) on pg 398, equation 13.27. They use $\sigma$ in the denominator because it is historically correct, however standard practice has misdefined d as using $s$, and so I follow through with that error in the above calculation. I have randomized both within participants (i.e. a participants RT may be sampled more than once) and across subjects (participants may be sampled more than once) such that even if participant 1 is sampled twice, their Mean RT in both samples is unlikely to be exactly equal. For each randomized/resampled dataset I recalculate d. In this case $N_{sim} = 10000$. What I'm observing is a trend for the observed value of Cohen's d to be typically closer to the 97.5th percentile of than to 2.5th percentile of simulated observed values. It also tends to be closer to 0 than the median of the bootstrap (by 5% to 10% of the density of the simulated distribution). What can account for this (keeping in mind the magnitude of the effect I'm observing)? Is it due to it being 'easier' upon resampling to obtain more extreme variances than those observed relative to the extremity of means upon resampling? Might this be a reflection of data that has been overly massaged/selectively trimmed? Is this resampling approach the same as a bootstrap? If not, what else must be done to come up with a CI?


title: Which is the current testing methodology for pool-based active learning?

text: I have seen many papers using 10-fold CV ("10-pool" CV), but I think the accuracy obtained this way can be, sometimes, optimistically incorrect since at any given time step **t** the sum of distinct queries of all folds might be greater than **t**. **t** is normally used in plots as the number of queries (seen instances). Another problem is the impossibility to know how many queries were used, expect by averaging across the folds. Another approach commonly used is the Holdout 50/50, 70/30 etc. Hold-out does not have those specific problems, but is inadequate to small datasets.


title: Approximating a matrix as the outer product of two vectors

text: I have a two-dimensional data set that I would like to model as the outer product of two vectors (or two unit vectors times a constant). I need to find a way to a) figure out what these two vectors are and b) how well the data fits the model. The caveat is that I need to do this in the presence of possible missing values. I know I can use PCA to find a rank-one matrix approximating my data as well as how good that approximation is, if the matrix is complete. I'm not sure how to do it with missing values. This seems to me like some kind of degenerate form of collaborative filtering, with only one parameter per row/column. I think there has to be an easier way to solve it, though.


title: Appropriate GLM when response variable is proportion, but not binomial

text: The response variable I'm dealing with is the proportion of a total area that is suitable habitat for a species of interest. So although the response variable is bounded between 0 and 1, my intuition is that it wouldn't be appropriate to call it binomial since the numerator and denominator of the proportion are non-integer. The beta distribution comes to mind, but I'm uncertain of the appropriate link function and whether there tools in R to deal with a beta. Some background on my eventual goal: I'll likely be pursuing a conditional autoregressive model to account for spatial autocorrelation. I'll be treating space as 1-dimensional since I'm dealing with a river system, and so each observational unit only has two neighbors: one upstream and downstream. I'll be working in R and JAGS/BUGS if I decide to go Bayesian.


title: Forward Sampling for HMM

text: Is it reasonable to use **forward sampling** to compute the probability of P(X_1=x_1, ..., X_N=x_n) in an HMM where is the observation variable? Is the forward sampling algorithm related to the forward-algorithm for HMM?


title: How do Restricted Boltzmann Machines work?

text: I don't have a strong math background, so I am kind of stuck trying to understand how RBMs work. On wikipedia page about RBMs, I am able to understand the formulas in isolation, however, I'm not sure how these formulas are connected together. More specifically, isn't there cyclic dependency between energy function, probability functions & sampling process? In the training algorithm, it states to use Gibbs sampling together with other stuff in order to build/train W. However, sampling depends on P, and P by definition depends on E, and E depends on W, so in my noob eyes then this is cyclic dependency. Could someone please clarify with a short example how I could break this "loop"? Thanks http://en.wikipedia.org/wiki/Restricted_Boltzmann_machine


title: Moments of the two-parameter generalized Pareto distribution (GPD) needed

text: In this thread the first two moments of the two-parameter GPD are given, where the distribution might be defined as $G(y)= \begin{cases} 1-\left(1+ \frac{\xi y}{\beta} \right)^{-\frac{1}{\xi}} & \xi \neq 0 \\ 1-\exp\left(-\frac{y}{\beta}\right) & \xi=0 \end{cases}$ Now I require the formula for the first four moments and the conditions under which they exist / are finite. Wikipedia didn't help since only other types of the GPD are discussed here. The same holds for this page.


title: How to interpret glmer results (variance, correlation and ICC)

text: I'm a beginner in statistics and I have to run multilevel logistic regressions. I am confused with the results as they differ from logistic regression with just one level. I don't know how to interpret the variance and correlation of the random variables. And I wonder how to compute the ICC. For example : I have a dependent variable about the protection friendship ties give to individuals (1 is for individuals who can rely a lot on their friends, 0 is for the others). There are 50 geographic clusters of respondant and one random variable which is a factor about the social situation of the neighborhood. Upper/middle class is the reference, the other modalities are working class and underprivileged neighborhoods. I get these results :               > summary(RLM3)     Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']      Family: binomial  ( logit )     Formula: Arp ~ Densite2 + Sexe + Age + Etudes + pcs1 + Enfants + Origine3 +      Sante + Religion + LPO + Sexe * Enfants + Rev + (1 + Strate |           Quartier)        Data: LPE     Weights: PONDERATION     Control: glmerControl(optimizer = "bobyqa")               AIC      BIC   logLik deviance df.resid        3389.9   3538.3  -1669.9   3339.9     2778           Scaled residuals:          Min      1Q  Median      3Q     Max      -3.2216 -0.7573 -0.3601  0.8794  2.7833           Random effects:      Groups   Name           Variance Std.Dev. Corr             Neighb. (Intercept)     0.2021   0.4495                             Working Cl.    0.2021   0.4495   -1.00                     Underpriv.     0.2021   0.4495   -1.00  1.00     Number of obs: 2803, groups:  Neigh., 50          Fixed effects:      The differences with the "call" part is due to the fact I translated some words. I think I understand the relation between the random intercept and the random slope for linear regressions but it is more difficult for logistics ones. I guess that when the correlation is positive, I can conclude that the type of neighborhood (social context) has a positive impact on the protectiveness of friendship ties, and conversely. But how do I quantify that ? Moreover, I find it odd to get correlation of 1 or -1 and nothing more intermediate. As for the ICC I am puzzled because I have seen a post about lmer regression that indicates that intraclass correlation can be computed by dividing the variance of the random intercept by the variance of the random intercept, plus the variance the random variables, plus the residuals. But there are no residuals in the results of a glmer. I have read in a book that ICC must be computed by dividing the random intercept variance by the random intercept variance plus 2.36 (pi²/3). But in another book, 2.36 was replaced by the inter-group variance (the first level variance I guess). What is the good solution ? I hope these questions are not too confused. Thank you for your attention !


title: Are pupils in a school independent and identically distributed when compared to a national population?

text: I am interested in the way in which data are analysed in education. In England, pupils' reported test results are compared to national test results using z-tests. These results are then used to judge the effectiveness of schools. My hypothesis is that students in a school are not IID, and that therefore analyses using tests based on the Central Limit Theorem which seeks to judge a 'school effect' rather than a 'pupil effect' is not justified. Is it reasonable to assume that students and therefore their test scores are independent and identically distributed?


title: How to make a good color intensity scale?

text: I am by no means good in statistics, but I think I have come to the right place. My question is simple: My problem consists of comparing the population of several states in a small country, but some states have a population of 3000,000 and some a population of 2,000.   I am painting it on a map, and the "intensity" of the color depends on how the population of every state compares to the population of the whole country. The problem is that the states with a lot of population are shown with really intense colors and the small states barely have any color. Is there an easy way to "normalize" or make the data comparable? I dont know if I am explaining myself properly but I hope some one can help me. Please comment if my question is not clear and I will clarify. Thank you for your help!


title: The positive stable distribution in R

text: Positive stable distributions are described by four parameters: the skewness parameter $\beta\in[-1,1]$, the scale parameter $\sigma>0$, the location parameter $\mu\in(-\infty,\infty)$, and the so-called index parameter $\alpha\in(0,2]$. When $\beta$ is zero the distribution is symmetric around $\mu$, when it is positive (resp. negative) the distribution is skewed to the right (resp. to the left). Stable distributions allow fat tails when $\alpha$ decreases. When $\alpha$ is strictly less than one and $\beta=1$ the support of the distribution restricts to $(\mu,\infty)$. The density function only has a closed-form expression for some particular combinations of values for the parameters. When $\mu=0$, $\alpha<1$, $\beta=1$, and $\sigma=\alpha$ it is (see formula (4.4) here): $f(y) = -\frac{1}{\pi y} \sum_{k=1}^{\infty} \frac{\Gamma(k\alpha+1)}{k!} (-y^{-\alpha})^k \sin(\alpha k \pi)$ It has infinite mean and variance. **Question** I would like to use that density in R. I use               > alpha <- ...     > dstable(y, alpha=alpha, beta=1, gamma=alpha, delta=0, pm=1)      where the dstable function comes with the fBasics package. Can you confirm this is the right way to compute that density in R? Thank you in advance! **EDIT** One reason why I am suspicious is that, in the output, the value of _delta_ is different from that in the input. Example:               > library(fBasics)     > alpha <- 0.4     > dstable(4, alpha=alpha, beta=1, gamma=alpha, delta=0, pm=1)     [1] 0.02700602     attr(,"control")        dist alpha beta gamma    delta pm     stable   0.4    1   0.4 0.290617  1


title: Updating a set of estimated forecasts

text: Suppose I have some stochastic process $X_t$. At each time $t$, I receive an estimated probability distribution for $x_t$, followed by an observation $x_t$. After receiving a set of observations ${x_1, \ldots, x_n}$, I want to go back and re-estimate the probability distribution for each $x_t$, $1 \le t \le n$. What are some ways of doing this? Let's make the assumption that the estimated probability distributions are "pretty good", so that we do want to use the information they contain; what other assumptions do I need in order to make this an interesting/tractable problem? For example, some (potentially totally misguided) thoughts I had were:   1. Suppose $X_t$ is a "first-order continuous" model: $x_t$ depends only on $x_{t-1}$, and if $x_{t-1}, x_{s-1}$ are fairly close, then the probability distributions for $x_t, x_s$ should also be fairly close. Then to revise my estimate of the probability distribution for $x_t$, I could use a kind of kernel density bootstrap sample and mix it with the original estimate: take all points $x_s$ such that $x_{s-1}$ is close to $x_{t-1}$, use these points $x_s$ in some weighted fashion to form a kernel density estimate $K$, and average $K$ with the original estimated probability distribution to get a new estimated probability distribution for $x_t$.   2. Suppose $X_t$ is a linear Gaussian model (so that it satisfies the assumptions of a standard Kalman filter). I could try forming a different kind of bootstrap sample: for each $x_t$, use a window around its estimated probability distribution to sample a new observation $x'_t$, and then use Kalman smoothing on the observations ${x'_t}$ to get a set of smoothed observations $\tilde{x_t}$. Use these bootstrap samples to generate new estimated probability distributions. This is a mostly random question a friend asked me a while ago that I just thought of again, so I don't have a particular application (and hence why I don't have a concrete set of assumptions on the process).


title: Test equivalence of non-nested models

text: Let's say $y$ is a linear function of $x$ and a dummy $d$. My hypothesis is that $d$ itself is like a hedonistic index of a vector of other variables, $Z$. I have support for this in a $MANOVA$ of $Z$ (i.e. $z_1$, $z_2$, ..., $z_n$) on $d$. Is there any way to test the _equivalence_ of these two models: Model 1: $y = b_0 + b_1 \cdot x + b_2\cdot d + e_1$ Model 2: $y = g_0 + Z\cdot G + e_2$ where $G$ is the column vector of parameters.


title: How to calculate pseudo-$R^2$ from R's logistic regression?

text: Christopher Manning's writeup on logistic regression in R shows a logistic regression in R as follows:               ced.logr <- glm(ced.del ~ cat + follows + factor(class),        family=binomial)      Some output:               > summary(ced.logr)     Call:     glm(formula = ced.del ~ cat + follows + factor(class),         family = binomial("logit"))     Deviance Residuals:     Min            1Q    Median       3Q      Max     -3.24384 -1.34325   0.04954  1.01488  6.40094          Coefficients:                   Estimate Std. Error z value Pr(>|z|)     (Intercept)   -1.31827    0.12221 -10.787 < 2e-16     catd          -0.16931    0.10032  -1.688 0.091459     catm           0.17858    0.08952   1.995 0.046053     catn           0.66672    0.09651   6.908 4.91e-12     catv          -0.76754    0.21844  -3.514 0.000442     followsP       0.95255    0.07400  12.872 < 2e-16     followsV       0.53408    0.05660   9.436 < 2e-16     factor(class)2 1.27045    0.10320  12.310 < 2e-16     factor(class)3 1.04805    0.10355  10.122 < 2e-16     factor(class)4 1.37425    0.10155  13.532 < 2e-16     (Dispersion parameter for binomial family taken to be 1)     Null deviance: 958.66 on 51 degrees of freedom     Residual deviance: 198.63 on 42 degrees of freedom     AIC: 446.10     Number of Fisher Scoring iterations: 4      He then goes into some detail about how to interpret coefficients, compare different models, and so on. Quite useful. However, how much variance does the model account for? A Stata page on logistic regression says: > Technically, R2 cannot be computed the same way in logistic regression as it > is in OLS regression. The pseudo-R2, in logistic regression, is defined as > (1 - L1)/L0, where L0 represents the log likelihood for the "constant-only" > model and L1 is the log likelihood for the full model with constant and > predictors. I understand that at the high level. The constant-only model would be without any of the parameters (only the intercept term). Log likelihood is a measure of how closely the parameters fit the data. In fact, Manning sort of hints that the deviance might be -2 log L. Perhaps null deviance is constant-only and residual deviance is -2 log L of the model? However, I'm not crystal clear on it. Can someone verify how one actually computes the pseudo-R^2 in R using this example?


title: Covariance and independence?

text: I read from my textbook that $\text{cov}(X,Y)=0$ does not guarantee X and Y are independent. But if they are independent, their covariance must be 0. I could not think of any proper example yet; could someone provide one?


title: Generating random samples from a custom distribution

text: I am trying to generate random samples from a custom pdf using R. My pdf is: $$f_{X}(x) = \frac{3}{2} (1-x^2), 0 \le x \le 1$$ I generated uniform samples and then tried to transform it to my custom distribution. I did this by finding the cdf of my distribution ($F_{X}(x)$) and setting it to the uniform sample ($u$) and solving for $x$. $$ F_{X}(x) = \Pr[X \le x] = \int_{0}^{x} \frac{3}{2} (1-y^2) dy = \frac{3}{2} (x - \frac{x^3}{3}) $$ To generate a random sample with the above distribution, get a uniform sample $u \in[0,1]$ and solve for $x$ in $$\frac{3}{2} (x - \frac{x^3}{3}) = u $$ I implemented it in `R` and I don't get the expected distribution. Can anyone point out the flaw in my understanding?               nsamples <- 1000;     x <- runif(nsamples);          f <- function(x, u) {        return(3/2*(x-x^3/3) - u);     }          z <- c();     for (i in 1:nsamples) {       # find the root within (0,1)        r <- uniroot(f, c(0,1), tol = 0.0001, u = x[i])$root;       z <- c(z, r);     }


title: Calculate how big a variable is in relation to others

text: This will sound weird, but I'm trying to calculate something that I don't know the name. I want to calculate how big/low a variable is in relation to the others. Like Google webmaster tools says "your site loads in 3s, this is slower than 59% of sites". How can I calculate this? For example, the points a player has. He has 4440 points and he is ranked on #437 out of 16801 players. The average of points in 976 and the top player has 30076 points.


title: Regression to the mean puzzle

text: In the "Regression to the Mean" chapter of "Thinking, Fast and Slow" by Daniel Kahneman, an example is given and the reader is asked to forecast the sales of individual stores given the overall sales forecast and the sales numbers from the previous year. For instance(the book's example has 4 stores, I use 2 here for simplicity):               Store    2011    2012     1        100      ?     2        500      ?     Total    600     660      The naive forecast would be 110 and 550 for stores 1 and 2, 10% increase for each. However, the author claims this naive approach is wrong. It is more likely for the poorer-performing store to increase more than 10%, and the better-performing store to increase(or even decrease) by less than 10%. So perhaps a forecast of 115(15% increase) and 535(7% increase) would be "more correct" than the naive forecast. What I don't understand is how we can conclude that sales of 100 of store 1 is necessarily the poorer-performing store? Perhaps, due to location differences, the true time-series means of stores 1 and 2 are 10 and 550, and store 1 had a super year in 2011, and store 2 had a disastrous year in 2011. Then wouldn't it make sense to forecast a decrease for store 1 and increase for store 2? I know that time series information was not given in the original example, but I am under the impression that "regression to the mean" refers to the cross- sectional mean and therefore time-series information does not matter. What am I misunderstanding?


title: Propagation of Poisson Confidence Intervals

text: I have a set of measurement with 95% Poisson confidence intervals, and I would like to subtract and multiply them and propagate the error. For example, I measured the copy number of a piece of DNA in a biological sample as 0.0268, 95% CI = 0.308-0.227. I know there are systematic biases in my measurement, ie issues with sensitivity and specificity, which I independently determined. For example, to measure sensitivity I spiked in a known concentration of copies of DNA, 467 copies per ul, and I measured 251 copies per ul, 95% CI 255-248, so arround 54%. To measure specificity I used a known negative DNA sample, so expected 0 copy number but measured 0.00162, 95% CI 0.0023-0.000936.   I would like to multiply the measured copy number by the specificity, which would give a value of ~.496, then subtract the specificity, which would yield a final copy number estimate of ~.494. How can I propagate the errors in this case? Thanks for any help!


title: Using ccf() to find correlation between two data frames offset by nine months

text: This is the next step from a question I asked earlier. I've got two data frames: one focused on birth data, and one focused on winter weather events. The aim of my project is to discover whether there exists a simple correlation between extreme winter weather events (i.e. winter storms) and a spike in births nine months later (due to people getting stuck in doors during the storms). There are several areas with which I'm struggling:   1. I need to factor out less extreme winter weather events from combined.weather.birth$EVENT_TYPE. The factors currently included are "Frost/Freeze", "Hail", "Heavy Snow", "Ice Storm", "Winter Storm", "Winter Weather", and "Blizzard". **Of these, I wish to exclude frost/freeze and hail.**   2. I'm having difficulty running the cff() function on this data. As described above, I want to discover and analyze potential correlations between these data sets. I'm only comparing data in Massachusetts from years 2007-2011. Here is what I've tried so far:               correlation1 <- ccf(birth.data$DATE, combined.weather.birth$DATE, lag.max = NULL, type="correlation", plot=TRUE)          correlation2 <- ccf(birth.data$DATE, combined.weather.birth$DATE+combined.weather$EVENT_TYPE, lag.max = NULL, type="correlation", plot=TRUE)      **I need to offset this data by nine months, to account for pregnancy after the winter weather events. Any tips?** Here is some information on the data I'm working with:               str(combined.weather.birth) <-         'data.frame':   966 obs. of  8 variables:      $ EVENT_ID       : int  9620 9619 9623 5391 13835 13845 13844 13847 13846 13836 ...          $ STATE          : Factor w/ 1 level "MASSACHUSETTS": 1 1 1 1 1 1 1 1 1 1 ...      $ YEAR           : int  2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...          $ MONTH_NAME     : Factor w/ 12 levels "April","August",..: 5 5 5 5 4 4 4 4 4 4 ...      $ EVENT_TYPE     : Factor w/ 7 levels "Frost/Freeze",..: 6 6 4 4 5 5 5 5 5 5 ...          $ INJURIES_DIRECT: int  0 0 0 1 0 0 0 0 0 0 ...      $ DEATHS_DIRECT  : int  0 0 0 0 0 0 0 0 0 0 ...          $ DATE           : POSIXct, format: "2007-01-01" "2007-01-01" "2007-01-01" "2007-01-01" ...          str(birth.data) <-         'data.frame':   60 obs. of  4 variables:      $ YEAR       : int  2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...          $ MONTH      : Factor w/ 12 levels "April","August",..: 5 4 8 1 9 7 6 2 12 11 ...      $ BIRTH_TOTAL: num  6250 5833 6570 6227 6858 ...          $ DATE       : POSIXct, format: "2007-01-01" "2007-02-01" "2007-03-01" "2007-04-01" ..      EDIT: I should add that I'm not married to using cff() here. If there is a better function for finding the specified correlation, I am open to learning about it. I've read a bit about cor(), but that doesn't seem appropriate here since it's designed to only work with matrices.


title: Is there a more elegent solution for exponentially weighted-mean in R?

text: OK so I have a working code to do what I want but I am new to R and feel like my solution is very clunky and there is likely a more efficient way to get the same result. I have a set of data that contain records with a date, city, and quality score. For each city and each date, I want to create an exponentially weighted-average of that city's prior quality scores. (The first date for each city has no previous data, so it will always be NA.) My solution expands the original dataframe and calculates the means in a single pass. It does work but I view it as a bit inefficient and expanding the dataframe could be problematic with large data. As I am new to R, I would love to hear other approaches that get the same end result. Here is my code:               #load in the data     mydata <- structure(list(date = structure(c(15513, 15507, 15476, 15439, 15442, 15435, 15419, 15410, 15508, 15464, 15461), class = "Date"), city = structure(c(2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 1L, 1L, 1L), .Label = c("MIA", "PHX", "POR"), class = "factor"), quality = c(3.5, 6.7, 5.2, 3.9, 2.6, 7.7, 2.4, 4.7, 3.5, 2.6, 1.8)), .Names = c("date", "city", "quality"), row.names = c(NA, -11L), class = "data.frame")          #fix the date data     mydata$date <- as.Date(mydata$date, "%d%b%Y")          #create a unique id for each record     mydata$id <- paste(mydata$city, as.numeric(mydata$date), sep='_')          #expand the data     mydata_expanded <- expand.grid(x=mydata$id, y=mydata$id)     mydata_expanded$cityX <- substr(mydata_expanded$x,start=1,stop=3)     mydata_expanded$cityY <- substr(mydata_expanded$y,start=1,stop=3)     mydata_expanded <- mydata_expanded[mydata_expanded$cityX==mydata_expanded$cityY,]     mydata_expanded$dateX <- substr(mydata_expanded$x,start=5,stop=10)     mydata_expanded$dateY <- substr(mydata_expanded$y,start=5,stop=10)     mydata_expanded <- merge(mydata_expanded, mydata, by.x=c("y"), by.y=c("id"))     mydata_expanded$dateX <- as.Date(as.numeric(mydata_expanded$dateX), origin="1970-01-01")     mydata_expanded$dateY <- as.Date(as.numeric(mydata_expanded$dateY), origin="1970-01-01")     mydata_expanded <- mydata_expanded[mydata_expanded[,5]>mydata_expanded[,6],]     mydata_expanded <- mydata_expanded[order(-as.numeric(mydata_expanded$dateX)),]          #create the variables needed to calculate the exponentially weighted-means     mydata_expanded$dayssince <- mydata_expanded$dateX - mydata_expanded$date         mydata_expanded$w <- exp(-as.numeric(mydata_expanded$dayssince)/60)         mydata_expanded$wQuality <- mydata_expanded$quality * mydata_expanded$w     sum_w <- aggregate(mydata_expanded$w,by=data.frame(mydata_expanded$dateX, mydata_expanded$city),sum,na.rm=T);         sum_wQuality <- aggregate(mydata_expanded$wQuality,by=data.frame(mydata_expanded$dateX, mydata_expanded$city),sum,na.rm=T)     mydata_quality <- merge(sum_w,sum_wQuality,by=c("mydata_expanded.dateX", "mydata_expanded.city"))          #calculate the actual exponentially weighted-mean     mydata_quality$ewQuality <- mydata_quality$x.y/mydata_quality$x.x          #merge the data back to the original dataframe and cleanup     mydata <- merge(mydata,mydata_quality, by.x=c("date","city"), by.y=c("mydata_expanded.dateX","mydata_expanded.city"),all.x=TRUE)     mydata$id <- NULL         mydata$x.x <- NULL     mydata$x.y <- NULL         mydata$ewQuality <- signif(mydata$ewQuality,digits=4)     mydata <- mydata[with(mydata, order(city)), ]


title: statistical error on difference f(x) - x

text: Suppose I've a random variable $X$ and a sample of it with size $N$. I count how many element of the sample fall in a specific range $a<x<b$ and I found $n$ entries, so if I use Poisson statistics I can say that the error on $n$ is $\sqrt{n}$. Now suppose I've a function $f$, quite close to the unit function, and I define a new random variable $X'=f(X)$. Now the number of entries in $a<x<b$ is $n'=n+\delta n$. I don't have an analitical from of $f$. I only now $n$ and $n'$. What is the error on $n'-n$? How to take into account the correlation between them? More difficult: suppose $f$ depends on a parameter $\lambda$ with error $\delta\lambda$.


title: What is a dimensionless indicator?

text: What is a dimensionless indicator? I assume it refers to the fact that it is consistent in the way that it doesn't vary disproportionately at the variation of its factors.


title: Fitting moments to a distribution

text: I have data on the first to fourth moments of a continuous random variable and I am trying to find what distribution best fits the data. Wikipedia has a list of about 20 distributions that could fit the data and I would like to try each of them to see which is the best. Some of these distributions have one parameter, some have two or three. Is there some test I can use to compare the fit of each distribution? Or at least a test to compare distributions with the same number of parameters?


title: How do I use math to predict the next number in the series?

text: Here's a series of data I'm observing:               1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1      How do I use math to predict whether the next number in the series will be a 1 or a 0?


title: Link between forward and inverse regression ($\text{E}(X|Y)$ and $\text{E}(Y|X)$ ;$ \text{var}[\text{E}(X|Y)]$ and $\text{var}[\text{E}(Y|X)]$)

text: In a multivariate context, that is with at least X or Y being a random vector, are there formulae or theorems that link (even remotely) the forward and inverse regression, $\text{E}(X|Y)$ and $\text{E}(Y|X)$ ? Or alternatively $\text{E}(X_i|Y)$ and $\text{E}(Y|X_i)$, where $X_i$ is a given component of $X$. Similarly, can something be said about $\text{var}[\text{E}(X|Y)]$ and $\text{var}[\text{E}(Y|X)]$ ? Or alternatively $\text{var}[\text{E}(X_i|Y)]$ and $\text{var}[\text{E}(Y|X_i)]$, where $X_i$ is a given component of $X$. Both $X$ and $Y$ are supposed to be random. The typical application case is that $X$ follows a chosen distribution (for instance a multigaussian) and $Y$ is a function of $X$. The link function is too complicated to be used or even unknown (it is typically a numerical model). The distribution of $X$, however, is controlled. Thus, one should formulate hypothesis on $X$ if needed. *Note : this question is related to dimension reduction trhough the sliced inverse regression technique. It is also related to sensitivity analysis, as $\text{var}[\text{E}(Y|X)]$ is the numerator of Sobol' index (surprisingly, the tag `sensitivity-analysis` does not exist -- if someone with sufficient reputation read this, he might want to create it).*


title: How to draw a side-by-side plot mentioned in "Graphical Display as an Aid to Analysis"

text: Emerson, J. D. (1991) Graphical Display as an Aid to Analysis, in _Fundamentals of Exploratory Analysis of Variance_ (eds D. C. Hoaglin, F. Mosteller and J. W. Tukey), John Wiley & Sons, Inc., Hoboken, NJ, USA. doi: 10.1002/9780470316832.ch8 There describes side-by-side plots of effects and residuals... I think it's a nice idea for exploratory analysis, but don't know how to draw it.. (The usual side-by-side boxplot is good for one-way analysis but when it gets to two way, it's too complicated to see main effects and interactions...) Does anyone know how? Here's some example from the excellent book "A First Course in Design and Analysis of Experiments", which is available on web!!! ![one-way](http://i.stack.imgur.com/roTYy.jpg) ![more complicated setting](http://i.stack.imgur.com/4mQVG.jpg)


title: Removing a 'trend in variance' from a time series

text: I'm looking at a time series which has a very strong daily cycle in it. However, on top of having a daily cycle in the actual values of the time series, it also has a very strong daily variance cycle. I am wondering if I can meaningfully remove the 'variance trend' from my time series. I can calculate the 'average variance' for a given hour and then divide by it, but I'm not sure I'm doing anything useful if I do this? The large variance does correlate with large values in the daily cycle too, but taking a log doesn't seem to help much. EDIT: EXTRA INFORMATION As pointed out in the comments, this problem is rather under-defined. Unfortunately, this is a symptom I am also dealing with and so I can't fix that. However, I will try to give some more information. The **goal** is to have a model of the time series from which realisations can be drawn. The time series data is a measure of concentration over time, with clear daily and yearly cycles and a linear trend. There are many covariate time series which at least follow the same daily and yearly cycles. I am currently just trying to explore and learn different time series techniques. I wasn't hoping for a 'solution' but rather for a bit of direction as I hadn't been able to find much information about trends in the variance.


title: Method to jointly estimate Logit and OLS models

text: I am working on the following two models. Logit: Binary outcome $y$ ($y=1$ if paid more than the item value to win the item, $y=0$ otherwise) $$ y=b_0+x_1b_1+x_2b_2+x_3b_3. $$ OLS: Continuous dependent variable $z$ (amount of overpayment) $$ z=b_0+x_1b_1+x_2b_2+x_3b_3+e. $$ My questions:   1. Do I need to estimate these two models jointly (instead of separately) since the dependent variables seem to be related?   2. If joint estimation is a must, which method (or package in Stata) should I use? Seemingly unrelated regression `suest`?   3. Any other thoughts in terms of methodology?


title: How to handle high dimensional feature vector in probability graph model?

text: I was doing some NLP related stuff which involves training a hidden Markov model, and use the model to segment sentences. For every sentence, I translate the tokens into feature vectors. The features are manually picked by me, and I can only think of 20 features temporarily. All of the features are binary. So an example feature vector in my case would be like:               [0, 0, 1, 1, 0…]       When training the HMM (supervised learning with maximum likelihood estimation), I convert the binary feature vector to integer, and use the integers as the observations in HMM. My question is what if I have more dimension of features, say 100 features, then it might not be possible or efficient to map feature vectors to integers when training the HMM. What if the features are not binary. What are some common solutions/best practice to these problems? P.S. I did some research and find that one solution is to use K-Means to find clusters of feature vectors and consider each cluster as a single observation in HMM (also known as vector quantization). But I don't think it is the best solution.


title: How do I calculate confidence intervals for a linear combination of coefficients of a logistic regression?

text: I want to test a series of hypotheses on linear combinations of coefficients for logistic regression on categorical count data implemented as a GLM model. I know how to do linear combinations under normality, and I know how to do basic confidence intervals and significance tests for logistic models. But I just realized with some embarrassment that knowing these two things does not add up to knowing how to construct significance tests and confidence intervals on linear combinations of logistic regression coefficients. **How do you construct significance tests and confidence intervals on linear combinations of logistic regression coefficients?**


title: Robust estimation of mean with small sample size

text: Working with a 2D sensor. A row consists of 4096 samples, of which there are 16 _reference_ samples (the rest are _active_ ). The reference samples do not measure the signal, but are used to estimate the unique measurement offset that affects the entire row. The measurements must be corrected by this row offset before being further processed. This row offset varies between successive rows, and also over time. Each sample is independently affected by measurement noise which seems to follow a Gaussian distribution. So a given measurement $M$ is given by the signal at that location $S$ plus the row offset $O$ plus some Gaussian noise $G$, thus: $$M_{xy} = S_{xy} + O_{y} + G_{xy}$$ The typical procedure is to use the mean of these reference measurements to estimate this offset and subtract it from the active measurements (the remaining samples in that row) to obtain a corrected measurement: $$O_{y} = \Sigma M_{ref}/N$$ when $N=16, S=0$ (in the reference samples). Of course, when $S=0$, $M=O_{y} + G_{xy}$ so we still need to filter out the noise to get the offset. Currently, using the mean we can get the noise (stddev) quite low, but there is still sufficient noise that there are noticeable artefacts in the results. So we are looking for a more robust way (than mean or median) of estimating the offset from these samples. If it helps, during the calibration procedure, we are able to obtain "zero" level readings from the entire row. Thus under lab conditions (and _not_ in the field), we can read out the zero level from the entire row. So the question is:   * How can you robustly estimate the mean value of a small ($N=16$) noisy set of samples that are representative of a much larger population? and the bonus question:   * Is there a way to exploit the "zero" measurements available only during calibration? This answer seems closest to this question, but so far the methods described don't seem to yield a significant improvement:   * Crash course in robust mean estimation


title: How does the Kinect use random forests?

text: I read on this site that apparently the Kinect uses the random forests algorithm for machine learning in some way. Can anyone explain what it uses random forests for, and how their approach works?


title: Time series analysis in Python

text: I am a beginner to time-series analysis. I have the model below; y is sales of product and x is tweet-rate: $y_t=ay_{t-1}+by_{t-2}+...+cy_{t-m}+dx_t+ex_{t-1}+...+fx_{t-n}$   1. What is this model called? I guess it's called an AR model but I am not sure since the dependent variable y is on R.H.S as well.    2. How do I fix the lag period, $m$ and $n$? Can $x$ and $y$ have different lags?   3. How can I use Python to build this model and also predict the sales for $t+1\ldots t+n$? Any solution for this without using rpy.


title: Generalization for interpreting a test?

text: At my college I am doing a lot of testing, such as testing for autocorrelation, White , DW, $t$, $F$ etc. However, I always have to look up, when to reject and when to accept for example $H_0$ Hence, I was wondering if there is a generalization for interpreting a statistical test such as: When $value$ is greater than $value$ always _reject_ H0. From the point of view of your professional experience, is there a pattern in all statistical tests?


title: Variance decomposition using ANOVA

text: Say you have a random variable $X$ (e.g., kilometers driven). Getting its variance is straightforward. But what if you want to say, $A$ percent of the variance in $X$ is due to $\text{Var}(X)$ for female drivers and $B$ percent is the rest, that is, $\text{Var}(X)$ for male drivers? $A + B$ should be 100 percent. Is this possible? Are there assumptions to be made to simplify things? Independence of female and male drivers? [Question was also asked here http://mathoverflow.net/questions/88185/variance-decomposition-anova but I realized it's best to do it at Stats instead.]


title: Can p-values for Pearson's correlation test be computed just from correlation coefficient and sample size?

text: **Background:** I read one article where authors report Pearson correlation 0.754 from sample size 878. Resulting p-value for correlation test is "two star" significant (i.e. p < 0.01). However, I think that with such a large sample size, corresponding p-value should be less than 0.001 (i.e. three star significant).   * **Can p-values for this test be computed just from Pearson correlation coefficient and sample size?**   * **If yes, how can this be done in R?**


title: Does anybody use star-schema databases to collect and organize their data?

text: I've been reading about the star schema (or dimensional) database structure, which puts all measurements main in 'facts' table, and all context for those measurements in 'dimension' table linked to the facts table (I'm doing a horrible job explaining this), so that you can query any measurement by conditioning on the dimension. When I read about this, it struck me as a perfect way to aggregate different csv files with different columns into one store... for example, stock, commodities, and options data in various files go be shoved into a star-schema db with a facts table containing the open/close/high/low/volume, with foreign- key links to dimensions such as time, company, exchange, and 'file pedigree' so that we know where each entry comes from. That would allow you to query for, say, first quarter chinese mining companies with a single, fairly simple sql statement. However, all of the literature on this topic seems to be for large scale 'enterprise' deployments and not small data analysis or research jobs. Does anyone have any experience organizing their data this way? Is there are simpler way?


title: How to evaluate goodness of fit for negative binomial regression

text: I'm trying to fit a model estimating waiting time using negative binomial regression, but I'm not sure how to assess the goodness of fit for my model. I would like to compare the negative binomial model to a Poisson model. I have approximately $4,000$ data points. Any suggestions? Thanks!


title: Can longitudinal and cross-sectional studies be "controlled"?

text: I have heard it said that longitudinal and cross-sectional studies are forms of observational studies. **Why can't longitudinal and cross-sectional studies be considered "controlled studies"?** By their definitions, I think they can.


title: Loss for Kernel Ridge Regression

text: Is $||Y-X\beta||_2^2 + \lambda\beta^T K\beta$ , the standard loss-function in kernel ridge regression, or is it different? Also, is the gaussian kernel a standard choice used for the kernel, in practice? If not, which kernels are used more often than not? Also, is $\lambda$ the only parameter to be tuned via cross-validation or is the kernel parameter like $\sigma$ in a gaussian kernel, also tuned via cross validation in practice? Please confirm and/or correct my understanding of Kernel ridge regression!


title: Clustering-based Identification

text: I am trying to find a way to merge data that correspond to users and have the form $user_i, property_1, property_2, .\dots, property_n$, where $n\geq50$, $i\geq 10^6$ We might have missing data in some cases (i.e. $property_k$ might be missing for some users), or we might have two or more users (lines in the database) that correspond to the same user, and the only difference is some small variation in some of their properties. Specifically, I am interested in finding what is the cluster that a user belongs to, using some machine learning algorithm, etc. Each cluster will have all the instances of a user (all the lines that correspond to him). For that, we need of course some distance metric, which will depend on the model, I guess. Any ideas on what is the best model for such a system?


title: Variation explained by single variable

text: I’m trying to find a way to measure how much a single variable ‘summarizes’ a full set of continuous variables. For instance, in a PCA the first principal component will explain a certain percentage of the total variability in a multivariate set. So, how can I obtain a similar measurement for a pre- existing (untransformed) variable? For instance, how much does altitude (single variable) explains overall climatic variability (i.e. multiple variables: mean precipitation; mean temperature…)? I am particularly interested in a measurement that is directly comparable to the variance explained by a PCA. Let me know if you have any suggestions.


title: Why is a likelihood-ratio test distributed chi-squared?

text: Why is the test statistic of a likelihood ratio test distributed chi-squared? $2(\ln \text{ L}_{\rm alt\ model} - \ln \text{ L}_{\rm null\ model} ) \sim \chi^{2}_{df_{\rm alt}-df_{\rm null}}$


title: What is the procedure for "bootstrap validation" (a.k.a. "resampling cross-validation")?

text: "Bootstrap validation"/"resampling cross-validation" is new to me, but was discussed by the answer to this question. I gather it involves 2 types of data: the real data and simulated data, where a given set of simulated data is generated from the real data by resampling-with-replacement until the simulated data has the same size as the real data. I can think of two approaches to using such data types: (1) fit the model once, evaluate it many times on many simulated data sets; (2) fit the model many times using each of many simulated data sets, each time evaluate it against the real data. Which (if either) is best?


title: Interpretation of intercept in Helmert coding

text: I have a general, presumably simple, question, but I couldn't find a conclusive answer so far. Assume I have a simple case of a General Linear Model with one categorical predictor variable, that has 3 levels. This corresponds to a one-way ANOVA with 3 groups. The Linear Model would then be: $Y_i = \beta_0 + \beta_1X_1 + \beta_2X_2$ I use the following standard Helmert contrast codes (I, II) for the three groups A, B, C:                    A    B    C     I   -2    1    1     II   0   -1    1      Based on these contrast codes, **will the intercept of the linear model represent the mean of all group means? More generally, does this always apply to Helmert contrast coded variables,** regardless of the number of predictor variables or variable levels? Obviously, the intercept represents the predicted value, when all predictor variables (in this case $X_1$ and $X_2$) equal zero, but I can't figure out if that represents the mean of the group means.


title: R linear regression categorical variable "hidden" value

text: This is just an example that I have come across several times, so I don't have any sample data. Running a linear regression model in R:               a.lm = lm(Y ~ x1 + x2)      `x1` is a continuous variable. `x2` is categorical and has three values e.g. "Low", "Medium" and "High". However the output given by R would be something like:               summary(a.lm)                   Estimate Std. Error t value Pr(>|t|)         (Intercept)   0.521     0.20       1.446   0.19             x1            -0.61     0.11       1.451   0.17     x2Low         -0.78     0.22       -2.34   0.005     x2Medium      -0.56     0.45       -2.34   0.005      I understand that R introduces some sort of dummy coding on such factors (`x2` being a factor). I'm just wondering, how do I interpret the `x2` value "High"? For example, what effect does "High" `x2`s have on the response variable in the example given here? I've seen examples of this elsewhere (e.g. here) but haven't found an explanation I could understand.


title: Extending the Hellinger Distance to multivariate distributions

text: The hellinger distance for a univariate distribution is $$ \ H(x) = 1 - \int\ \sqrt {f(x)g(x)} dx $$ I wish to use it for a bivariate distribution, by extending it to this form $$ \ H(x) = 1 - \int\ \sqrt {f(x, y)g(x, y)} dx dy $$ Can someone comment on whether this is still a valid Hellinger distance?


title: Negative parameters of beta distribution

text: I have a set of observations of credit loss data, where the mean is 37% and variance 25%. Now, I have to find the distribution and the base assumption is it will follow a beta distribution. the issue is that my alpha and beta derived from mean and variance is being estimated at -0.025012 and -0.042588. I dont understand what to do with the negative values of alpha and beta. The formula which I am using to calculate alpha is mean*(((mean*(1-var))/var)-1) and beta is calculated by (1-mean)(((mean*(1-mean))/var)-1). Please do let me know how can I solve the problem.


title: What does the logit value actually mean?

text: I have a logit model that comes up with a number between 0 and 1 for many cases, but how can we interprete this? Lets take a case with a logit of 0.20 Can we assert that there is 20% probability that a case belongs to group B vs group A? is that the correct way of interpreting the logit value?


title: Is it ok to bin residuals before examining them?

text: I'm analyzing the residuals from a regression model fit to a dataset that covers several years worth of data. I want to report the sum of the residuals from that model, by year, as a measure of how the overall error for each year changes over time. Is this an acceptable way to report residuals? Here is an example of my calculation:               library(plyr)     set.seed(1)     Year <- c(  rep(2007,12),                 rep(2008,12),                 rep(2009,12),                 rep(2010,12),                 rep(2011,12))     n <- length(Year)     index <- 1:n     A <- (index+5*rnorm(n))^2+5*rnorm(n)     P <- predict(lm(A~Year),data.frame(A,Year))          residual <- A-P     residual2 <- ddply(data.frame(Year,residual),'Year',sum)          plot(2007+(index-7)/12,residual)     lines(residual2$Year,residual2$residual,col=2)          par(mfrow=c(2,1))     plot(index,residual)     plot(residual2)


title: Obtaining an estimator for z given an estimator for log z

text: As per gung's advice in Getting the equation from R's lm when using a product, I am starting a new thread for this question. I have a model $\widehat{\log z} = a + bx + cy + dxy$ for random variables $x,y,z$ and coefficients $a,b,c,d$. I now want a model for $z$ instead of $\log z$. I was going to simply do $\hat{z} = \exp(a + bx + cy + dxy)$ but I have been cautioned against this. Why is this incorrect, and what is the correct way to proceed (if there is one)?


title: How to compute logarithm of zero determinant of covariance matrix?

text: I am performing Box ´s test for equality of covariance matrices and I need to compute log(det(S1)) and log(det(S2)), where S1, S2 are sample covariance matrices with nearly zero determinant. How can I improve S1 and S2 in order to get non-zero determinant? I have 24-dimensional data with 14 observations. Thank you


I am studying a complex system. My goal is to understand the impact of a spreading accident (which saturates the filter (sand and plant) and creates puddles). Puddles decrease the efficiency of the system (depollution of effluent, and moisture content). The impact on the efficiency will decrease with time (more or less depending on the season and meteorological situation). I want to develop a model of how long the incident will impact on the pilote. For that I think I should use time series analysis. My question is:   * Can I use time series analysis on short time scales (e.g., with 5 to 30 measurements?   * or can I compare the effluent volume with the parameter that I want to study, with a decalage which should not be regular because of its depending with the effluent volume injected?