BJM fits a backward joint model of multivariate longitudinal
outcomes and time-to-event data, and uses it to make dynamic
predictions: given a patient’s longitudinal history up to some
prediction_time, predict their risk of an event (and, if
desired, their future biomarker values) conditional on survival to that
point.
Fitting and prediction is a four-step pipeline:
survivalSub() fits the marginal survival sub-model
(and, optionally, a competing-risks sub-model).longitudinalSub() fits the longitudinal sub-model(s),
one per biomarker.dynamicPrediction() combines the two fits to predict a
patient’s event risk over a future horizon.dynamicPredictionBio() additionally predicts a future
value of one specific biomarker.This vignette walks through all four steps on the pbc3
data set that ships with the package, then covers two practical
questions that come up once the mechanics are working: how to choose the
bandcount1/ bandcount2/bandcount3
numerical-integration tuning parameters, and what happens when the
pipeline functions are called with malformed inputs.
pbc3 is in long format: one row per longitudinal
measurement, with baseline covariates (age,
sex, …) repeated on every row for a given patient, and
event-time variables (years, status3, …) that
are also constant within patient.
survivalSub()survivalSub() needs one row per patient, so we first
drop the repeated longitudinal rows. form_marginal_surv is
a standard survival::Surv() formula for the overall event.
form_conditional_cr is optional: supply it when there is
more than one competing event type and you want to predict
event-specific risk (here, status4 distinguishes the two
causes among patients who had some event).
data_survival_fitting <- pbc3[!duplicated(pbc3$id), ]
survival_fit_all <- survivalSub(
data_survival_fitting,
form_marginal_surv = Surv(years, status3) ~ age + sex,
form_conditional_cr = status4 ~ years + age + sex
)
survival_fit_all
#>
#> Call:
#> survivalSub(form_marginal_surv = Surv(years, status3) ~ age + sex,
#> form_conditional_cr = status4 ~ years + age + sex)
#>
#> Data Descriptives:
#> Number of subjects : 312
#> Number of events : 169
#> Cause-1 events (CR model) : 29
#> Cause-2 events (CR model) : 140
#>
#> =================================================================
#> Marginal Survival Sub-model [Cox PH]
#> -----------------------------------------------------------------
#> Formula: Surv(years, status3) ~ age + sex
#>
#> Coef exp(Coef) SE z p-value
#> age 0.020411 1.020621 0.007584 2.691 0.00712 **
#> sex -0.497707 0.607923 0.207078 -2.403 0.01624 *
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> n = 312, events = 169
#> Concordance = 0.570 (se = 0.0226)
#> Likelihood ratio = 14.20 on 2 df, p = 0.0008263
#> Wald test = 15.29 on 2 df, p = 0.0004776
#> Score (logrank) = 15.50 on 2 df, p = 0.0004299
#>
#> =================================================================
#> Conditional Competing-Risks Sub-model [Logistic GLM]
#> -----------------------------------------------------------------
#> Formula: status4 ~ years + age + sex
#>
#> Coef SE z p-value
#> (Intercept) 5.65622 1.70478 3.318 0.000907 ***
#> years -0.02898 0.08919 -0.325 0.745242
#> age -0.15305 0.03076 -4.976 6.5e-07 ***
#> sex 0.02109 0.77150 0.027 0.978192
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Null deviance : 154.94 on 168 df
#> Residual deviance : 115.05 on 165 df
#> AIC : 123.05
#> =================================================================longitudinalSub()Each biomarker gets its own fixed-effects formula and random-effects
formula, supplied as same-length, same-order lists.
data_fit_all is a matching list of fitting data sets (one
per biomarker) — or a single data.frame that is reused for
every biomarker, if all biomarkers are fit on the same data. Here we fit
two biomarkers, serBilir and albumin, on the
complete-case subset.
long_sub_fixed <- list(
"serBilir" = serBilir ~ year + age + sex + (years) + (years) * year,
"albumin" = albumin ~ year + age + sex + (years) + (years) * year
)
long_sub_random <- list(
"serBilir" = ~ year | id,
"albumin" = ~ year | id
)
data_fit_all <- list(pbc3[pbc3$status3 == 1, ], pbc3[pbc3$status3 == 1, ])
long_fit_all <- longitudinalSub(data_fit_all, long_sub_fixed, long_sub_random)
long_fit_all
#>
#> Call:
#> longitudinalSub(M = 2 longitudinal outcomes)
#>
#> =================================================================
#> Longitudinal Sub-model(s)
#> -----------------------------------------------------------------
#> [1] serBilir
#> Formula: serBilir ~ year + age + sex + (years) + (years) * year
#>
#> Value SE t p-value
#> (Intercept) 2.424554 0.384604 6.304 5.13e-10 ***
#> year 0.442242 0.042957 10.295 < 2e-16 ***
#> age -0.011675 0.005785 -2.018 0.0439 *
#> sex 0.078449 0.167140 0.469 0.6390
#> years -0.204130 0.021557 -9.469 < 2e-16 ***
#> year:years -0.027852 0.006112 -4.557 6.12e-06 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual std. error (sigma): 0.4186
#> n (subjects) = 169, N (observations) = 872
#> Log-likelihood: -772.87
#> AIC: 1565.74, BIC: 1613.45
#>
#> -----------------------------------------------------------------
#> [2] albumin
#> Formula: albumin ~ year + age + sex + (years) + (years) * year
#>
#> Value SE t p-value
#> (Intercept) 3.578828 0.150261 23.817 < 2e-16 ***
#> year -0.219541 0.026452 -8.300 5.36e-16 ***
#> age -0.006625 0.002261 -2.930 0.00350 **
#> sex -0.170236 0.064981 -2.620 0.00899 **
#> years 0.077939 0.008732 8.926 < 2e-16 ***
#> year:years 0.010955 0.003585 3.056 0.00233 **
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual std. error (sigma): 0.3121
#> n (subjects) = 169, N (observations) = 872
#> Log-likelihood: -386.29
#> AIC: 792.58, BIC: 840.29
#>
#> =================================================================
#> Multivariate Random-Effects Covariance Matrix (D)
#> -----------------------------------------------------------------
#> (Intercept)_1 year_1 (Intercept)_2 year_2
#> (Intercept)_1 0.6012 -0.0231 -0.0473 -0.0070
#> year_1 -0.0231 0.0193 0.0121 -0.0050
#> (Intercept)_2 -0.0473 0.0121 0.0722 -0.0045
#> year_2 -0.0070 -0.0050 -0.0045 0.0050
#> =================================================================dynamicPrediction()To predict for a specific patient, build a
data_predict_all list (one data.frame per
biomarker, mirroring data_fit_all) containing only that
patient’s measurements up to prediction_time —
later measurements would not be available yet in a real prediction
setting. Below we predict, for patient 2, the risk of each event type
within one year of year = 5.
survival_variable_all/survival_trans_function
describe how the raw event-time variable is transformed for the
integration grid; see ?dynamicPrediction for details.
survival_variable_all <- list("Tyears1", "Tyears2", "Tyears3", "Tyears4")
survival_trans_function <- list(
fun1 = function(x) abs(x - 1),
fun2 = function(x) abs(x - 3),
fun3 = function(x) abs(x - 5),
fun4 = function(x) abs(x - 7)
)
data_raw_predict <- pbc3[pbc3$id == 2, ]
data_predict_all <- list(data_raw_predict, data_raw_predict)
risk <- dynamicPrediction(
data_predict_all, long_fit_all, survival_fit_all,
prediction_time = 5, horizon = 1, time_variable = "year",
survival_variable_all, survival_trans_function,
bandcount1 = 10, bandcount2 = 20
)
risk
#>
#> =================================================================
#> Dynamic Prediction - Event Risk
#> -----------------------------------------------------------------
#> Competing risks : Yes
#> Subjects : 1
#> -----------------------------------------------------------------
#>
#> Subject Cause 1 Risk Cause 2 Risk Total Risk
#> S1 0.0271 0.0012 0.0283
#>
#> =================================================================survival_variable_all/survival_trans_function
almost always follow the same convention shown above: variables named
"Tyears1", "Tyears2", …, each the absolute
distance from a fixed cut point. survivalTrans() builds
exactly that pair from a plain vector of cut points, so you do not have
to hand-write two matching parallel lists:
trans <- survivalTrans(c(1, 3, 5, 7))
identical(trans$survival_variable_all, survival_variable_all)
#> [1] TRUE
trans$survival_trans_function[[1]](2)
#> [1] 1risk_prob_1 and risk_prob_2 are the
predicted probabilities of experiencing each of the two competing event
types within the one-year horizon, conditional on the patient’s
longitudinal history and survival to prediction_time.
(risk_prob_2 is NULL whenever
survival_fit_all was fit without
form_conditional_cr.)
dynamicPredictionBio()dynamicPredictionBio() answers a different question: not
whether an event happens, but what a specific biomarker’s value
is likely to be at prediction_time + horizon, conditional
on survival. bio_i selects the biomarker by its position in
long_fit_all (1 = serBilir
here).
bio_pred <- dynamicPredictionBio(
bio_i = 1, data_predict_all, long_fit_all, survival_fit_all,
prediction_time = 5, horizon = 1, time_variable = "year",
survival_variable_all, survival_trans_function,
bandcount2 = 20, bandcount3 = 50
)
bio_pred$Y_predict
#> [1] 1.03628Y_predict is the MAP (most likely) predicted value;
Y_density/Y_all give the full predicted
density over a grid of candidate values, which is what
predictPlot() visualizes (see
?predictPlot).
bandcount1, bandcount2,
bandcount3dynamicPrediction() and
dynamicPredictionBio() estimate their outputs by numerical
integration over patient-specific time and biomarker grids. The
bandcount* arguments control how fine those grids are:
bandcount1
(dynamicPrediction() only) is the number of grid points
spanning the prediction window itself, from prediction_time
to prediction_time + horizon. This is the numerator of the
risk probability.bandcount2 is the number of grid
points spanning [prediction_time, upper_bound], where
upper_bound is set internally to twice the longest observed
survival/censoring time among at-risk patients. This approximates
“integrating out to infinity” for the denominator that normalizes the
risk probability. A wider follow-up range needs a larger
bandcount2 to keep the grid spacing comparable.bandcount3
(dynamicPredictionBio() only) is the number of points in
the candidate-biomarker-value grid (Y_all) used to build
the predicted density and locate its mode (Y_predict). It
controls the resolution of the density curve, not a time integral.Larger values give more accurate, smoother results at the cost of
more computation. There is no universal “correct” value because it
depends on how quickly the underlying hazard and biomarker trajectories
change and on the length of follow-up in your data — so the practical
approach is a convergence check: run the prediction
once with the defaults, once with every bandcount* doubled,
and confirm the results barely move.
risk_default <- dynamicPrediction(
data_predict_all, long_fit_all, survival_fit_all,
prediction_time = 5, horizon = 1, time_variable = "year",
survival_variable_all, survival_trans_function,
bandcount1 = 10, bandcount2 = 20
)
risk_doubled <- dynamicPrediction(
data_predict_all, long_fit_all, survival_fit_all,
prediction_time = 5, horizon = 1, time_variable = "year",
survival_variable_all, survival_trans_function,
bandcount1 = 20, bandcount2 = 40
)
abs(risk_default$risk_prob_1 - risk_doubled$risk_prob_1)
#> [1] 0.003174651
abs(risk_default$risk_prob_2 - risk_doubled$risk_prob_2)
#> [1] 0.0001377631If doubling the bandcount* values changes the result by
more than you can tolerate, keep doubling until it doesn’t; if it barely
changes anything (as above), the smaller, cheaper value is fine to use.
The same check applies to bandcount3 for
dynamicPredictionBio().
The pipeline functions validate their arguments before doing any
model fitting or numerical integration, so common mistakes fail fast
with a message that names the offending argument — instead of a cryptic
error from deep inside model.matrix() or list indexing. For
example, passing a single data.frame instead of a list to
dynamicPrediction():
dynamicPrediction(
data_predict_all[[1]], long_fit_all, survival_fit_all,
prediction_time = 5, horizon = 1, time_variable = "year",
survival_variable_all, survival_trans_function,
bandcount1 = 10, bandcount2 = 20
)
#>
#> =================================================================
#> Dynamic Prediction - Event Risk
#> -----------------------------------------------------------------
#> Competing risks : Yes
#> Subjects : 1
#> -----------------------------------------------------------------
#>
#> Subject Cause 1 Risk Cause 2 Risk Total Risk
#> S1 0.0271 0.0012 0.0283
#>
#> =================================================================or a formula referencing a column that doesn’t exist:
?predictPlot and ?riskPlot visualize the
outputs of
dynamicPredictionBio()/dynamicPrediction() for
a single patient.?cmtPlot plots observed longitudinal trajectories
stratified by eventual outcome, useful for checking whether a biomarker
looks informative before fitting.