这是方差组分估计方法的第五篇博客，也是最后一篇，介绍 AI-REML 算法。

AI-REML

首先我们看 Newton-Raphson 算法，其迭代公式为： \[ \boldsymbol{\theta}^{(t+1)}=\boldsymbol{\theta}^{(t)}-\left(\mathbf{H}^{(t)}\right)^{-1} \frac{\partial L}{\partial \boldsymbol{\theta}} \mid \boldsymbol{\theta}^{(t)} \] 其中 \(L\) 是需要最大化的原函数；\(\frac{\partial L}{\partial \boldsymbol{\theta}}\) 是 \(L\) 函数关于参数 \(\boldsymbol{\theta}\) 的一阶偏导数向量； \(\mathbf{H}\) 是 \(L\) 函数关于参数 \(\boldsymbol{\theta}\) 的二阶偏导数矩阵，即 Hessian 矩阵，举例如下 \[ \boldsymbol{H}=\left[\begin{array}{cccc} \frac{\partial^{2} L}{\partial \theta_{1}^{2}} & \frac{\partial^{2} L}{\partial \theta_{1} \partial \theta_{2}} & \cdots & \frac{\partial^{2} L}{\partial \theta_{1} \partial \theta_{k}} \\ \frac{\partial^{2} L}{\partial \theta_{1} \partial \theta_{2}} & \frac{\partial^{2} L}{\partial \theta_{2}^{2}} & \cdots & \frac{\partial^{2} L}{\partial \theta_{2} \partial \theta_{k}} \\ \vdots & \vdots & & \vdots \\ \frac{\partial^{2} L}{\partial \theta_{1} \partial \theta_{k}} & \frac{\partial^{2} L}{\partial \theta_{2} \partial \theta_{k}} & \cdots & \frac{\partial^{2} L}{\partial \theta_{k}^{2}} \end{array}\right] \] 由于计算 \(\mathbf{H}\) 矩阵的计算量较大，scoring 方法采用以下迭代公式，即采用 \(\mathbf{F}\) 矩阵替换 \(\mathbf{H}\) 矩阵。 \(\mathbf{F}\) 矩阵是 \(\mathbf{H}\) 矩阵的数学期望，即 \(\boldsymbol{F}=E(\boldsymbol{H})\) ； \(- \mathbf{F}\) 矩阵则称为 Fisher 信息矩阵 \[ \boldsymbol{\theta}^{(t+1)}=\boldsymbol{\theta}^{(t)}-\left(\mathbf{F}^{(t)}\right)^{-1} \frac{\partial L}{\partial \boldsymbol{\theta}} \mid \boldsymbol{\theta}^{(t)} \] 所谓平均信息算法，就是将期望信息 (\(- \mathbf{F}\)) 和观察信息 (\(- \mathbf{H}\)) 平均起来，得到一个平均信息矩阵 \[ A \boldsymbol{I}=\frac{-\boldsymbol{F}+(-\boldsymbol{H})}{2} \] 此时相应的迭代公式为 \[ \boldsymbol{\theta}^{(t+1)}=\boldsymbol{\theta}^{(t)}+\left(\mathbf{AI}^{(t)}\right)^{-1} \frac{\partial L}{\partial \boldsymbol{\theta}} \mid \boldsymbol{\theta}^{(t)} \] 在前面的推导中，我们已知 REML 方法的似然函数为 \[ \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} \] 之前我们推导证明： \[ \frac{\partial^2 l_{R}}{\partial \sigma_i^2 \partial \sigma_j^2} =\frac{1}{2} \operatorname{tr}\left(\mathbf{P } \mathbf{V}_j \mathbf{P } \mathbf{V}_i\right)-\mathbf{y}^{\prime} \mathbf{P} \mathbf{V}_j \mathbf{P} \mathbf{V}_i \mathbf{P y} \] 举个例子，假设模型中随机效应只有加性效应（\(\mathbf{G} = \mathbf{A} \sigma^{2}_{a}\) ）和残差 (\(\mathbf{R} = \mathbf{I} \sigma^{2}_{e}\))，此时 \(\mathbf{H}\) 矩阵为 \[ \begin{aligned} \mathbf{H} &=\left[\begin{array}{cc} \frac{\partial^{2} L}{\partial\left(\sigma_{a}^{2}\right)^{2}} & \frac{\partial \cdot}{\partial \sigma_{a}^{2} \partial \sigma_{e}^{2}} \\ \frac{\partial^{2} L}{\partial \sigma_{\mathrm{a}}^{2} \partial \sigma_{e}^{2}} & \frac{\partial^{2} L}{\partial\left(\sigma_{e}^{2}\right)^{2}} \end{array}\right] \\ &\\ &=\frac{1}{2}\left[\begin{array}{cc} \operatorname{tr}\left(\mathbf{P} \mathbf{V}_{a} \mathbf{P} \mathbf{V}_{a}\right)-2 \mathbf{y}^{\prime} \mathbf{P} \mathbf{V}_{a} \mathbf{P} \mathbf{V}_{a} \mathbf{P} \mathbf{y} & \operatorname{tr}\left(\mathbf{P} \mathbf{V}_{a} \mathbf{P}\right)-2 \mathbf{y}^{\prime} \mathbf{P} \mathbf{V}_{a} \mathbf{P P} \mathbf{y} \\ \operatorname{tr}\left(\mathbf{P} \mathbf{V}_{a} \mathbf{P}\right)-2 \mathbf{y}^{\prime} \mathbf{P} \mathbf{V}_{a} \mathbf{P P} \mathbf{y} & \operatorname{tr}(\mathbf{P P})-2 \mathbf{y}^{\prime} \mathbf{P P P} \mathbf{y} \end{array}\right] \end{aligned} \] 其中 \(\mathbf{V}_{a} = \mathbf{ZAZ'}\) 。

同样地，我们之前推导得到： \[ -\mathrm{E}\left(\frac{\partial^2 l_{\mathbf{R}}}{\partial \sigma_i^2 \sigma_j^2}\right) =\frac{1}{2} \operatorname{tr}\left(\mathbf{P} \mathbf{V}_j \mathbf{P V}_i \right) \] 采用上面的例子，此时 \(\mathbf{F}\) 矩阵为 \[ \mathbf{F}=\mathrm{E}(\mathbf{H})=-\frac{1}{2}\left[\begin{array}{cc} \operatorname{tr}\left(\mathbf{P} \mathbf{V}_{a} \mathbf{P} \mathbf{V} _{a}\right) & \operatorname{tr}\left(\mathbf{P} \mathbf{V}_{a} \mathbf{P}\right) \\ \operatorname{tr}(\mathbf{P V}_{a} \mathbf{P}) & \operatorname{tr}(\mathbf{P P}) \end{array}\right] \] 此时平均信息矩阵为 \[ \begin{aligned} \mathbf{AI}_{ij} &= \frac{-\mathbf{F}_{ij}+(-\mathbf{H}_{ij})}{2} \\ &= \frac{1}{2} \left(\frac{1}{2} \operatorname{tr}\left(\mathbf{P} \mathbf{V}_j \mathbf{P V}_i \right) -\frac{1}{2} \operatorname{tr}\left(\mathbf{P } \mathbf{V}_j \mathbf{P } \mathbf{V}_i\right)+\mathbf{y}^{\prime} \mathbf{P} \mathbf{V}_j \mathbf{P} \mathbf{V}_i \mathbf{P y} \right) \\ &= \frac{1}{2} \mathbf{y}^{\prime} \mathbf{P} \mathbf{V}_j \mathbf{P} \mathbf{V}_i \mathbf{P y} \\ \end{aligned} \] 采用上面的例子，此时 \(\mathbf{AI}\) 矩阵为 \[ \mathbf{AI}=\frac{-\mathbf{F}+(-\mathbf{H})}{2} = \frac{1}{2}\left[\begin{array}{cc} \mathbf{y}^{\prime} \mathbf{P} \mathbf{V}_{a} \mathbf{P} \mathbf{V}_{a} \mathbf{P} \mathbf{y} & \mathbf{y}^{\prime} \mathbf{P} \mathbf{V}_{a} \mathbf{P P} \mathbf{y} \\ \mathbf{y}^{\prime} \mathbf{P} \mathbf{V}_{a} \mathbf{P P} \mathbf{y} & \mathbf{y}^{\prime} \mathbf{P P P} \mathbf{y} \end{array}\right] \] 此时在 AI 算法中，不涉及迹函数的运算，计算难度相对低一点。

如果我们设 \(\mathbf{q}_{\kappa_{i}} = \mathbf{\dot{V}}_{i}\mathbf{Py}\) ，定义 \(\mathbf{Q} = \left[ \mathbf{q}_{\kappa_{1}} \cdots \mathbf{q}_{\kappa_{k}}\right]\) ，则我们可以得到 \[ \mathbf{AI} = \frac{1}{2} \mathbf{Q}^{\prime}\mathbf{P}\mathbf{Q} \] 对于某一个随机效应的方差组分 \(\gamma_{ij}\) ，我们有 \(\mathbf{\dot{V}}_{i j} = \mathbf{Z}_{i}\mathbf{\dot{G}}_{ij}\mathbf{Z}_{i}^{\prime}\) ，则有 $$ \[\begin{aligned} \mathbf{q}_{\gamma_{ij}} &= \mathbf{Z}_{i}\mathbf{\dot{G}}_{ij}\mathbf{Z}_{i}^{\prime} \mathbf{Py} \\ &= \mathbf{Z}_{i}\mathbf{\dot{G}}_{ij} \mathbf{G}_{i}^{-1} \mathbf{G}_{i} \mathbf{Z}_{i}^{\prime} \mathbf{Py} \\ &= \mathbf{Z}_{i}\mathbf{\dot{G}}_{ij} \mathbf{G}_{i}^{-1} \mathbf{u}_{i} \quad \because \mathbf{u}_{i} = \mathbf{G}_{i}\mathbf{Z}_{i}^{\prime} \mathbf{Py} \\ \end{aligned}\] \[ 对于残差的方差组分 $\phi_{ij}$ ，我们有 $\mathbf{\dot{V}}_{ij} = \mathbf{\dot{\Sigma}}_{ij}$ ，则有 \] \[\begin{aligned} \mathbf{q}_{\phi_{ij}} &= \mathbf{\dot{\Sigma}}_{ij} \mathbf{Py} \\ &= \mathbf{\dot{\Sigma}}_{ij} \mathbf{R^{-1}e} \quad \because \mathbf{Py} = \mathbf{R^{-1}e} \\ &= \mathbf{\dot{\Sigma}}_{ij} \mathbf{\Sigma^{-1}e} \end{aligned}\]

$$ 对于中间的 \(\mathbf{P}\) 矩阵，我们采用公式 \(\mathbf{P} = \mathbf{R^{-1}-R^{-1}WC^{-1}W^{\prime}R^{-1}}\) 进行计算。

我们也可以通过高斯消元得到 \(\mathbf{AI}\) 矩阵，首先我们构建增广矩阵如下 (用 \(\mathbf{Q}\) 代替 \(\mathbf{y}\) ) \[ \left[\begin{array}{ll} \mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{X} & \mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{Z} & \mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{Q} \\ \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{X} & \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Z}+\mathbf{G}^{-1} & \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Q} \\ \mathbf{Q}^{\prime} \mathbf{R}^{-1} \mathbf{X} & \mathbf{Q}^{\prime} \mathbf{R}^{-1} \mathbf{Z} & \mathbf{Q}^{\prime} \mathbf{R}^{-1} \mathbf{Q} \end{array}\right] = \left[\begin{array}{c} \mathbf{C} & \mathbf{a} \\ \mathbf{a}^{\prime} & \mathbf{Q}^{\prime} \mathbf{R}^{-1} \mathbf{Q} \end{array}\right] \] 将 \(\mathbf{C}\) 吸收进入 \(\mathbf{Q}^{\prime} \mathbf{R}^{-1} \mathbf{Q}\) 得到 \[ \begin{aligned} &\mathbf{Q}^{\prime} \mathbf{R}^{-1} \mathbf{Q} - \mathbf{a}^{\prime}\mathbf{C}^{-1}\mathbf{a}\\ &= \mathbf{Q}^{\prime} \mathbf{R}^{-1} \mathbf{Q} - \left[\begin{array}{ll} \mathbf{Q}^{\prime} \mathbf{R}^{-1} \mathbf{X} & \mathbf{Q}^{\prime} \mathbf{R}^{-1} \mathbf{Z} \end{array}\right] \mathbf{C}^{-1}\left[\begin{array}{ll} \mathbf{X}^{\prime} \mathbf{R}^{-1} \mathbf{Q} \\ \mathbf{Z}^{\prime} \mathbf{R}^{-1} \mathbf{Q} \\ \end{array}\right] \\ &= \mathbf{Q}^{\prime} \mathbf{R}^{-1} \mathbf{Q} - \mathbf{Q}^{\prime} \mathbf{R^{-1}XC^{-1}X^{\prime}R^{-1}} \mathbf{Q} - \mathbf{Q}^{\prime} \mathbf{R^{-1}ZC^{-1}Z^{\prime}R^{-1}} \mathbf{Q} \\ &= \mathbf{Q}^{\prime}( \mathbf{R}^{-1} - \mathbf{R^{-1}XC^{-1}X^{\prime}R^{-1}} -\mathbf{R^{-1}ZC^{-1}Z^{\prime}R^{-1}})\mathbf{Q} \\ &= \mathbf{Q}^{\prime}( \mathbf{R}^{-1} - \mathbf{R^{-1}WC^{-1}W^{\prime}R^{-1}} )\mathbf{Q} \quad \because \mathbf{W} = \left[\begin{array}{l} \mathbf{X} & \mathbf{Z} \end{array}\right] \\ &= \mathbf{Q}^{\prime} \mathbf{P} \mathbf{Q} \end{aligned} \]

