`vignettes/Dynamic_Predictions.Rmd`

`Dynamic_Predictions.Rmd`

Based on the general framework of joint models presented earlier, we are interested in deriving cumulative risk probabilities for a new subject \(j\) that has survived up to time point \(t\) and has provided longitudinal measurements \(\mathcal Y_{kj}(t) = \{ y_{kj}(t_{jl}); 0 \leq t_{jl} \leq t, l = 1, \ldots, n_j, k = 1, \ldots, K\}\), with \(K\) denoting the number of longitudinal outcomes. The probabilities of interest are \[\begin{array}{l} \pi_j(u \mid t) = \mbox{Pr}\{T_j^* \leq u \mid T_j^* > t, \mathcal Y_j(t), \mathcal D_n\}\\\\ = \displaystyle 1 - \int\int \frac{S(u \mid b_j, \theta)}{S(t \mid b_j, \theta)} \; p\{b_j \mid T_j^* > t, \mathcal Y_j(t), \theta\} \; p(\theta \mid \mathcal D_n) \; db_j d\theta, \end{array}\] where \(S(\cdot)\) denotes the survival function conditional on the random effects, and \(\mathcal Y_j(t) = \{\mathcal Y_{1j}(t), \ldots, \mathcal Y_{Kj}(t)\}\). Combining the three terms in the integrand we can device a Monte Carlo scheme to obtain estimates of these probabilities, namely,

Sample a value \(\tilde \theta\) from the posterior of the parameters \([\theta \mid \mathcal D_n]\).

Sample a value \(\tilde b_j\) from the posterior of the random effects \([b_j \mid T_j^* > t, \mathcal Y_j(t), \tilde \theta]\).

Compute the ratio of survival probabilities \(S(u \mid \tilde b_j, \tilde \theta) \Big / S(t \mid \tilde b_j, \tilde \theta)\).

Replicating these steps \(L\) times, we can estimate the conditional cumulative risk probabilities by \[1 - \frac{1}{L} \sum_{l=1}^L \frac{S(u \mid \tilde b_j^{(l)}, \tilde \theta^{(l)})}{S(t \mid \tilde b_j^{(l)}, \tilde \theta^{(l)})},\] and their standard error by calculating the standard deviation across the Monte Carlo samples.

We will illustrate the calculation of dynamic predictions using package **JMbayes2** from a trivariate joint model fitted to the PBC dataset for the longitudinal outcomes `serBilir`

(continuous), `prothrombin`

time (continuous) and `ascites`

(dichotomous). We start by fitting the univariate mixed models. For the two continuous outcomes, we allow for nonlinear subject-specific time effects using natural cubic splines. For `ascites`

, we postulate linear subject-specific profiles for the log odds. The code is:

```
fm1 <- lme(log(serBilir) ~ ns(year, 3) * sex, data = pbc2,
random = ~ ns(year, 3) | id, control = lmeControl(opt = 'optim'))
fm2 <- lme(prothrombin ~ ns(year, 2) * sex, data = pbc2,
random = ~ ns(year, 2) | id, control = lmeControl(opt = 'optim'))
fm3 <- mixed_model(ascites ~ year * sex, data = pbc2,
random = ~ year | id, family = binomial())
```

Following, we fit the Cox model for the time to either transplantation or death. The first line defines the composite event indicator, and the second one fits the Cox model in which we have also included the baseline covariates `drug`

and `age`

. The code is:

```
pbc2.id$event <- as.numeric(pbc2.id$status != "alive")
CoxFit <- coxph(Surv(years, event) ~ drug + age, data = pbc2.id)
```

The joint model is fitted with the following call to `jm()`

:

We want to calculate predictions for the longitudinal and survival outcomes for Patients 25 and 93. As a first step, we extract the data of these patients and store them in the data.frame `ND`

with the code:

```
t0 <- 5
ND <- pbc2[pbc2$id %in% c(25, 93), ]
ND <- ND[ND$year < t0, ]
ND$status2 <- 0
ND$years <- t0
```

We will only use the first five years of follow-up (line three), and further we specify that the patients were event-free up to this time point (lines four and five).

We start with predictions for the longitudinal outcomes. These are produced by the `predict()`

method for class `jm`

objects, and follow the same lines as the procedure described above for cumulative risk probabilities. The only difference is in Step 3, where instead of calculating the cumulative risk we calculate the predicted values for the longitudinal outcomes. There are two options controlled by the `type_pred`

