Joint Modeling of Longitudinal and Time-to-Event Data in R under a Bayesian Approach

Description

This package fits shared parameter models for the joint modeling of normal longitudinal responses and event times under a Bayesian approach. Various options for the survival model and the association structure are provided.

Details

The package has a single model-fitting function called jointModelBayes, which accepts as main arguments a linear mixed effects object fit returned by function lme() of package nlme, and a Cox model object fit returned by function coxph() of package survival. The survMod argument of specifies the type of survival submodel to be fitted; available options are a relative risk model with a Weibull baseline hazard (default) and a relative risk model with a B-spline approximation of the log baseline risk function. In addition, the param specifies the association structure between the longitudinal and survival processes; available options are: "td-value" which is the classic formulation used in Wulfsohn and Tsiatis (1997); "td-extra" which is a user-defined, possibly time-dependent, term based on the specification of the extraForm argument of jointModelBayes. This could be used to include terms, such as the time-dependent slope (i.e., the derivative of the subject-specific linear predictor of the linear mixed model) and the time-dependent cumulative effect (i.e., the integral of the subject-specific linear predictor of the linear mixed model); "td-both" which is the combination of the previous two parameterizations, i.e., the current value and the user-specified terms are included in the linear predictor of the relative risk model; and "shared-RE" where only the random effects of the linear mixed model are included in the linear predictor of the survival submodel.

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

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

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

Rizopoulos, D. (2016). The R package JMbayes for fitting joint models for longitudinal and time-to-event data using MCMC. Journal of Statistical Software 72(7), 1–45. doi:10.18637/jss.v072.i07.

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. and Ghosh, P. (2011) A Bayesian semiparametric multivariate joint model for multiple longitudinal outcomes and a time-to-event. Statistics in Medicine 30, 1366–1380.

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

jointModelBayes, survfitJM, aucJM, dynCJM, prederrJM, predict.JMbayes, logLik.JMbayes

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, 
    integrand.fun = NULL, ...)

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

Arguments

Value

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

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

Individualized Predictions from Linear Mixed Models

Description

Calculates subject-specific predictions for new subjects from a linear mixed model.

Usage

IndvPred_lme(lmeObject, newdata, timeVar, times = NULL, M = 200L,
    interval = c("confidence", "prediction"), all_times = FALSE,
    level = 0.95, return_data = FALSE, seed = 1L)
    
extract_lmeComponents(lmeObject, timeVar)

Arguments

Value

If return_data = TRUE, a the data frame newdata with extra rows for the time points at which predictions were calculated, and extra columns with the predictions and the limits of the pointwise confidence intervals.

If return_data = FALSE, a list with components

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

See Also

predict.JMbayes

Examples

## Not run: 
# linear mixed model fit
fitLME <- lme(log(serBilir) ~ drug * ns(year, 2), data = subset(pbc2, id != 2), 
        random = ~ ns(year, 2) | id)
        
DF <- IndvPred_lme(fitLME, newdata = subset(pbc2, id == 2), timeVar = "year", 
    M = 500, return_data = TRUE)
    
require(lattice)
xyplot(pred + low + upp ~ year | id, data = DF,
    type = "l", col = c(2,1,1), lty = c(1,2,2), lwd = 2,
    ylab = "Average log serum Bilirubin")
    
# extract_lmeComponents() extract the required components from the lme object
# that are required to calculate the predictions; this is a light weight version of
# the object, e.g.,
fitLME_light <- extract_lmeComponents(fitLME, timeVar = "year")

DF <- IndvPred_lme(fitLME_light, newdata = subset(pbc2, id == 2), timeVar = "year", 
    M = 500, return_data = TRUE)

## End(Not run)

Fitted JMbayes Object

Description

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

Value

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

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

See Also

jointModelBayes

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.

Anova Method for Fitted Joint Models

Description

Comparison of (non)nested joint models using information criteria.

Usage

## S3 method for class 'JMbayes'
anova(object, ...)

Arguments

Value

A data frame with rows the different models, and columns the number or parameters in each model, the the log pseudo marginal likelihood value, the deviance information criterion value, and the pD value.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

See Also

jointModelBayes

Examples

## Not run: 
# composite event indicator
pbc2$status2 <- as.numeric(pbc2$status != "alive")
pbc2.id$status2 <- as.numeric(pbc2.id$status != "alive")

# linear mixed model with natural cubic splines for the time
# effect
lmeFit.pbc1 <- lme(log(serBilir) ~ ns(year, 2), data = pbc2,
                   random = ~ ns(year, 2) | id, method = "ML")

# Cox regression model with baseline covariates
coxFit.pbc1 <- coxph(Surv(years, status2) ~ drug * age, data = pbc2.id, x = TRUE)

# the standard joint model fit with only the m_i(t) term in 
# the linear predictor of the survival submodel
jointFit.pbc1 <- jointModelBayes(lmeFit.pbc1, coxFit.pbc1, timeVar = "year")

# we include the time-dependent slopes term
dForm <- list(fixed = ~ 0 + dns(year, 2), random = ~ 0 + dns(year, 2), 
              indFixed = 2:3, indRandom = 2:3)

jointFit.pbc2 <- update(jointFit.pbc1, param = "td-both", extraForm = dForm)

# we include the cumulative effect of the marker
iForm <- list(fixed = ~ 0 + year + ins(year, 2), random = ~ 0 + year + ins(year, 2), 
              indFixed = 1:3, indRandom = 1:3)

jointFit.pbc3 <- update(jointFit.pbc1, param = "td-extra", extraForm = iForm)


# we compare the three models
anova(jointFit.pbc1, jointFit.pbc2, jointFit.pbc3)

## End(Not run)

