幾霜::残日録::2008/01/18 (金)

　今やっている計算が全て無駄に。あーあぁーーーー。

　すげぇぇぇぇーーー。しかも作者は小笠原さんやないか! Zaurusを使い倒していた頃には小笠原さんのソフトにはお世話になったなぁ。

　WindowsではVisual C++ 2008 Express Editionを、MacOS XではXcode 3.0をダウンロードしてインストールします。

　最適化コンパイルしてバイナリを作成するには、まず自分のPCのCPUが何なのかを知らなくてはなりません。Windowsなら「コントロールパネル」の「システム」の項を、MacOS Xなら「このMacについて」を見て下さい。CPUが分かったら、SSE / SSE2 / SSE3 / AltiVecなどの拡張命令のサポート状況を検索して調べましょう。WindowsならCPU-ZとかCrystalCPUIDなどのソフトで調べることもできます。SSSE3とかSSE4というのも既にありますが上記のコンパイラではサポートされていません。

　拡張命令のサポートに応じて以下のようにコンパイルオプション(CFLAGS)を変更します。Windows用のMakefileは以前書いたエントリにあります。

Windows + Visual C++ 2008 Express Edition
without SSE
/Ox /MD /DWIN_VERSION
with SSE
/Ox /MD /arch:SSE /DWIN_VERSION
with SSE2
/Ox /MD /arch:SSE2 /DWIN_VERSION

MacOS X + Xcode 3.0
without SSE 32bit Intel
-O2 -frerun-loop-opt -fomit-frame-pointer -arch i386
with SSE 32bit Intel
-O2 -frerun-loop-opt -fomit-frame-pointer -arch i386 -mfpmath=sse -msse
with SSE2 32bit Intel
-O2 -frerun-loop-opt -fomit-frame-pointer -arch i386 -mfpmath=sse -msse -msse2
with SSE3 32bit Intel
-O2 -frerun-loop-opt -fomit-frame-pointer -arch i386 -mfpmath=sse -msse -msse2 -msse3
(64bitでは-arch x86_64)
(Core Duo / Core 2 Duoでは-mtune=nocona)
without AltiVec 32bit PowerPC
-O2 -frerun-loop-opt -fomit-frame-pointer -arch ppc
with AltiVec 32bit PowerPC
-O2 -frerun-loop-opt -fomit-frame-pointer -arch ppc -faltivec -maltivec
(64bitでは-arch ppc64)
(G3では-mtune=G3)
(G4では-mtune=G4)
(G5では-mtune=G5)

　MacOS XではMakefileの

USEREADLINE ?= yes

という行のyesをnoに書き換えて下さい。readlineはカレントディレクトリ内のファイル名補完やコマンド名補完やコマンドヒストリ機能を加えるもので、大変便利なのですがこれらを有効にする方法はめんどいのでこの辺りを見て下さい。こちらのページを参考にすればWindowsでも上記の機能が使える実行可能ファイルを作成することができます。

　Makefileができたら、Windowsなら「Visual Studio 2008 コマンド プロンプト」を起動してソースコードとMakefileのあるディレクトリ(たとえばC:\Documents and Settings\ユーザー名\My Documents\mrbayes-3.1.2)に以下のコマンドを打って移動し、nmakeでコンパイルします。

cd "C:\Documents and Settings\ユーザー名\My Documents\mrbayes-3.1.2"
nmake

　MacOS XではTerminalを起動して同様にソースのあるディレクトリ(~/mrbayes-3.1.2)に移動し、makeします。

cd ~/mrbayes-3.1.2
make

　これで実行可能ファイルができます。あとは通常のMrBayesと同様に利用できます。ただし、コンパイルオプションに起因する問題がないとは断言できません。できたバイナリの内容を保証するものではありませんのでご利用は自己の責任においてお願いいたします。

　MrBayesによる系統推定の際にもう一つ問題となるのが、「果たしてこれで十分だろうか?」「もっとMCMCを走らせるべきではないか?」ということでしょう。

　実は、MCMCからのサンプリング数が十分か否かを知る方法もTracerには用意されています。

　まず、TracerにMrBayesが吐き出した*.pファイルを読み込ませます(File → Import Trace File)。複数ある場合は読み込みを繰り返して全て読み込ませて下さい。ファイルの読み込みが終わったら、左上のペインで各ファイルのBurn-Inの値を適切に設定して下さい。Burn-Inを設定したら、「Combined」をクリックします。すると、左下のペインに各パラメータの平均値(Mean)と有効サンプルサイズ(ESS)が表示されます。ESSが100未満のものは赤で表示されます。このESSが一つでも赤で表示されていたら、まだサンプルが不足していると考えて下さい。

　Burn-Inを適切に設定していないと、誤ってESSを小さく評価してしまいますのでご注意下さい。前のエントリの上側でESSが赤く表示されているのがその例です。Burn-Inを収束してからの値に設定するとちゃんと黒くなります。また、このBurn-Inは、MrBayesとは設定方法が異なります。MrBayesのBurn-InはBurn-Inする「サンプル数」を指定しますが、Tracerでは「MCMCの世代数」を指定します。つまり、SampleFreq=100であれば、MrBayesでの設定値の100倍がTracerでの設定値になります。SampleFreq=1000なら1000倍です。

