#asreml 平均信息REML

2019-07-29 本文已影响3人董八七

Gilmour AR (2019) Average information residual maximum likelihood in practice. J Anim Breed Genet 136:262–272. doi: 10.1111/jbg.12398

实际中的平均信息残差最大似然

摘要

Gilmour、Thompson和Cullis（Biometrics，1995，511440）提出了线性混合模型中有效方差参数估计的平均信息残差最大似然（REML）算法。该论文专门处理传统的方差分量模型，但该算法很快被应用于更通用的模型，并在几个REML包中实现，包括ASREML（Gilmour等人，Biometrics，2015，511440）。本文概述了与这些更通用模型相关的理论，描述了在拟合这些模型时遇到的主要问题，以及如何在ASREML软件中解决这些问题。所涉及的问题是算法实现的基本步骤，将参数保持在参数空间内、最大化稀疏性、通过使用因子分析结构避免与非结构化方差矩阵相关的问题，并处理基于标记的关系矩阵与当前工作中的奇异问题。
关键词 asreml，echidna软件，估算，混合模型，reml，方差分量

1引言

Patterson和Thompson（1971）在田间试验中恢复区组间信息的背景下设计了剩余最大似然（REML）方法。他们使用基于Fisher评分（期望信息）的迭代算法来估计方差参数。该算法被证明对较大的问题和那些方差参数多于几个的问题具有计算限制。这导致其他人使用自由导数（df）和期望最大化（em）方法，但尽管他们可以处理更大的模型，参数的数量仍然是一个问题。
Hofer（1998）对用于计算REML估计的不同算法进行了综合比较，Thompson（2009）对其进行了更新，如表1所示。这些比较显示了EM方法比DF方法的预期改进。比较还表明，大多数二阶微分方法在相对较少的迭代中收敛。
Gleeson、Gilmour和Cullis受雇于新南威尔士州农业部门担任生物统计师。Gleeson和Cullis（1987）开发了用于田间试验分析的空间方法，由Cullis和Gleeson（1991）扩展到二维。Gilmour为集团开发了软件，包括程序Reg（Gilmour，1992a）和Twod（Gilmour，1992b）。他1992年给Cullis编写了一个定制的df-reml程序，用于分析一个大型的回顾性基因型x环境数据集（Cullis，Thomson，Fisher，Gilmour，&Thompson，1996a，1996b）。当汤普森在那年晚些时候访问时，他提出了平均信息（ai）算法，该算法很快得到了实施，并且运行良好。有两种基于信息的算法可用：费希尔评分（预期信息）和牛顿拉斐逊（观测信息）。两者都涉及复杂的trace项的计算，使得它们不适合于大问题。然而，对于传统的方差模型，两个信息矩阵的平均值导致了trace项的取消和二次型的信息矩阵。汤普森（2019）详细讨论了AI REML算法的发展。
Johnson和Thompson（1995）以及Gilmour、Thompson和Cullis（1995）提出了线性混合模型中REML有效方差分量估计的AI算法。这些论文专门讨论了传统的方差分量模型。他们使用混合模型方程，并利用这些方程的稀疏性。该算法很快被应用于更通用的模型，例如多变量数据（Jensen、Mantysaari、Madsen和Thompson，1997年）。基于AI的软件开发使得REML的应用得到了快速扩展。
本文重点介绍了ASREML软件中AI算法的实现（Gilmour、Cullis、Welham和Thompson，1999；Gilmour、Gogel、Cullis和Thompson，2006、2008；Gilmour、Gogel、Cullis、Welham和Thompson，2002；Gilmour等，2015）。ASREML是在1996年根据以前项目的代码开发的，并利用了Arthur Gilmour、Brian Cullis、Robin Thompson和Sue Welham（Thompson，2019）的专业知识。作为合作的一部分，Asreml内核由Sue Welham（Payne等人，1997年）移植到Genstat，并与昆士兰初级工业部的David Butler共享，以开发SPlus（后来的R）用户界面（Butler、Cullis、Gilmour、Gogel和Thompson，2018年）。开发商的雇主，新南威尔士州农业和Rothasted Research，于2002年与VSN International（VSNI）签订合同，将ASREML商业化，并分别于2010年和2012年将其在ASREML的知识产权出售给VSNI。
Asreml是由热衷于在混合模型中应用新思想的生物测定学家开发的，用于农业研究和生产数据分析。Gilmour编写了计算机代码，他的代码开发人员工作：将可访问性扩展到其他平台，改进语法和文档，验证asreml的方法。他们和许多其他同事提供了定义专门选项的想法，以扩展ASREML的范围并记录ASREML的应用程序。因此，它处理各种各样的模型，其中方差结构是参数化的方差矩阵的两个或三个乘积之和，以适应特定的应用。例如，多元分析的剩余结构通常是单位/受试者同一性的直接产物和性状的方差矩阵；遗传成分是性状的遗传方差矩阵和基因组的遗传关系矩阵的直接产物。类型；具有相关基因型的一系列田间试验的剩余结构通常涉及一个单独的方差矩阵，用于每个试验的剩余，该矩阵具有自回归行相关矩阵和自回归列相关矩阵的直积。
ASREML在许多科学研究中被证明是有用的，在科学网中引用了2600多篇，其中乳制品和动物科学期刊超过1000篇，遗传学期刊超过700篇，食品科学、林业、生态学、植物科学和农学期刊各有200多篇。这部分是因为可以拟合的方差模型范围广泛且灵活，包括残差模型，部分是因为所使用的方法适合混合模型方程稀疏的大数据集的遗传分析。然而，在过去的十年里，已经出现了向更密集的基因组模型的转变，这导致了对asreml内核的返工。
基于混合模型方程、AI算法和利用稀疏性的其他大型问题的流行包包括DMU（Madsen和Jensen，2013）、Wombat（Meyer，2018）和Remlf90（Misztal等人，2016）。最近，Matilanen、Mantysaari、Lidauer、Stranden和Thompson（2013年）、Loh等人（2015年）和Stewart等人（2012年）表明，使用抽样方法可以减少计算工作量。对于需要密集基因组矩阵的数据，Lee和van der Werf（2006年）和Yang、Lee、Goddard和Visscher（2011年）讨论了使用加权最小二乘而不是混合模型方程的方法和软件，并在其情况下具有计算优势。
本文介绍了ASREML软件中一些关键思想的个人观点。汤普森（2019）总结了一些改进的计算机算法、改进的统计技术以及与学生合作进行的研究，这些研究加强了ASREML的理论基础和性能。在用户指南（Gilmour等人，2015年）和参考文献中可以找到对大多数想法的更彻底的处理。

