{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Kaplan-Meier Estimation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Run this notebook on Colab](https://colab.research.google.com/github/AllenDowney/SurvivalAnalysisPython/blob/master/02_kaplan_meier.ipynb)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook introduces Kaplan-Meier estimation, a way to estimate a hazard function when the dataset includes both complete and incomplete cases.\n", "To demonstrate, I'll use a small set of hypothetical data. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Dog adoption data\n", "\n", "Suppose you are investigating the time it takes for dogs to get adopted from a shelter. You visit a shelter every week for 10 weeks, and record the arrival time for each dog and the adoption time for each dog that was adopted.\n", "\n", "Here's what the data might look like. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
startendstatus
0051
1121
2261
3290
4490
5681
6790
\n", "
" ], "text/plain": [ " start end status\n", "0 0 5 1\n", "1 1 2 1\n", "2 2 6 1\n", "3 2 9 0\n", "4 4 9 0\n", "5 6 8 1\n", "6 7 9 0" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obs = pd.DataFrame()\n", "\n", "obs['start'] = 0,1,2,2,4,6,7\n", "obs['end'] = 5,2,6,9,9,8,9\n", "obs['status'] = 1,1,1,0,0,1,0\n", "\n", "obs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This `DataFrame` contains one row for each dog and three columns:\n", "\n", "* `start`: arrival time, in weeks since the beginning of the study\n", "\n", "* `end`: adoption date, for dogs that were adopted, or `9` for dogs that had not been adopted at the end of the study\n", "\n", "* `status`: `1` for dogs that were adopted; `0` for dogs that were not." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting lifelines\n", "\n", "The following function visualizes the data." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def plot_lifelines(obs):\n", " \"\"\"Plot a line for each observation.\n", " \n", " obs: DataFrame\n", " \"\"\"\n", " for y, row in obs.iterrows():\n", " start = row['start']\n", " end = row['end']\n", " status = row['status']\n", " \n", " if status == 0:\n", " # ongoing\n", " plt.hlines(y, start, end, color='C0')\n", " else:\n", " # complete\n", " plt.hlines(y, start, end, color='C1')\n", " plt.plot(end, y, marker='o', color='C1')\n", " \n", " plt.xlabel('Time (weeks)')\n", " plt.ylabel('Dog index')\n", " plt.gca().invert_yaxis()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here are the results:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAATf0lEQVR4nO3df7BndX3f8efLZZWV8CMVdCyLrlrF2s0iejXSTRwBY6gEnURdtRgrmclqJwoYMxrtFBk66UwSEyKpOm4VNGrVLUojYhBqNGtay3AXcfmxkFCCcoOWS1sXTHdggXf/+J4b7m7u7n65e889937u8zFz5/s9n3Pu+bz3O7uv+ezne87npKqQJLXnCUMXIEnqhwEvSY0y4CWpUQa8JDXKgJekRh02dAGzHXvssbVu3bqhy5CkZWP79u33VdVxc+1bUgG/bt06Jicnhy5DkpaNJN/f3z6naCSpUQa8JDXKgJekRhnwktQoA16SGtVrwCc5I8ntSe5I8tt99iUtWTu2wsXr4cJjRq87tg5dkVaI3i6TTLIK+AjwC8AUcH2Sr1TVrX31KS05O7bClefCnt2j7V13j7YBNmwari6tCH1eB/9S4I6quhMgyReA1wL9BPxlZ/ZyWumQTF0Pjzy4d9ue3fCNiwx49a7PKZrjgbtnbU91bXtJsjnJZJLJ6enpHsuRBrBvuM/YNbW4dWhF6nMEnzna/sHTRapqC7AFYGJiYv5PHznnqnn/qtSbi9ePpmX2dfTaxa9FK06fI/gp4IRZ22uBe3rsT1p6Tr8AVq/Zu231mlG71LM+A/564LlJnpXkicCbgK/02J+09GzYBGddAkefAGT0etYlzr9rUfQ2RVNVDyd5J/B1YBVwaVXd0ld/0pK1YZOBrkH0uppkVX0N+FqffUiS5uadrJLUKANekhplwEtSowx4SWqUAS9JjTLgJalRBrwkNcqAl6RGGfCS1CgDXpIaZcBLUqMMeElqlAEvSY0y4CWpUQa8JDXKgJekRhnwktQoA16SGmXAS1KjDHhJapQBL0mNMuAlqVEGvCQ1yoCXpEb1FvBJLk1yb5Kb++pD+9ixFS5eDxceM3rdsXXoiiQNqM8R/KeAM3o8v2bbsRWuPBd23Q3U6PXKcw15aQU7rK8TV9W2JOv6Ov+SddmZw/Q7dT088uDebXt2wzcugg2bhqlJ0qAGn4NPsjnJZJLJ6enpoctZvvYN9xm7pha3DklLRm8j+HFV1RZgC8DExEQNXM6hO+eqYfq9eH03PbOPo9cufi2SloTBR/BaIKdfAKvX7N22es2oXdKKZMC3YsMmOOsSOPoEIKPXsy5x/l1awXqboknyeeAVwLFJpoAPVtUn++pPjMLcQJfU6fMqmjf3dW5J0sE5RSNJjTLgJalRBrwkNcqAl6RGGfCS1CgDXpIaZcBLUqMMeElqlAEvSY0y4CWpUQa8JDXKgJekRhnwktQoA16SGmXAS1KjDHhJapQBL0mNMuAlqVEGvCQ1yoCXpEYZ8JLUKANekhplwEtSowx4SWpUbwGf5IQk30yyM8ktSc7rqy9JY9ixFS5eDxceM3rdsXXoitSzw3o898PAe6rqhiRHAtuTXFtVt/bYp6S57NgKV54Le3aPtnfdPdoG2LBpuLrUq94Cvqp+CPywe/9Akp3A8YAB37fLzhy6Ai01U9fDIw/u3bZnN3zjIgO+YYsyB59kHXAycN0c+zYnmUwyOT09vRjlSCvPvuE+Y9fU4tahRdXnFA0ASX4K+BJwflXdv+/+qtoCbAGYmJiovutZEc65augKtNRcvH40LbOvo9cufi1aNL2O4JOsZhTun6uqL/fZl6QDOP0CWL1m77bVa0btalafV9EE+CSws6r+sK9+JI1hwyY46xI4+gQgo9ezLnH+vXF9TtFsBH4VuCnJjV3bB6rqaz32KWl/Nmwy0FeYPq+i+UsgfZ1fknRg3skqSY0y4CWpUQa8JDXqoAGf5KlztJ3YTzmSpIUyzgj+20n+/qv3JO8BruivJEnSQhjnKppXAFuSvAF4GrATeGmfRUmSDt1BR/DdomFXA6cA64A/qaqf9FyXJOkQHXQEn+RaRqtCrgfWApcm2VZVv9V3cZKk+RtnDv4jVfXWqvpxVd3MaCS/q+e6JEmHaJwpmv+S5OeSnNM1/TTw2X7LkiQdqnEuk/wg8D7g/V3TEzHgJWnJG2eK5peB1wB/B1BV9wBH9lmUJOnQjRPwD1VVAQWQ5Ih+S5IkLYRxAn5rko8DxyT5deC/Av+x37IkSYfqoJdJVtWHkvwCcD9wInBBVV3be2WSpEMy1nrwXaAb6pK0jOw34JM8QDfvPpeqOqqXiiRJC2K/AV9VRwIkuQj4EfAZRk9oOhuvopGkJW+cL1l/sao+WlUPVNX9VfUx4HV9FyZJOjTjBPwjSc5OsirJE5KcDTzSd2GSpEMzTsD/S2AT8L+6nzd0bZKkJWycyyTvAl7bfymSpIU0znLBxwG/zmgt+L8/vqp+rb+yJEmHapzr4P8U+DajO1ide5ekZWKcgH9yVb3v8Z44yeHANuBJXT+XV9UHH+95JEnzM86XrF9N8up5nPtB4LSqOgl4IXBGkpfN4zySpHkYZwR/HvCBJA8Cexjd7FQHu5O1W4Fy5tmtq7uf/d4Zq4Xzxo9/Z+gSJD0OX3z7Kb2cd5wnOh1ZVU+oqjVVdVS3PdYyBd218zcC9wLXVtV1cxyzOclkksnp6enH/QeQJM0to4H2HDuS51fVbUleNNf+qrph7E6SY4ArgHd1z3Wd08TERE1OTo57Wkla8ZJsr6qJufYdaIrmN4HNwB/Msa+A08YtoKp+nORbwBnAfgNekrRwDrTY2Obu9dT5nLi7fn5PF+5rgFcCvzuvKiVJj9tY68HP09OBTydZxWiuf2tVfbXH/iRJs/QW8FW1Azi5r/NLkg5snOvgJUnL0Dhr0cx1Fc0u4PtV9fDClyRJWgjjTNF8FHgRsIPRTU7ru/dPSfKOqrqmx/okSfM0zhTNXcDJVTVRVS9mNK9+M6OrYn6vx9okSYdgnIB/flXdMrNRVbcyCvw7+ytLknSoxpmiuT3Jx4AvdNtvBP4qyZMYrU0jSVqCxhnBvw24AzgfeDdwZ9e2B5jXTVCSpP6N88i+3Un+GLiG0RIFt1fVzMj9J/v/TUnSkMa5TPIVwKcZfdka4IQk/6qqtvVamSTpkIwzB/8HwKuq6naAJM8DPg+8uM/CJEmHZpw5+NUz4Q5QVX/F6OEdkqQlbJwR/GSSTwKf6bbPBrb3V5IkaSGME/D/GvgN4FxGc/DbGN3dKklawsa5iubBJJ8BPlNVPlNPkpaJ/c7BZ+TCJPcBtzG64Wk6yQWLV54kab4O9CXr+cBG4CVV9ZSq+kfAzwIbk7x7MYqTJM3fgQL+rcCbq+pvZhq69Wfe0u2TJC1hBwr41VV1376N3Ty8l0lK0hJ3oIB/aJ77JElLwIGuojkpyf1ztAc4vKd6JEkLZL8BX1WrFrMQSdLC8qHbktQoA16SGmXAS1Kjeg/4JKuSfDfJV/vuS5L0mMUYwZ8H7FyEfiRJs4yzmuS8JVkLnAn8DvCbffYl7c8bP/6doUuQDuiLbz+ll/P2PYL/I+C9wKP7OyDJ5iSTSSanp12sUpIWSm8j+CS/BNxbVdu757rOqaq2AFsAJiYmqq96tHL1NTqSlro+R/AbgdckuQv4AnBaks/22J8kaZbeAr6q3l9Va6tqHfAm4M+r6i199SdJ2pvXwUtSo3q9imZGVX0L+NZi9CVJGnEEL0mNMuAlqVEGvCQ1yoCXpEYZ8JLUKANekhplwEtSowx4SWqUAS9JjTLgJalRBrwkNcqAl6RGGfCS1CgDXpIaZcBLUqMMeElqlAEvSY0y4CWpUQa8JDXKgJekRhnwktQoA16SGmXAS1KjDHhJatRhfZ48yV3AA8AjwMNVNdFnf5I0lh1b4RsXwa4pOHotnH4BbNg0dFULrteA75xaVfctQj+SdHA7tsKV58Ke3aPtXXePtqG5kF+MgJe0lFx25tAVDGvqenjkwb3b9uwejegbC/i+5+ALuCbJ9iSb5zogyeYkk0kmp6eney5H0oq3b7jP2DW1uHUsgr5H8Bur6p4kTwWuTXJbVW2bfUBVbQG2AExMTFTP9Ug656qhKxjWxetH0zL7Onrt4tfSs15H8FV1T/d6L3AF8NI++5Okgzr9Ali9Zu+21WtG7Y3pLeCTHJHkyJn3wKuAm/vqT5LGsmETnHUJHH0CkNHrWZc0N/8O/U7RPA24IslMP/+pqq7usT9JGs+GTU0G+r56C/iquhM4qa/zS5IOzDtZJalRBrwkNcqAl6RGGfCS1CgDXpIaZcBLUqMMeElqlAEvSY0y4CWpUQa8JDXKgJekRhnwktQoA16SGmXAS1KjDHhJapQBL0mNMuAlqVEGvCQ1yoCXpEYZ8JLUKANekhplwEtSowx4SWqUAS9Jjeo14JMck+TyJLcl2ZnklD77kyQ95rCez/9h4Oqqen2SJwJP7rk/SVKnt4BPchTwcuBtAFX1EPBQX/1JWl7e+PHvDF3CkvHFt/czudHnFM2zgWngsiTfTfKJJEfse1CSzUkmk0xOT0/3WI4krSypqn5OnEwA/wPYWFXXJfkwcH9V/dv9/c7ExERNTk72Uo8ktSjJ9qqamGtfnyP4KWCqqq7rti8HXtRjf5KkWXoL+Kr6EXB3khO7ptOBW/vqT5K0t76vonkX8LnuCpo7gXN67k+S1Ok14KvqRmDOuSFJUr+8k1WSGmXAS1KjDHhJapQBL0mN6u1Gp/lIMg18f56/fixw3wKWs5z5WezNz2Nvfh6PaeGzeGZVHTfXjiUV8IciyeT+7uZaafws9ubnsTc/j8e0/lk4RSNJjTLgJalRLQX8lqELWEL8LPbm57E3P4/HNP1ZNDMHL0naW0sjeEnSLAa8JDVq2Qd8kjOS3J7kjiS/PXQ9Q0pyQpJvdg84vyXJeUPXNLQkq7onin116FqGluSYJJcnua37O9LPc+KWiSTv7v6d3Jzk80kOH7qmhbasAz7JKuAjwL8AXgC8OckLhq1qUA8D76mqfwq8DPiNFf55AJwH7By6iCXiw8DVVfV84CRW8OeS5HjgXGCiqtYDq4A3DVvVwlvWAQ+8FLijqu7sHur9BeC1A9c0mKr6YVXd0L1/gNE/4OOHrWo4SdYCZwKfGLqWoSU5Cng58EmAqnqoqn48aFHDOwxYk+Qw4MnAPQPXs+CWe8AfD9w9a3uKFRxosyVZB5wMXHeQQ1v2R8B7gUcHrmMpeDYwDVzWTVl9IskRQxc1lKr6W+BDwA+AHwK7quqaYataeMs94DNH24q/7jPJTwFfAs6vqvuHrmcISX4JuLeqtg9dyxJxGKNnIn+sqk4G/g5Ysd9ZJflpRv/bfxbwj4Ejkrxl2KoW3nIP+CnghFnba2nwv1mPR5LVjML9c1X15aHrGdBG4DVJ7mI0dXdaks8OW9KgpoCpqpr5H93ljAJ/pXol8DdVNV1Ve4AvA/984JoW3HIP+OuB5yZ5Vvfc1zcBXxm4psEkCaM51p1V9YdD1zOkqnp/Va2tqnWM/l78eVU1N0IbV1X9CLg7yYld0+nArQOWNLQfAC9L8uTu383pNPilc98P3e5VVT2c5J3A1xl9C35pVd0ycFlD2gj8KnBTkhu7tg9U1deGK0lLyLuAz3WDoTuBcwauZzBVdV2Sy4EbGF199l0aXLbApQokqVHLfYpGkrQfBrwkNcqAl6RGGfCS1CgDXpIaZcBryUvylCQ3dj8/SvK33fufJPloT32en+StPZ37bUn+w5jH/kyST/VRh9q3rK+D18pQVf8beCFAkguBn1TVh/rqr1t86tdYAnd6VtVNSdYmeUZV/WDoerS8OILXspXkFTPrvCe5MMmnk1yT5K4kv5Lk95LclOTqbgkHkrw4yV8k2Z7k60mePsepTwNu6G6ke2qS7d3vnpSkkjyj2/6f3Z2QxyX5UpLru5+N3f4jklzatX03yT9Y6TTJmUm+k+TYJG/o1ib/XpJtsw67kgaXslX/DHi15DmMlgd+LfBZ4JtV9TPAbuDMLuT/GHh9Vb0YuBT4nTnOsxHYDlBV9wKHd8vt/jwwCfx8kmcyWszs/zFaZ/3iqnoJ8DoeW5743zBaIuElwKnA789ewTHJLzNa8OvVVXUfcAHwi1V1EvCaWfVMdn1Lj4tTNGrJn1XVniQ3MVq64uqu/SZgHXAisB64drT8CKsYLRW7r6ez97ok/51R6L8c+PfAGYxWMv12t/+VwAu6cwIcleRI4FWMFjz7ra79cOAZ3ftTgQngVbNW/PxvwKeSbGW0+NWMexmteCg9Lga8WvIgQFU9mmRPPbYOx6OM/q4HuKWqDvaout2MwnjGtxmNoJ8J/CnwPkbLUs88BvAJwClVtXv2SbpFrF5XVbfv0/6zjNaCeTbwPEYjdKrqHd2+M4Ebk7yw+/7h8K4m6XFxikYrye3AcTPPIk2yOsk/m+O4ncA/mbW9DXgL8NdV9Sjwf4BXMxpxA1wDvHPm4CQv7N5+HXhXF/QkOXnWOb8P/ArwJzM1JHlOVV1XVRcA9/HYUtjPA26e159YK5oBrxWje6zj64HfTfI94EbmXgP8zxhNx8z83l3d25kvPv8S+HFV/d9u+1xgIsmOJLcC7+ja/x2wGtiR5OZue3Y9twNnA/85yXMYzdHf1B27Dfhed+ipwFXz+kNrRXM1SWkOSa4A3ltVfz1wHU8C/gL4uap6eMhatPwY8NIcugdjPK2qth304H7reC5wfFV9a8g6tDwZ8JLUKOfgJalRBrwkNcqAl6RGGfCS1CgDXpIa9f8BpGIEi1s0kOYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plot_lifelines(obs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each line represents the time a dog spends at the shelter. Each dot represents an adoption.\n", "We can see, for example:\n", "\n", "* The dog with index 0 arrived during week 0, and was adopted during week 5.\n", "\n", "* The dog with index 3 arrived during week 2, and had not been adopted at the end of week 9.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimating survival\n", "\n", "Now suppose we want to know the distribution of \"survival time\" from arrival to adoption.\n", "For the dogs that were adopted, we have all the data we need. \n", "For the others, we have only partial information: if a dog hasn't been adopted yet, we don't know when it will be, but we can put a lower bound on it.\n", "\n", "When we have a mixture of complete and incomplete observations -- adopted and unadopted dogs -- we can't compute the Survival function directly.\n", "Instead, we have to work backwards: we estimate the hazard function first, then use it to compute the survival function, CDF, and PMF.\n", "\n", "Specifically, we'll use Kaplan-Meier estimation, which is based on two key ideas.\n", "\n", "The first idea is that we can ignore the arrival time in the observed data, and consider only the durations. In effect, we can take the actual lifelines and shift them so they all start at 0, like this:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "tags": [ "fill-in" ] }, "outputs": [], "source": [ "duration = obs['end'] - obs['start']" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "tags": [ "fill-in" ] }, "outputs": [], "source": [ "shifted = obs.copy()\n", "shifted['start'] = 0\n", "shifted['end'] = duration" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "tags": [ "fill-in" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEGCAYAAABvtY4XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVWElEQVR4nO3df7RdZX3n8feHECVSSGZEXYyhE239MZ1MBHvBoWktP6qlgxFn1UYZrCOdGqdLRWu7rDpr0GFW23E6SmWs1siPWsQfGSytKFWplQm6lOEGmfAj6GIoDhlQLmvGAF0RQvzOH3tfuYn33pzc3H3POZv3a62sc/dzznmeL+GuT57znL2fnapCktQ/hw27AElSNwx4SeopA16SesqAl6SeMuAlqacOH3YBMx1zzDG1Zs2aYZchSWNj27ZtD1TV02Z7bqQCfs2aNUxOTg67DEkaG0m+M9dzLtFIUk8Z8JLUUwa8JPWUAS9JPWXAS1JPdRrwSc5I8q0kdyZ5R5djSerI9i1w4Vp4z6rmcfuWYVekAXV2mmSSZcCfAC8BdgI3JvlsVd3e1ZiSFtn2LXD1ebBnd3O8657mGGDdxuHVpYF0eR78ScCdVXUXQJJPAWcB3QT8ZWd20q30hLbzRtj7yL5te3bDly8w4MdAl0s0zwTumXG8s23bR5JNSSaTTE5NTXVYjqSDtn+4T9u1c2nr0IJ0OYPPLG0/dneRqtoMbAaYmJhY+N1Hzv38gt8qaQ4Xrm2WZfa3cvXS16KD1uUMfidw3Izj1cC9HY4nabGdfj4sX7Fv2/IVTbtGXpcBfyPwnCTPSvIk4NXAZzscT9JiW7cRNlwEK48D0jxuuMj19zHR2RJNVT2W5E3AF4FlwKVVdVtX40nqyLqNBvqY6nQ3yaq6BrimyzEkSbPzSlZJ6ikDXpJ6yoCXpJ4y4CWppwx4SeopA16SesqAl6SeMuAlqacMeEnqKQNeknrKgJeknjLgJamnDHhJ6ikDXpJ6yoCXpJ4y4CWppwx4SeopA16SesqAl6SeMuAlqacMeEnqKQNeknrKgJeknjLgJamnOgv4JJcmuT/JrV2NMXa2b4EL18J7VjWP27cMuyJJPdblDP7PgDM67H+8bN8CV58Hu+4Bqnm8+jxDXlJnDu+q46rammRNV/3/mMvOXLKhFmTnjbD3kX3b9uyGL18A6zYOpyZJvTb0Nfgkm5JMJpmcmpoadjnd2T/cp+3aubR1SHrC6GwGP6iq2gxsBpiYmKgFd3Tu5xerpG5cuLZdntnPytVLX4ukJ4Shz+CfME4/H5av2Ldt+YqmXZI6YMAvlXUbYcNFsPI4IM3jhotcf5fUmc6WaJJ8EjgFOCbJTuDdVXVJV+ONhXUbDXRJS6bLs2jO7qpvSdKBuUQjST1lwEtSTxnwktRTBrwk9ZQBL0k9ZcBLUk8Z8JLUUwa8JPWUAS9JPWXAS1JPGfCS1FMGvCT1lAEvST1lwEtSTxnwktRTBrwk9ZQBL0k9ZcBLUk8Z8JLUUwa8JPWUAS9JPWXAS1JPGfCS1FMGvCT1VGcBn+S4JF9JsiPJbUne0tVY0tjZvgUuXAvvWdU8bt8y7IrUQ4d32PdjwO9U1U1JjgK2Jbm2qm7vcExp9G3fAlefB3t2N8e77mmOAdZtHF5d6p3OAr6q7gPua39+KMkO4JlANwF/2ZmddCstup03wt5H9m3bsxu+fIEBr0W1JGvwSdYAJwA3zPLcpiSTSSanpqaWohxpuPYP92m7di5tHeq9LpdoAEjyE8BngLdW1YP7P19Vm4HNABMTE7Xggc79/ILfKi2pC9c2yzL7W7l66WtRr3U6g0+ynCbcr6iqv+hyLGlsnH4+LF+xb9vyFU27tIi6PIsmwCXAjqp6f1fjSGNn3UbYcBGsPA5I87jhItfftei6XKJZD/w6cEuSm9u2d1XVNR2OKY2HdRsNdHWuy7Novgqkq/4lSfPzSlZJ6ikDXpJ6yoCXpJ46YMAnefosbc/rphxJ0mIZZAZ/fZIffd2f5HeAq7orSZK0GAY5i+YUYHOSXwOeAewATuqyKEnSoTvgDL7dNOwLwMnAGuDPq+rhjuuSJB2iA87gk1xLsyvkWmA1cGmSrVX1u10XJ0lauEHW4P+kql5bVd+vqltpZvK7Oq5LknSIBlmi+cskP5/k3LbpHwAf77YsSdKhGuQ0yXcDvwe8s216Ega8JI28QZZo/iXwcuDvAarqXuCoLouSJB26QQL+0aoqoACSHNltSZKkxTBIwG9J8hFgVZLXA38DfLTbsiRJh+qAp0lW1X9J8hLgQeB5wPlVdW3nlUmSDslA+8G3gW6oS9IYmTPgkzxEu+4+m6o6upOKJEmLYs6Ar6qjAJJcAHwXuJzmDk3n4Fk0kjTyBvmS9Zer6kNV9VBVPVhVHwZ+tevCJEmHZpCA35vknCTLkhyW5Bxgb9eFSZIOzSAB/6+AjcD32j+/1rZJkkbYIKdJ3g2c1X0pkqTFNMh2wU8DXk+zF/yPXl9Vv9FdWZKkQzXIefB/BVxPcwWra++SNCYGCfinVNXvHWzHSY4AtgJPbse5sqrefbD9SJIWZpAvWT+X5F8soO9HgNOq6gXA8cAZSf75AvqRJC3AIDP4twDvSvIIsIfmYqc60JWs7Q6U0/duXd7+mfPK2EP1qo98vauuJalTn37DyZ30O8gdnY6qqsOqakVVHd0eD7RNQXvu/M3A/cC1VXXDLK/ZlGQyyeTU1NRB/wdIkmaXZqI9yxPJ86vqjiQvnO35qrpp4EGSVcBVwJvb+7rOamJioiYnJwftVpKe8JJsq6qJ2Z6bb4nmbcAm4H2zPFfAaYMWUFXfT3IdcAYwZ8BLkhbPfJuNbWofT11Ix+3583vacF8B/BLw3gVVKUk6aAPtB79AxwIfS7KMZq1/S1V9rsPxJEkzdBbwVbUdOKGr/iVJ8xvkPHhJ0hgaZC+a2c6i2QV8p6oeW/ySJEmLYZAlmg8BLwS201zktLb9+alJ/m1VfanD+iRJCzTIEs3dwAlVNVFVP0uzrn4rzVkx/7nD2iRJh2CQgH9+Vd02fVBVt9ME/l3dlSVJOlSDLNF8K8mHgU+1x68Cvp3kyTR700iSRtAgM/jXAXcCbwV+G7irbdsDLOgiKElS9wa5Zd/uJP8V+BLNFgXfqqrpmfvDc79TkjRMg5wmeQrwMZovWwMcl+RfV9XWTiuTJB2SQdbg3we8tKq+BZDkucAngZ/tsjBJ0qEZZA1++XS4A1TVt2lu3iFJGmGDzOAnk1wCXN4enwNs664kSdJiGCTgfwt4I3AezRr8VpqrWyVJI2yQs2geSXI5cHlVeU89SRoTc67Bp/GeJA8Ad9Bc8DSV5PylK0+StFDzfcn6VmA9cGJVPbWq/iHwImB9kt9eiuIkSQs3X8C/Fji7qv5uuqHdf+Y17XOSpBE2X8Avr6oH9m9s1+E9TVKSRtx8Af/oAp+TJI2A+c6ieUGSB2dpD3BER/VIkhbJnAFfVcuWshBJ0uLyptuS1FMGvCT1lAEvST3VecAnWZbkm0k+1/VYkqTHLcUM/i3AjiUYR5I0wyC7SS5YktXAmcDvA2/rcqxXfeTrXXYvPeF9+g0nD7sEHaSuZ/B/DLwd+OFcL0iyKclkksmpKTerlKTF0tkMPsnLgPuralt7X9dZVdVmYDPAxMRELXQ8ZxeStK8uZ/DrgZcnuRv4FHBako93OJ4kaYbOAr6q3llVq6tqDfBq4G+r6jVdjSdJ2pfnwUtST3V6Fs20qroOuG4pxpIkNZzBS1JPGfCS1FMGvCT1lAEvST1lwEtSTxnwktRTBrwk9ZQBL0k9ZcBLUk8Z8JLUUwa8JPWUAS9JPWXAS1JPGfCS1FMGvCT1lAEvST1lwEtSTxnwktRTBrwk9ZQBL0k9ZcBLUk8Z8JLUUwa8JPWUAS9JPXV4l50nuRt4CNgLPFZVE12Op0W2fQt8+QLYtRNWrobTz4d1G4ddlaQBdRrwrVOr6oElGEeLafsWuPo82LO7Od51T3MMhrw0JpYi4JfGZWcOu4J+2Xkj7H1k37Y9u5sZvQEvjYWu1+AL+FKSbUk2zfaCJJuSTCaZnJqa6rgcDWz/cJ+2a+fS1iFpwbqewa+vqnuTPB24NskdVbV15guqajOwGWBiYqIWPNK5nz+kQrWfC9c2yzL7W7l66WuRtCCdzuCr6t728X7gKuCkLsfTIjr9fFi+Yt+25SuadkljobOAT3JkkqOmfwZeCtza1XhaZOs2woaLYOVxQJrHDRe5/i6NkS6XaJ4BXJVkepxPVNUXOhxPi23dRgNdGmOdBXxV3QW8oKv+JUnz80pWSeopA16SesqAl6SeMuAlqacMeEnqKQNeknrKgJeknjLgJamnDHhJ6ikDXpJ6yoCXpJ4y4CWppwx4SeopA16SesqAl6SeMuAlqacMeEnqKQNeknrKgJeknjLgJamnDHhJ6ikDXpJ6yoCXpJ4y4CWppzoN+CSrklyZ5I4kO5Kc3OV4kqTHHd5x/x8AvlBVr0zyJOApHY8nSWp1FvBJjgZeDLwOoKoeBR7tarxXfeTrXXX9hPfpN/jBSxpHXS7RPBuYAi5L8s0kFyc5cv8XJdmUZDLJ5NTUVIflSNITS6qqm46TCeAbwPqquiHJB4AHq+rfz/WeiYmJmpyc7KQeSeqjJNuqamK257qcwe8EdlbVDe3xlcALOxxPkjRDZwFfVd8F7knyvLbpdOD2rsaTJO2r67No3gxc0Z5BcxdwbsfjSZJanQZ8Vd0MzLo2JEnqlleySlJPGfCS1FMGvCT1lAEvST3V2YVOC5FkCvjOAt9+DPDAIpbTpXGqFcar3nGqFcar3nGqFcar3kOp9R9X1dNme2KkAv5QJJmc62quUTNOtcJ41TtOtcJ41TtOtcJ41dtVrS7RSFJPGfCS1FN9CvjNwy7gIIxTrTBe9Y5TrTBe9Y5TrTBe9XZSa2/W4CVJ++rTDF6SNIMBL0k9NfYBn+SMJN9KcmeSdwy7nvkkuTTJ/UluHXYtB5LkuCRfaW+WfluStwy7pvkkOSLJ/0jyP9t6/8OwazqQJMvau519bti1HEiSu5PckuTmJCN9V54kq5JcmeSO9vd3ZO85meR57d/p9J8Hk7x10fof5zX4JMuAbwMvobnByI3A2VU1kvvOJ3kx8DDw51W1dtj1zCfJscCxVXVTkqOAbcArRvjvNsCRVfVwkuXAV4G3VNU3hlzanJK8jWa31aOr6mXDrmc+Se4GJqpq5C8cSvIx4PqqurjdqvwpVfX9IZd1QG2e/R/gRVW10As+9zHuM/iTgDur6q72pt6fAs4ack1zqqqtwP8ddh2DqKr7quqm9ueHgB3AM4db1dyq8XB7uLz9M7KzlySrgTOBi4ddS58kORp4MXAJQFU9Og7h3jod+F+LFe4w/gH/TOCeGcc7GeEQGldJ1gAnADcc4KVD1S553AzcD1w743aRo+iPgbcDPxxyHYMq4EtJtiXZNOxi5vFsYAq4rF3+ujjJkcMuakCvBj65mB2Oe8BnlraRnbWNoyQ/AXwGeGtVPTjseuZTVXur6nhgNXBSkpFcBkvyMuD+qto27FoOwvqqeiHwK8Ab2+XGUXQ4zb2fP1xVJwB/D4z0d3MA7VLSy4H/tpj9jnvA7wSOm3G8Grh3SLX0TruW/Rngiqr6i2HXM6j2I/l1wBnDrWRO64GXt+vanwJOS/Lx4ZY0v6q6t328H7iKZnl0FO0Eds749HYlTeCPul8Bbqqq7y1mp+Me8DcCz0nyrPZfwFcDnx1yTb3Qfml5CbCjqt4/7HoOJMnTkqxqf14B/BJwx1CLmkNVvbOqVlfVGprf2b+tqtcMuaw5JTmy/aKddrnjpcBInglWVd8F7knyvLbpdGAkTwzYz9ks8vIMdH/T7U5V1WNJ3gR8EVgGXFpVtw25rDkl+SRwCnBMkp3Au6vqkuFWNaf1wK8Dt7Tr2gDvqqprhlfSvI4FPtaeiXAYsKWqRv70wzHxDOCq5t98Dgc+UVVfGG5J83ozcEU76bsLOHfI9cwryVNozgR8w6L3Pc6nSUqS5jbuSzSSpDkY8JLUUwa8JPWUAS9JPWXAS1JPGfAaqiR72130bmt3gnxbkkX7vUzyuiT/aMbxxUl+ZpH6fkWS8xejr1n6PmXQXSbbawBG+bRFDclYnwevXtjdbi9AkqcDnwBWAu8etIMky6pq7xxPv47mopzpKzF/81CK3c/baS4vH6qqmkpyX5L1VfW1Ydej0eEMXiOjvQx+E/CmNF6X5IPTzyf5XJJT2p8fTnJBkhuAk5Ocn+TGJLcm2dy+/5U02/Fe0X5KWJHkuiQTbR9nt3uc35rkvTPGeTjJ77efKL6R5Bn715rkucAjVfVAu8nZXe2Yq5L8cHqvliTXJ/np9mrQS9sav5nkrPb5ZUn+qG3fnuTHLnZJcmL7nmcn+cU8vnf4N6evMAX+EjhnEf43qEcMeI2UqrqL5vfy6Qd46ZHArVX1oqr6KvDBqjqx3Wd/BfCyqroSmATOqarjq2r39JvbZZv3AqcBxwMnJnnFjL6/UVUvALYCr59l/PXA9HbKe2nuS/AzwM/T7J3/C0meDKyuqjuBf0ezJcGJwKnAH7WX/f8bYFfbfiLw+iTPmlHnzwF/CpzV/t38LvDG9lPPLwDT/02T7bH0Iwa8RtFsu4Tuby/NRmjTTk1yQ5JbaEL7nx7g/ScC11XVVFU9BlxBs484wKPA9Pr3NmDNLO8/lmZb2mnXt+9/MfCHNEF/Is1+SdDs3/KOdtuH64AjgJ9s21/btt8APBV4TvuefwJsBjZU1f9u274GvD/JecCqtnZotkj+0XcNEhjwGjFJnk0T3vcDj7Hv7+gRM37+wfS6e5IjgA8Br6yqfwZ8dL/XzjrUPM/tqcf38NjL7N9V7d5vjOtpZtAnAdcAq2j2Hdo6Y7xfbT9JHF9VP1lVO9r2N89of1ZVfal9z33AD2j24gegqv4T8Js0n1K+keT57VNH8PhsXgIMeI2QJE+jWY74YBuwdwPHJzksyXHMvUXtdNA+kGb/+lfOeO4h4Kgffws3AL+Y5Jh2g7Kzgf9+EOXuAH56v/5+DvhhVf0AuJlm86jr2+e/CLy53aWTJCfMaP+tNFszk+S5efwGFd+nuevTH8z47uGnquqWqnovzbLMdMA/lxHd4VHD41k0GrYV7fLEcpoZ++XA9PbEXwP+DriFJrxumq2Dqvp+ko+2r7ubx5dFAP4M+NMku4GTZ7znviTvBL5CM4u+pqr+6iDq3gq8L0na2wU+kuQeYPoesNfT/KNxS3v8H2nu4rS9Dfm7gZfR3LJvDXBT2z4FvGJGnd9LsgH46yS/Abwmyak0nyxuB/66fempwOcPon49AbibpLRAST4AXF1VfzMCtWyl+SL2/w27Fo0Ol2ikhfsD4CnDLqJd2nq/4a79OYOXpJ5yBi9JPWXAS1JPGfCS1FMGvCT1lAEvST31/wGK9tuQ8I2UkgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plot_lifelines(shifted)\n", "plt.xlabel('Duration (weeks)');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that the x-axis in this figure is duration, not time.\n", "\n", "The second key idea is that we can estimate the hazard function by considering:\n", "\n", "* The number of dogs adopted at each duration, divided by\n", "\n", "* The number of dogs \"at risk\" at each duration, where \"at risk\" means that they *could* be adopted.\n", "\n", "For example:\n", "\n", "* At duration 1, there is 1 adoption out of 7 dogs at risk, so the hazard rate is `1/7`.\n", "\n", "* At duration 2, there is 1 adoption out of 6 dogs at risk, so the hazard rate is `1/6`.\n", "\n", "* At duration 4, there is 1 adoption out of 4 dogs at risk, so the hazard rate is `1/4`.\n", "\n", "And so on. Now let's see how that works computationally." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Computing \"at risk\"\n", "\n", "For each observed duration, we would like to compute the number of dogs that were at risk.\n", "Here are the unique durations, in order:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "tags": [ "fill-in" ] }, "outputs": [ { "data": { "text/plain": [ "array([1, 2, 4, 5, 7])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ts = duration.unique()\n", "ts.sort()\n", "ts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To compute the number of dogs at risk, we can loop through `ts` and count the number of dogs where `t` is less than or equal to `end`." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "scrolled": true, "tags": [ "fill-in" ] }, "outputs": [ { "data": { "text/plain": [ "1 7\n", "2 6\n", "4 4\n", "5 3\n", "7 1\n", "dtype: int64" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "at_risk = pd.Series(0, index=ts)\n", "\n", "for t in ts:\n", " k = (t <= shifted['end'])\n", " at_risk[t] = k.sum()\n", " \n", "at_risk" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you don't like mixing for loops with array operations, we can do the same computation using mesh grids." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(5, 7)" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "E, T = np.meshgrid(shifted['end'], ts)\n", "T.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The results are arrays with one row for each value of `t` and one column for each dog.\n", "Now we can use comparison operators to compare all values of `t` to all values of `end` at the same time." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([7, 6, 4, 3, 1])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "at_risk = (T <= E).sum(axis=1)\n", "at_risk" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The result is an array with the number of dogs at risk for each value of `t`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Estimating the hazard function\n", "\n", "Now, to compute the hazard function, we need to know the number of dogs adopted at each value of `t`." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "tags": [ "fill-in" ] }, "outputs": [ { "data": { "text/plain": [ "1 1\n", "2 1\n", "4 1\n", "5 1\n", "7 0\n", "dtype: int64" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adopted = pd.Series(0, index=ts)\n", "\n", "for t in ts:\n", " k = (shifted['status'] == 1) & (t == shifted['end'])\n", " adopted[t] = k.sum()\n", " \n", "adopted" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or here's the same computation with array operations:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 5., 1., 4., nan, nan, 2., nan])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adopt_times = np.where(shifted['status'], shifted['end'], np.nan)\n", "adopt_times" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(5, 7)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A, T = np.meshgrid(adopt_times, ts)\n", "T.shape" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1, 1, 1, 1, 0])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "adopted = (T == A).sum(axis=1)\n", "adopted" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the next step, it will be easier to see what we're doing if we put the results in a table." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "tags": [ "fill-in" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
adoptedat_risk
117
216
414
513
701
\n", "
" ], "text/plain": [ " adopted at_risk\n", "1 1 7\n", "2 1 6\n", "4 1 4\n", "5 1 3\n", "7 0 1" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d = dict(adopted=adopted, \n", " at_risk=at_risk)\n", "df = pd.DataFrame(d, index=ts)\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, the hazard function is the ratio of `adopted` and `at_risk`:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "tags": [ "fill-in" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
adoptedat_riskhazard
1170.142857
2160.166667
4140.250000
5130.333333
7010.000000
\n", "
" ], "text/plain": [ " adopted at_risk hazard\n", "1 1 7 0.142857\n", "2 1 6 0.166667\n", "4 1 4 0.250000\n", "5 1 3 0.333333\n", "7 0 1 0.000000" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df['hazard'] = df['adopted'] / df['at_risk']\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Working backwards\n", "\n", "Given the hazard function, we can work backwards to compute the survival curve.\n", "\n", "The hazard function is the probability of being adopted at each duration, so its complement is the probability of *not* being adopted. \n", "\n", "In order to survive past `t`, a dog has to *not* be adopted at all durations up to and including `t`.\n", "\n", "So the survival function is the cumulative product of the complement of the hazard function." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "tags": [ "fill-in" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
adoptedat_riskhazardsurv
1170.1428570.857143
2160.1666670.714286
4140.2500000.535714
5130.3333330.357143
7010.0000000.357143
\n", "
" ], "text/plain": [ " adopted at_risk hazard surv\n", "1 1 7 0.142857 0.857143\n", "2 1 6 0.166667 0.714286\n", "4 1 4 0.250000 0.535714\n", "5 1 3 0.333333 0.357143\n", "7 0 1 0.000000 0.357143" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df['surv'] = (1 - df['hazard']).cumprod()\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The CDF is the complement of the survival function." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "tags": [ "fill-in" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
adoptedat_riskhazardsurvcdf
1170.1428570.8571430.142857
2160.1666670.7142860.285714
4140.2500000.5357140.464286
5130.3333330.3571430.642857
7010.0000000.3571430.642857
\n", "
" ], "text/plain": [ " adopted at_risk hazard surv cdf\n", "1 1 7 0.142857 0.857143 0.142857\n", "2 1 6 0.166667 0.714286 0.285714\n", "4 1 4 0.250000 0.535714 0.464286\n", "5 1 3 0.333333 0.357143 0.642857\n", "7 0 1 0.000000 0.357143 0.642857" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df['cdf'] = 1 - df['surv']\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And the PMF is the difference between adjacent elements of the CDF." ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "tags": [ "fill-in" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
adoptedat_riskhazardsurvcdfpmf
1170.1428570.8571430.1428570.142857
2160.1666670.7142860.2857140.142857
4140.2500000.5357140.4642860.178571
5130.3333330.3571430.6428570.178571
7010.0000000.3571430.6428570.000000
\n", "
" ], "text/plain": [ " adopted at_risk hazard surv cdf pmf\n", "1 1 7 0.142857 0.857143 0.142857 0.142857\n", "2 1 6 0.166667 0.714286 0.285714 0.142857\n", "4 1 4 0.250000 0.535714 0.464286 0.178571\n", "5 1 3 0.333333 0.357143 0.642857 0.178571\n", "7 0 1 0.000000 0.357143 0.642857 0.000000" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df['pmf'] = np.diff(df['cdf'], prepend=0)\n", "df" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## lifelines\n", "\n", "Kaplan-Meier estimation is available in a library called `lifelines`.\n", "First I'll import it and create a `KaplanMeierFitter`." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "# If we're running in Colab, install lifelines\n", "\n", "import sys\n", "IN_COLAB = 'google.colab' in sys.modules\n", "\n", "if IN_COLAB:\n", " !pip install lifelines" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "from lifelines import KaplanMeierFitter\n", "kmf = KaplanMeierFitter()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we need two sequences, the durations, including complete and ongoing cases." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 5\n", "1 1\n", "2 4\n", "3 7\n", "4 5\n", "5 2\n", "6 2\n", "Name: end, dtype: int64" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "T = shifted['end']\n", "T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And an event flag that indicates whether a case is complete." ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "E = shifted['status']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `fit` method does the Kaplan-Meier estimation." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/downey/anaconda3/envs/SurvivalAnalysisPython/lib/python3.7/site-packages/lifelines/fitters/kaplan_meier_fitter.py:262: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version. Convert to a numpy array before indexing instead.\n", " self.confidence_interval_ = self._bounds(cumulative_sq_[:, None], alpha, ci_labels)\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "kmf.fit(T, E)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now the `Fitter` object contains the estimated survival function." ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
KM_estimate
timeline
0.01.000000
1.00.857143
2.00.714286
4.00.535714
5.00.357143
7.00.357143
\n", "
" ], "text/plain": [ " KM_estimate\n", "timeline \n", "0.0 1.000000\n", "1.0 0.857143\n", "2.0 0.714286\n", "4.0 0.535714\n", "5.0 0.357143\n", "7.0 0.357143" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "kmf.survival_function_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`timelines` includes an element at `t=0`, but other than that it is identical to what we computed (except for floating-point error)." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1.6653345369377348e-16" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "max(abs(kmf.survival_function_['KM_estimate'] - df['surv']).dropna())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`lifelines` also computes a confidence interval for the survival function." ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
KM_estimate_lower_0.95KM_estimate_upper_0.95
0.01.0000001.000000
1.00.3340540.978561
2.00.2581540.919797
4.00.1319880.824997
5.00.0519770.698713
7.00.0519770.698713
\n", "
" ], "text/plain": [ " KM_estimate_lower_0.95 KM_estimate_upper_0.95\n", "0.0 1.000000 1.000000\n", "1.0 0.334054 0.978561\n", "2.0 0.258154 0.919797\n", "4.0 0.131988 0.824997\n", "5.0 0.051977 0.698713\n", "7.0 0.051977 0.698713" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ci = kmf.confidence_interval_survival_function_\n", "ci" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With such a small dataset, the CI is pretty wide." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1pklEQVR4nO3deXzcZbX48c+ZfZ9kJkm3pLRCoRSQAkGWYlX2YmlxB7n398MF5ApckCuKigLCVRFZL6L0J+DGBRGkFGQTZG3ZWih0SQst3dI2+74v8/z+mJmYplkmy2S283698iKTmXznkCZz5jnP85xHjDEopZTKXZZUB6CUUiq1NBEopVSO00SglFI5ThOBUkrlOE0ESimV42ypDmC0CgoKzKxZs1IdhlJKZZQ1a9bUGGMKB7sv4xLBrFmzWL16darDUEqpjCIiO4a6T0tDSimV4zQRKKVUjtNEoJRSOS7j5giUUumnu7ub8vJyOjo6Uh1KznO5XBQXF2O32xP+Hk0ESqlxKy8vx+/3M2vWLEQk1eHkLGMMtbW1lJeXM3v27IS/L2mlIRG5T0SqRGT9EPeLiNwpIltE5H0ROTpZsSilkqujo4NwOKxJIMVEhHA4POqRWTLnCH4PnDnM/YuAObGPi4DfJDEWpVSSaRJID2P5d0haIjDGvALUDfOQpcAfTdQbQJ6ITEtWPLUtnVz/xAaaO7qT9RRKKZWRUjlHMAPY1e92eexrewc+UEQuIjpqYObMmWN6spVba/nDqu089d5urjltJiccWIDdbsdms2Gz2bBarVitVn1Xo5TKOalcPjrYK+6gp+QYY5YZY0qNMaWFhYPukB7RkiOn8/t/+zhiernssa1cu3wd6zd9wMaNG1m3bh1r167lnXfeYd26dWzevJnt27ezd+9e6urqaGpqoq2tja6uLnp7e8f0/Eqp5PL5fH2fP/XUU8yZM4edO3dy3XXXISJs2bKl7/7bbrsNEZmwLgXLly9n48aNfbd/8pOf8Pzzz4/7ug0NDdx9993jvs5IUjkiKAdK+t0uBvYk8wkXzJ3Or06r4C9lHTyxuYn3q7v53sIpHFzgAqIz7r29vXR2dtLe3k5vby/GGAae4ma1WnE4HDidzr4Ph8PRN7qIf+joQqnJ98ILL3DZZZfx3HPP9VUQjjjiCB566CGuueYaAB555BHmzZs3Yc+5fPlyFi9e3HfNn/70pxNy3Xgi+Pa3vz0h1xtKKhPBCuBSEXkIOA5oNMbsVxaaSFarlakF+Xzz6FZOOMDHra9VcvmT5Zx3ZD5fPTKEzSJ9L+LDiUQi9Pb20t7eTktLy36jBGMMIoLVat0nWTidzn3KUfGSlCYMlU2uf2IDG/c0Teg1500PcO3Zh434uFdffZULL7yQp556igMPPLDv6+eccw6PP/4411xzDR999BHBYHDEdfbPPfcc1157LZ2dnRx44IHcf//9+Hw+rr76alasWIHNZuP000/n85//PCtWrODll1/mxhtv5NFHH+WGG25g8eLFfPGLX2TWrFl89atf5cUXX6S7u5tly5bxgx/8gC1btnDVVVdx8cUX09LSwtKlS6mvr6e7u5sbb7yRpUuXcvXVV7N161bmz5/Paaedxs0338zNN9/Mww8/TGdnJ5/73Oe4/vrrx/3zTVoiEJEHgU8DBSJSDlwL2AGMMb8FngLOArYAbcDXkhVLf+FwmIaGBo6Z4eeec2Zy95s1PLC2njd3tfG9hVM4IM8x4jUsFgsWy/BVNWMMkUiESCRCS0sLjY2NRCKRQR/rcDj6RhgOhwOXy4XNZsNut+NyuUZ8LqUUdHZ2snTpUl566SXmzp27z32BQICSkhLWr1/P448/zle+8hXuv//+Ia9VU1PDjTfeyPPPP4/X6+Wmm27i1ltv5dJLL+Wxxx5j06ZNiAgNDQ3k5eWxZMmSvhf+wZSUlPD666/zne98hwsuuICVK1fS0dHBYYcdxsUXX4zL5eKxxx4jEAhQU1PD8ccfz5IlS/jFL37B+vXrWbt2LRBNTh9++CFvvfUWxhiWLFnCK6+8wsKFC8f1s0taIjDGnDfC/Qa4JFnPP5R4HdEYg89p5XsLp3DiTC93rqrikhW7uODoEJ+bl4fVMr536fERgdVqHfadRzxhdHd309HR0VeOit9ntVrJy8sjPz8fr9eLwzFyolIqlRJ5554MdrudE088kXvvvZc77rhjv/vPPfdcHnroIZ599lleeOGFYRPBG2+8wcaNG1mwYAEAXV1dnHDCCQQCAVwuF9/85jf57Gc/y+LFixOKbcmSJUC0RNXS0oLf78fv9+NyuWhoaMDr9fLDH/6QV155BYvFwu7du6msrNzvOs899xzPPfccRx11FAAtLS18+OGH6ZsI0pXdbsfn89HV1YXT6QTgpFk+Dpvi4o5V1fy/t2t5Y2cr3/3kFKb6E9+iPVb9E8ZgIpEITU1N1NXVYYzB7XYTDofx+/14PB4dLSgVY7FYePjhhzn11FP52c9+xg9/+MN97j/77LO56qqrKC0tJRAIDHstYwynnXYaDz744H73vfXWW7zwwgs89NBD3HXXXfzzn/8cMbb4a43FYun7PH67p6eHBx54gOrqatasWYPdbmfWrFmDbgozxvCDH/yAb33rWyM+52jk5KtIQUEBnZ2d+3wt323j2pOn8l8nFbG1rouLl+/k6Q8a95sonmwWiwW3243f7+/75d2zZw9lZWW8++67bNmyhdra2v3+f5TKRR6PhyeffJIHHniAe++9d5/73G43N910Ez/60Y9GvM7xxx/PypUr+1YatbW18cEHH/SVec866yxuv/32vpKN3++nubl5zHE3NjZSVFSE3W7nxRdfZMeOHYNe94wzzuC+++6jpaUFgN27d1NVVTXm543LuREBRH+4gxERTp8TYP40N7e8VsXtK6tZtaOVKxYUEfakx4/Kbrf3lZoikQitra3U19cD0Xcd4XCYQCCAx+MZcpShVDYLhUI888wzLFy4kIKCgn3uO/fccxO6RmFhIb///e8577zz+t5k3Xjjjfj9fpYuXUpHRwfGGG677ba+61544YXceeedPPLII6OO+fzzz+fss8+mtLSU+fPn981xhMNhFixYwOGHH86iRYu4+eabKSsr44QTTgCipe4///nPFBUVjfo5+5NUv+MdrdLSUjMRa383bNiAMWbI+n3EGFaUNXLv6lqcNuGyEwr51OzBE0i66OnpoaOjg0gkgogQCAQIhUL4fD6cTqeuTlJJU1ZWxqGHHprqMFTMYP8eIrLGGFM62OPT421uCoTDYcrLy4dMBBYRzpmXx9HTPfzq1Up+9lIlq3a0cskJhQSc6flO22az7TMZ3t7ezrZt24DoyqRQKEQwGMTj8Yy4RFYplTty9tUgGAyya9euER83M8/BbZ8t5i/v1/PntXW8X9HOlScVcWyxdxKiHDsRweVy4XJFN8v19PRQXV1NRUUFIoLP5yMcDuPz+XC5XDpaUDnruOOO22+O7U9/+hNHHHFEiiKafDmbCFwuFw6Hg56enhHfHVstwlfnhzi22MPNr1ZxzT/2ctbBAS76RAFue2bMt/ffKGeMoaurix07dmCMwWazEQqFyMvLw+v16mhBjUl8I2WmefPNN1MdwoQaS7k/Z//iRYSCggIqKir26VEynDkFLu46u5g/vlvHI+sbeGdPG9/95BSOmOpOcrQTS0T6djoD9Pb2Ultb27f6wOfzEQqF8Pv9uN3ujPzjVpPL5XJRW1urZxKkWPxgmnglIFE5O1kM0c0YZWVlI64pHsy6inZuea2SiuYevnB4Hv/3qBAOW2aMDoZjjKG7u5vOzs6+DW35+fnk5+fj8Xh0Q5salB5VmT6GOqpyuMninE4EkUiEtWvX4nK5xrTUsr07wrK3a3hqcxMH5Dn43sIpHBR2jvyNGSTehK+npweIrtMOhUIEAgHcbrduaFMqQ2giGMaOHTuora3F6x375O/b5a3c+loVjR29nD8/xLkfzx93i4p0ZIyhp6eHzs5OIpEIVquVYDBIKBTC4/Hss2NSKZVedPnoMPLz88e9M+/YYi/3fG4mv36jmj++W8ebu1q5auEUSoLZVUYRkf02tDU3N1NXFz2Izul0EggE+kYLuhpJqcyQ8yOC3t5e1q5dO2F9e17e1sz/vF5NZ4/hG6VhlhwaxJIjL4Y9PT19h/cYY7BYLPj9foLBIF6vt6+rqlJq8umIYBjx7p7Nzc243eNf/fOp2X4On+Lm9pVV/ObNGlbtaOW/PlnEFF/yG9il2sCzHCKRCB0dHTQ3N/ctaXO5XASDwb4VSQ6HQ0cNSqVYzo8IAOrr69myZcuYVg8NxRjDMx82cc+bNYjAxccVcvpB/px+0YvPMfQ/8tNms/WVkzwej05AK5UkOiIYgc/nQ0QmdEOMiLDo4CBHTYu2qLj1tapYA7tC8t25+WMfOMcA0dJcS0sL9fX1fT9/r9dLMBjE5/PhdrtHPElKKTU+OiKI2bRp0z5nFEykiDE8tqGB+9+pw2MTLjuxiE/OSmwTW66J72Po6urqa57ncDh0ElqpcdLlowmoqalh+/btQ7aongg7Grr45SuVbKnt5JQD/Xz7uAJ8adrALp3oJLRS46eloQT4/f6kH0JzQJ6DOxYX8+B7dfzve/W8t7eNK0+awjEzPEl93kynk9BKJZeOCPrZsGEDkUhkUtoofFDTwS9fqWRXYzdnzw3yzdIwrgxpYJduhpuEDgaDuN1unYRWOU9HBAkqKChg165dk5IIDi5w8eslJfz+nVoe29DImt1tfPeTRRw2JbMa2KWDoSahm5ub95uEjndY1Ulopf5FE0E/gUBgUs8odtosfOsThRxf4uWW16r47tO7+dLhefzbUWEcVi1rjIfVat1nX0h8EnrPnj06Ca3UAJoI+nG5XDidzoTOKJhIR07z8JulM1n2Vg1/WdfAW+VtfG/hFD4W0t49EyX+wt9/tNfT00N9fT01NTU6Ca1ymv6W9yMihMNhKisrJ/0FwOuw8J2TijhhppfbV1Vx2RO7+PejQnzp8OxsYJcOdBJaqShNBAPk5eWxZ8+elD3/8TO93FM0k/95vYr719Txxs42rvpkETOyrIFdOrJYLPsc2BOfhI4f8Qk6Ca2ykyaCAeIHu/f29o7pjIKJEHRZ+dGnp/LiRy38+o1q/mPFLi4sLWDx3IC+G51EI01Cx+kktMp0mggGEBFCodC4zyiYiDhOPtDPx6e6uXVlFXe9Uc2qnS1cedIUCr36z5YqOgmtspHuIxhEU1MTmzdvntAmdONhjOHvm5tY9nYNNovw7eMKOOXA3G5gl850J7RKR7qPYJS8Xi8Wi4VIJJIW9V8RYfHcIEdPjzawu/nVKlbtbOU/Tywiz6UtKtKNTkKrTKMjgiFs3bp1ws4omEi9EcOjGxr44zu1eB1WrlgQXWmkMofuhFapoCOCMQiHw/tMCKYLq0X48hH5HDvDwy9freS6F/Zy2kF+/uO4ArwOHR1kgrFMQvt8Pk0MKmk0EQzB54u2iZ7IMwom0uyQkzsXl/DA2jr+sq6e9/a281+fLGL+NG1gl4lGmoR2uVwccMABaTNvpbKLvsUYQnyo3tXVlepQhmS3ChccE+bWs4qxW4XvP7OH37xRTUdPJNWhqXGKrz7y+Xx9rU82bdrERx99lNa/kyozJTURiMiZIrJZRLaIyNWD3B8UkSdE5D0R2SAiX0tmPKMVCoUy4o/u0CIXdy8tYemhQZaXNXLJil1squ5IdVhqAjmdTgKBAA0NDaxbt47KykoiEU34amIkLRGIiBX4NbAImAecJyLzBjzsEmCjMeZI4NPALSKSNltoJ7sJ3Xi4bBa+fXwhvzhjOh3dhu/8vZw/vFNLd29mxK9GFu+g6na72bFjBxs3bqS5uTnVYakskMwRwSeALcaYj4wxXcBDwNIBjzGAX6JFeB9QB/QkMaZRcTgceDyejBgVxB013cM955Rw8sf8/O979VzxZDnb6ztTHZaaQFarlWAwSCQSYdOmTWzbti2jfkdV+klmIpgB7Op3uzz2tf7uAg4F9gDrgMuNMfuNd0XkIhFZLSKrq6urkxXvoAoKCujszKwXUp/TylULp3DtyVOpbuvh0hW7+Ou6enojOjrIJk6nE7/fT11dHevXr6eqqkrLRWpMkpkIBltqM/CV6AxgLTAdmA/cJSL7LYswxiwzxpQaY0oLCwsnOs5hZVJ5aKATD/Cx7JwSji328rvVtVz1zG72NnenOiw1gUQEn8+Hy+Vix44dlJWV0dLSkuqwVIZJZiIoB0r63S4m+s6/v68BfzNRW4BtwNwkxjRq/c8oyER5bhs/OXkq3/1kEdvqurh4+U6e2tyYsclNDc5qtRIIBOjp6aGsrIwdO3bQ3a1JXyUmmYngbWCOiMyOTQCfC6wY8JidwCkAIjIFOAT4KIkxjZqIUFBQQEdH5q7CERFOOyjAPeeUMLfQxR2rqvnxP/ZS25aZyU0NzeVy4ff7qampYd26dX2H7ig1nKQlAmNMD3Ap8CxQBjxsjNkgIheLyMWxh90AnCgi64AXgO8bY2qSFdNYBYPBrPhjKvLZ+fkZ0/n2cQW8X9HORY/t5KWPdNVJtomXi5xOJ9u2bWPTpk20tramOiyVxrTXUAKMMaxduxan05myMwomWnljFze/Wsmm6k4WzvJx2QmFBLSBXVbq6Oigq6uLKVOmMH36dO18mqOG6zWkO4sTkA3loYGKgw5uPauYC44OsWpnC99avpO3dum7xmzkcrnw+XxUVVWxbt06amtrs2KEqyaOJoIE5eXl9XWKzBZWi3DekSHuXFxCwGnlx8/v5faVVbR16xLEbBM/E8HhcLB161Y2b95Me3t7qsNSaUITQYL6n1GQbQ4MO/mfJSV86fA8nvmgif9YvpN1FfoikY1sNhvBYJCOjg7Wr19PeXl5xq6IUxNHE0GCLBYL+fn5Gbe5LFEOq/DNYwu45awZiAhXPb2bZW/V0KUN7LKS2+3G5/NRUVHB+vXrqa+v13JRDtNEMAqhUCjr12YfNsXNb5aWcNYhAR7d0MAlT5TzYU32zI2of4mXi2w2Gx9++CEffvihlotylCaCUfD5fIhI1r9zctst/OeJRdx42jRau3q5/Mly/ry2jh5tUZGV7HY7wWCQ1tZWNmzYwJ49e7JuPkwNTxPBKMTPKMjW8tBAxxZ7ueecmSyc7eNP79bxnb+Xs7NBm5tlK4/Hg9frZc+ePaxfv56Ghoasf9OjojQRjFI4HM6pTo9+p5WrPzWVH316KhXN3VyyYhd/29BARF8gslK8XGSxWPjggw/YunVrVi2bVoPTRDBKfr8/J8pDAy2c7eOec2Zy1HQ397xVw/ef2UOFNrDLWg6Hg0AgQHNzM+vXr2fv3r1aLspimghGKX5GQbZPGg8m5LFx/SnTuHJBER/WdPAfj+/k2Q+bci4p5goRwePx4PF4KC8vZ8OGDTQ1NaU6LJUEmgjGIBPPKJgoIsIZBwf47TkzOSjs5NbXqrjuhb3UaQO7rBXvbCoibNq0ia1bt+bs73+2GrHpiIgUAhcCs/o/3hjz9eSFld4y+YyCiTLVb+emM2ewfGMj962p5cLHdrJ4bpDFc4MUerWXTTZyOBzY7XYaGxupr6+npKSEwsJCLBZ9P5npEvmLfRx4FXge0CIh0ZOhnE4n3d3d2O32VIeTMhYRPn9YHqUzPNy3ppaH19Xz8Lp6FhzgY+mhQQ6f4iJ6CqnKFvFzk3t7e9m5cydVVVXMmjULv9+f6tDUOIzYfVRE1hpj5k9OOCNLRffRwezdu5c9e/bg8/lSHUraqGju5slNjTz9QRMtXREODDlYOi+PT8/24bTpu8Zs1NnZSUdHB4WFhcyYMQOHw5HqkNQQhus+mkgiuBFYZYx5KhnBjVa6JILW1lY2btxIILDfyZo5r6M7wj8/aubxjY1sb+gi4LSw6JAgiw8JUOTL3RFUtjLG0NbWhjGGkpISCgoKtFyUhsabCJoBL9AFxJfKGGNMSl4B0yURZOMZBRPNGMN7Fe2s2NjI67EW1wsO8LL00DwtG2Wh3t5eWltbsVqt+m+bJOFwmJkzZ47pe4dLBCPOERhjtPg3iPgZBdXV1Xi93lSHk5ZEhPnTPMyf5tmnbPTq9tZo2ejQPD79MS0bZYv46qJs7NCbDrq6upLWCyqhE8pEZAmwMHbzJWPMk0mJJgHpMiIAaG5uZvPmzTpRNgodPRH+uXVA2ejgAIvnBrVspNQwurq6cDgcHHLIIWP6/nGNCETkF8CxwAOxL10uIicZY64eUzRZpP8ZBVoTTYzLZuGsQ4IsOjjA+xXtPL6xkb+ub+Cv6xs4caaXc+Zp2UipyZbI8tGzgPnGmAiAiPwBeBfI+UQQP6OgoaEBj8eT6nAyiohw5DQPR07zUNnSzRNl0bLRazta+VisbPQZLRspNSkS/SvL6/d5MAlxZKxQKKQnPI3TFJ+dbx5bwANfmcUVCwoxBm5bWcX5D2/n3tU1VLXkXjsPpSZTIiOCnwPvisiLgBCdK/hBUqPKIF6vt68JnZYzxsdls7Do4CBnzgmwrrKD5RsbeGR99OPEmV6WzsvjCC0bKTXhElk19KCIvER0nkCA7xtjKpIdWKaIn1HQ3t6Oy+VKdThZQUT4+FQ3H5/qprLlX6uNXtvRyux8B0vnBfnMx/y4tGyk1IQY8i9JRObG/ns0MA0oB3YB02NfUzG5dkbBZJris/ON0gL+/OVZfGdBEQC3r6zm37RspNSEGW5EcCVwEXDLIPcZ4OSkRJSB4stHtTyUPC6bhTMPDnDGHD/rKjt4vF/Z6ITYaiMtGyk1NkMmAmPMRbFPFxlj9jmiSES0BtKPw+HA5/PR3d2tvVaSrH/ZqCpWNnrqgyZWatlIqTFL5K9lVYJfy2m5fEZBqhT57Hy9tIAHYmUjkX+VjX73dg2VWjZSKiFDjghEZCowA3CLyFFEJ4oBAoAumh/A7/fn/BkFqeLsVzZaX9nB42UNPLoh+nFCiZel84J8fKpby0ZKDWG4OYIzgAuAYqLzBPG/oibgh8kNK/O4XC5cLlfOn1GQSiLCEVPdHNGvbPT0B02s3BktGy05NMjJB2rZSKmBEuk++gVjzKOTFM+I0qnX0EB6RkH66eyJ8NJHLSwva+Cjui58jn/1Nprq14StMkcyew0l8tboGBHJ63ex/NgZBWqAYDConRfTjNNm4YyDA9y9pIRfLZrBUdPdPLqhga89uoPrX9jL2r1tWtJTOS+RncWLjDF9pSBjTL2InAVck7ywMpPb7cZms9Hb26tnFKSZ/cpGm5t4enMjq3a2MisvutpIy0YqVyXyW28VEWf8hoi4Aecwj+8jImeKyGYR2SIigzapE5FPi8haEdkgIi8nFnZ6EhEKCwuT1jNcTYwin52vHxPmz1+exZUnFWG1wB2rqjn/L9HVRhXNutpI5ZZERgR/Bl4QkfuJbiT7OvCHkb5JRKzAr4HTiO5KfltEVhhjNvZ7TB5wN3CmMWaniBSN/n8hvQSDQfbu3ZvqMFQCnDYLZ8wJcPpBfjZUdbB8Y2PfaqPjY6uNjtTVRioHJNJr6Jcisg44hejKoRuMMc8mcO1PAFuMMR8BiMhDwFJgY7/HfBX4mzFmZ+y5qkYZf9rxer1YrVY9oyCDiAiHT3Fz+JRo2ejvm5t4amDZ6GN+XHb991TZKZERAcaYp4GnR3ntGUR7E8WVA8cNeMzBgD3W1M4P3GGM+eMonyet6BkFma3IZ+drx4Q5/8h8XtrWwuMbG7ljVTX3rq7lzIMDnK2rjVQWSuSEss8DNwFFREcEQmKH1w82nh64PMMGHEN0tOEGXheRN4wxHwyI4SKifY/GfHDzZAqFQtTW1qY6DDUODpuF0+cEOC1WNnp8YyN/29DA3zY0cFyJl3MODXLkNC0bqeyQyIjgl8DZxpiyUV67HCjpd7sY2DPIY2qMMa1Aq4i8AhwJ7JMIjDHLgGUQ3Ucwyjgmnc/n0zMKskT/slF1aw9/39TIU5sbeX1nKwfkOVh6aJBTDtSykcpsifz2Vo4hCQC8DcwRkdki4gDOBVYMeMzjwCdFxCYiHqKlo7E8V1qxWq0Eg0HtPZRlCr02LoitNvqvk4qwW4Q7X6/m/Ie3s+wtXW2kMlciI4LVIvIXYDnQ98pmjPnbcN9kjOkRkUuBZwErcJ8xZoOIXBy7/7fGmDIReQZ4H4gAvzPGrB/b/0p6CYVCNDQ06GE1Wah/2WhjVQePlzXy2MZo2ej4mV6WHhpkvpaNVAZJJBEEgDbg9H5fM8CwiQDAGPMU8NSAr/12wO2bgZsTiCOjxJvQaXkoe4kIh01xc1j/stEHTVo2UhlnxF5D6Sadew0NVFZWRnd3N05nQvvvVBbo6onw8rYWlpc1sqW2E58juldhyaG62kiNTzJ7DSWyaii+kWwfxpivjymaHFJQUMCOHTs0EeQQh83CaXMCnDpI2ei4Ei/nzNOykUo/iZSGnuz3uQv4HPuv/lGDCAQC2tAsR/UvG9W09vD3zY38fXMTbzzbysxY2ehULRupNDHq0pCIWIDnjTEpObM4k0pDAOvWrUNE9IwCFS0bbY9uUvuwthNvv7LRNC0bqRGktDQ0iDlA+u/qShOFhYXs3r1bE4GKlo0OCnDqgX7KqqOb1B7f2MBjGxo4rsTD0nl5HKVlI5UCicwRNLPvHEEF8P2kRZRlAoEAu3btGvmBKmeICPOK3MwrcnNRWw9PboqXjfYwM2hnybw8Tj3Qj1vLRmqSDHdm8QJjzEqg0BjTMYkxZRW3243D4aCnpwebbSwDMJXNwh4b//foMOcdGeKVbc0s39jIXa9Xc/+aWs6Y4+fsuXlMD+hoUiXXcK9MdxLtA7QKOHpywsk+IkI4HKayslKPsFRDcliFUw8KcEq8bFTWyOMbG3lsQyMnzPRy+YIi8lx62JFKjuESQXds6WixiNw58E5jzH8mL6zskpeXp2cUqITsUzY6Nlo2emR9A1c8uYufnjqdmXmOVIeostBwRcjFRNtDtANrBvlQCfJ4PH1nFCiVqHjZ6OZFM2jvNnzn7+Ws3dOW6rBUFhpyRGCMqQEeEpEyY8x7kxhT1rFYLIRCIerr6/WMAjVqcwtd3LG4mB8/v5cfPreHyxcUccackbrAK5W4EZclaBKYGPn5+fT09KQ6DJWhpvrt3HbWDD4+zc2tr1Vx35paIrpZUU0QXZ82SeJnFGh5SI2Vz2nlxtOms+jgAH95v56fv1RJZ4/+Pqnx00QwSaxWK3l5eXR1daU6FJXBbBbh8hML+WZpmFe2t/D9Z/bQ0NGb6rBUhhtuH8GVw32jMebWiQ8nu4XDYerr6/WMAjUuIsKXjshnqt/OL1+p5PIndnHDabqiSI3dcCMC/wgfapT6H2Gp1Hh9cpaPmxfNoLPXcIWuKFLjMNyqoesnM5BcYLfb8fl8dHV1aWtqNSH6VhT9I7ai6MQizjhYVxSp0Umk15AL+AZwGNE21ICeRzBWBQUFbNu2TROBmjBTfHZu++wM/vvFCm5dWcXu5m4uODqERZvXqQQlMln8J2AqcAbwMlAMNCczqGzm92tVTU08r8PKT0+bzlmHRFcU/UxXFKlRSCQRHGSM+THQaoz5A/BZ4IjkhpW9nE4nbreb7u7uVIeisozNIvznCYVceGyY17a38P1ndtPQrntX1MgSSQTxV6wGETkcCAKzkhZRDigoKKCjQxu6qoknInzx8Hx+fPJUPqrr4j+fLGdHgy5ZVsNLJBEsE5F84MfACmAjcFNSo8pywWBQVw6ppFpwQHRFUXdvtEfRO7qiSA0jkURwvzGm3hjzsjHmY8aYImPMPUmPLIu5XK6+MwqUSpZDCl3cvriYQo+Na57bw9MfNKY6JJWmEkkE20RkmYicInqG3oQQES0PqUkxxWfn1s/OYP40N7evrObe1TXao0jtJ5FEcAjwPHAJsF1E7hKRk5IbVvYLBoPad0hNCq/Dyg2xFUUPr2vgZy9V6IoitY9Euo+2G2MeNsZ8HpgPBIguI1XjED+joLdX+8So5LPus6Kolaue3k29rihSMQk1nRORT4nI3cA7RDeVfTmpUeUAi8VCOBzW8pCaNPEVRT85eSrb67u4/Mlyttd3pjoslQZGTAQisg24AngVONwY82VjzKPJDiwX5Ofn64hATboTD/Dxq7PiK4p2s2a3rijKdYmMCI40xnzOGPOgMaY16RHlEK/Xi8Vi0bkCNekOLoj2KJris3HNP/bw1GZdUZTLhmtD/T1jzC+B/xaR/ZYZ6OH14xc/o6C5uRm3253qcFSOKfLZueWsYn72UgV3rKpmT1M3Xy8Na4+iHDRc07my2H9XT0YguSoUClFXV6eJQKWE12Hhp6dO4+43q/nr+gb2Nndz1cIpuGx6ZlUuGa4N9ROxT983xrw7SfHknP5nFOg2DZUKVotw6fGFzAg4WPZWDdVP7+a6U6YR8ozYnFhliUTS/q0isklEbhCRw5IeUY7pf0aBUqkiInz+sLzoiqKGLq74u64oyiWJ7CP4DPBpoJpo36F1InJNIhcXkTNFZLOIbBGRq4d53LEi0isiX0w08GxSUFCgiUClhRMP8HHLIl1RlGsSKgQaYyqMMXcCFwNrgZ+M9D0iYgV+DSwC5gHnici8IR53E/Bs4mFnF7/fr03oVNqYoyuKck4i+wgOFZHrRGQ9cBewiujhNCP5BLDFGPORMaYLeAhYOsjjLgMeBaoSDzu7OJ1OPB6PjgpU2ijy2bn1s8UcM93DHauq+X9va4+ibJZQ91GgHjjdGPMpY8xvjDGJvGjPAHb1u10e+1ofEZkBfA747XAXEpGLRGS1iKyurq5O4KkzT0FBAZ2dWpNV6cNjt3D9qdM4e26QR9Y3cOM/K+jQHkVZadhEECvbbDXG3GGM2TPKaw+2BGbgW4rbge8bY4bdXmuMWWaMKTXGlBYWFo4yjMwQCAS0PKTSjtUiXHJ8Ad/6RAGrdkZ7FNW1aY+ibDNsIoi9QIdFxDGGa5cDJf1uFwMDk0kp8JCIbAe+CNwtIueM4bkynsvlwul06hkFKu3EVxRde8o0djRoj6JslEhpaAewUkR+LCJXxj8S+L63gTkiMjuWSM4lesJZH2PMbGPMLGPMLOAR4NvGmOWj+1/IDiJCOBzW8pBKWyfM9HLLohn0RKKnnq3erR1nskUiiWAP8GTssf5+H8MyxvQAlxJdDVQGPGyM2SAiF4vIxWMPOXvl5eVpEzqV1v61osjOj/+xlyc36YqibCCZVpcuLS01q1dnZ9cLYwxr167F6XRitVpTHY5SQ2rrjvDzlyp4q7yNLxyWxzdKw1gtujM+mbq6unA4HBxyyCFj+n4RWWOMKR3svhH3kIvIi+w/yYsx5uQxRaOGJCKEQiFqa2vxer2pDkepIXnsFq47ZRq/fauGRzdEexR9f+EUXHbtUZSJEmkm8t1+n7uALwA6o5kk+fn5VFXl7JYKlUGiK4oKmRGw89s3a/ju07u5/tRphLVHUcYZ8V/MGLNmwJdWiogeVZkk/c8osFj03ZVKf+fMy2Oqz87PX67g8ifLWTw3gAy6elyNR29vL4dN9TDGytCwEikNhfrdtADHAFMnPhQFekaBykzHz/Ryy1nFXP/CXu5fU5fqcLLWl49IzjnBiYzh1hCdIxCiJaFtwDeSEIuKCYfD1NfXpzoMpUbloLCT33/xAHoimbUAJVN0dXXhdjmTcu1ESkOzk/LMakg+nw9AzyhQGcdqEV09lCQSsWBL0s92yCJ0rDX01H63/4+IPC4idw4oF6kJZrPZCAQC2oROKTUphpuNvAfoAhCRhcAvgD8CjcCy5IeW20KhkCYCpdSkGK40ZDXGxGd9vgIsM8Y8CjwqImuTHlmOCwQCWCwWmpubcblc2O32VIeklMpSw40IrCISTxSnAP/sd58uFE4yh8PBYYcdRklJCcYYmpqaaG1t1RYUSqkJN9wL+oPAyyJSA7QDrwKIyEFEy0MqyZxOJ1OmTKGoqIj29nbq6uqorq6mp6cHm82Gy+XSvQZKqXEbMhEYY/5bRF4ApgHPmX81JbIQPVVMTRIRwePx4PF4mD59Oi0tLdTU1FBfX08kEsHpdOJwOHSFkVJqTIYt8Rhj3hjkax8kLxw1EovFQiAQIBAI0NPTQ1NTE1VVVTQ3NyMiOp+glBo1rfVnMJvNRigUIhQK0dnZSUNDA1VVVTQ1NWG1WnG5XNrFVCk1Ik0EWaL/fEJbWxv19fU6n6CUSogmgiwjIni9Xrxe7z7zCXV1dRhjdD5BKbUfTQRZrP98wsyZM2lsbKS6ulrnE5RS+9BEkCNsNhvhcLjvXOSGhgYqKyt1PkEppYkgF+l8glKqP00EOWyw+YTq6uq+FtgOh0PnE5TKAZoIFLD//gSdT1Aqd2giUPvR+QSlcosmAjWs4eYT7HY7TqdT5xOUynCaCFRCdD5BqeyliUCN2kjzCTabTUcKSmUQTQRqXPrPJ3R1ddHa2kptbS2NjY1EIhGdU1AqA2giUBMmXh7Kz8+nt7eX1tZWGhsbqamp6TtQR1cfKZV+NBGopLBarX3lo+LiYtrb22lqaqK2tpampiYgmjicTqfOKyiVYpoIVNL1P1hn6tSpdHZ20tLS0pcUjDFaQlIqhTQRqEnndDpxOp2Ew2F6e3tpaWmhoaGB2tpaIpEIIoLT6dQSklKTRBOBSimr1UowGCQYDDJz5kza2tpoamqipqamr4SkrbOVSi5NBCpt9N+rMG3aNDo7O2lubqa2tpbm5maMMdoUT6kk0ESg0la8hFRQUEBPTw8tLS3U19dTX1/fV0JyuVzYbPprrNR4JPVtlYicKSKbRWSLiFw9yP3ni8j7sY9VInJkMuNRmctms5GXl8fs2bOZP38+c+fOZerUqfT29tLU1ERzczOdnZ0YY1IdqlIZJ2lvpUTECvwaOA0oB94WkRXGmI39HrYN+JQxpl5EFgHLgOOSFZPKDhaLBZ/Ph8/nY/r06X0lpJqaGlpaWjDGaB8kpUYhmWPqTwBbjDEfAYjIQ8BSoC8RGGNW9Xv8G0BxEuNRWSheHnK5XBQWFtLd3U1LSwt1dXU0NDQQiUSwWCxaQlJqGMn8y5gB7Op3u5zh3+1/A3h6sDtE5CLgIoCZM2dOVHwqC9ntdvLz88nPzycSidDa2kpDQwN1dXW0tbXtszRVVyEpFZXMRDDYX9mgBVwR+QzRRHDSYPcbY5YRLRtRWlqqRWCVEIvFgt/vx+/3U1xcTEdHR9/u5ngJKd4WQ0tIKpclMxGUAyX9bhcDewY+SEQ+DvwOWGSMqU1iPCqHiQhutxu3282UKVPo6urap4QE9JWQdHezyjXJTARvA3NEZDawGzgX+Gr/B4jITOBvwL8bYz5IYixK7cPhcBAKhQiFQn0N8uK7m3t6erBYLLq7WeWMpCUCY0yPiFwKPAtYgfuMMRtE5OLY/b8FfgKEgbtj9doeY0xpsmJSajD9G+SVlJTQ3t5OY2NjXy8kEelbhaTzCiobSaatuy4tLTWrV69OdRgqR8SXptbV1fU1yIsfvKMlJDWZurq6cDgcHHLIIWP6fhFZM9QbbV1Pp9QwBu5ubm1tpb6+nrq6Onp7e7WEpLKCJgKlEmSz2foa5B1wwAG0tbXtV0LSs5tVJtJEoNQY9G+QN336dDo6OmhpaaGmpobm5mYA3d2sMoYmAqUmQHx382AN8uIlJN3drNKV/lYqNcHiDfLy8vL2KyHp7maVjjQRKJVEAxvkdXR07HPGAmgJSaWeJgKlJkn/3c1FRUWDNsjTs5tVKmgiUCpF+jfI6+3tpa2trW93c3d3d18JyeFwpDpUleU0ESiVBqxW6z4N8trb2/vOWIif3ex2u3W/gkoKTQRKpRkRwePx4PF4+hrkNTQ0UFFRQVNTE1arFbfbrXMKasJoIlAqzTkcDoqKiigsLKStrY2amhpqamqIRCJ9O5+VGg9NBEpliP6b2IqLi2lsbKSyspKmpiYsFgtut1snmdWYaCJQKgNZrda+NtodHR3U1dVRWVlJT08Pdrsdl8ulexRUwjQRKJXhXC4X06dPZ+rUqbS0tFBdXU19fX3ffTrBrEaiiUCpLGGxWPrOVeju7qa+vn6f0pHH49EJZjUoTQRKZSG73b7PBHNtbS3V1dV9E8zaIVX1p4lAqSzWf4J5xowZNDY2UlVVRXNzc99OZ22Ep/Q3QKkcMXCCub6+noqKCtra2nSCOcdpIlAqB7lcLqZNm8bUqVNpbm7eZ4JZ21rkHk0ESuUwEdlngrmxsZG9e/fqBHOO0USglAKiE8wFBQWEw2Ha2tqoq6ujurqa3t5eHA4HTqdTS0dZShOBUmofA4/hbG5upqKiQieYs5j+ayqlhmS1WvtOW+vs7KSuro6qqira2tqw2Wy4XC4tHWUBTQRKqYQ4nc6+Ceb4Dua6ujqMMbhcLp1gzmCaCJRSoyIifWcnlJSU7NP8zmq16lxCkkQikaRdWxOBUmrM+k8wt7e3U1dX1zdKUBPP4/Ek5bqaCJRS49b/MJ3i4uJUh6NGSWd5lFIqx2kiUEqpHKeJQCmlcpwmAqWUynGaCJRSKsdpIlBKqRyniUAppXKcJgKllMpxkmk7AEWkGtgxxm8vAGomMJxky6R4MylWyKx4MylWyKx4MylWGF+8BxhjCge7I+MSwXiIyGpjTGmq40hUJsWbSbFCZsWbSbFCZsWbSbFC8uLV0pBSSuU4TQRKKZXjci0RLEt1AKOUSfFmUqyQWfFmUqyQWfFmUqyQpHhzao5AKaXU/nJtRKCUUmoATQRKKZXjciYRiMiZIrJZRLaIyNWpjmc4InKfiFSJyPpUxzISESkRkRdFpExENojI5amOaSgi4hKRt0TkvVis16c6pkSIiFVE3hWRJ1Mdy3BEZLuIrBORtSKyOtXxjERE8kTkERHZFPv9PSHVMQ1GRA6J/UzjH00icsWEPkcuzBGIiBX4ADgNKAfeBs4zxmxMaWBDEJGFQAvwR2PM4amOZzgiMg2YZox5R0T8wBrgnHT82Ur0IF2vMaZFROzAa8Dlxpg3UhzasETkSqAUCBhjFqc6nqGIyHag1BiTERu0ROQPwKvGmN+JiAPwGGMaUhzWsGKvZbuB44wxY91Yu59cGRF8AthijPnIGNMFPAQsTXFMQzLGvALUpTqORBhj9hpj3ol93gyUATNSG9XgTFRL7KY99pHW74REpBj4LPC7VMeSTUQkACwE7gUwxnSlexKIOQXYOpFJAHInEcwAdvW7XU6avlhlMhGZBRwFvJniUIYUK7OsBaqAfxhj0jbWmNuB7wGRFMeRCAM8JyJrROSiVAczgo8B1cD9sbLb70TEm+qgEnAu8OBEXzRXEoEM8rW0fieYaUTEBzwKXGGMaUp1PEMxxvQaY+YDxcAnRCRtS28ishioMsasSXUsCVpgjDkaWARcEitxpisbcDTwG2PMUUArkO5zhw5gCfDXib52riSCcqCk3+1iYE+KYsk6sXr7o8ADxpi/pTqeRMTKAC8BZ6Y2kmEtAJbEau8PASeLyJ9TG9LQjDF7Yv+tAh4jWpJNV+VAeb8R4SNEE0M6WwS8Y4ypnOgL50oieBuYIyKzY1n1XGBFimPKCrEJ2HuBMmPMramOZzgiUigiebHP3cCpwKaUBjUMY8wPjDHFxphZRH9n/2mM+bcUhzUoEfHGFgsQK7GcDqTtqjdjTAWwS0QOiX3pFCDtFjgMcB5JKAtBdHiU9YwxPSJyKfAsYAXuM8ZsSHFYQxKRB4FPAwUiUg5ca4y5N7VRDWkB8O/AuljtHeCHxpinUhfSkKYBf4itvLAADxtj0npJZgaZAjwWfV+ADfhfY8wzqQ1pRJcBD8TeHH4EfC3F8QxJRDxEVz1+KynXz4Xlo0oppYaWK6UhpZRSQ9BEoJRSOU4TgVJK5ThNBEopleM0ESilVI7TRKByVqz75Ldjn08XkUcm6LrXich3Y5//VEROnYjrKpUsunxU5axYb6QnJ7rDq4hcB7QYY341kddVKll0RKBy2S+AA2M93v8aP/9BRC4QkeUi8oSIbBORS0XkylhzsjdEJBR73IEi8kysydqrIjJ34BOIyO9F5Iuxz7eLyPUi8k6sb//c2Ne9sTMo3o49R9p2xlXZSROBymVXE23pOx+4asB9hwNfJdov57+BtlhzsteB/xN7zDLgMmPMMcB3gbsTeM6aWGO238S+B+BHRNtHHAt8Brg5QzphqiyREy0mlBqDF2PnKzSLSCPwROzr64CPx7qtngj8NdZWAcCZwHXjTfnWAJ+PfX460eZy8cTgAmYSPdtBqaTTRKDU4Dr7fR7pdztC9O/GAjTERhNjuW4v//r7E+ALxpjNYwtVqfHR0pDKZc2AfyzfGDtzYZuIfAmiXVhF5MgxxvEscFmskysictQYr6PUmGgiUDnLGFMLrIxNEt88hkucD3xDRN4DNjD2409vIHps5vuxWG4Y43WUGhNdPqqUUjlORwRKKZXjNBEopVSO00SglFI5ThOBUkrlOE0ESimV4zQRKKVUjtNEoJRSOe7/A9MroQhdgGA6AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ts = ci.index\n", "low, high = np.transpose(ci.values)\n", "\n", "plt.fill_between(ts, low, high, color='gray', alpha=0.3)\n", "kmf.survival_function_.plot(ax=plt.gca())\n", "plt.ylabel('Survival function');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Part of [Survival Analysis in Python](https://allendowney.github.io/SurvivalAnalysisPython/)\n", "\n", "Allen B. Downey\n", "\n", "[Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0)](https://creativecommons.org/licenses/by-nc-sa/4.0/)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Tags", "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.7.11" } }, "nbformat": 4, "nbformat_minor": 2 }