在前面的章节中，我们提到了 BLP 方法需要已知固定效应的值，而这在实践中往往不可能实现，由此引出 BLUP 方法，这应该也是本系列的最后一篇。

BLUP

我们先看 BLUP 方法的一般推导，假设我们想要预测一个随机变量 \(w\) ，我们只知道其均值为 \(\mathbf{k}^{\prime} \boldsymbol{\beta}\) ，方差为 \(v\) ，其与 \(\mathbf{y}^{\prime}\) 的协方差为 \(Cov(w, \mathbf{y}^{\prime}) = \mathbf{c}^{\prime}\) ，我们该如何取预测 \(w\) 呢？

一种可行的方案是我们找到一个关于 \(\mathbf{y}\) 的线性函数，其均值为 \(\mathbf{k}^{\prime} \boldsymbol{\beta}\) （即无偏），并且其在所有线性无偏估计值中具有最小的预测误差方差，这种方法我们称为最佳线性无偏预测法 (BLUP) 。

假设预测函数为 \(\mathbf{a'y}\) ，其期望为 \(E(\mathbf{a}^{\prime} \mathbf{y})=\mathbf{a}^{\prime} \mathbf{X} \boldsymbol{\beta}\) ，我们需要使之为 \(\mathbf{k}^{\prime} \boldsymbol{\beta}\) ，即 \(\mathbf{a}^{\prime} \mathbf{X} \boldsymbol{\beta} = \mathbf{k}^{\prime} \boldsymbol{\beta}\) 对于所有可能的 \(\boldsymbol{\beta}\) 恒成立，因此存在 \[ \mathbf{a}^{\prime} \mathbf{X}=\mathbf{k}^{\prime} \] 预测误差方差为下式，其中 \(v = Var(w)\) \[ Var(\mathbf{a'y}-w) = \mathbf{a'Va} - 2 \mathbf{a'c}+v \] 构建拉格朗日函数如下 \[ \mathbf{L} = 0.5 (\mathbf{a'Va} - 2 \mathbf{a'c}+v) + (\mathbf{a}^{\prime} \mathbf{X}-\mathbf{k}^{\prime}) \boldsymbol{\theta} \] 分别对 \(\mathbf{a}\) 和 \(\boldsymbol{\theta}\) 求偏导，使之为 0 ，得到 \[ \begin{aligned} &\partial \mathbf{L} / \partial \mathbf{a}= \mathbf{V} \mathbf{a} - \mathbf{c} +\mathbf{X} \boldsymbol{\theta}= \mathbf{0} \\ &\partial \mathbf{L} / \partial \boldsymbol{\theta}=\mathbf{X’a-k} = \mathbf{0} \end{aligned} \] 整理一下，我们得到方程组如下 \[ \left(\begin{array}{ll} \mathbf{V} & \mathbf{X} \\ \mathbf{X}^{\prime} & 0 \end{array}\right)\left(\begin{array}{l} \mathbf{a} \\ \boldsymbol{\theta} \end{array}\right)=\left(\begin{array}{l} \mathbf{c} \\ \mathbf{k} \end{array}\right) \] 这个方程组和广义最小二乘的方程组很相似。下面我们开始求解，对第一个方程进行变换，得到 $$ \[\begin{aligned} \mathbf{V} \mathbf{a}+\mathbf{X} \boldsymbol{\theta} &= \mathbf{c} \\ \mathbf{a} &= \mathbf{V^{-1}(c-X \boldsymbol{\theta}}) \\ \mathbf{X'a} &= \mathbf{X'V^{-1}(c-X \boldsymbol{\theta}}) \\ \mathbf{k} &= \mathbf{X'V^{-1}(c-X \boldsymbol{\theta}}) \\ \end{aligned}\]

\[ 因此我们得到了一个只关于 $\boldsymbol{\theta}$ 的方程组 \] = $$ 我们需要先证明这个方程组相容，我们只需要证明对于任何一个广义逆 \(\mathbf{ (X'V^{-1}X)^{-} }\) ，\(\mathbf{(X'V^{-1}X) (X'V^{-1}X)^{-}(X'V^{-1}c-k)=X'V^{-1}c-k}\) 均成立即可。假设 \(\mathbf{k}^{\prime} \boldsymbol{\beta}\) 可估计，则必然存在某个 \(\mathbf{a_{0}}\) 向量使得 \(\mathbf{k = X'a_{0}}\) ，因此上式可以写成 \(\mathbf{(X'V^{-1}X) (X'V^{-1}X)^{-}(X'V^{-1}c-X'a_{0})=X'V^{-1}c-X'a_{0}}\) 。我们设这里的 \(\mathbf{A=V^{-1/2}X}\) ，这里的 \(\mathbf{V^{-1/2}}\) 是 \(\mathbf{V^{1/2}}\) 的逆矩阵，易得 \(\mathbf{(X'V^{-1}X) (X'V^{-1}X)^{-}X'} = \mathbf{A'A(A'A)^{-}A'V^{1/2}} = \mathbf{A'V^{1/2}} = \mathbf{X'V^{-1/2} V^{1/2}} = \mathbf{X'}\) ，因此得证 \(\mathbf{(X'V^{-1}X) (X'V^{-1}X)^{-}(X'V^{-1}c-X'a_{0})=X'V^{-1}c-X'a_{0}}\) ，即我们证明这个方程组相容。

因此 \(\boldsymbol{\theta}\) 的一个解为 \[ \mathbf{(X'V^{-1}X)^{-}(X'V^{-1}c-k)} \] 我们将其带入到第一个方程 \(\mathbf{V} \mathbf{a}+\mathbf{X} \boldsymbol{\theta} = \mathbf{c}\) ，得到 \[ \begin{aligned} \mathbf{a}&=-\mathbf{V}^{-1} \mathbf{X} \boldsymbol{\theta}+\mathbf{V}^{-1} \mathbf{c} \\ &=\mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{k}-\mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{c}+\mathbf{V}^{-1} \mathbf{c} \end{aligned} \] 此时我们的预测值为 （这里需要设定 \(\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-}\) 同样是一个对称矩阵 ） \[ \mathbf{a}^{\prime} \mathbf{y}=\mathbf{k}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{y}+\mathbf{c}^{\prime} \mathbf{V}^{-1}\left[\mathbf{y}-\mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{y}\right] \] 由于 \(\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{y}=\boldsymbol{\beta}^{o}\) ，所以预测值可以写为 \[ \hat{w} = \mathbf{k}^{\prime} \boldsymbol{\beta}^{o}+\mathbf{c}^{\prime} \mathbf{V}^{-1}\left(\mathbf{y}-\mathbf{X} \boldsymbol{\beta}^{o}\right) \] 我们注意到，如果 \(\mathbf{k}^{\prime} \boldsymbol{\beta}=0\) 或者如果 \(\boldsymbol{\beta}\) 已知，那么此时的预测值就是 \(\mathbf{c}^{\prime} \mathbf{V}^{-1}\left(\mathbf{y}-\mathbf{X} \boldsymbol{\beta}^{o}\right)\) ，就是选择指数法的预测公式。因此，BLUP 就是用 \(\boldsymbol{\beta}^{0}\) 替换 \(\boldsymbol{\beta}\) 的 BLP 。

从元素推广到向量，得到 \[ \mathbf{\hat{w}} = \mathbf{K}^{\prime} \boldsymbol{\beta}^{o}+\mathbf{C}^{\prime} \mathbf{V}^{-1}\left(\mathbf{y}-\mathbf{X} \boldsymbol{\beta}^{o}\right) \] 其中，\(\mathbf{C} = Cov(\mathbf{y}, \mathbf{w}^{\prime})\) 。

方差及协方差

我们现在来看一些有用的方差和协方差，我们有 \[ \mathbf{\hat{w}}=\left(\mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} +\mathbf{C}^{\prime} \mathbf{V}^{-1}-\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \right) \mathbf{y} \] 假设 \(\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-}\) 是一个对称矩阵，并且是一个自反广义逆矩阵。因此我们有 $$ \[\begin{aligned} \operatorname{Var}(\hat{\mathbf{w}}) =& \left(\mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} +\mathbf{C}^{\prime} \mathbf{V}^{-1}-\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \right) \mathbf{V} \left(\mathbf{V^{-1}X(X'V^{-1}X)^{-}K + V^{-1}C - V^{-1}X(X'V^{-1}X)^{-}X'V^{-1}C}\right) \\ =& \left(\mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} +\mathbf{C}^{\prime} -\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \right) \left(\mathbf{V^{-1}X(X'V^{-1}X)^{-}K + V^{-1}C - V^{-1}X(X'V^{-1}X)^{-}X'V^{-1}C}\right) \\ =& \mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{K} + \mathbf{C}^{\prime}\mathbf{V^{-1}X(X'V^{-1}X)^{-}K} - \mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-}\mathbf{K} \\ & \quad +\mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V^{-1}C} + \mathbf{C'V^{-1}C} - \mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V^{-1}C} \\ & \quad -\mathbf{K'}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-}\mathbf{X'V^{-1}C}- \mathbf{C' V^{-1}X(X'V^{-1}X)^{-}X'V^{-1}C} + \mathbf{C'V^{-1}X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-}\mathbf{X'V^{-1}C} \\ =& \mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{K}+\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{C} \\ &-\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{C} \\ \end{aligned}\] \[ 协方差 \] \[\begin{aligned} \operatorname{Cov}\left(\hat{\mathbf{w}}, \mathbf{w}^{\prime}\right) =& \left(\mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} +\mathbf{C}^{\prime} \mathbf{V}^{-1}-\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \right) \operatorname{Cov}\left(\mathbf{y}, \mathbf{w}^{\prime}\right) \\ =& \left(\mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} +\mathbf{C}^{\prime} \mathbf{V}^{-1}-\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \right) \mathbf{C} \\ =& \mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{C}+\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{C} \\ &-\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{C} \\ =& \operatorname{Var}(\hat{\mathbf{w}}) \\ \end{aligned}\] \[ 预测误差方差为 \] \[\begin{aligned} \operatorname{Var}(\hat{\mathbf{w}}-\mathbf{w}) =& \operatorname{Var}(\mathbf{\hat{w}})-\operatorname{Cov}\left(\hat{\mathbf{w}}, \mathbf{w}^{\prime}\right)-\operatorname{Cov}\left(\mathbf{w}, \hat{\mathbf{w}}^{\prime}\right)+\operatorname{Var}(\mathbf{w}) \\ =& \mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{K}-\mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{C} \\ &-\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{K}-\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{C} \\ &+\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{C} +\mathbf{G} \end{aligned}\]

$$

BLUP 的性质

1. 在所有的线性无偏预测值中，BLUP 最大化了 \(\hat{w}_{i}\) 和 \(w_{i}\) 的相关
2. BLUP 具有不变性，即 \(\text{BLUP of }\mathbf{m^{\prime}w} = \mathbf{m^{\prime}\hat{w}}\)
3. 在正态分布假设下，我们有
   • \(E(\mathbf{w | \hat{w}}) = \mathbf{\hat{w}}\)
   • $Var( - ) = Var( | ) $
   • \(\mathbf{\hat{w}}\) 是最大似然估计值，并且是 \(\mathbf{w}\) 的条件均值的最佳线性无偏估计值
   • \(\mathbf{\hat{w}}\) 元素正确排序的可能性最大

混合模型方程组

上面的BLUP计算公式中涉及到 \(\mathbf{V}^{-1}\) 的计算，当观测值较多时， \(\mathbf{V}^{-1}\) 计算困难甚至不可计算，为此 Henderson 提出了求解 \(\mathbf{\hat{w}}\) 的另一种方法， 这个方法原来被 Henderson 称为最大似然法，但是实际上这个方法比最大似然方法更加具有普遍性，因为它其实不需要正态假设。这个方法就是混合模型方程组 (Mixed Model Equation, MME) .

对于一般的混合线性模型 \[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta}+ \mathbf{Zu}+ \mathbf{e} \]

均值和方差，协方差如下

\[ \begin{gathered} E(\mathbf{u})=\mathbf{0}, E(\mathbf{e})=\mathbf{0}, E(\mathbf{y})=\mathbf{X} \beta \\ \operatorname{Var}(\mathbf{u})=\mathbf{G}, \operatorname{Var}(\mathbf{e})=\mathbf{R}, \operatorname{Cov}\left(\mathbf{u}, \mathbf{e}^{\prime}\right)=\mathbf{0}, \operatorname{Var}(\mathbf{y})=\mathbf{Z G Z}^{\prime}+\mathbf{R} \end{gathered} \]

我们沿着 Henderson 一开始推出 MME 的思路，我们假设 \(\mathbf{u}\) 和 \(\mathbf{e}\) 服从正态分布， \(\mathbf{y}\) 和 \(\mathbf{u}\) 服从联合多变量高斯分布 \[ \mathbf{u} \sim N(\mathbf{0}, \mathbf{G}) \\ \mathbf{e} \sim N(\mathbf{0}, \mathbf{R}) \] 那么 \(\mathbf{y}\) 和 \(\mathbf{u}\) 的联合密度函数为 \[ f(\mathbf{y}, \mathbf{u})=f_{1}(\mathbf{y} \mid \mathbf{u}) f_{2}(\mathbf{u}) \] 根据多元正态分布，我们知道多元正态分布的条件分布 也是正态分布，其均值和协方差矩阵满足公式 \[ \begin{aligned} E(\mathbf{y} \mid \mathbf{x}) &=\boldsymbol{\mu}_{y}+\mathbf{\Sigma}_{y x} \Sigma_{x x}^{-1}\left(\mathbf{x}-\boldsymbol{\mu}_{x}\right), \\ \operatorname{cov}(\mathbf{y} \mid \mathbf{x}) &=\mathbf{\Sigma}_{y y}-\mathbf{\Sigma}_{y x} \mathbf{\Sigma}_{x x}^{-1} \mathbf{\Sigma}_{x y} . \end{aligned} \] 因此，我们有 $$ \[\begin{aligned} E(\mathbf{y} \mid \mathbf{u}) &=\boldsymbol{\mu}_{y}+\mathbf{\Sigma}_\mathbf{y u'} \Sigma_\mathbf{u u'}^{-1}\left(\mathbf{u}-E(\mathbf{u})\right) \\ &= \mathbf{X} \boldsymbol{\beta} + \mathbf{ZG G^{-1}u} \\ &= \mathbf{X} \boldsymbol{\beta} + \mathbf{Zu} \\ \operatorname{cov}(\mathbf{y} \mid \mathbf{u}) &=\mathbf{\Sigma}_{y y'}-\mathbf{\Sigma}_{y u'} \mathbf{\Sigma}_{u u'}^{-1} \mathbf{\Sigma}_{u y'} \\ &=\mathbf{ZGZ'+R} - \mathbf{(ZG)(G^{-1})(GZ')} \\ &=\mathbf{R} \\ \end{aligned}\] \[ 因此，我们有 \] N( + , ) \ \[ 根据多元正态分布的概率密度公式，我们有 \] \[\begin{aligned} &f_{1}(\mathbf{y} \mid \mathbf{u})=C_{1} \exp \left\{-\frac{1}{2}(\mathbf{y}-\mathbf{X} \beta-\mathbf{Z u})^{\prime} \mathbf{R}^{-1}(\mathbf{y}-\mathbf{X} \beta-\mathbf{Z u})\right\} \\ &f_{2}(\mathbf{u})=C_{2} \exp \left\{-\frac{1}{2} \mathbf{u}^{\prime} \mathbf{G}^{-1} \mathbf{u}\right\} \end{aligned}\]

