Shared Frailty for Recurrent Event Clusters

Why Put a Shared Frailty in a Recurrent-Event Model?

People or other units in the same clinic, family, or community may experience recurrent events at persistently different rates. A shared frailty model represents one source of that clustering by assigning a common, unobserved multiplicative factor to the event intensities of members of each cluster. It therefore models heterogeneity as well as within-cluster dependence.

The interpretation depends on choices made before fitting the model. In particular, the analyst must specify the time scale, which events count as recurrences, who is at risk at each time, and which event-history features are included in the intensity. Holding those choices fixed, the frailty describes residual cluster-level variation in the conditional event rate; it does not change the definition of the recurrent-event estimand.

Definition: Let $N_{ij}(t)$ count recurrent events by time $t$ for member $j$ of cluster $i$, and let $Y_{ij}(t)$ indicate whether that member is under observation and at risk immediately before $t$. A shared frailty recurrent-event intensity model is $$ \lambda_{ij}(t\mid Z_i,\mathcal{H}_{ij}(t),X_{ij}(t)) = Y_{ij}(t)Z_i\lambda_0(t)\exp\{\beta^\mathsf{T}X_{ij}(t)\}. $$ Here $\lambda_{ij}$ is the conditional event intensity, $\mathcal{H}_{ij}(t)$ is the observed event history up to $t$, $X_{ij}(t)$ is a covariate vector, $\lambda_0(t)$ is a baseline intensity, and $Z_i$ is the cluster’s shared frailty. The intensity has units of events per unit time while at risk.

Assumptions and Interpretation

A common specification takes the frailties to be independent across clusters and independent of the modeled covariates, with $E(Z_i)=1$ and $\operatorname{Var}(Z_i)=\theta$. A gamma choice is $Z_i\sim\operatorname{Gamma}(1/\theta,1/\theta)$, using shape and rate parameters, when $\theta>0$; the limiting case $\theta=0$ means $Z_i=1$ with no frailty variation. The mean-one constraint separates the overall intensity scale, represented by $\lambda_0(t)$, from the relative cluster multiplier.

Conditional on the frailties and modeled covariates, clusters are independent. For the usual shared-frailty interpretation, members’ event processes are conditionally independent given their shared frailty and covariates. Their marginal processes are dependent because they share $Z_i$. Censoring and loss to follow-up must also be non-informative conditional on the information represented in the model; otherwise, observed risk time can be selectively truncated.

The model is a conditional intensity specification. It does not automatically assert that recurrent events are independent or that each person’s events form a Poisson process. Such a Poisson interpretation requires additional conditions, such as an intensity without event-history dependence and suitable fixed covariates. If event order or time since the previous event affects the rate, that dependence must be represented in the intensity while preserving the intended estimand.

For two covariate values $x_1$ and $x_0$, at the same time on the chosen time scale, and for the same frailty value, the conditional intensity ratio is $$ \frac{\lambda_{ij}(t\mid Z_i,x_1)}{\lambda_{ij}(t\mid Z_i,x_0)} = \exp\{\beta^\mathsf{T}(x_1-x_0)\}. $$ The shared frailty cancels. Thus, for a binary covariate with coefficient $\beta$, $\exp(\beta)$ is a conditional recurrent-event intensity ratio. It compares rates among units at risk at the same modeled time, conditional on their shared frailty and other covariates. It is not automatically a marginal rate ratio, a ratio of cumulative event counts, or a causal effect.

Worked Example: Interpreting a Covariate Effect at a Fixed Time Scale

Suppose a hypothetical study follows patients in clinics, using years since enrollment as the time scale. The model has baseline intensity $0.08$ events per patient-year, and treatment has estimated coefficient $\hat{\beta}=\log(1.5)$. Consider two otherwise comparable patients in the same clinic, both at risk at year 2, with shared frailty $Z_i=1.6$. Assume the fitted intensity includes no event-history term that differs between them.

For the untreated patient, $X=0$, so the conditional intensity is $$ 1.6(0.08)\exp\{0\}=0.128 $$ events per patient-year.

