gmvjoint (version 0.3.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,
  b.density = c("normal", "t"),
  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.

b.density

character string imposing the density for the current and proposal for each draw of the random effects \(\boldsymbol{b}^{(l)}, l=1,\dots,\code{nsim}\). Options are b.density="normal" i.e. using the assumption that the conditional density of the random effects is

$$\boldsymbol{b}_i^{(l)}|\boldsymbol{Y}_i(u);\boldsymbol{\Omega}^{(l)}\sim N(\hat{\boldsymbol{b}}_i^{(l)},\hat{\boldsymbol{\Sigma}}_i^{(l)})$$

(i.e. in line with proposal by Bernhardt et al. (2015)), or in keeping with other literature surrounding dynamic predictions (e.g. Rizopoulos (2011)) in imposing the b.density="t" distribution.

scale

if b.density = "t" then this 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.

df

if b.density = "t" then this numeric denotes the degrees of freedom of the proposed \(t\) distribution on the random effects. df=4 is suggested.

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 either the (multivariate) Normal, or \(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)
# 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, b.density = 'normal',
                 scale = 0.18)
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, b.density = 'normal',
                    scale = 0.66)
bi.preds
plot(bi.preds) # Appears to level-off quicker; perhaps indicative of this id's albumin levels.
# }