Three problems LASSO solves

1. Variable selection. LASSO shrinks unimportant coefficients to exactly zero, giving you a sparse model automatically. You do not need to specify which predictors to test.

2. Overfitting. By penalizing the total size of the coefficients, LASSO trades a small increase in bias for a large decrease in variance. The resulting model predicts new data better than OLS, especially when \(p\) (predictors) is large relative to \(n\) (observations).

3. Multicollinearity. When predictors are correlated (temperature and dissolved oxygen, nitrate and phosphate), OLS coefficients become unstable and inflated. LASSO stabilizes them by forcing the model to choose among correlated predictors rather than splitting the effect across all of them.

Why absolute value matters

The absolute value (\(|\beta_j|\)) is what makes LASSO special. Because the absolute value function has a corner at zero, the optimization can land exactly at \(\beta_j = 0\) for weak predictors. This is automatic variable selection. Ridge regression (Section 9) uses \(\beta_j^2\) instead, which shrinks coefficients toward zero but never reaches it. The geometric reason involves the shape of the constraint region – LASSO’s diamond has corners on the axes; ridge’s circle does not.

The Lagrangian form

The penalized objective is equivalent to a constrained optimization:

\[\min_\beta \; \text{RSS}(\beta) \quad \text{subject to} \quad \sum_{j=1}^{p} |\beta_j| \leq t\]

where \(t\) is a budget on the total coefficient size. Larger \(\lambda\) corresponds to smaller \(t\). The constraint boundary is a \(p\)-dimensional diamond (cross-polytope), and the OLS solution is a point that may lie inside or outside this diamond. If it lies outside, the LASSO solution is the point on the diamond boundary closest to the OLS solution. Because the diamond has corners on the coordinate axes, this point often has one or more coordinates exactly equal to zero.

Standardization

LASSO penalizes \(|\beta_j|\), so the scale of each predictor matters. A coefficient for temperature in Celsius will be numerically different from one for temperature in Fahrenheit, but the penalty treats them the same. To make the penalty fair, glmnet standardizes all predictors to have mean zero and unit variance before fitting. The reported coefficients are then back-transformed to the original scale.

Symptoms of overfitting

cat("Number of predictors:", length(coef(ols_full)) - 1, "\n")

#> Number of predictors: 14

cat("Number of observations:", nrow(wq), "\n")

#> Number of observations: 80

cat("Ratio n/p:", round(nrow(wq) / (length(coef(ols_full)) - 1), 1), "\n")

#> Ratio n/p: 5.7

cat("R-squared:", round(summary(ols_full)$r.squared, 3), "\n")

#> R-squared: 0.735

cat("Adjusted R-squared:", round(summary(ols_full)$adj.r.squared, 3), "\n")

#> Adjusted R-squared: 0.678

With 80 observations and 14 predictors, the ratio \(n/p \approx 5.7\). The model appears to fit well (high \(R^2\)), but much of this is overfitting. Several problems are visible in the summary:

Inflated standard errors. Correlated predictors (alkalinity and hardness, turbidity and TSS) inflate coefficient standard errors, making it hard to tell which predictors are truly important.

Spurious significance. With 14 predictors tested simultaneously, some will appear significant by chance alone (\(14 \times 0.05 = 0.7\) expected false positives at \(\alpha = 0.05\)).

Unstable coefficients. Remove a few data points and refit; the coefficients will change substantially. This instability is a hallmark of models that are too complex for the available data.

set.seed(123)
n_boot <- 200
boot_coefs <- matrix(NA, n_boot, p)
colnames(boot_coefs) <- names(coef(ols_full))[-1]

for (i in 1:n_boot) {
  idx <- sample(n, replace = TRUE)
  fit <- lm(EPT_index ~ ., data = wq[idx, ])
  boot_coefs[i, ] <- coef(fit)[-1]
}

# Standardize for display
boot_std <- scale(boot_coefs)

boot_long <- as_tibble(boot_std) |>
  pivot_longer(everything(), names_to = "predictor", values_to = "coef") |>
  mutate(
    true_driver = predictor %in% c("DO_mgl", "canopy_cover",
                                    "temp_C", "turbidity"),
    predictor = fct_reorder(predictor, coef, .fun = function(x) median(abs(x)),
                            .desc = TRUE)
  )

