"""
Supplementary File 1
Reproducible analysis for the manuscript:
"Modelling Mortality Risk in Malaria Patients Using Logistic Regression Model:
A Case Study of Nakuru Level 6 Hospital"

This script reproduces the reported analysis from the de-identified dataset:
  1. Descriptive checks and variable coding
  2. Variance inflation factors (VIFs)
  3. Box-Tidwell functional-form assessment
  4. Original five-predictor linear logistic regression model
  5. Corrected three-knot restricted cubic spline (RCS) model
  6. Linear-vs-RCS likelihood-ratio comparison and nonlinear component tests
  7. Hosmer-Lemeshow grouped goodness-of-fit statistic
  8. Formal Riley-style sample-size/shrinkage calculation
  9. Clinically interpretable RCS odds-ratio contrasts with 95% CIs
 10. Bootstrap optimism correction (1,000 resamples)
 11. Repeated 20 x 10-fold stratified CV that repeats functional-form assessment,
     knot estimation, model fitting, and Youden threshold selection inside training data
 12. Ridge-penalized RCS sensitivity analysis with nested cross-validation
 13. Threshold-specific confusion matrices and Wilson 95% CIs
 14. Smooth calibration, ROC-comparison, RCS effect, and decision-curve figures

Expected input file (same directory as this script unless DATA_PATH is edited):
    nakuru_malaria_dataset_2020_2026.csv

Outputs are written to:
    supplementary_outputs/

The random seed is fixed for reproducibility.
"""

from pathlib import Path
import sys
import warnings

import numpy as np
import pandas as pd
import scipy
from scipy.stats import chi2
import statsmodels
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.proportion import proportion_confint
from statsmodels.nonparametric.smoothers_lowess import lowess
import sklearn
from sklearn.model_selection import (
    RepeatedStratifiedKFold,
    StratifiedKFold,
    GridSearchCV,
)
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (
    roc_auc_score,
    brier_score_loss,
    roc_curve,
    confusion_matrix,
)
import matplotlib
import matplotlib.pyplot as plt

warnings.filterwarnings("ignore")

# -----------------------------------------------------------------------------
# Reproducibility settings
# -----------------------------------------------------------------------------
SEED = 20260818
BOOTSTRAP_REPS = 1000
CV_REPEATS = 20
CV_FOLDS = 10
BT_ALPHA = 0.05
C_GRID = np.logspace(-3, 3, 25)

HERE = Path(__file__).resolve().parent
DATA_PATH = HERE / "nakuru_malaria_dataset_2020_2026.csv"
OUT = HERE / "supplementary_outputs"
OUT.mkdir(exist_ok=True)

# -----------------------------------------------------------------------------
# Utility functions
# -----------------------------------------------------------------------------
def expit(x):
    x = np.asarray(x, dtype=float)
    return 1.0 / (1.0 + np.exp(-x))


def rcs_basis(x, knots):
    """Three-knot restricted cubic spline basis: linear + one nonlinear term."""
    x = np.asarray(x, dtype=float)
    k1, k2, k3 = np.asarray(knots, dtype=float)
    nonlinear = (
        np.maximum(x - k1, 0) ** 3
        - np.maximum(x - k2, 0) ** 3 * (k3 - k1) / (k3 - k2)
        + np.maximum(x - k3, 0) ** 3 * (k2 - k1) / (k3 - k2)
    ) / ((k3 - k1) ** 2)
    return np.column_stack([x, nonlinear])


def knot_rule(d):
    """10th, 50th and 90th percentiles, calculated in the supplied data only."""
    return {
        "Age": np.quantile(d["Age"], [0.10, 0.50, 0.90]),
        "Hb": np.quantile(d["Haemoglobin_g_dL"], [0.10, 0.50, 0.90]),
        "Platelet": np.quantile(d["Platelet_Count_10^9_L"], [0.10, 0.50, 0.90]),
        "Parasite": np.quantile(d["parasite_1000"], [0.10, 0.50, 0.90]),
    }


def linear_design(d):
    return pd.DataFrame(
        {
            "Age": d["Age"].to_numpy(),
            "Male": d["male"].to_numpy(),
            "Parasite_per_1000": d["parasite_1000"].to_numpy(),
            "Hb": d["Haemoglobin_g_dL"].to_numpy(),
            "Platelet": d["Platelet_Count_10^9_L"].to_numpy(),
        },
        index=d.index,
    )


def fixed_rcs_design(d, knots):
    """Corrected manuscript model: Age/Hb/Platelet spline; parasite linear; sex binary."""
    a = rcs_basis(d["Age"], knots["Age"])
    h = rcs_basis(d["Haemoglobin_g_dL"], knots["Hb"])
    p = rcs_basis(d["Platelet_Count_10^9_L"], knots["Platelet"])
    return pd.DataFrame(
        {
            "Age_linear": a[:, 0],
            "Age_nonlinear": a[:, 1],
            "Male": d["male"].to_numpy(),
            "Parasite_per_1000": d["parasite_1000"].to_numpy(),
            "Hb_linear": h[:, 0],
            "Hb_nonlinear": h[:, 1],
            "Platelet_linear": p[:, 0],
            "Platelet_nonlinear": p[:, 1],
        },
        index=d.index,
    )


def box_tidwell(d):
    """
    Joint Box-Tidwell assessment for continuous predictors.

    For age, Age*ln(Age+1) is used because the dataset includes age=0.
    All other continuous predictors are strictly positive and use x*ln(x).
    The test is fitted jointly with the original five predictors.
    """
    age = d["Age"].astype(float)
    parasite = d["parasite_1000"].astype(float)
    hb = d["Haemoglobin_g_dL"].astype(float)
    platelet = d["Platelet_Count_10^9_L"].astype(float)

    if (parasite <= 0).any() or (hb <= 0).any() or (platelet <= 0).any():
        raise ValueError("Box-Tidwell requires positive parasite, Hb and platelet values.")

    X = pd.DataFrame(
        {
            "Age": age,
            "Male": d["male"].to_numpy(),
            "Parasite_per_1000": parasite,
            "Hb": hb,
            "Platelet": platelet,
            "BT_Age": age * np.log(age + 1.0),
            "BT_Parasite": parasite * np.log(parasite),
            "BT_Hb": hb * np.log(hb),
            "BT_Platelet": platelet * np.log(platelet),
        },
        index=d.index,
    )
    fit = sm.Logit(d["y"], sm.add_constant(X, has_constant="add")).fit(
        disp=False, method="newton", maxiter=500
    )
    result = pd.DataFrame(
        {
            "Predictor": ["Age", "Parasite density", "Haemoglobin", "Platelet count"],
            "Interaction_term": ["BT_Age", "BT_Parasite", "BT_Hb", "BT_Platelet"],
            "Coefficient": [
                fit.params["BT_Age"],
                fit.params["BT_Parasite"],
                fit.params["BT_Hb"],
                fit.params["BT_Platelet"],
            ],
            "SE": [
                fit.bse["BT_Age"],
                fit.bse["BT_Parasite"],
                fit.bse["BT_Hb"],
                fit.bse["BT_Platelet"],
            ],
            "p_value": [
                fit.pvalues["BT_Age"],
                fit.pvalues["BT_Parasite"],
                fit.pvalues["BT_Hb"],
                fit.pvalues["BT_Platelet"],
            ],
        }
    )
    return result, fit


