title: Fitting mixed effects logistic regression with random effects

text: I have a data frame of 134 observations, 9 independent variables, and a binary, categorical response; please see its structure below:               > str(df)     'data.frame':   134 obs. of  10 variables:      $ trip     : Factor w/ 11 levels "1","2","3","4",..: 2 2 2 2 2 2 2 2 2 2 ...      $ person   : Factor w/ 16 levels "1","2","3","4",..: 14 7 2 14 7 14 2 14 14 7 ...      $ f1       : Factor w/ 2 levels "1","2": 1 2 2 1 2 1 2 1 1 2 ...      $ f2       : num  146 163 179 190 194 172 172 175 240 214 ...      $ f3       : num  37 35 42 38 37 41 37 52 38 37 ...      $ condition: Factor w/ 4 levels "1","2","3","4": 1 1 4 1 4 2 2 1 2 1 ...      $ f4       : num  80 45 114 110 60 121 56 140 93 68 ...      $ f5       : num  60 136 128 136 165 118 220 90 177 240 ...      $ f6       : num  15.2 15.2 15.2 15.3 15.3 15.2 15.2 15.5 15.6 15.6 ...      $ mortality: Factor w/ 2 levels "1","2": 1 1 2 1 2 1 2 1 1 1 ...      For this particular study, we went out on several sampling trips where each person in our company attempted to capture animals with certain types of equipment. In addition to factors `f1 - f6`, we visually assessed the degree of physical injury done to the animal by each type of equipment and ranked it on a scale from 1-4. Then we tracked the animal after returning it to its environment to see if it died from being captured. I am attempting to fit a binary logistic regression with `mortality` as the binary response (1 = no, 2 = yes) and `f1,f2,f3,f4,f5,f6, and condition` as explanatory variables. Since I am posting on SO, it can only mean that I am having issues including the random effects of `trip` and `person`. More specifically, I am trying to incorporate these random effect terms with respect to `condition` since each `trip` and `person` could have introduced bias when `condition` was finalized. When we look at the frequency of animals caught and released on each `trip` and `condition`, for example, we see that each level of the factor `condition` is not present on each sampling `trip`:               > xtabs(~condition + trip, df)              trip     condition  1  2  3  4  5  6  7  8  9 10 11             1 17  9  4  4  1  2  4  4  4  4  4             2  5  9  6  4  0  3  2  8  8  0  0             3  1  0  1  4  1  5  3  0  0  5  0             4  1  3  2  3  1  0  1  0  1  0  0      This information above would make me believe that my `condition` data are nested within `trip`. However, since our data relied on the opportunistic capture of these animals, each `trip` had the same chance of catching an animal of each `condition` level (which would be a crossed-over design). So... **Question 1:** Are my random-effects of a nested or crossed over design? Depending on how individuals answer Question 1, I need to properly represent that random effect in an R package capable of handling mixed effects logistic regression. I am using the `ordinal` package since my final goal is to see which independent variables have the greatest influence on `mortality` so I can use them as predictors on animals we can't track upon release. Moreover, the `ordinal` package has the ability to predict mortality probabilities:               > library(ordinal)     > mod <- clmm(mortality ~ f1 + f2 + f3 + f4 + f5 + f6 + condition + randomEffects,                         data = df, Hess = TRUE)      Where `randomEffects` would be some combination of `trip`, `person`, and/or `condition` with respect to the intercept such as `(1|condition:trip)` for a `condition` being nested in `trip` or `(1|trip` if it was a crossed over. **Question 2:** How to properly represent the type of random effects in R syntax depending on if `condition`, `trip`, and/or `person` are needed. I understand the latter portion of this question is more relevant on StackOverflow but the the entirety of the question relates to properly setting up this mixed effects logistic regression. Please ask questions and I wouldn't mind suggestions about using other R packages. Thank you for your time.


title: How to decide which decision tree classifier to use?

text: I am confused about which decision tree algorithm in weka to use for my application. I have 5 real input variables and 2 classes. In various online tutorials J48 (C 4.5) seems to be the algorithm of choice. Are there any rules of thumb / tips / tricks to decide which tree algorithm should be used ? ![Weka decision tree choices](http://i.stack.imgur.com/MFIGx.jpg)


title: How to calculate the probability of matching sock from a drawer?

text: I'm not sure how to approach this problem from homework. The only equation I have for calculating probability is if each outcome in the sample space has an equal opportunity. The below isn't equal opportunity right? A drawer contains 6 red socks, 4 green socks, and 2 black socks. Two socks are chosen at random. What is the probability that they match?


title: How to calculate Prob > chi2 in R to test model fit of conditional logistic regression

text: I used the clogit function (from the survival package) to run a conditional logistic regression in R with a big dataset of 1:M matched pairs with n=300368964 and number of events= 39995.               model <- clogit(Alliance ~ OVB + CVC + BVB + strata(Strata), method="exact")          I received following results:                                coef  exp(coef)   se(coef)       z Pr(>|z|)         OVB        -0.0498174  0.9514031  0.0166275  -2.996  0.00273 **      BVB         0.0277405  1.0281289  0.0304956   0.910  0.36300         CVC         1.1709851  3.2251683  0.1089709  10.746  < 2e-16 ***     EarlyStage -1.3215824  0.2667129  0.0205851 -64.201  < 2e-16 ***     AvgVCSize   0.0087976  1.0088364  0.0002035  43.224  < 2e-16 ***     NumberVC    0.0643579  1.0664740  0.0034502  18.653  < 2e-16 ***     ---     Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Rsquare= 0   (max possible= 0.001 )     Likelihood ratio test= 6511  on 6 df,   p=0     Wald test            = 6471  on 6 df,   p=0     Score (logrank) test = 6801  on 6 df,   p=0      Since Rsquare equals 0 and the test ratios seems very high, I tried to plot the results to check whether the model fits. But I wasn't able to plot it properly. I would online many papers which use the ratio Prob > chi2 = 0 from Stata as test ratio to proof the model fit. How could I calculate this ratio in R? Are there any other ways I could check the model fit of my clogit results? I would appreciate any help. Thanks you very much in advance.


title: Is it possible to center orthogonal polynomials in multiple regression

text: I have a regression model that looks like the one below. $$ Y = \beta_0 + \beta_1T + \beta_2T^2 + \beta_3T^3 + \beta_4D + \beta_5D*T + \beta_6D*T^2 + \beta_7D*T^3 $$ Where `T` is a time variable, which starts from 0 - the starting point, and `T^2` and `T^3` are the natural polynomials of `T`; `D` is a dummy variable. `T`, `T^2`, and `T^3` are highly collinear. To reduce collinearity, one possible option is to use orthogonal polynomials. However, one issue I have with using orthogonal polynomials is that the intercept of the orthogonal polynomial regression is no longer the predicted value of `Y` at `T = 0` (and when `D = 0`). I crucially need the intercept to indicate the mean value of `Y` at `T = 0` because of theoretical importance. ### So my questions are: ### 1). Is it possible to center orthogonal polynomials in some way (e.g., centering around the minimum value??) so that the intercept will be equivalent to the mean value of `Y` when T= 0` But I am not sure how this can be properly done given that there are a quadratic term and a cubic term involved. ### 2). If the answer to the previous question is "No", I wonder if there are any alternative ways in which I can deal with multicollinearity while making the intercept mean what I want it to mean - namely B0 = mean of Y when T = 0 (and D = 0).


title: Alternatives to a train-test split with a small data set

text: Let's say I have a small data set (~50 observations) and 50 or so explanatory variables, and I want to develop a linear model for the sake of prediction. First, I want to construct a correlation matrix or a variable selection technique to select a smaller number of predictors to include within the model. However, I'm wondering what would be the appropriate action when attempting to train the model and checking whether it has predictive power? train-test data split? cross validation? what would be the best could of action with a small dataset and large number of explanatory variables? Tagged w/ R because I'll be using R to conduct analysis.


title: EM on product of multinomials

text: I have the following conditional density: $$ P(x | \theta, \pi) = \prod_{i=1}^I \prod_{j=1}^J t_{ij}! \prod_{k=1}^K \frac{1}{x_{ijk}!}(\sum_{l=1}^L \theta_{il} \pi_{jkl})^{x_{ijk}} $$ Here, $x$ is the data (observed). $\theta$ and $\pi$ are latent parameters that are unobserved. Finally, $t_{ij}$ is a constant that is calculated from the data. My overall goal is to infer the latent parameters ($\theta$ and $\pi$) given the data ($x$). Basically, $P(x | \theta, \pi)$ is a product of I*J multinomials. The parameters of each multinomial are controlled by $\theta$ and $\pi$. I'll briefly describe these variables. First, $\theta_i$ is a probability vector of length L. In other words, for all $i$, we have: $$ \sum_{l=1}^L \theta_{il} = 1 $$ The $\pi$ variables are indicator variables. For a given $j$ and $k$, they tell you which $\theta_{il}$'s to include in the sum. Intuitively, each multinomial is built from $K$ total parts. For each part $k$, we iterate through all possible $\theta_{il}$ and select which ones to include. There are a few simple constraints on the $\pi$ variables. For a given $j$ and $l$, there is only one $k$ for which $\pi_{jkl}$ is non-zero. It follows that: $$ \sum_{k=1}^K \sum_{l=1}^L \theta_{il} \pi_{jkl} = 1 $$ I hope this explanation is ok. Anyway, my goal is to infer the latent parameters $\theta$ and $\pi$. I was hoping I could do this with EM. My problem is that given $\pi$ and $x$, I don't know how to update the $\theta$ variables. I tried maximizing $P(x | \theta, \pi)$ over $\theta$, but I am running into trouble because the $\theta_{il}$'s are all dependent (because they must sum to 1). This makes calculating the derivative rather difficult. I was hoping that someone could give me advice or point me in the right direction. Also, if people can recommend other solutions, I would really appreciate it. Right now, I am looking into sampling from the posterior distribution of the parameters (given the data), but I need to do more reading to fully understand how this works. I can place a Dirichlet prior on $\theta$ and a discrete uniform prior on $\pi$, if it helps. Thanks!


title: Regression modeling of differences-in-differences in Stata

text: I'm having some trouble fitting a model for differences-in-differences regression in Stata. Specifically, I'm confused as to the `B3` coefficient in the below:               y = B0 + B1p + B2trt + B3(p*trt) + u      Here, `p` and `trt` are dummy variables where `p = 0` for pre-treatment and `1` for post-treatment, and `trt = 0` for control group and `1` for treatment group. `p * trt` is an interaction term indicating the marginal impact of the treatment group over the experimental period. My question is whether this should be a simple dummy variable equal to zero if either `p / trt = 0` and `1` if both `p / trt = 1`?


title: Particle filters and loopy belief propagation

text: I want to implement a loopy belief propagation algorithm for factor graphs with continuous variables and messages represented using particles, that is vectors of samples for an empirical distributions. All the works about particle filters that I have found consider either Markov Random Fields with unary and binary terms only or temporal sequences (Hidden Markov Models). In the case of a general factor graph, one needs to compute expressions like $$ m_{g_p\rightarrow z_q}(z_q)=\int_{z_1}\int_{z_2}\int_{z_3}g_p(z_1,z_2,z_3,z_q) m_{z_1\rightarrow g_p}(z_1)m_{z_2\rightarrow g_p}(z_2)m_{z_3\rightarrow g_p}(z_3) {\rm d}z_1{\rm d}z_2{\rm d}z_3. $$ This is the expression for the message going from the factor node $g_p$ to the variable node $z_q$ where $g_p(z_1,z_2,z_3,z_q)$ is function with known closed for expression and the messages $m_{z_i\rightarrow g_p}(z_i)$ are defined by particle distribution. All the equations I have see have only one of these particle terms; what I would like like to understand is the proper way to extend importance sampling and in general particle/MCMC algorithm to the case where you have products of particles. I know that $m_{z_i\rightarrow g_p}(z_i)$ can be written as a sum of weighted Diracs, however I am not sure about how to compute the products. Any pointers to relevant literature are highly appreciated.


title: Sparse PCA/Dictionary learning when the features are extremely sparse?

text: I am trying to do sparse PCA/dictionary learning, that is decompose a matrix $X\approx UV$ where the loading matrix $V$ is sparse, usually enforced with an $\ell_1$ penalty (the difference between sparse PCA and dictionary learning being whether the inner dimension of $U$ and $V$ is greater than the smaller dimension of $X$). In my data the columns of $X$ are extremely large (order $10^8$) but extremely sparse (order $10^2$). Are there specialized algorithms for the case of sparse data? Online is a plus.


title: Is there an effect size for Kappa's?

text: I am staring a project on comparing standard ways of creating a classifier with some heuristic methods. The heuristic methods should result in a faster training for the classifier but should result in a classifier with less accuracy - the usual tradeoff. I am running the comparison on some 50 datasets, and I want to report the difference on the accuracy by using the heuristic. 1) I was first planning on using Cohen's h on the difference of the accuracy. $$ h = 2 \arcsin(ac1) - 2 \arcsin(ac2) $$ EDIT: I found two conflicting formulas for the Cohen h. The formula above is in book Statistical Strategies for Small Sample Research edited by Rick H. Hoyle The second formula is $$ h = 2 \arcsin(\sqrt{ac1}) - 2 \arcsin(\sqrt{ac2}) $$ and the source seems to be Cohen himself at this link http://people.ucalgary.ca/~ramsay/cohen-effect-size-h-arcsin- transformation.htm. END EDIT (the intuition on the Cohen h is that it is much harder to improve the accuracy of a classifier from 0.90 to 0.95 than to improve another classifier from 0.70 to 0.75, that is, the effect size of the first change is larger than that of the second. Cohen's h theory actually talks about variance - the variance in the 0.90 is smaller than the variance around 0.7.) 2) But it hit me that to be even more fair, I must take into consideration how hard it is to have the accuracy of 0.9 or 0.70 to begin with. This is where the Kappa enters. Kappa measures how much above a "guessing" classifier is my accuracy (ac) $$ \kappa = {ac - guess\over{1-guess}} $$ where guess is the accuracy of guessing classifier that is aware of the proportion of the classes in the data $$ guess = pos*pos + neg*neg $$ where pos and neg are the proportion of positive and negative examples. 3) Here my first question. Is there an effect size measure for Kappa? Should I make the arcsin transformation or just use the average difference in kappa? 4) I also realized that I can use the variance of the accuracy or kappa for each classifier. I was originally planing in using a hold-out cross- validation, but if there is a use for the variance, I can use any other cross- validation that gives me an estimate of the variance (repeated hold-out, k-fold, bootstrap). In this case I could use a Cohen d to report the difference in accuracy (or Kappa). Does this make sense? Any thoughts? Thanks


