Chapter 20 REML Estimation of Variance Components

Consider the GLM: \[ y = X\beta + \epsilon, \text{ where } \epsilon \sim N(0, \Sigma) \] and \(\Sigma\) is an \(n\times n\) positive definite variance matrix which depends on an unknown parameter vector \(\gamma\).

Residual/Restricted Maximum Likelihood (REML): an approach that produces unbiased estimators for these special cases and produces less biased estimates than ML in general.

REML Method

  1. Find \(n - \text{rank}(X) = n - r\) linearly independent vectors \(a_1, \ldots, a_{n-r}\) such that \(a_i'X = 0'\) for all \(i = 1, \ldots, n-r\).

  2. Find the MLE of \(\gamma\) using \(w_1 \equiv a_1'y, \ldots, w_{n-r} \equiv a_{n-r}'y\) as data. \[ A = [a_1, \ldots, a_{n-r}] \quad w = \begin{bmatrix} w_1 \\ \vdots \\ w_{n-r} \end{bmatrix} = \begin{bmatrix} a_1'y \\ \vdots \\ a_{n-r}'y \end{bmatrix} = A'y \]

    • If \(a'X = 0'\), \(a'y\) is known as error contrast. Thus, \(w_1, \ldots, w_{n-r}\) comprise a set of \(n-r\) error constrasts

    • There exist \(n-r\) linearly independent rows of \(I - P_X\) to get \(a_1, \ldots, a_{n-r}\).

    • The error constrasts \(w_1, \ldots, w_{n-r}\) will be a subset of the elements of the residual vector \((I-P_X)y = y - \hat y\)

Note that \(w = A'y = A'\epsilon \sim N(0, A'\Sigma A)\) and the distribution of \(w\) depends only on \(\gamma\).

The log likelihood function in this case is \[ \ell(\gamma\mid w) = -\frac{1}{2}\log |A'\Sigma A| - \frac{1}{2}w'(A'\Sigma A)^{-1}w - \frac{n-r}{2}\log(2\pi) \] An MLE of \(\gamma\) can be found in the general case using numerical methods to obtain the REML estimate of \(\gamma\).

In 611, we show that every set of \(n-r\) linearly independent error contrasts yields the same REML estimator of \(\gamma\).

In any case, once a REML estimate of \(\gamma\) has been obtained, the BLUE of an estimate \(C\beta\) can be approximated by \(C\hat\beta_{\hat\Sigma} = C(X'\hat \Sigma^{-1} X)^{-}X' \hat\Sigma^{-1} y\), where \(\hat\Sigma\) is \(\Sigma\) with \(\hat\gamma\) (the REML estimate of \(\gamma\)) in place of \(\gamma\).

Example in Chapter 19
lme(SeedlingWeight ̃Genotype,random= ̃1|Tray,data=d)

Linear mixed-effects model fit by REML
  Data: d
  Log-restricted-likelihood: -123.5705
  Fixed: SeedlingWeight  ̃ Genotype
(Intercept)   Genotype2
  15.288838   -3.550201
Random effects:
 Formula:  ̃1 | Tray
        (Intercept) Residual
StdDev:    3.414856  1.88223
Number of Observations: 56
Number of Groups: 8


lmer(SeedlingWeight ̃Genotype+(1|Tray),data=d)

Linear mixed model fit by REML [’lmerMod’]
Formula: SeedlingWeight  ̃ Genotype + (1 | Tray)
   Data: d
REML criterion at convergence: 247.1411
Random effects:
 Groups   Name        Std.Dev.
 Tray   (Intercept)   3.415
 Residual             1.882
Number of obs: 56, groups: Tray, 8
Fixed Effects:
(Intercept)    Genotype2
      15.29        -3.55

In SAS, REML is also the default for SAS proc mixed.