Regularized Regression: LASSO, Ridge, and Elastic Net

What regularized regression is and why it matters in RWE

Ordinary least squares (OLS) and maximum-likelihood estimation find the coefficients that minimize the residual objective with no constraint on coefficient magnitude. In high-dimensional real-world evidence settings — where a claims analyst confronts thousands of ICD-10 diagnosis codes, CPT procedure codes, and NDC drug classes as candidate predictors — unconstrained estimators break down in two ways. First, when the number of predictors p approaches or exceeds the number of observations n, OLS has no unique solution and fits noise perfectly, producing wildly inflated coefficient variance. Second, even when n > p, OLS coefficients have high variance under multicollinearity — and ICD code families (diabetes with renal complications, diabetes without) are almost by definition collinear. Regularized regression adds a penalty term to the objective that trades a controlled amount of bias for a substantial reduction in variance, yielding estimates that generalize better to held-out data and remain estimable when p >> n.

The three penalties: geometry and intuition

All three methods minimize the same penalized objective: minimize { sum of squared (or deviance) residuals + lambda * Penalty(beta) }

The penalty type determines the geometry of the solution:

• Ridge (L2 penalty): Penalty = sum of beta_j squared. The L2 feasible region is a circle in two dimensions (sphere in higher). The optimization finds where the OLS loss ellipsoid first touches this sphere, and because the sphere has no corners, the solution almost never lands exactly at zero. Every coefficient is shrunk proportionally toward zero but retained. Ridge excels when all predictors carry some signal and when correlated code groups (all diabetes-related diagnoses, all hypertension-related procedure codes) should be kept together rather than having one arbitrarily selected and the rest discarded. Ridge also uniquely defines the coefficient vector even when the OLS design matrix is singular — critical for p > n claims analyses.

• LASSO (L1 penalty, Least Absolute Shrinkage and Selection Operator): Penalty = sum of |beta_j|. The L1 feasible region is a diamond in two dimensions (cross-polytope in higher dimensions) with sharp corners at the coordinate axes. The OLS loss ellipsoid most often first contacts the L1 region at one of these corners, where all but one or a few coordinates are zero. LASSO simultaneously estimates and selects: it produces a sparse model with many exact zeros. In a claims cost model with 4,000 candidate codes, LASSO at a cross-validated lambda typically retains 15 to 30 codes with non-zero coefficients — a compact, communicable predictive fingerprint. The sparsity comes at a cost: with strongly correlated predictors (two nearly identical complication codes), LASSO arbitrarily selects one and zeros the other, making the selected set unstable across bootstrap replicates.

• Elastic net (L1 plus L2 mixture): Penalty = alpha * sum(|beta_j|) + (1-alpha) * sum(beta_j squared). The elastic net mixes LASSO sparsity with ridge shrinkage of the selected group. When correlated code families compete for inclusion — four closely related diabetes complication codes where the science implies all four should be present — elastic net tends to keep the group together while zeroing truly irrelevant codes. The mixing parameter alpha controls the L1 fraction: alpha = 1 is pure LASSO, alpha = 0 is pure ridge. A value of alpha between 0.5 and 0.8 is a common default for claims code analyses with hierarchically organized ICD families. A further extension, the group lasso, enforces sparsity at the group level — all members of a pre-specified ICD-10 chapter or ATC drug class either enter or leave together — which is clinically natural for hierarchically organized codes.

Standardization is mandatory before penalizing

The penalty term treats all coefficients on the same footing. A coefficient of 1.0 on age (measured in decades, range 3 to 8) receives the same penalty as a coefficient of 1.0 on a binary ever/never flag (range 0 to 1). Without standardization, predictors measured on large scales are penalized less and predictors on small scales are penalized more — an artifact of units, not biology or confounding structure. Both glmnet (R) and sklearn (Python) standardize predictors internally by default and return coefficients on the original scale; always verify this behavior is active and has not been silently disabled in the implementation.

Lambda selection by cross-validation