title: Effective cutoff frequency for LOESS smoothing

text: I need to vary the window width (span) of a LOESS smoothing filter to select a given effective lowpass filter frequency. When I was using Savitzky-Golay for smoothing I found this technical report which related the window size to frequency, but have not found anything similar online for LOESS. The LOESS values will be similar to the S-G ones, but is there any more specific information available (for the Tricube or any other weight function)? Alternatively is there a calculation (dependent on the weight function) which can convert the S-G cutoff frequency to the LOESS?


title: Reporting $R^2$

text: Following my question here, **when is it appropriate (or inappropriate) to report r squared for a bivariate linear correlation.** As explained in the earlier question, r = 0.74 which results in r squared of 55%. The size of r is an excellent outcome but when r squared is taken into account (especially when 45% of the variance in the dependent variable remains unexplained), it sort of dilutes the message. As pointed out in the answer to the earlier question, I do understand that it is not readily possible for a single factor to explain all the variation. I think all statistical textbooks that I have consulted (unless I am consulting the wrong books!) state that r squared is a better way of understanding r (or the effect size). I do understand though that r of at least 0.71 is needed to get r squared of 50%. **Questions**   1. Is reporting both r and r squared a good practice (in social science research)?   2. When can I not report r squared? I have asked Q2 in order to to avoid a situation like the above, where my non statistical target audience may focus more on 45% of the unexplained variance, rather than on the 55% of the explained variance.


title: Weighted least squares with error correlation

text: I was wondering if errors have a nonzero correlation, is weighted least squares appropriate in this case?


title: How many features can be used for classification?

text: Asked a similar question the other day without an answer Link. I think maybe the question there is too big. Here I want to ask a specific one:   * 2 Class labels (Binary classification labelled with either 1 or 0);   * 50 Cases (20 labelled as 1 and 30 labelled as 0);   * 18 Features extracted (some physical quantities);   * 12 Features remained after feature selection; How many features can be used for classification to avoid overfitting? Or there is no limit? Some paper mentioned that the number of features should be less than the 1/3 of smallest labelled data. For my example, it should less than (1/3)*20 that is around 7 that means at most 6 or 7 features can be used. Am I right? Thanks a lot. A.


title: If Chebyshev's upper bound gives same value as the actual probability calculation, what can we conclude?

text: As an example, if Chebyshev tells $P(|X-\mu|\geq k\sigma)\leq 0.25$ and the actual probability for $k=2$ is also $0.25$.


title: Why bootstrap to calculate the standard error?

text: Can you please tell me the advantage of bootstrapping in the example below:               sampleOne   <- function(x) sample(x, replace = TRUE)     sampleMany  <- function(x, n) replicate(n, sampleOne(x), simplify = FALSE)     listMeans <- function(x, n) lapply(sampleMany(x, n), mean)     bootData <- function(x,n) do.call(rbind, listMeans(x,n))     sampleSize <- 100000     numBoots <- 1000          # Left Skewed distribution # shape1 = a and shape2 = b     set.seed(400)     popSkewLeft <- rbeta(sampleSize, shape1 = 5, shape2 = 1)     hist(popSkewLeft)     skewLeftbootData <- bootData(popSkewLeft, numBoots)     (populationMean <- mean(popSkewLeft))# Mean = a/(a+b)  = (5)/(5+1) = 0.8333333     (bootMean <- mean(skewLeftbootData))          (populationSd <- sd(popSkewLeft)) #sd = sqrt(ab/((a+b)^2 (a+b+1))) = sqrt((5*1)/((5+1)^2*(5+1+1))) = 0.140859     (bootSd <- sd(skewLeftbootData) * sqrt(sampleSize))      I created a left skewed population as can be seen from the above code. I also calculated the sample standard deviation. the population standard deviation was calculated using the beta distribution equation. The sample standard deviation and the booted standard error * the square root of the sample size are almost the same. As a matter of fact the sample standard deviation is closer to the population parameter.


title: Proof that $\mathrm{E}(s^2) = \sigma^2 \cdot N/(N-1)$

text: What's the derivation for expected value for sample variance for a sample taken from simple random sampling without replacement, i.e., how do we show that $$\mathrm{E}(s^2) = \sigma^2 \frac{N}{N-1}$$ Is my assumption this only applies to SRS samples without replacements correct? On a related note, how is the sample variance from two numbers, $x_1, x_2$ equal to $((x_1-x_2)^2)/2$ The same notation is used in this webpage/doc, http://www.utdallas.edu/~serfling/handouts/ExpValSampVar.pdf


title: Best way to get started with and learn R

text: I have tried several times to "go it on my own" - but with limited success. I am a casual SPSS user and some SAS experience. Would appreciate a pointer or two from someone who has similar background and now uses R.


title: Examples for One class SVM in R

text: I am trying to do one-class SVM in R. I have been trying to use e1071/ksvm kernlab package. But I am not sure if I am doing it correctly. Is there any working example for one-class SVM in R? Also,   * I am giving a big matrix of predictors as X. Since its supposed to be one-class, is the assumption that all training data I gave forms 'positive' class? If so, we don't have to give the labels 'Y'?   * The predicted labels given as output are True/False. So I am assuming, True is 'positive' class. Edit: Attaching sample code. Here I sampled 60% of 'TRUE' class and I tested on the full data set.               library(e1071)     library(caret)          data(iris)          iris$SpeciesClass[iris$Species=="versicolor"] <- "TRUE"     iris$SpeciesClass[iris$Species!="versicolor"] <- "FALSE"     trainPositive<-subset(iris,SpeciesClass=="TRUE")     inTrain<-createDataPartition(1:nrow(trainPositive),p=0.6,list=FALSE)     trainpredictors<-iris[inTrain,1:4]     testpredictors<-iris[,1:4]     testLabels<-iris[,6]          svm.model<-svm(trainpredictors,y=NULL,                    type='one-classification',                    nu=0.5,                    scale=TRUE,                    kernel="radial")     svm.pred<-predict(svm.model,testpredictors)     confusionMatrixTable<-table(Predicted=svm.pred,Reference=testLabels)     confusionMatrix(confusionMatrixTable,positive='TRUE')


title: Paired t-test to support increase in student reading scores from fall to spring

text: I did a pre-test, post-test method and I tried to perform a t-test to show how reading scores increased from fall to spring. I’m not sure if my numbers are accurate. What throws me off is the NULL Hypothesis. According to my calculation using Excel it is $7.37 \times 10^{-27}$ or 0.00000000000000000000000000737 (for a paired, 2 tail t-test). Is this realistic to get? My understanding is that you have to achieve <.05 to reject the null hypothesis. And my calculated t-tests are really small. Does this mean I can reject the null hypothesis and say that what I did to increase reading scores was effective? Also, some students who took a test in fall withdrew from my class mid-year and didn’t take a spring test. Do I leave those cells blank or do I fill them in with zeros? Or do I just pretend that those kids never existed? How much would that throw off the results? Thanks.


title: How do I convert the variances explained by a principal component into an eigenvalue?

