Solved – Mixed Models: How to derive Henderson’s mixed-model equations

linear algebramatrixmixed model

In the context of best linear unbiased predictors (BLUP), Henderson specified the mixed-model equations (see Henderson (1950): Estimation of Genetic Parameters. Annals of Mathematical Statistics, 21, 309-310). Let us assume the following mixed effects model:

$y = X\beta+Zu+e$

where $y$ is a vector of n observable random variables, $\beta$ is a vector of $p$ fixed effects, $X$ and $Z$ are known matrices, and $u$ and $e$ re vectors of $q$ and $n$ random effects such that $E(u) = 0$ and $E(e) = 0$ and

$ Var \begin{bmatrix}
u \\
e \\
\end{bmatrix} =
\begin{bmatrix}
G & 0 \\
0 & R \\
\end{bmatrix}\sigma^2$

where $G$ and $R$ are known positive definite matrices and $\sigma^2$ is a positive constant.

According to Henderson (1950), the BLUP estimates of $\hat {\beta}$ of $\beta$ and $\hat {u}$ of $u$ are defined as solutions to the following system of equation:

$X'R^{-1}X\hat {\beta}+X'R^{-1}Z\hat {u} = X'R^{-1}y$

$Z'R^{-1}X\hat {\beta}+(Z'R^{-1}Z + G^{-1})\hat {u} = Z'R^{-1}y$

(Also see: Robinson (1991): That BLUP is a good thing: the estimation of random effects (with discussion). Statistical Science, 6:15–51).

I have not found any derivation of this solution but assume that he approached it as follows:

$(y – X\beta – Zu)'V^{-1}(y – X\beta – Zu)$

where $V = R + ZGZ'$. Hence the solutions should therefore be

$X'V^{-1}X\hat {\beta} + X'V^{-1}Z\hat {u} = X'V^{-1}y$

$Z'V^{-1}X\hat {\beta} + Z'V^{-1}Z\hat {u} = Z'V^{-1}y$.

We also know that $V^{-1} = R^{-1} – R^{-1}Z(G^{-1}+Z'R^{-1}Z)Z'R^{-1}$.

However, ho to proceed to arrive at the mixed-model equations?

Best Answer

One approach is to form the log-likelihood and differentiate this with respect to the random effects $\mathbf{u}$ and set this equal to zero, then repeat, but differentiate with respect to the fixed effects $\boldsymbol{\beta}$.

With the usual normality assumptions we have:

$$ \begin{align*} \mathbf{y|u} &\sim \mathcal{N}\mathbf{(X\beta + Zu, R)} \\ \mathbf{u} &\sim \mathcal{N}(\mathbf{0, G}) \end{align*} $$ where $\mathbf{y}$ is the response vector, $\mathbf{u}$ and $\boldsymbol{\beta}$ are the random effects and fixed effects coefficient vectors $\mathbf{X}$ and $\mathbf{Z}$ are model matrices for the fixed effects and random effects respectively. The log-likelihood is then:

$$ -2\log L(\boldsymbol{\beta,}\mathbf{u}) = \log|\mathbf{R}|+(\mathbf{y - X\boldsymbol{\beta} - Zu})'\mathbf{R}^{-1}(\mathbf{y - X\boldsymbol{\beta} - Zu}) +\log|\mathbf{G}|+\mathbf{u'G^{-1}u} $$ Differentiating with respect to the random and fixed effects: $$ \begin{align*} \frac{\partial \log L}{\partial \mathbf{u}} &= \mathbf{Z'R^{-1}}(\mathbf{y - X\boldsymbol{\beta} - Zu}) - \mathbf{G^{-1}u} \\ \frac{\partial \log L}{\partial \boldsymbol{\beta}} &= \mathbf{X'R^{-1}}(\mathbf{y - X\boldsymbol{\beta} - Zu}) \end{align*} $$ After setting these both equal to zero, with some minor re-arranging, we obtain Henderson's mixed model equations:

$$ \begin{align*} \mathbf{Z'R^{-1}}\mathbf{y} &= \mathbf{Z'R^{-1}X\boldsymbol{\beta}} + \mathbf{u(Z'R^{-1}Z+G^{-1})} \\ \mathbf{X'R^{-1}}\mathbf{y} &= \mathbf{X'R^{-1}X\boldsymbol{\beta}} + \mathbf{X'R^{-1}Zu} \end{align*} $$

Related Solutions

Solved – Model Matrices for Mixed Effects Models

Creating the $J_i$ matrix entails producing 3 levels (309, 330 and 371) each one with 10 observations or measurements (nrow(sleepstudy[sleepstudy$Subject==309,]) [1] 10). Following the code in the original link in the OP:

f <- gl(3,10) Ji<-t(as(f,Class="sparseMatrix"))

enter image description here

Building the $X_i$ matrix can be aided by using the function getME as a reference:

library(lme4) sleepstudy <- sleepstudy[sleepstudy$Subject %in% c(309, 330, 371), ] rownames(sleepstudy) <- NULL fm1<-lmer(Reaction~Days+(Days|Subject), sleepstudy)

Xi <- getME(fm1,"mmList")

enter image description here

Since we will need the transpose, and the object Xi is not a matrix, the t(Xi) can be built as:

t_Xi <- rbind(c(rep(1,30)),c(rep(0:9,3)))

$Z_i$ is calculated as $Z_i= (J_i^{T}∗X_i^{T})^⊤$:

Zi<-t(KhatriRao(t_Ji,t_Xi)):

enter image description here

This corresponds to equation (6) in the original paper:

$$Z_i = (J_i^T*X_i^T)^T=\begin{bmatrix}J_{i1}^T\otimes X_{i1}^T\\J_{i2}^T\otimes X_{i2}^T\\\vdots\\J_{in}^T\otimes X_{in}^T\end{bmatrix}$$

And to see this we can instead play with truncated $J_i^T$ and $X_i^T$ matrices by imagining that instead of 9 measurements and a baseline (0), there is only 1 measurement (and a baseline). The resulting matrices would be:

$J_i^T=\left[\begin{smallmatrix}1&1&0&0&0&0\\0&0&1&1&0&0\\0&0&0&0&1&1\end{smallmatrix}\right]$ and $X_i^T=\left[\begin{smallmatrix}1&1&1&1&1&1\\0&1&0&1&0&1\end{smallmatrix}\right]$.

And

$J_i^T*X_i^T=\left[ \begin{smallmatrix} \left(\begin{smallmatrix}1\\0\\0\end{smallmatrix}\right)\otimes\left(\begin{smallmatrix}1\\0\end{smallmatrix}\right)&& \left(\begin{smallmatrix}1\\0\\0\end{smallmatrix}\right)\otimes\left(\begin{smallmatrix}1\\1\end{smallmatrix}\right)&& \left(\begin{smallmatrix}0\\1\\0\end{smallmatrix}\right)\otimes\left(\begin{smallmatrix}1\\0\end{smallmatrix}\right)&& \left(\begin{smallmatrix}0\\1\\0\end{smallmatrix}\right)\otimes\left(\begin{smallmatrix}1\\1\end{smallmatrix}\right)&& \left(\begin{smallmatrix}0\\0\\1\end{smallmatrix}\right)\otimes\left(\begin{smallmatrix}1\\0\end{smallmatrix}\right)&& \left(\begin{smallmatrix}0\\0\\1\end{smallmatrix}\right)\otimes\left(\begin{smallmatrix}1\\1\end{smallmatrix}\right) \end{smallmatrix} \right]$

$ \small= \begin{bmatrix}J_{i1}^T\otimes X_{i1}^T && J_{i2}^T \otimes X_{i2}^T && J_{i3}^T\otimes X_{i3}^T && J_{i4}^T\otimes X_{i4}^T && J_{i5}^T\otimes X_{i5}^T && J_{i6}^T\otimes X_{i6}^T\end{bmatrix}$

$=\left[\begin{smallmatrix}1&1&0&0&0&0\\0&1&0&0&0&0\\0&0&1&1&0&0\\0&0&0&1&0&0\\0&0&0&0&1&1\\0&0&0&0&0&1\end{smallmatrix}\right]$. Which transposed and extended would result in $Z_i=\left[\begin{smallmatrix}1&0&0&0&0&0\\1&1&0&0&0&0\\1&2&0&0&0&0\\\vdots\\\\0&0&1&0&0&0\\0&0&1&1&0&0\\0&0&1&2&0&0\\\vdots\\\\0&0&0&0&1&0\\0&0&0&0&1&1\\0&0&0&0&1&2\\\\\vdots\end{smallmatrix}\right]$.

Extracting the $b$ vector of random effects coefficients can be done with the function:

b <- getME(fm1,"b")

[1,] -44.1573839
[2,]  -2.4118590
[3,]  32.8633489
[4,]  -0.3998801
[5,]  11.2940350
[6,]   2.8117392

If we add these values to the fixed-effects of the call fm1<-lmer(Reaction~Days+(Days|Subject), sleepstudy) we get the intercepts:

205.3016 for 309; 282.3223 for 330; and 260.7530 for 371

and the slopes:

2.407141 for 309; 4.419120 for 330; and 7.630739 for 371

values consistent with:

library(lattice)
xyplot(Reaction ~ Days | Subject, groups = Subject, data = sleepstudy, 
       pch=19, lwd=2, type=c('p','r'))

enter image description here

$Zb$ can be calculated as as.matrix(Zi)%*%b.

Solved – The form of the Log-Likelihood Function in Mixed Linear Models

In standard maximum likelihood theory the likelihood function is based on observed data. Because the random effects are unobserved random variables, you need to integrate them out. When you do that you obtain the log-likelihood function you have written in Case 1.

The random effects are often estimated in a second step using empirical Bayes methodology. In particular, they are estimated via the posterior distribution

$$p(u \mid y; \hat \theta) = \frac{p(y \mid u; \hat \theta) \, p(u; \hat \theta)}{p(y; \hat \theta)},$$

where $p(\cdot)$ denotes the corresponding probability density functions, and $\hat \theta$ the maximum likelihood estimates obtained from the log-likelihood function of Case 1. Note that this is a whole distribution. As estimates of the random effects we use the mode of this distribution with respect to $u$. And here is where the joint likelihood of Case 2 comes into play. That is, the mode of the posterior distribution of the random effects is the mode of the joint density of Case 2 (the numerator in the equation above) because the normalizing constant in the denominator is not a function of $u$.

Best Answer

Related Solutions

Solved – Model Matrices for Mixed Effects Models

Solved – The form of the Log-Likelihood Function in Mixed Linear Models

Related Question