目次

• はじめに
• 簡単にいうと…
• 事後確率
• マルコフ連鎖モンテカルロ法
• MrBayes による系統解析の仕組み
• 問題点と今後の展望
• 参考文献・サイト

最近ではベイズ法を用いて（MrBayes を用いて）系統解析を行う論文も大分増えており， 方法論としてはほぼ定着した感があります。 しかしその原理や考え方についてはあまり広く理解されていないようでもあります。 そこで，筆者が幾つかの文献や資料を読んで理解した範囲で，その考え方を紹介しておきたいと思います。 筆者は数学（統計学）を専門にしているわけではありませんので，数学的に厳密な解説ではないかもしれませんが， その程度の解説と思って読んでいただければいいと思います。 正確・詳細な解説については末尾の参考文献・サイトを当たってみてください。 なお，間違いなどがありましたら，ご指摘のほどよろしくお願い致します。

最尤法では通常，発見的探索法（heuristic search）という方法などで，多数の系統樹について尤度を計算し， その中から尤度が最も高いものを選択します。 ところがベイズ推定においてはこのような探索法は使えません（各系統樹について事後確率を計算することが出来ないため）。 代わりにプログラムを用いたベイズ推論では，マルコフ連鎖モンテカルロ（MCMC）法なる方法を用います。 MCMC 法ではある種の規則に従って大量の系統樹を作ります。その中で様々な単系統群が出現する頻度を求めると， これがその単系統群の事後確率の近似になります。 こうして求めた事後確率が 50% 以上の単系統群をまとめた系統樹を作るのがベイズ法です。

その他，最尤法による系統樹とベイズ法による系統樹の大きな違いとして，系統樹の統計的な評価法があります。 得られた系統樹上の各枝の単系統性を評価するために， 最尤法ではブートストラップ値（通常 % 表記）が用いられるのに対して， ベイズ法でブートストラップ値を用いるのは現実的でなく， 事後確率（通常 % 表記はしない）によって各枝の単系統性を評価します。

また，最尤法と対比してみます。最尤法で系統樹を測る尺度となっている尤度ですが， これはある樹形（仮定）に沿って遺伝子が進化したときに，現在の遺伝子配列のセットが実現する確率を表します。 遺伝子の進化には無数の可能性が考えられますから，現在の配列データが実現する確率は極めて低いものになります。 一方で，ベイズ法における事後確率は，尤度と逆の発想になっています。尤度が「仮定→データ」の確率だとすると， 事後確率は「データ→仮定」の確率，すなわち手元にある遺伝子の配列， アラインメントのデータの原因が，ある系統樹である確率と言えます。

具体的な例を挙げて考えて見ましょう。 ここでは系統樹の作成ではなく，植物の同定を最尤法，あるいはベイズ推論を用いてすることを考えます。 さて，ある公園で見つけた花 X の名前を調べるとします。花 X は以下のような特徴を持っていました （長さについては四捨五入した値とします）。

単純に考えるため，花 X は種（しゅ）P，Q，R のいずれかであるとします。 種 P〜R については各形質を持っている確率だけがわかっています。 例えば，種 P の花弁の長さが 3 cm である確率は 10 %，4 cm の確率は 15 %，5cm の確率は 50 % ・・・， といった感じです。問題になるのは，花 X と同じ形質を持っている確率だけなので，以下にその確率の表を示します。 （花弁の長さであれば，5 cm である確率，つまり 50 %）。

最下段に示した尤度はこの場合，各種について全ての形質が花 X と同じである確率で， 各形質の確率を積算することで求められます（種 P の場合は 50 % × 30 % × 30 % × 100 % × 50 % ＝ 2.3 %）。 つまり，種 P の 2.3 % が花 X と区別できず，種 Q の 8.7 % が花 X と区別できない，ということです。 ちなみに種 R は花弁の枚数で確実に花 X と区別できるため，尤度が 0 % になっています。 図にすると下図の赤い部分の割合が（花 X が各種である）尤度ということになります。

ここで，尤度が最も高いのは種 Q の場合ですから，最尤法の考え方に従うと， 花 X は種 Q と考えるのが最も尤もらしい（もっともらしい）推論として採用されることになります。 これは一つの考え方ですが，場合によっては最尤法の考え方がふさわしくない場合があります。

