Intro

The articles Flexible Bayesian Additive Joint Models with an Application to Type 1 Diabetes Research (Köhler et al. 2017) and Nonlinear Association Structures in Flexible Bayesian Additive Joint Models (Köhler, Umlauf, and Greven 2018) present a flexible framework for estimating joint models for longitudinal and survival data using MCMC. This vignette provides some example code how to set up joint models using bamlss.

This vignette explains how to estimate joint models with the following hazard rate model: \[ h_{i}(t) = \exp\left[\eta_{i}(t)\right] = \exp [\eta_{\lambda i}(t) + \eta_{\gamma i} + \eta_{\alpha i}(\eta_{\mu i}(t), t) ] \] The longitudinal response \(y_{ij}\) at time points \(t_{ij}\) is modeled as \[ y_{ij}=\eta_{\mu i}(t_{ij}) + \varepsilon_{ij} \]
with independent normal errors \(\varepsilon_{ij} \sim N(0, \exp[\eta_{\sigma_i}(t_{ij})]^2)\).

In this model the predictors \(\eta_k = \left( \eta^{\lambda}, \eta^{\gamma}, \eta^{\alpha}, \eta^{\mu}, \eta^{\sigma} \right)\) are a composition of \[ \eta_k=\eta_k\left(\mathbf{x}; \boldsymbol{\beta}_k\right) = f_k(\mathbf{x}) = f_{1k}(\mathbf{x}_{1k})+\ldots+f_{J_kk}(\mathbf{x}_{J_kk}), \] where functions \(f_{jk}( \cdot )\), \(j = 1, \ldots, J_k\) are unspecified functions of generic covariates \(\mathbf{x}\), e.g., that are modeled by splines. There are two types of parameters in this model, time-varying (\(\lambda\), \(\alpha\) and \(\mu\)) and time-constant (\(\gamma\) and \(\sigma\)) parameters. Note that parameter \(\lambda\) includes the baseline hazard effect.

For more details on model specifications and methodology see Köhler et al. (2017) and Köhler, Umlauf, and Greven (2018).

Mayo Clinic Primary Biliary Cirrhosis Data

This example is taken from the JMbayes package (Rizopoulos 2016). The serum bilirubin level is one of the most important body markers associated with the PBC liver disease. Here, we try to investigate how the mean bilirubin level influences survival. Therefore, we decompose the time-varying effect in parameter \(\mu\) into a main effect and a functional random intercept, to account for individual deviations.

The data can be loaded with

data("pbc2", package = "JMbayes")
head(pbc2)
##   id    years status      drug      age    sex      year ascites
## 1  1  1.09517   dead D-penicil 58.76684 female 0.0000000     Yes
## 2  1  1.09517   dead D-penicil 58.76684 female 0.5256817     Yes
## 3  2 14.15234  alive D-penicil 56.44782 female 0.0000000      No
## 4  2 14.15234  alive D-penicil 56.44782 female 0.4983025      No
## 5  2 14.15234  alive D-penicil 56.44782 female 0.9993429      No
## 6  2 14.15234  alive D-penicil 56.44782 female 2.1027270      No
##   hepatomegaly spiders                   edema serBilir serChol albumin
## 1          Yes     Yes edema despite diuretics     14.5     261    2.60
## 2          Yes     Yes edema despite diuretics     21.3      NA    2.94
## 3          Yes     Yes                No edema      1.1     302    4.14
## 4          Yes     Yes                No edema      0.8      NA    3.60
## 5          Yes     Yes                No edema      1.0      NA    3.55
## 6          Yes     Yes                No edema      1.9      NA    3.92
##   alkaline  SGOT platelets prothrombin histologic status2
## 1     1718 138.0       190        12.2          4       1
## 2     1612   6.2       183        11.2          4       1
## 3     7395 113.5       221        10.6          3       0
## 4     2107 139.5       188        11.0          3       0
## 5     1711 144.2       161        11.6          3       0
## 6     1365 144.2       122        10.6          3       0

and contains several potentially relevant covariates explaining the disease. Note that for the survival part, the time variable is years and for the longitudinal part year in this data set. To set up the joint model we need to specify formulas for each parameter (\(\lambda\), \(\gamma\), \(\mu\), \(\sigma\), \(\alpha\)) of the model.

## Set up the model formula including
## functional random intercepts using ti().
f <- list(
  Surv2(years, status2, obs = log(serBilir)) ~ s(years,k=20),
  gamma ~ s(age,k=20) + drug + sex,
  mu ~ ti(id,bs="re") + 
       ti(year,k=20) + 
       ti(id,year,bs=c("re","cr"),k=c(nlevels(pbc2$id),8)),
  sigma ~ 1,
  alpha ~ s(years,k=20),
  dalpha ~ -1
)