\[ 其中 $C_{1}$ 和 $C_{2}$ 为常数，如下 \] C_{1}=(2 )^{-}|R|{-}, C_{2}=(2 )^{-}|G|{-} \[ 于是 \] f(, )=C {-(- -)^{} ^{-1}(- -)- ^{} ^{-1} } $$ 其中 \(C=C_{1} \times C_{2}\) 。

取对数得到 \[ \log f(\mathbf{y}, \mathbf{u})= \log (C) -\frac{1}{2}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta}-\mathbf{Z u})^{\prime} \mathbf{R}^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta} -\mathbf{Z u})-\frac{1}{2} \mathbf{u}^{\prime} \mathbf{G}^{-1} \mathbf{u} \] 根据最大似然的思想，我们要找一组参数 \(\boldsymbol{\hat{\beta}}\) 和 \(\mathbf{\hat{u}}\) 使得 \(\log f(\mathbf{y}, \mathbf{u})\) 最大。我们分别针对 \(\boldsymbol{\beta}\) 和 \(\mathbf{u}\) 求偏导，得到 \[ \begin{aligned} &\frac{\partial \log f(\mathbf{y}, \mathbf{u})}{\partial \boldsymbol{\beta}} =\mathbf{X'R^{-1}}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta}-\mathbf{Z u}) =\mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{y}-\mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{X} \beta-\mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{Z u} = \mathbf{0} \\ &\frac{\partial \log f(\mathbf{y}, \mathbf{u})}{\partial \mathbf{u}} =\mathbf{Z'R^{-1}}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta}-\mathbf{Z u}) -\mathbf{G}^{-1} \mathbf{u} =\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{y}-\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{X} \beta-\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z u}-\mathbf{G}^{-1} \mathbf{u} = \mathbf{0} \end{aligned} \] 整理一下，可得 \[ \begin{aligned} &\mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{X} \boldsymbol{\hat{\beta}}+\mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{Z} \hat{\mathbf{u}}=\mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{y} \\ &\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{X} \boldsymbol{\hat{\beta}}+\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right) \hat{\mathbf{u}}=\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{y} \end{aligned} \] 我们也可以整理成方程组的形式 \[ \left[\begin{array}{cc} \mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{X} & \mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{Z} \\ \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{X} & \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1} \end{array}\right]\left[\begin{array}{l} \boldsymbol{\hat{\beta}} \\ \hat{\mathbf{u}} \end{array}\right]=\left[\begin{array}{l} \mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{y} \\ \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{y} \end{array}\right] \] 此方程组就称为混合模型方程组。

