title: Multigroup Longitudinal Item Response WinBUGS OpenBUGS

text: I was able to fit the Longitudinal IRT model in Winbugs for an ordinal response by extending the BUGS code I took from the paper by Curtis in JSS http://www.jstatsoft.org/v36/c01/paper/ However, I am having difficulty introducing a group effect into the code. The data are 7-category ordinal responses coded 1-7, and measured at multiple time points for each subject. I would like to compare 3 groups of subjects on their latent variable (ability parameter). Here is the BUGs code (I am using R2WinBUGS) without the group effect and it runs as expected:               lgrmar1 <- function(){             #  Graded response longitudinal model             #  Unconstrained - discrimination parameters alpha freely estimated         for (t in 1:T){         for (i in 1:n){             for (j in 1:p){                 Y[i, j, t] ~ dcat(prob[i, j, t, 1:K[j]])                      for (k in 1:(K[j] - 1)) {                 logit(P[i, j, t, k]) <- kappa[j, k] - alpha[j]*theta[i, t]                  }                 P[i, j, t, K[j]] <- 1.0             }                   for (j in 1:p){                 prob[i, j, t, 1] <- P[i, j, t, 1]                 for (k in 2:K[j]) {                     prob[i, j, t, k] <- P[i, j, t, k] - P[i, j, t, k - 1]                 }             }         }     }          for (i in 1:n){         theta[i, 1:T] ~ dmnorm(mu.theta[], Pr.theta[,])     }          # Prior for mu.theta     mu.theta[1] <- 0.0     for (t in 2:T){         mu.theta[t] ~ dnorm(m.mu.theta, pr.mu.theta)     }          pr.mu.theta <- pow(s.mu.theta, -2)          # AR(1) structure for Sigma.theta     sigsq.theta <- 1.0     Sigma.theta[1, 1] <- sigsq.theta     for (t in 2:T){         Sigma.theta[t, t] <- sigsq.theta         for (j in 1:(t - 1)){             Sigma.theta[t, j] <- sigsq.theta*pow(rho, t - j)             Sigma.theta[j, t] <- Sigma.theta[t, j]         }     }          Pr.theta[1:T, 1:T] <- inverse(Sigma.theta[,])     rho ~ dunif(-1.0, 1.0)          # Priors on item parameters     for (j in 1:p){         alpha[j] ~ dnorm(m.alpha, pr.alpha) %_% I(0, )     }     pr.alpha <- pow(s.alpha, -2)          # thresholds need to be ranked     for (j in 1:p){         for (k in 1:K[j] - 1){             kappa.star [j, k] ~ dnorm(m.kappa, pr.kappa)             kappa[j, k] <- ranked(kappa.star[j, 1:(K[j] - 1)], k)         }     }     pr.kappa <- pow(s.kappa, -2)            }          Constants are:      time      T <- dim(Y)[3]      n.variables      p <- dim(Y)[2]      sample size      n <- dim(Y)[1]       K - max level for ordinal variable for each item      K <- apply(apply(Y, c(2,3), max), 1, max)          Prior parameters:      m.alpha <- 1.0      s.alpha <- 2.5      m.kappa <- 0.0     s.kappa <- 2.5      m.mu.theta <- 0      s.mu.theta <- 1      * * * My question is, how can I code a group effect on the latent variable scale? Should I model the group effect directly in the likelihood statement (say, logit(P[i, j, t, k]) <\- kappa[j, k] - alpha[j]*theta[i, t] + coef1*grp1 + coef2*grp2, there are 3 groups), or should the group effect enter though the latent variable theta? For the latter case, what would be the modification to the code?


title: How/Why does resampling from "any" distribution lead to a normal distribution?

text: I was performing some Monte-Carlos on historical data and irrespective of the distribution of the data I would always get a normal distribution owing to resampling with replacement. That made it easy for me to predict with 95% confidence what the expected value of that 'variable' would be. So far so good and so cool! No matter what the historical distribution of the variable looked like, resampling and estimating future probability of occurrence always seemed to follow a normal distribution. Now, normal distribution is not so normal in practice. So what's the phenomenon that leads to a normal distribution? Is there a mathematical proof for it? I'm sure it has something to do with the central limit theorem but I'm quite baffled and intrigued at the beauty of producing a normal distribution when resampling with replacement. I may be incorrect but is this true in general? Irrespective of my historical distribution (whether, beta, poisson, binomial, random etc.) I keep getting a normal distribution on resampling. Any help on the mathematics underpinning this phenomenon would be helpful.


title: Convergence of value iteration

