Fitting Shared Frailty Models for Clustered Event Histories

Why Fit a Shared Frailty Model?

Individuals in the same clinic, family, center, or other cluster may have more similar event risks than individuals in different clusters. A Cox model that treats their risks as unrelated can miss this dependence. A shared frailty model represents cluster-level heterogeneity with a latent random multiplier: individuals in the same cluster share the same frailty, while different clusters have different frailties.

The model is useful when the scientific aim includes covariate effects conditional on this cluster-level risk and when a random-effects representation is appropriate. It does not replace careful definition of the event, time origin, censoring, or risk sets. As in ordinary Cox regression, those choices must match the study design.

Definition: For individual $i$ in cluster $j$, a shared gamma frailty Cox model specifies the event hazard conditional on the cluster frailty $Z_j$ as $$ \lambda_{ij}(t\mid Z_j,\mathbf{x}_{ij}(t)) = Z_j\lambda_0(t)\exp\{\boldsymbol{\beta}^{\mathsf T}\mathbf{x}_{ij}(t)\}. $$ Here $\lambda_0(t)$ is an unspecified baseline hazard, $\mathbf{x}_{ij}(t)$ is the covariate vector, and $Z_j$ is shared by every individual in cluster $j$. The frailties are independent across clusters and have a gamma distribution with mean 1 and variance $\theta$. Equivalently, when $\theta>0$, $Z_j$ has shape $1/\theta$ and scale $\theta$.

Conditional on the frailty, the model is a proportional hazards model: a one-unit increase in a covariate with coefficient $\beta_k$, holding other covariates and $Z_j$ fixed, multiplies the hazard by $\exp(\beta_k)$. This is a conditional hazard ratio. It is not generally the same as a population-average hazard ratio after averaging over the frailty distribution.

Data Structure and Assumptions

For a single event per person, the usual input has one row per person, a follow-up time, an event indicator coded 1 for the event and 0 for right censoring, measured covariates, and a cluster identifier. The cluster identifier must identify the group whose members share frailty; it is not automatically the individual identifier. If people can experience recurrent events, the data may instead have multiple start-stop rows per person. The row structure and event indicator must then represent the intended risk process.

Before fitting, check that follow-up times and event indicators are valid, cluster labels are consistent, and covariates are coded as intended. For start-stop data, verify that intervals have positive length, do not overlap improperly for an individual, and include only times when that individual is at risk. The fitted model also relies on substantive assumptions:

The hazard is proportional in the covariates conditional on frailty, with the chosen functional forms adequate for continuous and categorical covariates.
A single multiplicative frailty is shared within each specified cluster, and frailties are independent between clusters with the stated gamma distribution.
Clusters provide independent units of information; dependence within a cluster is represented by the shared frailty.
Censoring is conditionally non-informative for the event process, given the modeled covariates and frailty. The event and censoring times must be represented correctly.
Covariates are measured and included appropriately. A frailty term does not correct confounding or repair a misspecified time scale.

Worked Example: Fit a Model for One Event per Person

Suppose a hypothetical follow-up study records one event time or censoring time for each participant, with participants recruited from 18 clinics. The data frame `dat` contains `time` in months, `event` coded 1 for the event and 0 for censoring, `treatment` coded 1 for the active treatment and 0 for control, `age` in years, and `clinic` as the clinic identifier.

First check that each participant contributes one record, that `time` is positive, that `event` has only the values 0 and 1, and that `clinic` correctly identifies the 18 clusters. These checks support the proposed data layout; they do not establish proportional hazards or non-informative censoring.

A direct fit using the `survival` package is:

library(survival)
fit <- coxph(Surv(time, event) ~ treatment + age +
frailty(clinic, distribution = "gamma", method = "em"),
data = dat, ties = "efron")
summary(fit)

`Surv(time, event)` tells `coxph` that each row has an observed event or right-censoring time. The `frailty()` term specifies a gamma frailty shared by participants with the same clinic value. The Efron option specifies how tied event times are handled. Review the fitted output for the covariate coefficients, their standard errors and tests, and the estimated frailty variance. This command fits the specified model; it does not itself verify the model assumptions.

What the Frailty Term Changes

