8
$\begingroup$

I'm attempting to predict vegetation productivity based on climatic and land use variables (the latter are categorical). I found that there is a multicollinearity problem between the predictors (especially land use) as seen from the Variance Inflation Factor (VIF of the Ordinary Least Squares Regression).

Although my knowledge of lasso regression is basic, I assume lasso regression might solve the multicollinearity problem and also select variables that are driving the system. I appreciate an R code for estimating the standardized beta coefficients for the predictors or approaches on how to proceed.

Variable           Coeff.  Std Coeff.  VIF    Std Error    t      P  Value 
Constant          -0.228   0            0      0.086       -2.644  0.008  
Precipitation      <.001   0.151       2.688   <.001        8.541  0.0  
Solar Rad          0.002   0.343       2.836   <.001        18.939 <.001  
Temp              -0.116  -1.604       28.12   0.004       -28.11  0.0  
Water Stress       0.881   0.391       2.352   0.037        23.7   <.001  
Vapor Pressure     0.135   1.382       30.49   0.006        23.259 0.0    
  1               -0.103   -0.109      52.086  0.074       -1.398  0.162    
  2               -0.14    -0.048      6.49    0.079       -1.761  0.078   
  3               -0.11    -0.048      10.007  0.077       -1.42   0.156    
  4               -0.104   -0.234      236.288 0.073       -1.416  0.157    
  5               -0.097   -0.242      285.244 0.073       -1.331  0.183    
  6               -0.104   -0.09       35.067  0.074       -1.406  0.16    
  8               -0.119   -0.261      221.361 0.073       -1.629  0.103 
ELEVATION          <.001   -0.115      3.917   <.001       -5.381  <.001
Condition Number: 59.833 
Mean of Correlation Matrix: 0.221 1st    
Eigenvalue divided by m: 0.328
Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ You might want to look into the package glmnet. The problem with lasso, however, is that some of your categorical predictors might be left out while others might be retained. Are you okay with that? $\endgroup$ Commented Dec 30, 2015 at 19:16
  • $\begingroup$ @JohnK, I'm ok with some categorical predictors being left out. But i am new to R, and i was looking for a sample R code for glmnet. $\endgroup$ Commented Dec 30, 2015 at 19:20
  • $\begingroup$ Sure, I'll write a toy example in R. $\endgroup$ Commented Dec 30, 2015 at 19:24
  • 2
    $\begingroup$ Note that, regardless of some of the answers, requests for R code are off-topic here. But the longer answers pick up on your implication that you are interested in statistical advice and that content keeps the thread within forum guidelines. For any future questions please study advice in the Help Center on software-specific questions. (I edited the title to remove reference to R code.) $\endgroup$ Commented Dec 31, 2015 at 2:33

3 Answers 3

Reset to default
16
$\begingroup$

Be sure to install and load the glmnet package.

install.packages("glmnet")
library(glmnet)

First you need to form a matrix with all your predictors, we call that matrix $\mathbf{X}$. I have done this for three variables I have created but since you have more you will need to change the predictor matrix accordingly.

    set.seed(1)
    x1 <- rnorm(30)
    x2 <- rnorm(30)
    x3 <- rnorm(30)
    X <- matrix( c(x1, x2, x3), byrow = F, ncol = 3)

Then we need a response as well. This needs to be a vector as you know so let's form a linear combination of the predictors and corrupt it with some noise.

y <- 3 + 4*x1 + 3*x2 + 5*x3 + rnorm(30)

We have everything we need now to begin. Here is our first attempt with glmnet.

    fit <-glmnet(x = X, y = y, alpha = 1) 
# different values of alpha return different estimators, alpha = 1 is the lasso.
    plot(fit, xvar = "lambda")

which should produce

enter image description here

To interpet this plot, recall the optimization problem the lasso solves

$$\hat{\boldsymbol{\beta}}_{LASSO} = \min_{\beta} \left( \mathbf{y} - \mathbf{X} \boldsymbol{\beta} \right) ^{\prime} \left( \mathbf{y} - \mathbf{X} \boldsymbol{\beta} \right) + \lambda \sum_j \left|\beta_j \right| $$

so $\lambda$ is the penalty or the Lagrange multiplier if you prefer and is always positive. Setting $\lambda = 0 $ yields the familiar minimization of squared residuals while for greater values, some of the coefficients will be set to zero. As $\lambda \to \infty$, all the coefficents will be set to zero.

This is exactly what this plot shows then, the coefficient path for different values of lambda. For reasons beyond me the creators of this package have opted to present lambda on the log scale, thus values between zero and one are now negative.

Pick the value of lambda you like and you can extract the coefficients with the command

coef(fit, s = 0.3) # s is the value of lambda

To locate the point on the plot, simply do

log(0.3)

