Fitting Bayesian Joint Models With MCMC

Why Use MCMC for a Joint Model?

A shared-random-effects joint model links a subject’s longitudinal trajectory to their event risk through latent subject-specific effects. Maximum-likelihood fitting integrates over those random effects, often using numerical approximations. Bayesian fitting instead assigns prior distributions to model parameters and uses Markov chain Monte Carlo (MCMC) to approximate their joint posterior distribution. This gives direct posterior summaries for an association parameter, while requiring care about prior choices and whether the simulation has adequately explored the posterior.

The focus here is fitting and assessing MCMC for a shared-random-effects model. The event association will be through the marker’s current latent value. Other association structures require corresponding changes to the model specification; the interpretation of an association coefficient always depends on the quantity it multiplies.

Model and Assumptions

For subject \(i\), let \(Y_{ij}\) be the marker measured at time \(t_{ij}\), \(T_i\) the observed event or censoring time, and \(\delta_i\) the event indicator, equal to 1 for an observed event and 0 for right censoring. Let \(b_i\) be a vector of subject-specific random effects. A Gaussian longitudinal submodel and a current-value event submodel can be written as follows:

Here \(m_i(t)\) is the latent marker trajectory, \(w_i\) contains event-submodel covariates, and \(h_0(t)\) is the baseline hazard. The association parameter \(\eta\) describes how the current latent marker value relates to the hazard, conditional on the other model terms. For a marker difference of \(\Delta\) units, the corresponding conditional hazard ratio is \(\exp(\eta\Delta)\). Thus, the units and scaling of the marker matter.

The usual assumptions include independent subjects; the specified longitudinal and event distributions conditional on the random effects and covariates; and a censoring process that is noninformative conditional on the modeled information. The longitudinal measurement errors are conditionally independent with the stated variance. The visit schedule is treated as fixed or noninformative under the model. A piecewise-constant baseline hazard is one practical choice: it assigns a constant hazard \(\lambda_k\) to each prespecified time interval \(k\).

These assumptions are substantive, not consequences of running MCMC. For example, if the event process still depends on a subject-specific feature after conditioning on \(b_i\) and the covariates, the stated event model may not capture that dependence. MCMC can sample from the posterior of a specified model; it does not make a misspecified model correct.

Specify Priors Before Sampling

A Bayesian fit combines the likelihood with priors to define the posterior. One reasonable starting specification uses normal priors for regression coefficients and the association parameter, a positive prior for the residual standard deviation, and a structured prior for the random-effects covariance matrix. For instance:

A half-\(t\) distribution is a \(t\) distribution restricted to positive values. For \(D\), a common strategy separates random-effect standard deviations from their correlation matrix, assigns half-\(t\) priors to the standard deviations, and uses an LKJ prior for the correlation matrix. For a piecewise-constant baseline hazard, one may place normal priors on \(\log(\lambda_k)\), ensuring each \(\lambda_k\) is positive. These are examples of prior families, not universal defaults: scales should be chosen with the units and plausible ranges of the problem in mind.

Continuous predictors are often centered and scaled before fitting, with the same transformation used consistently in both submodels. A prior on \(\eta\) then has a more interpretable scale, and sampling may be better behaved. Prior predictive checks—simulating trajectories and event times from the priors before fitting—can reveal implausible implications, such as nearly all events occurring immediately or trajectories far outside plausible ranges.

From the Joint Likelihood to Posterior Draws

Conditional on \(b_i\), the subject’s event contribution is the event hazard raised to the event indicator, multiplied by the survival probability through the observed time. If \(H_i(t\mid b_i)\) is the cumulative hazard, that contribution is \(h_i(T_i\mid b_i)^{\delta_i}\exp\{-H_i(T_i\mid b_i)\}\). Multiplying it by the longitudinal densities and the random-effects density gives the subject’s conditional contribution.

