Regularized Regression: LASSO, Ridge, and Elastic Net

In plain language

Regularized regression (LASSO, ridge, and elastic net) is a family of statistical models that deliberately shrink coefficient estimates toward zero — and in the LASSO case, zero them out entirely — to prevent the model from overfitting when there are many more candidate predictors than can reliably be estimated from the data. In health research with claims data, this means a model can start with thousands of diagnosis and procedure codes and automatically find the handful that best predict costs, hospitalizations, or treatment choice, without the analyst having to hand-pick covariates. The key trade-off is that the shrunken coefficients are intentionally biased and cannot be directly interpreted as causal effects, only as predictive weights.

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

• LASSO (L1 penalty, Least Absolute Shrinkage and Selection Operator): Penalty = sum of

• Elastic net (L1 plus L2 mixture): Penalty = alpha sum(|beta_j|) + (1-alpha)

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,

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 to use

• High-dimensional prediction with tens to thousands of candidate predictors from a claims
• Propensity score construction from thousands of code candidates as an alternative to
• Confounder screening: reduce 4,000 codes to a manageable set of 20 to 50 using LASSO
• Elastic net for correlated code families (ICD-10 chapters, ATC drug classes) where the
• Nuisance model fitting inside double-ML or TMLE as one candidate in a super learner,
• Penalized Cox regression for time-to-event outcomes in linked registries when many candidate

When NOT to use

• Small p where domain knowledge should select covariates: if a pre-specified set of 8 to 12
• Primary causal effect estimation with naive SEs: never report a LASSO-selected and refitted
• Instruments or confounders chosen purely by predictive criteria: LASSO selects on outcome
• When the complete covariate set must be retained for regulatory or scientific reasons: a
• When LASSO selection is the deliverable but collinearity is severe: use elastic net or group

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

• Standardize all predictors before fitting. Each predictor is centered to mean zero and scaled to unit variance so the penalty treats them equally regardless of original units. Both glmnet (R) and sklearn (Python) perform this internally by default. The OLS coefficients in the table are already expressed on this standardized scale.

• Ridge with lambda = 4: the ridge formula shrinks each OLS coefficient by the factor 1 / (1 + lambda). With lambda = 4.0: prior_insulin_use is penalized to 5.0 (1.0/(1.0+4.0)) = 5.0 0.2 = 1.0; prior_ed_visit_count is penalized to 3.0 (1.0/(1.0+4.0)) = 3.0 0.2 = 0.6; comorbidity_flag is penalized to 0.5 (1.0/(1.0+4.0)) = 0.5 0.2 = 0.1. All three coefficients are shrunk toward zero, but none reach exactly zero — ridge retains every predictor in the model.

• LASSO with lambda = 1: the soft-threshold formula is beta_lasso = sign(ols_beta) * max(|ols_beta| - lambda, 0). With lambda = 1.0: prior_insulin_use has ols_beta = 5.0, which exceeds lambda, so the LASSO coefficient is 5.0 - 1.0 = 4.0; prior_ed_visit_count has ols_beta = 3.0, which exceeds lambda, so the LASSO coefficient is 3.0 - 1.0 = 2.0; comorbidity_flag has ols_beta = 0.5, which is less than lambda = 1.0, so its coefficient is zeroed out and the predictor is excluded from the model entirely.

• LASSO performs variable selection; ridge does not. With lambda = 1, LASSO keeps 2 of 3 predictors and zeros the weakest (comorbidity_flag). In a full analysis with 500 candidate codes, a cross-validated lambda would typically keep 15 to 30 codes. Choose lambda.1se for the most parsimonious model; lambda.min for best predictive accuracy. Pre-specify this choice in the SAP.

• Caution on interpretation: the retained LASSO coefficients of 4.0 and 2.0 are penalized predictive weights, not unbiased causal effects. They are deliberately biased toward zero by design. Their naive standard errors — from refitting OLS on just these two selected predictors — would be invalid (anticonservative) because they ignore the selection step. Treat these as weights in a scoring model, not as adjusted effect estimates. For causal inference, apply post-double-selection or use double-ML (see parent entry).

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.

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])