「線形混合モデルにおける分散成分の推定」への注釈

はじめに

このページは、私が2007年に公開したweb記事「線形混合モデルにおける分散成分の推定」に関する注釈です。このオリジナルの記事は、当時、分散成分の推定に関する話題が日本語web上にほとんど存在しなかったことを踏まえ、主に研究室の学生に向けて、REMLの導入として書いたものです。また、オリジナルの記事では量的遺伝学、とくに家畜の育種への応用を前提としていますが、これは私の専門を反映したものです。

このページは、その当時の記事に注釈や引用を加えたものです。オリジナルの記事の本文を引用し、それに対してコメントを付加します。2019年現在、多くの情報がweb上に公開され、日本語で読める記事も公開当時と比べて増えました。そこで、オリジナル記事を補足し、また、新たに利用できるようになった情報源を紹介します。なお、このページでもREMLの導入が主眼ですので、厳密な記述は専門書と関連論文にゆずります。

はじめに、分散成分の推定に関して、web上で読める日本語の書籍と論文について紹介します。

もしREMLに関する初等的なテキストを求めているなら、以下の論文が参考になるかもしれません。

量的遺伝学者と統計学者によるREMLの総説論文のいくつかを以下に示します。惜しむらくは、多くがオープンアクセスではないことです。

分散成分の推定と歴史

REMLの導出と解説については，別のページに詳しくまとめたので，そちらを参照してください。 また，1990年以前の分散成分の推定に関する詳細な歴史は，SearleらのVariance Components第2章に詳しく述べられています。

このREMLの導出ページは完成していません。 ここで紹介したSearle, Casella, McCulloch の Variance Components は、線形混合モデルにおける分散成分の推定に関する網羅的な書物です。

分散分析法とその拡張

統計学の教科書には，必ず分散の計算法が登場します。 しかし，その方法（なにやら2乗したものを足していって，最後にn-1で割るもの）は，様々な要因が関与するような測定値には利用できません。 そこで，Sir Ronald A. Fisherが開発した分散分析法（ANOVA；Analysis of Variance）の登場です。 分散分析法を用いると，測定値に影響する要因を分割し，各要因に起因する分散を計算することができます。 しかし，この方法には欠点があり，各要因グループ内のデータ数がそろっている場合（釣り合い型データ；balanced data）しか分析が出来ないようになっています。 つまり，十分に計画立てて実験を行うことを前提とした分析法なのです。 上でも述べたように，実際の家畜の測定値は，条件をそろえることが出来ず，グループ内のデータ数もはらばらです（不釣り合い型データ；unbalanced data）。 このような事情から，分散分析を不釣り合い型データに応用する試みが数多く行われました。

上記の「なにやら2乗したものを・・・」による値は、不偏分散として知られています。 これは、各観察値が平均からどのくらい離れているかを表したものです。 そして、どの観察値も期待値が同じである（すべての観察値が同一の群に属する）という暗黙の仮定があります。 さて、ここで疑問が生じます。 もし観察値が2つの群から得られていて、各群の平均が異なるならば、分散はどのように計算すべきでしょうか。 結論から言えば、いくつかの分散統計量が定義できるので、それらのうち、もっとも目的にかなった値を採用する必要があるということです。

1つの可能性は、分散の計算式に平均値の違いを反映させることです。 つまり、A群の観察値にはA群の平均を、B群の観察値にはB群の平均を当てはめ、全部の観察値を使って分散を計算します。 これは、A群とB群の平均を同じに揃えた上で得た分散であり、群内分散と呼ばれます。 もう1つの可能性は、各群の平均の違いを無視し、観察値全体で1つの平均値だけを想定し、不偏分散を得ることです。 このようにして得た値は、全分散と呼ばれます。 もし、A群とB群の間に平均の違いがあれば、全分散は群内分散より大きくなります（平均の違いが分散に反映されます）。 言い換えれば、A群とB群の間に平均の違いがなければ、全分散と群内分散は同じになるはずです。 そしてこれが、分散分析の原理につながります。

もし、全分散 ＞ 群内分散が成り立つなら、この差は「群間の平均の違い」により引き起こされています。 そしてこの差を群間分散とすれば、ある条件下で、以下の式が成り立ちます。 $$ 全分散 = 群間分散 + 群内分散 $$ 群間分散は、この例を使っておおざっぱに言えば、平均Aと平均Bを新たな観察値とみなして得た分散に相当します。 上で述べたように、平均Aと平均Bに差がなければ、群間分散はゼロになります。

さて、分散は平方和の関数ですので、上記の線形の関係は一般には成り立ちません。 ただし、釣り合い型のデータに対しては、この関係が成立します。 釣り合い型のデータとは、この例ではA群とB群内の観察値数が同じであるようなデータです。 実験計画法に沿えば、最終的には釣り合い型データが得られます。 そして、各分散は上記のようにきれいに分割され、ここから各処置間の平均値に関する検定ができます。 これが分散分析法です。