text: Why is the termination condition of the value-iteration algorithm ( example http://aima-java.googlecode.com/svn/trunk/aima- core/src/main/java/aima/core/probability/mdp/search/ValueIteration.java ) as it is? In the MDP (Markov Decision Process) we have $||U_{i+1}-U_i||< \text{error}\cdot(1-\gamma)/\gamma$, where $U_i$ is a vector of utilities   $U_{i+1}$ is the vector of updated utilities   $\text{error}$ is the error bound used in the algorithm   $\gamma$ is the discount factor used in the algorithm   * Where does "$\text{error}\cdot(1-\gamma)/\gamma$" come from?    * Is the term "$/\gamma$" because every step is discounted by $\gamma$? But then what about $\text{error}\cdot(1-\gamma)$?    * And how big must $\text{error}$ be?


title: Rate of convergence for SLLN

text: I am interested in writing a non-asymptotic rate of convergence for SLLN as a function of number of samples. From the literature I've read so far, CLT provides an asymptotic convergence rate of $(1/\sqrt N)$ for SLLN. Also, Berry-Esseen provides a non-asymptotic bound in terms of c.d.f's. $$|F_N(x) - \Phi(x)| \le \frac{C\mathbb{E}(|X|^3)}{\sigma^3\sqrt N}$$ Is there a Berry-Esseen like statement to bound the difference between sample mean and the expected value of the underlying distribution as a function of N (number of samples)?


title: Implementing a vector autoregressive model where only one variable is of interest

text: I have a dataset of 100 different time series and I am trying to forecast only one of them. However I think that the 99 other time series influence the one that I am interested in so I use a VAR model so that time series that I am interested in is a linear combination of its own lagged values and the lagged values of the other 99 time series'. The intuition is that the added information from the other time series will improve the forecast performance (i.e. it will be better than the forecast performance of a simple AR model where only the lagged values of the time series of interest are included). So my first question is whether this intuition is correct. My second question is if there are any known methods to determine which time-series will improve forecast performance of the time series of interest and which will not (similar to how feature selection works).


title: VARMAX model in R

text: Is there a function in R that estimates the VARMAX model? There is one for a VARX (MTS package), but I didn't find one that works with the MA part also...


title: R predict with "prediction" option

text: I'm trying to match the prediction interval of the predict.lm() function in R using the formula found in this discussion : Obtaining a formula for prediction limits in a linear model I'm using a student's quantile in my interval but in the end it's far larger from the one given by predict(). Is there any specific calculation in the predict function, I tried to look at the code but couldn't find any answer. The formula looks ok as I found exactly the same from others source. My R code :               airquality_clean <- na.omit(airquality)     attach(airquality_clean)          #Model estimation     model_1 <- lm(Ozone ~., data = airquality_clean)          #Unbias variance of the residuals     sigma_2 <- sum(model_1$residuals**2)/(dim(airquality_clean)[1]-dim(airquality_clean)[2])          #New observation     new <- data.frame(Solar.R=200,Wind=10,Temp=70,Day=1,Month=3)          #Calculated prediction interval     sigma <- sqrt(sigma_2*(1 + as.matrix(new)%*%solve(as.matrix(t(airquality_clean[,-1]))%*%as.matrix(airquality_clean[,-1]))%*%as.matrix(t(new))))     qt <- qt(0.995, df = dim(airquality_clean)[1]-dim(airquality_clean)[2])     int_pred_t <- cbind(predict(model_1, new)-(qt*sigma),predict(model_1, new)+(qt*sigma))     int_pred_t               [,1]     [,2]     [1,] -22.59931 95.82563          #R prediction interval     predict(model_1, new, interval="predict", level=0.99)}            fit       lwr      upr     1 36.61316 -21.12916 94.35548      I'm not too far but it's not the same results. If I use a p value from a normal distribution and not a student I'm even closer. Thank you.


title: Compare Frequencies

text: I have a continuous predictor variable that I have divided into bins (at 0%, 0 to 25%, ...). For each of these bins, I have shown you the histogram below of Number of Companies that are in that range and the number of bankruptcies. You'd expect bankruptcies to rise as asset growth reaches 100% (and also when asset growth is near 0%). An overall chi-square test of independence on these bins revealed that bankruptcies and asset growth% bins are not independent. How do I prove that higher ( **and/or** lower) asset growth are strongly related to bankruptcies? Of course, the hit rate for asset growth% is the highest. I used Chi-Sq on each bin (df=1) and expected bankruptcy = (avg hit rate) x bin size but this may be incorrect. Do I perform univariate logistic regression 6 times and compare model results? (all prob are=0 for the asset growth=100% bin using LR/glmfit). I use **MATLAB**. ![Histogram for Number of Companies & Bankruptcies over asset growth %](http://i.stack.imgur.com/MVqvu.jpg)


title: Why are diagnostics based on residuals?

text: In simple linear regression one often wants to verify if certain assumptions are met to be able to do inference (e.g. residuals are normally distributed). Is it reasonable to check the assumptions by checking if the fitted values are normally distributed?


title: Coefficient estimated with a binary predictor $\in \{0,1\}$, but making predictions with values between $0$ and $1$ - is this OK?

text: Let's say I have a variable $x_d$ that, in the estimation data, is a simple indicator ($x_d \in \left\\{0,1\right\\}$). I estimate a coefficient for it, $\beta_d$, along with several other coefficients for variables that could be continuous, more dummies, etc. Now someone wants to use these coefficients in a forecast. For their input data, they would like to treat $x_d$ not as an indicator, but rather as a continuous variable. That is, while the estimation data represented this quantity as an either/or, the forecast data represents it more as a matter of degree. In this forecast data, I would hope that my colleague would constrain the value of $x_d$ between 0 and 1, but we have not yet discussed that. ### Edit $x_d$ describes whether a loan came from a very specific lender. So in reality, the value really is only 0 or 1. In the forecast, my colleague wants to treat that variable as a continuous measure of similarity to that lender, i.e. how much like that lender is a loan’s source, with 0 meaning not at all, and 1 meaning the same lender. Is it a bad thing to alter the meaning in that way? I’m especially concerned about how it affects the meanings of the other coefficients, which were estimated alongside an indicator and are now being used alongside a continuous variable. Are there any issues when a coefficient is estimated with an indicator/dummy, but then used in a forecast with a continuous variable?


title: How does eigenvalues work with binary data in redundancy analysis?

text: I am using the vegan package in R to do a redundancy analysis (RDA, a part of canonical correlation analysis). My response data is binary and my explanatory variables contains 0, 0.5 an 1. I get quite low eigenvalues (~0.05) and my question is how is the binary data affecting the eigenvalue? Will the variability always be 'badly' explained?


title: I'm not sure that statsmodels is predicting out-of-sample

text: I am fitting an ARMA model to my data and here is my code               import statsmodels.tsa.arima_model as ari          model=ari.ARMA(pivoted['price'],(2,1))     ar_res=model.fit()     preds=ar_res.predict(100,400)      What I want is to train the ARMA model up to the 100th data point and then test out-of-sample on the 100-400th data points. But I don't think that is what's happening. I think this is because the model is fit on the entire data set, so how can the predictions be out-of-sample? It says here that out-of- sample prediction is possible. Can someone please help me understand this? OK here is my new code               model=ari.ARMA(pivoted['price'][:100],order=(2,1))     ar_res=model.fit()     preds=ar_res.predict(100,400)      however it return the error               AssertionError: Index length did not match values


title: Measure theory interpretation - textbook advice

text: I have seen the other questions about measure theory book advice, and I don't think any of them fit what I am looking for. The vast majority of measure theory textbooks are (naturally) math-based and discuss different measures. What I am looking for is a measure theory book that focuses more on the _interpretation_ and how properties of things like the variance/standard deviation would change if you modify their underlying measures upon which they're based. In other words, a book that discusses interpretation, how to judge metrics (e.g. is this sensitive to outliers, is it robust, etc.) and how to make a metric , etc. For example, the accepted centrality measure is the mean, not mode or median. So if I would like to know how the performance of the variance/standard deviation changes if I replaced the mean with the median, how would I do that? In other words, how would I know to conclude that "it doesn't make much difference", if I replace the mean with the median in the variance calculation? Thank you for your help. **EDIT** : Is this question really a decision theory question?


title: Can you use the chi-squared test with tables of real numbers?

text: Consider a $N \times 2$ matrix (MATLAB notation):               M = [1.2, 3;           1.4, 2;           1.8, 1;           2.0, 2];      That is a $4 \times 2$ matrix where the first column consists of reals and the second of integers. My question is simply: does it make sense to run a $\chi^2$ test on this and if so how would you go about doing it? I am inclined to say no as tallying frequencies for reals (I have here truncated to a single decimal point, but this could be 5 or 6 decimal points) does not make sense. Am I way off, on or in between?


title: SPSS Date Correlation

text: I have 2 sets of data: one is a list of all days in 2012 in date format (`DatesA`); the other is a list of dates in 2012 on which certain events occurred (`DatesB`). Some dates in `DatesB` appear more than once because multiple events occurred on those dates. I want to see how closely the dates in `DatesB` correlate to any given single date.   What is the process for achieving this?


title: Time series modeling with independent variables

text: Let's say that I have data depicting the number of museum visits per day. My challenge is to understand how certain external (exogenous?) variables, e.g., weather and advertizing affect the number of daily visits. Moreover, I need to predict how the number of visits will change when the external variables are changed. My data has some quite obvious seasonal fluctuations where weekends and certain holidays attract more visitors. A yearly seasonal length can also be assumed to exist. It has been quite some time since I last worked with time series, and need therefore to be pointed in the right direction. What is the approach I should use to model this data? Pre-processing steps? Model type? Typical pitfalls? Does anybody know of worked examples I could get inspiration from?


title: How to calculate this kind of kappa?

text: A colleague and I are trying to annotate an alignment task. We have two documents. The first document is the original copy and the other document is a modified version. The other annotator and I are trying to align the sentences in the second copy to the original. A sentence in the second document could be the same as in the original; a sentence could also be added; and a sentence in the original could also be deleted. Also, one sentence could be aligned to multiple sentences in the original copy. And multiple sentences could be aligned to one sentence. ======= Here is an example: The first document has 10 sentences: Sent-1, Sent-2, Sent-3, Sent-4, ... , Sent-10. The second document has 9 sentences: Sent-1, Sent-2, Sent-3, ... , Sent-9. My annotation is: 1, 2, 3, 4, ADD, 5, 6, 8, 10 It shows that I think Sent-1, 2, 3,4 are aligned to 1,2,3,4 in the original sentence, Sent-5 is added, Sent-6 is aligned to 5, and Sent-7 to 6, Sent-8 to 8, Sent-9 to 10. The other unaligned ones are deleted. The other annotator has: 1, 2, 3, 4, 5, 5, 6, 9, 10 This shows that the other annotator has a different opinion on sentence 5, and sentence 8. He thinks 5 should be aligned to 5 and 8 should be aligned to 9. So what is the kappa here?


title: Difference between homoscedastic and homogenic-variance data

text: I have a doubt regarding two tests: **Breusch–Pagan test** , to detect heteroscedasticity in a series, and **Bartlett's test** , to test for equal variances for samples from $k$ populations. What are the difference betwen those two tests?   * Are not those tests highly related?   * How could a homoscedastic data give no homogeneity in the variance?


title: Testing statistical software

text: What techniques/approaches are useful in testing statistical software? I'm particularly interested in programs that do parametric estimation using maximum likelihood. Comparing results to those from other programs or published sources is not always possible since most of the time when I write a program of my own it is because the computation I need is not already implemented in an existing system. I am not insisting on methods which can guarantee correctness. I would be happy with techniques that can catch some fraction of errors.


title: ANOVA when assumption of independence is questionable/violated

text: This is a tricky one to phrase, but I am wondering how to interpret the assumption of independence requirement for ANOVA in situations in which the overall 'space' is fixed. For example: Imagine you have a competition for a total prize fund of $1M between several different individuals each of whom can be assigned to a personality type A, B and C. The hypothesis is that one of these types is better than the others at securing money in the competition. so you run the competition, find out how much each person in each group has gained in the competition and want to test the hypothesis by comparing the means of each group. A simple one-way ANOVA seems the obvious choice. But in this situation none of the observations is independent of the others, as with a fixed sum available a gain by one player implies a loss by another. And this is as true within groups as it is between them. With only three players it's obviously ridiculous to do any sort of statistical analysis. With an infinite number of players the dependence of one person's outcome on the others is infinitely small which suggests that an ANOVA would be in principle fine. And somewhere in between is the problem. Is it ok to test for differences in means with ANOVA in situations like these? Having thought it through I think the answer is yes but only if you keep an eye on how the situation will manifest itself in the final test statistics. Let's assume that normality and equality of variance assumptions can be met, but how dependent of each other do events have to be before reaching for some other technique? And what might be a better approach in that case?


title: On how to formulate and apply maximum likelihood

text: I have just delved into the basics of maximum likelihood estimation and expectation maximization. The latter is really difficult to follow and I am having a tough time in figuring how I can apply the EM method for parameter estimation.


title: Covariance matrix specification in multivariate probit

text: Im having trouble with a multivariate probit model with partial observability/sample selection (written in GAUSS). In this model there is a probit at each of multiple stages, and only one of the two outcomes for each stage goes on to a new stage. In other words, its basically a Heckman model, where you observe an outcome equation (wages) but only for a subgroup passing the selection equation (only those who decide to work). My problem is with specifying the covariance matrix. It basically needs to satisfy 3 things: 1) must be positive definite 2) all diagonal entries are = 1 for identification 3) off-diagonal entries are thus correlations, and must be in [-1,1]. To ensure the cov matrix is positive definite (PD), people usually just parameterize a lower triangular Cholesky-type matrix, C. Then define the cov matrix=CC', which is guaranteed PD. However, in this decomposition C, the diagonals take the form sqrt(1-x). These values often wind up imaginary in my estimation, even if I brute force them back into good values (presumably since the canned optimization does its thing after my code evaluates). Secondly, I am unable to restrict the off-diagonal entries within [-1,1], so I wind up with values of magnitude >1, yet they are supposed to be correlations. Does anyone have any advice or experience with these issues? Much appreciated.


title: How to get confidence intervals for specific indirect effects in Amos?

text: Say I have an SEM model with 1 predictor variable (IV), 2 mediators (MV1, MV2) and 1 dependent variable (DV). Amos reports combined indirect effects for an IV on a DV. So the combined indirect effect would be the sum of the component indirect effects (i.e., the IV-MV1-DV indirect effect and the IV-MV2-DV) indirect effect. Amos makes it easy to get bootstrap confidence intervals on combined indirect effects. However, hypotheses often concern component indirect effects (e.g., that MV1 mediates the effect of IV on DV). **How can I get confidence intervals on component indirect effects using Amos?** I'd be interested in both asymptotic and bootstrap confidence intervals. Could I just use the values and standard errors for the individual coefficients (i.e., MV1<-IV, and DV<-MV1)?


This question has come up in another thread that I started so I thought I would get more people's opinions on it. My question is **Is the residual, e, an estimator of the error, $\epsilon$?** The reason that I ask is the following. In OLS, the variance of the residuals, $\frac{\text{RSS}}{(n - K )}$, is known as the variance of the regression (where RSS is the residual sum of squares). Similarly the square root of this variance, $\sqrt\frac{\text{RSS}}{(n - K )}$, is the standard error of the regression. Given the fact that the square root of the variance, $\frac{\text{RSS}}{(n - K )}$,is a standard error, it must mean that this variance is the variance of an estimator. We already know that it is the variance of the residuals, therefore, the residual is an estimator?? (I assume of $\epsilon$) Thoughts??


title: Multivariate Multiple Regression with Repeated Dependent Measures [MATLAB]

text: [Originally posted in Mathematics before I realized this was a better venue] I have a dataset with $p$ predictors for $i$ items (so multiple regression). For each of $s$ subjects, I have $r$ repeated observations of $v$ dependent variables (so it's a multivariate problem). I wish to find out for each individual how much variance the $p$ predictors account for across all $i$ items. The formula for calculating the parameter coefficients (aka betas) (using MATLAB notation) is `Betas=pinv(X'*X)*X'*Y` where $X$ is the $[i \times p]$ matrix of predictors, and $Y$ is the matrix of observations. Unfortunately, my $Y$ is a 3D matrix of size $[i \times r \times v]$. The matrix multiplication requires a 2D matrix; however I am not sure about which is the correct approach:   1. Concatenate the $r$ observations of $v$ measurements together so that e.g. $X=[i \times p], Y=[i \times (r \cdot v)]$   2. Concatenate the $r$ observations of $i$ items together, and replicate the predictor matrix $r$ times so that e.g., $X=[(r\cdot i) \times p], Y=[(r\cdot i) \times v]$   3. Average the $r$ observations of $v$ measurements together so that $X=[i \times p]$ and $Y=[i \times v]$ The second option multiplies the predictors which, by virtue of artificially increasing the number of explanatory variables, seemed likely to inflate the goodness-of-fit. The first option multiplies the dependent variables, which then makes it impossible to use the parameters to predict a set of $v$ values. A final option I was considering was to simply average the $r$ observations, but then you lose data. Are any of these the appropriate course of action?


