这是方差组分估计方法的第四篇博客，介绍 EM-REML 算法。

混合线性模型的一般式子

我们将混合线性模型写为 \[ \mathbf{y} = \mathbf{X}\boldsymbol{\beta}+\mathbf{Zu}+\mathbf{e} \] 其中 \(\mathbf{u}\) 和 \(\mathbf{e}\) 的联合分布为 \[ \left[\begin{array}{l} \mathbf{u} \\ \mathbf{e} \end{array}\right] \sim \mathrm{N}\left(\left[\begin{array}{l} \mathbf{0} \\ \mathbf{0} \end{array}\right], \theta\left[\begin{array}{cc} \mathbf{G} & \mathbf{0} \\ \mathbf{0} & \mathbf{R} \end{array}\right]\right) \] 其中 \(\mathbf{G} = \mathbf{G}(\boldsymbol{\gamma})\) ， \(\mathbf{R} = \sigma^{2} \mathbf{\Sigma}\) ，\(\mathbf{\Sigma} = \mathbf{\Sigma}(\boldsymbol{\phi})\) 。这里我们假定这三个矩阵均为正定矩阵。

这里的 \(\boldsymbol{\gamma}\) 和 \(\boldsymbol{\phi}\) 是与随机效应 \(\mathbf{u}\) 和 \(\mathbf{e}\) 相关的方差组分的向量，而标量参数 \(\sigma^{2}\) 和 \(\theta\) 是与残差和全部方差组分的缩放系数。因此我们需要估计的方差参数为 \(\theta\) 和 \(\boldsymbol{\kappa} = (\boldsymbol{\gamma}^{\prime}, \sigma^{2}, \boldsymbol{\phi}^{\prime} )^{\prime}\) 。

因为 \(\mathbf{y}\) 的分布为 \[ \mathbf{y} \sim \mathrm{N}(\mathbf{X}\boldsymbol{\beta}, \theta\mathbf{H}) \] 其中 \(\mathbf{H} = \mathbf{ZGZ^{\prime}}+\mathbf{R}\) 是一个正定矩阵。

这里的标量参数 \(\sigma^{2}\) 和 \(\theta\) 均是用于缩放的参数，我们发现如果我们同时包含这两个参数，那么这两个参数的值便均无法确定，因为此时 \(\mathrm{Var}(\mathbf{u}) = \theta \mathbf{G}(\boldsymbol{\gamma})\) ，\(\mathrm{Var}(\mathbf{e}) = \theta \sigma^{2} \mathbf{\Sigma}(\boldsymbol{\phi})\) 。实际上我们至少会将其中一个参数设定为 1 。

按照育种模型的一般思路，这里我将 \(\sigma^{2}\) 和 \(\theta\) 均设为 1， 此时 \(\mathrm{Var}(\mathbf{y}) = \mathbf{V} = \mathbf{ZGZ^{\prime}}+\mathbf{R} = \mathbf{ZG(\boldsymbol{\gamma})Z^{\prime}}+\mathbf{\Sigma}(\boldsymbol{\phi})\) ，我们需要估计的参数为 \(\boldsymbol{\kappa} = (\boldsymbol{\gamma}^{\prime}, \boldsymbol{\phi}^{\prime} )^{\prime}\) 。

