`vignettes/JMbayes2.Rmd`

`JMbayes2.Rmd`

The function that fits joint models in **JMbayes2** is
called `jm()`

. It has three required arguments,
`Surv_object`

a Cox model fitted by `coxph()`

or
an Accelerated Failure time model fitted by `survreg()`

,
`Mixed_objects`

a single or a list of mixed models fitted
either by the `lme()`

or `mixed_model()`

functions, and `time_var`

a character string indicating the
name of the time variable in the specification of the mixed-effects
models. We will illustrate the basic use of the package in the PBC
dataset. We start by fitting a Cox model for the composite event
transplantation or death, including sex as a baseline covariate:

```
pbc2.id$status2 <- as.numeric(pbc2.id$status != 'alive')
CoxFit <- coxph(Surv(years, status2) ~ sex, data = pbc2.id)
```

Our aim is to assess the strength of the association between the risk
of the composite event and the levels of serum bilirubin that has been
collected during follow-up. We will describe the patient-specific
profiles over time for this biomarker using a linear mixed model, with
fixed-effects, time, sex, and their interaction, and as random effects
random intercepts and random slopes. The syntax to fit this model with
`lme()`

is:

`fm1 <- lme(log(serBilir) ~ year * sex, data = pbc2, random = ~ year | id)`

The joint model that links the survival and longitudinal submodels is
fitted with the following call to the `jm()`

function:

```
jointFit1 <- jm(CoxFit, fm1, time_var = "year")
summary(jointFit1)
#>
#> Call:
#> jm(Surv_object = CoxFit, Mixed_objects = fm1, time_var = "year")
#>
#> Data Descriptives:
#> Number of Groups: 312 Number of events: 169 (54.2%)
#> Number of Observations:
#> log(serBilir): 1945
#>
#> DIC WAIC LPML
#> marginal 4353.451 5259.707 -3116.756
#> conditional 3535.508 3355.363 -1913.528
#>
#> Random-effects covariance matrix:
#>
#> StdDev Corr
#> (Intr) 1.0014 (Intr)
#> year 0.1822 0.4102
#>
#> Survival Outcome:
#> Mean StDev 2.5% 97.5% P Rhat
#> sexfemale -0.1698 0.2695 -0.6668 0.3911 0.5127 1.0043
#> value(log(serBilir)) 1.2415 0.0881 1.0722 1.4207 0.0000 1.0192
#>
#> Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
#> Mean StDev 2.5% 97.5% P Rhat
#> (Intercept) 0.7235 0.1695 0.3902 1.0564 0.0000 1.0000
#> year 0.2660 0.0374 0.1937 0.3414 0.0000 1.0021
#> sexfemale -0.2629 0.1796 -0.6149 0.0878 0.1482 0.9998
#> year:sexfemale -0.0875 0.0401 -0.1663 -0.0088 0.0296 1.0027
#> sigma 0.3470 0.0066 0.3344 0.3601 0.0000 1.0102
#>
#> MCMC summary:
#> chains: 3
#> iterations per chain: 3500
#> burn-in per chain: 500
#> thinning: 1
#> time: 16 sec
```

The output of the `summary()`

method provides some
descriptive statistics of the sample at hand, followed by some fit
statistics based on the marginal (random effects are integrated out
using the Laplace approximation) and conditional on the random effects
log-likelihood functions, followed by the estimated variance-covariance
matrix for the random effects, followed by the estimates for the
survival submodel, followed by the estimates for the longitudinal
submodel(s), and finally some information for the MCMC fitting
algorithm.

By default, `jm()`

adds the subject-specific linear
predictor of the mixed model as a time-varying covariate in the survival
relative risk model. In the output this is named as
`value(log(serBilir))`

to denote that, by default, the
current value functional form is used. That is, we assume that the
instantaneous risk of an event at a specific time \(t\) is associated with the value of the
linear predictor of the longitudinal outcome at the same time point
\(t\).

Standard MCMC diagnostics are available to evaluate convergence. For
example, the traceplot for the association coefficient
`value(log(serBilir))`

is produced with the following
syntax:

`ggtraceplot(jointFit1, "alphas")`

and the density plot with the call:

`ggdensityplot(jointFit1, "alphas")`

To fit a joint model with multiple longitudinal outcomes, we simply
provide a list of mixed models as the second argument of
`jm()`

. In the following example, we extend the joint model
we fitted above by also including the prothrombin time and the log odds
of the presence or not of ascites as time-varying covariates in the
relative risk model for the composite event. Because this model is more
complex, we increase the number of MCMC iterations, the number of
burn-in iterations and the thinning per chain, using the corresponding
control arguments:

```
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)
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: 169 (54.2%)
#> Number of Observations:
#> log(serBilir): 1945
#> prothrombin: 1945
#> ascites: 1885
#>
#> DIC WAIC LPML
#> marginal 11486.28 14328.44 -7488.443
#> conditional 12588.69 12331.25 -6678.692
#>
#> Random-effects covariance matrix:
#>
#> StdDev Corr
#> (Intr) 0.9899 (Intr) year (Intr) year (Intr)
#> year 0.1908 0.4363
#> (Intr) 0.7784 0.5624 0.4824
#> year 0.3329 0.4163 0.3647 0.0526
#> (Intr) 2.9882 0.5596 0.5294 0.6014 0.2829
#> year 0.4737 0.4288 0.6581 0.2332 0.4775 -0.0170
#>
#> Survival Outcome:
#> Mean StDev 2.5% 97.5% P Rhat
#> sexfemale -0.6673 0.3512 -1.3667 0.0108 0.0547 1.0041
#> value(log(serBilir)) 0.5086 0.1703 0.1587 0.8298 0.0010 1.0266
#> value(prothrombin) 0.0170 0.1164 -0.2287 0.2316 0.8420 1.0754
#> value(ascites) 0.5712 0.1226 0.3418 0.8255 0.0000 1.0881
#>
#> Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
#> Mean StDev 2.5% 97.5% P Rhat
#> (Intercept) 0.7052 0.1633 0.3772 1.0198 0.0003 1.0007
#> year 0.2769 0.0334 0.2124 0.3421 0.0000 1.0093
#> sexfemale -0.2474 0.1731 -0.5829 0.0924 0.1597 1.0012
#> year:sexfemale -0.0838 0.0351 -0.1520 -0.0152 0.0187 1.0062
#> sigma 0.3467 0.0065 0.3344 0.3598 0.0000 0.9997
#>
#> Longitudinal Outcome: prothrombin (family = gaussian, link = identity)
#> Mean StDev 2.5% 97.5% P Rhat
#> (Intercept) 10.9461 0.1683 10.6171 11.2746 0.0000 1.0041
#> year 0.2500 0.0731 0.1062 0.3927 0.0007 1.0089
#> sexfemale -0.4065 0.1797 -0.7588 -0.0502 0.0257 1.0025
#> year:sexfemale 0.0516 0.0770 -0.0994 0.2029 0.5137 1.0018
#> sigma 1.0549 0.0202 1.0152 1.0944 0.0000 1.0084
#>
#> Longitudinal Outcome: ascites (family = binomial, link = logit)
#> Mean StDev 2.5% 97.5% P Rhat
#> (Intercept) -4.5582 0.6825 -5.9670 -3.2900 0.0000 1.0803
#> year 0.6363 0.0729 0.5082 0.7917 0.0000 1.1819
#> sexfemale -0.4823 0.6556 -1.7851 0.7774 0.4583 1.0090
#>
#> MCMC summary:
#> chains: 3
#> iterations per chain: 12000
#> burn-in per chain: 2000
#> thinning: 5
#> time: 1.9 min
```

The output for the survival submodel contains now the estimated
coefficients for `value(prothrombin)`

and
`value(ascites)`

, and parameter estimates for all three
longitudinal submodels.

As mentioned above, the default call to `jm()`

includes
the subject-specific linear predictors of the mixed-effects models as
time-varying covariates in the relative risk model. However, this is
just one of the many possibilities we have to link the longitudinal and
survival outcomes. The argument `functional_forms`

of
`jm()`

provides additional options. Based on previous
experience, two extra functional forms are provided, namely, the
time-varying slope and the time-varying *normalized*
area/cumulative-effect. The time-varying slope is the first order
derivative of the subject-specific linear predictor of the mixed-effect
model with respect to the (follow-up) time variable. The time-varying
*normalized* area/cumulative-effect is the integral of the
subject-specific linear predictor of the mixed-effect model from zero to
the current (follow-up) time \(t\)
divided by \(t\). The integral is the
area under the subject-specific longitudinal profile; by dividing the
integral by \(t\) we obtain the average
of the subject-specific longitudinal profile over the corresponding
period \((0, t)\).

To illustrate how the `functional_forms`

argument can be
used to specify these functional forms, we update the joint model
`jointFit2`

by including the time-varying slope of log serum
bilirubin instead of the value, and also the interaction of this slope
with sex, and for prothrombin we include the normalized cumulative
effect. For ascites, we keep the value functional form. The
corresponding syntax to fit this model is:

```
fForms <- list(
"log(serBilir)" = ~ slope(log(serBilir)) + slope(log(serBilir)):sex,
"prothrombin" = ~ area(prothrombin)
)
jointFit3 <- update(jointFit2, functional_forms = fForms)
summary(jointFit3)
#>
#> Call:
#> jm(Surv_object = CoxFit, Mixed_objects = list(fm1, fm2, fm3),
#> time_var = "year", functional_forms = fForms, n_iter = 12000L,
#> n_burnin = 2000L, n_thin = 5L)
#>
#> Data Descriptives:
#> Number of Groups: 312 Number of events: 169 (54.2%)
#> Number of Observations:
#> log(serBilir): 1945
#> prothrombin: 1945
#> ascites: 1885
#>
#> DIC WAIC LPML
#> marginal 11508.40 12656.39 -6945.466
#> conditional 12434.01 12209.69 -6622.497
#>
#> Random-effects covariance matrix:
#>
#> StdDev Corr
#> (Intr) 0.9810 (Intr) year (Intr) year (Intr)
#> year 0.1938 0.4614
#> (Intr) 0.7721 0.5605 0.4947
#> year 0.3337 0.4291 0.3760 0.0721
#> (Intr) 2.8525 0.6004 0.5173 0.6168 0.2777
#> year 0.4760 0.4828 0.7075 0.2792 0.5148 0.0565
#>
#> Survival Outcome:
#> Mean StDev 2.5% 97.5% P Rhat
#> sexfemale -0.0470 0.8043 -1.5439 1.6300 0.9350 1.1080
#> slope(log(serBilir)) 4.1270 2.1588 0.0073 8.6628 0.0497 1.0880
#> slope(log(serBilir)):sexfemale -2.6037 2.3269 -7.5388 1.6883 0.2610 1.1260
#> area(prothrombin) -0.0678 0.2993 -0.7745 0.4262 0.9253 1.1826
#> value(ascites) 0.8323 0.1950 0.4994 1.2811 0.0000 1.3058
#>
#> Longitudinal Outcome: log(serBilir) (family = gaussian, link = identity)
#> Mean StDev 2.5% 97.5% P Rhat
#> (Intercept) 0.6713 0.1612 0.3542 0.9790 0.0000 1.0034
#> year 0.2722 0.0329 0.2074 0.3370 0.0000 1.0007
#> sexfemale -0.2115 0.1708 -0.5457 0.1231 0.2227 1.0046
#> year:sexfemale -0.0732 0.0338 -0.1382 -0.0061 0.0333 1.0016
#> sigma 0.3471 0.0067 0.3343 0.3607 0.0000 1.0028
#>
#> Longitudinal Outcome: prothrombin (family = gaussian, link = identity)
#> Mean StDev 2.5% 97.5% P Rhat
#> (Intercept) 10.9269 0.1676 10.6019 11.2533 0.0000 1.0013
#> year 0.2397 0.0739 0.1005 0.3882 0.0020 1.0042
#> sexfemale -0.3852 0.1766 -0.7338 -0.0423 0.0280 1.0023
#> year:sexfemale 0.0666 0.0765 -0.0843 0.2153 0.3837 1.0041
#> sigma 1.0552 0.0206 1.0154 1.0963 0.0000 1.0026
#>
#> Longitudinal Outcome: ascites (family = binomial, link = logit)
#> Mean StDev 2.5% 97.5% P Rhat
#> (Intercept) -4.5810 0.7075 -6.1089 -3.3105 0.0000 1.0272
#> year 0.6508 0.0856 0.5182 0.8492 0.0000 1.1133
#> sexfemale -0.3335 0.6528 -1.6229 1.0001 0.6017 1.0015
#>
#> MCMC summary:
#> chains: 3
#> iterations per chain: 12000
#> burn-in per chain: 2000
#> thinning: 5
#> time: 2 min
```

As seen above, the `functional_forms`

argument is a named
list with elements corresponding to the longitudinal outcomes. If a
longitudinal outcome is not specified in this list, then the default,
value functional form, is used for that outcome. Each element of the
list should be a one-sided R formula in which the functions
`value()`

, `slope()`

and `area()`

can
be used. Interaction terms between the functional forms and other
(baseline) covariates are also allowed.

When multiple longitudinal outcomes are considered with possibly
different functional forms per outcome, we require to fit a relative
risk model containing several terms. Moreover, it is often of scientific
interest to select which terms/functional-forms per longitudinal outcome
are more strongly associated with the risk of the event of interest. To
facilitate this selection, `jm()`

provides the option to
penalize the regression coefficients using shrinkage priors. As an
example, we refit `jointFit3`

by assuming a Horseshoe prior
for the `alphas`

coefficients (i.e., the coefficients of the
longitudinal outcomes in the relative risk model):

```
jointFit4 <- update(jointFit3, priors = list("penalty_alphas" = "horseshoe"))
cbind("un-penalized" = unlist(coef(jointFit3)),
"penalized" = unlist(coef(jointFit4)))
#> un-penalized penalized
#> gammas.Mean -0.04701237 -0.57666619
#> association.slope(log(serBilir)) 4.12697941 2.60645712
#> association.slope(log(serBilir)):sexfemale -2.60368357 -0.63050964
#> association.area(prothrombin) -0.06778051 0.00798849
#> association.value(ascites) 0.83234894 0.70108773
```

Apart from the Horseshoe prior, the ridge prior is also provided.