Lambda (the penalty strength) is the critical tuning parameter. Too small: the penalty is negligible and the estimator reverts to OLS. Too large: every coefficient is shrunk to zero. In practice, lambda is chosen by k-fold cross-validation across a fine grid of values: fit on k-1 folds, predict the held-out fold, compute the loss (mean squared error for linear outcomes, binomial deviance for logistic). Two conventional choices from the resulting cross-validation curve:

• lambda.min: the lambda minimizing mean CV error — lowest bias, maximum predictive accuracy.
• lambda.1se: the largest lambda within one standard error of the minimum — sparser model, slightly higher CV error, more parsimonious. Preferred in RWE contexts where parsimony and communication matter more than the last fraction of AUC or R-squared.

The choice between lambda.min and lambda.1se should be pre-specified in the SAP or protocol and justified scientifically, not chosen post-hoc to maximize a preferred model size.

RWE high-dimensional reality: claims code spaces and p >> n settings

The primary motivation for regularized regression in RWE is the claims code landscape. A standard 365-day lookback window in a Medicare FFS or commercial database generates roughly 3,000 to 6,000 unique ICD-10 codes, 1,500 to 3,000 CPT codes, and 800 to 1,500 NDC drug classes per cohort. OLS is not identified at these dimensions; even moderate-sized cohorts face p >> n for subgroup analyses. The coordinate descent algorithm underlying glmnet cycles through predictors without requiring matrix inversion, handles both n > p and p >> n equally, and scales to millions of observations with tens of thousands of predictors on standard hardware.

In the high-dimensional propensity score (hdPS) pipeline, the empirical Bayes selection step plays a role conceptually similar to LASSO: data-adaptively identifying which of thousands of codes most predict treatment assignment. Full LASSO logistic regression on all candidate codes is increasingly used as a direct alternative or complement to hdPS for propensity score construction, with lambda chosen by cross-validation.

For biomarker panels (genomics, proteomics, pharmacogenomics) and linked registry-claims datasets, p >> n is routine. Ridge excels when all biomarkers are expected to contribute (polygenic scores, pathway-level analyses); LASSO when a small number of truly predictive biomarkers is plausible; elastic net when correlated gene families or pathway blocks should be grouped together.

Interpreting the output

A LASSO claims-cost model fit at a cross-validation-chosen lambda retains 18 of 4,000 candidate codes. The coefficient on "prior insulin use" is 0.35 on the log-cost scale.

Formal interpretation

Penalized coefficients are deliberately biased toward zero — this is the mechanism, not a failure mode. The value 0.35 is a shrunken predictive weight, not an unbiased adjusted effect estimate. It is not entitled to a naive confidence interval: the standard errors from a refitted unpenalized OLS on only the selected 18 codes are invalid. They ignore the selection step (which consumed degrees of freedom not reflected in the refitted model) and consistently understate uncertainty — naive SEs after LASSO are known to be anticonservative. Reporting them as if they were unpenalized OLS standard errors is incorrect. Post-selection inference is an active research area (selective inference, PoSI) requiring specialized software; it is not solved by simply refitting OLS on the LASSO-selected set and using the resulting standard errors.

Selection is also unstable: bootstrap the analysis cohort 200 times and the set of 18 codes changes with every replicate. Prior insulin use may appear in 70 to 80 percent of replicates as a stable predictor; the marginal code at position 18 may appear in only 20 to 30 percent. The specific coefficient 0.35 shifts in each bootstrap replicate, reflecting both sampling variation and selection variation.

Practical interpretation

Statistically correct plain English: the model found prior insulin use one of the strongest cost predictors among the 18 codes retained at the selected lambda. Treat the 18 selected codes as a predictive fingerprint — a compact feature set for scoring new patients' cost risk — not as a list of cost drivers in any causal sense. The selection of 18 codes does not imply these are the 18 most important biological or clinical determinants of cost. If the scientific question is whether insulin use causes higher costs compared with an alternative therapy, this LASSO fit does not answer it and should not be reported as if it does.

The inference warning: regularization optimizes prediction, not causal inference

Regularized regression is designed to minimize prediction error. Using LASSO to select confounders and then refitting unpenalized regression on only the selected set is a two-step procedure with invalid inference. Naive SEs from the refitted model are anticonservative and the point estimate for the treatment effect is inconsistent — the selection step induces a form of bias at the variable level that the refitted standard errors cannot recover.