title: Distribution fitting problem

text: I have a real data set `275, 13, 147, 23, 181, 30, 65, 10, 300, 173, 106, 300, 300, 212, 300, 300, 300, 2, 261, 293, 88, 247, 28, 143, 300, 23, 300, 80, 245, 266` and try to fit a Weibull distribution using Maximum Likelihood. I used `Mathematica` for analysis and `Mathematica` gives me the parameter estimates of Weibull distribution as $\alpha=1.07484, \beta=171.242,\gamma=9.46641$ where $\alpha,\beta$ and $\gamma$ are the shape, scale and location parameters, respectively. Now my two questions are:   1. If I find CDF of Weibull distribution using estimated parameters for minimum value of data set answer is `Zero` which is obvious but for maximum value the answer become `0.828829` Why it is not `One`?    2. How I find the estimated value using these parameter estimates? Note. Question 1 problem also occurs for other distributions as well; here the Weibull distribution is taken as an example. `Mathematica` code's and result are here


title: How to summarize and compare non-linear relationships?

text: I have data on the percent of organic matter in lake sediments from 0 cm (i.e., the sediment - water interface) down to 9 cm for approximately 25 lakes. In each lake 2 cores were taken from each location so I have 2 replicate measures of organic matter percentage at each sediment depth for each lake. I am interested in comparing how lakes differ in the relationship between percent organic matter and sediment depth (i.e., slope). In some lakes the relationship between percent organic matter and sediment depth appears linear but in other cases the relationship is more complex (see examples below). My initial thoughts were to fit linear relationships where appropriate either to the whole curve or to a subset of the curve if it was "mainly" linear and only compare those lakes where a significant linear relationship was found. However I am unhappy with this approach in that it requires eliminating data for no other reason than they do not fit the linear model and it ignores potentially interesting information about the relationship between percent organic matter and sediment depth. What would be a good way to summarize and compare the curves from different lakes? Thank you Example curves: In all cases the y-axis is the percent organic matter in the sediment and the x-axis is the sediment depth where 0 = the sedi ment-water interface. A nice linear example: ![A nice linear example](http://i.stack.imgur.com/TxzGA.png) 2 non-linear examples: ![nonlinear 1](http://i.stack.imgur.com/6M4N8.png) ![enonlinear 2](http://i.stack.imgur.com/ESuxB.png) An example with no obvious relationship: ![no relationship](http://i.stack.imgur.com/4Cczp.png)


title: What is the difference between (universal) kriging and spatial autoregressive models?

text: As part of a course on missing observations in social/survey statistics I am trying to explore existing methods of predicting either point pattern or polygon data. I got quite confused by all the different terms and formalisations used in **geostatistics** (more on the kriging side) on the one hand and **spatial statistics/econometrics** (more on the regression side) on the other. So far, I understood that both **spatial autoregression** and **kriging**   * seek spatial prediction,    * rely on a variance-covariance (or variogram) Matrix - in other words, they rely on spatial autocorrelation,   * and, as far as _universal_ kriging is concerned, assume an underlying trend plus random residuals. However, kriging seems to be used primarily for point patterns with continuous regions while spatial autoregression involves aggregations of phenomena into polygons. * * * So here are my questions:   * Is the underlying spatial structure the only difference between the two?   * Are there equivalent mathematical model formulations for universal kriging and spatial autoregression?   * Would area-to-area universal kriging be exactly the same thing as prediction models using spatial error autoregression?


title: Getting the log-likelihood per fold in LOO Cross Validation for HMM in Matlab

text: I am running Leave-one-out Cross validations for deciding the optimal number of hidden states in an HMM. At each iteration I am getting one model and with forward algorithm I estimate the probability of the testing data based on the model. I also need the log-likelihod (LL) of the model at each fold in order to get the sum of all folds and finally have a goodness of fit index. How can I calculate LL per fold using matlab? I am asking whether there is a specific function for the calculation of it? And last question, I get one model per fold how can I get the global model? By getting the average of probabilities?


title: Probabilistic interpretation of regression for justifying squared loss function

text: I was reading Andrew Ng's CS229 lecture notes (page 12) about justifying squared loss risk as a means of estimating regressions parameters. Andres explains that we first need to assume that the target function $y^{(i)}$ can be written as: $$ y^{(i)} = \theta^Tx^{(i)} + \epsilon^{(i)}$$ where $e^{(i)}$ is the error term that captures unmodeled effects and random noise. Further assume that this noise is distributed as $\epsilon^{(i)} \sim \mathcal{N}(0, \sigma^2)$. Thus: $$p(e^{(i)}) = \frac{1}{\sqrt{2\pi\sigma}}exp \left( \frac{-(e^{(i)})^2}{2\sigma^2} \right) $$ Thus we can see that the error term is a function of $y^{(i)}$, $x^{(i)}$ and $\theta$ as in: $$e^{(i)} = f(y^{(i)}, x^{(i)}; \theta) = y^{(i)} - \theta^Tx^{(i)}$$ thus we can substitute to the above equation for $e^{(i)}$ $$p(y^{(i)} - \theta^Tx^{(i)}) = \frac{1}{\sqrt{2\pi\sigma}}exp \left( \frac{-(y^{(i)} - \theta^Tx^{(i)})^2}{2\sigma^2} \right)$$ Now we know that: $p(e^{(i)}) = p(y^{(i)} - \theta^Tx^{(i)}) = p(f(y^{(i)}, x^{(i)}; \theta))$ Which is a function of the random variables $x^{(i)}$ and $y^{(i)}$ (and the non random variable $\theta$). Andrew then favors $x^{(i)}$ as being the conditioning variable and says: $p(e^{(i)}) = p(y^{(i)} \mid x^{(i)})$ However, I can't seem to justify why he says that using f alone. I.e. using $f(y^{(i)}, x^{(i)}; \theta)$ doesn't seem to favor any variable to be conditioned (or at least in that form). The problem I have with his conclusion is that, with only the equation: $$\frac{1}{\sqrt{2\pi\sigma}}exp \left( \frac{-(e^{(i)})^2}{2\sigma^2} \right)$$ I couldn't really see why $p(e^{(i)}) = p(x^{(i)} \mid y^{(i)})$. The Normal equation seems to be symmetrical in $x^{(i)}$ and $y^{(i)}$. Why favor $p(y^{(i)} \mid x^{(i)})$ and not $p(x^{(i)} \mid y^{(i)})$. Furthermore, if its a supervised learning situation, we would get both pairs $(x^{(i)}, y^{(i)})$, right? Its not like we get one first and then the other. Basically, I am just trying to understand why $p(y^{(i)} \mid x^{(i)})$ and not $p(x^{(i)} \mid y^{(i)})$.


title: how can I include a variable ranging from 0 to 1 in a continuous linear model

text: Imagine that I am trying to estimate house prices with the following model: $Y_p=1 + \beta_{surface} * X_{surface} + \beta_{New York}*X_{New York} + \beta_{Boston}*X_{Boston} $ With $X_{Boston}$ a dummy variable if in Boston In other words, I assume that the house prices are linearly correlated with surface area, and that there is different mean price in Boston as in New York. Now I would like to include a house that is between Boston and New York. I could set the both $X_{New York}$ and $X_{Boston}$ to 0.5. But I am a bit stuck here, as now the predictor is not a dummy anymore but a continuous variable. What should I do to include all houses spread in between New York and Boston and still know the difference in mean price between the two cities?


title: AR terms and independent variable as regressors

text: After trying several models with my data, R^2 and p values are showing my model looks like below. ACF plot tells me AR term is significant. Insights into data tells me change in 'x' would have impact(say y is A/C usage and x is temperature).                 y(t)=a+b*y(t-1)+c*{x(t)-x(t-1)}+error      Question: Is the above model statistically meaningful? If so, what is the category of this "linear model" called? I came across AR models, MA models and ARIMA/ARMA models, but havent seen a model where AR term and diff(independent variable) exists together.


title: Changing the variable you omit due to the dummy variable trap

text: I've got a quick question. The practice exam on which I have the question is **linked below** and its the **first question** As you can see, **column (3) and (4)** on the table are both the same except the dummy variable they omit is different. We've been asked to work out the coefficients and standard errors for the new regression which changes the omitted variable and I've worked out all the new coefficients but I don't know how to work out the **new standard errors** of the 4 location dummy variables. Practice Exam - The relevant q is Question 1


title: How can we choose which model generalizes better?

text: I got asked this question today: > How could we reasonably safely conclude that, for a given dataset, in > practice, one classifier generalizes better than another? I could not give a concise answer, and thought maybe you guys would like to chip in. I understand obviously that there is still much in terms of research in this area that needs to be done, however it I think its still worth an attempted answer.


title: Compute Shannon entropy between every row of a large, sparse matrix

text: I have a sparse, binary matrix of user (rows) and items (columns). Each element of this matrix is either 0 or 1:               require(Matrix)     set.seed(42)     rows <- 10     cols <- 2     non_zero <- 20     x <-  sparseMatrix(i = sample(1:rows, non_zero, replace=TRUE),                        j = sample(1:cols, non_zero, replace=TRUE),                        x=1)     colnames(x) <- letters[1:ncol(x)]     x <- sign(x)     print(x)     10 x 2 sparse Matrix of class "dgCMatrix"           a b      [1,] . .      [2,] 1 1      [3,] . 1      [4,] . .      [5,] 1 1      [6,] . 1      [7,] 1 .      [8,] 1 1      [9,] . 1     [10,] 1 1      (In reality I have about 1e7 row, 1e6 columns and about 1e8 elements in this matrix) I would like to compute Shannon Entropy between every column of this Matrix (as described here). I would like to do this using sparse matrix operations, as using dense matrix operations will be extremely inefficient on a matrix of this size. It's pretty easy to calculate the various frequencies tables need, e.g.: A and B counts:               a_and_b <- crossprod(x)     print(a_and_b)     2 x 2 sparse Matrix of class "dsCMatrix"       a b     a 5 4     b 4 7      B, not A counts and A, not B counts:               b_not_a <- crossprod(sign(a_and_b), Diagonal(x=diag(a_and_b)))     a_not_b <- t(b_not_a)          a_not_b <- a_not_b - a_and_b     b_not_a <- b_not_a - a_and_b          print(a_not_b)     2 x 2 sparse Matrix of class "dgCMatrix"       a b     a 0 1     b 3 0          print(b_not_a)     2 x 2 sparse Matrix of class "dgCMatrix"       a b     a 0 3     b 1 0      And finally, not A not B counts:               not_a_b <- crossprod(       sign(a_and_b),        Diagonal(n=ncol(a_and_b), x=nrow(x))       ) - a_not_b - b_not_a - a_and_b     print(not_a_b)     2 x 2 sparse Matrix of class "dgCMatrix"       a b     a 5 2     b 2 3      I can compute these 4 contingency matrices (A and B, A not B, B not A, not A and not B) extremely fast, even on huge matrices. How can I use them to calculate the Shannon entropy between each item (column) of my original matrix `x`? The final output should look something like this:                    a    b     a 0.43 0.89     b 0.96 0.64


title: Confidence regions for non-normal trivariate distributions

text: We are analyzing haplotype diversity within and between six populations or demes. We came up with a measure that compares the similarity of every possible pair of haplotypes within a deme or between two demes. Hence, for each comparison between two demes we have three measures: _Hi_ (similarity within deme _i_ ), _Hj_ (the same for deme _j_ ) and _Hij_ (similarity between deme _i_ and _j_ ). Observed results are plotted in figure 1 (there are 15 possible comparisons between six demes). Note that _Hi_ and _Hj_ are independent from each other, but _Hij_ depends on both! For high values of _Hi_ and _Hj_ (high similarity of haplotypes) _Hij_ can be both high or low (if the haplotypes differ a lot between demes but not within them). For low values of _Hi_ and _Hj_ , _Hij_ can only be low. ![1- Observed results, 2- 50 generations, 3- 150 generations](http://i.stack.imgur.com/mQL0S.png) We developed a software to simulate spatially the generation of haplotype diversity under gradients of selection in different scenarios with variable demographic parameters. For each combination of parameters we run 1000 simulations and, for each simulation, we computed the three measures above ( _Hi_ , _Hj_ and _Hij_ ) for several subsets of demes (not necessarily six demes). These results provide a cloud of three-dimensional points, which are represented in figures 2 and 3 (for different generations). Our main goal is to know if the observed results are “plausible” or “likely” under any of the demographic scenarios we used. We first consider a centroid- based approach (comparing the distances of the observed points to the centroid of the cloud with the distances of every simulated point to the centroid), but the asymmetry of the clouds ruled it out. Right now, we are looking for some form of estimating confidence regions for three dimensional distributions, and came across the paper A New Distribution-Free Approach to Constructing the Confidence Region for Multiple Parameters. We would appreciate some advice on this issue. Is this the right analysis or are there any better approaches? Sincerely yours


title: Decision Trees - Test & Training Data Difference

text: In decision tree modelling what is the difference between Training & Test Data ? i am having some difficulty in locating the meaning of these two words and i hope you guys can help me out. Could you maybe also give me some simple examples along with your definitions , just so i can practically see and understand the difference between Training & Test data. I tried to post this on stack overflow but i was recommended that i post here as you guys might be more knowledgeable in the field of Decision Trees. Hope you can help ,and much thanks.


title: is Shapiro Wilk test insensitive on the tails?

text: Is Shapiro Wilk test insensitive on the tails of sample distribution? I've read such statement on a paper but I cannot figure out why from its test statistic W. Could someone help me understand this?


title: Comparing tournament results to chance

text: Basically what I want to achieve is to take loads of data from a certain game's tournaments and compare it to chance, i.e., I want to see whether a "coin flipping tournament" would have similar results. Or in other words, whether any given player has ~50% chance to win against any other given player, and at the start all of them have equal chance to win a tournament (due to equal skill level, too much randomness in the game, or other reasons, doesn't matter). We have players that became champions once, more than once, or never. I assume that would happen in a coin flipping tournament as well. My understanding of statistics is rudimentary at best, so first question is, is this even possible to check whether we have a real skill differences and can predict some players to be more likely to win in the future? Or do we only have pure 50/50 chance? If yes, what would be the method? Or perhaps, what do I need to read up on before starting it?


title: Clarification in definition of Holt-Winters algorithm

text: The algorithm mentioned on Wikipedia has a line $$s_t = \alpha \frac{x_t}{c_{t-L}} + (1-\alpha) (s_{t-1} + b_{t-1})$$. For $1 < t < L$ how should we interpret $c_{t-L}$? On that interval, should we replace it with $c_t$?


title: Comparing odds ratios of continuous and discrete variables

text: I need to compare the ability of two methods to predict an event with a binary response. Each method produces a score, where the higher score indicates 1 and a lower score indicates 0. I am looking to calculate the odds ratio (in univariable and multivariable analysis) to compare these methods. Here is the break down of method data: Method 1: Probability between 0 and 1. Method 2: Score between 0 and 12 In univariable analysis, the OR for method 1 will be VERY large, but for method 2 would be quite small. I think that because the two methods have different scales the odds ratios become so vastly different. The unit increase for method one should be 0.1, whereas for method two is 1. Possible solutions are to binarize both methods, by picking a cutoff, then everything is on the 0 and 1 scale and the OR are comparable, however you loose information when you do this. Also, it can be trickey to pick an objective cutoff. Are there any suggestions for comparing the odds ratios of continous and discrete variables?


title: Statistical model for my problem

text: I'm running a benchmark to find the efficiency of a my computer. There are $p$ control variables say, $x_1,x_2,...,x_p$ and one output variable $Y$. For example, every time I run an experiment I change the control variables and see the output (computer efficiency). So basically, for this problem I used regression to model this problem $\\{Y_i,\, x_{i1}, \ldots, x_{ip}\\}_{i=1}^n$. I used regression also because of the nature of the experiment i.e. If I'm conducting the same experiment in another computer I want to "fit" that computer w.r.t. to the values I have from my previous computer - which is fitting a line problem. Now, in the new setting I have the same control variables but two outcomes $Y_{i1},Y_{i2}$ which are throughput and bandwidth. I'm trying to vary the control variables so as to maximize profit i.e. get $high$ throughput and consume $less$ bandwidth. How should I go about modeling this experiment? I've knowledge of only basic linear algebra which is why I modelled first one using regression "easily". I would appreciate any help to put me thinking in the right direction.


title: Multivariate regression or PCA to reduce response variables?

text: I hope the title is self-explanatory, but essentially I want to know which method is better: does it make sense to use a PCA to reduce a number of response Y variables and then conduct a univariate multiple regression OR conduct a multivariate multiple regression using all response variables as they are? My variables are all continuous. I will be conducted the analysis in R. Thanks in advance!


title: Using GLM where one predictor is linear and one is logarithmic

text: I am fitting a GLM to binomial data with two independent variables, `a` and `b`. I have tried to fit a binomial GLM with the form `y~a+b+a*b`, so that I have the contributions of the two variables as well as an interaction term. This works OK, but in the data the effect of `b` is very gradual in comparison to `a`. The best fit I've tried by far is simply `y~a+log(b+1)+a*log(b+1)` (the +1 is because b=0 for some datapoints) with a binomial link function. Is it possible to use this as a GLM? It's just that I've not seen this done when I've searched online so I don't know if there's some good reason that people only seem to use polynomials. There is a question that seems related in a previous post Interpret interactions and logarithms in linear regression, but this was using a log-log function, and used `lm()` instead of `glm()`; so I'm not at all confident of the relation.


title: Estimating number of unique people

text: Assume from a prior experiment with have a known truth table of misclassification of logins on an individual basis. $$ \begin{array}{c|lcr} \text{Truth}/\text{Observed} & \text{Al (M)} & \text{Ben (M)} & \text{Claire (F)} & \cdots\\\ \hline \text{Al (M)} & 10 & 2 & 4 & \cdots\\\ \text{Ben (M)} & 5 & 20 & 6 & \cdots\\\ \text{Claire (F)} & 3 & 1 & 30 & \cdots\\\ \vdots & \vdots & \vdots & \vdots & \ddots \\\ \end{array} $$ (Example shown here for 3x3 for illustration purposes, actual table is significantly larger. Demographics can be Age+Gender as well.) Is there a way to use the calibration table of truth to help infer the estimated actual number of unique individuals from another larger study where only observed unique individuals by demographic and observed number of logins by demographic are known? Example: $$ \begin{array}{c|lcr} \text{Observed} & \text{Unique Indv.} & \text{Total Logins} \\\ \hline \text{M} & 257 & 5588 \\\ \text{F} & 603 & 4412 \\\ \vdots & \vdots & \vdots \\\ \end{array} $$ $$ \begin{array}{c|lcr} \text{Estimated Truth} & \text{Unique Indv.} \\\ \hline \text{M} & ? \\\ \text{F} & ? \\\ \vdots & \vdots \\\ \end{array} $$


title: contingency table and noise

text: I have a contingency table and I want to use it to evaluate my clustering method's results (e.g. RandIndex). I have also noises in my data which are labeled as NOISE and they obviously should not be clustered. My question is: Should I add one row in the contingency table? or ignore the noise. Moreover, should I consider the number of noise points for the summation of the rows and cols. For example I have the following contingency table: ![enter image description here](http://i.stack.imgur.com/dSwud.png) My synthetic data has a class label (NoiseA) which I know that it is noise. My clustering method also finds some data points as noise which I called it as: NoiseC. I want to know that should I consider the noise data points either calculated or actual data for calculating the TP/TN/FP/FN in order to calculate RandINdex and other evaluations criteria?


title: Is there a way to compute daily percentiles (median and 95th) using 24 hourly percentiles

text: Basically I have max, count, median and 95th percentiles for every hour for a data stream already pre-computed from an hourly workflow. Can I use these 24 snapshots to get an approximate median and 95th percentile for the day? Is there a way / formula or technique to do this? EDIT: Sorry i should have given the information on distribution. This is a data related web page performance timing. I am assuming it as approximately normal.


title: Creating a cluster analysis on multiple variables

text: I am working on creating a cluster analysis for some very basic data in r for Windows [Version 6.1.76]. The groups themselves are countries and then I have 2 column with continuous numerical variables. I have applied a Ward Hierachical Method to the data               # Applying Ward Hierarchical Clustering     d = dist(conversion_set, method="euclidean")     fit = hclust(d, method="ward")      But I don't feel this represents what I am really trying to get to as it is just taking into account the first variable and disregarding the second. Is there a way to include both variables into the clustering calculations? My data looks similar to this Country - Var 1 - Var 2 US - 10 - 20 Canada - 5 - 30 ....


title: R Sales Modeling with Non-Replenishing Inventory

text: I'm working on modeling sales for tickets for a specific event over time. The issue that I'm having is that the inventory does not replenish, so I can't have my model predicting over what I actually have available. So for example, if I have 100 tickets available 30 days out and 20 sell, the remaining pool I have to sell from is only 80 tickets, so I can't predict to sell 90 tickets the next day because the original 100 tickets aren't all available any more. Is there a way to take this into consideration when building a model?


title: Dummy coding for contrasts: 0,1 vs. 1,-1

text: I'm seeking your help in understanding the difference between two different contrasts for dichotomous variables. On this page: http://www.psychstat.missouristate.edu/multibook/mlt08.htm under "Dichotomous Predictor Variables", there are two ways to code dichotomous predictors: using the contrast 0,1 or the contrast 1,-1. I kind of understand the distinction here (0,1 is dummy coding and 1,-1 adds to one group and subtracts from the other) but don't understand which to use in my regression. For example, if I have two dichotomous predictors, gender (m/f) and athlete (y/n), I could use contrasts 0,1 on both or 1,-1 on both. What would be the interpretation of a main effect or an interaction effect when using the two different contrasts? Does it depend on whether my cells are of different sizes?


title: Am I interpreting this coefficient correctly?

text: I'm trying to get caught up in some notes, and the topic is pooled, cross- sectional data. If we have the regression: $$\log{wage_i}=\beta_0+\beta_1\text{married}_i+\delta_0\text{yr10}_i+\delta_1\text{married}_i\times \text{yr10}_i+\varepsilon_i$$ Where $\text{married}$ is a dummy variable which $=1$ if a person is married and $0$ if not, and $\text{yr10}$ is a dummy variable indicating whether the observation was recorded in the year 2010, as opposed to being recorded in the year 2000. There's a blank spot in the notes that I didn't get, and it asks: How do we interpret $\delta_1$? I think it's correct to say that $\delta_1$ measures how the return on whether or not someones married has changed in the past $10$ years. Is that correct? Or am I wrong? I read a little bit in my textbook and I think that's what it's saying.


title: Interpretation of dummy-variable when the dependent variable is in percent and log

text: I have the following estimations of a fixed effects model: Open=β0+ β1Pol_r+ β2Infr+ β3Nat_re+ δ4SACU+ɤi+εit Open_ln=β0+ β1Pol_r+ β2Infr+ β3Nat_re+ δ4SACU+ɤi+εit Where Open is the openness of a country measured as ((imports+exports)/GDP)*100, Open_ln is the natural logarithm of Open. SACU is a dummy-variable taking the value 1 for countries becoming members of a free trade area in a particular year. How do I interpret the dummy-variable when Open is in level-form? How do I interpret the dummy variabel when Open is in log-form? I get the coefficient SACU 45.37 in the first case, and the coefficient SACU 0.566 in the second case.


title: Means of determining monotonicity of missing pattern other than eyeballing? And how monotone is monotone for the purpose of multiple imputation?

text: I have a large dataset with a large number of variables. Missingness as high as ~25% in some variables, many vars with no missing. Judgment of the monotonicity of the missing pattern is important for the selection of a multiple imputation method. I have seen some sources talk about "monotone or near-monotone" patterns. But how monotone is monotone, and how do you know, other than by looking? Every text I have seen recommends eyeballing the patterns. But is there not some means of statistically summarizing monotonicity (e.g., sum of belongingness the main missing pattern controlling for belongingness to other patterns)? I have not seen anything like this in any text, but maybe I'm just missing it. If I organize the data matrix to maximize the number of missing measurements clustered in a monotone pattern, I have about 35% of measurements in the pattern. (Other measurements display something of a "fan" pattern similar to that seen here. How would I decide whether this were monotonic "enough" to assume it for the purpose of MI? I know that I could also do MI in two steps in order to isolated the monotonic pattern, then run MI again with a monotone method. As fun as this sounds, I'm trying to determine whether it's reasonable to do it in one fell swoop.


title: How to split the dataset for cross validation, learning curve, and final evaluation?

text: **What is an appropriate strategy for splitting the dataset?** _I ask for feedback on the following approach_ (not on the individual parameters like `test_size` or `n_iter`, but if I used `X`, `y`, `X_train`, `y_train`, `X_test`, and `y_test` appropriately and if the sequence makes sense): (extending this example from the scikit-learn documentation) ### 1\. Load the dataset               from sklearn.datasets import load_digits     digits = load_digits()     X, y = digits.data, digits.target      ### 2\. Split into training and test set (e.g., 80/20)               from sklearn.cross_validation import train_test_split     X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)      ### 3\. Choose estimator               from sklearn.svm import SVC     estimator = SVC(kernel='linear')      ### 4\. Choose cross-validation iterator               from sklearn.cross_validation import ShuffleSplit     cv = ShuffleSplit(X_train.shape[0], n_iter=10, test_size=0.2, random_state=0)      ### 5\. Tune the hyperparameters applying the cross-validation iterator on the **training set**               from sklearn.grid_search import GridSearchCV     import numpy as np     gammas = np.logspace(-6, -1, 10)     classifier = GridSearchCV(estimator=estimator, cv=cv, param_grid=dict(gamma=gammas))     classifier.fit(X_train, y_train)      ### 6\. Debug algorithm with learning curve `X_train` is randomly split into a training and a test set 10 times (`n_iter=10`). Each point on the training-score curve is the average of 10 scores where the model was trained and evaluated on the first _i_ training examples. Each point on the cross-validation score curve is the average of 10 scores where the model was trained on the first _i_ training examples and evaluated on all examples of the test set.               from sklearn.learning_curve import learning_curve     title = 'Learning Curves (SVM, linear kernel, $\gamma=%.6f$)' %classifier.best_estimator_.gamma     estimator = SVC(kernel='linear', gamma=classifier.best_estimator_.gamma)     plot_learning_curve(estimator, title, X_train, y_train, cv=cv)     plt.show()      ![Learning curve](http://i.stack.imgur.com/XJ4nN.png) plot_learning_curve() can be found in the current dev version of scikit-learn (0.15-git). ### 7\. Final evaluation on the test set               classifier.score(X_test, y_test)      ### 7a. Test over-fitting in model selection with nested cross-validation (using the whole dataset)               from sklearn.cross_validation import cross_val_score     cross_val_score(classifier, X, y)      **Additional question:** _Does it make sense to replace step 7 by nested cross-validation? Or should nested cv be seen as complementary to step 7_ (the code seems to work with k-fold cross validation in scikit-learn, but not with shuffle & split. So `cv` needs to be changed above to make the code work) ### 8\. Train final model on whole dataset               classifier.fit(X, y)      EDIT: I now agree with cbeleites that step 7a doesn't make much sense in this sequence. So I wouldn't adopt that.


title: Simultaneous Probability

text: There are three drivers and two roads, X and Y. The probability of choosing road X is p for all drivers, and the probability of choosing road Y is (1-p) for all drivers. To calculate the probability that all drivers choose road X simultaneous is p^3. However, what would 3p imply? What does it mean intuitively?


title: Is there an advantage to squaring dissimilarities when using Ward clustering?

text: Is there a reason to prefer squaring or not squaring the dissimilarities when clustering with Ward's method? The question is motivated by the following statement in the documentation for R's `hclust()` function: > Two different algorithms are found in the literature for Ward clustering. > The one used by option "`ward.D`" (equivalent to the only Ward option > "`ward`" in R versions <= 3.0.3) _does not_ implement Ward's (1963) > clustering criterion, whereas option "`ward.D2`" implements that criterion > (Murtagh and Legendre 2013). With the latter, the dissimilarities are > _squared_ before cluster updating. Does squaring improve the algorithm?


title: Post hoc test in a 2x3 mixed design ANOVA using SPSS?

text: I have two groups of 10 participants who were assessed three times during an experiment. To test for differences between groups and across the three assessments, I ran a 2x3 mixed design ANOVA with `group` (control, experimental), `time` (first, second, three), and `group x time`. Both `time` and `group` resulted significant, besides there was a significant interaction `group x time`. I don't know very well how to proceed to further check for the differences between the three times of assessments, also respect to group membership. In fact, at the beginning I only specified in the options of the ANOVA to compare all the main effects, using the Bonferroni's correction. However, then I realized that this way they were compared the differences in time of the total sample, without group distinction, am I right? Therefore, I searched a lot on the internet to find a possible solution, but with scarce results. I only found 2 cases similar to mine, but their solutions are opposite!   1. In an article, after the mixed design, the authors ran 2 repeated measures ANOVA as a post-hoc, one for each group of subjects. This way, the two groups are analysed separately without any correction, am I right?   2. In a guide on the internet, they say to add manually in the SPSS syntax `COMPARE(time) ADJ(BONFERRONI)`, just after `/EMMEANS=TABLES(newgroup*time)`, while running the mixed ANOVA. This way, the three times are compared separately for each group, with Bonferroni correction, am I right? What do you think? Which would be the correct way to proceed?


title: Need help with converting or analyzing uneven time intervals

text: I have next to no idea what I'm doing. I have unevenly sampled poll data over an 11 year time period and a database of events. I wish to show whether or not increases in the number of these events correlate to increases in the percentage of affirmative poll responses (the poll's responses are ranked "very," "not very," "not too" and "not at all," I'm lumping the first two together as affirmative and the latter as negative). Dates of the polls: http://puu.sh/1rj6l How would one go about converting this unevenly collected poll data to a normal time series? As you can see, there are some points where the polls are very close together and some where they are a year or so apart. Is it viable to convert it like this? The final interval doesn't much matter - what I'm doing is summing the number of events since the last interval and seeing how changes in that number correlate to changes in people's responses. I have Excel and SPSS at my disposal.


title: Fraction of variance unexplained and R-squared in linear and non-linear regression

text: I have a non-linear model of the following form: $y = a*x^b$ I can fit it using logarithms and a linear model or directly with a non-linear model. First approach, logarithms and linear model:               lmfit <- lm(log(y)~log(x))      Second approach, non-linear model:               nlsfit <- nls(y~a*x^b, start=list(a=200, b=1.6))      In the first case I can simply get the $R^2$ value from the linear model or calculate it myself by:               rsquared <- var(fitted(lmfit)) / var(log(y))      In the second case there is no $R^2$ value generated, but I can obtain one $pseudoR^2$ value myself by:               pseudorsquared <- var(fitted(nlsfit)) / var(y)      In a linear model I can calculate the fraction of variance unexplained by simply doing $1-R^2$. I have read that this is not applicable to non-linear regressions. I would like to know if there is an equivalent version of this measure, so that I can compare both regressions and use the best one. As an extra information, I would like to add that this is a regression of physical variables, and that the non-linear approach is providing more close- to-literature results for the coeficients, whereas the linear approach gives better statistical performance ($R^2$, bias, etc.).


title: Comparing treatment means with lots of zeros

text: I am trying to compare the means of two treatments on a continuous variable with a lot of zeros in it. I've tried a log(n+1) transformation but that did not get me to a normal distribution. Any suggestions?


title: Plotting binomial proportions on a box & whiskers plot using R

text: Apologies in advance - I'm a bit of a newbie with R but I've been Google-ing this for weeks and I haven't found a straight answer. **The problem** Basically, I'm looking to plot proportions with their confidence intervals on a box and whiskers plot using the R environment. i.e. Plot the results of binom.test(100,1000) and binom.test(125,1000) What is the best package (I don't have a preference either way) and way to go about doing this?


title: Understanding expectation maximization for simple 2 linear mixture case

text: I would appreciate some help getting some EM stuff straight. So, say I generate data in R as follows:               N       <- 100     epsilon <- rnorm(N)     X       <-  10*runif(N)     beta.0  <- 10     beta.1  <-  3     sigma   <- 2     Y       <-  beta.0 + beta.1 * X + sigma * epsilon     epsilon2 <- rnorm(N)     X2 <- 10*runif(N)     Y2 <-  3 - X2 + 0.25 * epsilon2     Y.mix <- c(Y, Y2)     X.mix <- c(X, X2)      Now, in expectation maximization, in the first step, I have some prior probability, say 0.5, of the data being from either one or the other distribution. So, using EM I know I can estimate the mean and variance of the two mixtures. From looking at a density plot, it seems like the means are at about -2 and 30 for the data I simulated. But, at what stage in EM do I back out the betas? I want to recover the slope, intercept, and sd deviation parameters for the 2 regression-type equations. Thanks for an clarification.


title: How should I deal with features with time related values in Bayesian network?

text: I would like to apply Bayesian network on some data. However, some of the variables are related to time. E.g. Number of time he/she visit library. As the value can be defined as Total number of this person visits, Average weekly number of he/she visit, Daily number of visits, or even Number of visits in the last 3 days etc... I do not want to include all these different variations in my network Then how should I treat feature like this properly? Many thanks for your help


title: Expected number of runs of length n in coin toss

text: I'm writing a program for an assignment for my Probability and Statistics course. This program performs several millions of simulated coin tosses, and then analyses various properties of the generated sequence. One of the things it does is count how many runs does the generated sequence contain, and their lengths. One result for 10.000.000 tosses is:               Number of runs of length 1 is: 2500986     Number of runs of length 2 is: 1246647     Number of runs of length 3 is: 625656     Number of runs of length 4 is: 311689     Number of runs of length 5 is: 156673     Number of runs of length 6 is: 78464     Number of runs of length 7 is: 39253     Number of runs of length 8 is: 19547     Number of runs of length 9 is: 9866     Number of runs of length 10 or more is: 9818      Now, my problem is following: I need to compare the result acquired by simulation to theoretical expectation of sorts. Basically, I need that same result as above, calculated on paper. A pattern in the results obviously exists, but I just cant figure out a formula. Anyone, help please?


title: How to compute Confidence Interval associated to a Binomial proportion's increase?

text: Not sure if you can help me but here's the problem. I have two binomial proportions A and B (95% CI) - $A = 2\% \pm 0.2\%,\quad B = 3\% \pm 0.2\%.$ In other words, B's proportion is 50% higher than A's. Business people tend to prefer speaking about this difference in terms of a percentage increase. In other words, B's proportion is 50% higher than A's. How do I go about calculating a Confidence Interval for the 50% increase of B's proportion over A. I've done some research into Fieller’s formula which is said to provide exact percent effect. However, besides not understanding Fieller's theorem, I don't know how I would go about calculating this in R. Can anyone help? **Edit:** E.g., if the CI for those intervals is 95%, the CI for the 50% increase = 98.52%.


title: Geostatistical analysis using spatial.exp in WinBugs

text: I have a logistic regression model in GeoBugs to estimate predictors of prevalence of a disease. Can anyone tell me if there is a simple way to determine the lower and upper bounds of phi for spatial.exp. This is part of the model:               Infected[i] ~ dbin(p[i], tested[i])     logit(p[i]) <- alpha + beta1*covar1[i] + beta2*covar2[i]  + u[i]      mu[i] <- 0     u[1:N] ~ spatial.exp(mu[], x[], y[], tau, phi, 1)      This how I have specified the priors:               phi ~ dunif(0.25, 20)      At the moment my selection has been arbitrary. The GeoBugs manual only goes as far as saying that the bounds of the distribution are determined by the data. What I would like to know is there a dummies version of the best way to select the values.


title: Elementary Probability Questions

text: Although I have read the related book chapters, I still have some confusion, so I am asking the following elementary questions. 1) What are the following equal to? (Some of them might be meaningless)   * $\sum_{x_1} p(x_1 | x_2) = ?$ (I think this is equal to $1$ for a setting of $x_2$, e.g. $\sum_{x_1} p(x_1 | x_2 = 1) = 1$)   * $\sum_{x_2} p(x_1 | x_2) = ?$   * $\sum_{x_1} p(x_1 | x_2, x_3) = ?$ (I think this is equal to $1$ for a setting of $x_2$ and $x_3$)   * $\sum_{x_2} p(x_1 | x_2, x_3) = ?$   * $p(x_1|x_3)p(x_2|x_3)p(x_3)$ = ? 2) Assuming each of $x_1$, $x_2$ and $x_3$ can take three different values, $1$, $2$ and $3$, how do we calculate the number of free parameters for the following?   * $p(x_1)$ (I think this is $3-1=2$)   * $p(x_1, x_2, x_3)$ (I think this is $3^3=27$)   * $p(x_1|x_2)$ (I think this is $(3-1)*3=6$)   * $p(x_1|x_2, x_3)$ (I think this is $(3-1)*3*3=18)$   * $p(x_1, x_2 | x_3)$ (I think it is $(3*3-1)*3$ 3) Are the following statements true?   * $p(x_1)$ is a probability distribution.   * $p(x_1, x_2)$ is a probability distribution.   * $p(x_1 | x_2)$ is a collection of probability distributions.   * $p(x_1 | x_2 = 1)$ is a probability distribution.