Time-Dependent ROCs and AUCs for Joint Models

Description

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

Usage

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

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

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

rocJM(object, newdata, Tstart, ...)

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

## S3 method for class 'mvJMbayes'
rocJM(object, newdata, Tstart, Thoriz = NULL, 
    Dt = NULL, idVar = "id", M = 100, ...)
    
predict_eventTime(object, newdata, cut_points, ...)

## S3 method for class 'mvJMbayes'
predict_eventTime(object, newdata, cut_points, 
    idVar = "id", M = 500L, low_percentile = 0.025, ...)

find_thresholds(object, newdata, Dt, ...)

## S3 method for class 'mvJMbayes'
find_thresholds(object, newdata, Dt, idVar = "id", 
    M = 200L, variability_threshold = NULL,
    n_cores =  max(1, parallel::detectCores() - 2), ...)

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. For the subjects who due to censoring we do not know if they are comparable, they contribute in the AUC with the probability that they would have been comparable.

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. (2016). The R package JMbayes for fitting joint models for longitudinal and time-to-event data using MCMC. Journal of Statistical Software 72(7), 1–45. doi:10.18637/jss.v072.i07.

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

survfitJM, dynCJM, jointModelBayes

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 <- jointModelBayes(lmeFit, survFit, timeVar = "year")

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

plot(rocJM(jointFit, pbc2, Tstart = 5, Thoriz = 8))

## End(Not run)

Combines Predictions for Bayesian Model Averaging

Description

Combines estimated survival probabilities or predictions for longitudinal responses.

Usage

bma.combine(..., JMlis = NULL, weights = NULL)

Arguments

Value

an object of class survfit.JMbayes or predict.JMbayes.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2016). The R package JMbayes for fitting joint models for longitudinal and time-to-event data using MCMC. Journal of Statistical Software 72(7), 1–45. doi:10.18637/jss.v072.i07.

Rizopoulos, D., Hatfield, L., Carlin, B. and Takkenberg, J. (2014). Combining dynamic predictions from joint models for longitudinal and time-to-event data using Bayesian model averaging. Journal of the American Statistical Association 109, 1385–1397.

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 two joint models using splines for the subject-specific 
# longitudinal trajectories and a spline-approximated baseline
# risk function; the first one with the current value parameterization
# and the other with the shared random effects parameterization
lmeFit <- lme(log(serBilir) ~ ns(year, 2), data = pbc2,
              random = ~ ns(year, 2) | id)
survFit <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)

jointFit1 <- jointModelBayes(lmeFit, survFit, timeVar = "year")
jointFit2 <- update(jointFit1, param = "shared-RE")

# we compute survival probabilities for Subject 2 with 
# different weights
ND <- pbc2[pbc2$id == 2, ] # the data of Subject 2
survPreds1 <- survfitJM(jointFit1, newdata = ND, weight = 0.4)
survPreds2 <- survfitJM(jointFit2, newdata = ND, weight = 0.6)

survPreds.bma <- bma.combine(survPreds1, survPreds2)
survPreds.bma
plot(survPreds.bma)

## End(Not run)

Estimated Coefficients and Confidence Intervals for Joint Models

Description

Extracts estimated coefficients and confidence intervals from fitted joint models.

Usage

## S3 method for class 'JMbayes'
coef(object, process = c("Longitudinal", "Event"), ...)

## S3 method for class 'JMbayes'
fixef(object, process = c("Longitudinal", "Event"), ...)

## S3 method for class 'JMbayes'
confint(object, parm = c("all", "Longitudinal", "Event"), ...)

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 or confidence intervals for the fitted model.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

See Also

ranef.JMbayes, jointModelBayes

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 <- jointModelBayes(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)

Dynamic Information

Description

Using the available longitudinal information up to a starting time point, this function computes an estimate of the cross-entropy function based on joint models.

Usage

cvDCL(object, newdata, Tstart, idVar = "id", M = 300L, seed = 123L)

Arguments

Details

This function calculates an estimate of the cross-entropy at any time point t (given in Tstart) that can be used to identify the joint model that best predicts future event giver survival up to t.

Value

A list of class aucJM with components:

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

See Also

survfitJM, dynCJM, jointModelBayes

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 <- jointModelBayes(lmeFit, survFit, timeVar = "year")

cvDCL(jointFit, Tstart = 5)

## End(Not run)

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 'JMbayes'
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=JMbayes.)

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.

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. (2016). The R package JMbayes for fitting joint models for longitudinal and time-to-event data using MCMC. Journal of Statistical Software 72(7), 1–45. doi:10.18637/jss.v072.i07.

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

survfitJM, aucJM, jointModelBayes

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, 2), data = pbc2,
              random = ~ ns(year, 2) | id)
survFit <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
jointFit <- jointModelBayes(lmeFit, survFit, timeVar = "year")

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

## End(Not run)

Dynamic Information of an Extra Longitudinal Measurement

Description

Using the available longitudinal information up to a particular time point, this function computes an estimate of the information we again by obtaining an extra longitudinal measurement.

Usage

dynInfo(object, newdata, Dt, K = 5, M = 500, idVar = "id", 
    simulateFun = function (eta, scale) rnorm(length(eta), eta, scale), 
    seed = 1L)

Arguments

Details

This functions computes the following posterior predictive distribution

E_{Y} [ E_{T^* | Y} (\log p (T_j^* \mid T_j^* > u, \{ Y_j(t), y_j(u) \}, D_n \bigr )) ],

where T_j^* denotes the time-to-event for subject j for whom we wish to plan the next visit, Y_j(t) the available longitudinal measurements of this subject up to time t, y_j(u) the future longitudinal measurement we wish to plan at time u > t, and D_n the data set that was used to fit the joint model.