2方法

2.1统计模型

ASReml适合一般的混合模型：
$\boldsymbol{y}=\boldsymbol{X} \boldsymbol{\beta}+\boldsymbol{Z} \boldsymbol{u}+\boldsymbol{\epsilon}$
其中y是观测向量，X是固定效应 $\beta$ 的设计矩阵，Z是随机效应 $u$ 的设计矩阵， $u \sim(\mathbf{0}, \boldsymbol{G}(\phi))$ ， $\boldsymbol{\epsilon} \sim(\boldsymbol{0}, \boldsymbol{R}(\boldsymbol{\theta}))$ 是残差向量， $\theta$ 是用于构造R的方差参数向量， $\phi$ 是用于构造G的方差参数向量。当R与传统混合模型中一样为 $\sigma^{2} \boldsymbol{I}$ 时， $\boldsymbol{G}(\phi)$ 可以相对于 $\sigma^{2}$ 进行估计，即所谓的 GAMMA ( $\gamma$ )尺度，其中 $\phi$ 中的比例参数定义为比率。这在早期的方案中很重要，以获得方差分量的REML估计，因为它减少了计算量，因此计算成本很高。
ASREML的设计是为了拟合R和G的更复杂的形式。首先，R和G可以分别表示为s和b部分的直和：
$\boldsymbol{R}=\oplus_{j=1}^{s} \boldsymbol{R}_{j}=\left[\begin{array}{ccccc}{\boldsymbol{R}_{1}} & {0} & {\ldots} & {0} & {0} \\ {0} & {\boldsymbol{R}_{2}} & {\ldots} & {0} & {0} \\ {\vdots} & {\vdots} & {\ddots} & {\vdots} & {\vdots} \\ {0} & {0} & {\ldots} & {\boldsymbol{R}_{s-1}} & {0} \\ {0} & {0} & {\ldots} & {0} & {\boldsymbol{R}_{s}}\end{array}\right]$
，
$\boldsymbol{G}=\oplus_{i=1}^{b} \boldsymbol{G}_{i}=\left[\begin{array}{ccccc}{\boldsymbol{G}_{1}} & {0} & {\ldots} & {0} & {0} \\ {0} & {\boldsymbol{G}_{2}} & {\ldots} & {0} & {0} \\ {\vdots} & {\vdots} & {\ddots} & {\vdots} & {\vdots} \\ {0} & {0} & {\ldots} & {\boldsymbol{G}_{b-1}} & {0} \\ {0} & {0} & {\ldots} & {0} & {\boldsymbol{G}_{b}}\end{array}\right]$
。
与每个部分相关的效应与其他部分的效应不相关。R的一般性有助于分析一系列试验，每个试验中可能有不同数量的限期和不同的方差参数。参数b与嵌入在u中的独立效应集数量有关。
其次，一组效应的方差结构可以定义为最多三个方差矩阵的直积。因此，多变量分析中残差的方差结构通常被定义为 $\boldsymbol{I} \otimes \boldsymbol{\Sigma}_{R}$ 对于将性状按顺序排列到units中的残差，以及 $\boldsymbol{\Sigma}_{G} \otimes \boldsymbol{A}$ 对于跨位点和/或跨性状相关的遗传效应。
组成矩阵可以有多种形式，包括同一矩阵、对角矩阵、遗传关系矩阵、时空相关矩阵、因子分析矩阵和一般（非结构化）矩阵。遗传关系矩阵可以是基于系谱的或基于标记的。用户甚至可以定义自己的结构（例如Jaffric、White、Thompson、Hill和Visscher，2002），因为人们不能总是预测用户的独创性。
定义设计矩阵也有灵活性。三次平滑样条（Verbyla、Cullis、Kenward和Welham，1999）和勒让德多项式通过基函数拟合，用户可以提供其他基函数。特定模型项的设计矩阵可以使用and()函数求和；这用于在双列实验中估计一般的结合能力、减少的动物模型以及估计遗传竞争效应（Bijma、Muir、Ellen、Wolf和van Arendonk，2007）。at()设计函数允许在模型中使用性状-/地点-特定项。例如，我们可能拟合某些地点的区组效应，而不是其他地点的区组效应，或者一个性状而不是另一个性状的协变量。str(...)函数允许方差结构根据需要跨多个模型项来估计加性和母本遗传效应之间的遗传协方差。
一个重要的发展是简化了方差结构的指定方式。最初的方法需要列出定义线性模型的术语（X和Z的分量），并分别为Z中的每个术语指定G的形式，这不是一个简单的有方差 $\sigma_{t}^{2} \boldsymbol{I}$ 的随机因子。这种方法是一种逻辑上的思考问题的方法，但从用户角度很难运行。它现在被Butler等人（2018）设计的功能语法所取代，用于在R中运行。每一个方差结构都由一个函数来表示，例如，我们不需要编写模型项Trait.animal，然后为Trait指定一个因子分析1结构和一个基于系谱的关系矩阵，我们只需编写xfa1(Trait).nrm(animal)。虽然Asreml和AsremlR具有相同的功能，但特定的函数名并不都相同。
这种函数方法也适用于R结构，其中方差结构的参数是数据中的变量，用于识别数据中的部分，并定义每个部分中小区的空间排列。例如，residual sat(Trial).ar1v(Row).ar1(Column)为数据文件中的每个测试中的残差建立了自回归相关结构，利用数据变量trial、row和Column. ar1v(.) (ar1(.))中的编码指定自回归方差（相关）矩阵；asreml为每个试验拟合一个方差和两个相关参数。
ASREML所涵盖的一般混合模型的特殊情况包括：