text: I'm using SigmaPlot to run PCA on various measurements that should all correspond to size of an animal. Running PCA using a covariance matrix (instead of a correlation matrix, since all of the measurements are measured in the same units) yielded really, really small eigenvalues for the first three PCs, despite relatively high proportion of variance explained. I've read that proportion of variance explained should correspond to eigenvalue, but I'm not sure how (and I can't seem to find out). Given the proportion of variance explained, how do I calculate the eigenvalues (i.e., confirm that SigmaPlot is not making some sort of error). Thanks!


title: Two sample $t$ test in R: where did it go wrong?

text: I am doing two sampled t test on two species of fish.               t.test(fish1,fish2)      It runs for a year so there are 365 values for each fish, i.e. $(1,0,1,0,0,2...)$ and $(0,0,01,21,0,5...)$. Yet the $df =368.586, t=2.2205, p=0.02699$, and the 95% confidence interval: $(0.1454645, 2.3970013)$. The sample means are               Mean of x    mean of y     1.5452055.   0.2739726      Considering the df is larger than the amount of values I have I think it's gone wrong somewhere.


title: How to get only positive values when imputing data?

text: Suppose age is normally distributed with mean 20 and standard deviation 5. How do you ensure that you get only positive values when you sample age from this distribution? I am trying to impute missing data. I am assuming that it is MAR (missing at random) since the missingness depends on gender. The missing ages come from both males and females. Thus I am sampling from two normal distributions (one for males and one for females).


title: How to model non-linear time courses with nlme