Value

A list with components:

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

See Also

survfitJM, jointModelBayes

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 <- jointModelBayes(lmeFit, survFit, timeVar = "year")

dynInfo(jointFit, newdata = pbc2[pbc2$id == 2, ], Dt = 5)[[1]]

## End(Not run)

Fitted Values and Residuals for Joint Models

Description

Calculates fitted values for joint models.

Usage

## S3 method for class 'JMbayes'
fitted(object, 
    process = c("Longitudinal", "longitudinal", "Event", "event"), 
    type = c("Marginal", "marginal", "Subject", "subject"), nullY = FALSE, ...)

## S3 method for class 'JMbayes'
residuals(object, 
    process = c("Longitudinal", "longitudinal", "Event", "event"), 
    type = c("Marginal", "marginal", "Subject", "subject", 
             "Martingale", "martingale", "nullMartingale", "nullmartingale"),
    standardized = FALSE, ...)

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}, where \hat{\beta} and \hat{b} denote the corresponding posterior means for the fixed and random effects. The corresponding residuals are calculated by subtracting the fitted values from the observed data y. If type = "Subject" and standardized = TRUE, the residuals are divided by the estimated residual standard error.

For process = "Event" function fitted() calculates the cumulative hazard function at each time point a longitudinal measurement has been recorded. If nullY = TRUE, then the cumulative hazard is calculated without the contribution of the longitudinal process. Function residuals() calculates the martingales residuals or the martingale residuals without the contribution of the longitudinal process when type = "nullMartingale".

Value

a numeric vector of fitted values or residuals.

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: 
lmeFit <- lme(log(serBilir) ~ ns(year, 2), data = pbc2,
              random = ~ ns(year, 2) | id)
survFit <- coxph(Surv(years, status2) ~ 1, data = pbc2.id, x = TRUE)
jointFit <- jointModelBayes(lmeFit, survFit, timeVar = "year")

fitted(jointFit, process = "Event")
residuals(jointFit, type = "Subject", standardized = TRUE)

## End(Not run)

The Generalized Student's t Distribution

Description

Density, distribution function, quantile function and random generation for the generalized Student's t distribution.

Usage

dgt(x, mu = 0, sigma = 1, df = stop("no df arg"), log = FALSE)

pgt(q, mu = 0, sigma = 1, df = stop("no df arg")) 

qgt(p, mu = 0, sigma = 1, df = stop("no df arg"))

rgt(n, mu = 0, sigma = 1, df = stop("no df arg"))

Arguments

Value

dgt gives the density, pgt gives the distribution function, qgt gives the quantile function, and rgt generates random deviates.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

Examples

## Not run: 
x <- rnorm(10, mean = 10, sd = 3)
dgt(x, mu = 10, sigma = 3, df = 4)
rgt(10, mu = 10, sigma = 3, df = 4)

## End(Not run)

Joint Models for Longitudinal and Time-to-Event Data

Description

Fits shared parameter joint models for longitudinal and survival outcomes under a Bayesian approach.

Usage

jointModelBayes(lmeObject, survObject, timeVar,  
    param = c("td-value", "td-extra", "td-both", "shared-betasRE", "shared-RE"), 
    extraForm = NULL, baseHaz = c("P-splines", "regression-splines"), 
    transFun = NULL, densLong = NULL, lag = 0, df.RE = NULL, 
    estimateWeightFun = FALSE, weightFun = NULL, init = NULL, 
    priors = NULL, scales = NULL, control = list(), ...)

Arguments

Details

Function jointModelBayes fits joint models for longitudinal and survival data under a Bayesian approach. For the longitudinal responses a linear mixed effects model represented by the lmeObject is assumed, unless the user specifies his own probability density function using argument densLong. For the survival times, let w_i denote the vector of baseline covariates in survObject, with associated parameter vector \gamma, m_i(t) the subject-specific linear predictor at time point t as defined by the mixed model (i.e., m_i(t) equals the fixed-effects part + random-effects part of the mixed effects model for sample unit i), m_i'(t) denotes an extra user-defined term (based on the specification of argument extraForm) to be included in the linear predictor of the survival submodel, and \alpha and \alpha_e vector of association parameters for m_i(t) and m_i'(t), respectively. Then, jointModelBayes assumes a relative risk model of the form

