这是方差组分估计方法的第二篇博客，介绍ML和REML两种方法及一些很早之前的算法。

ML

在最大似然方法中，我们需要新增分布假设，一般我们都采用正态假设，因此模型假设为 (这里假设所有随机效应的协方差矩阵均为对角矩阵) \[ \mathbf{y} \text { is } N_{n}(\mathbf{X} \boldsymbol{\beta}, \mathbf{V}), \quad \text { where } \quad \mathbf{V}=\sum_{i=1}^{m} \sigma_{i}^{2} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime}+\sigma^{2} \mathbf{I}_{n} \] 其中 \(\mathbf{X}\) 是一个 \(n \times p\) 并且秩为 \(r \leq p\) 的矩阵。\(\boldsymbol{\mathbf{V}}\) 是一个 \(n \times n\) 的正定矩阵。为了方便书写，我们记 \(\sigma_{0}^{2}=\sigma^{2}\) ， \(\mathbf{Z}_{0}=\mathbf{I}_{n}\) ，因此协方差矩阵变成 (假设随机向量之间互不相关) \[ \mathbf{V}=\sum_{i=0}^{m} \sigma_{i}^{2} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \]

随机向量 \(\mathbf{y}\) 的联合分布密度函数，即似然函数为 \[ L(\mathbf{y})=\frac{1}{(2 \pi)^{n / 2}|\mathbf{V}|^{1 / 2}} \exp \left\{-0.5(\mathbf{y}-\mathbf{X} \beta)^{\prime} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X} \beta)\right\} \] 其对数为 \[ l=-0.5 n \times \ln (2 \pi)-0.5 \ln |\mathbf{V}|-0.5(\mathbf{y}-\mathbf{X} \beta)^{\prime} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X} \beta) \] 对他求未知参数 \(\boldsymbol{\boldsymbol{\theta}}^{\prime}=\left[\begin{array}{lllll}\boldsymbol{\beta}^{\prime} & \sigma_{0}^{2} & \sigma_{1}^{2} & \cdots & \sigma_{m}^{2}\end{array}\right]'\) 的偏导数如下（证明略） \[ \begin{aligned} &\frac{\partial l}{\partial \boldsymbol{\beta}}=-\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X} \boldsymbol{\beta}-\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{y}\right) \\ &\frac{\partial l}{\partial \sigma_{i}^{2}}=-\frac{1}{2} \operatorname{tr}\left(\mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \sigma_{i}^{2}}\right)+\frac{1}{2}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta})^{\prime} \mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \sigma_{i}^{2}} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta}) \\ &\left(\because \frac{\partial}{\partial \sigma_{i}^{2}} \ln |\mathbf{V}|=\operatorname{tr}\left(\mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \sigma_{i}^{2}}\right), \quad \frac{\partial}{\partial \sigma_{i}^{2}} \mathbf{V}^{-1}=-\mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \sigma_{i}^{2}} \mathbf{V}^{-1}\right) \quad \text{参考 linear models in statistics } \end{aligned} \] 令上面两式为0，我们设 \(\mathbf{V}_{i}= \frac{\partial \mathbf{V}}{\partial \sigma_{i}^{2}}=\mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime}\) ，因此我们有 \[ \begin{aligned} \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X} \hat{\boldsymbol{\beta}} &=\mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{y} \\ \operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{V}_{i}\right) &=(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{V}_{i} \hat{\mathbf{V}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}}) \\ &=\mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{V}_{i} \hat{\mathbf{P}} {\mathbf{y}} \end{aligned} \] 其中 \(\hat{\mathbf{P}}=\hat{\mathbf{V}}^{-1}-\hat{\mathbf{V}}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1}\) ，因为下式成立 \[ \begin{aligned} \hat{\mathbf{P}} \mathbf{y} &=\left(\hat{\mathbf{V}}^{-1}-\hat{\mathbf{V}}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1}\right) \mathbf{y} \\ &=\hat{\mathbf{V}}^{-1}\left(\mathbf{y}-\mathbf{X}\left(\mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{y}\right) \\ &=\hat{\mathbf{V}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}}) \end{aligned} \] 因为
\[ \operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{V}_{i}\right)=\operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{V}_{i} \hat{\mathbf{V}}^{-1} \hat{\mathbf{V}}\right)=\sum_{j} \operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{V}_{i} \hat{\mathbf{V}}^{-1} \mathbf{V}_{j} \hat{\sigma}_{j}^{2}\right) \] 故有 ML 的估计方程为 \[ \begin{aligned} &\sum_{j} \left\{\operatorname{tr} \hat{\mathbf{V}}^{-1} \mathbf{V}_{i} \hat{\mathbf{V}}^{-1} \mathbf{V}_{j}\right\}\left\{\hat{\sigma}_{j}^{2}\right\}=\left\{\mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{V}_{i} \hat{\mathbf{P}} \mathbf{y}\right\} \\ &i, j=0,1, \cdots, m \end{aligned} \] 这个式子还可以换一种形式，对于左手项： $$ \[\begin{aligned} & \sum_{j} \left\{\operatorname{tr} \left( \hat{\mathbf{V}}^{-1} \mathbf{V}_{i} \hat{\mathbf{V}}^{-1} \mathbf{V}_{j} \right) \right\}\left\{\hat{\sigma}_{j}^{2}\right\} \\ &=\sum_{j} \left\{\operatorname{tr} \left( \hat{\mathbf{V}}^{-1} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{Z}_{j}\mathbf{Z}_{j}^{\prime} \right)\right\}\left\{\hat{\sigma}_{j}^{2}\right\} \\ & = \sum_{j} \left\{\operatorname{tr} \left( \mathbf{Z}_{j}^{\prime} (\hat{\mathbf{V}}^{-1} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{Z}_{j})\right)\right\}\left\{\hat{\sigma}_{j}^{2}\right\} \\ & = \sum_{j} \left\{\operatorname{tr}\left( (\mathbf{Z}_{i}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{Z}_{j})^{\prime} (\mathbf{Z}_{i}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{Z}_{j})\right)\right\}\left\{\hat{\sigma}_{j}^{2}\right\} \\ & = \sum_{j} \left\{ \mathrm{sesq} (\mathbf{Z}_{i}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{Z}_{j})\right\}\left\{\hat{\sigma}_{j}^{2}\right\} \\ \end{aligned}\]

$$ 这里的 \(\mathrm{sesq}(\mathbf{A})\) 就是矩阵 \(\mathbf{A}\) 所有元素的平方和。

同理右手项可以写为 \[ \mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{V}_{i} \hat{\mathbf{P}} \mathbf{y} = \mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \hat{\mathbf{P}} \mathbf{y} = \mathrm{sesq}(\mathbf{Z}_{i}^{\prime} \hat{\mathbf{P}} \mathbf{y}) \] 因此，上面的式子可以写为 \[ \begin{aligned} &\sum_{j} \left\{ \mathrm{sesq} (\mathbf{Z}_{i}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{Z}_{j})\right\}\left\{\hat{\sigma}_{j}^{2}\right\}=\mathrm{sesq}(\mathbf{Z}_{i}^{\prime} \hat{\mathbf{P}} \mathbf{y}) \\ &i, j=0,1, \cdots, m \end{aligned} \] 很显然这个式子不是 \(\boldsymbol{\sigma}^{2}\) 的线性方程组，左手项的 \(\hat{\mathbf{V}}^{-1}\) 和右手项的 \(\hat{\mathbf{P}}\) 均包含\(\boldsymbol{\sigma}^{2}\) 。根据此式子，我们显然找不到闭式解，还需要采用迭代的方式求解。我们可以设置 \(\boldsymbol{\hat{\sigma}}^{2}\) 的一个初始值，然后得到 \(\hat{\mathbf{V}}^{-1}\) 和 \(\hat{\mathbf{P}}\) ，此时上面的方程组就变成了线性方程组，就可以求解得到一组新的 \(\boldsymbol{\hat{\sigma}}^{2}\) 。

似然法的一个缺点式，当数据量比较小时，其估计值时有偏的。而且随着固定效应水平数目的增加和随机效应水平数目的减少，其偏差会愈加严重。

估计值的方差

由于 ML 估计值必须通过迭代求解得到，因此不可能求得估计值的准确的抽样方差，但当样本数目较大时，可以求得估计值的近似方差（称为大样本方差）。这个方差是借助参数的信息矩阵来求得的，参数 \(\boldsymbol{\theta}\) 的信息矩阵用 \(\mathbf{I}(\boldsymbol{\theta})\) 表示，其定义为 \[ \mathbf{I}(\boldsymbol{\theta}) = - \operatorname{E}\left\{ \frac{\partial^{2} l}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^{\prime} } \right\} \] 其为二阶导的数学期望的负值。

Searle (1970) 证明，当样本数目足够大时，参数 \(\boldsymbol{\theta}\) 的 ML 估计值 \(\boldsymbol{\hat{\theta}}\) 的方差近似为： \[ \operatorname{var} (\boldsymbol{\hat{\theta}}) \approx \left[\mathbf{I}(\boldsymbol{\theta}) \right]^{-1} \] 注意这式子本身并不需要 ML 估计值 \(\boldsymbol{\hat{\theta}}\) 。

ML 估计值的信息矩阵写为下式，其中 \(l_{\boldsymbol{\beta} \boldsymbol{\beta}}\) 表示 \(\partial^{2} l/\partial \boldsymbol{\beta} \partial \boldsymbol{\beta}^{\prime}\) ，其余式子类似。 \[ I\left[\begin{array}{c} \boldsymbol{\beta} \\ \boldsymbol{\sigma}^2 \end{array}\right]=-E\left[\begin{array}{cc} l_{\boldsymbol{\beta} \boldsymbol{\beta}} & l_{\boldsymbol{\beta} \boldsymbol{\sigma}^2} \\ l_{\boldsymbol{\sigma}^2 \boldsymbol{\beta}} & l_{\boldsymbol{\sigma}^2 \boldsymbol{\sigma}^2} \end{array}\right] \] 因为我们有 \[ \frac{\partial l}{\partial \boldsymbol{\beta}}=-\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X} \boldsymbol{\beta}-\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{y}\right) \]

\[ \frac{\partial}{\partial \sigma_{i}^{2}} \mathbf{V}^{-1}=-\mathbf{V}^{-1} \frac{\partial \mathbf{V}}{\partial \sigma_{i}^{2}} \mathbf{V}^{-1} = -\mathbf{V}^{-1} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime} \mathbf{V}^{-1} \]

因此，推导得到 \[ l_{\boldsymbol{\beta} \boldsymbol{\beta}} = -\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X} \]

\[ l_{\boldsymbol{\beta} \boldsymbol{\sigma}_{i}^2} = -\mathbf{X}^{\prime}\mathbf{V}^{-1}\mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime}\mathbf{V}^{-1}(\mathbf{y} - \mathbf{X} \boldsymbol{\beta}) \]

同时我们有 \[ \frac{\partial l}{\partial \sigma_{i}^{2}}=-\frac{1}{2} \operatorname{tr}\left(\mathbf{V}^{-1} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime}\right)+\frac{1}{2}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta})^{\prime} \mathbf{V}^{-1} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta}) \] 因此，推导得到 \[ \begin{aligned} l_{\boldsymbol{\sigma}_{i}^2 \boldsymbol{\sigma}_{j}^2} &= \frac{1}{2}\operatorname{tr}(\mathbf{V}^{-1} \mathbf{Z}_{j}\mathbf{Z}_{j}^{\prime}\mathbf{V}^{-1} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime}) - \frac{1}{2}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta})^{\prime} \mathbf{V}^{-1} \mathbf{Z}_{j}\mathbf{Z}_{j}^{\prime} \mathbf{V}^{-1} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta}) \\ & \quad - \frac{1}{2}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta})^{\prime} \mathbf{V}^{-1} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime} \mathbf{V}^{-1} \mathbf{Z}_{j}\mathbf{Z}_{j}^{\prime} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta}) \\ &= \frac{1}{2}\operatorname{tr}(\mathbf{V}^{-1} \mathbf{Z}_{j}\mathbf{Z}_{j}^{\prime}\mathbf{V}^{-1} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime}) - (\mathbf{y}-\mathbf{X} \boldsymbol{\beta})^{\prime} \mathbf{V}^{-1} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime} \mathbf{V}^{-1} \mathbf{Z}_{j}\mathbf{Z}_{j}^{\prime} \mathbf{V}^{-1}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta}) \\ \end{aligned} \] 然后再计算期望，我们已知 \(\operatorname{E}(\mathbf{y}) = \mathbf{X} \boldsymbol{\beta}\) ，即 \(\operatorname{E}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta}) = \mathbf{0}\)。我们再利用一个公式 \(\operatorname{E}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta})^{\prime} \mathbf{T} (\mathbf{y}-\mathbf{X} \boldsymbol{\beta}) = \operatorname{tr} (\mathbf{TV})\) , 其中 \(\mathbf{T}\) 是一个非随机的矩阵。