text: I have some time course data which plotted looks like this: ![enter image description here](http://i.stack.imgur.com/W2YE7.png) I am trying to fir a model to it using `nlme.lme()`. I am interfacing with R via RPy, and using the following code:               nlme.lme(r.formula('Pupil~CoI*Time'),random=r.formula('~1|ID'),correlation=nlme.corCAR1(),data=dfr)      so my formula is `Pupil~CoI*Time` my random effects are `~1|ID` and my correlation structure is `corCAR1()` My **questions** are:   * How do I tweak my formula (and random effects?) to obtain a non-linear time course as seen in my data? Currently I am getting values for first degree functions - which is not helpful. I have no theoretical basis on which to specifically add square roots or exponents and so on to the formula, all I know is linear is not enough.   * How do I decide what correlation structure to use to account for autocorrelation? I chose this because the description from the ETHZ R manual mentioned the word autocorrelation, but sadly I do not understand much from those brief descriptions.


title: Is it important to convert “integer” variables (with 0 or 1 values) to factors?

text: I am working on a high-dimensional dataset (1776 variables). When I read the csv file, R loads variables (with 0 or 1 values) as class of "integer". Is it important to convert these variables to factors before building predictive models (e.g gbm, randomForest, SVM, etc.)? What if I want to use PCA to reduce the dimension, is it correct to keep these variables as "integer"? (since PCA does not work with factor variables). Cheers


title: Using more measurements of lower quality or just the one with the best quality?

text: I am currently trying to analyze some pupillometry data. The device I use gives me X and Y measurements of the pupil area (from which I compute the total area via $A = \pi \times \left(\frac{d}{2}\right)^2$ ). The X-Y plot of my data looks like this: ![enter image description here](http://i.stack.imgur.com/0ey0a.png) In gray you can see my raw data points. As you see, though this effect does not seem to strongly affect participant means, there is a slight deviation from the regression line in the upper right quadrant of the figure. **My question is** : in such a situation, would $d=X$ or $d=\frac{X+Y}{2}$ give me a more reliable measurement?


title: Filtering with HMM

text: I want to use HMM for filtering, i.e. to find $p(x_t|y_{1:t})$. I see that the forward algorithm calculates the forward variable as a joint probability; $\alpha_t(i) = p(y_{1:t},x_t=S_i|\lambda)$, meaning joint of observations so far and the current state is $S_i$, given HMM parameters $\lambda$. To find $p(x_t|y_{1:t})$, it looks fine to sum all alphas at time t to find the marginal $p(y_{1:t}|\lambda) = \sum_{i=1}^{N}\alpha_t(i)$. Now, if I normalize $\alpha_t(i)$ for each of N possible states with this marginal, the resulting distribution should be $p(x_t|y_{1:t})$. Does that sound fine?


title: Handling redundant factor variable levels for linear regressions in R

text: Say I have two factor variables, X and Y, each with 3 levels. However, X==3 if and only if Y==3, while such a connection doesn't hold for X,Y==1,2. In this case, while X and Y are not redundant, my design matrix, when expanded, will contain the same column twice (0's everywhere, except 1's where X==Y==3). This is not good. Theoretically, there is no problem: Just remove one of the redundant columns. However, I don't want to go ahead and re-implement the entire linear regression just for this case. Is there some technical or theoretical trick that will still allow me to use R's `lm`?


title: Can a different indexing approach improve the relevance and efficiency of search engine results?

text: I am creating a search engine, with the corpus consisting of websites crawled through a webcrawler (Apache Nutch). I need the query searches to be both fast and relevant. So far, I have been trying to develop a Latent Semantic Indexing system, but I recently read an article that told me that LSA had not shown to be useful to them. Basically, LSA uses a k-rank approximation of the singular value decomposition of a matrix whose dimensions are the documents and a weighted term frequency (IDF, TF-IDF, Binary, and Raw). Does anyone know a much better data-mining algorithm that I could use? Through your own personal experience? I only have about 2 weeks left for this project, and I need to make sure I can end it right.


title: Shouldn't the beta estimates always be the same using lm and lmer?

text: Can anyone tell me under what conditions the beta estimates differ between lm and lmer with a random intercept? I came across a situation where the fixed effect differed considerably. I thought the std errors should change but the fixed effects should remain unchanged. The difference does not seem to be due to having different cluster sizes or having a large number of clusters with a single observation. I cannot supply the data but have concocted a simplified example below. In this case the correlation within clusters should be negligible.               library(lme4)     x=c(rep(0,10),rep(1,10))     y=rnorm(length(x),mean=3100,sd=400)-200*x     m=c(1,2,3,4,4,6,7,8,8,10,11,12,13,14,15,16,17,18,20,20)     summary(lm(y~x))     summary(lmer(y~x+(1|m)))


title: How to know the stochastic gradient descent is converging

text: How can I know if the stochastic gradient descent algorithm is converging. I cannot plot my objective function and take its average for lets say every 1000 iterations to see the trend. My objective function itself is huge, so it takes lot of time to compute the objective function itself. In this situation, how to determine if the algorithm is converging. I know what the parameter values should be. But I see the algorithm is fluctuating around a value which is smaller than the actual one. My function has a global maxima. So it should reach that maxima isn't it? The function is concave


title: What test should I choose if I want to see if two groups are different from each other in many categories?

text: I am currently analysing data of a big webshop with over 31,000 vendors selling goods, i.e. each vendor is selling items over the webshop. The cool thing is that we know the gender of these vendors (=individuals) and the categories of their products. Now, we want to see if there are categorical gender differences, i.e. we would like to test, if males sell more/less items than females of a certain category and if these differences are significant. Here is a graph of the data: ![http://i.stack.imgur.com/sitxf.png.](http://i.stack.imgur.com/I8srB.png) The question I have now is which statistical test is appropriate for this kind of analysis? So far I was using a chi-squared test to see if the distributions over all categories are significantly different between males and females. Now I would like to know: which statistical test is valid to use to show me for instance that the two distributions between males, females of the category "Apparel" is also significant?


title: Differences between Stata & SPSS in multiple regression outcomes

text: I'm a graduate student working on a thesis whose aim is to test complementarity between 2 practices (CI and INNO, according to theory of supermodularity). To do so, I calculated two new variables (each as a mean of 6 other variables) corresponding to the 2 practices, then derived from them 4 binary variables based on the combination of the two (that is, I used medians of the two variables as a threshold, then created four categories like 11,10,01,00: `HH,LH,HL,LL`). I took these four binaries, (plus two continuous variables and three other dummies as controls) and ran a multiple regression on `Y` (operational performance), without a constant term. I used both SPSS and Stata, but obtained two very different outcomes. I don't understand why!!! Moreover, correlation matrices gave the same results. Here are all the outcomes: **SPSS**               Variabili escluse a,b     Modello Beta In t   Sig.    Correlazioni parziali   Statistiche di collinearità                         Tolleranza  VIF Tolleranza minima     1   LL  .c  .   .   .   ,000    .   ,000     a. Variabile dipendente: PERF_medio     b. Regressione lineare che passa per l'origine     c. Predittori nel modello : industryT, industryM, industryE, HL, LH, HH,                                  log_age, log_size              Coefficienti non standardizzati                 B     STD ERROR DEV         LH       ,188     ,110         HL       ,156     ,115         HH       ,467     ,102         log_size   ,008   ,040         log_age   -,003   ,039         industryE 3,416   ,244         industryM 3,487   ,255         industryT 3,551   ,246      As you can see, the `LL` variable is excluded; Stata doesn't exclude this. Moreover multicollinearity tests give positive results, whereas VIF performed with Stata after regression do not. **STATA:**               regress PERF_medio LL LH HL HH log_size log_age industryE industryM industryT, noconstant     note: industryM omitted because of collinearity                  Source  SS       df     MS               Number of obs    =     190                                                      F(  8,   182)  =   1253.59             Model   2654.3682   8   331.796025           Prob > F           =    0.0000             Residual    48.1709704  182 .264675662   R-squared  =    0.9822                                                      Adj R-squared  =    0.9814             Total   2702.53917  190 14.2238904       Root MSE           =    .51447              PERF_medio |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]             -------------+----------------------------------------------------------------               LL |   3.486898   .2548016    13.68   0.000     2.984153    3.989643               LH |   3.675309   .2746033    13.38   0.000     3.133493    4.217124               HL |   3.643261   .2646596    13.77   0.000     3.121066    4.165457               HH |   3.954367   .2753474    14.36   0.000     3.411083     4.49765         log_size |    .008443   .0400578     0.21   0.833    -.0705944    .0874803          log_age |  -.0027773   .0393044    -0.07   0.944    -.0803283    .0747737         industryE |  -.0706952   .0935333    -0.76   0.451    -.2552442    .1138538         industryM |          0  (omitted)         industryT |   .0639002   .0960992     0.66   0.507    -.1257117     .253512              estat vif, uncentered              Variable |       VIF       1/VIF            -------------+----------------------         log_size |     42.17    0.023711               HH |     16.90    0.059170               LL |     15.21    0.065753          log_age |     12.70    0.078726               LH |     10.26    0.097499               HL |      8.73    0.114505         industryT |      2.23    0.447811         industryE |      2.21    0.451551         -------------+----------------------             Mean VIF |     13.80      Thus, SPSS excludes `LL`, but Stata does not. What's wrong? Why are coefficients so different?


title: How to justify using a gamma distribution for biological image analysis

text: I am trying to characterize objects within the nucleus of a cell. These objects are brighter than the background, and when I remove the background and map the frequency distribution of the objects I know that the distribution fits a gamma or negative binomial distribution. The appropriate way to choose a model is understand the basis of a model and logically apply it to your situation. Could anybody give me insight or give me direction to use either model?


title: Model averaging binomial glm with a categorical explanatory term

text: I am a using model averaging approach (using the MuMIn package) with a binomial GLM. One of my explanatory variable is a categorical factor 'treatment' with 5 levels (T1,T2,T3,T4 and T5). I also have a 2 binary explanatory variables Gender and Employed. I am also including all 2 way interactions .               glm(y ~ Gender*Employed+Gender*Treatment+Treatment*Employed,data=data1,family="binomial")      I am aware that because their is an interaction I must center my binary variables. I am not so sure what to do with my categorical variable though. Should I make the treatment variable into 5 dummy variables (0/1) and center these by subtracting the mean? Then include them in the model as below? In the model below T5 would be my reference so is not included in the model.               glm(y ~ Gender*Employed+Gender*T1+T1*Employed+Gender*T2+T2*Employed+Gender*T3             +T3*EmployedGender*T4+T4*Employed,             data=data1,family="binomial")       If I do this I am aware that the top models may not include all of the levels of the Treatment. Does this still make sense? Also this model compares all levels T5 to the other levels. If I wanted to know how all levels of treatment differed from each other would I need to make a total of 5 models, each time including all but one level of Treatment?


title: How to calculate the bounds for the probability of exceedance while eliciting from experts?

text: I intend to elicit rank correlation using probability of exceedance from domain expert to generate conditional probability table for a Bayesian Network. The probability of exceedance to be elicited is contrained by the values estimated previously. I need to know the lower & upper bounds for the exceedance probability for subsequent elicitations. How does one go about calculating this? Thank you


title: How can I analytically prove that randomly dividing an amount results in an exponential distribution (of e.g. income and wealth)?

text: In this current article in SCIENCE the following is being proposed: > Suppose you randomly divide 500 million in income among 10,000 people. > There's only one way to give everyone an equal, 50,000 share. So if you're > doling out earnings randomly, equality is extremely unlikely. But there are > countless ways to give a few people a lot of cash and many people a little > or nothing. In fact, given all the ways you could divvy out income, most of > them produce an exponential distribution of income. I have done this with the following R code which seems to reaffirm the result:               library(MASS)          w <- 500000000 #wealth     p <- 10000 #people          d <- diff(c(0,sort(runif(p-1,max=w)),w)) #wealth-distribution     h <- hist(d, col="red", main="Exponential decline", freq = FALSE, breaks = 45, xlim = c(0, quantile(d, 0.99)))          fit <- fitdistr(d,"exponential")     curve(dexp(x, rate = fit$estimate), col = "black", type="p", pch=16, add = TRUE)      ![enter image description here](http://i.stack.imgur.com/L9Qw4.png) **My question**   How can I analytically prove that the resulting distribution is indeed exponential? **Addendum**   Thank you for your answers and comments. I have thought about the problem and came up with the following intuitive reasoning. Basically the following happens (Beware: oversimplification ahead): You kind of go along the amount and toss a (biased) coin. Every time you get e.g. heads you divide the amount. You distribute the resulting partitions. In the discrete case the coin tossing follows a binomial distribution, the partitions are geometrically distributed. The continuous analogues are the poisson distribution and the exponential distribution respectively! (By the same reasoning it also becomes intuitively clear why the geometrical and the exponential distribution have the property of memorylessness - because the coin doesn't have a memory either).


title: Why does auto.arima() give negative output?

text: This is the dataset on which I am working currently, which is production data. Data:               > test.ts             Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct    Nov    Dec     1990                                                            0.0   10.8  180.0  418.2     1991  561.9  517.9  531.3  448.1  254.9   49.0    3.2    0.0    0.0   10.4  207.7  526.2     1992  597.2  581.5  596.4  518.4  378.3  209.9   32.1    0.0    0.0    7.9  166.7  571.7     1993  650.4  578.5  565.7  280.5   35.3    0.7    0.0    0.0    0.0   39.5  289.2  638.9     1994  643.8  533.8  410.9  159.3    0.0    0.0    0.0    0.0    0.0   38.3  322.8  684.9     1995  695.6  665.8  640.2  415.4  113.0   20.7   12.1    0.0    0.0   13.6  316.3  677.5     1996  754.5  683.4  719.6    0.0    0.0    0.0    0.0    0.0    0.0    0.0    0.0  678.0     1997  774.5  808.1  847.9  677.8  208.7    9.9    0.0    0.0    0.0    5.2  296.4  794.9     1998  952.0  873.1  732.0  264.6    3.9    0.0    0.0    0.0    0.0    0.0  245.8  833.0     1999  843.5  812.3  708.3  275.2   10.8    0.0    0.0    0.0    0.0    5.4  300.4  884.9     2000  949.0  898.6  892.7  474.7  130.0   19.8    0.0    0.0    0.0    8.5  367.2 1000.8     2001 1092.1  987.7  864.3  392.8   41.8    0.0    0.0    0.0    0.0    7.0  425.0  968.0     2002  983.6  925.0 1018.0  696.0  209.0   26.0    0.0    0.0    0.0    0.0   63.0  823.0     2003 1066.0  930.0  929.0 1071.0  614.0  125.0   29.0    0.0    0.0    0.5  300.0 1005.8     2004 1043.0 1051.9  863.2  279.1    8.0    0.0    0.0    0.0    0.0   67.8  597.1 1120.3     2005 1087.9 1015.9  855.7  292.9    0.8    0.0    0.0    0.0    0.0   78.3  683.3 1139.2     2006 1185.5 1162.1 1131.3  386.9   16.4    1.2    0.0    0.0    0.0    7.1  728.5 1493.0     2007 1572.9 1341.0 1652.9 1279.3  386.4   14.3    0.0    0.0    0.0    0.0  102.5 1570.1     2008 1864.7 1786.7 1523.9  422.7   48.1    0.8    0.0    0.0    0.0    0.0  192.4 1556.9     2009 1260.8  763.8  284.1    6.1    0.0    0.0    0.0    0.0    0.0    0.0   73.8 1495.6     2010 1280.8 1248.8  887.2  185.6    7.3    0.0    0.0    0.0    0.0    0.8  182.0 1524.9     2011 1461.5 1497.7 1111.5  108.6    0.0    0.0    0.0    0.0    0.0    2.9  519.3 1652.5     2012 1552.5 1563.2 1380.4  295.2    7.7    0.0    0.0    0.0    0.0    0.1  225.0 1677.6     2013 1686.2 1420.0 1691.0  795.0    0.0        I used auto.arima() from forecast package. Code:               ARIMAfit <- auto.arima(test.ts)     test.ar <- forecast(ARIMAfit, level=70, h=12)      Following is the output I got Output:               > test.ar              Point Forecast       Lo 70     Hi 70     Jun 2013      -4.429870  -186.37952  177.5198     Jul 2013      -4.429870  -261.74553  252.8858     Aug 2013      -4.429870  -319.57590  310.7162     Sep 2013      -4.429870  -368.32916  359.4694     Oct 2013      -3.416802  -410.26858  403.4350     Nov 2013     296.121405  -149.56239  741.8052     Dec 2013    1505.197792  1023.80428 1986.5913     Jan 2014    1477.195886   962.56457 1991.8272     Feb 2014    1327.574562   781.72562 1873.4235     Mar 2014    1423.251183   847.87588 1998.6265     Apr 2014     550.206881   -53.25183 1153.6656     May 2014      -1.892754  -632.18482  628.3993      Questions:   1. Why does the output shows negative value, when there has been no negative value in the historic data? The data is of production, which cannot be negative.   2. Is there any other model class which handles "zero" values appropriately? Kindly help.


title: A basic voting system question

text: I am sorry for the vague title - its because I don't know what I am talking about.   Imagine a committee of people in charge of hiring a CEO. A committee member can vote "No hire" (+0), "Maybe"(+1), "Hire"(+2) for a potential CEO candidate. Each CEO is scored based on the votes, and the CEO with the highest score is hired. Now, imagine a situation where two candidates are tied... both of them have a score of 5. However, candidate1 had 2Hires and 1Maybe (ie, 2+2+1), whereas candidate2 had 5Maybes (ie, 1+1+1+1+1). Here is my question:   * Intuitively, I think we can say that we are more confident about candidate1 eventho the scores are tied. What specific domain of statistics am I dealing with here (so that I can investigate on my own)   * In statistics, how does one go about resolving such situations? And I am not just merely talking about resolving a tie... but rather, how to resolve this in a generic manner? Thank you.


title: Can composite scores of Likert scale be obtained?

text: I am using a Likert scale for a group of 3 questions.   1. Strongly Agree   2. Agree   3. Undecided   4. Disagree   5. Strongly Disagree I have 3 questions with these possible responses under the category "Effectiveness". Can the responses in the 3 questions be combined to give a meaningful response? For example: the responses in Q1, Q2, Q3 are 2, 3, 5. Can the three responses be combined in some way? Can a composite score be obtained of these?


title: Graphing a Probability Curve for a Logit Model With Multiple Predictors

text: I have the following probability function: $$\text{Prob} = \frac{1}{1 + e^{-z}}$$ where $$z = B_0 + B_1X_1 + \dots + B_nX_n.$$ My model looks like $$\Pr(Y=1) = \frac{1}{1 + \exp\left(-[-3.92 + 0.014\times(\text{bid})]\right)}$$ This is visualized via a probability curve which looks like the one below. ![enter image description here](http://i.stack.imgur.com/JjQzf.png) Now, I am consider adding a couple variables to my original regression equation. So let's say I add gender (categorical: F and M) and age (categorical: < 25 and > 26) into the model. So I end up with: $$\Pr(Y=1) = \frac{1}{1 + \exp\left(-[-3.92 + 0.014\times(\text{bid}) + 0.25\times(\text{gender}) + 0.15\times(\text{age})]\right)}$$ Now, in R I can generate a similar probability curve which will tell me the probability of Y=1 when accounting for all three predictors. Where I'm lost is I want to find the probabilities for every possible permutation of these variation. So when bid = 1, gender = M, and age is >= 26, what is the probability that Y = 1. Similarly, when bid = 2, gender = F, and age is >= 26, what is the probability that Y = 1. I want to generate a probability curve which will allow me to visualize this. Can anyone help? I may be completely misunderstanding what kind of information one can glean from a logit model, but please tell me if I am misunderstanding the theory also.


title: How do I compare the accuracy of two measurement devices when one of them is the reference?

text: I'm currently trying to compare the measurement accuracy of two devices. However, I am using one of the devices as a "gold standard" to say something about the accuracy of the other device. The quantity I am trying to measure is the distance (depth) to a point in a scene (image). I am using a laser scanner which gives me (x,y,z) coordinates for a point and I am using optical geometry to give me (x,y,z) coordinates for the same point. My base (or reference) measurement is the laser scanner measurement. I am using this as a gold standard. I now have a corresponding measurement from the optical geometry based measurement system. I have corresponding measurements for 1000's of points. I average the error in the difference between these and get a mean error. This error, however, neglects the fact that the laser scanner too is a measurement system and has its own error. It comes with an accuracy of +/- 5mm for an object at 10m range. How do I incorporate this into an accuracy metric for the optical geometry based system?


title: Logistic regression weights of uncorrelated predictors

text: I am a bit puzzled about the behavior of uncorrelated predictors in logistic regression. As in OLS, I thought that if two predictors (`rv1` and `rv2`) are uncorrelated, then the regression weights of `rv1` will not change from a regression that only includes `rv1` to one that includes `rv1` and `rv2`. However, it seems to be the case that this is not true in logistic regression and coefficients change between the two regression models, even if the predictors are uncorrelated. I have pasted some R syntax below that demonstrates this behavior. Why is this the case and how do the regression weights from the two regressions (the one with only rv1 and the other one with rv1 and rv2) relate to each other? Is there a way to know what the regression weight of rv1 will be if one knows the regression weight of rv1 in the regression that includes both predictors? Thanks! P.S. This post is crossposted at another unrelated stat answer site.               library(MASS)          #generate lots of data (a little bit weird data handling, I know)     n <- 10000     rdta <- as.data.frame(mvrnorm(n=n,c(0,0),matrix(c(1,0,0,1),2,2),empirical=TRUE))     names(rdta) <- c("rv1","rv2")          #confirm that preds are uncorrelated     cov(rdta$rv1,rdta$rv2)          rv1 <- rdta$rv1         rv2 <- rdta$rv2          rv1ry <- 1     rv2ry <- 1          #generate binary data from known regression coefficients     ylinp <- (1 / (1+exp(-(-1 + rv1*rv1ry + rv2*rv2ry))))     y <- rbinom(n,1,ylinp)      glm(y~rv1+rv2,family=binomial(link='logit'))     glm(y~rv1,family=binomial(link='logit'))     glm(y~rv2,family=binomial(link='logit'))          #confirm that OLS regression works as expected (regression weights do not change)     rv1y <- .222     rv2y <- .333     y <- rv1y * rv1 + rv2y * rv2 + rnorm(n,0,.5)     lm(y~rv1+rv2)     lm(y~rv1)     lm(y~rv2)         I am not sure if it is expected to also paste relevant output here, but here goes: OLS results > lm(y~rv1+rv2)               Call:     lm(formula = y ~ rv1 + rv2)          Coefficients:     (Intercept)          rv1          rv2        0.001096     0.220051     0.333072        > lm(y~rv1)               Call:     lm(formula = y ~ rv1)          Coefficients:     (Intercept)          rv1        0.001096     0.220051        > lm(y~rv2)               Call:     lm(formula = y ~ rv2)          Coefficients:     (Intercept)          rv2        0.001096     0.333072        Logistic regression results > glm(ry~rv1+rv2,family=binomial(link='logit'))               Call:  glm(formula = ry ~ rv1 + rv2, family = binomial(link = "logit"))          Coefficients:     (Intercept)          rv1          rv2            -1.001        1.916        2.469        > glm(ry~rv1,family=binomial(link='logit'))               Call:  glm(formula = ry ~ rv1, family = binomial(link = "logit"))          Coefficients:     (Intercept)          rv1           -0.5495       1.0535        > glm(ry~rv2,family=binomial(link='logit'))               Call:  glm(formula = ry ~ rv2, family = binomial(link = "logit"))          Coefficients:     (Intercept)          rv2       -0.6538       1.6140


title: Origin of "5$\sigma$" threshold for accepting evidence in particle physics?

text: News reports say that CERN will announce tomorrow that the Higgs boson has been experimentally detected with 5$\sigma$ evidence. According to that article: > 5$\sigma$ equates to a 99.99994% chance that the data the CMS and ATLAS > detectors are seeing aren’t just random noise — and a 0.00006% chance that > they’ve been hoodwinked; 5$\sigma$ is the necessary certainty for something > to be officially labeled a scientific “discovery.” This isn't super rigorous, but it seems to say that physicists use standard "hypothesis testing" statistical methodology, setting $\alpha$ to $0.0000006$, which corresponds to $z=5$ (two-tailed)? Or is there some other meaning? In much of science, of course, setting alpha to 0.05 is done routinely. This would be equivalent to "two-$\sigma$" evidence, although I've never heard of it being called that. Are there other fields (besides particle physics) where a much stricter definition of alpha is standard? Anyone know a reference for how the five-$\sigma$ rule got accepted by particle physics? **Update:** I'm asking this question for a simple reason. My book Intuitive Biostatistics (like most stats books) has a section that explains how arbitrary the usual "P<0.05" rule is. I'd like to add this example of a scientific field where a much (much!) smaller value of $\alpha$ is considered necessary. But if the example is actually more complicated, with use of Bayesian methods (as some comments below suggest), then it wouldn't be quite apt or would require a lot more explanation.


title: What are some common measures of variation in an interval?

text: Context: I'm working on an app that "listens" to breathing/snoring and "diagnoses" apnea by taking note of pauses between breaths. For the most part this works well, but it's (rather amazingly) difficult, in the app, to distinguish snoring/breathing from talking/music (as with a TV playing in the room). I've pretty well exhausted signal processing techniques to distinguish this, so I'm trying some basic statistical measures. I find that simply taking the normalized standard deviation of the breath/snore interval is fairly effective (values below roughly 0.35 are "good"), but I'm wondering if there might be another measure of variation, other than standard deviation, that would produce better results (a crisper dividing line). Basically, I want to distinguish "snore, snore, snore, pause, snore" from a random string of intervals with the same mean (such as you would get with voice), keeping in mind that snoring/breathing is not exactly clockwork regular. So what other (reasonably simple) measures of variation in a train of numbers are there? **Update:** It occurs to me that a measure of skewness might be useful here. Music/voice would presumably tend to have less skew than breathing/snoring. How does one measure skewness (particularly with some sort of moving average)? **Update 2:** I've placed some sample interval lists out in Dropbox:   * Voice: http://dl.dropbox.com/u/72976325/intervalData_holm.csv   * Music: http://dl.dropbox.com/u/72976325/intervalData_epica.csv   * Snoring sample #1: http://dl.dropbox.com/u/72976325/intervalData_P1005.csv   * Snoring sample #2: http://dl.dropbox.com/u/72976325/intervalData_P1009.csv


title: How to run a logistic regression with multiple related dependent variables?

text: **Background** In a small online study I asked people (n=100) which products they would purchase. The choice set contained 20 products for which they had to indicate whether they would buy a product or not (= DV, yes/no format). Respondents could indicate yes for any number of products. The IVs were product and presentation characteristics. **Data format** resp_id (in rows), buy_prod_1, buy_prod_2 ... and IVs in columns **Previous analyses** I did a “simple” logistic regression in Stata. For this I reshaped my data so that there is one line for each respondent and each product (100*20 lines, buy_prod and IVs in columns). In the logistic regression buy_prod was my DV. **Problem** I do not know how to integrate the fact that the product choices for each respondent are not independent, i.e. how to combine the different DVs in the right way in one model.


title: Relation between mean of the hypergeometric distribution and binomial

text: The mean of the hypergeometric distribution is: $n \frac{K}{N}$ where:   * $n$ is the number of draws   * $K$ is the number of successes    * N is the size of the finite population. As the population size, $N$, goes to infinity, I would expect to have a situation with "replacement", and have the mean to converge to the mean of the Binomial distribution, but instead the mean converges to 0 (see the formula). Why is this the case?


title: Data mining techniques in Obama's campaign

text: I came across this article about the data mining team in Obama's reelection campaign. Unfortunately, the article is very fuzzy about the actual machinery of the statistical algorithms. However, it sounded as if the general techniques are known in social and political sciences. Since this is not my area of expertise, can anyone point me at (overview) literature about these kind of techniques?


title: What are the advantages of Multiple Kernel Learning (MKL) methods?

text: Multiple Kernel Learning methods aim to construct a kernel model where the kernel is a linear combination of fixed base kernels. Learning the kernel then consists of learning the weighting coefficients for each base kernel, rather than optimising the kernel parameters of a single kernel. The disadvantages of multiple kernel learning seem to be that they are less interpretable and computationally expensive (as to evaluate the model output you need to evaluate all of the base kernels). So if similar performance can be achieved by simply optimising a single kernel, what are the advantages of MKL?


title: Meaning of intercept in multinomial regression with binary predictors?

text: I am doing a multinomial regression and trying to interpret the results: In the basic model there is only one binary predictor variable (0 = high risk scenario, 1 = low risk scenario), the dependent variable has 3 categories (strategy 1,2 or 3). The output shows that the model is significant and most of the logits as well. However, I am wondering how to interpret the (significant) logits of the intercepts. This is how I am interpreting it: "The likelihood to choose strategy 1 over 2 in a risk-free scenario (the intercept?) is higher/lower, keeping the rest constant." It this the correct way to interpret it? Thanks in advance.


title: How to assess the quality of a return period estimate?

text: **Background** I have a time series of structural loads, which are measured forces on a moored ocean buoy, and I need to obtain the return period value so that the structure can be designed to withstand the max load expected in a storm. The return period is estimated from the measured time series in the following way: First, I choose a threshold and identify the peaks above this threshold. Extreme events are expected to have a distribution known as the Generalized Pareto Distribution (GPD), so I fit the peaks to that distribution. Then the GPD fit is inverted to give the load associated with the probability of a chosen time period (3 hours). This load is called the return period value. I am using a program called WAFO, which takes extreme value analysis routines from S-Plus, to do this analysis. **Problem** The problem is how to assess the quality of the return period estimate. WAFO produces quantile-quantile plots that show how well the peaks match the GPD, a p-value to check the quality of the fit, or confidence bounds on the return period estimate. But these 3 diagnostics seem to conflict. Sometimes the q-q plot shows that many peaks don't fit the GPD (bad), but the p-value for the fit is high (good) and the confidence bounds are narrow (good). Sometimes the confidence bounds are very wide (bad) but the q-q plot shows the peaks lining up with the GPD (good). My gut sense is that the quality might be poor because we only measured data for a limited duration and there is likely measurement noise. **Question** : How can I tell if the return period estimate is good vs bad?


title: Load forecast model with temperature data

text: I would like to estimate a regression model of the type: $ load_t = seasonality_t + trend_t + \beta * temperature_t, $ and I have load data and temperature on high frequency (hourly data). My impression is that the temperature does not influence load at the same frequency as I measure load (i.e. a change in temperature for 2 or 3 hours does not imply a change in load immediately). This is how I understand the application of the concept of coherency as in http://eprints.nuim.ie/1968/1/JR_C81dfisf.pdf I tried to look at average temperature per day and use this in the regression but the results were not satisfying. How can we incorporate temperature into a practical model in the best way? Any hints? Good references?


title: How to best visualize differences in many proportions across three groups?

text: I'm trying to visually compare how three different news publications cover different topics (determined through an LDA topic model). I have two related methods for doing so, but have received lots of feedback from colleagues that this isn't very intuitive. I'm hoping someone out there has a better idea for visualizing this. In the first graph, I show the proportions of each topic in each publication, like so: ![Proportions for all topics and publications](http://i.stack.imgur.com/sZ7XC.png) This is pretty straightforward and intuitive to almost everyone I've talked to. However, it's difficult to see the differences between the publications. Which newspaper covers which topic more? To get at this, I graphed the difference between the publication with the highest and second highest proportion of topics, colored by the publication with the highest. Like this: ![Difference between first and second highest topics](http://i.stack.imgur.com/mgB7p.png) So, the huge bar for football, for example, is really the distance between al- Ahram English and Daily News Egypt (#2 in football coverage), and it is colored red because Al-Ahram is #1. Similarly, trials is green because Egypt Independent has the highest proportion, and the bar size is the distance between Egypt Independent and Daily News Egypt (#2 again). The fact that I have to explain that all in two paragraphs is a pretty sure sign that the graph fails the self-sufficiency test. It's hard to tell what's really going on by just looking at it. Any general suggestions about how to visually highlight the dominant publication for each topic in a more intuitive way? **Edit: Data to play with:** Here's `dput` output from R, as well as a CSV file. **Edit 2:** Here's a preliminary dot plot version, with the diameters of the dots proportional to the proportion of the topic in the corpus (which is how the topics were originally sorted). Though I still need to tweak it a little more, it feels a lot more intuitive than what I was doing before. Thanks everyone! ![Dot plot](http://i.stack.imgur.com/2pc2N.png)


title: Regression with percentage points

text: Suppose I want to find out how growth in average income affects the growth in income of the poor. For this I have 60 countries and for each country the income of the poor and average income for 6 years. I am unsure on how to perform a regression of this. I have been calculating the average annual growth in percentage for average income and income for the poor. From what I understand this gives me percentage points change instead of percentage change, thus it would be a simple linear regression. Yet I have seen some tutorials that log the percentage points, creating a log- log regression. Which of these two methods are correct? EDIT: Example values: Country A:   * Average income for the poor for 4 years: 1990:150, 1991: 260, 1992:300, 1993: 400   * Average income for 4 years: 1990:340, 1991:600, 1992:710, 1993:1000 Compound annual growth poor: 38.67%   Compound annual growth all: 43.28% Say I do this for 60 countries, so i will get for average annual growth of income of the poor (in percent): 21.10 , 16.23 , 12.34 , 23.25 , 11.12 , 13.55 etc and for average annual growth if income of all: 43.13 , 35.31 , 56.31 etc and I do a regression on them and get say y = 2.433 + 0.6 x . So if income of all increases one unit, the income of the poor increases 0.6 units. Is it correct to say than that a 1 percent change increase in average income of all will result in 0.6 percent change average income of the poor? Or is it percentage points change?


title: Significant variables, but no improvement on predictive ability of model

text: When running binary logistic regression in SPSS, I obtain first a measure of the predictive ability of a model without independent variables (base model). This is then compared with that of the model containing independent variables. If there is no difference between the predictive ability of the models, is any significance shown by the variables therefore irrelevant or does it contain some value?


title: How do we interpret model priors if models intersect?

text: Let $\Theta$ be some parameter space and $\Theta_1,\Theta_2,\dots,\Theta_s \subset \Theta$ be parameter subsets which represent competing models in some model selection procedure. There are no assumptions about containment or disjointness of the $\Theta_i$ made. The fundamental postulate of Bayesian model selection is that posterior probability of model $\Theta_i$ is computed according to Bayes rule $$P(\Theta_i | D) \propto P(\Theta_i) \int_{\Theta_i} L (\theta | D) dQ_i(\theta) $$ where $D$ stands for the data, $L$ the likelihood, $Q_i$ the conditional prior distribution of model $\Theta_i$, and $P(\Theta_i)$ the prior. My question is: **What do we mean by the prior probability $P(\Theta_i)$?** The prior should express our a priori belief of $\Theta_i$ being the right model in some sense, so one might come up with something like $\text{Prob}(\theta_0 \in \Theta_i)$, where $\theta_0$ is the true parameter generating the data (assuming one for now). However, this definition does not work because the true parameter may well be contained in several models. Another way to ask the question: **Assume I know the true parameter $\theta_0$, which prior should I choose?** Do we have to let go of the notion of an unknown true parameter? Is modeling with overlapping sets $\Theta_i$ inconsistent?


title: How to compare two independent KS tests?

text: I have two vectors of observed data `obs1` and `obs2`, and two vectors of expected data `exp1` and `exp2`. All the vectors have different sizes.   I performed two independent one-sided Kolmogorov-Smirnov (KS) tests: `ks.test(obs1,exp1, alternative='greater')` and `ks.test(obs2,exp2,alternative='greater')`, to assess if `obs1` and `obs2` were significantly greater than the respective expected vectors `exp1` and `exp2`.   Now I am interested to see which of the two `obs` cases shows a distribution significantly shifted towards lower values, that is 'greater' in terms of KS test.   Basically the idea is to compare `obs1` and `obs2` taking into account their 'differences' from their respective `exp` distributions (I tried two-tailed KS test between `obs1` and `obs2` but in this way I don't consider their 'deviations' from `exp1` and `exp2` respectively). Does somebody have any suggestion? Thanks everyone!


title: Formula for "adjusted mean" - Is this correct?

text: I am working on a software project involving statistical analysis of myoelectrical signals (EMG). One client requested the calculation of "adjusted average" for a given signal interval. His description is the following: $$\text {AdjustedMean} = \text{Mean} + \frac{\text{StandardDeviation}^2}{\text{Mean}}$$ He sent me this explanation about it: > "The adjusted mean is used to correct for situations with high standard > deviation, such as those found with Electromyography." I googled a lot, but couldn't find any reference to such formula. Most references to "adjusted mean" are related do ANCOVA, and then the mean is adjusted to some other variable, not relative to the standard deviation of the same variable. (EDIT, as per WHuber request) My client has these other considerations about what he intends to do with the adjusted average: > "The adjusted mean presents a better linear correlation with the exerted > muscular force (measured by a dynamometer) than the unadjusted average. This > is more evident for stronger contractions / forces, because then more motor > units are recruited and thus the electrical values oscillate more widely". Some other paper I've read DEFINES the EMG amplitude envelope as the (moving?) Standard Deviation itself, a somewhat different definition of the more common rectification + smoothing.


title: validation of a Zero Adjusted Gamma model

text: I am using a Zero Adjusted Gamma regression (two-part ZAGA) model to estimate the effect of a psychometric categorical factor ("expected recovery", with 3 levels) on cost associated with treatment after whiplash injury. I have about 15 covariates/degrees of freedom that are relevant to the cost outcome and potential confounders. The ZAGA two-part model fits a logistic regression for the probability of having zero cost, and a gamma model (using log link) for those with cost outcome > 0. I understand that there are limitations with this type of models for the objective of interest, but for various reasons we need to stick with it for now. For the logistic part, we have only 19 observations with zero cost, so we cannot use many covariates. It also unclear how the measured covariates are associated with the probability of having zero costs (i.e. deciding not to pursue any further treatment options). So we ended up using a "two-step" stepwise regression approach, including in the logistic model only the two covariates that had significant univariate association with I(cost=0). On the other hand, we were able to include all covariates in the gamma model (we have more than 300 non zero cost observations). No two covariates had correlation > 0.8 so we decided to keep them all in the model. We would like to properly validate our model, keeping in mind that the main objective is the effect of (potentially modifiable) factor of perception of "expected recovery" and 1-yr health care cost (vs. prediction). It is unclear to me, i) whether we could have done a better job in the model selection process, ii) what would be the appropriate ways to validate such two-part model and potentially apply shrinkage to tackle possible over- fitting. Any ideas/suggestions/comments would be greatly appreciated. Thanks, Nicholas


title: Standard deviation exceeds maximum value

text: I am trying to plot data that only has two possible values – `1` (present) or `0` (absent) – clearly not a normal distribution. The max value could be 1, but my standard deviation bars go to 1.3, which is impossible. I tried increasing the sample size, but the max column value + standard deviation are still 1.3. Is it OK to manipulate my graph and cap my column bars off at 1 since it is not possible to have a sample value larger than 1?


title: Longitudinal paired analysis

text: I need your advice regarding a complicated design. I am testing data regarding new eye drops that suppose to decrease some quantitative measure. For every subject, one eye is randomly assigned with the drops, while the other is a control, and prior applying the drops, this measure (Y from now on) is measures on both eyes. From here I have calculated the difference, which had a mean near 0, with SD of 1.5, and the points looked unbiased around the mean. Now the drops were applied, and the measure was taken again on both eyes, and the difference was calculated. Now it was far from 0; the treatment works. After some days, another measure was taken, showing that the effect is decreasing. And then there was another visit of measuring. I want to find the time when the treatment isn't working anymore. How should I do it? I have for every subject something like 4 visits including the baseline, and I have the difference between the eyes on each visit. I can't simply say that below 1.5 SD still counts, or any other cutoff; I need some test. Which test / model should I use here? Thank you!


title: How to interpret canonical correlation analysis

text: I have done a canonical correlation analysis using the American Community Survey Dataset. The analysis is done between `Ancestry` and `Educational Attainment` variables. The values for canonical correlations are: {0.8140, 0.5716, 0.4708, 0.3946, 0.1465, 0.1365, 0.0409}. The values for multivariate tests of significance for the first canonical function:               Statistic        Value    F-Value    df1    df2    P-Value     Wilk's Lambda    0.143    35.373     98    19754   0.000     Pillai's Trace   1.408    54.358     98    21896   0.000     Hotelling Trace  2.961    22.779     98    21842   0.000     Roy's Root       0.662      The redundancy analysis for the first canonical function: For `Educational Attainment` variables               Average Squared Loading      Canonical R^2      Redundancy Index        2.119                    0.666               1.411      For `Ancestry` variables               Average Squared Loading      Canonical R^2      Redundancy Index        2.1514                    0.666               1.432      Please help me interpret these results. I understand the p-values are good (as they are very low); _F_ -values seems good as well (as they are high), but I'm not sure how to interpret the values of Wilk's Lambda, Pillai's Trace, the Hotelling–Lawley Trace, and Roy's Largest Root. According to my understanding of the analysis, the multivariate tests are statistically significant, but not practically significant according to redundancy analysis. Any help (a little explanation of the analysis) would be greatly appreciated. Thanks.


title: $I(X:Y|Z=z) \leq I(X:Y|Z)$?

text: $X,Y$ and $Z$ are discrete random variables. $I$ is the mutual information. The question is in the title. If the inequality is true, how would you show it? Thanks. My intuition: having more information on the conditional variable $Z$ ($e.g.$ $Z=z$) decreases the correlations between $X$ and $Y$, as visualized by the decrease in the intersection surface.


title: How to simulate a etas spatio-temporal model?

text: My question is how simulate the spatio-temporal ETAS point process (a point process used for earthquake prediction) defined by a conditional intensity function. More precisely: let $\mathbf{V}$ be a spatio temporal point process, $\mathbf{s}$ be a vector of dimension 2 containing the coordinates, $t$ be the temporal index and $m$ indicates a threshold level for the earthquake. In the ETAS model the intensity (conditional) function is \begin{equation} \lambda(\mathbf{s},t,m|\mathcal{H}(t))= \lambda(\mathbf{s},t|\mathcal{H}(t))\Gamma(m) = \left(\mu(\mathbf{s})+ \rho \sum_{i:t_i<t}\alpha(m_i)g(t-t_i)f(||\mathbf{s}-\mathbf{s}_i || | m_i)\right)\Gamma(m) \end{equation} where $\mathcal{H}(t)$ is the history of earthquake until time $t$, $t$ excluded. $\mu(\mathbf{s})$ is the background intensity and could be function on coordinates and some covariates: \begin{equation} \mu(\mathbf{s}) = \mathbf{Z} \boldsymbol{\eta} \end{equation} where $\mathbf{Z}$ is the matrix of covariate (possible with $\mathbf{s}$) and $\boldsymbol{\eta}$ is the vector of regressor coefficients. $\alpha(m)$ is the expected number of events triggered from an event of magnitude $m$: \begin{equation} \alpha(m) = \exp(\psi (m-\bar{m})) \end{equation} where $\bar{m}$ is a threshold magnitude (events over this threshold are not take into consideration in the model). $g(t-t_i)$ is the probability density function of the occurrence times of the triggered events: \begin{equation} g(t-t_i) = \frac{\beta_2-1}{\beta_1}\left(1+\frac{t}{\beta_1}\right)^{-\beta_2} \end{equation} The locations distribution of the triggered events we use is \begin{equation} f(||\mathbf{s}-\mathbf{s}_i || | m_i) = \frac{\xi-1}{\pi \delta^2\exp(\omega (m-\bar{m}))} \exp\left( 1+ \frac{(\mathbf{s}-\mathbf{s}_i)^{\prime}(\mathbf{s}-\mathbf{s}_i)}{\delta^2\exp(\omega (m-\bar{m}))} \right)^{- \xi} \end{equation} other choice are possible. Finally \begin{equation} \Gamma(m) = \gamma \exp (\gamma (m-\bar{m})) \end{equation} Then the log-likelihood is \begin{equation} \log L(\mathbf{V}| \mathcal{H}(t_1)) = \sum_{i=1}^n \lambda(\mathbf{s}_i,t_i,m_i|\mathcal{H}(t_i))-\int_{ \mathcal{T}_B} \int_{ \mathcal{S}_B} \int_{\mathcal{M}_b} \lambda\left( \mathbf{s}, t ,m| \mathcal{H}(t) \right) d \mathbf{s} d t d m \end{equation}


title: Verifying representativeness of a sample, after simple random sampling

text: I just used the standard formula to determine the sample size of a sample to match the mean of a pupulation with an error margin of 3 percentage points and with a 90% probability. I know would like to check that the mean of my sample is actually within the 3% margin of error of the mean in the population. To do so I am taking the population and appending the sample, identifying observations in the sample with a dummy called sample. Then I run the following regresion: reg var1 sample And I am testing wheter the coefficient for the dummy variable sample is less than 3 percentage points. To do so I am using a one sided t test. But this is only ok when the coefficient of the dummy variable "sample" has the "right" sign. I would like to test the joing hypoithesis that the coefficient of the dummy variable sample is greater than -3 and lower than 3. Does anyone have any insight on how to do this test? Thanks in advance.


title: Inverse Differencing and ARIMA Model Equivalence

text: I've developed a ARIMA model with exogenous variable. Before fitting the model, I made every time series stationary by differencing (each variable had a different order of integration). For simplicity, let's say there was only one exogenous variable. Finally I ended up with a model. Keeping in mind that all my data are stationary to create the model, I received the following coefficient estimates:               ARIMA(1,0,0) with zero mean               Coefficients:            ar1      X           -0.6     -0.002      There were no MA terms. Note that Y was second-differenced, and X was first differenced. So, my formula is really like: $Y_t'' = -0.6 Y_{t-1}'' - 0.002 X'$ Where the primes ( $'$ ) represent the level of differencing. I'd like an equivalent representation of the model using the original non- differenced variables. So I perform the following: $(Y_t - Y_{t-1}) - (Y_{t-1} - Y_{t-2}) = -0.6 [(Y_{t-1} - Y_{t-2}) - (Y_{t-2} - Y_{t-3})] - 0.002(X_{t} - X_{t-1}) $ $Y_{t} - 2Y_{t-1} + Y_{t-2} = -0.6Y_{t-1} + 1.2Y_{t-2} - 0.6Y_{t-3} - 0.002X_{t} + 0.002X_{t-1}$ $Y_{t} = 1.4Y_{t-1} - 0.2Y_{t-2} - 0.6Y_{t-3} - 0.002X_{t} + 0.002X_{t-1}$ So, I have 2 questions. 1) Is my interpretation correct? Seems straight-forward, but I just want to double check. Is there a quicker way to perform these conversions, besides the tedious expanding & combining I performed above? Some of the models I'm working with are quite more complicated. Perhaps using the backshift operator would help? I don't have much practice with algebra involving the backshift operator, and am not sure whether it would help do what I'd like. 2) How do I perform a similar procedure when there are MA terms? I know that an MA term is equivalent to infinite AR terms. So, is it still be possible to get some representation like $Y_t = ...$? For example, suppose my fit produced ARIMA(1,0,1), with an MA coefficient estimate of 0.3. Thank you!