不完全区组分析，其中固定区组对比会删除处理信息，因为它们与处理不正交；区组效应假定为随机和独立的，G是对角矩阵，R是（尺度）恒等式，例如，具有三个区组因子（mu为模型截距）的品种试验的模型的默认方差结构正式规定为：

yield~mu+variety !r idv(Replicate)+idv(RowBlock)+idv(COlumnBlock)
residual idv(units)

空间分析（Gilmour、Cullis和Verbyla，1997）使用R来说明残差（例如，自回归行和自回归列、matern函数）之间的相关性，指定为

yield~mu+variety 
residual ar1v(Row).ar1(Column)

在G的定义中使用基于系谱（A）或标记（ $\boldsymbol{K}=\boldsymbol{M} \boldsymbol{M}^\prime$ ）的关系矩阵的遗传分析。例如，

birthwt~mu+BirthType+Parity+Flock !r nrm(LambId)+idv(DamId)+idv(Litter)
residual idv(units)

重复测量和多变量分析，包括G和R中的 $\boldsymbol{\Sigma} \otimes \boldsymbol{K}$ （或 $\boldsymbol{K} \otimes \boldsymbol{\Sigma}$ ）结构，其中 $\boldsymbol{\Sigma}$ 是结构化或非结构化（性状/时间/位置之间）方差矩阵，K是已知的关系矩阵。
多环境试验：我们对每个试验都有单独的 $\boldsymbol{R}_{i}$ 矩阵，并且遗传关系在不同的环境中有不同的表达，建模为跨环境协方差矩阵和遗传关系矩阵的直接乘积。例如，假设12个YrLoc中64个实验的18432个产量代表5132个基因型的谱系，一个可能的模型（mv拟合缺失值，因此空间网格是完整的）是

yield~mu+YrLoc+mv !r xfal(YrLoc).nrm(Geno)
residual sat(expt).ar1v(Col).ar1(Row)

2.2参数估算

对于给定的G和R，通过求解混合模型方程得到 $\beta$ 和 $u$ 的解：
$\left[\begin{array}{cc}{X^{\prime} \boldsymbol{R}^{-1} \boldsymbol{X}} & {\boldsymbol{X}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{Z}} \\ {\boldsymbol{Z}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{X}} & {\boldsymbol{Z}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{Z}+\boldsymbol{G}^{-1}}\end{array}\right]\left[\begin{array}{l}{\hat{\boldsymbol{\beta}}} \\ {\tilde{\boldsymbol{u}}}\end{array}\right]=\left[\begin{array}{c}{\boldsymbol{X}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{y}} \\ {\boldsymbol{Z}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{y}}\end{array}\right] \quad (1.1)$
其中 $\hat{\beta}$ 表示固定效应 $\beta$ 的估计值， $\tilde{u}$ 表示基于已知方差参数的随机效应 $u$ 的预测值。
方差参数通过使用AI算法最大化REML（Patterson和Thompson，1971）对数似然函数来估计（Gilmour等人，1995）。这需要方差参数的得分函数（一阶导数）和相应的AI矩阵（二阶导数）。
除常数外，REML对数似然函数可写为：
$\begin{aligned} l &=-\frac{1}{2}\left(\log \left|X^{\prime} \boldsymbol{H}^{-1} \boldsymbol{X}\right|+\log |\boldsymbol{H}|+v \log \sigma^{2}+\boldsymbol{y}^{\prime} \boldsymbol{P} \boldsymbol{y} / \sigma^{2}\right) \\ &=-\frac{1}{2}\left(\log |\boldsymbol{C}|+\log |\boldsymbol{R}|+\log |\boldsymbol{G}|+v \log \sigma^{2}+\boldsymbol{y}^{\prime} \boldsymbol{P} y / \sigma^{2}\right) \end{aligned}$
$\boldsymbol{H}=\boldsymbol{R}+\boldsymbol{Z} \boldsymbol{G} \boldsymbol{Z}^\prime$ ， $\boldsymbol{P}=\boldsymbol{H}^{-1}-\boldsymbol{H}^{-1} \boldsymbol{X}\left(\boldsymbol{X}^{\prime} \boldsymbol{H}^{-1} \boldsymbol{X}\right)^{-1} \boldsymbol{X}^{\prime} \boldsymbol{H}^{-1}$ ， $v$ 是n‐rank(X)的自由度，C是（1.1）混合模型方程的系数矩阵。如果R和G在方差尺度上参数化， $\sigma^{2}$ 的值为1；如果相对于剩余方差参数化， $\hat{\sigma}^{2}=\boldsymbol{y}^{\prime} \boldsymbol{P} \boldsymbol{y} / v$ 。
如果 $\boldsymbol{G}_{i}$ 是 $\boldsymbol{G}(\phi)$ 对 $\phi_{i}(i=1 : g)$ 的导数，并且相应的工作变量计算为 $\boldsymbol{y}_{i}=\boldsymbol{Z} \boldsymbol{G}_{i} \boldsymbol{G}^{-1} \tilde{\boldsymbol{u}}$ ，则得分函数可以写为 $-\left(\operatorname{tr}\left[\boldsymbol{G}_{i} \boldsymbol{G}^{-1}\right]-\operatorname{tr}\left[\boldsymbol{C}_{i} \boldsymbol{C}^{-1[Z Z]}\right]-\boldsymbol{y}^{\prime} \boldsymbol{P} \boldsymbol{y}_{i}\right) / 2$ ，其中 $\boldsymbol{C}_{i}=-\partial \boldsymbol{C} / \partial \phi_{i}=-\partial \boldsymbol{G}^{-1} / \partial \phi_{i}=\boldsymbol{G}^{-1} \boldsymbol{G}_{i} \boldsymbol{G}^{-1}$ ， $\boldsymbol{C}^{-1[Z Z]}$ 是C对应的倒数块矩阵至 $\boldsymbol{Z}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{Z}+\boldsymbol{G}^{-1}$ 。同样，如果 $\boldsymbol{R}_{j}$ 是 $\boldsymbol{R}(\boldsymbol{\theta})$ 对 $\theta_{j}(j=1:r)$ 的导数，并且相应的工作变量计算为 $y_{j}=R_{j} R^{-1} \tilde{\epsilon}$ ，其中 $\tilde{\boldsymbol{\epsilon}}=\boldsymbol{y}-\boldsymbol{X} \hat{\boldsymbol{\beta}}-Z \tilde{\boldsymbol{u}}$ ，则得分函数可以写成 $-\left(\operatorname{tr}\left[\boldsymbol{R}_{j} \boldsymbol{R}^{-1}\right]-\operatorname{tr}\left[\boldsymbol{C}_{j} \boldsymbol{C}^{-1}\right]-\boldsymbol{y}^{\prime} \boldsymbol{P} \boldsymbol{y}_{j}\right) / 2$ ，其中 $\boldsymbol{C}_{j}=-\partial \boldsymbol{C} / \partial \boldsymbol{\theta}_{j}$ 是平方和叉积之和（SSP）使用 $\boldsymbol{R}^{-1} \boldsymbol{R}_{j} \boldsymbol{R}^{-1}$ 形成的矩阵。AI矩阵的计算公式为 $\mathcal{I}=\boldsymbol{Y}^{\prime} \boldsymbol{P} \boldsymbol{Y} / 2$ ，其中 $Y=\left[Y_{G} Y_{R}\right]$ ， $Y_{G}$ 为g列 $y_{i}$ 的矩阵， $Y_{R}$ 为r列 $y_{j}$ 的矩阵。如果参数是相对于剩余方差参数化的，则类似的公式成立。

