Integrating Gamma Frailty in Cluster Likelihoods

Why Integrate Frailty Out?

A shared frailty represents unobserved cluster-level variation in event risk. Conditional on a cluster’s frailty, the event processes may be modeled using familiar survival likelihoods. But when that frailty is unobserved, the likelihood for the observed data must average the conditional likelihood over the frailty distribution. For gamma frailty, this integration has a closed form when the frailty multiplies the event intensity.

The calculation is useful because it turns a conditional likelihood into an observed-data likelihood for the cluster. Its form also makes clear what information from the observed event history enters the frailty integration: the total number of events and the cluster’s cumulative exposure to its baseline risk.

Definition: Let $U_j$ be the unobserved frailty shared by members of cluster $j$. In this tutorial, $U_j$ has a gamma distribution with mean 1 and variance $\theta>0$, parameterized by shape $1/\theta$ and rate $1/\theta$. Thus its density is $$ f(u)=\frac{(1/\theta)^{1/\theta}}{\Gamma(1/\theta)} u^{1/\theta-1}\exp(-u/\theta), \qquad u>0. $$ The parameter $\theta$ measures between-cluster frailty variation; as $\theta$ tends to zero, $U_j$ tends to the constant 1.

Writing the Conditional Cluster Likelihood

Suppose the event intensity for member $i$ of cluster $j$, conditional on the shared frailty, is $$ \lambda_{ij}(t\mid U_j)=U_jY_{ij}(t)\lambda_0(t)\exp\{\beta^\mathsf{T}X_{ij}(t)\}. $$ Here $Y_{ij}(t)$ indicates whether member $i$ is at risk, $\lambda_0(t)$ is the baseline hazard, and $X_{ij}(t)$ is a covariate vector. This form also accommodates recurrent events when the intensity and risk indicator are defined for the event history being modeled. The essential condition for the integration below is that the same $U_j$ multiplies each member’s intensity and does not otherwise enter the likelihood.

Assume the observed event processes are conditionally independent within the cluster given $U_j$ and the modeled covariates and histories. Let $D_j$ be the total number of observed events in cluster $j$. Define $A_j$ as the summed cumulative baseline-risk exposure over its members:

A_j=\sum_i\int Y_{ij}(t)\lambda_0(t) \exp\{\beta^\mathsf{T}X_{ij}(t)\}\,dt.

For a fully specified conditional event-process likelihood, the event intensity factors contribute one factor of $U_j$ per event, while the integrated intensity contributes $\exp(-U_jA_j)$. Therefore the cluster likelihood conditional on $U_j=u$ can be written

L_j(u)=C_j u^{D_j}\exp(-uA_j),

where $C_j$ collects the event-time factors that do not involve $u$, such as baseline hazards and covariate terms evaluated at event times. It can also include other terms independent of $u$. This is a full conditional likelihood contribution, not a partial likelihood with the baseline hazard automatically removed.

The setup assumes that the event observation and censoring scheme is correctly represented by the likelihood and that censoring is ignorable for this likelihood, for example, conditionally independent of the event process given the modeled variables and frailty. The gamma distribution must describe the same shared frailty that multiplies the intensities. If members have distinct frailties, or the frailty enters the intensity in a different way, the expression $C_j u^{D_j}\exp(-uA_j)$ need not apply.

Integrating the Gamma Frailty

The observed-data contribution is the conditional likelihood averaged over the gamma density. For $D_j$ a nonnegative integer and $A_j\geq 0$, this gives

\begin{aligned} L_j &=\int_0^\infty L_j(u)f(u)\,du\\ &=C_j\frac{(1/\theta)^{1/\theta}}{\Gamma(1/\theta)} \int_0^\infty u^{D_j+1/\theta-1} \exp\{-(A_j+1/\theta)u\}\,du. \end{aligned}

For $a>0$ and $b>0$, the gamma integral is $\int_0^\infty u^{a-1}e^{-bu}\,du=\Gamma(a)/b^a$. Applying it with $a=D_j+1/\theta$ and $b=A_j+1/\theta$ yields the closed form:

Result: Under the shared gamma frailty model and likelihood assumptions above, the integrated cluster contribution is

