Survival data (2)

class: center, middle, inverse, title-slide

.title[
# Survival data (2)
]
.author[
### Yue Jiang
]
.date[
### STA 440 / Spring 25
]

---

### Estimation

Suppose we had a series of *non-censored* observations `$t_1, t_2, \cdots, t_n$`, 
and we thought that these survival times were i.i.d. and came from some
distribution with density `$f(t)$`.

.question[
How might we estimate the parameter(s) `$\theta$` of that distribution?
]

---

### Review: maximum likelihood estimation

`\begin{align*}
\mathcal{L}(\theta | T) &= f(t_1, t_2, \cdots, t_n | \theta)\\
&= f(t_1 | \theta)f(t_2 | \theta) \cdots f(t_n | \theta)\\
&= \prod_{i = 1}^n f(t_i|\theta),
\end{align*}`

which we can often maximize in closed form for many familiar distributions (or 
numerically).

---

### Review: maximum likelihood estimation

For instance, if we thought `$T \stackrel{i.i.d}{\sim} Exp(\lambda)$`, then

`\begin{align*}
\mathcal{L}(\lambda | T) &= \lambda^n \exp \left(-\lambda\sum_{i = 1}^n t_n \right)\\
\log \mathcal{L}(\lambda | T) &= n\log(\lambda) - \lambda \sum_{i = 1}^n t_i,\\
\hat{\lambda}_{MLE} &= \frac{n}{\sum_{i = 1}^n t_i}
\end{align*}`

---

### Estimation for censored data

How might we perform maximum likelihood estimation for *censored* data? What
would the likelihood look like? Suppose we have `$n$` i.i.d. observations with the
same `$f(t)$`, `$S(t)$`, and hazard `$\lambda(t)$`. Consider what might happen at time
`$t_i$`:

- Case 1: Individual experiences event at `$t_i: \delta_i = 1$`. In 
this case, they contribute the density at `$t_i$`, which is `$f(t_i)$`.

.question[
What would their contribution to the likelihood be if they were censored?
]

- Case 2: Individual is still alive at `$t_i: \delta_i = 0$`. In this
case, we know that their survival time is greater than `$t_i$`, which happens with
probability `$S(t_i)$`.

So, this individual's contribution to the likelihood is

`\begin{align*}
f(t_i)^{\delta_i}S(t_i)^{1 - \delta_i}
\end{align*}`

---

### Estimation for censored data

Note that

`\begin{align*}
\lambda(t) &= \lim_{\Delta t \to 0^+} \frac{P(t \le T < t + \Delta t | T \ge t)}{\Delta t}\\
&= \lim_{\Delta t \to 0^+} \frac{P(t \le T < t + \Delta t \cap T \ge t)/P(T \ge t)}{\Delta t}\\
&= \lim_{\Delta t \to 0^+} \frac{P(t \le T < t + \Delta t)}{P(T \ge t)\Delta t}\\
&= \lim_{\Delta t \to 0^+} \frac{f(t)\Delta t}{P(T \ge t)\Delta t}\\
&= \frac{f(t)}{S(t)}, \\
f(t) &= \lambda(t)S(t)
\end{align*}`

---

### Estimation for censored data

An individual's contribution to the likelihood is thus

`\begin{align*}
&\mathrel{\phantom{=}} f(t_i)^{\delta_i}S(t_i)^{1 - \delta_i}\\
&= \lambda(t_i)^{\delta_i}S(t_i)^{\delta_i}S(t_i)^{1 - \delta_i}\\
&= \lambda(t_i)^{\delta_i}S(t_i)
\end{align*}`

For our exponential example with `$T \stackrel{i.i.d.}{\sim} Exp(\lambda)$`, since
the hazard is given by `$\lambda$` and the survival function is given by 
`$e^{-\lambda t}$`, we would thus have

`\begin{align*}
\mathcal{L}(\lambda | T) &= \prod_{i = 1}^n \lambda^{\delta_i}\exp(-\lambda t_i),
\end{align*}`

which we could then maximize using familiar methods.

---

### Covariates?

.question[
How might we adjust for potential covariates that might be related to our
outcome? Is there any way to create a predictive model for survival time?
]

---

### Some candidate models

Could we use the following model for survival outcomes?

`\begin{align*}
T_i = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \cdots + \beta_px_{pi} + \epsilon_i
\end{align*}`

How about this model?