def vif_table(d):
    X = sm.add_constant(linear_design(d), has_constant="add")
    out = []
    for i, col in enumerate(X.columns):
        if col == "const":
            continue
        out.append((col, variance_inflation_factor(X.to_numpy(), i)))
    return pd.DataFrame(out, columns=["Predictor", "VIF"])


def calibration_stats(y, p):
    """Calibration-in-the-large (intercept with slope fixed at 1) and slope."""
    y = np.asarray(y, dtype=int)
    p = np.clip(np.asarray(p, dtype=float), 1e-8, 1 - 1e-8)
    lp = np.log(p / (1 - p))
    slope_model = sm.Logit(y, sm.add_constant(lp)).fit(disp=False, maxiter=200)
    cil_model = sm.GLM(
        y,
        np.ones((len(y), 1)),
        family=sm.families.Binomial(),
        offset=lp,
    ).fit()
    return float(cil_model.params[0]), float(slope_model.params[1])


def hosmer_lemeshow(y, p, groups=10):
    """Hosmer-Lemeshow statistic based on quantile groups of fitted risk."""
    tmp = pd.DataFrame({"y": np.asarray(y), "p": np.asarray(p)})
    tmp["group"] = pd.qcut(tmp["p"], q=groups, labels=False, duplicates="drop")
    tab = tmp.groupby("group", observed=True).agg(
        n=("y", "size"), observed_deaths=("y", "sum"), expected_deaths=("p", "sum")
    )
    e = tab["expected_deaths"].to_numpy()
    n = tab["n"].to_numpy()
    o = tab["observed_deaths"].to_numpy()
    stat = np.sum((o - e) ** 2 / (e * (1.0 - e / n)))
    df_hl = len(tab) - 2
    p_value = chi2.sf(stat, df_hl)
    tab["observed_survivals"] = tab["n"] - tab["observed_deaths"]
    tab["expected_survivals"] = tab["n"] - tab["expected_deaths"]
    return float(stat), int(df_hl), float(p_value), tab.reset_index()


def youden_threshold(y, p):
    fpr, tpr, thresholds = roc_curve(y, p)
    return float(thresholds[np.argmax(tpr - fpr)])


def classification_metrics(y, p, threshold):
    y = np.asarray(y, dtype=int)
    p = np.asarray(p, dtype=float)
    pred = (p >= threshold).astype(int)
    tn, fp, fn, tp = confusion_matrix(y, pred, labels=[0, 1]).ravel()

    sensitivity = tp / (tp + fn)
    specificity = tn / (tn + fp)
    ppv = tp / (tp + fp) if (tp + fp) else np.nan
    npv = tn / (tn + fn) if (tn + fn) else np.nan
    accuracy = (tp + tn) / len(y)

    return {
        "Threshold": threshold,
        "TN": int(tn),
        "FP": int(fp),
        "FN": int(fn),
        "TP": int(tp),
        "Sensitivity": sensitivity,
        "Sensitivity_low": proportion_confint(tp, tp + fn, method="wilson")[0],
        "Sensitivity_high": proportion_confint(tp, tp + fn, method="wilson")[1],
        "Specificity": specificity,
        "Specificity_low": proportion_confint(tn, tn + fp, method="wilson")[0],
        "Specificity_high": proportion_confint(tn, tn + fp, method="wilson")[1],
        "PPV": ppv,
        "PPV_low": proportion_confint(tp, tp + fp, method="wilson")[0],
        "PPV_high": proportion_confint(tp, tp + fp, method="wilson")[1],
        "NPV": npv,
        "NPV_low": proportion_confint(tn, tn + fn, method="wilson")[0],
        "NPV_high": proportion_confint(tn, tn + fn, method="wilson")[1],
        "Accuracy": accuracy,
    }


def riley_shrinkage_sample_size(predictor_parameters, prevalence, anticipated_cs_r2, target_shrinkage=0.90):
    """
    Riley et al. shrinkage criterion for a binary-outcome prediction model.
    n = p / ((S - 1) * ln(1 - R2_CS/S))
    """
    p = float(predictor_parameters)
    S = float(target_shrinkage)
    R2 = float(anticipated_cs_r2)
    n = p / ((S - 1.0) * np.log(1.0 - R2 / S))
    events = n * prevalence
    return float(n), float(events)


def spline_or_contrast(model, knots, variable, value, reference):
    """Adjusted OR contrast and Wald 95% CI for one spline-modeled predictor."""
    names = list(model.params.index)
    d = np.zeros(len(names), dtype=float)

    if variable == "Age":
        bx = rcs_basis([value], knots["Age"])[0]
        br = rcs_basis([reference], knots["Age"])[0]
        d[names.index("Age_linear")] = bx[0] - br[0]
        d[names.index("Age_nonlinear")] = bx[1] - br[1]
    elif variable == "Haemoglobin":
        bx = rcs_basis([value], knots["Hb"])[0]
        br = rcs_basis([reference], knots["Hb"])[0]
        d[names.index("Hb_linear")] = bx[0] - br[0]
        d[names.index("Hb_nonlinear")] = bx[1] - br[1]
    elif variable == "Platelet":
        bx = rcs_basis([value], knots["Platelet"])[0]
        br = rcs_basis([reference], knots["Platelet"])[0]
        d[names.index("Platelet_linear")] = bx[0] - br[0]
        d[names.index("Platelet_nonlinear")] = bx[1] - br[1]
    else:
        raise ValueError("variable must be Age, Haemoglobin, or Platelet")

    beta = model.params.to_numpy()
    cov = model.cov_params().to_numpy()
    log_or = float(d @ beta)
    se = float(np.sqrt(d @ cov @ d))
    return {
        "Variable": variable,
        "Value": value,
        "Reference": reference,
        "OR": np.exp(log_or),
        "CI_low": np.exp(log_or - 1.96 * se),
        "CI_high": np.exp(log_or + 1.96 * se),
        "log_OR": log_or,
        "SE_log_OR": se,
    }