L_j =C_j\frac{\Gamma(D_j+1/\theta)}{\Gamma(1/\theta)} \frac{\theta^{D_j}}{(1+\theta A_j)^{D_j+1/\theta}}.

This expression is the conditional likelihood integrated over the unobserved frailty. The data from the cluster enter the frailty-dependent part through $D_j$ and $A_j$.

An equivalent form follows from the gamma Laplace transform and its derivatives:

\mathbb{E}\!\left[U_j^{D_j}e^{-A_jU_j}\right] =\frac{\Gamma(D_j+1/\theta)}{\Gamma(1/\theta)} \frac{\theta^{D_j}}{(1+\theta A_j)^{D_j+1/\theta}}.

This expectation is the frailty-dependent factor; multiplying it by $C_j$ gives the likelihood contribution. For $D_j=0$, the gamma-function ratio and $\theta^{D_j}$ reduce to 1, so the factor is $(1+\theta A_j)^{-1/\theta}$. As $\theta$ tends to zero, the factor tends to $e^{-A_j}$, the value obtained when frailty is fixed at 1. These limiting cases are useful checks on algebra and software.

Worked Example: Two Events in One Cluster

Consider a hypothetical cluster with $D_j=2$ events and cumulative exposure $A_j=1.5$. Suppose the frailty variance is $\theta=0.5$, and the event-time factors collected in $C_j$ equal $0.12$. These values are illustrative, not from a real study. The gamma shape and rate are both $1/\theta=2$. The integrated frailty factor is

\begin{aligned} \mathbb{E}[U_j^2e^{-1.5U_j}] &=(0.5)^2\frac{\Gamma(2+2)}{\Gamma(2)} (1+0.5(1.5))^{-(2+2)}\\ &=0.25\frac{6}{1}(1.75)^{-4}\\ &=\frac{1.5}{9.3789}\\ &\approx 0.1599. \end{aligned}

Thus the integrated likelihood contribution is $L_j=0.12(0.1599)\approx 0.01919$. This is a likelihood contribution, not a probability that the cluster belongs to a particular risk group. Its numerical size depends on the event-time factors and the chosen data representation as well as on the frailty factor.

Worked Example: A Cluster With No Observed Events

Suppose, hypothetically, that a cluster has no observed events, cumulative exposure $A_j=2.4$, and frailty variance $\theta=0.25$. Let $C_j=1$, since there are no event-time factors in this simplified example. With $D_j=0$, the frailty-integrated contribution is

\begin{aligned} L_j &=(1+\theta A_j)^{-1/\theta}\\ &=(1+0.25(2.4))^{-4}\\ &=1.6^{-4}\\ &=\frac{1}{6.5536}\\ &\approx 0.1526. \end{aligned}

This likelihood contribution concerns observing zero events over the specified exposure under the model. It is not the probability that a cluster’s frailty is low. In particular, even a cluster with no events can contribute information about the frailty variance through its observed exposure and outcome.

What the Closed Form Reveals

The gamma integration replaces an unobserved cluster effect with a specific dependence of the likelihood on the cluster’s event count and cumulative exposure. The event count appears in both $\Gamma(D_j+1/\theta)$ and the power of $\theta$; exposure appears in $1+\theta A_j$. The constant $C_j$ still matters for the full likelihood, but it does not affect the integration over $U_j$ when it is independent of frailty.

A short justification for why the event count is the relevant exponent is to inspect the product of event intensities. Each event contributes an intensity proportional to $u$. Multiplying across $D_j$ events gives $u^{D_j}$. The integrated intensity over all members and follow-up is $uA_j$, giving the survival factor $\exp(-uA_j)$. These two pieces combine before integration, which is why it is important to form the entire conditional cluster likelihood before averaging over frailty.

Worked Example: Three Events and the Zero-Frailty-Variance Check

Take hypothetical values $D_j=3$, $A_j=2$, $\theta=0.2$, and $C_j=0.05$. The shape and rate are both $5$. The frailty factor is