2.3算法

主要要求是求解（1.1）并找到得分计算所需的 $\boldsymbol{C}^{-1}$ 元素。通过递归完成，将（1.1）作为
$\boldsymbol{C}_{i} \boldsymbol{x}_{i}=\boldsymbol{r}_{i} \text { as }\left(\begin{array}{cc}{\boldsymbol{A}_{i}} & {\boldsymbol{b}_{i}} \\ {\boldsymbol{b}_{i}^{\prime}} & {m_{i}}\end{array}\right)\left(\begin{array}{c}{\boldsymbol{x}_{i-1}} \\ {x_{i}}\end{array}\right)=\left(\begin{array}{c}{\boldsymbol{e}_{i}} \\ {f_{i}}\end{array}\right) \quad(1.2)$
给定 $C_{i-1} x_{i-1}=r_{i-1}$ 吸收 $x_{i}$ ，其中 $\boldsymbol{C}_{i-1}=\boldsymbol{A}_{i}-\boldsymbol{b}_{i} \boldsymbol{b}_{i}^{\prime} / m_{i}$ ， $\boldsymbol{r}_{i-1}=\boldsymbol{e}_{i}-\boldsymbol{b}_{i} f_{i} / m_{i}$ 。只对 $b_{i}$ 的非零元素进行计算。
使用（1.2）递归地将 $\boldsymbol{C} \boldsymbol{x}=\boldsymbol{r}$ 的方程转换为 $\boldsymbol{C}_{A} \boldsymbol{x}=\boldsymbol{r}_{A}$ ，其中 $C_{A}$ 是下三角矩阵， $i$ 行为 $\boldsymbol{b}_{i}^{\prime} \quad m_{i} \quad \mathbf{0}^{\prime}$ ， $\boldsymbol{r}_{A}$ 的第i个元素为 $f_{i}$ 。
同样， $\boldsymbol{C}_{i}^{-1}$ 可以用分区形式编写，用于递归计算，如下所示
$\left(\begin{array}{cc}{\boldsymbol{A}_{i}} & {\boldsymbol{b}_{i}} \\ {\boldsymbol{b}_{i}^{\prime}} & {m_{i}}\end{array}\right)^{-1}=\left(\begin{array}{cc}{\boldsymbol{C}_{i-1}^{-1}} & {-\boldsymbol{C}_{i-1}^{-1} \boldsymbol{b}_{i} / m_{i}} \\ {-\boldsymbol{b}_{i}^{\prime} \boldsymbol{C}_{i-1}^{-1} / m_{i}} & {m_{i}^{-1}+\boldsymbol{b}_{i}^{\prime} \boldsymbol{C}_{i-1}^{-1} \boldsymbol{b}_{i} m_{i}^{-2}}\end{array}\right) \quad(1.3)$
计算分数所需的 $\boldsymbol{C}_{i}^{-1}$ 的所有元素都对应于C具有非零值的位置，因此我们只计算 $C_{A}$ 中非零值的值。这两个操作都是在原处按行执行的。
给定方差参数的初始值，迭代过程更新这些参数的主要步骤是：

