title: GLM model validation

text: When performing model validation and we are dropping the least significant explanatory variables until we find the optimum model where all remaining variables are significant, how does one go about dropping a particular "unimportant level" within a nominal variable? In my output below I have run a Neg. bin. GLM inorder to test abundance against the various explanatory variables listed. (as.factor(JJ)2) is the most obviously insignificant variable that I would like to drop but because it is a factor with three levels I don't know how to do this. (Of course I will have to rejig the data to unmask baseline values that are within the intercept also) Thanks!               Coefficients:                 Estimate Std. Error z value Pr(>|z|)         (Intercept)     3.351414   0.945861   3.543 0.000395 ***     DD              0.039417   0.011373   3.466 0.000528 ***     EE             -0.001497   0.001015  -1.475 0.140093         PP             -0.138303   0.094287  -1.467 0.142423         as.factor(BB)2 -0.901877   0.365796  -2.466 0.013682 *       as.factor(JJ)1  0.495985   0.198394   2.500 0.012420 *       as.factor(JJ)2 -0.276830   0.374576  -0.739 0.459878        as.factor(VV)2 -0.499569   0.174710  -2.859 0.004244 **      as.factor(VV)3 -0.166065   0.175430  -0.947 0.343834         ---     Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1           (Dispersion parameter for Negative Binomial(4.8369) family taken to be 1)              Null deviance: 97.392  on 46  degrees of freedom     Residual deviance: 47.883  on 38  degrees of freedom     AIC: 401.37


title: The distribution of a sufficient statistic

text: If I understand correctly, a distribution in the exponential family... $$\underline X\sim f_{\underline\theta}(\underline x) = exp\\{\sum\limits_{i}\eta_i(\underline\theta)T_i(\underline x)-B(\underline\theta)\\}~h(\underline x)$$ ...where $\underline\eta(\centerdot)$ is a (possibly vector valued) parameter, $\underline T(\centerdot)$ is a corresponding vector of sufficient statistics, $B(\theta)$ is the log partition, and $h(\underline x)$ is the base measure and $\underline T(\underline x)$ itself has the following distribution: $$\underline T(\underline x) =\underline t\sim g_{\underline\theta}(\underline t) = exp\\{\sum\limits_{i}\eta_i(\underline\theta)t_i-B(\underline\theta)\\}~h^*(\underline t)$$ ...where $h^*(\centerdot)$ is not necessarily the same function as the one from $f_{\underline\theta}(\underline x)$ above. My questions are:   * Are the other two functions ($B(\underline\theta)$ and $\underline\eta(\underline\theta)$) identical to the ones from $f_{\underline\theta}(\underline x)$?   * What is the general strategy for finding $h^*(\centerdot)$?


title: How to fit an ARIMAX-model with R?

