遺伝型値

メンデル遺伝形質を例にして

この記事では、量的形質がどのように親から子に受け継がれるかのモデルを考えてみます。 例として、ある植物（2倍体）の花の色を取り上げます。 この形質は，完全なメンデル遺伝形質だと仮定します。 つまり、この花の色は染色体上の1箇所の遺伝子座によって完全にコントロールされていて，環境の影響はありません。 これは量的形質ではありませんが、量的形質の遺伝継承を考える上での簡単なモデルなのです。

遺伝子型値

この遺伝子は、花の色素量を調節します。 その遺伝子が色素を作らないなら、花の色は白になります。 色素が十分に作れるときは、赤い花になります。 色素量がそれらの中間であれば、花の色はピンクになります。 色素の発現量を5段階のスコア（1が白で5が赤）で表せば、花の色は数値で表現できます。

色素量を決める遺伝子は2種類あり、それぞれAとaと表記します。 Aは色素を作る遺伝子、aは色素を作らない遺伝子です。 これらの遺伝子をアリル（または対立遺伝子）と呼びます。 個体の持つアリルの組み合わせ（遺伝型）は、AA、Aa、aaのいずれかになります。 それぞれの遺伝型の発現値（遺伝型値）は、以下のようになるでしょう。

念のため書き下すと、AAの遺伝型値は5、Aaは3、aaは1です。 さて、これらの遺伝型値には明らかな傾向があります。 遺伝型がaa -> Aa -> AAと変化するにつれ、つまりアリルAを1つ多く持つたびに、数値が2ずつ大きくなっています。 グラフで示すと、遺伝型とその値は以下のように直線上に並びます。

このような傾向があるとき、アリルと遺伝型値との関連は相加的（加算的）であるといいます。 そして、この条件下で展開される議論は、相加的モデルに従うといいます。

平均値とアリルの効果

さて、アリルA単体では、遺伝型値をどのくらい増やす効果があるでしょうか。 アリルAは数値を2増やす効果があるように見えますが、そうではありません。 アリルAを1個多く持つのは、アリルaを1個減らすのと引き換えですから、これも考えなければなりません。 これは、遺伝型値を平均からの差で表現するとはっきりします。

たとえば、これら3つの遺伝型が同じ割合で集団に存在するならば、これらの遺伝型値の平均は3です。 上の表は、以下のように表現できます。

この表現は、以下のことを意味します。

• どの個体も、等しく遺伝平均値3を発現する能力を持っている。
• ただし、アリルAを1つ持つ個体は発現値が+1になり、アリルaを1つ持つ個体は-1になる。
• つまり、遺伝型値は、遺伝平均値3とアリル効果の合計値である。

アリルAは遺伝型値に+1の貢献を、アリルaは-1の貢献をする効果があると言えます。 どちらの効果も、符号が違うだけで絶対値は同じです。 これは相加的モデルを想定しているからです。

相加的モデル

ここまで検討したことをまとめます。

この例では、遺伝型値がそのまま表現型の値に反映されます。 式で表せば、以下の関係が成り立っています。

表現型値 = 遺伝型値

相加的モデルが仮定できれば、上の遺伝型値は以下のように分割されます。

表現型値 = 遺伝型平均値 + アリル効果の合計値

このアリル効果の合計値は、平均をゼロとしたプラス・マイナスの数値で表現されます。

相加的モデルでの遺伝継承

さて、上で示した相加的モデルものとで、遺伝型値がどのように親から子に受け継がれるかを考えてみます。 親が取りうるパターンは、AA x AA, aa x aa, AA x Aa, AA x aa, Aa x Aaの5通りです。 それぞれの交配パターンと、生まれる子の組み合わせを以下に示しました。 この図には各個体の遺伝型値も掲載してあります。

この図から、以下の事実が読み取れます。 親がAA x AAまたはaa x aaの場合、子は親と全く同じ遺伝型をとります。 つまり、親と子の遺伝型値は同じです。 両親の遺伝型が異なる場合、子の遺伝子型はさまざまに分離します。 中には両親よりも大きな値を持つもの、あるいは小さな値を持つものが現れます。 そうであっても、子の遺伝型値の平均をとると、両親平均と同じになります。 たとえばAa x Aaの交配では、子はAA:Aa:aa=25%:50%:25%のようにばらつきますが、その遺伝型平均は3になり、両親平均値と同じです。

上で示した関係は、遺伝型値ではなくてアリル効果の合計値でも成り立ちます。 つまり、上の図の遺伝型を、AAで+2に、Aaで0に、aaで-2に置き換えても、子のアリル効果の平均値は両親平均に等しいことがわかります。

相加的モデルの意味するところ

メンデルの法則が厳密に当てはまる状況ならば、上で示した遺伝継承の規則性は、以下のようにまとめることができます。

