Joint Modeling of Longitudinal and Time-to-Event Data in R

Description

This package fits shared parameter models for the joint modeling of normal longitudinal responses and event times under a maximum likelihood approach. Various options for the survival model and optimization/integration algorithms are provided.

Details

The package has a single model-fitting function called jointModel, which accepts as main arguments a linear mixed effects object fit returned by function lme() of package nlme, and a survival object fit returned by either function coxph() or function survreg() of package survival. In addition, the method argument of jointModel() specifies the type of the survival submodel to be fitted and the type of the numerical integration technique; available options are:

"Cox-PH-GH"

the time-dependent version of a proportional hazards model with unspecified baseline hazard function. The Gauss-Hermite integration rule is used to approximate the required integrals. (This option corresponds to the joint model proposed by Wulfsohn and Tsiatis, 1997)

"weibull-PH-GH"

the Weibull model under the proportional hazards formulation. The Gauss-Hermite integration rule is used to approximate the required integrals.

"weibull-AFT-GH"

the Weibull model under the accelerated failure time formulation. The Gauss-Hermite integration rule is used to approximate the required integrals.

"piecewise-PH-GH"

a proportional hazards model with a piecewise constant baseline risk function. The Gauss-Hermite integration rule is used to approximate the required integrals.

"spline-PH-GH"

a proportional hazards model, in which the log baseline hazard is approximated using B-splines. The Gauss-Hermite integration rule is used to approximate the required integrals.

"ch-Laplace"

an additive log cumulative hazard model, in which the log cumulative baseline hazard is approximated using B-splines. A fully exponential Laplace approximation method is used to approximate the required integrals (Rizopoulos et al., 2009).

For all the above mentioned options (except the last one), a pseudo-adaptive Gauss-Hermite integration rule is also available (Rizopoulos, 2012b). This is much faster than the classical Gauss-Hermite rule, and in several simulations it has been shown to perform equally well (though its performance is still under investigation).

The package also offers several utility functions that can extract useful information from fitted joint models. The most important of those are included in the See also section below.

Author(s)

Dimitris Rizopoulos

Maintainer: Dimitris Rizopoulos <d.rizopoulos@erasmusmc.nl>

References

Henderson, R., Diggle, P. and Dobson, A. (2000) Joint modelling of longitudinal measurements and event time data. Biostatistics 1, 465–480.

Rizopoulos, D. (2012a) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

Rizopoulos, D. (2012b) Fast fitting of joint models for longitudinal and event time data using a pseudo-adaptive Gaussian quadrature rule. Computational Statistics and Data Analysis 56, 491–501.

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

Rizopoulos, D. (2010) JM: An R package for the joint modelling of longitudinal and time-to-event data. Journal of Statistical Software 35 (9), 1–33. doi:10.18637/jss.v035.i09

Rizopoulos, D., Verbeke, G. and Lesaffre, E. (2009) Fully exponential Laplace approximation for the joint modelling of survival and longitudinal data. Journal of the Royal Statistical Society, Series B 71, 637–654.

Rizopoulos, D., Verbeke, G. and Molenberghs, G. (2010) Multiple-imputation-based residuals and diagnostic plots for joint models of longitudinal and survival outcomes. Biometrics 66, 20–29.

Tsiatis, A. and Davidian, M. (2004) Joint modeling of longitudinal and time-to-event data: an overview. Statistica Sinica 14, 809–834.

Wulfsohn, M. and Tsiatis, A. (1997) A joint model for survival and longitudinal data measured with error. Biometrics 53, 330–339.

See Also

jointModel, survfitJM, rocJM, aucJM, dynCJM, prederrJM, predict

Didanosine versus Zalcitabine in HIV Patients

Description

A randomized clinical trial in which both longitudinal and survival data were collected to compare the efficacy and safety of two antiretroviral drugs in treating patients who had failed or were intolerant of zidovudine (AZT) therapy.

Format

A data frame with 1408 observations on the following 9 variables.

patient

patients identifier; in total there are 467 patients.

Time

the time to death or censoring.

death

a numeric vector with 0 denoting censoring and 1 death.

CD4

the CD4 cells count.

obstime

the time points at which the CD4 cells count was recorded.

drug

a factor with levels ddC denoting zalcitabine and ddI denoting didanosine.

gender

a factor with levels female and male.

prevOI

a factor with levels AIDS denoting previous opportunistic infection (AIDS diagnosis) at study entry, and noAIDS denoting no previous infection.

AZT

a factor with levels intolerance and failure denoting AZT intolerance and AZT failure, respectively.

Note

The data frame aids.id contains the first CD4 cell count measurement for each patient. This data frame is used to fit the survival model.

References

Goldman, A., Carlin, B., Crane, L., Launer, C., Korvick, J., Deyton, L. and Abrams, D. (1996) Response of CD4+ and clinical consequences to treatment using ddI or ddC in patients with advanced HIV infection. Journal of Acquired Immune Deficiency Syndromes and Human Retrovirology 11, 161–169.

Guo, X. and Carlin, B. (2004) Separate and joint modeling of longitudinal and event time data using standard computer packages. The American Statistician 58, 16–24.

Examples

summary(aids.id)

Anova Method for Fitted Joint Models

Description

Produces marginal Wald tests or Performs a likelihood ratio test between two nested joint models.

Usage

## S3 method for class 'jointModel'
anova(object, object2, test = TRUE, 
   process = c("both", "Longitudinal", "Event"), L = NULL, ...)

Arguments

Value

An object of class aov.jointModel with components,

Warning

The code minimally checks whether the models are nested! The user is responsible to supply nested models in order the LRT to be valid.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

Rizopoulos, D. (2010) JM: An R Package for the Joint Modelling of Longitudinal and Time-to-Event Data. Journal of Statistical Software 35 (9), 1–33. doi:10.18637/jss.v035.i09

See Also

jointModel

Examples

## Not run: 
# linear mixed model fit without treatment effect
fitLME.null <- lme(sqrt(CD4) ~ obstime, 
    random = ~ 1 | patient, data = aids)
# cox model fit without treatment effect
fitCOX.null <- coxph(Surv(Time, death) ~ 1, 
    data = aids.id, x = TRUE)
# joint model fit without treatment effect
fitJOINT.null <- jointModel(fitLME.null, fitCOX.null, 
    timeVar = "obstime", method = "weibull-PH-aGH")

# linear mixed model fit with treatment effect
fitLME.alt <- lme(sqrt(CD4) ~ obstime * drug - drug, 
    random = ~ 1 | patient, data = aids)
# cox model fit with treatment effect
fitCOX.alt <- coxph(Surv(Time, death) ~ drug, 
    data = aids.id, x = TRUE)
# joint model fit with treatment effect
fitJOINT.alt <- jointModel(fitLME.alt, fitCOX.alt, timeVar = "obstime", 
    method = "weibull-PH-aGH")

# likelihood ratio test for treatment effect
anova(fitJOINT.null, fitJOINT.alt)

## End(Not run)

Time-Dependent AUCs for Joint Models

Description

Using the available longitudinal information up to a starting time point, this function computes an estimate of the prediction error of survival at a horizon time point based on joint models.

Usage

aucJM(object, newdata, Tstart, ...)

## S3 method for class 'jointModel'
aucJM(object, newdata, Tstart, Thoriz = NULL, 
    Dt = NULL, idVar = "id", simulate = FALSE, M = 100, ...)

Arguments

Details

Based on a fitted joint model (represented by object) and using the data supplied in argument newdata, this function computes the following estimate of the AUC:

\mbox{AUC}(t, \Delta t) = \mbox{Pr} \bigl [ \pi_i(t + \Delta t \mid t) < \pi_j(t + \Delta t \mid t) \mid \{ T_i^* \in (t, t + \Delta t] \} \cap \{ T_j^* > t + \Delta t \} \bigr ],

