Gamma Frailty and Its Laplace Transform

Why the Gamma Frailty Transform Is Useful

A frailty variable represents latent multiplicative risk: conditional on its value, an individual or cluster has a proportional hazards model, but the frailty itself is unobserved. To move from that conditional model to quantities averaged over the frailty distribution, we need to integrate over the latent variable. The gamma distribution is especially useful for this calculation because its Laplace transform has a simple closed form.

The goal here is to make that calculation transparent. We will parameterize a positive frailty to have mean one, derive its density and variance, obtain its Laplace transform, and use the transform to integrate a conditional proportional hazards model. This integration depends on assumptions about the frailty and the conditional model; the algebra alone does not make those assumptions true.

Mean-One Gamma Frailty

Definition: For $\theta>0$, a mean-one gamma frailty $U$ with variance parameter $\theta$ has a gamma distribution with shape $1/\theta$ and rate $1/\theta$. Its density is $$ f_U(u)=\frac{u^{1/\theta-1}\exp(-u/\theta)}{\Gamma(1/\theta)\theta^{1/\theta}}, \qquad u>0, $$ where $\Gamma(a)=\int_0^\infty v^{a-1}\exp(-v)\,dv$ is the gamma function. Thus $E(U)=1$ and $\operatorname{Var}(U)=\theta$. When $\theta=0$, the mean-one frailty is defined by the limiting point mass $U=1$.

The rate convention matters. A gamma distribution with shape $k$ and rate $r$ has density $r^k u^{k-1}\exp(-ru)/\Gamma(k)$, mean $k/r$, and variance $k/r^2$. Setting $k=r=1/\theta$ gives mean one and variance $\theta$. Equivalently, in the shape-scale convention, the shape is $1/\theta$ and the scale is $\theta$. The rate and scale are reciprocals, not interchangeable labels.

The parameter $\theta$ controls dispersion, not the mean. Smaller positive $\theta$ gives a distribution more concentrated around one; larger $\theta$ gives greater variation in latent risk. At $\theta=0$, there is no frailty variation: every unit has multiplier one.

Worked Example: Writing the Density and Checking Its Moments

Suppose a hypothetical model uses $\theta=0.25$. Find the gamma shape and rate, write the density, and verify its mean and variance.

The shape and rate are both $1/\theta=1/0.25=4$. Therefore $$ f_U(u)=\frac{4^4}{\Gamma(4)}u^{4-1}\exp(-4u) =\frac{256}{6}u^3\exp(-4u),\qquad u>0. $$ Here $\Gamma(4)=3!=6$, so the density coefficient is approximately $42.67$.

Using the shape-rate moment formulas, $E(U)=4/4=1$, and $\operatorname{Var}(U)=4/4^2=0.25$. These agree with the intended mean-one, variance-$\theta$ parameterization. The density is defined only for positive $u$; its value at a point is a density, not a probability that $U$ equals that point.

Deriving the Laplace Transform

For $s\geq 0$, the Laplace transform of a nonnegative random variable $U$ is $L_U(s)=E\{\exp(-sU)\}$. For gamma frailty, substitute the density into this expectation. The integral is finite because the exponent has a negative rate $s+1/\theta>0$.

\begin{aligned} L_U(s) &=\int_0^\infty \exp(-su) \frac{u^{1/\theta-1}\exp(-u/\theta)} {\Gamma(1/\theta)\theta^{1/\theta}}\,du\\ &=\frac{1}{\Gamma(1/\theta)\theta^{1/\theta}} \int_0^\infty u^{1/\theta-1} \exp\{-(s+1/\theta)u\}\,du\\ &=\frac{1}{\Gamma(1/\theta)\theta^{1/\theta}} \frac{\Gamma(1/\theta)}{(s+1/\theta)^{1/\theta}}\\ &=(1+\theta s)^{-1/\theta}. \end{aligned}

The third line uses the gamma integral identity $\int_0^\infty u^{a-1}\exp(-bu)\,du=\Gamma(a)/b^a$ for $a,b>0$. Substituting $a=1/\theta$ and $b=s+1/\theta$ yields the stated result. At $\theta=0$, the point-mass definition instead gives $L_U(s)=\exp(-s)$, which is also the limit of $(1+\theta s)^{-1/\theta}$ as $\theta$ tends to zero.

Worked Example: Evaluating a Gamma Frailty Transform

Let $U$ have mean one and variance $\theta=0.5$. Evaluate $E\{\exp(-2U)\}$ using the Laplace transform.

Here $s=2$. Substitution gives $$ L_U(2)=(1+0.5\times 2)^{-1/0.5} =2^{-2}=\frac{1}{4}=0.2500. $$ Thus the expectation of $\exp(-2U)$ over the hypothetical frailty distribution is $0.25$. The calculation integrates over all possible positive frailty values; it is not the value of $\exp(-2u)$ for one particular unit.

Using the Transform to Marginalize a Proportional Hazards Model

Suppose that, conditional on a frailty value $U=u$, a unit with fixed covariates $x$ has hazard $\lambda(t\mid u,x)=u\lambda_0(t)\exp(x^\mathsf{T}\beta)$. Here $\lambda_0(t)$ is a nonnegative baseline hazard, $\beta$ is a coefficient vector, and the conditional cumulative hazard is $u\exp(x^\mathsf{T}\beta)\Lambda_0(t)$, where $\Lambda_0(t)=\int_0^t\lambda_0(v)\,dv$. The corresponding conditional survival is the exponential of minus that cumulative hazard.

