L1-Penalized Logistic Regression: Performance vs Complexity

In applied data science, strong predictive performance is only one part of a successful model. In many settings, especially in regulated industries, models are expected to support decisions in a transparent and reliable way. A model that is marginally more accurate but substantially more complex can be difficult to interpret, harder to maintain, and less stable when deployed in new environments.

Sparse models offer a practical alternative. They are easier to validate, simpler to explain, and often more robust across datasets. At the same time, overly simple models risk missing important signal. The central question is how to balance predictive performance with model complexity in a principled way.

Start here

Article Overview

L1-penalized logistic regression drives coefficients to exactly zero, so the penalty strength, set through C, decides how many survive. Dropping coefficients can cost performance, and this article measures that exchange on heart-disease data: 918 patients, 11 clinical features, and a tolerance rule that trades a fixed margin of validation area under the precision-recall curve (AUPRC) for a sparser model.

L1 selects, L2 does not. L1 shrinks coefficients to exactly zero while L2 only shrinks them toward it, so the penalty strength, set through C, sets complexity, measured here as the count of non-zero coefficients.
Five figures. Three interactive panels compare the L1 and L2 coefficient paths, apply a tolerance rule to the performance-complexity curve, and track how often each feature survives resampling. Two more plot the C sweep and the held-out precision-recall curve (PRC) and receiver operating characteristic (ROC) curves.
The code path. Eleven Python blocks build the pipeline, sweep C over 200 log-spaced values under nested cross-validation, and score the best-AUPRC model against the sparsest models within 1%, 5%, and 10% of its validation score.

Read every result below with this in mind. A zero coefficient indicates that a predictor was not selected in this fit, and does not definitively determine whether it carries signal. Correlated predictors trade places from one resample to the next, so the selection frequencies across resamples carry more weight than any single fitted model.

L1 regularization as a mechanism for sparsity

L1-regularized logistic regression provides a direct way to control model complexity. By penalizing the absolute values of the coefficients, the method encourages many coefficients to shrink exactly to zero. This effectively performs feature selection during training.

As the strength of regularization increases, more coefficients are set to zero and the model becomes simpler. As regularization weakens, more features enter the model, capturing additional structure in the data. This creates a natural continuum of models ranging from highly sparse to more complex.

This behavior is particularly useful in domains where interpretability matters. A model that relies on a small set of meaningful predictors is easier to audit and more aligned with real-world constraints such as missing data or measurement variability.

Figure 1 lets you vary the penalty strength and compare the L1 and L2 penalties directly, so you can watch coefficients reach exactly zero as the penalty grows.

Figure 1

The L1 penalty shrinks coefficients to exactly zero

Drag the slider to change the penalty strength, then switch the penalty type. Under L1 the coefficients drop out one by one as the penalty grows. Under L2 they shrink toward zero but never reach it, so no feature is ever removed.

Penalty parameter CC = 0.130

Penalty type

Penalty L1 (lasso) at C = 0.130: 5 of 9 coefficients are non-zero. The rest have been set to exactly zero and dropped from the model.

Worst concave points1.09
Worst perimeter0.79
Worst texture0.45
Mean concavity0.27
Worst area0.09
Worst smoothness0 (removed)
Mean texture0 (removed)
SE compactness0 (removed)
Mean fractal dimension0 (removed)

Illustrative coefficient paths for a breast-cancer diagnosis model built from cell-nucleus measurements, a separate example from the heart-disease case study later in the article. The values are drawn to show the shape of the paths. Only an exact zero drops a feature from the model, so a very small coefficient under L2 still counts as one the model uses.

Reading the figure. The count of non-zero coefficients is the measure of model complexity used throughout this article. Smaller C means a stronger penalty and a sparser model, because C is the inverse of the penalty strength.

Defining model complexity

In this setting, model complexity can be defined in a straightforward way as the number of non-zero coefficients. This provides a clear and interpretable measure of how many features are actively used. Tracking this quantity across different regularization strengths allows for a structured comparison of models. It also shifts the focus away from performance metrics alone and introduces a second axis for evaluation.

The performance–complexity curve

A practical way to study the trade-off is to construct a performance–complexity curve. This involves training a sequence of models across a grid of regularization strengths and evaluating each model on held-out data.

For each model, two quantities are recorded:
• Predictive performance, such as area under the curve (AUC) or AUPRC
• Model complexity, defined as the number of non-zero coefficients

When performance is plotted against complexity, a consistent pattern often emerges. Performance improves quickly as the most informative predictors are included. After this initial phase, gains become smaller as additional variables contribute less signal. This plateau region is where model selection becomes most meaningful.

Figure 2, in the next section, plots this curve and lets you apply the selection rule directly.

Diminishing returns and model selection

The presence of diminishing returns creates an opportunity to favor simpler models without sacrificing much performance. Instead of selecting the single best-performing model, a more robust strategy is to identify a set of near-optimal models and then choose the simplest among them.

A common approach follows three steps:
• Identify the model with the highest validation performance
• Define a tolerance range, such as within 5 percent of the maximum
• Select the model within this range that has the fewest non-zero coefficients

This method formalizes the trade-off and avoids overfitting to small and potentially unstable differences in performance.

Figure 2 turns this three-step rule into an interactive curve. Choose a tolerance to see which model the rule selects and how much validation performance that choice gives up.

Figure 2

Performance rises quickly, then plateaus

Every point is a model with a given number of non-zero coefficients. Pick a tolerance to select the sparsest model whose validation AUPRC stays within that margin of the best model, or switch to performance-first selection to keep the highest score.

Selection rule

Tolerance below best AUPRC

Selection rule

Parsimony-aware

Selected complexity

4 of 20

Selected AUPRC

0.955

Cost vs best

−0.035 vs best

Parsimony-aware (within 5%): the sparsest model whose AUPRC stays above the threshold uses 4 coefficients at AUPRC 0.955, costing 0.035 AUPRC versus the best model.

Illustrative curve for the same breast-cancer diagnosis model, a separate example from the heart-disease case study later in the article. It shows the rise-then-plateau shape this trade-off produces. The selected model is simply the smallest one that satisfies the rule. A gap smaller than the shaded interval sits inside the resampling spread.

Reading the figure. Because the curve flattens, a much simpler model can come close to the best model at a small cost in AUPRC. Higher AUPRC is better, and complexity is the number of non-zero coefficients.

A practical detail: stability of feature selection

One important consideration is that feature selection under L1 regularization can be unstable, particularly when predictors are correlated. Small changes in the data can lead to different subsets of selected features, even when overall model performance remains similar.

A useful best practice is to evaluate feature selection stability using resampling methods such as cross-validation or bootstrapping. Features that are consistently selected across folds are more likely to represent true signal rather than noise. This perspective shifts attention from a single fitted model to patterns that persist across multiple samples of the data.

Figure 3 illustrates this variability by resampling the data and tracking how often each feature keeps a non-zero coefficient.

Figure 3

Which features are selected can change from sample to sample

Each bar is how often a feature keeps a non-zero coefficient across resampled folds. Draw a new resample to watch the correlated pair trade places. The frequency bars stay fixed because they summarize all 20 resamples.

Resampling

Resample 1 of 20

This resample keeps Mean perimeter (Mean radius dropped). Consistently selected features (long bars) are the stable signal. The correlated pair trades places from one resample to the next while the frequency bars stay fixed (they summarize all 20 resamples).

Illustrative resampling results for the same breast-cancer diagnosis model, a separate example from the heart-disease case study later in the article. A high selection frequency means the model chooses a feature consistently, which is evidence of a stable signal, and evidence of causation comes from study design.

Reading the figure. L1 tends to keep one feature from a correlated pair and drop the other, and the choice varies across samples. Features selected in almost every resample are the more trustworthy part of the model.

Interpreting results in context

Model selection rarely depends on metrics alone. In practice, a slightly less accurate model may be preferable if it is easier to use and interpret.

For example, simpler models may be favored when they:
• Rely on features that are routinely available
• Align with domain knowledge
• Are easier to communicate to stakeholders
• Generalize more reliably across settings

In clinical applications, these factors often determine whether a model is adopted in practice.

A practical decision framework

The worked example that follows implements steps 1 to 6. It does not run the step 7 stability check. The steps turn the trade-off into a repeatable routine.

1
Fit a grid of models. Train L1-penalized logistic regression across a range of regularization strengths C.
2
Estimate performance out of sample. Use nested cross-validation and score each candidate with AUPRC.
3
Record complexity. Track the average number of non-zero coefficients for each candidate.
4
Find the best score. Identify the model with the highest cross-validated AUPRC.
5
Set a tolerance. Decide how much performance you are willing to give up, for example staying within 5% of the best AUPRC.
6
Take the simplest qualifier. Among models within the tolerance, choose the one with the fewest non-zero coefficients.
7
Check stability. Confirm the selected features recur across resamples before you trust the interpretation.

Example: Evaluating the trade-off in L1-penalized logistic regression

LASSO Sparsity. How L1-regularized logistic regression balances predictive performance against model complexity, using nested cross-validation to find the sparsest model that still scores near the best.

Hover any card to explore

Overview

We run LASSO (L1-regularized) logistic regression on a clinical heart-disease dataset, then trace the trade-off between predictive performance and model complexity, counting complexity as the number of non-zero coefficients.

Key questions explored

How does predictive performance (AUPRC) change as regularization is relaxed and more features are included?
Where does the performance–complexity curve plateau, the point of diminishing returns?
Can we identify a sparser model that retains most of the predictive value of the globally best model?

Methodology

Nested cross-validation: An outer loop estimates generalization performance; an inner loop selects the optimal regularization strength C.
Repeated stratified folds: 100 repetitions × 5 outer folds produce stable performance estimates and 95% resampling intervals.
Four model variants are compared: the globally best model (Model A) and the sparsest models within 1% (Model B), 5% (Model C), and 10% (Model D) of maximum validation performance.

Dataset

The Heart Failure Prediction Dataset from Kaggle contains 918 patients with 11 clinical and demographic features and a binary outcome indicating the presence of heart disease. It has no empty cells, though Cholesterol is recorded as 0 for 172 of the 918 patients, and resting blood pressure for one, physiologically impossible readings that stand in for missing measurements. The notebook uses the columns as recorded, so the cholesterol coefficient partly reflects that sentinel group. The dataset remains well-suited for illustrating regularization trade-offs.

1. Import dependencies

Python

import os
from collections import defaultdict

# Numerical computing and data manipulation
import numpy as np
import pandas as pd

# Parallel processing across outer CV folds
from joblib import Parallel, delayed

# Scikit-learn: pipeline construction and preprocessing
from sklearn.base import clone
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.experimental import enable_iterative_imputer  # noqa: F401 — required to activate IterativeImputer
from sklearn.impute import IterativeImputer, SimpleImputer

# Scikit-learn: model, cross-validation, and evaluation metrics
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV, StratifiedKFold, RepeatedStratifiedKFold
from sklearn.metrics import (
    make_scorer,
    average_precision_score,
    roc_auc_score,
    roc_curve,
    precision_recall_curve,
)