argument, namely predictions at the scale of the response/outcome (default) or at the linear predictor level. The `type`

argument controls if the predictions will be for the mean subject (i.e., including only the fixed effects) or subject-specific including both the fixed and random effects. In the `newdata`

argument we provide the available measurements of the two patients. This will be used to sample their random effects at Step 2 presented above. This is done with a Metropolis-Hastings algorithm that runs for `n_mcmc`

iterations; all iterations but the last one are discarded as burn-in. Finally, argument `n_samples`

corresponds to the value of \(L\) defined above and specifies the number of Monte Carlo samples:

`predLong1 <- predict(jointFit, newdata = ND, return_newdata = TRUE)`

Argument `return_newdata`

specifies that the predictions are returned as extra columns of the `newdata`

data.frame. By default the 95% credible intervals are also included. Using the `plot()`

method for objects returned by `predict.jm(..., return_newdata = TRUE)`

, we can display the predictions. With the following code we do that for the first longitudinal outcome:

`plot(predLong1)`

When we want to calculate predictions for other, future time points, we can accordingly specify the `times`

argument. In the following example, we calculate predictions from time `t0`

to time 12:

```
predLong2 <- predict(jointFit, newdata = ND,
times = seq(t0, 12, length.out = 51),
return_newdata = TRUE)
```

We show these predictions for the second outcome and the second patient (i.e., Patient 93). This is achieved by suitably specifying the `outcomes`

and `subject`

arguments of the `plot()`

method:

`plot(predLong2, outcomes = 2, subject = 93)`

We continue with the predictions for the event outcome. To let `predict()`

know that we want the cumulative risk probabilities, we specify `process = "event"`

:

```
predSurv <- predict(jointFit, newdata = ND, process = "event",
times = seq(t0, 12, length.out = 51),
return_newdata = TRUE)
```

The predictions are included again as extra columns in the corresponding data.frame. To depict the predictions of both the longitudinal and survival outcomes combined, we provide both objects to the `plot()`

method:

`plot(predLong2, predSurv)`

Again by default, the plot is for the predictions of the first subject (i.e., Patient 25) and for the first longitudinal outcome (i.e., `log(serBilir)`

). However, the `plot()`

method has a series of arguments that allows users to customize the plot. We illustrate some of these capabilities with the following figure. First, we specify that we want to depict all three outcomes using `outcomes = 1:3`

(note: a max of three outcomes can be simultaneously displayed). Next, we specify via the `subject`

argument that we want to show the predictions of Patient 93. Note, that for serum bilirubin we used the log transformation in the specification of the linear mixed model. Hence, we receive predictions on the transformed scale. To show predictions on the original scale, we use the `fun_long`

argument. Because we have three outcomes, this needs to be a list of three functions. The first one, corresponding to serum bilirubin is the `exp()`

and for the other two the `identity()`

because we do not wish to transform the predictions. Analogously, we also have the `fun_event`

argument to transform the predictions for the event outcome, and in the example below we set that we want to obtain survival probabilities. Using the arguments `bg`

, `col_points`

, `col_line_long`

, `col_line_event`

, `fill_CI_long`

, and `fill_CI_event`

we have changed the appearance of the plot to a dark theme. Finally, the `pos_ylab_long`

specifies the relative positive of the y-axis labels for the three longitudinal outcomes.

```
cols <- c('#F25C78', '#D973B5', '#F28322')
plot(predLong2, predSurv, outcomes = 1:3, subject = 93,
fun_long = list(exp, identity, identity),
fun_event = function (x) 1 - x,
ylab_event = "Survival Probabilities",
ylab_long = c("Serum Bilirubin", "Prothrombin", "Ascites"),
bg = '#132743', col_points = cols, col_line_long = cols,
col_line_event = '#F7F7FF', col_axis = "white",
fill_CI_long = c("#F25C7880", "#D973B580", "#F2832280"),
fill_CI_event = "#F7F7FF80",
pos_ylab_long = c(1.9, 1.9, 0.08))
```

We evaluate the discriminative capability of the model using ROC methodology. We calculate the components of the ROC curve using information up to year five, and we are interested in events occurring within a three-year window. That is discriminating between patients who will get the event in the interval `(t0, t0 + Dt]`

, (i.e., in our case \(T_j \in (5, 8]\)) from patients who will survive at least 8 years (i.e., \(T_j > 8\)). The calculations are performed with the following call to `tvROC()`

:

```
pbc2$event <- as.numeric(pbc2$status != "alive")
roc <- tvROC(jointFit, newdata = pbc2, Tstart = t0, Dt = 3)
roc
#>
#> Time-dependent Sensitivity and Specificity for the Joint Model jointFit
#>
#> At time: 8
#> Using information up to time: 5 (202 subjects still at risk)
#>
#> cut-off SN SP qSN qSP
#> 1 0.06 0.02078117 1.000000000 0.01590944 1.000000000
#> 2 0.08 0.03349069 0.997475870 0.02382560 0.745062856
#> 3 0.11 0.03349069 0.990977279 0.01892020 0.392470137
#> 4 0.12 0.07505303 0.990977279 0.05157722 0.635482082
#> 5 0.13 0.09284276 0.990041808 0.06507264 0.664740858
#> 6 0.14 0.09284276 0.983543217 0.06027814 0.525104264
#> 7 0.16 0.11362393 0.977044626 0.07229033 0.484779394
#> 8 0.18 0.13440510 0.977044626 0.08932204 0.536301455
#> 9 0.19 0.17596745 0.977044626 0.12392329 0.613584546
#> 10 0.21 0.19674862 0.977044626 0.14149853 0.643308811
#> 11 0.22 0.21752979 0.977044626 0.15926074 0.668786753
#> 12 0.24 0.23831097 0.977044626 0.17721292 0.690867636
#> 13 0.26 0.25909214 0.977044626 0.19535813 0.710188409
#> 14 0.30 0.27987331 0.977044626 0.21369951 0.727236150
#> 15 0.32 0.30065448 0.970546036 0.22804484 0.686857855
#> 16 0.33 0.33383151 0.967923817 0.25654124 0.691457671
#> 17 0.34 0.35461268 0.961425227 0.27168582 0.661212860
#> 18 0.36 0.37539385 0.961425227 0.29117729 0.675328991
#> 19 0.37 0.39617502 0.961425227 0.31088901 0.688315831
#> 20 0.41 0.41695620 0.961425227 0.33082472 0.700303684
#> 21 0.42 0.42133847 0.956297041 0.33205926 0.673037119
#> 22 0.45 0.45299606 0.946701080 0.35758839 0.641120604
#> 23 0.49 0.45299606 0.933703899 0.35003061 0.581505860
#> 24 0.52 0.45832648 0.928872212 0.35255591 0.564611722
#> 25 0.55 0.45832648 0.922373621 0.34870207 0.538807986
#> 26 0.56 0.47910765 0.922373621 0.36993859 0.551984901
#> 27 0.57 0.47910765 0.915875031 0.36614304 0.527965451
#> 28 0.58 0.50830685 0.912008888 0.39437796 0.532257775
#> 29 0.59 0.51481328 0.901046367 0.39501409 0.500284152
#> 30 0.60 0.53949478 0.895767469 0.41861217 0.498691159
#> 31 0.62 0.58105713 0.895767469 0.46438949 0.521477925
#> 32 0.64 0.58105713 0.882770288 0.45752269 0.485208654
#> 33 0.65 0.60183830 0.882770288 0.48110540 0.496161662
#> 34 0.66 0.60183830 0.876271697 0.47773595 0.479310058
#> 35 0.68 0.62261947 0.863274516 0.49515983 0.458459847
#> 36 0.70 0.64340065 0.856775925 0.51655658 0.454127111
#> 37 0.71 0.64340065 0.843778744 0.50997912 0.426241932
#> 38 0.72 0.66418182 0.824282972 0.52562747 0.398392549
#> 39 0.73 0.68832103 0.812335901 0.55029177 0.388578257 *
#> 40 0.74 0.69240384 0.807114069 0.55298976 0.381540347
#> 41 0.76 0.69756538 0.802229569 0.55730584 0.375787076
#> 42 0.77 0.71011565 0.773661277 0.55972451 0.337385186
#> 43 0.78 0.71011565 0.760664095 0.55300275 0.319072497
#> 44 0.79 0.71447043 0.742530130 0.54939864 0.297176843
#> 45 0.80 0.71920892 0.737513338 0.55338741 0.293035923
#> 46 0.81 0.76338342 0.718834427 0.60822500 0.290056928
#> 47 0.82 0.76960102 0.694784406 0.60558820 0.266029166
#> 48 0.83 0.79114076 0.682024440 0.63313421 0.261711919
#> 49 0.84 0.82209010 0.639714047 0.66413272 0.233918289
#> 50 0.85 0.85156712 0.609940418 0.70313424 0.219880835
#> 51 0.86 0.85156712 0.583946055 0.69089235 0.199591459
#> 52 0.87 0.85615472 0.565884896 0.69088567 0.188043973
#> 53 0.88 0.85852204 0.547129423 0.68595003 0.175857148
#> 54 0.89 0.85927969 0.495377624 0.65752405 0.143414118
#> 55 0.90 0.88006086 0.456386081 0.68121439 0.128491970
#> 56 0.91 0.88352483 0.424976363 0.66861993 0.113322331
#> 57 0.92 0.90903674 0.380965605 0.70834003 0.100396059
#> 58 0.93 0.93078132 0.335776742 0.74577866 0.087257836
#> 59 0.94 0.93238287 0.258294483 0.68235672 0.057707457
#> 60 0.95 0.93432300 0.174419514 0.55777484 0.030422917
#> 61 0.96 0.97815379 0.090647216 0.70580441 0.017704479
#> 62 0.97 1.00000000 0.025994362 1.00000000 0.006317481
#> 63 0.98 1.00000000 0.006498591 1.00000000 0.001555798
```