h_i(t) = h_0(t) \exp\{\gamma^\top w_i + f(m_i(t), m_i'(t), b_i; \alpha, \alpha_e)\},

where the baseline risk function is approximated using splines, i.e.,

\log h_0(t) = \sum_k \tilde\gamma_k B(t; \lambda),

with B(.) denoting a B-spline basis function, \lambda a vector of knots, and \tilde \gamma_k the corresponding splines coefficients (\tilde \gamma correspond to Bs.gammas above). Argument baseHaz specifies whether a penalized- or regression-spline-approximation is employed. For the former the P-splines approach of Eilers and Marx (1996) is used, namely the prior for \tilde \gamma is taken to be proportional to

p(\tilde \gamma) \propto \exp \Bigl(- \frac{\tau_{bs}}{2} \tilde \gamma^\top \Delta^\top \Delta \tilde \gamma \Bigr),

where \Delta denotes the differences matrix (the order of the differences is set by the control argument diff).

Function f(m_i(t), m_i'(t), b_i; \alpha, \alpha_d) specifies the association structure between the two processes. In particular, for param = "td-value",

f(m_i(t), m_i'(t), b_i; \alpha, \alpha_d) = f_1(m_i(t)) \alpha,

for param = "td-extra",

f(m_i(t), m_i'(t), b_i; \alpha, \alpha_d) = f_2(m_i'(t)) \alpha_e,

for param = "td-both",

f(m_i(t), m_i'(t), b_i; \alpha, \alpha_d) = f_1(m_i(t)) \alpha + f_2(m_i'(t)) \alpha_e,

for param = "shared-RE",

f(m_i(t), m_i'(t), b_i; \alpha, \alpha_d) = \alpha^\top b_i,

and for param = "shared-betasRE",

f(m_i(t), m_i'(t), b_i; \alpha, \alpha_d) = \alpha^\top (\beta^* + b_i),

where f_1(.) and f_2(.) denote possible transformation functions, b_i denotes the vector of random effects for the ith subject and \beta^* the fixed effects that correspond to the random effects.

Value

See JMbayesObject for the components of the fit.

Note

1. The lmeObject argument should represent a mixed model object without any special structure in the random-effects covariance matrix (i.e., no use of pdMats).

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. Furthermore, the scale of the time variable (e.g., days, months, years) should be the same in both the lmeObject and survObject objects.

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 "AssoctE" that corresponds to the parameter \alpha_e when parameterization is "td-extra" or "td-both" (see Details).

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. (2016). The R package JMbayes for fitting joint models for longitudinal and time-to-event data using MCMC. Journal of Statistical Software 72(7), 1–45. doi:10.18637/jss.v072.i07.

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.

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

coef.JMbayes, ranef.JMbayes, logLik.JMbayes, survfitJM, aucJM, dynCJM, prederrJM, predict.JMbayes

Examples

## Not run: 
# A joint model for the AIDS dataset:
# First we fit the linear mixed model for the longitudinal measurements of
# sqrt CD4 cell counts
lmeFit.aids <- lme(CD4 ~ obstime * drug, random = ~ obstime | patient, data = aids)
# next we fit the Cox model for the time to death (including the 'x = TRUE' argument)
survFit.aids <- coxph(Surv(Time, death) ~ drug, data = aids.id, x = TRUE)

# the corresponding joint model is fitted by (the default is to assume 
# the current value parameterization)
jointFit.aids <- jointModelBayes(lmeFit.aids, survFit.aids, timeVar = "obstime")
summary(jointFit.aids)

# A joint model for the PBC dataset:
# We first fit the linear mixed and Cox models. In the first we include 
# splines to model flexibly the subject-specific longitudinal trajectories
lmeFit.pbc <- lme(log(serBilir) ~ ns(year, 2),
    random = list(id = pdDiag(form = ~ ns(year, 2))), data = pbc2)
survFit.pbc <- coxph(Surv(years, status2) ~ 1, data = pbc2.id, x = TRUE)

# the corresponding joint model is fitted by:
jointFit.pbc <- jointModelBayes(lmeFit.pbc, survFit.pbc, timeVar = "year", 
    baseHaz = "regression-splines")
summary(jointFit.pbc)

# we update the joint model fitted for the PBC dataset by including
# the time-dependent slopes term. To achieve this we need to define 
# the 'extraForm' argument, in which we use function dns() to numerically
# compute the derivative of the natural cubic spline. In addition, we increase
# the number of MCMC iterations to 35000
dform = list(fixed = ~ 0 + dns(year, 2), random = ~ 0 + dns(year, 2),
    indFixed = 2:3, indRandom = 2:3)
jointFit.pbc2 <- update(jointFit.pbc, param = "td-both", extraForm = dform,
    n.iter = 35000)
summary(jointFit.pbc2)

# we fit the same model with the shared random effects formulation
jointFit.pbc3 <- update(jointFit.pbc, param = "shared-betasRE")
summary(jointFit.pbc3)

# a joint model for left censored longitudinal data
# we create artificial left censoring in serum bilirubin
pbc2$CensInd <- as.numeric(pbc2$serBilir <= 0.8)
pbc2$serBilir2 <- pbc2$serBilir
pbc2$serBilir2[pbc2$CensInd == 1] <- 0.8

censdLong <- function (y, eta.y, scale, log = FALSE, data) {
    log.f <- dnorm(x = y, mean = eta.y, sd = scale, log = TRUE)
    log.F <- pnorm(q = y, mean = eta.y, sd = scale, log.p = TRUE)
    ind <- data$CensInd
    if (log) {
        (1 - ind) * log.f + ind * log.F
    } else {
        exp((1 - ind) * log.f + ind * log.F)
    }
}
lmeFit.pbc2 <- lme(log(serBilir2) ~ ns(year, 2), data = pbc2,
                   random = ~ ns(year, 2) | id, method = "ML")
jointFit.pbc4 <- jointModelBayes(lmeFit.pbc2, survFit.pbc, timeVar = "year",
                                  densLong = censdLong)

summary(jointFit.pbc4)

# a joint model for a binary outcome
pbc2$serBilirD <- 1 * (pbc2$serBilir > 1.8)
lmeFit.pbc3 <- glmmPQL(serBilirD ~ year, random = ~ year | id, 
	family = binomial, data = pbc2)

dLongBin <- function (y, eta.y, scale, log = FALSE, data) {
    dbinom(x = y, size = 1, prob = plogis(eta.y), log = log)
}

jointFit.pbc5 <- jointModelBayes(lmeFit.pbc3, survFit.pbc, timeVar = "year", 
	densLong = dLongBin)

summary(jointFit.pbc5)


# create start-stop counting process variables
pbc <- pbc2[c("id", "serBilir", "drug", "year", "years",
              "status2", "spiders")]
pbc$start <- pbc$year
splitID <- split(pbc[c("start", "years")], pbc$id)
pbc$stop <- unlist(lapply(splitID,
                          function (d) c(d$start[-1], d$years[1]) ))
pbc$event <- with(pbc, ave(status2, id,
                           FUN = function (x) c(rep(0, length(x)-1), x[1])))
pbc <- pbc[!is.na(pbc$spiders), ]

# left-truncation
pbc <- pbc[pbc$start != 0, ] 

lmeFit.pbc <- lme(log(serBilir) ~ drug * ns(year, 2),
                  random = ~ ns(year, 2) | id, data = pbc)

tdCox.pbc <- coxph(Surv(start, stop, event) ~ drug * spiders + cluster(id),
                   data = pbc, x = TRUE, model = TRUE)

jointFit.pbc6 <- jointModelBayes(lmeFit.pbc, tdCox.pbc, timeVar = "year")

summary(jointFit.pbc6)

## End(Not run)

Log-Likelihood for Joint Models

Description

Computes the log-likelihood for a fitted joint model.

Usage

## S3 method for class 'JMbayes'
logLik(object, thetas, b, priors = TRUE, marginal.b = TRUE, 
    marginal.thetas = FALSE, full.Laplace = FALSE, useModes = TRUE, ...)

Arguments

Details

Let y_i denote the vectors of longitudinal responses, T_i the observed event time, and \delta_i the event indicator for subject i (i = 1, \ldots, n). Let also p(y_i | b_i; \theta) denote the probability density function (pdf) for the linear mixed model, p(T_i, \delta_i | b_i; \theta) the pdf for the survival submodel, and p(b_i; \theta) the multivariate normal pdf for the random effects, where \theta denotes the full parameter vector. Then, if priors = TRUE, and marginal.b = TRUE, function logLik() computes

\log \int p(y_i | b_i; \theta) p(T_i, \delta_i | b_i; \theta) p(b_i; \theta) db_i + \log p(\theta),

where p(\theta) denotes the prior distribution for the parameters. If priors = FALSE the prior is excluded from the computation, i.e.,

\log \int p(y_i | b_i; \theta) p(T_i, \delta_i | b_i; \theta) p(b_i; \theta) db_i,

and when marginal.b = FALSE, then the conditional on the random effects log-likelihood is computed, i.e.,

\log p(y_i | b_i; \theta) + \log p(T_i, \delta_i | b_i; \theta) + \log p(b_i; \theta) + \log p(\theta),

when priors = TRUE and

\log p(y_i | b_i; \theta) + \log p(T_i, \delta_i | b_i; \theta) + \log p(b_i; \theta),

when priors = FALSE.

Value

a numeric scalar of class logLik with the value of the log-likelihood. It also has the attributes df the number of parameter (excluding the random effects), and nobs the number of subjects.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D., Hatfield, L., Carlin, B. and Takkenberg, J. (2014). Combining dynamic predictions from joint models for longitudinal and time-to-event data using Bayesian model averaging. Journal of the American Statistical Association. to appear.

See Also

jointModelBayes

Examples

## Not run: 
lmeFit <- lme(log(serBilir) ~ ns(year, 2), data = pbc2, 
    random = ~ ns(year, 2) | id)
survFit <- coxph(Surv(years, status2) ~ 1, data = pbc2.id, x = TRUE)

jointFit <- jointModelBayes(lmeFit, survFit, timeVar = "year")

logLik(jointFit)
logLik(jointFit, priors = FALSE)
logLik(jointFit, marginal.b = FALSE)

## End(Not run)

Calculates Marginal Subject-specific Log-Likelihood Contributions

Description

This function computes marginal subject-specific log-likelihood contributions based on a fitted joint model. The marginalization is done with respect to both the random effects and the parameters using a Laplace approximation.

Usage

marglogLik(object, newdata, idVar = "id", method = "BFGS", control = NULL)

Arguments

Value

a numeric vector of marginal log-likelihood contributions.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

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 a joint model using splines for the subject-specific 
# longitudinal trajectories and a spline-approximated baseline
# risk function
lmeFit <- lme(log(serBilir) ~ ns(year, 2), data = pbc2,
    random = ~ ns(year, 2) | id)
survFit <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)

jointFit <- jointModelBayes(lmeFit, survFit, timeVar = "year")

# we compute marginal log-likelihood contribution for Subject 2
ND <- pbc2[pbc2$id == 2, ] # the data of Subject 2
marglogLik(jointFit, newdata = ND)

## End(Not run)

Multivariate Joint Models for Longitudinal and Time-to-Event Data

Description

Fits multivariate shared parameter joint models for longitudinal and survival outcomes under a Bayesian approach.

Usage

mvJointModelBayes(mvglmerObject, survObject, timeVar,
    Formulas = list(NULL), Interactions = list(NULL),
    transFuns = NULL, priors = NULL, multiState = FALSE, 
    data_MultiState = NULL, idVar_MultiState = "id", 
    control = NULL, ...)

Arguments

Details

The mathematical details regarding the definition of the multivariate joint model, and the capabilities of the package can be found in the vignette in the doc directory.

Value

A list of class mvJMbayes with components:

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. (2016). The R package JMbayes for fitting joint models for longitudinal and time-to-event data using MCMC. Journal of Statistical Software 72(7), 1–45. doi:10.18637/jss.v072.i07.

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.

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

mvglmer, jointModelBayes

Examples

## Not run: 
pbc2.id$Time <- pbc2.id$years
pbc2.id$event <- as.numeric(pbc2.id$status != "alive")

##########################################################################################

##############################################
# Univariate joint model for serum bilirubin #
##############################################

# [1] Fit the mixed model using mvglmer(). The main arguments of this function are:
# 'formulas' a list of lme4-like formulas (a formular per outcome),
# 'data' a data.frame that contains all the variables specified in 'formulas' (NAs allowed),
# 'families' a list of family objects specifying the type of each outcome (currently only
# gaussian, binomial and poisson are allowed).
MixedModelFit1 <- mvglmer(list(log(serBilir) ~ year + (year | id)), data = pbc2,
                          families = list(gaussian))

# [2] Fit a Cox model, specifying the baseline covariates to be included in the joint
# model; you need to set argument 'model' to TRUE.
CoxFit1 <- coxph(Surv(Time, event) ~ drug + age, data = pbc2.id, model = TRUE)

# [3] The basic joint model is fitted using a call to mvJointModelBayes(), which is very
# similar to JointModelBayes(), i.e.,
JMFit1 <- mvJointModelBayes(MixedModelFit1, CoxFit1, timeVar = "year")
summary(JMFit1)
plot(JMFit1)

##########################################################################################

#########################################################
# Bivariate joint model for serum bilirubin and spiders #
#########################################################

MixedModelFit2 <- mvglmer(list(log(serBilir) ~ year + (year | id),
                               spiders ~ year + (1 | id)), data = pbc2,
                          families = list(gaussian, binomial))

CoxFit2 <- coxph(Surv(Time, event) ~ drug + age, data = pbc2.id, model = TRUE)

JMFit2 <- mvJointModelBayes(MixedModelFit2, CoxFit2, timeVar = "year")
summary(JMFit2)
plot(JMFit2)

##########################################################################################

#######################
# slopes & area terms #
#######################

# We extend model 'JMFit2' by including the value and slope term for bilirubin, and
# the area term for spiders (in the log-odds scale). To include these terms into the model
# we specify the 'Formulas' argument. This is specified in a similar manner as the
# 'derivForms' argument of jointModelBayes(). In particular, it should be a list of lists.
# Each component of the outer list should have as name the name of the corresponding
# outcome variable. Then in the inner list we need to specify four components, namely,
# 'fixed' & 'random' R formulas specifying the fixed and random effects part of the term
# to be included, and 'indFixed' & 'indRandom' integer indicices specifying which of the
# original fixed and random effects are involved in the claculation of the new term. In
# the inner list you can also optionally specify a name for the term you want to include.
# Notes: (1) For terms not specified in the 'Formulas' list, the default value functional
# form is used. (2) If for a particular outcome you want to include both the value
# functional form and an extra term, then you need to specify that in the 'Formulas'
# using two entries. To include the value functional form you only need to set the
# corresponding to 'value', and for the second term to specify the inner list. See
# example below on how to include the value and slope for serum bilirubin (for example,
# if the list below the entry '"log(serBilir)" = "value"' was not give, then only the
# slope term would have been included in the survival submodel).

Forms <- list("log(serBilir)" = "value",
              "log(serBilir)" = list(fixed = ~ 1, random = ~ 1,
                                     indFixed = 2, indRandom = 2, name = "slope"),
              "spiders" = list(fixed = ~ 0 + year + I(year^2/2), random = ~ 0 + year,
                               indFixed = 1:2, indRandom = 1, name = "area"))

JMFit3 <- update(JMFit2, Formulas = Forms)
summary(JMFit3)
plot(JMFit3)

##########################################################################################

#####################
# Interaction terms #
#####################

# We further extend the previous model by including the interactions terms between the
# terms specified in 'Formulas' and 'drug'. The names specified in the list that defined
# the interaction factors should match with the names of the output from 'JMFit3'; the
# only exception is with regard to the 'value' functional form. See specification below
# (to include the interaction of the value term of 'log(serBilir)' with 'drug', in the
# list we can either specify as name of the corresponding formula 'log(serBilir)_value'
# or just 'log(serBilir)'):

Ints <- list("log(serBilir)" = ~ drug, "log(serBilir)_slope" = ~ drug,
             "spiders_area" = ~ drug)

# because of the many association parameters we have, we place a shrinkage prior on the
# alpha coefficients. In particular, if we have K association parameters, we assume that
# alpha_k ~ N(0, tau * phi_k), k = 1, ..., K. The precision parameters tau and phi_k are
# given Gamma priors. Precision tau is global shrinkage parameter, and phi_k a specific
# per alpha coefficient.
JMFit4 <- update(JMFit3, Interactions = Ints, priors = list(shrink_alphas = TRUE))
summary(JMFit4)
plot(JMFit4)

##########################################################################################

########################
# Time-varying effects #
########################

# We allow the association parameter to vary with time; the time-varying coefficients are
# approximated with P-splines
JMFit_tveffect <- mvJointModelBayes(MixedModelFit1, CoxFit1, timeVar = "year",
                    Interactions = list("log(serBilir)_value" = ~ 0 + tve(Time, df = 8)))

plot(JMFit_tveffect, "tv_effect")


##########################################################################################

############################
# Interval censoring terms #
############################

# create artificial interval censoring in the PBC data set
pbc2$status2 <- as.numeric(pbc2$status != "alive")
pbc2.id$status2 <- as.numeric(pbc2.id$status != "alive")
pbc2$status3 <- as.character(pbc2$status)
ff <- function (x) {
    out <- if (x[1L] %in% c('dead', 'transplanted')) 'dead' else 
        switch(sample(1:3, 1), '1' = "alive", '2' = "left", '3' = "interval")
    rep(out, length(x))
}
pbc2$status3 <- unlist(with(pbc2, lapply(split(status3, id), ff)), use.names = FALSE)
pbc2$status3 <- unname(with(pbc2, sapply(status3, function (x) 
    switch(x, 'dead' = 1, 'alive' = 0, 'left' = 2, 'interval' = 3))))
pbc2$yearsU <- as.numeric(NA)
pbc2$yearsU[pbc2$status3 == 3] <- pbc2$years[pbc2$status3 == 3] + 
    runif(sum(pbc2$status3 == 3), 0, 4)
pbc2.id <- pbc2[!duplicated(pbc2$id), ]

# next we fit a weibull model for interval censored data
survFit <- survreg(Surv(years, yearsU, status3, type = "interval") ~ drug + age, 
                   data = pbc2.id, model = TRUE)

# next we fit the joint model (we use 'MixedModelFit1' from above)
JMFit_intcens <- mvJointModelBayes(MixedModelFit1, survFit, timeVar = "year")
summary(JMFit_intcens)

## End(Not run)

Multivariate Mixed Models

Description

Fits multivariate mixed models under a Bayesian approach using JAGS.

Usage

mvglmer(formulas, data, families, engine = c("JAGS", "STAN"), 
    overdispersion = FALSE, priors = NULL, init = NULL, control = NULL, ...)

Arguments

Details

This function creates a JAGS program representing a multivariate mixed effects that is run with JAGS using the jagsUI package. Currently only Gaussian, Bernoulli and Poisson longitudinal outcomes can be handled.

Value

A list of class mvglmer with components:

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

See Also

mvJointModelBayes, jointModelBayes

Examples

## Not run: 
MixedModelFit <- mvglmer(list(log(serBilir) ~ year + (year | id),
                               spiders ~ year + (1 | id)), data = pbc2,
                          families = list(gaussian, binomial))

summary(MixedModelFit)
plot(MixedModelFit)

## End(Not run)

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.

MCMC Diagnostics for Joint Models

Description

Produces MCMC diagnostics plots.

Usage

## S3 method for class 'JMbayes'
plot(x, which = c("trace", "autocorr", "density", "CPO", "weightFun"), 
    param = c("betas", "sigma", "D", "gammas", "alphas", "Dalphas", 
        "shapes", "Bs.gammas", "tauBs"), ask = TRUE, max.t = NULL, 
    from = 0, ...)

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.

See Also

jointModelBayes

Examples

## Not run: 
# linear mixed model fit
fitLME <- lme(log(serBilir) ~ drug * year, random = ~ 1 | id, data = pbc2)

# survival regression fit
fitSURV <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)

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

plot(fitJOINT)

## End(Not run)

Plot Method for survfit.JMbayes and survfit.mvJMbayes Objects

Description

Produces plots of conditional probabilities of survival.

Usage

## S3 method for class 'survfit.JMbayes'
plot(x, estimator = c("both", "mean", "median"), 
    which = NULL, fun = NULL, invlink = 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, xlim = NULL)
    
## S3 method for class 'survfit.mvJMbayes'
plot(x, split = c(1, 1), which_subjects = NULL, 
    which_outcomes = NULL, surv_in_all = TRUE, include.y = TRUE, fun = NULL,
    abline = NULL,
    main = NULL, xlab = "Time", ylab = NULL, zlab = "Event-Free Probability",
    include_CI = TRUE, fill_area_CI = TRUE, col_points = "black", 
    pch_points = 1, col_lines = "red", col_lines_CI = "black", 
    col_fill_CI = "lightgrey", lwd_lines = 2, lty_lines_CI = 2, 
    cex_xlab = 1, cex_ylab = 1, cex_zlab = 1, cex_main = 1, 
    cex_axis = 1, ...)

Arguments

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2016). The R package JMbayes for fitting joint models for longitudinal and time-to-event data using MCMC. Journal of Statistical Software 72(7), 1–45. doi:10.18637/jss.v072.i07.

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

survfitJM

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, 2), data = pbc2,
    random = ~ ns(year, 2) | id)
survFit <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
jointFit <- jointModelBayes(lmeFit, survFit, timeVar = "year")

# we will compute survival probabilities for Subject 2 in a dynamic manner, 
# i.e., after each longitudinal measurement is recorded
ND <- pbc2[pbc2$id == 2, ] # the data of Subject 2
survPreds <- vector("list", nrow(ND))
for (i in 1:nrow(ND)) {
    survPreds[[i]] <- survfitJM(jointFit, newdata = ND[1:i, ])
}

# the default call to the plot method using the first three
# longitudinal measurements
plot(survPreds[[3]])

# we produce the corresponding plot
par(mfrow = c(2, 2), oma = c(0, 2, 0, 2))
for (i in c(1,3,5,7)) {
    plot(survPreds[[i]], estimator = "median", conf.int = TRUE,
        include.y = TRUE, main = paste("Follow-up time:",
            round(survPreds[[i]]$last.time, 1)), ylab = "", ylab2 = "")
}
mtext("log serum bilirubin", side = 2, line = -1, outer = TRUE)
mtext("Survival Probability", side = 4, line = -1, outer = TRUE)

## End(Not run)

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 'JMbayes'
prederrJM(object, newdata, Tstart, Thoriz, 
    lossFun = c("square", "absolute"), interval = FALSE, idVar = "id", 
    simulate = FALSE, M = 100, ...)

## S3 method for class 'mvJMbayes'
prederrJM(object, newdata, Tstart, Thoriz, 
    lossFun = c("square", "absolute"), interval = FALSE, idVar = "id", 
    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. (2016). The R package JMbayes for fitting joint models for longitudinal and time-to-event data using MCMC. Journal of Statistical Software 72(7), 1–45. doi:10.18637/jss.v072.i07.

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

survfitJM, aucJM, dynCJM, jointModelBayes

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, 2), data = pbc2,
    random = ~ ns(year, 2) | id)
