00_ML_HCV_SHORT_TX.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c8b9924b",
   "metadata": {},
   "source": [
    "# __Machine learning prediction of short duration HCV treatment outcome__\n",
    "Exploration of  machine learning (ML) to develop a clinical tool that predicts the outcomes of short duration direct-acting antiviral (DAA) therapy in patients with Hepatitis C Virus (HCV). Determine whether an individual treated with short duration DAAs will achieve cure using baseline demographic, clinical, and treatment-related variables factors."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de43d195",
   "metadata": {},
   "source": [
    "# 1.0 __Data processing__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ce401c4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Import libraries for data manipulation\n",
    "import numpy as np\n",
    "from numpy import mean, std, argmax\n",
    "import pandas as pd\n",
    "from collections import defaultdict\n",
    "import scipy.stats as stats \n",
    "from scipy.stats import mannwhitneyu, loguniform  \n",
    "from scipy.signal import savgol_filter\n",
    "import statsmodels.api as sm\n",
    "from sklearn.experimental import enable_iterative_imputer\n",
    "from sklearn.impute import IterativeImputer\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "# Import libraries for data tables and visualization \n",
    "from tableone import TableOne\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import colormaps\n",
    "import seaborn as sns\n",
    "import os\n",
    "%matplotlib inline\n",
    "\n",
    "# Import libraries for the machine learning\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.utils import resample\n",
    "from imblearn.over_sampling import BorderlineSMOTE\n",
    "from sklearn.model_selection import cross_val_predict, KFold, RandomizedSearchCV #or GridSearchCV \n",
    "from sklearn.metrics import (confusion_matrix, accuracy_score, f1_score, precision_score, average_precision_score, recall_score, \n",
    "                             roc_auc_score, roc_curve, auc, precision_recall_curve,  matthews_corrcoef, brier_score_loss)\n",
    "from sklearn.calibration import calibration_curve, CalibratedClassifierCV\n",
    "import xgboost as xgb\n",
    "import shap\n",
    "\n",
    "# Suppress warnings\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=FutureWarning)  \n",
    "warnings.filterwarnings(\"ignore\", category=UserWarning)    \n",
    "warnings.filterwarnings(\"ignore\", category=RuntimeWarning)  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7c001b7",
   "metadata": {},
   "source": [
    "Import data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82e512f3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df=pd.read_csv(\"ML_DAA_short.csv\")\n",
    "# set id as index\n",
    "df.set_index('sid', inplace=True)\n",
    "# check missing data\n",
    "missing_data = df.isnull().sum()\n",
    "print(missing_data[missing_data != 0]) # <5% missing data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2245443",
   "metadata": {},
   "source": [
    "Impute small number of missing liver enzyme values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d555c7a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# use multiple imputation by chained equations (MICE) - best for missing data  with complex relationships between features\n",
    "imputer = IterativeImputer()\n",
    "df[['fibro', 'alt', 'ast', 'plt']] = imputer.fit_transform(df[['fibro', 'alt', 'ast', 'plt']])\n",
    "# CIRRHOSIS\n",
    "df['cirrhosis'] = (df['fibro'] > 12.5).astype(int)\n",
    "## AST-ALT RATIO\n",
    "df['ast_alt_ratio'] = df['ast'] / df['alt']\n",
    "# APRI \n",
    "ULN_AST = 40 # Upper Limit of Normal AST (ULN_AST)\n",
    "df['apri'] = (df['ast'] / ULN_AST) * 100 / df['plt']\n",
    "# view scores\n",
    "df[['ast','alt', 'plt', 'ast_alt_ratio', 'apri']].describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bf522d3",
   "metadata": {},
   "source": [
    "Standardise virology"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdc7a6dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# GENOTYPE\n",
    "# Create the binary variables for GT_1 and GT_3\n",
    "df['GT_1'] = np.where(df['genotype'] == 1, 1, 0)\n",
    "df['GT_3'] = np.where(df['genotype'] == 3, 1, 0)\n",
    "# apply the same log10 scale to normalise\n",
    "df['log_rna'] = df['rna', 'duration', 'age'].apply(np.log10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95acf332",
   "metadata": {},
   "source": [
    "DAA classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a0988f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the new variable DAA class\n",
    "df['PI'] = np.where(df['daa'].str.contains('GLE|GRZ|PrOD', regex=True), 1, 0)\n",
    "df['NS5A'] = np.where(df['daa'].str.contains('LDV|DCV|PrOD|PIB|ELB|VEL', regex=True), 1, 0)\n",
    "df['NS5B'] = np.where(df['daa'].str.contains('SOF|PrOD', regex=True), 1, 0)\n",
    "# Define conditions and corresponding classifications\n",
    "conditions = [\n",
    "    (df['PI'] == 1) & (df['NS5A'] == 1) & (df['NS5B'] == 1),\n",
    "    (df['PI'] == 1) & (df['NS5A'] == 1),\n",
    "    (df['NS5A'] == 1) & (df['NS5B'] == 1)\n",
    "]\n",
    "choices = ['PI+NS5A+NS5B', 'PI+NS5A', 'NS5A+NS5B']\n",
    "# Use numpy.select to classify\n",
    "df['daa_class'] = np.select(conditions, choices, default='Other')\n",
    "df['daa_class'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd84556c",
   "metadata": {},
   "source": [
    "Clinical Trial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96dcbbba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the new variable 'clin_trial' vs. cohort\n",
    "df['clin_trial'] = np.where(df['study'].str.contains('TARGET|REACT|SMART', regex=True), 1, 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bc6c2c1",
   "metadata": {},
   "source": [
    "# 2.0 Descriptives"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37f2aac6",
   "metadata": {},
   "source": [
    "## 2.1 Descriptive Table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66c616e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create table one\n",
    "columns = ['age', 'sex',  'white', 'idu', 'prison', 'gbm',  'hiv', 'clin_trial', 'acute', 'genotype','log_rna', \n",
    "           'alt', 'ast', 'plt', 'ast_alt_ratio', 'apri', 'fibro', 'cirrhosis', 'daa_class', 'duration', ]\n",
    "# identify variables to be tabulated\n",
    "categorical = ['sex', 'white', 'idu', 'prison', 'gbm',  'hiv', 'clin_trial', 'acute',  'genotype', 'cirrhosis',  'daa_class']\n",
    "# identify variables for median calcluations\n",
    "nonnormal = ['age', 'log_rna', 'alt', 'ast', 'plt', 'ast_alt_ratio', 'apri', 'fibro',  'duration']\n",
    "# Create a TableOne instance, stratified by the 'fail' column\n",
    "t1 = TableOne(data=df, columns=columns, categorical=categorical, nonnormal=nonnormal, pval=False, missing=False, \n",
    "              row_percent=False, decimals=1, display_all=False, groupby='fail')\n",
    "# Print the table \n",
    "print(t1.tabulate(tablefmt=\"github\"))\n",
    "# Convert the TableOne object to a DataFrame\n",
    "t1_df = t1.tableone.reset_index()\n",
    "# Export the DataFrame to a CSV file\n",
    "t1_df.to_csv('python_output/MT1_demo_table.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd99b9f8",
   "metadata": {},
   "source": [
    "## 2.2 Violin plots\n",
    "\n",
    "Violin plots illustrate the distribution of HCV RNA levels and liver function markers (e.g., AST-ALT ratio, APRI, ALT, and AST) stratified by HCV phase (chronic vs. acute). Median values and variability within each group are shown, with p-values calculated using the Mann-Whitney U test to assess differences between groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d622b36",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of rows and columns for the subplot grid\n",
    "nrows = 2\n",
    "ncols = 2\n",
    "# Adjust these values to change the overall size of the plot\n",
    "fig_width = 9\n",
    "fig_height = 6\n",
    "\n",
    "# Map acute infection to readable labels\n",
    "df['HCV Phase'] = df['acute'].map({0: 'Chronic', 1: 'Acute'})\n",
    "# List of columns to plot\n",
    "columns_to_plot = [\"log_rna\", \"ast_alt_ratio\", \"apri\", \"fibro\"]\n",
    "# Mapping of original column names to the desired axis labels\n",
    "axis_labels = {\n",
    "    \"log_rna\": \"HCV RNA (log10 IU/mL)\",\n",
    "    \"ast_alt_ratio\": \"AST-ALT Ratio\",\n",
    "    \"apri\": \"APRI Score\",\n",
    "    \"fibro\": \"Fibrosis Score (kPa)\"\n",
    "}\n",
    "# List of categorical variables to loop through\n",
    "x_variables = ['HCV Phase']\n",
    "# Loop over categorical variables\n",
    "for x_var in x_variables:\n",
    "    # Create a figure with a specified size\n",
    "    fig, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(fig_width, fig_height))\n",
    "    # Flatten the axes array for easy indexing\n",
    "    axes = axes.flatten()\n",
    "    # Plot each column on its respective axes for the current x variable\n",
    "    for i, col in enumerate(columns_to_plot):\n",
    "        if i < len(axes):\n",
    "            # Create the violin plot\n",
    "            ax = sns.violinplot(x=x_var, y=col, data=df, palette=\"cool\", ax=axes[i], inner=\"box\", scale=\"width\")\n",
    "\n",
    "            # Adjust the transparency of the violin plot\n",
    "            for violin in ax.collections:  # Access the collection of polygons in the plot\n",
    "                if isinstance(violin, matplotlib.collections.PolyCollection):  # Only adjust polygons\n",
    "                    violin.set_alpha(0.7)\n",
    "            # Perform a statistical test (Mann-Whitney U test for non-parametric data)\n",
    "            acute_values = df[df['acute'] == 1][col]\n",
    "            chronic_values = df[df['acute'] == 0][col]\n",
    "            stat, p_value = stats.mannwhitneyu(acute_values, chronic_values)\n",
    "            # Format the p-value for display\n",
    "            if p_value < 0.001:\n",
    "                p_text = \"p<0.001\"\n",
    "            else:\n",
    "                p_text = f\"p={p_value:.3f}\"\n",
    "            # Set the title with the p-value\n",
    "            axes[i].set_title(p_text)\n",
    "            # Set y-axis labels for each plot\n",
    "            axes[i].set_ylabel(axis_labels[col])\n",
    "    # Add a shared x-axis label for the figure\n",
    "    fig.text(0.5, 0.04, \"HCV Phase\", ha='center', fontsize=12)\n",
    "    # Adjust layout for better display\n",
    "    plt.tight_layout(rect=[0, 0.05, 1, 1]) \n",
    "    # Save the figure\n",
    "    output_path = \"python_output/SF2_violin_plot.png\"\n",
    "    plt.savefig(output_path, dpi=300, bbox_inches=\"tight\")\n",
    "    print(f\"Figure saved to {output_path}\")\n",
    "    # Show the plot\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "896a3559",
   "metadata": {},
   "source": [
    "## 2.2 Scatter Plots \n",
    "Scatter plots depict correlations between liver function markers, with Pearson correlation coefficients (r) to quantify strength and direction of linear relationships (r=-1 indicates perfect negative correlation; r=0 indicates no correlation; r=1 indicates perfect positive correlation). Data points are color-coded by treatment response: blue dots represent sustained virological response, and purple dots represent treatment failure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1ed5bf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Set figure style\n",
    "sns.set_style('whitegrid')\n",
    "sns.set_context(\"paper\", font_scale=1.1)\n",
    "plt.rcParams['text.color'] = 'black'\n",
    "plt.rcParams['axes.labelcolor'] = 'black'\n",
    "plt.rcParams['xtick.color'] = 'black'\n",
    "plt.rcParams['ytick.color'] = 'black'\n",
    "plt.rcParams['axes.titlecolor'] = 'black'\n",
    "# Create a new column 'SVR' based on the 'fail' column\n",
    "df['SVR'] = df['fail'].map({0: 'Yes', 1: 'No'})  # 0 -> Yes (SVR), 1 -> No (Treatment failure)\n",
    "# Mapping of original column names to the desired axis labels\n",
    "axis_labels = {\n",
    "    \"apri\": \"APRI Score\",\n",
    "    \"ast_alt_ratio\": \"AST-ALT Ratio\",\n",
    "    \"fibro\": \"Fibrosis Score (kPa)\",\n",
    "    \"Treatment Failure\": \"Treatment Failure\"\n",
    "}\n",
    "# Calculate correlation coefficients\n",
    "pearson_corr_apri_ast_alt = df[['apri', 'ast_alt_ratio']].corr().iloc[0, 1]\n",
    "pearson_corr_apri_fibro = df[['apri', 'fibro']].corr().iloc[0, 1]\n",
    "pearson_corr_ast_alt_fibro = df[['ast_alt_ratio', 'fibro']].corr().iloc[0, 1]\n",
    "pearson_corr_ast_alt = df[['ast', 'alt']].corr().iloc[0, 1]\n",
    "# Create a figure with 2x2 subplots\n",
    "fig, axs = plt.subplots(2, 2, figsize=(9, 6))  # axs is a 2D array\n",
    "# Flatten the 2D array to easily access each subplot\n",
    "axs = axs.flatten()\n",
    "# Custom hue order and mapping\n",
    "hue_order = ['No', 'Yes']\n",
    "# Plot 0: ALT vs AST\n",
    "scatter0= sns.scatterplot(x='alt', y='ast', hue='SVR', data=df, ax=axs[0], \n",
    "                           palette='cool_r', alpha=0.6, hue_order=hue_order)\n",
    "axs[0].set_title(f'(r = {pearson_corr_ast_alt:.2f})')\n",
    "axs[0].set_xlabel(\"ALT\")  # Corrected x-axis label\n",
    "axs[0].set_ylabel(\"AST\")  # Corrected y-axis label\n",
    "# Plot 1: APRI vs AST/ALT Ratio\n",
    "scatter1 = sns.scatterplot(x='apri', y='ast_alt_ratio', hue='SVR', data=df, ax=axs[1], \n",
    "                           palette='cool_r', alpha=0.6, hue_order=hue_order)\n",
    "axs[1].set_title(f'(r = {pearson_corr_apri_ast_alt:.2f})')\n",
    "axs[1].set_xlabel(axis_labels['apri'])  # Rename x-axis\n",
    "axs[1].set_ylabel(axis_labels['ast_alt_ratio'])  # Rename y-axis\n",
    "# Plot 2: APRI vs Fibroscan\n",
    "scatter2 = sns.scatterplot(x='apri', y='fibro', hue='SVR', data=df, ax=axs[2], \n",
    "                           palette='cool_r', alpha=0.6, hue_order=hue_order)\n",
    "axs[2].set_title(f'(r = {pearson_corr_apri_fibro:.2f})')\n",
    "axs[2].set_xlabel(axis_labels['fibro'])  # Rename x-axis\n",
    "axs[2].set_ylabel(axis_labels['apri'])  # Rename y-axis\n",
    "# Plot 3: AST/ALT Ratio vs Fibroscan\n",
    "scatter3 = sns.scatterplot(x='ast_alt_ratio', y='fibro', hue='SVR', data=df, ax=axs[3], \n",
    "                           palette='cool_r', alpha=0.6, hue_order=hue_order)\n",
    "axs[3].set_title(f'(r = {pearson_corr_ast_alt_fibro:.2f})')\n",
    "axs[3].set_xlabel(axis_labels['ast_alt_ratio'])  # Rename x-axis\n",
    "axs[3].set_ylabel(axis_labels['fibro'])  # Rename y-axis\n",
    "# Adjust the layout to avoid overlap\n",
    "plt.tight_layout()\n",
    "# Save the figure\n",
    "output_path = \"python_output/SF1_scatter_plot.png\"\n",
    "plt.savefig(output_path, dpi=300, bbox_inches=\"tight\")\n",
    "print(f\"Figure saved to {output_path}\")\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91c88135",
   "metadata": {},
   "source": [
    "# 3.0 Logistic Regression"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81d37bd9",
   "metadata": {},
   "source": [
    "## 3.1 Data perparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7f86ec7",
   "metadata": {},
   "outputs": [],
   "source": [
    "df1 = df[['log_rna', 'duration', 'ast_alt_ratio', 'daa_class', 'apri', 'GT_3', 'fibro', 'age', 'white', 'sex', 'idu', 'fail']].copy()\n",
    "# get dummies for categorical variables\n",
    "df1 = pd.get_dummies(df1, columns=['daa_class'], drop_first=True)\n",
    "# Convert boolean columns to integers\n",
    "df1[['daa_class_PI+NS5A', 'daa_class_PI+NS5A+NS5B']] = \\\n",
    "    df1[['daa_class_PI+NS5A', 'daa_class_PI+NS5A+NS5B']].astype(int)\n",
    "# Check for missing values and fill if necessary\n",
    "df1 = df1.fillna(0)\n",
    "df1.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb534f81",
   "metadata": {},
   "source": [
    "## 3.2 Odds ratios"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e0eb372",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define feature matrix (X) and target variable (y)\n",
    "X = df1.drop(['fail'], axis=1)\n",
    "y = df1['fail']\n",
    "# Add constant column\n",
    "X = sm.add_constant(X)\n",
    "# Fit logistic regression model\n",
    "logit_model = sm.Logit(y, X)\n",
    "result = logit_model.fit()\n",
    "# Calculate odds ratios, confidence intervals, and p-values\n",
    "odds_ratios = np.exp(result.params)  # ORs\n",
    "conf = np.exp(result.conf_int())     # Confidence intervals\n",
    "p_values = result.pvalues            # p-values\n",
    "# Combine results into a single DataFrame\n",
    "summary_df = conf.copy()\n",
    "summary_df.columns = ['2.5%', '97.5%']  # Rename CI columns\n",
    "summary_df['OR'] = odds_ratios\n",
    "summary_df['p-value'] = p_values\n",
    "# Display the summary DataFrame\n",
    "print(\"Odds Ratios, Confidence Intervals, and p-values:\")\n",
    "print(summary_df)\n",
    "# Save the results to a CSV\n",
    "summary_df.to_csv('python_output/ST2_log_regression.csv', index=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fc22274",
   "metadata": {},
   "source": [
    "## 3.3 Data conversion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6ed0e9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df2 = df[['log_rna', 'duration', 'ast_alt_ratio', 'daa_class', 'GT_3', 'fail']].copy()\n",
    "df2 = pd.get_dummies(df2, columns=['daa_class'], drop_first=False)\n",
    "df2 = df2.drop(['daa_class_PI+NS5A+NS5B'], axis=1)\n",
    "# Convert boolean columns to integers\n",
    "df2[['daa_class_PI+NS5A', 'daa_class_NS5A+NS5B']] = df2[['daa_class_PI+NS5A', 'daa_class_NS5A+NS5B']].astype(int)\n",
    "df2.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92745ef3",
   "metadata": {},
   "source": [
    "## 3.4 Define metrics functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1374e63c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1) Set seed\n",
    "seed = 23\n",
    "np.random.seed(seed)\n",
    "# 2) Helper: calculate_metrics\n",
    "def calculate_metrics(y_true, y_pred, y_proba):\n",
    "    tn, fp, fn, tp = confusion_matrix(y_true, y_pred).ravel()\n",
    "    precision = tp / (tp + fp) if (tp + fp) > 0 else 0\n",
    "    recall = tp / (tp + fn) if (tp + fn) > 0 else 0\n",
    "    fmi = (precision * recall)**0.5 # Calculate Fowlkes-Mallows Index\n",
    "    return {\n",
    "        'accuracy':  accuracy_score(y_true, y_pred),\n",
    "        'recall':    recall,\n",
    "        'specificity': tn / (tn + fp) if (tn + fp) > 0 else 0,\n",
    "        'precision': precision,\n",
    "        'npv':       tn / (tn + fn) if (tn + fn) > 0 else 0,\n",
    "        'f1':        f1_score(y_true, y_pred),\n",
    "        'mcc':       matthews_corrcoef(y_true, y_pred),\n",
    "        'fmi':       fmi \n",
    "    }\n",
    "# 3) Aggregate threshold-based metrics\n",
    "def aggregate_metrics(metrics_list):\n",
    "    metrics_summary = {\n",
    "        metric: np.mean([fold[metric] for fold in metrics_list])\n",
    "        for metric in metrics_list[0]\n",
    "    }\n",
    "    metrics_std = {\n",
    "        metric: np.std([fold[metric] for fold in metrics_list])\n",
    "        for metric in metrics_list[0]\n",
    "    }\n",
    "    return metrics_summary, metrics_std\n",
    "# 3) Brier aggregator (if you want calibration curves)\n",
    "def aggregate_calibration_metrics(brier_scores_list, calibration_prob_true_list, calibration_prob_pred_list, num_points=100):\n",
    "    mean_brier = np.mean(brier_scores_list)\n",
    "    std_brier = np.std(brier_scores_list)\n",
    "\n",
    "    interpolated_true = []\n",
    "    interpolated_pred = []\n",
    "    for prob_true, prob_pred in zip(calibration_prob_true_list, calibration_prob_pred_list):\n",
    "        interp_pred = np.linspace(0, 1, num_points)\n",
    "        interp_true = np.interp(interp_pred, prob_pred, prob_true)\n",
    "        interpolated_true.append(interp_true)\n",
    "        interpolated_pred.append(interp_pred)\n",
    "    mean_prob_true = np.mean(np.array(interpolated_true), axis=0)\n",
    "    mean_prob_pred = np.mean(np.array(interpolated_pred), axis=0)\n",
    "    return mean_brier, std_brier, mean_prob_true, mean_prob_pred\n",
    "\n",
    "# 4) Overconfidence penalty (optional)\n",
    "def calculate_overconfidence_penalty(y_true, y_proba_1d):\n",
    "    penalties = [\n",
    "        (p - y)**2 * (2 if p > 0.85 or p < 0.15 else 1)\n",
    "        for p, y in zip(y_proba_1d, y_true)\n",
    "    ]\n",
    "    return np.mean(penalties)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2200e42",
   "metadata": {},
   "source": [
    "## 3.5 Nested cross validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0cccb930",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 0) Lists for storing final metrics across outer folds\n",
    "fpr_list, tpr_list, roc_auc_list = [], [], []\n",
    "precision_list, recall_list, pr_auc_list = [], [], []\n",
    "calibration_prob_true_list, calibration_prob_pred_list = [], []\n",
    "interpolated_prob_true_list, interpolated_prob_pred_list = [], []\n",
    "overconfidence_penalties, brier_scores = [], []\n",
    "default_metrics = []\n",
    "# Define a fixed grid for predicted probabilities\n",
    "fixed_grid = np.linspace(0, 1, 100)\n",
    "# 1) Prepare data\n",
    "X = df2.drop('fail', axis=1)\n",
    "y = df2['fail']\n",
    "# 2) Define outer CV\n",
    "outer_cv = KFold(n_splits=5, shuffle=True, random_state=seed)  # Fixed seed\n",
    "# 3) Outer CV loop\n",
    "for fold, (train_ix, test_ix) in enumerate(outer_cv.split(X)):\n",
    "    X_train, X_test = X.iloc[train_ix], X.iloc[test_ix]\n",
    "    y_train, y_test = y.iloc[train_ix], y.iloc[test_ix]\n",
    "    # 4) SMOTE on the train fold\n",
    "    smote = BorderlineSMOTE(random_state=seed)\n",
    "    X_train_res, y_train_res = smote.fit_resample(X_train, y_train)\n",
    "    # 5) Inner CV for hyperparameter tuning\n",
    "    cv_inner = KFold(n_splits=3, shuffle=True, random_state=seed)\n",
    "    # Define the logistic regression model (L2 penalty)\n",
    "    model = LogisticRegression(penalty='l2', fit_intercept=True, solver='lbfgs', max_iter=100, C=1.0, class_weight=None)\n",
    "    # Define hyperparameter search space (excluding l1_ratio)\n",
    "    space = {\n",
    "        'C': np.logspace(-4, 4, 20),\n",
    "        'max_iter': [50, 100, 250, 500, 750, 1000],\n",
    "        'class_weight': [None, 'balanced'],\n",
    "        'fit_intercept': [True, False],\n",
    "        'tol': np.logspace(-6, -1, 6),\n",
    "        'warm_start': [False] \n",
    "    }\n",
    "    # 6) Fit on inner train fold\n",
    "    log_model = model.fit(X_train_res, y_train_res)\n",
    "    # 7) Evaluate on outer test fold\n",
    "    # Predict probabilities\n",
    "    y_proba = log_model.predict_proba(X_test)\n",
    "    # Apply threshold\n",
    "    threshold = 0.5\n",
    "    yhat = (y_proba[:, 1] >= threshold).astype(int)\n",
    "    # ROC data for this fold\n",
    "    fpr, tpr, _ = roc_curve(y_test, y_proba[:, 1])\n",
    "    roc_auc = auc(fpr, tpr)\n",
    "    fpr_list.append(fpr)\n",
    "    tpr_list.append(tpr)\n",
    "    roc_auc_list.append(roc_auc)\n",
    "    # PRC data for this fold\n",
    "    precision, recall, _ = precision_recall_curve(y_test, y_proba[:, 1])\n",
    "    pr_auc = average_precision_score(y_test, y_proba[:, 1])\n",
    "    precision_list.append(precision)\n",
    "    recall_list.append(recall)\n",
    "    pr_auc_list.append(pr_auc)\n",
    "    # Performance metrics at 0.5 threshold\n",
    "    y_pred_0_5 = (y_proba[:, 1] >= 0.5).astype(int)\n",
    "    default_metrics.append(calculate_metrics(y_test, y_pred_0_5, y_proba))\n",
    "    # ii) Calibration\n",
    "    prob_true, prob_pred = calibration_curve(y_test, y_proba[:, 1], n_bins=11, strategy='uniform')\n",
    "    calibration_prob_true_list.append(prob_true)\n",
    "    calibration_prob_pred_list.append(prob_pred)\n",
    "    fold_brier_score = brier_score_loss(y_test, y_proba[:, 1])\n",
    "    brier_scores.append(fold_brier_score)\n",
    "    # Interpolation for consistent grid\n",
    "    interp_prob_true = np.interp(fixed_grid, prob_pred, prob_true, left=0, right=1)\n",
    "    interpolated_prob_true_list.append(interp_prob_true)\n",
    "    interpolated_prob_pred_list.append(fixed_grid)\n",
    "    # iii) Overconfidence\n",
    "    overconfidence_penalties.append(calculate_overconfidence_penalty(y_test, y_proba[:, 1]))\n",
    "# Aggregate calibration data\n",
    "mean_brier, std_brier, mean_prob_true,  mean_prob_pred = aggregate_calibration_metrics(\n",
    "    brier_scores, calibration_prob_true_list, calibration_prob_pred_list\n",
    ")\n",
    "# Aggregate threshold-based metrics\n",
    "log_metrics_summary, log_metrics_std= aggregate_metrics(default_metrics)\n",
    "# Print results\n",
    "print(f\"AUROC: {np.mean(roc_auc_list):.4f} ± {np.std(roc_auc_list):.4f}\")\n",
    "print(f\"AUPRC: {np.mean(pr_auc_list):.4f} ± {np.std(pr_auc_list):.4f}\")\n",
    "print(f\"Brier Score: {mean_brier:.4f} ± {std_brier:.4f}\")\n",
    "print(f\"Overconfidence Penalty: {np.mean(overconfidence_penalties):.4f} ± {np.std(overconfidence_penalties):.4f}\")\n",
    "print(\"\\nMetrics at 0.5 Threshold:\")\n",
    "for m, val in log_metrics_summary.items():\n",
    "    print(f\"{m.capitalize()}: {val:.4f} ± {log_metrics_std[m]:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61ac14a7",
   "metadata": {},
   "source": [
    "## 3.6 Define smoothing function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b68c5382",
   "metadata": {},
   "outputs": [],
   "source": [
    "def smooth_curve_savgol(y, window_length=7, polyorder=1, mode='interp', cval=0.0):\n",
    "    \"\"\"\n",
    "    Apply Savitzky–Golay smoothing to array y with robust handling of edge cases.\n",
    "     Parameters:\n",
    "    - y (array-like): Input array to smooth.\n",
    "    - window_length (int): Length of the filter window; must be odd and <= len(y).\n",
    "    - polyorder (int): Order of the polynomial used to fit the samples.\n",
    "    - mode (str): Padding mode for edges (default: 'interp').\n",
    "    - cval (float): Value to use for padding when mode='constant' (default: 0.0).\n",
    "    \"\"\"\n",
    "    # Validate and adjust window_length\n",
    "    if window_length > len(y):\n",
    "        window_length = len(y) if len(y) % 2 == 1 else len(y) - 1\n",
    "        window_length = max(3, window_length)\n",
    "    # Validate polyorder\n",
    "    if polyorder >= window_length:\n",
    "        raise ValueError(f\"polyorder ({polyorder}) must be less than window_length ({window_length}).\")\n",
    "    # Apply smoothing\n",
    "    return savgol_filter(y, window_length=window_length, polyorder=polyorder, mode=mode, cval=cval)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aba1e6e1",
   "metadata": {},
   "source": [
    "## 3.7 Visualise performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b652f09",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_horizontal_panel(\n",
    "    fpr_list, tpr_list, roc_auc_list, precision_list, recall_list, pr_auc_list,\n",
    "    calibration_prob_true_list, calibration_prob_pred_list,\n",
    "    mean_prob_true, mean_prob_pred, brier_scores,\n",
    "    colormap='cool_r', save_path=None, smoothing_function=None):\n",
    "    \"\"\"\n",
    "    Plot ROC, PR, and Calibration curves in a horizontal layout.\n",
    "    \"\"\"\n",
    "    # Ensure colormap is valid\n",
    "    cmap = plt.get_cmap(colormap)\n",
    "    colors = cmap(np.linspace(0, 1, len(fpr_list)))\n",
    "    # Create the subplots (1 row, 3 columns)\n",
    "    fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 5))\n",
    "    # ROC Curve (Panel A)\n",
    "    for i, (fpr, tpr, color) in enumerate(zip(fpr_list, tpr_list, colors)):\n",
    "        ax1.plot(fpr, tpr, lw=1.5, color=color, label=f'Fold {i + 1} (AUROC={roc_auc_list[i]:.3f})')\n",
    "    mean_fpr = np.linspace(0, 1, 100)\n",
    "    mean_tpr = np.mean([np.interp(mean_fpr, fpr, tpr) for fpr, tpr in zip(fpr_list, tpr_list)], axis=0)\n",
    "    ax1.plot(mean_fpr, mean_tpr, color='black', lw=2, label='Average AUROC')\n",
    "    ax1.plot([0, 1], [0, 1], 'k--', label='No Discrimination')\n",
    "    ax1.set_xlabel(\"False Positive Rate\")\n",
    "    ax1.set_ylabel(\"True Positive Rate\")\n",
    "    ax1.text(-0.1, 1.05, \"A\", transform=ax1.transAxes, fontsize=14, fontweight='bold', va='top')\n",
    "    ax1.legend(loc=\"lower right\", framealpha=0.2)\n",
    "    ax1.set_xlim(-0.05, 1.05)\n",
    "    ax1.set_ylim(-0.05, 1.05)\n",
    "    # PR Curve (Panel B)\n",
    "    for i, (precision, recall, color) in enumerate(zip(precision_list, recall_list, colors)):\n",
    "        ax2.plot(recall, precision, lw=1.5, color=color, label=f'Fold {i + 1} (AUPRC={pr_auc_list[i]:.3f})')\n",
    "    mean_recall = np.linspace(0, 1, 100)\n",
    "    mean_precision = np.mean([np.interp(mean_recall, recall[::-1], precision[::-1]) for recall, precision in zip(recall_list, precision_list)], axis=0)\n",
    "    ax2.plot(mean_recall, mean_precision, color='black', lw=2, label='Average AUPRC')\n",
    "    ax2.set_xlabel(\"Recall (Sensitivity)\")\n",
    "    ax2.set_ylabel(\"Precision\")\n",
    "    ax2.text(-0.1, 1.05, \"B\", transform=ax2.transAxes, fontsize=14, fontweight='bold', va='top')\n",
    "    ax2.legend(loc=\"lower right\", framealpha=0.2)\n",
    "    ax2.set_xlim(-0.05, 1.05)\n",
    "    ax2.set_ylim(-0.05, 1.05)\n",
    "    # Calibration Curve (Panel C)\n",
    "    for i, (prob_true, prob_pred, color) in enumerate(zip(calibration_prob_true_list, calibration_prob_pred_list, colors)):\n",
    "        idx_sorted = np.argsort(prob_pred)\n",
    "        prob_pred_sorted = np.array(prob_pred)[idx_sorted]\n",
    "        prob_true_sorted = np.array(prob_true)[idx_sorted]\n",
    "        # Smoothed calibration curves\n",
    "        if smoothing_function is not None:\n",
    "            prob_true_smoothed = smoothing_function(prob_true_sorted)\n",
    "            ax3.plot(prob_pred_sorted, prob_true_smoothed, lw=1.5, color=color, label=f'Fold {i + 1} (Brier={brier_scores[i]:.4f})')\n",
    "        else:\n",
    "            ax3.plot(prob_pred_sorted, prob_true_sorted, lw=1.5, color=color, label=f'Fold {i + 1} (Brier={brier_scores[i]:.4f})')\n",
    "        # Scatter plot for raw data\n",
    "        ax3.scatter(prob_pred_sorted, prob_true_sorted, color=color, s=10, alpha=0.5, label='_nolegend_')\n",
    "    # Average Calibration Curve\n",
    "    if smoothing_function is not None:\n",
    "        mean_prob_true_smoothed = smoothing_function(mean_prob_true)\n",
    "        ax3.plot(mean_prob_pred, mean_prob_true_smoothed, color='black', lw=2, label='Average Calibration')\n",
    "    else:\n",
    "        ax3.plot(mean_prob_pred, mean_prob_true, color='black', lw=2, label='Average Calibration')\n",
    "    ax3.plot([0, 1], [0, 1], 'k--', label='Perfect Calibration')\n",
    "    ax3.set_xlabel(\"Predicted Probability\")\n",
    "    ax3.set_ylabel(\"Observed Probability\")\n",
    "    ax3.text(-0.1, 1.05, \"C\", transform=ax3.transAxes, fontsize=14, fontweight='bold', va='top')\n",
    "    ax3.legend(loc=\"upper left\", framealpha=0.2)\n",
    "    ax3.set_xlim(-0.05, 1.05)\n",
    "    ax3.set_ylim(-0.05, 1.05)\n",
    "    # Finalize and show/save\n",
    "    plt.tight_layout()\n",
    "    if save_path:\n",
    "        print(f\"Saving figure to {save_path}\")\n",
    "        plt.savefig(save_path, dpi=300, bbox_inches=\"tight\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9db680e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_horizontal_panel(\n",
    "    fpr_list, tpr_list, roc_auc_list, precision_list, recall_list, pr_auc_list,\n",
    "    calibration_prob_true_list, calibration_prob_pred_list,\n",
    "    mean_prob_true, mean_prob_pred, brier_scores, colormap='cool_r', save_path='python_output/SF3_logistic_perform_calib.png',\n",
    "    smoothing_function=lambda y: smooth_curve_savgol(y, window_length=5, polyorder=1)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "588d1297",
   "metadata": {},
   "source": [
    "# 4.0 Elastic Net"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6e7d42d",
   "metadata": {},
   "source": [
    "## 4.1 Define threshold functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab8085b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#1) Quadratic cost function\n",
    "def quadratic_cost(y_true, y_proba, penalty_fp=1, penalty_fn=15, penalty_tp=0.1):\n",
    "    costs = []\n",
    "    thresholds = np.linspace(0, 1, 100)\n",
    "    for t in thresholds:\n",
    "        y_pred = (y_proba >= t).astype(int)\n",
    "        tn, fp, fn, tp = confusion_matrix(y_true, y_pred).ravel()\n",
    "        fp_cost = penalty_fp * np.sum((1 - y_proba[(y_true == 0) & (y_pred == 1)]) ** 2)\n",
    "        fn_cost = penalty_fn * np.sum((y_proba[(y_true == 1) & (y_pred == 0)]) ** 2)\n",
    "        tp_cost = penalty_tp * np.sum((1 - y_proba[(y_true == 1) & (y_pred == 1)]) ** 2)  # Add TP penalty\n",
    "        cost = fp_cost + fn_cost + tp_cost\n",
    "        costs.append(cost)\n",
    "    optimal_threshold = thresholds[np.argmin(costs)]\n",
    "    return optimal_threshold, costs\n",
    "# 6) Utility maximization function\n",
    "def utility_maximization(y_true, y_proba, utility_tp=5, utility_tn=10, cost_fp=10, cost_fn=20):\n",
    "    utilities = []\n",
    "    thresholds = np.linspace(0, 1, 100)\n",
    "    for t in thresholds:\n",
    "        y_pred = (y_proba >= t).astype(int)\n",
    "        tn, fp, fn, tp = confusion_matrix(y_true, y_pred).ravel()\n",
    "        utility = (\n",
    "            utility_tp * tp + utility_tn * tn -  # Reward for correct predictions\n",
    "            cost_fp * fp - cost_fn * fn         # Penalty for incorrect predictions\n",
    "        )\n",
    "        utilities.append(utility)\n",
    "    optimal_threshold = thresholds[np.argmax(utilities)]\n",
    "    return optimal_threshold, utilities"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e95327e1",
   "metadata": {},
   "source": [
    "## 4.2 Nested cross validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "655585d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 0) Lists for storing final metrics across outer folds\n",
    "fpr_list, tpr_list, roc_auc_list = [], [], []\n",
    "precision_list, recall_list, pr_auc_list = [], [], []\n",
    "calibration_prob_true_list, calibration_prob_pred_list = [], []\n",
    "interpolated_prob_true_list, interpolated_prob_pred_list = [], []\n",
    "overconfidence_penalties, brier_scores, brier_optimized_metrics = [], [], []\n",
    "predictions_list, youden_metrics, default_metrics = [], [], []\n",
    "quadratic_costs, quadratic_metrics_list = [], []\n",
    "utilities, utility_metrics_list = [], []\n",
    "\n",
    "# Define a fixed grid for predicted probabilities\n",
    "fixed_grid = np.linspace(0, 1, 100)\n",
    "# 1) Prepare data\n",
    "X = df2.drop('fail', axis=1)\n",
    "y = df2['fail']\n",
    "# 2) Define outer CV\n",
    "outer_cv = KFold(n_splits=5, shuffle=True, random_state=seed)  # Fixed seed\n",
    "# 3) Outer CV loop\n",
    "for fold, (train_ix, test_ix) in enumerate(outer_cv.split(X)):\n",
    "    X_train, X_test = X.iloc[train_ix], X.iloc[test_ix]\n",
    "    y_train, y_test = y.iloc[train_ix], y.iloc[test_ix]\n",
    "    # 4) SMOTE on the train fold\n",
    "    smote = BorderlineSMOTE(random_state=seed)\n",
    "    X_train_res, y_train_res = smote.fit_resample(X_train, y_train)\n",
    "    # 5A) Inner CV for hyperparameter tuning\n",
    "    cv_inner = KFold(n_splits=3, shuffle=True, random_state=seed)\n",
    "    logistic = LogisticRegression(random_state=seed, solver='saga', penalty='elasticnet')\n",
    "    # 5B) Define hyperparameter search space\n",
    "    space = {\n",
    "        'C': np.logspace(-4, 4, 20),\n",
    "        'max_iter': [50, 100, 250, 500, 750, 1000],\n",
    "        'class_weight': [None, 'balanced'],\n",
    "        'fit_intercept': [True, False],\n",
    "        'l1_ratio': np.linspace(0.01, 0.1, 10),  # Only applicable if penalty='elasticnet'\n",
    "        'tol': np.logspace(-6, -1, 6),\n",
    "        'warm_start': [False]\n",
    "    }\n",
    "    # 5C) Randomized search\n",
    "    search = RandomizedSearchCV(\n",
    "        logistic, space, scoring='f1_weighted', cv=cv_inner,refit=True, random_state=seed, n_iter=500\n",
    "    )\n",
    "    # 6A) Fit on inner train fold\n",
    "    result = search.fit(X_train_res, y_train_res)\n",
    "    # 6B) Get best model\n",
    "    best_model = result.best_estimator_\n",
    "    # 7) Predictions on outer test fold\n",
    "    y_proba = best_model.predict_proba(X_test)[:, 1]\n",
    "    thresholds = np.linspace(0, 1, 100)\n",
    "    # i) ROC & PR\n",
    "    fpr, tpr, thresholds = roc_curve(y_test, y_proba)\n",
    "    precision, recall, thresholds_pr = precision_recall_curve(y_test, y_proba)\n",
    "    roc_auc_list.append(auc(fpr, tpr))\n",
    "    pr_auc_list.append(average_precision_score(y_test, y_proba))\n",
    "    fpr_list.append(fpr)\n",
    "    tpr_list.append(tpr)\n",
    "    precision_list.append(precision)\n",
    "    recall_list.append(recall)\n",
    "    # Performance metrics at 0.5 threshold\n",
    "    y_pred_0_5 = (y_proba >= 0.5).astype(int)\n",
    "    default_metrics.append(calculate_metrics(y_test, y_pred_0_5, y_proba))\n",
    "    # ii) Calibration\n",
    "    prob_true, prob_pred = calibration_curve(y_test, y_proba, n_bins=11, strategy='uniform')\n",
    "    calibration_prob_true_list.append(prob_true)\n",
    "    calibration_prob_pred_list.append(prob_pred)\n",
    "    fold_brier_score = brier_score_loss(y_test, y_proba)\n",
    "    brier_scores.append(fold_brier_score)\n",
    "    # Interpolation for consistent grid\n",
    "    interp_prob_true = np.interp(fixed_grid, prob_pred, prob_true, left=0, right=1)\n",
    "    interpolated_prob_true_list.append(interp_prob_true)\n",
    "    interpolated_prob_pred_list.append(fixed_grid)\n",
    "    # iii) Overconfidence\n",
    "    overconfidence_penalties.append(calculate_overconfidence_penalty(y_test, y_proba))\n",
    "    # iv) Youden threshold\n",
    "    youden_index = tpr - fpr\n",
    "    max_youden_idx = np.argmax(youden_index)\n",
    "    optimal_youden_threshold = thresholds[max_youden_idx]\n",
    "    yhat_youden = (y_proba >= optimal_youden_threshold).astype(int)\n",
    "    youden_metrics.append(calculate_metrics(y_test, yhat_youden, y_proba))\n",
    "    # v) Brier-optimized threshold\n",
    "    brier_threshold_scores = [\n",
    "        brier_score_loss(y_test, (y_proba >= thr).astype(int))\n",
    "        for thr in thresholds\n",
    "    ]\n",
    "    best_brier_idx = np.argmin(brier_threshold_scores)\n",
    "    optimal_brier_threshold = thresholds[best_brier_idx]\n",
    "    yhat_brier_optimized = (y_proba >= optimal_brier_threshold).astype(int)\n",
    "    brier_optimized_metrics.append(calculate_metrics(y_test, yhat_brier_optimized, y_proba))\n",
    "    # a) Quadratic Cost Function\n",
    "    optimal_quadratic_threshold, quadratic_cost_values = quadratic_cost(y_test, y_proba)\n",
    "    quadratic_costs.append(quadratic_cost_values)\n",
    "    yhat_quadratic = (y_proba >= optimal_quadratic_threshold).astype(int)\n",
    "    quadratic_metrics = calculate_metrics(y_test, yhat_quadratic, y_proba)\n",
    "    quadratic_metrics_list.append(quadratic_metrics)\n",
    "    # c) Utility Maximization Function\n",
    "    optimal_utility_threshold, utility_values = utility_maximization(y_test, y_proba)\n",
    "    utilities.append(utility_values)\n",
    "    yhat_utility = (y_proba >= optimal_utility_threshold).astype(int)\n",
    "    utility_metrics = calculate_metrics(y_test, yhat_utility, y_proba)\n",
    "    utility_metrics_list.append(utility_metrics)\n",
    "    # Store predictions\n",
    "    temp_df = pd.DataFrame({\n",
    "        'fold': fold + 1,\n",
    "        'sid': X_test.index,  \n",
    "        'y_true': y_test.values,\n",
    "        'y_proba': y_proba,\n",
    "        'y_pred': y_pred_0_5\n",
    "    })\n",
    "    predictions_list.append(temp_df)\n",
    "predictions_df = pd.concat(predictions_list, ignore_index=True)\n",
    "# Aggregate calibration data\n",
    "mean_brier, std_brier, mean_prob_true,  mean_prob_pred = aggregate_calibration_metrics(\n",
    "    brier_scores, calibration_prob_true_list, calibration_prob_pred_list\n",
    ")\n",
    "# Aggregate threshold-based metrics\n",
    "youden_summary, youden_std = aggregate_metrics(youden_metrics)\n",
    "default_summary, default_std = aggregate_metrics(default_metrics)\n",
    "brier_optimized_summary, brier_optimized_std = aggregate_metrics(brier_optimized_metrics)\n",
    "# Aggregate metrics for the new cost functions\n",
    "quadratic_summary, quadratic_std = aggregate_metrics(quadratic_metrics_list)\n",
    "utility_summary, utility_std = aggregate_metrics(utility_metrics_list)\n",
    "# Print results\n",
    "print(f\"AUROC: {np.mean(roc_auc_list):.4f} ± {np.std(roc_auc_list):.4f}\")\n",
    "print(f\"AUPRC: {np.mean(pr_auc_list):.4f} ± {np.std(pr_auc_list):.4f}\")\n",
    "print(f\"Brier Score: {mean_brier:.4f} ± {std_brier:.4f}\")\n",
    "print(f\"Overconfidence Penalty: {np.mean(overconfidence_penalties):.4f} ± {np.std(overconfidence_penalties):.4f}\")\n",
    "print(\"\\nMetrics at 0.5 Threshold:\")\n",
    "for m, val in default_summary.items():\n",
    "    print(f\"{m.capitalize()}: {val:.4f} ± {default_std[m]:.4f}\")\n",
    "print(\"\\nYouden Index Threshold Metrics:\")\n",
    "for m, val in youden_summary.items():\n",
    "    print(f\"{m.capitalize()}: {val:.4f} ± {youden_std[m]:.4f}\")\n",
    "print(\"\\nBrier Optimized Threshold Metrics:\")\n",
    "for m, val in brier_optimized_summary.items():\n",
    "    print(f\"{m.capitalize()}: {val:.4f} ± {brier_optimized_std[m]:.4f}\")\n",
    "# Print the new metrics\n",
    "print(\"\\nQuadratic Cost Threshold Metrics:\")\n",
    "for m, val in quadratic_summary.items():\n",
    "    print(f\"{m.capitalize()}: {val:.4f} ± {quadratic_std[m]:.4f}\")\n",
    "print(\"\\nUtility Maximization Threshold Metrics:\")\n",
    "for m, val in utility_summary.items():\n",
    "    print(f\"{m.capitalize()}: {val:.4f} ± {utility_std[m]:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ebd1b33",
   "metadata": {},
   "source": [
    "## 4.3 Overconfident cases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da82676d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define your filter: y_proba > 0.8 or < 0.1, and wrong prediction\n",
    "mask = (\n",
    "    ((predictions_df['y_proba'] > 0.85) | (predictions_df['y_proba'] < 0.15)) &\n",
    "    (predictions_df['y_pred'] != predictions_df['y_true'])\n",
    ")\n",
    "bad_cases = predictions_df.loc[mask].copy()\n",
    "# Sort by fold, etc.\n",
    "bad_cases.sort_values(by=['fold', 'sid'], inplace=True)\n",
    "# Merge to get all the data the index in df, reset it to a column\n",
    "df_reset = df.reset_index()\n",
    "# Perform an inner merge to retain only matched rows\n",
    "merged_df = bad_cases.merge(df_reset, how='inner', on='sid')\n",
    "# If you want to set 'sid' back as the index after the merge\n",
    "merged_df = merged_df.set_index('sid')\n",
    "# Export\n",
    "merged_df.to_csv('python_output/TX1_overconfident_cases.csv', index=False)\n",
    "merged_df.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49aab75a",
   "metadata": {},
   "source": [
    "## 4.4 Export metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "219ec943",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_results_table():\n",
    "    # Dynamically fetch the metrics keys from one of the summary dictionaries\n",
    "    metrics = list(youden_summary.keys())\n",
    "    methods = {\n",
    "        \"Logistic\": (log_metrics_summary, log_metrics_std),\n",
    "        \"Default\": (default_summary, default_std),\n",
    "        \"Youden\": (youden_summary, youden_std),\n",
    "        \"Quadratic\": (quadratic_summary, quadratic_std),\n",
    "        \"Utility\": (utility_summary, utility_std),\n",
    "        \"Brier\": (brier_optimized_summary, brier_optimized_std),\n",
    "    }\n",
    "    # Initialize the DataFrame\n",
    "    metrics_data = {\"Metric\": metrics}\n",
    "    # Add columns dynamically\n",
    "    for method, (mean_data, sd_data) in methods.items():\n",
    "        metrics_data[f\"{method} (Mean)\"] = [mean_data[metric] for metric in metrics]\n",
    "        metrics_data[f\"{method} (SD)\"] = [sd_data[metric] for metric in metrics]\n",
    "    # Create DataFrame\n",
    "    results_table = pd.DataFrame(metrics_data)\n",
    "    return results_table\n",
    "\n",
    "# Generate the results table\n",
    "results_table = create_results_table()\n",
    "# Format with 2 decimal places\n",
    "results_table = results_table.applymap(lambda x: f\"{x:.2f}\" if isinstance(x, float) else x)\n",
    "# Export the table to a CSV file\n",
    "results_table.to_csv(\"python_output/MT2_perform_metrics.csv\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1337ee7",
   "metadata": {},
   "source": [
    "## 4.5 Best model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a57d27f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06b72561",
   "metadata": {},
   "source": [
    "## 4.6 Visualise performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebd63429",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_horizontal_panel(\n",
    "    fpr_list, tpr_list, roc_auc_list, precision_list, recall_list, pr_auc_list,\n",
    "    calibration_prob_true_list, calibration_prob_pred_list,\n",
    "    mean_prob_true, mean_prob_pred, brier_scores, colormap='cool_r', save_path='python_output/F1_perform_calib.png',\n",
    "    smoothing_function=lambda y: smooth_curve_savgol(y, window_length=5, polyorder=1)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9f630fa",
   "metadata": {},
   "source": [
    "## 4.7 Visualise feature importance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f46ee7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ensure the output directory exists\n",
    "output_dir = \"python_output\"\n",
    "os.makedirs(output_dir, exist_ok=True)\n",
    "# Declare feature vector and target variable\n",
    "X = df2.drop(['fail'], axis=1)\n",
    "y = df2['fail']\n",
    "best_model.fit(X, y)\n",
    "# Rename columns for better readability\n",
    "rename_dict = {\n",
    "    'log_rna': 'HCV RNA',\n",
    "    'ast_alt_ratio': 'AST-ALT ratio',\n",
    "    'GT_3': 'Genotype 3',\n",
    "    'duration': 'DAA duration',\n",
    "    'daa_class_NS5A+NS5B': 'NS5A+NS5B inhibitor',\n",
    "    'daa_class_PI+NS5A': 'PI+NS5A inhibitor',\n",
    "    'daa_class_PI+NS5A+NS5B': 'PI+NS5A+NS5B inhibitor',\n",
    "}\n",
    "X_renamed = X.rename(columns=rename_dict)\n",
    "# Create SHAP explainer and calculate SHAP values\n",
    "explainer = shap.Explainer(best_model, X_renamed)\n",
    "shap_values = explainer(X_renamed)\n",
    "# Create a Matplotlib figure explicitly\n",
    "fig = plt.figure(figsize=(10, 6))  # Adjust the size as needed\n",
    "# Plot the SHAP summary plot\n",
    "shap.summary_plot(\n",
    "    shap_values.values,\n",
    "    X_renamed,\n",
    "    cmap=colormaps['cool'],\n",
    "    show=True  # Prevent the plot from being displayed immediately\n",
    ")\n",
    "# Save the figure in high resolution\n",
    "plt.close(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba75ceeb",
   "metadata": {},
   "source": [
    "## 4.8 Visualise HCV RNA cutoffs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52900fea",
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_names = X.columns.tolist()  # Use the DataFrame's column names\n",
    "# Fixed values for other features\n",
    "medi_ast_alt_ratio = df['ast_alt_ratio'].median()\n",
    "fixed_values = {\n",
    "    'ast_alt_ratio': medi_ast_alt_ratio,\n",
    "    'daa_class_PI+NS5A': 0,\n",
    "    'daa_class_NS5A+NS5B': 0,\n",
    "    'duration': 28,\n",
    "    'GT_3': 0\n",
    "}\n",
    "# Generate a range of log_rna values\n",
    "log_rna_values = np.linspace(0, 9, num=1000)\n",
    "# Initialize data storage for probabilities\n",
    "probabilities_by_group = {}\n",
    "cutoff_statistics = []\n",
    "# Populate probabilities_by_group\n",
    "for duration in [28, 42]:\n",
    "    for gt3 in [0, 1]:\n",
    "        group_key = (duration, gt3)\n",
    "        probabilities_by_group[group_key] = {}\n",
    "        for daa_class in [\"NS5A+NS5B\", \"PI+NS5A\"]:\n",
    "            # Update fixed values\n",
    "            fixed_values.update({\n",
    "                'duration': duration,\n",
    "                'GT_3': gt3,\n",
    "                'daa_class_NS5A+NS5B': 1 if daa_class == \"NS5A+NS5B\" else 0,\n",
    "                'daa_class_PI+NS5A': 1 if daa_class == \"PI+NS5A\" else 0\n",
    "            })\n",
    "            # Create DataFrame for predictions\n",
    "            dfx = pd.DataFrame(fixed_values, index=range(len(log_rna_values)))\n",
    "            dfx['log_rna'] = log_rna_values\n",
    "            dfx = dfx[feature_names]\n",
    "            # Predict probabilities\n",
    "            probabilities = best_model.predict_proba(dfx)[:, 1]\n",
    "            probabilities_by_group[group_key][daa_class] = probabilities\n",
    "# Simulate bootstrapping\n",
    "bootstrapped_probs_by_group = defaultdict(list)\n",
    "for (duration, gt3), daa_classes in probabilities_by_group.items():\n",
    "    for daa_class, probs in daa_classes.items():\n",
    "        for _ in range(100):  # Simulate 1000 bootstraps\n",
    "            bootstrapped_probs = probs + np.random.normal(0, 0.02, size=probs.shape)\n",
    "            bootstrapped_probs_by_group[(duration, gt3, daa_class)].append(bootstrapped_probs)\n",
    "# Plotting and statistics calculation\n",
    "fig, axes = plt.subplots(2, 2, figsize=(10, 7), sharey=False)\n",
    "plot_positions = [(28, 1), (28, 0), (42, 1), (42, 0)]\n",
    "subplot_labels = [\"A\", \"B\", \"C\", \"D\"]\n",
    "subplot_descriptions = [\n",
    "    \"GT3 28-days duration\", \n",
    "    \"GT3 42-days duration\", \n",
    "    \"non-GT3 28-days duration\", \n",
    "    \"non-GT3 42-days duration\"\n",
    "]\n",
    "# Define unique colors for each subplot\n",
    "subplot_colors = {\n",
    "    \"A\": [plt.get_cmap('cool')(0.30), plt.get_cmap('cool')(0.65)],\n",
    "    \"B\": [plt.get_cmap('cool')(0.25), plt.get_cmap('cool')(0.60)],\n",
    "    \"C\": [plt.get_cmap('cool')(0.20), plt.get_cmap('cool')(0.55)],\n",
    "    \"D\": [plt.get_cmap('cool')(0.15), plt.get_cmap('cool')(0.50)],\n",
    "}\n",
    "cutoff_statistics = []  # Reset statistics collection\n",
    "for (duration, gt3), ax, label, description in zip(plot_positions, axes.flatten(), subplot_labels, subplot_descriptions):\n",
    "    # Add subplot label in the top-left\n",
    "    ax.text(-0.07, 1.10, label, transform=ax.transAxes, fontsize=12, fontweight='bold', va='top')\n",
    "    # Add description in the bottom-left\n",
    "    ax.text(\n",
    "        0.575, 0.02, description, \n",
    "        transform=ax.transAxes, \n",
    "        fontsize=10, \n",
    "        fontweight='semibold', \n",
    "        va='bottom', \n",
    "        ha='left',\n",
    "        color='black'\n",
    "    )\n",
    "    colors = subplot_colors[label]  # Get unique colors for this subplot\n",
    "    for i, (daa_class, probs) in enumerate(probabilities_by_group[(duration, gt3)].items()):\n",
    "        mean_probs = np.mean(bootstrapped_probs_by_group[(duration, gt3, daa_class)], axis=0)\n",
    "        lower_ci = np.percentile(bootstrapped_probs_by_group[(duration, gt3, daa_class)], 2.5, axis=0)\n",
    "        upper_ci = np.percentile(bootstrapped_probs_by_group[(duration, gt3, daa_class)], 97.5, axis=0)\n",
    "        # Plot mean curve and CI ribbon\n",
    "        ax.plot(log_rna_values, mean_probs, color=colors[i], lw=1.5, label=daa_class)\n",
    "        ax.fill_between(log_rna_values, lower_ci, upper_ci, color=colors[i], alpha=0.6)\n",
    "        # Find the cutoff at threshold 0.5\n",
    "        crossing_mean = np.where(mean_probs >= 0.5)[0]\n",
    "        crossing_lower = np.where(lower_ci >= 0.5)[0]\n",
    "        crossing_upper = np.where(upper_ci >= 0.5)[0]\n",
    "        mean_cutoff = log_rna_values[crossing_mean[0]] if len(crossing_mean) > 0 else np.nan\n",
    "        lower_cutoff = log_rna_values[crossing_lower[0]] if len(crossing_lower) > 0 else np.nan\n",
    "        upper_cutoff = log_rna_values[crossing_upper[0]] if len(crossing_upper) > 0 else np.nan\n",
    "        # Store statistics\n",
    "        cutoff_statistics.append({\n",
    "            \"Duration (days)\": duration,\n",
    "            \"GT3\": gt3,\n",
    "            \"DAA Class\": daa_class,\n",
    "            \"Mean Cutoff\": mean_cutoff,\n",
    "            \"95% CI Lower\": upper_cutoff,\n",
    "            \"95% CI Upper\": lower_cutoff\n",
    "        })\n",
    "    ax.axhline(0.5, color='black', linestyle='--', label='Threshold = 0.5')\n",
    "    ax.set_xlim([0.3, 8.7])\n",
    "    ax.set_ylim([0.0, 1.0])\n",
    "    ax.grid(True, alpha=0.5)\n",
    "    ax.legend(loc='upper left', framealpha=0.8)\n",
    "    # Add individual x and y labels\n",
    "    ax.set_xlabel(\"HCV RNA (Log10 IU/mL)\", fontsize=10)\n",
    "    ax.set_ylabel(\"Probability of treatment failure\", fontsize=10)\n",
    "plt.tight_layout()\n",
    "plt.savefig('python_output/F3_HCV_RNA_cutoffs.png', dpi=300)\n",
    "plt.show()\n",
    "# Export cutoff statistics\n",
    "cutoff_df = pd.DataFrame(cutoff_statistics)\n",
    "cutoff_df.to_csv(\"python_output/TX2_HCV_RNA_cutoffs.csv\", index=False)\n",
    "print(\"\\nCutoff Statistics:\")\n",
    "print(cutoff_df)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}