　以下の図の上側がESSが小さいパラメータのヒストグラムで、下側がESSが100を超えているヒストグラムです。ESSが100を超えるものは基本的に一つ(かそれ以上)の山型になります。上側のようにどかんとでかい頻度があるようなのはESSが小さくなります。

　前のエントリの通り、MrBayesによる計算の最中でも*.pファイルを読み込むことはできますので、計算中にときどきTracerで読み込んで収束判定やサンプリング量が十分かを判断しながら系統推定すると良いでしょう。

　MrBayesでは、収束判定のために「Average standard deviation of split frequencies」という値を出力しますが、この値だけに頼っていてはいけません。必ず、各runの尤度のプロットを比較して下さい。でないとしばしば収束していないのに収束していると誤った判断を下してしまいます。

　上の画像が具体例です。Rambautらの開発したTracerというJavaプログラムで各runの尤度の値をプロットしたものです。横軸はMCMCの世代です。上の例では約1/4が過ぎたところで両runの尤度が収束し、ASDSFもその少し後から0.01以下になっています。これはうまく収束している場合です。

　下の例は、一方のrun(青線)が開始してすぐに最尤系統樹付近に収束していますが、もう一方のrun(黒線)は、ずっとふらふら別のところを漂っています。しかし、真ん中辺りの220万世代付近で実はASDSFは0.009台になっています。一般的には収束していると考えられるラインです。これは、ASDSFが「両runが樹形空間内で同じ辺りを探索しているのか」しか表していないために間違った収束判定をしてしまう例です。つまり、黒線のrunはパラメータ空間内では全然違うところを探索しているのです。黒線の方はパラメータ空間内で局所最適解に捕まっているために尤度がこのようになってしまっていると考えられます。もちろん、青線のrunも局所最適解に捕まっていないという保証はありません。つまり、よりマシな局所最適解に捕まっているだけという可能性があります。

　このような誤判定では、得られる樹形は結局同じなんじゃないか?と思われるかもしれませんが、事後確率が不当に低くなったりします。当然、樹形自体も誤っている可能性も決して小さくはありません。また、枝長の分布はほぼ信用できないと考えて間違いないでしょう。

　MrBayesで系統推定をしていると、しばしば収束したと思ったらまたASDSFが大きくなったり、尤度が悪化したりすることがあります。しかし、尤度が最高になってASDSFも非常に小さくなった後での悪化は、むしろその樹形や尤度に捕まる方が正常なわけで、それを避けてサンプリングした場合、それは「作為」に他なりません。十分に注意して下さい。

　また、ASDSFはDiagnFreqの値によって計算する頻度が決定されます。DiagnFreqのデフォルト値は1000ですが、多遺伝子や多OTUのデータを用いた大規模解析では不十分です。大規模解析ではかなり長い周期で尤度・樹形がふらつくため、デフォルト値のままではいつまでたっても収束していないように見えてしまうことがしばしばあります(逆に誤って収束していると判定してしまうこともあります)。実際にはDiagnFreq=10000くらいにしたりするとASDSFが低い値になることがあります。そのような場合、SampleFreq(MCMCから樹形・パラメータをサンプリングする頻度)をデフォルトの100(100世代ごとに1回サンプリングという意味)よりも大きな値にする必要があるでしょう。要するに、ASDSFの値はDiagnFreqの値に依存することがあるので、DiagnFreqには「十分に」大きな値を設定する必要がある、ということです。「十分」な値がいくつになるのかはデータに依存するので私にも分かりません。

　まとめましょう。

ASDSFの値だけでなく少なくとも尤度も見て収束しているか否かを考えるべし。
ASDSFの値はDiagnFreqの値によっては参考にならない場合がある。その場合はDiagnFreqの値を大きくすべし。同時にSampleFreqの値も大きくすることが望ましい。
ASDSFおよび尤度の両面から見て一度収束したのに、再びこれらの指標が悪化するのは問題無い。むしろそれが正当である(はずだ)。そのまま解析すべし。

ということです。

　ちなみに、計算途中であっても*.pファイルを読むことはできますので、Tracerで上記のプロットを行わせることは可能です。Tracerの使い方はBEAST公式サイトのWikiを見て下さい。

追記 - 10:42:42
　一応上記の表示方法だけ書いておきます。まず、File → Import Trace Fileで*.pファイルを全て読み込ませた上で、左下のペインでLnLを選択しておきます。次に、右上のタブをTraceにして左上のペインでCtrlを押しながらファイル名をクリックすることで複数のTrace Fileを選択します。この時点で複数のプロットが出ますが色分けはされていません。ここで右下のColour byをTrace Fileに設定すると色分けされます。Legendを適当に選択すると、それぞれの色に対応するTrace Fileがわかるように凡例が表示されます。