title: Imputation using caret package vs imputation without it

text: In the ` caret ` package, there is the option to impute missing data using decision trees. I already imputed my data by sampling the missing data from a distribution. Also, I determined the important variables to include in the model instead of relying on the ` caret ` package. In general, is it better to clean the data by yourself instead of relying on the ` caret ` package?


title: How to understand the log marginal likelihood of a Gaussian Process?

text: I'm trying to understand Gaussian Processes. Could anyone tell me:   1. Why we need to use the log marginal likelihood?   2. Why using log, the marginal likelihood can be decomposed to 3 terms (including a fitting term and a penalty term)?


title: In Gaussian Processes, how to understand the hyper-parameters optimization?

text: I know that in a GP, hyper-parameters are optimized by maximizing the marginal likelihood. Could anyone explain this method to me please? Thanks for your help in advance.


title: Two clocks problem and measurements

text: > If I have one clock, I know what the time is. > > If I have two clocks, I now am unsure. Although simplified, the problem is related to a more complex system I am testing. If I have the following measurements from two systems:               Attribute   System 1   System 2     A           0          0     B           10         100     C           90         100     D           100        110     E           900        1000      then I'd like a way to determine some value of error between the systems (assuming I cannot distinguish which system, if any, is correct) So   * attribute A values are considered no-error or both systems correct   * attribute B values are considered a larger error than A, C, D and E values   * attributes C, D and E values are have the same error (or similar error)  Also, assume the attributes are independent. My understanding of statistics is low, and I can only recall some things about relative error, but iirc I require knowing whether one system is correct to determine the relative error of the other system. How does one go about resolving this issue?