我们下面证明混合模型方程组与根据 BLUP 的定义直接得到的 \(\boldsymbol{\hat{\beta}}\) 和 \(\mathbf{\hat{u}}\) 等价，也就是说混合模型方程组实际不需要正态假设。

证明等价

由混合模型方程组的第二个方程得 \[ \mathbf{u}=\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \mathbf{R}^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\hat{\beta}}) \] 这里需要证明 \(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\) 可逆，由于 \(\mathbf{G}\) 阵是正定矩阵，因此其逆矩阵 \(\mathbf{G^{-1}}\) 也是正定矩阵，易证 \(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}\) 为半正定矩阵，因此 \(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\) 也是正定矩阵，因此其可逆。

将该式带入到第一个方程得到 \[ \begin{aligned} &\mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{X} \boldsymbol{\hat{\beta}}-\mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{Z}\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{X} \boldsymbol{\hat{\beta}} \\ &=\mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{y}-\mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{Z}\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{y} \end{aligned} \] 变换一下得到 $$ \[\begin{aligned} \mathbf{X'}\left(\mathbf{R}^{-1}-\mathbf{R}^{-1} \mathbf{Z}\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \mathbf{R}^{-1}\right)\mathbf{X}\boldsymbol{\hat{\beta}} &= \mathbf{X'}\left(\mathbf{R}^{-1}-\mathbf{R}^{-1} \mathbf{Z}\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \mathbf{R}^{-1}\right)\mathbf{y} \\ \mathbf{X'}\mathbf{V^{-1}}\mathbf{X}\boldsymbol{\hat{\beta}} &= \mathbf{X'}\mathbf{V^{-1}}\mathbf{y} \\ \end{aligned}\]

