title: Why is p for 8 times heads out of 21 flips not 8/21?

text: Flipping a fair coin 21 times. The probability of getting 8 heads in any order is $$p = \frac{21!}{8!(21-8)!}0.5^8(1-0.5)^{21-8} = 0.097$$ I get that the Binomial formula takes into account the different sequences in which the heads can occur during the 21 flips. But in any case, I would at least had expected the Binomial formula to give a value a bit closer to $0.5$. The question is: What am I then calculating by $8/21 = 0.38$?


title: What is the best ensemble sampler for highly correlated parameter space?

text: I have a likelihood that I want to estimate the free parameters for it and I am using MCMC to estimate the parameters. Two of the free parameters are positions and I defined uniform priors and one has exponential and the last one lognormal prior. I have written an MCMC code with **pymc** and used `AdaptiveMetropolis` as ensemble sampler and I am running my MCMC code on simulated data where I know the final answer. I used `Raftery-Lewis` test in pymc to see how many chains I need to run in order to get convergence and based on that I chose the number of steps. But after running with the estimated steps, I haven't obtained the input results. My guess is that since my parameters are correlated, it hasn't converged. **My question** is that what is the best step method for correlated parameter space? How could I test which ensemble sampler works the best for my case?


title: Multiple regression with many zero values in DV, skewed data

text: I am not an expert, but I am trying to run a fixed effect multiple regression model in STATA. My data is panel data and has around 12000 observations, 348 districts over 36 months. My dependent variable is continuous, but there are several months were there is the value zero (nothing was installed/bought these months). The data is spread out pretty evenly across my months, but I still have around 9000 observations with the value zero. So, my data is skewed and I originally wanted to run a natural log with it, but my histogram looks awful, with a huge spike. Should I run the regression without logging it? If I log it or not, will my results be flawed? I have read about Zero-Inflated Poisson Regression as an option, but I am not really dealing with counts. Also, I feel like the zero values carry some significance and that I should not exclude these values. I would be grateful for any help, thank you.


title: Adjusting for too many variables

text: Suppose that you have a linear regression model with outcome variable $Y$ and explanatory variables $X_{1}, \dots, X_{10}$. Suppose $X_1$ is the main effect. Also suppose $X_{2}, \dots, X_{10}$ are associated with the outcome and the explanatory variable based on subject matter expertise. Should you include all the variables $X_{2}, \dots , X_{10}$ in the adjusted model?


title: Repeated measures logistic regression

text: I have a dataset that features a binary outcome, a binary predictor, and an unordered factor with 7 levels, and 120 subjects. Each of the 120 subjects were asked a binary question on seven issues, hence the 7-level unordered factor. Therefore, I have a total of 840 observations. Each of the 840 observations has a binary predictor, and a binary outcome. I want to regress the binary outcome, on the binary predictor of course. The outcome is slightly complicated. I will use mac vs. PCs as an example. People are asked, "Do you think macs are better or do you think PCs are better?" That is recorded as a 1, for preferring macs, and 0 for preferring PCs. Then they are asked "Is your opinion based on objective truth, or subjective opinion?" In other words, do you think it is objective (1) or subjective (0). Then the last question is "In 20 years, will people think macs are better or will people think PCs are better?" Again, 1 will be coded for macs. A new variable is then created. The variable is 1 if their current opinion matches what they think the future opinion will be. Previous pilot studies found that people tend to believe their stance on issues is what the future will be (...I know, a pretty intuitive finding...people tend to think they are right). The outcome is that variable, whether their current opinion agrees with what they think the future will be. So the point is to see if how people view things as subjective vs. objective, predicts whether they think their vision of the future is in line with their beliefs. The mac vs PC is only an example of an issue. In reality, there are seven issues, ranging from gun control, to abortion, etc., and each of those seven issue ranges in "objectivity" given by how the participants answered the objectivity question.   * `Participant ID`: 120 Factors    * `Objectivity`: Binary predictor of 0s and 1s    * `Outcome`: Binary outcome of 0s and 1s    * `Issue`: 7-factor level variable  So I tried something like this:                 # Dummy coding for all participants     contrasts(data$Participant) <- contr.treatment(120)  #$       # Deviance coding to compare to the overall mean     contrasts(data$Issue)       <- contr.sum(7)      model <- glmer(Outcome ~ Objectivity*Issue + (Objectivity*Issue|Participant),                     family=binomial, data=data)      This gives a bunch of warnings. I essentially want to see if my predictor predicts the outcome, based on the seven issue, while adding a random effect for the participant ID.


title: Alternative to visualize hyper-ellipsoid defined by variance-covariance matrix

text: Let say at the beginning I have two variables $a_{1}$ and $a_{2}$ of the same type then I would have the variance-covariance matrix defined as: $\sum = \begin{bmatrix} \sigma^{2}_{a_{1}} & \sigma_{a_{1},a_{2}} \\ \sigma_{a_{2},a_{1}} & \sigma^{2}_{a_{2}} \end{bmatrix}$, where diagonals are indication of variances and off-diagonals as covariances. I know that it is possible to visualize the confidence region by an error ellipse which its parameters would be defined as:   * the center of the ellipse is the mean (or expected value) of two variables,    * the orientation computed as $\alpha = \frac{1}{2}\tan^{-1}(\frac{2\sigma_{{a_{1}},a_{2}}}{\sigma^{2}_{a_{1}}-\sigma^{2}_{a_{2}}})$    * the axes length computed as squared root of eigenvalues of the covariance matrix.    * etc. Now if I have let say, $n$ variables, where in this way I will have variance- covariance matirx of size $n\times n$: $\sum = \begin{bmatrix} \sigma^{2}_{a_{1}} & \sigma_{a_{1},a_{2}} &\dots & \sigma_{a_{1},a_{n}}\\ \sigma_{a_{2},a_{1}} & \sigma^{2}_{a_{2}} & \dots & \sigma_{a_{2},a_{n}} \\ \vdots & \vdots & \ddots & \vdots\\ \sigma_{a_{n},a_{1}} & \sigma_{a_{n},a_{2}} & \dots & \sigma^{2}_{a_{n}} \end{bmatrix}$, the confidence region would be a hyper-ellipsoid (with dimension of $n$ ) which is not possible to visualize. My question: is there any alternative way to visualize or map this region in 2-dimension?


title: Difference-in-Differences Estimator for Logistic Regressions

text: I have a pre-post intervention study with four groups: 1) Pre-Intervention Control, 2) Pre-Intervention Treatment, 3) Post-Intervention Control, and 4) Post-Intervention Treatment. The outcome is a binary variable. There four other predictor variables. I am reanalyzing a previous study. In the previous study, they used a difference-in-differences estimator in a logistic regression, while controlling for the four predictors. With the indicators for Treatment and Time, the model is: $$\mbox{logit}(Pr(y=1 | \mbox{Time}, \mbox{Treat}))=\alpha_0\mbox{Time}+\alpha_1\mbox{Treat}+\alpha_2\mbox{Time}\cdot\mbox{Treat}+\beta x$$ However, reviewing the literature surrounding DiD estimators, it appears that using the DiD estimators in a logit regression (any regression with a nonlinear link function), results in the common trend assumption be violated. You can think of how there can't possibly be a common trend on the predicted outcome because it has support between 0 and 1. In addition, depending on where baseline was, differences on the index value (the stuff inside the logit function), can result in different margins on the probability of y. If baseline was somewhere in the middle, small difference on the index value would drastically change predicted probability, while if baseline started high, the differences would be minimal. So what is the practical solution here? How should I go about reanalyzing this data? What can I do, practically, to still be able to draw a causal conclusion based on the DiD estimator and the already established study design? Any ideas?


title: Classifier feature importance

text: If I train a GNB/LDA/kNN/other classifier I would like to know, in the model built, how important are features to classify or which feature(s) drives the classifier. For example in SVM models the importance of the feature is sometimes evaluated looking at weights magnitude, but for non-linear and generative models is more difficult to extract weights. Do you know a book/paper that could be useful to solve my question?


title: Possible outcomes of Approximate profile-likelihood estimator (APLE) for spatial autocorrelation

text: I've been working with spatial autocorrelation for a while and now I'm trying to move from more traditional estimators such as Moran's I or Geary's C to the new APLE estimator. I read Li's papers on APLE and also the reference on R-project. But I am still missing something very basic, as I still didn't understand something: the results of the APLE statistic is said to have a closed form but I didn't figure which one it is. For instance, while Moran's I can gom from -1 to 1, I couldn't figure which are the possible boundaries of APLE statistics. More interesting, I've been running tests in my older data, using the spdep package in R, and it seems to me that APLE can go up to 2. Any one has any tip on that? Please, any insights here would be very much appreciated.


title: p-value for hypothesis test with given correlation, sample size

text: I'm writing a statistical tool in Java. Here I'm calculation pearson correlation values (r). My next step is to validate if my correlation values are good enough to prove that my to variables are correlated. So I have two approaches in mind: a) since I have the sample size = n and the correlation = r already calculated, I could use the these values to calculate the p-value and validate if `p<alpha` (level of significance) But for this approach I'm missing the formula for calculating the p-value b) Since I'm planning to use a fixed alpha, I was thinking of precalculating the minimum r_min for given n. (e.g n=1 -> n=10000) and then check the calculated r against the precalculated r_min. Does anyone know about a good solution for this Problem? I already found this but here as far as I see it, the library needs to recalculates the r with would take the performance down.


title: Quantile regression with one predictor

text: Is there any closed formula for quantile regression with only one predictor? **Motivation**   I need to implement in SQL median regression with one predictor. It is quite easy to implement OLS with one predictor, because there is simple formulat for coefficients, but I have no idea how whether there is a formula for quantile regression.


title: Non-intuitive answer from a Poisson regression

text: I am not very familiar with Poisson regression, so I think I must have made a mistake in the below analysis: I am studying the effects of smoking on lung cancer rates. The dataset is provided here. The variable `smoking_status` is defined: > smoking status: coded 1 = doesn't smoke, 2 = smokes cigars or pipe only, 3 = > smokes cigarrettes and cigar or pipe, and 4 = smokes cigarrettes only, I modified the data a bit and made two new categorical variables: pipe/cigar smoker and cigarette smoker, to replace smoking_status. So a smoking status 1 maps to (0,0), 2 maps to (1,0), 3 maps to (1,1) etc. I also added a constant column to my dataset. This is all I did to the data. I then performed Poisson regression on this dataset, using a exponential link function. My hopes was that the coefficients of the two new variables would be positive, but instead only cigarette_smoker is positive. The confidence intervals do not contain positive points either. Have I analysed the data incorrectly, or is my data just wrong? **EDIT** The output (it's from a Python library Statsmodels )                            Generalized Linear Model Regression Results     Dep. Variable:                      y   No. Observations:                   36     Model:                            GLM   Df Residuals:                       31     Model Family:                 Poisson   Df Model:                            4     Link Function:                    log   Scale:                             1.0     Method:                          IRLS   Log-Likelihood:                -815.93     Date:                Thu, 31 Jan 2013   Deviance:                       1391.8     Time:                        13:19:32   Pearson chi2:                 1.22e+03     No. Iterations:                     7                          coef    std err          t      P>|t|      [95.0% Conf. Int.]     ------------------------------------------------------------------------------     x1             0.2596      0.006     44.097      0.000         0.248     0.271     x2            -0.1850      0.024     -7.775      0.000        -0.232    -0.138     x3             0.5327      0.031     17.101      0.000         0.472     0.594     x4             0.0004   7.95e-06     54.637      0.000         0.000     0.000     const          2.9593      0.046     63.903      0.000         2.869     3.050      The variables in order are age, smoke_cigar (0,1), smoke_cigarettes (0,1), population (in hundred of thousands), constant_term. Some example data: array([[ 2., 0., 0., 359., 1.], [ 4., 0., 1., 3270., 1.]]) with target deaths [ 22., 514.] respectively.