使用用当前参数值（ $\boldsymbol{\theta}$ ）计算得出的 $\boldsymbol{R}^{-1}$ 形成SSP矩阵。形成的矩阵实际上是
$\left[\begin{array}{lll}{\boldsymbol{y}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{y}} & {\boldsymbol{y}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{X}} & {\boldsymbol{y}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{Z}} \\ {\boldsymbol{X}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{y}} & {\boldsymbol{X}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{X}} & {\boldsymbol{X}^{\prime} \boldsymbol{R}^{-1} \mathbf{Z}} \\ {\boldsymbol{Z}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{y}} & {\boldsymbol{Z}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{X}} & {\boldsymbol{Z}^{\prime} \boldsymbol{R}^{-1} \mathbf{Z}}\end{array}\right]$
但有一半以稀疏的存储形式按行存储下三角。
将从当前参数值（ $\phi$ ）形成的 $\boldsymbol{G}^{-1}$ 与 $\boldsymbol{Z}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{Z}$ 相加，
确定解方程的顺序（见下文讨论），该方程将保持高度稀疏性（除非结构与前一次迭代保持不变，当重复使用前一个顺序时），并将ssp矩阵放入该分析顺序中，
确定InFill模式，即，当较低的行（按新顺序）被吸收时，插入单元格后每行的相邻列表是必需的。也就是说，如果 $\left[\begin{array}{lll}{a^{\prime}} & {b} & {i}\end{array}\right]$ 是当前行（i）的邻接列表，则通过将其与a合并来更新b行的列表。邻接列表保留下来以便在后续迭代中重用。由于C的因子是 $\boldsymbol{D}_{A}^{1 / 2} \boldsymbol{C}_{A}$ ，其中 $D_{A}$ 是 $C_{A}$ 的对角线，此操作相当于David Johnson在其早期的aireml程序中使用的符号因子分解（Johnson&Thompson，1995年），
吸收方程，将C转换为 $C_{A}$ ，同时根据结果就地调整y列， $\boldsymbol{y}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{y}$ 变为 $\boldsymbol{y}^{\prime} \boldsymbol{P} \boldsymbol{y}$ （Gilmour等人，1995），
反解并计算残差，
形成工作变量（Y）及其交叉积（ $\boldsymbol{Y}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{Y}$ 、 $\boldsymbol{Y}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{X}$ 和 $\boldsymbol{Y}^{\prime} \boldsymbol{R}^{-1} \boldsymbol{Z}$ ），
使用 $C_{A}$ ，吸收工作变量交叉积以获得AI矩阵（ $\boldsymbol{Y}^{\prime} P \boldsymbol{Y}$ ）。如果需要预测值（Gilmour、Cullis、Welham、Gogel和Thompson，2004），给定预测设计矩阵D，这些都是在最后的迭代中使用 $C_{A}$ 来形成 $D^{\prime} C^{-1} X R^{-1} y$ 和 $\boldsymbol{D}^{\prime} \boldsymbol{C}^{-1} \boldsymbol{D}$ 。
通过将 $C_{A}$ 转换为 $\boldsymbol{C}^{-1}$ 来完成系数矩阵的稀疏逆：只计算 $C_{A}$ 的稀疏模式中存在的单元，
将稀疏逆矩阵转换为模型顺序，
计算方差分量的得分，
更新方差参数，
收敛性测试；如果未达到最大迭代次数且上次更新受到限制（见下一节），或者可能性变化大于0.01，则执行另一次迭代。
在ASREML4.2之前，没有执行符号分解，步骤5是使用链接列表执行的，根据需要插入单元格。在asreml 4.2中，步骤5使用索引（顺序存储邻接）列表，其中包含插入空间。处理长链接列表要比处理序列列表慢得多，而且这种更改将吸收某些大型（密集）作业所需的典型时间减半。最近，通过颠倒循环的顺序重新映射矩阵，使步骤10变得更快；对于一些大型密集作业，这将执行时间缩短了一半。此外，在ASREML4.2的步骤9中已经实现了并行处理，在四核PC上减少了70%的执行时间。

2.4参数空间

