BEAST2による分岐年代推定②

Tracerによる収束の判定

図1

BEASTによる解析が終了すると、logファイル（*.log）とtree（*.trees）ファイルができています。解析結果が収束しているかを確認するためには、logファイルをTracerに読み込ませる必要があります。Tracerを起動したら「+」をクリックするか、「File」>「Import Trace File」でlogファイルを選択します。
　
　図1がTracer の画面です（表示されるパラメータは解析の条件によって異なります）。この解析では200,000,000回の試行を行っています（図1の「States」）。また、MCMCでは、試行の最初の方は数値が安定しないため結果から省きます。これが「Burn-In」で、この例では全試行の10%にあたる20,000,000回が指定されています。
　

図2

続いて各パラメータが収束しているか見ていきます。例として「birthRate」を見てみます。左側のパネルで「birthRate」をクリックし、右側のパネルで「Trace」タブを表示してください（図2）。表示されたグラフの横軸は試行回数で、縦軸は各試行での推定値です。グラフの左端の白くなっている部分はBurn-Inによって省かれたところです。このグラフでは、Burn-Inによって削られて以降、推定値が一定の範囲に収まっているのがわかるかと思います。次に左側のパネルで「ESS」（有効サンプルサイズ）を見ると、「birthRate」では1714と大きな値となっています。この値はパラメータの収束に関する基準で、200を超えるようにすべきとされています。以上より、「birthRate」は十分に収束していると判断できます。

図3

　悪い例を見てみます。これは「prior」 のパラメータですが、グラフを見ると試行により推定値が一定の範囲に収まっていないことがわかると思います（図3）。ESSの値も61と低く、収束しているとは判断できません。このような場合、試行回数をさらに増やしたほうが良いでしょう。試行回数を増やしてもESSが改善されない場合には、パラメータの初期値などを変更したほうが良いかもしれません。


Treeファイルの要約

全てのパラメータの収束が十分であると判断されたら、treeファイルの要約を行います。BEASTの解析で得られたtreeファイルには各試行の結果が含まれています。すなわち、先の例でいえば、200,000,000 / 1000（どのくらいの頻度でファイルに記録するかという値。xmlファイル作成時に変更可能。）= 200,000本分の系統樹と分岐年代推定値が記録されています。当然一つ一つ見ていくわけにはいかないので、treeファイルの内容を要約します。これにはTreeAnnotatorを用います（BEASTと同じフォルダに入っています）。

図4

TreeAnnotatorを起動します（図4）。「Input Tree File」でtreeファイルを選択します。「Output Tree File」では出力するファイル名を入力します。「Burnin percentage」では何%削るかを指定します。これはTracerにおいて収束が確認された際のBurninに合わせます。200,000,000回の試行、20,000,000回のBurninで収束しているならば、10と入力してください。「Target tree type」はどの樹形で出力するかを指定します。基本的にはそのままで大丈夫でしょう。「Node heights」は「Mean heights」 にしておけば良いと思います。入力ができたら「Run」をクリックします。要約が完了したら出力ファイルをFigTreeに読み込ませて結果を確認してください。以上がBEASTによる分岐年代推定の基本的な方法となります。

2021/12/14 追記

「prior」が極端に低い場合の対応について

ほとんどのパラメータは収束しているのに「prior」と「posterior」のESSが低いままのことがしばしばあります。「posterior」は「prior」と「likelihood」に関係してくるので、実質「prior」のみの問題と考えて良いでしょう。グラフを見ても収束の兆しが見えず、試行回数を増やしても改善されない場合、推定するパラメータ数が多すぎることに起因している可能性があります。対処方法としては、とにかく推定するパラメータ数を減らすことです。例えば、Tree priorとしてBirth-Death modelを使っているならYule modelへ変更したり、インプット配列のパーティション数が多いのならばPartitionFinderなどを用いてパーティション数を減らす、また、塩基置換モデルとしてGTRなどの複雑なモデルを使用している場合、よりシンプルなHKYモデルにするなどすることで収束しやすくなります。