title: Why do I get equal AIC, BIC and log likelihood for different models in LME framework?

text: I have two LME models with the same interaction, one containing both main effects and one containing only one main effect, say : $$ H\\_CE = Season + Crownlevel + Season:Crownlevel , random = 1|CollectorID $$ and $$ H\\_CE = Season + Season:Crownlevel , random = 1|CollectorID $$ There are 4 levels of each, and every combination of Season, Crownlevel and CollectorID The AIC, BIC and log likelihood of both models are completely equal. Given the formula for AIC being $$ \mathit{AIC} = 2k - 2\ln(L)\ $$ one would expect this to be different, even if the likelihoods are exactly the same. In the end, they have a different number of parameters. Or so I thought... Trying this toy example in R :               library(nlme)          Season <- rep(as.factor(rep(letters[1:4],each=4)),4)     Crownlevel <-rep(as.factor(rep(letters[11:14],4)),4)     CollectorID <-rep(letters[20:23],each=16)     X <-  model.matrix(~Season+Crownlevel+Season:Crownlevel)     B <- c(1,1,-2,2,0.3,0.4,0.4,2,3,1,-2,-3,-4,2,1,2)     H_CE <- X %*% B + rnorm(16*4)     KBM <- data.frame(Season,Crownlevel,H_CE,CollectorID)          model1 <- lme(H_CE~Season+Crownlevel+Season:Crownlevel,data=KBM,            method="ML",random=~1|CollectorID)     model1e <- lme(H_CE~Season+Season:Crownlevel,data=KBM,            method="ML",random=~1|CollectorID)      I get :               anova(model1,model1e)             Model df      AIC      BIC    logLik     model1      1 18 174.1834 213.0433 -69.09168     model1e     2 18 174.1834 213.0433 -69.09168      What am I missing here? Why are the numbers completely equal? It has to do something with the model specification, but I can't really see what. The model specification in itself is faulty, I know that. But I can't explain what makes it return a different set of parameters, but exactly the same residuals, likelihood and degrees of freedom :               > all.equal(residuals(model1),residuals(model1e))     [1] TRUE      * * * As fabians rightfully pointed out, both models are perfectly equivalent. Yet, I fail to see why in the AIC calculation the same value for the number of parameters `k` is used. The `k` in AIC uses the df, which explains everything.


title: Occam's razor (when is it appropriate to add another free parameter?)

text: So if I fit data to a function you can almost always decrease $\chi_{\nu}^2$ by adding more free parameters. However, this becomes ridiculous if you are fitting a 100-order polynomial to a straight line data set. There is obviously some point where you stop adding parameters! Is there any quantitative means to find this point other than intuition? I've heard about the F-test but not sure on how to apply it! Thanks


title: Aggregating two categorical variables

text: I have 2 categorical variables that I want to add into one single variable entitled `UnhealthyConsumption`: > FastFoodConsumption (1=never, 2=once a month, 3= 1-2 days/week, 4= 3-4/week, > 5= 4-5/week, 6= everyday), > > SnacksConsumption (1=never, 2=once a month, 3= less than once a week, 4= 1-2 > days/week, 5= 3-4/week, 6= 4-5/week, 7= everyday). I understand that if I would do a maximum frequency between the two might not be the best solution? I am interested in frequency of consumption. What would be the best way to aggregate these? Thank you


title: Why does the McNemar's test use chi-square and not the normal distribution?

text: I just noticed how the non exact McNemar's test uses the chi square asymptotic distribution. But since the exact test (for the two case table) relies on the binomial distribution, how come it is not common to suggest the normal approximation to the binomial distribution? Thanks.


title: Why goodness of fit via deviance and chisq is poor for poisson regression (glm)?

text: Simulate a poisson model:               set.seed(1)     predictor <- rnorm(100000, 2.5, 0.5)     # describe(predictor)     lam <- 0.98 * predictor     # describe(lam)          rp <- function(lambda){rpois(1, lambda)}     vrp <- Vectorize(rp)          response <- vrp(lam)     # describe(response)          fit <- glm(response ~ 1, offset=log(predictor), family="poisson")     summary(fit)      Regression:               Call:     glm(formula = response ~ 1, family = "poisson", offset = log(predictor))          Deviance Residuals:          Min       1Q   Median       3Q      Max       -2.8188  -0.8334  -0.1146   0.5777   4.3871            Coefficients:                  Estimate Std. Error z value Pr(>|z|)         (Intercept) -0.017707   0.002018  -8.773   <2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          (Dispersion parameter for poisson family taken to be 1)              Null deviance: 111501  on 99999  degrees of freedom     Residual deviance: 111501  on 99999  degrees of freedom     AIC: 361334          Number of Fisher Scoring iterations: 5          > pchisq(111501, 99999)     [1] 1      This above should be an exact fit, I don't get why `pchisq` is giving me so extreme number. Reference link: When someone says residual deviance/df should ~ 1 for a Poisson model, how approximate is approximate? ### Update I found this http://pj.freefaculty.org/guides/stat/Regression- GLM/GLM2-SigTests/ And it suggests to use pearson's residual instead of deviance to perform a goodness of fit... and it works much better.               ssr <- sum(residuals(fit, type="pearson")^2)     pchisq(ssr, 99999)


title: Undetermined power & retrospective study

text: My question is about a retrospective medical study in which no significant difference was found between two different treatments. What exactly is meant by "the lack of a significant difference does not imply equivalence because of the undetermined power of statistical tests"?


title: Laplace transform and density

text: It is true that the Laplace transform of a (positive) random variable characterises that random variable, just like its density? ($L_X(z) = E(exp(-Xz))$)


title: Performance indicator for imputation method using R

text: I have a 10 imputed air quality data using 7 different method of single imputation. can someone show how to compute MSE,RMSE,PA and MAE for my data, without manually type the equation in r.   * Prediction Accuracy (PA)   * The Coefficient of Determination   * Mean Absolute Error (MAE)   * Mean-Squared Error.(RMSE) Thank you.


title: Singular Value Decomposition and PCR

text: Can anyone guide me to understand the relation between Singular Value Decomposition (SVD) and Principle Component Regression (PCR)? I know that we can construct the principle components (PCs) using SVD, but I don't understand if the loadings which are defined in SVD are the same with the loadings in PCR. If $X = (x1,...,xn)$, SVD of $X = UDV^T$, the $U$ are the PCs of unit length and the columns of $V$ are the corresponding loadings of the PCs. In PCR, we have: $X=Pt_i+e_i$ where $T=(t_1,...,t_k)$ is the latent component matrix and $P$ is the matrix of loadings. My questions: Is the $V$ in SVD equal to the $P$ in PCR? I mean, can we use the same notation for both, $V=P$? Also, is $U=T$? (Note that, $n$: is the number of observations, $p$= number of predictors and $k$ = number of components.)


title: Expected value of latent utility in logistic regression

text: I'm looking for an analytical expression for the expected value of the latent utility in a logistic regression. **Setup:** There are two choices indexed by $i \in \\{0, 1\\}$ with associated utilities $u_i = w_i + \epsilon_i$, with the $\epsilon_i$ i.i.d. Gumbel(0, 1), aka Type I Extreme Value with location 0 and scale 1. The $w_i$ are constants. (Same setup as http://en.wikipedia.org/wiki/Logistic_regression#As_a_two- way_latent-variable_model) Let $Z \equiv \max\\{u_0, u_1\\}$, i.e. $Z$ is the utility you get when choosing the preferred alternative. **Question:** Is there a closed-form expression for the expected value of $Z$? If so, what is it, and how is it derived? **Additional Info:** If you know the answer, skip this section :-) The $\epsilon_i$ have CDF $F(x) = e^{-e^{-x}}$ and expected value $\gamma \approx 0.5772$, the Euler–Mascheroni constant -- here I'm copying from http://en.wikipedia.org/wiki/Gumbel_distribution, in case this is helpful to anyone reading this question. I did some googling, and I found something saying $\mathbb{E}[Z] = \ln(\Sigma_i(e^{w_i})) +$ constant. Is this correct? I wrote some R code to investigate, and it seems like $\mathbb{E}[Z] = \ln(\Sigma_i(e^{w_i})) + \gamma \approx \ln(\Sigma_i(e^{w_i})) + 0.5772$. Is this correct? Here's the code:               # Binary choice     # w_i + \epsilon_i with \epsilon_i ~ Gumbel     # Is Emax{w_i + \epsilon_i} = log(sum(exp(w_i))) + \gamma?          rgumbel <- function(n) {       # Slower than Runuran, but avoids installing a new package       u <- runif(n)       return(-log(-log(u)))     }     library(Runuran)  # Gumbel     rgumbel <- function(n) {       distribution <- udgumbel()       generator <- pinvd.new(distribution)       return(ur(generator, n))     }          gamma <- 0.57722 # Euler–Mascheroni constant     GetError <- function(w=c(5, 10), n=100000) {       epsilon <- matrix(nrow=n, ncol=2, data=rgumbel(2*n))       utility <- pmax(epsilon[, 1] + w[1],                       epsilon[, 2] + w[2])       emax <- mean(utility)       emax.analytical <- log(sum(exp(w))) + gamma  # Correct?       return(emax - emax.analytical)     }      The GetError function returns the difference between a sample mean and my candidate analytical expression. I played around and tended to see values very close to zero. Here's my very feeble attempt at getting an analytical expression for $\mathbb{E}[Z]$: $G(z) \equiv Pr[Z \leq z] = Pr[w_0 + \epsilon_0 \leq z] Pr[w_1 + \epsilon_1 \leq z | w_0 + \epsilon_0 \leq z]$ $G(z) = F(z - w_0) F(z - w_1) = e ^ {-e^{-z + w_0} -e^{-z + w_1}}$ $g(z) \equiv G^\prime(z) = e ^ {-e^{-z + w_0} -e^{-z + w_1}} (e^{-z + w_0} + e^{-z + w_1})$ $\mathbb{E}[Z] = \int z g(z) dz$ So, assuming there are no mistakes up to here, I need to calculate: $\int z e ^ {-e^{-z + w_0} -e^{-z + w_1}} e^{-z + w_0} dz$ plus a similar second term, but I'm not sure how to deal with this integral. Many thanks for any help, Adrian


title: Calculating probability in questionnaire answers

text: I'm designing a piece of research and want to quote the probability of two unique responses matching. 10 questions, for each answerer is required to select a,b or c. Therefore, if selecting randomly there's a 33.33% chance of ansering the same. What would be the overall probability of answering all 10 questions the same? is it 33.33 X 33.33 x 33.33 X 33.33 x33.33 X 33.33 x33.33 X 33.33 x33.33 X 33.33? How would one express this as a percentage chance?


title: Is this kernel proper?

text: ![enter image description here](http://i.stack.imgur.com/Ecpme.png) I think it is proper as it follows the rule that K(x1,x2) = f(x1)f(x2) Also they are both a function of only one of the points in the kernel respectively. Finally the product of two valid kernels is hence a proper kernel. I was wondering whether my reasoning was correct?


title: Seasonal exponential smoothing without trend

text: Why multiplicative property exists only for the exponential smoothing with seasonality and trend (Winter's additive and Winter's multiplicative models) and not for the exponential smoothing with only seasonality (without trend)? As I have seen in forecasting software, the exponential smoothing model with only seasonality only exists in an additive form. Why not exist in multiplicative form also? Is there a technical problem with this form (i.e. it does not exist either in the literature) or is it a deficiency of the software (i.e. it has not be implemented although it could be implemented)?


title: Suitable plot for 5 dimensional feature vectors?

text: I have a list of personality scores obtained from 100 people, based on the Big-Five personality test. Each person has one score for each of the five assessed traits. I put these scores into a 5 dimensional feature vector, so each person's personality is now represented by this 5 element vector. I would like to find out how personality traits correlate in each person. For example, if someone scores high on extroversion, how does he score on openness? And I would like to find out if I can generalise this and if people's scores tend to cluster or are more dispersed. My long term goal is to build a personality recogniser, and I would like to train a model for each personality cluster (if they exist). How can I approach this? In 3d space I would make a 3d scatter plot of the data and see how it is distributed. But how do I go on about doing it for 5 dimensions? If you know of any statistical measures that could help, I am open for suggestions (I'm not a statistician and don't know much about this). Thanks in advance!


title: What algorithm could be used to predict consumables usage given data from past purchases?

text: Thinking about a supposedly simple but interesting problem, I'd like to write some code to forecast consumable I'll need in the near future given the full history of my previous purchases. I'm sure this kind of problem has some more generic and well studied definition (someone suggested this is related to some concepts in ERP systems and the like). The data I have is the full history of previous purchases. Let's say I'm looking at paper supplies, my data looks like (date, sheets):               2007-05-10   500     2007-11-11  1000     2007-12-18  1000     2008-03-25   500     2008-05-28  2000     2008-10-31  1500     2009-03-20  1500     2009-06-30  1000     2009-09-29   500     2009-12-16  1500     2010-05-31   500     2010-06-30   500     2010-09-30  1500     2011-05-31  1000      it's not 'sampled' at regular intervals, so I think it doesn't qualify as a Time Series data. I have no data on the actual stock levels every time. I'd like to use this simple and limited data to predict how much paper I'll need in (for example) 3,6,12 months. So far I came to know that what I'm looking for is called Extrapolation and not much more :) What algorithm could be used in such a situation? And what algorithm, if different from the previous one, could also take advantage of some more data points giving the current supply levels (e.g., if I know that on date X I had Y sheets of paper left)? Please feel free to edit the question, title and tags if you know a better terminology for this. EDIT: for what it's worth, I will be trying to code this in python. I know there are lots of libraries that implement more or less any algorithm out there. In this question I'd like to explore the concepts and the techniques that could be used, with the actual implementation to be left as an exercise to the reader.