survFit <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
jointFit <- jointModelBayes(lmeFit, survFit, timeVar = "year")

# 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 'JMbayes'
predict(object, newdata, type = c("Marginal", "Subject"),
    interval = c("none", "confidence", "prediction"), level = 0.95, idVar = "id", 
    FtTimes = NULL, last.time = NULL, LeftTrunc_var = NULL, 
    M = 300, returnData = FALSE, scale = 1.6, 
    weight = rep(1, nrow(newdata)), invlink = NULL, seed = 1, ...)

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", then it calculates the confidence intervals based on the percentiles of the MCMC sample for \beta.

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.

Note

The user is responsible to appropriately set the invlink argument when a user-specified mixed effects model has been fitted.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2016). The R package JMbayes for fitting joint models for longitudinal and time-to-event data using MCMC. Journal of Statistical Software 72(7), 1–45. doi:10.18637/jss.v072.i07.

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

See Also

survfitJM.JMbayes, jointModelBayes

Examples

## Not run: 
# linear mixed model fit
fitLME <- lme(log(serBilir) ~ drug * year, data = pbc2,
    random = ~ year | id)
# survival regression fit
fitSURV <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, 
    x = TRUE)
# joint model fit, under the (default) Weibull model
fitJOINT <- jointModelBayes(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)
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.