AI算法通过以下方式更新参数值：
$\boldsymbol{\kappa}_{i+1}=\boldsymbol{\kappa}_{i}+\mathcal{I}^{-1} \boldsymbol{s}$
$\boldsymbol{\kappa}_{i}=\left(\begin{array}{ll}{\boldsymbol{\phi}_{i}^{\prime}} & {\boldsymbol{\theta}_{i}^{\prime}}\end{array}\right)^{\prime}$ 是用于计算ai矩阵 $\mathcal{I}$ 的前一组参数值， $s$ 是分数向量。此更新有两个问题。第一个问题是，如果 $\boldsymbol{\kappa}_{i}$ 的值不接近reml解，则更新（ $\mathcal{I}^{-1} \boldsymbol{s}$ ）可能过多，因此算法无法收敛到reml解；它可能会找到局部解或没有解。我们可以通过将更新与 $\boldsymbol{\kappa}_{i}$ 中的值进行比较来评估这一点，如果判断为过度，则按0.5或0.1的比例缩放更新。
第二个问题是更新后的参数值可能不在所需的参数空间中。典型的情况是当相关在[-1,1]外时，方差分量是负的，或者方差分量矩阵是负定的。如果所拟合的模型不适合数据，当前的方差参数不接近REML解的值，或者仅仅因为采样的变化，则很可能发生这种情况。例如，一个传统的裂区分析可能有两个误差方差被标记为ErrorA（主区层）和ErrorB（子区层），这些方差是独立的。ErrorA应该大于ErrorB，但并非必须如此。就方差分量而言，在平衡情况下，如果ErrorB的预期均方为 $\sigma^{2}$ ，则ErrorA的预期均方为 $\sigma^{2}+k \sigma_{A}^{2}$ 。如果ErrorA小于误差B，则 $\sigma_{A}^{2}$ 为负。
然而，在REML估计中，将简单方差分量（如 $\sigma_{A}^{2}$ ）变为非负的情况并不少见。然后必须调整REML算法，因为应用约束意味着REML的可能性不最大。
如果在收缩上述所有更新之后，某些更新的参数值仍在可接受的参数空间之外，asreml会将这些值设置为更接近边界的可接受值****，并且如果问题再次出现，asreml最终会将参数固定为靠近边界的值。此外，如果某个校正的任何特定更新大于0.3，且方差参数更新大于因子15，则更新将被限制为当前迭代的这些限制。对于更新形式不为正定的不规则矩阵（通过逆矩阵进行测试），当似然函数的最大值与一个非正定矩阵重合时，该矩阵通过一个EM步更新，但该过程很慢地收敛到边界附近的矩阵。如果修改了任何参数，asreml将根据受约束的参数值重新计算其他参数的更新。分区 $\mathcal{I}=\left(\begin{array}{cc}{\mathcal{I}_{u u}} & {\mathcal{I}_{u f}} \\ {\mathcal{I}_{f u}} & {\mathcal{I}_{f f}}\end{array}\right)$ 和 $\boldsymbol{s}=\left(\begin{array}{c}{\boldsymbol{s}_{u}} \\ {\boldsymbol{s}_{f}}\end{array}\right)$ ，下标 $f$ 指的是约束参数，无约束参数的更新计算为 $\mathcal{I}_{u u}^{-1}\left(s_{u}-\mathcal{I}_{u f} \boldsymbol{v}\right)$ ，其中 $v$ 指应用于约束参数的实际更新。

2.5最大化稀疏度

处理大型模型最昂贵的步骤是步骤5和9；步骤3通常也很昂贵。由于AI算法只需要 $\boldsymbol{C}^{-1}$ 中与C的非零元素相对应的 $\boldsymbol{C}^{-1}$ 元素来计算 $\operatorname{tr}\left(\boldsymbol{C}_{i} \boldsymbol{C}^{-1}\right)$ 和 $\operatorname{tr}\left(\boldsymbol{C}_{j} \boldsymbol{C}^{-1}\right)$ ，因此通过保持 $\boldsymbol{C}_{A}$ 稀疏和只计算存在于 $\boldsymbol{C}_{A}$ 稀疏模式中的 $\boldsymbol{C}^{-1}$ 元素可以节省大量的精力和内存。因此，asreml在解这些方程之前对它们进行重新排序，以保持 $\boldsymbol{C}_{A}$ 稀疏，从而生成一个稀疏的 $\boldsymbol{C}^{-1}$ 。计算工作量大致与（ $\boldsymbol{C}_{A}$ 的平均行长）的平方乘以C中的行数成正比。因此， $\boldsymbol{C}_{A}$ 的大小至关重要。确定适当的订单是一个公认的问题。
ASREML中的算法随着时间的推移而发展。最初的方法（ns）只是选择在每个步骤中吸收剩余的最稀疏的行。这是通过处理列与集合相同的行来加速的（方法gns）。这些方法给出了相似的顺序。该算法在ASREML4中得到了进一步的改进，以首先清除密集块（通常是没有数据的K矩阵行），即使它们不是短的方法dgns。
这是在asreml中实现的一个昂贵的排序过程，因为它模拟了总吸收过程（模式，而不是实际的矩阵），通常在每个步骤中评估几个备选方案。
行之间的联系（通常在文献中称为边缘）来自模型拟合（交叉积），来自空间模型（相邻小区）产生的联系，以及遗传关系产生的联系。拟合模型的微小变化会对优化的加工顺序产生很大的影响。
Ducrocq和Druet（2003年）和Meyer（2005年）推荐了Metis库（Karypis，2011年）。这些算法基于图的多级K-way划分。Metis是使用初始粗化，分割，然后不粗化，导致多层，嵌套解剖涉及自上而下和自下而上的顺序。Meyer（2005）在她的软件中评估了Metis，并将其与一系列其他已发表的算法进行了比较，其中包括多重最小度常规genmmd（George&Liu，1989年）和从MUMPS包中提取的各种改进（Amestoy、Duff、L'Excellent和Koster，2001年）。她总结道，Metis是生成排序，它需要大约一半的最小度数算法执行时间，尽管时间随示例的特定特性略有不同。
Gilmour和Thompson（2006）比较了Metis与当时在ASREML中使用的算法，结果汇总在表2中。第一个例子是对26875只动物的数据进行六特征多变量SIRE模型分析，拟合31745个方程和39个方差参数。第二个例子是437632个数据记录的多元随机回归，拟合536288个方程和116个方差参数。Metis是一种比asreml程序更快地找到顺序的程序，但是它的填充更大，处理时间也更长，因此在多次迭代中没有任何净优势。该表已更新为使用ASREML3.0（建于2011年）、4.1和4.2在四核64位笔记本电脑（4）和八核PC（8）上运行的第一个示例的重新运行。新时代的一个重要组成部分是更快的硬件，因为填充值几乎没有改变。在asreml 3中，获得排序所花费的时间明显得到了很大的改进，在本例中，当前asreml中的时间得到了进一步的改进。这一改进的一部分来自于在行长度都很长的后一阶段使用一个更简单的规则。4.2结果表明，在算法中插入步骤4，并在步骤9中使用并行处理，可以获得处理增益。尽管如此，在许多大型工作中，步骤3仍然是处理时间的重要组成部分，排序问题将作为Echidna项目的一部分重新讨论（见下文）。