$$ 其中有 \(\mathbf{V}^{-1}=\mathbf{R}^{-1}-\mathbf{R}^{-1} \mathbf{Z}\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \mathbf{R}^{-1}\) ，推导可见混合线性模型一之广义最小二乘 。

因此，混合模型方程组得到的 \(\boldsymbol{\hat{\beta}}\) 就是广义最小二乘估计值。

然后我们再需要证明 \(\mathbf{\hat{u}} = \mathbf{C}^{\prime} \mathbf{V}^{-1}\left(\mathbf{y}-\mathbf{X} \boldsymbol{\hat{\beta}}\right)\) ，即等于选择指数法估计值。 \[ \begin{aligned} \mathbf{u} &=\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \mathbf{R}^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\hat{\beta}}) \\ &=\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{V} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X} \widehat{\beta}) \\ &=\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \mathbf{R}^{-1}\left(\mathbf{Z G Z}^{\prime}+\mathbf{R}\right) \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\hat{\beta}}) \\ &=\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1}\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right) \mathbf{G Z}^{\prime} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\hat{\beta}}) \\ &=\mathbf{G Z}^{\prime} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\hat{\beta}}) \\ &=\mathbf{C}^{\prime} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\hat{\beta}}) \quad( \because \mathbf{C}=\mathbf{Z G}) \end{aligned} \] 得证。

因此我们证明了混合模型方程组与根据 BLUP 的定义直接得到的 \(\boldsymbol{\hat{\beta}}\) 和 \(\mathbf{\hat{u}}\) 等价，所以在家畜育种中，混合模型方程组已成了 BLUP 法的同义词。

但是这种从最大似然法推导出来混合模型方程组的方式，和 BLUP 的定义似乎缺少概念上的直接联系，下面我们尝试直接从 BLUP 的定义出发，直接推导得到混合模型方程组。

重新推导

