{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22660227",
   "metadata": {
    "lines_to_next_cell": 0,
    "name": "setup",
    "tags": [
     "remove_input"
    ]
   },
   "outputs": [],
   "source": [
    "knitr::opts_chunk$set(error = TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40d5521d",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b3e8d85",
   "metadata": {},
   "source": [
    "# Survival Analysis\n",
    "\n",
    "In this lab, we perform survival analyses on three separate data sets. In Section 11.8.1 we analyze the `BrainCancer` data that was  first described in Section 11.3.\n",
    " In Section 11.8.2, we examine the `Publication` data from Section 11.5.4. Finally, Section 11.8.3 explores a simulated call center data set.\n",
    "\n",
    "## Brain Cancer Data\n",
    "We begin with the `BrainCancer` data set, which is part of the  `ISLR2` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "497a9e23",
   "metadata": {
    "name": "chunk1"
   },
   "outputs": [],
   "source": [
    "library(ISLR2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d96d5011",
   "metadata": {},
   "source": [
    "The rows index the 88 patients, while the columns contain the 8 predictors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44f4dee0",
   "metadata": {
    "name": "chunk2"
   },
   "outputs": [],
   "source": [
    "names(BrainCancer)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22b1813c",
   "metadata": {},
   "source": [
    "We first briefly examine the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3cc1eff",
   "metadata": {
    "name": "chunk3"
   },
   "outputs": [],
   "source": [
    "attach(BrainCancer)\n",
    "table(sex)\n",
    "table(diagnosis)\n",
    "table(status)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef3774d6",
   "metadata": {},
   "source": [
    "Before beginning an analysis, it is important to know how the `status` variable has been coded.  Most software, including `R`, uses the convention that `status = 1` indicates an uncensored observation,  and `status = 0` indicates a censored observation. But some scientists might use the opposite coding. For the `BrainCancer` data set $35$ patients died before the end of the study.\n",
    "\n",
    "To begin the analysis, we re-create  the Kaplan-Meier survival curve  shown in  Figure 11.2, using the `survfit()` function within the `R` `survival` library. Here `time` corresponds to $y_i$, the time to the $i$th event (either censoring or death)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f40ac817",
   "metadata": {
    "name": "chunk4"
   },
   "outputs": [],
   "source": [
    "library(survival)\n",
    "fit.surv <- survfit(Surv(time, status) ~ 1)\n",
    "plot(fit.surv, xlab = \"Months\",\n",
    "    ylab = \"Estimated Probability of Survival\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92184a20",
   "metadata": {},
   "source": [
    "Next we create Kaplan-Meier survival curves that are stratified by `sex`, in order to reproduce Figure 11.3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9da8f52a",
   "metadata": {
    "lines_to_next_cell": 2,
    "name": "chunk5"
   },
   "outputs": [],
   "source": [
    "fit.sex <- survfit(Surv(time, status) ~ sex)\n",
    "plot(fit.sex, xlab = \"Months\",\n",
    "    ylab = \"Estimated Probability of Survival\", col = c(2,4))\n",
    "legend(\"bottomleft\", levels(sex), col = c(2,4), lty = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00cf7223",
   "metadata": {},
   "source": [
    "As discussed in Section 11.4, we can perform a log-rank test to compare the survival of males to females, using the `survdiff()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6925849d",
   "metadata": {
    "name": "chunk6"
   },
   "outputs": [],
   "source": [
    "logrank.test <- survdiff(Surv(time, status) ~ sex)\n",
    "logrank.test"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6644ec43",
   "metadata": {},
   "source": [
    "The resulting $p$-value is $0.23$, indicating  no evidence of a difference in survival between the two sexes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d718d54",
   "metadata": {},
   "source": [
    "Next, we fit  Cox proportional hazards models using the `coxph()`  function.\n",
    "To begin, we consider a model that uses `sex` as the only predictor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b529a4b6",
   "metadata": {
    "name": "chunk7"
   },
   "outputs": [],
   "source": [
    "fit.cox <- coxph(Surv(time, status) ~ sex)\n",
    "summary(fit.cox)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb4f0b8a",
   "metadata": {},
   "source": [
    "Note that the values of the likelihood ratio, Wald, and score tests have been rounded.  It is possible to display additional digits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b610f27c",
   "metadata": {
    "name": "chunk8"
   },
   "outputs": [],
   "source": [
    "summary(fit.cox)$logtest[1]\n",
    "summary(fit.cox)$waldtest[1]\n",
    "summary(fit.cox)$sctest[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cc741d9",
   "metadata": {},
   "source": [
    "Regardless of which test we use, we see that there is no clear evidence for a difference in survival between males and females."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b67113bd",
   "metadata": {
    "name": "chunk9"
   },
   "outputs": [],
   "source": [
    "logrank.test$chisq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b01fd57",
   "metadata": {},
   "source": [
    "As we learned in this chapter, the score test from the Cox model is exactly equal to the log rank test statistic!\n",
    "\n",
    "Now we fit a  model that makes use of additional predictors. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3ef29b6",
   "metadata": {
    "name": "chunk10"
   },
   "outputs": [],
   "source": [
    "fit.all <- coxph(\n",
    "Surv(time, status) ~ sex + diagnosis + loc + ki + gtv +\n",
    "   stereo)\n",
    "fit.all"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4279f4ca",
   "metadata": {},
   "source": [
    "The `diagnosis` variable has been coded so that the baseline\n",
    "corresponds to meningioma. The results indicate that the risk associated with HG glioma is more than eight times (i.e. $e^{2.15}=8.62$) the risk associated with meningioma. In other words, after adjusting for the other predictors, patients with HG glioma have much worse survival\n",
    "compared to those with meningioma.  In addition, larger values of the Karnofsky index, `ki`, are associated with lower risk, i.e. longer survival.\n",
    "\n",
    "Finally, we plot survival curves for each diagnosis category,  adjusting for the other predictors.\n",
    "To make these plots, we set the values of the other predictors equal\n",
    "to the mean for quantitative variables, and the modal value for\n",
    "factors. We first create a data frame with four rows, one for each\n",
    "level of diagnosis. The `survfit()` function will produce a curve for each of the rows in this data frame,\n",
    "and one call to `plot()` will display them all in the same plot.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba945922",
   "metadata": {
    "lines_to_next_cell": 0,
    "name": "chunk11"
   },
   "outputs": [],
   "source": [
    "modaldata <- data.frame(\n",
    "     diagnosis = levels(diagnosis),\n",
    "     sex = rep(\"Female\", 4),\n",
    "     loc = rep(\"Supratentorial\", 4),\n",
    "     ki = rep(mean(ki), 4),\n",
    "     gtv = rep(mean(gtv), 4),\n",
    "     stereo = rep(\"SRT\", 4)\n",
    "     )\n",
    "survplots <- survfit(fit.all, newdata = modaldata)\n",
    "plot(survplots, xlab = \"Months\",\n",
    "    ylab = \"Survival Probability\", col = 2:5)\n",
    "legend(\"bottomleft\", levels(diagnosis), col = 2:5, lty = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87dcfcdf",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f20de24",
   "metadata": {},
   "source": [
    "## Publication Data\n",
    "\n",
    "The `Publication` data presented in Section 11.5.4 can be found in the `ISLR2` library. We first reproduce Figure 11.5 by plotting the Kaplan-Meier curves stratified on the `posres` variable, which records whether the study had a positive or negative result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e9570f4",
   "metadata": {
    "lines_to_next_cell": 2,
    "name": "chunk12"
   },
   "outputs": [],
   "source": [
    "fit.posres <- survfit(\n",
    "    Surv(time, status) ~ posres, data = Publication\n",
    "  )\n",
    "plot(fit.posres, xlab = \"Months\",\n",
    "    ylab = \"Probability of Not Being Published\", col = 3:4)\n",
    "legend(\"topright\", c(\"Negative Result\", \"Positive Result\"),\n",
    "    col = 3:4, lty = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5856c8f3",
   "metadata": {},
   "source": [
    "As discussed previously, the $p$-values from fitting Cox's proportional hazards model to the `posres` variable are quite large, providing no evidence of a difference in time-to-publication between studies with positive versus negative results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0533651d",
   "metadata": {
    "lines_to_next_cell": 2,
    "name": "chunk13"
   },
   "outputs": [],
   "source": [
    "fit.pub <- coxph(Surv(time, status) ~ posres,\n",
    "    data = Publication)\n",
    "fit.pub"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ef39682",
   "metadata": {},
   "source": [
    "As expected, the log-rank test provides an identical conclusion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5f741e5",
   "metadata": {
    "lines_to_next_cell": 2,
    "name": "chunk14"
   },
   "outputs": [],
   "source": [
    "logrank.test <- survdiff(Surv(time, status) ~ posres,\n",
    "    data = Publication)\n",
    "logrank.test"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f429ef11",
   "metadata": {},
   "source": [
    "However, the results change dramatically when we include other  predictors in the model. Here we have excluded the funding mechanism variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40e25769",
   "metadata": {
    "name": "chunk15"
   },
   "outputs": [],
   "source": [
    "fit.pub2 <- coxph(Surv(time, status) ~ . - mech,\n",
    "    data = Publication)\n",
    "fit.pub2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99c006b7",
   "metadata": {},
   "source": [
    "We see that there are a number of  statistically significant variables, including whether the trial focused on a clinical endpoint, the impact of the study, and whether the study had positive or negative results.\n",
    "\n",
    "## Call Center Data\n",
    "\n",
    "In this section, we will simulate survival data using the `sim.survdata()` function, which is part of the `coxed` library. Our simulated data will represent the observed wait times (in seconds) for 2,000 customers who have phoned  a call center.\n",
    "In this context, censoring occurs if a customer hangs up before his or her call is answered.\n",
    "\n",
    "There are three covariates: `Operators` (the number of call center operators available at the time of the call, which can range from $5$ to $15$), `Center` (either A, B, or C), and `Time` of day (Morning, Afternoon, or Evening). We generate data for these covariates so that all possibilities are equally likely: for instance, morning, afternoon and evening calls are equally likely, and any number of operators from $5$ to $15$ is equally likely."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b33c7c91",
   "metadata": {
    "name": "chunk16"
   },
   "outputs": [],
   "source": [
    "set.seed(4)\n",
    "N <- 2000\n",
    "Operators <- sample(5:15, N, replace = T)\n",
    "Center <- sample(c(\"A\", \"B\", \"C\"), N, replace = T)\n",
    "Time <- sample(c(\"Morn.\", \"After.\", \"Even.\"), N, replace = T)\n",
    "X <- model.matrix( ~ Operators + Center + Time)[, -1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d172deaa",
   "metadata": {},
   "source": [
    "It is worthwhile to take a peek at the design matrix `X`,  so that we can be sure that we understand how the\n",
    "variables have been coded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e51c552",
   "metadata": {
    "name": "chunk17"
   },
   "outputs": [],
   "source": [
    "X[1:5, ]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "051bcd70",
   "metadata": {},
   "source": [
    "Next,  we specify the coefficients and the hazard function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6ca61d2",
   "metadata": {
    "name": "chunk18"
   },
   "outputs": [],
   "source": [
    "true.beta <- c(0.04, -0.3, 0, 0.2, -0.2)\n",
    "h.fn <- function(x) return(0.00001 * x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d827ceb",
   "metadata": {},
   "source": [
    "Here, we have set the coefficient associated with `Operators` to equal $0.04$; in other words, each additional operator leads to a $e^{0.04}=1.041$-fold increase in the\n",
    "  \"risk\" that the call will be answered,\n",
    "given the `Center` and `Time` covariates. This makes sense: the greater the number of operators at hand, the shorter the wait time! The coefficient associated with `Center = B`\n",
    "is $-0.3$, and `Center = A` is treated as the baseline. This means that the risk of a call being\n",
    "answered at Center B is $0.74$ times the risk that it will be answered at Center A; in other words,\n",
    "the wait times are a bit longer at Center B.\n",
    "\n",
    "We are now ready to  generate data under the  Cox proportional hazards model. The `sim.survdata()` function allows us to specify the maximum possible failure time, which in this case corresponds to the longest possible wait time for a customer; we set this to equal $1{,}000$ seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d60561ad",
   "metadata": {
    "lines_to_next_cell": 0,
    "name": "chunk19"
   },
   "outputs": [],
   "source": [
    "library(coxed)\n",
    "queuing <- sim.survdata(N = N, T = 1000, X = X,\n",
    "    beta = true.beta, hazard.fun = h.fn)\n",
    "names(queuing)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a77ebed",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "24dbbfe8",
   "metadata": {},
   "source": [
    "The \"observed\" data is stored in `queuing$data`, with `y` corresponding to the event time and `failed` an indicator of whether the call was answered (`failed = T`) or the customer hung up before the call was answered (`failed = F`). We see that almost $90\\%$ of calls were answered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "511496b3",
   "metadata": {
    "lines_to_next_cell": 0,
    "name": "chunk20"
   },
   "outputs": [],
   "source": [
    "head(queuing$data)\n",
    "mean(queuing$data$failed)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "061f89d9",
   "metadata": {},
   "source": [
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29e3e744",
   "metadata": {},
   "source": [
    "We now plot  Kaplan-Meier survival curves. First, we stratify by `Center`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "807ed968",
   "metadata": {
    "name": "chunk21"
   },
   "outputs": [],
   "source": [
    "par(mfrow = c(1, 2))\n",
    "fit.Center <- survfit(Surv(y, failed) ~ Center,\n",
    "    data = queuing$data)\n",
    "plot(fit.Center, xlab = \"Seconds\",\n",
    "    ylab = \"Probability of Still Being on Hold\",\n",
    "    col = c(2, 4, 5))\n",
    "legend(\"topright\",\n",
    "     c(\"Call Center A\", \"Call Center B\", \"Call Center C\"),\n",
    "     col = c(2, 4, 5), lty = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "daff5cc6",
   "metadata": {},
   "source": [
    "Next, we stratify by `Time`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec47f2fd",
   "metadata": {
    "name": "chunk22"
   },
   "outputs": [],
   "source": [
    "fit.Time <- survfit(Surv(y, failed) ~ Time,\n",
    "   data = queuing$data)\n",
    "plot(fit.Time, xlab = \"Seconds\",\n",
    "    ylab = \"Probability of Still Being on Hold\",\n",
    "    col = c(2, 4, 5))\n",
    "legend(\"topright\", c(\"Morning\", \"Afternoon\", \"Evening\"),\n",
    "    col = c(5, 2, 4), lty = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f4a671f",
   "metadata": {},
   "source": [
    "It seems that calls at Call Center B\n",
    "  take longer to be answered than calls at Centers A and C. Similarly, it appears that wait times are longest in the morning and shortest in the evening hours. We can use a log-rank test to determine whether these differences are statistically significant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "266ccf00",
   "metadata": {
    "lines_to_next_cell": 0,
    "name": "chunk23"
   },
   "outputs": [],
   "source": [
    "survdiff(Surv(y, failed) ~ Center, data = queuing$data)\n",
    "survdiff(Surv(y, failed) ~ Time, data = queuing$data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07f0be29",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "db72ddaa",
   "metadata": {},
   "source": [
    "We find that  differences between centers are highly significant, as are differences between times of day."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52790932",
   "metadata": {},
   "source": [
    "Finally, we fit Cox's proportional hazards model to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9c88320",
   "metadata": {
    "name": "chunk24"
   },
   "outputs": [],
   "source": [
    "fit.queuing <- coxph(Surv(y, failed) ~ .,\n",
    "    data = queuing$data)\n",
    "fit.queuing"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a410178e",
   "metadata": {},
   "source": [
    "The $p$-values for `Center = B`, `Time = Even.` and `Time = Morn.` are very small. It is also clear that the hazard --- that is, the instantaneous risk that a call will be answered --- increases with the number of operators. Since we generated the data ourselves, we know that the true coefficients for `Operators`, `Center = B`, `Center = C`, `Time = Even.` and `Time = Morn.` are  $0.04$, $-0.3$, $0$, $0.2$, and $-0.2$, respectively. The coefficient estimates resulting from the Cox model are fairly accurate.\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "name,tags,-all"
  },
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}