The model is estimated by first running a backfitting algorithm, to find appropriate starting values, afterwards the MCMC sampler is started.

## Set the seed for reproducibility.
set.seed(123)

## Estimate model.
b <- bamlss(f, data = pbc2, family = "jm",
  timevar = "year", idvar = "id")

The model summary gives

## 
## Call:
## bamlss(formula = f, family = "jm", data = pbc2, timevar = "year", 
##     idvar = "id")
## ---
## Family: jm 
## Link function: lambda = log, gamma = log, mu = identity, sigma = log, alpha = identity, dalpha = identity
## *---
## Formula lambda:
## ---
## Surv2(years, status2, obs = log(serBilir)) ~ s(years, k = 20)
## -
## Smooth terms:
##                     Mean      2.5%       50%     97.5% parameters
## s(years).tau21 5.060e+00 2.095e-04 2.770e-01 4.387e+01      0.000
## s(years).edf   2.395e+00 1.001e+00 1.769e+00 6.330e+00      0.103
## s(years).alpha 8.921e-01 2.108e-01 9.867e-01 1.000e+00         NA
## ---
## Formula gamma:
## ---
## gamma ~ s(age, k = 20) + drug + sex
## -
## Parametric coefficients:
##                   Mean     2.5%      50%    97.5% parameters
## (Intercept)   -4.33809 -5.19152 -4.34750 -3.54898     -2.134
## drugD-penicil -0.13823 -0.49208 -0.13589  0.19728     -0.081
## sexfemale     -0.03222 -0.54328 -0.03723  0.48846     -0.568
## -
## Acceptance probabilty:
##         Mean   2.5%    50% 97.5%
## alpha 0.8499 0.0771 0.9707     1
## -
## Smooth terms:
##                   Mean      2.5%       50%     97.5% parameters
## s(age).tau21 0.6067465 0.0001197 0.0147995 8.8790071      0.000
## s(age).edf   1.5123370 1.0006645 1.1297232 4.9789590      0.608
## s(age).alpha 0.9179284 0.2268342 0.9931830 1.0000000         NA
## ---
## Formula mu:
## ---
## mu ~ ti(id, bs = "re") + ti(year, k = 20) + ti(id, year, bs = c("re", 
##     "cr"), k = c(nlevels(pbc2$id), 8))
## -
## Parametric coefficients:
##               Mean   2.5%    50%  97.5% parameters
## (Intercept) 0.7109 0.5619 0.6258 1.0593      0.603
## -
## Acceptance probabilty:
##         Mean   2.5%    50% 97.5%
## alpha 0.9978 0.9850 1.0000     1
## -
## Smooth terms:
##                         Mean       2.5%        50%      97.5% parameters
## ti(id).tau21       4.027e-01  7.969e-02  7.969e-02  1.376e+00      0.000
## ti(id).edf        -2.253e+01 -2.670e+01 -2.265e+01 -1.382e+01      0.000
## ti(id).alpha       6.804e-02  6.040e-24  1.482e-05  1.000e+00         NA
## ti(year).tau21     1.246e-04  2.762e-05  9.255e-05  4.187e-04      0.000
## ti(year).edf      -7.169e-01 -1.356e+00 -6.795e-01 -2.440e-01      0.041
## ti(year).alpha     9.863e-01  8.844e-01  9.996e-01  1.000e+00         NA
## ti(id,year).tau21  1.821e+02  3.785e+01  2.160e+02  2.910e+02      0.000
## ti(id,year).tau22  9.708e-05  3.451e-05  7.889e-05  2.038e-04      0.000
## ti(id,year).edf   -2.121e+01 -3.244e+01 -2.400e+01 -8.790e+00      0.000
## ti(id,year).alpha  2.698e-01  1.115e-18  3.000e-02  1.000e+00         NA
## ---
## Formula sigma:
## ---
## sigma ~ 1
## -
## Parametric coefficients:
##                Mean    2.5%     50%   97.5% parameters
## (Intercept) -0.3689 -0.4954 -0.3296 -0.2375      0.104
## -
## Acceptance probabilty:
##             Mean       2.5%        50% 97.5%
## alpha  6.962e-01 2.570e-143  9.920e-01     1
## ---
## Formula alpha:
## ---
## alpha ~ s(years, k = 20)
## -
## Parametric coefficients:
##              Mean  2.5%   50% 97.5% parameters
## (Intercept) 1.661 1.116 1.666 2.139          0
## -
## Acceptance probabilty:
##         Mean   2.5%    50% 97.5%
## alpha 0.9473 0.4908 0.9972     1
## -
## Smooth terms:
##                     Mean      2.5%       50%     97.5% parameters
## s(years).tau21 1.095e+01 3.477e-02 4.575e+00 5.909e+01      0.708
## s(years).edf   4.686e+00 1.510e+00 4.536e+00 7.718e+00      0.000
## s(years).alpha 6.828e-01 1.393e-10 8.595e-01 1.000e+00         NA
## ---
## Formula dalpha:
## ---
## dalpha ~ -1
## ---
## Sampler summary:
## -
## DIC = 5081.882 logLik = -2298.779 logPost = 4514.064
## pd = 484.3245 runtime = 10075.83
## ---
## Optimizer summary:
## -
## converged = 0 logPost = -33648.29 time = 174.306

