\[ \large\require{color}\mathcal{p}\left(y, t,\colorbox{#F5842070}{$b$}\right)=\mathcal{p}\left( y\mid\colorbox{#F5842070}{$b$}\right)\mathcal{p}\left( t\mid\colorbox{#F5842070}{$b$}\right)\mathcal{p}\left(\colorbox{#F5842070}{$b$}\right), \]

\[ \large\colorbox{#F5842070}{$b$} \sim \mathcal{N} \left(0, D\right). \]

• \(\colorbox{#F5842070}{$b$}\) explain all interrelationships.

• Mathematically convenient, but computationally intensive.

Estimation:

• Bayesian paradigm;

• Metropolis-Hastings and Robbins–Monro algorithms;

Implementation:

• MCMC in C++;

• Parallel sampling of random effects.

• Parallel MCMC chains.

\[{\require{color} \begin{cases} y_i(t)= \colorbox{#F5842070}{$\boldsymbol x_i(t)^\top\boldsymbol\beta + \boldsymbol z_i(t)^\top \mathbf b_i$} + \varepsilon_i(t) = \colorbox{#F5842070}{$\eta_i(t)$} + \varepsilon(t) & \text{Longitudinal process}\\ h_i(t)= h_0(t)\exp\left\{ \boldsymbol w_i(t)^\top \boldsymbol \gamma + \colorbox{#F5842070}{$\eta_i(t)$} \alpha \right\} & \text{Event process}\\ \end{cases} } \]

\[ \mathbf b_i \sim \mathcal{N} \left(\boldsymbol 0, \mathbf D \right), \qquad \varepsilon_i(t) \sim \mathcal{N} \left(0, \sigma^2\right), \qquad i=1,\dots,n \text{ individuals}. \]

\[{\require{color} \begin{cases} y_i(t)= \colorbox{#FFFFFF}{$\boldsymbol x_i(t)^\top\boldsymbol\beta + \boldsymbol z_i(t)^\top \mathbf b_i$} + \varepsilon_i(t) = \colorbox{#FFFFFF}{$\eta_i(t)$} + \varepsilon(t) & \text{Longitudinal process}\\ h_i(t)= h_0(t)\exp\left\{ \boldsymbol w_i(t)^\top \boldsymbol \gamma + \colorbox{#FFFFFF}{$\eta_i(t)$} \alpha \right\} & \text{Event process}\\ \end{cases} } \]

\[ \mathbf b_i \sim \mathcal{N} \left(\boldsymbol 0, \mathbf D \right), \qquad \varepsilon_i(t) \sim \mathcal{N} \left(0, \sigma^2\right), \qquad i=1,\dots,n \text{ individuals}. \]

library(JMbayes2)

long_fit <- lme(y ~ time + group, long_data, ~ time | id)

\[{\require{color} \begin{cases} y_i(t)= \colorbox{#FFFFFF}{$\boldsymbol x_i(t)^\top\boldsymbol\beta + \boldsymbol z_i(t)^\top \mathbf b_i$} + \varepsilon_i(t) = \colorbox{#FFFFFF}{$\eta_i(t)$} + \varepsilon(t) & \text{Longitudinal process}\\ h_i(t)= h_0(t)\exp\left\{ \boldsymbol w_i(t)^\top \boldsymbol \gamma + \colorbox{#FFFFFF}{$\eta_i(t)$} \alpha \right\} & \text{Event process}\\ \end{cases} } \]

\[ \mathbf b_i \sim \mathcal{N} \left(\boldsymbol 0, \mathbf D \right), \qquad \varepsilon_i(t) \sim \mathcal{N} \left(0, \sigma^2\right), \qquad i=1,\dots,n \text{ individuals}. \]

library(JMbayes2)

long_fit <- lme(y ~ time + group, long_data, ~ time | id)

ter_fit  <- coxph(Surv(tstop, status) ~ group, ter_data)

\[{\require{color} \begin{cases} y_i(t)= \colorbox{#FFFFFF}{$\boldsymbol x_i(t)^\top\boldsymbol\beta + \boldsymbol z_i(t)^\top \mathbf b_i$} + \varepsilon_i(t) = \colorbox{#FFFFFF}{$\eta_i(t)$} + \varepsilon(t) & \text{Longitudinal process}\\ h_i(t)= h_0(t)\exp\left\{ \boldsymbol w_i(t)^\top \boldsymbol \gamma + \colorbox{#FFFFFF}{$\eta_i(t)$} \alpha \right\} & \text{Event process}\\ \end{cases} } \]

\[ \mathbf b_i \sim \mathcal{N} \left(\boldsymbol 0, \mathbf D \right), \qquad \varepsilon_i(t) \sim \mathcal{N} \left(0, \sigma^2\right), \qquad i=1,\dots,n \text{ individuals}. \]

library(JMbayes2)

long_fit <- lme(y ~ time + group, long_data, ~ time | id)

ter_fit  <- coxph(Surv(tstop, status) ~ group, ter_data)

jm_fit   <- jm(ter_fit, long_fit, 'time')

summary(jm_fit)

 ## Call:
 ## jm(Surv_object = ter_fit, Mixed_objects = long_fit, 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    4400.430 5288.656 -4037.448
 ## conditional 4571.205 6412.573 -5511.117
 ## 
 ## Random-effects covariance matrix:
 ##                     
 ##        StdDev   Corr
 ## (Intr) 1.0100 (Intr)
 ## year   0.1823 0.3714
 ## 
 ## Survival Outcome:
 ##                         Mean  StDev    2.5%  97.5%      P   Rhat
 ## sexfemale            -0.2747 0.3945 -1.0233 0.5175 0.4882 1.0138
 ## value(log(serBilir))  1.2126 0.1100  0.9932 1.4357 0.0000 1.0082
 ## 
 ## Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
 ##                Mean  StDev    2.5%  97.5%      P   Rhat
 ## (Intercept)  0.6525 0.3062  0.0490 1.2563 0.0384 1.0051
 ## year         0.1834 0.0293  0.1257 0.2399 0.0000 1.0017
 ## sexfemale   -0.1803 0.3175 -0.8024 0.4480 0.5556 1.0054
 ## sigma        0.3864 0.0176  0.3586 0.4299 0.0000 1.0084
 ## 
 ## MCMC summary:
 ## chains: 3 
 ## iterations per chain: 3500 
 ## burn-in per chain: 500 
 ## thinning: 1 
 ## time: 35 sec

traceplot(jm_fit, "alphas")

\[\require{color} \begin{cases} y_i(t)= \colorbox{#D36EAB70}{$\color{#515151} \eta_i(t)$} + \varepsilon_i(t) \\ h_i(t)= h_0(t)\exp\left\{ w_i(t)^\top \gamma + \stackrel{}{\boxed{\;\mathcal{f}\left\{\colorbox{#D36EAB70}{$\color{#515151}\eta_i(t)$}\right\}\;}} \alpha \right\}\\ \end{cases}, \]

\[ \varepsilon_i(t) \sim \mathcal{N} \left(0, \sigma^2\right), \qquad i=1,\dots,n \text{ individuals}. \]

jm_fit <- update(jm_fit, functional_forms = ~ value(y))

\[\require{color}\DeclareMathOperator{\logit}{logit} \begin{cases} y_{\colorbox{#F5842070}{$k_i$}}(t)= \mathcal{g}^{-1}_k\left\{\eta_{k_i}(t)\right\} \\ \\ h_i(t)= h_0(t)\exp\left\{w_i(t)^\top \gamma + \colorbox{#F5842070}{$\sum_{k=1}^{p}$}\mathcal{f}_k\left\{\eta_{k_i}(t)\right\} \alpha_k \right\} \\ \end{cases}, \]

\[ \mathbf b \sim \mathcal{N} \left(\boldsymbol 0, \mathbf D \right), \qquad i=1,\dots,n \text{ individuals}, \qquad \colorbox{#F5842070}{$k=1,\dots,p \text{ biomarkers.}$} \]

library(JMbayes2)

long_fit  <- lme(y ~ time + group, long_data, ~ time | id)

ter_fit   <- coxph(Surv(tstop, status) ~ group, ter_data)

\[\require{color}\DeclareMathOperator{\logit}{logit} \begin{cases} y_{\colorbox{#FFFFFF}{$k_i$}}(t)= \mathcal{g}^{-1}_k\left\{\eta_{k_i}(t)\right\} \\ \\ h_i(t)= h_0(t)\exp\left\{w_i(t)^\top \gamma + \colorbox{#FFFFFF}{$\sum_{k=1}^{p}$}\mathcal{f}_k\left\{\eta_{k_i}(t)\right\} \alpha_k \right\} \\ \end{cases}, \]

\[ \mathbf b \sim \mathcal{N} \left(\boldsymbol 0, \mathbf D \right), \qquad i=1,\dots,n \text{ individuals}, \qquad \colorbox{#FFFFFF}{$k=1,\dots,p \text{ biomarkers.}$} \]

library(JMbayes2)

long_fit  <- lme(y ~ time + group, long_data, ~ time | id)

ter_fit   <- coxph(Surv(tstop, status) ~ group, ter_data)

long_fit2 <- mixed_model(y2 ~ time, ~ time | id, long_data, binomial() )

\[\require{color}\DeclareMathOperator{\logit}{logit} \begin{cases} y_{\colorbox{#FFFFFF}{$k_i$}}(t)= \mathcal{g}^{-1}_k\left\{\eta_{k_i}(t)\right\} \\ \\ h_i(t)= h_0(t)\exp\left\{w_i(t)^\top \gamma + \colorbox{#FFFFFF}{$\sum_{k=1}^{p}$}\mathcal{f}_k\left\{\eta_{k_i}(t)\right\} \alpha_k \right\} \\ \end{cases}, \]

\[ \mathbf b \sim \mathcal{N} \left(\boldsymbol 0, \mathbf D \right), \qquad i=1,\dots,n \text{ individuals}, \qquad \colorbox{#FFFFFF}{$k=1,\dots,p \text{ biomarkers.}$} \]

library(JMbayes2)

long_fit  <- lme(y ~ time + group, long_data, ~ time | id)

ter_fit   <- coxph(Surv(tstop, status) ~ group, ter_data)

long_fit2 <- mixed_model(y2 ~ time, ~ time | id, long_data, binomial() )

jm_fit    <- jm(ter_fit, list(long_fit, long_fit2), 'time')

Various distributions:

Univariate and Multivariate Joint Models and Non Gaussian Mixed Models vignettes.

\[\require{color} \begin{cases} y_{k_i}(t)= \mathcal{g}\left\{\eta^{-1}_{k_i}(t)\right\} \\ \\ h_{\colorbox{#D36EAB70}{$j_i$}}(t)= \,h_{\colorbox{#D36EAB70}{$j_0$}}(t)\exp\left\{ w_{j_i}(t)^\top \gamma_j \,+ \sum_{k=1}^{p}\mathcal{f}_{jk}\left\{\eta_{k_i}(t)\right\} \alpha_{jk} \,\colorbox{#F5842070}{$\left(+\,\upsilon_{j_i}\right)$} \right\} \\ \end{cases} \]

\[ \small \begin{pmatrix} \mathbf b \\ \colorbox{#F5842070}{$\upsilon_i$}\end{pmatrix} \sim \mathcal{N} \left(\mathbf 0, \begin{pmatrix} \mathbf D & 0 \\ & \sigma^2_F \\ \end{pmatrix}\right), \]

\[ \quad i=1,\dots,n \text{ individuals}, \qquad k=1,\dots,p \text{ biomarkers}, \qquad \colorbox{#D36EAB70}{$j=1,\dots,m \text{ events.}$} \]

Competing Risks and Multi-state Processes vignettes.

\[\require{color} \begin{cases} y_{k_i}(t)= \mathcal{g}\left\{\eta^{-1}_{k_i}(t)\right\} \\ \\ h_{T_i}(t)= \,h_{T_0}(t)\exp\left\{ w_{T_i}(t)^\top \gamma_T \,+ \sum_{k=1}^{p}\mathcal{f}_k\left\{\eta_{k_i}(t)\right\} \alpha_{T_k} \,+ \colorbox{#F5842070}{$\upsilon_i$} \alpha_{F} \right\} & \text{Terminal}\\ \\ \colorbox{#EBEBEB}{$h_{R_i}(t)= h_{R_0}(t)\exp\left\{ w_{R_i}(t)^\top \gamma_R + \sum_{k=1}^{p}\mathcal{f}_k\left\{\eta_{k_i}(t)\right\} \alpha_{R_k} + \colorbox{#F5842070}{$\upsilon_i$} \right\}$} & \text{Recurrent}\\ \end{cases} \]

\[ \begin{pmatrix} \mathbf b \\ \colorbox{#F5842070}{$\upsilon_i$}\end{pmatrix} \sim \mathcal{N} \left(\mathbf 0, \begin{pmatrix} \mathbf D & 0 \\ & \sigma^2_F \\ \end{pmatrix}\right),\]

\[ i=1,\dots,n \text{ individuals}, \qquad k=1,\dots,p \text{ biomarkers}. \]

 ##   id tstart tstop status strata
  ## 1  1      0     2      1    rec 
 ## 2  1      3     4      1    rec
 ## 3  1      5     7      0    rec
 ## 4  1      0     7      0    ter

surv <- rc_setup(rc_data, trm_data) 

library(JMbayes2)

long_fit  <- lme(y ~ time, long, ~ time | id)

surv_fit  <- coxph(Surv(tstart, tstop, status) ~ sex:strata(process), surv)

library(JMbayes2)

long_fit  <- lme(y ~ time, long, ~ time | id)

surv_fit  <- coxph(Surv(tstart, tstop, status) ~ sex:strata(process), surv)

jm_fit    <- jm(surv_fit, long_fit, 'time',

                functional_forms = ~ value(y):process,
                
                recurrent = 'gap')

summary(jm_fit)

 ## 
 ## Call:
 ## jm(Surv_object = ter_fit, Mixed_objects = long_fit, time_var = "time", 
 ##     recurrent = "gap", functional_forms = ~value(y):strata)
 ## 
 ## Data Descriptives:
 ## Number of Groups: 500        Number of events: 2250 (81.6%)
 ## Number of Observations:
 ##   y: 3075
 ## 
 ##                  DIC     WAIC      LPML
 ## marginal    13501.29 13337.81 -7021.504
 ## conditional 14644.81 14287.74 -7842.206
 ## 
 ## Random-effects covariance matrix:
 ##                      
 ##        StdDev   Corr 
 ## (Intr) 0.6737 (Intr) 
 ## time   0.4474 -0.0461
 ## 
 ## Frailty standard deviation:
 ##                 Mean   2.5%  97.5%
 ## sigma_frailty 0.6107 0.5249 0.7006
 ## 
 ## Survival Outcome:
 ##                           Mean  StDev   2.5%  97.5%     P   Rhat
 ## group:strata(strata)Rec 0.5354 0.1194 0.2977 0.7661 0e+00 1.0077
 ## group:strata(strata)Ter 0.4452 0.1193 0.2120 0.7023 2e-04 1.0127
 ## value(y):strataRec      0.4821 0.0309 0.4281 0.5467 0e+00 1.5571
 ## value(y):strataTer      0.4031 0.0454 0.3265 0.4964 0e+00 1.4754
 ## frailty                 0.6524 0.1249 0.4164 0.9093 0e+00 1.0072
 ## 
 ## Longitudinal Outcome: y (family = gaussian, link = identity)
 ##               Mean  StDev   2.5%  97.5% P   Rhat
 ## (Intercept) 6.9147 0.1160 6.6809 7.1357 0 1.0008
 ## time        0.2867 0.0619 0.1685 0.4110 0 1.0006
 ## sigma       0.7154 0.0123 0.6923 0.7411 0 1.0002
 ## 
 ## MCMC summary:
 ## chains: 3 
 ## iterations per chain: 3500 
 ## burn-in per chain: 500 
 ## thinning: 1 
 ## time: 3.5 min

Recurrent Events vignette.

pred_long <- predict(jm_fit, new_data, process = 'longitudinal',

                     times = seq(tstart, 12))

pred_long <- predict(jm_fit, new_data, process = 'longitudinal',

                     times = seq(tstart, 12))

pred_ter <- predict(jm_fit, new_data, process = 'event',

                    times = seq(tstart, 12))

pred_ter <- predict(jm_fit, new_data, process = 'event',

                    times = seq(tstart, 12))

• Predictive Accuracy

Dynamic Predictions vignette.

• Shrinkage priors.

• Combine competing risks and recurrent events.

• Dynamic predicitons:

  • Joint models with recurrent events;

  • Multistate joint models.