Transforming Random-Effects Integrals to Standard Normal Coordinates

Why Change the Integration Coordinates?

A GLMM marginal likelihood integrates a conditional likelihood over the distribution of its random effects. For a Gaussian random effect, the original integral is weighted by a Gaussian density whose mean and covariance may vary with the model parameters. Re-expressing that integral using standard-normal coordinates separates the standard Gaussian weight from the model-specific scale and correlation. This is an exact change of variables, not an approximation: it changes how the integral is written, not its value.

That distinction matters when evaluating likelihood integrals numerically. A standard-normal representation gives a common reference distribution even when the random-effect variance changes. But the transformation is only correct if the scale or matrix factor is applied to the argument of the conditional likelihood and its Jacobian is handled correctly when transforming a density-weighted integral.

Definition: Let \(b\) be a \(d\)-dimensional Gaussian random effect with mean vector \(\mu\) and positive-definite covariance matrix \(\Sigma\). Choose a nonsingular matrix \(A\) such that \(\Sigma=AA^\mathsf{T}\), and let \(z\sim N_d(0,I_d)\), where \(I_d\) is the \(d\)-dimensional identity matrix. Then \(b=\mu+Az\) has distribution \(N_d(\mu,\Sigma)\). The vector \(z\) is the standard-normal coordinate for \(b\).

The factor \(A\) expresses the location and covariance scale of \(b\). In one dimension, \(A\) is just the positive standard deviation. In multiple dimensions, the factor also maps the independent coordinates of \(z\) to the covariance structure in \(b\). This tutorial assumes \(\Sigma\) is positive definite, so that such a factor can be chosen nonsingular and the change of variables is invertible.

The Change of Variables and Its Jacobian

The affine transformation is \(b=\mu+Az\), with inverse \(z=A^{-1}(b-\mu)\). Its derivative matrix is \(A\), so the absolute Jacobian determinant is \(|\det A|\). In differential-volume notation, \(db=|\det A|\,dz\). The absolute value is essential: a volume element is nonnegative even if a chosen factor has negative determinant.

The Gaussian density of \(b\) is

f_b(b)=\frac{1}{(2\pi)^{d/2}|\Sigma|^{1/2}} \exp\!\left[-\frac{1}{2}(b-\mu)^\mathsf{T}\Sigma^{-1}(b-\mu)\right].

Because \(\Sigma=AA^\mathsf{T}\), we have \(|\Sigma|^{1/2}=|\det A|\). Substituting \(b=\mu+Az\) also gives \((b-\mu)^\mathsf{T}\Sigma^{-1}(b-\mu)=z^\mathsf{T}z\). Therefore,

f_b(\mu+Az)=\frac{\phi_d(z)}{|\det A|}, \qquad \phi_d(z)=\frac{1}{(2\pi)^{d/2}}\exp\!\left(-\frac{1}{2}z^\mathsf{T}z\right).

For any integrable function \(h\), the ordinary change-of-variables formula now gives

\int_{\mathbb{R}^d} h(b)f_b(b)\,db = \int_{\mathbb{R}^d} h(\mu+Az)\phi_d(z)\,dz.

The density contributes a factor \(1/|\det A|\), and the volume element contributes \(|\det A|\); they cancel. This cancellation is why a Gaussian expectation can be evaluated as an expectation under a standard normal vector, with the function evaluated at \(\mu+Az\). It does not mean the Jacobian can be omitted when transforming an unweighted integral or a density itself.

Key result: If \(L(b)\) is the conditional likelihood of the observed responses given \(b\), the Gaussian marginal likelihood contribution is \(\int L(b)f_b(b)\,db=\int L(\mu+Az)\phi_d(z)\,dz\). The equality is exact when the stated Gaussian model holds and the integrals exist.

In one dimension, if \(b\sim N(\mu,\sigma^2)\) with \(\sigma>0\), use \(b=\mu+\sigma z\), so \(db=\sigma\,dz\). The density in the original coordinate is \(f_b(b)=\phi(z)/\sigma\), and thus

\int_{-\infty}^{\infty} h(b)f_b(b)\,db = \int_{-\infty}^{\infty} h(\mu+\sigma z)\phi(z)\,dz.

The positive standard deviation is the scale factor. Using the variance \(\sigma^2\) in place of \(\sigma\) would give the wrong transformation and Jacobian.