def nonlinear_component_test(df, X_rcs, full_model, column):
    reduced_X = X_rcs.drop(columns=[column])
    reduced = sm.Logit(df["y"], sm.add_constant(reduced_X, has_constant="add")).fit(
        disp=False, method="newton", maxiter=500
    )
    stat = 2.0 * (full_model.llf - reduced.llf)
    return float(stat), float(chi2.sf(stat, 1))


def net_benefit(y, p, threshold):
    pred = np.asarray(p) >= threshold
    y = np.asarray(y, dtype=int)
    tp = np.sum(pred & (y == 1))
    fp = np.sum(pred & (y == 0))
    n = len(y)
    return tp / n - (fp / n) * threshold / (1.0 - threshold)


# -----------------------------------------------------------------------------
# Dynamic functional-form design used INSIDE repeated CV
# -----------------------------------------------------------------------------
def dynamic_train_design(train, bt_pvalues, alpha=BT_ALPHA):
    """
    Repeat the functional-form decision within a training sample.
    Any continuous predictor with Box-Tidwell p<alpha is represented with a
    three-knot RCS using knots estimated from that training sample.
    """
    X = pd.DataFrame(index=train.index)
    meta = {}

    # Age
    if bt_pvalues["Age"] < alpha:
        k = np.quantile(train["Age"], [0.10, 0.50, 0.90])
        b = rcs_basis(train["Age"], k)
        X["Age_linear"] = b[:, 0]
        X["Age_nonlinear"] = b[:, 1]
        meta["Age"] = k
    else:
        X["Age"] = train["Age"].to_numpy()
        meta["Age"] = None

    X["Male"] = train["male"].to_numpy()

    # Parasite density
    if bt_pvalues["Parasite density"] < alpha:
        k = np.quantile(train["parasite_1000"], [0.10, 0.50, 0.90])
        b = rcs_basis(train["parasite_1000"], k)
        X["Parasite_linear"] = b[:, 0]
        X["Parasite_nonlinear"] = b[:, 1]
        meta["Parasite density"] = k
    else:
        X["Parasite_per_1000"] = train["parasite_1000"].to_numpy()
        meta["Parasite density"] = None

    # Haemoglobin
    if bt_pvalues["Haemoglobin"] < alpha:
        k = np.quantile(train["Haemoglobin_g_dL"], [0.10, 0.50, 0.90])
        b = rcs_basis(train["Haemoglobin_g_dL"], k)
        X["Hb_linear"] = b[:, 0]
        X["Hb_nonlinear"] = b[:, 1]
        meta["Haemoglobin"] = k
    else:
        X["Hb"] = train["Haemoglobin_g_dL"].to_numpy()
        meta["Haemoglobin"] = None

    # Platelets
    if bt_pvalues["Platelet count"] < alpha:
        k = np.quantile(train["Platelet_Count_10^9_L"], [0.10, 0.50, 0.90])
        b = rcs_basis(train["Platelet_Count_10^9_L"], k)
        X["Platelet_linear"] = b[:, 0]
        X["Platelet_nonlinear"] = b[:, 1]
        meta["Platelet count"] = k
    else:
        X["Platelet"] = train["Platelet_Count_10^9_L"].to_numpy()
        meta["Platelet count"] = None

    return X, meta


def dynamic_test_design(test, meta):
    X = pd.DataFrame(index=test.index)

    if meta["Age"] is not None:
        b = rcs_basis(test["Age"], meta["Age"])
        X["Age_linear"] = b[:, 0]
        X["Age_nonlinear"] = b[:, 1]
    else:
        X["Age"] = test["Age"].to_numpy()

    X["Male"] = test["male"].to_numpy()

    if meta["Parasite density"] is not None:
        b = rcs_basis(test["parasite_1000"], meta["Parasite density"])
        X["Parasite_linear"] = b[:, 0]
        X["Parasite_nonlinear"] = b[:, 1]
    else:
        X["Parasite_per_1000"] = test["parasite_1000"].to_numpy()

    if meta["Haemoglobin"] is not None:
        b = rcs_basis(test["Haemoglobin_g_dL"], meta["Haemoglobin"])
        X["Hb_linear"] = b[:, 0]
        X["Hb_nonlinear"] = b[:, 1]
    else:
        X["Hb"] = test["Haemoglobin_g_dL"].to_numpy()

    if meta["Platelet count"] is not None:
        b = rcs_basis(test["Platelet_Count_10^9_L"], meta["Platelet count"])
        X["Platelet_linear"] = b[:, 0]
        X["Platelet_nonlinear"] = b[:, 1]
    else:
        X["Platelet"] = test["Platelet_Count_10^9_L"].to_numpy()

    return X


# -----------------------------------------------------------------------------
# Read and validate dataset
# -----------------------------------------------------------------------------
required = [
    "Patient_ID",
    "Year",
    "Age",
    "Gender",
    "Haemoglobin_g_dL",
    "Parasite_Density_per_uL",
    "Platelet_Count_10^9_L",
    "Outcome",
]

if not DATA_PATH.exists():
    raise FileNotFoundError(
        f"Dataset not found: {DATA_PATH}\n"
        "Place nakuru_malaria_dataset_2020_2026.csv in the same folder as this script."
    )

df = pd.read_csv(DATA_PATH)
missing_cols = [c for c in required if c not in df.columns]
if missing_cols:
    raise ValueError(f"Missing required columns: {missing_cols}")

if df[required].isna().any().any():
    raise ValueError("The submitted analytical dataset contains missing values in required variables.")

if df["Patient_ID"].duplicated().any():
    raise ValueError("Patient_ID is not unique; one record per patient is required.")

if not set(df["Outcome"].unique()).issubset({"Survived", "Died"}):
    raise ValueError("Unexpected values in Outcome.")
if not set(df["Gender"].unique()).issubset({"Female", "Male"}):
    raise ValueError("Unexpected values in Gender.")

