JMbayes2

Extended Joint Models under the Bayesian Approach

Dimitris Rizopoulos, P. Afonso, G. Papageorgiou

Department of Biostatistics

d.rizopoulos@erasmusmc.nl drizopoulos drizopoulos

1 / 26

Joint Models:
Current Standard

1 / 26

A Bit of History

2 / 26

A Bit of History

Joint models have been around for many years...

2 / 26

A Bit of History

Joint models have been around for many years...

Two parallel trajectories

Event Process

endogenous time-varying covariates

Longitudinal Outcomes

non-random (MNAR) dropout

2 / 26

JMs and Extensions

Basic assumption: Conditional Independence
The random effects

b

explain all interrelationships.

Rizopoulos (2012)

3 / 26

JMs and Extensions

Basic assumption: Conditional Independence
The random effects

b

explain all interrelationships.

Rizopoulos (2012)

$Y$ : Longitudinal process
$T$ : Event process

3 / 26

JMs and Extensions

Basic assumption: Conditional Independence
The random effects

b

explain all interrelationships.

Rizopoulos (2012)

$Y$ : Longitudinal process
$T$ : Event process

$\begin{aligned} p {Y, T ∣ b} = p {Y ∣ b} p {T ∣ b} \end{aligned}$

3 / 26

JMs and Extensions

DynPreds

4 / 26

JMs and Extensions

Longitudinal Process

multiple outcomes
continuous, discrete, left-censored
flexible models (splines, etc.)

Event Process

competing risks
multi-state processes
recurrent events
multivariate

5 / 26

JMs and Extensions

Conditional Independence: Rizopoulos (2012)

mathematically convenient
computionally intensive

6 / 26

JMs and Extensions

Conditional Independence: Rizopoulos (2012)

mathematically convenient
computionally intensive

Maximum Likelihood:

requires high-dimensional numerical integration

Bayesian:

requires sampling high-dimensional random effects

6 / 26

Theory vs Practice

As with many other statistical techniques,
it took some years for JMs software to appear

7 / 26

Theory vs Practice

As with many other statistical techniques,
it took some years for JMs software to appear

Nowadays, given the rise of the R ecosystem a number of packages is available

7 / 26

Available Software

Maximum Likelihood:

Bayesian:

8 / 26

Available Software

Still no software package exists to satisfy all needs

9 / 26

JMbayes2

10 / 26

A New Package

⇨ Learning from the lessons of the past we decided to create a new package

11 / 26

A New Package

⇨ Learning from the lessons of the past we decided to create a new package

Our aims

versatile
user-friendly
fast
complete

cover (almost) all extensions
straightforward syntax
C++ & computational tricks
support functions

11 / 26

Under the Hood

Work horse: Metropolis-Hastings algorithm
⇨ adaptive optimal scaling using the Robbins–Monro process

12 / 26

Under the Hood

Work horse: Metropolis-Hastings algorithm
⇨ adaptive optimal scaling using the Robbins–Monro process

Computational aspects
⇨ parallel sampling of random effects
⇨ info from separate fits

12 / 26

Under the Hood

Implementation
⇨ MCMC written in C++
⇨ avoiding double computations

13 / 26

An Example

library("JMbayes2")
fm1 <- lme(log(serBilir) ~ year * sex, data = pbc2, random = ~ year | id)
CoxFit <- coxph(Surv(years, status2) ~ sex, data = pbc2.id)
jointFit1 <- jm(CoxFit, fm1, time_var = "year")

14 / 26

An Example

library("JMbayes2")
fm1 <- lme(log(serBilir) ~ year * sex, data = pbc2, random = ~ year | id)
CoxFit <- coxph(Surv(years, status2) ~ sex, data = pbc2.id)
jointFit1 <- jm(CoxFit, fm1, time_var = "year")

Check convergence

ggtraceplot(jointFit1, "alphas")

14 / 26

Traceplot

15 / 26

Model Summary

summary(jointFit1)