with i and j denote a randomly selected pair of subjects, and \pi_i(t + \Delta t \mid t) and \pi_j(t + \Delta t \mid t) denote the conditional survival probabilities calculated by survfitJM for these two subjects, for different time windows \Delta t specified by argument Dt using the longitudinal information recorded up to time t = Tstart.

The estimate of \mbox{AUC}(t, \Delta t) provided by aucJM() is in the spirit of Harrell's c-index, that is for the comparable subjects (i.e., the ones whose observed event times can be ordered), we compare their dynamic survival probabilities calculated by survfitJM. As with Harrell's c-index, \mbox{AUC}(t, \Delta t) does not fully take into account censoring, and therefore is expected to mask the true discriminative capability of the joint model under heavy censoring.

Value

A list of class aucJM with components:

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Antolini, L., Boracchi, P., and Biganzoli, E. (2005). A time-dependent discrimination index for survival data. Statistics in Medicine 24, 3927–3944.

Harrell, F., Kerry, L. and Mark, D. (1996). Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics in Medicine 15, 361–387.

Heagerty, P. and Zheng, Y. (2005). Survival model predictive accuracy and ROC curves. Biometrics 61, 92–105.

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

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

Rizopoulos, D., Murawska, M., Andrinopoulou, E.-R., Lesaffre, E. and Takkenberg, J. (2013). Dynamic predictions with time-dependent covariates in survival analysis: A comparison between joint modeling and landmarking. under preparation.

See Also

survfitJM, dynCJM, jointModel

Examples

## Not run: 
# we construct the composite event indicator (transplantation or death)
pbc2$status2 <- as.numeric(pbc2$status != "alive")
pbc2.id$status2 <- as.numeric(pbc2.id$status != "alive")

# we fit the joint model using splines for the subject-specific 
# longitudinal trajectories and a spline-approximated baseline
# risk function
lmeFit <- lme(log(serBilir) ~ ns(year, 3),
    random = list(id = pdDiag(form = ~ ns(year, 3))), data = pbc2)
survFit <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
jointFit <- jointModel(lmeFit, survFit, timeVar = "year", 
    method = "piecewise-PH-aGH")

# AUC using data up to year 5 with horizon at year 8
aucJM(jointFit, pbc2, Tstart = 5, Thoriz = 8)

## End(Not run)

Estimated Coefficients for Joint Models

Description

Extracts estimated coefficients from fitted joint models.

Usage

## S3 method for class 'jointModel'
coef(object, process = c("Longitudinal", "Event"), 
    include.splineCoefs = FALSE, ...)
## S3 method for class 'jointModel'
fixef(object, process = c("Longitudinal", "Event"), 
    include.splineCoefs = FALSE, ...)

Arguments

Details

When process = "Event" both methods return the same output. However, for process = "Longitudinal", the coef() method returns the subject-specific coefficients, whereas fixef() only the fixed effects.

Value

A numeric vector or a matrix of the estimated parameters for the fitted model.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

See Also

ranef.jointModel

Examples

## Not run: 
# linear mixed model fit
fitLME <- lme(sqrt(CD4) ~ obstime * drug - drug, 
    random = ~ 1 | patient, data = aids)
# cox model fit
fitCOX <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)

# joint model fit
fitJOINT <- jointModel(fitLME, fitCOX, 
    timeVar = "obstime")

# fixed effects for the longitudinal process
fixef(fitJOINT)

# fixed effects + random effects estimates for the longitudinal 
# process
coef(fitJOINT)

# fixed effects for the event process
fixef(fitJOINT, process = "Event")
coef(fitJOINT, process = "Event")

## End(Not run)

Transform Competing Risks Data in Long Format

Description

In a competing risks setting this function expands the data frame with a single row per subject to the a data frame in long format in which each subject has as many rows as the number of competing events.

Usage

crLong(data, statusVar, censLevel,
    nameStrata = "strata", nameStatus = "status2")

Arguments

Value

A data frame in the long format with multiple rows per subject.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

Putter, H., Fiocco, M., and Geskus, R. (2007). Tutorial in biostatistics: Competing risks and multi-state models. Statistics in Medicine 26, 2389–2430.

Examples

head(crLong(pbc2.id, "status", "alive"))

Derivatives and Integrals of B-splines and Natural Cubic splines

Description

Numerical derivatives and integrals of functions bs() and ns() at their first argument.

Usage

dns(x, df = NULL, knots = NULL, intercept = FALSE, 
    Boundary.knots = range(x), eps = 1e-03)

dbs(x, df = NULL, knots = NULL, intercept = FALSE, 
    Boundary.knots = range(x), eps = 1e-03)

ins(x, df = NULL, knots = NULL, intercept = FALSE, 
    Boundary.knots = range(x), from = 0, weight.fun = NULL, ...)

ibs(x, df = NULL, knots = NULL, intercept = FALSE, 
    Boundary.knots = range(x), from = 0, weight.fun = NULL, ...)

Arguments

Value

an object of class dns, dbs, ins or ibs.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

Examples

x <- rnorm(10)
dns(x, df = 4)
ins(x, df = 4)

A Dynamic Discrimination Index for Joint Models

Description

This function computes a dynamic discrimination index for joint models based on a weighted average of time-dependent AUCs.

Usage

dynCJM(object, newdata, Dt, ...)

## S3 method for class 'jointModel'
dynCJM(object, newdata, Dt, idVar = "id", t.max = NULL, 
    simulate = FALSE, M = 100, weightFun = NULL, ...)

Arguments

Details

