# Why am I getting such poor results from this example of LASSO-regression?

I'm modelling the number of people in a population working, married, being ill, studying, having kids etc. I call one combination of these Boolean variables a "state".

I'm estimating transition probabilities between these states from one year to the next year. I want to use a logistic regression model to estimate these transition probabilities. Since there are a lot of variables in the data set, I want to try to choose variables automatically and I'm trying LASSO-regression at the moment.

Below is a made up example where I try to model the probability of working next year.

The usual glm-function in R works fine and gives estimations of the true parameters that are quite close (although the true influence of the variable married is actually 0). However, when I try LASSO-regression I get results that are far away from the true parameters. Furthermore, I thought that the LASSO-regression would pick up that the variable marriage does not have any influence, but keep the rest of the parameters. Why doesn’t it work?

library(glmnet)
library(boot)

seed <- 1

# Create the combinations of variables
set_working <- 0:1
set_married <- 0:1
set_age <- seq(from=20, to=60, by=1)

data_prob <- expand.grid(set_age, set_working,set_married)
colnames(data_prob) <- c("age","working","married")

# Create true underlying prob-model
intercept_parameter <- -3
age_parameter <- 1/20
working_parameter <- 4
married_parameter <- 0
true_parameters <- c(intercept_parameter,age_parameter,working_parameter,married_parameter)

data_prob$prob_working_next_year <- inv.logit(intercept_parameter + data_prob$age * age_parameter +
data_prob$working * working_parameter + data_prob$married * married_parameter)

# Generate random observations from the model
set.seed(seed)
data_prob$obs_working_next_year <- rbinom(n = nrow(data_prob), size = 1, prob=data_prob$prob_working_next_year)

# Fit glm-model using all the variables
fit_glm_all_variables <- glm(obs_working_next_year ~ age + working + married  ,
family = binomial,
data = data_prob)
# The result is reasonably close to the true parameters
coefficients(fit_glm_all_variables)

# Fit glm-model using all the variables and the LASSO-tecnique
x_matrix <- as.matrix(cbind(data_prob$age,data_prob$working,data_prob$married)) fit_glm_lasso <- cv.glmnet(x=x_matrix, y=data_prob$obs_working_next_year,
family = "binomial",
alpha=1)

# The result is far away from the true parameters
coefficients(fit_glm_lasso)
# 4 x 1 sparse Matrix of class "dgCMatrix"
# 1
# (Intercept) -0.5932176
# V1           .
# V2           2.6637460
# V3           .

# Get the whole sequence of parameters from using different lambda-values
fit_glm_lasso2 <- glmnet(x=x_matrix,
y=data_prob$obs_working_next_year, family = "binomial", alpha=1) coefficients(fit_glm_lasso2) # Get the whole sequence of parameters from using different lambda-values # This time explicately setting the glmnet-function to do the # standardization fit_glm_lasso3 <- glmnet(x=x_matrix, y=data_prob$obs_working_next_year,
family = "binomial",
alpha=1,
standardize = TRUE,
intercept = TRUE,
standardize.response = FALSE)
coefficients(fit_glm_lasso3)

# The "best" set of parameters chosen by the cv.glmnet-function
# When setting the glmnet-function to do the
# standardization
fit_glm_lasso4 <- cv.glmnet(x=x_matrix,
y=data_prob$obs_working_next_year, family = "binomial", alpha=1, standardize = TRUE, intercept = TRUE, standardize.response = FALSE) coefficients(fit_glm_lasso4) # 4 x 1 sparse Matrix of class "dgCMatrix" # 1 # (Intercept) -0.6273374 # V1 . # V2 2.7786131 # V3 . # >  • It's nice to not post code like rm(list = ls()) unless it's crucial for your example. No one wants to copy/paste your example and run that line accidentally. Commented Jul 19, 2017 at 16:09 • I'm unable to run your example at the moment, but a couple of clarifying questions: is there a particular reason that you are ending your investigation with cv.glmnet, and are using cv.glmnet on the whole data set? Are you certain that you want deviance (the default type.measure for cv.glmnet) rather than auc or classification error? What do you consider "quite close" in the coefficient estimates, and how much worse were the cv.glmnet coefficient estimates? Commented Jul 19, 2017 at 20:09 • @Gregor, you're right. I'll remove it. Commented Jul 20, 2017 at 6:45 • @Upper_Case, since this is a bit of a theoretical test I haven't divided the data in a traningset and a testingset. Since I decided the true parameters I figured the dividing of the dataset is not needed. Maybe in the real application though. I'm not certain that I want deviance, but why shouldn't the "presets" not work reasonably well in this case? The true values are -3, 0.05, 4 and 0. The coefficients(fit_glm_lasso) were -0.59, 0, 2.66, 0. Commented Jul 20, 2017 at 6:48 ## 1 Answer One thing is that your data isn't scaled; this means that age will have an artificially low coefficient compared to married. It is also why age is what drops out first, despite its true influence being higher than married. A good diagnostic is to plot the path for all lambda: fit_glm_lasso <- glmnet(x=x_matrix, y=data_prob$obs_working_next_year,
family = "binomial",
alpha=1)
plot(fit_glm_lasso)


where the scaling problem is apparent. If you scale your data first, then the lasso will drop out married first, as predicted:

x_scaled = scale(x_matrix)
fit_glm_lasso <- glmnet(x=x_scaled,
y=data_prob\$obs_working_next_year,
family = "binomial",
alpha=1)
plot(fit_glm_lasso)


The coefficients need to be reverse-transformed (not familiar with R enough to do this without a hacky solution), but they are ultimately much closer to your true model, especially if you take the point on the path where married drops out, before the optimum CV value is found.

• I didn't know that the x-values needed scaling in LASSO, but after reading a bit about it I found out that there are settings in glmnet that will do the scaling for me. I'm adding this in the code example above. There isn't that big of a difference though. But maybe I have misunderstood the cv.glmnet-function. I thought that it made the choice of lambda for me according to some criteria. But looking at the sequence of coefficents(fit_glm_lasso3) there are sets of estimates that are far better than the one that the cv.glmnet-function did chose. Commented Jul 20, 2017 at 6:54