For the treated patient, $X=1$, so the conditional intensity is $$ 1.6(0.08)\exp\{\log(1.5)\} =1.6(0.08)(1.5)=0.192 $$ events per patient-year. The ratio is $0.192/0.128=1.5$. At this time and for this shared frailty, the fitted treatment-associated intensity is 50% higher. This comparison does not claim that treatment causes the difference, and it does not compare patients in different risk sets or on a different time scale.

What the Frailty Variance Adds

The frailty variance $\theta$ quantifies heterogeneity in the multiplicative cluster factor on the model’s conditional-intensity scale. Since $E(Z_i)=1$, a cluster with $Z_i=1$ has the baseline cluster multiplier; values above or below 1 scale its members’ intensities up or down. For any mean-one frailty distribution, the coefficient of variation of $Z_i$ is $\sqrt{\theta}$. In a gamma model, larger $\theta$ also means a more dispersed distribution of cluster multipliers.

This is a model-based description of latent heterogeneity, not a direct measurement of all differences between clusters. Observed cluster characteristics should be included as covariates when appropriate. If frailty is associated with omitted covariates, interpreting it as independent residual heterogeneity can be misleading.

Worked Example: Translating Frailty Variance to Relative Spread

Suppose a hypothetical recurrent-event analysis uses a mean-one gamma frailty with estimated variance $\theta=0.25$. The gamma shape and rate are both $1/\theta=4$, and the standard deviation of the frailty multiplier is $$ \sqrt{\operatorname{Var}(Z_i)}=\sqrt{0.25}=0.5. $$

Its coefficient of variation is $0.5/1=0.5$, or 50%. This describes the standard deviation of the latent cluster multipliers relative to their mean. For example, multipliers of $0.7$ and $1.4$, if they occurred, would scale all modeled member intensities in their respective clusters by those factors, holding time, covariates, and risk status fixed. The estimate does not mean that every cluster’s rate is exactly 50% above or below the average, nor does it give an individual patient’s event probability.

Why Shared Frailty Creates Dependence

A short derivation makes the dependence mechanism explicit. Consider a special case in which two members of a cluster have conditionally independent Poisson event counts over fixed observation windows, with conditional means $Z_i a$ and $Z_i b$, where $a$ and $b$ are their respective expected counts when $Z_i=1$. Given $Z_i$, their count covariance is zero. Marginally, however, their conditional means vary together with the same frailty.

\begin{aligned} \operatorname{Cov}(N_{i1},N_{i2}) &=E\{\operatorname{Cov}(N_{i1},N_{i2}\mid Z_i)\} +\operatorname{Cov}\{E(N_{i1}\mid Z_i),E(N_{i2}\mid Z_i)\}\\ &=0+\operatorname{Cov}(Z_i a,Z_i b)\\ &=ab\operatorname{Var}(Z_i)=ab\theta. \end{aligned}

For positive $\theta$, this covariance is positive: clusters with larger frailty tend to have more events in both members. This calculation is exact under the stated conditional Poisson assumptions. For a more general recurrent-event intensity, dependence still arises from the shared latent multiplier, but the Poisson count formula need not apply.

Worked Example: A Gamma-Frailty Cluster Count

Suppose two hypothetical patients in one clinic each contribute 2 years at risk. Assume a constant baseline intensity of $0.20$ events per patient-year, no covariate effects, and conditionally independent homogeneous Poisson event processes given the clinic frailty. The combined conditional mean count is $Z_i A$, where $$ A=2(2)(0.20)=0.80. $$ Assume $Z_i\sim\operatorname{Gamma}(2,2)$, so its mean is 1 and variance is $1/2=0.5$. Three events are observed across the two patients.

Conditional on $Z_i=z$, the total count has a Poisson probability. Integrating over the gamma frailty gives the marginal count probability $$ P(N_i=3)= \frac{\Gamma(2+3)}{\Gamma(2)3!} \left(\frac{2}{2+0.80}\right)^2 \left(\frac{0.80}{2+0.80}\right)^3. $$

