Likelihood Integration for Lognormal Frailty

Why integrate over the random effect?

A cluster-level frailty makes event rates within a cluster conditionally related: given the frailty, the event process has a specified intensity, but the frailty itself is unobserved. The likelihood for observed data must therefore account for uncertainty about that latent effect. For a lognormal frailty, this means integrating the conditional likelihood against a Gaussian density. The resulting integral is generally not available in closed form, so numerical quadrature is a practical way to evaluate it.

The integration is over the random effect, not over the observed event times. Once the conditional likelihood has been written, the likelihood contribution for a cluster is its expectation over the random-effect distribution. Keeping this distinction clear helps avoid treating an unobserved frailty as if it were an observed covariate or multiplying conditional likelihoods evaluated at an arbitrary frailty value.

Definition: For cluster \(i\), let \(N_i\) be its observed event count and let \(A_i\) be its cumulative baseline intensity over the observation period, including any modeled covariate multiplier. Suppose that, conditional on a shared frailty \(Z_i\), the event process is Poisson with intensity multiplied by \(Z_i\). Its likelihood contribution can be written \(L_i(Z_i)=C_i Z_i^{N_i}\exp(-A_iZ_i)\), where \(C_i\) is a factor that does not depend on \(Z_i\). Under mean-one lognormal frailty, write \(Z_i=\exp(B_i)\), where \(B_i\sim N(-\sigma^2/2,\sigma^2)\). The marginal cluster likelihood contribution is \(L_i=\mathbb{E}\{L_i(\exp(B_i))\}\).

From conditional likelihood to a Gaussian integral

The mean-one convention is useful because it separates the overall intensity scale from the amount of heterogeneity. If \(B_i\sim N(-\sigma^2/2,\sigma^2)\), then the lognormal moment formula gives \(\mathbb{E}[\exp(B_i)]=\exp(-\sigma^2/2+\sigma^2/2)=1\). The variance parameter \(\sigma^2\) controls variation in frailty; when \(\sigma^2=0\), the frailty is identically one.

Substitute \(Z_i=\exp(B_i)\) into the conditional likelihood and integrate with respect to the density of \(B_i\). Equivalently, set \(B_i=-\sigma^2/2+\sigma U_i\), where \(U_i\sim N(0,1)\). The latter form exposes a standard normal expectation that can be approximated with familiar Gaussian quadrature.

L_i=C_i\int_{-\infty}^{\infty} \exp\{N_i b-A_i\exp(b)\} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{(b+\sigma^2/2)^2}{2\sigma^2}\right\}\,db

L_i=C_i\,\mathbb{E}\left[ \exp\left\{N_i\left(-\frac{\sigma^2}{2}+\sigma U_i\right) -A_i\exp\left(-\frac{\sigma^2}{2}+\sigma U_i\right)\right\} \right],\qquad U_i\sim N(0,1).

For independent clusters, conditional independence given the cluster frailties and independence of the frailties across clusters imply that the marginal likelihood is the product of these integrated contributions: \(\mathcal{L}=\prod_i L_i\). If the model includes event-time factors, they belong in \(C_i\) when they do not depend on \(B_i\); if the event intensity has a different frailty-dependent form, the conditional likelihood kernel must be changed before integration.

Quadrature rule: For a smooth function \(g\), the three-node Gauss–Hermite rule for a standard normal variable is \(\mathbb{E}\{g(U)\}\approx \frac{1}{6}g(-\sqrt{3})+\frac{2}{3}g(0)+\frac{1}{6}g(\sqrt{3})\). For a lognormal frailty random effect, evaluate the integrand at \(b_q=-\sigma^2/2+\sigma u_q\), where \(u_q\) are these three standard-normal nodes. The weights are \(1/6,2/3,1/6\), respectively.

Worked examples

Worked Example: Building the Cluster Contribution

Suppose a hypothetical cluster has two observed events and cumulative baseline intensity \(A_i=1.4\). Ignore factors independent of frailty, so \(C_i=1\), and suppose \(B_i\sim N(-0.125,0.25)\). Write its marginal likelihood as a one-dimensional integral.

Here \(N_i=2\), \(\sigma^2=0.25\), and \(\sigma=0.5\). Conditional on \(B_i=b\), the likelihood kernel is \(\exp\{2b-1.4\exp(b)\}\). The density of \(B_i\) has mean \(-0.125\) and standard deviation \(0.5\). Therefore,