根据 BLUP 定义，对于单个变量 \(w\) ，我们有方程组 \[ \left(\begin{array}{ll} \mathbf{V} & \mathbf{X} \\ \mathbf{X}^{\prime} & 0 \end{array}\right)\left(\begin{array}{l} \mathbf{a} \\ \boldsymbol{\theta} \end{array}\right)=\left(\begin{array}{l} \mathbf{c} \\ \mathbf{k} \end{array}\right) \] 推广至向量 \(\mathbf{w}\) ，我们有下面的方程组 (注意这里的 \(\boldsymbol{\theta}\) 是一个矩阵，上面是一个向量 )。矩阵 \(\mathbf{A}\) 满足 \(\mathbf{\hat{w} = A'y}\) 。 \[ \left(\begin{array}{ll} \mathbf{V} & \mathbf{X} \\ \mathbf{X}^{\prime} & 0 \end{array}\right)\left(\begin{array}{l} \mathbf{A} \\ \boldsymbol{\theta} \end{array}\right)=\left(\begin{array}{l} \mathbf{C} \\ \mathbf{K} \end{array}\right) \] 将 \(\mathbf{V = ZGZ' + R}\) 和 \(\mathbf{C=ZG}\) 带入第一个方程，得到 \[ \mathbf{(ZGZ' + R)A} + \mathbf{X}\boldsymbol{\theta} = \mathbf{CG} \\ \] 或者写成 \[ \mathbf{RA} + \mathbf{ZG(Z'A-I)} + \mathbf{X}\boldsymbol{\theta} = \mathbf{0} \] 令 \(\mathbf{T=G(Z'A-I)}\) ，则 \(\mathbf{Z'A-G^{-1}T=I}\) ，于是有方程组 \[ \mathbf{RA} + \mathbf{ZT} + \mathbf{X}\boldsymbol{\theta} = \mathbf{0}\\ \mathbf{Z'A-G^{-1}T=I} \\ \mathbf{X'A=K} \] 由其中的第一个方程可得 \[ \mathbf{A} = -\mathbf{R^{-1}X\boldsymbol{\theta} - R^{-1}ZT} \] 带入到第二方程，得 \[ \begin{aligned} -\mathbf{Z'R^{-1}X\boldsymbol{\theta} - Z'R^{-1}ZT-G^{-1}T} &= \mathbf{I} \end{aligned} \] 带入到第三方程，得 \[ -\mathbf{X'R^{-1}X\boldsymbol{\theta} - X’R^{-1}ZT} = \mathbf{K} \] 整理成方程组得格式为 \[ \left[\begin{array}{ll} \mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{X} & \mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{Z} \\ \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{X} & \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1} \end{array}\right]\left[\begin{array}{c} \boldsymbol{\theta} \\ \mathbf{T} \end{array}\right]=\left[\begin{array}{c} -\mathbf{K} \\ -\mathbf{I} \end{array}\right] \] 前面前面的推导，我们知道这个方程组应该是相容的。我们将上式系数矩阵的一个广义逆表示为（假设这个广义逆也是一个对称矩阵，满足 \(\mathbf{C_{12}} = \mathbf{C_{21}}^{\prime}\) ，\(\mathbf{C_{11}}\) 和 \(\mathbf{C_{22}}\) 也是对称矩阵） \[ \left[\begin{array}{ll} \mathbf{C}_{11} & \mathbf{C}_{12} \\ \mathbf{C}_{21} & \mathbf{C}_{22} \end{array}\right] \] 则该方程组的解为 \[ \left[\begin{array}{l} \boldsymbol{\theta} \\ \mathbf{T} \end{array}\right]=\left[\begin{array}{ll} \mathbf{C}_{11} & \mathbf{C}_{12} \\ \mathbf{C}_{21} & \mathbf{C}_{22} \end{array}\right]\left[\begin{array}{l} -\mathbf{K} \\ -\mathbf{I} \end{array}\right] \] 于是 $$ \[\begin{aligned} \mathbf{A} &= \left[\begin{array}{l} \mathbf{R^{-1}X} & \mathbf{R^{-1}Z} \end{array}\right] \left[\begin{array}{l} -\boldsymbol{\theta} \\ -\mathbf{T} \end{array}\right] \\ \\ &= \left[\begin{array}{l} \mathbf{R^{-1}X} & \mathbf{R^{-1}Z} \end{array}\right] \left[\begin{array}{ll} \mathbf{C}_{11} & \mathbf{C}_{12} \\ \mathbf{C}_{21} & \mathbf{C}_{22} \end{array}\right]\left[\begin{array}{l} \mathbf{K} \\ \mathbf{I} \end{array}\right] \end{aligned}\] \[ 因此 \] \[\begin{aligned} \mathbf{\hat{w}} &= \mathbf{A'y} \\ &= \left[\begin{array}{l} \mathbf{K'} & \mathbf{I} \end{array}\right] \left[\begin{array}{ll} \mathbf{C}_{11} & \mathbf{C}_{12} \\ \mathbf{C}_{21} & \mathbf{C}_{22} \end{array}\right] \left[\begin{array}{l} \mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{y} \\ \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{y} \end{array}\right] \end{aligned}\]

\[ 同时根据上面的推导，我们有 \] = ^{} +^{} ^{-1}(- )=^{} + = \[ 因此 \] =\[ 即 \] =$$ 证明完毕。

性质

假设 \(\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-}\) 是一个对称矩阵，并且是一个自反广义逆矩阵。

固定效应估计值

固定效应估计值就是一个广义最小二乘估计值，公式如下，其性质不再赘述。 \[ \boldsymbol{\hat{\beta}}=\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{y} \]

随机效应估计值

随机效应计算公式为 \[ \begin{aligned} \mathbf{\hat{u}} &=\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \mathbf{R}^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\hat{\beta}}) \\ &=\mathbf{G Z}^{\prime} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\hat{\beta}}) \\ \end{aligned} \] 根据广义最小二乘估计值的性质，我们易得 \(\mathbf{X} \boldsymbol{\hat{\beta}}\) 是一个可估计函数，即 \(\mathbf{X} \boldsymbol{\hat{\beta}}\) 的值唯一，因此 \(\mathbf{\hat{u}}\) 也是唯一的。

首先看期望 \[ \begin{aligned} E(\hat{\mathbf{u}}) &\left.=\mathbf{G Z}^{\prime} \mathbf{V}^{-1}\left(\mathbf{X} \boldsymbol{\beta}-\mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X} \boldsymbol{\beta}\right)\right) \\ &=\mathbf{G Z}^{\prime} \mathbf{V}^{-1}(\mathbf{X} \boldsymbol{\beta}-\mathbf{X} \boldsymbol{\beta})=0 \end{aligned} \] 其中存在 \(\mathbf{X} \left( \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X} = \mathbf{X}\) ，证明如下，我们设 \(\mathbf{A=V^{-1/2}X}\) ，这里的 \(\mathbf{V^{-1/2}}\) 是 \(\mathbf{V^{1/2}}\) 的逆矩阵，易得 \[ \mathbf{X} \left( \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X} = \mathbf{V^{1/2}A(A'A)^{-}A'A} = \mathbf{V^{1/2}A} = \mathbf{X} \] 第二，我们看随机效应的协方差矩阵，根据上面的 BLUP的一般推导，我们有 \[ \operatorname{Var}(\hat{\mathbf{w}}) = \mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{K}+\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{C} -\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{C} \\ \] 这里 \(\mathbf{u}\) 均值为 \(\mathbf{0}\) ，即 \(\mathbf{K = 0}\) ，同时有 \(\mathbf{C=ZG}\) 带入得到 \[ \begin{aligned} \operatorname{Var}(\hat{\mathbf{u}}) &= \mathbf{GZ'} \mathbf{V}^{-1} \mathbf{ZG} -\mathbf{GZ'} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{ZG} \\ & =\mathbf{G} \mathbf{Z}^{\prime} \mathbf{V}^{-1}\left(\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1}\right) \mathbf{Z G} \end{aligned} \] 根据上面的 BLUP的一般推导，我们有 \[ \operatorname{Cov}\left(\hat{\mathbf{u}}, \mathbf{u}^{\prime}\right) = \operatorname{Var}(\hat{\mathbf{u}}) \] 固定效应的可估计函数与随机效应估计值的协方差矩阵为 $$ \[\begin{aligned} \operatorname{Cov}\left(\mathbf{T} \boldsymbol{\hat{\beta}}, \mathbf{\hat{u}}^{\prime}\right) &= \mathbf{T}\operatorname{Cov}\left(\boldsymbol{\hat{\beta}}, \mathbf{\hat{u}}^{\prime}\right) \\ &= \mathbf{T}\operatorname{Cov}\left(\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{y}, \left( \mathbf{G Z}^{\prime} \mathbf{V}^{-1}(\mathbf{I}-\mathbf{X} \left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} ) \mathbf{y} \right)^{\prime}\right) \\ &=\mathbf{T} \left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{V} \left( \mathbf{I} - \mathbf{V^{-1}X\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} X'} \right) \mathbf{V^{-1}ZG} \\ &=\mathbf{T} \left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \left( \mathbf{I} - \mathbf{V^{-1}X\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} X'} \right) \mathbf{V^{-1}ZG} \\ &=\mathbf{T} \left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V^{-1}ZG} - \mathbf{T} \left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V^{-1}ZG} \\ &=\mathbf{0} \end{aligned}\] \[ 根据通式 \] \[\begin{aligned} \operatorname{Var}(\hat{\mathbf{w}}-\mathbf{w}) =& \operatorname{Var}(\mathbf{\hat{w}})-\operatorname{Cov}\left(\hat{\mathbf{w}}, \mathbf{w}^{\prime}\right)-\operatorname{Cov}\left(\mathbf{w}, \hat{\mathbf{w}}^{\prime}\right)+\operatorname{Var}(\mathbf{w}) \\ =& \mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{K}-\mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{C} \\ &-\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{K}-\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{C} \\ &+\mathbf{C}^{\prime} \mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{C} +\mathbf{G} \end{aligned}\] \[ 我们有 \] \[\begin{aligned} \operatorname{Var}(\hat{\mathbf{u}}-\mathbf{u}) &= -\mathbf{GZ'V^{-1}ZG} + \mathbf{GZ'}\mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{ZG} + \mathbf{G} \\ &=\mathbf{G}-\mathbf{G Z}^{\prime} \mathbf{V}^{-1}\left(\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1}\right) \mathbf{Z G} \\ \end{aligned}\]

