Estimating GLS Variance Components by ML and REML

Why Estimate Variance Components?

In generalized least squares, the covariance matrix determines how observations are weighted and how their dependence is represented. When observations in the same group share a random intercept, their covariance depends on two unknown quantities: the variance of the shared group effect and the variance of observation-specific noise. Estimating these variance components is therefore part of fitting the covariance model, not a separate choice of weights made in advance.

Maximum likelihood (ML) and restricted maximum likelihood (REML) use the same covariance model but account differently for estimating fixed effects. In a balanced random-intercept model, the distinction has a transparent expression: ML estimates a between-group variance using \(G\) dimensions, whereas REML uses \(G-1\) group contrasts after removing the overall mean. That difference can affect the fitted covariance, particularly when the number of groups is modest.

Model and Covariance Structure

In vector notation, \(y=X\beta+Zb+\varepsilon\), so after integrating out the random intercepts, \(y\sim N(X\beta,V)\), where \(V=\tau^2ZZ'+\sigma^2I\). The matrix \(Z\) assigns each observation to its group. Two different observations in the same group have covariance \(\tau^2\); an observation has variance \(\tau^2+\sigma^2\). Observations in different groups have covariance zero under this model.

The likelihood-based calculations below use the balanced, intercept-only special case: each of \(G\) groups has \(m\) observations, and the only fixed effect is a common intercept \(\beta\). Write \(N=Gm\), and define \(\lambda=\sigma^2+m\tau^2\). This is the covariance eigenvalue in a group-constant direction; it is not the variance of an individual observation. The component constraint \(\tau^2\geq0\) is equivalent to \(\lambda\geq\sigma^2\).

Within-Group and Between-Group Variation

Let \(\bar y_g\) be the mean in group \(g\), and let \(\bar y\) be the overall mean. The data separate into deviations within groups and deviations among group means:

The factor \(m\) in \(SS_B\) matters: this is the sum of squares in the original \(N\)-observation space, not just the sum of squared deviations of the \(G\) group means. There are \(q=G(m-1)\) within-group directions, each with variance \(\sigma^2\), and \(G\) group-constant directions, each with variance \(\lambda\). One of those latter directions is the overall mean direction.

ML and REML Estimates

The ML procedure maximizes the likelihood of the observations, estimating \(\beta\) along with the covariance parameters. For the balanced intercept-only model, profiling out \(\beta\) gives \(\hat\beta=\bar y\). For \(m\geq2\), when \(SS_W>0\) and \(SS_B>0\), differentiating the displayed ML log likelihood with respect to each gives the finite unconstrained interior estimates

When \(q=G(m-1)>0\) and \(SS_W=0\), the likelihood is unbounded as \(\sigma^2\downarrow0\), so no finite ML or REML estimate exists under the strict constraint \(\sigma^2>0\). The denominator \(G\), rather than \(G-1\), in the second expression arises because the ordinary likelihood includes the determinant contribution from the overall mean direction, even though the fitted intercept makes its residual quadratic contribution zero. If these estimates satisfy \(\hat\lambda_{\mathrm{ML,unc}}\geq\hat\sigma^2_{\mathrm{ML,unc}}\), the corresponding variance-component estimates are \(\hat\tau^2_{\mathrm{ML}}=(\hat\lambda_{\mathrm{ML,unc}}-\hat\sigma^2_{\mathrm{ML,unc}})/m\) and \(\hat\sigma^2_{\mathrm{ML}}=\hat\sigma^2_{\mathrm{ML,unc}}\).

REML maximizes a likelihood formed from contrasts that eliminate the fixed effects. In this balanced case it removes the overall mean direction from the likelihood, leaving \(G-1\) between-group contrasts. Thus the displayed REML formula for \(\hat\lambda\) requires \(G\geq2\); if \(G=1\), there are no between-group contrasts and REML provides no information to estimate \(\tau^2\). For \(G\geq2\), when \(SS_W>0\) and \(SS_B>0\), the finite unconstrained interior estimates are

The different denominator is not a correction that can be applied to every ML estimate after the fact. It follows from the restricted likelihood's removal of the fixed-effect direction. In models with a general design matrix, the corresponding likelihoods are expressed using \(V\), \(X\), and the residual-forming projection; the balanced formulas above are a particularly clear special case.

Boundary Estimates and the Nonnegative Constraint

An unconstrained calculation can imply \(\hat\tau^2<0\), which is not an admissible variance. The fit must then be maximized subject to \(\tau^2\geq0\). In the balanced model, the boundary is \(\tau^2=0\), or equivalently \(\lambda=\sigma^2\). At that boundary, all observation directions have the same variance. The ML estimate is the residual sum of squares divided by \(N\), and the REML estimate is divided by \(N-1\):

These expressions can also be obtained by setting \(\lambda=\sigma^2\) in the respective likelihoods and maximizing over the common variance. For ML the log-variance coefficient is \(q+G=N\); for REML it is \(q+G-1=N-1\). If the unconstrained estimates violate \(\lambda\geq\sigma^2\), the constrained optimum in this balanced model is at that boundary.

From Components to a Covariance Matrix

Why the ML–REML Distinction Matters

ML and REML answer closely related but not identical likelihood questions. ML uses the likelihood of the full data and estimates fixed effects within that likelihood. REML uses information in contrasts that do not depend on the fixed effects being estimated. In the balanced random-intercept example, REML's between-group denominator is \(G-1\), reflecting removal of the overall mean direction; ML's denominator is \(G\). Neither calculation changes the assumed random-intercept structure.

REML likelihood values depend on the fixed-effect design \(X\), so they are generally not used to compare models with different fixed-effect specifications. ML likelihoods are used for such comparisons because they are based on the full data likelihood under each model. This distinction concerns likelihood comparison, not a claim that one method is universally preferable for every purpose.

Common Mistakes

• Using the REML denominator in an ML calculation. In this balanced model, the between-group denominator is \(G\) for ML and \(G-1\) for REML. State which likelihood is being maximized before calculating.
• Forgetting the factor \(m\) in \(SS_B\). The between-group sum of squares used with \(\lambda=\sigma^2+m\tau^2\) is \(m\sum_g(\bar y_g-\bar y)^2\).
• Reporting a negative variance component. An unconstrained difference \((\hat\lambda-\hat\sigma^2)/m<0\) is not an admissible estimate. Refit under \(\tau^2\geq0\), including the boundary.
• Calling \(\lambda\) the individual variance. Here \(\lambda\) is a group-constant-direction covariance eigenvalue. The individual variance is \(\sigma^2+\tau^2\).
• Interpreting a zero boundary estimate as proof of no group variation. It is the constrained likelihood estimate for the observed data, not proof that the underlying population component equals zero.