(Note: The following contain some math formulas, which are better viewed in the pdf version of the manual accessible at https://cran.r-project.org/package=JM.)

Function dynC computes the following discrimination index

\mbox{C}_{dyn}^{\Delta t} = \int_0^{t_{max}} \mbox{AUC}(t, \Delta t) \, \mbox{Pr} \{ {\cal E}(t, \Delta t) \} \; dt \Big / \int_0^{t_{max}} \mbox{Pr} \{ {\cal E}(t, \Delta t) \} \; dt,

where

\mbox{AUC}(t, \Delta t) = \mbox{Pr} \bigl [ \pi_i(t + \Delta t \mid t) < \pi_j(t + \Delta t \mid t) \mid \{ T_i^* \in (t, t + \Delta t] \} \cap \{ T_j^* > t + \Delta t \} \bigr ],

and

{\cal E}(t, \Delta t) = \bigl [ \{ T_i^* \in (t, t + \Delta t] \} \cap \{ T_j^* > t + \Delta t \} \bigr ],

with i and j denote a randomly selected pair subjects, and \pi_i(t + \Delta t \mid t) and \pi_j(t + \Delta t \mid t) denote the conditional survival probabilities calculated by survfitJM for these two subjects, for different time windows \Delta t specified by argument Dt. The upper limit of integral in specified by argument t.max. The integrals in the numerator and denominator are approximated using a 15-point Gauss-Kronrod quadrature rule.

Index \mbox{C}_{dyn}^{\Delta t} is in the spirit of Harrell's c-index, that is for the comparable subjects (i.e., the ones whose observed event times can be ordered), we compare their dynamic survival probabilities calculated by survfitJM. As with Harrell's c-index, \mbox{C}_{dyn}^{\Delta t} does not take into account censoring, and therefore is expected to mask the true discriminative capability of the joint model under heavy censoring.

Value

A list of class dynCJM with components:

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Antolini, L., Boracchi, P., and Biganzoli, E. (2005). A time-dependent discrimination index for survival data. Statistics in Medicine 24, 3927–3944.

Harrell, F., Kerry, L. and Mark, D. (1996). Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics in Medicine 15, 361–387.

Heagerty, P. and Zheng, Y. (2005). Survival model predictive accuracy and ROC curves. Biometrics 61, 92–105.

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

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

Rizopoulos, D., Murawska, M., Andrinopoulou, E.-R., Lesaffre, E. and Takkenberg, J. (2013). Dynamic predictions with time-dependent covariates in survival analysis: A comparison between joint modeling and landmarking. under preparation.

See Also

survfitJM, aucJM, jointModel

Examples

## Not run: 
# we construct the composite event indicator (transplantation or death)
pbc2$status2 <- as.numeric(pbc2$status != "alive")
pbc2.id$status2 <- as.numeric(pbc2.id$status != "alive")

# we fit the joint model using splines for the subject-specific 
# longitudinal trajectories and a spline-approximated baseline
# risk function
lmeFit <- lme(log(serBilir) ~ ns(year, 3),
    random = list(id = pdDiag(form = ~ ns(year, 3))), data = pbc2)
survFit <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
jointFit <- jointModel(lmeFit, survFit, timeVar = "year", 
    method = "piecewise-PH-aGH")

# dynamic discrimination index up to year 10 using a two-year interval 
dynCJM(jointFit, pbc2, Dt = 2, t.max = 10)

## End(Not run)

Fitted Values for Joint Models

Description

Calculates fitted values for joint models.

Usage

## S3 method for class 'jointModel'
fitted(object, process = c("Longitudinal", "Event"), 
type = c("Marginal", "Subject", "EventTime", "Slope"), scale = c("survival", 
"cumulative-Hazard", "log-cumulative-Hazard"), M = 200, ...)

Arguments

Details

For process = "Longitudinal", let X denote the design matrix for the fixed effects \beta, and Z the design matrix for the random effects b. Then for type = "Marginal" the fitted values are X \hat{\beta}, whereas for type = "Subject" they are X \hat{\beta} + Z \hat{b}. For type = "EventTime" is the same as type = "Subject" but evaluated at the observed event times. Finally, type == "Slope" returns Xs \hat{\beta} + Zs \hat{b} where Xs and Zs denote the fixed- and random-effects design matrices corresponding to the slope term which is specified in the derivForm argument of jointModel.

For process = "Event" and type = "Subject" the linear predictor conditional on the random effects estimates is calculated for each sample unit. Depending on the value of the scale argument the fitted survival function, cumulative hazard function or log cumulative hazard function is returned. For type = "Marginal", random effects values for each sample unit are simulated M times from a normal distribution with zero mean and covariance matrix the estimated covariance matrix for the random effects. The marginal survival function for the ith sample unit is approximated by

S_i(t) = \int S_i(t | b_i) p(b_i) db_i \approx (1/M) \sum_{m = 1}^M S_i(t | b_{im}),

where p(b_i) denotes the normal probability density function, and b_{im} the mth simulated value for the random effect of the ith sample unit. The cumulative hazard and log cumulative hazard functions are calculated as H_i(t) = - \log S_i(t) and \log H_i(t) = \log \{ - \log S_i(t)\}, respectively.

Value

a numeric vector of fitted values.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

Rizopoulos, D. (2010) JM: An R Package for the Joint Modelling of Longitudinal and Time-to-Event Data. Journal of Statistical Software 35 (9), 1–33. doi:10.18637/jss.v035.i09

See Also

residuals.jointModel

Examples

## Not run: 
# linear mixed model fit
fitLME <- lme(log(serBilir) ~ drug * year, 
    random = ~ 1 | id, data = pbc2)
# survival regression fit
fitSURV <- survreg(Surv(years, status2) ~ drug, 
    data = pbc2.id, x = TRUE)
# joint model fit, under the (default) Weibull model
fitJOINT <- jointModel(fitLME, fitSURV, timeVar = "year")

# fitted for the longitudinal process
head(cbind(
    "Marg" = fitted(fitJOINT), 
    "Subj" = fitted(fitJOINT, type = "Subject")
))

# fitted for the event process - survival function
head(cbind(
    "Marg" = fitted(fitJOINT, process = "Ev"), 
    "Subj" = fitted(fitJOINT, process = "Ev", type = "Subject")
))

# fitted for the event process - cumulative hazard function
head(cbind(
    "Marg" = fitted(fitJOINT, process = "Ev", 
        scale = "cumulative-Hazard"), 
    "Subj" = fitted(fitJOINT, process = "Ev", type = "Subject", 
        scale = "cumulative-Hazard")
))

## End(Not run)

Joint Models for Longitudinal and Survival Data

Description

This function fits shared parameter models for the joint modelling of normal longitudinal responses and time-to-event data under a maximum likelihood approach. Various options for the survival model are available.

Usage

jointModel(lmeObject, survObject, timeVar, 
  parameterization = c("value", "slope", "both"),
  method = c("weibull-PH-aGH", "weibull-PH-GH", "weibull-AFT-aGH", 
    "weibull-AFT-GH", "piecewise-PH-aGH", "piecewise-PH-GH", 
    "Cox-PH-aGH", "Cox-PH-GH", "spline-PH-aGH", "spline-PH-GH", 
    "ch-Laplace"),
  interFact = NULL, derivForm = NULL, lag = 0, scaleWB = NULL,
  CompRisk = FALSE, init = NULL, control = list(), ...)

Arguments

Details

Function jointModel fits joint models for longitudinal and survival data (more detailed information about the formulation of these models can be found in Rizopoulos (2010)). For the longitudinal responses the linear mixed effects model represented by the lmeObject is assumed. For the survival times let w_i denote the vector of baseline covariates in survObject, with associated parameter vector \gamma, m_i(t) the value of the longitudinal outcome at time point t as approximated by the linear mixed model (i.e., m_i(t) equals the fixed-effects part + random-effects part of the linear mixed effects model for sample unit i), \alpha the association parameter for m_i(t), m_i'(t) the derivative of m_i(t) with respect to t, and \alpha_d the association parameter for m_i'(t). Then, for method = "weibull-AFT-GH" a time-dependent Weibull model under the accelerated failure time formulation is assumed. For method = "weibull-PH-GH" a time-dependent relative risk model is postulated with a Weibull baseline risk function. For method = "piecewise-PH-GH" a time-dependent relative risk model is postulated with a piecewise constant baseline risk function. For method = "spline-PH-GH" a time-dependent relative risk model is assumed in which the log baseline risk function is approximated using B-splines. For method = "ch-Laplace" an additive model on the log cumulative hazard scale is assumed (see Rizopoulos et al., 2009 for more info). Finally, for method = "Cox-PH-GH" a time-dependent relative risk model is assumed where the baseline risk function is left unspecified (Wulfsohn and Tsiatis, 1997). For all these options the linear predictor for the survival submodel is written as

\eta = \gamma^\top w_i + \alpha m_i\{max(t-k, 0)\},

when parameterization = "value",

\eta = \gamma^\top w_i + \alpha_s m_i'\{max(t-k, 0)\},

when parameterization = "slope", and

\eta = \gamma^\top w_i + \alpha m_i\{max(t-k, 0)\} + \alpha_s m_i'\{max(t-k, 0)\},

when parameterization = "both", where in all the above the value of k is specified by the lag argument and m_i'(t) = d m_i(t) / dt. If interFact is specified, then m_i\{max(t-k, 0)\} and/or m_i'\{max(t-k, 0)\} are multiplied with the design matrices derived from the formulas supplied as the first two arguments of interFact, respectively. In this case \alpha and/or \alpha_s become vectors of association parameters.

For method = "spline-PH-GH" it is also allowed to include stratification factors. These should be included in the specification of the survObject using function strata(). Note that in this case survObject must only be a 'coxph' object.

For all survival models except for the time-dependent proportional hazards model, the optimization algorithm starts with EM iterations, and if convergence is not achieved, it switches to quasi-Newton iterations (i.e., BFGS in optim() or nlminb(), depending on the value of the optimizer control argument). For method = "Cox-PH-GH" only the EM algorithm is used. During the EM iterations, convergence is declared if either of the following two conditions is satisfied: (i) L(\theta^{it}) - L(\theta^{it - 1}) < tol_3 \{ | L(\theta^{it - 1}) | + tol_3 \} , or (ii) \max \{ | \theta^{it} - \theta^{it - 1} | / ( | \theta^{it - 1} | + tol_1) \} < tol_2, where \theta^{it} and \theta^{it - 1} is the vector of parameter values at the current and previous iterations, respectively, and L(.) is the log-likelihood function. The values for tol_1, tol_2 and tol_3 are specified via the control argument. During the quasi-Newton iterations, the default convergence criteria of either optim() or nlminb() are used.

The required integrals are approximated using the standard Gauss-Hermite quadrature rule when the chosen option for the method argument contains the string "GH", and the (pseudo) adaptive Gauss-Hermite rule when the chosen option for the method argument contains the string "aGH". For method = "ch-Laplace" the fully exponential Laplace approximation described in Rizopoulos et al. (2009) is used. The (pseudo) adaptive Gauss-Hermite and the Laplace approximation are particularly useful when high-dimensional random effects vectors are considered (e.g., when modelling nonlinear subject-specific trajectories with splines or high-order polynomials).

Value

See jointModelObject for the components of the fit.

Note

1. The lmeObject argument should represent a linear mixed model object with a simple random-effects structure, i.e., only the pdDiag() class is currently allowed.

2. The lmeObject object should not contain any within-group correlation structure (i.e., correlation argument of lme()) or within-group heteroscedasticity structure (i.e., weights argument of lme()).

3. It is assumed that the linear mixed effects model lmeObject and the survival model survObject have been fitted to the same subjects. Moreover, it is assumed that the ordering of the subjects is the same for both lmeObject and survObject, i.e., that the first line in the data frame containing the event times corresponds to the first set of lines identified by the grouping variable in the data frame containing the repeated measurements, and so on.

4. In the print and summary generic functions for class jointModel, the estimated coefficients (and standard errors for the summary generic) for the event process are augmented with the element "Assoct" that corresponds to the association parameter \alpha and the element "Assoct.s" that corresponds to the parameter \alpha_s when parameterization is "slope" or "both" (see Details).

5. The standard errors returned by the summary generic function for class jointModel when method = "Cox-PH-GH" are based on the profile score vector (i.e., given the NPMLE for the unspecified baseline hazard). Hsieh et al. (2006) have noted that these standard errors are underestimated.

6. As it is the case for all types of mixed models that require numerical integration, it is advisable (especially in difficult datasets) to check the stability of the maximum likelihood estimates with an increasing number of Gauss-Hermite quadrature points.

7. It is assumed that the scale of the time variable (e.g., days, months years) is the same in both lmeObject and survObject.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Henderson, R., Diggle, P. and Dobson, A. (2000) Joint modelling of longitudinal measurements and event time data. Biostatistics 1, 465–480.

Hsieh, F., Tseng, Y.-K. and Wang, J.-L. (2006) Joint modeling of survival and longitudinal data: Likelihood approach revisited. Biometrics 62, 1037–1043.

Rizopoulos, D. (2012a) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

Rizopoulos, D. (2012b) Fast fitting of joint models for longitudinal and event time data using a pseudo-adaptive Gaussian quadrature rule. Computational Statistics and Data Analysis 56, 491–501.

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

Rizopoulos, D. (2010) JM: An R package for the joint modelling of longitudinal and time-to-event data. Journal of Statistical Software 35 (9), 1–33. doi:10.18637/jss.v035.i09

Rizopoulos, D., Verbeke, G. and Lesaffre, E. (2009) Fully exponential Laplace approximations for the joint modelling of survival and longitudinal data. Journal of the Royal Statistical Society, Series B 71, 637–654.

Rizopoulos, D., Verbeke, G. and Molenberghs, G. (2010) Multiple-imputation-based residuals and diagnostic plots for joint models of longitudinal and survival outcomes. Biometrics 66, 20–29.

Tsiatis, A. and Davidian, M. (2004) Joint modeling of longitudinal and time-to-event data: an overview. Statistica Sinica 14, 809–834.

Wulfsohn, M. and Tsiatis, A. (1997) A joint model for survival and longitudinal data measured with error. Biometrics 53, 330–339.

See Also

jointModelObject, anova.jointModel, coef.jointModel, fixef.jointModel, ranef.jointModel, fitted.jointModel, residuals.jointModel, plot.jointModel, survfitJM, rocJM, dynCJM, aucJM, prederrJM

Examples

## Not run: 
# linear mixed model fit (random intercepts)
fitLME <- lme(log(serBilir) ~ drug * year, random = ~ 1 | id, data = pbc2)
# survival regression fit
fitSURV <- survreg(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
# joint model fit, under the (default) Weibull model
fitJOINT <- jointModel(fitLME, fitSURV, timeVar = "year")
fitJOINT
summary(fitJOINT)

# linear mixed model fit (random intercepts + random slopes)
fitLME <- lme(log(serBilir) ~ drug * year, random = ~ year | id, data = pbc2)
# survival regression fit
fitSURV <- survreg(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
# joint model fit, under the (default) Weibull model
fitJOINT <- jointModel(fitLME, fitSURV, timeVar = "year")
fitJOINT
summary(fitJOINT)

# we also include an interaction term of log(serBilir) with drug
fitJOINT <- jointModel(fitLME, fitSURV, timeVar = "year",
    interFact = list(value = ~ drug, data = pbc2.id))
fitJOINT
summary(fitJOINT)


# a joint model in which the risk for and event depends both on the true value of
# marker and the true value of the slope of the longitudinal trajectory
lmeFit <- lme(sqrt(CD4) ~ obstime * drug, random = ~ obstime | patient, data = aids)
coxFit <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)

# to fit this model we need to specify the 'derivForm' argument, which is a list
# with first component the derivative of the fixed-effects formula of 'lmeFit' with
# respect to 'obstime', second component the indicator of which fixed-effects 
# coefficients correspond to the previous defined formula, third component the 
# derivative of the random-effects formula of 'lmeFit' with respect to 'obstime', 
# and fourth component the indicator of which random-effects correspond to the 
# previous defined formula
dForm <- list(fixed = ~ 1 + drug, indFixed = c(2, 4), random = ~ 1, indRandom = 2)
jointModel(lmeFit, coxFit, timeVar = "obstime", method = "spline-PH-aGH",
  parameterization = "both", derivForm = dForm)


# Competing Risks joint model
# we first expand the PBC dataset in the competing risks long format
# with two competing risks being death and transplantation
pbc2.idCR <- crLong(pbc2.id, "status", "alive")

# we fit the linear mixed model as before
lmeFit.pbc <- lme(log(serBilir) ~ drug * ns(year, 3), 
    random = list(id = pdDiag(form = ~ ns(year, 3))), data = pbc2)

# however, for the survival model we need to use the data in the long
# format, and include the competing risks indicator as a stratification
# factor. We also take interactions of the baseline covariates with the
# stratification factor in order to allow the effect of these covariates
# to be different for each competing risk
coxCRFit.pbc <- coxph(Surv(years, status2) ~ (drug + sex)*strata + strata(strata), 
    data = pbc2.idCR, x = TRUE)

# the corresponding joint model is fitted simply by including the above
# two submodels as main arguments, setting argument CompRisk to TRUE, 
# and choosing as method = "spline-PH-aGH". Similarly as above, we also 
# include strata as an interaction factor to allow serum bilirubin to 
# have a different effect for each of the two competing risks
jmCRFit.pbc <- jointModel(lmeFit.pbc, coxCRFit.pbc, timeVar = "year", 
    method = "spline-PH-aGH", 
    interFact = list(value = ~ strata, data = pbc2.idCR), 
    CompRisk = TRUE)
summary(jmCRFit.pbc)

# linear mixed model fit
fitLME <- lme(sqrt(CD4) ~ obstime * drug - drug, 
    random = ~ 1 | patient, data = aids)
# cox model fit
fitCOX <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)
# joint model fit with a spline-approximated baseline hazard function
fitJOINT <- jointModel(fitLME, fitCOX, 
    timeVar = "obstime", method = "spline-PH-aGH")
fitJOINT
summary(fitJOINT)

## End(Not run)

Fitted jointModel Object

Description

An object returned by the jointModel function, inheriting from class jointModel and representing a fitted joint model for longitudinal and time-to-event data. Objects of this class have methods for the generic functions anova, coef, fitted, fixed.effects, logLik, plot, print, random.effects, residuals, summary, and vcov.

Value

The following components must be included in a legitimate jointModel object.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

See Also

jointModel

Mayo Clinic Primary Biliary Cirrhosis Data

Description

Followup of 312 randomised patients with primary biliary cirrhosis, a rare autoimmune liver disease, at Mayo Clinic.

Format

A data frame with 1945 observations on the following 20 variables.

id

patients identifier; in total there are 312 patients.

years

number of years between registration and the earlier of death, transplantion, 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 gm/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.

Note

The data frame pbc2.id contains the first measurement for each patient. This data frame is used to fit the survival model.

References

Fleming, T. and Harrington, D. (1991) Counting Processes and Survival Analysis. Wiley, New York.

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

Examples

summary(pbc2.id)

Proportional Hazards Models with Piecewise Constant Baseline Hazard Function

Description

Based on a fitted Cox model this function fits the corresponding relative risk model with a piecewise constant baseline hazard using the Poisson regression equivalence

Usage

piecewiseExp.ph(coxObject, knots = NULL, length.knots = 6)

Arguments

Value

an object of class glm.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

Examples

coxFit <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)
piecewiseExp.ph(coxFit)