The full posterior is proportional to the product of these contributions over subjects, multiplied by the priors for all parameters. This follows directly from Bayes’ rule: posterior density is proportional to likelihood times prior density. The random effects are unknown latent quantities and can be sampled along with the population-level parameters.

The computational reason MCMC is useful is that the marginal likelihood would require integrating over all subjects’ random effects and, often, other parameters. Instead, a sampler generates dependent draws from the joint posterior of the parameters and latent effects. After initial draws used for warm-up, retained draws approximate the posterior distribution if the chains have reached and adequately explored their stationary distribution. This is an approximation that must be assessed, not a guarantee supplied by the algorithm.

Hamiltonian Monte Carlo, including its No-U-Turn Sampler (NUTS) variant, is widely used for continuous hierarchical models. It proposes moves using posterior gradients and can explore correlated parameter spaces more efficiently than simple random-walk proposals. Other samplers may combine Gibbs updates for parameters with convenient conditional distributions and Metropolis updates for the rest. The appropriate implementation depends on the software and model; whatever the algorithm, retain multiple chains initialized at dispersed plausible values and examine the resulting draws.

Trace Plots, Effective Sample Size, and Precision

A trace plot displays sampled values against iteration. For key quantities—especially \(\eta\), important variance components, and selected latent effects—look for chains that overlap and fluctuate across a stable range without persistent trends or long periods trapped in distinct regions. A trace plot is a visual diagnostic: a reassuring appearance does not prove convergence, and a troubling appearance is reason to investigate rather than to report an unqualified summary.

Because successive MCMC draws are dependent, the retained draw count is not the same as the number of independent draws. For a scalar quantity with \(M\) retained draws and autocorrelation \(\rho_k\) at lag \(k\), an approximate effective sample size is:

In practice, software estimates this using a suitable truncation or positive-sequence method. A small ESS means the chain contains less information for estimating that quantity than the raw draw count suggests. The approximate Monte Carlo standard error for a posterior mean is the posterior standard deviation divided by \(\sqrt{ESS}\). It measures simulation precision, not uncertainty about the parameter itself.

Use multiple-chain convergence summaries as a complement to plots. The split-\(\hat R\) statistic compares variation within and between chains after splitting chains into halves; values close to 1 indicate that those sources of variation are similar. A value such as 1.01 is often used as a practical screening target, not as proof of convergence. Examine ESS and Monte Carlo error as well. Diagnostics should cover the parameters that matter to the scientific question, not just an easy-to-sample regression coefficient.

Practical Fitting Sequence and Alternatives

MCMC and maximum likelihood are different inferential approaches, not interchangeable computational settings. Maximum likelihood targets parameter values that maximize the likelihood, often using quadrature or approximations to integrate random effects. Bayesian MCMC targets a posterior distribution that also reflects prior information, and can provide direct uncertainty summaries for latent quantities. Bayesian credible intervals have a posterior probability interpretation conditional on the model and prior; frequentist confidence intervals have a different repeated-sampling interpretation. Neither approach removes the need to check the model’s assumptions or the stability of the computation.

Common Pitfalls and Limitations

• Reading a coefficient without its scale. The hazard ratio is \(\exp(\eta\Delta)\), not simply \(\eta\). State the marker difference, units, and conditioning context.
• Treating many draws as proof of convergence. Autocorrelation can make the effective sample size much smaller than the retained count. Review traces, multiple-chain summaries, ESS, and Monte Carlo error.
• Assuming a good trace validates the model. Sampling diagnostics address computation, not whether the longitudinal distribution, event association, censoring assumption, or baseline hazard is appropriate.
• Ignoring prior influence. With limited information, posterior results can depend materially on prior scales. Check prior implications and consider whether conclusions change under scientifically defensible alternatives.

MCMC can be costly for large or highly correlated joint models, and diagnostics cannot establish that every posterior region has been explored. Posterior summaries should therefore be reported with adequate simulation diagnostics and in context. If chains show persistent trends, do not present the resulting summaries as reliable merely because the software completed.