In this post I provide a quickstart Python code based on the one provided in https://challengedata.ens.fr/challenges/162. In particular, I use Python package survivalist, that does Survival analysis with Machine Learning and uncertainty quantification (introduced in https://thierrymoudiki.github.io/blog/2024/12/15/python/agnostic-survival-analysis).

Data Challenge : Leukemia Risk Prediction

GOAL OF THE CHALLENGE and WHY IT IS IMPORTANT:

The goal of the challenge is to predict disease risk for patients with blood cancer, in the context of specific subtypes of adult myeloid leukemias.

The risk is measured through the overall survival of patients, i.e. the duration of survival from the diagnosis of the blood cancer to the time of death or last follow-up.

Estimating the prognosis of patients is critical for an optimal clinical management.
For exemple, patients with low risk-disease will be offered supportive care to improve blood counts and quality of life, while patients with high-risk disease will be considered for hematopoietic stem cell transplantion.

The performance metric used in the challenge is the IPCW-C-Index.

THE DATASETS

The training set is made of 3,323 patients.

The test set is made of 1,193 patients.

For each patient, you have acces to CLINICAL data and MOLECULAR data.

The details of the data are as follows:

• OUTCOME:
  • OS_YEARS = Overall survival time in years
  • OS_STATUS = 1 (death) , 0 (alive at the last follow-up)