分散分析法は、一般に群間の平均値に関する検定（F検定）を行う方法として知られています。 しかし、その過程で、群間分散と群内分散を計算することになります（平均平方和がそれに相当します）。 そして、それらの分散は、線形混合モデルを当てはめた場合に仮定される分散成分の線形関数であると期待できます。 もし平均値の差に興味がなく、単に分散成分のみに興味があるならば、途中経過の値から分散成分の推定値を得ることができるわけです。

データが不釣り合いの場合、上記の等式は成り立たず、各要因による分散を一意に得ることはできません。 つまり、群間分散と群内分散は、ad hocな方法でしか得られないのです。 次節にあるHendersonの方法は、この分散成分をできる限り合理的に計算しようとするものでした。 残念ながら、不釣合い型のデータに関して、分散分析を（検定ではなく）分散成分の推定に利用する方法は一般的ではなく、そのため、これを説明する教科書は多くありません。 この事情が、この方法を謎めいたものにしていた一つの理由です。

Hendersonが1953年に発表した論文には，不釣り合い型データに分散分析様の方法を応用し，分散成分を推定するための3つの方法が載っています。 彼の提案した方法は，非常に柔軟性が高かったため，画期的であるとして多方面から激賞されました。 しかしながら，分散分析とその拡張法（Hendersonの方法を含む）は，前提として固定効果モデルを仮定し，変量効果を固定効果と見なして分析する方法でした。 この方法は，データによっては強引で無理のあるものでした。 しかしながら，高速なコンピュータを簡単に利用できなかった時代，たとえ稚拙な方法であっても（分散の推定値が負の値になることも珍しくありませんでした），多くのデータに関して応用したとしても計算時間の短いHendersonの方法が利用され続けていました。 家畜育種学の論文を調べてみると，1980年代末まで頻繁に利用されていたことが分かります。 その当時は，分散成分の推定には固定効果モデルを，BLUPの計算には混合モデルを利用していたのです。 現在でも，SASのProc GLMでは，HendersonのMethod III（方法3）による分散の推定値を計算することが出来ます。

不釣合い型データは、フィールド研究など実験計画法が利用できない場面、あるいは事故等で欠損値が生じた場合に現れます。 そして不釣合い型データしか得られない研究課題でこそ、分散成分の推定が問題になりがちです。 Hendersonの方法が、その後長く使われたのも、このような背景があってのことです。

Hendersonの方法（分散分析法を含む）は、二次形式に基づく方法と呼ばれます。 分散の計算に必要な2乗の項（二次式）の総和は平方和と呼ばれており、より一般にはベクトル $\mathbf{y}$ と対称行列 $\mathbf{Q}$ の二次形式 $\mathbf{y}'\mathbf{Qy}$ として表現できることに由来します。 ここで観察ベクトル$\mathbf{y}$を確率変数とみなせば、二次形式の期待値$E(\mathbf{y}'\mathbf{Qy})$は、分散成分の線形関数で表現できます。 いちどこの関係を導いておけば、この等式の$\mathbf{y}$に実際の観察値を代入し、分散推定値を計算することができます。 そして、この推定量は、その期待値が真のパラメータ（分散）に等しいという性質を持っており、不偏推定量です。

特に引用されたHendersonの方法3は、複数の効果を含む線形モデルに基づきます。 まず、全部の効果を含むモデルを当てはめて、ある種の平方和を計算し、次に目的の効果を取り除いた縮小モデルで同種の平方和を計算します。 モデルの縮小に伴う平方和の減少量に注目していることから、この平方和の差をreduction（リダクション）と呼びます (Searle, 1968)。 そして、このリダクションの期待値が、目的の効果による分散成分の関数であることが導けます。

統計ソフトウエアのRでも、Hendersonの平方和を得ることができます。 興味のある方はType III Sum of Squares（タイプ3の平方和）で検索してみてください。 これは、もともとSASで不釣り合い型データから計算された一意でない平方和を識別するために使われていた用語です。 Type IIIは、Hendersonの方法3に相当する、リダクションを表す平方和です。

さて、この方法のメリットはいくつかあります。 上で紹介されている通り、1990年代以前の限られた計算資源でも計算が可能であったことのほか、あらゆる固定効果モデルにも応用できたこと、観察値の分布に対する仮定が必要ないこと（ただし分散推定値のSEの導出には正規分布を前提とする必要がある）、いくつかのソフトウエアパッケージが入手できたことがあります（SASとLSML76/LSMLMWが代表です）。 一方、メリットを上回るデメリットや疑問もあります。

  1. 推定量がパラメータ空間内におさまる保証のないこと。
  2. この推定量は不偏だが、不偏だけを備えた統計量がどれほど魅力的であるかということ。
  3. 不釣合い型データに対して、無限数の分散分割方法が存在するため、どれが合理的かを決定づける方法がないこと。
  4. 観察値が選抜の結果として得られている場合、分散の推定値に深刻な偏りが生じ得ること。
  5. 固定効果モデルを前提にしているのに、分散成分の推定の際にはその効果をランダム効果として扱っていることに一貫性がないこと。

