Shared Frailty

Shared Frailty in Clustered Survival Analysis

Patients treated at the same hospital tend to have more similar survival times than patients drawn at random, because they share influences that no recorded covariate captures. A shared frailty model represents those influences by multiplying the hazard of every member of a cluster by one unobserved factor, so a cluster with a high value carries a higher risk for all of its members. The variance of that factor measures how much clusters differ after case mix is allowed for, and the same variance governs how strongly event times within a cluster are correlated. Health economists meet the method in multicentre trials and in recurrent event data, where the cluster is one patient and the events are, for example, repeated serious infections or readmissions. Its most practical consequence for decision modelling is the gap between the hazard ratio for patients in the same cluster and the hazard ratio seen across a whole population.

Why clustered survival data need a shared term

A standard Cox proportional hazards model treats every event time as independent once covariates are included. Balan and Putter note that this assumption might not hold when observations come from the same centre. In their EORTC 10854 breast cancer illustration, 2687 patients from 15 centres, ranging in size between 6 and 880 patients, had Kaplan-Meier curves by centre that were hard to compare, partly because of case mix and partly because chance is difficult to disentangle from systematic differences between centres.

Two broad responses exist. A marginal model keeps the ordinary Cox fit and corrects its standard errors for clustering with a sandwich estimator, in the spirit of generalized estimating equations. A shared frailty model instead writes the clustering into the hazard, which lets the differences between clusters be quantified and gives covariate effects an interpretation for individuals within a cluster. The underlying idea comes from demography, where Vaupel, Manton and Stallard showed that ignoring heterogeneity in individual frailty distorts life table estimates, including potential gains in life expectancy from health interventions.

The shared frailty hazard and its gamma form

The model gives each cluster one frailty value and treats event times within the cluster as independent once that value is known. All dependence between members therefore runs through the shared factor.

$$h_{ij}(t \mid Z_i) = Z_i \, h_0(t) \exp\left(\beta^{\mathrm{T}} x_{ij}\right)$$

where $h_{ij}(t \mid Z_i)$ is the hazard at time $t$ for individual $j$ in cluster $i$, $Z_i$ is the positive frailty shared by every member of cluster $i$, $h_0(t)$ is the baseline hazard, $x_{ij}$ is the vector of observed covariates and $\beta$ is the vector of regression coefficients. The clusters are treated as a random sample from a population of clusters, and by convention the frailty has mean 1.

The gamma distribution is the most common choice, being the simplest and best understood, and averaging over the frailty gives the marginal, or population-averaged, survival function.

$$\bar{S}(t \mid x) = E\left[\exp\{-Z H_0(t) e^{\beta^{\mathrm{T}} x}\}\right] = \left(1 + \xi \, H_0(t) \, e^{\beta^{\mathrm{T}} x}\right)^{-1/\xi}$$

where $\bar{S}(t \mid x)$ is survival averaged over the frailty distribution for a person with covariates $x$, $H_0(t)$ is the baseline cumulative hazard and $\xi$ is the variance of the gamma frailty. For two members of one cluster the joint survival is the same expectation taken with the sum of their two cumulative hazards in place of one, which is what makes their event times positively dependent. Shared frailty can only produce positive dependence, a restriction that correlated frailty models relax.

What the frailty variance tells an analyst

The variance $\xi$ is the headline result. A value of 0 returns the ordinary Cox model, and larger values mean wider differences between clusters and stronger correlation within them. In the EORTC illustration, a gamma frailty model estimated a centre frailty variance of 0.122, with a profile likelihood 95% confidence interval of 0.048 to 0.324, and the likelihood ratio test still judged the adjusted differences in mortality between centres statistically significant. A log-normal frailty fitted to the same data gave a slightly larger variance (0.132) and almost identical regression coefficients.

Testing whether the variance is zero needs care because zero lies on the boundary of the parameter space. Balan and Putter describe a likelihood ratio test based on an equal mixture of chi-squared distributions with 0 and 1 degree of freedom, and suggest likelihood-based confidence intervals because symmetric intervals for a positive variance are not very meaningful.

The size of the variance depends on which covariates are in the model. In the chronic granulomatous disease trial of recurrent serious infections, a gamma frailty with only sex and treatment as covariates gave a variance of 0.821 in Balan and Putter's frailtyEM software paper, while their tutorial, adding age, pattern of inheritance and corticosteroid use, reported 0.555 (95% confidence interval 0.067 to 1.449). Observed covariates explain part of the heterogeneity, and the frailty absorbs the rest.

