Posterior Predictive Criteria for Joint Model Comparison

Why the Comparison Unit Matters

Bayesian information criteria summarize different aspects of a fitted model. DIC penalizes a model’s deviance according to an estimated effective number of parameters; WAIC and leave-one-out (LOO) criteria target predictive performance. In a joint longitudinal and time-to-event model, however, a score is only meaningful after specifying what is being predicted. Are we predicting the full record of a new subject, another measurement for an existing subject, or an event outcome conditional on that subject’s observed longitudinal history? These are different prediction tasks and can require different likelihood contributions.

For a new-subject prediction, the natural comparison unit is usually the subject: that subject’s longitudinal observations and event or censoring information are a dependent record linked through shared random effects. Treating every measurement and event record as independent comparison units can change the question and produce misleadingly precise comparisons.

Definition: Let $y_i$ denote subject $i$’s complete observed record, including their longitudinal measurements and their event-time or censoring data. Let $\psi$ denote population-level parameters and $b_i$ the subject-specific random effects. The marginal subject contribution is $$ p(y_i\mid\psi)=\int p(y_i\mid b_i,\psi)\,p(b_i\mid\psi)\,db_i. $$ For a shared-random-effects joint model, the conditional record likelihood inside this integral is the joint likelihood of all components for subject $i$, not a product of separately marginalized components. This contribution is appropriate when the prediction target is a new subject’s full record.

The integral expresses the dependence structure: conditional on $b_i$ and $\psi$, the longitudinal and survival components may have a specified factorization, but they share $b_i$. The integral over $b_i$ must therefore be applied to their joint conditional likelihood. For prediction about an already observed subject, conditioning on information from that subject may be appropriate instead. State the prediction target before selecting the factorization.

DIC: Deviance and Effective Complexity

For observed data $y$, define the deviance, up to an additive constant that is the same across candidate models, as $D(\theta)=-2\log p(y\mid\theta)$, where $\theta$ denotes the parameters used in that likelihood. DIC combines the posterior mean deviance with an effective parameter count. In the common posterior-mean form, $\bar D$ is the posterior mean of $D(\theta)$, $\hat\theta$ is the posterior mean of $\theta$, and $D(\hat\theta)$ is the deviance evaluated at that mean.

p_D=\bar D-D(\hat\theta), \qquad \operatorname{DIC}=\bar D+p_D=2\bar D-D(\hat\theta).

Smaller DIC indicates a preferred trade-off between fit and effective complexity among the models being compared, provided the same data and likelihood convention are used. Unlike a conventional parameter count, $p_D$ is an estimate that can be small, negative, or unstable in difficult posterior distributions. In a hierarchical joint model, a DIC based on the likelihood conditional on subject random effects is not automatically comparable to one based on the marginal subject likelihood. These choices represent different treatments of latent effects and must be reported.

Worked Example: Calculating DIC from Posterior Deviances

Suppose four draws from a hypothetical fitted joint model have deviances $118,122,120,$ and $124$. Their average is $$ \bar D=\frac{118+122+120+124}{4}=\frac{484}{4}=121. $$ Suppose the deviance evaluated at the posterior mean of the model parameters is $D(\hat\theta)=116$. Then $$ p_D=121-116=5,\qquad \operatorname{DIC}=121+5=126. $$ Equivalently, $2(121)-116=126$. If a second model fitted to the same observations has DIC $130$ under the same likelihood convention, the first model has the lower DIC by 4 points. This is a relative criterion, not a test or a measure of the probability that either model is true.

WAIC: Pointwise Predictive Accuracy

WAIC begins by factoring the data likelihood into contributions $p(y_i\mid\theta)$ for the chosen units $i$. For a new-subject target, use the marginal subject contributions $p(y_i\mid\psi)$, with latent random effects integrated out. Given posterior draws $\theta^{(1)},\ldots,\theta^{(S)}$, the log pointwise predictive density (lppd) and WAIC effective complexity $p_{\mathrm{WAIC}}$ are

