mvJointModelBayes: Multivariate Joint Models for Longitudinal and Time-to-Event...
In JMbayes: Joint Modeling of Longitudinal and Time-to-Event Data under a Bayesian Approach

Description Usage Arguments Details Value Author(s) References See Also Examples

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, ...)

`mvglmerObject`	an object of class 'mvglmer' fitted by function `mvglmer()`.
`survObject`	an object of class 'coxph' fitted by function `coxph()` or 'survreg' fitted by function `survreg()` from package survival.
`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 `fixed` a formula representing the fixed-effects part of the user-defined term, `indFixed` a numeric vector indicating which fixed effects of `mvglmerObject` are involved in the user-defined term, `random` a formula representing the random-effects part of the user-defined term, and `indRamdom` a numeric vector indicating which random effects of `mvglmerObject` are involved in the user-defined term. See Examples.
`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 `"identity"` (identity function), `"expit"` (logistic transformation), `"exp"`, `"log"`, `"log2"`, `"log10"` and `"sqrt"`.
`priors`	a named list of user-specified prior parameters: mean_Bs_gammas the prior mean vector of the normal prior for the B-splines coefficients used to approximate the baseline hazard. Tau_Bs_gammas the prior precision matrix of the normal prior for the B-splines coefficients used to approximate the baseline hazard. mean_gammas the prior mean vector of the normal prior for the regression coefficients of baseline covariates. Tau_gammas the prior precision matrix of the normal prior for the regression coefficients of baseline covariates. mean_alphas the prior mean vector of the normal prior for the association parameters. Tau_alphas the prior mean vector of the normal prior for the association parameters. A_tau_Bs_gammas the prior shape parameter of the Gamma prior for the precision parameter of the penalty term for the B-splines coefficients for the baseline hazard. B_tau_Bs_gammas the prior rate parameter of the Gamma prior for the precision parameter of the penalty term for the B-splines coefficients for the baseline hazard. shrink_gammas logical; should the regression coefficients for the baseline covariates be shrinked. A_tau_gammas the prior shape parameter of the Gamma prior for the precision parameter of the global penalty term baseline regression coefficients. Only relevant when `shrink_gammas = TRUE`. B_tau_gammas the prior rate parameter of the Gamma prior for the precision parameter of the penalty term for the baseline regression coefficients. Only relevant when `shrink_gammas = TRUE`. A_phi_gammas the prior shape parameter of the Gamma prior for the precision parameters of each baseline regression coefficients. Only relevant when `shrink_gammas = TRUE`. B_phi_gammas the prior rate parameter of the Gamma prior for the precision parameters of each baseline regression coefficients. Only relevant when `shrink_gammas = TRUE`. shrink_alphas logical; should the association parameters be shrinked. A_tau_alphas the prior shape parameter of the Gamma prior for the precision parameter of the global penalty term for the association parameters. Only relevant when `shrink_alphas = TRUE`. B_tau_alphas the prior rate parameter of the Gamma prior for the precision parameter of the penalty term for the association parameters. Only relevant when `shrink_alphas = TRUE`. A_phi_alphas the prior shape parameter of the Gamma prior for the precision parameters of each association parameter. Only relevant when `shrink_alphas = TRUE`. B_phi_alphas the prior rate parameter of the Gamma prior for the precision parameters of each association parameter. Only relevant when `shrink_alphas = TRUE`.
`multiState`	logical; if `TRUE` then a joint model for longitudinal and multi-state survival data is fitted.
`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 `trans` indicating the transition to which each row corresponds to, must be included in the data.frame. Function `msprep()` from package mstate can be used to easily convert datasets from wide format (one row per subject) to long format while including the necessary `trans` column. For more information and examples see the documentation for function `msprep()` from package mstate.
`idVar_MultiState`	A character string indicating the id variable in `data_MultiState` for the multi-state model.
`control`	a list of control values with components: n_iter integer specifying the total number of iterations after burn in; default is 300. n_burnin integer specifying how many of iterations to discard as burn-in; default is 1000. n_thin integer specifying the thinning of the chains; default is 300. n_block integer specifying the number of block iterations in which the acceptance rates are checked; default is 100. target_acc a numeric scalar denoting the target acceptance rate; default is 0.234. knots a numeric vector of knots positions for the spline approximation of the log baseline risk function; default is `NULL`, which means that the knots are calculated based on the percentiles of the observed event times. ObsTimes.knots logical; if `TRUE` (default), the knots are set using the percentiles of the observed event times (i.e., including both true events and censored observations). If `FALSE`, the knots are set based on the percentiles of the true event times alone. lng.in.kn a numeric scalar indicating the number of knots to use (based on the percentiles); default is 15. ordSpline an integer setting the order of the spline function. This is the number of coefficients in each piecewise polynomial segment, thus a cubic spline has order 4; default is 4. diff an integer setting the order of the differences in the calculation of the penalty term for the penalized baseline hazard; default is 2. seed an integer setting the random seed; default is 1. n_cores an integer specifying the number of cores to use. Default is the the available cores minus one. GQsurv a character string specifying the type of Gaussian quadrature to be used. Options are "GaussKronrod" (default) and "GaussLegendre". GQsurv.k an integer specifying the number of quadrature points; default is 15.
`...`	options passed to the `control` argument.

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 `components` element of `mvglmerObject`.
`Data`	a list of data used to fit the model.
`control`	a copy of the `control` values used in the fit.
`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 d.rizopoulos@erasmusmc.nl

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.

mvglmer, jointModelBayes

## 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)