Ridge and Lasso Exercises in R: 18 Practice Problems with Solutions

ex_1_1 <- model.matrix(medv ~ ., Boston)[, -1]
dim(ex_1_1)
#> [1] 506  13
head(ex_1_1[, 1:5], 3)
#>      crim zn indus chas   nox
#> 1 0.00632 18  2.31    0 0.538
#> 2 0.02731  0  7.07    0 0.469
#> 3 0.02729  0  7.07    0 0.469

Explanation: model.matrix() does the dummy-coding work even when all your predictors are numeric, so it is the safest converter. The [, -1] drop is critical because glmnet adds its own intercept and a duplicate constant column would be flagged as singular. Once you have this matrix, every later glmnet() and cv.glmnet() call reuses it.

x <- model.matrix(medv ~ ., Boston)[, -1]
y <- Boston$medv
set.seed(1)
ex_1_2 <- glmnet(x, y, alpha = 0)
class(ex_1_2)
#> [1] "elnet"  "glmnet"
dim(coef(ex_1_2))
#> [1]  14 100

Explanation: alpha = 0 selects the L2 penalty, which is ridge. glmnet automatically picks a sequence of 100 lambdas (returned as ex_1_2$lambda) ranging from one large enough to zero every coefficient down to the data-driven minimum. The 14 coefficient rows are the 13 predictors plus the intercept. Set a seed even though the path itself is deterministic; later cross-validation steps depend on a stable RNG state.

set.seed(1)
ex_1_3 <- glmnet(x, y, alpha = 1)
ex_1_3$lambda[1:3]
#> [1] 6.7849 6.1822 5.6331
sum(coef(ex_1_3, s = 0.5)[-1, ] == 0)
#> [1] 4

Explanation: The L1 penalty is non-differentiable at zero, which is exactly why coordinate descent can pin individual coefficients there. The [-1, ] slice strips the intercept row before counting zeros. At s = 0.5 four predictors out of thirteen are removed; raise the lambda and more drop out until eventually only the intercept remains.

plot(ex_1_3, xvar = "lambda", label = TRUE)
ex_2_1 <- ex_1_3
head(ex_2_1$beta[, 80:82], 3)
#>            s79         s80         s81
#> crim -0.105014 -0.105657 -0.106246
#> zn    0.046194  0.046408  0.046604
#> indus 0.000000  0.000000  0.000000

Explanation: The xvar = "lambda" argument flips the x-axis to log(lambda), which is the most readable scale because lambdas span orders of magnitude. The labels at the right edge identify each variable's coefficient at the smallest penalty. Variables that hug zero across the whole sweep are noise candidates; the survivors at high lambda are the strongest signals. For lasso on Boston, lstat (low-status share of population) and rm (rooms per dwelling) dominate.

ex_2_2 <- as.numeric(coef(ex_1_3, s = 0.1))[-1]
names(ex_2_2) <- rownames(coef(ex_1_3, s = 0.1))[-1]
length(ex_2_2)
#> [1] 13
round(ex_2_2[1:5], 4)
#>     crim       zn    indus     chas      nox
#>  -0.0966   0.0408   0.0000   2.6837 -16.6443

Explanation: coef.glmnet returns a column of a sparse dgCMatrix, which is memory-efficient but awkward to slice. as.numeric() materialises it to a dense numeric vector, then [-1] removes the intercept. Re-attaching names from the sparse object preserves the predictor labels. Always supply an explicit s value; the default is the entire lambda sequence and yields a 14-by-100 matrix.

ridge_b <- as.numeric(coef(ex_1_2, s = 0.5))[-1]
lasso_b <- as.numeric(coef(ex_1_3, s = 0.5))[-1]
ex_2_3 <- tibble::tibble(
  predictor = colnames(x),
  ridge     = round(ridge_b, 2),
  lasso     = round(lasso_b, 2)
) |> arrange(desc(abs(lasso)))
ex_2_3
#> # A tibble: 13 x 3
#>    predictor   ridge   lasso
#>    <chr>       <dbl>   <dbl>
#>  1 nox       -8.32   -3.11  
#> ...