证明过程如下，将 \(\mathbf{X} \boldsymbol{\beta}\) 改成为 \(\boldsymbol{\mu}\) $$ \[\begin{aligned} &\operatorname{E}(\mathbf{y}-\boldsymbol{\mu})^{\prime} \mathbf{T} (\mathbf{y}-\boldsymbol{\mu}) \\ &= \operatorname{E}(\mathbf{y}^{\prime}\mathbf{T}-\boldsymbol{\mu}^{\prime}\mathbf{T}) (\mathbf{y}-\boldsymbol{\mu}) \\ &= \operatorname{E}(\mathbf{y}^{\prime}\mathbf{T}\mathbf{y}-\boldsymbol{\mu}^{\prime}\mathbf{T}\mathbf{y} - \mathbf{y}^{\prime}\mathbf{T}\boldsymbol{\mu} + \boldsymbol{\mu}^{\prime}\mathbf{T} \boldsymbol{\mu}) \\ &= \operatorname{E}(\mathbf{y}^{\prime}\mathbf{T}\mathbf{y})-\operatorname{E}(\boldsymbol{\mu}^{\prime}\mathbf{T}\mathbf{y}) - \operatorname{E}(\mathbf{y}^{\prime}\mathbf{T}\boldsymbol{\mu}) + \operatorname{E}(\boldsymbol{\mu}^{\prime}\mathbf{T} \boldsymbol{\mu}) \\ &= \operatorname{E}(\mathbf{y}^{\prime}\mathbf{T}\mathbf{y})-\boldsymbol{\mu}^{\prime}\mathbf{T}\operatorname{E}(\mathbf{y}) - \operatorname{E}(\mathbf{y}^{\prime})\mathbf{T}\boldsymbol{\mu} + \boldsymbol{\mu}^{\prime}\mathbf{T} \boldsymbol{\mu} \\ &= \operatorname{E}(\mathbf{y}^{\prime}\mathbf{T}\mathbf{y})-\boldsymbol{\mu}^{\prime}\mathbf{T} \boldsymbol{\mu} \\ &= \operatorname{tr} (\mathbf{TV})+\boldsymbol{\mu}^{\prime}\mathbf{T} \boldsymbol{\mu}-\boldsymbol{\mu}^{\prime}\mathbf{T} \boldsymbol{\mu} \quad \because \operatorname{E}(\mathbf{y}^{\prime}\mathbf{T}\mathbf{y}) = \operatorname{tr} (\mathbf{TV})+\boldsymbol{\mu}^{\prime}\mathbf{T} \boldsymbol{\mu} \\ &= \operatorname{tr} (\mathbf{TV}) \\ \end{aligned}\] \[ 因此我们得到： \] \[\begin{aligned} -E l_{\boldsymbol{\beta} \boldsymbol{\beta}} & =E\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)=\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}, \\ -E l_{\boldsymbol{\beta} \sigma_i^2} & =\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{Z}_i \mathbf{Z}_i^{\prime} \mathbf{V}^{-1} E(\mathbf{y}-\mathbf{X} \boldsymbol{\beta})=\mathbf{0} \\ -E l_{\sigma_i^2 \sigma_j^2} & =-\frac{1}{2} \operatorname{tr}\left(\mathbf{V}^{-1} \mathbf{Z}_i \mathbf{Z}_i^{\prime} \mathbf{V}^{-1} \mathbf{Z}_j \mathbf{Z}_j^{\prime}\right)+\operatorname{tr}\left(\mathbf{V}^{-1} \mathbf{Z}_i \mathbf{Z}_i^{\prime} \mathbf{V}^{-1} \mathbf{Z}_j \mathbf{Z}_j^{\prime} \mathbf{V}^{-1} \mathbf{V}\right) \\ & =\frac{1}{2} \operatorname{tr}\left(\mathbf{V}^{-1} \mathbf{Z}_i \mathbf{Z}_i^{\prime} \mathbf{V}^{-1} \mathbf{Z}_j \mathbf{Z}_j^{\prime}\right) \end{aligned}\] \[ 因此信息矩阵为 \] =. \[ 因此估计值的协方差矩阵为 \] \[\begin{aligned} \operatorname{var}\left[\begin{array}{c} \boldsymbol{\beta} \\ \tilde{\sigma}^2 \end{array}\right] & \simeq\left(\mathbf{I}\left[\begin{array}{c} \boldsymbol{\beta} \\ \boldsymbol{\sigma}^2 \end{array}\right]\right)^{-1} \\ & =\left[\begin{array}{cc} \left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-1} & \mathbf{0} \\ \mathbf{0} & 2\left[\left\{\mathrm{~m} \operatorname{tr}\left(\mathbf{V}^{-1} \mathbf{Z}_i \mathbf{Z}_i^{\prime} \mathbf{V}^{-1} \mathbf{Z}_j \mathbf{Z}_j^{\prime}\right)\right\}_{i, j=0}^r\right]^{-1} \end{array}\right], \end{aligned}\]

$$

The Hartley-Rao form

Hartley 和 Rao 在 1967 年提出了基于 \(\mathbf{H}\) 矩阵的似然函数，其中 \(\mathbf{H}\) 矩阵定义为 \[ \mathbf{V} = \mathbf{H} \sigma_{e}^{2} \qquad \text{with} \qquad \mathbf{V}^{-1} = \mathbf{H}^{-1}/ \sigma_{e}^{2} \] 因此 \(\mathbf{H}\) 矩阵和 \(\mathbf{V}\) 矩阵具有相同的形式，不过 \(\mathbf{V}\) 矩阵中的 \(\sigma_{e}^{2}\) 变成了1，而 \(\sigma_{i}^{2}\) 变成了 \(\gamma_{i} = \sigma_{i}^{2}/\sigma_{e}^{2}\) ，这意味着此时需要估计的参数为 \(\boldsymbol{\beta}\) ，\(\sigma_{e}^{2}\) 和 \(\gamma_{1}, \cdots, \gamma_{r}\) (这里的 \(r\) 就是上面的 \(m\) ，随机效应数目)，此时我们将 \(\sigma_{e}^{2}\) 单独拎了出来。

我们从之前的公式 \(\operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime}\right) =\mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \hat{\mathbf{P}} {\mathbf{y}}\) 出发进行推导，首先，当 \(i=0\) 时，此时 \(\sigma_{0}^{2} = \sigma_{e}^{2}\)，\(\mathbf{Z}_{0} = \mathbf{I}\) ，因此我们有 \[ \operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \right) =\mathbf{y}^{\prime} \hat{\mathbf{P}}^{2} {\mathbf{y}} = (\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})^{\prime} \hat{\mathbf{V}}^{-2}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}}) \] 代入 \(\mathbf{V}^{-1} = \mathbf{H}^{-1}/ \sigma_{e}^{2}\) ，得到 \[ \operatorname{tr}\left(\hat{\mathbf{H}}^{-1} \right) = (\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})^{\prime} \hat{\mathbf{H}}^{-2}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})/ \hat{\sigma}_{e}^{2} \] 因此 \[ \hat{\sigma}_{e}^{2} = \frac{(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})^{\prime} \hat{\mathbf{H}}^{-2}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})}{\operatorname{tr}\left(\hat{\mathbf{H}}^{-1} \right)} \] 然后我们继续看其它方差组分 (\(i=1,\cdots,r\)) 的式子，对于公式 \(\operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime}\right) =\mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \hat{\mathbf{P}} {\mathbf{y}}\) ，左右乘以 \(\hat{\sigma}_{i}^{2}\) 并对 \(i=1,\cdots,r\) 的式子求和得到： \[ \operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \sum_{i=1}^{r} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \hat{\sigma}_{i}^{2} \right) =\mathbf{y}^{\prime} \hat{\mathbf{P}} \left(\sum_{i=1}^{r} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \hat{\sigma}_{i}^{2} \right) \hat{\mathbf{P}} {\mathbf{y}} \] 代入\(\mathbf{V}^{-1} = \mathbf{H}^{-1}/ \sigma_{e}^{2}\) ， \(\hat{\mathbf{P}} \mathbf{y} =\hat{\mathbf{V}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})\) 以及 \(\sum_{i=1}^{r} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \hat{\sigma}_{i}^{2} = \mathbf{V} - \sigma_{e}^{2}\mathbf{I}= (\mathbf{H} - \mathbf{I})\sigma_{e}^{2}\) ，得到 \[ \operatorname{tr}\left((\mathbf{\hat{H}}^{-1}/\hat{\sigma}_{e}^{2}) (\mathbf{\hat{H}} - \mathbf{I})\hat{\sigma}_{e}^{2} \right) =(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})^{\prime} \mathbf{\hat{H}}^{-1}/ \hat{\sigma}_{e}^{2} (\mathbf{\hat{H}} - \mathbf{I})\hat{\sigma}_{e}^{2} \mathbf{\hat{H}}^{-1}/ \hat{\sigma}_{e}^{2}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}}) \] 化简得到 \[ \operatorname{tr}\left(\mathbf{I} - \mathbf{\hat{H}}^{-1} \right) =(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})^{\prime} (\mathbf{\hat{H}}^{-1} - \mathbf{\hat{H}}^{-2}) (\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})/\hat{\sigma}_{e}^{2} \] 等价于 \[ N\hat{\sigma}_{e}^{2} = (\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})^{\prime} \mathbf{\hat{H}}^{-1} (\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}}) + \operatorname{tr}\left(\mathbf{\hat{H}}^{-1} \right) \left[ \hat{\sigma}_{e}^{2} - \frac{(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})^{\prime} \mathbf{\hat{H}}^{-2} (\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})}{\operatorname{tr}\left(\mathbf{\hat{H}}^{-1} \right)} \right] \] 根据上面的推导公式，方括号内的值为 0，因此我们得到 \(\hat{\sigma}_{e}^{2}\) 的计算公式 \[ \hat{\sigma}_{e}^{2} = \frac{(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})^{\prime} \mathbf{\hat{H}}^{-1} (\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})}{N} \] 对于公式 \(\operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime}\right) =\mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \hat{\mathbf{P}} {\mathbf{y}}\) ，化简得到 \[ \operatorname{tr}\left(\hat{\mathbf{H}}^{-1} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime}\right) =(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})^{\prime} \mathbf{\hat{H}}^{-1} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \mathbf{\hat{H}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})/\hat{\sigma}_{e}^{2} \] 汇总一下，Hartley 和 Rao 方法中设计的方程组为 $$ \[\begin{aligned} \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X} \hat{\boldsymbol{\beta}} &=\mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{y} \\ \hat{\sigma}_{e}^{2} &= \frac{(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})^{\prime} \mathbf{\hat{H}}^{-1} (\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})}{N} \\ &\text{and} \\ \operatorname{tr}\left(\hat{\mathbf{H}}^{-1} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime}\right) &=(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})^{\prime} \mathbf{\hat{H}}^{-1} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \mathbf{\hat{H}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})/\hat{\sigma}_{e}^{2} \\ &i=1, \cdots, r \end{aligned}\]

$$

EM 算法

混合详细模型如下： \[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \mathbf{Z}\mathbf{u} + \mathbf{e} \] 对于公式 \(\operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{V}_{i}\right)=\mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{V}_{i} \hat{\mathbf{P}} {\mathbf{y}}\) ，我们在等式两边同乘以 \(\hat{\sigma}^{2}_{i}\) ，并求总和，得到： \[ \begin{aligned} \operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \sum \mathbf{\hat{V}}_{i} \hat{\sigma}_{i}^{2}\right) &=\mathbf{y}^{\prime} \hat{\mathbf{P}} \sum \mathbf{V}_{i} \hat{\sigma}_{i}^{2} \hat{\mathbf{P}}\mathbf{y} \\ N &=\mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{\hat{V}} \hat{\mathbf{P}}\mathbf{y} \\ \end{aligned} \] 这里 \(N\) 为 \(\mathbf{V}\) 矩阵维度，即表型数目。

我们再看右手项，在之前 ML 的推导中，我们知道 \(\hat{\mathbf{P}} \mathbf{y} =\hat{\mathbf{V}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})\) ，因此我们有 $$ \[\begin{aligned} \mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{\hat{V}} \hat{\mathbf{P}}\mathbf{y} &=(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})^{\prime} \hat{\mathbf{V}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})\\ &=\mathbf{y}^{\prime} \hat{\mathbf{V}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}}) - \hat{\boldsymbol{\beta}}^{\prime} \left( \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{y} - \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X} \hat{\boldsymbol{\beta}} \right)\\ &=\mathbf{y}^{\prime} \hat{\mathbf{V}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}}) \quad\left(\because \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X} \hat{\boldsymbol{\beta}}=\mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{y}\right)\\ &=\mathbf{y}^{\prime}\left[\hat{\mathbf{R}}^{-1}-\hat{\mathbf{R}}^{-1} \mathbf{Z}\left(\mathbf{Z}^{\prime} \hat{\mathbf{R}}^{-1} \mathbf{Z}+\hat{\mathbf{G}}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \hat{\mathbf{R}}^{-1}\right](\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})\\ &=\mathbf{y}^{\prime} \hat{\mathbf{R}}^{-1}\left[(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})-\mathbf{Z}\left(\mathbf{Z}^{\prime} \hat{\mathbf{R}}^{-1} \mathbf{Z}+\hat{\mathbf{G}}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \hat{\mathbf{R}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})\right]\\ &=\mathbf{y}^{\prime} \hat{\mathbf{R}}^{-1}[(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})-\mathbf{Z} \hat{\mathbf{u}}]\\ &\left(\because \hat{\mathbf{u}}=\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1}\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{y}-\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{X} \hat{\boldsymbol{\beta}}\right)\right)\\ &=\frac{1}{\hat{\sigma}_{0}^{2}} \mathbf{y}^{\prime}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}}-\mathbf{Z} \hat{\mathbf{u}}) \quad \because \text{假设为单性状模型} \end{aligned}\]

\[ 因此，我们有 \] _{0}^{2}={}(- - ) /N $$ 其它方差组分公式的证明过程如下