The principled fix is post-double-selection (Belloni, Chernozhukov, and Hansen): run LASSO of the outcome on all candidate confounders; run LASSO of the exposure on all candidate confounders; take the union of the two selected sets; refit unpenalized regression on this union plus the exposure. This produces root-n-consistent, asymptotically normal estimates of the treatment effect even when thousands of candidate confounders are screened. An equivalent and more general approach is double/debiased ML (see parent entry predictive-and-causal-ml-models-rwe), which incorporates this logic in a cross-fit framework with formal efficiency guarantees.

A second critical rule: never penalize the exposure coefficient itself when the goal is causal effect estimation. The exposure must enter the model without penalty; only the potential confounders receive shrinkage. Penalizing the exposure biases the causal estimate toward the null — a form of dilution bias invisible to standard fit metrics.

Pros, cons, and trade-offs

Pros: Estimable when p approaches or exceeds n — the only linear framework that works in truly high-dimensional claims data without separate dimension reduction. Automatic variable selection (LASSO, elastic net) reduces the model to a compact set of predictors — a communication advantage in HEOR submissions. Ridge eliminates instability under multicollinearity; correlated code families are retained collectively rather than having arbitrary members selected. The coordinate descent algorithm (glmnet) is fast enough for millions of observations and tens of thousands of predictors. Extends to other outcomes via penalized GLMs: logistic (readmission, treatment failure), Poisson (utilization counts), Cox proportional hazards (time-to-event), multinomial. Lambda chosen by cross-validation is principled and reproducible.

Cons: Coefficients are biased by construction; they cannot be reported as unbiased adjusted effects without post-double-selection or a causal ML wrapper. Naive confidence intervals after LASSO selection are invalid and anticonservative. LASSO selection is unstable with strongly correlated predictors; elastic net or group lasso is needed when correlated groups must be represented collectively. The lambda.min vs lambda.1se choice materially affects the selected set and should be pre-specified. SAS does not provide a full cross-validated regularization path equivalent to glmnet; PROC GLMSELECT LASSO is available for moderate dimensions but requires additional calibration steps.

When NOT to use

• Small p where domain knowledge should select covariates: if a pre-specified set of 8 to 12 clinically motivated confounders is available, fit them in an unpenalized model; adding LASSO introduces bias without the high-dimensional benefit.
• Primary causal effect estimation with naive SEs: never report a LASSO-selected and refitted coefficient as an unbiased causal estimate without post-double-selection or an equivalent causal ML approach.
• Instruments or confounders chosen purely by predictive criteria: LASSO selects on outcome prediction; a code highly predictive of the outcome but connected to the exposure only through an unmeasured common cause is a collider — selecting it can induce M-bias or amplify unmeasured confounding. Data-adaptive selection cannot substitute for a causal DAG (see bias discussion in predictive-and-causal-ml-models-rwe).
• When the complete covariate set must be retained for regulatory or scientific reasons: a pre-specified confirmatory analysis with a small number of adjudicated covariates should use unpenalized regression.
• When LASSO selection is the deliverable but collinearity is severe: use elastic net or group lasso instead; pure LASSO with correlated predictors produces an arbitrary sparse representation of the correlation structure, not a scientifically meaningful feature set.

Decision diagram

flowchart TD
  Q["High-dimensional covariate space?<br/>(p > 50 or p > n)"] -->|yes| Pen["Use penalized regression"]
  Q -->|no| OLS["Use OLS / unpenalized logistic<br/>with pre-specified covariates"]
  Pen --> Goal["Primary goal?"]
  Goal -->|prediction / scoring| Pred["Tune lambda by CV<br/>cv.glmnet / LassoCV"]
  Goal -->|causal effect| Causal["Apply post-double-selection<br/>or double-ML wrapper<br/>NEVER penalize exposure"]
  Pred --> Corr["Correlated code families?"]
  Corr -->|yes, keep group| Enet["Elastic net<br/>alpha 0.5-0.8"]
  Corr -->|no, want sparsity| Lasso["LASSO<br/>alpha = 1"]
  Corr -->|keep all predictors| Ridge["Ridge<br/>alpha = 0"]
  Lasso --> Lam["Lambda choice: pre-specify<br/>lambda.min (max accuracy)<br/>lambda.1se (parsimonious)"]
  Enet --> Lam
  Ridge --> Lam
  Lam --> Stab["Bootstrap stability check<br/>200 replicates, track selection frequency"]
  Stab --> Interp["Report as predictive weights<br/>NOT unbiased causal effects<br/>Naive CIs after LASSO are INVALID"]

