---
title: Turbo Regression
date: 2025-12-21T12:23:41Z
modified: 2026-07-28T07:44:57Z
permalink: "https://www.micheledpierri.com/2025/12/21/turbo-regression/"
type: post
status: publish
excerpt: ""
wpid: 2399
categories:
  - Data Analysis
  - Statistics
tags:
  - Data Analysis
  - Statistics
  - Data Science
featured_image: "https://www.micheledpierri.com/wp-content/uploads/2025/12/turbo_regression_.png"
featured_image_alt: A vintage pastoral scene of a young boy repairing a mechanical cart in a rural village, surrounded by children, weathered homes, and warm sepia light.
timestamp: 2026-07-28T07:44:57Z
---

## **How to Put a Turbo on Regression Models**

## Bootstrap Validation and Cubic Splines for Medical Data

_A practical guide for data scientists and clinical researchers working with limited medical datasets_

Last updated: December 2025

Author: Michele D. Pierri

Reading time: 15–20 minutes

---

## Introduction — Why regression often disappoints in medicine

[Regression](https://www.micheledpierri.com/wp-content/uploads/wp-mfa-exports/page/linear-regression-in-statistics.md) models are still the backbone of medical research.

Risk scores, prognostic models, outcome prediction after surgery, ICU mortality, disease progression — all of these rely heavily on regression.

Yet, anyone who has tried to deploy a regression model outside the dataset it was built on knows a frustrating truth:

**models that look excellent on paper often perform worse in real patients.**

There are two recurring reasons for this failure.

First, models tend to be **overconfident**. They are evaluated on the same patients they were trained on, which leads to performance estimates that are systematically too optimistic.

Second, models are often **too rigid**. Continuous clinical variables — age, creatinine, hemoglobin, lactate — are forced into linear relationships that do not reflect physiology.

Bootstrap validation and cubic splines in regression address these two issues directly.

Not by replacing regression, but by **making it more honest and more realistic**.

---

## The first problem: optimism in clinical models

Imagine building a logistic regression model to predict 30-day mortality after cardiac surgery.

You use 18 preoperative variables, fit the model, and obtain an AUC of 0.85.

This number feels reassuring.

But where does it come from?

Almost always, it comes from evaluating the model on the **same dataset used to estimate the coefficients**. The model has already “seen” every patient. It has implicitly learned not only the signal, but also the noise.

In clinical terms, this is like testing a diagnostic score on the same cohort that was used to define it.

The result is not wrong — but it is **optimistic**.

The real question clinicians care about is different:

> How will this model behave on the next patient?

---

## Bootstrap validation: simulating future patients

Bootstrap validation answers this question without requiring an external cohort.

Suppose you have a dataset of 500 cardiac surgery patients.

A bootstrap sample is obtained by randomly drawing 500 patients **with replacement** from this dataset.

Some patients will appear multiple times.

Others will not appear at all.

On average:

- about 63% of patients are included at least once (this derives from the probability (1-1/n)^n → e^(-1) ≈ 0.632 as n grows large),
- about 37% are left out.

Those left-out patients are called **out-of-bag (OOB)** observations.

They play a crucial role.

Each bootstrap sample represents a plausible alternative reality in which the same study was conducted, but patient inclusion happened slightly differently.

---

## Where optimism is measured

For each bootstrap sample, we do something very specific:

- We fit the regression model on the bootstrap sample.
- We evaluate its performance on: 
    - the bootstrap sample itself (patients the model has “seen”),
    - the OOB patients (patients the model has _never_ seen).

The difference between these two performances is **optimism**.

In a mortality model, for example, we may observe:

- AUC = 0.87 on bootstrap data,
- AUC = 0.81 on OOB data.

The optimism for that bootstrap iteration is 0.06.

Repeating this process hundreds or thousands of times produces a distribution of optimism estimates.

Their average represents how much our apparent performance is inflated.

---

## Correcting clinical performance estimates

Once optimism is estimated, correction is straightforward.

If the apparent AUC of the original model is 0.85 and the average optimism is 0.05, the optimism-corrected AUC is:

> 0.85 − 0.05 = 0.80

This corrected value is a much better approximation of what we should expect in new patients.

This approach is now standard in high-quality clinical prediction modeling and is strongly recommended over split-sample validation when datasets are limited.

![infopraphic: bootstrap and optimism](https://www.micheledpierri.com/wp-content/uploads/2025/12/Bootstrap-683x1024.png)

---

## The second problem: linearity does not reflect physiology

Even a perfectly validated regression model can still be misleading if its structure is wrong.

Consider age as a predictor of mortality.

Is the risk increase from 40 to 50 the same as from 80 to 90?

Clinically, clearly not.

The same applies to creatinine, lactate, hemoglobin, or blood pressure.

These relationships are **nonlinear**, often with thresholds, plateaus, or acceleration zones.

Forcing them into a straight line is convenient, but biologically implausible.

---

## Cubic splines: letting data bend smoothly

Cubic splines allow regression models to adapt to nonlinear patterns while remaining smooth and interpretable.

Instead of fitting one equation across the entire range of a variable, splines divide the range into intervals separated by **knots**.

Within each interval, a cubic polynomial is fitted, but all pieces are constrained to join smoothly.

This ensures:

- continuity of the curve,
- continuity of the first and second derivatives,
- no abrupt changes in slope.

From a clinical perspective, this means:

- no artificial jumps in risk,
- no arbitrary cutoffs,
- a continuous physiological interpretation.

---

## A clinical example: creatinine and mortality

Creatinine often has a weak association with outcome at low values, but risk increases sharply beyond certain thresholds.

A linear term cannot capture this behavior.

A cubic spline can.

Instead of deciding _a priori_ where risk changes, the spline allows the data to reveal the shape of the relationship — smoothly.

![infographic: understanding cubic spline](https://www.micheledpierri.com/wp-content/uploads/2025/12/Cubic_splines-683x1024.png)

---

## Why splines and bootstrap belong together

Cubic splines increase model flexibility.

Flexibility increases the risk of overfitting.

Bootstrap validation controls that risk.

Together, they form a powerful and principled modeling strategy:

- splines improve biological realism,
- bootstrap ensures honest performance estimation.

This combination is common in high-impact clinical models, even if not always explicitly stated.

---

## Common mistakes in medical applications

The most frequent error is using splines without validation.

Flexible models _must_ be validated more carefully, not less.

Another common mistake is replacing continuous variables with categories.

This throws away information and introduces artificial thresholds.

Splines typically outperform categorization in terms of both accuracy and interpretability.

---

## Limitations and When to Use Alternatives

Bootstrap validation and splines are powerful, but not universal solutions.

**Computational cost**: Bootstrap with 1000 iterations can be slow with large datasets or complex models. If computation is a bottleneck, consider k-fold cross-validation as a faster alternative.

**Very small samples**: With fewer than 100 observations, bootstrap may be unstable. Leave-one-out cross-validation or penalized regression (ridge, lasso) may be more appropriate.

**Rare outcomes**: When the outcome occurs in less than 5-10% of cases, stratified cross-validation ensures adequate representation in each fold. Standard bootstrap may occasionally produce samples with very few events.

**Interpretation requirements**: While spline curves are interpretable, if you need simple risk scores for bedside use (e.g., “add 3 points if age > 65”), linear models with categorization may be more practical despite statistical drawbacks.

The methods described here are most valuable when:

- you have moderate sample sizes (100-5000 observations),
- you need realistic performance estimates,
- biological plausibility matters,
- the model will be used for individual predictions rather than simple screening.

---

# A fully reproducible example: bootstrap and cubic splines on simulated medical data

## Step 1 — Create a synthetic clinical dataset

- **Age**
- **Creatinine**
- **Hemoglobin**
- **Outcome** with _nonlinear_ risk (clinically plausible)

import numpy as np
import pandas as pd

np.random.seed(42)

n = 600

age = np.random.normal(70, 8, n).clip(40, 90)
creatinine = np.random.lognormal(mean=0.2, sigma=0.4, size=n).clip(0.5, 5)
hemoglobin = np.random.normal(13, 1.5, n).clip(8, 18)

# Nonlinear true risk function (unknown to the model)
logit = (
    0.04 * (age - 65)
    + 0.8 * np.maximum(creatinine - 1.2, 0) ** 1.5
    - 0.25 * (hemoglobin - 13)
)

prob = 1 / (1 + np.exp(-logit))
mortality = np.random.binomial(1, prob)

data = pd.DataFrame({
    "age": age,
    "creatinine": creatinine,
    "hemoglobin": hemoglobin,
    "death_30d": mortality
})

data.head()
```
<span class="line"><span style="color: #FF79C6">import</span><span style="color: #F8F8F2"> numpy </span><span style="color: #FF79C6">as</span><span style="color: #F8F8F2"> np</span></span>
<span class="line"><span style="color: #FF79C6">import</span><span style="color: #F8F8F2"> pandas </span><span style="color: #FF79C6">as</span><span style="color: #F8F8F2"> pd</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">np.random.seed(</span><span style="color: #BD93F9">42</span><span style="color: #F8F8F2">)</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">n </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> </span><span style="color: #BD93F9">600</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">age </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> np.random.normal(</span><span style="color: #BD93F9">70</span><span style="color: #F8F8F2">, </span><span style="color: #BD93F9">8</span><span style="color: #F8F8F2">, n).clip(</span><span style="color: #BD93F9">40</span><span style="color: #F8F8F2">, </span><span style="color: #BD93F9">90</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #F8F8F2">creatinine </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> np.random.lognormal(</span><span style="color: #FFB86C; font-style: italic">mean</span><span style="color: #FF79C6">=</span><span style="color: #BD93F9">0.2</span><span style="color: #F8F8F2">, </span><span style="color: #FFB86C; font-style: italic">sigma</span><span style="color: #FF79C6">=</span><span style="color: #BD93F9">0.4</span><span style="color: #F8F8F2">, </span><span style="color: #FFB86C; font-style: italic">size</span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2">n).clip(</span><span style="color: #BD93F9">0.5</span><span style="color: #F8F8F2">, </span><span style="color: #BD93F9">5</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #F8F8F2">hemoglobin </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> np.random.normal(</span><span style="color: #BD93F9">13</span><span style="color: #F8F8F2">, </span><span style="color: #BD93F9">1.5</span><span style="color: #F8F8F2">, n).clip(</span><span style="color: #BD93F9">8</span><span style="color: #F8F8F2">, </span><span style="color: #BD93F9">18</span><span style="color: #F8F8F2">)</span></span>
<span class="line"></span>
<span class="line"><span style="color: #6272A4"># Nonlinear true risk function (unknown to the model)</span></span>
<span class="line"><span style="color: #F8F8F2">logit </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> (</span></span>
<span class="line"><span style="color: #F8F8F2">    </span><span style="color: #BD93F9">0.04</span><span style="color: #F8F8F2"> </span><span style="color: #FF79C6">*</span><span style="color: #F8F8F2"> (age </span><span style="color: #FF79C6">-</span><span style="color: #F8F8F2"> </span><span style="color: #BD93F9">65</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #F8F8F2">    </span><span style="color: #FF79C6">+</span><span style="color: #F8F8F2"> </span><span style="color: #BD93F9">0.8</span><span style="color: #F8F8F2"> </span><span style="color: #FF79C6">*</span><span style="color: #F8F8F2"> np.maximum(creatinine </span><span style="color: #FF79C6">-</span><span style="color: #F8F8F2"> </span><span style="color: #BD93F9">1.2</span><span style="color: #F8F8F2">, </span><span style="color: #BD93F9">0</span><span style="color: #F8F8F2">) </span><span style="color: #FF79C6">**</span><span style="color: #F8F8F2"> </span><span style="color: #BD93F9">1.5</span></span>
<span class="line"><span style="color: #F8F8F2">    </span><span style="color: #FF79C6">-</span><span style="color: #F8F8F2"> </span><span style="color: #BD93F9">0.25</span><span style="color: #F8F8F2"> </span><span style="color: #FF79C6">*</span><span style="color: #F8F8F2"> (hemoglobin </span><span style="color: #FF79C6">-</span><span style="color: #F8F8F2"> </span><span style="color: #BD93F9">13</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #F8F8F2">)</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">prob </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> </span><span style="color: #BD93F9">1</span><span style="color: #F8F8F2"> </span><span style="color: #FF79C6">/</span><span style="color: #F8F8F2"> (</span><span style="color: #BD93F9">1</span><span style="color: #F8F8F2"> </span><span style="color: #FF79C6">+</span><span style="color: #F8F8F2"> np.exp(</span><span style="color: #FF79C6">-</span><span style="color: #F8F8F2">logit))</span></span>
<span class="line"><span style="color: #F8F8F2">mortality </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> np.random.binomial(</span><span style="color: #BD93F9">1</span><span style="color: #F8F8F2">, prob)</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">data </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> pd.DataFrame({</span></span>
<span class="line"><span style="color: #F8F8F2">    </span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">age</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">: age,</span></span>
<span class="line"><span style="color: #F8F8F2">    </span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">creatinine</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">: creatinine,</span></span>
<span class="line"><span style="color: #F8F8F2">    </span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">hemoglobin</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">: hemoglobin,</span></span>
<span class="line"><span style="color: #F8F8F2">    </span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">death_30d</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">: mortality</span></span>
<span class="line"><span style="color: #F8F8F2">})</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">data.head()</span></span>
<span class="line"></span>
```

👉 This dataset is **nonlinear by design**, just like clinical reality.

---

## Step 2 — Apparent performance of a simple logistic regression

Let’s start with a **classic linear model**.

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score

X = data[["age", "creatinine", "hemoglobin"]]
y = data["death_30d"]

model = LogisticRegression(max_iter=1000)
model.fit(X, y)

apparent_auc = roc_auc_score(y, model.predict_proba(X)[:, 1])
print(f"Apparent AUC: {apparent_auc:.3f}")
```
<span class="line"><span style="color: #FF79C6">from</span><span style="color: #F8F8F2"> sklearn.linear_model </span><span style="color: #FF79C6">import</span><span style="color: #F8F8F2"> LogisticRegression</span></span>
<span class="line"><span style="color: #FF79C6">from</span><span style="color: #F8F8F2"> sklearn.metrics </span><span style="color: #FF79C6">import</span><span style="color: #F8F8F2"> roc_auc_score</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">X </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> data[[</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">age</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">, </span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">creatinine</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">, </span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">hemoglobin</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">]]</span></span>
<span class="line"><span style="color: #F8F8F2">y </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> data[</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">death_30d</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">]</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">model </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> LogisticRegression(</span><span style="color: #FFB86C; font-style: italic">max_iter</span><span style="color: #FF79C6">=</span><span style="color: #BD93F9">1000</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #F8F8F2">model.fit(X, y)</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">apparent_auc </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> roc_auc_score(y, model.predict_proba(X)[:, </span><span style="color: #BD93F9">1</span><span style="color: #F8F8F2">])</span></span>
<span class="line"><span style="color: #8BE9FD">print</span><span style="color: #F8F8F2">(</span><span style="color: #FF79C6">f</span><span style="color: #F1FA8C">"Apparent AUC: </span><span style="color: #BD93F9">{</span><span style="color: #F8F8F2">apparent_auc</span><span style="color: #FF79C6">:.3f</span><span style="color: #BD93F9">}</span><span style="color: #F1FA8C">"</span><span style="color: #F8F8F2">)</span></span>
<span class="line"></span>
```

Result: Apparent AUC: 0.661![🔍](blob:https://www.micheledpierri.com/c3ddffd8-5f2e-4e6f-a95e-702d9d121bdb)

This AUC is optimistic: the model is evaluated on the same patients.

---

## Step 3 — Bootstrap optimism correction (fully reproducible)

Now let’s estimate the optimism.

from sklearn.utils import resample

n_boot = 500
optimism = []

for i in range(n_boot):
    boot_idx = resample(np.arange(len(data)), replace=True)
    oob_idx = np.setdiff1d(np.arange(len(data)), boot_idx)

    if len(oob_idx) < 30:
        continue

    X_boot, y_boot = X.iloc[boot_idx], y.iloc[boot_idx]
    X_oob, y_oob = X.iloc[oob_idx], y.iloc[oob_idx]

    model.fit(X_boot, y_boot)

    auc_boot = roc_auc_score(y_boot, model.predict_proba(X_boot)[:, 1])
    auc_oob = roc_auc_score(y_oob, model.predict_proba(X_oob)[:, 1])

    optimism.append(auc_boot - auc_oob)

mean_optimism = np.mean(optimism)
corrected_auc = apparent_auc - mean_optimism

print(f"Mean optimism: {mean_optimism:.3f}")
print(f"Optimism-corrected AUC: {corrected_auc:.3f}")
```
<span class="line"><span style="color: #FF79C6">from</span><span style="color: #F8F8F2"> sklearn.utils </span><span style="color: #FF79C6">import</span><span style="color: #F8F8F2"> resample</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">n_boot </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> </span><span style="color: #BD93F9">500</span></span>
<span class="line"><span style="color: #F8F8F2">optimism </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> []</span></span>
<span class="line"></span>
<span class="line"><span style="color: #FF79C6">for</span><span style="color: #F8F8F2"> i </span><span style="color: #FF79C6">in</span><span style="color: #F8F8F2"> </span><span style="color: #8BE9FD">range</span><span style="color: #F8F8F2">(n_boot):</span></span>
<span class="line"><span style="color: #F8F8F2">    boot_idx </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> resample(np.arange(</span><span style="color: #8BE9FD">len</span><span style="color: #F8F8F2">(data)), </span><span style="color: #FFB86C; font-style: italic">replace</span><span style="color: #FF79C6">=</span><span style="color: #BD93F9">True</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #F8F8F2">    oob_idx </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> np.setdiff1d(np.arange(</span><span style="color: #8BE9FD">len</span><span style="color: #F8F8F2">(data)), boot_idx)</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">    </span><span style="color: #FF79C6">if</span><span style="color: #F8F8F2"> </span><span style="color: #8BE9FD">len</span><span style="color: #F8F8F2">(oob_idx) </span><span style="color: #FF79C6"><</span><span style="color: #F8F8F2"> </span><span style="color: #BD93F9">30</span><span style="color: #F8F8F2">:</span></span>
<span class="line"><span style="color: #F8F8F2">        </span><span style="color: #FF79C6">continue</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">    X_boot, y_boot </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> X.iloc[boot_idx], y.iloc[boot_idx]</span></span>
<span class="line"><span style="color: #F8F8F2">    X_oob, y_oob </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> X.iloc[oob_idx], y.iloc[oob_idx]</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">    model.fit(X_boot, y_boot)</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">    auc_boot </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> roc_auc_score(y_boot, model.predict_proba(X_boot)[:, </span><span style="color: #BD93F9">1</span><span style="color: #F8F8F2">])</span></span>
<span class="line"><span style="color: #F8F8F2">    auc_oob </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> roc_auc_score(y_oob, model.predict_proba(X_oob)[:, </span><span style="color: #BD93F9">1</span><span style="color: #F8F8F2">])</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">    optimism.append(auc_boot </span><span style="color: #FF79C6">-</span><span style="color: #F8F8F2"> auc_oob)</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">mean_optimism </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> np.mean(optimism)</span></span>
<span class="line"><span style="color: #F8F8F2">corrected_auc </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> apparent_auc </span><span style="color: #FF79C6">-</span><span style="color: #F8F8F2"> mean_optimism</span></span>
<span class="line"></span>
<span class="line"><span style="color: #8BE9FD">print</span><span style="color: #F8F8F2">(</span><span style="color: #FF79C6">f</span><span style="color: #F1FA8C">"Mean optimism: </span><span style="color: #BD93F9">{</span><span style="color: #F8F8F2">mean_optimism</span><span style="color: #FF79C6">:.3f</span><span style="color: #BD93F9">}</span><span style="color: #F1FA8C">"</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #8BE9FD">print</span><span style="color: #F8F8F2">(</span><span style="color: #FF79C6">f</span><span style="color: #F1FA8C">"Optimism-corrected AUC: </span><span style="color: #BD93F9">{</span><span style="color: #F8F8F2">corrected_auc</span><span style="color: #FF79C6">:.3f</span><span style="color: #BD93F9">}</span><span style="color: #F1FA8C">"</span><span style="color: #F8F8F2">)</span></span>
<span class="line"></span>
```

Result: Mean optimism: 0.014 Optimism-corrected AUC: 0.647

---

## Step 4 — Why linear regression is inadequate here

Now let’s look at the true (simulated) relationship between creatinine and risk.

import matplotlib.pyplot as plt

plt.scatter(data["creatinine"], prob, alpha=0.3)
plt.xlabel("Creatinine (mg/dL)")
plt.ylabel("True mortality risk")
plt.title("True nonlinear relationship (unknown to the model)")
plt.show()
```
<span class="line"><span style="color: #FF79C6">import</span><span style="color: #F8F8F2"> matplotlib.pyplot </span><span style="color: #FF79C6">as</span><span style="color: #F8F8F2"> plt</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">plt.scatter(data[</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">creatinine</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">], prob, </span><span style="color: #FFB86C; font-style: italic">alpha</span><span style="color: #FF79C6">=</span><span style="color: #BD93F9">0.3</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #F8F8F2">plt.xlabel(</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">Creatinine (mg/dL)</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #F8F8F2">plt.ylabel(</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">True mortality risk</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #F8F8F2">plt.title(</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">True nonlinear relationship (unknown to the model)</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #F8F8F2">plt.show()</span></span>
<span class="line"></span>
```

![👉](blob:https://www.micheledpierri.com/984d40d3-694b-40b0-885e-8b6a347df663) No clinician would expect a linear relationship.

Yet the previous model forces it.

---

## Step 5 — Logistic regression with cubic splines

Now we introduce cubic splines.

import statsmodels.api as sm
from patsy import dmatrix

spline_creatinine = dmatrix(
    "bs(creatinine, df=4, include_intercept=False)",
    data,
    return_type="dataframe"
)

X_spline = pd.concat([
    data[["age", "hemoglobin"]],
    spline_creatinine
], axis=1)

X_spline = sm.add_constant(X_spline)

model_spline = sm.Logit(y, X_spline).fit(disp=False)

pred_spline = model_spline.predict(X_spline)
spline_auc = roc_auc_score(y, pred_spline)

print(f"Spline model apparent AUC: {spline_auc:.3f}")
```
<span class="line"><span style="color: #FF79C6">import</span><span style="color: #F8F8F2"> statsmodels.api </span><span style="color: #FF79C6">as</span><span style="color: #F8F8F2"> sm</span></span>
<span class="line"><span style="color: #FF79C6">from</span><span style="color: #F8F8F2"> patsy </span><span style="color: #FF79C6">import</span><span style="color: #F8F8F2"> dmatrix</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">spline_creatinine </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> dmatrix(</span></span>
<span class="line"><span style="color: #F8F8F2">    </span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">bs(creatinine, df=4, include_intercept=False)</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">,</span></span>
<span class="line"><span style="color: #F8F8F2">    data,</span></span>
<span class="line"><span style="color: #F8F8F2">    </span><span style="color: #FFB86C; font-style: italic">return_type</span><span style="color: #FF79C6">=</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">dataframe</span><span style="color: #E9F284">"</span></span>
<span class="line"><span style="color: #F8F8F2">)</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">X_spline </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> pd.concat([</span></span>
<span class="line"><span style="color: #F8F8F2">    data[[</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">age</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">, </span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">hemoglobin</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">]],</span></span>
<span class="line"><span style="color: #F8F8F2">    spline_creatinine</span></span>
<span class="line"><span style="color: #F8F8F2">], </span><span style="color: #FFB86C; font-style: italic">axis</span><span style="color: #FF79C6">=</span><span style="color: #BD93F9">1</span><span style="color: #F8F8F2">)</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">X_spline </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> sm.add_constant(X_spline)</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">model_spline </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> sm.Logit(y, X_spline).fit(</span><span style="color: #FFB86C; font-style: italic">disp</span><span style="color: #FF79C6">=</span><span style="color: #BD93F9">False</span><span style="color: #F8F8F2">)</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">pred_spline </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> model_spline.predict(X_spline)</span></span>
<span class="line"><span style="color: #F8F8F2">spline_auc </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> roc_auc_score(y, pred_spline)</span></span>
<span class="line"></span>
<span class="line"><span style="color: #8BE9FD">print</span><span style="color: #F8F8F2">(</span><span style="color: #FF79C6">f</span><span style="color: #F1FA8C">"Spline model apparent AUC: </span><span style="color: #BD93F9">{</span><span style="color: #F8F8F2">spline_auc</span><span style="color: #FF79C6">:.3f</span><span style="color: #BD93F9">}</span><span style="color: #F1FA8C">"</span><span style="color: #F8F8F2">)</span></span>
<span class="line"></span>
```

Result: Spline model apparent AUC: 0.669

The model now **adapts to physiology**, not the other way around.

---

## Step 6 — Bootstrap validation of the spline model

Flexibility must be paid for with rigorous validation.

optimism_spline = []

for i in range(n_boot):
    boot_idx = resample(np.arange(len(data)), replace=True)
    oob_idx = np.setdiff1d(np.arange(len(data)), boot_idx)

    if len(oob_idx) < 30:
        continue

    data_boot = data.iloc[boot_idx]
    data_oob = data.iloc[oob_idx]

    Xb = sm.add_constant(pd.concat([
        data_boot[["age", "hemoglobin"]],
        dmatrix("bs(creatinine, df=4, include_intercept=False)",
                data_boot, return_type="dataframe")
    ], axis=1))

    Xo = sm.add_constant(pd.concat([
        data_oob[["age", "hemoglobin"]],
        dmatrix("bs(creatinine, df=4, include_intercept=False)",
                data_oob, return_type="dataframe")
    ], axis=1))

    yb = data_boot["death_30d"]
    yo = data_oob["death_30d"]

    try:
        m = sm.Logit(yb, Xb).fit(disp=False)
        auc_b = roc_auc_score(yb, m.predict(Xb))
        auc_o = roc_auc_score(yo, m.predict(Xo))
        optimism_spline.append(auc_b - auc_o)
    except:
        continue

mean_opt_spline = np.mean(optimism_spline)
corrected_auc_spline = spline_auc - mean_opt_spline

print(f"Spline optimism: {mean_opt_spline:.3f}")
print(f"Spline corrected AUC: {corrected_auc_spline:.3f}")
```
<span class="line"><span style="color: #F8F8F2">optimism_spline </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> []</span></span>
<span class="line"></span>
<span class="line"><span style="color: #FF79C6">for</span><span style="color: #F8F8F2"> i </span><span style="color: #FF79C6">in</span><span style="color: #F8F8F2"> </span><span style="color: #8BE9FD">range</span><span style="color: #F8F8F2">(n_boot):</span></span>
<span class="line"><span style="color: #F8F8F2">    boot_idx </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> resample(np.arange(</span><span style="color: #8BE9FD">len</span><span style="color: #F8F8F2">(data)), </span><span style="color: #FFB86C; font-style: italic">replace</span><span style="color: #FF79C6">=</span><span style="color: #BD93F9">True</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #F8F8F2">    oob_idx </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> np.setdiff1d(np.arange(</span><span style="color: #8BE9FD">len</span><span style="color: #F8F8F2">(data)), boot_idx)</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">    </span><span style="color: #FF79C6">if</span><span style="color: #F8F8F2"> </span><span style="color: #8BE9FD">len</span><span style="color: #F8F8F2">(oob_idx) </span><span style="color: #FF79C6"><</span><span style="color: #F8F8F2"> </span><span style="color: #BD93F9">30</span><span style="color: #F8F8F2">:</span></span>
<span class="line"><span style="color: #F8F8F2">        </span><span style="color: #FF79C6">continue</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">    data_boot </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> data.iloc[boot_idx]</span></span>
<span class="line"><span style="color: #F8F8F2">    data_oob </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> data.iloc[oob_idx]</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">    Xb </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> sm.add_constant(pd.concat([</span></span>
<span class="line"><span style="color: #F8F8F2">        data_boot[[</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">age</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">, </span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">hemoglobin</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">]],</span></span>
<span class="line"><span style="color: #F8F8F2">        dmatrix(</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">bs(creatinine, df=4, include_intercept=False)</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">,</span></span>
<span class="line"><span style="color: #F8F8F2">                data_boot, </span><span style="color: #FFB86C; font-style: italic">return_type</span><span style="color: #FF79C6">=</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">dataframe</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #F8F8F2">    ], </span><span style="color: #FFB86C; font-style: italic">axis</span><span style="color: #FF79C6">=</span><span style="color: #BD93F9">1</span><span style="color: #F8F8F2">))</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">    Xo </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> sm.add_constant(pd.concat([</span></span>
<span class="line"><span style="color: #F8F8F2">        data_oob[[</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">age</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">, </span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">hemoglobin</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">]],</span></span>
<span class="line"><span style="color: #F8F8F2">        dmatrix(</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">bs(creatinine, df=4, include_intercept=False)</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">,</span></span>
<span class="line"><span style="color: #F8F8F2">                data_oob, </span><span style="color: #FFB86C; font-style: italic">return_type</span><span style="color: #FF79C6">=</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">dataframe</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #F8F8F2">    ], </span><span style="color: #FFB86C; font-style: italic">axis</span><span style="color: #FF79C6">=</span><span style="color: #BD93F9">1</span><span style="color: #F8F8F2">))</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">    yb </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> data_boot[</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">death_30d</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">]</span></span>
<span class="line"><span style="color: #F8F8F2">    yo </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> data_oob[</span><span style="color: #E9F284">"</span><span style="color: #F1FA8C">death_30d</span><span style="color: #E9F284">"</span><span style="color: #F8F8F2">]</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">    </span><span style="color: #FF79C6">try</span><span style="color: #F8F8F2">:</span></span>
<span class="line"><span style="color: #F8F8F2">        m </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> sm.Logit(yb, Xb).fit(</span><span style="color: #FFB86C; font-style: italic">disp</span><span style="color: #FF79C6">=</span><span style="color: #BD93F9">False</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #F8F8F2">        auc_b </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> roc_auc_score(yb, m.predict(Xb))</span></span>
<span class="line"><span style="color: #F8F8F2">        auc_o </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> roc_auc_score(yo, m.predict(Xo))</span></span>
<span class="line"><span style="color: #F8F8F2">        optimism_spline.append(auc_b </span><span style="color: #FF79C6">-</span><span style="color: #F8F8F2"> auc_o)</span></span>
<span class="line"><span style="color: #F8F8F2">    </span><span style="color: #FF79C6">except</span><span style="color: #F8F8F2">:</span></span>
<span class="line"><span style="color: #F8F8F2">        </span><span style="color: #FF79C6">continue</span></span>
<span class="line"></span>
<span class="line"><span style="color: #F8F8F2">mean_opt_spline </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> np.mean(optimism_spline)</span></span>
<span class="line"><span style="color: #F8F8F2">corrected_auc_spline </span><span style="color: #FF79C6">=</span><span style="color: #F8F8F2"> spline_auc </span><span style="color: #FF79C6">-</span><span style="color: #F8F8F2"> mean_opt_spline</span></span>
<span class="line"></span>
<span class="line"><span style="color: #8BE9FD">print</span><span style="color: #F8F8F2">(</span><span style="color: #FF79C6">f</span><span style="color: #F1FA8C">"Spline optimism: </span><span style="color: #BD93F9">{</span><span style="color: #F8F8F2">mean_opt_spline</span><span style="color: #FF79C6">:.3f</span><span style="color: #BD93F9">}</span><span style="color: #F1FA8C">"</span><span style="color: #F8F8F2">)</span></span>
<span class="line"><span style="color: #8BE9FD">print</span><span style="color: #F8F8F2">(</span><span style="color: #FF79C6">f</span><span style="color: #F1FA8C">"Spline corrected AUC: </span><span style="color: #BD93F9">{</span><span style="color: #F8F8F2">corrected_auc_spline</span><span style="color: #FF79C6">:.3f</span><span style="color: #BD93F9">}</span><span style="color: #F1FA8C">"</span><span style="color: #F8F8F2">)</span></span>
<span class="line"></span>
```

Result: Spline optimism: 0.030 Spline corrected AUC: 0.639

---

## Conclusion — A modern view of regression in medicine

Regression is not outdated.

It is underutilized.

Bootstrap validation teaches regression humility.

Cubic splines give it flexibility.

Together, they transform regression from a blunt instrument into a refined clinical modeling tool.

If regression is the engine of medical prediction,

**bootstrap and splines are what finally put the turbo on it.**

---

## References and Further Reading

[**Harrell FE.** _Regression Modeling Strategies: With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis._ 2nd ed. Springer, 2015.](https://nibmehub.com/opac-service/pdf/read/Regression%20Modeling%20Strategies-%202nd%20edition-%202015.pdf)

> The definitive reference on modern regression techniques in biostatistics, including extensive coverage of splines and bootstrap validation.

[**Steyerberg EW.** _Clinical Prediction Models: A Practical Approach to Development, Validation, and Updating._ 2nd ed. Springer, 2019](https://www.clinicalpredictionmodels.org/).

> Comprehensive guide to building and validating prediction models in medicine, with emphasis on optimism correction.

[**Efron B, Tibshirani RJ.** _An Introduction to the Bootstrap._ Chapman & Hall, 1993.](https://www.taylorfrancis.com/books/mono/10.1201/9780429246593/introduction-bootstrap-bradley-efron-tibshirani)

> The foundational text on bootstrap methods by their inventor.

[**Collins GS, Reitsma JB, Altman DG, Moons KGM.** Transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD): the TRIPOD Statement. _BMJ_ 2015;350:g7594.](https://pubmed.ncbi.nlm.nih.gov/25569120/)

> Essential guidelines for reporting clinical prediction models, emphasizing proper validation.