首先我们有 \[ \operatorname{var}(\mathbf{e}) = \mathbf{R} = \sigma_{e}^{2}\mathbf{I}_{N} \quad \text{and} \quad \operatorname{var}(\mathbf{u}) = \mathbf{G} = \left\{_{d} \sigma_{i}^{2} \mathbf{I}_{q_{i}} \right\}_{i=1}^{r} \] 这里我们进一步提出两个矩阵 \(\mathbf{W}\) 和 \(\mathbf{F}_{ii}\) 矩阵，第一个矩阵如下，其中 \(q_{.} = \sum_{i=1}^{r} q_{i}\) ，即所有随机效应水平总数。 \[ \mathbf{W} = (\mathbf{I+Z'R^{-1}ZG})^{-1} = \left\{ \mathbf{W}_{ij} \right\}_{i,j=1}^{r} \quad \text{with} \quad \mathbf{W}^{-1} - \mathbf{I}_{q_{.}} = \mathbf{Z'R^{-1}ZG} \] 第二个矩阵是 \(\mathbf{G}\) 矩阵的变体， \(\mathbf{F}_{ii}\) 是将 \(\mathbf{G}\) 矩阵中的所有 \(\sigma_{i}^{2}\) 替换为 1，而其它 \(\sigma_{j}^{2}\) (\(j \neq i\) ) 替换为0 （注意这里的 \(\mathbf{G}\) 矩阵是一个对角矩阵），因此我没有 \[ \mathbf{G}\mathbf{F}_{ii} = \sigma_{i}^{2}\mathbf{F}_{ii} \quad \text{and} \quad \mathbf{F}_{ii} = \mathbf{G}\mathbf{F}_{ii}/\sigma_{i}^{2} \] 举个例子，假设 \(q_{1}=2, q_{2}=3,q_{3}=4\) ，那么我们有 \[ \mathbf{F_{22}}=\left[\begin{array}{ccc} \mathbf{0}_{2 \times 2} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{1}_{3 \times 3} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0}_{4 \times 4} \end{array}\right] \text {. } \] 根据 \(\mathbf{V}^{-1}\) 公式，我们可以推导出 \[ \begin{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} \\ &= \mathbf{R}^{-1}-\mathbf{R}^{-1} \mathbf{Z}\left(\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z} \mathbf{G} + \mathbf{I} \right) \mathbf{G}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \mathbf{R}^{-1} \\ &= \mathbf{R}^{-1}-\mathbf{R}^{-1} \mathbf{Z}\mathbf{G}\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z} \mathbf{G} + \mathbf{I} \right)^{-1} \mathbf{Z}^{\prime} \mathbf{R}^{-1} \\ &= \mathbf{R}^{-1}-\mathbf{R}^{-1} \mathbf{Z}\mathbf{G}\mathbf{W} \mathbf{Z}^{\prime} \mathbf{R}^{-1} \\ \end{aligned} \] 对于其它方程组分我们有 \(\operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime} \right)=\mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{V}_{i} \hat{\mathbf{P}} {\mathbf{y}}\) 的左手项，我们有 $$ \[\begin{aligned} \operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime} \right) &= \operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \left( \mathbf{Z} \mathbf{F}_{ii}\right) \left( \mathbf{F}_{ii} \mathbf{Z}^{\prime} \right)\right) \\ &= \operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{Z} \mathbf{F}_{ii}^{2} \mathbf{Z}^{\prime}\right) \\ &= \operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{Z} \mathbf{F}_{ii} \mathbf{Z}^{\prime}\right) \\ &= \operatorname{tr}\left(\mathbf{Z}^{\prime}\hat{\mathbf{V}}^{-1} \mathbf{Z} \mathbf{F}_{ii} \right) \\ &= \operatorname{tr}\left(\mathbf{Z}^{\prime}\left[\mathbf{\hat{R}}^{-1}-\mathbf{\hat{R}}^{-1} \mathbf{Z}\mathbf{\hat{G}}\mathbf{W} \mathbf{Z}^{\prime} \mathbf{\hat{R}}^{-1} \right] \mathbf{Z} \mathbf{F}_{ii} \right) \\ &= \operatorname{tr}\left(\left[\mathbf{Z}^{\prime} \mathbf{\hat{R}}^{-1} \mathbf{Z} \mathbf{\hat{G}} -\mathbf{Z}^{\prime} \mathbf{\hat{R}}^{-1} \mathbf{Z}\mathbf{\hat{G}}\mathbf{W} \mathbf{Z}^{\prime} \mathbf{\hat{R}}^{-1} \mathbf{Z}\mathbf{\hat{G}} \right] \mathbf{F}_{ii}/\sigma_{i}^{2} \right) \quad \because \mathbf{F}_{ii} = \mathbf{G}\mathbf{F}_{ii}/\sigma_{i}^{2} \\ &= \operatorname{tr}\left(\left[ \mathbf{W}^{-1} - \mathbf{I} -\left( \mathbf{W}^{-1} - \mathbf{I} \right)\mathbf{W} \left( \mathbf{W}^{-1} - \mathbf{I} \right) \right] \mathbf{F}_{ii}/\sigma_{i}^{2} \right) \\ &= \operatorname{tr}\left(\left( \mathbf{I} - \mathbf{W} \right) \mathbf{F}_{ii}/\sigma_{i}^{2} \right) \\ &= \operatorname{tr}\left( \mathbf{F}_{ii} - \mathbf{W} \mathbf{F}_{ii} \right) /\sigma_{i}^{2} \\ &= \frac{ \mathbf{q}_{i} - \operatorname{tr} \left(\mathbf{W}_{ii}\right)} {\sigma_{i}^{2}} \\ \end{aligned}\]

$$ 其中 \(\mathbf{W}_{ii}\) 是 \(\mathbf{W}\) 矩阵中第 \((i,i)\) 个子矩阵。

右手项可以写为 $$ \[\begin{aligned} \mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \hat{\mathbf{P}} {\mathbf{y}} & = \mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{Z} \mathbf{F}_{ii} \mathbf{F}_{ii}^{\prime} \mathbf{Z}^{\prime} \hat{\mathbf{P}} {\mathbf{y}} \\ & = \mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{Z} \hat{\mathbf{G}} \mathbf{F}_{ii} \mathbf{F}_{ii}^{\prime} \hat{\mathbf{G}}^{\prime} \mathbf{Z}^{\prime} \hat{\mathbf{P}} {\mathbf{y}} / \sigma_{i}^{4} \\ & = \mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{Z} \hat{\mathbf{G}} \mathbf{F}_{ii} \hat{\mathbf{G}} \mathbf{Z}^{\prime} \hat{\mathbf{P}} {\mathbf{y}} / \sigma_{i}^{4} \\ \end{aligned}\]

$$ 因为 \(\hat{\mathbf{G}} \mathbf{Z}^{\prime} \hat{\mathbf{P}} {\mathbf{y}} = \hat{\mathbf{G}} \mathbf{Z}^{\prime} \hat{\mathbf{V}}^{-1} \hat{\mathbf{V}} \hat{\mathbf{P}} {\mathbf{y}} = \hat{\mathbf{G}} \mathbf{Z}^{\prime} \hat{\mathbf{V}}^{-1} (\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}}) = \hat{\mathbf{u}}\)

因此右手项可以继续写为 \[ \mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \hat{\mathbf{P}} {\mathbf{y}} = \frac{\hat{\mathbf{u}}^{\prime} \mathbf{F}_{ii} \hat{\mathbf{u}}} {\sigma_{i}^{4}} = \frac{\hat{\mathbf{u}}^{\prime}_{i} \hat{\mathbf{u}}_{i}} {\sigma_{i}^{4}} \] 因此我们将 \(\operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime} \right)=\mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{V}_{i} \hat{\mathbf{P}} {\mathbf{y}}\) 改写为下式 \[ \frac{ \mathbf{q}_{i} - \operatorname{tr} \left(\mathbf{W}_{ii}\right)} {\sigma_{i}^{2}} = \frac{\hat{\mathbf{u}}^{\prime}_{i} \hat{\mathbf{u}}_{i}} {\sigma_{i}^{4}} \] 因此我们得到了两个迭代公式 \[ \hat{\sigma}_{i}^{2}=\frac{\hat{\mathbf{u}}_{i}^{\prime} \hat{\mathbf{u}}_{i}+ \hat{\sigma}_{i}^{2} \operatorname{tr} \left(\mathbf{W}_{ii}\right) } {q_{i}} \] 或者 \[ \hat{\sigma}_{i}^{2}=\frac{\hat{\mathbf{u}}_{i}^{\prime} \hat{\mathbf{u}}_{i}}{q_{i}- \operatorname{tr} \left(\mathbf{W}_{ii}\right) } \] 这里 \(q_{i}\) 是第 \(i\) 个随机效应的水平数目；\(\hat{\mathbf{u}}_{i}\) 是第 \(i\) 个随机效应的预测值。

同时使用公式 \[ \hat{\sigma}_{0}^{2}=\mathbf{y}^{\prime}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}}-\mathbf{Z} \hat{\mathbf{u}}) /N \]

额外证明

1. \(\mathbf{W} = (\mathbf{I+Z'R^{-1}ZG})^{-1}\) 矩阵存在，因为 \(\mathbf{I+Z'R^{-1}ZG}\) 的行列式等于 \(|\mathbf{R}^{-1}||\mathbf{V}| \neq 0\) 。

   不会证明。

2. \(\operatorname{tr}\left(\mathbf{W}_{ii} \right) > 0\)

   证明：

   因为 \[ \mathbf{G}^{1/2}\mathbf{W}^{-1} = \left(\mathbf{I} + \mathbf{G}^{1/2}\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z} \mathbf{G}^{1/2} \right) \mathbf{G}^{1/2} \] 其逆矩阵为 \[ \mathbf{W}\mathbf{G}^{-1/2} =\mathbf{G}^{-1/2} \left(\mathbf{I} + \mathbf{G}^{1/2}\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z} \mathbf{G}^{1/2} \right)^{-1} \] 左右均乘以 \(\mathbf{G}^{1/2}\) 得到 \[ \mathbf{G}^{1/2}\mathbf{W} = \left(\mathbf{I} + \mathbf{G}^{1/2}\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z} \mathbf{G}^{1/2} \right)^{-1} \mathbf{G}^{1/2} = \left(\mathbf{I} + \mathbf{T}^{\prime} \mathbf{T} \right)^{-1} \mathbf{G}^{1/2} \] 这里 \(\mathbf{T} = \mathbf{R}^{-1/2} \mathbf{Z} \mathbf{G}^{1/2}\) 。

   其中 \(\mathbf{I} + \mathbf{T}^{\prime} \mathbf{T}\) 是一个正定矩阵，因为其二次型 \(\mathbf{x}^{\prime} \left( \mathbf{I} + \mathbf{T}^{\prime} \mathbf{T} \right) \mathbf{x} = \mathbf{x}^{\prime}\mathbf{x} + (\mathbf{T} \mathbf{x})^{\prime} (\mathbf{T} \mathbf{x}) > 0\) ，因此其逆矩阵 \((\mathbf{I} + \mathbf{T}^{\prime} \mathbf{T})^{\prime}\) 也是正定矩阵，其对角线元素均大于 0。因为 \(\mathbf{G}^{1/2}\) 是一个对角矩阵，根据上面的公式 \(\mathbf{G}^{1/2}\mathbf{W} = \left(\mathbf{I} + \mathbf{T}^{\prime} \mathbf{T} \right)^{-1} \mathbf{G}^{1/2}\) ，因此 \(\mathbf{W}\) 矩阵的对角线元素也均大于 0， 因此 \(\mathbf{W}\) 矩阵正定，因此 \(\operatorname{tr}\left(\mathbf{W}_{ii} \right) > 0\) 。

3. \(q_{i}- \operatorname{tr} \left(\mathbf{W}_{ii}\right) \geq 0\)

   证明：根据上面的公式，我们有 \[ \frac{ \mathbf{q}_{i} - \operatorname{tr} \left(\mathbf{W}_{ii}\right)} {\sigma_{i}^{2}} = \operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime} \right) = \operatorname{tr}\left( \mathbf{Z}_{i}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{Z}_{i} \right) = \operatorname{tr}\left( \mathbf{Z}_{i}^{\prime} \hat{\mathbf{V}}^{-1} \hat{\mathbf{V}} \hat{\mathbf{V}}^{-1} \mathbf{Z}_{i} \right) = \operatorname{var}( \mathbf{Z}_{i}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{y}) \geq 0 \] 因此我们可以看出每一轮迭代产生的 \(\hat{\sigma}_{i}^{2}\) 均为正数，即均在参数空间中。（不知道如何证明 \(\hat{\sigma}_{0}^{2}=\mathbf{y}^{\prime}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}}-\mathbf{Z} \hat{\mathbf{u}}) /N\) 大于 0 ）

REML

限制性最大似然 (restricted (residual) maximum likelihood, REML) ，由 Patterson 和 Thompson 在 1971 年提出。我们强调 REML 方法的原因是在标准的线性模型中，样本方差 \(s^{2}\) 也是 REML 的估计值。同时，REML 是一种一般的方法，比如 REML 在不平衡数据中也可以使用。在某些平衡数据中，REML 可能有解析解 (闭式解)。REML 同时是最佳二次无偏估计值。

REML 的思想是对数据 \(\mathbf{K'y}\) 进行最大似然估计，而不是 \(\mathbf{y}\) ，这里 \(\mathbf{K}^{\prime}\) 是我们人为挑选的，以使得 \(\mathbf{K'y}\) 的分布只包含方差组分，而不包含 \(\boldsymbol{\beta}\) 。为了实现这一点，矩阵 \(\mathbf{K}^{\prime}\) 需要满足 \(\mathbf{K' X}=\mathbf{O}\) ，因此 \(E(\mathbf{K' y})=\mathbf{K' X}\boldsymbol{\beta}=\mathbf{0}\) 。为了方便，我们同时要求 \(\mathbf{K}^{\prime}\) 是一个满秩矩阵。我们同时要求 \(\mathbf{K'y}\) 尽可能包含最多的关于方程组分的信息，因此 \(\mathbf{K}^{\prime}\) 必须有最大的行数。

