Frailty Model (Survival Analysis)

Frailty model is a random-effects extension of the proportional hazards model in biostatistics, demography, and genetic epidemiology that incorporates an unobservable, non-negative random variable—termed frailty—acting multiplicatively on the baseline hazard rate. Introduced to biodemography by James W. Vaupel, Kenneth G. Manton, and Eric Stallard in 1979 to explain why cohort mortality curves decelerate at advanced ages relative to individual biological aging, and extended to multivariate clustered failure-time data by David G. Clayton (1978) and Philip Hougaard (1986), frailty modeling serves two distinct statistical roles. In univariate frailty models, the latent variable captures unobserved individual heterogeneity (omitted genetic, physiological, or environmental risk factors) that triggers selective survivorship ("weeding out" of high-risk individuals). In shared and correlated frailty models, a common random effect induces dependence among clustered or recurrent event times—such as monozygotic twins, patients treated within the same surgical hospital, or repeated malaria infections within a single child.

Individual Conditional Hazard vs. Marginal Population Survival

Let individual i (i = 1, …, n) possess an observed covariate vector x_i and an unobserved, time-invariant positive random frailty Z_i > 0 drawn from a probability density g(z) normalized to unit expectation, E[Z_i] = 1, and finite variance Var(Z_i) = θ ≥ 0. Conditional on Z_i = z_i and x_i, the instantaneous hazard rate h(t | x_i, z_i) follows a multiplicative proportional hazards specification:

where h₀(t) is the baseline hazard for a reference individual with x = 0 and unit frailty z = 1, and H_cond(t | x_i) = ∫₀^t h₀(u) exp(x_i^⊤ β) du = H₀(t) exp(x_i^⊤ β) is the conditional cumulative hazard. An individual with z_i = 2.5 experiences 2.5 times the instantaneous hazard of an average peer (z = 1) with identical measured covariates x_i.

Because Z_i is unobserved, conditional survival S(t | x_i, z_i) = exp[−z_i H_cond(t | x_i)] must be integrated over the density g(z) to obtain the observable marginal (population-averaged) survival function S_marg(t | x_i), which is evaluated directly by the Kaplan–Meier estimator:

where ℒ_Z(s) = E[exp(−sZ)] is the Laplace transform of the frailty distribution g(z) evaluated at argument s = H_cond(t | x_i). Differentiating −ln S_marg(t | x_i) with respect to t yields the marginal population hazard h_marg(t | x_i) = E[Z | T ≥ t, x_i] × h_cond(t | x_i).

Gamma Frailty and Artifactual Hazard Deceleration

The most widely applied frailty specification is the one-parameter Gamma distribution Z ∼ Gamma(k = 1/θ, λ = 1/θ) with mean E[Z] = 1, variance Var(Z) = θ, and Laplace transform ℒ_Z(s) = (1 + θ s)^−1/θ. Substituting s = H_cond(t | x) yields closed-form expressions for the marginal survival function, the mean frailty among survivors at time t, and the marginal cohort hazard:

This equation reveals the selective survivorship paradox identified by Vaupel, Manton, and Stallard (1979): as time t progresses, high-frailty individuals (z > 1) fail early and exit the risk set, causing the average frailty among surviving cohort members E[Z | T ≥ t] to decline monotonically. Even if every individual's biological hazard h₀(t) rises exponentially with age (a Gompertz law), the observed population hazard h_marg(t) plateaus toward an asymptote in old age. Furthermore, the marginal hazard ratio between an exposed group (x = 1) and an unexposed group (x = 0) shrinks toward 1.0 as t increases, explaining why omitting balanced prognostic covariates in a randomized clinical trial biases standard Cox hazard ratio estimates toward the null (in contrast to the accelerated failure time model).

Different parametric choices for g(z) produce distinct patterns of late-life hazard convergence and multivariate cluster dependence, summarized via the Hougaard power-variance family:

As Philip Hougaard (1986) demonstrated, the Positive Stable distribution with stability index α ∈ (0, 1) (where E[exp(−sZ)] = exp(−s^α)) is the only frailty distribution under which proportional hazards at the individual conditional level (h(t | x, z)) also preserves proportional hazards at the marginal population level, shrinking the marginal log-hazard ratio by the constant attenuation factor α: β_marg = α β_cond.

Shared Frailty in Clustered Survival and Recurrent Events

When observations are grouped into G natural clusters i = 1, …, G with n_i members j = 1, …, n_i (such as patients nested within G clinical centers, siblings within families, or sequential asthma exacerbations within patient i), the shared frailty model assigns a common cluster-level frailty Z_i to all members of cluster i:

Integrating out Z_i ∼ Gamma(1/θ, 1/θ) across a bivariate cluster (j = 1, 2) yields the Clayton–Oakes survival copula (Clayton 1978; Oakes 1989):

In this shared gamma frailty copula, the frailty variance θ relates directly to Kendall's rank correlation coefficient τ between clustered event times via τ = θ / (θ + 2). A likelihood-ratio test of H₀ : θ = 0 (compared against a 50:50 mixture of χ²₀ and χ²₁ distributions because θ = 0 lies on the boundary of the parameter space θ ≥ 0) tests for statistically significant intra-cluster correlation or recurrent-event overdispersion relative to an Andersen–Gill Poisson process.

Estimation Algorithms: EM, Penalized Partial Likelihood, and Laplace Approximation

Fitting shared frailty models with semi-parametric baseline hazards (estimated via the Nelson–Aalen estimator step function) proceeds via three classical computational strategies:

1. Expectation–Maximization (EM) Algorithm (Gill 1985; Nielsen et al. 1992): In the E-step, the posterior expectation of each cluster's frailty given its D_i observed events, E[Z_i | data] = (1/θ + D_i) / [1/θ + ∑_j H₀(y_ij) exp(x_ij^⊤ β)], is computed analytically for gamma frailty; in the M-step, β is updated via a standard Cox partial likelihood with offset ln *Ê*[Z_i].
2. Penalized Partial Likelihood (PPL) (Therneau and Grambsch 2000): Writing random effects on the log scale u_i = ln Z_i transforms both gamma and Gaussian frailty estimation into maximizing a Cox partial log-likelihood penalized by the log-density of u, solved rapidly via Newton–Raphson iteration in clinical survival software.

1. Vaupel, J. W., Manton, K. G., & Stallard, E. (1979). "The impact of heterogeneity in individual frailty on the dynamics of mortality." Demography, 16(3), 439–454. https://doi.org/10.2307/2061224
2. Clayton, D. G. (1978). "A model for association in bivariate life tables and its application in epidemiological studies of familial tendency in chronic disease incidence." Biometrika, 65(1), 141–151. https://doi.org/10.1093/biomet/65.1.141
3. Hougaard, P. (1986). "A class of multivariate failure time distributions." Biometrika, 73(3), 671–678. https://doi.org/10.1093/biomet/73.3.671
4. Oakes, D. (1989). "Bivariate survival models induced by frailties." Journal of the American Statistical Association, 84(406), 487–493. https://doi.org/10.1080/01621459.1989.10478795
5. Wienke, A. (2010). Frailty Models in Survival Analysis. Chapman and Hall/CRC. https://doi.org/10.1201/9781420073911