ggplot(boot_long, aes(predictor, coef, fill = true_driver)) +
  geom_boxplot(alpha = 0.7, outlier.size = 0.5, linewidth = 0.4) +
  geom_hline(yintercept = 0, linewidth = 0.4, colour = "#7a7a72") +
  scale_fill_manual(values = c("TRUE" = "#4a7c5e", "FALSE" = "#c8832a"),
                    labels = c("Noise", "True driver"),
                    name = NULL) +
  coord_flip() +
  labs(x = NULL, y = "Standardized coefficient (bootstrap distribution)",
       title = "OLS Coefficient Instability: 200 Bootstrap Samples") +
  theme_minimal(base_size = 11) +
  theme(
    plot.title      = element_text(family = "serif", face = "bold",
                                    colour = "#1e3a2f"),
    legend.position = "bottom"
  )

Coefficient instability under bootstrap resampling. Each boxplot shows the distribution of an OLS coefficient across 200 bootstrap samples. Truly important predictors (DO, canopy, temperature, turbidity) should have distributions concentrated away from zero. Noise predictors should span zero.

Preparing the data

# Predictor matrix (no intercept, no response)
x <- as.matrix(wq |> select(-EPT_index))

# Response vector
y <- wq$EPT_index

cat("Dimensions of x:", nrow(x), "x", ncol(x), "\n")

#> Dimensions of x: 80 x 14

cat("Length of y:", length(y), "\n")

#> Length of y: 80

The coefficient path

When you call glmnet(), it fits the model across a sequence of 100 lambda values (by default), from the largest lambda that sets all coefficients to zero down to a small lambda that is close to OLS. The result is a coefficient path showing how each predictor enters the model as the penalty relaxes.

lasso_fit <- glmnet(x, y, alpha = 1)  # alpha = 1 is LASSO

plot(lasso_fit, xvar = "lambda", label = TRUE,
     col = rep(c("#1e3a2f", "#4a7c5e", "#8aad96", "#c8832a",
                 "#7a3a2f", "#2d5a3d", "#e8a84a"), 2),
     lwd = 1.5)
title("LASSO Coefficient Path", line = 2.5,
      font.main = 2, col.main = "#1e3a2f")
mtext("Numbers at left identify each predictor", side = 1, line = 4,
      cex = 0.8, col = "#7a7a72")

LASSO coefficient path. Each colored line is one predictor. As lambda decreases (left to right), predictors enter the model one at a time. The four true drivers (DO, canopy cover, temperature, turbidity) enter first and have the largest coefficients. Noise predictors remain near zero or enter late.

lambda.min vs. lambda.1se

cat("lambda.min:", round(cv_fit$lambda.min, 4), "\n")

#> lambda.min: 0.1162

cat("lambda.1se:", round(cv_fit$lambda.1se, 4), "\n")

#> lambda.1se: 0.5147

cv.glmnet() reports two special lambda values:

lambda.min is the lambda that minimizes the cross-validated MSE. This gives the best prediction accuracy, but may retain some noise predictors.

lambda.1se is the largest lambda whose MSE is within one standard error of the minimum. This is a more parsimonious model – it trades a tiny bit of prediction accuracy for a simpler, more interpretable model. For environmental science, where the goal is usually understanding (which variables matter?) rather than pure prediction, lambda.1se is the better default.

Extracting coefficients

# Coefficients at lambda.1se
coefs_1se <- coef(cv_fit, s = "lambda.1se")
coefs_min <- coef(cv_fit, s = "lambda.min")

# Convert to tidy data frame
tidy_coefs <- function(coef_obj, label) {
  cf <- as.matrix(coef_obj)
  tibble(
    predictor   = rownames(cf),
    coefficient = cf[, 1],
    model       = label
  ) |>
    filter(predictor != "(Intercept)")
}

coef_df <- bind_rows(
  tidy_coefs(coefs_1se, "lambda.1se"),
  tidy_coefs(coefs_min, "lambda.min")
) |>
  mutate(selected = coefficient != 0)

# lambda.1se results
coef_1se_display <- coef_df |>
  filter(model == "lambda.1se") |>
  arrange(desc(abs(coefficient))) |>
  mutate(
    status = ifelse(selected, "Selected", "Dropped"),
    true_driver = predictor %in% c("DO_mgl", "canopy_cover",
                                    "temp_C", "turbidity")
  )