# Spline smoothing for visualising noisy C-vs-performance curves
from scipy.interpolate import UnivariateSpline

# Interactive plotting
import plotly.graph_objects as go
from plotly.subplots import make_subplots

# Dataset acquisition from Kaggle
import kagglehub
from kagglehub import KaggleDatasetAdapter

import os
from collections import defaultdict

# Numerical computing and data manipulation
import numpy as np
import pandas as pd

# Parallel processing across outer CV folds
from joblib import Parallel, delayed

# Scikit-learn: pipeline construction and preprocessing
from sklearn.base import clone
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.experimental import enable_iterative_imputer  # noqa: F401 — required to activate IterativeImputer
from sklearn.impute import IterativeImputer, SimpleImputer

# Scikit-learn: model, cross-validation, and evaluation metrics
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import GridSearchCV, StratifiedKFold, RepeatedStratifiedKFold
from sklearn.metrics import (
    make_scorer,
    average_precision_score,
    roc_auc_score,
    roc_curve,
    precision_recall_curve,
)

# Spline smoothing for visualising noisy C-vs-performance curves
from scipy.interpolate import UnivariateSpline

# Interactive plotting
import plotly.graph_objects as go
from plotly.subplots import make_subplots

# Dataset acquisition from Kaggle
import kagglehub
from kagglehub import KaggleDatasetAdapter

2. Load data and define constants

We load the Heart Failure Prediction Dataset. If not already cached locally, it is downloaded from Kaggle via kagglehub. Column names are renamed to be human-readable.

Dataset characteristics

Property	Value
Samples	918 patients
Features	11 (mix of numerical and categorical)
Target	Heart Disease (1 = present, 0 = absent)

Feature list: Age, Sex, Chest Pain Type, Resting Blood Pressure, Cholesterol, Fasting Blood Sugar, Resting electrocardiogram (ECG), Max Heart Rate, Exercise Angina, Old Peak (ST depression), ST Slope.

Key constants

Constant	Value	Purpose
`N_TRIALS`	100	Outer CV repetitions. More trials pin down the resampling interval more precisely
`N_INNER_SPLITS`	5	Inner folds for hyperparameter selection
`N_OUTER_SPLITS`	5	Outer folds for out-of-sample performance estimation
`RANDOM_STATE`	42	Global reproducibility seed
`CATEGORICAL_THRESHOLD`	10	Features with ≤10 unique values are treated as categorical

Runtime note: With N_TRIALS=100, the nested CV run will take several minutes depending on hardware. To explore the code quickly, reduce N_TRIALS to 5-10. Resampling intervals will be far less reliable but the workflow is identical.

Python

DATA_FILE = "heart.csv"

# Download from Kaggle if not already cached locally
if not os.path.exists(DATA_FILE):
    df = kagglehub.load_dataset(
        KaggleDatasetAdapter.PANDAS,
        "fedesoriano/heart-failure-prediction",
        "heart.csv",
    )
    df.to_csv(DATA_FILE, index=False)
else:
    df = pd.read_csv(DATA_FILE)

# Rename columns to human-readable form for display and interpretation
df.rename(
    columns={
        "ChestPainType": "Chest Pain Type",
        "RestingBP": "Resting Blood Pressure",
        "FastingBS": "Fasting Blood Sugar",
        "RestingECG": "Resting ECG",
        "MaxHR": "Max Heart Rate",
        "ExerciseAngina": "Exercise Angina",
        "Oldpeak": "Old Peak",
        "ST_Slope": "ST Slope",
        "HeartDisease": "Heart Disease",
    },
    inplace=True,
)

TARGET_COLUMN = "Heart Disease"
CATEGORICAL_THRESHOLD = 10  # Columns with ≤10 unique values are treated as categorical
N_TRIALS = 100              # Number of outer CV repetitions (reduce to 5–10 for quick exploration)
N_INNER_SPLITS = 5          # Folds for inner hyperparameter selection
N_OUTER_SPLITS = 5          # Folds for outer performance estimation
RANDOM_STATE = 42

np.random.seed(RANDOM_STATE)

# Automatically classify each feature as categorical or numerical.
# Binary flags (e.g. Fasting Blood Sugar: 0/1) are treated as categorical even though
# they are stored as integers, because their cardinality is low.
categorical_features = [
    col for col in df.columns
    if col != TARGET_COLUMN and (
        df[col].dtype == "object"
        or df[col].dtype == "bool"
        or (df[col].dtype in ("int64", "float64") and df[col].nunique() <= CATEGORICAL_THRESHOLD)
    )
]
numerical_features = [col for col in df.columns if col not in categorical_features + [TARGET_COLUMN]]

# Split into feature matrix and binary target vector
X, y = df.drop(columns=TARGET_COLUMN), df[TARGET_COLUMN]

print(f"Dataset shape: {X.shape[0]} samples, {X.shape[1]} features")
print(f"Target prevalence: {y.mean():.1%} positive (Heart Disease)")
print(f"Categorical features ({len(categorical_features)}): {categorical_features}")
print(f"Numerical features  ({len(numerical_features)}): {numerical_features}")

DATA_FILE = "heart.csv"

# Download from Kaggle if not already cached locally
if not os.path.exists(DATA_FILE):
    df = kagglehub.load_dataset(
        KaggleDatasetAdapter.PANDAS,
        "fedesoriano/heart-failure-prediction",
        "heart.csv",
    )
    df.to_csv(DATA_FILE, index=False)
else:
    df = pd.read_csv(DATA_FILE)

# Rename columns to human-readable form for display and interpretation
df.rename(
    columns={
        "ChestPainType": "Chest Pain Type",
        "RestingBP": "Resting Blood Pressure",
        "FastingBS": "Fasting Blood Sugar",
        "RestingECG": "Resting ECG",
        "MaxHR": "Max Heart Rate",
        "ExerciseAngina": "Exercise Angina",
        "Oldpeak": "Old Peak",
        "ST_Slope": "ST Slope",
        "HeartDisease": "Heart Disease",
    },
    inplace=True,
)

TARGET_COLUMN = "Heart Disease"
CATEGORICAL_THRESHOLD = 10  # Columns with ≤10 unique values are treated as categorical
N_TRIALS = 100              # Number of outer CV repetitions (reduce to 5–10 for quick exploration)
N_INNER_SPLITS = 5          # Folds for inner hyperparameter selection
N_OUTER_SPLITS = 5          # Folds for outer performance estimation
RANDOM_STATE = 42

np.random.seed(RANDOM_STATE)

# Automatically classify each feature as categorical or numerical.
# Binary flags (e.g. Fasting Blood Sugar: 0/1) are treated as categorical even though
# they are stored as integers, because their cardinality is low.
categorical_features = [
    col for col in df.columns
    if col != TARGET_COLUMN and (
        df[col].dtype == "object"
        or df[col].dtype == "bool"
        or (df[col].dtype in ("int64", "float64") and df[col].nunique() <= CATEGORICAL_THRESHOLD)
    )
]
numerical_features = [col for col in df.columns if col not in categorical_features + [TARGET_COLUMN]]

# Split into feature matrix and binary target vector
X, y = df.drop(columns=TARGET_COLUMN), df[TARGET_COLUMN]

print(f"Dataset shape: {X.shape[0]} samples, {X.shape[1]} features")
print(f"Target prevalence: {y.mean():.1%} positive (Heart Disease)")
print(f"Categorical features ({len(categorical_features)}): {categorical_features}")
print(f"Numerical features  ({len(numerical_features)}): {numerical_features}")

Dataset shape: 918 samples, 11 features
Target prevalence: 55.3% positive (Heart Disease)
Categorical features (6): ['Sex', 'Chest Pain Type', 'Fasting Blood Sugar', 'Resting ECG', 'Exercise Angina', 'ST Slope']
Numerical features  (5): ['Age', 'Resting Blood Pressure', 'Cholesterol', 'Max Heart Rate', 'Old Peak']

3. Preprocessing pipeline

We construct a scikit-learn Pipeline that chains preprocessing with the classifier. Encapsulating everything in a single pipeline is essential for correct nested cross-validation. All preprocessing steps (imputation, scaling, encoding) are fitted only on training data within each fold, preventing data leakage.

Numerical features

Iterative imputation (IterativeImputer): Models each feature with missing values as a function of all other features, iterating to convergence. This is more principled than mean imputation when missingness is related to other observed variables, the missing-at-random case. It does not recover missingness that depends on the unobserved value itself, and unlike the categorical branch below it adds no indicator that the value was absent.
Standardization (StandardScaler): Centers and scales to zero mean and unit variance. Required for L1 regularization to penalize all coefficients on the same scale.

Categorical features

Constant imputation: Missing values are filled with a dedicated "missing" category rather than the most frequent value, preserving the information that the value was absent.
One-hot encoding (OneHotEncoder, drop="first"): Encodes categories as binary indicator variables. Dropping the first level avoids perfect multicollinearity (the dummy variable trap).

Classifier

L1-regularized logistic regression fitted with the SAGA solver (liblinear also supports an L1 penalty, but SAGA is the one that scales well). The key hyperparameter is C (the inverse of the regularization strength λ):

Small C → strong regularization → most coefficients shrink to exactly zero → sparse, interpretable model
Large C → weak regularization → many non-zero coefficients → complex, potentially overfitted model

Python

model_pipeline = Pipeline([
    ("preprocessing", ColumnTransformer([
        # Numerical branch: impute then standardise
        ("numeric", Pipeline([
            ("imputer", IterativeImputer(max_iter=10, sample_posterior=False, random_state=RANDOM_STATE)),
            ("scaler", StandardScaler()),
        ]), numerical_features),
        # Categorical branch: fill missing as its own category, then one-hot encode
        ("categorical", Pipeline([
            ("imputer", SimpleImputer(strategy="constant", fill_value="missing")),
            ("encoder", OneHotEncoder(drop="first", handle_unknown="ignore")),
        ]), categorical_features),
    ])),
    # L1 (LASSO) logistic regression; smaller C → stronger sparsity penalty
    ("classifier", LogisticRegression(
        solver="saga",       # Only solver supporting L1 at scale
        penalty="l1",
        max_iter=10_000,     # Generous iteration budget to ensure convergence across the full C grid
        random_state=RANDOM_STATE,
    )),
])

model_pipeline = Pipeline([
    ("preprocessing", ColumnTransformer([
        # Numerical branch: impute then standardise
        ("numeric", Pipeline([
            ("imputer", IterativeImputer(max_iter=10, sample_posterior=False, random_state=RANDOM_STATE)),
            ("scaler", StandardScaler()),
        ]), numerical_features),
        # Categorical branch: fill missing as its own category, then one-hot encode
        ("categorical", Pipeline([
            ("imputer", SimpleImputer(strategy="constant", fill_value="missing")),
            ("encoder", OneHotEncoder(drop="first", handle_unknown="ignore")),
        ]), categorical_features),
    ])),
    # L1 (LASSO) logistic regression; smaller C → stronger sparsity penalty
    ("classifier", LogisticRegression(
        solver="saga",       # Only solver supporting L1 at scale
        penalty="l1",
        max_iter=10_000,     # Generous iteration budget to ensure convergence across the full C grid
        random_state=RANDOM_STATE,
    )),
])