定理： \(\mathbf{K}^{\prime}\) 是一个满秩矩阵，满足 \(\mathbf{K' X}=\mathbf{O}\) ，但是同时拥有最大的行数，因此它是一个 \((n-r) \times n\) 的矩阵 (\(r\) 为 \(\mathbf{X}\) 的秩)。进一步地说，\(\mathbf{K}^{\prime}\) 矩阵的形式必须为 \(\mathbf{K}^{\prime}=\mathbf{C}(\mathbf{I}-\mathbf{H})=\mathbf{C}\left[\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{X}\right)^{-} \mathbf{X}^{\prime}\right]\) ，这里 \(\mathbf{C}\) 是一个 \((n-r) \times n\) 的满秩矩阵。

证明：\(\mathbf{K}^{\prime}\) 的行 \(\mathbf{k}_{i}^{\prime}\) 满足等式 \(\mathbf{k}_{i}^{\prime} \mathbf{X}=\mathbf{0}^{\prime}\) ，转置得到 \(\mathbf{X}^{\prime} \mathbf{k}_{i}=\mathbf{0}\) （注意这里出现了符号重叠问题，这里的 \(\mathbf{k}_{i}\) 是第 \(i\) 行元素组成的列向量）。根据定理， 如果线性方程组 \(Ax = c\) 是相容的，那么所有可能的解为 \(\mathbf{x}=\mathbf{A}^{-} \mathbf{c}+\left(\mathbf{I}-\mathbf{A}^{-} \mathbf{A}\right) \mathbf{h}\) ，这里 \(\mathbf{A}^{-}\) 是某个特定的广义逆， \(h\) 为所有可能的值组成的向量。因此，这里的方程组的解为 \(\mathbf{k}_{i}=\left(\mathbf{I}-\mathbf{X}^{-'} \mathbf{X}^{\prime}\right) \mathbf{c}\) （书中有笔误，书里是 \(\mathbf{k}_{i}=\left(\mathbf{I}-\mathbf{X}^{-} \mathbf{X}\right) \mathbf{c}\) ，这里又有符号重叠问题，这里的 \(\mathbf{c}\) 就是上面定理中的 \(\mathbf{h}\) ，为所有可能的 \(n \times 1\) 的向量） ，因此 \(\mathbf{k}_{i}^{\prime}=\mathbf{c}^{\prime} \left(\mathbf{I}-\mathbf{X}\mathbf{X}^{-}\right)\) ，也就是说 \(\mathbf{K}^{\prime}\) 的行 \(\mathbf{k}_{i}^{\prime}\) 是 \(\mathbf{I}-\mathbf{X}\mathbf{X}^{-}\) 的行的所有可能的线性组合。

我们知道 \(\operatorname{rank}\left(\mathbf{X}\mathbf{X}^{-}\right)=\operatorname{rank}(\mathbf{X})=r\) ，并且 \(\left(\mathbf{I}-\mathbf{X}\mathbf{X}^{-}\right)^{2} = \mathbf{I}-2\mathbf{X} \mathbf{X}^{-}+\mathbf{X} \mathbf{X}^{-}= \mathbf{I}-\mathbf{X} \mathbf{X}^{-}\) ，因此 \(\mathbf{I}-\mathbf{X}\mathbf{X}^{-}\) 是一个幂等矩阵，因此 \(\operatorname{rank}\left(\mathbf{I}-\mathbf{X}\mathbf{X}^{-}\right)=\operatorname{tr}\left(\mathbf{I}-\mathbf{X}\mathbf{X}^{-}\right)=\operatorname{tr}(\mathbf{I})-\operatorname{tr}\left(\mathbf{X}\mathbf{X}^{-}\right)=n-r\) ，因此 \(\mathbf{I}-\mathbf{X}\mathbf{X}^{-}\) 的行的张成空间的维度为 \(n-r\) ，因此最多有 \(n-r\) 个线性无关的向量 \(\mathbf{k}_{i}^{\prime}\) ，也就是说 \(\mathbf{K}^{\prime}\) 矩阵的行数最大为 \(n-r\) 。

因为 \(\mathbf{k}_{i}^{\prime}=\mathbf{c}^{\prime} \left(\mathbf{I}-\mathbf{X}\mathbf{X}^{-}\right)\) ，那么一定存在一个 \((n-r) \times n\) 的满秩矩阵 \(\mathbf{C}\) ，使得 \(\mathbf{K}^{\prime}\) 矩阵可以表示为 \(\mathbf{K}^{\prime}=\mathbf{C}\left(\mathbf{I}-\mathbf{X}\mathbf{X}^{-}\right)\) ，其中 \(\mathbf{K}^{\prime}\) 矩阵的每一行均是 \(\mathbf{I}-\mathbf{X}\mathbf{X}^{-}\) 的行的线性组合，并且由于 \(\mathbf{K}^{\prime}\) 矩阵满秩， 因此选择 \(\mathbf{C}\) 矩阵的标准是必须使得 \(\mathbf{K}^{\prime}\) 矩阵的行之间线性无关。根据广义逆的性质，我们知道 \(\left(\mathbf{X}^{\prime} \mathbf{X}\right)^{-} \mathbf{X}^{\prime}\) 是 \(\mathbf{X}\) 的一个广义逆，因此我们可以这里的 \(\mathbf{X}^{-}\) 可以选择为 \(\left(\mathbf{X}^{\prime} \mathbf{X}\right)^{-} \mathbf{X}^{\prime}\) ，因此此时 \(\mathbf{K}^{\prime}\) 矩阵可以表示为 \(\mathbf{C}(\mathbf{I}-\mathbf{H})=\mathbf{C}\left[\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{X}\right)^{-} \mathbf{X}^{\prime}\right]\) 。

满足要求的 \(\mathbf{K}^{\prime}\) 矩阵有无数个，但是这不影响我们使用。同时我们注意到校正固定效应的残差为 \(\hat{\boldsymbol{\varepsilon}} = (\mathbf{I}-\mathbf{H}) \mathbf{y}\) ，因此

\(\mathbf{K' y}=\mathbf{C}(\mathbf{I}-\mathbf{H}) \mathbf{y} = \mathbf{C} \hat{\boldsymbol{\varepsilon}}\) 的每个元素都是所有残差的线性组合，这也是 residual maximum likelihood 这个名称的由来。

而 \(\mathbf{K' y}\) 的分布见下面的定理。

定理：在混合线性模型中，假设 \(\mathbf{y} \sim N_{n}(\mathbf{X} \boldsymbol{\beta}, \mathbf{V})\) ，其中 \(\mathbf{V}=\sum_{i=0}^{m} \sigma_{i}^{2} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime}\) ，那么 \[ \mathbf{K' y} \sim N_{n-r}\left(\mathbf{0}, \mathbf{K' \mathbf{V}} \mathbf{K}\right) \] 证明很简单，因为 \(\mathbf{y}\) 满足正态，所以 \(\mathbf{K' y}\) 也满足正态，同时根据 \(\mathbf{K' X}=\mathbf{0}\) 易得 \(\mathbf{K' y}\) 的均值为 \(\mathbf{0}\) , 其方差为 \(\mathbf{K^{\prime} \mathbf{V}} \mathbf{K}\) ，得证。

因此 \(\mathbf{K^{\prime} y}\) 的分布只与 \(m+1\) 个方差组分有关。为了估计方差组分，REML 的下一步是针对这些方差组分最大化 \(\mathbf{K^{\prime} y}\) 的似然值。

对似然函数求偏导，使之为零，我们得到下面的定理。

定理：在上面的模型中， \(m+1\) 个方差组分 \(\sigma_{0}^{2}, \ldots, \sigma_{m}^{2}\) 的估计值满足下面的方程组，其中 \(i=0, \ldots, m\) \[ \operatorname{tr}\left[\mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K^{\prime} Z}_{i} \mathbf{Z}_{i}^{\prime}\right]=\mathbf{y}^{\prime} \mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K^{\prime} Z}_{i} \mathbf{Z}_{i}^{\prime} \mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K}^{\prime} \mathbf{y} \] 证明：因为 \(E(\mathbf{K^{\prime} y})=\mathbf{0}\) ，\(\mathbf{K^{\prime} y}\) 的对数似然值为 \[ \begin{aligned} \ln L\left(\sigma_{0}^{2}, \ldots, \sigma_{m}^{2}\right)=& -\frac{n-r}{2} \ln (2 \pi)-\frac{1}{2} \ln \left|\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right|-\frac{1}{2} \mathbf{y}^{\prime} \mathbf{K} \left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K} \right)^{-1} \mathbf{K^{\prime} y} \\ =& -\frac{n-r}{2} \ln (2 \pi)-\frac{1}{2} \ln \left|\mathbf{K}^{\prime}\left(\sum_{i=0}^{m} \sigma_{i}^{2} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime}\right) \mathbf{K} \right| \\ &-\frac{1}{2} \mathbf{y}^{\prime} \mathbf{K} \left[\mathbf{K}^{\prime}\left(\sum_{i=0}^{m} \sigma_{i}^{2} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime}\right) \mathbf{K}\right]^{-1} \mathbf{K^{\prime} y} \end{aligned} \] 根据矩阵求导的性质，假设 \(\mathbf{A}\) 是一个 \(n \times n\) 的非奇异矩阵，其元素 \(a_{ij}\) 是关于标量 \(x\) 的函数。我们定义 \(\partial \mathbf{A} / \partial x\) 是一个 \(n \times n\) 的矩阵，其元素为 \(\partial a_{i j} / \partial x\) 。那么存在 \[ \frac{\partial \mathbf{A}^{-1}}{\partial x}=-\mathbf{A}^{-1} \frac{\partial \mathbf{A}}{\partial x} \mathbf{A}^{-1} \] 和 \[ \frac{\partial \log |\mathbf{A}|}{\partial x}=\operatorname{tr}\left(\mathbf{A}^{-1} \frac{\partial \mathbf{A}}{\partial x}\right) \] 利用这两条性质，我们求 \(\ln L\left(\sigma_{0}^{2}, \ldots, \sigma_{m}^{2}\right)\) 对每一个 \(\sigma_{i}^{2}\) 的偏导数，得到 \[ \begin{aligned} \frac{\partial}{\partial \sigma_{i}^{2}} \ln L\left(\sigma_{0}^{2}, \ldots, \sigma_{m}^{2}\right)=&-\frac{1}{2} \operatorname{tr}\left(\left(\mathbf{K^{\prime} \mathbf{V}} \mathbf{K}\right)^{-1}\left[\frac{\partial}{\partial \sigma_{i}^{2}}\left(\mathbf{K^{\prime} \mathbf{V}} \mathbf{K}\right)\right]\right) \\ &+\frac{1}{2} \mathbf{y}^{\prime} \mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1}\left[\frac{\partial}{\partial \sigma_{i}^{2}}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)\right]\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K^{\prime} y} \\ =&-\frac{1}{2} \operatorname{tr}\left[\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K^{\prime} Z}_{i} \mathbf{Z}_{i}^{\prime} \mathbf{K}\right] \\ &+\frac{1}{2} \mathbf{y}^{\prime} \mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K^{\prime} Z}_{i} \mathbf{Z}_{i}^{\prime} \mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K}^{\prime} \mathbf{y} \\ =&-\frac{1}{2} \operatorname{tr}\left[\mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K^{\prime} Z}_{i} \mathbf{Z}_{i}^{\prime}\right] \\ &+\frac{1}{2} \mathbf{y}^{\prime} \mathbf{K} \left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K^{\prime} Z}_{i} \mathbf{Z}_{i}^{\prime} \mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K}^{\prime} \mathbf{y} \end{aligned} \] 将这个式子设为0，得到上面的结果。

证毕。

有一点很有意思，根据二次型的期望公式 \(E\left(\mathbf{y}^{\prime} \mathbf{A y}\right)=\operatorname{tr}(\mathbf{A} \mathbf{V})+\boldsymbol{\mu}^{\prime} \mathbf{A} \boldsymbol{\mu}\) ，上面式子的右手项的期望 (这里设期望公式里的 \(\mathbf{y}\) 是 \(\mathbf{K^{\prime}y}\) ，\(\mathbf{A} = \left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K} \right)^{-1} \mathbf{K^{\prime} Z}_{i} \mathbf{Z}_{i}^{\prime} \mathbf{K} \left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1}\) ) 正好是左手项。

SEARLE (1979) 证明 \[ \mathbf{K} \left(\mathbf{K}^{\prime} \mathbf{V K}\right)^{-1} \mathbf{K}^{\prime}=\mathbf{P} \] 证明如下：

\(\mathbf{K K}^{+}=\mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{K}\right)^{-1} \mathbf{K}^{\prime}\) 和 \(\mathbf{X X}^{+}=\mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{X}\right)^{+} \mathbf{X}^{\prime}\) 是对称幂等矩阵，对称证明略，幂等证明如下，\((\mathbf{K K}^{+})^{2} = \mathbf{K K}^{+} \mathbf{K K}^{+} = \mathbf{K K}^{+}\) 。我们已知 \(\mathbf{K}^{\prime} \mathbf{X}=\mathbf{0}\) ，因此 \(\mathbf{K K}^{+} \mathbf{X}=\mathbf{0}\) 并且 \(\mathbf{X X}^{+} \mathbf{K}=\mathbf{0}\) 。因此

\(\mathbf{T}=\mathbf{I}-\mathbf{X} \mathbf{X}^{+}-\mathbf{K} \mathbf{K}^{+}\) 是一个对称幂等矩阵，幂等证明如下（其实直接根据对称幂等矩阵的性质即可证明，如果 \(\mathbf{A}\) 对称幂等，那么 \(\mathbf{I-A}\) 对称幂等 ）： \[ \begin{aligned} &\left( \mathbf{I}-\mathbf{X} \mathbf{X}^{+}-\mathbf{K} \mathbf{K}^{+} \right)\left( \mathbf{I}-\mathbf{X} \mathbf{X}^{+}-\mathbf{K} \mathbf{K}^{+} \right) \\ &= \mathbf{I}-\mathbf{X} \mathbf{X}^{+}-\mathbf{K} \mathbf{K}^{+} -\mathbf{X} \mathbf{X}^{+} + \mathbf{X} \mathbf{X}^{+} \mathbf{X} \mathbf{X}^{+} +\mathbf{K} \mathbf{K}^{+} \mathbf{X} \mathbf{X}^{+} - \mathbf{K} \mathbf{K}^{+} + \mathbf{X} \mathbf{X}^{+} \mathbf{K} \mathbf{K}^{+} + \mathbf{K} \mathbf{K}^{+} \mathbf{K} \mathbf{K}^{+} \\ &= \mathbf{I}-\mathbf{X} \mathbf{X}^{+}-\mathbf{K} \mathbf{K}^{+} -\mathbf{X} \mathbf{X}^{+} + \mathbf{X} \mathbf{X}^{+} - \mathbf{K} \mathbf{K}^{+} + \mathbf{K} \mathbf{K}^{+} \\ &= \mathbf{I}-\mathbf{X} \mathbf{X}^{+}-\mathbf{K} \mathbf{K}^{+} \\ \end{aligned} \] 因此 \[ \begin{aligned} \operatorname{tr}\left(\mathbf{T T}^{\prime}\right)=\operatorname{tr}\left(\mathbf{T}^{2}\right)=\operatorname{tr}(\mathbf{T}) &=\operatorname{tr}(\mathbf{I})-\operatorname{tr}\left(\mathbf{X X}^{+}\right)-\operatorname{tr}\left(\mathbf{K K}^{+}\right) \\ &=N-r_{\mathbf{X}}-r_{\mathbf{K}} \\ &=N-r_{\mathbf{X}}-\left(N-r_{\mathbf{X}}\right) \\ &=0 \end{aligned} \] 根据幂等矩阵的性质，我们知道幂等矩阵的秩等于迹，因此 \(\operatorname{rank}(\mathbf{T}) = 0\) ，说明 \(\mathbf{T}=\mathbf{0}\) （只有零矩阵的秩为 0）。因此，我们有 \[ \mathbf{I}-\mathbf{X X}^{+}=\mathbf{K K}^{+} \] 因为 \(\mathbf{V}\) 是一个正定矩阵，因此总是存在 \(\mathbf{V}^{1/2}\) 矩阵使得 \(\mathbf{V} = (\mathbf{V}^{1/2})^{2}\) 。那么由于 \(\mathbf{K}^{\prime} \mathbf{X}=\mathbf{0}\) ，那么我们知道 \(\left(\mathbf{V}^{1/2} \mathbf{K}\right)^{\prime} \mathbf{V}^{-1/2} \mathbf{X}=\mathbf{0}\) ，之前适应 \(\mathbf{K}\) 和 \(\mathbf{X}\) 的结果均适用于 \(\mathbf{V}^{1/2} \mathbf{K}\) 和 \(\mathbf{V}^{-1/2} \mathbf{X}\) 。我们替换 \(\mathbf{I}-\mathbf{X X}^{+}=\mathbf{K K}^{+}\) 中的MP逆矩阵，我们得到 \(\mathbf{I}-\mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{X}\right)^{+} \mathbf{X}^{\prime} =\mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{K}\right)^{-1} \mathbf{K}^{\prime}\) 。将 \(\mathbf{V}^{1/2} \mathbf{K}\) 和 \(\mathbf{V}^{-1/2} \mathbf{X}\) 带入得到 \[ \mathbf{I}-\mathbf{V}^{-1/2} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{+} \mathbf{X}^{\prime} \mathbf{V}^{-1/2}=\mathbf{V}^{1/2} \mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K}^{\prime} \mathbf{V}^{1/2} \] 因此 （注意下面的 \(\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-}\) 是 MP 逆 ） \[ \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} \] 得证。