\operatorname{lppd} =\sum_{i=1}^{n}\log\left\{\frac{1}{S} \sum_{s=1}^{S}p(y_i\mid\theta^{(s)})\right\}, \qquad p_{\mathrm{WAIC}} =\sum_{i=1}^{n}\operatorname{Var}_{s} \left\{\log p(y_i\mid\theta^{(s)})\right\},

where the variance is the sample variance across posterior draws. The information criterion is $\operatorname{WAIC}=-2(\operatorname{lppd}-p_{\mathrm{WAIC}})$; smaller values indicate better estimated predictive performance. Equivalently, the estimated expected log predictive density is $\operatorname{elpd}_{\mathrm{WAIC}}=\operatorname{lppd}-p_{\mathrm{WAIC}}$, for which larger is better. Calculate from unrounded terms and round only the final reported result.

Worked Example: Computing Subject-Level WAIC

Consider three posterior draws of the log marginal likelihood contributions for each of two hypothetical subjects. The values are $(-2.0,-2.4,-1.8)$ for subject 1 and $(-1.5,-1.7,-1.4)$ for subject 2.

For subject 1, the log of the average likelihood is $$ \log\left\{\frac{e^{-2.0}+e^{-2.4}+e^{-1.8}}{3}\right\} \approx -2.037. $$ For subject 2, it is $$ \log\left\{\frac{e^{-1.5}+e^{-1.7}+e^{-1.4}}{3}\right\} \approx -1.526. $$ Using the original values without rounding the subject contributions first gives $\operatorname{lppd}\approx-3.562453$, or $-3.562$ to three decimal places.

For subject 1 the mean log likelihood is $-2.0667$, and its sample variance is approximately $0.09333$. For subject 2 the mean is $-1.5333$, and its sample variance is approximately $0.02333$. Thus $p_{\mathrm{WAIC}}\approx0.116667$. Substituting the unrounded quantities gives $$ \operatorname{WAIC} =-2(-3.562453-0.116667) \approx 7.358239, $$ or $7.358$ to three decimal places. This score is on the negative-twice-log-predictive-density scale; by itself, it does not establish that the model is adequate.

Leave-One-Out Prediction

In leave-one-out cross-validation, each unit is held out in turn. For subject $i$, the target predictive density is $$ p(y_i\mid y_{-i})=\int p(y_i\mid\theta)\,p(\theta\mid y_{-i})\,d\theta, $$ where $y_{-i}$ denotes the records for all other subjects. Summing the log predictive densities gives $\operatorname{elpd}_{\mathrm{LOO}}=\sum_i\log p(y_i\mid y_{-i})$; larger values indicate better predictive performance. The corresponding information-criterion scale is $\operatorname{LOOIC}=-2\operatorname{elpd}_{\mathrm{LOO}}$, so smaller is better. WAIC and LOO are closely related large-sample predictive criteria, but they are not identical calculations.

Exact refitting after removing each subject can be costly. Pareto-smoothed importance sampling LOO (PSIS-LOO) approximates these held-out predictions using draws from the full-data posterior and diagnoses the stability of that approximation. Its Pareto $k$ diagnostic signals influential observations or unreliable importance sampling; values above roughly $0.7$ warrant concern, and values near or above $1$ indicate especially serious instability. These are practical warning thresholds, not guarantees. When diagnostics are poor, refitting for the affected units or using another validation design may be necessary.

Worked Example: Comparing Two Models by LOO

Suppose three held-out hypothetical subjects have LOO predictive densities $0.08,0.12,0.05$ under model A and $0.10,0.10,0.06$ under model B. The sums of log predictive densities are $$ \operatorname{elpd}_{A} =\log(0.08)+\log(0.12)+\log(0.05) \approx -7.642, $$ $$ \operatorname{elpd}_{B} =\log(0.10)+\log(0.10)+\log(0.06) \approx -7.419. $$ Model B’s estimated advantage is $\Delta= -7.419-(-7.642)=0.223$ log-score units. The per-subject differences, B minus A, are approximately $0.2231,-0.1823,0.1823$. Their sample standard deviation is about $0.2233$, giving an approximate standard error for the summed difference of $$ SE(\Delta)=\sqrt{3}(0.2233)\approx0.387. $$ The estimated advantage is small relative to its uncertainty, so these three units do not provide clear evidence that B predicts better. This uncertainty calculation is only illustrative with such a small number of units; practical comparisons should report the unit-level differences and their uncertainty, not just the two totals.