Worked example

Scenario

A HEOR analyst builds a log-cost prediction model for 500 commercial insurance members who initiated a diabetes drug. After constructing a 365-day lookback feature matrix with hundreds of candidate codes, the analyst uses three illustrative predictors to demonstrate ridge and LASSO penalization. The OLS fit gives the coefficients shown in the table below (all on the standardized scale: mean 0, unit variance per predictor). The analyst compares ridge (lambda = 4) and LASSO (lambda = 1) to see which predictors survive penalization and by how much.

Dataset

OLS log-cost model coefficients on the standardized predictor scale. Three predictors shown for illustration; a real analysis would have hundreds of code-level features.

Steps

Result

Ridge (lambda = 4): prior_insulin_use = 1.0, prior_ed_visit_count = 0.6, comorbidity_flag = 0.1. All 3 predictors retained, each shrunk proportionally. LASSO (lambda = 1): prior_insulin_use = 4.0, prior_ed_visit_count = 2.0, comorbidity_flag = 0.0 (zeroed out and excluded from the model). 2 of 3 predictors selected. Ridge shrinks coefficients proportionally toward zero; LASSO performs automatic variable selection by zeroing the weakest coefficient below the penalty threshold.

Trade-offs

Runnable example

import numpy as np
import pandas as pd
from sklearn.linear_model import LassoCV, RidgeCV, ElasticNetCV, LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.model_selection import GridSearchCV

# ── Input: df with columns [person_id, outcome, ...baseline covariates] ──
# outcome is continuous (log-cost) for LASSO/ridge, binary (0/1) for logistic.
covar_cols = [c for c in df.columns if c not in ("person_id", "outcome")]
X = df[covar_cols].values
y = df["outcome"].values

# ── 1. LASSO (L1) with cross-validated lambda — continuous outcome ──
# Pipeline ensures StandardScaler is fit only on training folds (no leakage).
lasso_pipe = Pipeline([
    ("scaler", StandardScaler()),          # mandatory: standardize before penalizing
    ("lasso", LassoCV(cv=10, max_iter=10_000, random_state=42))
])
lasso_pipe.fit(X, y)
lasso_m = lasso_pipe.named_steps["lasso"]
print(f"LASSO chosen alpha (lambda): {lasso_m.alpha_:.5f}")
coefs = pd.Series(lasso_m.coef_, index=covar_cols)
selected = coefs[coefs != 0].sort_values(key=abs, ascending=False)
print(f"Selected predictors: {len(selected)} of {len(coefs)}")
print(selected.head(20))   # top 20 by absolute magnitude

# ── 2. Ridge (L2) — all predictors retained, fights multicollinearity ──
ridge_pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("ridge", RidgeCV(alphas=np.logspace(-3, 4, 100), cv=10))
])
ridge_pipe.fit(X, y)
ridge_m = ridge_pipe.named_steps["ridge"]
print(f"\nRidge chosen alpha (lambda): {ridge_m.alpha_:.5f}")
# All coefficients are non-zero; ridge retains every predictor.

# ── 3. Elastic net — L1+L2 mixture for correlated code families ──
enet_pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("enet", ElasticNetCV(
        l1_ratio=[0.1, 0.5, 0.7, 0.9, 0.95, 1.0],
        cv=10, max_iter=10_000, random_state=42
    ))
])
enet_pipe.fit(X, y)
enet_m = enet_pipe.named_steps["enet"]
print(f"\nElastic net: l1_ratio = {enet_m.l1_ratio_:.2f}, alpha = {enet_m.alpha_:.5f}")

