Marginal Survival From a Frailty Model

From Conditional Survival to Population Survival

A frailty model describes event risk conditional on a latent, positive risk multiplier. Its conditional survival function is not, by itself, the survival curve for the population: individuals differ in frailty, and the population curve averages over that distribution. The averaging step has a useful general form. It turns the cumulative hazard and covariate risk into the argument of the frailty distribution's Laplace transform.

This tutorial focuses on that transformation and its calculation for gamma and lognormal frailty. The result is a population-average survival probability at a specified time and covariate value. It does not require estimating a model from observed data, and it should not be confused with a mean cumulative function for recurrent-event burden.

Definition: Let \(Z\) be a positive frailty with \(E(Z)=1\), and let \(X=x\) be a fixed covariate value. Suppose the conditional hazard is \(\lambda(t\mid Z,x)=Z\exp(\beta^\mathsf{T}x)\lambda_0(t)\). Write \(H_0(t)=\int_0^t\lambda_0(s)\,ds\) for the cumulative baseline hazard and \(r(x)=\exp(\beta^\mathsf{T}x)\) for the covariate risk multiplier. The marginal survival is the survival probability averaged over the frailty distribution: \(S_m(t\mid x)=E_Z\{S(t\mid Z,x)\}\).

The model is conditional on the frailty: for a subject with frailty \(Z=z\), the hazard at time \(t\) is multiplied by both \(z\) and \(r(x)\). The baseline hazard describes the hazard when \(z=1\) and \(x=0\). The mean-one convention fixes the scale of frailty relative to the baseline hazard; the frailty variance describes heterogeneity around that mean.

For the formula below, assume that \(Z\) is positive with a specified distribution, that its distribution is the same at the covariate value \(x\), and that the conditional hazard has the stated multiplicative form. We also assume the cumulative hazard exists over the time interval considered. No particular frailty distribution is required for the general result.

The Laplace-Transform Formula

Given \(Z=z\), the cumulative hazard through time \(t\) is \(z r(x)H_0(t)\). The standard relationship between a hazard and survival therefore gives \(S(t\mid z,x)=\exp\{-z r(x)H_0(t)\}\). Averaging this conditional survival over the frailty distribution yields the marginal survival.

\begin{aligned} S_m(t\mid x) &=E_Z\left[\exp\{-Zr(x)H_0(t)\}\right]\\ &=\mathcal{L}_Z\{r(x)H_0(t)\}, \end{aligned}

Here \(\mathcal{L}_Z(u)=E\{\exp(-uZ)\}\), for \(u\geq 0\), is the Laplace transform of the frailty distribution. The argument \(u=r(x)H_0(t)\) is dimensionless: the cumulative hazard is dimensionless, as is the relative-risk multiplier. This formula is the central result. The conditional survival is averaged, not the conditional hazard, and the expectation is taken over the frailty distribution at the population level.

Key result: Under the stated conditional hazard model, \(S_m(t\mid x)=\mathcal{L}_Z\{r(x)H_0(t)\}\). The formula requires the full frailty distribution, not only its mean and variance. Different distributions with the same mean and variance can give different marginal survival values.

The formula has basic checks. At time zero, \(H_0(0)=0\), so \(S_m(0\mid x)=\mathcal{L}_Z(0)=1\). Since \(Z\geq 0\), increasing cumulative hazard cannot increase survival. If \(H_0(t)\) tends to infinity and \(P(Z=0)=0\), then \(\exp\{-Zr(x)H_0(t)\}\) tends to zero for every positive \(Z\), and marginal survival tends to zero as well.

Gamma Frailty

For mean-one gamma frailty with variance \(\theta>0\), use a gamma distribution with shape \(1/\theta\) and rate \(1/\theta\). Its Laplace transform has a closed form. For \(u\geq0\), integration of the gamma density gives \(\mathcal{L}_Z(u)=(1+\theta u)^{-1/\theta}\). Thus the marginal survival can be calculated directly from the cumulative baseline hazard, covariate risk, and frailty variance.