## 
## Call:
## jm(Surv_object = CoxFit, Mixed_objects = fm1, time_var = "year")
## 
## Data Descriptives:
## Number of Groups: 312        Number of events: 140 (44.9%)
## Number of Observations:
##   log(serBilir): 1945
## 
##                  DIC     WAIC      LPML
## marginal    4255.243 5193.189 -2927.695
## conditional 3438.104 3256.032 -1840.069
## 
## Random-effects covariance matrix:
##                     
##        StdDev   Corr
## (Intr) 1.0053 (Intr)
## year   0.1810 0.3976
## 
## Survival Outcome:
##                         Mean  StDev    2.5%  97.5%      P   Rhat
## sexfemale            -0.2626 0.2867 -0.7847 0.3157 0.3689 1.0020
## value(log(serBilir))  1.2475 0.0947  1.0613 1.4390 0.0000 1.0057
## 
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
##                   Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept)     0.7223 0.1720  0.3826  1.0622 0.0000 1.0000
## year            0.2638 0.0376  0.1912  0.3395 0.0000 1.0055
## sexfemale      -0.2586 0.1824 -0.6156  0.0971 0.1620 1.0000
## year:sexfemale -0.0899 0.0402 -0.1705 -0.0120 0.0229 1.0047
## sigma           0.3469 0.0068  0.3337  0.3607 0.0000 1.0207
## 
## MCMC summary:
## chains: 3 
## iterations per chain: 3500 
## burn-in per chain: 500 
## thinning: 1 
## time: 18 sec

16 / 26

Multivariate

fm2 <- lme(prothrombin ~ year * sex, data = pbc2, random = ~ year | id)
fm3 <- mixed_model(ascites ~ year + sex, data = pbc2,
                   random = ~ year | id, family = binomial())
jointFit2 <- jm(CoxFit, list(fm1, fm2, fm3), time_var = "year",
                n_iter = 12000L, n_burnin = 2000L, n_thin = 5L)

17 / 26

Multivariate

summary(jointFit2)

## 
## Call:
## jm(Surv_object = CoxFit, Mixed_objects = list(fm1, fm2, fm3), 
##     time_var = "year", n_iter = 12000L, n_burnin = 2000L, n_thin = 5L)
## 
## Data Descriptives:
## Number of Groups: 312        Number of events: 140 (44.9%)
## Number of Observations:
##   log(serBilir): 1945
##   prothrombin: 1945
##   ascites: 1885
## 
##                  DIC      WAIC       LPML
## marginal    11501.71 592782.35 -39636.864
## conditional 12491.57  12237.35  -6618.324
## 
## Random-effects covariance matrix:
##                                                  
##        StdDev   Corr                             
## (Intr) 0.9913 (Intr)   year (Intr)   year (Intr) 
## year   0.1899 0.4258                             
## (Intr) 0.7774 0.5618 0.4965                      
## year   0.3310 0.4145 0.3559 0.0587               
## (Intr) 2.9639 0.5629 0.5285 0.6224 0.2768        
## year   0.4601 0.3954 0.6588 0.2551 0.4505 -0.0095
## 
## Survival Outcome:
##                         Mean  StDev    2.5%   97.5%      P   Rhat
## sexfemale            -0.8052 0.3825 -1.5633 -0.0734 0.0347 1.0019
## value(log(serBilir))  0.4671 0.1765  0.1134  0.7965 0.0117 1.0737
## value(prothrombin)    0.0463 0.1518 -0.2836  0.3178 0.6847 1.0408
## value(ascites)        0.6305 0.1388  0.3947  0.9281 0.0000 1.0245
## 
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
##                   Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept)     0.7014 0.1681  0.3689  1.0278 0.0003 1.0008
## year            0.2689 0.0340  0.2017  0.3365 0.0000 1.0000
## sexfemale      -0.2395 0.1784 -0.5871  0.1087 0.1840 1.0006
## year:sexfemale -0.0794 0.0360 -0.1508 -0.0082 0.0280 1.0031
## sigma           0.3466 0.0066  0.3340  0.3597 0.0000 1.0011
## 
## Longitudinal Outcome: prothrombin (family = gaussian, link = identity)
##                   Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept)    10.9424 0.1693 10.6140 11.2733 0.0000 1.0018
## year            0.2436 0.0741  0.0980  0.3915 0.0003 1.0002
## sexfemale      -0.4000 0.1810 -0.7613 -0.0475 0.0240 1.0016
## year:sexfemale  0.0549 0.0781 -0.0992  0.2103 0.4783 1.0013
## sigma           1.0551 0.0200  1.0164  1.0952 0.0000 1.0085
## 
## Longitudinal Outcome: ascites (family = binomial, link = logit)
##                Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept) -4.5003 0.7081 -5.9580 -3.1942 0.0000 1.0087
## year         0.6081 0.0732  0.4700  0.7582 0.0000 1.0402
## sexfemale   -0.5043 0.6911 -1.8435  0.8308 0.4613 1.0045
## 
## MCMC summary:
## chains: 3 
## iterations per chain: 12000 
## burn-in per chain: 2000 
## thinning: 5 
## time: 2.1 min

18 / 26

Functional Forms

jointFit3 <- update(jointFit2, 
  functional_forms = ~ slope(log(serBilir)) + slope(log(serBilir)):sex + 
    area(prothrombin))

