title: Why do ANOVAs and GLM (negative binomial model) give different results for interaction effects?

text: In my study, children repeat sentences and I count the errors (DV = error rate). Children are divided into two groups (Group factor), and there are two different types of sentences (Condition factor). I am particularly interested in the Group by Condition interaction, which is theoretically important in my field (developmental psycholinguistics). When I run a by-subjects ANOVA, the "gold standard" in my field, I get a significant Group by Condition interaction. The assumptions of the ANOVA pretty much hold (normal distribution of residuals according to a QQ plot, homogeneity of variance (Levene's), Sphericity (Mauchly's)). The data themselves are not normally distributed, but I gather it's the distribution of the residuals that really matters. I'm also interested in applying GLM approaches. The raw data (i.e. both groups combined not averaged in any way) are strongly right-ward skewed and are well- approximated by a negative binomial model (I've used STATA nbvargr routine). When a negative binomial regression is run the residuals show a normalish distribution according to a QQ plot (better than the poisson and gaussian models) though not nearly as good as the ANOVA. I also get a significant interaction but it is in the OPPOSITE direction to that of the ANOVA (a negative as opposed to positive coefficient). I think this is because the ANOVA assumes an additive model whereas the negative binomial assumes a multiplicative model (as it employs logs). My question is, which analysis is "right", i.e. which best describes the underlying processes which result in the observed data? Though it has to be said that the ANOVA fit is excellent, there is a substantial loss of data as they are averaged for each cell of the ANOVA design (participant by group by condition = 50 observations per cell). From this perspective it's no wonder that the fit of the negative binomial model is not quite as good, as the model has to account for 50 times the data. In addition, the negative binomial model attempts to model the data generating process, i.e. a poisson process, whereas I don't think the same thing can be said for the ANOVA. (NB the issue of how to interpret interactions, and whether to use an additive or multiplicative model is hot issue in Psychiatry (Kendler, K. S., & Gardner, C. O. (2010). Interpretation of interactions: guide for the perplexed. The British Journal of Psychiatry, 197(3), 170–171. doi:10.1192/bjp.bp.110.081331)) Here's some code showing how the interaction effect is negative for the negative binomial model, but positive for the ANOVA.               ** ld = error rate, DV     ** limp = group (0 = control group, 1 = group of interest)     ** hard = experimental condition (0 = "easy" condition, 1 = "hard" condition)     ** morphemes = length of sentence in morphemes, reqd as an offset for     ***negative binomial          . nbreg ld limp##hard, robust cluster(id) offset(morphemes)          Fitting Poisson model:      [model-fitting procedure omitted for simplicity's sake]               Negative binomial regression                      Number of obs   =       3397     Dispersion           = mean                       Wald chi2(3)    =     116.39     Log pseudolikelihood = -7574.8427                 Prob > chi2     =     0.0000                                          (Std. Err. adjusted for 34 clusters in id_num)     ------------------------------------------------------------------------------                  |               Robust      ld_morph_br |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]     -------------+----------------------------------------------------------------           1.limp |   2.235774   .2752258     8.12   0.000     1.696341    2.775206           1.hard |   .7013264   .1065947     6.58   0.000     .4924047    .9102482                  |        limp#hard |             1 1  |  -.3988122    .123109    -3.24   0.001    -.6401013   -.1575231                  |            _cons |  -10.25841   .2499824   -41.04   0.000    -10.74837   -9.768454        morphemes |   (offset)     -------------+----------------------------------------------------------------         /lnalpha |    .696452   .0583161                      .5821544    .8107495     -------------+----------------------------------------------------------------            alpha |   2.006621   .1170184                       1.78989    2.249593     ------------------------------------------------------------------------------      Now dataset is collapsed to obtain the mean error rate per participant for each cell of the ANOVA design               . collapse (mean) ld_morph_br, by(id_num limp hard)          .      . anova ld_morph_br limp / id|limp hard limp#hard, repeated(hard)                                     Number of obs =      68     R-squared     =  0.9819                                Root MSE      = .323407     Adj R-squared =  0.9621                            Source |  Partial SS    df       MS           F     Prob > F                  ------------+----------------------------------------------------                        Model |  181.662443    35  5.19035552      49.62     0.0000                              |                         limp |  126.729506     1  126.729506      83.17     0.0000                  id_num|limp |  48.7592906    32  1.52372783                     ------------+----------------------------------------------------                         hard |  5.62123672     1  5.62123672      53.74     0.0000                    limp#hard |  .552409685     1  .552409685       5.28     0.0282                              |                     Residual |   3.3469382    32  .104591819                     ------------+----------------------------------------------------                        Total |  185.009381    67  2.76133405                  Between-subjects error term:  id_num|limp                          Levels:  34        (32 df)          Lowest b.s.e. variable:  id_num          Covariance pooled over:  limp      (for repeated variable)          Repeated variable: hard                                               Huynh-Feldt epsilon        =  1.0323                                               *Huynh-Feldt epsilon reset to 1.0000                                               Greenhouse-Geisser epsilon =  1.0000                                               Box's conservative epsilon =  1.0000                                                      ------------ Prob > F ------------                       Source |     df      F    Regular    H-F      G-G      Box                  ------------+----------------------------------------------------                         hard |      1    53.74   0.0000   0.0000   0.0000   0.0000                    limp#hard |      1     5.28   0.0282   0.0282   0.0282   0.0282                     Residual |     32                  -----------------------------------------------------------------      Command "regress" used to obtain coefficients               . regress                Source |       SS       df       MS              Number of obs =      68     -------------+------------------------------           F( 35,    32) =   49.62            Model |  181.662443    35  5.19035552           Prob > F      =  0.0000         Residual |   3.3469382    32  .104591819           R-squared     =  0.9819     -------------+------------------------------           Adj R-squared =  0.9621            Total |  185.009381    67  2.76133405           Root MSE      =  .32341          ------------------------------------------------------------------------------      ld_morph_br |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]     -------------+----------------------------------------------------------------           1.limp |   3.139737   .3327826     9.43   0.000     2.461881    3.817593                  |      id_num#limp |                  | [omitted for simplicity's sake]                  |           1.hard |   .3947685   .1109275     3.56   0.001     .1688164    .6207205                  |        limp#hard |             1 1  |   .3605257   .1568752     2.30   0.028     .0409813    .6800701                  |            _cons |   .9226158   .2353129     3.92   0.000     .4432992    1.401932     ------------------------------------------------------------------------------      Now we have a positive coefficient


title: Ordinal dependent variable with continuous independent variables

text: I have an ordinal dependent variable, named D, which varies from very small, small, medium, big, to very big. This variable depends on the independent variables X, V, which are continuous variables. Could I perform an ordinal logistic regression analysis in the former data in order to obtain the relationship between D and X, V? Could I obtain in this way the D variable as percentages?


title: Predicting the behavior of a model based on another

text: I'm trying to find a way of predicting various models from one "perfect model", using EXCEL. i.e. If I assume that all models should behave like the original one, for which I have all the data, I would like to find out the missing data based on the "perfect model". So the problem is as follows: I'm selling the same TV model to 10 different clients and over the span of 1 year I record how many complaints each client made. So the data looks like this:   Client ID-#of Complaints   Client1- 3   Client2- 17   Client3- 0   Client4- 7   Client5- 6   Client6- 20   Client7- 1   Client8- 7   Client9-2   Client10- 5   Now If a have the data of 10 other clients for half a year, can I predict what the # of complaints will be by the end of the year?   Client ID-#of Complaints   Client1-13   Client2- 7   Client3-10   Client4- 3   Client5- 4   Client6- 2   Client7- 9   Client8- 7   Client9-1   Client10- 2   Any advice on how I should go about this using excel would be greatly appreciated   Thanks


title: Using decision trees to make a binary decision

text: I have a button that I can press or not press, a binary target that I would like to be 1 as often as possible, and a bunch of features. I also have a bunch of (feature, button choice, target) data, in which the decision to press the button was made randomly (and thus independently of the other features). I'd like to train something to decide, given a feature vector, whether or not it's a good idea to press the button. The obvious thing that occurs to me is to make IsButton a feature, train a decision tree/forest. Then, given a feature vector, press the button if and only if the score of the decision tree with IsButton = 1 is greater than the score of the decision tree with IsButton = 0. This produces somewhat sensible results, but it's also very noisy, and I feel like there should be a better way of deciding whether or not to press the button. The binary target is 1 about half of the time, and in aggregate, pressing the button only moves the average value of the target by about 1%. So I'm subtracting two mid-sized numbers to try to estimate a very small difference. IsButton ends up being very low in feature importance, which is fine - IsButton probably isn't the most important factor in terms of influencing the target, I just want to figure out the small effect of pressing the button. But since much of my tree doesn't involve the IsButton feature, I'm spending most of my time building branches, while important for estimating the probability of the target, turn out to have no importance for the actual decision as to whether to press the button. Is there a better training method that reflects the fact that in the end I don't care about the absolute value of the target, but only the effect of pushing the button on the target?


title: Bias of autocorrelation function in finite sample for ARMA processes

text: I remember reading that, when estimating the autocorrelation function of a univariate ARMA time series using finite samples, the estimate is biased and specifically the lag-1 ACF is negatively biased in the case of iid observations. However, I can't recall the source and I am therefore not sure this statement is true. Does anyone know whether this is true or not? Thanks.


title: Two different formulas for AICc

text: Wikipedia's page on AIC gives a formula for the AICc, a "corrected" version of the AIC that helps to avoid overfitting when the sample size is small relative to the number of parameters in the models being considered. Wikipedia's formula is consistent with Burnham & Anderson (2002), with the correction term $2K(K+1)/(n-K-1)$, where $K$ is the number of parameters and $n$ is the sample size. However, in every source I have found prior to Burnham & Anderson (2002), the correction term is stated as $2(K+1)(K+2)/(n-K-2)$. For example, see Anderson et. al. (1994). Why are the two formulas different? Was there a mistake in the original derivation that was later corrected, or did the set of assumptions change at some point? I have not been able to find any explicit mention of the difference between these two formulas.


title: Standard Error of MLE

text: I have a model and a vector of parameters that I estimate using ML. Since I'm doing a MonteCarlo I can simulate as many times as I want, i.e. data is not real data. I was just wondering how I can compute the Standard Errors of my estimator and what their significance is.


title: Calculate error for monthly exchange rates

text: I'm working on a project that will convert $USD to other currencies. The system will use average monthly exchange rates instead of daily because of pre-existing limitations. I need to calculate an estimate of the amount of error that will result from using monthly average. My idea is to   1. Download historical daily exchange rates for several years   2. Calculate the monthly average for each month   3. Subtract the monthly average from each daily amount that month   4. Find the standard deviation for the differences I would report the standard deviation as the expected error amount. Is this a valid approach?


title: How true is this slide on deep learning?

text: I was listening to a talk and saw this slide: ![enter image description here](http://i.stack.imgur.com/qg54F.jpg) How true is it?


title: using cluster information in multiple imputation

text: i need to impute a dataset all categorical variables before doing analysis. I can just simply **impute** with mode of all data or a variable. I belief that better option will be to **classify** the subjects (observations) using some short of **clustering algorithm** and then use this information to **impute** the data. The following is small data (though the real data is really big) and the idea.               md <- data.frame(V1 = c("AA", "AA", "AA", NA, "AB"), V2 = c("AB", "AB", "BB", "BB", "BB"), V3 = c("BB", NA, "BB", "BB", "BB"), V4 = c("AA", "AA", "AA", "AA", "AA"),       V5=c(rep("AB", 5)), V6 = c("BB", "BB", "AB", "AB", NA), V7 = c("AB", "AB", "BB", "BB", "BB"))           md         V1 V2   V3 V4 V5   V6   V7     1   AA AB   BB AA   AB   BB AB     2   AA AB <NA> AA   AB   BB AB     3   AA BB   BB AA   AB   AB BB     4 <NA> BB   BB AA   AB   AB BB     5   AB BB   BB AA   AB <NA> BB      In visual observation, the sample 1 and 2 are more similar in one cluster while 3,4,5 are in second cluster. The side is hand drawn dendogram (I could not due HC because of missing values). Now I would like to impute all missing values based on the similarity. For example, missing value in column 1 is most likely to be AA as sample 4 is more similar to 3 than 1, 2, or 5. Similarly the 5 column missing value is BB as the **neighbor** in the cluster is also BB. Similarly column 6 value should be AB as the closest similar has AB. and so on. ![enter image description here](http://i.stack.imgur.com/DIcPA.jpg) Thus completed data would look like: ![enter image description here](http://i.stack.imgur.com/xa5j3.png) How can we perform this ?


title: Are there any well-understood circumstances when we should expect cross-validation to be indifferent to $\alpha$ in elastic net tuning?

text: I am fitting a lot of elastic net models, simultaneously tuning $\lambda$ and $\alpha$. I am often coming to the following conclusions:   * Cross-validation error is much more sensitive to changes in $\lambda$ than $\alpha$   * In fact, the differences in cross-validation error are negligible between $L_1$, $L_2$, and $\alpha$ = 0.5 Without looking at my data or knowing anything about it, even the family of GLM I'm using, are there any generalizable results about the conditions when this would be the case? Thank you.


title: Plotting a boxplot against multiple factors in R with ggplot2

text: I have a function that plots (I think) various conditions within a set of motor racing timing data using boxplot:               library(RCurl)     gsqAPI = function(key,query,gid=0){ return( read.csv( paste( sep="",'http://spreadsheets.google.com/tq?', 'tqx=out:csv','&tq=', curlEscape(query), '&key=', key, '&gid=', curlEscape(gid) ) ) ) }     pd=gsqAPI('0AmbQbL4Lrd61dHVNemlLLWNaZ1NzX3JhaS1DYURTZVE','select A,C,E,G',gid='6')          boxplot(Time~Stint*Session*DriverNum,data=subset(pd, (DriverNum==1 | DriverNum==2)&Time<110),las=2,ylab="Time(s)",ylim=c(80,110),xlab="FP Stint:Session:DriverNum",main='F1 2011 Italy FP Times (Red Bull)')      What's the equivalent way of doing this using ggplot? (That is, what's the ggplot equivalent of the boxplot command above?) As a second part to the question, how would I generate a separate chart using ggplot that displays the boxplot overlaid with a plot of the raw data points within each condition?


title: RBM - clamping the visible units during reconstruction

text: I am currently attempting to code an implementation of Deep NN for image recognition using stacked RBMs (by adding a linear classifier at the "top" of the last RBM in the stack). I am at the very early stage of learning RBMs. I have read several papers and some powerpoint presentation from the Toronto University Machine Learning group. From this material, I understood when computing (vi*hj)0 - (vi*hj)n (to update weights) I would first set the visible units to a data vector and use that to compute the states of the hidden layer (this would be used for (vi*hj)0) then use the newly computed hidden units to create a reconstruction of the visible units and use that reconstruction in a loop to calculate (vi*hj)n by updating the hidden layers with the reconstruction, then updating the reconstruction, and continuing n times updating the hidden layer with the newest reconstruction of the visible layer. Recently I came by this : http://halfbakedmaker.org/?p=11748 Here the author clamps the training vector to the visible units not only to calculate the first iteration of updates, but through all the process. His loop is also short, which is fine I guess, I read about such a shortcut in one of Hinton´s powerpoint files. My questions are, if there is a correct way of doing this, which one is it? Clamping the visible units for the first iteration only and then using the newest reconstruction or clamping the visible unit through the whole process? If both ways are correct, are there performance (learning) differences? Thanks.


title: How can I weight the independent variables in a regression?

text: If I have a regressive equation of the form: $Y = aX_1 + bX_2$, $X_1$ and $X_2$ are predictor variables of size (n x 1). I have an array of weights that I would like to apply to $X_1$ and $X_2$, named $w_1$ and $w_2$ also of size (n x 1) that I would like to apply here, such that the estimates of $a$ and $b$ are based on these weights. $w_1 + w_2 = 1$ (this is based on my own constraints), and essentially describes the amount of influence that $X_1$ and $X_2$ have on $Y$ at any given point in time. For instance, if $w_1 = 0.9$ and $w_2 = 0.1$, I want to weigh the corresponding values of $X_1$ and $X_2$ accordingly. I was reading up on weighted least squares, but that is a technique used to weight the errors of the fit, rather than the predictor variables themselves. **How can I weigh my variables?**


title: Type of Sampling

text: I am conducting a research project, and have recruited participants through the organisation that has come on board with the project. What kind of sampling is this? At first I thought it was convenience sampling because the organisation is participating, but now I am not so sure. Participants will be recruited via email (through the organisation's email base) and invited to participate. If anyone could help me out would be much appreciated.


title: Issues with using Expectation Maximization algorithm

text: I was using the EM algorithm to maximize a partially observed likelihood. However, I have certain doubts. Normally, the algorithm works fine. I could print the value of the log likelihood of the observed data that I am trying to minimize and the values keeps on going down and converge. However, at times, I could see that it goes down for a while and then starts increasing. What could be the reason for that? As far as the code is concerned, I am pretty sure it's fine. Any suggestions?


title: will a correlation between a percentage variable and a numerical variable work?

text: I am trying to conduct a simple correlation analysis between a group of variable, however I have a number of these variables in percentage change (or percentage yield) form and the others in strictly numerical prices. I was wondering whether it is possible to find the correlation coefficients between a percentage variable and a numerical variable? I have tried with the data as it is as well as with all variables in some percentage return format with varying results, hence my current confusion! The data I am working with relates to the prices of stocks, gold etc and the percentage yield/change of bonds and CPI %change. Any help would be greatly appreciated!


title: How to do mean comparison with trending data?

text: A client with a website has recently released a tablet app. The theory is that after installing the app a customer is less likely to log on to the website as they will use the app instead. For each customer that installed the app i have gathered data that shows the number of times they logged-in to the website in the two weeks before they installed the app, and the number of times they logged-in in the two weeks after app installation. The average for the after case is indeed less, and i've run a t-test in R which shows the difference in means is statistically significant. However, the data is not normal by a long shot; and wikipedia pointed to a wilcox test which i've also run in R and which also says the difference is significant. first question: is the (paired) wilcoxon signed rank test an appropriate test with non-normal data? But, if the analysis done so far is ok (please point out if it is not; i did a stats course 15 years ago and have spent the last 2 days re-learning the concepts), i have a problem. The interval of interest (mid-june to mid-july) includes the financial year rollover. So the trend of logons across all our customers is generally down over the interval. More people tend to logon just before the end of financial year to do what-not. There are hiccups in the trend, but it is is for sure that more customers logged on in last half of june than first half of july. So i'm worried that the statistically significant lower average for the after- case may just be a manifestation of this trend, and not significant at all. How do i take the down trend into account in the analysis?


title: One way Anova Test validation

text: Let me give you an example to show what i am trying to. Let's say 50 patients started a cure for a disease and it took 2 months. So, I want to compare the first month's tumor sizes to the second month's tumor sizes for each of 50 patients. And i should be able to tell the patients tumor sizes are stable/increase/decrease over the cure period. I have 50 patients' data that describe the effect of treatment in the first month and the second month. For example, the data looks like               patient-id  first-month  The size of tumor volume     1           1                 90     1           1                 89     1           1                 88     2           1                 90     2           1                 90     2           1                 91     1           2                 75     1           2                 80     1           2                 83     2           2                 98     2           2                 95     2           2                 97     ..          ..                 ..     ..          ..                 ..     50          2                 67      Here if the first month is 1, that means it is in the first month, 2 means it is in the second month of the treatment. And the IDs are just number to represents each patient. There are more than 3 measurements for each patient, but i have shown only 3 of them for brevity. How to compare the first month of the treatment to the second month of the treatment for the disease. I used A one way Anova to compare the first month to the second month by using the "The size of tumor volume" as a dependent variable and "first-month" as a factor in SPSS. The results show there is a significant difference between the first month of the tumor size and second month of the cancer size for all patients. But, i did not include patient-ids in SPSS, should I use it? Should I compare each patient's first month data and second month of the treatment period separately? So, instead of having one Anova test result, should i have 50 (number of patients) Anova results to say that the differences is significant?


title: Polya urn waiting time distributions

text: From an urn containing 1 ball each of n different colours, a ball is drawn, its colour noted and then it is returned to the urn along with k balls of the same colour. What is the distribution of the number of balls of each colour in a multicolour urn (after a certain number of draws)? Additionally, what is the distribution of wait times before drawing the first ball of a given colour in the multicolour urn model? EDIT in simulations i've run, i've found that the number of balls of each colour is described by the negative binomial, and the time to the first ball for each colour given by the Lomax (pareto type II) distribution. Ideally I'd like to find a proof or a reference to confirm these results.


title: Compute expectation of a random variable given the density function

text: The probability density function is given by: $f(t)= (8-t)^2 /9$ for $5 \le t \le 8$ Compute the mean daily CPU time. Hence state the mean of a new variable $W=T+12$ hours. For this question, based on the information given, I would think it to be a continuous uniform distribution. I know that the mean is $b+a$. Also, for the other part, I'm wondering how to make $T$ the subject of the formula so that $T=W-12$ hours and then put that in wherever $T$ is? Side note, can anyone tell me how to insert actual symbols into the posts (eg. the actual 'beta' symbol and not 'b')?


title: R: estVar and cor

text: I have a problem when estimating a multivariate Regression Model. I know that `cor(y1, y2)` is strongly positive but when I am trying to estimate a regression model like               x <- lm(cbind(y1, y2) ~ X)      then `estvar(x)` returns a slightly negative covariance. Is that possible and why? I can´t explain it to myself.


title: Computing sample size for having at least one individual of each level for a categorical variable?

text: Imagine a categorical variable with 10 levels. I would like to know how many individual I need to sample from the whole population to have a .95 or .99 probability of having at least one individual of each category for that variable?


title: percentage of variance

text: I try to implement PCA in MATLAB using the Fisher iris data. I am having a problem finding out the percentage of variance in Principal Component Analysis. I get the following result:               pc =              0.3614    0.6566   -0.5820    0.3155        -0.0845    0.7302    0.5979   -0.3197         0.8567   -0.1734    0.0762   -0.4798         0.3583   -0.0755    0.5458    0.7537               latent =              4.2282         0.2427         0.0782         0.0238               ans =              0.9246         0.9777         0.9948         1.0000      Now how do you choose the percentage of variance? (In the "hald" dataset it is 98%, but I also don't understand why 98%.)


title: Methods to determine cut-off between groups

text: I am working with a dataset of diseased vs non-diseased patients. I'm looking for a cut-off value for determining the status of a new patient. This is intended for lab technicians and MDs to allow for quick diagnosis. It's got to be practical, so please spare me from discussing whether it's the right thing to do—I already know it's not a favorable approach. To be as thorough as possible I want to use multiple methods to determine the cut-off. My hope is that these methods will yield similar values. So far I've used **logistic regression (ROC)** and **mixture modeling**. So my question: what other techniques can be used to determine the cut-off? 1/10/2014 EDIT: To be clear - I actually have 3 variables to be classified in the dataset, however only one of them is practically useful. The other two are highly correlated with one another, and not really correlated with the outcome. If I have to use more than one, I will use a ratio, so no need for multivariate classifiers. The goal is to be simple - just one cut-off value.


title: Likert scales: analyzing scales or items?

text: I am writing an essay based on a survey I made. First I set my hypothesis = Who is an impulsive buyer tend to be loyal to shop X. Than using "marketing scales book" I set different scales for each element:   * Impulsive buying tendency: I tend to buy things without thinking, I buy on impulse.   * Store loyalty: I love spending my time at store X, I am a loyal customer of store X, I feel in touch with store X. Than I asked how they agreed 1-7 for each one. Now I have to test my hypothesis… but witch is the best way? Is there a method to test all 2 scales representing IMPULSIVENESS VS 3 scales representing LOYALTY one shot, as a group so it will be fully representative or I have to test every single scale of impulsiveness vs. single scale of loyalty and than, I made a mean?


title: Using the covariance matrix to calculate correlations

text: I'm a bit lost here. I have a matrix of response variables, $y$, and I fit a model to account for a number of predictor variables, say $x_1$, $x_2$ and $x_3$:               lm.m <- lm( y ~ x1 + x2 + x3 )      I can then extract the estimated covariance matrix of the residuals (if I understand correctly), for example using the `estVar` function from R:               c <- estVar( lm.m )      What should I do to turn this into a correlation matrix? Divide each cell by the product of square roots variances of the residuals? Or the original variables? I'm confused. What I'm trying to achieve is to calculate the correlation matrix of the residuals from the model. Something that could be achieved by               y2 <- apply( y, 2, function( yy ) lm( yy ~ x1 + x2 + x3 )$residuals )     c <- cor( y2 )


title: Back-transformed confidence intervals

text: Having come across this discussion I'm raising the question on the back- transformed confidence intervals conventions. According to this article the nominal coverage back-transformed CI for the mean of a log-normal random variable is: $\ UCL(X)= \exp\left(Y+\frac{\text{var}(Y)}{2}+z\sqrt{\frac{\text{var}(Y)}{n}+\frac{\text{var}(Y)^2}{2(n-1)}}\right)$ $\ LCL(X)= \exp\left(Y+\frac{\text{var}(Y)}{2}-z\sqrt{\frac{\text{var}(Y)}{n}+\frac{\text{var}(Y)^2}{2(n-1)}}\right)$ /and not the naive $\exp((Y)+z\sqrt{\text{var}(Y)})$/ Now, what are such CIs for the following transformations:   1. $\sqrt{x}$ and $x^{1/3}$   2. $\text{arcsin}(\sqrt{x})$   3. $\log(\frac{x}{1-x})$   4. $1/x$ How about the tolerance interval for the random variable itself (I mean a single sample value randomly drawn from the population)? Is there the same issue with the back-transformed intervals, or will they have the nominal coverage?


title: Ongoing chance of drawing a particular number from a bag of numbers

text: I'm sure this question has been asked many times over applied in different ways so perhaps there is a duplicate question that can simply be pointed to. Let's say there is a bag with 10 numbers in it (1 through 10). Assume that after each time a number is drawn, that number is put back into the bag so there are always 10 numbers to choose from. The chance of me drawing a specific number is obviously always going to be 10% for any given attempt. But how would I calculate the chance of drawing a particular number for x number of attempts in this scenario? EDIT: Just to make sure I'm being completely clear here let me retell the scenario in a slightly different way. Suppose there are x number of bags and each one contains 10 numbers (1 through 10). What is my chance of drawing a 5 from one of the bags if I draw from each bag exactly once?


title: Aggregate stock return predictive panel analysis

text: I am trying to run a panel OLS predictive regression on stock returns for 7 countries using 8 macro economic variables. I have already done the simple linear regression analysis that has highlighted the presence or absence of predictive power for each variable per country, and i am aware of the various problems of sample bias and data mining, which can be addressed with a robust bootstrap procedure. However, I have been advised to run a panel analysis to look for a common break date in the data accross countries, and i would appreciate advise on how to go about this. I have run the pooled OLS regression and i went further to run the Fixed effects regression, which has increased the standard errors a bit and both tests confirm the individual predictive ability of the same variables (interest rate, tbills and business confidence index) as the earlier simple linear regression. My main problem is how to go about identifying these breaks both in the simple OLS and Panel OLS. My understanding is that if there are breaks then an out of sample technique would allow the parameters in the predictive regression to vary over time. I am not so great statistically, but any pointers in the right direction would be highly appreciated


title: How can I find distribution from mean and variance

text: If I know the mean and variance of a discrete random variable $\in \\{0,1,\dots \\}$, How can I find the distribution?


title: How to calculate relative error when the true value is zero?

text: How do I calculate relative error when the true value is zero? Say I have $x_{true} = 0$ and $x_{test}$. If I define relative error as: $\text{relative error} = \frac{x_{true}-x_{test}}{x_{true}}$ Then the relative error is always undefined. If instead I use the definition: $\text{relative error} = \frac{x_{true}-x_{test}}{x_{test}}$ Then the relative error is always 100%. Both methods seem useless. Is there another alternative?


title: Computing the expected value using the sample mean gives poor convergence?

text: This is for Project Euler problem #371: http://projecteuler.net/problem=371 My solution was to generate a huge random pool of numbers between 0 and 999 inclusive, repeatedly ask the question 'How many numbers do I need to pull into a holding area before I get a pair that adds to 1000?', and then average the resulting list(to get the expected number of plates). The problem is that hand-optimizing my program lets me crunch about 10 million samples of the second step(so a total pool of 400 million random numbers are used), but I'm not even getting 2 digits of accuracy behind the decimal(I'm using exact integer arithmetic):   * 39.6638047   * 39.6791329 There could be an error in my program otherwise, but I'm not asking you to debug my program. I'm sure there's a better way to do this in principle than brute-forcing it with the Law of Large Numbers, so I wanted to ask:   * Should I be expecting more than 2 digits of precision with my approach?   * If someone handed you a huge number of random samples, how would you efficiently compute their expected value to a high precision?


title: Bootstrap intervals from a sample not obtained with 'boot'

text: I have a bootstrap sample of a parameter and I want to calculate the confidence intervals that 'boot.ci' from the R package 'boot' produces: "norm","basic","perc","bca"and "stud; but I do not have an object of class 'boot', I just have the bootstrap sample. Do you know another package that I could use for this purpose? Thanks in advance.


title: Linear combination of two dependent multivariate normal random variables

text: Suppose we have two vectors of random variables, both are normal, i.e., $X \sim N(\mu_X, \Sigma_X)$ and $Y \sim N(\mu_Y, \Sigma_Y)$. We are interested in the distribution of their linear combination $Z = A X + B Y + C$, where $A$ and $B$ are matrices, $C$ is a vector. If $X$ and $Y$ are independent, $Z \sim N(A \mu_X + B \mu_Y + C, A \Sigma_X A^T + B \Sigma_Y B^T)$. The question is in the dependent case, assuming that we know the correlation of any pair $(X_i, Y_i)$. Thank you. Best wishes, Ivan


title: Statistical Practices Using Sparse Data: Methods for Approximating Standard Deviation

text: Suppose I know that for a discrete, non-negative r.v. $X$ that $X | X \geq 1$ has $\mu = 3.3$ while when $X \geq 0$ has $\mu = 2.1$. That is, the subset of the population that already has a value $1$ or more has a proportionally higher mean. I would like to assume that $X$ is more likely to evaluate to some value $x \leq 5$ so that the distributions of $X$ and $X | X \geq\ 1$ could roughly be modeled by a gamma distribution. Is there a way that I can estimate what the standard deviation for both distributions? What is the most information that can be gained from these small pieces?


title: Classification accuracy

text: I am classifying text based on news headlines and I am achieving accuracy up to approx 80%. I want to improve it more. But issue is that when I calculate the same with synonyms using the code below:               Doc = actxserver('Word.Application')         X = cellfun(@(word) invoke( Doc,'SynonymInfo',word), words, 'UniformOutput', false);     Synonyms = cellfun(@(X) get(X,'MeaningList'), X, 'UniformOutput', false);     Synonyms = cellfun(@(X) [words{X}; Synonyms{X}], num2cell(1:numel(words)), 'UniformOutput', false);      My accuracy falls very badly below and reach 40% or less. Why is this happening?


title: Truncated SVD matrix reconstruction: what is the meaning of the real values?

text: Im my algorithm, I am working with Singular Value Decomposition (SVD). I have an input matrix $A_{in} \in \\{0,1\\}^{(m * n)} $, made by $n$ rows and $m$ colums. All the entries are 0 or 1. I decompose it in $A = U * \Sigma * V^{T}$ After choosing a proper truncation level $k$, I construct an output matrix $A_{out} \in \mathbb{R}$, this way: $ U_{k} * \Sigma_{k} * \; V^{T}_{k} = A_{out}$ All the $A_{out}$ entries are real valued. **What is the meaning of each real value?** Is it the likelihood that the element is near to 1.0 in the reconstruction? Or are the real values just to create an order, from the less likely to the most likely?


title: How to interpret the coefficients and plot the relationship between IV_i and DV?

text: I'm using Tobit regression and I'm not sure how to interpret the coefficients. I made sample data with random function in Excel.               DV: Y( the value ranges from 0 to 100 at fist but is censored at 20 and 80)     IV_1: X1(the value ranges from 0 to 5 in interger value)     IV_2: X2(the vales is Y*rand() in excel)      As explained DV is censored as follows: ![enter image description here](http://i.stack.imgur.com/vOZp6.jpg) and I conducted Tobit regression in R. Here's the result:                               Estimate Std. Error   z value        pvals     (Intercept):1 12.0819679 2.27337253  5.314557 1.319982e-07     (Intercept):2  3.1529486 0.04336364 72.709498 0.000000e+00     X1             7.5977990 0.65988825 11.513766 6.956904e-29     X2             0.8100315 0.06801476 11.909643 1.149151e-30          Efron's R-square = 0.518107      With This, I want visualize how independent variable, namely X1, affects Y according to X1's value.(note that it ranges from 0 ~ 5 in integer) Thus I plotted like this: ![enter image description here](http://i.stack.imgur.com/zRqkY.jpg) It is draw from the coefficient of X1(7.5977990) and multiplication of each value of the X1 and the coefficient(0.00, 7.60, 15.20, 22.79, 30.39, 37.99). I guess this is valid approach for OLS because each coefficient of IV directly affects Y. However, in non-linear case like Tobit regression, I'm not sure it is okay. If there is any advice or other suggestions, please let me know. Thank you!


title: Why is the Pearson correlation 1 when only two data values are available?

text: I am trying to obtain a Pearson correlation between 6 different variables (represented by columns in the matrix below) with two datapoints each (rows). This is the matrix:                    scer       bay      par       mik      glab       lac     var1 2.2273444 2.0923416 2.044007 1.7664921 1.3832924 2.4294228     var2 0.3000878 0.2792936 0.286928 0.3246768 0.4946222 0.3083171       When I apply the standard R code for correlation:               cor(mat)      I obtain the following result:                    scer bay par mik glab lac     scer    1   1   1   1    1   1     bay     1   1   1   1    1   1     par     1   1   1   1    1   1     mik     1   1   1   1    1   1     glab    1   1   1   1    1   1     lac     1   1   1   1    1   1      If I add another two rows to the original matrix:                               scer       bay       par       mik      glab       lac     var1    2.2273444 2.0923416 2.0440068 1.7664921 1.3832924 2.4294228     var2    0.3000878 0.2792936 0.2869280 0.3246768 0.4946222 0.3083171     var3    1.1399738 1.2899311 1.1071462 1.0180361 1.4507592 2.4078977     var4    0.7107440 0.6415944 0.7197905 0.7357125 0.4571745 0.3173547      and re-execute the above code with the new matrix, I obtain a more familiar result:                         scer       bay       par       mik      glab       lac     scer 1.0000000 0.9895959 0.9991065 0.9967358 0.7860344 0.8246286     bay  0.9895959 1.0000000 0.9916464 0.9890492 0.8647974 0.8958393     par  0.9991065 0.9916464 1.0000000 0.9991332 0.7928330 0.8310776     mik  0.9967358 0.9890492 0.9991332 1.0000000 0.7845007 0.8235245     glab 0.7860344 0.8647974 0.7928330 0.7845007 1.0000000 0.9978420     lac  0.8246286 0.8958393 0.8310776 0.8235245 0.9978420 1.0000000      Why does the correlation function return a matrix of 1s if I use 2 values?


title: Calculating the size of Type 1 error, Type 2 error and power of the test

text: Let $X$ have a binomial distribution with parameter $n=5$ and $P\in [p:p=\frac{1}{4},\frac{1}{2}]$. The null hypothesis $H_{0}:P=\frac{1}{4}$ is rejected, and The alternative hypothesis $H_a:P=\frac{1}{2}$ is accepted. If the observed value of $X_1$, a random sample of size one, is less than or equal to $3$. Find the size of `Type 1 error`, `Type 2 error` and `power of the test`. I have no idea to solve the question. I only know that ### Size of a Type 1 error = Pr[rejecting$H_0|H_0 $is true] ### Size of a Type 2 error = Pr[not rejecting$H_0|H_0 $is False] The sign $|$ denotes "given that".


title: How do programs like BUGS/JAGS automatically determine conditional distributions for Gibbs sampling?

text: Seems like full conditionals are often quite difficult to derive, yet programs like JAGS and BUGS derive them automatically. Can someone explain how they algorithmically generate full conditionals for any arbitrary model specification?


title: Weight matrix definition in weighted least square method

text: On page 220 of Kutner, the author defined the weight matrix in weighted least square method as $\sigma_i^2 = \sigma^2 / w_i$, which is problematic because $\sigma^2$ is unknown. Even more confusing is that the estimate $MSE_w$ involves $w_i$ which is in turn defined by $\sigma^2$, forming a circular argument. Am I missing something here? ![enter image description here](http://i.stack.imgur.com/DAJuv.png)


title: Relating a sum to the standard deviation

text: I have a summation that looks like $$\sum_{j=0}^{n-1} S_j(v_{j} - \bar{v})$$ The $S_j$ terms all have values in $[0,1]$. $v_j$ is a data value and $\bar{v}$ is the mean of all $v_j$. This expression is clearly related to the standard deviation. If $\sigma=0$ than all $v_j=\bar{v}$ so the result is $0$. But what I would like to do is relate the change in value of this expression to $\sigma$ for all values of $\sigma$. Is this possible? I have tried without luck but think that perhaps there is some property of $\sigma$ which I do not know which might help here?


title: 2 sided vs 1 sided test

text: I am asking this questions straight from "Data Analysis and Statistical Inference " course by Dr. Mine Çetinkaya-Rundel on coursera! Please right click on the image and one in a new tab to make it bigger.![enter image description here](http://i.stack.imgur.com/agIA3.jpg) In below example, can we do one sided hypothesis test? i want to find out if Streptomycin is effective or not. Based upon histogram, only 1 out of of 100 observations has chance of having >0.2 probability. 0.2 = (51/55)-(38/52). For some reason the professor is doing two sided test. I want to perform one sided test as i just dont want to say that both treatments are different, but i also want to prove that Streptomycin is better.


title: Which error is displayed in an error ellipse?

text: I have some bivariate data and I have calculated the error ellipse in the following way: I have first calculated the covariance matrix and then to obtain the radii of the ellipse I have taken the eigenvalues and then their squared root. Does the length of each of these radii correspond to one standard deviation in the direction of the semi-principal axes?


title: SVM basic theory?

text: I have some questions about SVM:   1. In SVM there is a nonlinear and linear SVM. What is the difference between them?   2. To do classification in SVM, we will find the linearly separable boundary (hyperplane). And if we cannot find it, we must project our input space into a features space. What does that feature space mean?   I read on the internet that some explain that the features space is a scalar, and we get that scalar by inner product. For example, I have an input space matrix $A=[\\{(1,2); 0\\},\\{(0,6);1\\}]$; where 1, 2, 0, 6 are the features and 0, 1 are the labels. How can I project my input space A into feature space?   3. Why do most people say that with kernel tricks, we actually don't project our data into a higher space? But we project it into a low space (because of that kernel tricks = scalar)?   4. There are many kernel types: linear, polynomial, RBF and sigmoid. How can we determine which kernel type we have to use for our data?   5. On this webpage, we can see that with polynomial and RBF kernel, we don't find the linearly separable boundary / hyperplane, but there is a curve, a group (for RBF). But as mentioned before, SVM finds the linear hyperplane? Why are they different? I appreciate all answers.


title: Ways the likelihood ratio test may fail

text: I have a question related to: Why is a likelihood-ratio test distributed chi- squared? On point 2 of @StasK's answer, he states: > The theorem assumes that all the relevant derivatives are non-zero. This can > be challenged with some nonlinear problems and/or parameterizations, and/or > situations when a parameter is not identified under the null.   * What derivatives are these?    * In the case of mixture distributions, when a parameter is not identified under the null, is it possible for the Information matrix to not be invertible? How?


title: How to compare two non-normally distributed samples with very different sizes? (Mann-Whitney vs Randomization/Bootstrap)

text: Perhaps this is a very basic question, but I didn't find yet a simple solution for this simple problem: I want to compare two samples (say X and Y) for a continuous variable which is **non-normally** distributed and test if X and Y are significantly different. The sample size of **X** is **N=81** and **Y** is **N=5110** , so they are quite unbalanced. My first attempt was to use the Mann-Whitney (i.e. Wilcoxon Signed Rank test). However, I am bothered with this huge difference in sample sizes. I thought that some kind of randomization or bootstrap method is a good alternative, but I am not sure if my approach makes sense. My idea was to get 1000 random samples of size 81 from Y and X and then use the Mann-Whitney to compare both distributions. The empirical _p-value_ would be the proportion of tests with _p-value_ < 0.05. I "R", I've implemented as follows:               X = data1 # sample size 81     Y = data2 # sample size 5510     R = 1000     alpha = numeric(R)          for(i in 1:R) {         group1 = sample(X, replace=TRUE)         group2 = sample(Y, size=81, replace=TRUE)         alpha[i] = wilcox.test(group1, group2)$p.value     }      Empirical _p-value_ would be the proportion of _p-values_ < 0.05:               mean(alpha < 0.05)      Does this approach make sense? How can I do this hypothesis testing correctly?


title: How and why do normalization and feature scaling work?

text: I see that lots of machine learning algorithms work better with mean cancellation and covariance equalization. For example, Neural Networks tend to converge faster, and K-Means generally gives better clustering with pre- processed features. I do not see the intuition behind these pre-processing steps lead to improved performance. Can someone explain this me?


title: Whether to use robust linear regression or bootstrapping when there is heteroscedasticity?

text: I have a dataset where I need to do linear regression. Unfortunately there is a problem with heteroscedasticity. I´ve rerun the analysis using robust regression with the HC3 estimator for the variance and also done bootstrapping with the bootcov function in Hmisc for R. The results are quite close. What is generally recommended?


title: How would YOU compute IMDB movie rating?

text: I'm doing this only for learning purposes. I've no intentions of reversing the methods of IMDB. I asked myself I owned IMDB or similar website. How would I compute the movie rating? All I can think of is _Arithmetic Mean_ For a movie data provided below computation would be > (38591*10 + 27994*9 + 32732*8 + 17864*7 + 7361*6 + 2965*5 + 1562*4 + 1073*3 > + 891*2 + 3401*1) / 134434 = 8.17055953 ![alt text](http://i.stack.imgur.com/bWeRr.jpg) * * * My rating 8.17055953 which is close to IMDB's Arithmetic Mean.   1. Whats wrong with my score? Why is it not ideal (because IMDB didn't use it)?   2. If you had to compute. How would you have done it? What factors you would consider? **Note:** I'm not asking about the ranking mechanism for the Top/Bottom listing (that is shown here : http://www.imdb.com/chart/top). This question is about how IMDb turns the arithmetic mean into the weighted mean. Take "Manos: The Hands of Fate" — Unweighted mean 2.1, weighted mean 1.5. (Or "The Godfather", another example of ballot stuffing.)


title: How to analyze data from two-tiered experiment

text: I have an experiment where each subject undergoes M sessions of treatment drug & control drug. I have N subjects. How do I 'match/pair' my data together? In particular, lets say that d(n,m) is my performance metric for patient n and session m (which can be either treatment or control). I want to know the unbiased estimator for the difference between treatment and control drugs, or r=d(:,m={treatment}) - d(:,m={control}). How can I get this and the associated error bounds. Thanks!


title: How does "correlation among repeated measures" work for repeated measures power analysis in g*power?

text: I am using g*power to calculate the a priori power for a 2 group repeated measures anova with a continuous dependent variable with 3 time points. I think that higher values for the "correlation among rep measures" should indicate a better test (i.e., high test-retest correltaion) and thus yield better power and a lower total sample size; however the opposite seems to be true: Test family: F tests; Statistical test: ANOVA: repeated measures, between factors Effect size f: 0.25 alpha error probability: 0.05 Power: 0.8 Number of groups: 2 Number of measurements: 3 Correlation among rep measures: 0.8 (a reasonable value for these neuropsychological assessments of cognitive function) -> total sample size: 112 Correlation among rep measures: 0.5 (a poor test would have a value this low) -> 86 What am I not understanding? Thanks! Henry


title: Understanding productivity or expenses over time without falling victim to stochastic interruptions

text: Help me here, please. Perhaps before even giving me an answer you may need to help me ask the question. I have never learned about time series analysis and do not know if that is indeed what I need. I have never learned about time smoothed averages and do not know if that is indeed what I need. My statistics background: I have 12 credits in biostatistics (multiple linear regression, multiple logistic regression, survival analysis, multi-factorial anova but never repeated measures anova). **So please look at my scenarios below. What are the buzzwords I should be searching for and can you suggest a resource to learn what I need to learn?** I want to look at several different data sets for totally different purposes but common to all of them is that there are dates as one variable. So a couple of examples spring to mind: clinical productivity over time (as in how many surgeries or how many office visits) or electric bill over time (as in money paid to electricity company per month). For both of the above the near universal way to do it is to create a spreadsheet of month or quarter in one column and in the other column would be something such as electricity payment or number of patients seen in the clinic. However, counting per month leads to a lot of noise that has no meaning. For instance, if I usually pay the electricity bill on the 28th of every month but on one occasion I forget and so I only pay it 5 days later on the 3rd of the next month then one month will appear as if there was zero expense and the next month will show ginormous expense. Since one has the actual dates of payment why would one purposefully throw away the very granular data by boxing it into expenses by calendar month. Similarly if I am out of town for 6 days at a conference then that month will appear to be very unproductive and if those 6 days fell near the end of the month, the next month will be uncharacteristicaly busy since there will be a whole waiting list of people who wanted to see me but had to wait till I returned. Then of course there are the obvious seasonal variations. Air conditioners use a lot of electricity so obviously one has to adjust for summer heat. Billions of children are referred to me for recurrent acute otitis media in the winter and hardly any in the summer and early fall. No child of school age gets scheduled for elective surgery in the first 6 weeks that schools return following the long summer vacation. Seasonality is just one independent variable that affects the dependent variable. There must be other independent variables some of which can be guessed and others that are not known. A whole bunch of different issues crop up when looking at enrollment in a longstanding clinical study. What branch of statistics lets us look at this over time by simply looking at events and their actual dates but without creating artificial boxes (months/quarters/years) that do not really exist. I thought of making the weighted average count for any event. For instance number of patients seen this week is equal to 0.5*nr seen this week + 0.25*nr seen last week + 0.25*nr seen next week. I want to learn more about this. What buzzwords should I be searching for?


title: Fitting GLMMs for binomial data produces Error and Warning messages when using different fitting procedures in R.

text: I am trying to fit GLMM's to my data using the glmer function available in R's lme4 package. The data is available under the name block1and2 at: https://onedrive.live.com/redir?resid=1B727FC7180E87DF%21119 I keep getting error or warning messages. Can anybody help me get rid of them. There is probably a problem with the low number of Virus positive samples.Can these problems be resolved using GLMMs or should I think of other statistical tests? N.B. this is a similar problem I had posted about concerning a similar data set and similar analyses (After trying various optimzers, model simplification running more iterations, when fitting GLMMs, R still produces warning messages) Tried fitting the model: `Line.glmer<-glmer(Virus_DNA~Line+(1|Block/Day_of_Analyses),family=binomial,data=data)` Only to get the error message `Error: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate` When running the code, with more iterations: `Line.glmer<-glmer(Virus_DNA~Line+(1|Block/Day_of_Analyses),family=binomial,data=data,control=glmerControl(optCtrl=list(maxfun=1e9)))` I get the same message. When running the code with the optimizer bobyqa: `Line.glmer<-glmer(Virus_DNA~Line+(1|Block/Day_of_Analyses),family=binomial,data=data,control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=1e9)))` I get the warning messages: `Warning messages: 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00163126 (tol = 0.001, component 5) 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge: degenerate Hessian with 2 negative eigenvalues` When running the code with the optimizer Nelder_Mead: 'Line.Nelder_Mead.glmer<-glmer(Virus_DNA~Line+(1|Block/Day_of_Analyses),family=binomial,data=data,control=glmerControl(optimizer="Nelder_Mead",optCtrl=list(maxfun=1e9)))' I get the error message: `Error: (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate` I have also tried the optimizer "nlminb" from the package "optimx". With the code `Line.nlminb.glmer<-glmer(Virus_DNA~Line+(1|Block/Day_of_Analyses),family=binomial,data=data, control=glmerControl(optimizer="optimx",optCtrl=list(method="nlminb")))` Getting the warning messages: `Warning messages: 1: In optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings. This can cause poor performance in optimization. It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA. 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00184741 (tol = 0.001, component 5) 3: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge: degenerate Hessian with 2 negative eigenvalues` Similarly with the optimizer "nlminb" from the package "optimx": `Line.LBFGSB.glmer<-glmer(Virus_DNA~Line+(1|Block/Day_of_Analyses),family=binomial,data=data, control=glmerControl(optimizer="optimx",optCtrl=list(method="L-BFGS-B")))` Produces the warning messages: `Warning messages: 1: In optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower, : Parameters or bounds appear to have different scalings. This can cause poor performance in optimization. It is important for derivative free methods like BOBYQA, UOBYQA, NEWUOA. 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.00330343 (tol = 0.001, component 5) 3: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: large eigenvalue ratio - Rescale variables?`


title: KKT versus unconstrained formulation of lasso regression

text: L1 penalized regression (aka lasso) is presented in two formulations. Let the two objective functions be $$ Q_1 = \frac{1}{2}||Y - X\beta||_2^2 \\\ Q_2 =\frac{1}{2}||Y - X\beta||_2^2 + \lambda ||\beta||_1. $$ Then the two different formulations are $$ \text{argmin}_\beta \; Q_1 $$ subject to $$ ||\beta||_1 \leq t, $$ and, equivalently $$ \text{argmin}_\beta \; Q_2. $$ Using the Karush-Kuhn-Tucker (KKT) conditions, it's easy to see how the stationarity condition for the first formulation is equivalent to taking the gradient of the second formulation and setting it equal to 0. What I can not find, nor figure out, is how the complementary slackness condition for the first formulation, $\lambda\left(||\beta||_1 - t\right) = 0$, is guaranteed to be fulfilled by the solution to the second formulation.


title: Set of expected values

text: I have a doubt, I'd like to calculate a the expected value of a set, let suppose I have a set of n points, every point $x_i$ has a probability $p_i$, the expected value of this set is the sum of all ($p_i * x_i$). However, I'd like to select k points, what is the expected value to this k points? Trying to clarify a lit bit more, I'll give a example, let suppose we have a set of the heights of the women: $X$= {69, 69, 66, 68, 71, 65, 67, 66, 66, 67, 70, 72}, we can see the expected value as: $\sum_{i}x_i p_i$ = (69 + 69 + 66 + 68 + 71 + 65 + 67 + 66 + 66 + 67 + 70 + 72)/12 = 67.9 here $p_i$ = 1/12 67.9 is the expected value whether I get a woman randomly, however I'd like to know what expected value is whether I get k women randomly.


title: Adding proportions as features to logistic regression

text: I'm setting up a logistic regression that models a probability $\mathbb{P}_t$ of certain events. The probability is changing over time and I want to add proportions of past observations that estimate this probability. Suppose I have a proportion $q = \frac{\text{#positives}}{\text{#total}}=\frac{b}{c}$ that estimates $\mathbb{P}_t$ if times are relatively close. How do I add it as a feature in a logistic regression model? **My guess** Because the logistic regression model is $$ \log \frac{p}{1-p} = x \cdot \beta,$$ 1) my intuition says add $\log q$ and $\log(1-q)$. Then our logistic regression model would clearly encompass the baseline model of $p=q$. In the case where the coefficient of $\log q$ is $+1$ and the coefficient of $\log(1-q)$ is $-1$ and other coefficients in $\beta$ are $0$. 2) Another option, remembering that $q=\frac{b}{c}$, is to express \begin{align} \log{q}&= \log{b}-\log{c}\\\ \log{1-q}&=\log \frac{c-b}{c}=\log{(c-b)} - \log{c}. \end{align} Here I would be tempted to add $\log{(b + 1)}, \log{(c + 1)}, \log{(c - b + 1)}.$ as features.


title: Significance testing between two sets of Pearson's correlation coefficients

text: I have been reading through questions here to find an answer to my question but didn't find a satisfactory one. Here is my scenario: I have a system computing an algorithm on each of 100 units (so basically, the system will generate a 100 values each corresponds to one unit). I have two different algorithms to predict the value the system will generate for a given unit. I randomly generated 120 subsets of the 100 units. For a prediction algorithm, I computed the predictions of each subset then computed Pearson's correlation coefficient between actual and predicted values for each subset. Eventually, I got 120 correlations r for each of the prediction algorithms. My questions are:   1. Is it meaningful to compute the average correlation for each prediction algorithm over the 120 correlations? I am using Fisher's Z to compute the average as explained here. So can I safely say that for these 100 units, algorithm X managed to predict value by correlation of = average correlation?   2. I need to compare the significance of difference in predictive power of the two algorithms in such scenario, how can I do this having 120 correlation values per algorithm? can I normally apply a t-test in such case? Thanks in advance,


title: Defining a simple linear regression of reaction times against stimulus number

text: I have reaction times for each of nine stimuli (the stimuli being the numbers 1 to 9), all 9 variables being arranged as columns, and with subjects across rows. I would like to regress those scores (Y) against their respective numbers (X), but have trouble understanding how to define this regression in my stats package (Statistica), whose help suggests that each of the two variables (X and Y) should be selected from the variables list. X in my case would just be the numbers 1 to 9, whereas Y would be each of the 9 columns (variables). How should this simple regression be defined? Should it be a special type of regression given that the X variable is really categorical rather than continuous? Thanks!


title: Multi Factor analysis with FactoMineR

text: I have a database with different variables, both categorical and numerical. I wanted to analyse these using MultiFactor Analysis with FactoMineR, with the main idea to obtain which of them are the most descriptive, and the correlation between them. The problem is that the resulting graphs only performance separate graph ones for categorical variables and others for numerical variables. Is there any method or a better R package that I can use to obtain common results?


title: Can I perform Cox regression on left truncated records?

text: Is a Cox proportional hazards regression model appropriate for records which are **left truncated**? I'm developing a model that predicts risk of hospitalization during a two year time period. Some members of my study population have been hospitalized prior to the study period (left censoring) based on records prior to the study period. No problem, there are exisitng techniques to adjust for left censoring. However, there are also left truncated members for whom we have no information whether they have been hospitalized or not prior to the study. Some of these subjects float through the 2 year period without being hospitalized. Is there any way to retain the left truncated subjects in the Cox regression model? Should I consider another methodology, such as logistic regression? Edit: To clarify, here are some record examples for clarification. I apologize for the formatting - I don't know how to create tables in the editor. As per stnd methodology for a recurring event Cox regression, there are multiple records for member 1 (CT0001) associated with each interval ending in an IP admission or the end of the time period. The Start and End columns display the beginning and end dates of each interval, or . if missing. Member 2 (CT0002) didn't go to the hospital during the study period. But we still want member 2 in the dataset, since we're predicting risk of hospitalization for the entire population based on demographic variables even if there was no IP visit.               **MemberID    IP_admit      IP_disch    Start       End         Censoring/Truncation?**     ----------     CT0001        10/1/2010     10/20/2010  .           .           .     CT0001        3/1/2011      3/15/2011   .           3/1/2011    Left censored     CT0001        7/1/2011      7/5/2011    3/16/2011   7/1/2011    .     CT0001        11/1/2011     11/15/2011  7/6/2011    11/1/2011   .     CT0001        .             .           11/16/2011  .           Right censored     CT0002        .             .           .           .           Left truncation


title: VAR model with zero coefficients

text: ll, I'm working with a bivariate time series $(X_{t},Y_{t})$. Looking at the two time series separately, $X_{t}$ appears to be white noise. This is supported by looking at the empirical ACF and PACF of $X_{t}$ and by performing the Ljung-Box-Test for various lags (very high p-values). My model for $X_{t}$ would therefore be: $$X_{t} = \mu + \epsilon_{t}$$ with $\epsilon_{t} \sim WN(\mu, \sigma_{\epsilon})$. Looking at the empirical ACF and PACF of the other time series $Y_{t}$ I get significant values of both ACF and PACF even for relatively high lags. This could suggest an ARMA model but using AIC and a maximum model order of $p,q = 5$ points to an AR(1) model, like so: $$Y_{t} = \alpha_{1}Y_{t-1}+\epsilon^{y}_{t}$$ Finally, I had a look at the empirical cross-correlation function of $Y_{t}$ and $X_{t}$. It shows a significant value for a lag of $h = -1$, that is a significant correlation between $Y_{t}$ and $X_{t+1}$. It would therefore seem appropriate to account for this by letting $X_{t}$ depend on $Y_{t-1}$, yielding something like $$X_{t} = \mu + \beta_{1}Y_{t-1}+ \epsilon_{t}$$ I'm wondering how I would go about fitting such a model (preferably using R). Would I have to use a package for VAR models and specify the $X_{t-1}$ coefficient for $Y_{t}$ as zero? Thanks!


title: Interaction effect 2x2

text: I have a little question about a 2x2 design that keeps on confusing me. If I have this kind of 2x2 design where I can have these combinations in the conditions. No/No, Yes/No, No/Yes, Yes/Yes. The interaction effect is the effect of the two IVs, right? As the Yes/Yes condition consist of both IVs, is the Yes/Yes condition here the same as the interaction effect of A and B? ![2x2](http://i.stack.imgur.com/9eGfa.png)


title: Approximation for Beta distribution when alpha is less than 10

text: I know that we can approximate Beta distribution to Normal distribution when the values of alpha and beta are large numbers. In my problem alpha lies between 1 and 10, beta is always greater than 1000. Is there any distribution (for such a combination of alpha and beta) that we can approximate it to, with a simpler closed form so that it will be easy to calculate CDF and InverseCDF? (I have tried Kumaraswamy distribution but did not get good results due to the large values of beta.) **EDIT** "easy"- in terms of computation time in JAVA/PHP. (I need an approximate distribution since it takes lot of time to compute the InverseCDF of Beta distribution(0.05 ms in java is also very slow for my needs). For alpha greater than 10 i approximated Beta distribution to Normal distribution and computed InverseCDF like  this  which included just math. So i did not use any other library functions like Normal distribution  or  Beta distribution, which reduced the processing time drastically.)


title: Normalization in SVM

text: I have applied `libsvm` with a linear kernel to a set of instances and I have obtained a 68 % success:               instance1 : f11, f12, f13, f14     instance2 : f21, f22, f23, f24     instance3 : f31, f32, f33, f34     instance4 : f41, f42, f43, f44     ..............................     instanceN : fN1, fN2, fN3, fN4      Taking the same set of instances but multiplying each value (`f11*1000 ... fN4*1000`) I have obtained a 90% of success. However, I the multiplied by 1000 is normalized the percentatge of success turn to 68%. It seems something related with normalization but I don't know which it is the reason.


title: Does Fisher's factorization theorem provide the pdf of the sufficient statistic?

text: From Wikipedia > Fisher's factorization theorem or factorization criterion provides a > convenient characterization of a sufficient statistic. If the probability > density function is $ƒ_θ(x)$, then $T$ is sufficient for $θ$ if and only if > nonnegative functions $g$ and $h$ can be found such that $$ f_\theta(x)=h(x) > \, g_\theta(T(x)), \,\\! $$ i.e. the density $ƒ$ can be factored into a > product such that one factor, $h$, does not depend on $θ$ and the other > factor, which does depend on $θ$, depends on $x$ only through $T(x)$. I was wondering if $\frac{g_\theta(t)}{c}$, where $c := \int g_\theta(t) dt$, is the pdf of $T(X)$ when the pdf of $X$ is $ƒ_θ(x)$? Thanks and regards!


title: Ancillary statistic not containing information about sample distribution?

text: From a note by Jun Shao > If $V(X)$ is a nontrivial ancillary statistic, then $σ(V(X)) ⊂ σ(X)$ is a > nontrivial σ-field that **does not contain any information about** $P$. I was wondering in what sense "the σ-field of statistic $V(X)$ does not contain any information about the distribution $P$ of $X$" means?   1. I learned that σ-field itself represents information, but not sure about what it means by σ-field of $V(X)$ not containing information about distribution of $X$.   2. Does "information" here mean the fisher information? Thanks and regards


title: Finding temporal patterns that predict "life length"

text: I'm doing research on system usage. It is an online system that provide certain information to users, and I want to predict for how long users use the system. I already noticed that if people don't use it for 2 weeks, then it is unlikely that they will use it again. So my "dependent variable" is the length of time between user registration and the point in which they do not do anything for two weeks. I have a log of every visit to the system and every action the user did. I already did some "aggregate" analysis (e.g. average monthly usage, etc.) and created a random forest on a bunch of such variables. I now look for a more "temporal" approach, that is taking into account when certain types of actions occurred and what are the relations among those times. I am not an expert on this but I am sure that there is plenty of practical and theoretical work that has been done. So any reference or suggestion (or leading questions) that will direct me will be highly appreciated.


title: How to check the features which are selected by LASSO

text: I am using LASSO (glmnet) to do feature selection. However, how can I check which features are selected?


title: Does it make sense to add a quadratic term but not the linear term to a model?

text: I have a (mixed) model in which one of my predictors should a priori only be quadratically related to the predictor (due to the experimental manipulation). Hence, I would like to add only the quadratic term to the model. Two things keep me from doing so:   1. I think I read somehwere that you should always include the lower order polynomial when fitting higher order polynomials. I forgot where I found it and in the literature I looked at (e.g., Faraway, 2002; Fox, 2002) I cannot find a good explanation.   2. When I add both, the linear and quadratic term, both are significant. When I add only one of them, they are not significant. However, a linear relation of predictor and data is not interpretable. The context of my question is specifically a mixed-model using `lme4`, but I would like to get answers that could explain why it is or why it is not okay to inlcude a higher order polynomial and not the lower order polynomial. If necessary I can provide the data.


title: How do I interpret a variable in a regression if it is squared, but I don't include the linear variable?

text: I have a regression model with 3 different variables. One variable is only significant if I square it, and not include the linear version i.e. $Y = \beta_1X_1 + \beta_2X_2 + \beta_3X_3^2$ How should I interpret the $X_3^2$ variable?


title: With K-Fold cross-validation, do you average the model?

text: When performing K-fold cross validation, I understand that you obtain the accuracy metrics by pointing all the folds except one at that one fold and make predictions, and then repeat this process K times. You can then run accuracy metrics on all of your instances (precision, recall, % classified correctly), which ought to be the same as if you calculated them each time and then averaged the result (correct me if I'm wrong.) The end result you want is a model. Do you average the models obtained to make your set of K predictions to end up with the model that has the accuracy metrics obtained by the above method? Apologies if something is really, really wrong or just patently obvious -- just correct it. Thanks


title: Clustering of documents that are very different in number of words

text: I have a corpus of 643 documents with different sizes and my goal is to cluster them according their topics and label each cluster with semantic name for its main topic. I have tired different approaches for clustering including:   1. Clustering based on cosine distance between the tf-idf vector of weights representing words in each document.   2. Clustering based on finding the tf-idf of words in each document, then use LSA (Latent Semantic Analysis) to find an approximate representation of the document in the LSA.. Then I forumulate a number of query terms. Each query represent a topic and rank documents according to their similarity to the different queries.   3. Clustering based on finding the tf-idf of words in each document. Then use PCA to find the principal components and then choose the first 10 components (contributing to more than 97% of variance between documents) and perform clustering based on it. I use Hierarchical agglomerative Clustering and I employ the Silhouette value to get an indication for the optimal cut level of the resultant dendrogram and therefore the number of clusters. However, I don't think that the clustering results really represent documents belonging to different topics. Most of time I end with having one very big cluster (> 200 documents) and other few clusters with small number of documents each. I have noticed that the documents are different in their size significantly, The quantiles for the number of words / document are as follow: 25% of documents have 6 or more words. 50% of documents have 12 or more words 75% of documents have 21 or more words 3 documents have more than 200 words. (The one with maximum words have 608 words). Do you think that the big variation for number of words can be a reason for a problem in clustering. ?? and If Yes, would you please suggest possible methods to solve such problem.


title: burn in for Metropolis Hastings MCMC

text: I was wondering if there is a principled way to figure out how many samples to discard during the MH-MCMC burn-in stage. So, as I understand it, the initial samples can introduce bias in the computation of expectation in the estimated posterior distribution especially if they are really far from the mean. I was wondering, say I do a MAP estimation and try to find some local mode of the distribution to initialise my MH-MCMC. Would that mean I do not have to do the burn in (assuming this mode is in some high density region of the distribution)?


title: How to compute sample size to compare two diagnostic tests

text: I will be performing two diagnostic tests (one is the gold standard, one is novel) on the same subject aiming to establish sensitivity, specificity, PPV and NPV. What formula may be used to compute for a good sample size? What is better way to establish correlation of the two tests? I understand a Bland-Altman plot may not be used since my tests check for two different parameters (and diagnosis lies on the interpretation of the parameters)?


title: When does Naive Bayes perform better than SVM?

text: In a small text classification problem I was looking at, Naive Bayes has been exhibiting a performance similar to or greater than an SVM and I was very confused. I was wondering what factors decide the triumph of one algorithm over the other. Are there situations where there is no point in using Naive Bayes over SVMs? Can someone shed light on this?


title: Time Series Econometrics: Cointegration methods for series with mixed degree of integration

text: I am performing a time series analysis on dataset ranging from 1974-2008. I have performed Augmented Dickey Fuller tests and Phillips Perron tests to check the stationarity/ order of integration of my series. These tests both confirmed that all series are I(1) with the exception of GDP per capita which is I(2). The country in question is Ireland. For the purposes of my analysis I had been planning on using the ARDL Bounds Testing approach to cointegration (Pesaran & Shin, 1999) but this method is not possible given the presence of an I(2) variable. My question is: Does anyone know of any approach to cointegration which allows for the mix of I(1) and I(2) variables? I am unsure as to how to proceed with my analysis given this issue and thus far have not come across any solutions. Any help at all would be greatly appreciated.


title: Does talking about a probability of future events make you a Bayesian?

text: Okay, so I sent a confession of love to a girl detailing several reasons and coming up with a 90% probability that I will be with her together for my whole life. But I then thought that I was supposed to be a frequentist and not a Bayesian, because I cannot start up 100 samples of me and her in a virtual machine and build prediction intervals for breakup time. Is this true? I always thought that probability is utterly useless if it cannot allow you to, say, claim that a coin you never flipped has a probability 1/6. On the other hand it seems wrong to suppose, say, that my probability of being forever with my beloved would decrease due to a revelation that she is a spy for the NSA. What exactly is the difference between frequentists and Bayesians on probabilities of unknown events? My frequentist stat teacher simply handwaves about it and say that we just call it a 90% "confidence" without daring to call it a probability and violating Þe Olde Frequentist School. That seems so unrigorous, and besides when calculating my 90% "confidence", I used formulas with, guess what, $Pr(...)$ in them.


title: Convergence in probability and $L_2$ for normal random variables

text: In an answer here: Convergence of identically distributed normal random variables, the following lemma is mentioned: > **Lemma** : Let $X_1, X_2, \ldots$ be a sequence of zero-mean normal random > variables defined on the same space with variances $\sigma^2_n$. Then, $X_n > \to X_\infty$ both in probability and in $L_2$ (and hence also in > distribution) if and only if $\sigma^2_n \to \sigma^2 < \infty$. In this > case, the limit $X_\infty$ is also normally distributed with mean zero and > variance $\sigma^2$. I'd like to find out more about this lemma. Does it have a name?


title: Generating Random Vectors with Arbitrary Marginal Distributions via NORTA

text: When generating random variates from different marginal distributions using the NORTA (Normal-to-Anything) method, as described in Cario & Nelson 2007, why is $\varrho$ required? > To adjust this correlation, we can restrict attention to adjusting $\ > E[X_iX_j]$, since > > $$ Corr[X_i; X_j ] = \frac{E[X_iX_j]-E[X_i]E[X_j]}{\sqrt{Var[X_i]Var[X_j]}} > $$ > > and $\ E[X_i]$, $\ E[X_i]$, $\ Var[X_i]$ and $\ Var[X_i]$ are fixed by $\ > F_{X_i} $ and $\ F_{X_j} $. Then, since $\ (Z_i, Z_j)$ has a standard > bivariate normal distribution with correlation $\ Corr[Z_i; Z_j ] = (Z_i, > Z_j)$, we have > > $$ E[X_iX_j] = E \Bigl\\{ F^{-1}_{X_i}[\Phi(Z_i)]F^{-1}_{X_j}[\Phi(Z_j)] > \Bigl\\} $$ $$ = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} > F^{-1}_{X_i}[\Phi(Z_i)]F^{-1}_{X_j}[\Phi(Z_j)]\varphi_{\varrho > z(i,j)}(z_i,z_j)dz_idz_j $$ > > where $\ \varphi_{\varrho z(i,j)}(z_i,z_j) $ is the standard bivariate > normal probability density function (pdf) with correlation $\ \varrho z(i,j) > $. **We are only interested in distributions for which this expectation > exists** Isn't it possible to create the desired variates by feeding the required parameters straight into a Gaussian copula? For instance, a -0.5 correlation with U(10) and Expo(10). MATLAB example below.               corrZ = [1 -0.5;-0.5 1];     size = 200;     Z = mvnrnd([0 0], corrZ, size);      X(:,1) = expinv( normcdf(Z(:,1)), 10);      X(:,2) = unidinv(normcdf(Z(:,2)), 10);      X     mean(X)     scatter(X(:,2),X(:,1))      Does $\varrho$ preserve the variances and means of the marginals? What am I misunderstanding? One other question: what does it mean when the value found for $\varrho$ lies outside (-1,1)? I have some marginals that make the root finding for the $\varrho$ unsolvable. Why does this occur? I can clarify further if required.


title: Approximate function observed through noise

text: I have a function $f : [a,b] \rightarrow \mathbb{R}$ that I can observe through some noise, i.e. I can only directly measure $f(x) + \eta$ where $\eta$ is some random noise with mean 0. I have a large number of noisy observations of this function, $f_i$, at _different_ arguments $x_i$. What is a good way to approximate the function for any $x \in [a,b]$? If I had several observations of $f(x_0)$ at the _same_ $x_0$ then I could average over these to estimate the true $f(x_0)$. But while I have a relatively large number of observations, each one is at a different argument $x_i$. What is the best course of action in this case? * * * I have no background in statistics so excuse me if this is a basic question. What I'd like t know is if there's some sane theory to handle these situation (which would be much better than myself coming up with some ad-hoc solution such as moving averages). The names of some methods I could use would be sufficient, I can look up the details. Note: I would like to use the approximated $f$ mainly for things such as optimization and finding zero crossings. * * * My data is similar to this artificially generated dataset (Mathematica code):               n = 500;     data = {#, Sin[#] + RandomReal[.1 {-1, 1}]} & /@ RandomReal[{0, 2 Pi}, n];          ListPlot[data]      ![enter image description here](http://i.stack.imgur.com/7qoc5.png)


title: What is "the smallest solution of an equation" mean?

text: Suppose you have an quadratic equation (n) of `q`, what's the meaning of this statement "...This implies that q is the smallest solution to equation (n)"? What is the property of `smallest solution` in an equation?


title: Using Wilcoxon-Mann-Whitney test for comparing two population of different sizes

text: I am using Wilcoxon-Mann-Whitney test for comparing two populations. Unfortunately, sizes of my population are different; one has size 100000 and the other 6000. Can I use this test to compare these populations? If not, is there any alternative test? Thanks!!


title: Multilevel exploratory factor analysis

text: We are creating a scale to measure the 'strength of tie' between two people. We conducted a survey where each person was asked to name a number of other people they know. They were then asked to respond to a scale concerning each person. As a result we have a 23-item scale with 10 or more responses from each respondent. This is a brand new scale, so we want to know the structure. We expect it to be multi-dimensional. Can you suggest a direction I can take for this? I use R primarily, so that would be preferred, but if I need to use a different software package I can figure it out.


title: Poisson Approximation to Binomial

text: I have a Poisson approximation to binomial question, posted below. I'm not too sure if I'm using the proper formula: $$P(x) = e^{-np}(np)^x/x!$$ . Of anyone can tell me if I'm doing this right, that would be great. Again these are just practice questions not homework. > Q: A salesperson has found that the probability of making a sale on a > particular product manufactured by him or her company is .05. If the > salesperson contacts 140 potential customers, what is the probability he or > she will sell at least 2 of these products? Use and justify Poisson > approximation to Binomial. What I'm doing: $$P(x) = e^{-140 (.05)}(140*.05)^2/2!$$ .


title: Error when running glmnet in multinomial

text: **The problem mentioned in this question is fixed in version 1.7.3 of the R package glmnet.** I'm having some problems running glmnet with family=multinomial, and was wondering has encountered something similar or might be able to tell me what I'm doing wrong. When I put my own dummy data in, the error "Error in apply(nz, 1, median) : dim(X) must have a positive length" gets reported when I run `cv.glmnet`, which apart from saying "it didn't work" wasn't hugely informative to me.               y=rep(1:3,20) #=> 60 element vector     set.seed(1011)     x=matrix(y+rnorm(20*3*10,sd=0.4),nrow=60) # 60*10 element matrix     glm = glmnet(x,y,family="multinomial")   #=> returns without error     crossval = cv.glmnet(x,y,family="multinomial")   #=> Error in apply(nz, 1, median) : dim(X) must have a positive length     crossval = cv.glmnet(x,y,family="multinomial",type.measure="class")   #=> Error in apply(nz, 1, median) : dim(X) must have a positive length     crossval = cv.glmnet(x,y,family="multinomial",type.measure="mae")   #=> Error in apply(nz, 1, median) : dim(X) must have a positive length     cvglm = cv.glmnet(x,y,family="multinomial",lambda=2)   #=> Error in apply(nz, 1, median) : dim(X) must have a positive length      Here's a visual description of the problem I was trying to get glmnet to solve, if that helps:               my_colours = c('red','green','blue')     plot(x[,1],x[,2],col=my_colours[y])      I'm able to run the example code from the package docs, which makes me suspcious that I'm either misunderstanding something or that there is a bug in glmnet.               library(glmnet)     set.seed(10101)     n=1000;p=30     x=matrix(rnorm(n*p),n,p) #=> 1000*30 element matrix     beta3=matrix(rnorm(30),10,3)     beta3=rbind(beta3,matrix(0,p-10,3))     f3=x%*% beta3     p3=exp(f3)     p3=p3/apply(p3,1,sum)     g3=rmult(p3) #=> 1000 element vector     set.seed(10101)     cvfit=cv.glmnet(x,g3,family="multinomial")      This is using R version 2.13.1 (2011-07-08) and glmnet 1.7.1, though I can generate the same problem on R 2.14.1. Any ideas people?


title: Calculate variance explained by each predictor in multiple regression using R

text: I have run a multiple regression in which the model as a whole is significant and explains about 13% of the variance. However, I need to find the amount of variance explained by each significant predictor. How can I do this using R? Here's some sample data and code:               D = data.frame(         dv = c( 0.75, 1.00, 1.00, 0.75, 0.50, 0.75, 1.00, 1.00, 0.75, 0.50 ),         iv1 = c( 0.75, 1.00, 1.00, 0.75, 0.75, 1.00, 0.50, 0.50, 0.75, 0.25 ),         iv2 = c( 0.882, 0.867, 0.900, 0.333, 0.875, 0.500, 0.882, 0.875, 0.778, 0.867 ),         iv3 = c( 1.000, 0.067, 1.000, 0.933, 0.875, 0.500, 0.588, 0.875, 1.000, 0.467 ),         iv4 = c( 0.889, 1.000, 0.905, 0.938, 0.833, 0.882, 0.444, 0.588, 0.895, 0.812 ),         iv5 = c( 18, 16, 21, 16, 18, 17, 18, 17, 19, 16 ) )     fit = lm( dv ~ iv1 + iv2 + iv3 + iv4 + iv5, data=D )     summary( fit )      Here's the output with my actual data:               Call: lm(formula = posttestScore ~ pretestScore + probCategorySame +          probDataRelated + practiceAccuracy + practiceNumTrials, data = D)          Residuals:         Min      1Q  Median      3Q     Max      -0.6881 -0.1185  0.0516  0.1359  0.3690           Coefficients:                       Estimate Std. Error t value Pr(>|t|)      (Intercept)        0.77364    0.10603    7.30  8.5e-13 ***      iv1                0.29267    0.03091    9.47  < 2e-16 ***      iv2                0.06354    0.02456    2.59   0.0099 **      iv3                0.00553    0.02637    0.21   0.8340      iv4               -0.02642    0.06505   -0.41   0.6847      iv5               -0.00941    0.00501   -1.88   0.0607 .       --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1          Residual standard error: 0.18 on 665 degrees of freedom      Multiple R-squared:  0.13,      Adjusted R-squared:  0.123      F-statistic: 19.8 on 5 and 665 DF,  p-value: <2e-16      This question has been answered here, but the accepted answer only addresses uncorrelated predictors, and while there is an additional response that addresses correlated predictors, it only provides a general hint, not a specific solution. I would like to know what to do if my predictors are correlated.


title: How to calculate variance explained by an environmental variable in the model?

text: I have done some survival analysis, and used the following equation to calculate variance explained by environmental variables in the survival models: (DEV.-DEVenv)/(DEV.-DEVtime) where DEV. is the deviance of the model assuming constant survival rates, DEVenv is the deviance of the model assuming a linear effect of an environmental variable, and DEVtime is the deviance of the model assuming time-dependent survival rates. I should calculate variance explained by the variables temperature and salinity in the following model: Model 1: location + temperature + salinity, deviance = 318.2 In addition, I have models Model 2 (constant survival rate): location, deviance = 328.5 Model 3 (time- dependent survival rate): location + time, deviance = 318.7 However, if I calculate variance explained by temperature + salinity in the model 1, I get the following result: (328.5-318.2)/(328.5-318.7) = 1.05 (i.e. more than 100%!). Have I calculated something wrong to get this result? Should I consider temperature and salinity separately and how can one do this when both variables are in the same model?


title: Plotting a categorical response as a function of a continuous predictor using R

text: I would like to plot the relationship between a binary categorical response variable and a continuous predictor to study its shape. The goal is to prep a logistic regression. This image may clarify: ![Example Plot](http://i.stack.imgur.com/pUGv5.png) I have access to Minitab and R and would greatly appreciate any insight on how to recreate this histogram or alternatives that may do just as well.


title: matlab princomp latent

text: I have a feature vector containing 17 features. I need to take only features containing 98% variance. So i did a princomp and printed out latent and got this                   0.9082         0.9824         0.9999         0.9999         1.0000         1.0000         1.0000         1.0000         1.0000         1.0000         1.0000         1.0000         1.0000         1.0000         1.0000         1.0000         1.0000      The first value -> 0.9082 indicates the variance value of the first principle component am i correct ? Now if i need to select featurs containing 98% variance do i select only                   0.9824         0.9999         0.9999      What does 1.0 mean? does this mean 100% variance ? Could somebody please clear my confusions. PCA code               [COEFF,SCORE,latent] = princomp(X);     cumsum(latent)./sum(latent)


title: Asymptotics of 0-1 classification loss

text: I am interested in training a simple binary linear classifier. That is, I will find a vector of weights $\bf w$ such that I can predict the class of new example by the sign of $f(x) = w^T x$. I will find the weights in the standard way, by minimizing the empirical risk of some loss function $\sum_i L(w^T x^{(i)}, y^(i))$ Many different loss functions have been suggested to solve this problem, such as the hinge loss $L(w^T x, y) = \max(0,1-y\ w^T x)$, the "logistic" loss $L(w^T x, y) = \log(1 + \exp(-y\ w^T x))$, and the "0-1" loss $L(w^T x, y) = I(y = \text{sign}( w^T x))$. The 0-1 loss is rarely used in practice due to non-convexity. However, suppose that we were able to minimize it. If the data is "noisy" or non-separable, we can expect the 0-1 loss to outperform the others (in terms of classification rate) given infinite data, since it is robust to outliers. My question is this: does the 0-1 loss have a convergence _rate_ equal to the others? (Edit: Here by convergence rate, I mean in the traditional statistical sense of the expected squared difference between the recovered weights and the true weights, though I would be interested in other notions of convergence as well.) I would be interested in an answer in any setting, but the simplest might be to assume that the data _is_ in fact separable. My intuition tells me that here the logistic or hinge loss will outperform the 0-1 loss, but I can't quantify this. Of course, if there were a result for a different classifier than linear, that would be of interest as well. Thanks for any insight!


title: How do I train a (logistic?) regression in R using L1 loss function?

text: I can train a logistic regression in `R` using               glm(y ~ x, family=binomial(logit)))      but, IIUC, this optimizes for log likelihood. Is there a way to train the model using the linear ($L_1$) loss function (which in this case is the same as the total variation distance)? I.e., given a numeric vector $x$ and a bit (logical) vector $y$, I want to construct a monotonic (in fact, increasing) function $f$ such that $\sum |f(x)-y|$ is minimized. See also   * How do I train a logistic regression in R using L1 loss function?


title: Is chi-squared the right method to compare time periods?

text: We manufacture foobars. In July, 91% of foobars were defect-free, but in August that figure was 89%. Would chi-squared be the right method to determine if the difference of 2% between July and August is significant? I read somewhere that chi-squared should not be used to compare same variable from different time periods, as it would not qualify as independent variable in this case. **Update** On prompt from @Glen_b, I dug out the passage where I got this idea from: "What is a p-value anyway?" by Andrew Vickers, page 196. In answer to a review question, the book says: _"All of the common statistical tests assume that the data are independent. Applying these tests to non- independent data is a very common error. An obvious example is repeat observations."_ The example data follows, in which sales over a week in two different stores are split into daily sales, so that instead of comparing two data points (weekly sales in each store) we can compare 14 data points (daily sales in each store). The erroneous assumption here is that we have 14 independent data points, while in fact we have two sets of non-independent data points, thus t-test or chi-squared cannot be applied. The author defines independence as: _"...two variables are independent if information about one gives you no information about the other."_ He then says this about the given example: _"...if I tell you Monday's sales figures, you can take a reasonable guess at Tuesday's"._ Therefore, Tuesday sales are not independent from Monday sales. Similarly, in my example August measurements are not independent of July measurements.


title: Balanced Incomplete Block Design(BIBD)

text: If treatment $1$ appears in $r$ blocks and there are $(k-1)$ other treatments in each of those blocks, there are $r(k-1)$ observations in a block containing treatment $1$. But the following Table does not follow the above statement but the table is said to be BIBD: $$Block$$ $$ \begin{array}{l} 1 & 2 &3 &4 \\\ \hline A & A &A & B \\\ B & B &C & C \\\ C & D &D & D \\\ \end{array} $$ here, $k=3$ , $r=3$ But each block doesn't contain $r(k-1)=3(3-1)=6$ observations rather each block is containing $3$ observations. _Where am i doing the mistake?_