---
title: "Statistical Thinking Tutorial"
author: "Week 6"
date: "29/08/2021"
output:
  pdf_document: default
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## What's on deck for this week

This week we are going to look at censored data. We will

- A quick check in with QQ-plots
- Learn how to fit survival data with the `survival` package 
- A closer look at the Kaplan-Meier curve 
- A better look at medians
- Fitting parametric survival models


# An amuse-bouche 


a. Simulate samples of size $n=$ 30, 100, 500 from these distributions
    i. Lognormal(4, 2)
    ii. Gamma(3, 3)
b. Make a QQ-plot of each these samples comparing to . Explain how closely the samples, of different sizes, appears to match the theoretical distribution.
c. By drawing 100 samples of size 30 from the Gamma(3,3), visualise how variable the QQ-plot can be.

## Dealing with right-censored data in R

Right censored data is a common type of "time to event" data where the event
is only observed for some of the observations. An example would be when the 
event is a death due to acute myelogenous leukemia. Some patients in the sample
will die before the end of the study, while some (hopefully most!) will still
have not died when the study ends.

In class, we talked about how important^[We talked about the statistical reasons
for doing this, but when dealing with serious data about death, we are bound both
morally and ethically to use our data as best we can.] it is to correctly account 
for both the events that we have observed and the ones that will occur, 
unobserved, in the future.