Plot Diagnostics for Joint Models

Description

Produces a variety of plots for fitted joint models.

Usage

## S3 method for class 'jointModel'
plot(x, which = 1:4, caption = c("Residuals vs Fitted", 
  "Normal Q-Q", "Marginal Survival", "Marginal Cumulative Hazard", 
  "Marginal log Cumulative Hazard", "Baseline Hazard", 
  "Cumulative Baseline Hazard", "Subject-specific Survival", 
  "Subject-specific Cumulative Hazard", 
  "Subject-specific log Cumulative Hazard"), survTimes = NULL, 
  main = "", 
  ask = prod(par("mfcol")) < length(which) && dev.interactive(), 
  ..., ids = NULL, add.smooth = getOption("add.smooth"), 
  add.qqline = TRUE, add.KM = FALSE, cex.caption = 1, return = FALSE)

Arguments

Note

The plots of the baseline hazard and the cumulative baseline hazard are only produced when the joint model has been fitted using method = "Cox-PH-GH".

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

Rizopoulos, D. (2010) JM: An R package for the joint modelling of longitudinal and time-to-event data. Journal of Statistical Software 35 (9), 1–33. doi:10.18637/jss.v035.i09

See Also

jointModel

Examples

## Not run: 
# linear mixed model fit
fitLME <- lme(log(serBilir) ~ drug * year, random = ~ 1 | id, data = pbc2)
# survival regression fit
fitSURV <- survreg(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
# joint model fit, under the (default) Weibull model
fitJOINT <- jointModel(fitLME, fitSURV, timeVar = "year")

plot(fitJOINT, 3, add.KM = TRUE, col = "red", lwd = 2)

par(mfrow = c(2, 2))
plot(fitJOINT)

## End(Not run)

Plot Method for rocJM Objects

Description

Produces plots of ROC curves and the corresponding areas under the curve.

Usage

## S3 method for class 'rocJM'
plot(x, which = NULL, type = c("ROC", "AUC"), 
  ndt = "all", main = NULL, caption = NULL, xlab = NULL, 
  ylab = NULL, ask = NULL, legend = FALSE, lx = NULL, ly = NULL, 
  lty = NULL, col = NULL, cex.caption = 0.8, cex.axis = NULL, 
  cex.lab = NULL, cex.main = NULL, ...)

Arguments

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

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

See Also

rocJM

Examples

## Not run: 
fitLME <- lme(sqrt(CD4) ~ obstime + obstime:(drug + AZT + prevOI + gender), 
    random = ~ obstime | patient, data = aids)
fitSURV <- coxph(Surv(Time, death) ~ drug + AZT + prevOI + gender, 
    data = aids.id, x = TRUE)
fit.aids <- jointModel(fitLME, fitSURV, timeVar = "obstime", 
    method = "piecewise-PH-aGH")

ND <- aids[aids$patient == "7", ]
roc <- rocJM(fit.aids, dt = c(2, 4, 8), ND, idVar = "patient")
plot(roc, lwd = 2, legend = TRUE)
plot(roc, type = "AUC")

## End(Not run)

Plot Method for survfitJM Objects

Description

Produces plots of conditional probabilities of survival.

Usage

## S3 method for class 'survfitJM'
plot(x, estimator = c("both", "mean", "median"), 
    which = NULL, fun = NULL, conf.int = FALSE, 
    fill.area = FALSE, col.area = "grey", col.abline = "black", col.points = "black",
    add.last.time.axis.tick = FALSE, include.y = FALSE, main = NULL, 
    xlab = NULL, ylab = NULL, ylab2 = NULL, lty = NULL, col = NULL, 
    lwd = NULL, pch = NULL, ask = NULL, legend = FALSE, ...,
    cex.axis.z = 1, cex.lab.z = 1)

Arguments

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

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

Rizopoulos, D. (2010) JM: An R Package for the Joint Modelling of Longitudinal and Time-to-Event Data. Journal of Statistical Software 35 (9), 1–33. doi:10.18637/jss.v035.i09

See Also

survfitJM

Examples

# linear mixed model fit
fitLME <- lme(sqrt(CD4) ~ obstime + obstime:drug, 
    random = ~ 1 | patient, data = aids)
# cox model fit
fitCOX <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)

