Gaussian Random-Effects Density and Covariance Parameterization

Why the Random-Effects Density Matters

In a GLMM, the conditional response model is combined with a distribution for the cluster-specific random effects. For cluster $i$, integrating the conditional likelihood over that distribution gives its marginal likelihood contribution. The random-effects density is therefore not just a convenient description of between-cluster variation: it is part of the integrand, and its scale, determinant, and coordinate system all matter in computation.

This tutorial focuses on a $q$-dimensional Gaussian random-effects vector. The response model supplies a conditional likelihood $L_i(b_i)=f(y_i\mid b_i)$, where $b_i$ is the random-effects vector for cluster $i$. The density $f(b_i)$ weights each possible vector in the integral $\int L_i(b_i)f(b_i)\,db_i$. The same probability model can be written in different coordinates, but the density and integration measure must be transformed together.

Definition: Suppose the $q$-vector of random effects for cluster $i$ satisfies $b_i\sim N_q(0,\Sigma)$, where $\Sigma$ is a symmetric positive-definite $q\times q$ covariance matrix. Its density with respect to $q$-dimensional volume $db_i$ is $$ f(b_i;\Sigma)=(2\pi)^{-q/2}|\Sigma|^{-1/2} \exp\left(-\frac{1}{2}b_i^\mathsf{T}\Sigma^{-1}b_i\right). $$

The quadratic form $b_i^\mathsf{T}\Sigma^{-1}b_i$ is nonnegative. Consequently, the exponent is nonpositive, and the density is largest at the origin and decreases as the random-effects vector moves away in the geometry determined by $\Sigma$. The determinant factor sets the density's overall scale: a distribution spread over a larger volume has a lower peak, all else equal. The density formula requires a nonsingular covariance matrix; positive definiteness guarantees that $|\Sigma|>0$ and that $\Sigma^{-1}$ exists.

Standard Deviations, Correlations, and Covariances

For two random effects, a useful parameterization separates their marginal scales from their association: $$ \Sigma= \begin{pmatrix} \tau_1^2 & \rho\tau_1\tau_2\\ \rho\tau_1\tau_2 & \tau_2^2 \end{pmatrix}, \qquad \tau_1>0,\quad \tau_2>0,\quad -1<\rho<1. $$ Here $\tau_1$ and $\tau_2$ are the standard deviations, and $\rho$ is the correlation. The off-diagonal covariance is $\rho\tau_1\tau_2$, not $\rho$ alone. In general, $\Sigma=D R D$, where $D$ is diagonal with the random-effects standard deviations and $R$ is their correlation matrix.

For the bivariate case, $|\Sigma|=\tau_1^2\tau_2^2(1-\rho^2)$. Substitution into the multivariate density gives $$ f(b_i;\tau_1,\tau_2,\rho)= \frac{1}{2\pi\tau_1\tau_2\sqrt{1-\rho^2}} \exp\left\{-\frac{1}{2(1-\rho^2)} \left[ \left(\frac{b_{i1}}{\tau_1}\right)^2 -2\rho\left(\frac{b_{i1}}{\tau_1}\right)\left(\frac{b_{i2}}{\tau_2}\right) +\left(\frac{b_{i2}}{\tau_2}\right)^2 \right]\right\}. $$

This form makes the scales and correlation explicit. It also shows why the correlation cannot be set to $1$ or $-1$ in a nonsingular bivariate Gaussian density: then $1-\rho^2=0$, so the covariance matrix is singular and the displayed density is undefined.

Worked Example: Evaluating a Bivariate Random-Effects Density

Suppose a hypothetical random intercept and random slope have standard deviations $\tau_1=2$ and $\tau_2=1$, with correlation $\rho=0.6$. Evaluate their joint density at $b_i=(1,-0.5)^\mathsf{T}$. The covariance matrix is $$ \Sigma= \begin{pmatrix} 4 & 1.2\\ 1.2 & 1 \end{pmatrix}. $$

Its determinant is $4(1)-(1.2)^2=2.56$, which is positive. Its leading diagonal entry is positive as well, so this symmetric $2\times2$ matrix is positive definite. The inverse is $$ \Sigma^{-1}=\frac{1}{2.56} \begin{pmatrix} 1 & -1.2\\ -1.2 & 4 \end{pmatrix}. $$

The quadratic form is $$ b_i^\mathsf{T}\Sigma^{-1}b_i =\frac{1}{2.56}\left[1^2-2(1.2)(1)(-0.5)+4(-0.5)^2\right] =\frac{3.2}{2.56}=1.25. $$

Thus the density at this point is $$ f(b_i;\Sigma)=\frac{1}{2\pi\sqrt{2.56}}\exp(-1.25/2) =\frac{1}{3.2\pi}\exp(-0.625) \approx 0.05324. $$

The value is a density, not the probability that the random-effects vector equals exactly $(1,-0.5)^\mathsf{T}$; a continuous random vector assigns probability zero to any single point. The density is measured per unit of area in the intercept-slope coordinate system.

Changing Coordinates to Standard Normal Variables

