| Title: | Joint Modeling of Longitudinal and Time-to-Event Data under a Bayesian Approach |
|---|---|
| Description: | Shared parameter models for the joint modeling of longitudinal and time-to-event data using MCMC; Dimitris Rizopoulos (2016) <doi:10.18637/jss.v072.i07>. |
| Authors: | Dimitris Rizopoulos <[email protected]> |
| Maintainer: | Dimitris Rizopoulos <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.8-86 |
| Built: | 2024-11-07 03:26:56 UTC |
| Source: | https://github.com/drizopoulos/jmbayes |
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.
A data frame with 1408 observations on the following 9 variables.
patientpatients identifier; in total there are 467 patients.
Timethe time to death or censoring.
deatha numeric vector with 0 denoting censoring and 1 death.
CD4the CD4 cells count.
obstimethe time points at which the CD4 cells count was recorded.
druga factor with levels ddC denoting zalcitabine and ddI denoting didanosine.
gendera factor with levels female and male.
prevOIa factor with levels AIDS denoting previous opportunistic infection (AIDS
diagnosis) at study entry, and noAIDS denoting no previous infection.
AZTa factor with levels intolerance and failure denoting AZT intolerance and
AZT failure, respectively.
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.
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.
Comparison of (non)nested joint models using information criteria.
## S3 method for class 'JMbayes' anova(object, ...)## S3 method for class 'JMbayes' anova(object, ...)
object, ...
|
objects inheriting from class |
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.
Dimitris Rizopoulos [email protected]
## 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)## 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)
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.
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), ...)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), ...)
object |
an object inheriting from class |
newdata |
a data frame that contains the longitudinal and covariate information for the subjects for which prediction
of survival probabilities is required. The names of the variables in this data frame must be the same as in the data frames that
were used to fit the linear mixed effects model (using |
Tstart |
numeric scalar denoting the time point up to which longitudinal information is to be used to derive predictions. |
Thoriz |
numeric scalar denoting the time point for which a prediction of the survival status is of interest;
|
Dt |
numeric scalar denoting the length of the time interval of prediction; either |
idVar |
the name of the variable in |
simulate |
logical; if |
M |
a numeric scalar denoting the number of Monte Carlo samples; see |
cut_points |
a numeric matrix with first column time-points followed by other columns of optimal cut-points from an ROC curve. |
variability_threshold |
numeric value denoting the treshold in the spread of the posterior distribution calculated from the 2.5% percentile to the median. Default is the 25% percentile of the event times distribution. |
low_percentile |
a numeric value indicating the percentile based on which it will be judged whether the spread of the posterior predictive distribution is too large. |
n_cores |
an integer indicating the number of cores to use for parallel computing. |
... |
additional arguments; currently none is used. |
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:
with and denote a randomly selected pair of subjects, and
and denote the conditional survival probabilities calculated by
survfitJM for these two subjects, for different time windows specified by argument Dt using
the longitudinal information recorded up to time t = Tstart.
The estimate of provided by aucJM() is in the spirit of Harrell's
-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.
A list of class aucJM with components:
auc |
a numeric scalar denoting the estimated prediction error. |
Tstart |
a copy of the |
Thoriz |
a copy of the |
nr |
a numeric scalar denoting the number of subjects at risk at time |
classObject |
the class of |
nameObject |
the name of |
Dimitris Rizopoulos [email protected]
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.
survfitJM, dynCJM, jointModelBayes
## 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)## 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 estimated survival probabilities or predictions for longitudinal responses.
bma.combine(..., JMlis = NULL, weights = NULL)bma.combine(..., JMlis = NULL, weights = NULL)
... |
objects inheriting from class |
JMlis |
a list of |
weights |
a numeric vector of weights to be applied in each object. |
an object of class survfit.JMbayes or predict.JMbayes.
Dimitris Rizopoulos [email protected]
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.
## 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)## 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)
Extracts estimated coefficients and confidence intervals from fitted joint models.
## 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"), ...)## 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"), ...)
object |
an object inheriting from class |
process |
for which submodel (i.e., linear mixed model or survival model) to extract the estimated coefficients. |
parm |
for which submodel (i.e., linear mixed model or survival model) to extract credible intervals. |
... |
additional arguments; currently none is used. |
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.
A numeric vector or a matrix of the estimated parameters or confidence intervals for the fitted model.
Dimitris Rizopoulos [email protected]
ranef.JMbayes, jointModelBayes
## 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)## 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)
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.
cvDCL(object, newdata, Tstart, idVar = "id", M = 300L, seed = 123L)cvDCL(object, newdata, Tstart, idVar = "id", M = 300L, seed = 123L)
object |
an object inheriting from class |
newdata |
a data frame that contains the longitudinal and covariate information for the subjects for which prediction
of survival probabilities is required. The names of the variables in this data frame must be the same as in the data frames that
were used to fit the linear mixed effects model (using |
Tstart |
a numeric scalar indicating at which time to compute the cross-entropy. |
idVar |
the name of the variable in |
M |
a numeric scalar denoting the number of Monte Carlo samples. |
seed |
a numeric scalar |
This function calculates an estimate of the cross-entropy at any time point (given in Tstart) that can be
used to identify the joint model that best predicts future event giver survival up to .
A list of class aucJM with components:
auc |
a numeric scalar denoting the estimated prediction error. |
Tstart |
a copy of the |
Thoriz |
a copy of the |
nr |
a numeric scalar denoting the number of subjects at risk at time |
classObject |
the class of |
nameObject |
the name of |
Dimitris Rizopoulos [email protected]
survfitJM, dynCJM, jointModelBayes
## 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)## 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)
Numerical derivatives and integrals of functions bs() and ns() at their first argument.
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, ...)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, ...)
x, df, knots, intercept, Boundary.knots
|
see the help pages of functions |
eps |
a numeric scalar denoting the step length for the central difference approximation, which calculates the derivative. |
from |
a numeric scalar denoting the lower limit of the integral. |
weight.fun |
a function to be applied as weights. |
integrand.fun |
a function to be applied in the integrand. |
... |
extra arguments passed to |
an object of class dns, dbs, ins or ibs.
Dimitris Rizopoulos [email protected]
This function computes a dynamic discrimination index for joint models based on a weighted average of time-dependent AUCs.
dynCJM(object, newdata, Dt, ...) ## S3 method for class 'JMbayes' dynCJM(object, newdata, Dt, idVar = "id", t.max = NULL, simulate = FALSE, M = 100, weightFun = NULL, ...)dynCJM(object, newdata, Dt, ...) ## S3 method for class 'JMbayes' dynCJM(object, newdata, Dt, idVar = "id", t.max = NULL, simulate = FALSE, M = 100, weightFun = NULL, ...)
object |
an object inheriting from class |
newdata |
a data frame that contains the longitudinal and covariate information for the subjects for which prediction
of survival probabilities is required. The names of the variables in this data frame must be the same as in the data frames that
were used to fit the linear mixed effects model (using |
Dt |
a numeric scalar denoting the time frame within which the occurrence of events is of interest. |
idVar |
the name of the variable in |
t.max |
a numeric scalar denoting the time maximum follow-up time up to which the dynamic discrimination index is to be calculated.
If |
simulate |
logical; if |
M |
a numeric scalar denoting the number of Monte Carlo samples; see |
weightFun |
a function of two arguments the first denoting time and the second the length of the time frame of interest, i.e., |
... |
additional arguments; currently none is used. |
(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
where
and
with and denote a randomly selected pair subjects, and
and denote the conditional survival probabilities calculated by
survfitJM for these two subjects, for different time windows 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.
A list of class dynCJM with components:
dynC |
a numeric scalar denoting the dynamic discrimination index. |
times |
a numeric vector of time points at which the AUC was calculated. |
AUCs |
a numeric vector of the estimated AUCs at the aforementioned time points. |
weights |
a numeric vector of the estimated weights at the aforementioned time points. |
t.max |
a copy of the |
Dt |
a copy of the |
classObject |
the class of |
nameObject |
the name of |
Dimitris Rizopoulos [email protected]
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.
survfitJM, aucJM, jointModelBayes
## 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)## 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)
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.
dynInfo(object, newdata, Dt, K = 5, M = 500, idVar = "id", simulateFun = function (eta, scale) rnorm(length(eta), eta, scale), seed = 1L)dynInfo(object, newdata, Dt, K = 5, M = 500, idVar = "id", simulateFun = function (eta, scale) rnorm(length(eta), eta, scale), seed = 1L)
object |
an object inheriting from class |
newdata |
a data frame that contains the longitudinal and covariate information for
the subject for whom we wish to plan the next measurement. The names of the variables in
this data frame must be the same as in the data frames that were used to fit the linear
mixed effects model (using |
Dt |
numeric scalar denoting the length of the time interval to search for the
optimal time point of the next measurement, i.e., the interval is |
K |
numeric scalar denoting the number of time points to cosider in the interval
|
idVar |
the name of the variable in |
simulateFun |
a function based on which longitudinal measurement can be simulated.
This should have as a main argument the variable |
M |
a numeric scalar denoting the number of Monte Carlo samples. |
seed |
a numeric scalar |
This functions computes the following posterior predictive distribution
where denotes the time-to-event for subject for
whom we wish to plan the next visit, the available longitudinal measurements
of this subject up to time , the future longitudinal measurement we
wish to plan at time , and the data set that was used to fit the
joint model.
A list with components:
summary |
a numeric matrix with first column the time points at which the longitudinal measurement is hypothetically taken, second column the estimated information we gain by obtaining the measurement, and third column the estimated cumulative risk of an event up to the particular time point denoted in the first column. |
full.results |
a numeric matrix with columns representing the time points, rows the Monte Carlo samples, and entries the value of log posterio predictive density. |
Dimitris Rizopoulos [email protected]
## 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)## 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)
Calculates fitted values for joint models.
## 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, ...)## 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, ...)
object |
an object inheriting from class |
process |
for which model (i.e., linear mixed model or survival model) to calculate fitted values or residuals. |
type |
what type of fitted values or residuals to calculate. See Details. |
nullY |
logical; if |
standardized |
logical; if |
... |
additional arguments; currently none is used. |
For process = "Longitudinal", let denote the design matrix for the fixed effects , and
the design matrix for the random effects . Then for type = "Marginal" the fitted values are
whereas for type = "Subject" they are , where
and 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 . 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".
a numeric vector of fitted values or residuals.
Dimitris Rizopoulos [email protected]
Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.
## 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)## 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)
Density, distribution function, quantile function and random generation for the generalized Student's t distribution.
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"))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"))
x, q
|
vector of quantiles. |
p |
vector of probabilities. |
n |
a numeric scalar denoting the number of observations. |
mu |
a vector of means. |
sigma |
a vector of standard deviations. |
log |
logical; if |
df |
a numeric scalar denoting the degrees of freedom. |
dgt gives the density, pgt gives the distribution function,
qgt gives the quantile function, and rgt generates random deviates.
Dimitris Rizopoulos [email protected]
## 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)## 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)
Calculates subject-specific predictions for new subjects from a linear mixed model.
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)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)
lmeObject |
an object inheriting from class |
newdata |
a data frame in which to look for variables with which to predict. |
timeVar |
a character string specifying the time variable in the linear mixed model. |
interval |
a character string indicating what type of intervals should be computed. |
all_times |
logical; should predictions be calculated at all |
level |
a numeric scalar denoting the tolerance/confidence level. |
times |
a numeric vector denoting the time points for which we wish to compute the
subject-specific predictions after the last available measurement provided in
|
M |
numeric scalar denoting the number of Monte Carlo samples. See Details. |
return_data |
logical; if |
seed |
numeric scalar, the random seed used to produce the results. |
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
times_to_pred |
time points at which predictions were calculated. |
predicted_y |
the predictions. |
low |
the lower limits of the pointwise confidence intervals. |
upp |
the upper limits of the pointwise confidence intervals. |
Dimitris Rizopoulos [email protected]
## 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)## 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)
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.
| Package: | JMbayes |
| Type: | Package |
| Version: | 0.8-86 |
| Date: | 2020-08-11 |
| License: | GPL (>=2) |
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.
Dimitris Rizopoulos
Maintainer: Dimitris Rizopoulos <[email protected]>
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.
jointModelBayes,
survfitJM,
aucJM,
dynCJM,
prederrJM,
predict.JMbayes,
logLik.JMbayes
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.
The following components must be included in a legitimate JMbayes object.
mcmc |
a list with the MCMC samples for each parameter (except from the random effects if control argument |
postMeans |
a list with posterior means. |
postModes |
a list with posterior modes calculated using kernel desnisty estimation. |
postVarsRE |
a list with the posterior variance-covariance matrix for the random effects of each subject. |
StErr |
a list with posterior standard errors. |
EffectiveSize |
a list with effective sample sizes. |
StDev |
a list with posterior standard deviations. |
CIs |
a list with 95% credible intervals. |
vcov |
the variance-covariance matrix of the model's parameters based. |
pD |
the pD value. |
DIC |
the deviance information criterion value. |
CPO |
the conditional predictive ordinate value. |
LPML |
the log pseudo marginal likelihood value. |
time |
the time used to fit the model. |
scales |
a list with scaling constants in the Metropolis algorithm. |
Covs |
a list with the covariance matrices of the proposals in the Metropolis algorithm. |
acceptRates |
a list of acceptance rates. |
x |
a list with the design matrices for the longitudinal and event processes. |
y |
a list with the response vectors for the longitudinal and event processes. |
Data |
a list of data frames with the data used to fit the models. |
Terms |
a list of terms objects for the various parts of the joint model. |
Funs |
a list of functions used for the various parts of the joint model. |
Forms |
a list of formulas for the two submodels. |
timeVar |
the value of the |
control |
the value of the |
densLongCheck |
a logical indicating whether a scale parameter is required in the longitudinal submodel. |
param |
the value of the |
priors |
a list with the specification of the prior distributions for the model's parameters. This has the same components as
the |
baseHaz |
the value of the |
df.RE |
the value of the |
call |
the matched call. |
Dimitris Rizopoulos [email protected]
Fits shared parameter joint models for longitudinal and survival outcomes under a Bayesian approach.
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(), ...)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(), ...)
lmeObject |
an object of class 'lme' fitted by function |
survObject |
an object of class 'coxph' fitted by function |
timeVar |
a character string indicating the time variable in the mixed effects model. |
param |
a character string specifying the type of association structure between the longitudinal and survival processes. See Details. |
extraForm |
a list with components |
baseHaz |
a character string specifying the type of the baseline hazard function. See Details. |
transFun |
a function or a named list with elements |
densLong |
a function with arguments |
lag |
a numeric scalar denoting a lag effect in the time-dependent covariate represented by the mixed model; default is 0. |
df.RE |
a numeric scalar denoting the number of degrees of freedom for the Student's- |
estimateWeightFun |
logical; experimental, not in use yet. |
weightFun |
a weight function; experimental, not in use yet. |
init |
a named list of user-specified initial values:
When this list of initial values does not contain some of these components or contains components not of the appropriate length, then the default initial values are used instead. |
priors |
a named list of user-specified prior parameters:
|
scales |
a named list with names as in |
control |
a list of control values with components:
|
... |
options passed to the |
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 denote the vector of baseline covariates in
survObject, with associated parameter vector , the subject-specific linear predictor at time point
as defined by the mixed model (i.e., equals the fixed-effects part + random-effects part of the mixed
effects model for sample unit ), 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 and vector of
association parameters for and , respectively. Then, jointModelBayes
assumes a relative risk model of the form
where the baseline risk function is approximated using splines, i.e.,
with denoting a B-spline basis function, a vector of knots, and the corresponding
splines coefficients ( 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 is taken to be proportional to
where denotes the
differences matrix (the order of the differences is set by the control argument diff).
Function specifies the association structure between the two processes. In particular, for
param = "td-value",
for param = "td-extra",
for param = "td-both",
for param = "shared-RE",
and for param = "shared-betasRE",
where and denote possible
transformation functions, denotes the vector of random effects for
the th subject and the fixed effects that correspond to the random effects.
See JMbayesObject for the components of the fit.
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 and the element "AssoctE" that corresponds to the parameter
when parameterization is "td-extra" or "td-both" (see Details).
Dimitris Rizopoulos [email protected]
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.
coef.JMbayes,
ranef.JMbayes,
logLik.JMbayes,
survfitJM,
aucJM,
dynCJM,
prederrJM,
predict.JMbayes
## 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)## 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)
Computes the log-likelihood for a fitted joint model.
## S3 method for class 'JMbayes' logLik(object, thetas, b, priors = TRUE, marginal.b = TRUE, marginal.thetas = FALSE, full.Laplace = FALSE, useModes = TRUE, ...)## S3 method for class 'JMbayes' logLik(object, thetas, b, priors = TRUE, marginal.b = TRUE, marginal.thetas = FALSE, full.Laplace = FALSE, useModes = TRUE, ...)
object |
an object inheriting from class |
thetas |
a list with values for the joint model's parameters. This should have the same structure as
the |
b |
a numeric matrix with random effects value. This should have the same structure as
the |
priors |
logical, if |
marginal.b |
logical, if |
marginal.thetas |
logical, if |
full.Laplace |
logical, if |
useModes |
logical, if |
... |
extra arguments; currently none is used. |
Let denote the vectors of longitudinal responses, the observed event time, and
the event indicator for subject (). Let also denote the probability
density function (pdf) for the linear mixed model, the pdf for the survival submodel, and
the multivariate normal pdf for the random effects, where denotes the full parameter vector. Then,
if priors = TRUE, and marginal.b = TRUE, function logLik() computes
where denotes the prior distribution for the parameters. If priors = FALSE the prior is excluded from the
computation, i.e.,
and when
marginal.b = FALSE, then the conditional on the random effects log-likelihood is computed, i.e.,
when
priors = TRUE and
when priors = FALSE.
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.
Dimitris Rizopoulos [email protected]
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.
## 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)## 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)
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.
marglogLik(object, newdata, idVar = "id", method = "BFGS", control = NULL)marglogLik(object, newdata, idVar = "id", method = "BFGS", control = NULL)
object |
an object inheriting from class |
newdata |
a data frame that contains the longitudinal and covariate information for the subjects for which prediction
of survival probabilities is required. The names of the variables in this data frame must be the same as in the data frames that
were used to fit the linear mixed effects model (using |
idVar |
the name of the variable in |
method |
the |
control |
the |
a numeric vector of marginal log-likelihood contributions.
Dimitris Rizopoulos [email protected]
## 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)## 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)
Fits multivariate mixed models under a Bayesian approach using JAGS.
mvglmer(formulas, data, families, engine = c("JAGS", "STAN"), overdispersion = FALSE, priors = NULL, init = NULL, control = NULL, ...)mvglmer(formulas, data, families, engine = c("JAGS", "STAN"), overdispersion = FALSE, priors = NULL, init = NULL, control = NULL, ...)
formulas |
a list of R formulas representing the mixed models; these should be lme4-type formulas. |
data |
a data.frame that contains all the variable to be used when fitting the multivariate mixed model. |
families |
a list of families objects correspond to each outcome. |
engine |
a character string indicating whether to use JAGS or STAN to fit the model. |
overdispersion |
logical; for Poisson outcomes, should an overdispersion parameter be included. |
priors |
a named list of user-specified prior parameters:
|
init |
a list of initial values. |
control |
a list of control values with components:
|
... |
options passed to the |
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.
A list of class mvglmer with components:
mcmc |
a list with the MCMC samples for each parameter. |
components |
a list with design matrices and responses vectors extracted by applying
the |
data |
a copy of |
control |
a copy of the |
mcmc.info |
a list with information over the MCMC (i.e., time it took, iterations, etc.). |
DIC |
the DIC value for the fitted model. |
pD |
the pD valu for the fitted model. |
Rhat |
a list with the Rhat convergence diagnostics for each parameter. |
priors |
a copy of the priors used. |
postMeans |
a list with posterior means. |
postModes |
a list with posterior modes calculated using kernel desnisty estimation. |
EffectiveSize |
a list with effective sample sizes. |
StErr |
a list with posterior standard errors. |
StDev |
a list with posterior standard deviations. |
CIs |
a list with 95% credible intervals. |
Pvalues |
a list of tail probabilities for containg the zero value. |
call |
the matched call. |
Dimitris Rizopoulos [email protected]
mvJointModelBayes,
jointModelBayes
## 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)## 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)
Fits multivariate shared parameter joint models for longitudinal and survival outcomes under a Bayesian approach.
mvJointModelBayes(mvglmerObject, survObject, timeVar, Formulas = list(NULL), Interactions = list(NULL), transFuns = NULL, priors = NULL, multiState = FALSE, data_MultiState = NULL, idVar_MultiState = "id", control = NULL, ...)mvJointModelBayes(mvglmerObject, survObject, timeVar, Formulas = list(NULL), Interactions = list(NULL), transFuns = NULL, priors = NULL, multiState = FALSE, data_MultiState = NULL, idVar_MultiState = "id", control = NULL, ...)
mvglmerObject |
an object of class 'mvglmer' fitted by function
|
survObject |
an object of class 'coxph' fitted by function |
timeVar |
a character string indicating the time variable in the multivariate mixed effects model. |
Formulas |
a list of lists. Each inner list should have components
|
Interactions |
a list specifying interaction terms for the components of the longitudinal outcomes that are included in the survival submodel. See Examples. |
transFuns |
a character vector providing transformations of the linear predictors of
the mixed models that enter in the linear predictor of the relative risk model.
Currently available options are |
priors |
a named list of user-specified prior parameters:
|
multiState |
logical; if |
data_MultiState |
A data.frame that contains all the variables which were used to fit the multi-state model. This data.frame should be in long format and include one row for each transition for which a subject is at risk. A column called |
idVar_MultiState |
A character string indicating the id variable in |
control |
a list of control values with components:
|
... |
options passed to the |
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.
A list of class mvJMbayes with components:
mcmc |
a list with the MCMC samples for each parameter. |
components |
a copy of the |
Data |
a list of data used to fit the model. |
control |
a copy of the |
mcmc.info |
a list with information over the MCMC (i.e., time it took, iterations, etc.). |
priors |
a copy of the priors used. |
postwMeans |
a list with posterior weighted means. |
postMeans |
a list with posterior means. |
postModes |
a list with posterior modes calculated using kernel desnisty estimation. |
EffectiveSize |
a list with effective sample sizes. |
StErr |
a list with posterior standard errors. |
StDev |
a list with posterior standard deviations. |
CIs |
a list with 95% credible intervals. |
Pvalues |
a list of tail probabilities for containg the zero value. |
call |
the matched call. |
Dimitris Rizopoulos [email protected]
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.
## 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)## 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)
Followup of 312 randomised patients with primary biliary cirrhosis, a rare autoimmune liver disease, at Mayo Clinic.
A data frame with 1945 observations on the following 20 variables.
idpatients identifier; in total there are 312 patients.
yearsnumber of years between registration and the earlier of death, transplantion, or study analysis time.
statusa factor with levels alive, transplanted and dead.
druga factor with levels placebo and D-penicil.
ageat registration in years.
sexa factor with levels male and female.
yearnumber of years between enrollment and this visit date, remaining values on the line of data refer to this visit.
ascitesa factor with levels No and Yes.
hepatomegalya factor with levels No and Yes.
spidersa factor with levels No and Yes.
edemaa 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).
serBilirserum bilirubin in mg/dl.
serCholserum cholesterol in mg/dl.
albuminalbumin in gm/dl.
alkalinealkaline phosphatase in U/liter.
SGOTSGOT in U/ml.
plateletsplatelets per cubic ml / 1000.
prothrombinprothrombin time in seconds.
histologichistologic stage of disease.
status2a numeric vector with the value 1 denoting if the patient was dead, and 0 if the patient was alive or transplanted.
The data frame pbc2.id contains the first measurement for each patient. This data frame is used to
fit the survival model.
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.
Produces MCMC diagnostics plots.
## 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, ...)## 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, ...)
x |
an object inheriting from class |
which |
which types of plots to produce. |
param |
for which parameter to produce the MCMC diagnostic plots; default is for all parameters. |
ask |
logical, if |
max.t |
numeric scalar; up to which time point to plot the weight function, default is up to the third quantile of the observed event times. |
from |
numeric scalar; from which time point to start plotting the weight function, default is zero. |
... |
additional arguments; currently none is used. |
Dimitris Rizopoulos [email protected]
Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.
## 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)## 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)
Produces plots of conditional probabilities of survival.
## 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, ...)## 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, ...)
x |
an object inheriting from class |
estimator |
character string specifying, whether to include in the plot the mean of
the conditional probabilities of survival, the median or both. The mean and median are
taken as estimates of these conditional probabilities over the M replications of the
Monte Carlo scheme described in |
which |
an integer or character vector specifying for which subjects to produce the
plot. If a character vector, then is should contain a subset of the values of the
|
which_subjects |
an integer vector specifying for which subjects to produce the plot. |
split |
a integer vector of length 2 indicating in how many panels to construct, i.e., number of rows and number of columns. |
which_outcomes |
integer vector indicating which longitudinal outcomes to include in the plot. |
surv_in_all |
logical; should the survival function be included in all panels. |
fun |
a vectorized function defining a transformation of the survival curve. For
example, with |
abline |
a list with arguments to |
invlink |
a function to transform the fitted values of the longitudinal outcome. |
conf.int, include_CI
|
logical; if |
fill.area, fill_area_CI
|
logical; if |
col.area, col_fill_CI
|
the color of the area defined by the confidence interval of the survival function. |
col.abline, col.points, col_points, col_lines, col_lines_CI
|
the color for the
vertical line and the points when |
add.last.time.axis.tick |
logical; if |
include.y |
logical; if |
main |
a character string specifying the title in the plot. |
xlab |
a character string specifying the x-axis label in the plot. |
ylab |
a character string specifying the y-axis label in the plot. |
ylab2 |
a character string specifying the y-axis label in the plot,
when |
zlab |
a character string specifying the z-axis (vertical right-hand side) label in the plot. |
lty, lty_lines_CI
|
what types of lines to use. |
col |
which colors to use. |
lwd, lwd_lines
|
the thickness of the lines. |
pch, pch_points
|
the type of points to use. |
ask |
logical; if |
legend |
logical; if |
cex.axis.z, cex.lab.z
|
the par |
cex_xlab, cex_ylab, cex_zlab, cex_main, cex_axis
|
the par |
xlim |
the par |
... |
extra graphical parameters passed to |
Dimitris Rizopoulos [email protected]
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.
## 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)## 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)
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.
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, ...)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, ...)
object |
an object inheriting from class |
newdata |
a data frame that contains the longitudinal and covariate information for the subjects for which prediction
of survival probabilities is required. The names of the variables in this data frame must be the same as in the data frames that
were used to fit the linear mixed effects model (using |
Tstart |
numeric scalar denoting the time point up to which longitudinal information is to be used to derive predictions. |
Thoriz |
numeric scalar denoting the time point for which a prediction of the survival status is of interest; |
lossFun |
either the options |
interval |
logical; if |
idVar |
the name of the variable in |
simulate |
logical; if |
M |
a numeric scalar denoting the number of Monte Carlo samples; see |
... |
additional arguments; currently none is used. |
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:
where denotes the number of subjects at risk at time Tstart, denotes the available
longitudinal measurements up to time , denotes the observed event time for subject , is the event indicator,
is the starting time point Tstart up to which the longitudinal information is used, and is the horizon time point Thoriz.
Function is the loss function that can be the absolute value (i.e., ), the squared value (i.e., ),
or a user-specified function. The probabilities are calculated by survfitJM.
When interval is set to TRUE, then function prederrJM computes the integrated prediction error in the interval
(Tstart, Thoriz) defined as
where
with denoting
the Kaplan-Meier estimator of the censoring time distribution.
A list of class prederrJM with components:
prederr |
a numeric scalar denoting the estimated prediction error. |
nr |
a numeric scalar denoting the number of subjects at risk at time |
Tstart |
a copy of the |
Thoriz |
a copy of the |
interval |
a copy of the |
classObject |
the class of |
nameObject |
the name of |
lossFun |
a copy of the |
Dimitris Rizopoulos [email protected]
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.
survfitJM, aucJM, dynCJM, jointModelBayes
## 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)## 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)
Calculates predicted values for the longitudinal part of a joint model.
## 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, ...)## 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, ...)
object |
an object inheriting from class |
newdata |
a data frame in which to look for variables with which to predict. |
type |
a character string indicating the type of predictions to compute, marginal or subject-specific. See Details. |
interval |
a character string indicating what type of intervals should be computed. |
level |
a numeric scalar denoting the tolerance/confidence level. |
idVar |
a character string indicating the name of the variable in
|
FtTimes |
a list with components numeric vectors denoting the time points
for which we wish to compute subject-specific predictions after the last
available measurement provided in |
last.time |
a numeric vector. This specifies the known time at which each of the subjects in |
LeftTrunc_var |
character string indicating the name of the variable in |
M |
numeric scalar denoting the number of Monte Carlo samples. See Details. |
returnData |
logical; if |
scale |
a numeric value setting the scaling of the covariance matrix of the empirical Bayes estimates in the Metropolis step during the Monte Carlo sampling. |
weight |
a numeric vector of weights to be applied to the predictions of each subject. |
invlink |
a function to tranform the linear predictor of the mixed model to fitted means;
relevant when the user has specified her own density for the longitudinal outcome in
|
seed |
numeric scalar, the random seed used to produce the results. |
... |
additional arguments; currently none is used. |
When type = "Marginal", this function computes predicted values for the
fixed-effects part of the longitudinal submodel. In particular,
let denote the fixed-effects design matrix calculated using
newdata. The predict() calculates ,
and if interval = "confidence", then it calculates the confidence intervals
based on the percentiles of the MCMC sample for .
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" , whereas
for interval = "prediction" is a random vector from a normal
distribution with mean and standard deviation
. Based on this Monte Carlo simulation scheme we take as
estimate of the average of the M estimates
from each Monte Carlo sample. Confidence intervals are constructed using the
percentiles of from the Monte Carlo samples.
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.
The user is responsible to appropriately set the invlink argument when a user-specified
mixed effects model has been fitted.
Dimitris Rizopoulos [email protected]
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.
survfitJM.JMbayes, jointModelBayes
## 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)## 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)
A randomized trial on 488 liver cirrhosis patients
Two data frames with the following variable.
idpatients identifier; in total there are 467 patients.
proprothrobin measurements.
timefor data frame prothro the time points at which the prothrobin measurements were taken;
for data frame prothros the time to death or censoring.
deatha numeric vector with 0 denoting censoring and 1 death.
treatrandomized treatment; a factor with levels "placebo" and "prednisone".
http://www.gllamm.org/books/readme.html#14.6,
Andersen, P. K., Borgan, O., Gill, R. D. and Keiding, N. (1993). Statistical Models Based on Counting Processes. New York: Springer.
Extracts the random effects estimates from a fitted joint model.
## S3 method for class 'JMbayes' ranef(object, postVar = FALSE, ...)## S3 method for class 'JMbayes' ranef(object, postVar = FALSE, ...)
object |
an object inheriting from class |
postVar |
logical; if |
... |
additional arguments; currently none is used. |
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".
Dimitris Rizopoulos [email protected]
Rizopoulos, D. (2012) Joint Models for Longitudinal and Time-to-Event Data: with Applications in R. Boca Raton: Chapman and Hall/CRC.
## 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)## 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)
This function loads the shiny package and runs the application for calculating dynamic predictions using package JMbayes.
runDynPred(type = c("JM", "lme"))runDynPred(type = c("JM", "lme"))
type |
character string indicating whether dynamic predictions are based on joint models or mixed models alone. |
No value returned.
Dimitris Rizopoulos [email protected]
## Not run: runDynPred() ## End(Not run)## Not run: runDynPred() ## End(Not run)
This function computes the conditional probability of surviving later times than the last observed time for which a longitudinal measurement was available.
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, ...)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, ...)
object |
an object inheriting from class |
newdata |
a data frame that contains the longitudinal and covariate information for the subjects for which prediction
of survival probabilities is required. The names of the variables in this data frame must be the same as in the data frames that
were used to fit the linear mixed effects model (using |
type |
character string indicating what to compute, i.e., survival probabilities or the log conditional density. |
idVar |
the name of the variable in |
simulate |
logical; if |
survTimes |
a numeric vector of times for which prediction survival probabilities are to be computed. |
last.time |
a numeric vector or character string. This specifies the known time at which each of the subjects in |
LeftTrunc_var |
character string indicating the name of the variable in |
M |
integer denoting how many Monte Carlo samples to use – see Details. |
CI.levels |
a numeric vector of length two that specifies which quantiles to use for the calculation of confidence interval for the predicted probabilities – see Details. |
log |
logical, should results be returned in the log scale. |
scale |
a numeric scalar that controls the acceptance rate of the Metropolis-Hastings algorithm – see Details. |
weight |
a numeric vector of weights to be applied to the predictions of each subject. |
init.b |
a numeric matrix of initial values for the random effects. These are used in the optimization procedure that finds the mode of the posterior distribution described in Step 2 below. |
seed |
numeric scalar, the random seed used to produce the results. |
... |
additional arguments; currently none is used. |
Based on a fitted joint model (represented by object), and a history of longitudinal responses
and a covariates vector (stored in
newdata), this function provides estimates of , i.e., the conditional probability of surviving time given that subject , with covariate information
, has survived up to time and has provided longitudinal the measurements .
To estimate and if simulate = TRUE, a
Monte Carlo procedure is followed with the following steps:
Take randomly a realization, say from the MCMC sample of posterior of the joint model represented by object.
Simulate random effects values, say , from their posterior distribution given survival up to time ,
the vector of longitudinal responses and . This is achieved using a Metropolis-Hastings algorithm with
independent proposals from a properly centered and scaled multivariate distribution. The scale argument controls the
acceptance rate for this algorithm.
Using and , compute .
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
where denotes the posterior means for the parameters,
and denotes the posterior means for the random effects.
A list of class survfit.JMbayes with components:
summaries |
a list with elements numeric matrices with numeric summaries of the predicted probabilities for each subject. |
survTimes |
a copy of the |
last.time |
a numeric vector with the time of the last available longitudinal measurement of each subject. |
obs.times |
a list with elements numeric vectors denoting the timings of the longitudinal measurements for each subject. |
y |
a list with elements numeric vectors denoting the longitudinal responses for each subject. |
full.results |
a list with elements numeric matrices with predicted probabilities for each subject in each replication of the Monte Carlo scheme described above. |
success.rate |
a numeric vector with the success rates of the Metropolis-Hastings algorithm described above for each subject. |
scale |
a copy of the |
Dimitris Rizopoulos [email protected]
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.
plot.survfit.JMbayes, predict.JMbayes,
aucJM, dynCJM, prederrJM, jointModelBayes
## 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)## 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)
A B-spline expansion of the input variables to be used for a time-varying effect in the specification of joint model.
tve(x, df = NULL, knots = NULL, ord = 3)tve(x, df = NULL, knots = NULL, ord = 3)
x |
a numeric input variable. |
df |
integer denoting the degrees of freedom. |
knots |
a numeric vector of knots. |
ord |
an integer denoting the order of the spline. |
an object of class tve.
Dimitris Rizopoulos [email protected]
produces a LaTeX table with the results of a joint model using package xtable.
## 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, ...)## 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, ...)
x |
an object inheriting from class |
caption |
the |
label |
the |
align |
the |
digits |
the |
display |
the |
which |
a character string indicating which results to include in the LaTeX table. Options are all results, the results of longitudinal submodel or the results of the survival submodel. |
varNames.Long |
a character vector of the variable names for the longitudinal submodel. |
varNames.Event |
a character vector of the variable names for the survival submodel. |
... |
additional arguments; currently none is used. |
A LaTeX code chunk with the results of the joint modeling analysis.
Dimitris Rizopoulos [email protected]
## 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)## 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)