Here, the estimated parameters from the backfitting step are quite different from the posterior mean calculated from the MCMC samples. This might indicate that the default number of iterations, n.iter = 1200, is not enough to obtain convergence of the MCMC chains. To check this issue we look at some traceplots.

plot(b, which = "samples")

For convenience, we only show the traceplots of the intercept for parameter \(\gamma\) and \(\mu\). The traceplot for parameter \(\mu\) clearly shows an upwards trend of the MCMC chain. To overcome this issue, we estimate the model again and increase the burnin-in phase and the thinning parameter.

set.seed(321)

## Run model
m <- bamlss(f, data = pbc2, family = "jm",
  timevar = "year", idvar = "id", optimizer = FALSE,
  jm.start = coef(b), n.iter = 12000, burnin = 2000, thin = 10)

Note, executing the code will take some time because of the large number of parameters for the functional random intercept terms (\(> 2500\)). Another solution to speed up estimation is to use multiple cores for sampling. This can be done using the bamlss() wrapper by setting, e.g., cores = 4. Then, internally the mclapply() function of the parallel (R Core Team 2019) package is called. Now, the samples from model m look much better.

plot(m, which = "samples", lag = 100)

However, there is still some autocorrelation for parameter \(\mu\), which could be further reduced by, e.g., increasing the thinning parameter.

The estimated effects can be plotted instantly, e.g., for the time-varying parameters \(\lambda\) and \(\alpha\):

plot(m, model = c("lambda", "alpha"))

The left plot is the estimated baseline hazard effect, the right plot the estimated association parameter. Both effects seem to be quite linear, however, according the 95% credible intervals the effect for parameter alpha is not significant.

The age effect for the time-constant parameter \(\gamma\) of the risk model is plotted with:

plot(m, model = "gamma")

Again, the effect of age on the logarithmic risk scale is basically linear and increasing. The estimated effects for the bodymarker model in the \(\mu\) parameter are plotted in similar fashion:

plot(m, model = "mu")

The left plot shows the estimated random effect for each individual, the middle plot the estimated functional random intercepts and the right plot the corresponding main effect of year, which is linearly increasing.

Another feature that is implemented is the plot of estimated survival probabilities for single individuals.

jm_survplot(m, id = 21)

The upper plot shows the predicted survival probability up to 14 years of individual 21. The lower plot is the predicted effect of time in the longitudinal model \(\mu\). Here, the dots represent the observed bilirubin level of this individual.

References

Köhler, Meike, Nikolaus Umlauf, and Sonja Greven. 2018. “Nonlinear Association Structures in Flexible Bayesian Additive Joint Models.” Statistics in Medicine 37 (30): 4771–88. doi:10.1002/sim.7967.

Köhler, Meike, Nikolaus Umlauf, Andreas Beyerlein, Christiane Winkler, Anette-Gabriele Ziegler, and Sonja Greven. 2017. “Flexible Bayesian Additive Joint Models with an Application to Type 1 Diabetes Research.” Biometrical Journal 59 (6): 1144–65. doi:10.1002/bimj.201600224.

R Core Team. 2019. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Rizopoulos, Dimitris. 2016. “The R Package JMbayes for Fitting Joint Models for Longitudinal and Time-to-Event Data Using Mcmc.” Journal of Statistical Software 72 (7): 1–45. doi:10.18637/jss.v072.i07.

Umlauf, Nikolaus, Nadja Klein, Achim Zeileis, and Thorsten Simon. 2019. bamlss: Bayesian Additive Models for Location Scale and Shape (and Beyond). https://CRAN.R-project.org/package=bamlss.