Package {BJM}

Recode ordinal biomarker responses onto the training categories

Description

Shared helper for predictRisk, dynamicPredictionBio, and dynamicPredictionBioAll. The copula densities turn an ordinal response into its category code with as.numeric() and index the fitted thresholds with it, and the candidate categories of an ordinal bio_i are read off levels(). Both used to come from the prediction data itself, so any difference from the training categories – unused levels dropped (e.g. by droplevels() or subsetting and re-creating the factor), a different level order, or a character column – silently bracketed the latent score between the wrong thresholds (in one check, a risk of 0.53 instead of 0.10). Each ordinal response is now rebuilt as an ordered factor on the categories the clmm() fit was estimated with, matching by label; a value that is not one of those categories is an error. Missing values stay missing (drop_missing_longitudinal() removes them afterwards).

Usage

align_ordinal_levels(data_predict_all, long_fit_all)

Value

data_predict_all, with each ordinal response recoded.

Set a patient's future covariates

Description

For each row of future_rows (one per simulated time), the covariates of the last row of future_covariates (for this patient, if it has an id column) at or before that time; a covariate it does not set keeps its value in future_rows.

Usage

apply_future_covariates(
  future_rows,
  future_covariates,
  time_variable,
  id,
  patient_id
)

Call a single survival_trans_function element and validate its output

Description

assert_survival_trans() only probes each transform once, at a single representative time value, before the prediction grid runs. That catches a transform that is broken everywhere (wrong return type/length, or throws), but not one that only misbehaves away from the probe point – e.g. log(x - 10), which is fine near the probe value but returns NaN once the internal prediction grid (which can range up to 2 * max(observed survival time)) reaches x <= 10. This helper wraps every actual call site inside conditionalYT()/ conditionalYTBio()/conditionalYDT()/conditionalYDTBio() so a bad value is caught immediately, with a clear error, instead of silently corrupting a data.frame column or surfacing later as a cryptic "replacement has ... rows, data has ..." error.

Usage

apply_survival_trans(fun, x, surv_i)

Assert that every element of a list is a formula

Description

Shared input-validation helper for long_sub_fixed/ long_sub_random-style arguments, after they have been normalized to a list (a bare formula is wrapped in list() by the caller before calling this).

Usage

assert_all_formulas(x, arg_name)

Assert that a bandcount argument is a positive number or "auto"

Description

Shared input-validation helper for the bandcount1/ bandcount2/bandcount3 arguments of predictRisk(), dynamicPredictionBio(), predictPlot(), and riskPlot(), which now accept either an explicit positive number (the original behavior) or the literal string "auto" to have the value chosen automatically (see auto_tune_bandcount()).

Usage

assert_bandcount(x, arg_name)

Assert that biomarker_type is a valid per-biomarker type vector

Description