上記の4.は、家畜育種への応用においては致命的な欠点です。 家畜育種では、すぐれた個体を選抜して（あるいは劣った個体を淘汰して）後代を得て、それを繰り返すことで形質の改善を図ります。 もし、ある個体が繰り返し観察されるなら、最初の観察値が得られた段階で能力の低いものは淘汰されてしまい、優れたものだけが次の観察値を得ることになります。 したがって、このようにして得られた観察値は、選抜（淘汰）の影響下にあります。 分散分析様の二次形式に基づく推定値は、選抜による分散成分の偏りを排除することができません（尤度に基づく方法は、この偏りを大幅に軽減できます）。

また、最後の点に関して、この推定に用いる方程式は正規方程式であって、混合モデル方程式ではありません。 つまり、Hendersonの方法を含む分散分析法ベースの手法は、現在私達がよく知っている混合モデル（ランダムモデル）を用いたものではありません。 たとえば、ランダム効果の水準間に共分散構造がある場合、それを考慮することができません。

正当に混合モデルを当てはめ、そのBLUPを含む二次形式の期待値と取ると、その推定量はより合理的なものになるはずです。 後述するように、このようにして得た推定量はMINQUEと呼ばれ、REML推定量の近似になっています。 そして、この近似を繰り返せば、最終的にはREML推定量が得られます（いわゆるEMアルゴリズムです）。 この二次形式の期待値をとる方法は、1970年代から1980年代にかけて、量的遺伝学者の間でさかんに研究されました。 そしてこれがREML推定量を与えるものだとわかると、研究の興味は、REML推定値の効率的な計算方法へと移っていったのです。

尤度に基づく推定

一般の線形混合モデルに関する分散成分の推定法は，1967年のHartley and J. N. K. Raoが最尤法（ML；Maximum Likelihood）を一般化するまで待たねばなりませんでした。 彼らは，混合モデルのことを分散成分モデル（variance components model）と呼んでいました。 彼らの方法は，分散分析に基づく手法よりも格段に優れていましたが，以下の点から，あまり利用されませんでした。

  1. 計算量がとても多かったこと
  2. 推定値に深刻な偏りが生じる可能性があったこと

1971年にPatterson and Thompsonは，最尤法を修正し，偏りのない推定値を求める方法を発表しました。 これが現在では制限付き最尤法，通称REML（Restricted Maximum Likelihood）と呼ばれている方法です。 ヨーロッパを中心に，残差最尤法（Residual Maximum Likelihood）と呼ばれることもあります。 この方法は，これまでの方法の中で最も好ましい推定値を生み出すことが分かりました。 しかし，家畜育種などの分野のように，データ数が非常に多い場合には，計算時間がかかりすぎるため，後年まで実用にはなりませんでした。 ちなみに，REMLの発表者であるRobin Thompsonは統計遺伝学者であり，家畜育種学の応用論文にも頻繁に名を連ねています。

現在，家畜育種における分散成分の推定法として，REMLが標準的な手法になっています。 これは，コンピュータの高速化・主記憶領域の大型化のみならず，計算アルゴリズムの改善による部分が大きいと言えます。 計算手順は，後の章にて解説します。

最尤法（ML）は、現在ではきわめて広範囲に用いられている統計的方法のひとつです。 線形混合モデルの一般的な仮定は、観測ベクトル$\mathbf{y}$に対して多変量正規分布を当てはめ、$var(\mathbf{y})$は分散成分の線形関数とすることです。 この仮定はMLでもREMLでも同じです。 MLもREMLも、その推定量には各種の漸近的性質が保証されていて、観察ベクトルに対して推定値が一意に決まることから、分散分析様の方法よりも好ましいとみなされています。 ところが、線形混合モデルの分散成分の推定に関して、MLよりもREMLの方が広く用いられています。 これにはいくつかの理由があります。

たとえば、上の記述の 2. はデータ構造によっては深刻な問題になります。 端的な例として、たとえば1群だけを想定した固定効果モデル（$y_i=\mu+e_i$）の残差分散は、REMLであれば不偏分散（分母が$N-1$）が得られますが、MLでは分母が$N$になります。 これは、より多くの固定効果を含めるとき、REMLは残差分散における自由度の減少を説明できるのに対し、MLでは自由度が補正されず値が過小に推定される結果になります。 これは、最尤法が固定効果と分散成分を同時に推定するのに対し、REMLでは$\mathbf{y}$を線形変換して固定効果を消去し、分散成分だけを推定するためです（誤差対比に対する最尤法）。 この原理に関して、プロファイル尤度で説明するほうが通りが良いかもしれません。 より詳しい説明を求めるなら、Harville (1977) かSearle et al. (1992)を参照してみてください。

