Random-Effect Correlation in Approximate Likelihoods

Why Correlation Changes an Approximate Likelihood

In an NLME model, random effects can represent related subject-specific features, such as clearance and volume. If those effects are correlated, a subject with a higher-than-typical value of one may tend to have a higher or lower value of the other. This relationship matters even when the model’s predictions are locally linear: it changes which combinations of effects are likely, how uncertainty propagates to observations, and the determinant terms that normalize a likelihood.

The key is to carry the full random-effect covariance matrix through each calculation. Replacing it with its diagonal entries is not a harmless simplification: it changes the prior penalty, the induced observation covariance, and often the curvature determinant. The formulas below use a subject with \(q\) random effects and \(m\) observations, additive Gaussian residual error with covariance \(R\), and a positive-definite random-effect covariance matrix \(\Omega\). The focus is the role of correlation, not the choice of residual-error model.

Covariance, Inverse, and Determinant

For two random effects with standard deviations \(\omega_1,\omega_2\) and correlation \(\rho\), the covariance matrix is

A valid nonsingular covariance matrix requires \(\omega_1>0\), \(\omega_2>0\), and \(|\rho|<1\). Its inverse is

The off-diagonal terms in the inverse mean that the prior penalty is not generally a sum of separate penalties for each random effect. For a vector \(\eta=(\eta_1,\eta_2)^\mathsf{T}\), that penalty is

Thus, the relative plausibility of a pair of effects depends on their joint direction, not only on their separate magnitudes. The determinant also changes with \(\rho\): as the absolute correlation approaches one, the bivariate distribution becomes concentrated near a line, and \(|\Omega|\) approaches zero.

FO: Correlation Propagates Into Observation Covariance

Let \(y\) be the vector of observations and \(f(\eta)\) the model prediction vector at the fixed effects under consideration. In the first-order (FO) approximation, linearize the prediction around the population random-effect value \(\eta=0\): $$ f(\eta)\approx f(0)+G\eta, $$ where \(G\) is the \(m\times q\) derivative matrix at zero. With additive Gaussian residual error of covariance \(R\), the linearized observation model is Gaussian. Integrating over the random effects gives mean \(f(0)\) and covariance \(S=R+G\Omega G^\mathsf{T}\). In particular, the off-diagonal elements of \(\Omega\) contribute to the covariance among observations through \(G\Omega G^\mathsf{T}\).

This result follows because \(G\eta\) has covariance \(G\Omega G^\mathsf{T}\), and independent residual error contributes covariance \(R\). Independent Gaussian terms add their covariance matrices. The resulting FO subject log likelihood is

Both the determinant and inverse of \(S\) matter: the determinant normalizes the density, while the inverse weights the residual vector according to its covariance. Equivalently, the same Gaussian integral can be expressed using the random-effect precision. For fixed \(G\) and \(R\), define $$ H=\Omega^{-1}+G^\mathsf{T}R^{-1}G. $$ The determinant identity $$ |R+G\Omega G^\mathsf{T}|=|R|\,|\Omega|\,|H| $$ shows explicitly how the prior determinant and inverse enter a precision-based calculation. This identity follows from \(|I+AB|=|I+BA|\), applied to the matrices in the covariance sum.

FOCE and Laplace: The Prior Remains Correlated

FOCE uses a subject-specific conditional mode, denoted \(\hat\eta\), and linearizes the prediction at that point. Write \(J\) for the derivative matrix at \(\hat\eta\). The local prediction is \(f(\hat\eta)+J(\eta-\hat\eta)\). To see how the correlated prior enters, combine this linearized likelihood with the same prior \(N(0,\Omega)\). The prior is still centered at zero; centering the linearization at the mode does not replace the prior by a distribution centered there.

Under the stated additive Gaussian residual model, integrating this linearized likelihood against the prior gives the local Gaussian expression with observation covariance \(R+J\Omega J^\mathsf{T}\), and linearized mean \(f(\hat\eta)-J\hat\eta\). The prior precision contributes \(\Omega^{-1}\) to the quadratic form. Implementations and FOCE variants can differ in how they form the approximate subject objective, so the precise convention must be identified; in all cases, a correlated random-effect prior requires its full covariance structure. A conditional linearization should not be confused with Laplace’s approximation, which uses local curvature of the log integrand at its mode.

For Laplace, let \(h(\eta)=\log p(y\mid\eta)+\log p(\eta)\), and suppose \(h\) has an interior mode \(\hat\eta\) with negative-definite Hessian. Define the positive-definite curvature matrix \(H=-\nabla^2h(\hat\eta)\). A second-order Taylor expansion gives \(h(\eta)\approx h(\hat\eta)-\frac12(\eta-\hat\eta)^\mathsf{T}H(\eta-\hat\eta)\). Integrating this local Gaussian yields

Because the log prior contains \(-\frac12\eta^\mathsf{T}\Omega^{-1}\eta-\frac12\log|\Omega|\), the prior contributes \(\Omega^{-1}\) to \(H\) when the likelihood curvature is added, as well as the separate normalizing term \(-\frac12\log|\Omega|\) inside \(h(\hat\eta)\). Under a locally linear Gaussian likelihood with derivative \(J\), \(H=\Omega^{-1}+J^\mathsf{T}R^{-1}J\). With nonlinear predictions, the likelihood curvature can also include second-derivative terms; it is not generally just \(J^\mathsf{T}R^{-1}J\).

Checks, Alternatives, and Limitations

A practical check is to verify that \(\Omega\) is symmetric positive definite before evaluating densities or approximations. For the bivariate case, check positive variances and \(|\rho|<1\). In higher dimensions, positive diagonal entries alone are not sufficient; all eigenvalues must be positive. Cholesky factorization is a standard numerical way to check and work with positive-definite covariance matrices without explicitly forming an inverse.

When correlations are near \(1\) or \(-1\), \(|\Omega|\) can be very small and the inverse can have large entries. This can cause unstable determinant calculations, difficult optimization, or a nearly singular curvature matrix. A failed factorization, extreme parameter estimates, or sensitivity to small changes in starting values are warnings to investigate—not evidence that correlation should simply be discarded. Compute log determinants using a Cholesky factor, inspect the estimated correlations and eigenvalues, and verify that covariance and precision calculations use the same parameterization.

The main alternatives differ in where the local approximation is made. FO linearizes at \(\eta=0\); FOCE linearizes at a subject-specific conditional mode; Laplace uses a second-order approximation to the log integrand at that mode. Each must account for the same correlated prior, but their treatment of the nonlinear likelihood differs. None removes the need to check model fit or guarantees that an approximation is accurate when the integrand is strongly non-Gaussian, has multiple important modes, or is poorly identified. Those are limitations of the local approximation, not problems solved merely by adding a correlation parameter.