R/JMfit.R

In sa4khan/AFTjmr: Recurrent event analysis and joint modeling using accelerated failure time (AFT) models

########################################################################
#' Fit joint model
#' @description Bayesian fit of a joint model. The time-to-event process is described using an AFT model,
#'       and the longitudinal process is modeled using \code{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.
#' @keywords AFT model, JAGS, Joint modeling, MCMC
#' @param surv.fit an object of class \code{coxph} representing the Cox PH fit. 
#'        It is required to specify \code{x=TRUE} in \code{coxph}.
#' @param lme.fit an object of class \code{lme} (from package \code{nlme}) 
#'        representing the linear mixed-effects model fit.
#' @param surv.model the AFT model to be used to describe the event process. Available options are
#'       \code{"weibull"}, \code{"lnormal"}, \code{"llogistic"}, \code{"eweibull"} and \code{"ggamma"} 
#'       for Weibull, log-normal, log-logistic, exponentiated Weibull and generalized gamma
#'       distributions, respectively. The parameterization of these distributions is described in
#'       \code{\link{survreg.aft}}.
#' @param rand.model a character string indicating the form of the random component in \code{lme}.
#'       Available options are \code{"ns"} (default) and \code{"simple"}. Use \code{"ns"} if the
#'       random component is specified using a function of time, 
#'       expressed in terms of polynomials or splines (i.e., \code{poly} or \code{ns}; see Example 2), and use
#'       \code{"simple"} if the random component is specified by the simple linear form 
#'       \eqn{b_{0i}+b_{1i}s_{ij}}, where \eqn{s_{ij}} is the \eqn{j}th measurement occasion
#'       for individual \eqn{i}. For the simple linear form \eqn{b_{0i}+b_{1i}s_{ij}},
#'       \code{"simple"} is more efficient (\code{"ns"} also works).
#' @param timevar a character string indicating the time variable in the linear mixed effects model.
#' @param n.chain number of Markov chains (default is 2).
#' @param n.adapt number of iterations for adaptation in JAGS (default is 5000).
#' @param n.burn number of iterations to be discarded as burn-in from a chain (default is 15000).
#' @param n.iter number of iterations to monitor (default is 100000).
#' @param n.thin thinning interval for monitors (default is 5).
#' @param control a list of control values (see Details).
#' @details The approach of Henderson et al. (2000) has been implemented in JAGS (JAGS must be
#'    installed to use this function). The \code{coxph} and
#'    \code{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 \code{control}).
#'
#'     Henderson et al. (2000) proposed to model the joint distribution of longitudinal measurements and
#'     event time for subject \eqn{i} via a \emph{latent} zero-mean bivariate Gaussian process \eqn{\{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 \eqn{y_{ij}} at time \eqn{s_{ij}} by the relationship
#'     \eqn{y_{ij}=\mu_i(s_{ij})+U_i(s_{ij})+\epsilon_{ij}}, where \eqn{\mu_i(s_{ij})} is the mean response,
#'     \eqn{U_i(s_{ij})} incorporates subject-specific random effects, and \eqn{\epsilon_{ij}} ~ \eqn{N(0,\sigma^2)} 
#'     is a sequence of mutually independent measurement errors. We assume that the mean response at time \eqn{s}
#'     is characterized by a linear model \eqn{\mu_i(s)=x_i'(s)\alpha}, where \eqn{x_i(s)} is a vector of 
#'     covariates (possibly time-dependent) and \eqn{\alpha} is the corresponding vector of
#'     regression coefficients (fixed effects). For \eqn{U_i(s)}, we assume a linear random effects model
#'     \eqn{U_i(s)=w_i'(s)b_i}, where \eqn{w_i(s)} is a vector of covariates and \eqn{b_i} ~ \eqn{N(0,\Sigma_b)} 
#'     is the corresponding vector of random effects. Note that the \code{length} of the fixed-effects coefficient
#'     vector \eqn{\alpha} 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 \eqn{b_i}
#'     for individual \eqn{i} must be > 1. 
#'
#'     The event intensity process at time \eqn{t} can be expressed as 
#'     \eqn{\lambda_i(t)=\lambda_0[g_i(t)]}exp\eqn{[-z_i'\beta-V_i(t)]}, where \eqn{g_i(t)=\int_0^t}exp\eqn{[-z_i'\beta-V_i(u)]du}, 
#'     \eqn{z_i} is a vector of baseline covariates, \eqn{\beta} is the corresponding vector of regression
#'     coefficients, and \eqn{V_i(t)} is specified in a similar way to \eqn{U_i(t)}.
#'
#'     In our implementation, dependence between the longitudinal and the time-to-event sub-models is captured 
#'     through \eqn{V_i(t)= \phi U_i(t)}
#'     so that \eqn{\phi} 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.
#'
#' \code{control}: a list of control values with components:
#'     \itemize{ 
#'     \item \code{mean.alpha}: \eqn{\mu_\alpha} for prior \eqn{\alpha} ~ \eqn{N(\mu_\alpha, \Sigma_\alpha)}, where
#'           \eqn{\alpha} is the vector of fixed-effects coefficients of the longitudinal model;      
#'           default is \code{rep(0, p1)}, where \code{p1} is the \code{length} of \eqn{\alpha} (including an intercept).
#'     \item \code{inv.cov.alpha}: inverse of \eqn{\Sigma_\alpha};      
#'           default is \code{diag(p1)*1e-5}, where \code{p1} is the \code{length} of \eqn{\alpha} (including an intercept).
#'     \item \code{mean.beta}: \eqn{\mu_\beta} for prior \eqn{\beta^*} ~ \eqn{N(\mu_\beta, \Sigma_\beta)}, where
#'           \eqn{\beta^*=(\beta,\phi)}, \eqn{\beta} is the vector of coefficients corresponding to the baseline covariates 
#'           of the survival model, and \eqn{\phi} is the association parameter.      
#'           default is \code{rep(0, (p2 + 1))}, where \code{p2} is the \code{length} of \eqn{\beta} (no intercept).
#'     \item \code{inv.cov.beta}: inverse of \eqn{\Sigma_\beta}; default is \code{diag(p2 + 1)*1e-5}, where
#'           \code{p2} is the \code{length} of \eqn{\beta} (no intercept)
#'     \item \code{inv.scale.wishart}: a \code{p3} x \code{p3} positive definite matrix \eqn{A} for prior \eqn{\Sigma_b^{-1}} ~ \eqn{W(A, df)}, 
#'            where \eqn{\Sigma_b} is the covariance matrix of the random effect vector \eqn{b_i}, \code{p3} is the \code{length} of
#'            \eqn{b_i}, and \eqn{A} is the inverse-scale matrix of the Wishart distribution; default is \code{diag(p3)*0.1}.
#'     \item \code{df.wishart}: df for \eqn{\Sigma_b^{-1}} ~ \eqn{W(A, df)}; default is \code{p3 + 1}.
#'     \item \code{inv.sigma2.prior}: a vector (shape, rate) for prior \eqn{\sigma^{-2}} ~ G(shape, rate),
#'            where \eqn{\sigma^2} is the variance of the longitudinal response and G stands for gamma distribution;
#'            default is \code{c(0.1,0.1)}.
#'     \item \code{rho.prior}: a vector (shape, rate) for prior \eqn{\rho^{-1}} ~ G(shape, rate), where \eqn{\rho} is the
#'            inverse-scale parameter of the time-to-event model; default is \code{c(1,0.0005)}.
#'     \item \code{kappa.prior}: a vector (shape, rate) for prior \eqn{\kappa^{2}} ~ G(shape, rate), where \eqn{\kappa} is the
#'           shape parameter of the time-to-event model; default is \code{c(0.1,0.1)}.
#'     \item \code{gam.prior}: a vector (shape, rate) for prior \eqn{\gamma^{2}} ~ G(shape, rate), where
#'           \eqn{\gamma} is the shape parameter of the time-to-event model (for the exponentiated
#'            Weibull and generalized gamma distributions); default is \code{c(0.1,0.1)}.
#'     \item \code{init.alpha}: initial value for \eqn{\alpha} (a vector of length \code{p1}).
#'     \item \code{init.beta} initial value for \eqn{\beta^*} (a vector of length \code{p2} + 1).
#'     \item \code{init.inv.sigma2}: initial value for \eqn{\sigma^{-2}} (a scalar).
#'     \item \code{init.inv.cov.rand}: initial value for \eqn{\Sigma^{-1}} (a \code{p3} x \code{p3} positive denite matrix).
#'     \item \code{init.rho}: initial value for \eqn{\rho} (a scalar).
#'     \item \code{init.kappa}: initial value for \eqn{\kappa} (a scalar).
#'     \item \code{init.gam}: initial value for \eqn{\gamma} (a scalar).
#'     \item \code{init.reffects}: initial values for \eqn{b_i}'s (an n x \code{p3} matrix).
#'     \item \code{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).
#'          }
#' @return \code{survival model:} the survival model used to describe the event process.
#' @return \code{Association:} the association structure linking the two processes.
#' @return \code{longitudinal model:} LME.
#' @return \code{MCMC output:} MCMC output, including posterior summaries of the parameters.
#' @return \code{time:} computational time.
#' @references Cox DR and Oakes D, Analysis of survival data, Chapman and Hall/CRC, 1984.
#' @references 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.  
#' @references Henderson R, Diggle P, and Dobson A, Joint modelling of longitudinal measurements and event time data, 
#'           Biostatistics, 1: 465-480, 2000.
#' @references Plummer M, JAGS Version 4.3.0 user manual, 2017, 
#'           \url{http://people.stat.sc.edu/hansont/stat740/jags_user_manual.pdf}. 
#' @author Shahedul Khan <khan@math.usask.ca>
#' @seealso \code{\link{jm.DIC}}, \code{\link{jm.summary}}, \code{\link{jm.WAIC}}
#' @examples 
#' # 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)
#' @export