Explanation: At identical lambda the L2 ridge keeps every coefficient non-zero but pushes them toward each other, while the L1 lasso pins the weakest signals to exactly zero (tax, indus, age). The interesting cell is nox: ridge keeps a large negative coefficient, lasso shrinks it sharply because it shares signal with dis and rad. This kind of side-by-side is the most readable way to show stakeholders the practical difference between the two penalties.

set.seed(42)
ex_3_1 <- cv.glmnet(x, y, alpha = 1, nfolds = 10)
ex_3_1$lambda.min
#> [1] 0.02212
ex_3_1$lambda.1se
#> [1] 0.4339
ex_3_1$nzero[match(c(ex_3_1$lambda.min, ex_3_1$lambda.1se), ex_3_1$lambda)]
#> s52 s17
#>  12   8

Explanation: cv.glmnet runs the full path under k-fold CV and returns two principled lambdas. lambda.min is the empirical winner on average held-out MSE; lambda.1se is the largest lambda whose CV error sits within one standard error of the minimum, giving you a sparser, more conservative model. Reporting both is standard practice; choosing lambda.1se for production deployments is a good default when interpretability matters.

set.seed(99)
train_idx <- sample(nrow(Boston), 0.7 * nrow(Boston))
x_tr <- x[train_idx, ];  y_tr <- y[train_idx]
x_te <- x[-train_idx, ]; y_te <- y[-train_idx]

set.seed(42)
cvfit <- cv.glmnet(x_tr, y_tr, alpha = 1, nfolds = 10)
pred_min <- predict(cvfit, newx = x_te, s = "lambda.min")
pred_1se <- predict(cvfit, newx = x_te, s = "lambda.1se")

ex_3_2 <- tibble::tibble(
  lambda_choice = c("lambda.min", "lambda.1se"),
  mse = c(mean((y_te - pred_min)^2), mean((y_te - pred_1se)^2))
)
ex_3_2
#> # A tibble: 2 x 2
#>   lambda_choice    mse
#>   <chr>          <dbl>
#> 1 lambda.min      24.7
#> 2 lambda.1se      26.4

Explanation: The holdout is unseen by every CV fold, so its MSE is an unbiased estimate of generalisation error. lambda.min usually wins on the holdout when the sample is large; lambda.1se wins more often when n is small or noisy because it absorbs less variance. Always rerun cv.glmnet on the train slice, never on the full set, or you leak label information into the lambda choice.

set.seed(7)
ex_3_3 <- cv.glmnet(x, y, alpha = 1, nfolds = 10, type.measure = "mae")
ex_3_3$name
#>          mae
#> "Mean Absolute Error"
ex_3_3$lambda.min
#> [1] 0.0221
ex_3_3$lambda.1se
#> [1] 0.418

Explanation: type.measure accepts "mse" (default for Gaussian), "mae", "deviance", "class" (for classification), "auc" (binary), and "C" (Cox). MAE is robust to outliers because it weights every error linearly instead of quadratically. The chosen lambdas tend to be similar to MSE-based ones on clean data but diverge sharply when the response has fat tails or extreme observations.

beta_1se <- coef(ex_3_1, s = "lambda.1se")[-1, ]
ex_4_1 <- sum(beta_1se != 0)
ex_4_1
#> [1] 8

Explanation: The [-1, ] again strips the intercept so it does not get counted as a predictor. Comparing strictly with != 0 is safe because lasso enforces exact zeros (not near-zero), unlike regularizers that just shrink. Eight survivors out of thirteen on Boston is a typical sparsity level: the model retains the strongest neighbourhood signals (rm, lstat, dis, ptratio) while pruning the redundant ones.

set.seed(123)
n <- 200; p <- 50
X_sim <- matrix(rnorm(n * p), n, p)
true_beta <- c(3, -2, 1.5, -1, 0.8, rep(0, 45))
y_sim <- X_sim %*% true_beta + rnorm(n)

