EM Estimation With Latent Frailties

Why Treat Frailty as Missing Data?

A gamma frailty model represents unobserved differences in event propensity by a positive latent variable. For a cluster of recurrent-event observations, a shared frailty makes the cluster’s events conditionally more or less frequent together. The likelihood can be integrated over the frailty distribution, as in direct likelihood methods, but this integration can be difficult when model parameters must be updated repeatedly. The EM algorithm instead treats each cluster’s frailty as missing data: it calculates posterior frailty moments under the current parameters, then maximizes the expected complete-data log likelihood.

This tutorial uses a shared-frailty recurrent-event model with a nonparametric baseline cumulative hazard. The likelihood is a conditional counting-process likelihood; it is not a claim that recurrent events are independent after marginalizing over frailty.

Model and assumptions: For independent clusters $i=1,\ldots,n$, let $N_i(t)$ count recurrent events and let $Y_i(t)$ indicate that cluster $i$ is under observation and at risk. Given a cluster-specific frailty $Z_i$ and the observed covariate process $X_i(t)$, assume the event intensity is $$ \lambda_i(t\mid Z_i)=Y_i(t)Z_i\exp\{\beta^\top X_i(t)\}\lambda_0(t). $$ The frailties are independent across clusters and $Z_i\sim \operatorname{Gamma}(r,r)$, using shape and rate parameterization. Thus $E(Z_i)=1$, $\operatorname{Var}(Z_i)=1/r=\theta$, and $r=1/\theta$. Assume the conditional event-process likelihood is valid for the specified intensity, covariates are predictable, and censoring is noninformative conditional on the modeled information. The model does not include a separate likelihood for the censoring process.

Write $d_i=N_i(\tau_i)$ for the number of observed events in cluster $i$ through its observation end $\tau_i$, and define its cumulative exposure, at current regression and baseline-hazard parameters, by $$ A_i(\beta,\Lambda_0)=\int_0^{\tau_i}Y_i(t)\exp\{\beta^\top X_i(t)\}\,d\Lambda_0(t). $$ Ignoring terms that do not depend on the unknown parameters, the conditional likelihood contribution given $Z_i$ is proportional to $Z_i^{d_i}\exp(-Z_iA_i)$ times the usual event and covariate terms. This representation is the basis of the E-step.

The E-Step: Posterior Frailty Moments

At an iteration with current frailty shape $r$, the prior density is proportional to $z^{r-1}\exp(-rz)$. Multiplying it by the cluster likelihood terms involving $z$ gives a gamma posterior. Let $m_i$ denote its posterior mean and $\ell_i$ its posterior expected log frailty.

Posterior update: Given the observed events and exposure, $Z_i$ has a gamma posterior with shape $a_i=r+d_i$ and rate $b_i=r+A_i$. Therefore, $$ m_i=E(Z_i\mid \text{data})=\frac{r+d_i}{r+A_i}, \qquad \ell_i=E(\log Z_i\mid \text{data})=\psi(r+d_i)-\log(r+A_i), $$ where $\psi$ is the digamma function, the derivative of $\log\Gamma(\cdot)$.

Worked Example: Calculating Posterior Frailty Moments

Suppose the current frailty variance is $\theta=0.5$, so the current shape and rate are $r=1/\theta=2$. For two hypothetical clusters, suppose the current model gives $(d_1,A_1)=(2,1.5)$ and $(d_2,A_2)=(0,0.5)$.

For cluster 1, the posterior shape is $2+2=4$ and the rate is $2+1.5=3.5$. Hence $$ m_1=\frac{4}{3.5}=1.143. $$ For integer shape 4, $\psi(4)=1.2561$, so $$ \ell_1=\psi(4)-\log(3.5)=1.2561-1.2528=0.00335. $$

For cluster 2, the posterior shape is $2+0=2$ and the rate is $2+0.5=2.5$. Thus $$ m_2=\frac{2}{2.5}=0.8000, \qquad \ell_2=\psi(2)-\log(2.5)=0.4228-0.9163=-0.4935. $$

The event-free cluster has a posterior mean below the prior mean of 1, while the cluster with two events has a posterior mean above 1. These are posterior summaries for the E-step, not fixed estimates of individual frailties. The calculation conditions on the stated event counts and exposures and uses the current model parameters.

The M-Step: Updating Regression and Baseline Hazard