4. Hyperparameter grid and nested cross-validation setup

Regularization grid

We search C over 200 log-spaced values from 10^-2.1 ≈ 0.008 (strong L1, very sparse) to 10^0.9 ≈ 7.9 (weak L1, dense). Log-spacing is essential: the transition from zero to non-zero coefficients happens multiplicatively, so equal spacing on a log scale samples the critical sparsity-performance trade-off region densely.

Scoring: AUPRC

The Area Under the Precision-Recall Curve (AUPRC) is the primary metric. It is most informative when the positive class is rare, because AUROC itself does not change with prevalence, yet when negatives are abundant a large number of false positives barely moves the false positive rate, so a high AUROC can hide poor precision while AUPRC stays sensitive to performance on the positive class. Prevalence here is substantial (~55%), and the class balance does not by itself call for AUPRC. It is the primary metric because the positive class is the one of clinical interest, and its no-skill baseline sits at the prevalence, 0.553, in place of the familiar 0.5.

Nested cross-validation structure

Outer loop: RepeatedStratifiedKFold (100 repeats × 5 folds = 500 outer folds)
│   Provides an unbiased estimate of generalization performance.
│   Repetition over 100 trials yields a stable distribution for computing 95% resampling intervals.
│
└── Inner loop: StratifiedKFold (5 folds, used by GridSearchCV)
        Selects the best C on held-out inner validation data.
        Stratification preserves class balance in every fold.

Nested cross-validation (CV) avoids the optimistic bias that arises when the same data are used for both hyperparameter selection and performance estimation (a common mistake). The inner loop never sees the outer loop’s test data.

Python

# Log-spaced C grid: 200 values from ~0.008 (strong L1, sparse) to ~7.9 (weak L1, dense)
param_grid = {"classifier__C": np.logspace(-2.1, 0.9, 200)}

# AUPRC scorer using predicted probabilities as continuous scores.
# `response_method="predict_proba"` is required in scikit-learn >= 1.4;
# the deprecated `needs_threshold=True` was removed in 1.6 and causes a TypeError.
scoring_fn = make_scorer(average_precision_score, response_method="predict_proba")

# Outer CV: 100 repetitions × 5 folds = 500 outer evaluations for stable resampling intervals
outer_cv = RepeatedStratifiedKFold(n_splits=N_OUTER_SPLITS, n_repeats=N_TRIALS, random_state=RANDOM_STATE)

# Inner CV: standard 5-fold stratified split for GridSearchCV hyperparameter selection
inner_cv = StratifiedKFold(n_splits=N_INNER_SPLITS, shuffle=True, random_state=RANDOM_STATE)

# Log-spaced C grid: 200 values from ~0.008 (strong L1, sparse) to ~7.9 (weak L1, dense)
param_grid = {"classifier__C": np.logspace(-2.1, 0.9, 200)}

# AUPRC scorer using predicted probabilities as continuous scores.
# `response_method="predict_proba"` is required in scikit-learn >= 1.4;
# the deprecated `needs_threshold=True` was removed in 1.6 and causes a TypeError.
scoring_fn = make_scorer(average_precision_score, response_method="predict_proba")

# Outer CV: 100 repetitions × 5 folds = 500 outer evaluations for stable resampling intervals
outer_cv = RepeatedStratifiedKFold(n_splits=N_OUTER_SPLITS, n_repeats=N_TRIALS, random_state=RANDOM_STATE)

# Inner CV: standard 5-fold stratified split for GridSearchCV hyperparameter selection
inner_cv = StratifiedKFold(n_splits=N_INNER_SPLITS, shuffle=True, random_state=RANDOM_STATE)

5. Helper functions for nested cross-validation

Two functions power the nested CV loop:

`average_nonzero_coefficients`

Re-fits the pipeline on each inner training split and counts the coefficients whose absolute value exceeds a small tolerance. Averaging across folds accounts for L1 instability. When predictors are correlated, small perturbations in the data can cause the model to switch between equivalent sparse solutions. The average provides a stable measure of expected model complexity.

`evaluate_outer_fold`

The main workhorse for a single outer fold. It:

Runs GridSearchCV on the inner splits to find the C maximizing AUPRC.
Applies the diminishing-returns selection rule. For each tolerance (90%, 95%, 99%), finds the smallest C (i.e. the sparsest model) whose inner AUPRC meets the threshold.
Refits each candidate model on the full outer training set and generates held-out predictions.
Records the full inner CV score and complexity curve for later aggregation.

The function is designed to be called in parallel via joblib.Parallel, with each outer fold running independently on a separate CPU core.

Python

def average_nonzero_coefficients(pipeline, C_val, X_train, y_train, inner_splits, tolerance=1e-6):
    """
    Estimate the average number of non-zero coefficients for a given C.

    Re-fits the pipeline on each inner training split and counts coefficients
    with |value| > tolerance. Averaging across folds accounts for L1 instability:
    correlated predictors can produce different sparse solutions across folds even
    when overall performance is unchanged.

    Parameters
    ----------
    pipeline : sklearn Pipeline
        Full preprocessing + classifier pipeline (cloned before fitting).
    C_val : float
        Regularisation strength to evaluate.
    X_train : pd.DataFrame
        Outer-fold training features.
    y_train : pd.Series
        Outer-fold training labels.
    inner_splits : list of (train_idx, val_idx)
        Pre-computed inner cross-validation splits (integer index arrays into X_train).
    tolerance : float
        Threshold below which a coefficient is considered zero (default: 1e-6).

    Returns
    -------
    float
        Mean number of non-zero coefficients across inner folds.
    """
    nonzero_counts = []
    for train_idx, _ in inner_splits:
        # Clone the pipeline to avoid state contamination between folds
        model = clone(pipeline).set_params(classifier__C=C_val).fit(
            X_train.iloc[train_idx], y_train.iloc[train_idx]
        )
        nonzero = (np.abs(model.named_steps["classifier"].coef_) > tolerance).sum()
        nonzero_counts.append(nonzero)
    return float(np.mean(nonzero_counts))


def evaluate_outer_fold(train_idx, test_idx, split_id):
    """
    Execute nested CV for one outer fold.

    For a given outer train/test split:
    1. Run inner GridSearchCV to identify the C that maximises AUPRC.
    2. Apply the diminishing-returns rule: for each tolerance band (90%, 95%, 99%),
       select the *smallest* C (sparsest model) whose inner AUPRC exceeds the threshold.
    3. Refit each candidate model on the full outer training set and predict on the test set.
    4. Return inner CV scores and complexity counts for all 200 candidate C values.

    Parameters
    ----------
    train_idx, test_idx : array-like of int
        Indices into the global X/y arrays for the outer fold.
    split_id : int
        Global split index (0 to N_TRIALS * N_OUTER_SPLITS - 1), used to recover
        the trial number and fold number.

    Returns
    -------
    dict with keys: trial, outer_fold, inner_records, y_true,
                    y90, y95, y99, y_best (predicted probability scores).
    """
    # Recover trial and fold indices from the flat split_id
    trial_idx = split_id // N_OUTER_SPLITS
    outer_fold_idx = split_id % N_OUTER_SPLITS

    X_train, X_test = X.iloc[train_idx], X.iloc[test_idx]
    y_train, y_test = y.iloc[train_idx], y.iloc[test_idx]

    # Pre-compute inner splits once; reused by GridSearchCV and non-zero coefficient counting
    inner_splits = list(inner_cv.split(X_train, y_train))

    # Inner hyperparameter search across the 200-point C grid.
    # refit=False: we manually refit with the chosen C on the full outer training set below.
    grid_search = GridSearchCV(
        model_pipeline,
        param_grid,
        scoring=scoring_fn,
        cv=inner_splits,
        n_jobs=1,          # Parallelism is handled at the outer fold level
        refit=False,
        return_train_score=False,
    ).fit(X_train, y_train)

    candidate_Cs = grid_search.cv_results_["param_classifier__C"].astype(float)
    # mean_inner_cv_auprc: mean AUPRC across inner validation folds — model selection signal, NOT outer generalisation
    inner_cv_auprc = grid_search.cv_results_["mean_test_score"]
    max_inner_auprc = inner_cv_auprc.max()
    best_C_overall = candidate_Cs[inner_cv_auprc.argmax()]

    def select_best_C(threshold_ratio):
        """Return the sparsest (smallest) C within threshold_ratio of the peak inner AUPRC."""
        eligible = candidate_Cs[inner_cv_auprc >= threshold_ratio * max_inner_auprc]
        return eligible.min() if eligible.size else best_C_overall

    # Identify the sparsest C for each performance tolerance band
    C_90, C_95, C_99 = select_best_C(0.90), select_best_C(0.95), select_best_C(0.99)

    def predict_with_C(C_val):
        """Refit on the full outer training set with a given C and return test-set probabilities."""
        model = clone(model_pipeline).set_params(classifier__C=C_val).fit(X_train, y_train)
        return model.predict_proba(X_test)[:, 1]

    # Count non-zero coefficients for every C value, averaged over inner folds.
    # This builds the full complexity profile used in the performance–complexity curve.
    nonzero_coefficients = {
        c: average_nonzero_coefficients(model_pipeline, c, X_train, y_train, inner_splits)
        for c in np.unique(candidate_Cs)
    }

    # Assemble one record per C value for later aggregation across trials and folds
    inner_results = [
        {
            "trial": trial_idx,
            "outer_fold": outer_fold_idx,
            "C": c,
            "mean_inner_cv_auprc": s,
            "nonzero": nonzero_coefficients[c],
        }
        for c, s in zip(candidate_Cs, inner_cv_auprc)
    ]

    return {
        "trial": trial_idx,
        "outer_fold": outer_fold_idx,
        "inner_records": inner_results,
        "y_true": y_test.values,
        "y90": predict_with_C(C_90),      # Sparsest model within 10% of max AUPRC (Model D)
        "y95": predict_with_C(C_95),      # Sparsest model within  5% of max AUPRC (Model C)
        "y99": predict_with_C(C_99),      # Sparsest model within  1% of max AUPRC (Model B)
        "y_best": predict_with_C(best_C_overall),  # Globally best C by inner AUPRC  (Model A)
    }

def average_nonzero_coefficients(pipeline, C_val, X_train, y_train, inner_splits, tolerance=1e-6):
    """
    Estimate the average number of non-zero coefficients for a given C.

    Re-fits the pipeline on each inner training split and counts coefficients
    with |value| > tolerance. Averaging across folds accounts for L1 instability:
    correlated predictors can produce different sparse solutions across folds even
    when overall performance is unchanged.

    Parameters
    ----------
    pipeline : sklearn Pipeline
        Full preprocessing + classifier pipeline (cloned before fitting).
    C_val : float
        Regularisation strength to evaluate.
    X_train : pd.DataFrame
        Outer-fold training features.
    y_train : pd.Series
        Outer-fold training labels.
    inner_splits : list of (train_idx, val_idx)
        Pre-computed inner cross-validation splits (integer index arrays into X_train).
    tolerance : float
        Threshold below which a coefficient is considered zero (default: 1e-6).

    Returns
    -------
    float
        Mean number of non-zero coefficients across inner folds.
    """
    nonzero_counts = []
    for train_idx, _ in inner_splits:
        # Clone the pipeline to avoid state contamination between folds
        model = clone(pipeline).set_params(classifier__C=C_val).fit(
            X_train.iloc[train_idx], y_train.iloc[train_idx]
        )
        nonzero = (np.abs(model.named_steps["classifier"].coef_) > tolerance).sum()
        nonzero_counts.append(nonzero)
    return float(np.mean(nonzero_counts))