In the first line we define the event indicator as we did in the `pbc2.id`

data.frame. The cut-point with the asterisk on the right maximizes the Youden’s index. To depict the ROC curve, we use the corresponding `plot()`

method:

The area under the ROC curve is calculated with the `tvAUC()`

function:

```
tvAUC(roc)
#>
#> Time-dependent AUC for the Joint Model jointFit
#>
#> Estimated AUC: 0.8055
#> At time: 8
#> Using information up to time: 5 (202 subjects still at risk)
```

This function either accepts an object of class `tvROC`

or of class `jm`

. In the latter case, the user must also provide the `newdata`

, `Tstart`

and `Dt`

or `Thoriz`

arguments. Here we have used the same dataset as the one to fit the model, but, in principle, discrimination could be (better) assessed in another dataset.

To assess the accuracy of the predictions we produce a calibration plot:

`calibration_plot(jointFit, newdata = pbc2, Tstart = t0, Dt = 3)`

The syntax of the `calibration_plot()`

function is almost identical to that of `tvROC()`

. The kernel density estimation is of the estimated probabilities \(\pi_j(t + \Delta t \mid t) = \pi_j(8 \mid 5)\) for all individuals at risk at year `t0`

in the data frame provided in the `newdata`

argument. Using the `calibration_metrics()`

function we can also calculate metrics for the accuracy of predictions:

```
calibration_metrics(jointFit, pbc2, Tstart = 5, Dt = 3)
#> ICI E50 E90
#> 0.04380478 0.02027983 0.11170119
```

The ICI is the mean absolute difference between the observed and predicted probabilities, E50 is the median absolute difference, and E90 is the 90% percentile of the absolute differences. Finally, we calculate the Brier score as an overall measure of predictive performance. This is computed with the `tvBrier()`

function:

```
tvBrier(jointFit, newdata = pbc2, Tstart = t0, Dt = 3)
#>
#> Prediction Error for the Joint Model jointFit
#>
#> Estimated Brier score: 0.1328
#> At time: 8
#> Using information up to time: 5 (202 subjects still at risk)
```

**Notes:**

- To obtain valid estimates of the predictive accuracy measures (i.e., time-varying sensitivity, specificity, and Brier score) we need to account for censoring. A popular method to achieve this is via inverse probability of censoring weighting. For this approach to be valid, we need the model for the weights to be correctly specified. In standard survival analysis, this is achieved either using the Kaplan-Meier estimator or a Cox model for the censoring distribution. However, in the settings where joint models are used, it is often the case that the censoring mechanism may depend on the history of the longitudinal outcomes in a complex manner. This is especially the case when we consider multiple longitudinal outcomes in the analysis. Also, these outcomes may be recorded at different time points per patient and have missing data. Because of these reasons, in these settings, Kaplan-Meier-based or Cox-based censoring weights may be difficult to derive or be biased. The functions in
**JMbayes2**that calculate the predictive accuracy measures use joint-model-based weights to account for censoring. These weights allow censoring to depend in any possible manner on the history of the longitudinal outcomes. However, they require that the model is appropriately calibrated. - The calibration curve, produced by
`calibration_plot()`

, and the calibration metrics, produced by`calibration_metrics())`

, are calculated using the procedure described in Austin et al., 2020.