title: SVM decision boundary conditions : derivation problem

text: I was trying to understand the derivation of SVM decision boundary. Suppose my decision boundary is `y-x-1=0`. Now in the book it was written that `(y-x-1) = +/-(1)`, are the equation of the maximum margins. But my question is why are we choosing a distance of 1 from the decision boundary ? How do we even know that we would be able to hit the support vectors just at the distance of +/-(1) from the decision boundary ? Moreover how can we say, that if we are able to hit a support vector at a distance of +1 from decision boundary then it will be there exactly at -1 (and not at -3, -5, -10... etc for that matter) ? i.e. how are we so sure that at exactly +/-(1) "on both sides" from decision boundary we will be able to hit support vectors.


title: Can a narrow confidence interval around a non-significant effect provide evidence for the null?

text: It is obviously fallacious to assume that failure to reject the null implies that the null is true. But in a case where the null is not rejected and the corresponding confidence interval (CI) is narrow and centered around 0, does this not provide evidence **for** the null? I am of two minds: Yes, in practice this would provide evidence that the effect is more or less 0. However, in a strict hypothesis-testing framework, it seems that null effects are simply unusable for inference, as are their corresponding CIs. So what is the meaning of a CI when its point estimate is non-significant? Is it also unusable for inference or can it be used as in the preceding example to quantify evidence for the null? Answers with scholarly references are encouraged.