df["y"] = (df["Outcome"] == "Died").astype(int)
df["male"] = (df["Gender"] == "Male").astype(int)
df["parasite_1000"] = df["Parasite_Density_per_uL"] / 1000.0

# Software versions
versions = pd.DataFrame(
    {
        "Software": [
            "Python", "numpy", "pandas", "scipy", "statsmodels",
            "scikit-learn", "matplotlib"
        ],
        "Version": [
            sys.version.split()[0], np.__version__, pd.__version__, scipy.__version__,
            statsmodels.__version__, sklearn.__version__, matplotlib.__version__
        ],
    }
)
versions.to_csv(OUT / "software_versions.csv", index=False)

# -----------------------------------------------------------------------------
# Descriptive and diagnostic checks
# -----------------------------------------------------------------------------
summary = pd.DataFrame(
    {
        "Statistic": ["N", "Deaths", "Mortality prevalence"],
        "Value": [len(df), int(df["y"].sum()), df["y"].mean()],
    }
)
summary.to_csv(OUT / "sample_summary.csv", index=False)

vifs = vif_table(df)
vifs.to_csv(OUT / "vif_results.csv", index=False)

bt_table, bt_model = box_tidwell(df)
bt_table.to_csv(OUT / "box_tidwell_results.csv", index=False)

# -----------------------------------------------------------------------------
# Full-sample original linear and corrected fixed RCS models
# -----------------------------------------------------------------------------
K = knot_rule(df)
knots_table = pd.DataFrame(
    {
        "Predictor": ["Age", "Haemoglobin", "Platelet count"],
        "Knot_10th": [K["Age"][0], K["Hb"][0], K["Platelet"][0]],
        "Knot_50th": [K["Age"][1], K["Hb"][1], K["Platelet"][1]],
        "Knot_90th": [K["Age"][2], K["Hb"][2], K["Platelet"][2]],
    }
)
knots_table.to_csv(OUT / "rcs_knots.csv", index=False)

X_linear = linear_design(df)
X_rcs = fixed_rcs_design(df, K)

linear_model = sm.Logit(df["y"], sm.add_constant(X_linear, has_constant="add")).fit(
    disp=False, method="newton", maxiter=500
)
rcs_model = sm.Logit(df["y"], sm.add_constant(X_rcs, has_constant="add")).fit(
    disp=False, method="newton", maxiter=500
)

# Complete coefficient tables
def model_table(model):
    ci = model.conf_int()
    return pd.DataFrame(
        {
            "Term": model.params.index,
            "Coefficient": model.params.to_numpy(),
            "SE": model.bse.to_numpy(),
            "CI_low_beta": ci[0].to_numpy(),
            "CI_high_beta": ci[1].to_numpy(),
            "p_value": model.pvalues.to_numpy(),
        }
    )

linear_coef = model_table(linear_model)
# Add AORs to original linear table
linear_coef["AOR"] = np.exp(linear_coef["Coefficient"])
linear_coef["AOR_CI_low"] = np.exp(linear_coef["CI_low_beta"])
linear_coef["AOR_CI_high"] = np.exp(linear_coef["CI_high_beta"])
linear_coef.loc[linear_coef["Term"] == "const", ["AOR", "AOR_CI_low", "AOR_CI_high"]] = np.nan
linear_coef.to_csv(OUT / "linear_model_coefficients.csv", index=False)

rcs_coef = model_table(rcs_model)
rcs_coef.to_csv(OUT / "rcs_model_coefficients.csv", index=False)

p_linear_full = np.asarray(
    linear_model.predict(sm.add_constant(X_linear, has_constant="add"))
)
p_rcs_full = np.asarray(
    rcs_model.predict(sm.add_constant(X_rcs, has_constant="add"))
)

# Model comparison
lr_stat = 2.0 * (rcs_model.llf - linear_model.llf)
model_comparison = pd.DataFrame(
    {
        "Measure": [
            "Clinical predictors", "Predictor parameters", "Log-likelihood", "AIC",
            "McFadden pseudo-R2", "Apparent AUC", "Apparent Brier score"
        ],
        "Original_linear": [
            5, 5, linear_model.llf, linear_model.aic, linear_model.prsquared,
            roc_auc_score(df["y"], p_linear_full), brier_score_loss(df["y"], p_linear_full)
        ],
        "Corrected_RCS": [
            5, 8, rcs_model.llf, rcs_model.aic, rcs_model.prsquared,
            roc_auc_score(df["y"], p_rcs_full), brier_score_loss(df["y"], p_rcs_full)
        ],
    }
)
model_comparison.to_csv(OUT / "model_comparison.csv", index=False)

nested_lr = pd.DataFrame(
    {
        "Comparison": ["Corrected RCS vs original linear"],
        "Chi_square": [lr_stat],
        "df": [3],
        "p_value": [chi2.sf(lr_stat, 3)],
    }
)
nested_lr.to_csv(OUT / "nested_likelihood_ratio_test.csv", index=False)

nonlinear_tests = []
for label, col in [
    ("Age", "Age_nonlinear"),
    ("Haemoglobin", "Hb_nonlinear"),
    ("Platelet count", "Platelet_nonlinear"),
]:
    stat, pval = nonlinear_component_test(df, X_rcs, rcs_model, col)
    nonlinear_tests.append((label, stat, 1, pval))
pd.DataFrame(
    nonlinear_tests, columns=["Predictor", "Chi_square", "df", "p_value"]
).to_csv(OUT / "nonlinear_component_tests.csv", index=False)

# Hosmer-Lemeshow for corrected RCS model
hl_stat, hl_df, hl_p, hl_groups = hosmer_lemeshow(df["y"], p_rcs_full, groups=10)
hl_groups.to_csv(OUT / "hosmer_lemeshow_groups.csv", index=False)
pd.DataFrame(
    {"Chi_square": [hl_stat], "df": [hl_df], "p_value": [hl_p]}
).to_csv(OUT / "hosmer_lemeshow_test.csv", index=False)

# Formal sample-size / shrinkage assessment
required_n, required_events = riley_shrinkage_sample_size(
    predictor_parameters=8,
    prevalence=float(df["y"].mean()),
    anticipated_cs_r2=0.05,
    target_shrinkage=0.90,
)
sample_size_table = pd.DataFrame(
    {
        "Quantity": [
            "Predictor parameters", "Observed prevalence", "Anticipated Cox-Snell R2",
            "Target shrinkage", "Required sample size", "Required expected deaths",
            "Observed sample size", "Observed deaths", "Observed events per predictor parameter"
        ],
        "Value": [
            8, df["y"].mean(), 0.05, 0.90, required_n, required_events,
            len(df), df["y"].sum(), df["y"].sum() / 8.0
        ],
    }
)
sample_size_table.to_csv(OUT / "formal_sample_size_assessment.csv", index=False)

