Simulation-Based Calibration for Joint Models

Why Calibrate a Joint-Model Computation?

A Bayesian joint model can be difficult to fit: it combines repeated measurements, event-time information, latent random effects, and often a flexible baseline hazard. Chain diagnostics can reveal some computational problems, but they do not by themselves check whether the fitting algorithm is returning draws from the posterior distribution specified by the model. Simulation-based calibration (SBC) provides a direct test of that computation under controlled conditions: generate parameters and data from the model, refit the model, and compare the generating parameter with posterior draws.

SBC asks whether an inference procedure is calibrated when the model is true and its parameters are generated from the stated prior. It is not a test that the model is a good description of an external population. That distinction matters especially in a joint model: a flat rank histogram supports the calibration of the implementation under the simulation design, but does not establish that the longitudinal or event-time assumptions are appropriate for observed data.

Joint-Model Setup and Calibration Assumptions

Consider subjects \(i=1,\ldots,n\), with repeated measurements \(Y_{ij}\) at times \(t_{ij}\), and an event time observed subject to right censoring. A simple shared-random-effects model has a longitudinal mean \(m_i(t)=\beta_0+\beta_1t+b_{0i}+b_{1i}t\), with \(Y_{ij}\mid b_i,\theta\sim N(m_i(t_{ij}),\sigma_y^2)\) and \(b_i=(b_{0i},b_{1i})^\mathsf{T}\sim N(0,D)\). An example of a linked event hazard is \(h_i(t\mid b_i,\theta)=h_0(t)\exp\{\gamma^\mathsf{T}w_i+\alpha m_i(t)\}\). Here \(\theta\) denotes the population parameters, including the longitudinal coefficients, association parameter \(\alpha\), hazard parameters, residual variance, and random-effects covariance parameters.

For event time \(T_i\) and censoring time \(C_i\), the observed time is \(U_i=\min(T_i,C_i)\), and the event indicator is \(\delta_i=I(T_i\le C_i)\). Conditional on the modeled covariates and random effects, the event contribution to the likelihood is \(h_i(U_i)^{\delta_i}\exp\{-\int_0^{U_i}h_i(s)\,ds\}\). The longitudinal and event contributions are combined and the random effects integrated over or sampled as part of the joint model. In an SBC exercise, the data-generating mechanism, including the measurement schedule and censoring process, must match the assumptions used in fitting. In particular, censoring must be generated consistently with the conditional noninformative-censoring assumption used by the analysis.

The benchmark relies on several conditions. The generating parameter must actually be drawn from the prior used in fitting; data must be simulated from the stated likelihood and censoring mechanism; and the fitted posterior draws must represent the target posterior. For the simplest rank argument, posterior draws are conditionally independent draws from the exact posterior. Practical MCMC draws are often autocorrelated, so inadequate mixing or too few effective draws can distort rank behavior. SBC should therefore be used alongside chain diagnostics, not instead of them.

Why the Rank Should Be Uniform

For a scalar continuous parameter \(\theta_k\), draw its generating value \(\theta_k^{(0)}\), simulate data \(y\), and obtain \(L\) posterior draws \(\theta_k^{(1)},\ldots,\theta_k^{(L)}\). Define the rank \(R_k=\sum_{\ell=1}^{L}I(\theta_k^{(\ell)}<\theta_k^{(0)})\). It can take any integer from 0 through \(L\). The generating value is included conceptually in the comparison, so there are \(L+1\) possible positions.

The justification follows from the Bayesian joint distribution. Before seeing the simulated data, \(\theta_k^{(0)}\) is drawn from its prior and \(y\) is drawn from the likelihood given the full generating parameter. Conditional on \(y\), the generating parameter has the posterior distribution by Bayes' rule. Thus, conditional on \(y\), \(\theta_k^{(0)}\) and \(L\) independent posterior draws are identically distributed and exchangeable. With a continuous parameter, ties have probability zero, and each of the \(L+1\) possible order positions is equally likely. Averaging over simulated data preserves that uniform rank distribution.