title: How to generate random data that conforms to a given mean and upper / lower endpoints?

text: I'm in need of generating a set of random numbers (between 100 and 300 would suffice), which conform to a given mean, and also fall within certain upper and lower endpoints. Using the statistics software that I have (Minitab), I can find options to   1. generate random data with a given mean but no endpoints, and    2. generate random data with endpoints and a given mode (not mean).  I think the problem I'm having is that it doesn't seem to be a normal distribution; for instance, one of the data sets should give a mean of 2.90, but range between 0 and 30.66. It is for a paper I'm writing, wherein the concern is with being able to run statistical analyses on the data in question; in other words, once I generate the data set, I have to do ANOVA and other tests on it. (I'm fine with doing all that; I just can't figure out how to generate the data.)


title: Interpretation of three-way interaction in an output of lmer

text: I am using lme4 to run a three-way interaction model. I have three independent variables: animal (rat, lion, dog), color (red, green, blue) and sex (male, female). The baselines are as follows: Animal: dog Color: blue Sex: female My ultimate aim is to find out differences between rat VS dog, and lion VS dog. Then I interpret the output in this way:               1) rat = rat:blue:female     2) lion = lion:blue:female     3) red = dog:red:female     4) green = dog:green:female     5) male = dog:blue:male     6) rat:red = rat:red:female     7) lion:red = lion:red:female     8) rat:green = rat:green:female     9) lion:green = lion:green:female     10) rat:male = rat:blue:male     11) lion:male = lion:blue:male     12) rat:red:male = rat:red:male     13) lion:red:male = lion:red:male     14) rat:green:male = rat:green:male     15) rat:green:male = rat:green:male      So according to this link: How to interpret 2-way and 3-way interaction in lmer? It seems like all parameters will have to involve these three variables as I wrote them after "=" with the appearance of baselines if they are not in the parameter. Then my question is: as my goal is to find out differences between rat VS dog, and lion VS dog, I am not sure if I can directly compare rat VS dog in the first one in female [(1) rat = rat:blue:female] with rat VS dog in the tenth [(10) rat:male = rat:blue:male] as from my understanding, these two parameters are different from each other only in sex of participants. What do you think?


