Chapter_11_sec_11.8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 11.8 Lab: Survival Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from lifelines import KaplanMeierFitter\n",
    "from lifelines.statistics import logrank_test\n",
    "from lifelines import CoxPHFitter\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" \n",
    "I am also new to this topic, so let us learn those concept together. Feedbacks are welcome.\n",
    "Survival analysis is a statistical method used to estimate the survival function of a population.\n",
    "These arise in the analysis of a unique kind of outcome variable: the analysis time until an event occurs.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 11.8.1 Brain Cancer Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BrainCancer = pd.read_csv('data/BrainCancer.csv', header=0)\n",
    "# use some options in .describe() to get a quick overview of the data\n",
    "BrainCancer.describe(include = 'object')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BrainCancer.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# I did a bit of google search and found the package lifelines.\n",
    "# % pip install lifelines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a kmf object\n",
    "kmf = KaplanMeierFitter() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fit the data into the model\n",
    "kmf.fit(BrainCancer.time, BrainCancer.status,label='Kaplan Meier Estimate')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create an estimate\n",
    "kmf.plot(ci_show=True) ## ci_show is meant for Confidence interval, which is the shaded area in the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kmf1 = KaplanMeierFitter() ## instantiate the class to create an object\n",
    "\n",
    "## two Cohorts are compared. Cohort 1. Female; Cohort 2. Male \n",
    "groups = BrainCancer['sex'] \n",
    "T = BrainCancer.time\n",
    "E = BrainCancer.status  \n",
    "i1 = (groups == 'Female')     \n",
    "i2 = (groups == 'Male')     \n",
    "\n",
    "\n",
    "## fit the model for 1st cohort\n",
    "kmf1.fit(T[i1], E[i1], label='Female')\n",
    "a1 = kmf1.plot()\n",
    "\n",
    "## fit the model for 2nd cohort\n",
    "kmf1.fit(T[i2], E[i2], label='Male')\n",
    "kmf1.plot(ax=a1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can perform a log-rank test to compare the survival of males to females,\n",
    "results=logrank_test(T[i1],T[i2],event_observed_A=E[i1], event_observed_B=E[i2])\n",
    "results.print_summary()\n",
    "\"\"\" \n",
    "The resulting p-value is 0.23, indicating no evidence of a difference in sur- vival between the two sexes.\n",
    "This also can be seen from the overlapping the confidence intervals. \n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_dummy = pd.get_dummies(BrainCancer, drop_first=True)\n",
    "df_dummy.head()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use Cox Proportional Hazards model\n",
    "cph1 = CoxPHFitter()   \n",
    "cph1.fit(df_dummy[['status', 'sex_Male', 'time']], 'time', event_col='status')   \n",
    "cph1.print_summary()\n",
    "\"\"\" \n",
    "I was not able to use the optional 'formula' in the fit() function due to the compatibility with new version, instead I sliced the dataframe to only contains 3 cols ('status', 'sex_Male', 'time')\n",
    "Here the p value is at 0.233. Regardless of which test we use, we see that there is no clear evidence for a difference in survival between males and females.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use Cox Proportional Hazards model and use more features \n",
    "cph2 = CoxPHFitter()   \n",
    "cph2.fit(df_dummy, 'time', event_col='status')   \n",
    "cph2.print_summary()\n",
    "\"\"\" \n",
    "after adjusting for the other predictors, larger values of the Karnofsky index, ki, are associated with lower risk, \n",
    "i.e. longer survival.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 11.8.2 Publication Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Publication = pd.read_csv('data/Publication.csv', header=0)\n",
    "Publication.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kmf1 = KaplanMeierFitter() ## instantiate the class to create an object\n",
    "\n",
    "## two Cohorts are compared. Cohort 1. posres==0 ; Cohort 2. posres==1\n",
    "groups = Publication['posres'] \n",
    "T = Publication.time\n",
    "E = Publication.status  \n",
    "i1 = (groups == 0)     \n",
    "i2 = (groups == 1)     \n",
    "\n",
    "## fit the model for 1st cohort\n",
    "kmf1.fit(T[i1], E[i1], label='Positive Results')\n",
    "a1 = kmf1.plot()\n",
    "\n",
    "## fit the model for 2nd cohort\n",
    "kmf1.fit(T[i2], E[i2], label='Negative Results')\n",
    "kmf1.plot(ax=a1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can perform a log-rank test to compare the survival of males to females,\n",
    "results=logrank_test(T[i1],T[i2],event_observed_A=E[i1], event_observed_B=E[i2])\n",
    "results.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_dummy = pd.get_dummies(Publication, drop_first=True)\n",
    "df_dummy.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use Cox Proportional Hazards model\n",
    "cph3 = CoxPHFitter()   \n",
    "cph3.fit(df_dummy[['status', 'posres', 'time']], 'time', event_col='status')   \n",
    "cph3.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cph4 = CoxPHFitter()   \n",
    "cph4.fit(df_dummy[['status', 'posres', 'time', 'multi', 'clinend','sampsize', 'budget', 'impact']], 'time', event_col='status')   \n",
    "cph4.print_summary()\n",
    "\"\"\" \n",
    "After we control for other features, posres becomes an important factor (well, at least significant).\n",
    "We see that there are a number of statistically significant variables, \n",
    "including whether the trial focused on a clinical endpoint (clinend), the impact of the study(impact),\n",
    "and whether the study had positive or negative results (posres).\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 11.8.3 Call Center Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1)\n",
    "N = 2000\n",
    "Operators = np.random.choice(range(5, 16), N)\n",
    "Center = np.random.choice([\"A\", \"B\", \"C\"], N)\n",
    "Time = np.random.choice([\"Morn.\", \"After.\", \"Even.\"], N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we generate a similar random data set\n",
    "X_pre = pd.DataFrame({\"Operators\": Operators, \"Center\": Center, \"Time\": Time})\n",
    "X = pd.get_dummies(X_pre, drop_first=True)\n",
    "X.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "true_coeff = np.array([0.04, -0.3, 0, 0.2, -0.2])\n",
    "# well, I was not able to fully following the simulation in the book. \n",
    "# I think the highlevel idea is to use those coefficients to generate a dataset and show the model fit could \n",
    "# sucessfully recover the coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this simulation code is wrong. And I will come back and fix it. \n",
    "X['y'] = 350*np.exp(np.sum(-X*true_coeff,axis=1).tolist() + np.random.normal(0, 0.005, N))\n",
    "X['answered'] = np.where( X['y'] < 300 , 1, 0)\n",
    "X.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use Cox Proportional Hazards model\n",
    "cph5 = CoxPHFitter()   \n",
    "cph5.fit(X, 'y', event_col='answered')   \n",
    "cph5.print_summary()\n",
    "\"\"\" \n",
    "Since the simulation is wrong, so the summary is not correct. \n",
    "But at least we can see the p-values for Operatator, Center = B, Time = Even. and Time = Morn are very small, \n",
    "and they are directly related to the ground truth coefficients.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# End of Chapter 11"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "4548a0e672c5b3a287feee7b2962606840aa548749d1830ef724408652b0c250"
  },
  "kernelspec": {
   "display_name": "Python 2.7.16 64-bit ('base': conda)",
   "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}