Shared input-validation helper for the biomarker_type argument of longitudinalSub: when supplied explicitly (as opposed to being auto-detected from each biomarker's response column), it must be a character vector with exactly one entry per longitudinal outcome, and every entry must be one of "continuous" or "ordinal".

Usage

assert_biomarker_type(
  biomarker_type,
  M,
  valid_types = c("continuous", "ordinal")
)

Assert that an object is the output of a specific BJM fitting function

Description

Shared input-validation helper: checks the S3 class tag attached by survivalSub()/longitudinalSub(), so passing the wrong object (or the arguments in the wrong order) fails immediately with a clear message instead of deep inside the prediction code.

Usage

assert_class(x, expected_class, arg_name, expected_source)

Assert that an object is a non-empty data.frame

Description

Shared input-validation helper: raises a clear error instead of letting a malformed argument fail deep inside model-fitting code with a cryptic message.

Usage

assert_data_frame(x, arg_name)

Assert that an object is a list of data.frame objects of the expected length

Description

Shared input-validation helper for the data_fit_all/ data_predict_all arguments, which must be a list with one data.frame per longitudinal outcome. When allow_bare_df = TRUE, a single bare data.frame is also accepted, matching the auto-repeat convenience some functions apply for a single data.frame shared across every outcome.

Usage

assert_data_list(x, arg_name, n_expected, allow_bare_df = FALSE)

Assert that an index is a valid, in-range biomarker position

Description

Shared input-validation helper for bio_i: catches an out-of-range or non-integer value before it becomes a "subscript out of bounds" error from indexing into data_predict_all/lfit.

Usage

assert_index(x, max_value, arg_name, context)

Check that the survival time enters the event-type model linearly

Description

Used by survivalSub(): at prediction time the event-type model form_conditional_cr is evaluated at every point of the survival-time integration grid by multiplying the survival time's coefficient by the grid value. That is only correct when the survival time appears as a plain main effect (status ~ time + ...); a transformation (log(time), I(time^2)) or an interaction (time:age) would silently be evaluated at the wrong value, so it is rejected here. Leaving the survival time out entirely is allowed (event type then does not depend on event time).

Usage

assert_linear_time_term(form_conditional_cr, survival_variable)

Arguments

Assert that an optional (Suggests-only) package is installed

Description

Shared input-validation helper for functions that depend on a package listed in Suggests rather than Imports (so that most users – and CRAN's own checks – are unaffected by a large, optional dependency they never call). Every call site into that package must still be fully namespaced (pkg::fun()), never @importFrom; this helper just turns a missing package into a clear, actionable error instead of a cryptic "could not find function" failure deep inside the calling code.

Usage

assert_package_installed(pkg, context)

Assert that a list of longitudinalSub.BJM fits is poolable with Rubin's rules

Description

Shared input-validation helper for poolLongitudinalSub: checks that long_fit_all_list is a list of at least two longitudinalSub.BJM objects (the S3 class longitudinalSub attaches to its return value), and that every one of them was fit for the same set of biomarkers with the same fixed-effect coefficients – otherwise "the mean of the estimates" would be averaging unrelated quantities across completions.

Usage

assert_poolable_longitudinal_fits(
  long_fit_all_list,
  arg_name = "long_fit_all_list"
)

Assert that an object is a single positive integer

Description

Shared input-validation helper for n_cores: catches a non-integer, zero, negative, or non-scalar value before it reaches parallel::mclapply()'s own (less informative) mc.cores validation.

Usage

assert_positive_integer(x, arg_name)

Assert that an object is a single, non-missing numeric value

Description

Assert that an object is a single, non-missing numeric value

Usage

assert_scalar_numeric(x, arg_name, positive = FALSE)

Assert that an object is a single, non-missing character string

Description

Assert that an object is a single, non-missing character string

Usage

assert_string(x, arg_name)

Assert that survival_variable_all/survival_trans_function are consistent

Description

Shared input-validation helper for predictRisk(), dynamicPredictionBio(), predictPlot(), and riskPlot(): the two arguments must have matching length, and every transform must be a function. When probe_value is supplied, every transform is also test-called once on it, and must return a single, finite, non-missing numeric value. Without this, a transform that throws an error, or returns a character value, a length != 1 vector, or a non-finite value (e.g. log(x) evaluated at x <= 0), would only surface deep inside the per-patient prediction grid built by conditionalYT()/ conditionalYDT()/conditionalYTBio()/conditionalYDTBio() – as a cryptic error, or, worse, as silently corrupted data with no error at all. The probe is a single call per transform, so it is cheap even though the same transform is later called many times inside the prediction grid.

Usage

assert_survival_trans(
  survival_variable_all,
  survival_trans_function,
  probe_value = NULL
)

Check that a longitudinal sub-model uses the event time

Description

Check that a longitudinal sub-model uses the event time

Usage

assert_uses_event_time(formula_i, i, survival_variable, survival_variable_all)

Arguments

Value

Invisibly TRUE; an error naming biomarker i otherwise.

Assert that a formula's variables are all present in a data.frame

Description

Shared input-validation helper: catches missing columns before they surface as an opaque error from deep inside coxph(), lme(), or model.matrix().

Usage

assert_vars_in_data(vars, data, source_name, data_name)

Auto-select "auto" bandcount arguments by doubling until convergence

Description

Shared implementation backing bandcount1/ bandcount2/bandcount3 = "auto" support in predictRisk()/dynamicPredictionBio(), and the bandcount pre-resolution done once, up front, by predictPlot()/ riskPlot() (so their internal horizon/landmark loops do not repeat the auto-tuning search on every iteration).

Starts every argument named in auto_names at its entry in bandcount_auto_start(), doubles all of them together each round, and compares consecutive results with max_relative_diff() until the largest relative change drops below tol, or max_rounds extra doublings have been tried – a hard cap, so this never loops indefinitely: at most max_rounds + 1 calls to predict_fun (the default max_rounds = 2 means at most 3 calls). If the cap is hit without converging, a warning is issued and the result/bandcount at the largest value tried is returned anyway, rather than erroring, so automated pipelines are not interrupted.

Usage

auto_tune_bandcount(
  predict_fun,
  args,
  auto_names,
  tol = 0.01,
  max_rounds = 2,
  multiplier = 2
)

Arguments

Value

A list with result (predict_fun's return value at the resolved bandcount) and bandcount (a named list of the resolved numeric bandcount value(s), one per element of auto_names).

Auto-select "auto" bandcount3 for a single biomarker's per-marker step

Description

dynamicPredictionBioAll()'s counterpart to auto_tune_bandcount(): bandcount3 only controls compute_bio_marker_step()'s own candidate-value grid (Y_all), not the shared step, so it is tuned per biomarker by re-calling compute_bio_marker_step() directly against an already-computed shared object (from compute_bio_shared_step()), rather than re-running the whole dynamicPredictionBio() pipeline – which would recompute the shared denominator on every doubling round, for every biomarker, exactly the redundant work dynamicPredictionBioAll() is meant to avoid. Same doubling-until-stable check as auto_tune_bandcount() (compares Y_predict via max_relative_diff(), capped at max_rounds doublings, warns instead of erroring if not converged by then).

Usage

auto_tune_marker_bandcount3(
  shared,
  bio_i,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  horizon,
  time_variable,
  survival_variable_all,
  survival_trans_function,
  tol = 0.01,
  max_rounds = 2,
  multiplier = 2
)

Arguments

Value

A list with result (compute_bio_marker_step()'s return value at the resolved bandcount3) and bandcount3 (the resolved numeric value).

Starting values for "auto" bandcount doubling

Description

The built-in starting point auto_tune_bandcount() doubles from for each bandcount argument. bandcount1/ bandcount2/bandcount3 used to default to fixed numbers (10, 40, and 300 respectively, across predictRisk()/dynamicPredictionBio()); those same numbers are reused here as the starting point for auto-tuning, so that the first call auto_tune_bandcount() makes matches what a caller relying on the old fixed defaults would have gotten. This cannot instead be read off formals(predict_fun), because that default is now the literal string "auto" itself. The bandcount3 start was lowered from 300 to 100 when its grid started covering only where the predictive densities are not negligible (see continuous_value_grid()): 100 points there are finer than 300 were over the old fixed, wide range.

Usage

bandcount_auto_start

Format

An object of class numeric of length 3.

Barnard & Rubin (1999) small-sample pooled degrees of freedom

Description

Internal helper for rubin_pool_scalar. Interpolates between the classical Rubin (1987) degrees of freedom (df_old, which ignores the complete-data sample size and can therefore exceed it) and the complete-data degrees of freedom (dfcom), weighted by how much of the total variance is due to missingness (lambda).

Usage

barnard_rubin_df(m, lambda, dfcom)

Ordinal box probabilities by adaptive Gauss-Hermite quadrature

Description

For each column k, \log P(lower < m + l u + e \le upper) with u \sim N(0, 1) and e \sim N(0, I) independent: \log \int \prod_j P_j(u) \phi(u) du, where P_j is a univariate normal probability. The integrand is log-concave in u; Newton's method finds its mode and the curvature there, which place an n_nodes-point Gauss-Hermite rule, all columns at once. With 20 nodes the relative error was below 10^{-5} in every check.

Usage

box_prob_quadrature(m, lower, upper, l, n_nodes = 20)

Arguments

Value

A vector of log probabilities.

Reshape recurrent-event data into gap-time episodes for the backward joint model

Description

Experimental. Prepares recurrent-event data for the backward joint model by a gap-time (renewal) reformulation: it splits each subject's event history into episodes – the interval from one restart (the previous event, or the discharge from the stay it began) to the next event, with the clock reset to 0 at the start of each – and re-clocks any accompanying longitudinal measurements onto the same per-episode timeline. The result plugs into the existing survivalSub()/ longitudinalSub()/predictRisk() pipeline, unchanged, to predict for a subject currently at risk in an open episode whether their next event (within a horizon) is a recurrence or a terminal event; predictRecurrentCount builds on the same data to predict the number of recurrences in a window. Unlike the single-event and competing-risks methods the rest of the package rests on (Shen and Li, 2021; Li et al., 2026 – see citation("BJM")), this extension is unpublished, and has been checked only in simulation studies, not on real data (see the "Recurrent events" article on the package website).

Usage

buildRecurrentEpisodes(
  data_events,
  data_long = NULL,
  id_variable = "id",
  time_variable = "time",
  status_variable = "status",
  time_variable_long = NULL,
  episode_variable = "episode",
  gap_variable = "gap",
  time_in_episode_variable = "time_in_episode",
  event_any_variable = "event_any",
  event_type_variable = "event_type",
  restart_variable = NULL,
  lookback = 0,
  episode_id_variable = "episode_id"
)

Arguments

Details

What counts as an episode. Episode 1 runs from time 0 to a subject's first row in data_events; episode k > 1 runs from the end of episode k - 1 to that subject's k-th row. A subject contributes no further episodes after a row with status %in% c(0, 2) (administrative censoring or a terminal event); any later rows for that id are ignored, with a message saying how many. status == 1 rows (recurrences) always start a new episode.

Feeding the result into survivalSub(). Fit the gap-time recurrence hazard and the recurrence-vs-terminal-event type model exactly as survivalSub() already fits a marginal hazard and a competing-risks type model, using episode_data in place of one-row-per-subject data, Surv(gap, event_any) ~ covariates in place of form_marginal_surv, and event_type ~ gap + covariates in place of form_conditional_cr (event_type is 0 for a recurrence, 1 for a terminal event – so predictRisk()'s risk_prob_1/risk_prob_2 come out as the recurrence/terminal risk respectively). Administratively-censored episodes get event_any = 0 and event_type = NA, exactly like status3 = 0/status4 = NA for a censored subject in pbc3. The episodes of one subject are then treated as independent: the coefficients are still consistent for the gap-time model, but the standard errors that summary() reports ignore the correlation between a subject's episodes and are too small.

Feeding the result into longitudinalSub()/predictRisk(). gap plays the same dual role years plays for pbc3: it is both the Surv() time in episode_data and the backward covariate swept over candidate values inside predictRisk(), so it must appear in long_sub_fixed (e.g. biomarker ~ time_in_episode + gap + gap * time_in_episode + covariates), and the training data passed to longitudinalSub() must be filtered to event_any == 1 rows of long_data – only completed episodes have an observed gap – exactly as the package's own examples filter pbc3 to status3 == 1 before fitting longitudinalSub(). time_in_episode (elapsed time since the current episode began) plays the role year plays for pbc3. To predict from "now" for a subject currently at risk in an open episode, add a trailing status = 0 row to data_events at the current data-cutoff time, and call predictRisk() with prediction_time set to that subject's elapsed time_in_episode at the cutoff.

Random effects: group by episode, not by patient. Fit the biomarker model with the random effects grouped by episode_id (random = ~ time_in_episode | episode_id), so that each episode is its own "subject", as it already is in the survival sub-model. A patient's biomarker level generally differs from one episode to the next (the underlying trajectory keeps moving), so a single random effect per patient is the wrong unit; in a simulation in which the biomarker's current level drives the risk, grouping by patient gave a clearly worse AUC than grouping by episode. predictRisk() and predictRecurrentCount() then take one episode_id per prediction, so give each prediction its own id.

The open episode's gap is NA. For an episode that has not ended (event_any == 0) – typically the one a subject is currently at risk in – the gap column of long_data is NA: its final length is exactly what is being predicted. predictRisk() reads a missing survival-time variable as "event time not yet known, still at risk", so that episode's rows can be passed straight to it, with no placeholder to set or override (episode_data keeps the real elapsed time of such an episode, which is what the Cox fit needs as a censoring time). Training rows are unaffected, since longitudinalSub() is fit on event_any == 1 rows only. Note that predictRisk() also expects the outcome columns it was fit with (the time and event-type variables) to exist in the prediction data; long_data already carries event_type (NA for an open episode).

Episodes that start at discharge (restart_variable, lookback). By default an episode starts at the previous event (time 0 for the first). For hospitalization data that is not the at-risk clock: after a readmission the patient is in hospital, not at risk, until discharge. Give data_events a numeric column – passed as restart_variable – holding, for each recurrence row, the time the process restarts (the discharge ending that stay); the next episode then starts there and its gap is measured from it. Time 0 is the first episode's start (e.g. the index discharge). Biomarkers measured during the stay, or before the index discharge, precede the episode they predict: lookback (default 0, none) keeps measurements up to that long before an episode's start – never reaching back past the previous event, so a stay's own labs are not mixed into the earlier at-risk period – and gives them negative time_in_episode. Use Inf to keep the whole stay (and, for the first episode, everything before time 0). Negative times are fine for longitudinalSub() and predictRisk().

Death, or the end of follow-up, during a hospital stay. With restart_variable, a terminal event or censoring row whose time is at or before the restart time of the readmission before it happened while the patient was in hospital, not at risk. That readmission then ends the subject's process: no further episode is opened, the row is dropped with a message, and so a death in hospital is not counted as an out-of-hospital terminal event (predictRecurrentCount()'s p_death is the probability of dying after discharge). Record such a stay's restart time as the time of death or of the end of follow-up, or later, when the discharge never happened. A data cutoff during a stay gives no open episode to predict for.

Covariates that change over time. Each episode's covariates come from the row of data_events that ends it, as in counting-process data, where each interval carries the covariate values in force during it. A covariate that changes over time must therefore be recorded on each row as its value at the start of that row's episode (age at the discharge that started it, the number of earlier admissions, the length of the previous stay), and the open episode's row in the same way. A value measured at the event itself – age at readmission, say – depends on how long the episode lasted, which is what the survival sub-model explains, and so biases it; a value measured at the time of the prediction is not what the model was fit with. Biomarkers that change within an episode belong in data_long. Covariates built from earlier episodes can also be added to episode_data afterwards, e.g. the length of the previous gap:

ed <- episodes$episode_data
ed$previous_gap <- ave(ed$gap, ed$id, FUN = function(g) c(0, g[-length(g)]))

Value

A list with elements:

episode_data

One row per subject per episode: id_variable, episode_variable (episode number, starting at 1), gap_variable (that episode's gap time), event_any_variable (0/1: did the episode end in an event of either type, i.e. is gap an observed event time rather than an administrative censoring time), event_type_variable (0 = recurrence, 1 = terminal event, NA when event_any == 0), and every other column of data_events.

long_data

NULL if data_long was not supplied; otherwise data_long, in the same bare-data.frame-or-list shape it was supplied in, with episode_variable, gap_variable, time_in_episode_variable, event_any_variable, and event_type_variable columns added (and out-of-range rows dropped, as described above). Any other column of episode_data (baseline or event-time-varying covariates, e.g. age/sex) is also carried onto matching rows here, the way every row of pbc3's long-format table already repeats age/sex/years/status3 – but only into a column data_long does not already have of that name.

Examples

## A minimal example: reshaping alone (no model fitting).
## Subject 1: two recurrences, then a terminal event; subject 2: censored.
data_events <- data.frame(
  id     = c(1, 1, 1, 2),
  time   = c(1.2, 2.5, 4.0, 3.1),
  status = c(1,   1,   2,   0),
  age    = c(58,  58,  58,  61)
)
episodes <- buildRecurrentEpisodes(data_events)
episodes$episode_data

## Hospitalizations: each episode starts at DISCHARGE, not at the readmission.
## `time` = readmission (or death/censoring) day; `discharge` = day the stay
## it started ends. Labs drawn in the stay (and before the index discharge at
## day 0) are kept with negative time_in_episode via lookback = Inf.
events <- data.frame(id = 1, time = c(20, 60, 90), status = c(1, 1, 0),
                     discharge = c(26, 66, NA))
labs <- data.frame(id = 1, day = c(-3, 0, 5, 21, 26, 40, 61, 66, 80))
hf <- buildRecurrentEpisodes(events, labs, time_variable_long = "day",
                             restart_variable = "discharge", lookback = Inf)
hf$episode_data   # gaps measured from each discharge: 20, 34, 24
hf$long_data[, c("day", "episode", "time_in_episode", "gap")]


## A full pipeline on simulated recurrent-event + longitudinal data:
## build episodes -> survivalSub() -> longitudinalSub() -> predictRisk() and
## predictRecurrentCount().
set.seed(1)
n <- 60
data_events <- do.call(rbind, lapply(seq_len(n), function(id) {
  k <- sample(1:3, 1)
  gaps <- rexp(k, rate = 1)
  data.frame(
    id = id, time = cumsum(gaps),
    status = c(rep(1, k - 1), sample(c(0, 2), 1, prob = c(0.6, 0.4))),
    age = round(50 + stats::rnorm(1, sd = 8)), sex = id %% 2
  )
}))

## subject 9999 is still at risk right now, 0.7 time units into their
## 2nd episode (1 prior recurrence, no resolution yet) -- this is who
## we will predict for below.
data_events <- rbind(data_events, data.frame(
  id = 9999, time = c(0.8, 1.6), status = c(1, 0), age = 55, sex = 1
))

data_long <- do.call(rbind, lapply(c(seq_len(n), 9999), function(id) {
  maxt <- max(data_events$time[data_events$id == id]) + 0.2
  tt <- seq(0, maxt, by = 0.25)
  data.frame(id = id, year = tt, biomarker = rnorm(length(tt)))
}))

built <- buildRecurrentEpisodes(data_events, data_long, time_variable_long = "year")

## gap-time recurrence hazard + recurrence-vs-terminal-event type model
survival_fit_recur <- survivalSub(
  built$episode_data,
  Surv(gap, event_any) ~ age + sex + episode,
  event_type ~ gap + age + sex + episode
)

## longitudinal sub-model, trained on completed episodes only, with the
## random effects grouped by episode (see Details)
long_fit_recur <- longitudinalSub(
  built$long_data[built$long_data$event_any == 1, ],
  biomarker ~ time_in_episode + age + sex + gap + gap * time_in_episode,
  ~ time_in_episode | episode_id
)

## subject 9999's open episode: its gap is NA, which predictRisk() reads as
## "still at risk" -- nothing to override
open <- built$long_data[built$long_data$id == 9999 & built$long_data$event_any == 0, ]
landmark <- max(open$time_in_episode)

## the risk that their *next* event, within 1 time unit, is a recurrence /
## a terminal event
risk <- predictRisk(
  open, long_fit_recur, survival_fit_recur,
  prediction_time = landmark, horizon = 1, time_variable = "time_in_episode",
  survival_variable_all = list(), survival_trans_function = list(),
  bandcount1 = 10, bandcount2 = 10
)
risk$risk_prob_1 # P(next event within horizon is a recurrence)
risk$risk_prob_2 # P(next event within horizon is the terminal event)

## the number of recurrences over the next 3 time units
predictRecurrentCount(
  open, long_fit_recur, survival_fit_recur,
  prediction_time = landmark, horizon = 3, time_variable = "time_in_episode",
  n_sim = 100, seed = 1
)

Build the per-patient longitudinal design matrices

Description

Shared helper for conditionalYT and conditionalYDT: constructs the longitudinal outcome matrix, fixed-effect parameter matrix, and random-effect variance-covariance pieces used by the conditional density quadratic form, for a single patient.

Usage

build_conditional_design(
  rep_num_i_list,
  data_num_i_list,
  lfit,
  Sigma,
  sigma.longitudinal,
  time_variable,
  n_longitudinal,
  long_sub_random
)

Value

A list with longitudinal_all_matrix, parameter_matrix, and cov_fac, the patient's factored covariance from cov_factor().

Build the per-patient longitudinal design pieces for a mixed continuous/ordinal (Gaussian-copula) joint model

Description

Copula-aware counterpart to build_conditional_design() (see R/conditionalDesign.R), for use whenever long_fit_all$biomarker_type contains at least one "ordinal" biomarker (see longitudinalSubCopula()). Continuous biomarkers' rows are handled identically to build_conditional_design(): their exact observed value is recorded. For an ordinal biomarker's rows, the exact underlying latent Gaussian score is unobserved – only the pair of cumulative-link thresholds (long_fit_all$thresholds[[m]]) bracketing the observed category is known (the classical Albert & Chib probit data-augmentation representation, same convention as longitudinalSubVarCopula() in R/copulaLongitudinal.R: latent score Z = eta + epsilon, epsilon ~ N(0, 1), category k iff full_alpha[k] < Z <= full_alpha[k+1] with full_alpha = c(-Inf, alpha_m, Inf), thresholds used directly with no sign flip).

The random-effects design (A_i) and Sigma_all construction is unchanged in form from build_conditional_design() – marker type only affects sigma.longitudinal (already fixed at 1 for ordinal markers by the caller, the probit identification constraint; see longitudinalSubCopula()), which flows in exactly the same way for both continuous and ordinal markers.

Usage

build_conditional_design_copula(
  rep_num_i_list,
  data_num_i_list,
  lfit,
  Sigma,
  sigma.longitudinal,
  time_variable,
  n_longitudinal,
  biomarker_type,
  long_sub_fixed,
  thresholds,
  long_sub_random
)

Value

A list with y_all_vec (stacked observed values; NA at ordinal-marker rows), alpha_lower_vec/alpha_upper_vec (stacked threshold bounds; NA at continuous-marker rows), row_marker_type (stacked per-row "continuous"/"ordinal" label), parameter_matrix (block-diagonal fixed-effect coefficients, lfit[[i]]$coefficients$fixed for continuous markers, lfit[[i]]$beta – no intercept – for ordinal markers), and Sigma_all.

Build the conditional diffusion denoising network

Description

Internal helper for imputeLongitudinal's method = "diffusion" backend. A single MLP predicts the noise in the target entries of a row from: the noisy target entries (zero elsewhere), the conditioning entries' values (zero elsewhere), the conditioning mask, the covariates and an embedding of the diffusion step – a tabular version of CSDI's conditional denoiser (Tashiro et al., 2021).

Usage

build_diffusion_module(n_cov, M, hidden_units, time_emb_dim)

Arguments

Value

An instantiated torch::nn_module.

Build the longitudinal outcome matrix for all candidate biomarker values

Description

Shared helper for conditionalYTBio and conditionalYDTBio: for each candidate value in Y_all, assembles the row of observed longitudinal outcomes across biomarkers, substituting the candidate value for the biomarker being predicted.

Usage

build_longitudinal_matrix_bio(
  data_num_i_list,
  lfit,
  bio_i,
  Y_select_all,
  n_longitudinal,
  Y_all
)

Value

A matrix with length(Y_all) rows, one per candidate value.

Build the MIWAE encoder/decoder network

Description

Internal helper for imputeLongitudinal. The encoder input is (standardized covariates, zero-filled standardized biomarkers, missingness mask); the decoder input is (latent draw, standardized covariates), matching the "condition the generative model on fully-observed covariates too" design used throughout.

Usage

build_miwae_module(n_cov, M, hidden_units, latent_dim)

Plot calibration of dynamic risk predictions

Description

Checks whether the predicted risks are numerically right, not just well ranked. At each landmark s in prediction_time, every subject in data_predict_all still event-free at s gets a predicted risk from predictRisk of an event in (s, s + horizon], using only their longitudinal history up to s. The subjects are split into n_groups groups of about equal size by predicted risk, and for each group the mean predicted risk (horizontal axis) is plotted against the observed risk (vertical axis): one minus the Kaplan–Meier estimate at s + horizon or, with competing risks, the Aalen–Johansen cumulative incidence of that event type, both of which account for censoring, with 95\ intervals. Points on the diagonal mean the predictions are well calibrated; points below it mean the model overestimates the risk, above it that it underestimates it.

Evaluating on the data the sub-models were fit on gives an optimistic (apparent) calibration; pass held-out validation data in data_predict_all when available.

For an interval-censored fit (experimental), the observed risk of a group is a nonparametric (Turnbull-type) estimate from the subjects' intervals, left-truncated at each subject's last visit up to the landmark – with competing risks, the cumulative incidence estimate of Hudgens, Satten and Longini (2001) – so it estimates the same quantity predictRisk() predicts, without using the model. Within an interval the timing of the event is not identified, so the estimate is less precise than with exact times, and no confidence interval is drawn. In simulations it was unbiased when the visits were at most about horizon apart, but with visits further apart than horizon the window mostly falls inside the subjects' intervals and the observed risk of the higher-risk groups was much too high (up to +0.27 with visits every 2–6 years and a 2-year window). An observed risk computed from the model's own posterior was more precise, but followed the predictions of a deliberately miscalibrated model instead of the truth, so it is not offered.

Usage

calibrationPlot(
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  horizon,
  time_variable,
  survival_variable_all,
  survival_trans_function,
  n_groups = 10,
  bandcount1 = "auto",
  bandcount2 = "auto"
)

Arguments

Value

A ggplot object with one panel per landmark (and, with competing risks, per event type). Its data element is a data.frame with one row per risk group, and columns landmark, cause, group, predicted, observed, lower, upper, n (subjects in the group) and n_cases (events of that type observed in the window).

Examples

data(pbc3)
survival_fit_all <- survivalSub(pbc3[!duplicated(pbc3$id), ],
                                Surv(years, status3) ~ age + sex, NULL)
long_sub_fixed <- list("long1" = serBilir ~ year + age + sex + years)
long_sub_random <- list("long1" = ~ year | id)
long_fit_all <- longitudinalSub(list(pbc3[pbc3$status3 == 1, ]),
                                long_sub_fixed, long_sub_random)
trans <- survivalTrans(c(1, 3, 5, 7))

# apparent calibration on the fitting data: 3-year landmark, 2-year window
calibrationPlot(pbc3, long_fit_all, survival_fit_all,
  prediction_time = 3, horizon = 2, time_variable = "year",
  trans$survival_variable_all, trans$survival_trans_function, n_groups = 5)

Predicted and observed risk per risk group at one landmark

Description

Predicted and observed risk per risk group at one landmark

Usage

calibration_groups(risk, time, status, cause, k, end, n_groups)

Arguments

Value

A data.frame with one row per group: group, predicted (mean predicted risk), observed, lower, upper (Kaplan–Meier/Aalen–Johansen estimate at end with 95\

Observed risk per risk group for interval-censored outcomes

Description

The interval-censored counterpart of calibration_groups(): within each group of predicted risk, the observed risk of an event of type k in (s, e] for a subject known to be event-free at their last visit V is (F_k(e) - F_k(s)) / (1 - F(V)) with F from npmle_interval(), averaged over the group's subjects – the same quantity predictRisk() predicts. No confidence interval is given.

Usage

calibration_groups_interval(risk, L, R, cause, V, k, s, end, n_groups)

Arguments

Value

As calibration_groups().

Check whether bandcount1/bandcount2/bandcount3 are large enough

Description

predictRisk()/predictLongitudinal() approximate the integrals behind the predicted risk probabilities (and, for predictLongitudinal(), the predicted biomarker density) with a finite grid, controlled by bandcount1/bandcount2/ bandcount3. There is no universal correct value: too few grid points silently bias the answer, and too many just cost more time, and the right value depends on the data (e.g. how wide the follow-up range is). Rather than guess, or auto-loop until some tolerance is met – which multiplies runtime unpredictably, especially for predictLongitudinal()'s nested bandcount2 x bandcount3 grid – this runs the prediction once at the bandcount value(s) you supply, once more with those value(s) scaled up, and reports the largest relative change between the two, so you can see whether you have already converged or need to increase the checked bandcount(s) and re-run. It costs exactly 2 calls to predict_fun, regardless of how many times you invoke it.

Usage

checkBandcountConvergence(
  predict_fun,
  ...,
  bandcount_args,
  multiplier = 2,
  tol = 0.01
)

Arguments

Value

An object of class "checkBandcountConvergence.BJM", a named list with elements:

base_bandcount

The bandcount value(s) checked, as supplied in bandcount_args.

scaled_bandcount

base_bandcount scaled by multiplier.

tol

The relative-change tolerance used to decide converged.

by_field

A named numeric vector giving the largest relative change, per comparable output field (e.g. risk_prob_1, Y_predict), between the base and scaled call.

max_rel_diff

The largest value in by_field, or NA if there was no comparable output (e.g. horizon <= 0).

converged

TRUE if max_rel_diff < tol, FALSE if not, NA if there was nothing to compare.

base_result

The full result of calling predict_fun at base_bandcount.

scaled_result

The full result of calling predict_fun at scaled_bandcount.

Printing the object summarizes these fields.

Examples

data(pbc3)

data_survival_fitting = pbc3[!duplicated(pbc3$id), ]
survival_fit_all = survivalSub(data_survival_fitting, Surv(years, status3) ~ age + sex, NULL)

long_sub_fixed = serBilir ~ year + age + sex + (years) + (years) * year
long_sub_random = ~ year | id
long_fit_all = longitudinalSub(pbc3[pbc3$status3 == 1, ], long_sub_fixed, long_sub_random)

trans = survivalTrans(c(1, 3, 5, 7))

data_predict_all = pbc3[pbc3$id == 2 & pbc3$year <= 3, ]

# Check bandcount1/bandcount2 together: are they both already large enough?
check = checkBandcountConvergence(
  predictRisk, data_predict_all, long_fit_all, survival_fit_all,
  prediction_time = 3, horizon = 3, time_variable = "year",
  trans$survival_variable_all, trans$survival_trans_function,
  bandcount_args = list(bandcount1 = 10, bandcount2 = 10)
)
print(check)

Validate future_covariates

Description

Checks that it is a data frame with a non-missing time column, that its id column (if any) names patients being simulated, and that it does not set an outcome. Returns it ordered by time.

Usage

check_future_covariates(
  future_covariates,
  time_variable,
  id,
  patient_ids,
  outcome_variables
)

Check the event-time transformations against the training data

Description

The longitudinal sub-models were fit on data whose survival_variable_all columns the user computed from the event time; at prediction time the same columns are recomputed with survival_trans_function. The two are passed as separate lists, so a function out of order (or a different one) used to give wrong predictions silently. This applies each function to the event times stored with every fit (lme()'s data, clmm()'s model frame) and stops if the result differs from the column the model was fit on.

Usage

check_survival_trans_training(
  long_fit_all,
  survival_fit_all,
  survival_variable_all,
  survival_trans_function
)

Arguments

Value

TRUE, invisibly; an error on a mismatch.

Plot cumulative incidence functions

Description

Plots the nonparametric cumulative incidence of each event type over time, estimated by the Aalen–Johansen estimator (via survfit on a multi-state Surv outcome). With a single event type this equals one minus the Kaplan–Meier survival curve; with competing risks each event type gets its own curve, correctly accounting for the other event types instead of treating them as censoring. An optional group_variable draws separate curves per group (e.g. treatment arm).

Usage

cifPlot(
  data_plot_all,
  survival_variable,
  event_type_variable,
  id_variable = "id",
  censor_value = 0,
  event_labels = NULL,
  group_variable = NULL,
  conf_int = TRUE
)

Arguments

Value

A ggplot object.

Examples

data(pbc3)

# competing risks: death (1) and transplantation (2); 0 = censored
cifPlot(pbc3, survival_variable = "years", event_type_variable = "status5",
  event_labels = c("1" = "Death", "2" = "Transplant"))

# the same events coded as in status4 (censored subjects are NA)
cifPlot(pbc3, survival_variable = "years", event_type_variable = "status4",
  censor_value = NA, event_labels = c("0" = "Death", "1" = "Transplant"))

# single event type, by treatment arm
cifPlot(pbc3, survival_variable = "years", event_type_variable = "status2",
  group_variable = "drug", conf_int = FALSE)

Legend labels for the event types in cifPlot()

Description

Legend labels for the event types in cifPlot()

Usage

cif_event_labels(event_values, event_labels)

Arguments

Value

A character vector of labels, one per event type.

Expand cumulative incidence confidence limits into a step-shaped ribbon

Description

geom_ribbon() joins points with straight lines; to follow the step curves, each jump time gets two rows, one with the previous limits and one with the new limits.

Usage

cif_step_ribbon(curves)

Arguments

Value

A data.frame with columns time, lower, upper, hue, curve.

Normalize a risk-probability ratio into a valid probability

Description

Shared helper for predictRisk and dynamicPredictionBio: divides summed predicted-event mass by summed total mass and clamps the result to [0, 1].

Usage

clamp_risk_prob(numerator_sum, denominator_sum)

Value

A numeric vector of risk probabilities in [0, 1].

Plot conditional mean trajectories (CMT)

Description

This function generates the Conditional Mean Trajectories (CMT) plot: the mean biomarker trajectory of the patients whose event happened in the window (condi_time2event - event_window, condi_time2event]. The time from the first measurement up to condi_time2event is cut into consecutive intervals of length interval_time; each point is the mean of those patients' biomarker values measured within one interval, plotted at the interval's end, and a straight line is fitted through the points. With competing risks there is one trajectory per event type.

All times – time_variable, survival_variable, condi_time2event, event_window and interval_time – are in the units of the data (years, months, days, ...), and the time axis shows time_variable on that scale, labelled with its name.

Only patients with an observed event enter: with competing risks, those whose event_type_variable equals one of event_type; without, those whose status_variable is not 0. Without competing risks and with status_variable = NULL, patients are selected by their survival time alone, so censored ones are included too, with a warning. An interval without any measurement is left out rather than filled in.

Usage

cmtPlot(
  data_plot_all,
  condi_time2event,
  event_type_variable,
  event_type,
  bio_variable,
  time_variable,
  survival_variable,
  interval_time = 1/12,
  id_variable = "id",
  status_variable = NULL,
  event_window = 1
)

Arguments

Value

Conditional mean trajectories plot.

Examples

# example without competing risks

data(pbc3)

pbc.cmt <- cmtPlot(data_plot_all = pbc3, condi_time2event = 5,
   event_type_variable = NULL, event_type = NULL,
   bio_variable = "serBilir", time_variable = "year",
   survival_variable = "years",
   interval_time = 1/12, status_variable = "status3"
)

pbc.cmt


# example with competing risks

data(pbc3)

pbc.cmt.cr <- cmtPlot(pbc3, condi_time2event = 5,
   event_type_variable = 'status4', event_type = c("0", "1"),
   bio_variable = "albumin", time_variable = "year",
   survival_variable = "years",
   interval_time = 1/4
)

pbc.cmt.cr

A patient's history rows usable for conditioning

Description

The rows of one biomarker's data with the biomarker and every covariate of its sub-model observed (the outcome columns, which are set at each candidate event time, excepted). May have no rows.

Usage

complete_history_rows(d, i, long_fit_all, survival_variable, other_outcomes)

Bio_i-specific computation for dynamicPredictionBio()

Description

Internal helper for the part of dynamicPredictionBio()'s pipeline that does depend on which biomarker is being predicted: building the candidate-value grid Y_all, evaluating the numerator conditional density on it, and assembling/normalizing Y_density and Y_predict. Takes the output of compute_bio_shared_step() so the bio_i-independent pieces are not recomputed.

Usage

compute_bio_marker_step(
  shared,
  bio_i,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  horizon,
  time_variable,
  survival_variable_all,
  survival_trans_function,
  bandcount3
)

Arguments

Value

A list with Y_predict, Y_density, Y_all (matching dynamicPredictionBio()'s return value fields). For an ordinal bio_i, Y_all/Y_predict hold integer category codes (1:K, in threshold order) rather than a numeric grid – see Details below – and Y_all additionally carries a "category_labels" attribute with the matching level-label strings.

Bio_i-independent shared computation for dynamicPredictionBio()

Description

Internal helper factoring out the part of dynamicPredictionBio()'s pipeline that does not depend on which biomarker (bio_i) is being predicted: restricting to at-risk patients, building the survival-side integration grid over T > prediction_time + horizon out to upper_bound (the prediction is conditional on survival through the time the biomarker is predicted at), dispatching to the Gaussian-copula conditional-density variants when needed (see predictRisk()), and evaluating the denominator conditional density (conditionalYT/conditionalYDT). dynamicPredictionBioAll() computes this once and reuses it across every requested biomarker, instead of recomputing it once per biomarker (including, under the copula path, its mvtnorm::pmvnorm() calls) the way calling dynamicPredictionBio() once per biomarker would.

Usage

compute_bio_shared_step(
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  horizon,
  time_variable,
  survival_variable_all,
  survival_trans_function,
  bandcount2
)

Arguments

Value

A list with the pieces compute_bio_marker_step() needs: data_predict_all (at-risk-filtered), survival_variable, predict.time.infinity, S_T_all_infinity, has_cr, D_T_all_infinity (NULL if no competing risk), f_y_D_all_infinity, and the dispatched conditionalYTBio_fun/conditionalYDTBio_fun functions.

conditional distribution of D|T

Description

This function computes the conditional probability density function of competing risk event type D, given the survival time T.

Usage

conditionalDT(data_predict_all, long_fit_all, survival_fit_all, l_i)

Arguments

Value

Probability matrices of competing risk event type D conditional on survival outcome T.

conditional distribution of Y|D, T, if with competing risk

Description

This function computes the conditional probability density function of longitudinal variable Y, given the survival time T with competing risk D.

Usage

conditionalYDT(
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  l_i,
  survival_variable,
  time_variable,
  survival_variable_all,
  survival_trans_function
)

Arguments

Value

The output is a list containing probability matrices of log densities (so that they neither overflow nor underflow; see predictRisk() for how they are exponentiated). In the presence of competing risks, this list includes two elements; otherwise, it contains only one element. Each element within the list is a probability matrix, with the number of rows (l_i) corresponding to specific time points and columns representing different patients. Every matrix element represents the conditional probability derived from the conditional distribution of longitudinal variable Y given the survival time T with competing risk D for a particular patient at a specific time point.

conditional distribution of Y|D, T, if with competing risk

Description

This function computes the conditional probability density function of longitudinal variable Y, given the survival time T with competing risk D.

Usage

conditionalYDTBio(
  Y_all,
  time_new,
  bio_i,
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  l_i,
  survival_variable,
  time_variable,
  survival_variable_all,
  survival_trans_function
)

Arguments

Value

The output is a list containing probability matrices of log densities (so that they neither overflow nor underflow; see predictRisk() for how they are exponentiated). In the presence of competing risks, this list includes two elements; otherwise, it contains only one element. Each element within the list is a probability matrix, with the number of rows (l_i) corresponding to specific time points and columns representing different patients. Every matrix element represents the conditional probability derived from the conditional distribution of longitudinal variable Y given the survival time T with competing risk D for a particular patient at a specific time point.

Conditional distribution of Y|D,T for a mixed continuous/ordinal (Gaussian-copula) joint model, if with competing risk – biomarker-value prediction

Description

Copula-aware counterpart to conditionalYDTBio(), used by dynamicPredictionBio() whenever long_fit_all$biomarker_type contains at least one "ordinal" biomarker and a competing-risk survival_fit_all is supplied. Both a continuous and an ordinal bio_i target are supported: for an ordinal bio_i, Y_all is expected to already be the vector of candidate category labels (the fitted factor's levels(), in threshold order) rather than a numeric grid – see compute_bio_marker_step(), which builds that vector and translates the result back into integer category codes for its caller. See conditionalYTBioCopula() for why no branching is actually needed here – build_conditional_design_copula() already looks up each row of biomarker bio_i generically via long_fit_all$biomarker_type[bio_i], whether that row is historical or the candidate row assigned below – and for why this function evaluates mvtnorm::pmvnorm() once per patient/l_i/ Y_all candidate/event-type branch, rather than one vectorized dmvnorm() call per patient/l_i as conditionalYDTBio() does.

Usage

conditionalYDTBioCopula(
  Y_all,
  time_new,
  bio_i,
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  l_i,
  survival_variable,
  time_variable,
  survival_variable_all,
  survival_trans_function
)

Arguments

Conditional distribution of Y|D,T for a mixed continuous/ordinal (Gaussian-copula) joint model, if with competing risk

Description

Copula-aware counterpart to conditionalYDT(), used by predictRisk() whenever long_fit_all$biomarker_type contains at least one "ordinal" biomarker (see longitudinalSubCopula()). See conditionalYTCopula() for the mixed continuous/ordinal density/probability construction shared by both event-type branches (w0/w1) here.

Usage

conditionalYDTCopula(
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  l_i,
  survival_variable,
  time_variable,
  survival_variable_all,
  survival_trans_function
)

Arguments

conditional distribution of Y|T, if no competing risk;

Description

This function computes the conditional probability density function of longitudinal variable Y, given the survival time T without competing risk D.

Usage

conditionalYT(
  data_predict_all,
  long_fit_all,
  l_i,
  survival_variable,
  time_variable,
  survival_variable_all,
  survival_trans_function
)

Arguments

Value

The output is a list containing probability matrices of log densities (so that they neither overflow nor underflow; see predictRisk() for how they are exponentiated). In the presence of competing risks, this list includes two elements; otherwise, it contains only one element. Each element within the list is a probability matrix, with the number of rows (l_i) corresponding to specific time points and columns representing different patients. Every matrix element represents the conditional probability derived from the conditional distribution of longitudinal variable Y given the survival time T without competing risk D for a particular patient at a specific time point.

conditional distribution of Y|T, if no competing risk;

Description

This function computes the conditional probability density function of longitudinal variable Y, given the survival time T without competing risk D.

Usage

conditionalYTBio(
  Y_all,
  time_new,
  bio_i,
  data_predict_all,
  long_fit_all,
  l_i,
  survival_variable,
  time_variable,
  survival_variable_all,
  survival_trans_function
)

Arguments

Value

The output is a list containing probability matrices of log densities (so that they neither overflow nor underflow; see predictRisk() for how they are exponentiated). In the presence of competing risks, this list includes two elements; otherwise, it contains only one element. Each element within the list is a probability matrix, with the number of rows (l_i) corresponding to specific time points and columns representing different patients. Every matrix element represents the conditional probability derived from the conditional distribution of longitudinal variable Y given the survival time T without competing risk D for a particular patient at a specific time point.

Conditional distribution of Y|T for a mixed continuous/ordinal (Gaussian-copula) joint model, if no competing risk – biomarker-value prediction

Description

Copula-aware counterpart to conditionalYTBio(), used by dynamicPredictionBio() whenever long_fit_all$biomarker_type contains at least one "ordinal" biomarker (see longitudinalSubCopula()). Both a continuous and an ordinal bio_i target are supported: for an ordinal bio_i, Y_all is expected to already be the vector of candidate category labels (the fitted factor's levels(), in threshold order) rather than a numeric grid – see compute_bio_marker_step(), which builds that vector and translates the result back into integer category codes for its caller. No branching is actually needed here: build_conditional_design_copula() already looks up each row of biomarker bio_i (historical and the candidate row assigned below) generically via long_fit_all$biomarker_type[bio_i], bracketing an ordinal row's latent score between its category's cumulative-link thresholds instead of matching it to an exact value – exactly the same mechanism already used for every other ordinal biomarker's observed history in the joint density. Assigning a candidate label into data_it_Y's ordinal factor column (below) therefore evaluates the same mixed_density_prob_copula() box probability as any other ordinal row would, with no separate code path required.

Unlike conditionalYTBio() – which can evaluate mvtnorm::dmvnorm() once per l_i across every Y_all candidate in one vectorized call, because the multivariate-normal density does not need re-normalizing per candidate – here each candidate value in Y_all changes the conditional block's mean/covariance (see mixed_density_prob_copula()) and therefore requires its own mvtnorm::pmvnorm() evaluation. This makes the copula path in this function O(patients * length(l_i) * length(Y_all)) Monte-Carlo box probability evaluations, rather than one dmvnorm() call per patient/l_i – keep bandcount3 modest for mixed fits with a continuous bio_i (an ordinal bio_i's Y_all length is fixed at its category count, not controlled by bandcount3 at all).

Usage

conditionalYTBioCopula(
  Y_all,
  time_new,
  bio_i,
  data_predict_all,
  long_fit_all,
  l_i,
  survival_variable,
  time_variable,
  survival_variable_all,
  survival_trans_function
)

Arguments

Conditional distribution of Y|T for a mixed continuous/ordinal (Gaussian-copula) joint model, if no competing risk

Description

Copula-aware counterpart to conditionalYT(), used by predictRisk() whenever long_fit_all$biomarker_type contains at least one "ordinal" biomarker (see longitudinalSubCopula()). Continuous biomarkers contribute an exact Gaussian density factor, exactly as conditionalYT() computes for every biomarker; ordinal biomarkers instead contribute a Gaussian-copula box probability (via mixed_density_prob_copula()), since only the cumulative-link category – not the exact underlying latent score – is observed for them. All-continuous fits are unaffected: they are still routed to conditionalYT() by predictRisk(), this function is never called for them, and mixed_density_prob_copula() reduces to (a constant multiple of) conditionalYT()'s own computation when there are no ordinal markers – see that function's documentation.

Usage

conditionalYTCopula(
  data_predict_all,
  long_fit_all,
  l_i,
  survival_variable,
  time_variable,
  survival_variable_all,
  survival_trans_function
)

Arguments

Log densities of an empty longitudinal history

Description

Stand-ins for conditionalYT()/conditionalYDT() for patients without any measurement: f(Y \mid T) = 1, i.e. a log density of 0 at every time (and event type), so predictRisk() gives the survival sub-model's risk.

Usage

conditionalYT_none(data_predict_all, long_fit_all, l_i, ...)

conditionalYDT_none(data_predict_all, long_fit_all, survival_fit_all, l_i, ...)

Value

A list of one (conditionalYT_none()) or two (conditionalYDT_none()) zero matrices, one row per time in l_i and one column per patient.

Per-patient pieces of the marginal survival model

Description

Shared by integration_upper_bound() and prepare_infinity_grid(): each at-risk patient's Cox linear predictor (reference = "zero", NA if a covariate is missing), and the patients grouped by the baseline cumulative hazard they use (one group, or one per stratum of a stratified model).

Usage

conditional_survival_setup(data_predict_all, long_fit_all, survival_fit_all)

Value

A list with lp, groups (each a list with cum_basehaz and patients, indices into lp), and last_time (the last time in the data survivalSub() was fit on).

Candidate-value grid for a continuous biomarker

Description

Helper for compute_bio_marker_step(): the grid of candidate values (Y_all) a continuous biomarker's predictive density is tabulated on, shared by all patients predicted. Built in two passes:

1. A coarse grid over a deliberately wide range – 5 spans either side of the patients' observed values, the span being their range but at least the biomarker's SD in the training data – with a step of the biomarker's residual SD \sigma, at least 30 and at most 100 points. A predicted future measurement includes its measurement error, so every predictive density has an SD of at least \sigma: within half a step of its mode it is still above 88\ rel_tol, so no patient's density falls between coarse points.

2. bandcount3 + 1 equally spaced points over the part of the coarse grid where some patient's density exceeds rel_tol times that patient's maximum, widened by 2\sigma (and at least one coarse step) on either side.

The wide range used to be the final grid, so most of the bandcount3 points (about 90\ every density was practically 0. If the coarse pass finds no usable density, the wide range is returned with bandcount3 + 1 points, as before.

Usage

continuous_value_grid(
  shared,
  bio_i,
  bio_i_name,
  long_fit_all,
  bandcount3,
  density_fun,
  rel_tol = 1e-06
)

Arguments

Value

A numeric vector of candidate values.

Copula EM with adaptive Gauss-Hermite quadrature

Description

Given the ordinal random effects b_o, the continuous ones b_c have a Gaussian posterior with mean a + B b_o and covariance A_c, and the continuous data contribute a Gaussian factor in b_o; only b_o is integrated numerically. For each subject the posterior mode of b_o (Newton, warm-started from the previous iteration) and the curvature there place a product Gauss-Hermite grid, which gives E[b b^T \mid y] block by block; their average over subjects is the next D. Uses 7 nodes per dimension up to two ordinal random effects and 5 beyond. The EM is accelerated by SQUAREM (Varadhan and Roland, 2008), falling back to plain EM steps where an extrapolated D is not positive definite; convergence is a change in D, scaled by \sqrt{D_{jj} D_{kk}} so that it does not depend on the biomarkers' units, below tol.em over one accelerated cycle (three EM updates).

Usage

copula_D_aghq(D, setup, tol.em = 1e-04, max.iter = 500, verbose = FALSE)

Arguments

Value

A list with D.

Copula EM by stochastic approximation (SAEM)

Description

Fallback of longitudinalSubVarCopula() for many ordinal random effects. Each iteration draws every ordinal latent score from its truncated normal given the subject's current random effects (one dimension at a time, as they are independent given b), computes E[b b^T \mid y, z] exactly, draws new random effects from their Gaussian posterior, and updates D by stochastic approximation: step 1 for n_burnin iterations, then 1/k, averaging the next n_average. The draws use a fixed seed, with the caller's random number state restored, so a fit is reproducible.

Usage

copula_D_saem(D, setup, n_burnin = 100, n_average = 200, verbose = FALSE)

Arguments

Value

A list with D.

Per-subject pieces of the copula E-step that do not depend on D

Description

For every subject: the ordinal rows' fixed-effect linear predictor, the thresholds bracketing each observed category and the random-effects design of those rows; and the continuous rows' residuals from the fixed effects, their random-effects design and residual precisions. Also the positions ro/rc of the ordinal and continuous biomarkers' random effects in D.

Usage

copula_estep_setup(thetaLong, l, biomarker_type, thresholds)

Arguments

Value

A list with subjects (one entry per subject), ro and rc.

A patient's copula log densities at every grid point

Description

Shared by conditionalYTCopula(), conditionalYDTCopula(), conditionalYTBioCopula() and conditionalYDTBioCopula(). The patient's covariance, observed values and threshold bounds do not depend on the event time, so the design (build_conditional_design_copula()) is built once per candidate value of the predicted biomarker instead of once per grid point, and the means at all grid points come from copula_mean_grid(). All densities are evaluated at once.

Predicting a continuous biomarker with the Monte Carlo box probability (two or more ordinal random effects, see mixed_density_copula_prep()) needs it at every grid point for every candidate value, which made this the slowest prediction in BJM. The candidate value only shifts the ordinal latent scores' conditional mean, so the log box probability is a smooth function of it: it is computed at n_exact candidate values spread over the grid, and a cubic spline gives it at the others, while the continuous biomarkers' density is computed exactly at every value. In checks the difference from computing every value was no larger than the Monte Carlo error itself (around 10^{-4}), at 4-8 times the speed. With one ordinal random effect, the quadrature is fast and every value is computed exactly.

Usage

copula_logdens_grid(
  data_num_i_list,
  rep_num_i_list,
  long_fit_all,
  l_i,
  survival_variable,
  time_variable,
  survival_variable_all,
  survival_trans_function,
  sigma.longitudinal,
  event_type_variable = NULL,
  bio = NULL,
  n_exact = 15
)

Arguments

Value

A list with one element per event type (1, then 0; just one without competing risks), each a length(l_i) x (number of candidate values, 1 without bio) matrix of log densities.

A biomarker's mean at every grid point, continuous or ordinal

Description

For the copula density functions: fixed_mean_grid() for a continuous biomarker; for an ordinal one, the linear predictor X(t)\beta of its clmm() fit (no intercept; the thresholds take its place) from long_sub_fixed[[i]], with all grid points stacked into one model.matrix() call. An ordinal formula with a term whose values depend on the other rows (poly(), splines::ns()/bs(), factor(), scale(), cut(), which longitudinalSub() warns about) is evaluated one grid point at a time, as before, so that its value does not change.

Usage

copula_mean_grid(
  i,
  data_i,
  long_fit_all,
  l_i,
  survival_variable,
  survival_variable_all,
  survival_trans_function,
  fixed_values = list()
)

Value

A nrow(data_i) x length(l_i) matrix.

Factor one patient's joint longitudinal covariance

Description

The stacked (across every biomarker) observation vector of one patient has covariance Sigma_all = A %*% D %*% t(A) + diag(r_diag): A is the N x r random-effects design (random_effects_design()), D the r x r random-effects covariance (Sigma_fit) and r_diag each observation's residual variance. This factors it once so that cov_logdens() can evaluate the Gaussian log-density at any number of residual vectors.

With method = "woodbury", Sigma_all is never formed. Writing D = G %*% t(G) and B = A %*% G, the matrix determinant lemma and the Woodbury identity give

\log|\Sigma| = \sum \log r_j + \log|K|,\quad e'\Sigma^{-1}e = e'R^{-1}e - u'K^{-1}u,

with K = I + t(B) %*% R^{-1} %*% B (r x r) and u = t(B) %*% R^{-1} %*% e. Factoring costs O(N r^2 + r^3) instead of O(N^3), each density O(N r) instead of O(N^2), and memory O(N r) instead of O(N^2). D^{-1} is never needed: G is t(chol(D)), or, when D is only positive semi-definite (a random-effect variance collapsing to 0 during EM), its eigenvectors scaled by the square roots of the non-negative eigenvalues, which still reproduces D exactly.

With method = "dense", Sigma_all is formed and Cholesky factored, as mvtnorm::dmvnorm() does. "auto" picks "woodbury" only when every residual variance is positive and N >= 3 r: in timings the two break even around N = 2.5 r (the Woodbury steps carry larger constants), and dense is faster below. Both give the same density up to floating-point rounding.

Usage

cov_factor(A, D, r_diag, method = c("auto", "woodbury", "dense"))

Arguments

Value

An object for cov_logdens().

Gaussian log-density from a factored covariance

Description

Evaluates the N-variate normal log-density with covariance factored by cov_factor() at each column of a residual matrix (observation minus mean).

Usage

cov_logdens(fac, E)

Arguments

Value

A length-m vector of log-densities. If the covariance is not positive definite ("dense" only), -Inf, or Inf for a zero residual, as mvtnorm::dmvnorm() returns.

Log densities of every candidate outcome vector at every grid point

Description

For conditionalYTBio()/conditionalYDTBio(): the log density of each row of x (one per candidate value of the predicted biomarker) with mean mu[, k], for every grid point k, in one cov_logdens() call.

Usage

cov_logdens_grid(fac, x, mu)

Arguments

Value

A nrow(x) x ncol(mu) matrix.

Baseline cumulative hazard at a set of time points

Description

Helper for marginalT and the integration-grid helpers: for each element of l_i, the Breslow cumulative hazard as the right-continuous step function it is – the value at the last tabulated time <= l_i, or 0 before the first. Past the table's last time it continues from the last tabulated value with the slope of the least-squares line through the whole table (a linear extrapolation of the cumulative hazard, i.e. a constant hazard). It used to switch to that line itself, which does not pass through the last tabulated value: the cumulative hazard jumped there, and when it jumped down (a convex cumulative hazard, i.e. an increasing hazard) the interval spanning the last time got a negative probability mass. The line used to be tabulated on a fixed 0.005 grid out to twice the integration upper limit and then looked up, which assumed time was measured in years (with time in days the table ran to millions of rows) and failed with "wrong sign in 'by'" when the upper limit was below the last training time; it is now evaluated directly. The table used to be read at the nearest tabulated time, which for a point just before an event time took that event's jump early.

Usage

cumulative_baseline_at(cum_basehaz, l_i)

Arguments

Value

A numeric vector, one cumulative hazard per element of l_i.

Default covariate update for the episodes after the open one

Description

Adds k to the episode counter, but never past cap: the largest episode number in the data the survival sub-model was fit on. Beyond it the model's effect of the counter is an extrapolation, and a hazard that rises with each earlier recurrence, run out to episode 30, multiplies the hazard many hundredfold, so paths drag on through ever shorter episodes and the counts are inflated.

Usage

default_next_covariates(episode_variable, cap = Inf)

Mode of a density tabulated on a grid

Description

Helper for compute_bio_marker_step(): the grid point where density peaks, refined (when refine = TRUE) by fitting a parabola through the log-density at that point and its two neighbours. For a smooth, near-normal predictive density this recovers the mode to well within one grid step, instead of snapping to the grid: the unrefined mode moved in whole grid steps as bandcount3 changed, which for a predicted value near 0 is a large relative change, so "auto" tuning of bandcount3 reported non-convergence even on a fine grid. No refinement at the grid's ends, for an ordinal biomarker (whose grid is its categories), or if a neighbour's density is 0.

Usage

density_mode(Y_all, density, refine = TRUE)

Arguments

Value

A single number, or NA if density has no finite maximum.

Build a linear variance-preserving diffusion (DDPM) noise schedule

Description

Internal helper for imputeLongitudinal's method = "diffusion" backend. Returns the standard Ho et al. (2020) linear beta schedule and the derived alpha/ alpha_bar sequences, as plain numeric vectors (not torch tensors) so diffusion_impute can index them by a plain R integer timestep inside its reverse-sampling loop without repeated tensor<->R round trips.

Usage

diffusion_beta_schedule(diffusion_steps, beta_start = 1e-04, beta_end = 0.02)

Draw conditional reverse-diffusion completions

Description

Internal helper for imputeLongitudinal's method = "diffusion" backend. The missing entries of every row are the targets and all its observed entries the conditions: starting from pure noise on the targets, the DDPM reverse process (Ho et al., 2020) is run with the network conditioned on the observed values at every step. Observed entries are returned unchanged.

Usage

diffusion_impute(
  net,
  X_cov,
  X_resp,
  X_mask,
  resp_mean,
  resp_sd,
  schedule,
  diffusion_steps,
  time_emb_dim,
  n_imputations
)

Arguments

Value

A list of n_imputations completed response matrices on the original scale.

Drop longitudinal measurements taken after the prediction time

Description

Shared helper for predictRisk, dynamicPredictionBio, and dynamicPredictionBioAll: a dynamic prediction at prediction_time may only condition on the history observed up to prediction_time, so rows of data_predict_all whose time_variable is later than that are removed, with a warning saying how many. (predictPlot/riskPlot already truncate per landmark time before calling these, so they never trigger the warning.) Rows with a missing time_variable are kept, as before. It is an error if no row is left: there is no history to condition on. With allow_empty = TRUE (predictRisk(), which predicts a patient without measurements from the survival sub-model alone) it is an error only if data_predict_all had no row to begin with.

Usage

drop_after_prediction_time(
  data_predict_all,
  time_variable,
  prediction_time,
  allow_empty = FALSE
)

Value

data_predict_all, filtered per element.

Drop prediction rows with a missing biomarker value or covariate

Description

Shared helper for predictRisk, dynamicPredictionBio, and dynamicPredictionBioAll: a patient's conditional density only involves the measurements actually observed, so for each biomarker, rows of its data_predict_all element with a missing response, or a missing covariate/time/ID used by that biomarker's long_sub_fixed/long_sub_random formulas, are removed. Previously such a row put an NA into the stacked outcome vector (or misaligned the design matrix) and made that patient's whole prediction NA. The outcome variables (survival_variable: the survival time and, with competing risks, the event type – see outcome_variables()) and survival_variable_all are not checked: they are overwritten with each integration grid point / event type before use, and are unknown (typically NA) for a patient still at risk. The event type used to be checked, so with the event type in long_sub_fixed every at-risk patient was dropped. A patient left with no rows for some biomarker cannot be predicted; a warning names them, and it is an error if that leaves no patient at all.

With no_measurement_ok = TRUE (predictRisk()), a patient left with no rows for any biomarker is not reported: it is predicted from the survival sub-model alone (see prepare_risk_data()). It is then an error only if no patient has measurements of every biomarker, none is left without any measurement, and n_other (the patients without any row up to prediction_time, already removed) is 0.

Usage

drop_missing_longitudinal(
  data_predict_all,
  long_fit_all,
  survival_variable,
  survival_variable_all,
  no_measurement_ok = FALSE,
  n_other = 0
)

Value

data_predict_all, filtered per element.

Predict a future biomarker value from fitted sub-models, for a single biomarker

Description

Internal single-biomarker engine behind predictLongitudinal – call predictLongitudinal() directly instead (it dispatches here automatically when bio_i names exactly one biomarker, and to dynamicPredictionBioAll otherwise). Kept as a separate internal function – rather than folded into predictLongitudinal() – because predictPlot and checkBandcountConvergence also call it directly for a single biomarker at a time.

Companion to predictRisk, using the same fitted longitudinal (longitudinalSub) and survival (survivalSub) sub-models, but instead of an event-risk probability this returns a predictive density for a future value of one chosen biomarker (bio_i) at prediction_time + horizon, conditional on the subject's observed longitudinal history up to prediction_time and on survival (being event-free) through prediction_time + horizon, the time the biomarker is predicted at. The density (Y_density, evaluated over a candidate-value grid Y_all) is obtained by integrating the biomarker's predictive distribution against the survival sub-model's hazard, using the subject's empirical-Bayes random-effects update from their observed history; its mode (Y_predict) is reported as the point prediction. As in predictRisk, the survival-side integrals are evaluated on numerical grids controlled by bandcount2, and the density itself is evaluated on a grid controlled by bandcount3; see Details.

bio_i may refer to either a continuous or an ordinal biomarker (see longitudinalSubCopula); for an ordinal bio_i, Y_all and Y_predict hold integer category codes (1:K, in the fitted factor's levels()/threshold order) rather than a numeric grid, with Y_all additionally carrying a "category_labels" attribute giving the matching level-label strings; bandcount3 is ignored in that case, since the candidate grid is fixed at the biomarker's category count (see Details).

The prediction is conditional on the longitudinal history observed up to prediction_time: rows of data_predict_all whose time_variable is later than prediction_time are dropped, with a warning, before predicting.

Usage

dynamicPredictionBio(
  bio_i,
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  horizon,
  time_variable,
  survival_variable_all,
  survival_trans_function,
  bandcount2 = "auto",
  bandcount3 = "auto"
)

Arguments

Details

There is no universal correct value for bandcount2/bandcount3: as a practical check, double both and confirm the results barely change; if they do, keep doubling. By default (bandcount2 = "auto", bandcount3 = "auto"), this doubling check is done for you: starting from small built-in values, both are doubled together, and the result is compared to the previous round, until the largest relative change in Y_predict drops below 1%, or 2 doublings have been tried (so at most 3 calls' worth of work). If it still has not converged by then, a warning reports this and the result at the largest value tried is returned anyway (not an error), so this never silently loops for an unbounded amount of time. Pass an explicit number for either argument to skip auto-tuning it and use a fixed value instead (as in previous package versions), or call checkBandcountConvergence() directly for more control over the tolerance and doubling count. See also vignette("BJM-intro", package = "BJM") for a worked example.

Value

An object of class "predictLongitudinal.BJM", a named list with elements:

Y_predict

A vector, one entry per at-risk patient (named by patient id), giving the MAP (most likely) predicted value of biomarker bio_i at prediction_time + horizon. For an ordinal bio_i, this is an integer category code (see Details), not a raw value.

Y_density

A probability matrix whose rows correspond to the candidate biomarker values in Y_all and whose columns correspond to individual patients; each entry is the dynamically predicted density of the biomarker taking that value. For an ordinal bio_i, each row is instead that category's predicted probability.

Y_all

The grid of candidate biomarker values used to build Y_density. For an ordinal bio_i, this is the vector of integer category codes 1:K, carrying a "category_labels" attribute with the matching level-label strings (see Details).

Predict multiple future biomarker values from fitted sub-models

Description

Internal multi-biomarker engine behind predictLongitudinal – call predictLongitudinal() directly instead (it dispatches here automatically when bio_i names more than one biomarker, or is left NULL).

Batch counterpart to dynamicPredictionBio: predicts every requested biomarker in bio_i from the same fitted longitudinal/survival sub-models and the same data_predict_all, but computes the bio_i-independent part of the pipeline (restricting to at-risk patients, the survival-side integration grid, and the denominator conditional density – see compute_bio_shared_step()) only once and reuses it across every biomarker, instead of recomputing it once per biomarker the way calling dynamicPredictionBio() in a loop would. This matters most under the Gaussian-copula path (see longitudinalSubCopula), where that denominator involves mvtnorm::pmvnorm() Monte-Carlo evaluations that are otherwise the dominant cost of a multi-biomarker prediction run.

bandcount2 (controlling the shared survival-integration grid) and bandcount3 (controlling each biomarker's own candidate-value grid) are auto-tuned separately when left at their default "auto": if bandcount2 = "auto", it is resolved once, using a single representative biomarker (the first one in bio_i), via the same doubling-until-stable check dynamicPredictionBio uses; the shared step is then built once at that resolved value. If bandcount3 = "auto", it is then resolved independently for every biomarker in bio_i (their candidate-value grids need not converge at the same resolution), reusing the once-computed shared step for every doubling round rather than rebuilding it – except for any ordinal biomarker, for which bandcount3 tuning is always skipped (its candidate grid is fixed at its category count; see dynamicPredictionBio) and "bandcount3" is recorded as NA for that biomarker. See Details in dynamicPredictionBio for the general auto-bandcount rationale.

Usage

dynamicPredictionBioAll(
  bio_i = NULL,
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  horizon,
  time_variable,
  survival_variable_all,
  survival_trans_function,
  bandcount2 = "auto",
  bandcount3 = "auto"
)

Arguments

Value

A named list of "predictLongitudinal.BJM" objects (see dynamicPredictionBio), one per requested biomarker, named by that biomarker's response-variable name; with attributes "bandcount2" (the single resolved/used bandcount2) and "bandcount3" (a named numeric vector of the resolved/used bandcount3 for each biomarker, or NA for an ordinal biomarker). Classed "predictLongitudinalAll.BJM".

Interval edges holding equal conditional survival mass

Description

Helper for prepare_infinity_grid(): n_intervals + 1 edges from prediction_time to upper_bound such that the average, over at-risk patients, of S(t \mid x) / S(s \mid x) drops by the same amount across every interval. Falls back to equally spaced edges if no patient has a usable linear predictor.

Usage

equal_mass_edges(
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  upper_bound,
  n_intervals
)

Value

A strictly increasing numeric vector of length n_intervals + 1.

Posterior distribution of the event time at a landmark

Description

For the model-based performance measures: each patient's distribution of T (and event type) given their history up to the landmark s and T > V, their last visit – the same quantity predictRisk() integrates – tabulated on the cells of predictRisk()'s grids: (V, s] (interval_gap_grid()) and (s, upper] (prepare_infinity_grid()).

Usage

event_time_posterior(
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  s,
  time_variable,
  survival_variable_all,
  survival_trans_function,
  n_inf = 100,
  n_gap = 20
)

Arguments

Value

A named list (patient ids) of lists with edges (cell boundaries, from V) and mass (cells x event types, summing to 1).

Vertical line at the end of a patient's follow-up

Description

Helper for riskPlot(): a vertical line at the patient's survival time, labelled in a legend – solid, coloured by event type (as the risk curves: black for the first, red for the second), for an observed event, and dashed grey for a censored follow-up. Nothing is drawn when the time or the status is unknown, as for a new patient or an interval-censored fit. The line used to be drawn the same way (red, with no legend) for a censored patient as for an event.

Usage

event_time_vline(
  data_predict_all_pre,
  survival_fit_all,
  event_type_variable = NULL
)

Arguments

Value

A list of ggplot2 components (possibly empty).

Exponentiate log densities after a per-patient shift

Description

Exponentiate log densities after a per-patient shift

Usage

exp_shifted(log_density, shift)

Arguments

Value

exp(log_density - shift), column-wise.

Slope of the extrapolated baseline cumulative hazard

Description

The constant hazard cumulative_baseline_at() uses past the last tabulated time: the slope of the least-squares line through the whole table, or 0 if that slope is negative or undefined (a table with a single time), so the cumulative hazard never decreases.

Usage

extrapolation_slope(cum_basehaz)

Arguments

Value

A single non-negative number.

Write event times (and their transforms) into longitudinal data

Description

Write event times (and their transforms) into longitudinal data

Usage

fill_event_time(
  data_fit_all,
  id,
  T_values,
  event_time,
  survival_variable_all = NULL,
  survival_trans_function = NULL
)

Fit the backward joint model to interval-censored event times

Description

Experimental. With interval censoring nobody's event time T is known exactly – only L < T \le R – so the longitudinal sub-model f(Y \mid T), which has T as a covariate, cannot be fit on the subjects with an observed event the way it is for right-censored data. fitIntervalBJM() fills in T by stochastic EM, then by multiple imputation:

1. fit the marginal survival model to the intervals (survivalSub() with Surv(L, R, type = "interval2")); it does not depend on T, so it is fit once;

2. start every interval-censored T at its interval midpoint and fit longitudinalSub();

3. draw each subject's T from f(T \mid L < T \le R, Y, X) \propto f(Y \mid T, X) f(T \mid X) on (L, R] (imputeEventTime()), refit longitudinalSub(), and repeat for n_burnin iterations;

4. keep iterating, saving every n_thin-th draw and its longitudinalSub() fit until n_imputations are saved.

Step 3 uses the backward model's own factorization: f(Y \mid T) on a grid of T is exactly what dynamic prediction already computes (conditionalYT()), so each draw costs about one predictRisk() call per subject.

For a subject whose event was detected, only the longitudinal measurements up to the last visit without the event (L) are used: the one at the detecting visit R (and any later one) is taken after the event, which f(Y \mid T) does not describe, and they are left out with a message. As in the rest of the interval-censoring support, the biomarker visits are taken to be the inspection visits, so no measurement lies between L and R.

By default only subjects whose event was observed (finite R) enter the longitudinal fit, mirroring the complete-case fit used for right-censored data, which is valid when censoring is independent of (T, Y) given the covariates. include_right_censored = TRUE also imputes T \in (L, \infty) for right-censored subjects and uses them; their draws then rely on the survival model's extrapolated tail.

The saved draws come from one chain, so successive draws are correlated; saving only every n_thin-th iteration makes them closer to independent. They are drawn with the model parameters at their point estimates, not from their posterior, so the between-imputation variance in pooled can still be somewhat understated.

For prediction, long_fit_all combines the imputations: its fixed effects are the Rubin's-rules pooled estimates and its residual and random-effects variances the averages over the imputations (see poolLongitudinalSub), so a prediction does not depend on which single imputation happened to come first. With competing risks the event-type model of survival_fit_all likewise has its coefficients averaged over the imputations.

Usage

fitIntervalBJM(
  data_survival,
  data_fit_all,
  form_marginal_surv,
  event_time,
  long_sub_fixed,
  long_sub_random,
  time_variable,
  survival_variable_all = NULL,
  survival_trans_function = NULL,
  form_conditional_cr = NULL,
  n_burnin = 10,
  n_imputations = 5,
  n_thin = 5,
  n_grid = 40,
  baseline = c("spline", "piecewise"),
  df = NULL,
  include_right_censored = FALSE,
  seed = NULL
)

Arguments

Value

An object of class "fitIntervalBJM", a list with survival_fit_all (with competing risks, its event-type model has the coefficients averaged over the imputations, while the rest of that glm object, e.g. its standard errors, is the first imputation's fit; survival_fit_all_list has one fit per imputation); long_fit_all_list (one longitudinalSub() fit per imputation); long_fit_all (the imputations combined, see Description, for use with predictRisk() and friends); pooled (poolLongitudinalSub of the fits, when n_imputations > 1); imputed_T (a subjects x imputations matrix); trace (each iteration's fixed effects, to check the chain has settled); and settings (the iteration settings). print() shows a short summary.

Examples

# simulated study: visits every 1-3 years, an event known only to lie
# between the last visit without it (L) and the first visit with it (R)
set.seed(1)
n <- 60
x <- rnorm(n)
event_time <- (-log(runif(n)) / (0.05 * exp(0.5 * x)))^(1 / 1.5)
surv <- data.frame(id = seq_len(n), x = x, L = NA_real_, R = NA_real_)
long <- list()
for (i in seq_len(n)) {
  visits <- c(0, cumsum(runif(30, 1, 3)))
  visits <- visits[visits < 12]
  if (event_time[i] > max(visits)) {
    surv$L[i] <- max(visits)
    seen <- visits
  } else {
    k <- findInterval(event_time[i], visits)
    surv$L[i] <- visits[k]
    surv$R[i] <- visits[k + 1]
    seen <- visits[visits <= event_time[i]]
  }
  long[[i]] <- data.frame(id = i, year = seen, x = x[i],
                          y = 1 + 0.3 * seen - 0.2 * event_time[i] + 0.4 * x[i] +
                            rnorm(1, 0, 0.5) + rnorm(length(seen), 0, 0.3))
}
long <- do.call(rbind, long)

fit <- fitIntervalBJM(surv, long, Surv(L, R, type = "interval2") ~ x,
                      event_time = "Tev", long_sub_fixed = y ~ year + Tev + x,
                      long_sub_random = ~ year | id, time_variable = "year",
                      n_burnin = 2, n_imputations = 2, n_thin = 1, n_grid = 15,
                      seed = 1)
fit$pooled

# the fit plugs into the prediction functions; Tev is unknown for a new patient
history <- long[long$id == 1 & long$year <= 2, ]
history$Tev <- NA_real_
predictRisk(history, fit$long_fit_all, fit$survival_fit_all,
            prediction_time = 2, horizon = 2, time_variable = "year",
            survival_variable_all = NULL, survival_trans_function = NULL,
            bandcount1 = 10, bandcount2 = 20)

Event-type model for an interval-censored outcome

Description

The competing-risks part of survivalSub(): a logistic model for the event type among subjects whose event was observed (finite R). The event type is taken as known once the event is detected, but the event time is not, so a model that uses it needs the event_time column filled for those subjects (fitIntervalBJM() fills it with imputed times and refits this model every iteration).

Usage

fit_event_type_interval(
  data,
  form_marginal_surv,
  form_conditional_cr,
  event_time
)

Fit the marginal model for a single ordinal/categorical biomarker

Description

Fits a cumulative link mixed model (probit link) via ordinal::clmm(), used as the marginal model for a categorical biomarker in the Gaussian-copula extension of longitudinalSub() (see longitudinalSubCopula()). The probit link is used specifically (rather than the more common logit) because it integrates directly with the latent-Gaussian-score representation the copula extension relies on: the underlying continuous latent variable implied by a probit cumulative link model is, by construction, Gaussian, which is exactly the representation needed to fold a categorical biomarker into the same random-effects covariance structure as the continuous biomarkers.

This follows the classical two-stage IFM (Inference Functions for Margins; Joe 2005) strategy already used throughout BJM: fit each biomarker's marginal model separately first (here, exactly as nlme::lme() does for continuous biomarkers), then re-estimate the joint random-effects covariance across all biomarkers afterwards (see longitudinalSubVarCopula()).

Usage

fit_marginal_ordinal(fixed_formula, random_formula, data)

Value

A fitted ordinal::clmm object.

A biomarker's fixed-effects mean at every integration grid point

Description

The longitudinal mean X(t)\beta depends on the event time t only through the survival-time column and its survival_variable_all transformations. This evaluates it for all rows of data_i (all patients) at all grid points l_i with one model.frame()/model.matrix() call on the rows stacked once per grid point, instead of one call per patient and grid point, whose fixed overhead dominated the run time of the prediction functions. The terms, factor levels and contrasts are those of the fit, as in predict.lme(), so the design matrix is the one the model was fit with, including data-dependent bases such as poly() and splines::ns(). Transformed terms such as log(years) are evaluated at each grid point. Rows are processed in chunks of at most max_rows stacked rows.

Usage

fixed_mean_grid(
  data_i,
  lfit_i,
  xlev_i,
  l_i,
  survival_variable,
  survival_variable_all,
  survival_trans_function,
  fixed_values = list(),
  beta = lfit_i$coefficients$fixed,
  max_rows = 2e+05
)

Arguments

Value

A nrow(data_i) x length(l_i) matrix; a row with a missing covariate is NA.

Gauss-Hermite quadrature rule for the standard normal

Description

Nodes and weights (summing to 1) of the k-point rule for E f(X), X \sim N(0, 1), by the Golub-Welsch eigenvalue method.

Usage

gauss_hermite_normal(k)

Arguments

Value

A list with x (nodes) and w (weights).

Gaussian log density from a Cholesky factor

Description

Zero-mean multivariate normal log density of each column of E, with covariance R^T R given by its upper Cholesky factor R.

Usage

gaussian_logdens_chol(R, E)

What a subject's other visits say about each biomarker

Description

Inputs for the deep imputation models (history = TRUE in imputeLongitudinal): without them each visit is completed from that visit alone, so an imputed value does not carry the subject's own level, which shrinks the estimated variation between subjects. For each visit and biomarker: the mean of the subject's other observed values of the biomarker, how many there are (log1p), and the value at the nearest other observed visit and its distance in time (log1p). The visit's own value is left out, so the model cannot simply copy it while training. A subject with no other observed value gets the biomarker's overall mean, a count of 0 and the largest distance.

Usage

history_features(aligned, id, time_variable, resp_mat, mask)

Value

A matrix with four columns per biomarker, one row per visit.

Proportional hazards model for interval-censored event times

Description

Internal fitter behind survivalSub() when the outcome is interval-censored (Surv(L, R, type = "interval2")). coxph cannot fit such data, so the marginal model is a proportional hazards model H(t \mid x) = H_0(t) e^{x\beta} with a parametric baseline, fit by maximum likelihood. Each subject contributes S(L) - S(R) for L < T \le R, S(L) when right-censored (R = \infty), 1 - S(R) when left-censored, and h(T) S(T) for an exactly observed T. Two baselines:

"spline" (default)

Royston–Parmar: \log H_0(t) is a natural cubic spline in \log t with df degrees of freedom (df - 1 internal knots at quantiles of the log interval endpoints, boundary knots at their extremes). df = 1 is a Weibull model.

"piecewise"

a piecewise-constant hazard on df pieces cut at quantiles of the interval endpoints.

The spline is the default because interval-censored data say little about the shape of the hazard within a visit interval, and that shape is exactly what imputeEventTime() draws T from: the piecewise-constant hazard is flat there, which in simulations biased the longitudinal sub-model fit by fitIntervalBJM() when the true hazard was increasing, while the smooth spline carries the trend across neighbouring intervals into each one.

The baseline cumulative hazard is then tabulated on a fine grid in the same hazard/time layout as basehaz, so the prediction code can use it exactly as it uses the Breslow estimate of a right-censored fit (see survival_accessors).

Usage

icphFit(
  formula,
  data,
  event_time,
  baseline = c("spline", "piecewise"),
  df = 3,
  lower_tail = c("spline", "weibull")
)

Arguments

Value

An object of class "icph.BJM".

Baselines for icphFit()

Description

Each returns a list with n_par, cumhaz(t, g) (baseline cumulative hazard, 0 at t <= 0, Inf at Inf), loghaz(t, g) (log baseline hazard, for exactly observed times), monotone(g) (whether cumhaz is increasing), start(rate0) (parameters of a constant hazard rate0) and knot_times.

Usage

rp_baseline(endpoints, df, tail_slope = NULL)

piecewise_baseline(endpoints, df)

Arguments

Functions

• rp_baseline(): Royston–Parmar spline for \log H_0 in \log t. Parameters: intercept, then one per basis column; the basis columns are scaled by their standard deviation over endpoints so the optimizer sees comparable scales.

• piecewise_baseline(): Piecewise-constant hazard; parameters are the log hazards of the pieces.

Linear predictor of an interval-censored PH fit

Description

Linear predictor of an interval-censored PH fit

Usage

icph_lp(fit, newdata)

Arguments

Stratum of each patient for a stratified interval-censored fit

Description

Stratum of each patient for a stratified interval-censored fit

Usage

icph_strata(fit, newdata)

Arguments

Value

Character labels in the fit's strata() format, or NULL when the fit is not stratified.

Covariates of the imputation model

Description

The right-hand-side variables of long_sub_fixed other than the subject id, including the visit time time_variable when a formula uses it. The visit time is also the key that aligns rows across biomarkers, and it used to be left out of the model along with the id, so values were imputed without knowing when the measurement was taken (pulling a trending biomarker towards its overall level).

Usage

imputation_covariates(covariate_vars, id)

Arguments

Value

Character vector of covariate names.

Draw interval-censored event times given the biomarker history

Description

Experimental. One draw of each subject's event time from f(T \mid L < T \le R, Y, X) \propto f(Y \mid T, X) f(T \mid X): f(Y \mid T) from the longitudinal sub-model on a grid of n_grid intervals tiling (L, R] (conditionalYT()), times the survival model's probability of each grid interval (marginalT()); an interval is drawn with those weights and T uniformly within it. A right-censored subject (R = \infty) uses the same equal-probability grid from L out to the survival model's tail as predictRisk(). A subject with L = R keeps T = L. Used by fitIntervalBJM.

Usage

imputeEventTime(
  data_fit_all,
  long_fit_all,
  survival_fit_all,
  bounds,
  time_variable,
  survival_variable_all = NULL,
  survival_trans_function = NULL,
  n_grid = 40,
  n_draws = 1
)

Arguments

Value

A named numeric vector of draws, one per row of bounds; with n_draws > 1, a matrix with one column per draw.

Impute missing longitudinal biomarker values

Description

Fills in missing biomarker values before longitudinalSub() is fit, assuming they are missing at random (MAR): whether a value is missing depends only on what is observed.

It is often not needed. longitudinalSub() fits each biomarker's linear mixed model by maximum likelihood on all of that biomarker's observed values, which is already valid under MAR; only the joint covariance of the random effects across biomarkers, Sigma_fit, is estimated from the subjects with some value of every biomarker. In simulations on pbc3 (see the "Missing data imputation" article on the package website), no imputation method made the predictions closer to those from the complete data than this complete-case analysis when visits were missing values at random; with half of the subjects never measured for a biomarker, multiple imputation with method = "mixed" was slightly closer. Imputation can also help when the missingness depends on variables outside the model, if they are added to the formulas used here.

Three methods are available, selected with method:

• "mixed" (default) draws each subject's missing values from the multivariate linear mixed model longitudinalSub() fits, given everything observed for that subject – all visits and all biomarkers – with the model's parameters drawn from an approximation to their posterior for each completion (proper multiple imputation; see mixed_model_imputations()). It needs no extra package, keeps the variation between subjects, and was the only method in the simulations that left the random-effects variances and the event-time coefficients unbiased.

• "miwae" fits a MIWAE (Missing data Importance-Weighted AutoEncoder; Mattei & Frellsen, 2019) – a variational autoencoder whose training objective only scores reconstruction on the observed entries – jointly across all supplied biomarkers, and uses self-normalized importance resampling from the fitted decoder to draw one or more plausible completions of the missing cells.

• "diffusion" fits a conditional denoising diffusion probabilistic model (DDPM; Ho, Jain & Abbeel, 2020) over the same per-row biomarker vector, trained as in CSDI (Tashiro et al., 2021): in every mini-batch, each row's observed entries are split at random into targets and conditions, and the network learns to denoise the targets given the conditions' values, the covariates and the observation mask, so the loss only ever involves observed entries. Missing cells are then filled by running the reverse diffusion with the missing entries as targets and all observed ones as conditions; observed entries are returned unchanged.

The two deep generative methods need the torch package. They complete each visit from that visit's covariates and observed biomarkers (and, with history = TRUE, a summary of the subject's other visits), and in the simulations they shrank the estimated variation between subjects, which every prediction uses; with history = TRUE they also biased the event-time coefficients when subjects lacked a biomarker entirely. The completed data.frame(s) of any method can be passed directly to longitudinalSub() in place of the original data_fit_all.

Usage

imputeLongitudinal(
  data_fit_all,
  long_sub_fixed,
  long_sub_random,
  time_variable,
  n_imputations = 5,
  latent_dim = 8,
  hidden_units = c(64, 32),
  epochs = 300,
  importance_samples = 20,
  diffusion_steps = 100,
  method = c("mixed", "miwae", "diffusion"),
  impute = c("single", "multiple"),
  history = FALSE,
  seed = NULL
)

Arguments

Details

Only the biomarker (response) columns named on the left-hand side of long_sub_fixed are imputed, and they must be numeric (continuous): an ordinal biomarker cannot be imputed with this function, and gives an error. Every covariate referenced on the right-hand side of long_sub_fixed/long_sub_random is assumed fully observed (an error is raised if any covariate has missing values). The deep generative methods model p(biomarkers | covariates, missingness mask) row by row (one row per subject-visit; the covariates include the visit time time_variable when long_sub_fixed uses it), and use long_sub_random only for the id variable; "mixed" uses the full random-effects structure.

Rows are aligned across biomarkers by (id, time_variable): when data_fit_all is a list of different data.frames with different visit schedules per biomarker (rather than one data.frame shared by every biomarker), the union of every (id, time_variable) pair seen in any of them is used, so a visit at which one biomarker was measured but another was not still gets a value imputed for the missing one.

Because the deep generative methods need enough data to fit reliably; the diagnostics element of the return value also reports a classical single regression imputation (each biomarker regressed on the covariates, fit on its observed rows) as a lightweight baseline, so the two can be compared on a given dataset before trusting the deep model's completions.

Value

A named list with elements:

data_fit_all

A list with one completed data frame per biomarker (also when data_fit_all was a single data frame), with every missing biomarker cell filled. Each holds one row per (id, time_variable) pair seen for any biomarker, and only the columns its model uses: the id, time_variable, the biomarker and the covariates of its long_sub_fixed formula – other columns of the input are not carried over. Ready to pass straight to longitudinalSub with the same formulas.

data_fit_all_list

When impute = "multiple", a list of n_imputations completed versions of data_fit_all, each from its own draw; NULL when impute = "single".

diagnostics

A named list, one element per biomarker (its deep_generative entry describes the completions of whichever method was used), each with the number of missing cells and observed-vs-imputed mean/sd for both the deep-generative and classical-baseline completions.

References

Mattei, P.-A. and Frellsen, J. (2019). MIWAE: Deep Generative Modelling and Imputation of Incomplete Data Sets. Proceedings of the 36th International Conference on Machine Learning, PMLR 97:4413-4423.

Ho, J., Jain, A., and Abbeel, P. (2020). Denoising Diffusion Probabilistic Models. Advances in Neural Information Processing Systems, 33:6840-6851.

Tashiro, Y., Song, J., Song, Y., and Ermon, S. (2021). CSDI: Conditional Score-based Diffusion Models for Probabilistic Time Series Imputation. Advances in Neural Information Processing Systems, 34:24804-24816.

Examples

data(pbc3)
data_fit_all <- pbc3[pbc3$status3 == 1, ]

# pbc3 has no missing serBilir/albumin values; delete some at random
# (MCAR, a special case of MAR) to have something to impute
set.seed(1)
n <- nrow(data_fit_all)
data_fit_all$serBilir[sample.int(n, floor(0.1 * n))] <- NA
data_fit_all$albumin[sample.int(n, floor(0.1 * n))] <- NA

long_sub_fixed <- list(
  "long1" = serBilir ~ year + age + sex + years,
  "long2" = albumin ~ year + age + sex + years)
long_sub_random <- list("long1" = ~ year | id, "long2" = ~ year | id)

# multiple imputation from the mixed model, pooled with Rubin's rules
imputed <- imputeLongitudinal(data_fit_all, long_sub_fixed, long_sub_random,
                              time_variable = "year", impute = "multiple", seed = 1)
fits <- lapply(imputed$data_fit_all_list, longitudinalSub,
               long_sub_fixed = long_sub_fixed, long_sub_random = long_sub_random)
pooled <- poolLongitudinalSub(fits)


# a deep generative method instead (needs torch)
if (requireNamespace("torch", quietly = TRUE) && torch::torch_is_installed()) {
  imputed_miwae <- imputeLongitudinal(data_fit_all, long_sub_fixed, long_sub_random,
                                      time_variable = "year", method = "miwae",
                                      epochs = 50, seed = 1)
}

Upper limit of the prediction-to-infinity integration grid

Description

Shared helper for predictRisk and dynamicPredictionBio. The denominator of a dynamic prediction integrates over every event time after prediction_time, out to infinity; the grid has to stop somewhere, and it should stop where the probability left beyond it is negligible. So the upper limit is the earliest time by which every at-risk patient's model-based conditional survival S(t \mid x) / S(s \mid x), s = prediction_time, has dropped below tail_prob, using the same (linearly extrapolated, per-stratum) baseline cumulative hazard as marginalT(). It never falls below min_upper.

This replaces twice the largest observed survival time of the patients being predicted, which used each patient's own future outcome (not available at prediction_time), stopped the integral early – and inflated the risk – for a patient whose event came soon after prediction_time, and failed when that time was missing.

Beyond the last time in the data survivalSub() was fit on, the baseline hazard is an extrapolation; if it is so flat that the tail criterion is not met by max_multiple times that last time, the bound is capped there, with a warning if some patient's conditional survival at the cap is still above warn_prob.

Usage

integration_upper_bound(
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  min_upper = prediction_time,
  tail_prob = 1e-04,
  max_multiple = 20,
  warn_prob = 0.01
)

Arguments

Value

A single number.

Spline interpolation of a log box probability over candidate values

Description

For copula_logdens_grid(): a cubic spline (end conditions from the end points' own cubic, not zero curvature, which a log probability, curving like a parabola in its tails, does not have) through the finite log probabilities at the exactly computed values x, evaluated at x_out; outside the range of the finite ones (where the probability underflowed to 0) the log probability is -Inf.

Usage

interpolate_log_box(x, y, x_out)

Value

A vector, one value per element of x_out.

Interval endpoints from an interval-type Surv object

Description

L < T \le R for every subject: R = \infty when right-censored, L = 0 when left-censored, L = R when T is observed exactly.

Usage

interval_bounds(y)

Arguments

Value

A list with numeric vectors L and R.

Integration grid between the last negative visit and the prediction time

Description

Interval-censored fits only. With interval censoring an event is only detected at the next visit, so at prediction_time s a patient is known to be event-free only up to their last visit V \le s (taken to be the last time in their biomarker history up to s), not up to s itself. Conditioning on T > V instead of T > s adds the stretch (V, s] to the denominator of every prediction. This builds, for each patient, n equal intervals tiling their own (V, s] – the stretch differs per patient, so the grid does too – and the survival model's probability of each.

Returns NULL for a right-censored fit (an event is observed when it happens, so V = s) and when every patient was seen at s, so callers keep their right-censored code path unchanged.

Usage

interval_gap_grid(
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  time_variable,
  n
)

Arguments

Value

NULL, or a list with last_visit (V, one per patient), mids and S (n x patients matrices of interval midpoints and probabilities; zero probability for a patient with V = s), and eval(fun), which calls fun(data_j, l_i) once per patient with a gap (data_j the patient's rows, l_i their midpoints) and binds the resulting one-column results into an n x patients matrix (filled with fill, -Inf by default for log densities, for patients without a gap). fun may also return a (nested) list of such columns, in which case eval() returns the matching list of matrices.

Is the outcome of a survival formula interval-censored?

Description

Is the outcome of a survival formula interval-censored?

Usage

is_interval_surv(form_marginal_surv, data)

Predicted risks and observed outcomes at each landmark time

Description

Shared by performancePlot() and calibrationPlot(): validates their common arguments and, for each landmark s, calls predictRisk() on the subjects event-free at s, with their history up to s and their outcome hidden.

Usage

landmark_predictions(
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  horizon,
  time_variable,
  survival_variable_all,
  survival_trans_function,
  bandcount1,
  bandcount2,
  posterior = FALSE
)

Arguments

Value

A list with one element per landmark that had subjects to predict, each a list with landmark, outcome (see performance_outcome(), rows matching the predictions) and risks (a list with one vector of predicted risks per event type).

Rank-one factor of the ordinal latent scores' covariance

Description

Writes cov as l l^T + I when cov - I has rank at most one: a single ordinal random effect, as seen through the patient's ordinal rows. Returns NULL (use pmvnorm()) otherwise, including when cov - I is not positive semi-definite (the latent residual variance is then not 1).

Usage

latent_factor(cov)

Value

The vector l (all zero for rank zero), or NULL.

Fit nlme::lme(), retrying with more optimizer iterations

Description

lme() with opt = "optim" stops after msMaxIter (default 50) optimizer iterations and then fails with "optim problem, convergence error code = 1". That happens when a random-effect variance is close to zero – e.g. a random slope in long_sub_random for data with little slope variation, more often with few measurements per subject. The call is evaluated as given first, so a fit that succeeds is exactly what it was before (call included); only if it fails to converge is it refit with msMaxIter = 1000, with a warning. If that fails too, the error suggests a simpler random-effects formula. Errors that are not about convergence are passed on unchanged, without a retry.

Usage

lme_with_retry(call_expr, label)

Arguments

Value

The lme fit.

Seed R's RNG until the calling function returns

Description

Calls set.seed(seed) and registers, in the calling function's frame, an exit handler that puts back the .Random.seed that existed before (or removes it if there was none), so a seed argument does not change the caller's random number stream.

Usage

local_r_seed(seed, frame = parent.frame())

Arguments

Fit a multivariate longitudinal sub-model

Description

Fits one linear mixed-effects model per longitudinal biomarker separately (via lme), keeping each biomarker's own fixed-effects and residual-variance estimates from that separate fit. The M separate fits are then combined into a single multivariate model by re-estimating the full random-effects variance-covariance matrix jointly across all M biomarkers, via an EM algorithm initialized from the block-diagonal covariance implied by the separate fits – so correlation between biomarkers' random effects is captured, even though the fixed effects and residual variances themselves are not re-estimated jointly and remain exactly what each biomarker's own lme() fit produced. This lets you assess the effect of predictor variables on several longitudinal outcomes at once while accounting for both population-level (fixed) effects and subject-level (random) variability, and for how the outcomes covary within a subject. The result is one of the two sub-models – together with survivalSub – that predictRisk/dynamicPredictionBio combine to produce dynamic risk/biomarker predictions.

Usage

longitudinalSub(
  data_fit_all,
  long_sub_fixed,
  long_sub_random,
  biomarker_type = NULL
)

Arguments

Value

An object of class "longitudinalSub.BJM", a named list with elements:

lfit

A list of fitted univariate linear mixed models, one per longitudinal outcome, each obtained via lme.

Sigma_fit

The estimated variance-covariance matrix of the random effects in the multivariate linear mixed model.

long_sub_fixed

The long_sub_fixed argument, as supplied.

long_sub_random

The long_sub_random argument, as supplied.

xlevels

A list, one element per longitudinal outcome, of the factor levels observed in the full training data for that outcome. Used internally at prediction time so that poly()/ splines::ns()/splines::bs()/factor() terms in long_sub_fixed reuse the basis/contrasts fit at training time instead of recomputing one from a small, patient-specific slice of data.

Examples

long_sub_fixed = list(
  "long1" = serBilir ~ year + age + sex +  (years) + (years) * year,  
  "long2" = prothrombin ~ year + age + sex + (years) + (years) * year,  
  "long3" = albumin ~ year + age + age * year + sex + (years) + (years) * year,  
  "long4" = alkaline ~ year + age + sex + (years) + (years) * year, 
  "long5" = SGOT ~ year + age + sex + (years) + (years) * year, 
  "long6" = platelets ~ year + age + sex + (years)  + (years) * year)

long_sub_random =list(
  "long1" =  ~ year| id,   
  "long2" =  ~ year| id,    
  "long3" =  ~ year| id,    
  "long4" =  ~ year| id,    
  "long5" =  ~ year| id,    
  "long6" =  ~ year| id)

# Complete case analysis
data_fit_all = list()
for(i in seq_len(length(long_sub_fixed))){
  data_fit_all[[i]] = pbc3[pbc3$status3 == 1, ]
}

# fitting longitudinal submodel
long_fit_all = longitudinalSub(data_fit_all, long_sub_fixed, long_sub_random)

# poly() in its default orthogonal mode, splines::ns()/bs(), and
# factor() are safe in a continuous biomarker's long_sub_fixed: the
# terms/xlevels/contrasts fit on the full training data are cached and
# reused when estimating Sigma_fit and at prediction time, instead of
# being recomputed from each patient's small per-prediction slice.
long_fit_poly = longitudinalSub(
  pbc3[pbc3$status3 == 1, ],
  serBilir ~ year + poly(age, 2) + factor(sex) + years,
  ~ year | id)

# A few more nonlinear-term styles. None of these would trigger the
# warning even in long_sub_random, because their basis doesn't depend on the surrounding
# data at all -- there's simply nothing that could disagree between the
# full training data and the small per-patient slice used at prediction
# time.
data_fit_extra_terms = pbc3[pbc3$status3 == 1, ]
# A transformed-time column can also be built once with survivalTrans()
# and merged onto the fitting data ahead of time -- by the time
# long_sub_fixed sees it, it is just an ordinary numeric column, same as
# any of the other terms below.
data_fit_extra_terms$Tyears1 = survivalTrans(c(1, 3, 5, 7))$survival_trans_function[[1]](
  data_fit_extra_terms$year)

long_fit_nonlinear = list(
  longitudinalSub(data_fit_extra_terms, serBilir ~ year + I(year^2) + age + sex,
                   ~ year | id),
  longitudinalSub(data_fit_extra_terms, serBilir ~ year + I(year^2) * age + sex,
                   ~ year | id),
  longitudinalSub(data_fit_extra_terms, serBilir ~ year + log(year + 1) + age + sex,
                   ~ year | id),
  longitudinalSub(data_fit_extra_terms, serBilir ~ year + sqrt(year) + age + sex,
                   ~ year | id),
  longitudinalSub(data_fit_extra_terms, serBilir ~ poly(year, 2, raw = TRUE) + age + sex,
                   ~ year | id),
  longitudinalSub(data_fit_extra_terms, serBilir ~ year + Tyears1 + age + sex,
                   ~ year | id)
)

Fit a multivariate longitudinal sub-model with mixed continuous/ordinal biomarkers (Gaussian-copula extension)

Description

Internal counterpart to longitudinalSubGaussian(), used by the longitudinalSub dispatcher whenever at least one biomarker is categorical/ordinal. Each biomarker is still fit on its own first, exactly following longitudinalSubGaussian()'s IFM strategy: continuous biomarkers via nlme::lme() (byte-identical code path to longitudinalSubGaussian()), ordinal biomarkers via fit_marginal_ordinal() (ordinal::clmm() with a probit link). The resulting per-biomarker fixed effects, thresholds, and random-effects covariance blocks seed the same D/beta/sigma2 bookkeeping longitudinalSubGaussian() builds, which is then handed to longitudinalSubVarCopula() (rather than longitudinalSubVar()) to jointly re-estimate D across all biomarkers, including the ordinal ones' latent Gaussian scores.

Usage

longitudinalSubCopula(
  data_fit_all,
  long_sub_fixed,
  long_sub_random,
  biomarker_type,
  M
)

Fit a multivariate longitudinal sub-model (all-continuous biomarkers)

Description

Internal workhorse for longitudinalSub when every biomarker is continuous. This is the original longitudinalSub implementation, kept verbatim and called with the caller's pristine, unmodified arguments, so that the all-continuous case behaves exactly as it always has – see longitudinalSub for the argument documentation, and longitudinalSubCopula() for the mixed continuous/ordinal extension used when at least one biomarker is categorical.

Usage

longitudinalSubGaussian(data_fit_all, long_sub_fixed, long_sub_random)

Variance-covariance matrix Reference: package "lmm" and package "joineRML" function mvlme.

Description

The EM is accelerated by SQUAREM (Varadhan and Roland, 2008), and stops once the relative change in D over an accelerated cycle (three EM updates) is at most tol.em, or after max.iter EM updates with a warning.

Usage

longitudinalSubVar(thetaLong, l, tol.em, verbose, max.iter = 1000)

Joint random-effects covariance with ordinal biomarkers (Gaussian copula)

Description

Re-estimates the shared random-effects covariance D across all biomarkers by maximum likelihood, holding each biomarker's fixed effects, residual variance and (for an ordinal biomarker) thresholds at their separate-fit values, as longitudinalSubVar() does for continuous ones. An ordinal observation is the category of a latent score Z = \eta + \epsilon, \epsilon \sim N(0, 1) (probit link, Albert & Chib), observed only through the thresholds that bracket it, so the posterior of the random effects given the data is not Gaussian and the EM's E[b b^T \mid y] has no closed form. It is computed by

• copula_D_aghq(), an EM with adaptive Gauss-Hermite quadrature over the ordinal biomarkers' random effects (the continuous biomarkers' random effects are integrated out exactly), when those number at most max_quadrature_dim; deterministic, with the usual convergence check;

• otherwise copula_D_saem(), a stochastic-approximation EM (SAEM) that draws the latent scores, whose cost does not grow with that dimension.

Both target the same maximum-likelihood estimate, but with the default number of iterations SAEM stopped short of it: for two ordinal biomarkers with a random intercept and slope each (four ordinal random effects), the random-effects variances were off by about 5\ quadrature (5^4 = 625 nodes per subject) removes; hence the default max_quadrature_dim = 4. The earlier version plugged in the mean of each latent score given its category, which ignored the score's own variance and underestimated the ordinal biomarkers' random-effects variances (by about 40\ for a random slope in simulations).

Usage

longitudinalSubVarCopula(
  thetaLong,
  l,
  biomarker_type,
  thresholds,
  tol.em = 1e-04,
  max.iter = 500,
  max_quadrature_dim = 4,
  verbose = FALSE
)

Arguments

Value

A list with the estimate D.

Marginal distribution of T

Description

Marginal distribution of T

Usage

marginalT(data_predict_all, long_fit_all, survival_fit_all, l_i)

Arguments

Value

Marginal density probability of the survival variable T is represented as a probability matrix. In this matrix, the rows (l_i) are aligned with specific time points, while the columns correspond to individual patients. Each entry in the matrix denotes the marginal density probability of survival for a given patient at a particular time point.

Point estimates, standard errors and complete-data df of one fit

Description

Helper for poolLongitudinalSub(). For a continuous biomarker's nlme::lme() fit: the fixed effects, their standard errors, and lme()'s own df per coefficient (fixDF$X) – a between-subject covariate such as age has far fewer than the within-subject time terms, so applying the intercept's df to every coefficient would overstate the df of between-subject effects. For an ordinal biomarker's ordinal::clmm() fit: the thresholds followed by the slopes, their Wald standard errors, and an infinite complete-data df (clmm()'s inference is asymptotic). Previously only lme() fits were handled, and pooling a fit with an ordinal biomarker failed.

Usage

marker_estimates(fit)

Arguments

Value

A list with named vectors estimate, std_error and dfcom.

Compare two prediction results' plain numeric-vector fields

Description

Shared comparison logic for checkBandcountConvergence() and the "auto" bandcount support in predictRisk()/ dynamicPredictionBio(): only plain numeric vectors (no dim) that have the same length in both results are compared. This naturally skips fields whose size is itself controlled by the bandcount being varied (e.g. dynamicPredictionBio()'s Y_density matrix and Y_all grid, whose resolution is exactly what bandcount3 sets), while still comparing the actual per-patient estimates derived from them (risk_prob_1/risk_prob_2, Y_predict).

Usage

max_relative_diff(result_a, result_b)

Value

A list with max (the largest relative change across all comparable fields, or NA if none were comparable) and by_field (a named numeric vector, one entry per comparable field).

Draw K importance-weighted latent samples and decode them

Description

Internal helper shared by train_miwae and miwae_impute: encodes each row once, draws K samples from q(z|x_o) (reusing the same encoder mean/logvar for all K, as in the original MIWAE), decodes all of them, and returns everything needed to compute the importance-weighted log-likelihood terms.

Usage

miwae_forward(net, X_cov, X_resp, X_mask, K)

Draw self-normalized-importance-resampled completions from a fitted MIWAE

Description

Internal helper for imputeLongitudinal: implements the MIWAE multiple-imputation procedure (Mattei & Frellsen, 2019, Section 3.2) – draw K decoder samples per row, importance-weight them by their (masked-to-observed-entries) likelihood, then resample one index per completion and use that draw's decoder mean plus Gaussian noise as the imputed value. Each completion uses its own K samples, so completions are independent draws. Already-observed cells are left untouched by the caller; this only needs to return sensible values for missing cells.

Usage

miwae_impute(net, X_cov, X_resp, X_mask, resp_mean, resp_sd, K, n_imputations)

mixed_density_prob_copula() for many means under one covariance

Description

Factors the parts of mixed_density_prob_copula() that depend only on the covariance – the inverse and log-determinant of the continuous block, and the conditional covariance of the ordinal block – once, and returns a function that evaluates the log density for many means, observed values and threshold bounds at once (one per column), e.g. a patient's density at every integration grid point.

The ordinal box probability is computed by adaptive Gauss-Hermite quadrature (box_prob_quadrature()) when the patient's ordinal rows involve a single ordinal random effect (e.g. one ordinal biomarker with a random intercept). Given the continuous values, the ordinal latent scores then have covariance l l^T + I, l l^T from the random effect and I from the latent residuals, so they are independent given the random effect, and the probability is a one-dimensional integral of a product of univariate normal probabilities. It used to be computed by mvtnorm::pmvnorm()'s Monte Carlo algorithm, which dominated the run time of every prediction with an ordinal biomarker, differed between runs (relative error around 10^{-4}), and needed set.seed() to be reproducible. With two or more ordinal random effects pmvnorm() is still used: a product quadrature grid then needs so many nodes to match its accuracy (a random intercept with a large variance makes the integrand lopsided) that it is no faster. It is also used when the covariance does not have this structure, and when box = "pmvnorm".

Usage

mixed_density_copula_prep(
  Sigma_all,
  cc_idx,
  oo_idx,
  box = c("quadrature", "pmvnorm")
)

Arguments

Value

A function (mu_full, y_all_vec, alpha_lower_vec, alpha_upper_vec) returning the log densities, with attribute "monte_carlo": whether pmvnorm() is used. mu_full is a vector or a matrix with one column per evaluation; the other three may be vectors (the same for every column) or matrices of the same shape.

Joint density x probability for a mixed continuous/ordinal stacked vector

Description

Given one patient's full stacked mean vector and covariance matrix (across all jointly-fit biomarkers' repeated measurements), computes the exact Gaussian density at the observed values for the continuous-marker rows, times the Gaussian-copula probability that the unobserved ordinal-marker latent scores fall in the box implied by their observed categories' cumulative-link thresholds – conditional on the continuous rows' observed values (via the standard multivariate-normal conditioning/Schur-complement formula). This is the quantity conditionalYT()/conditionalYDT() compute directly as a single Gaussian density when every marker is continuous (see the quadratic-form comment in R/conditionalYT.R); setting oo_idx = integer(0) here reduces to exactly that computation (the full (2*pi)^{-n/2} normalizing constant is kept, not dropped, even though it would cancel in predictRisk()'s own ratio – dynamicPredictionBio() instead divides this function's output by mvtnorm::dmvnorm()-normalized values from conditionalYTBioCopula()/conditionalYDTBioCopula(), so the constant must be correct in absolute terms, not just consistent within one call, for Y_density to remain a genuine, correctly-scaled conditional density – see R/conditionalYTBio.R's analogous all-continuous computation, which keeps its own constant for the same reason). Returns the log of that density/probability, so that it neither overflows nor underflows (see predictRisk() for how it is exponentiated); -Inf if the ordinal box probability is 0.

Usage

mixed_density_prob_copula(
  Sigma_all,
  mu_full,
  y_all_vec,
  alpha_lower_vec,
  alpha_upper_vec,
  cc_idx,
  oo_idx
)

Arguments

Multiple imputation from the multivariate linear mixed model

Description

method = "mixed" in imputeLongitudinal. Fits the model longitudinalSub fits – each biomarker's linear mixed model, with the random effects of all biomarkers jointly normal with covariance Sigma_fit – and draws each subject's missing values from their normal distribution given everything observed for that subject: all of its visits and all biomarkers. For each completion the model's parameters are first drawn from an approximation to their posterior, so that the completions also reflect the uncertainty in the fit (proper multiple imputation): the fixed effects from their estimated sampling distribution, each residual variance from a scaled inverse chi-squared distribution, and Sigma_fit from an inverse Wishart distribution with as many degrees of freedom as subjects.

Usage

mixed_model_imputations(
  aligned,
  data_fit_all,
  long_sub_fixed,
  long_sub_random,
  id,
  resp_mat,
  mask,
  n_imputations
)

Value

A list of n_imputations matrices like resp_mat, with the missing cells filled.

Variables a biomarker's sub-model uses

Description

The rows longitudinalSub() keeps for the joint EM must be the rows lme() was fit on, i.e. complete in every variable of both the fixed- and the random-effects formula (including the id). Only the fixed-effects variables used to be checked, so a missing value in a variable used only in long_sub_random left a row in the fixed-effects design but not in the random-effects one, and failed with "arguments imply differing number of rows".

Usage

model_vars(fixed_formula, random_formula)

Arguments

Value

A character vector of variable names.

Number of subjects the random-effects covariance is estimated from

Description

The EM update for Sigma_fit averages each subject's E[b b^T] over the subjects that are actually in the fit – those with at least one complete observation of every biomarker (yi has one element per such subject). It used to divide by the number of subjects in the first biomarker's raw data instead, which underestimated Sigma_fit (badly, since EM compounds it over iterations) whenever some subject was excluded, e.g. one never measured for some biomarker.

Usage

n_subjects(yi)

Arguments

Value

The number of subjects.

Convert an nlme-style random-effects formula to an lme4-style bar term

Description

longitudinalSub()'s long_sub_random argument uses nlme::lme()'s ~ terms | group formula convention. ordinal::clmm() (used to fit the marginal model for a categorical/ ordinal biomarker – see fit_marginal_ordinal()) instead expects random effects written as an lme4-style (terms | group) bar term embedded directly in the model formula. This helper translates one into the other so that the exact same long_sub_random argument drives both nlme::lme() (continuous biomarkers) and ordinal::clmm() (ordinal biomarkers), with no second random-effects argument for the user to keep in sync.

Usage

nlme_random_to_lme4_bars(random_formula)

Value

A length-1 character string, e.g. "(year | id)".

Nonparametric estimate of the event-time distribution from interval data

Description

Self-consistency (EM) estimate in the spirit of Turnbull (1976) for interval-censored event times (L, R], left-truncated at V (each subject is only in the sample because T > V), with an optional known event type per detected event (the competing-risks estimate of Hudgens, Satten & Longini, 2001, Biometrics 57:74). Mass sits on the cells between consecutive distinct time points (all L, R, V and extra); within an observed interval it is not identified, and the EM's choice there is what the estimate reports.

Usage

npmle_interval(
  L,
  R,
  V,
  cause = NULL,
  extra = NULL,
  tol = 1e-09,
  max_iter = 5000
)

Arguments

Value

A list with edges and mass (cells x event types, the last cell running to Inf).

Latent-score means of an ordinal biomarker at every candidate event time

Description

X(l)\beta for an ordinal biomarker's rows, with the design built the way copula_mean_grid() builds it (the columns of lfit[[i]]$beta, which has no intercept: the thresholds play that role).

Usage

ordinal_mean_on_grid(
  data_i,
  i,
  long_fit_all,
  xlev_i,
  l_unique,
  survival_variable,
  survival_variable_all,
  survival_trans_function
)

Value

A matrix with one row per row of data_i and one column per element of l_unique.

Legend label for event type k in performancePlot()

Description

Legend label for event type k in performancePlot()

Usage

outcome_cause_label(survival_fit_all, k)

Outcome variables of the survival sub-model

Description

The survival-time variable of survival_fit_all and, with competing risks, its event-type variable (the response of form_conditional_cr). Both are unknown for a patient still at risk at prediction_time; the prediction functions integrate over them.

Usage

outcome_variables(survival_fit_all)

Arguments

Value

A character vector of variable names.

Apply a function to every patient, on one core or several

Description

Used by the copula density functions, whose cost is per patient. With getOption("BJM.n_cores") above 1 (set by the n_cores argument of predictRisk() and predictLongitudinal()) the patients are spread over that many cores with parallel::mclapply() (serially on Windows). Each patient is computed with its own fixed random-number seed, derived from its id (patient_seed()), and the caller's random-number state is restored afterwards, so the Monte Carlo box probabilities still used with two or more ordinal random effects give the same result serially, in parallel, and on every run, without set.seed().

Usage

patient_lapply(ids, FUN)

Arguments

Value

A list, one element per id.

Per-patient shift for exponentiating log densities

Description

The conditional-density helpers (conditionalYT() and relatives) return log densities, one column per patient. Every quantity a prediction needs is a ratio of sums of those densities within one patient, so any per-patient constant cancels: this returns, for each patient, the largest finite log density across all the matrices given, to subtract before exponentiating (see exp_shifted()). Without it the densities – whose scale is the determinant of a covariance matrix that grows with the number of observations and with the biomarkers' units – underflowed to 0 or overflowed to Inf.

Usage

patient_log_shift(...)

Arguments

Value

A numeric vector, one shift per patient (0 for a patient with no finite value).

One patient's stacked fixed-effects means at every grid point

Description

fixed_mean_grid() for each biomarker's rows of one patient (data_num_i_list), stacked in biomarker order as in build_conditional_design().

Usage

patient_mean_grid(
  data_num_i_list,
  long_fit_all,
  l_i,
  survival_variable,
  survival_variable_all,
  survival_trans_function,
  fixed_values = list()
)

Value

A matrix, one row per stacked observation and one column per grid point in l_i.

A fixed random-number seed for a patient

Description

Derived from the printed id, so that a patient gets the same seed whichever other patients are predicted alongside it.

Usage

patient_seed(id)

Value

An integer.

Mayo Clinic primary biliary cirrhosis data used as example code

Description

The dataset originates from the Mayo Clinic trial on primary biliary cirrhosis (PBC) of the liver, carried out from 1974 to 1984. It includes data from 424 PBC patients who were referred to the Mayo Clinic within this decade and met the eligibility requirements for a randomized placebo-controlled trial of D-penicillamine. However, only the initial 312 cases from the dataset were enrolled in the randomized trial. Thus, the dataset specifically pertains to these 312 patients, for whom the data is largely complete.

Usage

data(pbc3)

Format

A data frame with 1945 observations on the following 27 variables:

id

patients identifier; in total there are 312 patients.

years

number of years between registration and the earlier of death, transplantation, or study analysis time.

status

a factor with levels alive, transplanted and dead.

drug

a factor with levels placebo and D-penicil.

age

at registration in years.

sex

a factor with levels male and female.

year

number of years between enrollment and this visit date, remaining values on the line of data refer to this visit.

ascites

a factor with levels No and Yes.

hepatomegaly

a factor with levels No and Yes.

spiders

a factor with levels No and Yes.

edema

a factor with levels No edema (i.e. no edema and no diuretic therapy for edema), edema no diuretics (i.e. edema present without diuretics, or edema resolved by diuretics), and edema despite diuretics (i.e. edema despite diuretic therapy).

serBilir

serum bilirubin in mg/dl.

serChol

serum cholesterol in mg/dl.

albumin

albumin in mg/dl.

alkaline

alkaline phosphatase in U/liter.

SGOT

SGOT in U/ml.

platelets

platelets per cubic ml/1000.

prothrombin

prothrombin time in seconds.

histologic

histologic stage of disease.

status2

a numeric vector with the value 1 denoting if the patient was dead, and 0 if the patient was alive or transplanted.

status3

a numeric vector with the value 1 denoting if the patient was dead or transplanted, and 0 if the patient was alive.

status4

a numeric vector with the value 1 denoting if the patient was transplanted, and 0 if the patient was dead. Used for competing risks.

status5

a numeric vector with the value 2 if the patient was transplanted, 1 denoting if the patient was dead, and 0 if the patient was alive. Used for competing risks with censored.

Tyears1

a numeric vector with a transformed value of time-to-event outcome.

Tyears2

a numeric vector with a transformed value of time-to-event outcome.

Tyears3

a numeric vector with a transformed value of time-to-event outcome.

Tyears4

a numeric vector with a transformed value of time-to-event outcome.

Source

pbc.

References

Fleming T, Harrington D. Counting Processes and Survival Analysis. 1991; New York: Wiley.

Therneau T, Grambsch P. Modeling Survival Data: Extending the Cox Model. 2000; New York: Springer-Verlag.

Inverse of a covariance matrix that may be close to singular

Description

By the eigendecomposition, with eigenvalues floored at 1e-10 times the largest, so a covariance estimate approaching the boundary (e.g. a random-effect variance tending to 0) still gives a finite, positive-definite precision.

Usage

pd_inverse(V)

Arguments

Value

Its (regularized) inverse.

Plot predictive performance across landmark times

Description

Evaluates how well the backward joint model's dynamic risk predictions discriminate and are calibrated, at each of several landmark times. For each landmark s in prediction_time, every subject in data_predict_all still event-free at s gets a predicted risk from predictRisk, using only their longitudinal history up to s, of an event in (s, s + horizon]; this is compared with what actually happened, by

AUC

the time-dependent (cumulative/dynamic) area under the ROC curve: the probability that a subject with an event in the window has a higher predicted risk than one still event-free at s + horizon. 0.5 is no better than chance, 1 is perfect.

Brier score

the mean squared difference between the predicted risk and the observed event indicator in the window. Lower is better.

Subjects censored within the window are handled by inverse probability of censoring weighting (IPCW), with the censoring distribution among the subjects at risk at s estimated by Kaplan–Meier. With competing risks, both measures are computed separately for each event type: the cases are the subjects with that event type in the window, and for the Brier score a subject with the other event type in the window counts as not having the event.

Evaluating on the data the sub-models were fit on gives an optimistic (apparent) performance; pass held-out validation data in data_predict_all when available.

Usage

performancePlot(
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  horizon,
  time_variable,
  survival_variable_all,
  survival_trans_function,
  bandcount1 = "auto",
  bandcount2 = "auto",
  interval_method = c("model", "ipcw")
)

Arguments

Value

A ggplot object with one panel per measure, the landmark time on the horizontal axis. Its data element is a data.frame with one row per landmark, measure (and event type), and columns landmark, measure, cause, value, n_at_risk (subjects predicted) and n_cases (events of that type observed in the window).

Interval-censored fits (experimental)

The event is then only known to lie between the last negative visit L and the visit R that detected it. A subject is at risk at s when followed up past s without a detected event, and the predicted risk is that of predictRisk() (conditional on being event-free at the last visit). Following Yang, Rizopoulos, Newcomb and Erler (2026, Biometrical Journal 68:e70108), interval_method chooses how subjects whose interval straddles s or s + horizon are handled:

• "model" (default): every subject at risk counts, as a case or a control with the probability the fitted model gives their event time given their observed interval. The model only splits the interval observed for each subject, so in our simulations (visits every 1–3 years, 2-year window) both measures were within about 0.003 of their true values, also with a survival sub-model that left out a covariate. The model does grade itself, though: for a deliberately miscalibrated model the Brier score came out slightly too good (by about 0.003), so judge calibration with calibrationPlot, which does not use the model.

• "ipcw": only subjects whose status is certain are used – cases detected inside the window after a negative visit at or after s, controls with a negative visit at or after s + horizon – weighted by the inverse Kaplan–Meier probability of still being followed up. It does not use the model, but these weights do not account for how likely the status of a subject is to be certain, which depends on the visit schedule: in the same simulations only about 5 of 60 cases were certain, the Brier score came out at about 0.03 instead of 0.13, and the AUC was about 0.01 too high with three times the error of "model". With visits every 3–9 months the Brier score was still about 25% too low. Use it as a model-free check of the AUC.

Examples

data(pbc3)
survival_fit_all <- survivalSub(pbc3[!duplicated(pbc3$id), ],
                                Surv(years, status3) ~ age + sex, NULL)
long_sub_fixed <- list("long1" = serBilir ~ year + age + sex + years)
long_sub_random <- list("long1" = ~ year | id)
long_fit_all <- longitudinalSub(list(pbc3[pbc3$status3 == 1, ]),
                                long_sub_fixed, long_sub_random)
trans <- survivalTrans(c(1, 3, 5, 7))

# apparent performance on the fitting data, 2-year window
performancePlot(pbc3, long_fit_all, survival_fit_all,
  prediction_time = c(1, 3, 5), horizon = 2, time_variable = "year",
  trans$survival_variable_all, trans$survival_trans_function)

IPCW time-dependent AUC and Brier score at one landmark

Description

IPCW time-dependent AUC and Brier score at one landmark

Usage

performance_metrics(risk, time, status, cause, k, s, horizon)

Arguments

Value

A list with auc, brier, n_at_risk, n_cases.

IPCW AUC and Brier score for interval-censored outcomes at one landmark

Description

Uses only subjects whose status in (s, e] is certain: cases of type k (last negative visit L \ge s, event detected at R \le e), the other event type detected the same way (a non-case in the Brier score), and controls (L \ge e). They are weighted by the inverse probability of still being followed up – at R for an event, at e for a control – from the Kaplan–Meier estimate of the end of follow-up among the subjects at risk at s, as in Yang et al. (2026), Biometrical Journal 68:e70108. Subjects whose interval straddles s or e get weight 0.

Usage

performance_metrics_interval_ipcw(risk, L, R, followup, cause, k, s, horizon)

Arguments

Value

As performance_metrics(), plus n_uncertain.

Model-based AUC and Brier score for interval-censored outcomes

Description

Every subject at risk counts, as a case with probability p_case and a control with probability p_control (window_status_probs()): the expected squared error for the Brier score, and case/control pairs of different subjects weighted by p_i q_j for the AUC (Yang et al., 2026). These weights come from the model being evaluated, so the measures lean optimistic when that model is misspecified.

Usage

performance_metrics_interval_model(risk, p_case, p_control)

Arguments

Value

As performance_metrics(), with n_cases the expected number of cases.

One row per subject with the observed outcome, for performancePlot()

Description

One row per subject with the observed outcome, for performancePlot()

Usage

performance_outcome(data, survival_fit_all, id_variable)

Arguments

Value

A data.frame with columns id, time, status (1 = event, 0 = censored) and cause (1 or 2 for the event type matching risk_prob_1/risk_prob_2 of predictRisk(); NA without competing risks or when censored).

Observed outcome per subject for an interval-censored fit

Description

The interval-censored counterpart of performance_outcome(): each subject's event is known only to lie in (L, R] (R = Inf if no event was detected), with its type when detected.

Usage

performance_outcome_interval(data, survival_fit_all, id_variable)

Arguments

Value

A data.frame with columns id, L, R, followup (end of follow-up: R after an event, L otherwise) and cause (1 or 2, matching risk_prob_1/risk_prob_2; NA without competing risks or without an event).

Cumulative hazard of a piecewise-constant hazard

Description

Cumulative hazard of a piecewise-constant hazard

Usage

piecewise_cumhaz(t, lambda, cuts)

Arguments

Plot a fitted longitudinal sub-model

Description

Diagnostic plots for the per-biomarker mixed models and the shared random-effects covariance of a longitudinalSub fit, one panel per biomarker:

"residuals"

Standardized (Pearson) residuals against fitted values, with a LOESS smoother; a trend or a funnel shape suggests a misspecified mean or non-constant variance.

"qq"

Normal Q-Q plot of the standardized residuals.

"ranef"

Normal Q-Q plot of each predicted random effect (the model assumes they are normal).

"corr"

Heat map of the correlation matrix of the multivariate random-effects covariance Sigma_fit, i.e. how the subject-level deviations of the different biomarkers move together.

Residual plots are drawn for continuous biomarkers only; an ordinal biomarker (fit by ordinal::clmm()) has no residuals on the response scale and is left out of them with a message.

Usage

## S3 method for class 'longitudinalSub.BJM'
plot(x, which = c("residuals", "qq", "ranef", "corr"), ...)

Arguments

Value

A ggplot object.

Examples

data(pbc3)
long_sub_fixed <- list("long1" = serBilir ~ year + age + sex + years,
                       "long2" = albumin ~ year + age + sex + years)
long_sub_random <- list("long1" = ~ year | id, "long2" = ~ year | id)
data_fit <- pbc3[pbc3$status3 == 1, ]
long_fit_all <- longitudinalSub(list(data_fit, data_fit), long_sub_fixed, long_sub_random)

plot(long_fit_all)
plot(long_fit_all, which = "qq")
plot(long_fit_all, which = "ranef")
plot(long_fit_all, which = "corr")

Plot predicted biomarker distributions

Description

Plots the predicted distribution of the biomarker at prediction_time + horizon returned by predictLongitudinal, one curve per patient, with a dashed line at the patient's point prediction (the most likely value, Y_predict). For an ordinal biomarker, plots each patient's predicted category probabilities as bars instead.

Usage

## S3 method for class 'predictLongitudinal.BJM'
plot(x, subject = NULL, ...)

Arguments

Value

A ggplot object.

Examples

data(pbc3)
survival_fit_all <- survivalSub(pbc3[!duplicated(pbc3$id), ],
                                Surv(years, status3) ~ age + sex, NULL)
long_sub_fixed <- list("long1" = serBilir ~ year + age + sex + years)
long_sub_random <- list("long1" = ~ year | id)
long_fit_all <- longitudinalSub(list(pbc3[pbc3$status3 == 1, ]),
                                long_sub_fixed, long_sub_random)

data_predict <- pbc3[pbc3$id %in% c(2, 4, 5) & pbc3$year <= 3, ]
data_predict$years <- NA
trans <- survivalTrans(c(1, 3, 5, 7))
pred <- predictLongitudinal(data_predict, long_fit_all, survival_fit_all,
  prediction_time = 3, horizon = 1, time_variable = "year",
  trans$survival_variable_all, trans$survival_trans_function, bio_i = 1)
plot(pred)

Plot simulated futures

Description

Plots the draws of simulateTrajectory, pooling every draw of the patients chosen by id (all of them by default: one patient's draws picture that patient's predictive distribution, one draw each of a synthetic cohort pictures the cohort):

"trajectory"

One panel per biomarker. For a continuous biomarker, n_paths of the simulated trajectories as thin lines, with the median and the central level interval of the draws at each time, and the measurements conditioned on as points. For an ordinal biomarker, the share of the draws in each category at each time. With truncate = TRUE in simulateTrajectory() the summaries at a time describe the draws still event-free then.

"event"

The cumulative incidence of the event (of each event type, with competing risks) after prediction_time, computed from the drawn event times.

Usage

## S3 method for class 'simulateTrajectory.BJM'
plot(
  x,
  which = c("trajectory", "event"),
  id = NULL,
  n_paths = 30,
  level = 0.9,
  ...
)

Arguments

Value

A ggplot object.

Examples

data(pbc3)
survival_fit_all <- survivalSub(pbc3[!duplicated(pbc3$id), ],
                                Surv(years, status3) ~ age + sex, NULL)
long_fit_all <- longitudinalSub(list(pbc3[pbc3$status3 == 1, ], pbc3[pbc3$status3 == 1, ]),
                                list(serBilir ~ year + age + sex + years,
                                     albumin ~ year + age + sex + years),
                                list(~ year | id, ~ year | id))
history <- pbc3[pbc3$id == 2 & pbc3$year <= 3, ]
sims <- simulateTrajectory(history, long_fit_all, survival_fit_all,
                           prediction_time = 3, times = seq(3, 8, by = 0.5),
                           time_variable = "year", survival_variable_all = list(),
                           survival_trans_function = list(), n_sim = 200, seed = 1)
plot(sims)
plot(sims, which = "event")

Plot a fitted survival sub-model

Description

"forest"

Forest plot of the hazard ratios (with 95\ confidence intervals) of the Cox marginal survival model and, when the fit has competing risks, of the odds ratios (with 95\ intervals) of the logistic event-type model, on a log scale with a reference line at 1.

"basehaz"

The Cox model's baseline cumulative hazard (covariates at 0, not centered), one step curve per stratum.

Usage

## S3 method for class 'survivalSub.BJM'
plot(x, which = c("forest", "basehaz"), ...)

Arguments

Value

A ggplot object.

Examples

data(pbc3)
data_survival_fitting <- pbc3[!duplicated(pbc3$id), ]
survival_fit_all <- survivalSub(data_survival_fitting,
                                Surv(years, status3) ~ age + sex,
                                status4 ~ years + age + sex)
plot(survival_fit_all)
plot(survival_fit_all, which = "basehaz")

Pool longitudinal sub-model fits across multiple imputations with Rubin's rules

Description

When imputeLongitudinal is run with impute = "multiple", it returns n_imputations independently completed versions of data_fit_all in data_fit_all_list, instead of a single completed dataset filled with the across-draw mean (impute = "single"). Fitting longitudinalSub once per completion and combining the resulting fixed-effect estimates with Rubin's rules (Rubin, 1987) – rather than averaging into one completed dataset up front – is the standard way multiple imputation propagates the extra uncertainty from not knowing the true missing values into the final standard errors; impute = "single" does not do this (see the "Missing data imputation" article on the package website). poolLongitudinalSub() performs that combination step.

Usage

poolLongitudinalSub(long_fit_all_list)

Arguments

Details

For each fixed-effect coefficient, across the m fits in long_fit_all_list:

• the pooled estimate is the mean of the m per-completion estimates (\bar{Q});

• the pooled variance (T) adds the average within-imputation variance (\bar{U}, from each fit's own standard error) to the between-imputation variance (B, the sample variance of the m estimates) inflated by a factor of 1 + 1/m: T = \bar{U} + (1 + 1/m) B;

• the pooled degrees of freedom use the Barnard & Rubin (1999) adjustment, which interpolates between the classical Rubin (1987) degrees of freedom (based only on m and the fraction of missing information) and the complete-data degrees of freedom reported by each nlme::lme() fit. This adjustment matters most at the small sample sizes typical of this package's use case, where the classical, unadjusted degrees of freedom can otherwise come out larger than the complete-data degrees of freedom, which is not sensible.

A large fraction of missing information (fmi, close to 1) for a coefficient means the missing data cost a lot of precision for that coefficient specifically, even if the biomarker's overall missingness rate is modest.

Value

An object of class "poolLongitudinalSub.BJM", a named list with elements:

pooled

A list, one data.frame per longitudinal outcome, with one row per fixed-effect coefficient and columns estimate, std_error, df, statistic, p_value, riv (relative increase in variance due to missingness), and fmi (fraction of missing information).

m

The number of completed-data fits pooled.

long_sub_fixed

The long_sub_fixed used by every fit (taken from the first element of long_fit_all_list).

long_fit_all

A longitudinalSub.BJM object usable for prediction (e.g. as long_fit_all in predictRisk): the first fit, with each biomarker's fixed effects replaced by their pooled estimates, its residual variance by the average across completions, and Sigma_fit by the average of the m fits' Sigma_fit (Rubin's rules point estimates). Averaging the m fits' predictions instead is a common alternative.

References

Rubin, D. B. (1987). Multiple Imputation for Nonresponse in Surveys. John Wiley & Sons.

Barnard, J. and Rubin, D. B. (1999). Small-Sample Degrees of Freedom with Multiple Imputation. Biometrika, 86(4):948-955.

Examples

if (requireNamespace("torch", quietly = TRUE) && torch::torch_is_installed()) {
  data(pbc3)
  data_fit_all <- pbc3[pbc3$status3 == 1, ]

  set.seed(1)
  n <- nrow(data_fit_all)
  data_fit_all$serBilir[sample.int(n, floor(0.1 * n))] <- NA
  data_fit_all$albumin[sample.int(n, floor(0.1 * n))] <- NA

  long_sub_fixed <- list(
    "long1" = serBilir ~ year + age + sex + years,
    "long2" = albumin ~ year + age + sex + years)
  long_sub_random <- list("long1" = ~ year | id, "long2" = ~ year | id)

  imputed <- imputeLongitudinal(data_fit_all, long_sub_fixed, long_sub_random,
                                 time_variable = "year", n_imputations = 5,
                                 impute = "multiple", epochs = 50, seed = 1)

  long_fit_all_list <- lapply(imputed$data_fit_all_list, function(d) {
    longitudinalSub(d, long_sub_fixed, long_sub_random)
  })
  pooled <- poolLongitudinalSub(long_fit_all_list)
  pooled
}

Average the event-type model over the imputations

Description

With competing risks fitIntervalBJM() refits the event-type glm on every imputation, since it uses the event time. For prediction the coefficients are averaged over the imputations (the Rubin's-rules point estimate); the rest of the returned glm object is the first imputation's fit.

Usage

pool_event_type_fit(survival_fit_all_list)

Arguments

Value

The first fit, with the averaged event-type coefficients.

A prediction-ready fit from Rubin's-rules point estimates

Description

Helper for poolLongitudinalSub(): see its long_fit_all return element.

Usage

pooled_longitudinal_fit(long_fit_all_list, pooled)

Arguments

Value

A longitudinalSub.BJM object.

Dynamic prediction function for future longitudinal outcomes

Description

Companion to predictRisk, using the same fitted longitudinal (longitudinalSub) and survival (survivalSub) sub-models, but instead of an event-risk probability this returns a predictive density for a future value of one or more chosen biomarkers (bio_i) at prediction_time + horizon, conditional on the subject's observed longitudinal history up to prediction_time and on survival (being event-free) through prediction_time + horizon, the time the biomarker is predicted at.

This is the single, unified entry point for biomarker-value prediction:

• If bio_i names exactly one biomarker, this predicts just that biomarker and returns a single "predictLongitudinal.BJM" object.

• If bio_i names more than one biomarker, or is left at its default NULL (meaning every biomarker in long_fit_all), this instead predicts every requested biomarker and returns a named list of "predictLongitudinal.BJM" objects (classed "predictLongitudinalAll.BJM"), computing the bio_i-independent part of the pipeline (restricting to at-risk patients, the survival-side integration grid, and the denominator conditional density) only once and reusing it across every biomarker, instead of recomputing it once per biomarker. This matters most under the Gaussian-copula path (see longitudinalSubCopula), where that denominator involves mvtnorm::pmvnorm() Monte-Carlo evaluations that are otherwise the dominant cost of a multi-biomarker prediction run.

See Value for the exact return shape in each case.

bio_i may refer to either continuous or ordinal biomarkers (see longitudinalSubCopula), and a single call may mix both. For an ordinal biomarker, Y_all/Y_predict hold integer category codes (1:K, in the fitted factor's levels()/threshold order) rather than a numeric grid, with Y_all additionally carrying a "category_labels" attribute giving the matching level-label strings; bandcount3 is ignored for that biomarker, since its candidate grid is fixed at its category count (see Details).

The prediction is conditional on the longitudinal history observed up to prediction_time: rows of data_predict_all whose time_variable is later than prediction_time are dropped, with a warning, before predicting.

bandcount2 (controlling the shared survival-integration grid) and bandcount3 (controlling each biomarker's own candidate-value grid) are auto-tuned separately when left at their default "auto", via the doubling-until-stable check described in Details. When bio_i names more than one biomarker, bandcount2 is resolved once using a single representative biomarker (the first one in bio_i), and bandcount3 is then resolved independently for every requested biomarker (their candidate-value grids need not converge at the same resolution) – except for any ordinal biomarker, for which bandcount3 tuning is always skipped (see above) and "bandcount3" is recorded as NA for that biomarker.

Usage

predictLongitudinal(
  bio_i = NULL,
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  horizon,
  time_variable,
  survival_variable_all,
  survival_trans_function,
  bandcount2 = "auto",
  bandcount3 = "auto",
  n_cores = 1
)

Arguments

Details

There is no universal correct value for bandcount2/bandcount3: as a practical check, double both and confirm the results barely change; if they do, keep doubling. By default (bandcount2 = "auto", bandcount3 = "auto"), this doubling check is done for you: starting from small built-in values, both are doubled together, and the result is compared to the previous round, until the largest relative change in Y_predict drops below 1%, or 2 doublings have been tried (so at most 3 calls' worth of work). If it still has not converged by then, a warning reports this and the result at the largest value tried is returned anyway (not an error), so this never silently loops for an unbounded amount of time. Pass an explicit number for either argument to skip auto-tuning it and use a fixed value instead, or call checkBandcountConvergence() directly for more control over the tolerance and doubling count. See also vignette("BJM-intro", package = "BJM") for a worked example.

The integration grid, and the candidate-value grid of a continuous biomarker (which spans the observed values of every patient predicted), are shared by all patients predicted in one call, so a patient's value can differ slightly, within the grids' discretization error, depending on which other patients are predicted alongside it.

The denominator integrates over every event time after prediction_time, including times beyond the last follow-up time in the data survivalSub() was fit on. There, the baseline cumulative hazard is extrapolated linearly (a constant hazard), and the longitudinal sub-model's mean is evaluated at event times it was never fit on. When much of an at-risk patient's survival probability lies beyond the last follow-up time (e.g. a prediction late in follow-up), the prediction depends on this extrapolation; in pbc3 it moved the risks checked by less than one percentage point.

For a fit with an ordinal biomarker, a patient's ordinal measurements contribute the probability of their observed categories. With a single ordinal random effect (e.g. one ordinal biomarker with a random intercept) it is computed by Gauss-Hermite quadrature, and is exact; with more, by Monte Carlo (mvtnorm::pmvnorm()), with an error around the fourth significant digit. Each patient's Monte Carlo draws use a fixed seed derived from the patient's id, so the results are the same on every run and for any n_cores, without set.seed(), and the caller's random-number state is left unchanged.

Value

If bio_i names exactly one biomarker: an object of class "predictLongitudinal.BJM", a named list with elements:

Y_predict

A vector, one entry per at-risk patient (named by patient id), giving the MAP (most likely) predicted value of the biomarker at prediction_time + horizon. For an ordinal biomarker, this is an integer category code (see Details), not a raw value.

Y_density

A probability matrix whose rows correspond to the candidate biomarker values in Y_all and whose columns correspond to individual patients; each entry is the dynamically predicted density of the biomarker taking that value. For an ordinal biomarker, each row is instead that category's predicted probability.

Y_all

The grid of candidate biomarker values used to build Y_density. For an ordinal biomarker, this is the vector of integer category codes 1:K, carrying a "category_labels" attribute with the matching level-label strings (see Details).

Otherwise (bio_i names more than one biomarker, or is NULL): a named list of such objects, one per requested biomarker, named by that biomarker's response-variable name; with attributes "bandcount2" (the single resolved/used bandcount2) and "bandcount3" (a named numeric vector of the resolved/used bandcount3 for each biomarker, or NA for an ordinal biomarker). Classed "predictLongitudinalAll.BJM".

Examples

data(pbc3)

data_survival_fitting =  pbc3[!duplicated(pbc3$id), ]

form_marginal_surv = Surv(years, status3) ~ age + sex
form_conditional_cr = NULL

survival_fit_all = survivalSub(data_survival_fitting, form_marginal_surv,
                               form_conditional_cr)

long_sub_fixed = list(
  "long1" = serBilir ~ year + age + sex +  (years) + (years) * year,
  "long2" = prothrombin ~ year + age + sex + (years) + (years) * year,
  "long3" = albumin ~ year + age + age * year + sex + (years) + (years) * year,
  "long4" = alkaline ~ year + age + sex + (years) + (years) * year,
  "long5" = SGOT ~ year + age + sex + (years) + (years) * year,
  "long6" = platelets ~ year + age + sex + (years)  + (years) * year)

long_sub_random =list(
  "long1" =  ~ year| id,
  "long2" =  ~ year| id,
  "long3" =  ~ year| id,
  "long4" =  ~ year| id,
  "long5" =  ~ year| id,
  "long6" =  ~ year| id)

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)}
)

# Complete case analysis
data_fit_all = list()
for(i in seq_len(length(long_sub_fixed))){
  data_fit_all[[i]] = pbc3[pbc3$status3 == 1, ]
}

# fitting longitudinal submodel
long_fit_all = longitudinalSub(data_fit_all, long_sub_fixed, long_sub_random)

i_PID = 2
data.raw.predict.1 = pbc3[pbc3$id == i_PID, ]

data_predict_all = list()
for(i in seq_len(length(long_sub_fixed))){
  data_predict_all[[i]] = data.raw.predict.1[data.raw.predict.1$year <= 3,]
}

# A single biomarker -> a single predictLongitudinal.BJM result
Y_predict = predictLongitudinal(bio_i = 1, data_predict_all, long_fit_all,
                                survival_fit_all, prediction_time = 3,
                                horizon = 3, time_variable = "year",
                                survival_variable_all, survival_trans_function,
                                bandcount2 = 40, bandcount3 = 400)

# Multiple biomarkers (or bio_i = NULL for all of them) -> a named list
Y_predict_all = predictLongitudinal(bio_i = NULL, data_predict_all, long_fit_all,
                                    survival_fit_all, prediction_time = 3,
                                    horizon = 3, time_variable = "year",
                                    survival_variable_all, survival_trans_function,
                                    bandcount2 = 40, bandcount3 = 400)

Plot risk and future biomarker predictions across a horizon sweep

Description

Fixes the landmark time at prediction_time and sweeps forward over a range of prediction horizons, calling predictRisk (and, when bio_pred is supplied, predictLongitudinal) at each horizon to trace out how predicted event risk (and, optionally, the predicted density/point-forecast of biomarker bio_pred) evolves the further out the prediction window extends from a single fixed landmark. This contrasts with riskPlot, which instead fixes the horizon and sweeps over a set of landmark times.

Usage

predictPlot(
  data_predict_all_one,
  long_fit_all,
  survival_fit_all,
  prediction_time = 4,
  horizon = seq(0, 3, 0.5),
  time_variable,
  survival_variable_all,
  survival_trans_function,
  bandcount1 = "auto",
  bandcount2 = "auto",
  bandcount3 = "auto",
  bio_his = 1,
  bio_pred = 1,
  density = 1,
  n_cores = 1
)

Arguments

Value

Plot of risk and future biomarker with density using dynamic prediction.

Examples

data(pbc3)

data_survival_fitting =  pbc3[!duplicated(pbc3$id), ]

form_marginal_surv = Surv(years, status3) ~ age + sex
form_conditional_cr = NULL

survival_fit_all = survivalSub(data_survival_fitting, form_marginal_surv, 
                               form_conditional_cr)

long_sub_fixed = list(
  "long1" = serBilir ~ year + age + sex +  (years) + (years) * year,  
  "long2" = prothrombin ~ year + age + sex + (years) + (years) * year,  
  "long3" = albumin ~ year + age + age * year + sex + (years) + (years) * year,  
  "long4" = alkaline ~ year + age + sex + (years) + (years) * year, 
  "long5" = SGOT ~ year + age + sex + (years) + (years) * year, 
  "long6" = platelets ~ year + age + sex + (years)  + (years) * year)

long_sub_random =list(
  "long1" =  ~ year| id,   
  "long2" =  ~ year| id,    
  "long3" =  ~ year| id,    
  "long4" =  ~ year| id,    
  "long5" =  ~ year| id,    
  "long6" =  ~ year| id)

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)}
)

# Complete case analysis
data_fit_all = list()
for(i in seq_len(length(long_sub_fixed))){
  data_fit_all[[i]] = pbc3[pbc3$status3 == 1, ]
}

# fitting longitudinal submodel
long_fit_all = longitudinalSub(data_fit_all, long_sub_fixed, long_sub_random)

i_PID = 2
data.raw.predict.plot = pbc3[pbc3$id == i_PID, ]
data_predict_all_one = list(data.raw.predict.plot, data.raw.predict.plot, data.raw.predict.plot,
                            data.raw.predict.plot, data.raw.predict.plot, data.raw.predict.plot)

# plot biomarker 1 history,  predict future biomarker

predictPlot(data_predict_all_one, long_fit_all, survival_fit_all,
            prediction_time = 5, bio_his = 1, bio_pred = 1,
            horizon = seq(0.5, 3.0, 0.5), time_variable = "year",
            survival_variable_all, survival_trans_function,
           bandcount1 = 10, bandcount2 = 10, bandcount3 = 200)
       

    

Predict the number of recurrent events in a window, by rolling episodes forward

Description

Experimental. For a patient at the start of (or partway through) an episode of a gap-time recurrent-event model fit with buildRecurrentEpisodes, predicts the distribution of the number of recurrences in (prediction_time, prediction_time + horizon], together with the probability of dying in it, by Monte Carlo: every patient is simulated forward n_sim times, episode after episode, until the window closes or the patient dies. predictRisk() answers "is the next event a recurrence or a terminal event, and when"; this answers "how many recurrences".

Usage

predictRecurrentCount(
  data_predict,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  horizon,
  time_variable,
  survival_variable_all = list(),
  survival_trans_function = list(),
  episode_variable = "episode",
  next_covariates = NULL,
  stay_lengths = 0,
  n_sim = 200,
  max_stages = 10,
  bandcount2 = 100,
  seed = NULL
)

Arguments

Details

Each simulated path goes as follows.

1. The current (open) episode's event time and type are drawn with simulateTrajectory, from their posterior given the patient's biomarker history in that episode – the same distribution predictRisk integrates, so the first step uses the biomarkers exactly as predictRisk() does.

2. If the event is a terminal event inside the window, the path ends (the patient dies). If it is a recurrence inside the window, the count goes up by one; a stay length is then drawn (from stay_lengths) and the next episode starts at discharge. If the event (or the discharge) falls after the window closes, the path ends.

3. Each later episode is drawn the same way, but with no biomarker history – the patient has not been observed in it yet – so from the model's prior given the patient's covariates, with the episode counter advanced by next_covariates. Repeat until no path is still running, or max_stages episodes have been drawn.

Limitation. Only the first step uses the patient's biomarkers. Later episodes depend on the patient only through their covariates (and the episode counter), so any patient-level tendency to recur that is not in those covariates – a frailty – is not carried forward, and the predicted spread of the counts (e.g. the chance of two or more recurrences) will be too small when such a tendency is present. Putting the number of previous episodes, or the length of the previous gap, among the survival sub-model's covariates (and updating them in next_covariates) recovers some of it.

The stay length is not part of the backward joint model, so it has to be supplied. The default, 0, means a recurrence takes no time (a tumour recurrence, say); for hospitalizations pass observed stay lengths (they are resampled with replacement).

Value

An object of class "predictRecurrentCount.BJM", a list with

summary

One row per patient: mean_count (expected number of recurrences in the window), p_count_0, p_count_1, p_count_2, p_count_3plus, p_ge1, p_ge2, and p_death (probability of dying in the window, outside hospital: a death during a stay is not part of the model, see buildRecurrentEpisodes).

counts

The raw draws: a patients x n_sim matrix of recurrence counts.

died

A patients x n_sim logical matrix: died in the window.

horizon, prediction_time, n_sim

As supplied.

stages

Most episodes drawn on any path.

capped

TRUE if some path was still running when max_stages stopped it (a warning, giving their share, is also given): its later recurrences are not counted, so the counts are slightly low.

capped_fraction

The share of all simulated paths stopped by max_stages (0 if none).

See Also

buildRecurrentEpisodes, predictRisk, simulateTrajectory

Examples

## Simulated recurrent events and a biomarker. Subject 9999 is still at
## risk in their 2nd episode, and is who we predict for.
set.seed(1)
n <- 60
data_events <- do.call(rbind, lapply(seq_len(n), function(id) {
  k <- sample(1:3, 1)
  data.frame(
    id = id, time = cumsum(rexp(k, rate = 1)),
    status = c(rep(1, k - 1), sample(c(0, 2), 1, prob = c(0.6, 0.4))),
    age = round(50 + rnorm(1, sd = 8)), sex = id %% 2
  )
}))
data_events <- rbind(data_events, data.frame(
  id = 9999, time = c(0.8, 1.6), status = c(1, 0), age = 55, sex = 1
))
data_long <- do.call(rbind, lapply(c(seq_len(n), 9999), function(id) {
  tt <- seq(0, max(data_events$time[data_events$id == id]) + 0.2, by = 0.25)
  data.frame(id = id, year = tt, biomarker = rnorm(length(tt)))
}))
built <- buildRecurrentEpisodes(data_events, data_long, time_variable_long = "year")

survival_fit_recur <- survivalSub(
  built$episode_data,
  Surv(gap, event_any) ~ age + sex + episode,
  event_type ~ gap + age + sex + episode
)
long_fit_recur <- longitudinalSub(
  built$long_data[built$long_data$event_any == 1, ],
  biomarker ~ time_in_episode + age + sex + gap + gap * time_in_episode,
  ~ time_in_episode | episode_id
)

## The number of recurrences over the next 3 time units
open <- built$long_data[built$long_data$id == 9999 & built$long_data$event_any == 0, ]
counts <- predictRecurrentCount(
  open, long_fit_recur, survival_fit_recur,
  prediction_time = max(open$time_in_episode), horizon = 3,
  time_variable = "time_in_episode", n_sim = 50, bandcount2 = 30, seed = 1
)
counts

Dynamic prediction function for future event risk

Description

Combines a fitted longitudinal sub-model (longitudinalSub) and survival sub-model (survivalSub) into a backward joint model prediction: for each at-risk subject at prediction_time, the subject's observed longitudinal history up to that time is used to update their individual random effects (empirical Bayes), which are then integrated against the survival sub-model's hazard to give the predicted probability of experiencing the event within (prediction_time, prediction_time + horizon], conditional on being event-free at prediction_time. The numerator (event within the window) is evaluated on a grid over the window whose resolution is controlled by bandcount1; the denominator (survival to prediction_time, used for normalization) is that same window integral plus the integral beyond the window, evaluated on a grid controlled by bandcount2; see Details.

The prediction is conditional on the longitudinal history observed up to prediction_time: rows of data_predict_all whose time_variable is later than prediction_time are dropped, with a warning, before predicting.

A patient with no biomarker measurement up to prediction_time (rows whose biomarker values are all NA, or only rows after prediction_time, e.g. a prediction_time before the first visit) has no history to condition on. The longitudinal factor f(Y \mid T) of the model is then 1, and the patient's risk is that of the survival sub-model alone, given the covariates in their earliest row; a warning names these patients. The other patients' predictions do not change. A patient with measurements of some biomarkers but none of another is still not predicted (NA, with a warning).

Usage

predictRisk(
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  horizon,
  time_variable,
  survival_variable_all,
  survival_trans_function,
  bandcount1 = "auto",
  bandcount2 = "auto",
  n_cores = 1
)

Arguments

Details

There is no universal correct value for bandcount1/bandcount2: as a practical check, double both and confirm the resulting risk probabilities barely change; if they do, keep doubling. By default (bandcount1 = "auto", bandcount2 = "auto"), this doubling check is done for you: starting from small built-in values, both are doubled together, and the result is compared to the previous round, until the largest relative change in the risk probabilities drops below 1%, or 2 doublings have been tried (so at most 3 calls' worth of work). If it still has not converged by then, a warning reports this and the result at the largest value tried is returned anyway (not an error), so this never silently loops for an unbounded amount of time. Pass an explicit number for either argument to skip auto-tuning it and use a fixed value instead (as in previous package versions), or call checkBandcountConvergence() directly for more control over the tolerance and doubling count. See also vignette("BJM-intro", package = "BJM") for a worked example.

The denominator integrates over every event time after prediction_time, including times beyond the last follow-up time in the data survivalSub() was fit on. There, the baseline cumulative hazard is extrapolated linearly (a constant hazard), and the longitudinal sub-model's mean is evaluated at event times it was never fit on. When much of an at-risk patient's survival probability lies beyond the last follow-up time (e.g. a prediction late in follow-up), the prediction depends on this extrapolation; in pbc3 it moved the risks checked by less than one percentage point.

For a fit with an ordinal biomarker, a patient's ordinal measurements contribute the probability of their observed categories. With a single ordinal random effect (e.g. one ordinal biomarker with a random intercept) it is computed by Gauss-Hermite quadrature, and is exact; with more, by Monte Carlo (mvtnorm::pmvnorm()), with an error around the fourth significant digit. Each patient's Monte Carlo draws use a fixed seed derived from the patient's id, so the results are the same on every run and for any n_cores, without set.seed(), and the caller's random-number state is left unchanged.

Value

An object of class "predictRisk.BJM", a named list with elements:

risk_prob_1

A vector of dynamically predicted probabilities, one per at-risk patient, of experiencing the (first) event within the prediction horizon, named by patient id. All 0 when horizon = 0.

risk_prob_2

When survival_fit_all was fit with competing risks, a vector of dynamically predicted probabilities, one per at-risk patient, of experiencing the competing event within the prediction horizon, named by patient id (all 0 when horizon = 0). NULL when there is no competing risk.

prob_undetected_1, prob_undetected_2

Only for an interval-censored survival_fit_all (see below): the probability that the (first / competing) event already happened between the patient's last visit and prediction_time, undetected. prob_undetected_2 is NULL without competing risks.

Interval-censored fits. An event is then only detected at a visit, so at prediction_time s a patient is known to be event-free only up to their last visit V (the last time in their history up to s). risk_prob_1 is then P(s < T \le s + horizon \mid T > V, history) and prob_undetected_1 is P(V < T \le s \mid T > V, history) (with competing risks, each _1/_2 value is that probability jointly with the event type); with V = s they reduce to the usual conditioning on T > s. Only patients still at risk at prediction_time are predicted: a patient whose recorded survival time is before prediction_time is left out (one with a missing survival time is kept). The names show which patients each value belongs to. The numerical integration grids are shared by all patients predicted in one call, so a patient's value can differ slightly (within the grid's discretization error, which shrinks as bandcount1/bandcount2 grow) depending on which other patients are predicted alongside it.

Examples

data(pbc3)

data_survival_fitting =  pbc3[!duplicated(pbc3$id), ]

form_marginal_surv = Surv(years, status3) ~ age + sex
form_conditional_cr = NULL

survival_fit_all = survivalSub(data_survival_fitting, form_marginal_surv, 
                               form_conditional_cr)

long_sub_fixed = list(
  "long1" = serBilir ~ year + age + sex +  (years) + (years) * year,  
  "long2" = prothrombin ~ year + age + sex + (years) + (years) * year,  
  "long3" = albumin ~ year + age + age * year + sex + (years) + (years) * year,  
  "long4" = alkaline ~ year + age + sex + (years) + (years) * year, 
  "long5" = SGOT ~ year + age + sex + (years) + (years) * year, 
  "long6" = platelets ~ year + age + sex + (years)  + (years) * year)

long_sub_random =list(
  "long1" =  ~ year| id,   
  "long2" =  ~ year| id,    
  "long3" =  ~ year| id,    
  "long4" =  ~ year| id,    
  "long5" =  ~ year| id,    
  "long6" =  ~ year| id)

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)}
)

# Complete case analysis
data_fit_all = list()
for(i in seq_len(length(long_sub_fixed))){
  data_fit_all[[i]] = pbc3[pbc3$status3 == 1, ]
}

# fitting longitudinal submodel
long_fit_all = longitudinalSub(data_fit_all, long_sub_fixed, long_sub_random)

i_PID = 2
data.raw.predict.1 = pbc3[pbc3$id == i_PID, ]

data_predict_all = list()
for(i in seq_len(length(long_sub_fixed))){
  data_predict_all[[i]] = data.raw.predict.1[data.raw.predict.1$year <= 3,]
}

# predict risk probability
risk.prob = predictRisk(data_predict_all, long_fit_all, survival_fit_all,
                              prediction_time = 3,
                              horizon = 3, time_variable = "year",
                              survival_variable_all, survival_trans_function,
                              bandcount1 = 10, bandcount2 = 10)

# poly() in its default orthogonal mode, splines::ns()/bs(), and factor()
# in a continuous biomarker's long_sub_fixed produce correct dynamic
# predictions (see ?longitudinalSub) -- including
# when a patient has only a single longitudinal observation to condition
# on -- because the basis/contrasts fit on the full training data are
# cached (via each biomarker's terms object and long_fit_all$xlevels) and
# reused here at prediction time, instead of being recomputed from that
# patient's small per-prediction-time slice:
long_sub_fixed_nonlinear = list(
  "long1" = serBilir ~ poly(year, 2) + age + sex + years,
  "long2" = albumin ~ year + age + sex + years)
long_sub_random_nonlinear = list("long1" = ~ year | id, "long2" = ~ year | id)
long_fit_nonlinear = longitudinalSub(pbc3[pbc3$status3 == 1, ],
                                     long_sub_fixed_nonlinear,
                                     long_sub_random_nonlinear)

data_predict_normal = data.raw.predict.1[data.raw.predict.1$year <= 3, ]
data_predict_sparse = data.raw.predict.1[1, ]

risk.prob.normal = predictRisk(data_predict_normal, long_fit_nonlinear,
                                     survival_fit_all, prediction_time = 3,
                                     horizon = 3, time_variable = "year",
                                     survival_variable_all, survival_trans_function,
                                     bandcount1 = 10, bandcount2 = 10)
risk.prob.sparse = predictRisk(data_predict_sparse, long_fit_nonlinear,
                                     survival_fit_all, prediction_time = 3,
                                     horizon = 3, time_variable = "year",
                                     survival_variable_all, survival_trans_function,
                                     bandcount1 = 10, bandcount2 = 10)
# both give a sane, non-degenerate risk_prob_1 (not 0, no error)
risk.prob.normal$risk_prob_1
risk.prob.sparse$risk_prob_1

Patient ids, in the order predictions are returned

Description

The density helpers (marginalT(), conditionalYT(), ...) all produce one column per patient, in order of first appearance in data_predict_all[[1]]; this returns those ids, as character, to name the returned predictions with.

Usage

prediction_patient_ids(data_predict_all, long_fit_all)

Value

A character vector.

Standard deviation of each patient's predictive distribution

Description

Helper for max_relative_diff(): from a dynamicPredictionBio()-style result (Y_density, one column per patient, tabulated on Y_all), the standard deviation of each patient's predicted biomarker distribution. 0 if the result has no density.

Usage

predictive_sd(result)

Arguments

Value

A numeric vector, one per patient, or 0.

Build the prediction-to-infinity integration grid and marginal survival

Description

Shared helper for predictRisk and dynamicPredictionBio: builds the numerical-integration grid from prediction_time out to upper_bound, and evaluates the marginal survival function S(T) over it.

The bandcount2 + 1 intervals are spaced by probability, not by time: each holds an equal share of the at-risk patients' (averaged) model-based conditional survival mass beyond prediction_time. The upper bound sits far out in the tail (see integration_upper_bound()), so equally spaced intervals would spend most of the grid where there is almost no mass and converge more slowly as bandcount2 grows. The first interval starts exactly at prediction_time; the previous equally spaced grid started half an interval earlier, counting mass from before the prediction time. Each grid point is its interval's midpoint.

Usage

prepare_infinity_grid(
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  upper_bound,
  bandcount2
)

Value

A list with predict.time.infinity (the grid points), predict.time.infinity.1 (the interval edges), and S_T_all_infinity.

Prediction data for predictRisk(), keeping patients without measurements

Description

Without any biomarker measurement up to prediction_time, the longitudinal factor f(Y \mid T) of the backward joint model is 1, and the predicted risk is that of the survival sub-model alone given the patient's covariates. predictRisk() predicts such patients that way instead of refusing them. This runs the usual drop_after_prediction_time(), align_ordinal_levels() and drop_missing_longitudinal(), then puts back one row per patient left without any measurement: their earliest row, with the biomarker set to NA and time_variable set to their last visit up to prediction_time (min(0, prediction_time) if there is none), which an interval-censored fit takes as the last time they were known to be event-free. The covariates come from that row: the survival sub-model's covariates are baseline values, one row per patient. A warning names these patients. split_no_measurement() separates them again; running this on its own output returns the same data, as "auto" bandcount tuning does.

A patient with measurements of some biomarkers but not of others is still reported as unpredictable, as before.

Usage

prepare_risk_data(
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  time_variable,
  survival_variable_all
)

Value

data_predict_all, filtered, with the added rows in the position of each patient's earliest row.

Print a fitIntervalBJM fit

Description

A short summary of a fitIntervalBJM fit: the imputation settings, the survival sub-model's hazard ratios and the longitudinal sub-model's fixed effects (pooled over the imputations, with the fraction of missing information). The components themselves are in the returned list (e.g. x$pooled, x$trace).

Usage

## S3 method for class 'fitIntervalBJM'
print(x, digits = 4, ...)

Arguments

Value

x, invisibly.

Print method for longitudinalSub.BJM objects

Description

Automatically called when you type long_fit_all at the console.

Usage

## S3 method for class 'longitudinalSub.BJM'
print(x, digits = 4, ...)

Arguments

Value

Invisibly returns x.

Print method for poolLongitudinalSub.BJM objects

Description

Automatically called when you type pooled or print(pooled) at the console. Displays a JMbayes2-style formatted summary of the Rubin's-rules-pooled longitudinal sub-model(s).

Usage

## S3 method for class 'poolLongitudinalSub.BJM'
print(x, digits = 4, ...)

Arguments

Value

Invisibly returns x.

Examples

if (requireNamespace("torch", quietly = TRUE) && torch::torch_is_installed()) {
  data(pbc3)
  data_fit_all <- pbc3[pbc3$status3 == 1, ]

  set.seed(1)
  n <- nrow(data_fit_all)
  data_fit_all$serBilir[sample.int(n, floor(0.1 * n))] <- NA
  data_fit_all$albumin[sample.int(n, floor(0.1 * n))] <- NA

  long_sub_fixed <- list(
    "long1" = serBilir ~ year + age + sex + years,
    "long2" = albumin ~ year + age + sex + years)
  long_sub_random <- list("long1" = ~ year | id, "long2" = ~ year | id)

  imputed <- imputeLongitudinal(data_fit_all, long_sub_fixed, long_sub_random,
                                 time_variable = "year", n_imputations = 5,
                                 impute = "multiple", epochs = 50, seed = 1)

  long_fit_all_list <- lapply(imputed$data_fit_all_list, function(d) {
    longitudinalSub(d, long_sub_fixed, long_sub_random)
  })
  pooled <- poolLongitudinalSub(long_fit_all_list)
  pooled   # triggers print.poolLongitudinalSub.BJM automatically
}

Print method for predictLongitudinal.BJM objects

Description

Automatically called when you type the result of predictLongitudinal() (with a single bio_i) at the console.

Usage

## S3 method for class 'predictLongitudinal.BJM'
print(
  x,
  bio_i = NULL,
  long_fit_all = NULL,
  prediction_time = NULL,
  horizon = NULL,
  subject_ids = NULL,
  digits = 4,
  ...
)

Arguments

Value

Invisibly returns x.

Print method for predictLongitudinalAll.BJM objects

Description

Called when you print the result of predictLongitudinal() with several biomarkers (bio_i naming more than one, or NULL): prints each biomarker's prediction in turn, under its name.

Usage

## S3 method for class 'predictLongitudinalAll.BJM'
print(x, ...)

Arguments

Value

Invisibly returns x.

Print method for predictRisk.BJM objects

Description

Automatically called when you type the result of predictRisk() at the console.

Usage

## S3 method for class 'predictRisk.BJM'
print(
  x,
  prediction_time = NULL,
  horizon = NULL,
  subject_ids = NULL,
  digits = 4,
  ...
)

Arguments

Value

Invisibly returns x.

Print method for survivalSub.BJM objects

Description

Automatically called when you type survival_fit_all or print(survival_fit_all) at the console. Displays a JMbayes2-style formatted summary of the survival sub-model.

Usage

## S3 method for class 'survivalSub.BJM'
print(x, digits = 4, ...)

Arguments

Value

Invisibly returns x.

Examples

data(pbc3)
data_survival_fitting <- pbc3[!duplicated(pbc3$id), ]
form_marginal_surv  <- Surv(years, status3) ~ age + sex
form_conditional_cr <- status4 ~ years + age + sex
survival_fit_all  <- survivalSub(data_survival_fitting,
                                 form_marginal_surv, form_conditional_cr)
survival_fit_all   # triggers print.survivalSub.BJM automatically

Log probability of an ordinal category and its derivatives

Description

\log P(lo < \eta + \epsilon \le hi) for \epsilon \sim N(0, 1), with its first and second derivatives in \eta; evaluated on the side of the mean where the interval lies, so a category far in a tail keeps its precision.

Usage

probit_box_derivs(eta, lo, hi)

Arguments

Value

A list with logP, g (first derivative) and h (second derivative, \le 0).

Factor one patient's joint longitudinal covariance for conditionalYTBio()/conditionalYDTBio()

Description

Factor one patient's joint longitudinal covariance for conditionalYTBio()/conditionalYDTBio()

Usage

process_variance(
  num_i,
  time_new,
  bio_i,
  data_predict_all,
  long_fit_all,
  time_variable
)

Value

The result of cov_factor(), or NA if the patient has no rows for some biomarker.

Random-effects design matrix for one patient, across all biomarkers

Description

Shared helper for build_conditional_design(), build_conditional_design_copula(), and process_variance(): builds each biomarker's random-effects design from the left-hand side of its own long_sub_random formula (e.g. ~ 1, ~ year, ~ year + I(year^2)) and stacks them block-diagonally, matching the block-diagonal Sigma_fit built by longitudinalSub(). Previously this assumed every biomarker had either a random intercept only, or a random intercept plus a slope on time_variable, so a mix of the two (or any other random-effects formula) failed with "non-conformable arguments".

Usage

random_effects_design(data_num_i_list, long_sub_random, Sigma)

Arguments

Value

A matrix with one row per observation (stacked across biomarkers) and one column per random effect.

Resolve the per-biomarker type used by longitudinalSub()

Description

If biomarker_type is supplied explicitly, it always wins (after validation). Otherwise each biomarker's type is auto-detected from its own response column in data_fit_all: a factor (or ordered factor) response is treated as "ordinal"; anything else is treated as "continuous". This mirrors exactly the priority rule requested for longitudinalSub: explicit biomarker_type overrides auto-detection, and auto-detection is only consulted when biomarker_type is NULL.

Usage

resolve_biomarker_type(biomarker_type, data_fit_all, long_sub_fixed, M)

Plot predicted risk across a sweep of landmark times

Description

Fixes the prediction horizon and sweeps backward/forward over a set of landmark times, calling predictRisk at each landmark in prediction_time to trace out how predicted event risk over that fixed-length window changes as the landmark moves – i.e. as more (or less) of the subject's longitudinal history is used to update their risk. prediction_time may be a vector of landmark times to use directly; a single value, which is expanded to three landmarks (prediction_time, 1.5 * prediction_time, 2 * prediction_time); or NULL, which uses every observed longitudinal measurement time (across subjects, from the first biomarker's data) as a landmark. This contrasts with predictPlot, which instead fixes the landmark and sweeps over a range of horizons. The optional bio_i argument only selects which biomarker's observed trajectory is overlaid on the plot for visual reference – it does not affect the risk computation itself (unlike predictPlot's bio_pred, which drives a biomarker density prediction).

When the data give the subject's survival time and status, a vertical line marks the end of follow-up: solid for an observed event (coloured by event type with competing risks), dashed grey for a censoring.

Usage

riskPlot(
  data_predict_all_pre,
  long_fit_all,
  survival_fit_all,
  prediction_time = NULL,
  bio_i = NULL,
  horizon,
  time_variable,
  survival_variable_all,
  survival_trans_function,
  bandcount1 = "auto",
  bandcount2 = "auto",
  n_cores = 1
)

Arguments

Value

Plot of risk using dynamic prediction.

Examples

data(pbc3)
survival_fit_all <- survivalSub(pbc3[!duplicated(pbc3$id), ],
                                Surv(years, status3) ~ age + sex, NULL)
long_fit_all <- longitudinalSub(pbc3[pbc3$status3 == 1, ],
                                list(serBilir ~ year + age + sex + years,
                                     albumin ~ year + age + sex + years),
                                list(~ year | id, ~ year | id))

# one patient's whole history; risks within 2 years of landmarks 2, 3, 4
patient <- pbc3[pbc3$id == 2, ]
riskPlot(patient, long_fit_all, survival_fit_all, prediction_time = c(2, 3, 4),
         bio_i = 1, horizon = 2, time_variable = "year",
         survival_variable_all = list(), survival_trans_function = list(),
         bandcount1 = 10, bandcount2 = 20)

Biomarker axis with the risk probabilities on a secondary axis

Description

Helper for riskPlot() and predictPlot(), which draw risks on the biomarker's axis multiplied by scale_prob. The secondary (risk) axis is labelled only within [0, 1]: where the biomarker goes below 0 it used to show negative "probabilities". For an ordinal biomarker the primary axis is labelled with its categories.

Usage

risk_axis_scale(scale_prob, category_labels = NULL)

Arguments

Value

A ggplot2 scale.

Restricted cubic spline basis (Royston–Parmar)

Description

Columns x and, for each internal knot k_j, (x - k_j)_+^3 - \lambda_j (x - k_{min})_+^3 - (1 - \lambda_j) (x - k_{max})_+^3, \lambda_j = (k_{max} - k_j) / (k_{max} - k_{min}): a cubic spline that is linear beyond the boundary knots. With deriv = TRUE, the derivative of each column in x.

Usage

rp_basis(x, knots, deriv = FALSE)

Arguments

Draw from a truncated multivariate normal distribution

Description

n draws of N(\mu, V) restricted to the box (lower, upper]: by rejection from batches of untruncated proposals, and, for any draws still missing once max_proposals proposals have been spent (a box of small probability), by Gibbs sampling, one coordinate at a time from its univariate truncated normal conditional distribution, after burn_in sweeps and keeping every thin-th one.

Usage

rtmvnorm_box(
  n,
  mean,
  sigma,
  lower,
  upper,
  max_proposals = 20000,
  burn_in = 50,
  thin = 5
)

Value

A matrix with one column per draw.

Draw from a univariate truncated normal distribution

Description

By inversion, on the side of the mean where the truncation interval lies (so a box far in the upper tail does not lose precision to pnorm() rounding to 1).

Usage

rtnorm_one(mean, sd, lower, upper)

Vectorized draws from truncated standard-variance normals

Description

N(mean, 1) truncated to (lower, upper], by inversion on the side of the mean where the interval lies.

Usage

rtnorm_vec(mean, lower, upper)

Arguments

Value

A numeric vector of draws.

Combine one coefficient's per-completion estimates with Rubin's rules

Description

Internal helper for poolLongitudinalSub: implements the scalar-parameter pooling rules of Rubin (1987, Ch. 3) and the Barnard & Rubin (1999) small-sample degrees-of-freedom adjustment, for a single fixed-effect coefficient's m per-completion estimates/standard errors.

Usage

rubin_pool_scalar(estimates, std_errors, dfcom)

Arguments

Filter longitudinal data down to a single patient

Description

Shared helper for conditionalYT and conditionalYDT: extracts each biomarker's rows for one patient ID, replicating the shared data frame across biomarkers when only one was supplied.

Usage

select_patient_longitudinal_data(
  data.long,
  num,
  num_i,
  n_longitudinal,
  time_variable
)

Value

A list with rep_num_i_list and data_num_i_list, or NULL if any biomarker has zero rows for this patient (caller should skip the patient in that case).

Filter longitudinal data down to a single patient, substituting the candidate biomarker prediction value

Description

Shared helper for conditionalYTBio and conditionalYDTBio: extracts each biomarker's rows for one patient ID, and for the biomarker being predicted (bio_i), appends a row at time_new with each candidate value in Y_all substituted in turn.

Usage

select_patient_longitudinal_data_bio(
  data.long,
  num,
  num_i,
  n_longitudinal,
  time_variable,
  bio_i,
  time_new,
  Y_all,
  long_fit_all
)

Value

A list with rep_num_i_list, data_num_i_list, and Y_select_all (a matrix of the substituted biomarker values, one column per element of Y_all), or NULL if any biomarker has zero rows for this patient (caller should skip the patient in that case).

Set a patient's survival-time columns to one integration grid point

Description

Sets the survival-time column and each survival_variable_all column (through its survival_trans_function) to their values at event time l.

Usage

set_survival_columns(
  df,
  survival_variable,
  l,
  survival_variable_all,
  survival_trans_function
)

Value

df, updated.

Simulate future trajectories from the backward joint model

Description

Draws n_sim complete futures for each patient from the fitted backward joint model, conditional on the patient's longitudinal history up to prediction_time and on being event-free at prediction_time: an event time (and, with competing risks, an event type), and the values of every biomarker at the requested times. Where predictRisk and predictLongitudinal summarize the predictive distribution, this returns samples from it, so it can be used to generate synthetic patients or to propagate the prediction into any other quantity (e.g. a patient's "digital twin" under the fitted model).

Each draw follows the model's factorization f(T, D)\, f(Y \mid T, D):

1. (T, D) is drawn from its posterior given the history, proportional to f(Y_{hist} \mid T, D)\, f(D \mid T)\, f(T) restricted to T > prediction_time, evaluated on the same integration grid as predictRisk's denominator;

2. the random effects b are drawn from their Gaussian posterior given the history and the drawn (T, D);

3. the biomarkers at times are drawn from the longitudinal sub-model given (T, D, b), adding residual error.

A patient with no non-missing biomarker measurement up to prediction_time (e.g. a single baseline row with the biomarkers set to NA) is simulated from the model's prior given their baseline covariates, i.e. as a new synthetic patient.

Usage

simulateTrajectory(
  data_predict_all,
  long_fit_all,
  survival_fit_all,
  prediction_time,
  times,
  time_variable,
  survival_variable_all,
  survival_trans_function,
  n_sim = 100,
  bandcount2 = 100,
  max_event_time = NULL,
  truncate = TRUE,
  future_covariates = NULL,
  seed = NULL
)

Arguments