Shared Effects, Diagnostics, and Model Choice

The consequence of shared random effects is visible even in a simple illustration. Suppose a random effect has two possible values, each with prior probability $0.5$. For a particular subject, suppose the conditional likelihood of the observed longitudinal data is $0.9$ at the first value and $0.1$ at the second; the conditional likelihood of the observed survival data is $0.1$ and $0.9$, respectively. The marginal joint likelihood is $$ 0.5(0.9\times0.1)+0.5(0.1\times0.9)=0.09. $$ But multiplying the separately marginalized likelihoods gives $$ \{0.5(0.9)+0.5(0.1)\} \{0.5(0.1)+0.5(0.9)\}=0.5\times0.5=0.25, $$ which is not the joint marginal likelihood. Marginalizing each component separately discards their dependence through the shared effect. This matters directly when constructing subject-level likelihood contributions for WAIC or LOO.

A conditional likelihood that includes the held-out subject’s estimated random effect can also use information from that subject’s own record. It may suit a target such as predicting a later outcome for an existing subject conditional on earlier measurements, if the validation scheme is designed accordingly. It is generally not the same as predicting a new subject’s entire record. Grouping each subject’s complete record into one holdout prevents the subject’s remaining measurements from leaking into prediction of that subject’s other held-out components.

Specify the target.
State whether prediction concerns a new subject’s full record, a future measurement for an observed subject, or another clearly defined unit.

Construct the matching likelihood contribution.
For a new subject, combine their longitudinal and survival contributions conditional on their shared random effects, then integrate over those effects.

Calculate and diagnose the criterion.
Use the same units and observations across candidate models. Inspect posterior computation and, for PSIS-LOO, Pareto $k$; do not interpret an unstable score as a reliable ranking.

Report uncertainty and context.
Give the score difference, its uncertainty when available, the prediction target, and any influential units or limitations.

DIC is computationally simple but can behave poorly with skewed, multimodal, or weakly identified posteriors, and its effective parameter count depends on how latent effects enter the deviance. WAIC uses the full posterior distribution and is pointwise predictive, but its approximation can be unstable for highly influential units. LOO more directly mimics held-out prediction, although exact refitting can be expensive and importance-sampling approximations can fail. None of these criteria replaces substantive model assessment, posterior predictive checks, or careful sensitivity analysis. They compare predictive performance under a specified scoring unit; they do not show that a model is scientifically correct or that a coefficient is causal.

Common pitfall: Counting every longitudinal measurement and the event record as independent WAIC or LOO units, then interpreting the resulting score as new-subject prediction, ignores the dependence induced by shared random effects. Instead, define the prediction target first and use a likelihood factorization that respects it. When the target is a new subject, hold out and score the subject’s joint record.

Check Your Understanding

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

A model has $\bar D=84$ and $D(\hat\theta)=78$. What are $p_D$ and DIC?
For new-subject prediction in a shared-random-effects joint model, should the longitudinal and survival likelihoods be integrated over the random effects separately or jointly?
Model A has $\operatorname{elpd}_{\mathrm{LOO}}=-50$, and model B has $-47$. Which has the better estimated predictive score, and by how much?
Why can a conditional likelihood using a held-out subject’s random effect give an overly optimistic assessment of prediction for a new subject?
What should be considered if a PSIS-LOO fit has a Pareto $k$ value above roughly $0.7$?

Solutions. (1) $p_D=84-78=6$, so $\operatorname{DIC}=84+6=90$. (2) Integrate the joint conditional likelihood jointly over the shared effects; separately integrating components generally changes the joint probability. (3) Model B has the larger elpd, with an estimated advantage of $3$ log-score units; uncertainty is needed to judge how compelling that difference is. (4) The held-out subject’s data may inform their estimated random effect, so the assessment can use information unavailable for a new subject. (5) Treat the importance-sampling approximation as potentially unreliable; examine influential units and consider exact refitting or another validation approach.