set.seed(7)
cv_sim <- cv.glmnet(X_sim, y_sim, alpha = 1)
beta_min <- coef(cv_sim, s = "lambda.min")[-1, ]
ex_4_2 <- which(beta_min != 0)
ex_4_2
#> [1]  1  2  3  4  5
length(ex_4_2)
#> [1] 5

Explanation: This kind of recovery experiment is the cleanest way to see why lasso earns its name: with truly sparse ground truth and uncorrelated predictors, lasso reliably identifies the right variables once n exceeds a multiple of the active-set size. lambda.1se would prune more aggressively (sometimes dropping the smallest true signal); lambda.min is more permissive and often includes a noise variable or two. Repeat with different seeds to see the variability in selected sets.

beta_1se <- coef(ex_3_1, s = "lambda.1se")[-1, ]
selected <- names(beta_1se)[beta_1se != 0]
form <- as.formula(paste("medv ~", paste(selected, collapse = " + ")))
ex_4_3 <- lm(form, data = Boston)
coef(ex_4_3)
#> (Intercept)         crim         chas          nox           rm          dis      ptratio        black        lstat
#>     32.9890      -0.0786       2.7197     -16.6612       4.5263      -1.2864      -0.9534       0.0094      -0.5395

Explanation: The two-step "select with lasso, then refit with OLS" recipe is sometimes called the relaxed lasso (Meinshausen, 2007). It removes the shrinkage bias on the chosen coefficients while keeping lasso's feature selection. Be aware: standard errors from summary(ex_4_3) are too small because they ignore the random selection step. Use selectiveInference::fixedLassoInf or sample-splitting if you need honest p-values.

fit_alpha <- function(a) {
  set.seed(42)
  fit <- cv.glmnet(x, y, alpha = a, nfolds = 10)
  min(fit$cvm)
}
ex_5_1 <- tibble::tibble(
  alpha      = c(0, 0.5, 1),
  cv_mse_min = sapply(c(0, 0.5, 1), fit_alpha)
)
ex_5_1
#> # A tibble: 3 x 2
#>   alpha cv_mse_min
#>   <dbl>      <dbl>
#> 1   0         24.4
#> 2   0.5       23.7
#> 3   1         23.6

Explanation: On clean low-correlation data lasso usually beats ridge slightly because variable selection is genuinely useful. Elastic net at alpha = 0.5 lands between the two; on Boston the differences are small (under 1 MSE point on a response that ranges 5 to 50). Always reset the same seed before each cv.glmnet call so the fold assignments match across alphas; otherwise you cannot tell whether a difference reflects the penalty or the random folds.

alphas <- seq(0, 1, by = 0.1)
results <- lapply(alphas, function(a) {
  set.seed(42)
  fit <- cv.glmnet(x, y, alpha = a, nfolds = 10)
  data.frame(alpha = a, cv_mse_min = min(fit$cvm), lambda_min = fit$lambda.min)
})
ex_5_2 <- dplyr::bind_rows(results) |> tibble::as_tibble()
ex_5_2
#> # A tibble: 11 x 3
#>    alpha cv_mse_min lambda_min
#>    <dbl>      <dbl>      <dbl>
#>  1   0         24.4    0.692
#> ...
ex_5_2$alpha[which.min(ex_5_2$cv_mse_min)]
#> [1] 1

Explanation: A grid over alpha is the standard way to tune elastic net because there is no closed-form optimum for the mix. Eleven points is enough resolution for most problems; finer grids rarely change the winner. The caret::train or tidymodels::tune_grid wrappers automate the same loop with nested CV. On Boston the search confirms lasso, but on highly correlated designs the optimum often lands between 0.1 and 0.5.

set.seed(2)
n <- 200
ch1 <- rnorm(n)
ch2 <- ch1 + rnorm(n, sd = 0.1)
ch3 <- ch1 + rnorm(n, sd = 0.1)
noise <- matrix(rnorm(n * 2), n, 2)
Xc <- cbind(ch1, ch2, ch3, noise1 = noise[, 1], noise2 = noise[, 2])
yc <- 1 * ch1 + 1 * ch2 + 1 * ch3 + rnorm(n)