事後確率の考え方で上の例のを考えて見ましょう。ここでは，花 X が見つかった公園に， 種 P，Q，R がどれくらいの個体数あるのかが事前にわかっていたものとします（下図）。

この例の場合，今示した頻度の割合を事前確率と呼びます。 それぞれの種に事前確率をかけた値は，全体の中で例えば種 P であって花 X と同じ特徴を持った個体の割合になります。 言い換えると，公園の中で 100 本花があったときに，花 X と同じ特徴を持った種 P の花は，1.8 本存在すると言えます （2.3 % × 78 % = 1.8 % なので）。同様に他の 2 種についても計算してみましょう（下図）。

この結果，花 X と同じ特徴を持った個体数は，種 Q よりも 種 P の方が大きいことがわかります（1.8 % > 0.17 %）。 もう少し考えると，花 X が種 P である確率が計算できることがわかります。 すなわち図の赤い面積の合計の中で，種 P に含まれている部分がどれくらいの割合なのかを計算すればいいわけです。 これは以下の計算式で計算され，事後確率と呼ばれる値になります。

同様に他の 2 種についても事後確率を計算することが出来ます。

この結果から，花 X は 91 % の確率で種 P に同定され，種 Q の確率も 8.6 % ながら存在するが，種 R の確率はない， ということが言えます。これがベイズ推論の考え方です。 なお事前確率が全て等しい場合（種 P，Q，R ともに 33 % の場合）， ベイズ推論での結果と最尤法で選ばれる結果は厳密に一致します。この例の場合，事前確率が等しければ，

となり，やはり花 X は種 Q である確率が一番高くなります。

さて，このベイズ推論を系統樹の推定に導入するとどうなるでしょうか？ Huelsenbeck et al. (2001) に示された式を引用すると，あるデータの下でのある樹形が正しい確率（Pr[Tree | Data]）は， 以下の形式で示されます。

ここで，Pr[Data | Tree] はその樹形の下での尤度で，Pr[Tree] はその樹形の事前確率になります （通常の系統解析では各樹形が出現する確率は平等と考え，全て同じ値にしています）。

さて，ここで問題になるのが Pr[Data] の値です。これは上で示した花の同定の例と同様に， 「事前確率×尤度」の総和になります。総和と書きましたが，これは花の場合には 3 種の候補についての和を求めただけで済みました。 しかし系統樹の場合には考えられる全ての樹形について計算しなければなりません。 しかも各樹形について枝の長さについて積分したり（連続値なので和ではありません）， 塩基置換速度などのパラメータの値についても全て考慮しなければなりません。 当然ながらそのような値は計算できません。

これが，系統樹などに事前確率を導入するときの問題となっていました。 これを解決するために現在用いられているのが MCMC 法になります。

マルコフ連鎖モンテカルロ（MCMC）法とは，簡潔に言えばマルコフ連鎖に対してモンテカルロ法を適用する方法です。 といってもこれでは意味が通じないでしょうから，もう少し踏み込んで解説してみましょう。 まずモンテカルロ法というのは，ランダムな試行を繰り返して目的の値を近似的に測定する方法です。 例えば，円の面積 M を近似的に求めることを考えましょう（計算を簡単にするため，円の 1/4 で考えます）。 以下のような扇形と正方形を考えます。 この正方形の中にランダムに無数の点を発生させます（X 軸の値と Y 軸の値を乱数で決めます）。 大量の点のうち，扇形の中に入っている点（）の割合を数えれば, 扇形の面積の近似になると考えられます。

実際に計算して見ますと，理論値では M = 0.785 に対してモンテカルロ法（n = 5000；n は試行の数） で計算した値は M = 0.778 と，比較的近い値が与えられています。 点の数を増やせば増やすほど，この値は真の値に近づくと考えられます。

通常のモンテカルロ法と異なり，MCMC 法では一回ごとの試行が独立ではなく， ある試行の結果（の確率分布）が前の試行の結果に影響されるような現象について，モンテカルロ法を行います。 例えば，4 人でゲームをしている状況を考えてください。あるゲームで勝ったプレイヤーは次のゲームで親となり， 勝ちやすくなるとしましょう。これを続けていくような試行をマルコフ連鎖と呼びます。 これを無限に繰り返した場合，それぞれのプレイヤーの勝率がどうなるのかを MCMC 法で求めることが出来ます （注：この仮定の場合，MCMC 法を使わずとも解析的に勝率を計算することが出来ます）。