The expected complete-data log likelihood replaces the unobserved $Z_i$ and $\log Z_i$ contributions by $m_i$ and $\ell_i$. For fixed posterior moments from the E-step, the terms involving $\beta$ and $\Lambda_0$ have the form $$ Q_{\beta,\Lambda_0} = \sum_i\int_0^{\tau_i}\{\beta^\top X_i(t)+\log d\Lambda_0(t)\}\,dN_i(t) - \sum_i m_i A_i(\beta,\Lambda_0). $$ The frailty-prior terms are handled separately to update $r$. In the M-step, the posterior weights $m_i$ are held fixed; they are not recalculated during that maximization.

Represent the nonparametric baseline hazard by jumps at distinct observed event times $t_k$. At fixed $\beta$, maximizing $Q$ over the jump at $t_k$ gives the update below. The numerator counts events at that time, while the denominator is a risk-set sum weighted by posterior mean frailty and the covariate multiplier.

Baseline jump update: For $dN_i(t_k)$ events in cluster $i$ at $t_k$, $$ \Delta\widehat{\Lambda}_0(t_k) = \frac{\sum_i dN_i(t_k)} {\sum_i m_iY_i(t_k)\exp\{\beta^\top X_i(t_k)\}}. $$ At fixed baseline jumps, the regression update maximizes $Q_{\beta,\Lambda_0}$. Its score equation is $$ \sum_i\int_0^{\tau_i}X_i(t)\,dN_i(t) - \sum_i m_i\int_0^{\tau_i}Y_i(t)\exp\{\beta^\top X_i(t)\}X_i(t)\,d\Lambda_0(t) =0. $$

In general, the regression coefficients and baseline jumps are updated together by maximizing $Q$, or by iterating conditional maximizations until the expected complete-data log likelihood no longer increases. If a computational routine only takes steps that increase $Q$, rather than finding its exact maximum, the procedure is a generalized EM algorithm.

Worked Example: Updating One Baseline-Hazard Jump

At a hypothetical event time $t_k$, suppose three clusters are at risk. Their current posterior frailty means are $1.2$, $0.8$, and $1.0$, and their current covariate multipliers $\exp\{\beta^\top X_i(t_k)\}$ are $1$, $2$, and $1.5$. Suppose clusters 1 and 3 each have one event at $t_k$, and cluster 2 has none. For this update, the covariates, risk indicators, regression coefficients, and posterior means are treated as fixed.

The numerator is the total number of events at $t_k$, namely $1+0+1=2$. The denominator is $$ 1.2(1)+0.8(2)+1.0(1.5)=1.2+1.6+1.5=4.3. $$ Therefore, $$ \Delta\widehat{\Lambda}_0(t_k)=\frac{2}{4.3}=0.4651. $$

This is a baseline cumulative-hazard jump, not a hazard rate per unit time. Its value depends on the risk set at $t_k$ and on the current posterior frailty and regression updates. The update uses two observed events in the numerator; it does not multiply each event by its cluster’s posterior mean.

Updating the Frailty Variance

The gamma prior also contributes to the expected complete-data log likelihood. Let $r=1/\theta$. For one cluster, its log density, apart from no omitted parameter-dependent terms, is $$ r\log r-\log\Gamma(r)+(r-1)\log Z_i-rZ_i. $$ After taking the E-step expectation and summing over clusters, the part of $Q$ involving $r$ is $$ Q_r=n\{r\log r-\log\Gamma(r)\}+(r-1)\sum_i\ell_i-r\sum_i m_i. $$ Differentiating with respect to $r$, and setting the derivative to zero, gives $$ \log r+1-\psi(r)+\frac{1}{n}\sum_i\ell_i-\frac{1}{n}\sum_i m_i=0. $$ The E-step values in this equation are held fixed while maximizing over $r>0$. There is generally no closed-form solution, so a one-dimensional numerical maximization or root-finding procedure is used; then the updated variance is $\theta=1/r$.

Worked Example: Solving the Frailty-Variance Update

Use the two clusters from the first example. Their E-step summaries give $$ \bar m=\frac{1.1429+0.8000}{2}=0.9714, \qquad \bar\ell=\frac{0.00335-0.4935}{2}=-0.2451. $$ The shape update solves $$ f(r)=\log r+1-\psi(r)+\bar\ell-\bar m=0. $$ Substitution gives $\bar\ell-\bar m=-1.2165$, so a numerical solver finds $r$ approximately $2.46$. The updated frailty variance is consequently $$ \widehat{\theta}_{\mathrm{new}}=\frac{1}{2.46}\approx 0.406. $$

