Posterior Predictive Checks for Joint Models

Why Check Replicated Data?

A joint model connects repeated biomarker measurements to an event process, often through shared subject-specific random effects or a latent biomarker trajectory. After fitting the model, a central diagnostic question is whether data generated under its fitted structure can resemble the data actually observed. A model may estimate parameters precisely yet fail to reproduce features that matter, such as the spread of individual slopes, the number of early events, or the relationship between biomarker trajectories and event timing.

A posterior predictive check answers this question by generating replicated data sets from posterior draws and comparing summaries of those replicas with the corresponding summaries of the observed data. These checks are descriptive diagnostics: they can reveal mismatches, but agreement does not prove that the model is correct.

Definition: Let \(D\) denote the observed longitudinal and event-time data, and let \(\theta\) collect model parameters. A posterior predictive replicate \(D^{\mathrm{rep}}\) is generated by first drawing \(\theta^{(s)}\) from the posterior distribution \(p(\theta\mid D)\), then generating data from the model conditional on that draw. A discrepancy statistic \(T(D)\) is a summary chosen to highlight a feature of the data; a posterior predictive check compares \(T(D)\) with \(T(D^{\mathrm{rep}})\) across replicated data sets.

What Must Be Replicated?

For subject \(i\), write the observed biomarker measurements as \(Y_{ij}\) at visit times \(t_{ij}\), and let \(T_i\) be the observed follow-up time and \(\Delta_i\) the event indicator, with \(\Delta_i=1\) for an observed event and \(0\) for censoring. A common shared-random-effects structure has a longitudinal model such as

Y_{ij}=m_i(t_{ij};\theta,b_i)+\epsilon_{ij},

and an event hazard such as

h_i(t\mid\theta,b_i)=h_0(t;\theta)\exp\{x_i^\mathsf{T}\gamma+\alpha m_i(t;\theta,b_i)\}.

Here \(b_i\) is the subject-specific random effect, \(m_i(t;\theta,b_i)\) is the subject's latent biomarker trajectory, \(h_0(t;\theta)\) is the baseline hazard, \(x_i\) contains event-model covariates, and \(\alpha\) links the trajectory to the event hazard. This is an illustrative model form; the check must use the actual longitudinal and event submodels that were fitted.

For a check conditional on the observed subjects, draw \((\theta^{(s)},b_1^{(s)},\ldots,b_n^{(s)})\) from their joint posterior, generate new measurement errors and event times, and retain the observed covariates and visit schedule. This asks whether the model can reproduce the data for these subjects, given the information their measurements and outcomes have already supplied about their random effects. For a check aimed at new subjects, instead draw new random effects from the fitted population distribution. The two checks answer different questions: conditioning on estimated subject effects can make individual trajectories look good while failing to expose a poor model for the population distribution of those effects.

Event replication must also respect the censoring design. With fixed administrative censoring times \(C_i\), simulate a latent event time \(T_i^*\), then set \(T_i^{\mathrm{rep}}=\min(T_i^*,C_i)\) and \(\Delta_i^{\mathrm{rep}}=I(T_i^*\leq C_i)\). If censoring times are random and are to be replicated, a censoring model is needed. Treating observed follow-up times as event times, or comparing event times without accounting for censoring, does not reproduce the observed-data process.

Choose the predictive target.
Decide whether the check is conditional on observed subjects and their posterior random effects, or concerns new subjects drawn from the population distribution.

Draw from the fitted posterior.
For each retained posterior draw, use a coherent draw of the model parameters and, for conditional checks, the subject random effects.

Generate both data components.
Simulate biomarker measurements at the chosen visit times and event outcomes from the event model using the same draw, then apply the censoring design.

Compare relevant summaries.
Calculate the same trajectory and event-pattern summaries for observed and replicated data, and look for systematic discrepancies across draws.

Worked Checks

Worked Example: Comparing Replicated Trajectories

Hypothetical data. Suppose three subjects have visits at times \(0,1,2\) months. For one posterior draw, the fitted trajectories are \(m_i(t)=b_{0i}+b_{1i}t\), with \((b_{01},b_{11})=(100,2)\), \((b_{02},b_{12})=(110,-1)\), and \((b_{03},b_{13})=(90,3)\), in mg/dL and mg/dL per month. Let the residual standard deviation be \(4\) mg/dL. The observed measurements are \((101,103,105)\), \((109,106,109)\), and \((92,95,95)\) mg/dL.