# Clinically interpretable spline contrasts
contrast_specs = [
    ("Age", 50, 44), ("Age", 60, 44), ("Age", 70, 44), ("Age", 80, 44),
    ("Haemoglobin", 5, 10.2), ("Haemoglobin", 7, 10.2),
    ("Haemoglobin", 8, 10.2), ("Haemoglobin", 9, 10.2),
    ("Platelet", 75, 198.5), ("Platelet", 100, 198.5),
    ("Platelet", 150, 198.5), ("Platelet", 300, 198.5),
]
contrast_table = pd.DataFrame(
    [spline_or_contrast(rcs_model, K, *spec) for spec in contrast_specs]
)
contrast_table.to_csv(OUT / "spline_OR_contrasts.csv", index=False)

# -----------------------------------------------------------------------------
# Original development-sample thresholds (original linear model)
# -----------------------------------------------------------------------------
orig_youden = youden_threshold(df["y"], p_linear_full)
original_threshold_table = pd.DataFrame(
    [
        classification_metrics(df["y"], p_linear_full, 0.50),
        classification_metrics(df["y"], p_linear_full, orig_youden),
    ]
)
original_threshold_table.to_csv(OUT / "original_linear_thresholds.csv", index=False)

# -----------------------------------------------------------------------------
# Bootstrap optimism correction for the FIXED corrected RCS model
# (conditional on the corrected full-data specification; knots re-estimated)
# -----------------------------------------------------------------------------
rng = np.random.default_rng(SEED)
auc_apparent = roc_auc_score(df["y"], p_rcs_full)
brier_apparent = brier_score_loss(df["y"], p_rcs_full)
boot_rows = []

for b in range(BOOTSTRAP_REPS):
    idx = rng.integers(0, len(df), len(df))
    db = df.iloc[idx].reset_index(drop=True)
    kb = knot_rule(db)
    Xb = fixed_rcs_design(db, kb)
    Xo = fixed_rcs_design(df, kb)
    mb = sm.Logit(db["y"], sm.add_constant(Xb, has_constant="add")).fit(
        disp=False, method="newton", maxiter=200
    )
    p_b = np.asarray(mb.predict(sm.add_constant(Xb, has_constant="add")))
    p_o = np.asarray(mb.predict(sm.add_constant(Xo, has_constant="add")))

    auc_b = roc_auc_score(db["y"], p_b)
    auc_o = roc_auc_score(df["y"], p_o)
    br_b = brier_score_loss(db["y"], p_b)
    br_o = brier_score_loss(df["y"], p_o)
    cal_i, cal_s = calibration_stats(df["y"], p_o)

    boot_rows.append((auc_b, auc_o, br_b, br_o, cal_i, cal_s))

bootstrap = pd.DataFrame(
    boot_rows,
    columns=[
        "AUC_bootstrap_apparent", "AUC_test_original", "Brier_bootstrap_apparent",
        "Brier_test_original", "Calibration_intercept", "Calibration_slope"
    ],
)
bootstrap.to_csv(OUT / "bootstrap_1000_replications.csv", index=False)

auc_corrected_each = auc_apparent - (
    bootstrap["AUC_bootstrap_apparent"] - bootstrap["AUC_test_original"]
)
brier_corrected_each = brier_apparent - (
    bootstrap["Brier_bootstrap_apparent"] - bootstrap["Brier_test_original"]
)

bootstrap_summary = pd.DataFrame(
    {
        "Measure": [
            "Apparent AUC", "Optimism-corrected AUC", "Apparent Brier score",
            "Optimism-corrected Brier score", "Calibration intercept", "Calibration slope"
        ],
        "Estimate": [
            auc_apparent,
            auc_corrected_each.mean(),
            brier_apparent,
            brier_corrected_each.mean(),
            bootstrap["Calibration_intercept"].mean(),
            bootstrap["Calibration_slope"].mean(),
        ],
        "CI_low": [
            np.nan,
            np.quantile(auc_corrected_each, 0.025),
            np.nan,
            np.quantile(brier_corrected_each, 0.025),
            np.quantile(bootstrap["Calibration_intercept"], 0.025),
            np.quantile(bootstrap["Calibration_slope"], 0.025),
        ],
        "CI_high": [
            np.nan,
            np.quantile(auc_corrected_each, 0.975),
            np.nan,
            np.quantile(brier_corrected_each, 0.975),
            np.quantile(bootstrap["Calibration_intercept"], 0.975),
            np.quantile(bootstrap["Calibration_slope"], 0.975),
        ],
    }
)
bootstrap_summary.to_csv(OUT / "bootstrap_validation_summary.csv", index=False)

# -----------------------------------------------------------------------------
# Repeated 20 x 10-fold CV repeating ALL training-data modeling decisions:
# Box-Tidwell -> nonlinear/linear choice -> knots -> fit -> Youden threshold
# -----------------------------------------------------------------------------
rkf = RepeatedStratifiedKFold(
    n_splits=CV_FOLDS, n_repeats=CV_REPEATS, random_state=SEED
)

cv_pred_sum = np.zeros(len(df))
cv_pred_count = np.zeros(len(df))
linear_pred_sum = np.zeros(len(df))
linear_pred_count = np.zeros(len(df))
cv_fold_rows = []
selection_rows = []

