Shrinkage with Correlated Random Effects

Why One Parameter’s Data Can Inform Another

In a nonlinear mixed-effects model, an individual’s observations may be informative about one random effect but contain little direct information about another. If the population random effects are correlated, however, the data can still shift both empirical Bayes estimates (EBEs). The shift is not evidence that the observations directly measured both parameters: it arises because the model says that individuals with one random effect tend, on average, to have a related value of the other.

This distinction matters when interpreting individual estimates. A nonzero EBE for a weakly informed parameter may reflect a genuine cross-parameter update under the fitted covariance model, rather than a strong individual likelihood contribution for that parameter. The calculations below isolate this mechanism in Gaussian settings before discussing its use as an approximation in nonlinear models.

Definition: Let \(\boldsymbol{\eta}_i=(\eta_{i1},\eta_{i2})^\mathsf{T}\) be the two random effects for individual \(i\), with population distribution \(\boldsymbol{\eta}_i\sim N(\mathbf{0},\Omega)\), where \[ \Omega = \begin{pmatrix} \omega_{11} & \omega_{12}\\ \omega_{12} & \omega_{22} \end{pmatrix}. \] The off-diagonal entry \(\omega_{12}\) is the covariance between the random effects. The correlation is \(\rho=\omega_{12}/\sqrt{\omega_{11}\omega_{22}}\). An empirical Bayes estimate is an individual-specific estimate based on that individual’s data and the population model; in a Gaussian model, the conditional mean and conditional mode coincide.

The Gaussian Update and Its Cross-Parameter Effect

Assume that the random effects have a zero-mean multivariate normal distribution with positive-definite covariance matrix \(\Omega\). For a linear Gaussian observation model, write \(\mathbf{y}_i=H\boldsymbol{\eta}_i+\mathbf{e}_i\), where \(H\) describes which combinations of random effects enter the observations and \(\mathbf{e}_i\sim N(\mathbf{0},R)\). Assume \(R\) is positive definite and that the observation error is independent of the random effects. These assumptions make the conditional distribution of \(\boldsymbol{\eta}_i\) given \(\mathbf{y}_i\) multivariate normal.

The conditional covariance and mean are

C=(\Omega^{-1}+H^\mathsf{T}R^{-1}H)^{-1}, \qquad \widehat{\boldsymbol{\eta}}_i =C H^\mathsf{T}R^{-1}\mathbf{y}_i =\Omega H^\mathsf{T}(H\Omega H^\mathsf{T}+R)^{-1}\mathbf{y}_i.

The mean is the Gaussian EBE. The covariance describes uncertainty remaining in the random effects conditional on the observations and the assumed population parameters. To see why the formula holds, multiply the Gaussian likelihood by the Gaussian prior and collect the terms that depend on \(\boldsymbol{\eta}_i\).

\begin{aligned} \log p(\boldsymbol{\eta}_i\mid\mathbf{y}_i) &= \text{constant} -\frac{1}{2}\boldsymbol{\eta}_i^\mathsf{T}\Omega^{-1}\boldsymbol{\eta}_i -\frac{1}{2}(\mathbf{y}_i-H\boldsymbol{\eta}_i)^\mathsf{T} R^{-1}(\mathbf{y}_i-H\boldsymbol{\eta}_i)\\ &= \text{constant} -\frac{1}{2}\boldsymbol{\eta}_i^\mathsf{T} (\Omega^{-1}+H^\mathsf{T}R^{-1}H)\boldsymbol{\eta}_i +\boldsymbol{\eta}_i^\mathsf{T}H^\mathsf{T}R^{-1}\mathbf{y}_i. \end{aligned}

Completing the square gives a normal distribution with covariance \(C\) and mean \(C H^\mathsf{T}R^{-1}\mathbf{y}_i\). The equivalent expression using \(\Omega H^\mathsf{T}\) follows from the matrix identity for this Gaussian update. In particular, if an observation depends only on \(\eta_{i1}\), the update for \(\eta_{i2}\) depends on \(\operatorname{Cov}(\eta_{i2},y_i)\). When \(\omega_{12}\ne 0\), this covariance need not be zero.

Worked Example: An Observation on the First Random Effect

Suppose the hypothetical random effects are log-scale deviations, so both are dimensionless, with standard deviations \(0.4\) and \(0.6\), and correlation \(0.5\). An individual has one noisy summary measurement \(z=0.30\) of \(\eta_1\), with measurement error standard deviation \(0.2\). Use the model \(z=\eta_1+e\), where \(e\) is independent normal error. There is no direct observation of \(\eta_2\).

The prior covariance is

\Omega= \begin{pmatrix} 0.4^2 & 0.5(0.4)(0.6)\\ 0.5(0.4)(0.6) & 0.6^2 \end{pmatrix} = \begin{pmatrix} 0.16 & 0.12\\ 0.12 & 0.36 \end{pmatrix}.