# ── 4. LASSO logistic for binary outcomes (readmission, event within horizon) ──
# In sklearn, C = 1/lambda; smaller C = more regularization.
# CRITICAL: exposure variable must NOT be penalized — use a separate unpenalized step
# or manually set penalty_factor equivalent (not natively in sklearn pipeline;
# consider using glmnet via rpy2 for the penalty.factor argument).
Cs = np.logspace(-4, 2, 50)
log_pipe = Pipeline([
    ("scaler", StandardScaler()),
    ("logit", LogisticRegression(penalty="l1", solver="saga", max_iter=5_000))
])
gs = GridSearchCV(log_pipe, {"logit__C": Cs}, cv=10, scoring="neg_log_loss")
gs.fit(X, y)
best_C = gs.best_params_["logit__C"]
print(f"\nLASSO logistic: best lambda = {1/best_C:.5f} (C = {best_C:.5f})")
coefs_logit = pd.Series(
    gs.best_estimator_.named_steps["logit"].coef_[0], index=covar_cols
)
print(f"Selected: {(coefs_logit != 0).sum()} predictors")

# ── 5. Bootstrap stability (200 replicates) ──
B = 200
rng = np.random.default_rng(0)
selections = np.zeros((len(covar_cols), B), dtype=int)
for b in range(B):
    idx = rng.choice(len(X), size=len(X), replace=True)
    pipe_b = Pipeline([("sc", StandardScaler()),
                       ("ls", LassoCV(cv=5, max_iter=5_000, random_state=b))])
    pipe_b.fit(X[idx], y[idx])
    selections[:, b] = (pipe_b.named_steps["ls"].coef_ != 0).astype(int)
stability = pd.Series(selections.mean(axis=1), index=covar_cols)
print("\nTop 10 most stable predictors (fraction of replicates selected):")
print(stability.sort_values(ascending=False).head(10))

library(glmnet)

# ── Prepare data: X must be a numeric matrix ──
# For claims data convert factor predictors via model.matrix (removes intercept).
X <- model.matrix(~ . - 1, data = df[, covariate_cols])
y <- df[["log_cost"]]      # continuous outcome; use 0/1 for binomial, Surv() for cox

# ── 1. LASSO (alpha = 1): variable selection + shrinkage ──
set.seed(42)
fit_lasso <- cv.glmnet(X, y, alpha = 1, family = "gaussian",
                       nfolds = 10, standardize = TRUE)
cat("LASSO lambda.min:", round(fit_lasso$lambda.min, 5), "\n")
cat("LASSO lambda.1se:", round(fit_lasso$lambda.1se, 5), "\n")

# Coefficients at lambda.1se (parsimonious; pre-specify in SAP)
coef_1se  <- coef(fit_lasso, s = "lambda.1se")
coef_min  <- coef(fit_lasso, s = "lambda.min")
n_1se <- sum(coef_1se[-1, 1] != 0)   # exclude intercept
n_min <- sum(coef_min[-1, 1] != 0)
cat("Predictors at lambda.1se:", n_1se, "| at lambda.min:", n_min, "\n")

# Top selected predictors at lambda.1se, sorted by |coefficient|
sel_df <- data.frame(
  predictor = rownames(coef_1se)[-1],
  coef      = coef_1se[-1, 1]
)
sel_df <- sel_df[sel_df$coef != 0, ]
sel_df <- sel_df[order(abs(sel_df$coef), decreasing = TRUE), ]
print(head(sel_df, 20))

# ── 2. Ridge (alpha = 0): shrinkage only, all predictors retained ──
fit_ridge <- cv.glmnet(X, y, alpha = 0, family = "gaussian", nfolds = 10)
cat("\nRidge lambda.min:", round(fit_ridge$lambda.min, 5), "\n")

# ── 3. Elastic net: search over alpha values ──
alphas    <- c(0.2, 0.5, 0.8, 0.9, 1.0)
cv_errors <- sapply(alphas, function(a) {
  fit <- cv.glmnet(X, y, alpha = a, family = "gaussian", nfolds = 10)
  min(fit$cvm)   # minimum CV mean squared error
})
best_alpha <- alphas[which.min(cv_errors)]
fit_enet   <- cv.glmnet(X, y, alpha = best_alpha, family = "gaussian", nfolds = 10)
cat("Elastic net best alpha:", best_alpha, "\n")