A common integration coordinate is a standard-normal vector $z_i\sim N_q(0,I_q)$. Let $A$ be a nonsingular matrix satisfying $AA^\mathsf{T}=\Sigma$, and set $b_i=Az_i$. Then $b_i$ has mean zero and covariance $A I_q A^\mathsf{T}=\Sigma$. A Cholesky factor is one standard choice for $A$, though the density representation does not depend on choosing that particular factor.

Change-of-variables result: If $b_i=Az_i$ for nonsingular $A$, then $db_i=|A|\,dz_i$, and $$ f(b_i;\Sigma)\,db_i =(2\pi)^{-q/2}\exp\left(-\frac{1}{2}z_i^\mathsf{T}z_i\right)\,dz_i. $$ Thus the random-effects integral can be written either in $b_i$-coordinates or in standard-normal $z_i$-coordinates, provided the transformed likelihood is evaluated at $b_i=Az_i$.

To verify the result, use $\Sigma=AA^\mathsf{T}$, so $|\Sigma|=|A|^2$ and $\Sigma^{-1}=A^{-\mathsf{T}}A^{-1}$. With $b_i=Az_i$, $$ b_i^\mathsf{T}\Sigma^{-1}b_i =z_i^\mathsf{T}A^\mathsf{T}A^{-\mathsf{T}}A^{-1}Az_i =z_i^\mathsf{T}z_i. $$

The original density's normalizing factor is $|\Sigma|^{-1/2}=|A|^{-1}$, while the volume element contributes $|A|$. They cancel, leaving the standard-normal density. This cancellation is essential: evaluating $f(Az_i;\Sigma)$ but omitting the Jacobian would not preserve the integral.

Worked Example: Transforming a Random-Effects Integral

Consider a hypothetical one-dimensional random effect $b\sim N(0,4)$, and a conditional response density $f(y\mid b)$. Write $b=2z$, where $z\sim N(0,1)$. The marginal density of an observed response value $y$ is $$ \int_{-\infty}^{\infty} f(y\mid b)\frac{1}{2\sqrt{2\pi}}\exp\left(-\frac{b^2}{8}\right)\,db. $$

Under $b=2z$, the derivative is $db/dz=2$, so $db=2\,dz$. Substituting both the transformed random effect and the volume element gives $$ \int_{-\infty}^{\infty} f(y\mid 2z) \frac{1}{2\sqrt{2\pi}}\exp\left(-\frac{(2z)^2}{8}\right)(2)\,dz = \int_{-\infty}^{\infty} f(y\mid 2z) \frac{1}{\sqrt{2\pi}}\exp\left(-\frac{z^2}{2}\right)\,dz. $$

The two integrals are equal because the factor $2$ from the Jacobian cancels the $1/2$ in the original density's normalizing constant. The likelihood must still be evaluated at $b=2z$. Integrating $f(y\mid z)$ against a standard-normal density instead would describe a different random-effects scale.

Parameterizing the Covariance Matrix

The density is mathematically determined by $\Sigma$, but a statistical computation needs parameters that keep $\Sigma$ in the valid covariance-matrix space. Directly estimating every covariance entry is possible if positive definiteness is enforced, but unconstrained entries can produce an invalid matrix. In the bivariate case, a parameterization by $\tau_1>0$, $\tau_2>0$, and $-1<\rho<1$ builds validity into the allowed parameter values.

For more dimensions, one may parameterize a Cholesky factor $A$ with positive diagonal entries and set $\Sigma=AA^\mathsf{T}$. This guarantees positive definiteness. Another option is to work with standard deviations and a valid correlation matrix, but its correlations must jointly define a positive-definite matrix; checking each pair separately is not enough when $q>2$. Computational implementations often use unconstrained parameters internally, transforming them to positive standard deviations and valid correlations. Such transformations are choices of coordinates for optimization or integration, not changes to the intended probability model.

Worked Example: Checking and Repairing an Invalid Covariance

Suppose a hypothetical bivariate covariance estimate has variances $1$ and $1$, and covariance $1.1$: $$ \Sigma_{\text{candidate}}= \begin{pmatrix} 1 & 1.1\\ 1.1 & 1 \end{pmatrix}. $$

The implied correlation is $1.1/(\sqrt{1}\sqrt{1})=1.1$, outside the required interval $(-1,1)$. Equivalently, the determinant is $1-(1.1)^2=-0.21$, so the matrix is not positive definite. It cannot be the covariance matrix of a nonsingular bivariate Gaussian random-effects distribution, and substituting it into the Gaussian density would not give a valid density.

A valid alternative with the same marginal variances and correlation $0.8$ is $$ \Sigma_{\text{valid}}= \begin{pmatrix} 1 & 0.8\\ 0.8 & 1 \end{pmatrix}, \qquad |\Sigma_{\text{valid}}|=1-0.8^2=0.36>0. $$

This repaired matrix is positive definite, but it is not a correction justified by the data on its own: it is an example of a valid parameter value. In estimation, the fitting method must enforce positive definiteness rather than silently treating an invalid candidate as a permissible covariance matrix.