假设我们有 \(s\) 个随机效应，即 \(\mathbf{u} = (\mathbf{u}_{1}^{\prime}, \mathbf{u}_{2}^{\prime}, \cdots, \mathbf{u}_{s}^{\prime})^{\prime}\) ，因此此时随机效应的协方差矩阵可以写为 \[ \mathbf{G} = \oplus_{i=1}^{s}\mathbf{G}_{i} = \left[\begin{array}{ccccc} \mathbf{G}_1 & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{G}_2 & \cdots & \mathbf{0} & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \mathbf{0} & \mathbf{0} & \ldots & \mathbf{G}_{s-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{0} & \mathbf{G}_s \end{array}\right] \] 此时的方差组分 \(\boldsymbol{\gamma}\) 也可以进行相应的分块，例如 \(\mathbf{G}_{i} = \mathbf{G}_{i}(\boldsymbol{\gamma}_{i})\) ，因此我们有 \(\boldsymbol{\gamma} = (\boldsymbol{\gamma}_{1}^{\prime}, \boldsymbol{\gamma}_{2}^{\prime}, \cdots, \boldsymbol{\gamma}_{s}^{\prime})^{\prime}\) 。

但是对于残差，我们通常还是认为所有数据来自与一个数据集，即残差是均质的，其协方差矩阵不用分块。

推导 C 逆矩阵

假设 \(\mathbf{X}\) 矩阵满秩，则混合线性模型的系数矩阵 \(\mathbf{C}\) 可逆。

我们将原始的系数矩阵分块为 \[ \begin{aligned} \mathbf{C} &= \left[\begin{array}{c} \mathbf{C}_{11} & \mathbf{C}_{12} \\ \mathbf{C}_{21} & \mathbf{C}_{22} \\ \end{array}\right] \\ &= \left[\begin{array}{c} \mathbf{X^{\prime}R^{-1}X} & \mathbf{X^{\prime}R^{-1}Z} \\ \mathbf{Z^{\prime}R^{-1}X} & \mathbf{Z^{\prime}R^{-1}Z+G^{-1}} \\ \end{array}\right] \\ \end{aligned} \] 根据舒尔补公式，系数矩阵逆矩阵可以写为 \[ \begin{aligned} \mathbf{C}^{-1} &= \left[\begin{array}{c} \mathbf{C}^{11} & \mathbf{C}^{12} \\ \mathbf{C}^{21} & \mathbf{C}^{22} \\ \end{array}\right] \\ &= \left[\begin{array}{c} (\mathbf{C}_{11} - \mathbf{C}_{12}\mathbf{C}_{22}^{-1}\mathbf{C}_{21})^{-1} & -(\mathbf{C}_{11} - \mathbf{C}_{12}\mathbf{C}_{22}^{-1}\mathbf{C}_{21})^{-1}\mathbf{C}_{12}\mathbf{C}_{22}^{-1} \\ -\mathbf{C}_{22}^{-1}\mathbf{C}_{21}(\mathbf{C}_{11} - \mathbf{C}_{12}\mathbf{C}_{22}^{-1}\mathbf{C}_{21})^{-1} & \mathbf{C}_{22}^{-1}+\mathbf{C}_{22}^{-1}\mathbf{C}_{21}(\mathbf{C}_{11} - \mathbf{C}_{12}\mathbf{C}_{22}^{-1}\mathbf{C}_{21})^{-1}\mathbf{C}_{12}\mathbf{C}_{22}^{-1} \\ \end{array}\right] \\ &= \left[\begin{array}{c} \mathbf{C}^{11} & -\mathbf{C}^{11}\mathbf{C}_{12}\mathbf{C}_{22}^{-1} \\ -\mathbf{C}_{22}^{-1}\mathbf{C}_{21}\mathbf{C}^{11} & \mathbf{C}_{22}^{-1}+\mathbf{C}_{22}^{-1}\mathbf{C}_{21}\mathbf{C}^{11}\mathbf{C}_{12}\mathbf{C}_{22}^{-1} \\ \end{array}\right] \\ \end{aligned} \] 在之前的推导中，我们已经得到 \[ \mathbf{C^{11}} = (\mathbf{C}_{11} - \mathbf{C}_{12}\mathbf{C}_{22}^{-1}\mathbf{C}_{21})^{-1} = (\mathbf{X^{\prime}V^{-1}X})^{-1} \] 根据附录1，我们有 \(\mathbf{(Z^{\prime}R^{-1}Z+G^{-1})^{-1}Z^{\prime}R^{-1}} = \mathbf{GZ^{\prime}(R+ZGZ^{\prime})^{-1}}= \mathbf{GZ^{\prime}V^{-1}}\) ，因此其它分块子矩阵推导如下 \[ \begin{aligned} \mathbf{C}^{12} &= -\mathbf{C}^{11}\mathbf{C}_{12}\mathbf{C}_{22}^{-1} \\ &= -(\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X^{\prime}R^{-1}Z}(\mathbf{Z^{\prime}R^{-1}Z+G^{-1}})^{-1} \\ &= -(\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime}\mathbf{V^{-1}ZG} \\ \mathbf{C}^{21} &= (\mathbf{C}^{12})^{\prime} \\ &= -\mathbf{GZ^{\prime}V^{-1}} \mathbf{X} (\mathbf{X^{\prime}V^{-1}X})^{-1} \\ \mathbf{C}^{22} &= \mathbf{C}_{22}^{-1}+\mathbf{C}_{22}^{-1}\mathbf{C}_{21}\mathbf{C}^{11}\mathbf{C}_{12}\mathbf{C}_{22}^{-1} \\ &= (\mathbf{Z^{\prime}R^{-1}Z+G^{-1}})^{-1}+(\mathbf{Z^{\prime}R^{-1}Z+G^{-1}})^{-1}\mathbf{Z^{\prime}R^{-1}X}(\mathbf{X^{\prime}V^{-1}X})^{-1}\mathbf{X^{\prime}R^{-1}Z} (\mathbf{Z^{\prime}R^{-1}Z+G^{-1}})^{-1} \\ &=(\mathbf{Z^{\prime}R^{-1}Z+G^{-1}})^{-1}+ \mathbf{GZ^{\prime}V^{-1}}\mathbf{X}(\mathbf{X^{\prime}V^{-1}X})^{-1}\mathbf{X^{\prime}} \mathbf{V^{-1}ZG} \\ &=(\mathbf{Z^{\prime}R^{-1}Z+G^{-1}})^{-1}+ \mathbf{GZ^{\prime}}\mathbf{(V^{-1}-P)} \mathbf{ZG} \quad \because \mathbf{P}=\mathbf{V}^{-1}-\mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \\ &=(\mathbf{Z^{\prime}R^{-1}Z+G^{-1}})^{-1}+ \mathbf{GZ^{\prime}}\mathbf{V^{-1}} \mathbf{ZG}- \mathbf{GZ^{\prime}}\mathbf{P} \mathbf{ZG} \\ &=(\mathbf{Z^{\prime}R^{-1}Z+G^{-1}})^{-1}+ \mathbf{(Z^{\prime}R^{-1}Z+G^{-1})^{-1}Z^{\prime}R^{-1}} \mathbf{ZG}- \mathbf{GZ^{\prime}}\mathbf{P} \mathbf{ZG} \\ &=(\mathbf{Z^{\prime}R^{-1}Z+G^{-1}})^{-1}(\mathbf{I}+\mathbf{Z^{\prime}R^{-1}} \mathbf{ZG})- \mathbf{GZ^{\prime}}\mathbf{P} \mathbf{ZG} \\ &=(\mathbf{Z^{\prime}R^{-1}Z+G^{-1}})^{-1}(\mathbf{G}^{-1}+\mathbf{Z^{\prime}R^{-1}} \mathbf{Z})\mathbf{G}- \mathbf{GZ^{\prime}}\mathbf{P} \mathbf{ZG} \\ &=\mathbf{G}- \mathbf{GZ^{\prime}}\mathbf{P} \mathbf{ZG} \\ \end{aligned} \] 因此我们得到 $$ \[\begin{aligned} \mathbf{C}^{-1} &= \left[\begin{array}{c} \mathbf{C}^{11} & \mathbf{C}^{12} \\ \mathbf{C}^{21} & \mathbf{C}^{22} \\ \end{array}\right] \\ &= \left[\begin{array}{c} (\mathbf{X^{\prime}V^{-1}X})^{-1} & -(\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime}\mathbf{V^{-1}ZG} \\ -\mathbf{GZ^{\prime}V^{-1}} \mathbf{X} (\mathbf{X^{\prime}V^{-1}X})^{-1} & \mathbf{G}- \mathbf{GZ^{\prime}}\mathbf{P} \mathbf{ZG} \\ \end{array}\right] \\ \end{aligned}\]

$$

推导 P 矩阵

之前我们得到了两个 P 矩阵的式子，如下（\(\mathbf{K}\) 为满足 \(\mathbf{K'X}=0\) 的一个矩阵） \[ \mathbf{P}=\mathbf{V}^{-1}-\mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1}=\mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K}^{\prime} \] 这里我们推导第三个式子 \[ \mathbf{P} = \mathbf{R^{-1}-R^{-1}WC^{-1}W^{\prime}R^{-1}} \] 其中 \(\mathbf{W} = \left[\begin{array}{l} \mathbf{X} & \mathbf{Z} \end{array}\right]\) ，\(\mathbf{C}\) 为系数矩阵。

证明过程

根据下面的附录1，我们有 \[ \mathbf{R^{-1}Z(Z^{\prime}R^{-1}Z+G^{-1})^{-1}} = \mathbf{(R+ZGZ^{\prime})^{-1}ZG} \] 我们可以将其写为 \[ \mathbf{R^{-1}ZM} = \mathbf{V^{-1}ZG} \] 其中 \(\mathbf{M}=\mathbf{(Z^{\prime}R^{-1}Z+G^{-1})^{-1}}\) 。

我们也可以用 \(\mathbf{M}\) 来表示 \(\mathbf{V}^{-1}\) \[ \begin{aligned} \mathbf{V^{-1}} &= \mathbf{R^{-1}-R^{-1}Z(Z^{\prime}R^{-1}Z+G^{-1})^{-1}Z^{\prime}R^{-1}} \\ &=\mathbf{R^{-1}-R^{-1}ZMZ^{\prime}R^{-1}} \end{aligned} \] 我们同样可以用 \(\mathbf{M}\) 来表示 \(\mathbf{C}^{-1}\) 矩阵的元素 $$ \[\begin{aligned} \mathbf{C}^{12} &= -(\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime}\mathbf{V^{-1}ZG} \\ &= -(\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime}\mathbf{R^{-1}ZM} \\ \mathbf{C}^{21} &= -\mathbf{MZ^{\prime}R^{-1}}\mathbf{X}(\mathbf{X^{\prime}V^{-1}X})^{-1} \\ \mathbf{C}^{22} &=(\mathbf{Z^{\prime}R^{-1}Z+G^{-1}})^{-1}+ \mathbf{GZ^{\prime}V^{-1}}\mathbf{X}(\mathbf{X^{\prime}V^{-1}X})^{-1}\mathbf{X^{\prime}} \mathbf{V^{-1}ZG} \\ &=\mathbf{M}+ \mathbf{MZ^{\prime}R^{-1}}\mathbf{X}(\mathbf{X^{\prime}V^{-1}X})^{-1}\mathbf{X^{\prime}} \mathbf{R^{-1}ZM} \\ \end{aligned}\] \[ 因此系数矩阵逆矩阵可以表示为 \] \[\begin{aligned} \mathbf{C}^{-1} &= \left[\begin{array}{c} (\mathbf{X^{\prime}V^{-1}X})^{-1} & -(\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime}\mathbf{R^{-1}ZM} \\ -\mathbf{MZ^{\prime}R^{-1}}\mathbf{X}(\mathbf{X^{\prime}V^{-1}X})^{-1} & \mathbf{M}+ \mathbf{MZ^{\prime}R^{-1}}\mathbf{X}(\mathbf{X^{\prime}V^{-1}X})^{-1}\mathbf{X^{\prime}} \mathbf{R^{-1}ZM} \\ \end{array}\right] \\ \end{aligned}\] \[ 因此，我们有 \] \[\begin{aligned} \mathbf{WC^{-1}W^{\prime} } &= \left[\begin{array}{l} \mathbf{X} & \mathbf{Z} \end{array}\right] \left[\begin{array}{c} (\mathbf{X^{\prime}V^{-1}X})^{-1} & -(\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime}\mathbf{R^{-1}ZM} \\ -\mathbf{MZ^{\prime}R^{-1}}\mathbf{X}(\mathbf{X^{\prime}V^{-1}X})^{-1} & \mathbf{M}+ \mathbf{MZ^{\prime}R^{-1}}\mathbf{X}(\mathbf{X^{\prime}V^{-1}X})^{-1}\mathbf{X^{\prime}} \mathbf{R^{-1}ZM} \\ \end{array}\right] \left[\begin{array}{l} \mathbf{X}^{\prime} \\ \mathbf{Z}^{\prime} \end{array}\right] \\ &= \mathbf{X} (\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime} - \mathbf{X} (\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime}\mathbf{R^{-1}ZM}\mathbf{Z}^{\prime} - \mathbf{Z} \mathbf{MZ^{\prime}R^{-1}}\mathbf{X}(\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime} \\ & \quad + \mathbf{Z}\mathbf{M} \mathbf{Z}^{\prime} + \mathbf{Z} \mathbf{MZ^{\prime}R^{-1}}\mathbf{X}(\mathbf{X^{\prime}V^{-1}X})^{-1}\mathbf{X^{\prime}} \mathbf{R^{-1}ZM} \mathbf{Z}^{\prime} \end{aligned}\] \[ 因此 \] \[\begin{aligned} &\mathbf{R^{-1}-R^{-1}WC^{-1}W^{\prime}R^{-1}} \\ &= \mathbf{R^{-1}} -\mathbf{R^{-1}} \mathbf{X} (\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime} \mathbf{R^{-1}} + \mathbf{R^{-1}} \mathbf{X} (\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime}\mathbf{R^{-1}ZM}\mathbf{Z}^{\prime} \mathbf{R^{-1}} + \mathbf{R^{-1}} \mathbf{Z} \mathbf{MZ^{\prime}R^{-1}}\mathbf{X}(\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime} \mathbf{R^{-1}} \\ & \quad - \mathbf{R^{-1}} \mathbf{Z}\mathbf{M} \mathbf{Z}^{\prime} \mathbf{R^{-1}} - \mathbf{R^{-1}} \mathbf{Z} \mathbf{MZ^{\prime}R^{-1}}\mathbf{X}(\mathbf{X^{\prime}V^{-1}X})^{-1}\mathbf{X^{\prime}} \mathbf{R^{-1}ZM} \mathbf{Z}^{\prime} \mathbf{R^{-1}} \\ &= \mathbf{R^{-1}} - \mathbf{R^{-1}} \mathbf{Z}\mathbf{M} \mathbf{Z}^{\prime} \mathbf{R^{-1}} - (\mathbf{R^{-1}}-\mathbf{R^{-1}} \mathbf{Z} \mathbf{MZ^{\prime}R^{-1}} ) \mathbf{X} (\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime} \mathbf{R^{-1}}\\ &+ (\mathbf{R^{-1}}-\mathbf{R^{-1}} \mathbf{Z} \mathbf{MZ^{\prime}R^{-1}}) \mathbf{X} (\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime}\mathbf{R^{-1}ZM}\mathbf{Z}^{\prime} \mathbf{R^{-1}} \\ &= \mathbf{R^{-1}} - \mathbf{R^{-1}} \mathbf{Z}\mathbf{M} \mathbf{Z}^{\prime} \mathbf{R^{-1}} \\ & - (\mathbf{R^{-1}}-\mathbf{R^{-1}} \mathbf{Z} \mathbf{MZ^{\prime}R^{-1}} ) \mathbf{X} (\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime} (\mathbf{R^{-1}} -\mathbf{R^{-1}ZM}\mathbf{Z}^{\prime} \mathbf{R^{-1}} )\\ &= \mathbf{V}^{-1} - \mathbf{V}^{-1} \mathbf{X} (\mathbf{X^{\prime}V^{-1}X})^{-1} \mathbf{X}^{\prime} \mathbf{V}^{-1} \\ &= \mathbf{P} \end{aligned}\]

$$ 证明完毕

EM-REML

在之前的DF-REML 推导中，我们得到 REML 的一个似然函数的式子如下： \[ l(\mathbf{y}_{2}) = -\frac{1}{2}\left( \ln |\mathbf{V}| + \ln |\mathbf{X^{\prime}V^{-1}X}|+\mathbf{y}^{\prime} \mathbf{P} \mathbf{y} \right) \] 对于某个参数 \(\boldsymbol{\kappa}_{i}\) 进行求导，如下： \[ \begin{aligned} U(\kappa_{i}) &= \frac{\partial{l_{R}(\mathbf{y}_{2})}}{\partial{\kappa_{i}}} \\ &= -\frac{1}{2}\left( \frac{\partial\ln |\mathbf{V}|}{\partial{\kappa_{i}}} + \frac{\partial\ln |\mathbf{X^{\prime}V^{-1}X}|}{\partial{\kappa_{i}}}+\frac{\partial\mathbf{y}^{\prime} \mathbf{P} \mathbf{y} }{\partial{\kappa_{i}}} \right) \end{aligned} \] 我们定义 \(\dot{\mathbf{V}}_{i} = \frac{\partial \mathbf{V}}{\partial \kappa_{i}}\) ，并且使用下面的求导公式 \[ \begin{aligned} \frac{\partial \log |\mathbf{V}|}{\partial \kappa_i} & =\operatorname{tr}\left(\mathbf{V}^{-1} \dot{\mathbf{V}}_i\right) \\ \frac{\partial \mathbf{V}^{-1}}{\partial \kappa_i} & =-\mathbf{V}^{-1} \dot{\mathbf{V}}_i \mathbf{V}^{-1} \end{aligned} \] 因此我们得到 \[ \begin{aligned} & \frac{\partial \log |\mathbf{V}|}{\partial \kappa_i}+\frac{\partial \log \left|\boldsymbol{X}^{\prime} \mathbf{V}^{-1} \boldsymbol{X}\right|}{\partial \kappa_i} \\ & =\operatorname{tr}\left(\mathbf{V}^{-1} \dot{\mathbf{V}}_i\right)-\operatorname{tr}\left(\left(\boldsymbol{X}^{\prime} \mathbf{V}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\prime} \mathbf{V}^{-1} \dot{\mathbf{V}}_i \mathbf{V}^{-1} \boldsymbol{X}\right) \\ & =\operatorname{tr}\left(\mathbf{V}^{-1} \dot{\mathbf{V}}_i-\mathbf{V}^{-1} \boldsymbol{X}\left(\boldsymbol{X}^{\prime} \mathbf{V}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\prime} \mathbf{V}^{-1} \dot{\mathbf{V}}_i\right) \\ & =\operatorname{tr}\left(\boldsymbol{P} \dot{\mathbf{V}}_i\right) \end{aligned} \]

\[ \begin{aligned} \frac{\partial \boldsymbol{y}^{\prime} \boldsymbol{P} \boldsymbol{y}}{\partial \kappa_i} & =\boldsymbol{y}^{\prime} \frac{\partial \boldsymbol{P}}{\partial \kappa_i} \boldsymbol{y} \\ & =-\boldsymbol{y}^{\prime} \boldsymbol{P} \dot{\boldsymbol{V}}_i \boldsymbol{P} \boldsymbol{y} \end{aligned} \]

因此我们得到 REML score 式子： \[ U(\kappa_{i}) = - \frac{1}{2} \Big( \operatorname{tr}\left(\mathbf{P} \dot{\mathbf{V}}_{i}\right)-\mathbf{y}^{\prime} \mathbf{P} \dot{\mathbf{V}}_{i} \mathbf{P} \mathbf{y} \Big) \] 下面我们要考虑 REML score 对于 \(\boldsymbol{\gamma}\) 和 \(\boldsymbol{\phi}\) 的具体式子。

对于一个与 \(\mathbf{G}_{i}\) 相关的方差组分 \(\gamma_{ij}\) ，\(\mathbf{V}_{ij}\) 推导如下 \[ \begin{aligned} \boldsymbol{V}_{i j} & =\frac{\partial \boldsymbol{V}}{\partial \gamma_{i j}} \\ & =\frac{\partial \boldsymbol{Z G Z}^{\prime}}{\partial \gamma_{i j}} \\ & =\frac{\partial}{\partial \gamma_{i j}}\left(\boldsymbol{Z}_1 \boldsymbol{G}_1 \boldsymbol{Z}_1^{\prime}+\ldots+\boldsymbol{Z}_i \boldsymbol{G}_i \boldsymbol{Z}_i^{\prime}+\ldots \boldsymbol{Z}_s \boldsymbol{G}_s \boldsymbol{Z}_s^{\prime}\right) \\ & =\boldsymbol{Z}_i \frac{\partial \boldsymbol{G}_i}{\partial \gamma_{i j}} \boldsymbol{Z}_i^{\prime} \\ & =\boldsymbol{Z}_i \dot{\boldsymbol{G}}_{i j} \boldsymbol{Z}_i^{\prime} \end{aligned} \] 其中 \(\dot{\mathbf{G}}_{i j} = \frac{\partial \mathbf{G}_i}{\partial \gamma_{i j}}\) 。

如果模型中只有一个随机组分，则有： $$ \[\begin{aligned} \boldsymbol{V}_{i} & =\frac{\partial \boldsymbol{V}}{\partial \gamma_{i}} \\ & =\boldsymbol{Z} \frac{\partial \boldsymbol{G}_i}{\partial \gamma_{i}} \boldsymbol{Z}^{\prime} \\ & =\boldsymbol{Z} \dot{\boldsymbol{G}}_{i} \boldsymbol{Z}^{\prime} \end{aligned}\] \[ 对于残差的方差组分 $\phi_{i}$ 则有 \] \[\begin{aligned} &\frac{\partial \boldsymbol{V}}{\partial \phi_{i}} \\ & = \frac{\partial \boldsymbol{\Sigma}}{\partial \phi_{i}} \\ & = \dot{\boldsymbol{\Sigma}}_{i} \end{aligned}\]

\[ 因此我们得到对于 $\gamma_{ij}$ 的 REML score 式子为： \] U(_{ij}) = - ( ( i {i j} _i^{})-{} i {i j} _i^{} ) $$ 我们定义 \(\mathbf{S}_{i}\) 是一个 \((p+q) \times q_{i}\) 的矩阵（\(p,q\) 为所有固定因子和所有随机因子的水平数目，\(q_{i}\) 是第 \(i\) 个随机因子的水平数目 ），并且满足 \(\mathbf{W}\mathbf{S}_{i}=\mathbf{Z}_{i}\) ，因此 \(\mathbf{S}_{i}\) 中只有对应 \(\mathbf{Z}_{i}\) 在 \(\mathbf{W}\) 中的位置存在一个单位矩阵，其它元素全是 0 （因为矩阵向量乘积可以理解为矩阵的列的线性组合，因此 \(\mathbf{S}_{i}\) 每一列只有一个对应 \(\mathbf{Z}_{i}\) 的列号的元素为1，其它为0 ）。

我们将系数矩阵 \(\mathbf{C}\) 写为 \[ \mathbf{C} = \mathbf{W^{\prime}R^{-1}W} + \mathbf{G}^{*} \] 其中 \[ \mathbf{G}^{*} = \left[\begin{array}{c} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{G}^{-1} \\ \end{array}\right] \] 我们对上面式子进一步推导，得到 $$ \[\begin{aligned} &\operatorname{tr}\left( \mathbf{PZ_{i}\dot{\mathbf{G}}_{i j}Z_{i}^{\prime}}\right) \\ &= \operatorname{tr}\left( \mathbf{\dot{G}_{i j}Z_{i}^{\prime}PZ_{i}}\right) \\ &= \operatorname{tr}\left( \mathbf{\dot{G}_{i j}S_{i}^{\prime}W^{\prime}PWS_{i}}\right) \\ &= \operatorname{tr}\left( \mathbf{\dot{G}_{i j}S_{i}^{\prime}W^{\prime}(R^{-1}-R^{-1}WC^{-1}W^{\prime}R^{-1})WS_{i}}\right) \\ &= \operatorname{tr}\left( \mathbf{\dot{G}_{i j}S_{i}^{\prime}(W^{\prime}R^{-1}W-W^{\prime}R^{-1}WC^{-1}W^{\prime}R^{-1}W)S_{i}}\right) \\ &= \operatorname{tr}\left( \mathbf{\dot{G}_{i j}S_{i}^{\prime}W^{\prime}R^{-1}WC^{-1}(C-W^{\prime}R^{-1}W)S_{i}}\right) \\ &= \operatorname{tr}\left( \mathbf{\dot{G}_{i j}S_{i}^{\prime}(C-G^{*})C^{-1}G^{*}S_{i}}\right) \\ &= \operatorname{tr}\left( \mathbf{\dot{G}_{i j}S_{i}^{\prime}G^{*}S_{i}}\right) - \operatorname{tr}\left( \mathbf{\dot{G}_{i j}S_{i}^{\prime}G^{*}C^{-1}G^{*}S_{i}}\right) \\ \end{aligned}\] \[ 其中 $\mathbf{G^{*}C^{-1}G^{*}}$ 可以写为： \] \[\begin{aligned} \mathbf{G^{*}C^{-1}G^{*}} &= \left[\begin{array}{c} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{G}^{-1} \\ \end{array}\right] \left[\begin{array}{c} \mathbf{C}^{XX} & \mathbf{C}^{XZ} \\ \mathbf{C}^{ZX} & \mathbf{C}^{ZZ} \\ \end{array}\right] \left[\begin{array}{c} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{G}^{-1} \\ \end{array}\right] \\ &= \left[\begin{array}{c} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{G}^{-1} \mathbf{C}^{ZZ} \mathbf{G}^{-1} \\ \end{array}\right] \end{aligned}\]

\[ 这里 $\mathbf{C}^{ZZ}$ 就是$\mathbf{C}^{-1}$ 中对应着随机效应的部分，由于我们可以有多个随机效应 ，因此 $\mathbf{C}^{ZZ}$ 可以进一步分块为（其中 $s$ 为随机效应数目）： \] ^{ZZ} = $$ 其中 \(\mathbf{C}^{Z_{i}Z_{i}}\) 就是对应着第 \(i\) 个随机效应 \(\mathbf{u}_{i}\) 的部分。

同样的道理，我们可以将 \(\mathbf{G}\) 矩阵分块为： \[ \mathbf{G} = \left[\begin{array}{c} \mathbf{G}_{1} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{G}_{2} & \cdots & \mathbf{0} \\ \vdots & \vdots & \ddots &\vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{G}_{s} \\ \end{array}\right] \] 因此我们得到 \[ \mathbf{G}^{-1}\mathbf{C}^{ZZ}\mathbf{G}^{-1} = \left[\begin{array}{c} \mathbf{G}_{1}^{-1}\mathbf{C}^{Z_{1}Z_{1}}\mathbf{G}_{1}^{-1} & \mathbf{G}_{1}^{-1}\mathbf{C}^{Z_{1}Z_{2}}\mathbf{G}_{2}^{-1} & \cdots & \mathbf{G}_{1}^{-1}\mathbf{C}^{Z_{1}Z_{s}}\mathbf{G}_{s}^{-1} \\ \mathbf{G}_{2}^{-1}\mathbf{C}^{Z_{2}Z_{1}}\mathbf{G}_{1}^{-1} & \mathbf{G}_{2}^{-1}\mathbf{C}^{Z_{2}Z_{2}}\mathbf{G}_{2}^{-1} & \cdots & \mathbf{G}_{2}^{-1}\mathbf{C}^{Z_{2}Z_{s}}\mathbf{G}_{s}^{-1} \\ \vdots & \vdots & \ddots &\vdots \\ \mathbf{G}_{s}^{-1}\mathbf{C}^{Z_{s}Z_{1}}\mathbf{G}_{1}^{-1} & \mathbf{G}_{s}^{-1}\mathbf{C}^{Z_{s}Z_{2}}\mathbf{G}_{2}^{-1} & \cdots & \mathbf{G}_{s}^{-1}\mathbf{C}^{Z_{s}Z_{s}}\mathbf{G}_{s}^{-1} \\ \end{array}\right] \] 因为 \(\mathbf{S}_{i}\) 中只有对应 \(\mathbf{Z}_{i}\) 在 \(\mathbf{W}\) 中的位置存在一个单位矩阵，其它元素全是 0 ，因此我们得到 \[ \mathbf{S}_{i}^{\prime}\mathbf{G^{*}C^{-1}G^{*}}\mathbf{S}_{i} = \mathbf{G}_{i}^{-1}\mathbf{C}^{Z_{i}Z_{i}}\mathbf{G}_{i}^{-1} \]

\[ \mathbf{S}_{i}^{\prime}\mathbf{G}^{*}\mathbf{S}_{i} = \mathbf{G}_{i}^{-1} \]

因此，\(\operatorname{tr}\left( \mathbf{PZ_{i}\dot{\mathbf{G}}_{i j}Z_{i}^{\prime}}\right)\) 可以进一步推导得到 $$ \[\begin{aligned} &\operatorname{tr}\left( \mathbf{PZ_{i}\dot{\mathbf{G}}_{i j}Z_{i}^{\prime}}\right)\\ &= \operatorname{tr}\left( \mathbf{\dot{G}_{i j}S_{i}^{\prime}G^{*}S_{i}}\right) - \operatorname{tr}\left( \mathbf{\dot{G}_{i j}S_{i}^{\prime}G^{*}C^{-1}G^{*}S_{i}}\right) \\ &= \operatorname{tr}\left( \mathbf{\dot{G}_{i j}G^{-1}_{i}}\right) - \operatorname{tr}\left( \mathbf{\dot{G}_{i j}}\mathbf{G}_{i}^{-1}\mathbf{C}^{Z_{i}Z_{i}}\mathbf{G}_{i}^{-1}\right) \\ \end{aligned}\] \[ 对第二个式子推导前，我们先用 $\mathbf{P}$ 得到随机效应 $\mathbf{u}$ 的式子，根据BLUP方法，我们有 \] = \[ 同时我们有 \] = \[ 因此，我们得到 \] \[\begin{aligned} \mathbf{u}&=\mathbf{GZ^{\prime}Py} \\ &=\left[\begin{array}{c} \mathbf{G}_{1} & \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{G}_{2} & \cdots & \mathbf{0} \\ \vdots & \vdots & \ddots &\vdots \\ \mathbf{0} & \mathbf{0} & \cdots & \mathbf{G}_{s} \\ \end{array}\right] \left[\begin{array}{c} \mathbf{Z}_{1}^{\prime} \\ \mathbf{Z}_{2}^{\prime} \\ \vdots \\ \mathbf{Z}_{s}^{\prime} \\ \end{array}\right]\mathbf{Py} \\ &=\left[\begin{array}{c} \mathbf{G}_{1}\mathbf{Z}_{1}^{\prime} \mathbf{Py} \\ \mathbf{G}_{2}\mathbf{Z}_{2}^{\prime} \mathbf{Py} \\ \vdots \\ \mathbf{G}_{s}\mathbf{Z}_{s}^{\prime}\mathbf{Py} \\ \end{array}\right] \\ \end{aligned}\] \[ 我们同时得到 \] {i} = {i}_{i}^{} \[ 因此，对右手项推导得到 \] \[\begin{aligned} &\mathbf{y}^{\prime} \mathbf{P} \mathbf{Z}_{i} \dot{\mathbf{G}}_{i j} \mathbf{Z}_{i}^{\prime} \mathbf{P} \mathbf{y} \\ &= \mathbf{y}^{\prime} \mathbf{P} \mathbf{Z}_{i} \mathbf{G}_{i} \mathbf{G}_{i}^{-1} \dot{\mathbf{G}}_{i j} \mathbf{G}_{i}^{-1} \mathbf{G}_{i} \mathbf{Z}_{i}^{\prime} \mathbf{P} \mathbf{y} \\ &= \mathbf{u}_{i}^{\prime} \mathbf{G}_{i}^{-1} \dot{\mathbf{G}}_{i j} \mathbf{G}_{i}^{-1} \mathbf{u}_{i}\\ \end{aligned}\] \[ 因此，我们得到对于 $\gamma_{ij}$ 的 REML score 式子可以写为： \] U({ij}) = - ( ( ) - ( ^{{Z_{i}Z_{i}}_{i}}{-1}{i}^{-1})-_{i}{} {i}^{-1} {i j} {i}^{-1} {i} ) \[ 如果模型中只有一个随机效应，则有 \] U({i}) = - ( ( ) - ( ^{ZZ}{-1}^{-1})-{} ^{-1} {i} ^{-1} ) \[ 对于残差的方差组分 $\phi_{i}$ ，我们将 $\dot{\mathbf{V}}_{i} = \dot{\mathbf{\Sigma}}_{i}$ 代入 $U(\kappa_{i}) = - \frac{1}{2} \Big( \operatorname{tr}\left(\mathbf{P} \dot{\mathbf{V}}_{i}\right)-\mathbf{y}^{\prime} \mathbf{P} \dot{\mathbf{V}}_{i} \mathbf{P} \mathbf{y} \Big)$ 得到 \] U({i}) = - ( ( {i})-^{} _{i} ) \[ 我们对其进行推导得到： \] \[\begin{aligned} &\operatorname{tr}\left(\mathbf{P} \dot{\mathbf{\Sigma}}_{i}\right) \\ &= \operatorname{tr} \left( \left(\mathbf{R^{-1}-R^{-1}WC^{-1}W^{\prime}R^{-1}} \right) \dot{\mathbf{\Sigma}}_{i}\right) \\ &= \operatorname{tr} \left( \left(\mathbf{\Sigma^{-1}-\Sigma^{-1}WC^{-1}W^{\prime}\Sigma^{-1}} \right) \dot{\mathbf{\Sigma}}_{i}\right) \\ &= \operatorname{tr} \left(\mathbf{\Sigma^{-1}\dot{\mathbf{\Sigma}}_{i}-\Sigma^{-1}WC^{-1}W^{\prime}\Sigma^{-1}}\dot{\mathbf{\Sigma}}_{i} \right) \\ &\\ \end{aligned}\] \[ 对于第二个式子，首先我们推导 $\mathbf{Py}=\mathbf{R^{-1}e}$ ，如下 \] \[\begin{aligned} &\mathbf{Py} \\ &= \mathbf{V^{-1}}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}) \\ &= \left( \mathbf{R^{-1}-R^{-1}Z\left(Z^{\prime}R^{-1}Z+G^{-1}\right)^{-1}Z^{\prime}R^{-1}} \right)(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}) \\ &= \mathbf{R^{-1}}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta})- \mathbf{R^{-1}Z\left(Z^{\prime}R^{-1}Z+G^{-1}\right)^{-1}Z^{\prime}R^{-1}}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}) \\ &= \mathbf{R^{-1}}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta})- \mathbf{R^{-1}Z}\mathbf{u} \quad \because \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{\beta}) \\ &= \mathbf{R^{-1}}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta} -\mathbf{Z}\mathbf{u}) \\ &= \mathbf{R^{-1}}\mathbf{e} \\ &= \mathbf{\Sigma^{-1}}\mathbf{e} \\ \end{aligned}\]

