1 Random effects in hierarchical models

1.1 Definition and motivation

Random effects are latent components in hierarchical (multilevel) statistical models that account for systematic differences among groups, subjects, or clusters. Instead of treating group-specific deviations as fixed parameters to be estimated independently, random effects assume these deviations arise from a common probability distribution. This approach supports partial pooling (also called shrinkage), improves generalization to new groups, and reduces the number of free parameters.

In practice, random effects describe heterogeneity that is expected to vary by grouping structure—such as participants in a study, classrooms in an education dataset, or clinics in a health registry. The “covariance of random effects” refers to how different latent deviations co-move, for example whether groups with higher baseline outcomes also tend to respond more strongly to a predictor.

1.2 Common model formulations

A standard formulation expresses the outcome as a combination of fixed effects, random effects, and residual error. For a linear mixed-effects model, one common form is \[ y_{ij} = x_{ij}^\top \beta + z_{ij}^\top b_j + \varepsilon_{ij}, \] where \(i\) indexes observations within group \(j\), \(x_{ij}\) are fixed-effect covariates with coefficients \(\beta\), and \(z_{ij}\) are design terms that load random effects \(b_j\) for each group. Typically, \[ b_j \sim \mathcal{N}(0, G), \quad \varepsilon_{ij} \sim \mathcal{N}(0, R), \] with \(G\) capturing covariance among random effects and \(R\) capturing residual covariance (often diagonal, but not always).

Generalized mixed models extend the idea to non-Gaussian outcomes through link functions and latent random effects, with the covariance of random effects still playing a central role in how the model distributes variability across groups.

1.3 Random intercepts and random slopes

Random intercepts allow each group to have its own baseline level. If a model includes a random intercept, group \(j\) has a deviation \(b_{0j}\) from the overall intercept.

Random slopes allow the effect of a predictor to vary by group. If the slope for a covariate \(x\) is random, group \(j\) has a deviation \(b_{1j}\) from the global slope. A common two-part random-effects vector is \[ b_j = \begin{pmatrix} b_{0j} \\ b_{1j} \end{pmatrix}, \] with covariance matrix \[ G = \begin{pmatrix} \operatorname{Var}(b_{0j}) & \operatorname{Cov}(b_{0j}, b_{1j})\\ \operatorname{Cov}(b_{0j}, b_{1j}) & \operatorname{Var}(b_{1j}) \end{pmatrix}. \] The covariance term \(\operatorname{Cov}(b_{0j}, b_{1j})\) determines whether groups with a higher-than-average baseline also tend to have a larger or smaller slope.

1.4 Distributional assumptions (e.g., normality)

Most classical covariance formulations assume random effects follow a multivariate normal distribution. Under that assumption, the covariance matrix \(G\) fully characterizes second-order properties such as variances and covariances of random effects. Normality also simplifies likelihood construction and enables closed-form expressions for certain marginal and conditional moments in Gaussian mixed models.

In non-Gaussian settings, normality is often treated as a working assumption for the latent random effects distribution. When the true random-effects distribution departs substantially from normality, estimated covariance parameters may still be informative about variability patterns, but inference can become sensitive and may require robust or alternative modeling strategies.

2 Covariance structure of random effects

2.1 Covariance matrix of random effects

The covariance structure is typically summarized by a covariance matrix \(G\) for the random-effects vector \(b_j\). Each entry in \(G\) corresponds to the covariance between two random-effect components within the same group.

2.1.1 Parameters and dimensionality

If a group-specific random-effects vector has dimension \(q\), then \(G\) is a \(q \times q\) symmetric matrix. The number of distinct covariance parameters is \[ \frac{q(q+1)}{2}, \] which grows quickly with the number of random-effect terms (e.g., multiple random slopes and intercepts).

2.1.1.1 Positive semidefiniteness and identifiability

Because \(G\) is a covariance matrix, it must be positive semidefinite: all variances must be nonnegative and the matrix must represent a feasible joint variability structure. Estimation routines enforce this through parameterizations (e.g., Cholesky factorization).