2.6因子分析模型

非结构化（完全参数化）方差矩阵约束为正（半）定的估计常会出现问题。针对这种情况，参数扩展EM算法（Liu、Rubin和Wu，1998）在ASREML中实现，但当最大似然值与参数空间外的参数值相对应时，不会很快产生可接受的结果。Meyer（2019）提出了一种方法，在多分析的背景下，寻求将遗传方差矩阵向表型方差矩阵的协方差结构弯曲。另一种方法是使用更简洁的参数化。因子分析方差结构由 $\boldsymbol{\Sigma}=\boldsymbol{\Gamma} \boldsymbol{\Gamma}^{\prime}+\boldsymbol{\Psi}$ 给出， $\boldsymbol{\Gamma}$ 中的列被称为荷载，将 $k$ 个变量与 $f$ 个潜在因子联系起来； $\boldsymbol{\Psi}$ 是k个特定方差的对角矩阵。这些因子类似于通常被约束为正交的主要成分。这一模型可以理解为复合对称模型的延伸，即Genotype + Site.Genotype。当 $\boldsymbol{\Gamma}$ 有一列的所有负荷相等时（即Genotype方差分量的平方根），且所有特定方差相等时（即Site.Genotype方差分量），因子分析模型等价于复合对称。该结构已证明是一种有效的完全参数化公式的替代方案，而这种公式往往过于参数化。因子分析结构很容易通过将特定的方差限制为零或正来确定正（半）定。
Thompson、Cullis、Smith和Gilmour（2003）提出了一个编写方差的模型的稀疏实现，它的逆矩阵是 $\left(\begin{array}{cc}{\Psi+\Gamma \Gamma^{\prime}} & {-\Gamma} \\ {-\Gamma^{\prime}} & {I}\end{array}\right)^{-1}=\left(\begin{array}{cc}{\Psi^{-1}} & {\Psi^{-1} \Gamma} \\ {\Gamma^{\prime} \Psi^{-1}} & {I+\Gamma^{\prime} \Psi^{-1} \Gamma}\end{array}\right)$ 。为了使用这个扩展公式，需要扩展带有k列的设计矩阵，以包含f个零的列用于因子。当某些特定的方差为零时，汤普森等人（2003）展示如何将 $\Psi$ 零行折叠成因子行，以便可以估计所有负载。
该模型的优点是：（a）如果需要，可以先用1个因子，然后用2个因子，再加上更多的因子；（b） $\Psi$ 的元素可能为零，从而得到一个可能的奇异矩阵（相关性为1）；（c）虽然有时是脆弱的，但事实分析模型比完全参数更具鲁棒性。在这些情况下，检验模型；（d）在大多数情况下，当k很大时，协方差的大部分只由少数几个因素获得，从而形成一个节俭的模型；（e）稀疏公式的运行速度大大快于具有稠密逆的方差结构。
Ganesaligam（2013）提出了一种更为稀疏的因子分析模型。它包括将负载拟合为一个完全降秩的FA模型（所有特定的方差均为零），并在单独的对角项（RR + DIAG, rr1(Site).id(Genotype) + diag(Site).id(Genotype)）中拟合特定的方差。她在对菜籽油多环境试验的分析中发现，对于两个涉及加性遗传方差效应的模型，RR+DIAG参数化比XFA参数化快10和3倍。
表3报告了ASREML4.2中关于多环境试验分析的一些最新时间比较。该模型适用于59个评价3755个基因型的环境中的AR-AR空间变异。根据需要在环境中拟合额外的随机区组因子。遗传模型 RR + D为rrk(Environment).id(Genotype) + diag(Environment).id(Genotype)和XFA为xfak(Environment).id(Genotype)，其中k取1:6。XFA1模型报告了37608个数据记录、544598个方程和482个方差分量。在这种情况下，对于有两个以上因素的情况，RR+D明显比XFA慢。这可能是排序问题，需要进一步调查。
Tolhurst、Mathews、Smith和Cullis（2019）的多环境试验数据集的进一步分析表明，与表3的结果有些相似。在这项分析中，遗传变异被建模为独立的或基因相关的基因型。对于加性和基因组关系矩阵，秩1模型的rr+diag参数化比xfa参数化快8倍。对于加性模型下的4阶模型参数化，其优势随着阶数的增加而减小，但对于4阶基因组方差模型，RR+XFA参数化仍然具有优势。
FA模型有两个问题。（a）当存在多个荷载时， $\mathcal{I}$ 矩阵中出现奇异。历史上，我们将一些加载设置为零，但实际上，最好将加载保持在接近正交的位置，但不要更新本来设置为零的加载；自从采用这种方法以来，作者对这些模型的收敛性问题较少；（b）与加载相关的 $\mathcal{I}$ 矩阵行NGS可以是非条件的，加载参数有时是强（负）相关的。这可以通过使用岭回归技术、测试不良状况以及将与载荷相关的 $\mathcal{I}$ 对角线元素膨胀1%或更多来克服。

2.7标记模型