The frailty variance $\theta$ quantifies the spread of the cluster multipliers around their mean of 1. When $\theta=0$, the frailty distribution collapses to $Z_j=1$, yielding the ordinary Cox model without cluster-level frailty. A positive estimate indicates estimated heterogeneity between clusters under the gamma model; it is not a measure of the fraction of outcome variation “explained” by the cluster.

The distinction between conditional and marginal effects follows directly from integrating over the frailty. Let $H_0(t)=\int_0^t\lambda_0(u)\,du$, and consider a fixed covariate vector with $\eta=\boldsymbol{\beta}^{\mathsf T}\mathbf{x}$. Conditional survival is $\exp\{-Z_jH_0(t)e^\eta\}$. For a gamma variable with mean 1 and variance $\theta$, its Laplace transform is $E\{\exp(-sZ_j)\}=(1+\theta s)^{-1/\theta}$. Therefore:

S_{\mathrm{marg}}(t\mid\mathbf{x}) = E_{Z}\left[\exp\{-ZH_0(t)e^\eta\}\right] = \left\{1+\theta H_0(t)e^\eta\right\}^{-1/\theta}.

This marginal survival expression is not generally the survival function of a proportional hazards model with the same coefficient. Consequently, the coefficients reported from the shared frailty Cox fit describe covariate associations conditional on the latent cluster frailty. Do not label them population-average effects.

Worked Example: Calculate and Interpret a Conditional Hazard Ratio

Suppose a hypothetical fitted model reports a treatment coefficient of $\hat{\beta}_{\text{treatment}}=-0.357$, with estimated standard error 0.160, and an estimated gamma frailty variance of 0.420. Treatment is coded 1 for active treatment and 0 for control.

The conditional hazard ratio is obtained by exponentiating the coefficient:

\widehat{HR}_{\text{conditional}} = \exp(\hat{\beta}_{\text{treatment}}) = \exp(-0.357) = 0.700.

Under the fitted model, and at the same age and frailty value, the active-treatment group has an estimated event hazard 0.700 times that of the control group. Equivalently, the estimated conditional hazard is 30.0% lower. This is an association conditional on the model and its covariates, not by itself evidence that treatment caused a reduction.

For a Wald interval, the coefficient-scale endpoints are $-0.357\pm1.96(0.160)$, or $-0.6706$ and $-0.0434$. Exponentiating gives an approximate 95% confidence interval of $(0.511,\ 0.958)$ for the conditional hazard ratio. In repeated sampling under the model, this interval procedure has 95% coverage; it is not a 95% probability statement about this particular interval. The frailty variance estimate is 0.420 in the model’s variance parameterization; its implications depend on the assumed gamma frailty model.

Fitting Recurrent Event Histories

A person with recurrent events may contribute multiple at-risk intervals. If the intended event process is an Andersen-Gill counting-process formulation, one row per interval can be supplied using start-stop syntax. The frailty term still identifies the level at which the random multiplier is shared. For example, if individuals are nested in hospitals and each individual can contribute several event intervals, use the hospital identifier in `frailty()` to model a hospital-level shared effect, not the person identifier.

This choice encodes a substantive dependence structure. If person-level dependence is also important, a single hospital frailty does not automatically represent it. Likewise, changing the time scale or risk-set construction changes the model being fitted. Use the event-history representation justified by the study question and design.

Worked Example: Fit Cluster Frailty with Start-Stop Data

Consider a hypothetical recurrent-event data frame `longdat` with `start` and `stop` times in days, `event` coded 1 if an event occurs at the interval end and 0 otherwise, a time-varying `exposure`, baseline `age`, `person`, and `hospital`. Each person belongs to one hospital and contributes non-overlapping intervals while at risk.

After checking that each interval satisfies `start < stop`, that records are ordered correctly within person, and that event indicators occur only as defined by the chosen counting-process setup, fit:

fit_recurrent <- coxph(Surv(start, stop, event) ~ exposure + age +
frailty(hospital, distribution = "gamma", method = "em"),
data = longdat, ties = "efron")
summary(fit_recurrent)