text: I have four different time series of hourly measurements:   1. The heat consumption inside a house   2. The temperature outside the house   3. The solar radiation   4. The wind speed I want to be able to predict the heat consumption inside the house. There is a clear seasonal trend, both on a yearly basis, and on a daily basis. Since there is a clear correlation between the different series, I want to fit them using an ARIMAX-model. This can be done in R, using the function arimax from the package TSA. I tried to read the documentation on this function, and to read up on transfer functions, but so far, my code:               regParams = ts.union(ts(dayy))     transferParams = ts.union(ts(temp))     model10 = arimax(heat,order=c(2,1,1),seasonal=list(order=c(0,1,1),period=24),xreg=regParams,xtransf=transferParams,transfer=list(c(1,1))     pred10 = predict(model10, newxreg=regParams)      gives me: ![enter image description here](http://i.stack.imgur.com/UuJGj.png) where the black line is the actual measured data, and the green line is my fitted model in comparison. Not only is it not a good model, but clearly something is wrong. I will admit that my knowledge of ARIMAX-models and transfer functions is limited. In the function arimax(), (as far as I have understood), xtransf is the exogenous time series which I want to use (using transfer functions) to predict my main time series. But what is the difference between xreg and xtransf really? More generally, what have I done wrong? I would like to be able to get a better fit than the one achieved from lm(heat ~ temp*radi*wind*time). **Edits:** Based on some of the comments, I removed transfer, and added xreg instead:               regParams = ts.union(ts(dayy), ts(temp), ts(time))     model10 = arimax(heat,order=c(2,1,1),seasonal=list(order=c(0,1,1),period=24),xreg=regParams)      where dayy is the "number day of the year", and time is the hour of the day. Temp is again the temperature outside. This gives me the following result: ![enter image description here](http://i.stack.imgur.com/GDuzQ.png) which is better, but not nearly what I expected to see.


title: Analysis for cross-subjects correlation

text: Suppose that an investigator has collected some longitudinal data from each of the $n$ subjects, and that a correlation coefficient is computed on the data between any pair of subjects, leading to totally $n(n-1)/2$ correlation coefficients. After converting those correlation coefficients into $z$-values through Fisher's transformation, we can denote the $z$-values in the following matrix, $$ Z=\begin{bmatrix} * \ z_{12} \ z_{13} \ ... \ z_{1n}\\\ z_{12} \ * \ z_{23} \ ... \ z_{2n} \\\ z_{13} \ z_{23} \ * \ ... \ z_{3n} \\\ ... \\\ z_{1n} \ z_{2n} \ z_{3n} \ ... \ * \end{bmatrix}$$ The statistical goal here is to test the null hypothesis about the cross- subjects correlation at the population level. In the literature I have seen the following approach: use those $n(n-1)/2$ $z$-values (upper or lower triangle in the above matrix $Z$), and run a one-sample $t$-test. This seems troubling to me because those $n(n-1)/2$ $z$-values are not fully independent from each other (some of them are from the same subjects). In other words, the test would be inflated because the associated degrees of freedom are assigned with $n(n-1)/2-1$ even though there are only $n$ subjects. Another approach is to obtain the average $z$-value on each row in $Z$. That is, $i$th subject is represented with an average $z$-value $Z_i$ ($i=1,2,...,n$): $$ \begin{bmatrix} * \ z_{12} \ z_{13} \ ... \ z_{1n}\\\ z_{12} \ * \ z_{23} \ ... \ z_{2n} \\\ z_{13} \ z_{23} \ * \ ... \ z_{3n} \\\ ... \\\ z_{1n} \ z_{2n} \ z_{3n} \ ... \ * \end{bmatrix} \begin{matrix} Z_1 \\\ Z_2 \\\ Z_3 \\\ ... \\\ Z_n \end{matrix} $$ Then we can run a one-sample $t$-test on $Z_1, Z_2, ...Z_n$. The problem with this approach is that $Z_1, Z_2, ...Z_n$ are not independent with each other. Yet another alternative is to obtain the average $z$-value on each row in $Z$, but only in the lower or upper triangle to avoid the dependence issue. That is, the $i$th subject is represented with an average $z$-value $Y_i$ ($i=2,3...,n$) as below: $$ \begin{bmatrix} * \ - \ - \ ... \ -\\\ z_{12} \ * \ - \ ... \ - \\\ z_{13} \ z_{23} \ * \ ... \ - \\\\... \\\ z_{1n} \ z_{2n} \ z_{3n} \ ... \ * \end{bmatrix} \begin{matrix} - \\\ Y_2 \\\ Y_3 \\\ ... \\\ Y_n \end{matrix} $$ Then a one-sample $t$-test can be performed on $Y_2, Y_3, ...Y_n$. My intuition is that $Y_2, Y_3, ...Y_n$ are independent with each other, but the problem with this approach is that I'm not so sure whether it is a legitimate test for the original hypothesis about the cross-subjects correlations at the population level. Any suggestion or comment? Thanks!


title: How to check if data needs to be seasonally differenced in MATLAB

text: I am forecasting data using ARIMA. I would like to know if there is any test to check if seasonal differencing is needed in the ARIMA model. I know R uses Canova Hansen test, but does MATLAB provide anything similar to that? Also, is there an easy way of predicting if seasonally differenced model would give better results than the one without seasonal differencing? I would like to take the decision before fitting the model and estimating the data. Currently I am using BIC to check the performance and evaluate the best model. But the process is really slow. So I would like to know beforehand if i should seasonally difference the data or not.


title: Bounded expectation implied bounded conditional or vice versa?

text: If $\mathrm{E}\left(X\right)<\infty$ does that imply $\mathrm{E}\left(X|Y\right)<\infty$? How about vice versa? I'm thinking if we condition on an event (say $Y>2$) then if we have $\mathrm{E}\left(X\right)<\infty$ we can somehow use the theorem of total expectation to bound the other terms. Any thoughts? I've been searching for some kind of list of useful rules about expectations so please let me know if you have such a thing. Thank you! **Edit:** After reading what fellow users have said I believe I should ask whether we can bound expectations given events (for example $E(X\mid Y>2)$) if we know that $\mathrm{E}\left(X\right)<\infty$.


title: Why is an L1-estimator a special case of an M-estimator?

text: An L1 estimator minimizes the absolute difference between $Y$ and $x\beta$, whereas an M-estimator minimizes some function of the residuals. But the objective function needs to be differentiable right? Then why is an L1-estimator a special case of M-estimator?


title: Convert Cross Correlation to Probability value

text: How can I convert the [-1;1] bound cross correlation value to a probability [0;1]. Cross correlation measures the similarity of two signals / images A,B, I aim at expressing this in the probability that signal A is equal to signal B. Essentially, cross correlation has been used in the field of bio-image processing as a similarity measure but I was wondering how cross correlation could be transformed into a probability so that I could use various probability based methods for analysis. Is applying the square to the CC valid? Somebody have some references on that?


title: EM algorithm decreases!

text: I have used the Bayes Net Toolbox to build a small network, which consists of 3 nodes and is shown below. ![enter image description here](http://i.stack.imgur.com/KkC1t.png) Node 1 is a Bernoulli random variable, node 2 is a Gaussian random variable and node 3 is a softmax random variable with 3 possible values. The data is incomplete, so I use the EM algorithm to estimate the parameters. But the log- likelihood decreases from the beginning and stops with the last two iterations being equal. Does anyone know the possible solution to this situation?


title: Effect size interpretation

text: If r between X and Y is 0.74, then r-squared is 55%. According to many 'rules of thumb' and others, the effect size (in this case r) is considered 'moderate' or 'medium'. As such, only 55% of the variance in Y can be explained by X. The remaining 45% remain unexplained. **Question 1: Given the above, how good is r when it can only explain 55% of the variance?** I know this depends on the context of the study but 45% of unexplained variance, by any measure, is huge. I am afraid this is the number (not 55%) that my (non-statistical) target audience will remember at the end of the day. In other words, I may get a question such as "what is the point of the research when 45% of the result remain unexplained?" **Question 2: How should I answer this question?** Notwithstanding the above, I think 0.74 is a good and strong result!


title: Degrees of freedom for Chi-squared test

text: I am facing the following dilemma. I am aware of how to handle the one-sided Chi-distribution, but I am falling victim to how to handle degrees of freedom. Let me clarify with an example what I mean. I have the following obseverd and expected values:               [Observed Data]          #Periods      CountryI   CountryII     CountryIII     #(1900-1950)     100      150            20     #(1951-2000)     59       160            50          [Expected DATA]          #Periods   country I     Country II       CountryIII             #(1900-1950)  118.4         52                40     #(1951-2000)   80.5         90                25      My question is: Since this is a one sided-Chi square test, are the degrees of freedom counted by the formula: (columns-1)(rows-1), in which case I would have $(6-1)(2-1) = 5$? Or is that really just `country1 country2 country3` that matters, so that d.f. would be 3-1=2? Because d.f. is usually defined as the terms for the chi squared = 6, where we usually subtract 1 from it. Please help me out with this one.


title: Unable to understand fuzzy logic example

text: I shall be grateful for an explanation of how to apply fuzzy logic based reasoning for the following example: I have 3 crisp sets of features = color, shape and texture. The objects are red ball, red box, green cup and blue cup. Rule can be: R1: If color is red and shape is circular and texture is shiny then object is red ball. R2 : If color is red and shape is square and texture is coarse then object is red box. When the rule has AND then I should use the min operator, when there is an OR I should apply max operator. What I don't understand is how do I fuzzify the features and then make a decision that it is red ball or red box? Say the crisp value of color (intensity) is 0.45, crisp value of shape is 0.2 and crisp value of texture is 0.6. These values should pass through R1 and R2, so when the crisp numerical input pass through the membership set of red, circular and shiny and coarse how do I say that the input color is red and not something else (how to decide)? Same for the other sets? Do we assume that the membership value of color input is HIGHEST (max operator) in the Red fuzzy set and so we assume that the color is RED? I should have membership function for each of the features color, shape and texture. Will color have a membership function for red,green and blue set? Same for shape (circular, square) and texture (shiny and coarse). I have followed the examples in this tutorial but I am unable to understand how to apply for the case of color, shape and texture. Thank you. Explanation with any other example will also be highly appreciated.


title: How to describe study design looking at effect of a training program on two age groups  and how do I calculate required sample size?

text: I am conducting an experiment on comparing the training effects on young and older adults (meaning it will be 2 independent groups). This also means it will be a repeated measure (within groups) because I will test the score before and after training. I think I should use mixed 2 way ANOVA. (Am I right?) However I will also record the score before and after training in 2 different conditions. Is there a test to incorporate this variable as well? And how would I estimate the sample size if I already have the mean and standard deviation from a similar study using that design? Thank you in advance.


title: How to analyze a 34 item standardised survey purportedly measuring 6 dimensions using SPSS?

text: I am currently analyzing a 34-item standardized survey on Physics expectations. I gathered data from 104 female respondents from a certain school. The survey's supporting literature said that the survey measures 6 specific dimensions, with each dimension corresponding to 5-6 items. My professor asked me to conduct SPSS analysis on the data aside from the usual mean and standard deviation. Now, I am in the dilemma of identifying the correct statistical analysis to apply using the SPSS. Can anyone here please help me with this? What analyses can I apply? This is a project for a Statistics class, and we have no specific questions to answer. We are just asked to make meaning from the gathered data out of various statistical methods. Thank you very much for the response.


title: What is the advantage of transforming variables from nominal to ordinal/numerical when it reduces variance explained in CatPCA?

text: ### Context I have a dataset of 8 categorical variables. And I want to apply Categorical Principal Component Analysis (CatPCA). Before doing that, I have been advised to look at the transformation plots of all these variables after transforming them into nominal variables. This shows that some variables should be transformed into ordinal variables (when non- decreasing) and some into numeric (linear trend). Now when doing the CatPCA with all nominal variables and comparing this to the CatPCA with newly transformed variables, there is a slight decrease in variance explained in the dependent variable. ### Question   * So if variance explained decreases after the transformations, what is the advantage of the transformation?


title: How to find duplicate cases in a large dataset?

text: I have a large dataset which is made up of responses from 600 shoppers on 15 variables. When I tally the responses, it comes to 603. I know there were 600 respondents, so I must have recorded the responses of three people twice. Instead of checking each questionnaire against the dataset, can I run a test in Excel that would point out the double entries.


title: Computing Cohen's Kappa variance (and standard errors)

text: The Kappa ($\kappa$) statistic was introduced in 1960 by Cohen [1] to measure agreement between two raters. Its variance, however, had been a source of contradictions for quite a some time. My question is about which is the best variance calculation to be used with large samples. I am a inclined to believe the one tested and verified by Fleiss [2] would be the right choice, but this does not seems to be the only published one which seems to be correct (and used thoroughly fairly recent literature). Right now I have two concrete ways to compute its asymptotic large sample variance:   * The corrected method published by Fleiss, Cohen and Everitt [2];   * The delta method which can be found in the book by Colgaton, 2009 [4] (page 106). To illustrate some of this confusion, here is a quote by Fleiss, Cohen and Everitt [2], emphasis mine: > Many human endeavors have been cursed with repeated failures before final > success is achieved. The scaling of Mount Everest is one example. The > discovery of the Northwest Passage is a second. **The derivation of a > correct standard error for kappa is a third**. So, here is a small summary of what happened:   * 1960: Cohen publishes his paper "A coefficient of agreement for nominal scales" [1] introducing his chance-corrected measure of agreement between two raters called $\kappa$. However, he publishes incorrect formulas for the variance calculations.   * 1968: Everitt attempts to correct them, but his formulas were incorrect as well.   * 1969: Fleiss, Cohen and Everitt publish the correct formulas in the paper "Large Sample Standard Errors Of Kappa and Weighted Kappa" [2].   * 1971: Fleiss publishes another $\kappa$ statistic (but a different one) under the same name, with incorrect formulas for the variances.   * 1979: Fleiss Nee and Landis publish the corrected formulas for Fleiss' $\kappa$. At first, consider the following notation. This notation implies the summation operator should be applied to all elements in the dimension over which the dot is placed: $\ \ \ p_{i.} = \displaystyle\sum_{j=1}^{k} p_{ij}$ $\ \ \ p_{.j} = \displaystyle\sum_{i=1}^{k} p_{ij}$ Now, one can compute Kappa as: $\ \ \ \hat\kappa = \displaystyle\frac{p_o-p_c}{1-p_e}$ In which $\ \ \ p_o = \displaystyle\sum_{i=1}^{k} p_{ii}$ is the observed agreement, and $\ \ \ p_c = \displaystyle\sum_{i=1}^{k} p_{i.} p_{.i}$ is the chance agreement. So far, the correct variance calculation for Cohen's $\kappa$ is given by: $\ \ \ \newcommand{\var}{\mathrm{var}}\widehat{\var}(\hat{\kappa}) = \frac{1}{N(1-p_c)^4} \\{ \displaystyle\sum_{i=1}^{k} p_{ii}[(1-p_o) - (p_{.i} + p_{i.})(1-p_o)]^2 \\\ \ \ \ + (1-p_o)^2 \displaystyle\sum_{i=1}^{k} \displaystyle\sum_{j=1 \atop i\not=j}^{k} p_{ij} (p_{.i} + p_{j.})^2 - (p_op_c-2p_c+p_o)^2 \\} $ and under the null hypothesis it is given by: $\ \ \ \widehat{\var}(\hat{\kappa}) = \frac{1}{N(1-p_c)^2} \\{ \displaystyle\sum_{i=1}^{k} p_{.i}p_{i.} [1- (p_{.i} + p_{i.})^2] + \displaystyle\sum_{i=1}^{k} \displaystyle\sum_{j=1, i\not=j}^{k} p_{.i}p_{j.}(p_{.i} + p_{j.})^2 - p_c^2 \\} $ Congalton's method seems to be based on the delta method for obtaining variances (Agresti, 1990; Agresti, 2002); however I am not sure on what the delta method is or why it has to be used. The $\kappa$ variance, under this method, is given by: $\ \ \ \widehat{\var}(\hat{\kappa}) = \frac{1}{n} \\{ \frac{\theta_1 (1-\theta_1)}{(1-\theta_2)^2} + \frac{2(1-\theta_1)(2\theta_1\theta_2-\theta_3)}{(1-\theta_2)^3} + \frac{(1-\theta_1)^2(\theta_4-4\theta_2^2)}{(1-\theta_2)^4} \\} $ in which $\ \ \ \theta_1 = \frac{1}{n} \displaystyle\sum_{i=1}^{k} n_{ii}$ $\ \ \ \theta_2 = \frac{1}{n^2} \displaystyle\sum_{i=1}^{k} n_{i+}n_{+i}$ $\ \ \ \theta_3 = \frac{1}{n^2} \displaystyle\sum_{i=1}^{k} n_{ii}(n_{i+} + n_{+i})$ $\ \ \ \theta_4 = \frac{1}{n^3} \displaystyle\sum_{i=1}^{k} \displaystyle\sum_{j=1}^{k} n_{ij}(n_{j+} + n_{+i})^2$ (Congalton uses a $+$ subscript rather than a $.$, but it seems to mean the same thing. In addition, I am supposing that $n_{ij}$ should be a counting matrix, i.e. the confusion matrix before being divided by the number of samples as related by the formula $p_{ij} = \frac{n_{ij}}{\mathrm{samples}}$) Another weird part is that Colgaton's book seems to refer the original paper by Cohen, but does not seems to cite the corrections to the Kappa variance published by Fleiss et al, not until he goes on to discuss weighted Kappa. Perhaps his first publication was written when the true formula for kappa was still lost in confusion? Is somebody able to explain why those differences? Or why would someone use the delta method variance instead of the corrected version by Fleiss? [1]: Fleiss, Joseph L.; Cohen, Jacob; Everitt, B. S.; Large sample standard errors of kappa and weighted kappa. Psychological Bulletin, Vol 72(5), Nov 1969, 323-327. doi: 10.1037/h0028106 [2]: Cohen, Jacob (1960). A coefficient of agreement for nominal scales. Educational and Psychological Measurement 20 (1): 37–46. DOI:10.1177/001316446002000104. [3]: Alan Agresti, Categorical Data Analysis, 2nd edition. John Wiley and Sons, 2002. [4]: Russell G. Congalton and Green, K.; Assessing the Accuracy of Remotely Sensed Data: Principles and Practices, 2nd edition. 2009.


title: How to test the correlation between two variables in 50+ data subsets for statistical significance?

text: I have performed 61 Pearson correlation analyses between Variable A (buying ice cream) and Variable B (buying yoghurt) in my dataset (n=550). Each correlation is based on a different subset of the data (e.g., separate correlations for males and females; separate correlations for different age groups). Fortunately (or unfortunately!), for all 61 correlations, my p value is less than 0.0007431 (which is the correction of 0.05 using Sidax). I interpret this to mean that all the findings are significant, which in my cases is that there is a relationship between Variable A (people who buy ice cream) and Variable B (people who buy yoghurt). ### Question   * Is it normal to get all significant correlations when analysing the same correlation across a large number of subsets of data?   * Is it an acceptable results (i.e., "is not questionable" etc.)?  I can assure that my procedures for data collection and analysis are robust.


title: Best practices for measuring and avoiding overfitting?

text: I am developing automated trading systems for the stock market. The big challenge has been overfitting. Can your recommend some resources describing methods for measuring and avoiding overfitting? I started with training/validation sets, but the validation set always gets tainted. Also, the time series data is always changing because the market is always changing. How do you measure this and determine the likelihood of consistent results on unseen data? Thanks.


title: How to choose between nested or split-plot analysis of variance?

text: I am having serious problems understanding the type of design I have, I would appreciate any help you can give me to understand it: I have 324 data points, and four factors:   * `epoch` factor, with three levels.   * altitudinal `floor` factor, with three levels.   * `site` factor, with nine levels.   * `time` factor, with three levels. I have four points at each site, I have nine sites, these nine sites are divided in three, three sites correspond to the highest altitudinal floor, three to the middle altitudinal floor and three to the lower altitudinal floor, that give me a total of 36 points, in each point I took data three different times of day, giving me a total of 108 data points. And I repeated these procedure (the exactly same 36 points) in three epochs of the year. In other words: Each epoch have 108 data points, these 108 are divided in three (the levels of altitudinal floor factor), 36 to each altitudinal floor, these 36 are divided in three (three levels of the site factor), 12 to each site, and these 12 are divided in three (the levels of time factor) each time have 4 data points. I want know which anova is appropriate for test if there are differences at each factor: a nested anova:               aov(y ~ (epoch/floor/site/time), data=data)      or a split-plot anova:               aov(y ~ floor*site*time + Error(epoch/floor/site/time), data=data)      From what I've read, I think that is a nested design, but I dont understand why it is not a split-plot. Can anyone help me?


title: Modelling "gaps" in a Poisson regression

text: I am conducting an analysis into offending in relation to a time-dependent exposure - onset of heroin use, over a persons life. I am using poisson GEE with offending aggregate to a persons age-year. Many people have experienced time in prison, which means they can't experience the outcome, but are still experiencing the exposure. Including the reduced time (year minus proportion of year spent in prison) as an offset does not seem to be appropriate, as the exposure still applies. Is there any other modelling approaches you can recommend? Many thanks


title: Why are most standard goodness of fit tests based only on continuous distributions?

text: I tried to search info regarding this fact but I don't really understand why most of the standard goodness of fit tests (e.g. Kolmogorov-Smirnov, Anderson- Darling, a part of the Chi-square test, perhaps!) work only with continuous distributions. Can someone help me? Thank You.


title: What test is appropriate for this binned data?

text: I have collected questions from 10th grade students about genetic engineering before and after learning about the subject and classified the questions based on their cognitive level.. Before learning about the subject, I collected a total of 45 student questions. I classified 28 of those as "input" (basic) questions, 8 questions as "processing" (intermediate) and 9 as "output" (advanced). After learning about the subject, I collected from the same students 28 questions, distributed as follows: 11 "input", 7 "processing" and 10 "output". What test is appropriate to see if the distribution of student questions among the three categories is different? I tried conducting a prop.test but that only distinguishes between two categories (success and failure). Thanks.


title: Modelling both mean and dispersion of count data

text: I have a model of the following form: $P(Y \mid X) = \,D(\mu,\sigma^2) ~~\text{where}$ $\mu = f(X) ~~\text{and}~~ \sigma^2=g(X)$ where $y$ is the response vector of count data, $X$ is the predictor matrix, $f()$ and $g()$ are linear inverse link functions (I don't mind what they are), and $D$ is some probability distribution (maybe Poisson with overdispersion). Does this model have a name? How would you estimate it? Is there a package in R that can do it? I understand that usually one would use negative binomial or Poisson GLM, but here I have heteroscedasticity that depends on $X$ not in the same way as the expected value of the response (that is $f()$ and $g()$ are not the same). Thank you for your help!


title: How to check whether maximum likelihood estimation optimizer has converged in R?

text: I got AIC values of all models to identify the best model using R language. As I heard, best model produce the smallest AIC value, but maximum likelihood estimation procedure optimizer should converge. How can I check whether maximum likelihood estimation procedure optimizer has converged or not in R language?


title: ezANOVA package in R: right method for my experimental design?

text: I have a problem with the `ezANOVA` package using for rm ANOVA or a problem with my statistic design of the experiment (or both ;-)) My file looks like this,               dput(head(reg,20))          structure(list(sample = c(13L, 14L, 15L, 13L, 14L, 15L, 13L,      14L, 15L, 13L, 14L, 15L, 13L, 14L, 15L, 13L, 14L, 15L, 1L, 2L     ), days = c(7L, 7L, 7L, 14L, 14L, 14L, 21L, 21L, 21L, 28L, 28L,      28L, 35L, 35L, 35L, 42L, 42L, 42L, 7L, 7L), treat = structure(c(1L,      1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,      1L, 2L, 2L), .Label = c("asynchron", "control", "synchron"), class = "factor"),      repl = c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L,      1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L), disp = structure(c(2L, 2L,      2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,      2L, 2L, 2L), .Label = c("high", "low", "medium"), class = "factor"),      gbiov = c(205088809L, 236043290L, 182258079L, 356710654L,      519036032L, 423556500L, 265916788L, 260495312L, 424304021L,      456717288L, 412530454L, 495216491L, 360012083L, 246070292L,      351623075L, 99075488L, 141565279L, 109303015L, 157719021L,      275901505L), rich = c(15L, 15L, 15L, 12L, 12L, 10L, 10L,      10L, 10L, 8L, 7L, 9L, 9L, 9L, 8L, 6L, 8L, 8L, 15L, 15L),      evenness = c(0.810893754, 0.766130198, 0.82401103, 0.422734873,      0.365797654, 0.406629574, 0.317890623, 0.434262493, 0.475385016,      0.311819454, 0.307131166, 0.321269955, 0.365280751, 0.505205296,      0.43643584, 0.484220445, 0.455471176, 0.423730082, 0.811921742,      0.833306048)), .Names = c("sample", "days", "treat", "repl",      "disp", "gbiov", "rich", "evenness"), row.names = c(NA, 20L), class = "data.frame")      Again, my dataframe looks like (column sorted)               sample   days  treat    repl   disp   gbiov richness evenness      I have measured totalbiovolume(gbiov), richness and evenness over six weeks with a weekly sampling interval. I have 2x2 factorail design, where I have three treatments (treat): Synchron, asynchron and control and three different dispersal rate (disp) low, medium and high. So I decided to use a rm ANOVA by using `ezANOVA`, because I took a sample every week from the same erlenmeyer flask. My within SSw is (treatment and disp) and between (days?). I used repl. column as the specifying case/Ss identifier for the replicates of each combination ($2\times 2$ design). **My input:**               test<-ezANOVA(data=reg, dv=.(gbiov), wid=.(repl), within=.(treat, disp, days))      **Error message:**               Warnung: Converting "repl" to factor for ANOVA.     Warnung: "days" will be treated as numeric.     Fehler in ezANOVA_main(data = data, dv = dv, wid = wid, within = within,  :      One or more cells is missing data. Try using ezDesign() to check your data.      My problem in my dataframe is, that for the combination (asynchron $\times$ low) I have only three replicates (one is missing during the experiment periode), I got the advice to use `ezDesign` to check my dataframe and after I did so, I add the missing replicates by reseting the missing value to zero for each time point, but still I get the same error message, except when I add days to the within-SSw and omit between-SSb. But I don´t know if that is correct? Perhaps, the SS identifier makes no sence?! Or my design is wrong and I have no within and between conditions Or it is a problem of one missing replicate Thanks in advance for solving the problem Best, Nils


title: Compare R-squared from two different Random Forest models

text: I'm using the randomForest package in R to develop a random forest model to try to explain a continuous outcome in a "wide" dataset with more predictors than samples. Specifically, I'm fitting one RF model allowing the procedure to select from a set of ~75 predictor variables that I think are important. I'm testing how well that model predicts the actual outcome for a reserved testing set, using the approach posted here previously, namely, ![](http://mathurl.com/435zb2a.png) ... or in R:               1 - sum((y-predicted)^2)/sum((y-mean(y))^2)      But now I have an additional ~25 predictor variables that I can add. When using the set of ~100 predictors, the R² is higher. I want to test this statistically, in other words, when using the set of ~100 predictors, does the model test _significantly_ better in testing data than the model fit using ~75 predictors. I.e., is the R² from testing the RF model fit on the full dataset significantly higher than the R² from testing the RF model on the reduced dataset. This is important for me to test, because this is pilot data, and getting those extra 25 predictors was expensive, and I need to know whether I should pay to measure those predictors in a larger follow-up study. I'm trying to think of some kind of resampling/permutation approach but nothing comes to mind. Any help is appreciated.


title: Correct use of partial derivatives? (Example: polynomial regression)

text: [update] It seems I had my question's original title misfocused a bit; concerning the use of partial derivatives I found some explanation/confirmation in wikipedia's partial derivative There it is stated -and explained with a 3-d-picture of a surface-, that indeed the three partial derivatives suffice to find the/a common extremum. In the occurrence of my problem I seem to have had a misleading intuition, based on the observation of non-independence of the x and x-square-data and felt thus cautious (too much, apparently).   Full exposition of my proceeding/reasoning and source of question is here _I adapted the subject of my question; if this seems inappropriate please feel free to roll it back. The original question follws below_ [end update] * * * I'm just trying to reinvent the wheel - I want to understand (and implement for an example) the computation of the polynomial regression. For the simplest nonlinear approach let's use the estimated model $$\hat y_k=a+b x_k + c x_k^2 \qquad \text{(for k=1 to n)}$$ with the minimizing criterion $$\sum (y_k - \hat y_k)^2 =\sum (y_k - (a+bx_k+cx_k^2))^2= \text{Min}$$ I have expressed all the occurrng and needed variances and covariances in a matrix, cofactored with $a,a^2,b,b^2,c \text{ and }c^2$ . To find the minimum I must use the derivatives; however I'm not much familiar with partial derivatives. So **my question is** : do I get the global minimum if I solve my matrix- expression setting the partial derivatives resp to $a,b,c$ _separately_ (and find the then possible common solution)? Or do I need to involve some "joint derivative" or "common derivative" where, for instance something like $\partial ab$, thus _a_ and _b_ , are simultanously used in one derivative?


title: A construct as an independent variable?

text: I know that you can have a construct as a dependent variable: MANOVA or MANCOVA. Is it possible for me to have a construct as an independent variable? I have a construct that consists of 3 variables that I want to use as my independent variable and I have one variable that I want to use as a dependent variable.


title: Excel pivot table with multiple total columns

text: Is there a way in Excel to create a pivot table with multiple column totals? Here's column names to demonstrate what I'm trying to do: Jan | Feb | Mar | Apr | May | Jun | **6 Month total** | Jul | Aug | Sep | Nov | Dec | **6 Month total** | **12 Month total** The **12 Month total** column is of course the regular Grand Total. The two 6 Month Total columns are the ones I'm looking to create.


title: Confidence interval of slope in linear regression

text: When computing a confidence interval of slope in linear regression, should you use the z- or t-statistic?


title: ANOVA in lmerTest for in-between study with repeated measurements

text: I have a bit of trouble understanding how to formulate a correct formula in R using the lmerTest-library. Let’s assume there is an 2x2 in-between design with repeated measurements. Those repeated measurements can also be categorized in two ways: the position and the difficulty. Therefrom I got the following data table:               proband;between1;between2;position;difficulty;score     1      ;0       ;0       ;0       ;0         ;12     1      ;0       ;0       ;1       ;0         ;13     1      ;0       ;0       ;2       ;1         ;16     1      ;0       ;0       ;3       ;1         ;12     2      ;1       ;0       ;0       ;1         ;11     2      ;1       ;0       ;1       ;0         ;14     2      ;1       ;0       ;2       ;0         ;12     2      ;1       ;0       ;3       ;1         ;16     ...      between1 can be 0 and 1   between2 can be 0 and 1   position can be 0, 1, 2 and 3   difficulty can be 0 and 1 Is it correct, if I now use the following R script to analyse if there are any main or interaction effects?               zf.data=read.csv(file="data.csv", sep=";")     attach(zf.data)     library(lmerTest)     model = lmerTest::lmer(score ~ between1*between2*position*difficulty + (1 | proband), data=zf.data)     anova(model)      I am quite new to statistics, in case this question is too simple I will withdraw it and be deeply sorry.


title: Combinatorial Probability

text: Let there be an exam consisting of **q questions** each of which having **a alternatives** where there's only one correct answer to each question. Assumptions:   * The questions are independent; that is, the responses do not interfere with the responses to other questions.   * The probability of each alternative is equal. The probability of the simultaneous occurrence of **k identical wrong answers (IWA)** to a single question is as follows: $$P(k) = \left( \frac{1}{a^2}\right)^k$$ I need to adopt this to a situation where there are **q questions** with **a alternatives** under those assumptions. And to find $P(\geq k)$, actually.


title: Computation of a binomial probability: nan result

text: I am computing binomial probabilities for many parameters. I am using the library scipy.stats:               b = scipy.stats.binom(n,p)     probability = b.pmf(k)      I have noticed that the result is `nan` when both n and p are 0. This raises errors in my program, so I replaced nan by 0, since the probability is 0. Is this a good replacement in your opinion?


title: Is this variable ordinal or interval ?

text: For my research paper, I want to measure headache frequency and I asked participants " Since last month how many times you have had headche ? and they gave answers like 2 times, 3, 5, 10, 14 times etc.. I am confused with the type of this data, Shall I treat it as ordinal or interval ? Thanks alot in advance


title: Is the bootstrap estimate of the mean biased when a single extreme value is in the sample?

text: My sample includes $n$ random observations, while $n-1$ of these observations are in the range (0-1) there is also one observation that gets very high value. For example, a sample of prices where $n-1$ prices are lower then 1 dollar and one price in the sample is 1000 dollars. **Is the bootstrap estimate of the mean biased when a single extreme value such as above is in the sample?** It seems to me that it should depend on $n$ and $B$, where $B$ is the number of sub-samples.


title: Unbiased estimation of the mean of a normal distribution when the parameter space is restricted

text: Let $X_1, . . . , X_n$ be an iid sample from $\text{N}(\theta,\sigma^2)$, in which $\sigma^2$ is known and $\theta$ can have the values $0,±1,±2, . . .$, and let $\bar{X}_n^*=\max\\{k\in\mathbb{Z}:k\leq \bar{X}_n\\}$ be the largest integer smaller or equal than $\bar{X}_n$.   1. Show that $\bar{X}_n^*$ is an unbiased estimator for the restricted parameter $\theta$. That is, prove that $\text{E}_\theta\left[\bar{X}_n^*\right]=\theta$, for every integer $\theta$.   2. There exist positive constants $a$ and $b$ such that, for all sufficiently large $n$, $\text{Var}_\theta\left[\bar{X}_n^*\right] \leq a\cdot e^{−b\,n}$, for every integer $θ$.


title: Detecting initial trend or outliers

text: In my test procedure I sequentially take 10 measurements of a recently perturbed physical system, and I often find the first few (between 0 and 4) measurements can be inaccurate because the system has not settled/equilibrated. The initial non-equilibrated measurements may just have less accuracy or show a distinct trend, while the remainder should not have a trend and just show random variation. I am currently eyeballing the data to discard any initial poor data but want to have a more rigourous and automated method for selecting the "good" data. I am aware (eg here) of various outlier tests (Chauvenet, Grubbs, Pierce, Generalized ESD) but because of "masking" and trending don't think they will give reliable results. The experimental conditions are fairly well controlled so there are unlikely to be other genuinely sporadic outliers appearing. Is there a more appropriate statistical test (than say Generalised ESD) for my procedure? _Update:_ To give some more details I have run about 15 procedures so far and while some datasets have all good values most need 2 or 3 measurements dropped, so it seems reasonable to query the first 4. The problem description above is somewhat simplified in that more than one number is determined for each measurement, ie there are several system parameters measured. Only 1 or 2 parameters show non-equilibrium effects though so there is limited scope for cross-checking. The aim is to get reliable values for parameters of an unchanging system and therefore the model I am trying to fit is a series of constant values with some measurement error. The measuring is time-consuming but automated and can currently run overnight, the analysis is performed afterwards. It is possible to increase the number of measurements but this obviously takes more time. From the previous results I have a good idea of the expected variation in measurements at equilibrium, so I could use that information in addition to the within-dataset variation to assess outliers. Currently I only use that information to determine if the final error is acceptable. Six good data-points is sufficient to get an adequate estimate of the parameters and to detect any obvious problems with the process, so just dropping the first 4 is an acceptable solution. It seems from the discussion this is the best solution, unless the Dixon's test variant can be used. Where can I read more about that?


title: Is matching only for treatment effects with selection bias?

text: Propensity score matching (and other matching techniques) are used, as far as I have seen, exclusively for identifying causal effects of a treatment (intervention) and particularly where there is a suspicion of selection bias. However I do not see why they might not be useful for other cases of comparing two groups within a sample. E.g. I am looking at differences between rural-to- urban migrants and urban natives in education outcomes. It makes little sense to see this as a 'treatment' effect (although it might if I was comparing rural-to-urban migrants with _rural_ natives), and causation is not exactly the issue. I want to understand whether there is a difference that can be explained by migration status separately from that explained by differences between migrants and urban natives in terms of wealth, parents' education and other observed variables. This would usually be done using linear regression and ordinary least squares. But wouldn't there be some advantages to using a matching technique here? Specifically, doesn't OLS apply the same linear relationships among covariates to the whole sample, while matching could better allow for variation in those relationships? In this type of case, would matching perhaps be comparable to regression with lots of interaction effects allowed among the covariates? (I am aware of this earlier question but everything I come across on matching seems to assume you are trying to evaluate causal effects of a treatment. It would be great to hear about any counter-examples to my claim that matching is only ever used for treatment effects, if there are any.)


title: Looking for significant variability of categorical data in columns of data frame

text: I have a data frame with each column contains categorical data: 0,0.5,1,-1                          T116 T111 T11 T109 T108 T107 T105 T104 T55  T52 T58  T1     cg21277995  1.0  0.5 0.5  0.5    1  0.5  1.0  1.0 1.0  0.5 1.0 0.5     cg01324372  0.0  0.5 0.0  0.5    1  1.0  0.5  1.0 1.0  0.5 1.0 0.5     cg13836638  1.0  0.5 0.0  0.5    1  1.0  0.5  1.0 1.0  1.0 1.0 0.5     cg00002719  0.5  0.0 0.0  0.5    1  0.5  0.5  1.0 1.0  0.5 1.0 0.5     cg15804767  0.5  0.5 0.5  0.5    1  1.0  0.5  1.0 0.5  0.5 1.0 0.5     cg17228900 -1.0  0.5 0.5  0.5   -1  0.5  0.5  0.5 0.5 -1.0 0.5 0.0      I want to get only the columns that have a significant variability among the others. Is there any test to calculate the variability in each column compared with the others? Probably permutation tests? I want to obtain a p-value for each column and then select the best ones Thanks in advance


title: Consequences from comparison of two models

text: I have time series data for two firms relating cost to quantity produced and I model the two firms' cost functions separately using OLS: Firm 1:〖cost〗_1t= β_1+β_2 〖qty〗_1t+ β_3 〖qty〗_1t^2+ β_4 〖qty〗_1t^3+ ε_1t Firm 2:〖cost〗_2t= β_1+β_2 〖qty〗_2t+ β_3 〖qty〗_2t^2+ β_4 〖qty〗_2t^3+ ε_2t It is believed the firms have the same regression coefficients but doubted the error variances are the same (and a Goldfeld-Quandt test shows the latter). What are the consequences for the coefficient estimates derived above if: (i) the 2 firms actually have different variances but the same coefficients and; (ii) the 2 firms have different coefficients but the same variance? Note, the two separate regressions for the firms _appear_ to show different coefficients.


title: Learning classifier system frameworks

text: For some time already I am studiyng and applying machine learning methods. Recently for some particular problem, where methods like SVM, RF, neural nets etc. failed I got limited success with genetic programming, which was quite new area for me in practice at the moment. But right now it looks like, I really need to implement something like Learning classifier system. So the question is, are there any existing frameworks for LCS implementation ? To be more precise I am looking for frameworks on any language from the list: Java, Python, R, Matlab. Thanks. P.S.: I am still not sure if this is a correct place to ask my question, should I go with it to CrossValidated instead?


title: Why doesn't the Cramer-Rao lower bound apply?

text: Let $X_1, X_2, \dots, X_n$ be a sample of i.i.d. random variables, with density $$f_\theta=\frac{2}{3\theta}\left(1-\frac{x}{3\theta}\right) $$ for $0 < x < 3\theta$. And $f_\theta=0$ if $ x < 0$ or $ x>3\theta$ Let $\hat{\theta}=\overline{X}$ be an estimate for $\theta$ I showed that $\hat\theta$ is an unbiased estimator for $\theta$ and it's a consistent estimator. My question is: Why doesn't the Cramer-Rao lower bound apply to unbiased estimates of for this distribution?


title: Hamiltonian Monte Carlo and discrete parameter spaces

text: I've just started building models in **stan** ; to build familiarity with the tool, I'm working through some of the exercises in Bayesian Data Analysis (2nd ed.). The Waterbuck exercise supposes that the data $n \sim \text{binomial}(N, \theta)$, with $(N, \theta)$ unknown. Since Hamiltonian Monte Carlo doesn't permit discrete parameters, I've declared $N$ as a real $\in [72, \infty)$ and coded a real-valued binomial distribution using the `lbeta` function. A histogram of the results looks virtually identical to what I found by computing the posterior density directly. However, I'm concerned that there may be some subtle reasons that I should not trust these results in general; since the real-valued inference on $N$ assigns positive probability to non- integer values, we know that these values are impossible, as fractional waterbuck don't exist in reality. On the other hand, the results appear to be fine, so the simplification would appear to have no effect on inference in this case. Are there any guiding principles or rules of thumb for modeling in this way, or is this method of "promoting" a discrete parameter to a real bad practice?


title: Predicting values for a Naive Bayes in R

text: I have the following data set and am trying to predict on the naive bayes values               > head(dat)        lead_id                      email owns_home credit_type buyer_returned     1 12073402              b@hotmail.com         0        GOOD            \\0     2 12073514                 mo@cox.net         0        GOOD            \\0     3 12073541                v@yahoo.com         0   EXCELLENT            \\0     4 12073591              rob@yahoo.com         1   EXCELLENT           \001     5 12073653                c@ymail.com         0   EXCELLENT            \\0     6 12073768                r@yahoo.com         0   EXCELLENT            \\0      So I run the following command to run NB and to predict on the desired values.               m <- naiveBayes(buyer_returned ~ owns_home + credit_type, dat)     m          > predict(m, dat[1:10, -5])     Error in FUN(1:2[[2L]], ...) : subscript out of bounds     > predict(m, dat[1:10, -5], type="raw")     Error in FUN(1:2[[2L]], ...) : subscript out of bounds          > predict(m, dat[1:10, 5])     Error in log(sapply(attribs, function(v) { :        Non-numeric argument to mathematical function     > predict(m, dat[1:10, 5], type="raw")     Error in log(sapply(attribs, function(v) { :        Non-numeric argument to mathematical function      I'm not sure why I get the subscript out of bounds error. What is the issue. Here's info that might be useful:               > str(dat)     'data.frame':   50000 obs. of  5 variables:      $ lead_id       : int  12073402 12073514 12073541 12073591 12073653 12073768 12073809 12073818 12073897 12073943 ...          $ email         : Factor w/ 49755 levels "0011werbicki1100@gmail.com",..: 6231 32482 47424 38883 6784 37311 33483 37525 41427 35459 ...      $ owns_home     : int  0 0 0 1 0 0 0 0 0 0 ...          $ credit_type   : Factor w/ 5 levels "\\N","EXCELLENT",..: 4 4 2 2 2 2 2 4 2 4 ...      $ buyer_returned: Factor w/ 2 levels "\001","\\0": 2 2 2 1 2 2 2 2 2 2 ...         > str(m)         List of 4          $ apriori: 'table' int [1:2(1d)] 2865 47135       ..- attr(*, "dimnames")=List of 1       .. ..$ Y: chr [1:2] "\001" "\\0"          $ tables :List of 2       ..$ owns_home  : num [1:2, 1:2] 0.358 0.268 0.479 0.443           .. ..- attr(*, "dimnames")=List of 2           .. .. ..$ Y        : chr [1:2] "\001" "\\0"       .. .. ..$ owns_home: NULL           ..$ credit_type: table [1:2, 1:5] 0.1888 0.086 0.3218 0.4358 0.0932 ...       .. ..- attr(*, "dimnames")=List of 2       .. .. ..$ Y          : chr [1:2] "\001" "\\0"           .. .. ..$ credit_type: chr [1:5] "\\N" "EXCELLENT" "FAIR" "GOOD" ...      $ levels : chr [1:2] "\001" "\\0"          $ call   : language naiveBayes.default(x = X, y = Y, laplace = laplace)      - attr(*, "class")= chr "naiveBayes"      I am trying to get the predicted values for each lead_id.


title: Distribution of a centered standardized sample

text: Foreword : this is **not** homework, but a real problem : in a Bayesian model comparison context, I am trying to work out the correct prior density of mixed-model parameters which are, for computational reasons, computed as the centered and standardized means of other "raw" parameters, which are themselves given prior distributions. See Gelman A, Hill J. _Data Analysis Using Regression and Multilevel/Hierarchical Models_. 1st ed. Cambridge University Press; 2006. This textbook uses this trick (or variants thereof) quite often, and various simulations have convinced me of its efficacy, at least in selected cases. I currently do not have on hand any reference material (textbook, etc...) for this question, whose answer is quite probably a standard result. Could someone lend me a hand (if only by pointing me to a Web-accessible reference, which my Google attempts have been unsuccessful at referring me to) ? Let $X=X_1,\ldots X_N$ be a sample (i. e. vector of $N$ _independent_ random variables sampled from the same) **normal** distribution ${\textrm N}(\mu_x, \sigma^2_x)$. Let $m={{\sum_{i=1}^n X_i} \over N}$ and $s^2={{\sum (X_i)^2} \over {N-1}}$ the unbiased estimators of the mean and variance of this distribution. What is the (joint) distribution of the _components_ $Y_i$ of $Y={{X-m} \over s}$ ? I am sorely tempted to answer that they are $N$ independent normal standard random variables (i. e. distributed as ${\textrm N}(0, 1)$), and my intuition is seconded by various simulation results, but I have trouble **_proving_** it. What I know (or prove to my satisfaction) can be summarized as follows :   * $m$ has a ${\textrm N}(\mu_x, {{\sigma^2_x} \over N})$ distribution ; therefore, the common mean of the $Y_i$ is 0 (trivial from the definition of the mean).    * The common variance of the $Y_i$ is 1 (a trifle less trivial but can be worked out from the definition).   * ${s^2} \over {\sigma^2_x}$ has a $\chi^2_{N-1}$ distribution (classical result that I remember having worked out).   * Any textbook will state that $m \over s$ has a $t_{N-1}(0, 1)$ distribution ( _Not_ trivial, but I remember having worked it out, the crucial point being the independence of $m$ and $s$. I think that this historically is the origin of the distinction of the $t$ distribution, BTW...).   * Obviously, $Y_i={{X_i} \over s}-{m \over s}$. However, neither $m$ nor $s$ are independent of the $X_i$. Therefore   * the distribution of ${X_i} \over s$ is unknown (to me), and   * the distribution of the sum is _not_ the convolution of the distribution of its terms. Where am I goofing ?


title: How would you do Bayesian ANOVA and regression in R?

text: I have a fairly simple dataset consisting of one independent variable, one dependent variable, and a categorical variable. I have plenty of experience running frequentist tests like `aov()` and `lm()`, but I cannot figure out how to perform their bayesian equivalents in R. I would like to run a bayesian linear regression on the first two variables and a bayesian analysis of variance using the categorical variable as the groupings, but I cannot find any simple examples on how to do this with R. Can someone provide a basic example for both? Additionally, what exactly are the output statistics created by bayesian analysis and what do they express? I am not very well-versed in stats, but the consensus seems to be that using basic tests with p-values is now thought to be somewhat misguided, and I am trying to keep up. Regards.


title: Distribution of ratio of sample means from two independent normal variables?

text: ## The Question We have a sample of size $N$ with mean $\bar{x}$ and SD $\bar{\sigma_x}$ from a random variable $X \sim \mathcal{N} (\mu, \sigma^2)$ We have a sample of size $M$ with mean $\bar{y}$ and SD $\bar{\sigma_y}$ from a random variable $Y \sim \mathcal{N} (c\mu, c^2\sigma^2)$ We wish to find estimates of $\mu$ and $c$, along with the distributions of those estimates ## My Progress So Far   * $\bar{x}$ is one obvious estimate of $\mu$, and we know it has a t-distribution, but using only this statistic ignores the information about $\mu$ contained in $\bar{y}$   * $\frac{\bar{y}}{\bar{x}}$ would give us an estimate of $c$, but what distribution would it have? $\frac{Y}{X}$ has a Cauchy distribution, but what is the analogous distribution when using the ratio of sample means? Said another way: Normal : t-distribution :: Cauchy : ???   * Once we have an estimate of $c$, we could divide $\bar{y}$ by that estimate to get another estimate of $\mu$, and thereby extract the additional info about $\mu$ contained in $\bar{y}$. But what is the distribution of that estimate, and how to combine it with our $\bar{x}$ estimate? Things seem to be getting confusing since now we'd have a t-distribution random variable divided by our sample-equivalent-of-a-Cauchy random variable.... Is there a more straightforward way to do this?


title: ARIMA(0,2,2) model fitting

text: My data is below.After taking differences 2 times,it looks stationary. Only lag 2 is significant in both pacf and acf plot.EACF says it is ARIMA(0,2,2) and I think it can be. However,if it is, aren't there supposed to be tailing off in PACF and significant both 1 and 2 lag. I put some plots and model fitting result below. I am waiting for your advices.Thanks is adavance               [1]    78.148736    82.362301    85.723792    89.179080    93.161329    98.285106       [7]   102.859345   105.121553   109.301252   114.655529   119.834085   123.862782      [13]   127.814803   130.247296   131.758083   134.446371   137.915851   138.865710      [19]   139.498699   140.453177   142.105410   143.506534   142.464967   138.959921      [25]   135.560157   132.511408   128.042383   123.033841   117.357823   111.048266      [31]   104.973723    99.239756    93.275708    90.735442    86.725773    81.298565      [37]    76.319429    72.158145    70.384700    68.749347    66.504599    62.839408      [43]    58.803546    55.084801    50.034986    45.158016    42.119715    38.714206      [49]    34.212532    28.673336    23.231000    18.985425    14.896459    10.089654      [55]     4.871620    -2.948507   -11.121487   -16.665217   -21.012469   -25.845888      [61]   -32.918252   -40.058004   -46.531557   -52.461981   -58.746629   -65.520068      [67]   -73.166780   -81.246150   -90.916097  -100.950862  -110.556203  -120.576848      [73]  -130.998991  -140.491158  -150.041804  -158.967854  -167.749404  -177.381233      [79]  -185.706326  -192.954778  -200.089766  -207.319277  -215.431529  -224.885963      [85]  -234.760472  -243.693291  -251.378333  -259.255769  -267.695742  -277.253981      [91]  -288.460978  -299.293036  -309.631924  -319.873367  -330.624362  -342.764577      [97]  -354.945638  -365.290975  -376.070030  -387.484814  -399.084482  -410.684334     [103]  -422.114317  -432.606447  -443.512131  -455.033873  -465.441883  -475.172921     [109]  -486.433274  -497.249294  -508.283388  -521.443677  -535.586560  -548.697855     [115]  -561.976309  -576.901046  -591.117945  -604.964464  -619.652964  -633.832902     [121]  -648.036244  -662.790856  -678.667738  -694.466336  -709.636897  -723.648894     [127]  -736.047369  -747.682971  -758.410251  -768.075198  -777.187845  -787.093211     [133]  -799.103186  -809.792309  -819.243285  -829.071917  -837.829408  -845.610798     [139]  -854.877214  -865.713652  -876.429800  -886.308464  -895.490517  -903.948036     [145]  -914.319994  -924.675138  -933.085928  -942.196491  -952.717216  -962.231467     [151]  -970.993981  -979.975584  -989.348858  -999.394030 -1009.710203 -1017.895786     [157] -1025.066235 -1032.983167 -1041.579235 -1051.556946 -1060.286981 -1069.899993     [163] -1081.721708 -1092.839751 -1101.202291 -1109.546476 -1118.355565 -1124.851188     [169] -1131.410932 -1138.830048 -1147.041922 -1156.586296 -1165.934825 -1177.124854     [175] -1188.691257 -1201.133894 -1214.156490 -1227.369868 -1241.058276 -1255.250397     [181] -1270.291862 -1286.125954 -1302.427868 -1318.572693 -1332.972092 -1347.901450     [187] -1363.525918 -1378.926688 -1393.967150 -1409.059311 -1421.825466 -1434.386226     [193] -1446.977634 -1458.660748 -1470.548458 -1482.422500 -1492.986396 -1504.392771     [199] -1515.468191 -1525.285415                    arima(x = mydata, order = c(0, 2, 2))          Coefficients:              ma1     ma2           0.1936  -0.327     s.e.  0.0652   0.061          sigma^2 estimated as 0.9952:  log likelihood = -280.63,  aic = 565.26               arima(x = mydata, order = c(0, 2, 2), fixed = c(0, NA))          Coefficients:           ma1      ma2             0  -0.3204     s.e.    0   0.0601          sigma^2 estimated as 1.041:  log likelihood = -285.03,  aic = 572.05      ![acf after twice diff](http://i.stack.imgur.com/fvzbl.png) ![pacf after twice diff](http://i.stack.imgur.com/RktKa.png)


title: Determining MLE of a distribution with sufficient statistics

text: Please pardon me for asking a simple question but I find this field rather very hard. Anyway, here it goes:   1. We are trying to find MLE for a distribution. Now for a given distribution, its MLE is a number that is difficult to estimate. So what we do then is to use Sufficient Statistics for that particular distribution to find MLE over that distribution. Am I correct here? I will be absolutely grateful for any help here.


title: When to use likelihood ratio test vs. incremental F-test in nested nonlinear model selection?

text: I have two nested nonlinear models and I want to know which provides a better fit to some data. I see descriptions of both the likelihood ratio test and of the incremental F-test (also called the nested F-test or the extra sum-of- squares F-test) but I cannot find any discussion of when you would use one or the other or their relative advantages/disadvantages. Can anyone enlighten me?


title: Mean Square Error definition for symmetric models

text: I want to estimate the parameters of the sum of two complex sinusoids embedded in Gaussian noise: $y[n]=a_1 {\rm e}^{\,j n \theta_1} + a_2 {\rm e}^{\,j n \theta_2} + w[n]$ where $w[n]\sim {\cal N}(0,\sigma^2)$ and $n=1,\ldots,N$. Say that I want to compare the performance of different estimators (Prony, Maximum Likelihood, etc.), so I run a Monte Carlo simulation. For each realization of the noise term $w[n]$ I will obtain the estimates $\hat{a}_1$,$\hat{a}_2$,$\hat{\theta}_1$,$\hat{\theta}_2$. The problem now is that if I want to compute, say, the Mean Square Error for the amplitude estimates $\mbox{MSE}(\hat{a}_1,\hat{a}_2)=E\left[(a_1-\hat{a}_1)^2+(a_2-\hat{a}_2)^2\right]$, I need to associate each of the values returned by the estimator to the correct "ground truth" parameter $a_1$ or $a_2$. If the noise variance is low and $|a_i-\hat{a}_i|$ is not too small this is not hard; otherwise it does not seem to be an easy task. Any suggestions?


title: High dimensional parametric skew distribution

text: What are the best researched options for a parametric multivariate probability distribution that is able to be skewed, given curvature, and possibly even multiple modes? I have done some research into unified skew normal distributions and unified skew elliptical distributions, but the available research is very lean, and I've yet to find any method for actually estimating such distributions given data which is purported to follow it.


title: Canonical correlation and prediction

text: Is it possible to predict values of variables in a matrix Y from values of variables in a matrix X, ie. using canonical correlations since some variables are correlated among these two matrix? X and Y contain variables for the same individuals. If possible, I would like an example using R. Edit: the aim is to fill missing values in the matrix Y from correlated variables that are in matrix X.


title: Deviation from Gaussian smoothed curve and relation of its derivatives to the results

text: ![enter image description here](http://i.stack.imgur.com/wZwEX.png) I was plotting a deviation of my data points from a Gaussian smoothed curve depending on different properties of this Gaussian curve that smooths the series. So I had a series of:               y - result     g0 - value of gaussian smooth     g1 - first derivative of the gaussian curve     g2 - 2nd derivative     g3 - 3rd derivative      I was experimenting with different relations and I came to this plot:               plot(data[,"g1"] + 100*data[,"g2"] + 5000*data[,"g3"], data[,"y"]-data[,"g0"])     this is plotting Distane(g0,y) <- g1 + 100*g2 + 5000*g3      It gave me a fantastic, inspiring correlation. While the coefficients 100 and 5000 were selected experimentally to find the slimmest cloud of points, I'm unable to identify the reason for these particular coefficients. It would be a real breakthrough for me to calculate them from the data themselves rather than based on my subjective visual identification. It isn't said that this is the best or optimal setting. Any ideas or explanation for these settings? Any theoretical reasoning behind this phenomenon? ![coefficients 10 and 2000](http://i.stack.imgur.com/3WIGU.png)   For instance, for coefficients 10 and 2000, this will not look much different from just a simple plot of `Distance(g0, Y) <- g1`. So `Distance(g0, Y) <- g1 + 10*g2 + 2000*g3` and `Distance(g0, Y) <- g1` look pretty much the same.


title: How to read result from Kolmogorov-Smirnov one-sample test in SPSS?

text: I'm newbie and bit lost... Can someone help me read this result from SPSS 1-sample Kolmogorov-Smirnov test (poisson distribution).                                                       SUBKEY1   SUBKEY2  SUBKEY3  SUBKEY4    SUBKEY5          N                                128       128       128      128    128          Poisson Parametera,,b     Mean          .4609     .4609    .4922     .5156  .4922          Most Extreme Differences  Absolute  .092      .092     .103      .113   .103                                Positive  .079      .079     .088      .095   .088                                Negative     -.092     -.092    -.103     -.113      -.103          Kolmogorov-Smirnov Z                1.037    1.037    1.171     1.276   1.171          Asymp. Sig. (2-tailed)              .233      .233     .129      .077   .129      From the result above, I define (in simple words):               Ho = satisfied     Ha = not satisfied      but I don't know how to specify the level significance for this test and I don't know whether have to read Asymp. Sig. (2-tailed) (p-value) or Kolmogorov-Smirnov Z (D value) to be compare with critical value...


title: Approximate a poisson binomial distribution with a binomial distribution?

text: I have samples of Bernoulli distributed variable that are neither the first nor the second i in iid. My goal is to model their sum. I got from Wikipedia that I can use the poisson binomial distribution to make up for one of the i's, but then I have to keep all the inidividual probabilities. It would probably also be possible to throw the central limit theorem against it somehow to model it as a Gaussian, but I wonder if I can do better. Are there any results on how well a binomial distribution fits the sum of non identically non independently distributed Bernoulli samples. Especially if I can get some bounds on the accuracy wrt the correlation of the samples or something like that.


title: How to interpret the divergence of Fisher information expectation?

text: Consider translated Weibull distribution with probability density function: $$ f(x ; k, \lambda, \theta) = \frac{k}{\lambda} \left( \frac{x-\theta}{\lambda} \right)^{k-1} \exp\left( - \left(\frac{x-\theta}{\lambda} \right)^k \right) \chi_{x \ge \theta}(x) $$ Consider $ \mathbb{E}\left( -\partial_\theta^2 \log f(x; k, \lambda, \theta) \right) = \mathbb{E}\left( \frac{k-1}{(x-\theta)^2 } \left(1 + \left( \frac{x-\theta}{\lambda} \right)^k \right) \right)$. Notice that the expectation is only convergent for $k > 2$. **Question** : What is the significance, or interpretation, of infinite matrix element of the Fisher information matrix ? Does it mean that the maximum likelihood estimator for $\theta$ parameter is not asymptotically normal ? Thank you for any light shed on the subject.


title: Conditionals for Gibbs sampling of relational clusters

text: I'm trying to implement a Gibbs sampler, but I'm having trouble to find some of the conditionals of this model. ## Model We have $A$ actors, $K$ classes or clusters, and a matrix $\phi$ that determines the probabilities of two classes interacting. The generative process would be like the following: First, every actor chooses a class from a set of available classes. Then, every actor makes a set of $N$ interactions with other actors. To create every interaction, the user first chooses a class to interact with, and then an actor in that class.   1. Generate probabilities of every class: $$\pi \sim Dirichlet(\alpha_1, \alpha_2...\alpha_R)$$   2. Assign a class to every actor: $$z_i \sim Categorical(\pi)$$   3. Generate a relationships pattern for every class: $$\phi_{k} \sim Dirichlet(\beta)$$   4. For every actor $a$, repeat $N$ times: 5.1 Choose a class to interact with: $$r_{an} \sim Categorical(\phi_{z_i})$$ 5.2 Interact with an actor $x$ in that class: $$u_{an} \sim Uniform(p=\frac{1}{C_{r_{an}}}) $$ where $C$ counts the the number of actors in a given class. The uniform distribution means that, given a class choice, the probability of any actor from this class to be chosen is $\frac{1}{C_{r_{an}}}$. The actors out of this class have probability 0. Here is the graphical model: ![Probabilistic Graphical Model](http://i.stack.imgur.com/g4MGk.png) I want to sample the class assignments $z_j$ and find the probability of interactions between classes $\phi_k$. This is the joint probability distribution (I leave the brackets to show the plates in the graph model): $$ \prod_{i=1}^R\pi_{i}^{\alpha_i -1} \prod_{a=1}^{A} \left[ \prod_{i=1}^{R} \pi_{i}^{z_{ai}} \prod_{n=1}^{N} \left[ \prod_{k=1}^{R} \phi_{z_ak}^{r_{ank}} \frac{1}{C_{r_{an}}} \right] \right] \prod_{k=1}^{R} \left[ \prod_{i=1}^R \phi_{k_i}^{\beta_i-1} \right] $$ which, after a couple of simple steps, as the two $\phi$ correspond to a Dirichlet and a Categorical, can be converted to ( **option A** ): $$ \prod_{i=1}^{R}\pi_{i}^{C_i + \alpha_i -1} \prod_{a=1}^{A} \prod_{n=1}^{N} \prod_{k=1}^{R} \phi_{z_ak}^{r_{nk}} \frac{1}{C_{r_{an}}} \prod_{k=1}^{R} \prod_{i=1}^R \phi_{k_i}^{\beta_i-1} $$ where $C_i$ indicates the number of authors in class $i$. This equation can also be expressed as: $$ \prod_{i=1}^{R}\pi_{i}^{C_i + \alpha_i -1} \prod_{k=1}^{R} \prod_{l=1}^{R} \phi_{z_k,l}^{Q_{kl}} \prod_{a=1}^{A} \prod_{n=1}^{N} \frac{1}{C_{r_{an}}} \prod_{k=1}^{R} \prod_{i=1}^R \phi_{k_i}^{\beta_i-1} $$ where $Q_{kl}$ counts the number of interactions between a user of class $k$ to a user of class $j$. If we merge the two $\phi$ (easy because, again, they come from conjugates) we would have another expression of the joint distribution ( **option B** ): $$ \prod_{i=1}^{R}\pi_{i}^{C_i + \alpha_i -1} \prod_{a=1}^{A} \prod_{n=1}^{N} \frac{1}{C_{r_{an}}} \prod_{k=1}^{R} \prod_{i=1}^R \phi_{k_i}^{Q_{kl}+ \beta_i-1} $$ ## Question: I tried two ways to get to the conditionals: **Option A** : If I go on with the first expression, I can integrate out $\pi$ since this is a Dirichlet-like factor. Then, I would do the same with $r_{an}$, as they are $A\times N$ variables and I am not interested on them: $$ P(\boldsymbol{u}, \boldsymbol{z}, \boldsymbol{\phi} | \alpha, \beta) \propto \frac{ \Gamma(\sum_{r=1}^{R} C_r + \alpha_r)} {\prod_{r=1}^{R} \Gamma(C_r + \alpha_r)} \prod_{k=1}^{R} ( \sum_{x=1}^R \phi_{k,x} \frac{1}{C_{x}} )^{C_k N} \prod_{k=1}^{R} \prod_{i=1}^R \phi_{k_i}^{\beta_i-1} $$ This gives me the following conditional of $z_j$: $$ P(z_j | P(\boldsymbol{u}^{(-j)}, \boldsymbol{z}^{(-j)}, \boldsymbol{\phi} | \alpha, \beta) \propto \frac{ \prod_{k=1}^{R} ( \sum_{r_{an}=1}^R \phi_{k,r_{an}} \frac{1}{C_{r_{an}}} )^{C_k N} } { (C_x + \alpha_x -1) \prod_{k=1}^{R} ( \sum_{r_{an}=1}^R \phi_{k,r_{an}} \frac{1}{C_{r_{an}}^{(-j)}} )^{C_k^{(-j)} N} } $$ However, I can't get the conditional of $\phi$. **Option B** : If I go on with the expression B of the joint distribution, I am not able to get the conditional of $z$ because of the counts $Q_{ij}$. Is this model overcomplicated, or I am missing something?


title: Cross-validation and feature selection of a multivariate regression

text: I've been trying to create a multivariate regression model to fit my training data into the prediction of a value. I've put my data into a matrix `X` with `m x n` where `m` is the number of instances and `n` the number of features/predictors. My label vector is then `m x 1`. This is my code to predict the theta values, or parameters.               theta_matrix =  pinv(X'*X)*X'*y_label;      Now I want to slip the data into train and test, and by researching I've found that `cross-validation` in 10-fold can be a good option. If I do so, wouldn't it get me 10 sets of parameters theta? So what to choose from then? And about feature selection, I've found that `stepwise` can be a good choice, but I think it does not take into account that features can be correlated. Any alternative?


title: Expected value of coin tossed twice

text: Here's the problem: A (Laplace) coin with sides $0$ and $1$ is tossed twice, each time independently from the other. Let $f:=$ maximum of both results $g:=$ sum of both results. What is the expected value and variance of $f$ and $g$? I think $E[g]=2*\frac{1}{4} + 1*\frac{1}{2}+0*\frac{1}{2}+0*\frac{1}{4}=1$, but I'm not sure how to go about calculating the rest. Any guidance would be appreciated!


title: Randomly sample bounds from many multi-dimensional points

text: Let $\mathbb{M}$ be an $m\times n$ matrix, with $m < n$. My data consists of a large list of points $\mathbf{x}_i\in\mathbb{R}^n$, $i=1,...,N$, where each point satisfies $\mathbb{M}\cdot\mathbf{x}_i = 0$ and $\mathbf{A} < \mathbf{x} < \mathbf{B}$, where vector inequality denotes all component-wise inequality, and $\mathbf{A},\mathbf{B}$ are some vector bounds. For each vector $\mathbf{b}\in\mathbb{R}^n$, let $\mu(\mathbf{b})$ denote the number of points $\mathbf{x}_i$ that satisfy $\mathbf{x}_i \leq \mathbf{b}$. Obviously $\mu(\mathbf{A}) = 0$ and $\mu(\mathbf{B}) = N$. Naively taking a random vector $\mathbf{b}$ inside $\mathbf{A} < \mathbf{b} < \mathbf{B}$ will surely yield very low values of $\mu(\mathbf{b})$. Similarly, $\mu(\mathbf{x}_i) = 1$ for most points $\mathbf{x}_i$ I've tried. So here's my problem. I need a way to randomly sample vectors $\mathbf{b}$ such that all values of $\mu$ are more or less equally represented.


title: How to analyse means between two similar scales?

text: In a questionnaire I conducted for my dissertation, i asked respondents to rate the frequency (5 point scale from never to always) to which they resorted to certain unethical behaviour during their work. I also asked them to rate their acceptance of such behaviour using a reliable scale (5 point scale from Strongly Disagree to Strongly Agree). Now, as expected, the mean scores of the self-admitted unethical behaviour is lower than the mean of the behaviour acceptance scale. Is there an appropriate statistical test to use in this case? What test is there to prove that the behaviour acceptance scale is a good measure/representative of actual behaviour? Is correlation enough?


title: Basic algebra question regarding properties of the summation operator

text: I admit that my algebra is rusty. Can someone please explain how the middle part of the top equation is transformed to the middle part of the bottom equation (both are highlighted in black). ![equation](http://imgur.com/DqgVr26.jpg)


title: Calculating 90% confidence intervals in regression to identify outlying data

text: I have been searching for the answer for a problem with confidence intervals for a long time now, so I hope someone can help! The two psychological papers below performed the same analysis, regressing irregular word reading onto nonword reading for controls, and plotting these regressions with 90% confidence intervals, before superimposing dyslexic participants to find outliers. However, one paper has straight line intervals, the other curved intervals. Why? (I have provided links to the papers so you can observe the plots). I have scoured both textbooks and online resources, and the only descriptions of confidence intervals I have found are for confidence intervals of the regression line, which estimate an area in which there is a given probability that the regression line will fall. This is not what the studies below used, however. They used what may, I think, be termed ‘prediction intervals’ and describe an area in which a given percentage of the data will fall (which are then, in this kind of research, used to denote an area of 'normal' performance on a task). Can anyone help with these distinctions, and/or provide a calculation for this type of interval? The curious thing is that SPSS/PASW will plot these confidence intervals for me, but offers no explanation of where they come from, nor can I edit the provided plot. Any help, no matter how small, will be greatly appreciated! References: Manis, F. R., Seidenberg, M. S., Doi, L. M., McBride-Chang, C. & Petersen, A. (1996). On the bases of two subtypes of developmental dyslexia. Cognition: 58, 157-195. LINK: http://lcnl.wisc.edu/publications/archive/45.pdf Stanovich, K.E., Siegel, L.S. & Gottardo, A. (1997). Converging evidence for phonological and surface subtypes of dyslexia. Journal of Educational Psychology: 89, 114-127. LINK: http://web.mac.com/kstanovich/Site/Research_on_Reading_files/RdJEdPsy97.pdf


title: How should I deal with the consequences of proportional hazards violations in log-rank (and related) tests?

text: I have (censored) time-to event data for subjects in four groups. I would like to do something like a logrank test, but the survival curves do not satisfy the proportional hazards assumption. I think I have heard that the consequence of a prop. hazards violation is loss of power to identify differences between survival curves, however, in the case of my study, the logrank tests _do_ show significant differences between the groups. The proportional hazards violation takes the following form: one group has a relatively larger probability of an event early in the observation period, and another group has a relatively larger probability of an event late in the observation period. I believe that the G-$\rho$ family of tests (for example in the `survdiff` function in R's survival package) can be parameterized such that the earlier or later portion of the observation period is more heavily weighted. However, in this case, the different groups would "do better" (the event being studied is a good thing) in different specifications for the test. I would like to know several things:   1. if the logrank function _does_ find significant difference despite the presence of a proportional hazards violation, can we interpret this as indicating a true (overall or average) difference between the curves? or does this violation mean the test's results are totally meaningless?   2. is there a principled way of describing survival times in cases like mine? I would ideally like to be able to report an overall hazard ratio (I know -- this would lack external validity with the non-proportional hazards and censored observations, but would be useful in describing the experiment), as well as give information about which groups were more likely to have events at which times. I could choose a break point in the middle of the observation period and just do separate tests before and after (assuming these subsets of the data did satisfy the proportional hazards assumption), but the choice of such a point feels somewhat ad hoc. Related discussion: This thread discusses alternatives to the logrank test, but doesn't consider my issue: What are the pros and cons of using the logrank vs. the Mantel- Haenszel method for computing the Hazard Ratio in survival analysis?


title: Including or discarding a variable based on its frequency

text: I am trying to dummy code a variable sector which have 20 levels. I am trying to combine some of levels. However I got this situation, if I have a say 2000 observations, how do I decide on the minimum frequency of a particular level that I should consider before making it a dummy variable, as variable with low frequency not add much information. How do I decide on what minimum % of total number it should be, in order to take into modeling?


title: Using and understanding the zelig package in R

text: Is any one here familiar with an R package called Zelig? I have a data frame like this:               IQ   AGE     80   50     100  18     90   25      etc. What I need to do is build a model of IQ given AGE, I am running these commands: `z.out <- zelig(IQ~AGE,data=df,model="ls")` this runs the what-if given age 110, what would be the IQ `x.out <- setx(z.out, AGE=110)` This is a simulation model where given the age 110, after running 1 million runs of simulation, what would be the IQ with 95% confidence interval. `s.out <- sim(z.out,x.out, num=1000000, level=95)` I have a hard time understanding from what pool of data the `sim()` function draws the numbers. I read though the docs, but they are written for Ph.D. students, if not more advanced readers. I have asked the Zelig creators this question multiple times but they are directing me to the docs which I read multiple times, with no luck. However, one of the person that works with Zelig sent me this email: > Suppose that you fit $$\text{IQ} = a + \text{Age} * b + e$$ Then you get a > table of regression coefficients where a=50, b=2, and their standard errors > are something like $\text{s.e.}(a)=\sqrt{10}$ and $\text{s.e.}(b)=1$. These > are all hypothetical examples. In maximum likelihood estimation, this > regression output is another way of saying that $a$ and $b$ are distributed > bivariate normal with means $[50,2]$ and there's a variance-covariance > matrix that looks something like this (all numbers are made up): > $$\begin{array}{cc} 10 & cov(a,b) \\\ cov(a,b) & 1 \\\ \end{array}$$ So, the > variance of $a$ is 10, the variance of $b$ is 1, and their covariance is > $cov(a,b)$. It won't be shown in your regression table, but Zelig remembers > it for you. Let's pretend it's 3. This variance-covariance matrix is the > inverse of the Hessian I mentioned earlier. Don't worry about it. For this > example, you need only remember that $\text{mean}(a,b) = [50,2]$ and > $cov(a,b)=\begin{array}{cc}10&3\\\3&1\end{array}$. For purposes of plain > text email, I'm representing matrices with columns separated by commas and > rows by semicolons. In addition, suppose that the error term $e$ is > distributed with mean 0 and s.d.=1. Now, one way to predict what IQ you > might get for somebody aged 88, based on this regression table, is exactly > what you would expect: you simply calculate 50 + 88 * 2 = 226. This is your > point estimate. The 95% confidence interval around this point estimate is a > function of the standard errors of the coefficient estimates of $a=50$ and > $b=2$, and the exact formula for that is in any econometrics textbook. > Simulation makes it unnecessary to dig up that textbook. Instead, for 1000 > rounds, `sim()` will come up with 1000 different pairs of $(a,b)$ estimates > drawn from the bivariate normal with mean=[50,2] and > cov=$\begin{array}{cc}10&3\\\3&1\end{array}$. One such pair might be > $(47,1.5)$; another might be $(52,3)$; yet another might be $(10,5)$. > Whatever they are, `sim()` plugs them into the formula and gives you 1000 > different estimates for the IQ. Their average is your point estimate. If you > stack them from lowest to highest, the ends of the 95% confidence interval > are the top 25th value and the bottom 25th. That's it. That's all that > `sim()` does. Given the above explanation, can anybody tell me in lay terms, what numbers `sim()` is picking? How are those numbers in pool generated? I would greatly appreciate if anyone brings some light into this.


title: Relationships between t distributions and normal distributions

text: In Gelman's Bayesian Data Analysis: > The t distribution is **the marginal posterior distribution for the normal > mean with unknown variance** and **conjugate prior distribution** and can be > interpreted as **a mixture of normals with common mean and variances that > follow an inverse-gamma distribution.**   * What precisely (in terms of mathematical formulas) does it mean by "the t distribution is the marginal posterior distribution for the normal mean with unknown variance"?   * In the table on Wikipedia, I didn't find that the t distribution is a conjugate prior distribution for a normal distribution. So where does "conjugate prior distribution" come from?


title: how to implement linear or non linear regression for 3d position estimation?

text: I am a beginner in Machine Learning. For my project I need a regression algorithm that can estimate the 3D position of a device based on some constraints (moreover inputs). I know how to implement regression algorithms (using gradient descent) for scalar outputs. But in this case the output is a position vector. Can someone provide me any material or book chapter. Any quick theoretical reference? Regards.


title: What is the difference between time-varying coefficients and time-varying variables?

text: In survival analysis, what is the difference between time-varying coefficients and time-varying variables? Suppose we have the following Cox-Proportional Hazards Model: $$\log[h(t)] = \log[h_{0}(t)]+ \beta_{1}x_1(t)$$ vs $$\log[h(t)] = \log[h_{0}(t)]+ \beta_{1}(t)x_1$$ In the first model, $x_1$ is a time varying covariate. In the second model, $\beta_1$ is changing in time. How do we interpret this?


title: Is it possible to have an F-test with a $p<.05$ even if there is an $R^2<.1$ for a least squares analysis in JMP?

text: I've been working with JMP, and I have found that often I get a significant p value but a relatively small r. Is this supposed to happen?


title: Using inverse of cube In linear model

text: What's the formula for a equation that can produce the continuum from the red to green lines in this graph below? I can easily get anywhere from the green line to the blue with $$y = B_0 + B_1x + B_2x^3$$ but getting anywhere from the blue line to the red line requires an inverse of the cube or x^(1/3) which doesn't work for negative numbers. Is there any continuous formula that can produce the red line? Alternatively, how do I use imaginary numbers in a linear regression?   ![enter image description here](http://i.stack.imgur.com/z8T0c.png) Here's some data if anyone can fit a model to it.                 x, y      -10,  2.154434525      -9,  2.080083671      -8,  1.999999861      -7,  1.912931059      -6,  1.817120484      -5,  1.709975855      -4,  1.587400979      -3,  1.442249517      -2,  1.259921021      -1,  1       0,  0       1, -1       2, -1.259921021       3, -1.442249517       4, -1.587400979       5, -1.709975855       6, -1.817120484       7, -1.912931059       8, -1.999999861       9, -2.080083671      10, -2.154434525


title: Bootstrapped optimism with low prevalence

text: My dataset is not quite large enough to be split into training and testing sets, so I am using bootstrapped optimism to account for overfitting when reporting model performance. The cases in my data set have a low prevalence, such that in some iterations of sampling they're missed. To account for this I've increased the number of times I am sampling and consider set the optimism, for those instances that didn't sample enough (or any) cases for the model to converge, as NA. **I wanted to get a sense for whether that made sense, or if anyone has any experience doing a stratified sampling approach for bootstrapped optimism.** I am using Frank Harrell's 1996 stats medicine paper (Tutorial in biostatistics: Multivariable Prognostic models..") as a guide. I know that RMS has the validate function that can do this for lm and cph models, but I am using neither of these.


title: Alternative to Chi-Square, etc for categorical variables

text: I am trying again to explain my problem, and I have more concrete results. I have tried every potentially relevant test I can find, and none seem to work quite right. The question is illustrated here: ![enter image description here](http://i.stack.imgur.com/fAk1q.jpg) The actual data are protein sequences. I have two categories of specimens, shown by the black and grey circles, and many positions of genes, shown by the colored squares. I want to find which of the gene positions correlate with the specimen categories. For example in row 1, there is perfect correspondence. Nearly all metrics pick up this signal. In row 5, the genes are exactly evenly distributed between the traits, and again most metrics get this. However, in case 4, there is zero predictive power for that trait, (every position has to have one or the other value, so expected value should be 1), but most tests report a high value. I have a script which calculates metrics for test cases and actual data, and some of these results are illustrated below. The data set has been rotated so the rows are now columns, and the two categories are indicated by the colors of the letters, with two black rows and 4 grey rows. ![enter image description here](http://i.stack.imgur.com/wMAp1.png) Yates-corrected Chi-square, Mutual Information, Fisher's exact test, Concordance (with some assumptions), etc, all fail on case 4, and some of them give case 2 an equal value to case 1, or score case 3 higher than case 2. I came up with a metric (MyChi, the rust-colored line above) which is an attempt at a one-tailed Chi-square (only adds values higher than expected, and not lower), normalized by the number of traits present, and not counting the singleton traits. It works pretty well (the columns in the second figure are ordered by how well they "should" score), but I am leery of just making up a test. The actual data also have a limitation that there will always be lots of zeroes and ones in a contingency table, so many of the standard tests are not recommended. Although there are many rows per sample (and I would like to figure out a way to sum the scores from all rows into one metric), there might only be 20 values for each row. Does anyone have a recommendation for a test to address this analysis?


title: Should I bootstrap at the cluster level or the individual level?

text: I have a survival model with patients nested in hospitals that includes a random-effect for the hospitals. The random effect is gamma-distributed, and I am trying to report the 'relevance' of this term on a scale that is easily understood. I have found the following references which use the Median Hazard Ratio (a bit like the Median Odds Ratio), and calculated this. Bengtsson T, Dribe M: Historical Methods 43:15, 2010 However, now I wish to report the uncertainty associated with this estimate using the bootstrap. The data is survival data, and hence there are multiple observations per patient, and multiple patients per hospital. It seems obvious that I need to cluster the patient observations when re-sampling. But I don't know if I should cluster the hospitals too (i.e. resample hospitals, rather than patients? I am wondering if the answer depends on the parameter of interest, and so would be different if the target was something that was relevant at the patient level rather than the hospital level? I have listed the stata code below in case that helps.               cap program drop est_mhr     program define est_mhr, rclass     stcox patient_var1 patient_var2 ///         , shared(hospital) ///         noshow     local twoinvtheta2 = 2 / (e(theta)^2)     local mhr = exp(sqrt(2*e(theta))*invF(`twoinvtheta2',`twoinvtheta2',0.75))     return scalar mhr = `mhr'     end          bootstrap r(mhr), reps(50) cluster(hospital): est_mhr


title: Output probabilities of binary support vector machine classifier in Matlab R2014a

text: I’m using SVM for classification of my binary output problem. I want probability of belonging to every class. How can I obtain it? For instance suppose this is our structure:               X = [1 4 1 5 3 ; 4 5 4 3 1 ; 6 7 5 5 5];     %% (Every row is a sample and every column is a feature)     Y=[1 ; 0 ; 0 ; 1 ; 1];      After designing SVM model I want check an outsample data with it. For example:               Out_sample_data = [3 3 4 0 3 ; 7 5 7 1 2 ; 6 7 3 3 4];      Besides out sample data classes, I want find probability of belonging to every class. How can I do this with MATLAB R2014a ? Thanks. Ps. in Matlab R2014a we have some new functions like "fitSVMPosterior" that I think we can use it for this purpose.


title: Regression Analysis on Segments (Subsets of Dimensions)

text: My data has subsets of 8 categorical dimensions that are mapped to a real valued dependent variable. So if my dimensions are X1, X2,.., X8 my data might look like this: (X1=x11, X3=x31, X8=x81) -> r1 (X4=x42) -> r2 (X1=x13, X5=x53, X7=x73) -> r3 ... So, the difference between this scenario, and your run of the mill regression, is that the dimensions are incomplete. I am wondering what is a usual way to proceed in this type of a scenario. One thing that comes to mind is to add the missing dimensions with a new dummy categorical variable, and then train the model. But I am not sure this would be the correct way to proceed. Any help is much appreciated.


title: Testing for membership in a distribution with one single value

text: I came across an interesting little quirk when I was coding up a hypothesis test in Python. I was wondering why very time I tested some values for normalcy using the Anderson–Darling test, I was getting a positive. Turns out my loop was just testing over one single value (oops)! But then I got to thinking, I wonder if this is generally the case -- can a single real- valued number be shown to belong to every distribution that has real-valued numbers as it's support? Is there a trivial result that states this somewhere in the statistical literature?


title: Using cross correlation to infer dependence, can it be done?

text: I have a very particular question, I have seen a similar one here, but my knowledge is too limited to make use of it. I will try to explain myself as clearly as possible... Wish me luck! I have a sequence of ( _a priori_ ) probabilities of a binary variable, we could say $X_i \sim Ber(p_i)$, i.e. $P(X_i = 1) = 1 - P(X_i = 0)=p_i$. I know that $X_i$ are not independent, but I don't know exactly $P(X_i|X_{j\neq i})$ nor $P(X_i|X_{j_1\neq i},...,X_{j_n\neq i})$ (I don't know anything about their dependence), and I want to know $$P(\sum_{j=-K}^{j=K} X_{i+j} =0)$$ What kind of knowledge/hypothesis do you think I need in order to approximate this more efficiently? Assuming that they are independent does not work well enough. What would you do if you found a similar situation? I am not an expert at all in this matters, and the more I learn about statistics the less I know! For example, I can see that the probabilities form small "triangle" shapes, so maybe something like $P(X_i = 1 | X_{i-1} = 1, X_{i-2} = 0) = P(X_i = 0 | X_{i-1} = 0, X_{i-2} = 1) $ can help? If so, is there any way to use this? As I said, I have no idea... Thank you very much for your help! edit: I think the title is not very good, but I don't know how to explain it better... One more question, can I say $Y_i = X_i - X_{i-1}$ and try to see $P(X_i|Y_{i-1})$?


title: Gaussian process estimation

text: The stochastic process $(X_t)_{t\in T}$ is called Gaussian if for all $t_1,\dots,t_k\in T$, for all $k$, the joint distribution of $X_{t_1},\dots,X_{t_k}$ is multivariate normal. The process is completely characterized by its mean $$\mu(t) = E[X_t]$$ and covariance functions $$\sigma(s,t) = Cov[X_s,X_t]$$. Given a centered (0-mean) Gaussian process, is it possible to estimate its covariance function?


title: Can I compare results from the Kruskal-Wallis test with different numbers of inputs?

text: I have repeatedly sampled 9 populations and calculated the error from a "true" estimate. I have done this where I have taken progressively less samples (A-G). So for example in A I have taken 10,000 samples each from 9 different populations, calculated a statistic for each population and compared them to the true value, for B I have taken 5,000 samples and so on. I would like to know at what level of sampling the error becomes significantly different from the error in strategy A.               A = [ 0.001, 0.062,  0.014, 0.03,  0.087, 0.071,  0.005,  0.018,  0.058]     B = [ 0.001, 0.083,  0.041, 0.048, 0.153, 0.072,  0.035,  0.007,  0.143]     C = [ 0.002, 0.008,  0.04,  0.019, 0.083, 0.004,  0.011,  0.008,  0.132]     D = [ 0.039, 0.262,  0.089, 0.082, 0.06,  0.013,  0.013,  0.049,  0.083]     E = [ 0.11,  0.079,  0.032, 0.127, 1.037, 0.039,  0.063,  0.245,  0.772]     F = [ 0.262, 0.784,  0.095, 0.438, 0.937, 0.226,  0.269,  0.626,  1.161]     G = [ 1.003, 2.173,  0.545, 5.417, 0.944, 1.958,  1.795,  3.053,  5.168]      Using `scipy.stats.kruskal(A, B)` returns a significance value of `0.33` therefore there is no difference between the two strategies. I keep adding strategies e.g. `scipy.stats.kruskal(A, B, C)` until I get a significant answer (p<0.05) which occurs when I use A - F. Why question is, is the Kruskal-Wallis test the right test to use and can I compare the answers from each test? Or would it be better to use some other test, maybe the t-test between pairs? And if the KW test is OK is the above approach valid?


title: Interpreting the results of a pca

text: I applied pca on r using `prcomp`. I would like to use the main PCs to reduce the size of my problem. In order to do so, I want to express the values of my samples on these principal components. As far as I understand:   1. The rotation matrix expresses the weighting of the variables found: it returns, for each variable, the weighting of the each sample on this variable. This would be the values of the samples on my new features and would replace my old data matrix.   2. `pr.result$x` displays the coefficients of the linear combination of my new features with respect to the old ones Question 1: Are the weightings expressed by my rotation matrix the same as the linear combination of the original data matrix taking the coefficients of `pr.result$x`? i.e If I take all the coefficients of the samples on the old variables and perform the linear combination with the coefficients of `pr.result$x`, would I get the same result as the one given by my rotation matrix? Question 2: When performing my pca, I scaled my data (`scale=TRUE` in `prcomp`). Is there an easy way to apply the inverse transform to the rotation matrix?