具体的な数値を入れて考えて見ましょう。プレイヤー A，B，C，D がいたとして， それぞれゲームに上手下手があるとします。また親のときは，そうでないときより 20 % 勝率があがります。 勝率は以下のように決まっているとしましょう（なお，最初に親になる確率は平等に 25 % ずつとします） （注：親の時の確率は太字で示してあります）。

このような場合，充分に試行を行った後の A の勝率はどのようになるでしょうか？ これを MCMC 法で計算してみます。 まずは今記述したような条件でのシミュレーションを延々と繰り返します。 試行数が少ないうちは，初めに誰が親だったかの影響が残りますが， そのうちに最初に誰が親だったのかは関係なくなります。これを定常状態に達したと表現しますが， この，定常状態に達して以降に A が親になった回数を数え，その割合を計算すれば目的の確率が計算できます。 実際にシミュレーションを走らせた結果を下図に示してみます。表示してあるのは親になったプレイヤーのみです。

さらにここから，定常状態（ここでは 1001 回目〜 5000 回目をカウントしました） でそれぞれのプレイヤーが親になる確率を求めました。理論値も参考に示しておきましょう。

このように MCMC 法の結果は理論値と比較的近い値になっています。

と，これが MCMC 法の基本的な考え方ですが，この方法は事後確率の計算にも利用することができます。 なぜ、MCMC 法を事後確率の推定に用いることができるのかについては当然ながら数学的な背景があるようですが、 ここでは簡単な説明にとどめることにします（というか、詳細は私が理解していません。ごめんなさい）。 事後確率を求める場合には、Metropolis-Hastings のアルゴリズム、なるものを用います。 これは、ある試行（世代とも言う）の結果から次の試行の結果（の確率）を求めるルールです。 ある世代の結果、いくつかのパラメータ（例えば樹形や枝長）が得られたとします。 このとき次の世代を求めるために、これらのパラメータをランダムに変化させます。 変化させた結果と変化する前の値について、尤度、事前確率、互いの値から次の世代に相手の値が導かれる確率、 をそれぞれ求め、これらの値から公式に従って、新しい値を受理するか棄却するかの確率を決定します。 受理する場合は　「新しい値」＝「次世代の値」　となります。棄却する場合には　「もとの値」＝「次世代の値」　 となります。 このように作られたマルコフ連鎖について MCMC 法を適用すると、事後確率が求められるという仕組みになっています。

MrBayes などの系統解析プログラムにおいては， 系統樹の樹形，枝の長さ，塩基置換速度に関するパラメータなどを求めることになります。 最初はランダムな系統樹からスタートし、世代が 1 個進むごとに，系統樹を一部変化させます。 この時に上記の確率に基づいて新しい系統樹を次の世代に採用するか，それとも棄却するかを決定します。 棄却する場合には前の世代の系統樹をそのまま次の世代に採用することになります。

なお系統樹の場合の尤度とは，最初に述べた「ある樹形（仮定）に沿って遺伝子が進化したときに， 現在の遺伝子配列のセットが実現する確率」です。従って，これを計算するには塩基などの進化のモデルが必要です。 つまり，どの塩基間でどの程度の確率で置換が起きるかの推定を求めなければならないのです。 具体的な値については MCMC で同時に求めるのですが，パラメータの数については事前に決める必要があります。

さて，このようなマルコフ連鎖をやはり延々と繰り返していきます。 これも同様に連鎖が定常状態に達するのまで続けて， 定常状態にある大量のサンプルから樹形などのデータをサンプリングしていきます。 すると，ここでの樹形の出現頻度がその樹形の事後確率になるそうです。

これは各単系統群についても言えるので，実際には出現頻度（事後確率）が 50 % を超える単系統群（枝）を求めます。 そのような枝を集めて 50 % 合意樹を描くことで，最も事後確率の高い系統樹（に近い系統樹）を描くことが出来ます。 簡単に言うと，これがベイズ法によって系統樹を描く方法となります。ちなみに事前確率が等しい場合， これは最尤系統樹の近似にもなっていると考えられます。