This data is so common that there is a packing in R that deals with it 
_that you you don't even have to install_. It's part of the basic R distribution.
This package is, creatively, names `survival`^[Survival can deal with more 
complex censoring than just right censoring. It does a bunch of stuff that 
we are not going to look at in this class. It's a real industrial-strength thing.].

In order to deal with survival data, we need to keep track of two things:

1. The observed time
2. An indicator of whether this was an event or if we just stopped measuring (a 
censoring). You can check the help for other options, but this is usally coded as 
a `1` if it's an event and a `0` if it is not.

We can see this data structure in the `aml` dataset (the aforementioned leukemia) 
from the `survival` package.

```{r, echo = TRUE, eval = TRUE, message = FALSE}
library(survival)
library(tidyverse)
head(aml)
dim(aml)
```

This data contains an extra factor (`x`) that indicates what type of chemotherapy
regime was used. We are going to use that to filter the data into subpopulations.

One of the key functions in the survival package is `Surv`, which lets us
put the time and event information together.

```{r, echo = TRUE, eval = TRUE}
Surv(aml$time, aml$status) %>% head()
## If you don't understand it, you can view on the next code
view(cbind(aml, Surv(aml$time, aml$status)))
```

Notice that this separates the times into numbers (9, 13) and cenored times (13+), 
where the + indicates that the event occurs at some point after 13.

## Making a survival curve

To fit survival data, we use the `survival::survfit()` function.

```{r, echo = TRUE, eval = TRUE}
fit <- survfit(Surv(time, status) ~ 1, data = aml)
```

Let's break that down:

- The `Surv(time, status)` term is there to link the time to event with the event
indicator
- The tilde (`~`) is a special character in R that we will see more of when we start
moving into regression, but it exists to separate the data (the left of `~`) from 
the "predictors" on the right. For today, we are just using `~ 1`, which says "there
are no predictors".
- The `data = aml` tells the function which data to fit the survival curve on (and
where to look for the columns `time` and `status`)

The fit object contains a bunch of stuff, but the easiest thing we can do is look
at what it prints out as a default
```{r, echo = TRUE, eval = TRUE}
fit
```

This gives the number of observations, the number of observed events and an estimate
with 95% confidence interval for the median survival time. We could've also gotten
this median with the commnad 

```{r, echo = TRUE, eval = TRUE}
median(Surv(aml$time, aml$status))
```

(The `50` is labelling the median as the 50% quantile. Is it terrible output? Yes.)

We can use the `fit` object to plot the Kaplan-Meier estimate of the survival 
curve. There are a few ways to do this, but we will use the `survminer::ggsurvplot()`
function. To use this, you will need to install the `survminer` package.

```{r, echo = TRUE, eval = FALSE}
install.packages("survminer")
```

We can then make the plot.

```{r echo = TRUE, eval = TRUE}
library(survminer)
ggsurvplot(fit)
```
By default, this plots a confidence interval and you can turn this off with the 
option `conf.int = FALSE`.


### Question
1. You can plot the Kaplan-Meier curves for the sub-populations defined by the 
covariate `x`, by replacing `~ 1` with `~ x` in the `survfit` call. Look at the 
resulting medians and Kaplan-Meier curves (use the same call to `ggsurvplot`). 
Write a few sentences commenting on the similarities and differences.

## A closer look at the Kaplan-Meier curve

So how did `ggsurvplot` make that plot? If we look a little closer at the `fit` 
object, we can work it out .

```{r echo = TRUE, eval  =TRUE}
glimpse(fit)

```

These are all the things that are computed when you call `survfit`.

- `time` is the set of unique times in the data ($t_j$ in the lectures)
- `n.event` is the number of events that were observed at each time ($d_j$ in 
the lectures)
- `n.censor` is the number of censored observations at each time
- `n.risk` is the number at risk at each time ($n_j$ in the lectures)
- `surv` is the Kaplan-Meier estimate of the survival time.

We can check that we can compute `surv` from `n.event` and `n.risk`:

```{r, echo = TRUE, eval = TRUE}
cumprod(1 - fit$n.event / fit$n.risk) - fit$surv
```

### Questions
- What does it mean that the above expression is (almost) zero?
- Why isn't it _exactly_ zero?
- What does `cumprod` do? (Also look up what `cumsum` does)
- How can you compute `n.censor` from this data?

## Computing the median

Finally, we can use the estimated survival curve to find the median. This is not
unique except for in the very unusual situation where  `fit$surv` is exactly `0.5`
at one time.  When there is no unique median, we have a number of options, two
popular ones are:

- Select the average of the two closest times with `fit$surv` on either side of 0.5
- Select the time with `fit$surv` closest to 0.5.

For big enough data sets, both of these should give about the same value, but 
for small data sets the second option (which is what R uses) is slightly more reliable
as it's possible that one of the two closest times has a survival estimate that is quite
far from 0.5.

Let's look at this!
```{r, echo = TRUE, eval = TRUE}
summary(fit)
```
The "closest" median is `27`, while the "average" median is `r (23+27)/2`.

### Questions

- Write functions to find the median using both methods. They should take a `fit` 
as an input.

- Consider the `aml` data when `x == "Maintained"`. Compute the Kaplan-Meier curve
for this subset of the data. Which of the two medians is most appropriate? 

- Repeat the exercise for `x == Nonmaintained`.

- Look at the `diabetic` data (load it with `data(diabetic)`). Why is the median 
of this data not defined?

## Fitting a parametric survival distribution

So far, we have only talked about non-parametric estimates of the survival curve.
What if we think that the uncensored data came from a particular distribution? 
Not every distribution can be used to fit survival data in R, but three of the 
most useful options are

- Exponential
- Weibull (_WARNING_ the parameterisation of the Weibull that `survival` uses is different to the standard one in R. Don't even ask.)
- Log-normal.

We fit these models using the `survival::survreg` function

```{r, echo = TRUE, eval = TRUE}
fit2 <- survreg(Surv(time, status) ~ 1, data = aml, dist = "exponential")
fit2
```
What do the paraemters mean? The `(Intercept)` term is the _log_ mean of the survival 
distribution (in this case the exponential). So the single unknown parameter in the
exponential distribution can be computed as

```{r eval = TRUE, echo = TRUE}

lambda <- exp(-coef(fit2))

```

We can then look up what the median of an exponentially distributed random variable is.

![]("exponential.png")

Using this we can compare the parametric estimate of the median to the Kaplan-Meier
estimate.
```{r, eval = TRUE, echo = TRUE}
log(2) / lambda
median(Surv(aml$time, aml$status))$quantile
```

We can also compare the survival curves.

```{r, echo = TRUE, eval = TRUE}
km_curve <- ggsurvplot(fit, conf.int = FALSE)
km_curve$plot + stat_function(fun = pexp, 
                              colour = "blue", 
                              line_type = "dashed",
                              args = list(rate = lambda,
                                          lower.tail = FALSE)
                              )
```

The resulting survival curve is a prettty good match with the empirical estimate
for early times, but after the median, the two estimates diverge wildly.

### Question

- Fit the data using a Weibull distribution. And compare the median and the survival curve and comment on your findings.

The following code will be useful for translating between the `survival` package's
parameterisation of the Weibull distribution and the standard R parameterisation.

```{r, echo = TRUE, eval = FALSE}
## fit3 is the weibull fit

pweibull_shape <- 1 / fit3$scale
pweibull_scale <- exp(coef(fit3))
```