kable(coef_1se_display |>
        select(predictor, coefficient, status, true_driver),
      digits = 4,
      col.names = c("Predictor", "LASSO Coefficient", "Status",
                     "True Driver?"),
      caption = "LASSO coefficients at lambda.1se. Predictors with coefficient = 0 have been dropped from the model. The 'True Driver' column shows whether LASSO correctly identified the variables that actually drive EPT in this simulated dataset.") %>%
  kable_styling(full_width = FALSE, font_size = 14) %>%
  column_spec(3, bold = TRUE)

coef_1se_display |>
  filter(selected) |>
  ggplot(aes(reorder(predictor, abs(coefficient)), coefficient,
             fill = true_driver)) +
  geom_col(alpha = 0.8, width = 0.6, colour = "#1e3a2f") +
  geom_hline(yintercept = 0, linewidth = 0.4) +
  scale_fill_manual(values = c("TRUE" = "#4a7c5e", "FALSE" = "#c8832a"),
                    labels = c("Noise", "True driver"), name = NULL) +
  coord_flip() +
  labs(x = NULL, y = "LASSO coefficient (lambda.1se)",
       title = "Selected Predictors of EPT Index") +
  theme_minimal(base_size = 12) +
  theme(
    plot.title      = element_text(family = "serif", face = "bold",
                                    colour = "#1e3a2f"),
    legend.position = "bottom"
  )

LASSO coefficients at lambda.1se. Only selected predictors have nonzero coefficients. Green bars are true drivers; orange bars are noise predictors. A successful LASSO selects all four true drivers and excludes most noise.

Comparing lambda.min and lambda.1se

compare_df <- coef_df |>
  select(predictor, coefficient, model) |>
  pivot_wider(names_from = model, values_from = coefficient) |>
  mutate(across(c(lambda.1se, lambda.min), ~ round(.x, 4)))

kable(compare_df,
      col.names = c("Predictor", "lambda.1se", "lambda.min"),
      caption = "Coefficients at both lambda values. lambda.min typically retains more predictors than lambda.1se.") %>%
  kable_styling(full_width = FALSE, font_size = 14)

The family

When to use each

Ridge (\(\alpha = 0\)) is preferred when you believe most predictors are relevant and you just want to stabilize coefficients. It handles multicollinearity better than LASSO because it does not force a choice among correlated predictors – it splits the effect equally.

LASSO (\(\alpha = 1\)) is preferred when you believe the true model is sparse – most predictors are noise and only a few matter. This is the most common scenario in environmental variable selection.

Elastic net (\(0 < \alpha < 1\)) is preferred when correlated predictors exist in groups and you want to select or exclude the whole group together. A common choice is \(\alpha = 0.5\). If LASSO’s selection feels arbitrary among correlated predictors (keeping nitrate but dropping phosphate, or vice versa, depending on the random seed), elastic net will tend to keep both or drop both.

# Fit ridge and elastic net for comparison
ridge_fit <- glmnet(x, y, alpha = 0)
enet_fit  <- glmnet(x, y, alpha = 0.5)

par(mfrow = c(1, 3), mar = c(4.5, 4, 3, 1))

plot(lasso_fit, xvar = "lambda", main = "LASSO (alpha = 1)",
     col = "#1e3a2f", lwd = 1.2)
plot(enet_fit, xvar = "lambda", main = "Elastic Net (alpha = 0.5)",
     col = "#c8832a", lwd = 1.2)
plot(ridge_fit, xvar = "lambda", main = "Ridge (alpha = 0)",
     col = "#4a7c5e", lwd = 1.2)

Coefficient paths for LASSO (left), elastic net (center), and ridge (right). LASSO produces exact zeros; ridge never does; elastic net is intermediate. The number at the top of each panel shows how many nonzero coefficients remain at each lambda.

par(mfrow = c(1, 1))

Prediction accuracy

# Predicted values at lambda.1se
y_pred <- predict(cv_fit, newx = x, s = "lambda.1se")[, 1]

# In-sample metrics
residuals <- y - y_pred
rmse <- sqrt(mean(residuals^2))
r2   <- 1 - sum(residuals^2) / sum((y - mean(y))^2)

cat("LASSO (lambda.1se):\n")