ベイズ法を用いた基本的な系統解析方法は前章で紹介したとおりですが， 実際の解析には MrBayes というプログラムが最もよく用いられています（Huelsenbeck & Ronquist, 2001; Ronquist & Huelsenbeck, 2003）。 MrBayes では MCMC 法にもう一工夫が凝らされています。ここではその，Metropolis Coupled MCMC（MCMCMC，あるいは (MC)³；ちなみに Metropolis はこの方法を原形を開発した研究者の名前です）について説明を補足します。 また，定常状態に達していることを評価する方法や，定常状態に達して以降のサンプリングが充分かどうかも， 最新のバージョンでは判定できるようになりました。これについても解説したいと思います。

MCMC 法で系統解析を行う場合，尤度がある程度高くなると， その付近の系統樹，「局所的な最適解」にマルコフ連鎖がつかまってしまう場合があります。 これでは正しい定常状態に達することが出来ませんので，系統樹を導くことが出来ません。 そこで局所的な最適解につかまることを防ぐために導入されたのが，MCMCMC 法です。

MCMCMC 法では複数の（しばしば 4 本の）マルコフ連鎖を同時に走らせます。 この内 1 本は cold chain と呼ばれ，これが通常の MCMC に相当します。サンプリングはこの cold chain から行います。 cold chain 以外の他の連鎖は heated chains と呼ばれ，それぞれ異なる程度で cold chain より激しく系統樹が撹乱されます。そして各世代ごとに（あるいは数世代ごとに）ランダムに 2 つの連鎖の状態が入れ替えられます（下図）。

この MCMCMC 法を用いることによって，マルコフ連鎖は局所的な最適解にとらわれにくくなり， 特に heated chains によって効率的に系統樹の可能性が探索されることになります。 実際に MrBayes を使用する場合には，撹乱の程度や連鎖の本数なども変更することが出来ます。 また，定常状態におけるサンプリングの頻度は，図では毎世代サンプリングしているように描きましたが， 実際には設定可能で，通常は 100 世代ごとにサンプリングしているようです。

さて，サンプリングを行うときには，マルコフ連鎖が定常状態になっているかどうかが非常に重要な意味を持ちます。 これを理論的に計算することは出来ないため，MrBayes 3.0b4 以前のバージョンでは尤度の変化を見ていました。 尤度（の対数値）のグラフは普通以下のようになります。

今までのやり方では，このグラフやあるいは尤度の値を見て尤度が一定の値の近くで余り変化をしなくなれば， マルコフ連鎖が定常状態に達した（収束した）として，この状態のマルコフ連鎖からサンプリングを行ってきました。 ところがこれは実は収束を正確には再現していないことが最近になって示されてきています（Mossel & Vigoda, 2005： ベイズ法の解析の問題点も参照）。 そこで，収束をより客観的，定量的に測定する方法が必要でした。 このことは以前からも考えられており，MrBayes 3.1 以降では収束を判定するために， Average Standard Deviation of Split Frequencies（以下，ASDSF）という値を用いることが出来るようになりました。

ASDSF はどのような考え方から出てきた値なのでしょうか？ この概念の理解に先立って，単系統群（および側系統群）を区画分けして表現する方法について触れておく必要があります。 無根系統樹は任意の枝の上を境目にして，常に二つに分けることが出来ます。 これに併せて系統樹の末端の分類単位も二つに仕分けすることが出来ます（下図；太線を境に (A,B,C)と(D,E,F,G,H)）。

この仕分けを文字で表記すると，次のようになります。

左三つのピリオドは，A〜C が同じ区画に含まれることを意味し，右の＊は残りがもう一つの区画を作ることを示します。 同様に，この系統樹は 13 個の区画分け（つまり枝）によって記述できます。

記述方法はさておき，MCMC 法ではこのような区画分けをサンプリングし，それぞれの出現頻度を求めます。 マルコフ連鎖が定常状態に達すると，この各区画分けの出現頻度が事後確率に従うと考えられます。