其中 \(\hat{\mathbf{P}}=\hat{\mathbf{V}}^{-1}-\hat{\mathbf{V}}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1}\) ，将该式带入，得到 \[ \operatorname{tr}\left(\hat{\mathbf{P}} \mathbf{V}_{i}\right)=\mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{V}_{i} \hat{\mathbf{P}} {\mathbf{y}} \] 注意最大似然法的估计公式为 \(\operatorname{tr}\left(\hat{\mathbf{V}}^{-1} \mathbf{V}_{i}\right) =\mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{V}_{i} \hat{\mathbf{P}} {\mathbf{y}}\) ，因此 REML 方法相比与 ML 方法就是将左手项的 \(\mathbf{\hat{V}}^{-1}\) 替换为 \(\mathbf{\hat{P}}\) 。

我们有 $$ \[\begin{aligned} \mathbf{PVP} &= \left(\mathbf{V}^{-1}-\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}-\mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \right) \\ &= \left(\mathbf{I}-\mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \right) \left(\mathbf{V}^{-1}-\mathbf{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \right) \\ &= \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{V}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \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{X}\left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{V}^{-1} \\ &= \mathbf{P} \quad \because \text{设} \left(\mathbf{X}^{\prime} \mathbf{V}^{-1} \mathbf{X}\right)^{-} \text{是自反广义逆矩阵} \\ \end{aligned}\] \[ 因此，我们有 \] \[\begin{aligned} \operatorname{tr}\left(\hat{\mathbf{P}} \mathbf{V}_{i}\right) &=\operatorname{tr}\left(\hat{\mathbf{P}} \hat{\mathbf{V}} \hat{\mathbf{P}} \mathbf{V}_{i}\right) \\ &=\operatorname{tr}\left(\hat{\mathbf{P}} \mathbf{V}_{i} \mathbf{P} \hat{\mathbf{V}}\right) \\ &=\operatorname{tr}\left(\hat{\mathbf{P}} \mathbf{V}_{i} \hat{\mathbf{P}} \sum_{j=0}^{m} \mathbf{V}_{j}\right) \hat{\sigma}_{j}^{2} \\ &=\sum_{j=0}^{m} \operatorname{tr}\left(\hat{\mathbf{P}} \mathbf{V}_{i} \hat{\mathbf{P}} \mathbf{V}_{j}\right) \hat{\sigma}_{j}^{2} \end{aligned}\]

\[ 因而，对于 $i,j = 0,1,\cdots, m$ ，我们有 \] {j=0}^{m}{( {i} )}{{j}^{2}}=\left{{} {i} } $$ 这就是 REML 的估计方程，它也必须通过迭代求解。

根据这个定理，我们可以得到 \(m+1\) 个方程，而我们总共有 \(m+1\) 个未知数。在某些情况下，这些方程组可以进一步简化从而得到闭式解。但是在大多数情况下，这些方程组是无法直接求解的。

我们一般会用一些迭代的方法来进行估计 (Rao 1997 pp. 104 – 105, McCulloch and Searle 2001, pp. 265 – 269)。我们注意到上面定理中的 \(m +1\) 个方程组成的方程组可以写为 \[ \mathbf{M} \boldsymbol{\sigma}=\mathbf{q} \] 这里 \(\boldsymbol{\sigma}=\left(\sigma_{0}^{2}, \sigma_{1}^{2}, \cdots, \sigma_{m}^{2}\right)^{\prime}\) ，这里 \(\mathbf{M}\) 是一个非奇异的 \((m+1) \times(m+1)\) 的矩阵，其中的元素 \(M_{ij}=\operatorname{tr}\left(\hat{\mathbf{P}} \mathbf{V}_{i} \hat{\mathbf{P}}{\mathbf{V}_{j}}\right)\) ，然后 \(\mathbf{q}\) 是一个 \((m+1) \times 1\) 的向量，其中的元素 \(q_{i} = \mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{V}_{i} \hat{\mathbf{P}} \mathbf{y}\) 。

我们注意到 \(\mathbf{M} \boldsymbol{\sigma}=\mathbf{q}\) 这个式子貌似更加复杂，因为 \(\mathbf{M}\) 和 \(\mathbf{q}\) 中都含有未知数 \(\boldsymbol{\sigma}\) 。然而，这个方程组对于我们从一个初始值 \(\boldsymbol{\sigma}_{(1)}\) 逐步迭代来估计的方式而言很有用。我们在第 \(t\) 步可以用 \(\boldsymbol{\sigma}_{(t)}\) 来计算 \(\mathbf{M}_{(t)}\) 和 \(\mathbf{q}_{(t)}\) ，那么 \(\boldsymbol{\sigma}_{(t+1)}=\mathbf{M}_{(t)}^{-1} \mathbf{q}_{(t)}\) ，然后一直迭代直到收敛。

估计值的方差

首先我们证明 \[ \begin{aligned} \frac{\partial \mathbf{P}}{\partial \sigma_i^2} & =\frac{\partial}{\partial \sigma_i^2} \mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K}^{\prime} \\ & =\mathbf{K} \left( \frac{\partial}{\partial \sigma_i^2} \left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \right) \mathbf{K}^{\prime} \quad \because \frac{\partial \mathbf{P} \mathbf{A} \mathbf{Q} }{\partial x} = \mathbf{P} \frac{\partial \mathbf{A}}{\partial x} \mathbf{Q} \\ & =-\mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K} \right)^{-1}\frac{\partial \mathbf{K}^{\prime} \mathbf{V} \mathbf{K}}{\partial \sigma_i^2}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K}^{\prime} \quad \because \frac{\partial \mathbf{A}^{-1}}{\partial x}=-\mathbf{A}^{-1} \frac{\partial \mathbf{A}}{\partial x} \mathbf{A}^{-1}\\ & =-\mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K}^{\prime} \frac{\partial \mathbf{V}}{\partial \sigma_i^2} \mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V K}\right)^{-1} \mathbf{K}^{\prime} \\ & =-\mathbf{P} \frac{\partial \mathbf{V}}{\partial \sigma_i^2} \mathbf{P}\\ &=-\mathbf{P} \mathbf{Z}_i \mathbf{Z}_i^{\prime} \mathbf{P}\\ \end{aligned} \] 我们已知 $$ \[\begin{aligned} \frac{\partial l_{\mathbf{R}}}{\partial \sigma_i^2} &=-\frac{1}{2} \operatorname{tr}\left[\mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K^{\prime} Z}_{i} \mathbf{Z}_{i}^{\prime}\right] \\ &+\frac{1}{2} \mathbf{y}^{\prime} \mathbf{K} \left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K^{\prime} Z}_{i} \mathbf{Z}_{i}^{\prime} \mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{V} \mathbf{K}\right)^{-1} \mathbf{K}^{\prime} \mathbf{y}\\ & =-\frac{1}{2} \operatorname{tr}\left(\mathbf{P}\mathbf{Z}_i \mathbf{Z}_i^{\prime}\right)+\frac{1}{2} \mathbf{y}^{\prime} \mathbf{P} \mathbf{Z}_i \mathbf{Z}_i^{\prime} \mathbf{P y} \end{aligned}\] \[ 求二阶导得到l \] \[\begin{aligned} \frac{\partial^2 l_{R}}{\partial \sigma_i^2 \partial \sigma_j^2} & =\frac{1}{2} \operatorname{tr}\left(\mathbf{P Z}_j \mathbf{Z}_j^{\prime} \mathbf{P}\mathbf{Z}_i \mathbf{Z}_i^{\prime}\right)-\frac{1}{2} \mathbf{y}^{\prime} \mathbf{P} \mathbf{Z}_j\mathbf{Z}_j^{\prime} \mathbf{P} \mathbf{Z}_i \mathbf{Z}_i^{\prime} \mathbf{P y}-\frac{1}{2} \mathbf{y}^{\prime} \mathbf{P} \mathbf{Z}_i \mathbf{Z}_i^{\prime} \mathbf{P }\mathbf{Z}_j \mathbf{Z}_j^{\prime} \mathbf{P y} \\ & =\frac{1}{2} \operatorname{tr}\left(\mathbf{P Z}_j \mathbf{Z}_j^{\prime} \mathbf{P } \mathbf{Z}_i\mathbf{Z}_i^{\prime}\right)-\mathbf{y}^{\prime} \mathbf{P} \mathbf{Z}_j \mathbf{Z}_j^{\prime} \mathbf{P} \mathbf{Z}_i \mathbf{Z}_i^{\prime} \mathbf{P y} \end{aligned}\]

$$ 这里利用了两个性质：\(\frac{\partial \operatorname{tr} (\mathbf{A} \mathbf{B}) }{\partial x} = \operatorname{tr} (\mathbf{A}\frac{\partial \mathbf{B} }{\partial x}) + \operatorname{tr} (\frac{\partial \mathbf{A} }{\partial x}\mathbf{B})\) 和 \(\frac{\partial \mathbf{A}_{1} \mathbf{A}_{2} \cdots \mathbf{A}_{n} }{\partial x} = (\frac{\partial \mathbf{A}_{1}}{\partial x} \mathbf{A}_{2} \cdots \mathbf{A}_{n})+\cdots+( \mathbf{A}_{1} \mathbf{A}_{2} \cdots \frac{\partial \mathbf{A}_{n}}{\partial x})\) ，证明过程参考 Linear models in statistics 。

再求期望得到信息矩阵， 这里利用二次型的期望公式 \(E\left(\mathbf{y}^{\prime} \mathbf{A y}\right)=\operatorname{tr}(\mathbf{A} \mathbf{V})+\boldsymbol{\mu}^{\prime} \mathbf{A} \boldsymbol{\mu} .\) \[ \begin{aligned} & -E\left(\frac{\partial^2 l_{\mathbf{R}}}{\partial \sigma_i^2 \sigma_j^2}\right)\\ &=-\frac{1}{2} \operatorname{tr}\left(\mathbf{P Z}_j \mathbf{Z}_j^{\prime} \mathbf{P Z_{ i }} \mathbf{Z}_i^{\prime}\right)+\operatorname{tr}\left(\mathbf{P Z}_j \mathbf{Z}_j^{\prime} \mathbf{P }\mathbf{Z}_i \mathbf{Z}_i^{\prime} \mathbf{P V}\right)-\boldsymbol{\beta}^{\prime} \mathbf{X}^{\prime} \mathbf{P } \mathbf{Z}_j\mathbf{Z}_j^{\prime} \mathbf{P Z _ { i }} \mathbf{Z}_i^{\prime} \mathbf{P X} \boldsymbol{\beta} \\ & =-\frac{1}{2} \operatorname{tr}\left(\mathbf{P Z}_j \mathbf{Z}_j^{\prime} \mathbf{P Z _ { i }} \mathbf{Z}_i^{\prime}\right)+\operatorname{tr}\left(\mathbf{Z}_j \mathbf{Z}_j^{\prime} \mathbf{P Z}_i \mathbf{Z}_i^{\prime} \mathbf{P V P}\right)+\mathbf{0} \quad \text {, because } \mathbf{P X}= \mathbf{K} \left(\mathbf{K}^{\prime} \mathbf{V K}\right)^{-1} (\mathbf{K}^{\prime} \mathbf{X})=\mathbf{0} \text {, } \\ & =\frac{1}{2} \operatorname{tr}\left(\mathbf{P Z}_j \mathbf{Z}_j^{\prime} \mathbf{P Z}_i \mathbf{Z}_i^{\prime}\right) \quad \text {, because } \mathbf{P V P}=\mathbf{P} . \\ & \end{aligned} \] 因此估计值的协方差矩阵为 \[ \begin{aligned} \operatorname{var}\left(\tilde{\sigma}_{\mathrm{REML}}^2\right) & \simeq 2\left[\left\{_{\mathrm{m}} \operatorname{tr}\left(\mathbf{P} \mathbf{Z}_i \mathbf{Z}_i^{\prime} \mathbf{P} \mathbf{Z}_j \mathbf{Z}_j^{\prime}\right)\right\}_{i, j=0}^r\right]^{-1} \\ & \simeq 2\left[\left\{_{\mathrm{m}} \operatorname{sesq}\left(\mathbf{Z}_i^{\prime} \mathbf{P} \mathbf{Z}_j\right)\right\}_{i, j=0}^r\right]^{-1}, \end{aligned} \] 这和 ML 的式子格式相同，只是将 \(\mathbf{V}^{-1}\) 替换成了 \(\mathbf{P}\) 。

EM 算法

