Why the Frailty Distribution Matters
A frailty model represents unobserved differences in risk by multiplying a subject’s or cluster’s conditional event rate by a positive random effect. The frailty distribution is more than a mathematical convenience: it determines how individual risks mix together when frailty is unobserved. Consequently, two models with the same mean frailty and similar estimated variance can imply different marginal survival, dependence among members of a cluster, and predicted risks.
The question here is not which frailty distribution is universally correct. It is whether conclusions are robust to plausible choices of distribution while keeping the model’s other structure comparable. Interpret comparisons in light of the data, the fitting procedure, and the estimand of interest.
The comparison assumes the same event-process specification, covariate structure, baseline-hazard treatment, cluster definition, and censoring assumptions across candidate fits. In particular, clusters are independent, and event histories are conditionally independent across members given the shared frailty and modeled covariates, as required by the chosen model. Censoring must be noninformative conditional on the modeled information, or its mechanism must be handled in the likelihood. A candidate frailty distribution must be positive and properly normalized to mean one; its variance parameter must be defined consistently. If these assumptions are inadequate, a distribution comparison cannot repair the model.
The gamma, lognormal, and inverse Gaussian choices below all have mean one. For \(\theta>0\), their parameterizations and transforms are:
| Frailty distribution | Mean-one parameterization | Laplace transform \(L(s)\) |
|---|---|---|
| Gamma | Shape \(1/\theta\), scale \(\theta\) | \((1+\theta s)^{-1/\theta}\) |
| Lognormal | \(\log Z\sim N(-\sigma^2/2,\sigma^2)\), with \(\theta=\exp(\sigma^2)-1\) | \(E[\exp\{-s\exp(U)\}]\), where \(U\sim N(-\sigma^2/2,\sigma^2)\); typically evaluated numerically |
| Inverse Gaussian | Mean \(1\), shape parameter \(1/\theta\) | \(\exp\{(1/\theta)[1-\sqrt{1+2\theta s}]\}\) |
At \(\theta=0\), each model is interpreted by its limiting no-frailty case \(Z=1\). For the lognormal distribution, matching a variance \(\theta\) requires \(\sigma^2=\log(1+\theta)\). All three distributions therefore have the same mean and variance when compared at the same \(\theta\), but they need not have the same skewness or tail behavior.
How Distributional Shape Changes Marginal Quantities
Conditional survival through cumulative exposure \(H\) is \(\exp(-ZH)\). Averaging over the frailty distribution gives marginal survival \(S_m(H)=L(H)\); therefore the marginal event risk by that time is \(1-L(H)\). The transform, not just the first two moments, determines these quantities.
The marginal hazard also changes with time, even when the conditional hazard is proportional. Differentiating marginal survival with respect to \(H\) gives the marginal hazard as the conditional baseline-and-covariate multiplier times \(-L'(H)/L(H)\). This ratio depends on the entire frailty distribution. Thus conditional coefficients, marginal hazards, and predicted risks should not be treated as interchangeable summaries.
For two members of the same cluster, let \(Y_j(t)\) indicate that member \(j\) has had no event by time \(t\). Conditional on their shared \(Z\), their survival indicators are independent under the model. Their marginal covariance is \(L(H_1+H_2)-L(H_1)L(H_2)\), where \(H_1,H_2\) are their respective cumulative exposures. A different frailty distribution can therefore change the induced dependence even when the mean and variance match.
For \(H_1,H_2>0\), this covariance is nonnegative: as \(Z\) increases, each conditional survival probability \(\exp(-ZH_j)\) decreases, so the two conditional survival probabilities vary in the same direction across frailty values. More generally, the exact dependence and its magnitude are governed by the transform. At larger cumulative exposure, the factor \(\exp(-HZ)\) suppresses contributions from large \(Z\), so marginal survival can become especially sensitive to the distribution’s lower tail, including behavior near zero. Differences in distribution shape can matter even when mean and variance agree.
Worked Comparisons
Worked Example: Same Mean and Variance, Different Dependence
Suppose two members share a mean-one frailty with variance \(\theta=0.5\), and each has cumulative exposure \(H=1\). Compare gamma and inverse Gaussian frailties. The common model assumptions are that the cluster shares one frailty and that the two members’ survival outcomes are conditionally independent given it.
For the gamma model, \(L(1)=(1+0.5)^{-2}=1/2.25=0.4444\), and \(L(2)=(1+0.5\times2)^{-2}=2^{-2}=0.2500\). Thus the covariance of the two no-event indicators is \(0.2500-(0.4444)^2\approx0.05251\) (using unrounded values, \(0.25-(4/9)^2=0.05247\)).
For the inverse Gaussian model, \(L(1)=\exp\{2[1-\sqrt{2}]\}\approx0.4367\), and \(L(2)=\exp\{2[1-\sqrt{3}]\}\approx0.2313\). Using unrounded transform values, the covariance is approximately \(0.04055\). Using the displayed rounded values gives \(0.2313-(0.4367)^2\approx0.04059\); the small difference is due to rounding. Both covariances are positive, but the implied dependence differs despite identical frailty means and variances.
Worked Example: Comparing Three Hypothetical Fits
Consider hypothetical fits to the same clustered recurrent-event data. Each uses the same covariates and baseline-hazard specification, and each has five estimated parameters in the likelihood criterion. Suppose the maximized log likelihoods are \(-412.6\) for gamma, \(-410.9\) for lognormal, and \(-411.7\) for inverse Gaussian. The fitted frailty variances are \(0.46\), \(0.61\), and \(0.40\), respectively.
For each model, \(\mathrm{AIC}=-2\ell+2k\). Hence the gamma AIC is \(-2(-412.6)+2(5)=835.2\); the lognormal AIC is \(-2(-410.9)+10=831.8\); and the inverse Gaussian AIC is \(-2(-411.7)+10=833.4\). Among these candidates, lognormal has the lowest AIC, by 3.4 relative to gamma and 1.6 relative to inverse Gaussian. This is relative evidence under the specified candidate set, not proof that the lognormal mixing law is true.
Suppose each fitted model also reports a two-year event-risk prediction for the same covariate profile: 0.286 under gamma, 0.301 under lognormal, and 0.279 under inverse Gaussian. These are hypothetical model outputs, with each model’s own fitted parameters and baseline. The lognormal prediction is 0.015 higher than the gamma prediction, or 1.5 percentage points. Whether this difference matters depends on the decision context and its uncertainty. The variance estimates alone would not reveal this prediction spread.
Worked Example: Translating a Transform Into Risk
At a specified covariate profile and time, suppose the fitted conditional cumulative exposure is \(H=1\), and compare gamma and inverse Gaussian frailties with \(\theta=0.5\). For gamma frailty, the marginal event risk is \(1-L(1)=1-(1+0.5)^{-2}=1-0.4444=0.5556\).
For inverse Gaussian frailty, it is \(1-\exp\{2[1-\sqrt{2}]\}\approx1-0.4367=0.5633\). The inverse Gaussian model gives an event risk about \(0.5633-0.5556=0.0077\), or 0.77 percentage points, higher at this exposure. The comparison holds \(H\) and \(\theta\) fixed to isolate distributional shape; fitted model comparisons generally also reflect changes in estimated baseline and other parameters.
Fitting, Diagnostics, and Practical Comparisons
A sensitivity analysis should change the frailty distribution while holding the rest of the analysis as comparable as possible. Fit each candidate with the same observations, covariates, time scale, censoring treatment, baseline-hazard flexibility, and cluster structure. Use the appropriate likelihood for each distribution; numerical integration or adaptive quadrature must be sufficiently accurate. Check convergence, gradients, parameter boundaries, and stability from reasonable starting values before interpreting a model comparison.
Compare more than one output. Useful summaries include the fitted frailty variance and its uncertainty, likelihood-based fit criteria where applicable, conditional covariate effects, marginal survival or event risks at prespecified times, and the dependence measures relevant to the scientific question. A variance estimate is standardized across the mean-one models, but it does not fully summarize their different shapes. If a model’s predicted risks shift materially while its variance changes little, the transform’s shape is influencing the scientific conclusion.
Information criteria such as AIC are useful for comparing candidate likelihood fits to the same data and response structure. They penalize fit for the number of estimated parameters, but do not measure absolute calibration and do not establish that the selected distribution is correct. Likelihood-ratio reference distributions also require care when comparing models that are not nested; do not apply a standard nested-model test merely because two distributions are alternative candidates. Where uncertainty in a difference in predictions matters, report intervals or use a suitable resampling or model-based uncertainty procedure that respects the clustering and refitting.
Use the same data, covariates, event structure, baseline specification, and censoring assumptions for each frailty distribution.
Check numerical integration, convergence, boundary estimates, and whether the estimated likelihood is stable to starting values and quadrature settings.
Assess fit criteria and uncertainty, then compare estimated heterogeneity, induced dependence, and risks at scientifically meaningful times and covariate profiles.
State which conclusions remain stable, which change, and what the candidate distributions cannot establish about the true frailty law.
Residual and calibration diagnostics can reveal lack of fit, but they may have limited power to distinguish frailty distributions, particularly with few clusters or little follow-up. Check calibration over time and across risk groups, and inspect whether predictions are systematically too high or too low. Warnings include unstable variance estimates, large changes with quadrature settings, a variance estimate near zero, sparse cluster information, or important changes in prediction under plausible distribution choices. A flexible baseline or omitted covariates can also masquerade as frailty-distribution sensitivity.
Interpretation and Common Pitfall
The gamma model has a convenient closed-form transform; the lognormal transform generally needs numerical evaluation; and the inverse Gaussian has a different transform and tail shape. These are computational and distributional distinctions, not a ranking from simple to correct. Each is a parametric assumption about unobserved heterogeneity. Stronger tail behavior can affect some predictions and dependence summaries, but the direction and size of change depend on the cumulative exposure and fitted parameters.
A common mistake is to conclude that similar estimated frailty variances imply equivalent fitted models. That reasoning overlooks the fact that \(L(H)\) is not determined by mean and variance alone. The correction is to compare the quantities used for inference and prediction under each fitted model, with uncertainty, and to report any consequential sensitivity. Conversely, choosing the model with the smallest AIC does not remove uncertainty about the frailty distribution or validate the rest of the model assumptions.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- Why can two mean-one frailty distributions with the same variance imply different marginal survival?
- For gamma frailty with \(\theta=0.5\) and \(H=1\), calculate marginal event risk.
- In a shared-frailty model, what does positive covariance between two members’ no-event indicators represent?
- Why is the lowest AIC not proof that a particular frailty distribution is true?
Solutions. (1) Marginal survival is \(L(H)=E\{\exp(-HZ)\}\), which depends on the distribution’s full shape, not only its first two moments. (2) \(L(1)=(1+0.5)^{-2}=0.4444\), so risk is \(1-0.4444=0.5556\), or 55.56% at the specified exposure. (3) Members in clusters with smaller shared frailty tend to have higher probabilities of no event, while those with larger frailty tend to have lower probabilities; this shared variation induces positive marginal association. (4) AIC ranks relative expected information loss among the specified fitted candidates, under its assumptions and approximations. It does not establish absolute fit or show that the candidate set contains the true frailty distribution.