Here \(H=(1\ \ 0)\) and \(R=(0.2^2)=(0.04)\). The measurement variance is \(\operatorname{Var}(z)=0.16+0.04=0.20\). The covariances of the two random effects with \(z\) are \(0.16\) and \(0.12\), respectively. Therefore,

\widehat{\boldsymbol{\eta}} = \frac{1}{0.20} \begin{pmatrix} 0.16\\ 0.12 \end{pmatrix}(0.30) = \begin{pmatrix} 0.24\\ 0.18 \end{pmatrix}.

The data shift the estimate of \(\eta_1\) to \(0.24\), and also shift the estimate of \(\eta_2\) to \(0.18\), despite no direct measurement of \(\eta_2\). The second estimate is positive because the fitted population correlation is positive and the observed first effect is positive. For comparison, the conditional covariance is \[ C=\Omega- \frac{1}{0.20} \begin{pmatrix}0.16\\0.12\end{pmatrix} \begin{pmatrix}0.16&0.12\end{pmatrix} = \begin{pmatrix}0.032&0.024\\0.024&0.288\end{pmatrix}. \] Thus uncertainty about \(\eta_2\) decreases from its prior variance \(0.36\) to \(0.288\), but remains substantial because it was not directly measured.

Direction and Size Depend on the Covariance

For the one-measurement setup, let \(z=\eta_1+e\), with \(\operatorname{Var}(e)=\tau^2\). The joint Gaussian formula reduces to

E(\eta_2\mid z) = \frac{\operatorname{Cov}(\eta_2,z)}{\operatorname{Var}(z)}z = \frac{\omega_{12}}{\omega_{11}+\tau^2}z.

This makes the mechanism explicit. The sign of the shift in \(\eta_2\) follows the product of the sign of \(z\) and the sign of \(\omega_{12}\). Its magnitude is reduced when measurement error is large relative to the variation in \(\eta_1\). If \(\omega_{12}=0\), the measurement does not shift the conditional mean of \(\eta_2\) in this setup. That does not mean every model with a zero prior covariance has uncoupled EBEs: the likelihood itself may contain information linking the two effects, for example through a nonzero off-diagonal entry in its local information matrix.

Worked Example: A Negative Correlation Reverses the Shift

Consider hypothetical log-scale random effects with standard deviations \(0.3\) and \(0.5\), and correlation \(-0.6\). A noisy measurement of \(\eta_1\) is \(z=0.20\), with error standard deviation \(0.10\). As before, the observation model is \(z=\eta_1+e\), with independent normal error, and there is no direct measurement of \(\eta_2\).

The covariance is \(\omega_{12}=(-0.6)(0.3)(0.5)=-0.09\), and the variance of \(z\) is \(0.3^2+0.1^2=0.10\). Consequently,

E(\eta_2\mid z) = \frac{-0.09}{0.10}(0.20) =-0.18.

Although the observed first effect is positive, the conditional estimate of the second effect is negative. This is the direction implied by the negative population correlation. It is a model-based conditional estimate, not proof that the individual’s unobserved parameter is negative.

When Both Parameters Have Data

If observations contain information about both random effects, their estimates can still influence one another. The relative size of that influence depends on both the population covariance and the individual likelihood information. The following calculation compares a joint update with the separate updates that would result from a diagonal covariance.

Worked Example: Comparing Joint and Separate Updates

Suppose the hypothetical random-effect covariance is the matrix from the first example. The individual has two noisy summaries \(\mathbf{y}=(0.30,-0.10)^\mathsf{T}\), one for each effect, with independent error variances \(0.04\) and \(0.09\). Thus \(H\) is the identity matrix and \(R=\operatorname{diag}(0.04,0.09)\). The assumptions are a bivariate normal prior and independent normal measurement errors with those stated variances.

The posterior precision and precision-weighted data are

\Omega^{-1}= \begin{pmatrix} 8.333&-2.778\\ -2.778&3.704 \end{pmatrix}, \qquad R^{-1}= \begin{pmatrix} 25&0\\ 0&11.111 \end{pmatrix}, \qquad R^{-1}\mathbf{y}= \begin{pmatrix} 7.500\\ -1.111 \end{pmatrix}.

Adding the precisions gives \[ \Omega^{-1}+R^{-1} = \begin{pmatrix} 33.333&-2.778\\ -2.778&14.815 \end{pmatrix}. \] The determinant is approximately \(33.333(14.815)-(-2.778)^2=486.111\). Inverting this matrix and multiplying by the precision-weighted data gives

\widehat{\boldsymbol{\eta}} = \frac{1}{486.111} \begin{pmatrix} 14.815&2.778\\ 2.778&33.333 \end{pmatrix} \begin{pmatrix} 7.500\\ -1.111 \end{pmatrix} \approx \begin{pmatrix} 0.222\\ -0.0333 \end{pmatrix}.

For comparison, if the covariance were diagonal with variances \(0.16\) and \(0.36\), each parameter would update independently. The separate estimates would be \(0.16/(0.16+0.04)(0.30)=0.240\) and \(0.36/(0.36+0.09)(-0.10)=-0.0800\). With the correlated covariance, the first estimate is slightly smaller and the second is closer to zero. Both observations contribute, but the joint prior changes how their information is combined.