\begin{aligned} \mathbb{E}[U_j^3e^{-2U_j}] &=(0.2)^3\frac{\Gamma(8)}{\Gamma(5)} (1+0.2(2))^{-8}\\ &=0.008(7\cdot6\cdot5)(1.4)^{-8}\\ &=\frac{1.68}{14.7579}\\ &\approx 0.1138. \end{aligned}

The integrated cluster contribution is $0.05(0.1138)\approx 0.005692$. As a check, if frailty variance were instead zero, the frailty-dependent term would be $1^3e^{-2}=e^{-2}\approx0.1353$, giving a contribution of about $0.006767$. The two values differ because the first model averages over heterogeneous frailties; the calculation does not imply that frailty necessarily lowers every cluster’s likelihood contribution.

Assumptions, Diagnostics, and Alternatives

The closed form is convenient, but it does not make the model assumptions automatic. Before using it, check the following features of the likelihood construction:

One shared multiplier: The same $U_j$ must multiply every modeled intensity in cluster $j$. If some members or event types have different frailties, the cluster likelihood generally has a different form.
Correct event count and exposure: $D_j$ must count the events represented by the likelihood, and $A_j$ must sum the corresponding at-risk intensity exposures. Omitting a member, event, or risk period changes the integral.
Conditional independence: The modeled processes must have the assumed conditional likelihood given frailty and covariates. Residual dependence beyond the shared frailty is not accounted for by this integration.
Appropriate censoring treatment: Informative censoring not represented in the likelihood can distort the contribution. Censoring is not made ignorable merely by adding a frailty term.
Frailty distribution and independence: The gamma density must be an appropriate model for $U_j$, typically assumed independent of the included covariates. The closed form depends on this distributional choice.

A practical diagnostic is to calculate $D_j$ and $A_j$ independently from the event and at-risk records for a few clusters, then verify that the integrated formula reduces to $C_je^{-A_j}$ as $\theta$ approaches zero. Check also that $A_j\geq0$, that the gamma-function arguments are positive, and that computed log-likelihood contributions are finite. For numerical work, it is often safer to evaluate the logarithm of the contribution, using log-gamma calculations, than to multiply very small factors directly.

The main alternatives answer different modeling or computational needs. Treating $U_j$ as fixed at 1 removes cluster heterogeneity; it is not the same as integrating over a gamma distribution. Keeping a conditional likelihood in which $U_j$ is known is appropriate only when the frailty is observed. For this gamma model and likelihood structure, analytic integration avoids numerical integration. With a different frailty distribution or a likelihood that does not reduce to $u^{D_j}e^{-uA_j}$, numerical integration may be needed. In every case, integrating a frailty changes the likelihood, not the underlying observed event records.

A common mistake is to integrate each member’s likelihood over a separate frailty distribution and then multiply the results, even though the model specifies one shared frailty. That calculation treats the members’ frailties as independent. The correct shared-frailty calculation first multiplies the members’ conditional likelihoods given the same $u$, combines their event counts and exposures, and then integrates once over $u$. Another error is to replace $A_j$ by the number of people in a cluster: exposure is an integrated intensity, not a cluster size.

Key takeaway: When a shared gamma frailty multiplies every event intensity, the conditional cluster likelihood has frailty-dependent part $u^{D_j}e^{-uA_j}$. Integrating this part against the mean-one gamma density gives a closed-form contribution determined by the total event count, cumulative exposure, and frailty variance.

Check Your Understanding

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

For $D_j=0$, $A_j=1$, and $\theta=1$, what is the frailty-integrated factor when $C_j=1$?
A cluster has two members, each with one event. Under a shared frailty, what is $D_j$, and how many times should the frailty density be integrated?
In one sentence, explain why the conditional likelihood includes $u^{D_j}$.
What does the integrated contribution reduce to as $\theta$ tends to zero?

Solutions. (1) The no-event formula gives $(1+\theta A_j)^{-1/\theta}=(1+1)^{-1}=0.5$. (2) $D_j=2$, and the combined conditional likelihood is integrated once over the one shared frailty. (3) Each of the $D_j$ event intensities contributes one multiplicative factor of $u$, so their product contains $u^{D_j}$. (4) The frailty-dependent factor tends to $e^{-A_j}$, as if $U_j=1$ for every cluster.