Cholesky Transformation for Multivariate Quadrature

Why Transform the Random Effects?

In a joint model, the subject-specific longitudinal and event-time contributions are evaluated conditional on latent random effects. The marginal likelihood contribution for a subject is obtained by integrating that conditional contribution over the random-effects distribution. When the random effects are correlated, the integral is multivariate and its Gaussian density contains the covariance matrix and its inverse. A Cholesky transformation removes those terms from the integration measure, so standard-normal quadrature can be used directly.

The transformation is a change of coordinates, not a new statistical model or an approximation by itself. It preserves the integral exactly when the change of variables and its Jacobian are handled correctly. Numerical quadrature is applied afterward, and its accuracy still depends on the integrand and the integration rule.

Definition: Let \(b_i\) be a \(q\)-dimensional random-effects vector with \(b_i\sim N_q(0,D)\), where \(D\) is symmetric positive definite. Its Cholesky factor is the unique lower-triangular matrix \(L\) with positive diagonal entries such that \(D=LL^\mathsf{T}\). Setting \(b_i=Lz_i\), where \(z_i\sim N_q(0,I_q)\), represents the correlated random effects as a linear transformation of independent standard-normal variables.

Here \(i\) indexes subjects, \(q\) is the number of random effects, and \(I_q\) is the \(q\)-dimensional identity matrix. The Cholesky factor carries the variances and correlations in \(D\); the transformed vector \(z_i\) has independent components, each with mean zero and variance one. If the random-effects mean is \(m\) rather than zero, the corresponding transformation is \(b_i=m+Lz_i\).

Deriving the Standard-Normal Integral

Write \(f_i(b_i;\theta)\) for the subject-specific conditional likelihood contribution, with all observed longitudinal and event-time data for subject \(i\) held fixed. The marginal contribution under the zero-mean Gaussian random-effects model is

L_i(\theta)=\int_{\mathbb{R}^q} f_i(b_i;\theta)\, (2\pi)^{-q/2}|D|^{-1/2} \exp\!\left(-\frac{1}{2}b_i^\mathsf{T}D^{-1}b_i\right)\,db_i.

Because \(D=LL^\mathsf{T}\), we have \(D^{-1}=(L^\mathsf{T})^{-1}L^{-1}\), \(|D|=(\det L)^2\), and \(\det L>0\). Under \(b_i=Lz_i\), the quadratic form becomes \(b_i^\mathsf{T}D^{-1}b_i=z_i^\mathsf{T}z_i\), while the volume element changes by \(db_i=(\det L)\,dz_i\). Thus the factor \(|D|^{-1/2}\) cancels the Jacobian \(\det L\). The result is

L_i(\theta)=\int_{\mathbb{R}^q} f_i(Lz_i;\theta)\, (2\pi)^{-q/2}\exp\!\left(-\frac{1}{2}z_i^\mathsf{T}z_i\right)\,dz_i =E_{z_i}\!\left[f_i(Lz_i;\theta)\right], \qquad z_i\sim N_q(0,I_q).

This derivation is the central justification: the transformed integral is an expectation under a standard multivariate normal distribution. It does not require \(f_i\) to be Gaussian. For joint models, \(f_i\) may combine a longitudinal density and an event-time contribution, including the appropriate contribution for censoring or delayed entry; those details are contained in the conditional likelihood.

Worked Example: Factoring a Correlated Covariance Matrix

Suppose a hypothetical subject has a random intercept and random slope with covariance matrix

D=\begin{pmatrix}4&1.2\\1.2&1\end{pmatrix}.

To obtain a lower-triangular factor, write \(L=\begin{pmatrix}a&0\\c&d\end{pmatrix}\). Then \(LL^\mathsf{T}=\begin{pmatrix}a^2&ac\\ac&c^2+d^2\end{pmatrix}\). Matching entries gives \(a=2\), \(c=1.2/2=0.6\), and \(d=\sqrt{1-0.6^2}=\sqrt{0.64}=0.8\). Therefore

L=\begin{pmatrix}2&0\\0.6&0.8\end{pmatrix}, \qquad LL^\mathsf{T}= \begin{pmatrix}4&1.2\\1.2&0.36+0.64\end{pmatrix}=D.

For transformed coordinates \(z=(0.5,-1)^\mathsf{T}\), the corresponding random effects are \(b=Lz=(2(0.5),\,0.6(0.5)+0.8(-1))^\mathsf{T}=(1,-0.5)^\mathsf{T}\). Thus quadrature can evaluate the conditional contribution at \(b=(1,-0.5)^\mathsf{T}\) by evaluating the transformed integrand at the standard-normal point \(z=(0.5,-1)^\mathsf{T}\). The covariance and correlation are represented by \(L\), rather than by altering the standard-normal nodes.