def evaluate_outer_fold(train_idx, test_idx, split_id):
    """
    Execute nested CV for one outer fold.

    For a given outer train/test split:
    1. Run inner GridSearchCV to identify the C that maximises AUPRC.
    2. Apply the diminishing-returns rule: for each tolerance band (90%, 95%, 99%),
       select the *smallest* C (sparsest model) whose inner AUPRC exceeds the threshold.
    3. Refit each candidate model on the full outer training set and predict on the test set.
    4. Return inner CV scores and complexity counts for all 200 candidate C values.

    Parameters
    ----------
    train_idx, test_idx : array-like of int
        Indices into the global X/y arrays for the outer fold.
    split_id : int
        Global split index (0 to N_TRIALS * N_OUTER_SPLITS - 1), used to recover
        the trial number and fold number.

    Returns
    -------
    dict with keys: trial, outer_fold, inner_records, y_true,
                    y90, y95, y99, y_best (predicted probability scores).
    """
    # Recover trial and fold indices from the flat split_id
    trial_idx = split_id // N_OUTER_SPLITS
    outer_fold_idx = split_id % N_OUTER_SPLITS

    X_train, X_test = X.iloc[train_idx], X.iloc[test_idx]
    y_train, y_test = y.iloc[train_idx], y.iloc[test_idx]

    # Pre-compute inner splits once; reused by GridSearchCV and non-zero coefficient counting
    inner_splits = list(inner_cv.split(X_train, y_train))

    # Inner hyperparameter search across the 200-point C grid.
    # refit=False: we manually refit with the chosen C on the full outer training set below.
    grid_search = GridSearchCV(
        model_pipeline,
        param_grid,
        scoring=scoring_fn,
        cv=inner_splits,
        n_jobs=1,          # Parallelism is handled at the outer fold level
        refit=False,
        return_train_score=False,
    ).fit(X_train, y_train)

    candidate_Cs = grid_search.cv_results_["param_classifier__C"].astype(float)
    # mean_inner_cv_auprc: mean AUPRC across inner validation folds — model selection signal, NOT outer generalisation
    inner_cv_auprc = grid_search.cv_results_["mean_test_score"]
    max_inner_auprc = inner_cv_auprc.max()
    best_C_overall = candidate_Cs[inner_cv_auprc.argmax()]

    def select_best_C(threshold_ratio):
        """Return the sparsest (smallest) C within threshold_ratio of the peak inner AUPRC."""
        eligible = candidate_Cs[inner_cv_auprc >= threshold_ratio * max_inner_auprc]
        return eligible.min() if eligible.size else best_C_overall

    # Identify the sparsest C for each performance tolerance band
    C_90, C_95, C_99 = select_best_C(0.90), select_best_C(0.95), select_best_C(0.99)

    def predict_with_C(C_val):
        """Refit on the full outer training set with a given C and return test-set probabilities."""
        model = clone(model_pipeline).set_params(classifier__C=C_val).fit(X_train, y_train)
        return model.predict_proba(X_test)[:, 1]

    # Count non-zero coefficients for every C value, averaged over inner folds.
    # This builds the full complexity profile used in the performance–complexity curve.
    nonzero_coefficients = {
        c: average_nonzero_coefficients(model_pipeline, c, X_train, y_train, inner_splits)
        for c in np.unique(candidate_Cs)
    }

    # Assemble one record per C value for later aggregation across trials and folds
    inner_results = [
        {
            "trial": trial_idx,
            "outer_fold": outer_fold_idx,
            "C": c,
            "mean_inner_cv_auprc": s,
            "nonzero": nonzero_coefficients[c],
        }
        for c, s in zip(candidate_Cs, inner_cv_auprc)
    ]

    return {
        "trial": trial_idx,
        "outer_fold": outer_fold_idx,
        "inner_records": inner_results,
        "y_true": y_test.values,
        "y90": predict_with_C(C_90),      # Sparsest model within 10% of max AUPRC (Model D)
        "y95": predict_with_C(C_95),      # Sparsest model within  5% of max AUPRC (Model C)
        "y99": predict_with_C(C_99),      # Sparsest model within  1% of max AUPRC (Model B)
        "y_best": predict_with_C(best_C_overall),  # Globally best C by inner AUPRC  (Model A)
    }

6. Run nested cross-validation

We dispatch all 500 outer folds (100 trials × 5 folds) in parallel across all available CPU cores. Each job is independent: it receives its own train/test index arrays and returns predictions and inner CV records without shared mutable state.

What happens inside the pipeline?
Each outer fold call:

Fits 200 × 5 = 1,000 models for the inner GridSearchCV

Fits an additional 200 × 5 = 1,000 models to count non-zero coefficients

Fits 4 final models (one per candidate C: best, 99%, 95%, 90%) on the full outer training set

Across 500 outer folds, this amounts to 1,002,000 model fits. This remains feasible in practice due to the small dataset size, the parallel execution, and SAGA’s efficient convergence.

Python

cv_results = Parallel(n_jobs=-1, verbose=1)(
    delayed(evaluate_outer_fold)(train_idx, test_idx, idx)
    for idx, (train_idx, test_idx) in enumerate(outer_cv.split(X, y))
)

cv_results = Parallel(n_jobs=-1, verbose=1)(
    delayed(evaluate_outer_fold)(train_idx, test_idx, idx)
    for idx, (train_idx, test_idx) in enumerate(outer_cv.split(X, y))
)

[Parallel(n_jobs=-1)]: Using backend LokyBackend with 32 concurrent workers.
[Parallel(n_jobs=-1)]: Done 136 tasks      | elapsed:  4.6min
[Parallel(n_jobs=-1)]: Done 386 tasks      | elapsed: 12.5min
[Parallel(n_jobs=-1)]: Done 500 out of 500 | elapsed: 15.8min finished

7. Aggregate inner CV results

We aggregate the raw inner CV records to construct the performance–complexity curve, the central diagnostic tool for the sparsity trade-off analysis.

Aggregation steps

Trial-level summary: For each (trial, C) pair, average AUPRC and non-zero coefficient count across the 5 outer folds of that trial. This produces one performance estimate per C per trial.
Global summary: Across all 100 trials, compute the mean AUPRC and mean non-zero count for each C, plus 2.5th and 97.5th percentiles to form 95% resampling intervals.
Model selection: Four models are identified from the global summary:
- Model A (best_C_max_auprc): The C with the highest mean AUPRC, the globally best-performing model.
- Model B (best_C_within_1): Sparsest model within 1% of Model A’s AUPRC.
- Model C (best_C_within_5): Sparsest model within 5% of Model A’s AUPRC.
- Model D (best_C_within_10): Sparsest model within 10% of Model A’s AUPRC.

Models B–D represent the diminishing-returns candidates. These are simpler models whose mean inner-CV validation AUPRC stays within 1%, 5%, or 10% of the maximum observed. The rule constrains that validation curve. It does not bound the held-out gap.

Python

# Flatten all inner CV records from all outer folds into a single DataFrame
inner_cv_df = pd.concat(pd.DataFrame(fold["inner_records"]) for fold in cv_results)

# Step 1: Trial-level summary — average over the 5 outer folds within each trial
# This gives one AUPRC and one complexity estimate per (trial, C) combination
trial_summary_df = (
    inner_cv_df.groupby(["trial", "C"])
    .agg(
        mean_inner_cv_auprc=("mean_inner_cv_auprc", "mean"),
        mean_nonzero=("nonzero", "mean"),
    )
    .reset_index()
)

# Step 2: Global summary — aggregate over all 100 trials per C value
# The 2.5/97.5 percentiles are empirical resampling intervals (not classical CIs):
# they describe the range of inner-CV AUPRC estimates across 100 random data splits.
# Percentiles form the 95% resampling interval
global_summary_df = (
    trial_summary_df.groupby("C")
    .agg(
        mean_auprc=("mean_inner_cv_auprc", "mean"),
        auprc_ci_low=("mean_inner_cv_auprc", lambda x: np.percentile(x, 2.5)),
        auprc_ci_high=("mean_inner_cv_auprc", lambda x: np.percentile(x, 97.5)),
        mean_nonzero=("mean_nonzero", "mean"),
        nonzero_ci_low=("mean_nonzero", lambda x: np.percentile(x, 2.5)),
        nonzero_ci_high=("mean_nonzero", lambda x: np.percentile(x, 97.5)),
    )
    .reset_index()
)

# Step 3: Model selection — identify the best C and the sparsest Cs within tolerance bands
best_C_max_auprc = global_summary_df.loc[global_summary_df["mean_auprc"].idxmax(), "C"]

def select_best_C_within(threshold):
    """Select the sparsest C whose mean AUPRC is within `threshold` of the global maximum."""
    filtered = global_summary_df[
        global_summary_df["mean_auprc"] >= threshold * global_summary_df["mean_auprc"].max()
    ]
    return filtered.loc[filtered["mean_nonzero"].idxmin(), "C"]

best_C_within_1  = select_best_C_within(0.99)   # Within 1% of max AUPRC
best_C_within_5  = select_best_C_within(0.95)   # Within 5% of max AUPRC
best_C_within_10 = select_best_C_within(0.90)   # Within 10% of max AUPRC

print(f"Model A — Best AUPRC:           C = {best_C_max_auprc:.4f}")
print(f"Model B — Sparsest within  1%:  C = {best_C_within_1:.4f}")
print(f"Model C — Sparsest within  5%:  C = {best_C_within_5:.4f}")
print(f"Model D — Sparsest within 10%:  C = {best_C_within_10:.4f}")

# Flatten all inner CV records from all outer folds into a single DataFrame
inner_cv_df = pd.concat(pd.DataFrame(fold["inner_records"]) for fold in cv_results)

# Step 1: Trial-level summary — average over the 5 outer folds within each trial
# This gives one AUPRC and one complexity estimate per (trial, C) combination
trial_summary_df = (
    inner_cv_df.groupby(["trial", "C"])
    .agg(
        mean_inner_cv_auprc=("mean_inner_cv_auprc", "mean"),
        mean_nonzero=("nonzero", "mean"),
    )
    .reset_index()
)

# Step 2: Global summary — aggregate over all 100 trials per C value
# The 2.5/97.5 percentiles are empirical resampling intervals (not classical CIs):
# they describe the range of inner-CV AUPRC estimates across 100 random data splits.
# Percentiles form the 95% resampling interval
global_summary_df = (
    trial_summary_df.groupby("C")
    .agg(
        mean_auprc=("mean_inner_cv_auprc", "mean"),
        auprc_ci_low=("mean_inner_cv_auprc", lambda x: np.percentile(x, 2.5)),
        auprc_ci_high=("mean_inner_cv_auprc", lambda x: np.percentile(x, 97.5)),
        mean_nonzero=("mean_nonzero", "mean"),
        nonzero_ci_low=("mean_nonzero", lambda x: np.percentile(x, 2.5)),
        nonzero_ci_high=("mean_nonzero", lambda x: np.percentile(x, 97.5)),
    )
    .reset_index()
)