• 相加的モデルでは、子の遺伝型値（アリル効果の合計値）は、メンデルの法則に従って分離する。
• 子の遺伝型値の平均は、親の遺伝型値の平均値と同じである。

あとから検討しますが、親の世代の遺伝平均値と子の世代の遺伝遺伝平均値は同じである、という事実も間接的に得られます。

この例では、たった1つの遺伝子座だけを考えましたが、これを複数の遺伝子座が関与する場合に拡張することができます。 たとえば、2つの遺伝子座があり、それぞれA/aとB/bのアリルが独立に表現型に関与するならば、上の相加的モデルは以下のようになります。

表現型値 = (遺伝型平均値1 + アリル効果の合計値1) + (遺伝型平均値2 + アリル効果の合計値2)

これは結局、以下のようにまとめることができます。

表現型値 = 遺伝型平均値 + アリル効果の合計値

このように複数の遺伝子座の効果をまとめても、上記の規則性は常に当てはまります。

相加的モデルの仮定は、動植物の改良にとって非常に便利です。 両親のアリル効果合計値を推定できれば、子の遺伝型値を予測することができます。 つまり、大きな数値をもつ親を選んで交配すれば、子も平均して大きな数値を持つことになるということです。 そして、子の遺伝型値がどのくらいばらつくのかも、ある程度は予測できるでしょう。 これは、交配によって遺伝型値の低い個体がどのくらい生じるかを見積もる（リスクを評価する）のに有用です。

一般の場合

さて、ここまでは、特別に作られた例に基づいて議論してきました。 遺伝型値が相加的であり、すべての遺伝型が等しく集団内に存在し、さらにいくつかの理想的な条件もありました。 この議論を一般化するならば、以下の2つの問題を解決する必要があります。

問題1: 平均値とアリル効果

ある集団にはAA個体しか存在しないとしましょう。 この場合、上で示した遺伝型値の平均は5になり、各遺伝型値は以下のようになるでしょう。

• AA: 平均値5 + 0
• Aa: 平均値5 - 2
• aa: 平均値5 - 4

さて、この状況ではAaとaaの個体は存在しないので、それらの遺伝型について議論する意味はありません。 そして、AAの全個体の遺伝型値は、平均値と同じ5です。 アリルAを持っていても、平均値からプラスにもマイナスにもなりません。 つまり、アリルAの効果は0です。

最初の例で示した集団ではアリルAの効果は+1なのに、この極端な集団ではアリルAの効果は0です。 つまり、アリルの効果は、集団内にAA、Aa、aaがどの割合で存在しているかによって異なるということです。 これは、アリル効果が「メリット」であると考えれば、より明確になります。 つまりアリル効果は、そのアリルを持つことで、集団平均よりもどのくらい表現型値が大きく（小さく）なるかを表しています。 みんなが同じアリルを持ってしまえば、そのアリルは平凡なものになり、そのメリットは消滅します。

問題2: 優性の存在

メンデルの法則では、アリルの片方が優性であることが想定されています。 つまり、AaがAAと同じ発現をするということです。 この場合、遺伝型値は以下のようになります。

• AA: 5
• Aa: 5
• aa: 1

この遺伝型値は相加的ではありません。 この状況では、アリルAが増えれば単純に値がどのくらい増えるかを議論することはできません。 しかし、相加的モデルは育種において非常に有用なので、この状況でもそれが利用できる方法を考える必要があります。

次の記事では、この相加的モデルを一般的なケースに拡張します。

メンデルの法則とモデル

メンデルの法則は，遺伝の原理を簡潔にあらわしたものです。 複雑な遺伝現象を単純化して，できるだけシンプルなルールだけで説明しようとしています。 このように，複雑な現象を人間が理解できるよう，その本質を抜き出して単純化したものをモデル（model）といいます。 メンデルの法則は，まさに遺伝現象のモデルなのです。

モデルを当てはめることは，何らかの仮定をおき，いくつかの事実を無視することです。 そのため，現実の事象をおおよそうまく捉えることはできるでしょうが，完璧ではなく，ときに例外が生じます。 そして，その例外を説明するため，モデルはさらに拡張されてゆきます。

メンデルの法則は，遺伝現象を確率で説明するモデルです。 このモデルのうち，もっとも単純なものは，2倍体の生物を仮定して以下のように記述できます。