19 / 26

Functional Forms

jointFit3 <- update(jointFit2, 
  functional_forms = ~ slope(log(serBilir)) + slope(log(serBilir)):sex + 
    area(prothrombin))

Slope $h_{i} (t) = h_{0} (t) \exp {α \frac{d η_{i} (t)}{d t}}$

(normalized) Area $h_{i} (t) = h_{0} (t) \exp {α \frac{1}{t} \int_{0}^{t} η_{i} (s) d s}$

η_{i} (t) = x_{i} (t)^{⊤} β + z_{i} (t)^{⊤} b_{i}

19 / 26

Functional Forms

summary(jointFit3)

## 
## Call:
## jm(Surv_object = CoxFit, Mixed_objects = list(fm1, fm2, fm3), 
##     time_var = "year", functional_forms = ~slope(log(serBilir)) + 
##         slope(log(serBilir)):sex + area(prothrombin), n_iter = 12000L, ...
## 
## Data Descriptives:
## Number of Groups: 312        Number of events: 140 (44.9%)
## Number of Observations:
##   log(serBilir): 1945
##   prothrombin: 1945
##   ascites: 1885
## 
##                  DIC     WAIC      LPML
## marginal    11407.35 12593.53 -6801.446
## conditional 12464.61 12210.27 -6597.434
## 
## Random-effects covariance matrix:
##                                                 
##        StdDev   Corr                            
## (Intr) 0.9839 (Intr)   year (Intr)   year (Intr)
## year   0.1920 0.4627                            
## (Intr) 0.7730 0.5623 0.5129                     
## year   0.3329 0.4331 0.3694 0.0773              
## (Intr) 2.9335 0.5821 0.5332 0.6110 0.2939       
## year   0.4724 0.4700 0.6992 0.3232 0.4553 0.0607
## 
## Survival Outcome:
##                                   Mean  StDev    2.5%  97.5%      P   Rhat
## sexfemale                      -0.2101 0.7811 -1.7096 1.3984 0.7723 1.0581
## slope(log(serBilir))            4.1106 2.0171 -0.0720 7.8839 0.0543 1.0085
## slope(log(serBilir)):sexfemale -2.2855 2.2898 -7.1533 2.0401 0.3127 1.0604
## area(prothrombin)               0.1355 0.2434 -0.3670 0.5848 0.5390 1.1560
## value(ascites)                  0.7797 0.1901  0.4654 1.1858 0.0000 1.2316
## 
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
##                   Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept)     0.6770 0.1658  0.3506  0.9979 0.0003 1.0005
## year            0.2642 0.0338  0.1977  0.3311 0.0000 1.0004
## sexfemale      -0.2145 0.1763 -0.5628  0.1313 0.2197 1.0005
## year:sexfemale -0.0688 0.0351 -0.1364 -0.0001 0.0500 1.0010
## sigma           0.3474 0.0065  0.3350  0.3600 0.0000 1.0001
## 
## Longitudinal Outcome: prothrombin (family = gaussian, link = identity)
##                   Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept)    10.9271 0.1716 10.5903 11.2658 0.0000 1.0016
## year            0.2387 0.0744  0.0969  0.3896 0.0017 1.0025
## sexfemale      -0.3811 0.1812 -0.7357 -0.0270 0.0357 1.0014
## year:sexfemale  0.0639 0.0781 -0.0880  0.2200 0.4093 1.0027
## sigma           1.0551 0.0203  1.0157  1.0956 0.0000 1.0102
## 
## Longitudinal Outcome: ascites (family = binomial, link = logit)
##                Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept) -4.6078 0.7255 -6.1072 -3.2645 0.0000 1.0750
## year         0.6355 0.0768  0.5003  0.7970 0.0000 1.1650
## sexfemale   -0.3989 0.6978 -1.7520  0.9894 0.5607 1.0013
## 
## MCMC summary:
## chains: 3 
## iterations per chain: 12000 
## burn-in per chain: 2000 
## thinning: 5 
## time: 2.3 min

20 / 26

Competing Risks

Bring the data in the competing risks long format

pbc2.idCR <- crisk_setup(pbc2.id, statusVar = "status", censLevel = "alive", 
                    nameStrata = "CR")
pbc2.idCR[pbc2.idCR$id %in% c(1, 2, 5), c("id", "years", "status", "status2", "CR")]

##     id     years       status status2           CR
## 1    1  1.095170         dead       1         dead
## 1.1  1  1.095170         dead       0 transplanted
## 2    2 14.152338        alive       0         dead
## 2.1  2 14.152338        alive       0 transplanted
## 5    5  4.120578 transplanted       0         dead
## 5.1  5  4.120578 transplanted       1 transplanted