• CLINICAL DATA, with one line per patient:

  • ID = unique identifier per patient
  • CENTER = clinical center
  • BM_BLAST = Bone marrow blasts in % (blasts are abnormal blood cells)
  • WBC = White Blood Cell count in Giga/L
  • ANC = Absolute Neutrophil count in Giga/L
  • MONOCYTES = Monocyte count in Giga/L
  • HB = Hemoglobin in g/dL
  • PLT = Platelets coutn in Giga/L
  • CYTOGENETICS = A description of the karyotype observed in the blood cells of the patients, measured by a cytogeneticist. Cytogenetics is the science of chromosomes. A karyotype is performed from the blood tumoral cells. The convention for notation is ISCN (https://en.wikipedia.org/wiki/International_System_for_Human_Cytogenomic_Nomenclature). Cytogenetic notation are: https://en.wikipedia.org/wiki/Cytogenetic_notation. Note that a karyotype can be normal or abnornal. The notation 46,XX denotes a normal karyotype in females (23 pairs of chromosomes including 2 chromosomes X) and 46,XY in males (23 pairs of chromosomes inclusing 1 chromosme X and 1 chromsome Y). A common abnormality in the blood cancerous cells might be for exemple a loss of chromosome 7 (monosomy 7, or -7), which is typically asssociated with higher risk disease

• GENE MOLECULAR DATA, with one line per patient per somatic mutation. Mutations are detected from the sequencing of the blood tumoral cells.
  We call somatic (= acquired) mutations the mutations that are found in the tumoral cells but not in other cells of the body.

  • ID = unique identifier per patient
  • CHR START END = position of the mutation on the human genome
  • REF ALT = reference and alternate (=mutant) nucleotide
  • GENE = the affected gene
  • PROTEIN_CHANGE = the consequence of the mutation on the protei that is expressed by a given gene
  • EFFECT = a broad categorization of the mutation consequences on a given gene.
  • VAF = Variant Allele Fraction = it represents the proportion of cells with the deleterious mutations.

!pip install survivalist --upgrade --no-cache-dir --verbose

!pip install scikit-survival

!pip install nnetsauce

!curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh -s -- -y

import os
os.environ['PATH'] = f"/root/.cargo/bin:{os.environ['PATH']}"

!echo $PATH
!rustc --version
!cargo --version

!pip install glmnetforpython --verbose --upgrade --no-cache-dir

!pip install git+https://github.com/Techtonique/genbooster.git --upgrade --no-cache-dir

import numpy as np
import pandas as pd
import glmnetforpython as glmnet
from matplotlib import pyplot as plt
from sklearn.utils.discovery import all_estimators
from sklearn.datasets import load_iris, load_breast_cancer, load_wine, load_digits
from sklearn.linear_model import Ridge, RidgeCV
from sklearn.tree import ExtraTreeRegressor
from sklearn.model_selection import train_test_split
from genbooster.genboosterregressor import BoosterRegressor
from genbooster.randombagregressor import RandomBagRegressor
from sklearn.metrics import mean_squared_error
from tqdm import tqdm
from sklearn.utils.discovery import all_estimators
from time import time

from survivalist.nonparametric import kaplan_meier_estimator
from survivalist.datasets import load_whas500, load_gbsg2, load_veterans_lung_cancer
from survivalist.custom import SurvivalCustom
from survivalist.custom import PISurvivalCustom
from survivalist.ensemble import GradientBoostingSurvivalAnalysis

import matplotlib.pyplot as plt
from sklearn.tree import plot_tree
from sklearn.model_selection import train_test_split
from sksurv.ensemble import RandomSurvivalForest
from sksurv.linear_model import CoxPHSurvivalAnalysis
from sksurv.metrics import concordance_index_censored , concordance_index_ipcw
from sklearn.impute import SimpleImputer
from sksurv.util import Surv

# Clinical Data
df = pd.read_csv("./clinical_train.csv")
df_eval = pd.read_csv("./clinical_test.csv")

# Molecular Data
maf_df = pd.read_csv("./molecular_train.csv")
maf_eval = pd.read_csv("./molecular_test.csv")

target_df = pd.read_csv("./target_train.csv")
#target_df_test = pd.read_csv("./target_test.csv")

# Preview the data
df.head()

Step 1: Data Preparation (clinical data only)

For survival analysis, we’ll format the dataset so that OS_YEARS represents the time variable and OS_STATUS represents the event indicator.

# Drop rows where 'OS_YEARS' is NaN if conversion caused any issues
target_df.dropna(subset=['OS_YEARS', 'OS_STATUS'], inplace=True)

# Check the data types to ensure 'OS_STATUS' is boolean and 'OS_YEARS' is numeric
print(target_df[['OS_STATUS', 'OS_YEARS']].dtypes)

# Contarget_dfvert 'OS_YEARS' to numeric if it isn’t already
target_df['OS_YEARS'] = pd.to_numeric(target_df['OS_YEARS'], errors='coerce')

# Ensure 'OS_STATUS' is boolean
target_df['OS_STATUS'] = target_df['OS_STATUS'].astype(bool)

# Select features
features = ['BM_BLAST', 'HB', 'PLT']
target = ['OS_YEARS', 'OS_STATUS']

# Create the survival data format
X = df.loc[df['ID'].isin(target_df['ID']), features]
y = Surv.from_dataframe('OS_STATUS', 'OS_YEARS', target_df)

OS_STATUS    float64
OS_YEARS     float64
dtype: object

Step 2: Splitting the Dataset

We’ll split the data into training and testing sets to evaluate the model’s performance.

# Split the data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# Survival-aware imputation for missing values
imputer = SimpleImputer(strategy="median")
X_train[['BM_BLAST', 'HB', 'PLT']] = imputer.fit_transform(X_train[['BM_BLAST', 'HB', 'PLT']])
X_test[['BM_BLAST', 'HB', 'PLT']] = imputer.transform(X_test[['BM_BLAST', 'HB', 'PLT']])

display(X_train.head())
display(y_train)

array([(False, 1.91780822), ( True, 1.28219178), ( True, 1.49041096), ...,
       (False, 8.63561644), (False, 0.47671233), (False, 1.29041096)],
      dtype=[('OS_STATUS', '?'), ('OS_YEARS', '<f8')])

Step 3: Cox Proportional Hazards Model

\[h(t | X) = h_0(t) \exp(\beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p)\]

where \(h_0(t)\) is the baseline hazard function, and \(\beta\) values are coefficients for each covariate, representing the effect of $X$ on the hazard. Importantly, the proportional hazards assumption implies that the hazard ratios between individuals are constant over time. This approach effectively leverages both observed and censored survival times, making it a more suitable method for survival data compared to standard regression techniques that ignore censoring.

# Initialize and train the Cox Proportional Hazards model
cox = CoxPHSurvivalAnalysis()
cox.fit(X_train, y_train)

# Evaluate the model using Concordance Index IPCW
cox_cindex_train = concordance_index_ipcw(y_train, y_train, cox.predict(X_train), tau=7)[0]
cox_cindex_test = concordance_index_ipcw(y_train, y_test, cox.predict(X_test), tau=7)[0]
print(f"Cox Proportional Hazard Model Concordance Index IPCW on train: {cox_cindex_train:.5f}")
print(f"Cox Proportional Hazard Model Concordance Index IPCW on test: {cox_cindex_test:.5f}")

Cox Proportional Hazard Model Concordance Index IPCW on train: 0.66302
Cox Proportional Hazard Model Concordance Index IPCW on test: 0.66060

Step 4: other models

4 – 1 demo

import xgboost as xgb
import lightgbm as lgb

# Initialize and train the XGBoost model


event_time = [y[1] for y in y_test]
event_status = [y[0] for y in y_test]
km = kaplan_meier_estimator(event_status, event_time,
                            conf_type="log-log")
estimator = PISurvivalCustom(regr=xgb.XGBRegressor(),
                             type_pi="kde")

estimator.fit(X_train, y_train)

surv_funcs = estimator.predict_survival_function(X_test.iloc[:1])

for fn in surv_funcs.mean:
    plt.step(fn.x, fn(fn.x), where="post")
    plt.fill_between(fn.x, surv_funcs.lower[0].y, surv_funcs.upper[0].y, alpha=0.25, color="lightblue", step="post")
    plt.step(km[0], km[1], where="post", color="red", label="Kaplan-Meier")
    plt.fill_between(km[0], km[2][0], km[2][1], alpha=0.25, color="pink", step="post")
    plt.ylim(0, 1)
    plt.show()

# Evaluate the model using Concordance Index IPCW
cox_cindex_test = concordance_index_ipcw(y_train, y_test, estimator.predict(X_test.iloc[:]).mean, tau=7)[0]
print(f"Cox Proportional Hazard Model Concordance Index IPCW on test: {cox_cindex_test:.5f}")
cox_cindex_test = concordance_index_ipcw(y_train, y_test, estimator.predict(X_test.iloc[:]).lower, tau=7)[0]
print(f"Cox Proportional Hazard Model Concordance Index IPCW on test: {cox_cindex_test:.5f}")
cox_cindex_test = concordance_index_ipcw(y_train, y_test, estimator.predict(X_test.iloc[:]).upper, tau=7)[0]
print(f"Cox Proportional Hazard Model Concordance Index IPCW on test: {cox_cindex_test:.5f}")

Cox Proportional Hazard Model Concordance Index IPCW on test: 0.60130
Cox Proportional Hazard Model Concordance Index IPCW on test: 0.60106
Cox Proportional Hazard Model Concordance Index IPCW on test: 0.59588

import lightgbm as lgb

event_time = [y[1] for y in y_test]
event_status = [y[0] for y in y_test]
km = kaplan_meier_estimator(event_status, event_time,
                            conf_type="log-log")
estimator = PISurvivalCustom(regr=lgb.LGBMRegressor(verbose=0),
                             type_pi="kde")

estimator.fit(X_train, y_train)

surv_funcs = estimator.predict_survival_function(X_test.iloc[:1])

for fn in surv_funcs.mean:
    plt.step(fn.x, fn(fn.x), where="post")
    plt.fill_between(fn.x, surv_funcs.lower[0].y, surv_funcs.upper[0].y, alpha=0.25, color="lightblue", step="post")
    plt.step(km[0], km[1], where="post", color="red", label="Kaplan-Meier")
    plt.fill_between(km[0], km[2][0], km[2][1], alpha=0.25, color="pink", step="post")
    plt.ylim(0, 1)
    plt.show()

# Evaluate the model using Concordance Index IPCW
cox_cindex_test = concordance_index_ipcw(y_train, y_test, estimator.predict(X_test.iloc[:]).mean, tau=7)[0]
print(f"Cox Proportional Hazard Model Concordance Index IPCW on test: {cox_cindex_test:.5f}")
cox_cindex_test = concordance_index_ipcw(y_train, y_test, estimator.predict(X_test.iloc[:]).lower, tau=7)[0]
print(f"Cox Proportional Hazard Model Concordance Index IPCW on test: {cox_cindex_test:.5f}")
cox_cindex_test = concordance_index_ipcw(y_train, y_test, estimator.predict(X_test.iloc[:]).upper, tau=7)[0]
print(f"Cox Proportional Hazard Model Concordance Index IPCW on test: {cox_cindex_test:.5f}")

Cox Proportional Hazard Model Concordance Index IPCW on test: 0.61745
Cox Proportional Hazard Model Concordance Index IPCW on test: 0.62040
Cox Proportional Hazard Model Concordance Index IPCW on test: 0.61757

4 – 2 models galore

# prompt: loop on scikit-learn regressors
import nnetsauce as ns
import pandas as pd
import xgboost as xgb

from functools import partial
from sklearn.utils.discovery import all_estimators
from sklearn.metrics import mean_squared_error
from tqdm import tqdm
from sksurv.metrics import concordance_index_ipcw
from sksurv.linear_model import CoxPHSurvivalAnalysis
from sksurv.util import Surv
from sklearn.impute import SimpleImputer


# Get all regressors from scikit-learn
regressors = [est for est in all_estimators()  if 'Regressor' in est[0]]

# Append xgb.XGBRegressor and lgb.LGBMRegressor as (name, class) tuples
regressors += [('XGBRegressor', xgb.XGBRegressor),
 ('LGBMRegressor', partial(lgb.LGBMRegressor, verbose=0))]

results = []

for name, Regressor in tqdm(regressors):

    print("\n\n ----- base learner", name)

    try:

      # Initialize and train the model
      estimator = PISurvivalCustom(regr=Regressor(), type_pi="kde")
      estimator.fit(X_train, y_train)
      # Make predictions and evaluate the model
      y_pred = estimator.predict(X_test.iloc[:])
      c_index = concordance_index_ipcw(y_train, y_test, y_pred.mean, tau=7)[0]
      c_index_upper = concordance_index_ipcw(y_train, y_test, y_pred.upper, tau=7)[0]
      c_index_lower = concordance_index_ipcw(y_train, y_test, y_pred.lower, tau=7)[0]
      print("\n c_index", c_index)
      results.append([name, c_index, c_index_lower, c_index_upper])

      # Initialize and train the model
      estimator = PISurvivalCustom(regr=ns.CustomRegressor(Regressor()), type_pi="kde")
      estimator.fit(X_train, y_train)
      # Make predictions and evaluate the model
      y_pred = estimator.predict(X_test.iloc[:])
      c_index = concordance_index_ipcw(y_train, y_test, y_pred.mean, tau=7)[0]
      c_index_upper = concordance_index_ipcw(y_train, y_test, y_pred.upper, tau=7)[0]
      c_index_lower = concordance_index_ipcw(y_train, y_test, y_pred.lower, tau=7)[0]
      print("\n c_index", c_index)
      results.append(["custom" + name, c_index, c_index_lower, c_index_upper])

      # Initialize and train the model
      estimator = PISurvivalCustom(regr=RandomBagRegressor(Regressor()), type_pi="kde")
      estimator.fit(X_train, y_train)
      # Make predictions and evaluate the model
      y_pred = estimator.predict(X_test.iloc[:])
      c_index = concordance_index_ipcw(y_train, y_test, y_pred.mean, tau=7)[0]
      c_index_upper = concordance_index_ipcw(y_train, y_test, y_pred.upper, tau=7)[0]
      c_index_lower = concordance_index_ipcw(y_train, y_test, y_pred.lower, tau=7)[0]
      print("\n c_index", c_index)
      results.append(["bagging" + name, c_index, c_index_lower, c_index_upper])

    except Exception as e:
      continue

pd.options.display.float_format = '{:.5f}'.format
results_df = pd.DataFrame(results, columns=['Regressor', 'Concordance Index IPCW', 'lower bound', 'upper bound'])
results_df.drop(columns=['lower bound', 'upper bound'], inplace=True)
results_df.sort_values(by='Concordance Index IPCW', ascending=False, inplace=True)
results_df

results_df.shape

(60, 2)