Estimation and the choice of frailty distribution

With a parametric baseline hazard the marginal likelihood can be maximised directly by maximum likelihood estimation. With an unspecified Cox baseline, Balan and Putter describe the expectation maximisation algorithm and the penalised likelihood method, the latter implemented in the coxph function of the R survival package for gamma and log-normal frailties. The expectation step yields an empirical Bayes estimate of each cluster's frailty, which for the gamma case has a simple form.

$$\hat{z}_i = \frac{1/\hat{\xi} + N_i}{1/\hat{\xi} + \hat{E}_i} = w_i \frac{N_i}{\hat{E}_i} + (1 - w_i), \qquad w_i = \frac{\hat{E}_i}{1/\hat{\xi} + \hat{E}_i}$$

where $\hat{z}_i$ is the estimated frailty of cluster $i$, $N_i$ is its observed number of events, $\hat{E}_i$ is its expected number of events given case mix and follow-up (the sum of members' fitted cumulative hazards), $\hat{\xi}$ is the estimated frailty variance and $w_i$ is the weight given to the cluster's own observed to expected ratio. Small clusters have small $\hat{E}_i$ and are pulled further towards 1.

The distribution matters more than its convenience suggests. Hougaard notes that the gamma assumption implies dependence that is most important for late events, and Balan and Putter describe the positive stable frailty as inducing early dependence. In multicentre trials, Glidden and Vittinghoff found substantial advantages for the gamma frailty model, which appeared minimally affected when the true frailty was not gamma. Balan and Putter add that choosing between frailty distributions is notoriously difficult and that the choice is usually made by convenience.

Worked example with hospital frailties and a decision model

All numbers in this example are illustrative. A gamma frailty model for mortality after a surgical procedure has been fitted across many hospitals, with an estimated frailty variance of 0.25, so the inverse variance $\hat{\xi}^{-1}$ equals 4.

1. A small hospital. Hospital A had 3 deaths against 1.5 expected, a raw observed to expected ratio of 2.00. Its frailty estimate is (4 + 3) / (4 + 1.5), or 7 divided by 5.5, which is 1.27, with weight 1.5 / 5.5 = 0.27 on its own ratio. As a check, 0.27 × 2.00 + 0.73 = 1.27.

2. A large hospital. Hospital B had 40 deaths against 32 expected, a raw ratio of 1.25. Its frailty estimate is (4 + 40) / (4 + 32), or 44 divided by 36, which is 1.22, with weight 32 / 36 = 0.89. Hospital A's apparent doubling of risk shrinks to 27% above average, while Hospital B's estimate barely moves.

A second illustrative model is fitted to a multicentre trial with a shared gamma frailty variance of 0.5, a baseline hazard of 0.2 deaths per patient-year at frailty 1, so $H_0(t) = 0.2t$, and a conditional hazard ratio of 0.6 for a new treatment against standard care.

3. Marginal survival. Under standard care, survival at 5 years is (1 + 0.5 × 1.0) to the power minus 2, which is 1 / 2.25 = 0.444, and at 10 years (1 + 0.5 × 2.0) to the power minus 2 = 0.25. Under the new treatment the cumulative hazard is multiplied by 0.6, giving 1 / 1.69 = 0.592 at 5 years and 1 / 2.56 = 0.391 at 10 years.

4. Marginal hazard ratio. The ratio of marginal hazards is 0.6 × (1 + 0.5 × 0.2t) / (1 + 0.5 × 0.6 × 0.2t). It is 0.60 at the start, 0.6 × 1.5 / 1.3 = 0.692 at 5 years and 0.6 × 2.0 / 1.6 = 0.75 at 10 years. The population-level effect weakens over time because frail patients leave the standard care risk set faster.

5. Restricted mean survival. For survival of the form (1 + ct) to the power minus 2, the area up to $T$ equals $T/(1 + cT)$. With $c$ = 0.1 for standard care and 0.06 for the new treatment, 10-year restricted mean survival time is 10 / 2 = 5.00 years and 10 / 1.6 = 6.25 years, a gain of 1.25 years.

6. The common error. An analyst who applies the conditional hazard ratio of 0.6 as a constant ratio to the marginal standard care curve obtains (1 + 0.1t) to the power minus 1.2. Its 10-year area is 1 minus 2 to the power minus 0.2, divided by 0.02, which is (1 minus 0.8706) / 0.02 = 6.47 years, and 10-year survival is 0.25 raised to the power 0.6, which is 0.435, against the correct 0.391. The life-year gain becomes 1.47 years, overstating the correct 1.25 by 0.22 years, or about 18%, and so would any QALY gain built on it.

Conditional and marginal effects in economic evaluation

A coefficient from a shared frailty model compares two patients with the same frailty, that is, two patients in the same cluster. Balan and Putter show that for gamma frailty the marginal hazard ratio is pulled towards 1 as cumulative hazard grows, while the hazard ratios of an ordinary Cox model have a marginal interpretation, averaged over unmeasured covariates such as centre effects. In their EORTC table the coefficient for node-positive disease was 0.97 in the Cox model and 0.79 in the gamma frailty model, close to the 0.77 from fixed centre effects. This shift reflects adjustment for centre and runs opposite to the attenuation just described. The robust standard error of the Cox coefficient (0.26) was more than twice the model-based one (0.11).

A cost-effectiveness model for a national population needs survival averaged across the centres and settings where the technology will be used. In a simulation study for the NICE Decision Support Unit's Technical Support Document 21 (TSD 21), Rutherford and colleagues generated survival with an unobserved individual-level frailty term and stated that the marginal survival function was the quantity of interest, adding that the shape of the marginal hazard can be very different from the conditional hazards. A partitioned survival model or state transition model fed with frailty-model outputs should therefore integrate over the frailty distribution, as in step 3, instead of plugging a conditional hazard ratio into a population curve, as in step 6. The same distinction affects extrapolation, since a population hazard can fall even when every cluster's conditional hazard is rising.

Shared frailty beside neighbouring methods

The univariate frailty model gives each individual their own frailty. With no clustering, Balan and Putter judge it to involve a large degree of speculation, because frailty and a time-varying covariate effect cannot easily be told apart. The shared version is identifiable through within-cluster correlation, and the authors judge random effects particularly sensible when clusters can be seen as a random sample from a larger population. Fixed centre effects can be considered when clusters are large, but their estimates for small clusters are highly variable and would benefit from shrinkage. In meta-analysis, the random effects model applies the same hierarchical idea to studies, and a mixed effects model uses random effects for clustered outcomes that are not times to event. Cluster randomised designs raise a related design issue through the intraclass correlation.

Limitations and common misreadings

Most errors with shared frailty come from forgetting that the model describes patients within a cluster and rests on assumptions about the frailty that the data can only partly check. Several can change a cost-effectiveness result.

1. Reading the frailty coefficient as a population effect. It is a within-cluster effect, and a decision model needs the marginal survival curve.
2. Treating the variance as fixed. It can shrink as covariates are added and depends on the assumed distribution.
3. Ignoring the conditional proportional hazards check. Balan and Putter suggest Schoenfeld residuals with the log frailties as an offset, most sensible when clusters exceed about 10 subjects. With small clusters, frailty and non-proportional hazards are confounded.
4. Overlooking informative censoring. Valid inference needs censoring that does not depend on the frailty.
5. Over-reading small clusters. Empirical Bayes shrinkage tempers them but does not create information that is not in the data.

Sources

• Balan TA, Putter H. A tutorial on frailty models. Statistical Methods in Medical Research. 2020;29(11):3424-3454. https://doi.org/10.1177/0962280220921889
• Balan TA, Putter H. frailtyEM: an R package for estimating semiparametric shared frailty models. Journal of Statistical Software. 2019;90(7):1-29. https://doi.org/10.18637/jss.v090.i07
• Glidden DV, Vittinghoff E. Modelling clustered survival data from multicentre clinical trials. Statistics in Medicine. 2004;23(3):369-388. https://doi.org/10.1002/sim.1599
• Hougaard P. Frailty models for survival data. Lifetime Data Analysis. 1995;1(3):255-273. https://doi.org/10.1007/BF00985760
• Rutherford MJ, Lambert PC, Sweeting MJ, Pennington B, Crowther MJ, Abrams KR, Latimer NR. NICE DSU Technical Support Document 21: Flexible methods for survival analysis. Sheffield: Decision Support Unit, ScHARR, University of Sheffield; 2020 (updated March 2022). Section 5.4. https://doi.org/10.15131/shef.data.33849118
• Vaupel JW, Manton KG, Stallard E. The impact of heterogeneity in individual frailty on the dynamics of mortality. Demography. 1979;16(3):439-454. https://doi.org/10.2307/2061224