$$

系数矩阵广义逆矩阵的性质

假设混合模型方程组的系数矩阵的任意一个对称的广义逆矩阵为 \[ \left[\begin{array}{ll} \mathbf{C}_{11} & \mathbf{C}_{12} \\ \mathbf{C}_{12}^{\prime} & \mathbf{C}_{22} \end{array}\right] \] 则存在下式 \[ \begin{aligned} &\operatorname{Var}(\boldsymbol{\hat{\beta}})=\mathbf{ C}_{11} \\ &\operatorname{Var}(\hat{\mathbf{u}})=\operatorname{Cov}(\hat{\mathbf{u}}, \mathbf{u})=\mathbf{G}-\mathbf{C}_{22} \\ &\operatorname{Var}(\hat{\mathbf{u}}-\mathbf{u})=\mathbf{C}_{22} \\ &\operatorname{Cov}\left( \boldsymbol{\hat{\beta}}, \quad \hat{\mathbf{u}}^{\prime}\right)=\mathbf{0} \\ &\operatorname{Cov}( \boldsymbol{\hat{\beta}}, \quad \mathbf{u^{\prime}})=-\mathbf{ C}_{12} \\ &\operatorname{Cov}( \boldsymbol{\hat{\beta}}, \quad \hat{\mathbf{u}}^{\prime}-\mathbf{u}^{\prime})=\mathbf{ C}_{12} \end{aligned} \] 下面的证明是基于系数矩阵满秩的前提下，即此时的广义逆就是逆矩阵存在下式成立。 $$

= $$ 这是 Henderson 在 1975 年提出的，我们有