S_m(t\mid x) = \left\{1+\theta r(x)H_0(t)\right\}^{-1/\theta}.

When \(\theta\) approaches zero, the mean-one gamma distribution approaches a point mass at \(Z=1\), and the expression approaches the no-frailty survival \(\exp\{-r(x)H_0(t)\}\). This is a useful limiting check. The parameter \(\theta\) is the variance of \(Z\), not the standard deviation.

Worked Example: Gamma-Frailty Survival at a Fixed Time

Suppose a hypothetical subject has covariate risk multiplier \(r(x)=1.5\), cumulative baseline hazard \(H_0(t)=0.8\), and the model uses mean-one gamma frailty with variance \(\theta=0.5\). Calculate the marginal survival at time \(t\).

First form the Laplace-transform argument: \(u=r(x)H_0(t)=1.5(0.8)=1.2\). The gamma formula gives

\begin{aligned} S_m(t\mid x) &=(1+\theta u)^{-1/\theta}\\ &=(1+0.5(1.2))^{-1/0.5}\\ &=1.6^{-2}\\ &=\frac{1}{2.56}=0.390625. \end{aligned}

The marginal survival is approximately \(0.391\), or 39.1%, at the specified time and covariate value. This is an average over the gamma frailty distribution. It is not the survival probability for every individual, nor is it a statement that 39.1% of a particular observed sample must survive.

Lognormal Frailty

For mean-one lognormal frailty, let \(W\sim N(0,\sigma^2)\) and define \(Z=\exp(W-\sigma^2/2)\). Then \(E(Z)=1\), and the frailty variance is \(\exp(\sigma^2)-1\). The shift by \(-\sigma^2/2\) is necessary: without it, the mean of \(Z\) would be \(\exp(\sigma^2/2)\), not one.

The lognormal Laplace transform is an expectation over a normal random variable. In general it has no closed form in elementary functions, so numerical integration is used. With \(Y\sim N(0,1)\), write \(W=\sigma Y\). Then

\mathcal{L}_Z(u) = \int_{-\infty}^{\infty} \exp\left[-u\exp\left(\sigma y-\frac{\sigma^2}{2}\right)\right] \frac{\exp(-y^2/2)}{\sqrt{2\pi}}\,dy.

The integrand is bounded between zero and the standard normal density, making numerical evaluation straightforward with suitable quadrature. In practice, increase the number of quadrature points or otherwise check numerical convergence, especially when the frailty variance or \(u\) is large. The answer should lie in \([0,1]\); as \(u\) increases, the transform must not increase.

Worked Example: Numerical Lognormal-Frailty Survival

Consider a hypothetical model with \(r(x)=1\), \(H_0(t)=1\), and mean-one lognormal frailty with \(\sigma^2=0.5\). The variance of \(Z\) is \(\exp(0.5)-1\), approximately \(0.6487\). Calculate the marginal survival using numerical integration.

Here \(u=r(x)H_0(t)=1\), and \(Z=\exp(\sqrt{0.5}Y-0.25)\). The required expectation is

S_m(t\mid x) = \int_{-\infty}^{\infty} \exp\left[-\exp\left(\sqrt{0.5}\,y-0.25\right)\right] \frac{\exp(-y^2/2)}{\sqrt{2\pi}}\,dy.

For example, five-point Gauss–Hermite quadrature approximates this integral by summing the integrand's exponential term at the quadrature nodes with their normal-expectation weights. The resulting value is about \(0.449\) (a higher-order quadrature calculation gives a value close to \(0.448\)). Thus the estimated population-average survival at this time is about 44.8%, to three significant figures. The normal density in the integral is essential: simply evaluating the exponential at \(Z=E(Z)=1\) would give \(\exp(-1)\approx0.3679\), a conditional value rather than the required average.

What the Average Does—and Does Not—Mean

A frailty model averages conditional survival probabilities, not frailty values before calculating survival. In general, \(E\{\exp(-uZ)\}\) is not equal to \(\exp\{-uE(Z)\}\). The latter expression substitutes the mean frailty into the conditional survival model and discards heterogeneity.