for fold_id, (tr, te) in enumerate(rkf.split(df, df["y"])):
    repeat_id = fold_id // CV_FOLDS + 1
    within_repeat_fold = fold_id % CV_FOLDS + 1
    train = df.iloc[tr].copy()
    test = df.iloc[te].copy()

    # Repeat Box-Tidwell functional-form assessment in training data only
    bt_train, _ = box_tidwell(train)
    bt_p = dict(zip(bt_train["Predictor"], bt_train["p_value"]))

    # Build model from the training-data functional-form decisions
    Xtr, meta = dynamic_train_design(train, bt_p, alpha=BT_ALPHA)
    Xte = dynamic_test_design(test, meta)
    model = sm.Logit(train["y"], sm.add_constant(Xtr, has_constant="add")).fit(
        disp=False, method="newton", maxiter=200
    )
    p_tr = np.asarray(model.predict(sm.add_constant(Xtr, has_constant="add")))
    p_te = np.asarray(model.predict(sm.add_constant(Xte, has_constant="add")))

    # Repeat threshold selection in training data only
    threshold = youden_threshold(train["y"], p_tr)
    cm = classification_metrics(test["y"], p_te, threshold)

    cv_pred_sum[te] += p_te
    cv_pred_count[te] += 1

    # Cross-validated original linear comparator using same splits
    Xlin_tr = linear_design(train)
    Xlin_te = linear_design(test)
    lin = sm.Logit(train["y"], sm.add_constant(Xlin_tr, has_constant="add")).fit(
        disp=False, method="newton", maxiter=200
    )
    p_lin_te = np.asarray(lin.predict(sm.add_constant(Xlin_te, has_constant="add")))
    linear_pred_sum[te] += p_lin_te
    linear_pred_count[te] += 1

    row = {
        "Repeat": repeat_id,
        "Fold": within_repeat_fold,
        "Youden_threshold_training": threshold,
        **cm,
    }
    cv_fold_rows.append(row)

    selection_rows.append(
        {
            "Repeat": repeat_id,
            "Fold": within_repeat_fold,
            "Age_nonlinear": int(bt_p["Age"] < BT_ALPHA),
            "Parasite_nonlinear": int(bt_p["Parasite density"] < BT_ALPHA),
            "Haemoglobin_nonlinear": int(bt_p["Haemoglobin"] < BT_ALPHA),
            "Platelet_nonlinear": int(bt_p["Platelet count"] < BT_ALPHA),
            "BT_p_Age": bt_p["Age"],
            "BT_p_Parasite": bt_p["Parasite density"],
            "BT_p_Haemoglobin": bt_p["Haemoglobin"],
            "BT_p_Platelet": bt_p["Platelet count"],
        }
    )

cv_folds = pd.DataFrame(cv_fold_rows)
cv_folds.to_csv(OUT / "repeated_cv_fold_results.csv", index=False)
selection_df = pd.DataFrame(selection_rows)
selection_df.to_csv(OUT / "repeated_cv_functional_form_selection.csv", index=False)

p_cv = cv_pred_sum / cv_pred_count
p_linear_cv = linear_pred_sum / linear_pred_count

cv_cil, cv_slope = calibration_stats(df["y"], p_cv)
cv_summary = pd.DataFrame(
    {
        "Measure": [
            "Repeated-CV AUC", "Repeated-CV Brier score", "Calibration intercept",
            "Calibration slope", "Median training-fold Youden threshold",
            "IQR low", "IQR high", "Mean held-out sensitivity", "Mean held-out specificity",
            "Mean held-out PPV", "Mean held-out NPV", "Mean held-out accuracy"
        ],
        "Value": [
            roc_auc_score(df["y"], p_cv),
            brier_score_loss(df["y"], p_cv),
            cv_cil,
            cv_slope,
            cv_folds["Youden_threshold_training"].median(),
            cv_folds["Youden_threshold_training"].quantile(0.25),
            cv_folds["Youden_threshold_training"].quantile(0.75),
            cv_folds["Sensitivity"].mean(),
            cv_folds["Specificity"].mean(),
            cv_folds["PPV"].mean(),
            cv_folds["NPV"].mean(),
            cv_folds["Accuracy"].mean(),
        ],
    }
)
cv_summary.to_csv(OUT / "repeated_cv_summary.csv", index=False)

selection_summary = pd.DataFrame(
    {
        "Predictor": ["Age", "Parasite density", "Haemoglobin", "Platelet count"],
        "Proportion_selected_nonlinear": [
            selection_df["Age_nonlinear"].mean(),
            selection_df["Parasite_nonlinear"].mean(),
            selection_df["Haemoglobin_nonlinear"].mean(),
            selection_df["Platelet_nonlinear"].mean(),
        ],
    }
)
selection_summary.to_csv(OUT / "repeated_cv_functional_form_selection_summary.csv", index=False)

# Fixed 5% and 50% thresholds on averaged held-out predictions
cv_thresholds = pd.DataFrame(
    [
        classification_metrics(df["y"], p_cv, 0.05),
        classification_metrics(df["y"], p_cv, 0.50),
    ]
)
cv_thresholds.to_csv(OUT / "cross_validated_fixed_thresholds.csv", index=False)

# -----------------------------------------------------------------------------
# Ridge-penalized corrected RCS model: nested 10-fold CV
# -----------------------------------------------------------------------------
ridge_pipe = Pipeline(
    [
        ("scale", StandardScaler()),
        ("model", LogisticRegression(solver="lbfgs", penalty="l2", max_iter=5000)),
    ]
)

outer = StratifiedKFold(n_splits=10, shuffle=True, random_state=SEED)
p_ridge = np.zeros(len(df))
selected_C = []

for i, (tr, te) in enumerate(outer.split(df, df["y"])):
    train = df.iloc[tr].copy()
    test = df.iloc[te].copy()
    k = knot_rule(train)
    Xtr = fixed_rcs_design(train, k)
    Xte = fixed_rcs_design(test, k)

    inner = StratifiedKFold(n_splits=5, shuffle=True, random_state=SEED + i)
    search = GridSearchCV(
        ridge_pipe,
        {"model__C": C_GRID},
        scoring="neg_log_loss",
        cv=inner,
        n_jobs=1,
    )
    search.fit(Xtr, train["y"])
    p_ridge[te] = search.predict_proba(Xte)[:, 1]
    selected_C.append(search.best_params_["model__C"])

ridge_cil, ridge_slope = calibration_stats(df["y"], p_ridge)
ridge_summary = pd.DataFrame(
    {
        "Measure": ["Median selected C", "AUC", "Brier score", "Calibration intercept", "Calibration slope"],
        "Value": [
            np.median(selected_C), roc_auc_score(df["y"], p_ridge),
            brier_score_loss(df["y"], p_ridge), ridge_cil, ridge_slope
        ],
    }
)
ridge_summary.to_csv(OUT / "ridge_nested_cv_summary.csv", index=False)
pd.DataFrame({"Outer_fold": np.arange(1, 11), "Selected_C": selected_C}).to_csv(
    OUT / "ridge_selected_C_by_outer_fold.csv", index=False
)

