Cholesky Whitening for Correlated Random Effects

Why whiten correlated random effects?

A Gaussian random-effects vector can have components that vary together. Numerical integration is often more convenient when the integration variables are independent standard normals, because standard-normal quadrature rules can then be used directly. Cholesky whitening makes this change of coordinates without changing the integral being evaluated. It changes the coordinates, not the statistical model or the dimension of the integral.

Suppose \(y\) denotes the observed responses, \(b\) is a \(q\)-dimensional Gaussian random-effects vector, and \(\beta\) denotes fixed-effect parameters. Conditional on \(b\), the responses follow the specified GLMM response model. The random effects have mean zero and covariance matrix \(\Sigma\). The marginal likelihood contribution is obtained by integrating the conditional likelihood over the distribution of \(b\). The discussion below concerns how to transform that integral when \(\Sigma\) is correlated.

Positive definiteness matters: it guarantees an ordinary nonsingular Cholesky factor and an invertible change of variables. A zero-mean Gaussian vector with a positive-definite covariance has density \( \phi_q(b;0,\Sigma)=(2\pi)^{-q/2}|\Sigma|^{-1/2}\exp\{-b^\mathsf{T}\Sigma^{-1}b/2\} \). The factor \(L\) is not generally diagonal. Consequently, \(b=Lz\) mixes the independent coordinates of \(z\) to produce the intended covariance among the components of \(b\).

How the change of coordinates preserves the integral

Write \(g(b)=p(y\mid b;\beta)\) for the conditional likelihood at a given random-effects vector. For fixed model parameters, the marginal likelihood contribution is \( \int_{\mathbb R^q}g(b)\phi_q(b;0,\Sigma)\,db \). Under \(b=Lz\), the Jacobian determinant is \(|\det L|\), and \(\Sigma^{-1}=L^{-\mathsf{T}}L^{-1}\). These identities give the transformed integral exactly.

For the cancellation in the second line, \(|\Sigma|=|LL^\mathsf{T}|=(\det L)^2\). Since the diagonal entries of \(L\) are positive, \(|\det L|=\det L=|\Sigma|^{1/2}\). The factor from the change-of-variables Jacobian therefore cancels the covariance determinant in the Gaussian density. The exponent becomes \(z^\mathsf{T}z\), since \(b^\mathsf{T}\Sigma^{-1}b=z^\mathsf{T}z\). These steps establish that the transformed variables are independent standard normals and that no extra determinant multiplier remains in the expectation.

The equality is exact; numerical quadrature is used only to approximate the resulting expectation. In a GLMM, the integrand evaluated at a quadrature point \(z\) is the conditional likelihood at \(b=Lz\). If the covariance depends on model parameters, its Cholesky factor must be calculated from the current covariance matrix when that matrix changes.

Using the whitened coordinates for quadrature

Once the integral is written as \(\mathbb E\{g(LZ)\}\), quadrature points are generated in standard-normal coordinates. For each point, first compute \(b=Lz\); then evaluate the conditional likelihood using that \(b\). The Gaussian density weighting is already represented by the standard-normal quadrature weights. Multiplying by the original correlated density again, or adding a further determinant factor, would count the density incorrectly.

The following example uses a two-point standard-normal quadrature rule with nodes \(-1\) and \(1\), each with weight \(1/2\). Its two-dimensional product rule has four points, each with weight \(1/4\). This is a deliberately small illustrative rule; it shows the coordinate transformation and arithmetic, not a guarantee of high accuracy.

Numerical checks and failure modes

Whitening is algebraically straightforward, but implementation errors can produce plausible-looking yet incorrect likelihoods. The most important checks concern the factorization, the direction of the transformation, and the numerical stability of the covariance matrix.

• Check the factor orientation. With a lower-triangular factor, verify \(\Sigma=LL^\mathsf{T}\), not \(L^\mathsf{T}L\). The two products generally differ.
• Check the direction of the map. To generate random effects from standard-normal nodes use \(b=Lz\). To whiten a given \(b\), solve \(Lz=b\); do not multiply by \(L\).
• Check covariance validity. A positive-definite covariance must have positive Cholesky diagonal entries. Failure of the factorization can indicate an invalid covariance matrix, a matrix that is only positive semidefinite, or numerical roundoff near the boundary.
• Check the numerical residual. After factorization, examine the size of \(LL^\mathsf{T}-\Sigma\). When converting a point back, check that \(Lz-b\) is near zero relative to the scale of the entries.
• Check quadrature convergence separately. A correct whitening transformation does not make a low-order quadrature rule accurate automatically. Compare marginal likelihood estimates across increasing orders or appropriate integration rules.

A covariance matrix with a very high correlation can be positive definite but nearly singular. Then the Cholesky factor has a small diagonal entry, and solving for whitened coordinates can magnify rounding errors. For example, if \( \Sigma=\begin{pmatrix}1&0.99\\0.99&1\end{pmatrix} \), then \( L=\begin{pmatrix}1&0\\0.99&\sqrt{0.0199}\end{pmatrix} \), with \(\sqrt{0.0199}\approx0.1411\). For \(b=(0.5,0.4)^\mathsf{T}\), solving \(Lz=b\) gives \(z_1=0.5\) and \(z_2=(0.4-0.99(0.5))/0.1411\approx-0.6734\). Substitution returns a second component close to \(0.4\), but the division by the relatively small diagonal entry makes accurate arithmetic important.

How this approach compares with alternatives

One alternative is to integrate directly in the original \(b\)-coordinates, retaining the correlated Gaussian density explicitly. That is mathematically valid, but a standard independent-normal quadrature rule no longer applies without incorporating the covariance and density carefully. Cholesky whitening makes the standard-normal weight explicit and leaves the conditional likelihood as the function evaluated at transformed nodes.

An eigenvalue factorization can also express a covariance as \(AA^\mathsf{T}\), followed by \(b=Az\). For a positive-definite matrix it provides a valid transformation, although the factor need not be triangular. Cholesky is a conventional choice with a unique lower-triangular factor when the diagonal is required to be positive. If a covariance is exactly rank deficient, an ordinary invertible Cholesky transformation is unavailable; a lower-dimensional representation using only the nonzero covariance directions may be appropriate, but it changes the dimension of integration and requires handling that rank reduction deliberately.

Cholesky whitening does not reduce the number of integration dimensions, remove dependence from the conditional likelihood, or ensure that a chosen quadrature rule is accurate. Tensor-product and sparse-grid quadrature remain alternatives for approximating the standard-normal expectation; the transformation supplies the appropriate coordinates, not a cure for high-dimensional cost. The chosen ordering of random-effect components can affect the triangular factor and therefore the arrangement of finite quadrature points, even though it does not change the exact integral. For an approximate calculation, checking sensitivity to quadrature order—and, when relevant, coordinate ordering—is useful.