逆に定常状態に達しているかどうかの判定に，この事実を用いようとするのが ASDSF の考え方です。 どのような初期条件で何度 MCMC を走らせようと，同じ配列情報と進化モデルのもとでは， 最終的な事後確率の分布は同一であるはずです。ということは，つまり様々な区画分けの出現頻度が等しくなるはずなのです。 このことを逆に利用して，二つ以上の MCMCMC で区画分けの出現頻度が一致するようになれば， 両方の MCMCMC が定常状態にあると考えてよいわけです。

MrBayes 3.1 以降では二つ以上の MCMCMC を同時に走らせることができるようになりました。 これらの MCMCMC について適当な世代ごとに任意の burn-in（まだ定常状態に達していない世代と仮定して， 計算から切り捨てる部分）を設定し，それぞれの MCMCMC における様々な区画分けの出現頻度の標準偏差を求めます。 これを平均したものが ASDSF に相当します。

仮に burn-in 以降になっても定常状態に達していなかったとすると， 各 MCMCMC ごとに区画分けの出現頻度は異なってくるため，標準偏差が大きくなり，ASDSF は大きくなります。 一方で，設定した burn-in より後で MCMCMC が定常状態に達していたならば， 全ての区画分けの出現頻度はどの MCMCMC でも等しいはずで，標準偏差は全て 0 に近い値をとるはずです。 当然その平均である ASDSF の値もほぼ 0 になります。 このように，ASDSF が 0 に近づくことを指標にして，MCMCMC の収束が客観的に測定できるようになったのです。 具体的な数値としては ASDSF が 0.01 を切るのが一つの目標値と見ると良いでしょう（MrBayes 3.1 のマニュアル）。

もう一つ解析を行う際に問題になるのは，定常状態において一体どれくらいの数のサンプリングを行えばよいか， ということです。これを判断する指標として用いられるのが PSRF（Potential Scale Reduction Factor）という値です。

PSRF については MrBayes の作者もあまり詳細な解説を与えておらず， 単にパラメータの値が充分に収束しているのかを判断する基準とだけ紹介されています。 この値は，MCMC からのサンプリングが充分であれば 1 に近づくとされており， PSRF 値が 1 から外れている場合にはサンプリングを増やす必要があります （マルコフ連鎖からのサンプリング密度を増やしたり，長く連鎖を実行したりします）。

PSRF 値はマルコフ連鎖から取り出した結果が Student の t 分布に従っているかどうかの指標のようですが， 正確な式などについては PSRF 値を導入した論文（Gelman & Rubin, 1992）などを参照してください。なお，MrBayes 3.1 の解説によると，系統解析の場合は PSRF 値を用いるための条件が完全には満たされていないとのことで， 収束の厳密な指標にはならないため，あくまで粗い推定に用いるべきだとされています。

MCMCMC の収束の評価については，以下のようにまとめられます。

ここまでの解説では，ベイズ法にまつわる問題については触れてきませんでしたが， 実際の系統解析の場合には様々な問題点が知られてきています。 問題点を理解することでより適切にベイズ法を使用することができるようになるでしょう。

ベイズ法について指摘されている最も深刻な問題点は，MCMCMC を用いた場合に事後確率が過大評価されるという可能性です。 Suzuki et al. (2005) によるシミュレーション研究ではこの点が詳細に議論されています。 彼らは 4 種の仮想上の遺伝子配列を構築，これらの系統関係を最尤法，近隣結合法， ベイズ法で解析し，前 2 者についてはブートストラップ値で，ベイズ法については事後確率で信頼性の評価を行いました。 この結果，ブートストラップ値による評価が適切であったのに対して，事後確率による評価は過大評価となっていました。 この一因として，ブートストラップ値では遺伝子配列が有限の長さしか用いられないために起こる統計的誤差を考慮しており， 事後確率ではこれは考慮されていないことが指摘できます。

しかしながらベイズ法の事後確率では，次善の系統樹やそれ以外の系統樹と比べて， 得られた系統樹がどの程度優れているのかも判断することが可能です。一方で最尤法を用いた場合は， ブートストラップ値の比較とは独立に，尤度比検定を個々の（問題の）系統樹について行わなければなりません。 ですから事後確率とブートストラップ値は全く異なる統計的誤差を扱っているとも言えると思います。 そのためシミュレーションの行い方がそもそも適切であるのか，まだ検証の余地があるかもしれません。