# -----------------------------------------------------------------------------
# Figures
# -----------------------------------------------------------------------------
# 1. Original linear-model forest plot
terms_for_forest = ["Age", "Male", "Parasite_per_1000", "Hb", "Platelet"]
labels_for_forest = [
    "Age (per year)", "Sex (Male vs Female)", "Parasite density (per 1,000/µL)",
    "Haemoglobin (per g/dL)", "Platelet count (per ×10⁹/L)"
]
lt = linear_coef.set_index("Term")
ors = lt.loc[terms_for_forest, "AOR"].to_numpy()
lo = lt.loc[terms_for_forest, "AOR_CI_low"].to_numpy()
hi = lt.loc[terms_for_forest, "AOR_CI_high"].to_numpy()
ypos = np.arange(len(terms_for_forest))[::-1]
fig, ax = plt.subplots(figsize=(6.2, 4.3))
ax.errorbar(ors, ypos, xerr=[ors - lo, hi - ors], fmt="o", capsize=3)
ax.axvline(1.0, linestyle=":", linewidth=1)
ax.set_xscale("log")
ax.set_yticks(ypos)
ax.set_yticklabels(labels_for_forest)
ax.set_xlabel("Adjusted odds ratio (95% CI)")
fig.tight_layout()
fig.savefig(OUT / "Figure_1_Original_Linear_Forest.png", dpi=300, bbox_inches="tight")
plt.close(fig)

# Fit a full-data ridge model (tuned internally) for dashed spline sensitivity curves
inner_full = StratifiedKFold(n_splits=5, shuffle=True, random_state=SEED)
search_full = GridSearchCV(
    ridge_pipe, {"model__C": C_GRID}, scoring="neg_log_loss", cv=inner_full, n_jobs=1
)
search_full.fit(X_rcs, df["y"])

# Helper to get unpenalized OR curve + CI using beta covariance
def unpenalized_or_curve(variable, grid, ref):
    rows = [spline_or_contrast(rcs_model, K, variable, float(x), ref) for x in grid]
    d = pd.DataFrame(rows)
    return d["OR"].to_numpy(), d["CI_low"].to_numpy(), d["CI_high"].to_numpy()

# Ridge OR curve by constructing full design rows; other predictors cancel in contrasts.
def ridge_or_curve(variable, grid, ref):
    base = pd.DataFrame(
        {
            "Age": np.repeat(float(df["Age"].median()), len(grid)),
            "male": np.repeat(0, len(grid)),
            "parasite_1000": np.repeat(float(df["parasite_1000"].median()), len(grid)),
            "Haemoglobin_g_dL": np.repeat(float(df["Haemoglobin_g_dL"].median()), len(grid)),
            "Platelet_Count_10^9_L": np.repeat(float(df["Platelet_Count_10^9_L"].median()), len(grid)),
        }
    )
    if variable == "Age":
        base["Age"] = grid
    elif variable == "Haemoglobin":
        base["Haemoglobin_g_dL"] = grid
    elif variable == "Platelet":
        base["Platelet_Count_10^9_L"] = grid

    Xg = fixed_rcs_design(base, K)

    ref_df = base.iloc[[0]].copy()
    if variable == "Age":
        ref_df["Age"] = ref
    elif variable == "Haemoglobin":
        ref_df["Haemoglobin_g_dL"] = ref
    elif variable == "Platelet":
        ref_df["Platelet_Count_10^9_L"] = ref
    Xref = fixed_rcs_design(ref_df, K)

    # Pipeline decision_function is logit for binary logistic regression.
    eta_g = search_full.best_estimator_.decision_function(Xg)
    eta_ref = float(search_full.best_estimator_.decision_function(Xref)[0])
    return np.exp(eta_g - eta_ref)

plot_specs = [
    ("Age", np.linspace(df["Age"].min(), df["Age"].max(), 250), 44.0, K["Age"], "Age (years)", "Figure_2_Age_RCS.png"),
    ("Haemoglobin", np.linspace(df["Haemoglobin_g_dL"].min(), df["Haemoglobin_g_dL"].max(), 250), 10.2, K["Hb"], "Haemoglobin (g/dL)", "Figure_3_Haemoglobin_RCS.png"),
    ("Platelet", np.linspace(df["Platelet_Count_10^9_L"].min(), df["Platelet_Count_10^9_L"].max(), 250), 198.5, K["Platelet"], "Platelet count (×10⁹/L)", "Figure_4_Platelet_RCS.png"),
]

for variable, grid, ref, knots, xlabel, filename in plot_specs:
    est, low, high = unpenalized_or_curve(variable, grid, ref)
    ridge_curve = ridge_or_curve(variable, grid, ref)
    fig, ax = plt.subplots(figsize=(6.3, 4.4))
    ax.plot(grid, est, linewidth=2, label="Restricted cubic spline")
    ax.fill_between(grid, low, high, alpha=0.20, label="95% CI")
    ax.plot(grid, ridge_curve, linestyle="--", linewidth=1.5, label="Ridge-penalized spline")
    ax.axhline(1.0, linestyle=":", linewidth=1)
    for knot in knots:
        ax.axvline(knot, linestyle=":", linewidth=0.8)
    ax.set_yscale("log")
    ax.set_xlabel(xlabel)
    ax.set_ylabel("Adjusted odds ratio (log scale)")
    ax.legend(frameon=False)
    fig.tight_layout()
    fig.savefig(OUT / filename, dpi=300, bbox_inches="tight")
    plt.close(fig)

# Smooth calibration curve from repeated-CV predictions
smooth = lowess(df["y"].to_numpy(), p_cv, frac=0.40, it=0, return_sorted=True)
smooth[:, 1] = np.clip(smooth[:, 1], 0, 1)
fig, ax = plt.subplots(figsize=(6.2, 4.2))
ax.plot([0, 0.8], [0, 0.8], linestyle=":", linewidth=1.2, label="Ideal")
ax.plot(smooth[:, 0], smooth[:, 1], linewidth=2, label="Cross-validated smooth calibration")
ax.set_xlim(0, 0.8)
ax.set_ylim(0, 0.8)
ax.set_xlabel("Predicted mortality probability")
ax.set_ylabel("Observed mortality probability")
ax.legend(frameon=False)
fig.tight_layout()
fig.savefig(OUT / "Figure_5_Smooth_CrossValidated_Calibration.png", dpi=300, bbox_inches="tight")
plt.close(fig)