`\begin{align*}
\log(T_i) = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \cdots + \beta_px_{pi} + \epsilon_i
\end{align*}`

---

### An accelerated failure time model

An .vocab[accelerated failure time] (AFT) model assumes

`\begin{align*}
\log(T_i) = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \cdots + \beta_px_{pi} + \epsilon_i
\end{align*}`

where `$\epsilon_i$` are commonly assumed to be i.i.d. and follow some specified
distribution.

.question[
How might we "know" what distribution `$\epsilon$` takes on?
]

---

### An accelerated failure time model

`\begin{align*}
\log(T_i) = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \cdots + \beta_px_{pi} + \epsilon_i
\end{align*}`

There is a one-to-one relationship between the distribution of `$T$` and the
assumed error distribution in the AFT model.

For instance, assuming `$\epsilon_i$` has a normal distribution implies `$T_i$` 
has a log-normal distribution.

---

### An accelerated failure time model

Note that we can also write the AFT model as

`\begin{align*}
\log(T_i) &= \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \cdots + \beta_px_{pi} + \epsilon_i\\
T_i &= \exp\left(\beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \cdots + \beta_px_{pi} \right)e^{\epsilon_i}
\end{align*}`

.question[
How might we interpret coefficients in the AFT model?
]

---

### An accelerated failure time model

`\begin{align*}
T_i &= \exp\left(\beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \cdots + \beta_px_{pi} \right)e^{\epsilon_i}\\
&= e^{\beta_0}e^{\beta_1x_{1i}}e^{\beta_2x_{2i}}\cdots e^{\beta_px_{pi}}e^{\epsilon_i}
\end{align*}`

Covariates in an AFT model have a multiplicative effect on *time*.

For instance,
if `$\beta_k = 0.4$`, then `$\exp(\beta_k) \approx 1.5$`. Holding all else equal,
an individual with covariate `$x_k$` one unit greater than another is expected
to survive approximately 1.5 times longer than the other. The probability that
they have survived to time `$1.5t$` is the same as the probability that the other
individual is alive has survived to time `$t$`.

---

### Estimation in survival models

Remember that an individual's contribution to the likelihood is

`\begin{align*}
\lambda(t)^{\delta_i}S(t),
\end{align*}`

and so the full likelihood function is given by

`\begin{align*}
\mathcal{L}(\theta | Y) = \prod_{i = 1}^n \lambda(t_i)^{\delta_i}S(t_i).
\end{align*}`

This likelihood function holds regardless of what type of model is being
estimated (for independently right-censored data).

---

### Fitting an AFT model

Recall that our model specifies a functional form between `$Y$`, `$\beta$`, and 
`$\mathbf{X}$`. For the AFT model,

`\begin{align*}
\log(T_i) = \beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \cdots + \beta_px_{pi} + \epsilon_i,
\end{align*}`

Suppose we think `$T$` follows an exponential distribution such that
`$T \sim Exp(\exp(\mathbf{X}^T\beta))$`; that is, `$\epsilon$` has a standard 
extreme value distribution: `$f(u) = \exp(u - \exp(u))$`.

Then we have:

`\begin{align*}
S(t | \mathbf{X}) &= \exp(-t\exp(\mathbf{X}^T\beta))\\
\lambda(t | \mathbf{X}) &= \exp(\mathbf{X}^T\beta)
\end{align*}`

---

### Fitting an AFT model

Thus we can maximize the likelihood for this model

`\begin{align*}
&\mathrel{\phantom{=}} \prod_{i = 1}^n \lambda(t_i | \mathbf{X}_i)^{\delta_i}S(t_i | \mathbf{X}_i)\\
&= \prod_{i = 1}^n \exp(\mathbf{X}_i^T\beta)^{\delta_i}\exp(-t\exp(\mathbf{X}_i^T\beta))
\end{align*}`

with respect to the regression parameters of interest `$\beta$`. Different 
assumptions re: error distributions (and thus of `$T$`) result in different hazard 
and survival functions being plugged into our generic likelihood.

Note that there is no closed form for this likelihood function with respect to 
`$\beta$`; choose your favorite numerical method (we'll see Newton-Raphson once
again in model output on the following slides).

---

### Fitting an AFT model

```r
library(survival)
aft_e <- survreg(Surv(DWHFDAYS, DWHF) ~ TRTMT + EJFPER + PREVMI, 
                 data = dig, dist = "exponential")
summary(aft_e)
```

```
## 
## Call:
## survreg(formula = Surv(DWHFDAYS, DWHF) ~ TRTMT + EJFPER + PREVMI, 
##     data = dig, dist = "exponential")
##                Value Std. Error     z       p
## (Intercept)  6.62850    0.07495 88.43 < 2e-16
## TRTMT        0.29856    0.04167  7.17 7.7e-13
## EJFPER       0.04037    0.00239 16.92 < 2e-16
## PREVMI      -0.04978    0.04364 -1.14    0.25
## 
## Scale fixed at 1 
## 
## Exponential distribution
## Loglik(model)= -20511.5   Loglik(intercept only)= -20680.3
## 	Chisq= 337.51 on 3 degrees of freedom, p= 7.6e-73 
## Number of Newton-Raphson Iterations: 5 
## n= 6799
```

---

### Fitting an AFT model

```r
library(survival)
aft_w <- survreg(Surv(DWHFDAYS, DWHF) ~ TRTMT + EJFPER + PREVMI, 
                 data = dig, dist = "weibull")
summary(aft_w)
```

```
## 
## Call:
## survreg(formula = Surv(DWHFDAYS, DWHF) ~ TRTMT + EJFPER + PREVMI, 
##     data = dig, dist = "weibull")
##                Value Std. Error     z       p
## (Intercept)  6.58614    0.10266 64.15 < 2e-16
## TRTMT        0.39910    0.05768  6.92 4.5e-12
## EJFPER       0.05263    0.00337 15.61 < 2e-16
## PREVMI      -0.06661    0.06001 -1.11    0.27
## Log(scale)   0.31844    0.01891 16.84 < 2e-16
## 
## Scale= 1.37 
## 
## Weibull distribution
## Loglik(model)= -20351.2   Loglik(intercept only)= -20504.9
## 	Chisq= 307.43 on 3 degrees of freedom, p= 2.5e-66 
## Number of Newton-Raphson Iterations: 5 
## n= 6799
```

---

### Fitting an AFT model

```r
library(survival)
aft_n <- survreg(Surv(DWHFDAYS, DWHF) ~ TRTMT + EJFPER + PREVMI, 
                 data = dig, dist = "lognormal")
summary(aft_n)
```

```
## 
## Call:
## survreg(formula = Surv(DWHFDAYS, DWHF) ~ TRTMT + EJFPER + PREVMI, 
##     data = dig, dist = "lognormal")
##                Value Std. Error     z       p
## (Intercept)  6.06283    0.12297 49.30 < 2e-16
## TRTMT        0.55329    0.06741  8.21 2.2e-16
## EJFPER       0.05932    0.00385 15.42 < 2e-16
## PREVMI      -0.11991    0.07059 -1.70   0.089
## Log(scale)   0.82037    0.01636 50.14 < 2e-16
## 
## Scale= 2.27 
## 
## Log Normal distribution
## Loglik(model)= -20343.2   Loglik(intercept only)= -20501
## 	Chisq= 315.72 on 3 degrees of freedom, p= 3.9e-68 
## Number of Newton-Raphson Iterations: 3 
## n= 6799
```

---

### Fitting an AFT model

```r
library(survival)
aft_l <- survreg(Surv(DWHFDAYS, DWHF) ~ TRTMT + EJFPER + PREVMI, 
                 data = dig, dist = "loglogistic")
summary(aft_l)
```

```
## 
## Call:
## survreg(formula = Surv(DWHFDAYS, DWHF) ~ TRTMT + EJFPER + PREVMI, 
##     data = dig, dist = "loglogistic")
##                Value Std. Error     z       p
## (Intercept)  6.04825    0.11365 53.22 < 2e-16
## TRTMT        0.48004    0.06243  7.69 1.5e-14
## EJFPER       0.05570    0.00359 15.50 < 2e-16
## PREVMI      -0.08625    0.06517 -1.32    0.19
## Log(scale)   0.19883    0.01835 10.84 < 2e-16
## 
## Scale= 1.22 
## 
## Log logistic distribution
## Loglik(model)= -20342   Loglik(intercept only)= -20498.1
## 	Chisq= 312.28 on 3 degrees of freedom, p= 2.2e-67 
## Number of Newton-Raphson Iterations: 4 
## n= 6799
```

---

### How well are we fitting?

