Dataset used in the examples
We use the Mayo Clinic’s primary biliary cirrhosis (PBC) dataset in
the examples below. The dataset contains 312 individuals with primary
biliary cirrhosis who participated in a randomised placebo controlled
trial of D-penicillamine conducted at the Mayo Clinic between 1974 and
1984 [19]. However, to ensure the examples run quickly, we use a small
random subset of just 40 patients from the full data.
These example data are contained in two separate data frames. The
first data frame contains multiple-row per patient longitudinal
biomarker information, as shown in
id age sex trt year logBili albumin platelet
1 1 58.76523 f 1 0.0000000 2.67414865 2.60 190
2 1 58.76523 f 1 0.5256674 3.05870707 2.94 183
3 2 56.44627 f 1 0.0000000 0.09531018 4.14 221
4 2 56.44627 f 1 0.4982888 -0.22314355 3.60 188
5 2 56.44627 f 1 0.9993155 0.00000000 3.55 161
6 2 56.44627 f 1 2.1026694 0.64185389 3.92 122
while the second data frame contains single-row per patient survival
information, as shown in
id age sex trt futimeYears status death
1 1 58.76523 f 1 1.095140 2 1
3 2 56.44627 f 1 14.151951 0 0
12 3 70.07255 m 1 2.770705 2 1
16 4 54.74059 f 1 5.270363 2 1
23 5 38.10541 f 0 4.120465 1 0
29 6 66.25873 f 0 6.852841 2 1
The variables included across the two datasets can be defined as
follows:
age in yearsalbumin serum albumin (g/dl)logBili logarithm of serum bilirubindeath indicator of death at endpointfutimeYears time (in years) between baseline and the
earliest of death, transplantion or censoringid numeric ID unique to each individualplatelet platelet countsex gender (m = male, f = female)status status at endpoint (0 = censored, 1 =
transplant, 2 = dead)trt binary treatment code (0 = placebo, 1 =
D-penicillamine)year time (in years) of the longitudinal measurements,
taken as time since baseline)
A description of the example datasets can be found by accessing the
following help documentation:
help("datasets", package = "rstanarm")
Fitting the models
Univariate joint model (current value association structure)
In this example we fit a simple univariate joint model, with one
normally distributed longitudinal marker, an association structure based
on the current value of the linear predictor, and B-splines baseline
hazard. To fit the model we use the joint (longitudinal and
time-to-event) modelling function in the rstanarm
package: stan_jm. When calling stan_jm we
must, at a minimum, specify a formula object for each of the
longitudinal and event submodels (through the arguments
formulaLong and formulaEvent), the data frames
which contain the variables for each of the the longitudinal and event
submodels (through the arguments dataLong and
dataEvent), and the name of the variable representing time
in the longitudinal submodel (through the argument
time_var).
The formula for the longitudinal submodel is specified using the
lme4 package formula style. That is
y ~ x + (random_effects | grouping_factor). In this example
we specify that log serum bilirubin (logBili) follows a
subject-specific linear trajectory. To do this we include a fixed
intercept and fixed slope (year), as well as a random
intercept and random slope for each subject id
((year | id)).
The formula for the event submodel is specified using the
survival package formula style. That is, the outcome of
the left of the ~ needs to be of the format
Surv(event_time, event_indicator) for single row per
individual data, or
Surv(start_time, stop_time, event_indicator) for multiple
row per individual data. The latter allows for exogenous time-varying
covariates to be included in the event submodel. In this example we
assume that the log hazard of death is linearly related to gender
(sex) and an indicator of treatment with D-penicillamine
(trt).
library(rstanarm)
mod1 <- stan_jm(formulaLong = logBili ~ sex + trt + year + (year | id),
dataLong = pbcLong,
formulaEvent = survival::Surv(futimeYears, death) ~ sex + trt,
dataEvent = pbcSurv,
time_var = "year",
chains = 1, refresh = 2000, seed = 12345)
Fitting a univariate joint model.
Please note the warmup may be much slower than later iterations!
SAMPLING FOR MODEL 'jm' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000308 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 3.08 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 21.097 seconds (Warm-up)
Chain 1: 25.727 seconds (Sampling)
Chain 1: 46.824 seconds (Total)
Chain 1:
The argument refresh = 2000 was specified so that Stan
didn’t provide us with excessive progress updates whilst fitting the
model. However, if you are fitting a model that will take several
minutes or hours to fit, then you may wish to request progress updates
quite regularly, for example setting refresh = 20 for every
20 iterations (by default the refresh argument is set to 1/10th of the
total number of iterations).
The fitted model is returned as an object of the S3 class
stanjm. We have a variety of methods and post-estimation
functions available for this class, including: print,
summary, plot, fixef,
ranef, coef, VarCorr,
posterior_interval, update, and more. Here, we
will examine the most basic output for the fitted joint model by typing
print(mod1):
stan_jm
formula (Long1): logBili ~ sex + trt + year + (year | id)
family (Long1): gaussian [identity]
formula (Event): survival::Surv(futimeYears, death) ~ sex + trt
baseline hazard: bs
assoc: etavalue (Long1)
------
Longitudinal submodel: logBili
Median MAD_SD
(Intercept) 0.304 0.558
sexf 0.477 0.523
trt -0.127 0.370
year 0.212 0.043
sigma 0.354 0.016
Event submodel:
Median MAD_SD exp(Median)
(Intercept) -3.179 0.621 0.042
sexf -0.336 0.620 0.715
trt -0.720 0.415 0.487
Long1|etavalue 1.341 0.238 3.821
b-splines-coef1 -0.836 1.035 NA
b-splines-coef2 0.518 0.881 NA
b-splines-coef3 -1.797 1.196 NA
b-splines-coef4 0.299 1.652 NA
b-splines-coef5 -0.061 1.665 NA
b-splines-coef6 -0.820 1.665 NA
Group-level error terms:
Groups Name Std.Dev. Corr
id Long1|(Intercept) 1.2954
Long1|year 0.1921 0.52
Num. levels: id 40
Sample avg. posterior predictive distribution
of longitudinal outcomes:
Median MAD_SD
Long1|mean_PPD 0.585 0.032
------
For info on the priors used see help('prior_summary.stanreg').
The “Long1|etavalue” row under “Event submodel” is our \(\alpha_{mq}\) parameter (\(m = 1\), \(q =
1\)). The estimated median of tells us that for each one unit
increase in an individual’s underlying level of log serum bilirubin,
their estimated log hazard of death increases by some amount. The mean
absolute deviation (MAD) is provided as a more robust estimate of the
standard deviation of the posterior distribution. In this case the
MAD_SD for the association parameter indicates there is quite large
uncertainty around the estimated association between log serum bilirubin
and risk of death (recall this is a small dataset).
If we wanted some slightly more detailed output for each of the model
parameters, as well as further details regarding the model estimation
(for example computation time, number of longitudinal observations,
number of individuals, type of baseline hazard, etc) we can instead use
the summary method:
summary(mod1, probs = c(.025,.975))
Model Info:
function: stan_jm
formula (Long1): logBili ~ sex + trt + year + (year | id)
family (Long1): gaussian [identity]
formula (Event): survival::Surv(futimeYears, death) ~ sex + trt
baseline hazard: bs
assoc: etavalue (Long1)
algorithm: sampling
priors: see help('prior_summary')
sample: 1000 (posterior sample size)
num obs: 304 (Long1)
num subjects: 40
num events: 29 (72.5%)
groups: id (40)
runtime: 0.8 mins
Estimates:
mean sd 2.5%
Long1|(Intercept) 0.294 0.586 -0.926
Long1|sexf 0.499 0.560 -0.529
Long1|trt -0.144 0.386 -0.973
Long1|year 0.214 0.042 0.136
Long1|sigma 0.354 0.017 0.324
Long1|mean_PPD 0.586 0.030 0.527
Event|(Intercept) -3.191 0.622 -4.477
Event|sexf -0.299 0.620 -1.407
Event|trt -0.740 0.445 -1.663
Event|b-splines-coef1 -0.949 1.065 -3.263
Event|b-splines-coef2 0.462 0.908 -1.410
Event|b-splines-coef3 -1.807 1.186 -4.100
Event|b-splines-coef4 0.388 1.567 -2.682
Event|b-splines-coef5 -0.077 1.647 -3.410
Event|b-splines-coef6 -1.014 1.692 -4.884
Assoc|Long1|etavalue 1.347 0.237 0.915
Sigma[id:Long1|(Intercept),Long1|(Intercept)] 1.678 0.438 1.012
Sigma[id:Long1|year,Long1|(Intercept)] 0.129 0.076 0.018
Sigma[id:Long1|year,Long1|year] 0.037 0.020 0.014
log-posterior -329.512 9.997 -349.426
97.5%
Long1|(Intercept) 1.367
Long1|sexf 1.627
Long1|trt 0.582
Long1|year 0.299
Long1|sigma 0.387
Long1|mean_PPD 0.647
Event|(Intercept) -2.068
Event|sexf 1.036
Event|trt 0.128
Event|b-splines-coef1 0.882
Event|b-splines-coef2 2.114
Event|b-splines-coef3 0.370
Event|b-splines-coef4 3.494
Event|b-splines-coef5 2.989
Event|b-splines-coef6 1.810
Assoc|Long1|etavalue 1.829
Sigma[id:Long1|(Intercept),Long1|(Intercept)] 2.703
Sigma[id:Long1|year,Long1|(Intercept)] 0.300
Sigma[id:Long1|year,Long1|year] 0.085
log-posterior -311.470
Diagnostics:
mcse Rhat n_eff
Long1|(Intercept) 0.029 0.999 411
Long1|sexf 0.028 0.999 407
Long1|trt 0.023 1.001 281
Long1|year 0.002 0.999 288
Long1|sigma 0.001 1.002 953
Long1|mean_PPD 0.001 1.001 1165
Event|(Intercept) 0.019 1.001 1131
Event|sexf 0.018 1.002 1156
Event|trt 0.013 0.999 1149
Event|b-splines-coef1 0.040 0.999 703
Event|b-splines-coef2 0.029 1.000 1008
Event|b-splines-coef3 0.046 0.999 660
Event|b-splines-coef4 0.065 0.999 589
Event|b-splines-coef5 0.064 0.999 654
Event|b-splines-coef6 0.061 0.999 767
Assoc|Long1|etavalue 0.007 0.999 1074
Sigma[id:Long1|(Intercept),Long1|(Intercept)] 0.025 0.999 306
Sigma[id:Long1|year,Long1|(Intercept)] 0.006 1.001 182
Sigma[id:Long1|year,Long1|year] 0.001 1.003 214
log-posterior 0.708 1.002 199
For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
The easiest way to extract the correlation matrix for the random
effects (aside from viewing the print output) is to use the
VarCorr function (modelled on the VarCorr
function from the lme4 package). If you wish to extract
the variances and covariances (instead of the standard deviations and
correlations) then you can type the following to return a data frame
with all of the relevant information:
as.data.frame(VarCorr(mod1))
grp var1 var2 vcov sdcor
1 id Long1|(Intercept) <NA> 1.67802152 1.2953847
2 id Long1|year <NA> 0.03691896 0.1921431
3 id Long1|(Intercept) Long1|year 0.12937327 0.5197818
Univariate joint model (current value and current slope association
structure)
In the previous example we were fitting a shared parameter joint
model which assumed that the log hazard of the event (in this case the
log hazard of death) at time t was linearly related to the
subject-specific expected value of the longitudinal marker (in this case
the expected value of log serum bilirubin) also at time t. This
is the default association structure, although it could be explicitly
specified by setting the assoc = "etavalue" argument.
However, let’s suppose we believe that the log hazard of death is
actually related to both the current value of log serum
bilirubin and the current rate of change in log serum
bilirubin. To estimate this joint model we need to indicate that we want
to also include the subject-specific slope (at time t) from the
longitudinal submodel as part of the association structure. We do this
by setting the assoc argument equal to a character vector
c("etavalue", "etaslope") which indicates our desired
association structure:
mod2 <- stan_jm(formulaLong = logBili ~ sex + trt + year + (year | id),
dataLong = pbcLong,
formulaEvent = survival::Surv(futimeYears, death) ~ sex + trt,
dataEvent = pbcSurv,
assoc = c("etavalue", "etaslope"),
time_var = "year",
chains = 1, refresh = 2000, seed = 12345)
In this example the subject-specific slope is actually constant
across time t since we have a linear trajectory. Note however
that we could still use the "etaslope" association
structure even if we had a non-linear subject specific trajectory (for
example modelled using cubic splines or polynomials).
Multivariate joint model (current value association structures)
Suppose instead that we were interested in two repeatedly
measured clinical biomarkers, log serum bilirubin and serum albumin, and
their association with the risk of death. We may wish to model these two
biomarkers, allowing for the correlation between them, and estimating
their respective associations with the log hazard of death. We will fit
a linear mixed effects submodel (identity link, normal distribution) for
each biomarker with a patient-specific intercept and linear slope but no
other covariates. In the event submodel we will include gender
(sex) and treatment (trt) as baseline
covariates. Each biomarker is assumed to be associated with the log
hazard of death at time \(t\) via it’s
expected value at time \(t\) (i.e. a
current value association structure).
The model we are going to fit can therefore be specified as:
\[
y_{im}(t_{ijm}) \sim N(\mu_{im}(t_{ijm}), \sigma_m)
\]
\[
\eta_{im}(t) = \mu_{im}(t) = \beta_{0m} + \beta_{1m} t + b_{0mi} +
b_{1mi} t
\]
\[
h_i(t) = h_0(t; \boldsymbol{\omega}) \exp(\gamma_1 w_{1i} + \gamma_2
w_{2i} + \alpha_{1i} \mu_{i1}(t) + \alpha_{2i} \mu_{i2}(t))
\]
where \(t\) is time in years, and
\(w_{1i}\) and \(w_{2i}\) are, respectively, the gender and
treatment indicators for individual \(i\).
(Note that due to the very small sample size, the clinical findings
from this analysis should not to be overinterpreted!).
mod3 <- stan_jm(
formulaLong = list(
logBili ~ sex + trt + year + (year | id),
albumin ~ sex + trt + year + (year | id)),
formulaEvent = survival::Surv(futimeYears, death) ~ sex + trt,
dataLong = pbcLong, dataEvent = pbcSurv,
time_var = "year",
chains = 1, refresh = 2000, seed = 12345)
Fitting a multivariate joint model.
Please note the warmup may be much slower than later iterations!
SAMPLING FOR MODEL 'jm' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.000368 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 3.68 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 45.961 seconds (Warm-up)
Chain 1: 43.684 seconds (Sampling)
Chain 1: 89.645 seconds (Total)
Chain 1:
We can now examine the output from the fitted model, for
example
stan_jm
formula (Long1): logBili ~ sex + trt + year + (year | id)
family (Long1): gaussian [identity]
formula (Long2): albumin ~ sex + trt + year + (year | id)
family (Long2): gaussian [identity]
formula (Event): survival::Surv(futimeYears, death) ~ sex + trt
baseline hazard: bs
assoc: etavalue (Long1), etavalue (Long2)
------
Longitudinal submodel 1: logBili
Median MAD_SD
(Intercept) 0.166 0.576
sexf 0.585 0.561
trt -0.050 0.370
year 0.218 0.043
sigma 0.355 0.017
Longitudinal submodel 2: albumin
Median MAD_SD
(Intercept) 3.497 0.267
sexf 0.027 0.237
trt -0.008 0.183
year -0.154 0.023
sigma 0.291 0.014
Event submodel:
Median MAD_SD exp(Median)
(Intercept) 6.838 2.753 932.668
sexf -0.043 0.691 0.958
trt -0.494 0.450 0.610
Long1|etavalue 0.801 0.299 2.228
Long2|etavalue -3.081 0.881 0.046
b-splines-coef1 -1.064 1.126 NA
b-splines-coef2 0.502 0.896 NA
b-splines-coef3 -2.688 1.388 NA
b-splines-coef4 -0.481 1.808 NA
b-splines-coef5 -1.288 1.808 NA
b-splines-coef6 -2.691 1.942 NA
Group-level error terms:
Groups Name Std.Dev. Corr
id Long1|(Intercept) 1.25984
Long1|year 0.18631 0.48
Long2|(Intercept) 0.52341 -0.65 -0.49
Long2|year 0.09583 -0.56 -0.80 0.44
Num. levels: id 40
Sample avg. posterior predictive distribution
of longitudinal outcomes:
Median MAD_SD
Long1|mean_PPD 0.588 0.029
Long2|mean_PPD 3.343 0.024
------
For info on the priors used see help('prior_summary.stanreg').
or we can examine the summary output for the association parameters
alone:
summary(mod3, pars = "assoc")
Model Info:
function: stan_jm
formula (Long1): logBili ~ sex + trt + year + (year | id)
family (Long1): gaussian [identity]
formula (Long2): albumin ~ sex + trt + year + (year | id)
family (Long2): gaussian [identity]
formula (Event): survival::Surv(futimeYears, death) ~ sex + trt
baseline hazard: bs
assoc: etavalue (Long1), etavalue (Long2)
algorithm: sampling
priors: see help('prior_summary')
sample: 1000 (posterior sample size)
num obs: 304 (Long1), 304 (Long2)
num subjects: 40
num events: 29 (72.5%)
groups: id (40)
runtime: 1.5 mins
Estimates:
mean sd 2.5% 25% 50% 75% 97.5%
Assoc|Long1|etavalue 0.806 0.301 0.221 0.602 0.801 1.010 1.416
Assoc|Long2|etavalue -3.157 0.939 -5.258 -3.737 -3.081 -2.515 -1.546
Diagnostics:
mcse Rhat n_eff
Assoc|Long1|etavalue 0.008 1.001 1514
Assoc|Long2|etavalue 0.033 1.001 815
For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
Posterior predictions
We can also access the range of post-estimation functions (described
in the stan_jm and related help documentation; see for
example help(posterior_traj) or
help(posterior_survfit)).
Predicted individual-specific longitudinal trajectory for in-sample
individuals
Predicted individual-specific biomarker values can be obtained using
either the posterior_traj or posterior_predict
function. The posterior_traj is preferable, because it can
be used to obtain the biomarker values at a series of evenly spaced time
points between baseline and the individual’s event or censoring time by
using the default interpolate = TRUE option. Whereas, the
posterior_predict function only provides the predicted
biomarker values at the observed time points, or the time points in the
new data. Predicting the biomarker values at a series of evenly spaced
time points can be convenient because they can be easily used for
plotting the longitudinal trajectory. Moreover, by default the
posterior_traj returns a data frame with variables
corresponding to the individual ID, the time, the predicted mean
biomarker value, the limits for the 95% credible interval
(i.e. uncertainty interval for the predicted mean biomarker value), and
limits for the 95% prediction interval (i.e. uncertainty interval for a
predicted biomarker data point), where the level for the uncertainty
intervals can be changed via the prob argument. Conversely,
the posterior_predict function returns an \(S\) by \(N\) matrix of predictions where \(S\) is the number of posterior draws and
\(N\) is the number of prediction time
points (note that this return type can also be obtained for
posterior_traj by specifying the argument
return_matrix = TRUE).
As an example, let’s plot the predicted individual-specific
longitudinal trajectories for each of the two biomarkers (log serum
bilirubin and serum albumin) in the multivariate joint model estimated
above. We will do this for three individuals (IDs 6, 7 and 8) who were
included in the model estimation.
Here are the plots for log serum bilirubin:
p1 <- posterior_traj(mod3, m = 1, ids = 6:8)
pp1 <- plot(p1, plot_observed = TRUE)
pp1
and here are the plots for serum albumin:
p2 <- posterior_traj(mod3, m = 2, ids = 6:8)
pp2 <- plot(p2, plot_observed = TRUE)
pp2
The m argument specifies which biomarker we want to
predict for (only relevant for a multivariate joint model). The
ids argument is optional, and specifies a subset of
individuals for whom we want to predict. In the plotting method, the
plot_observed = TRUE specifies that we want to include the
observed biomarker values in the plot of the longitudinal
trajectory.
If we wanted to extrapolate the trajectory forward from the event or
censoring time for each individual, then this can be easily achieved by
specifying extrapolate = TRUE in the
posterior_traj call. For example, here is the plot for log
serum bilirubin with extrapolation:
p3 <- posterior_traj(mod3, m = 1, ids = 6:8, extrapolate = TRUE)
pp3 <- plot(p3, plot_observed = TRUE, vline = TRUE)
pp3
and for serum albumin with extrapolation:
p4 <- posterior_traj(mod3, m = 2, ids = 6:8, extrapolate = TRUE)
pp4 <- plot(p4, plot_observed = TRUE, vline = TRUE)
pp4
Here, we included the vline = TRUE which adds a vertical
dashed line at the timing of the individual’s event or censoring time.
The interpolation and extrapolation of the biomarker trajectory can be
further controlled through the control argument to the
posterior_traj function; for example, we could specify the
number of time points at which to predict, the distance by which to
extrapolate, and so on.
We could customize these plots further, for example, by using any of
the ggplot2 functionality or using the additional
arguments described in help(plot.predict.stanjm).
Predicted individual-specific survival curves for in-sample
individuals
Predicted individual-specific survival probabilities and/or survival
curves can be obtained using the posterior_survfit
function. The function by default returns a data frame with the
individual ID, the time, and the predicted survival probability
(posterior mean and limits for the 95% credible interval). The
uncertainty level for the credible interval can be changed via the
prob argument. By default, individual-specific survival
probabilities are calculated conditional on the individual’s
last known survival time. When we are predicting survival probabilities
for individuals that were used in the estimation of the model
(i.e. in-sample individuals, where no new covariate data is provided),
then the individual’s “last known survival time” will be their event or
censoring time. (Note that if we wanted didn’t want to condition on the
individual’s last known survival time, then we could specify
condition = FALSE, but we probably wouldn’t want to do this
unless we were calculating marginal or standardised survival
probabilities, which are discussed later).
The default argument extrapolate = TRUE specifies that
the individual-specific conditional survival probabilities will be
calculated at evenly spaced time points between the individual’s last
known survival time and the maximum follow up time that was observed in
the estimation sample. The behaviour of the extrapolation can be further
controlled via the control argument. If we were to specify
extrapolate = FALSE then the survival probabilities would
only be calculated at one time point, which could be specified in the
times argument (or otherwise would default to the
individual’s last known survival time).
As an example, let plot the predicted individual-specific conditional
survival curve for the same three individual’s that were used in the
previous example. The predicted survival curve will be obtained under
the multivariate joint model estimated above.
p5 <- posterior_survfit(mod3, ids = 6:8)
pp5 <- plot(p5)
pp5
We could customize the plot further, for example, by using any of the
ggplot2 functionality or using the additional arguments
described in help(plot.survfit.stanjm).
Combined plot of longitudinal trajectories and survival curves
The package also provides a convenience plotting function, which
combines plots of the individual-specific longitudinal trajectories, and
the individual-specific survival function. We can demonstrate this by
replotting the predictions for the three individuals in the previous
example:
plot_stack_jm(yplot = list(pp3, pp4), survplot = pp5)
Here we can see the strong relationship between the underlying values
of the biomarkers and mortality. Patient 8 who, relative to
patients 6 and 7, has a higher underlying
value for log serum bilirubin and a lower underlying value for serum
albumin at the end of their follow up has a far worse predicted
probability of survival.
Predicted individual-specific longitudinal trajectory and survival
curve for out-of-sample individuals (i.e. dynamic predictions)
Let us take an individual from our training data, in this case the
individual with subject ID value 8. However, we will
pretend this individual was not a member of our training data and rather
that they are a new individual for whom we have obtained new biomarker
measurements. Our goal is to obtain predictions for the longitudinal
trajectory for this individual, and their conditional survival curve,
given that we know they are conditional on their biomarker measurements
we currently have available.
First, let’s extract the data for subject 8 and then
rename their subject ID value so that they appear to be an individual
who was not included in our training dataset:
ndL <- pbcLong[pbcLong$id == 8, , drop = FALSE]
ndE <- pbcSurv[pbcSurv$id == 8, , drop = FALSE]
ndL$id <- paste0("new_patient")
ndE$id <- paste0("new_patient")
Note that we have both the longitudinal data and event data for this
new individual. We require data for both submodels because we are going
to generate dynamic predictions that require drawing new
individual-specific parameters (i.e. random effects) for this individual
conditional on their observed data. That means we need to evaluate the
likelihood for the full joint model and that requires both the
longitudinal and event data (note however that the status indicator
death will be ignored, since it is assumed that the
individual we are predicting for is still alive at the time we wish to
generate the predictions).
Now we can pass this data to the posterior_traj function
in the same way as for the in-sample individuals, except we will now
specify the newdataLong and newdataEvent
arguments. We will also specify the last_time argument so
that the function knows which variable in the event data specifies the
individual’s last known survival time (the default behaviour is to use
the time of the last biomarker measurement).
Our predictions for this new individual for the log serum bilirubin
trajectory can be obtained using:
p6 <- posterior_traj(mod3, m = 1,
newdataLong = ndL,
newdataEvent = ndE,
last_time = "futimeYears")
Drawing new random effects for 1 individuals. Monitoring progress:
| | | 0% | |======================================================================| 100%
pp6 <- plot(p6, plot_observed = TRUE, vline = TRUE)
pp6
and for the serum albumin trajectory:
p7 <- posterior_traj(mod3, m = 2,
newdataLong = ndL,
newdataEvent = ndE,
last_time = "futimeYears")
Drawing new random effects for 1 individuals. Monitoring progress:
| | | 0% | |======================================================================| 100%
pp7 <- plot(p7, plot_observed = TRUE, vline = TRUE)
pp7
For the conditional survival probabilities we use similar
information, provided to the posterior_survfit
function:
p8 <- posterior_survfit(mod3,
newdataLong = ndL,
newdataEvent = ndE,
last_time = "futimeYears")
Drawing new random effects for 1 individuals. Monitoring progress:
| | | 0% | |======================================================================| 100%
We can then use the plot_stack_jm function, as we saw in
a previous example, to stack the plots of the longitudinal trajectory
and the conditional survival curve:
plot_stack_jm(yplot = list(pp6, pp7), survplot = pp8)
Here we see that the predicted longitudinal trajectories and
conditional survival curve for this individual, obtained using the
dynamic predictions approach, are similar to the predictions we obtained
when we used their individual-specific parameters from the original
model estimation. This is because in both situations we are conditioning
on the same outcome data.
Side note: We can even compare the estimated
individual specific parameters obtained under the two approaches. For
example, here is the posterior mean for the estimated
individual-specific parameters for individual 8 from the
fitted model:
c(ranef(mod3)[["Long1"]][["id"]][8,],
ranef(mod3)[["Long2"]][["id"]][8,])
$`(Intercept)`
[1] -1.688937
$year
[1] 0.1197032
$`(Intercept)`
[1] 0.4621096
$year
[1] -0.05095125
and here is the mean of the draws for the individual-specific
parameters for individual 8 under the dynamic predictions
approach:
colMeans(attr(p6, "b_new"))
b[Long1|(Intercept) id:new_patient] b[Long1|year id:new_patient]
-1.69392611 0.11030577
b[Long2|(Intercept) id:new_patient] b[Long2|year id:new_patient]
0.46042200 -0.04324049
Predicted population-level longitudinal trajectory
Suppose we wanted to predict the longitudinal trajectory for each of
the biomarkers, marginalising over the distribution of the
individual-specific parameters. To do this, we can pass a new data frame
with the covariate values we want to use in the predictions. Here, we
will demonstrate this by obtaining the predicted trajectory for log
serum bilirubin, under the multivariate joint model that was estimated
previously. Our prediction data will require the variables
year, sex and trt, since these
were the covariates used in the longitudinal submodel.
We will predict the value of log serum bilirubin at years 0 through
10, for each combination of sex and trt. We
also need to include the id variable in our prediction data
because this is relevant to the longitudinal submodel. Since we want to
marginalise over the individual-specific parameters
(i.e. individual-level random effects) we need to note two things:
First, the values for the id variable must
not match any individual used in the model estimation. Here, we
use the following id values: "male_notrt",
"female_notrt", "male_trt", and
"female_trt", since each individual in our prediction data
represents a different combination of sex and
trt. However, we could have given the individuals any
id value just as long as is didn’t match an individual who
was used in the model estimation
Second, we need to specify the argument
dynamic = FALSE when calling posterior_traj.
This specifies that we do not want to draw new individual-specific
parameters conditional on outcome data observed up to some time \(t\). Instead, we want predictions that
marginalise over the distribution of individual-specific parameters and
are therefore conditional only on covariates and not
conditional on outcome data for the new individuals.
Here is our prediction data:
ndL <- expand.grid(year = seq(0, 10, 1),
sex = c("m", "f"),
trt = 0:1)
ndL$id <- rep(c("male_notrt", "female_notrt",
"male_trt", "female_trt"), each = 11)
ndL <- ndL[, c(4,1,2,3)]
str(ndL)
'data.frame': 44 obs. of 4 variables:
$ id : chr "male_notrt" "male_notrt" "male_notrt" "male_notrt" ...
$ year: num 0 1 2 3 4 5 6 7 8 9 ...
$ sex : Factor w/ 2 levels "m","f": 1 1 1 1 1 1 1 1 1 1 ...
$ trt : int 0 0 0 0 0 0 0 0 0 0 ...
And to predict the marginal longitudinal trajectory for log serum
bilirubin under each covariate profile and plot it we can type:
p1 <- posterior_traj(mod3, m = 1, newdataLong = ndL, dynamic = FALSE)
plot(p1) + ggplot2::coord_cartesian(ylim = c(-10,15))
Because we are marginalising over the distribution of the
individual-specific parameters, we are incorporating all the variation
related to between-individual differences, and therefore the prediction
interval is wide (shown by the shaded area around the marginal
longitudinal trajectory). The magnitude of the effects of both
sex and trt are relatively small compared to
the population-level effect of year and the
between-individual variation in the intercept and slope. For example,
here are the point estimates for the population-level effects of
sex, trt, and year:
(Intercept) sexf trt year
0.16627636 0.58536144 -0.04983978 0.21776053
and here are the standard deviations for the individual-level random
effects:
Groups Name Std.Dev. Corr
id Long1|(Intercept) 1.259840
Long1|year 0.186312 0.479
Long2|(Intercept) 0.523409 -0.655 -0.489
Long2|year 0.095831 -0.559 -0.805 0.444
This shows us that the point estimates for the population-level
effects of sex and trt are 0.57 and -0.10,
respectively, whereas the standard deviation for the individual-specific
intercept and slope parameters are 1.24 and 0.19; hence, any differences
due to the population-level effects of gender and treatment
(i.e. differences in the black line across the four panels of the plot)
are swamped by the width of the uncertainty intervals (i.e. the grey
shaded areas).
Standardised survival curves
In this example we show how a standardised survival curve can be
obtained, where the \(i =
1,...,N^{pred}\) individuals used in generating the standardised
survival curve are the same individuals that were used in estimating the
model. We will obtain the survival curve for the multivariate joint
model estimated in an earlier example (mod3). The
standardise = TRUE argument to
posterior_survfit specifies that we want to obtain
individual-specific predictions of the survival curve and then average
these. Because, in practical terms, we need to obtain survival
probabilities at time \(t\) for each
individual and then average them we want to explicitly specify the
values of \(t\) we want to use (and the
same values of \(t\) will be used for
individuals). We specify the values of \(t\) to use via the times
argument; here we will predict the standardised survival curve at time 0
and then for convenience we can just specify
extrapolate = TRUE (which is the default anyway) which will
mean we automatically predict at 10 evenly spaced time points between 0
and the maximum event or censoring time.
p1 <- posterior_survfit(mod3, standardise = TRUE, times = 0)
head(p1) # data frame with standardised survival probabilities
year survpred ci_lb ci_ub
1 0.0000 1.0000 1.0000 1.0000
2 1.0154 0.8324 0.7721 0.8923
3 2.0307 0.7355 0.6726 0.7901
4 3.0461 0.6855 0.6248 0.7450
5 4.0614 0.6418 0.5838 0.7104
6 5.0768 0.5973 0.5446 0.6535
plot(p1) # plot the standardised survival curve