歴史的には、上で言及した通り、REML相当の推定値を経験的に得る方法が発見され、それが合理的であると考えられた、という事実も重要です。 これは、前述の二次形式の期待値に基づく方法を拡張し、固定効果モデルではなく混合効果モデルを仮定し、正規方程式のかわりに混合モデル方程式を利用することで導けます。 最初に初期の分散成分を与えて混合モデル方程式を解き、得たBLUPとBLUPを用いて分散成分を推定するのです。 分散成分は、固定効果のBLUPとランダム効果のBLUPを含む二次形式で表現されます（ただし分散分析由来のものよりも複雑になります）。 この方法は、一般にMINQUE（Minimum Norm Quadratic Unbiased Estimation）と呼ばれ、多変量正規分布が仮定できる状況ではMIVQUE（Minimum Variance Quadratic Unbiased Estimation）といいます (Rao, 1971a,b; LaMotte 1973)。 さて、その更新された分散成分を使って方程式を作り直し、それを解いて分散成分を更新し・・・を繰り返すと、やがて分散成分はある値に収束します。 こうして得られた推定値は反復MIVQUE推定値と呼ばれ、もしその値が妥当であれば、これはREML推定値と一致します (Harville, 1977)。 そして、この方法そのものが、EMアルゴリズムによるREML推定に相当しています (Dempster 1977)。

分散成分の推定にREMLが広く受け入れられたのは、以上の複数の要因によるものです。

ベイズ統計学に基づく方法

1950年代に始まったベイズ統計学の隆盛により，家畜育種にベイズ統計学を導入する試みがなされました。 1970年代には，David A. Harvilleによって，分散のMLならびにREML推定量に関するベイズ流解釈が発表されています。 しかしながら，実際のデータから分散成分を推定するためにベイズ統計学を利用する実務家は，多くありませんでした。 それは，ベイズ統計学において分散成分の推定に必要な「積分」の計算が不可能であったこと，旧来の方法に勝るメリットが明確でなかったことが挙げられます。

ここで示したHarvilleの論文は、つまり、ベイズ統計学におけるある条件下での分散成分の推定結果は、REMLと同じであることを証明しています。厳密には、REML統計量が、分散成分に対して無情報事前分布を仮定することで得られた周辺事後分布（尤度）のモードと一致することを示しています。以下に日本語の解説があります。

nshi.jp

オリジナルの原稿にある「積分」は、同時事後分布の周辺化を意味しています。 つまり、分散成分に任意の共役事前分布を仮定した場合、周辺事後分布を導くには積分が必須ですが、その結果は解析的に得ることができません。 この事実は、分散成分の事後周辺分布は数値的にしか得られず、マルコフ連鎖モンテカルロ法（MCMC）が広く応用される前、その計算は極めて複雑であったことを意味します (Gianola nad Foulley, 1990)。

ところが，1984にGeman and GemanがGibbs Samplingと呼ばれる手法を定式化したことにより，直接積分を実行すること無く，最終的な結果が数値的に得られるようになりました。 1990年代に入り，Gibbs Samplingは爆発的に普及し始めます。 家畜育種の文脈でGibbs Samplingが初めて応用されたのは，1993年のことです。 その後，徐々に複雑なモデルへと応用が進み，現在ではREMLと並んで，よく利用される方法となっています。REMLよりもGibbs Samplingが利用されるのは，以下のような場面に多く見られます。 主に計算上の理由であることが多いのです。

  1. 方程式のサイズが非常に大きい場合（REMLでは主記憶容量と計算時間の面から，現実的ではないことが多い）
  2. 数学モデルが非常に複雑な場合（REMLは収束しにくい）

Gibbs samplingが家畜育種で分散成分の推定に本格的に応用されたのは、Wang et al. の1993年の研究が最初です。その後、多形質モデル、閾値モデル、ランダム回帰モデル、母性効果モデルなど、量的遺伝学において重要なモデルに速やかに適用されてゆきました。

これが書かれた2007年当時と2019年現在とで大きく異なるのは、ベイズ統計学がさまざなま自然科学分野で広く受け入れられる方法になったことと、Rをはじめとするソフトウエアが容易に利用できるようになったことです。テキストが充実し、ツールも利用できるようになった今、Gibbs samplingとREMLのどちらを選ぶかは、分析者の考え方によります。

さて、上の記述は、実務家の視点からのもので、これに関して2019年現在でもさほど変わるところはありません。量的遺伝学（特に家畜育種）の応用において、計算コストと推定値の安定性を重視して手法を選びたいのであれば、Misztal (2008)が参考になるかもしれません。この論文には、REMLの計算手法についての助言もあります。

プログラマの視点では、Gibbs samplingはREMLに比べて、プログラミングが容易ということも特筆すべき点です。 これは、Gibbs samplingのもっとも典型的な実装が、方程式の反復解法であるガウス・ザイデル(Gauss-Seidel)法とほぼ同じことによります（Sorensen and Gianola, 2002）。 線形方程式の反復解法のプログラムがすでに存在すれば、あとは適切な乱数生成ルーチンを用意するだけで、大きな困難なく単純なGibbs samplerを実装できます。