21 / 26

Competing Risks

Fit a competing risks Cox model

CoxFit_CR <- coxph(Surv(years, status2) ~ (age + drug) * strata(CR),
                     data = pbc2.idCR)

22 / 26

Competing Risks

Fit the competing risks joint model

jFit_CR <- jm(CoxFit_CR, fm1, time_var = "year", 
              functional_forms = ~ value(log(serBilir)):CR, 
              n_iter = 25000L, n_burnin = 5000L, n_thin = 5L)

23 / 26

Competing Risks

summary(jFit_CR)

## 
## Call:
## jm(Surv_object = CoxFit_CR, Mixed_objects = fm1, time_var = "year", 
##     functional_forms = ~value(log(serBilir)):CR, n_iter = 25000L, 
##     n_burnin = 5000L, n_thin = 5L)
## 
## Data Descriptives:
## Number of Groups: 312        Number of events: 169 (27.1%)
## Number of Observations:
##   log(serBilir): 1945
## 
##                  DIC     WAIC      LPML
## marginal    4435.774 5313.425 -3238.973
## conditional 3624.102 3446.399 -2000.840
## 
## Random-effects covariance matrix:
##                     
##        StdDev   Corr
## (Intr) 1.0020 (Intr)
## year   0.1817 0.3973
## 
## Survival Outcome:
##                                        Mean  StDev    2.5%   97.5%      P
## age                                 -0.0681 0.0252 -0.1160 -0.0191 0.0035
## drugD-penicil                       -0.2892 0.3936 -1.0726  0.4564 0.4737
## age:strata(CR)dead                   0.1300 0.0248  0.0808  0.1757 0.0000
## drugD-penicil:strata(CR)dead         0.2697 0.4208 -0.5363  1.1023 0.5288
## value(log(serBilir)):CRtransplanted  1.1675 0.2083  0.7802  1.5974 0.0000
## value(log(serBilir)):CRdead          1.3682 0.1079  1.1643  1.5861 0.0000
##                                       Rhat
## age                                 1.0588
## drugD-penicil                       1.0069
## age:strata(CR)dead                  1.0744
## drugD-penicil:strata(CR)dead        1.0077
## value(log(serBilir)):CRtransplanted 1.0152
## value(log(serBilir)):CRdead         1.0031
## 
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
##                   Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept)     0.7182 0.1717  0.3821  1.0565 0.0000 0.9999
## year            0.2624 0.0372  0.1898  0.3356 0.0000 1.0000
## sexfemale      -0.2591 0.1835 -0.6227  0.0989 0.1607 0.9998
## year:sexfemale -0.0845 0.0395 -0.1628 -0.0079 0.0305 1.0000
## sigma           0.3472 0.0067  0.3344  0.3605 0.0000 0.9998
## 
## MCMC summary:
## chains: 3 
## iterations per chain: 25000 
## burn-in per chain: 5000 
## thinning: 5 
## time: 4.9 min

24 / 26

More Features

× Using the compare_jm() function we can compare fitted joint models using DIC, WAIC and LPML


  fit1 <- jm(CoxModel1, MixedModels1, ...)
  fit2 <- jm(CoxModel2, MixedModels2, ...)
  
  compare_jm(fit1, fit2)

25 / 26

More about JMbayes2 in the dedicated website

https://drizopoulos.github.io/JMbayes2/

26 / 26

Thank you for your attention!

d.rizopoulos@erasmusmc.nl drizopoulos drizopoulos http://www.drizopoulos.com/

26 / 26

JMbayes2

Extended Joint Models under the Bayesian Approach

Dimitris Rizopoulos, P. Afonso, G. Papageorgiou

Department of Biostatistics

d.rizopoulos@erasmusmc.nl drizopoulos drizopoulos

1 / 26

Joint Models:
Current Standard

1 / 26

A Bit of History

2 / 26

A Bit of History

Joint models have been around for many years...

2 / 26

A Bit of History

Joint models have been around for many years...

Two parallel trajectories

Event Process

endogenous time-varying covariates

Longitudinal Outcomes

non-random (MNAR) dropout

2 / 26

JMs and Extensions

Basic assumption: Conditional Independence
The random effects

b

explain all interrelationships.

Rizopoulos (2012)

3 / 26

JMs and Extensions

Basic assumption: Conditional Independence
The random effects

b

explain all interrelationships.

Rizopoulos (2012)

$Y$ : Longitudinal process
$T$ : Event process

3 / 26

JMs and Extensions

