Package {BJM}

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 dynamicPrediction(), 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 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)

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 dynamicPrediction(), 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
)

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 dynamicPrediction()/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).

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 dynamicPrediction()/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.

Usage

bandcount_auto_start

Format

An object of class numeric of length 3.

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
)

Value

A list with longitudinal_all_matrix, parameter_matrix, Sigma_all, det_Var_cov_estep, Sigma_all_solve, and long_sigma_long.

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.

Check whether bandcount1/bandcount2/bandcount3 are large enough

Description

dynamicPrediction()/dynamicPredictionBio() approximate the integrals behind the predicted risk probabilities (and, for dynamicPredictionBio(), 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 dynamicPredictionBio()'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(
  dynamicPrediction, 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)

Normalize a risk-probability ratio into a valid probability

Description

Shared helper for dynamicPrediction 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. All patients in this plot experience events at the same time point, specified by condi_time2event. Several evenly spaced time points between the baseline and condi_time2event are selected for plotting. Each point is calculated using the mean value of all patients' biomarker values at that time point. The interval between two time points is defined by interval_time

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

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
)

pbc.cmt


# example with competing risks

data(pbc3)

data_plot_all = pbc3[!is.na(pbc3$status4),]

pbc.cmt.cr <- cmtPlot(data_plot_all, 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

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. 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. 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|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. 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. 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.

Dynamic prediction function

Description

The time values in the prediction data subset must be less than the specified prediction_time which is the prediction time. The time points for longitudinal repeated measurements must not surpass the prediction time.

Usage

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

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.

Value

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

risk_prob_1

A vector of dynamically predicted probabilities, one per patient, of experiencing the (first) event within the prediction horizon. 0 when horizon <= 0.

risk_prob_2

When survival_fit_all was fit with competing risks, a vector of dynamically predicted probabilities, one per patient, of experiencing the competing event within the prediction horizon. NULL when there is no competing risk, or when horizon <= 0.

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 = dynamicPrediction(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)

Dynamic prediction function for future biomarker

Description

The time values in the prediction data subset must be less than the specified prediction_time which is the prediction time. The time points for longitudinal repeated measurements must not surpass the prediction time.

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 "dynamicPredictionBio.BJM", a named list with elements:

Y_predict

A vector, one entry per patient, giving the MAP (most likely) predicted value of biomarker bio_i at prediction_time + horizon.

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.

Y_all

The grid of candidate biomarker values used to build Y_density.

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,]
}

Y_predict = dynamicPredictionBio(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)

The process involves estimating parameters for a multivariate linear mixed-effects model, which simultaneously analyzes multiple dependent variables that may be correlated. This approach incorporates both fixed effects, which are consistent across the population, and random effects, accounting for variations within groups or subjects. By fitting this model, one can assess the influence of predictor variables on several longitudinal outcomes while considering the inherent variability in the data due to random effects.

Description

The process involves estimating parameters for a multivariate linear mixed-effects model, which simultaneously analyzes multiple dependent variables that may be correlated. This approach incorporates both fixed effects, which are consistent across the population, and random effects, accounting for variations within groups or subjects. By fitting this model, one can assess the influence of predictor variables on several longitudinal outcomes while considering the inherent variability in the data due to random effects.

Usage

longitudinalSub(data_fit_all, long_sub_fixed, long_sub_random)

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)

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)

# poly() in its default orthogonal mode, splines::ns()/bs(), and
# factor() trigger a warning (see the long_sub_fixed argument above),
# but are still safe to use: the terms/xlevels/contrasts fit on the
# full training data are cached and reused 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)

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

Description

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

Usage

longitudinalSubVar(thetaLong, l, tol.em, verbose)

Marginal distribution of T

Description

Marginal distribution of T

Usage

marginalT(data_predict_all, long_fit_all, survival_fit_all, l_i, upper_bound)

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.

Compare two prediction results' plain numeric-vector fields

Description

Shared comparison logic for checkBandcountConvergence() and the "auto" bandcount support in dynamicPrediction()/ 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).

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.

Plot of risk and future biomarker with density using dynamic prediction

Description

This function gives the risk and biomarker prediction plot.

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)
       

    

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

Description

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

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, predict.time.infinity.1, and S_T_all_infinity.

Print method for dynamicPrediction.BJM objects

Description

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

Usage

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

Arguments

Value

Invisibly returns x.

Print method for dynamicPredictionBio.BJM objects

Description

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

Usage

## S3 method for class 'dynamicPredictionBio.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 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 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

Print both sub-models of a fitted backward joint model

Description

Convenience function that prints the longitudinal and survival sub-model summaries together. Unlike print.longitudinalSub.BJM and print.survivalSub.BJM, this does not dispatch on a single "BJM"-classed object, because longitudinalSub and survivalSub are fit and returned separately; it simply prints both fit objects you already have.