混合详细模型如下： \[ \mathbf{y} = \mathbf{X} \boldsymbol{\beta} + \mathbf{Z}\mathbf{u} + \mathbf{e} \] 对于公式 \(\operatorname{tr}\left(\hat{\mathbf{P}} \mathbf{V}_{i}\right)=\mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{V}_{i} \hat{\mathbf{P}} {\mathbf{y}}\) ，我们在等式两边同乘以 \(\hat{\sigma}^{2}_{i}\) ，并求总和，得到： \[ \begin{aligned} \operatorname{tr}\left(\hat{\mathbf{P}} \sum \mathbf{\hat{V}}_{i} \hat{\sigma}_{i}^{2}\right) &=\mathbf{y}^{\prime} \hat{\mathbf{P}} \sum \mathbf{V}_{i} \hat{\sigma}_{i}^{2} \hat{\mathbf{P}}\mathbf{y} \\ \operatorname{tr}\left(\hat{\mathbf{P}} \mathbf{\hat{V}}\right) &=\mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{\hat{V}} \hat{\mathbf{P}}\mathbf{y} \\ \end{aligned} \] 我们先看左手项 (下面的 \(N\) 是表型数目，\(r\) 为 \(\mathbf{X}\) 的秩) $$ \[\begin{aligned} \operatorname{tr}(\hat{\mathbf{P}} \hat{\mathbf{V}}) &=\operatorname{tr}\left[\left(\hat{\mathbf{V}}^{-1}-\hat{\mathbf{V}}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X}\right) \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1}\right) \hat{\mathbf{V}}\right] \\ &=\operatorname{tr}(\mathbf{I})-\operatorname{tr}\left[ \hat{\mathbf{V}}^{-1} \mathbf{X} \left(\mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \right] \\ &=\operatorname{tr}(\mathbf{I})-\operatorname{tr}\left[\left(\mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X}\right] \quad \because \left(\mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X} \text{ 是一个对称幂等矩阵} \\ &=\operatorname{tr}(\mathbf{I})-\operatorname{Rank}\left[\left(\mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X}\right] \\ &=N-r \end{aligned}\] \[ 我们再看右手项，在之前 ML 的推导中，我们知道 $\hat{\mathbf{P}} \mathbf{y} =\hat{\mathbf{V}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})$ ，因此我们有 \] \[\begin{aligned} \mathbf{y}^{\prime} \hat{\mathbf{P}} \mathbf{\hat{V}} \hat{\mathbf{P}}\mathbf{y} &=(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})^{\prime} \hat{\mathbf{V}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})\\ &=\mathbf{y}^{\prime} \hat{\mathbf{V}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}}) - \hat{\boldsymbol{\beta}}^{\prime} \left( \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{y} - \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X} \hat{\boldsymbol{\beta}} \right)\\ &=\mathbf{y}^{\prime} \hat{\mathbf{V}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}}) \quad\left(\because \mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{X} \hat{\boldsymbol{\beta}}=\mathbf{X}^{\prime} \hat{\mathbf{V}}^{-1} \mathbf{y}\right)\\ &=\mathbf{y}^{\prime}\left[\hat{\mathbf{R}}^{-1}-\hat{\mathbf{R}}^{-1} \mathbf{Z}\left(\mathbf{Z}^{\prime} \hat{\mathbf{R}}^{-1} \mathbf{Z}+\hat{\mathbf{G}}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \hat{\mathbf{R}}^{-1}\right](\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})\\ &=\mathbf{y}^{\prime} \hat{\mathbf{R}}^{-1}\left[(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})-\mathbf{Z}\left(\mathbf{Z}^{\prime} \hat{\mathbf{R}}^{-1} \mathbf{Z}+\hat{\mathbf{G}}^{-1}\right)^{-1} \mathbf{Z}^{\prime} \hat{\mathbf{R}}^{-1}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})\right]\\ &=\mathbf{y}^{\prime} \hat{\mathbf{R}}^{-1}[(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}})-\mathbf{Z} \hat{\mathbf{u}}]\\ &\left(\because \hat{\mathbf{u}}=\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1}\right)^{-1}\left(\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{y}-\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{X} \hat{\boldsymbol{\beta}}\right)\right)\\ &=\frac{1}{\hat{\sigma}_{0}^{2}} \mathbf{y}^{\prime}(\mathbf{y}-\mathbf{X} \hat{\boldsymbol{\beta}}-\mathbf{Z} \hat{\mathbf{u}}) \quad \because \text{假设为单性状模型} \end{aligned}\]

\[ 因此，我们有 \] _{0}^{2}=({} -^{} -^{} ) /(N-r) $$ 对于某个随机因子的方差组分 \(\hat{\sigma}_{i}^{2}\) ，其迭代公式推导如下。

首先我们需要证明一个引理，因为 \(\mathbf{V} = \mathbf{Z}\mathbf{G}\mathbf{Z}^{\prime}+\mathbf{R}\) ，如果我们将 \(\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}=\mathbf{0}\) ，将新的式子定义为 \(\mathbf{S}\) ，得到： \[ \mathbf{S} = \mathbf{R}^{-1}-\mathbf{R}^{-1} \mathbf{X}\left(\mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{X}\right)^{-} \mathbf{X}^{\prime} \mathbf{R}^{-1} \] 那么我们有下面的引理。

引理： \[ \mathbf{P}=\mathbf{S-S Z\left(D^{-1}+Z^{\prime} S Z\right)^{-1} Z^{\prime} S} \] 证明：因为我们有 \(\mathbf{P} = \mathbf{K} \left(\mathbf{K}^{\prime} \mathbf{V K}\right)^{-1} \mathbf{K}^{\prime}\) ，代入 \(\mathbf{Z}=\mathbf{0}\) 得到 \[ \mathbf{S} = \mathbf{K} \left(\mathbf{K}^{\prime} \mathbf{R K}\right)^{-1} \mathbf{K}^{\prime} \] 因此我们有 \[ \begin{aligned} \mathbf{P}= & \mathbf{K}\left[\mathbf{K}^{\prime}(\mathbf{Z G Z^{\prime}}+\mathbf{R}) \mathbf{K}\right]^{-1} \mathbf{K}^{\prime}=\mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{R K}+\mathbf{K}^{\prime} \mathbf{Z G Z^{\prime}} \mathbf{K}\right)^{-1} \mathbf{K}^{\prime} \\ = & \mathbf{K}\left\{\left(\mathbf{K}^{\prime} \mathbf{R K}\right)^{-1}-\left(\mathbf{K}^{\prime} \mathbf{R K}\right)^{-1} \mathbf{K}^{\prime} \mathbf{Z G}\left[\mathbf{I}+\mathbf{Z}^{\prime} \mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{R K}\right)^{-1} \mathbf{K}^{\prime} \mathbf{Z G}\right]^{-1}\right. \\ & \left.\times \mathbf{Z}^{\prime} \mathbf{K}\left(\mathbf{K}^{\prime} \mathbf{R K}\right)^{-1}\right\} \mathbf{K}^{\prime}, \text { 使用第二个舒尔补公式 }\\ = & \mathbf{S}-\mathbf{S Z G}\left(\mathbf{I}+\mathbf{Z}^{\prime} \mathbf{S Z G}\right)^{-1} \mathbf{Z}^{\prime} \mathbf{S} \\ = & \mathbf{S}-\mathbf{S Z}\left(\mathbf{G}^{-1}+\mathbf{Z}^{\prime} \mathbf{S Z}\right)^{-1} \mathbf{Z}^{\prime} \mathbf{S} . \end{aligned} \] 注意这个公式和 \(\mathbf{V}^{-1}\) 有着相似的格式，只是其中的 \(\mathbf{R}^{-1}\) 替换成了 \(\mathbf{S}\) \[ \mathbf{V}^{-1} = \mathbf{R}^{-1}-\mathbf{R^{-1} Z}\left(\mathbf{G}^{-1}+\mathbf{Z}^{\prime} \mathbf{R^{-1} Z}\right)^{-1} \mathbf{Z}^{\prime} \mathbf{R}^{-1} \] 舒尔补公式说明如下，如果一个非奇异分块矩阵中，如果 \(\mathbf{A}\) 可逆，则其逆矩阵可以表示为 \[ \left[\begin{array}{ll} \mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D} \end{array}\right]^{-1}=\left[\begin{array}{cc} \mathbf{A}^{-1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{array}\right]+\left[\begin{array}{c} -\mathbf{A}^{-1} \mathbf{B} \\ \mathbf{I} \end{array}\right]\left(\mathbf{D}-\mathbf{C A}^{-1} \mathbf{B}\right)^{-1}\left[-\mathbf{C A}^{-1} \quad \mathbf{I}\right] \] 其中矩阵 \(\mathbf{D}-\mathbf{C A}^{-1} \mathbf{B}\) 称为 \(\mathbf{A}\) 的舒尔补，Marsaglia 和 Styan (1974a,b) 给出了两个重要公式，分别为（第二个式子可以理解为将 \(\mathbf{A}\) 替换成了 \(\mathbf{-A}\) ）。 \[ \left(\mathbf{D}-\mathbf{C A}^{-1} \mathbf{B}\right)^{-1}=\mathbf{D}^{-1}+\mathbf{D}^{-1} \mathbf{C}\left(\mathbf{A}-\mathbf{B D}^{-1} \mathbf{C}\right)^{-1} \mathbf{B D}^{-1}, \]

\[ \left(\mathbf{D}+\mathbf{C A}^{-1} \mathbf{B}\right)^{-1}=\mathbf{D}^{-1}-\mathbf{D}^{-1} \mathbf{C}\left(\mathbf{A}+\mathbf{B D}^{-1} \mathbf{C}\right)^{-1} \mathbf{B D}^{-1} \text {, } \]

我们再看公式 \(\operatorname{tr}\left(\mathbf{P} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime} \right)=\mathbf{y}^{\prime} \mathbf{P} \mathbf{V}_{i} \mathbf{P} {\mathbf{y}}\) 的左手项 \(\operatorname{tr}\left(\mathbf{P} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime} \right)\) ，之前在 ML 的 EM 算法推导中，我们得到 \[ \operatorname{tr}\left(\mathbf{V}^{-1} \mathbf{Z}_i \mathbf{Z}_i^{\prime}\right)=\left[q_i-\operatorname{tr}\left(\mathbf{W}_{ii}\right)\right] / \sigma_i^2 \] 这里 \(\mathbf{W}_{ii}\) 是 \(\mathbf{W}=\left(\mathbf{I}+\mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z} \mathbf{G}\right)^{-1}\) 的第 \((i,i)\) 个子矩阵。

因为 \(\mathbf{P}\) 和 \(\mathbf{V}^{-1}\) 的差别就是将 \(\mathbf{R}^{-1}\) 替换成了 \(\mathbf{S}\)， 因此我们得到： \[ \operatorname{tr}\left(\mathbf{P Z}_i \mathbf{Z}_i^{\prime}\right)=\left[q_i-\operatorname{tr}\left(\mathbf{T}_{ii}\right)\right] / \sigma_i^2 \] 这里 \(\mathbf{T}_{ii}\) 是 \(\mathbf{W}\) 中将 \(\mathbf{R}^{-1}\) 替换成了 \(\mathbf{S}\) 的子矩阵，我们标记为 \[ \mathbf{T}=\left(\mathbf{I}+\mathbf{Z}^{\prime} \mathbf{S} \mathbf{Z} \mathbf{G}\right)^{-1} \] 因为 REML 公式和 ML 公式的右手项不变，因此我们得到： \[ \hat{\sigma}_{i}^{2}=\frac{\hat{\mathbf{u}}_{i}^{\prime} \hat{\mathbf{u}}_{i}+ \hat{\sigma}_{i}^{2} \operatorname{tr} \left(\mathbf{T}_{ii}\right) } {q_{i}} \] 或者 \[ \hat{\sigma}_{i}^{2}=\frac{\hat{\mathbf{u}}_{i}^{\prime} \hat{\mathbf{u}}_{i}}{q_{i}- \operatorname{tr} \left(\mathbf{T}_{ii}\right) } \] 通过上面这两个公式我们得到了 REML 的 EM 算法。在计算时，我们先给出一组初始的 \(\hat{\sigma}_{i}^{2}\) 值，由混合模型方程组求出相应的 \(\hat{\boldsymbol{\beta}}\) 和 \(\hat{\mathbf{u}}\) ，以及相应的系数矩阵逆矩阵元素；再由上面的两个公式求出一组新的 \(\hat{\sigma}_{i}^{2}\) 值，如此迭代下去，知道两次迭代得到的 \(\hat{\sigma}_{i}^{2}\) 值得差异小于给定阈值时，迭代收敛。

而且根据上面两个公式，可以得到 \(\hat{\sigma}_{0}^{2}\) 和 \(\hat{\sigma}_{i}^{2}\) 都比大于0，因此 EM 算法一定满足参数空间的要求。

EM 算法

EM 算法的全名为 expectation-maximization ，因为它在计算条件期望值和最大化简化的似然函数值间不断跳转。EM算法只会计算估计值，而不会计算估计值的方差。

EM 算法用于当数据增加时会简化问题难度的最大似然估计，应用 EM 算法的核心在于决定什么是完整数据（观测数据+未观测到的数据）。实际数据（观测数据）在 EM 算法中一般称为不完整的数据。因此在方差组分估计中，我们将表型 \(\mathbf{y}\) 视为不完整的数据，完整数据为表型 \(\mathbf{y}\) 加上未观测的随机效应 \(\mathbf{u}_{i}\) (\(i=1,2,\cdots,r\)) 。

如果我们知道随机效应 \(\mathbf{u}_{i}\) 的值，我们可以用其平方和均值估计随机效应的方差，该式为基于完整数据的最大似然估计值。 \[ \hat{\sigma}_{i}^{2} = \mathbf{u}_{i}^{\prime}\mathbf{u}_{i}/q_{i} \] 但是，实际上我们并不知道随机效应 \(\mathbf{u}_{i}\) 的值，但是 EM 算法给了我们一种基于 \(\mathbf{y}\) 估计 \(\mathbf{u}_{i}\) 的方法。首先我们采用一组参数的初始值，我们之后基于 \(\mathbf{y}\) 计算 \(\mathbf{u}_{i}^{\prime}\mathbf{u}_{i}\) 的条件期望（E 步, expectation step），然后利用其条件期望值使用最大似然法得到一组新的参数估计值（M 步， maximization step）。然后不断循环这两步直到收敛。

EM 算法的一个重要特征是，因为 EM 算法对完整数据执行最大似然估计，因此每一次迭代得到的估计值均在参数空间中。

联合分布