For this draw, use hypothetical replicated residuals \((1,-2,0.5)\), \((-3,0,2)\), and \((2,1,-1)\) mg/dL. Adding each residual to its trajectory gives replicated measurements \((101,100,104.5)\), \((107,109,110)\), and \((92,94,95)\) mg/dL.

The observed mean at time \(0\) is \((101+109+92)/3=100.7\) mg/dL; at time \(1\) it is \((103+106+95)/3=101.3\) mg/dL; and at time \(2\) it is \((105+109+95)/3=103.0\) mg/dL. The corresponding replicated means are \((101+107+92)/3=100.0\), \((100+109+94)/3=101.0\), and \((104.5+110+95)/3=103.2\) mg/dL. These means are similar for this draw. However, individual trajectories also matter: the observed endpoint slopes are \(2.0\), \(0\), and \(1.5\) mg/dL per month, with mean \(1.17\); the replicated endpoint slopes are \(1.75\), \(1.5\), and \(1.5\), with mean \(1.58\) mg/dL per month. One replicate is not enough to assess fit. Across many draws, compare the distribution of these summaries and the spread of subject-specific slopes with their observed values.

Worked Example: Replicating Event Outcomes Under Administrative Censoring

Hypothetical setup. Consider three subjects followed for at most 12 months. For a particular posterior draw, suppose their hazards are constant at \(0.10\), \(0.05\), and \(0.20\) events per month. These may differ because their covariates or shared random effects differ. With a constant hazard \(h_i\), the event time has an exponential distribution, and inverse-transform simulation gives \(T_i^*=-\log(U_i)/h_i\) for \(U_i\) drawn uniformly between 0 and 1. Use hypothetical uniforms \(0.20\), \(0.50\), and \(0.80\).

The simulated latent event times are \(-\log(0.20)/0.10=16.09\) months, \(-\log(0.50)/0.05=13.86\) months, and \(-\log(0.80)/0.20=1.116\) months. Applying the 12-month administrative limit gives observed replicated follow-up times of \(12\), \(12\), and \(1.116\) months, with event indicators \(0,0,1\). Thus, one of the three subjects has a replicated event by month 12.

If the observed data had two events by month 12, this particular replicate would have fewer events than observed. That difference alone does not establish a mismatch: the relevant comparison is the observed event count against its distribution across many posterior replicates. The censoring rule is part of the calculation; counting the first two latent event times as events would be wrong because both occur after the administrative censoring time.

Discrepancy Statistics and Their Interpretation

A discrepancy should target a scientifically or statistically important feature and be calculated the same way for the observed data and each replicate. For longitudinal data, useful summaries include the mean trajectory at selected times, between-subject variation in slopes, within-subject residual spread, or the frequency of unusually large changes between visits. For the event component, summaries can include the number of events by a prespecified time, the distribution of observed follow-up times with censoring represented, or event counts within covariate groups. A joint check can examine a feature connecting the two components, such as whether subjects with larger estimated biomarker slopes also tend to have earlier events.

For a scalar discrepancy, a posterior predictive tail area can be estimated from \(S\) posterior draws as

\widehat{p}_{\mathrm{pp}}=\frac{1}{S}\sum_{s=1}^S I\{T(D^{\mathrm{rep},(s)})\geq T(D)\}.

This is the fraction of replicated discrepancies at least as large as the observed discrepancy. It is a descriptive tail-area summary, not a classical p-value: the observed data were used to form the posterior, and posterior predictive tail areas are not generally uniformly distributed under a true model. A value near zero or one can flag a feature that the model often fails to reproduce in one direction, but it is not the probability that the model is true or false.

Worked Example: Assessing an Event-Count Tail Area

Hypothetical data. In a cohort of 40 subjects, 8 events are observed by month 12. From \(S=1{,}000\) posterior predictive replicates, suppose 72 have at least 8 events by month 12. For the discrepancy \(T(D)=\) number of events by month 12, the estimated upper-tail area is

\widehat{p}_{\mathrm{pp}}=\frac{72}{1{,}000}=0.072.

The observed event count is in the upper 7.2% of the replicated counts. This suggests checking whether the fitted model tends to generate fewer events than observed, but the result is not a formal rejection threshold. The approximate Monte Carlo standard error of this estimated proportion is

