Lasso Regression – Resolving Null Returns in GAMLasso from PLSMSelect Package

generalized-additive-modellasso

I was fitting a generalized additive model (GAM) with LASSO penalty to my data using gamlasso, as in the following (here I use simulated data, but the idea is similar):

set.seed(1)
Sex <- sample(c(rep(1, 40), rep(0, 60)))
set.seed(2)
Group <- sample(c(rep(1, 30), rep(0, 70)))
set.seed(3)
Ed <- sample(c(rep(1, 20), rep(0, 80)))
Age <- sample(c(rep(7, 15), rep(8, 20), rep(9, 30), rep(10, 10), rep(11, 15), rep(12, 10)))

sample_data <-data.frame(Age = Age,
                         Sex = Sex,
                         Ed = Ed,
                         Group = Group,
                         Perf1 = runif(100, -200, 100),
                         Perf2 = runif(100, 20, 200),
                         Perf3 = runif(100, 0.03, 0.3)
                         )

sample_data$x <- model.matrix( ~ Age + Sex + Ed,
                        data=sample_data)[,-1]

gamlasso_fit <- gamlasso(Group ~ x + s(Perf1) + s(Perf2) + s(Perf3),
                         data = sample_data, 
                         smooth.penalty = "l1", 
                         linear.penalty = "l1", 
                         family = "binomial")

gamlasso_fit$cv.glmnet
summary(gamlasso_fit, s = lambda.1se)

Which gives me:

> gamlasso_fit$cv.glmnet
NULL
> summary(gamlasso_fit, s = lambda.1se)
$lasso
NULL

$gam

Family: binomial 
Link function: logit 

Formula:
Group ~ x + s(Perf1) + s(Perf2) + s(Perf3)

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept) -1.88730    1.63831  -1.152   0.2493  
xAge         0.09693    0.17289   0.561   0.5750  
xSex         0.71345    0.53603   1.331   0.1832  
xEd         -1.71468    0.88280  -1.942   0.0521 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
           edf Ref.df Chi.sq p-value
s(Perf1) 4.381  5.363  4.858   0.494
s(Perf2) 1.000  1.000  0.009   0.926
s(Perf3) 4.734  5.740  5.770   0.395

R-sq.(adj) =   0.12   Deviance explained = 19.4%
UBRE = 0.26659  Scale est. = 1         n = 100

My actual model is giving me the same kind of result with the lasso part ($lasso) missing in the output (see summary and gamlasso_fit$cv.glmnet), which should be part of the output due to my model specification. My question is whether I am specifying anything wrong here. I followed the vignette of the plsmselect package and their examples, however I think I might be missing something here. I am using a GAM (binomial) with LASSO penalty since non-linearity is found between the log-odds of the outcome (Group) and some of the predictors in a glm (with/without LASSO penalty).

Best Answer

You should change all instances of x (lowercase) to X (uppercase). So your code will be

set.seed(1)
Sex <- sample(c(rep(1, 40), rep(0, 60)))
set.seed(2)
Group <- sample(c(rep(1, 30), rep(0, 70)))
set.seed(3)
Ed <- sample(c(rep(1, 20), rep(0, 80)))
Age <- sample(c(rep(7, 15), rep(8, 20), rep(9, 30), rep(10, 10), rep(11, 15), rep(12, 10)))

sample_data <-data.frame(Age = Age,
                         Sex = Sex,
                         Ed = Ed,
                         Group = Group,
                         Perf1 = runif(100, -200, 100),
                         Perf2 = runif(100, 20, 200),
                         Perf3 = runif(100, 0.03, 0.3)
                         )

sample_data$X <- model.matrix( ~ Age + Sex + Ed,
                        data=sample_data)[,-1]

gamlasso_fit <- gamlasso(Group ~ X + s(Perf1) + s(Perf2) + s(Perf3),
                         data = sample_data, 
                         smooth.penalty = "l1", 
                         linear.penalty = "l1", 
                         family = "binomial")

is.null(gamlasso_fit$cv.glmnet) # FALSE