Usage

printBJM(long_fit_all, survival_fit_all, digits = 4)

Arguments

Value

Invisibly returns a named list with both fit objects.

Construct variance

Description

Construct variance

Usage

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

Plot of risk using dynamic prediction

Description

This function gives the risk prediction plot.

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.

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

Restrict prediction data to patients still at risk

Description

Shared helper for dynamicPrediction and dynamicPredictionBio: drops rows whose survival-time variable is below prediction_time from every biomarker's data frame.

Usage

subset_at_risk(data_predict_all, survival_variable, prediction_time)

Value

data_predict_all, filtered in place per element.

Summary method for dynamicPrediction.BJM objects

Description

Like print but also shows mean, SD, and range of predicted risks.

Usage

## S3 method for class 'dynamicPrediction.BJM'
summary(
  object,
  prediction_time = NULL,
  horizon = NULL,
  subject_ids = NULL,
  digits = 4,
  ...
)

Arguments

Value

Invisibly returns object.

Summary method for dynamicPredictionBio.BJM objects

Description

Like print but also shows distribution-level summaries across subjects.

Usage

## S3 method for class 'dynamicPredictionBio.BJM'
summary(
  object,
  bio_i = NULL,
  long_fit_all = NULL,
  prediction_time = NULL,
  horizon = NULL,
  subject_ids = NULL,
  digits = 4,
  ...
)

Arguments

Value

Invisibly returns object.

Summary method for longitudinalSub.BJM objects

Description

Like print but adds per-outcome random-effects variance components and the full correlation matrix of D.

Usage

## S3 method for class 'longitudinalSub.BJM'
summary(object, digits = 4, ...)

Arguments

Value

Invisibly returns a list of per-outcome summary.lme objects.

Summary method for survivalSub.BJM objects

Description

Called via summary(survival_fit_all). Returns (and prints) an extended summary including baseline hazard range, BIC, and McFadden R2 for the competing-risks GLM.

Usage

## S3 method for class 'survivalSub.BJM'
summary(object, digits = 4, ...)

Arguments

Value

Invisibly returns a list with components cox_summary and (if competing risks) glm_summary.

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)
summary(survival_fit_all)

Fitting survival sub-model

Description

Fitting survival sub-model

Usage

survivalSub(data_survival_fitting, form_marginal_surv, form_conditional_cr)

Arguments

Value

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

coxph_fit

The fitted coxph marginal survival model.

form_marginal_surv

The form_marginal_surv formula, as supplied.

glm_fit

The fitted glm competing-risks (event type) model, or NULL if form_conditional_cr was not supplied.

form_conditional_cr

The form_conditional_cr formula, as supplied (or NULL).

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)

Build a survival-time transform basis from cut points

Description

dynamicPrediction(), dynamicPredictionBio(), predictPlot(), and riskPlot() all take a pair of arguments, survival_variable_all/survival_trans_function, that describe transformed basis variables of the (remaining) survival time; these can optionally be referenced in long_sub_fixed formulas to let the longitudinal sub-model depend flexibly on time-to-event. In every example in this package, that pair follows the same convention: variables named "Tyears1", "Tyears2", ... , each defined as the absolute distance from a fixed cut point (function(x) abs(x - k)). survivalTrans() builds exactly that pair from a plain vector of cut points, so you do not have to hand-write two parallel lists of matching names and closures. You remain free to construct survival_variable_all/survival_trans_function by hand for any other transform.

Usage

survivalTrans(cut_points, prefix = "Tyears")

Arguments

Value

A named list with elements survival_variable_all and survival_trans_function, in the format expected by dynamicPrediction(), dynamicPredictionBio(), predictPlot(), and riskPlot().

Examples

trans <- survivalTrans(c(1, 3, 5, 7))
trans$survival_variable_all
trans$survival_trans_function[[1]](2)

# equivalent to hand-writing:
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)
)

Warn about formula terms whose basis is recomputed from whatever data they are given

Description

poly() (in its default orthogonal mode), splines::ns()/splines::bs(), and factor() compute their basis/contrasts from whatever data is passed to model.matrix(). BJM's prediction functions rebuild the design matrix from a small, patient-specific slice of data at every point on the internal prediction grid, which is not the data the model was fit on, so the basis recomputed at prediction time silently does not match the one used at fitting time (or, with too few distinct values, model.matrix() fails outright). poly(..., raw = TRUE), I(x^2), log(), sqrt(), and similar terms that do not depend on the surrounding data are unaffected and are not flagged.

Usage

warn_unsafe_formula_terms(formula_list, arg_name)