This argument applies to scalar quantities whose induced distribution is continuous, including parameters and scalar functions of parameters. It does not require the quantity to be a basic coefficient: for example, a hazard at a specified time can be evaluated for the generating parameter and for each posterior draw, then ranked. For quantities with ties, such as discrete-valued functions of parameters or rounded output, a standard remedy is to randomize the rank uniformly over the positions occupied by tied values; otherwise, a nonuniform rank histogram may be an artifact of tie handling.

Running SBC in Practice

The full hierarchy matters. In a shared-random-effects model, the random effects are generated from their population distribution, and the longitudinal and event outcomes are then simulated conditional on them. Generating the longitudinal data from random effects but simulating event times independently of those effects would test a different model than the one being fitted. Similarly, if the actual analysis includes right censoring, the synthetic data must include censoring handled in the same way.

Summarizing Rank Histograms

For \(S\) simulation replicates and \(L\) posterior draws per replicate, there are \(L+1\) rank categories. Under calibration, the expected count in each category is \(S/(L+1)\). A histogram that is approximately level is consistent with calibration, subject to sampling variation. A systematic slope can indicate a location bias: ranks near one end occur more often when generating values tend to lie below or above their posterior draws. Excess counts at both extremes can indicate posterior draws that are too concentrated relative to the generating values; excess counts near the middle can indicate draws that are too dispersed. These patterns are clues, not unique diagnoses, and can also arise from poor mixing or implementation errors.

For a formal summary, a Pearson statistic can compare observed counts \(O_r\) with expected counts \(E_r=S/(L+1)\): \(\sum_{r=0}^{L}(O_r-E_r)^2/E_r\). Its usual chi-square approximation is more credible when expected counts are not small, and its interpretation must account for the many parameters or derived quantities examined. A histogram is often more informative than a single p-value: a non-significant test is not proof of calibration, while a large departure should prompt investigation of both the model code and the posterior computation.

Ranking Derived Quantities and Diagnosing Failure

A model parameter can be poorly calibrated even when a familiar coefficient appears stable, and a derived quantity can reveal errors in transformations or parameter combinations. For example, calculate a subject-specific hazard at a fixed time from each posterior draw and from the generating values, then rank those quantities across replicates. The calculation must use the same definition in simulation and fitting, including units and covariates.

When rank plots look unusual, first verify the simulation and fitting code agree on parameterization, units, random-effects structure, censoring, and likelihood contributions. Then examine chain behavior, effective sample size, and Monte Carlo error for the ranks. A rank trend can result from biased posterior computation; extreme-heavy ranks can reflect under-dispersed draws, but can also arise from insufficiently mixed chains. Rank histograms for individual quantities are marginal checks: they do not guarantee that the joint dependence among parameters is calibrated.

SBC differs from chain diagnostics and simulation-recovery summaries. Chain diagnostics assess whether sampled chains have mixed and agree; SBC assesses whether the overall procedure recovers the known prior-generated value in rank across repeated model-generated data sets. A typical recovery exercise may compare estimates with generating values, but SBC uses a posterior-rank property that does not depend on choosing a point estimator or a nominal interval. Conversely, SBC is computationally intensive and only tests cases drawn from the chosen prior and design. It cannot establish robustness to misspecification or guarantee performance in regions the prior rarely generates.

Common Pitfalls and Limitations

• Using the wrong generating distribution. If parameters are generated from a convenient distribution different from the fitting prior, the uniform-rank result does not apply. Draw from the actual prior hierarchy.
• Simulating a simplified data structure. Omitting shared random effects, censoring, or the longitudinal-event dependence changes the target. Simulate the full joint data-generating process used in the analysis.
• Treating autocorrelated draws as independent. The exchangeability justification assumes exact independent posterior draws. MCMC autocorrelation can affect rank histograms, so run adequate chains and assess effective sample size rather than relying on nominal draw counts.
• Interpreting a flat histogram as proof of model fit. SBC checks computation under the assumed model; it does not show that the model describes an external cohort adequately.
• Over-reading individual bins or p-values. Random variation is expected. Use enough independent replicates, inspect the pattern, and consider how many targets were tested.