混合线性模型理论的第二篇，如何估计随机效应的值（假设不存在固定效应或者已知固定效应）。上一篇见混合线性模型一之广义最小二乘 。

代价函数

这里运行机器学习的思想，不同的预测函数的好坏需要通过代价函数来计算比较，或者说我们需要通过最小化代价函数来确定一个合适的预测函数，因此首要的问题是确定一个合适的代价函数，也就是说我们如何定义“最佳”。

本节内容主要参考 Henderson 1973年的论文。

最佳预测

假设 \(\mathbf{y}\) 是观测值向量，我们想要预测 \(\mathbf{w}\) 向量， \(\mathbf{w}\) 向量的第 \(i\) 个元素的预测值为 \[ \hat{w}_{i} = f(\mathbf{y}) \] 这里的 \(f\) 是某个关于观测值的函数，不一定要是线性函数。

此时我们的问题是，我们如何定义“最佳”的预测函数？ 一种方式是我们想尽可能地减少预测误差 (prediction error) ，即 \(\hat{w}_{i} - w_{i}\) ，那么我们如何去衡量预测误差的大小呢？从理论上，预测误差的大小可以用下式表示（等价于 \(L_{p}\) 范数） \[ E|\hat{w}_{i} - w_{i}|^{p} \] 即预测误差绝对值的 \(p\) 次方的期望值。那么 \(p\) 应该选多少呢？最普遍的选择是 2 ，因为 \(L_2\) 范数计算简便，并且有其它优良性质。

我们将最小化 \(E|\hat{w}_{i} - w_{i}|^{2}\) 的方法称为最小均方误差预测 (Minimum Mean Square Error Prediction) ，有时也称为最佳预测 (Best Prediction, BP) 。我们可以证明此时的估计值等价于下式 (Cochran , 1951 或 Rao 1965 )
\[ \hat{w}_{i} = E(w_{i} | \mathbf{y} ) \] 也就是说等同于给定 $ $ 的 \(w_{i}\) 的条件期望。