The `Surv(start, stop, event)` response defines the start-stop risk intervals. The coefficient for `exposure` is conditional on the hospital frailty and the other modeled covariates. Because there are repeated records per person, the analyst must ensure the risk intervals correctly encode the event process; merely adding a frailty term does not fix duplicated or incorrectly assigned risk time. The frailty is hospital-level, so this specification does not separately introduce a person-level random effect.

Estimation, Diagnostics, and Alternatives

In `coxph`, the frailty term is estimated along with the covariate effects using the software’s frailty estimation procedure; it is not an extra fixed-effect indicator for each cluster. The exact printed presentation can depend on the `survival` package version and fitting options, so inspect the model output and record the variance estimate and the estimation method used. A small estimated variance does not prove that clusters are irrelevant, particularly when the number of clusters or events is limited.

Check the fit at several levels. Confirm that the fit converged and that the variance estimate is finite and interpretable in the specified parameterization. Look for unusually large standard errors, unstable coefficients, or warnings, and check whether a small number of clusters dominate the event information. Examine the proportional-hazards assumption and continuous-covariate functional forms using appropriate residual-based or graphical diagnostics. Compare plausible time-scale and covariate specifications when scientifically justified. These checks diagnose features of the model; they do not test the gamma distribution directly or establish that the frailty is independent of the covariates.

The main alternatives answer different questions:

Approach	Cluster dependence representation	Effect target or emphasis
Ordinary Cox model	No explicit cluster dependence term	Conditional PH coefficients under a model without latent cluster effects
Shared gamma frailty Cox model	One latent gamma multiplier shared within each cluster	Covariate effects conditional on frailty, plus an estimated frailty variance
Marginal Cox model with robust variance	Uses a working Cox model and adjusts uncertainty for within-cluster dependence	Population-average or marginal estimating target, depending on the formulation and conditions

A robust variance adjustment and a frailty model are not interchangeable. Robust variance primarily changes uncertainty estimates for a marginal estimating approach; it does not estimate a gamma frailty distribution. Conversely, a shared frailty model makes distributional assumptions and yields conditional effects. The number of independent clusters matters for inference, and few clusters can make variance estimation and standard errors unreliable. Do not add a robust cluster variance mechanically to a frailty fit without checking the software’s supported model and the intended inferential target.

Common Mistakes and Limitations

Using the wrong grouping variable. A patient ID does not represent a clinic-level shared frailty. Use the identifier for the scientific cluster whose members are assumed to share the latent multiplier.
Calling the coefficient marginal. The fitted covariate hazard ratio is conditional on frailty. Marginal hazards can change over time after averaging over frailty.
Assuming frailty repairs every dependence problem. A single gamma effect does not automatically account for multiple nested dependence levels, time-varying cluster effects, or misspecified event histories.
Treating the estimated variance as a causal cluster effect. It describes latent heterogeneity under the model; it does not identify a causal effect of belonging to a particular cluster.
Ignoring sparse cluster information. A large number of individuals concentrated in only a few clusters may still provide limited information about the between-cluster frailty distribution.

The model is most useful when its cluster structure and conditional interpretation match the scientific question. Its conclusions remain dependent on the gamma frailty assumption, the proportional hazards specification, the censoring assumptions, and correct event-history construction.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

In a model with people nested in clinics, which identifier belongs in the shared frailty term if the assumed latent effect is clinic-level?
A treatment coefficient is $-0.223$. Calculate its conditional hazard ratio and give a contextual interpretation.
What does a frailty variance estimate of zero correspond to in this model?
Why is a robust-variance Cox analysis not equivalent to fitting a shared gamma frailty model?

Solutions. (1) Use the clinic identifier, because the model assumes people in the same clinic share the same latent multiplier. (2) $\exp(-0.223)=0.800$ to three significant figures. Conditional on the frailty and other covariates, the treatment group’s estimated hazard is 0.800 times the comparison group’s hazard, or 20.0% lower. (3) It is the limiting no-frailty case, $Z_j=1$ for all clusters, reducing to an ordinary Cox model without cluster-level latent heterogeneity. (4) The robust approach adjusts uncertainty for clustering in a marginal estimating framework; it does not estimate the gamma frailty distribution or the same conditional effect.