title: Can RBMs be used for feature selection / reduction?

text: I have a data set that's ~ 150R X 2000C and was curious if an RBM is appropriate in situation with this type of imbalance. It's a microarray and I'm looking at a 0/1 classification problem. I'd be curious if there are any coded examples of using this especially in R/Python/Java.


title: Multi-peak Gaussian fit in R

text: I have a large data set composed of several "independent" data frames like this one                    Tiempo UT1 UT2 UT3 UT4 UT5 UT6 UT7 UT8 UT9     3082  616.2   0   4   0   1   1   0   0   0   0     3083  616.4   0   1   0   0   1   0   0   0   0     3084  616.6   0   0   0   0   1   0   0   0   0     3085  616.8   0   1   0   1   2   0   0   0   0     3086  617.0   0   0   0   1   0   0   0   0   0     3087  617.2   0   0   0   0   4   0   0   0   0     3088  617.4   0   1   0   1   0   0   0   0   0     3089  617.6   0   1   0   1   1   0   0   0   0     3090  617.8   0   0   0   1   1   0   0   0   0     3091  618.0   0   1   0   1   0   0   0   0   0     3092  618.2   0   0   0   2   1   0   0   0   0     3093  618.4   1   1   0   0   0   0   0   0   0     3094  618.6   0   0   0  16   4   0   0   6   5     3095  618.8   0  11   0  23   1   0   0   6   4     3096  619.0   1  16   0  23   1   0   0   0   1     3097  619.2   0   6   0   9   0   0   0   1   0     3098  619.4   0   0   1   4   1   0   0   0   1     3099  619.6   0  18   9  16   0   1   1  15   3     3100  619.8   0  21  23  25   0  21   8  22  16     3101  620.0   0  28   3  17   1  10   3   8   2     3102  620.2   0   3   2   7   0  29   0   5   5     3103  620.4   0   2   1   4   0  23   1   2   3     3104  620.6   0   0   0   3   0  19   0   1   1     3105  620.8   0   0   0   3   0  13   0   0   2     3106  621.0   0   0   0   3   1  11   0   1   1     3107  621.2   0   0   0   1   0   8   1   1   1     3108  621.4   0   0   0   2   0  10   0   2   2     3109  621.6   0   0   0   2   0   5   0   3   2     3110  621.8   0   0   0   2   0   7   1   2   1          UT10 UT11 UT12 UT13 UT14     3082    5    0    0    0    3     3083    1    0    0    0    2     3084    0    0    0    0    1     3085    0    0    0    0    1     3086    0    0    0    0    2     3087    5    0    0    0    2     3088   13    0    0    0    1     3089    2    0    0    0    1     3090    0    0    0    0    0     3091    0    0    0    0    2     3092    1    0    0    0    1     3093    1    1    1    0    0     3094    4    2    0    0    0     3095    6    7    1    1    1     3096   12    1    0    0    3     3097    4    2    0    1    5     3098   10    3    0    0    3     3099   10    1   12    0    1     3100   15   13   21   15    0     3101   11    6    6   14    2     3102    1    3    4    2    0     3103    8    2    2    1    0     3104    1    2    3    1    0     3105    2    4    1    2    0     3106    2    3    1    2    0     3107    0    0    1    2    0     3108    1    0    1    3    0     3109    0    2    1    2    0     3110    0    2    1    4    0      Here's an example of the plot I use to analyse my data (`UT4` vs. `Tiempo`): ![](http://i.stack.imgur.com/GG7Pl.png) My goal is to fit a multi-peak Gaussian of every column `UT`$_i$ in order to get the parameters for a generic UT and use it for a further statistical analysis. I've found some ideas here using `ksmooth` fitting multiple peaks to a dataset and extracting individual peak information in R, but the result I got was a unimodal fit of my data. I hope I've made myself clear with the question. Thank you in advance.