Worked Transformations

Worked Example: A Scalar Gaussian Expectation

Suppose a hypothetical random intercept satisfies \(b\sim N(0,1)\), and consider the integral \(I=E[\exp(-b^2/2)]\). This simple function isolates the coordinate transformation; the same change applies when the function is a conditional GLMM likelihood.

Here \(\sigma=1\), so \(b=z\) and \(db=dz\). The transformed integral is

I=\int_{-\infty}^{\infty}e^{-z^2/2}\phi(z)\,dz =\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}e^{-z^2}\,dz =\frac{1}{\sqrt{2}}.

For the last equality, use \(\int_{-\infty}^{\infty}e^{-z^2}\,dz=\sqrt{\pi}\); hence \(I=\sqrt{\pi}/\sqrt{2\pi}=1/\sqrt{2}\approx0.7071\). The transformed expression retains the standard-normal density as its integration weight.

Worked Example: A Bivariate Scale-and-Correlation Transformation

Let a hypothetical two-dimensional random effect have mean \(\mu=(0,0)^\mathsf{T}\), and choose \(A=\begin{pmatrix}1&0\\0.5&2\end{pmatrix}\). Its covariance is \(\Sigma=AA^\mathsf{T}=\begin{pmatrix}1&0.5\\0.5&4.25\end{pmatrix}\), and \(\det(A)=2\). Thus \(|\Sigma|=4\) and \(|\Sigma|^{1/2}=2\), as required.

Take the standard-normal coordinate \(z=(0.5,1)^\mathsf{T}\). The corresponding random effect is

b=Az= \begin{pmatrix}1&0\\0.5&2\end{pmatrix} \begin{pmatrix}0.5\\1\end{pmatrix} = \begin{pmatrix}0.5\\2.25\end{pmatrix}.

Since \(z^\mathsf{T}z=0.5^2+1^2=1.25\), its standard bivariate normal density is \(\phi_2(z)=(2\pi)^{-1}e^{-0.625}\approx0.08519\). The density at the corresponding \(b\) is therefore \(f_b(b)=\phi_2(z)/2\approx0.04259\). Multiplying this density by the Jacobian \(2\) recovers \(0.08519\). In an integral of \(h(b)f_b(b)\), the two factors cancel, leaving \(h(Az)\phi_2(z)\) in standard-normal coordinates.

Worked Example: Transforming a Random-Intercept Likelihood

Consider five hypothetical Bernoulli responses in one cluster, with three successes and two failures. Given a random intercept \(b\), let the success probability be \(p(b)=1/(1+e^{-b})\). The conditional likelihood for this particular response sequence is \(L(b)=p(b)^3[1-p(b)]^2\). Suppose \(b\sim N(0,0.64)\), so \(\sigma=0.8\).

The cluster's marginal likelihood contribution is \(\int L(b)f_b(b)\,db\). Set \(b=0.8z\), for which \(db=0.8\,dz\) and \(f_b(0.8z)=\phi(z)/0.8\). The transformed expression is

\int_{-\infty}^{\infty}L(b)f_b(b)\,db = \int_{-\infty}^{\infty} \left[\frac{1}{1+e^{-0.8z}}\right]^3 \left[1-\frac{1}{1+e^{-0.8z}}\right]^2 \phi(z)\,dz.

At the illustrative coordinate \(z=0.5\), \(b=0.8(0.5)=0.4\) and \(p(0.4)\approx0.5987\). Thus \(L(0.4)\approx0.5987^3(0.4013)^2\approx0.03456\). Also \(\phi(0.5)\approx0.3521\), so the transformed integrand at this coordinate is approximately \(0.03456(0.3521)=0.01217\). In the original coordinate, \(f_b(0.4)=0.3521/0.8\approx0.4401\); multiplying by \(L(0.4)\) gives approximately \(0.01521\) per unit of \(b\). Multiplying by \(db/dz=0.8\) gives \(0.01217\) per unit of \(z\), matching the transformed integrand. This pointwise check verifies the Jacobian; it is not, by itself, an evaluation of the full marginal likelihood.

Assumptions, Checks, and Failure Modes