The root is a maximum of $Q_r$: its derivative with respect to $r$ is $f(r)$, and $$ f'(r)=\frac{1}{r}-\psi_1(r)<0, $$ where $\psi_1$ is the trigamma function and $\psi_1(r)>1/r$ for $r>0$. Thus the score decreases with $r$, and an interior root is unique. This update is conditional on the E-step summaries; it does not use newly calculated moments until the next E-step.

Why the EM Steps Increase the Likelihood

The EM principle can be justified through the conditional distribution of the missing data. Let $\eta$ denote all model parameters, and let $q_{\text{old}}(Z)$ be the conditional distribution of the frailties given the observed data under the current parameter value $\eta_{\text{old}}$. For any candidate $\eta$, the observed log likelihood satisfies $$ \ell_{\mathrm{obs}}(\eta)-\ell_{\mathrm{obs}}(\eta_{\mathrm{old}}) = Q(\eta\mid\eta_{\mathrm{old}})-Q(\eta_{\mathrm{old}}\mid\eta_{\mathrm{old}}) + \operatorname{KL}\{q_{\text{old}}\|q_{\eta}\}. $$ Here $q_\eta$ is the frailty distribution conditional on the observed data at $\eta$, and $\operatorname{KL}$ is the nonnegative Kullback–Leibler divergence. Therefore, if the M-step raises $Q$, the observed log likelihood cannot decrease. An exact M-step maximizes $Q$; a generalized EM step need only increase it. This guarantee does not establish that the algorithm finds the global maximum.

Implementation, Alternatives, and Pitfalls

An iteration consists of recalculating $A_i$, $m_i$, and $\ell_i$ using the current parameters; maximizing the expected complete-data objective; and checking convergence. In practice, monitor the observed log likelihood as well as changes in parameter estimates. A substantial decrease in the observed log likelihood suggests an implementation error, an inaccurate numerical maximization, or a mismatch between the likelihood and the assumed observation process.

Check the posterior rate. The rate is $r+A_i$, not $r-A_i$. Greater exposure with no corresponding event evidence should pull the posterior frailty mean downward.
Use both posterior moments. The regression and baseline updates use $m_i$; the gamma prior update also requires $\ell_i$. Substituting $\log(m_i)$ for $\ell_i$ is generally incorrect.
Keep E-step weights fixed in the M-step. Recomputing posterior moments while optimizing the same $Q$ function changes the objective and is not the specified EM update.
Track numerical behavior. Poor starting values, an inaccurate root solve for $r$, or nearly vanishing frailty variance can cause slow or unstable iterations. Verify that each M-step increases $Q$, and verify the observed likelihood independently when possible.

Direct marginal likelihood maximization, including adaptive quadrature, integrates over frailty and optimizes the observed likelihood without these EM iterations. It may converge faster in some settings, but requires accurate likelihood integration and derivatives. EM replaces that integration in the update with posterior moments and often gives transparent updates, but can converge slowly near a maximum. Neither approach repairs an incorrectly specified event-process likelihood or noninformative-censoring assumption. The introductory material on gamma frailty and its Laplace transform gives the distributional foundation; the EM construction here uses the cluster-level conditional likelihood.

Key takeaway: For a mean-one gamma shared-frailty model, the E-step has a closed-form gamma posterior with shape $r+d_i$ and rate $r+A_i$. Its mean weights the event-process M-step, its expected log enters the frailty-prior update, and the updated variance is obtained by solving a one-dimensional digamma equation.

Check Your Understanding

Try each question, then compare your reasoning with the brief solutions.

With current $r=3$, a cluster has $d_i=1$ event and $A_i=2$. What are its posterior shape and rate?
For a fixed regression parameter, why does a cluster’s posterior mean $m_i$ appear in the denominator of a baseline-jump update?
What posterior summaries are required for the frailty-shape update, and why is the posterior mean alone insufficient?
If a numerical M-step raises $Q$ but the observed log likelihood decreases, what should be checked?

Solutions. (1) The posterior is gamma with shape $3+1=4$ and rate $3+2=5$. (2) The expected complete-data exposure term is $-m_iA_i$, so the risk contribution is weighted by the conditional expected frailty. (3) Both $m_i=E(Z_i\mid\text{data})$ and $\ell_i=E(\log Z_i\mid\text{data})$ are needed; the gamma prior log density contains both $Z_i$ and $\log Z_i$. (4) Check that the likelihood, E-step posterior moments, and M-step objective correspond to the same model, and that the likelihood and M-step calculations are numerically accurate.