title: What do error bounds in forecasting represent?

text: What do error bounds actually mean in forecasting timeseries? For example, when I get a forecast I get the 85% and 95% high and low error bounds. I can also set my own error bounds to be calculated but am unsure of what an error bound in terms of timeseries forecasting means.


title: Literature on this kind of bias

text: Suppose I'm trying to investigate the lifespan distribution of light bulbs. The catch is that I can only observe each bulb at most $T$ time units. So if the bulb doesn't blow before $T$ I will not know what the value of the lifespan would be for that specific bulb, only that its lifespan was larger than $T$. So, the questions:   1. Is there a known name for this kind of bias?   2. Are there standard references to techniques to deal with this problem?


title: Why do we use k-means instead of other algorithms?

text: I researched about k-means and these are what I got: k-means is one of the simplest algorithm which uses unsupervised learning method to solve known clustering issues. It works really well with large datasets. However, there are also drawbacks of K-Means which are:   * Strong sensitivity to outliers and noise   * Doesn't work well with non-circular cluster shape -- number of cluster and initial seed value need to be specified beforehand   * Low capability to pass the local optimum. Is there anything great about k-means, because it seems that the drawbacks are beyond the good things about k-means. Please teach me.


title: Do we multiply the derivative of the sigmoid function when computing the errors of the weights connecting the last hidden layer to the output layer?

text: See I was under the impression that when we are calculating the error with respect to the weights connecting the last hidden layer with the output layer we are suppose to get the error of the output layer multiply it with the derivative of sigmoid, the current weights and the outputs of the last hidden layer. But if you go through the Andrew NG's course in Coursera you will see that he doesn't multiply the sigmoid derivative for the last layer!! This is even mentioned in the pdf he provided page 9. Also I find that when I do multiply it with a sigmoid derivative I get multiple peaks in the regularization vs training error plot, the NN converges to very high values of training errors in few iterations on the MNSIT handwritten digits dataset. Without the sigmoid derivative the regularization vs error is a smooth singular peak plot and goes on decreasing till 1000 iterations! Am I missing something here ? Cause from the calculus I learnt it seems obvious that sigmoid derivative needs to be multiplied !


title: Is there a Bayesian equivalent to a Wald-Wolfowitz runs test?