# ── 4. Penalized logistic for binary outcomes ──
# CRITICAL: protect the exposure from penalization via penalty.factor.
# Assume column 1 of X is the exposure (treat); all others are confounders.
pf <- rep(1, ncol(X))
pf[1] <- 0   # exposure enters unpenalized
fit_bin <- cv.glmnet(X, df$event, alpha = 1, family = "binomial",
                     nfolds = 10, penalty.factor = pf)

# ── 5. Penalized Cox for time-to-event outcomes ──
library(survival)
fit_cox <- cv.glmnet(X, Surv(df$time, df$event_indicator),
                     alpha = 1, family = "cox", nfolds = 10)
cat("\nCox LASSO selected at lambda.1se:",
    sum(coef(fit_cox, s = "lambda.1se")[, 1] != 0), "predictors\n")

# ── 6. Bootstrap selection stability ──
B     <- 200
n     <- nrow(X)
boot  <- matrix(0, nrow = ncol(X), ncol = B)
for (b in seq_len(B)) {
  idx    <- sample(n, replace = TRUE)
  fit_b  <- cv.glmnet(X[idx, ], y[idx], alpha = 1, nfolds = 5)
  boot[, b] <- as.integer(coef(fit_b, s = "lambda.1se")[-1, 1] != 0)
}
stability <- rowMeans(boot)
names(stability) <- colnames(X)
cat("\nTop stable predictors (fraction selected across", B, "bootstrap replicates):\n")
print(sort(stability, decreasing = TRUE)[1:10])

/* ── Standardize predictors before penalization ── */
proc standard data=work.analytic out=work.std_data mean=0 std=1;
  var &covariates;   /* space-delimited list of baseline predictor names */
run;

/* ── 1. LASSO for continuous outcome (log-cost) ── */
/* SELECTION=LASSO: L1 penalty path.                                                   */
/* STOP=CV CHOOSE=CV: stop and choose by 10-fold cross-validation error.               */
proc glmselect data=work.std_data plots=all;
  model log_cost = &covariates / selection=lasso
                                 stop=cv(method=random(seed=42) groups=10)
                                 choose=cv;
  output out=work.lasso_scored predicted=pred_lasso;
run;
/* Inspect the ANOVA table and parameter estimates for the selected model.              */
/* Predictors with non-zero estimates at the chosen step are selected by LASSO.         */

/* ── 2. Elastic net (L1 + L2) ── */
/* L2= specifies the ridge penalty weight added to the LASSO path.                     */
proc glmselect data=work.std_data;
  model log_cost = &covariates / selection=elasticnet(l2=0.1)
                                 stop=cv(method=random(seed=42) groups=10)
                                 choose=cv;
run;

/* ── 3. LASSO logistic for binary outcomes ── */
/* Available in SAS 9.4 TS1M5 and later. SLENTRY controls the entry threshold.        */
proc logistic data=work.std_data;
  model readmission(event='1') = &covariates / selection=lasso slentry=0.05;
run;
/* Note: SAS LASSO logistic does not use a CV criterion equivalent to cv.glmnet;       */
/* for a full CV-tuned path on large claims data, prefer R glmnet or Python sklearn.   */

/* ── 4. Score held-out data with selected model ── */
/* After PROC GLMSELECT, use PROC PLM or the OUTPUT statement from the model step      */
/* to score new patients. The selected (non-zero) predictors and their coefficients    */
/* are embedded in the parameter estimates output dataset for deployment.              */
proc plm restore=work.lasso_model;
  score data=work.holdout out=work.holdout_scored predicted=pred_lasso;
run;

Citations

1. [1]Tibshirani R. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society Series B (Methodological). 1996;58(1):267-288.
2. [2]Zou H, Hastie T. Regularization and variable selection via the elastic net. Journal of the Royal Statistical Society Series B (Statistical Methodology). 2005;67(2):301-320.

1. [3]Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software. 2010;33(1):1-22.

1. [4]Simon N, Friedman J, Hastie T, Tibshirani R. Regularization paths for Cox's proportional hazards model via coordinate descent. Journal of Statistical Software. 2011;39(5):1-13.