一般にはブートストラップ値の 50 % を系統樹に示す場合，事後確率では 90 %（0.90） 以上の値を示す慣習があるようですので，全く異なる指標と考えて取り扱えば問題ないと考えます。

最尤法など他の系統推定法と，ベイズ法では系統樹の構築の方法が根本的に異なっていますから， いずれかに特定のバイアスがかかる場合もあるでしょう。そのようなバイアスを回避する目的では， ベイズ法と他の方法を併用して整合性が取れるかどうかを比較することが重要となるでしょう。

さて，他の問題として最近指摘されているのは，heterotachy というバイアスを回避できるのかどうかという点です。 heterotachy については雑記中に幾つか記事を書いていますので，そちらも参照してください （最尤法か最節約法か，続報， 続報 2，続報 3）。 Heterotachy とは，サイトや領域ごとに，進化速度の系統樹上での変化が異なっている場合を指します。 これは特に，複数のドメインをもつタンパク質の場合や，複数タンパク質をつなげて解析する場合などに問題になります。

もう一点，技術的な問題について一点だけ触れて起きます。 実際に解析を行う場合にはマルコフ連鎖が中々収束しない場合が出てくるようです。 このような場合には系統樹を得ることが出来ませんから，何らかの対策を打つ必要が出てきます。 MrBayes 3.1 のマニュアルに与えられている解決法としては，まず解析を長くすることが挙げられています。

次に MCMC の際に系統樹（やパラメータ）を撹乱する方法（proposal: "prop" コマンドで変更可）が修正 できます。解析の後に，cold chain において新しい状態が受理された確率が以下のような形で示されています。

これは 10 % から 70 % の間にあることが望ましいとされており，前述の "prop" コマンドを用いて proposal を変更することで調整することが出来るでしょう（詳細はマニュアルの Appendix も参照）。 ただし不適切な調整を行うと，収束の可能性を減らすため，慎重に行う必要があるそうです。 目安としては，各世代につき一つのパラメータ変更が行われることが望ましいとされています。

複数のマルコフ連鎖を走らせた場合，隣接する連鎖の間で受理の確率が低すぎる場合も問題があり， MCMCMC の温度（temperature）を変更する必要があります。これはやはり解析の後に表示される以下の表の， 赤丸で強調した部分の値を見ることで判断できます。

ここでは 1 が cold chain を指し，残りは 2，3，4 の順に heat の程度が上がります。 従って上記の赤丸の数値が低い場合は，温度（撹乱の程度）が近いもの同士の間での受理の確率なので， この値が低いと MCMCMC の効率に影響すると考えられます。これは，temperature を下げることによって調整できます （"mcmc temp=**" あるいは "mcmcp temp=**" で設定できます。 また，コンピュータの能力に余裕があれば連鎖数を増やすことも一つの手だそうです。

収束を早めるための手として，他にはランダムな樹形からマルコフ連鎖を始めるのではなく， （例えば他の解析法で描いた）最初の樹形にまともな系統樹を用いる方法も用意されています （初めに "usertree = (A, B, (C, D))" のような形式で樹形を読み込ませ， マルコフ連鎖の設定時に "mcmcp startingtree=user" として利用します）。 しかし複数の MCMCMC を比較する場合（ASDSF を用いたい場合），同じ系統樹から始めては意味がありません。 そのため初めに設定した系統樹を撹乱しておくという方法も用意されています（"mcmcp nperts=<撹乱の回数>"）。

現時点でははブートストラップ値に相当する系統樹の信頼度の指標が得られないことが， ベイズ法の欠点の一つとなっていますが，新しい方法が発明されれば， ブートストラップ値と同様にデータの統計的誤差を判断する指標が使えるようになるかもしれません。

その他，系統樹の事前確率については現在考慮されていませんが， ランダムな系統分岐が起こったようなケースなのか，特定の系統で多様化が進んだケースなのか， 事前に考えることも出来るでしょう。 このような系統パターンの進化モデルも充分な理論的な研究がなされれば， ベイズ法に導入していくことも出来るのではないでしょうか。

などなど，まだまだ系統解析に用いられてから歳が浅いだけに，今後の発展性にも多くの期待が持てるでしょう。 一方で，未熟な方法であることも事実ですから，常に他の解析法と併用し， 極度におかしなことが起こっていないようチェックすることも必要でしょう。