因为我们需要计算给定 \(\mathbf{y}\) 条件下 \(\mathbf{u}_{i}^{\prime}\mathbf{u}_{i}\) 的条件期望，因此我们需要 \(\mathbf{y}\) 和 \(\mathbf{u} = \left[ \mathbf{u}_{1}^{\prime}, \mathbf{u}_{2}^{\prime}, \cdots, \mathbf{u}_{r}^{\prime} \right]^{\prime}\) 的联合分布。

我们有下面的式子 \[ \begin{aligned} \mathbf{y} &= \mathbf{X} \boldsymbol{\beta} + \sum_{i=1}^{r} \mathbf{Z}_{i} \mathbf{u}_{i} + \mathbf{e} \\ &= \mathbf{X} \boldsymbol{\beta} + \sum_{i=0}^{r} \mathbf{Z}_{i} \mathbf{u}_{i} \end{aligned} \] 其中 \(\mathbf{u}_{0} = \mathbf{e}\) , \(q_{0} = N\) , \(\mathbf{Z}_{0} = \mathbf{I}_{n}\) 。

我们有 \[ \operatorname{cov}(\mathbf{y}, \mathbf{u}_{j}^{\prime}) = \operatorname{cov} \left( \mathbf{X} \boldsymbol{\beta} + \sum_{i=0}^{r} \mathbf{Z}_{i} \mathbf{u}_{i}, \mathbf{u}_{j}^{\prime} \right) = \mathbf{Z}_{j} \operatorname{cov} (\mathbf{u}_{j}, \mathbf{u}_{j}^{\prime}) = \sigma_{j}^{2} \mathbf{Z}_{j} \] 同时我们有 \[ \mathbf{V} = \operatorname{cov} (\mathbf{y}) = \sum_{i=0}^{r} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \sigma_{i}^{2} = \sum_{i=1}^{r} \mathbf{Z}_{i} \mathbf{Z}_{i}^{\prime} \sigma_{i}^{2} + \sigma_{0}^{2} \mathbf{I}_{N} \] 因此 \(\mathbf{y}\) 和 \(\mathbf{u}\) 的联合分布为 \(N(\boldsymbol{\mu}, \boldsymbol{\Sigma})\) ，其中 \[ \boldsymbol{\mu} =\left[ \begin{array} {c} \mathbf{X}\boldsymbol{\beta} \\ \mathbf{0} \\ \vdots \\ \mathbf{0}\\ \end{array}\right] \]

\[ \boldsymbol{\Sigma} =\left[ \begin{array} {cc} \mathbf{V} & \left\{_{r}\sigma_{i}^{2}\mathbf{Z}_{i} \right\}_{i=1}^{r} \\ \left\{_{c}\sigma_{i}^{2}\mathbf{Z}_{i}^{\prime} \right\}_{i=1}^{r} & \left\{_{d} \sigma_{i}^{2}\mathbf{I}_{q_{i}} \right\}_{i=1}^{r}\\ \end{array}\right] \]

因此我们有概率密度函数 \[ f_{\mathbf{y}, \boldsymbol{\mu}}(\mathbf{y}, \boldsymbol{\mu})= (2\pi)^{-\frac{1}{2} \sum_{i=0}^{r}q_{i}} |\mathbf{\Sigma}|^{-\frac{1}{2}} \exp \left\{- \frac{1}{2} \mathbf{Q} \right\} \] 其中 \[ \mathbf{Q} = \left[ (\mathbf{y-X}\boldsymbol{\beta})^{\prime} \quad \mathbf{u}_{1}^{\prime} \quad \mathbf{u}_{2}^{\prime} \quad \cdots \quad \quad \mathbf{u}_{r}^{\prime} \right] \mathbf{\Sigma}^{-1} \left[\begin{array}{c} (\mathbf{y-X}\boldsymbol{\beta}) \\ \mathbf{u}_{1} \\ \mathbf{u}_{2} \\ \vdots \\ \mathbf{u}_{r} \end{array} \right] \] 根据舒尔补，我们有公式 $$ | \[\begin{array}{c} \mathbf{A} & \mathbf{B}\\ \mathbf{C} & \mathbf{D}\\ \end{array}\] | = ||| \[ 我们推导得到 \] \[\begin{aligned} |\mathbf{\Sigma}| &= |\left\{_{d} \sigma_{i}^{2}\mathbf{I}_{q_{i}} \right\}_{i=1}^{r}|\mathbf{V} - \left\{_{r}\sigma_{i}^{2}\mathbf{Z}_{i} \right\}_{i=1}^{r} \left\{_{d} (\sigma_{i}^{2})^{-1}\mathbf{I}_{q_{i}} \right\}_{i=1}^{r} \left\{_{c}\sigma_{i}^{2}\mathbf{Z}_{i}^{\prime} \right\}_{i=1}^{r} | \\ &=\prod_{i=1}^{r}(\sigma_{i}^{2})^{q_{i}} |\mathbf{V} - \sum_{i=1}^{r}\sigma_{i}^{2}\mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime}| \\ &=\prod_{i=1}^{r}(\sigma_{i}^{2})^{q_{i}} |\sigma_{0}^{2}\mathbf{I}_{N}| \\ &=\prod_{i=1}^{r}(\sigma_{i}^{2})^{q_{i}} (\sigma_{0}^{2})^{N} \\ &=\prod_{i=0}^{r}(\sigma_{i}^{2})^{q_{i}} \\ \end{aligned}\] \[ 同时根据舒尔补公式 \] ^{-1} = + (-^{-1} )^{-1}\[ 我们得到 \] ^{-1} = + {0}^{-2}{N} \[ 因此 \] \[\begin{aligned} \mathbf{Q} &= \left[ (\mathbf{y-X}\boldsymbol{\beta})^{\prime} \quad \mathbf{u}^{\prime} \right] \mathbf{\Sigma}^{-1} \left[\begin{array}{c} (\mathbf{y-X}\boldsymbol{\beta}) \\ \mathbf{u} \\ \end{array} \right] \\ &= \sum_{i=1}^{r} \frac{\mathbf{u}_{i}^{\prime}\mathbf{u}_{i}}{\sigma_{i}^{2}} + (\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i})^{\prime} (\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i})/\sigma_{0}^{2} \\ \end{aligned}\] \[ 因此，将上面的式子代入基于完全数据的对数似然函数（概率密度函数），得到 \] \[\begin{aligned} l & = - \frac{1}{2} \left( \sum_{i=0}^{r} q_{i} \right) \ln 2 \pi- \frac{1}{2} \sum_{i=0}^{r} q_{i} \ln \sigma_{i}^{2} - \frac{1}{2}\sum_{i=1}^{r} \frac{\mathbf{u}_{i}^{\prime}\mathbf{u}_{i}}{\sigma_{i}^{2}} \\ & \quad - \frac{1}{2} (\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i})^{\prime} (\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i})/\sigma_{0}^{2} \\ & = - \frac{1}{2} \left( \sum_{i=0}^{r} q_{i} \right) \ln 2 \pi- \frac{1}{2} \sum_{i=0}^{r} q_{i} \ln \sigma_{i}^{2} - \frac{1}{2}\sum_{i=0}^{r} \frac{\mathbf{u}_{i}^{\prime}\mathbf{u}_{i}}{\sigma_{i}^{2}} \\ \end{aligned}\]

$$ 其中第二步推导是因为 \(\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i} = \mathbf{e} = \mathbf{u}_{0}\) 。

因此根据完全数据的对数似然函数，对 \(\sigma_{i}^{2}\) 求偏导，得到 \[ \frac{\partial l}{\partial \sigma_{i}^{2}} = - \frac{1}{2} \left( \frac{q_{i}}{\sigma_{i}^{2}} - \frac{\mathbf{u}_{i}^{\prime}\mathbf{u}_{i}}{\sigma_{i}^{4}} \right) \] 使该式为 0 ，我们可以轻松得到 ML 估计值 \[ \hat{\sigma}_{i}^{2} = \mathbf{u}_{i}^{\prime}\mathbf{u}_{i}/q_{i} \] 对 \(\boldsymbol{\beta}\) 求偏导，计算过程如下 \[ \begin{aligned} &\mathrm{d}(\mathrm{tr}((\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i})^{\prime} (\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i}))) \\ &=\mathrm{tr}(\mathrm{d}((\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i})^{\prime} (\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i}))) \\ &=\mathrm{tr}\left(-\mathrm{d}(\boldsymbol{\beta}^{\prime}\mathbf{X}^{\prime}) (\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i})-(\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i})^{\prime}\mathrm{d}(\mathbf{X} \boldsymbol{\beta}) \right) \\ &=-\mathrm{tr}\left(\mathrm{d}(\boldsymbol{\beta}^{\prime}) \mathbf{X}^{\prime} (\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i})\right)-\mathrm{tr}\left((\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i})^{\prime}\mathbf{X}\mathrm{d}( \boldsymbol{\beta}) \right) \\ &=-\mathrm{tr}\left((\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i})^{\prime}\mathbf{X}\mathrm{d}( \boldsymbol{\beta}) \right)-\mathrm{tr}\left((\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i})^{\prime}\mathbf{X}\mathrm{d}( \boldsymbol{\beta}) \right) \\ &=-2\mathrm{tr}\left((\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i})^{\prime}\mathbf{X}\mathrm{d}( \boldsymbol{\beta}) \right) \\ \end{aligned} \] 因此 \[ \frac{\partial l}{\partial \boldsymbol{\beta}} = \frac{1}{\sigma_{0}^{2}} \mathbf{X}^{\prime} (\mathbf{y-X}\boldsymbol{\beta}-\sum_{i=1}^{r}\mathbf{Z}_{i}\mathbf{u}_{i}) \] 使之为 \(\mathbf{0}\) ，得到 \[ \boldsymbol{\hat{\beta}} = \mathbf{(X'X)^{-}X'}\left( \mathbf{y} - \sum_{i=1}^{r} \mathbf{Z}_{i} \mathbf{u}_{i} \right) \] 即 \[ \mathbf{X}\boldsymbol{\hat{\beta}} = \mathbf{X(X'X)^{-}X'}\left( \mathbf{y} - \sum_{i=1}^{r} \mathbf{Z}_{i} \mathbf{u}_{i} \right) \] 这个公式等价于从 \(\mathbf{y}\) 中剔除除了残差之外的随机效应，再应用普通最小二乘法得到的估计值。

为了完善 EM 算法的迭代步骤，我们还需要知道 \(\mathbf{u}_{i}^{\prime}\mathbf{u}_{i}\) 和 \(\mathbf{y} - \sum_{i=1}^{r} \mathbf{Z}_{i} \mathbf{u}_{i}\) 在给定 \(\mathbf{y}\) 时的条件期望。我们可以直接使用多变量正态分布的公式 \[ \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} \] 进行推导得到 \[ \mathbf{u}_{i} | \mathbf{y} \sim \mathcal{N}\left[ \sigma_{i}^{2} \mathbf{Z}_{i}^{\prime}\mathbf{V}^{-1}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}), \sigma_{i}^{2}\mathbf{I}_{q_{i}} - \sigma_{i}^{4} \mathbf{Z}_{i}^{\prime}\mathbf{V}^{-1} \mathbf{Z}_{i}\right] \] 因此我们得到 \[ \mathrm{E}(\mathbf{u}_{i} | \mathbf{y}) = \sigma_{i}^{2} \mathbf{Z}_{i}^{\prime}\mathbf{V}^{-1}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}) \] 根据二次型期望公式 \(\mathrm{E}\left(\mathbf{y}^{\prime} \mathbf{A y}\right)=\operatorname{tr}(\mathbf{A} \Sigma)+\boldsymbol{\mu}^{\prime} \mathbf{A} \boldsymbol{\mu}\) ，我们可以推导得到 \[ \mathrm{E}(\mathbf{u}_{i}^{\prime}\mathbf{u}_{i}|\mathbf{y}) = \sigma_{i}^{4}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta})^{\prime}\mathbf{V}^{-1}\mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime}\mathbf{V}^{-1}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}) + \mathrm{tr}(\sigma_{i}^{2}\mathbf{I}_{q_{i}} - \sigma_{i}^{4} \mathbf{Z}_{i}^{\prime}\mathbf{V}^{-1} \mathbf{Z}_{i}) \]

ML 的 EM 算法1

我们现在可以形成 ML 的 EM 算法的正式描述，其中上标 \(m\) 为第 m 次迭代得到的估计值

step 0: 给定一组初始值 \(\boldsymbol{\beta}^{(0)}\) 和 \(\mathbf{\sigma}^{2(0)}\) ，\(m=0\) 。

step 1 (E 步): 计算 \(\mathbf{u}_{i}^{\prime}\mathbf{u}_{i}\) 和 \(\mathbf{y} - \sum_{i=1}^{r} \mathbf{Z}_{i} \mathbf{u}_{i}\) 这两个统计量的条件期望，将它们标记为 \(t_{i}\) 和 \(\mathbf{s}\) 。 \[ \begin{aligned} \hat{t}_{i}^{(m)} &= \mathrm{E}(\mathbf{u}_{i}^{\prime}\mathbf{u}_{i}|\mathbf{y})|_{\boldsymbol{\beta}=\boldsymbol{\beta}^{(m)}, \mathbf{\sigma}^{2}=\mathbf{\sigma}^{2(m)}} \\ &=\sigma_{i}^{4(m)}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}^{(m)})^{\prime}(\mathbf{V}^{(m)})^{-1}\mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime}(\mathbf{V}^{(m)})^{-1}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}^{(m)}) \\ & \quad + \mathrm{tr}(\sigma_{i}^{2(m)}\mathbf{I}_{q_{i}} - \sigma_{i}^{4(m)} \mathbf{Z}_{i}^{\prime}(\mathbf{V}^{(m)})^{-1} \mathbf{Z}_{i}) \\ \end{aligned} \]