# Step 3: Model selection — identify the best C and the sparsest Cs within tolerance bands
best_C_max_auprc = global_summary_df.loc[global_summary_df["mean_auprc"].idxmax(), "C"]

def select_best_C_within(threshold):
    """Select the sparsest C whose mean AUPRC is within `threshold` of the global maximum."""
    filtered = global_summary_df[
        global_summary_df["mean_auprc"] >= threshold * global_summary_df["mean_auprc"].max()
    ]
    return filtered.loc[filtered["mean_nonzero"].idxmin(), "C"]

best_C_within_1  = select_best_C_within(0.99)   # Within 1% of max AUPRC
best_C_within_5  = select_best_C_within(0.95)   # Within 5% of max AUPRC
best_C_within_10 = select_best_C_within(0.90)   # Within 10% of max AUPRC

print(f"Model A — Best AUPRC:           C = {best_C_max_auprc:.4f}")
print(f"Model B — Sparsest within  1%:  C = {best_C_within_1:.4f}")
print(f"Model C — Sparsest within  5%:  C = {best_C_within_5:.4f}")
print(f"Model D — Sparsest within 10%:  C = {best_C_within_10:.4f}")

Model A — Best AUPRC:           C = 6.2298
Model B — Sparsest within  1%:  C = 0.1073
Model C — Sparsest within  5%:  C = 0.0241
Model D — Sparsest within 10%:  C = 0.0116

8. Smooth the performance–complexity curves

The raw aggregated curves are noisy due to finite sample variability and the discrete nature of sparsity (the number of non-zero coefficients changes in steps). We apply smoothing splines on the log₁₀(C) scale to produce cleaner visualizations.

Key choices

Log-scale smoothing: Because C spans three orders of magnitude, fitting splines on log₁₀(C) distributes smoothing effort evenly across the range rather than over-smoothing the dense low-C region.
Smoothing factor s: A small value (s = 0.01) for AUPRC preserves the subtle plateau shape. A larger value (s = 1) for the non-zero coefficient count smooths the discrete staircase pattern.
Clipping: AUPRC values are clipped to [0, 1] to prevent the spline from producing physically meaningless overshoot near the boundary.

Interval ordering enforcement

Smoothing the mean and each interval bound independently does not guarantee that the smoothed bounds remain on the correct side of the smoothed mean. If the interval is narrow relative to the spline’s flexibility, the upper bound can dip below the mean or the lower bound can rise above it, producing interval bands that cross the mean curve. After smoothing, we therefore clip each bound to enforce lower ≤ mean ≤ upper pointwise.

Python

def smooth_spline_logx(x, y, s_val):
    """Fit a smoothing spline on log10(x) and evaluate it at the same x positions."""
    return UnivariateSpline(np.log10(x), y, s=s_val)(np.log10(x))

C_values = global_summary_df["C"].to_numpy()

# Smooth AUPRC mean and interval bounds (tight smoothing to preserve the plateau shape)
smooth_auprc       = np.clip(smooth_spline_logx(C_values, global_summary_df["mean_auprc"],    s_val=0.01), 0, 1)
smooth_auprc_upper = np.clip(smooth_spline_logx(C_values, global_summary_df["auprc_ci_high"], s_val=0.01), 0, 1)
smooth_auprc_lower = np.clip(smooth_spline_logx(C_values, global_summary_df["auprc_ci_low"],  s_val=0.01), 0, 1)

# Smooth non-zero coefficient mean and interval bounds (looser smoothing to remove the staircase)
smooth_nonzero       = smooth_spline_logx(C_values, global_summary_df["mean_nonzero"],    s_val=1)
smooth_nonzero_upper = smooth_spline_logx(C_values, global_summary_df["nonzero_ci_high"], s_val=1)
smooth_nonzero_lower = smooth_spline_logx(C_values, global_summary_df["nonzero_ci_low"],  s_val=1)
# Enforce ordering: independently smoothed bounds can cross the mean curve due to spline overfitting.
# Clipping guarantees lower <= mean <= upper everywhere, preserving valid interval geometry.
smooth_auprc_upper = np.clip(smooth_auprc_upper, smooth_auprc, 1.0)
smooth_auprc_lower = np.clip(smooth_auprc_lower, 0.0, smooth_auprc)
smooth_nonzero_upper = np.maximum(smooth_nonzero_upper, smooth_nonzero)
smooth_nonzero_lower = np.minimum(smooth_nonzero_lower, smooth_nonzero)

def smooth_spline_logx(x, y, s_val):
    """Fit a smoothing spline on log10(x) and evaluate it at the same x positions."""
    return UnivariateSpline(np.log10(x), y, s=s_val)(np.log10(x))

C_values = global_summary_df["C"].to_numpy()

# Smooth AUPRC mean and interval bounds (tight smoothing to preserve the plateau shape)
smooth_auprc       = np.clip(smooth_spline_logx(C_values, global_summary_df["mean_auprc"],    s_val=0.01), 0, 1)
smooth_auprc_upper = np.clip(smooth_spline_logx(C_values, global_summary_df["auprc_ci_high"], s_val=0.01), 0, 1)
smooth_auprc_lower = np.clip(smooth_spline_logx(C_values, global_summary_df["auprc_ci_low"],  s_val=0.01), 0, 1)

# Smooth non-zero coefficient mean and interval bounds (looser smoothing to remove the staircase)
smooth_nonzero       = smooth_spline_logx(C_values, global_summary_df["mean_nonzero"],    s_val=1)
smooth_nonzero_upper = smooth_spline_logx(C_values, global_summary_df["nonzero_ci_high"], s_val=1)
smooth_nonzero_lower = smooth_spline_logx(C_values, global_summary_df["nonzero_ci_low"],  s_val=1)
# Enforce ordering: independently smoothed bounds can cross the mean curve due to spline overfitting.
# Clipping guarantees lower <= mean <= upper everywhere, preserving valid interval geometry.
smooth_auprc_upper = np.clip(smooth_auprc_upper, smooth_auprc, 1.0)
smooth_auprc_lower = np.clip(smooth_auprc_lower, 0.0, smooth_auprc)
smooth_nonzero_upper = np.maximum(smooth_nonzero_upper, smooth_nonzero)
smooth_nonzero_lower = np.minimum(smooth_nonzero_lower, smooth_nonzero)

9. Hyperparameter selection plot

This plot is the core visualization of the performance–complexity trade-off.

Left axis (blue): Mean inner-CV AUPRC (± 95% empirical interval) as a function of C. These scores come from the inner validation folds of the nested CV (models trained on ~64% of the data, evaluated on ~16%). They are the model-selection signal used to choose C within each outer fold, and are slightly pessimistic compared to models trained on the full outer training set (~80% of the data).
Right axis (purple): Mean number of non-zero coefficients (± 95% empirical interval). This rises with C as more features are admitted, though not strictly, since a feature can leave the L1 solution as the penalty weakens, which is part of why the raw count curve is a noisy staircase.
Empirical 95% interval: The shaded bands show the 2.5th–97.5th percentile range of the 100 trial-level estimates. This is a resampling interval. It describes where 95% of inner-CV estimates fall across different random data splits. A classical confidence interval around the population parameter is a different object.
Annotated models (A–D): Each marker shows where the four candidate models sit on the curve. Models B–D sit where the tolerance rule places them: C and D are by construction the sparsest fits within 5% and 10% of the maximum validation score, so their position follows from the tolerance, which is no demonstration that sparsity is free. The held-out scores confirm the cost is real for C and D.

What to look for:

The plateau onset: the C at which AUPRC levels off and additional features contribute negligible gain.
The vertical gap between Model A and Model D: the performance cost of moving from densest to sparsest.
The width of the interval band: a wide band indicates high variability across splits. Small AUPRC differences between candidate models may not be meaningful.

Python

# Truncate all plotted curves to the region where the mean number of retained features is at least 1
mask = smooth_nonzero >= 1

C_values_trunc = C_values[mask]
smooth_auprc_trunc = smooth_auprc[mask]
smooth_auprc_upper_trunc = smooth_auprc_upper[mask]
smooth_auprc_lower_trunc = smooth_auprc_lower[mask]

smooth_nonzero_trunc = smooth_nonzero[mask]
smooth_nonzero_upper_trunc = smooth_nonzero_upper[mask]
smooth_nonzero_lower_trunc = smooth_nonzero_lower[mask]

fig_hp = go.Figure()

# Plot AUPRC curve
fig_hp.add_trace(go.Scatter(
    x=C_values_trunc,
    y=smooth_auprc_trunc,
    mode="lines",
    name="Mean AUPRC",
    yaxis="y1",
    hovertemplate="<b>C</b>: %{x:.3f}<br><b>Mean AUPRC</b>: %{y:.3f}<extra></extra>"
))
fig_hp.add_trace(go.Scatter(
    x=np.r_[C_values_trunc, C_values_trunc[::-1]],
    y=np.r_[smooth_auprc_upper_trunc, smooth_auprc_lower_trunc[::-1]],
    fill="toself",
    fillcolor="rgba(0,0,255,0.12)",
    line=dict(color="rgba(0,0,0,0)"),
    hoverinfo="skip",
    name="95% CI (AUPRC)",
    yaxis="y1"
))

# Plot non-zero coefficients curve
fig_hp.add_trace(go.Scatter(
    x=C_values_trunc,
    y=smooth_nonzero_trunc,
    mode="lines",
    name="Mean N of Features Retained",
    yaxis="y2",
    line=dict(color="purple", dash="dot"),
    hovertemplate="<b>C</b>: %{x:.3f}<br><b>Features</b>: %{y:.3f}<extra></extra>"
))
fig_hp.add_trace(go.Scatter(
    x=np.r_[C_values_trunc, C_values_trunc[::-1]],
    y=np.r_[smooth_nonzero_upper_trunc, smooth_nonzero_lower_trunc[::-1]],
    fill="toself",
    fillcolor="rgba(128,0,128,0.12)",
    line=dict(color="rgba(0,0,0,0)"),
    hoverinfo="skip",
    name="95% CI (N of Features)",
    yaxis="y2"
))

# Annotate selected C values on both curves
selected_Cs = [
    (best_C_max_auprc, "darkorange",   "Model A"),
    (best_C_within_1,  "forestgreen",  "Model B"),
    (best_C_within_5,  "mediumorchid", "Model C"),
    (best_C_within_10, "steelblue",    "Model D"),
]