证明 (来自于 条件方差和预测 | 其重要属性与 5+ 示例) ：这里证明根据一个随机变量的观测值来预测另一个随机变量，假设我们有一个随机变量 X 的观测值，我们需要预测随机变量 Y ，证明最好的预测值为 \(E(Y|X)\) 。假设 \(g(X)\) 是任意一个预测函数，因此我们要证明 \[ E\left[(Y-g(X))^{2}\right] \geq E\left[(Y-E[Y \mid X])^{2}\right] \] 我们有 \[ \begin{aligned} E\left[(Y-g(X))^{2} \mid X\right]=& E\left[(Y-E[Y \mid X]+E[Y \mid X]-g(X))^{2} \mid X\right] \\ =& E\left[(Y-E[Y \mid X])^{2} \mid X\right] \\ &+E\left[(E[Y \mid X]-g(X))^{2} \mid X\right] \\ &+2 E[(Y-E[Y \mid X])(E[Y \mid X]-g(X)) \mid X] \end{aligned} \] 在给定 X 的前提下，\(\mathrm{E}[\mathrm{Y} \mid \mathrm{X}]-\mathrm{g}(\mathrm{X})\) 是一个常数，因此上式的最后一项为 0 ，证明如下 \[ \begin{aligned} &E[(Y-E[Y \mid X])(E[Y \mid X]-g(X)) \mid X] \\ &=(E[Y \mid X]-g(X)) E[Y-E[Y \mid X] \mid X] \\ &=(E[Y \mid X]-g(X))(E[Y \mid X]-E[Y \mid X]) \\ &=0 \end{aligned} \] 因此下式成立 \[ E\left[(Y-g(X))^{2} \mid X\right]= E\left[(Y-E[Y \mid X])^{2} \mid X\right] +E\left[(E[Y \mid X]-g(X))^{2} \mid X\right] \\ \] 其中右手第二项是一个非负数，因此下式成立 \[ E\left[(Y-g(X))^{2}\right] \geq E\left[(Y-E[Y \mid X])^{2}\right. \] 得证。

因此我们必须先已知 \([\mathbf{y}, \mathbf{w}]\) 的联合分布，然后再推导条件分布。

BP 的一些性质如下：

1. \(E(\hat{w}_{i}) = E(w_{i})\)
2. \(Var(\hat{\mathbf{w}}- \mathbf{w} ) = \text{ average conditional } Var(\mathbf{w | y})\) 。这个没看懂，感觉怪怪的。
3. 在所有可能的预测值中，\(\hat{w}_{i} = E(w_{i} | \mathbf{y} )\) 最大化了 \(\hat{w}_{i}\) 和 \(w_{i}\) 的相关

最佳线性预测

通常在实践中，我们不知道 \([\mathbf{y}, \mathbf{w}]\) 的联合分布 ，我们就也没有办法计算条件分布。那么此时我们怎么办呢？我们可以将预测值的范围限制为线性预测值，这样我们可以极大地简化问题。我们称这种方法为 最佳线性预测 (Best Linear Prediction, BLP)，也就是育种值估计方法中的选择指数法 (Selection Index) 。

我们先看一般 最佳线性预测的一般推导，假设我们的观测值向量为 \(\mathbf{y}\) ，我们需要预测的随机变量是 \(w\) （注意这里我们采用单个变量），我们采用的线性预测函数为 \(\hat{w} = \mathbf{a'y} + b\) ，这里 \(\mathbf{a'}\) 是一个向量，\(b\) 是一个标量。我们的目标是找到一组 \(\mathbf{a'}\) 和 \(b\)，使得 \(E(\hat{w} - w)^{2}\) 最小。

均值，协方差等设定如下 \[ \begin{aligned} E(w) &=\gamma, \\ E(\mathbf{y}) &= \boldsymbol{\alpha}, \\ \operatorname{Cov}(\mathbf{y}, w) &=\mathbf{c}, \text { and } \\ \operatorname{Var}(\mathbf{y}) &=\mathbf{V} \end{aligned} \] 那么 $$ \[\begin{aligned} E\left(\mathbf{a}^{\prime} \mathbf{y}+b-w\right)^{2} =&E\left((\mathbf{a}^{\prime} \mathbf{y})^{2}+b^{2}+w^{2}+2(\mathbf{a}^{\prime} \mathbf{y})b-2(\mathbf{a}^{\prime} \mathbf{y})w-2bw\right)\\ =&E((\mathbf{a}^{\prime} \mathbf{y})^{2})+b^{2}+E(w^{2})+2E((\mathbf{a}^{\prime} \mathbf{y}))b-2E((\mathbf{a}^{\prime} \mathbf{y})w)-2bE(w)\\ =& \mathbf{a}^{\prime} \mathbf{V a} + (\mathbf{a}^{\prime} \boldsymbol{\alpha})^2 + b^{2} + \operatorname{Var}(w)+\gamma^{2}+2 \mathbf{a}^{\prime} \boldsymbol{\alpha} b - 2 \mathbf{a}^{\prime} (\mathbf{c}+ \boldsymbol{\alpha} \gamma) -2b \gamma\\ =& \mathbf{a}^{\prime} \mathbf{V a} + \mathbf{a}^{\prime} \boldsymbol{\alpha} \boldsymbol{\alpha}^{\prime} \mathbf{a} + b^{2} + \operatorname{Var}(w)+\gamma^{2}+2 b \mathbf{a}^{\prime} \boldsymbol{\alpha} - 2 \mathbf{a}^{\prime} \mathbf{c}- 2 \gamma \mathbf{a}^{\prime} \boldsymbol{\alpha} -2b \gamma\\ \end{aligned}\] \[ 分别对 $\mathbf{a'}$ 和 $b$ 求偏导，使之为 0，得到 \] \[\begin{aligned} &\partial L_{2} / \partial \mathbf{a}= 2\mathbf{V} \mathbf{a}+ 2\boldsymbol{\alpha \alpha^{\prime}} \mathbf{a}+ 2b \boldsymbol{\alpha} - 2 \mathbf{c} - 2 \gamma \boldsymbol{\alpha} = \mathbf{0} \\ &\partial L_{2} / \partial b=2b+2\mathbf{a}^{\prime} \boldsymbol{\alpha} - 2 \gamma= 0 \end{aligned}\] \[ 整理得到下面的方程组 \] ( \[\begin{array}{cc} \mathbf{V}+\boldsymbol{\alpha \alpha^{\prime}} & \boldsymbol{\alpha} \\ \boldsymbol{\alpha^{\prime}} & 1 \end{array}\] )( \[\begin{array}{l} \mathbf{a} \\ b \end{array}\] )=( \[\begin{array}{c} \mathbf{c}+\boldsymbol{\alpha} \gamma \\ \gamma \end{array}\] ) \[ 根据第二个方程，我们得到 $b = \gamma - \mathbf{a}^{\prime} \boldsymbol{\alpha} = \gamma - \boldsymbol{\alpha}^{\prime} \mathbf{a}$ ，带入到第一个方程得到 \] \[\begin{aligned} (\mathbf{V}+\boldsymbol{\alpha \alpha^{\prime}}) \mathbf{a} + (\gamma - \boldsymbol{\alpha}^{\prime} \mathbf{a}) \boldsymbol{\alpha} - \mathbf{c} - \gamma \boldsymbol{\alpha} &= \mathbf{0} \\ (\mathbf{V}+\boldsymbol{\alpha \alpha^{\prime}}) \mathbf{a} - (\boldsymbol{\alpha}^{\prime} \mathbf{a})\boldsymbol{\alpha} - \mathbf{c} &= \mathbf{0} \\ (\mathbf{V}+\boldsymbol{\alpha \alpha^{\prime}}) \mathbf{a} - \boldsymbol{\alpha} (\boldsymbol{\alpha}^{\prime} \mathbf{a}) - \mathbf{c} &= \mathbf{0} \\ \mathbf{Va} &= \mathbf{c} \\ \end{aligned}\]

\[ 因此我们得到解为 \] = , b = -^{} ^{-1} \[ 因此，预测值为 \] =+^{} ^{-1}(-) \[ 从单个元素推广至向量，得到 \] =+^{} ^{-1}(-) $$ 其中，\(\boldsymbol{\gamma} = E(\mathbf{w})\) ， \(\mathbf{C} = Cov(\mathbf{y}, \mathbf{w}^{\prime})\) 。

当 \(\mathbf{y}, w\) 的共同分布为多元正态分布时，BLP 就是 \(E(w \mid y)\) ，即此时 BLP = BP 。

均值及协方差

首先我们看期望，BLP是一个无偏估计值，即 \[ E(\hat{w})=E(w) \] 证明： \[ \begin{aligned} E(\hat{w}) &=E\left[\gamma+\mathbf{c}^{\prime} \mathbf{V}^{-1}(\mathbf{y}-\boldsymbol{\alpha})\right] \\ &=\gamma+\mathbf{c}^{\prime} \mathbf{V}^{-1}(\boldsymbol{\alpha}-\boldsymbol{\alpha}) \\ &=\gamma=E(w) . \end{aligned} \] 其协方差矩阵为 \[ \operatorname{Var}(\hat{w})=\operatorname{Var}\left(\mathbf{c}^{\prime} \mathbf{V}^{-1} \mathbf{y}\right)=\mathbf{c}^{\prime} \mathbf{V}^{-1} \mathbf{V} \mathbf{V}^{-1} \mathbf{c}=\mathbf{c}^{\prime} \mathbf{V}^{-1} \mathbf{c} \] 其与真实值的协方差矩阵为 \[ \operatorname{Cov}(\hat{w}, w)=\mathbf{c}^{\prime} \mathbf{V}^{-1} \operatorname{Cov}(\mathbf{y}, w)=\mathbf{c}^{\prime} \mathbf{V}^{-1} \mathbf{c}=\operatorname{Var}(\hat{w}) \] 误差方差为 \[ \begin{aligned} \operatorname{Var}(\hat{w}-w)&=\operatorname{Var}(w)+\operatorname{Var}(\hat{w})-2\operatorname{Cov}(\hat{w}, w)\\ &=\operatorname{Var}(w)+\operatorname{Var}(\hat{w})-2\operatorname{Var}(\hat{w})\\ &=\operatorname{Var}(w)-\operatorname{Var}(\hat{w}) \end{aligned} \] 下面我们证明 BLP 最大化了 \(r_{\hat{w} w}\) ，根据相关系数公式，我们有一般公式 $$ \[\begin{aligned} r_{\hat{w} w} &= \frac{Cov(\hat{w} ,w)}{\left(Var(\hat{w} ) Var(w)\right)^{0.5}} \\ &= \frac{Cov(\mathbf{a'y} + b ,w)}{\left(Var(\mathbf{a'y} + b) Var(w)\right)^{0.5}} \\ &= \frac{\mathbf{a'c}}{\left(\mathbf{a'Va} Var(w)\right)^{0.5}} \\ \end{aligned}\]

\[ 取对数（容易看出 $r_{\hat{w} w} > 0$ ），得到 \] r=^{} -0.5 (w) \[ 对 $\mathbf{a}$ 微分，使之为 0 ，得到 \] r / = - = \[ 我们得到 \] = \[ 注意到两边的分母均为标量，因此 $\mathbf{Va}$ 和 $\mathbf{c}$ 这两个向量平行，即满足下式，其中 $k$ 是一个标量 \] = k \[ 因此，我们得到 $\mathbf{a} = k \mathbf{V^{-1}c}$ ，带入到上式中，得到 \] = $$ 因此 \(k\) 可以取任意值，不妨设 \(k=1\) ，此时 \(\mathbf{a} = \mathbf{V^{-1}c}\) ，即为 BLP 估计值，因此我们证明了 BLP 最大化了 \(r_{\hat{w} w}\) 。

选择指数法

下面的内容来自于张沅老师的 《畜禽育种中的线性模型》，选择指数法就是BLP在育种中的应用。

选择指数法需要满足三个条件或者假设：

1. 不存在系统环境效应，或者所有系统环境效应均已知
2. 被选择的个体间不存在固定的遗传差异
3. 误差和育种值的方差，协方差矩阵均已知

其实前两个条件可以合并为不存在固定效应或者已知固定效应，也就是说此时的模型对应为随机模型。

根据定义，选择指数的公式为 \[ I = \sum b_{i}(y_{i} - \mu_{i}) \] 其中：

\(I\) = 某个体一个性状的选择指数

\(b_{i}\) = 与该个体有关的第 \(i\) 个信息的加权系数（偏回归系数）

\(y_{i}\) = 第 \(i\) 个信息来源的表型值

\(u_{i}\) = 与该表型值对应的总体均数

上式可以用矩阵形式表示： \[ I=\mathbf{b'd} \] 其中：

\(\mathbf{b}\) = 加权系数向量

$ $ ，即各信息来源的表型值与各自的群体平均数离差向量

我们确定了选择指数的公式之后，我们需要找一组合适的加权系数，设该个体的真实育种值为 \(A\) ，则代价函数为 \(E(I-A)^{2}\) ，根据协方差公式，我们有 $ Var(I-A) = E(I-A)^{2} - (E(I-A))^{2} = E(I-A)^{2}$ ，因此等价于需要一个 \(\mathbf{b}\) 向量使得误差方差 $ Var(I-A) $ 最小，我们有 \[ \begin{aligned} \operatorname{Var}(I-A) &=\operatorname{Var}(I)+\operatorname{Var}(A)-2 \operatorname{Cov}(I, A) \\ &=\operatorname{Var}\left(\mathbf{b}^{\prime} \mathbf{d}\right)+\sigma_{A}^{2}-2 \operatorname{Cov}\left(\mathbf{b}^{\prime} \mathbf{d}, A\right) \\ &=\mathbf{b}^{\prime} \operatorname{\mathbf{Vb}}+\sigma_{A}^{2}-2 \mathbf{b}^{\prime} \mathbf{c} \end{aligned} \] 其中 \[ \mathbf{V}=\operatorname{Var}(\mathbf{d})=\operatorname{Var}(\mathbf{y}-\boldsymbol{\mu})=\operatorname{Var}(\mathbf{y}) \]

\[ \mathbf{c}=\operatorname{Cov}(\mathbf{d}, \mathrm{A})=\operatorname{Cov}(\mathbf{y}, \mathrm{A}) \]

求上式关于 \(\mathbf{b}\) 的最小值，我们求偏导数，使之为 \(\mathbf{0}\) 得到 \[ \frac{\partial}{\partial b} \operatorname{Var}(I-A)=2 \mathbf{Vb}-2 \mathbf{ c}=\mathbf{0} \] 我们得到 \(\mathbf{Vb = c}\) ，因此我们有 \[ \mathbf{b = V^{-1} c} \] 如果我们根据单一表型来预测育种值，即此时 \(\mathbf{y}\) 是一个标量，那么有 \[ b = \mathbf{ V^{-1} c } = (Cov(y))^{-1} Cov(y,A) = \frac{Cov(y,A)}{Cov(y)} \] 根据数量遗传学理论，我们知道两个个体的加性遗传相关等于两个个体的育种值之间的相关，即 \[ a_{XY} = \frac{Cov(A_{X},A_{Y})}{\sigma_{X}\sigma_{Y}} \] 这里由于所有个体属于同一群体，因此我们有 \(\sigma_{X} = \sigma_{Y} = \sigma_{A}\) ，因此 \(a_{XY} = \frac{Cov(A_{X},A_{Y})}{\sigma_{A}^{2}}\) 。假设我们的目标个体称为 \(X\) ，因此\(\mathbf{c}\) 向量中的元素 $c_{i} =(y_{i}, ) =(, ) =(, ) = a_{iX} _{A}^{2} $ ，即等于相应的加性遗传相关乘以加性遗传方差。设 \(\mathbf{a}\) 为个体与各信息来源的个体间加性遗传相关的向量，因此我们有 \(\mathbf{c = a \sigma_{A}^{2}}\) ，因此 \(\mathbf{b}\) 的计算公式可以写为 \[ \mathbf{b = V^{-1} a \sigma_{A}^{2}} \] 因此，选择指数的基本公式可以写为： \[ I = \mathbf{c'V^{-1}(y-\boldsymbol{\mu})} = \sigma_{A}^{2} \mathbf{a'V^{-1}(y-\boldsymbol{\mu})} \]

从单个元素推广至向量，得到 \[ \mathbf{I}=\mathbf{C}^{\prime} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta})=\mathbf{GZ'} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}) \] 其中，\(\mathbf{C}^{\prime} = Cov(\boldsymbol{\mu}, \mathbf{y}^{\prime}) = \mathbf{GZ'}\)

选择指数估计的准确性的高低，直接取决于信息来源的种类和数量。而准确度的高低又直接关系到将获得的遗传进展的大小，所以这是一个非常重要的指标。准确性通常用真实育种值与估计育种值的相关系数来表示，即： \[ r_{AI} = \frac{Cov(A,I)}{\sigma_{A} \sigma_{I}} \] 因为 \(\sigma_{I}^{2}=\mathbf{b}^{\prime} \mathbf{V b} = \mathbf{b}^{\prime} \mathbf{c},\quad \operatorname{Cov}(A, I)=\mathbf{b}^{\prime} \mathbf{c}=\sigma_{I}^{2}\) ，所以 \[ r_{AI} = \frac{ \sigma_{I}^{2}}{\sigma_{A} \sigma_{I}} = \frac{ \sigma_{I}}{\sigma_{A} } = \frac{\mathbf{(b'c)}^{0.5}}{\sigma_{A} } = \frac{(\mathbf{b'a }\sigma^{2}_{A})^{0.5}}{\sigma_{A} } = (\mathbf{b'a })^{0.5} \] BLP 的一些性质为

1. 在所有关于 \(\mathbf{y}\) 的线性组合中，BLP 最大化了 \(I\) 和 \(A\) 的相关
2. BLP 具有不变性，即 \(\text{BLP of }\mathbf{m^{\prime}w} = \mathbf{m^{\prime}\hat{w}}\)
3. 当 \([\mathbf{y}, \mathbf{w}]\) 服从多变量正态分布是，此时的 BLP 就是 BP

缺点

最主要的缺点就是需要已知 \(\mathbf{}\) \(\mathbf{y}\) 的均值，即 \(\mathbf{X} \boldsymbol{\beta}\) ，这在实际情况中基本是不可能实现的。

参考文献

1. Henderson，1973，SIRE EVALUATION AND GENETIC TRENDS
2. Henderson，《Applications_of_Linear_Models_in_Animal_Breeding》, 1984
3. 张沅，张勤，《畜禽育种中的线性模型》