量的遺伝学でのGibbs samplingの利用を想定したパッケージは、Rを中心に数多く公開されています。 そのRパッケージの中でも、BGLR (Bayesian Generalized Linear Regression)は最近開発されてよく利用されています。 また、大規模なフィールドデータをサポートする単独パッケージとして、MisztalらのGibbsf90、Van TasselとVan VleckのMTGSAM、MeyerのRRGIBBS、Madsenらのrjmcがあります。

REML推定値を計算するためのアルゴリズム

REML推定量は非線形方程式の解であるため，通常はNewton-Rhapson法またはそれに準じた方法（準Newton法）を利用します。 しかし，観測値数が非常に大きいとき，Newton法に必要なヘッセ行列（またはヘシアン；Hessian）を求めるための計算量が増大し，コンピュータの能力の限界を超える可能性が高まります。 家畜育種では，一般に多数のデータを利用するため，アルゴリズムの改善によって計算量を減らす工夫がなされてきました。 以下，実用化された順に紹介してゆきます。

現在、REMLを用いる多くのケースでは、計算量とアルゴリズムは問題にならないかもしれません。 データがコンパクトで、モデルもシンプルならば、最近のコンピュータで十分対処できるでしょう。 この場合、後述するDFまたは準ニュートン法を使えば、実用時間内でREML推定値が得られます。 家畜育種では、数多くの固定効果のほか、数万から数十万の個体をランダム効果とみなし、その共分散構造（血縁関係）をモデルに含めます。 したがって、REML推定値のための計算量は膨大になり、これはコンピュータが非力だった時代には大きな問題でした。 そのため、多くの方法が提案され、大規模なデータに対応する計算プログラムも多く書かれました。

オリジナルのテキストでは、実用化された順にと書いてありますが、正しくは量的遺伝学（特に家畜育種）の文脈において広く用いられるようになった順に紹介しています。ここで紹介しなかった方法にも、いくつもの提案があり、計算パッケージに実装されたものもあります（SAS Proc MixedのW transformationなど）。

なお、多変量の場合にcanonical transformationを当てはめるケース、Factor Analysisとの組み合わせ、罰則付き最尤法など、多様なモデルに拡張されていますが、それらをここで扱うことはしません。このあたりはKarin Meyerが幅広く研究しているので、彼女のWebサイトから業績をたどるのが確実だと思います。

導関数を使用しない方法（Derivative Free REML）

単にDF REML，DFアルゴリズム，デリバティブフリー法などとも呼ばれます。 動作の原理は簡単で，パラメータを試行錯誤的に変化させながら尤度を計算するだけです。 尤度が大きくなる方向へパラメータを動かしてゆき，尤度の値がほとんど変化しなくなった時点で検索を終了します。 この方法は，1986年ころにGraserとSmithによって提案され，後にKarin MeyerがDFREMLパッケージとして公開したため，広く使用されるようになりました。 1990年に入る頃には，Boldman and Van VleckがDFREMLを大幅に改良し，疎行列パッケージを組み込んで多変量分析にも対応させたMTDFREMLプログラムを開発しました。 DFREMLおよびMTDFREMLプログラムは，爆発的に普及し，家畜育種学で分散成分を推定する際には，必ず使用されていた時期がありました。

最尤法は、尤度関数を最大にする値を得る方法なので、いろいろな最適化手法を使うことができます。ここでDFと呼ばれているのは、Nelder-Mead法やPowell法などのヒューリスティックを用いて、尤度関数の最大値とそれに対応する分散成分値を得るものです。なお、Nelder-Mead法は、そもそも尤度関数の最適化問題のために開発されたので、まさにDF REMLにうってつけです。

Karin MeyerのDF REMLは発展的に消滅し、WOMBAT（ウォンバット）として生まれ変わりました。このプログラムでは、もはやDFアルゴリズムは使われていません。彼女のWebサイトからプログラムをダウンロードできます（リンクは後ほど紹介します）。BoldmanとVan VleckのMTDFREMLは、すでに十分古いソフトウエアですが、作者のひとりによると、いまでも年に何度か問い合わせがあるそうです。

この方法の欠点は，言うまでもなく，試行錯誤型の検索が上手くゆくとは限らないことです。 モデルが複雑になると，尤度の最大点（global maximum）ではなく極大点（local maximum）に落ち着いてしまうのです。 そのため，初期値を変えて何度も繰り返し実行する必要がありました。 MTDFREMLを使用している1990年代の論文を読むと，コールドスタート（cold start）なる用語が登場しますが，これは初期値を変えて再度実行するという意味です。 計算量は少ないですが，収束まで時間がかかり，収束値が妥当かどうか判断しにくいので，現在では使われなくなってきています。