and as you can see since lambda is quite close to zero there is not much shrinkage. Of course now we have to select one of these values of $\lambda$ and visual inspection is not good enough.

We do this by the crossvalidation function of glmnet. If you have never heard of crossvalidation, all you need to know is that it is a predictive criterion that evaluates the sample performance by splitting the sample into training and validation sets and choosing the value of lambda with which the error of prediction is minimal.

crossval <-  cv.glmnet(x = X, y = y)
plot(crossval)
penalty <- crossval$lambda.min #optimal lambda
penalty #minimal shrinkage
fit1 <-glmnet(x = X, y = y, alpha = 1, lambda = penalty ) #estimate the model with that
coef(fit1)

I have chosen to plot the crossvalidation results just so you can see how this method works. We finally estimate the lasso with the optimal CV parameter and extract the coefficients. You will notice in my example that the shrinkage is minimal. This occurs because there is no multicollinearity in the sample thanks to the naive generation. I am certain your results will be much much different.

Here is a webpage where the creators explain in detail how to use the package. I kept it simple because you are only interested in the lasso and no other estimators.

Hope this helps.

Improve this answer
$\endgroup$
3
  • $\begingroup$ This is greatly helpful. I have the Matrix X imported from Excel as: mydata <- read.csv("2377Dummy_Popds.csv"). Where the csv file consists of all the predictors and the response variable. How do i just import the both the response and predictor as X ( or mydata matrix) and apply glmnet? i.e. the response variable is one column in the mydata matrix) $\endgroup$ Commented Dec 30, 2015 at 20:45
  • $\begingroup$ @Kaleab If I were you, I would just import the entire file first and then select the appropriate columns, i.e. variables manually. If you google "import csv in R" a million results will come up, it's quite straightforward. $\endgroup$ Commented Dec 30, 2015 at 20:59
  • 1
    $\begingroup$ @JohnK This is so helpful! I have been searching for quite some time to find a "beginner's guide" to model selection via LASSO and this is by far the best one I've come across. Thank you! $\endgroup$ Commented Aug 9, 2017 at 22:41
2
$\begingroup$

You might be better off using ridge regression rather than LASSO, or devising a scale that combines some of the correlated variables in an intelligent way, based on your knowledge of the subject matter. LASSO will not select the variables that are "driving the system," just a subset that happens to work on this data sample and which could be a good deal different from the variables selected on a different sample.

Once you get glmnet working, it will be trivial to try ridge regression instead of LASSO. If your interest is in prediction then including all variables with ridge regression will tend to work better, particularly with so few predictors.

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ Lasso might not select the variables of interest but ridge will not select anything, just shrink. $\endgroup$ Commented Dec 31, 2015 at 0:16
  • $\begingroup$ @JohnK the problem is that the variables LASSO will select have no reason to be the variables "driving the system" as the OP implies. Try doing LASSO on multiple bootstraps of the same sample to see how variable the selection of variables can be, particularly with severe multicollinearity. $\endgroup$ Commented Dec 31, 2015 at 0:28
  • 1
    $\begingroup$ I agree that lasso is often sub-optimal but I gather that the OP is interested in performing variable selection so it's either the lasso or the elastic net which is a compromise between ridge and lasso. The point is that ridge will keep all variables. Anyway +1 for the recommendation. I believe in ridge regression. $\endgroup$ Commented Dec 31, 2015 at 0:32
  • $\begingroup$ @EdM the property u describe about variable selection 'could be a good deal different' given a new sample - is that why Lasso is called an unstable selector or does that mean something else? $\endgroup$ Commented Oct 20, 2021 at 1:48
  • 1
    $\begingroup$ @JohnVandivier yes. When there are multiple correlated predictors, the particular one(s) selected by LASSO can vary from sample to sample, so in that sense its predictor selection is “unstable.” Any of the selections might work OK for prediction, nevertheless. So although there can be “significance tests” for predictors selected by LASSO, the meaning of that “significance” might not be what you might naively expect. There is no assurance that the selected predictors are “driving the system” as the OP might hope. $\endgroup$ Commented Oct 20, 2021 at 2:23
0
$\begingroup$

The one thing that I hate about glmnet is that it doesn't use data frames as an argument. In order to use glmnet, a data frame needs to be transformed into a model matrix and a response vector. Most model functions, like lm and glm, allow for formulas and data frames.

The caret package is a wrapper for many modelling packages. Here's some sample code that uses caret to fit a glmnet model. Note that elastic net is a generalization of ridge and lasso regression.

https://stats.stackexchange.com/a/188780/24030

Improve this answer
$\endgroup$
1
  • 2
    $\begingroup$ It's quite easy to write a wrapper for glmnet that accepts a formula. A tiny investment can save you the hate. $\endgroup$ Commented Dec 31, 2015 at 6:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.