# ROC comparison
fig, ax = plt.subplots(figsize=(7.0, 5.2))
for label, pred in [
    ("Linear", p_linear_cv),
    ("RCS", p_cv),
    ("Ridge RCS", p_ridge),
]:
    fpr, tpr, _ = roc_curve(df["y"], pred)
    auc = roc_auc_score(df["y"], pred)
    ax.plot(fpr, tpr, linewidth=2, label=f"{label} (AUC = {auc:.3f})")
ax.plot([0, 1], [0, 1], linestyle=":", linewidth=1.2)
ax.set_xlabel("False positive rate (1 − specificity)")
ax.set_ylabel("True positive rate (sensitivity)")
ax.legend(frameon=False, loc="lower right")
fig.tight_layout()
fig.savefig(OUT / "Figure_6_CrossValidated_ROC_Comparison.png", dpi=300, bbox_inches="tight")
plt.close(fig)

# Decision curve using repeated-CV predictions
threshold_grid = np.linspace(0.01, 0.20, 100)
prevalence = df["y"].mean()
nb_model = np.array([net_benefit(df["y"], p_cv, t) for t in threshold_grid])
nb_all = prevalence - (1 - prevalence) * threshold_grid / (1 - threshold_grid)
nb_none = np.zeros_like(threshold_grid)

dca = pd.DataFrame(
    {
        "Threshold": threshold_grid,
        "Model_net_benefit": nb_model,
        "Classify_all_net_benefit": nb_all,
        "Classify_none_net_benefit": nb_none,
    }
)
dca.to_csv(OUT / "decision_curve_values.csv", index=False)

fig, ax = plt.subplots(figsize=(6.2, 4.2))
ax.plot(threshold_grid, nb_model, linewidth=2, label="RCS model")
ax.plot(threshold_grid, nb_all, linestyle="--", linewidth=1.5, label="Classify all")
ax.plot(threshold_grid, nb_none, linestyle=":", linewidth=1.5, label="Classify none")
ax.set_xlabel("Risk threshold")
ax.set_ylabel("Net benefit")
ax.legend(frameon=False)
fig.tight_layout()
fig.savefig(OUT / "Figure_7_Decision_Curve.png", dpi=300, bbox_inches="tight")
plt.close(fig)

# -----------------------------------------------------------------------------
# Plain-text manuscript-oriented summary
# -----------------------------------------------------------------------------
summary_lines = []
summary_lines.append("REPRODUCIBLE ANALYSIS SUMMARY")
summary_lines.append("=" * 80)
summary_lines.append(f"N = {len(df)}; deaths = {int(df['y'].sum())}; prevalence = {df['y'].mean():.4f}")
summary_lines.append("")
summary_lines.append("VIF range (five original predictors): " +
                     f"{vifs['VIF'].min():.3f} to {vifs['VIF'].max():.3f}")
summary_lines.append("Box-Tidwell p-values:")
for _, r in bt_table.iterrows():
    summary_lines.append(f"  {r['Predictor']}: p = {r['p_value']:.6g}")
summary_lines.append("")
summary_lines.append(
    f"Original linear model: LL={linear_model.llf:.3f}, AIC={linear_model.aic:.3f}, "
    f"McFadden R2={linear_model.prsquared:.4f}, AUC={roc_auc_score(df['y'], p_linear_full):.3f}, "
    f"Brier={brier_score_loss(df['y'], p_linear_full):.4f}"
)
summary_lines.append(
    f"Corrected RCS model: LL={rcs_model.llf:.3f}, AIC={rcs_model.aic:.3f}, "
    f"McFadden R2={rcs_model.prsquared:.4f}, AUC={roc_auc_score(df['y'], p_rcs_full):.3f}, "
    f"Brier={brier_score_loss(df['y'], p_rcs_full):.4f}"
)
summary_lines.append(
    f"Nested LR comparison: chi2(3)={lr_stat:.3f}, p={chi2.sf(lr_stat,3):.6g}"
)
summary_lines.append(
    f"Hosmer-Lemeshow corrected RCS: chi2({hl_df})={hl_stat:.3f}, p={hl_p:.3f}"
)
summary_lines.append(
    f"Formal shrinkage sample-size criterion: n={required_n:.1f}, expected deaths={required_events:.1f}; "
    f"observed n={len(df)}, deaths={int(df['y'].sum())}"
)
summary_lines.append("")
summary_lines.append("Bootstrap validation (1,000 resamples):")
for _, r in bootstrap_summary.iterrows():
    if pd.notna(r["CI_low"]):
        summary_lines.append(
            f"  {r['Measure']}: {r['Estimate']:.4f} (95% interval {r['CI_low']:.4f} to {r['CI_high']:.4f})"
        )
    else:
        summary_lines.append(f"  {r['Measure']}: {r['Estimate']:.4f}")
summary_lines.append("")
summary_lines.append("Repeated 20 x 10-fold CV with functional-form selection repeated inside training folds:")
summary_lines.append(f"  AUC={roc_auc_score(df['y'], p_cv):.4f}")
summary_lines.append(f"  Brier={brier_score_loss(df['y'], p_cv):.4f}")
summary_lines.append(f"  Calibration intercept={cv_cil:.4f}")
summary_lines.append(f"  Calibration slope={cv_slope:.4f}")
summary_lines.append(
    "  Nonlinear selection frequencies: " + ", ".join(
        f"{row.Predictor}={row.Proportion_selected_nonlinear:.3f}"
        for row in selection_summary.itertuples()
    )
)
summary_lines.append(
    f"  Training-fold Youden median={cv_folds['Youden_threshold_training'].median():.4f}; "
    f"IQR={cv_folds['Youden_threshold_training'].quantile(.25):.4f}–{cv_folds['Youden_threshold_training'].quantile(.75):.4f}"
)
summary_lines.append("")
summary_lines.append("Nested ridge RCS sensitivity analysis:")
summary_lines.append(f"  Median selected C={np.median(selected_C):.4g}")
summary_lines.append(f"  AUC={roc_auc_score(df['y'], p_ridge):.4f}")
summary_lines.append(f"  Brier={brier_score_loss(df['y'], p_ridge):.4f}")
summary_lines.append(f"  Calibration intercept={ridge_cil:.4f}")
summary_lines.append(f"  Calibration slope={ridge_slope:.4f}")

(OUT / "analysis_summary.txt").write_text("\n".join(summary_lines), encoding="utf-8")

print("\n".join(summary_lines))
print(f"\nAll tables and figures written to: {OUT}")