---
title: "Superlandmarking with static covariates"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Superlandmarking with static covariates}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup, message = FALSE}
library(landmaRk)
library(survival)
```

# Overview

This vignette fits both a basic landmarking model and a superlandmarking
model to a cohort of ovarian cancer patients, the dataset `ova`
(see `?ova`). The example is taken from Section 7.2 ("Obtaining dynamic
predictions by landmarking") of van Houwelingen HC, Putter H (2012),
*Dynamic Prediction in Clinical Survival Analysis*.

Landmarking produces a prediction from landmark time $s$ up to horizon
$s + w$ by fitting a survival (typically Cox) model to the risk set still at risk at $s$.
Optionally, the dataset is administratively censored at $s + w$.
$$
h(t \mid x, s, w) = h_0(t \mid s, w) \exp(x^\top \beta_{LM}), \qquad
s \le t \le s + w.
$$
Repeating this at a grid of landmark times $s_1, \dots, s_L$ gives landmark-specific
estimate of $\beta_{LM}$ per landmark (the "Crude" estimates below).

Superlandmarking instead pools every landmark-specific dataset into
a single stacked superdataset, and lets $\beta_{LM}$ vary smoothly with $s$
via a basis expansion rather than being re-estimated independently at each
landmark:
$$
h(t \mid x, s, w) = h_0(t \mid s, w) \exp\{x^\top \beta_{LM}(s)\}, \qquad
s \le t \le s + w, 
$$
$$
\beta_{LM}(s) = \sum_{j=1}^{m_b} \gamma_j f_j(s).
$$
Fitting this is a single Cox model stratified on $s$ (giving every landmark
its own baseline hazard $h_0(t \mid s, w)$) with covariate $\times f_j(s)$
interaction terms for $\hat\gamma$ — exactly what the `formula` argument of
`fit_superlandmark()` specifies directly, via `strata(LM)` and
`covariate:f_j(LM)` terms.

The stratified baseline $h_0(t \mid s, w)$ is only defined at the
landmark times actually used for fitting. Instead, the following equation lets a single
baseline hazard vary smoothly with $s$ too:
$$
h_0(t \mid s, w) = h_0(t) \exp\{\theta(s)\}, \qquad
\theta(s) = \sum_{j=1}^{m_h} \eta_j g_j(s), 
$$
for basis functions $g_j$ standardized so that $g_j(s_1) = 0$. 

The outcome of interest is overall survival. Following Section 7.2 of van Houwelingen and Putter,
the predictor is a single
prognostic index summarizing a patient's clinical risk factors at baseline, and the prediction window is
$w = 3$ years following each landmark time. The 
superlandmarking approach is used because the proportional hazard assumption
does not hold here.

# Example: ovarian cancer dataset

```{r str-ova}
data(ova)
str(ova)
```

`ova` contains 358 patients from two clinical trials of combination
chemotherapy for advanced ovarian cancer (Neijt et al. 1984, 1987). Every
column other than the outcome and patient identifier is a baseline covariate
recorded once, at diagnosis:

- `tyears`: time in years from diagnosis to death or last follow-up (the
  survival outcome's time).
- `d`: event indicator, `1` = death, `0` = censored (the survival outcome's
  status).
- `Karn`: Karnofsky performance score at diagnosis, recoded from 6
  ($\le 60$) to 10 (100); higher is better.
- `Broders`: Broders' histological tumour grade, a factor with levels
  `"1"`-`"4"` and `"unknown"`.
- `FIGO`: FIGO disease stage at diagnosis, a factor with levels `"III"` and
  `"IV"`.
- `Ascites`: presence of ascites at diagnosis, a factor with levels
  `"absent"`, `"present"`, `"unknown"`.
- `Diam`: diameter of the residual tumour after surgery, a factor with
  levels `"micr."`, `"<1cm"`, `"1-2cm"`, `"2-5cm"`, `">5cm"`.
- `id`: patient identifier.

There is no time-varying (longitudinal) covariate in this data set. 

# Static Cox model

```{r cox-static}
fit_baseline <- coxph(
  Surv(tyears, d) ~ FIGO + Diam + Broders + Ascites + Karn,
  data = ova
)
coef(summary(fit_baseline))
```

The fitted (centered)
linear predictor of this static Cox model (prognostic index) is used to inform the landmarking and superlandmarking models.

```{r pi}
lp <- predict(fit_baseline, type = "lp")
ova$PI <- lp - mean(lp)
c(mean = mean(ova$PI), sd = sd(ova$PI))
```

# Setting up the landmark analysis

`PI` is a single, time-fixed covariate (computed once from baseline data), so
this analysis does not require a longitudinal submodel. `LandmarkAnalysis()`
still requires `data_dynamic` to be a named list of data frames, so an empty
placeholder is supplied.

```{r setup-object}
x <- LandmarkAnalysis(
  data_static = ova,
  data_dynamic = list(
    none = data.frame(id = integer(0), time = numeric(0), value = numeric(0))
  ),
  event_indicator = "d",
  ids = "id",
  event_time = "tyears",
  times = "time",
  measurements = "value"
)
```

Following van Houwelingen and Putter (2012), landmark times are chosen on an equidistant grid from
$s_1 = 0$ to $s_L = 5$ years with distance 0.1, with a fixed prediction
window $w = 3$ years.

```{r risk-sets}
sL <- 5
w <- 3
landmarks <- seq(0, sL, by = 0.1)
x <- compute_risk_sets(x, landmarks)
```

# Crude landmarking

A separate Cox model is fitted at each landmark time, using `fit_survival()`.

```{r crude}
x <- fit_survival(
  x,
  formula = Surv(event_time, event_status) ~ PI,
  landmarks = landmarks,
  horizons = landmarks + w,
  method = "coxph"
)

# One fit per landmark, e.g. at s = 0:
coef(summary(x@survival_fits[["0-3"]]))
```

# Superlandmarking

A single Cox supermodel is fitted using `fit_superlandmark()` to the stacked superdataset, allowing the effect of `PI` to
vary smoothly with landmark time via the basis
$f_1(s) \equiv 1, f_2(s) = s/s_L, f_3(s) = (s/s_L)^2$ used in van Houwelingen
and Putter (2012) (Equation 7.2), with `strata(LM)` for the landmark-specific
baseline hazard.

```{r supermodel}
x <- fit_superlandmark(
  x,
  formula = Surv(event_time, event_status) ~
    PI + PI:I(LM / sL) + PI:I((LM / sL)^2) + strata(LM),
  landmarks = landmarks,
  w = w
)

fit_super <- x@superlandmark_fits[[length(x@superlandmark_fits)]]
summary(fit_super)$coefficients
```

# Displaying results

Hazard ratios estimated using both the crude landmarking and superlandmarking
approaches can be displayed using `plot_hazard_ratio()`.

```{r figure-7-9, fig.width = 6, fig.height = 4}
plot_hazard_ratio(x, covariate = "PI", landmarks = landmarks, w = w, type = "both") +
  ggplot2::coord_cartesian(ylim = c(-1, 1.5))
```