for C_val, color, label in selected_Cs:
    if len(C_values_trunc) > 0 and C_values_trunc.min() <= C_val <= C_values_trunc.max():
        idx = np.argmin(np.abs(C_values_trunc - C_val))

        # Marker on AUPRC curve
        fig_hp.add_trace(go.Scatter(
            x=[C_values_trunc[idx]],
            y=[smooth_auprc_trunc[idx]],
            mode="markers+text",
            marker=dict(size=12, color=color),
            text=[label],
            textposition="top center",
            name=label,
            yaxis="y1",
            hovertemplate="<b>C</b>: %{x:.3f}<br><b>AUPRC</b>: %{y:.3f}<extra></extra>",
            showlegend=True
        ))

        # Marker on purple feature-retention curve (same circle style)
        fig_hp.add_trace(go.Scatter(
            x=[C_values_trunc[idx]],
            y=[smooth_nonzero_trunc[idx]],
            mode="markers+text",
            marker=dict(size=12, color=color),
            text=[label],
            textposition="bottom center",
            name=f"{label} (Features)",
            yaxis="y2",
            hovertemplate="<b>C</b>: %{x:.3f}<br><b>Features</b>: %{y:.3f}<extra></extra>",
            showlegend=False
        ))

fig_hp.update_layout(
    title=dict(text="<b>LASSO Logistic Regression Hyperparameter Selection</b>", x=0.5),
    xaxis=dict(
        type="log",
        showgrid=False,
        title=dict(text="<b>C (1/λ)</b>", font=dict(size=14)),
    ),
    yaxis=dict(
        range=[0, 1],
        showgrid=False,
        title=dict(text="<b>Mean AUPRC</b>", font=dict(color="blue")),
        tickfont=dict(color="blue"),
    ),
    yaxis2=dict(
        range=[0, 17],
        overlaying="y",
        side="right",
        showgrid=False,
        title=dict(text="<b>N of Features Retained</b>", font=dict(color="purple")),
        tickfont=dict(color="purple"),
    ),
    legend=dict(
        orientation="h",
        yanchor="bottom",
        y=-0.35,
        xanchor="center",
        x=0.5,
        bgcolor="rgba(255,255,255,0.95)",
        bordercolor="LightGray",
        borderwidth=1
    ),
    template="plotly_white",
    width=900,
    height=520,
    margin=dict(l=60, r=60, t=60, b=10),
    font=dict(color="black"),
)

fig_hp.write_html("Hyperparameter_Selection_nt.html")
fig_hp.show()

# Truncate all plotted curves to the region where the mean number of retained features is at least 1
mask = smooth_nonzero >= 1

C_values_trunc = C_values[mask]
smooth_auprc_trunc = smooth_auprc[mask]
smooth_auprc_upper_trunc = smooth_auprc_upper[mask]
smooth_auprc_lower_trunc = smooth_auprc_lower[mask]

smooth_nonzero_trunc = smooth_nonzero[mask]
smooth_nonzero_upper_trunc = smooth_nonzero_upper[mask]
smooth_nonzero_lower_trunc = smooth_nonzero_lower[mask]

fig_hp = go.Figure()

# Plot AUPRC curve
fig_hp.add_trace(go.Scatter(
    x=C_values_trunc,
    y=smooth_auprc_trunc,
    mode="lines",
    name="Mean AUPRC",
    yaxis="y1",
    hovertemplate="<b>C</b>: %{x:.3f}<br><b>Mean AUPRC</b>: %{y:.3f}<extra></extra>"
))
fig_hp.add_trace(go.Scatter(
    x=np.r_[C_values_trunc, C_values_trunc[::-1]],
    y=np.r_[smooth_auprc_upper_trunc, smooth_auprc_lower_trunc[::-1]],
    fill="toself",
    fillcolor="rgba(0,0,255,0.12)",
    line=dict(color="rgba(0,0,0,0)"),
    hoverinfo="skip",
    name="95% CI (AUPRC)",
    yaxis="y1"
))

# Plot non-zero coefficients curve
fig_hp.add_trace(go.Scatter(
    x=C_values_trunc,
    y=smooth_nonzero_trunc,
    mode="lines",
    name="Mean N of Features Retained",
    yaxis="y2",
    line=dict(color="purple", dash="dot"),
    hovertemplate="<b>C</b>: %{x:.3f}<br><b>Features</b>: %{y:.3f}<extra></extra>"
))
fig_hp.add_trace(go.Scatter(
    x=np.r_[C_values_trunc, C_values_trunc[::-1]],
    y=np.r_[smooth_nonzero_upper_trunc, smooth_nonzero_lower_trunc[::-1]],
    fill="toself",
    fillcolor="rgba(128,0,128,0.12)",
    line=dict(color="rgba(0,0,0,0)"),
    hoverinfo="skip",
    name="95% CI (N of Features)",
    yaxis="y2"
))

# Annotate selected C values on both curves
selected_Cs = [
    (best_C_max_auprc, "darkorange",   "Model A"),
    (best_C_within_1,  "forestgreen",  "Model B"),
    (best_C_within_5,  "mediumorchid", "Model C"),
    (best_C_within_10, "steelblue",    "Model D"),
]

for C_val, color, label in selected_Cs:
    if len(C_values_trunc) > 0 and C_values_trunc.min() <= C_val <= C_values_trunc.max():
        idx = np.argmin(np.abs(C_values_trunc - C_val))

        # Marker on AUPRC curve
        fig_hp.add_trace(go.Scatter(
            x=[C_values_trunc[idx]],
            y=[smooth_auprc_trunc[idx]],
            mode="markers+text",
            marker=dict(size=12, color=color),
            text=[label],
            textposition="top center",
            name=label,
            yaxis="y1",
            hovertemplate="<b>C</b>: %{x:.3f}<br><b>AUPRC</b>: %{y:.3f}<extra></extra>",
            showlegend=True
        ))

        # Marker on purple feature-retention curve (same circle style)
        fig_hp.add_trace(go.Scatter(
            x=[C_values_trunc[idx]],
            y=[smooth_nonzero_trunc[idx]],
            mode="markers+text",
            marker=dict(size=12, color=color),
            text=[label],
            textposition="bottom center",
            name=f"{label} (Features)",
            yaxis="y2",
            hovertemplate="<b>C</b>: %{x:.3f}<br><b>Features</b>: %{y:.3f}<extra></extra>",
            showlegend=False
        ))

fig_hp.update_layout(
    title=dict(text="<b>LASSO Logistic Regression Hyperparameter Selection</b>", x=0.5),
    xaxis=dict(
        type="log",
        showgrid=False,
        title=dict(text="<b>C (1/λ)</b>", font=dict(size=14)),
    ),
    yaxis=dict(
        range=[0, 1],
        showgrid=False,
        title=dict(text="<b>Mean AUPRC</b>", font=dict(color="blue")),
        tickfont=dict(color="blue"),
    ),
    yaxis2=dict(
        range=[0, 17],
        overlaying="y",
        side="right",
        showgrid=False,
        title=dict(text="<b>N of Features Retained</b>", font=dict(color="purple")),
        tickfont=dict(color="purple"),
    ),
    legend=dict(
        orientation="h",
        yanchor="bottom",
        y=-0.35,
        xanchor="center",
        x=0.5,
        bgcolor="rgba(255,255,255,0.95)",
        bordercolor="LightGray",
        borderwidth=1
    ),
    template="plotly_white",
    width=900,
    height=520,
    margin=dict(l=60, r=60, t=60, b=10),
    font=dict(color="black"),
)

fig_hp.write_html("Hyperparameter_Selection_nt.html")
fig_hp.show()

Figure 4

AUPRC and complexity across C

Open this figure at full size in a new tab

Figure 4. Mean inner-fold AUPRC in blue on the left axis, mean count of non-zero coefficients in purple on the right, both against C on a log scale. Bands are 2.5th to 97.5th percentiles over the 100 trials. Every line is a smoothing spline on log C. Markers fix Models A to D on both curves. AUPRC climbs out of the sparse region and flattens, peaking at C = 6.230 (Model A) while the purple curve keeps rising. Models C and D are the sparsest fits within 5% and 10% of the peak, and both land below the plateau. The right axis counts one-hot encoded coefficients, which outnumber the 11 raw features.

10. Evaluate model performance: PRC and ROC curves

We now assess the out-of-sample performance of all four model variants using the held-out outer fold predictions accumulated over 100 trials.

How the curves are computed

For each trial, the five outer fold predictions are concatenated to form a single trial-level score vector. This gives one PRC and one ROC curve per trial, per model variant. The mean curve and 95% resampling interval are then computed across all 100 trials.

Precision-Recall Curve (PRC): Plots precision vs recall as the classification threshold varies. The AUPRC summarizes the area under this curve. A monotone correction is applied after interpolation onto a common recall grid, replacing precision by its running maximum at higher recall. Real precision-recall curves are not monotone, so this is the standard interpolated-precision convention for averaging curves across trials, imposed by the averaging itself.

ROC Curve: Plots the True Positive Rate (sensitivity) vs False Positive Rate (1 − specificity). The AUROC is reported as a secondary metric.

Why both metrics?

AUPRC is the primary metric because it focuses on the positive class, the class of clinical interest. Its advantage over AUROC is largest when positives are rare. With substantial positive class prevalence (~55%) here the class balance does not by itself call for it, and its no-skill baseline sits at the prevalence, 0.553, in place of the familiar 0.5. AUROC is included for completeness and comparability with the broader literature. In highly imbalanced datasets, AUROC can give an overly optimistic impression. Here, the two metrics are broadly consistent.

Python

# Common interpolation grids for aligning curves across trials
FPR_GRID    = np.linspace(0, 1, 100)
RECALL_GRID = np.linspace(0, 1, 100)

# Map each model variant to the key in the cv_results dict that holds its predictions
MODEL_VARIANTS = {
    "Model A": "y_best",  # Max AUPRC (no tolerance constraint)
    "Model B": "y99",     # Sparsest within 1% of max AUPRC
    "Model C": "y95",     # Sparsest within 5% of max AUPRC
    "Model D": "y90",     # Sparsest within 10% of max AUPRC
}
MODEL_COLORS = {
    "Model A": ("darkorange",   "rgba(255,165,0,0.12)"),
    "Model B": ("forestgreen",  "rgba(34,139,34,0.12)"),
    "Model C": ("mediumorchid", "rgba(186,85,211,0.12)"),
    "Model D": ("steelblue",    "rgba(70,130,180,0.12)"),
}
MODEL_LINESTYLES = {
    "Model A": "solid",
    "Model B": "dash",
    "Model C": "dot",
    "Model D": "dashdot",
}

# Group predictions by model variant and trial
trial_predictions = {model: defaultdict(list) for model in MODEL_VARIANTS}
for result in cv_results:
    trial = result["trial"]
    for model, key in MODEL_VARIANTS.items():
        trial_predictions[model][trial].append((result["y_true"], result[key]))