#> LASSO (lambda.1se):

cat("  RMSE:", round(rmse, 2), "\n")

#>   RMSE: 2.33

cat("  R²:  ", round(r2, 3), "\n")

#>   R²:   0.631

# Compare to OLS
y_pred_ols <- predict(ols_full)
rmse_ols   <- sqrt(mean((y - y_pred_ols)^2))
r2_ols     <- summary(ols_full)$r.squared

cat("\nOLS (all 14 predictors):\n")

#> 
#> OLS (all 14 predictors):

cat("  RMSE:", round(rmse_ols, 2), "\n")

#>   RMSE: 1.98

cat("  R²:  ", round(r2_ols, 3), "\n")

#>   R²:   0.735

Residual plot

resid_df <- tibble(fitted = y_pred, residual = residuals)

ggplot(resid_df, aes(fitted, residual)) +
  geom_hline(yintercept = 0, linewidth = 0.4, colour = "#7a7a72",
             linetype = "dashed") +
  geom_point(colour = "#4a7c5e", alpha = 0.6, size = 2) +
  geom_smooth(method = "loess", se = FALSE, colour = "#c8832a",
              linewidth = 0.8, span = 0.7) +
  labs(x = "Fitted values (EPT index)", y = "Residuals",
       title = "LASSO Residuals vs. Fitted Values") +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(family = "serif", face = "bold",
                               colour = "#1e3a2f")
  )

Residuals vs. fitted values for the LASSO model at lambda.1se. No systematic pattern should be visible; random scatter around zero indicates the model is adequate. Curvature would suggest a missing nonlinear term.

Selection stability via repeated CV

A single CV run depends on the random fold assignment. To assess stability, repeat the entire CV procedure many times and record which predictors are selected.

set.seed(2024)
n_repeats <- 100
selection_count <- rep(0, p)
names(selection_count) <- colnames(x)

for (i in 1:n_repeats) {
  cv_i <- cv.glmnet(x, y, alpha = 1, nfolds = 10)
  coefs_i <- coef(cv_i, s = "lambda.1se")[-1, 1]
  selection_count <- selection_count + (coefs_i != 0)
}

stability_df <- tibble(
  predictor = names(selection_count),
  frequency = selection_count / n_repeats,
  true_driver = predictor %in% c("DO_mgl", "canopy_cover",
                                  "temp_C", "turbidity")
) |>
  arrange(desc(frequency))

ggplot(stability_df, aes(reorder(predictor, frequency), frequency,
                          fill = true_driver)) +
  geom_col(alpha = 0.8, width = 0.6, colour = "#1e3a2f") +
  geom_hline(yintercept = 0.5, linetype = "dashed", colour = "#7a7a72") +
  scale_fill_manual(values = c("TRUE" = "#4a7c5e", "FALSE" = "#c8832a"),
                    labels = c("Noise", "True driver"), name = NULL) +
  scale_y_continuous(labels = percent_format()) +
  coord_flip() +
  labs(x = NULL, y = "Selection frequency (100 CV runs)",
       title = "LASSO Selection Stability") +
  theme_minimal(base_size = 11) +
  theme(
    plot.title      = element_text(family = "serif", face = "bold",
                                    colour = "#1e3a2f"),
    legend.position = "bottom"
  )

Selection frequency across 100 repeated cross-validation runs. A predictor selected in 90+ runs is a robust finding; one selected in 30-60 runs is uncertain.

What to include in a methods section

When you use LASSO in a paper or report, include:

1. The number of candidate predictors and the sample size (n/p ratio).
2. The software and version (glmnet R package, version X.X).
3. The alpha value (1 for LASSO, or the elastic net value if used).
4. How lambda was chosen (10-fold CV, lambda.1se or lambda.min).
5. The set of selected predictors with their coefficients.
6. A stability assessment (repeated CV or bootstrap selection frequency).
7. The CV-estimated prediction error (MSE or RMSE).

What to include in a results section

Common pitfalls

What’s next?

• Multiple Regression & GLMs – the OLS and AIC foundation that LASSO extends
• Mixed-Effects Models – when observations are grouped, combine LASSO variable selection with random effects
• Simple Linear Regression – review the OLS fundamentals that LASSO modifies
• Correlation Analysis – detecting the multicollinearity that motivates regularization