\[ 因此，我们得到 \] ^{} {i} = ^{}\mathbf{{-1}} {i} \[ 因此我们得到关于 $\phi_{i}$ 的 REML score 式子如下 \] U({i}) = - ( () -({i} )-^{}\mathbf{{-1}} _{i} ) $$

单性状模型的迭代公式

对于单性状模型（假设所有随机因子之间均不存在相关），我们将 \(U(\gamma_{ij})\) 设为 0， 推导一下，我们得到： \[ q_{i}/\sigma_{i}^{2} - \operatorname{tr}\left( \mathbf{C}^{Z_{i}Z_{i}}\mathbf{A}_{i}^{-1}\right)/\sigma_{i}^{4} = \mathbf{u}_{i}^{\prime} \mathbf{A}_{i}^{-1} \mathbf{u}_{i}/\sigma_{i}^{4} \] 整理一下，我们得到（其中 \(q_{i}\) 为 \(u_{i}\) 的水平数目，或者说是 \(\mathbf{Z}_{i}\) 的列数） \[ \hat{\sigma}_{i}^{2}=\left(\hat{\mathbf{u}}_{i}^{\prime} \mathbf{A}_{i}^{-1} \hat{\mathbf{u}}_{i}+\operatorname{tr}\left( \mathbf{C}^{Z_{i}Z_{i}}\mathbf{A}_{i}^{-1}\right)\right) / q_{i} \] 一个收敛更快的迭代公式为 \[ \hat{\sigma}_{i}^{2}=\frac{\hat{\mathbf{u}}_{i}^{\prime} \mathbf{A}_{i}^{-1} \hat{\mathbf{u}}_{i}}{q_{i}- \operatorname{tr}\left( \mathbf{C}^{Z_{i}Z_{i}}\mathbf{A}_{i}^{-1}\right)/\hat{\sigma}_{i}^{2} } \] 联合使用之前推导得到的残差计算公式(其中 \(N\) 是表型数目，\(r\) 为 \(\mathbf{X}\) 的秩，只适合单性状模型)，我们便得到了完整的 EM-REML 估计公式 \[ \hat{\sigma}_{0}^{2}=\left(\mathbf{y}^{\prime} \mathbf{y}-\mathbf{y}^{\prime} \mathbf{X} \hat{\boldsymbol{\beta}}-\mathbf{y}^{\prime} \mathbf{Z} \hat{\mathbf{u}}\right) /(N-r) \] 这里我们还有残差的另外两种推导公式，如下 \[ \begin{aligned} \sigma_{0}^{2} &= \frac{1}{n}(\mathbf{e}^{\prime}\mathbf{e}+ \operatorname{tr}(\mathbf{WC^{-1}W^{\prime}})) \\ \text{or}&\\ \sigma_{0}^{2} &= \frac{1}{n}\left(\sigma_{0}^{2}\left(p+q\right)- \sigma_{0}^{2} \operatorname{tr}\left(\mathbf{G^{-1}C^{ZZ}}\right) + \mathbf{e}^{\prime}\mathbf{e}\right) \end{aligned} \] 将 \(U(\phi_{i})\) 设为0，我们得到 \[ \frac{n}{\sigma_{0}^{2}}-\frac{\operatorname{tr}\left(\mathbf{WC^{-1}W^{\prime}} \right)}{\sigma_{0}^{4}} = \frac{\mathbf{e^{\prime}e}}{\sigma_{0}^{4}} \] 因此可以得到第一个式子，我们对其中 \(\operatorname{tr}\left(\mathbf{WC^{-1}W^{\prime}} \right)/\sigma_{0}^{2}\) 推导如下 \[ \begin{aligned} &\operatorname{tr}\left(\mathbf{WC^{-1}W^{\prime}} \right)/\sigma_{0}^{2} \\ &=\operatorname{tr}\left(\mathbf{W^{\prime}R^{-1}WC^{-1}} \right) \\ &=\operatorname{tr}\left(\mathbf{(C-G^{*})C^{-1}} \right) \quad \because \mathbf{C} = \mathbf{W^{\prime}R^{-1}W} + \mathbf{G}^{*} \\ &=\operatorname{tr}\left(\mathbf{I_{p+q}-G^{*}C^{-1}} \right)\\ &= p+q-\operatorname{tr}\left(\mathbf{G^{*}C^{-1}} \right)\\ \end{aligned} \] 对于 \(\mathbf{G^{*}C^{-1}}\) ，我们对其推导得到 \[ \begin{aligned} \mathbf{G^{*}C^{-1}} &= \left[\begin{array}{c} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{G}^{-1} \\ \end{array}\right] \left[\begin{array}{c} \mathbf{C}^{XX} & \mathbf{C}^{XZ} \\ \mathbf{C}^{ZX} & \mathbf{C}^{ZZ} \\ \end{array}\right] \\ &= \left[\begin{array}{c} \mathbf{0} & \mathbf{0} \\ \mathbf{G}^{-1}\mathbf{C}^{ZX} & \mathbf{G}^{-1}\mathbf{C}^{ZZ} \\ \end{array}\right] \\ \end{aligned} \] 因此，我们得到 \[ \operatorname{tr}(\mathbf{G^{*}C^{-1}}) = \operatorname{tr}(\mathbf{G^{-1}C^{ZZ}}) \] 因此我们得到 \[ \operatorname{tr}\left(\mathbf{WC^{-1}W^{\prime}} \right)/\sigma_{0}^{2} = p+q-\operatorname{tr}(\mathbf{G^{-1}C^{ZZ}}) \] 代入上式，我们得到 \[ \frac{n-p-q+\operatorname{tr}(\mathbf{G^{-1}C^{ZZ}})}{\sigma_{0}^{2}} = \frac{\mathbf{e^{\prime}e}}{\sigma_{0}^{4}} \] 因此可以得到第二个式子。

