title: GLMM Models in R

text: I want to see if there is a significant relationship between the mating (response variable) and several explanatory variables such as individual number (categorical), leg movement (continuous) and the reuse of individuals (categorical). My data looks like this:               Mating  ID  Replication Shaking     1       M1  R1          10     0       M1  R2          14     0       M2  R1          15     0       M2  R2          12     1       M3  R1          14     0       M3  R2          17     1       M4  R1          19     0       M4  R2          22     1       M5  R1          18     0       M5  R2          16     1       M6  R1          17      1 means mating happened, 0 means did not occur. I am trying to organize the data to apply a GLMM to the data, but have not been able to do it. I followed the model proposed by Crawley for binary response with pseudorreplication but it does not work. The script goes like this:               model2 <- lmer(Mating~Replication+(1|ID), family=binomial, method="PQL")      Can somebody tell me what part is wrong and how could I fix it? If you know a better method, that will be really welcomed as well.


title: Standard deviation of a ratio (percentage change)

text: I have 2 data sets. The first data set, let's call it $X$ has an average value of ($\bar X$) and standard deviation of ($STD_X$), the second set of data also has the average value of ($\bar Y$) and standard deviation of ($STD_Y$). I want to find out the standard error or standard deviation of a percentage change of data set 2 compared to data set 1. So I have $((\bar Y-\bar X)/\bar X)*100$. Now my question is, how do you take into account the standard deviations for this percentage value?


title: Time series modeling of choppy data

text: I'm trying to model 10 years of monthly time series data that is very choppy, and overall it has an upward trend. At first glance it looks like a strong seasonal series, however the test results indicate that it is definitely not seasonal. This is a pricing variable that I'm trying to model as a function of macroeconomic environment, such as interest rates and yield curves. I've tried linear OLS regression (proc reg), but I don't get a very good model with that. I've also tried autoregressive error models (proc autoreg), but it captures 7 lags of the error term as significant factors. I don't really want to include that many lag of the error term in the model. In addition, most of the macroeconomic variables become insignificant when I include all these error lags in the model. Any suggestions on modeling method/technique that could help me model this choppy data is really appreciated.


title: Finding similar items based on user likes and dislikes

text: I have data on content-based recommendations of movies and their attributes. Suppose a user likes x,y and z and also dislikes c and d movies. I want to predict movies that he will like based on his likes and dislikes. It is actually quite easy to find similar items if the only interaction were a movie like, because I would then be looking at the closest items to x movie by calculating distances via attributes. But it is quite confusing when interactions are multiple and based on both likes and dislikes. What is the right approach in that case ?


title: Kernel density estimation - application

text: What is the validity of using a _kernel-density-estimation_ to _compare model x observed data_? In other words, if the _KDE_ curve for the observed data looks like the _KDE_ for the model forecast, can I use this result as a quality measure for the model? If not, why? In this case, it's a weather forecast model - but I have some complex variables like _wind direction,_ for which a simple correlation will not help too much. For example, the wind going from 360 to 5 is a huge numeric decrease, but it's almost the same direction! I deeply appreciate any help on this, Thanks!


title: Error Analysis in ANN (modeling)