L_i=\int_{-\infty}^{\infty} \exp\{2b-1.4\exp(b)\} \frac{1}{\sqrt{2\pi(0.25)}} \exp\left\{-\frac{(b+0.125)^2}{2(0.25)}\right\}\,db.

Equivalently, this is the standard normal expectation of \(\exp\{2(-0.125+0.5U)-1.4\exp(-0.125+0.5U)\}\). This expression identifies precisely what must be integrated. It is not the likelihood conditional on setting \(B_i\) to its mean: because the likelihood kernel is nonlinear in \(b\), evaluating it at \(b=-0.125\) does not generally equal its average over the Gaussian distribution.

Worked Example: Three-Node Quadrature for Two Clusters

Consider a hypothetical pair of independent clusters under the same mean-one frailty distribution, \(B_i\sim N(-0.125,0.25)\). Cluster 1 has \(N_1=1\) and \(A_1=1\); cluster 2 has \(N_2=0\) and \(A_2=1\). Take \(C_1=C_2=1\). Approximate the joint marginal likelihood contribution with three-node quadrature.

The standard normal nodes are \(-\sqrt{3},0,\sqrt{3}\), approximately \(-1.7321,0,1.7321\). Transforming them to the random-effect scale gives \(b_1=-0.9910\), \(b_2=-0.1250\), and \(b_3=0.7410\). The weights remain \(1/6,2/3,1/6\).

For cluster 1, the kernel is \(g_1(b)=\exp\{b-\exp(b)\}\). At the transformed nodes its values are approximately \(0.2561, 0.3651,\) and \(0.2574\). Thus,

L_1\approx \frac{1}{6}(0.2561)+\frac{2}{3}(0.3651)+\frac{1}{6}(0.2574) \approx 0.3290.

Using the unrounded node evaluations, the quadrature contribution is approximately \(0.3290057\). The displayed four-decimal kernel values are for readability; retaining more precision in the calculation avoids small rounding discrepancies.

For cluster 2, the kernel is \(g_2(b)=\exp\{-\exp(b)\}\). Its values at the same nodes are approximately \(0.6899, 0.4137,\) and \(0.1227\). Hence,

L_2\approx \frac{1}{6}(0.6899)+\frac{2}{3}(0.4137)+\frac{1}{6}(0.1227) \approx 0.4113.

The clusters are independent after marginalizing their independent frailties, so the joint contribution is the product of the two integrated contributions. Using the unrounded quadrature contributions gives \(\mathcal{L}\approx 0.3290057\times L_2\approx 0.1353087\), or \(0.1353\) to four significant figures. Multiplying rounded intermediate values may shift the last displayed digit, so the unrounded values should be used for the final product.

Worked Example: Checking the No-Frailty Limit

Use cluster 1 from the previous example, with one event and \(A_1=1\), but set \(\sigma^2=0\). Then \(B_1=0\) and \(Z_1=1\) with probability one. The conditional likelihood is \(\exp\{1(0)-1\exp(0)\}=\exp(-1)\approx0.3679\).

This is the correct zero-variance limit of the integral: as \(\sigma^2\) tends to zero, the Gaussian distribution of \(B_1\) concentrates at zero, so averaging the continuous likelihood kernel approaches its value at zero. It also provides a basic implementation check. A quadrature routine with zero variance should return the no-frailty likelihood, not an undefined expression or a contribution evaluated at an unrelated node.

Why the integral is justified, and how to assess it

The marginalization follows directly from the law of total probability for a latent variable. Given the random effect \(B_i=b\), the data have likelihood \(L_i(b)\). Averaging that conditional likelihood over the probability distribution of \(B_i\) gives the probability model for the observed data alone. In density notation, this is the integral of the conditional data density times the random-effect density.

Specify the conditional model.
State the event intensity and form its likelihood conditional on the cluster frailty. Confirm that the frailty multiplier enters the event term and the cumulative intensity consistently.

Specify the random-effect distribution.
For mean-one lognormal frailty, use \(B_i\sim N(-\sigma^2/2,\sigma^2)\), and state that independent clusters have independent random effects when the model assumes this.