# Compute mean curves and 95% resampling intervals across trials for each model variant
model_metrics = {}
for model in MODEL_VARIANTS:
    tpr_trials       = []
    precision_trials = []
    auroc_scores     = []
    auprc_scores     = []

    for trial in range(N_TRIALS):
        # Concatenate all 5 outer folds for this trial into a single evaluation set
        y_true_all  = np.concatenate([yt for yt, _  in trial_predictions[model][trial]])
        y_score_all = np.concatenate([yp for _,  yp in trial_predictions[model][trial]])

        # ROC: interpolate onto the common FPR grid so curves can be averaged
        fpr, tpr, _ = roc_curve(y_true_all, y_score_all)
        tpr_trials.append(np.interp(FPR_GRID, fpr, tpr))

        # PRC: interpolate onto the common recall grid.
        # Sort by recall before interpolating (precision_recall_curve returns descending recall).
        # Apply the interpolated-precision convention (running max) so curves can be averaged
        # across trials after interpolation onto the common recall grid.
        prec, rec, _ = precision_recall_curve(y_true_all, y_score_all)
        order        = np.lexsort((prec, rec))  # tie-break: ascending prec -> last = max prec
        interp_prec  = np.interp(RECALL_GRID, rec[order], prec[order])
        interp_prec  = np.minimum(1.0, np.maximum.accumulate(interp_prec[::-1])[::-1])
        precision_trials.append(interp_prec)

        auroc_scores.append(roc_auc_score(y_true_all, y_score_all))
        auprc_scores.append(average_precision_score(y_true_all, y_score_all))

    tpr_trials       = np.stack(tpr_trials)
    precision_trials = np.stack(precision_trials)

    model_metrics[model] = {
        "mean_tpr":      tpr_trials.mean(axis=0),
        "tpr_ci_low":    np.percentile(tpr_trials, 2.5,  axis=0),
        "tpr_ci_high":   np.percentile(tpr_trials, 97.5, axis=0),
        "mean_prec":     precision_trials.mean(axis=0),
        "prec_ci_low":   np.percentile(precision_trials, 2.5,  axis=0),
        "prec_ci_high":  np.percentile(precision_trials, 97.5, axis=0),
        "auroc_mean":    np.mean(auroc_scores),
        "auroc_ci_low":  np.percentile(auroc_scores, 2.5),
        "auroc_ci_high": np.percentile(auroc_scores, 97.5),
        "auprc_mean":    np.mean(auprc_scores),
        "auprc_ci_low":  np.percentile(auprc_scores, 2.5),
        "auprc_ci_high": np.percentile(auprc_scores, 97.5),
    }

# Print a summary table for quick reference
print(f"{'Model':<10} {'AUPRC':<30} {'AUROC'}")
print("-" * 65)
for model, m in model_metrics.items():
    auprc_str = f"{m['auprc_mean']:.3f} [{m['auprc_ci_low']:.3f}–{m['auprc_ci_high']:.3f}]"
    auroc_str = f"{m['auroc_mean']:.3f} [{m['auroc_ci_low']:.3f}–{m['auroc_ci_high']:.3f}]"
    print(f"{model:<10} {auprc_str:<30} {auroc_str}")

# Common interpolation grids for aligning curves across trials
FPR_GRID    = np.linspace(0, 1, 100)
RECALL_GRID = np.linspace(0, 1, 100)

# Map each model variant to the key in the cv_results dict that holds its predictions
MODEL_VARIANTS = {
    "Model A": "y_best",  # Max AUPRC (no tolerance constraint)
    "Model B": "y99",     # Sparsest within 1% of max AUPRC
    "Model C": "y95",     # Sparsest within 5% of max AUPRC
    "Model D": "y90",     # Sparsest within 10% of max AUPRC
}
MODEL_COLORS = {
    "Model A": ("darkorange",   "rgba(255,165,0,0.12)"),
    "Model B": ("forestgreen",  "rgba(34,139,34,0.12)"),
    "Model C": ("mediumorchid", "rgba(186,85,211,0.12)"),
    "Model D": ("steelblue",    "rgba(70,130,180,0.12)"),
}
MODEL_LINESTYLES = {
    "Model A": "solid",
    "Model B": "dash",
    "Model C": "dot",
    "Model D": "dashdot",
}

# Group predictions by model variant and trial
trial_predictions = {model: defaultdict(list) for model in MODEL_VARIANTS}
for result in cv_results:
    trial = result["trial"]
    for model, key in MODEL_VARIANTS.items():
        trial_predictions[model][trial].append((result["y_true"], result[key]))

# Compute mean curves and 95% resampling intervals across trials for each model variant
model_metrics = {}
for model in MODEL_VARIANTS:
    tpr_trials       = []
    precision_trials = []
    auroc_scores     = []
    auprc_scores     = []

    for trial in range(N_TRIALS):
        # Concatenate all 5 outer folds for this trial into a single evaluation set
        y_true_all  = np.concatenate([yt for yt, _  in trial_predictions[model][trial]])
        y_score_all = np.concatenate([yp for _,  yp in trial_predictions[model][trial]])

        # ROC: interpolate onto the common FPR grid so curves can be averaged
        fpr, tpr, _ = roc_curve(y_true_all, y_score_all)
        tpr_trials.append(np.interp(FPR_GRID, fpr, tpr))

        # PRC: interpolate onto the common recall grid.
        # Sort by recall before interpolating (precision_recall_curve returns descending recall).
        # Apply the interpolated-precision convention (running max) so curves can be averaged
        # across trials after interpolation onto the common recall grid.
        prec, rec, _ = precision_recall_curve(y_true_all, y_score_all)
        order        = np.lexsort((prec, rec))  # tie-break: ascending prec -> last = max prec
        interp_prec  = np.interp(RECALL_GRID, rec[order], prec[order])
        interp_prec  = np.minimum(1.0, np.maximum.accumulate(interp_prec[::-1])[::-1])
        precision_trials.append(interp_prec)

        auroc_scores.append(roc_auc_score(y_true_all, y_score_all))
        auprc_scores.append(average_precision_score(y_true_all, y_score_all))

    tpr_trials       = np.stack(tpr_trials)
    precision_trials = np.stack(precision_trials)

    model_metrics[model] = {
        "mean_tpr":      tpr_trials.mean(axis=0),
        "tpr_ci_low":    np.percentile(tpr_trials, 2.5,  axis=0),
        "tpr_ci_high":   np.percentile(tpr_trials, 97.5, axis=0),
        "mean_prec":     precision_trials.mean(axis=0),
        "prec_ci_low":   np.percentile(precision_trials, 2.5,  axis=0),
        "prec_ci_high":  np.percentile(precision_trials, 97.5, axis=0),
        "auroc_mean":    np.mean(auroc_scores),
        "auroc_ci_low":  np.percentile(auroc_scores, 2.5),
        "auroc_ci_high": np.percentile(auroc_scores, 97.5),
        "auprc_mean":    np.mean(auprc_scores),
        "auprc_ci_low":  np.percentile(auprc_scores, 2.5),
        "auprc_ci_high": np.percentile(auprc_scores, 97.5),
    }

# Print a summary table for quick reference
print(f"{'Model':<10} {'AUPRC':<30} {'AUROC'}")
print("-" * 65)
for model, m in model_metrics.items():
    auprc_str = f"{m['auprc_mean']:.3f} [{m['auprc_ci_low']:.3f}–{m['auprc_ci_high']:.3f}]"
    auroc_str = f"{m['auroc_mean']:.3f} [{m['auroc_ci_low']:.3f}–{m['auroc_ci_high']:.3f}]"
    print(f"{model:<10} {auprc_str:<30} {auroc_str}")

Model      AUPRC                          AUROC
-----------------------------------------------------------------
Model A    0.931 [0.925–0.935]            0.924 [0.920–0.926]
Model B    0.925 [0.919–0.928]            0.919 [0.915–0.921]
Model C    0.891 [0.882–0.898]            0.881 [0.874–0.886]
Model D    0.859 [0.848–0.872]            0.836 [0.819–0.852]

11. Plot PRC and ROC curves

Below: mean PRC and ROC curves for all four model variants, with shaded 95% percentile bands.

PRC (left): the dotted gray line marks the no-skill baseline, where a classifier just predicts the prevalence every time. Anything sitting well above it is discriminating for real.
ROC (right): the dashed diagonal is the random-classifier baseline. On this dataset every model should sit well clear of it.

Models A and B overlap. C and D sit below both, and that gap is what their extra sparsity costs.

Python

fig = make_subplots(
    rows=1, cols=2,
    subplot_titles=("<b>Precision-Recall Curves (PRC)</b>", "<b>ROC Curves</b>"),
)

# Plot PRC curves
for model in MODEL_VARIANTS:
    color, fill = MODEL_COLORS[model]
    linestyle = MODEL_LINESTYLES[model]
    stats = model_metrics[model]
    label = f"{model}  (AUPRC {stats['auprc_mean']:.3f} [{stats['auprc_ci_low']:.3f}–{stats['auprc_ci_high']:.3f}])"

    fig.add_trace(go.Scatter(
        x=RECALL_GRID, y=stats["mean_prec"], mode="lines", name=label,
        line=dict(color=color, dash=linestyle), legendgroup="PRC",
        hovertemplate="<b>Recall</b>: %{x:.3f}<br><b>Precision</b>: %{y:.3f}<extra></extra>",
    ), row=1, col=1)
    fig.add_trace(go.Scatter(
        x=np.r_[RECALL_GRID, RECALL_GRID[::-1]],
        y=np.r_[stats["prec_ci_high"], stats["prec_ci_low"][::-1]],
        fill="toself", fillcolor=fill, line=dict(color="rgba(0,0,0,0)"),
        hoverinfo="skip", showlegend=False, legendgroup="PRC"
    ), row=1, col=1)

# Plot ROC curves
for model in MODEL_VARIANTS:
    color, fill = MODEL_COLORS[model]
    linestyle = MODEL_LINESTYLES[model]
    stats = model_metrics[model]
    label = f"{model}  (AUROC {stats['auroc_mean']:.3f} [{stats['auroc_ci_low']:.3f}–{stats['auroc_ci_high']:.3f}])"

    fig.add_trace(go.Scatter(
        x=FPR_GRID, y=stats["mean_tpr"], mode="lines", name=label,
        line=dict(color=color, dash=linestyle), legendgroup="ROC",
        hovertemplate="<b>FPR</b>: %{x:.3f}<br><b>TPR</b>: %{y:.3f}<extra></extra>",
    ), row=1, col=2)
    fig.add_trace(go.Scatter(
        x=np.r_[FPR_GRID, FPR_GRID[::-1]],
        y=np.r_[stats["tpr_ci_high"], stats["tpr_ci_low"][::-1]],
        fill="toself", fillcolor=fill, line=dict(color="rgba(0,0,0,0)"),
        hoverinfo="skip", showlegend=False, legendgroup="ROC"
    ), row=1, col=2)

# Add reference lines
fig.add_trace(go.Scatter(x=[0, 1], y=[0, 1], mode="lines",
    line=dict(color="gray", dash="dash"), name="Reference", legendgroup="ROC", hoverinfo="skip"), row=1, col=2)
fig.add_trace(go.Scatter(x=[0, 1], y=[y.mean()] * 2, mode="lines",
    line=dict(color="gray", dash="dot"), name="Prevalence", legendgroup="PRC", hoverinfo="skip"), row=1, col=1)