Identifiability refers to whether the data contain enough information to separate fixed-effect contributions from random-effect variance components and, when applicable, to distinguish among correlated random effects. Poor experimental design (such as insufficient replication within groups for certain random slopes) can make some covariance elements weakly identified or effectively non-estimable.

2.2 Off-diagonal elements and interpretation

Off-diagonal entries of \(G\) represent covariances between random effects. For two random effects \(b_{k}\) and \(b_{\ell}\), the covariance \(\operatorname{Cov}(b_k, b_\ell)\) indicates whether deviations in those components tend to occur together.

  • Positive covariance: groups showing larger-than-average values in one component also tend to show larger-than-average values in the other.
  • Negative covariance: larger deviations in one component are associated with smaller deviations in the other.
  • Near-zero covariance: deviations are largely uncorrelated at the random-effects level, though conditional relationships may still exist through shared covariates or residual structure.

Importantly, covariance does not imply deterministic dependence; it summarizes linear co-variation in the random-effects distribution.

2.3 Variance partitioning across hierarchy levels

In hierarchical models, covariance among observed outcomes arises from multiple sources: random effects at different grouping layers and residual variation within the lowest level. Variance partitioning expresses how much of the total variability is attributable to each layer. For example, in a two-level linear mixed model, outcomes from observations in the same group share the group-level random effects, while outcomes from different groups do not.

The random-effects covariance contributes directly to shared variability within groups and, depending on the design and model structure, can also affect marginal covariances between observations that load different random components.

2.4 Relationship to correlation

Covariance is scale-dependent, because it depends on the units of the random effects being modeled. Correlation expresses association standardized by the square roots of variances.

2.4.1 Converting covariance to correlation

For random effects \(b_k\) and \(b_\ell\), \[ \rho_{k\ell} = \frac{\operatorname{Cov}(b_k,b_\ell)}{\sqrt{\operatorname{Var}(b_k)\operatorname{Var}(b_\ell)}}. \] This conversion helps interpret the magnitude of association on a standardized [-1, 1] scale, making it easier to compare across different model specifications or across datasets with different measurement scales.

3 Implications for the marginal covariance of observations

3.1 From random effects covariance to outcome covariance

Even if residual errors are independent across observations, outcomes can be correlated because they share random effects. In linear mixed models, the covariance of two observations is driven by how their design vectors \(z_{ij}\) and \(z_{i'j'}\) overlap in group membership and random-effect loadings.

Conceptually, the random-effects covariance matrix \(G\) is “projected” into the space of observed outcomes through the random-effects design. Larger covariances among random-effect components tend to increase (or decrease) marginal association among outcomes when the corresponding random-effect components influence those outcomes jointly.

3.2 Within-group vs between-group dependence

Observations typically become correlated primarily within the same group because they share the same random-effects realization \(b_j\). Between-group covariances are often zero under the assumption that random effects for different groups are independent.

If the model uses multiple grouping structures, dependence may appear beyond a single cluster boundary, but in the simplest single-factor setting, shared group membership is the key driver of nonzero marginal covariance.

3.3 Cluster-level covariance patterns

Within a group, the covariance between two outcomes depends on:

  1. which random effects each observation loads (e.g., intercept-only vs intercept-and-slope),
  2. the magnitudes of the relevant covariates at each observation,
  3. the covariance structure \(G\) among random effects.

As a result, covariance patterns can vary across the within-group index \(i\) depending on covariate values. For example, if slopes are random and the slope deviation is correlated with the intercept deviation, the implied covariance among observations at different predictor levels can change sign or magnitude as predictor values vary.

3.4 Design-driven covariance (e.g., random slopes)

Random slopes introduce a dependence on the design matrix through the covariates. Two observations within the same cluster may have different implied correlation depending on their predictor values, because their random slope contributions scale with those values.

This means that the same estimated random-effects covariance matrix can produce different marginal covariance patterns under different study designs or coding choices for predictors, even when the underlying latent variability is unchanged.

4 Estimation of random-effects covariance

4.1 Likelihood-based estimation