Random Effects Estimates for Joint Models

Description

Extracts the random effects estimates from a fitted joint model.

Usage

## S3 method for class 'JMbayes'
ranef(object, 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.JMbayes, jointModelBayes

Examples

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

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

## End(Not run)

Shiny Application for Dynamic Predictions

Description

This function loads the shiny package and runs the application for calculating dynamic predictions using package JMbayes.

Usage

runDynPred(type = c("JM", "lme"))

Arguments

Value

No value returned.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

Examples

## Not run: 
runDynPred()

## End(Not run)

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 'JMbayes'
survfitJM(object, newdata, 
   type = c("SurvProb", "Density"), idVar = "id", 
   simulate = TRUE, survTimes = NULL, last.time = NULL, 
   LeftTrunc_var = NULL, M = 200L, 
   CI.levels = c(0.025, 0.975), log = FALSE, scale = 1.6, 
   weight = rep(1, nrow(newdata)), 
   init.b = NULL, seed = 1L, ...)

## S3 method for class 'mvJMbayes'
survfitJM(object, newdata, 
    survTimes = NULL, idVar = "id", last.time = NULL, 
    M = 200L, scale = 1.6, log = FALSE, 
    CI.levels = c(0.025, 0.975), seed = 1L, ...)

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:

Take randomly a realization, say \theta^* from the MCMC sample of posterior of the joint model represented by 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 percentiles (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 posterior means for the parameters, and \hat{b}_i denotes the posterior means for the random effects.

Value

A list of class survfit.JMbayes with components:

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

Rizopoulos, D. (2016). The R package JMbayes for fitting joint models for longitudinal and time-to-event data using MCMC. Journal of Statistical Software 72(7), 1–45. doi:10.18637/jss.v072.i07.

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

plot.survfit.JMbayes, predict.JMbayes, aucJM, dynCJM, prederrJM, jointModelBayes

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, 2), data = pbc2,
    random = ~ ns(year, 2) | id)