これらヒューリスティックの欠点は上記のとおりで、家畜育種で利用する人はほぼ皆無になりました。その理由は、コンピュータの能力が向上した上、導関数を用いるソフトウエアも利用可能になり、DF REMLを使わないほうが計算時間を短縮できるようになったためです。また、当てはめるモデルが年々複雑になり、多くのパラメータを含むモデルでは、DF REMLが収束に困難をきたすことも理由です。なお、Rのlme4はデフォルトでPowellのBOBYQA法を使っています。

1次導関数を使用する方法（EM REML）

EMとは，EMアルゴリズム（EM algorithm）を意味しています。 EMアルゴリズムは，1977年にDempsterらが発表した方法で，DF REMLよりも早い段階から計算式は知られていました。 しかしながら，EM REMLの計算式自体は極めて簡単であるにもかかわらず，その中には巨大な逆行列が含まれていたために計算がほぼ不可能であり，長らく実用化されてきませんでした。

ここで計算が不可能としているのは、1990年代初頭までの実務データへの応用です。Dempsterらの有名な1977年の論文では、REML推定量がEMアルゴリズムにより算出できることが明記されています。そして、その論文に付随した長いdiscussionではRobin ThompsonとShayle Searleが分散成分の推定に関するコメントを寄せており、EMアルゴリズムの収束の遅さと計算量の多さへの懸念が示されています。現在では後述するような工夫により、かなり大きなデータに対しても利用できます。また、導関数を使う場合のREMLの計算方法については、Meyer and Smith (1996)に詳しい記述があります。

1990年代に入り，Perez-Enciso，Misztal，Elzoはフリーな疎行列パッケージFSPAKを開発し，そこに逆行列を計算するサブルーチンを組み込みました。 これは，逆行列の全要素を計算するのではなく，元の行列で値が入っている部分の要素だけを計算するアルゴリズム（Takahashi algorithmと呼ばれます）を利用しています。 EM REMLでは，この形式の逆行列を利用しても推定値を得ることが可能でした。 DF REMLに取って代わる形で，現在でも使用されています。 FSPAKの発表以降，比較的簡単にプログラミングできるため，分析者が独自にソフトウエアを開発することも多いです。 また，FSPAKをFortran95でモジュール化したライブラリ（SPARSEM）も利用できます。 利用可能なパッケージには，MisztalらのREMLF90，MadsenらのDMUなどがあります。

この項目は私の研究に関連するため、すこし詳しく説明させてください。

尤度の1次導関数は、後述するquasi-Newton法（2次導関数に基づく方法）でも、勾配ベクトルの計算に必要とされます。Misztal and Perez-Enciso (1993)は、この勾配の算出には混合モデル方程式の係数行列の部分逆行列（疎逆行列）のみが必要で、これがTakahashiアルゴリズムにより現実的なコストで計算できることを示しました。この部分逆行列の計算は、係数行列に密構造が含まれるときに計算コストが跳ね上がります。このようなケースは、量的遺伝学であれば多形質モデル（分散共分散行列）あるいはゲノミックモデル（ゲノム関係行列）で頻繁に出現します。Takahashiアルゴリズムを発展させ、密構造をうまく扱えるようにした方法がいくつか提案されています。たとえば、マルチフロンタル法をはじめとしたsupernodal法を逆行列の計算に応用した、Inverse Multifrontal Methodがあります (Campbell, 1995)。この手法が量的遺伝学における諸問題の解決に有効であることを実証することが、私の研究課題の一つでした。これに取り組んでいた当時、そのアルゴリズム自体は公開されていたものの、既存のパッケージが存在しなかったため、自分で開発する必要がありました (Masuda et al., 2014)。この成果は、現在ではBLUPF90パッケージに組み込まれ、計算速度の改善に貢献しています (Masuda et al., 2015)。

なお、モンテカルロ法を使えば、この部分逆行列を数値的に得る（近似する）ことができます（García‐Cortés et al., 1992）。これを用いた方法はモンテカルロEMと呼ばれ、後の文献にも散見されます。

EM REMLの最大の欠点は，収束が極めて遅いことです。 現在のパラメータが収束値に近づくにつれて値の変動が小さくなるという性質があるため，推定値がほとんど動かなくなっても収束していない可能性があります。 また，複雑なモデルでは，小刻みに値が変動するので，相当回数反復させても値が固定しません。

EMアルゴリズムは、分散成分の初期値を与え、それを繰り返し更新し、その値が十分に動かなくなったところで反復を打ち切ります。 初期値がパラメータ空間内に収まっていれば、理論上、収束値もパラメータ空間内に収まります。 しかし、収束が遅く、さらには必ずしも尤度関数の最大値に収束するとは限りません。 EM REMLは、後述するAI REMLの初期値を得る目的のほか、データに対して複雑なモデルを当てはめたい場合にも利用されます。EM REMLを加速するための手法（Aitkenの方法）も提案されており、REMLF90がそれを実装しています。 また、PX-EM法はEM REMLの反復回数を減らす方法として有効です (Foulley and Van Dyk, 2000)。

2次導関数を使用する方法（AI REML）