Basic assumption: Conditional Independence
The random effects

b

explain all interrelationships.

Rizopoulos (2012)

$Y$ : Longitudinal process
$T$ : Event process

$\begin{aligned} p {Y, T ∣ b} = p {Y ∣ b} p {T ∣ b} \end{aligned}$

3 / 26

JMs and Extensions

DynPreds

4 / 26

JMs and Extensions

Longitudinal Process

multiple outcomes
continuous, discrete, left-censored
flexible models (splines, etc.)

Event Process

competing risks
multi-state processes
recurrent events
multivariate

5 / 26

JMs and Extensions

Conditional Independence: Rizopoulos (2012)

mathematically convenient
computionally intensive

6 / 26

JMs and Extensions

Conditional Independence: Rizopoulos (2012)

mathematically convenient
computionally intensive

Maximum Likelihood:

requires high-dimensional numerical integration

Bayesian:

requires sampling high-dimensional random effects

6 / 26

Theory vs Practice

As with many other statistical techniques,
it took some years for JMs software to appear

7 / 26

Theory vs Practice

As with many other statistical techniques,
it took some years for JMs software to appear

Nowadays, given the rise of the R ecosystem a number of packages is available

7 / 26

Available Software

Maximum Likelihood:

Bayesian:

8 / 26

Available Software

Still no software package exists to satisfy all needs

9 / 26

JMbayes2

10 / 26

A New Package

⇨ Learning from the lessons of the past we decided to create a new package

11 / 26

A New Package

⇨ Learning from the lessons of the past we decided to create a new package

Our aims

versatile
user-friendly
fast
complete

cover (almost) all extensions
straightforward syntax
C++ & computational tricks
support functions

11 / 26

Under the Hood

Work horse: Metropolis-Hastings algorithm
⇨ adaptive optimal scaling using the Robbins–Monro process

12 / 26

Under the Hood

Work horse: Metropolis-Hastings algorithm
⇨ adaptive optimal scaling using the Robbins–Monro process

Computational aspects
⇨ parallel sampling of random effects
⇨ info from separate fits

12 / 26

Under the Hood

Implementation
⇨ MCMC written in C++
⇨ avoiding double computations

13 / 26

An Example

library("JMbayes2")
fm1 <- lme(log(serBilir) ~ year * sex, data = pbc2, random = ~ year | id)
CoxFit <- coxph(Surv(years, status2) ~ sex, data = pbc2.id)
jointFit1 <- jm(CoxFit, fm1, time_var = "year")

14 / 26

An Example

library("JMbayes2")
fm1 <- lme(log(serBilir) ~ year * sex, data = pbc2, random = ~ year | id)
CoxFit <- coxph(Surv(years, status2) ~ sex, data = pbc2.id)
jointFit1 <- jm(CoxFit, fm1, time_var = "year")

Check convergence

ggtraceplot(jointFit1, "alphas")

14 / 26

Traceplot

15 / 26

Model Summary

summary(jointFit1)

## 
## Call:
## jm(Surv_object = CoxFit, Mixed_objects = fm1, time_var = "year")
## 
## Data Descriptives:
## Number of Groups: 312        Number of events: 140 (44.9%)
## Number of Observations:
##   log(serBilir): 1945
## 
##                  DIC     WAIC      LPML
## marginal    4255.243 5193.189 -2927.695
## conditional 3438.104 3256.032 -1840.069
## 
## Random-effects covariance matrix:
##                     
##        StdDev   Corr
## (Intr) 1.0053 (Intr)
## year   0.1810 0.3976
## 
## Survival Outcome:
##                         Mean  StDev    2.5%  97.5%      P   Rhat
## sexfemale            -0.2626 0.2867 -0.7847 0.3157 0.3689 1.0020
## value(log(serBilir))  1.2475 0.0947  1.0613 1.4390 0.0000 1.0057
## 
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
##                   Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept)     0.7223 0.1720  0.3826  1.0622 0.0000 1.0000
## year            0.2638 0.0376  0.1912  0.3395 0.0000 1.0055
## sexfemale      -0.2586 0.1824 -0.6156  0.0971 0.1620 1.0000
## year:sexfemale -0.0899 0.0402 -0.1705 -0.0120 0.0229 1.0047
## sigma           0.3469 0.0068  0.3337  0.3607 0.0000 1.0207
## 
## MCMC summary:
## chains: 3 
## iterations per chain: 3500 
## burn-in per chain: 500 
## thinning: 1 
## time: 18 sec

16 / 26

Multivariate

fm2 <- lme(prothrombin ~ year * sex, data = pbc2, random = ~ year | id)
fm3 <- mixed_model(ascites ~ year + sex, data = pbc2,
                   random = ~ year | id, family = binomial())