title: Comparing distributions as opposed to medians?

text: Assume a hypothetical scenario of two events. During event 1, I observe a set of values $[X_1, X_2,X_3,...,X_n]$. A physical phenomena occurs and this triggers event 2 for a short period of time and I end up observing another set of values $[Y_1, Y_2,Y_3,...,Y_n]$ during that time. I don't know anything about the underlying distributions. I have a number of $X-Y$ pairs and I want a method to detect when the Y values are generally larger than the $X$ values for further investigation. The pairs are not comparable - that is I cannot use the $X$ values from another subject to inform me about the distribution, however, in the null case I expect $X$ and $Y$ should be from the approximately the same distribution. My first attempt is to define event 2 to have some importance if the median of the observed values from event 2 is higher than that of the observed values in event 1. I have a few questions:   * Because median is known to be robust to outliers, what are the pitfalls of such a comparison i.e., when will this fail?    * If this is not good, then can I make use of something like KL-Divergence or JS-Divergence to understand the extent of divergence (and if it is beyond a threshold, conclude that event 2 is important)?   * Does comparing divergence have any advantages over the median? I think I want to detect significant change as an event. Any clarification will be greatly appreciated.


title: Variability of data

text: There is a student, who is taking only 2 classes in one semester. Class 1: Math, Class 2: Science. For each of the classes, every Thursday a quiz is taken. There are no other exams or hw; only the quizzes. At the end, if we look at his math performance, we can draw autocorrelation verses week lag; and find correlation length. In a word, we can say whether his math performance was consistent throughout the semester or not (assuming that, if the correlation length was greater than some value, we will call the performance to be consistent). We can do the same thing, for his science class. Now, we want to answer, how consistent was his all over performance?