\[ \begin{aligned} \mathbf{\hat{s}}^{(m)} &= \mathrm{E}(\mathbf{y} - \sum_{i=1}^{r} \mathbf{Z}_{i} \mathbf{u}_{i}|\mathbf{y})|_{\boldsymbol{\beta}=\boldsymbol{\beta}^{(m)}, \mathbf{\sigma}^{2}=\mathbf{\sigma}^{2(m)}} \\ &= \mathbf{y} - \sum_{i=1}^{r} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime} \sigma_{i}^{2(m)} (\mathbf{V}^{(m)})^{-1}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}^{(m)}) \\ &= \mathbf{y} - (\mathbf{V}^{(m)} - \sigma_{0}^{2(m)}\mathbf{I}) (\mathbf{V}^{(m)})^{-1}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}^{(m)}) \\ &= \mathbf{y} - (\mathbf{V}^{(m)} - \sigma_{0}^{2(m)}\mathbf{I}) (\mathbf{V}^{(m)})^{-1}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}^{(m)}) \\ &= \mathbf{X}\boldsymbol{\beta}^{(m)}+\sigma_{0}^{2(m)} (\mathbf{V}^{(m)})^{-1} (\mathbf{y}-\mathbf{X}\boldsymbol{\beta}^{(m)}) \\ \end{aligned} \]

step 2 (M 步): 基于完整数据，采用最大似然法，得到 \[ \sigma_{i}^{2(m+1)} = \hat{t}_{i}^{(m)}/q_{i}, \quad i=0,1,2,\cdots,r \]

\[ \mathbf{X}\boldsymbol{\hat{\beta}^{(m+1)}} = \mathbf{X(X'X)^{-}X'}\mathbf{\hat{s}}_{(m)} \]

step 3: 如果达到了收敛标准，则将 \(\boldsymbol{\hat{\sigma}}^{2} = \boldsymbol{\sigma}^{2(m+1)}\) ， \(\boldsymbol{\hat{\beta}} = \boldsymbol{\beta}^{(m+1)}\) ；不然 m 递增1，返回 step 1 。

ML 的 EM 算法2

如果我们使用 \(\mathbf{V}\) 的 ML 估计值 \(\mathbf{\hat{V}}\) ，则 \(\mathbf{X} \boldsymbol{\beta}\) 的 ML 估计值为 \(\mathbf{X(X'\hat{V}^{-1}X)^{-}X'\hat{V}^{-1}y}\) 。

Laird (1992) 建议在 EM 算法的迭代过程中不计算 \(\mathbf{\hat{s}}^{(m)}\) 和 \(\boldsymbol{\beta}^{(m)}\) ，而只是收敛时计算一次 \(\boldsymbol{\hat{\beta}}\) ，推导过程如下，如果上面式子中的 \(\mathbf{X} \boldsymbol{\beta}^{(m)}\) 均采用 \(\mathbf{X} \boldsymbol{\hat{\beta}}\) 的公式计算，即 \[ \mathbf{X} \boldsymbol{\beta}^{(m)}_{*} = \mathbf{X(X'(V^{(m)})^{-1}X)^{-}X'(V^{(m)})^{-1}y} \] 而在计算 \(\hat{t}_{i}^{(m)}\) 其中用到了 \((\mathbf{V}^{(m)})^{-1}(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}^{(m)})\) ，这正好是 \(\mathbf{P}^{(m)}\mathbf{y}\) ，其中 \(\mathbf{P}^{(m)}\) 为： \[ \mathbf{P^{(m)}} = (\mathbf{V^{(m)}})^{-1} - (\mathbf{V^{(m)}})^{-1} \mathbf{X(X'(V^{(m)})^{-1}X)^{-}X'(V^{(m)})^{-1}} \] 因此此时在每个迭代过程中，我们不需要得到 \(\boldsymbol{\beta}^{(m)}\) ，我们只需要计算 \(\mathbf{P}^{(m)}\mathbf{y}\) 。

这个算法其实严格来说不是 EM 算法，但是和 EM 算法差距很小，具体过程如下：

step 0: 给定一组初始 \(\mathbf{\sigma}^{2(0)}\) ，\(m=0\) 。

step 1 (E 步): 计算 \(\mathbf{u}_{i}^{\prime}\mathbf{u}_{i}\) 的条件期望，标记为 \(t_{i}\) 。 \[ \begin{aligned} \hat{t}_{i}^{(m)} &= \mathrm{E}(\mathbf{u}_{i}^{\prime}\mathbf{u}_{i}|\mathbf{y})|_{\boldsymbol{\beta}=\boldsymbol{\beta}^{(m)}, \mathbf{\sigma}^{2}=\mathbf{\sigma}^{2(m)}} \\ &=\sigma_{i}^{4(m)}\mathbf{y}^{\prime}\mathbf{P}^{(m)}\mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime}\mathbf{P}^{(m)}\mathbf{y} \\ & \quad + \mathrm{tr}(\sigma_{i}^{2(m)}\mathbf{I}_{q_{i}} - \sigma_{i}^{4(m)} \mathbf{Z}_{i}^{\prime}(\mathbf{V}^{(m)})^{-1} \mathbf{Z}_{i}) \\ \end{aligned} \] step 2 (M 步): 基于完整数据，采用最大似然法，得到 \[ \sigma_{i}^{2(m+1)} = \hat{t}_{i}^{(m)}/q_{i}, \quad i=0,1,2,\cdots,r \] step 3: 如果达到了收敛标准，则将 \(\boldsymbol{\hat{\sigma}}^{2} = \boldsymbol{\sigma}^{2(m+1)}\) ， \(\mathbf{X} \boldsymbol{\hat{\beta}} = \mathbf{X(X'\hat{V}^{-1}X)^{-}X'\hat{V}^{-1}y}\) ；不然 m 递增1，返回 step 1 。

EM 算法和 ML 方程组的等价性

我们现在考虑EM 算法和 ML 方程组的等价性，首先当 EM 算法收敛时，我们有 \(\boldsymbol{\hat{\beta}}=\boldsymbol{\beta}^{(m+1)}=\boldsymbol{\beta}^{(m)}\)， \(\boldsymbol{\hat{\sigma}^{2}}=\boldsymbol{\sigma}^{2(m+1)}=\boldsymbol{\sigma}^{2(m)}\) ，利用第一个等式，根据 ML 的 EM 算法1， 我们得到 \[ \mathbf{X}\boldsymbol{\hat{\beta}} = \mathbf{X(X'X)^{-}X'}\left[\mathbf{X}\boldsymbol{\hat{\beta}}+\hat{\sigma}_{0}^{2} (\mathbf{\hat{V}}^{-1} (\mathbf{y}-\mathbf{X}\boldsymbol{\hat{\beta}}) \right] \] 根据广义逆性质我们有 \(\mathbf{X} = \mathbf{X(X'X)^{-}X'X}\) ，因此化简得到 \[ \mathbf{X(X'X)^{-}X'} \mathbf{\hat{V}}^{-1}\mathbf{X}\boldsymbol{\hat{\beta}} = \mathbf{X(X'X)^{-}X'} \mathbf{\hat{V}}^{-1}\mathbf{y} \] 等式左右两侧乘以 \(\mathbf{X}^{\prime}\) , 根据广义逆性质 \(\mathbf{X'} = \mathbf{X'X(X'X)^{-}X'}\), 我们得到 \[ \mathbf{X'}\mathbf{\hat{V}}^{-1}\mathbf{X}\boldsymbol{\hat{\beta}} = \mathbf{X'} \mathbf{\hat{V}}^{-1}\mathbf{y} \] 这就是 ML 方程组中的一个式子。

类似地，对于 \(\hat{\sigma}_{i}^{2}\) ，我们有 \[ \begin{aligned} q_{i} \hat{\sigma}_{i}^{2} &= \hat{\sigma}_{i}^{4}\mathbf{y}^{\prime}\mathbf{\hat{P}}\mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime}\mathbf{\hat{P}}\mathbf{y} + \mathrm{tr}(\hat{\sigma}_{i}^{2}\mathbf{I}_{q_{i}} - \hat{\sigma}_{i}^{4} \mathbf{Z}_{i}^{\prime}(\mathbf{\hat{V}})^{-1} \mathbf{Z}_{i}) \\ &= \hat{\sigma}_{i}^{4}\mathbf{y}^{\prime}\mathbf{\hat{P}}\mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime}\mathbf{\hat{P}}\mathbf{y} + q_{i}\hat{\sigma}_{i}^{2} - \hat{\sigma}_{i}^{4} \mathrm{tr}( \mathbf{Z}_{i}^{\prime}(\mathbf{\hat{V}})^{-1} \mathbf{Z}_{i}) \\ \end{aligned} \] 只要 \(\hat{\sigma}_{i}^{2} \neq 0\) ，我们可以简单得到 \[ \mathrm{tr}( \mathbf{\hat{V}}^{-1} \mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime}) = \mathbf{y}^{\prime}\mathbf{\hat{P}}\mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime}\mathbf{\hat{P}}\mathbf{y} \] 这是 ML 方程组中的另外一个式子。

REML 的 EM 算法

我们只需要将 ML 的 EM 算法中的矩阵和向量进行以下替换： \[ \begin{aligned} &\text{replacing} \quad \mathbf{y} \quad \text{by} \quad \mathbf{K'y} \\ &\text{replacing} \quad \mathbf{X} \quad \text{by} \quad \mathbf{K'X=0} \\ &\text{replacing} \quad \mathbf{Z} \quad \text{by} \quad \mathbf{K'Z} \\ &\text{replacing} \quad \mathbf{V} \quad \text{by} \quad \mathbf{K'VK} \\ \end{aligned} \] 然后我们就可以从 ML 的 EM 算法2 推导得到 REML 的 EM 算法。

其中会应用到下面的式子 \[ \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} \] 我们先推导一下ML 的 EM 算法2式子中的 \(\mathbf{P}\) 矩阵，将其中的 \(\mathbf{X}\) 和 \(\mathbf{V}\) 矩阵进行替换，得到： \[ \begin{aligned} \mathbf{P} &=\mathbf{(K'VK)}^{-1}-\mathbf{(K'VK)}^{-1} \mathbf{(K'X)}\left(\mathbf{(K'X)}^{\prime} \mathbf{(K'VK)}^{-1} \mathbf{(K'X)}\right)^{-} \mathbf{(K'X)}^{\prime} \mathbf{(K'VK)}^{-1} \\ &=\mathbf{(K'VK)}^{-1} \quad \because \mathbf{K'X=0} \end{aligned} \] 因此ML 的 EM 算法2式子中的 \(\mathbf{Z_{i}^{\prime}Py = Z_{i}^{\prime}K(K'VK)^{-1}K'y }\) ，因此其在 REML 的算法中仍为 \(\mathbf{Z_{i}^{\prime}Py}\) 不变。

因此，我们得到 REML 的 EM 算法步骤如下：

step 0: 给定一组初始 \(\mathbf{\sigma}^{2(0)}\) ，\(m=0\) 。

step 1 (E 步): 计算 \(\mathbf{u}_{i}^{\prime}\mathbf{u}_{i}\) 的条件期望，标记为 \(t_{i}\) 。 \[ \begin{aligned} \hat{t}_{i}^{(m)} &= \mathrm{E}(\mathbf{u}_{i}^{\prime}\mathbf{u}_{i}|\mathbf{y})|_{\boldsymbol{\beta}=\boldsymbol{\beta}^{(m)}, \mathbf{\sigma}^{2}=\mathbf{\sigma}^{2(m)}} \\ &=\sigma_{i}^{4(m)}\mathbf{y}^{\prime}\mathbf{P}^{(m)}\mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime}\mathbf{P}^{(m)}\mathbf{y} \\ & \quad + \mathrm{tr}(\sigma_{i}^{2(m)}\mathbf{I}_{q_{i}} - \sigma_{i}^{4(m)} \mathbf{Z}_{i}^{\prime} \mathbf{K}(\mathbf{K}^{\prime}\mathbf{V}^{(m)}\mathbf{K})^{-1} \mathbf{K}^{\prime}\mathbf{Z}_{i}) \\ &=\sigma_{i}^{4(m)}\mathbf{y}^{\prime}\mathbf{P}^{(m)}\mathbf{Z}_{i}\mathbf{Z}_{i}^{\prime}\mathbf{P}^{(m)}\mathbf{y} + \mathrm{tr}(\sigma_{i}^{2(m)}\mathbf{I}_{q_{i}} - \sigma_{i}^{4(m)} \mathbf{Z}_{i}^{\prime} \mathbf{P}^{(m)}\mathbf{Z}_{i}) \\ \end{aligned} \] step 2 (M 步): 基于完整数据，采用最大似然法，得到 \[ \sigma_{i}^{2(m+1)} = \hat{t}_{i}^{(m)}/q_{i}, \quad i=0,1,2,\cdots,r \] step 3: 如果达到了收敛标准，则将 \(\boldsymbol{\hat{\sigma}}^{2} = \boldsymbol{\sigma}^{2(m+1)}\) ；不然 m 递增1，返回 step 1 。

我们可以看到 REML 的 EM 算法与 ML 的 EM 算法只有一点不同，就是将公式中的 \(\mathbf{V}^{-1}\) 替换成了 \(\mathbf{P}\) 。

参考文献

1. Rencher A C, Schaalje G B. Linear models in statistics[M]. John Wiley & Sons, 2008.
2. Patterson H D, Thompson R. Recovery of inter-block information when block sizes are unequal[J]. Biometrika, 1971, 58(3): 545-554.
3. Hofer A. Variance component estimation in animal breeding: a review[J]. Journal of Animal Breeding and Genetics, 1998, 115(1‐6): 247-265.
4. Misztal I. Reliable computing in estimation of variance components[J]. Journal of animal breeding and genetics, 2008, 125(6): 363-370.
5. 张沅, 张勤, 家畜育种. 畜禽育种中的线性模型[M]. 北京农业大学出版社, 1993.
6. Searle S R. Large sample variances of maximum likelihood estimators of variance components using unbalanced data[J]. Biometrics, 1970: 505-524.
7. Searle S R. Notes on variance components estimation-A detailed account of maximum likelihood and kindred methodology[J]. Technical Report BU-673-M, 1979.
8. Searle S R. An overview of variance component estimation[J]. Metrika, 1995, 42(1): 215-230.
9. Searle S R, Casella G, McCulloch C E. Variance components[M]. John Wiley & Sons, 2009.
10. Matsaglia G, PH Styan G. Equalities and inequalities for ranks of matrices[J]. Linear and multilinear Algebra, 1974, 2(3): 269-292.
11. Marsaglia G, Styan G P H. Rank conditions for generalized inverses of partitioned matrices[J]. Sankhyā: The Indian Journal of Statistics, Series A, 1974: 437-442.
12. Bard, Y. (1974). Nonlinear Parameter Estimation. Academic Press, New York.