jointFit2 <- jm(CoxFit, list(fm1, fm2, fm3), time_var = "year",
                n_iter = 12000L, n_burnin = 2000L, n_thin = 5L)

17 / 26

Multivariate

summary(jointFit2)

## 
## Call:
## jm(Surv_object = CoxFit, Mixed_objects = list(fm1, fm2, fm3), 
##     time_var = "year", n_iter = 12000L, n_burnin = 2000L, n_thin = 5L)
## 
## Data Descriptives:
## Number of Groups: 312        Number of events: 140 (44.9%)
## Number of Observations:
##   log(serBilir): 1945
##   prothrombin: 1945
##   ascites: 1885
## 
##                  DIC      WAIC       LPML
## marginal    11501.71 592782.35 -39636.864
## conditional 12491.57  12237.35  -6618.324
## 
## Random-effects covariance matrix:
##                                                  
##        StdDev   Corr                             
## (Intr) 0.9913 (Intr)   year (Intr)   year (Intr) 
## year   0.1899 0.4258                             
## (Intr) 0.7774 0.5618 0.4965                      
## year   0.3310 0.4145 0.3559 0.0587               
## (Intr) 2.9639 0.5629 0.5285 0.6224 0.2768        
## year   0.4601 0.3954 0.6588 0.2551 0.4505 -0.0095
## 
## Survival Outcome:
##                         Mean  StDev    2.5%   97.5%      P   Rhat
## sexfemale            -0.8052 0.3825 -1.5633 -0.0734 0.0347 1.0019
## value(log(serBilir))  0.4671 0.1765  0.1134  0.7965 0.0117 1.0737
## value(prothrombin)    0.0463 0.1518 -0.2836  0.3178 0.6847 1.0408
## value(ascites)        0.6305 0.1388  0.3947  0.9281 0.0000 1.0245
## 
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
##                   Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept)     0.7014 0.1681  0.3689  1.0278 0.0003 1.0008
## year            0.2689 0.0340  0.2017  0.3365 0.0000 1.0000
## sexfemale      -0.2395 0.1784 -0.5871  0.1087 0.1840 1.0006
## year:sexfemale -0.0794 0.0360 -0.1508 -0.0082 0.0280 1.0031
## sigma           0.3466 0.0066  0.3340  0.3597 0.0000 1.0011
## 
## Longitudinal Outcome: prothrombin (family = gaussian, link = identity)
##                   Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept)    10.9424 0.1693 10.6140 11.2733 0.0000 1.0018
## year            0.2436 0.0741  0.0980  0.3915 0.0003 1.0002
## sexfemale      -0.4000 0.1810 -0.7613 -0.0475 0.0240 1.0016
## year:sexfemale  0.0549 0.0781 -0.0992  0.2103 0.4783 1.0013
## sigma           1.0551 0.0200  1.0164  1.0952 0.0000 1.0085
## 
## Longitudinal Outcome: ascites (family = binomial, link = logit)
##                Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept) -4.5003 0.7081 -5.9580 -3.1942 0.0000 1.0087
## year         0.6081 0.0732  0.4700  0.7582 0.0000 1.0402
## sexfemale   -0.5043 0.6911 -1.8435  0.8308 0.4613 1.0045
## 
## MCMC summary:
## chains: 3 
## iterations per chain: 12000 
## burn-in per chain: 2000 
## thinning: 5 
## time: 2.1 min

18 / 26

Functional Forms

jointFit3 <- update(jointFit2, 
  functional_forms = ~ slope(log(serBilir)) + slope(log(serBilir)):sex + 
    area(prothrombin))

19 / 26

Functional Forms

jointFit3 <- update(jointFit2, 
  functional_forms = ~ slope(log(serBilir)) + slope(log(serBilir)):sex + 
    area(prothrombin))

Slope $h_{i} (t) = h_{0} (t) \exp {α \frac{d η_{i} (t)}{d t}}$

(normalized) Area $h_{i} (t) = h_{0} (t) \exp {α \frac{1}{t} \int_{0}^{t} η_{i} (s) d s}$

η_{i} (t) = x_{i} (t)^{⊤} β + z_{i} (t)^{⊤} b_{i}

19 / 26

Functional Forms

summary(jointFit3)