title: Does smoothing of Bayes classifier will increase precision?

text: I have implemented Bayes multinominal and Bernoulli's model and my question is does the smoothing have any impact of the performance of both models (Laplace’s law of succession or add one smoothing)?


title: Cholesky versus eigendecomposition for drawing samples from a multivariate normal distribution

text: I would like to draw a sample $\mathbf{x} \sim N\left(\mathbf{0}, \mathbf{\Sigma} \right)$. Wikipedia suggests either using a Cholesky or Eigendecomposition, i.e. $ \mathbf{\Sigma} = \mathbf{D}_1\mathbf{D}_1^T $ or $ \mathbf{\Sigma} = \mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^T $ And hence the sample can be drawn via: $ \mathbf{x} = \mathbf{D}_1 \mathbf{v} $ or $ \mathbf{x} = \mathbf{Q}\sqrt{\mathbf{\Lambda}} \mathbf{v} $ where $ \mathbf{v} \sim N\left(\mathbf{0}, \mathbf{I} \right) $ Wikipedia suggests that they are both equally good for generating samples, but the Cholesky method has the faster computation time. Is this true? Especially numerically when using a monte-carlo method, where the variances along the diagonals may differ by several orders of magnitude? Is there any formal analysis on this problem?


title: Linear Regression: Ordinal or dummy independent variables?