For Gaussian linear mixed models, estimation typically relies on maximizing a likelihood or a related objective function with respect to fixed effects and covariance parameters. The covariance parameters in \(G\) influence the marginal covariance of the outcomes, so they enter the likelihood through the implied covariance matrix of \(y\).

Likelihood-based procedures provide coherent estimates of covariance components under model assumptions, and they allow comparison of different covariance structures through likelihood ratio tests, information criteria, or other selection rules.

4.2 Restricted maximum likelihood (REML)

Restricted maximum likelihood (REML) is widely used to estimate variance components in mixed models. The key idea is to base estimation on a likelihood that accounts for the estimation of fixed effects, which reduces bias in variance component estimates that can occur under ordinary maximum likelihood in finite samples.

REML commonly yields more reliable estimates of random-effects variances and covariances, especially when the number of groups is not large or when fixed effects consume substantial degrees of freedom.

4.3 Bayesian estimation and priors

Bayesian approaches treat covariance parameters as random variables with prior distributions. Priors can help stabilize estimation when data provide limited information about certain covariance elements (e.g., weakly identified off-diagonal terms).

Common strategies include:

  • priors directly on covariance matrix entries or on decompositions such as standard deviations and correlations,
  • hierarchical priors that borrow strength across related parameters,
  • sensitivity analyses to assess how conclusions change under different prior choices.

Bayesian estimation also produces full posterior uncertainty for covariance and derived correlation quantities, which can be valuable for reporting and interpretation.

4.4 Model checking for covariance adequacy

After estimating a covariance structure, it is important to assess whether the implied covariance patterns match observed features of the data. While diagnostics vary by software and outcome type, general principles include checking adequacy of variability and correlation implied by the model.

4.4.1 Diagnostics for variance and correlation parameters

Diagnostics may include:

  • examining residual and random-effect predictive checks,
  • comparing observed within-group variability to model-implied variability,
  • assessing whether estimated covariance parameters produce plausible shrinkage behavior,
  • using sensitivity analyses across alternative covariance structures (e.g., diagonal vs full covariance for random effects).

If the model systematically under- or over-predicts within-group differences, the covariance structure—particularly correlations among random intercepts and slopes—may require refinement.

5 Identifiability, constraints, and special cases

5.1 When covariance is not estimable

Some covariance parameters may not be estimable due to limited replication, collinearity in design terms, or insufficient overlap of covariate values that inform the relationship between random effects. For instance, if a random slope is effectively unobserved within many groups (or if all groups experience nearly the same predictor range), the model cannot reliably distinguish slope variability or its covariance with other random effects.

In such cases, estimators may produce unstable results, large standard errors, or boundary estimates (e.g., variances driven toward zero).

5.2 Boundary issues (near-zero variances)

Variance and covariance parameters can lie near the boundary of the parameter space, particularly when true variability is small. Optimizers may then report:

  • estimated variances at or extremely close to zero,
  • correlations that appear unstable because they involve division by small standard deviations.

Boundary behavior also affects the interpretation of off-diagonal covariance parameters: an estimated covariance may be small because either the association is weak or because one component’s variance is near zero, making correlation interpretations less meaningful.

5.3 Simplified covariance structures

To improve stability and interpretability, practitioners sometimes restrict \(G\) to simpler patterns.

5.3.1 Diagonal vs full covariance

A diagonal covariance structure assumes random effects are uncorrelated at the group level, setting off-diagonal entries to zero. A full covariance structure estimates covariances between random effects. Full structures can capture important dependence, but they may require more data to estimate reliably.

Choosing between diagonal and full covariance typically involves trade-offs: parsimony and stability versus flexibility to represent correlated heterogeneity.

5.3.2 Exchangeable/compound-symmetry patterns

In some contexts, covariance structures are designed to reflect symmetry assumptions. Exchangeable (compound-symmetry) patterns impose that all pairwise covariances among a set of related random components are equal. While not universally appropriate, such constraints reduce the number of parameters and can be justified when the grouping mechanism and design imply similar variability and dependence across random-effect components.

5.4 Centering and scaling effects