\[ \left(\begin{array}{l} \boldsymbol{\hat{\beta}} \\ \hat{\mathbf{u}} \end{array}\right)=\left[\begin{array}{lll} \mathbf{C}_{11} & & \mathbf{C}_{12} \\ \mathbf{C}_{12}^{\prime} & & \mathbf{C}_{22} \end{array}\right]\left[\begin{array}{l} \mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{y} \\ \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{y} \end{array}\right]=\left(\begin{array}{l} \mathbf{Q}_{1}^{\prime} \\ \mathbf{Q}_{2}^{\prime} \end{array}\right) \mathbf{y} \] 其中，我们设 \[ \left(\begin{array}{l} \mathbf{Q}_{1}^{\prime} \\ \mathbf{Q}_{2}^{\prime} \end{array}\right) = \left[\begin{array}{lll} \mathbf{C}_{11} & & \mathbf{C}_{12} \\ \mathbf{C}_{12}^{\prime} & & \mathbf{C}_{22} \end{array}\right]\left[\begin{array}{l} \mathbf{X}^{\prime} \mathbf{R}^{-1} \\ \mathbf{Z}^{\prime} \mathbf{R}^{-1} \end{array}\right] \] 继续拆开，得到 \[ \begin{aligned} \mathbf{Q}_{1}^{\prime} &= \mathbf{C}_{11} \mathbf{X}^{\prime} \mathbf{R}^{-1} + \mathbf{C}_{12} \mathbf{Z}^{\prime} \mathbf{R}^{-1} \\ \mathbf{Q}_{2}^{\prime} &= \mathbf{C}_{12}^{\prime} \mathbf{X}^{\prime} \mathbf{R}^{-1} + \mathbf{C}_{22} \mathbf{Z}^{\prime} \mathbf{R}^{-1} \\ \end{aligned} \] 或者 \[ \begin{aligned} \mathbf{Q}_{1} &= \mathbf{R}^{-1}\mathbf{X} \mathbf{C}_{11} + \mathbf{R}^{-1}\mathbf{Z} \mathbf{C}_{12}^{\prime} \\ \mathbf{Q}_{2} &= \mathbf{R}^{-1}\mathbf{X} \mathbf{C}_{12} + \mathbf{R}^{-1}\mathbf{Z} \mathbf{C}_{22} \\ \end{aligned} \] 因此，我们有 \[ \begin{aligned} \boldsymbol{\hat{\beta}} &= \mathbf{Q}_{1}^{\prime} \mathbf{y} \\ \mathbf{\hat{u}} &= \mathbf{Q}_{2}^{\prime} \mathbf{y} \\ \end{aligned} \] 我们可以得到下面这个等式 \[ \begin{aligned} \left[\begin{array}{l} \mathbf{Q}_{1}^{\prime} \\ \mathbf{Q}_{2}^{\prime} \end{array}\right][\mathbf{X} \quad \vdots \quad \mathbf{Z}] &= \left[\begin{array}{lll} \mathbf{C}_{11} & & \mathbf{C}_{12} \\ \mathbf{C}_{12}^{\prime} & & \mathbf{C}_{22} \end{array}\right]\left[\begin{array}{l} \mathbf{X}^{\prime} \mathbf{R}^{-1} \\ \mathbf{Z}^{\prime} \mathbf{R}^{-1} \end{array}\right] [\mathbf{X} \quad \vdots \quad \mathbf{Z}] \\ &=\left[\begin{array}{cc} \mathbf{C}_{11} & \mathbf{C}_{12} \\ \mathbf{C}_{12}^{\prime} & \mathbf{C}_{22} \end{array}\right]\left(\left[\begin{array}{cc} \mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{X} & \mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{Z} \\ \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{X} & \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1} \end{array}\right]-\left[\begin{array}{cc} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{G}^{-1} \end{array}\right]\right) \\ &=\left(\begin{array}{cc} \mathbf{I} & \mathbf{0} \\ \mathbf{0} & \mathbf{I} \end{array}\right)-\left(\begin{array}{cc} \mathbf{0} & \mathbf{C}_{12} \mathbf{G}^{-1} \\ \mathbf{0} & \mathbf{C}_{22} \mathbf{G}^{-1} \end{array}\right) \end{aligned} \] 因此我们得到 \[ \begin{aligned} &\mathbf{Q_{1}^{\prime} X=I} \\ &\mathbf{Q_{1}^{\prime} Z=-C_{12} G^{-1}} \\ &\mathbf{Q_{2}^{\prime} X=0} \\ &\mathbf{Q_{2}^{\prime} Z=I-C_{22} G^{-1}} \end{aligned} \] 利用这些结果，我们得到 $$ \[\begin{aligned} \operatorname{Var}(\boldsymbol{\hat{\beta}}) &=\operatorname{Var}(\mathbf{Q}_{1}^{\prime} \mathbf{y})\\ &=\mathbf{Q}_{1}^{\prime}\left(\mathbf{R}+\mathbf{Z} \mathbf{G} \mathbf{Z}^{\prime}\right) \mathbf{Q}_{1} \\ &=\mathbf{Q}_{1}^{\prime}\mathbf{R}\mathbf{Q}_{1}+ \mathbf{Q}_{1}^{\prime} \mathbf{Z} \mathbf{G} \mathbf{Z}^{\prime} \mathbf{Q}_{1} \\ &=\mathbf{Q}_{1}^{\prime}\mathbf{R} \left( \mathbf{R}^{-1}\mathbf{X} \mathbf{C}_{11} + \mathbf{R}^{-1}\mathbf{Z} \mathbf{C}_{12}^{\prime} \right)+ \mathbf{Q}_{1}^{\prime} \mathbf{Z} \mathbf{G} \mathbf{Z}^{\prime} \mathbf{Q}_{1} \\ &=\mathbf{Q}_{1}^{\prime} \mathbf{X C}_{11}+\mathbf{Q}_{1}^{\prime} \mathbf{Z C}_{12}^{\prime}+(\mathbf{Q}_{1}^{\prime} \mathbf{Z}) \mathbf{ G } (\mathbf{Q}_{1}^{\prime} \mathbf{Z})^{\prime} \\ &=\mathbf{I C}_{11}-\mathbf{C}_{12} \mathbf{G}^{-1} \mathbf{C}_{12}^{\prime}+ (-\mathbf{C}_{12} \mathbf{ G }^{-1}) \mathbf{ G } (-\mathbf{C}_{12} \mathbf{ G }^{-1})^{\prime} \\ &=\mathbf{ C}_{11}-\mathbf{C}_{12} \mathbf{G}^{-1} \mathbf{C}_{12}^{\prime}+ \mathbf{C}_{12} \mathbf{G}^{-1} \mathbf{C}_{12}^{\prime} \\ &=\mathbf{C}_{11} \end{aligned}\]

$$

\[ \begin{aligned} \operatorname{Cov}\left(\boldsymbol{\hat{\beta}}, \hat{\mathbf{u}}^{\prime}\right) &=\operatorname{Cov}\left( \mathbf{Q}_{1}^{\prime} \mathbf{y}, (\mathbf{Q}_{2}^{\prime} \mathbf{y})^{\prime}\right) \\ &=\mathbf{Q}_{1}^{\prime}\left(\mathbf{R}+\mathbf{ZGZ}^{\prime}\right) \mathbf{Q}_{2} \\ &=\mathbf{Q}_{1}^{\prime}\mathbf{R}\mathbf{Q}_{2}+\mathbf{Q}_{1}^{\prime}\mathbf{ZGZ}^{\prime}\mathbf{Q}_{2} \\ &=\mathbf{Q}_{1}^{\prime}\mathbf{R} \left( \mathbf{R}^{-1}\mathbf{X} \mathbf{C}_{12} + \mathbf{R}^{-1}\mathbf{Z} \mathbf{C}_{22}\right) +(\mathbf{Q}_{1}^{\prime}\mathbf{Z})\mathbf{G}(\mathbf{Q}_{2}^{\prime}\mathbf{Z})^{\prime} \\ &=\mathbf{Q}_{1}{ }^{\prime} \mathbf{X C}_{12}+\mathbf{Q}_{1}{ }^{\prime} \mathbf{Z C}_{22}+(\mathbf{Q}_{1}^{\prime}\mathbf{Z})\mathbf{G}(\mathbf{Q}_{2}^{\prime}\mathbf{Z})^{\prime} \\ &=\mathbf{I C}_{12}-\mathbf{C}_{12} \mathbf{G}^{-1} \mathbf{C}_{22}-\mathbf{C}_{12} \mathbf{G}^{-1} \mathbf{G}\left(\mathbf{I}-\mathbf{G}^{-1} \mathbf{C}_{22}\right) \\ &= \mathbf{0} \end{aligned} \]