title: Average partial effects

text: I need to explain what average partial effects (APEs) are to a very general non-statistical audience (i.e. the APEs derived from a probit model). I have tried to define APEs using layman's terms but I find it hard, can somebody please help me with this?


title: Is it possible to estimate non-nested levels in a multilevel model in R?

text: I have a research data with 3 levels of agregation: level1 is the individual, level2 is house area and level3 is school. This means, the levels are not nested in the usual way (level1 within level2 within level3). In my case, level 1 lays within level 2 and also within level 3, but level 2 and level 3 are not within each other. Is is still possible to estimate a hlm regression with such configuration with lme funcion in R?


title: Median survival time for stratified Cox model

text: I have a dataset of 30,000 patients who are matched by surgical procedure. My main exposure of interest is if the patient got regional anesthesia during surgery or did not. I have 27 different surgical procedure types, so instead of including 26 dummy variables in my Cox regression model, I stratified on them. In addition to the main exposure of interest in the model (regional anesthesia/no regional anesthesia), I have 3 other covariates I'm controlling for...Is there a simple way to calculate median survival for the patients who got regional anesthesia and for those that did not? Or would I have to calculate median survival for these two groups for all 27 strata since they all have different baseline hazards?


title: Transformation of latitudes to include them a linear model

text: I have in a dataset a variable `latitude` that contains values such as `39° 37' 39" N` and `39deg 24min`. As far as I know, latitude is measured like above internationally. My question is: To process the above data to include them in a linear model for example, which is the proper way to express or format `latitude` values, keeping their scientific and logical aspect?


title: Binomial distribution where probability of success is dependent on another binomial distribution

text: How does one model the Binomial distribution where the probability of success is the result of another Binomial distribution. For example, say I make 10 coin tosses many times and record the number of heads ( _H_ ). Then for each set ( _i_ ) of 10 coin tosses I put _Hi_ black marbles, and 10- _Hi_ white marbles in a jar and make 50 draws with replacement. How would I model the distribution of the black marble draws taking into account their dependence on the previous Binomial distribution that generated their probability of success.


title: Distribution of correlation coefficients for uniform random variables

text: Let $n>1$, let $X$ be uniformly distributed on $[-\frac12,\frac12]$, and consider the sequence $X_1,\ldots,X_{n+1}$ of independent copies of $X$. R implements `pearson`, `kendall` and `spearman` measures of correlation in its function `cor`. So for each of these methods, if I compute               cor(X[1:n], X[2:(n+1)], method=...)      is there a known distribution for the resulting correlation coefficient? The `pearson` coefficient seems to be almost normally distributed, but for $n=10$ and method `spearman` it seems to be noticeably non-normal. Also, both `kendall` and `spearman` are discrete. I'm trying to figure out if in a long sequence $X_1,\ldots,X_{m(n+1)}$ of $m$ windows of size $n$, each window's consecutive random variables are independent. So I calculate the correlation coefficient for each window, and want to know its distribution under the hypothesis that they are all independent uniformly distributed as above. Are these distributions known? If not, are there any good approximations?


title: Sufficient statistic, specifics/intuition problems

text: I'm teaching myself some statistics for fun and I have some confusion regarding sufficient statistics. I'll write out my confusions in list format:   1. If a distribution has $n$ parameters then will it have $n$ sufficient statistics?   2. Is there any sort of direct correspondence between the sufficient statistics and the parameters? Or do the sufficient statistics just serve as a pool of "information" so that we can recreate the setting so we can calculate the same estimates for the parameters of the underlying distribution.    3. Do all distributions have sufficient statistics? ie. can the factorization theorem ever fail?   4. Using our sample of data, we assume a distribution that the data is most likely to be from and then can calculate estimates (e.g. the MLE) for the parameters for the distribution. Sufficient statistics are a way to be able to calculate the same estimates for the parameters without having to rely on the data itself, right?   5. Will all sets of sufficient statistics have a minimal sufficient statistic? This is the material which I am using to try to understand the topic matter: https://onlinecourses.science.psu.edu/stat414/node/283 From what I understand we have a factorization theorem which separates the joint distribution into two functions, but I do not understand how we are able to extract the sufficient statistic after factorizing the distribution into our functions.   1. The Poisson question given in this example had a clear factorization, but then it was stated that the sufficient statistics were the sample mean and the sample sum. How did we know that those were the sufficient statistics just by looking at the form of the first equation?    2. How is it possible to conduct the same MLE estimates using sufficient statistics if the second equation of the factorization result will sometimes depend on the data values $X_i$ themselves? For instance in the Poisson case the second function depended on the inverse of the product of the factorials of the data, and we would no longer have the data!   3. Why would the sample size $n$ not be a sufficient statistic, in relation to the Poisson example on the webpage? We would require $n$ to reconstruct certain parts of the first function so why is it not a sufficient statistic as well?


title: Determining if these random variables are independent

text: > Let $X$ and $Y$ be iid $N(0, \sigma^2)$ variables. > >   1. Show that $X^2+Y^2$ and $X/ \sqrt{X^2 + Y^2}$ are independent. **Hint** > : Compute the joint distribution of the first and the square of the second > expression. > Since $X \sim N(0, \sigma^2)$, we know that $X/\sigma \sim N(0,1)$. So $X^2/\sigma^2 \sim \chi^2(1)$. But this means that $X^2 = \sigma^2(X^2/\sigma^2) \sim \chi^2(\sigma^2)$. Of course, we have the same result for $Y$. Let $W_1 = X^2/(X^2 + Y^2)$ and $W_2 = X^2 + Y^2$. We get the joint distribution, $$f_{W_1, W_2}(w_1, w_2) = \left[ \left( \frac{1}{\Gamma(\sigma^2/2)2^{\sigma^2/2}} \right)^2 (w_1(1-w_1))^{\sigma^2/2 - 1} \right] \left[ (w_2)^{\sigma^2 - 2} e^{-w_2} \right]$$ There is a theorem in my textbook ( _Introduction to Mathematical Statistics_ By Hogg) which says > _Let the random variables $X_1$ and $X_2$ have supports $S_1$ and $S_2$ > respectively, and have the joint pdf $f(x_1, x_2)$. Then $X_1$ and $X_2$ are > independent if and only if $f(x_1, x_2)$ cna be written as a product of a > nonnegative function of $x_1$ and a nonnegative function of $x_2$. That is > $$f(x_1, x_2) \equiv g(x_1)h(x_2),$$ where $g(x_1) > 0$, $x_1 \in S_1$, zero > elsewhere, and $h(x_2) > 0$, $x_2 \in S_2$, zero elsewhere._ So, by this theorem, $W_1$ and $W_2$ are independent. (We know that $\left( \frac{1}{\Gamma(\sigma^2/2)2^{\sigma^2/2}} \right)^2 (w_1(1-w_1))^{\sigma^2/2 - 1}$ is positive because $w_1 = \frac{x^2}{x^2 + y^2} \leq 1$.) >   1. Show that $X+Y$ and $X-Y$ are independent. **Hint:** Use a standard > property of the normal distribution. > $X+Y \sim N(0,2)$ and $X-Y \sim N(0,2)$. Let $U_1 = X+Y$ and $U_2 = X-Y$. We have, $$f_{U_1, U_2}(u_1, ,u_2) = \left[ \left( \frac{1}{4 \pi} \right) e^{-u_1^2}{4} \right] \left[ e^{ (-u_2^2)/4} \right]$$ By the above theorem, $X+Y$ and $X-Y$ are independent. Do you think my answers are correct? Thanks you in advance


title: Can conditional Gaussian estimation problem solved by linear regression?

text: It just occurs to my mind about under what situation can the linear regression produce good result (valid). for example, suppose   y, $x_1$, $x_2$, ... $x_n$ follow multivariate Gaussian distribution, and now when $x_1$, $x_2$,... are observed and we are to estimate $y$, the most "reasonable" way is to just use the conditional Gaussian to estimate. However, from a different point of view, we can also solve it as a linear regression problem - suppose we have a group of samples (or previous observations) $x^{(1)}_1,x^{(1)}_2,... x^{(1)}_n, y^{(1)}$, and $x^{(2)}_1,x^{(2)}_2,... x^{(2)}_n, y^{(2)}$, ... Then we can forget about the correlations between $x_1, x_2, \ldots, x_n$, and then solve it as a linear regression problem. I wonder whether the results are the same. I doubt it because the second solution doesn't consider the correlations between $x$. Please advice. Thanks.


title: Ratio that accounts for different sample sizes

text: I'm trying to get a measure for a set of data that indicates someone's "spam score". Essentially, the higher the spam score, the more likely they are to be spammers. Right now, I'm measuring a person's spam score as a ratio: bad posts/total posts. If this ratio is high, then they're likely to be spammers. The difficulty I'm having is trying to compare different people - for example, a person with 6/8 bad posts is not AS likely to be as spammy as someone who has 600/800 bad posts (the person who's made 600 bad posts is clearly a spammer, but the one who's made 6 has not proven themselves to be a spammer to that extent). However, right now they are being assigned the same spam score of 6/8 = 0.75. Is there a way I can account for the size of the sample in my spam score?