## 
## Call:
## jm(Surv_object = CoxFit, Mixed_objects = list(fm1, fm2, fm3), 
##     time_var = "year", functional_forms = ~slope(log(serBilir)) + 
##         slope(log(serBilir)):sex + area(prothrombin), n_iter = 12000L, ...
## 
## Data Descriptives:
## Number of Groups: 312        Number of events: 140 (44.9%)
## Number of Observations:
##   log(serBilir): 1945
##   prothrombin: 1945
##   ascites: 1885
## 
##                  DIC     WAIC      LPML
## marginal    11407.35 12593.53 -6801.446
## conditional 12464.61 12210.27 -6597.434
## 
## Random-effects covariance matrix:
##                                                 
##        StdDev   Corr                            
## (Intr) 0.9839 (Intr)   year (Intr)   year (Intr)
## year   0.1920 0.4627                            
## (Intr) 0.7730 0.5623 0.5129                     
## year   0.3329 0.4331 0.3694 0.0773              
## (Intr) 2.9335 0.5821 0.5332 0.6110 0.2939       
## year   0.4724 0.4700 0.6992 0.3232 0.4553 0.0607
## 
## Survival Outcome:
##                                   Mean  StDev    2.5%  97.5%      P   Rhat
## sexfemale                      -0.2101 0.7811 -1.7096 1.3984 0.7723 1.0581
## slope(log(serBilir))            4.1106 2.0171 -0.0720 7.8839 0.0543 1.0085
## slope(log(serBilir)):sexfemale -2.2855 2.2898 -7.1533 2.0401 0.3127 1.0604
## area(prothrombin)               0.1355 0.2434 -0.3670 0.5848 0.5390 1.1560
## value(ascites)                  0.7797 0.1901  0.4654 1.1858 0.0000 1.2316
## 
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
##                   Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept)     0.6770 0.1658  0.3506  0.9979 0.0003 1.0005
## year            0.2642 0.0338  0.1977  0.3311 0.0000 1.0004
## sexfemale      -0.2145 0.1763 -0.5628  0.1313 0.2197 1.0005
## year:sexfemale -0.0688 0.0351 -0.1364 -0.0001 0.0500 1.0010
## sigma           0.3474 0.0065  0.3350  0.3600 0.0000 1.0001
## 
## Longitudinal Outcome: prothrombin (family = gaussian, link = identity)
##                   Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept)    10.9271 0.1716 10.5903 11.2658 0.0000 1.0016
## year            0.2387 0.0744  0.0969  0.3896 0.0017 1.0025
## sexfemale      -0.3811 0.1812 -0.7357 -0.0270 0.0357 1.0014
## year:sexfemale  0.0639 0.0781 -0.0880  0.2200 0.4093 1.0027
## sigma           1.0551 0.0203  1.0157  1.0956 0.0000 1.0102
## 
## Longitudinal Outcome: ascites (family = binomial, link = logit)
##                Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept) -4.6078 0.7255 -6.1072 -3.2645 0.0000 1.0750
## year         0.6355 0.0768  0.5003  0.7970 0.0000 1.1650
## sexfemale   -0.3989 0.6978 -1.7520  0.9894 0.5607 1.0013
## 
## MCMC summary:
## chains: 3 
## iterations per chain: 12000 
## burn-in per chain: 2000 
## thinning: 5 
## time: 2.3 min

20 / 26

Competing Risks

Bring the data in the competing risks long format

pbc2.idCR <- crisk_setup(pbc2.id, statusVar = "status", censLevel = "alive", 
                    nameStrata = "CR")
pbc2.idCR[pbc2.idCR$id %in% c(1, 2, 5), c("id", "years", "status", "status2", "CR")]

##     id     years       status status2           CR
## 1    1  1.095170         dead       1         dead
## 1.1  1  1.095170         dead       0 transplanted
## 2    2 14.152338        alive       0         dead
## 2.1  2 14.152338        alive       0 transplanted
## 5    5  4.120578 transplanted       0         dead
## 5.1  5  4.120578 transplanted       1 transplanted

21 / 26

Competing Risks

Fit a competing risks Cox model

CoxFit_CR <- coxph(Surv(years, status2) ~ (age + drug) * strata(CR),
                     data = pbc2.idCR)

22 / 26

Competing Risks

Fit the competing risks joint model

jFit_CR <- jm(CoxFit_CR, fm1, time_var = "year", 
              functional_forms = ~ value(log(serBilir)):CR, 
              n_iter = 25000L, n_burnin = 5000L, n_thin = 5L)

23 / 26

Competing Risks

summary(jFit_CR)