多性状模型的迭代公式

随机效应

对于多性状模型，假设只有一个随机因子（其实多个因子同理，只是这样不用写 \(\mathbf{G}\) 矩阵的下标了，符号更加清晰），按照性状排序， 则有 \[ \mathbf{G} = \mathbf{G}_{k} \otimes \mathbf{A} \] 其中 \(\mathbf{G}_{k}\) 就是相应的方差组分的矩阵形式，如下所示，即反过来 \(\boldsymbol{\gamma} = \mathrm{vec}(\mathbf{G}_{k})\) ，其中 \(k\) 为性状数目。 \[ \boldsymbol{G}_k=\left[\begin{array}{cccc} \gamma_{11} & \gamma_{12} & \cdots & \gamma_{1 k} \\ \gamma_{12} & \gamma_{22} & \cdots & \gamma_{2 k} \\ \vdots & \vdots & \ddots & \vdots \\ \gamma_{1 k} & \gamma_{2 k} & \cdots & \gamma_{k k} \end{array}\right] \] 对于随机因子 \(\mathbf{u}\) ，假设对于每个性状的随机因子水平数目均为 \(m\) ，则 \(\mathbf{u}\) 也可以表示为一个 \(m \times k\) 的矩阵 \(\mathbf{U}\) ，其中每一列代表一个性状的所有水平，如下 \[ \mathbf{U} =\left[\begin{array}{cccc} u_{11} & u_{12} & \cdots & u_{1 k} \\ u_{21} & u_{22} & \cdots & u_{2 k} \\ \vdots & \vdots & \ddots & \vdots \\ u_{m1} & u_{m2} & \cdots & u_{m k} \end{array}\right] = \left[\begin{array}{cccc} \mathbf{u}_{1} & \mathbf{u}_{2} & \cdots & \mathbf{u}_{ k} \\ \end{array}\right] \] 反过来则有 \(\mathbf{u} = \mathrm{vec} (\mathbf{U})\) 。

对于 \(\gamma_{ij}\) ，我们有公式 \[ U(\gamma_{ij}) = - \frac{1}{2} \Big( \operatorname{tr}\left( \mathbf{G^{-1}\dot{G}_{i j}}\right) - \operatorname{tr}\left( \mathbf{C}^{ZZ}\mathbf{G}^{-1}\mathbf{\dot{G}_{i j}}\mathbf{G}_{i}^{-1}\right)-\mathbf{u}^{\prime} \mathbf{G}^{-1} \dot{\mathbf{G}}_{i j} \mathbf{G}^{-1} \mathbf{u} \Big) \] 因为 \(\mathbf{\dot{G}_{i j}}\) 相应 \((m,n)\) 位置的元素 就是 \(\partial g_{mn} / \partial \gamma_{ij}\) ，所有易得 \(\mathbf{\dot{G}_{i j}} = \mathbf{\dot{G}_{k_{i j}}} \otimes \mathbf{A}\) ，其中 \(\mathbf{\dot{G}_{k_{i j}}} = \frac{\partial \mathbf{G}_{k}}{\partial \gamma_{ij}}\) 。

因此，我们推导得到。 $$ \[\begin{aligned} & \mathbf{G^{-1}\dot{G}_{i j}} \\ &= \left( \mathbf{G}_{k}^{-1} \otimes \mathbf{A}^{-1} \right) \left(\mathbf{\dot{G}}_{k_{i j}} \otimes \mathbf{A}\right) \\ &= \left( \mathbf{G}_{k}^{-1} \mathbf{\dot{G}}_{k_{i j}} \otimes \mathbf{I}_{m} \right) &\\ &\\ & \operatorname{tr} (\mathbf{G^{-1}\dot{G}_{i j}}) \\ &= \operatorname{tr}\left( \mathbf{G}_{k}^{-1} \mathbf{\dot{G}}_{k_{i j}} \otimes \mathbf{I}_{m} \right) \\ &= m \operatorname{tr}\left( \mathbf{G}_{k}^{-1} \mathbf{\dot{G}}_{k_{i j}} \right) \\ &\\ &\\ &\mathbf{G}^{-1} \dot{\mathbf{G}}_{i j} \mathbf{G}^{-1} \\ &= \left( \mathbf{G}_{k}^{-1} \otimes \mathbf{A}^{-1} \right) \left(\mathbf{\dot{G}}_{k_{i j}} \otimes \mathbf{A}\right) \left( \mathbf{G}_{k}^{-1} \otimes \mathbf{A}^{-1} \right)\\ &= \left( \mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} \otimes \mathbf{A}^{-1} \right) &\\ &\\ &\mathbf{u}^{\prime} \mathbf{G}^{-1} \dot{\mathbf{G}}_{i j} \mathbf{G}^{-1} \mathbf{u} \\ &= \mathrm{vec}(\mathbf{U})^{\prime} \left( \mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} \otimes \mathbf{A}^{-1} \right) \mathrm{vec}(\mathbf{U}) \\ &= \mathrm{tr}(\mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} \mathbf{U}^{\prime} \mathbf{A}^{-1} \mathbf{U}) \quad \because \mathrm{vec}(\mathbf{A})^{\prime} (\mathbf{B} \otimes \mathbf{C}) \mathrm{vec}(\mathbf{A}) = \mathrm{tr}(\mathbf{B^{\prime}A^{\prime}CA}) \\ &= \mathrm{tr}(\mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} \mathbf{\Delta} ) \\ \end{aligned}\]

$$

最后一个式子中，我们设 \(\mathbf{\Delta} = \mathbf{U}^{\prime} \mathbf{A}^{-1} \mathbf{U}\) 。

代入 \(U(\gamma_{ij})\) 得到 \[ \begin{aligned} U(\gamma_{ij}) &= - \frac{1}{2} \Big( m \operatorname{tr}\left( \mathbf{G}_{k}^{-1} \mathbf{\dot{G}}_{k_{i j}} \right) - \operatorname{tr}\left(\mathbf{C}^{ZZ} \left( \mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} \otimes \mathbf{A}^{-1} \right) \right)-\mathrm{tr}(\mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} \mathbf{\Delta} ) \Big)\\ \end{aligned} \] 我们需要对这个式子进一步简化，我们知道 \(\mathbf{G}_{k}\) 是一个对称矩阵，因此 \(\mathbf{G}_{k}^{-1}\) 也是一个对称矩阵， 我们将 \(\mathbf{G}_{k}^{-1}\) 中的元素表示为 \(\gamma^{ij}\) (\(i \leq j\) ) 。