survFit <- coxph(Surv(years, status2) ~ drug, data = pbc2.id, x = TRUE)
jointFit <- jointModelBayes(lmeFit, survFit, timeVar = "year")

# we will compute survival probabilities for Subject 2 in a dynamic manner, 
# i.e., after each longitudinal measurement is recorded
ND <- pbc2[pbc2$id == 2, ] # the data of Subject 2
survPreds <- vector("list", nrow(ND))
for (i in 1:nrow(ND)) {
    survPreds[[i]] <- survfitJM(jointFit, newdata = ND[1:i, ])
}
survPreds

###########################################################################

# Predictions from multivariate models

pbc2 <- pbc2[!is.na(pbc2$serChol), ]
pbc2.id <- pbc2[!duplicated(pbc2$id), ]
pbc2.id$Time <- pbc2.id$years
pbc2.id$event <- as.numeric(pbc2.id$status != "alive")

# Fit a trivariate joint model
MixedModelFit <- mvglmer(list(log(serBilir) ~ year + (year | id),
                              sqrt(serChol) ~ year + (year | id),
                              hepatomegaly ~ year + (year | id)), data = pbc2,
                         families = list(gaussian, gaussian, binomial), engine = "STAN")

CoxFit <- coxph(Surv(Time, event) ~ drug + age, data = pbc2.id, model = TRUE)