Worked Example: Comparing the Average With Frailty Equal to One

Use the gamma-model values from the earlier calculation: \(\theta=0.5\), \(r(x)=1.5\), and \(H_0(t)=0.8\), so \(u=1.2\). The marginal survival is \(0.390625\). For a subject with \(Z=1\), the conditional survival is

S(t\mid Z=1,x)=\exp(-1.2)\approx0.3010.

The two values answer different questions. Approximately 30.1% is the model-based survival for frailty fixed at one; approximately 39.1% is the mean of the conditional survival probabilities across the gamma frailty distribution. These calculations do not imply that a high-frailty individual has better survival. Rather, the population average gives greater weight to the survival contributions of lower-frailty individuals, who retain a larger conditional survival probability at this time.

For the same \(u=1\), a gamma frailty model with variance \(0.5\) gives \((1+0.5)^{-2}=0.4444\), while the lognormal example with \(\sigma^2=0.5\) gives about \(0.448\). The models have the same mean frailty and approximately the same variance, but their distributions are not identical, so their Laplace transforms need not agree. A difference in marginal survival between these models is a distributional consequence, not a numerical error by itself.

Assumptions, Diagnostics, and Limitations

The calculation is appropriate only when the conditional model and frailty distribution correspond to the target population and outcome. Check these points before interpreting a result:

Multiplicative conditional hazard. The frailty and covariate risk must multiply the baseline hazard as specified. If the conditional effect changes form over time, the simple argument \(r(x)H_0(t)\) may not apply.
Frailty distribution at the covariate value. The formula uses the frailty distribution being averaged over. If frailty distributions differ across covariate groups, each group requires its own expectation.
Correct scale and mean. Verify that frailty is parameterized with mean one and that the baseline hazard uses the same scale. For a lognormal frailty, omitting the mean-correcting shift changes the model.
Numerical stability. For lognormal frailty, check that the computed probability is between zero and one, decreases as \(u\) increases, and is stable when quadrature accuracy is increased.
Target interpretation. Marginal survival averages over the specified frailty distribution. It is not the conditional survival for a named frailty value, and it does not by itself describe how event times in a particular cluster are related.

The marginal formula also depends on the assumed frailty distribution. The gamma case is convenient because its transform is explicit; the lognormal case generally needs numerical integration. If the distribution is misspecified, the resulting marginal survival can differ even if the conditional hazard form is otherwise plausible. The formula is a model-based population quantity, not an empirical survival estimate independent of assumptions.

Common mistake: Substituting \(E(Z)=1\) into the conditional survival and reporting \(\exp\{-r(x)H_0(t)\}\) as the marginal survival. This calculates survival at frailty \(Z=1\), not the average over frailty. The correction is to calculate \(E_Z[\exp\{-Zr(x)H_0(t)\}]\), or equivalently evaluate the frailty Laplace transform at \(r(x)H_0(t)\).

Check Your Understanding

Use the definitions and formulas in this tutorial. The solutions follow the questions.

For a conditional hazard multiplier \(Zr(x)\), what quantity is used as the argument of the frailty Laplace transform?
A mean-one gamma frailty has variance \(0.25\), and \(r(x)H_0(t)=2\). Calculate marginal survival.
For mean-one lognormal frailty with \(\sigma^2=0.2\), write the expression for \(Z\) in terms of a standard normal random variable \(Y\).
Why is \(\exp\{-r(x)H_0(t)\}\) generally not the marginal survival when frailty is variable?

Solutions. (1) The argument is \(u=r(x)H_0(t)\). (2) \(S_m=(1+0.25(2))^{-1/0.25}=1.5^{-4}=1/5.0625\approx0.1975\). (3) \(Z=\exp(\sqrt{0.2}\,Y-0.1)\), which has mean one. (4) That expression is the conditional survival at \(Z=1\). Marginal survival averages \(\exp\{-Zr(x)H_0(t)\}\) over the full frailty distribution.