text: I am a beginner who needs some help reading some results from linear regressions. I am looking for factors that influence the location of civil conflict events. My dependent variable is the distance from the capital (km), a continuous variable. I have 4 control variables which are all continuous too. The variable (fightcap) I want to test is the fight capability of the rebels. This ordinal variable has three levels:   1 = low   2 = moderate   3 = high   However, I read in a couple of textbooks that to deal with ordinal predictors it is better to convert them to dummy variables. I thus created two dummies:   fightCap(low): 1 if fightcap = 1, 0 otherwise   fightCap(moderate): 1 if fightCap = 2, 0 otherwise Using JMP (the software I know best), I ran the two models in parallel: one (left) with the fightcap ordinal variable, the second (right) with both fightCap(Low) and fightcap(moderate) dummy variables. ![enter image description here](http://i.stack.imgur.com/Bu9Dt.png) I notice that the output is identical for both models, excepts regarding the estimates of the ordinal/dummy variables. **My questions are thus the following:**   1. How do I interpret the terms "fightcap[2-1]", "fightcap[3-2]" and their estimates? Do I read that a change of "fightcap" from 1 to 2 results in a 59.42 change in the dependent variable?   2. How do I interpret the fact that "fightcap[2-1]" has a significant estimate, but not "fightcap[3-2]"? Does it mean that only the variation from 1 to 2 explains the variation, but not the variation from 2 to 3?   3. How can I understand that "half" of my variable works if using the ordinal setup, but none of the dummies has a significant estimate? Thanks for your help! Damien


title: Relative Importance using Out-of-Sample Error

text: I'm currently running some regressions to understand the nature of price changes in a financial market. In this market, price changes can occur for two reasons: local price changes (which affect the specific user) and market-wide price changes (which affect all users). Let $P$ denote the price change, $M$ denote the market-wide price change, $L$ denote the local price change, and subscript $i$ denote trade number. Then, the equation: $P_i = M_i + L_i$ holds EXACTLY for each trade $i$. What I would like to understand, however, is how strongly the local effects $L$ influence the $P$ series. To do this, I have used an ordinary (linear) least squares regression to fit the coefficients $c$ and $k$ for the regression line $P = kL + c$, which is equivalent to assuming that the relationship for a give trade is given by $P_i = kL_i + E_i$, where $E_i$ is a noise term. At present, I perform this regression, then estimate the out-of-sample error of the regression (using 10-fold cross validation). I then repeat this whole process for a similar regression that includes only the market-wide term, i.e., $P = jM + d$. I again estimate the out-of-sample error of the regression using 10-fold cross validation. I then compare how important the $M$ and $L$ factors are in causing price changes by comparing the out-of-sample error of these two regressions. The reservation that I have about this method is that I already know the values of $E_i$ in the first regression (they are exactly equal to $M_i$ in each case), and I similarly know that the noise terms in the second regression are exactly equal to $L_i$ for each trade. Therefore, I am unsure about whether using a regression in which I "pretend" not to know them is really telling me anything. Is there a standard (frequentist) way to approach this type of problem (preferrably by comparing out-of-sample performance, rather than $R^2$)? Thanks!


title: Mixed-model or not? Case-study in R

text: I performed some observations of clownfishes. For each individual fish (each observation), I recorded the anemone species where it lives, the size of the anemonea and the size of the fish (and many other things but we'll forget them for the purpose of this question). In total I have 8 fishes and 7 anemones. It represents the total number of species present in the sampling area. I want to know if the **size of the anemonae** correlates with the **size of the fish**. But I have to decorrelate the **fish species** and the **anemonea species** from this analysis as both are maybe correlated with the **size of the fish**. What model should I run on R? May I run the following...               aov(fish.size ~ factor(fish.species) + factor(anemone.species) + anemone.size)      and record only the p.value for the `anemone.size`? Or should I use a mixed-model-effect?               lmer(fish.size ~ anemone.size + 1 | factor(fish.species) + 1 |         factor(anemone.species))      or something else?..


title: What diagnostics for random effects logistic regression?

text: I'm wondering what checks/diagnostics you usually calculate and report if you do a random effects logistic regression.   * C-statistic/ROC curve   * Check multicollinearity?   * Check heteroscedasticity?   * McKelvey & Zalvoina R-squared   * ... There is a lot of information about standard logistic regression and what you should check for, but I'm missing an overview about what applies to logistic regression on panel data. Thanks


title: expressing a univariate normal as a multivariate normal

text: I need to express a univariate normal as a multivariate normal to make certain calculations possible (for example: being able to divide two Gaussian distributions). So, my univariate normal is defined over only one random variable ex: $G(x)$ and now I want to define it over $G(x_1, x_2.....x_n)$. However, even though it is defined over $n$ random variables, it does not depend on the values of the other random variables. Can someone tell me how I can define the mean and covariance of such a distribution, so it is mathematically correct? Does the mean vector only have non-zero value at one location or should I try and describe $(x-\mu)$ as 0? What about the covariance matrix?


title: from z-score to p-value - how to do this for this hypothesis-set?

text: I know that if               H0: μ = 1.35       then if `H1: μ != 1.35` means p-value is 2P(Z ≥ |z|) else if               H1: μ > 1.35 then P(Z ≥ z) and H1: μ < 1.35 means p-value is P(Z ≤ z)      However, I am stumped about how to do it with following sets of hypothesis.               H0: μ ≥ 1.35      H1: μ < 1.35       ## p-value=??               H0: μ ≤ 1.35      H1: μ > 1.35       p-value=??


title: Clustered data WITHOUT multilevel / GEE model?

text: I have a data-set with around 700 observations from 12 centres. Although the clustering effect as tested in a random intercept model didn't seem significant, it seems more appropriate to use a multilevel / GEE model. But it may not be easy to do so as I want to test for multiple potential mediators. (to be specific, the model in my mind is a SEM without latent variables). So my question is: is it theoretically advisable to fit the data in a single- level path model using ML (Huber sandwich estimator, or bootstrap SE)? UPDATE: Thanks for those helping out. Just a update for what I did yesterday. I tried the following models:   1. Single-level SEM in R (lavaan package)   2. Single-level SEM in MPlus    3. Two-level SEM in Mplus (1) and (2) were fitted with ML (Huber / Bootstrap SE). (3) was fitted with Bayesian estimator as MLR did not converge. Basically all of the above gave me similar results but (1) and (2) gave me the chance to do multi-group analysis on top of the SEM and (3) seems to be more theoretically correct. I am not quite sure if the results from (1) and (2) can be trusted. A quick literature search told me that Huber SE will give a correct SE even if the model is specified, but the parameters estimated may be biased. Under this situation, will (1) and (2) provide biased parameter estimates?


title: Name of a tool to find out the best power to use in order to normalize a variable

text: It should be an easy question. I am looking for the name of a tool. It is used to normalize a variable with the best possible power. I think it uses an iterative process to find out this best power. What is the name of this tool? Thank you


title: Sampling from a distribution instead learning parameters by inference

text: I am not sure if this question/approach makes sense. I need to generate synthetic data based on the parameters that are learned from a model, feed the generated data back to the model and learn the parameters again (bootstrap) for calculating some p-values. Inference is very time consuming and if we want to do it thousands times, it is out of question. I thought that if I centre a distribution **D** based on the learned parameters and sample from that distribution, I can avoid doing inference and learning parameters again. Is it is statistically sound approach? Also is there any guideline or theory for choosing distribution **D**. The naive approach would be using a normal distribution with mean equal to our parameter. Of course depending on the chosen variance results would be very different.


title: What is the autocorrelation function of a time series arising from computing a moving standard deviation?

text: Say I have a time series of observations and I compute a measure of the variance of that time series as the standard deviation (SD) in a rolling window of width $w$ and that window is moved in single time steps over the series. Assume further that $w = \left \lceil{n/2}\right \rceil$, where $n$ is the number of observations, and that the window is right aligned; I have to observe $w = \left \lceil{n/2}\right \rceil$ values of the series before I start yielding moving window estimates of the SD of the time series. Is there an expected form for the ACF of the new time series of SD values? I presume the dependence on previous values will relate to the window with $w$, but is the ACF of such a series related to the ACF of an $\mathrm{MA}(w)$ process? ### Background I'm trying to think through the implications of deriving a time series of the variance of the original time series via moving windows. Having computed the derived series of SD values the next step that is commonly applied is to see if there is some trend in the derived series of SD values. As each value in the derived series depends to some extent on the previous values of the original series the values of the derived series are not independent. Thus a question that crops up frequently is how to account for that lack of independence. Such computations (the moving windows) are often done to time series to look for evidence of indicators (increasing variance, increasing AR(1) coefficient) of impending threshold response (so-called critical transitions).


title: Asymptotic Property of the Likelihood surface

text: I have a questions which I am not quite sure how to frame, so I apologies if it does not make sense, but I will try my best to make it interpretable. I have been running some network models in R, using the sna package and lnam function. By using this function I am implementing a Spatial auto-regressive disturbance model. When I run this code I get the error:               Warning in sqrt(diag(o$acvm)): NaNs produced      Where acvm is the the estimated asymptotic covariance matrix for the model parameters. This also gives 'NaN' p-values in the summary of the model. When researching a solution I came across this blog which seems to have a similar problem. The author of the blog suggests that this is because _"This is an asymptotic property of the likelihood surface, and for finite sample this approximation can be quite bad"_ I have asked for the blogger to elaborate - as you will see from the link - however have not had a reply (understandably as the link is quite old) and was hoping that someone here at Cross Validated could shed some light on the issue. Please do not hesitate to ask if anything needs clarifying , like I said I am still not sure how to frame the question, so it is likely that I have not divulged all the necessary information. Thanks, Sam


title: Problem understanding 'Percentage' when working out Confidence Interval

text: I have attempted to use on line tools to work out a population sample size for an experiment I will be running. However I do not really understand the "Percentage" part I have attempted to use this www.calculator.net Sample Size Calculator but I don't feel happy not understanding it even if most of the information I have read has said to leave it at 50 (the worst case scenario?). The description on the website does not seem to help at all > "Percentage: The percentage of a particular answer was chosen." and some other resources seem only slightly better. > "Your accuracy also depends on the percentage of your sample that picks a > particular answer..." My issue is, what answer, and what question?, what if all my survey questions are Likert types scales and the user has to answer them? Everyone who undertakes this survey has to complete it or their answers are not recorded, so is my percentage 100%? Any help is much appreciated


title: Vector Quantization of heavy tailed distribution

text: I'm generating with Monte Carlo simulation some stock price $X$.   Once I have the stock price sample, I want to cluster it with 100 points $\hat{X}$. My problem is that the error associate with my k-mean clustering $$ E \lVert X - \hat{X} \rVert^2, \tag{1} $$ doesn't decrease with the number of simulation $n$.   My cluster $\hat{X}$ is obtained with Lloyd algorithm and SciKit learn. I looked at the empirical pdf of X, obtained through Monte Carlo, and it seems to be very heavy tail.   Is this the reason why my k-mean error grows with $n$?   Do I need $n$ to be really really large? How can I quantize $X$ so that the error (1) is small?