JMFit <- mvJointModelBayes(MixedModelFit, CoxFit, timeVar = "year")

# We want survival probabilities for three subjects
ND <- pbc2[pbc2$id %in% c(2, 25, 81), ]

sprobs <- survfitJM(JMFit, ND)
sprobs

# Basic plot
plot(sprobs)

# split in a 2 rows 2 columns and include the survival function in 
# a separate panel; plot only the third & first subjects; change various defaults
plot(sprobs, split = c(3, 2), surv_in_all = FALSE, which_subjects = c(3, 1),
     lty_lines_CI = 3, col_lines = "blue", col_fill_CI = "red", 
     col_points = "pink", pch_points = 12)

###########################################################################

# run Shiny app
runDynPred()

## End(Not run)

Time-Varying Effects using P-splines

Description

A B-spline expansion of the input variables to be used for a time-varying effect in the specification of joint model.

Usage

tve(x, df = NULL, knots = NULL, ord = 3)

Arguments

Value

an object of class tve.

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

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 'JMbayes'
xtable(x, caption = NULL, label = NULL, align = NULL,
    digits = NULL, display = NULL, which = c("all", "Longitudinal", "Event"),
    varNames.Long = NULL, varNames.Event = NULL, ...)

Arguments

Value

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

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

See Also

jointModelBayes

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 <- jointModelBayes(lmeFit, survFit, timeVar = "time")
    
if (require("xtable")) {
    xtable:::xtable(jointFit, math.style.negative = TRUE)
}

## End(Not run)