# joint model fit
fitJOINT <- jointModel(fitLME, fitCOX, 
    timeVar = "obstime", method = "weibull-PH-aGH")

# sample of the patients who are still alive
ND <- aids[aids$patient == "141", ]
ss <- survfitJM(fitJOINT, newdata = ND, idVar = "patient", M = 50)
plot(ss)
plot(ss, include.y = TRUE, add.last.time.axis.tick = TRUE, legend = TRUE)

Prediction Errors for Joint Models

Description

Using the available longitudinal information up to a starting time point, this function computes an estimate of the prediction error of survival at a horizon time point based on joint models.

Usage

prederrJM(object, newdata, Tstart, Thoriz, ...)

## S3 method for class 'jointModel'
prederrJM(object, newdata, Tstart, Thoriz, 
    lossFun = c("absolute", "square"), interval = FALSE, idVar = "id", 
    simulate = FALSE, M = 100, ...)

Arguments

Details

Based on a fitted joint model (represented by object) and using the data supplied in argument newdata, this function computes the following estimate of the prediction:

PE(u | t) = \{R(t)\}^{-1} \sum_{i: T_i \geq s} I(T_i \geq u) L\{1 - Pr(T_i > u | T_i > t, \tilde{y}_i(t), x_i)\}

+ \delta_i I(T_i < u) L\{0 - Pr(T_i > u | T_i > t, \tilde{y}_i(t), x_i)\}

+ (1 - \delta_i) I(T_i < u) [S_i(u \mid T_i, \tilde{y}_i(t)) L\{1 - Pr(T_i > u | T_i > t, \tilde{y}_i(t), x_i)\}

+ \{1 - S_i(u \mid T_i, \tilde{y}_i(t))\} L\{0 - Pr(T_i > u | T_i > t, \tilde{y}_i(t), x_i)\}],