首先我们考虑 \(\operatorname{tr}\left( \mathbf{G}_{k}^{-1} \mathbf{\dot{G}}_{k_{i j}} \right)\) 这一项，当 \(i=j\) 时，矩阵 \(\mathbf{\dot{G}}_{k_{i j}}\) 只有 \((i,i)\) 位置元素为1，其它位置元素均为 0 。因此 \(\mathbf{G}_{k}^{-1} \mathbf{\dot{G}}_{k_{i j}}\) 只有第 \(i\) 列元素为 \(\mathbf{G}_{k}^{-1}\) 的第 \(i\) 列元素，其它元素为 0 。因此，我们得到 \[ \operatorname{tr}\left( \mathbf{G}_{k}^{-1} \mathbf{\dot{G}}_{k_{i j}} \right) = \gamma^{ii} \] 当 \(i \neq j\) 时, 矩阵 \(\mathbf{\dot{G}}_{k_{i j}}\) 只有 \((i,j)\) 和 \((j,i)\) 位置为 1，其它位置元素为 0 。因此 \(\mathbf{G}_{k}^{-1} \mathbf{\dot{G}}_{k_{i j}}\) 其第 \(i\) 列为 \(\mathbf{G}_{k}^{-1}\) 的第 \(j\) 列， 其第 \(j\) 列为 \(\mathbf{G}_{k}^{-1}\) 的第 \(i\) 列，其它位置元素为 0 。因此我们得到 \[ \operatorname{tr}\left( \mathbf{G}_{k}^{-1} \mathbf{\dot{G}}_{k_{i j}} \right) = \gamma^{ij} + \gamma^{ji} = 2 \gamma^{ij} \] 即我们得到 \[ \operatorname{tr}\left(\boldsymbol{G}_k^{-1} \dot{\boldsymbol{G}}_{k_{i j}}\right)=\left\{\begin{array}{cc} \gamma^{i i}, & i=j \\ 2 \gamma^{i j}, & i \neq j \end{array}\right. \] 我们现在考虑第三项，如下 \[ \begin{aligned} \mathrm{tr}(\mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} \mathbf{\Delta}) & = \mathrm{tr}(\mathbf{G}_{k}^{-1} \mathbf{\Delta}\mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}) \\ & = \mathrm{tr}(\mathbf{\Delta}^{*}\mathbf{\dot{G}}_{k_{i j}}) \\ \end{aligned} \] 其中我们设 \(\mathbf{\Delta}^{*} = \mathbf{G}_{k}^{-1} \mathbf{\Delta}\mathbf{G}_{k}^{-1} = \{ \delta_{ij}^{*}\}\) 。

