Rizopoulos 2016 - Dynamic Predictions and Bayesian Model Averaging

Summary

For a new patient with biomarker history $Y_j(t)$ and survival to $t$, the model gives the conditional survival probability ${\pi }_j(u\mid t)=P(T_j^{\ast }\ge u\mid T_j^{\ast }>t,Y_j(t),{\mathcal{D}}_n)$ and expected future biomarker ${\omega }_j(u\mid t)$, updated as new measurements arrive. Predictions from several association structures are combined with subject- and time-specific posterior model weights, computed by Laplace approximation.

Main Content

Dynamic predictions (Sec. 5.1)

$${\pi }_j(u\mid t)=\int P(T_j^{\ast }\ge u\mid T_j^{\ast }>t,Y_j(t),\theta )p(\theta \mid {\mathcal{D}}_n)d\theta ,{\omega }_j(u\mid t)=\int \mathbb{E}[y_j(u)\mid T_j^{\ast }>t,Y_j(t),\theta ]p(\theta \mid {\mathcal{D}}_n)d\theta .$$

Conditional survival via random effects

$$P(T_j^{\ast }\ge u\mid T_j^{\ast }>t,Y_j(t),\theta )=\int \frac{S_j(u\mid H_j(u,{\mathbf{b}}_j),\theta )}{S_j(t\mid H_j(t,{\mathbf{b}}_j),\theta )}p({\mathbf{b}}_j\mid T_j^{\ast }>t,Y_j(t),\theta )d{\mathbf{b}}_j,$$

and $\mathbb{E}[y_j(u)\mid \dots ]={\mathbf{x}}_j(u{)}^{\top }\mathbf{\beta }+{\mathbf{z}}_j(u{)}^{\top }{\bar{\mathbf{b}}}_j^{(t)}$ with ${\bar{\mathbf{b}}}_j^{(t)}=\int {\mathbf{b}}_{j}p({\mathbf{b}}_j\mid T_j^{\ast }>t,Y_j(t),\theta )d{\mathbf{b}}_j$. Estimated by Monte Carlo over posterior draws of $\theta$ and ${\mathbf{b}}_j$ (200 samples by default; survfitJM(), predict()); a first-order version ${\tilde{\pi }}_j(u\mid t)=S_j(u\mid H_j(u,{\hat{\mathbf{b}}}_j),\hat{\theta })/S_j(t\mid H_j(t,{\hat{\mathbf{b}}}_j),\hat{\theta })$ is available with simulate=FALSE. Conditional independence (Eqs. 4-5) makes the first factor tractable.

Patient 2 example. Using jointFit.pbc15 and last measurement at about 8.83 years, the event-free probability falls from 1.00 to 0.99 (8.92 y), 0.945 (9.25 y), 0.90 (9.58 y), 0.857 (9.90 y), 0.81 (10.23 y), 0.77 (10.55 y), 0.72 (10.88 y), 0.68 (11.21 y) (posterior mean, widening 95% intervals: at 11.21 y the interval is 0.27-0.93). predict(type="Subject", interval=...) gives future log bilirubin with intervals; a shiny app (runDynPred()) provides the same interactively.

Bayesian model averaging (Eq. 10)

For models $M_1..M_K$ (different association structures):

$$P(T_j^{\ast }>u\mid {\mathcal{D}}_j(t),{\mathcal{D}}_n)=\sum _{k=1}^{K}P(T_j^{\ast }>u\mid M_k,{\mathcal{D}}_j(t),{\mathcal{D}}_n)p(M_k\mid {\mathcal{D}}_j(t),{\mathcal{D}}_n),$$

with weights $p(M_k\mid {\mathcal{D}}_j(t),{\mathcal{D}}_n)=\frac{p({\mathcal{D}}_j(t)\mid M_k)p({\mathcal{D}}_n\mid M_k)p(M_k)}{{\sum }_{\ell }p({\mathcal{D}}_j(t)\mid M_{\ell })p({\mathcal{D}}_n\mid M_{\ell })p(M_{\ell })}$, where ${\mathcal{D}}_j(t)=\{T_j^{\ast }>t,Y_j(t),{\mathbf{w}}_j\}$. Marginal densities use Laplace approximations (first integrate random effects, then parameters; Tierney and Kadane 1986). The new patient informs only the weights, not $\theta$; priors $p(M_k)=1/K$. Weights are subject- and time-specific (Rizopoulos et al. 2014). In JMbayes: marglogLik(), logLik(marginal.thetas=TRUE), bma.combine().

BMA for Patient 2 (first five measurements)

Models pbc1, pbc12, pbc13, pbc14, pbc15. The weights are dominated by pbc14 (shared random effects). BMA event-free probabilities fall from 0.99 at 5.0 y to 0.74 at 7.9 y and 0.43 at 11.2 y (95% bands 0.09-0.75 at the end). Interpretation: with only five measurements the model that captures baseline level and slope through random effects explains this patient’s trajectory best; a different patient or later time could favour other structures.

import numpy as np
def bma_weights(log_pDj, log_pDn, log_prior):
    w = np.asarray(log_pDj) + np.asarray(log_pDn) + np.asarray(log_prior)
    w = np.exp(w - w.max()); return w / w.sum()       # stabilised softmax, as in the paper's R code

Connections

Model Comparison, Survival Analysis, Rizopoulos 2016 - Predictive Accuracy and Cross-Validation.

See Also

Rizopoulos 2016 - JMbayes Joint Models Overview