where R(t) denotes the number of subjects at risk at time t = Tstart, \tilde{y}_i(t) = \{y_i(s), 0 \leq s \leq t\} denotes the available longitudinal measurements up to time t, T_i denotes the observed event time for subject i, \delta_i is the event indicator, s is the starting time point Tstart up to which the longitudinal information is used, and u > s is the horizon time point Thoriz. Function L(.) is the loss function that can be the absolute value (i.e., L(x) = |x|), the squared value (i.e., L(x) = x^2), or a user-specified function. The probabilities Pr(T_i > u | T_i > t, \tilde{y}_i(t), x_i) are calculated by survfitJM.

When interval is set to TRUE, then function prederrJM computes the integrated prediction error in the interval (u,t) = (Tstart, Thoriz) defined as

IPE(u | t) = \sum_{i: t \leq T_i \leq u} w_i(T_i) PE(T_i | t),

where

w_i(T_i) = \frac{\delta_i G(T_i) / G(t)}{\sum_{i: t \leq T_i \leq u} \delta_i G(T_i) / G(t)},

with G(.) denoting the Kaplan-Meier estimator of the censoring time distribution.

Value

A list of class prederrJM with components:

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Henderson, R., Diggle, P. and Dobson, A. (2002). Identification and efficacy of longitudinal markers for survival. Biostatistics 3, 33–50.

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

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

Rizopoulos, D., Murawska, M., Andrinopoulou, E.-R., Lesaffre, E. and Takkenberg, J. (2013). Dynamic predictions with time-dependent covariates in survival analysis: A comparison between joint modeling and landmarking. under preparation.

See Also

survfitJM, aucJM, dynCJM, jointModel

Examples

## Not run: 
# we construct the composite event indicator (transplantation or death)
pbc2$status2 <- as.numeric(pbc2$status != "alive")
pbc2.id$status2 <- as.numeric(pbc2.id$status != "alive")

# we fit the joint model using splines for the subject-specific 
# longitudinal trajectories and a spline-approximated baseline
# risk function
lmeFit <- lme(log(serBilir) ~ ns(year, 3),
    random = list(id = pdDiag(form = ~ ns(year, 3))), data = pbc2)
survFit <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
jointFit <- jointModel(lmeFit, survFit, timeVar = "year", 
    method = "piecewise-PH-aGH")

# prediction error at year 10 using longitudinal data up to year 5 
prederrJM(jointFit, pbc2, Tstart = 5, Thoriz = 10)
prederrJM(jointFit, pbc2, Tstart = 5, Thoriz = 6.5, interval = TRUE)

## End(Not run)

Predictions for Joint Models

Description

Calculates predicted values for the longitudinal part of a joint model.

Usage

## S3 method for class 'jointModel'
predict(object, newdata, type = c("Marginal", "Subject"),
    interval = c("none", "confidence", "prediction"), level = 0.95, idVar = "id", 
    FtTimes = NULL, M = 300, returnData = FALSE, scale = 1.6, ...)

Arguments

Details

When type = "Marginal", this function computes predicted values for the fixed-effects part of the longitudinal submodel. In particular, let X denote the fixed-effects design matrix calculated using newdata. The predict() calculates \hat{y} = X \hat{\beta}, and if interval = "confidence", var(\hat{y}) = X V X^t, with V denoting the covariance matrix of \hat{\beta}. Confidence intervals are constructed under the normal approximation.

When type = "Subject", this functions computes subject-specific predictions for the longitudinal outcome based on the joint model. This accomplished with a Monte Carlo simulation scheme, similar to the one described in survfitJM. The only difference is in Step 3, where for interval = "confidence" y_i^* = X_i \beta^* + Z_i b_i^*, whereas for interval = "prediction" y_i^* is a random vector from a normal distribution with mean X_i \beta^* + Z_i b_i^* and standard deviation \sigma^*. Based on this Monte Carlo simulation scheme we take as estimate of \hat{y}_i the average of the M estimates y_i^* from each Monte Carlo sample. Confidence intervals are constructed using the percentiles of y_i^* from the Monte Carlo samples.

Value

If se.fit = FALSE a numeric vector of predicted values, otherwise a list with components pred the predicted values, se.fit the standard error for the fitted values, and low and upp the lower and upper limits of the confidence interval. If returnData = TRUE, it returns the data frame newdata with the previously mentioned components added.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

See Also

jointModel

Examples

## Not run: 
# linear mixed model fit
fitLME <- lme(log(serBilir) ~ drug * year, 
    random = ~ year | id, data = pbc2)
# survival regression fit
fitSURV <- survreg(Surv(years, status2) ~ drug, 
    data = pbc2.id, x = TRUE)
# joint model fit, under the (default) Weibull model
fitJOINT <- jointModel(fitLME, fitSURV, timeVar = "year")

DF <- with(pbc2, expand.grid(drug = levels(drug),
    year = seq(min(year), max(year), len = 100)))
Ps <- predict(fitJOINT, DF, interval = "confidence", return = TRUE)
require(lattice)
xyplot(pred + low + upp ~ year | drug, data = Ps,
    type = "l", col = c(2,1,1), lty = c(1,2,2), lwd = 2,
    ylab = "Average log serum Bilirubin")


# Subject-specific predictions
ND <- pbc2[pbc2$id == 2, ]
Ps.ss <- predict(fitJOINT, ND, type = "Subject",
  interval = "confidence", return = TRUE)
require(lattice)
xyplot(pred + low + upp ~ year | id, data = Ps.ss,
    type = "l", col = c(2,1,1), lty = c(1,2,2), lwd = 2,
    ylab = "Average log serum Bilirubin")

## End(Not run)

Prednisone versus Placebo in Liver Cirrhosis Patients

Description

A randomized trial on 488 liver cirrhosis patients

Format

Two data frames with the following variable.

id

patients identifier; in total there are 467 patients.

pro

prothrobin measurements.

time

for data frame prothro the time points at which the prothrobin measurements were taken; for data frame prothros the time to death or censoring.

death

a numeric vector with 0 denoting censoring and 1 death.

treat

randomized treatment; a factor with levels "placebo" and "prednisone".

Source

http://www.gllamm.org/books/readme.html#14.6,

References

Andersen, P. K., Borgan, O., Gill, R. D. and Keiding, N. (1993). Statistical Models Based on Counting Processes. New York: Springer.

Examples

summary(prothros)

Random Effects Estimates for Joint Models

Description

Extracts the random effects estimates from a fitted joint model.

Usage

## S3 method for class 'jointModel'
ranef(object, type = c("mean", "mode"), postVar = FALSE, ...)

Arguments

Value

a numeric matrix with rows denoting the individuals and columns the random effects (e.g., intercepts, slopes, etc.). If postVar = TRUE, the numeric matrix has an extra attribute “postVar".

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

See Also

coef.jointModel, fixef.jointModel

Examples

## Not run: 
# linear mixed model fit
fitLME <- lme(log(serBilir) ~ drug * year, random = ~ 1 | id, data = pbc2)
# survival regression fit
fitSURV <- survreg(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)

# joint model fit, under the (default) Weibull model
fitJOINT <- jointModel(fitLME, fitSURV, timeVar = "year")
ranef(fitJOINT)

## End(Not run)

Residuals for Joint Models

Description

Calculates residuals for joint models.

Usage

## S3 method for class 'jointModel'
residuals(object, process = c("Longitudinal", "Event"), 
  type = c("Marginal", "Subject", "stand-Marginal", 
  "stand-Subject", "Martingale", "nullMartingale", "CoxSnell", "AFT"), 
  MI = FALSE, M = 50, time.points = NULL, return.data = FALSE, 
  ...)

Arguments

Details

