What a Shared-Frailty Result Describes
A shared-frailty model allows members of the same cluster to have correlated event histories because they share an unobserved multiplicative risk factor. After fitting such a model, two results are easy to overstate: the estimated frailty variance and the predicted frailty for a particular cluster. The first summarizes residual heterogeneity across clusters under the fitted model; the second estimates a cluster’s relative risk given its observed history and the model assumptions. Neither is a direct measurement of a biological mechanism or a causal cluster effect.
The interpretation below uses a proportional-hazards shared-frailty model. Let \(j\) index clusters and \(i\) index individuals within cluster \(j\). Conditional on the shared frailty \(U_j\) and covariates \(X_{ij}\), the hazard for an individual is
Here \(h_0(t)\) is the baseline hazard, \(\beta\) is the vector of regression coefficients, and \(U_j\) is positive. A common convention is \(E(U_j)=1\), so the frailty scales risk relative to the model’s reference level rather than changing its overall scale. For a gamma frailty, write \(\operatorname{Var}(U_j)=\theta\). The frailties are typically assumed independent across clusters and independent of the modeled covariates; conditional on the frailties and covariates, the event processes are assumed independent across individuals. The event and censoring assumptions of the fitted survival model must also be suitable for the data.
Interpreting the Frailty Variance
A positive estimate of \(\theta\) indicates that the fitted model attributes some remaining cluster-to-cluster risk variation to a shared latent multiplier, after accounting for the included covariates and the specified baseline hazard. It does not say that a fraction \(\theta\) of events is caused by clustering, nor does it quantify the proportion of outcome variation explained in the way an \(R^2\) might in a linear model. Its numerical meaning depends on the frailty distribution and model specification.
The null hypothesis of no frailty variation is \(H_0:\theta=0\). This value lies on the boundary of the permissible parameter space \(\theta\geq 0\), so a routine Wald test based on an estimate divided by a standard error can be unreliable. Under standard regularity conditions for a single variance component, a likelihood-ratio statistic for this null has an asymptotic reference distribution that is a 50:50 mixture of a point mass at zero and a chi-square distribution with one degree of freedom. The conditions and finite-sample behavior matter; with few clusters or weak information, a parametric bootstrap or other model-appropriate calibration may be preferable.
This mixture adjustment is not a universal substitute for checking the software’s test and its assumptions. In particular, an ordinary chi-square p-value may not be the boundary-adjusted p-value. A confidence interval for \(\theta\) should respect the nonnegative parameter space; an interval that extends below zero is a warning that a symmetric Wald approximation is being used where it may be unsuitable.
Worked Example: Assessing Evidence of Residual Heterogeneity
Suppose a hypothetical analysis of 48 clinics reports a gamma frailty variance estimate \(\hat{\theta}=0.30\). The likelihood-ratio statistic comparing the fitted model with the no-frailty model is \(LR=4.20\). Assume the boundary-mixture approximation is appropriate for this illustration.
Under the mixture reference distribution, the p-value is half the upper-tail probability from a chi-square distribution with one degree of freedom:
At a prespecified significance level \(\alpha=0.05\), this is evidence against the no-frailty model, conditional on the model assumptions and the validity of the approximation. The result supports residual cluster heterogeneity as represented by the gamma frailty; it does not establish why clinics differ or prove that an unmeasured clinic characteristic caused the difference. If the ordinary chi-square tail probability \(0.0404\) had been reported instead, it would not be the boundary-mixture p-value in this example.
The estimate \(\hat{\theta}=0.30\) describes the fitted variance of the mean-one frailty distribution, not an absolute risk difference. Its interpretation should be accompanied by uncertainty, the frailty distribution, the number of independent clusters, and the model context. Statistical evidence and practical importance are separate questions.
Conditional Hazard Ratios and Marginal Risk
For two individuals with covariate vectors \(x_1\) and \(x_0\) in the same cluster, their hazards at a given time share the same \(U_j\). Their conditional hazard ratio is therefore
Thus, a coefficient-based hazard ratio compares individuals with different covariates conditional on the same frailty. It is not automatically a population-average hazard ratio. When frailty is integrated out, the distribution of frailty among those still event-free changes over time: high-frailty individuals tend to experience events earlier, leaving a more selected group at risk. Consequently, marginal hazards and hazard ratios need not retain the proportional form of the conditional model.
For example, with a single individual, gamma frailty of mean one and variance \(\theta>0\), and cumulative baseline hazard \(H_0(t)\), integrating over frailty gives marginal survival
Differentiating the negative logarithm of this survival function gives the marginal hazard
The denominator depends on time and covariates, so the marginal hazard ratio for two covariate values generally varies with time. This calculation illustrates why conditional regression coefficients from a frailty model should be labelled as conditional, not silently interpreted as population-average effects.
Worked Example: Interpreting a Conditional Hazard Ratio
In a hypothetical shared-frailty model, the fitted coefficient for a treatment indicator is \(\hat{\beta}=-0.40\). The model reports a treatment hazard ratio of \(\exp(-0.40)=0.670\). Suppose treatment and control individuals being compared belong to the same clinic, so they share the same clinic frailty.
Conditional on that clinic frailty and the other modeled covariates, the estimated treatment hazard is \(0.670\) times the control hazard at a given time under the proportional-hazards model. Equivalently, it is estimated to be 33.0% lower, since \(1-0.670=0.330\). This comparison is about the instantaneous event rate among individuals still at risk, not a 33.0% reduction in cumulative event probability.
If the same comparison is made across populations with different frailty distributions or after averaging over frailty, the marginal hazard ratio need not equal \(0.670\), and it may vary over time. The conditional estimate should therefore be reported with its conditioning interpretation and uncertainty.
Understanding Cluster-Level Frailty Predictions
A fitted model may provide cluster-specific frailty predictions, often described as empirical Bayes estimates or posterior means. These summarize how the observed event history for a cluster updates the assumed frailty distribution, while also accounting for uncertainty and typically shrinking noisy cluster estimates toward the population mean of one. A prediction above one indicates higher latent relative risk than the model’s reference level, conditional on the fitted model and that cluster’s history. It is not a directly observed cluster hazard ratio with no uncertainty.
A simple gamma model shows the source of shrinkage. Suppose, for illustration, that a cluster’s event count \(D_j\), conditional on its frailty, follows a Poisson model with mean \(U_j A_j\), where \(A_j\) is the cluster’s known cumulative exposure on the model’s hazard scale. With prior frailty \(U_j\sim\operatorname{Gamma}(1/\theta,1/\theta)\), using the shape-rate convention, the likelihood is proportional to \(U_j^{D_j}\exp(-A_jU_j)\). Multiplying this likelihood by the gamma density gives
This conjugate calculation is an illustration, not a general recipe for extracting predictions from every fitted survival model. Real event-history likelihoods, time-varying covariates, censoring, and parameter estimation can make the calculation more involved. Use the fitted model’s documented prediction method and report what summary it returns, such as a conditional mode or posterior mean.
Worked Example: A Shrunk Cluster Prediction
Consider a hypothetical gamma frailty model with \(\theta=0.25\). For one clinic, use the simplified count-and-exposure setup above, with \(D_j=6\) events and \(A_j=20\) units of cumulative hazard exposure. The prior shape and rate are both \(1/\theta=4\).
The posterior mean frailty is
In this illustrative model, the clinic’s posterior mean multiplier is about \(0.417\), below the population mean of one. The raw event count of six alone does not imply high frailty: the clinic had 20 units of cumulative exposure on the hazard scale, so its event count is not unexpectedly large relative to that exposure. The prediction is also pulled toward one by the prior information represented by the gamma distribution. In a full fitted model, uncertainty in the regression parameters and baseline hazard may contribute further uncertainty.
Diagnostics, Comparisons, and Limitations
A frailty estimate is only as informative as the model that produced it. Review the number of independent clusters and the distribution of cluster sizes: many individuals in only a few clusters do not provide the same information about between-cluster variation as many independent clusters. Check convergence, whether the variance estimate is near zero or very large, and whether results are sensitive to plausible modeling choices. Compare observed and model-predicted event patterns across clusters, and inspect whether a small number of clusters dominate the fitted heterogeneity. Cluster predictions with very wide uncertainty are weak evidence for ranking clusters.
A shared frailty model and a marginal model with robust variance answer different questions. A frailty model specifies a latent cluster-level distribution and gives conditional hazard effects; a marginal robust-variance analysis targets regression effects at the population level while adjusting uncertainty for within-cluster dependence. Robust standard errors do not, by themselves, estimate a frailty variance or provide cluster-specific latent-risk predictions. Neither approach automatically corrects misspecification of the mean model, informative cluster size, or dependence structures not represented by the analysis.
The frailty distribution is also an assumption. Gamma frailty is convenient and common, but a different distribution can yield different variance estimates or predictions. A nonzero \(\hat{\theta}\) is evidence of residual heterogeneity only relative to the included covariates, baseline hazard, event-history structure, and frailty distribution. It is not proof that the assumed distribution is true, that the clusters are exchangeable in reality, or that a particular omitted variable explains the variation.
Name the frailty distribution, its mean convention, and the meaning of its variance parameter.
For a test of zero variance, account for the boundary at zero and consider the number of independent clusters.
Describe a cluster estimate as model-based and history-informed, with shrinkage and uncertainty—not as an observed or causal effect.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- In a mean-one gamma frailty model, what does \(\theta=0.4\) represent, and what does it not represent?
- A treatment coefficient is \(-0.25\). Calculate and interpret its conditional hazard ratio.
- Why can a cluster prediction be below one even when the cluster has several observed events?
- Why is an ordinary chi-square reference distribution potentially inappropriate for a likelihood-ratio test of zero frailty variance?
Solutions. (1) It is the variance of the model’s positive, mean-one latent cluster multiplier; it is not the percentage of events due to clustering. (2) \(\exp(-0.25)=0.779\), so, conditional on the same frailty and other modeled covariates, the estimated treatment hazard is about 22.1% lower at a given time under proportional hazards. This is not a cumulative-risk reduction. (3) The prediction uses event history relative to exposure and other modeled information, and it is shrunk toward the population mean; an event count alone is not enough to determine relative risk. (4) Zero is on the boundary of the nonnegative variance parameter space, so under standard conditions the null likelihood-ratio distribution is a boundary mixture rather than the usual chi-square distribution.