当 \(i=j\) 时，矩阵 \(\mathbf{\Delta}^{*}\mathbf{\dot{G}}_{k_{i j}}\) 只有第 \(i\) 列元素为 \(\mathbf{\Delta}^{*}\) 的第 \(i\) 列元素 ，其它位置元素为 0 。因此，我们得到 \[ \mathrm{tr}(\mathbf{\Delta}^{*}\mathbf{\dot{G}}_{k_{i j}}) = \delta_{ii}^{*} \] 当 \(i \neq j\) 时, 矩阵 \(\mathbf{\Delta}^{*}\mathbf{\dot{G}}_{k_{i j}}\) 其第 \(i\) 列为 \(\mathbf{\Delta}^{*}\) 的第 \(j\) 列， 其第 \(j\) 列为 \(\mathbf{\Delta}^{*}\) 的第 \(i\) 列，其它位置元素为 0 。因此，我们得到 (因为 \(\mathbf{\Delta}^{*}\) 对称) \[ \mathrm{tr}(\mathbf{\Delta}^{*}\mathbf{\dot{G}}_{k_{i j}}) = \delta_{ij}^{*} + \delta_{ji}^{*} = 2 \delta_{ij}^{*} \] 综上，我们有 \[ \mathrm{tr}(\mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} \mathbf{\Delta})=\left\{\begin{array}{cc} \delta^{*}_{i i}, & i=j \\ 2 \delta^{*}_{i j}, & i \neq j \end{array}\right. \] 我们现在考虑第二项，我们将 \(\mathbf{C}^{ZZ}\) 按照性状分块为 \(k^{2}\) 个 \(m \times m\) 的子矩阵，如下 \[ \mathbf{C}^{ZZ} =\left[\begin{array}{cccc} \mathbf{C}^{ZZ}_{11} & \mathbf{C}^{ZZ}_{12} & \cdots & \mathbf{C}^{ZZ}_{1 k} \\ \mathbf{C}^{ZZ}_{21} & \mathbf{C}^{ZZ}_{22} & \cdots & \mathbf{C}^{ZZ}_{2 k} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{C}^{ZZ}_{k1} & \mathbf{C}^{ZZ}_{k2} & \cdots & \mathbf{C}^{ZZ}_{k k} \end{array}\right] \] 处于简化的目的，我们这里考虑 \(k=3\) 的情况（这个不影响推导出来的结果），根据矩阵乘法，当 \(i=j\) 时，我们有 \[ \mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} = \left[\begin{array}{cccc} \gamma^{1i}\gamma^{i1} & \gamma^{1i}\gamma^{i2} & \gamma^{1i}\gamma^{i3} \\ \gamma^{2i}\gamma^{i1} & \gamma^{2i}\gamma^{i2} & \gamma^{2i}\gamma^{i3} \\ \gamma^{3i}\gamma^{i1} & \gamma^{3i}\gamma^{i2} & \gamma^{3i}\gamma^{i3} \end{array}\right] \] 因此，我们得到 \[ \begin{aligned} & \mathbf{C}^{ZZ} \left( \mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} \otimes \mathbf{A}^{-1} \right) \\ &= \left[\begin{array}{cccc} \mathbf{C}^{ZZ}_{11} & \mathbf{C}^{ZZ}_{12} & \mathbf{C}^{ZZ}_{1 3} \\ \mathbf{C}^{ZZ}_{21} & \mathbf{C}^{ZZ}_{22} & \mathbf{C}^{ZZ}_{2 3} \\ \mathbf{C}^{ZZ}_{31} & \mathbf{C}^{ZZ}_{32} & \mathbf{C}^{ZZ}_{3 3} \end{array}\right] \left[\begin{array}{cccc} \gamma^{1i}\gamma^{i1} \mathbf{A}^{-1} & \gamma^{1i}\gamma^{i2} \mathbf{A}^{-1}& \gamma^{1i}\gamma^{i3} \mathbf{A}^{-1}\\ \gamma^{2i}\gamma^{i1} \mathbf{A}^{-1}& \gamma^{2i}\gamma^{i2}\mathbf{A}^{-1} & \gamma^{2i}\gamma^{i3} \mathbf{A}^{-1}\\ \gamma^{3i}\gamma^{i1} \mathbf{A}^{-1}& \gamma^{3i}\gamma^{i2} \mathbf{A}^{-1} & \gamma^{3i}\gamma^{i3} \mathbf{A}^{-1} \end{array}\right] \\ \end{aligned} \] 利用迹函数的性质，我们得到 \[ \begin{aligned} \operatorname{tr}\left(\mathbf{C}^{ZZ} \left( \mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} \otimes \mathbf{A}^{-1} \right) \right) &= \sum_{s=1}^{3}\sum_{t=1}^{3} \operatorname{tr} (\mathbf{C}^{ZZ}_{st} \gamma^{ti}\gamma^{is}\mathbf{A}^{-1}) \\ &= \sum_{s=1}^{3}\sum_{t=1}^{3}\gamma^{it}\gamma^{si} \operatorname{tr} (\mathbf{C}^{ZZ}_{st} \mathbf{A}^{-1}) \\ \end{aligned} \] 当 \(i \neq j\) 时, 我们有 \[ \mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} = \left[\begin{array}{cccc} \gamma^{1j}\gamma^{i1}+\gamma^{1i}\gamma^{j1} & \gamma^{1j}\gamma^{i2}+\gamma^{1i}\gamma^{j2} & \gamma^{1j}\gamma^{i3}+\gamma^{1i}\gamma^{j3} \\ \gamma^{2j}\gamma^{i1}+\gamma^{2i}\gamma^{j1} & \gamma^{2j}\gamma^{i2}+\gamma^{2i}\gamma^{j2} & \gamma^{2j}\gamma^{i3}+\gamma^{2i}\gamma^{j3} \\ \gamma^{3j}\gamma^{i1}+\gamma^{3i}\gamma^{j1} & \gamma^{3j}\gamma^{i2}+\gamma^{3i}\gamma^{j2} & \gamma^{3j}\gamma^{i3}+\gamma^{3i}\gamma^{j3} \end{array}\right] \] 因此 $$ \[\begin{aligned} & \mathbf{C}^{ZZ} \left( \mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} \otimes \mathbf{A}^{-1} \right) \\ &= \left[\begin{array}{cccc} \mathbf{C}^{ZZ}_{11} & \mathbf{C}^{ZZ}_{12} & \mathbf{C}^{ZZ}_{1 3} \\ \mathbf{C}^{ZZ}_{21} & \mathbf{C}^{ZZ}_{22} & \mathbf{C}^{ZZ}_{2 3} \\ \mathbf{C}^{ZZ}_{31} & \mathbf{C}^{ZZ}_{32} & \mathbf{C}^{ZZ}_{3 3} \end{array}\right] \left[\begin{array}{cccc} (\gamma^{1j}\gamma^{i1}+\gamma^{1i}\gamma^{j1})\mathbf{A}^{-1} & (\gamma^{1j}\gamma^{i2}+\gamma^{1i}\gamma^{j2})\mathbf{A}^{-1} & (\gamma^{1j}\gamma^{i3}+\gamma^{1i}\gamma^{j3})\mathbf{A}^{-1} \\ (\gamma^{2j}\gamma^{i1}+\gamma^{2i}\gamma^{j1})\mathbf{A}^{-1} & (\gamma^{2j}\gamma^{i2}+\gamma^{2i}\gamma^{j2})\mathbf{A}^{-1} & (\gamma^{2j}\gamma^{i3}+\gamma^{2i}\gamma^{j3})\mathbf{A}^{-1} \\ (\gamma^{3j}\gamma^{i1}+\gamma^{3i}\gamma^{j1})\mathbf{A}^{-1} & (\gamma^{3j}\gamma^{i2}+\gamma^{3i}\gamma^{j2})\mathbf{A}^{-1} & (\gamma^{3j}\gamma^{i3}+\gamma^{3i}\gamma^{j3})\mathbf{A}^{-1} \end{array}\right] \\ \end{aligned}\] \[ 因此我们得到 \] \[\begin{aligned} \operatorname{tr}\left(\mathbf{C}^{ZZ} \left( \mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} \otimes \mathbf{A}^{-1} \right) \right) &= \sum_{s=1}^{3}\sum_{t=1}^{3} \operatorname{tr} \left(\mathbf{C}^{ZZ}_{st} \left(\gamma^{tj}\gamma^{is}+\gamma^{ti}\gamma^{js}\right)\mathbf{A}^{-1}\right) \\ &= \sum_{s=1}^{3}\sum_{t=1}^{3} \left(\gamma^{tj}\gamma^{is}+\gamma^{ti}\gamma^{js}\right) \operatorname{tr} \left(\mathbf{C}^{ZZ}_{st} \mathbf{A}^{-1}\right) \\ &= \sum_{s=1}^{3}\sum_{t=1}^{3}\gamma^{tj}\gamma^{is}\operatorname{tr} \left(\mathbf{C}^{ZZ}_{st} \mathbf{A}^{-1}\right) + \sum_{s=1}^{3}\sum_{t=1}^{3}\gamma^{ti}\gamma^{js}\operatorname{tr} \left(\mathbf{C}^{ZZ}_{st} \mathbf{A}^{-1}\right) \end{aligned}\] \[ 我们对需要求和的第一项交换 $s$ 和 $t$ ，并且我们知道 $\gamma^{ij} = \gamma^{ji}$ , $\mathbf{C}^{ZZ}_{st} = (\mathbf{C}^{ZZ}_{ts})^{T}$ ，因此我们继续推导得到 \] \[\begin{aligned} \operatorname{tr}\left(\mathbf{C}^{ZZ} \left( \mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} \otimes \mathbf{A}^{-1} \right) \right) &= \sum_{s=1}^{3}\sum_{t=1}^{3}\gamma^{tj}\gamma^{is}\operatorname{tr} \left(\mathbf{C}^{ZZ}_{st} \mathbf{A}^{-1}\right) + \sum_{s=1}^{3}\sum_{t=1}^{3}\gamma^{ti}\gamma^{js}\operatorname{tr} \left(\mathbf{C}^{ZZ}_{st} \mathbf{A}^{-1}\right) \\ & = \sum_{t=1}^{3}\sum_{s=1}^{3}\gamma^{sj}\gamma^{it}\operatorname{tr} \left(\mathbf{C}^{ZZ}_{ts} \mathbf{A}^{-1}\right) + \sum_{s=1}^{3}\sum_{t=1}^{3}\gamma^{ti}\gamma^{js}\operatorname{tr} \left(\mathbf{A}^{-1}\mathbf{C}^{ZZ}_{ts} \right) \\ & = \sum_{s=1}^{3}\sum_{t=1}^{3}\gamma^{sj}\gamma^{it}\operatorname{tr} \left(\mathbf{C}^{ZZ}_{ts} \mathbf{A}^{-1}\right) + \sum_{s=1}^{3}\sum_{t=1}^{3}\gamma^{ti}\gamma^{js}\operatorname{tr} \left(\mathbf{C}^{ZZ}_{ts} \mathbf{A}^{-1}\right) \\ & = 2\sum_{s=1}^{3}\sum_{t=1}^{3}\gamma^{it}\gamma^{sj}\operatorname{tr} \left(\mathbf{C}^{ZZ}_{ts} \mathbf{A}^{-1}\right) \\ \end{aligned}\] \[ 综上，我们将 $k=3$ 改为一般的 $k$ ，我们得到 \] (^{ZZ} ( {k}^{-1}{k_{i j}}_{k}^{-1} ^{-1} ) )={ \[\begin{array}{cc} \sum_{s=1}^{k}\sum_{t=1}^{k}\gamma^{it}\gamma^{si} \operatorname{tr} (\mathbf{C}^{ZZ}_{st} \mathbf{A}^{-1}), & i=j \\ 2\sum_{s=1}^{k}\sum_{t=1}^{k}\gamma^{it}\gamma^{sj}\operatorname{tr} \left(\mathbf{C}^{ZZ}_{st} \mathbf{A}^{-1}\right), & i \neq j \end{array}\] . \[ 现在我们定义 $\mathbf{\Psi}$ 是一个 $k \times k$ 的矩阵，其元素为 $\Psi_{ij} = \operatorname{tr}(\mathbf{C}_{ij}^{ZZ}\mathbf{A}^{-1})$ ，根据矩阵乘法性质（假设 $\mathbf{A,B,C}$ 均为 $n \times n$ 的矩阵，则 $\mathbf{ABC} = \{\sum_{s=1}^{n}\sum_{t=1}^{n}a_{it}b_{ts}c_{sj}\}$ ），我们发现 \] \[\begin{aligned} \mathbf{G}_{k}^{-1} \mathbf{\Psi}\mathbf{G}_{k}^{-1} &= \left \{\sum_{s=1}^{k}\sum_{t=1}^{k}\gamma^{it}\operatorname{tr}(\mathbf{C}_{ts}^{ZZ}\mathbf{A}^{-1})\gamma^{sj} \right\}\\ &= \left \{\sum_{s=1}^{k}\sum_{t=1}^{k}\gamma^{it}\gamma^{sj}\operatorname{tr}(\mathbf{C}_{ts}^{ZZ}\mathbf{A}^{-1}) \right\}\\ \end{aligned}\]

$$ 我们定义 \(\psi_{ij}^{*}\) 是 \(\mathbf{\Psi}^{*} = \mathbf{G}_{k}^{-1} \mathbf{\Psi}\mathbf{G}_{k}^{-1}\) 的元素，即 \(\psi_{ij}^{*} = \left\{ \mathbf{G}_{k}^{-1} \mathbf{\Psi}\mathbf{G}_{k}^{-1} \right\}_{ij}\) 。

因此，我们得到 \[ \operatorname{tr}\left(\mathbf{C}^{ZZ} \left( \mathbf{G}_{k}^{-1}\mathbf{\dot{G}}_{k_{i j}}\mathbf{G}_{k}^{-1} \otimes \mathbf{A}^{-1} \right) \right)=\left\{\begin{array}{cc} \psi_{ii}^{*}, & i=j \\ 2\psi_{ij}^{*}, & i \neq j \end{array}\right. \] 将上面的式子代入 \(U(\gamma_{ij})\) ，进一步简化得到 \[ U(\gamma_{ij}) =\left\{\begin{array}{cc} - \frac{1}{2} \Big(m\gamma^{ii} - \psi_{ii}^{*}- \delta_{ii}^{*} \big), & i=j \\ - \frac{1}{2} \Big(2m\gamma^{ij} - 2\psi_{ij}^{*}- 2\delta_{ij}^{*} \big), & i \neq j \end{array}\right. \] 将之设为0， 对 \(\gamma^{ij}\) 求解得到 EM 推导公式 \[ m\gamma^{ij} = \psi_{ij}^{*} + \delta_{ij}^{*} \] 使用矩阵格式，我们得到 \[ \begin{aligned} m\mathbf{G}_{k}^{-1} &= \mathbf{\Psi}^{*} + \mathbf{\Delta}^{*} \\ &= \mathbf{G}_{k}^{-1} \mathbf{\Psi}\mathbf{G}_{k}^{-1} + \mathbf{G}_{k}^{-1} \mathbf{\Delta}\mathbf{G}_{k}^{-1} \end{aligned} \] 左右均乘以 \(\mathbf{G}_{k}^{-1}\) ，得到 \[ m\mathbf{G}_{k}= \mathbf{\Psi} + \mathbf{\Delta} \] 因此，我们得到 EM 推导公式 \[ m\gamma_{ij}^{(m+1)} = \psi_{ij}^{(m)} + \delta_{ij}^{(m)} \] 因为 \(\psi_{ij} = \operatorname{tr}(\mathbf{C}_{ij}^{ZZ}\mathbf{A}^{-1})\) ， \(\mathbf{\Delta} = \mathbf{U}^{\prime} \mathbf{A}^{-1} \mathbf{U}\) 可以进一步分块为： \[ \mathbf{\Delta} = \mathbf{U}^{\prime} \mathbf{A}^{-1} \mathbf{U} =\left[\begin{array}{cccc} \mathbf{u}_{1}^{\prime}\mathbf{A}^{-1}\mathbf{u}_{1} & \mathbf{u}_{1}^{\prime}\mathbf{A}^{-1}\mathbf{u}_{2} & \cdots & \mathbf{u}_{1}^{\prime}\mathbf{A}^{-1}\mathbf{u}_{k} \\ \mathbf{u}_{2}^{\prime}\mathbf{A}^{-1}\mathbf{u}_{1} & \mathbf{u}_{2}^{\prime}\mathbf{A}^{-1}\mathbf{u}_{2} & \cdots & \mathbf{u}_{2}^{\prime}\mathbf{A}^{-1}\mathbf{u}_{k} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{u}_{k}^{\prime}\mathbf{A}^{-1}\mathbf{u}_{1} & \mathbf{u}_{k}^{\prime}\mathbf{A}^{-1}\mathbf{u}_{2} & \cdots & \mathbf{u}_{k}^{\prime}\mathbf{A}^{-1}\mathbf{u}_{k} \\ \end{array}\right] \] 因此 \(\delta_{ij} = \mathbf{u}_{i}^{\prime}\mathbf{A}^{-1}\mathbf{u}_{j}\) 。因此上面的 EM 推导公式可以写为 \[ \gamma_{ij}^{(m+1)} = \frac{1}{m} \left(\operatorname{tr}(\mathbf{C}_{ij}^{ZZ^{(m)}}\mathbf{A}^{-1}) + \mathbf{u}_{i}^{(m)^{\prime}}\mathbf{A}^{-1}\mathbf{u}_{j}^{(m)} \right) \] 这个式子和单性状模型的推导公式正好具有相同的格式。

残差

按照性状排序， 此时我们有下式成立，其中 \(n\) 为表型行数，\(\mathbf{\Sigma}_{k}\) 为需要估计的协方差组分的矩阵。 \[ \mathbf{\Sigma} = \mathbf{\Sigma}_{k} \otimes \mathbf{I}_{n} \] \(\mathbf{\Sigma}_{k}\) 形式如下（这里我们使用 \(\phi_{ij}\) ，与上面的随机效应格式相同 ）。 \[ \boldsymbol{\Sigma}_k=\left[\begin{array}{cccc} \phi_{11} & \phi_{12} & \cdots & \phi_{1 k} \\ \phi_{12} & \phi_{22} & \cdots & \phi_{2 k} \\ \vdots & \vdots & \ddots & \vdots \\ \phi_{1 k} & \phi_{2 k} & \cdots & \phi_{k k} \end{array}\right] \] 对于随机因子 \(\mathbf{e}\) ，同上，\(\mathbf{e}\) 也可以表示为一个 \(n \times k\) 的矩阵 \(\mathbf{E}\) ，其中每一列代表一个性状的所有残差，如下 \[ \mathbf{E} =\left[\begin{array}{cccc} e_{11} & e_{12} & \cdots & e_{1 k} \\ e_{21} & e_{22} & \cdots & e_{2 k} \\ \vdots & \vdots & \ddots & \vdots \\ e_{n1} & e_{n2} & \cdots & e_{n k} \end{array}\right] = \left[\begin{array}{cccc} \mathbf{e}_{1} & \mathbf{e}_{2} & \cdots & \mathbf{e}_{ k} \\ \end{array}\right] \] 反过来则有 \(\mathbf{e} = \mathrm{vec} (\mathbf{E})\) 。

对于 \(\phi_{ij}\) ，我们有公式 \[ U(\phi_{ij}) = - \frac{1}{2} \Big( \operatorname{tr} \left(\mathbf{\Sigma^{-1}\dot{\mathbf{\Sigma}}_{ij}}\right) -\operatorname{tr}\left(\mathbf{C^{-1}W^{\prime}\Sigma^{-1}}\dot{\mathbf{\Sigma}}_{ij}\mathbf{\Sigma^{-1}W} \right)-\mathbf{e}^{\prime}\mathbf{\Sigma^{-1}} \dot{\mathbf{\Sigma}}_{ij}\mathbf{\Sigma^{-1}}\mathbf{e} \Big) \] 易得 \(\mathbf{\dot{\Sigma}_{i j}} = \mathbf{\dot{\Sigma}_{k_{i j}}} \otimes \mathbf{I}_{n}\) ，其中 \(\mathbf{\dot{\Sigma}_{k_{i j}}} = \frac{\partial \mathbf{\Sigma}_{k}}{\partial \phi_{ij}}\) 。

因此，我们推导得到 $$ \[\begin{aligned} & \mathbf{\Sigma^{-1}\dot{\Sigma}_{i j}} \\ &= \left( \mathbf{\Sigma}_{k}^{-1} \otimes \mathbf{I}^{-1} \right) \left(\mathbf{\dot{\Sigma}}_{k_{i j}} \otimes \mathbf{I}\right) \\ &= \left( \mathbf{\Sigma}_{k}^{-1} \mathbf{\dot{\Sigma}}_{k_{i j}} \otimes \mathbf{I}_{n} \right) &\\ &\\ & \operatorname{tr} (\mathbf{\Sigma^{-1}\dot{\Sigma}_{i j}}) \\ &= \operatorname{tr}\left( \mathbf{\Sigma}_{k}^{-1} \mathbf{\dot{\Sigma}}_{k_{i j}} \otimes \mathbf{I}_{n} \right) \\ &= n \operatorname{tr}\left( \mathbf{\Sigma}_{k}^{-1} \mathbf{\dot{\Sigma}}_{k_{i j}} \right) \\ &\\ &\\ &\mathbf{\Sigma}^{-1} \dot{\mathbf{\Sigma}}_{i j} \mathbf{\Sigma}^{-1} \\ &= \left( \mathbf{\Sigma}_{k}^{-1} \otimes \mathbf{I}^{-1} \right) \left(\mathbf{\dot{\Sigma}}_{k_{i j}} \otimes \mathbf{I}\right) \left( \mathbf{\Sigma}_{k}^{-1} \otimes \mathbf{I}^{-1} \right)\\ &= \left( \mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \otimes \mathbf{I} \right) &\\ &\\ &\mathbf{e}^{\prime} \mathbf{\Sigma}^{-1} \dot{\mathbf{\Sigma}}_{i j} \mathbf{\Sigma}^{-1} \mathbf{e} \\ &= \mathrm{vec}(\mathbf{E})^{\prime} \left( \mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \otimes \mathbf{I} \right) \mathrm{vec}(\mathbf{E}) \\ &= \mathrm{tr}(\mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \mathbf{E}^{\prime} \mathbf{E}) \quad \because \mathrm{vec}(\mathbf{A})^{\prime} (\mathbf{B} \otimes \mathbf{C}) \mathrm{vec}(\mathbf{A}) = \mathrm{tr}(\mathbf{B^{\prime}A^{\prime}CA}) \\ &= \mathrm{tr}(\mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \mathbf{\Delta} ) \\ \end{aligned}\]

$$ 最后一个式子中，我们设 \(\mathbf{\Delta} = \mathbf{E}^{\prime} \mathbf{E}\) 。

代入 \(U(\phi_{ij})\) 得到 \[ U(\phi_{ij}) = - \frac{1}{2} \Big( n \operatorname{tr}\left( \mathbf{\Sigma}_{k}^{-1} \mathbf{\dot{\Sigma}}_{k_{i j}} \right) -\operatorname{tr}\left(\mathbf{C^{-1}W^{\prime}\left( \mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \otimes \mathbf{I} \right) W} \right)-\mathrm{tr}(\mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \mathbf{\Delta} ) \Big) \] 我们需要对这个式子进一步简化，我们知道 \(\mathbf{\Sigma}_{k}\) 是一个对称矩阵，因此 \(\mathbf{\Sigma}_{k}^{-1}\) 也是一个对称矩阵， 我们将 \(\mathbf{\Sigma}_{k}^{-1}\) 中的元素表示为 \(\phi^{ij}\) (\(i \leq j\) ) 。

首先我们考虑第一项 \(\operatorname{tr}\left( \mathbf{\Sigma}_{k}^{-1} \mathbf{\dot{\Sigma}}_{k_{i j}} \right)\) ，同上，我们可以得到 \[ \operatorname{tr}\left(\boldsymbol{\Sigma}_k^{-1} \dot{\boldsymbol{\Sigma}}_{k_{i j}}\right)=\left\{\begin{array}{cc} \phi^{i i}, & i=j \\ 2 \phi^{i j}, & i \neq j \end{array}\right. \] 我们现在考虑第三项，如下 \[ \begin{aligned} \mathrm{tr}(\mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \mathbf{\Delta}) & = \mathrm{tr}(\mathbf{\Sigma}_{k}^{-1} \mathbf{\Delta}\mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}) \\ & = \mathrm{tr}(\mathbf{\Delta}^{*}\mathbf{\dot{\Sigma}}_{k_{i j}}) \\ \end{aligned} \] 其中我们设 \(\mathbf{\Delta}^{*} = \mathbf{\Sigma}_{k}^{-1} \mathbf{\Delta}\mathbf{\Sigma}_{k}^{-1} = \{ \delta_{ij}^{*}\}\) 。

当 \(i=j\) 时，矩阵 \(\mathbf{\Delta}^{*}\mathbf{\dot{\Sigma}}_{k_{i j}}\) 只有第 \(i\) 列元素为 \(\mathbf{\Delta}^{*}\) 的第 \(i\) 列元素 ，其它位置元素为 0 。因此，我们得到 \[ \mathrm{tr}(\mathbf{\Delta}^{*}\mathbf{\dot{\Sigma}}_{k_{i j}}) = \delta_{ii}^{*} \] 当 \(i \neq j\) 时, 矩阵 \(\mathbf{\Delta}^{*}\mathbf{\dot{\Sigma}}_{k_{i j}}\) 其第 \(i\) 列为 \(\mathbf{\Delta}^{*}\) 的第 \(j\) 列， 其第 \(j\) 列为 \(\mathbf{\Delta}^{*}\) 的第 \(i\) 列，其它位置元素为 0 。因此，我们得到 (因为 \(\mathbf{\Delta}^{*}\) 对称) \[ \mathrm{tr}(\mathbf{\Delta}^{*}\mathbf{\dot{\Sigma}}_{k_{i j}}) = \delta_{ij}^{*} + \delta_{ji}^{*} = 2 \delta_{ij}^{*} \] 综上，我们有 \[ \mathrm{tr}(\mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \mathbf{\Delta})=\left\{\begin{array}{cc} \delta^{*}_{i i}, & i=j \\ 2 \delta^{*}_{i j}, & i \neq j \end{array}\right. \] 对于第二项，首先我们对 \(\mathbf{WC^{-1}W^{\prime}}\) 进行推导如下 \[ \begin{aligned} &\mathbf{WC^{-1}W^{\prime}}\\ &= \left[\begin{array}{c} \mathbf{X} & \mathbf{Z} \end{array}\right] \left[\begin{array}{c} \mathbf{C}^{XX} & \mathbf{C}^{XZ} \\ \mathbf{C}^{ZX} & \mathbf{C}^{ZZ} \\ \end{array}\right] \left[\begin{array}{c} \mathbf{X}^{\prime} \\ \mathbf{Z}^{\prime} \end{array}\right] \\ &=\mathbf{X}\mathbf{C}^{XX}\mathbf{X}^{\prime}+\mathbf{X}\mathbf{C}^{XZ}\mathbf{Z}^{\prime}+\mathbf{Z}\mathbf{C}^{ZX}\mathbf{X}^{\prime}+\mathbf{Z}\mathbf{C}^{ZZ}\mathbf{Z}^{\prime} \\ \end{aligned} \] 因此，第二项推导得到 $$ \[\begin{aligned} &\operatorname{tr}\left(\mathbf{C^{-1}W^{\prime}\left( \mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \otimes \mathbf{I} \right) W} \right) \\ &= \operatorname{tr}\left(\mathbf{WC^{-1}W^{\prime}\left( \mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \otimes \mathbf{I} \right) } \right) \\ &= \operatorname{tr}\left(\mathbf{X}\mathbf{C}^{XX}\mathbf{X}^{\prime}\left( \mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \otimes \mathbf{I} \right) \right) + \operatorname{tr}\left(\mathbf{X}\mathbf{C}^{XZ}\mathbf{Z}^{\prime}\left( \mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \otimes \mathbf{I} \right) \right) + \operatorname{tr}\left(\mathbf{Z}\mathbf{C}^{ZX}\mathbf{X}^{\prime}\left( \mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \otimes \mathbf{I} \right) \right) + \operatorname{tr}\left(\mathbf{Z}\mathbf{C}^{ZZ}\mathbf{Z}^{\prime}\left( \mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \otimes \mathbf{I} \right) \right) \end{aligned}\] \[ 对于第二项的第一项进行推导（其余项同理可得），假设 $k=3$ ，则 $\mathbf{X}\mathbf{C}^{XX}\mathbf{X}^{\prime}$ 可以分解为 \] \[\begin{aligned} &\mathbf{X}\mathbf{C}^{XX}\mathbf{X}^{\prime} \\ &= \left[\begin{array}{cccc} \mathbf{X}_{1} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{X}_{2} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{X}_{3} \\ \end{array}\right] \left[\begin{array}{cccc} \mathbf{C}^{XX}_{11} & \mathbf{C}^{XX}_{12} & \mathbf{C}^{XX}_{13} \\ \mathbf{C}^{XX}_{21} & \mathbf{C}^{XX}_{22} & \mathbf{C}^{XX}_{23} \\ \mathbf{C}^{XX}_{31} & \mathbf{C}^{XX}_{32} & \mathbf{C}^{XX}_{33} \\ \end{array}\right] \left[\begin{array}{cccc} \mathbf{X}_{1}^{\prime} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{X}_{2}^{\prime} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{X}_{3}^{\prime} \\ \end{array}\right] \\ &= \left[\begin{array}{cccc} \mathbf{X}_{1}\mathbf{C}^{XX}_{11}\mathbf{X}_{1}^{\prime} & \mathbf{X}_{1}\mathbf{C}^{XX}_{12}\mathbf{X}_{2}^{\prime} & \mathbf{X}_{1}\mathbf{C}^{XX}_{13}\mathbf{X}_{3}^{\prime} \\ \mathbf{X}_{2}\mathbf{C}^{XX}_{21}\mathbf{X}_{1}^{\prime} & \mathbf{X}_{2}\mathbf{C}^{XX}_{22}\mathbf{X}_{2}^{\prime} & \mathbf{X}_{2}\mathbf{C}^{XX}_{23}\mathbf{X}_{3}^{\prime} \\ \mathbf{X}_{3}\mathbf{C}^{XX}_{31}\mathbf{X}_{1}^{\prime} & \mathbf{X}_{3}\mathbf{C}^{XX}_{32}\mathbf{X}_{2}^{\prime} & \mathbf{X}_{3}\mathbf{C}^{XX}_{33}\mathbf{X}_{3}^{\prime} \\ \end{array}\right] \end{aligned}\] \[ 类似与随机效应，我们可以得到（推导过程略） \] (^{XX}{}( {k}^{-1}{k_{i j}}_{k}^{-1} ) )={ \[\begin{array}{cc} \sum_{s=1}^{k}\sum_{t=1}^{k}\phi^{it}\phi^{si} \operatorname{tr} (\mathbf{X}_{s}\mathbf{C}^{XX}_{st}\mathbf{X}_{t}^{\prime} ), & i=j \\ \sum_{s=1}^{k}\sum_{t=1}^{k}\phi^{it}\phi^{sj} \left( \operatorname{tr} \left(\mathbf{X}_{s}\mathbf{C}^{XX}_{st}\mathbf{X}_{t}^{\prime}\right) + \operatorname{tr} \left(\mathbf{X}_{t}\mathbf{C}^{XX}_{ts}\mathbf{X}_{s}^{\prime}\right)\right) , & i \neq j \end{array}\] . \[ 因此，我们得到（同上，我们有 $\mathbf{C}^{AB}_{st} = (\mathbf{C}^{BA}_{ts})^{T}$ ） \] \[\begin{aligned} &\operatorname{tr}\left(\mathbf{C^{-1}W^{\prime}\left( \mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \otimes \mathbf{I} \right) W} \right)\\ &=\left\{\begin{array}{cc} \sum_{s=1}^{k}\sum_{t=1}^{k}\phi^{it}\phi^{si} \operatorname{tr} (\mathbf{X}_{s}\mathbf{C}^{XX}_{st}\mathbf{X}_{t}^{\prime} + \mathbf{X}_{s}\mathbf{C}^{XZ}_{st}\mathbf{Z}_{t}^{\prime} + \mathbf{Z}_{s}\mathbf{C}^{ZX}_{st}\mathbf{X}_{t}^{\prime} + \mathbf{Z}_{s}\mathbf{C}^{ZZ}_{st}\mathbf{Z}_{t}^{\prime} ), & i=j \\ \sum_{s=1}^{k}\sum_{t=1}^{k}\phi^{it}\phi^{sj} \left( \operatorname{tr} \left(\mathbf{X}_{s}\mathbf{C}^{XX}_{st}\mathbf{X}_{t}^{\prime} + \mathbf{X}_{s}\mathbf{C}^{XZ}_{st}\mathbf{Z}_{t}^{\prime} +\mathbf{Z}_{s}\mathbf{C}^{ZX}_{st}\mathbf{X}_{t}^{\prime} +\mathbf{Z}_{s}\mathbf{C}^{ZZ}_{st}\mathbf{Z}_{t}^{\prime} \right) + \operatorname{tr} \left(\mathbf{X}_{t}\mathbf{C}^{XX}_{ts}\mathbf{X}_{s}^{\prime}+\mathbf{X}_{t}\mathbf{C}^{XZ}_{ts}\mathbf{Z}_{s}^{\prime}+\mathbf{Z}_{t}\mathbf{C}^{ZX}_{ts}\mathbf{X}_{s}^{\prime}+\mathbf{Z}_{t}\mathbf{C}^{ZZ}_{ts}\mathbf{Z}_{s}^{\prime}\right) \right), & i \neq j \end{array}\right. \\ &=\left\{\begin{array}{cc} \sum_{s=1}^{k}\sum_{t=1}^{k}\phi^{it}\phi^{si} \operatorname{tr} (\mathbf{X}_{s}\mathbf{C}^{XX}_{st}\mathbf{X}_{t}^{\prime} + \mathbf{X}_{s}\mathbf{C}^{XZ}_{st}\mathbf{Z}_{t}^{\prime} + \mathbf{Z}_{s}\mathbf{C}^{ZX}_{st}\mathbf{X}_{t}^{\prime} + \mathbf{Z}_{s}\mathbf{C}^{ZZ}_{st}\mathbf{Z}_{t}^{\prime} ), & i=j \\ \sum_{s=1}^{k}\sum_{t=1}^{k}\phi^{it}\phi^{sj} 2 \operatorname{tr} \left(\mathbf{X}_{s}\mathbf{C}^{XX}_{st}\mathbf{X}_{t}^{\prime} + \mathbf{X}_{s}\mathbf{C}^{XZ}_{st}\mathbf{Z}_{t}^{\prime} +\mathbf{Z}_{s}\mathbf{C}^{ZX}_{st}\mathbf{X}_{t}^{\prime} +\mathbf{Z}_{s}\mathbf{C}^{ZZ}_{st}\mathbf{Z}_{t}^{\prime} \right), & i \neq j \end{array}\right. \\ \end{aligned}\] \[ 现在我们定义 $\mathbf{\Psi}$ 是一个 $k \times k$ 的矩阵，其元素为 $\psi_{ij} = \operatorname{tr} \left(\mathbf{X}_{i}\mathbf{C}^{XX}_{ij}\mathbf{X}_{j}^{\prime} + \mathbf{X}_{i}\mathbf{C}^{XZ}_{ij}\mathbf{Z}_{j}^{\prime} +\mathbf{Z}_{i}\mathbf{C}^{ZX}_{ij}\mathbf{X}_{j}^{\prime} +\mathbf{Z}_{i}\mathbf{C}^{ZZ}_{ij}\mathbf{Z}_{j}^{\prime} \right)$，根据矩阵乘法性质（假设 $\mathbf{A,B,C}$ 均为 $n \times n$ 的矩阵，则 $\mathbf{ABC} = \{\sum_{s=1}^{n}\sum_{t=1}^{n}a_{it}b_{ts}c_{sj}\}$ ），我们发现 \] \[\begin{aligned} \mathbf{\Sigma}_{k}^{-1} \mathbf{\Psi}\mathbf{\Sigma}_{k}^{-1} &= \left \{\sum_{s=1}^{k}\sum_{t=1}^{k}\phi^{it}\operatorname{tr} \left(\mathbf{X}_{t}\mathbf{C}^{XX}_{ts}\mathbf{X}_{s}^{\prime} + \mathbf{X}_{t}\mathbf{C}^{XZ}_{ts}\mathbf{Z}_{s}^{\prime} +\mathbf{Z}_{t}\mathbf{C}^{ZX}_{ts}\mathbf{X}_{s}^{\prime} +\mathbf{Z}_{t}\mathbf{C}^{ZZ}_{ts}\mathbf{Z}_{s}^{\prime} \right)\phi^{sj} \right\}\\ &= \left \{\sum_{s=1}^{k}\sum_{t=1}^{k}\phi^{it}\phi^{sj}\operatorname{tr} \left(\mathbf{X}_{t}\mathbf{C}^{XX}_{ts}\mathbf{X}_{s}^{\prime} + \mathbf{X}_{t}\mathbf{C}^{XZ}_{ts}\mathbf{Z}_{s}^{\prime} +\mathbf{Z}_{t}\mathbf{C}^{ZX}_{ts}\mathbf{X}_{s}^{\prime} +\mathbf{Z}_{t}\mathbf{C}^{ZZ}_{ts}\mathbf{Z}_{s}^{\prime} \right)\right\}\\ &= \left \{\sum_{s=1}^{k}\sum_{t=1}^{k}\phi^{it}\phi^{sj}\operatorname{tr} \left(\mathbf{X}_{s}\mathbf{C}^{XX}_{st}\mathbf{X}_{t}^{\prime} + \mathbf{X}_{s}\mathbf{C}^{XZ}_{st}\mathbf{Z}_{t}^{\prime} +\mathbf{Z}_{s}\mathbf{C}^{ZX}_{st}\mathbf{X}_{t}^{\prime} +\mathbf{Z}_{s}\mathbf{C}^{ZZ}_{st}\mathbf{Z}_{t}^{\prime} \right)\right\}\\ \end{aligned}\]

$$ 我们定义 \(\psi_{ij}^{*}\) 是 \(\mathbf{\Psi}^{*} = \mathbf{\Sigma}_{k}^{-1} \mathbf{\Psi}\mathbf{\Sigma}_{k}^{-1}\) 的元素，即 \(\psi_{ij}^{*} = \left\{ \mathbf{\Sigma}_{k}^{-1} \mathbf{\Psi}\mathbf{\Sigma}_{k}^{-1} \right\}_{ij}\) 。

因此，我们得到 \[ \operatorname{tr}\left(\mathbf{C^{-1}W^{\prime}\left( \mathbf{\Sigma}_{k}^{-1}\mathbf{\dot{\Sigma}}_{k_{i j}}\mathbf{\Sigma}_{k}^{-1} \otimes \mathbf{I} \right) W} \right)=\left\{\begin{array}{cc} \psi_{ii}^{*}, & i=j \\ 2\psi_{ij}^{*}, & i \neq j \end{array}\right. \] 将上面的式子代入 \(U(\phi_{ij})\) ，进一步简化得到 \[ U(\phi_{ij}) =\left\{\begin{array}{cc} - \frac{1}{2} \Big(n\phi^{ii} - \psi_{ii}^{*}- \delta_{ii}^{*} \big), & i=j \\ - \frac{1}{2} \Big(2n\phi^{ij} - 2\psi_{ij}^{*}- 2\delta_{ij}^{*} \big), & i \neq j \end{array}\right. \] 将之设为0， 对 \(\gamma^{ij}\) 求解得到 EM 推导公式 \[ n\phi^{ij} = \psi_{ij}^{*} + \delta_{ij}^{*} \] 使用矩阵格式，我们得到 \[ \begin{aligned} n\mathbf{\Sigma}_{k}^{-1} &= \mathbf{\Psi}^{*} + \mathbf{\Delta}^{*} \\ &= \mathbf{\Sigma}_{k}^{-1} \mathbf{\Psi}\mathbf{\Sigma}_{k}^{-1} + \mathbf{\Sigma}_{k}^{-1} \mathbf{\Delta}\mathbf{\Sigma}_{k}^{-1} \end{aligned} \] 左右均乘以 \(\mathbf{\Sigma}_{k}\) ，得到 \[ n\mathbf{\Sigma}_{k}= \mathbf{\Psi} + \mathbf{\Delta} \] 因此，我们得到 \[ n\phi_{ij} = \psi_{ij} + \delta_{ij} \] 因为 \(\psi_{ij} = \operatorname{tr} \left(\mathbf{X}_{i}\mathbf{C}^{XX}_{ij}\mathbf{X}_{j}^{\prime} + \mathbf{X}_{i}\mathbf{C}^{XZ}_{ij}\mathbf{Z}_{j}^{\prime} +\mathbf{Z}_{i}\mathbf{C}^{ZX}_{ij}\mathbf{X}_{j}^{\prime} +\mathbf{Z}_{i}\mathbf{C}^{ZZ}_{ij}\mathbf{Z}_{j}^{\prime} \right)\) ， \(\mathbf{\Delta} = \mathbf{E}^{\prime} \mathbf{E}\) 可以进一步分块为： \[ \mathbf{\Delta} = \mathbf{E}^{\prime} \mathbf{E} =\left[\begin{array}{cccc} \mathbf{e}_{1}^{\prime}\mathbf{e}_{1} & \mathbf{e}_{1}^{\prime}\mathbf{e}_{2} & \cdots & \mathbf{e}_{1}^{\prime}\mathbf{e}_{k} \\ \mathbf{e}_{2}^{\prime}\mathbf{e}_{1} & \mathbf{e}_{2}^{\prime}\mathbf{e}_{2} & \cdots & \mathbf{e}_{2}^{\prime}\mathbf{e}_{k} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{e}_{k}^{\prime}\mathbf{e}_{1} & \mathbf{e}_{k}^{\prime}\mathbf{e}_{2} & \cdots & \mathbf{e}_{k}^{\prime}\mathbf{e}_{k} \\ \end{array}\right] \] 因此 \(\delta_{ij} = \mathbf{e}_{i}^{\prime}\mathbf{e}_{j}\) 。因此上面的 EM 推导公式可以写为 \[ \phi_{ij} = \frac{1}{n} \left(\operatorname{tr} \left(\mathbf{X}_{i}\mathbf{C}^{XX}_{ij}\mathbf{X}_{j}^{\prime} + \mathbf{X}_{i}\mathbf{C}^{XZ}_{ij}\mathbf{Z}_{j}^{\prime} +\mathbf{Z}_{i}\mathbf{C}^{ZX}_{ij}\mathbf{X}_{j}^{\prime} +\mathbf{Z}_{i}\mathbf{C}^{ZZ}_{ij}\mathbf{Z}_{j}^{\prime} \right) + \mathbf{e}_{i}^{\prime}\mathbf{e}_{j} \right) \] 如果我们设 \(\mathbf{W}_{i} = \left[\begin{array}{c} \mathbf{X}_{i} & \mathbf{Z}_{i} \end{array}\right]\) , \(\mathbf{W}_{j} = \left[\begin{array}{c} \mathbf{X}_{j} & \mathbf{Z}_{j} \end{array}\right]\) , \(\mathbf{C}^{ij} = \left[\begin{array}{c} \mathbf{C}^{XX}_{ij} & \mathbf{C}^{XZ}_{ij} \\ \mathbf{C}^{ZX}_{ij} & \mathbf{C}^{ZZ}_{ij} \\ \end{array}\right]\) ，则我们有 \[ \psi_{ij} = \operatorname{tr} \left(\mathbf{X}_{i}\mathbf{C}^{XX}_{ij}\mathbf{X}_{j}^{\prime} + \mathbf{X}_{i}\mathbf{C}^{XZ}_{ij}\mathbf{Z}_{j}^{\prime} +\mathbf{Z}_{i}\mathbf{C}^{ZX}_{ij}\mathbf{X}_{j}^{\prime} +\mathbf{Z}_{i}\mathbf{C}^{ZZ}_{ij}\mathbf{Z}_{j}^{\prime} \right) = \operatorname{tr} \left( \mathbf{W}_{i} \mathbf{C}^{ij} \mathbf{W}_{j}^{\prime} \right) \] 则 \(\phi_{ij}\) 的 EM 推导公式可以写为 \[ \phi_{ij} = \frac{1}{n} \left(\operatorname{tr} \left( \mathbf{W}_{i} \mathbf{C}^{ij} \mathbf{W}_{j}^{\prime} \right) + \mathbf{e}_{i}^{\prime}\mathbf{e}_{j} \right) \]

附录1: V逆相关推导

根据舒尔补公式 \(\mathbf{(D+CA^{-1}B)^{-1}=D^{-1}-D^{-1}C(A+BD^{-1}C)^{-1}BD^{-1}}\) ， 我们可以得到 \[ \mathbf{V^{-1}} = \mathbf{(R+ZGZ^{\prime})^{-1}} = \mathbf{R^{-1}-R^{-1}Z(Z^{\prime}R^{-1}Z+G^{-1})^{-1}Z^{\prime}R^{-1}} \] 两边左乘 \(\mathbf{GZ^{\prime}}\) 得到 \[ \begin{aligned} & \quad \mathbf{GZ^{\prime}(R+ZGZ^{\prime})^{-1}} \\ &= \mathbf{GZ^{\prime}R^{-1}-GZ^{\prime}R^{-1}Z(Z^{\prime}R^{-1}Z+G^{-1})^{-1}Z^{\prime}R^{-1}}\\ &= \mathbf{GZ^{\prime}R^{-1}-G(Z^{\prime}R^{-1}Z+G^{-1})(Z^{\prime}R^{-1}Z+G^{-1})^{-1}Z^{\prime}R^{-1}+(Z^{\prime}R^{-1}Z+G^{-1})^{-1}Z^{\prime}R^{-1}} \\ &= \mathbf{GZ^{\prime}R^{-1}-GZ^{\prime}R^{-1}+(Z^{\prime}R^{-1}Z+G^{-1})^{-1}Z^{\prime}R^{-1}} \\ &= \mathbf{(Z^{\prime}R^{-1}Z+G^{-1})^{-1}Z^{\prime}R^{-1}} \\ \end{aligned} \] 因此我们得到 \[ \mathbf{(Z^{\prime}R^{-1}Z+G^{-1})^{-1}Z^{\prime}R^{-1}} = \mathbf{GZ^{\prime}(R+ZGZ^{\prime})^{-1}} \] 转置一下，得到 \[ \mathbf{R^{-1}Z(Z^{\prime}R^{-1}Z+G^{-1})^{-1}} = \mathbf{(R+ZGZ^{\prime})^{-1}ZG} \]

附录2：P 矩阵对方差组分的导数

我们有 \[ \frac{\partial \boldsymbol{P}}{\partial \kappa_i}=-\boldsymbol{P} \boldsymbol{V}_i \boldsymbol{P} \] 证明如下： \[ \begin{aligned} \frac{\partial \boldsymbol{P}}{\partial \kappa_i}= & \frac{\partial}{\partial \kappa_i}\left(\boldsymbol{V}^{-1}-\boldsymbol{V}^{-1} \boldsymbol{X}\left(\boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\prime} \boldsymbol{V}^{-1}\right) \\ = & \frac{\partial \boldsymbol{V}^{-1}}{\partial \kappa_i}-\frac{\partial \boldsymbol{V}^{-1}}{\partial \kappa_i} \boldsymbol{X}\left(\boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\prime} \boldsymbol{V}^{-1}-\boldsymbol{V}^{-1} \boldsymbol{X} \frac{\partial\left(\boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{X}\right)^{-1}}{\partial \kappa_i} \boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \\ & \quad-\boldsymbol{V}^{-1} \boldsymbol{X}\left(\boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\prime} \frac{\partial \boldsymbol{V}^{-1}}{\partial \kappa_i} \\ = & -\boldsymbol{V}^{-1} \boldsymbol{V}_i \boldsymbol{V}^{-1}+\boldsymbol{V}^{-1} \boldsymbol{V}_i \boldsymbol{V}^{-1} \boldsymbol{X}\left(\boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \\ & \quad-\boldsymbol{V}^{-1} \boldsymbol{X}\left(\boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{V}_i \boldsymbol{V}^{-1} \boldsymbol{X}\left(\boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X} \boldsymbol{V}^{-1} \\ & \quad+\boldsymbol{V}^{-1} \boldsymbol{X}\left(\boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{V}_i \boldsymbol{V}^{-1} \\ = & -\boldsymbol{V}^{-1} \boldsymbol{V}_i\left(\boldsymbol{V}^{-1}-\boldsymbol{V}^{-1} \boldsymbol{X}\left(\boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{X}^{-1} \boldsymbol{X}^{\prime} \boldsymbol{V}^{-1}\right)\right. \\ & \quad+\boldsymbol{V}^{-1} \boldsymbol{X}\left(\boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{V}_i\left(\boldsymbol{V}^{-1}-\boldsymbol{V}^{-1} \boldsymbol{X}\left(\boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\prime} \boldsymbol{V}^{-1}\right) \\ = & -\boldsymbol{V}^{-1} \boldsymbol{V}_i \boldsymbol{P}+\boldsymbol{V}^{-1} \boldsymbol{X}^{\prime}\left(\boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{V}_i \boldsymbol{P} \\ = & -\left(\boldsymbol{V}^{-1}-\boldsymbol{V}^{-1} \boldsymbol{X}\left(\boldsymbol{X}^{\prime} \boldsymbol{V}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\prime} \boldsymbol{V}^{-1}\right) \boldsymbol{V}_i \boldsymbol{P} \\ = & -\boldsymbol{P} \boldsymbol{V}_i \boldsymbol{P} \end{aligned} \]

参考文献

1. Knight E. Improved iterative schemes for REML estimation of variance parameters in linear mixed models[D]. , 2008.