1. ある形質には遺伝子が関与すると仮定する。
2. その遺伝子には2つのタイプがある。記号で書けば，Aまたはaである。これをアリルという。
3. それぞれのアリルは，機能を持っている。たとえば豆の色ならば，Aは黄色でaは緑色を発現するはたらきがある。
4. 各個体は2個のアリルを持っていると仮定し，その2個の組み合わせ（遺伝型）で遺伝子の発現が決まる。記号で書けば，遺伝型はAA，Aa，aaのいずれかである。
5. 異なるタイプのアリルを持つ個体，つまりAaは，Aまたはaのどちらか一方が発現する。発現するほうのアリルを優性（顕性），覆い隠されるほうのアリルを劣性（潜性）という。
6. 親は，自分の持つ2個のアリルのうちの片方だけを子に伝える。どのアリルを伝えるかはランダムである。もしその生物に父と母がいるなら，子は父と母からそれぞれ1個ずつアリルを受け取る。これを分離の法則という。
7. もしこれ以外の遺伝子が存在するなら，上記の仮定は，それぞれの遺伝子ごとに別個に当てはまる（独立である）。たとえば，Aとaのほかに，別の遺伝子Bとbがあるなら，上記のルールはA/aとB/bにそれぞれ当てはまる。これを独立の法則という。

メンデルの法則で最も重要なのは，分離の法則です。 親から子に遺伝子を伝えるプロセスは，完全にランダムだと仮定しているのです。 つまり，親のもつどちらのアリルが伝わるかは確定できませんが，その確率は50%ずつだと想定しています。 このおかげで，遺伝現象のモデルは，確率理論と統計学を用いてより精密化できるようになるのです。

さて，上記のモデルは，ごく単純な遺伝様式を持つ形質には当てはまりますが，一般の形質には当てはまりません。 つまり，その近似には限界があります。

1. これは常に正しいと仮定する。
2. 2種類以上のアリルが存在することがある。たとえばABO式血液型では，A，B，Oの3種類のアリルが存在する。
3. それらのアリルが機能を持っている必要はない。この場合，遺伝子ではなく遺伝マーカーと呼ぶほうが混乱が少ない。
4. その形質の表現型が遺伝子で完全にコントロールされているなら，この記述は正しい。もし，その表現型が遺伝と環境の両方で決まるなら，遺伝的発現は環境の影響で変化するかもしれない。また，性染色体上の遺伝子は必ずしもこの仮定には当てはまらない。
5. 初等的な遺伝学の教科書では，これを優性の法則と呼ぶが，これは例外が多すぎて法則と呼ぶのはふさわしくない。Aaの発現は形質によって異なり，必ずしもホモ接合体と同じにはならない。
6. その遺伝子が性染色体上にあれば，この仮定は必ずしも成り立たない。
7. 異なる遺伝子が同染色体上に連鎖していれば，この法則は成り立たない。また，遺伝子間に相互作用があるなら，遺伝子の発現も独立であるとはいえない。

このように，多数の例外があるにもかかわらず，メンデルの法則は遺伝の基本原則として教科書には必ず採用されています。 この法則にいくつかの拡張を加えれば，簡潔さを保ちつつ，例外にもいくらか対応できます。 そして，これが実際に細胞内の染色体の挙動と合致し，遺伝物質の継承のモデルとしても優れています。 また，メンデルの法則は，遺伝継承だけでなく，遺伝的変異（表現型の個体差・多様性）を説明するための基礎にもなるのです。

量的遺伝学

量的遺伝学とは，量的形質を対象とする遺伝学です。端的には，ある量的形質について，それがどのように親から子に遺伝継承されるか，ある集団における遺伝的な個体差（遺伝的多様性）はどのくらいあるか，を扱う学問です。

量的形質は，遺伝要因と環境要因の両方が表現型に影響する形質です。加えて，表現型に関与する遺伝子は1個だけでなく，多数が関与すると想定される形質でもあります。多数の要因が表現型に影響するところから，複合形質と呼ばれることもあります。生物の身長（体長）や体重（重量），植物の子実の数，さらには生物の寿命（死亡までの日数）など，かなり多くの形質がこの性質を満たすと言えるでしょう。結果として，量的形質の多くの表現型は，連続的な数値をとります。なお，表現型が数値でなくても，上記の性質を満たすものであれば，それは量的形質です。たとえば疾病に罹患したかどうかを考えると，その表現型は「罹患した」か「罹患しないか」の二者択一であって数値ではありません。しかし、多数の要因が関与すると想定できるならば，それは量的形質なのです。

一方，上記の性質を満たさない形質は，質的形質とよばれることがあります。初等的なメンデル遺伝学の例としてあげられる形質（エンドウの豆の色など）は，多くが質的形質です。これらは，環境の影響を受けず，1つないしは少数個の遺伝子で完全にコントロールされており，表現型は離散的なカテゴリで表示されます。

