jmreg.aft: Fit joint model
In sa4khan/AFTjmr: Recurrent event analysis and joint modeling using accelerated failure time (AFT) models

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

Bayesian fit of a joint model. The time-to-event process is described using an AFT model, and the longitudinal process is modeled using lme. It is assumed that the association between the two submodels is induced by the longitudinal value. This version does not allow association induced by the random intercepts and/or slopes. JAGS (Plummer, 2017) must be installed to use this function.

1
2
3

jmreg.aft(surv.fit, lme.fit, surv.model = "weibull", rand.model = "ns",
  timevar, n.chain = 2, n.adapt = 5000, n.burn = 15000,
  n.iter = 1e+05, n.thin = 5, control = list())

`surv.fit`	an object of class `coxph` representing the Cox PH fit. It is required to specify `x=TRUE` in `coxph`.
`lme.fit`	an object of class `lme` (from package `nlme`) representing the linear mixed-effects model fit.
`surv.model`	the AFT model to be used to describe the event process. Available options are `"weibull"`, `"lnormal"`, `"llogistic"`, `"eweibull"` and `"ggamma"` for Weibull, log-normal, log-logistic, exponentiated Weibull and generalized gamma distributions, respectively. The parameterization of these distributions is described in `survreg.aft`.
`rand.model`	a character string indicating the form of the random component in `lme`. Available options are `"ns"` (default) and `"simple"`. Use `"ns"` if the random component is specified using a function of time, expressed in terms of polynomials or splines (i.e., `poly` or `ns`; see Example 2), and use `"simple"` if the random component is specified by the simple linear form b_{0i}+b_{1i}s_{ij}, where s_{ij} is the jth measurement occasion for individual i. For the simple linear form b_{0i}+b_{1i}s_{ij}, `"simple"` is more efficient (`"ns"` also works).
`timevar`	a character string indicating the time variable in the linear mixed effects model.
`n.chain`	number of Markov chains (default is 2).
`n.adapt`	number of iterations for adaptation in JAGS (default is 5000).
`n.burn`	number of iterations to be discarded as burn-in from a chain (default is 15000).
`n.iter`	number of iterations to monitor (default is 100000).
`n.thin`	thinning interval for monitors (default is 5).
`control`	a list of control values (see Details).

The approach of Henderson et al. (2000) has been implemented in JAGS (JAGS must be installed to use this function). The coxph and lme fits are used to organize the data to be used in JAGS. The event process is modeled by an AFT model (available options are Weibull, log-normal, log-logistic, exponentiated Weibull and generalized gamma distributions), and the longitudinal process is characterized by the the linear mixed effects model. The formulation of an AFT model for time-dependent covariates is described in Cox and Oakes (1984). Gauss-Kronrod or Gauss-Legendre quadrature rule is used to approaximate the integral in the definition of the survival function (see control).

Henderson et al. (2000) proposed to model the joint distribution of longitudinal measurements and event time for subject i via a latent zero-mean bivariate Gaussian process \{U_i,V_i\}. This latent process is assumed to generate two linked sub-models, representing the longitudinal and the event time processes.

We model the longitudinal response y_{ij} at time s_{ij} by the relationship y_{ij}=μ_i(s_{ij})+U_i(s_{ij})+ε_{ij}, where μ_i(s_{ij}) is the mean response, U_i(s_{ij}) incorporates subject-specific random effects, and ε_{ij} ~ N(0,σ^2) is a sequence of mutually independent measurement errors. We assume that the mean response at time s is characterized by a linear model μ_i(s)=x_i'(s)α, where x_i(s) is a vector of covariates (possibly time-dependent) and α is the corresponding vector of regression coefficients (fixed effects). For U_i(s), we assume a linear random effects model U_i(s)=w_i'(s)b_i, where w_i(s) is a vector of covariates and b_i ~ N(0,Σ_b) is the corresponding vector of random effects. Note that the length of the fixed-effects coefficient vector α must be > 1. For example, if there is an intercept, there must be at least one covariate/slope. Similarly, the length of the random-effects coefficient vector b_i for individual i must be > 1.

The event intensity process at time t can be expressed as λ_i(t)=λ_0[g_i(t)]exp[-z_i'β-V_i(t)], where g_i(t)=\int_0^texp[-z_i'β-V_i(u)]du, z_i is a vector of baseline covariates, β is the corresponding vector of regression coefficients, and V_i(t) is specified in a similar way to U_i(t).

In our implementation, dependence between the longitudinal and the time-to-event sub-models is captured through V_i(t)= φ U_i(t) so that φ is a measure of association induced by the fitted longitudinal values. The current version does not allow association induced by the random intercepts and/or slopes.

control: a list of control values with components:

mean.alpha: μ_α for prior α ~ N(μ_α, Σ_α), where α is the vector of fixed-effects coefficients of the longitudinal model; default is rep(0, p1), where p1 is the length of α (including an intercept).
inv.cov.alpha: inverse of Σ_α; default is diag(p1)*1e-5, where p1 is the length of α (including an intercept).
mean.beta: μ_β for prior β^* ~ N(μ_β, Σ_β), where β^*=(β,φ), β is the vector of coefficients corresponding to the baseline covariates of the survival model, and φ is the association parameter. default is rep(0, (p2 + 1)), where p2 is the length of β (no intercept).
inv.cov.beta: inverse of Σ_β; default is diag(p2 + 1)*1e-5, where p2 is the length of β (no intercept)
inv.scale.wishart: a p3 x p3 positive definite matrix A for prior Σ_b^{-1} ~ W(A, df), where Σ_b is the covariance matrix of the random effect vector b_i, p3 is the length of b_i, and A is the inverse-scale matrix of the Wishart distribution; default is diag(p3)*0.1.
df.wishart: df for Σ_b^{-1} ~ W(A, df); default is p3 + 1.
inv.sigma2.prior: a vector (shape, rate) for prior σ^{-2} ~ G(shape, rate), where σ^2 is the variance of the longitudinal response and G stands for gamma distribution; default is c(0.1,0.1).
rho.prior: a vector (shape, rate) for prior ρ^{-1} ~ G(shape, rate), where ρ is the inverse-scale parameter of the time-to-event model; default is c(1,0.0005).
kappa.prior: a vector (shape, rate) for prior κ^{2} ~ G(shape, rate), where κ is the shape parameter of the time-to-event model; default is c(0.1,0.1).
gam.prior: a vector (shape, rate) for prior γ^{2} ~ G(shape, rate), where γ is the shape parameter of the time-to-event model (for the exponentiated Weibull and generalized gamma distributions); default is c(0.1,0.1).
init.alpha: initial value for α (a vector of length p1).
init.beta initial value for β^* (a vector of length p2 + 1).
init.inv.sigma2: initial value for σ^{-2} (a scalar).
init.inv.cov.rand: initial value for Σ^{-1} (a p3 x p3 positive denite matrix).
init.rho: initial value for ρ (a scalar).
init.kappa: initial value for κ (a scalar).
init.gam: initial value for γ (a scalar).
init.reffects: initial values for b_i's (an n x p3 matrix).
quad.points: number of quadrature points to approximate the integral involved in the calculation of the survival function; available options are 3, 5, 7, 15 and 21 (default is 5).

survival model: the survival model used to describe the event process.

Association: the association structure linking the two processes.

longitudinal model: LME.

MCMC output: MCMC output, including posterior summaries of the parameters.

time: computational time.

Shahedul Khan <khan@math.usask.ca>

Cox DR and Oakes D, Analysis of survival data, Chapman and Hall/CRC, 1984.

Guo X and Carlin BP, Separate and joint modeling of longitudinal and event time data using standard computer packages, The American Statistician, 58: 16-24, 2004.

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

Plummer M, JAGS Version 4.3.0 user manual, 2017, http://people.stat.sc.edu/hansont/stat740/jags_user_manual.pdf.

jm.DIC, jm.summary, jm.WAIC

# Example 1: AIDS data from package JM
  library(JM)
  data(aids.id)
  data(aids)
  surv.fit<-coxph(Surv(Time,death)~drug+gender+prevOI+AZT,
          data=aids.id,x=TRUE)
  lme.fit<-lme(CD4~obstime+obstime:drug+gender+prevOI+AZT,
         random=~obstime|patient,data=aids)
# weibull time-to-event process
  jmfit.w<-jmreg.aft(surv.fit=surv.fit,lme.fit=lme.fit,
         surv.model="weibull",rand.model="simple",timevar="obstime",
         n.chain=2,n.adapt = 5000, n.burn = 15000, 
         n.iter = 150000, n.thin = 5)
# Computatinoal time: 25.35 minutes in a desktop with intel(R) 
# core(TM)i7-4800, 2.69 GH, 4 cores, 32 GB RAM.
  jm.summary(jmfit.w)
  jm.DIC(jmfit.w)
  jm.WAIC(jmfit.w)
# generalized gamma time-to-event process 
  jmfit.gg<-jmreg.aft(surv.fit=surv.fit,lme.fit=lme.fit,
          surv.model="ggamma",rand.model="simple",timevar="obstime",
          n.chain=2,n.adapt = 5000, n.burn = 15000, 
          n.iter = 150000, n.thin = 5)
# Computatinoal time: 25.35 minutes in a desktop with intel(R) 
# core(TM)i7-6700, 3.40 GH, 4 cores, 40 GB RAM. 
  jm.summary(jmfit.gg)
# 
# Example 2: pbc data
  data(pbc.long)
  data(pbc.surv)
  agec<-pbc.surv$age-mean(pbc.surv$age)
  pbc.surv0<-data.frame(pbc.surv,agec=agec)
# use natural splines
  lme.fit<-lme(log(bilirubin)~drug+ns(futime,2),data=pbc.long,
         random=~ns(futime,2)|id)
  surv.fit<-coxph(Surv(st,status2)~drug*agec,data=pbc.surv0,x=TRUE)
# use rand.model="ns"
  jmfit.ll<-jmreg.aft(surv.fit=surv.fit,lme.fit=lme.fit,
          surv.model="llogistic",rand.model="ns",timevar="futime",
          n.chain=2,n.adapt = 5000, n.burn = 15000, 
          n.iter = 150000, n.thin = 5)
# Computatinoal time: 25.97 minutes in a desktop with intel(R) 
# core(TM)i7-6700, 3.40 GH, 4 cores, 40 GB RAM. 
  jm.summary(jmfit.ll)