{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "EZqu9a_ChWYv"
},
"source": [
"# Cox Proportional Hazards and Random Survival Forests\n",
"\n",
"Welcome to the final assignment in Course 2! In this assignment you'll develop risk models using survival data and a combination of linear and non-linear techniques. We'll be using a dataset with survival data of patients with Primary Biliary Cirrhosis (pbc). PBC is a progressive disease of the liver caused by a buildup of bile within the liver (cholestasis) that results in damage to the small bile ducts that drain bile from the liver. Our goal will be to understand the effects of different factors on the survival times of the patients. Along the way you'll learn about the following topics: \n",
"\n",
"- Cox Proportional Hazards\n",
" - Data Preprocessing for Cox Models.\n",
"- Random Survival Forests\n",
" - Permutation Methods for Interpretation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Outline\n",
"\n",
"- [1. Import Packages](#1)\n",
"- [2. Load the Dataset](#2)\n",
"- [3. Explore the Dataset](#3)\n",
"- [4. Cox Proportional Hazards](#4)\n",
" - [Exercise 1](#Ex-1)\n",
"- [5. Fitting and Interpreting a Cox Model](#5)\n",
"- [6. Hazard ratio](#3)\n",
" - [Exercise 2](#Ex-2)\n",
"- [7. Harrell's C-Index](#7)\n",
" - [Exercise 3](#Ex-3)\n",
"- [8. Random Survival Forests](#8)\n",
"- [9. Permutation Method for Interpretation](#9)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "IH0ukiNS3zG-"
},
"source": [
"<a name='1'></a>\n",
"## 1. Import Packages\n",
"\n",
"We'll first import all the packages that we need for this assignment. \n",
"\n",
"- `sklearn` is one of the most popular machine learning libraries.\n",
"- `numpy` is the fundamental package for scientific computing in python.\n",
"- `pandas` is what we'll use to manipulate our data.\n",
"- `matplotlib` is a plotting library.\n",
"- `lifelines` is an open-source survival analysis library."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "0JHzRJaQi_nU"
},
"outputs": [],
"source": [
"import sklearn\n",
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from lifelines import CoxPHFitter\n",
"from lifelines.utils import concordance_index as cindex\n",
"from sklearn.model_selection import train_test_split\n",
"\n",
"from util import load_data"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "vZMwq0VfW5TW"
},
"source": [
"<a name='2'></a>\n",
"## 2. Load the Dataset\n",
"\n",
"Run the next cell to load the data."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"df = load_data()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a name='3'></a>\n",
"## 3. Explore the Dataset\n",
"\n",
"In the lecture videos `time` was in months, however in this assignment, `time` will be converted into years. Also notice that we have assigned a numeric value to `sex`, where `female = 0` and `male = 1`.\n",
"\n",
"Next, familiarize yourself with the data and the shape of it. "
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 241
},
"colab_type": "code",
"id": "T1a_aHGmXT_C",
"outputId": "1bbcf6d9-f293-49f4-963a-827c8e79813b"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(258, 19)\n"
]
},
{
"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>time</th>\n",
" <th>status</th>\n",
" <th>trt</th>\n",
" <th>age</th>\n",
" <th>sex</th>\n",
" <th>ascites</th>\n",
" <th>hepato</th>\n",
" <th>spiders</th>\n",
" <th>edema</th>\n",
" <th>bili</th>\n",
" <th>chol</th>\n",
" <th>albumin</th>\n",
" <th>copper</th>\n",
" <th>alk.phos</th>\n",
" <th>ast</th>\n",
" <th>trig</th>\n",
" <th>platelet</th>\n",
" <th>protime</th>\n",
" <th>stage</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>1.095890</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>58.765229</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>14.5</td>\n",
" <td>261.0</td>\n",
" <td>2.60</td>\n",
" <td>156.0</td>\n",
" <td>1718.0</td>\n",
" <td>137.95</td>\n",
" <td>172.0</td>\n",
" <td>190.0</td>\n",
" <td>12.2</td>\n",
" <td>4.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>12.328767</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>56.446270</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>1.1</td>\n",
" <td>302.0</td>\n",
" <td>4.14</td>\n",
" <td>54.0</td>\n",
" <td>7394.8</td>\n",
" <td>113.52</td>\n",
" <td>88.0</td>\n",
" <td>221.0</td>\n",
" <td>10.6</td>\n",
" <td>3.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>2.772603</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>70.072553</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.5</td>\n",
" <td>1.4</td>\n",
" <td>176.0</td>\n",
" <td>3.48</td>\n",
" <td>210.0</td>\n",
" <td>516.0</td>\n",
" <td>96.10</td>\n",
" <td>55.0</td>\n",
" <td>151.0</td>\n",
" <td>12.0</td>\n",
" <td>4.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>5.273973</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>54.740589</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.5</td>\n",
" <td>1.8</td>\n",
" <td>244.0</td>\n",
" <td>2.54</td>\n",
" <td>64.0</td>\n",
" <td>6121.8</td>\n",
" <td>60.63</td>\n",
" <td>92.0</td>\n",
" <td>183.0</td>\n",
" <td>10.3</td>\n",
" <td>4.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>6</th>\n",
" <td>5.019178</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>55.534565</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>322.0</td>\n",
" <td>4.09</td>\n",
" <td>52.0</td>\n",
" <td>824.0</td>\n",
" <td>60.45</td>\n",
" <td>213.0</td>\n",
" <td>204.0</td>\n",
" <td>9.7</td>\n",
" <td>3.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" time status trt age sex ascites hepato spiders edema \\\n",
"0 1.095890 1.0 0.0 58.765229 0.0 1.0 1.0 1.0 1.0 \n",
"1 12.328767 0.0 0.0 56.446270 0.0 0.0 1.0 1.0 0.0 \n",
"2 2.772603 1.0 0.0 70.072553 1.0 0.0 0.0 0.0 0.5 \n",
"3 5.273973 1.0 0.0 54.740589 0.0 0.0 1.0 1.0 0.5 \n",
"6 5.019178 0.0 1.0 55.534565 0.0 0.0 1.0 0.0 0.0 \n",
"\n",
" bili chol albumin copper alk.phos ast trig platelet protime \\\n",
"0 14.5 261.0 2.60 156.0 1718.0 137.95 172.0 190.0 12.2 \n",
"1 1.1 302.0 4.14 54.0 7394.8 113.52 88.0 221.0 10.6 \n",
"2 1.4 176.0 3.48 210.0 516.0 96.10 55.0 151.0 12.0 \n",
"3 1.8 244.0 2.54 64.0 6121.8 60.63 92.0 183.0 10.3 \n",
"6 1.0 322.0 4.09 52.0 824.0 60.45 213.0 204.0 9.7 \n",
"\n",
" stage \n",
"0 4.0 \n",
"1 3.0 \n",
"2 4.0 \n",
"3 4.0 \n",
"6 3.0 "
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"print(df.shape)\n",
"\n",
"# df.head() only outputs the top few rows\n",
"df.head()"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "Zy5BmjCV-Uo2"
},
"source": [
"Take a minute to examine particular cases."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 357
},
"colab_type": "code",
"id": "01I3ChzL-T-f",
"outputId": "68e209dc-7a44-434b-d44c-4a1e817ee6ca"
},
"outputs": [
{
"data": {
"text/plain": [
"time 11.175342\n",
"status 1.000000\n",
"trt 0.000000\n",
"age 44.520192\n",
"sex 1.000000\n",
"ascites 0.000000\n",
"hepato 1.000000\n",
"spiders 0.000000\n",
"edema 0.000000\n",
"bili 2.100000\n",
"chol 456.000000\n",
"albumin 4.000000\n",
"copper 124.000000\n",
"alk.phos 5719.000000\n",
"ast 221.880000\n",
"trig 230.000000\n",
"platelet 70.000000\n",
"protime 9.900000\n",
"stage 2.000000\n",
"Name: 23, dtype: float64"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"i = 20\n",
"df.iloc[i, :]"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "pYZKl_9Tk2vS"
},
"source": [
"Now, split your dataset into train, validation and test set using 60/20/20 split. "
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "V4HJSZaMk1xG"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Total number of patients: 258\n",
"Total number of patients in training set: 154\n",
"Total number of patients in validation set: 52\n",
"Total number of patients in test set: 52\n"
]
}
],
"source": [
"np.random.seed(0)\n",
"df_dev, df_test = train_test_split(df, test_size = 0.2)\n",
"df_train, df_val = train_test_split(df_dev, test_size = 0.25)\n",
"\n",
"print(\"Total number of patients:\", df.shape[0])\n",
"print(\"Total number of patients in training set:\", df_train.shape[0])\n",
"print(\"Total number of patients in validation set:\", df_val.shape[0])\n",
"print(\"Total number of patients in test set:\", df_test.shape[0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Before proceeding to modeling, let's normalize the continuous covariates to make sure they're on the same scale. Again, we should normalize the test data using statistics from the train data."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"continuous_columns = ['age', 'bili', 'chol', 'albumin', 'copper', 'alk.phos', 'ast', 'trig', 'platelet', 'protime']\n",
"mean = df_train.loc[:, continuous_columns].mean()\n",
"std = df_train.loc[:, continuous_columns].std()\n",
"df_train.loc[:, continuous_columns] = (df_train.loc[:, continuous_columns] - mean) / std\n",
"df_val.loc[:, continuous_columns] = (df_val.loc[:, continuous_columns] - mean) / std\n",
"df_test.loc[:, continuous_columns] = (df_test.loc[:, continuous_columns] - mean) / std"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's check the summary statistics on our training dataset to make sure it's standardized."
]
},
{
"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</th>\n",
" <th>bili</th>\n",
" <th>chol</th>\n",
" <th>albumin</th>\n",
" <th>copper</th>\n",
" <th>alk.phos</th>\n",
" <th>ast</th>\n",
" <th>trig</th>\n",
" <th>platelet</th>\n",
" <th>protime</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>count</th>\n",
" <td>1.540000e+02</td>\n",
" <td>1.540000e+02</td>\n",
" <td>1.540000e+02</td>\n",
" <td>1.540000e+02</td>\n",
" <td>1.540000e+02</td>\n",
" <td>1.540000e+02</td>\n",
" <td>1.540000e+02</td>\n",
" <td>1.540000e+02</td>\n",
" <td>1.540000e+02</td>\n",
" <td>1.540000e+02</td>\n",
" </tr>\n",
" <tr>\n",
" <th>mean</th>\n",
" <td>9.833404e-16</td>\n",
" <td>-3.258577e-16</td>\n",
" <td>1.153478e-16</td>\n",
" <td>1.153478e-16</td>\n",
" <td>5.767392e-18</td>\n",
" <td>1.326500e-16</td>\n",
" <td>-1.263059e-15</td>\n",
" <td>8.074349e-17</td>\n",
" <td>2.018587e-17</td>\n",
" <td>1.291896e-14</td>\n",
" </tr>\n",
" <tr>\n",
" <th>std</th>\n",
" <td>1.000000e+00</td>\n",
" <td>1.000000e+00</td>\n",
" <td>1.000000e+00</td>\n",
" <td>1.000000e+00</td>\n",
" <td>1.000000e+00</td>\n",
" <td>1.000000e+00</td>\n",
" <td>1.000000e+00</td>\n",
" <td>1.000000e+00</td>\n",
" <td>1.000000e+00</td>\n",
" <td>1.000000e+00</td>\n",
" </tr>\n",
" <tr>\n",
" <th>min</th>\n",
" <td>-2.304107e+00</td>\n",
" <td>-5.735172e-01</td>\n",
" <td>-1.115330e+00</td>\n",
" <td>-3.738104e+00</td>\n",
" <td>-9.856552e-01</td>\n",
" <td>-7.882167e-01</td>\n",
" <td>-1.489281e+00</td>\n",
" <td>-1.226674e+00</td>\n",
" <td>-2.058899e+00</td>\n",
" <td>-1.735556e+00</td>\n",
" </tr>\n",
" <tr>\n",
" <th>25%</th>\n",
" <td>-6.535035e-01</td>\n",
" <td>-4.895812e-01</td>\n",
" <td>-5.186963e-01</td>\n",
" <td>-5.697976e-01</td>\n",
" <td>-6.470611e-01</td>\n",
" <td>-5.186471e-01</td>\n",
" <td>-8.353982e-01</td>\n",
" <td>-6.884514e-01</td>\n",
" <td>-6.399831e-01</td>\n",
" <td>-7.382590e-01</td>\n",
" </tr>\n",
" <tr>\n",
" <th>50%</th>\n",
" <td>-6.443852e-03</td>\n",
" <td>-3.846612e-01</td>\n",
" <td>-2.576693e-01</td>\n",
" <td>5.663556e-02</td>\n",
" <td>-3.140636e-01</td>\n",
" <td>-3.416086e-01</td>\n",
" <td>-2.260984e-01</td>\n",
" <td>-2.495932e-01</td>\n",
" <td>-4.100373e-02</td>\n",
" <td>-1.398807e-01</td>\n",
" </tr>\n",
" <tr>\n",
" <th>75%</th>\n",
" <td>5.724289e-01</td>\n",
" <td>2.977275e-02</td>\n",
" <td>1.798617e-01</td>\n",
" <td>6.890921e-01</td>\n",
" <td>3.435366e-01</td>\n",
" <td>-4.620597e-03</td>\n",
" <td>6.061159e-01</td>\n",
" <td>3.755727e-01</td>\n",
" <td>6.617988e-01</td>\n",
" <td>3.587680e-01</td>\n",
" </tr>\n",
" <tr>\n",
" <th>max</th>\n",
" <td>2.654276e+00</td>\n",
" <td>5.239050e+00</td>\n",
" <td>6.243146e+00</td>\n",
" <td>2.140730e+00</td>\n",
" <td>5.495204e+00</td>\n",
" <td>4.869263e+00</td>\n",
" <td>3.058176e+00</td>\n",
" <td>5.165751e+00</td>\n",
" <td>3.190823e+00</td>\n",
" <td>4.447687e+00</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" age bili chol albumin copper \\\n",
"count 1.540000e+02 1.540000e+02 1.540000e+02 1.540000e+02 1.540000e+02 \n",
"mean 9.833404e-16 -3.258577e-16 1.153478e-16 1.153478e-16 5.767392e-18 \n",
"std 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 \n",
"min -2.304107e+00 -5.735172e-01 -1.115330e+00 -3.738104e+00 -9.856552e-01 \n",
"25% -6.535035e-01 -4.895812e-01 -5.186963e-01 -5.697976e-01 -6.470611e-01 \n",
"50% -6.443852e-03 -3.846612e-01 -2.576693e-01 5.663556e-02 -3.140636e-01 \n",
"75% 5.724289e-01 2.977275e-02 1.798617e-01 6.890921e-01 3.435366e-01 \n",
"max 2.654276e+00 5.239050e+00 6.243146e+00 2.140730e+00 5.495204e+00 \n",
"\n",
" alk.phos ast trig platelet protime \n",
"count 1.540000e+02 1.540000e+02 1.540000e+02 1.540000e+02 1.540000e+02 \n",
"mean 1.326500e-16 -1.263059e-15 8.074349e-17 2.018587e-17 1.291896e-14 \n",
"std 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 \n",
"min -7.882167e-01 -1.489281e+00 -1.226674e+00 -2.058899e+00 -1.735556e+00 \n",
"25% -5.186471e-01 -8.353982e-01 -6.884514e-01 -6.399831e-01 -7.382590e-01 \n",
"50% -3.416086e-01 -2.260984e-01 -2.495932e-01 -4.100373e-02 -1.398807e-01 \n",
"75% -4.620597e-03 6.061159e-01 3.755727e-01 6.617988e-01 3.587680e-01 \n",
"max 4.869263e+00 3.058176e+00 5.165751e+00 3.190823e+00 4.447687e+00 "
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"df_train.loc[:, continuous_columns].describe()"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "BX3woHz-jit1"
},
"source": [
"<a name='4'></a>\n",
"## 4. Cox Proportional Hazards\n",
"\n",
"Our goal is to build a risk score using the survival data that we have. We'll begin by fitting a Cox Proportional Hazards model to your data.\n",
"\n",
"Recall that the Cox Proportional Hazards model describes the hazard for an individual $i$ at time $t$ as \n",
"\n",
"$$\n",
"\\lambda(t, x) = \\lambda_0(t)e^{\\theta^T X_i}\n",
"$$\n",
"\n",
"The $\\lambda_0$ term is a baseline hazard and incorporates the risk over time, and the other term incorporates the risk due to the individual's covariates. After fitting the model, we can rank individuals using the person-dependent risk term $e^{\\theta^T X_i}$. \n",
"\n",
"Categorical variables cannot be used in a regression model as they are. In order to use them, conversion to a series of variables is required.\n",
"\n",
"Since our data has a mix of categorical (`stage`) and continuous (`wblc`) variables, before we proceed further we need to do some data engineering. To tackle the issue at hand we'll be using the `Dummy Coding` technique. In order to use Cox Proportional Hazards, we will have to turn the categorical data into one hot features so that we can fit our Cox model. Luckily, Pandas has a built-in function called `get_dummies` that will make it easier for us to implement our function. It turns categorical features into multiple binary features.\n",
"\n",
"<img src=\"1-hot-encode.png\" style=\"padding-top: 5px;width: 60%;left: 0px;margin-left: 150px;margin-right: 0px;\">\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a name='Ex-1'></a>\n",
"### Exercise 1\n",
"In the cell below, implement the `to_one_hot(...)` function."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<details> \n",
"<summary>\n",
" <font size=\"3\" color=\"darkgreen\"><b>Hints</b></font>\n",
"</summary>\n",
"<p>\n",
"<ul>\n",
" <li>Remember to drop the first dummy for each each category to avoid convergence issues when fitting the proportional hazards model.</li>\n",
" <li> Check out the <a href=\"https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.get_dummies.html\" > get_dummies() </a> documentation. </li>\n",
" <li>Use <code>dtype=np.float64</code>.</li>\n",
"</ul>\n",
"</p>"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "VMzvx0xF_C3I"
},
"outputs": [],
"source": [
"# UNQ_C1 (UNIQUE CELL IDENTIFIER, DO NOT EDIT)\n",
"def to_one_hot(dataframe, columns):\n",
" '''\n",
" Convert columns in dataframe to one-hot encoding.\n",
" Args:\n",
" dataframe (dataframe): pandas dataframe containing covariates\n",
" columns (list of strings): list categorical column names to one hot encode\n",
" Returns:\n",
" one_hot_df (dataframe): dataframe with categorical columns encoded\n",
" as binary variables\n",
" '''\n",
" \n",
" ### START CODE HERE (REPLACE INSTANCES OF 'None' with your code) ###\n",
" \n",
" one_hot_df = pd.get_dummies(dataframe, columns = columns, drop_first = True, dtype=np.float64)\n",
" \n",
" ### END CODE HERE ###\n",
" \n",
" return one_hot_df"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "rM2tIzvG_ifc"
},
"source": [
"Now we'll use the function you coded to transform the training, validation, and test sets."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "SGZfLeup_fUL"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"['time', 'status', 'trt', 'age', 'sex', 'ascites', 'hepato', 'spiders', 'bili', 'chol', 'albumin', 'copper', 'alk.phos', 'ast', 'trig', 'platelet', 'protime', 'edema_0.5', 'edema_1.0', 'stage_2.0', 'stage_3.0', 'stage_4.0']\n",
"There are 22 columns\n"
]
}
],
"source": [
"# List of categorical columns\n",
"to_encode = ['edema', 'stage']\n",
"\n",
"one_hot_train = to_one_hot(df_train, to_encode)\n",
"one_hot_val = to_one_hot(df_val, to_encode)\n",
"one_hot_test = to_one_hot(df_test, to_encode)\n",
"\n",
"print(one_hot_val.columns.tolist())\n",
"print(f\"There are {len(one_hot_val.columns)} columns\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Expected output\n",
"```Python\n",
"['time', 'status', 'trt', 'age', 'sex', 'ascites', 'hepato', 'spiders', 'bili', 'chol', 'albumin', 'copper', 'alk.phos', 'ast', 'trig', 'platelet', 'protime', 'edema_0.5', 'edema_1.0', 'stage_2.0', 'stage_3.0', 'stage_4.0']\n",
"There are 22 columns\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Look for new features\n",
"Now, let's take a peek at one of the transformed data sets. Do you notice any new features?"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 241
},
"colab_type": "code",
"id": "w8EG8A9gXcpu",
"outputId": "384d9ade-2c96-4979-d3b7-da2b8e50f2e0"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(154, 22)\n"
]
},
{
"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>time</th>\n",
" <th>status</th>\n",
" <th>trt</th>\n",
" <th>age</th>\n",
" <th>sex</th>\n",
" <th>ascites</th>\n",
" <th>hepato</th>\n",
" <th>spiders</th>\n",
" <th>bili</th>\n",
" <th>chol</th>\n",
" <th>...</th>\n",
" <th>alk.phos</th>\n",
" <th>ast</th>\n",
" <th>trig</th>\n",
" <th>platelet</th>\n",
" <th>protime</th>\n",
" <th>edema_0.5</th>\n",
" <th>edema_1.0</th>\n",
" <th>stage_2.0</th>\n",
" <th>stage_3.0</th>\n",
" <th>stage_4.0</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>279</th>\n",
" <td>3.868493</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>-0.414654</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>-0.300725</td>\n",
" <td>-0.096081</td>\n",
" <td>...</td>\n",
" <td>0.167937</td>\n",
" <td>0.401418</td>\n",
" <td>0.330031</td>\n",
" <td>0.219885</td>\n",
" <td>-1.137178</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>137</th>\n",
" <td>3.553425</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>0.069681</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>0.895363</td>\n",
" <td>0.406085</td>\n",
" <td>...</td>\n",
" <td>0.101665</td>\n",
" <td>0.472367</td>\n",
" <td>1.621764</td>\n",
" <td>-0.120868</td>\n",
" <td>-0.239610</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>249</th>\n",
" <td>4.846575</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>-0.924494</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>-0.510565</td>\n",
" <td>-0.225352</td>\n",
" <td>...</td>\n",
" <td>0.245463</td>\n",
" <td>1.899020</td>\n",
" <td>-0.580807</td>\n",
" <td>0.422207</td>\n",
" <td>0.159309</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>266</th>\n",
" <td>0.490411</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>1.938314</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>0.748475</td>\n",
" <td>-0.608191</td>\n",
" <td>...</td>\n",
" <td>-0.650254</td>\n",
" <td>-0.288898</td>\n",
" <td>-0.481443</td>\n",
" <td>-0.727833</td>\n",
" <td>1.356065</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>12.328767</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.563645</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>1.0</td>\n",
" <td>-0.405645</td>\n",
" <td>-0.210436</td>\n",
" <td>...</td>\n",
" <td>2.173526</td>\n",
" <td>-0.144699</td>\n",
" <td>-0.531125</td>\n",
" <td>-0.450972</td>\n",
" <td>-0.139881</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>0.0</td>\n",
" <td>1.0</td>\n",
" <td>0.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>5 rows × 22 columns</p>\n",
"</div>"
],
"text/plain": [
" time status trt age sex ascites hepato spiders \\\n",
"279 3.868493 0.0 0.0 -0.414654 0.0 0.0 0.0 0.0 \n",
"137 3.553425 1.0 0.0 0.069681 1.0 0.0 1.0 0.0 \n",
"249 4.846575 0.0 1.0 -0.924494 0.0 0.0 1.0 0.0 \n",
"266 0.490411 1.0 0.0 1.938314 0.0 1.0 1.0 1.0 \n",
"1 12.328767 0.0 0.0 0.563645 0.0 0.0 1.0 1.0 \n",
"\n",
" bili chol ... alk.phos ast trig platelet \\\n",
"279 -0.300725 -0.096081 ... 0.167937 0.401418 0.330031 0.219885 \n",
"137 0.895363 0.406085 ... 0.101665 0.472367 1.621764 -0.120868 \n",
"249 -0.510565 -0.225352 ... 0.245463 1.899020 -0.580807 0.422207 \n",
"266 0.748475 -0.608191 ... -0.650254 -0.288898 -0.481443 -0.727833 \n",
"1 -0.405645 -0.210436 ... 2.173526 -0.144699 -0.531125 -0.450972 \n",
"\n",
" protime edema_0.5 edema_1.0 stage_2.0 stage_3.0 stage_4.0 \n",
"279 -1.137178 0.0 0.0 0.0 1.0 0.0 \n",
"137 -0.239610 0.0 0.0 0.0 1.0 0.0 \n",
"249 0.159309 0.0 0.0 0.0 0.0 1.0 \n",
"266 1.356065 0.0 1.0 0.0 0.0 1.0 \n",
"1 -0.139881 0.0 0.0 0.0 1.0 0.0 \n",
"\n",
"[5 rows x 22 columns]"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"print(one_hot_train.shape)\n",
"one_hot_train.head()"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "hNxuymLwyjqM"
},
"source": [
"<a name='5'></a>\n",
"## 5. Fitting and Interpreting a Cox Model"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "ygiFcUKcAFQk"
},
"source": [
"Run the following cell to fit your Cox Proportional Hazards model using the `lifelines` package."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
},
"colab_type": "code",
"id": "dDCS7p3xjbXB",
"outputId": "41b12f82-8b35-43e1-d2a9-05258ac50b20"
},
"outputs": [
{
"data": {
"text/plain": [
"<lifelines.CoxPHFitter: fitted with 154 total observations, 90 right-censored observations>"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"cph = CoxPHFitter()\n",
"cph.fit(one_hot_train, duration_col = 'time', event_col = 'status', step_size=0.1)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "5MUITR0QANDH"
},
"source": [
"You can use `cph.print_summary()` to view the coefficients associated with each covariate as well as confidence intervals. "
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 903
},
"colab_type": "code",
"id": "fH5AZs8vjcEv",
"outputId": "5429f7d5-5669-431f-a014-cf609c90997f"
},
"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",
" <tbody>\n",
" <tr>\n",
" <th>model</th>\n",
" <td>lifelines.CoxPHFitter</td>\n",
" </tr>\n",
" <tr>\n",
" <th>duration col</th>\n",
" <td>'time'</td>\n",
" </tr>\n",
" <tr>\n",
" <th>event col</th>\n",
" <td>'status'</td>\n",
" </tr>\n",
" <tr>\n",
" <th>number of observations</th>\n",
" <td>154</td>\n",
" </tr>\n",
" <tr>\n",
" <th>number of events observed</th>\n",
" <td>64</td>\n",
" </tr>\n",
" <tr>\n",
" <th>partial log-likelihood</th>\n",
" <td>-230.82</td>\n",
" </tr>\n",
" <tr>\n",
" <th>time fit was run</th>\n",
" <td>2020-06-26 05:28:02 UTC</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div><table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>coef</th>\n",
" <th>exp(coef)</th>\n",
" <th>se(coef)</th>\n",
" <th>coef lower 95%</th>\n",
" <th>coef upper 95%</th>\n",
" <th>exp(coef) lower 95%</th>\n",
" <th>exp(coef) upper 95%</th>\n",
" <th>z</th>\n",
" <th>p</th>\n",
" <th>-log2(p)</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>trt</th>\n",
" <td>-0.22</td>\n",
" <td>0.80</td>\n",
" <td>0.30</td>\n",
" <td>-0.82</td>\n",
" <td>0.37</td>\n",
" <td>0.44</td>\n",
" <td>1.45</td>\n",
" <td>-0.73</td>\n",
" <td>0.46</td>\n",
" <td>1.11</td>\n",
" </tr>\n",
" <tr>\n",
" <th>age</th>\n",
" <td>0.23</td>\n",
" <td>1.26</td>\n",
" <td>0.19</td>\n",
" <td>-0.13</td>\n",
" <td>0.60</td>\n",
" <td>0.88</td>\n",
" <td>1.82</td>\n",
" <td>1.26</td>\n",
" <td>0.21</td>\n",
" <td>2.27</td>\n",
" </tr>\n",
" <tr>\n",
" <th>sex</th>\n",
" <td>0.34</td>\n",
" <td>1.41</td>\n",
" <td>0.40</td>\n",
" <td>-0.45</td>\n",
" <td>1.14</td>\n",
" <td>0.64</td>\n",
" <td>3.11</td>\n",
" <td>0.84</td>\n",
" <td>0.40</td>\n",
" <td>1.33</td>\n",
" </tr>\n",
" <tr>\n",
" <th>ascites</th>\n",
" <td>-0.10</td>\n",
" <td>0.91</td>\n",
" <td>0.56</td>\n",
" <td>-1.20</td>\n",
" <td>1.01</td>\n",
" <td>0.30</td>\n",
" <td>2.75</td>\n",
" <td>-0.17</td>\n",
" <td>0.86</td>\n",
" <td>0.21</td>\n",
" </tr>\n",
" <tr>\n",
" <th>hepato</th>\n",
" <td>0.31</td>\n",
" <td>1.36</td>\n",
" <td>0.38</td>\n",
" <td>-0.44</td>\n",
" <td>1.06</td>\n",
" <td>0.64</td>\n",
" <td>2.89</td>\n",
" <td>0.81</td>\n",
" <td>0.42</td>\n",
" <td>1.26</td>\n",
" </tr>\n",
" <tr>\n",
" <th>spiders</th>\n",
" <td>-0.18</td>\n",
" <td>0.83</td>\n",
" <td>0.38</td>\n",
" <td>-0.94</td>\n",
" <td>0.57</td>\n",
" <td>0.39</td>\n",
" <td>1.77</td>\n",
" <td>-0.47</td>\n",
" <td>0.64</td>\n",
" <td>0.66</td>\n",
" </tr>\n",
" <tr>\n",
" <th>bili</th>\n",
" <td>0.05</td>\n",
" <td>1.05</td>\n",
" <td>0.18</td>\n",
" <td>-0.29</td>\n",
" <td>0.39</td>\n",
" <td>0.75</td>\n",
" <td>1.48</td>\n",
" <td>0.29</td>\n",
" <td>0.77</td>\n",
" <td>0.37</td>\n",
" </tr>\n",
" <tr>\n",
" <th>chol</th>\n",
" <td>0.19</td>\n",
" <td>1.20</td>\n",
" <td>0.15</td>\n",
" <td>-0.10</td>\n",
" <td>0.47</td>\n",
" <td>0.91</td>\n",
" <td>1.60</td>\n",
" <td>1.28</td>\n",
" <td>0.20</td>\n",
" <td>2.33</td>\n",
" </tr>\n",
" <tr>\n",
" <th>albumin</th>\n",
" <td>-0.40</td>\n",
" <td>0.67</td>\n",
" <td>0.18</td>\n",
" <td>-0.75</td>\n",
" <td>-0.06</td>\n",
" <td>0.47</td>\n",
" <td>0.94</td>\n",
" <td>-2.28</td>\n",
" <td>0.02</td>\n",
" <td>5.46</td>\n",
" </tr>\n",
" <tr>\n",
" <th>copper</th>\n",
" <td>0.30</td>\n",
" <td>1.35</td>\n",
" <td>0.16</td>\n",
" <td>-0.01</td>\n",
" <td>0.61</td>\n",
" <td>0.99</td>\n",
" <td>1.84</td>\n",
" <td>1.91</td>\n",
" <td>0.06</td>\n",
" <td>4.14</td>\n",
" </tr>\n",
" <tr>\n",
" <th>alk.phos</th>\n",
" <td>-0.22</td>\n",
" <td>0.80</td>\n",
" <td>0.14</td>\n",
" <td>-0.49</td>\n",
" <td>0.05</td>\n",
" <td>0.61</td>\n",
" <td>1.05</td>\n",
" <td>-1.62</td>\n",
" <td>0.11</td>\n",
" <td>3.24</td>\n",
" </tr>\n",
" <tr>\n",
" <th>ast</th>\n",
" <td>0.21</td>\n",
" <td>1.24</td>\n",
" <td>0.16</td>\n",
" <td>-0.10</td>\n",
" <td>0.53</td>\n",
" <td>0.91</td>\n",
" <td>1.69</td>\n",
" <td>1.34</td>\n",
" <td>0.18</td>\n",
" <td>2.48</td>\n",
" </tr>\n",
" <tr>\n",
" <th>trig</th>\n",
" <td>0.20</td>\n",
" <td>1.23</td>\n",
" <td>0.16</td>\n",
" <td>-0.11</td>\n",
" <td>0.52</td>\n",
" <td>0.89</td>\n",
" <td>1.68</td>\n",
" <td>1.27</td>\n",
" <td>0.21</td>\n",
" <td>2.28</td>\n",
" </tr>\n",
" <tr>\n",
" <th>platelet</th>\n",
" <td>0.14</td>\n",
" <td>1.15</td>\n",
" <td>0.15</td>\n",
" <td>-0.16</td>\n",
" <td>0.43</td>\n",
" <td>0.86</td>\n",
" <td>1.54</td>\n",
" <td>0.92</td>\n",
" <td>0.36</td>\n",
" <td>1.48</td>\n",
" </tr>\n",
" <tr>\n",
" <th>protime</th>\n",
" <td>0.36</td>\n",
" <td>1.43</td>\n",
" <td>0.17</td>\n",
" <td>0.03</td>\n",
" <td>0.69</td>\n",
" <td>1.03</td>\n",
" <td>1.99</td>\n",
" <td>2.15</td>\n",
" <td>0.03</td>\n",
" <td>4.97</td>\n",
" </tr>\n",
" <tr>\n",
" <th>edema_0.5</th>\n",
" <td>1.24</td>\n",
" <td>3.47</td>\n",
" <td>0.46</td>\n",
" <td>0.35</td>\n",
" <td>2.14</td>\n",
" <td>1.42</td>\n",
" <td>8.50</td>\n",
" <td>2.72</td>\n",
" <td>0.01</td>\n",
" <td>7.28</td>\n",
" </tr>\n",
" <tr>\n",
" <th>edema_1.0</th>\n",
" <td>2.02</td>\n",
" <td>7.51</td>\n",
" <td>0.60</td>\n",
" <td>0.84</td>\n",
" <td>3.20</td>\n",
" <td>2.31</td>\n",
" <td>24.43</td>\n",
" <td>3.35</td>\n",
" <td><0.005</td>\n",
" <td>10.28</td>\n",
" </tr>\n",
" <tr>\n",
" <th>stage_2.0</th>\n",
" <td>1.21</td>\n",
" <td>3.35</td>\n",
" <td>1.08</td>\n",
" <td>-0.92</td>\n",
" <td>3.33</td>\n",
" <td>0.40</td>\n",
" <td>28.06</td>\n",
" <td>1.11</td>\n",
" <td>0.27</td>\n",
" <td>1.91</td>\n",
" </tr>\n",
" <tr>\n",
" <th>stage_3.0</th>\n",
" <td>1.18</td>\n",
" <td>3.27</td>\n",
" <td>1.09</td>\n",
" <td>-0.96</td>\n",
" <td>3.33</td>\n",
" <td>0.38</td>\n",
" <td>27.86</td>\n",
" <td>1.08</td>\n",
" <td>0.28</td>\n",
" <td>1.84</td>\n",
" </tr>\n",
" <tr>\n",
" <th>stage_4.0</th>\n",
" <td>1.41</td>\n",
" <td>4.10</td>\n",
" <td>1.15</td>\n",
" <td>-0.85</td>\n",
" <td>3.67</td>\n",
" <td>0.43</td>\n",
" <td>39.43</td>\n",
" <td>1.22</td>\n",
" <td>0.22</td>\n",
" <td>2.18</td>\n",
" </tr>\n",
" </tbody>\n",
"</table><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",
" <tbody>\n",
" <tr>\n",
" <th>Concordance</th>\n",
" <td>0.83</td>\n",
" </tr>\n",
" <tr>\n",
" <th>Log-likelihood ratio test</th>\n",
" <td>97.63 on 20 df, -log2(p)=38.13</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"cph.print_summary()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Question:**\n",
"\n",
"- According to the model, was treatment `trt` beneficial? \n",
"- What was its associated hazard ratio? \n",
" - Note that the hazard ratio is how much an incremental increase in the feature variable changes the hazard."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<details> \n",
"<summary>\n",
" <font size=\"3\" color=\"darkgreen\"><b>Check your answer!</b></font>\n",
"</summary>\n",
"<p>\n",
"<ul>\n",
"<ul>\n",
" <li>You should see that the treatment (trt) was beneficial because it has a negative impact on the hazard (the coefficient is negative, and exp(coef) is less than 1).</li>\n",
" <li>The associated hazard ratio is ~0.8, because this is the exp(coef) of treatment.</li>\n",
"</ul>\n",
"</p>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can compare the predicted survival curves for treatment variables. Run the next cell to plot survival curves using the `plot_covariate_groups()` function. \n",
"- The y-axis is th survival rate\n",
"- The x-axis is time"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 282
},
"colab_type": "code",
"id": "Uxl0icyBS4Dr",
"outputId": "5fa08369-e89e-424f-f9f0-60cf7a1cfbcd"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3df1yUZb7/8deV2qrlr6TSFfyVRpp6sDDWbHMq3DQNKy2IbBPb6FuaWdlJ9+SPrE062jnqZj+0Fd08rrZqpVZqbqKZZlBLK6KWaSmupmhqairG9f1jYBxggAFmGGZ4Px8PHnHPfXPPZ9I+XVz35/pcxlqLiIgEvwsCHYCIiPiGErqISIhQQhcRCRFK6CIiIUIJXUQkRNQN1BuHhYXZtm3bBurtRUSC0hdffJFrrb3U07mAJfS2bduSkZERqLcXEQlKxpjvSzunKRcRkRChhC4iEiKU0EVEQkTA5tBFpGry8vLIycnh9OnTgQ5F/KB+/fqEh4dTr149r39GCV0kSOXk5NCoUSPatm2LMSbQ4YgPWWs5fPgwOTk5tGvXzuufK3fKxRgzxxhz0BiTVcp5Y4yZYYzZaYz5lzHmmgrELSKVdPr0aZo3b65kHoKMMTRv3rzCv315M4c+F+hbxvl+QMeCr2TgtQpFICKVpmQeuirzZ1tuQrfWrgeOlHHJQOCv1ukzoKkxpmV599225Z+Mffhu7yMVEZEy+aLKpRWw1+04p+C1EowxycaYDGNMxi+/5BP5wwoyn7+eHu2bkHLPlWz++8s+CEdEqsvRo0d59dVXK3yuNLt37yYmJoYOHToQHx/P2bNnfRFmrVGtZYvW2lnW2mhrbXS7NuH0uO5a17nwc3vpmvkcPdo30chdJEiUlrTPnTtXqYT+zDPP8MQTT7Bz506aNWvGX/7yF1+FWiv4osplHxDhdhxe8FqZLmp2OVf/cQMA6eNg899fJnvrYs79/DlRh9/n2MzfMfDPXzAyti13XduC3J/OsvTAr0l+cZ4PQhYRXxgzZgzffvstUVFR1KtXj/r169OsWTO2b9/ONddc4zrXp08fpkyZUua9rLV8/PHHLFiwAIAHHniAiRMn8sgjj1THRwkJvkjoy4ARxpiFQAxwzFq7v6I3ibn7Kbj7KR5pMo72x9aw7YfjnDx7jj0/nmLr/mOkrkhnU84vJHfMha6DITrJB6GLhIbnlm8l+9/HfXrPzr9uzITbry7zmpSUFLKyssjMzCQtLY3+/fuTlZVFu3bt+O6771znAH766Sd++9vferzPggULuOyyy2jatCl16zrTUnh4OPv2lTs2FDflJnRjzN8ABxBmjMkBJgD1AKy1rwMfALcBO4FTQJUy7cXXP8TUzNsAaP8MbML5dbzpOO5ovZcDW9fD1vW02LJYiV2khrnuuutKrZtu1KiRK7l7kpub66+wao1yE7q19t5yzltguK8CSoxpTWJM6xKvL4iaw3uZ+/jzi/cQZo7x5cUb2Puv9dw/ZDTPDuhA7NVhJW+mhC+1RHkj6epy0UUXlXquvBF6p06dOHr0KOfOnaNu3brk5OTQqpXH+gopRdCsFC1M9H/4bjhfn/wVY/P28pv8tZw8u4vvjpzk8911+M+3d/DcHR3pHdmcHV+sY93Sj0m+a/H5myjBi/hUo0aN+Omnn7w6V94IHeCmm25i8eLFJCQkMG/ePAYOHOjTeENd0CT0Qm9OHgPAgs17WJB5N+07wwu7j5B38AB5+TOZmPd7Ljvbnb2f59Do53+TDGTl/MSIv2wksatbgldyF6my5s2b06tXL7p06UKDBg24/PLLPZ7r169fuQ9FAV566SUSEhJ49tln6d69Ow8++KA/ww85xjljUv2io6Otrza4WLB5D+9lFn14suadBdQLa0Pv3/bi2L5vOfjun3ixz8W0aFyfUambSOxaj+S7bj7/A0rwEmS2bdtGp06dAh2G+JGnP2NjzBfW2mhP1wfdCN0TT/PuC6JauZJ8k1ZXsP22qfw/4Mq8XE42O8IVnS4iffdRnn57O6/cfIYu3ztLKJXURSRYhURC96R4kj8/ir+Er2/9bx4EOuYdIDt/Jn+/uAfH984n67lHNecuIkErZBN6ce4J3j253/TUTObtPsKCT9ZwbbNTJAPrdhxmwsKMonPuxSnZi0gNU2sSujtPo/fHMm/jk4guxJ/txsG8f/JLy7NEda/Pmq25vLBiJ2899B9EXNKAlVsOkbL0Sz1gFZEap1Ym9OISY1rDc0+75twvu7I72Y0m83LLxjx0xUnIeAHumQcREbByJVH7XyYuxrI88wf279pK8vcbYIumakQksJTQCxQftce/sYns/ceZTWMuv3cyoz/IwdlIsgk4JvFxVCtm/v1eerbrDW3qnL/R9xucX6CkLiLVSptEl2JgVCs6t2zs8dzm3Uf44ztbONPuBjren8LSJg/imHeSY3ctgAHTnBetGAWp/Yt+ZaRW4ycQ8T9ft8995ZVX6NChA8YYtQKoBI3QS1FaCwI4/1B1Mzcyee0+7mx6yHVu0bcNSflbUx75TWOS27j9kEbuEoIKk/ajjz5a5HX39rnFz5WlV69eDBgwAIfD4eNIa4eQWFgUKAs27+GP72wBIKbdJQB8+8m7tP7Vaeb/eTJLly7lzTff5KOPPqLeV/Odo3aANjecv4nm26WSasLCooSEBN577z0iIyM9ts8tPOdN+1x3bdu2JSMjg7AwDz2aapFaubAoUApH8O6rVHPDbyQXeHdL0V8X52YZUuY2ZPsEtz20C0ftWzyURirRS0V8OAYObPHtPVt0hX4pZV7iy/a5nTt39m38tZASehV5KoH84ztb+OM7W4hpdwWX3zuZIXMy+PaTbVxwSVtIep/Zs2ezaNEi1qRM85zMleglSFWlfa5UnRK6j3katYNz5E74jSzYvOf8i9FJvPr5zyxfvpwPP/wQgOnTp/OPj//NsqRfl7y55uGlNOWMpKtLVdrnaoRedUrofuCxt0yRkXsXmt/dhfg3NtFw95GSN2gaAUnLmDp1Kps2bWLJkiUApIwYzCXfv09ym8VK6FIj+Lp9rlSNyharSWJMa168s6vr4SlA9v7jnOpwi2t0DvD444+zbNkyj/eYuyaLjw82c86VqhRSagD3FrlPP/201+dKM2PGDMLDw8nJyaFbt2784Q9/8EfYIUtVLgFUuHjJU737wKhWnssmM1LPz60f2OJ8cJX0vp8jlZqoJlS5iH+pyiWIDIzyvL3W5t1H2Lz7SIl5eID8L7PoFv4bxo0b5xyhF47WPdEDVJFaRQk9gErdP9XDhh3gTPS5G75kf0RTxgHjP61DeF4Tklt4uHlZlTLFKfGLhAQl9BqorET/R0bTqmAe/u0NX9OzZ0+Sk1IZO3Yshw8fZtasWQCMHtKXK+vsK7pa1ZPC2mUldJGgp4QeRBJjWhcZuW/fvr3Ua1dkfMctt9wCSTPLvqmnaRuN2EWCkhJ6EMref5z4NzYB5x+eTp48ucg1hcl+1Chnu4Fp05xNw4YPH06DBg2YOnUqAMnLThDd2G3aRiN2kaClhB5k3B+kZu8/DlBqEzFvrN/2A83vvBOSCv6HoBG7SNBS2WIQK6vsEcooffQgKSmJiIgIJsW1K/ogVaWRNVagyxa/++47BgwYQFZWll/un5aWxtSpU1mxYgXLli0jOzubMWPG+OW9qiojI4O//vWvzJgxo8I/O3fuXDIyMnjllVdKnFPZYi1SWtkjlF36WOZ9opOKjsaLj9g1WpcAiIuLIy4uLqAxnDt3jrp1PafM6OhooqM95thqpZWiQSwxpjWLHu7p8av4qtTSZO8/znuZ+0hNTWXSpEkMGTKEq666ynU+4W+5zNrexHlwYIt3ZZBSa5w7d4777ruPTp06MXjwYE6dOgXApEmT6NGjB126dCE5OZnCmYAZM2bQuXNnunXrRkJCAgAnT55k2LBhXHfddXTv3p333nuvxPvMnTuXESNGADB06FBGjhzJ9ddfT/v27Vm8+PzfySlTptCjRw+6devGhAkTStznl19+YejQoXTp0oWuXbvyv//7vwA4HA4KZwxyc3Np27at633j4uK4+eabueWWW0hISOD998//tjp06FAWL15MWloaAwYMID8/n7Zt23L06FHXNR07duSHH35g+fLlxMTE0L17d2JjY/nhhx8q/e+9VNbagHxde+21VgLvntc32i4TVtp7Xt9o73l9o732vv+0g5OfdJ2PjIy0kydPdh7Muc3aFyOc/yz+lT4nQJ+g9srOzi5y3Lt3b5uammqttfbs2bO2d+/e9q233rLWWnvy5Enbu3dvu3DhQmuttUePHrW9e/e2S5YssdZae+jQIdu7d2+7bNkya621+/fvL/f9d+/ebQG7YcMGa621SUlJdsqUKdZaaw8fPuy6bsiQIa77tmzZ0p4+fdpaa+2PP/5orbV27Nixrjh//PFH27FjR3vixAm7du1a279/f2uttampqXb48OHWWmsfeOABO3jwYPvLL7/YrVu32iuuuMJaa+2qVavsQw89ZPPz8+0vv/xi+/fvb9etW1ck5oyMDBsbG+s6Loyhd+/eNj093fXvok2bNq73bdWqlevzLF261P7+97+31lp75swZGx4ebk+dOlUk1pEjR9o5c5z/PXz22Wf2lltusdZae+TIEZufn2+ttXb27Nn2ySefLPHZiiv+Z2yttUCGLSWvaoReyxXfau9Mh5u54JrBruPt27czZswYBg0axNTMi53z6cVp5F5rRURE0KtXLwCGDBnChg3ObqBr164lJiaGrl278vHHH7N161YAunXrxn333cf8+fNd0xerV68mJSWFqKgoHA4Hp0+fZs+ePZ7fsMAdd9zBBRdcQOfOnV0j3dWrV7N69Wq6d+/ONddcw/bt2/nmm2+K/Fz79u3ZtWsXjz32GCtXrqRxY8/Pn9z16dOHSy5x/rbbr18/1q5dy5kzZ/jwww+58cYbadCgQZHr4+PjWbRoEQALFy4kPj4egJycHG699Va6du3KlClTXP9OfElz6LWcp82xSxXRA5JGl3xd8+w1Qlpamuv7evXqFTlu2LBhkeMmTZoUOQ4LCyty3KKFp+XHJRljShyfPn2aRx99lIyMDCIiIpg4cSKnT58G4P3332f9+vUsX76cP/3pT2zZsgVrLUuWLCEyMrLIvcqakvjVr37l+t4WTOdYaxk7diwPP/xwqT/XrFkzvvrqK1atWsXrr7/O22+/zZw5c6hbty75+fkArlgLubcErl+/Pg6Hg1WrVrFo0SLXtJG7nj17snPnTg4dOsS7777Ls88+C8Bjjz3Gk08+SVxcHGlpaUycOLHUOCtLI3QpobDOvfBrweY9LFmyhNGjRxMXF8f06dNd1/br149Xt19yfuSu0XqtsmfPHjZtcg4CFixYwA033OBKiGFhYZw4ccI1x52fn8/evXu56aabeOmllzh27BgnTpzg1ltv5c9//rMrMf/zn/+sVCy33norc+bM4cSJEwDs27ePgwcPFrkmNzeX/Px8Bg0axAsvvMCXX34JOLe8++KLLwCKzMl7Eh8fT2pqKp988gl9+/Ytcd4Yw5133smTTz5Jp06daN68OQDHjh2jVStnAcK8efMq9RnLoxG6FFG8csarWvc2PSGpYCNgNQyrVSIjI5k5cybDhg2jc+fOPPLIIzRs2JCHHnqILl260KJFC3r06AE4H0gOGTKEY8eOYa1l5MiRNG3alHHjxjFq1Ci6detGfn4+7dq1Y8WKFRWO5Xe/+x3btm2jZ8+eAFx88cXMnz+fyy67zHXNvn37SEpKco3GCxfkjR49mnvuuYdZs2bRv38pf3fd3uf+++9n4MCBXHjhhR6viY+Pp0ePHsydO9f12sSJE7n77rtp1qwZN998M7t3767wZyyP6tClTMVr3cutbXdv71ucatp9KtB16OJ/qkMXn/JmZWpsbCw5OTnOdgPF69jdlTZqFxGfUEKXMrk/NC0crRd/cHq0xbXc0KMP4KznhaIP6IrwNB2jaRgRn9BDUfFa8RLHQmc63MxPbW8EnAsthg4dSl5eHg6Hg/nz5wNw6tQpHK98e36RUiE9RK2SQE2Ziv9V5s9WI3TxWml92t1H7EOHDgUgLy+vxHXH8hvQ5LYxUFCXC5T9EFUj9zLVr1+fw4cP07x58xLlgxLcrLUcPnyY+vXrV+jnlNDFJzy19C1eB+2xHK3r4JKvgdr4eqFwM+VDhw4FOhTxg/r16xMeHl6hn1FClyrztqXvsWPHSElJoUePHtx1113k5uYyePQ8EhMTSU5OLnqxHqCWq169erRr1y7QYUgN4lVCN8b0BaYDdYA3rbUpxc63BuYBTQuuGWOt/cDHsUoN5c2DU4BvP3mXU9nrXXXJ4Jyaadmypecbl1XP7k5TMyKAFwndGFMHmAn0AXKAdGPMMmtttttlzwJvW2tfM8Z0Bj4A2vohXqnhymrpe6bDzXT97R3cdZdz4UdYWBiffvqp54tLm4opTlMzIi7ejNCvA3Zaa3cBGGMWAgMB94RugcLyhybAv30ZpASP0h6cguc+MQcOHCAlJYW+ffvSt29f9u7dy/3338+zzz5LrDcLkDQ1I+LiTUJvBex1O84BYopdMxFYbYx5DLgIiPV0I2NMMpAM0Lp15bdNk+BVfDrGOQ2z2WNPDK8VTs1o6kVqOV/Vod8LzLXWhgO3AW8ZY0rc21o7y1obba2NvvTSS3301hIsPNWxn+lwM10f/h9XQo+IiCAtLc3V5rRcXQc72wmonl3Eq4S+D4hwOw4veM3dg8DbANbaTUB9IMwXAUro8LTDUmn7oQ4bNozhw4cDsGPHDhwOBxs3bgQgKysLh8NBenq6c0Se9L7nPu0itYw3Uy7pQEdjTDuciTwBSCx2zR7gFmCuMaYTzoSu4liptMTExBL9sT3JzMwEIMrfAYkEgXITurX2nDFmBLAKZ0niHGvtVmPMJJxbIS0DngJmG2OewPmAdKjVmmTxkucyx65cXN9ZMRMZGVlkkVKXLl1cx67eMQ9c5H2ZI2i+XUKSV3XoBTXlHxR7bbzb99lAL9+GJrVBaWWOXvVhB6ZNm+b85lwFNkVQqaOEKK0UlYDypj9MWaKiokhPT8fxdCkrTj1RqaOEKCV0qbFKW3FaXJcLcwkPD+f666/3/uYqdZQQpIQuNVJZK07dZe8/Di3DWFTQptcrhatQNfUiIUZb0ElQKxzBP97VkpKSwpQpU4iMjGTdunVMmDCBOXPm0L59e9asWcMLL7zAW2+9RUREQRVu4dSLtsSTIKIt6CTkZWVlcfz4ca+uXblyJQB9oWhljKZfJMhphC5BzX0T63I3sC7gKnWc+sD51aXawFqChEboErIK59q9LXMEWLhwofObFi3Oj8hV+SIhQHuKSlArbCdQWgsBT1q0aOEsdXQ4mDVrlh+jE6leGqFLrRUdHc1dd911/oWyVppqfl2CgBK6hAxv69YBBkb9B1On3n7+hbI21FB5owQJJXQJCd7WrcP5+fb6+zJISUnho48+okl0UukJW/PrEiSU0CUklLVTUnGFo/jc3FwaNmwIwKJFi5zn4uM9/9D3GyAjVaN0qdH0UFRqreTkZNLS0mjSpAmvvfYar732mucLC6djtIGG1HAaoUutVDjfXli7/sEHH5R+cXSSkrkEBY3QpdYp3Aove/9x3st0br7VsGFDli5disPhIC8vD4C5c+e6FiGJBAON0KXWKZxv97YixkUbaEgNp6X/Umu5tw0ASm0dMHv2bPjuUx7q4OWuimojIH6kpf8iHriXOpbVOqCwAuahP63x7sYqc5QAUUKXWsu91LGs6Zc1a7xM5CIBpoQu4oVXX30VgEcffdS7H1BbXgkAVbmIeGH58uUsX77cu4u7DnbOoYMzsavkUaqJRugiBcrqBXP/xDe8XomKexsBzadLNdIIXYTztemeFNarT58+nbi4ONfrU6dOZdCgQa7jlJQUEhIS/B6rSGk0Qheh7F4wFa5XL059YKSaqA5dpByFCX3Rwz29uv75558HYNy4cc5EvmIUtLlBdeniE2XVoWvKRcQLm3cfYcHmPV5du2PHDnbs2OE8iE5yJnORaqApF5FyDIxqxebdR3gvc59XD0bnz58PwPjx49m7dy+pN/o7QhEnJXSRciTGtHY18aq0ivSBqS6qjw85SugifjJp0iTnNxmpgQ3EE22rF5KU0EW8VLE9S52NvsaOHQvA5Mk17IFoTfttQXxCCV3EC5XZszQxpjWHDx/2V0giJSihi3ihMnuWAsyaNctfIYmUoLJFET8bPXo0o0ePDnQYJRUueJKQoYQu4geF8+3xb2ziq+8O8vPPPwc6pKK08XVI0pSLiI8V3zijc59HmOnlKtNqo42vQ5JG6CI+lhjTmkUP92TRwz1dDb9GjRrFqFGjXNcMHz68yDRMcnKyqyIGICkpifHjx7uOhwwZ4mopAJCQkEBKSorreNCgQUydOtV1HBcXx/Tp013H/fr1c/V0B4iNjWVWmncrXyV4KKGL1EI5OTkYE+goxNfUnEvEjyra2KtapfY/v6F1cVpFWmOpOZdIAFWksVe1ct9ZyZ12WQpaXj0UNcb0BaYDdYA3rbUpHq65B5gIWOAra22iD+MUCUoVbexVXRwOBwBpaWklT2oVadAqN6EbY+oAM4E+QA6QboxZZq3NdrumIzAW6GWt/dEYc5m/AhYJJj5p7OUHQ4cODXQI4gfejNCvA3Zaa3cBGGMWAgOBbLdrHgJmWmt/BLDWHvR1oCLiO0OHDiUvL4+8vDzq1asX6HDER7yZQ28F7HU7zil4zd2VwJXGmE+NMZ8VTNGUYIxJNsZkGGMyDh06VLmIRcQn+vTpQ58+fQIdhviQrxYW1QU6Ag4gHFhvjOlqrT3qfpG1dhYwC5xVLj56b5EaryKdGqtLZ8cdXN8hzPNJT/3bVflS43mT0PcBEW7H4QWvucsBNltr84DdxpivcSb4dJ9EKRLEKtKpsbpk7z9O55YxDBnioZyysC2AO/VPDwreJPR0oKMxph3ORJ4AFK9geRe4F0g1xoThnILZ5ctARYJVRTo1Vpf4NzZx7uxpxo8fz9VXX018fDzHjh1j4MCBjBw5kruKb2itypegUO4curX2HDACWAVsA9621m41xkwyxsQVXLYKOGyMyQbWAk9ba9UIWqQG+37zStavXx/oMMSHtFJUpBYqawVrbm4uAGFhbvPrhSP04iN3qXZaKSoiJZS2gnXw4MEMHlzKPHpqf+eX+qjXSGqfK1ILlbWC9amnnir5A+4PSvWAtMZSQhephcpawXr77beXfDE66XwC1wPSGktTLiJSxIEDBxg1ahQrV64EYO/evTgcDu2PGgSU0EWkiGXLlpGZmVnktbp169K+ffsARSTe0pSLiBSRnJxMcnKy6zgiIoI1a9YEMCLxlhK6SC1WWkuCgVGtijws3bXLuU7QNUr31BqgOLUKqHZK6CK1VGktCbL3HwcoktCHDRsGFPRP99QaoDhVwgSEErpILVVaSwJPI/bnnnvu/IF7xUtpVAkTEHooKiLl6t27Ny1atMDhcLBx48ZAhyOlUEIXEQkRSugi4pXIyEjS0tJo3LgxWVlZgQ5HPNAcuoiU4Kn6pbDyZcSIEUApG0xLQCmhi0gRnqpf3CtfpkyZUt0hiZeU0EWkCE/VL+6j9R49enh3o+K16qpL9zvNoYtIhWRmZpKQkOBqD5Ceno7D4XDNq2/cuBHHK98ya3uT8z90YAtsWRyIcGsVjdBFpEI+//xzDhw4UOY17br1pPeYMRAZ6XxBdenVQjsWiUi5ytrhyCva8chnytqxSCN0EfGKe+VL8V4vxa1btw5wLkiS6qOELiLlcq988dTrpbgJEyYAKm2sbkroIlIu98oXT71eipszZ46/QxIPlNBFxOc8bobhTctdUHljFSihi4jPFW6IERsb63zBm5a7oLa7VaSELiI+98ILL3D11VcTGxvLypUrSUmZR2JiYpGdkDxSeWOVKKGLiM8lJibSr18/13FUVBRxcXEBjKh2UEIXkQorbeu687py8b8tiRHQt29f+vbtW22x1WZK6CJSIaVtXeeueGnj8uXLAbj99tv9F5gooYtIxZS2dZ274qP3l19+GVBC9zcldBHxu8WL1ZirOiihi4jfhYWFBTqEWkHtc0XE75YuXcrSpUsDHUbIU0IXEb+bMWMG6enpACxatAiHw8GpU6cAmD9/Pg6Hg7y8vECGGBI05SIifpeYmEh8fHy5183dkAPAUC0UrRT1QxcRn4t/YxPZ+4/TuWVj12vltdwFcFzVHIC07Yf9Gl8wUz90EalWxWvVvWm5C5D2zG/8FlNtoIQuIj5XvFbdm5a7UnV6KCoiNcbsdXuYvW5PoMMIWhqhi0iNsejz/QA8FOA4gpUSuojUGGuejgl0CEHNq4RujOkLTAfqAG9aa1NKuW4QsBjoYa1VCYuIVJy3Oxt5qxbtgFRuQjfG1AFmAn2AHCDdGLPMWptd7LpGwOPAZn8EKiLBrfyWu1Dno3y6mQaMuc1Hb1rLdkDyZoR+HbDTWrsLwBizEBgIZBe77nngJeBpn0YoIkHPm5a7AMu/ymXtry5hzN/f980b17IdkLxJ6K2AvW7HOUCRiS5jzDVAhLX2fWNMqQndGJMMJAO0bl12PaqIhA5vWu4CxDOjGqIJXVUuWzTGXAD8D/BUeddaa2dZa6OttdGXXnppVd9aRETceDNC3wdEuB2HF7xWqBHQBUgzxgC0AJYZY+L0YFREKuLrfyxyfvNwz8AGEqS8SejpQEdjTDuciTwBSCw8aa09BriaHRtj0oDRSuYiUlEHd3wR6BCCWrkJ3Vp7zhgzAliFs2xxjrV2qzFmEpBhrV3m7yBFpHa44dH/DnQIQc2rOnRr7QfAB8VeG1/KtY6qhyUiIhWllaIiUmNsX73A+Y3m0CtFzblEpMY4vDuLw7uzAEhJSSEhIcF17vnnn2fIkCGu4/Hjx5OUdH7B0NixY0lOTnYdjx49muFvZVVD1DWHRugiUmP0evjFQIcQ1LRjkYjUGIWtARb5asoltT98vwEGTAuZ5f9l7VikKRcRCVmjVp9l1MrTsGVxoEOpFkroIhK6Lr8aGv060FFUG82hi8XBuoEAAAfPSURBVEjImjZtGqR+E+gwqo1G6CIiIUIJXURC1vDhw2tV6aKmXEQkZDVo0AAurBPoMKqNErqIhKypU6dC6rZAh1FtlNBFpEbxZqu6iphx5jSXN6rvs/vVZJpDF5EaY2BUKzq3bOyz+30063lGvLXVZ/er6TRCF5Eaw9ut6rzV6Z0wmvxce9Jc7fmkIlLrdLvzEZ44vCvQYVQbJXQRCWnPLvmGppwglf5Vu1HXwTW+H4wSuoiEtLzmVxJR/0DVbnJgi/OfSugiIoFz0aD/ZlJVuzemVnF0X01U5SIiEiKU0EUkpH02Z2KRnY5CmaZcRCSkNbq8DZGREYEOo1oooYtISLu6fxLjasmm05pyEREJEUroIhLSVkx7hp594gIdRrVQQheRkDUwqhUXXt6ekxe3AmDQoEHODowF4uLimD59uuu4X79+vPrqq67j2NhYZs+eXX0BV5Hm0EUkZCXGtOa9ex8OdBjVxlhrA/LG0dHRNiMjIyDvLSK1R2Er3kVVeDDquKo5AGnbD/skpqowxnxhrY32dE4jdBGRcgztFR7oELyihC4iUo6hNyihi4iEhLxz+fBDFvVqeE8XVbmIiJSjz+u76TP/VKDDKJdG6CIi5fjDE886v6kJPWGGmVJPKaGLiJQjWJp7KaGLSMjL3n/cVb5YGefOnmZAt1+T1DvSh1H5nhK6iIS0gVGtqnyPlS8/zicX1iXp6y99EJH/KKGLSEhLjGlNYkzrKt2jZ8bdPorGv5TQRUTK0To6NtAheEUJXUSkHGd/PhHoELyiOnQRkXJ8+uozrJ/2OAC5ubk4HA6WL18OwIEDB3A4HKxcuRKAvXv34nA4WLNmDQC7du3C4XCwbt06AHbs2IHD4WDjxo0AZGVl4XA4SE9PByAzMxOHw0FmZiYA6enpOBwOsrKyyo3Tq4RujOlrjNlhjNlpjBnj4fyTxphsY8y/jDH/MMa08ea+IiLBoPV1fbiq7/2BDqNc5XZbNMbUAb4G+gA5QDpwr7U22+2am4DN1tpTxphHAIe1Nr6s+6rboogEC190bPSVsrotejNCvw7Yaa3dZa09CywEBrpfYK1da60tXBf7GRAcnWxEREKINwm9FbDX7Tin4LXSPAh86OmEMSbZGJNhjMk4dOiQ91GKiEi5fPpQ1BgzBIgGpng6b62dZa2NttZGX3rppb58axGRWs+bssV9QITbcXjBa0UYY2KB/wJ6W2vP+CY8ERHxljcj9HSgozGmnTHmQiABWOZ+gTGmO/AGEGetPej7MEVEpDzlJnRr7TlgBLAK2Aa8ba3daoyZZIyJK7hsCnAx8HdjTKYxZlkptxMRET/xaqWotfYD4INir413+z441sWKiFRSVTs2Vgct/RcRKYcvOjZWh3IXFvmLFhaJiFRcVRcWiYhIEFBCFxEJEUroIiIhQgldRCREKKGLiIQIJXQRkRChhC4iEiKU0EVEQkTAFhYZYw4B3wfkzasmDMgNdBB+EsqfDfT5glkofzao2OdrY6312H88YAk9WBljMkpbpRXsQvmzgT5fMAvlzwa++3yachERCRFK6CIiIUIJveJmBToAPwrlzwb6fMEslD8b+OjzaQ5dRCREaIQuIhIilNBFREKEEroXjDERxpi1xphsY8xWY8zjgY7JH4wxdYwx/zTGrAh0LL5mjGlqjFlsjNlujNlmjOkZ6Jh8xRjzRMHfyyxjzN+MMfUDHVNVGGPmGGMOGmOy3F67xBjzkTHmm4J/NgtkjFVRyuebUvB381/GmHeMMU0rc28ldO+cA56y1nYGfgMMN8Z0DnBM/vA4zo3AQ9F0YKW19irgPwiRz2mMaQWMBKKttV2AOkBCYKOqsrlA32KvjQH+Ya3tCPyj4DhYzaXk5/sI6GKt7QZ8DYytzI2V0L1grd1vrf2y4PufcCaD4Nhk0EvGmHCgP/BmoGPxNWNME+BG4C8A1tqz1tqjgY3Kp+oCDYwxdYGGwL8DHE+VWGvXA0eKvTwQmFfw/TzgjmoNyoc8fT5r7Wpr7bmCw8+A8MrcWwm9gowxbYHuwObARuJz04D/BPIDHYgftAMOAakFU0pvGmMuCnRQvmCt3QdMBfYA+4Fj1trVgY3KLy631u4v+P4AcHkgg/GzYcCHlflBJfQKMMZcDCwBRllrjwc6Hl8xxgwADlprvwh0LH5SF7gGeM1a2x04SXD/yu5SMJc8EOf/tH4NXGSMGRLYqPzLOmutQ7Le2hjzXzineP+vMj+vhO4lY0w9nMn8/6y1SwMdj4/1AuKMMd8BC4GbjTHzAxuST+UAOdbawt+qFuNM8KEgFthtrT1krc0DlgLXBzgmf/jBGNMSoOCfBwMcj88ZY4YCA4D7bCUXCCmhe8EYY3DOv26z1v5PoOPxNWvtWGttuLW2Lc4Hah9ba0NmlGetPQDsNcZEFrx0C5AdwJB8aQ/wG2NMw4K/p7cQIg98i1kGPFDw/QPAewGMxeeMMX1xTnnGWWtPVfY+Suje6QXcj3PkmlnwdVugg5IKeQz4P2PMv4Ao4MUAx+MTBb91LAa+BLbg/G86qJfJG2P+BmwCIo0xOcaYB4EUoI8x5hucv5WkBDLGqijl870CNAI+Ksgvr1fq3lr6LyISGjRCFxEJEUroIiIhQgldRCREKKGLiIQIJXQRkRChhC4iEiKU0EVEQsT/B7YjD2mqmcoIAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"cph.plot_covariate_groups('trt', values=[0, 1]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice how the group without treatment has a lower survival rate at all times (the x-axis is time) compared to the treatment group."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a name='6'></a>\n",
"## 6. Hazard Ratio\n",
"\n",
"Recall from the lecture videos that the Hazard Ratio between two patients was the likelihood of one patient (e.g smoker) being more at risk than the other (e.g non-smoker).\n",
"$$\n",
"\\frac{\\lambda_{smoker}(t)}{\\lambda_{nonsmoker}(t)} = e^{\\theta (X_{smoker} - X_{nonsmoker})^T}\n",
"$$\n",
"\n",
"Where\n",
"\n",
"$$\n",
"\\lambda_{smoker}(t) = \\lambda_0(t)e^{\\theta X_{smoker}^T}\n",
"$$\n",
"and\n",
"$$\n",
"\\lambda_{nonsmoker}(t) = \\lambda_0(t)e^{\\theta X_{nonsmoker}^T} \\\\\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a name='Ex-2'></a>\n",
"### Exercise 2\n",
"In the cell below, write a function to compute the hazard ratio between two individuals given the cox model's coefficients."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<details> \n",
"<summary>\n",
" <font size=\"3\" color=\"darkgreen\"><b>Hints</b></font>\n",
"</summary>\n",
"<p>\n",
"<ul>\n",
" <li>use numpy.dot</li>\n",
" <li>use nump.exp</li>\n",
"</ul>\n",
"</p>\n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "WbBmxbeDA3k1"
},
"outputs": [],
"source": [
"# UNQ_C2 (UNIQUE CELL IDENTIFIER, DO NOT EDIT)\n",
"def hazard_ratio(case_1, case_2, cox_params):\n",
" '''\n",
" Return the hazard ratio of case_1 : case_2 using\n",
" the coefficients of the cox model.\n",
" \n",
" Args:\n",
" case_1 (np.array): (1 x d) array of covariates\n",
" case_2 (np.array): (1 x d) array of covariates\n",
" model (np.array): (1 x d) array of cox model coefficients\n",
" Returns:\n",
" hazard_ratio (float): hazard ratio of case_1 : case_2\n",
" '''\n",
" \n",
" ### START CODE HERE (REPLACE INSTANCES OF 'None' with your code) ###\n",
" \n",
" hr = np.exp(cox_params.dot((case_1 - case_2).T))\n",
" \n",
" ### END CODE HERE ###\n",
" \n",
" return hr"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "zbDQUxE6CcA3"
},
"source": [
"Now, evaluate it on the following pair of indivduals: `i = 1` and `j = 5`"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "7flsvTRXCgqO"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"15.029017732492221\n"
]
}
],
"source": [
"i = 1\n",
"case_1 = one_hot_train.iloc[i, :].drop(['time', 'status'])\n",
"\n",
"j = 5\n",
"case_2 = one_hot_train.iloc[j, :].drop(['time', 'status'])\n",
"\n",
"print(hazard_ratio(case_1.values, case_2.values, cph.params_.values))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Expected Output:\n",
"```CPP\n",
"15.029017732492221\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Question:** \n",
"\n",
"Is `case_1` or `case_2` at greater risk? "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<details> \n",
"<summary>\n",
" <font size=\"3\" color=\"darkgreen\"><b>Check your answer!</b></font>\n",
"</summary>\n",
"<p>\n",
"<ul>\n",
"<ul>\n",
" Important! The following answer only applies if you picked i = 1 and j = 5\n",
" <li>You should see that `case_1` is at higher risk.</li>\n",
" <li>The hazard ratio of case 1 / case 2 is greater than 1, so case 1 had a higher hazard relative to case 2</li>\n",
"</ul>\n",
"</p>"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Inspect different pairs, and see if you can figure out which patient is more at risk."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 374
},
"colab_type": "code",
"id": "g2PZ3sGvCs0K",
"outputId": "59336868-d421-4645-d88e-76a8a8cffc9f"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Case 1\n",
"\n",
" trt 0.000000\n",
"age 0.563645\n",
"sex 0.000000\n",
"ascites 0.000000\n",
"hepato 1.000000\n",
"spiders 1.000000\n",
"bili -0.405645\n",
"chol -0.210436\n",
"albumin 1.514297\n",
"copper -0.481961\n",
"alk.phos 2.173526\n",
"ast -0.144699\n",
"trig -0.531125\n",
"platelet -0.450972\n",
"protime -0.139881\n",
"edema_0.5 0.000000\n",
"edema_1.0 0.000000\n",
"stage_2.0 0.000000\n",
"stage_3.0 1.000000\n",
"stage_4.0 0.000000\n",
"Name: 1, dtype: float64 \n",
"\n",
"Case 2\n",
"\n",
" trt 0.000000\n",
"age 0.463447\n",
"sex 0.000000\n",
"ascites 0.000000\n",
"hepato 1.000000\n",
"spiders 0.000000\n",
"bili -0.489581\n",
"chol -0.309875\n",
"albumin -1.232371\n",
"copper -0.504348\n",
"alk.phos 2.870427\n",
"ast -0.936261\n",
"trig -0.150229\n",
"platelet 3.190823\n",
"protime -0.139881\n",
"edema_0.5 0.000000\n",
"edema_1.0 0.000000\n",
"stage_2.0 0.000000\n",
"stage_3.0 0.000000\n",
"stage_4.0 1.000000\n",
"Name: 38, dtype: float64 \n",
"\n",
"Hazard Ratio: 0.1780450006997129\n"
]
}
],
"source": [
"i = 4\n",
"case_1 = one_hot_train.iloc[i, :].drop(['time', 'status'])\n",
"\n",
"j = 7\n",
"case_2 = one_hot_train.iloc[j, :].drop(['time', 'status'])\n",
"\n",
"print(\"Case 1\\n\\n\", case_1, \"\\n\")\n",
"print(\"Case 2\\n\\n\", case_2, \"\\n\")\n",
"print(\"Hazard Ratio:\", hazard_ratio(case_1.values, case_2.values, cph.params_.values))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<details> \n",
"<summary>\n",
" <font size=\"3\" color=\"darkgreen\"><b>Check your answer!</b></font>\n",
"</summary>\n",
"<p>\n",
"<ul>\n",
"<ul>\n",
" Important! The following answer only applies if you picked i = 4 and j = 7\n",
" <li>You should see that `case_2` is at higher risk.</li>\n",
" <li>The hazard ratio of case 1 / case 2 is less than 1, so case 2 had a higher hazard relative to case 1</li>\n",
"</ul>\n",
"</p>"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "KUa6r-KOyySp"
},
"source": [
"<a name='7'></a>\n",
"## 7. Harrell's C-index"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "woQAtSmRXrgr"
},
"source": [
"To evaluate how good our model is performing, we will write our own version of the C-index. Similar to the week 1 case, C-index in the survival context is the probability that, given a randomly selected pair of individuals, the one who died sooner has a higher risk score. \n",
"\n",
"However, we need to take into account censoring. Imagine a pair of patients, $A$ and $B$. \n",
"\n",
"#### Scenario 1\n",
"- A was censored at time $t_A$ \n",
"- B died at $t_B$\n",
"- $t_A < t_B$. \n",
"\n",
"Because of censoring, we can't say whether $A$ or $B$ should have a higher risk score. \n",
"\n",
"#### Scenario 2\n",
"Now imagine that $t_A > t_B$.\n",
"\n",
"- A was censored at time $t_A$ \n",
"- B died at $t_B$\n",
"- $t_A > t_B$\n",
"\n",
"Now we can definitively say that $B$ should have a higher risk score than $A$, since we know for a fact that $A$ lived longer. \n",
"\n",
"Therefore, when we compute our C-index\n",
"- We should only consider pairs where at most one person is censored\n",
"- If they are censored, then their censored time should occur *after* the other person's time of death. \n",
"\n",
"The metric we get if we use this rule is called **Harrel's C-index**.\n",
"\n",
"Note that in this case, being censored at time $t$ means that the true death time was some time AFTER time $t$ and not at $t$. \n",
"- Therefore if $t_A = t_B$ and A was censored:\n",
" - Then $A$ actually lived longer than $B$. \n",
" - This will effect how you deal with ties in the exercise below!\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a name='Ex-3'></a>\n",
"### Exercise 3\n",
"Fill in the function below to compute Harrel's C-index."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<details> \n",
"<summary>\n",
" <font size=\"3\" color=\"darkgreen\"><b>Hints</b></font>\n",
"</summary>\n",
"<p>\n",
"<ul>\n",
" <li>If you get a division by zero error, consider checking how you count when a pair is permissible (in the case where one patient is censored and the other is not censored).</li>\n",
"</ul>\n",
"</p>"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"# UNQ_C3 (UNIQUE CELL IDENTIFIER, DO NOT EDIT)\n",
"def harrell_c(y_true, scores, event):\n",
" '''\n",
" Compute Harrel C-index given true event/censoring times,\n",
" model output, and event indicators.\n",
" \n",
" Args:\n",
" y_true (array): array of true event times\n",
" scores (array): model risk scores\n",
" event (array): indicator, 1 if event occurred at that index, 0 for censorship\n",
" Returns:\n",
" result (float): C-index metric\n",
" '''\n",
" \n",
" n = len(y_true)\n",
" assert (len(scores) == n and len(event) == n)\n",
" \n",
" concordant = 0.0\n",
" permissible = 0.0\n",
" ties = 0.0\n",
" \n",
" result = 0.0\n",
" \n",
" ### START CODE HERE (REPLACE INSTANCES OF 'None' and 'pass' with your code) ###\n",
" \n",
" # use double for loop to go through cases\n",
" for i in range(n):\n",
" # set lower bound on j to avoid double counting\n",
" for j in range(i+1, n):\n",
" \n",
" # check if at most one is censored\n",
" if event[i] == 1 or event[j] == 1:\n",
" \n",
" # check if neither are censored\n",
" if event[i] == 1 and event[j] == 1:\n",
" \n",
" permissible += 1.0\n",
" \n",
" # check if scores are tied\n",
" if scores[i] == scores[j]:\n",
" ties += 1.0\n",
" \n",
" # check for concordant\n",
" elif y_true[i] < y_true[j] and scores[i] > scores[j]:\n",
" concordant += 1.0\n",
" elif y_true[i] > y_true[j] and scores[i] < scores[j]:\n",
" concordant += 1.0\n",
" \n",
" # check if one is censored\n",
" elif event[i] != event[j]:\n",
" \n",
" # get censored index\n",
" censored = j\n",
" uncensored = i\n",
" \n",
" if event[i] == 0:\n",
" censored = i\n",
" uncensored = j\n",
" \n",
" # check if permissible\n",
" # Note: in this case, we are assuming that censored at a time\n",
" # means that you did NOT die at that time. That is, if you\n",
" # live until time 30 and have event = 0, then you lived THROUGH\n",
" # time 30.\n",
" if y_true[uncensored] <= y_true[censored]:\n",
" permissible += 1.0\n",
" \n",
" # check if scores are tied\n",
" if scores[uncensored] == scores[censored]:\n",
" # update ties \n",
" ties += 1.0\n",
" \n",
" # check if scores are concordant \n",
" if scores[uncensored] > scores[censored]:\n",
" concordant += 1.0\n",
" \n",
" # set result to c-index computed from number of concordant pairs,\n",
" # number of ties, and number of permissible pairs (REPLACE 0 with your code) \n",
" result = (concordant + 0.5*ties) / permissible\n",
" \n",
" ### END CODE HERE ###\n",
" \n",
" return result "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can test your function on the following test cases:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Case 1\n",
"Expected: 1.0, Output: 1.0\n",
"\n",
"Case 2\n",
"Expected: 0.0, Output: 0.0\n",
"\n",
"Case 3\n",
"Expected: 1.0, Output: 1.0\n",
"\n",
"Case 4\n",
"Expected: 0.75, Output: 0.75\n",
"\n",
"Case 5\n",
"Expected: 0.583, Output: 0.5833333333333334\n",
"\n",
"Case 6\n",
"Expected: 1.0 , Output:1.0000\n"
]
}
],
"source": [
"y_true = [30, 12, 84, 9]\n",
"\n",
"# Case 1\n",
"event = [1, 1, 1, 1]\n",
"scores = [0.5, 0.9, 0.1, 1.0]\n",
"print(\"Case 1\")\n",
"print(\"Expected: 1.0, Output: {}\".format(harrell_c(y_true, scores, event)))\n",
"\n",
"# Case 2\n",
"scores = [0.9, 0.5, 1.0, 0.1]\n",
"print(\"\\nCase 2\")\n",
"print(\"Expected: 0.0, Output: {}\".format(harrell_c(y_true, scores, event)))\n",
"\n",
"# Case 3\n",
"event = [1, 0, 1, 1]\n",
"scores = [0.5, 0.9, 0.1, 1.0]\n",
"print(\"\\nCase 3\")\n",
"print(\"Expected: 1.0, Output: {}\".format(harrell_c(y_true, scores, event)))\n",
"\n",
"# Case 4\n",
"y_true = [30, 30, 20, 20]\n",
"event = [1, 0, 1, 0]\n",
"scores = [10, 5, 15, 20]\n",
"print(\"\\nCase 4\")\n",
"print(\"Expected: 0.75, Output: {}\".format(harrell_c(y_true, scores, event)))\n",
"\n",
"# Case 5\n",
"y_true = list(reversed([30, 30, 30, 20, 20]))\n",
"event = [0, 1, 0, 1, 0]\n",
"scores = list(reversed([15, 10, 5, 15, 20]))\n",
"print(\"\\nCase 5\")\n",
"print(\"Expected: 0.583, Output: {}\".format(harrell_c(y_true, scores, event)))\n",
"\n",
"# Case 6\n",
"y_true = [10,10]\n",
"event = [0,1]\n",
"scores = [4,5]\n",
"print(\"\\nCase 6\")\n",
"print(f\"Expected: 1.0 , Output:{harrell_c(y_true, scores, event):.4f}\")"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "CtQVe4pAn8ic"
},
"source": [
"Now use the Harrell's C-index function to evaluate the cox model on our data sets."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
},
"colab_type": "code",
"id": "8nzHc_Qbn7dM",
"outputId": "bc2f960d-16e5-46b2-a41f-695892c311c7"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Train: 0.8265139116202946\n",
"Val: 0.8544776119402985\n",
"Test: 0.8478543563068921\n"
]
}
],
"source": [
"# Train\n",
"scores = cph.predict_partial_hazard(one_hot_train)\n",
"cox_train_scores = harrell_c(one_hot_train['time'].values, scores.values, one_hot_train['status'].values)\n",
"# Validation\n",
"scores = cph.predict_partial_hazard(one_hot_val)\n",
"cox_val_scores = harrell_c(one_hot_val['time'].values, scores.values, one_hot_val['status'].values)\n",
"# Test\n",
"scores = cph.predict_partial_hazard(one_hot_test)\n",
"cox_test_scores = harrell_c(one_hot_test['time'].values, scores.values, one_hot_test['status'].values)\n",
"\n",
"print(\"Train:\", cox_train_scores)\n",
"print(\"Val:\", cox_val_scores)\n",
"print(\"Test:\", cox_test_scores)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What do these values tell us ?"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "AuNjR_wNkpWz"
},
"source": [
"<a name='8'></a>\n",
"## 8. Random Survival Forests\n",
"\n",
"This performed well, but you have a hunch you can squeeze out better performance by using a machine learning approach. You decide to use a Random Survival Forest. To do this, you can use the `RandomForestSRC` package in R. To call R function from Python, we'll use the `r2py` package. Run the following cell to import the necessary requirements. \n"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 530
},
"colab_type": "code",
"id": "ZgSy-Dj6kquK",
"outputId": "4aa5d2fa-30f4-4328-ae29-a2ff05223e22"
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"R[write to console]: Loading required package: ggplot2\n",
"\n"
]
}
],
"source": [
"%load_ext rpy2.ipython\n",
"%R require(ggplot2)\n",
"\n",
"from rpy2.robjects.packages import importr\n",
"# import R's \"base\" package\n",
"base = importr('base')\n",
"\n",
"# import R's \"utils\" package\n",
"utils = importr('utils')\n",
"\n",
"# import rpy2's package module\n",
"import rpy2.robjects.packages as rpackages\n",
"\n",
"forest = rpackages.importr('randomForestSRC', lib_loc='R')\n",
"\n",
"from rpy2 import robjects as ro\n",
"R = ro.r\n",
"\n",
"from rpy2.robjects import pandas2ri\n",
"pandas2ri.activate()"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "LXBgqQBfuMA5"
},
"source": [
"Instead of encoding our categories as binary features, we can use the original dataframe since trees deal well with raw categorical data (can you think why this might be?).\n",
"\n",
"Run the code cell below to build your forest."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"colab": {},
"colab_type": "code",
"id": "B-pio4o4mdVJ"
},
"outputs": [],
"source": [
"model = forest.rfsrc(ro.Formula('Surv(time, status) ~ .'), data=df_train, ntree=300, nodedepth=5, seed=-1)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 289
},
"colab_type": "code",
"id": "zZfcUvJ3nL04",
"outputId": "27d00bd8-ea33-4c1b-f5d7-ca2c73ead721"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" Sample size: 154\n",
" Number of deaths: 64\n",
" Number of trees: 300\n",
" Forest terminal node size: 15\n",
" Average no. of terminal nodes: 6.54\n",
"No. of variables tried at each split: 5\n",
" Total no. of variables: 17\n",
" Resampling used to grow trees: swor\n",
" Resample size used to grow trees: 97\n",
" Analysis: RSF\n",
" Family: surv\n",
" Splitting rule: logrank *random*\n",
" Number of random split points: 10\n",
" Error rate: 19.07%\n",
"\n",
"\n"
]
}
],
"source": [
"print(model)"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "9Mwzm55H-QKV"
},
"source": [
"Finally, let's evaluate on our validation and test sets, and compare it with our Cox model."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
},
"colab_type": "code",
"id": "vfl4LbGfpbKp",
"outputId": "13f8b560-e171-41e9-f6dc-cf0468a3f786"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Cox Model Validation Score: 0.8544776119402985\n",
"Survival Forest Validation Score: 0.8296019900497512\n"
]
}
],
"source": [
"result = R.predict(model, newdata=df_val)\n",
"scores = np.array(result.rx('predicted')[0])\n",
"\n",
"print(\"Cox Model Validation Score:\", cox_val_scores)\n",
"print(\"Survival Forest Validation Score:\", harrell_c(df_val['time'].values, scores, df_val['status'].values))"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
},
"colab_type": "code",
"id": "uhqSQJhrplSG",
"outputId": "752c266e-0234-45c5-d53f-554e2ff17a5a"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Cox Model Test Score: 0.8478543563068921\n",
"Survival Forest Validation Score: 0.8621586475942783\n"
]
}
],
"source": [
"result = R.predict(model, newdata=df_test)\n",
"scores = np.array(result.rx('predicted')[0])\n",
"\n",
"print(\"Cox Model Test Score:\", cox_test_scores)\n",
"print(\"Survival Forest Validation Score:\", harrell_c(df_test['time'].values, scores, df_test['status'].values))"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "Gp_SgUXreAWn"
},
"source": [
"Your random forest model should be outperforming the Cox model slightly. Let's dig deeper to see how they differ."
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "ZtPMPaSli8GB"
},
"source": [
"<a name='9'></a>\n",
"## 9. Permutation Method for Interpretation\n",
"\n",
"We'll dig a bit deeper into interpretation methods for forests a bit later, but for now just know that random surival forests come with their own built in variable importance feature. The method is referred to as VIMP, and for the purpose of this section you should just know that higher absolute value of the VIMP means that the variable generally has a larger effect on the model outcome.\n",
"\n",
"Run the next cell to compute and plot VIMP for the random survival forest."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 281
},
"colab_type": "code",
"id": "u7M4_N_d-YJu",
"outputId": "7e1830cb-4b67-444f-8ba5-d49d3ff2f172"
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEICAYAAAB4YQKYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZwdZZn28d+VgIQtCZCogMQWBNmJ2KAsKiI6CI4ioCggImoGQXjxHVQcmREcF5R3xo1xiQyyBDfCgA5oQGVRUAwJCVmIoJIogiJhCUtYQnK9f9TTcjjp7vRytg7X9/PpT9epeuqp+5x0+u6nnqq7ZJuIiIjhGtXuACIiYu2QhBIREQ2RhBIREQ2RhBIREQ2RhBIREQ2RhBIREQ2RhBLPeZI+J+mUAba9TtL7G3z8MyRNa2SfzSBpiaQDGtDPSZI+34iYorMkoURbSJoh6VO9rH+rpL9KWkfS+ZI+XdZ3SbKkOXXtJ0h6StKSmnVLJD0u6VFJ95Z+NuojjonAMcA3G/oGm2SkJJ81+BZwlKTntzuQaKwklGiXC4CjJalu/buBi20/3cd+G0jaueb1kcDiXtr9o+2NgN2BbuD0Pvo7Fvix7ccHHHkMi+0ngJ9QJfJYiyShRLtcDmwGvLpnhaRNgDcDF/az30XAe2peH9Nfe9t3U/3y2rmPJm8Crq+NQdIVku6T9GBZflHdPttIminpYUk/lLRp2XeMpGmS7pf0kKSbJb2gbNtC0o8kPSDp95I+0FswkvaT9Oe6dUskHSDpQOBfgCPK6OvWsn2cpP+W9BdJd0v6tKTRvfS9RRm5bVqz7uWSlkpaV9I2kq4p8S+VdLGk8X3E+ffRY29xl2NdWj7HxZJOruviOuDg3vqOkSsJJdqijAh+wLP/Sn0H8Fvbt/az6zTgnZJGS9oR2Aj4TV+NJW0FHATM6aPJLsDtNa9HAd8GXgxMAh4Hzqnb5xjgOGBz4GngK2X9e4BxwFZUyfL4sj/A94A/A1sAhwOflbR/P+9zNbZnAJ8Fvm97I9u7lU3nlzheCrwceCOw2jyP7XuAXwOH1aw+EphuewUg4HMlxh3K+zhjMDECSBoF/C9wK7Al8HrgFEn/UNNsEbBbL7vHCJaEEu10AXC4pDHl9TFlXX/+TJUADijtL+qj3eWSHgJuoBqBfLaPduOBR3pe2L7f9qW2l9t+BPgM8Nq6fS6yvcD2Y8C/Au8oI4IVVInkpbZX2p5t++GS1PYBPmb7CdtzgXNpwCmfMgI6CDjF9mO2/wZ8EXhnH7t8B3hX2Vel3XfKe/+97Z/aftL2fcB/9vLeB2IPYKLtT9l+yvadVPMmtTE9QpV8Yy2yTrsDiOcu2zdIWgocIulmYE/g0AHseiHV3MfeVKfMtuulzSG2fzaAvh4ENu55IWkDql/IBwKblNUbSxpte2V5fVfN/n8E1gUmUCW3rYDvlVNF04BPUP3F/0BJULX7dQ8gvjV5cTn+X2qmo0bVxVjrUuCrkjan+txWAb+EvyenL1N9phuXfh4cYkxblITeY3TPcYqNgWVD6Ds6WEYo0W4XUv2lfjRwle17B7DPpVTn3++0/adhHn8ez05I/wy8DHil7bHAa8r62osHtqpZnkQ1Mllqe4XtM23vSJXs3kz13u4BNpW0cd1+d/cSz2PABj0vyshnYs32+vLgdwFPAhNsjy9fY23v1Nubtf0gcDVwBNXpru/5mZLjny3971Le+9F177vPOIEX1sW0uCae8bY3tn1QTZsdqE6JxVokCSXa7UKq01cfYM2nuwAop5r2p5d5giH4Mc8+rbMx1bzHQ2Xy+pO97HO0pB3LaOZTVHMQKyW9TtIuJQk8TJVoVtm+C/gV8Lkycb8r8D6qEUy9O4Axkg6WtC7V1Wnr1Wy/F+gq8xTY/gtVgvgPSWMljSqT6/2dqvoOVaI7vCzXvvdHgWWStgQ+0k8fc4GDJG0q6YVA7X08M4FHJH1M0vplvmtnSXvUtHkt1cUSsRZJQom2sr2E6pfthsCPBrHfLNt/aEAIF1L9Yly/vP4SsD6wFLgJmNHLPhdRTYT/FRgD9FzB9EJgOlUyWUQ1d9Mzx/MuoItqtHIZ8MneTsnZXgacQDXHcjfVSKD2qq9Lyvf7Jd1Slo8BngfcRnWKajrVBQN9+RGwLfDXugsgzqS6zHoZcCXwP/30cRHVCGMJVUL7fs17WEk1OptMdUn30vJ+xkF1NRzVvM+A/oCIkUN5wFY810n6LPA3219qdyzPBZJOAray/dF2xxKNlYQSERENkVNeERHREEkoERHREEkoERHREGv9jY0TJkxwV1dXu8OIiBhRZs+evdT2xDW3fMZan1C6urqYNWtWu8OIiBhRJP1xsPvklFdERDREEkpERDREEkpERDREEkpERDREEkpERDREEkpERDREEkpERDREEkpERDTEWp9Q5t+9jK7TrqTrtCvbHUpExFptrU8oERHRGi1NKJJOKY9NjYiItUyrRyinAEkoERFroaYlFEkbSrpS0q2SFkj6JLAFcK2ka0ubr0uaJWmhpDNr9j1I0m8lzZb0FUlX1PR5nqSZkuZIemuz4o+IiMFpZrXhA4F7bB8MIGkc8F7gdbaXljafsP2ApNHAzyXtCtwBfBN4je3Fkr5b0+cngGtsHydpPDBT0s9sP1Z7YElTgCkAo8cOqvpyREQMUTNPec0H3iDp85JebXtZL23eIekWYA6wE7AjsD1wp+3FpU1tQnkjcJqkucB1wBhgUn2ntqfa7rbdPXqDcY17RxER0aemjVBs3yFpd+Ag4NOSfl67XdJLgFOBPWw/KOl8qgTRHwGH2b69GTFHRMTQNXMOZQtgue1pwNnA7sAjwMalyVjgMWCZpBcAbyrrbwe2ltRVXh9R0+1VwEmSVI7x8mbFHxERg9PMOZRdgLMlrQJWAB8E9gJmSLrH9uskzQF+C9wF3Ahg+3FJJ5R2jwE31/T578CXgHmSRgGLgTc38T1ERMQAyXa7Y1iNpI1sP1pGIv8F/M72F4fSV3d3t/MI4IiIwZE023b3YPbp1DvlP1Am3hcC46iu+oqIiA7WkSOURlpv8229+Xu+1O4w+rXkrIPbHUJExLOsTSOUiIgYYVpdy+vYcvVXz+tzJe3YyhgiIqI5Gp5Qyl3vfTmWqvwKALbfb/u2RscQERGtN6iEIqmr1Ni6WNIiSdMlbSBpSbkj/hbg7ZImS7pJ0jxJl0naRNLhQDdwsaS5ktaXdJ2k7tL3o5LOLnW9fiZpz7L9TklvKW1GlzY3l77/qeGfSEREDMlQRigvA75mewfgYeCEsv5+27vb/h5wIfAx27tSlWD5pO3pwCzgKNuTbT9e1++GVHW6dqK6AfLTwBuAtwGfKm3eByyzvQewB9XVYC+pD1DSlFJ0ctbK5b1VfImIiEYbSkK5y/aNZXkasG9Z/j78vQjkeNvXl/UXAK8ZQL9PATPK8nzgetsrynJXWf9G4JhySfFvgM2Abes7Si2viIjWG8qd8vXXGfe8fqy+4SCt8DPXMK8CngSwvUpST5wCTrJ91TCPFRERDTaUEcokSXuV5SOBG2o3lqrCD0p6dVn1bqBntFJby2sorgI+KGldAEnbSdpwGP1FRESDDGWEcjtwoqTzgNuArwMn1bV5D/CN8rjfO6megwJwfln/OFVdr8E6l+r01y2lLMt9wCH97bDLluOYlRsHIyKablB3ypcKwFfY3rlZATVaanlFRAxe7pSPiIi2GVRCsb1kJI1OAObfvYyu065sdxgREWu9ho1Qam9S7KfNKWVepWV9RUREa7T6lNcpQKOSQCP7ioiIYRp0Qumr/Epdm6+XO9UXSjqzrDuZqo7XtZKuLeveKOnXkm6RdImkjXo53mpteusrIiLaa6gjlL7Kr/T4RLk6YFfgtZJ2tf0V4B7gdeXxvxOA04EDbO9OVZbl/9Z20leb+r7qg0vplYiI1hvqM+Xry6+cXLf9HZKmlP43B3YE5tW1eVVZf2N1SwnPA349hDarsT0VmArVA7YG9pYiImI4hppQ+iq/QinWeCqwh+0HJZ0PjOmlDwE/tf2ufo4zkDYREdEBhnrKq7/yK2Op6notk/QC4E0122pLr9wE7CPppQCSNpS0Xd1x+msz3DIuERHRQENNKD3lVxYBm1CVXwHA9q3AHOC3wHeAG2v2mwrMkHSt7fuoHrj1XUnzqE5lbV97kDW0+XtfQ3wPERHRQIMqvQIjr/xKSq9ERAxeSq9ERETbDHpS3vYSYESMTiIionXW+hFKanlFRLRG2xOKpPGS6m+MrN3+q1bGExERQ9P2hAKMZ/U77el57K/tvVseUUREDNpQb2xspLOAbSTNBVYATwAPUl0evJ2kR21vJGkUcA6wP3BXaXue7eltijsiImp0QkI5DdjZ9mRJ+wFXlteL69odSvX43x2B5wOLgPN667CUfZkCMHrsxOZEHRERz9IJp7zqzewlmQDsC1xie5XtvwJ93tBoe6rtbtvdozcY17RAIyLiGZ2YUB5rdwARETF4nZBQBlqT60bgMEmjSo2w/ZoaVUREDErb51Bs3y/pRkkLgMeBe/toeinweuA2qkn5W4A87CQiokMMupZXO0nayPajkjYDZgL7lPmUPqWWV0TE4A2lllfbRyiDdIWk8VQP2vr3NSWTiIhonRGVUGzvN9h91lR6ZclZBw8npIiIKDphUj4iItYCIzKhSDpW0hbtjiMiIp4xIhMK1VMck1AiIjpIR82hSLoc2AoYA3wZ+O/y1Q2YqtTKXeX1xZIeB/ay/Xh7Io6IiB4dlVCA42w/IGl94GZgNrBlz+OGJY23/ZCkDwGn2u71euDU8oqIaL1OO+V1sqRbgZuoRirPA7aW9FVJBwIPD6ST1PKKiGi9jkkopdLwAVSnsHYD5gDrAbsB1wHHA+e2K76IiOhfJ53yGgc8aHu5pO2BVwETgFG2L5V0OzCttB1o/a+IiGiRTkooM4DjJS0Cbqc67bUlcF15uBbAx8v384FvZFI+IqJzjKhaXkORWl4REYM3lFpeHTOHEhERI1snnfJqit5qeaV+V0RE42WEEhERDTHkhCJpiaQJkrrKw7GG1c9Q94+IiM6QEUpERDTEgBKKpMslzZa0sJQ16avd1pLmSNqjbn2XpN9KuljSIknTJW1Q0+QkSbdIml/uQUHSpuW48yTdJGnXsv61kuaWrzmScj9KREQHGOgI5Tjbr6AqynhyeQTvs0h6GdVz34+1fXMvfbwM+JrtHahKqJxQs22p7d2BrwOnlnVnAnNs7wr8C3BhWX8qcKLtycCrqZ5DXx/LFEmzJM1auTyPnY+IaIWBJpT6Glvb1m2fCPwQOMr2rX30cZftG8vyNGDfmm3/U77PBrrK8r7ARQC2rwE2kzQWuBH4T0knA+NtP11/oNTyiohovTUmlD5qbI2pa7YM+BPPThL16u+grH39ZPm+kjVcymz7LOD9wPrAjT2nyCIior0GMkLprcZWvaeAtwHHSDqyj34mSdqrLB8J3LCG4/4SOAr+ntSW2n5Y0ja259v+PFWJ+ySUiIgOMJAbG3ursbUa249JejPwU0mPArOAc20fVJrcDpwo6TzgNqr5kv6cAZwnaR6wHHhPWX+KpNcBq4CFwE/662SXLccxKzcyRkQ0XUtqeUnqAq7oeVBWK6WWV0TE4KWWV0REtE1LannZXgK0fHQCvdfy6k3qe0VEDE9GKBER0RAjPqFIGt3uGCIiogEJRdIxpTzKrZIuKmVWrinrfi5pUml3vqRvlDvY7yhXhCHpWEk/lHSdpN9J+mRN30dLmlnKrHyzJ3lIelTSf5SbLffqNbCIiGipYc2hSNoJOB3Y2/ZSSZsCFwAX2L5A0nHAV4BDyi5dwJ7ANsC1kl5a1u9JNceyHLhZ0pXAY8ARwD62V0j6GtV9KRcCGwK/sf3Pw4k/IiIaZ7iT8vsDl9heCmD7gXLz4qFl+0XAF2ra/8D2KuB3ku7kmZsSf2r7fgBJ/0N1x/3TwCuoEgxUd8b/rbRfSVU3rFelgOUUgNFjJw7zLUZExEC0+omNfZVf6W29qEY6H++lnydsr+zzIPZUYCrAeptv2/wbbSIiYthzKNcAb++pPlxOef0KeGfZfhRVCZUeb5c0StI2wNZUd88DvKGUq1+f6vTYjcDPgcMlPb+nb0kvHma8ERHRJMMaodheKOkzwPWSVlIVjjwJ+LakjwD3Ae+t2eVPwExgLHC87SfK6ayZVKewXgRMsz0LQNLpwNWSRgErgBOBPw4n5oiIaI6WlF6B6iovqvIr0+vWHwt02/5QM46b0isREYOX0isREdE2LZuUt31sH+vPB85vVRwREdEcrb7Kq+UGWsurGVIfLCKeS4Z1ykvSEkkTyt3xCxoVVN0xzpW0YzP6joiIxun4EYrt97c7hoiIWLMBj1AkXS5ptqSF5U70eutIuljSIknTJW1Q9lsiaUJZ7pZ0XVk+Q9IFkn4p6Y+SDpX0BUnzJc2QtG5pd52k7rL8qKTPlLphN0l6wXA/gIiIaIzBnPI6zvYrgG7g5J6bGWu8DPia7R2Ah4ETBtDnNlTlW94CTAOutb0L8DjQ2wTEhsBNtncDfgF8oLdOJU0pRShnrVy+bABhRETEcA0moZxcqvveBGwFbFu3/S7bN5blaVT1uNbkJ7ZXAPOB0VTPr6e87uql/VPAFWV5dh9tsD3Vdrft7tEbjBtAGBERMVwDmkORtB9wALCX7eXltNWYumZ91el6mmcSV/0+TwLYXiVphZ+5y3JVH7HVtlk50PgjIqL5BjpCGQc8WJLJ9sCremkzqVQaBjgSuKEsL6GqGgxw2FADjYiIzjbQhDKDatJ9EXAW1WmvercDJ5Y2mwBfL+vPBL4saRbVqCIiItZCLavl1S6p5RURMXip5RUREW2z1k9q91d6JaVRIiIaJyOUiIhoiI5JKJLOl3T4INo3rX5YREQMXscklIiIGNnallAkHSNpXqnLdVFZ/RpJv5J0Z89oRZWzJS0odb6OaFfMERHRt7ZMykvaCTgd2Nv2UkmbAv8JbE5VsmV74EfAdOBQYDKwGzABuFnSL9bQ/xRgCsDosROb9TYiIqJGu0Yo+wOX2F4KYPuBsv5y26ts3wb0VBLeF/iu7ZW27wWuB/bor/PU8oqIaL1Om0N5smZZbYsiIiIGrV0J5Rrg7T0l8Mspr778EjhC0mhJE4HXADNbEGNERAxCW+ZQbC+U9BngekkrgTn9NL8M2Au4laqC8Udt/1VSV9MDjYiIAUstr4iIWE1qeUVERNs8p2t5RUSsrdpRqzAjlIiIaIiOSih91eeSdK6kHcvyEkkTyvKvWh1jRET0bkSc8rL9/j7W793qWCIioncdNUIp1pF0saRFkqZL2kDSdZJWu9pA0qPtCDAiIlbXiQnlZcDXbO8APAycMNgOJE2RNEvSrJXLlzU8wIiIWF0nJpS7bN9YlqdR1fIalNTyiohovU5MKPV3Wq7dd15GRKwlOjGhTJK0V1k+ErihncFERMTAdGJCuR04UdIiYBPg622OJyIiBiC1vCIiYjWp5RUREW0zIm5sHI7U8oo1aUfNo4i1UdNHKJKOlXROs48TERHtlVNeERHREMNOKJKOljRT0lxJ3yyP6n2vpDskzQT2qWk7UdKlkm4uX/uU9WdIukDSLyX9UdKhkr4gab6kGZLWLe3+rey3QNJUSXnufEREhxhWQpG0A3AEsI/tycBK4GjgTKpEsi+wY80uXwa+aHsP4DDg3Jpt2wD7A2+hukP+Wtu7AI8DPSe5z7G9h+2dgfWBNw8n/oiIaJzhTsq/HngFcHMZLKwP7A1cZ/s+AEnfB7Yr7Q8AdqwZWIyVtFFZ/ontFZLmA6OBGWX9fKCrLL9O0keBDYBNgYXA/9YHJWkKMAVg9NiJw3yLERExEMNNKAIusP3xv6+QDgEO7aP9KOBVtp94VidVgnkSwPYqSSv8zA0yq6gqEI8BvgZ0275L0hnAmN4OYnsqMBVgvc23XbtvtImI6BDDnUP5OXC4pOcDSNoUmAO8VtJmZe7j7TXtrwZO6nkhafIgjtWTPJaWUc3hw4o8IiIaalgjFNu3SToduFrSKGAFcCJwBvBr4CFgbs0uJwP/JWleOfYvgOMHeKyHJH0LWAD8Fbh5OLFHRERjpfRKRESsJqVXIiKibZJQIiKiIZJQIiKiIVpRy+t4Scf0sr5L0oJmHz8iIlqj6dWGbX+jEf1IWsf2043oKyIiGm9IIxRJG0q6UtKtpa7WEZKW1NTfminppaXtGZJOLcuvKPvcSnV5cU9/oyWdXep0zZP0T2X9fqW+14+A23o77vA/goiIaIShnvI6ELjH9m6lrlZPmZRlpf7WOcCXetnv28BJtnerW/++su8ewB7AByS9pGzbHfg/trfr57jPImmKpFmSZt13331DfIsRETEYQ00o84E3SPq8pFfbXlbWf7fm+161O0gaD4y3/Yuy6qKazW8EjpE0F/gNsBmwbdk20/biNRz3WWxPtd1tu3vixNTyiohohSElFNt3UI0c5gOflvRvPZtqmw2iS1GNXCaXr5fYvrpse2wAx42IiDYb6hzKFsBy29OAs6l+yUNVyr7n+69r97H9EPCQpH3LqqNqNl8FfLDmuSfbSdpwEMeNiIg2G+pVXrsAZ0taRVW/64PAdGCTUqfrSeBdvez3XuA8SaYqFNnjXKoS9beUh2bdBxwywONGREQHaFgtL0lLqErLL21Ihw2SWl4REYOXWl4REdE2Dbux0XZXo/pqpPl393ohWERENFhGKBER0RDDSijNrsdV+j+yWf1HRETjdPoIpQtIQomIGAEakVBGS/qWpIWSrpa0vqRtJM2QNLvU4toeQNL5kr5RyqLcIenNZX1XaXdL+dq79H0W8GpJcyV9WNIYSd8u9cLmSHpdA+KPiIgGaMSk/LbAu2x/QNIPgMOo7jc53vbvJL0S+Bqwf2nfBewJbANcW4pI/g14g+0nJG1LVbqlGzgNONV2T+L5Z8C2dylJ6mpJ29l+ojYgSVOAKQCjx6b0SkREKzQioSy2Pbcsz6ZKGHsDl1T3KAKwXk37H9heBfxO0p3A9sBi4BxJk4GVwHZ9HGtf4KsAtn8r6Y+l7bzaRranAlMB1tt828bcaBMREf1qREJ5smZ5JfAC4CHbk/toX/8L3sCHgXuB3ahOwz1Rv1NERHS2ZkzKPwwslvR2AFVqy9W/XdIoSdsAWwO3A+OAv5SRy7uB0aXtI8DGNfv+klIDTNJ2wKSyf0REtFmzrvI6CnhfeZDWQuCtNdv+BMwEfkI1z/IE1RzLe0r77XmmwvA8YGV5oNaHS7tRkuYD3weOtV07QoqIiDZpWC2vAR1MOh+4wvb0Vh0ztbwiIgYvtbwiIqJtGlbLayBsH9vK40FVy6vrtCsH3H7JWQc3MZqIiLVXRigREdEQbU0okraQNL0sT5Z0UDvjiYiIoWtrQrF9j+3Dy8vJQBJKRMQINeyEIunyUrNroaQpkkaXml0LSs2tD5d2L5X0s3IJ8C2l3ldXafc84FPAEaVu1xGSNpR0nqSZpW7XW0s/O5V1cyXNK6VaIiKizRoxKX+c7QckrQ/cTFV+ZUvbOwNIGl/aXQycZfsySWOoktnzAWw/JenfqB4h/KGy32eBa2wfV/qYKelnwPHAl21fXBLRaOqklldEROs14pTXyeWGxJuArYDnAVtL+qqkA4GHJW1MlWQuA7D9hO3la+j3jcBpkuYC1wFjqO6M/zXwL5I+BrzY9uP1O9qearvbdvfoDcY14C1GRMSaDPcBW/sBBwB72d4NmENVCHI3qiRwPHDuULsHDrM9uXxNsr3I9neAtwCPAz+WtH//3URERCsMd4QyDnjQ9vJSTv5VwARglO1LgdOB3W0/AvxZ0iEAktaTtEFdX/V1u64CTlIpWSzp5eX71sCdtr8C/BDYdZjvISIiGmC4CWUGsI6kRVQPw7oJ2BK4rpyqmgZ8vLR9N9XpsXnAr4AX1vV1LbBjz6Q88O/AusA8SQvLa4B3AAtK/zsDFw7zPURERAO0tJZXO6SWV0TE4KWWV0REtM1an1Dm372s3SFERDwnrPUJJSIiWiMJJSIiGiIJJSIiGqLtCaXU7Lqy1PhaUOp4vULS9aVG2FWSNpe0jqSby82USPqcpM+0OfyIiCha+oCtPhwI3GP7YABJ46ieN/9W2/eVe1I+U2p6HQtMl3RS2e+VvXWYWl4REa3XCQllPvAfkj4PXAE8SHXD4k/LTfKjgb8A2F4o6aLSbi/bT/XWoe2pwFSA9Tbfdu2+0SYiokO0PaHYvkPS7lTPQvk0cA2w0PZefeyyC/AQpVJxRER0hk6YQ9kCWG57GnA21WmsiZL2KtvXlbRTWT4U2BR4DfDVmtL4ERHRZm0foVCNOM6WtApYAXwQeBr4SplPWQf4kqR7qeqFvd72XZLOAb4MvKffzrdM+fqIiFZoe0KxfRVVZeF6r+ll3XY1+32laUFFRMSgtf2UV0RErB3aPkJptvl3L6PrtCub1v+Ssw5uWt8RESNJRigREdEQSSgREdEQbU8oki4vJVYWljvckfQ+SXdIminpW+WKLiRNlHRpKcFys6R92ht9RET06IQ5lONsPyBpfeBmSVcC/wrsTvWc+WuAW0vbLwNftH2DpElUV4ftUN9hSq9ERLReJySUkyW9rSxvRfXs+ettPwAg6RKeuVz4AKrnzvfsO1bSRrYfre0wpVciIlqvrQmlVA4+gKou13JJ1wG/pZdRRzEKeJXtJ1oTYUREDFS751DGAQ+WZLI98CpgQ+C1kjaRtA5wWE37q4GTel5ImtzSaCMiok/tTigzgHUkLaIqq3ITcDfwWWAmcCOwBOh5MPzJQLekeZJuA45vecQREdEr2Z03xdAzL1JGKJcB59m+bCh9dXd3e9asWY0NMCJiLSdptu3uwezT7hFKX86QNBdYACwGLm9zPBERsQadcJXXamyf2u4YIiJicDp1hBIRESNMRyYUSeMlnTDYbRER0T4dmVCA8cBqSaNM0ve6LSIi2qsj51CoLiHepkzMrwCeAB4Etgduqdn2U9sfaV+YERHRo1MTymnAzrYnl7vpryyvF0vq6tnW1861tbwmTZrU/GgjIqJjT3nVm2l78dNZZ3IAAAU9SURBVEAb255qu9t298SJKQ4ZEdEKIyWhPNbuACIion+dmlAeATYewraIiGiTjpxDsX2/pBslLQAeB+7tY9tPMikfEdEZOjKhANg+cijbIiKiPTr1lFdERIwwSSgREdEQSSgREdEQSSgREdEQSSgREdEQSSgREdEQSSgREdEQSSgREdEQSSgREdEQst3uGJpK0iPA7e2OYw0mAEvbHcQAjIQ4E2NjjIQYYWTEOVJjfLHtQZVr79jSKw10u+3udgfRH0mzOj1GGBlxJsbGGAkxwsiI87kUY055RUREQyShREREQzwXEsrUdgcwACMhRhgZcSbGxhgJMcLIiPM5E+NaPykfERGt8VwYoURERAskoUREREOM6IQi6UBJt0v6vaTTetm+nqTvl+2/kdRVs+3jZf3tkv6h02KUtJmkayU9KumcZsU3zBjfIGm2pPnl+/4dGOOekuaWr1slva1ZMQ4nzprtk8q/+amdFqOkLkmP13ye3+i0GMu2XSX9WtLC8rM5ppNilHRUzWc4V9IqSZObEeMw41xX0gXlM1wk6eNrPJjtEfkFjAb+AGwNPA+4Fdixrs0JwDfK8juB75flHUv79YCXlH5Gd1iMGwL7AscD53To5/hyYIuyvDNwdwfGuAGwTlneHPhbz+tOirNm+3TgEuDUTosR6AIWNOtnsUExrgPMA3YrrzfrtP/bdW12Af7QoZ/lkcD3yvIGwBKgq7/jjeQRyp7A723fafsp4HvAW+vavBW4oCxPB14vSWX992w/aXsx8PvSX8fEaPsx2zcATzQhrkbFOMf2PWX9QmB9Set1WIzLbT9d1o8BmnkVynB+JpF0CLCY6rPsyBhbZDgxvhGYZ/tWANv3217ZYTHWelfZt1mGE6eBDSWtA6wPPAU83N/BRnJC2RK4q+b1n8u6XtuUXyrLqP5iGci+7Y6xVRoV42HALbaf7LQYJb1S0kJgPnB8TYLpmDglbQR8DDizSbENO8ay7SWS5ki6XtKrOzDG7QBLukrSLZI+2oEx1joC+G6TYnxWDMVg4pwOPAb8BfgT8P9sP9DfwZ4LpVeiySTtBHye6q/DjmP7N8BOknYALpD0E9vNHvkN1hnAF20/2trBwKD8BZhk+35JrwAul7ST7X7/am2xdahOFe8BLAd+Lmm27Z+3N6zVSXolsNz2gnbH0oc9gZXAFsAmwC8l/cz2nX3tMJJHKHcDW9W8flFZ12ubMmwbB9w/wH3bHWOrDCtGSS8CLgOOsf2HToyxh+1FwKNU8z2dFucrgS9IWgKcAvyLpA91UozlFPH9ALZnU52b366TYqT6C/wXtpfaXg78GNi9w2Ls8U6aOzp5VgzFYOI8Ephhe4XtvwE3Av3X+2rWZFCzv6j+ErmTalK9Z7Jpp7o2J/LsyaYflOWdePak/J00Z+JuyDHWbD+W5k7KD+dzHF/aH9rB/9Yv4ZlJ+RcD9wATOi3OujZn0LxJ+eF8lhN7/p9QTfLeDWzaYTFuAtxCuRgD+BlwcCfFWF6PKp/f1s34d27QZ/kx4NtleUPgNmDXfo/XzDfT7C/gIOAOqr+UPlHWfQp4S1keQ3XFzO+BmbX/eMAnyn63A2/q0BiXAA9Q/VX9Z+quzmh3jMDpVOdY59Z8Pb/DYnw31ST33PKL5pBO/Zms6eMMmpRQhvlZHlb3Wf5jp8VYth1d4lwAfKFDY9wPuKmZP4sN+PfeqKxfSJVMPrKmY6X0SkRENMRInkOJiIgOkoQSERENkYQSERENkYQSERENkYQSERENkYQSERENkYQSEREN8f8BxC7X3bcCzPgAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"vimps = np.array(forest.vimp(model).rx('importance')[0])\n",
"\n",
"y = np.arange(len(vimps))\n",
"plt.barh(y, np.abs(vimps))\n",
"plt.yticks(y, df_train.drop(['time', 'status'], axis=1).columns)\n",
"plt.title(\"VIMP (absolute value)\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"colab_type": "text",
"id": "2YGhK2xwjkiA"
},
"source": [
"### Question:\n",
"\n",
"How does the variable importance compare to that of the Cox model? Which variable is important in both models? Which variable is important in the random survival forest but not in the Cox model? You should see that `edema` is important in both the random survival forest and the Cox model. You should also see that `bili` is important in the random survival forest but not the Cox model ."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Congratulations!\n",
"\n",
"You've finished the last assignment in course 2! Take a minute to look back at the analysis you've done over the last four assignments. You've done a great job!"
]
}
],
"metadata": {
"accelerator": "GPU",
"colab": {
"collapsed_sections": [],
"include_colab_link": true,
"name": "C2M4_Assignment.ipynb",
"provenance": []
},
"coursera": {
"schema_names": [
"AI4MC2-4"
]
},
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}