The R Package JMbayes for Fitting Joint Models for Longitudinal and Time-to-Event Data using MCMC

Joint models for longitudinal and time-to-event data constitute an attractive modeling framework that has received a lot of interest in the recent years. This paper presents the capabilities of the R package JMbayes for fitting these models under a Bayesian approach using Markon chain Monte Carlo algorithms. JMbayes can fit a wide range of joint models, including among others joint models for continuous and categorical longitudinal responses, and provides several options for modeling the association structure between the two outcomes. In addition, this package can be used to derive dynamic predictions for both outcomes, and offers several tools to validate these predictions in terms of discrimination and calibration. All these features are illustrated using a real data example on patients with primary biliary cirrhosis.


Introduction
Joint models for longitudinal and time-to-event data constitute an attractive modeling paradigm that currently enjoys great interest in the statistics and medical literature (Rizopoulos and Lesaffre 2014;Rizopoulos 2012;Tsiatis and Davidian 2004). These models are utilized in follow-up studies where interest is in associating a longitudinal response with an event time outcome. In general, there are mainly two settings in which such type of models are required. First, when one is interested in measuring the strength of the association between the hazard of an event and a time-varying covariate, then we should pay special attention to the attributes of the covariate process. In particular, when this is an endogenous time-varying covariate (Kalbfleisch and Prentice 2002, Section 6.3), standard methods, such as the timedependent Cox model (Therneau and Grambsch 2000), are not optimal for measuring this association. Standard examples of endogenous covariates are covariates, which are measured on the sample units themselves, for instance, biomarkers or other parameters measured on patients during follow-up. The important feature of such covariates is that their existence and/or future path is directly related to the event status. By postulating a model for the joint distribution of the covariate and event processes we explicitly acknowledge this link, and hence we obtain a more accurate estimate for their association. The second case in which joint models are of use is when one needs to account for incomplete data. More specifically, when the probability of missingness depends on unobserved longitudinal responses, then in order to obtain valid inferences we need to postulate a model for the joint distribution of the arXiv:1404.7625v1 [stat.CO] 30 Apr 2014 longitudinal and missingness processes (Little and Rubin 2002;Molenberghs and Kenward 2007). In this context, three main frameworks have been proposed to define such joint distributions, namely, selection, pattern mixture and shared parameter models. The majority of the models that have been proposed in the literature under these frameworks have focused on standard designs assuming a fixed set of time points at which subjects are expected to provide measurements. Nonetheless, in reality, subjects often do not adhere to the posited study schedule and may skip visits and dropout from the study at random time points. Even though in many of those occasions information on the exact dropout time is available, the typical convention in selection and pattern mixture modeling has been to ignore this feature and coerce measurements back to discrete follow-up times. Another alternative that makes better use of the available data is to acknowledge that dropout occurs in continuous time and consider it as a time-to-event outcome.
Following the increasing in these models, currently there are several software implementations available to fit them. The R package JM Rizopoulos (2014Rizopoulos ( , 2012Rizopoulos ( , 2010 fits joint models for a continuous longitudinal outcome and an event time process under maximum likelihood. Several types of association structures are supported and the package also allows to fit joint models with competing risk survival data. In addition, JM can be used to calculate dynamic predictions for either of the two outcomes. The R package joineR (Philipson, Sousa, Diggle, Williamson, Kolamunnage-Dona, and Henderson 2012) similarly fits joint models for a continuous longitudinal outcome and a time-to-event, following the formulation of Henderson, Diggle, and Dobson (2000). In addition, the stjm package for STATA (Crowther 2013) implements joint modeling of a normal longitudinal response and a time-to-event using maximum likelihood, with emphasis on parametric time-to-event models. The implementation of joint models in SAS and WinBUGS has been discussed by Guo and Carlin (2004). Finally, contrary to the previous software implementations, function Jointlcmm() from the R package lcmm (Proust-Lima, Philipps, Diakite, and Liquet 2013) fits joint latent class models for a continuous longitudinal outcome and a survival outcome using maximum likelihood; these models postulate that the association between the two processes is captured by categorical random effects (i.e., latent classes).
In this paper we introduce the R package JMbayes that fits joint models under a Bayesian approach. JMbayes can fit a wide range of joint models, including among others joint models for continuous and categorical longitudinal responses. It provides several options for modeling the association structure between the two outcomes, with the possibility of different terms from the longitudinal submodel entering the linear predictor of the survival submodel, allowing also for general transformation of these terms. In addition, the package provides extensive capabilities to derive dynamic predictions for both outcomes, it allows to combine predictions from different models using innovative Bayesian model averaging techniques, and facilitates the utilization of these predictions in practice using a web interface. Moreover, it offers several tools to quantify the quality of these predictions in terms of discrimination and calibration, and code is provided for their validation. The rest of the paper is organised as follows. Section 2 presents a short review of the underlying methodological framework behind joint models. Section 3 gives the details behind the implementation of joint models in package JMbayes, and Section 4 illustrates in detail the use of the package in a real dataset on patients with primary biliary cirrhosis. Finally, Section 5 presents how dynamic predictions for the longitudinal and event time outcomes are defined, and how they can be calculated and validated with the package.

Theoretical framework
Let D n = {T i , δ i , y i ; i = 1, . . . , n} denote a sample from the target population, where T * i denotes the true event time for the i-th subject , C i the censoring time, T i = min(T * i , C i ) the corresponding observed event time, and δ i = I(T * i ≤ C i ) the event indicator, with I(·) being the indicator function that takes the value 1 when T * i ≤ C i , and 0 otherwise. In addition, we let y i denote the n i × 1 longitudinal response vector for the i-th subject, with element y il denoting the value of the longitudinal outcome taken at time point t il , l = 1, . . . , n i .
To accommodate different types of longitudinal responses in a unified framework, we postulate a generalized linear mixed effects model. In particular, the conditional distribution of y i given a vector of random effects b i is assumed to be a member of the exponential family, with linear predictor given by g E{y where g(·) denotes a known one-to-one monotonic link function, and y i (t) denotes the value of the longitudinal outcome for the i-th subject at time point t, x i (t) and z i (t) denote the timedependent design vectors for the fixed-effects β and for the random effects b i , respectively. The random effects are assumed to follow a multivariate normal distribution with mean zero and variance-covariance matrix D. For the survival process, we assume that the risk for an event depends on a function of the subject-specific linear predictor η i (t). More specifically, we have where H i (t) = {η i (s), 0 ≤ s < t} denotes the history of the underlying longitudinal process up to t, h 0 (·) denotes the baseline hazard function, w i is a vector of baseline covariates with corresponding regression coefficients γ. Parameter vector α quantifies the association between features of the marker process up to time t and the hazard for an event at the same time point. Various options for the form of function f (·) are presented in Section 4.3. To complete the specification of the survival process we need to make appropriate assumptions for the baseline hazard function h 0 (·). To model this function, while still allowing for flexibility, we use a B-splines approach. In particular, the logarithm of the baseline hazard function is expressed as where B q (t, v) denotes the q-th basis function of a B-spline with knots v 1 , . . . , v Q and γ h 0 the vector of spline coefficients. Increasing the number of knots Q increases the flexibility in approximating log h 0 (·); however, we should balance bias and variance and avoid over-fitting. A standard rule of thumb is to keep the total number of parameters, including the parameters in the linear predictor in (2) and in the model for h 0 (·), between 1/10 and 1/20 of the total number of events in the sample (Harrell 2001, Section 4.4). After the number of knots has been decided, their location can be based on percentiles of the observed event times T i or of the true event times {T i : T * i ≤ C i , i = 1, . . . , n} in order to allow for more flexibility in the region of greatest density. A standard alternative approach that avoids the task of choosing the appropriate number and position of the knots is to include a relatively high number of knots (e.g., 15 to 20) and appropriately penalize the B-spline regression coefficients γ h 0 for smoothness (Eilers and Marx 1996).
Under the Bayesian approach, estimation of joint model's parameters proceeds using Markov chain Monte Carlo (MCMC) algorithms. The expression for the posterior distribution of the model parameters is derived under the assumptions that given the random effects, both the longitudinal and event time process are assumed independent, and the longitudinal responses of each subject are assumed independent. Formally we have, where θ denotes the full parameter vector, and p(·) denotes an appropriate probability density function. Under these assumptions the posterior distribution is analogous to: where with ψ il (b i ) and ϕ denoting the natural and dispersion parameters in the exponential family, respectively, c(·), a(·), and d(·) are known functions specifying the member of the exponential family, and for the survival part with h i (·) given by (2). The integral in the definition of the survival function does not have a closed-form solution, and thus a numerical method must be employed for its evaluation. Standard options are the Gauss-Kronrod and Gauss-Legendre quadrature rules.
For the parameters θ we take standard prior distributions. In particular, for the vector of fixed effects of the longitudinal submodel β, for the regression parameters of the survival model γ, for the vector of spline coefficients for the baseline hazard γ h 0 , and for the association parameter α we use independent univariate diffuse normal priors. The penalized version of the B-spline approximation to the baseline hazard can be fitted by specifying for γ h 0 the improper prior (Lang and Brezger 2004): where τ h is the smoothing parameter that takes a Gamma(1, 0.005) hyper-prior in order to ensure a proper posterior for γ h 0 , K = ∆ r ∆ r , where ∆ r denotes r-th difference penalty matrix, and ρ(K) denotes the rank of K. For the covariance matrix of the random effects we assume an inverse Wishart prior, and when fitting a joint model with a normally distributed longitudinal outcome, we take an inverse-Gamma prior for the variance of the error terms σ 2 . More details regarding Bayesian estimation of joint models can be found in Ibrahim, Chen, and Sinha (2001, Chapter 7) and Brown, Ibrahim, and DeGruttola (2005).
3. The R package JMbayes

Design
In many regards the design of package JMbayes is similar to the one of package JM for fitting joint models under maximum likelihood. In particular, JMbayes has a basic modelfitting function called jointModelBayes(), which accepts as main arguments a linear mixed effects object fit as returned by functions lme() of package nlme (Pinheiro, Bates, De-bRoy, Sarkar, and R Development Core Team 2014) or from function glmmPQL() from package MASS (Venables and Ripley 2002), and a survival object fit as returned by function coxph() of package survival (Therneau and Lumley 2014). The final required argument is timeVar, a character string that denotes the name of the time variable in the mixed model. By default jointModelBayes() fits joint models with a linear mixed effects submodel for a continuous longitudinal outcome, and a relative risk submodel of the form (2 i.e., the risk for an event at time t is associated with the subjectspecific mean of the longitudinal outcome at the same time point. Joint models for other types of longitudinal outcomes can be fitted by appropriately specifying argument densLong, and arguments param, extraForm and transFun can be used to add extra terms involving components of the longitudinal process and possibly transform these terms. A detailed account on the use of these arguments, with examples, is given in Section 4. The baseline hazard is by default approximated using penalized B-splines; regression splines can be instead invoked by appropriately setting argument baseHaz. The number and position of the knots can be controlled via the lng.in.kn and knots control arguments. The former defines the number of knots to use (by default placed at equally spaced percentiles of the observed event times), whereas argument knots can be invoked to specify knots at specific positions. The type of numerical integration algorithm used to approximate the survival function (7) is specified with the control argument GQsurv with options "GaussKronrod" (default) and "GaussLegendre", while the number of quadrature points is specified using the control argument GQsurv.k (for the Gauss-Kronrod rule only 7 or 15 can be specified). The fitting process can be further appropriately tweaked using a series of extra control arguments explained in the following section and in the help file of jointModelBayes(). In addition, the default values of the parameters of the prior distributions can be altered using the priors argument, and analogously the default initial values using argument init.

Implementation details
The MCMC algorithm that samples from the posterior conditional distributions of the parameters and the random effects is implemented by the internal function MCMCfit(). For the majority of the posterior conditionals random walk Metropolis is used, with exceptions for the precision parameter of the error terms distribution when a linear mixed model is used for the longitudinal outcome in which case slice sampling is used, and for the random effects precision matrix D −1 in which case when the random effects are assumed normally distributed the posterior conditional is a Wishart distribution (if argument df.RE of jointModelBayes() is not NULL the distribution of the random effects is assumed to be a Student's-t distribution with df.RE degrees of freedom; in this case the random effects precision matrix is updated with a Metropolis-Hastings algorithm). The implementation behind MCMCfit() takes full advantage of the separately fitted mixed effects and Cox models in order to appropriately define the covariance matrix of the normal proposal distributions for the random walk Metropolis algorithm. In particular, for β and b i these covariance matrices are taken from the mixed model, whereas for the regression coefficients in the linear predictor of the survival submodel and the B-spline coefficients γ h 0 a two-stage approach is employed, where a time-dependent Cox model is fitted using the mixed model to compute f {η i (t), b i , α}. These proposal distributions are tuned during an adaptive phase of n.adapt iterations (default 3000), where every n.batch iterations (default 100) the acceptance rate of the algorithms are checked. Following a burn-in period of n.burnin iterations (default 3000) is performed, and after these iterations the algorithm continues to run for an extra of n.iter iterations (default 20000). The chains are thinned according to the n.thin argument (default is to keep 2000 iterations for each parameter).
From the two schools of running MCMC algorithms, namely the 'one long chain school' and the 'multiple shorter chains school', JMbayes implements the former. Users who wish to check convergence using multiple chains can still do it but with a bit of extra programming. More specifically, they could call jointModelBayes() with different initial values (by appropriately specifying argument init), and following they could extract component mcmc from the fitted models, which is the list of simulated values for each parameter. These lists could subsequently be processed using the coda package (Plummer, Best, Cowles, and Vines 2006) and perform these diagnostic tests.

The basic joint model
We will illustrate the capabilities of package JMbayes using the primary biliary cirrhosis (PBC) data collected by the Mayo Clinic from 1974 to 1984 (Murtaugh, Dickson, Van Dam, Malincho, Grambsch, Langworthy, and Gips 1994). PBC is a chronic, fatal, but rare liver disease characterized by inflammatory destruction of the small bile ducts within the liver, which eventually leads to cirrhosis of the liver. Patients with PBC have abnormalities in several blood tests, such as elevated levels of serum bilirubin. For our analysis we will consider 312 patients who have been randomized to D-penicillamine and 154 placebo. During follow-up several biomarkers associated with PBC have been collected for these patients. Here we focus on serum bilirubin levels, which is considered one of the most important ones associated with disease progression. Patients had on average 6.2 measurements (std. deviation 3.8 measurements), with a total of 1945 observations. In package JMbayes the PBC data are available in the data frames pbc2 and pbc2.id containing the longitudinal and survival information, respectively (i.e., the former is in the long format while the latter contains a single row per patient).
We start by loading packages JMbayes and lattice (Sarkar 2008) and defining the indicator status2 for the composite event, namely transplantation or death: R> library("JMbayes") R> library("lattice") R> pbc2$status2 <-as.numeric(pbc2$status != "alive") R> pbc2.id$status2 <-as.numeric(pbc2.id$status != "alive") Descriptive plots for the survival and longitudinal outcomes are presented in Figures 1 and 2 that depict the Kaplan-Meier estimate of transplantation-free survival for the two treatment groups, and the sample subject-specific longitudinal trajectories for patients with and without and endpoint, respectively.
R> sfit <-survfit(Surv(years, status2)~drug, data = pbc2.id) R> plot(sfit, lty = 1:2, lwd = 2, col = 1:2, mark.time = FALSE, + xlab = "Time (years)", ylab = "Transplantation-free Survival") R> legend("topright", levels(pbc2.id$drug), lty = 1:2, col = 1:2, lwd = 2, + cex = 1.3, bty = "n") R> pbc2$status2f <-factor(pbc2$status2, levels = 0:1, + labels = c("alive", "transplanted/dead")) R> xyplot(log(serBilir)~year | status2f, group = id, data = pbc2, type = "l", + col = 1, xlab = "Time (years)", ylab = "log(serum Bilirubin)") We continue by separately fitting a linear mixed model for the longitudinal and a Cox model for the survival one. Careful investigation of the shapes of the log serum bilirubin profiles indicates that for many individuals these seem to be nonlinear. Hence, to allow for flexibility in the specification of these profiles we include natural cubic splines in both the fixed-and random-effects parts of the mixed model. This model can be fitted using the following call to functions lme() and ns() (the latter from package splines): R> lmeFit.pbc1 <-lme(log(serBilir)~ns(year, 2), data = pbc2, + random =~ns(year, 2) | id) Analogously, in the Cox model we control for treatment and age, and also allow for their interaction: R> coxFit.pbc1 <-coxph(Surv(years, status2)~drug * age, data = pbc2.id, In the call to coxph() argument x is set to TRUE such that the design matrix is also included in the resulting model object. Using as main arguments the lmeFit.pbc1 and coxFit.pbc1 objects, the corresponding joint model is fitted using the code: R> jointFit.pbc1 <-jointModelBayes(lmeFit.pbc1, coxFit.pbc1, timeVar = "year", + n.iter = 30000) R> summary(jointFit.pbc1) Call: jointModelBayes(lmeObject = lmeFit.pbc1, survObject = coxFit.pbc1,   As explained earlier, argument timeVar is a character string that specifies the name of the time variable in the mixed model (the scale of time (e.g., days, months, years) in both the mixed and Cox models must be the same). In addition, using the control argument n.iter we specified that after adaption and burn-in, the MCMC should run for 30000 iterations. The default call to jointModelBayes() includes in the linear predictor of the relative risk model the subject-specific linear predictor of the mixed model η i (t), which in this case represents the average subject-specific log serum bilirubin level. The output of the summary() method is rather self-explanatory and contains model summary statistics, namely LPML (the log pseudo marginal likelihood value), DIC (deviance information criterion), and pD (the effective number of parameters component of DIC), posterior means for all parameters, and standard errors (effective sample size estimated using time series methodology), standard deviations, 95% credibility intervals and tail probabilities for all regression coefficients in the two submodels. The association parameter α is denoted in the output as Assoct. The tail probabilities, under the column with the heading P, are calculated as 2 × min{Pr(θ > 0), Pr(θ < 0)}, with θ denoting here the corresponding regression coefficient from the longitudinal or the survival submodel. The results suggest that serum bilirubin is strongly associated with the risk for the composite event, with a doubling of serum bilirubin levels, resulting in a 2.7-fold (95% CI: 2.4; 3.1) increase of the risk. In the appendix we show how the plot() method can be used produce diagnostic plots for investigating the convergence of the MCMC.

Extended joint models
The previous section showed how the basic joint model for a continuous normally distributed longitudinal outcome and a time-to-event can be fitted in JMbayes. In this section we will illustrate how joint models for other types of longitudinal responses may be fitted using function jointModelBayes() by suitably specifying argument densLong. In particular, this argument accepts a function that calculates the probability density function (and its natural logarithm) of the longitudinal outcome, with arguments y denoting the vector of longitudinal responses y, eta.y the subject-specific linear predictor η i (t), scale a potential scale parameter (e.g., the standard deviation of the error terms), log a logical denoting whether logarithm of the density is computed, and data a data frame that contains variables that are potentially required in the definition of densLong. To better illustrate the use of this function, we present here three examples of joint models with more elaborate longitudinal outcomes. We start with an extension of model jointFit.pbc1 that allows for a more heavier-tailed error distribution, that is, Function dgt() of package JMbayes calculates the probability density function of the generalized Student's t distribution (i.e., a Student's t with mean parameter mu and scale parameter sigma). Supplying this function in the densLong fits the corresponding joint model: R> jointFit.pbc2 <-jointModelBayes(lmeFit.pbc1, coxFit.pbc1, timeVar = "year", + densLong = dLongST) R> summary(jointFit.pbc2) . . . We observe some slight changes in the regression coefficients of both submodels, where a doubling of serum bilirubin levels is now associated with a 2.6-fold (95% CI: 2.3; 2.9) increase of the risk for the composite event.
Following we illustrate the use of densLong for fitting a joint model with a dichotomous (binary) longitudinal outcome. Since in the PBC data there was no binary biomarker recorded during follow-up, we artificially create one by dichotomizing serum bilirubin at the threshold value of 1.8 mg/dL. To fit the corresponding joint model we need first to fit a mixed effects logistic regression for the longitudinal binary outcome using function glmmPQL() from package MASS, the syntax is R> pbc2$serBilirD <-as.numeric(pbc2$serBilir > 1.8) R> lmeFit.pbc2 <-glmmPQL(serBilirD~year, random =~year | id, + family = binomial, data = pbc2) As for continuous longitudinal outcomes, this mixed effects model object is merely used to extract the required data (response vector, design matrices for fixed and random effects), and starting values for the parameters and random effects. The definition of densLong and the call to jointModelBayes() take the form: . . .
As we have already seen, the default parameterization posits that the subject-specific linear predictor η i (t) from the mixed model is included as a time-varying covariate in the relative risk model. This means that, in this case, the estimate of the association parameter α = 0.2 denotes the log hazard ratio for a unit increase in the log odds of having serum bilirubin above 1.8 mg/dL. The flexibility that the user has in defining her own density function for the longitudinal outcome is evident, for example, we can easily fit a mixed effects probit regression instead of using the logit link by defining densLong as R> dLongBin <-function (y, eta.y, scale, log = FALSE, data) { + dbinom(x = y, size = 1, prob = pnorm(eta.y), log = log) + } As a final example, we illustrate how densLong can be utilized to fit joint models with censored longitudinal data (detection limit problem) by making use of extra variables in the data frame containing the longitudinal information. Similarly to the previous example, the biomarkers collected in PBC study were not subject to detection limits, and therefore we again artificially create a censored version of serum bilirubin with values below the threshold value of 0.8 mg/dL set equal to the detection limit of 0.8 mg/dL. The code creating the censored longitudinal response vector is: In addition to the censored version of serum bilirubin we have also included in the data frame pbc2 the censoring indicator CensInd. We again assume a normal error distribution for the logarithm of serum bilirubin but in the definition of the corresponding density function we need to account for censoring, that is for observations above the detection limit we use the density function whereas for observations below this limit we use the cumulative distribution function. The definition of the censored density becomes: R> 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 + log.dens <-(1 -ind) * log.f + ind * log.F + if (log) log.dens else exp(log.dens) + } Note that the censoring indicator is extracted from the data argument of censdLong(). Again in order to fit the joint model, we first need to fit the linear mixed model for the censored response variable serBilir2 and following supply this object and the censdLong() to jointModelBayes(), i.e., R> lmeFit.pbc3 <-lme(log(serBilir2)~ns(year, 2), data = pbc2, + random =~ns(year, 2) | id) R> jointFit.pbc4 <-jointModelBayes(lmeFit.pbc3, coxFit.pbc1, timeVar = "year", + densLong = censdLong) R> summary(jointFit.pbc4) . . . We observe that the estimate of the association parameter α is relatively close to the estimate obtained in model jointFit.pbc1 that was based on the original (uncensored) version of serum bilirubin.

Association structures
The joint models we fitted in Sections 4.1 and 4.2 assumed that the hazard for an event at any time t is associated with the current underlying value of the biomarker at the same time point, denoted as η i (t), and the strength of this association is measured by parameter α. Even though under this formulation parameter α enjoys a clear interpretation, it is not realistic to expect that it will always be the most appropriate in expressing the correct relationship between the two processes. In general, there could be other characteristics of the subjects' longitudinal profiles that are more strongly predictive for the risk of an event; for example, the rate of increase/decrease of the biomarker's levels or a suitable summary of the whole longitudinal trajectory, among others. In this section we illustrate how such association structures could be postulated and fitted with jointModelBayes().
We start with the parameterization proposed by Ye, Lin, and Taylor (2008), Brown (2009) and Rizopoulos (2012) that posits that the risk depends on both the current true value of the trajectory and its slope at time t. More specifically, the relative risk survival submodel takes the form, where The interpretation of parameter α 1 remains the same as in the standard parameterization. Parameter α 2 measures the association between the slope of the true longitudinal trajectory at time t and the risk for an event at the same time point, provided that η i (t) remains constant. To fit the joint model with the extra slope term in the relative risk submodel we need to specify the param and extraForm arguments of jointModelBayes(). The first one is a character string with options "td-value" (default) that denotes that only the current value term η i (t) is included, "td-extra" which means that only the extra, user-defined, term is included, and "td-both" which means that both η i (t) and the user-defined terms are included. The exact definition of the extra term is provided via the argument extraForm which is a list with four components, namely * "fixed" an R formula specifying the fixed-effects part of the extra term, * "random" an R formula specifying the random-effects part of the extra term, * "indFixed" an integer vector denoting which of the fixed effects of the original mixed model are encountered in the definition of the extra term, and * "indRandom" an integer vector denoting which of the random effects of the original mixed model are encountered in the definition of the extra term.
For example, to include the slope term η i (t) under the linear mixed model lmeFit.pbc1, this list takes the form: R> dForm <-list(fixed =~0 + dns(year, 2), random =~0 + dns(year, 2), + indFixed = 2:3, indRandom = 2:3) Function dns() computes numerically (with a central difference approximation) the derivative of a natural cubic spline as calculated by function ns() (there is also a similar function dbs() that computes numerically the derivative of a cubic spline as calculated by function bs()). The corresponding joint model is fitted with the code: We observe that both the current level and the current slope of the longitudinal profile are strongly associated with the risk for the composite event. For patients with the same treatment and age at baseline, and who have the same underlying level of serum bilirubin at time t, if serum bilirubin has increased by 50% within a year then the corresponding hazard ratio is 2.9 (95% CI: 1.8; 4.5).
A common characteristic of the two parameterizations we have seen so far is that the risk for an event at any time t is assumed to be associated with features of the longitudinal trajectory at the same time point (i.e., current value η i (t) and current slope η i (t)). However, this assumption may not always be appropriate, and we may benefit from allowing the risk to depend on a more elaborate function of the history of the time-varying covariate (Sylvestre and Abrahamowicz 2009). In the context of joint models, one option to account for the cumulative effect of the longitudinal outcome is to include in the linear predictor of the relative risk submodel the integral of the longitudinal trajectory from baseline up to time t (Brown 2009;Rizopoulos 2012). More specifically, the survival submodel takes the form where for any particular time point t, α measures the strength of the association between the risk for an event at time point t and the area under the longitudinal trajectory up to the same time t, with the area under the longitudinal trajectory taken as a summary of the whole marker history H i (t) = {m i (s), 0 ≤ s < t}. To fit a joint model with this term in the linear predictor of the survival submodel, we need first again to appropriately define the formulas that calculate its fixed-effects and random-effects parts. Similarly to including the slope term, the list with these formulas takes the form R> iForm <-list(fixed =~0 + year + ins(year, 2), + random =~0 + year + ins(year, 2), + indFixed = 1:3, indRandom = 1:3) where function ins() calculates numerically (using the Gauss-Kronrod rule) the integral of function ns(). The corresponding joint model is fitted by supplying this list in the extraForm argument and by also setting in argument param that we only want to include the extra term in the linear predictor of the survival submodel: To explicitly denote that in the relative risk submodel we only want to include the user-defined integral term, we have set argument param to "td-extra". Similarly to the previous results we observe that the area under the longitudinal profile of log serum bilirubin is strongly associated with the risk for an event, with a unit increase corresponding to a 1.3-fold (95% CI: 1.2; 1.3) increase of the risk.
The final type of association structure we consider assumes that only the random effects are shared between the two processes, namely or potentially the corresponding fixed effects may also be included, i.e., with β b denoting the fixed effects that correspond to the random effects. This type of parameterization is more meaningful when a simple random-intercepts and random-slopes structure is assumed for the longitudinal submodel, in which case the random effects express subjectspecific deviations from the average intercept and average slope. Under this setting this parameterization postulates that patients who have a lower/higher level for the longitudinal outcome at baseline (i.e., intercept) or who show a steeper increase/decrease in their longitudinal trajectories (i.e., slope) are more likely to experience the event. In that respect, this formulation shares also similarities with the time-dependent slopes formulation (8). A joint model with a relative risk model of the form (9) can be fitted by setting argument param to "shared-RE" in the call to jointModelBayes(), whereas formulation (10)  The results suggest that both the baseline levels of the underlying log serum bilirubin (i.e., parameter Assoct:(Intercept)) as well as the longitudinal evolution of the marker (i.e., parameters Assoct:ns(year, 2)1 and Assoct:ns(year, 2)2) are strongly related to the hazard of the composite event.

Transformation functions
The previous section illustrated several options for the definition of function f (·) in (2) for studying which features of the longitudinal process are associated with the event of interest. Yet another set of options for function f (·) would be to consider adding interaction or nonlinear terms for the components of the longitudinal outcome that are included in the linear predictor of the relative risk model. Such options are provided in jointModelBayes() by suitably specifying argument transFun. This should be a function (or a list of two functions) with arguments x denoting the term from the longitudinal model, and data a data frame that contains other variables that potentially should be included in the calculation. When a single function is provided, then this function is applied to the current value term η i (t) and potentially also to the extra term provided by the user if param was set to "td-both". If a list is provided, then this should be a named list with components "value" and "extra" providing separate functions for the current value and the user-defined terms, respectively. We illustrate how these transformation functions can be used in practice by extending model jointFit.pbc12, which included the current value term η i (t) and the current slope term η i (t), by including the quadratic effect of η i (t) and the interaction of η i (t) with the randomized treatment, i.e., To fit the corresponding joint model we first define the two transformation functions as: R> tf1 <-function (x, data) { + cbind(x, "^2" = x*x) + } R> tf2 <-function (x, data) { + cbind(x, "D-penicil" = x * (data$drug == 'D-penicil')) + } Following we update the call to jointFit.pbc12 and supply the list of the two functions in argument transFun, R> jointFit.pbc15 <-update(jointFit.pbc12, + transFun = list("value" = tf1, "extra" = tf2)) R> summary ( There is some weak evidence that the effect of η i (t) could be nonlinear, but clearly the association between η i (t) and the hazard does not seem to be different between the two treatment groups.

Supporting functions
Several supporting functions are available in the package that extract or calculate useful statistics based on the fitted joint model. In particular, function jointModelBayes() return objects of class "JMbayes", for which there are S3 methods defined for several of the standard generic functions in R. The most important are enlisted below: Functions coef() and fixef() extract the estimated coefficients for the two submodels from a fitted joint model. For the survival process both provide the same output, but for the longitudinal model, the former returns the subject-specific regression coefficients (i.e., the fixed effects plus their corresponding random effects estimates), whereas the latter only returns the estimated fixed effects.
Function ranef() extracts the empirical Bayes estimates for the random effects for each subject. The function also extracts estimates for the dispersion matrix of the posterior of the random effects using argument postVar.
Function vcov() extracts the estimated variance-covariance matrix of the parameters from the MCMC sample.
Functions fitted() and residuals() compute several kind of fitted values and residuals, respectively, for the two outcomes. For the longitudinal outcome the fitted() method computes the marginal Xβ and subject-specific Xβ + Zb i fitted values, whereβ and b i denote the posterior means of the fixed and random effects, whereas for the survival outcome it computes the cumulative hazard function for every subject and every time point a longitudinal measurement was collected. Analogously, method residuals() calculates the marginal and subject-specific residuals for the longitudinal outcome, and the martingale residuals for the survival outcome.
Function anova() can be used to compare joint models on the basis of the DIC, pD and LPML values.
Function plot() produces diagnostic plots for the MCMC sample, including trace plots, auto-correlation plots and kernel density estimation plots.
Function predict() produces subject-specific and marginal predictions for the longitudinal outcome, while function survfitJM() produces subject-specific predictions for the event time outcome, along with associated confidence or prediction confidence intervals. The use of these functions is explained in detail and illustrated in Section 5.
Function logLik() calculates the log-likelihood value for the posterior means of the parameters and the random effects, and can be also used to obtain the marginal log-likelihood (integrated over the parameters and random effects) using the Laplace approximation.
Function xtable() returns the L A T E X code to produce the table of posterior means, posterior standard deviations, and 95% credibility intervals from a fitted joint model. This is a method for the generic function xtable() from package xtable (Dahl 2014).

Definitions and estimation
In recent years there has been increasing interest in medical research towards personalized medicine. In particular, physicians would like to tailor decision making on the characteristics of individuals patients with aim to optimize medical care. In the same sense, patients who are informed about their individual health risk often decide to adjust their lifestyles to mitigate it. In this context it is often of interest to utilize results from tests performed on patients on a regular basis to derive medically-relevant summary measures, such as survival probabilities. Joint models constitute a valuable tool that can be used to derive such probabilities and also provide predictions for future biomarker levels. More specifically, under the Bayesian specification of the joint model, presented in Section 2, we can derive subjectspecific predictions for either the survival or longitudinal outcomes (Yu, Taylor, and Sandler 2008;Rizopoulos 2011Rizopoulos , 2012Taylor, Park, Ankerst, Proust-Lima, Williams, Kestin, Bae, Pickles, and Sandler 2013). To put it more formally, based on a joint model fitted in a sample D n = {T i , δ i , y i ; i = 1, . . . , n} from the target population, we are interested in deriving predictions for a new subject j from the same population that has provided a set of longitudinal measurements Y j (t) = {y j (t jl ); 0 ≤ t jl ≤ t, l = 1, . . . , n j }, and has a vector of baseline covariates w j . The fact that biomarker measurements have been recorded up to t, implies that subject j was event-free up to this time point, and therefore it is more relevant to focus on conditional subject-specific predictions, given survival up to t. In particular, for any time u > t we are interested in the probability that subject j will survive at least up to u, i.e., Similarly, for the longitudinal outcome we are interested in the predicted longitudinal response at u, i.e., The time-dynamic nature of both π j (u | t) and ω j (u | t) is evident from the fact that when new information is recorded for subject j at time t > t, we can update these predictions to obtain π j (u | t ) and ω j (u | t ), and therefore proceed in a time-dynamic manner.
Under the joint modeling framework of Section 2, estimation of either π j (u | t) or ω j (u | t) is based on the corresponding posterior predictive distributions, namely for the survival outcome, and analogously for the longitudinal one. The calculation of the first part of each integrand takes full advantage of the conditional independence assumptions (4) and (5). In particular, we observe that the first term of the integrand of π j (u | t) can be rewritten by noting that: whereas for ω j (u | t) we similarly have: Combining these equations with the MCMC sample from the posterior distribution of the parameters for the original data D n , we can devise a simple simulation scheme to obtain Monte Carlo estimates of π j (u | t) and ω j (u | t). More details can be found in Yu et al. (2008), Rizopoulos (2011Rizopoulos ( , 2012, and Taylor et al. (2013).
In package JMbayes these subject-specific predictions for the survival and longitudinal outcomes can be calculated using functions survfitJM() and predict(), respectively. As an illustration we show how these functions can be utilized to derive predictions for Patient 2 from the PBC dataset using joint model jointFit.pbc15. We first extract the data of this patient in a separate data frame

R> ND <-pbc2[pbc2$id == 2, ]
Function survfitJM() has two required arguments, the joint model object based on which predictions will be calculated, and the data frame with the available longitudinal data and baseline information. For Patient 2 estimates of π j (u | t) are calculated with the code: R> sfit.pbc15 <-survfitJM(jointFit.pbc15, newdata = ND) R> sfit.pbc15  (Time), quantile(Time, 0.9) + 0.01, length.out = 35) with Time denoting the observed event times variable. The user may override these default time points and specify her own using argument survTimes. In the output of survfitJM() we obtain as estimates of π j (u | t) the mean and median over the Monte Carlo samples along with the 95% pointwise confidence intervals. If only point estimates are of interest, survfitJM() provides the option (by setting argument simulate to FALSE) to use a first order estimator of π j (u | t) calculated as: whereθ denotes here the posterior means of the model parameters, andb j the mode of the posterior density p(b j | T * j > t, Y j (t),θ) with respect to b j . The corresponding plot() method for objects created by survfitJM() produces the figure of estimated conditional survival probabilities; for Patient 2 this is depicted in Figure 3. By setting logical argument include.y to TRUE, the fitted longitudinal profile is also included in the plot, i.e., R> plot(sfit.pbc15, estimator = "mean", include.y = TRUE, + conf.int = TRUE, fill.area = TRUE, col.area = "lightgrey") Argument estimator specifies whether the "mean" or the "median" over the 200 Monte Carlo samples should be used as an estimate of π j (u | t), and in addition arguments conf.int, fill.area and col.area control the appearance of the 95% confidence intervals. In a similar manner, predictions for the longitudinal outcome are calculated by the predict() function. For example, predictions of future log serum bilirubin levels for Patient 2 are produced with the code: R> Ps.pbc15 <-predict(jointFit.pbc15, newdata = ND, type = "Subject", + interval = "confidence", return = TRUE) Argument type specifies if subject-specific or marginal predictions are to be computed 1 , argument interval specifies the type of interval to compute (i.e., confidence or prediction), and by setting argument return to TRUE, predict() returns the data frame supplied in the required argument newdata having as extra columns the corresponding predictions and the limits of the confidence/prediction interval. This option facilitates plotting these predictions by a simple call to xyplot(), i.e.,

Web interface using shiny
To facilitate the use of package JMbayes for deriving individualized predictions, a web interface has been written using using package shiny (RStudio and Inc. 2014). This is available in the demo folder of the package and can be invoked with the code (assuming that JMbayes has been installed in the default library): R> library("shiny") R> runApp(file.path(.Library, "JMbayes/demo")) With this interface users may load an R workspace with the fitted joint model(s), following load the data of the new subject, and subsequently obtain dynamic estimates of π j (u | t) and ω j (u | t) (i.e., an estimate after each longitudinal measurement). Several additional options are provided to calculate predictions based on different joint models (if the R workspace contains more than one models), to obtain estimates at specific horizon times, and to extract the dataset with the estimated conditional survival probabilities.

Bayesian model averaging
Section 4.3 demonstrated that there are several choices to link the longitudinal and event time outcomes. When faced with this problem, the common practice in prognostic modeling is to base predictions on a single model that has been selected based on an automatic algorithm, such as, backward, forward or stepwise selection, or on likelihood-based information criteria, such as, AIC, BIC, DIC and their variants. However, what is often neglected in this procedure is the issue of model uncertainty. For example, if we choose a model using any of these criteria, say DIC, we usually treat it as the true model, even if there could be more than one models with DIC values of similar magnitude. In addition, when it comes to using a model for deriving predictions, we implicitly make the assumption that this model is adequate for all future patients. This seldom will be true in clinical practice. In our setting, a joint model with a specific formulation of the association structure may produce more accurate predictions for subjects with specific longitudinal profiles, while other models with other association structures may produce better predictions for subjects whose profiles have other characteristics. Here we follow another approach and we explicitly take into account model uncertainty by combining predictions under different association structures using Bayesian model averaging (BMA) (Hoeting, Madigan, Raftery, and Volinsky 1999;Rizopoulos, Hatfield, Carlin, and Takkenberg 2014).
We focus here on dynamic BMA predictions of survival probabilities. BMA predictions for the longitudinal outcome can be produced with similar methodology. Following the definitions of Section 5.1, we assume that we have available data D n = {T i , δ i , y i ; i = 1, . . . , n} based on which we fit M 1 , . . . , M K joint models with different association structures. Interest is in calculating predictions for a new subject j from the same population who has provided a set of longitudinal measurements Y j (t), and has a vector of baseline covariates w j . We let D j (t) = {T * j > t, Y j (t), w j } denote the available data for this subject. The model-averaged probability of subject j surviving time u > t, given her survival up to t is given by the expression: The first term in the right-hand side of (11) denotes the model-specific survival probabilities, derived in Section 5.1, and the second term denotes the posterior weights of each of the competing joint models. The unique characteristic of these weights is that they depend on the observed data of subject j, in contrast to classic applications of BMA where the model weights depend only on D n and are the same for all subjects. This means that, in our case, the model weights are both subject-and time-dependent, and therefore, for different subjects, and even for the same subject but at different times points, different models may have higher posterior probabilities . Hence, this framework is capable of better tailoring predictions to each subject than standard prognostic models, because at any time point we base risk assessments on the models that are more probable to describe the association between the observed longitudinal trajectory of a subject and the risk for an event.
For the calculation of the model weights we observe that these are written as : and p(D n | M k ) is defined analogously. The likelihood part p(D n | θ k ) is based on (6), and similarly p(D j (t) | θ k ) equals Thus, the subject-specific information in the model weights at time t comes from the available longitudinal measurements Y j (t) but also from the fact that this subject has survived up to t. We should note that the new subject j does not contribute any information about θ k (i.e., we do not refit the models using the data of this subject), the information for the parameters only comes from the original dataset in which the joint models have been fitted via the posterior distribution p(θ k | D n , M k ). A priori we assume that all models are equally probable, i.e., p(M k ) = 1/K, for all k = 1, . . . , K. Closed-form expressions for the marginal densities p(D n | M k ) and p(D j (t) | M k ) are obtained by means of Laplace approximations (Tierney and Kadane 1986) performed in two-steps, namely, first integrating out the random effects and then the parameters.
In package JMbayes BMA predictions for either the survival or longitudinal outcome can be calculated using function bma.combine(). This function accepts a series or a list of objects returned by either survfitJM() or predict() and a vector of posterior model weights, and returns a single object of the same class as the input objects with the combined predictions. We illustrate how this function can be used to produce the BMA prediction of π j (u | t) using the first five measurements of Patient 2 from the PBC dataset based on joint models jointFit.pbc1, jointFit.pbc12, jointFit.pbc13, jointFit.pbc14, and jointFit.pbc15. We start by computing the posterior model weights. As seen above for the calculation these weights we need to compute the marginal densities p(D n | M k ) and p(D j (t) | M k ). The former is obtained using the logLik() method for JMbayes objects, and the latter using function marglogLik(). The following code illustrates how this can be achieved: Argument newdata of marglogLik() is used to provide the available data D j (t) of the j-th subject, whereas argument marginal.thetas is invoked in order the logLik() method to compute the marginal log-likelihood. As just mentioned, we should stress that marglogLik() and logLik() compute log p(D j (t) | M k ) and log p(D n | M k ), respectively. Hence, to calculate the weights we need to transform them back to the original scale, i.e.,

Predictive accuracy
The assessment of the predictive performance of time-to-event models has received a lot of attention in the statistical literature. In general two main lines have emerged, namely one focusing on calibration, i.e., how well the model predicts the observed data (Schemper and Henderson 2000;Gerds and Schumacher 2006) and a second on focusing on discrimination, i.e., how well can the model discriminate between patients that had the event from patients that did not (Harrell, Kerry, and Mark 1996;Pencina, D'Agostino, D'Agostino, and Vasan 2008). In the following we present discrimination and calibration measures suitably adapted to the dynamic prediction setting and their implementation in JMbayes.

Discrimination
To measure the discriminative capability of a longitudinal marker we focus on a time interval of medical relevance within which the occurrence of events is of interest. In this setting, a useful property of the model would be to successfully discriminate between patients who are going to experience the event within this time frame from patients who will not. To put this formally, as before, we assume that we have collected longitudinal measurements Y j (t) = {y j (t jl ); 0 ≤ t jl ≤ t, l = 1, . . . , n j } up to time point t for subject j. We are interested in events occurring in the medically-relevant time frame (t, t + ∆t] within which the physician can take an action to improve the survival chance of the patient. Under the assumed model and the methodology presented in Section 5.1, we can define a prediction rule using π j (t + ∆t | t) that takes into account the available longitudinal measurements Y j (t). In particular, for any value c in [0, 1] we can term subject j as a case if π j (t + ∆t | t) ≤ c (i.e., occurrence of the event) and analogously as a control if π j (t + ∆t | t) > c. Thus, in this context, we define sensitivity and specificity as and Pr π j (t + ∆t | t) > c | T * j > t + ∆t , respectively. For a randomly chosen pair of subjects {i, j}, in which both subjects have provided measurements up to time t, the discriminative capability of the assumed model can be assessed by the area under the receiver operating characteristic curve (AUC), which is obtained for varying c and equals, that is, if subject i experiences the event within the relevant time frame whereas subject j does not, then we would expect the assumed model to assign higher probability of surviving longer than t+∆t for the subject who did not experience the event. To summarize the discriminating power of the assumed model over the whole follow-up period, we need to take into account that the number of subjects contributing to the comparison of the fitted π i (t + ∆t | t) with the observed data is not the same for all time points t. Following an approach similar to Antolini, Boracchi, and Biganzoli (2005) and Heagerty and Zheng (2005), we can utilize a weighted average of AUCs, i.e., where E(t) = {T * i ∈ (t, t + ∆t]} ∩ {T * j > t + ∆t} , and Pr{E(t)} denotes the probability that a random pair is comparable at t. We can call C ∆t dyn a dynamic concordance index since it summarizes the concordance probabilities over the follow-up period. Note also that AUC(t, ∆t) and as a result also C ∆t dyn depend on the length ∆t of the time interval of interest, which implies that different models may exhibit different discrimination power for different ∆t.
For the estimation of AUC(t, ∆t) and C ∆t dyn we need to take care of two issues, namely, the calculation of the integrals in the definition of (12) and censoring. For the former we use the 15-point Gauss-Kronrod quadrature rule. Estimation of AUC(t, ∆t) is directly based on its definition, namely by appropriately counting the concordant pairs of subjects. More specifically, we have A UC(t, ∆t) = A UC 1 (t, ∆t) + A UC 2 (t, ∆t).
A UC 1 (t, ∆t) refers to the pairs of subjects who are comparable (i.e, their observed event times can be ordered), where i, j = 1, . . . , n with i = j. For such comparable subjects i and j, we can estimate and compare their survival probabilities π i (t+∆t | t) and π j (t+∆t | t), based on the methodology presented in Section 5.1. This leads to a natural estimator for AUC 1 (t, ∆t) as the proportion of concordant subjects out of the set of comparable subjects at time t: where I(·) denotes the indicator function. Analogously, A UC 2 (t, ∆t) refers to the pairs of subjects who due to censoring cannot be compared, namely with again i, j = 1, . . . , n with i = j. Concordant subjects in this set contribute to the overall AUC appropriately weighted with the probability that they would be comparable, i.e., , withν i (t + ∆t | T i ) = 1 −π i (t + ∆t | T i ) being the probability that subject i who survived up to time T i will have the event before t + ∆t.
Having estimated AUC(t, ∆t), the next step in estimating C ∆t dyn is to obtain estimates for the weights Pr{E(t)}. We observe that these can be rewritten as where the simplification in the second line comes from the independence of subjects i and j, and S(·) here denotes the marginal survival function. In practice calculation of C ∆t dyn is restricted into a follow-up interval [0, t max ] where we have information. Let t 1 , . . . , t 15 denote the re-scaled abscissas of the Gauss-Kronrod rule in the interval [0, t max ] with corresponding weights 1 , . . . , 15 . We combine the estimates A UC(t k , ∆t), k = 1, . . . , 15 with the estimates of the weights Pr{E(t)} to obtain where Pr{E(t k )} = S(t k ) − S(t k + ∆t) S(t k + ∆t), with S(·) denoting here the Kaplan-Meier estimate of the marginal survival function S(·).
The AUC(t, ∆t) and the dynamic discrimination index can be calculated for joint models fitted by jointModelBayes() using functions aucJM() and dynCJM(), respectively. We illustrate their use based again on joint model jointFit.pbc15. The basic call to aucJM() requires the user to provide the fitted joint model object, the data frame upon which the AUC is to be calculated, the time point t (argument Tstart) up to which longitudinal measurements are to be used and the length of the time window ∆t (argument Dt) 2 : R> auc.pbc15 <-aucJM(jointFit.pbc15, newdata = pbc2, Tstart = 5, Dt = 2) R> auc.pbc15 Time-dependent AUC for the Joint Model jointFit.pbc15 Estimated AUC: 0.842 At time: 7 Using information up to time: 5 (202 subjects still at risk) We observe that using the first five year longitudinal measurements, serum bilirubin exhibits nice discrimination capabilities for patients who are to die within a two-year time frame. To investigate if this is also the case during the whole follow-up period, we calculate the dynamic discrimination index for the same time window. The syntax of dynCJM() is (almost) identical to the one of aucJM(), i.e., R> dynC.pbc15 <-dynCJM(jointFit.pbc15, newdata = pbc2, Dt = 2) R> dynC.pbc15 Dynamic Discrimination Index for the Joint Model jointFit.pbc15 The estimate of C ∆t=2 dyn is almost identical to the one of AUC(t = 5, ∆t = 2) indicating that serum bilirubin can discriminate well between patients during follow-up.

Prediction error
The assessment of the accuracy of predictions of survival models is typically based on the expected error of predicting future events. In our setting, and again taking into account the dynamic nature of the longitudinal outcome, it is of interest to predict the occurrence of events at u > t given the information we have recorded up to time t. This gives rise to expected prediction error: where N i (t) = I(T * i > t) is the event status at time t, L(·) denotes a loss function, such as the absolute or square loss, and the expectation is taken with respect to the distribution of the event times. An estimate of PE(u | t) that accounts for censoring has been proposed by Henderson, Diggle, and Dobson (2002): where n(t) denotes the number of subjects at risk at time t. The first two terms in the sum correspond to patients who were alive after time u and dead before u, respectively; the third term corresponds to patients who were censored in the interval [t, u]. Using the longitudinal information up to time t, PE(u | t) measures the predictive accuracy at the specific time point u. Alternatively, we could summarize the error of prediction in a specific interval of interest, say [t, u], by calculating a weighted average of {PE(s | t), t < s < u} that corrects for censoring, similarly to C ∆t dyn . An estimator of this type for the integrated prediction error has been suggested by Schemper and Henderson (2000), which adapted to our time-dynamic setting takes the form I PE(u | t) = i:t≤T i ≤u δ i S C (t)/ S C (T i ) PE(T i | t) i:t≤T i ≤u δ i S C (t)/ S C (T i ) , where S C (·) denotes the Kaplan-Meier estimator of the censoring time distribution.
Both PE and IPE can be calculated for joint models fitted by jointModelBayes() using function prederrJM(). This has a similar syntax as function aucJM(), and requires a fitted joint model, a data frame based on which the prediction error will be calculated, and the time points t (argument Tstart) and u (argument Thoriz) that denotes up to which time point to use the longitudinal information and at which time point to make the prediction, respectively. For model jointFit.pbc15 using the biomarker information during the first five years of follow-up the estimated prediction error at year seven is R> pe.pbc15 <-prederrJM(jointFit.pbc15, pbc2, Tstart = 5, Thoriz = 7) R> pe.pbc15 Prediction Error for the Joint Model jointFit.pbc15 Estimated prediction error: 0.107 At time: 7 Using information up to time: 5 (202 subjects still at risk) Loss function: square By default the loss function is the square one (i.e., L(x) = x 2 ), but the user may specify the absolute loss or define her own loss function using argument lossFun. The integrated prediction error can be simply calculated by setting logical argument interval to TRUE in the call to prederrJM(); for example, for the same joint model and in the interval [5,9] the IPE is calculated with the code: R> ipe.pbc15 <-prederrJM(jointFit.pbc15, pbc2, Tstart = 5, + Thoriz = 9, interval = TRUE) R> ipe.pbc15 Prediction Error for the Joint Model jointFit.pbc15 Estimated prediction error: 0.0907 In the time interval: [5,9] Using information up to time: 5 (202 subjects still at risk) Loss function: square

Future plans
In this paper we have illustrated the capabilities of package JMbayes for fitting joint models for longitudinal and time-to-event data under a Bayesian approach. As we have seen, the current version of the package provides several options for fitting different types of joint models, but nonetheless several extensions are planned in the future to further expand on what is currently available. These include among others: • The consideration of multiple longitudinal outcomes, while allowing for the various association structures we have presented in Sections 4.3 and 4.4.
• Handling of exogenous time-varying covariates by supplying a time-dependent Cox model as an argument to jointModelBayes().
• Extend functionality in the survival submodel to handle, competing risks, recurrent events, and left-and interval-censored event time data.
• Update dynamic predictions to handle the aforementioned extensions.