title: Predict values with multiple linear regression

text: I am learning `multiple linear regression`. I know for `simple linear regression` if I have the following `X` and `Y` values, then how to predict `Y'`. Like if I have the following datas:               X     Y   Y'     1     0   ?     0     1   ?     0     1   ?     0     0   ?      then I can calculate `a` and `b` and get `Y'` from the equation `y = a + bx` as shown here http://www.easycalculation.com/statistics/learn-regression.php. Now, I have multiple `independent variables` then what is the formula to calculate the `Y'` using multiple linear regression? My current datas are:               X1   X2   X3   X4   X5   Y  Y'      1    0    0    1    0    1  ?      0    1    0    1    0    0  ?      0    0    0    0    1    0  ?      1    0    1    0    0    1  ?      Can anybody show me the formula to predict `Y'`? I have searched a lot but those shows me scatter plots and analysis and all. I just want to get the `Y'` values. I didn't get any tutorial like I got for linear regression (link provided above) where it is shown clearly how to calculate `Y'` values. Edit: So, according to the answers my `X` and `Y` matrix are:                X = 1   1    0   0    1    0              1   0    1   0    1    0              1   0    0   0    0    1             1   1    0   1    0    0             Y = 1          0          0          1      and then I calculate `betas` with help of `(X'X)^-1 X'Y` equation. Then put the `x1`, `x2`, `x3`, values to predict `y1`, `y2`, `y3` and so on - right?


title: How to get p value and confidence intervals for nls functions?

text: I have 2 questions. 1) How can I have p.value for my 2 functions? My hypothesis is that I have a correlation between my function and my data. 2) How can I have a confidence intervals for my 2 functions?               library(ggplot2)     g <- function (x, a,b,c) a * (1-exp(-(x-c)/abs(b)))     X1 <- c(129.08,109.92,85.83,37.72)     Y1 <- c(0.7,0.5,0.39,-1.36)     dt1 <- data.frame(x1=X1,y1=Y1)     model1 <- nls(Y1 ~ g(X1, a, b, c),                start = list(a=0.5, b=60, c=50),control=nls.control(maxiter = 200))          ggplot(data = dt1,aes(x = x1, y = y1)) +           theme_bw() + geom_point() +           geom_smooth(data=dt1, method="nls", formula=y~g(x, a, b, c),            se=F, start=list(a=0.5, b=60, c=50))               f <- function (x, a, b, c) a*(x^2)+b*x+c        X2 <- c(589.62,457.92,370.16,295.98,243.99,199.07,159.91,142.63,     124.15, 101.98, 87.93, 83.16, 82.2, 74.48, 47.68, 37.51, 31,     27.9, 21.24,18.28)     Y2 <- c(0.22,0.37,0.49,0.65,0.81,0.83,1,0.81,0.65,0.44,0.55,0.63,     0.65,0.55,0.37,0.32,0.27,0.22,0.17,0.14)     dt2 <- data.frame(x2=X2,y2=Y2)     model2 <- nls(Y2 ~ f(X2, a, b, c),                 start = list(a=-1, b=3, c=0),control=nls.control(maxiter = 200))     ggplot(data = dt2,aes(x = x2, y = y2)) +            theme_bw() + geom_point() +            geom_smooth(data=dt2, method="nls", formula=y~f(x, a, b, c),            se=F, start=list(a=-1, b=3, c=0))      Thank you in advance


title: Conditions on Normal Random Variable

text: Let $X$ be a multivariate normal random variable with known mean and covariance, and $Y$ be another random variable. We know for sure $\|X−Y\|_2\le C$ where $C$ is a known constant. A few questions: Is it possible for $Y$ to be normally distributed? Can we say anything about its mean and/or covariance? If we know the mean of $Y$, can we say something about its covariance?


title: Divide-and-conquer approach for hierarchical clustering

text: I have a huge data set (33K), each represented as a bit-vector of 275-dimensions. basically my data set can be represented as a 33000 x 275 matrix. I want to cluster these bit-vectors. I have tried single link hierarchical clustering on a small data set, 3000 x 275, the result is promising. I know that single link hierarchical clustering algorithm is not scalable as the time complexity is $O(n^2)$. I am planning to apply divide-and-conquer approach, i.e., divide the dataset into chunks of equal size and cluster each chunks individually and finally merge the clustered chunks based on distance (if: $d(C1,C2)< t$; then: merge $C1$ and $C2$). The time complexity for my new approach is $O(p)O(1) + O(pq)$, where $p$ is number of chunks and $q$ is the average number of clusters in each chunk.   Note: I assume that when hierarchical clustering is applied, each chunck will take same amount of time and its constant for all chunks, thus $O(n^2)$ will become $O(1)$. I want to know, whether the above mentioned clustering approach is feasible and efficient. or is there any logic flaws in applying divide-and-conquer approach for clustering


title: Asymptotic normality in linear mixed-effects models

text: Estimates for the fixed effects in a mixed-effects model are known to be asymptotically normal in distribution (under mild assumptions). Are there any results that quantify the deviation from normality (in general, in simple cases, or just simulation studies)?


title: Relative predictive power of predictors used in time series models like kalman filter

text: How would we measure the predictive power of predictors in time series models. For e.g. in linear regression we have the magnitude and direction of the regression co-efficients and their p-values. Is there any measure like that to evaluate the performance of predictors in kalman filter?


title: How to calculate predicted hazard rates from a Cox PH model?

text: I have the following Cox PH model                (Time, Event) ~ X + Y + Z      I would like to get the predicted hazard **rates** (i am talking about hazard rates **NOT** hazard ratios) given specific values of X Y Z. I know the muhaz package can calculate the observed hazard rates, but I am interested in the predicted model. Is there a way to do this in R?


title: Variance Inequality for Random Vectors

text: I know that if X and Y are random scalar variables, then: \begin{align*} \mathrm{Var}(X+Y) & = \mathrm{Var}(X) + \mathrm{Var}(Y) + 2\mathrm{Corr}(X,Y)\sqrt{\mathrm{Var}(X)\mathrm{Var}(Y)} \\ & \le \left(\sqrt{\mathrm{Var}(X)} + \sqrt{\mathrm{Var}(Y)}\right)^2 \end{align*} I would like to know if the same inequality holds for random vectors, i.e.: $X , Y \in \mathcal{R^{n}}$ \begin{align*} \mathrm{Var}(X+Y) & = \mathrm{Var}(X) + \mathrm{Var}(Y) + \mathrm{Cov}(X,Y)+ \mathrm{Cov}(X,Y)' \\ & \le \left(\sqrt{\mathrm{Var}(X)} + \sqrt{\mathrm{Var}(Y)}\right)\left(\sqrt{\mathrm{Var}(X)} + \sqrt{\mathrm{Var}(Y)}\right)' \end{align*} Where **square root of a variance (positive definite) matrix is its Cholesky decomposition**. And the **inequality means that the difference is semidefinite positive.** As an extra hypothesis, I have that both $\mathrm{Var}(X)$ and $\mathrm{Var}(Y)$ are **diagonal matrices**. I have been looking for Cauchy Schwartz inequalities for variance matrices, and I found this: \begin{equation} \mathrm{Var}(Y) \ge \mathrm{Cov}(Y,X) \mathrm{Var}(X)^{-1} \mathrm{Cov}(X,Y) \end{equation} But I don't know how to give any use to it. Thanks a lot!


title: Good introduction into different kinds of entropy

text: I'm looking for a book or online resource that explains different kinds of entropy such as Sample Entropy and Shannon Entropy and their advantages and disadvantages. Can someone point me in the right direction?


title: Principal component analysis and variance normalization

text: I want to transform my data $X$ such that the variances will be one and the covariances will be zero. Furthermore the means should be zero. I know I will get there by doing Z-standardization and PCA-transformation, but in which order should I do them? I should add that the composed transformation should have the form $\mathbf{x} \mapsto W\mathbf{x} + \mathbf{b}$. Is there a method similar to PCA which does exactly both these transformations and gives me a formula of the form above?


title: How does gentle boosting differ from AdaBoost?

text: There is a variant of boosting called gentleboost. How does gentle boosting differ from the better-known AdaBoost?


title: characteristic roots of AR(p)

text: In Tsay's Financial Time Series, for a AR(p) process, > In the time series literature, **inverses** of the two solutions (to the > characteristic equation of AR(2)) are referred to as the **characteristic > roots of the AR(2) model**. From other areas of mathematics where the concept of characteristic equations also exists, e.g. difference equations, differential equation, eigenvalues of matrices, isn't the characteristic roots define as the solutions to the characteristic equations? If yes, what the book says is very unusual, and confusing. I hope this is not widely used in time series.


title: How to gently introduce epidemiologists/public health coworkers to advanced predictive modeling?

text: Coming from a social science and epidemiology background, my coworkers were trained on least squares regression, logistic regression, and survival analysis. They like to see 95% confidence intervals and p-values with the parameter coefficents, and are distrustful of more current predictive tools such as neural networks, CART, bagging & boosting, as well as penalized regression techniques.


title: Assuming training data as set of target funtions

text: In this 10th slide of http://www.cs.cmu.edu/afs/cs.cmu.edu/project/theo-20/www/mlbook/ch6.pdf presentation The Training data set $D$ is assumed as the set of target function . Actually $D = \left\{(x_i,C(x_i))\right\}$ but its assumed as $D=\left\{C(x_i)\right\}$ for i=1 to m Why its assumed like that and in 14th slide of same one ,its given like if Training data $D = \left\{(x_i,d_i)\right\}$ then $p(D/h) = \prod_{i=1}^{m}(p(\frac{d_i}{h})) $ ,but its possible $\iff$ $D= \cap_{i=1}^{m} d_1$ How it has to be interpreted ?


title: Calculating posterior of difference given posterior of two means

text: I am using R and `MCMCpack` to do a Bayesian analysis of some data that I have. I have generated posterior distributions (`postDist`) on the means of some parameters (`y1,y2,y3`) using `MCMCregress` (`postDist <\- MCMCRegress( x ~ y + z ,...)`). Now, I would like to take those posterior distributions on the means and generate a posterior distribution on the difference between the means. Is that a reasonable thing to do in a Bayesian analysis, and if so, how do you do it (either in theory or in R)?


title: bootstrapping for a parameter estimation

text: At the moment I'm studying bootstrap and I want to do some examples to get a feeling for it. Suppose I generate the following data:               set.seed(10)     data <- rnorm(1000,1,3)     number <- min(data[data <= 0.9])     number      so `number.min` will contain the minimum number of the data in the interval $(-\infty,0.9]$ which is equal to $-8.036491$. If we know there were again 1000 generated numbers using a normal distribution with standard deviation 3. Suppose the `number.min` is this time $a:=-12.2385$, is it still plausible to assume that the expected value is $1$? We denote `number.min` as a random variable $X$. My idea was to perform a hypothesis test. So $H_0: \mu = 1$. The problem here is we do not know the exact distribution of $X$ and I don't want to assume normal distribution. My idea was to use bootstrap. Under the assumption $H_0$, I can generate as many outcomes of $X$ as I want. So if I'm right, I'm interested in $$P(X\in [x_l,x_u]|\mu =1)=1-\alpha$$ when choosing $\alpha=0.01$ I want to find $x_u,x_l$. If $a\in[x_l,x_u]$ I would argue that $\mu=1$ is plausible, right? **1. Question:** This is theoretical nature: Is my idea correct, i.e. finding $x_l,x_u$ is this whole procedure from a statistical viewpoint ok? Assuming this is right, we can go over and implement this. First we generate enough outcomes of $X$                f <- function(n){     data <- rnorm(n,1,3)     number <- min(data[data <= 0.9])     return(number)     }     data <- replicate(500,f(1000))      This leads to 500 realisations of $X$. Now we can generate our boostrap samples               bst.smpl <- boot(data,?,1000)      **2. Question:** the boot function needs an argument "statistics", but I don't understand which it is in my case? Last I would just compute the confidence intervals               ci <- boot.ci(bst.smpl,conf=0.99,type="basic",index=1)      **3.Question:** Is there any easier way to decide if we are still using $\mu=1$?


title: Resources for learning about the Statistical Analysis of Financial Data

text: I realize that the statistical analysis of financial data is a huge topic, but that is exactly why it is necessary for me to ask my question as I try to break into the world of financial analysis. As at this point I know next to nothing about the subject, the results of my google searches are overwhelming. Many of the matches advocate learning specialized tools or the R programming language. While I will learn these when they are necessary, I'm first interested in books, articles or any other resources that explain modern methods of statistical analysis specifically for financial data. I assume there are a number of different wildly varied methods for analyzing data, so ideally I'm seeking an overview of the various methods that are practically applicable. I'd like something that utilizes real world examples that a beginner is capable of grasping but that aren't overly simplistic. What are some good resources for learning bout the statistical analysis of financial data?


title: AIC guidelines in model selection

text: I typically use BIC as my understanding is that it values parsimony more strongly than does AIC. However, I have decided to use a more comprehensive approach now and would like to use AIC as well. I know that Raftery (1995) presented nice guidelines for BIC differences: 0-2 is weak, 2-4 is positive evidence for one model being better, etc. I looked in textbooks and they seem strange on AIC (it looks like a larger difference is weak and a smaller difference in AIC means one model is better). This goes against what I know I have been taught. My understanding is that you want lower AIC. Does anyone know if Raftery's guidelines extend to AIC as well, or where I might cite some guidelines for "strength of evidence" for one model vs. another? And yes, cutoffs are not great (I kind of find them irritating) but they are helpful when comparing different kinds of evidence.


title: Determine if difference in class distribution is statistically significant

text: I have a dataset of some observations with class attribute of values 0 and 1. The dataset is quite unbalanced (class 1 – 15%, class 0 – 85%). Further this dataset consists of 5 years, and the distribution of class values varies between years (let's say class 1: 12%, 13%, 15%, 15%, 14%). I would like to know if this difference in class distribution is statistically significant between years. The results should be a matrix where each year is compared with each other year. I think I can use a vector of class values for one year and compare it with a vector of class values for another year with a t-test. Is this a correct approach? And is there a way how to get full matrix, so I don't need to run multiple t-tests (if t-test is the correct choice)? Any comment, help, link, or probably R code snippet would be very appreciated.


title: Weighing data during distribution fits

text: I currently have multiple sets of data in the form of a histogram. Using Matlab, the all look to be very well described by the negative binomial distribution. However, the far right tail not modeled as well as I'd like by the distribution. Plotted on a log-scale, the right tail in the actual data has a different slope than what the negative binomial distribution predicts, most of the time it predicts a smaller tail than reality. My question is then: is it possible to weigh the tail data so that the fitting pays more attention to it? If not, is there another method to solving this? I was thinking of trying to combine the negative binomial fit with some other curve that better represents the tail. At the same time I am looking into whether there is some reason for the tail from the data itself. I am rather new to statistics, so any guidance on how to implement your solution in Matlab would be much appreciated.


title: What to do when exploratory factor analysis results are different for complete-cases and imputed data?

text: I have a hundred items that I'm performing EFA on, with around 370 complete cases. Using parallel analysis to determine the number of factors to extract, EFA gave 9 factors, all of which make theoretical sense. However, I have 280 incomplete cases, which are missing answers for up to 10% of the items. Missing data is not at random, with a couple of items missing responses for 20-40% of all respondents. Using the `R` package `missForest`, I imputed the missing data and performed factor analysis on the resulting imputed dataset. This gave me 10 factors, around 5 of which are the same as the previous analysis and 5 'new' factors. What is the protocol when the results of imputed data analysis are different to those of a complete cases analysis? Any advice is much appreciated.


title: Singular Value Decomposition more columns than rows

text: I am confused by the singular value decomposition of a matrix. This may just be a misunderstanding of what singular value decomposition does so please be gentle with me. If I do a singular value decomposition on $X$ with $m$ rows and $n$ columns such that I get $X=U\Sigma V$. $V^*$ is supposed to be a nxn square matrix (see e.g. Wikipedia. However, for a matrix x e.g. size (40,100) I get in Julia (and also in R):               x = randn(40, 100)     xsvd = svdfact(x)     size(xsvd.Vt)     (40,100)      I am expecting (100,100). However, for               x = randn(100, 40)     xsvd = svdfact(x)     size(xsvd.Vt)     (40,40)      I get what I expect. Can someone please explain to me what is going on here? And possibly point me to somewhere where I can read up on the fundamentals?


title: How to create a model with dummy coding to put into log regression

text: How do I use information from a chi square table to build a statistical model. To be specific, I am using Likert-style survey information that measures student engagement. I would like to see if students with a specific intervention are more engaged than students who do not have that intervention. If they are, how does it affect their outcomes? I'm checking intervention>engagement>outcome. I do not know how to turn the results of the chi square test into a model that has new values based on salience. I will compute new variables by adding items/means together. The model will have be coded to include with the log regression measurement of the intervention which has all categorical variables.


title: Help me in solving question on sufficient statistics

text: Let X1 and X2 be a random sample of size 2 from the exponential distribution with mean θ. Let T=X1+X2 be a sufficient statistic for θ. Let θ^=X^2 (based on the 2nd observation) be an estimator of θ. Consider an estimator θ∗=E(θ^/T^=t). Show that Var(θ∗) For solving above question I have to find the conditional distribution of X2/X1+X2.but the difficulty is in solving their joint distribution.Please give me any hint for solving it.


title: Modeling empirical probability densities

text: Say we have a classic regression problem in wich we have a numeric outcome along with some predictors, both numeric and categoric. In a typical prediction problem we use to estimate some parameter like the mean of the outcome , but what if we want to estimate the "statistically different" densities based on that predictors? I know two possible ways of doing this, first would be a bayesian approach, and other could be some parametric model like AFT in survival analysis where we estimate the parameters of a family of distributions (logistic, lognormal etc.) restricted in some kind of problems where like that where the formulation allow it. The question is wether exists some general approximation to estimate that in a nonparametrical way (supose we have two-modal or some other no theoretical distributions)and if it could be a non-bayesian approach or we necessarely have to take a bayesian aproach to the problem. Thanks in advance,


title: Variable transformation in a ANOVA and multivariate analyses

text: I want to run an ANOVA, and the multiple correlation and regression analyses. Some of my variables to be included in these analyses are normally distributed and some are not. I log transformed the not normally distributed variables which helped to normalize them. My question is whether I can use both the log transformed and not transformed variables in the same analysis? For example, can I have 2 variables transformed and 3 not transformed and use them all in the multiple regression; can I have the dependent variable transformed and the independent variable not transformed, and use it in ANOVA? Thanks for help.


title: Understand the parameters for multiple regression - question based on the notations in a given algorithm

text: I've set up a regression model but am not sure if I'm doing it right. I'm using multiple regression to help do multi-class classification. So far I feel like I've understood the theory, but I'm confused in implementation/interpretation. My confusion is between two possible interpretations of the regression part of a Gentleboost algorithm, which I will show below. Let me explain my problem. My data, $X$, is 2D, with $n$ number of observations, and $m$ number of features per observation. Now, the response vector is a 2D array too. In other words, given an observation $x_i$, the response vector doesn't merely tell me what the $y$ value is, but it gives me an array of values, one value for each class. For example, if my data can be classified as being from one of two possible classes, then the response vector might be like this: $y = (0,1)$. So that's telling me that the value for the first class is 0, and for the second class it's 1. Now, on to my regression function, $g(x)$. Recall that for each observation I have $m$ features, in other words, I need $m$ regression parameters in the regression function. So that means I would store the parameters in a 1D array that holds $m$ values. To illustrate this, say my data $X$ is as follows, with 10 observations and 3 features each. And the observations can easily be separated into two different classes. So the number of classes for $X$ below is 2: > [1,2,3] > > [1,2,3] > > [4,5,6] > > [1,2,3] > > [4,5,6] > > [1,2,3] > > [1,2,3] > > [4,5,6] > > [4,5,6] > > [1,2,3] So this was $X$, a 10x3 array. This means that my regression function $g(x)$ should have 3 parameters, let it be a 3x1 array. Now because my response vector is 1x2 array, not not merely a single value, then if I put in an observation into $g(x)$, it should also give out a response _vector_ , with one value per class, i.e. a 1x2 array. But in order for $g(x)$ to return a 1x2 array for me, it needs to have more than 3 parameters. It has to have 3 parameters for _each class_ , that means in this case of two classes, it has to be 3x2 array. (I'm not sure if it really _has_ to, I'm just explaining according to my understanding and implementation of this using `statsmodels` in python. If I'm wrong, please correct me). So $g(x)$ was 3x2. Let's say it looks like this: > [a, d] > > [b, e] > > [c, f] And let's say that getting this 3x2 array of parameters above was called the "fitting" step. However, recall that I needed a $g(x)$ for _each class_. So to store all the regression parameters for my data, given that I have two classes, I decided that I need a 2x3x2 array, which is a 3D array. Like this: > [[a, d] > > [b, e] > > [c, f]] > > **,** > > [[p, s] > > [q, t] > > [r, u]] Is this true? Is this the correct understanding of **`Step 2(a)(ii)`** in this algorithm below? (In the algorithm below, the calculated working response is $z$, not $y$). ![enter image description here](http://i.stack.imgur.com/d5C6P.png) So my confusion is whether $g(x)$ should be:   1. A 2x3x2 array, as I mentioned above. (And if that's the case, then when I input test observations, $g(x)$ gives output of a 2x2 array, and how is that resulting array supposed to be interpretted?)   2. Or, instead, should I stick to the 3x2 $g(x)$ that I got in the fitting step, and consider each column to be $g_j(x)$? In other words, given that my 3x2 $g(x)$ was  > [a, d] > > [b, e] > > [c, f] , should I now consider that the first column is $g_0(x)$, and the second column is $g_1(x)$? This would mean that the parameters for the regression function for the first class are `a, b, c`, and for the second class they're `d, e, f`. Is this the correct way to implement/understand the above step **`Step 2(a)(ii)`** in the algorithm? If I do this, then I believe that giving test data to $g(x)$ would produce a 1x2 array, which makes more sense as the response vector.


title: power analysis when mean and standard deviation are given

text: I have been asked to perform a sample adequacy test (i.e., I suppose, to verify what is the minimum difference that can be quantified between two grous) where a new sample (of which mu, sigma and sample size are given) is compared with a reference sample. I was thinking about t-test power analysis but I am also given the standard deviation and the sample size of the other group. What is the most appropriate approch to handle this problem? Possibly in R....


title: What does p-value mean in R?

text: I want to apply a t-test on two independent samples, and decide whether sample_2 is greater than sample_1. Below are the results for two different `t.test`s I ran on the samples : default parameters : t.test(sample_1,sample_2)                       Welch Two Sample t-test          data:  sample_1 and sample_2     t = -1.8795, df = 1121.445, p-value = 0.06043     alternative hypothesis: true difference in means is not equal to 0     95 percent confidence interval:      -4.45823402  0.09580241     sample estimates:     mean of x mean of y       45.08116  47.26238       alternative "greater" :               > t.test(sample_1,sample_2,alternative="greater")                  Welch Two Sample t-test          data:  sample_1 and sample_2     t = -1.8795, df = 1121.445, p-value = 0.9698     alternative hypothesis: true difference in means is greater than 0     95 percent confidence interval:      -4.091665       Inf     sample estimates:     mean of x mean of y       45.08116  47.26238       alternative "less" :               t.test(sample_1,sample_2,alternative="less" , paired=FALSE)                  Welch Two Sample t-test          data:  sample_1 and sample_2     t = -1.8795, df = 1121.445, p-value = 0.03022     alternative hypothesis: true difference in means is less than 0     95 percent confidence interval:            -Inf -0.2707661     sample estimates:     mean of x mean of y       45.08116  47.26238       can someone please explain me what are the meaning of p-values in these two tests? the first p-value 0.06, is this the probability that they are equal, or that they are **not** equal? And the second value 0.97, does that mean sample_2 > sample_1 , or vice versa : sample_2 <= sample_1 ? This is a very basic question but I cannot understand from other resources. Any help would be appreciated, thanks.


title: Do coefficients of logistic regression have a meaning?

text: I have a binary classification problem from several features. Do the coefficients of a (regularized) logistic regression have an interpretable meaning? I thought they could indicate the size of influence, given the features are normalized beforehand. However, in my problem the coefficients seem to depend sensitively on the features I select. Even the sign of the coefficients changes with different feature sets chosen as input. Does it make sense to examine the value of the coefficients and what is the correct way to find the most meaningful coefficients and _state their meaning in words_? Are some fitted models and their sign of the coefficients wrong - even if when they sort-of fit the data? (The highest correlation that I have between features is only 0.25, but that certainly plays a role?)


title: Polynomial Contrasts for Regression

text: I cannot understand the usage of polynomial contrasts in regression fitting. In particular, I am referring to an encoding used by `R` in order to express an interval variable (ordinal variable with equally spaced levels), described at this page. In the example of that page, if I understood correctly, R fits a model for an interval variable, returning some coefficients which weights its linear, quadratic, or cubic trend. Hence, the fitted model should be: > $write$ ~ $52.7870 + 14.2587*X - 0.9680*X^2 - 0.1554*X^3$, where $X$ should take values $1$, $2$, $3$, or $4$ according to the different level of the interval variable. Is this correct? And, if so, what was the purpose of polynomial contrasts? Thanks!


title: simulate ANOVA data to follow a specific eta-squared

text: I am simulating RM-ANOVA data with three treatment conditions (one-way within subjects designs). I would like to be able to add an effect to the data so that eta-squared has a specific value. Imagine a more simple case for two independent groups and Cohen's d. Here, if I know both standard deviations in the two groups, I can have a look at the formula of Cohen' d to see that it's possible to multiply the mean of both standard deviations (the denominator) with the desired effect size, say 0.4, and I get the difference both means should have (numerator) and add this to one of the groups. see here I am wondering if that's also possible for eta squared (or some of its variants, partial etc.). It is somewhat more difficult because one has to work with sums of squares instead of mean differences. Do you have any suggestions on how to achieve this?


title: Showing characteristic function is infinitely differentiable

text: Suppose that X is a real-valued random variable. Suppose that there exists a constant M > 0 such that the support of X lies entirely in the interval $[-M,M]$. Let $\phi$ denote the characteristic function of X. Show that $\phi$ is infinitely differentiable. If infinitely differentiable is equivalent to absolutely continuous, then $$\int_{-M}^{M}|\phi(t)dt < \infty$$


title: Best practice to understand a dataset

text: My supervisor once told me that before I run any classifier or do anything with a dataset I must fully understand the data and make sure that the data is clean and correct. My questions: What are the best practices to understand a dataset (high dimensional with numerical and nominal attributes)? Practices to make sure the dataset is clean? Practices to make sure the dataset doesn't have wrong values or so?


title: What is the added value of using a 2-part model over an OLS model on a subsample when estimating health care expenditures?

text: I want to estimate the impact of having health insurance on health care expenditures and have a hunch that it is best to use a two-part model to first estimate the probability of using any health care and then estimating the amount spent by those who used health care. However, I am not clear on the advantage of using a two-part model over just running ordinary least squares regression on the subsample of people who used any health care. Any insight into the advantages of the 2-part model over OLS on the subsample of users, and the differences in these two approaches would be much appreciated!


title: Is this a correct statement of the probability chain rule?

text: Everywhere I see the chain rule written in this order: $$ \Pr (A,B,C) = \Pr (A|B,C) \Pr (B|C) \Pr (C). $$ Is this the same as $$ \Pr (A) \Pr (B|A) \Pr (C|B,A), $$ where the conditioning is done with different variables?


title: Likelihood ratio test to compare two predictions

text: I have two predictions from two different types of methods. "predictedHousePrices1" is a continuous variable and the output of a prediction from a RandomForest model, "predictedHousePrices2" is the output of the predict() function a RidgeRegression model. I'd like to compare which one better explains the variability in the real data. I'm wondering if a likelihood ratio test is the best way to do this, for example:               m1 <- lm(realHousePrices~predictedHousePrices1) # R^2 = 0.25     m2 <- lm(realHousePrices~predictedHousePrices2) # R^2 = 0.30      Is it correct to use a likelihood ratio test to check if an $R^2$ of .25 is "significantly" more than 0.30, for example:               m3 <- lm(realHousePrices~predictedHousePrices1+predictedHousePrices2)     library("epicalc")     lrtest(m3, m1)      Or is there a better way to do this?


title: Why is normalizing seasonal adjustment factors done additively?

text: One of the methods for dealing with seasonal influence is to establish a multiplicative factor for each season within a "year". For example, this happens with exponential smoothing models of type $(*,*,M)$. These adjustment factors need to be normalized - if there are $m$ seasons in a "year", then the seasonality normally represents only $m-1$ degrees of freedom. The standard normalization requirement for adjustment factors $f_1, \dotsc, f_m$ seems to be to require that $\sum_{i=1}^m f_i = m$ - that is, the arithmetic mean of the factors is $1$. This is certainly the case in Hyndman et al., "Forecasting wtih Exponential Smoothing" (section 8.2). I would have expected the requirement that the geometric mean be $1$, i.e., $\prod_{i=1}^m f_i = 1$. This because the standard for _additive_ seasonal adjustment is that $\sum_{i=1}^m f_i = 0$; these two are equivalent under taking logarithms. Is there a good mathematical reason for the choice, or is it just a convention?


title: Variance on the sum of predicted values from a mixed effect model on a timeseries

text: I have a mixed effect model (in fact a generalized additive mixed model) that gives me predictions for a timeseries. To counter the autocorrelation, I use a corCAR1 model, given the fact I have missing data. The data is supposed to give me a total load, so I need to sum over the whole prediction interval. But I should also get an estimate of the standard error on that total load. If all predictions would be independent, this could be easily solved by : $Var(\sum^{n}_{i=1}E[X_i]) = \sum^{n}_{i=1}Var(E[X_i])$ with $Var(E[X_i]) = SE(E[X_i])^2$ Problem is, the predicted values are coming from a model, and the original data has autocorrelation. The whole problem leads to following questions :   1. Am I correct in assuming that the SE on the calculated predictions can be interpreted as the root of the variance on the expected value of that prediction? I tend to interprete the predictions as "mean predictions", and hence sum a whole set of means.   2. How do I incorporate the autocorrelation in this problem, or can I safely assume that it wouldn't influence the results too much? This is an example in R. My real dataset has about 34.000 measurements, so scaleability is a problem. That's the reason why I model the autocorrelation within each month, otherwise the calculations aren't possible any more. It's not the most correct solution, but the most correct one isn't feasible.               set.seed(12)     require(mgcv)          Data <- data.frame(         dates = seq(as.Date("2011-1-1"),as.Date("2011-12-31"),by="day")     )          Data <- within(Data,{     X <- abs(rnorm(nrow(Data),3))     Y <- 2*X + X^2 + scale(Data$dates)^2     month <- as.POSIXlt(dates)$mon+1     mday <- as.POSIXlt(dates)$mday     })          model <- gamm(Y~s(X)+s(as.numeric(dates)),correlation=corCAR1(form=~mday|month),data=Data)          preds <- predict(model$gam,se=T)          Total <- sum(preds$fit)      Edit : Lesson to learn : first go through all the samples in all the help files before panicking. In the help files of predict.gam, I can find :               #########################################################     ## now get variance of sum of predictions using lpmatrix     #########################################################          Xp <- predict(b,newd,type="lpmatrix")           ## Xp %*% coef(b) yields vector of predictions          a <- rep(1,31)     Xs <- t(a) %*% Xp ## Xs %*% coef(b) gives sum of predictions     var.sum <- Xs %*% b$Vp %*% t(Xs)      Which seems to be close to what I want to do. This still doesn't tell me exactly how it is done. I could get as far as the fact that it's based on the linear predictor matrix. Any insights are still welcome.


title: test for trend in ordinal response variable

text: I have an ordinal response variable and multiple explanatory variables-both continuous and categorical. What I'd like to do is check for any univariate trend between the the different explanatory variables and the response variable. is it reasonable to just use a logistic ordinal model for each explanatory variable as long as the proportional odds assumption is fulfilled?


title: Should factor loadings be dominated by items' ranges of answer options?

text: If an exploratory factor analysis is done with some 1-5 agreement items and some 0/1 "choose all that apply" items, theoretically how much of a spurious tendency would there be for the 1-5 items to load on one or two factors and the 0/1 items to load on a separate set of one or two factors? (I've heard arguments both for and against the idea that correlations tend to be much higher among items of a similar scale. My own experimentation/simulation has not found much of an effect.)


title: Formula for autocorrelation in R vs. Excel

text: I am trying to figure out how R computes lag-k autocorrelation (apparently, it is the same formula used by Minitab and SAS), so that I can compare it to using Excel's CORREL function applied to the series and its k-lagged version. R and Excel (using CORREL) give slightly different autocorrelation values. I'd also be interested to find out whether one computation is more correct than the other.


title: Probability Higher level

text: Alice's school is planning to take some students from her class on a field trip. Alice is really excited about it. There are a total of $S$ students in her class. But due to budget constraints, the school is planning to take only $N$ students for the trip. These students will be picked randomly. And each student has equal chance of being picked. Alice's friend circle has M students including her. Though she is excited about the field trip, she will enjoy it only if there are at least $K$ of her friends with her on the trip. She is wondering what are the chances of that happening. She needs your help. Tell her the probability that she will enjoy given that she goes on the trip.


title: How to normalize data to let each feature lie between [-1,1]?

text: I have samples with each sample has n features, how to normalize these features to let feature values lie between interval [-1,1], please give a formula.


title: Tukey's HSD after a randomized block design in R

text: I want to run a Tukey's HSD on my data after performing a three factor anova with block as a fourth random factor. The anova works fine in R, using the following code,               model<-aov(dv~f1*f2*f3+Error(block/(f1*f2*f3)), data=*.txt)      where dv is my dependent variable, f1, 2 and 3 are the three fixed factors and *.txt is some text file which has the data to be analysed. When I attempt then, to do a Post hoc (Tukey's HSD) using the package agricolae, the following error shows up,               out<-HSD.test(model,"f1", group=TRUE) #### used this command          "Error in as.data.frame.default(x[[i]], optional = TRUE, stringsAsFactors = stringsAsFactors) :        cannot coerce class "c("aovlist", "listof")" to a data.frame"          out<-HSD.test(model,"f1*f2*f3", group=TRUE) ### used this as command as well          "Error in as.data.frame.default(x[[i]], optional = TRUE, stringsAsFactors = stringsAsFactors) :        cannot coerce class "c("aovlist", "listof")" to a data.frame"      Can someone help me with this please? Sorry about the earlier error.


title: Compressing an integer series

text: I have a set of integer series $S_1$, $S_2$, ... $S_n$. Each series has 3600 data points. Each data point is a positive integer. Each data point is stored as an an `unsigned int` requiring 4 bytes. So, storing the entire series requires 4 * 3600 bytes. This is too large a size and these series need to be compressed. The data is queried in the following way and (any compression mechanism should minimise the error on every query).   1. A few series are selected, say, $S_p, S_q$ and $ S_r$ (Not that any number of series might be selected, I am using 3 as just an example)   2. The selected s eries are summed to get $S_T$ such that, $S_T[i] = S_p[i] + S_q[i] + S_r[i]$, for $i = 1$ to $3600$    3. The minimum value and the maximum value of $S_T$ are computed. It is possible that the peaks in the series are aligned. So each of the series cannot be thought as as independent. Thoughts?


title: Minimum distance estimator

text: I'm trying to replicate the paper of Blundell et al. (2008) to separate permanent and transitory shocks on income on a panel dataset. He solves the non-linear system of equations using Chamberlain's minimum distance estimator (it is shown in the appendix of the paper), but I haven't found a library, nor in R, nor in STATA, that does the trick. Can someone help me?


title: How to examine the change of event sequences

text: Let's assume we have a sequence of events $x_1, x_2, ...,x_n$ and each event can be described as a categorical variable from domain $\{A, B, C...\}$. The time interval between two consecutive events is not fixed. An example could be a researcher's publication record. Each publication is an event and the domain is the set of all fields, {Statistics, Data Mining, ...}. Of course, the time interval between two consecutive publications is not fixed. My question is, are there any math/stat/data mining/.. tools to analyze these kind of data, so we could answer following questions?   * **Given a sequence, we can tell whether there is a significant change in the value of events**. For example, in the beginning, the researcher's publications are mostly about statistics, and for some reasons the interest has changed and more machine learning papers are published.   * **Given a group of sequences, we can tell whether there is a significant change in the value of events in many of them**. I guess the answer of this question is based on the answer of first question.   * **Given a sequence, we can tell whether there is a significant change in the time interval**. For example, a researcher publishes 1 article per month in the beginning, but publishes 1 article per year after he got tenure:)   * **And there may be more questions....** A simple method is to draw a timeline chart for each sequence. But it would not be acceptable if we want to draw sound conclusions from the data. Thank you!


title: code needed for p values in GLMM

text: I need to get p values for the fixed effects in the following GLMM's I ran. Does anyone know of code that I can run that will give me the p values I need? At the moment the output from the ANOVA only gives me one p value and I believe I need a separate p value for each of the fixed effects in the models. Thanks in advance. Code is as follows - For GLMM 1 I ran this code -               m1<-lmer(step~Depth*threshold+(1|ind))     m2<-lmer(step~(1|ind))     anova(m1,m2)      For GLMM 2 I ran this code -               m2<-lmer(PDBA~step*threshold+Depth*threshold+(1|ind))     m3<-update(m2,~.-step*threshold)     anova(m2,m3)      and this one:               m2<-lmer(PDBA~step*threshold+Depth*threshold+(1|ind))     m4<-update(m2,~.-Depth*threshold)     anova(m2,m4)      When I ran GLMM 1 code this is what I got:               m2: step ~ (1 | ind)     m1: step ~ Depth * threshold + (1 | ind)     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)        m2  3 373235 373259 -186615   373229                                 m1  8 373225 373290 -186605   373209 19.767      5   0.001382 **     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1      summary               > summary(m1)     Linear mixed model fit by REML ['lmerMod']     Formula: step ~ Depth * threshold + (1 | ind)           REML criterion at convergence: 373184           Random effects:      Groups   Name        Variance Std.Dev.      ind      (Intercept) 196519   443.3         Residual             469370   685.1        Number of obs: 23473, groups: ind, 11          Fixed effects:                      Estimate Std. Error t value       (Intercept)      160.95895  134.80279   1.194     Depth              0.06438    0.44777   0.144     threshold2        51.18065   17.62222   2.904     threshold3         1.47733   21.43879   0.069     Depth:threshold2  -1.23654    0.60029  -2.060     Depth:threshold3  -0.09587    0.65088  -0.147          Correlation of Fixed Effects:                 (Intr) Depth  thrsh2 thrsh3 Dpth:2     Depth       -0.094                                 threshold2  -0.090  0.712                          threshold3  -0.075  0.588  0.567                   Dpth:thrsh2  0.071 -0.737 -0.745 -0.435            Dpth:thrsh3  0.064 -0.674 -0.490 -0.857  0.502      OUTPUT FROM SUGGESTED CODE BY SETH (IN COMMENTS)               Models:     m6: step ~ Depth + threshold + (1 | ind)     m5: step ~ Depth + threshold + Depth:threshold + (1 | ind)        Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)       m6  6 373227 373275 -186607   373215                                m5  8 373225 373290 -186605   373209 5.2901      2      0.071 .     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


title: How to draw an interaction plot with confidence intervals?

text: My attempts:   1. I couldn't get confidence intervals in `interaction.plot()`   2. and on the other hand `plotmeans()` from package 'gplot' wouldn't display two graphs. Furthermore, I couldn't impose two `plotmeans()` graphs one on top of the other because by default the axis are different.   3. I had some success using `plotCI()` from package 'gplot' and superimposing two graphs but still the match of the axis wasn't perfect. Any advice on how to make an interaction plot with confidence intervals? Either by one function, or advice on how to superimpose `plotmeans()` or `plotCI()` graphs. code sample               br=structure(list(tangle = c(140L, 50L, 40L, 140L, 90L, 70L, 110L,      150L, 150L, 110L, 110L, 50L, 90L, 140L, 110L, 50L, 60L, 40L,      40L, 130L, 120L, 140L, 70L, 50L, 140L, 120L, 130L, 50L, 40L,      80L, 140L, 100L, 60L, 70L, 50L, 60L, 60L, 130L, 40L, 130L, 100L,      70L, 110L, 80L, 120L, 110L, 40L, 100L, 40L, 60L, 120L, 120L,      70L, 80L, 130L, 60L, 100L, 100L, 60L, 70L, 90L, 100L, 140L, 70L,      100L, 90L, 130L, 70L, 130L, 40L, 80L, 130L, 150L, 110L, 120L,      140L, 90L, 60L, 90L, 80L, 120L, 150L, 90L, 150L, 50L, 50L, 100L,      150L, 80L, 90L, 110L, 150L, 150L, 120L, 80L, 80L), gtangles = c(141L,      58L, 44L, 154L, 120L, 90L, 128L, 147L, 147L, 120L, 127L, 66L,      118L, 141L, 111L, 59L, 72L, 45L, 52L, 144L, 139L, 143L, 73L,       59L, 148L, 141L, 135L, 63L, 51L, 88L, 147L, 110L, 68L, 78L, 63L,      64L, 70L, 133L, 49L, 129L, 100L, 78L, 128L, 91L, 121L, 109L,      48L, 113L, 50L, 68L, 135L, 120L, 85L, 97L, 136L, 59L, 112L, 103L,      62L, 87L, 92L, 116L, 141L, 70L, 121L, 92L, 137L, 85L, 117L, 51L,      84L, 128L, 162L, 102L, 127L, 151L, 115L, 57L, 93L, 92L, 117L,      140L, 95L, 159L, 57L, 65L, 130L, 152L, 90L, 117L, 116L, 147L,      140L, 116L, 98L, 95L), up = 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)), .Names = c("tangle", "gtangles", "up"     ), class = "data.frame", row.names = c(NA, -96L))          plotmeans2 <- function(br, alph) {     dt=br;   tmp   <- split(br$gtangles, br$tangle);        means <- sapply(tmp, mean);  stdev <- sqrt(sapply(tmp, var));       n <- sapply(tmp,length);       ciw   <- qt(alph, n) * stdev / sqrt(n)     plotCI(x=means, uiw=ciw, col="black", barcol="blue", lwd=1,ylim=c(40,150),  xlim=c(1,12));      par(new=TRUE) dt= subset(br,up==1);        tmp   <- split(dt$gtangles, dt$tangle);       means <- sapply(tmp, mean);       stdev <- sqrt(sapply(tmp, var));       n <- sapply(tmp,length);      ciw  <- qt(0.95, n) * stdev / sqrt(n)     plotCI(x=means, uiw=ciw, type='l',col="black", barcol="red", lwd=1,ylim=c(40,150), xlim=c(1,12),pch='+');     abline(v=6);abline(h=90);abline(30,10); par(new=TRUE);     dt=subset(br,up==-1);        tmp <- split(dt$gtangles, dt$tangle);       means <- sapply(tmp, mean);       stdev <- sqrt(sapply(tmp, var));       n <- sapply(tmp,length);      ciw <- qt(0.95, n) * stdev / sqrt(n)     plotCI(x=means, uiw=ciw, type='l', col="black", barcol="blue",   lwd=1,ylim=c(40,150), xlim=c(1,12),pch='-');abline(v=6);abline(h=90);     abline(30,10);     }          plotmeans2(br,.95)


title: Is there a library that implements a fast maximum rank correlation estimator?

text: I'm looking for an implementation of a fast maximum rank correlation (MRC) estimator. This will be applied to large-ish sparse matrices (~100,000 by 10,000) in a text-mining application. I'm working in python and R, so it would be nice to find something in those languages. Failing that, I could probably convert code from some other language. Any suggestions? Quick note: the best algorithm I've seen described is in Abrevaya (1999), which runs in $nlog(n)$ time. Wang (2007) has an "IMO" algorithm that he claims runs in $n^2log(n)$ -- much worse. Note: I've cross-posted this question on stackoverflow because I'm not sure which community is the right place to ask.


title: Time series clustering

text: I have many time series in this format 1 column in which I have date (d/m/yr) format and many columns that represent different time series like here:               DATE         TS1     TS2     TS3 ...     24/03/2003   0.00    0.00    ...     17/04/2003  -0.05    1.46     11/05/2003   0.46   -3.86     04/06/2003  -2.21   -1.08     28/06/2003  -1.18   -2.16     22/07/2003   0.00    0.23      With R, how I can group the time series that show similar trends?


title: Probability of a certain run of dice rolls

text: OK, so here's a simple problem: Suppose I roll a D6 one hundred times. What is the probability of _any_ number appearing _exactly_ 50 times? Now, the probability of rolling 50 ones is presumably given by the number of 100-roll sequences that contain exactly 50 ones divided by the total number of possible sequences. But I could roll 50 twos instead. Which suggests that you need to add up all six ways to meet the criteria... except that these cases are not mutually-exclusive. It's possible, for example, to roll 50 threes and 50 sixes. Freakishly unlikely, but possible. So... um... how do I attack this?


title: How to test for variance between conditions

text: I have some Nanostring results that depict the fold change in expression after various conditions are applied. The people I am working with would like to know how to see what conditions give the largest change. For example which genes have the largest variance between all of the conditions and which genes have the largest variance between only a few of the conditions, i.e., between condition 1 and 2. What I have done was a simple test of variance across the 4 conditions, using the var.s function in Excel Here is a link to the Excel function var.s. Here is a small example:               Gene Name   Tconv WT/Treg WT    Tconv DEL/Treg WT   Tconv WT/Treg DEL   Tconv DEL/Treg DEL     ACVR2A      1                    1.07               -1.02                1.01     AHR         1                    1.48                1.26                1.72     AchE        1                    1.51               -1.46               -3.32     Adrb2       1                   -1.13                1.04               -1.19     Arg1        1                    1.1                -1.61               -1.14     Arnt1       1                    1.07               -1.03                1.08     Arnt2       1                    1.53                1.3                 1.3      All of the fold change data can be found on pastbin here. The question I have is, is this even the best test to do for this, and if I wanted to then look for the largest variation between condition 1 and 2 for example, as well as all of the other combinations can I still use the same test?


title: How to interpret the overall statistically significant effect of a MANOVA (In words)?

text: I have the following analysis: I ran a MANOVA with two DVs (both empathy scales) and an IV (gender). Now the $p$-value (.03) tells me that there is a significant difference between men and women. However, I do not know which is higher on which variable. When I then look at the ANOVA table (given with SPSS) output, it tells me that there is no significant difference (all $p$-values larger than .05). How do I interpret this? Had men and women score differently or not? Also, my supervisor asked me to look at the means, however they are only different by maybe one point, so I cannot see there which one should be larger (and the MANOVA only tells me the joint effect of the two DVs together). When I ran a t-test to have a look at the individual means, they were also non-significant. I am a bit lost what to do with this and would highly appreciate some help in regards how to put this interpretation into words.


title: why isn't the the marginal distribution needed when using a conjugate prior?

text: What is a good explanation as to why you wouldn't have to integrate to find the posterior when you use a conjugate prior. Most examples (like for instance: http://www.youtube.com/watch?v=0XD6C_MQXXE) multiply the prior times the likelihood and end with 'look how nice the (unscaled) posterior is of the same distribution as the prior and we didn't need to integrate'. But they never calculate the marginal distribution and divide by it. That kind of makes the conjugate look like a diversion from the fact that the integration problem still isn't solved at all..so why isn't the marginal distribution needed?


title: $a_n$ consistency and other consistency

text: From Jun Shao's Mathematical Statistics > Definition 2.10 (Consistency of point estimators). Let $X = (X_1 , ..., > X_n)$ be a sample from $P ∈ \mathcal P$ and $T_n(X)$ be a point estimator of > $θ$ for every $n$. > > Let $\{a_n\}$ be a sequence of positive constants diverging to $∞$. $T_n(X)$ > is called $a_n$-consistent for $θ$ if $a_n[T_n(X) − θ] = O_p(1)$ w.r.t. any > $P∈\mathcal P$. $O_p(\cdot)$ is defined here. When $a_n := \sqrt{n}$, the consistency is called root $n$ consistency. Do one of $a_n$ consistency and weak consistency imply the other? if not, how about $a_n := \sqrt{n}$? Do one of $a_n$ consistency and strong consistency imply the other? if not, how about $a_n := \sqrt{n}$? Do one of $a_n$ consistency and $L^2$ consistency imply the other? if not, how about $a_n := \sqrt{n}$? Thanks and regards!


title: Parametric segmented or piecewise regression with heteroscedastic errors

text: I am fitting longitudinal data with an increase in variance over time. The standard physiological model is a bi- or tri-linear model with variable breakpoints. The estimated parameters are used to describe rate of growth and time to maturity. Is there an efficient method to account for within-segment heteroscedastic errors in a parametric piecewise regression model?


title: Rescaling columns of bootstrapped matrix when using matrix multiplication to compute correlation coefficients

text: If I have $n$ samples and want to compute the pairwise correlation coefficients between $l$ and $m$ variables, I can construct the matrices $A_{l,n}$ and $B_{n,m}$. If I then subtract the mean and divide by the standard deviation from the rows of $A$ and columns of $B$, the matrix product $AB$ will contain the pairwise Pearson correlation coefficients of my $m$ variables with my $l$ variables. My question is this: If I want to bootstrap this procedure to get quantiles for each correlation coefficient, is it necessary to rescale the matrices $A$ and $B$ after I sample with replacement the rows of $A$ and columns of $B$, or can I scale once and resample? My intuition is that the rescaling for every bootstrap iteration is necessary.


title: Hierarchical regression with partially nested variables

text: I would like to do a logistic regression of vaccination coverage of, for example, $HPV$. Let's assume that we have 50 states ($S$). 30 states have an urban part and a country part (urbanity variable $U$). The other 20 have only a country area. 10 of the states have population areas where individuals speaks only Spanish and other areas only English ($L$). The linguistic region correlates only slightly with urbanity. For example: $S$: 1,.., 50; $U$: urban/country; $L$: English/Spanish.   * S(1)-U(urban)-L(english)   * S(1)-U(country)-L(english)   * S(1)-U(urban)-L(spanish)   * S(1)-U(country)-L(spanish)   * S(2)-U(urban)-L(english)   * S(2)-U(urban)-L(spanish)   * S(3)-U(urban)-L(english)   * S(3)-U(country)-L(english)   * S(4)-U(urban)-L(english)   * S(4)-U(country)-L(spanish) ... where the number indicate a state. The situation is that $U$ is nested in $S$ but not each $S$ has both levels of $U$. This is also the case for the linguistic variable. If $U$ would be nested completely in $S$ and $L$ completely in $S$, the model would look like (in R-terminology for lmer): `y ~ (1|S) +(1|S:U) + (1|S:L)` where $y$: 0/1 vaccinated Is it possible to estimate odds ratio using this model even if variables $U$ and $L$ are not completely nested? Or Should be used another model, for example: `y ~ (1|S) + U + L` where`(1|S)` is random variable and `U` and `L` are fixed variables. With fixed levels $U$ and $L$ the nesting character would be lost and the variables would be characterstics of the individuals,that is, the aggregation of individuals would not be considered. Problems which arise using fixed variables: would give significant results or narrow CI due to the high number of individuals. I would be very happy if someone could give me some hints how to cope with partially nested variables.


title: Residual plot for nonlinear regression

text: I have a couple of questions regarding performance of nonlinear regression models.   1. Are the residuals from a nonlinear regression model supposed to be randomly distributed too (as in linear regression)?   2. I am comparing two nonlinear regression models (non-nested). What model performance indicators can I use for the purpose? My understanding was to compare RSS, RSE, residual plots, autocorrelation and residual normality plots.   3. The models have different Y axes. So I tried standardizing variables and then re-estimated RSS, RSE. Is that the right approach? What other model performance measures should be estimated / plotted with standardized variables?


title: What makes an econometric model robust?

text: I was reading a paper on robustness (http://econ.ucsb.edu/~doug/245a/Papers/Robustness%20Checks.pdf) and they say: "To determine whether one has estimated effects of interest, $\beta$; or only predictive coefficients, $\hat{\beta}$ one can check or test robustness by dropping or adding covariates." What does a model being robust mean to you? Is this the only way to consider it in an econometric sense?


title: Kalman Filter with Linear Constraints

text: A question on this topic has been asked before: Combining a linear Kalman Filter with additional linear constraints? and I checked out some of the references given: http://academic.csuohio.edu/simond/pubs/IETKalman.pdf I will probably use estimate projection. However I have a couple of questions.   1. For the log-likelihood calculation, (i.e. when fitting the model) should I use the original or projected estimates? I presume the former?    2. When forecasting should I use the original or projected estimates, in this case I think the latter? Could someone please confirm which to use in each case and ideally why? Lastly for the prefect measurement method shown in the link above: Why is the number of restrictions forced to be less than the number of variables? I would have thought mathematically they could also be equal to the number of variables?


title: What are Quasi correlations?

text: Someone know what quasi correlations are in a DCC GARCH output. I am familiar with the Pearson correlation, always between -1(perfect negative) and +1(perfect positive). In my output I have a quasi correlation of 1,6 between 2 indices. The command that I runned: mgarch dcc (var1 var2 = , noconstant), arch(1) garch(1) constraints(1 2)


title: pseudo R-squared for model estimated with maximum likelihood

text: I want to compute a pseudo-$R^2$ for a model whose parameter estimation was based on maximum likelihood (function `likfit()`, package `geoR`, R software). I tried to compute the $R^2$ proposed by Maddala (1983) which compares the maximized likelihood for the model without any predictor and the maximized likelihood for the model with all predictors. I got a really low value (0.01%). Did I miss something? Are there other $R^2$s which are more appropriate than the $R^2$ of Maddala?


title: Wolfram Mathematica, MATLAB or something else?

text: I have only a little experience in math calculation tools. So please don't be strict if my question is not well defined :) Which tool (Wolfram Mathematica, MATLAB or something else) is better to use for statistics? Where can I find quickstart guides and tutorials for it?


title: MLEs of Poisson lambda values

text: I have sample data that I expect to contain values from at least several Poisson distributions (set around various lambda values). Some of these lambda values are nicely spaced, leading to what are visually obvious distinct distributions. My goal is to set confidence intervals (bins) around the most believable of these distributions so I can parse the data in a defensible way. So I'm thinking about using the conditional Chi-squared dispersion test once I have each lambda estimated... But I'm wondering what would be the best way to choose these values for lambda other than just simple visual inspection, which is where I'm stuck right now.


title: Two simple or one complex model, BIC and likelihood

text: I have a set of data points with a total number of Nt. I know a priori that the data comes from two distinct processes (distributions). I am trying to find the optimal model parameters together with the optimal partitioning for the dataset; such that N1+N2=Nt. I want to use the Bayesian Information Criterion, my question is which of the following BIC expressions is the correct one to optimize: As an example I will assume the problem of bi-linear fit of a dataset ![enter image description here](http://i.stack.imgur.com/xOywo.jpg) BIC= LogLikelihood - 0.5*k*log(N) [k: free parameters, N: number of datapoints] 1) Two separate linear fit models: LL(M1) - log(N1) + LL(M2) - log(N2) (k=2 for linear fit) 2) One single bi-linear fit model: LL(M)-2.5*log(Nt) (k=5; 4 (2 linear fits) + 1 (split point))


title: Kaplan-Meier, survival analysis and plots in R

text: What is the best package to to do some survival analysis and plots in R? I have tried some tutorials but I couldn't find a definite answer. TIA


title: Statistics and causal inference?

text: In his 1984 paper "Statistics and Causal Inference", Paul Holland raised one of the most fundamental questions in statistics: > What can a statistical model say about causation? This led to his motto: > NO CAUSATION WITHOUT MANIPULATION which emphasized the importance of restrictions around experiments that consider causation. Andrew Gelman makes a similar point: > "To find out what happens when you change something, it is necessary to > change it."...There are things you learn from perturbing a system that > you'll never find out from any amount of passive observation. His ideas are summarized in this article. What considerations should be made when making a causal inference from a statistical model?


title: Sample size for SEM. What are the alternatives

text: I have 3 latent variables with 3 measures each. I understand the conventional sample size would be 15 participants per measure (i.e. 9 x 15 = 135). I have 85 participants and don't think I will have enough time to collect many more (maybe another 15 max). Literature I read suggest no less than 200 for SEM. So my question is, if I can't use SEM what would be valid alternatives?


title: A limit theorem for non-independent variables

text: Let $X_n$ be a sequence of identically distributed (e.g., binomial $B(1,1/2)$) random variables which are _not_ independent (say, for any $n$ and $m$, $corr(X_n,X_m)=c$). What can be said about the limit of their mean $\sum_n X_n/N$? Motivation: a reason price variations assumed to be Gaussian is that they are composed by allegedly independent actions of many traders. However, in reality the traders are not independent, because, e.g., they all read the same news. So, $X_n$ is best viewed as a "set", not an "ordered sequence".


title: Accounting for the presence of a variable in regression?

text: So, as an example, let's say that I have data on a small auction that I conducted.   Let's say I am selling a car and there are 5 bidders, and I have data on each of their bids, who won, who participated, etc. Now, let's say that I have 5 different cars, so the auction occurs 5 separate times. For example, in one auction, all 5 buyers participate while in another only 3 buyers participate. I am wondering how to model out this situation for expected revenue (continuous variable -> linear model), while accounting for which buyers were present. I am modeling the data for all 5 of the auction together, but want to account for the presence or absence of a buyer in an auction. Is this simply a matter of including a dummy variable for each buyer, where 0 is absent and 1 is present? Or is there some alternative solution like tree models or other advanced tools that can be utilized to examine this problem. EDIT 1: I have more than 5 "auctions", but I just specified that number to keep things simple. EDIT 2: I'd be interested in hearing about non dummy-variable solutions to this problem.


title: How should I report errors for aggregate data?

text: I wish to "predict" the number of bugs my company will cause using a linear regression. For any given piece of development, the prediction is not great ($r^2 = .2$). However, my real task is to predict, company-wide, how many bugs we'll produce in a week. My errors are roughly zero-mean, so in this aggregated data set my $r^2$ is closer to .8. How should I report this? My coworkers (and admittedly myself) are not terribly familiar with statistics, but if I say something like "20% of the variance is 'explained' by the model" that makes it sound as if there is a long way to go. If 80% is explained, that seems quite good. Am I unfairly increasing my perceived correctness if I report the higher $r^2$ from the aggregate model?


title: Multiple regression with constraints on coefficients

text: I am performing a returns analysis. The idea is to regress a time-series of returns on the returns of various asset classes. The beta coefficients must be constrained such that sum of the coefficients is 1 and no coefficient is less than 0 or greater than 1. These beta coefficients can then be interpreted as explaining what % of returns are explained by exposure to the various asset classes. Are there any packages in R that let me setup the above regression and benefit from the attendant reporting on model fit statistics? Or do I need to do some homework on setting up constrained least squares optimization in R (please provide any references to recommended R packages)?


title: Conditional model using function tslm in R package forecast

text: I would like to use tslm with data that has intraday seasonality and a different pattern on business days and on non-business days. If data.ts is my time series then I would like to use something like               tslm(data.ts~season|businesss.dummy)      Thus I want to model season given that the dummy for this hour is True or False. I don't want               tslm(data.ts~season + businesss.dummy)      as this would just give a parallel shift on business days. I know that I can subset the data before applying the model and thus get business day data and non-business day data only but can I achieve this aim more elegantly using the right formula in tslm? Thanks!


title: Finding the PMF of conditional probability, poisson process. Don't understand where $10^6$ goes

text: "Customers arrive at a bank according to a Poisson process with rate 6 per hour. State (together with a proof) clearly the (conditional) probability mass function of the numbers of customers arrived during the first 20 minutes, given that 10 customers have arrived during the first hour." I did this: $$ P(X_{\frac{1}{3}} = x | X_1 = 10) $$ From the definition of joint probability, we know $$ P(X_{\frac{1}{3}} = x , X_1 = 10) = P(X_{\frac{1}{3}} = x | X_1 = 10) P(X_1 = 10)$$ Re-aranging gives: $$ P(X_{\frac{1}{3}} = x | X_1 = 10) = \frac{ P(X_{\frac{1}{3}} = x , X_1 = 10)} {P(X_1 = 10)}$$ Let's first calculate $P(X_{\frac{1}{3}} = x , X_1 = 10)$: This becomes: $$ P(X_{\frac{1}{3}} = x , X_{\frac{2}{3}}= y) = P(X_{\frac{1}{3}} = x) P(X_{\frac{2}{3}} = y) \hspace{1cm} y = 10 - x $$ Using the Poisson formula, you get: $$P(X_{\frac{1}{3}} = x) = \frac{e^{-2}2^x}{x!} \hspace 2cm P(X_{\frac{2}{3}}= y) = \frac{e^{-4}4^y}{y!}$$ And so $$ P(X_{\frac{1}{3}} = x , X_{\frac{2}{3}}= y) = e^{-6} \frac{2^x4^y}{x!y!} $$ Also, using the Poisson formula, we get: $$ P(X_1 = 10) = \frac{e^{-6}6^{10}}{10!} $$ So $$\frac{ P(X_{\frac{1}{3}} = x , X_1 = 10)} {P(X_1 = 10)} = \frac{e^{-6} \frac{2^x4^y}{x!y!}}{\frac{e^{-6}6^{10}}{10!}} $$ Which I get to be $$ \frac{10! 2^x 4^y}{x!y! 6^{10}}$$ But in the answers, it says it should be $$ \binom{10}{x} \left( \frac{1}{3} \right)^x \left( \frac{2}{3} \right)^{10 - x} $$ Why is this?


title: Max steps in lars

text: I don't know how many steps are necessary for `lars()` to select the variables till the algorithm proceeds to the saturated fit (especially using the `lars` package in `R`)? Can anyone give me a hint?


title: Relation between different ML models

text: Is there a paper or book that talks about the mathematical relation between different machine learning models - how they are different and how they can (sometimes) be equivalent? For instance logistic regression can be viewed as a single-layer perceptron. Other examples are here


title: Grouping trials decreases standard error?

text: My apologies is this is too rudimentary to be asked here. If it does not belong here, could someone recommend a more appropriate place to ask? A little context. I am in a Senior Physics Lab class in college. Our first "experiment" involved rolling dice and calculating the probability of rolling each face. (The purpose is purely for practice in applying and understand statistical analysis.) I chose to roll two dice separately with the intention of testing whether the two dice are identical. I know now I should have taken more samples, but I am not able to collect more data at this time. I originally recorded my data in trials of 20 rolls, and I took 15 trials for each die. I then calculated the counts of each face within each trial and divided by 20 to obtain an estimate of the probability in each trial. I then averaged these probabilities over all the trials and calculated the standard error (calculated as the standard deviation of the sample). However, this produced very large standard error, from about 40% to 60% for each face. Here is the data for Die 1.               Face  Avg. Probability  Error  Percent Error     1     0.13              0.08   59.6     2     0.20              0.08   39.0     3     0.18              0.09   47.9     4     0.13              0.08   64.8     5     0.18              0.10   57.1     6     0.18              0.10   55.4      However, I discovered that by grouping the trials into larger sets, the error decreased. That is, I grouped trials 1 through 5, 6 through 10, and 11 through 15, and I calculated the probability of obtaining each face as the count of the face within the "grouped trial" divided by 100. This produced the following values for Die 1:               Face  Avg. Probability  Error  Percent Error     1     0.13              0.03   26.6     2     0.20              0.01   5.0     3     0.18              0.04   20.0     4     0.13              0.06   43.5     5     0.18              0.06   32.9     6     0.18              0.02   11.1      Have I definitely done something wrong, or is the change in the averages and the decrease in error explainable? Would I be doing something to invalidate my results by using the larger "group trial" samples? I can provide my original data in CSV format if desired, and of course, if I am saying things that make no sense, I am open to correction. **Edit:** I've found the mistake in my calculations with the groups of 20 rolls, and the means are now identical as expected. However, this hasn't resolved the question of why the error is so much larger when I group my data into a greater number of samples, which was really the question I meant to ask. Yes, I understand that I could simply consider this problem using other means, but I'm still wondering why the error/standard deviation decreases. I think that question applies more generally than to this specific probability question. Or is the approach of trying to estimate the probability this way entirely invalid?


title: Using a random variable to represent a distribution?

text: The necessary context is, I want to model the world as states, and a state is a frequency distribution of people's opinions. So I wonder if I can use a random variable, say $T$, to model the state the world is in and hence the frequency distribution. Thanks.


title: Comparing regression coefficients of same model across different data sets

text: I'm evaluating two (2) refrigerants (gases) that were used in the same refrigeration system. I have saturated suction temperature ($S$), condensing temperature ($D$), and amperage ($Y$) data for the evaluation. There are two (2) sets of data; 1st refrigerant ($R_1$) & 2nd refrigerant ($R_2$). I'm using a non-linear, multivariate ($S$ & $D$), model for the regression analyses; polynomial of the 3rd degree. I would like determine how much less/more amperage on average, as a percentage, is being drawn by the second refrigerant. Or, some similar metric as a performance comparison. First thought was :   1. Determine the model to use... $Y = b_0 + b_1S + b_2D + b_3SD + b_4S^2 + b_5D^2 + b_6S^2D + b_7D^2S + b_8D^3 + b_9S^3$   2. Derive coefficients ($b_i$) from the baseline data ($R_1$).   3. Using those coefficients, for each $S$ & $D$ in the $R_2$ data set, calculate each expected amp draw ($\hat{Y}$) and then average.   4. Compare the $\hat{Y}$ average to the actual average amp draw ($Y_2$) of the $R_2$ data.   5. $\text{percent (%) change} = (Y_2 - \hat{Y}) / \hat{Y}$ However, since the 2nd refrigerant has slightly different thermal properties & small changes were made to the refrigeration system (TXV & superheat adjustments) I don't believe this 'baseline comparison method' is accurate. Next thought was to do two (2) separate regression analyses:   * $Y_1 = a_{0} + a_{1}S_1 + a_{2}D_1 + a_{3}S_1D_1 + a_{4}S_1^2 + a_{5}D_1^2 + a_{6}S_1^2D_1 + a_{7}D_1^2S_1 + a_{8}D_1^3 + a_{9}S_1^3$   * $Y_2 = b_{0} + b_{1}S_2 + b_{2}D_2 + b_{3}S_2D_2 + b_{4}S_2^2 + b_{5}D_2^2 + b_{6}S_2^2D_2 + b_{7}D_2^2S_2 + b_{8}D_2^3 + b_{9}S_2^3$ and then, for saturated suction temp ($S$), compare coefficients ($a_{1}$ vs $b_{1}$)... $\text{% change} = (b_{1} - a_{1}) / a_{1}$ However, again, these coefficients should be weighted differently... Therefore, the results would be skewed. I believe I could use a z-test to determine how differently weighted the coefficients are, but I'm not sure I fully understand the meaning of the output: $z = (a_{1} - b_{1}) / \sqrt{SE_{a_{1}}^2 + SE_{b_{1}}^2 )}$. But, that still wouldn't give me a performance metric, which is the overall objective... Thanks in advance for any guidance...


title: Path analysis of empirically evauluated model and sample size

text: I am new to the Path Analysis and ESM. I have been doing reading on this method, but still feel very much confused. I have a empirically evaluated failure model having 45 variables connected with each other in a cause and effect model. In my models, there are total 40 failure path. The longest paths has as many as 7 variables. After evaluating failure model, I have measured each variable through a questionnaire survey and have received total 50 replies. Now, I am wondering if this sample size is enough to perform Path Analysis to evaluate the "strongest failure paths"? or to perform SEM if I am attaching 5 latent veriables to my model? One of my friend suggested to use PLS techniques by using software 'SmartPLS" as it is less data sensitive.


title: What summary statistics to use with categorical or qualitative variables?

text: Just to clarify, when I mean summary statistics, I refer to the Mean, Median Quartile ranges, Variance, Standard Deviation. When summarising a univariate which is **categorical or qualitative** , considering both **Nominal** and **Ordinal** cases, does it make sense to find its mean, median, quartile ranges, variance, and standard deviation? If so is it different than if you were summarising a continuous variable, and how?


title: Entropy based on euclidian distances between datapoints / clusters centers?

text: Is it possible to define any useful entropy or conditional entropy which is based on the distance between datapoint(s) and cluster center(s), instead of basing on the number of points assigned to cluster like it is defined in for example to compute the v-measure ? I am not implying that I want an equivalent to v-measure, but I just wonder if it is possible and maybe useful to define a conditional entropy based on distances. So, What I'm looking for is a kind of conditional entropy based on distances, which allow to have an idea about the "homogeneity" and "completeness" of clusters, with respect to distances.


title: How to incorporate per-observation emission priors into HMM

text: I have a two-state HMM, in which my belief in emission probabilities depends on the observation. Basically, in addition to the two vectors of emission spectra (one for each state), I also have two vectors of priors that depend on parameters other than state "history". In other words, they have nothing to do with transition probabilities. I'm wondering what's the best way to incorporate them into the HMM. Using the priors just to define the initial states doesn't seem informative enough. On the other hand, using posteriors p(State)*p(Obs|State) rather than likelihoods p(Obs|State) as my emission spectra is probably also not a good idea as any sort of bias will be amplified in the HMM. Is there a better way to do this?


title: How to interpret correlations with negative numbers in SPSS

text: I'm trying to establish a bivariate Pearson correlation between two groups of variables in SPSS, however one of the groups has positive decimal numbers and the other negative decimal numbers. The results show a significant negative correlation between the two groups. If the negative numbers were positive instead this analysis would show a significant positive correlation. What Im wondering is, how do I interpret these results? Can I accurately say that there's a negative correlation? Thanks.


title: Posterior distribution and computation of probability of a future event

text: I am a beginner in statistics, and am self-studying from "Information Theory, Inference, and Learning Algorithms" by David MacKay. I've hit a wall with one of the questions, and was wondering if any of you could be so kind as to point me in the right direction. Admittedly, my title may be wrong or totally misleading as well - if you found it was, please let me know what would have been a better title for this question! The example which the questions are based on is as follows: > Bill tosses a bent coin N times, obtaining a sequence of heads and tails. We > assume that hte coin has a probability fH of coming up heads; we do not know > fH. If nH heads have occurred in N tosses, what is the probability > distribution of fH? What is the probability that the N+1th outcome will be a > head, given nH heads in N tosses. The question is (this is Exercise 2.8 on p. 30): > Assuming a uniform prior on fH, P(fH) = 1, solve the problem in the example > (given above). Sketch the posterior distribution of fH and compute the > probability that the N+1th outcome will be a head for > >   1. N = 3 and nH = 0; >   2. N = 3 and nH = 2; >   3. N = 10 and nH = 3; >   4. N = 300 and nH = 29; > He also states: > You will find the beta integral useful: (Since I'm new I guess I can only > give a link to the image - sorry about that): > http://i.stack.imgur.com/7MxQt.png Now I'm guessing that my difficulty has a lot to do with the fact that I'm new to statistics, and still learning the parlance so to speak. So let me break my question into a few pieces:   1. When he says "Sketch the posterior distribution of fH," what exactly does he mean by that? Does he mean it literally? (Yes...this is how ignorant I am...)   2. Since he doesn't give us fH, is he expecting us to come up with it on our own? Would it be something like (fH + 1)/(N + 2)?   3. First of all, how do I read the beta integral (what exactly is meant by Fa, Fb, pa) and secondly, how is it useful for this problem? Like I said, I'm a beginner, and trying to grok as much of what he said, but I think at this point I need some help. Thanks in advance!


title: Statistical significance of ordered binary vectors

text: I have a program to predict some values for people. For validation, I keep track of whether the prediction is correct or not, which gives me a binary vector with a length of about 600. To test if it is significant above a random predictor, I created 100 random binary vectors of predictions for the same values. This gives me 100 random vectors and 1 real vector. I want to test significance between the 1 real vector and the set of 100 random prediction vectors.   * Is there a good test to do this?    * Would a good methodology be to calc significance between two vectors at and then average p values?   * Would it be better to do something with distances, for example by calculating the Hamming distance between the real vector and the 100 random vectors, to get 100 distance values? If I do that, how can I know if that is significant or not? By checking skewness from normal?   * I just want to see if my predictions are greater than chance. Random predictions, averaged from 100 trials, are 0.33 accurate. Is there an easier way to test significance?


title: Time series forecasting using genetic algorithms

text: I am a beginner in the field of forecasting. I wish to know which are the best tools that can be used for forecasting future values in a time series using genetic algorithms. Are there any tools in MATLAB dedicated to forecasting using genetic algorithms similar to nntoolbox is dedicated to neural networks? Help is highly appreciated.


title: When will PCA be equivalent to ICA?

text: $X = AS$ where $A$ is my mixing matrix and each column of $S$ represents my sources. $X$ is the data I observe. If the columns of $S$ are independent and Gaussian, will the components of PCA be extremely similar to that of ICA? Is this the only requirement for the two methods to coincide? Can someone provide an example of this being true when the $cov(X)$ isn't diagonal?


title: Testing between two competing linear models with different lagged independent variables

text: What would be an appropriate statistical test to determine the best model for two competing linear models? Both models use the same independent variables (IVs); however, some independent variables are measured at different times. In the first model, the IVs lag the dependent variable (DV) by one period of time, and in the second model, the IVs lag the DV by a longer period of time. I'm interested in determining which lagged period is more appropriate to use. The DV for both models is measured at the same time period. In symbols: $$ y_{i,t} = x_{i,t-1} + y_{i,t-1} + z_{i,t-1} + \epsilon_i $$ versus $$ y_{i,t} = x_{i,t-2} + y_{i,t-2} + z_{i,t-2} + \epsilon_i $$ Are AIC, BIC, etc. appropriate for this?