Centering and scaling covariates can influence the practical meaning and numerical stability of random intercept and random slope terms. Because the intercept and slope may be correlated in \(G\), changes in the predictor’s coding alter how the model partitions variation between these components.

A common consequence is that covariance between intercept and slope can change even if the implied relationship between baseline and predictor effect in the data remains similar. Centering (e.g., centering time or a continuous predictor) can improve interpretability by aligning the intercept with a meaningful reference value and can reduce correlation among estimated random effects in some parameterizations.

6 Extensions and variants

6.1 Multivariate random effects

Multivariate random-effects models introduce multiple correlated response variables or multiple latent effects that are linked by shared random-effects components. In these settings, the covariance structure must account for correlations across dimensions—either across different outcomes, different random-effect components, or both.

The covariance of random effects becomes a block matrix, where each block corresponds to covariances among components across outcomes and within a grouping factor.

6.2 Crossed random effects (multiple grouping factors)

Crossed designs arise when observations belong to combinations of multiple grouping factors rather than a single nested hierarchy. Examples include patients measured under multiple therapists or items rated by multiple users.

Crossed random effects require covariance modeling for each grouping factor and, when appropriate, allow correlations between random effects across factors. Covariance interpretation depends on whether random effects are assumed independent across factors or whether the model estimates dependencies.

6.3 Heteroscedastic and correlated residual models

While the covariance of random effects addresses between-group or subject-level heterogeneity, residual covariance models capture within-group correlation or changing variance across conditions. Heteroscedastic models allow the residual variance to differ across observations or predictor levels, and correlated residual models introduce dependence among residual errors.

These extensions interact with random-effects covariance: distinguishing random-effect-driven dependence from residual-driven dependence can be crucial. Without adequate residual structure, some covariance may be incorrectly attributed to random effects.

6.4 Non-Gaussian random effects considerations

When random effects are non-Gaussian, the covariance matrix still summarizes second-order dependence, but higher-order behavior may differ from the Gaussian case. Estimation can become more complex, especially for nonlinear link functions.

Approaches to address non-Gaussianity include flexible distributions, mixture models, or approximate inference methods. The goal is to represent realistic latent variability while retaining an interpretable covariance structure for co-variation among random effects.

7 Practical workflow and interpretation

7.1 Specifying the random-effects term

Model specification typically begins by choosing which grouping factors warrant random effects and which effects within groups should vary (intercepts, slopes, or both). The random-effects design matrices determine which covariates enter the random part and therefore which covariance elements in \(G\) become identifiable.

A common workflow uses incremental modeling: start with a simpler random-effects structure, then add correlations or additional random slopes if justified by theory, diagnostics, or improvement in predictive fit.

7.2 Interpreting estimated covariance and correlation

After estimation, the covariance matrix \(G\) is interpreted in the context of the random-effects design. Variances quantify group-to-group variability in each random effect, while covariances indicate whether deviations co-occur.

Because covariance depends on scale, many practitioners report both covariance and the derived correlations. Correlation is especially helpful to understand whether baseline heterogeneity and slope heterogeneity move together, but it should be interpreted alongside the magnitude of the underlying variances.

7.3 Reporting conventions in statistical output

Reporting often includes:

  • the estimated random-effects covariance matrix \(G\) or its components (variances and covariances),
  • derived correlations among random intercepts and slopes,
  • uncertainty measures such as standard errors or credible intervals,
  • the residual variance structure (since it affects marginal covariances).

Clear notation—indicating which random effects are included and the order of terms in covariance matrices—helps avoid ambiguity when comparing results across studies.

7.4 Common pitfalls and troubleshooting

Common issues include:

  • overparameterization, where a full covariance structure is attempted with insufficient data for certain covariance elements;
  • misinterpretation of correlation when associated variances are near zero;
  • confusing random-effects covariance with residual covariance, especially in repeated-measures contexts;
  • neglecting the influence of centering/scaling on intercept-slope covariance interpretation.

Troubleshooting often involves comparing alternative covariance structures (diagonal vs full; simpler parametrizations), checking whether random slopes are supported by within-group variation, and conducting model diagnostics focused on within-group dependence patterns.