I was working through the lab on ridge regression and LASSO in ISLR and I came across a strange behavior in the cv.glmnet function. When I followed the lab as written I got the following   set.seed(1) train &lt;- sample(1:nrow(x), nrow(x)/2) test &lt;- (-train) y.test &lt;- y[test] set.seed(1) cv.out &lt;- cv.glmnet(x[train,], y[train], lambda=grid, alpha=0) plot(cv.out) bestlam &lt;- cv.out$lambda.min bestlam [1] 231.013   For my own benefit I tried it using a different seed (8675309) and got back a different result. Any combination of setting the seeds resulted in different answers. I am assuming this has to do with how the 10-folds are changed with the different seeds, however the different lambda.min can vary so much I am concerned the package might not be stable. Am I missing something?
121663
I am using cv.glmnet to find predictors. The setup I use is as follows:  lassoResults&lt;-cv.glmnet(x=countDiffs,y=responseDiffs,alpha=1,nfolds=cvfold) bestlambda&lt;-lassoResults$lambda.min  results&lt;-predict(lassoResults,s=bestlambda,type="coefficients")  choicePred&lt;-rownames(results)[which(results !=0)]   To make sure the results are reproducible I set.seed(1). The results are highly variable. I ran the exact same code 100 times to see how variable the results were. In the 98/100 runs had one particular predictor always selected (sometimes just on its own); other predictors were selected (co-efficient was non-zero) usually 50/100 times.  So it tells me that each time the cross-validation is running it's going to probably select a different best lambda, because of the initial randomization of the folds matter. Others have seen this problem () but there isn't a suggested solution.  I am thinking that maybe that one which shows up 98/100 is probably pretty highly correlated to all the others? The results do stabilize if I just run LOOCV ($\text{fold-size} = n$), but I am curious why they are so variable when $\text{nfold} &lt; n$.
39980
I'm looking at a regression model where a very large number of possible explanatory variables are being evaluated, and a small number are finally chosen via the lasso method of variable selection. The $\lambda$ tuning parameter in the lasso is chosen by looking at cross-validation forecast performance, which is pretty standard.   However, when I take the list of chosen variables and just run OLS on them, many turn out to be statistically insignificant. That may be perfectly fine if they are jointly significant and the forecast performance is superior to other models (in addition, there would be a question of what the t-test means when you have already screened the variables in a separate step, but I'm leaving that aside).   I'm curious though whether it makes sense to look at statistical significance of individual variables in a model chosen by lasso using CV forecast performance to select the tuning parameter. The problem is that lasso ends up selecting various dummy variables that are only true on small segments of the population and which are insignificant in OLS, and there is a natural question as to whether the model should be judgmentally simplified.
94017
I am interested in utilizing caret for making inferences on a particular data set. Is it possible to do the following:   produce coefficients of a glmnet model I trained in caret. I would like to use glmnet because of the inherent feature selection as I do not believe glm has it? other than the ROC metric, is there another metric that I can utilize to asses fit of the model? Such as adjusted $R^2$?   The purpose of this analysis is to derive some inference on the effects of particular variables, rather than for prediction. I just like the caret package because it's been easy to work with thus far using matrices.
93747
In a data set I have, all of the results of the residual vs. predicted values  appear as parallel lines as pictured below. I interpret this according to my training as a violation of homoscedasticity which makes running a multiple regression untenable.      Here is some context about the data sets, type of analysis used, the design of the study, etc.   The type of analysis I conducted was multiple regression. In SPSS - Analyze > Regression > Linear. The independent (predictor) and dependent (outcome) variables are all ordinal data (though I'm treating them as "scale" in the analysis). Both the independent and dependent variables represent satisfaction scores on a 7-point scale. These satisfaction scores were collected from the same participants at the same point in time (one survey). There are 4 independent variables that represent participants' satisfaction with several more granular dimensions or aspects of a service (e.g., easy to use, effective in meeting needs, etc.) The dependent variable is satisfaction overall with the service.  More about the other results from testing the multiple regression assumptions: (1) The data set passed the "independence of observations" test as measured by a Durbin-Watson statistics which was approximately 2; (2) the independent and dependent variables are approximately linearly related which was based on an inspection of the partial plots produced by the regression analysis.  There are around 20 highly influential data points out of approximately 450 participants as measured by studentized deleted residuals of over +/- 3.    I've tried transforming the data using a log transformation and this did not at all help the apparent heteroscedasticity problem. Plotting the log transformed residuals vs. predicted values yielded the same general pattern in the chart.  The other technique I tried was outlined in this paper  which described providing a regression model with heteroscedasticity-corrected error residuals. The beta coefficients for each predictor variable were essentially the same as the non-heteroscedasticity corrected model, as well as the overall R^2 value for the model.   Some of the predictor variables are highly correlated with one another, but around .6 - .8 - so there may be a slight multicollinearity problem - however, the tolerance and VIF values are fine.  More intuitively, a lot of the participants provided similar satisfaction ratings for all of the predictor variables and the dependent variable. Quite a few on the 7-point scale answered things like " 7 7 7" for the predictors and "7" for the outcome or "1 1 1" and "1". However, this was not always the case.    Looking for some guidance around potential reasons for this, alternative ways to analyze the model, etc. Thanks in advance.
52604
Using R on some data and trying to see whether or not my data is heteroscedastic, I've found two implementations of the Breusch-Pagan test,  (package lmtest) and  (package car). However, these produce different results. What is the difference between the two? When should you choose to use one or the other?  &gt; model &lt;- lm(y ~ x) &gt; bp &lt;- bptest(model) &gt; bp studentized Breusch-Pagan test  data:  model  BP = 3.3596, df = 1, p-value = 0.06681  &gt; ncvTest(model) Non-constant Variance Score Test  Variance formula: ~ fitted.values  Chisquare = 3.858704    Df = 1     p = 0.04948855    These example shows that according to the tests, my data is in one case heteroscedastic and in the other case homoscedastic. I did find  so bptest might be studentized and ncvTest might not be, however, what does this mean then?
10268
I have the following output:   Generalized linear mixed model fit by the Laplace approximation  Formula: aph.remain ~ sMFS2 +sAG2 +sSHDI2 +sbare +season +crop +(1|landscape)    AIC   BIC    logLik deviance  4062  4093  -2022   4044  Random effects: Groups    Name        Variance Std.Dev. landscape (Intercept) 0.82453  0.90804  Number of obs: 239, groups: landscape, 45  Fixed effects:             Estimate Std. Error z value Pr(&gt;|z|)     (Intercept)  2.65120    0.14051  18.868   &lt;2e-16      sMFS2        0.26922    0.17594   1.530   0.1260     sAG2         0.09268    0.14529   0.638   0.5235     sSHDI2       0.28345    0.17177   1.650   0.0989   sbare        0.41388    0.02976  13.907   &lt;2e-16  seasonlate  -0.50165    0.02729 -18.384   &lt;2e-16  cropforage   0.79000    0.06724  11.748   &lt;2e-16  cropsoy      0.76507    0.04920  15.551   &lt;2e-16   Correlation of Fixed Effects:            (Intr) sMFS2  sAG2   sSHDI2 sbare  sesnlt crpfrg sMFS2      -0.016                                           sAG2        0.006 -0.342                                    sSHDI2     -0.025  0.588 -0.169                             sbare      -0.113 -0.002  0.010  0.004                      seasonlate -0.034  0.005 -0.004  0.001 -0.283               cropforage -0.161 -0.005  0.012 -0.004  0.791 -0.231        cropsoy    -0.175 -0.022  0.013  0.013  0.404 -0.164  0.557   All of my continuous variables (denoted by a small s before the variable name) are standardized (z-scores). season is a categorical variable with 2 levels (early and late), and crop is a categorical variable with 3 levels (corn, forage, and soy).  This correlation of fixed effects matrix is really confusing me, because all of the correlations have the opposite sign that they do when I look at the simple regressions of pairs of variables. i.e., the correlation of fixed effects matrix suggests a strong positive correlation between cropforage and sbare, when in fact there is a very strong NEGATIVE correlation between these variables - forage crops tended to have much less bare ground compared to corn and soy crops. Pairs of continuous variables have the same issue, the correlation of fixed effects matrix says everything is the opposite of what it should be ... Could this just be due to the complexity of the model (not being a simple regression)? Could it have something to do with the fact that the variables are standardized?   Thanks.
72652
When explaining LASSO regression, the diagram of a diamond and circle is often used. It is said that because the shape of the constraint in LASSO is a diamond, the least squares solution obtained might touch the corner of the diamond such that it leads to a shrinkage of some variable. However, in ridge regression, because it is a circle, it will often not touch the axis. I could not understand why it cannot touch the axis or maybe have a lower probability than LASSO to shrink certain parameters. On top of that, why do LASSO and ridge have lower variance than ordinary least squares? The above is my understanding of ridge and LASSO and I might be wrong. Can someone help me understand why these two regression methods have lower variance?
83480
My data is described here   I am trying to see the effect of an interaction using lmer so my base case is:  my_null.model &lt;- lmer(value ~ Condition+Scenario+                               (1|Player)+(1|Trial), data = my, REML=FALSE)  my.model &lt;- lmer(value ~ Condition*Scenario+                               (1|Player)+(1|Trial), data = my, REML=FALSE)   Running the anova gives me significant results, but when I try to account for random slope ((1+Scenario|Player)) the model fails with this error:     Warning messages:  1: In commonArgs(par, fn, control, environment()) :    maxfun &lt; 10 * length(par)^2 is not recommended.  2: In optwrap(optimizer, devfun, getStart(start, rho$lower, rho$pp),  :   convergence code 1 from bobyqa: bobyqa -- maximum number of function evaluations exceeded  3: In commonArgs(par, fn, control, environment()) :   maxfun &lt; 10 * length(par)^2 is not recommended.  4: In optwrap(optimizer, devfun, opt$par, lower = rho$lower, control = control,  :    convergence code 1 from bobyqa: bobyqa -- maximum number of function evaluations exceeded  5: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :    Model failed to converge with max|grad| = 36.9306 (tol = 0.002)  6: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :    Model failed to converge: degenerate  Hessian with 1 negative eigenvalues   Alternatively if it fails to converge after a lot of iterations (I set it to 100 000) and I am getting the same results after 50k and 100k it means that it is very close to the actual value, just it does not reach it. So can I report my results like this?  Note that when I set the iterations so high I get only these warnings:   Warning messages: 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :  Model failed to converge with max|grad| = 43.4951 (tol = 0.002) 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :  Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
105690
I have two different implementations of ridge in MATLAB. One is simply    $\mathbf x = (\mathbf{A}'\mathbf{A}+\mathbf{I}\lambda)^{-1}\mathbf{A}'\mathbf b$  (as seen on ), with $\mathbf{I}$ being the identity matrix of size columns($\mathbf{A}$) $\times$ columns($\mathbf{A}$), and  I'm simply calling Matlab's "ridge" with   x = ridge(A, b, lambda)    My problem is that both return different results. (1) returns the results that I want (I know this by comparing results with other people) but why does (2) not return the same results?  My matrix $\mathbf A$ is sparse, it's filled with 1% 1's and 99% 0's. Some columns contain almost no 1's. The biggest difference seems to be that the coefficient for those columns with very few 1's are very close to 0 in (1), but can be quite far from 0 in (2)  Does anyone have any idea why it's different and how I can modify the call in (2) to give the same results as (1)?
5580
I would like to test the difference in response of two variables to one predictor. Here is a minimal reproducible example.   library(nlme)  ## gls is used in the application; lm would suffice for this example m.set &lt;- gls(Sepal.Length ~ Petal.Width, data = iris,                 subset = Species == "setosa") m.vir &lt;- gls(Sepal.Length ~ Petal.Width, data = iris,                 subset = Species == "virginica") m.ver &lt;- gls(Sepal.Length ~ Petal.Width, data = iris,                 subset = Species == "versicolor")   I can see that the slope coefficients are different:  m.set$coefficients (Intercept) Petal.Width    4.7771775   0.9301727 m.vir$coefficients (Intercept) Petal.Width    5.2694172   0.6508306  m.ver$coefficients (Intercept) Petal.Width     4.044640    1.426365    I have three questions:   How can I test the difference between slopes? How can I test the difference between residual variances? What is a simple, effective way to present these comparisons?    A related question, , suggests re-running the model with a dummy variable to differentiate the slopes, are there options that would allow the use of independent data sets?
78291
I'm trying to fit generalized linear models to some sets of count data that might or might not be overdispersed. The two canonical distributions that apply here are the Poisson and Negative Binomial (Negbin), with E.V. $\mu$ and variance   $Var_P = \mu $  $Var_{NB} = \mu + \frac{\mu^2}{\theta} $  which can be fitted in R using glm(..,family=poisson) and glm.nb(...), respectively. There is also the quasipoisson family, which in my understanding is an adjusted Poisson with the same E.V. and variance  $Var_{QP} = \phi\mu$,  i.e. falling somewhere in between Poisson and Negbin. The main problem with the quasipoisson family is that there is no corresponding likelihood for it, and hence a lot of extremely useful statistical tests and fit measures (AIC, LR etcetera) are unavailable.   If you compare the QP and Negbin variances, you might notice that you could equate them by putting $\phi = 1 + \frac{\mu}{\theta}$. Continuing on this logic, you could try to express the quasipoisson distribution as a special case of the Negbin:  $QP\,(\mu,\phi) = NB\,(\mu,\theta = \frac{\mu}{\phi-1})$,  i.e. a Negbin with $\theta$ linearly dependent on $\mu$. I tried to verify this idea by generating a random sequence of numbers according to the above formula and fitting it with glm:  #fix parameters  phi = 3 a = 1/50 b = 3 x = 1:100  #generating points according to an exp-linear curve #this way the default log-link recovers the same parameters for comparison  mu = exp(a*x+b)  y = rnbinom(n = length(mu), mu = mu, size = mu/(phi-1)) #random negbin generator  #fit a generalized linear model y = f(x)   glmQP = glm(y~x, family=quasipoisson) #quasipoisson glmNB = glm.nb(y~x) #negative binomial  &gt; glmQP  Call:  glm(formula = y ~ x, family = quasipoisson)  Coefficients: (Intercept)            x       3.11257      0.01854   (Dispersion parameter for quasipoisson family taken to be 3.613573)  Degrees of Freedom: 99 Total (i.e. Null);  98 Residual Null Deviance:      2097  Residual Deviance: 356.8    AIC: NA  &gt; glmNB  Call:  glm.nb(formula = y ~ x, init.theta = 23.36389741, link = log)  Coefficients: (Intercept)            x       3.10182      0.01873    Degrees of Freedom: 99 Total (i.e. Null);  98 Residual Null Deviance:      578.1  Residual Deviance: 107.8    AIC: 824.7   Both fits reproduce the parameters, and the quasipoisson gives a 'reasonable' estimate for $\phi$. We can now also define an AIC value for the quasipoisson:  df = 3 # three model parameters: a,b, and phi phi.fit = 3.613573 #fitted phi value copied from summary(glmQP) mu.fit = glmQP$fitted.values   #dnbinom = negbin density, log=T returns log probabilities AIC = 2*df - 2*sum(dnbinom(y, mu=mu.fit, size = mu.fit/(phi.fit - 1), log=T)) &gt; AIC [1] 819.329   (I had to manually copy the fitted $\phi$ value from summary(glmQP), as I couldn't find it in the glmQP object)  Since $AIC_{QP} &lt; AIC_{NB}$, this would indicate that the quasipoisson is, unsurprisingly, the better fit; so at least $AIC_{QP}$ does what it should do, and hence it might be a reasonable definition for the AIC (and by extension, likelihood) of a quasipoisson. The big questions I am left with are then   Does this idea make sense? Is my verification based on circular reasoning? The main question for anyone that 'invents' something that seems to be missing from a well established topic: if this idea makes sense, why isn't it already implemented in glm ?   Edit: figure added
32113
We have comprehensive coverage on the topic of perfect separation in logistic regression. When it happens in R we usually see two warnings:     Warning messages:   1: glm.fit: algorithm did not converge   2: glm.fit: fitted probabilities numerically 0 or 1 occurred    Is it possible that the algorithm did converge and fitted probabilities numerically 0 or 1 occurred nonetheless?
43059
I hope it's okay to ask theoretically driven R questions here.  R has given me the following results from my 'tournament of models'.  All models are entirely distinct except from 3 basic control variables. It's clear to me that Model 3 (third column) is the best performing. My concern is that Model 4 has such as low $R^2$, despite having a low AIC. How would you theoretically explain why a model has low and AIC but also a low $R^2$?  Secondly, imagine I went to add the variables from Model 4 into each of the other models and note the reduction of AIC scores. If the reductions in AIC is different across each model, what does this suggest?  Should it be the case that the inserted variables do better when inserted into a very different thereotical model?    [For info I consulted this other thread but decided to create my own since I was still unclear of the distinction between goodness of fit and explanatory power.]
27355
Let's take the following example:  set.seed(342) x1 &lt;- runif(100) x2 &lt;- runif(100) y &lt;- x1+x2 + 2*x1*x2 + rnorm(100) fit &lt;- lm(y~x1*x2)   This creates a model of y based on x1 and x2, using a OLS regression.  If we wish to predict y for a given x_vec we could simply use the formula we get from the summary(fit).  However, what if we want to predict the lower and upper predictions of y? (for a given confidence level).  How then would we build the formula?
38410
In almost all of the analysis work that I've ever done I use:   set.seed(42)    It's an homage to Hitchhiker's Guide to the Galaxy. But I'm wondering if I'm creating bias by using the same seed over and over.
36352
I'm studying LASSO regression, in particular the choice of the optimal tuning parameter.  The glmnet package and the book "Elements of Statistical Learning" offer two possible tuning Parameters: The $\lambda$, that minimizes the average error, and the $\lambda$, selected by the "one-standard-error" rule.  I couldn't find any research or evidence on    the "one-standard-error" rule, and which $\lambda$ I should use for my LASSO-regression.   EDIT: From the textbook elements of statistical learning   "Often a “one-standard error” rule is used with cross-validation, in which we choose the most par- simonious model whose error is no more than one standard error above the error of the best model."
87211
Following is the plot of glmnet with default alpha(1, hence lasso) using mtcars data set in R with mpg as the DV and others as the predictor variables.  glmnet(as.matrix(mtcars[-1]), mtcars[,1])     What can we conclude from this plot regarding different variables, especially am, cyl and wt (red, black and light blue lines)? How would we phrase the output in a report to be published?  I thought of following:   wt is the most important predictor of mpg. It is negatively affecting the mpg. cyl is a weak negative predictor of mpg. am may be a positive predictor of mpg. Other variables are not robust predictors of mpg.    Thanks for your thoughts on this.  (Note: cyl is the black line which does not reach 0 till very close to it.)  Edit: Following is plot(mod, xvar='lambda') which shows x-axis in order reverse of above plot: