Lab 10: Survival Analysis, Part II

class: center, middle, inverse, title-slide

# Lab 10: Survival Analysis, Part II
## PHS2000A
### Christopher Boyer
### 2021-11-22

---

## Lab Overview

1. What's new?
2. Continuous time-to-event functions
3. Nonparametric estimation and inference
 - Kaplan Meier survival estimator
 - Log-rank test 
 - Other estimators
4. Model-based estimation and inference
 - Cox proportional hazards model
 - Semi-parametric partial likelihood
 - Hypothesis testing
 - Confidence intervals
 - Prediction

---

## What's new?

Last week we covered:
- Discrete time survival methods
- Grouped continuous time survival methods

Question: What is different this week?

---

## Data structure

<table>
 <thead>
  <tr>
   <th style="text-align:right;"> id </th>
   <th style="text-align:right;"> time </th>
   <th style="text-align:right;"> status </th>
   <th style="text-align:right;"> age </th>
   <th style="text-align:right;"> treat </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 59 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 72.33 </td>
   <td style="text-align:right;"> 1 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:right;"> 115 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 74.49 </td>
   <td style="text-align:right;"> 1 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:right;"> 156 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 66.47 </td>
   <td style="text-align:right;"> 1 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:right;"> 421 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 53.36 </td>
   <td style="text-align:right;"> 2 </td>
  </tr>
  <tr>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:right;"> 431 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:right;"> 50.34 </td>
   <td style="text-align:right;"> 1 </td>
  </tr>
</tbody>
</table>

---

## Continuous time-to-event functions
Recall, we've learned several ways to describe the distribution of a continuous time-to-event random variable `$T_i$`:
 - `$f(t)$`, the *probability density function*.
 - `$F(t)$`, the *cumulative distribution function*.
 - `$S(t)=1-F(t)$`, the *survivor function*.
 - `$h(t)=\frac{f(t)}{S(t)}$`, the instantaneous *hazard function*.
 - `$H(t) = \int_0^t h(u) du$`, the *cumulative hazard function*.

Where the hazard in continuous time is defined as
`$$h(t)=\lim_{\Delta t \rightarrow 0} \frac{P(t \le T < t+ \Delta t | T \ge t)}{\Delta t} = \frac{f(t)}{S(t)}$$`

---

## Continuous time-to-event functions
We mentioned that, for a continuous time to event random variable,
`$$h(t) = \frac{f(t)}{S(t)}$$`
It can also be shown that 
`$$h(t) = - \frac{d}{dt}\{ \mbox{log}[S(t)] \}$$`
and so
`$$S(t) = \mbox{exp}\{ - H(t) \}$$`

Note that the cumulative hazard can be obtained from the survivor function, since
`$$H(t) = -\mbox{log}[S(t)]$$`

---

## Kaplan-Meier estimator
The **Kaplan-Meier estimator** of the survival probability `$S(t) = \Pr(T>t)$` is:

`$$\hat{S}_{KM}(t) = \prod_{j: \tau_j \leq t} \left( 1 - \frac{d_j}{r_j} \right)$$`

where

- `$\{ \tau_1, \dots, \tau_K \}$` is the set of `$K$` distinct event times observed in the sample
 - `$d_j$` is the number of deaths at `$\tau_j$`
 - `$r_j$` is the number of individuals *"at risk"* right before the `$j$`th event time
 - `$c_j$` is the number of censored observations between the `$j$`th and `$(j+1)$`th death times. Censorings tied at `$\tau_j$` are included in `$c_j$`.
 - Note that the risk set `$r_j$` can be calculated as either `$r_{j-1} - d_{j-1} - c_{j-1}$` or alternatively as `$\sum_{l \geq j} (c_l + d_l)$`.

The estimator is **non-parametric**, i.e. we don't make any assumptions about the distribution of `$T$`. For `$\hat{S}_{KM}(t)$` to converge to true survival we just need our samples to be *independent* and censoring to be *non-informative*.

---

## Kaplan-Meier estimator

**Example:** We wish to estimate survival for a sample of 12 individuals with the following observed event times:

$$6^+, 6, 6, 6, 7, 9^+, 10^+, 10, 11^+, 13, 16, 17^+ $$
Note that times with `$^+$` are right-censored.

---

## Kaplan-Meier estimator

**Example:** Let's make a table with a row for every rank-ordered death or censoring time:

$$6^+, 6, 6, 6, 7, 9^+, 10^+, 10, 11^+, 13, 16, 17^+ $$

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> time $\tau_j$ </th>
   <th style="text-align:left;"> events $d_j$ </th>
   <th style="text-align:left;"> censored $c_j$ </th>
   <th style="text-align:left;"> n at risk $r_j$ </th>
   <th style="text-align:left;"> $1-(d_j/r_j)$ </th>
   <th style="text-align:left;"> $\hat{S}(\tau_j)$ </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> $6$ </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $7$ </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $9$ </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $10$ </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $11$ </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $13$ </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $16$ </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $17$ </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
  </tr>
</tbody>
</table>

---

## Kaplan-Meier estimator

**Example:** Let's make a table with a row for every rank-ordered death or censoring time:

$$6^+, 6, 6, 6, 7, 9^+, 10^+, 10, 11^+, 13, 16, 17^+ $$

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> time $\tau_j$ </th>
   <th style="text-align:left;"> events $d_j$ </th>
   <th style="text-align:left;"> censored $c_j$ </th>
   <th style="text-align:left;"> n at risk $r_j$ </th>
   <th style="text-align:left;"> $1-(d_j/r_j)$ </th>
   <th style="text-align:left;"> $\hat{S}(\tau_j)$ </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> $6$ </td>
   <td style="text-align:left;"> $3$ </td>
   <td style="text-align:left;"> $1$ </td>
   <td style="text-align:left;"> $12$ </td>
   <td style="text-align:left;"> $\frac{9}{12}=0.750$ </td>
   <td style="text-align:left;"> $0.750$ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $7$ </td>
   <td style="text-align:left;"> $1$ </td>
   <td style="text-align:left;"> $0$ </td>
   <td style="text-align:left;"> $8$ </td>
   <td style="text-align:left;"> $\frac{7}{8}=0.875$ </td>
   <td style="text-align:left;"> $0.656$ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $9$ </td>
   <td style="text-align:left;"> $0$ </td>
   <td style="text-align:left;"> $1$ </td>
   <td style="text-align:left;"> $7$ </td>
   <td style="text-align:left;"> $\frac{7}{7}=1.000$ </td>
   <td style="text-align:left;"> $0.656$ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $10$ </td>
   <td style="text-align:left;"> $1$ </td>
   <td style="text-align:left;"> $1$ </td>
   <td style="text-align:left;"> $6$ </td>
   <td style="text-align:left;"> $\frac{5}{6}=0.833$ </td>
   <td style="text-align:left;"> $0.547$ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $11$ </td>
   <td style="text-align:left;"> $0$ </td>
   <td style="text-align:left;"> $1$ </td>
   <td style="text-align:left;"> $4$ </td>
   <td style="text-align:left;"> $\frac{4}{4}=1.000$ </td>
   <td style="text-align:left;"> $0.547$ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $13$ </td>
   <td style="text-align:left;"> $1$ </td>
   <td style="text-align:left;"> $0$ </td>
   <td style="text-align:left;"> $3$ </td>
   <td style="text-align:left;"> $\frac{2}{3}=0.666$ </td>
   <td style="text-align:left;"> $0.365$ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $16$ </td>
   <td style="text-align:left;"> $1$ </td>
   <td style="text-align:left;"> $0$ </td>
   <td style="text-align:left;"> $2$ </td>
   <td style="text-align:left;"> $\frac{1}{2}=0.500$ </td>
   <td style="text-align:left;"> $0.182$ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $17$ </td>
   <td style="text-align:left;"> $0$ </td>
   <td style="text-align:left;"> $1$ </td>
   <td style="text-align:left;"> $1$ </td>
   <td style="text-align:left;"> $\frac{1}{1}=1.000$ </td>
   <td style="text-align:left;"> $0.182$ </td>
  </tr>