title: Is there any relation between Power Law and Negative Binomial distribution?

text: In a social experiment that I was conducting, I was trying to count the number of people each user contacted in a period of 10 days. The population size was 100 for the experiment. Based on the values that I calculated, I fit a negative binomial distribution to the data (the Q-Q plot is given below). Conventional wisdom says that most networks amongst humans follow a power law distribution. I am guessing that my population size is too small to make a full conclusion about anything but is there any kind of relation between a negative binomial distribution and a power law distribution? I am asking this because I read a few days back that Normal Distribution and Gamma distribution (whose discrete analogue is the negative binomial) have a special role in that many other distributions can be derived from the Gamma distribution. I am wondering if this is true even with the power law distribution. I am a beginner in statistics so kindly point me in the right direction if I am out of track. ![alt text](http://i.stack.imgur.com/MMla7.png)


title: Is there a mathematical model that distinguishes between volatility and trend?

text: Say we are studying Twitter hashtags over time. We monitor how popular they are day to day. Some hashtags may be volatile (i.e. "lunch", "Celtics", "Friday"). Their popularity rises and falls frequently. Some hashtags may be in the process of becoming unpopular (i.e. "Gulf oil spill", "Transformers 2", "Christine O'Donnell"). Is there a mathematical model that can distinguish between a hashtag that has temporarily fallen in popularity but is likely to go up in popularity later and a hashtag that is sinking and likely to stay sunk in popularity? thanks


title: Computing the decision boundary of a linear SVM model

text: Given the support vectors of a linear SVM, how can I compute the equation of the decision boundary?


title: Newey West standard errors give me more significance

text: I had a time series model with 5 time series variables and it's a model that's reputed in the literature for having autocorrelation problems. Why when I use standard OLS var-covar matrix only 2 variables are (highly) significant but with Newey West varcovar matrix 4 of the variables are (highly) significant (1% level)?


title: (Quantile regression) AR(1) variable in the design matrix

text: I'm not doing a pure QAR (quantile auto regression) but I do have a lagged dependent variable (AR(1)) as a predictor. I'm using the `quantreg` package in `R` to do it. I have two very closely related questions: (1) I use `quantreg` in `R` to do my QRs. Do I need to do anything fancy to be fine by putting an AR(1) term on the RHS? (2) Is it fine to use AR(1) term as a regressor to deal with autocorrelation in quantile regression as you can with OLS? Or is QR different in this regards?


title: Conditional distribution of $X$ given $|X|$

text: If a continuous random variable $X$ has a symmetric distribution around 0, what is the conditional distribution of $X$ given $|X|$?


title: Why are Bayesian methods widely considered particularly "convenient"?

text: Several times I have come across (in lectures, papers, etc.) informal/tangential remarks that have the following basic form: > Philosophical issues aside, Bayesian methods are _extremely_ convenient. ...with the implication being that they are more "convenient" (or "handy", "easy", etc.) than competing ("frequentist"?) methods1. As one can tell from the assertion's preamble, such remarks usually come from people who would not describe themselves as "Bayesians". (This is what I meant to suggest by the adverb "widely" in this post's title: i.e. _even_ among "non-Bayesians", Bayesian methods are considered particularly convenient.) What is it about Bayesian methods that make them so convenient? The best answer to this question I can come up with is that Bayesian methods obviate the need for making some otherwise rather arbitrary decisions about parameters, and whole models. Instead, the Bayesian formalism serves as a principled way to simultaneously consider entire ensembles of such parameters and models (weighted by their corresponding priors). Now, my grasp of "Bayesian statistics" is tenuous at best, so I don't put much faith in this explanation. I would love to see a more informed elaboration. * * * 1 One example I can cite is a recent essay by Gelman and Shalizi's, where the authors disavow the "usual story" (i.e. the conventional interpretation of Bayesian statistics) as their reason for using Bayesian methods, and say instead they use them because with them they can "fit rich enough models" (presumably richer than they could fit using non-Bayesian methods). Unfortunately, this is the only hard reference I can provide. The other ones were in the form of informal lecture remarks, or in past reading that I can no longer remember specifically.


title: How to write a mathematical equation for GLM model with gamma and gaussian distribution?

text: I am writing a paper and the following is a code that I wrote in R. The reason that I am struggling with this is because I tried hundreds of models with different variables and the following model below had the best model with the least MAPE. I added and multiplied some variables and I am not sure how my model should be expressed in a mathematical equation and would appreciate any help. Thank you. FYI, B_1 is a categorical variable and the rest are numerical variables.               model = glm(as.formula(y ~ (B_1+B_2)*B_3*B_4*B_5), data = train.set, family="Gamma")       Thank you!


title: Using the eigenvalues from PCA in k-nearest-neighbours

text: I'm quite new to this StackExchange, only been a lurker till now, but my StackOverflow fellows have said you'd be the best people to ask about this. Anyway, enough introduction. I'm using the weighted k-Nearest-Neighbours algorithm. My original data set has 37 features. I've looked into using PCA to reduce dimensionality, and I'm going to follow this method. For simplicity's sake, let us assume that two of the new features created account for 90% of the variance and I'm only going to use these two new features. Let us call them feature 1 and feature 2 ($f_1$, $f_2$). Let us say that $f_1$ accounts for 60% of the variance and $f_2%$ accounts for 30% of the variance. I know wish to select the weights ($w_1 , w_2$) for these two features. My initial intuition is that we could correlate the variance accounted for with the weight of the feature. Therefore, I would use a weight combination of $w_1 = 0.6$ and $w_2 = 0.3$ in my k-Nearest-Neighbours algorithm. I am well aware that there is much literature suggesting the best way to select weights would be to use a lattice type of method where we select different combinations of weights and then follow through with combination that yields the best results. I was just wondering if the intuition of weights being related to total variance accounted for was. Also, as my dataset actually requires 11 features to be used to account for 90% of the variance, I'd like to have a starting point for determining combinations of weights. **Summary:** When using PCA as a precursor for kNN, is it possible to base the weights of features in k-NN on the total variance said features accounts for in the data? Sorry if there are any formatting errors or if I'm breaking any protocols. Let me know if I have, and I will update the post.


title: Census data comparisons

text: What is the best way to compare population of children age 0-6 years old (head counts) across various regions? Population density (population per sq km) seems to make sense. A colleague mentioned looking at population per capita. Does that make sense? Any help is much appreciated. Tx!


title: Estimate area of 95% confidence ellipse in R

text: I have data from a PCA and would like to know the area of the 95% confidence ellipse of the first two PCs. Does anyone know how to estimate the area of a 95% confidence ellipse in R?


title: F test and categorical variables

text: I am trying to compute a regression. I have a continuous dependent variable (investment in millions of dollars)   One continuous independent variables   Five dichotomous indicators (1 / 0)   Four interactions.   Robust standard errors are being used. N=569. Is there any particular reason why the F-Test is not being worked out in Stata? Stata help says: > Your estimation results show an F or chi2 model statistic reported to be > missing. Stata has done that so as to not be misleading, not because there > is something necessarily wrong with your model The problem seems to be related to the interactive dummy variables, because when I remove them the F test is calculated. The regression output is the following               Linear regression                                        Number of obs =  580                                                              F(  9,   569) =  .                                                              Prob > F      =  .                                                              R-squared     =  0.2336                                                              Root MSE      =  1.2608                                                         Robust     lninv                         Coef.       Std. Err.  t      P>t    [95% Conf. Interval]          age                           0.1392668   0.0363307  3.83   0.000    .0679081  0.2106255     1.idicatorone                 0.4070426   0.2380546  1.71   0.088   -.0605304  0.8746155     1.indicatortwo                1.085003    0.1235029  8.79   0.000    .8424262  1.327581     1.indicatorthree              -0.4390229  0.1781991  -2.46  0.014   -.7890312 -0.0890145     1.fincrisis                   0.4526393   0.1376538  3.29   0.001    .1822677  0.7230108     1.bubble                      1.657902    0.1088199  15.24  0.000    1.444165  1.87164     indicatorone#indicatortwo                                                              1 1                       0.0810565   0.2113363  0.38   0.701   -.3340379  0.496151      bubble#indicatorthree                                                                 1 1                       -0.8225963  0.2688304  -3.06  0.002   -1.350617 -0.2945752     fincrisis#indicatorthree                                                               1 1                       -0.1068982  0.2141135  -0.5   0.618   -.5274474  0.313651      indicatorthree#c.age                                                              1                             -0.0103719  0.0450553  -0.23  0.818   -.0988669  0.0781231      _cons                        0.0378515   0.1878819  0.2    0.840   -.3311753  0.4068783


title: Bagging, boosting and stacking in machine learning

text: What's the similarities and differences beetween this 3 methods: bagging, boosting, stacking? Which is the best one? And why? Can you give me an example for each?


title: Understanding MANOVA in case of a single predictor

text: I'm trying to understand the statistical analysis I saw in a clinical study. They measured performance of 3 groups of subjects with a series of performance measures (A, B, C, ..., N). The objective was to find out if there was a difference between the performance of the 3 groups, and if yes, then in which ways they differed. Their analysis consisted of doing a MANOVA like this:               set.seed(100)     group <- rep(c(0,1), each=40)     A <- rnorm(80, 5, .5) + .1 * group     B <- rnorm(80, 9, .3) + .2 * group + .5 * A     C <- rnorm(80, 12, .3) + .2 * group + .7 * B     d.1 <- data.frame(A = A, B = B, C = C, group = group)          fit.manova <- manova(cbind(A, B, C) ~ group, d.1)     summary(fit.manova, test="Pillai")                    Df  Pillai approx F num Df den Df    Pr(>F)         group      1 0.19669   6.2027      3     76 0.0007949 ***     Residuals 78      When they showed they see significant differences, they continued to perform a series of ANOVA tests for each DV, i.e.:               summary(lm(A ~ group, d.1))     summary(lm(B ~ group, d.1))     summary(lm(C ~ group, d.1))      to find out which performance measures where different between the groups. What I'd like to know is:   1. Is this approach (MANOVA followed by series of ANOVA) justified? Are there strict assumptions before we could take this path?    2. If yes, should there be some kind of correction for the second step, i.e. series of ANOVAs for individual DVs (multiple-comparisons)? What kid of correction?   3. What is the recommended approach for problems like this with multiple DVs? EDIT 1: changed the text to include an example code. EDIT 2: updated the example. DVs are now correlated. EDIT 3: this is actually a very common situation in studies involving objective measures. Devices typically just spit out an array of measures, even if you don't specifically ask for them. For example, you have two groups, a control and patient group, and do a laboratory gait analysis. The gait analysis systems gives you 50 different gait measures. Your research question might be: do the two groups have similar gait? If not, in which ways do they differ?