## 
## Call:
## jm(Surv_object = CoxFit_CR, Mixed_objects = fm1, time_var = "year", 
##     functional_forms = ~value(log(serBilir)):CR, n_iter = 25000L, 
##     n_burnin = 5000L, n_thin = 5L)
## 
## Data Descriptives:
## Number of Groups: 312        Number of events: 169 (27.1%)
## Number of Observations:
##   log(serBilir): 1945
## 
##                  DIC     WAIC      LPML
## marginal    4435.774 5313.425 -3238.973
## conditional 3624.102 3446.399 -2000.840
## 
## Random-effects covariance matrix:
##                     
##        StdDev   Corr
## (Intr) 1.0020 (Intr)
## year   0.1817 0.3973
## 
## Survival Outcome:
##                                        Mean  StDev    2.5%   97.5%      P
## age                                 -0.0681 0.0252 -0.1160 -0.0191 0.0035
## drugD-penicil                       -0.2892 0.3936 -1.0726  0.4564 0.4737
## age:strata(CR)dead                   0.1300 0.0248  0.0808  0.1757 0.0000
## drugD-penicil:strata(CR)dead         0.2697 0.4208 -0.5363  1.1023 0.5288
## value(log(serBilir)):CRtransplanted  1.1675 0.2083  0.7802  1.5974 0.0000
## value(log(serBilir)):CRdead          1.3682 0.1079  1.1643  1.5861 0.0000
##                                       Rhat
## age                                 1.0588
## drugD-penicil                       1.0069
## age:strata(CR)dead                  1.0744
## drugD-penicil:strata(CR)dead        1.0077
## value(log(serBilir)):CRtransplanted 1.0152
## value(log(serBilir)):CRdead         1.0031
## 
## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
##                   Mean  StDev    2.5%   97.5%      P   Rhat
## (Intercept)     0.7182 0.1717  0.3821  1.0565 0.0000 0.9999
## year            0.2624 0.0372  0.1898  0.3356 0.0000 1.0000
## sexfemale      -0.2591 0.1835 -0.6227  0.0989 0.1607 0.9998
## year:sexfemale -0.0845 0.0395 -0.1628 -0.0079 0.0305 1.0000
## sigma           0.3472 0.0067  0.3344  0.3605 0.0000 0.9998
## 
## MCMC summary:
## chains: 3 
## iterations per chain: 25000 
## burn-in per chain: 5000 
## thinning: 5 
## time: 4.9 min

24 / 26

More Features

× Using the compare_jm() function we can compare fitted joint models using DIC, WAIC and LPML


  fit1 <- jm(CoxModel1, MixedModels1, ...)
  fit2 <- jm(CoxModel2, MixedModels2, ...)
  
  compare_jm(fit1, fit2)

25 / 26

More about JMbayes2 in the dedicated website

https://drizopoulos.github.io/JMbayes2/

26 / 26

Thank you for your attention!

d.rizopoulos@erasmusmc.nl drizopoulos drizopoulos http://www.drizopoulos.com/

26 / 26

↑, ←, Pg Up, k	Go to previous slide
↓, →, Pg Dn, Space, j	Go to next slide
Home	Go to first slide
End	Go to last slide
Number + Return	Go to specific slide
b / m / f	Toggle blackout / mirrored / fullscreen mode
c	Clone slideshow
p	Toggle presenter mode
t	Restart the presentation timer
?, h	Toggle this help
o	Tile View: Overview of Slides
Alt + f	Fit Slides to Screen

Joint Models:Current Standard

A Bit of History

A Bit of History

A Bit of History

JMs and Extensions

JMs and Extensions

JMs and Extensions

JMs and Extensions

JMs and Extensions

JMs and Extensions

JMs and Extensions

Theory vs Practice

Theory vs Practice

Available Software

Available Software

JMbayes2

A New Package

A New Package

Under the Hood

Under the Hood

Under the Hood

An Example

An Example

Traceplot

Model Summary

Multivariate

Multivariate

Functional Forms

Functional Forms

Functional Forms

Competing Risks

Competing Risks

Competing Risks

Competing Risks

More Features

Thank you for your attention!

Joint Models:Current Standard

Help

JMbayes2

Joint Models:Current Standard

A Bit of History

A Bit of History

A Bit of History

JMs and Extensions

JMs and Extensions

JMs and Extensions

JMs and Extensions

JMs and Extensions

JMs and Extensions

JMs and Extensions

Theory vs Practice

Theory vs Practice

Available Software

Available Software

JMbayes2

A New Package

A New Package

Under the Hood

Under the Hood

Under the Hood

An Example

An Example

Traceplot

Model Summary

Multivariate

Multivariate

Functional Forms

Functional Forms

Functional Forms

Competing Risks

Competing Risks

Competing Risks

Competing Risks

More Features

Thank you for your attention!

Joint Models:
Current Standard

Joint Models:
Current Standard

Joint Models:
Current Standard