\[ \begin{aligned} \operatorname{Cov}\left(\boldsymbol{\hat{\mathbf{\beta}}}, \mathbf{u}^{\prime}\right) &= \operatorname{Cov}\left(\mathbf{Q}_{1}^{\prime} \mathbf{y}, \mathbf{u}^{\prime}\right) \\ &=\mathbf{Q}_{1}^{\prime} \mathbf{ZG} \\ &=-\mathbf{C}_{12} \mathbf{G}^{-1} \mathbf{G} \\ &=-\mathbf{C}_{12} \end{aligned} \]

\[ \begin{aligned} \operatorname{Cov}\left(\boldsymbol{\hat{\beta}}, \mathbf{\hat{u}^{\prime}-\mathbf{u}^{\prime}}\right) &= \operatorname{Cov}\left(\boldsymbol{\hat{\beta}}, \mathbf{\hat{u}^{\prime}}\right) - \operatorname{Cov}\left(\boldsymbol{\hat{\beta}}, \mathbf{u}^{\prime}\right) \\ &=0-\left(-\mathbf{C}_{12}\right) \\ &=\mathbf{C}_{12} \end{aligned} \]

$$ \[\begin{aligned} \operatorname{Var}(\hat{\mathbf{u}}) &= \operatorname{Var}(\mathbf{Q}_{2}^{\prime} \mathbf{y}) \\ &=\mathbf{Q}_{2}^{\prime}\left(\mathbf{R}+\mathbf{Z G Z}^{\prime}\right) \mathbf{Q}_{2} \\ &=\mathbf{Q}_{2}^{\prime} \mathbf{R} \mathbf{Q}_{2} + (\mathbf{Q}_{2}^{\prime} \mathbf{Z}) \mathbf{G} (\mathbf{Q}_{2}^{\prime} \mathbf{Z})^{\prime} \\ &=\mathbf{Q}_{2}^{\prime} \mathbf{R} \left( \mathbf{R}^{-1}\mathbf{X} \mathbf{C}_{12} + \mathbf{R}^{-1}\mathbf{Z} \mathbf{C}_{22} \right) + \left( \mathbf{I-C_{22} G^{-1}} \right) \mathbf{G} \left( \mathbf{I- G^{-1} C_{22}} \right) \\ &=\mathbf{Q}_{2}^{\prime} \mathbf{X C}_{12}+\mathbf{Q}_{2}^{\prime} \mathbf{Z C}_{22} + \left( \mathbf{G-C_{22} } \right) \left( \mathbf{I- G^{-1} C_{22}} \right) \\ &=\mathbf{0}+\left(\mathbf{I}-\mathbf{C}_{22} \mathbf{G}^{-1}\right) \mathbf{C}_{22}+ \left( \mathbf{G-C_{22} } \right) \left( \mathbf{I- G^{-1} C_{22}} \right) \\ &=\mathbf{C}_{22}-\mathbf{C}_{22} \mathbf{G}^{-1}\mathbf{C}_{22} + \mathbf{G-C_{22} } - \mathbf{C}_{22} +\mathbf{C}_{22} \mathbf{G}^{-1}\mathbf{C}_{22}\\ &=\mathbf{G}-\mathbf{C}_{22} \end{aligned}\]

$$

\[ \begin{aligned} \operatorname{Cov}(\hat{\mathbf{u}}, \mathbf{u}) &=\operatorname{Cov}( \mathbf{Q}_{2}^{\prime} \mathbf{y} , \mathbf{u}) \\ &=\mathbf{Q}_{2}^{\prime} \mathbf{Z G} \\ &=\left(\mathbf{I}-\mathbf{C}_{22} \mathbf{G}^{-1}\right) \mathbf{G} \\ &=\mathbf{G}-\mathbf{C}_{22} \end{aligned} \]

\[ \begin{aligned} \operatorname{Var}(\hat{\mathbf{u}}-\mathbf{u}) &=\operatorname{Var}(\hat{\mathbf{u}}) - \operatorname{Cov}(\hat{\mathbf{u}}, \mathbf{u}) - \operatorname{Cov}(\mathbf{u}, \hat{\mathbf{u}}) + \operatorname{Var}(\mathbf{u}) \\ &=\mathbf{G}-\mathbf{C}_{22}-2\left(\mathbf{G}-\mathbf{C}_{22}\right)+\mathbf{G} \\ &=\mathbf{C}_{22} \end{aligned} \]

全部证明完毕。

简化方程组

在家畜育种的实际运用中，我们常有 \[ \operatorname{Var}(\mathbf{e})=\mathbf{R}=\mathbf{I} \sigma_{e}^{2}, \quad \operatorname{Var}(\mathbf{u})=\mathbf{G}=\mathbf{A} \sigma_{u}^{2} \] 此时，我们可以对混合模型方程组的两侧均消去 \(\mathbf{R^{-1}}\) ，得到一个简化版本 \[ \left[\begin{array}{ll} \mathbf{X}^{\prime} \mathbf{X} & \mathbf{X}^{\prime} \mathbf{Z} \\ \mathbf{Z}^{\prime} \mathbf{X} & \mathbf{Z}^{\prime} \mathbf{Z}+k \mathbf{A}^{-1} \end{array}\right]\left[\begin{array}{l} \boldsymbol{\hat{\beta}} \\ \mathbf{\hat{u}} \end{array}\right]=\left[\begin{array}{c} \mathbf{X^{\prime} y} \\ \mathbf{Z^{\prime} y} \end{array}\right] \] 其中 \[ k=\sigma_{e}^{2} / \sigma_{u}^{2} \]

优势

混合模型方程组中不涉及 \(\mathbf{V}\) 矩阵，在常规评估中， \(\mathbf{A^{-1}}\) 计算比 \(\mathbf{V^{-1}}\) 容易得多，当然这主要得益于 Henderson 提出了一种不用求逆直接构建 \(\mathbf{A^{-1}}\) 的方法。

参考文献

1. Henderson，1975，Best Linear Unbiased Estimation and Prediction under a Selection Model
2. Henderson，《Applications_of_Linear_Models_in_Animal_Breeding》, 1984
3. 张沅，张勤，《畜禽育种中的线性模型》