text: I am using Artificial Neural Networks (FFBP, RBF, GRN) and Multi Linear Regression techniques to model a real world case study. I have a small list of records (60 data items). I am using 48 of 60 as my training data set and the remaining 12 as my testing data set. I have two diagrams: one is based on actual observations for my test set (49 to 60), the other diagram is based on the estimations from the models for my test set (49 to 60). I have three questions:   1. What is the best error range for this models (short record data)?   2. How can I calculate the error amount that can be added and subtracted with observation to form the error range?   3. What is the level of confidence? I would appreciate if you can direct me to similar examples (or some tutorial) where I can learn to find the answers to my questions. ![enter image description here](http://i.stack.imgur.com/cChXx.jpg)


title: How would you deal with categorical data in a naive Bayesian classifier?

text: I've built a little naive Bayesian classifier that works with Boolean and real values. Boolean distributions are dealt with via Bernoulli distributions, while real valued data are dealt with kernel mixture estimators. I'm currently in the process of adding count data in. How would one deal with categorical data though, e.g. Monday, Tuesday, Wednesday, or Toyota, Honda, Ford? My initial thoughts are to assign a number to each category, treat it as a normal real value and round to the nearest integer category on prediction. That seems very wrong to me though.


title: Variance explained from factors for step-wise regression in R

text: I used a stepwise multiple regression to generate a model and I am trying to appropriately report the results by indicating the (additional) variance explained by each included factor. I'm not sure how to extract that information, however. Basically my code thus far is simply:               init <- lm(dep ~ fac1 + fac2 + fac3 + fac4,data=data)     final <- stepAIC(init, direction="both")      or               final <- step(init, direction="both")      In the case of my data, the final model is `final ~ fac1 + fac2 + fac3` At this point, I want to know how much variance the most important factor explains, how much additional variance is explained by adding a second factor and so on. I can get the coefficients, but I'm not sure how to extract the R/R^2 for each factor. Where is that reported? Thanks in advance.


title: Distribution of transformation

text: Suppose $X_1,\ldots,X_n$ are i.i.d. $\mathcal U(0,1)$. I am looking for the asymptotic distribution of $$T_n = \prod_{i=1}^n [e{X_i}]^{1/\sqrt{n}} \>.$$


title: How are standard errors affected in a multivariate regression?

text: Initially we have this regression: $$ \text{hourly wage} = 12.69 + 5.44\text{CollegeEducation} - 2.64\text{Female} $$ with standard errors of $0.21$ on $\text{CollegeEducation}$ and $0.20$ on $\text{Female}$. The variable $\text{Age}$ is added to the regression, which has a coefficient $0.29$, but the standards errors on college education and female stay the same. Dummies for four locations are then added, they have coefficients ranging from $-0.2$ to $0.6$, the standard errors on college education and female stay the same. The questions are:   * Why don't the standard errors change?   * I have read somewhere that added uncorrelated controls can increases accuracy and lower standard errors, is this true?


title: Find correlating hashtags

text: I have a couple of user profiles, that contain hashtags like this: user1: #sports, #bike, #enduro, #racing, #mountainbike user2: #sports, #bike, #racing, #mountainbike, #mtb (...) userN: #mtb I would need to find an algorithm, that constantly builds up correlating word groups out of all the hashtags in my database. In the end I would like to better match users by comparing not just the hashtags because they might all have different writings, I would like to compare if those hashtags belong to a certain group of words. I hope my question is understandable :)


title: Ordered Probit and categorical variables

text: I wanted to run a quick and easy (or so I thought) regression on some data I have but now I am starting to doubt whether or not the regression makes any sense. I have seen some similar questions but they don't really answer what I want to know. I want to run an ordered probit/logit with the dependent variable having five ordinal outcomes. My explanatory variables are: age, work experience, gender, position with-in the organisation, country of origin and so on. "Gender" will be a dummy variable and "Country of origin" will be coded using dummmies. "Age" is a categorical variable so it will have values such as 18-25 years, 25-32 years etc. Should I code each group as a dummy variable? "Position with-in the organisation" is categories such as employee, executive and so on. Should I code each of these as dummy variables? Work experience are in categories such as 1-2 years, 2-3 years and 5+ years etc. Should I code these as dummy variables? What I am interested in knowing is that when I only have dummy variables as explanatory variables, will I get meaningful results? For some reason I find it very unintuitive.


title: How to summarize and understand the reults of DBSCAN clustering on big data?

text: Many clustering algorithms can be used with big data, eg. versions of KMeans, DBSCAN based on Hadoop, etc. But, with k means we will get k centroids for k clusters and we can map them to the space and somehow understand the results. But what about density based algorithms like DBSCAN? In DBSCAN we'll get m clusters containing millions and billions of data points in case of big-data. **How do we try to understand the results of such clusterings?** Clustering more data is costly and also meaningless if we don't understand each cluster.


title: One-way ANOVA model for categorical data

text: I have data to fit a one-way ANOVA model to categorical data. I used a simple model first, getting significant coefficients, but with an $R^2$ of 0.2. I then tried a log model and got more significant coefficients, but the $R^2$ is still only 0.3. Is this normal for categorical data? Are there any other methods I should try to improve the model?


title: what I do with statically significant

text: I have a class with 52 students. Every student has a grade. I want to know how many of the student grade in compaire with a special number(for example 15) is statically significant?


title: What is the best test to estimate the correlation between binomial/categorical dataset?

text: I'm trying to analyze if there are correlations between binomial dataset. I have binomial data (presence/absence) of two variables in different periods and I need to know what is the best way to find some relationship. Data were taken annually from the same samples.


title: Modelling rating time series

text: Imagine a simple set up where where entities (e.g. a restaurant) are rated every year with 5 rating categories $c_t$: _excellent_ , _good_ , _medium_ , _poor_ and _bad_. Now assume that I have some underlying data $X^i_t=(x^{1i}_t,x^{2i}_t, ..., x^{pi}_t)$ at entity level $i$ on top of the rating history $c^i_t$. I am interested in a modelling approach to predict $c^i_t$ from $X^i_t$. In the past, I have randomly sampled each entity only once from the entire rating history (to minimise the correlation effects) and used an ordered logit model to estimate the coefficients of $X^i_t$. However, I believe I can make more efficient use of my data if I could use the entire data. Is there such a model? Could you please point to some literature?


title: How to select predictors in regression using backward method?

text: When is it appropriate to use a backward method in regression? I have read that it is permitted for exploratory model-building, but I have also read negative things about it. I am making a model about streaming. The 12 predictors I chose are based on theory. The enter method is not an option, because when I put all independent variables (sociodemographic and others) in the regression with the enter method, only one independent variable is significant. When I use this independent variable and add an independent variable (which was significant from backward selection), both are significant, which means that the enter method filters variables out of the regression that are significant.   1. Can I use the backward method to test all variables first, then use the significant variables to rerun the regression with an enter method, or must I be consistent and use the backward method again? Or is the backward method out of the question?   2. May I select independent variables that have a significant effect on the dependent variable according to the correlation matrix to run the regression with, in this case maybe, the enter method?


title: Considering correlation of independent variables in regression model

text: I have one dependent variable `C` and four independent variables:   * `S`   * `SD`   * `TPA`   * `VPT` The correlation matrix between them is the following                   TPA     VPT      S     TPA -       -       -     VPT -0.37   -       -     S   0.04    0.00    -     SD  0.05    -0.04   0.01      `TPA`and `VPT` have a moderate negative linear relationship. Now I want to create a regression model:               C ~ S + SD + VPT + TPA      Does it makes sense to consider the correlation between `VPT` and `TPA` in the model? If so, what should I add to it to consider the correlation?


title: Picking and permuting over randomized reference set from large amount of data

text: I have a question about different strategies for randomly picking up a reference set for statistical testing from a large amount of data (and permuting over the random picking/testing). For illustration purposes, we could assume that the data is e.g. time series data (could be something else as well), where different phenomena (either continuous or binary) has been collected densely over a long period of time. What I want to do is: Given a set of time periods (test set), see if the measured phenomena is enriched in the set when compared to the data (reference set). Practical example of this could be:   * There is a data set that contains temperature measurements measured every 1 second for the last 50 years   * Given a set of time periods (test set), ranging from a couple of seconds to a couple of days, test if the temperatures in the test set is statistically different from the temperatures in the general data.   * For the purpose of this question we can assume that the test itself is some basic statistical test such as t-test or linear regression. My naïve solution would be to calculate all permutations of test set time period length in the data, and compare the test set to those values. Because of the massive amounts of the data and permutations, this would not be computationally feasible (not to mention that the test set can contain time periods of greatly varying lengths). The second and more feasible option would be to randomly pick a reference set from the data, consisting of time periods of similar length as in the test set. It just sounds that doing this random selection of reference set once is not enough, the obvious choice to me is to perform the random selection and testing several times.. And then to the actual question, is this feasible at all, and if so, what should I actually do with the results from the repeated testing, assuming that the end result I want would be:   1. Measurement of statistical significance (i.e. p-value)   2. Measurement of effect size Yet again naïve solutions I can come up are:   * Take the mean of the p-value/t-score for each randomized test   * Set up some kind of threshold for the testing, e.g. p < 0.05 and calculate the proportion of tests that pass this criteria All help/ideas will be greatly appreciated! **TL;DR** _I have a large set of data, and want to test a set of given ranges against it through randomly selecting reference set multiple times. How do I combine results from these multiple tests into a single values._


title: Size or level $ \alpha $ in Neyman-Pearson

text: I found two different definitions of Neyman-Person lemma, one using the size $ \alpha $ and one the level $ \alpha $ of a test. I know the two concepts are related, but not synonymous. Is the lemma the same in both cases? Does it make any difference in the proof? Thank you!


title: What is the distribution of the sample variance of the Skellam distribution?

text: I want to estimate the parameter $\mu$ using the difference between two Poisson distributions with the same parameter, i.e. a Skellam distribution with $\mu_1=\mu_2 = \mu$. I can calculate the variance of the Skellam distribution as the average variance over multiple samples, however, depending on sample size and some good/bad luck, I have quite a bit of error in my estimate (red crosses in the figure below); the variance should be $2\mu$. If I knew what distribution I was looking for, I could fit that instead to hopefully get a more robust estimate. Why am I doing this if I appear to know $\mu$? In my application, I only know $\hat{\mu}=\beta_1\mu+\beta_0$, so I need a second estimate of $\mu$ to find out the betas. ![enter image description here](http://i.stack.imgur.com/sw2G0.png)


title: Computing ANOVA in R

text: Edit: I think the data was poorly given, varables coke, pepsi, and sprite were just given in 2 columns when I got it, but it has nothing to do with before and after treatment? here's the original questions since someone asked For the study, 2 males and 2 females are picked randomly from NYC, and a random sample of 2 males and 2 females are picked from LA. Each person is given 6 drinks (2 from each brand), where the order of the drinks for each person is randomized. Each person rates the drink’s taste, Y, on a scale of 60 to 100, where 100 is the best tasted. Please let me know if I can provide any other info. Pretty much I am asked to fit ANOVA. End of Edit# I would like to understand ANOVA in R, I am working in this data set.               dat <- structure(list(city = structure(c(2L, 2L, 2L, 2L, 1L, 1L, 1L,      1L), .Label = c("LA", "NY"), class = "factor"), sex = structure(c(2L,      2L, 1L, 1L, 2L, 2L, 1L, 1L), .Label = c("F", "M"), class = "factor"),          rater = 1:8, frit0 = c(77L, 72L, 78L, 84L, 77L, 78L, 83L,          81L), frit1 = c(76L, 71L, 78L, 84L, 76L, 77L, 86L, 84L),          utz0 = c(78L, 73L, 78L, 84L, 74L, 81L, 85L, 83L), utz1 = c(80L,          75L, 80L, 86L, 72L, 73L, 81L, 79L), weiss0 = c(78L, 76L,          84L, 90L, 78L, 79L, 88L, 86L), weiss1 = c(81L, 73L, 81L,          87L, 81L, 82L, 92L, 90L)), .Names = c("city", "sex", "rater",      "coke0", "coke1", "pepsi0", "pepsi1", "sprite0", "sprite1"),       class = "data.frame", row.names = c(NA, -8L))      I want to do an analysis of variance, check assumptions, and interpret results. I am sure this is done many times by the experts, so, what is the correct step to do it? I have tried               am1 <- aov(coke1 ~ coke0 + sex, data=dat2)     summary(am1)      Is this correct? Or even relevant? thanks in advance!


title: Computing the mode of data sampled from a continuous distribution

text: What are the best methods for fitting the 'mode' of data sampled from a continuous distribution? Since the mode is technically undefined (right?) for a continuous distribution, I'm really asking 'how do you find the most common value'? If you assume the parent distribution is gaussian, you could bin the data and find say the mode is the bin location with the greatest counts. However, how do you determine the bin size? Are there robust implementations available? (i.e., robust to outliers). I use `python`/`scipy`/`numpy`, but I can probably translate `R` without too much difficulty.


title: One sided confidence interval for hypothesis testing

text: I'm having a hard time understanding how to use confidence intervals for hypothesis testing. I've seen examples to reject a null hypothesis, but don't understand how to use a one-sided confidence interval to show significance. If you have two samples, x and y, how you would you use a one sided confidence interval to show $H_1:\mu_x > \mu_y$ is significant with 95%?


title: How can I perform a two-way ANOVA without replication in R?

text: Here is an example I am trying to copy from a text that I have. No idea how to proceed.               car.noise <- data.frame( speed = c("idle", "0-60mph", "over 60"), chrysler = c(41,65,76),     bmw = c(45,67,72), ford = c(44,66,76), chevy = c(45,66,77), subaru = c(46,76,64))          mcar.noise<- melt(car.noise, id.var="speed")          > mcar.noise          speed variable value     1     idle chrysler    41     2  0-60mph chrysler    65     3  over 60 chrysler    76     4     idle      bmw    45     5  0-60mph      bmw    67     6  over 60      bmw    72     7     idle     ford    44     8  0-60mph     ford    66     9  over 60     ford    76     10    idle    chevy    45     11 0-60mph    chevy    66     12 over 60    chevy    77     13    idle   subaru    46     14 0-60mph   subaru    76     15 over 60   subaru    64      I have been looking around and found anova, aov, anoava.lm and manova. One of those cases of info overload and now I am stuck. Any advice on which to use and how for a two-way anova without replication?


title: Scalability of Markov Clustering

text: I want to do graph clustering on a large dataset (A graph with 600,000 Nodes and tens of millions of edges). I read about Markov clustering. I saw this algorithm involved the calculation of a stochastic matrix and its powers. Therefore, I was wondering whether it could be used here? Otherwise, what algorithm do you suggest?


title: Two-sample bootstrap?

text: I have two independent samples of observations. From each sample I produce a statistic. Let's denote these as $\theta_1$ and $\theta_2$. I'd like to test the hypothesis that $H_0: \Theta_1=\Theta_2$, but I have these two constraints:   * There is no analytical estimate of the distributions of $\theta_1$ and $\theta_2$ (the statistic is a product of some (computationally expensive) algorithm that operates on each sample).   * Even under $H_0$, exchanging observations between the two samples in not sensible. Therefore, a permutation approach might reject $H_0$ erroneously. My current idea was to bootstrap $\theta_1$ and $\theta_2$ independently and then estimate the distribution of $\Theta_2-\Theta_1$ from these two bootstrapped distributions by means of convolution. Q1: Is this a valid approach? Q2: Any reason why not the extend this to jackknifing (instead of bootstrapping) as well? Q3: Any references to such 'two-samples' bootstrap? Q4: Any recommended alternatives?


title: Statistics for Area under the ROC curve

text: I have a question regarding statistical evaluation of the AUC. In their paper (http://www.jstor.org/stable/2531595), DeLong et al. describe a method to evaluate AUC curves. (Another good explanation can be found in the book "Statistics with Confidence: Confidence Intervals and Statistical Guidelines" by Altman et al.). As far as I understood, we compute the $\text{AUC}$ and the standard deviation $\sigma$ of the Kernel matrix. Assuming the normal distribution $\mathcal{N}(\text{AUC},\sigma)$ it is possible to compute confidence intervals. My question is about the normality assumption:   1. The $\text{AUC}$ usually lies in the interval $[0,1]$ but the interval for the normal distribtion is $(-Inf, Inf)$. Is this problem really negligible? (This problem e.g. is solved in `pROC` package by just restricting the CI to $[0,1]$)   2. The $Beta$ distribution is defined on the interval $[0,1]$ and has the shape parameters $\alpha$ and $\beta$. Can we estimate them given the data like we are able to do it for the AUC? To give an example: Given a vector `c(T,F,F,F,T,F,F,T,F,F)` the $\text{AUC} = 0.619$ and $\sigma = 0.237$ which results in 95% CI $(0.156, 1.083)$.               library(pROC)     temp.in <- c(T,F,F,F,T,F,F,T,F,F)     pROC::auc(pROC::roc(controls=which(temp.in), cases=which(!temp.in)))     pROC::ci.auc(pROC::roc(controls=which(temp.in), cases=which(!temp.in)))      Intead of using the normal distribution I would like to use the $Beta$ distribution. But how we can estimate $\alpha$ and $\beta$ for $Beta$ distribution given `c(T,F,F,F,T,F,F,T,F,F)`?


title: C parameter does not affect accuracy (LibSVM)

text: I'm working on LibSVM and I'm using linear kernel. But i have a problem with C parameter, it does not affect accuracy during test set! Is it possible? I send a picture with the distributions of the two classes and the separator line. This line did not change when I changed C. I have no idea why! :-( ![enter image description here](http://i.stack.imgur.com/B2v7T.png)


title: Selecting most important variable based on individual p-value vs. partial $R^2$

text: I'm trying to solve a problem where the goal is to find an association between children's cortisol values (`y`) against their mother's weekly cortisol averages (`x1` to `x6`) and gender (`z`). After conducting model selection strategies in R (all subset regression, backward elimination, etc.), the following two 'optimal' models were found:               y = a0 + a5*x5      and               y = a0 + a1*x1 + a3*x3 + a4*x4 + a5*x5 + a6*x6      For the second model, I found something interesting:   * x3: p-value = 0.16, partial $R^2$ = 37.1%   * x5: p-value = 0.04, partial $R^2$ = 5.5% (Let's ignore other variables for now – their p-values and partial $R^2$ values were between these two variables.) [Note: p-value refers to the probability / significance of the variable not being equal to zero; partial $R^2$ to the percent of variation that cannot be explained in a model that doesn't contain the variable] **Now to my question:** Why do I see `x5` being much more significant to the model than `x3`, yet dropping `x3` from the model will drop my $R^2$ values a lot (from around 20% to about 5%), but not so much for `x5`? Is the reason collinearity between all the variables in the model (which does exist)? Or is it something else? Also, the ultimate goal is to find the most important variable describing the response. Would I choose `x3` or `x5` in this case, and why? Or can such a choice be made?


title: Examples of costly consequences from improper use of statistical tools

text: I suspect that most users of statistical tools are ancillary users (folks who have had little to no formal training in statistics). It’s very tempting for researchers and other professionals to apply statistical methods to their data simply because they have seen it “done before” in peer-reviewed papers, grey literature, the web or at a conference. However, doing so without a clear understanding of the required assumptions and the statistical tool’s limitations can lead to erroneous results—errors often unacknowledged! I find that undergraduate students (particularly in the social and natural sciences) are either unaware of the statistical pitfalls or find these pitfalls inconsequential (the latter being most often the case). Though examples of improper use of statistical tools can be found in many introductory level text books, the web or StackExchange, I have a difficult time finding real-world examples that have had detrimental results (e.g. cost in $, lives impacted and careers lost). To that end, I am looking for real- world examples that highlight the misuse of statistical methods for which:   1. the statistical methods used are typically covered in **introductory stats courses** (i.e. inferential stats, regressions, etc…)    2. the end result have had **costly consequences** (dollars lost, lives impacted, careers shattered etc…)    3. the **data are readily available** for use as working examples in a course (the purpose is to have students work through real world examples that have had real world consequences.) One non-statistical example I like to bring up to students when discussing the importance of properly defining the units in a research project is the “metric mishap” that led to the loss of a $125M satellite! This usually invokes an :-o factor from the students and seems to have a lasting impression (at least throughout their short academic lives).


title: What is the probability of observing not more than N positive outcomes

text: I have the population which is split into three classes; A, B and C. We know that the number of observations in each class is Na, Nb and Nc respectively. Assuming that my observations are binary distributed how can I specify (with a formula) probability of observing not more than P positive outcomes.


title: Model performance in quantile modelling

text: I am using quantile regression (for example via `gbm` or `quantreg` in R) - not focusing on the median but instead an upper quantile (e.g. 75th). Coming from a predictive modeling background, I want to measure how well the model fits on a test set and be able to describe this to a business user. My question is how? In a typical setting with a continuous target I could do the following:   * Calculate the overall RMSE   * Decile the data set by the predicted value and compare the average actual to the average predicted in each decile.   * Etc. What can be done in this case, where there really is no actual value (i don't think at least) to compare the prediction to? Here is an example code:               install.packages("quantreg")     library(quantreg)          install.packages("gbm")     library(gbm)          data("barro")          trainIndx<-sample(1:nrow(barro),size=round(nrow(barro)*0.7),replace=FALSE)     train<-barro[trainIndx,]     valid<-barro[-trainIndx,]          modGBM<-gbm(y.net~., # formula                 data=train, # dataset                 distribution=list(name="quantile",alpha=0.75), # see the help for other choices                 n.trees=5000, # number of trees                 shrinkage=0.005, # shrinkage or learning rate,                 # 0.001 to 0.1 usually work                 interaction.depth=5, # 1: additive model, 2: two-way interactions, etc.                 bag.fraction = 0.5, # subsampling fraction, 0.5 is probably best                 train.fraction = 0.5, # fraction of data for training,                 # first train.fraction*N used for training                 n.minobsinnode = 10, # minimum total weight needed in each node                 cv.folds = 5, # do 3-fold cross-validation                 keep.data=TRUE, # keep a copy of the dataset with the object                 verbose=TRUE) # don’t print out progress          best.iter<-gbm.perf(modGBM,method="cv")          pred<-predict(modGBM,valid,best.iter)      Now what - since we don't observe the percentile of the conditional distribution? **Add:** I hypothesized several methods and I would like to know if they are correct and if there are better ones - also how to interpret the first:   1. Calculate the average value from the loss functions:              qregLoss<-function(actual, estimate,quantile)     {        (sum((actual-estimate)*(quantile-((actual-estimate)<0))))/length(actual)          }      This is the loss function for quantile regression - but how do we interpret the value?   2. Should we expect that if for example we are calculating the 75th percentile that on a test set, the predicted value should be greater than the actual value around 75% of the time? Are there other methods formal or heuristic to describe how well the model predicts new cases?


title: The ratio of the dependent normal rv and chi rv is?

text: (I hope this question may be too easy for some experts here). It is well-known that if $X$ is normal rv and $Y$ is a chi-square rv, then $Z=X/\sqrt{Y/n}$ follows student's $t$ distribution, where $n$ is the degrees of freedom of $Y$, and $X$ and $Y$ need to be independent. However, what comes if $X$ and $Y$ are dependent? More specifically, $X_1$ and $X_2$ are independent normal rvs, then $Z=X_1/\sqrt{X_1^2+X_2^2}$ would follow what kind of distribution?


title: Is sequential Bayesian updating an option when using MCMC?

text: I have an implementation of the Griddy Gibbs sampler, but my observations on which I'm conditioning model parameters are too many in number, thus the likelihood underflows quickly, even with a log transformation. Sequential updating holds if one is using conjugate priors, since the posterior can be found analytically, but would this be the same with MCMC, and in particular my sampler? What I'm not sure about is the loss that comes from the candidate distribution, or approximating the full conditional, which is my case. Every chunk of observation is bound to have some error, and I'm not sure if chaining chunks of observations can be trusted. In general, how can one implement a MCMC algorithm of this sort (non conjugate full conditionals for Gibbs) for say, tens of millions of observations from a model?


title: how to find out the likelihood of a model given data

text: If i have a non-stochastic model that predicts the following dataset: > [.2, .2] and the actual dataset found empirically (averaged over participants) is > [.3, .3] How would I determine the likelihood of the model? > p(Data|Model) What if the data was: > [.2, .2] Would the likelihood be 1.0? What is the procedure for calculating model likelihood?


title: Which statistical test?

text: I am designing my study, but I am a little stuck in which test I eventually should use. I have a between-subject design with 6 conditions (let's say A, B, C, D, E, F), with each having 6 responses (let's call these a, b, c, d, e, f), on a 7-point scale:               IV (condition) -> DV     A -> a, b, c, d, e, f     B -> a, b, c, d, e, f     C -> a, b, c, d, e, f     D -> a, b, c, d, e, f     E -> a, b, c, d, e, f     F -> a, b, c, d, e, f      So, each condition has 6 dependent variables. I don't want to test difference between the conditions, but only the effect of the condition on the different DV's. For example, when someone is in condition A, I would expect them to score higher on DV a than on DV b, c, d, e, and f. When someone is in condition B, I would expect them to score higher on another DV, and so forth. Later on, I expect also some moderator effects, but I guess that's only adding interaction terms to the model. Can I use separate tests for each condition? Or should I test it in one big model? What kind of test would be most suitable for this. Most tables that should help you choose a statistical test only show options with 1 DV.


title: Is it OK to do additive smoothing before applying Pearson's chi-square test for independence?

text: I'm concerned about treating my data as gold, especially in areas of low data support, so I would like to apply additive smoothing. I'm then doing several things with this data, and one of them is Pearson's chi-square test for independence. Is it still acceptable to first do the smoothing and then do the chi-square test? The former is very Bayesian and the latter is very frequentest, so at the least it feels somewhat strange. However, it seems like smoothing should curtail the need to worry much about the data support or the need for something like Fisher's exact test. My reasoning for wanting to use competing philosophies is because the audience really wants to see p-values, but the Bayesian inside of me is fighting to be let loose. As such I came up with this compromise. EDIT: To clarify how I envision the smoothing taking place, consider the following contingency table:                   X=0  X=1        +-------+     Y=0|  2|  2|        +---+---+     Y=1|  3|  3|        +---+---+      So without smoothing we have: $P(X=0)=\frac{5}{10}$, $P(X=1)=\frac{5}{10}$, $P(Y=0)=\frac{4}{10}$, $P(Y=1)=\frac{6}{10}$ With additive smoothing and say $\alpha=1$ we get $P(X=0)=\frac{6}{12}$, $P(X=1)=\frac{6}{12}$, $P(Y=0)=\frac{5}{12}$, $P(Y=1)=\frac{7}{12}$ I haven't yet decided how to handle the contingency table. One way would be to smooth each cell:                   X=0  X=1        +-------+     Y=0|  3|  3|        +---+---+     Y=1|  4|  4|        +---+---+


title: Finding precision of Monte Carlo simulation estimate

text: **Background** I am designing a Monte Carlo simulation that combines the outputs of series of models, and I want to be sure that the simulation will allow me to make reasonable claims about the probability of the simulated outcome and the precision of that probability estimate. The simulation will find the probability that a jury drawn from a specified community will convict a certain defendant. These are the steps of the simulation:   1. Using existing data, generate a logistic probability model ( _M_ ) by regressing “juror first ballot vote” on demographic predictors.   2. Use Monte Carlo methods to simulate 1,000 versions of _M_ (i.e., 1000 versions of the coefficients for the model parameters).   3. Select one of the 1,000 versions of the model ( _M i_).   4. Empanel 1,000 juries by randomly selecting 1,000 sets of 12 “jurors” from a “community” ( _C_ ) of individuals with specified demographic characteristic distributions.   5. Deterministically calculate the probability of a first ballot guilty vote for each juror using _M i_.   6. Render each "juror’s" probable vote into a determinate vote (based on whether it is greater or less than randomly selected value between 0-1).   7. Determine each "jury’s" “final vote” by using a model (derived from empirical data) of the probability a jury will convict, conditional on the proportion of jurors voting for conviction on the first ballot.   8. Store the proportion of guilty verdicts for the 1000 juries ( _PG i_).   9. Repeat steps 3-8 for each of the 1,000 simulated versions of _M_.   10. Calculate the mean value of _PG_ and report that as the point estimate of the probability of conviction in _C_.   11. Identify the 2.5 & 97.5 percentile values for _PG_ and report that as 0.95 confidence interval. I am currently using 1,000 jurors and 1,000 juries on the theory that 1,000 random draws from a probability distribution--demographic characteristics of _C_ or versions of _M_ --will fill out that distribution. **Questions** Will this allow me to accurately determine the precision of my estimate? If so, how many juries do I need to empanel for each _PG i_ calculation to cover _C_ 's probability distribution (so I avoid selection bias); may I use fewer than 1,000? Thank you so much for any help!


title: Meaning and use of `reltol` in `nnet` library in R

text: I'm trying to use the `nnet` library in R, and can't seem to work out how to use the `reltol` parameter. It says in the docs: > Stop if the optimizer is unable to reduce the fit criterion by a factor of > at least 1 - reltol. I assume this means that if `reltol = 0.5` then it will stop if it can't halve the error rate, and if `reltol = 0.1` it will stop if it can't reduce the error rate by 10%. Is that correct? The experiments that I've done don't seem to work like that...


title: How to check reliability of a given sampling technology?

text: Very basic question, I suspect. It's sort of the inverse of this one, and the same disclaimers about my own ignorance apply, but here goes: I have a sampling technology. I know the total population, and I know the number of samples I've done. What's the best way to calculate and express how confident I am that my results are correct? A few more points in case they make a difference:   * My test checks whether an item is in one of three states (and they are the only possible states).   * I'm assuming that the test is truly a random distribution.   * In general, my total population is > 1,000,000 and my sample is ~20% of that, but it can vary quite a bit. And a wrinkle (I hesitate to even put this in because I want to keep things simple):   * There's a slight error where my reported Total Population value might be slightly smaller than the actual total (by 5-10%). I can get the actual total if need be, but won't bother if it doesn't really make a difference. Rough and simple / clear is probably better than more accurate if there are different ways of expressing this. I.e., something like "we're 95% confident that this is correct" is probably better than "we're 95% confident that this is correct to within 3%". (Am I talking about P?) **UPDATE:** Hmm... I think my "wrinkle" may have introduced confusion. I'm not trying to find out how accurate the total is. I can get the actual total, but I have to use a more convoluted method. So my main question was about saying if I have 10,000 samples out of 1,000,000 items how confident can I be that my sampled distribution among the three states is correct. My secondary question (which is where the wrinkle comes in) is: given that the actual total is 5-10% larger than the one I'm using, how much of a difference does that make to my confidence? In other words, I'm saying 1,000,000 items, but it might actually be 1,100,000. Should I bother to go through the convoluted process to get the actual total, just to compute the confidence? To me, that seems unlikely to make a significant difference to my confidence level, but I thought I should check.


title: How to get standard error of a function (delta method vs. simulation)?

text: I want to get a confidence interval of a function of some parameters. for example, from the data I estimate parameters of Pareto. Now I want to get 95% CI for 90th quantile (it's a function of parameters of Pareto), so I would need standard error. I know delta method is one option. For simulation method, I am wondering if it is legitimate to simulate 1000 samples of size 50 from Pareto, calculate each of the 90th quantiles and take the standard deviation of the 1000 data. Is the standard deviation I get equivalent to standard error? Thanks for your help!


title: Differences between heavy tail and fat tail distributions

text: I thought heavy tail = fat tail, but some articles I read gave me a sense that they aren't. One of them says: heavy tail means the distribution have infinite jth moment for some integer j. Additionally all the dfs in the pot-domain of attraction of a Pareto df are heavy-tailed. If the density has a high central peak and long tails, then the kurtosis is typically large. A df with kurtosis larger than 3 is fat-tailed or leptokurtic. I still don't have a concrete distinction between these two (heavy tail vs. fat tail). Any thoughts or pointers to relevant articles would be appreciated.


title: Are there arithmetical tests analogues for graphical test?

text: I learned in statistics basic course that sometimes one can make reasoning about the data by making graphical plots, like testing outliers or P-P-plot or Q-Q-plot. But are there some test which measures the same things by numbers? Like what can one do if he or she have to analyze some data in a computer that don't have graphical desktop installed?


title: Interpreting coefficients of a dlog-dlog regression

text: I am running a dlog-dlog (difference of logarithm*) regression and I want to convert the coefficients into marginal effects. I know that it's different from a log-log regression, in which the coefficients directly give us the elasticities. How can we interpret the coefficients from a dlog-dlog regression? * For example dlog (Y)= a + b dlog(X)+ error term.


title: Test data distribution different than training data distribution

text: What could be some issues if the distribution of the test data is significantly different than the distribution of the training data and why is that? ## Edit 1 ## The training and testing sets came from different processes. I was thinking if I use a set of data, that is say normally distributed, to train and build a model (e.g. logistic regression) and then use that model to make predictions on a new set of data, say it is right skewed. What could be the issues here and how would it impact the predictions?


title: Help with R code

text: I am a newbie at R and I'm stuck at the following code               $y=c(2.5,6.0,6.0,7.5,8.0,8.0,16.0,6.0,5.0,6.0,28.0,5.0,9.5,6.0,4.5,10.0,14.0,3.0,4.5,5.5,3.0,3.5,6.0,2.0,3.0,4.0,6.0,5.0,6.5,5.0,10.0,6.0,18.0,4.5,20.0)     x2=c(650,2500,900,800,3070,2866,7500,800,800,650,2100,2000,2200,500,1500,3000,2200,350,1000,600,300,1500,2200,900,600,2000,800,950,1750,500,4400,600,5200,850,5000)     x1=c(16.083,48.350,33.650,45.600,62.267,73.2170,204.617,36.367,29.750,39.7500,192.667,43.050,65.000,44.133,26.933,72.250,98.417,78.650,17.417,32.567,15.950,27.900,47.633,17.933,18.683,26.217,34.433,28.567,50.500,20.950,85.583,32.3830,170.250,28.1000,159.8330)     library(MASS)      n=length(y)     X=matrix(nrow=n,ncol=2)     for(i in 1:n){X[i,1]=x1[i]}     for(i in 1:n){X[i,2]=x2[i]}     nbrun=0     gibbs=function(data,m01=0,m02=0,k01=0.1,k02=0.1,a0=0.1,L0=0.1,nburn=0,ndraw=5000)     {m0=c(m01,m02)     C0=matrix(nrow=2,ncol=2)     C0[1,1]=1/k01     C0[1,2]=0     C0[2,1]=0     C0[2,2]=1/k02     beta=mvrnorm(1,m0,C0)     omega=rgamma(1,a0,1)/L0     draws=matrix(ncol=3,nrow=ndraw)     while(-nburn < ndraw){ -nburn=-nburn+1     C1=solve(solve(C0)+omega*t(X)%*%X)     m1=C1%*%(solve(C0)%*%m0+omega*t(X)%*%y)     beta=mvrnorm(1,m1,C1)     a1=a0+n/2     L1=L0+t(y-X%*%beta)%*%(y-X%*%beta)/2     omega=rgamma(1,a1,1)/L1     if(-nburn>0){     draws[-nburn,1]=beta[1]     draws[-nburn,2]=beta[2]     draws[-nburn,3]=omega     }     }     return(draws)}$      whenever I run the code I get:"Error: object 'nburn' not found" and I don't understand why


title: Binary linear SVM in terms of structured SVM

text: On slide 32/85 of this tutorial on structured SVM learning, the author formulates binary SVM classification in terms of the structured SVM just introduced. Why is there a constant factor of $1/2$ in the feature map and inference? I can't see it making any difference to the result.


title: Statistical Analysis on mostly boolean values

text: So I have a large dataset, and I was wondering what the best way to conduct statistical analysis of it is. I'm very green in terms of statistical methods, but I learn quickly. Basically, each item has a couple attributes, and each attribute has several possibilities. Each item has their specific attribute set-up in terms of booleans for each possibility of each attribute (i.e. 1 means it has that possible version of that attribute, 0 means it doesn't). For some of the attributes, each item can only have one possible value, and for some it can have multiple. I'm using R through RStudio to conduct my analysis. Any help would be appreciated! Basically, I have a bunch of error messages generated by different customers. Each error has different attributes, which include server downtime, country of origin, the type of change, which specific item(s) caused the problem, and a couple others. I'm mostly trying to look for which applications/types of changes prompt downtime.


title: Knapsack problem with uncertain profits

text: I am trying to find the optimal strategy for a game where the goal is to pick a number of items where the profit is uncertain but the weights are set. Each round of picking is a specific point in time, and each round is partially independent from the past (think sports games). This is basically a version of the knapsack problem. (see Knapsack Problem) Here is the problem :   Maximize: $$\sum _{j=1}^{N}\sum _{j=1}^{\left | G_{i} \right |}p_{ij}x_{ij} +\sum _{j=1}^{\left | G_{s} \right |}p_{sj}x_{sj}$$ Subject to: $$ \begin{align} (\sum _{i=1}^{n}\sum _{j=1}^{\left | G_{i} \right |}w_{ij}x_{ij} +\sum _{j=1}^{\left | G_{s} \right |}w_{sj}x_{sj})&\leq c \\ \sum _{j=1}^{\left | G_{i} \right |}x_{ij} &= 1 \\ \sum _{j=1}^{\left | G_{s} \right |}x_{sj} &= k \\ x_{ij} \in \{0,1\}\\ x_{sj} \in \{0,1\} \end{align} $$ Notation:   _N_ the number of items   _n_ the number of groups   _c_ the capacity of a single knapsack   _G_ = {$G_{1},...,G_{n}$} the set of groups   ${\left | G_{i} \right |}$ the number of items in group $G_{i}$   ${\left | G_{s} \right |}$ the number of items in group $G_{s}$ (special group)   $p_{ij}$ the profit of item _j_ of group $G_{i}$   $w_{ij}$ the weight of item _j_ of group $G_{i}$   $p_{sj}$ the profit of item _j_ of group $G_{s}$   $w_{sj}$ the weight of item _j_ of group $G_{s}$   _k_ restraint of group $G_{s}$   So for example, say you have 100 items per group where each item has a profit and a weight. In this particular instance, the profit is uncertain and therefore you must use regression to come up with a reasonable estimate. So far, I am using linear regression paired with a k-NN algorithm to predict the profit of each item based on past observations. I also use a decision tree learner for a classification between items which have either a positive profit or a negative profit (again using past observations). I'm going to use GLPK to solve the knapsack problem. In fact, I've already coded the module and its ready to go. The problem is with the uncertainty of the profits. I'm having a tough time conceptualizing uncertainty and risk. Is there a way to segregate item choices according to risk? Perhaps I can look at the historical variance in profit value for each item? Is there any particular method that would be useful here? (Something akin to Monte Carlo, etc?)


title: Why does "explaining away" make intuitive sense?

text: I recently learned about a principle of probabilistic reasoning called "explaining away," and I am trying to grasp an intuition for it. Let me set up a scenario. Let $A$ be the event that an earthquake is occurring. Let event $B$ be the event that the jolly green giant is strolling around town. Let $C$ be the event that the ground is shaking. Let $A \perp\!\!\!\perp B$. As you see, either $A$ or $B$ can cause $C$. I use "explain away" reasoning, if $C$ occurs, one of $P(A)$ or $P(B)$ increases, but the other decreases since I don't need alternative reasons to explain why $C$ occurred. However, my current intuition tells me that both $P(A)$ and $P(B)$ should increase if $C$ occurs since $C$ occurring makes it more likely that any of the causes for $C$ occurred. How do I reconcile my current intuition with the idea of explaining away? How do I use explaining away to justify that $A$ and $B$ are conditionally dependent on $C$?


title: How to interpret log-log regression coefficients for other than 1 or 10 percent change?

text: I have read many threads here on how to interpret coefficients in a regression where the predictor and the dependent variable are log-transformed. Most give an answer for a one or ten percent change. However, I am not sure whether I am right about interpreting other changes. I have a beta coefficient of -0.5057. Given the standard interpretation (as elasticity) a one percent increase in the predictor variable (X) leads to an decrease of 0.5057 percent in the dependent variable (Y), alternatively a ten percent increase in X leads to a 5.057 percent decrease in Y. I found the formulas 1.1^(-0.5057) and exp(-0.5057(ln(1.1))) to calculate the affect of any change in X on Y. Are these formulas right because I get for a ten percent increase slightly different results than the 5.057 decrease. Furthermore, is it right that the percentage values for the changes in Y do not increase proportianally with increases in X? If a one percent increase in X leads to an decrease of 0.5057 percent in Y, it is wrong to assume that a 50 percent increase in X leads to (50 x 0.5057) decrease in Y? exp(-0.5057(ln(1.5)))= 0.8146 Given the difference between 1 and this value, it suggests that Y decreases by 18.5 percent. Am I right? The formulas mentioned above suggest so. Is that because of the nonlinearity in the model?


title: Calculating principal component scores after PC analysis

text: I am carrying out a study to find out meteorological patterns using daily met observations including around 30 met parameters (each day is a case with 30 variables). My methodology includes carrying out a PCA:   1. To reduce 30 variables to smaller number of PCs.   2. Find out the PC scores for all days in the ten years. As first I standardized all 30 variables using SPSS function `ANALYSE>>> DESCRIPTIVE`. Now that I have got 6 PCs explaining 80% of variance of original data, I have to calculate the PC scores. My questions are as under;   1. Should I use Component Matrix OR Rotated Component Matrix for PC score calculation?   2. Should I multiply the component/Rotated component matrix with original variable matrix (un-standardized) OR i should multiply it with standardized form of original variable matrix?


title: Good book about theoretical approach to statistics

text: When I took courses in theoretical statistics as an undergrad 10 years ago, we used _Modern Mathematical Statistics_ by Dudewicz and Mishra. I find myself referring back to the book now and am reminded some of the code examples are in assembly for an IBM 370. While quaint, I cannot help but feel this is somewhat dated. What high quality books exist of more recent vintage?


title: What value is conditional logistic regression if two cohorts are already matched on everything of interest?

text: My situation is this:   * I am comparing two cohorts   * I have matched the two cohorts on all of the factors that I am interested in   * The two cohorts are already balanced on all of the factors that I would include in a regression Questions:   * Is it appropriate to just report the univariate results?   * Would conditional logistic regression add anything to/change the results?


title: Error using rfe in caret package in R

text: I am doing some exploratory data analysis in the Heritage Health Prize , and have come across a weird error using R's caret package. In the dataset, I've created a dataframe counting how many times a patient (by `MemberID`) has visited a specific primary care physician (denoted by a number) in year 1 insurance claims. Let's say there are 57 such primary care physicians, and so our dataframe looks like this:               > head(cbind(right.a[,1],right.a[,107:111]))           MemberID 2448 91972 20893 30870 13281     1     24027523    0     0     0     0     0     2     98321677    2     0     0     3     0     3     33896267    0     1     0     0     1     4      5481382    0     0     2     0     0     5     69908341    0     0     0     0    17     6      8169352    0     0     0     0     0      As you can see, this is a very sparse matrix. Our outcome here is the number of days they spent in hospital the following year in year 2 (`DaysInHospital`). I'm using year 2 insurance claims and year 3 days in hospital as a validation set. I'm using rfe to do some recursive feature selection, and I'm getting the following error:               > subsets <- seq(12,57,by = 3)     > lmProfile <- rfe(right.a[,c(107:164)],                  log1p(right.a$DaysInHospital),                  right.b[,c(107:164)],                  log1p(right.b$DaysInHospital),                  sizes=subsets,                  metric="RMSE",                  maximize="FALSE",                  rfeControl=ctrl)     .     (snip)     .     External resampling iter:        10      Computing importance     Fitting subset size:   58      Computing importance     Computing importance     Error in y$pred : $ operator is invalid for atomic vectors      Am I missing something here?


title: What is the difference between hazard and crude ratio?

text: I've recently calculated a hazard ratio doing a multi variable Cox regression analysis (Kaplan Meier) for an assignment. Part of my feedback was to adjust the hazard ratio to a 'crude ratio'. How do I do that? I'm using SPSS 19/PASW.


title: Temporal correlation of a discrete signal

text: I have a physical signal (for example, `Temperature`, `Humidity`) which is sampled at discrete times. I wanted to know how to compute the temporal correlation between "a value of the signal" and "another group of values in the same signal, but acquired previously". I found this formula on the topic but I don't understand it at all. Another question linked to this precedent is to know if using of the temporal correlation could help me to get a simple and practical solution to my problem which is: I wanted to detect in real time if the new acquired sample from the physical phenomenon (e.g., `Temperature`) presents a variation (strong or weak) in comparison with the previously acquired samples. My primary goal is to detect in real time when there is variation in the signal and then increase the sampling frequency to have more samples and so have a higher resolution of the signal. This known as "Adaptive Sampling".


title: What is the best forecasting model for time series data?

text: I have a time series data (1 minute and sometimes 5 minute data) data I would like use `forecasting` package to forecast couple hours ahead. Here is my data:               dput(head(p,20))     structure(list(time = structure(c(1373889420, 1373889480, 1373889540,      1373889600, 1373889660, 1373889720, 1373889780, 1373889840, 1373889900,      1373889960, 1373890020, 1373890080, 1373890140, 1373890200, 1373890260,      1373890320, 1373890380, 1373890440, 1373890500, 1373890560), class = c("POSIXct",      "POSIXt"), tzone = "America/New_York"), cpu = c(2.25892, 2.04144,      5.04823333333333, 4.9947, 1.72982857142857, 4.82655, 3.6168625,      4.7357, 2.42683333333333, 3.62635, 5.02315714285714, 2.57147142857143,      7.16005, 2.34253333333333, 2.82315714285714, 5.17668, 2.2899375,      6.92, 5.172375, 4.63735), name = c("servers", "servers", "servers",      "servers", "servers", "servers", "servers", "servers", "servers",      "servers", "servers", "servers", "servers", "servers", "servers",      "servers", "servers", "servers", "servers", "servers")), .Names = c("time",      "cpu", "name"), row.names = c(1116L, 1411L, 123L, 226L, 1014L,      435L, 538L, 569L, 1081L, 342L, 74L, 865L, 178L, 890L, 281L, 166L,      1035L, 143L, 112L, 91L), class = "data.frame")          x.xts <- xts(p$cpu, p$time)     x.ts <- as.ts(x.xts)     x.ets <- ets(x.ts)     x.fore <- forecast(x.ets, h=120)     f<-data.frame(x.fore$mean)         DateTime<-tail(z,1)$time     f$DATE <- DateTime + 60 * (seq_len(nrow(f))-1)         colnames(f)<-c("cpu", "time")         f$name<-c("forecast")      I see that `cpu` is the same for all future data times:                cpu                time     name     1 6.020207 2013-07-15 11:57:00 forecast     2 6.020207 2013-07-15 11:58:00 forecast     3 6.020207 2013-07-15 11:59:00 forecast     4 6.020207 2013-07-15 12:00:00 forecast     5 6.020207 2013-07-15 12:01:00 forecast     6 6.020207 2013-07-15 12:02:00 forecast      Is there a better forecasting model besides `ets` for time series data?


title: relationship between ARMA and AR

text: I once heard some statements regarding the relationship between ARMA and AR process, such as               An average of severl lags of an autoregression forms an ARMA process      or               A weighted mixture of lags of an AR(P) model is ARMA      I am not very clear about these two statements, whether they are right? If they are right, how to get them? Thanks.


title: Calculating when significance will be reached in A/B test ($\chi^2$ test of the difference between two proportions)

text: I'm a marketer trying to understand the math involved in A/B testing. I've come across a tool that estimates how many more samples are required to reach significance for a test that's already begun. When you enter your test results thus far, it tells you how many more samples are required to reach the 90% confidence level. If I enter the following data into the form:                          Visitors    Goals     Group A:        200       25     Group B:        190       26      The tool tells me that I need an estimated 8405 more samples to reach significance at the 90% confidence level. So far I've:   * Read through the explanation of sample size in a couple of different statistics books.   * Manipulation of power analysis equations on paper.   * Experimenting with `power.prop.test()` in `R`. I'm just not seeing how to calculate this. Can someone help?


title: Online feature preprocessing

text: Say we have a batch of data we use to learn an SVM, a Regressor, etc. And we use the learned model to do online classification, the learning is not done online, the online part is only the processing. How do we go about Normalization, Whitening, Rescaling, etc Since most of these techniques require a batch of data to work, and even then, if the features do not have the same normalization, it may not work as advertised. Thanks


title: How to fit ARMA(4,4) model, but some coefficients set to zero / fixed?

text: Suppose I have the following ACF and PACF (data: ![ap](http://i.stack.imgur.com/A7B44.png) I want to fit an ARMA-GARCH process. Currently I want to do the first step, specify the mean equation. The first model just uses a constant $\mu$, so no ARMA. In the second model I was thinking about a modified ARMA(1,1) or ARMA(4,4), I don't know what this is called. I want to only use the 4th lag order in the AR and MA part. So this is basically an ARMA(4,4) where the coefficients of the first three lags are set to zero. $r_t=\delta + \epsilon_t + \alpha_4 r_{t-4}+ b_4 \epsilon_{t-4}$ How can I fit this model in R? I tried                  arima(logloss, order=c(4,0,4),fixed=c(0,0,0,NA,0,0,0,NA,NA))      First of all: Is this correct? Second: Does this make sense? My output is the following: ![ts](http://i.stack.imgur.com/nfi5A.png) If I calculate the p-values via               # p-values     (1-pnorm(abs(aa$coef)/sqrt(diag(aa$var.coef))))/2      I get               > (1-pnorm(abs(aa$coef)/sqrt(diag(aa$var.coef))))/2              ar1          ar2          ar3          ar4          ma1          ma2      2.500000e-01 2.500000e-01 2.500000e-01 4.431378e-08 2.500000e-01 2.500000e-01               ma3          ma4    intercept      2.500000e-01 2.523225e-06 1.886732e-01      >       So can I say, that the both coefficients of the 4th lag order are highly significant, but the intercept is not significant, correct? So should I also fix it to zero? If I just fit a model with a mean, so no AR or MA, I get: ![ts2](http://i.stack.imgur.com/Mr8Gp.png) So the mean is also not significant. What should I do? Fit a GARCH without a mean equation? So no mean, no AR or MA part? EDIT: I played around with it and I found, that an ARIMA(5,0,5) with the first 3 lags fixed to zero and the mean fixed to zero seems to be approrpriate. The output is: ![ts4](http://i.stack.imgur.com/VxuTJ.png) The AIC is smaller than in case of the ARIMA(4,0,4) with mean fixed to zero and the residuals look ok. Are my model building steps correct?


title: R Packages for Panel GARCH?

text: Are there any packages that let me estimate panel GARCH models in R? I have looked extensively in Google but have not found anything.


title: Best method to analyse survey results with multiple choice of answers

text: I have designed a survey with multiple choice answers. Each question contains same set of answers   * Strongly agree    * Agree    * Disagree    * Strongly disagree   * Don't know There are 25 questions and questions are grouped into 5 areas. I need to analyze each area. What is the best approach to analyze the survey results which I collected from a number of people? Is it a good method to calculate the mean response from each area as I need area-wise analysis of questions? Can any one suggest a best option? I am not a statistician Thanks in advance


title: Understanding output of a GEE model in R - Naive S.E vs. Robust S.E?

text: Happy New Year to all of you. I have a question regarding the interpretation of the output I get by fitting a GEE model in R. Attached is the picture of the output I get: ![enter image description here](http://i.stack.imgur.com/3LhaW.png) 1) The first column, Estimate, is obviously the estimated $\beta$'s under the GEE model. 3) The third column seems to be the Z statistic calculated using the Naive S.E 5) Fifth column seems to be the Z statistic calculated using the Robust S.E 2) **Naive S.E** : I guess this one is the estimated S.E using the correlation structure I used in my GEE (here: corstr = "exchangeable") 4) **Robust S.E** : I guess is the S.E estimated using the empirical sandwich estimator. Is that right?


title: Chi-square test in finite mixture models

text: I was trying to figure out whether or not the distribution of a biomarker came from heterogeneous populations. Analyzing the data with `normalmixEM` in the `mixtools` package in R, I got two distributions, but the histogram was unimodal and seemed to be homogeneous. (Sorry, I cannot post an image because I'm new here.) So I conducted a chi-square test, `test.equality(y = biomarker)`, to see if it was from a homogeneous population but it gave a significant p-value following a warning of non-convergence. So my two related questions are:   1. How do you normally figure out if the distribution is two-component instead of just one?    2. Even if the analysis gives you two components, can the distribution still be from one population?


title: Confidence bands in case of fitting ARIMA in R?

text: I want to look at the acf and pacf of my data, to identify the model for my mean equation, so I want to fit an ARMA for my mean equation and later on model the conditional variance by a ARCH/GARCH (I know I have to do jointly model estimation). In the first step I want to look at the ACF and PACF for identifying, I used the standard Acf of the forecast library, but I noticed, that the confidence bands in these plots are given for testing randomness and not for fitting a ARMA. As wikipedia says: > Correlograms are also used in the model identification stage for fitting > ARIMA models. In this case, a moving average model is assumed for the data > and the following confidence bands should be generated: $ \pm z_{1-\alpha/2}\sqrt{\frac{1}{N}\left(1+2\sum_{i=1}^{k} y_i^2\right)} $ How can I get these confidence bands, which increase as the lag increases?


title: Why doesn't classifier accuracy increase linearly with k in k-fold validation?

text: I am experimenting with WEKA and I'm trying to understand the impact of k in cross-fold validation. It seems reasonable to me that the higher k is, the more accurate the classifier will be (with respect to the train/test data, that is). This doesn't seem to be the case when I am varying that parameter. Why not?


title: Comparing two different leagues of similar but not equal distributions around a standard deviation of error of a prediction from a rating system

text: This query ties a lot of my interests in rating sports teams together, because as I’ve mentioned before I do a version of the Kenneth Massey method (as per his 1997 thesis http://masseyratings.com/theory/massey97.pdf) of rating Sports Teams. I do it in Excel and usually do it by setting up a least squares based problem and either have the Solver application solve it, or setting up a conditional recursive loop (as this will provide two extra elements - average margin of victory and strength of schedule). This enables me to make projections of a predicted margin of victory in a match (by having the ratings for two teams and a Home Field/Court advantage). Because the standard deviation of the error of the prediction (against the actual score) is usually normally distributed (as per Stern “On the Probability of Winning a Football Game” American Statistician 45 (1991): 179-83) which can be virtually proved via a QQ plot amongst other things - it allows you to put percentages on projections (by taking the predicated margin of victory between the two teams from the ratings as a mean and the standard deviation as per the above). My query is this, just before Christmas in the Spanish Asbol Handball League Barcelona were 15 wins with 0 losses and a derived rating of 10.88. The standard deviation of the error about the prediction was 4.1. In the Russian Superleague, Chehovski were also 15 wins with 0 losses and a derived rating of 13.97. The standard deviation of the error about the prediction was 4.4. Who would be the better team based on the approach above? Can I or is there a way to compare two teams effectively from two different distributions? By this I am assuming that it would be unsound to assume that Chehovski is 3.09 points better on a neutral ground than Barcelona (e.g. 13.97 – 10.88). Is there a way to make a correction to adjust for the fact that the two leagues have a different distribution (I haven’t looked at these leagues over a couple of years yet, which I would normally do to get a more accurate standard deviation so currently I am assuming the standard deviations are different for the two different leagues and the 4.1 and 4.4 are approximately reflective after half a season). In this case I am just keen to figure out who is effectively the best handball team statistically (and how you can effectively separate two teams in different leagues with the same record). Kenneth Massey in the ratings part of his site doesn’t treat leagues individually (and also refined his rating formula from his thesis so that it is bias free and not affected by points scored in the same way as Dr Wesley Colley’s approach http://www.colleyrankings.com/method.html). I don’t have a great maths background (I just enjoy sports) so I would be keen for any nudges in the right direction as I'd like to find a mathematically sound way of doing this. Thanks,


title: Student T CDF - How to read the hypergeometric function formula?

text: I am trying to program out a couple of statistical formulas with python. I am not concerned with efficiency or the fact that there already are tools to calculate the values, but I am doing it to increase my understanding of statistics. I tried setting up the Student T CDF, but I have no clue how to even start programming this. On wikipedia it states the studentt cdf as follows: $$ \frac{1}{2}+x\Gamma\left(\frac{\nu+1}{2}\right)\cdot\frac{_2F_1\left(\frac{1}{2},\frac{\nu+1}{2},\frac{3}{2},-\frac{x^2}{\nu}\right)}{\sqrt{\pi\nu}\Gamma\left(\frac{\nu}{2}\right)} $$ It is noted that 2F1 is the Hypergeometric Function I have no idea what it does, how it works or what it looks like. I can't figure out how to program this function or even read it, since there must be something mysterious about the 2F1 notion that I don't understand. I tried reading up on the topic but I can only find very thick and old books that make me feel very stupid. The wikipedia link does not really help me, since I don't understand it. I don't have a good understanding of function theory and have never encountered gamma, beta and so on. My desire might be too reductive and limited: I think of a function as something that calculates a value for a given input, so f(x) = 2x can pre easily programmed and will return 2 for f(1), 4 for f(2) and so on. So my question is: **How do I even read this formula?** **Update** I am sorry to admit it, but I am really dependent on your patience here. I have really big problems understanding the formula's "code" to begin with. I feel that until I do, I will have no use for your answers - so please help me on a lower level, as if you were explaining it to your grandma. 2F1 - what is the 2 for, what is the F for, what is the 1 for? Where does the (q)n come from? 2F1(a,b;c;z)? Why are a and b separated with a comma, not a semicolon? How is 2F1 different from 3F8, or is there even such a thing? This is all hieroglyphs to me. How do I read it? My question is more based on a reading issue than on a logical or a lack of understanding. I feel like translating chinese without a dictionary!


title: How to assess whether a coin tossed 900 times and comes up heads 490 times is biased?

text: A coin is tossed 900 times and heads appeared 490 times. Does the result support the hypothesis that the coin is unbiased?


title: Neg Binomial and the Jeffreys' Prior

text: I'm trying to obtain the Jeffreys' prior for a negative binomial distribution. I can't see where I go wrong, so if someone could help point that out that would be appreciated. Okay, so the situation is this: I am to compare the prior distributions obtained using a binomial and a negative binomial, where (in both cases) there are $n$ trials and $m$ successes. I get the right answer for the binomial case, but not for the negative binomial. Let's call the Jeffreys' prior $\pi_J(\theta)$. Then, $$ \pi_J(\theta)\propto [I(\theta)]^{1/2}. $$ Under the regularity conditions (fulfilled as we are dealing with the exponential family), $$ I(\theta)=-E\left(\frac{\partial^2 \log L(\theta|x)}{\partial \theta^2}\right) $$ where for the negative binomial $n$ is $x$ in the above expression (the total number of successes $m$ is fixed, $n$ is not). The distribution --I think-- is $$ p(m|\theta)\propto\theta^m(1-\theta)^{n-m} $$ since $\theta$ is defined as the probability of success and $m$ is the number of successes. This is also the likelihood, since $m$ is a scalar and not a vector. Hence, $$ L(\theta|n)\propto\theta^m(1-\theta)^{n-m}\\ \log L(\theta|n)=m\log\theta +(n-m)\log (1-\theta)\\ \frac{\partial\log L(\theta|n)}{\partial \theta}=\frac{m}{\theta}-\frac{n-m}{1-\theta}\\ \frac{\partial^2\log L(\theta|n)}{\partial \theta^2}=-\frac{m}{\theta^2}-\frac{n-m}{(1-\theta)^2} $$ so the Fisher information is $$ I(\theta)=-E\left(\frac{\partial^2\log L(\theta|n)}{\partial \theta^2}\right)=\frac{m}{\theta^2}+\frac{E(n)-m}{(1-\theta)^2}=\frac{m}{\theta^2}+\frac{\frac{m\theta}{1-\theta}-m}{(1-\theta)^2}\\ =\frac{m(1-\theta)^2+\frac{m\theta^3}{(1-\theta)}-m\theta^2}{\theta^2(1-\theta)^2}=\frac{m(1-2\theta)+\frac{m\theta^3}{(1-\theta)}}{\theta^2(1-\theta)^2}\\ =\frac{m(1-2\theta)(1-\theta)+m\theta^3}{\theta^2(1-\theta)^3}=\frac{m(1-3\theta+2\theta^2+\theta^3)}{\theta^2(1-\theta)^3}\\ \propto\frac{1-3\theta+2\theta^2+\theta^3}{\theta^2(1-\theta)^3} $$ This, however, does not give me the correct answer. The correct answer is $$ \pi_J(\theta)\propto \frac{1}{\theta(1-\theta)^{1/2}} $$ which means that the Information I get should be $$ I(\theta)=\frac{1}{\theta^2(1-\theta)} $$ since the prior should be proportional to the square root of the information. Can anyone find any mistakes? I wouldn't be surprised if I screwed something up with the set up of the distribution (successes vs failures with their respective probabilities, etc). I used the expected value from Wikipedia and I know the correct answer from here (page 3).


title: Basic descriptive statistics on an irregular data set

text: I have not used my stats in a while and was interested in how the group would handle the following: I have a multi-year insect data set that has irregularities in trapping effort, trap type and area of trap size within a specific trap type. I am interested in obtaining measures of relative abundance, dominant species, species diversity and richness numbers. I was considering using a basic ANOVA and Shannon-Weiner Diversity index and looking at the data based on hours and area of trapping effort. Any thoughts...


title: Kolmogorov-Smirnov two-sample $p$-values

text: I am using the Kolmogorov-Smirnov two-sample test to compare distributions, and I noticed a $p$-value is frequently reported as the test statistic. How is this $p$-value determined? I know it's the probability of obtaining a result at least as large as the one obtained, but how is this $p$-value determined given this is a nonparametric test? That is, we can't assume Gaussian fluctuations in the distribution and compute the $p$-value using a $t$-test. Thanks!


title: C++ for statistics

text: I'd like to get into C++ for solving statistical problems. At the moment I should learn the basics, in particular things related to matrix computation (svd decomposition and the like), numerical derivatives, and so on. I hope somebody out there is willing to suggest source materials (both books and websites).


title: Smooth a circular/periodic time series

text: I have data for motor vehicle crashes by hour of the day. As you would expect, they are high in the middle of the day and peak at rush-hour. ggplot2's default geom_density smooths it out nicely A subset of the data, for drink-drive-related crashes, is high at either end of the day (evenings and early mornings) and highest at the extremes. But ggplot2's default geom_density still dips at the right-hand extreme. What to do about this? The aim is merely visualisation -- no need (is there?) for robust statistical analysis. ![Imgur](http://i.imgur.com/MyKHHVS.png)               x <- structure(list(hour = c(14, 1, 1, 9, 2, 11, 20, 5, 22, 13, 21,                              2, 22, 10, 18, 0, 2, 1, 2, 15, 20, 23, 17, 3, 3, 16, 19, 23,                              3, 4, 4, 22, 2, 21, 20, 1, 19, 18, 17, 23, 23, 3, 11, 4, 23,                              4, 7, 2, 3, 19, 2, 18, 3, 17, 1, 9, 19, 23, 9, 6, 2, 1, 23, 21,                              22, 22, 22, 20, 1, 21, 6, 2, 22, 23, 19, 17, 19, 3, 22, 21, 4,                              10, 17, 23, 3, 7, 19, 16, 2, 23, 4, 5, 1, 20, 7, 21, 19, 2, 21)                    , count = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,                                 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,                                 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,                                 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,                                 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,                                 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,                                 1L, 1L, 1L, 1L, 1L, 1L, 1L))               , .Names = c("hour", "count")               , row.names = c(8L, 9L, 10L, 29L, 33L, 48L, 51L, 55L, 69L, 72L, 97L, 108L, 113L,                                118L, 126L, 140L, 150L, 171L, 177L, 184L, 202L, 230L, 236L, 240L,                                242L, 261L, 262L, 280L, 284L, 286L, 287L, 301L, 318L, 322L, 372L,                                380L, 385L, 432L, 448L, 462L, 463L, 495L, 539L, 557L, 563L, 566L,                                570L, 577L, 599L, 605L, 609L, 615L, 617L, 624L, 663L, 673L, 679L,                                682L, 707L, 730L, 733L, 746L, 754L, 757L, 762L, 781L, 793L, 815L,                                817L, 823L, 826L, 856L, 864L, 869L, 877L, 895L, 899L, 918L, 929L,                                937L, 962L, 963L, 978L, 980L, 981L, 995L, 1004L, 1005L, 1007L,                                1008L, 1012L, 1015L, 1020L, 1027L, 1055L, 1060L, 1078L, 1079L,                                1084L)               , class = "data.frame")          ggplot(x, aes(hour)) +        geom_bar(binwidth = 1, position = "dodge", fill = "grey") +       geom_density() +        aes(y = ..count..) +       scale_x_continuous(breaks = seq(0,24,4))      Happy for anyone with better stats vocabulary to edit this question, especially the title and tags.


title: Definition of multivariate regression coefficient

text: I know that the regression coefficient of $Y$ and $X$ is defined as $$\beta(Y,X) = \frac{\mathrm{Cov}(X,Y)}{\mathrm{Var}(X)}$$ Does this expression also hold in a multivariate regression with $Y$, $X$ and $Z$ as additional explanatory variable?


title: Unscented Kalman filter-negative covariance matrix

text: I have recently started working on the unscented Kalman filter. I coded the numerically stable version (i.e., square root Kalman filter) and used MATLAB for implementing. In the final update step, where we update the state covariance matrix using `cholupdate` I get an error.               Pk = cholupdate(expected_S,K*Sy,'-')      Basically, `expected_S'*expected_S-(K*Sy)*(K*Sy)'` is not positive semi- definite (where `expected_P=expected_S'*expected_S;`) and so `cholupdate` returns an error. I cannot understand why this is happening in the square root implementation. (I have checked the derivations and coding didn't seem to find an error.) PS: I am currently testing linear models on UKF just to check if its optimal for linear models, but I keep getting the above error (for local linear trend model). Please help and tell what can be done to overcome it.


title: Is there a way to continue a R/JAGS MCMC chain that did not converge?

text: I'm building a series of hierarchical models using R and JAGS, linked using the R2jags library. The runs are fairly long -- from several hours to several days. I've had the sad experience of running some chains that did not converge. In that case, is there a way to extend the chain, rather than starting over?


title: How to define a distribution such that draws from it correlate with a draw from another pre-specified distribution?

text: How do I define the distribution of a random variable $Y$ such that a draw from $Y$ has correlation $\rho$ with $x_1$, where $x_1$ is a single draw from a distribution with cumulative distribution function $F_{X}(x)$?


title: Analyzing (hierarchical?) clustered pair-matched binary data

text: We teach supplementary lessons in nearly two dozen local schools, and have two data sets of approximately four hundred records each from pre-post tests given at these schools. Each record contains pre and post values (correct, incorrect) for questions on 12 topics as well as whether or not there was an intervention (lesson taught) relating to that topic. Our goal is of course to assess the impact of the lessons taught. The pre and post test responses are pair-matched by student, and clustered by school since the lessons taught vary by school, in addition to the other factors generally accepted as valid for clustering at this level i.e; similar socioeconomic background, school culture, common instructors, etc. Without clustering, McNemar is the accepted test for analysis of this data, and several authors have explored various modifications to McNemar to allow for clustering including:   * Methods for the Analysis of Pair-Matched Binary Data from School-Based Intervention Studies. Vaughan & Begg. 1999. doi: 10.3102/10769986024004367   * Analysis of clustered matched-pair data. Durkalski et al. 2003. doi: 10.1002/sim.1438   * Methods for the Statistical Analysis of Binary Data in Split-Cluster Designs. Donner, Klar, Zou. 2004. doi: 10.1111/j.0006-341X.2004.00247.x   * Adjustment to the McNemar’s Test for the Analysis of Clustered Matched-Pair Data. McCarthy. 2007. http://biostats.bepress.com/cobra/ps/art29/ (Free to download) I have subsequently experimented with the Durkalski method as documented in McCarthy, since it seems to be deemed rather robust, as well as being the simplest for me to understand and code. However, none of the documented methods fit our case exactly as they use the matched pairs for pre-post or control-treatment only, and treat the clusters as a single class. We actually have matched pairs of pre-post in multiple control & treatment clusters, but this latter level of information is not used and discarding the ~50% of our data points from the control groups seems sub-optimal. Is anyone aware of a technique designed to analyze this data configuration, or someone whom might be interested in exploring this area? Thanks in advance!


title: Data simulation in R using covariance method

text: I simulated my data successfully in R by applying the R codes below. However, I need to simulate errors separately for each different Y outcome by using covariance method in R. I don't know how to use covariance to simulate errors in R. Could you help me about this problem? Thanks               library(MASS)          ## var1(Verb.comp), var2(word.Id), var3(word.att) ~ N(0, 1)     ## r_12 = 0.5, r_23 = 0.6, r_13 = 0.5          covmat <- matrix(c(1, 0.5, 0.5,     0.5, 1, 0.7,     0.5, 0.7, 1), ncol = 3, byrow = TRUE)          ## this is what my variance-covariance matrix looks like     covmat          ## simulate data     df <- mvrnorm(1000, mu = c(0, 0, 0), Sigma = covmat)     head(df)     ## estimate pm-corelations     cor(df)     evalA <- matrix(c(1,1,1,0,1,1,0,0,1),3,3)     head(evalA)          head(mydat <- as.matrix(df)%*%evalA)     colnames(mydat) <- c("VC", "WI", "WA")     head(mydat)          err<-rnorm(1000,0,1)     head(err)          Y.vc <- 0.8 + 1*mydat[,'VC'] + 0.54*mydat[,'WI']+0.51*mydat[,'WA'] + err     Y.wi <- 0.8 + 0.54*mydat[,'VC'] + 1*mydat[,'WI']+0.75*mydat[,'WA'] + err     Y.wa <- 0.8 + 0.51mydat[,'VC']+ 0.75*mydat[,'WI']+1*mydat[,'WA'] + err


title: Show that the probability that the first return to state $1$ occurs at time $k$ is $(0.5)^{k-1}$

text: > Suppose that the chain is intitially in state $1$, i.e $P(X_0 = 1) = 1$. Let > $\tau$ denote the time of first return to state $1$, i.e > > $$\tau = \min\{n > 0: X_N = 1\}.$$ > > Show that > > $$P(\tau = k) = (0.5)^{k-1}, k = 2, 3, ...$$ State $1$ only communicates with state $4$. I have already (correctly) got that $P_{14} = 1, P_{44} = \frac{1}{2}, P_{41} = \frac{1}{2}$. So to do this question, what basically happens is that my process will first go from $1$ to $4$. It will then stay in $4$ for some time $k$. Then, it will either remain in $4$ or go back to $1$. The probability of this happening is $$P_{14} \times P_{44}^k \times P_{41}^k = 1 \times (0.5)^k \times (0.5)^k = (0.5)^{2k}$$ which isn't the right answer. Where have I gone wrong? EDIT: Also, the next part tells me that using this relation and the definition of recurrence, I need to verify that state $1$ is recurrent. In the answers, they say > We need to show that $P(\tau = \infty) = 1$. Observe that > > $$P(\tau < \infty) = \sum_{k = 2}^{\infty} P(\tau = k) = \sum_{k = > 2}^{\infty} (0.5)^{k-1} = \sum_{j = 1}^{\infty} (0.5)^j = \frac{0.5}{1 - > 0.5} = 1$$ How have they managed to do this. I get what we want to show, due to the definition of recurrence, but why have they then worked it out for $\tau < \infty$ and how have they gone between each of the summation signs to get $\frac{0.5}{1-0.5}$?


title: How to calculate sample size needed for comparing the "change from baseline" scores between two groups?

text: I'd like to know how to calculate the sample size needed for comparing the "change from baseline" scores between two groups? Thank you in advance for your consideration


title: Which R package to use to conduct a latent class growth analysis (LCGA) / growth mixture model (GMM)?

text: I am trying to perform a latent class growth analysis (LCGA) and/or growth mixture models (GMMs) in R. The data I am using is an increasing number of forks of git repositories (discrete variable, not categorical), as you can see in this dataset. I tried `lavaan`, which helped me fit a latent growth curve model, but not to identify latent classes. I also tried `poLCA`, which only works for categorical polytomous variables, hence did not suffice either. What is the most appropriate R package for performing a latent class growth analysis on discrete variable data? The analysis I want to do is similar to that in Qureshi & Fang (2010): _Qureshi, I., and Fang, Y. 2010. “Socialization in Open Source Software Projects: A Growth Mixture Modeling Approach,” Organizational Research Methods (14:1), pp. 208–238._


title: Quantifying QQ plot

text: The qq-plot can be used to visualize how similar two distributions are (e.g. visualizing the similarity of a distribution to a normal distribution, but also to compare two artibrary data distributions). Are there any statistics that generate a more objective, numerical measure that represent their similarity (preferably in a normalized (0 <= x <= 1) form)? The Gini coefficient is for example used in economics when working with Lorenz curves; is there something for QQ-plots? Best, A.


title: Fitting ratios in multiple regression formula

text: I would like to ask a (probably very simple) question with regards to multiple linear regression. I have an experimental formula in the form: $$ Y \sim \frac{a_0 \cdot X_0}{(a_1 \cdot X1) * (a_2 \cdot X_2)} $$ and I want to set up a linear model which fits the $a_i$ coefficients so as to maximize the goodness of fit, how should I proceed? I thought to decouple the dependent variables by means of logarithms but I am not sure if this the right way to go. A hint seems to be given in this question (answer by whuber): _"[...] marginal changes in the explanatory variables are interpreted in terms of multiplicative (percentage) changes in the dependent variable"_. Analytically: $$ \begin{aligned} Y & \sim e^{\log{\frac{a_0 X_0}{(a_1 X1) * (a_2 X_2)}}} \\ & \sim e^{\log{a_0 X_0} - [\log{a_1 X_1} + \log{a_2 X_2}]} \\ \log{Y} & \sim \log{a_0 X_0} - \log{a_1 X_1} - \log{a_2 X_2} \end{aligned} $$ Not really sure how to work it out. Thanks for any help.


title: Adding a quadratic variable to regression

text: In my textbook I often see a quadratic term in the regression. For example let's say I have the regression model: $\log(\mathrm{wage})=\beta_0+\beta_1\mathrm{educ}+\beta_2\mathrm{exp}+\beta_3\mathrm{exp}^2+\varepsilon$. Why is it so common to put the squared term in there? Here's a quote from my book: > It often makes sense to add quadratic terms of any significant variables to > a model. Why is that? And if this is true, then if I'm studying endogeneity and a quadratic variable isn't in there is it likely to be that case the the squared term is endogenous if it's omitted?


title: Metric spaces and the support of a random variable

text: Does the use of metric spaces to describe the support of a random variable provide any greater illumination? I ask this after reading about how metrics spaces have been used to unify the mathematical measure theoretic nature of probability and the physical intuition that most associate with probability. You can read my inspiration here: http://www.arsmathematica.net/archives/2009/02/14/complete-metric-spaces-and- the-interpretation-of-probability/


title: Confidence intervals for repeatability

text: I have calculated the repeatability of individuals' responses to a stimulus using the methodology of Lessells & Boag (1987) Auk 104:116, where repeatability r = among-groups variance component / (among-groups variance component + within-groups variance component). How do I assign confidence intervals to my estimate of r?


title: Factor analysis on mixed (continuous/ordinal/nominal) data?

text: What approaches are there to perform FA on data that is clearly ordinal (or nominal for that matter) by nature? Should the data be transformed our are there readily available `R` packages that can handle this format? What if the data is of a mixed nature, containing both numerical, ordinal and nominal data? The data is from a survey where subjects have answered questions of many types: yes/no; continuous; scales. My aim is to use FA as a method for analyzing the underlying factors. I do not yet know what factors I'm looking for. However, condensing the underlying factors into a manageable number of factors is important. EDIT: Also, can I approximate a survey question answered on the Likert-type scale as a continuous variable? Thank you.


title: Does the ghyp R package create different GHD random data than the fBasic version?

text: I have been comparing two packages that are calculating generalised hyperbolic density (GHD) functions as I intend to become more familiar with these. But the results are different for each package. How do the implementations differ?               # Comparing GHD fBasics and ghyp functions     set.seed(4242)     plot(density(rhyp(100,2,0,1,0),bw = .5),col="black",ylim=c(0,.6), main="")     set.seed(4242)     lines(density(rghyp(100, ghyp.ad(lambda=1, alpha=2, delta=1, mu=0, beta=0)), bw=.5), col="red")     # Comparing random samples     set.seed(4242)     fb_ghd <- rhyp(100,2,0,1,0)     set.seed(4242)     gh_ghd <- rghyp(100, ghyp.ad(lambda=1, alpha=2, delta=1, mu=0, beta=0))     tail(cbind(fb_ghd, gh_ghd))


title: Is there a covariance MLE which takes into account independence relationships?

text: In the extreme case where all of the components of an $M$-variate observation are pairwise independent from each other, a multivariate normal distribution can be decomposed into the product of $M$ univariate normal distributions. For example, $$p \left( X_{1},X_{2},X_{3};[\mu_{1},\mu_{2},\mu_{3}]^{T},\left[\begin{array}{ccc} \sigma_{1} & 0 & 0\\ 0 & \sigma_{2} & 0\\ 0 & 0 & \sigma_{3}\end{array}\right]\right)=\prod_{i=1}^{3}\ p(X_{i};\mu_{i},\sigma_{i}).$$ The sample complexity of a univariate distribution is less than that of a full-covariance $M$-variate normal distribution, and thus fewer observations are required to get a good approximation of this model, since only univariate distributions must be estimated. With zero-mean observations $\{x^{(1)},\cdots,x^{(N)}\}$, the covariance MLE for the case where _the covariance matrix is isotropic (i.e. all observation components -- denoted with subscripts -- are pairwise independent)_ is simply $\hat{\Sigma}=\mathrm{diag}(\hat{\sigma}_{1,}\hat{\sigma}_{2},\cdots,\hat{\sigma}_{M})$ where $\hat{\sigma}_{i}=\frac{1}{N}\sum_{j=1}^{N}\left(x_{i}^{(j)}\right)^{2}$. This pairwise independence is represented as zeros in the _precision matrix_ $\Sigma^{-1}$, which is diagonal when $\Sigma$ is diagonal. As mentioned above, this is identical to (i.e. gives the same result as) separately estimating $M$ univariate distributions, and taking their product. My question is: In the _in between_ case, where there is some (not fully connected) dependency graph $\mathcal{G}$ between observation components, what does the covariance MLE look like then? As a specific example, if we say that $X_{1},X_{2}\perp X_{3}$, then the precision matrix must always have the form $$\Sigma^{-1} = \left[\begin{array}{ccc} \lambda_{11} & \lambda_{12} & 0\\ \lambda_{21} & \lambda_{22} & 0\\ 0 & 0 & \lambda_{33}\end{array}\right].$$ Therefore, what would the decomposition, and the covariance (or inverse covariance) MLE look like for a normal distribution such as the following? : $$ p\left(X_{1},X_{2},X_{3};[\mu_{1},\mu_{2},\mu_{3}]^{T},\left[\begin{array}{ccc} \lambda_{11} & \lambda_{12} & 0\\ \lambda_{21} & \lambda_{22} & 0\\ 0 & 0 & \lambda_{33}\end{array}\right]^{-1}\right)$$ For a given number of observations, the MLE should give the same result as the decomposed version of the pdf represented by $\mathcal{G}$.


title: What happens with the beta-binomial distribution, when n approaches infinity?

text: Short question: What happens to the beta-binomial distribution, when n increases to infinity? Is there a count distribution arising like it's for the classical binomial distribution?


title: What questions would a statistician ask about analysis of E. coli outbreak?

text: You may have heard about the recent enterohaemorrhagic _E. coli_ (EHEC) outbreak in Germany.   What questions would a statistician ask about EHEC analysis ?   I'm thinking of Q+As between reporters / public officials ↔ non-experts, say teachers and engineers with a Diplom / Master's degree but at most a smattering of statistics. (Is a picture, a map of EHEC land showing various strains of EHEC and the coverage of various tests, possible ?) Monday 20 June: I thought that the EHEC outbreak would be an area where statistics really matters in the world at large: what's the evidence for various causes, how can these be communicated to the public ? So, starting a bounty.


title: t.test returns an error "data are essentially constant"

text: R version 3.1.1 (2014-07-10) -- "Sock it to Me"     > bl <- c(140, 138, 150, 148, 135)     > fu <- c(138, 136, 148, 146, 133)     > t.test(fu, bl, alternative = "two.sided", paired = TRUE)     Error in t.test.default(fu, bl, alternative = "two.sided", paired = TRUE) :      data are essentially constant      Then I change just a single character in my fu dataset:               > fu <- c(138, 136, 148, 146, 132)      and it runs...               > t.test(fu, bl, alternative = "two.sided", paired = TRUE)              Paired t-test      What am I missing here?


title: Plotting changes in a three-valued ordinal variable across two time points using R

text: I have a dataset with 4025 participants across two time points. I have scored them on a three-point categorical variable (`Unlikely, Possible, Probable`) at each time point. I would like to visualize the various patterns of change (e.g. going from `Unlikely` at T1 to `Possible` at T2 or going from `Possible` at T1 to `Unlikely` at T2). I would also like some way of representing the number of participants in each of these clusters on the graph (somehow weighing by `N` and representing this by line thickness, number of lines etc.). The data is currently on the form:               id1, id2, variable_t1, variable_t2     1     500        0            1     2     501        1            0     ...      Any suggestions for how to do this? I have tried using ggplot2 and geom_line and grouping by id, but this graph looks very messy. I am looking for something more along the lines of a clustergram, but am open to suggestions. **Update** : I recently discovered Parallel Sets which is very close to what I would like to achieve. The only downside to this program is that it allows for very little customization of the resulting plot (e.g. rotating plot, adding titles and axis labels, manually adjusting size etc.). Although this is possible with a bit of post processing of the png file that the program can export to. (Now, is there a way of achieving the same with R and ggplot?). **Solution** : thanks to a reposting of my question by Tal (of R-bloggers fame) here there is now a solution for this question using R and lattice.


title: In plain language, why is there no VIF for binary outcome regression models?

text: As far as I know, the variance inflation factor is not computed with pseudo-$R^{2}$ or generalized $R^{2}$ in binary outcome models (e.g. logistic regression). Are there other measures of multi-colinearity than VIF appropriate to such models? Why ought or ought not we think about multi-colinearity in such models?


title: DLM out of sample errors

text: I'm using the DLM package to estimate a multivariate time series, I wanna check the out of sample forecasting, by estimating the residuals for 1, 6, 12 months ahead forecast? How can I calculate the 6 and 12 months ahead forcast like the kalman filter does for 1 month ahead forecast? Thanks Since this question is too general I update it. My question was: I have a times series which go from 1970 to 1990, and I want to check if my model gives a good out of sample fit. In order to do so I divide my dataset in two parts and starting from January 1980 I calculate 1 month ahead forecast errors, by dlm (f). Than I want to calculate 12 months ahead forecast errors, so once my t is january 1980 then february 1980, and so on. I would like to know if there's a way to do so? Thanks Maybe is better to specify my question a little more, because I did a mistake, sorry. I estimate the model recursevely from 1970:1 to 1980:1 (dlm), , at t=1980:1 I estimate y(t+12) and I compare it with the real y(t+12), then I estimate y(t+12) but t=1980:2, and so on. I would like to know which is the way to do it automically? Cause I thought that i can ran a dlm and use the dlmForecast and change every time the dataset through the window command, but I don't think it's the right way. Maybe for (i in 1:10){ fit = dlmFilter((window(data, start=1, end=12+i),mod), dlmForcast(FIT, nahed=12)