To integrate this conditional model over $U$, assume the frailty follows the mean-one gamma distribution above and is independent of the covariates $x$. The frailty is drawn for the unit (or cluster, under a shared-frailty structure) and remains fixed over the time interval being modeled. Conditional on $U$ and $x$, the hazard specification must be appropriate for the event process. With $a(t,x)=\exp(x^\mathsf{T}\beta)\Lambda_0(t)$, taking the expectation of conditional survival gives

\begin{aligned} S_{\mathrm{marg}}(t\mid x) &=E_U\{\exp[-U a(t,x)]\}\\ &=L_U\{a(t,x)\}\\ &=\{1+\theta\exp(x^\mathsf{T}\beta)\Lambda_0(t)\}^{-1/\theta}. \end{aligned}

This is the direct use of the transform: the conditional survival has the form $\exp(-Us)$, so averaging it over $U$ is exactly evaluating $L_U(s)$. Differentiating the resulting marginal cumulative hazard, $-\log S_{\mathrm{marg}}(t\mid x)$, gives the marginal hazard

\lambda_{\mathrm{marg}}(t\mid x) = \frac{\lambda_0(t)\exp(x^\mathsf{T}\beta)} {1+\theta\exp(x^\mathsf{T}\beta)\Lambda_0(t)}.

To see the differentiation, write $-\log S_{\mathrm{marg}}=(1/\theta)\log(1+\theta a)$. Its time derivative is $a'(t,x)/(1+\theta a)$, and $a'(t,x)=\exp(x^\mathsf{T}\beta)\lambda_0(t)$. This marginal hazard is obtained from the marginal survival; it is not obtained by simply averaging the conditional hazard at baseline. In particular, replacing $U$ by $E(U)=1$ before integrating would discard the effect of frailty variation.

Worked Example: Marginalizing a Conditional Proportional Hazards Model

Consider a hypothetical model with constant baseline hazard $\lambda_0(t)=0.08$ per year, a covariate multiplier $\exp(x^\mathsf{T}\beta)=1.5$, frailty variance $\theta=0.5$, and time $t=3$ years. Calculate the marginal survival and marginal hazard at that time.

The baseline cumulative hazard is $\Lambda_0(3)=0.08\times3=0.24$ (unitless). Thus $a=1.5\times0.24=0.36$. The marginal survival is $$ S_{\mathrm{marg}}(3\mid x) =(1+0.5\times0.36)^{-2} =1.18^{-2} =\frac{1}{1.3924} \approx 0.7182. $$

The marginal hazard at year 3 is $$ \lambda_{\mathrm{marg}}(3\mid x) =\frac{0.08\times1.5}{1+0.5\times0.36} =\frac{0.12}{1.18} \approx 0.1017\ \text{per year}. $$ The calculation uses the same frailty variance and covariate multiplier in both quantities. The survival is a probability and has no time unit; the hazard has units of inverse years.

Assumptions, Checks, and Alternatives

The transform calculation is exact under its distributional and model assumptions, but it does not verify them. In an application, the modeler must justify the conditional hazard form and the role assigned to the frailty. In particular, the independence of $U$ and $x$ is part of the stated model; if latent risk and covariates are associated, this integration does not describe the same model. A shared frailty also represents a common latent multiplier for the members of a cluster, whereas an individual frailty represents a multiplier specific to a unit.

Useful algebraic checks include the support $u>0$, the moments $E(U)=1$ and $\operatorname{Var}(U)=\theta$, and the transform at $s=0$: $L_U(0)=1$. For $s\geq0$, the transform must be positive and no greater than one. As $\theta$ approaches zero, the frailty becomes constant at one and the transform approaches $\exp(-s)$. A formula that fails these checks often signals a rate-scale interchange or an incorrect exponent.

Gamma frailty is convenient because its transform yields a closed-form integral. Other frailty distributions, including lognormal and inverse Gaussian frailties, can also be used, but their transforms and resulting calculations differ. Numerical integration may be needed when a closed form is unavailable. A closed form is a computational advantage, not evidence that the gamma distribution is a good description of latent heterogeneity; conclusions can depend on the distributional choice.

Common pitfall: Treating $\theta$ as the gamma scale while also using $1/\theta$ as the rate produces an inconsistent density. Under the convention here, shape $=1/\theta$, rate $=1/\theta$, and scale $=\theta$. Also, $E(U)=1$ does not mean $U=1$ for every unit when $\theta>0$; setting $U=1$ before integrating removes the variation the frailty model was introduced to represent.

Check Your Understanding

Answer each question, then compare with the solutions below.

For $\theta=0.2$, what are the gamma shape and rate, and what are the mean and variance of $U$?
For mean-one gamma frailty with $\theta=0.25$, evaluate $L_U(3)$.
In a conditional proportional hazards model, set $a(t,x)=\exp(x^\mathsf{T}\beta)\Lambda_0(t)$. What expression results from averaging the conditional survival over gamma frailty?
Why is replacing $U$ with its mean before evaluating conditional survival generally not the same as integrating over $U$?

Solutions. (1) Shape and rate are both $1/0.2=5$; the mean is $5/5=1$, and the variance is $5/5^2=0.2$. (2) $L_U(3)=(1+0.25\times3)^{-4}=1.75^{-4}\approx0.1066$. (3) The average conditional survival is $L_U\{a(t,x)\}=[1+\theta a(t,x)]^{-1/\theta}$. (4) Conditional survival is a nonlinear function of $U$, so in general $E\{\exp(-aU)\}\ne\exp\{-aE(U)\}$. The transform performs the required expectation rather than substituting the mean.