How Parameterization Affects the Integration Density

For cluster $i$, the marginal likelihood contribution is $$ L_i^{\mathrm{marg}}=\int L_i(b_i)f(b_i;\Sigma)\,db_i. $$

In $z_i$-coordinates it is instead $$ L_i^{\mathrm{marg}}=\int L_i(Az_i)(2\pi)^{-q/2} \exp\left(-\frac{1}{2}z_i^\mathsf{T}z_i\right)\,dz_i. $$

These are the same integral. However, the integrand as a function of its coordinate argument differs: one is expressed in the original random-effects units, and the other in standard-normal units. A numerical integration rule therefore encounters different coordinate scales and shapes. The Jacobian is part of the transformation, not an optional adjustment. Likewise, the determinant term in the original Gaussian density cannot be discarded when evaluating the density in $b_i$-coordinates.

This distinction matters in likelihood computation because a density value depends on the coordinates and units used to define volume. Probability assigned to a region does not change when the region and density are transformed consistently. For example, a peak density can change when a variable is rescaled even though the corresponding probabilities are identical. Parameterization can also affect numerical stability: correlations near $\pm1$ make $1-\rho^2$ small, while very different standard deviations can create poorly scaled coordinates.

Assumptions, Diagnostics, and Common Pitfalls

The formulas above assume a Gaussian random-effects distribution with mean zero and a positive-definite covariance matrix. In the usual independent-cluster likelihood, the cluster random-effects vectors are modeled as independent draws from the specified distribution. The conditional response model is then combined with the Gaussian density through the integral; the Gaussian assumption concerns the random effects, not necessarily the response distribution. Conditional independence assumptions for the responses are those of the specified GLMM and are not implied by the covariance parameterization itself.

Before using a Gaussian density or a numerical integral, check the following:

Symmetry and positive definiteness. Confirm $\Sigma$ is symmetric and has positive eigenvalues, or use a parameterization that guarantees this property. A zero determinant means the ordinary full-dimensional density formula does not apply.
Scale and correlation consistency. Each diagonal element is a variance, so its square root is a standard deviation. Each correlation is covariance divided by the product of standard deviations and must lie between $-1$ and $1$.
Normalization and Jacobians. In $b_i$-coordinates, retain $|\Sigma|^{-1/2}$. After changing coordinates, include the determinant of the transformation in the volume element; verify the transformed density integrates to one.
Numerical conditioning. Very small eigenvalues, extreme scale differences, or correlations close to their boundaries can make inverse and determinant calculations unstable. These are warnings to inspect the parameterization and computation, not proof by themselves that the Gaussian model is wrong.

A common mistake is to write $b_i=Az_i$, replace the covariance density by a standard-normal density, and then integrate over $db_i$ without a Jacobian—or to integrate over $dz_i$ while leaving the likelihood as $L_i(b_i)$ rather than $L_i(Az_i)$. The correction is to transform both the integrand and its volume element: $L_i(b_i)f(b_i)\,db_i=L_i(Az_i)\phi_q(z_i)\,dz_i$, where $\phi_q$ is the standard $q$-variate normal density.

Direct covariance-entry parameterization is transparent but requires constraints; standard-deviation and correlation parameterization gives interpretable components but must preserve a valid correlation matrix; Cholesky parameterization enforces positive definiteness and is often computationally convenient. These are alternatives for representing the same covariance structure when translated correctly. None removes the need to specify the Gaussian random-effects assumption, and no coordinate change can compensate for a covariance model that is substantively inadequate.

Key takeaway: A Gaussian random-effects density is defined jointly by its covariance matrix, its determinant normalization, and its coordinate volume. Standard-deviation, correlation, covariance, and standard-normal representations agree only when their transformations—including the Jacobian—are handled consistently.

Check Your Understanding

Answer each question, then compare with the solutions below.

For a bivariate Gaussian random effect with standard deviations $3$ and $2$, and correlation $-0.5$, what is the covariance between the two components?
Why must the determinant of a nonsingular covariance matrix be positive in the Gaussian density formula?
If $b=3z$, where $z\sim N(0,1)$, what Jacobian factor converts integration with respect to $db$ to integration with respect to $dz$?
A fitted bivariate covariance matrix has variances $4$ and $9$, and covariance $7$. Is it positive definite? Show the check.
Does a change from $b$-coordinates to standard-normal $z$-coordinates change the marginal likelihood? State what must change along with the coordinates.

Solutions. (1) The covariance is $\rho\tau_1\tau_2=(-0.5)(3)(2)=-3$. (2) Positive definiteness gives $|\Sigma|>0$, which makes $|\Sigma|^{-1/2}$ real and finite and ensures the inverse in the quadratic form exists. (3) Since $db/dz=3$, $db=3\,dz$, so the Jacobian factor is $3$. (4) The determinant is $4(9)-7^2=36-49=-13<0$; the matrix is not positive definite. (5) No. The likelihood argument must become $L(Az)$, and the density and integration measure must be transformed together, including the Jacobian; then the integral is unchanged.