Why Change Coordinates?
In a nonlinear mixed-effects model, a subject’s likelihood conditional on that subject’s random effects is integrated over the random-effects distribution to obtain a marginal likelihood. When the random effects are correlated Gaussian variables, the integral is naturally written with a multivariate normal density. Re-expressing it in coordinates with independent standard-normal components can simplify the integral and make its geometry easier to work with.
The coordinate change does not approximate the likelihood or alter the random-effects distribution. It changes only the variables of integration. For that equality to hold, the transformation and its Jacobian must be handled together: omitting the Jacobian, or using the wrong matrix orientation, changes the integral.
The zero-mean convention is standard for random effects in this setting. More generally, if \(\eta_i\sim N(\mu,\Omega)\), use \(\eta_i=\mu+Lz_i\); the translation by \(\mu\) does not change the Jacobian. The covariance factorization requires \(\Omega\) to be positive definite for an ordinary, nonsingular \(q\)-dimensional transformation.
Deriving the Change of Variables
Write \(L\) for the lower-triangular Cholesky factor. Since \(\eta_i=Lz_i\), the derivative matrix of \(\eta_i\) with respect to \(z_i\) is \(L\). The multivariate change-of-variables rule therefore gives \(d\eta_i=|\det L|\,dz_i\). Because \(L\) has positive diagonal entries, \(\det L>0\), so \(|\det L|=\det L\). Also, \(\det(\Omega)=\det(L)\det(L^\mathsf{T})=(\det L)^2\), hence \(\det L=\sqrt{\det(\Omega)}\).
The density of \(\eta_i\) is the multivariate normal density. Substituting \(\eta_i=Lz_i\) into its quadratic form gives \(\eta_i^\mathsf{T}\Omega^{-1}\eta_i=z_i^\mathsf{T}z_i\), because \(\Omega^{-1}=L^{-\mathsf{T}}L^{-1}\). Together, the density and volume element become
where \(\phi_q(z_i)=(2\pi)^{-q/2}\exp(-z_i^\mathsf{T}z_i/2)\) is the standard \(q\)-dimensional normal density. The factor \(\sqrt{\det(\Omega)}\) from the Jacobian cancels the same factor in the denominator of the random-effects density. This cancellation is why the transformed integral has standard-normal measure, not an extra covariance-dependent multiplier.
For subject \(i\), let \(f_i(y_i\mid\eta_i,\theta)\) denote the conditional data density, with \(y_i\) the observed data and \(\theta\) the model’s fixed-effect and residual parameters. The marginal likelihood contribution is
The conditional model itself has not changed: wherever it previously used \(\eta_i\), it now uses \(Lz_i\). The equality follows from the change-of-variables rule, provided the original integral exists and the transformation is nonsingular.
Worked Example: Transforming Two Correlated Random Effects
Suppose the hypothetical random-effects covariance matrix is \(\Omega=\begin{pmatrix}4&1.2\\1.2&1\end{pmatrix}\). Its Cholesky factor is \(L=\begin{pmatrix}2&0\\0.6&0.8\end{pmatrix}\), since \(LL^\mathsf{T}=\begin{pmatrix}4&1.2\\1.2&0.36+0.64\end{pmatrix}=\Omega\).
For standard-normal coordinates \(z_i=(1,1)^\mathsf{T}\), the corresponding random effects are \[ \eta_i=Lz_i= \begin{pmatrix}2&0\\0.6&0.8\end{pmatrix} \begin{pmatrix}1\\1\end{pmatrix} = \begin{pmatrix}2\\1.4\end{pmatrix}. \] Thus a standard-normal draw with two independent unit-variance components maps to a pair of correlated effects with covariance \(\Omega\).
The Jacobian determinant is \(\det(L)=2(0.8)=1.6\), which also equals \(\sqrt{\det(\Omega)}\), because \(\det(\Omega)=4(1)-1.2^2=2.56\). At \(\eta_i=(2,1.4)^\mathsf{T}\), the Gaussian quadratic form is \(z_i^\mathsf{T}z_i=2\). The correlated random-effects density is therefore \[ p_{\eta}(\eta_i)=\frac{1}{2\pi(1.6)}\exp(-1) \approx 0.03659. \] Multiplying by the Jacobian gives \(p_{\eta}(\eta_i)\det(L)\approx 0.03659(1.6)\approx 0.05854\) using the displayed rounded density. Using the unrounded density, the product is approximately \(0.05855\). This agrees with the standard bivariate normal density at \(z_i=(1,1)^\mathsf{T}\): \[ \phi_2(z_i)=\frac{1}{2\pi}\exp(-1)\approx 0.05855. \] The slight difference between the two displayed products is rounding, not a difference in the transformed measure.
What the Transformed Integral Means
The integral in standard-normal coordinates is an expectation with respect to \(z_i\sim N(0,I_q)\). In particular, \(L_i(\theta,\Omega)=E\{f_i(y_i\mid Lz_i,\theta)\}\), where the expectation is over standard-normal \(z_i\). This notation highlights the two parts of the integrand: the standard-normal measure and the conditional likelihood evaluated at the mapped random effects. It does not itself specify how the expectation will be numerically evaluated.
The transformation can also be checked in the scalar case. If \(\eta_i\sim N(0,\omega^2)\) with \(\omega>0\), then \(L=\omega\), \(z_i=\eta_i/\omega\), and \(d\eta_i=\omega\,dz_i\). The density contributes \(1/\omega\); the volume element contributes \(\omega\). Their product is the standard-normal density times \(dz_i\), exactly as in the multivariate derivation.
Worked Example: A Nonlinear Integrand in Standard-Normal Coordinates
Suppose a hypothetical scalar random effect has distribution \(\eta_i\sim N(0,0.25)\), so its standard deviation is \(0.5\), and consider the nonlinear function \(g(\eta_i)=\exp(\eta_i)\). The Cholesky factor is the scalar \(L=0.5\), so \(\eta_i=0.5z_i\) for \(z_i\sim N(0,1)\). The transformed expectation is \[ E\{g(\eta_i)\} =\int_{-\infty}^{\infty}\exp(0.5z_i)\phi_1(z_i)\,dz_i. \]
Using \(\phi_1(z_i)=(2\pi)^{-1/2}\exp(-z_i^2/2)\), complete the square: \[ -\frac{z_i^2}{2}+0.5z_i =-\frac{(z_i-0.5)^2}{2}+\frac{0.125}{1}. \] Thus the integrand is \(\exp(0.125)\) times a normal density centered at \(0.5\) with variance \(1\). That density integrates to one, giving \[ E\{\exp(\eta_i)\}=\exp(0.125)\approx 1.133. \] This calculation shows that a nonlinear function is simply evaluated at the mapped value \(0.5z_i\); the transformation does not replace it with a linear approximation.
Determinant Check Through a Gaussian Integral
A second check comes from integrating a Gaussian-shaped function against the random-effects distribution. For positive-definite \(\Omega\), \(E\{\exp(-\eta_i^\mathsf{T}\eta_i/2)\}=\det(I_q+\Omega)^{-1/2}\). To justify the result, transform to standard-normal coordinates and combine the exponents: \[ -\frac{1}{2}z_i^\mathsf{T}z_i-\frac{1}{2}z_i^\mathsf{T}L^\mathsf{T}Lz_i =-\frac{1}{2}z_i^\mathsf{T}(I_q+L^\mathsf{T}L)z_i. \] The Gaussian integral of this quadratic form is \(\det(I_q+L^\mathsf{T}L)^{-1/2}\). The identity \(\det(I_q+L^\mathsf{T}L)=\det(I_q+LL^\mathsf{T})=\det(I_q+\Omega)\) yields the stated expression. The integral is finite because \(I_q+\Omega\) is positive definite.
Worked Example: Checking the Covariance Dependence
Use the hypothetical covariance matrix from the earlier example and let \(g(\eta_i)=\exp(-\eta_i^\mathsf{T}\eta_i/2)\). The formula gives \[ E\{g(\eta_i)\}=\frac{1}{\sqrt{\det(I_2+\Omega)}}. \] Here \(I_2+\Omega=\begin{pmatrix}5&1.2\\1.2&2\end{pmatrix}\), so \[ \det(I_2+\Omega)=5(2)-1.2^2=10-1.44=8.56, \qquad E\{g(\eta_i)\}=\frac{1}{\sqrt{8.56}}\approx 0.3418. \]
In standard-normal coordinates the same quantity is \(\int \exp(-z_i^\mathsf{T}L^\mathsf{T}Lz_i/2)\phi_2(z_i)\,dz_i\). Combining its exponents gives the determinant \(\det(I_2+L^\mathsf{T}L)\), which equals \(\det(I_2+\Omega)=8.56\). The result is unchanged by the coordinate transformation; the covariance dependence remains in the mapped integrand.
Practical Checks and Common Failure Modes
The transformed integral is useful only if the transformation is implemented consistently. The following checks catch errors before they are mistaken for features of a model or its likelihood:
- Check the factorization. Verify numerically that \(LL^\mathsf{T}\) reproduces \(\Omega\), including its off-diagonal entries.
- Check the direction of the map. With the stated lower-triangular convention, map standard coordinates by \(\eta_i=Lz_i\), and recover them by \(z_i=L^{-1}\eta_i\). Using \(L^\mathsf{T}z_i\) generally gives the wrong covariance.
- Check the measure. If transforming the density directly, include \(\det(L)\) in \(p_{\eta}(\eta_i)d\eta_i\). If the integral has already been rewritten with \(\phi_q(z_i)dz_i\), do not multiply by the Jacobian a second time.
- Check covariance empirically. Standard-normal draws should have sample means near zero and covariance near \(I_q\). Their mapped values \(Lz_i\) should have covariance near \(\Omega\), allowing for ordinary sampling variation.
- Watch for near-singularity. A positive-definite matrix can still be poorly conditioned. A nearly singular \(\Omega\) can make numerical factorization and the inverse map sensitive to small rounding errors.
If \(\Omega\) is singular, a full-dimensional Cholesky factor with positive diagonal does not exist, and the density with respect to ordinary \(q\)-dimensional volume is not defined in the usual way. A lower-dimensional representation on the supported subspace may be appropriate, but it is a different integration setup and must be stated explicitly.
Comparing Coordinate Choices
Direct integration in \(\eta_i\)-coordinates and integration in \(z_i\)-coordinates are mathematically equivalent when the Jacobian is included. Direct coordinates may be convenient when the model or a numerical routine is already expressed in the original random effects. Standard-normal coordinates separate the fixed covariance transformation from the common standard-normal measure, which can make integration routines easier to formulate consistently.
The Cholesky factor is a practical choice because it is triangular and has a positive diagonal for a positive-definite covariance matrix. Another valid square root \(A\) satisfying \(AA^\mathsf{T}=\Omega\) also gives \(\eta_i=Az_i\) with standard-normal \(z_i\); the determinant magnitude and distributional result are the same. Different valid roots can rotate or otherwise reorient the standard-normal coordinates, so their integrands need not look identical point by point even though their integrals agree. The Cholesky convention fixes one reproducible map.
Check Your Understanding
Answer each question, then compare with the solutions below.
- For \(L=\begin{pmatrix}2&0\\1&3\end{pmatrix}\), calculate \(\Omega=LL^\mathsf{T}\) and the Jacobian determinant.
- If \(\eta_i=Lz_i\), which covariance matrix does \(\eta_i\) have when \(z_i\sim N(0,I_q)\)?
- Why is there no extra factor of \(\det(L)\) in the final standard-normal integral once it is written with \(\phi_q(z_i)dz_i\)?
- What problem arises if one uses \(L^\mathsf{T}z_i\) instead of \(Lz_i\) under the lower-triangular Cholesky convention?
Solutions. (1) \(LL^\mathsf{T}=\begin{pmatrix}4&2\\2&10\end{pmatrix}\), and \(\det(L)=2(3)=6\). (2) Its covariance is \(L I_q L^\mathsf{T}=LL^\mathsf{T}=\Omega\). (3) The Jacobian factor has already canceled the \(\sqrt{\det(\Omega)}\) denominator in the Gaussian density; including it again would alter the integral. (4) In general, \(\operatorname{Cov}(L^\mathsf{T}z_i)=L^\mathsf{T}L\), which need not equal \(\Omega=LL^\mathsf{T}\), so it would produce the wrong random-effects covariance.