{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to Survival Analysis with scikit-survival"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**scikit-survival** is a Python module for [survival analysis](https://en.wikipedia.org/wiki/Survival_analysis) built on top of [scikit-learn](http://scikit-learn.org/). It allows doing survival analysis while utilizing the power of scikit-learn, e.g., for pre-processing or doing cross-validation.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "\n",
    "1. [What is Survival Analysis?](#What-is-Survival-Analysis?)\n",
    "2. [The Veterans' Administration Lung Cancer Trial](#The-Veterans%27-Administration-Lung-Cancer-Trial)\n",
    "3. [Survival Data](#Survival-Data)\n",
    "4. [The Survival Function](#The-Survival-Function)\n",
    "5. [Considering other variables by stratification](#Considering-other-variables-by-stratification)\n",
    "6. [Multivariate Survival Models](#Multivariate-Survival-Models)\n",
    "7. [Measuring the Performance of Survival Models](#Measuring-the-Performance-of-Survival-Models)\n",
    "8. [Feature Selection: Which Variable is Most Predictive?](#Feature-Selection:-Which-Variable-is-Most-Predictive?)\n",
    "9. [What's next?](#What%27s-next?)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is Survival Analysis?\n",
    "\n",
    "The objective in survival analysis — also referred to as reliability analysis in engineering — is to establish a connection between covariates and the time of an event. The name *survival analysis* originates from clinical research, where predicting the time to death, i.e., survival, is often the main objective. Survival analysis is a type of regression problem (one wants to predict a continuous value), but with a twist. It differs from traditional regression by the fact that parts of the training data can only be partially observed – they are *censored*.\n",
    "\n",
    "As an example, consider a clinical study, which investigates coronary heart disease and has been carried out over a 1 year period as in the figure below.\n",
    "\n",
    "![image censoring](https://k-d-w.org/clipboard/censoring.png)\n",
    "\n",
    "Patient A was lost to follow-up after three months with no recorded cardiovascular event, patient B experienced an event four and a half months after enrollment, patient D withdrew from the study two months after enrollment, and patient E did not experience any event before the study ended. Consequently, the exact time of a cardiovascular event could only be recorded for patients B and C; their records are *uncensored*. For the remaining patients it is unknown whether they did or did not experience an event after termination of the study. The only valid information that is available for patients A, D, and E is that they were event-free up to their last follow-up. Therefore, their records are *censored*.\n",
    "\n",
    "Formally, each patient record consists of a set of covariates $x \\in \\mathbb{R}^d$ , and the time $t>0$ when an event occurred or the time $c>0$ of censoring. Since censoring and experiencing and event are mutually exclusive, it is common to define an event indicator $\\delta \\in \\{0;1\\}$ and the observable survival time $y>0$. The observable time $y$ of a right censored sample is defined as\n",
    "\n",
    "$$\n",
    "y = \\min(t, c) = \n",
    "\\begin{cases} \n",
    "t & \\text{if } \\delta = 1 , \\\\ \n",
    "c & \\text{if } \\delta = 0 .\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Consequently, survival analysis demands for models that take this unique characteristic of such a dataset into account, some of which are showcased below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Veterans' Administration Lung Cancer Trial\n",
    "\n",
    "The Veterans' Administration Lung Cancer Trial is a randomized trial of two treatment regimens for lung cancer. The [data set](http://lib.stat.cmu.edu/datasets/veteran) (Kalbfleisch J. and Prentice R, (1980) The Statistical Analysis of Failure Time Data. New York: Wiley) consists of 137 patients and 8 variables, which are discribed below:\n",
    "\n",
    "- `Treatment`: denotes the type of lung cancer treatment; `standard` and `test` drug.\n",
    "- `Celltype`: denotes the type of cell involved; `squamous`, `small cell`, `adeno`, `large`.\n",
    "- `Karnofsky_score`: is the Karnofsky score.\n",
    "- `Diag`: is the time since diagnosis in months.\n",
    "- `Age`: is the age in years.\n",
    "- `Prior_Therapy`: denotes any prior therapy; `none` or `yes`.\n",
    "- `Status`: denotes the status of the patient as dead or alive; `dead` or `alive`.\n",
    "- `Survival_in_days`: is the survival time in days since the treatment.\n",
    "\n",
    "Our primary interest is studying whether there a subgroups that differ in survival and whether we can predict survival times."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Survival Data\n",
    "\n",
    "As described in the section *What is Survival Analysis?* above, survival times are subject to right-censoring, therefore, we need to consider an individual's status in addition to survival time. To be fully compatible with scikit-learn, `Status` and `Survival_in_days` need to be stored as a [structured array](https://docs.scipy.org/doc/numpy/user/basics.rec.html) with the first field indicating whether the actual survival time was observed or if was censored, and the second field denoting the observerd survival time, which corresponds to the time of death (if `Status == 'dead'`, $\\delta = 1$) or the last time that person was contacted (if `Status == 'alive'`, $\\delta = 0$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([( True,  72.), ( True, 411.), ( True, 228.), ( True, 126.),\n",
       "       ( True, 118.), ( True,  10.), ( True,  82.), ( True, 110.),\n",
       "       ( True, 314.), (False, 100.), ( True,  42.), ( True,   8.),\n",
       "       ( True, 144.), (False,  25.), ( True,  11.), ( True,  30.),\n",
       "       ( True, 384.), ( True,   4.), ( True,  54.), ( True,  13.),\n",
       "       (False, 123.), (False,  97.), ( True, 153.), ( True,  59.),\n",
       "       ( True, 117.), ( True,  16.), ( True, 151.), ( True,  22.),\n",
       "       ( True,  56.), ( True,  21.), ( True,  18.), ( True, 139.),\n",
       "       ( True,  20.), ( True,  31.), ( True,  52.), ( True, 287.),\n",
       "       ( True,  18.), ( True,  51.), ( True, 122.), ( True,  27.),\n",
       "       ( True,  54.), ( True,   7.), ( True,  63.), ( True, 392.),\n",
       "       ( True,  10.), ( True,   8.), ( True,  92.), ( True,  35.),\n",
       "       ( True, 117.), ( True, 132.), ( True,  12.), ( True, 162.),\n",
       "       ( True,   3.), ( True,  95.), ( True, 177.), ( True, 162.),\n",
       "       ( True, 216.), ( True, 553.), ( True, 278.), ( True,  12.),\n",
       "       ( True, 260.), ( True, 200.), ( True, 156.), (False, 182.),\n",
       "       ( True, 143.), ( True, 105.), ( True, 103.), ( True, 250.),\n",
       "       ( True, 100.), ( True, 999.), ( True, 112.), (False,  87.),\n",
       "       (False, 231.), ( True, 242.), ( True, 991.), ( True, 111.),\n",
       "       ( True,   1.), ( True, 587.), ( True, 389.), ( True,  33.),\n",
       "       ( True,  25.), ( True, 357.), ( True, 467.), ( True, 201.),\n",
       "       ( True,   1.), ( True,  30.), ( True,  44.), ( True, 283.),\n",
       "       ( True,  15.), ( True,  25.), (False, 103.), ( True,  21.),\n",
       "       ( True,  13.), ( True,  87.), ( True,   2.), ( True,  20.),\n",
       "       ( True,   7.), ( True,  24.), ( True,  99.), ( True,   8.),\n",
       "       ( True,  99.), ( True,  61.), ( True,  25.), ( True,  95.),\n",
       "       ( True,  80.), ( True,  51.), ( True,  29.), ( True,  24.),\n",
       "       ( True,  18.), (False,  83.), ( True,  31.), ( True,  51.),\n",
       "       ( True,  90.), ( True,  52.), ( True,  73.), ( True,   8.),\n",
       "       ( True,  36.), ( True,  48.), ( True,   7.), ( True, 140.),\n",
       "       ( True, 186.), ( True,  84.), ( True,  19.), ( True,  45.),\n",
       "       ( True,  80.), ( True,  52.), ( True, 164.), ( True,  19.),\n",
       "       ( True,  53.), ( True,  15.), ( True,  43.), ( True, 340.),\n",
       "       ( True, 133.), ( True, 111.), ( True, 231.), ( True, 378.),\n",
       "       ( True,  49.)],\n",
       "      dtype=[('Status', '?'), ('Survival_in_days', '<f8')])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sksurv.datasets import load_veterans_lung_cancer\n",
    "\n",
    "data_x, data_y = load_veterans_lung_cancer()\n",
    "data_y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily see that only a small number of survival times are right-censored (`Status` is `False`), i.e., most veteran's died during the study period (`Status` is `True`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Survival Function\n",
    "\n",
    "A key quantity in survival analysis is the so-called survival function, which relates time to the probability of surviving beyond a given time point.\n",
    "\n",
    "> Let $T$ denote a continuous non-negative random variable corresponding to a patient’s survival time. The survival function $S(t)$ returns the probability of survival beyond time $t$ and is defined as\n",
    "> $$ S(t) = P (T > t). $$\n",
    "\n",
    "If we observed the exact survival time of all subjects, i.e., everyone died before the study ended, the survival function at time $t$ can simply be estimated by the ratio of patients surviving beyond time $t$ and the total number of patients:\n",
    "\n",
    "$$\n",
    "\\hat{S}(t) = \\frac{ \\text{number of patients surviving beyond $t$} }{ \\text{total number of patients} }\n",
    "$$\n",
    "\n",
    "In the presence of censoring, this estimator cannot be used, because the numerator is not always defined. For instance, consider the following set of patients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Status</th>\n",
       "      <th>Survival_in_days</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>True</td>\n",
       "      <td>8.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>True</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>True</td>\n",
       "      <td>20.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>False</td>\n",
       "      <td>25.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>True</td>\n",
       "      <td>59.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Status  Survival_in_days\n",
       "1    True               8.0\n",
       "2    True              10.0\n",
       "3    True              20.0\n",
       "4   False              25.0\n",
       "5    True              59.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "pd.DataFrame.from_records(data_y[[11, 5, 32, 13, 23]], index=range(1, 6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the formula from above, we can compute $\\hat{S}(t=11) = \\frac{3}{5}$, but not $\\hat{S}(t=30)$, because we don't know whether the 4th patient is still alive at $t = 30$, all we know is that when we last checked at $t = 25$, the patient was still alive.\n",
    "\n",
    "An estimator, similar to the one above, that *is* valid if survival times are right-censored is the [Kaplan-Meier estimator](https://en.wikipedia.org/wiki/Kaplan%E2%80%93Meier_estimator)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'time $t$')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEMCAYAAAA4S+qsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHGBJREFUeJzt3X2UXXV97/H3Rx5bJEJM0JgwJHAjEi0CjpmgdBUENKAl2EsJpCwpoqEISq21F7gUkS7r40WrFzFZPkLlSRRJaSpairViGBMeTMgEykgEBoJEAyRFuSXyvX/sfZKdk5lk7z37PH9ea501Z+/zm32+ezbm6+9ZEYGZmVkeL2l1AGZm1jmcNMzMLDcnDTMzy81Jw8zMcnPSMDOz3Jw0zMwsNycNMzPLzUnDzMxyc9IwM7Pcdm11AFWbNGlSTJ8+vdVhmJl1lLvvvvtXETF5Z+W6LmlMnz6dFStWtDoMM7OOIumRPOXcPGVmZrk5aZiZWW4tTRqSvirpKUn3j/G5JH1e0rCklZKOaHaMZma2VatrGl8H5u7g8xOAmelrIXBVE2IyM7MxtDRpRMSPgA07KDIPuDoSdwH7SJrSnOjMzKxeq2saOzMVeCxzPJKe24akhZJWSFqxfv36pgVnZtZr2j1paJRz2201GBGLI6I/IvonT97pMGMzMyup3ZPGCLB/5nga8ESjvuyj/7Saj/7T6kZd3sys47V70lgCvCsdRTUHeDYi1jXqy4ae2MjQExsbdXkzs47X0hnhkq4DjgYmSRoBPgLsBhARXwKWAicCw8BvgLNaE6mZmUGLk0ZEnL6TzwM4r0nhmJnZTrR785SZmbURJ406g2s3cO3go60Ow8ysLTlpZMw7LJkCcst9j7c4EjOz9uSkkbFgoI+BGRNbHYaZWdty0jAzs9ycNEbhfg0zs9E5adRxv4aZ2dicNOq4X8PMbGxOGmZmlpuThpmZ5eakYWZmuTlpmJlZbk4aZmaWm5OGmZnl5qQxBk/wMzPbnpPGKDzBz8xsdE4ao/AEPzOz0TlpmJlZbk4aZmaWm5OGmZnl5qRhZma5OWnsgIfdmplta1xJQ9JeknapKph24mG3ZmbbK5Q0JL1E0gJJ/yzpKeABYJ2k1ZI+LWlmY8JsPg+7NTPbXtGaxh3AQcBFwCsjYv+I2A/4Q+Au4BOSzqg4RjMzaxO7Fix/XES8IOmAiHixdjIiNgDfBr4tabdKIzQzs7ZRqKYRES+kb2+u/0zSnLoyXcGd4WZmWxXt0zhV0ieAvSUdUtcJvrja0FrPneFmZtsq2qdxJzAE7AtcATwk6R5JtwK/rTq4VnNnuJnZtgr1aUTE48DVkn4eEXcCSJoIzCAZSWVmZl2sUNKQpEjcWTuXdoJvqC9TYYxmZtYmCg+5lfR+SX3Zk5J2l/QWSd8Azsx7MUlzJT0oaVjShaN83ifpDkn3Slop6cSC8ZqZWYWKDrmdC7wbuE7SDOAZ4PdIks/3gc9GxH15LpR2ol8JHA+MAMslLYmIoUyxS4AbI+IqSbOApcD0gjGbmVlFivZpPA98EfhiOh9jEvDbiHimxHfPBoYj4mEASdcD80g62rd8JTAhff8y4IkS32NmZhXZYfOUpFmS/jFzfLuk18KW+RhvBM6XNLvEd08FHsscj6Tnsi4DzpA0QlLLeP8YcS6UtELSivXr15cIxczM8thZn8btJE1ENdMiYjWApDcB1wB9wNclvbPgd2uUc/Ud6KcDX4+IacCJwDWStos5IhZHRH9E9E+ePLlgGGZmltfOksZbgY9ljjdm3r8L+FJELASOBv5Xwe8eAfbPHE9j++ans4EbASJiGbAnSZNYU3lWuJlZYodJIyJWRcSfZU4NSzpF0n7AycAtabmngD0KfvdyYKakGZJ2B04DltSVeRQ4FkDSISRJo6ntT54Vbma2VdEhtx8EzgEeB+6JiJ8ApJ3iLy1yoYjYDJwP3AasIRkltVrS5ZJOSot9CHivpJ8B1wF/3uw5IJ4Vbma2VdHRU08Cx0t6SXaVW+AYkmXTC4mIpSQd3Nlzl2beDwFvLnpdMzNrjFI799UlDCLi+2nfRtdyv4aZmfcIz8X9GmZmCSeNHNyvYWaWcNIwM7Pciq5yu4ntJ+BBMlEvImLCKJ+ZmVmXKDp6au9GBWJmZu2v6Cq3W0jaF5hJMuEOgIj4URVBmZlZeyqVNCS9B7iAZOmP+4A5wDLgLdWFZmZm7aZsR/gFJCvcPhIRxwCH0+TlPczMrPnKJo3n0701kLRHRDwAHFxdWO3JE/zMrNeVTRojkvYBvgv8QNItdPkGSZ7gZ2ZWsk8jImp7Z1wm6Q6SXfW+V1lUbWjBQJ8Thpn1vLId4R8EvhURIxHx7xXHZGZmbaps89QE4DZJ/yHpPEmvqDKoduZ+DTPrZWVXuf1oRLwWOA94FfDvkv610sjakPs1zKzXjXftqaeAJ4FfA/uNP5z25oULzazXlUoaks6V9EPgdpI9u98bEYdWGZiZmbWfssuIHAD8ZUTcV2UwZmbW3sr2aVzYywnDneFm1qsKJQ1JP05/bpK0MfPaJGljY0JsL+4MN7NeVnRp9KPSnz27RLon+ZlZLyvbEf5BSVOrDqaTDK3byPxFy9xMZWY9pWxH+ATg+5I2ANcDN0XEL6sLq73VmqiG1iUtcgsG+loZjplZ03hyXwkLBvq44ZwjmTXFu9uaWW/x5L5x8kgqM+slntw3Dh5JZWa9pnCfhiQB/Xhyn0dSmVnPKVzTiIgADu/1hJHlJioz6xVl+zSWSXpjpZF0KDdRmVkvKZs0jgHukvRzSSslrZK0ssrAOoVXvjWzXlJ2nsYJlUbRBWpNVJ6zYWbdrGxN48wxXoVImivpQUnDki4co8ypkoYkrZZ0bcl4G8pNVGbWK8rWNJ7LvN8TeAewpsgFJO0CXAkcD4wAyyUtiYihTJmZwEXAmyPiaUltORfEo6jMrFeUShoR8X+yx5I+AywpeJnZwHBEPJxe43pgHjCUKfNe4MqIeDr93qfKxGtmZtUY74zwmt8HDiz4O1OBxzLHI+m5rFcDr5Z0p6S7JM0dR4xmZjZOpWoaklYBkR7uAkwGLi96mVHORd3xrsBM4GhgGvAfkl4XEc/UxbMQWAjQ1+eOaDOzRinbp/GOzPvNwC8jYnPBa4wA+2eOpwFPjFLmroh4AVgr6UGSJLI8WygiFgOLAfr7++sTj5mZVaRs89RsYENEPAKcBdwo6YiC11gOzJQ0Q9LuwGls3y/yXZI5IUiaRNJc9XDJmM3MbJzKJo2/jYhNko4C3gZ8A7iqyAXSmsn5wG0kI69ujIjVki6XdFJa7Dbg15KGgDuAD0fEr0vGbGZm41S2eep36c+3A1dFxC2SLit6kYhYCiytO3dp5n0Af5W+zMysxcrWNB6XtAg4FVgqaY9xXKtreOFCM+t2Zf+hP5Wk6WhuOpJpIvDhyqLqQJ4Vbma9oOzkvt8A38kcrwPWVRVUJ/KscDPrBT3fpGRmZvk5aVTM/Rpm1s0KJQ1J16Q/L2hMOJ3N/Rpm1u2K1jTeIOkA4N2S9pU0MftqRICdpLYhk2sbZtatiiaNLwHfA14D3F33WlFtaJ3JtQ0z62aFkkZEfD4iDgG+GhEHRsSMzKvoKrddydu/mlk3Kzvk9lxJrwf+MD31o4joyT3CxzK0biPzFy0DktqHt4E1s25Qdmn0D5AsRV6bq/FNSYsj4guVRdbBak1UkIymGly7AcCJw8w6npLlnQr+krQSODIinkuP9wKWRcShFcdXWH9/f6xY0T7dK9cOPsrFN69iYMZEbjjnyFaHY2Y2Kkl3R0T/zsqVnachti5aSPp+tE2Vel6tj6PWXOVRVWbWycqucvs1YFDSzenxycBXqgmp+9Saq4bWbQTcTGVmnatUTSMiriDZfGkD8DRwVkR8rsrAusmCgT5uOOdIZk2Z0OpQzMzGpWxNg4i4B7inwljMzKzNee2pJvNscTPrZE4aTeTZ4mbW6UolDUnnS9q36mC6ndemMrNOV7am8UpguaQbJc2V5OG2OdVqGxffvMqJw8w6TtnRU5cAM0mG2f458JCkv5d0UIWxdaUFA338/Tv/AHAzlZl1ntJ9GpFMJX8yfW0G9gVukvSpimLrWl7U0Mw6Vdk+jQ9Iuhv4FHAn8AcRcS7wBuB/VhhfV3Pfhpl1mrI1jUnAn0TE2yLiWxHxAkBEvAi8o7LouphHUplZJyqbNPaIiEeyJyR9EiAi1ow7qh7gJioz60Rlk8bxo5w7YTyB9Co3UZlZJymUNCSdK2kVcLCklZnXWsCbMBXkJioz6zRF1566FvgX4OPAhZnzmyJiQ2VR9YgFA31OGGbWUQoljYh4FngWOL0x4ZiZWTsr2jz14/TnJkkbM69NkjY2JsTu534NM+sUhZJGRByV/tw7IiZkXntHhDeLKMHLiphZJ2npKrfpulUPShqWdOEOyp0iKSTtdP/aTuNlRcyskxTq05C0CQhG3w88itQ2JO0CXEkyfHeEZAHEJRExVFdub+ADwGCRWDuJO8TNrFMUbZ7aO9McVf8q2jw1GxiOiIcj4r+B64F5o5T7O5LlSp4veP2OM7RuI/MXLXMzlZm1rao6wjeW6AifCjyWOR5Jz2W/73Bg/4i4teC1O868w6Yya8oEhtZtdK3DzNpWVR3hE0rUNEZt4tryofQS4LPAh3Z6IWmhpBWSVqxfv75gGO1hwUAfN5xzJLOmeDyBmbWvVnaEjwD7Z46nAU9kjvcGXgf8UNIvgDnAktE6wyNicUT0R0T/5MmTGxiymVlvKzojHABJewLvA44iqR38GLgqIor0OywHZkqaATwOnAYsqH2YTiSclPnOHwJ/HRErysRsZmbjV7amcTXwWuALwP8FDgGuKXKBiNgMnA/cBqwBboyI1ZIul3RSybjMzKyBStU0gIMj4vWZ4zsk/azoRSJiKbC07tylY5Q9uuj1zcysWmVrGvdKmlM7kDRAsoOfVcBDb82sXRWd3LeKpA9jN+Bdkmr/qvUBQ2P+ouVWW1ZkaF0ygnnBQF8rwzEz20bR5ilv5dpgCwb6WDDQx/xFy1odipnZdooujb5li1dJ+wIzgT0zRR7Z7pfMzKxrlB1y+x7gApK5FfeRzKFYBrylutDMzKzdlO0IvwB4I/BIRBwDHA505lTsNuZ9Nsys3ZRNGs/XJvJJ2iMiHgAOri4sy+6z4ZFUZtYuys7TGJG0D/Bd4AeSnmbbJUBsnGqjpm6573GPpDKztlGqphER74yIZyLiMuBvga8AJ1cZmG27iKGbqsysHVS59lRLdwHsZvMOm8rg2g1cfPMqwDUOM2udlq09Zfl5S1gzaxctXXvK8qttCVtbYmTeYVNd4zCzpiubNO6VNCci7gKvPdUsXmLEzFrNa091EC8xYmat5rWnOlRtNJVrG2bWTEX3CH+k9gL2Af44fe2TXZfKGqvWTOVOcTNrtlKjpyRdAHwT2C99/aOk91cZmI1twUAfAzMmtjoMM+tBZTvCzwYGIuI5AEmfJFmw8AtVBWZmZu2n7DwNAb/LHP8uPWdN5B3+zKzZytY0vgYMSro5PT6ZZCkRa5Jav8bg2g0Mrt0AeAiumTVe2bWnrgDOAjYATwNnRcTnqgzMdqy2LpVniptZMxWuaUgSMC0i7gHuqT4kK6I2U9xDcM2sGQrXNCIiSJZEtzbhIbhm1ixlO8LvkvTGSiOx0mpDcN0xbmaNVrYj/BjgLyT9AniOZORURMShVQVmxXhdKjNrhrJJ44RKo7Bx87pUZtYMZZPGL9l+E6arqgrKzMzaU9mkcTWwia0zwE8n2YTpT6sIyszM2pM3YepCtQ7xGm/YZGZV8SZMXabWIV5TmzFeG47rBGJm41E2aQyw/SZMa2qbNHkUVevUOsRrrh18dEvC8MgqMxsvJXP1Cv6SdMCOPs+7t4akucA/ALsAX46IT9R9/lfAe4DNwHrg3Tu7dn9/f6xYsSLP1/ec+YuWMbRuI7OmTNhyzjUPMwOQdHdE9O+sXKmaRhUbLknaBbgSOB4YAZZLWhIR2W1j7wX6I+I3ks4FPgXMH+9396r6pivXPMysqLLNU1WYDQxHxMMAkq4H5pHZazwi7siUvws4o6kRdpn6pqv5i5Z5zSozK6SVSWMq8FjmeISkr2QsZwP/0tCIesy8w6YyuHYDF9+8apt1q9xkZWZjqSxpSHplRDxZ5FdGOTdqB4ukM4B+4I/G+HwhsBCgr8//2OVVSwzZhOEmKzPbkVId4aNeSPrniHh7gfJHApdFxNvS44sAIuLjdeWOI5lE+EcR8dTOruuO8PFxZ7lZb2poR/hoiiSM1HJgpqQZwOPAacCCbAFJhwOLgLl5EoaNnzvLzWxHSi2NLumTec7tSERsBs4HbgPWADdGxGpJl0s6KS32aeClwLck3SdpSZl4Lb/ajoC116wpE7Z0lpuZld1P4/hRzhVe+TYilkbEqyPioIj4WHru0ohYkr4/LiJeERGHpa+TdnxFq5o3eDKzrELNU+lcifcBB0lamflob+AnVQZm7aG2nWx2PSv3cZj1rqJ9GteSDHv9OHBh5vymiNhQWVTWVrL9HNm1rJw8zHpPoaQREc8Cz0r6DrAhIjZJugQ4QtLfRcS9DYnSWio7KbC2lpU7yM16U9m1p1ZGxKGSjiKpdXwGuDgidjQ5ryk85LY5PDTXrLvkHXJbtiP8d+nPtwNXRcQtwO4lr2UdaN5hU7dJGEPrNrqz3KwHlK1p3Eoyt+J44Ajgt8BP6zZmagnXNFpjtJoHuPZh1ikaPbnvVGAu8JmIeEbSFODDJa9lXaB+UiB4YqBZNyqbNH4L7EWyN/jlwG7AM1UFZZ2nfgVd2Fr7mL9omWscZl2ibJ/GF4E5JEkDYBPJ3hhmW9T6PdzfYdY9yiaNgYg4D3geICKexh3hVqe2JImXIjHrHmWbp15Id94LAEmTgRcri8q6ylj7doxWzk1YZu2tbNL4PHAzsJ+kjwGnAJdUFpV1ldH27ajnTnOzzlB6Pw1JrwGOJdlM6faIWFNlYGV5yG1nGmvIbpZrImaN0/D9NCLiAeCBsr9vljXakN0s10TM2kNlO/e1C9c0upMnD5o1VtN37jNrJE8eNGsPrmlYx9pZP4hrIWb5uaZhXW9H/SCuhZg1hpOGdazRli6pmb9o2ZYJhU4cZtUpOyPcrK3VaiEX37yK+YuWeTa6WUVc07CulJ1Q6KYqs+q4pmFdy2tfmVXPScO6Xq2pyivtmo2fm6es6y0Y6NvSTDV/0bLtPvfQXLP8nDSsJ4w1PNf9HWbFeHKf9bQ8CyXm4dqKdTpP7jPLYWcLJebh2or1Etc0zMapqtpKHq7RWKO4pmHWJFXUVvIYXLuBwbUbco0Cc3KxRnHSMBunHS1nUqVrBx/NlTDcXGaN5OYpsy7TzOYya72qapV5m6daOrlP0lxJD0oalnThKJ/vIemG9PNBSdObH6VZZ5l32FQnjB4xtG5j0yettqymIWkX4D+B44ERYDlwekQMZcq8Dzg0Iv5C0mnAOyNi/o6u65qGmfWK+lrleGodnVDTmA0MR8TDEfHfwPXAvLoy84BvpO9vAo6VpCbGaGbWtrK1ymbVOlrZET4VeCxzPAIMjFUmIjZLehZ4OfCrpkRoZtbGsoMwRlsipxFamTRGqzHUt5XlKYOkhcBCgL4+jxgxs94z61XN6cdqZdIYAfbPHE8DnhijzIikXYGXARvqLxQRi4HFkPRpNCRaM7M29pE/fm1TvqeVfRrLgZmSZkjaHTgNWFJXZglwZvr+FODfotvGCJuZdZCW1TTSPorzgduAXYCvRsRqSZcDKyJiCfAV4BpJwyQ1jNNaFa+ZmbV4RnhELAWW1p27NPP+eeBPmx2XmZmNzjv3mZlZbk4aZmaWm5OGmZnl5qRhZma5OWmYmVluXbc0uqT1wCMlf30SvbdEie+5N/iee8N47vmAiJi8s0JdlzTGQ9KKPKs8dhPfc2/wPfeGZtyzm6fMzCw3Jw0zM8vNSWNbi1sdQAv4nnuD77k3NPye3adhZma5uaZhZma5OWmkJM2V9KCkYUkXtjqeqkjaX9IdktZIWi3pgvT8REk/kPRQ+nPf9LwkfT79O6yUdERr76AcSbtIulfSrenxDEmD6f3ekC7Hj6Q90uPh9PPprYx7PCTtI+kmSQ+kz/vIbn7Okj6Y/jd9v6TrJO3Zjc9Z0lclPSXp/sy5ws9V0plp+YcknTnad+XhpEHyDwxwJXACMAs4XdKs1kZVmc3AhyLiEGAOcF56bxcCt0fETOD29BiSv8HM9LUQuKr5IVfiAmBN5viTwGfT+30aODs9fzbwdET8D+CzablO9Q/A9yLiNcDrSe6/K5+zpKnAB4D+iHgdyfYKp9Gdz/nrwNy6c4Weq6SJwEdIttSeDXyklmgKi4iefwFHArdlji8CLmp1XA2611uA44EHgSnpuSnAg+n7RcDpmfJbynXKi2QXyNuBtwC3kmwb/Ctg1/rnTbKfy5Hp+13Tcmr1PZS45wnA2vrYu/U5A1OBx4CJ6XO7FXhbtz5nYDpwf9nnCpwOLMqc36ZckZdrGonaf4A1I+m5rpJWyQ8HBoFXRMQ6gPTnfmmxbvhbfA74G+DF9PjlwDMRsTk9zt7TlvtNP382Ld9pDgTWA19Lm+W+LGkvuvQ5R8TjwGeAR4F1JM/tbrr/OdcUfa6VPW8njYRGOddVw8okvRT4NvCXEbFxR0VHOdcxfwtJ7wCeioi7s6dHKRo5PuskuwJHAFdFxOHAc2xtshhNR9932rQyD5gBvArYi6Rppl63PeedGes+K7t/J43ECLB/5nga8ESLYqmcpN1IEsY3I+I76elfSpqSfj4FeCo93+l/izcDJ0n6BXA9SRPV54B9JNV2qsze05b7TT9/GcnWwp1mBBiJiMH0+CaSJNKtz/k4YG1ErI+IF4DvAG+i+59zTdHnWtnzdtJILAdmpiMvdifpUFvS4pgqIUkke62viYgrMh8tAWojKM4k6euonX9XOgpjDvBsrRrcCSLiooiYFhHTSZ7jv0XEnwF3AKekxervt/Z3OCUt33H/DzQingQek3RweupYYIgufc4kzVJzJP1++t947X67+jlnFH2utwFvlbRvWkt7a3quuFZ38LTLCzgR+E/g58D/bnU8Fd7XUSTV0JXAfenrRJL23NuBh9KfE9PyIhlJ9nNgFcnolJbfR8l7Pxq4NX1/IPBTYBj4FrBHen7P9Hg4/fzAVsc9jvs9DFiRPuvvAvt283MGPgo8ANwPXAPs0Y3PGbiOpN/mBZIaw9llnivw7vT+h4GzysbjGeFmZpabm6fMzCw3Jw0zM8vNScPMzHJz0jAzs9ycNMzMLDcnDTMzy81Jw8zMcnPSMCso3bfifZnjnzToe6ZJmt+Ia5uV5aRhVtw+wJakERFvatD3HEuyfpRZ23DSMCvuE8BBku6T9GlJ/wXJ0vPprnlfTneT+6ak4yTdme6WNrt2AUlnSPppeo1F6UZgZD4/CrgCOCUtM6Opd2g2Bi8jYlZQui/JrZHsGIek/4qIl6bnh0n2LFlNshDmz0jWCjqJZL2fkyUdAnwK+JOIeEHSF4G7IuLquu/5HvDXEXE/Zm1i150XMbMC1kbEKgBJq0m25AxJq0h2X4Ok2ekNwPJkgVZ+j61LW2cdTLLzmlnbcNIwq9b/y7x/MXP8Ilv/9ybgGxFx0VgXkfRykmWtX2hIlGYluU/DrLhNwN7j+P3bSfoq9gOQNFHSAXVlZtBZmyJZj3DSMCsoIn4N3Jl2dn+6xO8PAZcA35e0EvgBMKWu2APApPQ7GjU6y6wwd4SbmVlurmmYmVluThpmZpabk4aZmeXmpGFmZrk5aZiZWW5OGmZmlpuThpmZ5eakYWZmuf1/M6tvEbU8V+YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from sksurv.nonparametric import kaplan_meier_estimator\n",
    "\n",
    "time, survival_prob = kaplan_meier_estimator(data_y[\"Status\"], data_y[\"Survival_in_days\"])\n",
    "plt.step(time, survival_prob, where=\"post\")\n",
    "plt.ylabel(\"est. probability of survival $\\hat{S}(t)$\")\n",
    "plt.xlabel(\"time $t$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimated curve is a step function, with steps occuring at time points where one or more patients died. From the plot we can see that most patients died in the first 200 days, as indicated by the steep slope of the estimated survival function in the first 200 days."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Considering other variables by stratification\n",
    "\n",
    "### Survival functions by treatment\n",
    "\n",
    "Patients enrolled in the Veterans' Administration Lung Cancer Trial were randomized to one of two treatments: `standard` and a new `test` drug. Next, let's have a look at how many patients underwent the standard treatment and how many received the new drug."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "standard    69\n",
       "test        68\n",
       "Name: Treatment, dtype: int64"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data_x[\"Treatment\"].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Roughly half the patients received the alternative treatment.\n",
    "\n",
    "The obvious questions to ask is:\n",
    "> *Is there any difference in survival between the two treatment groups?*\n",
    "\n",
    "As a first attempt, we can estimate the survival function in both treatment groups separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f55a1c4d9e8>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEMCAYAAAA4S+qsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmcFPWd//HXB0FAGaMcI8qIoEEjiIwwMIPyy6pEBVzF+CCiaDQEr0Rdo0nWMXFF3UeyuGa94klQE41nNFGieIU1ByIoIisBNCCKDkJQUMADAf38/qjqoWl6Zrpqevp8Px+PfkxXdXXVp6aAD9/b3B0REZFMtMt3ACIiUjyUNEREJGNKGiIikjElDRERyZiShoiIZExJQ0REMqakISIiGVPSEBGRjClpiIhIxtrnO4Bs6969u/fp0yffYYiIFJVXXnnlA3fv0dJxJZc0+vTpw7x58/IdhohIUTGzFZkcp+opERHJmJKGiIhkLK9Jw8zuMrM1Zvb3Jj43M7vJzJaZ2WtmNjjXMYqIyDb5btP4NXAzcE8Tn48G+oWvWuC28KdI2dmyZQsNDQ1s2rQp36FIEevUqRNVVVV06NAh1vfzmjTc/a9m1qeZQ8YC93iw6MccM9vdzPZy91U5CVCkgDQ0NFBRUUGfPn0ws3yHI0XI3Vm7di0NDQ307ds31jkKvU2jF/Bu0nZDuG87ZnaOmc0zs3nvv/9+zoITyaVNmzbRrVs3JQyJzczo1q1bq0qrhZ400v3t2GGpQXef6u417l7To0eL3YxFipYShrRWa/8MFXrSaAD2SdquAt5rq4td9cdFXPXHRW11ehGRolfoSWM6cEbYi6oOWN+W7RmL39vA4vc2tNXpRYra2rVrqa6uprq6mp49e9KrV6/G7c2bN7fq3HfddRerV6/OUqT5uc6XX37JlClTsna+rVu3svvuu7fqHNOmTeMHP/hBliIK5LvL7QPAi8CBZtZgZpPM7DwzOy88ZAawHFgG/Ar4fp5CFSl73bp1Y8GCBSxYsIDzzjuPiy++uHF75513BoKG1i+//DLyuZU0Wi/u7z6qvCYNdz/V3fdy9w7uXuXud7r77e5+e/i5u/v57r6/uw909zadH+TM9bdz7Ls3cP/cd9ryMiIlZdmyZRx88MGcd955DB48mFWrVvHUU08xfPhwBg8ezPjx4/nkk08AmDx5MkOHDm083t156KGHWLBgAePHj28stVRVVfHTn/6Uuro6hg4dyvz58znmmGPYf//9+dWvftV47SlTpjBs2DAOOeQQrr766u3imTRpEgMGDGD06NFs2rQp7XXiWrlyJSNGjKC6upqDDz6Y2bNnU19fz8aNG6muruaMM84A4Pjjj2fIkCEMGDCAadOmAdtKEPX19QwaNIjhw4ezZs0aAN58801qa2sZOnQoV155ZeP1NmzYwFFHHcXgwYM55JBDeOKJJ5r83U+bNo0DDjiAI444gjlz5sS+x6ZY0Ju1dNTU1Hjcuaf+edNI3vrgE67vdT0PnTs8y5GJtM6SJUs46KCDgKD9LdtVqf333o3Jxw/I6Ngrr7ySLl268KMf/Yhly5ZxwAEHMHfuXIYOHcqaNWv41re+xVNPPcUuu+zCz372M8yMn/zkJ6xbt46uXbvi7kyYMIEzzjiD0aNHM2LECG6++Waqq6sBqKqqYvLkyZx99tlceOGFzJo1i1mzZvHxxx8zaNAgVq9ezYwZM3jiiSe45ZZbcHfGjBnDFVdcQWVlJQcddBDz589n4MCBnHTSSZx88smccsopO1wn2ZQpU3jwwQd32H/kkUdy/fXXb7fvmmuuAeDSSy/liy++4LPPPqNTp050796djz76qPG4xP1++umn1NTU8MILL1BRUUGHDh2YMWMGo0eP5pJLLqGyspL6+nrGjBnD6aefzoQJE7jxxhuZPHkyH330EVu2bGHTpk1UVFSwZs0aDj/8cJYuXbrD776hoYHDDz+c+fPnU1FRwde//nXq6uq44YYbtos/+c9Sgpm94u41LT37fA/uKyh7VnTig48/z3cYIkVn//33Z+jQoQDMnj2bxYsXc9hhhwGwefNmRowYAcDMmTO59tpr2bRpEx988AFDhgxh9OjRac95wgknADBw4EC2bt3Krrvuyq677kq7du34+OOPefbZZ3nqqac49NBDAfj444/5xz/+QWVlJV/96lcZOHAgAEOGDOHtt99u8R7q6+upr6/P6H6HDh3Kueeey6ZNmzjxxBMZNGgQW7du3eG466+/nunTpwPBOJs333yT6upqOnfu3HjfQ4YM4W9/+xsAL774In/84x8B+Pa3v83kyZOBoOrp0ksvZdasWbRr1453332XDz74ANj+dz9nzhxGjhxJt27dADj55JN5553s1pwoaaQYsHkh+73zO+6f24sJtb3zHY5IWpmWCHJl1113bXzv7owaNYp77713u2M+/fRTLrjgAubPn0+vXr24/PLLmx0v0LFjRwDatWvX+D6xvXXrVtydyy+/nEmTJm33vWXLlm13/E477ZT2H/RUUUoaRx11FH/+85958sknOe2007jssssYP378dsf86U9/4q9//Stz5syhc+fOjBgxovF+E21AqfGZWdousffccw/r169n/vz5tG/fnqqqqsZzJf/uE+doS4Xeeyq3Bo4DYOxOs3l8wco8ByNSnA477DD+8pe/sHz5cgA++eQTli5dymeffUa7du3o3r07Gzdu5NFHH238TkVFBRs3box0nWOPPZY777yzsb2koaGh8X/fTWnuOvX19Y0N+8mv1IQBsGLFCnr27Mk555zDd77zHV599VXatw/+D55IAOvXr6dr16507tyZRYsW8fLLL7d4T3V1dTz88MMA3HfffY37169fT2VlJe3bt+e5555j5cr0/z7V1dUxc+ZM1q1bx+bNm3nkkUdavGZUKmkkq5kICx+hYtX6fEciUrT23HNP7rzzTsaPH9/Y2Pzzn/+c4447jjPPPJODDz6Yfffdl9rabdPITZw4kbPOOovOnTvz0ksvZXSdMWPG8Prrr1NXVwcECeH+++9v9jup10n+H38UM2fO5LrrrqNDhw506dKF3/72twBMmjSJQw45hJqaGqZOncrUqVMZNGgQX/va17a736bcdNNNnHbaaVx33XV885vfbNz/7W9/m+OPP56amhoGDx5Mv3790n6/qqqKyy+/nLq6Ovbee29qalpsoohMDeGp7j6ORavWc3W3a9UYLgUlXeOlSBytaQhX9ZSIiGRMSUNERDKmpJHGgM0LGfnpjHyHISJScJQ0UoU9qA7/7Pk8ByIiUniUNFLVTGTRzgPzHYWISEFS0hARkYwpaYhIRjQ1+jbLly9PO3o8E/meDbe1lDREJCOaGn0bJQ3ZwYDNC5n7u//JdxgiBa8cp0avr6/n+eefp7q6mptuuomtW7dyySWXNMaSmAY90ynUi4mmEUnj437fhEUL6bL0D8AP8x2OyI6eqofVC7N7zp4DYXS8/wEvXryYu+++m9tvv501a9YwZcoUZs6c2Tg1+o033shPfvITLrroIq666qrGqdGffvppxo8fzy9/+csdpizv06cPc+bM4cILL2TSpEnbTY1+9tlnM2PGDN555x3mzp3bODX67Nmzqays5I033uCBBx5onBr9scce45RTTkl7nYQoExZOmTKFm2++mcceewyAW2+9lcrKSl566SU+//xz6urqOOaYY3jggQc4/vjjt5tCfdiwYUybNo0FCxbE+l3nm5JGGrXf+iGLlv4h32GIFI1ymxo91bPPPsuSJUsak8769etZunRpxlOoFxMljWYM2LwQ5t0dTGQoUkhilgjaSrlNjZ7K3bn11lsZOXLkDp+1NIV6sVGbRhNe6Hxk8GZh9qcWFill5TA1eup5jj32WG699dbG5PTGG2/w2WefZTSFerFRSaMJM3cZw+GfPU9hLXUjUvjKYWr0Qw89lC+++IJBgwYxadIkzj//fN55553GtpLKykoef/zxjKZQv+eee2LFkC+aGr0J4+94kSvW/pgBe30FJj6ZhchEWkdTo0u2aGp0ERHJiVYlDTPb1cx2ylYwIiJS2CIlDTNrZ2YTzOxJM1sDvA6sMrNFZnatmaVfg7CYrZgV9KASKQClVp0sudfaP0NRSxrPA/sDlwE93X0fd68E/h8wB5hiZqe3KqICoh5UUkg6derE2rVrlTgkNndn7dq1dOrUKfY5ovae+oa7bzGzfd29cYIZd18HPAo8amYdYkdTYH7+zzrG7j2bPfMdiAhQVVVFQ0MD77//fr5DkSLWqVMnqqqqYn8/UtJw9y3h2z8Ag5M/M7M6d5+TdExRG1vdi7lvreODjz9nz4r4WVkkWzp06EDfvn3zHYaUuahtGieb2RSgwswOSmkEn5rd0PJrQm1vavt2zXcYIiIFJWr11AtAJ+As4DrgQDP7CHgP+CzLsRWEjZu28s+Nm1RFJSJC9OqplcA9Zvamu78AYGZdgb4EPalKytjqXrCSoIoq38GIiBSAqNVTBpBIGOH7de7+irt/knxMKZhQ25uKTpppRUQkIXKXWzO70Mx6J+80s53N7Cgz+w1wZqYnM7NRZvaGmS0zsx3mJDaz3mb2vJm9amavmdmYiPGKiEgWRf1v9Cjgu8ADZtYX+AjoTJB8ngWud/eMVhYJG9FvAY4GGoCXzWy6uy9OOuxy4GF3v83M+gMzgD4RYxYRkSyJ2qaxCbgVuDUcj9Ed+MzdP4px7WHAMndfDmBmDwJjgeSk4cBu4fuvEDS4i4hInjRbPWVm/c3st0nbM81sADSO2RgKXGBmw2JcuxfwbtJ2Q7gv2ZXA6WbWQFDKuLCJOM8xs3lmNk8Dn0RE2k5LbRozCaqIEqrcfRGAmR0G3Av0Bn5tZt+MeO10Deap8yOcCvza3auAMcC9ZrZDzO4+1d1r3L2mR48eEcNoWeMKfiIiZa6lpHEM8LOk7Q1J788Abnf3c4AjgEsjXrsB2Cdpu4odq58mAQ8DuPuLBGNEuke8Tqto/ikRkW2aTRruvtDdT0vatczMxplZJXAi8Hh43BqgY7pzNONloJ+Z9TWznYFTgOkpx7wDjAQws4MIkkZO659m7jKGOV8exD83Nr2WsYhIuYja5fZi4FxgJTDf3WcDhI3iXaKcyN23AhcAzwBLCHpJLTKzq83shPCwHwJnm9n/AQ8A3/EcT/E5tjpoZtlz3TxVUYlI2Yvae2o1cLSZtUue5RY4kmDa9EjcfQZBA3fyviuS3i8GDo963myaUNubqbOPpm79kqCKqmZiPsMREcmrWCv3pSQM3P3ZsG2jJM3cZQyLdh6Y7zBERPJOa4SLiEjGlDRERCRjShoZSkyRLiJSziI1hJvZRnYcgAfBQD13993SfFb0NEW6iEggau+pirYKpJBNqO3NopmaIl1EJHb1lJntYWbDzOzriVc2AytEGzdt5f657+Q7DBGRvImVNMzsLOCvBAPzrgp/Xpm9sApP9y7BgPfHF6zMcyQiIvkTt6RxEcEMtyvc/UjgUHI8vUeu7VnRSav4iUjZi5s0NoVra2BmHd39deDA7IVVuOa+tU5VVCJStuImjQYz2x14DHjOzB6nDBZIUhWViJS7WPUt7p5YO+NKM3ueYFW9p7MWVYHas6ITtX275jsMEZG8idsQfrGZVQG4+1/cfbq7b85uaAVoxSxGfjpDVVQiUrbiVk/tBjxjZn8zs/PNrPTHvA0cB8DYnWYDqqISkfIUd5bbq9x9AHA+sDfwFzP7U1YjKzQ1E2HfEaqiEpGy1tq5p9YAq4G1QGXrwykei1dtYPwdL6qaSkTKSqyGcDP7HjAe6AE8ApwdLphUFhKr+S1eFSyZPqG2dz7DERHJmbij1fYFfuDuC7IZTLGYUNubCbW9GX/Hi/kORUQkp+K2adSXa8Jgxazt1gpPVFOpqkpEykGkpGFms8KfG81sQ9Jro5ltaJsQC0jYg4qFjwBBNVX/vYLZ4Bev2qAeVSJS8qJOjT4i/FmWU6RTM7ExYcC2aipAVVUiUhZaM7ivV7aDERGRwtaawX3PltXgPhER0eA+ERHJnAb3xZHSgypBc1KJSKmL26bxPTP7MzAT6E4wuO+QbAZWsFJ6UCUkBvypB5WIlLLIg/vMzIAaynVwX0oPqoQJtb2VMESk5EUuabi7A4eWZcJItnoh3H1c2moqEZFSFbdN40UzG5rVSIrJwHHQc2CQONKUOkRESlXcpHEkMMfM3jSz18xsoZm9ls3AClrNRJj4ZJA4RETKSNwJC0dnNYpiluhJVTMx35GIiLS5uCWNM5t4RWJmo8zsDTNbZmb1TRxzspktNrNFZnZ/zHjbRhM9qURESlXcpPFJ0usLgpJHnygnMLOdgFvC7/YHTjWz/inH9AMuAw4PBxP+IGa8bSNczS+ZxmqISCmLOyL8f5JePwOOAKLORTUMWObuy919M/AgMDblmLOBW9z9w/C6a+LEmysaqyEipa61I8ITdgH2i/idXsC7SdsN7Jh4DgAOMLMXzGyOmY1qRYxtbkJtb60fLiIlLe5yrwsBDzd3Ilj29eqop0mzz1O22wP9CEoyVcDfzOxgd/8oJZ5zgHMAevfO/9KriSoqLQMrIqUmbknjX4Hjw9cxwN7ufnPEczQA+yRtVwHvpTnmcXff4u5vAW8QJJHtuPtUd69x95oePXpEDCMLkuaiUhWViJSyuEljGLDO3VcAE4GHzWxwxHO8DPQzs75mtjNwCjA95ZjHCMaEYGbdCaqrlseMuW2k9KBSFZWIlLK4SeM/3H2jmY0AjgV+A9wW5QTuvhW4AHgGWAI87O6LzOxqMzshPOwZYK2ZLQaeB37s7mtjxtw20vSgAq0dLiKlKe7gvi/Cn8cBt7n742Z2ZdSTuPsMYEbKviuS3jtwSfgqGokqKgiSB6D2DREpCXGTxkozuwP4BnCNmXUkez2xip7WDheRUhU3aZwMjAJ+4e4fmdlewI+zF1ZpSVRVQVAKUalDRIpVrKTh7p8Cv0/aXgWsylZQRamJOahUVSUipURVStnQzBxUE2p789C5w3no3OH032u3HAcmIpJdcaunJFkTq/mlo6oqESlmkUoaZnZv+POitgmntI2t7tVY2li8aoMGAIpI0Yla0hhiZvsC3zWze0iZCsTd12UtshKkXlUiUuyitmncDjwNfA14JeU1L7uhlT5Noy4ixSZS0nD3m9z9IOAud9/P3fsmvaLOclvWNEeViBSjuF1uv2dmg4D/F+76q7uXzxrhWTChtrcShogUnVhdbs3s34D7gMrwdZ+ZXZjNwMrF3LfWaX4qESkacbvcngXUuvsnAGZ2DfAi8MtsBVYOElVUc98K+g+o+62IFLq4g/uMbZMWEr5Pt6iSNCMx8E9TqYtIsYibNO4G5prZleHstnOAO7MWVRlSTyoRKQaxkoa7X0ew+NI64ENgorvfkM3Ayol6UolIsYg9jYi7zwfmZzGWsqWeVCJSLDRhoYiIZEwTFmbT6oVw93HB+4HjdpgmXUSk2MUdp3GBme2R7WCK2sBx0HNg8H71woxnvRURKSZxSxo9gZfNbD5wF/BMuJ53+aqZuK1kkShtiIiUmLi9py4H+hF0s/0OsNTMfm5m+2cxtuKWqKpKvObd3eJXEmttqOutiBSq1vSecjNbDawGtgJ7AI+Y2XPu/u/ZCrAoJVbyS1i9MPjZTBtHotutloQVkUJmcWqVwrmnzgQ+AKYBj7n7FjNrByx197yVOGpqanzevAKbpf3u44LE0XNgiw3kiXU2Hjp3eK6iExHBzF5x95qWjotb0ugOnOTuK5J3uvuXZvavMc9ZuhIljwxKHLCtmkrLwYpIoYk7TqNjasIIJy3E3Ze0OqpSUzMRJj65rXdVMxJLwmo5WBEpRHGTxtFp9o1uTSASSEximFhLXESkkESqnjKz7wHfB/Yzs+RFlyqAF7IZmIiIFJ6obRr3A08B/wXUJ+3f6O7rshaViIgUpKhrhK9397fd/VR3X5H0UsLI1IpZGY3ZEBEpRJGShpnNCn9uNLMNSa+NZrahbUIsIYleVJpiRESKVKTqKXcfEf6saJtwSlzNRCUMESlqeZ0a3cxGmdkbZrbMzOqbOW6cmbmZtTjwpCioikpEilTU6qmNSdVRqa9I1VNmthNwC0FX3f7AqWbWP81xFcC/AXOjnL9gRaii0hKwIlJoojaEV7j7buHP1FfUgQXDgGXuvtzdNwMPAmPTHPefwH8DmyKevzDVTIR9R7R4mJaAFZFClK2G8A0xGsJ7Ae8mbTeE+5Kvdyiwj7s/EfHcha+FKqoJtb2p7ds1hwGJiLQsakmjsSE8LHFs94p4bUt3icYPg8kPrwd+2OKJzM4xs3lmNu/999+PGEYeRKyi0nTpIlIo8tkQ3gDsk7RdBbyXtF0BHAz82czeBuqA6ekaw919qrvXuHtNjx492jDkLIlQRVXbtytz31qnaioRKQhxl3vtZGaXmNnvzexRM7vYzDpFPM3LQD8z62tmOwOnANMTH4YDCbu7ex937wPMAU5w9wKb97ztJOahUjWViBSKuCWNe4ABwC+Bm4GDgHujnMDdtwIXAM8AS4CH3X2RmV1tZifEjKu4qOutiBSZuOtpHOjug5K2nzez/4t6EnefAcxI2XdFE8ceEfX8BW3guCBpLHykxfU1YFv3W62vISL5FLek8aqZ1SU2zKwWzXIbTYbtGqDutyJSOKJ2uV0YToleC8w2s7fDRuoXga+3QXylL4MqKnW/FZFCEbV6Sku5ZlPEKqrEMrCAloIVkbyIOmFh4xKvZrYH0A9I7jW1YocvSdMiTGCYqKKCIHkAShoiknOxGsLN7CzgIoKxFQsIxlC8CByVvdAk2YTa3o1JIlHaEBHJtbgN4RcBQ4EV7n4kcChQBEOxS4cmMxSRfIibNDa5+yYAM+vo7q8DB2YvLGmOelOJSL7EHafRYGa7A48Bz5nZh2w/BYhEsXoh3H3ctu2B45ptGJ9Q21sJQ0TyIlbScPdvhm+vNLPnga8AT2ctqnKSmLwwYfXC4GcGvalERHItbkN4J+D7wAiCmWlnkedVAItWzcTtE0RyiUNEpMDErZ66B9hIMPcUwKkEc099KxtBiYhIYcrr3FMiIlJcNPdUEUuMEFfXWxHJlUglDTNbSNCG0QE4w8wS/1r1BhZnOTZpRqLbrUaHi0guae6pIpUYIa7R4SKSS1HXCF+ReAG7A8eHr92T56WS3NLocBHJlbjLvV4E3AdUhq/fmtmF2QxMMqPR4SKSS3F7T00Cat39EwAzu4ZgwsJfNvstyUzqCPFkKaPFNTpcRHIpbu8pA75I2v4i3CetNXAc9ByY/rPVC5ucSl1VVCKSC3FLGncDc83sD+H2icCd2QmpzKWOEE/WROljbHUv5r61jscXrFQvKhFpU3HnnrrOzP5MMI2IARPd/dVsBiaZUxWViORK5OopC+zj7vPd/SZ3v1EJI4cS7R1p1hXXYD8RaWuRSxru7mb2GDCkDeKR5iRmxE0zE64G+4lILsRtCJ9jZkOzGom0rGYiTHwybUP5hNrePHTucPrvtVseAhORchG3IfxI4Dwzexv4hKBdw939kGwFJi1oolvuFWvX80LnI4HhuY9JREpe3KQxOqtRSDSpCzcl6bNleQ4DEZFyEzdp/JMdF2G6LVtBSQua6Zb79s9H5DgYESknWoSpBPXZsjzjEeUiIlFoEaYS80LnI/l08xfssmp9477uXTqyZ0UnrT8uIq0WN2m8amZ17j4HtAhTIely2Nn8YsGYxu3FqzbQv+NuPDRxuNYfF5FWi5s0atlxEaYliUWa1IsqfxLrbCRovQ0Ryaa4SWNUNi5uZqOAG4GdgGnuPiXl80uAs4CtwPvAd7VuR3SJkeJXrF1Pny3Lt2ss796lI3sedrqqrEQkI3Hnnmr1P9xmthNwC3A00AC8bGbT3T152dhXgRp3/9TMvgf8NzC+tdcuJ4mR4kA4fmObTzd/QZcPlwQz5yppiEgG4pY0smEYsMzdlwOY2YPAWJLWGnf355OOnwOcntMIS8D21VXbD/gLSh8/ZkDuwxKRIpXPpNELeDdpu4GgraQpk4Cn2jSictXcok8J6qorImQxaZhZT3dfHeUrafZ5E+c+HagB/qWJz88BzgHo3VsT9UXxyObhsDOQrotugrrqikgomyWNO4EofTobgH2StquA91IPMrNvAD8F/sXdP093InefCkwFqKmpSZt4ZEdjq3vxOCdxNSc17tuui26CuuqKSChrScPdo/7L8jLQz8z6AiuBU4AJyQeY2aHAHcAod1+TlUClUWr3XGimi25yFZaqqkTKVqykYWbXuPulLe1rjrtvNbMLgGcIutze5e6LzOxqYJ67TweuBboAvzMzgHfc/YQ4MUvmEl10Ey7pehS1PcMNVVWJlLW4JY2jgdQEMTrNvma5+wxgRsq+K5LefyNmfBJTchddCBLIdRzGQ+f+MNihqiqRshYpaYRjJb4P7G9mryV9VAHMzmZgkh8aUS4izYla0rifoNvrfwH1Sfs3uvu6rEUlhS2TLrqp1A4iUhIiJQ13Xw+sN7PfA+vcfaOZXQ4MNrP/dPdX2yRKyavkNo6Rnw7mcF8Pq9bv2DW3KWoHESkZcds0/sPdf2dmI4BjgV8At9P84DwpQqltHDN3GcPMXcak75rbFLWDiJSMuEnji/DnccBt7v64mV2ZnZCkkKTrlgtq6xApV3GTxkozu4OgF9U1ZtYRaJe9sKTkxGkHSaV2EZG8i5s0TiaYHv0X7v6Rme0F/Dh7YUkxSB3PMba6V9pSCQPHtf5iahcRKQhxk8ZnwK4Ea4NfDXQAPspWUFL40o3nANInjZqJrf/HXu0iIgUhbtK4FfgSOIogaWwEHgWGZikuKXB5Gc+RropLVVYiORV7uVd3H2xmrwK4+4dmtnMW45IilFpdBc1UWUWVropLVVYiORc3aWwJV95zADPrQVDykDKVWl0FLVRZRZWuiktVViI5Fzdp3AT8Aag0s58B44DLsxaVFJ1IM+aKSNGKu0b4fWb2CjCSYDGlE919SVYjE8lEc1151d4hknWx19Nw99eB17MYi5Sg5HaOrLVvJDTXlVftHSJtIp9rhEuJS27nyGr7RkJzXXnUPSS0AAAIhUlEQVTV3iHSJpQ0pM0kt3OofUOkNChpSM6k65KbLOvVV1qiViTrlDQkJ9J1yU2W9eqr5PYOtW+IZI2ShuREU7PlJmS9+iq5vUPtGyJZo6Qh5SEbs+w2R9VfUiaUNKRgtNTmEcfY6l5MyMYsu81R9ZeUESUNKQgttXnE0dhOcm4WZtltjqq/pIwoaUhBaKnNIw518xXJPiUNKWltOiI9WVu3mUSlNhZpI0oaUrLafER6Qlu3mUSlNhZpQ0oaUrJyNiI9GysTZlMhlXik5LTLdwAiIlI8VNKQstFUl942bevIl0JrY5G2k+P2KyUNKQtNdelt07aOfCm0NhZpO3lov1LSkLLQVJfekuyWW2htLNJ28lCaVNKQstcWI9HbSklWpUnr5Hg257w2hJvZKDN7w8yWmVl9ms87mtlD4edzzaxP7qOUUja2uhf999ot32FkZPGqDTy+YGW+w5BCMnAc9BwYvF+9EBY+0uaXzFtJw8x2Am4BjgYagJfNbLq7L046bBLwobt/1cxOAa4Bxuc+WilVbTESva0US2lIcigPsznns6QxDFjm7svdfTPwIDA25ZixwG/C948AI83MchijiIgkyWfS6AW8m7TdEO5Le4y7bwXWA91yEp2IiOwgnw3h6UoMHuMYzOwc4ByA3r2Lo6pBJKr+exdH24vkSaJto43lM2k0APskbVcB7zVxTIOZtQe+AqxLPZG7TwWmAtTU1OyQVERKweTjB+Q7BClko6fk5DL5rJ56GehnZn3NbGfgFGB6yjHTgTPD9+OA/3V3JQURkTzJW0nD3bea2QXAM8BOwF3uvsjMrgbmuft04E7gXjNbRlDCOCVf8YqISJ4H97n7DGBGyr4rkt5vAr6V67hERCQ9zXIrIiIZU9IQEZGMKWmIiEjGlDRERCRjShoiIpIxK7VhD2b2PrAi5te7Ax9kMZxioHsuD7rn8tCae97X3Xu0dFDJJY3WMLN57l6T7zhySfdcHnTP5SEX96zqKRERyZiShoiIZExJY3tT8x1AHuiey4PuuTy0+T2rTUNERDKmkoaIiGRMSSNkZqPM7A0zW2Zm9fmOJ1vMbB8ze97MlpjZIjO7KNzf1cyeM7Ol4c89wv1mZjeFv4fXzGxwfu8gHjPbycxeNbMnwu2+ZjY3vN+Hwun4MbOO4fay8PM++Yy7NcxsdzN7xMxeD5/38FJ+zmZ2cfhn+u9m9oCZdSrF52xmd5nZGjP7e9K+yM/VzM4Mj19qZmemu1YmlDQI/oEBbgFGA/2BU82sf36jypqtwA/d/SCgDjg/vLd6YKa79wNmhtsQ/A76ha9zgNtyH3JWXAQsSdq+Brg+vN8PgUnh/knAh+7+VeD68LhidSPwtLt/DRhEcP8l+ZzNrBfwb0CNux9MsLzCKZTmc/41MCplX6TnamZdgclALTAMmJxINJG5e9m/gOHAM0nblwGX5TuuNrrXx4GjgTeAvcJ9ewFvhO/vAE5NOr7xuGJ5EawCORM4CniCYNngD4D2qc+bYD2X4eH79uFxlu97iHHPuwFvpcZeqs8Z6AW8C3QNn9sTwLGl+pyBPsDf4z5X4FTgjqT92x0X5aWSRiDxBzChIdxXUsIi+aHAXGBPd18FEP6sDA8rhd/FDcC/A1+G292Aj9x9a7idfE+N9xt+vj48vtjsB7wP3B1Wy00zs10p0efs7iuBXwDvAKsIntsrlP5zToj6XLP2vJU0ApZmX0l1KzOzLsCjwA/cfUNzh6bZVzS/CzP7V2CNu7+SvDvNoZ7BZ8WkPTAYuM3dDwU+YVuVRTpFfd9h1cpYoC+wN7ArQdVMqlJ7zi1p6j6zdv9KGoEGYJ+k7SrgvTzFknVm1oEgYdzn7r8Pd//TzPYKP98LWBPuL/bfxeHACWb2NvAgQRXVDcDuZpZYqTL5nhrvN/z8KwRLCxebBqDB3eeG248QJJFSfc7fAN5y9/fdfQvwe+AwSv85J0R9rll73koagZeBfmHPi50JGtSm5zmmrDAzI1hrfYm7X5f00XQg0YPiTIK2jsT+M8JeGHXA+kQxuBi4+2XuXuXufQie4/+6+2nA88C48LDU+038HsaFxxfd/0DdfTXwrpkdGO4aCSymRJ8zQbVUnZntEv4ZT9xvST/nJFGf6zPAMWa2R1hKOybcF12+G3gK5QWMAf4BvAn8NN/xZPG+RhAUQ18DFoSvMQT1uTOBpeHPruHxRtCT7E1gIUHvlLzfR8x7PwJ4Iny/H/ASsAz4HdAx3N8p3F4Wfr5fvuNuxf1WA/PCZ/0YsEcpP2fgKuB14O/AvUDHUnzOwAME7TZbCEoMk+I8V+C74f0vAybGjUcjwkVEJGOqnhIRkYwpaYiISMaUNEREJGNKGiIikjElDRERyZiShoiIZExJQ0REMqakIRJRuG7F95O2Z7fRdarMbHxbnFskLiUNkeh2BxqThrsf1kbXGUkwf5RIwVDSEIluCrC/mS0ws2vN7GMIpp4PV82bFq4md5+ZfcPMXghXSxuWOIGZnW5mL4XnuCNcCIykz0cA1wHjwmP65vQORZqgaUREIgrXJXnCgxXjMLOP3b1LuH8ZwZoliwgmwvw/grmCTiCY7+dEMzsI+G/gJHffYma3AnPc/Z6U6zwN/Mjd/45IgWjf8iEiEsFb7r4QwMwWESzJ6Wa2kGD1NQiqnYYALwcTtNKZbVNbJzuQYOU1kYKhpCGSXZ8nvf8yaftLtv19M+A37n5ZUycxs24E01pvaZMoRWJSm4ZIdBuBilZ8fyZBW0UlgJl1NbN9U47pS3EtiiRlQklDJCJ3Xwu8EDZ2Xxvj+4uBy4Fnzew14Dlgr5TDXge6h9doq95ZIpGpIVxERDKmkoaIiGRMSUNERDKmpCEiIhlT0hARkYwpaYiISMaUNEREJGNKGiIikjElDRERydj/B70Gnu8IqWDgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for treatment_type in (\"standard\", \"test\"):\n",
    "    mask_treat = data_x[\"Treatment\"] == treatment_type\n",
    "    time_treatment, survival_prob_treatment = kaplan_meier_estimator(\n",
    "        data_y[\"Status\"][mask_treat],\n",
    "        data_y[\"Survival_in_days\"][mask_treat])\n",
    "    \n",
    "    plt.step(time_treatment, survival_prob_treatment, where=\"post\",\n",
    "             label=\"Treatment = %s\" % treatment_type)\n",
    "\n",
    "plt.ylabel(\"est. probability of survival $\\hat{S}(t)$\")\n",
    "plt.xlabel(\"time $t$\")\n",
    "plt.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately, the results are inconclusive, because the difference between the two estimated survival functions is too small to confidently argue that the drug affects survival or not.\n",
    "\n",
    "*Sidenote: Visually comparing estimated survival curves in order to assess whether there is a difference in survival between groups is usually not recommended, because it is highly subjective. Statistical tests such as the [log-rank test](https://en.wikipedia.org/wiki/Log-rank_test) are usually more appropriate.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Survival functions by cell type\n",
    "\n",
    "Next, let's have a look at the cell type, which has been recorded as well, and repeat the analysis from above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f55a1bb97f0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEMCAYAAAA4S+qsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl4VFW28OHfAgIRSJBRQJRBIUyBCGFQUUGRIURAIdCCA1HwqtehbeepRfvaarettu2ICrReVAQnPm7EFgUFGoGgyJggjRDCJAIyGiGwvj9OVVlJKklVpSpVlaz3eepJzlDnrJOCrOxz9l5bVBVjjDHGHzUiHYAxxpjYYUnDGGOM3yxpGGOM8ZslDWOMMX6zpGGMMcZvljSMMcb4zZKGMcYYv1nSMMYY4zdLGsYYY/xWK9IBhFqTJk20TZs2kQ7DGGNiysqVK39S1abl7VflkkabNm3Izs6OdBjGGBNTRGSrP/tF9PaUiEwVkR9FZG0p20VEnheRTSKyWkR6VHaMxhhjfhPpZxrTgSFlbB8KtHe9bgBeroSYjDHGlCKiSUNVvwL2lbHLCOBNdXwNnCoiLSonOmOMMcVF+zON04FtXsv5rnU7IxOOMbHp+PHj5OfnU1BQEOlQTITFx8fTqlUr4uLignp/tCcN8bGuxAQgInIDzu0rzjzzzHDHZEzMyc/PJyEhgTZt2iDi67+VqQ5Ulb1795Kfn0/btm2DOkakn2mUJx84w2u5FbCj+E6qOkVVU1U1tWnTcnuMGVPtFBQU0LhxY0sY1ZyI0Lhx4wq1OKM9acwBrnH1ouoLHFBVuzVlTBAsYRio+L+DiN6eEpF3gP5AExHJBx4B4gBU9RUgC0gDNgFHgcywBvTJfc7XoU+G9TTGGBOrIt176kpVbaGqcaraSlXfUNVXXAkDV6+p/1bVs1Q1WVXDO2pv1xrnZYwxAXjuued48803w3b8hx9+mG7dupGSksKgQYPYscO5S79w4UIaNGhASkoKKSkpPPbYYwAcO3aMCy+8kMLCwpDHEu23p4wxJqoVFhYydepUxo0bF7Zz3H333axevZpVq1aRnp7uSQ4AF1xwAatWrWLVqlX88Y9/BKB27dpccsklzJw5M+SxWNIwxoTdkSNHGDZsGN27d6dr166eX2bz5s2jY8eO9OvXj9tuu4309HQAJk+ezNNPP+15f9euXdmyZQsAI0eOpGfPnnTp0oUpU6Z49qlfvz733nsvPXv2ZODAgSxfvpz+/fvTrl075syZAzgdAjIzM0lOTuacc85hwYIFAEyfPp1bbrnFc6z09HQWLlzIiRMnmDBhAl27diU5OZlnn322xLV98cUX9OjRg1q1nLv9/fv3595776V379506NCBRYsWVfjnl5iYWORn6c9ziZEjRzJjxowKn7u4aO9yW6mekv3kFOwh7V93cEF2U37dkENiejoNx46JdGjGhMyj/28d63ccDOkxO7dM5JHLupS6fd68ebRs2ZL/+7//A+DAgQMUFBQwadIkvvjiC84++2zGjh3r17mmTp1Ko0aN+OWXX+jVqxejRo2icePGHDlyhP79+/PUU09x+eWX89BDD/HZZ5+xfv16rr32WoYPH86LL74IwJo1a8jJyWHQoEFs3Lix1HOtWrWK7du3s3atU+no559/LrHPkiVL6NmzZ5F1hYWFLF++nKysLB599FHmz59fZPuhQ4e44IILfJ7z7bffpnPnziXWP/jgg7z55ps0aNDAk+wAli5dSvfu3WnZsiVPP/00Xbo4n0PXrl1ZsWJFqdcWLGtpeGt8NtmnxJO1w/nL4OiKFRycOzfCQRkT+5KTk5k/fz733nsvixYtokGDBuTk5NC2bVvat2+PiHDVVVf5daznn3+e7t2707dvX7Zt28b3338POLdkhgwZ4jnfRRddRFxcHMnJyZ5WyuLFi7n66qsB6NixI61bty4zabRr147Nmzdz6623Mm/evCJ/8bvt3LmT4l39r7jiCgB69uzpObe3hIQEzy2l4i9fCQPg8ccfZ9u2bYwfP54XXngBgB49erB161a+++47br31VkaOHOnZv2bNmtSuXZtDhw6Ven3BsJaGl3sve5Oc6akANH/gAX7dkBPhiIwJvbJaBOHSoUMHVq5cSVZWFvfffz+DBg1i+PDhpd5mqVWrFidPnvQsu8cVLFy4kPnz57N06VLq1q1L//79Pdvi4uI8x6tRowZ16tTxfO9+IKxaYmxwmedr2LAh3333HZ9++ikvvvgi7733HlOnTi3y3lNOOaXEuAf3uWvWrOnzYXQwLQ23cePGMWzYMB599NEiSSwtLY2bb76Zn376iSZNmgDw66+/Eh8fX+qxgmFJwxgTdjt27KBRo0ZcddVV1K9fn+nTp3PPPffwww8/8J///IezzjqLd955x7N/mzZtmOtq5X/zzTf88MMPgHNbq2HDhtStW5ecnBy+/vrrgOK48MILmTFjBhdffDEbN24kLy+PpKQkDh48yEsvvcTJkyfZvn07y5cvB+Cnn36idu3ajBo1irPOOosJEyaUOGanTp3YtGlTQHG4Wxr++v7772nfvj0Ac+bMoWPHjgDs2rWL0047DRFh+fLlnDx5ksaNGwOwd+9emjZtGnS5kNJY0vDy6P9bx5Fjzl8Fby/L4/wIx2NMVbFmzRruvvtuatSoQVxcHC+//DLx8fFMmTKFYcOG0aRJE/r16+d5djBq1CjefPNNUlJS6NWrFx06dABgyJAhvPLKK3Tr1o2kpCT69u0bUBw333wzN954I8nJydSqVYvp06dTp04dzj//fNq2bUtycjJdu3alRw9nFobt27eTmZnpaYU88cQTJY45dOhQzy2vcLnvvvvIzc2lRo0atG7dmldeeQWA2bNn8/LLL1OrVi1OOeUU3n33XU9ra8GCBaSlpYU8FimtuRarUlNTNdhJmB79f+tYt3M8J08qNeUV/rLYqcTe+q3w9b82pjJs2LCBTp06RTqMMi1cuJCnn37a08KIJZdffjl/+ctfPK2BaHDFFVfwxBNPkJSUVGKbr38PIrJSVVPLO649CPfyyGVdqFe7FjVqWLkFY4z/nnzySXbujJ4KR8eOHWPkyJE+E0ZFWdLwYUOdE+yv+RXg9KDaP/O9CEdkTNXXv3//mGxlACQlJXHhhRdGOgyP2rVrc80114Tl2JY0iknTegAcqLmcRNdAI+t2a4wxDksaxWRQn9RfCmh4Yi8Nx46hbq9ekQ7JGGOihiWN4pJHA5B4suTIT2OMqe4saRSXmsmRGvUiHYUxxkQlSxrGmJhVv359ALZs2ULXrl2DOsaECROYPXs24DyML63L/ujRo9m8eXNwgQZg+PDhRa5l1apV9O3bl5SUFFJTUz0DD+fOncsjjzwS9niKs6RhjDHlWLduHSdOnKBdu3ZhPc8HH3zgSYRu99xzD4888girVq3iscce45577gFg2LBhzJkzh6NHj4Y1puIsafihICeHrVdf4/Nl3XGNKV9ppdHbtGnDAw88wLnnnktqairffPMNgwcP5qyzzvKMej58+DCXXHIJPXr0IDk5mY8//rjMc504cYK77rqL5ORkunXrxj/+8Q8AVq5cyUUXXUTPnj0ZPHhwQOMqZsyYwYgRIzzL9evX58EHH/QUTty9e3egP5ISDh8+zDPPPMNDDz1UZL2IcPCgU5X4wIEDtGzZ0rM+Et2UrYxIKU6cVN5elsdQV7dbXwpynIKGVjrdxJRP7gv9DJXNk8ucJtlXaXS3M844g6VLl3LHHXcwYcIElixZQkFBAV26dOHGG28kPj6eDz/8kMTERH766Sf69u1bZrHDKVOm8MMPP/Dtt99Sq1Yt9u3bx/Hjx7n11lv5+OOPadq0KTNnzuTBBx8sUXywNEuWLOHKK6/0LB85coS+ffvy+OOPc8899/Daa6+V+GW/YMEC7rjjjhLHqlu3Lv/+979LrH/44Ye58847qVu3bpH1zz33HIMHD+auu+7i5MmTRd6bmprKokWLGDOm8n4HWdLwIa5mDU6cPMHHq7Yz7r/GlJoUtl4dnsEzxlQ1ycnJ3HXXXdx7772kp6cXqfA6fPhwzz6HDx8mISGBhIQE4uPj+fnnn6lXrx4PPPAAX331FTVq1GD79u3s3r2b5s2b+zzX/PnzufHGGz2TIjVq1Ii1a9eydu1aLr30UsBpjbRo0cLv+IuXP69du7ZnwqiePXvy2WeflXjPgAED/C5KuGrVKjZt2sSzzz5bopT6yy+/zLPPPsuoUaN47733uP766z3zczRr1swz9WtlsaThQ+2aNah9/BDnHs0Czo10OMaEVhktgnDxVRrdPTWpdwlz9/fu5cLCQmbMmMGePXtYuXIlcXFxtGnTpkQpcm+qWqIVoqp06dKFpUuXBhV/8fLn3mXYSyt/HkhLY+nSpaxcuZI2bdpQWFjIjz/+SP/+/Vm4cCH//Oc/+fvf/w5ARkYGEydO9LyvoKCAU045JahrCpY90/ClnvMXxfm/LChnR2OMP3bs2EHdunW56qqruOuuu/jmm2/8fu+BAwdo1qwZcXFxLFiwgK1bt5a5/6BBg3jllVc8v8j37dtHUlISe/bs8SSN48ePs27dOr9jCKb8ubulUfzl69bUTTfdxI4dO9iyZQuLFy+mQ4cOLFy4EICWLVvy5ZdfAs7Ust5FETdu3Bh0r7FgWUvDl4TmHNkfPcXHjIl1vkqj+2v8+PFcdtllpKamkpKS4plLojQTJ05k48aNdOvWjbi4OCZNmsQtt9zC7Nmzue222zhw4ACFhYX8/ve/90yNWp5hw4axcOFCBg4c6HfcofLaa69x++23U1hY6Ckn77ZgwQKf5drDyUqj+5A5L5Mjed/yyE/16fLA4lL3cz/TsNLpJtrFQmn0aPbLL78wYMAAlixZQs2aNSMdDgC7d+9m3LhxfP755wG/10qjh8mhgkLeXpYX6TCMMRF2yimn8Oijj7J9+/ZIh+KRl5fH3/72t0o/b4WShojUE5HoSLshtqHOCb5KOM7Hq8r+R+Iew2HjNYyp2gYPHsyZZ54Z6TA8evXqRUpKSqWfN6CkISI1RGSciPyfiPwI5AA7RWSdiPxVRKJn2qoKSGvnTJGY3eBEmfslpqcT37EjBTk5Vj7dGFMtBNrSWACcBdwPNFfVM1S1GXAB8DXwpIhcFeIYK11GhwxStU65+zUcO4bWb71JfDkP5owxpqoItPfUQFU9LiKtVfWke6Wq7gPeB94XkbiQRmiMMSZqBNTSUNXjrm8/LL5NRPoW2yfm1Tt5hEuOZkU6DGOMiRqBPtMYIyJPAgki0qnYQ/Appb0vJrkG+J1zYD5jX13K2FeXWk8qY8Js+vTp3HLLLWE9h6py8cUXe4oAhtq2bdsYMGAAnTp1okuXLp7R3ABjx44lJSWFlJQU2rRp43mQvWbNGiZMmBCWeEIt0NtTS4B4YCLwDJAkIj8DO4BfQhxbZCU059jB3dSt7eTF9Tudf2Dj+kRP7wljTOCysrLo3r07iYmJYTl+rVq1+Nvf/kaPHj04dOgQPXv25NJLL6Vz586e6r4Ad955Jw0aNACculv5+fnk5eVFVQ8tXwK9PbVdVd8ERqjqUFVtBwwEHgEuDkeAkfRDzUKeblnA6AH5dG5R9j+woytWWLdbY8owcuRIevbsSZcuXYqMap42bRodOnTgoosuYsmSJZ71e/bsYdSoUfTq1YtevXp5tk2ePJnrrruO/v37065dO55//nnPe5555hm6du1K165dee6553zG4V3mfMuWLXTq1IlJkybRpUsXBg0axC+/VOzv3xYtWtCjRw8AEhIS6NSpU4nxHarKe++9V6Ry7mWXXca7775boXNXhoBaGiIi6vB8sq6H4PuK7+Pn8YYAfwdqAq+r6pPFtp8J/BM41bXPfapaKQ8Z0tqlwa415J78Bb6bCvyx1H0T09M5umIFB+fOtTLpJuo9tfwpcvblhPSYHRt15N7e95a5z9SpU2nUqBG//PILvXr1YtSoURw7doxHHnmElStX0qBBAwYMGMA555wDwO23384dd9xBv379yMvLY/DgwWzYsAGAnJwcFixYwKFDh0hKSuKmm25i9erVTJs2jWXLlqGq9OnTh4suushzPLclS5bw6quvepa///573nnnHV577TXGjBnD+++/z1VXFe0EOmPGDP7617+WuKazzz7bM+ufL1u2bOHbb7+lT58+RdYvWrSI0047rUgdqdTUVJ588knPJEvRKtDbUwtE5H3gY1X13OAXkdpAP+BanG6508s7kOt5yIvApUA+sEJE5qjqeq/dHgLeU9WXRaQzkAW0CTDmoGR0yCDj4GEys/8MR/aUuW/DsWNsnIYx5Xj++ef58EOnD822bdv4/vvv2bVrF/379/eUHR87diwbN24EnBLn69f/9uvg4MGDHDp0CHBqQdWpU4c6derQrFkzdu/ezeLFi7n88supV68eAFdccQWLFi0qkTT27dtHQkKCZ7lt27aeZws9e/YsUZocnPpX48ePD+h6Dx8+zKhRo3juuedK3Ap75513irQyIDJlzoMRaNIYAlwHvCMibYGfgVNwbnP9C3hWVf0rIA+9gU2quhlARN4FRgDeSUMB90+7Ac6zk8qTmglrX/QsLvthH28vy7PnGiamldciCIeFCxcyf/58li5dSt26denfv7+n1HhpkymdPHmSpUuX+iz97V1C3V2a3N86erVq1eLkyZPUqFHD57F83Z4KtKVx/PhxRo0axfjx47niiiuKbCssLOSDDz5g5cqVRdZHosx5MAJ9plGgqi+p6vlAa+AS4BxVba2qkwJIGACnA9u8lvNd67xNBq4SkXycVsatvg4kIjeISLaIZO/ZU3arIFgjUpzQyisrYowp6cCBAzRs2JC6deuSk5PD119/DUCfPn1YuHAhe/fu5fjx48yaNcvznkGDBvHCCy94lsub0OjCCy/ko48+4ujRoxw5coQPP/ywyGRPbklJSWzevDmg+MePH++zzLmvhKGqXH/99XTq1Ik//OEPJbbPnz+fjh070qpVqyLrI1HmPBhlJg0R6Swi/+u1/LmIdAHPeIxewC0i0juIc/v686L4nwpXAtNVtRWQBrwlIiViVtUpqpqqqqnes2uF0rg+Z9KnbaOwHNuYqm7IkCEUFhbSrVs3Hn74Yfr27Qs4D40nT57Mueeey8CBAz0PkMG5nZWdnU23bt3o3LmzZ87w0vTo0YMJEybQu3dv+vTpw8SJE0vcmoLfypyHy5IlS3jrrbf44osvPN1rs7J+exT77rvvlrg1BU6Z82HDhoUtrpBR1VJfwE6gjddyrtf35wEHcMZnrAcuL+tYPo59LvCp1/L9wP3F9lkHnOG1vBloVtZxe/bsqaE0YVpPnTDNOeaYV/6tY175t8/9tlx1tW656uqQntuYUFm/fn2kQ4gaO3bs0IEDB0Y6jCIKCgq0T58+evz48Uo5n69/D0C2+vG7u7zbU4OAx72WvUfDXAO8oqo3AP2BQG+UrgDai0hb14P03wFziu2Th3MLDBHphDNGJDz3n4wx1UKLFi2YNGlS2Ab3BSMvL48nn3zSM695NCszQlVdA3h3GdgkIqOBr4CRwBWu/X4UkfIr/BU9dqGI3AJ8itOddqqqrhORx3Ay3hzgTuA1EbkD59bVBFdGNMaYoI0ZE11d49u3b1+k+200CzSt3QG8BbwDfKaq/wZwFSmsH+jJ1RlzkVVs3R+9vl8PnB/occPJelAZY6qzQHtP7VLVS4E6qprmtWkAzviMKs16UBljqrugZu5Tr7LoruV/uZ5tVGnWg8oYU93ZHOHGGGP8ZknDGFMp6tcP+LFnhX377bdMnDgxbMefMWMG3bp1o1u3bpx33nl89913AOTm5nrGaKSkpJCYmOgpoHjXXXfxxRdfhC2mcIv+/l3GmGrHMyagRsX+rv3zn//MQw89FKKoSmrbti1ffvklDRs25JNPPuGGG25g2bJlJCUleUawnzhxgtNPP53LL78cgFtvvZVJkyZx8cWxWRg80Cq3hyg5ahuc0d2qquEpUB9huRwjc16mU/mWVqzfeZCxry4FnIfj7p5UBTk5bL36mqDOkZiebhVyTbVw+PBhRowYwf79+zl+/Dj/8z//w4gRI9iyZQtDhw5lwIABLF26lI8++oj58+fz1FNP0bJlS9q3b0+dOnV44YUX2LNnDzfeeCN5eU7d1Oeee47zzy/a0fLQoUOsXr2a7t27A05J9by8PDZv3kxeXh6///3vue222yp0Leedd57n+759+5Kfn19in88//5yzzjqL1q1bA9C6dWv27t3Lrl27aN68eYXOHwkBJQ1VTSh/r6olTeuBQO6+XABGpDzq2eY9MVNienrQ5yjIccpUW9IwlWHXn//MrxtCWxq9TqeONH/gAb/2jY+P58MPPyQxMZGffvqJvn37Mnz4cMC5rTNt2jReeuklduzYwZ/+9Ce++eYbEhISuPjiiz0JoKyy6W7Z2dklajn5KqkeFxdXZJ+xY8eSm5tbIu4//OEPXHNN6X8UvvHGGwwdOrTEel9lQ3r06MGSJUsYNWpUGT+p6BT07SkRaQi0xxmlDYCqfhWKoKJJBvXJ0PpkNmoGOAnC3bJwtzbA+YUf7C/9YFsnxsQiVeWBBx7gq6++okaNGmzfvp3du3cDzl/h7rpUy5cv56KLLqJRI6fHYkZGRrll071Lnu/cuZPiteh8lVQvXjjQe3Y9fy1YsIA33niDxYsXF1l/7Ngx5syZwxNPPFFkfayUQfclqKQhIhOB24FWwCqgL7CUKjh7HwBbF0P98yAh9pqSxhTnb4sgXGbMmMGePXtYuXIlcXFxtGnTxlMm3T0XBlBmqfOyyqa7nXLKKZ7juvkqqV5coC2N1atXM3HiRD755BMaN25cZNsnn3xCjx49OO2004qsj5Uy6L4E+5TpdpwKt1tVdQBwDlW1JlTyaOdrORMxGWP8c+DAAZo1a0ZcXBwLFixg69atPvfr3bs3X375Jfv376ewsJD333/fs82fsumdOnVi06ZNAcc3c+ZMn2XQfSWMvLw8rrjiCt566y06dOhQYruvyZYgdsqg+xJs0ihQ1QIAEamjqjlAUujCiiKpmdC6X6SjMKbKGD9+PNnZ2aSmpjJjxgw6duzoc7/TTz+dBx54gD59+jBw4EA6d+5MgwYNAP/Kpnfs2JEDBw54ZvsLh8cee4y9e/dy8803k5KSQmpqqmfb0aNH+eyzz0pMwnT8+HE2bdpUZN9YEuwzjXwRORX4CPhMRPZT2bPqRQmrRWWMfw4fPgxAkyZNWLp0qc991q5dW2R53Lhx3HDDDRQWFnL55ZczaNAgzzH8efZw3XXXMXPmTCZOnMjkyZPLPFcwXn/9dV5//XWf2+rWrcvevXtLrJ87dy6jR4+OiYq2vgRbRuRyVf1ZVScDDwNv4FS9rVasFpUx4TV58mRSUlLo2rUrbdu2ZeTIwH7N3HTTTUWeY0SDwsJC7rzzzkiHEbRgH4TfAcxS1XxV/TLEMUWt3H25nvEaGR0yGNfnzJAljKMrVrB/5nvW7dYYL08//XSF3h8fH8/VV18domhCIyMjI9IhVEiwzzQSgU9FZJGI/LeInFbuO2JcmtYjqVESuftyydqcVf4bAuAe43Fw7tyQHtcYbzYVjYGK/zsI9vbUo6raBfhvoCXwpYjMr1AkUS6D+kwbMo2kRqF/3t9w7Bjq9uoV8uMa4xYfH8/evXstcVRzqsrevXuJj48vf+dSVPRJzI/ALmAv0KyCxzLGhEmrVq3Iz89nzx7rOl7dxcfHlxjQGIhgn2ncBIwFmgKzgUmuWfaqJe9aVL5416cyJhLi4uJo27ZtpMMwVUCwLY3WwO9VteSImmrG3YOqNN71qYwxJtYFlTRU9b5QBxKrvGtR+VJWC8QYY2JNQA/CRWSx6+shETno9TokIgfDE2KU2LoYsqcBkL07m1kbZ0U4IGOMqXwBJQ1V7ef6mqCqiV6vhKo6lwbwW/2pNbNdc2oQ8m63xhgTC4Lqcisid4hI2TfzqxKv+lMZHTJIPS02a8YYY0xFBfsgPBH4l4jsA94FZqvq7tCFVbWUNtOfMcbEGhvcF2YjUk6ncwvnzt36nQetTpUxJqbZ4L4wK22mP2OMiUXBPtO4SUQWAp8DTXAG93ULZWDGGGOiT8AtDRERIBUb3GeMMdVOwC0NdSqenWMJI/QKcnLYevU17J/5XqRDMcYYn4J9prFURHqp6oqQRlMNlFanqmdiJ5IbHKT56rXsPvQrvW1eDWNMFAo2aQwAbhSRLcARQHAaIfZcowxl1ala2a0/K7v1J+N/H6fe4V8rMSpjjPFfsEljaChOLiJDgL8DNYHXVfVJH/uMASYDCnynquNCce5IKK9OFcAns2Jz3mBjTPUQ7G+oa0tZ/5i/BxCRmsCLwKVAPrBCROZ4l1gXkfbA/cD5qrpfRKxbrzHGRFCw070e8XqdwGl5tAnwGL2BTaq6WVWP4YwsH1Fsn0nAi6q6H0BVfwwyXmOMMSEQbGn0v3kvi8jTwJwAD3M6sM1rOR/oU2yfDq7jL8G5hTVZVecFeB5jjDEhEqob6HWBdgG+R3ysKz6BcS2gPdAfaAUsEpGuqvpzkQOJ3ADcAHDmmZVT1yl3Xy6Z8zIBSGuXRkaHjJAdu/mePLZefY3PbYnp6TS0nlXGmAgJdrrXNfz2C74mzrSvfj/PcMkHzvBabgXs8LHP16p6HPhBRHJxkkiRrr6qOgWYApCamlo88YScuzw6OMkDCFnSWJPkNLYa+dhWkJMDYEnDGBMxwbY00r2+LwR2q2phgMdYAbQXkbbAduB3QPGeUR8BVwLTRaQJzu2qzcGFHAKuiZgyUjM9ScLd2ggVd9fbmf91bsnTl9L6MMaYyhLsg/DewD5V3QpkAu+JSI9ADuBKMrcAnwIbgPdUdZ2IPCYiw127fQrsFZH1wALgblXdG2TMFeM1EZMxxlRXwbY0HlbVWSLSDxgMPA28TMkH2WVS1Swgq9i6P3p9r8AfXK/ISs20hGGMqfaCbWmccH0dBrysqh8DtUMTkjHGmGgVbNLYLiKvAmOALBGpU4FjGWOMiRHB/qIfg/O8YYir+2sj4O6QRWWMMSYqBTvd61FV/UDPkFWcAAAV50lEQVRVv3ct71TVf4U2tNiRvTubWRtnVcq53OXTrYS6MSYS7JZSBbnHbGRtzipnz4pLTE8nvmNHwEkeB+fODfs5jTHGm5VUraCMDhmVkjDAGdTnHthnYzaMMZEQUEtDRN5yfb09POHEANcAP2OMqY4CvT3VU0RaA9eJSEMRaeT9CkeAUcUG+BljqrlAb0+9AszDKU64kqJFB5XAixbGFhvgZ4yp5gJqaajq86raCZiqqu1Uta3Xq2onjErknkf87WV5kQ7FGGOKCHY+jZtEpDtwgWvVV6q6OnRhRblda2DasN+WZTfZ8iuz/nUHGYOerdCh3fOIr995EKDM6WHd3W/drGy6MSbcgupyKyK3ATOAZq7XDBG5NZSBRa3k0dA8uciqNK0HQNaORRU+/Lg+ZzLzv86lc4vEMvfz7n4L1gXXGFM5gu1yOxHoo6pHAETkKWAp8I9QBRa1UjOdl5cMIGt6aqWG4d39FqwLrjGmcgQ7uE/4rWghru99zcRnjDGmCgm2pTENWCYiH7qWRwJvhCYkY4wx0SrYB+HPiMhCoB9OCyNTVb8NZWDGGGOiT9BlRFT1G+CbEMZiinF3vQWnV1VZPamMMaYyWO2pKOXuegv+db81xpjKYFVuo5S7660/3W/d3OM2rGS6MSZcgh2ncYuINAx1MCZ47nEbNl7DGBNOwbY0mgMrROQ9ERkiItbdNsIajh1D67feLDLgzxhjQi3YmfseAtrjdLOdAHwvIn8WkbNCGJsxxpgoU5HeUyoiu4BdQCHQEJgtIp+p6j2hCjCmFBxw5tooNmI8FLx7UrlZjypjTGULuvaUiKwE/gIsAZJV9SagJzAqhPHFjnpNna9hKJ0+IuX0Eg/D1+88yMertof8XMYYU5ZgWxpNgCtUdav3SlU9KSLpFQ8rBiU0hyN7nFlFQmxcnzNLtCiKtzqMMaYyBJs06hRPGCLylKreq6obQhBXTMrlGJmyG+aVvD2V1i6NjA4ZEYjKGGNCJ9jeU5f6WDe0IoHEurR2aSRR2+e23H25ZG3OquSIjDEm9AJqaYjITcDNQDsR8Z50KQHn2Ua1ldEhg4wl052FIdOKbMv00fIwxphYFOjtqbeBT4AngPu81h9S1X0hi8oYY0xUCihpqOoB4ABwZXjCMYHw1Q13gqtO1T0BPCi3rrvGGH8Fentqsar2E5FDFO0nJDhDN/wrkmQqzLugYUVYMURjTCACbWn0c31NCE84xl++uuECbF3s5O2Z/3WuX8exrrvGmEBEtMqtq25VrohsEpH7ythvtIioiFTuRNzGGGOKCPT2lPu2lK8ChQHdnhKRmsCLON1383EKIM5R1fXF9ksAbgOWBRJrtMndl0vmvMxKGa/hLpHulpieTsOxY8J6TmNM9RDo7alQ3pbqDWxS1c0AIvIuMAJYX2y/P+GUK7krhOeuVGnt0gAncQBhTRqJ6UUH5Bfk5ABY0jDGhERAt6dEZLHr6yEROVj8FeC5Twe2eS3nu9Z5n+8c4AxVLXOCCBG5QUSyRSR7z549AYYRYlsXO0ULvWR0yGDakGkkNUoK++ndJdLdLyuVbowJpUg+CPd5i8uzUaQG8CxO6fXy4poCTAFITU0NQ/UnPyWPdpLGmtlhqXQbLr667nqzLrnGGLdIzhGeD5zhtdwK2OG1nAB0BRa65nhqDswRkeGqml1pUQYiNTMsVW7Dqbyuu9Yl1xjjLaikISLxOOVE+uG0DhYDL6tqQQCHWQG0F5G2wHbgd8A490bXQMImXudcCNwVtQkjRpXWddfNuuQaY7wF2+X2TaAL8A/gBaAT8FYgB1DVQuAW4FNgA/Ceqq4TkcdEZHiQcUUHH881jDGmKgj29lSSqnb3Wl4gIt8FehBVzQKyiq37Yyn79g/0+BFRznON7N3ZzNo4q1LLpBfvguvNuuMaYwIRbEvjWxHp614QkT5U8yq3HqmZ0Lqfz03urreVWSY9MT291B5UBTk5HJxbZsc0Y4wpItDBfWtwnmHEAdeISJ5r05mUHF9hisnokFHp82o0HDum1JZEaa2P4pb9sK/Isw3rTWVM9RXo7anqOZVrNVa8d5X1pjKmegt0nIZnilcRaQi0B+K9dtla4k0mphXvXWW9qYyp3oLtcjsRuB1nbMUqoC+wFLg4dKEZY4yJNsE+CL8d6AVsVdUBwDlAhOt3GGOMCbdgk0aBeyCfiNRR1Rwg/IWVjDHGRFSw4zTyReRU4CPgMxHZT9ESIGbXGpg2rOT6+ifJPX6g0sqkl8fXGA4bu2GMKU1QSUNVL3d9O1lEFgANgHkhiyrWJY/2vX7XGtKat4UWSZVSJr08xcuog3+l1MsrcOiLddM1pmoIZe2piM4CGFVSM31XuZ02jAwgY8g0MudFvgqurzEc5Y3dCGZucuuma0zVEeztqTeBQzi1pwCuxKk9Fdl7LSbsyitw6It10zWm6oho7alqyf2so/5JSGge6WiMMSYgVnuqMiWPhubJTuI4Yj2UjTGxx2pPVSb3s45pw4DdkY7GGGMCZrWnjDHG+C2g21OqutX9Ak4FLnO9TvWuS2X8cOyIc5vKJmsyxsSQYLvc3g5MAj5wrfpfEZmiqv8o423GLXk0rH3RSRylTNYUSeEY8OfP2A4by2FM9Au299T1QB9VPQIgIk/hFCy0pOGP1Ez46SunpaGRDqaoYAf8lcWfsR02lsOY2BBs0hDghNfyCdc6E+OCGfBXHn/GdthYDmNiQ7BJYxqwTEQ+dC2PBN4ITUjVzNbFznONKLtFZYwxvgQ1TkNVnwEygX3AfiBTVZ8LZWDVQr2mztc1syMbhzHG+CngloaICNBKVb8Bvgl9SNVIQnNo3S/SURhjjN8CbmmoquKURDfGGFPNBPtM42sR6aWqK0IajanWrOS6MdEv2KQxALhRRLYAR3B6TqmqdgtVYNVB9u5sZtGIDOpHOpRyhXuyJiu5bkxsCDZpDA1pFNVQWrs0sndnkyVHyNDoThrhGLtRnJVcNyY2BJs0dlNyEqaXQxVUdZDRIYOszVnOAL8oF46xG8aY2GSTMBljjPGbTcIUDdwTMxWXPNoG/RljoopNwhRp9Zo6EzMVt2uNDfozxkSdYFsafSg5CdMG9yRN1osqAAnNIcNHeXRfLQ9Tgr/ddK1rrjGhEWzSGBKKk4vIEODvQE3gdVV9stj2PwATgUJgD3Cdzdth3Pztpmtdc40JnaCSRih+cYtITeBF4FIgH1ghInNU1Xva2G+BVFU9KiI3AX8Bxlb03NEkd18umfN8PLeQ3aRpvajuWVB87EYox234w99uutY115jQCfaZRij0Bjap6mZVPQa8C4zw3kFVF6jqUdfi10CrSo4xrNLapZHUKMnntlyOkfXrjqid2S8xPZ34jh09ywU5ORycOzeCERljKkOwt6dC4XRgm9dyPs6zktJcD3wS1ogqWUaHDDI6+G5LZM4aCgWbonJmPyg5dsPGbRhTPYSspSEizQN9i491PuexE5GrgFTgr6Vsv0FEskUke8+ePQGGEaUSmkN8g0hHYYwxRYTy9lSgkzDlA2d4LbcCdhTfSUQGAg8Cw1X1V18HUtUpqpqqqqlNmzYNMAxjjDH+ClnSUNVA+4iuANqLSFsRqQ38DpjjvYOInAO8ipMwfgxNpMYYY4IV1DMNEXlKVe8tb11ZVLVQRG4BPsXpcjtVVdeJyGNAtqrOwbkdVR+Y5cz9RJ6qDg8m5pi1dbGNFg+BssZz2BgOY/wX7IPwS4HiCWKoj3VlUtUsIKvYuj96fT8wyPiqhnpNofVpJde7ixxa0vBLWeM5bAyHMYEJKGm4xkrcDJwlIqu9NiUA/w5lYAYbLR4iZY3nsDEcxgQm0JbG2zjdXp8A7vNaf0hV94UsKgOUHPiX1i6t1C660SDcEzUZYyIvoKShqgeAAyLyAbBPVQ+JyENADxH5k6p+G5Yoq6G0dmlFlnP35QJEbdKojImajDGRF+wzjYdVdZaI9AMGA08Dr1D24DwTgOID/3yWGokiNlGTMdVDsF1uT7i+DgNeVtWPgdqhCcn4xT0Hx7RhUVtqxBhT9QTb0tguIq/i9KJ6SkTqENk6VtVL8ujfvreeVBXmb3n1sli3XVNdBJs0xuCUR39aVX8WkRbA3aELy5QpNfO3JGE9qSrE3/LqZbFuu6Y6CTZp/ALUw5kb/DEgDvg5VEEZU1n8La9eFuu2a6qTYG8pvQT0xUkaAIdw5sYwxhhThQU93auq9hCRbwFUdb+rfpQJI+9xG9E+ZsPN19iNYNmYD2MiL9iWxnHXzHsKICJNgZMhi8qU4D1hU+6+XLI2e1Vf2bo4KntQFZ+oqSJskidjokOwLY3ngQ+BZiLyODAaeChkUZkSvMdtFBmzkTzaSRpROFmTr7EbwbIxH8ZEh2DnCJ8hIiuBS3AmUxqpqhtCGpnxT2qmkzBMRIWi224oWRdgEy5BT/eqqjlATghjMSYmhaLbbihZF2ATTpGcI9yYKiEU3XZDKZpaPKbqsVHcxhhj/GYtjarCXYsqmtlsg8bEPEsaMarIXBv1T0KLZsBuANK0HhnUj1xwvoSgRlagYz5sXIcxoWdJIwYVn2uDhObOC9e8G43akTEkysZtVLAV5Gu+jrLYXB7GhIcljRhUfK4Nb9E+70awAh3zUd3HdURbF2ATPpXdvdqShjFVTLR1ATbhE4nu1ZY0jKlioq0LsAmfSLQmLWmYylOZPbx27XC+Zk+zHlvGhJCN0zCVI3k0NE+u3HMeO2IlVowJMWtpmMrhPdtgZfjimt+6+RpjQsZaGsYYY/xmLY0qqMjAv2JiZfKmUCj48Rhb397htDqqCBuwaCLNkkYVU2Lgn5fcfbkA1SJpJKanV7nbUzZg0fjiPSanMsZsWNKoYqrjwD9fGo4dQ8Oj/3QSR/O9zsoYr31V3QcsmpK8x+RU1pgNSxqm6koe/dv3Iah9ZUy08R6TU1ljNixpmKrLu8dWtFcANiZGRLT3lIgMEZFcEdkkIvf52F5HRGa6ti8TkTaVH6Uxxhi3iCUNEakJvAgMBToDV4pI52K7XQ/sV9WzgWeBpyo3SmOMMd4i2dLoDWxS1c2qegx4FxhRbJ8RwD9d388GLhERqcQYjTHGeIlk0jgd2Oa1nO9a53MfVS0EDgCNKyU6Y4wxJUTyQbivFoMGsQ8icgNwA8CZZ1p1z9J0bNQx0iFETmXXvQqDOp2q8ednytW5ZWKlnEdUS/wOrpwTi5wLTFbVwa7l+wFU9QmvfT517bNURGoBu4CmWkbQqampmp2dHd7gjTGmihGRlaqaWt5+kbw9tQJoLyJtRaQ28DtgTrF95gDXur4fDXxRVsIwxhgTXhG7PaWqhSJyC/ApUBOYqqrrROQxIFtV5wBvAG+JyCZgH05iMcYYEyERHdynqllAVrF1f/T6vgCo+oWSjDEmRlhpdGOMMX6zpGGMMcZvljSMMcb4zZKGMcYYv1nSMMYY47eIDe4LFxHZA2wN8u1NgJ9CGE4ssGuuHuyaq4eKXHNrVW1a3k5VLmlUhIhk+zMisiqxa64e7Jqrh8q4Zrs9ZYwxxm+WNIwxxvjNkkZRUyIdQATYNVcPds3VQ9iv2Z5pGGOM8Zu1NIwxxvjNkoaLiAwRkVwR2SQi90U6nlARkTNEZIGIbBCRdSJyu2t9IxH5TES+d31t6FovIvK86+ewWkR6RPYKgiMiNUXkWxGZ61puKyLLXNc701WOHxGp41re5NreJpJxV4SInCois0Ukx/V5n1uVP2cRucP1b3qtiLwjIvFV8XMWkaki8qOIrPVaF/DnKiLXuvb/XkSu9XUuf1jSwPkFA7wIDAU6A1eKSOfIRhUyhcCdqtoJ6Av8t+va7gM+V9X2wOeuZXB+Bu1drxuAlys/5JC4HdjgtfwU8KzrevcD17vWXw/sV9WzgWdd+8WqvwPzVLUj0B3n+qvk5ywipwO3Aamq2hVneoXfUTU/5+nAkGLrAvpcRaQR8AjQB+gNPOJONAFT1Wr/As4FPvVavh+4P9JxhelaPwYuBXKBFq51LYBc1/evAld67e/ZL1ZeQCvXf6SLgbk40wb/BNQq/nnjzOdyruv7Wq79JNLXEMQ1JwI/FI+9qn7OwOnANqCR63ObCwyuqp8z0AZYG+znClwJvOq1vsh+gbyspeFw/wN0y3etq1JcTfJzgGXAaaq6E8D1tZlrt6rws3gOuAc46VpuDPysqoWuZe9r8lyva/sB1/6xph2wB5jmui33uojUo4p+zqq6HXgayAN24nxuK6n6n7NboJ9ryD5vSxoO8bGuSnUrE5H6wPvA71X1YFm7+lgXMz8LEUkHflTVld6rfeyqfmyLJbWAHsDLqnoOcITfbln4EtPX7bq1MgJoC7QE6uHcmimuqn3O5SntOkN2/ZY0HPnAGV7LrYAdEYol5EQkDidhzFDVD1yrd4tIC9f2FsCPrvWx/rM4HxguIluAd3FuUT0HnCoi7pkqva/Jc72u7Q1wphaONflAvqoucy3PxkkiVfVzHgj8oKp7VPU48AFwHlX/c3YL9HMN2edtScOxAmjv6nlRG+eB2pwIxxQSIiI4c61vUNVnvDbNAdw9KK7FedbhXn+NqxdGX+CAuxkcC1T1flVtpaptcD7HL1R1PLAAGO3arfj1un8Oo137x9xfoKq6C9gmIkmuVZcA66minzPObam+IlLX9W/cfb1V+nP2Eujn+ikwSEQaulppg1zrAhfpBzzR8gLSgI3Af4AHIx1PCK+rH04zdDWwyvVKw7mf+znwvetrI9f+gtOT7D/AGpzeKRG/jiCvvT8w1/V9O2A5sAmYBdRxrY93LW9ybW8X6bgrcL0pQLbrs/4IaFiVP2fgUSAHWAu8BdSpip8z8A7Oc5vjOC2G64P5XIHrXNe/CcgMNh4bEW6MMcZvdnvKGGOM3yxpGGOM8ZslDWOMMX6zpGGMMcZvljSMMcb4zZKGMcYYv1nSMMYY4zdLGsYEyDVvxc1ey/8O03laicjYcBzbmGBZ0jAmcKcCnqShqueF6TyX4NSPMiZqWNIwJnBPAmeJyCoR+auIHAan9Lxr1rzXXbPJzRCRgSKyxDVbWm/3AUTkKhFZ7jrGq66JwPDa3g94Bhjt2qdtpV6hMaWwMiLGBMg1L8lcdWaMQ0QOq2p91/pNOHOWrMMphPkdTq2g4Tj1fkaKSCfgL8AVqnpcRF4CvlbVN4udZx5wl6quxZgoUav8XYwxAfhBVdcAiMg6nCk5VUTW4My+Bs5tp57ACqdAK6fwW2lrb0k4M68ZEzUsaRgTWr96fX/Sa/kkv/1/E+Cfqnp/aQcRkcY4Za2PhyVKY4JkzzSMCdwhIKEC7/8c51lFMwARaSQirYvt05bYmhTJVBOWNIwJkKruBZa4Hnb/NYj3rwceAv4lIquBz4AWxXbLAZq4zhGu3lnGBMwehBtjjPGbtTSMMcb4zZKGMcYYv1nSMMYY4zdLGsYYY/xmScMYY4zfLGkYY4zxmyUNY4wxfrOkYYwxxm//H8hislNNj8yVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for value in data_x[\"Celltype\"].unique():\n",
    "    mask = data_x[\"Celltype\"] == value\n",
    "    time_cell, survival_prob_cell = kaplan_meier_estimator(data_y[\"Status\"][mask],\n",
    "                                                           data_y[\"Survival_in_days\"][mask])\n",
    "    plt.step(time_cell, survival_prob_cell, where=\"post\",\n",
    "             label=\"%s (n = %d)\" % (value, mask.sum()))\n",
    "\n",
    "plt.ylabel(\"est. probability of survival $\\hat{S}(t)$\")\n",
    "plt.xlabel(\"time $t$\")\n",
    "plt.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, we observe a pronounced difference between two groups. Patients with *squamous* or *large* cells seem to have a better prognosis compared to patients with *small* or *adeno* cells."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multivariate Survival Models\n",
    "\n",
    "In the Kaplan-Meier approach used above, we estimated multiple survival curves by dividing the dataset into smaller sub-groups according to a variable. If we want to consider more than 1 or 2 variables, this approach quickly becomes infeasible, because subgroups will get very small. Instead, we can use a linear model, [Cox's proportional hazard's model](https://en.wikipedia.org/wiki/Proportional_hazards_model), to estimate the impact each variable has on survival.\n",
    "\n",
    "First however, we need to convert the categorical variables in the data set into numeric values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Age_in_years</th>\n",
       "      <th>Celltype=large</th>\n",
       "      <th>Celltype=smallcell</th>\n",
       "      <th>Celltype=squamous</th>\n",
       "      <th>Karnofsky_score</th>\n",
       "      <th>Months_from_Diagnosis</th>\n",
       "      <th>Prior_therapy=yes</th>\n",
       "      <th>Treatment=test</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>69.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>64.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>70.0</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>38.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>3.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>63.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>60.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>65.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>70.0</td>\n",
       "      <td>11.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Age_in_years  Celltype=large  Celltype=smallcell  Celltype=squamous  \\\n",
       "0          69.0             0.0                 0.0                1.0   \n",
       "1          64.0             0.0                 0.0                1.0   \n",
       "2          38.0             0.0                 0.0                1.0   \n",
       "3          63.0             0.0                 0.0                1.0   \n",
       "4          65.0             0.0                 0.0                1.0   \n",
       "\n",
       "   Karnofsky_score  Months_from_Diagnosis  Prior_therapy=yes  Treatment=test  \n",
       "0             60.0                    7.0                0.0             0.0  \n",
       "1             70.0                    5.0                1.0             0.0  \n",
       "2             60.0                    3.0                0.0             0.0  \n",
       "3             60.0                    9.0                1.0             0.0  \n",
       "4             70.0                   11.0                1.0             0.0  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sksurv.preprocessing import OneHotEncoder\n",
    "\n",
    "data_x_numeric = OneHotEncoder().fit_transform(data_x)\n",
    "data_x_numeric.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Survival models in **scikit-survival** follow the same rules as estimators in scikit-learn, i.e., they have a `fit` method, which expects a data matrix and a structered array of survival times and binary event indicators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CoxPHSurvivalAnalysis(alpha=0, n_iter=100, tol=1e-09, verbose=0)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sksurv.linear_model import CoxPHSurvivalAnalysis\n",
    "\n",
    "estimator = CoxPHSurvivalAnalysis()\n",
    "estimator.fit(data_x_numeric, data_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is a vector of coefficients, one for each variable, where each value corresponds to the [log hazard ratio](https://en.wikipedia.org/wiki/Hazard_ratio)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Age_in_years            -0.008549\n",
       "Celltype=large          -0.788672\n",
       "Celltype=smallcell      -0.331813\n",
       "Celltype=squamous       -1.188299\n",
       "Karnofsky_score         -0.032622\n",
       "Months_from_Diagnosis   -0.000092\n",
       "Prior_therapy=yes        0.072327\n",
       "Treatment=test           0.289936\n",
       "dtype: float64"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.Series(estimator.coef_, index=data_x_numeric.columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the fitted model, we can predict a patient-specific survival function, by passing an appropriate data matrix to the estimator's `predict_survival_function` method .\n",
    "\n",
    "First, let's create a set of four synthetic patients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Age_in_years</th>\n",
       "      <th>Celltype=large</th>\n",
       "      <th>Celltype=smallcell</th>\n",
       "      <th>Celltype=squamous</th>\n",
       "      <th>Karnofsky_score</th>\n",
       "      <th>Months_from_Diagnosis</th>\n",
       "      <th>Prior_therapy=yes</th>\n",
       "      <th>Treatment=test</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>65</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>60</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>65</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>60</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>65</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>60</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>65</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>60</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Age_in_years  Celltype=large  Celltype=smallcell  Celltype=squamous  \\\n",
       "1            65               0                   0                  1   \n",
       "2            65               0                   0                  1   \n",
       "3            65               0                   1                  0   \n",
       "4            65               0                   1                  0   \n",
       "\n",
       "   Karnofsky_score  Months_from_Diagnosis  Prior_therapy=yes  Treatment=test  \n",
       "1               60                      1                  0               1  \n",
       "2               60                      1                  0               0  \n",
       "3               60                      1                  0               0  \n",
       "4               60                      1                  0               1  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_new = pd.DataFrame.from_dict({\n",
    "    1: [65, 0, 0, 1, 60, 1, 0, 1],\n",
    "    2: [65, 0, 0, 1, 60, 1, 0, 0],\n",
    "    3: [65, 0, 1, 0, 60, 1, 0, 0],\n",
    "    4: [65, 0, 1, 0, 60, 1, 0, 1]},\n",
    "     columns=data_x_numeric.columns, orient='index')\n",
    "x_new"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similar to `kaplan_meier_estimator`, the `predict_survival_function` method returns a sequence of step functions, which we can plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f55a0c87c18>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEMCAYAAAA4S+qsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xu8VXWd//HXR0TwAglyVOKuw1VBtKPg6BSVGiJKNcJRTJ2TiVlO5pQzluaFyUprLJus9JEiWpni6MiPSC1zplGBOBiBICgq6FFELioo4fXz+2Ptddxs9t5nr7XXvr+fj8d+nLPWXnvt7z7b+PT9fr7fz9fcHRERkULsVukGiIhI7VDQEBGRgiloiIhIwRQ0RESkYAoaIiJSMAUNEREpmIKGiIgUTEFDREQKpqAhIiIF273SDUhanz59fPDgwZVuhohITVmyZMkmd2/q7Lq6CxqDBw+mra2t0s0QEakpZraukOs0PCUiIgVT0BARkYIpaIiISMHqLqchIvXpnXfeob29nR07dlS6KTWte/fu9O/fn65du8Z6fUWDhpndAkwGXnH3Q7M8b8D1wCRgO/BP7v54eVspItWgvb2dHj16MHjwYIJ/GiQqd2fz5s20t7czZMiQWPeo9PDUrcDEPM+fCAxNPWYAPytDm0SkCu3YsYP99ttPAaMIZsZ+++1XVG+tokHD3f8EbMlzyRTgNg8sBPY1s77laZ2IVBsFjOIV+zesdE+jM/2AF9KO21PndmJmM8yszczaNm7cGP/dfndJ8BARkayqPWhkC4m7bGru7je5e7O7Nzc1dbqgMadL1z7EpWsfiv16EalvV199NYcccghjxoxh7NixLFq0qKTvN2HChEiLlefMmcMhhxzCbrvtVrJFztU+e6odGJB23B94qVRvdvCf/kZXf4dFff6DcVO/Vqq3EZEatGDBAubNm8fjjz9Ot27d2LRpE2+//Xalm7WTQw89lHvuuYfzzjuvZO9R7T2NucBZFhgPvO7u60v1Zgdu7sL+G41xK2ZC26xSvY2I1KD169fTp08funXrBkCfPn348Ic/DMDMmTM58sgjOfTQQ5kxYwbuwYDIhAkTuOiii/joRz/KyJEjWbx4MZ/97GcZOnQol112GQBr165lxIgRnH322YwZM4ZTTz2V7du37/L+Dz74IEcffTRHHHEEU6dO5Y033tjlmpEjRzJ8+PBS/QmAyk+5vQOYAPQxs3bgCqArgLv/HJhPMN12DcGU29ZStqdL125s423m9Nibjz72Sw5oLunbiUhMV/2/Fax8aWui9xz14Z5ccfIhOZ8/4YQTmDlzJsOGDeO4446jpaWFj33sYwBccMEFXH755QCceeaZzJs3j5NPPhmAPfbYgz/96U9cf/31TJkyhSVLltC7d28OPvhgLrroIgBWr17NzTffzDHHHMPnP/95fvrTn/L1r3+94703bdrEt7/9bf7whz+w9957c80113Ddddd1vGc5VXr21Onu3tfdu7p7f3e/2d1/ngoYpGZNfdndD3b30e5e0kqE++25HwC/2buJTW+8Vcq3EpEas88++7BkyRJuuukmmpqaaGlp4dZbbwXg4YcfZty4cYwePZo//vGPrFixouN1p5xyCgCjR4/mkEMOoW/fvnTr1o2DDjqIF14I5vkMGDCAY445BoDPfe5zPPLIIzu998KFC1m5ciXHHHMMY8eOZfbs2axbV1B9wcRVe06jrJr2bGLvJ9Zy1BPvc0if5cEQlXobIlUnX4+glLp06cKECROYMGECo0ePZvbs2Zx22ml86Utfoq2tjQEDBnDllVfutA4iHM7abbfdOn4Pj999911g12mwmcfuzvHHH88dd9xRqo9WsGrPaZRVz8mTATj4qXeZ02NvWH53hVskItVi9erVPP300x3HS5cuZdCgQR0Bok+fPrzxxhvcfXf0fzeef/55FixYAMAdd9zBscceu9Pz48eP59FHH2XNmjUAbN++naeeeiruRymKgkaaXi3TePPQwUAwRMW6R5QQFxEA3njjDc4++2xGjRrFmDFjWLlyJVdeeSX77rsv5557LqNHj+bTn/40Rx55ZOR7jxw5ktmzZzNmzBi2bNnC+eefv9PzTU1N3HrrrZx++umMGTOG8ePHs2rVql3uc++999K/f38WLFjASSedxKc+9anYnzcXC7P89aK5udmLmZ+87syzWL7xSb4/ZU/+d8NfYdCx0PrbBFsoInE8+eSTjBw5stLNSNzatWuZPHkyTzzxRNneM9vf0syWuHtzZ69VTyOLg9e+wbhVu7Nij9HqbYiIpFHQyBDmNY5e+Tq/fPOo4KRyGyJSIoMHDy5rL6NYChoZerVM4/mDe9C1y27c8d4ng96GiIgACho5de1ijBvSm2073mXDNm36IiICCho5bXtnGwMHLwPggC1tymuIiKCgkVW4MnyTL2RZr+ODk8priIgoaGTTtGcTPbr2AOChvSZpFpWIANVfGv3iiy9mxIgRjBkzhs985jO89tpribdJQaMAj+758eAX9TZEGlZ6afRly5bxhz/8gQEDBnT+wjI6/vjjeeKJJ1i2bBnDhg3ju9/9buLvoaBRgO9sGM+G3p2ueRGROlYLpdFPOOEEdt89KCk4fvx42tvbE/87qGBhDgOf2cZhj21g/5P6sei5LWx64y0O6NG90s0SEQi2ZX55ebL3PHA0nPi9nE/XWmn0W265hZaWlqL/LJnU08giXODXf+FzdO21iHFDegdPKK8h0rBqqTT61Vdfze67784ZZ5yR5J8AUE8jq14t03hmzizYupb5z84HZvDLN4/iu12XB3kNlUsXqaw8PYJSqoXS6LNnz2bevHk89NBDu9wnCepp5JA+g2rK2H4frA5Xb0OkIdVCafT777+fa665hrlz57LXXntFbkchFDQ60bahrWOISrOoRBpXLZRGv+CCC9i2bRvHH388Y8eO5Ytf/GLsz5uLSqPnsO7Ms9j4t420ntxO8wHNbF83A4A79/h20NuY/CMNU4mUkUqjJ0el0Uukac8mmg/I+BuOPjX4qd6GiDSgooKGme1tZl2Saky12b54MYc9tqHjeNFzW/j1e58MNmZ6eTnMOkn5DREpSl2XRjez3cxsupn91sxeAVYB681shZl938yGlqaZ5RdOux25ZBMQJMMB7lv6YtDbODCVFJ/3VQUOEWkYUXsaDwMHA98ADnT3Ae6+P/APwELge2b2uYTbWBG9WqaxV1pCa/q4gR+s12huDbaAnfyj4FhDVSLSIKIGjePc/d+B1939/fCku29x9/9y938E7ky0hVVm0XNb+PWi54OD5tZgqErTcEWkQUQKGu7+TurXezOfM7PxGdfUhbCcCGQMUYWUGBeRBhI1pzHNzL4H9DCzkRlJ8JuSbVrlZeY1dhqiCqm3IdIwqr00+re+9a2Otp1wwgm89NJLibcp6vDUo8BKoBdwHfC0mT1uZvOAvyXduEoL9wvf9s425jw1p+P8TkNUoN6GSAOohdLoF198McuWLWPp0qVMnjyZmTNnJv4eUYenXnT324Ap7n6iux8EHAdcAXwi8dZVgXAXv6AGVY4hqrC3ISJ1qxZKo/fs2bPj9zfffLMktaciFSw0M/PAo+E5d98CbMm8JsE2VlTTnk1s/tvmjuPp4wZy39IXWbl+Ky03LmDK2H5MHzewgi0UaTzX/PkaVm3ZtYxGMUb0HsG/HfVvOZ+vldLol156Kbfddhsf+tCHePjhhxP7+4QiT7k1s382s53+lTSzPczsE2Y2Gzg7ueZVpylj+zGqb09Wrt+6c49DROpWrZRGv/rqq3nhhRc444wz+MlPfpL0nyFyafSJwOeBO8xsCPAasCdB8HkQ+KG7Ly30ZmY2Ebge6AL8wt2/l/H8QGA2sG/qmkvcfX7ENieibUMbc56aw9RhU5k+biDTxw2k5cYFHfmN6ZVolEiDytcjKKVaKI0emj59OieddBJXXXVV5M+ZT9Scxg53/6m7HwMMAj4JHO7ug9z93IgBowtwA3AiMAo43cxGZVx2GXCXux8OnAb8NEp7k5KZ1whlzW+ISF2qhdLo6e2bO3cuI0aMiNyWzuQNGmY2ysx+mXb8kJkdAh3rMY4ELjCzo2K891HAGnd/1t3fBn4DTMm4xoEws/MhIPn5YwXIWriQD6bgLnpuCxu27dC0W5E6Vgul0S+55BIOPfRQxowZw4MPPsj1118f+/Pmkrc0upmtB45297Wp49XuPjz1+98DvyNYAX4scKm777LoL8+9TwUmuvsXUsdnAuPc/YK0a/oSDHv1AvYmWJG+JN99kyqNHlp35lkAXHlGsCRl1sSdg8KvFz3PN+9dzjcPWMiM138czKJq/W1i7y8iAZVGT04pS6OfAFyddrw17fezgJ+7+wxgAhB1kDHbXLDMCHY6cKu79wcmAbeb2S5tNrMZZtZmZm0bN26M2IzihL2N72wYz4benf69RURqWt6g4e7L3T19Z/I1Znaqme0PfBq4L3XdK0C3bPfIox1IXxnTn12Hn84B7kq9xwKgO9AnSztvcvdmd29uamqK2IzOZZZIzxTmNja98ZaGqEQkkroujQ5cBJwHvAg87u6PAZhZV2CfiPdaDAw1syFmtgdBontuxjXPEyTbMbORBEGjrF2JzFIi2YS9jY7tYFUuXUTqVNTZUy+7+/FAN3eflPbUxwnKpke517vABcADwJMEs6RWmNlMMzslddnXgHPN7K/AHcA/lXvhYGaJ9Hwe2muSyqWLSF2Luk4DgPSy6KnjBwkS1lHvMx+Yn3Hu8rTfVwLHxGljKaSv1cipuTUIGOHOfhDUptJ+4iJSB7RHeIFyrdXIKtzZD7S7n4jUFQWNAu39xFrOeabAipbhzn7a3U+krlR7afTQD37wA8yMTZty52LjUtAoQCHJ8JzCCrjhcJV6HCI1qRZKowO88MIL/P73v2fgwNIUUo26CdM2M9ua5bHNzLZ2fofaFCUZnlU4XPXycvU4RGpULZRGB7jooou49tprS1IWHSImwt29R0laUUPyJcM7ihdmlkpvbg0eYWJcRIry8ne+w1tPJlsavdvIERz4zW/mfL4WSqPPnTuXfv36cdhhhyX6t0kXe3jKzHqZ2VFm9tHwkWTDqlG+ZHjBxQu1+E+kJlV7afTt27dz9dVXl2S3vnSxptya2ReACwlWcS8FxgMLqNPd+0JB4cIDsz4Xbs6U1+hTg6Cx/G5NwRUpQr4eQSlVc2n0Z555hueee66jl9He3s4RRxzBn//8Zw48MPu/W3HE7WlcSFDhdp27fxw4nDKv1K5JYVJcvQ2RmlPtpdFHjx7NK6+8wtq1a1m7di39+/fn8ccfTzRgQPygscPddwCYWTd3XwUMT65Z1S3Ma2QT5jVyGn1q8FNrN0RqSi2URi+HvKXRc77I7F6gFfgqwZDUq0DXjNIiFZF0afRQWCL9z986mZkLZtJ8QHPOMukA3/nM6Nx7h7fNCoKGyqiLFEyl0ZNTytLoWbn7Z9z9NXe/EvgWcDNB1du6N3XY1KwbMkGQ1/jOZ4KV4HnzG+EwlYhIjYkVNMzsIjPrD+Du/+vuc1O77zWM1VtW03p/6y7DVGHF24IotyHS8Oq9NHqoJ/CAmf2fmX3ZzA5IslHVavvixbx6511MOmgSw3sPZ/WW1TlrUa1cv5WWGxfQcuOC7DmOMLehxX4iBStzkeu6VOzfMO7w1FXufgjwZeDDwP+a2R+KakmVC0uJbJ03j6nDpjJr4iyG986e+58yth+j+gZbmy96bgvfvHf5roFDQ1QikXTv3p3NmzcrcBTB3dm8eTPdu3ePfY9Y6zTSvAK8DGwG9i/yXlWtV8s0ts6bt8v5bCvEp48b2JEED5Pj9y19MXtiPByi0roNkbz69+9Pe3s75d7Sud50796d/v37x3593MV95wMtQBNwN3Buau+LhjLpoEm0bWhj/rPzc+6xES76C4erpozt90Hw0GI/kYJ17dqVIUOGVLoZDS9uT2MQ8FV3X5pkY2rN1GFTC9pfIywxsnJ9UNOxI2iEGzaJiNSIuDmNSxo1YITJ8HT5FvtBECTuPO/ojjzHLjSLSkRqRNTS6I+kfmaWSK/r0uih9GR4aNJBwXrGgnb0y0YrxEWkhkQKGu5+bOpnD3fvmfbo4e45/m90/ci2r0a+xX7Z7FJmpLlVu/uJSM0oZnFfv6QbU+9ylk/X9FsRqRHFLO57sNEW9+XTWV4DPlgt3mlRQxGRKqXFfQmIktfIu1mT9hEXkSoXe+e+lIZZ3JcucwZVlLxGzt6G9hEXkRoQN6dxvpn9D/AQ0Idgcd+YJBtWrbLNoIoqa2+juTUok37g6KLaJyJSSpEX91mwD2EzDbq4L1c5kSgK2hpWRKQKRe5peFAt7PBGDBjpdqxaxbozz9ploV8USoiLSK2Jm9NYYGbR9zSsEz0nT6b7iBHsWLUqdq8jb0JcRKRKxQ0aHwcWmtkzZrbMzJab2bIkG1bNerVMY9Dtt9F9xIjY98i7WZNmUYlIlYpbsPDERFvRwMIhqp0q30IQOEDVb0WkqsTtaZyd49FwshUwLJRmUYlIrYkbNN5Me7xH0PMYHPUmZjbRzFab2RozuyTHNdPMbKWZrTCzX8dsb0kUO/020n7iIiJVINbwlLv/R/qxmf0AmBvlHmbWBbgBOB5oBxab2dz0zZzMbCjwDeAYd3/VzKpqAWES029FRGpJsSvCQ3sBB0V8zVHAGnd/1t3fBn4DTMm45lzgBnd/FcDdXym6pSVUSP2pbHJOvdU+GyJSZeKuCF+emjW1zMxWAKuB6yPeph/wQtpxe+pcumHAMDN71MwWmtnEOO0tte2LF3PGU0EnKOq+Gjmn3oYJcZUUEZEqErenMRk4OfU4Afiwu/8k4j0syznPON4dGApMAE4HfmFm++5yI7MZZtZmZm3l3nQ+zGsMXbw+0r4aoZy1qMJy6Zp+KyJVJG7QOArY4u7rgFbgLjM7IuI92oEBacf9gZeyXHOfu7/j7s8R9GiGZt7I3W9y92Z3b25qaorYjOJk25gpqry9DRUxFJEqEjdofMvdt5nZscCngNnAzyLeYzEw1MyGmNkewGnsmkz/b4KFhJhZH4Lhqmdjtrmkti9ezGGPbYj12pyzqNKn3yq/ISJVIG7QeC/18yTgZ+5+H7BHlBu4+7vABcADwJPAXe6+wsxmmtkpqcseADab2UrgYeBid98cs80lEw5RjVyyqTRvoPyGiFSJuEHjRTO7EZgGzDezbnHu5e7z3X2Yux/s7lenzl3u7nNTv7u7/4u7j3L30e7+m5jtLan0Iaq4M6gAVq7fSsuNC3adSRXmN9TbEJEKixs0phH0Aia6+2tAb+DixFpVo/bbcz8g+gwqCPIao/r2ZOX6rdmLGKq3ISJVIO52r9vd/R53fzp1vN7dH0y2abWnac8mmg9ojtXbmD5uIHeedzSj+vbMfoFmU4lIFUhqcZ+kRNkvPDLNphKRClPQSFiU/cIjUzFDEamwSEHDzG5P/bywNM2pH6u3rKb1/lZa72+NnRjPS0lxEamAqD2Nj5jZIODzZtbLzHqnP0rRwFo06aBJDO89HAiCR+JDVUqKi0iFRA0aPwfuB0YASzIebck2rXZNHTaVWRNnMWviLIb3Hh45MZ5z6m0oTIqLiJRZpKDh7j9295HALe5+kLsPSXtErXLbEKImxjudeisiUkFxp9yeb2aHmdkFqceYpBtWL8LEeJjj6KzH0enU23SafisiZRa3NPpXgF8B+6cevzKzf06yYfUkzHEkmt/Q9FsRqYC4U26/AIxLlfy4HBhPsGGSZBHmOMLkeCJUzFBEKiBu0DA+KFpI6vds+2NIqWkmlYiUUdygMQtYZGZXmtmVwELg5sRaVaO2L17Mq3feldj9cm4Dm07FDEWkjOImwq8j2HxpC/Aq0OruP0qyYbUmLI++dd68RO6Xc2OmbNTbEJEyiV1GxN0fT03Bvd7d/5Jko2pRoTv4FbpmI9yYqdM1G6B1GyJSNqo9VUYlX7OhKbgiUmIKGmUUrtmI0tsI12x0mt/QFFwRKYO46zQuMLNeSTemEcQpnV5QfkMVcEWkDOL2NA4EFpvZXWY20cw03TalsxlUUXsb8EF+o6DZVCIiJRR39tRlwFCCabb/BDxtZt8xs4MTbFvNKXQGVcl6GyIiJVbM7CkHXk493gV6AXeb2bUJta3mhDOoCu1tRBFpNpWISInErj1lZkuAa4FHgdHufj7wEeAfE2xfzUl6vUa6gmdThbOowodmU4lIQnaP+bo+wGfdfV36SXd/38wmF9+s2tWrZVpJAgYEvY3p4wbScuOC3BeFC/1CLy8Pfja3lqRNItJY4g5PdcsMGGZ2DYC7P1l0qxpE1M2ZChLOogof4TRc9TpEJAFxg8bxWc6dWExD6s2OVatYd+ZZrDvzrKz5jTjJ8HQFz6QK12+A1nCISNEiBQ0zO9/MlgPDzWxZ2uM5YFlpmlh7ek6eTPcRI4AgeGQbroqTDA9FmkmV3vPQGg4RKVLUnMavgd8B3wUuSTu/zd23JNaqGterZRq9WqYBsO7MszpmU4XnijV93MD4U2/DoarQ6FOV7xCRgkXdI/x1d1/r7qe7+7q0hwJGDp3NpipJXiOX9KEq0HCViEQWdXjqkdTPbWa2Ne2xzcy2lqaJtS1cuxHmONLzG8XmNSLLliTXPhwiEkHUnsaxqZ893L1n2qOHu/csTRNrX5jjyMxvFJPXAIpf6Kd9OEQkorjrNCSCMMex7syzErtnmAxfuT7o4E0fNzD6TZpbg4ChPIeIFCjq8NS2tOGozEfk4alUscPVZrbGzC7Jc92pZuZmFv//llexOHmNSGXT81GeQ0QiiDo81SNtOCrzEWl4ysy6ADcQrO8YBZxuZqOyXNcD+AqwKMr9a0WxeY2iCxlmy3OIiOSQVCJ8a4yexlHAGnd/1t3fBn4DTMly3b8T1LjaEfH+NSFOqfR0JSmbrtpVIpJDUonwnjES4f2AF9KO21PnOpjZ4cAAdy9NMacqUfHeRjoNV4lIHpXc7jXbxk3e8aTZbsAPga91eiOzGWbWZmZtGzduTLCJyctWNr3YWVSJlk3XtFwRySNuafTuZvYvZnaPmf2XmV1kZt0j3qYdGJB23B94Ke24B3Ao8D9mthYYD8zNlgx395vcvdndm5uamiI2o3yqomx6VJqWKyJp4vY0bgMOAf4T+AkwErg94j0WA0PNbIiZ7QGcBswNn0ytPu/j7oPdfTCwEDjF3dtitrniCt2kKY702VSJam6FQccme08RqVlxg8Zwdz/H3R9OPWYAw6LcwN3fBS4AHgCeBO5y9xVmNtPMTonZrqqXr7exestqWu9vLV9ZkShUXl1EiL+47y9mNt7dFwKY2TiCHfwicff5wPyMc5fnuHZCjHZWnXCTpswihmEyfPWW1UCQ56ga6Rs7aVMnkYYWdcrtcjNbBowDHjOztal8wwLgoyVoX13K1tuYOmwqsybOYnjv4ZVqVm4qry4iKVF7Gg29lWtSSrklbFmEQ1UqNyLScCIFjfQtXs2sFzAUSJ81tW6XF0nZhVNvQ1PG9otXmyqbcKhKw1QiDSlWTsPMvgBcSDBNdinBdNgFwCeSa5rEES70CxVV0DCb5tbgMeskFToUaUBxE+EXAkcCC93942Y2ArgquWY1hqR39IMgOKQHiJYbF+zU80is15GeHAf1PEQaRNygscPdd5gZZtbN3VeZWRVmcKtXz8mT2b54MVvnzUs0aGRK73kk2usIexyh9B6HiNStuEGj3cz2Bf4b+L2ZvcrOq7mlE+VKhqf3PNLzHCWROVwFGrISqTOxgoa7fyb165Vm9jDwIeD+xFrV4MJFfpMOmpT4eo2SJckzh6tAQ1YidShuIrw78CXgWIIig49Q2eKHNauci/xKmiTPHK4CDVmJ1KFK1p5qeOVe5BfWpwofRe/6V4hwyEqlR0TqQsVqT8kHBQwrJdF9OLIJ9+bQnhwidSNu0PiLmY0PD+LWnpLKCvfhKJmw/IhKj4jUjUg5DTNbTpDD6AqcZWbhuMZAYGXCbWsYpVivEUVJV5CHwo2clBQXqWlRexqTgZOBicAQ4GOpxxBAWc8Y8pVKj7tveBTh5k2hxDdxAm3kJFJHiqk9dRjwD6nD/3P3vybZsEaRa73GpIMm0bahjfnPzi9pmfTOVpBDAj2P5lYFDJE6EXfK7YXAucA9qVO/NLOb3P0/E2tZg8kcopo6bCrzn51f0jUb2ZR0Wm62xX/ptBBQpOrFXRF+DjDO3d8EMLNrCAoWKmjEkKukSCU2ZsrW8win5RYVOLIt/kunhYAiNSFu0DDgvbTj91LnJIZcQ1RTh01l6rCptN7f2pHfKPeOflPG9mPRc1u4b+mLxQ9R5QsIWggoUhPiBo1ZwCIzuzd1/Gng5mSaJJnKld/IZvq4gaVbx5FJw1ciVS/WOg13vw5oBbYArwKt7v6jJBvWiHasWsW6M8/i1Tvv2un81GFTaT6guSyzqXIJk+MtNy4ozQrycCFgLlogKFIVIvc0zMyA/u7+OPB48k1qTOHU2x2rVgHssmajkr2NkpVXT6fhK5GaYO4e/UVmS9z9IyVoT9Gam5u9ra2t0s2Ibd2ZZwEw6Pbbdnmu9f5WVm9ZzfDew8s2mypTOCV3VN+epVkEmEu4U2Bmb0RDViKJSP273tzZdXFzGgvN7Eh3Xxzz9RJDJWZTZQp7HSXrceSi0usiVSFu0Pg48EUzWwu8STBzyt19TFINk12lz6aqlHBKbsk3dMqk0usiVSFu0Dgx0VZITSpLzarO5JtxpaErkcTFrXK7AfhH4IfAdcBnU+ekTCo5kwrKVLOqM/lmXGm2lUhJxO1p3AZs44MV4KcTbMJU/kH2BlTJmVShbCvHyy7fjCsNXYmURNygMdzdD0s7ftjMVLAwIeF6jZ6TJ2ctlx7Wpao2mcNVUKEhKxEpGW3CVGV6Tp5M9xEj2LFqVdbSItUqc7gKKjRklS7cw0NEEhN3ncaTwHAgfROmJ4H3qfAsqlpfpxHKt14DqmPNRmcqtqYDgmAx76sw6Nhg90ARyavU6zQmxnydJKQa1mx0pmJrOkB7eIiUSNzaU+vyPQq9j5lNNLPVZrbGzC7J8vy/mNntzY1xAAAPOElEQVRKM1tmZg+Z2aA47a1HU4dNZdbEWQzvPbziM6lymT5uIHeed/Quw1ZlFU7J1TCVSCLi5jSKZmZdgBsI1nyMAk43s1EZl/0FaE4Nd90NXFveVla/sMdRjYnxdOGeHGUVTsnV9FuRxFQsaABHAWvc/Vl3fxv4DTAl/QJ3f9jdt6cOFwL9y9zGigp388unGirgdiYcpip7Ury5Nchn5KueKyKRxM1p7MLMDnT3lyO8pB/wQtpxOzAuz/XnAL/L8d4zgBkAAwfWx/TOcDe/l6+4gq3z5uWcfgvVsW4jn3BPjmxTctNpeq5I9UuypxF1E6ZsO/1lncplZp8DmoHvZ3ve3W9y92Z3b25qaorYjOrUq2UaB151FXsdeWSn02/D3ka4n3g19jiyTclNV/LpuZp+K5KIxHoa7h51CW47MCDtuD/wUuZFZnYccCnwMXd/K34La0+vlmn0apnGujPP6hiqytfbgOqdTZW5gjxTSVeUjz41CBrL71YtKpEixQoaZnaNu/9bZ+c6sRgYamZDgBeB04DpGfc8HLgRmOjur8Rpaz0Ih6q2zpuXM2hUQwXcYpVs+ErTb0USE7encTyQGSBOzHIuJ3d/18wuAB4AugC3uPsKM5sJtLn7XILhqH2AOcGGgTzv7qfEbHPN6tUyja3z5nXa26hl6bsDZpPIWo9cFXFVDVekYJGChpmdD3wJONjMlqU91QN4LOqbu/t8YH7GucvTfj8u6j3rVSG9jVpW8uGrbJs4gTZyEokoak/j1wQzmL4LpC/G2+buWxJrlewi7G10VsywnhVVEDFXRVxVwxWJJFLQcPfXgdfN7B5gi7tvM7PLgCPM7N/d/S8laaUAQW8Dgiq4QM6gEc6iClVrbaoosg1fVaQ8iUiDi5vT+Ja7zzGzY4FPAT8Afk7+dRZSpPTZVLmEs6hCbRvaOtZwhM/XYgDJNnxVkT08RBpc3KDxXurnScDP3P0+M7symSZJIXIlxcNZVKE5T83pCBjVOh23GPlmXGmxoEjy4gaNF83sRoJZVNeYWTcqW5KkoURJiqcHkbCcer0MXeWbcRVp6CrfPuOF0gwsaRBxg8Y0gvLoP3D318ysL3Bxcs2SfMKkeFSZQ1e13vPIN+Oq5cYFHUUS8waOXLOqotAMLGkgcYPG34C9CfYGnwl0BV5LqlFSmKgzqTKHrlrvb+0odFirgSOXKWP7sei5Ldy39MX8QSPfPuOF0gwsaSBxh5R+CownCBoA2wjKnEuZJLEtbK2UVY9j+riBjBvSu9LNEKk7cYPGOHf/MrADwN1fBfZIrFXSqV4t0xh0+210HzGioBLq2WQWOgwf1VjwMK4wUV72vTxE6lTc4al3UpsoOYCZNRHsDy5lVuxK8c6m6IbX1OLwVVm3m00imV4IJdylwuIGjR8D9wL7m9nVwKnAZYm1SgpWbF2qfFN0obaT5WGivOTrOZJIphdCCXepArGChrv/ysyWAJ8k2Bfj0+7+ZKItk4JF2bCpM/WaLC9oJlVcSSTTC6GEu1SB2Gsr3H2Vu9/g7j9RwKisKBs2RVUPyfKKbTcrUoe0IK9OpCfGk5QtWV5riXLNpBJJTmI790n9Sk+WpyfKay1BnqvkSE2VGyk04a6EuZSIgkYdChf9hYoto56e5wgT5bWWIM9VcqSmKuUWmnBXwlxKyNy90m1IVHNzs7e1tVW6GRXz6p137ZTT2LFqFd1HjGDQ7bcl+j5hHavhvYd3nKu1ngd8UCn3zvOOrnBLEhT2RFp/W9l2SE0xsyXu3tzZdepp1JmwfHpo3Zln7dTzSGrzpnqrYyUihVHQqHPhxk0QlFMPFwKmPx8niGSbmptZQRdqo/eRr7x6oWoqLyJSBAWNOpfe88gcugqDSHhdMTJ7HlAbSfN85dULVZV5kXKtUJfKK/OkB+U0Gtird97Fy1dcAcBeRx7ZcT6pIaz0pPnw3sOZNXFW0fesRlWXF2mbBcvvrnQrpBxeXg4Hjk4kf6WchnQqDAyZifP054oRDmHVy6ryfJIY4ipEQcNg5VqhLpVXgd6kgkaD6yxxDsX3PCYdNIm2DW3MXDAz78ryah3C6kwSQ1yFqMphMGk4Chqyk/TEOSST9wgDQb6Aka26bqZqDSr5dhBMUjl6MiKdUU5D8sqW90gq55Eus7pupnrPixSi5cYFrFy/lVF9e3Z6rWZzNYgE1+QopyGJyMx7pE/bTTJ4ZE7hzVTLU3qTUugwmIaxGkz6TLkyzKRS0JBOZZu2m23NRzalWkwInQ9p1VtAKXQYTMNYDSS9tEyZysdoeEpiyVzzkU2YC9nryCPLPqTVtiH4b+Dyoy+vq8BRiCjDWFL7OoYiixyq0vCUlFTmrKtsMnsl4euSkm9Ia85Tc5i5YGbHjK1663XkU67ZXFJ5lRiKVNCQkgkDS5hMD3cWjCtKbyV9xlY4jJV+vp6VazaXVF4lhiIrOjxlZhOB64EuwC/c/XsZz3cDbgM+AmwGWtx9bb57aniqOhUynJVP+lBXukICSdjrAGg+YNfedyP1QqS+pA9FXr75Yvrs040DvvJQrHtV/fCUmXUBbgCOB9qBxWY2191Xpl12DvCqu/+dmZ0GXAO0lL+1UqxChrPyyRZ0Ck3GHwX88MhJ/GrYK7s8V8j6kEIo8EglpA9Fbn/7PTa98RYHlPg9K9bTMLOjgSvd/VOp428AuPt30655IHXNAjPbHXgZaPI8jVZPo3EU2nvJ1UsB2Pi3jWz+2+ai2rHtnW08Mmo3Xp94VFH3SZKCWOOZf+3ZAEz619mxXl/1PQ2gH/BC2nE7MC7XNe7+rpm9DuwHbCpLC6WqFdp7yRdcmvZsomnPpqLasX3xYg55/n2ef3pFUfdJyrZ3tgELeaDr9yvdFCmjt95+jxf371Py96lk0LAs5zJ7EIVcg5nNAGYADByoBKDsrNihsc6EQWlEyd4hmiR6T1J7dtvN6L57l5K/TyWDRjswIO24P/BSjmvaU8NTHwK2ZN7I3W8CboJgeKokrRXJodRBKapBlW6A1LXdKvjei4GhZjbEzPYATgPmZlwzFzg79fupwB/z5TNERKS0KtbTSOUoLgAeIJhye4u7rzCzmUCbu88FbgZuN7M1BD2M0yrVXhERqfDiPnefD8zPOHd52u87AE0BERGpEpUcnhIRkRqjoCEiIgVT0BARkYIpaIiISMEUNEREpGB1twmTmW0E1sV8eR8ar0SJPnNj0GduDMV85kHu3mlNnboLGsUws7ZCCnbVE33mxqDP3BjK8Zk1PCUiIgVT0BARkYIpaOzspko3oAL0mRuDPnNjKPlnVk5DREQKpp6GiIgUTEEjxcwmmtlqM1tjZpdUuj1JMbMBZvawmT1pZivM7MLU+d5m9nszezr1s1fqvJnZj1N/h2VmdkRlP0E8ZtbFzP5iZvNSx0PMbFHq896ZKsePmXVLHa9JPT+4ku0uhpnta2Z3m9mq1Pd9dD1/z2Z2Ueq/6SfM7A4z616P37OZ3WJmr5jZE2nnIn+vZnZ26vqnzezsbO9VCAUNgn9ggBuAE4FRwOlmNqqyrUrMu8DX3H0kMB74cuqzXQI85O5DgYdSxxD8DYamHjOAn5W/yYm4EHgy7fga4Iepz/sqcE7q/DnAq+7+d8APU9fVquuB+919BHAYweevy+/ZzPoBXwGa3f1Qgu0VTqM+v+dbgYkZ5yJ9r2bWG7iCYEvto4ArwkATmbs3/AM4Gngg7fgbwDcq3a4Sfdb7gOOB1UDf1Lm+wOrU7zcCp6dd33FdrTwIdoF8CPgEMI9g2+BNwO6Z3zfBfi5Hp37fPXWdVfozxPjMPYHnMtter98z0A94Aeid+t7mAZ+q1+8ZGAw8Efd7BU4Hbkw7v9N1UR7qaQTC/wBD7alzdSXVJT8cWAQc4O7rAVI/909dVg9/ix8B/wq8nzreD3jN3d9NHad/po7Pm3r+9dT1teYgYCMwKzUs9wsz25s6/Z7d/UXgB8DzwHqC720J9f89h6J+r4l93woaActyrq6mlZnZPsB/AV919635Ls1yrmb+FmY2GXjF3Zekn85yqRfwXC3ZHTgC+Jm7Hw68yQdDFtnU9OdODa1MAYYAHwb2JhiayVRv33Nncn3OxD6/gkagHRiQdtwfeKlCbUmcmXUlCBi/cvd7Uqc3mFnf1PN9gVdS52v9b3EMcIqZrQV+QzBE9SNgXzMLd6pM/0wdnzf1/IcIthauNe1Au7svSh3fTRBE6vV7Pg54zt03uvs7wD3A31P/33Mo6vea2PetoBFYDAxNzbzYgyChNrfCbUqEmRnBXutPuvt1aU/NBcIZFGcT5DrC82elZmGMB14Pu8G1wN2/4e793X0wwff4R3c/A3gYODV1WebnDf8Op6aur7n/B+ruLwMvmNnw1KlPAiup0++ZYFhqvJntlfpvPPy8df09p4n6vT4AnGBmvVK9tBNS56KrdIKnWh7AJOAp4Bng0kq3J8HPdSxBN3QZsDT1mEQwnvsQ8HTqZ+/U9UYwk+wZYDnB7JSKf46Yn30CMC/1+0HAn4E1wBygW+p899TxmtTzB1W63UV83rFAW+q7/m+gVz1/z8BVwCrgCeB2oFs9fs/AHQR5m3cIegznxPlegc+nPv8aoDVue7QiXERECqbhKRERKZiChoiIFExBQ0RECqagISIiBVPQEBGRgiloiIhIwRQ0RESkYAoaIhGl9q34UtrxYyV6n/5m1lKKe4vEpaAhEt2+QEfQcPe/L9H7fJKgfpRI1VDQEInue8DBZrbUzL5vZm9AUHo+tWveL1K7yf3KzI4zs0dTu6UdFd7AzD5nZn9O3ePG1EZgpD1/LHAdcGrqmiFl/YQiOaiMiEhEqX1J5nmwYxxm9oa775M6v4Zgz5IVBIUw/0pQK+gUgno/nzazkcC1wGfd/R0z+ymw0N1vy3if+4Gvu/sTiFSJ3Tu/REQieM7dlwOY2QqCLTndzJYT7L4GwbDTR4DFQYFW9uSD0tbphhPsvCZSNRQ0RJL1Vtrv76cdv88H/3szYLa7fyPXTcxsP4Ky1u+UpJUiMSmnIRLdNqBHEa9/iCBXsT+AmfU2s0EZ1wyhtjZFkgahoCESkbtvBh5NJbu/H+P1K4HLgAfNbBnwe6BvxmWrgD6p9yjV7CyRyJQIFxGRgqmnISIiBVPQEBGRgiloiIhIwRQ0RESkYAoaIiJSMAUNEREpmIKGiIgUTEFDREQK9v8BN3sNwSm6GB0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pred_surv = estimator.predict_survival_function(x_new)\n",
    "for i, c in enumerate(pred_surv):\n",
    "    plt.step(c.x, c.y, where=\"post\", label=\"Sample %d\" % (i + 1))\n",
    "plt.ylabel(\"est. probability of survival $\\hat{S}(t)$\")\n",
    "plt.xlabel(\"time $t$\")\n",
    "plt.legend(loc=\"best\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Measuring the Performance of Survival Models\n",
    "\n",
    "Once we fit a survival model, we usually want to assess how well a model can actually predict survival. Our test data is usually subject to censoring too, therefore metrics like root mean squared error or correlation are unsuitable. Instead, we use generalization of the area under the receiver operating characteristic (ROC) curve called [Harrell's concordance index](https://pdfs.semanticscholar.org/7705/392f1068c76669de750c6d0da8144da3304d.pdf) or c-index.\n",
    "\n",
    "The interpretation is identical to the traditional area under the [ROC curve](https://en.wikipedia.org/wiki/Receiver_operating_characteristic) metric for binary classification:\n",
    "- a value of 0.5 denotes a random model,\n",
    "- a value of 1.0 denotes a perfect model,\n",
    "- a value of 0.0 denotes a perfectly wrong model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7362562471603816"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sksurv.metrics import concordance_index_censored\n",
    "\n",
    "prediction = estimator.predict(data_x_numeric)\n",
    "result = concordance_index_censored(data_y[\"Status\"], data_y[\"Survival_in_days\"], prediction)\n",
    "result[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or alternatively"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7362562471603816"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "estimator.score(data_x_numeric, data_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our model's c-index indicates that the model clearly performs better than random, but is also far from perfect."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Feature Selection: Which Variable is Most Predictive?\n",
    "\n",
    "The model above considered all available variables for prediction. Next, we want to investigate which single variable is the best risk predictor. Therefore, we fit a Cox model to each variable individually and record the c-index on the training set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Karnofsky_score          0.709280\n",
       "Celltype=smallcell       0.572581\n",
       "Celltype=large           0.561620\n",
       "Celltype=squamous        0.550545\n",
       "Treatment=test           0.525386\n",
       "Age_in_years             0.515107\n",
       "Months_from_Diagnosis    0.509030\n",
       "Prior_therapy=yes        0.494434\n",
       "dtype: float64"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def fit_and_score_features(X, y):\n",
    "    n_features = X.shape[1]\n",
    "    scores = np.empty(n_features)\n",
    "    m = CoxPHSurvivalAnalysis()\n",
    "    for j in range(n_features):\n",
    "        Xj = X[:, j:j+1]\n",
    "        m.fit(Xj, y)\n",
    "        scores[j] = m.score(Xj, y)\n",
    "    return scores\n",
    "\n",
    "scores = fit_and_score_features(data_x_numeric.values, data_y)\n",
    "pd.Series(scores, index=data_x_numeric.columns).sort_values(ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Karnofsky_score` is the best variable, whereas `Months_from_Diagnosis` and `Prior_therapy='yes'` have almost no predictive power on their own.\n",
    "\n",
    "Next, we want to build a parsimonious model by excluding irrelevant features. We could use the ranking from above, but would need to determine what the optimal cut-off should be. Luckily, scikit-learn has built-in support for performing grid search.\n",
    "\n",
    "First, we create a pipeline that puts all the parts together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.feature_selection import SelectKBest\n",
    "from sklearn.pipeline import Pipeline\n",
    "\n",
    "pipe = Pipeline([('encode', OneHotEncoder()),\n",
    "                 ('select', SelectKBest(fit_and_score_features, k=3)),\n",
    "                 ('model', CoxPHSurvivalAnalysis())])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to define the range of parameters we want to explore during grid search. Here, we want to optimize the parameter `k` of the `SelectKBest` class and allow `k` to vary from 1 feature to all 8 features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean_fit_time</th>\n",
       "      <th>std_fit_time</th>\n",
       "      <th>mean_score_time</th>\n",
       "      <th>std_score_time</th>\n",
       "      <th>param_select__k</th>\n",
       "      <th>params</th>\n",
       "      <th>split0_test_score</th>\n",
       "      <th>split1_test_score</th>\n",
       "      <th>split2_test_score</th>\n",
       "      <th>mean_test_score</th>\n",
       "      <th>std_test_score</th>\n",
       "      <th>rank_test_score</th>\n",
       "      <th>split0_train_score</th>\n",
       "      <th>split1_train_score</th>\n",
       "      <th>split2_train_score</th>\n",
       "      <th>mean_train_score</th>\n",
       "      <th>std_train_score</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.126878</td>\n",
       "      <td>0.008825</td>\n",
       "      <td>0.005426</td>\n",
       "      <td>0.000059</td>\n",
       "      <td>3</td>\n",
       "      <td>{'select__k': 3}</td>\n",
       "      <td>0.650628</td>\n",
       "      <td>0.718220</td>\n",
       "      <td>0.754649</td>\n",
       "      <td>0.707491</td>\n",
       "      <td>0.043067</td>\n",
       "      <td>1</td>\n",
       "      <td>0.765569</td>\n",
       "      <td>0.695647</td>\n",
       "      <td>0.714082</td>\n",
       "      <td>0.725099</td>\n",
       "      <td>0.029590</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.125702</td>\n",
       "      <td>0.008785</td>\n",
       "      <td>0.005293</td>\n",
       "      <td>0.000101</td>\n",
       "      <td>5</td>\n",
       "      <td>{'select__k': 5}</td>\n",
       "      <td>0.644874</td>\n",
       "      <td>0.738347</td>\n",
       "      <td>0.728822</td>\n",
       "      <td>0.703834</td>\n",
       "      <td>0.042098</td>\n",
       "      <td>2</td>\n",
       "      <td>0.783946</td>\n",
       "      <td>0.698207</td>\n",
       "      <td>0.718498</td>\n",
       "      <td>0.733551</td>\n",
       "      <td>0.036585</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.125528</td>\n",
       "      <td>0.008009</td>\n",
       "      <td>0.005279</td>\n",
       "      <td>0.000040</td>\n",
       "      <td>4</td>\n",
       "      <td>{'select__k': 4}</td>\n",
       "      <td>0.650628</td>\n",
       "      <td>0.719809</td>\n",
       "      <td>0.725723</td>\n",
       "      <td>0.698523</td>\n",
       "      <td>0.034138</td>\n",
       "      <td>3</td>\n",
       "      <td>0.768121</td>\n",
       "      <td>0.691037</td>\n",
       "      <td>0.707327</td>\n",
       "      <td>0.722162</td>\n",
       "      <td>0.033172</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.129095</td>\n",
       "      <td>0.011951</td>\n",
       "      <td>0.005358</td>\n",
       "      <td>0.000054</td>\n",
       "      <td>2</td>\n",
       "      <td>{'select__k': 2}</td>\n",
       "      <td>0.630753</td>\n",
       "      <td>0.717161</td>\n",
       "      <td>0.747934</td>\n",
       "      <td>0.698256</td>\n",
       "      <td>0.049604</td>\n",
       "      <td>4</td>\n",
       "      <td>0.758550</td>\n",
       "      <td>0.683611</td>\n",
       "      <td>0.705638</td>\n",
       "      <td>0.715933</td>\n",
       "      <td>0.031448</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.132264</td>\n",
       "      <td>0.016141</td>\n",
       "      <td>0.005395</td>\n",
       "      <td>0.000110</td>\n",
       "      <td>1</td>\n",
       "      <td>{'select__k': 1}</td>\n",
       "      <td>0.630753</td>\n",
       "      <td>0.715042</td>\n",
       "      <td>0.737087</td>\n",
       "      <td>0.693982</td>\n",
       "      <td>0.045843</td>\n",
       "      <td>5</td>\n",
       "      <td>0.744640</td>\n",
       "      <td>0.676697</td>\n",
       "      <td>0.695246</td>\n",
       "      <td>0.705527</td>\n",
       "      <td>0.028675</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.125061</td>\n",
       "      <td>0.008772</td>\n",
       "      <td>0.005272</td>\n",
       "      <td>0.000049</td>\n",
       "      <td>6</td>\n",
       "      <td>{'select__k': 6}</td>\n",
       "      <td>0.657427</td>\n",
       "      <td>0.669492</td>\n",
       "      <td>0.724690</td>\n",
       "      <td>0.683572</td>\n",
       "      <td>0.029179</td>\n",
       "      <td>6</td>\n",
       "      <td>0.783946</td>\n",
       "      <td>0.698848</td>\n",
       "      <td>0.716160</td>\n",
       "      <td>0.732985</td>\n",
       "      <td>0.036722</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.125344</td>\n",
       "      <td>0.009078</td>\n",
       "      <td>0.005225</td>\n",
       "      <td>0.000024</td>\n",
       "      <td>7</td>\n",
       "      <td>{'select__k': 7}</td>\n",
       "      <td>0.654812</td>\n",
       "      <td>0.659958</td>\n",
       "      <td>0.714876</td>\n",
       "      <td>0.676269</td>\n",
       "      <td>0.027083</td>\n",
       "      <td>7</td>\n",
       "      <td>0.788412</td>\n",
       "      <td>0.695519</td>\n",
       "      <td>0.712133</td>\n",
       "      <td>0.732021</td>\n",
       "      <td>0.040447</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.124911</td>\n",
       "      <td>0.009041</td>\n",
       "      <td>0.005294</td>\n",
       "      <td>0.000054</td>\n",
       "      <td>8</td>\n",
       "      <td>{'select__k': 8}</td>\n",
       "      <td>0.656904</td>\n",
       "      <td>0.653602</td>\n",
       "      <td>0.716942</td>\n",
       "      <td>0.675516</td>\n",
       "      <td>0.029004</td>\n",
       "      <td>8</td>\n",
       "      <td>0.786371</td>\n",
       "      <td>0.695006</td>\n",
       "      <td>0.713692</td>\n",
       "      <td>0.731690</td>\n",
       "      <td>0.039411</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   mean_fit_time  std_fit_time  mean_score_time  std_score_time  \\\n",
       "2       0.126878      0.008825         0.005426        0.000059   \n",
       "4       0.125702      0.008785         0.005293        0.000101   \n",
       "3       0.125528      0.008009         0.005279        0.000040   \n",
       "1       0.129095      0.011951         0.005358        0.000054   \n",
       "0       0.132264      0.016141         0.005395        0.000110   \n",
       "5       0.125061      0.008772         0.005272        0.000049   \n",
       "6       0.125344      0.009078         0.005225        0.000024   \n",
       "7       0.124911      0.009041         0.005294        0.000054   \n",
       "\n",
       "  param_select__k            params  split0_test_score  split1_test_score  \\\n",
       "2               3  {'select__k': 3}           0.650628           0.718220   \n",
       "4               5  {'select__k': 5}           0.644874           0.738347   \n",
       "3               4  {'select__k': 4}           0.650628           0.719809   \n",
       "1               2  {'select__k': 2}           0.630753           0.717161   \n",
       "0               1  {'select__k': 1}           0.630753           0.715042   \n",
       "5               6  {'select__k': 6}           0.657427           0.669492   \n",
       "6               7  {'select__k': 7}           0.654812           0.659958   \n",
       "7               8  {'select__k': 8}           0.656904           0.653602   \n",
       "\n",
       "   split2_test_score  mean_test_score  std_test_score  rank_test_score  \\\n",
       "2           0.754649         0.707491        0.043067                1   \n",
       "4           0.728822         0.703834        0.042098                2   \n",
       "3           0.725723         0.698523        0.034138                3   \n",
       "1           0.747934         0.698256        0.049604                4   \n",
       "0           0.737087         0.693982        0.045843                5   \n",
       "5           0.724690         0.683572        0.029179                6   \n",
       "6           0.714876         0.676269        0.027083                7   \n",
       "7           0.716942         0.675516        0.029004                8   \n",
       "\n",
       "   split0_train_score  split1_train_score  split2_train_score  \\\n",
       "2            0.765569            0.695647            0.714082   \n",
       "4            0.783946            0.698207            0.718498   \n",
       "3            0.768121            0.691037            0.707327   \n",
       "1            0.758550            0.683611            0.705638   \n",
       "0            0.744640            0.676697            0.695246   \n",
       "5            0.783946            0.698848            0.716160   \n",
       "6            0.788412            0.695519            0.712133   \n",
       "7            0.786371            0.695006            0.713692   \n",
       "\n",
       "   mean_train_score  std_train_score  \n",
       "2          0.725099         0.029590  \n",
       "4          0.733551         0.036585  \n",
       "3          0.722162         0.033172  \n",
       "1          0.715933         0.031448  \n",
       "0          0.705527         0.028675  \n",
       "5          0.732985         0.036722  \n",
       "6          0.732021         0.040447  \n",
       "7          0.731690         0.039411  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.model_selection import GridSearchCV\n",
    "\n",
    "param_grid = {'select__k': np.arange(1, data_x_numeric.shape[1] + 1)}\n",
    "gcv = GridSearchCV(pipe, param_grid, return_train_score=True, cv=3, iid=True)\n",
    "gcv.fit(data_x, data_y)\n",
    "\n",
    "pd.DataFrame(gcv.cv_results_).sort_values(by='mean_test_score', ascending=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results show that it is sufficient to select the 3 most predictive features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Celltype=large       -0.067277\n",
       "Celltype=smallcell    0.271007\n",
       "Karnofsky_score      -0.031285\n",
       "dtype: float64"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pipe.set_params(**gcv.best_params_)\n",
    "pipe.fit(data_x, data_y)\n",
    "\n",
    "encoder, transformer, final_estimator = [s[1] for s in pipe.steps]\n",
    "pd.Series(final_estimator.coef_, index=encoder.encoded_columns_[transformer.get_support()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## What's next?\n",
    "\n",
    "Cox's proportional hazards model is by far the most popular survival model, because once trained, it is easy to interpret. However, if prediction performance is the main objective, more sophisticated, non-linear or ensemble models might lead to better results. Check-out\n",
    "[this notebook](https://nbviewer.jupyter.org/github/sebp/scikit-survival/blob/master/examples/evaluating-survival-models.ipynb) for getting a better understanding on how to evaluate survival models,\n",
    "and [this notebook](https://nbviewer.jupyter.org/github/sebp/scikit-survival/blob/master/examples/survival-svm.ipynb) to learn more about Kernel Survival Support Vector Machines. The [API reference](https://scikit-survival.readthedocs.io/en/latest/api.html) contains a full list of models that are available within **scikit-survival**. In addition, you can use any unsupervised pre-processing method available with scikit-learn, for instance, you could perform dimensionality reduction using [Non-Negative Matrix Factorization (NMF)](http://scikit-learn.org/dev/modules/generated/sklearn.decomposition.NMF.html#sklearn.decomposition.NMF), before training a Cox model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}