Substituting and calculating, $$ P(N_i=3)=4(0.7143)^2(0.2857)^3\approx 0.0476. $$

Under these assumptions, the probability of exactly three total events in this clinic is about 0.0476. This is a marginal probability for the count in this particular two-patient, four-patient-year observation setup, not a probability that the frailty takes a particular value.

The extra variation can also be checked directly. The marginal mean count is $E(N_i)=A=0.80$, while the variance is $$ \operatorname{Var}(N_i)=E(Z_iA)+\operatorname{Var}(Z_iA) =A+\theta A^2=0.80+0.5(0.80)^2=1.12. $$

The variance exceeds the mean because clinics with high frailty contribute disproportionately many events. Without frailty, a Poisson count with mean 0.80 would have variance 0.80. The calculation depends on the homogeneous conditional Poisson assumptions and should not be transferred unchanged to a history-dependent event model.

Specification, Alternatives, and Failure Modes

Before fitting a shared frailty model, state the event process and time scale precisely. For example, time since enrollment and time since the previous event define different clocks and generally different baseline intensity functions and risk histories. An estimate from one clock cannot be interpreted as though it came from the other. Likewise, if the target concerns all recurrent events, a model that stops counting after the first event targets a different process.

The main alternative emphasized in the earlier material on Frailty Versus Cluster-Robust Standard Errors is a marginal recurrent-event analysis with cluster-robust standard errors. The robust approach can account for within-cluster correlation in uncertainty estimates without specifying a latent frailty distribution. A shared frailty model instead provides a conditional intensity interpretation and explicitly models residual cluster heterogeneity, but relies on its distributional and independence assumptions. Neither method makes the other’s estimand identical: comparisons require attention to the conditional or marginal target and to the specified event process.

Warning signs include a time scale that is not clearly defined, unmodeled predictors of both censoring and event occurrence, event-history effects omitted despite a scientific reason to expect them, or strong differences between clusters that the covariates and frailty distribution do not represent adequately. Model diagnostics should assess whether fitted event patterns are compatible with observed recurrent-event and cluster patterns; a fitted frailty variance alone does not establish adequate fit. Formal residual and calibration diagnostics are a separate topic.

Mistake: treating $Z_i$ as an observed cluster characteristic. Frailty is latent in this model. Correction: describe it as an unobserved multiplier, and do not assign a measured interpretation to it without additional data.
Mistake: interpreting $\exp(\beta)$ as a marginal effect. Correction: it is a conditional intensity ratio at the same time, covariates, risk status, and frailty value.
Mistake: assuming frailty fixes a wrong time scale or event definition. Correction: specify those features first; frailty models heterogeneity around the chosen process, not a substitute for defining it.
Mistake: concluding that a robust standard error and frailty model answer the same question. Correction: robust standard errors address uncertainty under clustering; the frailty model additionally imposes a latent heterogeneity structure and yields conditional effects.

Key takeaway: A shared frailty multiplies the recurrent-event intensity for every modeled member of a cluster. Its variance describes latent cluster heterogeneity, while the regression coefficients describe covariate contrasts conditional on that frailty and on the specified time scale, risk process, and event-history model.

Check Your Understanding

Answer each question, then compare with the short solutions below.

In a shared frailty model, what does $\exp(\beta)$ represent for a binary covariate?
For mean-one frailty variance $\theta=0.09$, what is the frailty coefficient of variation?
In the conditional Poisson count setting, take $a=1$, $b=2$, and $\theta=0.25$. What is the marginal covariance between the two counts?
Why does changing from time since enrollment to time since the previous event change the interpretation of the fitted model?

Solutions. (1) It is the conditional intensity ratio comparing covariate values 1 and 0, holding the frailty, time, risk status, and other modeled features fixed. (2) The coefficient of variation is $\sqrt{0.09}=0.30$, or 30%. (3) The covariance is $ab\theta=(1)(2)(0.25)=0.50$ expected-count units squared. (4) The clock changes the time coordinate, baseline intensity, and event-history/risk representation; therefore the model describes a different conditional event-rate contrast unless the estimand is explicitly held fixed by appropriate specification.