一组个体/动物/基因型的标记数据（M）通常在许多标记位置编码为0/1/2。如果个体多于标记，通常将标记作为一组具有共同方差的效应进行拟合。当标记物的数量超过个体数量时，一般情况下最好将个体与基于 $M M^{\prime}$ 的基因组关系矩阵（K）相匹配。对于后一种情况，Asreml有一个快速的Bayes A选项（Sun、Qu、Garrick、Dekkers和Fernando，2012），对于已知的方差分量，根据已知自由度的t分布，相对于较大的效果，该选项缩小较小的效应。当目的是识别与qtl相关的标记物并且进一步发展是可能的时候，这可能是有用的。
有时，K是奇异矩阵。ASREML将在求逆时检测到这一点，并在这种情况下使用一个特殊的扩展逆，其推导如下（Kelly，2013）。首先要注意的是，如果我们为每一行K添加一个空方程，并使用 $\left(\begin{array}{cc}{0} & {1} \\ {1} & {-K}\end{array}\right)$ 形式的扩展逆矩阵，就可以拟合一个等效模型。在这个扩展的逆矩阵中，在尽可能多地吸收K行并将结果进行排列以便于解释后，我们得到了
$\left(\begin{array}{ccc}{\mathbf{0}_{s s}} & {\mathbf{0}_{s d}} & {\boldsymbol{I}_{s s}} \\ {\mathbf{0}_{d s}} & {\boldsymbol{K}_{d d}^{-1}} & {-\boldsymbol{K}_{d d}^{-1} \boldsymbol{K}_{d s}} \\ {\boldsymbol{I}_{s s}} & {-\boldsymbol{K}_{s d} \boldsymbol{K}_{d d}^{-1}} & {\mathbf{0}_{s s}}\end{array}\right)$
其中下标s属于奇异行，d属于正定行。混合模型方程需要扩展以适应第一种情况下额外的s+d行，而在简化情况下则需要扩展s行。当集合s中存在个体的数据时，这种形式的k逆函数起作用，因此当它添加到交叉积（ $\mathbf{Z}^{\prime} \boldsymbol{Z}$ ）项中时，系统不再是奇异的。多余的s行可以看作是拉格朗日方程。
如果遗传方差矩阵是线性组合，比方L关系矩阵，例如 $\boldsymbol{G}=\Sigma_{l=1}^{L} \boldsymbol{\phi}_{l} \boldsymbol{K}_{l}$ （R.Pong Wong，per.随后，Lee&van der Werf，2006）也可以使用自己的函数来拟合模型。这种参数化的优点是在混合模型方程中只需要一组效应，而不需要L。事实上，将此参数化与上述扩展模型结合使用，可以更简化计算，因为不需要G的逆矩阵，并且可以更直接地计算分数项（R.Thompson）

2.8瓦尔德F统计

有时需要对固定效果进行测试，但在混合模型中并不简单。在固定效应模型中，我们可以用一个剩余方差和自由度来划分一个平方和，得到一个ftest。在这种情况下，存在着边缘化的挑战。在混合模型中，通常不能构造一个平方和，因为我们希望考虑协方差，而且通常我们没有一个唯一的残差方差。因此，我们将估计值与指定值（有时为零）进行比较，并通过估计值的渐近方差矩阵对所得平方和矩阵进行加权。这是一个Wald式统计，一个无量纲的量，相当于固定效应模型中的F检验。它有一个渐近的齐次方分布，Kenward和Roger（1997）提出了一个具有近似f分布的比例瓦尔德统计。
这个wald f统计量的值通常取决于模型中项的拟合顺序。基本的增量wald f统计量检验每个模型项的重要性，条件是前面的项，但忽略后面的项，也就是说，按照指定的顺序，只有在所有后面的项都不重要时才被认为是有效的；如果模型项的检验是无效的，则边缘到显著的高阶项。为了测试所有使用增量F统计的可能性，通常需要在多次运行中对模型项进行明智的重新排序。
asreml报告了第二个wald f统计量，称为最大条件增量（mci）f统计量。每个术语都是以所有其他术语为条件进行测试的，除了边缘术语。例如，拟合模型项mu+a+b+c+d+a.b+a.c+a.d+b.c+b.d+a.b.c+a.b.d可得到三组结果（主要因素、一阶相互作用和二阶相互作用）。最后一组术语的MCI F统计检验这些术语，就好像它们是最后一组术语一样。如果a.b.c不重要，则a.c和b.c的MCI测试对测试这些术语有效，因为它们与a.b.c相差无几。同样，如果a.b.d不重要，则a.d和b.d的MCI测试对测试这些术语有效，如果三者都不相互作用，则a.b的MCI测试有效。离子是重要的。主要影响的测试条件是所有不涉及该主要影响的条款。同时使用两个F统计数据，可以大大减少彻底测试所有术语所需的运行次数。边际性的定义不仅基于基本因素，而且还基于一个术语是否从另一个术语中移除了自由度。例如，在模型MU+区域+位置中，区域不能以位置为条件，因为位置嵌套在区域中。这是隐藏的边缘。
混合模型中wald f统计量的正式测试依赖于分母自由度的计算。Asreml使用Kenward和Roger（1997）的方法，使用数值导数计算分母自由度。但是，如果存在许多方差参数，这是很昂贵的，因为每个参数都需要额外的半迭代（步骤1-5）。

3结束语

ASREML的改进是一个持续的过程；最近的代码更改已经将非解析模型的执行时间缩短了70%以上。由于计算量与CA中的行数乘以其平均稀疏长度直接相关，因此希望通过替代模型公式、更快的计算方法和更优的方程排序，对大型复杂稀疏模型进行进一步的改进。
20多年来，ASREML中增加了新的选项，因此很难进一步适应处理大型密集模型。因此，作者开始编写一个新的方案（echidna，www.echidnams.org）利用多年积累的知识。人们希望这将扩展功能和效率相对于目前的ASREML在一个更好的组织方案，以满足基因组时代的挑战。由于vsni拥有大部分asreml知识产权，作者希望在echidna的写作中与vsni合作，使它可以免费提供给非商业研究人员和学生。Echidna的发展可供VSNI自行决定在ASREML中实施。