这是混合线性模型理论系列的第一篇，内容为概述及广义最小二乘（固定效应的估计）。

混合线性模型概述

前言

统计分析的目的是推断 (inference) 和预测 (prediction)，推断指估计参数的大小（比如确定哪些自变量有作用，哪些自变量没有作用），预测是对未来的观测值进行预测。

假设我们定义我们的观测值向量是一个 \(n\) 维的向量 \(\mathbf{y}\) ，一般来说我们将观测值向量视为对真实的或概念上的总体的一次随机抽样，即向量 \(\mathbf{y}\) 是一个随机向量，向量 \(\mathbf{y}\) 的每一个元素是随机变量 。但是我们很少知道总体的分布，一般我们可能会假设总体服从多元正态分布。

模型

混合线性模型表示如下 \[ \mathbf{y}=\mathbf{X} \boldsymbol{\beta}+\mathbf{Zu}+\mathbf{e} \] 其中，

\(\mathbf{X}\) 是一个已知的固定的 \(n \times p \quad (n > p)\) 的矩阵， 其秩为 \(r \leq p\) 。 \(\boldsymbol{\beta}\) 是一个固定的由未知参数组成的 \(p \times 1\) 向量 。 \(\mathbf{Z}\) 是一个已知的固定的 \(n \times q\) 矩阵。 \(\mathbf{u}\) 是一个 \(q \times 1\) 的随机向量。 \(\mathbf{e}\) 是由随机误差组成的 \(n \times 1\) 的随机向量

注意：这里的 \(\mathbf{y, u ,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} \] 这里的协方差矩阵 \(\mathbf{G, R}\) 可能是已知的，也可能是未知的。由于 \(\mathbf{u}\) 和 \(\mathbf{e}\) 均无法观测，因此我们一般不知道 \(\mathbf{G}\) 和 \(\mathbf{R}\) ，但是我们可以对这两个协方差矩阵做一些假设，从而减少估计的参数数目。我们一般假设 \(\mathbf{R}\) 的非对角线元素均为 0，对角线元素均相同，也就是不同样本具有相同的残差并且彼此的残差之间不相干，即假设 \(\mathbf{R} = \sigma_{e}^{2} \mathbf{I}\) 。

其中 \(\boldsymbol{\beta}\) 称为固定效应，因为其元素均为固定的不变的常数，只是我们不知道具体的值；\(\boldsymbol{\mu}\) 称为随机效应，因为它是从总体得到的一次随机抽样，其元素均为随机变量，我们无法完全控制其大小。

除了随机误差外，完全由固定效应组成的模型称为固定效应模型，或固定模型；除群体平均数 \(\mu\) 之外，完全由随机效应组成的模型称为随机效应模型，简称为随机模型。

其实所有的线性模型都可以表示为混合线性模型，固定模型和随机模型均可以视为混合线性模型的一个特例。这样一来我们可以通过这一个模型的通式来解决所有问题（线性问题）。

在混合线性模型中，我们需要估计固定效应 \(\boldsymbol{\beta}\) 和随机效应 \(\boldsymbol{\mu}\) ，我们将固定效应的估计称为估计，将随机效应的估计称为预测。

下面我们看如何通过广义最小二乘法来估计固定效应 \(\boldsymbol{\beta}\) 。

广义最小二乘

在普通最小二乘中，我们假设 \(\operatorname{Var}(\mathbf{y}) = \operatorname{Var}(\mathbf{e}) = \sigma_{e}^{2} \mathbf{I}\) ，此时得到的正规方程组 \(\mathbf{X^{T}X } \boldsymbol{\hat{\beta}} = \mathbf{X^{T}y}\) 及满秩情况下的唯一解 \(\boldsymbol{\beta}_{\mathrm{LS}}=\left(\boldsymbol{X}^{\mathrm{T}} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\mathrm{T}} \boldsymbol{y}\) 即为最佳线性无偏估计值 (best linear unbiased estimators, BLUE) 。这一部分内容可见 矩阵微分与正规方程组推导 。

但是如果 \(\mathbf{y}\) 的协方差矩阵并不满足普通最小二乘的假设，此时如果我们仍采用普通最小二乘方法，此时得到的参数估计值就不再是 BLUE ，此时我们需要对 \(\mathbf{y}\) 作相应的线性变换。