Using the Transformation in Quadrature

Let \(x_k\) and \(w_k\), for \(k=1,\ldots,K\), be nodes and weights for a one-dimensional standard-normal expectation rule, normalized so that the weights sum to one. The transformed multivariate expectation can be approximated by a product rule:

L_i(\theta)\approx \sum_{k_1=1}^{K}\cdots\sum_{k_q=1}^{K} \left(\prod_{r=1}^{q}w_{k_r}\right) f_i\!\left(L \begin{pmatrix}x_{k_1}\\ \vdots\\ x_{k_q}\end{pmatrix};\theta\right).

The vector of nodes is first assembled in standard-normal coordinates and then multiplied by \(L\). The conditional likelihood is evaluated at that resulting random-effects vector. Using nodes intended for a standard-normal expectation avoids separately coding the original Gaussian density and its determinant. If a quadrature convention instead supplies nodes and weights for a different Gaussian weight, its corresponding scaling and normalization must be applied consistently.

Worked Example: A Two-Point Rule and a Quadrature Diagnostic

Continue with \(D=\begin{pmatrix}4&1.2\\1.2&1\end{pmatrix}\) and \(L=\begin{pmatrix}2&0\\0.6&0.8\end{pmatrix}\). Consider the hypothetical integrand \(g(b)=\exp(-b^\mathsf{T}b/2)\). Use the two-point standard-normal rule with nodes \(x_1=-1,x_2=1\) and weights \(w_1=w_2=1/2\). The four equally weighted transformed node vectors are as follows.

At \(z=(1,1)^\mathsf{T}\), \(b=(2,1.4)^\mathsf{T}\), so \(b^\mathsf{T}b=2^2+1.4^2=5.96\) and \(g(b)=\exp(-2.98)\approx0.05079\). At \(z=(1,-1)^\mathsf{T}\), \(b=(2,-0.2)^\mathsf{T}\), so \(b^\mathsf{T}b=4+0.04=4.04\) and \(g(b)=\exp(-2.02)\approx0.1327\). The two remaining points are the negatives of these vectors, so they give the same respective values. The quadrature estimate is therefore

\frac{2\exp(-2.98)+2\exp(-2.02)}{4} =\frac{\exp(-2.98)+\exp(-2.02)}{2} \approx 0.09172.

For comparison, if \(b\sim N_2(0,D)\), direct Gaussian integration gives \(E[\exp(-b^\mathsf{T}b/2)]=|I_2+D|^{-1/2}\). Here \(|I_2+D|=\left|\begin{pmatrix}5&1.2\\1.2&2\end{pmatrix}\right|=10-1.44=8.56\), so the exact expectation is \(1/\sqrt{8.56}\approx0.3418\). The two-point quadrature estimate is far below this value. The transformation and its arithmetic are correct; the discrepancy shows that this low-order rule poorly approximates this integrand. In particular, its nodes do not include the region near zero where much of the integrand’s contribution lies. Increasing quadrature accuracy is a separate decision from applying the Cholesky transformation correctly.

Worked Example: Checking the Transformed Integral Analytically

As another hypothetical calculation, let \(D=\begin{pmatrix}1&0.6\\0.6&1\end{pmatrix}\), with \(L=\begin{pmatrix}1&0\\0.6&0.8\end{pmatrix}\), and take \(f(b)=\exp(t^\mathsf{T}b)\) for \(t=(0.5,-1)^\mathsf{T}\). The Gaussian moment-generating function gives the exact expectation \(E[\exp(t^\mathsf{T}b)]=\exp(t^\mathsf{T}Dt/2)\). First, \(t^\mathsf{T}Dt=0.25+1+2(0.5)(-1)(0.6)=0.65\), so the expectation is \(\exp(0.325)\approx1.384\).

In standard-normal coordinates, \(t^\mathsf{T}b=t^\mathsf{T}Lz=(L^\mathsf{T}t)^\mathsf{T}z\), where \(L^\mathsf{T}t=(-0.1,-0.8)^\mathsf{T}\). Since the components of \(z\) are independent standard normals,

E\!\left[\exp\!\left((-0.1)z_1-0.8z_2\right)\right] =\exp\!\left(\frac{1}{2}\{(-0.1)^2+(-0.8)^2\}\right) =\exp(0.325)\approx1.384.