In each of these past examples, we've assumed a distribution of survival times.
How can we check if we've made an appropriate assumption? Suppose we use the
following residuals:

`\begin{align*}
\widehat{e}_i &= \frac{\log(Y_i) - \mathbf{X}_i^T\widehat{\beta}}{\widehat{\phi}}
\end{align*}`

where `$\phi$` is the scale parameter from the model (don't worry about this).

.question[
What does `$\widehat{e}_i$` represent? Is it appropriate to use these in some way?
]

---

### How well are we fitting?

`$\widehat{e}_i$` examine scaled differences from `$Y_i$`, the *observed* data. 
Thus, these are *censored* residuals, which we must take into account.

One way to deal with this issue is to plot a Kaplan-Meier plot of these 
residuals against the survival function of the assumed distribution.

```r
resids <- (log(dig$DWHFDAYS) - aft_ll$linear.predictors) / 
  (aft_ll$scale)
m1 <- survfit(Surv(resids, DWHF) ~ 1, data = dig)
plot(m1, xlab = "AFT Residuals (Logistic model)", 
     ylab = "Survival Probability")
exp.x <- seq(min(resids), max(resids), length = 100)
exp.y <- plogis(exp.x, lower.tail = F) # F(t)
lines(exp.x, exp.y, col = "red", lwd = 2)
```

---

### How well are we fitting?

---

### Visualizing the AFT model

---

### A similar model for hazards?

---

### A proportional hazards model

`\begin{align*}
\lambda(t) &= \lambda_0(t)\exp\left(\beta_1x_1 + \beta_2x_2 + \cdots + \beta_px_p \right)
\end{align*}`

where `$\lambda_0$` is the .vocab[baseline hazard], which describes the hazard
function for an individual with `$\mathbf{X} = \mathbf{0}$`.

A .vocab[proportional hazards] model is a linear model for the log hazard ratio.

.question[
How might we interpret `$\beta$` terms in such a model? Why is there no intercept
term `$\beta_0$`?
]

---

### A proportional hazards model

`\begin{align*}
\lambda(t) &= \lambda_0(t)\exp\left(\beta_1x_1 + \beta_2x_2 + \cdots + \beta_px_p \right)
\end{align*}`

Covariates in a PH model have a multiplicative effect on *hazards*. For
instance, if `$\beta_k = 0.4$`, then `$\exp(\beta_k) \approx 1.5$`. Holding all else 
equal, an individual with covariate `$x_k$` one unit greater than another is 
expected to have approximately 1.5 times the *hazard* of the event than the
other. This multiplicative effect holds *regardless* of the time `$t$`.

.question[
Would you rather have a higher linear predictor or a lower linear predictor
in a PH model? How does this compare to the AFT model?
]

---

### A proportional hazards model

---

### Fitting a PH model

```r
library(eha)
ph_ln <- phreg(Surv(DWHFDAYS, DWHF) ~ TRTMT + EJFPER + PREVMI, 
               data = dig, dist = "lognormal")
ph_ln
```

```
## Call:
## phreg(formula = Surv(DWHFDAYS, DWHF) ~ TRTMT + EJFPER + PREVMI, 
##     data = dig, dist = "lognormal")
## 
## Covariate          W.mean      Coef Exp(Coef)  se(Coef)    Wald p
## (Intercept)                   4.480               1.109     0.000 
## TRTMT               0.518    -0.290     0.748     0.042     0.000 
## EJFPER             29.455    -0.038     0.963     0.002     0.000 
## PREVMI              0.645     0.048     1.049     0.044     0.271 
## 
## log(scale)                   15.421               2.721     0.000 
## log(shape)                   -1.358               0.123     0.000 
## 
## Events                    2332 
## Total time at risk        6092212 
## Max. log. likelihood      -20339 
## LR test statistic         303.79 
## Degrees of freedom        3 
## Overall p-value           0
```

.vocab[
How might we interpret the coefficient estimates here? How do they relate to
our previous models?
]

---

### The Cox proportional hazards model

`\begin{align*}
\lambda(t) &= \lambda_0(t)\exp\left(\beta_1x_1 + \beta_2x_2 + \cdots + \beta_px_p \right)
\end{align*}`

In a parametric survival model (such as ones we've seen), the survival times
are assumed to follow a specific distribution, which might be a strong
assumption.

Often times, we only care about the `$\beta$` terms and not so much the 
`$\lambda_0$`. Using the concept of partial likelihood, Cox (1972) found that we 
can "separate" inference for the `$\beta$` terms from `$\lambda_0$`.

The Cox model is a *semiparametric* survival model; `$\lambda_0(t)$` (*infinite*-
dimensional) is left completely unspecified, with no assumptions made on its 
shape. (A semi-parametric version of the AFT model also exists, but isn't very 
popular for reasons we won't get into).

---

### Estimating `$\beta$`s

.question[
How can we estimate `$\beta$` parameters if `$\lambda_0(t)$` is infinite-dimensional?
]

---

### Estimating `$\beta$`s

For now, assume that there are no tied failure times. Suppose we know that there
is one failure at time `$t_j$`. For this failure time, the contribution to the
"likelihood" is the probability that individual `$j$` fails, *given* that there is exactly one failure in the risk set at that time:

`\begin{align*}
P(\text{individual }j\text{ fails} | \text{one failure from }\mathcal{R}(t_j))
\end{align*}`

where `$\mathcal{R}(t_j) = \left\{i: T_i \ge t_j \right\}$` is the risk set at time
`$t_j$`.

---

### Estimating `$\beta$`s

Note that the probability

`\begin{align*}
P(\text{individual }j\text{ fails} | \text{one failure from }\mathcal{R}(t_j))
\end{align*}`

conditions on exactly one failure at time `$t_j$`; we don't know who it is! This
is equivalent to the probability

`\begin{align*}
&\mathrel{\phantom{=}} \frac{P(\text{individual }j\text{ fails} | \text{at risk at }t_j)}{\sum_{k \in \mathcal{R}(t_j)}P(\text{individual }k\text{ fails} | \text{at risk at }t_j)}\\
&= \frac{\lambda(t_j | \mathbf{X}_j)}{\sum_{k \in \mathcal{R}(t_j)}\lambda(t_j | \mathbf{X}_k)}
\end{align*}`

---

### Estimating `$\beta$`s

Under proportional hazards,
`$\lambda(t | \mathbf{X}) = \lambda_0(t)\exp(\mathbf{X}^T\beta)$`. Thus, for each
failure time `$t_j$`, we have

`\begin{align*}
&\mathrel{\phantom{=}} \frac{\lambda(t_j | \mathbf{X}_j)}{\sum_{k \in \mathcal{R}(t_j)}\lambda(t_j | \mathbf{X}_k)}\\
&= \frac{\lambda_0(t_j)\exp(\mathbf{X}_j^T\beta)}{\sum_{k \in \mathcal{R}(t_j)}\lambda_0(t_j)\exp(\mathbf{X}_k^T\beta)}\\
&= \frac{\exp(\mathbf{X}_j^T\beta)}{\sum_{k \in \mathcal{R}(t_j)}\exp(\mathbf{X}_k^T\beta)}
\end{align*}`

---

### Estimating `$\beta$`s

If there are `$J$` total failure times, the "likelihood" is thus the product over
all possible failure times:

`\begin{align*}
\mathcal{L}^\star(\beta | \mathbf{X}) &= \prod_{j = 1}^J \frac{\exp(\mathbf{X}_j^T\beta)}{\sum_{k \in \mathcal{R}(t_j)}\exp(\mathbf{X}_k^T\beta)}
\end{align*}`

Note that we can maximize this quantity with respect to `$\beta$` terms *without* 
having to estimate (or even specify) the baseline hazard `$\lambda_0(t)$`!

---

### The Cox proportional hazards model

`\begin{align*}
\lambda(t) &= \lambda_0(t)\exp\left(\beta_1x_1 + \beta_2x_2 + \cdots + \beta_px_p \right)
\end{align*}`

- Fewer assumptions than fully parametric models, but still requires PH 
assumption
- By far the most commonly used regression model for survival data
- Attractive interpretations using hazard ratios
- Can be extended to accommodate time-varying covariates, recurrent events, etc.

---

### Fitting the Cox PH model

```r
coxm1 <- coxph(Surv(DWHFDAYS, DWHF) ~ TRTMT + EJFPER + PREVMI, 
                    data = dig)
summary(coxm1)
```

```
## Call:
## coxph(formula = Surv(DWHFDAYS, DWHF) ~ TRTMT + EJFPER + PREVMI, 
##     data = dig)
## 
##   n= 6799, number of events= 2332 
## 
##             coef exp(coef)  se(coef)       z Pr(>|z|)    
## TRTMT  -0.289745  0.748454  0.041666  -6.954 3.55e-12 ***
## EJFPER -0.038061  0.962654  0.002381 -15.985  < 2e-16 ***
## PREVMI  0.047921  1.049087  0.043638   1.098    0.272    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##        exp(coef) exp(-coef) lower .95 upper .95
## TRTMT     0.7485     1.3361    0.6898    0.8121
## EJFPER    0.9627     1.0388    0.9582    0.9672
## PREVMI    1.0491     0.9532    0.9631    1.1428
## 
## Concordance= 0.605  (se = 0.006 )
## Likelihood ratio test= 304.4  on 3 df,   p=<2e-16
## Wald test            = 302.5  on 3 df,   p=<2e-16
## Score (logrank) test = 306.2  on 3 df,   p=<2e-16
```

.question[
How might we interpret these coefficients? How do they compare to our previous
models?
]

---

### Fitting the Cox PH model

One strength of the Cox model is that we can ignore estimation of `$\lambda_0$`
completely (it doesn't matter for valid inference on the `$\beta$`s).

If we want to estimate survival probabilities, then we must estimate the 
baseline hazard. A non-parametric estimate of the cumulative hazard (the 
.vocab[Breslow estimator]) is implemented by the `basehaz` function in the 
`survival` package by:

`\begin{align*}
\hat{\Lambda}_0(t) = \sum_{i:t_{(i)} < t} \frac{1}{\sum_{k \in \mathcal{R}_i} \exp(\mathbf{X}_k\beta)}
\end{align*}`

We can then estimate survival in the Cox model by:

`\begin{align*}
\hat{S}(t|\mathbf{X}) = \exp(-\exp(\mathbf{X}\beta) \hat{\Lambda}_0(t))
\end{align*}`

---

### Cox model diagnostics

Recall the inverse CDF result: if `$T_i$` has survival function `$S_i(t)$`, then
`$S_i(T_i)$` should have a uniform distribution on (0, 1) and `$\Lambda_i(T_i)$`
should follow Exp(1) distribution.

Thus, if the model is correct, then the estimated cumulative hazard `$\hat{\Lambda}$`
at observed times should be a censored sample from Exp(1). These `$\hat{\Lambda}_i(Y_i)$`
are known as .vocab[Cox-Snell residuals].

Plotting `$\log(-\log\hat{S}(Y_i))$` against `$\log(Y_i)$` should thus follow a 
straight line through the origin at a 45-degree angle.

Although Cox-Snell residuals can also be used for other models (e.g., checking
whether distribution specified in AFT model is appropriate), they're not too
useful in practice (for a variety of reasons).

---

### Cox model diagnostics

To assess PH assumption, we can examine .vocab[Schoenfeld residuals]. 
Intuitively, each Schoenfeld residual is the difference between the observed
covariate and the expected covariate for each observed failure, conditioned on
the risk set at that time.

In plotting Schoenfeld residuals vs. survival times, we expect to see them
randomly distributed around 0.

---

### Cox model diagnostics

.vocab[Martingale residuals] are based on the difference between observed number 
of events up until time `$Y_i$` and the expected count based on the fitted model.

.vocab[Deviance residuals] are a normalized transformation of the martingale
residuals that correct their skewness. They should be randomly distributed 
around mean 0 with a variance of 1.

In practice, these residuals are useful for finding potential outliers: negative
values "lived longer than expected" and positive values "died too soon."

---

### Cox model diagnostics

```r
library(survminer)
ggcoxdiagnostics(coxm1, type = "schoenfeld")
```

![](survival-2_files/figure-html/unnamed-chunk-14-1.png)

---

### Cox model diagnostics

```r
library(survminer)
ggcoxdiagnostics(coxm1, type = "deviance", linear.predictions = F)
```

![](survival-2_files/figure-html/unnamed-chunk-15-1.png)

---

### Bayesian survival analysis

As in all Bayesian models, we must specify our prior information for various
model components.

.question[
What pieces of the model should have priors specified on them? What types of 
prior might be reasonable?
]

- Intercept (baseline hazard)
- Regression "slopes"
- Additional distribution-specific parameters (e.g., shape parameters for 
Weibull)

We can extend models to more flexible situations, but these'll suffice for the
vast majority of time-independent regression situations you'll encounter.