Newton-Rhapson法はヘッセ行列の計算が困難でした。 ここで用いるヘッセ行列を観測情報行列（observed information matrix）としておきます。 一方，準Newton法として知られるFisher scoring法は，ヘッセ行列としてFisher情報行列を利用します。 しかしながら，線形混合モデルでは，Fisher情報行列も計算量が極めて多いのです。

Newton-Raphson法（準Newton法）は、典型的な反復解法です。分散成分の推定で問題になる計算量の多さは、この情報行列に含まれる行列の対角和に由来します。その行列を計算してから対角和を得る他ないため、計算を簡単にすることができないのです。

ところが，1994年にJensenらのグループMadsen, Jensen, ThompsonのグループとJohnson & Thompsonの両グループは，観測情報行列とFisher情報行列期待情報行列を足して2で割った行列は，比較的簡単に計算できることを発見しました。 この行列を平均情報行列（average information matrix）とし，これをヘッセ行列とする準Newton法がAI REMLです。 そのグループには，REMLを発表したThompsonも含まれており，翌年，彼は実際のデータを使用してAI REMLの動作を実証しています。 AI REMLは，従来の方法と比較して多少の計算量の増加はあるものの，劇的に少ない反復回数で収束することが知られています。 また，収束時の情報行列は，そのままパラメータの漸近分散として利用できるため，分散成分の標準誤差を容易に得ることができます。

平均情報行列は、正しくは、観測情報行列と期待情報行列（Fisher情報行列の符号を反転したもの；あるいは観測情報行列の期待値そのもの）の平均値です。 こうすることで、計算が困難な対角和を消去し、計算可能な二次形式の項だけを残すことができます。 また、オリジナルの記事の中で、平均情報行列を導出したグループへの言及が正しくありませんでした。 これは、1994年のWCGALP (World Congress on Genetics Applied to Livestock Production, Guelph, Canada)において、Madsen, Jensen, Thompsonによるグループと、JohnsonとThompsonのグループの両者によって発表されました。前者のグループは、一般的なモデルにおける導出方法を示し、後者はよりシンプルなモデルへの応用について述べています。両者とも、後に論文としてこれらの成果を公式に発表しています。また、同時期にGilmour, Thompson, CullisもAI REMLに関する論文を出しており、これはASRemlパッケージの基礎論文になっています。これらの発展について、明らかにRobin Thompsonが重要な役回りをしたことが伺えます。論文の詳細は、記事末尾の文献リストをにまとめてあります。

なお、日本語で「制限最尤法」「AIアルゴリズム」などの語でインターネット検索すると、新潟大学（当時）にいらした祝前先生のグループの研究がヒットします（例えば分散成分のAIアルゴリズムによる制限最尤推定のための一演算手法ならびにその特性）。これは、あるモデル下では、平均情報行列がある種の双一次形式と二次形式で表現され、これが計算上有利であることを示した研究です。

欠点としては，他のNewton法と同様に発散の危険が高いこと，EM REMLよりは計算量が多いので，コンピュータの能力に関する不安が存在していることです。 AI REML推定を行うパッケージは非常に多く，MisztalらのAIREMLF90，MeyerのWOMBAT，MadsenらのDMU，商用ではGilmourらのASREMLにおいて利用できます。

これらのパッケージへのリンクを紹介します。

なお、ASRemlはGilmourの手を離れ、VSN International社が完全にその知的財産権を取得しました。この商用化の過程で新規開発が滞るようになったようです。Gilmourはその後、ASRemlの機能を一から実装し直し、学生と研究者が自由に使えるソフトウエア Echidna を公開しています。

なお，AI REML以外にも，いくつかの準ニュートン法による分散成分の推定法が提案されています。それらの手法は，たとえばSASのMIXEDプロシジャやGroeneveldらのVCEにおいて応用されています。

SASのMixedプロシジャには、W-transformationと呼ばれる方法をREML向けに改良した手法を用いています (Hemmerle and Hartley, 1973; Corbeil and Searle, 1976)。 また、VCEソフトウエアは、AI REMLとは別の準ニュートン法を実装しています (Neumaier and Groeneveld, 1998)。 2019年現在、VCEの配布サイトに接続ができず、これを入手することができません。

