doc/user_guide/00-introduction.ipynb

{
 "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](https://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'-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's-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 described 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 are 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 observed 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": {
    "tags": []
   },
   "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 few 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": {
    "tags": []
   },
   "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": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'time $t$')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAERCAYAAABl3+CQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlsUlEQVR4nO3de1SUdf4H8PcImsniBVQExhuO4oCgIqSImcoKRi4qti6meeUQalnrVtp2ftaWm7qdXXPVZMmWNFPyuJvsmmCG1pFQLspqB5HlIBqDFhdRSQMRvr8/bGa5DjwPzzC39+sczuGZeXjm8+105uP39vmqhBACREREbehm7gCIiMiyMVEQEZFRTBRERGQUEwURERnFREFEREYxURARkVFmTRQrVqzAwIEDMWbMmFbfF0Jg7dq10Gg08Pf3x/nz57s4QiIiMmuiWLZsGVJTU9t8PyUlBYWFhSgsLERCQgJWrVrVhdERERFg5kQxdepUuLi4tPl+cnIylixZApVKhUmTJuHWrVu4ceNGF0ZIREQWPUdRWlqKwYMHG67VajVKS0vNGBERkf1xNHcAxrRWXUSlUrV4LSEhAQkJCQCAy5cvY/To0SaPjYjIlly9ehUVFRWtvmfRiUKtVqOkpMRwrdPp4OHh0eK+2NhYxMbGAgACAwORk5PTZTESEdmCwMDANt+z6KGnyMhI7Nu3D0IInD17Fn369IG7u7vJPu8P/85D3P5zJns+EZE1MmuPYuHChfjqq69QUVEBtVqNP/zhD6irqwMAxMXFISIiAseOHYNGo0GvXr2QmJho0nguXb+DOzV1Jv0MIiJrY9ZEcfDgQaPvq1Qq7Nq1q4uiISKi1lj00BMREZkfE0Uz+TeqcSDzO3OHQURkMZgoGpkzzhMA8PHZq+YNhIjIgjBRNPLMxCHQujubOwwiIovCREFEREYxURARkVFMFEREZBQTRSu48omI6H+YKJoJGdEfAJD8H1apJSICmChaCNW6QevuzFIeREQ/Y6IgIiKjmCiIiMgoJopmgke4mjsEIiKLwkRBRERGMVEQEZFRTBRGnCmqNHcIRERmx0RBRERGMVG0If9GNdLyfzB3GEREZsdE0Qr97uxviirMHAkRkfkxUbRCvzsb4DwFERETRSu4l4KI6H+YKIiIyCgmCiIiMoqJgoiIjGKiICIio5gojOBeCiIiJoo2cS8FEdFDTBRtaLyXgojInjFREBGRUUwUHcDd2URkz5go2sEJbSKyd51KFHfv3kV9fb1SsVgcTmgTEQGOUm5uaGhAUlISPvnkE2RnZ+ORRx5BbW0tBgwYgIiICMTGxmLkyJGmirVL6es9MUkQkb2T1KOYPn06ioqKsHnzZnz//fcoKSlBWVkZTp8+jUmTJmHDhg3Yv3+/qWI1K85TEJG9ktSj+PLLL9G9e3dcu3YN3br9L8e4uLhg/vz5mD9/Purq6hQPkoiIzEdSj6J79+4AgHnz5rV47+zZs03uISIi2yApURw6dAgbNmxAdXU18vPzm0xkx8bGSv7w1NRUeHt7Q6PRYMuWLS3ev337Nn71q19h7Nix8PX1RWJiouTPICKizpE09BQSEoKamhrs2bMH69atQ0FBAfr27QsPDw88+uijkj64vr4ea9aswYkTJ6BWqxEUFITIyEj4+PgY7tm1axd8fHzw73//G+Xl5fD29saiRYvQo0cPSZ9FRETySUoUnp6eWLJkCUaMGIGQkBAAwM2bN1FcXIzRo0dL+uCsrCxoNBp4eXkBAKKjo5GcnNwkUahUKlRXV0MIgR9//BEuLi5wdJQUMhERdZKkb10hBFQqlSFJAA8nsl1cXFrc057S0lIMHjzYcK1Wq5GZmdnknueffx6RkZHw8PBAdXU1Pv300yaT6HoJCQlISEgAAJSXl0tpEhERtUPy8tgdO3bgu+++a/L6/fv3cfLkSSxduhR79+7t0LOEEC1ea55gjh8/jnHjxuH69ev4z3/+g+effx537txp8XexsbHIyclBTk4OBgwYIKFF0nCJLBHZI0mJIjU1FQ4ODli4cCE8PDzg4+OD4cOHY+TIkTh48CB++9vfYtmyZR16llqtRklJieFap9PBw8OjyT2JiYmIioqCSqWCRqPB8OHDcfnyZSkhd4p+0x3LeBCRPZM09NSzZ0+sXr0aq1evRl1dHSoqKvDoo4+ib9++kj84KCgIhYWFKC4uhqenJ5KSknDgwIEm9wwZMgRpaWl4/PHH8cMPP6CgoMAwp9FVQkb0R/6NanxTVIFQrVuXfjYRkSUw2qO4dOkSFi9ebLgODQ1FXl4egIf7JbKzs7Fz505kZWVJ/mBHR0fs3LkT4eHh0Gq1WLBgAXx9fREfH4/4+HgAwP/93/8hIyMDfn5+CA0NxdatW9G/f3/Jn9UZzc+l4PATEdkboz2K0NBQnDlzxnCt0+ng6+sLAMjIyMCzzz6L3/zmN1i2bBn++Mc/troRz5iIiAhEREQ0eS0uLs7wu4eHB7744gtJzyQiImUZ7VF88cUXeP311w3XvXv3Nvy+b98+xMXFISEhAV999RW2bt1quiiJiMhsjCYKPz8/fPLJJ4ZrjUaDw4cPo6ysDEeOHMGcOXMAAAMHDkRtba1pIzUzTmgTkb2StOpp27Zt+Nvf/gZPT08EBARg8uTJAIC6ujr8+OOPJgnQ3IJHuPJcCiKya5JWPQ0aNAgnTpxAQ0NDk41vp06dwvTp0xUPzlKEat2YJIjIbsmqh9F8d3RYWBjCwsIUCYiIiCwLz8wmIiKjmChk4F4KIrInTBRERGSUpDkKZ2fnVivD6ivGtlawj4iIrJukRFFdXW2qOKyCfi9FqNbNMPykLxxIRGSrZJ8CVFVVhcLCQtTU1Bhemzp1qiJBWSIWByQieyUrUezZswfbt2+HTqfDuHHjcPbsWQQHB+PkyZNKx2cxuJeCiOyVrMns7du3Izs7G0OHDsWpU6eQm5tr0gODLEnzUh5cAUVEtk5WoujZsyd69uwJAKitrcXo0aNRUFCgaGCWiKU8iMgeyRp6UqvVuHXrFubOnYuZM2eiX79+LU6nsyXBI1xxpqiSw09EZJdkJYrPPvsMAPDmm29i+vTpuH37NmbNmqVoYNbkTFElVz8Rkc2SNfS0bds26HQ6AMATTzyByMhI9OjRQ9HAiIjIMshKFHfu3EF4eDgef/xx7Nq1Cz/8wHMaiIhslaxE8cYbbyAvLw+7du3C9evX8cQTT+CXv/yl0rFZFa5+IiJb1alaTwMHDsSgQYPg6uqKsrIypWIiIiILIitR7N69G9OmTUNoaCgqKirwwQcf4OLFi0rHZlEaT1bzWFQisieyVj1du3YN7733HsaNG6dwOJaPpTyIyN7IShRbtmxROg6rwb0URGRvJA09TZkyBcDDcuO9e/c2/OiviYjI9kjqUaSnpwNguXEAuFZ5D28dzUPIiP4cgiIimyZ7w11paanSsViNkBH9MdS1F65V3uMwFBHZPNkb7sLCwux2w12o1g0bZ/tiqGsvc4dCRGRy3HDXSY2XynLTHRHZIm64k6B54T+WHScie8ANd50QqnWD1t3Z3GEQEZmU5H0UQgjk5OTY7Ya71uiHn7j6iYhskeQehUqlQm5uLpPEzzj8RES2TtbQU3BwMLKzs5WOxSo0n6fg8BMR2TpZJTxOnTqF+Ph4DBs2DE5OThBCQKVS2eU8hV7j1U887Y6IbImsRJGSkqJ0HFaNhQKJyJbJShR79+5t9fWNGzdKek5qaipefPFF1NfXIyYmBhs2bGhxz1dffYWXXnoJdXV16N+/P77++ms5IZtU80KBPEObiGyJrETh5ORk+L2mpgZHjx6FVquV9Iz6+nqsWbMGJ06cgFqtRlBQECIjI+Hj42O459atW1i9ejVSU1MxZMgQq9qrwWRBRLZCVqL43e9+1+T65ZdfRmRkpKRnZGVlQaPRwMvLCwAQHR2N5OTkJoniwIEDiIqKwpAhQwA83OBHRERdq1M7s/Xu3buHK1euSPqb0tJSDB482HCtVqtbFBr873//i6qqKkybNg0TJkzAvn37lAi3y5wpqmRZDyKyerJ6FH5+flCpVAAeDiGVl5dLnp8QQrR4Tf9MvQcPHuDcuXNIS0vDTz/9hODgYEyaNAmjRo1qcl9CQgISEhIAAOXl5ZLiICIi42QliqNHj/7vAY6OcHNzg6OjtEep1WqUlJQYrnU6HTw8PFrc079/fzg5OcHJyQlTp07FhQsXWiSK2NhYxMbGAgACAwOlNkey4BGuknoKnK8gImsma+gpKysLLi4uGDp0KBITE7FgwQKcP39e0jOCgoJQWFiI4uJi3L9/H0lJSS3mOebMmYPTp0/jwYMHuHfvHjIzMyVPmlsKDkMRkbWSlSjefvttODs7Iz09HcePH8fSpUuxatUqSc9wdHTEzp07ER4eDq1WiwULFsDX1xfx8fGIj48HAGi1WsyaNQv+/v547LHHEBMTgzFjxsgJmYiIZFKJ1iYL2jF+/Hjk5ubitddeg5+fH5555hnDa+YWGBiInJwck39O897BW0fzkH+jGjFThhvddMchKCKyRMa+O2X1KDw9PfHcc8/h0KFDiIiIQG1tLRoaGjoVpLUJHuHa5EufxQGJyFbJShSHDh1CeHg4UlNT0bdvX9y8eRPvvvuu0rFZFRYHJCJbJWvVU69evRAVFWW4dnd3h7u7u2JB2TKugCIia6PIhjt7xi99IrJ1TBRmwGWyRGRNJCWKZ599FgCwfft2kwRjT5gsiMhaSEoU586dw7Vr1/D3v/8dVVVVuHnzZpMfe9V4+KnxAUZERLZA0mR2XFwcZs2ahStXrmDChAlN6jWpVCrJhQFtjf4Aoz3pxQDAQ4yIyCZI6lGsXbsW+fn5WLFiBa5cuYLi4mLDj70nCeBhYoiZMhwA91MQke2QtTx29+7duHDhAk6fPg0AmDp1Kvz9/RUNzFrpT7vTD0EZ61VwqSwRWQNZq57++te/YtGiRSgrK0NZWRkWLVqEHTt2KB2bVeEubSKyVbJqPfn7++PMmTOGI1Hv3r2L4OBgXLx4UfEApeqqWk+tabyS6a2jebhWeQ9DXXshZET/NnsW7FEQkSUw9t0pa+hJCAEHBwfDtYODQ6sHEdmzh72Kh0NQ+TeqAbQ+ua1PLkwYRGSpZCWK5cuXY+LEiZg3bx4A4MiRI1i5cqWigVm7UK0bQrVuSMv/AXvSi/FNUQVXQRGRVZKVKNatW4dp06YhPT0dQggkJiZi/PjxSsdmE/ST29cq7+Gto3lGh6GIiCyRrEQBAAEBAQgICFAyFpulH4a6VnkPAHsWRGRdWOupC4Rq3bBxti+GuvYydyhERJIxUVgI1n4iIkslK1Hs3LkTVVVVSsdCREQWSFai+P777xEUFIQFCxYgNTWVS2MlYNFAIrI2shLFpk2bUFhYiJUrV+Kjjz7CyJEj8fvf/x5FRUVKx2dV2tsLwR3bRGSNZM9RqFQqDBo0CIMGDYKjoyOqqqrw9NNP49VXX1UyPpuiP1e7rV7FmaJKzlUQkcWRXetpwoQJePXVVxESEoJvv/0Wu3fvxrlz5/CPf/xD6Rhtir5XsSe9mENQRGQVZO2jqKiowD//+U8MHTq0yevdunXD0aNHFQnMVun3UHC3NhFZC1k9itra2hZJYv369QAArVbb+ahsnH4IiojIGshKFCdOnGjxWkpKSqeDISIiyyMpUezevRt+fn4oKCiAv7+/4Wf48OE8uOhnUqrAGpvUJiKyFJLmKJ555hk8+eSTeO2117BlyxbD687OznBxcVE8OFumP1+b8xREZOkkJYo+ffqgT58+OHjwoKnisQnBI1zb7RXoq8q2hcekEpGlkDT0NGXKFAAPexC9e/c2/Oiv6X86+iWvLz/OpbJEZKkk9SjS09MBANXV1SYJxt6w/DgRWQNWjzWh4BGuRnsWLD9ORNZAUo/C2dkZKpWq1SKAKpUKd+7cUSwwIiKyDJISBYec5OnI5HZr9H/DSW0iMidFJrP1P9S29oahWH6ciCyVpETReDL7zp07LX6kSk1Nhbe3NzQaTZN9Gc1lZ2fDwcEBhw8flvwZ1qC98uPcgEdE5mS2yez6+nqsWbMGKSkpuHTpEg4ePIhLly61et/69esRHh5uhii7Bms/EZElk5Uoampq8Je//AVRUVGYP38+tm3bhpqaGknPyMrKgkajgZeXF3r06IHo6GgkJye3uG/Hjh2YP38+Bg4cKCdUm8GzKojIXGQliiVLliAvLw8vvPACnn/+eeTn5+PZZ5+V9IzS0lIMHjzYcK1Wq1FaWtrins8++wxxcXFywrQ63HxHRJZI1nkUBQUFuHDhguF6+vTpGDt2rKRntLXEtrGXXnoJW7duhYODg9FnJSQkICEhAQBQXl4uKY6u1tYKKG6+IyJLJStRjB8/HmfPnsWkSZMAAJmZmQgJCZH0DLVajZKSEsO1TqeDh4dHk3tycnIQHR0N4OFhSceOHYOjoyPmzp3b5L7Y2FjExsYCAAIDA6U2xyKEat0QqnXDW0fzzB0KEVETkhKFn58fVCoV6urqsG/fPgwZMgQA8N1338HHx0fSBwcFBaGwsBDFxcXw9PREUlISDhw40OSe4uJiw+/Lli3D7NmzWyQJayR3XwXAvRVE1PUkJQoljzl1dHTEzp07ER4ejvr6eqxYsQK+vr6Ij48HAJufl+hMsiAi6kqSEkXj40+rqqpQWFjYZLVT8+NR2xMREYGIiIgmr7WVID766CNJz7Zm+s13nKcgIksga9XTnj17MHXqVISHh+ONN95AeHg43nzzTYVDs0/tbb7TY2+EiLqKrESxfft2ZGdnY+jQoTh16hRyc3MxYMAApWOzea3NM+g337GkBxFZClmJomfPnujZsycAoLa2FqNHj0ZBQYGigdmzjvYqiIi6gqzlsWq1Grdu3cLcuXMxc+ZM9OvXr8XSVpKvvWNSiYi6kqxE8dlnnwEA3nzzTUyfPh23b9/GrFmzFA2M2selskTUFWQlipqaGrz//vtIT0+HSqXClClT0NDQoHRsdkH/Jc/JaSKyVGar9UTt62jtJyYZIjIls9V6IuNY+4mILIWsHoW+1pOenFpPZFyo1g0bZ/tiqGuvDi2VZa+CiEzFbLWeqGNCRvRH/o1qfFPEXgURmYfZaj1RU23VftIvldXPV4SM6N9mwuAqKCIyBdm1ni5cuIDTp08DAB5//HHOUZgQ5yuIyJxkl/BYtGgRysrKUFZWhsWLF2PHjh1Kx0Y/azxfQUTU1WStevrwww+RmZkJJycnAMD69esRHByMF154QdHg7E1H9lR0ZAiKiEhJshKFEKLJ8aQODg6tHm1KyuIQFBGZg6xEsXz5ckycOBHz5s0DABw5cgQrV65UNDBqicelEpE5yEoU69atw7Rp05Ceng4hBBITEzF+/HilYyOZmg9dcRUUEXWG5EQhhIBOp0NAQAACAgJMERN1AE/BI6KuInnVk0qlwty5c00QCum11wOQel7FmaJKww8RkVSylsdOmjQJ2dnZSsdCHaQ/BY+IqCvImqM4deoU4uPjMWzYMDg5OUEIAZVKhYsXLyodn91qa6d2Z3H+goikkpUoUlJSlI6DZOCeCiLqCrIShZubW4uDi1atWqV0bGSEUnsqzhRVsldBREbx4CIrJbUMORGRXDy4yMrpy5DvSS8GAFk9C85bEJExPLjIggWPcG33SztU64aYKcMBdHy5LBGRFCoho0iTVqtFQUFBk4OLtFotunXrZvbVT4GBgcjJyTHb55tKeyug3jqah2uV9zDUtZfik9vsYRDZPmPfnbKGnlJTUzsVECmPBQOJyFRkJYrGBxhR12ivBLkpCwZyZRSRfZM1R0FERPZDVo+C7E/jngx7F0T2RbEexffff6/Uo4iIyIIo1qNYuXIlPv/8c6UeR22QclyqntKroNi7ILIviiUKJgnLoF/9pJd/oxr5N6oNeyxYF4qIpJI19LR+/foOvUZdT1/aQ/8TM2W4oST5tcp7im/K41kXRLZP1oa7gIAAnD9/vslr/v7+FlFm3FY33LVG6pdz4015ekr2MDgMRWS9jH13SupR7N69G35+frh8+TL8/f0NP8OHD4efn5/kwFJTU+Ht7Q2NRoMtW7a0eP+TTz4xfMbkyZOb1Jci6UJG9G+SJEzRwyAi2yOpR3H79m1UVVXhtddea/LF7uzsDBcXF0kfXF9fj1GjRuHEiRNQq9UICgrCwYMH4ePjY7gnIyMDWq0W/fr1Q0pKCt58801kZmYafS57FB331tE85N+oRsyU4Sabt2Avg8g6KNaj6NOnD4YNG4aoqCi4uLhg6NCh+PjjjxETE4Pc3FxJQWVlZUGj0cDLyws9evRAdHQ0kpOTm9wzefJk9OvXD8DD41d1Op2kz7B1HSkaaIzUs7eJyD7Jmsx+++234ezsjPT0dBw/fhxLly5FXFycpGeUlpZi8ODBhmu1Wo3S0tI27//www/x5JNPtvpeQkICAgMDERgYiPLycklx2DP92dv65bT6HyXPtuBEN5H1k5UoHBwcADxcErtq1SrMmTMH9+/fl/SM1ka8VCpVq/eeOnUKH374IbZu3drq+7GxscjJyUFOTg4GDBggKQ5b0NleBectiMgYWfsoPD098dxzz+HEiRNYv349amtr0dDQIOkZarUaJSUlhmudTgcPD48W9128eBExMTFISUmBqyvHu9sSPMJV1r/c9cUE9fQro5TesMfCgkTWS1aP4tChQwgPD8fx48fRt29f3Lx5E++++66kZwQFBaGwsBDFxcW4f/8+kpKSEBkZ2eSe7777DlFRUfj4448xatQoOaGSROxhEFFzsnoUjz76KO7evYuDBw9i48aNqKurQ9++faV9sKMjdu7cifDwcNTX12PFihXw9fVFfHw8ACAuLg5vvfUWKisrsXr1asPf2MuKJjk6Ut6jPa31MPRncivRqzCGPQ4iyyRrw92qVavQrVs3nDx5Evn5+aiqqkJYWBiys7NNEaMk9rQ8ti1KTh6n5f+APenF0Lo7Y+NsX8We2xomCiLzUfyEu8zMTJw/fx7jx48HAPTr10/yZDZZh1CtG74pqmgyb8F6UUT2RVai6N69O+rr6w2rlMrLy9GtG89AslWNCw2a8qhVDk0RWSZZiWLt2rWYN28eysrK8Prrr+Pw4cPYtGmT0rGRhWg8b9F8VRR7F0S2T1aiWLRoESZMmIC0tDQIIXDkyBFotVqlYyOZmv/LW8k5i67qXbSmtXawl0FkerImsy0ZJ7ONUzJpNK9Ga47eBRMFkTIUn8wmAszbu9AzlviYRIiUwURhZ5QcljI2d6HHOQwi68dEQYpofgQrYL5eBhEpi3MUZLLqro3nMMzVs+DwE1HHcI6CjDLVKil9L4M9CyLrxl1yZDKhWjdsnO2Loa69DPWiuhrPwyDqPPYoqAWlexghI/oj/0Y1vilir4LIGjFRkMm1Vi+qNVwhRWSZOPRE7VJiQrj5ORfNmfrcCw4/EcnHHgV1SONkocRJes21tQ+jMfY4iMyDiYIsQmv7MBpTYuVUWwmOS2iJjGOiIMk627tojdweB3sZRKbHREFWgTu/icyHO7NJMV09Ydy8em1zXdnb4PAVWTvuzKYu0dqXpSmTh7F5DfY2iJTDREFWy9i8xltH8wy7wZksiDqHiYJMSt/L6Ophqa7eDd5V7eMQF5kDN9xRl+jqL7hQrRu07s6GlVLmqDNFZCvYoyCbZYvVazvac2HPg5TEREFdpr0vL6WHb/RzGJyvIOocJgqyGKZKJKxeS9Q5TBRk89qrXmuLu7tZBNG+mHqokYmCrEZn9mm0tefCluYviEyFiYKsmrF/STVOIm3tuehI1dqOsMVeCZEeEwXZtfaq1nYEeyVkbvp/FJlqCIqJgmxW8AjXdoem2qta2xFK9Uo6gj0XMgcmCrJpXbEzXIleSUew50LmwkRB1ElK9Eo6QkrPhT0PUhJLeJBdsIWdyu2dO65n6vPHyf6wR0F2o61kYS17Djrac+nKOROyDKbuQZo1UaSmpuLFF19EfX09YmJisGHDhibvCyHw4osv4tixY+jVqxc++ugjBAQEmClaslVK9DYsKdl01ZwJWYaumLsyW6Kor6/HmjVrcOLECajVagQFBSEyMhI+Pj6Ge1JSUlBYWIjCwkJkZmZi1apVyMzMNFfIRG3qqqGtjiSkrpozIcvQuAf57KRheGbiEMU/w2yJIisrCxqNBl5eXgCA6OhoJCcnN0kUycnJWLJkCVQqFSZNmoRbt27hxo0bcHd3N1fYRGbV0YRkST0cMq3GVZKT/1NqW4mitLQUgwcPNlyr1eoWvYXW7iktLWWiIGqHLUzeU8fo9wuZck7KbIlCCNHiNZVKJfkeAEhISEBCQgIA4PLlywgMDJQVU3l5OQYMGCDrb60V22wf2Gb7cKm8HIEfyGvz1atX23zPbIlCrVajpKTEcK3T6eDh4SH5HgCIjY1FbGxsp2MKDAxETk5Op59jTdhm+8A22wdTtdls+yiCgoJQWFiI4uJi3L9/H0lJSYiMjGxyT2RkJPbt2wchBM6ePYs+ffpw2ImIqIuZrUfh6OiInTt3Ijw8HPX19VixYgV8fX0RHx8PAIiLi0NERASOHTsGjUaDXr16ITEx0VzhEhHZLbPuo4iIiEBEREST1+Li4gy/q1Qq7Nq1q8viUWL4ytqwzfaBbbYPpmqzSrQ2Y0xERPQz1noiIiKjmCh+lpqaCm9vb2g0GmzZssXc4SimpKQE06dPh1arha+vL7Zv3w4AuHnzJmbOnImRI0di5syZqKqqMvzN5s2bodFo4O3tjePHj5sr9E6pr6/H+PHjMXv2bAC2314AuHXrFp5++mmMHj0aWq0WZ86csel2b9u2Db6+vhgzZgwWLlyImpoam2zvihUrMHDgQIwZM8bwmpx2njt3Dn5+ftBoNFi7dm2r2w/aJEg8ePBAeHl5iaKiIlFbWyv8/f1FXl6eucNSxPXr18W5c+eEEELcuXNHjBw5UuTl5YlXXnlFbN68WQghxObNm8Wrr74qhBAiLy9P+Pv7i5qaGnHlyhXh5eUlHjx4YLb45frzn/8sFi5cKJ566ikhhLD59gohxJIlS8QHH3wghBCitrZWVFVV2Wy7dTqdGDZsmLh3754QQohf//rXIjEx0Sbb+/XXX4tz584JX19fw2ty2hkUFCQyMjJEQ0ODmDVrljh27FiHY2CiEEJkZGSIsLAww/U777wj3nnnHTNGZDqRkZHiiy++EKNGjRLXr18XQjxMJqNGjRJCtGx7WFiYyMjIMEuscpWUlIgZM2aItLQ0Q6Kw5fYKIcTt27fFsGHDRENDQ5PXbbXdOp1OqNVqUVlZKerq6sRTTz0ljh8/brPtLS4ubpIopLbz+vXrwtvb2/D6gQMHRGxsbIc/n0NPaLtUiK25evUqcnNzMXHiRPzwww+GPSnu7u4oKysDYBv/LV566SX86U9/Qrdu//vf25bbCwBXrlzBgAEDsHz5cowfPx4xMTG4e/euzbbb09MTL7/8MoYMGQJ3d3f06dMHYWFhNtve5qS2s7S0FGq1usXrHcVEgY6XCrFmP/74I+bPn4/33nsPvXv3bvM+a/9vcfToUQwcOBATJkzo0P3W3l69Bw8e4Pz581i1ahVyc3Ph5ORkdK7N2ttdVVWF5ORkFBcX4/r167h79y7279/f5v3W3t6OaqudnW0/EwU6XirEWtXV1WH+/PlYtGgRoqKiAABubm64ceMGAODGjRsYOHAgAOv/b/HNN9/gX//6F4YNG4bo6GicPHkSixcvttn26qnVaqjVakycOBEA8PTTT+P8+fM22+4vv/wSw4cPx4ABA9C9e3dERUUhIyPDZtvbnNR2qtVq6HS6Fq93FBMFOlZOxFoJIbBy5UpotVqsW7fO8HpkZCT27t0LANi7dy/mzJljeD0pKQm1tbUoLi5GYWEhHnvsMbPELsfmzZuh0+lw9epVJCUlYcaMGdi/f7/Ntldv0KBBGDx4MAoKCgAAaWlp8PHxsdl2DxkyBGfPnsW9e/cghEBaWhq0Wq3Ntrc5qe10d3eHs7Mzzp49CyEE9u3bZ/ibDpE3tWJ7Pv/8czFy5Ejh5eUlNm3aZO5wFHP69GkBQPj5+YmxY8eKsWPHis8//1xUVFSIGTNmCI1GI2bMmCEqKysNf7Np0ybh5eUlRo0aJWllhKU5deqUYTLbHtqbm5srJkyYIPz8/MScOXPEzZs3bbrdGzduFN7e3sLX11csXrxY1NTU2GR7o6OjxaBBg4Sjo6Pw9PQUe/bskdXO7Oxs4evrK7y8vMSaNWtaLHwwhjuziYjIKA49ERGRUUwURERkFBMFEREZxURBRERGMVEQEZFRTBRERGQUEwURERnFREEk0a1bt/D+++8bridPnmySz9HpdPj0009N8mwiKZgoiCRqnigyMjJM8jlpaWk4f/68SZ5NJAUTBZFEGzZsQFFREcaNG4dXXnkFv/jFLwA8LOM+evRoxMTEYMyYMVi0aBG+/PJLhISEYOTIkcjKyjI8Y//+/Xjssccwbtw4PPfcc6ivr2/yGenp6Vi3bh0OHz6McePGobi4uEvbSNSEgiVJiOxC80NknJycDK87ODiIixcvivr6ehEQECCWL18uGhoaxJEjR8ScOXOEEEJcunRJzJ49W9y/f18IIcSqVavE3r17W3xOeHi4+Pbbb03fIKJ2OJo7URHZkuHDh8PPzw8A4Ovri9DQUKhUKvj5+eHq1asAHg4pnTt3DkFBQQCAn376yVAmurGCggJ4e3t3WexEbWGiIFLQI488Yvi9W7duhutu3brhwYMHAB6Wfl+6dCk2b97c5nMqKyvRp08fdO/e3bQBE3UA5yiIJHJ2dkZ1dbXsvw8NDcXhw4cNx1fevHkT165da3JPcXGxVR+sQ7aFiYJIIldXV4SEhGDMmDF45ZVXJP+9j48PNm3ahLCwMPj7+2PmzJmG08r0Ro8ejYqKCowZM8Zkq6qIOornURARkVHsURARkVFMFEREZBQTBRERGcVEQURERjFREBGRUUwURERkFBMFEREZxURBRERG/T8EDYr0dI50RAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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, conf_int = kaplan_meier_estimator(\n",
    "    data_y[\"Status\"], data_y[\"Survival_in_days\"], conf_type=\"log-log\"\n",
    ")\n",
    "plt.step(time, survival_prob, where=\"post\")\n",
    "plt.fill_between(time, conf_int[0], conf_int[1], alpha=0.25, step=\"post\")\n",
    "plt.ylim(0, 1)\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 occurring 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": {
    "tags": []
   },
   "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": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f3d7de38490>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAERCAYAAABl3+CQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAz+klEQVR4nO3deVRUV7o28KcEDEojihNCKYOoDBaigorGKCHOinFKTDTOl4hGb+JNWnN7pTN5o4nd1zbRDk1iG+1Evcavo9EEjEG0pREVNZoGNLQCsXAGRByY9/cHqQpDUXAOp6jp+a3Faqo4dc570i5ezn73frdKCCFARETUiDbmDoCIiCwbEwURERnFREFEREYxURARkVFMFEREZBQTBRERGWXWRLFo0SJ069YN/fv3N/hzIQRWrlwJf39/hISE4OzZs60cIRERmTVRLFiwAImJiY3+PCEhAdnZ2cjOzkZ8fDxiY2NbMToiIgLMnCieeOIJuLu7N/rz/fv3Y968eVCpVBg2bBju3r2L69evt2KERERk0TWK/Px89OzZU/9arVYjPz/fjBEREdkfR3MHYIyh7iIqlarBe/Hx8YiPjwcAXLx4EQEBAbKv+aCsEi6PWfR/FiIixeXm5uLOnTsGf2bRvxHVajWuXr2qf63VauHp6dnguJiYGMTExAAAwsLCkJ6eLvuaJy4XIKJ3Z9mfJyKyRmFhYY3+zKKHnqKjo7Fjxw4IIZCWlgY3Nzf06NHDZNd7+0AG3jmYgZ0nfzbZNYiIrI1Znyiee+45HD16FHfu3IFarcbbb7+NiooKAMDSpUsxceJEfPvtt/D390f79u2xbds2k8eUdb0Ef0vLxfNDe5n8WkRE1sCsiWLXrl1Gf65SqbBly5ZWigZ4c0ow0i/morLVrkhEZPksukZBRA1VVFRAq9WitLTU3KGQFXJ2doZarYaTk1OzP8NEUVvCGiwpy0Fc2+XmjoSoUVqtFq6urvDx8TE4C5CoMUIIFBQUQKvVwtfXt9mfs+hidqu78SP8qvPgWF5i7kiIGlVaWorOnTszSZBkKpUKnTt3lvw0ykRBZIWYJEguOf92mCjq0VRnof/dJE6RJWpEQUEBQkNDERoaCg8PD3h5eelfl5eXyzrn0aNHkZqaqnCkrX+dzz77DNeuXVPsfKNHj27RurCjR49i8uTJLY6DiaI2zUwAwFSHVOz/ga1CiAzp3LkzfvjhB/zwww9YunQpXnnlFf3rtm3borJS+rxBJgplVFVVmeS8LGbXFrYQD/75F7gUAw/vFZo7GiKrsWDBAri7u+PcuXMYNGgQli1bhuXLl+P27dto3749PvnkEwQEBODAgQNYu3YtysvL0blzZ3zxxRd49OgR4uLi4ODggM8//xwfffQRtm7dinbt2uHixYvIy8vDtm3bsH37dpw4cQJDhw7FZ599BgD47rvv8Oabb6KsrAy9e/fGtm3b8Jvf/AY+Pj6YP38+Dhw4gIqKCnz55ZdwdnZucJ2RI0fKut+qqiosXrwY6enpUKlUWLRoEXr27In09HTMmTMH7dq1w4kTJ7BhwwYcOHAAjx49wvDhw/GXv/wFKpUKo0ePxtChQ5GcnIy7d+9i69atGDlyJB49eoSFCxciMzMTgYGBePTokf6asbGxOH36NB49eoSZM2fi7bffBgD4+Phg0aJF+O677/DSSy+hY8eOePnll9GlSxcMGjSoxf/fAgCEjRk8eHCLPn//T0PFhbeHiskfHFQoIiJlZWZmmjsEvTfffFNs2LBBzJ8/X0yaNElUVlYKIYR48sknxU8//SSEECItLU1ERkYKIYQoLCwU1dXVQgghPvnkE7Fq1ao659GZP3++ePbZZ0V1dbXYt2+fcHV1FRcuXBBVVVVi0KBB4ty5c+L27dti5MiR4v79+0IIIdavXy/efvttIYQQ3t7e4sMPPxRCCLFlyxaxePFig9ep7ciRI2LAgAENviIiIhocm56eLp566in966KiIiGEEKNGjRKnT5/Wv19QUKD/fu7cueLrr7/WH6e792+++UZERUUJIYT44x//KBYuXCiEEOL8+fPCwcFBfz7duSorK8WoUaPE+fPn9ff6/vvvCyGEePTokVCr1eKnn34S1dXVYtasWWLSpEkN4jf0b8jY704+URiTcxzwlfcXB1FrePtABjKv3VP0nEGeHfDmlGDJn5s1axYcHBxw//59pKamYtasWfqflZWVAaiZ2vvss8/i+vXrKC8vNzpFc8qUKVCpVNBoNOjevTs0Gg0AIDg4GLm5udBqtcjMzMSIESMAAOXl5YiIiNB/fvr06QCAwYMH4+9//3uT8UdGRuKHH35o1r36+fnhypUrWLFiBSZNmoSxY8caPC45ORkffPABHj58iMLCQgQHB2PKlCkN4svNzQUA/OMf/8DKlSsBACEhIQgJCdGfa8+ePYiPj0dlZSWuX7+OzMxM/c+fffZZADVNUX19fdGnTx8AwNy5c/UNU1uCiYKIFOHi4gIAqK6uRseOHQ3+0l2xYgVWrVqF6OhoHD16FG+99Vaj53vssccAAG3atNF/r3tdWVkJBwcHjBkzptEOD7rPODg4NKtukpycjFdeeaXB++3bt29Q1+jUqRPOnz+PQ4cOYcuWLdizZw/++te/1jmmtLQUy5YtQ3p6Onr27Im33nqrzrTUxuIzNCspJycHf/jDH3D69Gl06tQJCxYsqHMu3X/7xj7fUkwUBmiqszCuIgnANHOHQmSUnL/8Ta1Dhw7w9fXFl19+iVmzZkEIgQsXLmDAgAEoLi6Gl5cXAGD79u36z7i6uuLePWlPRsOGDcPy5cvx73//G/7+/nj48CG0Wi369u3b6GeMXUfKE8WdO3fQtm1bzJgxA71798aCBQv05y8pqVmHpftF3qVLF9y/fx979+7FzJkzjZ73iSeewBdffIHIyEj861//woULFwAA9+7dg4uLC9zc3HDz5k0kJCRg9OjRDT4fEBCAnJwcXL58Gb17926yTVJzcdZTPcUewwEAoypNPwODyFZ98cUX2Lp1KwYMGIDg4GDs378fAPDWW29h1qxZGDlyJLp06aI/fsqUKfjqq68QGhqK48ePN+saXbt2xWeffYbnnnsOISEhGDZsGC5evGj0M3KuY0h+fj5Gjx6N0NBQLFiwAOvWrQNQU9RfunQpQkND8dhjj+E//uM/oNFo8PTTTyM8PLzJ88bGxuL+/fsICQnBBx98gCFDhgAABgwYgIEDByI4OBiLFi3SD7fV5+zsjPj4eEyaNAmPP/44vL29Zd9jbSohDOwOZMVauh9Fxj+/QXXSuwAAzQsbWKMgi5OVlYXAwEBzh0FWzNC/IWO/O/lE0ZQc+X91EBHZAiYKIiIyionCiAyFpx0SEVkjJgoiIjKKiaKeYM8OAGqmyOafTzJzNERE5sdEYcD1LjVTZHvcSeXwExHZPSYKA7wGROHHNpx+SGQI24zXeO+992R/1txdZqViojBCU52FTloOPxHVxjbjNZgoCMcca4af3G5whTZRUxYsWIBVq1YhMjISq1evxuXLlzF+/HgMHjwYI0eO1K+YPnDgAIYOHYqBAwfiqaeews2bN5Gbm4u4uDhs3LhRv2J6wYIFiI2NRWRkJPz8/HDs2DEsWrQIgYGB+nYZQE2b8YiICAwaNAizZs3C/fv3AdS03n7zzTcxaNAgaDQaXLx40eB15FqzZg0ePXqE0NBQzJkzBwDw+eefY8iQIQgNDcWLL76IqqoqVFVVYcGCBejfvz80Gg02btyIvXv36tuRh4aG1mklbrEa7StrpVraZlxc+Yf4V8pBMfmDg+LC20PF/T8NVSYwIoWwzbj524wLIYSLi4v++8zMTDF58mRRXl4uhBAiNjZWbN++vdntyFsb24ybAtuNk6VKWAPc+FHZc3pogAnrJX/MntqM15eUlIQzZ87o+zk9evQI3bp1w5QpU5rVjtzSMVHU5zsSwTgOgLOdiKSwpzbj9QkhMH/+fH1zwNqaakduDZgomiHj2j0EN/6HD5H5yPjL39Tsoc04ADg5OaGiogJOTk6IiorC1KlT8corr6Bbt24oLCxESUkJXFxcmmxHbg1YzG4uNgckajZbbzMOADExMQgJCcGcOXMQFBSEtWvXYuzYsQgJCcGYMWNw/fr1ZrUjt4ZiNtuMG5JzHFP23sN7j96FpjoL1wIXw3PINNYpyCKwzTi1FNuMK4hTZImImCiM+sO9KFxtx7/ciMi+MVE0YlRNvQ13y8wbBxGRubUoUTx48ABVVVVKxWJRJngDms6/vs64dg8nLheYLyCiWmystEitSM6/HUmJorq6Gjt37sSkSZPQrVs3BAQEoEePHggODsZrr72G7OxsyQFYJBatyYI5OzujoKCAyYIkE0KgoKAAzs7Okj4naR1FZGQknnrqKaxbtw79+/dHmzY1eaawsBDJyclYs2YNpk2bhrlz50oKwpI9qAQKS399feJyASJ6d278A0QmplarodVqcfv2bXOHQlbI2dkZarVa0mckJYrvv/8eTk5OyMvL0ycJAHB3d8eMGTMwY8YMVFRUSArAko3yAnCxpk7BYg5ZCicnJ6OtL4iUJun3n5OTEwBg2rRpDX6WlpZW5xhbMMEbcOHadSKyc5ISxZ49e7BmzRqUlJQgKyurTiE7JiZG8sUTExPRr18/+Pv7Y/36hq0IiouLMWXKFP3qzm3btkm+BhERtYykv5dHjBiB0tJSfPrpp1i1ahUuXbqEjh07wtPTE+3atZN04aqqKixfvhyHDx+GWq1GeHg4oqOjERQUpD9my5YtCAoKwoEDB3D79m3069cPc+bMQdu2bSVdi4iI5JOUKLy8vDBv3jz07t1b39q3sLAQOTk5CAgIkHThU6dOwd/fH35+fgCA2bNnY//+/XUShUqlQklJCYQQuH//Ptzd3eHoyLEgIqLWJGnoSTcdT5ckgJpC9uDBg/Uthps7ZS8/Px89e/bUv1ar1cjPz69zzEsvvYSsrCx4enpCo9Fg06ZNdYroOvHx8QgLC0NYWFirzAThegoisieSEkVkZCQ++ugj/Pzzz3XeLy8vx5EjRzB//vw6rYONMZRQVCpVndeHDh1CaGgorl27hh9++AEvvfSSwRbBMTExSE9PR3p6Orp27SrhjhoXPGKS/nvd3tkdbqQpcm4iImsiKVEkJibCwcEBzz33HDw9PREUFARfX1/06dMHu3btwiuvvFJnP1tj1Go1rl69qn+t1Wrh6elZ55ht27Zh+vTpUKlU8Pf3h6+vb5NthJXGxoBEZO8kDfg7Oztj2bJlWLZsGSoqKnDnzh20a9cOHTt2lHzh8PBwZGdnIycnB15eXti9ezd27txZ55hevXohKSkJI0eOxM2bN3Hp0iV9TaO1HHKKQlhpKrx+WXTX4UYa7nkM0w8/cfEdEdk6o08UmZmZdVZZR0VFISMjA0DNeonTp09j8+bNOHXqlOQLOzo6YvPmzRg3bhwCAwPxzDPPIDg4GHFxcYiLiwMAvPHGG0hNTYVGo0FUVBTef//9OpudtAZdc8Cej2qGn4iI7I3RjYt69OiBEydOwMfHBwDQr18/XLp0CQCQmpqKCRMm4Nlnn0VKSgr+53/+x+BCvNamyMZFvzhxuQAdbqQh+VgSXirfigedApEb9gbueQzTH8MnCiKyBbI3Lvruu+/wu9/9Tv+6Q4cO+u937NiBpUuXIj4+HkePHsX777+vULiW55BTFH5sY3hfihOXC/RfRES2yGii0Gg0+OKLL/Sv/f39sXfvXty6dQv79u3D1KlTAQDdunVDWRk3biAiskWSZj1t3LgRf/nLX+Dl5YVBgwZh+PCaGUEVFRW4f/++SQIkIiLzkjTrycPDA4cPH0Z1dXWdhW/JycmIjIxUPDgiIjI/Wf0w6q+OHjt2LMaOHatIQJbqQSVwpRg4lAeM8DB3NERErYfbLDTDKK+aduMPKoFj+Y0fx4I2EdkiJopmmOAN+Ln9ujcFW3kQkT1homhC7TUTAPBjAZCQZ6ZgiIjMQFKNwtXVtUHjPqCmwZ9KpTLYsM+WdHwMwMOa4acRTR5NRGQbJCWKkpISU8VhFdydAc0vC7F1w0/1nziIiGyN7F2AioqKkJ2djdLSUv17TzzxhCJBWYqI3p1ZoCYiuycrUXz66afYtGkTtFotQkNDkZaWhoiICBw5ckTp+IiIyMxkFbM3bdqE06dPw9vbG8nJyTh37pxiGwZZG86AIiJbJytRODs7w9nZGQBQVlaGgIAAfVdZW8PusERk72QlCrVajbt37+Lpp5/GmDFjMHXq1Aa709kil6IsjKtIMjpFljUNIrI1smoUX331FQDgrbfeQmRkJIqLizF+/HhFA7M0xR7D4VKUhakOqfgDonAsv2YhHhGRrZP1RLFx40ZotVoAwKhRoxAdHY22bdsqGpilKVJH4UGnwDpTZImI7IGsRHHv3j2MGzcOI0eOxJYtW3Dz5k2l47J4V4qBNamGh6C4kRER2RJZieLNN99ERkYGtmzZgmvXrmHUqFF46qmnlI7NIrkUZWGFSxL83GqSxbF8znwiItvWol5P3bp1g4eHBzp37oxbt24pFZPFKvao2ahpREUq1g+vaRRIRGTrZCWKjz/+GKNHj0ZUVBTu3LmDTz75BBcuXFA6Noujq1PUphuCeudgBpKy6g7BcQiKiGyBrFlPeXl5+NOf/oTQ0FCFw7Euo7x+/T6v4CGAO4gK7G62eIiITEFWoli/fr3ScVi0ex7DDNYhJnj/OkX21fT2rRwVEVHrkJQoHn/8caSkpDRoN24vbcblqj38xJXeRGRtJCWKlJQUAPbXbjyid2dk3DB+jGN5CSrburZOQERErUj2grv8fCObR9sg7jtBRPZKVo3i3r17GDt2LNzd3TF79mzMnDkT3buziNscHIYiImvDBXcyuBRloZM2qcH7WddLGkyRJSKydlxwJ5Fu0Z3bjdQ67+umyv7z8p3WDomIyKS44E4iQ4vugJppsoE9pBWzuRiPiKyB5BqFEALp6elccEdEZCckP1GoVCqcO3fO7pJERO/OCPbsoH/tXJIHn/R3DdYqiIhsiayhp4iICJw+fVrpWKxGscdwlLp6w7kkr0GtgojI1siaHpucnIy4uDj4+PjAxcVFvzLbnuoUReoo+KS/2+Jz6eoUnCpLRJZKVqJISEhQOg6rpZsqW6SOMncoREQmIStRbN++3eD7v//97yWdJzExEf/5n/+JqqoqLFmyBGvWrGlwzNGjR/Hyyy+joqICXbp0wbFjx+SErJhgzw7IuFbT00q3j7bbjVQmCiKyWbJqFC4uLvovBwcHJCQkIDc3V9I5qqqqsHz5ciQkJCAzMxO7du1CZmZmnWPu3r2LZcuW4euvv0ZGRga+/PJLOeGaTP2pso7lJbIX3XGqLBFZKllPFP/1X/9V5/Wrr76K6OhoSec4deoU/P394efnBwCYPXs29u/fj6CgIP0xO3fuxPTp09GrVy8ANQv8LNkoL+DHgppFd9yXgohsRYtWZus8fPgQV65ckfSZ/Px89OzZU/9arVY3aDT4008/oaioCKNHj8bgwYOxY8cOJcI1GTmL7oiILJ2sJwqNRqPfj6Kqqgq3b9+WXJ8QQjR4r/YeFwBQWVmJM2fOICkpCY8ePUJERASGDRuGvn371jkuPj4e8fHxAIDbt29LikNpjuX21YKdiGyfrERx8ODBX0/g6Iju3bvD0VHaqdRqNa5evap/rdVq4enp2eCYLl266OshTzzxBM6fP98gUcTExCAmJgYAEBYWJvV2Wqz+zKcrxTV7aI/o3YVDUERk9WQNPZ06dQru7u7w9vbGtm3b8Mwzz+Ds2bOSzhEeHo7s7Gzk5OSgvLwcu3fvblDnmDp1Ko4fP47Kyko8fPgQJ0+eRGBgwz5L5lS/SeAoL8DPrWYPbTYIJCJbICtRvPvuu3B1dUVKSgoOHTqE+fPnIzY2VtI5HB0dsXnzZowbNw6BgYF45plnEBwcjLi4OMTFxQEAAgMDMX78eISEhGDIkCFYsmQJ+vfvLydkRdVu5VF/5tMEb2D9cMC7s/Q9tE9cLuDsJyKyOCphqFjQhIEDB+LcuXN4/fXXodFo8Pzzz+vfM7ewsDCkp6eb7gI5xwFAv5YCgH6Fdm7YG/r31qTWDEH5udU8ZYwY2vwd8rhKm4ham7HfnbKeKLy8vPDiiy9iz549mDhxIsrKylBdXd2iIG2NbgjqSjFwzL52jSUiGyMrUezZswfjxo1DYmIiOnbsiMLCQmzYsEHp2CxasGeHOkNQ9emGoPzcWjEoIiITkDXrqX379pg+fbr+dY8ePdCjRw/FgrJoviP1w0/NdaUYeO/vaQCAoYG+Tc6EYqNAIrIkiiy4o8b30dYNQQE1CYMzoYjI2jBRKKCxfbSBX4egOAxFRNZKUqJ44YUXAACbNm0ySTDWKNizAzyHTAO6m3/aLhGRKUiqUZw5cwZ5eXn461//innz5jVow+Hu7q5ocLYor+Ah3jmYAQBcuU1EVkFSoli6dCnGjx+PK1euYPDgwXUShUqlktwY0N6M8gIqb9csxMsreAiAXWaJyPJJGnpauXIlsrKysGjRIly5cgU5OTn6LyaJpk3wBv4QVoLfTw6WtXKbiMgcZE2P/fjjj3H+/HkcP14zTfSJJ55ASEiIooHZA90mR3yqICJLJmvW04cffog5c+bg1q1buHXrFubMmYOPPvpI6dhs2ojeXQBwuiwRWT5ZTxSffvopTp48CRcXFwDA6tWrERERgRUrViganLVxaevQrOM63EhDVOAwJgkisgqyniiEEHBw+PWXooODg8GNiGyW70jFTpV1vQTvHMwwuM82u8kSkSWQ9USxcOFCDB06FNOmTQMA7Nu3D4sXL1Y0MHugG37Kul6zKx5rFURkiWQlilWrVmH06NFISUmBEALbtm3DwIEDlY7N5kUFdkdUYHf9ugoiIkskK1EAwKBBgzBo0CAlY7ErHW6k4Z7Hr3tUGJsBdeJyARsEEpHZsNeTXL4jG9Yqbv7LYGPApnAGFBFZMiYKpfiNBmC4MWBTogK7I7CHq8IBEREpQ1ai2Lx5M4qKipSOxbr1HS+5MWCHG2kmCoaISDmyEsWNGzcQHh6OZ555BomJifY1Nba+esNPziV58El/V9YQlDGcKktE5iIrUaxduxbZ2dlYvHgxPvvsM/Tp0wf//d//jcuXLysdn3XxG41SV284l+TJGoIiIrJEsmsUKpUKHh4e8PDwgKOjI4qKijBz5kz89re/VTI+69J3PFyiN6DU1dvckRARKUbW9NgPP/wQ27dvR5cuXbBkyRJs2LABTk5OqK6uRp8+ffDBBx8oHafV0Q1B6RR7DEeROqrBcb/WKVz1e1VwnwoisiSyEsWdO3fw97//Hd7edf9ybtOmDQ4ePKhIYNZMtzWqjnNJHgAYTBQ6NVNk7zS5T4WuTsF1FUTUWmQNPZWVlTVIEqtXrwYABAYGtjwqK1ekjkJu2Bv6r+YMRUUFduc+FURkkWQlisOHDzd4LyEhocXBWK16M5+CPTsg2LODmYIhIlKWpETx8ccfQ6PR4NKlSwgJCdF/+fr6cuOiJkiZNqurVRjqKKvDqbJE1Fok1Sief/55TJgwAa+//jrWr1+vf9/V1RXu7u6KB2crdDULJWsVREStRVKicHNzg5ubG3bt2mWqeGyKbvgpA1EoUkfVmQXVGHaUJSJLIylRPP7440hJSYGrqytUKpX+fSEEVCoV7t27p3iA9qJ2O4/aXWWJiMxNUqJISUkBAJSUlJgkGCIisjyy96Og5tPPgGrrgAflVeYNhohIIkmJQjfkZKgJoN0PPfmOBHKON3mYS1EWOmmTjBa0iYgsiaREwSGnFvIbDdz8F9xupDJREJHVkLSO4vHHHwdQ82TRoUOHBl/UhF/2rHBp66DIgjy2Hiei1iApUdQuZt+7d6/Bl1SJiYno168f/P3966zLqO/06dNwcHDA3r17JV+DiIhaxmxboVZVVWH58uVISEhAZmYmdu3ahczMTIPHrV69GuPGjTNDlCZy81/AT4mN/lg3VbY5K7SJiExNVqIoLS3F//7v/2L69OmYMWMGNm7ciNLSUknnOHXqFPz9/eHn54e2bdti9uzZ2L9/f4PjPvroI8yYMQPdunWTE6rl+WVvbVw5anT4aUTvLvDu3B55BQ/xz8t3jJ6Sw09EZEqyEsW8efOQkZGBFStW4KWXXkJWVhZeeOEFSefIz89Hz5499a/VajXy8/MbHPPVV19h6dKlcsK0TM3cW5vdZInIUshaR3Hp0iWcP39e/zoyMhIDBgyQdI7GptjW9vLLL+P999+Hg4OD0XPFx8cjPj4eAHD79m1JcZiNbvjpN8ObPpaIyIxkPVEMHDgQaWm/tpw4efIkRowYIekcarUaV69e1b/WarXw9PSsc0x6ejpmz54NHx8f7N27F8uWLcO+ffsanCsmJgbp6elIT09H165dpd2MOdQafmpMhxtp6HAjDY7lJci6XoJ/nkyr0+aDiKi1SHqi0Gg0UKlUqKiowI4dO9CrVy8AwM8//4ygoCBJFw4PD0d2djZycnLg5eWF3bt3Y+fOnXWOycnJ0X+/YMECTJ48GU8//bSk61ikvuP1SULfOPCa4Vljo7yAHwuAY/nABCP7H524XMBd74jIJCQlCiW3OXV0dMTmzZsxbtw4VFVVYdGiRQgODkZcXBwA2FZdojG64ae+4xs9ZIJ3TZIgIjIXSYmi9vanRUVFyM7OrjPbqf72qE2ZOHEiJk6cWOe9xhLEZ599JuncFu+XVdq4ctRootC5UgysSQWGBt7kHhVE1KpkFbM//fRTbNq0CVqtFqGhoUhLS0NERASOHDmidHzWpZn9ngDUGX4CaoagjA0/ATXJAlk5mNYph63IiajVyCpmb9q0CadPn4a3tzeSk5Nx7tw56ygiW6kJ3sD64YCfm/Hj2NKDiExBVqJwdnaGs7MzAKCsrAwBAQG4dOmSooHZjcIcIHGN0ZXaRETmJGvoSa1W4+7du3j66acxZswYdOrUqcHUVrvlO7Lmf5szBKWbJlv4y+yukKbXVLBWQUStTVai+OqrrwAAb731FiIjI1FcXIzx45suyFI9fcfXfCWuadbhhmoVBvWepEx8RESQmShKS0vx5z//GSkpKVCpVHj88cdRXV2tdGxUzwTvmq81qeaOhIjsiaxEMW/ePLi6umLFihUAgF27duGFF17Al19+qWhwVk03BKXT3NlQREQWxmy9nqihplZp1/ZjAZCQZ3y1NhGREmQlCl2vp2HDaubyy+n1RPI12dbD2NNL/ScdIqImmK3Xk91pajbUzX/pi9qd3Ica3VO7qbYe9Z9IlNh2lYjsl9l6PVEtummyAFCYA08AnkOmNWsIiojI1GT3ejp//jyOH6/563jkyJGsUbSEbpos0OypskRErUV2C485c+bg1q1buHXrFubOnYuPPvpI6djIFHKON/wiIjJCVjF769atOHnyJFxcXAAAq1evRkREhH66LFkWY0NYrF8QUVNkPVEIIepsT+rg4GBwa1MyLV07j4Q8c0dCRLZM1hPFwoULMXToUEybNg0AsG/fPixevFjRwGyWlFbkRtRp5wGupyAi05GVKFatWoXRo0cjJSUFQghs27YNAwcOVDo2+6Xb+e43jTcJVLSdR3MSF9dfENktyYlCCAGtVotBgwZh0KBBpojJvtXe+a4Z3WRbylD9gnULIqpNco1CpVLh6aefNkEoBKBmmmz3/pI+omvnQURkCrKK2cOGDcPp06eVjsV+KDiMo6tVGFupTUTUErJqFMnJyYiLi4OPjw9cXFwghIBKpcKFCxeUjs9+FeYg+MJ7AIAH5VX6t4s9htdp79FUOw/F1K5jsF5BZFdkJYqEhASl46Daarf0qMW5pGZ8yVAfKCW7ybJXFBHVJitRdO/evcHGRbGxsUrHZtuMTZOt3dIDQO4vv7h90t81eHiT3WSJiFpAVo1i3rx5yMjIwIoVK/DSSy8hKysLL7zwgtKx0S+a+ot+gjeg6dxKwRCR3eHGRVbGuSQPPunvNqhVEBGZCjcuMieJq7SLPWrWVTRWq9C19BjlZeIhKDkry1kAJ7JashLFyZMnG2xcFBgYqN/YiLOflBfs2QHwnIaMa1EGaxWmbOnBpoJE9k1WokhMTFQ6DpJINwSlEwsgti1wpS1wrGI4AA5LEZEyZCWK2hsYUQs1tUWqAbohKIOnq84DKgEmCiJSiqxEQeZTewjKkOokw1NozU6JDZJY5yAyCyYKapHa9QvWK4hsk2KJ4saNG/Dw8FDqdNQCvtV5UDWyOI/TaolIKlkL7gzhxkWtq7G/3o85Dkem8MaVYui/CktrfuZckge3G0psYEFE9kSxJ4pvvvlGqVNRC5T6ROHd/F+fGK4UA37tgPVhjbcAISIyRtYTxerVq5v1HkngO/LXrxaY4A2sH/7rl5+bQvFZgpzjDb+IyORkPVEcPnwY77//fp33EhISGrxHlqf++gtAuboFd8sjsk2Snig+/vhjaDQaXLx4ESEhIfovX19faDQayRdPTExEv3794O/vj/Xr1zf4+RdffKG/xvDhw+v0l6Lm07X22PZgODKr69YvHIpZtyAi4yQ9UTz//POYMGECXn/99Tq/2F1dXeHu7i7pwlVVVVi+fDkOHz4MtVqN8PBwREdHIygoSH+Mr68vjh07hk6dOiEhIQExMTE4efKkpOvYO11rDwA45BSFQ0516xe7274LPzPERUTWQ1KicHNzg5ubG6ZPnw53d3e4urpi7dq1OHv2LN544w0MHDiw2ec6deoU/P394edX82tq9uzZ2L9/f51EMXz4ryuQhw0bBq1WKyVcQk3NorG+T2tSATxq1XCIyArJqlG8++67mDVrFlJSUnDo0CG8+uqrWLp0qaS/9vPz89GzZ0/9a7VabfTzW7duxYQJEwz+LD4+HvHx8QCA27dvNzsGa6cb/zfWtK85DNUt6rPY9RfGCtpcyU2kCFmznhwcHADUTImNjY3F1KlTUV5eLukcQogG76lUKoPHJicnY+vWrY0Wy2NiYpCeno709HR07dpVUhz27v9VNKxb6NZd6LRk/UXGtXv6LyKyTrKeKLy8vPDiiy/i8OHDWL16NcrKylBdXS3pHGq1GlevXtW/1mq18PT0bHDchQsXsGTJEiQkJKBzZzvZxk3iPhVyjfICjiEKZ2B43YUO118Q2TdZiWLPnj1ITEzEq6++io4dO+L69evYsGGDpHOEh4cjOzsbOTk58PLywu7du7Fz5846x/z888+YPn06/va3v6Fv375yQrULwZ4dZP3Fbqh+saaRB4faw1MWOwxFRCYha+ipXbt2ePDgAXbt2gUAqKioQMeOHSWdw9HREZs3b8a4ceMQGBiIZ555BsHBwYiLi0NcXBwA4J133kFBQQGWLVuG0NBQhIWFNXFWMoVij+Eoda3JKFbVBoQL84gUoRKGigVNiI2NRZs2bXDkyBFkZWWhqKgIY8eOxenTp00RoyRhYWFIT083dxim08gvPaVqAGtSfxl+qrWiu/bWqrqnitywNxS5Xm0mXZzHwjaRUcZ+d8p6ojh58iS2bNkCZ2dnAECnTp0kF7NJWUr9kh3lVTdJXCkGjuUrcmoislKyahROTk6oqqrSz1K6ffs22rRRrBEtmVH9uoWhmkVzptPWx7oGkfWSlShWrlyJadOm4datW/jd736HvXv3Yu3atUrHRobI2Dq1pXQtQABgXMVwjKoGUFzzuuNjgLuz8c87l+QBABMFkZWSlSjmzJmDwYMHIykpCUII7Nu3D4GBgUrHRhIptQCvttotQIC6bUAMTaU1hNNriayb7P0oAgICEBAQoGQsZIGabAGiIFMsytPXbkz9BMZiOdkw7pltrYz9YrpmeZtIyalr1Mc6B5F5MFFQi9SuX9SeRltbscfwhm9KxDoHkfkwUZBstesXV34pbhtKFEXqqBb/gmedg8h8mChs0D2PYY3+rMONNMWuU7t+oXS9gogsBxOFDYroXbd54onLBWaKRFmm3Ma1xSytTQiL66QgJgpSTO16hU5jdQupDNU5mlO3qD2Tivt3E8nDRGFnDA1LKTEcVX+9BWC8biGVoToH6xZErYOJwg7UH4qqTalhKSkty4nIujBRkFUztj7DYuoXRFaOiYKMzpLSkTs81Zx1FnIZW59h9+suLK24TqZl4skLTBR2ztCwlFLDUc1dZyGXsfUZhp4yGmsRwiI3kXFMFGQyXGdBZBuYKKhZag9PKbloj4gsHxMFNWCqBXuG1lnUpnQNo3ahm4VtIvmYKKhVGFpnUZvSNYzahW67L2wTtRATBUnWnFlS9U1AmtEkoHQNo3ahmwvzyObpZrmZaPYTEwU1qTUW7JmTKTZMMhXO0CJzYKIgi9FUDUOOUV5ALJTZOMkY1kDIljFRUKtoarhqaOBNVF6+g0ojxziWl0i6pq7u8Xyvlm+cZAxrIGTrmCioRWoPS7VkGCoqsDuiArsbPUbqtFzd04kSGycZwxoI2bo25g6AiIgsG58oSDGmbAcCyGuRbspeU7WZugai19ahecf5jQb6jjdpKGQ/mCjIZpm615SOseaEZlGYU/O/TBSkECYKslmt1WvK1DWQ2po1PTZxjekDIbvCREEmpVSxuzHGZlOxJxWRMpgoyG40tk7DlLULsynM4ZOFPTFxTYqJglqNsRXegPJPHLWfNkZ5GX66MGXtwmz8Rps7AmpNrVCTYqIgu2BoT2/ARvfJ6DuehWx70gpPjkwUZBcaq2VUts345efBLTo/6yFky5goyGI0NTQlR3OGs/IKHuKdgxktuo5jeYs+3mzPBZXj+aC2rXMxsh66mlT4EiBsoeKnN+vK7MTERPTr1w/+/v5Yv359g58LIbBy5Ur4+/sjJCQEZ8+eNUOUZMtG9O4C787tzR1Gs1wpBvb/u8LcYZCl8RsNuPvWJIsf95rkEmZ7oqiqqsLy5ctx+PBhqNVqhIeHIzo6GkFBQfpjEhISkJ2djezsbJw8eRKxsbE4efKkuUImK9TU9Nzm9JiyFO8czACcnQDfiKYP1u1PQLZPV5MyYa3CbIni1KlT8Pf3h5+fHwBg9uzZ2L9/f51EsX//fsybNw8qlQrDhg3D3bt3cf36dfTo0cNcYZMVa2xoyxb21GjARBvYkIUy8R8GZksU+fn56Nmzp/61Wq1u8LRg6Jj8/HwmClKUKWojptDB2cncIZCl8h0JOHc02enNliiEEA3eU6lUko8BgPj4eMTHxwMALl68iLCwMFkx3b59G127dpX1WWvFe7Y+YZ9I/4y137McdnvPW+T9/svNzW30Z2ZLFGq1GlevXtW/1mq18PT0lHwMAMTExCAmJqbFMYWFhSE9Pb3F57EmvGf7wHu2D6a6Z7PNegoPD0d2djZycnJQXl6O3bt3Izo6us4x0dHR2LFjB4QQSEtLg5ubG4ediIhamdmeKBwdHbF582aMGzcOVVVVWLRoEYKDgxEXFwcAWLp0KSZOnIhvv/0W/v7+aN++PbZt22aucImI7JZZF9xNnDgREydOrPPe0qVL9d+rVCps2bKl1eJRYvjK2vCe7QPv2T6Y6p5VwlDFmIiI6BfcM5uIiIxiovhFU+1ErNXVq1cRGRmJwMBABAcHY9OmTQCAwsJCjBkzBn369MGYMWNQVFSk/8y6devg7++Pfv364dChQ+YKvUWqqqowcOBATJ48GYDt3y8A3L17FzNnzkRAQAACAwNx4sQJm77vjRs3Ijg4GP3798dzzz2H0tJSm7zfRYsWoVu3bujfv7/+PTn3eebMGWg0Gvj7+2PlypUGlx80SpCorKwUfn5+4vLly6KsrEyEhISIjIwMc4eliGvXrokzZ84IIYS4d++e6NOnj8jIyBCvvfaaWLdunRBCiHXr1onf/va3QgghMjIyREhIiCgtLRVXrlwRfn5+orKy0mzxy/XHP/5RPPfcc2LSpElCCGHz9yuEEPPmzROffPKJEEKIsrIyUVRUZLP3rdVqhY+Pj3j48KEQQohZs2aJbdu22eT9Hjt2TJw5c0YEBwfr35Nzn+Hh4SI1NVVUV1eL8ePHi2+//bbZMTBRCCFSU1PF2LFj9a/fe+898d5775kxItOJjo4W3333nejbt6+4du2aEKImmfTt21cI0fDex44dK1JTU80Sq1xXr14VTz75pEhKStInClu+XyGEKC4uFj4+PqK6urrO+7Z631qtVqjValFQUCAqKirEpEmTxKFDh2z2fnNycuokCqn3ee3aNdGvXz/9+zt37hQxMTHNvj6HntB4qxBbk5ubi3PnzmHo0KG4efOmfk1Kjx49cOvWLQC28d/i5ZdfxgcffIA2bX79523L9wsAV65cQdeuXbFw4UIMHDgQS5YswYMHD2z2vr28vPDqq6+iV69e6NGjB9zc3DB27Fibvd/6pN5nfn4+1Gp1g/ebi4kCzW8VYs3u37+PGTNm4E9/+hM6dOjQ6HHW/t/i4MGD6NatGwYPHtys4639fnUqKytx9uxZxMbG4ty5c3BxcTFaa7P2+y4qKsL+/fuRk5ODa9eu4cGDB/j8888bPd7a77e5GrvPlt4/EwWa3yrEWlVUVGDGjBmYM2cOpk+fDgDo3r07rl+/DgC4fv06unXrBsD6/1v885//xNdffw0fHx/Mnj0bR44cwdy5c232fnXUajXUajWGDh0KAJg5cybOnj1rs/f9/fffw9fXF127doWTkxOmT5+O1NRUm73f+qTep1qthlarbfB+czFRoHntRKyVEAKLFy9GYGAgVq1apX8/Ojoa27dvBwBs374dU6dO1b+/e/dulJWVIScnB9nZ2RgyZIhZYpdj3bp10Gq1yM3Nxe7du/Hkk0/i888/t9n71fHw8EDPnj1x6dIlAEBSUhKCgoJs9r579eqFtLQ0PHz4EEIIJCUlITAw0Gbvtz6p99mjRw+4uroiLS0NQgjs2LFD/5lmkVdasT3ffPON6NOnj/Dz8xNr1641dziKOX78uAAgNBqNGDBggBgwYID45ptvxJ07d8STTz4p/P39xZNPPikKCgr0n1m7dq3w8/MTffv2lTQzwtIkJyfri9n2cL/nzp0TgwcPFhqNRkydOlUUFhba9H3//ve/F/369RPBwcFi7ty5orS01Cbvd/bs2cLDw0M4OjoKLy8v8emnn8q6z9OnT4vg4GDh5+cnli9f3mDigzFcmU1EREZx6ImIiIxioiAiIqOYKIiIyCgmCiIiMoqJgoiIjGKiICIio5goiIjIKCYKIonu3r2LP//5z/rXw4cPN8l1tFot/u///s8k5yaSgomCSKL6iSI1NdUk10lKSsLZs2dNcm4iKZgoiCRas2YNLl++jNDQULz22mv4zW9+A6CmjXtAQACWLFmC/v37Y86cOfj+++8xYsQI9OnTB6dOndKf4/PPP8eQIUMQGhqKF198EVVVVXWukZKSglWrVmHv3r0IDQ1FTk5Oq94jUR0KtiQhsgv1N5FxcXHRv+/g4CAuXLggqqqqxKBBg8TChQtFdXW12Ldvn5g6daoQQojMzEwxefJkUV5eLoQQIjY2Vmzfvr3BdcaNGyd+/PFH098QURMczZ2oiGyJr68vNBoNACA4OBhRUVFQqVTQaDTIzc0FUDOkdObMGYSHhwMAHj16pG8TXdulS5fQr1+/VoudqDFMFEQKeuyxx/Tft2nTRv+6TZs2qKysBFDT+n3+/PlYt25do+cpKCiAm5sbnJycTBswUTOwRkEkkaurK0pKSmR/PioqCnv37tVvX1lYWIi8vLw6x+Tk5Fj1xjpkW5goiCTq3LkzRowYgf79++O1116T/PmgoCCsXbsWY8eORUhICMaMGaPfrUwnICAAd+7cQf/+/U02q4qoubgfBRERGcUnCiIiMoqJgoiIjGKiICIio5goiIjIKCYKIiIyiomCiIiMYqIgIiKjmCiIiMio/w92pSbw4fLAYwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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, conf_int = kaplan_meier_estimator(\n",
    "        data_y[\"Status\"][mask_treat],\n",
    "        data_y[\"Survival_in_days\"][mask_treat],\n",
    "        conf_type=\"log-log\",\n",
    "    )\n",
    "\n",
    "    plt.step(time_treatment, survival_prob_treatment, where=\"post\", label=f\"Treatment = {treatment_type}\")\n",
    "    plt.fill_between(time_treatment, conf_int[0], conf_int[1], alpha=0.25, step=\"post\")\n",
    "\n",
    "plt.ylim(0, 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": [
    "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://scikit-survival.readthedocs.io/en/latest/api/generated/sksurv.compare.compare_survival.html#sksurv.compare.compare_survival) *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": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f3d75d069e0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAERCAYAAABl3+CQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAABEwElEQVR4nO3deVhTd/Y/8HcCKIuAAiJItIAg+04EKipoBUWLC9raugxai6hdrVbbzlTbsdV2Zro4boO1Wr91a+0o/bnQurbiggoqVkQpgqwqi+yCBO7vj0yuhISQXBJCwnk9D89Dbm7uPZdaDvd+Pud8eAzDMCCEEEI6wNd2AIQQQno2ShSEEEIUokRBCCFEIUoUhBBCFKJEQQghRCFKFIQQQhTSaqJYsGABbG1t4e3tLfd9hmHwxhtvwMXFBb6+vsjIyOjmCAkhhGg1UcTHxyMlJaXD948dO4acnBzk5OQgKSkJixcv7sboCCGEAFpOFKNHj4aVlVWH7ycnJ2PevHng8XgIDQ1FVVUVSktLuzFCQgghPXqMori4GEOGDGFfCwQCFBcXazEiQgjpfQy1HYAi8rqL8Hg8mW1JSUlISkoCAGRnZ8Pd3V3jsRFCiD7Jz89HeXm53Pd6dKIQCAQoLCxkXxcVFWHw4MEy+yUkJCAhIQEAEBwcjCtXrnA+5+Wr2yEMeAUAUJ92CWYhIzgfixBCdEVwcHCH7/XoR0+xsbHYtWsXGIbBxYsXYWlpCXt7e22HRQghvYpW7yheeuklnDlzBuXl5RAIBPjoo4/Q3NwMAEhMTERMTAyOHj0KFxcXmJqaYseOHRqPKa+8HqLcCoQNs9b4uQghRBdoNVHs3btX4fs8Hg+bNm3qpmgAHFuFiFvH8GdtETDscwD0+IkQQnr0GIU22Nbko7U4VdthEKJWzc3NKCoqQmNjo7ZDIVpmbGwMgUAAIyMjpT9DiaKtievx8NYxbUdBiNoVFRXB3Nwcjo6OcmcOkt6BYRhUVFSgqKgITk5OSn+uRw9md7fPLn2GZRatOMqvw/1PP8WDdetQd/q0tsMipMsaGxthbW1NSaKX4/F4sLa2VvnOku4o2rluBAAMYgA0ZWcDEI9TdITGL4iuoCRBAG7/DuiOoo2VI1bCTzzpCnbvv4++VLhHCFHCV199hV27dmns+H/729/g6+sLf39/REVFoaSkBIC4SM7ExAT+/v7w9/dHYmIi+5nnnnsOjx49Usv5KVF04EJuhbZDIIToAJFIhG+//RYvv/yyxs6xYsUKZGZm4tq1a5g8eTI+/vhj9r1hw4bh2rVruHbtGrZu3cpunzt3LjZv3qyW81OiIIRoXH19PSZNmgQ/Pz94e3tj//79AICUlBS4u7sjPDwcb7zxBiZPngwAWLNmDf75z3+yn/f29kZ+fj4AYOrUqQgKCoKXlxfbugcA+vXrh5UrVyIoKAjPPfccLl26hIiICDg7O+Pnn38GIB6rmT9/Pnx8fBAQEIDT/xuD3LlzJ1577TX2WJMnT8aZM2fQ0tKC+Ph4eHt7w8fHB19++aXMtZ06dQqBgYEwNBQ/yY+IiMDKlSsxYsQIDB8+HGfPnu3yz8/CwkLqZ6nM46PY2NhOSxCURWMUbXz0/26iqUWcO0/eegDN/X1ASO+SkpKCwYMH48iRIwCA6upqNDY24tVXX8WpU6fg4uKCF198Ualjffvtt7CyssLjx48hFAoRFxcHa2tr1NfXIyIiAp999hmmTZuGv/71rzh+/DiysrLwl7/8BbGxsWxd1o0bN5CdnY2oqCjcuXOnw3Ndu3YNxcXF+OOPPwAAVVVVMvucO3cOQUFBUttEIhEuXbqEo0eP4qOPPsKJEyek3q+trcWoUaPknnPPnj3w9PSU2f7BBx9g165dsLS0ZBMcAOTl5SEgIAAWFhZYu3Yte9wBAwagqakJFRUVsLbuWgExJQo5GAY49UeWUolCMtBNg9pEV3z0/24iq6RGrcf0HGyB1c97dfi+j48Pli9fjpUrV2Ly5MkYNWoUrl27BicnJ7i6ugIA5syZI3WH0JENGzbg4MGDAIDCwkLk5OTA2toaffr0wYQJE9jz9e3bF0ZGRvDx8WHvRlJTU/H6668DANzd3fHMM88oTBTOzs64e/cuXn/9dUyaNAlRUVEy+5SWlsLDw0Nq2/Tp0wEAQUFB7LnbMjc3x7Vr1zq91rY++eQTfPLJJ1i3bh02btyIjz76CPb29igoKIC1tTXS09MxdepU3Lx5k70DsbW1RUlJSZcTBT16amP1817oa9AKyV1dQ5NIuwERoieGDx+O9PR0+Pj44L333mOfsXf0CMXQ0BCtra3sa8l0zjNnzuDEiRO4cOECrl+/joCAAPY9IyMj9nh8Ph99+/ZlvxeJxP8vy+tIreh8AwYMwPXr1xEREYFNmzZh4cKFMp81MTGRmW4qObeBgQF77rZqa2vZAej2X1lZWXJjlHj55Zfx008/seeRJIGgoCAMGzZMKvE1NjbCxMRE4fGUQXcUcmQbt6KfSToA8RTZutOn0S8yUstREaIeiv7y15SSkhJYWVlhzpw56NevH3bu3Il3330XeXl5yM3NxbBhw6Sepzs6OuLw4cMAgIyMDOTl5QEQP7IaMGAATE1NkZ2djYsXL6oUx+jRo7F7926MHTsWd+7cQUFBAdzc3FBTU4PNmzejtbUVxcXFuHRJ/KSgvLwcffr0QVxcHIYNG4b4+HiZY3p4eODPP/9UKQ5V7yhycnLYO6+ff/6ZXUqhrKwMVlZWMDAwwN27d5GTkwNnZ2cA4qR4//59ODo6qhSbPJQo2hnbxMN1IwZNJjfQ7B8Mg7w/UX/xIiUKQrrgxo0bWLFiBfh8PoyMjLBlyxYYGxsjKSkJkyZNgo2NDcLDw9mxgLi4OOzatQv+/v4QCoUYPnw4AGDChAnYunUrfH194ebmhtDQUJXiWLJkCRITE+Hj4wNDQ0Ps3LkTffv2xciRI+Hk5AQfHx94e3sjMDAQgHjxtPnz57N3G+vWrZM55sSJEzF37tyu/Hg6tWrVKty+fRt8Ph/PPPMMO7vp999/x4cffghDQ0MYGBhg69at7Kqh6enpCA0NZQfZu4LHdHQvpqO6uh7Fwy/8sNQMyOUNwX8Gvgnr78TTywa9957Cz9EYBenJbt26JfMcvac5c+YM/vnPf7J3Erpk2rRp+Pzzz9m/+nuCN998E7GxsRg3bpzMe/L+PSj63UljFHL0YxpgxVShoO5PVDbUoqFJhKySaoWfUVS9TQjRb+vXr0dpaam2w5Di7e0tN0lwQY+e2ikY6AvUX8QA1MK4Nl/b4RDSa0RERCAiIkLbYXDi5uYGNzc3bYch5dVXX1XbseiOop27dkLU8Uy1HQYhhPQYlCgIIYQoRIlCBVkl1Z2OVRBCiL6hRNEJfgutCEYI6d1oMLsDrQyQWmeHmQD4pcUwSdrAvveg79Mfm1loKFtjQetrE0L0Ed1RyGHAE5eWXKm3RZ2bA1rtHeTu96SgAPUqVoYSQtSnX79+AMTrMnh7e3M6Rnx8PA4cOABAPPOqo1qCGTNm4O7du9wCVUFsbKzUtRQUFCAyMhIBAQHw9fXF0aNHAYirsiW9rTSN7ijkMOQzaGHEPWNqvR1hMv45qfedBlsCAB7IqdIkhOifmzdvoqWlhW2PoSn//e9/2eQnsXbtWrzwwgtYvHgxsrKyEBMTg/z8fAwcOBD29vY4d+4cRo4cqdG46I6iA/2YBsQgVdthEKIXOlqPwtHREe+//z7CwsIQHByMjIwMREdHY9iwYWybirq6OowbNw6BgYHw8fFBcnKywnO1tLRg+fLl8PHxga+vL/79738DELe0GDNmDIKCghAdHa1Sgdzu3bsxZcoU9nW/fv3wwQcfwM/PD6GhoXjw4IGqPxIZdXV1+OKLL/DXv/5VajuPx0NNjbjbb3V1NQYPHsy+N3XqVOzevbvL5+4M3VHI0WhkBiNRI8biMoqgnspGQnqMY6uA+zfUe0w7H2Di+g7flrcehcSQIUNw4cIFvP3224iPj8e5c+fQ2NgILy8vJCYmwtjYGAcPHoSFhQXKy8sRGhqK2NjYDjvPJiUlIS8vD1evXoWhoSEqKyvR3NyM119/HcnJyRg4cCD279+PDz74AN9++61Sl3fu3Dm89NJL7Ov6+nqEhobik08+wbvvvott27bJ/II/ffo03n77bZljmZqa4vz58zLb//a3v+Gdd96Bqal0HdeaNWsQFRWFf//736ivr5da2yI4OFjmvJpAiUKOxj4W4DU+1nYYhOgNeetRSMTGxrL71NXVwdzcHObm5jA2NkZVVRXMzMzw/vvv4/fffwefz0dxcTEePHgAOzs7uec6ceIEEhMT2WZ4VlZW+OOPP/DHH39g/PjxAMR3Hfb29krHX1paioEDB7Kv+/Tpw67GFxQUhOPHj8t8JjIyUukOsdeuXcOff/6JL7/8Umb9ir179yI+Ph7vvPMOLly4gLlz5+KPP/4An89n15vQNEoUSiioE7cQHtrPpdN9ufZ8otlSpNso+MtfUyTrURw9ehTvvfceoqKi8OGHHwKA1LoRku8lr0UiEXbv3o2ysjKkp6fDyMgIjo6OMus/tMUwjMzdBsMw8PLywoULFzjF337NibZrX3S05oQqdxQXLlxAeno6HB0dIRKJ8PDhQ0RERODMmTPYvn07UlJSAABhYWFobGxEeXk5bG1t1bbeRGdojKKdBqunSxA+ZgyRWif/rxZCiPJKSkpgamqKOXPmYPny5cjIyFD6s9XV1bC1tYWRkRFOnz6Ne/fuKdw/KioKW7duZX95V1ZWws3NDWVlZWyiaG5uxs2bN5WOgcuaE5I7ivZf8h47LV68GCUlJcjPz0dqaiqGDx+OM2fOAACGDh2KkydPAhB3fW1sbGTvbu7cucN5tpcqKFF0INu4Fb+bN+NKva3Me5IK7YYmEZ4UFODBunWoa7OGLSFE2o0bNzBixAj4+/vjk08+Uem5+uzZs3HlyhUEBwdj9+7d7KI9HVm4cCGGDh0KX19f+Pn5Yc+ePejTpw8OHDiAlStXws/PD/7+/nJ/YXdk0qRJ7C/u7vavf/0L27Ztg5+fH1566SXs3LmTvZs5ffo0Jk2apPEYaD2Kdnamn0B24U84Up0K98c8VJWuRKKLeMZB+0dPhpfOwfLWdTwpKECfoUM7XbNCEXr0RDRJF9aj6MkeP36MyMhInDt3DgYGBtoOhzV69GgkJydjwIABKn2O1qNQg0AzD/g1d76faMRIDHrvPfQZOlTzQRFCtMbExAQfffQRiouLtR0Kq6ysDMuWLVM5SXBBg9mEEKKE6OhobYcgZeDAgZg6dWq3nKtLdxT19fVoaWlRVyw9Tj88liq6K6j7k50BpW60Qh4hpKdSKVG0trZiz549mDRpEmxtbeHu7g57e3t4eXlhxYoVyMnJ0VSc3a7RyAwAMJLJwNY/LbA5zxYXKs1k9qO244QQfadSooiMjERubi7WrVuH+/fvo7CwEA8fPsTZs2cRGhqKVatW4fvvv9dUrN2qsY8FGnnG6MsT3zEVN/bB1WrZREEIIfpOpTGKEydOwMjICPfu3QOf/zTHWFlZIS4uDnFxcWhuVmIUWEcY8Foh6FOHtwZk4otympVECOmdVLqjMDIyAgBMmzZN5r2L/2u3LdlHH/xpxMObAw2QZljQ6b5N2dlUS0FIF+3cuROvvfaaRs/BMAzGjh3LNtpTt8LCQkRGRsLDwwNeXl74+uuv2fdefPFF+Pv7w9/fH46OjvD39wcgrjOJj4/XSDzqoFKi+OGHH7Bq1SrU1tbi1q1bUgPZCQkJKp88JSUFbm5ucHFxwfr1sm0Fqqur8fzzz8PPzw9eXl7YsWOHyufgyttkGFyaGdw1BLL5isdezEJDAYDWpiBEBxw9ehR+fn6wsLDQyPENDQ3xr3/9C7du3cLFixexadMmZGVlAQD279/PVmjHxcVh+vTpAMR9roqKilBQ0PkfpdqgUqIYOXIkPD098ejRIyxbtgyurq4IDAzE5MmTVe430tLSgqVLl+LYsWPIysrC3r172R+mxKZNm+Dp6Ynr16/jzJkzeOedd/DkyROVzsNVoJkH3q11htuTJ+jTKt0gsP3MpwK3QLQ4uaChScRWbbf9IoSIW2IHBQXBy8sLSUlJ7PYdO3Zg+PDhGDNmDM6dO8duLysrQ1xcHIRCIYRCIfvemjVrsGDBAkRERMDZ2RkbNjxdffKLL76At7c3vL298dVXX8mNo23L8Pz8fHh4eODVV1+Fl5cXoqKi8Phx1xqC2tvbIzAwEABgbm4ODw8PmfoLhmHwww8/SHWkff7557Fv374unVtTVBqjcHBwwLx58zBs2DB2oYzKykrk5eV1Wlbf3qVLl+Di4sIuBDJr1iwkJyfD0/NpryUej4fa2lowDIO6ujpYWVmxHSG7Q6GJJ5r5T1e0Km7sg815tgiwrMfQfgo+SEgP9tmlz5Bdma3WY7pbuWPliJUK9/n2229hZWWFx48fQygUIi4uDk+ePMHq1auRnp4OS0tLdiU3AHjzzTfx9ttvIzw8HAUFBYiOjsatW7cAANnZ2Th9+jRqa2vh5uaGxYsXIzMzEzt27EBaWhoYhkFISAjGjBnDHk/i3Llz+M9//sO+zsnJwd69e7Ft2za88MIL+OmnnzBnzhypz+zevRv/+Mc/ZK7JxcWFXR1Pnvz8fFy9ehUhISFS28+ePYtBgwbB1dWV3RYcHIz169fj3XffVfhz1AaVfutKujK2XU3JysoKVlZWMvt0pri4GEOGDGFfCwQCpKWlSe3z2muvITY2FoMHD0ZtbS32798vNYgukZSUxP6FUlZWpsolKS3Asl4cd2MfAMCLGjkLIfprw4YNOHjwIADxc/ycnBzcv38fERERbJO7F198EXfu3AEgnjzT9ilDTU0NamtrAYh7L/Xt2xd9+/aFra0tHjx4gNTUVEybNg1mZuLZidOnT8fZs2dlEkVlZSXMzc3Z105OTuxYQVBQkEybb0Dcb2r27NkqXW9dXR3i4uLw1VdfyTzm2rt3r9TdBIBuaxnOhUqJIjIyEnFxcZgyZQqGtmlb8eTJE6SmpuK7775DZGSkUoMy8lpMtU8wv/zyC/z9/XHq1Cnk5uZi/PjxGDVqlMwPPSEhgR0jCQ4OVuWSlBZmVY8wq3pszpNtEkiILunsL39NOHPmDE6cOIELFy7A1NQUERERbNvujv6wbG1txYULF+Q+1m7bjlzS5lvZtnWGhoZobW1l/+hsfyx5j55UvaNobm5GXFwcZs+ezY5DSIhEIvz3v/9Fenq61PbuahnOhUpjFCkpKTAwMMBLL72EwYMHw9PTE05OTnB1dcXevXvZFaqUIRAIUFhYyL4uKiqSWuIPED+7nD59Ong8HlxcXODk5ITsbPXeMsvTttU4IaTrqqurMWDAAJiamiI7O5udJRkSEoIzZ86goqICzc3N+PHHH9nPREVFYePGjezrzhYBGj16NA4dOoSGhgbU19fj4MGDUgskSbi5ueHu3btyjtCx2bNny20ZLi9JMAyDV155BR4eHli2bJnM+ydOnIC7uzsEAoHU9u5qGc6FSncUxsbGWLJkCZYsWYLm5maUl5fDxMQE/fv3V/nEQqEQOTk5yMvLg4ODA/bt24c9e/ZI7SPpwz5q1Cg8ePAAt2/f1vji5l4OlgCAvEqNnoaQXmXChAnYunUrfH194ebmhtD/zRS0t7fHmjVrEBYWxg4CS2ZTbtiwAUuXLoWvry9EIhFGjx7NrqMtT2BgIOLj4zFihLjmaeHChTKPnYCnLcNdXDpfiIyLc+fO4f/+7//g4+PDPtL69NNPERMTAwDYt2+fzGMnoPtahnOhsM14VlYWPv30U7baety4cdiwYQO8vLwAAD///DMyMzMRFRXF/sdRxdGjR/HWW2+hpaUFCxYswAcffMD+Q0hMTERJSQni4+NRWloKhmGwatUqmUGm9rraZvzy/csAgLwbF2CWU4y9hocBANNMxHO7JY+e1vtIP/4ySdoAg7w/0TjtRYhGjERHPAdbdvgetRonmkJtxp8qLS3FvHnz5C5fqi1NTU0YM2YMUlNTu2XCjqptxhVGNG7cOKmlA4uKitgkcf78ecydOxcvvvgi4uPj8cknn8gtxFMkJiaGzbISiYmJ7PeDBw/Gr7/+qtIxNcW4Nh+N5o4dvt/sHwSDvD9hdC1dYaJQRNIYkBIGIZpjb2+PV199FTU1NRqrpVBVQUEB1q9f362zOlWhcIzi119/xQcffMC+bvtD3bVrFxITE5GUlIQzZ87gs88+01yUOkA0YiRanDRzK0sIUa8XXnihxyQJAHB1dUVERIS2w+iQwkTh4+OD3bt3s68lI/wPHz7EoUOH2KIVW1tbNDU1aTZSQgghWqHSrKcvv/wS//nPf+Dg4IDAwEA8++yzAMRTwerq6jQSYHcT2gnZ7+tdHbQYCSGE9AwqPRCzs7PD8ePHpeYgA+LR+sjISLUHp4+Uaekh7HQPQgjpPpxWuGtfHR0VFSXVu4UQQtrr16/7+94cOnQIH3/8scaO/8UXX8DT0xO+vr4YN24c7t27B0D8x7OkS6y/vz+MjY1x6NAhAOJ2Rbq2yFuXlkIlhBBNYBgGra2tXT7O559/jiVLlqghIvkCAgJw5coVZGZmYsaMGWyfpsjISLYo79SpUzA1NUVUVBQAYPHixfj88881FpMm9My5WDqgbQfZof3UO9tJ0frZNHWW6Lq6ujpMmTIFjx49QnNzM9auXYspU6YgPz8fEydORGRkJC5cuIBDhw5h165d2L17N4YMGQIbGxsEBQVh+fLlyM3NxdKlS1FWVgZTU1Ns27ZNpjHpnTt30LdvX9jY2AAA4uPjYWFhgStXruD+/fv4/PPPMWPGjC5dS9tH7qGhoXJX+Dxw4AAmTpwIU1NTAMCoUaMQHx8PkUjUY6fDtqcbUWpBg5UnTCvFDcn+NOLhP4YX4cXUA7Blu8gCkOkkyy8thknSBjlHVKzZP4hz/QUhusTY2BgHDx6EhYUFysvLERoaitjYWADA7du3sWPHDmzevBlXrlzBTz/9hKtXr0IkEiEwMBBBQUEAxP3dtm7dCldXV6SlpWHJkiU4deqU1HnOnTvHtvuWKC0tRWpqKrKzsxEbGys3UYwaNYptPtjWP//5Tzz33HMdXtf27dsxceJEme379u2TauXB5/Ph4uKC69evs9fT06mUKMzNzeU28JJ0jNXUilHa9FxDK2DKx+0+NWhlchBg+fQv+vadZJv9g8BlfT9+aTGMAEoUpFvc//RTNN1Sb8+0vh7usHv/faX2ZRgG77//Pn7//Xfw+XwUFxfjwYMHAIBnnnmGbe+RmpqKKVOmsI3ynn/+eQDiO5Lz589j5syZ7DHlTc8vLS1lu9JKTJ06FXw+H56enuw52zt79qxS19HW999/jytXruC3336TieHGjRuIjo6W2i7pFKuXiUJeltVHQjshbhafAADEtPZBfEk+XhC4osHoaRdZADKdZEUjRnL6Zc/lDoQQXbV7926UlZUhPT0dRkZGcHR0ZDvJSlqEA/I7TAPirrL9+/fvtEmgiYkJqqulZxm27RTb0fFVvaM4ceIEPvnkE/z2229SxwfEq4JOmzZNZonontwpVh7Oj54ePXqEnJwc9j8wIO7eqG8KBvrCtiYffVofo0HbwRCiBsr+5a8p1dXVsLW1hZGREU6fPs3OFGovPDwcixYtwnvvvQeRSIQjR47g1VdfhYWFBZycnPDjjz9i5syZYBgGmZmZ8PPzk/q8h4eH3DGDzqhyR3H16lUsWrQIKSkpsLWVXYJg7969WLduncz2O3fusO2QdAGnWU/ffPMNRo8ejejoaKxevRrR0dFYs2aNmkPTLkkX2bt2QlQY2QMA+C2Nij5CCFHC7NmzceXKFQQHB2P37t0dro4pFAoRGxsLPz8/TJ8+HcHBwbC0FP9/uXv3bmzfvh1+fn7w8vJCcnKyzOdHjx6Nq1evKr1OBRcrVqxAXV0dZs6cCX9/f3asBRCvbldYWIgxY8ZIfebBgwcwMTGBvb29xuJSN053FF9//TUuX76M0NBQnD59GtnZ2Vi9erW6YyOE6BFJ9wYbGxupZqNt/fHHH1Kvly9fjjVr1qChoQGjR4/GO++8A0C8Kl1KSorC85mamuK5557DyZMn8dxzz2Hnzp1y4+mKEydOdPieo6OjzFrZALBnzx4sWrSoy+fuTpzuKIyNjWFsbAxAPIjk7u6O27dvqzUwXZHbYIz9BaXddr76tEvsFyH6LiEhAf7+/ggMDERcXJzMLKbOvP/++2ho6FkPjfv374+//OUv2g5DJZzuKAQCAaqqqjB16lSMHz8eAwYMkFmdrjcIsKxHboMxrlab0RrahGhA+8XMVDVo0CCpx0E9wfz587Udgso4JQrJAulr1qxBZGQkqqurMWHCBLUG1hOV8muwo/UQfHiuCOZ5IcyqHlerzTr/oBIM8v6E4aVzEI0YiaySaoULHBFCSHfi9Ojpyy+/RFFREQBgzJgxiI2NRZ8+fdQaWE/zXEMr7FstcB/luMGot09Ls794LrXRtfRO9iSEO00O6hLdweXfAadEUVNTg+joaIwaNQqbNm3qsHBFnzxfz2BRYyjsYKP2Y9OiR0TTjI2NUVFRQcmil2MYBhUVFewYs7I4PXpavXo1Vq9ejczMTOzfvx9jxoyBQCBQOANAF0naeLSYiIto6l0d0FIlHrSvdXYV75TXedtwQrRNIBCgqKgIZWVl2g6FaJmxsTEEAoFKn+lSrydbW1vY2dnB2toaDx8+7MqheiQvB0vkVUpvMxA1oMXQVDsBEcKRkZERnJyctB0G0VGcEsWWLVuwf/9+lJWVYcaMGdi2bRs8PT3VHZvOKG7sg1U3atDXQH5JfrgNH+NtDdR+3vZTZKmzLCFEEzglinv37uGrr76Cv7+/msPRPQGW9Qrfz29ggPJWjSQKQgjpDpwSxfr169Udh85q2yRQ3roUq7OaOR2XpsgSQnoKlWY9hYeHAxC3G7ewsGC/JK/1UYOV+JGabU0+nO9fBgAUNeQis/K8NsMihJBuo1KiSE1NBSBuN15TU8N+SV7rq4KBvgCAoWWZ8DYZBgC4VZWhzZAIIaTbcC64k9fsSl/dtRPioYUjACDQzAMC02HaDYgQQroRpzGKmpoaREVFwcrKCrNmzcKMGTMwaNAgdcfWI3g5WAJ8MxgbGaCxuYXTMfIbGHasQlMzoAghRFM43VGsXr0aN2/exKZNm1BSUoIxY8YoXEtWFwnthErtV+vpilpP1w7fD7fhw9FUvHxsfgOD1PJWtcQnD3WVJYRoAhXccWAgErctNq3MYrcZ1xYDcmY9jbc1YO8gVJ0BlVXytOqbZkARQrSF0x3Fli1bEBERgXHjxqG8vBzbtm1DZmamumMjhBDSA6h8R8EwDK5cuUIFd4QQ0kuofEfB4/Fw9epVShIawC8thknSBhheOqftUAghhMXp0VNYWBguX76s7lh6tWb/ILTaO4BfWkzrUhBCehROg9mnT5/G1q1b4ejoCDMzMzAMAx6Pp/fjFMZGBnCyMQPKuR+j7VRZKf1GAOEjsPjXzbBsZiC/vSAhhHQ/Toni2LFj6o6jVwi34QOdTI9tbGEA8ChREEJ6DE6J4rvvvpO7/cMPP1TpOCkpKXjzzTfR0tKChQsXYtWqVTL7nDlzBm+99Raam5thY2OD3377jUvI3aKg7s8O3xvaz0VqqmxH7v/Kk7u97VRZCZoySwjpDpwShZmZGft9Y2MjDh8+DA8PD5WO0dLSgqVLl+L48eMQCAQQCoWIjY2VWteiqqoKS5YsQUpKCoYOHdrttRqSorvLpZ0/Uqt3dYDBfU1HRAgh3Y9TonjnnXekXi9fvhyxsbEqHePSpUtwcXGBs7MzAGDWrFlITk6WShR79uzB9OnTMXToUADiAj9tsqjIxcD8C1qNgRBCuhunWU/tNTQ04O7duyp9pri4GEOGDGFfCwQCmUaDd+7cwaNHjxAREYGgoCDs2rVLHeFyUikIBABYFVHXWEJI78LpjsLHxwc8nvhZektLC8rKylQen2AYRmab5JgSIpEI6enpOHnyJB4/foywsDCEhoZi+PDhUvslJSUhKSkJADS2eHyZY5hUkrjfXIFd5YcBAN4mwxBoptqjN0UcKovRJ2mD3Pea/YMgGjFSbecihJDOcEoUhw8ffnoAQ0MMGjQIhoaqHUogEKCwsJB9XVRUhMGDB8vsY2NjAzMzM5iZmWH06NG4fv26TKJISEhAQkICACA4OFjVy1GZZE0KQJwwAKgtUVx1DAAAOMl5j19aDCOAEgUhpFtxShSXLl3ChAkTYG5ujrVr1yIjIwN//etfERgYqPQxhEIhcnJykJeXBwcHB+zbtw979uyR2mfKlCl47bXXIBKJ8OTJE6SlpeHtt9/mErJaBZp5sIlBclcBAMa1+Wg0d+zSsS8OD8PF4WH4yNNI5j2TDu4yCCFEkzglir///e+YOXMmUlNT8csvv2D58uVYvHgx0tLSlD+xoSE2btyI6OhotLS0YMGCBfDy8sLWrVsBAImJifDw8MCECRPg6+sLPp+PhQsXwtvbm0vIeknelFkAEOVWqHScsGHW6giHEKKnOCUKAwNxLcCRI0ewePFiTJkyBWvWrFH5ODExMYiJiZHalpiYKPV6xYoVWLFiBZcw1cfeF1BiiqwikhqLoXJakaubYab0gLvIV/k7PUIIaY/TrCcHBwcsWrQIP/zwA2JiYtDU1ITWVs0tyEMIIUR7OCWKH374AdHR0UhJSUH//v1RWVmJf/zjH+qOTefUuzpoOwRCCFE7To+eTE1NMX36dPa1vb097O3t1RZUT2ZRkQvn+5dxV8mlUgkhRNeppeCut5AU3Q2vugknGzNxJ1kNyKplcPxhi0aOTQghqqJEoYIyxzDUWA+T2X7vSSky6m+p5RzhNuL/JKkddJmVLG5ECxwRQrqLSoli7ty5AICvv/5aI8HoIknx3R+Pc9VyvPG2BvA0l99BVrK4EQC1LnB0IbdC4RchpHdTaYwiPT0d9+7dw7fffot58+bJtOGwsrJSa3C6INDMQ21JojOiESPZqmwqviOEdBeVEkViYiImTJiAu3fvIigoSCpR8Hg8lRsD9kYFdX92Sy1FW+3rKtqiGgtCSGdUevT0xhtv4NatW1iwYAHu3r2LvLw89ouSBCGE6CdO02O3bNmC69ev4+zZswCA0aNHw9fXV62BEUII6Rk4JYoNGzYgKSmJraWYPXs2EhIS8Prrr6s1uJ5CaCcEHjfictUdmfecbMyA8qevW+ya0WDlKbNfW6L+ATLbFD0eIoQQbeKUKL755hukpaWxS6KuXLkSYWFhepso2jOtLoZb6ib2dX+TZgCAs0h9hXiSWorO1tgmhBBN45QoGIZhGwMC4iaB8hYi0keSorv2rhsBmTXX0U8NiSLcho+s2haklrcqTBSSmgoJTS1q1H6KLHWbJaR34ZQo5s+fj5CQEEybNg0AcOjQIbzyyitqDaynKnMMQ5ljmNQ2wb0ruF6dilN9Gai2crh8420NOiy4k2j2D0LbFStoUSNCiKZwShTLli1DREQEUlNTwTAMduzYgYAA2efuesVpFHBVdowCENdSFJWndms4bWsqAKqrIIRoDqdEAQCBgYEqrWjXm5hWZnU6oE0IIbqCc6Ig3N2uugoAcJMz+6m7yZttRUV4hJC2KFGogZONGfgFPPZ7AIC9Jfv+zWL5S5YSQogu4NQ9duPGjXj06JG6Y+nxhP2HQ9h/uEaOLe+v+PwGBquzmrE6q5najhNCtIZTorh//z6EQiFeeOEFpKSk9Jqpsd0p3IYPR1PxXUp+A9PpLChCCNEUTo+e1q5di7///e/49ddfsWPHDrz22mt44YUX8Morr2DYMNn1Gnql0kzAnntbk/G2BmwNxeqsZqU+I6mr0FQ9hQSX1uNUe0GI7uK8cBGPx4OdnR3s7OxgaGiIR48eYcaMGXj33XfVGR9RkmStCnWuU0EIIQDHRLFhwwYEBQXh3XffxciRI3Hjxg1s2bIF6enp+Omnn9QdI1GCaMRIPE54g13YiBBC1IXTo6fy8nL897//xTPPPCO1nc/n4/Dhw2oJjBBCSM/A6Y6iqalJJkmsXLkSAODh4dH1qPRFaSZQmgkv/j32ixBCdA2nRHH8+HGZbceOHetyMLrOoPkxBuZf0Mix206V1fSUWcPMDJkvQkjvpVKi2LJlC3x8fHD79m34+vqyX05OTr1q4SJ5tRTNffsBAKyK1P9Lte1UWQmaMksI6S4qjVG8/PLLmDhxIt577z2sX7+e3W5ubg4rKyu1B6dLmo0tYNRUByg3k1UlbafKSig7ZZYQQrpKpURhaWkJS0tL7N27V1Px9GxOo4C8sx2+nWPA4DVTERpyf5R5L6S/OyKsfTQZXY+mTO0F1VoQ0jOp9OgpPDwcgPgOwsLCgv2SvO7NQvq7w7WFJ/e9gsdlSKvKBiDuLGtamdWdoRFCSJeodEeRmipec6G2tlYjweiyCGsfLLr1OwDgduBMqfc+a3OH0bZpoNBO+i/o+nJxI8GsEmoiSAjpOThXZhNCCOkdVLqjMDc3B4/Hk9sEkMfjoaamRm2Bkc5Jpsy2tbhB/N9miwqD3eE2fIVrcxNCejeVEgU9cuo5wm34gBqmx+Y3MEB5KyUKQkiHVEoU4eHhSE1NZe8s2qM7ChWUZuJyaSaEAa+wm8xCRoi/OShb0NievCmzAGCSKv7v8pGnkVJhKDvNtqOiO1oNjxD9p9IYRdvB7JqaGpkvVaWkpMDNzQ0uLi5SdRntXb58GQYGBjhw4IDK5yCEENI1WlsKtaWlBUuXLsXx48chEAggFAoRGxsLT09Pmf1WrlyJ6OhoLUWqHgWPy/BZ7o/dUk8hWZdCQtPrU6gLrXNBSM/EadZTY2MjvvjiC0yfPh1xcXH48ssv0djYqNIxLl26BBcXFzg7O6NPnz6YNWsWkpOTZfb797//jbi4ONja2nIJVf2cRqn8kZD+7hhqMlCqnkJTJOtSSND6FISQruJ0RzFv3jyYm5vj9ddfBwDs3bsXc+fOxY8/ylYkd6S4uBhDhgxhXwsEAqSlpcnsc/DgQZw6dQqXL1/mEmqPEGHtgwhrH6l6Ck0RjRgpdffQ9s6CEEK44JQobt++jevXr7OvIyMj4efnp9IxOppi29Zbb72Fzz77DAYGimfkJCUlISkpCQBQVlamUhzqZlGRi4H5F1DmGKbcByQtQTjcqRBCSHfglCgCAgJw8eJFhIaGAgDS0tIwcqRqz8AFAgEKCwvZ10VFRRg8eLDUPleuXMGsWbMAiBdLOnr0KAwNDTF16lSp/RISEpCQkAAACA4OVvVy1KZSEAiLilxYFWUonyh6AHn1GG1RnQUhvZtKicLHxwc8Hg/Nzc3YtWsXhg4dCgAoKCiQGYTujFAoRE5ODvLy8uDg4IB9+/Zhz549Uvvk5eWx38fHx2Py5MkySaInKXMM00ibcU3qrB6D6iwIISolCnUuc2poaIiNGzciOjoaLS0tWLBgAby8vLB161YAQGJiotrOpQlt16S4XHVHrcdWVJug7kWEOqrHkKB25oQQlRJF2+VPHz16hJycHKnZTu2XR+1MTEwMYmJipLZ1lCB27typ0rEJIYSoB6cxim+++QZff/01ioqK4O/vj4sXLyIsLAynTp1Sd3w6SZUBbfZu5L4xu+12lbh7rFv/AI3Ep05t73C0UaWtbO0F1VsQwh2nRPH111/j8uXLCA0NxenTp5GdnY3Vq1erO7aeq/0MJckv+dLMTge0b9cX40zFDYVFd14OlrhZrL5W4+0L8NrSlWI8Qoj2cCq4MzY2hrGx+JdjU1MT3N3dcfv2bbUGpqvKHMNQYz1M7nsh/d0BQONFd221L8Bri4rxCCHK4HRHIRAIUFVVhalTp2L8+PEYMGCAzNRWIivC2qdbkwQgW4DXlrLFeFm10tNnabosIb0Lp0Rx8OBBAMCaNWsQGRmJ6upqTJgwQa2B6TqVC+96qHAbPoCn02dpuiwhvQ+nRNHY2IjNmzcjNTUVPB4P4eHhaG3t+toI+kJXC+/kaT99lqbLEtL7aK3Xkz4R2glx+f7TXlScCu9KM6VemlbWA3JmPbWdWaTumgpCCJFHa72eCCGE6AZOs54kvZ4kuPR6Ioo52ZhpOwRCCAGgxV5PvYFpdTHcUjdJbTMaYIxmYwstRdR7KSrMo2I8QhTTWq8nfVcpkK1SNq0uhpHpIOQyTd222l1n5BXjUREeIaQtzr2erl+/jrNnxWspjBo1isYo2ilzDJOZ8eSWugnjm1vRYCle7Q6AwkRxu+qqzDZ1tvVo9g+CUbtt/NJiGAEKE0WHbcmz0mQ2Mf36AQBGDrPBOI9BXYiWEKItnFt4bNu2DdOnTwcAzJkzBwkJCewsKNKxKc18uA+bqdRqd6aVWQCABivNPNaTV4zXWRFeZ23J5blX0QCgnBIFITqKU6LYvn070tLSYGYmHnBduXIlwsLCenWiENoJcbndFFdNa9+Erzumy3bWlrw9ka8XPj58U4MREUI0jdOsJ4ZhpJYnNTAwkLu0aa9j79vpLpIBbqPGmk73dbIxg5ONGbwcLNURHSGEcMLpjmL+/PkICQnBtGnTAACHDh3CK6+8otbA9JFkgFsyqE2znwghuoBToli2bBkiIiKQmpoKhmGwY8cOBAT0/LUTtE0ywC2eMivSdjiEEKIUlRMFwzAoKipCYGAgAgO7f6EaQggh3UvlMQoej4epU6dqIJTehS9qgml1CQbmX+h859JMmFZmwbQyCxb3L8Li/sXOP0MIIWrC6dFTaGgoLl++DKFQqO54eoVKQSBaKy+BL2rqkR1mNVGEd6+iodPZT9qqtVB2OVVFqLqb6DNOieL06dPYunUrHB0dYWZmBoZhwOPxkJnZvdNDdVWZYxgaWopgWl0CPFHfcSXTZbsyTZZrEZ4iI4fZAChXuA/VWhDSc3FKFMeOHVN3HKQTkiaBDbgPQPwXrDr+Em6PSxGeIoaZGYgGED1U8X6r61qBujqFSa593QghpHtwShSDBg2SWbho8eLF6o6NEEJID8Cp4G7evHm4efMmXn/9dbz22mu4desW5s6dq+7YegWLily4pW5SblCbEEK0gBYu0qLmvv1QYz0MFhW5AKDaoHbeWVjcF1d319iFaiI8QggBQAsXaVWzsQVuhy9FjfUwbYdCCCEd4nRHkZaWJrNwkYeHB7uwUa+e/WTvK7P+NSGE6DJOiSIlJUXdcegFod3/6koeN8q8d7nqjtrP5zX4f72inKTn8NeXq9ZEMKukWl0hdUmH61z8D1MgW4dB61wQonmcEkXbBYwId7fri3Gm4gbctB2IEuQV4QHqWw1P19e50MRU5a6gAkCiTpwSBem6kP7uuF1fjLSqbCzSdjCdkFeEB3S9EK8tZda5EPl6Sb2mdS4I6R6UKDTBaRSQd1bhLhHWPkiryu6mgLpGXhEe0LVCPEKI7uA060me+/fvq+tQRFWdJCVCCOkKtSUKWriIEEL0k9oePR05ckRdh+qVJBXabVUKAjsswpOZRXXfmP3WU8Vzew5WPEuqp8yKat8HilfXLHe7BPWGIkQ9ON1RrFy5UqltvZrTKKmXwv7DIew/XO6ulYJAmaI70+piWBVx6wJrFjJC6osQQrqC0x3F8ePH8dlnn0ltO3bsmMw2ohzJEqlttb+7IPIpqr1oW3dB9RaEcKfSHcWWLVvg4+OD7Oxs+Pr6sl9OTk7w8fFR+eQpKSlwc3ODi4sL1q9fL/P+7t272XM8++yzUv2lCAm34cPRlNfpfvcqGnAuV/F6GISQjql0R/Hyyy9j4sSJeO+996R+sZubm8PKykqlE7e0tGDp0qU4fvw4BAIBhEIhYmNj4en59Am7k5MTfvvtNwwYMADHjh1DQkIC0tLSVDpPT1fwuAyf5f4os93UVITxzXy4ayEmVbQtxFNX8Z2yOqu9kNRdUL0FIV2jUqKwtLSEpaUlpk+fDisrK5ibm2Pt2rXIyMjA3/72NwQEBCh9rEuXLsHFxQXOzs4AgFmzZiE5OVkqUTz77LPs96GhoSgqKlIlXO1rN04BALj6dBA6pH/HaSDHgAFf9BhTUjcpHNTWpraFeOosviOE9Cycxij+/ve/Y+bMmUhNTcUvv/yC5cuXIzExUaW/9ouLizFkyBD2tUAgUPj57du3Y+LEiXLfS0pKQlJSEgCgrKxM6Ri0LcLaBxHW8h/ZfXFzO1pb6mBaVgxAiRbkbRsRtu81df92V8KEZ7sHlPnl9eJvXAYALs+Jv9+XDL6oHqaVWR0ep8FK1flYhKue1lKEaJamW7ZwShQGBuLb/SNHjmDx4sWYMmUK1qxZo9IxGIaR2cbjyX/efPr0aWzfvh2pqaly309ISEBCQgIAIDg4WKU4upuw/3ClGgQ2G1ug2dgCDQ1dn8Fs5ttxN6n6zK4lEUKI/uM0PdbBwQGLFi3C/v37ERMTg6amJrS2qtbQTSAQoLCwkH1dVFSEwYMHy+yXmZmJhQsXIjk5GdbW1OiMEEK6G6dE8cMPPyA6Ohq//PIL+vfvj8rKSvzjH/9Q6RhCoRA5OTnIy8vDkydPsG/fPsTGxkrtU1BQgOnTp+P//u//MHy4/BoEneM0SrxmBSGE6AhOzzVMTExQX1+PvXv34sMPP0RzczP69++v2okNDbFx40ZER0ejpaUFCxYsgJeXF7Zu3QoASExMxMcff4yKigosWbKE/cyVK1e4hEx6uXsVDR3OfqIaC0IU45QolixZAj6fj1OnTuHDDz+Eubk54uLicPnyZZWOExMTg5iYGKltiYmJ7PfffPMNvvnmGy4h6hXT6mKV2nsQaSOH2QCQX0fRk9a0IKSn4rwUakZGBjsddsCAAXjy5IlaAyNilQLZfkWm1YpnQikaLO+ojYg+kvSAigYQPVT+PqvrWoG6ug77RfU01L+KaAOnRGFkZISWlhZ2llJZWRn4fLU1ou0dJOMUCtbXLnhchuUmAOxt2W0h/d2x6NbvGg5OMUcbM5ltD4wM8KSkDH1+Oiy13czfA/1CfJ9OqSWE6BxOieKNN97AtGnT8PDhQ3zwwQc4cOAA1q5dq+7YejV5xXgFj8U1Ij1xRTwzfw+ZbU9KxPH2CxEnRUU1FppCtRuEdB2nRDF79mwEBQXh5MmTYBgGhw4dgoeH7C8KIp/QTggAuHy/4zEdecV48lp9dJW8GgsutRX9QnzZhCDx4D/7OcdFCOk5OFdzubu7w929p3ci0gFKPIIihBBtojWzdVTbmVA0A4qQ3k3SskVTrTwoUeigtjOhOpsBRTqnaE0LZYXb8BV2siVEl1Gi0EFtFzrSlQWO5M2Uak8bM6PCbfhAuWrtZ9rLb2CA8lZKFERvUaIgvVpna1ooo6t3I4T0dJQoiAx1zYQihOgHShQ6RrIiXkh/9w7XsuhJnpSUKT1NltfcovB9xsMV8KO6CEK6GyUKHSIpwpMU3vX0RCGvCI+zhxXgAWBUTBRti/yo+I4QbihR9BTyWo+3q62QFOG1L7yzqMjFwPwLPW7mk7wiPEUUDWbz9iWrIyRCCAeUKHRcpSAQFhW5sCrK6HGJQlWKZkY9MBIPOA9qtw/1kCJE8yhR6LgyxzBYFelG51N9po5aDGUwBfLX1GiP1tgg6kQtXwnponAbPhxN5a/3rg33KhpwLlf++huEcEF3FFoktBMqbAzYk8ibMgvQtFlAPbUYyhL5enW6T0cr+RHCFd1R9GT2vrS+NiFE6+iOQk/IWy5VLkMTtZ7XuP6x0vuKBgRAZKPbA+6E9EaUKPSAvOVSexr+42IYAmpPFMr0kAJodhQhXUGJQkdJKrRZ7ZZL7agYT91rZjcqOUZhnLO5y+dSpcobeLoMKyGkayhR6IJ24xQhokrg3km5u+pK1baqVK3ybr8MKyGEO0oUWiZZFlWejmZERQyJQIShldz3NLFcqiJtZ0NpcgaUqlXevX0Z1nsVDTT7qRfRdN0MJQpC9MzIYTYAqI6it7hX0QCgnBIFIV3laGMG2FnKfS+rpLqbo9GscR6DqCq7F+mOO0dKFLqq7bhFu+aBPRX/cbFaBrWVOldDEwDAsPwCTcklpIuo4I50C9GAALSaOHTvSVubYPjoaveekxA9RHcUPZi8ge6e3PJD0cp4Ipuwbv3LvtV0P/gNJd12PkL0GSUK0ut5DpY/dqEKfRvnIKQtShR6SKYYr40XB4/BTPvwbo5IO5oqW1F4uAmtpv+bKtvniMbOZdIk0tix22qaNA1PJk7plnMRIkFjFPqgzcB2SH93DDUZKHe3gsdlOPqw5z66Uiczfw/0tdKvf9780mL0+e24tsMgPZCkbmZPWoFGjk93FHpGslyqPN1djKdN/UJ8YWOVCv7jYrSa9AUAiJ55FqKhkzRyvu549GSStEHj5yC6R1I3c6+iAcnXivFyyFC1n4MShY6RDHD35EHttjpax4ILVSu/RQMC2H/g/MfFMCw9pbFEQYi2SOpmNFlPQYlCXyizbkXhMaCPOeA0SvPxyJN3tltP13amlabrN9QxIN6ZB30NYWBsBJ9h1p3ueyG3QuPxkN5Dvx7iEkIIUTutJoqUlBS4ubnBxcUF69evl3mfYRi88cYbcHFxga+vLzIyMrQQJSGE9G5ae/TU0tKCpUuX4vjx4xAIBBAKhYiNjYWnpye7z7Fjx5CTk4OcnBykpaVh8eLFSEtL01bIPQqXrrNa18VHXmYqfr4+7RJwXzfam6hbmBKPp4j+0PSjRq3dUVy6dAkuLi5wdnZGnz59MGvWLCQnJ0vtk5ycjHnz5oHH4yE0NBRVVVUoLS3VUsSEENI7ae2Oori4GEOGDGFfCwQCmbsFefsUFxfD3t6+2+LURR3dbZj3Me/mSLTLLGQEgBHiFzsOttmmmwzMe9d/P6K8sGHWsDA20tjxtZYoGIaR2cbj8VTeBwCSkpKQlJQEAMjOzkZwcDCnmMrKyjBwoPxiNX0S/NenP5/ecs2sTcG6f80c/n3r/DVz0FuvOXgbt2vOz8/v8D2tJQqBQIDCwkL2dVFREQYPHqzyPgCQkJCAhISELscUHByMK1eudPk4uoSuuXega+4dNHXNWhujEAqFyMnJQV5eHp48eYJ9+/YhNjZWap/Y2Fjs2rULDMPg4sWLsLS0pMdOhBDSzbR2R2FoaIiNGzciOjoaLS0tWLBgAby8vLB161YAQGJiImJiYnD06FG4uLjA1NQUO3bs0Fa4hBDSa2m1MjsmJgYxMTFS2xITE9nveTweNm3a1G3xqOPxla6ha+4d6Jp7B01dM4+RN2JMCCGE/A+18CCEEKIQJYr/6aydiK4qLCxEZGQkPDw84OXlha+//hoAUFlZifHjx8PV1RXjx4/Ho0eP2M+sW7cOLi4ucHNzwy+//KKt0LukpaUFAQEBmDx5MgD9v14AqKqqwowZM+Du7g4PDw9cuHBBr6/7yy+/hJeXF7y9vfHSSy+hsbFRL693wYIFsLW1hbe3N7uNy3Wmp6fDx8cHLi4ueOONN+SWH3SIIYxIJGKcnZ2Z3NxcpqmpifH19WVu3ryp7bDUoqSkhElPT2cYhmFqamoYV1dX5ubNm8yKFSuYdevWMQzDMOvWrWPeffddhmEY5ubNm4yvry/T2NjI3L17l3F2dmZEIpHW4ufqX//6F/PSSy8xkyZNYhiG0fvrZRiGmTdvHrNt2zaGYRimqamJefTokd5ed1FREePo6Mg0NDQwDMMwM2fOZHbs2KGX1/vbb78x6enpjJeXF7uNy3UKhULm/PnzTGtrKzNhwgTm6NGjSsdAiYJhmPPnzzNRUVHs608//ZT59NNPtRiR5sTGxjK//vorM3z4cKakpIRhGHEyGT58OMMwstceFRXFnD9/XiuxclVYWMiMHTuWOXnyJJso9Pl6GYZhqqurGUdHR6a1tVVqu75ed1FRESMQCJiKigqmubmZmTRpEvPLL7/o7fXm5eVJJQpVr7OkpIRxc3Njt+/Zs4dJSEhQ+vz06AkdtwrRN/n5+bh69SpCQkLw4MEDtibF3t4eDx8+BKAfP4u33noLn3/+Ofj8p/+89fl6AeDu3bsYOHAg5s+fj4CAACxcuBD19fV6e90ODg5Yvnw5hg4dCnt7e1haWiIqKkpvr7c9Va+zuLgYAoFAZruyKFFA+VYhuqyurg5xcXH46quvYGFh0eF+uv6zOHz4MGxtbREUFKTU/rp+vRIikQgZGRlYvHgxrl69CjMzM4Vjbbp+3Y8ePUJycjLy8vJQUlKC+vp6fP/99x3ur+vXq6yOrrOr10+JAsq3CtFVzc3NiIuLw+zZszF9+nQAwKBBg9hOvKWlpbC1tQWg+z+Lc+fO4eeff4ajoyNmzZqFU6dOYc6cOXp7vRICgQACgQAhISEAgBkzZiAjI0Nvr/vEiRNwcnLCwIEDYWRkhOnTp+P8+fN6e73tqXqdAoEARUVFMtuVRYkCyrUT0VUMw+CVV16Bh4cHli1bxm6PjY3Fd999BwD47rvvMGXKFHb7vn370NTUhLy8POTk5GDECN3puLpu3ToUFRUhPz8f+/btw9ixY/H999/r7fVK2NnZYciQIbh9W7yu+MmTJ+Hp6am31z106FBcvHgRDQ0NYBgGJ0+ehIeHh95eb3uqXqe9vT3Mzc1x8eJFMAyDXbt2sZ9RCrehFf1z5MgRxtXVlXF2dmbWrl2r7XDU5uzZswwAxsfHh/Hz82P8/PyYI0eOMOXl5czYsWMZFxcXZuzYsUxFRQX7mbVr1zLOzs7M8OHDVZoZ0dOcPn2aHczuDdd79epVJigoiPHx8WGmTJnCVFZW6vV1f/jhh4ybmxvj5eXFzJkzh2lsbNTL6501axZjZ2fHGBoaMg4ODsw333zD6TovX77MeHl5Mc7OzszSpUtlJj4oQpXZhBBCFKJHT4QQQhSiREEIIUQhShSEEEIUokRBCCFEIUoUhBBCFKJEQQghRCFKFIQQQhSiREGIiqqqqrB582b29bPPPquR8xQVFWH//v0aOTYhqqBEQYiK2ieK8+fPa+Q8J0+eREZGhkaOTYgqKFEQoqJVq1YhNzcX/v7+WLFiBfr16wdA3Mbd3d0dCxcuhLe3N2bPno0TJ05g5MiRcHV1xaVLl9hjfP/99xgxYgT8/f2xaNEitLS0SJ0jNTUVy5Ytw4EDB+Dv74+8vLxuvUZCpKixJQkhvUL7RWTMzMzY7QYGBkxmZibT0tLCBAYGMvPnz2daW1uZQ4cOMVOmTGEYhmGysrKYyZMnM0+ePGEYhmEWL17MfPfddzLniY6OZm7cuKH5CyKkE4baTlSE6BMnJyf4+PgAALy8vDBu3DjweDz4+PggPz8fgPiRUnp6OoRCIQDg8ePHbJvotm7fvg03N7dui52QjlCiIESN+vbty37P5/PZ13w+HyKRCIC49ftf/vIXrFu3rsPjVFRUwNLSEkZGRpoNmBAl0BgFISoyNzdHbW0t58+PGzcOBw4cYJevrKysxL1796T2ycvL0+mFdYh+oURBiIqsra0xcuRIeHt7Y8WKFSp/3tPTE2vXrkVUVBR8fX0xfvx4drUyCXd3d5SXl8Pb21tjs6oIURatR0EIIUQhuqMghBCiECUKQgghClGiIIQQohAlCkIIIQpRoiCEEKIQJQpCCCEKUaIghBCiECUKQgghCv1/2wGnXkBmWC8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for value in data_x[\"Celltype\"].unique():\n",
    "    mask = data_x[\"Celltype\"] == value\n",
    "    time_cell, survival_prob_cell, conf_int = kaplan_meier_estimator(\n",
    "        data_y[\"Status\"][mask], data_y[\"Survival_in_days\"][mask], conf_type=\"log-log\"\n",
    "    )\n",
    "    plt.step(time_cell, survival_prob_cell, where=\"post\", label=f\"{value} (n = {mask.sum()})\")\n",
    "    plt.fill_between(time_cell, conf_int[0], conf_int[1], alpha=0.25, step=\"post\")\n",
    "\n",
    "plt.ylim(0, 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": [
    "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 structured array of survival times and binary event indicators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CoxPHSurvivalAnalysis()"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn import set_config\n",
    "from sksurv.linear_model import CoxPHSurvivalAnalysis\n",
    "\n",
    "set_config(display=\"text\")  # displays text representation of estimators\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",
    "    {\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",
    "    },\n",
    "    columns=data_x_numeric.columns,\n",
    "    orient=\"index\",\n",
    ")\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 0x7fd9aec07280>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEICAYAAABBBrPDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABJ+klEQVR4nO3dd3hUZfbA8e9JMpl0kpBAgAChhE4ogoKgolRFBBvY1rZ2UdfddXV/lsXCuvayVuyoa1krskhRUZRepPdO6ElISK/v7487CamTzGTSJufzPPfJzL133nvuiDm5977vecUYg1JKKVUVn4YOQCmlVOOmiUIppZRTmiiUUko5pYlCKaWUU5oolFJKOeXX0AF4WlRUlImLi2voMJRSqklZvXp1kjEmurJtXpco4uLiWLVqVUOHoZRSTYqI7Ktqm956Ukop5ZQmCqWUUk5polBKKeWU1z2jUEp5r/z8fBITE8nJyWnoUJqsgIAAYmNjsdlsNf5MgyYKEXkXuBA4ZozpU8l2AV4CLgCygOuNMWvqN0qlVGORmJhIaGgocXFxWL8elCuMMSQnJ5OYmEinTp1q/LmGvvX0PjDOyfbzgXjHcgvwej3EpJRqpHJycmjZsqUmCTeJCC1btnT5iqxBE4UxZhGQ4mSXicBMY1kGhItIm/qJTinVGGmSqB13vr+GvqKoTjvgQKn3iY51ZYjILSKySkRWHT9+3K0D5Wenc/jrh8jc/rNbn1dKKW/V2BNFZamvwgQaxpgZxphBxphB0dGVDiysVlJqEit2vUvK55fCnkVutaGU8n7Tp0+nd+/eJCQk0L9/f5YvX16nxxsxYoRLg4j/+9//0rt3b3x8fDw2+LixJ4pEoH2p97HAoTo5kggZP4Xz6olW8MEE8hPX1clhlFJN19KlS5k9ezZr1qxh/fr1/PDDD7Rv3776D9ajPn368NVXX3H22Wd7rM3GnihmAdeKZQiQZow5XBcHamkPY+Buw8U/QZYIx1d8XheHUUo1YYcPHyYqKgq73Q5AVFQUbdu2BeCxxx5j8ODB9OnTh1tuuYXi2UNHjBjBvffey9lnn03Pnj1ZuXIll1xyCfHx8Tz00EMA7N27lx49enDdddeRkJDAZZddRlZWVoXjz58/n6FDhzJw4EAuv/xyMjIyKuzTs2dPunfv7tHzbujusZ8AI4AoEUkE/gHYAIwxbwBzsLrG7sTqHntDXcXiHxFJdlQoYUnpPBMewUUpybStq4MppWrt0e82sfnQSY+22attGP+Y0LvK7WPGjOGxxx6jW7dujBo1iilTpnDOOecAMHXqVB555BEA/vCHPzB79mwmTJgAgL+/P4sWLeKll15i4sSJrF69msjISLp06cK9994LwLZt23jnnXcYNmwYN954I6+99hp//etfS46dlJTEE088wQ8//EBwcDBPPfUUzz//fMkx61JD93q60hjTxhhjM8bEGmPeMca84UgSOHo73WmM6WKM6WuMqdNqfxETJwGw6WQwvse/hCxnHbKUUs1NSEgIq1evZsaMGURHRzNlyhTef/99ABYuXMgZZ5xB3759+emnn9i0aVPJ5y666CIA+vbtS+/evWnTpg12u53OnTtz4IDVX6d9+/YMGzYMgGuuuYbffvutzLGXLVvG5s2bGTZsGP379+eDDz5g374q6/h5lI7MLqXVWSPZ/86H3PdlIbmX+MD//gyXv9/QYSmlKuHsL/+65Ovry4gRIxgxYgR9+/blgw8+4IorruCOO+5g1apVtG/fnmnTppUZq1B8q8rHx6fkdfH7goICoGK31fLvjTGMHj2aTz75pK5OrUqN/RlFvQrs2wcTFEBIDhSc8INNXzd0SEqpRmTbtm3s2LGj5P3atWvp2LFjSVKIiooiIyODL774wuW29+/fz9KlSwH45JNPGD58eJntQ4YMYfHixezcuROArKwstm/f7u6puEQTRSk+wcGYfz0AwM/+kdbKXQsbMCKlVGOSkZHBddddR69evUhISGDz5s1MmzaN8PBwbr75Zvr27cukSZMYPHiwy2337NmTDz74gISEBFJSUrj99tvLbI+Ojub999/nyiuvJCEhgSFDhrB169YK7Xz99dfExsaydOlSxo8fz9ixY90+32JS/GTeWwwaNMjUpu9wxrJlHLj+Bl6YEMeM4CUQGAH37/VcgEopt23ZsoWePXs2dBget3fvXi688EI2btxYL8er7HsUkdXGmEGV7a9XFOWIY4xfVLqdYyYcsk/AyncaNiillGpAmijK8e/YAYCOWblMyXvYWrnxywaMSCnl7eLi4urtasIdmijKsbVpQ5Y/nDT72OMTZq3ct7hhg1JKqQakiaISfj5Wr+FWLfL5sXCAtTJ5VwNGpJRSDUcTRSX8xBeA8/rl8m3hmdbKDa53d1NKKW+giaISPuKDAEksYU3gUGvluvof5KKUUo2BJorK5OYxfqWBtHQSMx1f0Yk9sGV2w8allGpwjb3M+H333UePHj1ISEjg4osvJjU1tdYxaKKoRNg4a3bWIwe2cvlpsVySO83asEq7ySrVnDWFMuOjR49m48aNrF+/nm7duvHkk0/Wuk1NFJUIHTUSgAIKaddS+N10tTYU5jdgVEqphtYUyoyPGTMGPz+rQ86QIUNITEys9XlrUcDK+FoPswftMPSfGITBhwMhCbTf+yukH4HQmAYOUCnF9w/AkQ2ebTOmL5z/ryo3N7Uy4++++y5Tpkyp9deiVxSVCB4yBIB+ewz/O/gmAOtzHcnhOc9OCKKUajqaUpnx6dOn4+fnx9VXX13r89Yrikr4hIQA4J9vyC3MBeCu9GsZHzDf2uHQ79B2QEOFp5QCp3/516WmUGb8gw8+YPbs2fz4448V2nGHXlFUQnx8CD7zTOy2QBYeWMirVw2kCB/m9H3J2mHnDw0boFKqQTSFMuNz587lqaeeYtasWQQFBbkcR2X0isKJjvuyEeNLcIvdALy3N5ILABLrdKI9pVQjlZGRwV133UVqaip+fn507dqVGTNmlCkzHhcXV6sy47feeivx8fFOy4zn5lp3Op544gm6detWZr+pU6eSm5vL6NGjASvBvPHGG26esUXLjFdh/823kPnrr/zxHl/G9LuML+cPpW+7Fvzn0PlgimDKR9BzggciVkrVlJYZ9wwtM+4hIY6eDC0DogDo3jqUJbuSybngRWuHz66BoqIGik4ppepPrRKFiASLOAojealhKzM4mXeSYLt1l+7uLb2hg6P+0z/bwKe171GglGrevKrMuIj4iMhVIvI/ETkGbAUOi8gmEXlGROLrJsz6FzLc6qYWm5jLgn0LePlKq5fT/M1HybrwdTjteijIga2z4Ujj/Q+slFK15eoVxUKgC/B3IMYY094Y0wo4C1gG/EtErvFwjA3CPy4Oe48etLCHEW4Pp0WgjbvOs0Zorz0ZAhNegskfWjvPuqsBI1VKqbrlaqIYZYx5HEgzxpTcoDfGpBhjvjTGXAp85tEIG1JREV03pJCam0pWfhZnxUcD8PuBVGt78cPsQ2u0DLlSymu5lCiMMcXFjr4uv01EhpTbp8mTwAAAItINu1J30SkqGIDj6bmOHQQuesV6/eUfoaiwIcJUSqk65eoziski8i8gVER6lnuQPcOzoTW88EmTALh4SREGQ3SonSB/X95fspfsPEdSGPgHiB9rvf65YUaKKqXqT2MvM/7www+XxDZmzBgOHTpU6xhcvfW0GNgMRADPAztEZI2IzAayax1NIxN+xRUYP1/8imDBvgUAxLW0rioe/e5UHRcufN76mbStvkNUStWjplBm/L777mP9+vWsXbuWCy+8kMcee6zWbbp66+mgMWYmMNEYc74xpjMwCvgHcF6to2lkRAS/iEgAbD42AN653hqP8unKAxQUOh7TtIiFqO6QV7EssFLKezSFMuNhYWElrzMzMz1S68mlEh4iIsayuHidMSYFSCm/T60jaySKv+L1SesBaNMikLG9WzNv01HOfnohY3rHMO2i3pCfDTsXNFygSjUzT614iq0pWz3aZo/IHtx/+v1Vbm8qZcYffPBBZs6cSYsWLVi4cGGtvxeXu8eKyF0i0qH0ShHxF5HzROQD4LpaR9WImIICzltnPcwudv+4Hozp1ZpDaTm8v2QvqVl50LpXA0aplKoPTaXM+PTp0zlw4ABXX301r7zySq3P29WigOOAG4FPRKQTkAoEYiWc+cALxpi1tY6qEfEJCiI3L4Ok7CTm753PmLgxdI4OYca1g3h+wXZe/nEHn686wC2tesL2uQ0drlLNhrO//OtSUygzXuyqq65i/PjxPProoy6fZ2muPqPIMca8ZowZBnQERgIDjDEdjTE3u5okRGSciGwTkZ0i8kAl21uIyHciss4x+vsGV9r3hODhwwmwW6V6j2UdK7Pt2qEdAfjP8v1QkGet1FHaSnmtplBmvHR8s2bNokePHi7HUp7TRCEivUTko1LvfxSR3lAyXmIwMFVETnf1wI6uta8C5wO9gCtFpPz9mzuBzcaYfsAI4DkR8Xf1WLViDJKSRmCOqZAookKsvwz2Jmfxa4Fj5ruVb9VreEqp+pORkcF1111Hr169SEhIYPPmzUybNq1MmfFJkybVqsx4QkICKSkpTsuMJyQkMGTIELZurfiM5oEHHqBPnz4kJCQwf/58XnrpJbfPt5jTMuMichgYaozZ63i/zRjT3fH6TOB7rJHYw4EHjTEVBuI5aXsoMM0YM9bx/u8AxpgnS+3zd6A9VsKIAxYA3UqPCi/PU2XGix177jmS33qbv1/ni71vbz67sOzA81nrDnH3J7/jRwE7A66FwTfB+Oc8dnyl1ClaZtwzPF1mfAwwvdT7k6VeXwu8YYy5BeuvfVdvGLYDDpR6n+hYV9orQE/gELABuKeyJCEit4jIKhFZdfz4cRfDcC5okPW9tQ1uQ6gttML2i/q1ZXBcBAX4UWALhZVvQ1ZKhf2UUqqpcpoojDEbjDGl62jvFJHLRKQVMAn41rHfMcBeSRPOVNa5t/zlzVhgLdAW6A+8IiJh5fbBGDPDGDPIGDMoOjraxTBqpn1G1ad3x7lWscCtNkeGfnskFHpNJROlVB3zqjLjwL3ArcBBYI0xZgmAiNiAEBfbSsS6rVQsFuvKobQbgK8cYzd2AnuA2j+ZcYGtnXWR0zoxk9zC3Er3Obd7K+JbhfBR6/usFSm7YYXXVTRRSjVTrvZ6OmKMGQ3YjTEXlNp0LlYJclesBOJFpJPjAfUVwKxy++zH6lmFiLQGugO7XTxOrdi7dgU/PwrFsPb4WtYcXVPlvif9WsJtjrGI8/4Psk/UU5RKKVV33JrhrvxzAmPMfMezClfaKACmAvOALcDnxphNInKbiNzm2O1x4EwR2QD8CNxvjElyJ+ZaKSigb5B1e2nG+sqvFPILi5iz4QjE9IHujhw67yFY+hps+qaeAlVKKc9zdcCdRxlj5gBzyq17o9TrQ1gP1BtcxG+boA8sPrS40u3RoXb2Jjtqs0x4CXYthLUfndqh0x4IiqyHSJVSyrNqNWd2c+HfpQtFJ1K5NvYSAHae2Flhn0Fxkfj7Or7OkFZw/164fx+M/Ie1bv1nkHYQvKcMllLNUmMvM17s2WefRURISqr9TRhNFDUQdNppAAzdbM1BkZRT+RefV1jqjpwtAALDIWGK9X7uA/BCL1ihA/KUaqqaQplxgAMHDrBgwQI6dOhQ/c414OrERekicrKSJV1ETlbfQtPU6s9WdcdQn8Aq98nJt5LI2uJpUou1aAdXfnZqJrzMsqO7lVJNR1MoMw5w77338vTTT3ukxDi4+IzCGFNxxFlz4PiyfcTKq/P2zmNImyFldhnRvRXvLd7L7HWH6N8+vOznu4+zfs6aCouegXPuB19bXUetlFc78s9/krvFs2XG7T17EPN//1fl9qZQZnzWrFm0a9eOfv36eex7cfvWk4hEiMjpInJ28eKxqBobH+trCv5hJQB234qD7wbHRQCQW1BldREIcww836HzVijVFDX2MuNZWVlMnz7dI7PaleZWrycRuQm4B2uQ3FpgCLAUL5zlDsA31LqQ8g1rQZBfEJn5mRX2CfL3IyzAj/0pTma5m/QazJxojbHoNhZ8fKveVynllLO//OtSYy4zvmvXLvbs2VNyNZGYmMjAgQNZsWIFMTExbp+zu1cU92BVjt1njDkXGAB4tshSIxPoeKCdX5TPNzu/4XhWxdPNKSjil+3H2Z9cRbLoPAJiB8OJPTD/4TqMVilVFxp7mfG+ffty7Ngx9u7dy969e4mNjWXNmjW1ShLgfqLIMcbkAIiI3RizFWvUtNeb3H0yAF/t+KrCtntHdQPg4tcWk1tQWHkDE1+zfi57FXb8UCcxKqXqRlMoM14XnJYZr/JDIl9j1WH6E9btphOArVxZjwbh6TLjxfZefQ1is5H/woNcPOtiurTowjeTvimzT0pmHsOf+omsvEJeu3ogF/RtU3lji1+GBQ9DeEe4Z13Jw3KllHNaZtwzPF1mvFLGmIuNManGmGnAw8A7WNVkvZcxZC1bRsesQAbHDGZX2i6y8sveYooM9ufd662/JF7+cUdlrViG3W39TN0Hybuq3k8ppRoBtxKFiNwrIrEAxphfjDGzjDF5ng2tcQnoZU2+l/X7WrpHWHfZnln1DLN3zy6z35DOLenVJoytR9JZtju56gYvcQy8e2MYFFRelVYp1Tx4W5nxYmHAPBH5VUTudFR29WoRV18FwJHHH+e8DucRbg/nmx3f8NSKpyrse9d5VgHBR77dyJUzlvGXz9dRVFTuFl/8GAhuBQU5kOO1YxWV8jh3bperU9z5/ty99fSoMaY31hSlbYFfRMSrn8z6d+gAIhSlpTE4ZjC/XvErE7tOJDU3lYMZB8vsO6J7K0b3ak14oD8HU7P5ck0iSZnlrhoCw+Gcv1mviwrq5ySUauICAgJITk7WZOEmYwzJyckEBAS49LnaVo89BhwBkoFWtWyrURNfXyKuuooTH39M7o4d2OPj6RlpPQxad2wd7UJOzeIa6O/LW9daz4TeX7yHad9t5tftSVx6WmzZRn0cX/9LCXDTD9DGcyMplfJGsbGxJCYm4ukpj5uTgIAAYmNjq9+xFHcH3N0OTAGigS+Am40xm91pqykJGnQaJz7+mIKkJOzx8QyOsR5cb0zeyAWdK+/wNbJna6Z9t5lpszbx1NytPD6pD2N7O/o095oISTusrrJpBzVRKFUNm81Gp06dGjqMZsfdZxQdgT8ZY3obY/7RHJIEgF+5+bg7h3cGYPau2ZXtDkC78EDuHhnPhP5tOZaey8aDaac2BkVCwuQ6iVUppTzF3WcUDxhj1no4liYj/+jRkteh/qGcyD3BtpRtle7r4yP8eXQ3/nlxXwAKyz/UdhQaZOv/6iRWpZSqLVfLjP/m+Fm+3LhXlxkv5tuyJQDZa34vWXdJV2syo2dXPVujNl77eRfpOfmnVrSyut2y+RvISvFInEop5UkuJQpjzHDHz1BjTFipJdQYE1Y3ITYe9k6dELsdKVXU6y+D/kJ8RDzLDi+jsKiKsh0Oo3paz/uX7io1vsLXz3o2kZcBy9+sk7iVUqo2ajPgrl31e3qf0kkCrAqPbYKtUh0n85xfVN0wzHoI997iveU2zLV+/vIvHXynlGp0ajPgbn5zGnBXzBQUcOLDDynMOFVq/Ox21lQc1d1+GtY1ij7twli6O5nPVu4/tcE/CCKtB+McXgdFTua0UEqpeqYD7lwU0M2qEJt/MLFk3YQu1ixW2QXZ1X7+1rO7APDqwl0l06cCcJ41JSLvjIb5D3ooWqWUqj23Z7hzaDYD7oq1vOmPAByd/s+SdUG2ILq06FKjz0/o15ZebcLYn5LFawt3ntrQ7XyrBHlwNGypurutUkrVN3efUdwuIj8DPwJRWAPuEjwZWGMVPHQoAHmO6Qvd8dZ11qjtVxbuPNVd1j8IBlwN+TmQtt/Jp5VSqn65nCjEmp9vEM1wwB2AT3Awgf37U3D4MIUZGW610S48kPAgG0UGdh0v18Zp150q7aGUUo2Ay4nCWNW4BjTnAXfBjikK0776muwNGwAoMAUs2Legxm0UD8C74b2V5BeWenjt46tFApVSjYq7zyiWiojrc/15CXu3eACO/vOfJN45FYBwe7hLbQyKi8BH4GBqNmnZpQbg5TsmZD+8zhOhKqVUrbmbKM4FlonILhFZLyIbRGS9JwNrzMLGjKHrzwsJu2gCBceOYfLyGNR6EDYfW43baBUawKMX9a64ocu51s8Vb8GBFR6KWCml3OfuzfDzPRpFE2SLicE3JBSAnM21e0Tz09ZjXDygHTZfH4juAX4B8PuHcGwL3PyjJ8JVSim3uXtFcV0VS7MSct55AGSvc+82UWSwNcr7b1+sZ0lxWY/ITvD3ROg6Sp9VKKUaBXcTRWappRDrCiPOQzE1GcXPKgock6jkF+WTW1jzEhzjE9rw3g3Wo568glIPtH1tIL6Q6/V1FpVSTYC7I7OfK7VMB0YALtd+EpFxIrJNRHaKyANV7DNCRNaKyCYR+cWdeOuKrZU1xjD122+xeg3DmqNrXGojOsRe+Ya8TEjZrfNpK6UaXG1HZhcLAjq78gER8QVexboa6QVcKSK9yu0TDrwGXOQoGXK5R6L1IJ+gIEx2DufEDANgc7J7zyuW704uu6LDGdbP1e9BfvWlQZRSqq64OzJ7g6O303oR2QRsA15ysZnTgZ3GmN3GmDzgU2BiuX2uAr4yxuwHMMYccyfeuhQyaiRFGRmEPvkuAIcyDrn0+bioYABOlp6jAqD9EOvngkdg98+1DVMppdzm7hXFhcAExzIGaGuMecXFNtoBpetgJFLx9lU3IEJEfhaR1SJybWUNicgtIrJKRFbV96Trre65B4Ci9ZuJsEeU3IKqqRC7H61C7Xy+KpH5m46c2tBtDPzRUWdx4XTYt8RTISullEvcTRSnAynGmH3ADcDnIjLQxTYq+41abp5Q/IDTgPHAWOBhEelW4UPGzDDGDDLGDIouN691XbO1a4c9visFR49iyysffs38aZR1Sr8fSC27oXUviB8LRzbCtu9rGalSSrnH3UTxsDEmXUSGY/0C/wB43cU2EoH2pd7HAuXv2yQCc40xmcaYJGAR0M/NmOtM6JixAASn5/HZts9Iy01z6fNXndEBgNd/3lX2FpR/MFz9uTWuYtW78PIAOLTWU2ErpVSNuJsoiidSGA+8boz5FvB3sY2VQLyIdBIRf+AKYFa5fb4FzhIRPxEJAs4AtrgZc52xtYkB4MYlAYDrzykARvawelAdSs3GKqdVyqh/QNeRVi+oo5tqF6xSSrnI3URxUETeBCYDc0TE7mpbxpgCYCowD+uX/+fGmE0icpuI3ObYZwswF1gPrADeNsZsdDPmOtPi4otBhI5BsQBMnj2ZrPwsl9q4qH9bAMa9+CtPzd1WduOQ22GsY/6L9Z/VOl6llHKFu4liMtYv+HHGmFQgErjP1UaMMXOMMd2MMV0c4zEwxrxhjHmj1D7PGGN6GWP6GGNedDPeOiW+vvh36UxkQARxYXEApOSkuNTGyJ6tmTahF5HB/mw/ml5xhzDHc34RSNqp06UqpeqNuwPusowxXxljdjjeHzbGzPdsaE2Pr48ftyTcAsBTK55y6bMhdj+uH9aJQJsvP22tpBewCLQdYHWVfeU0WDHDAxErpVT1PDXgTjmM7jgagA1JGygyrv/VH986BJ+qethe8hZcZo3X4OgGNyNUSinXaKLwFAP5Bw8S4BdAz8ieJOcks/TQUpeb6R4TalWRrUxUPPS51Hq9e1EtglVKqZpzKVGIyIeOn/fUTThNV+GJE+RstJ6z/2XQXwD47eBvzNk9h1VHVrnUVm5BEUVFTsZktB1gzaudvMvteJVSqqZcvaI4TUQ6AjeKSISIRJZe6iLApiLk3BElr+PC4hCEj7Z8xP2/3s+N824kMz+zRu3k5Fk9j1fvP1H1Tm0cQ0nmPehmtEopVXOuJoo3sLqr9gBWl1tc+7PZy/hFtkRs1gx3rYNb89Pkn5g1aRY3970Zg+FkDUuGX9C3DQBzNhxmy+EqPjP+eQiL1bIeSql64erYh5eNMT2Bd40xnY0xnUotLlWP9XZRgVF0atGJtiHW+IjlR5bX6HPRoXZ8fYT3Fu/lwa+reGDt4ws+PpCbBgV5ngpZKaUq5W732NtFpJ+ITHUsCZ4OzFsMbTsUgPl75/PympdJTE90un/n6BDWPDyaYV1bkl/o5DnFQEd9xMUvwqr3oPxobqWU8hB3y4zfDXwMtHIsH4vIXZ4MrCmq7Fd1hD2CmOAYlh1exlsb3uLbXd9W206LQBt2P1/nO0V2AcSqLDv7T/pgWylVZ9ztHnsTcIYx5hFjzCPAEOBmz4XV9JiCAsjPJ3dX2V/YQbYgFly2gDV/WIOP+Lg0tmLDwbSKdZ+K9bkEHjoKFzsG3iWucDd0pZRyyt1EIZwqDIjjtWsTMXiZgN69Acg/dNgj7eUXWgnlWLqTObj97BARZ73WkdpKqTribqJ4D1guItNEZBqwDHjHY1E1QbZ21kPrzN9+q3KfIlPEjPUzKCgqqLa98/tYvZ8e+mYjx07mVL1jhzOgw5lw6HdIc/78Qyml3OHuw+znsSYsSgFOADc01oJ99SWgRw8ATH7VvZAiA6yhJseyqp/RtX/7cHrEhLJg81HWOBtTARBmJSnmP1yzYJVSygVul/AwxqxxdJd9yRjzuyeDaop8AgPxjYhwus+fBv4JgJ/2/1Rte73ahvHClP4APDNvG7/tSKp650mvgX8oZB6HY1u1sqxSyqO01lM96tmyJwA/J/5co/3jWgZzYUIbdidl8usOJ3OB+9mhVQ/Y+yu8dgYse9UD0SqllEUThQeZwkJO/OcTinIqf6bQI7IH/aL7sfzwctYfX19te4H+vrxy1UDsfj78b8Nhbp65iv3JVUyIdPGbcPkHID6QXc2tKqWUcoG74yimiojz+yzNkK21NZ1p3u7dVe6TEG2NTfzP1v/UuN1LBsYSYvdjweajrNhbxYRILbtA70lWolBKKQ9y97dKDLBSRD4XkXEi0qy7xhaL/tOfADj20ktVjn/42+C/ERsSy/92/48D6Qdq1O4/L+7LW9cOAmD57mSPxKqUUjXlbq+nh4B4rC6x1wM7ROSfItLFg7E1Of6drXJXmb8souBY1T2buoRbX9OnWz+tcdvRoXYATubkO9+xqAB+fQ7ys2vctlJKOVObXk8GOOJYCoAI4AsRedpDsTU59k6diHnsUQCyVlQ9Uvrf5/0bu6+dRYmLqh55XU6AzZfurUPZfPgk02ZtYvfxjMp3jHGU3co46lLsSilVFbdrPYnIauBpYDHQ1xhzO3AacKkH42ty7F27ApD2TdU1nUSE3MJc9p7cy/70/TVu+7S4CDJzC3l/yV6+XXuo8p2G3G79nPM3+O8N1jLvQS0aqJRym7tXFFHAJcaYscaY/xpj8gGMMUXAhR6LrgkKGjiQgL59yVy8mKLMqicremLYEwDsP1nzRPHPi/uy5mFrTu7V+6ro2RSTAK37QMpuOLLB6jK79BXIqH6Qn1JKVcbdRGE3xuwrvUJEngIwxmypdVRNnE9AAACZS6ueMzsiwOo0dsePd9R49rvSjqVXUdYjpg/cvhjuWmUtI/5urf/9Q1g+AxJXu3wspVTz5m6iGF3JuvNrE4g3af3QQwBkLPqVE599TsYvv1TY58y2Z3Je+/MAyClwUsupEqN6tmb70Qy2HqnBrHnF5T1+ehy+vw9mNftq8EopF/m5srOI3A7cAXQWkdIjxkKxnlUowDciHGw2Uj//vGRdt1Ur8Q0JKXnv5+PH0LZD+elA9eU8youNCATg6bnbePf6wc537n4+3L8PigqtJHF8C2SWKgdiCwL/IJdjUEo1Hy4lCuA/wPfAk8ADpdanG2OqGAnW/NhataLbksUUZWWT+uUXJL38b0x2NpRKFLXxjwm9+G1nEtuOpNfsA4Hh1k//YOvZxTOlejH7h8BftoI91COxKaW8j0uJwhiTBqQBV9ZNON7DNzTUWkKsX8CZS5bQYuLESve99vtreXXkq8S1iKtR2yJCXkERB1OzyS8swuZbwzuI5/4d2p9+6v3+ZbDxC8g5qYlCKVUll55RiMhvjp/pInKy1JIuIjW4Yd78hJx3LgA527eTuXwFRbmnJiI6O/ZsRncczf70/exJ2+NSuxcmWPNVrE9M40BKFfWfyovsDKfffGrpdLa1PnmnS8dWSjUvLiUKY8xwx89QY0xYqSXUGBNWNyE2bT5BQSBCyjvvsv+660iZObNkW9uQttzU9ya32i0eqX3p60s46+mFHE5zYyR2ZCfrZ9J2t2JQSjUPrj6jUC7yi4yk83ezKDxxgn3XXkdRVg3/+q/GFYM7ENcymGW7k3lz0W7Scwpo08LFRqKtsucc3QhbZp9aH94B2iR4JE6lVNPnaq+ndMBQ+fzYxtWrChEZB7wE+AJvG2P+VcV+g7GmW51ijPnClWM0BsWjtXFSO9GVEdpglSA/t0crcgusqcunzdrEf24e4lpgtkDw9YfV71tLMf8Q+L+DrrWllPJarj7M9tgTTxHxBV7FGpORiFWNdpYxZnMl+z0FzPPUsRtMUREFx8tOQBQVGAXAtpRtbjV5Vnw0ANn5ha5/2B4C96wr21125Vvw+0duxaKU8k6eeph90o2H2acDO40xu40xecCnQGXdgu4CvgS8ogZF2hdfUpiaWvK+VVArIgMiWXd8ncsD7wCC7X6cFR/F7/tT+XxVzcqWlxHW1rrNVLyExIApgifbn1qW/Nv1dpVSXsNTD7PD3HiY3Q4o/Zst0bGuhIi0Ay4G3nCx7UYp4qqrACrUgGphb8H+9P3M3Dyzso9V65KB1te2YLMHKsb2uwKGToUB11gLAofW1r5dpVST1ZDToVX6nKPc+xeB+40xTu+riMgtIrJKRFYdP+5kbukGFtC3LwDHXnypzPq3Rr8FwMlc93oYXzwglp5twth9PINXF+7kUGot5qJo2QXGTodxT1pLUKQ11uLEvuo/q5TySu6WGQ8QkT+LyFci8qWI3CsiAS42kwi0L/U+FihfO3sQ8KmI7AUuA14TkUnlGzLGzDDGDDLGDIqOjnYxjPoTMuIcAAqTk8qsbx3cGoAPNn/AyTz3kkXnqGB2Hc/kmXnb+HSFaw/GnWrpeBC/Yobn2lRKNSnuXlHMBHoD/wZeAXoCH7rYxkogXkQ6iYg/cAUwq/QOxphOxpg4Y0wc8AVwhzHmGzdjbnB+EREEDhhA5pKlpH03u8y2sXFjAUjLSXOr7VeuGsCO6efj6yMUenLuias+B78Aq1aUUqpZcjdRdDfG/NEYs9Cx3AJ0c6UBY0wBMBWrN9MW4HNjzCYRuU1EbnMzrkYv/DJrXqfyFWXPibWuNv78y5/53+7/udyuiGDz9an0fl6t+PhYiWL5G/BYS2v5YIKnj6KUasTcTRS/i0hJp30ROQM3qscaY+YYY7oZY7oYY6Y71r1hjKnw8NoYc31THENRXvill2Jr356Ts2eTvWlTyfr+rfozInYE+07uY+GBhQ0YYSUufAHO+jMMuwda94ajm6r/jFLKa7jaPXaDo7z4GcASEdnreH6wFDi7DuLzSsHDzgQg7dtT06W2D23Pv0f+m9ZBrVlyaInbbRcUGV5duIuM3IJax1mizyUw8hFriR0Mhfmwd7E1CVJRkeeOo5RqlFy9orgQmACMAzoB5ziWTsB4z4bmvdpMm4YEBpI+fwGm3POEtNw00vNqWD68EsO7WgP4TmTm1SrGKvkHQ+5JeP8CePs82D63bo6jlGo0XB2ZXdJHUkQigHigdG8n7UNZQyYnh4LsbPL27MXeuVPJ+knxk/jPlv+43e6kAe34bWcS//f1BoL8fQFoHRbAtAm98fHxwBOMc+6HrqMg7SB8cxvkZdS+TaVUo+Zu99ibgEVYD6Ifdfyc5rmwvF/bp58GIP9Q+R7BtdO3XQsSYltwPD2XfclZrD2Qysyl+0jKzK3+wzXhH2yVJy+e12LXT9Zc3Kves+a1UEp5HXerx94DDAaWGWPOFZEeWAlD1ZBPqDXb3YGbbqL76lX4BAeXbMstzKXIFOEjrufx7jGhzJo6vOT9R8v28dA3G5m19hChAX70adeC3m1dLTNbicAI8AuEdZ9YC4CPLwy8tvZtK6UaFXd7PeUYY3IARMRujNkKdPdcWN4vZNgwQkaMAKAo51SNp+J6T+uPr6/sYy4rnrfiif9t4f4vN/DX/3qmXYIi4W+74b7dcMcya11hvmfaVko1Ku5eUSSKSDjwDbBARE5QcVS1ckJsNoLPPouMn38us350x9F8svUTdqXuoltEN4JsQbU6ztjeMax4cCQFhYaHv9nInuRMiopOPUAXscZguMU/yFqKHD2sNn1dcRKkdoMg4XI3o1dKNQZuJQpjzMWOl9NEZCHQAtDuLx5QnBimLZ3G/H3zeXP0m7Vus1Wo1d8g2O7H7uOZdP6/OSXbokLsLPrbCIL8azGHVUALiO4BR9ZbS7G8LNg8SxOFUk2cW78dHHWd7gCGYxXy+42GLTDYpCXeOZV2zz2LrV07ekb25PVRr/PympdJzU316HFuH9GFLtEhJe83HEzjhy1HScvOr12isAXAncsrrp91F+xY4H67SqlGoSFrPTV7wUOGEHTGGWSvXUvONmviIh/xYXi74UQHeb64Yc82YdwzKr5kGdWzFQCHUl2fB6PGspLh48lwbGvdHUMpVacarNaTAnvnzrT6231Vbs8rrKNBcw6tW1i3pDYedK8QYbXix1olP3bMgz2L6uYYSqk616C1nlTVcgtz2Zm6k6z8rDo7Rt92VjfZg6nZrNl/omQ5nFaL+SxK63khXP2lZ9pSSjUYl25Mi8gGrGcSNuBaESme+KADsLnKD6pqFRwrO9Nrt4huLD+8nNTc1Fr3fKqK3c8HH4EZi3YzY9HukvVhAX6snzbWswf77QXroXe/KZ5tVylV51x9gnlhnUTRjPmFhwOQtXIVEVdcUbI+PjwegPFfj+fD8z+kT1Qfjx87NMDGd3cN53j6qVHbX645yP/We7Cnc1AkDL4Z1n4MO3/QRKFUE1SbWk/9gLMcb381xqzzZGDNha1dO3yCgsjbuxdTUID4Wf9JzutwHnvS9vDepvdITE+sk0QBVBilvXrfCYoMnPX0TyXrbjunC1ef0dG9A4jA+GetJKGUapLc7R57D3Az8JVj1UciMsMY82+PRdaM+EVHk7NpE8nvvUfUzTcD0MLegoldJ/Lepve4b9F9nMw7yeTuk+s8lvEJbTicllMyKG/epiMs353ifqIoJj7W3Ntbvqt6n/hRMOWj2h1HKeVx7nae/yNwhjEmE0BEnsKak0IThRvav/kGu8adT+GJ1DLrO7XoxMNDHubJ5U+SmJ5YL7H0iAnj2cv7lbw/55mFzFp3iPvP70G78ED3Gx7zOOxfVvX2HQvgkF6UKtUYuZsoBCg9iXKhY51yg39cHAAp775L1G234hsWBlhjKiZ3n8yzq57lk62f8P3e73lxxIv0jupdb7F1iAxiX3IW7y/ew4Pje7nfUI/x1lKVzCTY+5v77Sul6oy7ieI9YLmIfO14Pwl4xyMRNVOhY8eSPm8ehampJYmi2F9O+wtrjq1hzp457EzdWa+J4v0bTqf7Q9+TW1BETr71t4G/r49n5rYo72QivODkWUyviTB2uuePq5Ryyt1aT8+LyM9YJTwEuMEY87snA2tuQs87l/R589j/x5toefNNREw+9TxiSo8pDGs3jDl75vBL4i9M7Dqx3uLy9RFCA/yYuXQfM5dafRkGdAjn6zuGefZAp11nPfiuyu6fddCeUg3E5UQhVqnRWGPMGmCN50NqnoLOOIMWEy8ifcEPZC5dWiZRALQJbgNQYerU+vD85P5sPWJNzzp/8xH2JdfBIMAOQ6ylKv+5Ak4e9PxxlVLVcjlRGGOMiHwDnOb5cJovW0wMbZ96il3rzid97jx4oex2Xx9fuoZ3ZeXRldww9wbuGnAXA1sPrJfYzu3RinN7FNeFymbbkXQ+XLYPu58PExLaEuiYcrXOndgHn15ddl2rnnDeQ/VzfKWaKXdLeCwTkcEejUQBYPKrnvxnUtdJ9Ijowaqjq1h6eGk9RnVKTIsAsvIKefibjfzti/X8tPVY9R/yhPjR0CIWUvacWvYvg1+fq5/jK9WMuZsozsVKFrtEZL2IbBARD02d1ryFnT8O8fevdNt1va/j7bFvA/DT/p8q3aeu3TGiC6seGsWXt58JQH5hUf0cePAf4Y4lZZdBN9TPsZVq5tzt9XS+R6NQLpMG6o0sIkSF2EnPsWa1e2XhTj5beQCbnw8Pj+9JfOvQ+g3IFMHilyqu9/WHfldCYHj9xqOUF3I3URyl4sRFr3sqqObO5OZa3WQddaDKG9F+BD8f+Jmvd3zNxfEXV7pPXWsbHsC43jGkZOaRnV/I0t3JLO/Vun4TRbhjtPiCRyrfbg+FAdfUXzxKeSl3E8VMIJ1TI7GvxJq4SOe8rCWfUGsMRfamTYQMq7wL6ti4sfx84Gfm7p3bYInC7ufLG3+w+jMcS8/h9Ok/kplbQHJGbpn9QgNs+PvV0eSHA/8AfS+zripKSz8C/x54ai5vpVStuJsouhtj+pV6v1BEtP6CBwQNHgTA8RdfIj/xIBFTKtZ3urDzhXyy5ROWHFqCMQZxNv6gHth8rETw5PdbefL7sjPZ9Wsfzrd3enjMRZmDV1JWpHhdym44ugla9XI+RkMp5ZS7ieJ3ERlijFkGOnGRJ/l36kTQkCHkbNxI6uefV5ooAAzWeIq8ojzsvvb6DLGCiGB/3rjmNI6ll51S9ZvfD3Iw1UOTILnCLwAQ69nF4pfgjmVWN1qllFvcTRRnUHHioi3FExsZYxI8El0z5BcRQcf33+PArbdRkJRU5X7ndTiPDUkbuHHujYQHhPPsOc8S6FeLon21NK5PTIV1mw+dZPvRDN5atJvLB8USHlR5by6PC4qE236FHfPhx8cgN71+jquUl3I3UYzzaBTKZWe2PZPVR1dzNOso6xPXcyjjEF3CuzR0WGV0bBlMRm4B0+dsIdDfl2uG1LJUuSti+kLG0fo7nlJezK2njMaYfc6WmrYjIuNEZJuI7BSRByrZfrVjnMZ6EVnimCypWTCFheRs2kRRTk6l23u17MXro17nloRbAPhw84f1GV6N3D6iC8v+PhKAwqL6Lz2COEaMvzMa/vfX+j++Ul6ijrqjVE9EfIFXscZk9AKuFJHydaz3AOc4bmU9Dsyo3ygbjq2NVdspb+9ep/sNiB4AwI7UHXUdkluKezy9snAnX6yunzk1SnQYAqOmQYv2cHxrtbsrpSrXYIkCOB3YaYzZbYzJAz4FypRFNcYsMcaccLxdBsTWc4wNJnCgVcdpz6SLyVq1qsr9Wge3Zmiboaw/vp4vtn9RX+HVWESQjRuHdSIjp4Dfdhyv34PbAmH4vRDeoX6Pq5SXcfcZRQUiEmOMOeLCR9oBB0q9T8R6SF6VPwLfV3HsW4BbADp08I5fCmHjxlKYnMSxZ58j/YcfsXftWuUAvEviL2Hp4aUsSlzEZd0uq99AqyEiPDKhFz9uPcqRkzksdFIbqkt0CB1aBnk+CGPgyHrIywL/OmhfKS/nsUSBNXGRkynMKqisY3ulN7JF5FysRDG8su3GmBk4bksNGjSoAW6Ge55PYCBh48dz7LnnSXn/fXyCAom+++5K9x3XaRxvb3iblUdWcsXsK7hn4D0MbTu0niN2LjTAj2W7U1i2O6XKfeJbhbDgz+d4/uD+wZCbBnPvh4t0tl6lXOWxRGGMcSVJgHUF0b7U+1jgUPmdRCQBeBs43xiT7H6ETY+tTRviF/3CzpGjyFy2HJ+gtwmfPLnCDHgAk7tPZlHiIn5J/IWVR1Y2ukQx88Yz2J9S9TwWL/2wne1HM+rm4BNfgee6Q87JumlfKS/nVqIQkaeMMfdXt64aK4F4EekEHASuAK4q12YH4CvgD8aY7e7E2tT5RUfj36UL2WvWkL1mDX6tWtHioosq7De5+2Qmd5/MgJkDGiDK6kUG+xMZXPU4ishgO6lZKTz49YYq9zkrPrrS8RrVCo2BqO6uf04pBbh/RTEaKJ8Uzq9kXZWMMQUiMhWYB/gC7xpjNonIbY7tbwCPAC2B1xxlKgqMMYPcjLnJ6vTVl+Tv38+usePIXLqs0kTR1A3oEM4v248xb1Plj7lOZhewLjHVvURR7NAa+HZq2XU+fjDsbojs7H67Snk5lxKFiNyOVTW2S7n5J0KBJa4e3BgzB5hTbt0bpV7fBNzkarveRkTwi7F+QRZlOB9lXGAKeGvDW9zY50ZC/EPqIzyPuGZIR6cD8m76YCWH0yofU1Ijnc6CrXNg54+n1pkiyDhiJYlhlT//UUq5fkXxH6yeR08CpQfIpRtjqn5KqWrNx27H3q0bGT//wo6zziZm2j8IHTmywn5ntTuLXw/+ypaULcQEWcklzB5GC3uL+g7Z43Ycy2Dkcz+XWdc9JpTXrq7BrLzjn7OW0vIy4Z9tqaIPhVLKwaVEYYxJA9JE5CsgxRiTLiIPAQNF5HFjzO91EqUCIOqOO8hctpTUTz8jc9nyShPFhC4T+PXgr9w478aSdXZfO79M+YVgW3B9hutRV57eAbut7Nzc246kM3ejKz2ylVLucPcZxcPGmP+KyHBgLPAs8AbOx0GoWgobN5awcWNJ/fQz0ufOJebB/6uwz3kdzuOZc54hv9Cae3vZ4WXM2jWL+XvnE+YfRnxEPB3Cmt5Yk5E9WzOyZ+sy656bv41XFu7kQCW9qWy+PsS0CKiv8JTyau4mikLHz/HA68aYb0VkmmdCUtXx79yZvN27yVq1iqBBZZ/t233tjIs7VbPR5mNj1q5ZPLLEmgWuZ2RPPp/web3GW1f8fX0wBs56emGl22f84TTG9K7Bw+/CPOs2lNtEB/Ipr+ZuojgoIm9i9X56SkTsNGw5kGYlbNw4kl57jaQ3Z9BhkPNOYGPjxtItohv5Rfk8v/p5tiRv4a31b5VsjwqMarBZ8mrrD0M7EhsZSGG5Ce7SsvN5fPZmkjPznDcgjn+yPz1hLbUx/nkY/MfataFUI+VuopiMVWr8WWNMqoi0Ae7zXFjKmei77yLjt9/I/PVXTFER4lN1jhYROodbXT97t+zNkkNLePn3l8vsc1bsWUQFRtVpzHUhPMifiwdULP91JC2Hx2dvZuHWY+QXFnHV6R3w863kO7IFwuXvQ+r+ittc8cOjkFrjoslKNTnuJopsIBhrruzHABuQ6qGYVA34BFu3OtJ//BG/qCiCBlQ/0O7ugXdze7/bS95/vfNrHl/2OHvS9jTJRFGVsEA/IoJszN98lPmbj9IvNpx+7cMr37m3B66mFv6z9m0o1Yi5e7voNWAIVqIASMcqGa7qSdiYMQAcvOtu9l15FXmJB2v0OZuvrWRpH2pVUNmZurPO4mwIQf5+rHpoNG9fa92WK2iIuTCU8iLuJoozjDF3AjkAjlLg9TTPpQIInzyZTl99SfSf/gTAiY8+crmN7pFWWYsXVr/A8E+HlyzXz73eg5E2DF8fweaYC+O2j1bz8o+Nc74OpZoCd2895TsmHjIAIhINFDn/iPIk8fUloFcvfEJDOf7ii6TNmkWrP9+L+Nc8X0fYI7hn4D0cyzpV+nv98fWsObqG5YeXl6zr1bIXof6hHo2/PvRvH841Qzowd+MRluxK4u6R8XV3sM3fwvFtddd+sV4Tof9V1e+nlAe5myheBr4GWonIdOAy4CGPRaVqzL99e0LOPZeMhQs5/OijtJ0+vcafFRFu6lu2Qsq7G99lU/Imbpp/av1l3S7jH0P/4bGY60uLQBtPTOpbd1Vpi/W93JrvIv1w3R4neRfkpGmiUPVOjHHv/q2I9ABGYs0r8aMxZosnA3PXoEGDzConM8J5o7zEg+waNQpb27YEn30Wkddcg71rV/faKsxjU/ImCousoTL3L7qfY9nH+GLCFyW3qpqaKW8uZfmeFB68oCc3n92Ei/99MAEK8+HGuQ0difJCIrK6qqKrbo99MMZsNca8aox5pbEkiebKP7YdEX/4A0X5eaR++hlps2e735avPwNaDWBQzCAGxQxieKw1V9Q3O7/xULT17/YRXQBYeyC1YQNRqonSQXJeIubB/6Pbr7+Cn59Ha9w9euajhPmH8dGWj+g/sz/9Z/bnpnlNq6DviO6t6NoqhLTsfDYeTCuzbD1ykiLtFaWUU56cClU1BgUFJL/5JtFT70RsNo80+cjQR9iWYj2o/e3gb2xM3sgvB37B39efQTGDsPl45jh1KdDmy287k7jw379V2PbkJX258vSmUP9KIHEVPNez+l2DW8KN86xpYJWqJU0UXsavTRsKDh8mZeZMfCMirXXR0YScVel04zUyNm4sY+PGAmAwbNmwhak/WRMAvTjiRUZ2rFjFtrF55aoBbDtSdi6P3IIi7vrkd05m5zdQVC4aOhXCa5DQUnbDvsWQcQwiO9V9XMrraaLwMm0encaBW27l2DPPllnfZcF8fMPC8AkJQXx9q/h09e7odwejOo7icMZh7v35Xt7b9B7z9s3D5mPj9n63ExtasaRGY9CxZTAdW5b96zorr6CBonFTtzHWUp21n1iJQikP0UThZULOPpuui37B5Fl/JafP/Z5jzz7HrtHWL5jg4cPp8PZbzppwyuZro3fL3nQK60T/6P6cyD1BUnYSBzMO0i2iG5O6Tiqzf5AtqNHfmkrKyGXHUeczBzojInSOCsbHRzwYlVKNh9vdYxur5tg91pnC9HROzp6Nyc8n9etvyNu9m4BevUq2+0VH0+65Z2v1PONEzgnO/uzsSrclRCfw8QUfu912XcrJL6T3P+ZR6IGH2feP61HSu6rBrf0EvrkN7l6rt55UjTnrHqtXFF7ONzSUiCutklx+MTGkfvppybb8o8fI/v13crZvJ7B3b7ePEREQwXPnPMfx7ONl1s/ZM4cdJ3bwrxX/4vre1xMTXIO5IepRgM2Xz28dyuG07Fq186dP15KaXU1J84bw3+utCrnKu/nZ4cIX6/SPAk0UzUjYmDElxQTBqjybeOdUTsz8kLyzzipZ798+lsB+/Vxqe0xcxXvndl87L695mY+3fEyofyi3JtyKn0/j+id3WscIIKJWbfzl83WeCcZT2p8OXUdDYW5DR6LqWl4W7F8KB1fXaaLQW0/NWFFmJtvPHIbJLfsLRQICaP/mmyXvA3r3wjckxK1jZORlcOYnZ2IwXN3zah44/YFaxdwYdX/oe3rEhDKkc8s6P9Y53aI5s6v3lIRXtXR8O7w6GC59B/peVqum9NaTqpRPcDBdF/5EYWpaybrUL78g5Z132X/ddSXrAgcOpOPMDxA/1/+5hPiH8M7Yd/jrL3/lm53fsOTQkir3HdJmCP93RsV5wBu7vu1asPFQGttq8UC8JnILivj9QKomClXvNFE0c36RkfhFRpa8j77nHkLPPRfjmF/0yCOPkL1mDQf/eh/Rd95hfaZ1a3zDwmp8jMExg7mz/52sOLKiyn02HN/Ad7u+I8RW9ZXL4JjBDG07tMbHrS9f3H5mvRznqreWkV9+3lel6oHeelJO5e7Zw+7zLyizToKC6Pz1V4jdji3GMw+o39v4Hi+teanK7YWmkLiwuApVbGOCYxrt2A1Pu+qtZexNyuTyQe2r3TcqxJ9rhnRERLvserV6uvWkiUJVK2frVvL27gXg2DPPkn/w1Gx6rR64n5DhwxGbDVuHDnX2i+lvi/7G93u+r7A+1BbK3MsqVlP1Ez+CbEF1EktD+ce3G/lgac3n5v7t/nOJjfCu70CVU5woorpDcLQ1X8kZt7jVlCYK5TEFyclkLllC/qHDHH/hhTLbwqdMIez8cWXW+XfogK1t21ofNy03ja0pW8us+27Xd3y769tK9xeE10a9xvB27pcuaaq+XJ3IX/67jkX3nUuHlpoovFp+Dnx1E2SdgKMboWVXuPlHt5rSRKHqROayZRSeOEFhejpHHql6YqPWDz9EiwkTXHquURNJ2UnM2zuPIlP2vn1mfiavrrWmcA/zD+PLi75sdGM46lJxorjqjA6EBTTuUfGq9ux+Ptw4vBMtvphiTWxVB4lCH2YrtwUPGVLyOmjQIAqTk8tsT3rjTTIXL+bo409w9PEnaP/224QMH+ax40cFRnF1z6srrDfG0MLegtm7ZrM+aT2jvxjNlO5TeGhI85iEMS4qmNAAP75YndjQoag6Zowhv9AQ3zqEC+vwOHpFoepUwYkTJN51F9mrVle7r09QEG2eeLzK7QG9euEfF1fjY+cU5PDJ1k94fvXzAFwafym3JtxaZh9fH19aBbWqcZtKNSY7j2Uw6vlfePnKAVy0fmqdXVFoolD1ImPRIrLXVj2COenNN6GwsNp2gkpdxQD42O1ETb2zwn5is2Hv1g3x8WHVkVXcMO+GKtt8YtgTXNi5dn+P+fq4X5FXKXcVJ4qxvVvzQPKDRPnlEDr1F7faarSJQkTGAS8BvsDbxph/ldsuju0XAFnA9caYNc7a1ETRNJmCAvL2Vd2jJ+3bWWStWgWlelVlr67+KsXeowcRV0zhSO8YtgacKLMtvyifx5dVfQXjiiu6X8Gl3S71SFvOdAjt4HW9uZT7UrPyGPncL5zIyuM9v3/RJiCXbg9WPV7JmUaZKETEF9gOjAYSgZXAlcaYzaX2uQC4CytRnAG8ZIw5w1m7miiaj6KsLLJWr8aUvxIpLOT4q69ScPgIhSdOJQex2QgZVXaSpcMnDnAsrgVZsZG4a+7uOeyLFrLtbjdRY+mBMCF+YrX75RXmcXXPq/ERne3Y29l97XSL6MaGp0ZhL0in+0OeTxQN+TD7dGCnMWY3gIh8CkwENpfaZyIw01jZbJmIhItIG2PM4foPVzU2PkFBhJQqZlha6EgrIeQfOULGokUkvfY6PgEB5G7bXrKPKSggZP9+QpbVLo4+tfu4y3JtX1e+wXGxVfzHn2E2xSnUOBneUuWfiu58pppj1Vt7jf18PRSfGLDnwwqBjEBhXfd2dHc9nGo1ZKJoBxwo9T4R66qhun3aAWUShYjcAtwC0KFDU5j7WNUXW0wMEZMnEzF5cqXb848eozAludJtNZW9cSMmv+6nU83buQux28vcfitR6s5AkSniaNZRCo2VJsTpbzlT+e/Hau40VPl7rqrPVRODk42eba/KTc4+U/W2Kr9bd+/UVPn9VbG+oJCT+3eRHxFC5Np9xISEu3fcajRkoqjs31r5b6Mm+2CMmQHMAOvWU+1DU82FrXUrbK1r1+spoGdPD0XjOW0aOgDVIAbWUbsNeQMzEShdtCYWOOTGPkoppepQQyaKlUC8iHQSEX/gCmBWuX1mAdeKZQiQps8nlFKqfjXYrSdjTIGITAXmYXWPfdcYs0lEbnNsfwOYg9XjaSdW99iqO8MrpZSqEw1awsMYMwcrGZRe90ap1waoOJpKKaVUvdFO1koppZzSRKGUUsopTRRKKaWc0kShlFLKKa+rHisix4GazxdZVhSQ5MFwmgI95+ZBz7l5qM05dzTGRFe2wesSRW2IyKqqimJ5Kz3n5kHPuXmoq3PWW09KKaWc0kShlFLKKU0UZc1o6AAagJ5z86Dn3DzUyTnrMwqllFJO6RWFUkoppzRRKKWUckoThYOIjBORbSKyU0QeaOh4PEVE2ovIQhHZIiKbROQex/pIEVkgIjscPyNKfebvju9hm4iMbbjo3SciviLyu4jMdrz36vMFcEwV/IWIbHX89x7qzectIvc6/k1vFJFPRCTAG89XRN4VkWMisrHUOpfPU0ROE5ENjm0vi1Q2VWIVjDHNfsEqc74L6Az4A+uAXg0dl4fOrQ0w0PE6FNgO9AKeBh5wrH8AeMrxupfj/O1AJ8f34tvQ5+HGef8Z+A8w2/Heq8/XcS4fADc5XvsD4d563lhTIu8BAh3vPweu98bzBc7GmrxuY6l1Lp8nsAIYijVz6PfA+TWNQa8oLKcDO40xu40xecCnwMQGjskjjDGHjTFrHK/TgS1Y/5NNxPrFguPnJMfricCnxphcY8werLlATq/XoGtJRGKB8cDbpVZ77fkCiEgY1i+UdwCMMXnGmFS8+7z9gEAR8QOCsGa/9LrzNcYsAlLKrXbpPEWkDRBmjFlqrKwxs9RnqqWJwtIOOFDqfaJjnVcRkThgALAcaG0cswU6fhZPHO0N38WLwN+AolLrvPl8wboaPg6857jl9raIBOOl522MOQg8C+wHDmPNfjkfLz3fSrh6nu0cr8uvrxFNFJbK7tV5Vb9hEQkBvgT+ZIw56WzXStY1me9CRC4EjhljVtf0I5WsazLnW4of1u2J140xA4BMrFsSVWnS5+24Jz8R6/ZKWyBYRK5x9pFK1jWZ83VBVedZq/PXRGFJBNqXeh+LdRnrFUTEhpUkPjbGfOVYfdRxOYrj5zHH+qb+XQwDLhKRvVi3EM8TkY/w3vMtlggkGmOWO95/gZU4vPW8RwF7jDHHjTH5wFfAmXjv+Zbn6nkmOl6XX18jmigsK4F4EekkIv7AFcCsBo7JIxw9G94Bthhjni+1aRZwneP1dcC3pdZfISJ2EekExGM9BGsSjDF/N8bEGmPisP47/mSMuQYvPd9ixpgjwAER6e5YNRLYjPee935giIgEOf6Nj8R6/uat51ueS+fpuD2VLiJDHN/XtaU+U72GfqLfWBbgAqweQbuABxs6Hg+e13CsS8z1wFrHcgHQEvgR2OH4GVnqMw86vodtuNAzorEtwAhO9XpqDufbH1jl+G/9DRDhzecNPApsBTYCH2L19PG68wU+wXoOk491ZfBHd84TGOT4rnYBr+CozFGTRUt4KKWUckpvPSmllHJKE4VSSimnNFEopZRyShOFUkoppzRRKKWUckoThVJKKac0USillHJKE4VSLnLM+3BHqfdL6ug4sSIypS7aVsoVmiiUcl04UJIojDFn1tFxRmLVa1KqQWmiUMp1/wK6iMhaEXlGRDLAKuPumF3ubcesax+LyCgRWeyYiaxk/gMRuUZEVjjaeFNEfEsfQESGA88Dlzn26VSvZ6hUKVrCQykXOeb1mG2M6eN4n2GMCXGs34k158cmrGKT67Bq81wE3GCMmSQiPbFmKLvEGJMvIq8By4wxM8sdZy7wV2PMRpRqQH4NHYBSXmaPMWYDgIhsAn40xhgR2QDEOfYZCZwGrHRMWxzIqTLRpXXHKuymVIPSRKGUZ+WWel1U6n0Rp/5/E+ADY8zfq2pERFpizdqWXydRKuUCfUahlOvSgdBafP5HrGcPrQBEJFJEOpbbpxNNe2Id5UU0USjlImNMMrDY8cD6GTc+vxl4CJgvIuuBBUCbcrttBaIcx6irXlVK1Yg+zFZKKeWUXlEopZRyShOFUkoppzRRKKWUckoThVJKKac0USillHJKE4VSSimnNFEopZRy6v8BgFovnXDniEIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "pred_surv = estimator.predict_survival_function(x_new)\n",
    "time_points = np.arange(1, 1000)\n",
    "for i, surv_func in enumerate(pred_surv):\n",
    "    plt.step(time_points, surv_func(time_points), where=\"post\", label=f\"Sample {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",
    "\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",
    "\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(\n",
    "    [\n",
    "        (\"encode\", OneHotEncoder()),\n",
    "        (\"select\", SelectKBest(fit_and_score_features, k=3)),\n",
    "        (\"model\", CoxPHSurvivalAnalysis()),\n",
    "    ]\n",
    ")"
   ]
  },
  {
   "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>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>4</th>\n",
       "      <td>5</td>\n",
       "      <td>{'select__k': 5}</td>\n",
       "      <td>0.716093</td>\n",
       "      <td>0.719862</td>\n",
       "      <td>0.716685</td>\n",
       "      <td>0.717547</td>\n",
       "      <td>0.001655</td>\n",
       "      <td>1</td>\n",
       "      <td>0.732087</td>\n",
       "      <td>0.742432</td>\n",
       "      <td>0.731710</td>\n",
       "      <td>0.735410</td>\n",
       "      <td>0.004968</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>{'select__k': 4}</td>\n",
       "      <td>0.697368</td>\n",
       "      <td>0.722332</td>\n",
       "      <td>0.727324</td>\n",
       "      <td>0.715675</td>\n",
       "      <td>0.013104</td>\n",
       "      <td>2</td>\n",
       "      <td>0.732477</td>\n",
       "      <td>0.743090</td>\n",
       "      <td>0.727138</td>\n",
       "      <td>0.734235</td>\n",
       "      <td>0.006630</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>8</td>\n",
       "      <td>{'select__k': 8}</td>\n",
       "      <td>0.706478</td>\n",
       "      <td>0.723320</td>\n",
       "      <td>0.716685</td>\n",
       "      <td>0.715494</td>\n",
       "      <td>0.006927</td>\n",
       "      <td>3</td>\n",
       "      <td>0.739356</td>\n",
       "      <td>0.746249</td>\n",
       "      <td>0.737519</td>\n",
       "      <td>0.741041</td>\n",
       "      <td>0.003758</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>6</td>\n",
       "      <td>{'select__k': 6}</td>\n",
       "      <td>0.704453</td>\n",
       "      <td>0.719368</td>\n",
       "      <td>0.716685</td>\n",
       "      <td>0.713502</td>\n",
       "      <td>0.006491</td>\n",
       "      <td>4</td>\n",
       "      <td>0.735722</td>\n",
       "      <td>0.747565</td>\n",
       "      <td>0.731710</td>\n",
       "      <td>0.738332</td>\n",
       "      <td>0.006731</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>7</td>\n",
       "      <td>{'select__k': 7}</td>\n",
       "      <td>0.700405</td>\n",
       "      <td>0.719368</td>\n",
       "      <td>0.720045</td>\n",
       "      <td>0.713272</td>\n",
       "      <td>0.009103</td>\n",
       "      <td>5</td>\n",
       "      <td>0.741173</td>\n",
       "      <td>0.742564</td>\n",
       "      <td>0.728621</td>\n",
       "      <td>0.737453</td>\n",
       "      <td>0.006271</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>{'select__k': 2}</td>\n",
       "      <td>0.699393</td>\n",
       "      <td>0.717885</td>\n",
       "      <td>0.718365</td>\n",
       "      <td>0.711881</td>\n",
       "      <td>0.008833</td>\n",
       "      <td>6</td>\n",
       "      <td>0.732087</td>\n",
       "      <td>0.727428</td>\n",
       "      <td>0.714409</td>\n",
       "      <td>0.724642</td>\n",
       "      <td>0.007481</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>{'select__k': 1}</td>\n",
       "      <td>0.698887</td>\n",
       "      <td>0.707510</td>\n",
       "      <td>0.712206</td>\n",
       "      <td>0.706201</td>\n",
       "      <td>0.005516</td>\n",
       "      <td>7</td>\n",
       "      <td>0.710670</td>\n",
       "      <td>0.714793</td>\n",
       "      <td>0.700445</td>\n",
       "      <td>0.708636</td>\n",
       "      <td>0.006032</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>{'select__k': 3}</td>\n",
       "      <td>0.708502</td>\n",
       "      <td>0.714427</td>\n",
       "      <td>0.694849</td>\n",
       "      <td>0.705926</td>\n",
       "      <td>0.008198</td>\n",
       "      <td>8</td>\n",
       "      <td>0.734034</td>\n",
       "      <td>0.722559</td>\n",
       "      <td>0.716634</td>\n",
       "      <td>0.724409</td>\n",
       "      <td>0.007223</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  param_select__k            params  split0_test_score  split1_test_score  \\\n",
       "4               5  {'select__k': 5}           0.716093           0.719862   \n",
       "3               4  {'select__k': 4}           0.697368           0.722332   \n",
       "7               8  {'select__k': 8}           0.706478           0.723320   \n",
       "5               6  {'select__k': 6}           0.704453           0.719368   \n",
       "6               7  {'select__k': 7}           0.700405           0.719368   \n",
       "1               2  {'select__k': 2}           0.699393           0.717885   \n",
       "0               1  {'select__k': 1}           0.698887           0.707510   \n",
       "2               3  {'select__k': 3}           0.708502           0.714427   \n",
       "\n",
       "   split2_test_score  mean_test_score  std_test_score  rank_test_score  \\\n",
       "4           0.716685         0.717547        0.001655                1   \n",
       "3           0.727324         0.715675        0.013104                2   \n",
       "7           0.716685         0.715494        0.006927                3   \n",
       "5           0.716685         0.713502        0.006491                4   \n",
       "6           0.720045         0.713272        0.009103                5   \n",
       "1           0.718365         0.711881        0.008833                6   \n",
       "0           0.712206         0.706201        0.005516                7   \n",
       "2           0.694849         0.705926        0.008198                8   \n",
       "\n",
       "   split0_train_score  split1_train_score  split2_train_score  \\\n",
       "4            0.732087            0.742432            0.731710   \n",
       "3            0.732477            0.743090            0.727138   \n",
       "7            0.739356            0.746249            0.737519   \n",
       "5            0.735722            0.747565            0.731710   \n",
       "6            0.741173            0.742564            0.728621   \n",
       "1            0.732087            0.727428            0.714409   \n",
       "0            0.710670            0.714793            0.700445   \n",
       "2            0.734034            0.722559            0.716634   \n",
       "\n",
       "   mean_train_score  std_train_score  \n",
       "4          0.735410         0.004968  \n",
       "3          0.734235         0.006630  \n",
       "7          0.741041         0.003758  \n",
       "5          0.738332         0.006731  \n",
       "6          0.737453         0.006271  \n",
       "1          0.724642         0.007481  \n",
       "0          0.708636         0.006032  \n",
       "2          0.724409         0.007223  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.model_selection import GridSearchCV, KFold\n",
    "\n",
    "param_grid = {\"select__k\": np.arange(1, data_x_numeric.shape[1] + 1)}\n",
    "cv = KFold(n_splits=3, random_state=1, shuffle=True)\n",
    "gcv = GridSearchCV(pipe, param_grid, return_train_score=True, cv=cv)\n",
    "gcv.fit(data_x, data_y)\n",
    "\n",
    "results = pd.DataFrame(gcv.cv_results_).sort_values(by=\"mean_test_score\", ascending=False)\n",
    "results.loc[:, ~results.columns.str.endswith(\"_time\")]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results show that it is sufficient to select the 5 most predictive features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Celltype=large       -0.754714\n",
       "Celltype=smallcell   -0.328059\n",
       "Celltype=squamous    -1.147673\n",
       "Karnofsky_score      -0.031112\n",
       "Treatment=test        0.257313\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": {},
   "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.\n",
    "Before you dive deeper into the various survival models, it is highly recommended reading [this notebook](evaluating-survival-models.ipynb) for getting a better understanding on how to evaluate survival models.\n",
    "The [User Guide](https://scikit-survival.readthedocs.io/en/latest/user_guide/index.html) is a good starting point to learn more about various models implemented in **scikit-survival**, and the [API reference](https://scikit-survival.readthedocs.io/en/latest/api/index.html) contains a full list of available classes and functions and their parameters. 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)](https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.NMF.html#sklearn.decomposition.NMF), before training a Cox model."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}