Frailty Model Diagnostics Using Residuals and Calibration

Why Frailty Models Need More Than a Convergence Check

A frailty model can converge and still describe the data poorly. Its baseline hazard may have the wrong time pattern, its frailty distribution may fail to capture cluster differences, or its predicted risks may be systematically too high or too low. Diagnostics ask whether these failures leave recognizable discrepancies between observed and model-implied event experience.

For recurrent events, observations within a cluster are dependent. A cluster with unusually many events can reflect a high latent frailty, a misspecified baseline hazard, or chance variation. A useful diagnostic therefore examines both when events occur and how event experience varies across clusters. Calibration adds a separate question: do clusters assigned a given predicted risk experience events at about that rate?

Model and Residual Definitions

Let \(N_{ij}(t)\) count recurrent events for individual \(j\) in cluster \(i\) by time \(t\), and let \(Y_{ij}(t)\) indicate that the individual is under observation and at risk just before \(t\). A shared-frailty proportional hazards intensity can be written

Here \(Z_i>0\) is shared by members of cluster \(i\), \(\mathcal{F}_{t-}\) is the observed history just before time \(t\), and the frailty is conventionally scaled so \(E(Z_i)=1\). The model assumes independent clusters, a correctly specified conditional intensity, predictable covariates and risk indicators, and censoring that is non-informative conditional on the model variables. The frailty distribution must also be specified for likelihood-based fitting and marginal prediction. Diagnostics assess these assumptions; they do not make them true.

With known parameters and the actual frailty, the compensator is the accumulated conditional intensity. Its difference from the event count is a martingale residual:

A residual based on a frailty estimate learned from the very same events it is meant to diagnose can be misleading: a high observed count tends to produce a high empirical Bayes frailty estimate, which then raises the fitted expected count. When feasible, calculate residuals for a later validation period using a frailty prediction based only on earlier history, or use leave-one-period-out predictions. For a new cluster, use the population frailty distribution rather than an estimate informed by that cluster's outcomes.

Why Cluster-Level Variation Matters

Frailty induces extra variation even when the conditional model is correct. For illustration, suppose that over a fixed interval the conditional event count \(N_i\), given \(Z_i\), is Poisson with mean \(Z_i\mu_i\), and \(Z_i\) has mean 1 and variance \(\theta\). Then

The second variance term is cluster-level variation. Thus, a set of cluster counts should not be judged as though all event counts were independent Poisson observations with variance \(\mu_i\). This calculation is exact under the stated conditional Poisson model; more general recurrent-event processes need not have that count distribution, but the dependence warning remains.

Residual Patterns and Deviance-Style Summaries

A single residual is rarely decisive. Plot cumulative residuals against time, and inspect cluster-level sums against fitted expected counts, covariates, cluster size, and follow-up. Under a well-specified intensity, the conditional martingale residual has mean zero over repeated realizations when the parameters and frailty are known. A systematic time trend in residual accumulation suggests that the model predicts event timing poorly. It may indicate a baseline-hazard problem, but time-varying covariate effects, changing risk-set composition, or informative censoring can create similar patterns.

For a simple count summary with observed count \(O_i\) and expected count \(E_i>0\), a signed Poisson-deviance-style transformation is

where the term \(O_i\log(O_i/E_i)\) is defined as zero when \(O_i=0\). This is the deviance residual for a Poisson count likelihood. For recurrent-event frailty data, it can be used as a descriptive transformation of observed and expected counts, but it is not generally an exact deviance residual from the full event-process likelihood. Its usual Poisson reference behavior also does not account for estimated frailties, parameter fitting, or extra cluster variation.

A useful qualitative distinction is the shape of the discrepancy. If many clusters show residual accumulation in the same time interval, the baseline hazard or a common time effect deserves scrutiny. If some clusters remain consistently above expectation while others remain below it across periods, the model may not capture the degree or structure of cluster heterogeneity. This distinction is suggestive rather than definitive: a wrongly specified baseline can interact with unequal follow-up, and an incorrectly specified frailty distribution can alter both the size and timing of residuals.

Calibration of Predicted Cluster Risk

Calibration compares predicted event risk with observed event frequency. Define the prediction target first: for example, the probability that a new cluster has at least one event in a specified period, or the expected number of events in that period. Predictions for a new cluster integrate over the population frailty distribution. Predictions for an observed cluster can use its earlier event history to update frailty, but must not use the validation outcome being assessed.

For a new cluster in a first-event setting, suppose the cumulative baseline-and-covariate hazard over the period is \(H\), and conditional survival given frailty is \(\exp(-Z H)\). The marginal risk of at least one event is \(1-E\{\exp(-ZH)\}\). For gamma frailty with mean 1 and variance \(\theta\), its Laplace transform gives

This formula assumes the stated gamma frailty model and the conditional hazard specification. It is not a generic risk formula for every frailty distribution. With recurrent events, the probability of at least one event can use the same expression when the conditional count process has probability of no events equal to \(\exp(-ZH)\); otherwise the event-process model must determine that probability.

Calibration can be summarized by comparing observed and predicted risk across groups of clusters ordered by predicted risk, or by comparing total observed events with the sum of predicted counts. Calibration-in-the-large asks whether predictions are systematically too high or too low overall; calibration across risk groups asks whether the model separates lower- and higher-risk clusters appropriately. Use cluster-level resampling or other methods that preserve dependence when quantifying uncertainty. A validation set with few clusters or few events can give a noisy calibration assessment, even if it contains many individuals within each cluster.

Separating Baseline-Hazard Problems from Cluster Heterogeneity

These are diagnostic clues, not a way to identify a unique cause from one plot. In particular, a flexible baseline hazard cannot repair unmodeled dependence, and a frailty term cannot repair a systematic time trend in the conditional event rate. Review risk-set definitions, censoring, covariate histories, and the prediction horizon before attributing a pattern to either component.

The main alternatives answer different questions. Martingale residuals emphasize signed discrepancies in event counts and their timing. Deviance-style residuals make count discrepancies more comparable in magnitude, but their reference interpretation is limited under frailty. Calibration directly assesses predicted risks or counts, usually on held-out clusters, but may conceal compensating errors: overprediction in one subgroup can offset underprediction in another. Use these tools together rather than treating one as a universal goodness-of-fit test.

Common Pitfalls and Practical Workflow

A frequent mistake is to compute an empirical Bayes frailty from all events in a cluster, use it to predict the same events, and conclude that small residuals establish good fit. The correction is to make the information boundary explicit: use earlier cluster history to predict a later period, or evaluate entirely new clusters using population-integrated predictions. Another mistake is to treat each individual or event as independent when calculating uncertainty. Shared frailty creates within-cluster dependence, so validation and uncertainty summaries should respect the cluster as the independent sampling unit.

Diagnostics have limits. Residuals depend on estimated parameters and, for observed clusters, on how frailty is predicted. Calibration can be unstable with few independent clusters, while good aggregate calibration does not guarantee correct individual predictions. No residual plot can establish that the frailty distribution is correct. The strongest assessment combines temporal patterns, cluster-level discrepancies, and out-of-sample calibration, then reports the population and prediction horizon to which each conclusion applies.