第二种计算方法

我们现在逐个计算 \(\mathbf{AI}\) 矩阵的元素 \(\mathbf{AI}_{ij}=\frac{1}{2} \mathbf{y}^{\prime} \mathbf{P} \mathbf{V}_j \mathbf{P} \mathbf{V}_i \mathbf{P y}=\frac{1}{2} \mathbf{q}_{\kappa_{j}}^{\prime} \mathbf{P} \mathbf{q}_{\kappa_{i}}=\frac{1}{2} \mathbf{q}_{\kappa_{i}}^{\prime} \mathbf{P} \mathbf{q}_{\kappa_{j}}\) ，这里我们先忽略 1/2 这个常数项，推导如下 $$ \[\begin{aligned} &\mathbf{q}_{\kappa_{i}}^{\prime} \mathbf{P} \mathbf{q}_{\kappa_{j}} \\ &= \mathbf{q}_{\kappa_{i}}^{\prime} \left( \mathbf{R^{-1}-R^{-1}WC^{-1}W^{\prime}R^{-1}} \right) \mathbf{q}_{\kappa_{j}} \\ &= \mathbf{q}_{\kappa_{i}}^{\prime}\mathbf{R^{-1}} \mathbf{q}_{\kappa_{j}} - \mathbf{q}_{\kappa_{i}}^{\prime} \mathbf{R^{-1}WC^{-1}W^{\prime}R^{-1}}\mathbf{q}_{\kappa_{j}} \\ &= \mathbf{q}_{\kappa_{i}}^{\prime}\mathbf{R^{-1}} \mathbf{q}_{\kappa_{j}} - \mathbf{t}^{\prime} \mathbf{W^{\prime}R^{-1}}\mathbf{q}_{\kappa_{j}} \\ \end{aligned}\]

$$ 这里我们设 \(\mathbf{t} = \mathbf{C^{-1}W^{\prime}R^{-1}}\mathbf{q}_{\kappa_{i}}\) 。

这里的第一项容易计算，第二项我们将 \(\mathbf{q}_{\kappa_{i}}\) 替换 \(\mathbf{y}\) 构建混合模型方程组 \(\mathbf{Ct=W^{\prime}R^{-1}\mathbf{q}_{\kappa_{i}}}\) ，此时 \(\mathbf{W^{\prime}R^{-1}\mathbf{q}_{\kappa_{i}}}\) 是右手项，\(\mathbf{t}\) 是相应的解。

参考文献

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

2. Gilmour A R, Thompson R, Cullis B R. Average information REML: an efficient algorithm for variance parameter estimation in linear mixed models[J]. Biometrics, 1995: 1440-1450.

3. Smith S P, Graser H U. Estimating variance components in a class of mixed models by restricted maximum likelihood[J]. Journal of Dairy Science, 1986, 69(4): 1156-1165.