text: I have a sequence of observations and I would like to determine if the observations in the sequence are mutually independent. Wald-Wolfowitz is a non-parametric test that can be used to check for randomness in the sequence. I think I have a general idea of how to formulate the problem, but don't want to sprain my brain if it's been done before. [My literature searches and Google searches didn't turn up anything, but I may not be asking the question properly.] Thanks in advance.


title: What is the correct SD to use to get a 95% CI for skewed data?

text: Let $X = [94, 10, 100, 100, 16, 14, 100, 100, 70, 88, 100, 100, 12, 100, 100, 58, 32, 100, 32, 36, 98, 0, 100, 100, 100]$ where $X$ are students' scores (between 0 and 100), and note many full marks! The Question is what statistics will best describe the data (note data is non- Gaussian) **Option 1** If I fit a Gaussian using maximum likelihood I will get sample mean = 70.4, and SD = 37.96, so a mean +/- 1 SD gives an interval from 32.43 to 108.36. Finally, If I fit a Gaussian to the data $X$ using `normfit` in `matlab(R)` and obtain a 95% confidence bound on the mean and standard deviation I will get $$ \begin{aligned} \mu &= 70.4 ; &CI_{95\%} = [54.73, 86.06] \\\ \sigma &= 37.96 ; &CI_{95\%} = [29.64, 52.80] \end{aligned} $$ **Option 2** On the other hand, what if I use left / right SD instead? I.e., to report two SD values, SD_left and SD_right where: $$ \begin{aligned} SD_{left} &= \sqrt{\frac{1}{N_{left}} * \sum(X*I(X<\mu) - \mu)^2} &= 49.94 \\\ SD_{right} &= \sqrt{\frac{1}{N_{right}} * \sum(X*I(X\ge\mu) - \mu)^2} &= 29.45 \end{aligned} $$ where $\mu=70.4$ is the mean, $N_{left} = \sum(I(X<\mu)) - 1 = 9$ is number of samples less than the mean (minus 1 to remove bias) and $I$ is the indicator function which gives 1 if its argument is true or else 0; $N_{right} = \sum(I(X\ge\mu)) - 1 = 14$ In this case the interval around the mean is [20.46, 99.85], instead of the previous result, [32.43, 108.36]. Which one shall I go for, 1 or 2?


title: Use of two-sample Kolmogorov-Smirnov test to evaluate similarities between two different distributions

text: $X_1, X_2, \dots X_n$ and $Y_1, Y_2, \dots Y_n; n = 1000$ are two samples of physical quantities coming from the application of two different mathematical models to some independent and identically distributed (iid) data. The mathematical model used to generate $Y_i%$ is a simplified version (tuned by the parameter $r$) of the mathematical model used to generate $X_i$. My goal is to find the value of $r$ that makes the empirical distribution $F_Y(y)$ to be as similar as possible to the distribution $F_X(x)$. I decided to use the Two-Sample Kolmogorov-Smirnov test (in R) for different value of $r$ ranging in a specific interval. Is this choice correct? I know that the null hypothesis for the K-S test is that the two distributions are the same. However, I know for sure that the two distributions are different because the two mathematical models are different. Is it correct to evaluate the best value of $r$ by looking at the p-value and the D statistic coming from the K-S test?


Is there any statistical test that is parametric and non-parametric? This question was asked by an interview panel. Is it valid question?


title: How do you perform the likelihood ratio test in SAS?

text: I want to compare model $A$ with 5 variables and model $B$ with 3 variables. I want to use the likelihood ratio test to do this. Note that this is logistic regression. How do I do this in SAS?


title: Statistical significance in Excel

text: I have data from a market research survey which covers behaviour of different customer groups. To be specific, I am looking at the number of channels that a customer will visit before they buy a product. I would like to know whether the value of the customer (as measured by the value of the product they buy) is a meaningful way of splitting the data. I have calculated the average number of channels visited for the whole sample (let's call it X) and then for customers who bought only the most expensive products (Y) and separately for customers who bought only the cheapest products (Z). I would now like to check whether Y and Z are statistically significant. i.e. what is the probability that in my whole population I would find a customer who has bought a cheap product having visited Y or more channels beforehand. There are a number of statistical functions in Excel that look relevant (`CONFIDENCE.NORM`, `CONFIDENCE.T`, `T.TEST`, etc...). Which one should I use?


title: theoretical concerns in logistic regression

text: I have a dataset with 260 patients. I aim to study factors associated the certain finding in magnetic resonance imaging. I use logistic regression with six predictors. Regression yields to several significant predictors. I have philosophical concerns however. Nagelkerke´s R is only 20%. My model is however significant compared to empty model. Only 7% of the positive outcomes can be predicted. However diagnostics show that just ten cases has normalized residuals outside 1.96 SD. Moreover leverage criteria of 3x mean value is met in 250 cases. However I have been told that I do not have look too deeply to this fact with low R2 (http://www.theanalysisfactor.com/small-r-squared/ I would like to think so. My aim is not to construct a PREDICTIVE model but instead to study RELATIONSHIP between predictors and outcome. Therefore low R2 can be tolerated as far as my model is better than empty model. Moreover this is a highly clinical issue and it is impossible to include all relevant predictors. These most likely contribute much of the variance. My predictors are those which can be measured in reasonable manner. This fact highlights more the impossibility to construct a PREDICTIVE model. Is it valid to draw any conclusions from my analysis with statistical findings stated above? And is there any basis to differentiate between predictive logistic regression versus logistic regression studying relationship?


title: Maximum likelihood estimator for variance in two linear models

text: I am learning MLE's at my inference class and this is a problem I came accross. Consider two simple linear models. $y_{1j}=\alpha _1+\beta_{1}x_{1j}+\epsilon_{1j}$ and   $y_{2j}=\alpha _2+\beta_{2}x_{2j}+\epsilon_{2j}$ , $ j=1,2,...,n>2$ where $ \epsilon_{ij}$~$N(0,\sigma^2)$ ,   $\epsilon$ are independent and identically distributed. **Obtain the maximum likelihood estimators of** $\sigma^2$. For this can I use only **one** equation of the two .   Can I find with only using $y_{1j}=\alpha _1+\beta_{1}x_{1j}+\epsilon_{1j}$ and would it be the same answer I would get had I used the other equation


title: Estimate average percentage error based on another gaussian measurement

text: I have a model where the error is proportional to the throughput. This is, the observations I got come from a measurement instrument that has some error and it measures material going through in little buckets. I am approximating the model as follows $$ x_i = x~(1 + \epsilon_X) \\\ \epsilon_X \sim N(\mu_X,\sigma_X) $$ I am also able to determine via another instrument, the total material before and after, so I am able to state (after simplifying) $$ \delta_i = \delta + \epsilon_D' \approx x_i \\\ \epsilon_D' = \epsilon_D^1 - \epsilon_D^2 \\\ \epsilon_D^{\\{1,2\\}} \sim N(\mu_D, \sigma_D) $$ The mean $\mu_D$ does not matter, as it cancels out in $\epsilon_D'$, but **the mean $\epsilon_X$, if not zero, is what I am after**. Things get complicated because the variance of $x_i$ depends on the actual throughput, so it seems I have heteroscedasticity in the model. Put other way, $\epsilon_X$ is given as a random percentage of the measured quantity. So, I have a number of $\delta_i$ readings, with random noise; and another set of $x_i$ readings with noise as above; they should match, but they don't due to noise. Actually, I am interested in the problem of finding out if the noise in $x_i$ is not zero-mean. Any pointers would be appreciated. **Update** For what is worth, in addition to _Alecos'_ answer, I tried also what _seanv507_ suggested: by averaging $(1/n)\sum_{i=1}^n(Ln(x_i)-Ln(\delta_i))$, and then reverting with $Exp(.)$. It is very close to the ratio model (as said $1+x \approx e^x$ for small _x_ ) with the additional _advantage_ that $x_i$ is is never negative under this model (if this is what one wants).


title: Calculating NBA shooting consistency

text: What would be the proper way to evaluate/determine an NBA player's 3-point shooting consistency? For instance, I have a player that shoots 37% from 3-point range and takes 200 attempts all year. I was considering taking the rolling average 3-point% of an arbitrary number of shots (say 20). Then using those averages to determine the standard deviation from the 37% mean. Using a rolling sample size of 20 shots only allows for a precision of 5% in shooting percentage, but I'm concerned that using too many shots won't reveal the inconsistencies in performance. Is there a better approach to determine consistency?


title: Probability of certain value of $\chi^2$ given a number of trials

text: I am calculating a result and an associated $\chi^2$ via an iterative process that makes use of 3 random numbers (randomly sampling a Gaussian) in an equation, and I select as my final result the one that has the lowest value of $\chi^2$. I am trying to calculate the probability of ending up with this value of $\chi^2$, given the number of iterations that I use (with the expectation being that by using more iterations I'll end up with a lower value for $\chi^2_{min}$ given simply that there will be more attempts. I am trying to do this like: $Probability = \frac{e^{-\chi^2_{min}}}{\sum_i^{N_{trials}}e^{-\chi^2_{i}}}$ So essentially the probability is the exponential of the lowest calculated $-\chi^2_{min}$, divided by the sum of all calculated exp($-\chi^2$) values. I expected that for more trials I'd get a lower value for $\chi^2_{min}$ and the probability would approach 1 (with the assumption that the value in the numerator would dominate the sum in the denominator, so that they are approximately equal). What I see is actually for more trials I get a lower value for the probability as the denominator grows due to the sum. I would be very grateful if someone could tell me where my mistake is (I'm assuming that because the denominator can (and does) grow larger than 1, more trials will simply lead to a lower value for the probability. It isn't clear to me how to fix this though and any insight would be very welcome! Thanks.


title: How to dampen forecast to improve accuracy?

text: According to Armstrong there is ample empirical evidence that dampening trends in uncertain and complex long term forecasting helps improve accuracy/reduce forecasting errors. What I'm not able to find is a formula or an implementation of dampening trends. I know there exists dampening trends in exponential smoothing which is not what I'm looking for. I was looking for a simple formula which would dampen recursively the forecasted values based on a **user provided dampening factor**. Lets suppose I have the following dataset,first 8 values are actual and the remainder are forecast, how would you dampen the forecasted values with the user provided dampening factors that dampens the trend recursively i.e., as the time progresses there should be more dampening as the uncertainits increase in long run? I have artificially provided _dampening values_ and dampened the forecasted values, is there a mechanical process formula to create the _dampening value_ ? Any suggestions/recommendations of articles/formula spreadsheet implementation of dampening trends are greatly appreciated. Best Regards ![enter image description here](http://i.stack.imgur.com/xpvz4.png) Edit: I would like to clarify what I'm trying to achive in terms of an illustrative example. Lets suppose if we have a trend forecast as shown below,                   stl.ap <-stl(AirPassengers,s.window="periodic")               ## extract Trend component              trend.ap <- stl.ap$time.series[,2]         plot(trend.ap)           ## Ignore Seasonal data since this is an example to determine      ## how to dampen a trend          linear.ap <- tslm(trend.ap~+ trend + I(trend^2))     forecast.ap <- forecast(linear.ap, h=36)          plot(forecast.ap)      ![enter image description here](http://i.stack.imgur.com/IRnTE.jpg) As can be seen above, this is an unrealistic projection of trend, is there a formula that recursively dampens this trend. I would like to dampen the trend heavily as the time horizon increases.


title: Avoid negative results in Holt Winters forecasting

text: I understood that Holt Winters forecasting may results in negative values due to trending. I did reduce trending component value, but still forecast values are negative territory. Our data set will never be in negative values (like electricity data set, which never falls below ZERO). What sort of post algorithm techniques we can apply to make this value as non- negative value? Any help would be appreciated. This implementation is in Java, so I can't use `ets()` package at this point of time (and I am assuming `ets()` also can't avoid negative values). ![enter image description here](http://i.stack.imgur.com/7DJ5z.png)


title: correlation of bootstrap sample means

text: Given a sample $\\{x_1,\dots,x_n\\}$, $z_1$ and $z_2$ are two bootstrap realizations of sample means, that is, $$z_1 = \frac{1}{n}\sum\\{x\in\text{bootstrap sample 1}\\}$$ $$z_2 = \frac{1}{n}\sum\\{x\in\text{bootstrap sample 2}\\}$$ , how to compute the $corr(z_1, z_2)$? **UPDATE** To make sure I understand the problem correctly, here is what I'm trying to solve: ![enter image description here](http://i.stack.imgur.com/1AN3N.png) Here is what I tried, as @probabilityislogic suggested, $$z = \frac{1}{n}\sum(x_ik_i) ;k_1+k_2+\cdots+k_n = n$$ , since this is a bootstrap sample, so $x_i$ are all constant and $k_i$ are random variables that $$(k_1,k_2,\cdots,k_n) \sim Mult(k_1k_2\cdots k_n|n, p_1,p_2,\cdots,p_n)$$ . I tried to compute the mean and variance like this, \begin{align*} E(z) &= \sum_k\begin{pmatrix}n\\\ k\end{pmatrix}p^kz\\\ &= \sum_k\begin{pmatrix}n\\\ k\end{pmatrix}p^k(\frac{1}{n}\sum_{i=1}^n x_ik_i)\\\ &= \frac{1}{n}\left[\sum_k\begin{pmatrix}n\\\ k\end{pmatrix}p^k(x_1k_1+\cdots+x_nk_n)\right]\\\ &= \frac{1}{n}[np_1x_1+\cdots+np_nx_n]\\\ &= p_1x_1+\cdots+p_nx_n \end{align*} , and since each $k_i$ could be any number in $[0,n]$, then $p_i = \frac{1}{n+1}$, right? If so, $$E(z) = \frac{n}{n+1}\bar x$$ . To compute covariance, \begin{align*} Cov(z_1, z_2) = \sum_{k_{z1}, k_{z2}}\left[\text{Pr}(k_{z1},k_{z2})(z_1 - E(z))(z_2 - E(z))\right] \end{align*} , where $\text{Pr}(k_{z1}, k_{z2})$ is the joint probability of $k$ for $z_1,z_2$. Am I going the right way here for covariance of $z_1$ and $z_2$? I can't seem to figure this out here, seems too complicated to me.


I have seen the term _deep autoencoders_ in a couple of articles such as Krizhevsky, Alex, and Geoffrey E. Hinton. "Using very deep autoencoders for content-based image retrieval." ESANN. 2011. What's the difference between autoencoders and deep autoencoders?


title: Test whether simple slopes are different from zero in 3-way interaction in multiple regression

text: I'm trying to test whether 4 different slopes from a 3-way interaction in multiple regression are significantly different from zero. The four lines are plotted at 2 levels of each of the 2 moderators (lo-lo, hi-lo, lo-hi, hi-hi). Here's how we could test the significance slopes in a 2-way interaction model. I would like to extend this method to test the significance of slopes in a 3-way interaction model in multiple regression. A 2-way interaction in multiple regression takes the following form:                y = a + b1(X) + b2(Z) + b3(X)(Z)      With a 2-way interaction, one can examine whether the slopes at various levels of the moderator, Z, are significantly different from zero using the following equations:                b1 at Z = b1 + b3(Z)      where b1 is the slope of the predicted effects of X on Y at any particular value of Z                SE(b1 at Z) = (var(b1) + (Z^2)(var(b3) + (2Z)(cov(b1,b3))^(1/2)      where var(b1) is the variance of the b1 regression coefficient, var(b3) is the variance of the b3 regression coefficient, cov(b1,b3) is the covariance between the b1, b3 regression coefficients                t = (b1 at Z)/SE(b1 at Z)      df = N - k - 1      where _N_ =sample size and _k_ =number of predictors One can then test the significance of each slope using a t-test. My question is, how can I extend this to the significance of slopes in a 3-way interaction model? See my example R syntax below:               set.seed(123)     predictor <- rnorm(1000, 10, 5)     moderator1 <- rnorm(1000, 100, 25)     moderator2 <- rnorm(1000, 50, 20)     outcome <- predictor*moderator1*moderator2*rnorm(20, 30)/10000     mydata <- data.frame(predictor, moderator1, moderator2, outcome)          model <- lm(outcome ~ predictor + moderator1 + moderator2 + predictor*moderator1 + predictor*moderator2 + moderator1*moderator2 + predictor*moderator1*moderator2, data=mydata)     plotData <- expand.grid(                         predictor = pretty(qnorm(pnorm(c(-1, 1)), mean = mean(mydata$predictor, na.rm = TRUE), sd = sd(mydata$predictor, na.rm = TRUE))),                         moderator1 = qnorm(pnorm(c(-1, 1)), mean = mean(mydata$moderator1, na.rm = TRUE), sd = sd(mydata$moderator1, na.rm = TRUE)),                         moderator2 = qnorm(pnorm(c(-1, 1)), mean = mean(mydata$moderator2, na.rm = TRUE), sd = sd(mydata$moderator2, na.rm = TRUE))                         )          plotData$outcome <- predict(model, newdata = plotData, level = 0)         plotData$mod1 <- factor(plotData$moderator1, labels = c("Lo mod1", "Hi mod1"))         plotData$mod2 <- factor(plotData$moderator2, labels = c("Lo mod2", "Hi mod2"))          mod1Lo_mod2Lo <- plotData[plotData$mod1=="Lo mod1" & plotData$mod2=="Lo mod2",]     mod1Hi_mod2Lo <- plotData[plotData$mod1=="Hi mod1" & plotData$mod2=="Lo mod2",]     mod1Lo_mod2Hi <- plotData[plotData$mod1=="Lo mod1" & plotData$mod2=="Hi mod2",]     mod1Hi_mod2Hi <- plotData[plotData$mod1=="Hi mod1" & plotData$mod2=="Hi mod2",]          #Generate Plot     plot(mod1Lo_mod2Lo$predictor, mod1Lo_mod2Lo$outcome, lty=1, lwd=2, type='l', xlab="predictor", ylab="outcome", ylim=c(min(plotData$outcome), max(plotData$outcome)))     lines(mod1Hi_mod2Lo$predictor, mod1Hi_mod2Lo$outcome, lty=2, lwd=2)     lines(mod1Lo_mod2Hi$predictor, mod1Lo_mod2Hi$outcome, lty=1, lwd=2, col="gray")     lines(mod1Hi_mod2Hi$predictor, mod1Hi_mod2Hi$outcome, lty=2, lwd=2, col="gray")     legend("topleft", legend=c("lo Mod1, lo Mod2","hi Mod1, lo Mod2","lo Mod1, hi Mod2","hi Mod1, hi Mod2"), lty=c(1,2,1,2), lwd=c(2,2,2,2), col=c("black","black","grey","grey"))      How can I test whether the slopes of each of these four lines, separately, is different from zero? Many thanks in advance!


title: What is this distribution?

text: I have a right-skewed distribution that has this qq-plot against a normal with the same mean and standard deviation: ![QQ-Plot](http://i.stack.imgur.com/xQCh1.jpg) The data are time-delays between two events. Also, this is the histogram of the data: ![enter image description here](http://i.stack.imgur.com/nUS3m.jpg) and this is the QQ-plot of log(data + 1), vs. normal of the same mean and standard deviation as the log(data+1) values: ![enter image description here](http://i.stack.imgur.com/keTZU.jpg) I have tried fitting Weibull, lognormal, and Pareto, without success. Any thoughts on what would be the right candidate?


title: Balancing servers and users

text: Looking for some direction here. I have a model for a project where I need to balance a bunch of users amongst different servers based on their usage. For example, the heaviest users get placed on the lightest server. Problem a: I have a multiple data parameters for each server such as response time, disc size etc. I need to use these parameters to come with a score for the server. Example 1 to 5. 5 Being the server with most usage. What method can I use to weight the different parameters to come up with a server score? Problem b: Using those scores , I need to split the users amongst these servers. I already have a weighting score for each user, say 1 to 5. 1 being a light user and 5 being a heavy user. So I'm basically need to balance the servers, by placing user with a score of 1 in server with a score of 5. I am trying to use Excel solver to accomplish this. Any suggestions?


title: cointegration - same thing as stationary residuals?

text: So I'm aware that cointegration means there is some linear combination of the set of variables that is stationary. So, if you do a regression and find stationary residuals, can you just immediately say "oh, yes these variables are cointegrated" ? Because you have then just found some linear combination that is stationary? My concern is the oft cited reasons/bad consequences for why you shouldn't do a regression with trending/undifferenced/non-stationary variables. I won't repeat them here, but the one caveat is always "oh but you can regress these non-stationary variables if they are cointegrated". Is this the same as saying "go ahead and try regressions with any non-stationary variables you want, and if the residuals are stationary you can say, well it was ok to do this regression after all because these variables are cointegrated"? It just seems too simple, or too good to be true. Thanks in advance.


title: P value calculation

text: I have three groups of experiments. For each experiment I am looking for the percentage of occurrence of case x. In the first group I have 15 experiments. The case x was seen 10.191% of the total time for 15 experiments. In the second group I have 6 experiments. The percentage of x is 1.564%. In the third group I have 3 experiments. The percentage of x is 0%. I want to show that occurrence of case x significantly decreased from group one to two and to three. Thus I want to calculate the p-values. Can anyone tell me how to do it? Edit: The number of measurements for the first group is 22568 and 10.191% of these measurements are case x. The second group has 1854 measurements (1.564% are x) and the third group has 1164 measurements (0% x)


title: Stata - Multiple rotated plots on graph (including distributions on sides of axes)

text: I would like to produce a single graph containing both: (1) a scatter plot (2) either histograms or kernel density functions of the Y and X variables to the left of the Y axis and below the X axis. I found a graph that does this in MATLAB -- I would just like to produce something similar in Stata: ![enter image description here](http://i.stack.imgur.com/9mAjq.png) That graph was produced using the following MATLAB code:               n = 1000;     rho = .7;     Z = mvnrnd([0 0], [1 rho; rho 1], n);     U = normcdf(Z);     X = [gaminv(U(:,1),2,1) tinv(U(:,2),5)];          [n1,ctr1] = hist(X(:,1),20);     [n2,ctr2] = hist(X(:,2),20);     subplot(2,2,2); plot(X(:,1),X(:,2),'.'); axis([0 12 -8 8]); h1 = gca;     title('1000 Simulated Dependent t and Gamma Values');     xlabel('X1 ~ Gamma(2,1)'); ylabel('X2 ~ t(5)');     subplot(2,2,4); bar(ctr1,-n1,1); axis([0 12 -max(n1)*1.1 0]); axis('off'); h2 = gca;     subplot(2,2,1); barh(ctr2,-n2,1); axis([-max(n2)*1.1 0 -8 8]); axis('off'); h3 = gca;     set(h1,'Position',[0.35 0.35 0.55 0.55]);     set(h2,'Position',[.35 .1 .55 .15]);     set(h3,'Position',[.1 .35 .15 .55]);     colormap([.8 .8 1]);


title: Where does p-value come from in McNemar's Test?

text: I was using this online calculator to work out McNemar's Test on this data set: ![enter image description here](http://i.imgur.com/DFsziOx.png) I get the results: > The two-tailed P value equals 0.1859 > > Chi squared equals 1.750 with 1 degrees of freedom. Now I understand how to calculate the Chi squared but I have no idea where this P value comes from. Also how does the P value relate to whether there is a significant difference between classifiers? If someone could explain this like a program in JavaScript or something it would be so much clearer. Every explanation I see throws jargon all over the place.


title: Can I compare two regression coefficients

text: I am comparing treatment outcome to two therapeutic treatments. Specifically, I am looking at how attachment moderates the relationship between therapeutic alliance and outcome. I hypothesize that the moderator of attachment will be stronger in one treatment versus the other. The only way I can think to test this is to do two regression models. 1) alliance and outcome to Treatment one, as moderated by attachment. 2) alliance and outcome to treatment two, as moderated by attachment. Then, compare the regression coefficients to see if they are statistically different. Is there a way to put this all in one model?


title: How to obtain t-value for significance of median

text: I collected analysts’ forecasts and calculated bias and accuracy such as: $bias=actual-forecast$ $accuracy=abs(actual-forecast)$ What I am wondering is now on how to calculate the t-statistic to determine that the median bias is significantly different from zero. I am using R and found the t.test function to calculate mean significance but did not find anything appropriate for the median. (in my case I need the one sample t-test)


title: Variable selection using mixed-models (lme4)

text: I want to model the infection rates in bees based on weather conditions. The weather variables are rolling means for different time periods and durations. Dependent data is infection levels gathered in March and the independent variables are the weather aggregates (e.g. from 30 day period from Jan1-Jan30, 90 day period from Dec1-Feb28), a few thousands of them and highly correlated. PCA techniques did not work since the infections are not so strongly related to weather. I have also tried Bayesian Model Averaging and Boosted Regression Trees, since variables could be selected based on variable importance they calculate. But since, my data is longitudinal and my apiaries have a fixed location, I think mixed-models are a good choice. Is there a way to do variable selection based on mixed-models? What I have done now is to   1\. run `glmer` for each of the independent variables separately,   2\. remove those variable whose p-values for fixed-effect estimates are below 0.05 (not sure if this is a right thing - if the estimate for a variable is not significant, that variable being the only one in the model, it is right to drop that variable, is it?)   3\. from the variables that are left over, test for correlation between the variables   4\. remove the variables that are highly correlated, giving preference to the variable that has the lowest AIC. Or should I at this stage, not worry about p-values of Intercepts and only focus on AIC (or BIC)? since some of the variables have high p-values but AICs lower than than those with low p-values. I have tried reading up a lot, and there is no one fool-proof solution for variable selection, but would like to know if there is anything inherently wrong with my method. As I am not a statistician, equations often look like beautiful Arabic calligraphy and there lies my dead-end.


title: Possible reason for failing to build a support vector machine

text: I was trying to build a classifier for a set of documents using a support vector machine. I choose to build the feature space using term occurrence. While experimenting, I found the following scenario: When removing stop words, the svm-based classifier was successfully built; otherwise, when keeping stop words, the SVM just could not be built and I got an error message “no support vector can be found”. I am very confused about this scenario. What might be the possible reason for this scenario?


title: Sparsity-inducing regularization for stochastic matrices

text: It is well-known (e.g. in the field of compressive sensing) that the $L_1$ norm is "sparsity-inducing," in the sense that if we minimize the functional (for fixed matrix $A$ and vector $\vec{b}$) $$f_{A,\vec{b}}(\vec{x})=\|A\vec{x}-\vec{b}\|_2^2+\lambda\|\vec{x}\|_1$$ for large enough $\lambda>0$, we're likely for many choices of $A$, $\vec{b}$, and $\lambda$ to have many exactly-zero entries in the resulting $\vec{x}$. But if we minimize $f_{A,\vec{b}}$ subject to the condition that the entries of $\vec{x}$ are positive and sum to $1$, then the $L_1$ term doesn't have any effect (because $\|\vec{x}\|_1=1$ by fiat). Is there an analogous $L_1$-type regularizer that works in this case to encourage that the resulting $\vec{x}$ is sparse?


title: How to compare pre- and post-tests when subject identifier is missing?

text: My colleague and I are professors and we conducted an experiment in which we would please like some advice on deciding which tests to use in SPSS. He taught 4 classes of students on 2 different occasions and I taught a different group of 4 classes of students on 2 different occasions. I was the experimental group and taught all of my sessions incorporating an educational game. He was the control group and did not use any games. We each gave the students a pre-test at the beginning of the first session and a post-test at the end of the second class. Our hypothesis is that students in the experimental (games) classes performed better on the post-test than students in the control group. Unfortunately we didn't think to assign each student a number so that we could figure out which pre and post-test belonged to who. So basically we have a ton of pre and post- tests divided by class but not by student. Is there a a way we could conduct statistical analyses for the groups instead of individuals to see if our hypothesis was concerned?


title: Over-fitting in the cross-validation

text: I have a response variable(y) and 20 independent variables (Xs). I want to select several Xs in the linear regression, but I'm not sure how many variables should be selected. To select the best number of variables, I use the sum of the squared residuals (Res) in the 10-fold cross-validation given N selected variables (N=2~20). The process is repeated 1000 times given each N. My idea is that Res should firstly decrease as more variable could explain y better and then it should increase as too many variables should lead to over- fitting. To my surprise, Res decrease continually as N increase(see the Figure). I don't know how to explain it. Is it mean that all 20 variables contribute to y, or over-fitting happened? **P.S.:** there are about 600 data points. The Res is calculated as the sum of the square of the difference between observed y and predicted y in each 10-fold cross-validation. ![enter image description here](http://i.stack.imgur.com/mT7v7.png)


title: When does the prediction of random effects matter?

text: In linear or generalized linear mixed effects models, random effects are incorporated to explain the within-unit correlation for repeated measures over time. In Bayesian modeling, conventional prior choice for random effects are Gaussian, which is lack of theoretical justification and is chosen only for mathematical convenience. Despite its correctness, research suggests that the fixed effects are robust to the misspecification of the distribution of random effects, i.e. the estimations of fixed effects is not much affected by the choice of prior. So the choice of normal prior seems ok as long as as the random effects are not of interest. My question then is when is the prediction (since they are random variables not parameters) of random effects of interest? When should we prioritize modeling them accurately? And, since random effects are latent, what methods we can use to check if our prediction is adequate?


title: Using control variables in experiments?

text: **Why would one want to control for any number of baseline covariates in a situation where the assignment to treatment group is random?** My understanding is that randomly assigning treatment should make the treatment variable strictly exogenous, creating a control group that can appropriately be considered as a counterfactual. The only exception I can think of is when sample sizes are small, and that random assignment can still produce unbalanced groups. Any thoughts are much appreciated. Thanks!


title: Fastest way to run ridge regression on large datasets where n>>p

text: Provided that you don't want to do any variable selection: **Is there any software which is faster than`glmnet` at _vanilla ridge regression_ for large datasets?**


title: maximum gap between order statistics of normally distributed random variables

text: Hello Cross Validated community, I am currently working on a not-that-easy problem involving order statistics. As I am unsure as to how I could solve it, I thought it might already possess a solution. So here I am, my questions is: do you know of a pre-existing solution/approximation to (or do you see an easy way of solving) the probability below? $$\pi = \Pr\left[G > c\right],$$ where $$G_1 = \max_{1 \leq j \leq K} \left(Z_{(j)} - Z_{(j - 1)}\right), Z_{(0)} = 0, Z_j \sim N(\mu, \sigma^2) \text{ and } Z_i \perp Z_j, \forall i \ne j.$$ Just to clarify things up, $Z_{(j)}$ is the $j$th order statistics of $Z_1, Z_2, \ldots, Z_K$. Also, if it helps, the problem can actually be reformulated this way (but I would prefer the one above): $$G_2 = \max_{1 \leq j \leq K} \left(\tfrac{Z_{(j)}}{Z_{(j - 1)}}\right), Z_{(0)} = 1, Z_j \sim N(\mu, \sigma^2) \text{ and } Z_i \perp Z_j, \forall i \ne j.$$ I understand that mathematically speaking, those expressions are not equivalent, but for $\mu \in (0, 1)$ and $\sigma$ small enough, it should not matter. Sorry if the exact same problem was posted before, sorry for my english and thanks in advance! Max //=======================================================================// **ETA** : OK, worked quite a bit on this _damned_ problem, and found that for $\mu \in (0, 1)$, e.g. $\mu = 0.1$, and $\sigma$ small enough, e.g. $\sigma = 0.01$, $G$ actually simplifies to $$G = \max_{1 \leq j \leq n} \left(Z_{(j)} - Z_{(j - 1)}\right) \approx Z_{(1)}.$$ Thus, the probability $\pi$ is easily computed $$\pi = \Pr\left[G > c\right] \approx \Pr\left[Z_{(1)} > c\right] = {\left(1 - \Phi\left(\tfrac{c - \mu}{\sigma}\right)\right)}^K$$ The critical value $c$ is also easily extracted $${\left(1 - \Phi\left(\tfrac{c - \mu}{\sigma}\right)\right)}^K = \alpha \quad \Leftrightarrow \quad c = \mu + \Phi^{-1} \left(1 - \alpha^{\tfrac{1}{K}}\right)\sigma.$$ All I got to do now (other than finding a more general solution to the above problem) is find values of $\mu$ and $\sigma$, which are context-dependent. If anyone feels I should pile on and add some more details, please let me know. I am, after all, new to this forum. Cheers!


title: Law of total variance as Pythagorean theorem

text: Assume $X$ and $Y$ have finite second moment. In the Hilbert space of random variables with second finite moment (with inner product of $T_1,T_2$ defined by $E(T_1T_2)$, $||T||^2=E(T^2)$), we may interpret $E(Y|X)$ as the projection of $Y$ onto the space of functions of $X$. We also know that Law of Total Variance reads $$Var(Y)=E(Var(Y|X)) + Var(E(Y|X))$$ Is there a way to interpret this law in terms of the geometric picture above? I have been told that the law is the same as Pythagorean Theorem for the right-angled triangle with sides $Y, E(Y|X), Y-E(Y|X)$. I understand why the triangle is right-angled, but not how the Pythagorean Theorem is capturing the Law of Total Variance.


title: Regression Predictions

text: I don't think I fully explained my last question because it was marked as a duplicate. I went to the other question, and read it through several times, and it was asking something very different. Trying again: I have two variables that both influence Y, and I calculated regression equations for each. Based on the slope coefficients, am I able to determine if one is a better predictor of Y than the other, or do I have to find the correlation now to do that?


title: "Critical value at 5%" or "critical value at 95%": what is the right way to say it?

text: I have a quite simple (and not particularly intriguing) question about the way to present a test. I was wondering what is the right expression (or more widespread) when I talk about critical values. In particular, if for example the rejection region is defined by the 95 percentile, is it preferred to say that I am using the critical value of the, say, t "at 95%" or "at 5%"? Is is correct in both ways? Is one "more correct" than the other? Do I incur in some ambiguity if I do not specify further? Are there some conventions? Do they vary based on the distribution I am referring? Thank you very much


title: Experimental design and mixed models

text: I want to test effect of 3 PH on larval development. I would like to know what is the best experimental design and statistical analysis. We can only use 3 compartments of sea water, each one with a different ph. In each compartment, larvae will grow in 4 differents cylinder. Each day, we will take off severals larvae in each cylinder, measured larvae and count proportion of dead larvae. Can I then analyse proportion of dead larvae per day with a GLMM with cylinder as random factor (and LMM for larvae biometry) ? Can I "use" different parents to produce larvae of cylinders and analyse a kind of genetic effect ? In this case, could I use mixed models with family nested in cylinder as random effect ? ![experiment design](http://i.stack.imgur.com/Us3GR.png) Thank you.


title: Why is there not always a collinearity problem when introducing interaction terms to a MLR or logistic model?

text: Collinearity is described as occurring when an independent variable is a linear function of another variable(s). Such that the variable in question can be consistently predicted given the other variable(s). So given an interaction term (x1*x2) that is simply a function of x1 and x2, how is there not always a collinearity problem? I think that the answer is somewhere along the lines that even though x1*x2 is a function of x1 and x2, there may not necessarily be a linear dependency among the variables. And this is the part that I am having trouble conceptualizing. Does anyone have any examples of this or an R code that I can run to understand this better? Thanks, WC


title: How can I combine data from 2 separate experiments?

text: I am a scientist and by no means a statistician. I have, I think, a very basic question for you all. I have performed the same experiment twice (analyzing levels of a certain protein) and am comparing levels in untreated vs treated animals. I performed a t-test with welch correction (because there was significant variance). I am observing a similar trend, but neither experiment alone reaches statistical significance. Is there a way to combine the two experiments to increase power? Although the trend is similar, the actual values are different between the experiments, which is, I think, increasing my standard error. Also, if anybody knows how to perform this analysis using graphpad prism, that would be ideal. Thanks so much in advance for any assistance you may offer. Please let me know if you need more information.


title: Random Forest: Important variables, Important values

text: I've been reading up on random forests and have come to a stumbling block in regard to their practical application (an probably fundamental understanding). This problem i in regard to important variables and their use in prediction of new outcomes. Take the iris example. We know that based on the random forest that both petal length and petal width are the most important values to attain when trying to classify an iris species. Now, theoretically, we can go outside, measure some iris petal dimensions, plug them into predict(0 as a new test observation and the RF will generate an answer. To me, this suggests that the RF has probabilities (or probability like values) assigned to each value of iris petal dimensions. My question is: how do you retrieve these values? Is there a way I can extract a matrix or decision tree which says if you have petal width x and petal length y you likely have z species of iris? In essence, I want to know if it's possible to use the outcome of RF to classify new observations without having to run those new observations back through predict function.


title: Kruskal–Wallis Test not significant but Mann Whitney U significant

text: I have a non-normalized dependent variable and an independent variable broken down into 4 groups. As such, I used the Kruskal–Wallis analysis to look for significant differences in the ranks of the groups. The data look like the following: \begin{array}{clc}\rm Group&\rm Size&\rm Means\\\\\hline 0&n=24&31.79\\\ 1&n=13&26.65\\\ 2&n=8 &15.94\\\ 3&n=10&30.30\\\ \rm Total& N=55\end{array} I get an asymptotic significance of .103. However, if I run a Mann–Whitney _U_ on the same data, I see a clear significance between the `0` and `2` groups:                  0    24  18.81   451.50        2    8   9.56    76.50     Total   32                Mann-Whitney U  40.500     Wilcoxon W  76.500     Z   -2.416          Asymp. Sig. (2-tailed)  .016     Exact Sig. [2*(1-tailed Sig.)]  .013b      The only things I could come up with were that:   * I may be running into some kind of confounding error with the Mann–Whitney _U_ test, and thus the significance is a mistake, or that    * The large difference in between the sample sizes is what is causing me to see significance with the Mann–Whitney, but not the K–W since it has a bit more "power" to tease out what is significant and what isn't. Any help or guidance would be appreciated. You all have been of great help already! Cheers.


title: Calculating effect size variance from adjusted means (or least squared means)

text: I have a question related to this post regarding the calculation of the variance of an effect size. In my specific case, I would like to calculate the effect size and its variance from adjusted means (i.e., least squared means) derived from a mixed model. The output of the adjusted means (from the R function 'lsmeans') provides the standard errors of the adjusted means. However, the calculation of the effect size variance calls for the standard deviation: $\frac{(SD_{E})^2}{n_{E}\bar{X}_{E}^2}$ + $\frac{(SD_{C})^2}{n_{C}\bar{X}_{C}^2}$ In the answer to the above post, the commenter explains that the standard error is essentially the same as the standard deviation for an estimated parameter. Since I am using estimated parameters from a mixed model, would it be appropriate to use the standard errors in the above formula instead of the standard deviation?


title: Unsupervised pre-training for Reinforcement Learning

text: Since the advent of many unsupervised learning methods, as a pretraining step for the main supervised task (mostly under the name of Deep Learning), it shouldn't be a farto ask, what is the current state of "pretraining and learning from unlabeled data, for Reinforcement Learning"? Any recent/old works on this? Any suggestions for future work?


title: How to transform data to normality?

text: We have financial some data (500-1000 samples), which is not normally distributed (well known fact from the literature). I have some ideas to do parametric transformations of this data (using some other data) to produce "adjusted" series. My goal is to find a transformation that makes the series normally distributed (with mean 0 and std deviation 1). What is the most appropriate statistic and corresponding test to optimize my parameters and determine if the outcome can be considered normally distributed? Please also point me to an implementation, ideally in C/C++ or java.


title: Approximating a binomial distribution with a mixture normal

text: This is purely a theoretical question (I legitimately can't think of a real application), but if you wanted to approximate a binomial distributed variable with a two-component mixture normal, is there a closed form for the answer? In general, what are the preferred approaches for answering questions like that? I was looking at moment matching, but the algebra got hairy pretty quickly.


title: If the residuals of each group show same tendency, what does it mean?

text: I'm doing my reseach paper, and there are 7 groups for 5 years. I'm trying to do "panel gls" in STATA.   It's a little embarrassing to say, but I think I still don't have enough konwledge to analize the results.   I'm using STATA and it shows that every coefficient are all good (zero p-value).   But when I draw plot of residuals in EXCEL, it shows same shape to each group, not scattered but it looks like a heart rate.   Then what does this result mean and what should I do to remedy this?


title: How to visualize concurrent function blocks in a tree

text: I am developing a online interactive tool using d3.js and I am trying to visualize runtimes of functions as rectangular blocks and their calls to other functions in a tree structure. I am able to visualize the sequential functional flow but stuck when a parent function calls two functions at the same time(ie Concurrency). The blocks are overlapping and I would like to visualize the blocks wthout overlap and also display the underlying tree for each block. Could anyone please suggest me some existing algorithms to achieve this visualization.


title: Under what conditions should Likert scales be used as ordinal or interval data?

text: Many studies in the social sciences use Likert scales. When is it appropriate to use Likert data as ordinal and when is it appropriate to use it as interval data?


title: How to analyze Five point Likert scale data using regression

text: My research is on "Effectiveness on Audit Committee Practices". I have collected the data using Five point Likert scale (1=Strongly Disagree to 5= Strongly Agree). Both Independent and Dependent variables are Five point Likert scale data type. Also each variable include 2 or 3 questions with Five point Likert scale. I want to run a regression. So, can I use Principal component method to get values for variables? then is the regression linear? What are the other types of techniques to analyze these data. is regression sufficient?


title: Should Likert scale measures be treated as interval, ratio, ordinal, or nominal?

text: Forgive me, I am forgetting my basic statistics.


title: Working with correlation coefficients

text: I have three Pearson correlation coefficients (.8978, .5676 and .7865) for three age groups (i.e. 21 to 30 years, 31 to 40 years and 41 to 50 years) whose behaviour I am studying in regard to their shopping habits versus weight gain. Can I say .8978 is the strongest relationship between shopping habits and weight gain? Based on the diffference in the coefficients, can I say that there is a difference in the shopping habits and weight gain of the three age groups? Finally, can I just add the the three coefficients and divide by three to come up with an average? All these are 'face value' interpretations. Are they acceptable or do I need to perform some sort of statistical analysis? If it is the latter (heaven forbid), can SPSS do it.


title: How to test for differences between two group means when the data is not normally distributed?

text: I'll eliminate all the biological details and experiments and quote just the problem at hand and what I have done statistically. I would like to know if its right, and if not, how to proceed. If the data (or my explanation) isn't clear enough, I'll try to explain better by editing. Suppose I have two groups/observations, X and Y, with size $N_x=215$ and $N_y=40$. I would like to know if the means of these two observations are equal. My first question is:   1. If the assumptions are satisfied, is it relevant to use a parametric two-sample t-test here? I ask this because from my understanding its usually applied when the size is small?   2. I plotted histograms of both X and Y and they were not normally distributed, one of the assumptions of a two-sample t-test. My confusion is that, I consider them to be two populations and that's why I checked for normal distribution. But then I am about to perform a two-SAMPLE t-test... **Is this right?**   3. From central limit theorem, I understand that if you perform sampling (with/without repetition depending on your population size) multiple times and compute the average of the samples each time, then it will be approximately normally distributed. And, the mean of this random variables will be a good estimate of the population mean. So, I decided to do this on both X and Y, 1000 times, and obtained samples, and I assigned a random variable to the mean of each sample. The plot was very much normally distributed. The mean of X and Y were 4.2 and 15.8 (which were the same as population +- 0.15) and the variance was 0.95 and 12.11.   I performed a t-test on these two observations (1000 data points each) with unequal variances, because they are very different (0.95 and 12.11). And the null hypothesis was rejected.   Does this make sense at all? Is this correct / meaningful approach or a two- sample z-test is sufficient or its totally wrong?   4. I also performed a non-parametric Wilcoxon test just to be sure (on original X and Y) and the null hypothesis was convincingly rejected there as well. In the event that my previous method was utterly wrong, I suppose doing a non-parametric test is good, except for statistical power maybe? In both cases, the means were significantly different. However, I would like to know if either or both the approaches are faulty/totally wrong and if so, what is the alternative?


title: can you use mutual information to determine how much one variable tells you about another

text: I wrote a Naive Bayes Classifier and want to be able to test how much information one variable gives about another. The idea is to use the ones that are the most orthogonal and avoid using highly redundant variables. My idea was to use the mutual information, but have found that it is not as clear cut as I first thought. For example, if my variables are X = { 1, 2, 3, 4 } and Y = { 9, 10, 11, 12 }, then the mutual information would seem to be 2. But if X = { 1, 2, 3, 4, 5, 6, 7, 8 } and Y = { 9, 10, 11, 12, 13, 14, 15, 16 }, then the mutual information is 3. In both cases, however, if these numbers are lists (like columns in a database or spreadsheet), then X = 1 will always be paired with Y = 9 etc., meaning knowing one will always tell you the value of the other. So if both cases have a 1 to 1 mapping between X and Y and the mutual information is different, how useful is the mutual information in determining how much you know about one variable by knowing the other? Would some quantity like $$\frac{H(X,Y)}{H(X) + H(Y)}$$ or $$\frac{H(X,Y)}{I(X,Y)}$$ work better? Does anyone have experience with determining how much one variable will tell you about another? Here is a summary of the two examples. Example 1   * X = { 1, 2, 3, 4 }   * Y = { 9, 10, 11, 12 }   * H(X) = 2   * H(Y) = 2   * H(X,Y) = 2   * I(X,Y) = 2   * H(X,Y) / (H(X) + H(Y)) = 0.5   * H(X,Y) / I(X,Y) = 1.0 Example 2   * X = { 1, 2, 3, 4, 5, 6, 7, 8 }   * Y = { 9, 10, 11, 12, 13, 14, 15, 16 }   * H(X) = 3   * H(Y) = 3    * H(X,Y) = 3   * I(X,Y) = 3   * H(X,Y) / (H(X) + H(Y)) = 0.5      * H(X,Y) / I(X,Y) = 1.0


title: What is the range of values that can be expected in the result of Principal Component Analysis (PCA)?

text: I want to normalize all of my preprocessing techniques between 0 and 1 so I want to know what the PCA range of values is so that I can apply a proper normalization to it. I applied PCA by using the result of SVD. For example, I saw that the SVD range of values was limited to positive values, so I feature scaled to normalize the values between 0 and 1. Is there any such range for PCA?


title: Does it make sense to perform a one-tailed Kolmogorov-Smirnov test?

text: Is it meaningful and possible to perform a one-tailed KS test? What would the null hypothesis of such a test be? Or is the KS test inherently a two-tailed test? I would benefit from an answer that helped me understand the distribution of _D_ (I am working through Massey's 1951 paper, and find the description challenging, for example are $D^{+}$ and $D^{-}$ the supremum and infimum of the differences of the non-absolute value of differences in empirical CDFs?). Followup question: how are $p$-values for $D^{+}$ and $D^{-}$ obtained? So many of the publications I am encountering are presenting tabled values, rather than CDF of $D_{n}$, $D^{+}$ and $D^{-}$. **Update:** I just discovered the related question What's the null hypothesis in a one-sided Kolmogorov-Smirnov test?, which I missed on my initial scan before writing this one.


title: What is max-sum / max-product variant of loopy BP computing?

text: In (Nowazin and Lampert, _Structured Learning and Prediction in Computer Vision_, p. 29.), they say that in the max-sum variant of loopy belief propagation, the "variable max-beliefs are no longer interpretable as marginals but instead $\mu_i(y_i)$ describes the maximum negative energy achievable when fixing the variable $Y_i = y_i$." I understand that to mean that if we set $Y_i = y_i$, the maximum negative energy that the MRF can achieve across all variable assignments to other nodes (with $Y_i = y_i$ being fixed) is $\mu_i(y_i)$. Is that correct? In the next paragraph, they say: to recover a joint minimum energy labeling, we select for each variable $Y_i$ the state $y_i \in \mathcal{Y}_i$ with the maximum max-belief, $y_i^* = \text{argmax}_{y_i \in \mathcal{Y}_i} \mu_i(y_i) \text{ } \forall i \in V$. How is this a joint minimum energy labelling? I don't see how this chooses the $y_i$ such that their pair-wise interaction terms are considered. It looks like it considers each on its own, selecting the $y_i$ such that the selection of the other variables _could_ result in a maximum given that selection of $y_i$, but then it goes on to select the other $y$ locally as well, without consideration of this $y_i$. Is this "joint minimum energy labelling" not a global minimum energy labelling? If not what is meant by "joint"?


title: Proof that the bootstrap realizations are iid

text: Let $x_1,x_2, \ldots, x_n $ be the values of a random sample. A bootstrap sample $x_1 ^* ,x_2^*, \ldots, x_n^* $ is a random sample of $x_1, x_2, \ldots, x_n$ drawn with replacement. My question is how can I show that the bootstrap sample values are iid with common cdf $\hat{F_n}$, the empirical cdf of $x_1,x_2, \ldots, x_n $? Thank you.


title: How to use regression analysis to set an optimal price

text: I am working on a side project with very small dataset where i am trying to figure out the optimal price i should set for a transaction fee (something like payPal). Currently i am using an arbitrary number for the earning which is a fixed cost. The data i have:               Transaction_completed | value_of_transaction | my_earning      yes                         $10                 $1      no                          $100                $1      .                             .                  .      .                             .                  .      Was hoping to maybe used a more statistical approach to set an optimal price for the fee but I was not able to find any examples out there. Any regression ideas would be greatly appreciated.


title: Region comparability over a matrix sample

text: We are given a list of rows `R`, with elements `R_i`, `i=1..n`, and a list of columns `C`, with elements `C_j`, `j=1..n`, which build matrix `A[R, C]`. Each element of the matrix is either 0 or 1. Region is an axis-parallel rectangle from matrix `A`, whose boundaries are defined by two `R_i` and two `C_j` elements. A region has an associated weight. `weight = a * halfPerimeter + b * frequency`:   * `a, b` are given positive constants,    * `halfPerimeter` is sum of the number of rows from `R` and the number of columns from `C` which intersect the region, and   * `frequency` is the number of 1-elements within the region. We need to construct a sample matrix `A'[R',C']`, where list `R'` is a sample from `R`, and and `C'` is a sample from `C`. Both `R'` and `C'` are of the same size `r`, and the process of their building is symmetric, so we illustrate the process on `R'`. `R'` consists of elements `R'_i'`, `i' = 1..r`, such that   * `0 < r < n`, and    * For every `i'`, there exist `i` such that `R'_i' = R_i`. We say that row `R'_i'` correspond to `R_i`.    * _Same sort order._ For two rows `R'_i1'`, `R'_i2'`, if `i1' < i2'`,   then for their corresponding rows it holds that `R_i1 < R_i2`.   To achieve this, we sort after sampling. A region `reg'` in the matrix `A'` corresponds to the region `reg` in the matrix `A`, iff all the region boundaries (two rows and two columns) from `reg'` correspond to those of `reg`. **The problem** boils down to:   * fill the matrix `A'`, and    * define a weight function on it such that:     * if `reg1'` corresponds to `reg1` (`reg2'` corresponds to `reg2`), and   `weight(reg1') < weight(reg2')`,   then `weight(reg1) < weight(reg2)`     * `weight(reg') < C`, where `C` is a (relatively large) constant. I appreciate any help, suggestions, ideas, or pointers. Thanks!


title: MANOVA for count data

text: I'm comparing two groups where each individual was tested on 9 response variables. These variables are a count with an upper limit of 2 (it's a count of missing copies, and max number of copies is 2). I want to test if, considering all response variables as a whole, there are differences between both groups. MANOVA comes to mind as my first option, but my response variables are clearly not adequate for MANOVA. Is non-parametric version of MANOVA the best approach? I've tried with _vegan_ package in R but I'm not sure what distance method should I use (most of them give me warnings). Here are the data:               "group" "Y1" "Y2" "Y3" "Y4" "Y5" "Y6" "Y7" "Y8" "Y9"     "1" "A" 2 2 0 0 0 0 0 0 0     "2" "A" 1 1 0 0 1 1 0 0 0     "3" "A" 0 0 0 0 0 0 0 1 0     "4" "A" 0 0 0 0 0 0 0 0 0     "5" "A" 0 0 0 0 0 0 0 0 0     "6" "A" 0 0 0 0 0 0 0 0 0     "7" "A" 0 0 0 0 0 0 0 0 0     "8" "A" 2 2 0 1 0 1 0 0 1     "9" "A" 2 2 0 0 0 0 0 0 0     "10" "A" 2 2 0 0 0 0 0 0 0     "11" "A" 2 2 0 0 0 0 0 0 1     "12" "A" 2 2 0 1 0 1 0 0 0     "13" "A" 2 2 1 0 0 0 0 0 0     "14" "A" 2 2 0 1 0 0 0 0 0     "15" "A" 1 1 1 1 0 0 0 0 0     "16" "A" 1 1 0 1 0 0 0 0 0     "17" "A" 1 2 0 0 0 0 0 0 0     "18" "A" 1 2 1 0 0 0 1 0 0     "19" "A" 1 1 0 1 0 0 0 0 0     "20" "A" 1 1 0 0 0 0 0 0 0     "21" "A" 1 1 0 0 0 0 0 0 0     "22" "A" 1 1 0 0 0 0 0 0 0     "23" "A" 1 1 1 0 0 0 0 0 0     "24" "A" 1 1 0 0 0 0 0 0 0     "25" "A" 0 0 0 0 0 0 0 0 0     "26" "A" 0 0 0 0 0 0 0 0 0     "27" "A" 0 0 1 0 0 0 0 0 0     "28" "A" 0 0 1 0 0 0 0 0 0     "29" "A" 0 0 0 0 0 0 0 0 0     "30" "A" 0 0 0 0 0 0 0 0 0     "31" "A" 0 0 0 0 0 0 0 0 0     "32" "A" 0 0 0 0 0 0 0 0 0     "33" "A" 0 0 0 0 0 0 0 0 0     "34" "A" 0 0 0 0 0 0 0 0 0     "35" "A" 0 0 0 0 0 0 0 0 0     "36" "A" 0 0 1 0 0 0 0 0 0     "37" "A" 0 0 0 0 1 0 0 0 0     "38" "A" 0 0 0 0 0 0 0 0 0     "39" "A" 0 0 0 0 0 0 0 0 0     "40" "A" 0 0 0 0 0 0 0 0 0     "41" "A" 0 0 1 0 0 0 0 0 0     "42" "A" 0 0 0 0 0 0 0 0 0     "43" "A" 0 0 0 0 0 0 0 0 0     "44" "A" 0 0 0 1 0 0 0 0 0     "45" "A" 0 0 1 0 0 0 0 0 0     "46" "A" 0 0 0 0 0 0 0 0 0     "47" "A" 0 0 0 0 0 1 0 0 0     "48" "A" 0 0 1 1 0 0 0 0 0     "49" "A" 0 0 1 0 0 0 0 0 0     "50" "A" 0 0 0 0 0 0 0 0 0     "51" "A" 0 0 0 0 0 0 0 0 0     "52" "A" 0 0 0 0 0 0 0 0 0     "53" "A" 0 0 1 0 0 0 0 0 0     "54" "A" 0 0 0 0 0 0 0 0 0     "55" "A" 0 0 1 1 0 0 0 0 0     "56" "A" 0 0 0 0 0 0 0 0 0     "57" "A" 0 0 1 1 0 0 0 0 0     "58" "A" 0 0 0 0 0 0 0 0 0     "59" "A" 0 0 0 0 0 0 0 0 0     "60" "A" 0 0 0 0 0 0 0 0 0     "61" "A" 0 0 0 0 0 0 0 0 0     "62" "A" 0 0 0 0 1 1 0 0 0     "63" "A" 0 0 0 0 1 2 0 0 1     "64" "A" 0 0 0 0 0 0 0 0 0     "65" "A" 0 0 0 0 0 0 0 0 1     "66" "A" 0 0 1 0 0 0 0 0 0     "67" "A" 0 0 0 0 0 0 0 0 0     "68" "A" 0 0 0 0 0 0 0 0 1     "69" "A" 0 0 0 1 1 0 0 0 0     "70" "A" 0 0 0 1 0 0 0 0 1     "71" "A" 0 0 0 0 1 1 0 0 0     "72" "A" 0 0 0 0 0 0 0 0 0     "73" "A" 0 0 0 0 0 0 0 0 0     "74" "A" 0 0 0 1 0 0 0 0 0     "75" "A" 0 0 0 0 0 1 0 0 0     "76" "A" 0 0 1 0 0 0 0 0 1     "77" "A" 0 0 0 0 0 0 0 0 1     "78" "A" 0 0 0 0 0 0 0 0 0     "79" "A" 0 0 0 0 0 0 0 0 0     "80" "A" 0 0 0 0 1 1 0 0 0     "81" "A" 0 0 0 0 1 1 0 0 0     "82" "A" 0 0 0 0 0 0 0 0 0     "83" "A" 2 2 0 1 0 0 0 0 0     "84" "A" 2 2 0 1 0 0 1 0 0     "85" "A" 2 2 0 0 0 0 0 0 0     "86" "A" 2 2 0 0 0 0 0 0 0     "87" "A" 2 2 0 1 0 0 0 0 0     "88" "A" 1 1 0 1 0 0 0 0 0     "89" "A" 1 1 0 1 0 0 2 0 0     "90" "A" 1 1 0 0 0 0 0 0 0     "91" "A" 1 1 0 1 0 0 0 0 0     "92" "A" 1 1 0 0 0 0 0 0 0     "93" "A" 1 1 0 0 0 0 0 0 0     "94" "A" 0 0 0 0 0 0 1 0 0     "95" "A" 0 0 0 0 0 0 0 0 0     "96" "A" 0 0 0 0 0 1 0 0 0     "97" "A" 0 0 0 0 0 1 0 0 0     "98" "A" 0 0 1 1 0 0 0 0 0     "99" "A" 0 0 0 0 0 0 0 0 0     "100" "A" 0 0 1 1 1 0 0 0 0     "101" "A" 0 0 1 1 0 0 0 0 0     "102" "A" 0 0 0 0 0 0 0 0 0     "103" "A" 0 0 0 0 0 0 0 0 0     "104" "A" 0 0 0 0 0 0 0 0 1     "105" "A" 0 0 0 0 0 0 0 0 0     "106" "A" 0 0 0 1 0 0 0 0 0     "107" "A" 0 0 0 0 0 1 0 0 0     "108" "A" 0 0 0 0 1 1 0 0 0     "109" "A" 2 2 0 0 0 0 0 0 0     "110" "A" 2 2 0 0 0 0 0 0 0     "111" "A" 1 1 0 0 0 0 0 0 0     "112" "A" 2 2 0 0 0 0 0 0 0     "113" "A" 0 0 0 0 0 0 0 0 0     "114" "A" 0 0 0 0 0 0 0 0 0     "115" "A" 0 0 0 0 0 0 0 0 0     "116" "A" 2 2 0 0 0 0 0 0 0     "117" "A" 2 2 0 0 0 0 0 0 0     "118" "A" 2 2 0 1 0 0 0 0 0     "119" "A" 1 1 0 0 0 0 0 0 0     "120" "A" 1 1 0 0 0 1 0 0 0     "121" "A" 1 1 0 0 0 0 0 0 0     "122" "A" 0 0 0 0 0 0 0 0 0     "123" "A" 0 0 0 0 0 0 1 0 0     "124" "A" 0 0 0 0 0 0 0 0 0     "125" "A" 0 0 0 0 0 0 0 0 0     "126" "B" 1 0 0 0 0 0 0 0 0     "127" "B" 0 0 1 0 1 1 0 0 1     "128" "B" 0 0 0 0 0 0 0 0 0     "129" "B" 0 0 0 0 0 0 0 0 0     "130" "B" 0 0 0 0 0 0 0 0 0     "131" "B" 0 0 0 0 0 0 0 0 0     "132" "B" 0 0 0 0 0 1 1 0 0     "133" "B" 0 0 0 0 0 0 0 0 0     "134" "B" 0 0 0 0 0 0 0 0 0     "135" "B" 2 1 0 1 0 0 0 0 1     "136" "B" 1 2 0 1 0 0 0 0 0     "137" "B" 1 1 0 1 0 0 0 0 0     "138" "B" 0 0 0 0 0 0 0 0 1     "139" "B" 0 0 0 0 0 0 0 0 0     "140" "B" 0 0 0 0 0 0 0 0 0     "141" "B" 0 0 0 0 0 0 0 0 0     "142" "B" 0 0 0 0 0 0 0 0 1     "143" "B" 0 0 0 0 1 0 0 0 0     "144" "B" 0 0 0 0 0 0 0 0 0     "145" "B" 1 1 0 0 0 0 0 0 0     "146" "B" 0 0 0 1 0 0 0 0 0     "147" "B" 0 0 0 0 0 0 0 0 0     "148" "B" 2 1 0 1 0 0 0 0 1     "149" "B" 2 2 0 0 0 0 1 0 0     "150" "B" 2 0 0 0 1 0 0 0 0     "151" "B" 1 1 0 1 0 0 0 0 0     "152" "B" 1 1 0 1 0 0 0 1 1     "153" "B" 0 0 0 0 0 0 0 0 0     "154" "B" 0 0 0 0 1 0 0 0 1     "155" "B" 0 0 0 0 0 0 0 0 0     "156" "B" 0 0 0 0 0 0 0 0 0     "157" "B" 0 0 0 0 0 0 0 0 0     "158" "B" 0 0 0 0 0 0 0 0 0     "159" "B" 0 0 0 0 0 0 0 0 0     "160" "B" 0 0 0 1 0 2 0 1 1     "161" "B" 0 0 0 0 0 0 0 0 1     "162" "B" 0 0 0 0 0 0 0 0 0     "163" "B" 0 0 0 0 0 0 0 0 0     "164" "B" 0 0 0 0 0 0 0 0 0     "165" "B" 0 2 0 0 1 0 0 0 1     "166" "B" 0 0 0 0 0 0 0 0 0     "167" "B" 0 0 0 0 0 0 0 0 1     "168" "B" 0 0 0 0 0 0 1 0 0     "169" "B" 0 0 0 0 0 0 1 1 0     "170" "B" 0 0 0 0 0 0 0 0 0     "171" "B" 0 0 0 0 0 0 0 0 1     "172" "B" 0 0 0 0 0 1 0 1 1     "173" "B" 0 0 0 0 0 0 0 0 0     "174" "B" 0 0 0 0 0 0 0 0 0     "175" "B" 0 0 0 0 0 0 0 0 1     "176" "B" 0 0 0 0 0 0 0 1 1     "177" "B" 0 0 0 0 0 0 0 0 1     "178" "B" 0 0 0 0 0 0 0 0 0     "179" "B" 0 0 0 0 0 0 0 0 2     "180" "B" 0 0 0 0 0 0 0 1 1     "181" "B" 2 1 0 1 1 0 0 0 1     "182" "B" 2 2 0 0 0 0 1 0 0     "183" "B" 1 1 1 1 0 0 0 0 1     "184" "B" 1 1 0 1 0 1 0 0 0     "185" "B" 1 1 0 0 0 0 0 0 0     "186" "B" 1 1 0 0 0 2 1 0 1     "187" "B" 1 1 0 1 0 0 0 0 1     "188" "B" 1 1 0 0 0 0 0 0 1     "189" "B" 0 0 0 1 0 0 0 1 1     "190" "B" 0 0 0 0 0 2 0 0 0     "191" "B" 0 0 0 0 0 0 0 0 2     "192" "B" 0 0 0 0 0 0 0 0 0     "193" "B" 0 0 0 0 0 0 1 0 0     "194" "B" 0 0 0 0 0 0 0 0 0     "195" "B" 0 0 0 0 0 0 0 0 1     "196" "B" 0 0 0 0 0 0 0 0 0     "197" "B" 0 0 0 0 0 0 0 0 1     "198" "B" 0 0 0 0 0 0 1 0 0     "199" "B" 0 0 0 0 0 0 0 0 0     "200" "B" 0 0 0 0 0 0 0 0 0     "201" "B" 0 0 0 0 0 0 0 0 2     "202" "B" 0 0 0 0 0 0 0 0 2     "203" "B" 0 0 0 1 2 0 0 0 0     "204" "B" 0 0 0 0 0 0 0 0 0     "205" "B" 0 0 0 0 0 0 0 0 0     "206" "B" 0 0 0 0 0 0 0 0 1     "207" "B" 0 0 0 0 0 0 0 0 0     "208" "B" 0 0 0 0 0 0 0 0 0     "209" "B" 0 0 0 0 0 0 0 0 0     "210" "B" 0 0 0 1 0 0 0 0 2     "211" "B" 2 0 0 0 0 0 0 0 1     "212" "B" 2 0 0 0 0 0 0 0 0     "213" "B" 2 2 0 0 0 0 0 0 0     "214" "B" 1 0 1 0 1 0 1 0 1     "215" "B" 0 0 2 0 0 0 0 0 1     "216" "B" 0 0 1 0 0 0 0 0 1     "217" "B" 0 0 0 1 0 0 0 0 0     "218" "B" 0 0 0 0 0 1 0 0 1     "219" "B" 0 0 0 0 0 0 0 0 1     "220" "B" 1 0 0 0 0 0 0 0 0     "221" "B" 1 1 0 1 1 0 0 0 1     "222" "B" 1 2 0 1 0 0 0 0 1     "223" "B" 1 2 1 0 0 1 0 0 1     "224" "B" 1 1 0 1 0 0 0 0 1     "225" "B" 0 0 0 0 0 0 0 0 1     "226" "B" 0 0 0 0 0 0 0 0 0     "227" "B" 0 0 0 0 0 0 0 0 0     "228" "B" 0 0 0 1 0 2 0 0 0     "229" "B" 0 0 1 0 1 0 0 0 1     "230" "B" 2 2 0 1 0 0 0 0 1     "231" "B" 0 0 0 0 0 1 0 0 0     "232" "B" 0 0 0 0 0 0 0 0 0     "233" "B" 0 0 0 0 0 0 0 0 1     "234" "B" 0 0 0 1 1 0 0 0 1     "235" "B" 2 2 0 1 0 0 0 0 0     "236" "B" 2 2 0 1 2 0 0 0 1     "237" "B" 1 1 0 0 0 0 0 0 0     "238" "B" 0 0 0 0 0 0 0 0 0     "239" "B" 0 0 0 0 0 0 0 0 0     "240" "B" 0 0 0 0 0 0 0 0 0     "241" "B" 0 0 0 0 0 0 0 0 0     "242" "B" 2 2 1 0 0 0 0 0 0     "243" "B" 0 0 0 0 0 0 0 0 0     "244" "B" 0 0 0 0 0 0 0 0 0     "245" "B" 0 0 0 0 0 0 0 0 0     "246" "B" 0 0 0 0 0 0 0 0 1     "247" "B" 0 0 0 0 0 0 0 0 1     "248" "B" 0 0 0 0 0 0 0 0 0     "249" "B" 0 0 0 0 0 0 0 0 0     "250" "B" 1 0 0 0 0 1 0 0 0     "251" "B" 1 1 0 0 0 0 1 0 0     "252" "B" 0 0 0 0 0 0 0 0 0     "253" "B" 0 0 0 0 0 0 1 0 0     "254" "B" 0 0 0 0 0 0 0 0 0     "255" "B" 0 0 0 0 0 0 0 1 0     "256" "B" 0 0 0 0 0 0 0 0 0     "257" "B" 0 0 0 0 0 0 0 0 0     "258" "B" 2 2 0 0 0 0 0 0 0     "259" "B" 1 0 1 0 0 0 0 0 0     "260" "B" 0 0 0 0 0 0 0 0 0     "261" "B" 0 0 0 0 0 0 0 0 0     "262" "B" 2 2 0 1 0 0 0 0 0     "263" "B" 1 1 0 0 0 0 0 0 0     "264" "B" 0 0 0 0 0 1 0 0 0     "265" "B" 0 0 0 0 0 0 0 0 0     "266" "B" 2 1 0 0 0 2 0 0 0     "267" "B" 0 0 0 0 0 0 0 0 0     "268" "B" 0 0 0 0 0 0 0 0 0     "269" "B" 2 2 0 0 0 2 0 1 0     "270" "B" 2 2 0 0 0 2 0 0 0     "271" "B" 0 0 0 0 0 0 0 0 0     "272" "B" 2 2 0 0 0 0 0 0 0     "273" "B" 0 0 0 0 0 0 0 0 0     "274" "B" 2 2 0 1 0 0 1 0 0     "275" "B" 2 2 0 1 0 0 0 0 0     "276" "B" 0 0 0 0 0 0 0 0 0     "277" "B" 1 0 0 0 0 0 0 0 0     "278" "B" 0 0 0 0 0 0 0 0 0     "279" "B" 0 0 0 0 0 0 1 0 0     "280" "B" 1 0 0 0 0 0 0 0 0


title: Replacing outliers with mean

text: This question was asked by my friend who is not internet savvy. I've no statistics background and I've been searching around internet for this question. The question is : is it possible to replace outliers with mean value? if it's possible, is there any book reference/journals to backup this statement?


title: If using Glmm with Gamma distribution do i need to transform my data to be between 0 and 1?

text: Says it all in the title really. When creating a Glmm with Gamma distribution do i need to transform my response variable data to be between 0 and 1? Thanks Jonathan


title: Mathematical basis for conditional probability

text: In order to get more fundamental in my understanding of probability I watched mathematicalmonk's lectures involving $\sigma$-algebras etc. - good. One of my main concerns was to better understand the basis for conditional probability: $P[A|B] = P[AB]/P[B]$. My question is simple How do we know that this quotient is indeed a probability measure?   After all, things like the ``Odds'' $P[A] / P[\bar{A}]$ is not a probability measure.


title: Does Mahalanobis distances have "significance" associated with them?

text: I have a "distance matrix". let's say a 6x6 distance matrix, each cell is the Mahalanobis distance of two "clusters" (or sets/groups of things in a multidimensional space), I want to "count" the number of actual clusters (significantly separated clusters) using this distance matrix, My question is **if I only know the Mahalanobis distance between two clusters can I say anything about how well they are separated?** More information: these distances come from a Gaussian Mixture model fitting 6 Gaussians to a dataset that can have 1-4 actual clusters in it.


title: Non-nested model selection

text: Both the likelihood ratio test and the AIC are tools for choosing between two models and both are based on the log-likelihood. But, why the likelihood ratio test can't be used to choose between two non- nested models while AIC can?


title: Random effects structure of unbalanced design

text: The random effects structure of the following experiment puzzles me. I am interested in the consumption (quantity) differences between two types (A and B) of cattle forage. Therefore, I selected four type A and four type B products for testing. An additional hypothesis is that the raw protein content of each product is implicated, so I would like to add it as a covariate in the analysis. Thus, beforehand the raw protein in all eight products was measured. To quantify consumption, eight separately housed calves will test each of these eight products in one day after which amount consumed is measured. Each calve feeds on a single product only once (i.e. 24 hours), whereas each day all calves will feed on a different product – essentially following a Latin square design. However, trials involving two products failed as there were insufficient amounts of them due to delivery troubles (one of type A, one of type B). Therefore, these had to be excluded in the analysis. As each calve might have its own consumption speed, I’d like to include calve as a random effect. The consumption could also fluctuate on a daily basis, so time would seem another random effect. Exploratory data analysis revealed that consumption meets normality. Therefore, a linear mixed model would be best suited. In R, the model formula I tried is: consumption ~ foodType + protein + (1 + foodType | calve) + (1 | day) As each calve has its own feeding rate (random intercept) and might have its own preference for food types (random slope) I included the first random clause. The second random intercept for day is because consumption might change over time. However, I am unsure if the random effects should be nested, as per day each calve may have a different consumption. But, I cannot distinguish this because each day the calve feeds on a different product, so daily consumption per calve per day and the food product offered cannot be distinguished if I am correct. So, would this be the correct random effect structure to analyse this dataset?


title: Ridge regression – Bayesian interpretation

text: I have heard that ridge regression can be derived as the mean of a posterior distribution, if the prior is adequately chosen. Is the intuition that the constraints as set on the regression coefficients by the prior (e.g. standard normal distributions around 0) are identical / replace the penalty set on the squared size of the coefficients? Does the prior have to be Gaussian for this equivalence to hold?


title: Logistic regression with probability of success close to one

text: Are there any caveats when logistic regression is used on a sample with average probability of success close to one (1.4M dataset, mean prob. of success = 0.975)?


title: How to check a similar question in a survey

text: I have a survey in which I asked, "Would you be ready for a t-shirt shop in St Jean Avenue? [yes/no/don't know]". Then I asked many other questions, ultimately including, "Would you be ready for a t-shirt shop in St Jean Avenue knowing that the staff working there is disabled (mentally or physically)? [yes/no/don't know]" I would like to know:   1. What's the best way to visualize the two questions?   2. What's the best statistical test to perform? I would like to know if the people saying yes to the question 1 say yes to the question 2 and the people saying no to question 1 says no to question 2.


title: How to estimate pmf without knowing the entire sample space?

text: Suppose we estimate the probability mass function $p$ from a iid training sample $\\{x_1,\dots,x_n\\}$. Then we want to use the estimated pmf $\hat{p}$ to estimate the probabilities $p(y_j)$ for a iid test sample $\\{y_1, \dots,y_m\\}$ which are supposed to come from the same distribution as the training sample. However, we don't know the entire sample space when we estimate the pmf from the training sample. Specifically, we anticipate that there may be some test outcomes, $y_j$'s, which have not been seen in the training example, i.e. $y_j \neq x_i, i \neq 1,\dots, n$, and furthermore, we don't know what $y_j$'s are and how many of them. My question is how we can estimate the pmf from the training sample to deal with such unknown testing outcomes? Note that, this is different from the case when we know the entire sample space, but not observe some outcomes. In this case, we can use Laplace smoothing to move probability mass to those with zero occurrences. The question arises from estimating the distributions for N-gram models in NLP. Thanks!


title: Problems with acceptance ratio of MCMC when using a Kalman filter and conjugate-prior for sampling

text: ### 1\. Model I am trying to build a MCMC estimation of the following model (simplified): $\log(P^{-1}(obs_t, \sigma)) = \log(Y_t) + \epsilon_t$ where $\epsilon_t \sim \mathcal{N}(0,\sigma_{e})$. $P^{-1}$ is the inverse of a pricing function, that connects an observed price to an unobserved variable, for which I assume the error $\epsilon$. Apart from the observed price it depends on a parameter $\sigma$ (I omit other parameters for simplicity). Further, the unobserved variable follows a Geometric Brownian Motion, i.e. $\log(Y_t) = \log(Y_{t-1}) + \mu - \frac{1}{2}\sigma_Y^2 + \theta_t$ where $\theta_t \sim \mathcal{N}(0,\sigma_Y)$. ### 2\. Estimation Approach _Sampling $Y_t$:_ I use the Forward Filtering Backward Sampling method to generate samples of all $Y_t$: After calculating $P^{-1}$ i run a Kalman filter on the $Y_t$ and draw samples from the Normal Distribution starting in $T$ afterwards treating the already drawn samples as additional observations. [1] _Samping $\sigma_Y$:_ Starting from an Inverse Gamma prior I apply the sampled $Y_t$ using the conjugate prior (standard Bayesian regression) and draw a sample from the posterior Inverse Gamma distribution. I then calculate the likelihood of the $Y^*_t$ implied by the sampled $\sigma_Y$ (as according to the first model equation) and specifiy a Metroplis-Hastings acceptance ratio, accepting the sample or not. ### 3\. Problem I am testing the estimation on a simulated model and experience the following problem:   * The parameter $\sigma_Y$ does not influence the volatilty of the $Y_t$ implied by the inverse pricing function $P^{-1}$ much   * Hence the sample from the Inverse Gamma distribution is quite near to the "true" value of the parameter    * However $\sigma_Y$ does influence the level of the $Y_t$ strongly   * Therefore the likelihood of the $Y^*_t$ is very small and the sample is not accepted I think the likelihood would be higher, if the variance used to draw the samples of $Y_t$ were higher in the Kalman filter. However the variance is implied by both $\sigma_Y$ and $\sigma_e$ so I see no way to influence it (value of $\sigma_e$ seems to be reasonable). ### 4\. Questions   * Is there any conceptual flaw in my approach which could lead to the problem?   * Is ther perhaps a technique to improve the acceptance ratio in such a case? I suppose moving from the intial value of $\sigma_Y$ to the "true" value more slowly would help. But I don't see how I could achieve that... So, if you are still in at this point thanks for your patience! Any ideas welcome... Best, Matt [1] Carter, Chris K., and Robert Kohn. "On Gibbs sampling for state space models." Biometrika 81.3 (1994): 541-553.


title: Predicting count data with ANN

text: I am trying to model a count dependent variable (0~60) with excessive zeros. The objective is prediction. I've tried some parametric statistical method such as Poisson Regression, Negative-Binomial Regression, Zero-Inflated and Hurdle model. All of them do not give good predicts (RMSE ~= 5) I also tried ANN and the result is pretty good (RMSE ~= 1). But it's just the standard ANN so far, which normally is applied to a continuous DV. I am just wondering that using ANN for this problem is "awkward" or not? What method should work best to predict count DV ? Thank you very much. Any help is appreciated ! Cheers, Leo


title: Countable intersection of almost sure events is also almost sure

text: Suppose that $\\{B_n\\}_{n\in\mathbb N}$ are almost sure events, i.e. $\mathbb P(B_n)=1, \forall n$. Then how do I show that their intersection is also a sure event, i.e. $\mathbb P\left(\cap_{n=1}^\infty B_n\right)=1$. Thanks!


title: How to pedict spatio-temporal events?

text: My data set consist of a large set of events. Each event occurs at an arbitrary location and lasts for an arbitrary but limited amount of time. Since I have no experience with spatio-temporal modeling, can you give my a starting point on how to build a model for such data which allows me to predict where events are likely to occur in the future and how long they will likely last? Is this possible after all?


title: Textbooks pertaining to creating models?

text: I'm not sure if a question like this has been asked yet on this website. I recently graduated as an undergraduate statistics major with a bit of a dense math background (I know some probability theory and stochastic calculus and took two semesters of real analysis and abstract algebra). I'm a bit disappointed, since despite all of the theory I have learned, I don't have too much background in terms of applying what I've learned. The only experience I have with applications is an introductory-level econometrics course, which strangely enough was not required for the statistics major (I took it back when I was still an actuarial major). I found that class to be a bore with the lack of mathematical rigor in it, and I asked the professor for a proof of the Gauss-Markov Theorem while I took the course. For someone with my mathematical background, what would you recommend as mathematically-oriented textbooks for learning about econometrics and issues that come with creating models (multicollinearity, heteroskedasticity, serial correlation, etc.)? (If my question is vaguely worded, let me know, and I can edit it to be more specific.)


title: How to predict extremes of a continuous response variable?

text: Suppose I have a response $Y_t \sim N(0,\sigma^2_Y)$ and features $X_{i,t} \sim N(0,\sigma^2_{X_i})$ for $i \in \\{1,...,100\\}$. $Y_t$ conforms to a _linear model_ $Y_t = a + \sum_{i=1}^{100} b_i X_{i,t} + e_t$ My objectives are:   * Build an accurate out of sample prediction model for $Y_t$ based on a regularised linear model with all $X_{i,t}$.   * Build a model with extra precise accuracy specifically at the extremes of $Y_t$. That is, out of sample errors are very small when $Y_t$ approaches very negative or very positive values (relative to $N(0,\sigma^2_Y)$).   * I am interested in practical solutions.   * I do not care about standard error unbiasedness/consistency. How do I achieve this in the context of regularised linear models? I realise I need some kind of asymmetrically weighted loss function. Has there been work on this? Alternatively perhaps, I could forecast $sign(Y_t)*|Y_t|^c$, where $c > 1$, in order to weight the extremes more heavily, and choose $c$ through causal cross-validation based on a utility function over the errors at the extremes (optimise $c$ based on $\sum_{t=T+1}^{T+T'} U(error(Y_t,e_t))$?


title: How to normalize data to 0-1 range?

text: I am lost in normalizing, could anyone guide me please. I have a minimum and maximum values, say -23.89 and 7.54990767, respectively. If I get a value of 5.6878 how can I scale this value on a scale of 0 to 1.


title: Normalizing difference between two real values to [0,1] interval

text: If I have two positive real numbers that can take on any value between 0 and some finite real number, how do I normalized the difference between these two numbers to [0,1] interval where 0 indicates no similarity and 1 indicates that these two numbers are the same? Thank you in advance.


title: How to check if a random sample is biased

text: A question on how to prove that differences in mean value are statictically significant, and not just random noise. I have a set of two observations, one of which I'm going to deliberately bias like this (in R):               m <- cbind( A=sample(0:10, 10, replace=TRUE), B=sample(1:11, 10, replace=TRUE) )      Now, even with such a small data set, the mean of B is usually (but not always) higher than that of A.               colMeans(m)      As I accumulate more measurements, the bias should be more clearly revealed. Assuming I didn't know about the bias beforehand, how can I _prove_ it mathematically (in R if possible)?


title: Specifying the validation dataset for liblinear

text: I am trying to use the liblinear logistic regression model with L2 regularization. I don't want the training data to be splitted for the cross validation. I want to specify my own validation set for the validation purpose. Can anyone let me know how to do it?