A Review of Mixed Effect Model 1 - REML & BLUE & BLUP

Posted on 2023-01-18 Edited on 2023-01-21 In Stat Model Word count in article: 1.1k Reading time ≈ 4 mins.

The 1st post to review the Mixed Effect Model.

Background

Consider the general linear model \[ y=X\beta+\epsilon \text{, where } \epsilon \sim N(0,\Sigma) \] and $\Sigma$ is an $n \times n$ positive definite variance matrix that depends on unknown parameters that are organized in a vector $\gamma$.

We know that maximum likelihood (ML) estimation of the parameter vector $\gamma$ (variance component vector) can be biased.

For example, suppose

$$ \Sigma = \sigma^2 \begin{bmatrix} 1 & \rho & \rho^2 \\ \rho & 1 & \rho\\ \rho^2 & \rho & 1 \end{bmatrix} $$

and

$$ \gamma = \begin{bmatrix} \sigma^2 \\ \rho \end{bmatrix} $$

where $\sigma^2 \gt 0$ and $\rho \in (-1,1)$.

More details could be found at STAT 510 Lecture Note MLE.

A series of lecture notes for STAT 510

Example of MLE Bias

For $\Sigma=\sigma^2I$, where $\gamma=\sigma^2$, the MLE of $\sigma^2$ is \[ \frac{(y-X\hat{\beta})^{'}(y-X\hat{\beta})}{n} \] with expectation \[ \frac{n-r}{n}\sigma^2. \] This MLE for $\sigma^2$ is criticized for "failing to account for the loss of degrees of freedom needed to estimate $\beta$."

REML

REML is an approach that produces unbiased estimators for these special cases and produces less biased estimates than ML in general.

Depending on whom you ask, REML stands for REsidual Maximum Likelihood or REstricted Maximum Likelihood.

For linear mixed effects models, the REML estimators of variance components produce the same estimates as the unbiased ANOVA-based estimators formed by taking appropriate linear combinations of mean squares when the latter are positive and data are balanced.

In any case, once a REML estimate of $\gamma$ (and thus $\Sigma$) has been obtained, the BLUE (Best Linear Unbiased Estimate) of an estimable $C\beta$ can be approximated by \[ C\hat{\beta}_\hat{\Sigma}=C(X^{'}\hat{\Sigma}^{-1}X)^{-1}X^{'}\hat{\Sigma}^{-1}y, \] where $\hat{\Sigma}$ is $\Sigma$ with $\hat{\gamma}$ (the REML estimate of $\gamma$) in place of $\gamma$.

A detailed REML method (algorithm) could be found at REML Method.

REML is the default method for lmer in R;
REML is also the default for PROC MIXED in SAS.

BLUE v.s. BLUP

Best Linear Unbiased Estimates (BLUEs) are the solutions (or estimates) associated with the fixed effects.

Best Linear Unbiased Predictions (BLUPs) are the solutions (but identified as predictions) associated with the random effects of a model.

But what does Best Linear Unbiased mean?

Best	Among all possible unbiased linear estimators these solutions have minimum variance
Linear	Solutions are formed from a linear combination of the observations
Unbiased	Expectations of these solutions are equal to their true values

The use of BLUPs to predict random effects was first described by C. R. Henderson. He developed the Mixed Model Equations (MME, equation [2]) for Linear Mixed Models (LMM) in order to calculate BLUPs of breeding values (or any random effect) and BLUEs of fixed effects.

Let’s have a look at how solutions are obtained from a linear mixed model. For this, we will consider the model written in matrix notation:

[1] $y=Xb+Zu+e$

where $y$, $b$, $u$, and $e$ are vectors of observations, fixed effects, random effects, and random residuals, respectively; and $X$ and $Z$ are design matrices connecting the observations to the effects.

The above model is the basis for most genetic analyses with what is known as the Animal Model, that once fitted provides us with an estimation of the breeding values for each of the ‘animals’ (individuals) considered in the model.

Considering the Animal Model on its matrix notation of above, we have $b$ corresponding to a series of nuisance fixed effects, for example contemporary group, age or replicate. However, what is more interesting to us is $u$, corresponding to the breeding values. As this factor is random, we have distributional assumptions, namely that $u \sim MVN(0,\sigma^2_aA)$ where the $A$ is the numerator relationship matrix that describes the additive genetic relationship between individuals, and $\sigma^2_a$ is the additive variance.

As the error or residual term, $e$, is also a random effect, we have their distributional assumptions of $e \sim MVN(0,\sigma^2_eI)$, with $I$ an identity matrix and $\sigma^2_e$ the residual variance (or mean square error).

As you can see, the central problem in predicting breeding values from observed phenotypic data is separating the genetic and environmental effects. This separation is clearly done by the Animal Model.

So, in order to get our solutions (i.e., BLUEs and BLUPs), we need to solve the following system of equations, known as Henderson’s MME:

[2]

$$ \begin{vmatrix} X'X & X'Z \\ Z'X & Z'Z+A^{-1}a \end{vmatrix} \begin{vmatrix} b \\ u \end{vmatrix}= \begin{vmatrix} X'y \\ Z'y \end{vmatrix} $$

All the elements we have previously defined, except for $a=\sigma^2_e/\sigma^2_a$, which is formed from variance components estimated in the first step when fitting a LMM. This system can involve hundreds or even thousands of equations, requiring some intense computational calculations in order to estimate our solutions.

A different way to see the above equations is to express the formulae for BLUEs as:

[3] $b=(X'V^{−1}X)^{−1}X'V^{−1}y$

and for the BLUPs (breeding values) as:

[4] $u=\sigma^2_aAZ'V^{−1}(y−Xb)$

where $V=\sigma^2_aZAZ'+\sigma^2_eI$

Interestingly, expression [3] is equivalent to the estimation of Weighted Least Squares (WLS) as described for linear regression or linear models. Here, the weights are defined as $W=V^{−1}$. This is relevant as the elements on $V$ (particularly the diagonals) will be associated with the amount of information, or weight, that each record has (in this case phenotypic response of an individual).

Also, the expression for the breeding values is extremely relevant. Note that [4] can be written as:

$u=H(y−Xb)$

where $H$ is $\sigma^2_aAZ'V^{-1}$ resembling the calculation of heritability with the genetic component in the numerator ($\sigma^2_aAZ'$) and the phenotypic variance ($V$) in the denominator as an inverse. Also, $y−Xb$, represents the phenotypic response ‘corrected’ for nuisance effects. Therefore, the above expression resembles the breeder’s equation:

$gi=h^2S_i=h^2(y_i−μ)$

where the index ii is used to identify the $i$th individual.

Therefore, fitting a LMM using an Animal Model is equivalent to that well-known expression, but in this case, the specification of matrices connecting pieces of data allows us to use all information and relationships between individuals and to have different weights for each of these pieces. Therefore, Henderson’s MME provide us with the best linear unbiased predictions of our breeding values!

Source: BLUEs, BLUPs and Breeding Values, what is the difference?