jmreg.aft<-function(surv.fit,lme.fit,
 surv.model="weibull",rand.model="ns",timevar,
 n.chain=2,n.adapt=5000,n.burn=15000,n.iter=100000,n.thin=5,
 control=list()){

# ---------- Checking -----------------------
if(!inherits(lme.fit, "lme")) stop("\nlme.fit must be from class lme.")
if(!inherits(surv.fit, "coxph")) stop("\nsurv.fit must be from class coxph.")
if(is.null(surv.fit$x)) stop("\nUse x = TRUE in coxph fit.")
if(!(surv.model=="weibull" || surv.model=="llogistic" || surv.model=="lnormal" || 
  surv.model=="eweibull" || surv.model=="ggamma"))
stop("\nThe survival model must be one of weibull, llogistic, lnormal, eweibull and ggamma")
if(!exists("timevar"))
stop("\ntimevar is missing")
# ----------------------------------------
if(rand.model=="ns"){
if(surv.model=="ggamma" || surv.model=="eweibull") {
return(bayes.jointM3(surv.fit=surv.fit,lme.fit=lme.fit,timevar=timevar,
 surv.model=surv.model,lme.model=rand.model,
 n.chain=n.chain,n.adapt=n.adapt,n.burn=n.burn,n.iter=n.iter,n.thin=n.thin,control=control))
}
if(surv.model=="llogistic" || surv.model=="weibull" || surv.model=="lnormal") {
return(bayes.jointM2(surv.fit=surv.fit,lme.fit=lme.fit,timevar=timevar,
 surv.model=surv.model,lme.model=rand.model,
 n.chain=n.chain,n.adapt=n.adapt,n.burn=n.burn,n.iter=n.iter,n.thin=n.thin,control=control))
}
}
if(rand.model=="simple"){
if(surv.model=="ggamma" || surv.model=="eweibull") {
return(bayes.jointM3.simple(surv.fit=surv.fit,lme.fit=lme.fit,timevar=timevar,
 surv.model=surv.model,lme.model=rand.model,
 n.chain=n.chain,n.adapt=n.adapt,n.burn=n.burn,n.iter=n.iter,n.thin=n.thin,control=control))
}
if(surv.model=="llogistic" || surv.model=="weibull" || surv.model=="lnormal") {
return(bayes.jointM2.simple(surv.fit=surv.fit,lme.fit=lme.fit,timevar=timevar,
 surv.model=surv.model,lme.model=rand.model,
 n.chain=n.chain,n.adapt=n.adapt,n.burn=n.burn,n.iter=n.iter,n.thin=n.thin,control=control))
}
}
}
#########################################################################