Integrate, then approximate if needed.
Write the Gaussian integral before applying a numerical rule. Transform standard-normal quadrature nodes onto the \(B_i\) scale and apply their weights to the full conditional likelihood kernel.

Check numerical stability.
Compare results using more quadrature nodes or a different numerical integration setting, and check limiting cases such as \(\sigma^2=0\). Material changes indicate that the initial approximation may be inadequate.

The conditional model must itself be appropriate. For the Poisson-process example, the assumed intensity and observation period must describe the event process, and \(A_i\) must correctly account for exposure, covariates, and the observation scheme. Conditional on the frailty and modeled covariates, the likelihood must match the assumed event process. If censoring or observation is informative in a way not represented in the likelihood, numerical integration cannot repair that modeling problem.

Quadrature is an approximation to an expectation, not a change to the statistical model. Three-node Gauss–Hermite quadrature is exact for polynomial functions up to a specified degree under its weighting rule, but the lognormal likelihood kernel is not generally such a polynomial. Accuracy should therefore be checked by increasing the number of nodes and comparing the resulting likelihoods or log-likelihoods. For large likelihoods or many clusters, computing sums of log contributions can also help prevent numerical underflow.

Alternatives, diagnostics, and limitations

A useful comparison is with gamma frailty. In some gamma-frailty models, integrating the conditional likelihood can be reduced to a closed-form expression using the gamma distribution's Laplace transform. The lognormal distribution has no corresponding elementary Laplace transform for this likelihood kernel, so numerical integration is typically needed. Direct integration on the \(B_i\) scale is natural because the Gaussian density is explicit; integrating on the \(Z_i\) scale is also possible, but requires the lognormal density and careful handling of its support \(Z_i>0\).

The single-cluster integral is one-dimensional, but the joint integral over many independent cluster effects would be high-dimensional if formed as one integral. The factorization into cluster contributions avoids that: integrate each cluster's one-dimensional conditional likelihood separately, then multiply the results (or add their log contributions). This factorization relies on the model's independence assumptions; correlated or nested random effects require an integral reflecting their joint distribution.

Signs that an evaluation may be going wrong include transformed nodes that omit the mean shift \(-\sigma^2/2\), quadrature weights that do not sum to one, a likelihood contribution outside its possible range in a probability model, or an estimate that changes substantially as nodes are added. A common implementation mistake is to apply standard-normal nodes directly as values of \(B_i\). Here the correct nodes are \(b_q=-\sigma^2/2+\sigma u_q\); using \(u_q\) without this transformation evaluates the integrand under the wrong random-effect distribution.

Numerical integration can be computationally demanding when a likelihood must be evaluated repeatedly across many parameter values. A fixed low-order rule may be inaccurate if the likelihood kernel is sharply concentrated or if the variance is large. Conversely, a very large number of nodes may waste computation for a smooth, well-behaved integrand. The appropriate choice is supported by convergence checks, not by assuming one node count works for every cluster and parameter setting.

Key takeaway: Under lognormal frailty, a marginal cluster likelihood is a Gaussian average of the conditional likelihood evaluated at the exponentiated random effect. Write that integral first; then transform the nodes, apply the weights, and verify the approximation.

Check Your Understanding

Answer each question, then compare with the solutions below.

For mean-one lognormal frailty with \(B\sim N(-\sigma^2/2,\sigma^2)\), show why \(\mathbb{E}[\exp(B)]=1\).
A cluster has \(N=0\), \(A=2\), and \(C=1\). Write its conditional likelihood kernel as a function of \(b\).
With \(\sigma=0.5\), what is the transformed random-effect node corresponding to the standard-normal node \(u=0\)?
Why is the product of cluster-specific marginal contributions appropriate for independent clusters with independent frailties?

Solutions. (1) The lognormal moment formula gives \(\mathbb{E}[\exp(B)]=\exp\{-\sigma^2/2+\sigma^2/2\}=1\). (2) Substituting \(N=0,A=2,C=1\) into \(C\exp\{Nb-A\exp(b)\}\) gives \(\exp\{-2\exp(b)\}\). (3) The transformed node is \(b=-\sigma^2/2+\sigma u=-0.125\). (4) Independence makes the joint distribution of cluster data and random effects factor into cluster-specific terms; integrating each independent effect therefore yields a product of its marginal likelihood contributions.