The transformation requires a correctly specified Gaussian random-effects distribution and, for the full-dimensional formula above, a positive-definite covariance matrix. The conditional likelihood must be a function of the random effects at the transformed value \(\mu+Az\), and its Gaussian-weighted integral must exist. These conditions concern the probability model and the validity of the change of variables; they do not guarantee that any particular numerical integration method will be accurate.

Several direct checks catch common implementation errors:

Check the covariance: multiply \(AA^\mathsf{T}\) and confirm that it equals the model's \(\Sigma\). A factor that instead satisfies \(A^\mathsf{T}A=\Sigma\) is not automatically suitable for the stated mapping.
Check the scale: in one dimension, the multiplier is \(\sigma\), not \(\sigma^2\). Verify that the transformed random effect has variance \(\sigma^2\).
Check the absolute determinant: the volume multiplier is \(|\det A|\), while the transformed density includes its reciprocal.
Check the function argument: evaluate the conditional likelihood at \(\mu+Az\), not at \(z\) unless the model truly has that form.
Check dimensions: \(z\) and \(b\) are both \(d\)-vectors, \(A\) is \(d\) by \(d\), and the Jacobian is a scalar.

A common mistake is to transform the random effect in the likelihood but leave the original-coordinate Gaussian density unchanged. For example, replacing \(b\) by \(\sigma z\) while retaining \(f_b(\sigma z)\) and integrating with respect to \(dz\) introduces an erroneous factor \(1/\sigma\). The correction is to include \(db=\sigma\,dz\); then \(f_b(\sigma z)\sigma=\phi(z)\). Equivalently, use the standard-normal expectation formula from the outset.

How This Transformation Relates to Other Approaches

Direct integration in \(b\)-coordinates retains the original Gaussian density and integrates \(L(b)f_b(b)\). Standard-normal coordinates express exactly the same integral as an expectation of \(L(\mu+Az)\). The first representation can be natural when the original density is already convenient; the second makes the Gaussian weight fixed while the covariance enters through the argument transformation.

Laplace approximation is different: it approximates an integral using local curvature near a mode. A coordinate change alone does not use a quadratic approximation and does not discard parts of the integrand. Likewise, numerical quadrature is a means of approximating the resulting integral, not part of the change-of-variables identity itself. Mapping to standard-normal coordinates does not ensure an easy integral, nor does it make a numerical approximation exact.

If the covariance is singular, a nonsingular square \(A\) does not exist, so the full-dimensional density and Jacobian formula stated here cannot be applied as written. Such a distribution is supported on a lower-dimensional subspace and needs a parameterization of that support. Also, different valid factors \(A\) can represent the same covariance; the identity remains valid for each nonsingular factor satisfying \(\Sigma=AA^\mathsf{T}\), though the numerical behavior of a later integration procedure may depend on the chosen coordinates.

Key takeaway: For \(b\sim N_d(\mu,\Sigma)\), write \(b=\mu+Az\) with \(\Sigma=AA^\mathsf{T}\). The Jacobian is \(|\det A|\), which cancels the corresponding Gaussian density factor, giving \(\int L(b)f_b(b)\,db=\int L(\mu+Az)\phi_d(z)\,dz\).

Check Your Understanding

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

If \(b\sim N(2,9)\), what standard-normal transformation should be used, and what is \(db/dz\)?
For a two-dimensional transformation with \(\det(A)=-3\), what is the absolute Jacobian determinant?
Suppose \(A=\begin{pmatrix}2&0\\0&1\end{pmatrix}\). What covariance matrix does it produce under \(\Sigma=AA^\mathsf{T}\)?
In a density-weighted integral, why does the factor \(|\det A|\) disappear from the final standard-normal expectation?
A calculation substitutes \(b=\mu+Az\) into \(L(b)\), but integrates \(L(\mu+Az)f_b(\mu+Az)\) directly with respect to \(dz\). What is missing?

Solutions. (1) Use \(b=2+3z\), because the standard deviation is \(\sqrt{9}=3\); thus \(db/dz=3\). (2) The absolute Jacobian is \(|-3|=3\). (3) The covariance is \(\begin{pmatrix}4&0\\0&1\end{pmatrix}\). (4) The transformed density contributes \(1/|\det A|\), while the volume element contributes \(|\det A|\); their product is one. (5) The integration variable change requires \(db=|\det A|\,dz\), so the Jacobian is missing. Including it converts the original Gaussian density to \(\phi_d(z)\).