R: Joint Models for Longitudinal and Time-to-Event Data

jointModelBayes {JMbayes}

R Documentation

Joint Models for Longitudinal and Time-to-Event Data

Description

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

Usage

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

Arguments

`lmeObject`	an object of class 'lme' fitted by function `lme()` from package nlme or by function `glmmPQL()` from package MASS.
`survObject`	an object of class 'coxph' fitted by function `coxph()` from package survival.
`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 `fixed` a formula representing the fixed-effects part of the user-defined term, `indFixed` a numeric vector indicating which fixed effects of `lmeObject` 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 `lmeObject` are involved in the user-defined term. Required only when `param = "td-extra"` or `param = "td-both"`. See Examples.
`baseHaz`	a character string specifying the type of the baseline hazard function. See Details.
`transFun`	a function or a named list with elements `value` and `extra` which should be functions. In either case the functions should always have two arguments, namely `x` and `data` (even when the second one is not needed). The purpose is to transform the `value` and/or `extra`, for example including an interaction term, a nonlinear function, etc.
`densLong`	a function with arguments `y`, `eta.y`, `scale`, `log`, and `data` that calculates the density of the longitudinal outcome. `y` denotes the longitudinal responses, `eta.y` the linear predictor that includes the fixed and random effects, `scale` a possible scale parameter (e.g., the measurement error standard deviation), `log` a logical argument that controls whether the density should be calculated in the log scale, and `data` a data frame which may be used to extract variables used in the definition of the density function (e.g., a censoring indicator for left censored longitudinal data).
`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-`t` random-effects distribution. If `NULL` the random effects are assumed to have a multivariate normal distribution.
`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: betas the vector of fixed effects for the linear mixed effects model. tau the precision parameter from the linear mixed effects model (i.e., `\tau = 1/\sigma^2` with `\sigma` denoting the error terms standard deviation). invD the inverse variance-covariance matrix of the random effects. b a matrix of random effects values. gammas the vector of baseline covariates for the survival model. alphas the association parameter(s). Dalphas the association parameter for the true slopes formulation. Bs.gammas the vector of spline coefficients for the spline-approximated baseline risk function. 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: priorMean.betas the prior mean vector of the normal prior for the fixed effects of the linear mixed effects model. priorTau.betas the prior precision matrix of the normal prior for the fixed effects of the linear mixed effects model. priorA.tau the prior shape parameter of the Gamma prior for the precision parameter of the linear mixed effects model. priorB.tau the prior rate parameter of the Gamma prior for the precision parameter of the linear mixed effects model. priorMean.gammas the prior mean vector of the normal prior for the regression coefficients of the survival model. priorTau.gammas the prior precision matrix of the normal prior for the regression coefficients of the survival model. priorMean.alphas the prior mean vector of the normal prior for the association parameter in the survival model. priorTau.alphas the prior precision matrix of the normal prior for the association parameter in the survival model. priorMean.Dalphas the prior mean vector of the normal prior for the slope association parameter in the survival model. priorTau.Dalphas the prior precision matrix of the normal prior for the slope association parameter in the survival model. priorMean.Bs.gammas the prior mean vector of the normal prior for the spline coefficients of the baseline risk function. priorTau.Bs.gammas the prior precision matrix of the normal prior for the spline coefficients of the baseline risk function. priorA.tauBs the prior shape parameter of the Gamma prior for the precision parameter of the penalty term when `baseHaz = "P-splines"`. priorB.tauBs the prior rate parameter of the Gamma prior for the precision parameter of the penal term when `baseHaz = "P-splines"`. priorR.D the prior precision matrix of the Wishart prior for the precision matrix of the random effects. priorK.D the degrees of freedom of the Wishart prior for the precision matrix of the random effects.
`scales`	a named list with names as in `init` specifying scaling constants for the proposal distributions in the Metropolis algorithm.
`control`	a list of control values with components: adapt logical default `FALSE`; should adaptive Metropolis be used. Currently experimental. n.iter integer specifying the total number of iterations; default is 20000. n.burnin integer specifying how many of `n.iter` to discard as burn-in; default is 3000. n.thin integer specifying the thinning of the chains; default is to set `n.thin` such that 2000 samples are kept. n.adapt integer specifying the number of iterations to use for adaptation; default is 3000. keepRE logical; if `TRUE` the MCMC samples for the random effect are kept in the output object. priorVar integer used as the prior precision in the normal prior for the fixed effects, the regression coefficients of the survival submodel, the association parameters, the extra association parameters, and in the spline coefficients; default is 100. 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 for the penalized spline baseline hazard and 5 for the regression spline baseline hazard. ordSpline a numeric scalar 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 a numeric scalar setting the order of the differences in the calculation of the penalty term for the penalized baseline hazard; default is 2. seed a numeric scalar setting the random seed; default is 1. verbose logical; if `TRUE` (default), a progress bar is shown in the console.
`...`	options passed to the `control` argument.

Details

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

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

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

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

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

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

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

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

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

for param = "td-extra",

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

for param = "td-both",

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

for param = "shared-RE",

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

and for param = "shared-betasRE",

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

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

Value

See JMbayesObject for the components of the fit.

Note

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

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

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

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

Author(s)

Dimitris Rizopoulos d.rizopoulos@erasmusmc.nl

References

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

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

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

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

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

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

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

Examples


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

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

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

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

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

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

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

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

summary(jointFit.pbc4)

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

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

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

summary(jointFit.pbc5)


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

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

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

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

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

summary(jointFit.pbc6)

## End(Not run)

[Package JMbayes version 0.8-85 Index]