The two calculations agree because \(\|L^\mathsf{T}t\|^2=t^\mathsf{T}LL^\mathsf{T}t=t^\mathsf{T}Dt\). This check illustrates that the transformation preserves the covariance structure even though the integration variables are independent.

Assumptions, Diagnostics, and Alternatives

For the ordinary Cholesky transformation, \(D\) must be symmetric positive definite. This ensures that a nonsingular lower-triangular \(L\) exists and that the Gaussian density has a valid determinant and inverse. In an estimated model, \(D\) should be checked for symmetry and positive definiteness at each parameter value where the integration is performed. Small numerical asymmetry can be corrected by replacing \(D\) with \((D+D^\mathsf{T})/2\), but a materially negative eigenvalue is not merely a rounding issue: it signals an invalid covariance estimate or computation.

Useful implementation checks include verifying that \(LL^\mathsf{T}\) reconstructs \(D\) to numerical tolerance, confirming that transformed standard-normal draws \(Lz\) have empirical covariance close to \(D\), and checking that the transformation uses \(Lz\), not \(L^\mathsf{T}z\). The latter error generally produces covariance \(L^\mathsf{T}L\), which need not equal \(D\). Also check the integration weights: standard-normal weights should sum to one, and applying an additional Gaussian density or determinant factor after using the expectation rule would count the density twice.

Cholesky whitening is particularly convenient when covariance parameters change during likelihood optimization: each current \(D\) is factored, and the same standard-normal nodes can be transformed by its current \(L\). It is not the only way to map a Gaussian vector to independent standard normals. An eigenvalue decomposition can also construct a square root, and may be useful for diagnosing nearly singular matrices. For a positive-definite matrix, however, the Cholesky factor is triangular and usually efficient. Direct integration in the original coordinates is mathematically equivalent, but requires careful handling of the correlated density, its inverse covariance, and its normalizing constant.

When a quadrature rule is adapted around a posterior mode, the proposal’s center and scale are part of that rule; if the proposal differs from the random-effects prior, its density must be accounted for in the integral. The Cholesky factor of \(D\) alone does not supply that adjustment. The transformation also does not reduce the dimension \(q\) or guarantee accurate quadrature. It improves the coordinates in which Gaussian integration is carried out, while approximation error and computational burden remain.

Common Pitfall

Common mistake: “The random effects are correlated, so I can use independent normal quadrature nodes directly as values of \(b_i\).” This drops the covariance structure. The nodes are values of \(z_i\), not \(b_i\); map each node vector through \(b_i=Lz_i\) before evaluating the conditional contribution. Conversely, if the integrand is written in \(z_i\)-coordinates, do not multiply by the original Gaussian density and determinant again: the Jacobian derivation has already converted the measure to the standard-normal density.

The practical importance is not limited to computational convenience. A wrong factor orientation, missing Jacobian, or duplicated density changes the marginal likelihood and can therefore affect parameter estimates and uncertainty summaries. A reliable implementation checks the factorization and transformation independently of the quadrature approximation, then assesses whether the chosen integration rule is sufficiently accurate for the likelihood at hand.

Check Your Understanding

Answer each question, then compare with the solutions.

If \(D=LL^\mathsf{T}\) and \(z\sim N_q(0,I_q)\), what is the covariance matrix of \(b=Lz\)?
For \(L=\begin{pmatrix}2&0\\0.6&0.8\end{pmatrix}\) and \(z=(-1,1)^\mathsf{T}\), calculate \(b\).
Why does the determinant factor in the Gaussian density cancel under the Cholesky change of variables?
If a transformed quadrature result is inaccurate, does that alone show that the Cholesky factor is wrong? Give one separate diagnostic for the factor and one for quadrature accuracy.

Solutions. (1) \(\operatorname{Cov}(b)=L\operatorname{Cov}(z)L^\mathsf{T}=LI_qL^\mathsf{T}=D\). (2) \(b=(2(-1),\,0.6(-1)+0.8(1))^\mathsf{T}=(-2,0.2)^\mathsf{T}\). (3) Since \(|D|=(\det L)^2\) and \(\det L>0\), \(|D|^{-1/2}=1/\det L\), which cancels the Jacobian \(db=(\det L)\,dz\). (4) No. Check whether \(LL^\mathsf{T}\) reconstructs \(D\) for the factor; assess quadrature accuracy by increasing the number of nodes and comparing the resulting integral estimates.