I try to understand lasso regression. So far, I do understand that it can be used to shrink regression coefficients in case of few events. The coefficients of some covariates are even shrunk to zero. By shrinking coefficients, the issue of overfitting is addressed. That is how I understand this paper: BMJ. 2016 Jun 8;353:i3235.
I want to fit a logistic regression model to predict a future event. Unfortunately, available data is sparse and we have only 40 events. I thought to use lasso regression instead of stepwise backward selection this time. However, I do not know how to get odds ratios with respective 95% CIs for the covariates retained in the lasso regression model?
I found a working example on this web page, which I modified for my question: http://www.sthda.com/english/articles/36-classification-methods-essentials/149-penalized-logistic-regression-essentials-in-r-ridge-lasso-and-elastic-net/
Here is my example code:
library(dplyr)
library(caret)
library(glmnet)
set.seed(123)
#####################################
# Load the data and remove NAs
#####################################
data("PimaIndiansDiabetes2", package = "mlbench")
PimaIndiansDiabetes2 <- na.omit(PimaIndiansDiabetes2)
# Inspect the data
sample_n(PimaIndiansDiabetes2, 3)
# Split the data into training and test set
set.seed(123)
training.samples <- PimaIndiansDiabetes2$diabetes %>%
createDataPartition(p = 0.8, list = FALSE)
train.data <- PimaIndiansDiabetes2[training.samples, ]
test.data <- PimaIndiansDiabetes2[-training.samples, ]
# Dummy code categorical predictor variables
x <- model.matrix(diabetes~., train.data)[,-1]
x
head(x)
# Convert the outcome (class) to a numerical variable
y <- ifelse(train.data$diabetes == "pos", 1, 0)
#####################################
# Find the best lambda using cross-validation
#####################################
cv.lasso <- cv.glmnet(x, y, alpha = 1, family = "binomial")
cv.lasso
#####################################
# Final model with lambda.min
#####################################
lasso.model <- glmnet(x, y, alpha = 1, family = "binomial", lambda = cv.lasso$lambda.min)
lasso.model
coef(cv.lasso, cv.lasso$lambda.min)
#####################################
# Standard model
#####################################
mylogit <- glm(diabetes ~ pregnant + glucose + pressure + triceps + insulin + mass + pedigree + age, data = train.data, family = "binomial")
summary(mylogit)
# lasso standard
# (Intercept) -7.57530592 -9.503717
# pregnant 0.02027797 0.045710
# glucose 0.03350748 0.042303
# pressure . -0.007004
# triceps 0.01349090 0.018578
# insulin . -0.001592
# mass 0.02172361 0.045017
# pedigree 0.57491081 0.968452
# age 0.03302719 0.042557
#####################################
# Get OR plus 95% CI from standard model
#####################################
exp(cbind(OR = coef(mylogit), confint(mylogit)))