An alternative approach to joint modelling is to use Bayesian MCMC methods. Such an approach has been considered previously by Faucett and Thomas (

1996), Brown et al. (

2005), Brown and Ibrahim (

2003), and Chi and Ibrahim (

2006) where tailor-made Gibbs sampling algorithms have been developed. In contrast, Guo and Carlin (

2004) use generic MCMC software such as WinBUGS (Spiegelhalter et al.,

2003) to implement a joint model. The attraction of using WinBUGS is that a joint model can be coded relatively simply, which applied researchers can then easily adapt to suit their needs. There is also considerable flexibility in model specification. For example, an individual's underlying intercept, slope or current value can all be associated with the time-to-event process. Second, standard errors of the model parameters can be obtained from posterior distributions without the need for asymptotic theory. Finally, in and out-of-sample predictions of both longitudinal and time-to-event processes can be obtained relatively easily using posterior predictive distributions.

However, a current limitation of implementing a joint model in WinBUGS is that, for the time-to-event process, a closed-form evaluation of the integrated hazard function is required in order to formulate the likelihood. For this reason, our Bayesian model is currently restricted to joint models for which the longitudinal process is specified by a linear mixed-effects model, so that *m*_{i}(*t*)=β_{11}+β_{12}*t*+*U*_{1i}+*U*_{2i}*t* while the hazard for the time-to-event process is modelled using a constant baseline hazard. The baseline hazard is multiplied by a function of the individual-specific intercept, slope and underlying measurement, as follows:

Under this specification the integrated hazard has the following closed-form expression:

More generally, to relax the assumption of a constant baseline hazard, a piecewise-constant baseline hazard can be specified and the integrated hazard evaluated within each time band separately. Furthermore, to allow for a wider range of hazards and non-linear longitudinal profiles, the trapezium rule could be used to approximate the integrated hazard function, an approach taken by Brown et al. (

2005). The WinBUGS code used to implement a constant baseline hazard joint model with linear longitudinal process is given as Supporting Information on the journal website. We use minimally informative Normal(0, 10

^{4}) priors for all fixed effect parameters (β

_{11}, β

_{12}, log(λ), α

_{1}, α

_{2} and α

_{3}), a Uniform(0, 100) prior for the within patient standard deviation, and a scaled inverse-Wishart prior for the between patient variance–covariance matrix, with degrees of freedom equal to one plus the number of random effects, as suggested previously for the specification of mixed-effect models (Gelman and Hill,

2007).

Inferences concerning future longitudinal values from a new individual (an out-of-sample prediction) can be obtained from the posterior predictive distribution. The existing data on

*m* patients can be summarised by

recalling that

**y**_{i}=

*y*_{i}(

*t*_{ij}),

*j*=1,…,

*n*_{i},

, and

. Suppose now that a set of longitudinal measurements are obtained from a new individual, labelled

*m*+1, together with information that they have survived up to time

*s*. This can be summarised as

. Then, given this data, the predictive distribution for a new observation

from this individual with random effects

and hyperparameters

is

Similarly, the predicted survival probability, for the time to failure

for the new individual, at time

*t* given survival up to time

*s* is

where

is the posterior distribution of the random effects for the new individual conditional on their data and the hyperparameters

**θ**. These quantities are easily obtained within the WinBUGS software, and as demonstrated in the code given as Supporting Information.