gmvjoint (version 0.4.0)

dynPred: Dynamic predictions for survival sub-model in a multivariate joint model.

Description

Calculates individualised conditional survival probabilities for subjects during a period of follow-up using a joint model fit along with requisite longitudinal process history.

Note that this function is largely designed for use within the ROC function which assesses discriminatory power of the joint model, however it does function by itself with proper use of its arguments.

Usage

dynPred(
  data,
  id,
  fit,
  u = NULL,
  nsim = 200,
  progress = TRUE,
  scale = NULL,
  df = NULL
)

Value

A list of class dynPred which consists of three items:

pi

A data.frame which contains each candidate failure time (supplied by u), with the mean, median and 2.5% and 97.5% quantiles of probability of survival until this failure time.

pi.raw

A matrix of with nsim rows and length(u) columns, each row represents the \(l\)th conditional survival probability of survival each u survival time. This is largely for debugging purposes.

MH.accept

The acceptance rate of the Metropolis-Hastings algorithm on the random effects.

Arguments

data

the data to which the original joint model was fit.

id

subject identifier, i.e. for which subject is the conditional survival probabilities desired?

fit

a joint model fit by the joint function.

u

a numeric vector of candidate follow-up times for which a survival probability should be calculated. Note that the first item u[1] denotes the start of the "window" and is dropped from calculations. If u=NULL (the default), then the probability of surviving all failure times after the id's final longitudinal time is calculated.

nsim

how many Monte Carlo simulations should be carried out? Defaults to nsim=200. First-order estimates are calculated if nsim=0.

progress

a logical, if progress=TRUE (the default) then a progress bar displaying the current percentage of simulations have been completed.

scale

numeric scales the variance-covariance parameter in the proposal distribution for the Metropolis-Hastings algorithm. Defaults to scale = NULL which doesn't scale the variance term at all. Users are encouraged to experiment with values here; this parameter controls the acceptance rate of the MH scheme.

df

numeric denotes the degrees of freedom of the proposed \(t\) distribution on the random effects; df=4 is suggested and is the default.

Author

James Murray ().

Details

Dynamic predictions for the time-to-event process based on information available on the subject's longitudinal process up to given time \(t\) are calculated by Monte Carlo simulation outlined in Rizopoulos (2011). For a subject last observed at time \(t\), the probability that they survive until future time \(u\) is

$$ Pr(T_i \ge u | T \ge t; \boldsymbol{Y}_i, \boldsymbol{b}_i; \boldsymbol{\Omega}) \approx \frac{S(u|\hat{\boldsymbol{b}}_i; \boldsymbol{\Omega})} {S(t|\hat{\boldsymbol{b}}_i; \boldsymbol{\Omega})} $$ where \(T_i\) is the true failure time for subject \(i\), \(\boldsymbol{Y}_i\) their longitudinal measurements up to time \(t\), and \(S()\) the survival function.

\(\boldsymbol{\Omega}\) is drawn from the multivariate normal distribution with mean \(\hat{\boldsymbol{\Omega}}\) and its variance taken from a fitted joint object. \(\hat{\boldsymbol{b}}\) is drawn from the \(t\) distribution by means of a Metropolis-Hastings algorithm with nsim iterations.

References

Bernhardt PW, Zhang D and Wang HJ. A fast EM Algorithm for Fitting Joint Models of a Binary Response to Multiple Longitudinal Covariates Subject to Detection Limits. Computational Statistics and Data Analysis 2015; 85; 37--53

Rizopoulos D. Dynamic predictions and prospective accuracy in joint models for longitudinal and time-to-event data. Biometrics 2011; 67: 819–829.

See Also

ROC and plot.dynPred.

Examples

# \donttest{
data(PBC)
PBC$serBilir <- log(PBC$serBilir)
# Focus in on id 81, who fails at around 7 years of follow-up. \code{dynPred} allows us to
# infer how the model believes their survival probability would've progressed (ignoring the
# true outcome at start time).
# Univariate -----------------------------------------------------------
long.formulas <- list(serBilir ~ drug * time + (1 + time|id))
surv.formula <- Surv(survtime, status) ~ drug
family <- list('gaussian')
fit <- joint(long.formulas, surv.formula, PBC, family)
preds <- dynPred(PBC, id = 81, fit = fit, u = NULL, nsim = 200,
                 scale = 2)
preds
plot(preds)
# Bivariate ------------------------------------------------------------
# Does introduction of albumin affect conditional survival probability?
long.formulas <- list(
  serBilir ~ drug * time + I(time^2) + (1 + time + I(time^2)|id),
  albumin ~ drug * time + (1 + time|id)
)
fit <- joint(long.formulas, surv.formula, data = PBC, family = list("gaussian", "gaussian"))
bi.preds <- dynPred(PBC, id = 81, fit = fit, u = NULL, nsim = 200, 
                    scale = fit$coeffs$D/sqrt(fit$ModelInfo$n))
bi.preds
plot(bi.preds) # Appears to level-off dramatically; perhaps indicative of this id's albumin
               # levels, or acceleration in serBilir trajectory around 8.5 years.
# }