set.seed(42)
fit_l <- cv.glmnet(Xc, yc, alpha = 1)
fit_e <- cv.glmnet(Xc, yc, alpha = 0.3)
ex_5_3 <- tibble::tibble(
  predictor = colnames(Xc),
  lasso     = round(as.numeric(coef(fit_l, s = "lambda.min"))[-1], 2),
  enet      = round(as.numeric(coef(fit_e, s = "lambda.min"))[-1], 2)
)
ex_5_3
#> # A tibble: 5 x 3
#>   predictor   lasso   enet
#> 1 ch1         2.10   0.92
#> 2 ch2         0      0.85
#> 3 ch3         0      0.81

Explanation: When several predictors carry the same signal, lasso arbitrarily picks one and zeros the rest because the L1 penalty is indifferent between them. Elastic net's L2 component breaks that tie by spreading the coefficient mass across the correlated group, which is what marketers usually want when channels are causally distinct but statistically twinned. This grouping effect is the main reason to choose elastic net over pure lasso for attribution and genomics work.

set.seed(11)
train_idx <- sample(nrow(Boston), 0.8 * nrow(Boston))
x_tr <- x[train_idx, ];  y_tr <- y[train_idx]
x_te <- x[-train_idx, ]; y_te <- y[-train_idx]

set.seed(42)
cvfit <- cv.glmnet(x_tr, y_tr, alpha = 1)
pred  <- predict(cvfit, newx = x_te, s = "lambda.min")
ex_6_1 <- sqrt(mean((y_te - pred)^2))
ex_6_1
#> [1] 4.95

Explanation: RMSE is in the same units as the response (here, $1000s of median home value), which makes it the easiest metric to discuss with non-technical stakeholders. The predict.cv.glmnet method recognises the string "lambda.min" so you do not have to dig out the numeric value. For a deployable model you would refit on the full data at the chosen lambda after CV; for benchmarking against alternatives, the train/test split shown here is enough.

set.seed(11)
train_idx <- sample(nrow(Boston), 0.8 * nrow(Boston))
train_df <- Boston[train_idx, ]
test_df  <- Boston[-train_idx, ]

ols_fit  <- lm(medv ~ ., data = train_df)
ols_pred <- predict(ols_fit, newdata = test_df)
ols_rmse <- sqrt(mean((test_df$medv - ols_pred)^2))

ex_6_2 <- tibble::tibble(
  model = c("OLS", "Lasso CV"),
  rmse  = c(round(ols_rmse, 2), 4.95)
)
ex_6_2
#> # A tibble: 2 x 2
#>   model      rmse
#>   <chr>     <dbl>
#> 1 OLS        5.03
#> 2 Lasso CV   4.95

Explanation: On Boston the gap between OLS and lasso is small because n (506) is roughly 40 times p (13), so OLS is not particularly overfitted. The picture flips when p is comparable to or larger than n: ridge and lasso then routinely cut RMSE by 20% or more. The lesson is that regularisation is most valuable when the predictor count is high relative to the sample size, not as a default for every regression.

fit_default <- glmnet(x, y, alpha = 1, lambda = 0.1)
b1 <- as.numeric(coef(fit_default))[-1]

x_scaled <- scale(x)
fit_manual <- glmnet(x_scaled, y, alpha = 1, lambda = 0.1, standardize = FALSE)
b2_scaled <- as.numeric(coef(fit_manual))[-1]
b2 <- b2_scaled / attr(x_scaled, "scaled:scale")

ex_6_3 <- max(abs(b1 - b2))
ex_6_3
#> [1] 1.2e-13

Explanation: glmnet always solves the optimisation on standardised columns internally, then back-transforms the coefficients to the original scale before returning them. Doing the standardisation by hand and dividing by the per-column standard deviation reproduces the exact same coefficient vector to numerical precision. This matters when you are wiring glmnet into a pipeline that already standardises (recipes::step_normalize); set standardize = FALSE to avoid double-scaling.