It isn't explicitly mentioned in the documentation to use an uppercase X but the package is actually hardcoded like that. In the next version of the package we will try to make it more flexible and add clearer documentation.

Related Solutions

Solved – Using LASSO from lars (or glmnet) package in R for variable selection

Using glmnet is really easy once you get the grasp of it thanks to its excellent vignette in http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html (you can also check the CRAN package page). As for the best lambda for glmnet, the rule of thumb is to use

cvfit <- glmnet::cv.glmnet(x, y)
coef(cvfit, s = "lambda.1se")

instead of lambda.min.

To do the same for lars you have to do it by hand. Here is my solution

cv <- lars::cv.lars(x, y, plot.it = FALSE, mode = "step")
idx <- which.max(cv$cv - cv$cv.error <= min(cv$cv))
coef(lars::lars(x, y))[idx,]

Bear in mind that this is not exactly the same, because this is stopping at a lasso knot (when a variable enters) instead of at any point.

Please note that glmnet is the preferred package now, it is actively maintained, more so than lars, and that there have been questions about glmnet vs lars answered before (algorithms used differ).

As for your question of using lasso to choose variables and then fit OLS, it is an ongoing debate. Google for OLS post Lasso and there are some papers discussing the topic. Even the authors of Elements of Statistical Learning admit it is possible.

Edit: Here is the code to reproduce more accurately what glmnet does in lars

  cv <- lars::cv.lars(x, y, plot.it = FALSE)
  ideal_l1_ratio <- cv$index[which.max(cv$cv - cv$cv.error <= min(cv$cv))]
  obj <- lars::lars(x, y)
  scaled_coefs <- scale(obj$beta, FALSE, 1 / obj$normx)
  l1 <- apply(X = scaled_coefs, MARGIN = 1, FUN = function(x) sum(abs(x)))
  coef(obj)[which.max(l1 / tail(l1, 1) > ideal_l1_ratio),]

Solved – implement lasso Regression in GAMs in R

It's not the Lasso, but it has the same effect; Marra & Wood (2011) proposed two options for this and describe several other ways to do feature selection in GAMs.

Both proposals involve modifying the usual ridge/wiggliness penalty(ies), which do not affect the perfectly smooth basis functions (the null space of the basis), such that the model can penalise the perfectly smooth parts of the basis and thus shrink terms out of the model entirely.

The two proposed options are:

For thin plate and cubic regression splines, Marra & Wood (2011) suggest using an eigendecomposition of the penalty matrix, where the perfectly smooth basis functions have an eigenvalue equal to zero (which is why they don't contribute to the usual wiggliness penalty). They propose that a small value be added to the zero eigenvalues such that they now are also affected by the wiggliness penalty.

The small value that is added implies that the wiggly parts of the basis will get shrunk more than the perfectly smooth parts.

In mgcv you can choose to use these types of smoother via the bs argument to s(), te() etc. For shrinkage thin plate splines use bs = 'ts' and for the cubic spline equivalent use bs = 'cs'.
The second option, the double penalty method and which has somewhat better performance, is to impose a second penalty that affect only the perfectly smooth terms. Now there are two penalties, one controlling the wiggliness and one controlling the smooth parts of the basis, which means that the linear terms in the model can be shrunk from the model also.

This implies that you want to shrink the wiggly and smooth parts of the basis in the same way.

You add this to the model in mgcv by using the select = TRUE argument to gam() etc. This means that you have to turn this method on for all smooth terms in the model, whereas the shrinkage smooth option above allows you to selectively decide which terms should be subject to feature selection. The double penalty approach also means that two smoothness parameters per smooth (instead of the usual one) need to be estimated/selected during fitting.

Marra, G., Wood, S.N., 2011. Practical variable selection for generalized additive models. Comput. Stat. Data Anal. 55, 2372–2387. https://doi.org/10.1016/j.csda.2011.02.004

Best Answer

Related Solutions

Solved – Using LASSO from lars (or glmnet) package in R for variable selection

Solved – implement lasso Regression in GAMs in R

Related Question