Using the Result in a Nonlinear Mixed-Effects Model

In an NLME model, the individual likelihood is generally not Gaussian in the random effects. The exact conditional distribution is proportional to the individual likelihood times the multivariate normal random-effects density. Its mode (often called the conditional mode) and mean need not be equal, and the Gaussian update above is not generally an exact formula for either one.

A local quadratic approximation explains how the covariance can still couple the estimates. If the log likelihood for one individual is approximated near \(\boldsymbol{\eta}=\mathbf{0}\) by \[ \ell(\boldsymbol{\eta})\approx \ell(\mathbf{0})+\mathbf{g}^\mathsf{T}\boldsymbol{\eta} -\tfrac12\boldsymbol{\eta}^\mathsf{T}J\boldsymbol{\eta}, \] where \(\mathbf{g}\) is the likelihood score and \(J\) is a positive-semidefinite local information approximation, adding the normal log prior gives a quadratic posterior log density. Its mode solves

(\Omega^{-1}+J)\widehat{\boldsymbol{\eta}}=\mathbf{g}, \qquad \widehat{\boldsymbol{\eta}}=(\Omega^{-1}+J)^{-1}\mathbf{g}.

The inverse covariance, or prior precision, and the likelihood information jointly determine the update. Off-diagonal terms can transmit a score for one parameter into the estimate of another. This approximation is most credible when the individual likelihood is reasonably represented by a quadratic near the mode; marked asymmetry, multiple modes, or strong nonlinearity can make it inaccurate. In such cases, the actual conditional distribution or the fitting method’s computed EBEs must be interpreted on their own terms.

Diagnostics, Alternatives, and Limitations

A correlated-random-effects model is not automatically preferable to a diagonal one. The covariance should be supported by the population data and estimated with adequate precision. Useful checks include whether the estimated covariance matrix is positive definite, whether correlations are extreme or poorly determined, and whether estimates or individual predictions change materially under a plausible diagonal-covariance sensitivity fit. An estimated correlation near \(1\) or \(-1\) can indicate a nearly singular covariance and unstable separation of the two random effects. Sparse individual data can also leave an EBE highly dependent on the population covariance estimate.

A diagonal covariance model is a simpler alternative: it rules out prior cross-parameter association, so information about one parameter does not shift another through the prior. It may be more stable when the data cannot estimate correlations reliably, but it can misrepresent the population distribution if meaningful correlations exist. A full covariance model allows that association, at the cost of estimating more parameters. In nonlinear likelihoods, the individual data may also couple the effects directly; that likelihood coupling is distinct from the coupling introduced by \(\Omega\).

Common pitfalls include treating a cross-updated EBE as if the individual data directly identified that parameter, interpreting a population correlation as a causal relationship, and assuming the sign of the EBE shift is always the sign of the correlation. The simple one-measurement formula shows that the observation’s sign matters too; with multiple observations and nonlinear likelihoods, the total update depends on the full information structure. Correlation is an association in the random-effects distribution, not a causal mechanism.

Key takeaway: Correlated random effects allow an individual’s data to shift an EBE for a parameter with little or no direct information. The direction and size of that shift depend on the covariance, the observation or likelihood information, and the assumptions used to calculate the conditional estimate.

Check Your Understanding

For each question, give the answer before reading the solution in the same item.

If \(\omega_{12}=0\) and the likelihood measures only \(\eta_1\), what is the Gaussian conditional mean of \(\eta_2\)? Solution: It remains at its prior mean, zero. In the one-measurement formula, \(E(\eta_2\mid z)=\omega_{12}z/(\omega_{11}+\tau^2)=0\).
In the one-measurement setup, suppose \(z>0\) and \(\omega_{12}<0\). What is the direction of the shift in \(E(\eta_2\mid z)\)? Solution: It is negative because the denominator \(\omega_{11}+\tau^2\) is positive, so the sign is the sign of \(\omega_{12}z\).
Why does a larger measurement-error variance generally reduce the cross-parameter update in the one-measurement formula? Solution: The update is \(\omega_{12}z/(\omega_{11}+\tau^2)\). For fixed \(z\), \(\omega_{12}\), and \(\omega_{11}\), increasing \(\tau^2\) increases the denominator and reduces the magnitude.
Does a positive EBE for \(\eta_2\) after observing only \(\eta_1\) establish that the individual’s true \(\eta_2\) is positive? Solution: No. It is a conditional estimate under the fitted population covariance and observation model; uncertainty remains, and the estimate is not the unobserved individual value known without error.
Name one reason the Gaussian update may not equal an EBE computed from a nonlinear model. Solution: The individual likelihood may be asymmetric or otherwise non-quadratic in the random effects. Then the posterior need not be Gaussian, and its conditional mean and mode need not equal the local Gaussian update.