我们考虑一般的情况，假设固定模型如下 \[ \mathbf{y} = \boldsymbol{X \beta} + \mathbf{e} \] 假设\(E(\mathbf{e}) = \mathbf{0}\)， \(\operatorname{Var}(\mathbf{y}) = \operatorname{Var}(\mathbf{e}) = \mathbf{V}=\sigma^{2} \mathbf{R}\) ，其中 \(\mathbf{R}\) 为一个 \(n \times n\) 的相关系数矩阵，这里 \(\mathbf{R}\) 是一个正定矩阵（ 因此 \(\mathbf{R}^{-1}\) 同样是正定矩阵）。根据正定矩阵的性质，一定存在一个非奇异矩阵 \(\mathbf{Q}\) 使得 \[ \mathbf{QQ^{\prime}= R^{-1}} \] 用 \(\mathbf{Q}\) 对 \(\mathbf{y}\) 作线性变换，有 \[ \begin{aligned} \mathbf{y}^{*} &=\mathbf{Q y} \\ &=\mathbf{Q} \mathbf{X} \boldsymbol{\beta}+\mathbf{Q} \mathbf{e} \\ &=\mathbf{X}^{*} \boldsymbol{\beta}+\mathbf{e}^{*} \end{aligned} \] 其中 \(\mathbf{X}^{*}=\mathbf{Q} \mathbf{X}, \quad \mathbf{e}^{*}=\mathbf{Q} \mathbf{e}\) ，则 \[ \operatorname{Var}\left(\mathbf{y}^{*}\right)=\operatorname{Var}(\mathbf{Q y})=\mathbf{Q} \operatorname{Var}(\mathbf{y}) \mathbf{Q}^{\prime}=\sigma^{2} \mathbf{Q} \mathbf{R} \mathbf{Q}^{\prime} \] 因为 \(\mathbf{QQ^{\prime}= R^{-1}}\) ，对该式左乘 \(\mathbf{Q}^{-1}\) ，右乘 \(\mathbf{Q}^{'-1}\) ，得到 \(\mathbf{Q}^{-1} \mathbf{R}^{-1} \mathbf{Q}^{\prime-1}=\mathbf{I}\) ，即 \(\mathbf{Q R Q}^{\prime}=\mathbf{I}\) ，于是 \[ \operatorname{Var}\left(\mathbf{y}^{*}\right)=\sigma^{2} \mathbf{I} \] 因此，我们知道线性转换后的模型 \(\mathbf{y}^{*}=\mathbf{X}^{*} \boldsymbol{\beta}+\mathrm{e}^{*}\) 符合普通最小二乘的假设，此时我们得到最小二乘方程组如下 \[ \mathbf{X}^{*'} \mathbf{X}^{*}\boldsymbol{\hat{\beta}}=\mathbf{X}^{*'} \mathbf{y}^{*} \] 将 \(\mathbf{X}^{*}=\mathbf{Q} \mathbf{X}\) ， $^{*} = $ 带入，得到 $$ \[\begin{aligned} &\mathbf{X}^{*'} \mathbf{X}^{*}\boldsymbol{\hat{\beta}}=\mathbf{X}^{*'} \mathbf{y}^{*} \\ &\mathbf{X'Q'} \mathbf{QX} \boldsymbol{\hat{\beta}}=\mathbf{X'Q'} \mathbf{Qy} \\ &\mathbf{X'R^{-1}} \mathbf{X} \boldsymbol{\hat{\beta}}=\mathbf{X'R^{-1}} \mathbf{y} \\ \end{aligned}\]

$$ \(\mathbf{X'R^{-1}} \mathbf{X} \boldsymbol{\hat{\beta}}=\mathbf{X'R^{-1}} \mathbf{y}\) 称为 广义最小二乘方程组。 当 \(\mathbf{X}\) 满秩时，我们得到唯一的广义最小二乘估计值 \(\boldsymbol{\hat{\beta}}= (\mathbf{X'R^{-1}} \mathbf{X})^{-1} \mathbf{X'R^{-1}} \mathbf{y}\) ，也就是最佳线性无偏估计值；当 \(\mathbf{X}\) 不满秩时，我们有无穷个解，某个解可以表示为 \(\boldsymbol{\hat{\beta}}= (\mathbf{X'R^{-1}} \mathbf{X})^{-} \mathbf{X'R^{-1}} \mathbf{y}\) ，其中 \((\mathbf{X'R^{-}} \mathbf{X})^{-}\) 为任意一个广义逆，此时若 \(\boldsymbol{q' \beta}\) 为可估计函数，则 \(\boldsymbol{q' \hat{\beta}}\) 是对 \(\boldsymbol{q' \beta}\) 的最佳线性无偏估计。

对上面的广义最小二乘方程组的左右两侧均乘以 \(1/\sigma^{2}\) ，我们得到该方程组的另一种形式。 \[ \mathbf{X'V^{-1}} \mathbf{X} \boldsymbol{\hat{\beta}}=\mathbf{X'V^{-1}} \mathbf{y} \] 上面我们推导了固定模型的情况，我们可以轻松推广至一般的混合线性模型，设混合线性模型表达式同上，我们可以将随机效应和随机误差效应视为一体，即令 \(\mathbf{Zu + e = e^{*}}\) ，此时模型写为 \[ \mathbf{y} = \boldsymbol{X \beta} + \mathbf{e}^{*} \] 此时，\(\operatorname{Var}(\mathbf{y}) = \operatorname{Var}(\mathbf{e^{*}}) =\mathbf{Z G Z}^{\prime}+\mathbf{R}\) ，即此时 \(\mathbf{V} =\mathbf{Z G Z}^{\prime}+\mathbf{R}\) ，同样有 \(\mathbf{X'V^{-1}} \mathbf{X} \boldsymbol{\hat{\beta}}=\mathbf{X'V^{-1}} \mathbf{y}\) 方程组成立。

另一种证明

我们接着延续高斯-马尔可夫定理的思想，从所有的线性无偏估计值中找到一个“最佳的”估计值，这里“最佳”的意思是指估计值具有最小的方差。

假设 \(\mathbf{k}^{\prime} \boldsymbol{\beta}\) 是一个可估计函数，即我们可以找到一个向量 \(\mathbf{a}\) 使得 \(\mathrm{E}\left(\mathbf{a}^{\prime} \mathbf{y}\right)=\mathbf{k}^{\prime} \boldsymbol{\beta}\) 。我们知道 \(E\left(\mathbf{a}^{\prime} \mathbf{y}\right)=\mathbf{a}^{\prime} \mathbf{X} \boldsymbol{\beta}\) ，因此 \(\mathbf{a}^{\prime} \mathbf{X} \boldsymbol{\beta} = \mathbf{k}^{\prime} \boldsymbol{\beta}\) 需要对所有可能的 \(\boldsymbol{\beta}\) 均成立，那么当且仅当下式成立时 \(\mathbf{k}^{\prime} \boldsymbol{\beta}\) 是一个可估计函数 \[ \mathbf{a}^{\prime} \mathbf{X}=\mathbf{k}^{\prime} \] 同时，我们有 \[ \operatorname{Var}\left(\mathbf{a}^{\prime} \mathbf{y}\right)=\mathbf{a}^{\prime}(\operatorname{Var}(\mathbf{y})) \mathbf{a}=\mathbf{a}^{\prime} \mathbf{V a} \] 这里我们假设 \(\mathbf{V}\) 矩阵已知。

因此我们是在满足约束条件 \(\mathbf{a}^{\prime} \mathbf{X}=\mathbf{k}^{\prime}\) 的前提下使得 \(\mathbf{a'Va}\) 的对角线元素最小，构建拉格朗日函数如下 \[ \mathbf{L}=\frac{1}{2}\mathbf{tr(a'Va) + (a'X-k')\boldsymbol{\theta}} \] 求偏导数，使之为 \(\mathbf{0}\) 得到 \[ \begin{aligned} &\partial \mathbf{L} / \partial \mathbf{a}= \mathbf{V} \mathbf{a}+\mathbf{X} \boldsymbol{\theta}= \mathbf{0} \\ &\partial \mathbf{L} / \partial \boldsymbol{\theta}=\mathbf{X’a-k} = \mathbf{0} \end{aligned} \] 将两个式子整理成一个方程组，得到 \[ \left(\begin{array}{cc} \mathbf{V} & \mathbf{X} \\ \mathbf{X}^{\prime} & \mathbf{0} \end{array}\right)\left(\begin{array}{l} \mathbf{a} \\ \boldsymbol{\theta} \end{array}\right)=\left(\begin{array}{l} \mathbf{0} \\ \mathbf{k} \end{array}\right) \] 对第一个方程进行变换可得 $$ \[\begin{aligned} \mathbf{V} \mathbf{a}+\mathbf{X} \boldsymbol{\theta} &= \mathbf{0} \\ \mathbf{a} &= \mathbf{-V^{-1}X \boldsymbol{\theta}} \\ \mathbf{X'a} &= \mathbf{-X'V^{-1}X \boldsymbol{\theta}} \\ \mathbf{k} &= \mathbf{-X'V^{-1}X \boldsymbol{\theta}} \\ \end{aligned}\]

\[ 因此我们得到了一个只关于 $\boldsymbol{\theta}$ 的方程组 \] = $$ 我们只需要证明对于任何一个广义逆 \(\mathbf{ (X'V^{-1}X)^{-} }\) ，\(\mathbf{(X'V^{-1}X) (X'V^{-1}X)^{-}k=k}\) 均成立，则该方程组相容。因为 \(\mathbf{k}^{\prime} \boldsymbol{\beta}\) 可估计，因此必然存在某个 \(\mathbf{a_{0}}\) 向量使得 \(\mathbf{k = X'a_{0}}\) ，因此上式可以写成 \(\mathbf{(X'V^{-1}X) (X'V^{-1}X)^{-}X'a_{0}=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'a_{0}=X'a_{0}}\) 。

因此 \(\boldsymbol{\theta}\) 的一个解为 \[ \mathbf{-(X'V^{-1}X)^{-}k} \] 带入到第一个方程 \(\mathbf{V} \mathbf{a}+\mathbf{X} \boldsymbol{\theta}= \mathbf{0}\) ，得到 \[ \mathbf{a} = \mathbf{V^{-1}X (X'V^{-1}X)^{-}k} \] 我们可以进一步证明无论 \(\mathbf{ (X'V^{-1}X)^{-} }\) 如何取值，\(\mathbf{a}\) 均不变。同样设 \(\mathbf{A=V^{-1/2}X}\) ，\(\mathbf{a} = \mathbf{V^{-1}X (X'V^{-1}X)^{-}k} = \mathbf{V^{-1}X (X'V^{-1}X)^{-}X'a_{0}}\) ，则 \(\mathbf{a = V^{-1/2} A(A'A)^{-}A'V^{1/2}a_{0}}\) 。根据广义逆的性质，无论 \((\mathbf{A}^{\prime} \mathbf{A})^{-}\) 取何值， \(\mathbf{A}\left(\mathbf{A}^{\prime} \mathbf{A}\right)^{-} \mathbf{A}^{\prime}\) 均保持不变，因此 \(\mathbf{a}\) 保持不变。

此时 \(\mathbf{k}^{\prime} \boldsymbol{\beta}\) 的 BLUE 为 \[ \mathbf{a'y} = \mathbf{k'(X'V^{-1}X)^{-}X'V^{-1}y} \] 如果我们设 \[ \boldsymbol{\beta^{0}} = \mathbf{(X'V^{-1}X)^{-}X'V^{-1}y} \] 即 \(\boldsymbol{\beta^{0}}\) 是下式的任意一个解，即 \(\boldsymbol{\beta^{0}}\) 是广义最小二乘方程组的一个解，其上标的 0 表示其是某一个解，而不是唯一解。 \[ \mathbf{(X'V^{-1}X)}\boldsymbol{\beta^{0}} = \mathbf{X'V^{-1}y} \] 因此，我们得到 \(\mathbf{k}^{\prime} \boldsymbol{\beta}\) 的 BLUE 值就是 \(\mathbf{k}^{\prime} \boldsymbol{\beta^{0}}\) 。

\(V^{-1}\) 的计算简化

在广义最小二乘估计时需要计算 \(\mathbf{V^{-1}}\) ，当表型值数据量很大时， \(\mathbf{V^{-1}}\) 计算困难甚至无法计算。Henderson 在 1959 年给出了 \(\mathbf{V^{-1}}\) 的一个计算公式，如下： \[ \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} \] 一般 \(\mathbf{R}^{-1}\) 很好计算，如果我们可以轻松计算 \(\mathbf{G^{-1}}\) 和 \(\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1}\) 的话 ，那么采用这个公式计算 \(\mathbf{V^{-1}}\) 就很有优势。在个体动物模型中貌似没法简单计算 \(\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1}\) ，所以这个计算公式现在好像没啥用了。

证明：我们只需要证明 \(\mathbf{(R+ZGZ')} \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{I}\) 即可，我们有 \[ \begin{aligned} &\mathbf{\left(R+Z G Z^{\prime}\right)\left[R^{-1}-R^{-1} Z\left(Z^{\prime} R^{-1} Z+G^{-1}\right)^{-1} Z^{\prime} R^{-1}\right]} \\ &=\mathbf{I+Z G Z^{\prime} R^{-1}-Z\left(Z^{\prime} R^{-1} Z+G^{-1}\right)^{-1} Z^{\prime} R^{-1}} \\ & \quad-\mathbf{Z G Z^{\prime} R^{-1} Z\left(Z^{\prime} R^{-1} Z+G^{-1}\right)^{-1} Z^{\prime} R^{-1}} \\ &=\mathbf{I+Z G Z^{\prime} R^{-1}-Z\left(I+G Z^{\prime} R^{-1} Z\right)\left(Z^{\prime} R^{-1} Z+G^{-1}\right)^{-1} Z^{\prime} R^{-1}} \\ &=\mathbf{I+Z G Z^{\prime} R^{-1}-Z G\left(G^{-1}+Z^{\prime} R^{-1} Z\right)\left(Z^{\prime} R^{-1} Z+G^{-1}\right)^{-1} Z^{\prime} R^{-1}} \\ &=\mathbf{ I+Z G Z^{\prime} R^{-1}-Z G Z^{\prime} R^{-1}} \\ &=\mathbf{I} \end{aligned} \] 得证。

其实这里可以推导出一个一般定理，如下：

对于矩阵 \(\mathbf{A}_{m \times m}, \mathbf{B}_{m \times n}, \mathbf{C}_{n \times n}, \mathbf{D}_{n \times m}\), 其中 \(\mathbf{A}, \mathbf{C}, \mathbf{C}^{-1}+\mathbf{D}\mathbf{A}^{-1}\mathbf{B}\) 可逆，则存在以下结果，证明过程同上。 \[ (\boldsymbol{A}+\boldsymbol{B C D})^{-1}=\boldsymbol{A}^{-1}-\boldsymbol{A}^{-1} \boldsymbol{B}\left(\boldsymbol{C}^{-1}+\boldsymbol{D} \boldsymbol{A}^{-1} \boldsymbol{B}\right)^{-1} \boldsymbol{D} \boldsymbol{A}^{-1} \]

BLUE的方差

假设我们有一个可估计的函数 \(\mathbf{K}^{\prime} \boldsymbol{\beta^{0}}\) ，那么其方差为 \[ \begin{aligned} \operatorname{Var}\left(\mathbf{K}^{\prime} \boldsymbol{\beta}^{o}\right) &=\operatorname{Var}\left[\mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{y}\right] \\ &=\mathbf{K}^{\prime}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{V} \mathbf{V}^{-1} \mathbf{X}\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{K} \quad\text { 注：这里应该需要假设 } \left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \text{是一个自反广义逆矩阵}\ \end{aligned} \] 由于 \(\mathbf{K}^{\prime} \boldsymbol{\beta^{0}}\) 可估计，因此一定可以找到一个矩阵 \(\mathbf{A}_{0}\) 使得 \(\mathbf{K = X'A_{0}}\) ，同上易证 $^{}({} ^{-1} )^{-} $ 具有不变性。

参考文献

1. Henderson，The Estimation of Environmental and Genetic Trends from Records Subject to Culling, 1959
2. Henderson，《Applications_of_Linear_Models_in_Animal_Breeding》, 1984
3. 张沅，张勤，《畜禽育种中的线性模型》