参考文献

  • Boldman, KG, Van Vleck, LD. 1991. Derivative-free restricted maximum likelihood estimation in animal models with a sparse matrix solver. J. Dairy Sci. 74:4337-4343.
  • Campbell, YE. 1995. Multifrontal algorithms for sparse inverse subsets and incomplete LU factorization. PhD thesis, University of Florida, Gainesville, FL.
  • Corbeil, RR, Searle, SR. 1976. Restricted maximum likelihood (REML) estimation of variance components in the mixed model. Technometrics 18:31-38.
  • Dempster, AP, Laird, NM, Rubin, DB. 1977. Maximum likelihood from incomplete data via the EM algorithm. J. Royal Stat. Soc. Series B (Methodological). 39:1-22.
  • Foulley, JL, Van Dyk, DA. 2000. The PX-EM algorithm for fast stable fitting of Henderson's mixed model. Genet. Sel. Evol. 32:143-163.
  • García‐Cortés, LA, Moreno, C., Varona, L., Altarriba, J. 1992. Variance component estimation by resampling. J. Anim. Breed. Genet. 109:1-6.
  • Gianola, D, Foulley, JL. 1990. Variance estimation from integrated likelihood (VEIL). Genet. Sel. Evol. 22:403-441.
  • Gilmour, AR, Thompson, R, Cullis, BR. 1995. Average information REML: An efficient algorithm for variance parameter estimation in linear mixed models. Biometrics 51, 1440-1450.
  • Graser, HU, Smith, SP, Tier, B. 1987. A derivative-free approach for estimating variance components in animal models by restricted maximum likelihood. J. Anim. Sci. 64:1362-1370.
  • Hemmerle, WJ, Hartley, HO. 1973. Computing maximum likelihood estimates for the mixed A. O. V. model using the W Transformation. Technometrics 15:819-831.
  • Henderson, CR. 1953. Estimation of variance and covariance component. Biometrics 9:226-252.
  • Hartley, HO, Rao, JN. 1967. Maximum-likelihood estimation for the mixed analysis of variance model. Biometrika 54:93-108.
  • Harville, DA. 1974. Bayesian inference for variance components using only error contrasts. Biometrika 61:383-385.
  • Harville, DA. 1977 Maximum likelihood approaches to variance component estimation and to related problems. J. Amer. Stat. Assoc., 72:320-338.
  • Jensen, J, Mäntysaari, EA, Madsen, P, Thompson, R. 1996-1997. Residual maximum likelihood estimation of (co)variance components in multivariate mixed linear models using average information. J. Indian Soc. Agr. Stat. 49:215-236.
  • Johnson, DL, Thompson, R. 1994. Restricted Maximum Likelihood estimation of variance components for univari-ate animal models using sparse matrix techniques and average information. Proc. 5th World Congr. Genet. Appl. Livest. Prod, University of Guelph, Guelph, Canada. 18, 410–413.
  • Johnson, DL and Thompson, R. 1995. Restricted maximum likelihood estimation of variance components for univariate animal models using sparse matrix techniques and average information. J. Dairy Sci. 78, 449-456.
  • LaMotte, LR. 1973. Quadratic estimation of variance components. Biometrics 29:311-330.
  • Madsen, P, Jensen, J, Thompson, R. 1994 Estimation of (co)variance components by REML in multivariate mixedlinear models using average of observed and expected information. Proc. 5th World Congr. Genet. Appl. Livest. Prod, University of Guelph, Guelph, Canada. Vol 22, 19–22.
  • Masuda, Y, Aguilar, I, Tsuruta, S, Misztal, I, 2015. Acceleration of sparse operations for average-information REML analyses with supernodal methods and sparse-storage refinements. J. Anim. Sci. 93:4670-4674.
  • Masuda, Y, Baba, T, Suzuki, M. 2014. Application of supernodal sparse factorization and inversion to the estimation of (co)variance components by residual maximum likelihood. J. Anim. Breed. Genet. 131:227-236.
  • Meyer, K, Smith, SP. 1996. Restricted Maximum Likelihood estimation for animal models using derivatives of the likelihood. Genet. Sel. Evol., 28:23-49.
  • Misztal I. 2008. Reliable computing in estimation of variance components. J. Anim. Breed. Genet. 125:363-370.
  • Misztal, I, Perez-Enciso, M. 1993. Sparse matrix inversion for restricted maximum likelihood estimation of variance components by expectation-maximization. J. Dairy Sci. 76:1479-1483.
  • Neumaier, A, Groeneveld, E. 1998. Restricted maximum likelihood estimation of covariances in sparse linear models. Genet. Sel. Evol. 30:3-26.
  • Patterson, HD, Thompson, R. 1971. Recovery of inter-block information when block sizes are unequal. Biometrika 58:545-554.
  • Rao, CR. 1971a. Estimation of variance and covariance components — MINQUE theory. J. Multivariate Anal. 1:257-275.
  • Rao, CR. 1971b. Minimum variance quadratic unbiased estimation of variance components. J. Multivariate Anal. 1:445-456.
  • Searle, SR. 1968. Another look at Henderson's methods of estimating variance components. Biometrics 24:749-778.
  • Searle, SR, Casella, G, McCulloch, C. 1992. Variance Components. John Wiley & Sons. NY.
  • Smith, SP, Graser, HU. 1986. Estimating Variance Components in a Class of Mixed Models by Restricted Maximum Likelihood. J. Dairy Sci., 69:1156-1165.
  • Sorensen, D, Gianola, D. 2002. Likelihood, Bayesian, and MCMC methods in quantitative genetics. Springer Science & Business Media. NY.
  • Wang, CS, Rutledge, JJ, Gianola D. 1993. Marginal inferences about variance components in a mixed linear model using Gibbs sampling. Genet. Sel. Evol. 25:41-62.