{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c7f4eb5a",
   "metadata": {},
   "source": [
    "\n",
    "# Chapter 11\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ae03bc9",
   "metadata": {},
   "source": [
    "# Lab: Survival Analysis\n",
    " In this lab, we perform survival analyses on three separate data\n",
    "sets. In  Section 11.8.1 we analyze the  `BrainCancer` \n",
    "data  that was first described in Section 11.3. In Section 11.8.2, we examine the  `Publication` \n",
    "data  from Section 11.5.4. Finally,  Section 11.8.3  explores\n",
    "a simulated call-center data set.\n",
    "\n",
    "We begin by importing some of our libraries at this top\n",
    "level. This makes the code more readable, as scanning the first few\n",
    "lines of the notebook tell us what libraries are used in this\n",
    "notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "91ac40fd",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:03.064349Z",
     "iopub.status.busy": "2023-07-26T05:18:03.064029Z",
     "iopub.status.idle": "2023-07-26T05:18:04.312445Z",
     "shell.execute_reply": "2023-07-26T05:18:04.312016Z"
    }
   },
   "outputs": [],
   "source": [
    "from matplotlib.pyplot import subplots\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from ISLP.models import ModelSpec as MS\n",
    "from ISLP import load_data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3dbcbbf",
   "metadata": {},
   "source": [
    "We  also collect the new imports\n",
    "needed for this lab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "99782418",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.315115Z",
     "iopub.status.busy": "2023-07-26T05:18:04.314859Z",
     "iopub.status.idle": "2023-07-26T05:18:04.379768Z",
     "shell.execute_reply": "2023-07-26T05:18:04.379303Z"
    }
   },
   "outputs": [],
   "source": [
    "from lifelines import \\\n",
    "     (KaplanMeierFitter,\n",
    "      CoxPHFitter)\n",
    "from lifelines.statistics import \\\n",
    "     (logrank_test,\n",
    "      multivariate_logrank_test)\n",
    "from ISLP.survival import sim_time\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c538d28",
   "metadata": {},
   "source": [
    "## Brain Cancer Data\n",
    "\n",
    "We begin with the `BrainCancer` data set, contained in the `ISLP` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3137149a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.382594Z",
     "iopub.status.busy": "2023-07-26T05:18:04.382386Z",
     "iopub.status.idle": "2023-07-26T05:18:04.389761Z",
     "shell.execute_reply": "2023-07-26T05:18:04.389375Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Index(['sex', 'diagnosis', 'loc', 'ki', 'gtv', 'stereo', 'status', 'time'], dtype='object')"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BrainCancer = load_data('BrainCancer')\n",
    "BrainCancer.columns\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e798f172",
   "metadata": {},
   "source": [
    "The rows index the 88 patients, while the 8 columns contain the predictors and outcome variables.\n",
    "We first briefly examine the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "45963c92",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.392289Z",
     "iopub.status.busy": "2023-07-26T05:18:04.392117Z",
     "iopub.status.idle": "2023-07-26T05:18:04.396285Z",
     "shell.execute_reply": "2023-07-26T05:18:04.395840Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Female    45\n",
       "Male      43\n",
       "Name: sex, dtype: int64"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BrainCancer['sex'].value_counts()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "73be61f6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.400201Z",
     "iopub.status.busy": "2023-07-26T05:18:04.399975Z",
     "iopub.status.idle": "2023-07-26T05:18:04.403967Z",
     "shell.execute_reply": "2023-07-26T05:18:04.403581Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Meningioma    42\n",
       "HG glioma     22\n",
       "Other         14\n",
       "LG glioma      9\n",
       "Name: diagnosis, dtype: int64"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BrainCancer['diagnosis'].value_counts()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "572f0b9e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.406232Z",
     "iopub.status.busy": "2023-07-26T05:18:04.406110Z",
     "iopub.status.idle": "2023-07-26T05:18:04.409773Z",
     "shell.execute_reply": "2023-07-26T05:18:04.409374Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0    53\n",
       "1    35\n",
       "Name: status, dtype: int64"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "BrainCancer['status'].value_counts()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbd132de",
   "metadata": {},
   "source": [
    "Before beginning an analysis, it is important to know how the\n",
    "`status` variable has been coded.  Most software\n",
    "uses the convention that a `status` of 1 indicates an\n",
    "uncensored observation (often death), and a `status` of 0 indicates a censored\n",
    "observation. But some scientists might use the opposite coding. For\n",
    "the  `BrainCancer`  data set 35 patients died before the end of\n",
    "the study, so we are using the conventional coding.\n",
    "\n",
    "To begin the analysis, we re-create  the Kaplan-Meier survival curve shown in Figure 11.2. The main\n",
    "package we will use for survival analysis\n",
    "is `lifelines`.\n",
    "The variable  `time`  corresponds to $y_i$, the time to the $i$th event (either censoring or\n",
    "death). The first argument to `km.fit` is the event time, and the\n",
    "second argument is the censoring variable, with a 1 indicating an observed\n",
    "failure time. The `plot()` method produces a survival curve with pointwise confidence\n",
    "intervals. By default, these are  90% confidence intervals, but this can be changed\n",
    "by setting the `alpha` argument to one minus the desired\n",
    "confidence level."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "92c39707",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.411840Z",
     "iopub.status.busy": "2023-07-26T05:18:04.411663Z",
     "iopub.status.idle": "2023-07-26T05:18:04.539774Z",
     "shell.execute_reply": "2023-07-26T05:18:04.539439Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='timeline'>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 25648 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots(figsize=(8,8))\n",
    "km = KaplanMeierFitter()\n",
    "km_brain = km.fit(BrainCancer['time'], BrainCancer['status'])\n",
    "km_brain.plot(label='Kaplan Meier estimate', ax=ax)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f037665b",
   "metadata": {},
   "source": [
    "Next we create Kaplan-Meier survival curves that are stratified by\n",
    "`sex`, in order to reproduce  Figure 11.3.\n",
    "We  do this using the `groupby()` method of  a  dataframe.\n",
    "This  method returns a generator that can\n",
    "be  iterated over in the `for` loop. In this case,\n",
    "the  items  in the `for` loop are  2-tuples  representing\n",
    "the groups: the first entry is the  value\n",
    "of the grouping column  `sex` while the  second  value\n",
    "is the dataframe consisting of all  rows in the\n",
    "dataframe matching that value of `sex`.\n",
    "We will want to use this data below\n",
    "in the log-rank test,  hence we store this\n",
    "information in the dictionary `by_sex`. Finally,\n",
    "we have also used the notion of\n",
    " *string interpolation* to automatically\n",
    "label the  different lines in the  plot. String\n",
    "interpolation is a powerful technique to format strings ---\n",
    "`Python` has many ways to facilitate such operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3fc7848c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.542265Z",
     "iopub.status.busy": "2023-07-26T05:18:04.542002Z",
     "iopub.status.idle": "2023-07-26T05:18:04.674613Z",
     "shell.execute_reply": "2023-07-26T05:18:04.674126Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 30140 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots(figsize=(8,8))\n",
    "by_sex = {}\n",
    "for sex, df in BrainCancer.groupby('sex'):\n",
    "    by_sex[sex] = df\n",
    "    km_sex = km.fit(df['time'], df['status'])\n",
    "    km_sex.plot(label='Sex=%s' % sex, ax=ax)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0c1a16a",
   "metadata": {},
   "source": [
    "As discussed in Section 11.4, we can perform a\n",
    "log-rank test to compare the survival of males to females. We use\n",
    "the `logrank_test()` function from the `lifelines.statistics` module.\n",
    "The first two arguments are the event times, with the second\n",
    "denoting the corresponding (optional) censoring indicators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "bf30d26f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.679625Z",
     "iopub.status.busy": "2023-07-26T05:18:04.679410Z",
     "iopub.status.idle": "2023-07-26T05:18:04.743501Z",
     "shell.execute_reply": "2023-07-26T05:18:04.743166Z"
    },
    "lines_to_next_cell": 2
   },
   "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>t_0</th>\n",
       "      <td>-1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>null_distribution</th>\n",
       "      <td>chi squared</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>degrees_of_freedom</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>test_name</th>\n",
       "      <td>logrank_test</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>test_statistic</th>\n",
       "      <th>p</th>\n",
       "      <th>-log2(p)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.44</td>\n",
       "      <td>0.23</td>\n",
       "      <td>2.12</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>"
      ],
      "text/latex": [
       "\\begin{tabular}{lrrr}\n",
       " & test_statistic & p & -log2(p) \\\\\n",
       "0 & 1.44 & 0.23 & 2.12 \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "<lifelines.StatisticalResult: logrank_test>\n",
       "               t_0 = -1\n",
       " null_distribution = chi squared\n",
       "degrees_of_freedom = 1\n",
       "         test_name = logrank_test\n",
       "\n",
       "---\n",
       " test_statistic    p  -log2(p)\n",
       "           1.44 0.23      2.12"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logrank_test(by_sex['Male']['time'],\n",
    "             by_sex['Female']['time'],\n",
    "             by_sex['Male']['status'],\n",
    "             by_sex['Female']['status'])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e270649c",
   "metadata": {},
   "source": [
    "The resulting $p$-value is $0.23$, indicating no evidence of a\n",
    "difference in survival between the two sexes.\n",
    "\n",
    "Next, we  use the `CoxPHFitter()`  estimator\n",
    "from `lifelines` to fit Cox proportional hazards models.\n",
    "To begin, we consider a model that uses  `sex`  as the only predictor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2ab78e07",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.746088Z",
     "iopub.status.busy": "2023-07-26T05:18:04.745814Z",
     "iopub.status.idle": "2023-07-26T05:18:04.773972Z",
     "shell.execute_reply": "2023-07-26T05:18:04.773498Z"
    }
   },
   "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>coef</th>\n",
       "      <th>se(coef)</th>\n",
       "      <th>p</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>covariate</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>sex[Male]</th>\n",
       "      <td>0.407668</td>\n",
       "      <td>0.342004</td>\n",
       "      <td>0.233262</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               coef  se(coef)         p\n",
       "covariate                              \n",
       "sex[Male]  0.407668  0.342004  0.233262"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "coxph = CoxPHFitter # shorthand\n",
    "sex_df = BrainCancer[['time', 'status', 'sex']]\n",
    "model_df = MS(['time', 'status', 'sex'],\n",
    "              intercept=False).fit_transform(sex_df)\n",
    "cox_fit = coxph().fit(model_df,\n",
    "                      'time',\n",
    "                      'status')\n",
    "cox_fit.summary[['coef', 'se(coef)', 'p']]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b58b93ae",
   "metadata": {},
   "source": [
    "The first argument to `fit` should be a data frame containing\n",
    "at least the event time (the second argument `time` in this case),\n",
    "as well as an optional censoring variable (the argument `status` in this case).\n",
    "Note also that the Cox model does not include an intercept, which is why\n",
    "we used the `intercept=False` argument to `ModelSpec` above.\n",
    "The `summary()` method delivers many columns; we chose to abbreviate its output here.\n",
    "It is possible to obtain the likelihood ratio test comparing this model to the one\n",
    "with no features as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "4716b7b0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.776775Z",
     "iopub.status.busy": "2023-07-26T05:18:04.776463Z",
     "iopub.status.idle": "2023-07-26T05:18:04.783362Z",
     "shell.execute_reply": "2023-07-26T05:18:04.782865Z"
    }
   },
   "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>null_distribution</th>\n",
       "      <td>chi squared</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>degrees_freedom</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>test_name</th>\n",
       "      <td>log-likelihood ratio test</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>test_statistic</th>\n",
       "      <th>p</th>\n",
       "      <th>-log2(p)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.44</td>\n",
       "      <td>0.23</td>\n",
       "      <td>2.12</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>"
      ],
      "text/latex": [
       "\\begin{tabular}{lrrr}\n",
       " & test_statistic & p & -log2(p) \\\\\n",
       "0 & 1.44 & 0.23 & 2.12 \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "<lifelines.StatisticalResult: log-likelihood ratio test>\n",
       "null_distribution = chi squared\n",
       "  degrees_freedom = 1\n",
       "        test_name = log-likelihood ratio test\n",
       "\n",
       "---\n",
       " test_statistic    p  -log2(p)\n",
       "           1.44 0.23      2.12"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cox_fit.log_likelihood_ratio_test()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2820f486",
   "metadata": {},
   "source": [
    "Regardless of which test we use, we see that there is no clear\n",
    "evidence for a difference in survival between males and females.  As\n",
    "we learned in this chapter, the score test from the Cox model is\n",
    "exactly equal to the log rank test statistic!\n",
    "\n",
    "Now we fit a  model that makes use of additional predictors. We first note\n",
    "that one of our `diagnosis` values is missing, hence\n",
    "we drop that observation before continuing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "c2767d88",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.785790Z",
     "iopub.status.busy": "2023-07-26T05:18:04.785586Z",
     "iopub.status.idle": "2023-07-26T05:18:04.824235Z",
     "shell.execute_reply": "2023-07-26T05:18:04.823739Z"
    }
   },
   "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>coef</th>\n",
       "      <th>se(coef)</th>\n",
       "      <th>p</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>covariate</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>sex[Male]</th>\n",
       "      <td>0.183748</td>\n",
       "      <td>0.360358</td>\n",
       "      <td>0.610119</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>diagnosis[LG glioma]</th>\n",
       "      <td>-1.239530</td>\n",
       "      <td>0.579555</td>\n",
       "      <td>0.032455</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>diagnosis[Meningioma]</th>\n",
       "      <td>-2.154566</td>\n",
       "      <td>0.450524</td>\n",
       "      <td>0.000002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>diagnosis[Other]</th>\n",
       "      <td>-1.268870</td>\n",
       "      <td>0.617672</td>\n",
       "      <td>0.039949</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>loc[Supratentorial]</th>\n",
       "      <td>0.441195</td>\n",
       "      <td>0.703669</td>\n",
       "      <td>0.530665</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>ki</th>\n",
       "      <td>-0.054955</td>\n",
       "      <td>0.018314</td>\n",
       "      <td>0.002693</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gtv</th>\n",
       "      <td>0.034293</td>\n",
       "      <td>0.022333</td>\n",
       "      <td>0.124661</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>stereo[SRT]</th>\n",
       "      <td>0.177778</td>\n",
       "      <td>0.601578</td>\n",
       "      <td>0.767597</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                           coef  se(coef)         p\n",
       "covariate                                          \n",
       "sex[Male]              0.183748  0.360358  0.610119\n",
       "diagnosis[LG glioma]  -1.239530  0.579555  0.032455\n",
       "diagnosis[Meningioma] -2.154566  0.450524  0.000002\n",
       "diagnosis[Other]      -1.268870  0.617672  0.039949\n",
       "loc[Supratentorial]    0.441195  0.703669  0.530665\n",
       "ki                    -0.054955  0.018314  0.002693\n",
       "gtv                    0.034293  0.022333  0.124661\n",
       "stereo[SRT]            0.177778  0.601578  0.767597"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cleaned = BrainCancer.dropna()\n",
    "all_MS = MS(cleaned.columns, intercept=False)\n",
    "all_df = all_MS.fit_transform(cleaned)\n",
    "fit_all = coxph().fit(all_df,\n",
    "                      'time',\n",
    "                      'status')\n",
    "fit_all.summary[['coef', 'se(coef)', 'p']]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eee4ab1f",
   "metadata": {},
   "source": [
    " The `diagnosis` variable has been coded so that the baseline\n",
    "corresponds to HG glioma. The results indicate that the risk associated with HG glioma\n",
    "is more than eight times (i.e. $e^{2.15}=8.62$) the risk associated\n",
    "with meningioma. In other words, after adjusting for the other\n",
    "predictors, patients with HG glioma have much worse survival compared\n",
    "to those with meningioma.  In addition, larger values of the Karnofsky\n",
    "index, `ki`, are associated with lower risk, i.e. longer survival.\n",
    "\n",
    "Finally, we plot estimated survival curves for each diagnosis category,\n",
    "adjusting for the other predictors.  To make these plots, we set the\n",
    "values of the other predictors equal to the mean for quantitative variables\n",
    "and equal to the mode for categorical. To do this, we use the\n",
    "`apply()` method across rows (i.e. `axis=0`) with a function\n",
    "`representative` that checks if a column is categorical\n",
    "or not."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "ede1d219",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.826852Z",
     "iopub.status.busy": "2023-07-26T05:18:04.826662Z",
     "iopub.status.idle": "2023-07-26T05:18:04.831013Z",
     "shell.execute_reply": "2023-07-26T05:18:04.830464Z"
    }
   },
   "outputs": [],
   "source": [
    "levels = cleaned['diagnosis'].unique()\n",
    "def representative(series):\n",
    "    if hasattr(series.dtype, 'categories'):\n",
    "        return pd.Series.mode(series)\n",
    "    else:\n",
    "        return series.mean()\n",
    "modal_data = cleaned.apply(representative, axis=0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1c307ae",
   "metadata": {},
   "source": [
    "We make four\n",
    "copies of the column means and assign the `diagnosis` column to be the four different\n",
    "diagnoses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "dc032a71",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.833522Z",
     "iopub.status.busy": "2023-07-26T05:18:04.833389Z",
     "iopub.status.idle": "2023-07-26T05:18:04.840549Z",
     "shell.execute_reply": "2023-07-26T05:18:04.840260Z"
    }
   },
   "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>sex</th>\n",
       "      <th>diagnosis</th>\n",
       "      <th>loc</th>\n",
       "      <th>ki</th>\n",
       "      <th>gtv</th>\n",
       "      <th>stereo</th>\n",
       "      <th>status</th>\n",
       "      <th>time</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Female</td>\n",
       "      <td>Meningioma</td>\n",
       "      <td>Supratentorial</td>\n",
       "      <td>80.91954</td>\n",
       "      <td>8.687011</td>\n",
       "      <td>SRT</td>\n",
       "      <td>0.402299</td>\n",
       "      <td>27.188621</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Female</td>\n",
       "      <td>HG glioma</td>\n",
       "      <td>Supratentorial</td>\n",
       "      <td>80.91954</td>\n",
       "      <td>8.687011</td>\n",
       "      <td>SRT</td>\n",
       "      <td>0.402299</td>\n",
       "      <td>27.188621</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Female</td>\n",
       "      <td>LG glioma</td>\n",
       "      <td>Supratentorial</td>\n",
       "      <td>80.91954</td>\n",
       "      <td>8.687011</td>\n",
       "      <td>SRT</td>\n",
       "      <td>0.402299</td>\n",
       "      <td>27.188621</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Female</td>\n",
       "      <td>Other</td>\n",
       "      <td>Supratentorial</td>\n",
       "      <td>80.91954</td>\n",
       "      <td>8.687011</td>\n",
       "      <td>SRT</td>\n",
       "      <td>0.402299</td>\n",
       "      <td>27.188621</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      sex   diagnosis             loc        ki       gtv stereo    status  \\\n",
       "0  Female  Meningioma  Supratentorial  80.91954  8.687011    SRT  0.402299   \n",
       "0  Female   HG glioma  Supratentorial  80.91954  8.687011    SRT  0.402299   \n",
       "0  Female   LG glioma  Supratentorial  80.91954  8.687011    SRT  0.402299   \n",
       "0  Female       Other  Supratentorial  80.91954  8.687011    SRT  0.402299   \n",
       "\n",
       "        time  \n",
       "0  27.188621  \n",
       "0  27.188621  \n",
       "0  27.188621  \n",
       "0  27.188621  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "modal_df = pd.DataFrame(\n",
    "              [modal_data.iloc[0] for _ in range(len(levels))])\n",
    "modal_df['diagnosis'] = levels\n",
    "modal_df\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84da2586",
   "metadata": {},
   "source": [
    "We then construct the model matrix based on the model specification `all_MS` used to fit\n",
    "the model, and name the rows according to the levels of `diagnosis`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e7c1fe43",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.843157Z",
     "iopub.status.busy": "2023-07-26T05:18:04.843004Z",
     "iopub.status.idle": "2023-07-26T05:18:04.852876Z",
     "shell.execute_reply": "2023-07-26T05:18:04.852545Z"
    }
   },
   "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>sex[Male]</th>\n",
       "      <th>diagnosis[LG glioma]</th>\n",
       "      <th>diagnosis[Meningioma]</th>\n",
       "      <th>diagnosis[Other]</th>\n",
       "      <th>loc[Supratentorial]</th>\n",
       "      <th>ki</th>\n",
       "      <th>gtv</th>\n",
       "      <th>stereo[SRT]</th>\n",
       "      <th>status</th>\n",
       "      <th>time</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Meningioma</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>80.91954</td>\n",
       "      <td>8.687011</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.402299</td>\n",
       "      <td>27.188621</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>HG glioma</th>\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",
       "      <td>80.91954</td>\n",
       "      <td>8.687011</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.402299</td>\n",
       "      <td>27.188621</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>LG glioma</th>\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>80.91954</td>\n",
       "      <td>8.687011</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.402299</td>\n",
       "      <td>27.188621</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Other</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>80.91954</td>\n",
       "      <td>8.687011</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.402299</td>\n",
       "      <td>27.188621</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            sex[Male]  diagnosis[LG glioma]  diagnosis[Meningioma]  \\\n",
       "Meningioma        0.0                   0.0                    1.0   \n",
       "HG glioma         0.0                   0.0                    0.0   \n",
       "LG glioma         0.0                   1.0                    0.0   \n",
       "Other             0.0                   0.0                    0.0   \n",
       "\n",
       "            diagnosis[Other]  loc[Supratentorial]        ki       gtv  \\\n",
       "Meningioma               0.0                  1.0  80.91954  8.687011   \n",
       "HG glioma                0.0                  1.0  80.91954  8.687011   \n",
       "LG glioma                0.0                  1.0  80.91954  8.687011   \n",
       "Other                    1.0                  1.0  80.91954  8.687011   \n",
       "\n",
       "            stereo[SRT]    status       time  \n",
       "Meningioma          1.0  0.402299  27.188621  \n",
       "HG glioma           1.0  0.402299  27.188621  \n",
       "LG glioma           1.0  0.402299  27.188621  \n",
       "Other               1.0  0.402299  27.188621  "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "modal_X = all_MS.transform(modal_df)\n",
    "modal_X.index = levels\n",
    "modal_X\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cfe1ec4",
   "metadata": {},
   "source": [
    "We can use the `predict_survival_function()` method to obtain the estimated survival function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f89fbed7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.855400Z",
     "iopub.status.busy": "2023-07-26T05:18:04.855197Z",
     "iopub.status.idle": "2023-07-26T05:18:04.864584Z",
     "shell.execute_reply": "2023-07-26T05:18:04.864194Z"
    },
    "lines_to_next_cell": 0
   },
   "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>Meningioma</th>\n",
       "      <th>HG glioma</th>\n",
       "      <th>LG glioma</th>\n",
       "      <th>Other</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.07</th>\n",
       "      <td>0.997947</td>\n",
       "      <td>0.982430</td>\n",
       "      <td>0.994881</td>\n",
       "      <td>0.995029</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.18</th>\n",
       "      <td>0.997947</td>\n",
       "      <td>0.982430</td>\n",
       "      <td>0.994881</td>\n",
       "      <td>0.995029</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.41</th>\n",
       "      <td>0.995679</td>\n",
       "      <td>0.963342</td>\n",
       "      <td>0.989245</td>\n",
       "      <td>0.989555</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.54</th>\n",
       "      <td>0.995679</td>\n",
       "      <td>0.963342</td>\n",
       "      <td>0.989245</td>\n",
       "      <td>0.989555</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2.03</th>\n",
       "      <td>0.995679</td>\n",
       "      <td>0.963342</td>\n",
       "      <td>0.989245</td>\n",
       "      <td>0.989555</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>65.02</th>\n",
       "      <td>0.688772</td>\n",
       "      <td>0.040136</td>\n",
       "      <td>0.394181</td>\n",
       "      <td>0.404936</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>67.38</th>\n",
       "      <td>0.688772</td>\n",
       "      <td>0.040136</td>\n",
       "      <td>0.394181</td>\n",
       "      <td>0.404936</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>73.74</th>\n",
       "      <td>0.688772</td>\n",
       "      <td>0.040136</td>\n",
       "      <td>0.394181</td>\n",
       "      <td>0.404936</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>78.75</th>\n",
       "      <td>0.688772</td>\n",
       "      <td>0.040136</td>\n",
       "      <td>0.394181</td>\n",
       "      <td>0.404936</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>82.56</th>\n",
       "      <td>0.688772</td>\n",
       "      <td>0.040136</td>\n",
       "      <td>0.394181</td>\n",
       "      <td>0.404936</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>85 rows × 4 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       Meningioma  HG glioma  LG glioma     Other\n",
       "0.07     0.997947   0.982430   0.994881  0.995029\n",
       "1.18     0.997947   0.982430   0.994881  0.995029\n",
       "1.41     0.995679   0.963342   0.989245  0.989555\n",
       "1.54     0.995679   0.963342   0.989245  0.989555\n",
       "2.03     0.995679   0.963342   0.989245  0.989555\n",
       "...           ...        ...        ...       ...\n",
       "65.02    0.688772   0.040136   0.394181  0.404936\n",
       "67.38    0.688772   0.040136   0.394181  0.404936\n",
       "73.74    0.688772   0.040136   0.394181  0.404936\n",
       "78.75    0.688772   0.040136   0.394181  0.404936\n",
       "82.56    0.688772   0.040136   0.394181  0.404936\n",
       "\n",
       "[85 rows x 4 columns]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predicted_survival = fit_all.predict_survival_function(modal_X)\n",
    "predicted_survival\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29afd641",
   "metadata": {},
   "source": [
    "This returns a data frame,\n",
    "whose plot methods yields the different survival curves. To avoid clutter in\n",
    "the plots, we do not display confidence intervals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "8f0329b4",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.867553Z",
     "iopub.status.busy": "2023-07-26T05:18:04.867105Z",
     "iopub.status.idle": "2023-07-26T05:18:04.987503Z",
     "shell.execute_reply": "2023-07-26T05:18:04.986681Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 47476 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots(figsize=(8, 8))\n",
    "predicted_survival.plot(ax=ax);\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12723ce5",
   "metadata": {},
   "source": [
    "## Publication Data\n",
    "The  `Publication`  data   presented in Section 11.5.4  can be\n",
    "found in the `ISLP` package.\n",
    "We first reproduce Figure 11.5  by plotting the Kaplan-Meier curves\n",
    "stratified on the  `posres`  variable, which records whether the\n",
    "study had a positive or negative result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "3045bfc0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:04.990450Z",
     "iopub.status.busy": "2023-07-26T05:18:04.990225Z",
     "iopub.status.idle": "2023-07-26T05:18:05.181343Z",
     "shell.execute_reply": "2023-07-26T05:18:05.178676Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 38840 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots(figsize=(8,8))\n",
    "Publication = load_data('Publication')\n",
    "by_result = {}\n",
    "for result, df in Publication.groupby('posres'):\n",
    "    by_result[result] = df\n",
    "    km_result = km.fit(df['time'], df['status'])\n",
    "    km_result.plot(label='Result=%d' % result, ax=ax)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fcb22f7",
   "metadata": {},
   "source": [
    "As discussed previously, the $p$-values from fitting Cox’s\n",
    "proportional hazards model to the `posres` variable are quite\n",
    "large, providing no evidence of a difference in time-to-publication\n",
    "between studies with positive versus negative results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "d070f716",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:05.184743Z",
     "iopub.status.busy": "2023-07-26T05:18:05.184569Z",
     "iopub.status.idle": "2023-07-26T05:18:05.222538Z",
     "shell.execute_reply": "2023-07-26T05:18:05.221996Z"
    },
    "lines_to_next_cell": 2
   },
   "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>coef</th>\n",
       "      <th>se(coef)</th>\n",
       "      <th>p</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>covariate</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>posres</th>\n",
       "      <td>0.148076</td>\n",
       "      <td>0.161625</td>\n",
       "      <td>0.359579</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               coef  se(coef)         p\n",
       "covariate                              \n",
       "posres     0.148076  0.161625  0.359579"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "posres_df = MS(['posres',\n",
    "                'time',\n",
    "                'status'],\n",
    "                intercept=False).fit_transform(Publication)\n",
    "posres_fit = coxph().fit(posres_df,\n",
    "                         'time',\n",
    "                         'status')\n",
    "posres_fit.summary[['coef', 'se(coef)', 'p']]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "513a55b1",
   "metadata": {},
   "source": [
    "However, the results change dramatically when we include other\n",
    "predictors in the model. Here we exclude the funding mechanism\n",
    "variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "2bbcdd0c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:05.225196Z",
     "iopub.status.busy": "2023-07-26T05:18:05.224969Z",
     "iopub.status.idle": "2023-07-26T05:18:05.270258Z",
     "shell.execute_reply": "2023-07-26T05:18:05.269721Z"
    }
   },
   "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>coef</th>\n",
       "      <th>se(coef)</th>\n",
       "      <th>p</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>covariate</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>posres</th>\n",
       "      <td>0.570773</td>\n",
       "      <td>0.175960</td>\n",
       "      <td>1.179610e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>multi</th>\n",
       "      <td>-0.040860</td>\n",
       "      <td>0.251194</td>\n",
       "      <td>8.707842e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>clinend</th>\n",
       "      <td>0.546183</td>\n",
       "      <td>0.262000</td>\n",
       "      <td>3.709944e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>sampsize</th>\n",
       "      <td>0.000005</td>\n",
       "      <td>0.000015</td>\n",
       "      <td>7.507005e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>budget</th>\n",
       "      <td>0.004386</td>\n",
       "      <td>0.002465</td>\n",
       "      <td>7.515984e-02</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>impact</th>\n",
       "      <td>0.058318</td>\n",
       "      <td>0.006676</td>\n",
       "      <td>2.426306e-18</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               coef  se(coef)             p\n",
       "covariate                                  \n",
       "posres     0.570773  0.175960  1.179610e-03\n",
       "multi     -0.040860  0.251194  8.707842e-01\n",
       "clinend    0.546183  0.262000  3.709944e-02\n",
       "sampsize   0.000005  0.000015  7.507005e-01\n",
       "budget     0.004386  0.002465  7.515984e-02\n",
       "impact     0.058318  0.006676  2.426306e-18"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = MS(Publication.columns.drop('mech'),\n",
    "           intercept=False)\n",
    "coxph().fit(model.fit_transform(Publication),\n",
    "            'time',\n",
    "            'status').summary[['coef', 'se(coef)', 'p']]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75bb8aa6",
   "metadata": {},
   "source": [
    "We see that there are a number of statistically significant variables,\n",
    "including whether the trial focused on a clinical endpoint, the impact\n",
    "of the study, and whether the study had positive or negative results."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfe236e5",
   "metadata": {},
   "source": [
    "## Call Center Data\n",
    "\n",
    "In this section, we will simulate survival data using the relationship\n",
    "between cumulative hazard and\n",
    "the survival function explored in Exercise 8.\n",
    "Our simulated data will represent the observed\n",
    "wait times (in seconds) for 2,000 customers who have phoned a call\n",
    "center.  In this context, censoring occurs if a customer hangs up\n",
    "before his or her call is answered.\n",
    "\n",
    "There are three covariates: `Operators` (the number of call\n",
    "center operators available at the time of the call, which can range\n",
    "from $5$ to $15$), `Center` (either A, B, or C), and\n",
    "`Time` of day (Morning, Afternoon, or Evening). We generate data\n",
    "for these covariates so that all possibilities are equally likely: for\n",
    "instance, morning, afternoon and evening calls are equally likely, and\n",
    "any number of operators from $5$ to $15$ is equally likely. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "b8ece43a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:05.272852Z",
     "iopub.status.busy": "2023-07-26T05:18:05.272613Z",
     "iopub.status.idle": "2023-07-26T05:18:05.277670Z",
     "shell.execute_reply": "2023-07-26T05:18:05.277169Z"
    }
   },
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(10)\n",
    "N = 2000\n",
    "Operators = rng.choice(np.arange(5, 16),\n",
    "                       N,\n",
    "                       replace=True)\n",
    "Center = rng.choice(['A', 'B', 'C'],\n",
    "                    N,\n",
    "                    replace=True)\n",
    "Time = rng.choice(['Morn.', 'After.', 'Even.'],\n",
    "                   N,\n",
    "                   replace=True)\n",
    "D = pd.DataFrame({'Operators': Operators,\n",
    "                  'Center': pd.Categorical(Center),\n",
    "                  'Time': pd.Categorical(Time)})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c93e44f3",
   "metadata": {},
   "source": [
    "We then build a model matrix (omitting the intercept)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "3e4f766f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:05.279951Z",
     "iopub.status.busy": "2023-07-26T05:18:05.279774Z",
     "iopub.status.idle": "2023-07-26T05:18:05.289519Z",
     "shell.execute_reply": "2023-07-26T05:18:05.288948Z"
    }
   },
   "outputs": [],
   "source": [
    "model = MS(['Operators',\n",
    "            'Center',\n",
    "            'Time'],\n",
    "           intercept=False)\n",
    "X = model.fit_transform(D)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cad1ed19",
   "metadata": {},
   "source": [
    "It is worthwhile to take a peek at the model matrix `X`, so\n",
    "that we can be sure that we understand how the variables have been coded. By default,\n",
    "the levels of categorical variables are sorted and, as usual, the first column of the one-hot encoding\n",
    "of the variable is dropped."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "72f42d14",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:05.292713Z",
     "iopub.status.busy": "2023-07-26T05:18:05.292424Z",
     "iopub.status.idle": "2023-07-26T05:18:05.299084Z",
     "shell.execute_reply": "2023-07-26T05:18:05.298414Z"
    }
   },
   "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>Operators</th>\n",
       "      <th>Center[B]</th>\n",
       "      <th>Center[C]</th>\n",
       "      <th>Time[Even.]</th>\n",
       "      <th>Time[Morn.]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>13</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>15</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>2</th>\n",
       "      <td>7</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>3</th>\n",
       "      <td>7</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>13</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Operators  Center[B]  Center[C]  Time[Even.]  Time[Morn.]\n",
       "0         13        0.0        1.0          0.0          0.0\n",
       "1         15        0.0        0.0          1.0          0.0\n",
       "2          7        1.0        0.0          0.0          1.0\n",
       "3          7        0.0        1.0          0.0          1.0\n",
       "4         13        0.0        1.0          1.0          0.0"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[:5]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38c40ae1",
   "metadata": {},
   "source": [
    "Next,  we specify the coefficients and the hazard function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "8b921536",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:05.301800Z",
     "iopub.status.busy": "2023-07-26T05:18:05.301580Z",
     "iopub.status.idle": "2023-07-26T05:18:05.340807Z",
     "shell.execute_reply": "2023-07-26T05:18:05.311843Z"
    }
   },
   "outputs": [],
   "source": [
    "true_beta = np.array([0.04, -0.3, 0, 0.2, -0.2])\n",
    "true_linpred = X.dot(true_beta)\n",
    "hazard = lambda t: 1e-5 * t\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0698ffd",
   "metadata": {},
   "source": [
    "Here, we have set the coefficient associated with `Operators` to\n",
    "equal $0.04$; in other words, each additional operator leads to a\n",
    "$e^{0.04}=1.041$-fold increase in the “risk” that the call will be\n",
    "answered, given the `Center` and `Time` covariates. This\n",
    "makes sense: the greater the number of operators at hand, the shorter\n",
    "the wait time! The coefficient associated with `Center == B` is\n",
    "$-0.3$, and `Center == A` is treated as the baseline. This means\n",
    "that the risk of a call being answered at Center B is 0.74 times the\n",
    "risk that it will be answered at Center A; in other words, the wait\n",
    "times are a bit longer at Center B.\n",
    "\n",
    "Recall from Section 2.3.7 the use of `lambda`\n",
    "for creating short functions on the fly.\n",
    "We use the function\n",
    "`sim_time()` from the `ISLP.survival` package. This function\n",
    "uses the relationship between the survival function\n",
    "and cumulative hazard $S(t) = \\exp(-H(t))$ and the specific\n",
    "form of the cumulative hazard function in the Cox model\n",
    "to simulate data based on values of the linear predictor\n",
    "`true_linpred` and the cumulative hazard. \n",
    " We need to provide the cumulative hazard function, which we do here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "96ce0f99",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:05.399214Z",
     "iopub.status.busy": "2023-07-26T05:18:05.386941Z",
     "iopub.status.idle": "2023-07-26T05:18:05.423513Z",
     "shell.execute_reply": "2023-07-26T05:18:05.417411Z"
    },
    "lines_to_next_cell": 0
   },
   "outputs": [],
   "source": [
    "cum_hazard = lambda t: 1e-5 * t**2 / 2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1956e4c2",
   "metadata": {},
   "source": [
    "We are now ready to generate data under the Cox proportional hazards\n",
    "model. We truncate the maximum time to 1000 seconds to keep\n",
    "simulated wait times reasonable. The function\n",
    "`sim_time()` takes a linear predictor,\n",
    "a cumulative hazard function and a\n",
    "random number generator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "63d78ff9",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:05.450050Z",
     "iopub.status.busy": "2023-07-26T05:18:05.449814Z",
     "iopub.status.idle": "2023-07-26T05:18:05.587170Z",
     "shell.execute_reply": "2023-07-26T05:18:05.586673Z"
    }
   },
   "outputs": [],
   "source": [
    "W = np.array([sim_time(l, cum_hazard, rng)\n",
    "              for l in true_linpred])\n",
    "D['Wait time'] = np.clip(W, 0, 1000)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "035e4ecf",
   "metadata": {},
   "source": [
    "We now simulate our censoring variable, for which we assume\n",
    "90% of calls were answered (`Failed==1`) before the\n",
    "customer hung up (`Failed==0`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "fe008dbf",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:05.590252Z",
     "iopub.status.busy": "2023-07-26T05:18:05.589902Z",
     "iopub.status.idle": "2023-07-26T05:18:05.596018Z",
     "shell.execute_reply": "2023-07-26T05:18:05.595604Z"
    },
    "lines_to_next_cell": 2
   },
   "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>Operators</th>\n",
       "      <th>Center</th>\n",
       "      <th>Time</th>\n",
       "      <th>Wait time</th>\n",
       "      <th>Failed</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>13</td>\n",
       "      <td>C</td>\n",
       "      <td>After.</td>\n",
       "      <td>525.064979</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>15</td>\n",
       "      <td>A</td>\n",
       "      <td>Even.</td>\n",
       "      <td>254.677835</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>7</td>\n",
       "      <td>B</td>\n",
       "      <td>Morn.</td>\n",
       "      <td>487.739224</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>7</td>\n",
       "      <td>C</td>\n",
       "      <td>Morn.</td>\n",
       "      <td>308.580292</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>13</td>\n",
       "      <td>C</td>\n",
       "      <td>Even.</td>\n",
       "      <td>154.174608</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Operators Center    Time   Wait time  Failed\n",
       "0         13      C  After.  525.064979       1\n",
       "1         15      A   Even.  254.677835       1\n",
       "2          7      B   Morn.  487.739224       1\n",
       "3          7      C   Morn.  308.580292       1\n",
       "4         13      C   Even.  154.174608       1"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D['Failed'] = rng.choice([1, 0],\n",
    "                         N,\n",
    "                         p=[0.9, 0.1])\n",
    "D[:5]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "c3a2bec7",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:05.598122Z",
     "iopub.status.busy": "2023-07-26T05:18:05.597969Z",
     "iopub.status.idle": "2023-07-26T05:18:05.601009Z",
     "shell.execute_reply": "2023-07-26T05:18:05.600579Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9075"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D['Failed'].mean()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "207937e5",
   "metadata": {},
   "source": [
    "We now plot  Kaplan-Meier survival curves. First, we stratify by `Center`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "2b27af56",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:05.603318Z",
     "iopub.status.busy": "2023-07-26T05:18:05.603142Z",
     "iopub.status.idle": "2023-07-26T05:18:05.827076Z",
     "shell.execute_reply": "2023-07-26T05:18:05.826151Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Probability of Still Being on Hold')"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 102144 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots(figsize=(8,8))\n",
    "by_center = {}\n",
    "for center, df in D.groupby('Center'):\n",
    "    by_center[center] = df\n",
    "    km_center = km.fit(df['Wait time'], df['Failed'])\n",
    "    km_center.plot(label='Center=%s' % center, ax=ax)\n",
    "ax.set_title(\"Probability of Still Being on Hold\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be6d37f7",
   "metadata": {},
   "source": [
    "Next, we stratify by `Time`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "9625598d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:05.830554Z",
     "iopub.status.busy": "2023-07-26T05:18:05.830059Z",
     "iopub.status.idle": "2023-07-26T05:18:06.008529Z",
     "shell.execute_reply": "2023-07-26T05:18:06.008079Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Probability of Still Being on Hold')"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 108004 characters]",
      "text/plain": [
       "<Figure size 800x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots(figsize=(8,8))\n",
    "by_time = {}\n",
    "for time, df in D.groupby('Time'):\n",
    "    by_time[time] = df\n",
    "    km_time = km.fit(df['Wait time'], df['Failed'])\n",
    "    km_time.plot(label='Time=%s' % time, ax=ax)\n",
    "ax.set_title(\"Probability of Still Being on Hold\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1408ebc0",
   "metadata": {},
   "source": [
    "It seems that calls at Call Center B take longer to be answered than\n",
    "calls at Centers A and C. Similarly, it appears that wait times are\n",
    "longest in the morning and shortest in the evening hours. We can use a\n",
    "log-rank test to determine whether these differences are statistically\n",
    "significant using the function `multivariate_logrank_test()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "75a744ef",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:06.012649Z",
     "iopub.status.busy": "2023-07-26T05:18:06.012037Z",
     "iopub.status.idle": "2023-07-26T05:18:06.033832Z",
     "shell.execute_reply": "2023-07-26T05:18:06.033320Z"
    },
    "lines_to_next_cell": 2
   },
   "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>t_0</th>\n",
       "      <td>-1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>null_distribution</th>\n",
       "      <td>chi squared</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>degrees_of_freedom</th>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>test_name</th>\n",
       "      <td>multivariate_logrank_test</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>test_statistic</th>\n",
       "      <th>p</th>\n",
       "      <th>-log2(p)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>20.30</td>\n",
       "      <td>&lt;0.005</td>\n",
       "      <td>14.65</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>"
      ],
      "text/latex": [
       "\\begin{tabular}{lrrr}\n",
       " & test_statistic & p & -log2(p) \\\\\n",
       "0 & 20.30 & 0.00 & 14.65 \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "<lifelines.StatisticalResult: multivariate_logrank_test>\n",
       "               t_0 = -1\n",
       " null_distribution = chi squared\n",
       "degrees_of_freedom = 2\n",
       "         test_name = multivariate_logrank_test\n",
       "\n",
       "---\n",
       " test_statistic      p  -log2(p)\n",
       "          20.30 <0.005     14.65"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "multivariate_logrank_test(D['Wait time'],\n",
    "                          D['Center'],\n",
    "                          D['Failed'])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be5055e4",
   "metadata": {},
   "source": [
    "Next, we consider the  effect of `Time`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "9badb3e3",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:06.036375Z",
     "iopub.status.busy": "2023-07-26T05:18:06.036170Z",
     "iopub.status.idle": "2023-07-26T05:18:06.055850Z",
     "shell.execute_reply": "2023-07-26T05:18:06.055183Z"
    },
    "lines_to_next_cell": 2
   },
   "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>t_0</th>\n",
       "      <td>-1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>null_distribution</th>\n",
       "      <td>chi squared</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>degrees_of_freedom</th>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>test_name</th>\n",
       "      <td>multivariate_logrank_test</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>test_statistic</th>\n",
       "      <th>p</th>\n",
       "      <th>-log2(p)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>49.90</td>\n",
       "      <td>&lt;0.005</td>\n",
       "      <td>35.99</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>"
      ],
      "text/latex": [
       "\\begin{tabular}{lrrr}\n",
       " & test_statistic & p & -log2(p) \\\\\n",
       "0 & 49.90 & 0.00 & 35.99 \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "<lifelines.StatisticalResult: multivariate_logrank_test>\n",
       "               t_0 = -1\n",
       " null_distribution = chi squared\n",
       "degrees_of_freedom = 2\n",
       "         test_name = multivariate_logrank_test\n",
       "\n",
       "---\n",
       " test_statistic      p  -log2(p)\n",
       "          49.90 <0.005     35.99"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "multivariate_logrank_test(D['Wait time'],\n",
    "                          D['Time'],\n",
    "                          D['Failed'])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64b2bc33",
   "metadata": {},
   "source": [
    "As in the case of a categorical variable with 2 levels, these\n",
    "results are similar to the likelihood ratio test\n",
    "from the Cox proportional hazards model. First, we\n",
    "look at the results for  `Center`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "026e9ff8",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:06.058057Z",
     "iopub.status.busy": "2023-07-26T05:18:06.057906Z",
     "iopub.status.idle": "2023-07-26T05:18:06.184883Z",
     "shell.execute_reply": "2023-07-26T05:18:06.184356Z"
    },
    "lines_to_next_cell": 2
   },
   "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>null_distribution</th>\n",
       "      <td>chi squared</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>degrees_freedom</th>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>test_name</th>\n",
       "      <td>log-likelihood ratio test</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>test_statistic</th>\n",
       "      <th>p</th>\n",
       "      <th>-log2(p)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>20.58</td>\n",
       "      <td>&lt;0.005</td>\n",
       "      <td>14.85</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>"
      ],
      "text/latex": [
       "\\begin{tabular}{lrrr}\n",
       " & test_statistic & p & -log2(p) \\\\\n",
       "0 & 20.58 & 0.00 & 14.85 \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "<lifelines.StatisticalResult: log-likelihood ratio test>\n",
       "null_distribution = chi squared\n",
       "  degrees_freedom = 2\n",
       "        test_name = log-likelihood ratio test\n",
       "\n",
       "---\n",
       " test_statistic      p  -log2(p)\n",
       "          20.58 <0.005     14.85"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = MS(['Wait time',\n",
    "        'Failed',\n",
    "        'Center'],\n",
    "        intercept=False).fit_transform(D)\n",
    "F = coxph().fit(X, 'Wait time', 'Failed')\n",
    "F.log_likelihood_ratio_test()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ed54fe0",
   "metadata": {},
   "source": [
    "Next, we look at the results for `Time`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "7cab3789",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:06.187428Z",
     "iopub.status.busy": "2023-07-26T05:18:06.187255Z",
     "iopub.status.idle": "2023-07-26T05:18:06.315832Z",
     "shell.execute_reply": "2023-07-26T05:18:06.314374Z"
    },
    "lines_to_next_cell": 2
   },
   "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>null_distribution</th>\n",
       "      <td>chi squared</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>degrees_freedom</th>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>test_name</th>\n",
       "      <td>log-likelihood ratio test</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>test_statistic</th>\n",
       "      <th>p</th>\n",
       "      <th>-log2(p)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>48.12</td>\n",
       "      <td>&lt;0.005</td>\n",
       "      <td>34.71</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>"
      ],
      "text/latex": [
       "\\begin{tabular}{lrrr}\n",
       " & test_statistic & p & -log2(p) \\\\\n",
       "0 & 48.12 & 0.00 & 34.71 \\\\\n",
       "\\end{tabular}\n"
      ],
      "text/plain": [
       "<lifelines.StatisticalResult: log-likelihood ratio test>\n",
       "null_distribution = chi squared\n",
       "  degrees_freedom = 2\n",
       "        test_name = log-likelihood ratio test\n",
       "\n",
       "---\n",
       " test_statistic      p  -log2(p)\n",
       "          48.12 <0.005     34.71"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = MS(['Wait time',\n",
    "        'Failed',\n",
    "        'Time'],\n",
    "       intercept=False).fit_transform(D)\n",
    "F = coxph().fit(X, 'Wait time', 'Failed')\n",
    "F.log_likelihood_ratio_test()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d250dc9",
   "metadata": {},
   "source": [
    "We find that differences between centers are highly significant, as\n",
    "are differences between times of day.\n",
    "\n",
    "Finally, we fit Cox's proportional hazards model to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "5cc4b898",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-07-26T05:18:06.320288Z",
     "iopub.status.busy": "2023-07-26T05:18:06.319261Z",
     "iopub.status.idle": "2023-07-26T05:18:06.585074Z",
     "shell.execute_reply": "2023-07-26T05:18:06.577564Z"
    },
    "lines_to_next_cell": 2
   },
   "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>coef</th>\n",
       "      <th>se(coef)</th>\n",
       "      <th>p</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>covariate</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>Operators</th>\n",
       "      <td>0.043934</td>\n",
       "      <td>0.007520</td>\n",
       "      <td>5.143589e-09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Center[B]</th>\n",
       "      <td>-0.236060</td>\n",
       "      <td>0.058113</td>\n",
       "      <td>4.864162e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Center[C]</th>\n",
       "      <td>0.012231</td>\n",
       "      <td>0.057518</td>\n",
       "      <td>8.316096e-01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Time[Even.]</th>\n",
       "      <td>0.268845</td>\n",
       "      <td>0.057797</td>\n",
       "      <td>3.294956e-06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Time[Morn.]</th>\n",
       "      <td>-0.148217</td>\n",
       "      <td>0.057334</td>\n",
       "      <td>9.733557e-03</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 coef  se(coef)             p\n",
       "covariate                                    \n",
       "Operators    0.043934  0.007520  5.143589e-09\n",
       "Center[B]   -0.236060  0.058113  4.864162e-05\n",
       "Center[C]    0.012231  0.057518  8.316096e-01\n",
       "Time[Even.]  0.268845  0.057797  3.294956e-06\n",
       "Time[Morn.] -0.148217  0.057334  9.733557e-03"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = MS(D.columns,\n",
    "       intercept=False).fit_transform(D)\n",
    "fit_queuing = coxph().fit(\n",
    "                  X,\n",
    "                 'Wait time',\n",
    "                 'Failed')\n",
    "fit_queuing.summary[['coef', 'se(coef)', 'p']]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bec9d61d",
   "metadata": {},
   "source": [
    "The $p$-values for Center B and evening time\n",
    "are very small. It is also clear that the\n",
    "hazard --- that is, the instantaneous risk that a call will be\n",
    "answered --- increases with the number of operators. Since we\n",
    "generated the data ourselves, we know that the true coefficients for\n",
    " `Operators`, `Center = B`, `Center = C`, \n",
    "`Time = Even.` and `Time = Morn.`   are $0.04$, $-0.3$,\n",
    "$0$,   $0.2$, and $-0.2$, respectively. The coefficient estimates\n",
    "from the fitted Cox model are fairly accurate.\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "main_language": "python",
   "notebook_metadata_filter": "-all"
  },
  "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.9.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