When process = "Longitudinal", residuals are calculated for the longitudinal outcomes. In particular, if type = "Marginal" these are e_{ij} = y_{ij} - x_{ij}^T \hat{\beta}, whereas for type = "Subject", e_{ij} = y_{ij} - x_{ij}^T \hat{\beta} - z_{ij}^T b_i, where i denotes the subject and j the measurement, y_{ij} the longitudinal responses, x_{ij}^T and z_{ij}^T the corresponding rows of the fixed and random effects design matrices, respectively, and \beta and b_i denote the fixed effects and random effects components. If type = "stand-Marginal" or type = "stand-Subject", the above defined residuals are divided by the estimated standard deviation of the corresponding error term. If MI = TRUE, multiple-imputation-based residuals are calculated for the longitudinal process; for more information regarding these residuals, check Rizopoulos et al. (2009).

When process = "Event", residuals are calculated for the survival outcome. Martingale residuals are available for all options for the survival submodel (for the different options of survival submodel, check the method argument of jointModel). when option type = "nullMartingale" is invoked, the martingale residuals are calculated with the coefficient(s) that correspond to the marker set to zero. Cox-Snell residuals (Cox and Snell, 1968) are available for the Weibull model and the additive log cumulative hazard model. AFT residuals are only available for the Weibull model.

Value

If MI = FALSE, a numeric vector of residual values. Otherwise a list with components:

Note

The multiple-imputation-based residuals are not available for joint models with method = "Cox-PH-GH".

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Cox, D. and Snell, E. (1968) A general definition of residuals. Journal of the Royal Statistical Society, Series B 30, 248–275.

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

Rizopoulos, D., Verbeke, G. and Molenberghs, G. (2010) Multiple-imputation-based residuals and diagnostic plots for joint models of longitudinal and survival outcomes. Biometrics 66, 20–29.

Rizopoulos, D. (2010) JM: An R Package for the Joint Modelling of Longitudinal and Time-to-Event Data. Journal of Statistical Software 35 (9), 1–33. doi:10.18637/jss.v035.i09

See Also

fitted.jointModel

Examples

## Not run: 
# linear mixed model fit
fitLME <- lme(sqrt(CD4) ~ obstime * drug - drug, 
    random = ~ 1 | patient, data = aids)
# cox model fit
fitCOX <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)

# joint model fit, under the additive log cumulative hazard model
fitJOINT <- jointModel(fitLME, fitCOX, 
    timeVar = "obstime")

# residuals for the longitudinal outcome
head(cbind(
    "Marginal" = residuals(fitJOINT),
    "std-Marginal" = residuals(fitJOINT, type = "stand-Marginal"),
    "Subject" = residuals(fitJOINT, type = "Subject"),
    "std-Subject" = residuals(fitJOINT, type = "stand-Subject")
))

# residuals for the survival outcome
head(cbind(
    "Martingale" = residuals(fitJOINT, process = "Event", type = "Martingale"),
    "CoxSnell" = residuals(fitJOINT, process = "Event", type = "CoxSnell")
))

## End(Not run)

Predictive Accuracy Measures for Longitudinal Markers under a Joint Modelling Framework

Description

It computes sensitivity, specificity, ROC and AUC measures for joint models.

Usage

rocJM(object, dt, data, idVar = "id", directionSmaller = NULL, cc = NULL, min.cc = NULL,
  max.cc = NULL, optThr = c("sens*spec", "youden"), 
  diffType = c("absolute", "relative"), abs.diff = 0, rel.diff = 1, 
  M = 300, burn.in = 100, scale = 1.6)

Arguments

Details