\sqrt{\frac{0.072(1-0.072)}{1{,}000}}=0.00817.

This measures simulation error in estimating the tail area, not uncertainty in the model parameters and not a measure of model fit. A histogram or interval display of all 1,000 replicated counts is more informative than the tail area alone: it shows the full spread and whether 8 is far from typical values.

Plots, Failure Modes, and Alternatives

Plotting observed summaries against replicated distributions makes patterns easier to detect than relying on a single number. For trajectories, overlay the observed mean trajectory with a pointwise envelope of replicated mean trajectories, and separately examine subject-level slopes or residual spread. For event outcomes, show the observed event count or event-time summary against the replicated distribution, preserving censoring in each replicate. Stratifying by an important covariate can reveal an overall fit that hides systematic underprediction in one group and overprediction in another.

A useful diagnostic should match the question. If the goal is population-level calibration, conditional checks that reuse each subject's posterior random effect may be too forgiving. If the goal is whether the model reproduces the observed cohort conditional on its covariates and design, checks for new subjects answer a different question. It is also important to generate longitudinal and event outcomes using the same posterior draw: drawing them independently can erase the dependence the joint model is intended to represent.

Posterior predictive checks complement rather than replace other diagnostics. Residual plots can focus attention on measurement-model structure, while posterior predictive checks assess whether complete replicated data sets reproduce selected features. Convergence and effective sample size diagnostics address whether the posterior simulation is reliable; they do not establish that the model can reproduce the data. Conversely, a plausible-looking predictive plot does not repair poorly mixing chains. Predictive checks also depend on which discrepancies are examined, so no finite set can establish overall model adequacy.

Common pitfall: “The observed trajectory falls inside the 95% posterior predictive envelope, so the model is correct.” An envelope only assesses a selected summary, and inclusion does not prove correctness. Examine multiple scientifically relevant features, note whether the envelope is pointwise or simultaneous, and check whether the model reproduces both longitudinal and event patterns. A conditional envelope may also understate problems in the population distribution of random effects.

Justification for Posterior Predictive Generation

The posterior predictive distribution follows directly from averaging the model's sampling distribution over uncertainty about its parameters. For a replicate \(D^{\mathrm{rep}}\), the law of total probability gives

p(D^{\mathrm{rep}}\mid D)=\int p(D^{\mathrm{rep}}\mid\theta)\,p(\theta\mid D)\,d\theta.

This identity justifies the simulation procedure: drawing \(\theta^{(s)}\) from \(p(\theta\mid D)\) and then drawing \(D^{\mathrm{rep},(s)}\) from \(p(D^{\mathrm{rep}}\mid\theta^{(s)})\) produces draws from the posterior predictive distribution. In a joint model, the conditional sampling distribution must generate both components together, including their shared random effects and the censoring mechanism. If the model conditions on observed subjects, the integral also includes their posterior random effects; for new-subject prediction, it instead includes random effects drawn from their population distribution.

Key takeaway: A posterior predictive check compares observed features with features generated by the fitted joint model, under a clearly stated predictive target and censoring design. It is most useful when the replicated data preserve the joint structure and the chosen summaries probe both biomarker trajectories and event-time patterns.

Check Your Understanding

Answer each question, then compare with the concise solutions below.

In a conditional check for the observed cohort, should you draw new random effects from the population distribution or use posterior draws of the subjects' random effects?
With administrative censoring at month 10, a simulated latent event time is 12 months. What are the replicated observed follow-up time and event indicator?
In 500 replicates, 25 event counts are at least as large as the observed count. Calculate the upper-tail posterior predictive area and interpret it.
Why can a conditional trajectory check look satisfactory even when the population random-effects distribution is misspecified?
Why should replicated longitudinal and event outcomes be generated from a coherent shared posterior draw?

Solutions. (1) Use posterior draws of the observed subjects' random effects; drawing new effects instead targets new-subject prediction. (2) The follow-up time is 10 months and the event indicator is 0, because the event occurs after censoring. (3) The area is \(25/500=0.050\); 5% of replicated event counts are at least as large as observed, a descriptive upper-tail summary rather than a formal p-value. (4) Conditioning on subject-specific effects already informed by the data can reproduce those subjects without adequately testing whether the population distribution generates realistic effects. (5) Separate draws can break the modeled dependence between biomarker trajectory and event risk, making a check unable to assess the joint structure.