# Final layout
fig.update_layout(
    title=dict(text="<b>Model Performance Comparison</b>", x=0.5),
    template="plotly_white", width=1200, height=700,
    font=dict(color="black"),
    legend=dict(orientation="h", x=0.5, y=-0.2, xanchor="center", yanchor="top",
                bgcolor="rgba(255,255,255,0.95)", bordercolor="LightGray", borderwidth=1),
)
fig.update_xaxes(title_text="<b>Recall</b>", row=1, col=1)
fig.update_yaxes(title_text="<b>Precision</b>", row=1, col=1)
fig.update_xaxes(title_text="<b>False Positive Rate</b>", row=1, col=2)
fig.update_yaxes(title_text="<b>True Positive Rate</b>", row=1, col=2)

fig.write_html("PRC_ROC_Comparison_nt.html")
fig.show()


fig = make_subplots(
    rows=1, cols=2,
    subplot_titles=("<b>Precision-Recall Curves (PRC)</b>", "<b>ROC Curves</b>"),
)

# Plot PRC curves
for model in MODEL_VARIANTS:
    color, fill = MODEL_COLORS[model]
    linestyle = MODEL_LINESTYLES[model]
    stats = model_metrics[model]
    label = f"{model}  (AUPRC {stats['auprc_mean']:.3f} [{stats['auprc_ci_low']:.3f}–{stats['auprc_ci_high']:.3f}])"

    fig.add_trace(go.Scatter(
        x=RECALL_GRID, y=stats["mean_prec"], mode="lines", name=label,
        line=dict(color=color, dash=linestyle), legendgroup="PRC",
        hovertemplate="<b>Recall</b>: %{x:.3f}<br><b>Precision</b>: %{y:.3f}<extra></extra>",
    ), row=1, col=1)
    fig.add_trace(go.Scatter(
        x=np.r_[RECALL_GRID, RECALL_GRID[::-1]],
        y=np.r_[stats["prec_ci_high"], stats["prec_ci_low"][::-1]],
        fill="toself", fillcolor=fill, line=dict(color="rgba(0,0,0,0)"),
        hoverinfo="skip", showlegend=False, legendgroup="PRC"
    ), row=1, col=1)

# Plot ROC curves
for model in MODEL_VARIANTS:
    color, fill = MODEL_COLORS[model]
    linestyle = MODEL_LINESTYLES[model]
    stats = model_metrics[model]
    label = f"{model}  (AUROC {stats['auroc_mean']:.3f} [{stats['auroc_ci_low']:.3f}–{stats['auroc_ci_high']:.3f}])"

    fig.add_trace(go.Scatter(
        x=FPR_GRID, y=stats["mean_tpr"], mode="lines", name=label,
        line=dict(color=color, dash=linestyle), legendgroup="ROC",
        hovertemplate="<b>FPR</b>: %{x:.3f}<br><b>TPR</b>: %{y:.3f}<extra></extra>",
    ), row=1, col=2)
    fig.add_trace(go.Scatter(
        x=np.r_[FPR_GRID, FPR_GRID[::-1]],
        y=np.r_[stats["tpr_ci_high"], stats["tpr_ci_low"][::-1]],
        fill="toself", fillcolor=fill, line=dict(color="rgba(0,0,0,0)"),
        hoverinfo="skip", showlegend=False, legendgroup="ROC"
    ), row=1, col=2)

# Add reference lines
fig.add_trace(go.Scatter(x=[0, 1], y=[0, 1], mode="lines",
    line=dict(color="gray", dash="dash"), name="Reference", legendgroup="ROC", hoverinfo="skip"), row=1, col=2)
fig.add_trace(go.Scatter(x=[0, 1], y=[y.mean()] * 2, mode="lines",
    line=dict(color="gray", dash="dot"), name="Prevalence", legendgroup="PRC", hoverinfo="skip"), row=1, col=1)

# Final layout
fig.update_layout(
    title=dict(text="<b>Model Performance Comparison</b>", x=0.5),
    template="plotly_white", width=1200, height=700,
    font=dict(color="black"),
    legend=dict(orientation="h", x=0.5, y=-0.2, xanchor="center", yanchor="top",
                bgcolor="rgba(255,255,255,0.95)", bordercolor="LightGray", borderwidth=1),
)
fig.update_xaxes(title_text="<b>Recall</b>", row=1, col=1)
fig.update_yaxes(title_text="<b>Precision</b>", row=1, col=1)
fig.update_xaxes(title_text="<b>False Positive Rate</b>", row=1, col=2)
fig.update_yaxes(title_text="<b>True Positive Rate</b>", row=1, col=2)

fig.write_html("PRC_ROC_Comparison_nt.html")
fig.show()

Figure 5

Held-out PRC and ROC by model

Open this figure at full size in a new tab

Figure 5. Two panels built from the outer-fold predictions pooled within each of the 100 trials, showing precision against recall on the left, with the dotted gray line at the 55.3% prevalence baseline, and true against false positive rate on the right, with the diagonal chance line. Every model carries a mean curve and a 2.5th to 97.5th percentile band. Model A reaches AUPRC 0.931 [0.925 to 0.935] against a no-skill floor of 0.553, and Model B overlaps it. Model C at 0.891 [0.882 to 0.898] and Model D at 0.859 [0.848 to 0.872] fall entirely below both, so Models A and B keep a measurable performance advantage over the sparsest fits.

Interpreting the PRC and ROC results

What to look for:

Curve overlap: Compare the shaded bands across both panels. Models A and B overlap, so these bands do not separate the two, but Models C and D sit below both over most of the curve and their AUPRC intervals fall entirely below, and the extra sparsity does cost measurable performance. Overlapping bands are consistent with equivalence without establishing it. Every trial scores all four models on the same splits, so a paired comparison is what separates A from B.
AUPRC vs. AUROC ordering: The tolerance rule on the validation AUPRC imposes the ordering (A ≥ B ≥ C ≥ D), but it need not survive on held-out data, and AUPRC and AUROC can rank the same models differently, making a reversal informative when it happens. Larger gaps in AUPRC than AUROC suggest that the models differ more in their handling of the positive class, the class of clinical interest.
Interval width vs. sample size: The interval width reflects both the intrinsic variability of the task (difficulty of the prediction problem) and the sample size. Narrow resampling intervals across 100 trials indicate that the performance estimates are stable and not artefacts of a particular train/test split.

Summary and conclusions

This notebook demonstrates a principled approach to the performance–complexity trade-off in LASSO logistic regression:

Model	Constraint	Implication
A	Maximum AUPRC	Best raw performance; most complex
B	Sparsest within 1%	Near-identical performance; slightly simpler
C	Sparsest within 5%	Small performance cost; meaningfully simpler
D	Sparsest within 10%	Modest performance cost; substantially simpler

Key takeaways aligned with the article:

The performance–complexity curve shows a rapid rise followed by a plateau. Most of the predictive value is captured by a small subset of features.
Diminishing returns allow the selection of sparser, more interpretable models, though only Model B has a band that overlaps Model A; Models C and D trade 0.040 and 0.072 AUPRC for their extra sparsity.
Nested cross-validation with repetition is essential for reliable performance estimation and for constructing resampling intervals that show how much the estimates move across random splits.
In clinical applications, the choice between Models A–D depends on context, since interpretability, data availability, and stakeholder trust often favor simpler models even when the performance difference is measurable.

Summary

Key Takeaways

L1 zeroes coefficients, L2 only shrinks them

The L1 penalty sets coefficients to exactly zero, so selection happens while the model is fitted. L2 shrinks coefficients toward zero without ever reaching it, and every predictor stays in. Only an exact zero drops a predictor from the model.

Complexity is the count of non-zero coefficients

Tracking that count across a grid of regularization strengths gives a second axis alongside performance. Smaller C means a stronger penalty and a sparser model, because C is the inverse of the penalty strength.

The curve rises quickly, then plateaus

Performance climbs as informative predictors enter, then flattens. The pooled validation curve peaks at C = 6.230, and Model A, refit at each fold’s own best value, reached a held-out AUPRC of 0.931. The sparsest model within 5% of that score, at C = 0.024, reached 0.891.

Someone has to pick the margin

Selecting the sparsest model within a set margin of the best validation score makes the trade-off explicit and avoids chasing small differences. The margin is a judgment about acceptable cost, and meeting it does not show the simpler model generalizes better.

Which features are kept can change

When predictors are correlated, small changes in the data change which subset L1 keeps while performance stays similar. Resampling and counting how often each feature survives moves the focus from one fitted model to the pattern that persists.

Nested cross-validation keeps selection out of the score

An inner loop picks C and an outer loop scores the chosen model on data it never saw, repeated over 100 trials of 5 outer folds. The shaded bands are percentile ranges over resampled splits of the same 918 patients, and no classical confidence interval is implied.

Data & License

Dataset

Heart Failure Prediction Dataset by fedesoriano (Kaggle, 2021), combining five UCI heart-disease databases over 11 shared clinical features (918 records). Original clinical data © its creators — Andras Janosi, M.D. (Budapest); William Steinbrunn, M.D. (Zurich); Matthias Pfisterer, M.D. (Basel); and Robert Detrano, M.D., Ph.D. (Long Beach / Cleveland); donor David W. Aha (UCI). Source: kaggle.com/datasets/fedesoriano/heart-failure-prediction.License: Open Database License (ODbL) v1.0 for the database; the contents remain © the original authors named above. Used with attribution. ODbL 1.0

Article

Code

Code examples in this article are licensed under the Common Public Attribution License Version 1.0 (CPAL-1.0), an OSI-approved copyleft license based on the Mozilla Public License 1.1. Initial Developer: Philip Sarajlic.

Modifications must be released in Source Code form under CPAL-1.0. Making the code usable by anyone other than you over a network is External Deployment under the license and is treated as distribution, so the Source Code must be made available to those users.

Full text: opensource.org/license/cpal-1-0 (SPDX identifier CPAL-1.0)

L1 regularization as a mechanism for sparsity

The L1 penalty shrinks coefficients to exactly zero

Defining model complexity

The performance–complexity curve

Diminishing returns and model selection

Performance rises quickly, then plateaus

A practical detail: stability of feature selection

Which features are selected can change from sample to sample

Interpreting results in context

A practical decision framework

Example: Evaluating the trade-off in L1-penalized logistic regression

LASSOSparsity

Overview

Key questions explored

Methodology

Dataset

1. Import dependencies

2. Load data and define constants

Dataset characteristics

Key constants

3. Preprocessing pipeline

Numerical features

Categorical features

Classifier

4. Hyperparameter grid and nested cross-validation setup

Regularization grid

Scoring: AUPRC

Nested cross-validation structure

5. Helper functions for nested cross-validation

average_nonzero_coefficients

evaluate_outer_fold

6. Run nested cross-validation

7. Aggregate inner CV results

Aggregation steps

8. Smooth the performance–complexity curves

Key choices

Interval ordering enforcement

9. Hyperparameter selection plot

AUPRC and complexity across C

10. Evaluate model performance: PRC and ROC curves

How the curves are computed

Why both metrics?

11. Plot PRC and ROC curves

Held-out PRC and ROC by model

Interpreting the PRC and ROC results

Summary and conclusions

Must Read

`average_nonzero_coefficients`

`evaluate_outer_fold`