Residual And Predictive Checks For Joint Models

Why Joint Models Need Joint Checks

A model for several endpoints must describe more than each endpoint’s average response. It must also describe how outcomes occur together. For example, a model may predict the observed proportion of patients with an adverse event and the observed distribution of a symptom score, yet still predict too few patients experiencing both outcomes. Endpoint-specific checks can expose marginal problems; joint predictive checks are needed to assess whether the model reproduces clinically relevant combinations.

This tutorial focuses on non-continuous endpoints, for which raw residuals such as observed minus predicted values are often discrete, bounded, or strongly dependent on the predicted mean. A residual plot that appears familiar from continuous-data analysis may therefore be difficult to interpret. The key is to compare diagnostics with what the fitted model itself predicts, including the dependence structure and the sampling design.

Endpoint-Specific Residuals

For subject \(i\) and endpoint \(j\), let \(Y_{ij}\) be the observed outcome and let \(F_{ij}(y)\) be its fitted predictive cumulative distribution function under the model and the relevant covariates. The appropriate predictive distribution depends on the diagnostic question. A population-level check generally uses predictions that integrate over random effects; a conditional check may instead condition on estimated or known subject-specific effects. State which is being used, because conditional predictions can look better simply by using information from the same subject’s outcomes.

For a discrete outcome, \(F_{ij}(Y_{ij})\) is not generally uniform when the model is correct: it has jumps at possible outcome values. A useful residual is the randomized quantile residual. Draw \(U_{ij}\) from a uniform distribution over the probability interval assigned to the observed outcome, then transform it through the standard normal quantile function.

Here \(F_{ij}(Y_{ij}^{-})\) is the probability of an outcome strictly below the observed value, and \(\Phi^{-1}\) is the standard normal quantile function. For an ordered categorical endpoint, the interval spans the probability assigned to the observed category. Under a correctly specified predictive distribution, randomized probability integral transform values are uniform, so the transformed residuals are standard normal in distribution. The randomization fills in the jumps; it does not add information to the data.

For binary outcomes, the residual therefore depends on whether the observed value is 0 or 1 and on the fitted event probability. For counts, the interval is defined by the cumulative probabilities up to and below the observed count. Residual-versus-time, residual-versus-predicted-risk, and residual-versus-exposure plots can reveal systematic patterns, changing dispersion, or a missing time trend. A histogram or normal quantile plot can also be informative, but neither should be judged against a continuous-data ideal without accounting for estimation and randomization.

Checking Cross-Endpoint Residual Association

After assessing each endpoint’s marginal fit, examine whether residuals from different endpoints remain associated. If the fitted joint model adequately represents the dependence relevant to the diagnostic, residual association should be compatible with the model’s own predictive variation. A sample correlation of endpoint-specific residuals is a useful summary for approximately continuous residual scales, but it is not a universal test of dependence. In particular, randomized residuals contain added random variation, and a correlation near zero does not rule out nonlinear or tail dependence.

A defensible comparison uses simulated reference data. Fit the model, simulate complete endpoint vectors under its fitted joint distribution using the same covariates and observation structure, and apply the same residual construction to each simulated data set. Compare the observed cross-endpoint residual correlation or another chosen association statistic with the simulated distribution. If parameters were estimated, a parametric bootstrap that refits the model in each replicate better reflects estimation uncertainty than simulation at fixed parameter estimates. A residual statistic calculated with the fitted parameters is not automatically distributed as it would be if those parameters were known.

For a more direct transformation-based check, a continuous correctly specified joint distribution can be transformed sequentially using conditional cumulative distributions. The resulting probability integral transform variables are independent uniforms under the model. Discrete outcomes require randomized conditional probability intervals; results can depend on endpoint ordering and the randomization. Thus, sequential transforms are useful diagnostics, but should be accompanied by joint predictive checks rather than treated as a definitive independence test.

Joint Predictive Checks

A joint predictive check begins by defining a discrepancy that captures a feature of the endpoint vector the model should reproduce. Examples include the number of patients with both binary events, the frequency of a high count occurring together with a binary event, or a correlation or tail-event rate. Simulate full vectors from the fitted joint distribution, preserving covariates, repeated-measure structure, and the relevant population or conditional prediction target. Recalculate the discrepancy in every simulated data set and compare its reference distribution with the observed value.

The check is most interpretable when the statistic reflects the scientific question and its direction is clear. “Both events” is different from “at least one event”; matching each marginal event rate does not guarantee that either joint frequency is correct. A central simulation interval can flag unusual values, but it is a diagnostic summary rather than a confidence interval for a model parameter. If many statistics, endpoints, subgroups, and times are checked, some will look unusual by chance, so patterns and scientific relevance matter.

Justification for Simulation-Based Reference Checks

The logic of a predictive check follows directly from the assumed joint distribution. Let \(T(\mathbf{Y})\) be a chosen statistic of all endpoints, and let \(\mathbf{Y}^{(b)}\) be a complete replicate generated under the fitted joint model. If the model and its parameters were known and correct, the observed vector and a simulated vector would have the same distribution conditional on the design. Therefore, \(T(\mathbf{Y})\) and \(T(\mathbf{Y}^{(b)})\) would have the same distribution as well. An observed statistic in the tail of simulated values signals a mismatch between the data and that model-based reference.

In practice, \(\theta\) is estimated, so this equality is only an ideal reference if the estimate is treated as fixed. Refitting simulated data can approximate the additional variability from estimation and model fitting. A Bayesian posterior predictive check instead averages predictions over a posterior distribution of parameters. These approaches answer related but not identical diagnostic questions; neither turns a check into a formal proof of adequacy.

Failure Modes, Alternatives, and Limitations

A common mistake is to calculate raw binary residuals, \(Y-\hat p\), and interpret their spread as though it should be constant across fitted risks. Their conditional variance is \(\hat p(1-\hat p)\), so the scale changes with the predicted probability. A second mistake is to conclude that a model has correct dependence because each endpoint has an acceptable residual plot. Marginal fit does not determine the joint distribution. The correction is to use endpoint-appropriate residuals and separately assess joint statistics against simulated data from the fitted model.

Other warning signs include residual trends over time or exposure, systematic underprediction in a covariate subgroup, excessive simulated-versus-observed zero counts, and joint event frequencies consistently outside predictive ranges. Discrete randomization can make a single plot appear different across random seeds; repeat the randomization or use stable summaries. Sparse joint outcomes can also leave checks with little power, so an unflagged diagnostic is not proof of adequate fit.

Alternatives depend on the diagnostic target. Endpoint-specific residual plots are easy to interpret for marginal structure, while randomized conditional transforms can probe the joint distribution more directly but can be ordering-sensitive for discrete outcomes. Joint predictive simulations are flexible and accommodate mixed endpoint types, but their conclusions depend on the chosen discrepancy and the model’s prediction target. Formal tests can quantify a prespecified departure, but may be poorly calibrated after model selection or with sparse data. All checks are complements: none replaces assessment of design, model assumptions, scientific plausibility, and parameter uncertainty.