notebook1_R.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dancing with censored data: How to survive with explainable survival analysis?\n",
    "**ML in PL Conference 2023**\n",
    "\n",
    "\n",
    "Mateusz Krzyziński, Mikołaj Spytek (**MI2.AI**)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Agenda\n",
    "\n",
    "- 9:00 - 9:15 - Introduction, technicalities\n",
    "- 9:15 - 10:00 - ML in survival analysis, creating models `(Python and R)`\n",
    "- 10:00 - 10:10 - *break*\n",
    "\n",
    "- 10:10 - 10:20 - XAI in survival analysis\n",
    "- 10:20 - 10:45 - SurvLIME `(Python and R)`\n",
    "- 10:45 - 11:30 - SurvSHAP(t) and its aggregations `(Python and R)`\n",
    "- 11:30 - 11:40 - *break*\n",
    "\n",
    "- 11:40 - 11:50 - creating explainers for any models `(R with Python via reticulate)`\n",
    "- 11:50 - 12:15 - global performance & global variable importance explanations `(R)`\n",
    "- 12:15 - 12:30 - local variable dependence explanations `(R)`\n",
    "- 12:30 - 12:45 - global variable dependence explanations `(R)`\n",
    "- 12:45 - 13:00 - summary and Q&A\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Technicalities\n",
    "### How to run R in Jupyter Notebook?\n",
    "1. Install R: https://cran.r-project.org/bin/. \n",
    "2. Use R terminal to install required packages: \n",
    "    ```R\n",
    "    install.packages(c('repr', 'IRdisplay', 'evaluate', 'crayon', 'pbdZMQ', 'devtools', 'uuid', 'digest'))\n",
    "    install.packages('IRkernel')\n",
    "    ```\n",
    "3. Make kernel available to Jupyter Notebook:\n",
    "    ```R\n",
    "    IRkernel::installspec()\n",
    "    ```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Installing and loading packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# install.packages(\"survival\")\n",
    "# install.packages(\"survex\")\n",
    "# install.packages(\"treeshap\")\n",
    "# install.packages(\"survminer\")\n",
    "# install.packages(\"ranger\")\n",
    "# install.packages(\"partykit\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "library(survival)\n",
    "library(survminer)\n",
    "library(survex)\n",
    "library(treeshap)\n",
    "library(partykit)\n",
    "library(ranger)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "options(jupyter.plot_scale=1, repr.plot.width = 6, repr.plot.height = 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Intro to survival analysis\n",
    "\n",
    "### What is survival analysis?\n",
    "- type of supervised learning task\n",
    "- also known as time-to-event analysis, reliability analysis, duration analysis, ... \n",
    "- **data modality:** mostly tabular, censored data (the most popular case is right censoring)\n",
    "\n",
    "<img src=\"images/censoring.png\" width=\"600\">\n",
    "\n",
    "- **output:** survival probability distribution (can be represented by survival function)\n",
    "\n",
    "<img src=\"images/survival_probs.png\" width=\"600\">\n",
    "\n",
    "\n",
    "### Mathematical background\n",
    "\n",
    "#### Data format\n",
    "- instance $i$: $(\\textbf{x}_i, y_i, \\delta_i)$\n",
    "- $\\textbf{x}_i$ is a vector of features\n",
    "-  $\\delta_i$ is the indicator of event of interest's occurrence (usually 0 is for censored data, 1 for event of interest's occurrence)\n",
    "-  $y_i$ stands for the observed time \n",
    "    - it can be either survival time $T_i$ when $\\delta_i = 1$ \n",
    "    - or censoring time $C_i$ when $\\delta_i = 0$\n",
    "- we want to model the conditional distribution of $T$, not $Y$ (but we don't observe $T$ directly!)\n",
    "\n",
    "\n",
    "#### Key quantities  \n",
    "- **survival function:** (probability of surviving past a given time point)\n",
    " $$S(t|\\mathbf{x}) = \\mathbb{P}(T > t|\\mathbf{X}=\\mathbf{x}) = 1 - \\mathbb{P}(T \\leq t|\\mathbf{X}=\\mathbf{x})$$\n",
    "\n",
    "\n",
    "- **hazard function:** (potential for an event to occur at infinitesimally small time interval $[t, \\Delta t]$)\n",
    "$$h(t|\\mathbf{x}) = \\lim_{\\Delta t \\to 0} \\frac{\\mathbb{P}(t \\leq T < t + \\Delta t|T \\geq t, \\mathbf{X}=\\mathbf{x})}{\\Delta t}$$\n",
    "\n",
    "- **cumulative hazard function:**  (total accumulated risk of an event up to time $t$)\n",
    "$$H(t|\\mathbf{x}) = \\int_0^t h(s|\\mathbf{x})ds$$\n",
    "\n",
    "- they are related by: $S(t|\\mathbf{x}) = \\exp(-H(t|\\mathbf{x}))$\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First look at the data & statistical models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "df <- read.csv(\"datasets/lung_dataset.csv\")\n",
    "df <- df[complete.cases(df),]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "dim(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "head(df, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data contain subjects with advanced lung cancer.\n",
    "\n",
    "- time: survival time (in days)\n",
    "- status: event status, 0=censored, 1=deceased\n",
    "- age: age (in years)\n",
    "- sex: 1=male, 2=female\n",
    "- ph.ecog: ECOG performance score as rated by the physician: \n",
    "    - 0=asymptomatic\n",
    "    - 1=symptomatic but completely ambulatory\n",
    "    - 2=in bed <50% of the day\n",
    "    - 3=in bed > 50% of the day but not bedbound\n",
    "    - 4=bedbound\n",
    "- ph.karno: Karnofsky performance score (0=bad, 100=good) rated by physician\n",
    "- pat.karno: Karnofsky performance score (0=bad, 100=good) rated by patient\n",
    "- meal.cal: calories consumed at meals\n",
    "- wt.loss: weight loss in last six months"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kaplan-Meier estimator \n",
    "- non-parametric estimator for survival function \n",
    "- estimates $S(t)$ only looking at the $(y_i, \\delta_i)$ pairs - do not take into account features $\\textbf{x}_i$\n",
    "- separate estimators can be fitted for each group of data instances with the same value of categorical feature (stratifying)\n",
    "\n",
    "Let's check how it looks like for our data, first for the small subset, then for the whole dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "tmp_df <- df[1:20, c(\"time\", \"status\")] \n",
    "tmp_df "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggsurvplot(survfit(Surv(time, status) ~ 1, data = tmp_df), \n",
    "    xlab = \"Time\", \n",
    "    ylab = \"Survival Probability\", \n",
    "    legend = \"none\",\n",
    "    conf.int = FALSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "ggsurvplot(survfit(Surv(time, status) ~ 1, data = df), \n",
    "    xlab = \"Time\", \n",
    "    ylab = \"Survival Probability\", \n",
    "    legend = \"none\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cox proportional hazards model \n",
    "- semi-parametric model for hazard function\n",
    "- assumes that the hazard function is a product of two parts: \n",
    "    - baseline hazard function $h_0(t)$ - does not depend on features $\\textbf{x}_i$\n",
    "    - exponential function of linear combination of features $\\textbf{x}_i$\n",
    "$$ h(t|\\mathbf{x}) = h_0(t) \\exp(\\beta_1 x_1 + \\beta_2 x_2 + \\dots + \\beta_p x_p)$$\n",
    "- good baseline model (~*linear regression of survival analysis*)\n",
    "- considered interpretable (exponentiated coefficients $\\exp(\\beta_j)$ can be interpreted as the hazard ratio)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "cph_model <- coxph(Surv(time, status) ~ ., data = df, x = TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "summary(cph_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "patient_A <- df[42, 3:9] \n",
    "patient_B <- patient_A\n",
    "patient_B[\"ph.ecog\"] <- 2  \n",
    "patients <- rbind(patient_A, patient_B)\n",
    "patients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "risks <- predict(cph_model, rbind(patient_A, patient_B), type = \"risk\")\n",
    "print(risks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "print(risks[2]/risks[1]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "times <- sort(unique(df$time))\n",
    "probs <- pec::predictSurvProb(cph_model, \n",
    "                                patients, \n",
    "                                times = times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "plot(times, probs[2,], type = \"l\", col = \"red\", \n",
    "    ylab = \"Survival Probability\", xlab = \"Time\")\n",
    "lines(times, probs[1,], col = \"blue\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Machine learning models\n",
    "\n",
    "### Survival trees\n",
    "- extension of decision trees for survival analysis\n",
    "- for each leaf node, based on the data instances in the node, one of the key quantities is calculated (e.g., cumulative hazard function)\n",
    "- there are many splitting rules but the most common is the log-rank splitting rule (based on the log-rank test comparing survival distributions of two samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "tree_model <- ctree(Surv(time, status) ~ ., data = df, )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "plot(tree_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random survival forest\n",
    "- works similarly to classical random forest but uses survival trees instead of classification or regression decision trees\n",
    "- predictions are based on the average of predictions of all trees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "rsf_model <- ranger(Surv(time, status) ~ ., data = df[complete.cases(df),])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "rsf_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "predictions <- predict(rsf_model, patients, type = \"response\")\n",
    "surv_probs <- predictions$survival\n",
    "times <- predict(rsf_model, patients, type = \"response\")$unique.death.times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "plot(times, surv_probs[2,], type = \"l\", col = \"red\", \n",
    "    ylab = \"Survival Probability\", xlab = \"Time\")\n",
    "lines(times, surv_probs[1,], col = \"blue\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"images/models.png\" width=\"1000\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Find more here:**\n",
    "- Wang, P., Li, Y., & Reddy, C. K. (2019). Machine learning for survival analysis: A survey. ACM Computing Surveys (CSUR), 51(6), 1-36. https://doi.org/10.1145/3214306\n",
    "- Wiegrebe, S., Kopper, P., Sonabend, R., & Bender, A. (2023). Deep Learning for Survival Analysis: A Review. arXiv preprint arXiv:2305.14961. https://arxiv.org/abs/2305.14961"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.3.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}