</tbody>
</table>

---

## Kaplan-Meier estimator

Note that:
 - `$\hat{S}(t)$` is 1 up to the first death time.
 - `$\hat{S}(t)$` only changes at death (failure) times. At times where there is only censoring, `$\hat{S}(t)$` stays the same as the preceding interval.
 - `$\hat{S}(t)$` only goes to 0 if the last event is a death. (In this example the last five events are censorings, so `$\hat{S}(t)$` stays the same as it was at the last failure at `$\tau_j=17$`).

---

## Kaplan-Meier estimator

---

## Confidence intervals

We can construct asymptotic point-wise confidence intervals for our survival estimates based on a normal approximation.
`$$\hat{S}_{KM}(t) \pm Z_{1-\alpha/2}\sqrt{\widehat{\text{Var}}[\hat{S}_{KM}(t)]}$$`
where the asypmtotic sampling variance can be calculated using Greenwood's formula, i.e. 
`$$\widehat{\mbox{Var}}(\hat{S}_{KM}(t)) = [\hat{S}_{KM}(t)]^2 \sum_{j:\tau_{j} \leq t} \frac{d_{j}}{(r_j - d_j)r_j}$$`
Note that this approach can sometimes yield values `$>1$` or `$<0$`. A better approach is to get a 95% confidence limit for `$\text{log}[-\text{log}(S(t))]$`.

Also note, one can also construct confidence *bands*, i.e. a band for entire function of `$\hat{S}_{KM}(t)$` such that, across repeated sampling, in `$100(1 - \alpha)$`% of samples the entire true function will lie within the band. It is usually *wider* than the point-wise intervals.

---

## Confidence intervals

---

## Alternative estimators

While we cannot estimate the continuous hazard function, there are two ways to estimate the **cumulative hazard function**.

- The Nelson-Aalen method
 
`$$\hat{H}_{NA}(t) = \sum_{t_j \leq t} \frac{d_j}{r_j}$$`

- The negative log survivor function method takes advantage of the relationship we learned between the cumulative hazard and the survivor function:
 
`$$\hat{H}_{KM}(t) = -\text{log} \left[\hat{S}_{KM}(t) \right]$$`

- The Nelson-Aalen method also leads us to another estimator of `$S(t)$` called the Fleming-Harrington estimator:

`$$\hat{S}_{FH}(t) = \text{exp}(-\hat{H}_{NA}(t))$$`

---

## Alternative estimators

---

## Nonparametric comparison of survival curves

Often we want to compare the survival estimates between two or more groups.
How might we compare the two estimated survival distributions, `$\hat{S}_{A}(t)$` and `$\hat{S}_{B}(t)$`?

If we were interested in comparing survival at a fixed time point `$t^*$`, this would be relatively straightforward. We could evaluate this using our data by checking whether the point-wise confidence intervals we calculated previously for the survival curves overlap at `$t^*$`, or equivalently we can check whether the 95% CI for the *difference in survival estimates* at time `$t^*$`
includes 0.

`$$\left[ \hat{S}_{A}(t^*) - \hat{S}_{B}(t^*) \right] \pm Z_{1-\alpha/2} \times \sqrt{Var(\hat{S}_{A}(t^*) - \hat{S}_{B}(t^*))}$$`

But often we want to know whether survival differs across the full distribution, i.e. are `$\hat{S}_{A}(t)$` and `$\hat{S}_{B}(t)$` meaningfully different at any point?

---

## Log-rank test
There are several non-parametric methods for comparing the survival distribution for two groups over the whole distribution. The most well known and widely used is the Mantel Haenszel **log-rank test**.

`$$H_0: S_A(t) = S_B(t)$$`
`$$H_1: S_A(t) \neq S_B(t)$$`

The log-rank test statistic is:

`$$\chi^{2}_{log-rank} = \frac{\left[ \sum^{K}_{j=1} (d_{Aj} - r_{Aj} \times d_{j}/r_{j})\right]^{2}} {\sum^{K}_{j=1} \frac{r_{Bj}r_{Aj} d_{j}(r_j - d_j)}{[r_j^2(r_j-1)]}}$$`
which has an asymptotic `$\chi^{2}$` distribution with 1 df under the null.

---

## Log-rank test

The log-rank test can be obtained by constructing a `$(2 \times 2)$` table at each
distinct death time and comparing the death rates between the two groups,
conditional on the number at risk in the group. The tables are then combined
using the Cochran-Mantel-Haenszel test.

Let `$\tau_1, \dots, \tau_K$` represent the ordered, distinct death times.
At the `$j$`th death time, we have the following table:

<table>
 <thead>
  <tr>
   <th style="text-align:left;"> Group </th>
   <th style="text-align:left;"> Died </th>
   <th style="text-align:left;"> Survived </th>
   <th style="text-align:left;"> Total </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> $A$ </td>
   <td style="text-align:left;"> $d_{Aj}$ </td>
   <td style="text-align:left;"> $r_{Aj} - d_{Aj}$ </td>
   <td style="text-align:left;"> $r_{Aj}$ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> $B$ </td>
   <td style="text-align:left;"> $d_{Bj}$ </td>
   <td style="text-align:left;"> $r_{Bj} - d_{Bj}$ </td>
   <td style="text-align:left;"> $r_{Bj}$ </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Total </td>
   <td style="text-align:left;"> $d_{j}$ </td>
   <td style="text-align:left;"> $r_{j} - d_{j}$ </td>
   <td style="text-align:left;"> $r_{j}$ </td>
  </tr>
</tbody>
</table>

where `$d_{0j}$` and `$d_{1j}$` are the number of deaths in groups 0 and 1, respectively,
at the `$j$`th death time, and `$r_{0j}$` and `$r_{1j}$` are the corresponding number at risk
at that time.

Assuming the tables are all independent, then this statistic will have an approximate `$\chi^{2}$` distribution with 1 df.

---

## Exercise

Load the `ovarian` dataset from the `survival` package in `R`.

1. Calculate the Kaplan-Meier estimates of survival by hand for the 13 women in the treatment group.
2. Confirm your results using the `survfit()` function in `R`.
3. Calculate the Flemming-Harrington estimates of survival by hand for the 13 women in the treatment group.
4. Confirm your results using the `survfit()` function in `R`.
5. Plot both estimates for treatment and control groups with 95% confidence intervals (use `log-log` variance).
6. Conduct a log-rank test comparing the survival estimates across treatment and control. How do you interpret the results?

---

## Exercise

```r
km <- survfit(Surv(futime, fustat) ~ rx, data = ovarian)
summary(km)
```

```
## Call: survfit(formula = Surv(futime, fustat) ~ rx, data = ovarian)
## 
##                 rx=1 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    59     13       1    0.923  0.0739        0.789        1.000
##   115     12       1    0.846  0.1001        0.671        1.000
##   156     11       1    0.769  0.1169        0.571        1.000
##   268     10       1    0.692  0.1280        0.482        0.995
##   329      9       1    0.615  0.1349        0.400        0.946
##   431      8       1    0.538  0.1383        0.326        0.891
##   638      5       1    0.431  0.1467        0.221        0.840
## 
##                 rx=2 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   353     13       1    0.923  0.0739        0.789        1.000
##   365     12       1    0.846  0.1001        0.671        1.000
##   464      9       1    0.752  0.1256        0.542        1.000
##   475      8       1    0.658  0.1407        0.433        1.000
##   563      7       1    0.564  0.1488        0.336        0.946
```

---

## Exercise

```r
fh <- survfit(Surv(futime, fustat) ~ rx, data = ovarian, stype = 2)
summary(fh)
```

```
## Call: survfit(formula = Surv(futime, fustat) ~ rx, data = ovarian, 
##     stype = 2)
## 
##                 rx=1 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    59     13       1    0.926  0.0712        0.796        1.000
##   115     12       1    0.852  0.0966        0.682        1.000
##   156     11       1    0.778  0.1131        0.585        1.000
##   268     10       1    0.704  0.1242        0.498        0.995
##   329      9       1    0.630  0.1313        0.419        0.948
##   431      8       1    0.556  0.1351        0.345        0.895
##   638      5       1    0.455  0.1433        0.246        0.843
## 
##                 rx=2 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   353     13       1    0.926  0.0712        0.796        1.000
##   365     12       1    0.852  0.0966        0.682        1.000
##   464      9       1    0.762  0.1210        0.558        1.000
##   475      8       1    0.673  0.1359        0.453        1.000
##   563      7       1    0.583  0.1443        0.359        0.947
```

---

## Exercise

```r
km <- survfit(Surv(futime, fustat) ~ rx, data = ovarian, conf.type = "log-log")
autoplot(km)
```

---

## Exercise

```r
fh <- survfit(Surv(futime, fustat) ~ rx, data = ovarian, conf.type = "log-log", stype = 2)
autoplot(fh) 
```

---

## Exercise

```r
lrank <- survdiff(Surv(futime, fustat) ~ rx, data = ovarian, rho = 0)
lrank
```

```
## Call:
## survdiff(formula = Surv(futime, fustat) ~ rx, data = ovarian, 
##     rho = 0)
## 
##       N Observed Expected (O-E)^2/E (O-E)^2/V
## rx=1 13        7     5.23     0.596      1.06
## rx=2 13        5     6.77     0.461      1.06
## 
##  Chisq= 1.1  on 1 degrees of freedom, p= 0.3
```

---

## Cox proportional hazards model
Cox defined the following model for the hazard at time `$T = t$` given covariates `$\mathbf{X}$`:

`$$h(t \mid \mathbf{x}) = h_0(t) \exp\{\mathbf{x}^T\beta\}$$`

where `$h_0(t)$` is the baseline hazard and `$e^{\beta_1x_1 +\beta_2x_2 + \dots + \beta_p x_p}$` is a multiplicative factor (the hazard ratio) by which we multiply the baseline hazard to reflect the effect of covariates.

The Cox model is probably the most popular way of modeling time to event data. Here, we focus on modeling the *log of the hazard ratio* and treat the baseline hazard as a *nuisance parameter*.

what's different: no distributional component.

---

## The partial likelihood
Under the proportional hazard assumption, `$h(t\mid\mathbf{x}) = h_0(t) e^{\mathbf{x}^T \beta}$`, so we get
`$$L^{partial}(\beta) = \prod_{j=1}^{K} \frac{h_0(\tau_j) e^{\mathbf{x}^T_j \beta}}
{\sum_{\ell \in \mathcal{R}(\tau_j)} h_0(\tau_j)e^{\mathbf{x}^T_{\ell}\beta }}
= \prod^{K}_{j=1} \frac{e^{\mathbf{x}^T_j \beta}}{\sum_{\ell \in \mathcal{R}(\tau_j)} e^{\mathbf{x}^T_{\ell} \beta}}$$`
Because each individual's contribution to the partial likelihood is computed when he or she experiences the event, only individuals with non-censored event times have an explicit term. However, everyone contributes indirectly to the denominator of the explicit contributions for any non-censored individual whose event time is less than or equal to their observed (or censored) event time.

---

## Estimation
As usual, we take the logarithm of the (partial) likelihood function and obtain:

$$
`\begin{aligned}
\ell(\beta) & = \log \left[ \prod^{K}_{j=1} \frac{e^{\mathbf{x}^T_{j} \beta}}{\sum_{\ell \in \mathcal{R}(\tau_j)} e^{\mathbf{x}^T_\ell \beta}} \right] \\
& = \sum^{K}_{j=1} \left[ \mathbf{x}^T_j \beta - \log \left( \sum_{\ell \in \mathcal{R}(\tau_j)} e^{\mathbf{x}^T_\ell \beta} \right) \right]
\end{aligned}`
$$
The partial likelihood score equations (i.e. the derivatives of the partial log likelihood with respect to `$\beta$`) are:
`$$U(\beta) = \frac{\partial}{\partial \beta} \ell(\beta) =
\sum^{K}_{j=1} \left \{ \mathbf{x}_{j} - \frac{\sum_{\ell \in \mathcal{R}(\tau_j)} x_\ell e^{\mathbf{x}^T_{j} \beta}}
{\sum_{\ell \in \mathcal{R}(\tau_j)} e^{\mathbf{x}^T_\ell \beta}} \right \}$$`
The maximum partial likelihood estimator (MPLE), `$\hat{\beta}$` can be found by solving `$U(\beta)=0$`. This is done numerically using similar methods as for GLMs.

---

## What is a semi-parametric model?

As it sounds a semi-parametric model is somewhere between a nonparametric and a parametric model. If a nonparametric model makes no assumptions about the shape of the distribution and a parametric model specifies it exactly, a semi-parametric model splits the difference by restricting to a certain class of distributions, which are indexed by parameters that are infinite dimensional and therefore extremely flexible.

The **Cox proportional hazard model** is a semi-parametric model, because:
 - it restricts to class of distributions where hazards are proportional.
 - but it allows the baseline hazard `$h_0(t)$` to be completely arbitrary... it could be anything (i.e. infinite dimensional).

This means:
 - The shape of the baseline hazard function is irrelevant.
 - Unlike glm's we have no formal distributional component (other than the proportional hazards assumption)
 - The downside is that the fitted model does not provide predicted values of the hazard, which sometimes we would like to know about. The best we can do is to use nonparametric strategies to "recover" survivor and cumulative hazard functions.

---

## Interpretation
The form of a proportional hazard model is
`$$h(t \mid \mathbf{x}) = h_0(t) e^{\beta_1x_1 + \beta_2x_2 + \dots + \beta_p x_p}$$`
which we can re-write as
`$$\mbox{log}\left[ \frac{h(t \mid \mathbf{x})}{h_0(t)} \right] = \mbox{log}\left[HR\right]= \beta_1x_1 + \beta_2x_2 + \dots + \beta_p x_p$$`
So this is a linear model for the *log of the hazard ratio*.

- How do we interpret `$\hat{\beta}_1$`, `$\hat{\beta}_2$`, etc.?
 - How do we interpret `$\mbox{exp}(\hat{\beta}_1)$`, `$\mbox{exp}(\hat{\beta}_2)$`,  etc.?

---

## Hypothesis testing
For each covariate of interest, consider the following hypothesis test:
`$$H_0 : \beta_j = 0 \iff H_0 :\mbox{HR}_j = 1$$`
`$$H_1 : \beta_j \neq 0 \iff H_1 :\mbox{HR}_j \neq 1$$`

A Wald test of the above hypothesis is constructed as
`$$Z = \frac{\hat{\beta}_j}{se(\hat{\beta}_j)} \quad or \quad \chi^2 = \left( \frac{\hat{\beta}_j}{se(\hat{\beta}_j)}\right)^2$$`
where the first statistic, `$Z$`, follows a normal distribution and the second, `$\chi^2$`, follows a Chi-quare distribution with 1 df. Different software will give different versions of the Wald test. For both, the p-values are given and don't depend on which form, `$Z$` or `$\chi^2$`, is provided.

We can also conduct likelihood ratio and score tests (the latter is equivalent to log-rank test).

---

## Confidence intervals
Again, we can construct asymptotic confidence intervals for `$\hat{\beta}$` based on a normal approximation on the log-scale.
`$$\hat{\beta} \pm 1.96se(\hat{\beta})$$`
Wc can exponentiate the endpoints to get a confidence interval for the hazard ratio.
`$$\left(e^{\hat{\beta}-1.96se(\hat{\beta})}, e^{\hat{\beta}+1.96se(\hat{\beta})}\right)$$`

---

## Prediction

The Cox PH model says that `$h_i(t \mid \mathbf{x}) = h_0(t) \exp(\mathbf{x}^T\beta)$`. What does this imply about the survivor function, `$S_i(t)$`, for the i-th individual with covariates `$\mathbf{x}_i$`?

For the baseline (reference) group, we have
`$$S_0(t) = e^{-\int_0^t h_0(u) du} = e^{-H_0(t)}$$`

For the i-th individual with covariates `$\mathbf{x}_i$`, we have

$$
`\begin{aligned}
S_i(t) & = e^{-\int_0^t h_i \; du}  \\
& = e^{-\int_0^t h_0(u) \exp(\mathbf{x}_i^T \beta) \; du} \\
& = e^{-\exp(\mathbf{x}_i^T \beta) \int_0^t h_0(u) \; du} \\
& = \left[ e^{-\int_0^t h_0(u) \; du} \right]^{\exp(\mathbf{x}^T_i \beta)} \quad \mbox{since }[e^b]^a = e^{ab} \\
& = \left[S_0(t)\right]^{\exp(\mathbf{x}_i^T \beta)} \\
& = \left[S_0(t)\right]^{HR_i}
\end{aligned}`
$$

---

## Prediction

So after estimating `$\hat{\beta}$`, we can obtain an estimate for `$S_i(t)$` as
`$$\hat{S}_i(t) = \left[S_0(t)\right]^{\widehat{HR}_i}$$`

How can we estimate the baseline survivor function, `$S_0(t)$`? We could use the KM estimator, but

Instead, we will use a baseline hazard estimator which takes advantage of the proportional hazards assumption to get a smoother estimate.
`$$\hat{S}_i(t) = \left[ \hat{S}_0(t) \right]^{\exp(\mathbf{x}^T_i \hat{\beta})}$$`
Using the above formula, we substitute `$\hat{\beta}$` based on fitting the Cox PH model and calculate `$\hat{S}_0(t)$` by one of the following approaches:
 - Breslow estimator
 - Kalbfleisch/Prentice estimator 
 
---

## Exercise

Load the `colon` dataset from the `survival` package in `R`.

1. Calculate the Kaplan-Meier estimates of survival for the treatment groups using the `survfit()` function in `R`. Does the PH assumption seem reasonable here?
2. Fit a Cox proportional hazards model for effect of treatment and adjusting for `age` and `sex`. Interpret the coefficients on the treatment terms.
3. Conduct Wald and likelihood ratio tests for the effect of treatment? What do you conclude?
4. Does the effect of treatment differ by `age` or `sex`? How do you know?
5. What is the predicted survival for a male aged 50 and female aged 50 under the `Lev+5FU` treatment at time 1500?

---

## Exercise

```r
km <- survfit(Surv(time, status) ~ rx, data = colon, conf.type = "log-log")
autoplot(km)
```

---

## Exercise

```r
cox <- coxph(Surv(time, status) ~ rx + sex + age, data = colon)
summary(cox)
```

```
## Call:
## coxph(formula = Surv(time, status) ~ rx + sex + age, data = colon)
## 
##   n= 1858, number of events= 920 
## 
##                coef exp(coef)  se(coef)      z Pr(>|z|)    
## rxLev     -0.018333  0.981834  0.076871 -0.238    0.811    
## rxLev+5FU -0.440908  0.643452  0.083961 -5.251 1.51e-07 ***
## sex       -0.050352  0.950894  0.066079 -0.762    0.446    
## age       -0.002034  0.997968  0.002805 -0.725    0.468    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##           exp(coef) exp(-coef) lower .95 upper .95
## rxLev        0.9818      1.019    0.8445    1.1415
## rxLev+5FU    0.6435      1.554    0.5458    0.7585
## sex          0.9509      1.052    0.8354    1.0824
## age          0.9980      1.002    0.9925    1.0035
## 
## Concordance= 0.548  (se = 0.01 )
## Likelihood ratio test= 36.35  on 4 df,   p=2e-07
## Wald test            = 34.24  on 4 df,   p=7e-07
## Score (logrank) test = 34.75  on 4 df,   p=5e-07
```

---

## Exercise

```r
cox <- coxph(Surv(time, status) ~ rx + sex + age, data = colon)
summary(cox)
```

---

## Exercise

```r
cox <- coxph(Surv(time, status) ~ rx + sex + age, data = colon)
cox.nested <- coxph(Surv(time, status) ~ sex + age, data = colon)
anova(cox, cox.nested)
```

```
## Analysis of Deviance Table
##  Cox model: response is  Surv(time, status)
##  Model 1: ~ rx + sex + age
##  Model 2: ~ sex + age
##    loglik  Chisq Df P(>|Chi|)    
## 1 -6587.8                        
## 2 -6605.4 35.332  2 2.127e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

---

## Exercise

```r
cox <- coxph(Surv(time, status) ~ rx * (sex + age), data = colon)
summary(cox)
```

```
## Call:
## coxph(formula = Surv(time, status) ~ rx * (sex + age), data = colon)
## 
##   n= 1858, number of events= 920 
## 
##                    coef exp(coef)  se(coef)      z Pr(>|z|)   
## rxLev         -0.150268  0.860478  0.415880 -0.361  0.71786   
## rxLev+5FU      0.701454  2.016684  0.416870  1.683  0.09244 . 
## sex            0.061860  1.063813  0.108206  0.572  0.56753   
## age            0.001671  1.001673  0.004602  0.363  0.71648   
## rxLev:sex      0.056175  1.057782  0.155271  0.362  0.71751   
## rxLev+5FU:sex -0.520538  0.594201  0.170831 -3.047  0.00231 **
## rxLev:age      0.001563  1.001565  0.006693  0.234  0.81532   
## rxLev+5FU:age -0.015332  0.984785  0.006839 -2.242  0.02497 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##               exp(coef) exp(-coef) lower .95 upper .95
## rxLev            0.8605     1.1621    0.3808    1.9442
## rxLev+5FU        2.0167     0.4959    0.8908    4.5654
## sex              1.0638     0.9400    0.8605    1.3151
## age              1.0017     0.9983    0.9927    1.0107
## rxLev:sex        1.0578     0.9454    0.7802    1.4340
## rxLev+5FU:sex    0.5942     1.6829    0.4251    0.8305
## rxLev:age        1.0016     0.9984    0.9885    1.0148
## rxLev+5FU:age    0.9848     1.0154    0.9717    0.9981
## 
## Concordance= 0.563  (se = 0.01 )
## Likelihood ratio test= 56.26  on 8 df,   p=3e-09
## Wald test            = 48.75  on 8 df,   p=7e-08
## Score (logrank) test = 50.32  on 8 df,   p=4e-08
```

---

## Exercise

```r
df <- data.frame(
  rx = factor("Lev+5FU", levels = levels(colon$rx)),
  sex = c(1,0), 
  age = 50
)

cox.pred <- survfit(coxph(Surv(time, status) ~ rx * (sex + age), data = colon), newdata = df)
summary(cox.pred)
```

```
## Call: survfit(formula = coxph(Surv(time, status) ~ rx * (sex + age), 
##     data = colon), newdata = df)
## 
##  time n.risk n.event survival1 survival2
##     8   1858       1     1.000     0.999
##     9   1857       1     0.999     0.999
##    19   1856       1     0.999     0.998
##    20   1855       1     0.999     0.998
##    23   1854       1     0.998     0.997
##    24   1852       1     0.998     0.997
##    28   1850       1     0.998     0.996
##    34   1849       1     0.997     0.996
##    35   1848       1     0.997     0.995
##    36   1847       1     0.997     0.995
##    38   1846       1     0.996     0.994
##    40   1845       1     0.996     0.993
##    43   1844       1     0.996     0.993
##    45   1843       2     0.995     0.992
##    49   1840       1     0.994     0.991
##    52   1839       1     0.994     0.991
##    56   1838       1     0.994     0.990
##    59   1836       1     0.993     0.990
##    62   1835       2     0.993     0.989
##    63   1833       1     0.992     0.988
##    68   1832       1     0.992     0.987
##    72   1831       2     0.991     0.986
##    77   1829       2     0.991     0.985
##    78   1827       1     0.990     0.985
##    79   1826       2     0.990     0.984
##    80   1824       3     0.989     0.982
##    85   1821       2     0.988     0.981
##    86   1819       2     0.987     0.980
##    88   1817       1     0.987     0.979
##    91   1816       2     0.986     0.978
##    93   1814       1     0.986     0.978
##    94   1813       1     0.985     0.977
##    98   1812       3     0.984     0.975
##    99   1809       2     0.984     0.974
##   100   1807       1     0.983     0.974
##   101   1806       2     0.983     0.973
##   102   1804       1     0.982     0.972
##   103   1803       1     0.982     0.972
##   105   1802       1     0.982     0.971
##   106   1801       1     0.981     0.971
##   108   1800       1     0.981     0.970
##   109   1799       1     0.981     0.969
##   111   1798       1     0.980     0.969
##   113   1797       4     0.979     0.967
##   116   1793       3     0.978     0.965
##   118   1790       1     0.977     0.965
##   119   1789       1     0.977     0.964
##   121   1788       2     0.976     0.963
##   122   1786       2     0.976     0.962
##   125   1784       1     0.975     0.961
##   127   1783       2     0.975     0.960
##   129   1781       1     0.974     0.960
##   131   1780       1     0.974     0.959
##   132   1779       1     0.974     0.959
##   133   1778       1     0.973     0.958
##   134   1777       1     0.973     0.957
##   136   1776       1     0.973     0.957
##   138   1775       1     0.972     0.956
##   139   1774       1     0.972     0.956
##   141   1773       2     0.971     0.955
##   143   1771       1     0.971     0.954
##   144   1770       1     0.970     0.954
##   145   1769       2     0.970     0.953
##   146   1767       2     0.969     0.951
##   147   1765       1     0.969     0.951
##   150   1764       1     0.968     0.950
##   154   1763       2     0.968     0.949
##   157   1761       3     0.967     0.948
##   160   1758       1     0.966     0.947
##   161   1757       2     0.965     0.946
##   164   1755       1     0.965     0.945
##   165   1754       4     0.964     0.943
##   166   1750       2     0.963     0.942
##   167   1748       1     0.963     0.942
##   168   1747       1     0.962     0.941
##   169   1746       1     0.962     0.940
##   171   1745       2     0.961     0.939
##   173   1743       3     0.960     0.938
##   174   1740       3     0.959     0.936
##   175   1737       1     0.959     0.936
##   176   1736       2     0.958     0.934
##   179   1734       1     0.958     0.934
##   181   1733       1     0.957     0.933
##   183   1732       2     0.957     0.932
##   185   1730       5     0.955     0.930
##   186   1725       2     0.954     0.928
##   187   1723       1     0.954     0.928
##   188   1722       1     0.953     0.927
##   189   1721       2     0.953     0.926
##   191   1719       3     0.952     0.925
##   196   1716       1     0.951     0.924
##   198   1715       1     0.951     0.924
##   199   1714       1     0.951     0.923
##   201   1713       2     0.950     0.922
##   203   1711       1     0.950     0.921
##   204   1710       1     0.949     0.921
##   205   1709       1     0.949     0.920
##   206   1708       1     0.948     0.920
##   208   1707       3     0.947     0.918
##   215   1704       3     0.946     0.916
##   216   1701       1     0.946     0.916
##   218   1700       4     0.945     0.914
##   219   1696       3     0.943     0.912
##   221   1692       1     0.943     0.911
##   222   1691       1     0.943     0.911
##   223   1690       1     0.942     0.910
##   224   1689       1     0.942     0.910
##   226   1688       1     0.942     0.909
##   227   1687       1     0.941     0.909
##   228   1686       1     0.941     0.908
##   229   1685       2     0.940     0.907
##   230   1683       5     0.938     0.904
##   232   1678       1     0.938     0.904
##   235   1677       1     0.938     0.903
##   237   1676       2     0.937     0.902
##   238   1674       3     0.936     0.901
##   241   1671       1     0.936     0.900
##   242   1670       2     0.935     0.899
##   243   1668       1     0.934     0.898
##   245   1667       2     0.934     0.897
##   246   1665       1     0.933     0.897
##   248   1664       1     0.933     0.896
##   250   1663       1     0.933     0.896
##   251   1662       1     0.932     0.895
##   252   1661       1     0.932     0.894
##   253   1660       2     0.931     0.893
##   255   1658       1     0.931     0.893
##   256   1657       2     0.930     0.892
##   257   1655       2     0.929     0.891
##   258   1653       1     0.929     0.890
##   259   1652       2     0.928     0.889
##   260   1650       2     0.928     0.888
##   261   1648       1     0.927     0.887
##   262   1647       1     0.927     0.887
##   263   1646       2     0.926     0.886
##   264   1644       2     0.925     0.885
##   269   1642       1     0.925     0.884
##   271   1641       3     0.924     0.882
##   273   1638       1     0.924     0.882
##   274   1637       2     0.923     0.881
##   275   1635       1     0.923     0.880
##   276   1634       3     0.921     0.879
##   279   1631       4     0.920     0.876
##   280   1627       1     0.920     0.876
##   283   1626       2     0.919     0.875
##   285   1624       1     0.919     0.874
##   286   1623       3     0.917     0.873
##   289   1620       1     0.917     0.872
##   290   1619       2     0.916     0.871
##   291   1617       1     0.916     0.870
##   293   1616       1     0.916     0.870
##   294   1615       2     0.915     0.869
##   296   1613       2     0.914     0.868
##   300   1611       1     0.914     0.867
##   302   1610       1     0.913     0.866
##   303   1609       1     0.913     0.866
##   304   1608       2     0.912     0.865
##   308   1606       1     0.912     0.864
##   311   1605       1     0.912     0.864
##   313   1604       2     0.911     0.863
##   314   1602       2     0.910     0.862
##   315   1600       2     0.909     0.860
##   316   1598       1     0.909     0.860
##   322   1597       3     0.908     0.858
##   323   1594       1     0.908     0.858
##   324   1593       1     0.907     0.857
##   325   1592       1     0.907     0.857
##   326   1591       1     0.906     0.856
##   328   1589       1     0.906     0.856
##   329   1588       1     0.906     0.855
##   330   1587       2     0.905     0.854
##   331   1585       1     0.905     0.853
##   333   1584       1     0.904     0.853
##   334   1583       1     0.904     0.852
##   335   1582       1     0.903     0.852
##   336   1581       2     0.903     0.851
##   337   1579       3     0.902     0.849
##   340   1576       1     0.901     0.848
##   341   1574       1     0.901     0.848
##   342   1573       1     0.901     0.847
##   343   1572       1     0.900     0.847
##   344   1571       1     0.900     0.846
##   348   1570       2     0.899     0.845
##   349   1568       2     0.898     0.844
##   352   1566       2     0.898     0.843
##   354   1564       1     0.897     0.842
##   355   1563       2     0.896     0.841
##   356   1560       3     0.895     0.840
##   360   1557       2     0.895     0.838
##   362   1555       2     0.894     0.837
##   363   1553       1     0.894     0.837
##   365   1552       2     0.893     0.836
##   366   1550       2     0.892     0.835
##   369   1548       1     0.892     0.834
##   370   1547       1     0.891     0.834
##   372   1546       2     0.891     0.832
##   374   1544       1     0.890     0.832
##   376   1543       1     0.890     0.831
##   378   1542       1     0.889     0.831
##   379   1541       1     0.889     0.830
##   380   1540       2     0.888     0.829
##   381   1538       1     0.888     0.829
##   382   1537       2     0.887     0.827
##   384   1535       3     0.886     0.826
##   386   1532       2     0.885     0.825
##   389   1530       2     0.885     0.824
##   390   1528       1     0.884     0.823
##   392   1527       1     0.884     0.823
##   393   1526       1     0.883     0.822
##   398   1525       1     0.883     0.821
##   400   1524       1     0.883     0.821
##   401   1523       1     0.882     0.820
##   402   1522       2     0.882     0.819
##   405   1520       1     0.881     0.819
##   406   1519       2     0.880     0.818
##   408   1517       1     0.880     0.817
##   409   1516       1     0.880     0.816
##   411   1515       2     0.879     0.815
##   413   1513       3     0.878     0.814
##   415   1510       2     0.877     0.813
##   417   1508       1     0.877     0.812
##   420   1507       1     0.876     0.812
##   421   1506       1     0.876     0.811
##   422   1504       2     0.875     0.810
##   428   1502       1     0.875     0.809
##   429   1501       1     0.874     0.809
##   430   1500       2     0.874     0.808
##   431   1498       1     0.873     0.807
##   433   1497       2     0.873     0.806
##   434   1495       1     0.872     0.805
##   435   1494       1     0.872     0.805
##   437   1493       2     0.871     0.804
##   438   1491       3     0.870     0.802
##   439   1488       2     0.869     0.801
##   440   1486       1     0.869     0.800
##   441   1485       1     0.868     0.800
##   443   1484       2     0.868     0.799
##   444   1482       1     0.867     0.798
##   448   1481       2     0.867     0.797
##   449   1479       1     0.866     0.797
##   454   1476       3     0.865     0.795
##   458   1473       2     0.864     0.794
##   459   1471       1     0.864     0.793
##   460   1470       1     0.863     0.793
##   461   1469       1     0.863     0.792
##   462   1468       1     0.863     0.792
##   464   1467       1     0.862     0.791
##   465   1466       3     0.861     0.789
##   466   1463       2     0.860     0.788
##   469   1461       1     0.860     0.788
##   472   1460       1     0.860     0.787
##   474   1459       2     0.859     0.786
##   475   1457       1     0.859     0.786
##   476   1456       1     0.858     0.785
##   480   1455       1     0.858     0.784
##   482   1454       1     0.857     0.784
##   484   1453       1     0.857     0.783
##   485   1452       2     0.856     0.782
##   486   1449       1     0.856     0.782
##   489   1447       1     0.855     0.781
##   490   1446       2     0.855     0.780
##   491   1444       2     0.854     0.779
##   493   1442       1     0.854     0.778
##   495   1441       1     0.853     0.778
##   496   1440       1     0.853     0.777
##   497   1439       1     0.852     0.777
##   498   1438       2     0.852     0.776
##   499   1436       4     0.850     0.773
##   503   1432       1     0.850     0.773
##   504   1431       1     0.849     0.772
##   505   1430       1     0.849     0.772
##   506   1429       1     0.849     0.771
##   510   1428       2     0.848     0.770
##   511   1426       1     0.847     0.770
##   512   1425       1     0.847     0.769
##   513   1424       1     0.847     0.768
##   522   1423       1     0.846     0.768
##   523   1422       1     0.846     0.767
##   525   1421       1     0.845     0.767
##   526   1420       1     0.845     0.766
##   527   1419       1     0.845     0.766
##   528   1418       1     0.844     0.765
##   529   1417       1     0.844     0.765
##   532   1416       2     0.843     0.763
##   534   1414       1     0.843     0.763
##   536   1413       1     0.842     0.762
##   537   1412       1     0.842     0.762
##   540   1411       1     0.842     0.761
##   542   1410       1     0.841     0.761
##   543   1409       2     0.840     0.760
##   546   1407       1     0.840     0.759
##   547   1406       1     0.840     0.758
##   548   1405       1     0.839     0.758
##   550   1404       1     0.839     0.757
##   553   1403       1     0.839     0.757
##   554   1402       2     0.838     0.756
##   555   1400       1     0.837     0.755
##   559   1399       1     0.837     0.755
##   560   1398       1     0.837     0.754
##   561   1397       1     0.836     0.753
##   563   1396       2     0.835     0.752
##   565   1394       1     0.835     0.752
##   569   1393       1     0.835     0.751
##   570   1392       1     0.834     0.751
##   573   1391       3     0.833     0.749
##   576   1388       2     0.832     0.748
##   577   1386       1     0.832     0.747
##   578   1385       3     0.831     0.746
##   580   1382       2     0.830     0.745
##   581   1380       1     0.830     0.744
##   582   1379       1     0.829     0.744
##   583   1378       2     0.828     0.742
##   587   1376       1     0.828     0.742
##   589   1375       1     0.828     0.741
##   591   1374       2     0.827     0.740
##   592   1372       1     0.826     0.740
##   593   1371       3     0.825     0.738
##   594   1368       2     0.824     0.737
##   595   1366       1     0.824     0.736
##   599   1365       2     0.823     0.735
##   601   1362       1     0.823     0.735
##   602   1361       3     0.822     0.733
##   603   1358       1     0.821     0.732
##   604   1357       1     0.821     0.732
##   608   1356       2     0.820     0.731
##   609   1354       1     0.820     0.730
##   612   1353       1     0.819     0.730
##   613   1352       1     0.819     0.729
##   614   1351       1     0.819     0.729
##   615   1350       1     0.818     0.728
##   616   1349       2     0.817     0.727
##   617   1347       1     0.817     0.726
##   622   1346       2     0.816     0.725
##   625   1344       1     0.816     0.725
##   628   1343       1     0.815     0.724
##   629   1341       1     0.815     0.724
##   632   1340       1     0.815     0.723
##   636   1339       1     0.814     0.722
##   638   1338       1     0.814     0.722
##   641   1337       1     0.813     0.721
##   642   1336       2     0.813     0.720
##   643   1334       2     0.812     0.719
##   647   1332       1     0.811     0.719
##   649   1331       1     0.811     0.718
##   653   1330       1     0.811     0.717
##   654   1329       1     0.810     0.717
##   657   1328       1     0.810     0.716
##   659   1327       2     0.809     0.715
##   663   1325       3     0.808     0.714
##   664   1322       1     0.807     0.713
##   665   1321       1     0.807     0.712
##   666   1319       1     0.807     0.712
##   668   1318       1     0.806     0.711
##   669   1317       1     0.806     0.711
##   670   1316       1     0.805     0.710
##   672   1315       1     0.805     0.710
##   673   1314       1     0.805     0.709
##   674   1313       1     0.804     0.709
##   675   1312       2     0.804     0.707
##   678   1310       1     0.803     0.707
##   680   1309       1     0.803     0.706
##   683   1308       1     0.802     0.706
##   684   1307       1     0.802     0.705
##   685   1306       1     0.802     0.705
##   686   1305       1     0.801     0.704
##   687   1304       1     0.801     0.704
##   692   1303       3     0.800     0.702
##   693   1300       1     0.799     0.701
##   696   1299       1     0.799     0.701
##   697   1298       1     0.798     0.700
##   700   1297       2     0.798     0.699
##   701   1295       1     0.797     0.699
##   702   1294       2     0.796     0.697
##   706   1292       1     0.796     0.697
##   708   1291       1     0.796     0.696
##   709   1290       2     0.795     0.695
##   711   1288       1     0.794     0.695
##   712   1287       2     0.793     0.694
##   716   1285       1     0.793     0.693
##   717   1284       2     0.792     0.692
##   718   1282       1     0.792     0.691
##   720   1281       1     0.791     0.691
##   721   1280       1     0.791     0.690
##   723   1279       1     0.791     0.690
##   726   1278       1     0.790     0.689
##   729   1277       1     0.790     0.689
##   730   1276       2     0.789     0.687
##   731   1274       1     0.789     0.687
##   735   1273       1     0.788     0.686
##   736   1272       1     0.788     0.686
##   739   1271       2     0.787     0.685
##   742   1269       1     0.787     0.684
##   743   1268       2     0.786     0.683
##   748   1266       1     0.785     0.682
##   751   1265       1     0.785     0.682
##   752   1264       1     0.785     0.681
##   753   1263       1     0.784     0.681
##   755   1262       1     0.784     0.680
##   758   1261       1     0.783     0.680
##   759   1260       2     0.783     0.679
##   760   1258       1     0.782     0.678
##   761   1257       1     0.782     0.677
##   764   1256       1     0.781     0.677
##   765   1255       1     0.781     0.676
##   766   1254       1     0.781     0.676
##   770   1253       1     0.780     0.675
##   772   1252       1     0.780     0.675
##   774   1251       2     0.779     0.674
##   775   1249       1     0.779     0.673
##   795   1248       1     0.778     0.672
##   797   1247       2     0.777     0.671
##   802   1245       2     0.776     0.670
##   803   1243       1     0.776     0.670
##   805   1242       1     0.776     0.669
##   806   1241       2     0.775     0.668
##   811   1239       1     0.774     0.667
##   827   1238       1     0.774     0.667
##   828   1237       1     0.774     0.666
##   832   1236       1     0.773     0.666
##   833   1235       2     0.772     0.665
##   835   1233       1     0.772     0.664
##   840   1232       1     0.772     0.663
##   844   1231       1     0.771     0.663
##   845   1229       1     0.771     0.662
##   846   1227       1     0.770     0.662
##   849   1226       1     0.770     0.661
##   851   1225       1     0.770     0.661
##   853   1224       1     0.769     0.660
##   854   1223       1     0.769     0.660
##   855   1222       1     0.768     0.659
##   858   1221       1     0.768     0.658
##   862   1220       1     0.767     0.658
##   863   1219       1     0.767     0.657
##   871   1218       1     0.767     0.657
##   874   1217       1     0.766     0.656
##   875   1216       1     0.766     0.656
##   883   1215       2     0.765     0.655
##   884   1213       1     0.765     0.654
##   885   1211       1     0.764     0.653
##   887   1210       3     0.763     0.652
##   890   1206       1     0.763     0.651
##   891   1205       1     0.762     0.651
##   900   1204       1     0.762     0.650
##   901   1203       1     0.761     0.650
##   902   1202       1     0.761     0.649
##   904   1201       1     0.760     0.648
##   905   1200       2     0.760     0.647
##   909   1198       1     0.759     0.647
##   911   1197       1     0.759     0.646
##   912   1196       1     0.758     0.646
##   916   1195       1     0.758     0.645
##   918   1194       1     0.758     0.644
##   922   1193       1     0.757     0.644
##   924   1192       1     0.757     0.643
##   928   1191       1     0.756     0.643
##   929   1190       1     0.756     0.642
##   930   1189       1     0.755     0.642
##   931   1188       1     0.755     0.641
##   934   1187       1     0.755     0.641
##   936   1186       2     0.754     0.639
##   938   1184       1     0.753     0.639
##   939   1182       1     0.753     0.638
##   940   1181       1     0.753     0.638
##   942   1180       1     0.752     0.637
##   944   1179       1     0.752     0.637
##   949   1178       1     0.751     0.636
##   952   1177       1     0.751     0.636
##   957   1176       1     0.750     0.635
##   959   1175       1     0.750     0.634
##   960   1174       1     0.750     0.634
##   961   1173       4     0.748     0.632
##   963   1169       1     0.747     0.631
##   966   1168       1     0.747     0.630
##   968   1167       2     0.746     0.629
##   969   1165       1     0.746     0.629
##   975   1164       1     0.745     0.628
##   976   1163       1     0.745     0.628
##   977   1162       1     0.745     0.627
##   986   1161       1     0.744     0.627
##   993   1160       1     0.744     0.626
##   997   1159       2     0.743     0.625
##  1013   1157       1     0.742     0.624
##  1018   1156       1     0.742     0.624
##  1020   1155       1     0.742     0.623
##  1021   1154       1     0.741     0.623
##  1022   1153       1     0.741     0.622
##  1024   1152       1     0.740     0.621
##  1025   1151       1     0.740     0.621
##  1026   1150       1     0.739     0.620
##  1029   1149       1     0.739     0.620
##  1031   1148       1     0.739     0.619
##  1032   1147       1     0.738     0.619
##  1034   1146       1     0.738     0.618
##  1037   1145       2     0.737     0.617
##  1041   1143       1     0.737     0.616
##  1042   1142       1     0.736     0.616
##  1046   1141       1     0.736     0.615
##  1048   1140       1     0.735     0.615
##  1052   1139       1     0.735     0.614
##  1055   1138       1     0.734     0.614
##  1057   1137       2     0.734     0.613
##  1061   1135       1     0.733     0.612
##  1070   1134       1     0.733     0.611
##  1079   1133       1     0.732     0.611
##  1081   1132       1     0.732     0.610
##  1083   1131       1     0.731     0.610
##  1089   1130       1     0.731     0.609
##  1092   1129       1     0.731     0.609
##  1101   1128       1     0.730     0.608
##  1103   1127       1     0.730     0.607
##  1105   1125       1     0.729     0.607
##  1106   1124       1     0.729     0.606
##  1108   1123       1     0.728     0.606
##  1112   1122       1     0.728     0.605
##  1114   1121       1     0.728     0.605
##  1117   1120       1     0.727     0.604
##  1122   1119       2     0.726     0.603
##  1130   1117       1     0.726     0.602
##  1133   1116       1     0.725     0.602
##  1134   1115       1     0.725     0.601
##  1135   1114       1     0.725     0.601
##  1136   1113       1     0.724     0.600
##  1138   1112       1     0.724     0.600
##  1139   1111       2     0.723     0.598
##  1142   1109       1     0.722     0.598
##  1145   1108       2     0.722     0.597
##  1151   1106       1     0.721     0.596
##  1154   1105       1     0.721     0.596
##  1159   1104       2     0.720     0.595
##  1161   1102       1     0.719     0.594
##  1166   1101       1     0.719     0.593
##  1178   1100       2     0.718     0.592
##  1183   1098       1     0.718     0.592
##  1186   1097       1     0.717     0.591
##  1191   1096       1     0.717     0.591
##  1193   1095       1     0.716     0.590
##  1195   1094       1     0.716     0.589
##  1198   1093       1     0.716     0.589
##  1201   1091       1     0.715     0.588
##  1207   1090       1     0.715     0.588
##  1209   1089       1     0.714     0.587
##  1211   1088       1     0.714     0.587
##  1212   1087       1     0.713     0.586
##  1215   1086       1     0.713     0.586
##  1216   1085       1     0.713     0.585
##  1219   1084       1     0.712     0.584
##  1230   1083       1     0.712     0.584
##  1233   1082       1     0.711     0.583
##  1236   1081       1     0.711     0.583
##  1237   1080       1     0.710     0.582
##  1246   1079       2     0.709     0.581
##  1252   1077       1     0.709     0.580
##  1262   1076       2     0.708     0.579
##  1272   1074       1     0.708     0.579
##  1273   1073       1     0.707     0.578
##  1274   1072       1     0.707     0.578
##  1275   1071       1     0.706     0.577
##  1276   1070       2     0.706     0.576
##  1277   1068       1     0.705     0.575
##  1279   1067       1     0.705     0.575
##  1290   1064       1     0.704     0.574
##  1295   1063       2     0.703     0.573
##  1298   1061       1     0.703     0.573
##  1302   1060       1     0.702     0.572
##  1304   1059       1     0.702     0.571
##  1306   1058       1     0.702     0.571
##  1313   1057       1     0.701     0.570
##  1314   1056       1     0.701     0.570
##  1323   1054       1     0.700     0.569
##  1325   1053       1     0.700     0.569
##  1327   1052       1     0.699     0.568
##  1329   1051       1     0.699     0.567
##  1353   1050       1     0.699     0.567
##  1363   1049       1     0.698     0.566
##  1365   1047       1     0.698     0.566
##  1375   1046       1     0.697     0.565
##  1387   1045       1     0.697     0.565
##  1388   1044       1     0.696     0.564
##  1399   1043       1     0.696     0.564
##  1405   1042       1     0.695     0.563
##  1424   1039       1     0.695     0.562
##  1432   1038       1     0.695     0.562
##  1434   1037       2     0.694     0.561
##  1436   1035       1     0.693     0.560
##  1437   1034       1     0.693     0.560
##  1439   1033       1     0.692     0.559
##  1446   1032       2     0.691     0.558
##  1447   1030       1     0.691     0.557
##  1455   1029       1     0.691     0.557
##  1466   1028       1     0.690     0.556
##  1471   1027       1     0.690     0.556
##  1475   1024       1     0.689     0.555
##  1482   1023       1     0.689     0.554
##  1488   1022       1     0.688     0.554
##  1495   1021       1     0.688     0.553
##  1509   1020       1     0.687     0.553
##  1511   1019       1     0.687     0.552
##  1521   1018       1     0.687     0.552
##  1530   1017       1     0.686     0.551
##  1535   1016       1     0.686     0.550
##  1539   1013       1     0.685     0.550
##  1540   1012       1     0.685     0.549
##  1548   1011       2     0.684     0.548
##  1550   1008       1     0.683     0.548
##  1551   1006       1     0.683     0.547
##  1561   1005       1     0.683     0.547
##  1564   1004       1     0.682     0.546
##  1568   1003       1     0.682     0.545
##  1589   1002       1     0.681     0.545
##  1606   1001       2     0.680     0.544
##  1607    999       1     0.680     0.543
##  1620    997       1     0.679     0.543
##  1637    996       1     0.679     0.542
##  1644    995       1     0.678     0.541
##  1647    994       1     0.678     0.541
##  1652    993       1     0.678     0.540
##  1656    992       1     0.677     0.540
##  1668    991       2     0.676     0.539
##  1671    989       1     0.676     0.538
##  1679    988       1     0.675     0.537
##  1687    987       1     0.675     0.537
##  1692    986       1     0.674     0.536
##  1709    985       1     0.674     0.536
##  1723    984       2     0.673     0.535
##  1743    982       1     0.673     0.534
##  1745    981       1     0.672     0.533
##  1749    980       1     0.672     0.533
##  1752    979       1     0.671     0.532
##  1759    978       1     0.671     0.532
##  1767    977       1     0.670     0.531
##  1768    976       1     0.670     0.531
##  1772    975       1     0.669     0.530
##  1783    974       1     0.669     0.529
##  1786    973       1     0.669     0.529
##  1788    972       1     0.668     0.528
##  1790    971       1     0.668     0.528
##  1798    968       1     0.667     0.527
##  1812    963       1     0.667     0.527
##  1818    958       1     0.666     0.526
##  1829    942       1     0.666     0.525
##  1831    941       1     0.665     0.525
##  1839    937       1     0.665     0.524
##  1850    936       1     0.664     0.524
##  1851    935       1     0.664     0.523
##  1856    928       1     0.663     0.522
##  1875    918       1     0.663     0.522
##  1876    915       1     0.662     0.521
##  1879    914       1     0.662     0.521
##  1884    911       1     0.661     0.520
##  1885    910       1     0.661     0.519
##  1895    907       1     0.661     0.519
##  1896    906       1     0.660     0.518
##  1907    898       1     0.660     0.518
##  1915    895       1     0.659     0.517
##  1918    894       1     0.659     0.516
##  1932    886       1     0.658     0.516
##  1950    880       1     0.658     0.515
##  1976    871       1     0.657     0.515
##  1981    864       1     0.657     0.514
##  1995    855       1     0.656     0.513
##  2012    840       1     0.656     0.513
##  2018    839       1     0.655     0.512
##  2021    834       1     0.655     0.511
##  2023    833       1     0.654     0.511
##  2028    828       1     0.653     0.510
##  2031    823       1     0.653     0.509
##  2035    822       1     0.652     0.509
##  2036    821       1     0.652     0.508
##  2052    814       1     0.651     0.508
##  2067    800       1     0.651     0.507
##  2074    795       1     0.650     0.506
##  2077    786       1     0.650     0.506
##  2079    785       1     0.649     0.505
##  2083    783       1     0.649     0.504
##  2085    780       1     0.648     0.504
##  2127    737       1     0.648     0.503
##  2128    736       1     0.647     0.502
##  2133    728       1     0.646     0.501
##  2148    711       1     0.646     0.501
##  2152    707       1     0.645     0.500
##  2171    671       2     0.644     0.498
##  2174    668       1     0.643     0.498
##  2197    618       1     0.643     0.497
##  2213    591       1     0.642     0.496
##  2231    571       1     0.641     0.495
##  2257    543       1     0.640     0.494
##  2284    513       1     0.639     0.493
##  2287    510       1     0.639     0.492
##  2288    509       1     0.638     0.491
##  2318    480       1     0.637     0.490
##  2351    446       1     0.636     0.489
##  2458    365       1     0.635     0.487
##  2482    346       1     0.634     0.486
##  2527    294       1     0.632     0.484
##  2542    276       1     0.631     0.482
##  2552    265       1     0.629     0.480
##  2593    235       1     0.627     0.478
##  2683    187       1     0.625     0.475
##  2695    181       1     0.623     0.473
##  2718    162       1     0.620     0.470
##  2725    154       1     0.617     0.466
##  2789    105       1     0.614     0.462
##  2910     59       1     0.606     0.453
```