このblogで私は，量的遺伝学にまつわる話題を，動物育種への応用を見据えて，できるだけ平易な言葉で記したいと考えています。そのため，中には厳密でないばかりか，過剰な単純化を図ることもあるでしょう。また，多くの話題は，私の専門である育種価の予測法と分散成分の推定法に関連します。動物の育種は，遺伝評価法，選抜法，交配様式，集団の遺伝的構成，育種計画の立案，遺伝子座の発見など，極めて多くの分野から構成されています。それらは量的遺伝学，集団遺伝学，統計学，計算機科学に深く結びついています。私の理解の限界から，このblogで紹介する話題には限りがあります。必要に応じて，この分野で定番とされている（そして現在でも容易に入手できる）専門書や論文を紹介することで，厳密で詳細な議論をするかわりにしたいと考えています。

遺伝育種学の専門家が日本語で書いた，量的遺伝学の総説は以下のWebサイトで読めます。前者は主として動物を，後者は植物を専門とする研究者によって書かれています。出版されてから時間が経過していますが，歴史的経緯についての記述は現在でも十分に参考になります。

• 野村哲郎. 2000. 集団・量的遺伝学の歴史と動物育種への寄与. 動物遺伝育種研究. 28:19-27.
• 鵜飼保雄. 2011. 統計遺伝学の歴史 —とくに量的形質の解析に関連して—. 計量生物学. 32:S1-17.

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

はじめに

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

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

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

• 久保川 達也. 2011. 線形混合モデルの理論と応用 -特に小地域推定を巡って- 「21 世紀の統計科学」第 III 巻 第4章. 日本統計学会HP版. http://ebsa.ism.ac.jp/ebooks/sp/sites/default/files/ebook/1881/pdf/vol3_ch4.pdf. 分散成分の推定よりもむしろ線形混合モデル全般を取り扱うテキストである。
• 南 美穂子. 1996. 制限付き最尤推定法(REML推定法). 応用統計学, 25:73-78. https://www.jstage.jst.go.jp/article/jappstat1971/25/2/25_2_73/_article/-char/ja. REMLの導出と性質について記述されたコンパクトなテキストである。

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

• Foulley, JL. 1993. A simple argument showing how to derive restricted maximum likelihood. J. Dairy Sci. 76:2320–2324. https://www.journalofdairyscience.org/article/S0022-0302(93)77569-4/fulltext. シンプルな1変量の例をもとに、MLからREMLを導いている。これは酪農科学の雑誌だが、分散成分の推定に関する重要な論文のいくつかは、家畜育種で盛んに研究されたこともあって、農畜産系の雑誌に掲載されていることがある。
• Walsh, B. 1999. REML Variance-Component Estimation. Course Note. http://nitro.biosci.arizona.edu/courses/EEB596-1999/handouts/REML.pdf. 量的遺伝学者として著名な著者の講義ノートである。このノートの拡充版は本人による著書 Genetics and Analysis of Quantitative Traits の第27章にある。

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

• Harville, DA. (1977). Maximum likelihood approaches to variance component estimation and to related problems. J. Amer. Stat. Assoc. 72:320-338. https://www.tandfonline.com/doi/abs/10.1080/01621459.1977.10480998#.VL2LrHY8qt8. REMLに関する初期の総説で、その性質に関する重要な説明がなされている。
• Hofer, A. 1998. Variance component estimation in animal breeding: a review. J Anim. Breed. 115:247-265. https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1439-0388.1998.tb00347.x. 分散成分の推定方法について、REMLに限らず幅広く取り上げ、その応用上の留意点にまで言及している。
• Thompson, R., Brotherstone, B., White, IMS. 2005. Estimation of quantitative genetic parameters. Philos. Trans. R. Soc. Lond. B Biol. Sci. 360:1469-1477. https://royalsocietypublishing.org/doi/abs/10.1098/rstb.2005.1676. 量的遺伝学者のRobin Thompsonによる遺伝的パラメータ（分散成分）の推定に関する総説である。
• Thompson, R. 2008. Estimation of quantitative genetic parameters. Philos. Trans. R. Soc. Lond. B Biol. Sci. 360:679-686. https://royalsocietypublishing.org/doi/pdf/10.1098/rspb.2007.1417. 線形混合モデルから分散成分の推定とその量的遺伝学の応用までを取り扱った論文である。

分散成分の推定と歴史

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において利用できます。

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

• BLUPF90
• WOMBAT
• DMU
• 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.

線形混合モデルにおける分散成分の推定

はじめに

このページは、私が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）と呼んでいました。 彼らの方法は，分散分析に基づく手法よりも格段に優れていましたが，以下の点から，あまり利用されませんでした。

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

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が利用されるのは，以下のような場面に多く見られます。 主に計算上の理由であることが多いのです。

1. 方程式のサイズが非常に大きい場合（REMLでは主記憶容量と計算時間の面から，現実的ではないことが多い）
2. 数学モデルが非常に複雑な場合（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において応用されています。