(Note: the following contain some math formulas, which are better viewed in the pdf version of the manual accessible at https://cran.r-project.org/package=JM.)

Assume that we have collected longitudinal measurements Y_i(t) = \{y_i(s); 0 \leq s \leq t\} up to time point t for subject i. We are interested in events occurring in the medically relevant time frame (t, t + \Delta t] within which the physician can take an action to improve the survival chance of the patient. Using an appropriate function of the marker history Y_i(t), we can define a prediction rule to discriminate between patients of high and low risk for an event. For instance, for in HIV infected patients, we could consider values of CD4 cell count smaller than a specific threshold as predictive for death. Since we are in a longitudinal context, we have the flexibility of determining which values of the longitudinal history of the patient will contribute to the specification of the prediction rule. That is, we could define a composite prediction rule that is not based only on the last available measurement but rather on the last two or last three measurements of a patient. Furthermore, it could be of relevance to consider different threshold values for each of these measurements, for instance, we could define as success the event that the pre-last CD4 cell count is c and the last one 0.5c, indicating that a 50% decrease is strongly indicative for death. Under this setting we define sensitivity and specificity as,

Pr \bigl \{ {\cal S}_i(t, k, c) \mid T_i^* > t, T_i^* \in (t, t + \Delta t] \bigr \},

and

Pr \bigl \{ {\cal F}_i(t, k, c) \mid T_i^* > t, T_i^* > t + \Delta t \bigr \},

respectively, where we term {\cal S}_i(t, k, c) = \{y_i(s) \leq c_s; k \leq s \leq t\} as success (i.e., occurrence of the event of interest), and {\cal F}_i(t, k, c) = \{y_i(s) > c_s; k \leq s \leq t\} as a failure, T_i^* denotes the time-to-event, and \Delta t the length of the medically relevant time window (specified by argument dt). The cut values for the marker c are specified by the cc, min.cc and max.cc arguments. Two types of composite prediction rules can be defined depending on the value of the diffType argument. Absolute prediction rules in which, between successive measurements there is an absolute difference of between the cut values, and relative prediction rules in which there is a relative difference between successive measurements of the marker. The lag values for these differences are defined by the abs.diff and rel.diff arguments. Some illustrative examples:

Ex1:

keeping the defaults we define a simple rule that is only based on the last available marker measurement.

Ex2:

to define a prediction rule that is based on the last two available measurements using the same cut values (e.g., if a patient had two successive measurements below a medically relevant threshold), we need to set abs.diff = c(0, 0).

Ex3:

to define a prediction rule that is based on the last two available measurements using a drop of 5 units between the cut values (e.g., the pre-last measurement is c and the last c-5), we need to set abs.diff = c(0, -5).

Ex4:

to define a prediction rule that is based on the last two available measurements using a drop of 20% units between the cut values (e.g., the pre-last measurement is c and the last 0.8c), we need to set diffType = "relative" and rel.diff = c(0, 0.8).

The estimation of the above defined probabilities is achieved with a Monte Carlo scheme similar to the one described in survfitJM. The number of Monte Carlo samples is defined by the M argument, and the burn-in iterations for the Metropolis-Hastings algorithm using the burn.in argument.

More details can be found in Rizopoulos (2011).

Value

An object of class rocJM is a list with components,

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Heagerty, P. and Zheng, Y. (2005). Survival model predictive accuracy and ROC curves. Biometrics 61, 92–105.

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

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

Rizopoulos, D. (2010) JM: An R package for the joint modelling of longitudinal and time-to-event data. Journal of Statistical Software 35 (9), 1–33. doi:10.18637/jss.v035.i09

Zheng, Y. and Heagerty, P. (2007). Prospective accuracy for longitudinal markers. Biometrics 63, 332–341.

See Also

plot.rocJM, survfitJM, dynCJM, aucJM, prederrJM, jointModel

Examples

## Not run: 
fitLME <- lme(sqrt(CD4) ~ obstime * (drug + AZT + prevOI + gender), 
    random = ~ obstime | patient, data = aids)
fitSURV <- coxph(Surv(Time, death) ~ drug + AZT + prevOI + gender, 
    data = aids.id, x = TRUE)
fit.aids <- jointModel(fitLME, fitSURV, timeVar = "obstime", 
    method = "piecewise-PH-aGH")

# the following will take some time to execute...
ND <- aids[aids$patient == "7", ]
roc <- rocJM(fit.aids, dt = c(2, 4, 8), ND, idVar = "patient")
roc

## End(Not run)

Simulate from Joint Models.

Description

simulate longitudinal responses and event times from joint models

Usage

simulateJM(nsim, nsub, thetas, times, formulas, Data = NULL, 
    method = c("weibull-PH", "weibull-AFT", "piecewise-PH", "spline-PH"), 
    lag = 0, censoring = "uniform", max.FUtime = NULL, seed = NULL, 
    return.ranef = FALSE)
## S3 method for class 'jointModel'
simulate(object, nsim, seed = NULL, times = NULL, 
    Data = NULL, ...)

Arguments

Value

A list of length nsim of data frames that contains the simulated responses for the longitudinal process "y", the simulated event times "Time", the event indicator "Event", and the subject identification number "id". If extra covariates were assumed, these are also included.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

See Also

jointModel

Examples

## Not run: 
prothro$t0 <- as.numeric(prothro$time == 0)
lmeFit <- lme(pro ~ treat * (time + t0), random = ~ time | id, data = prothro)
survFit <- coxph(Surv(Time, death) ~ treat, data = prothros, x = TRUE)
jointFit <- jointModel(lmeFit, survFit, timeVar = "time", 
    method = "weibull-PH-aGH")

newData <- simulate(jointFit, nsim = 1, times = seq(0, 11, len = 15))
newData

## End(Not run)

Summary Method for weibull.frailty Objects

Description

Summarizes the fit of a Weibull model with Gamma frailties

Usage

## S3 method for class 'weibull.frailty'
summary(object, sand.se = FALSE, ...)

Arguments

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

See Also

weibull.frailty

Examples

fit <- weibull.frailty(Surv(time, status) ~ age + sex, kidney)
summary(fit)
summary(fit, TRUE)

Prediction in Joint Models

Description

This function computes the conditional probability of surviving later times than the last observed time for which a longitudinal measurement was available.

Usage

survfitJM(object, newdata, ...)

## S3 method for class 'jointModel'
survfitJM(object, newdata, idVar = "id", simulate = TRUE, survTimes = NULL, 
	last.time = NULL, M = 200, CI.levels = c(0.025, 0.975), scale = 1.6, ...)

Arguments

Details

Based on a fitted joint model (represented by object), and a history of longitudinal responses \tilde{y}_i(t) = \{y_i(s), 0 \leq s \leq t\} and a covariates vector x_i (stored in newdata), this function provides estimates of Pr(T_i > u | T_i > t, \tilde{y}_i(t), x_i), i.e., the conditional probability of surviving time u given that subject i, with covariate information x_i, has survived up to time t and has provided longitudinal the measurements \tilde{y}_i(t).

To estimate Pr(T_i > u | T_i > t, \tilde{y}_i(t), x_i) and if simulate = TRUE, a Monte Carlo procedure is followed with the following steps:

Step 1:

Simulate new parameter values, say \theta^*, from N(\hat{\theta}, C(\hat{\theta})), where \hat{\theta} are the MLEs and C(\hat{\theta}) their large sample covariance matrix, which are extracted from object.

Step 2:

Simulate random effects values, say b_i^*, from their posterior distribution given survival up to time t, the vector of longitudinal responses \tilde{y}_i(t) and \theta^*. This is achieved using a Metropolis-Hastings algorithm with independent proposals from a properly centered and scaled multivariate t distribution. The scale argument controls the acceptance rate for this algorithm.

Step 3

Using \theta^* and b_i^*, compute Pr(T_i > u | T_i > t, b_i^*, x_i; \theta^*).

Step 4:

Repeat Steps 1-3 M times.

Based on the M estimates of the conditional probabilities, we compute useful summary statistics, such as their mean, median, and quantiles (to produce a confidence interval).

If simulate = FALSE, then survival probabilities are estimated using the formula

Pr(T_i > u | T_i > t, \hat{b}_i, x_i; \hat{\theta}),

where \hat{\theta} denotes the MLEs as above, and \hat{b}_i denotes the empirical Bayes estimates.

Value

A list of class survfitJM with components:

Note

Predicted probabilities are not computed for joint models with method = "ch-Laplace" and method = "Cox-PH-GH".

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

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

Rizopoulos, D. (2010) JM: An R Package for the Joint Modelling of Longitudinal and Time-to-Event Data. Journal of Statistical Software 35 (9), 1–33. doi:10.18637/jss.v035.i09

See Also

jointModel, plot.survfitJM

Examples

# linear mixed model fit
fitLME <- lme(sqrt(CD4) ~ obstime + obstime:drug, 
    random = ~ 1 | patient, data = aids)
# cox model fit
fitCOX <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)

# joint model fit
fitJOINT <- jointModel(fitLME, fitCOX, 
    timeVar = "obstime", method = "weibull-PH-aGH")

# sample of the patients who are still alive
ND <- aids[aids$patient == "141", ]
ss <- survfitJM(fitJOINT, newdata = ND, idVar = "patient", M = 50)
ss

Wald Test for Stratification Factors

Description

It performs a Wald test to test the hypothesis of equal spline coefficients among strata in the approximation of baseline risk function.

Usage

wald.strata(fit)

Arguments

Value

an object of class wald.strata with components:

Note

This test is valid when the same knots have been used across strata.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.

Examples

## Not run: 
fitLME <- lme(log(serBilir) ~ drug * year - drug, random = ~ year | id, 
    data = pbc2)
fitSURV <- coxph(Surv(years, status2) ~ drug + strata(hepatomegaly), 
    data = pbc2.id, x = TRUE)
fit.pbc <- jointModel(fitLME, fitSURV, timeVar = "year", method = "spline-PH-aGH")
wald.strata(fit.pbc)

## End(Not run)

Weibull Model with Gamma Frailties

Description

Fits a Weibull model with Gamma frailties for multivariate survival data under maximum likelihood

Usage

weibull.frailty(formula = formula(data), data = parent.frame(), 
    id = "id", subset, na.action, init, control = list())

Arguments

Details

The fitted model is defined as follows:

\lambda(t_i | \omega_i) = \lambda_0(t_i) \omega_i \exp(x_i^T \beta),

where i denotes the subject, \lambda(\cdot) denotes the hazard function, conditionally on the frailty \omega_i, x_i is a vector of covariates with corresponding regression coefficients \beta, and \lambda_0(\cdot) is the Weibull baseline hazard defined as \lambda_0(t) = shape * scale * t^{shape -1}. Finally, for the frailties we assume \omega_i \sim Gamma(\eta, \eta), with \eta^{-1} denoting the unknown variance of \omega_i's.

Value

an object of class weibull.frailty with components:

Note

weibull.frailty() currently supports only right-censored data.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

Examples

weibull.frailty(Surv(time, status) ~ age + sex, kidney)

xtable Method from Joint Models.

Description

produces a LaTeX table with the results of a joint model using package xtable.

Usage

## S3 method for class 'jointModel'
xtable(x, caption = NULL, label = NULL, align = NULL,
    digits = NULL, display = NULL, which = c("all", "Longitudinal", "Event"),
    varNames.Long = NULL, varNames.Event = NULL, p.values = TRUE, 
    digits.pval = 4, ...)

Arguments

Value

A LaTeX code chunk with the results of the joint modeling analysis.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

See Also

jointModel

Examples

## Not run: 
require(xtable)
prothro$t0 <- as.numeric(prothro$time == 0)
lmeFit <- lme(pro ~ treat * (time + t0), random = ~ time | id, data = prothro)
survFit <- coxph(Surv(Time, death) ~ treat, data = prothros, x = TRUE)
jointFit <- jointModel(lmeFit, survFit, timeVar = "time", 
    method = "weibull-PH-aGH")

xtable(jointFit, math.style.negative = TRUE)

## End(Not run)