はじめに

このページは、私が2007年に公開したwebサイト「線形混合モデルにおける分散成分の推定」のミラーです。これは、もともと私が大学院生として帯広畜産大学に在籍していた当時に作成したものです。このWebサイトは、私が帯広を離れてからも長らく公開されていましたが、このたび大学全体のWebサイトのリニューアルに伴って閲覧できなくなったようです。（追記：Slideshareに上がっているスライドも、もともと私が作成したものです。私がアップロードしたものではなく、それを入手したどなたかが公開したものと思われます。）

いくつかの誤りや記述ミスなどもありますが、オリジナルの本文そのままに、ここに改めて公開します。オリジナルの完全なミラーはInternel Archiveで閲覧することができます。

また、この本文の不足を補い、より専門的な内容を付け加えた記事を新たに書き起こしました（線形混合モデルにおける分散成分の推定への注釈）。興味のある方は、そちらもあわせてご覧ください。

分散成分の推定と歴史

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

分散分析法とその拡張

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

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

尤度に基づく推定

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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