----------------------------------------------
●最尤法でのモデル選択問題（ModelTestを用いて）
----------------------------------------------
分子系統学で最尤法を用いる際には，どのような分子進化モデルを仮定するかによって結果が左右されます．具体的には，塩基置換モデルやアミノ酸置換モデルをどのように経験的に選択するのかという問題です．統計学の世界では「モデル選択問題」はそれ自身が大きな研究領域となりつつあり，分子系統樹の構築を考えるときにもそれは関わりをもってきます．
ここでは塩基置換モデルに焦点を当ててモデル選択問題を考えてみます．一般時間逆転モデル（ＧＴＲ）をもっとも緩いモデル，Jukes-Cantorモデル（ＪＣ）をもっとも厳しいモデルとすると，現在用いられている塩基置換モデルはＧＴＲとＪＣを両極端とする階層のいずれかの場所に位置づけられます．
講義の中でデモした David Posada のソフトウェア「ModelTest」：
ModelTest Home Page
http://inbio.byu.edu/Faculty/kac/crandall_lab/modeltest.htm
は，これらの塩基置換モデルから「階層的尤度比検定（hLRT）」と「赤池情報量基準（AIC）」に基づいて自分のデータから絞りこむためのものです．
----------------
階層的尤度比検定
----------------
塩基置換モデル群が「階層的」に構造化されているとは，モデルのパラメーターに関する複雑性の間に順序関係があり，より単純なモデルＭ0を帰無仮説とし，より複雑なモデルＭ1を対立仮説として尤度比に基づく比較を行なうことになります．与えられた塩基配列データＤのもとでの最大尤度をそれぞれ
maxＬ(Ｍ0｜Ｄ) 帰無モデル（ｒ個のパラメーターを制約づける）
maxＬ(Ｍ1｜Ｄ) 対立モデル（ｒ個のパラメーターは無制約）
とするとき，尤度比：
λ＝maxＬ(Ｍ0｜Ｄ)／maxＬ(Ｍ1｜Ｄ)
の対数値δ＝−２logλは，漸近的にある非心χ二乗分布にしたがう統計量です．とくに，帰無モデルＭ0のもとでδは非心度ゼロのχ二乗分布に漸近的にしたがい，その自由度は両モデルの自由パラメーター数の差ｒとなります．この性質を用いて，モデル間の尤度比検定を行ないます．（尤度比の統計学的性質については，Stuart & Ord 1991, Kendall's Advanced Theory of Statistics, vol.2. §23.7参照）
ModelTest では，PAUP* を併用することにより，この尤度比検定を実行します．テストデータ（model3.nex）：
-----
#NEXUS
BEGIN DATA;
DIMENSIONS NTAX=10 NCHAR=30;
FORMAT DATATYPE=DNA ;
MATRIX
'Taxa A' GCAATCATCCAAATCGGTCAACTTAATAAA
'Taxa B' GCCATCATCCATAACGGTGAACTTTTAATG
'Taxa C' GCCATACTCCATAACGGTGAACTTGTAATA
'Taxa D' GCCAAACCCCATATCGTGCAACTTAATAAG
'Taxa E' GGCTATCCACCTTAAGTGTAAATTGTTGAT
'Taxa F' GGCTATCCAACTATAGTGCAACTTAATACA
'Taxa G' GGCTAGGCCAATAATATGAAACTTTTAATG
'Taxa H' GTCTAGGCCAAAAATATGAAACTTGTTATA
'Taxa I' GTCGAAGCAAAAATAGTGCAACTCAATAAA
'Taxa J' GTCGAAGCAAAAATAGTGAAACTCAATAGA
;
END;
-----
をまずPAUP*に取り込ませた後，ModelTestのコマンドファイル（modelblock3.nex）：
-----
#NEXUS
BEGIN PAUP;
log file= modelfit.log replace;
DSet distance=JC objective=ME base=equal rates=equal pinv=0
subst=all negbrlen=setzero;
NJ showtree=no breakties=random;
End;
[!***** BEGIN TESTING 56 MODELS OF EVOLUTION ***** ]
BEGIN PAUP;
Set criterion=likelihood;
[!** Model 1 of 56 * Calculating JC **]
lscores 1/ nst=1 base=equal rates=equal pinv=0 scorefile=model.scores replace;
[!** Model 2 of 56 * Calculating JC+I **]
lscores 1/ nst=1 base=equal rates=equal pinv=est scorefile=model.scores append;
...（中略）...
[!** Model 5 of 56 * Calculating F81 **]
lscores 1/ nst=1 base=est rates=equal pinv=0 scorefile=model.scores append;
...（中略）...
[!** Model 56 of 56 * Calculating GTR+I+G **]
lscores 1/ nst=6 base=est rmat=est rates=gamma shape=est pinv=est
scorefile=model.scores append;
-----
を PAUP* でオープンして実行させます．このコマンドファイルは，事前に取り込んである塩基配列データ（model3.nex）から求められた近隣結合樹（距離はＪＣで求める）の樹形に対する尤度を，ＪＣを出発点とするモデル階層に含まれる計５６モデルのそれぞれについて計算します．
PAUP*で計算された尤度値やその他のパラメーター推定値を含む出力ファイル（model3.scores）を ModelTest に入力し，尤度比検定を行なわせます．
-----
Testing models of evolution - Modeltest Version 3.06
(c) Copyright, 1998-2000 David Posada (dp47@email.byu.edu)
Department of Zoology, Brigham Young University
WIDB 574, Provo, UT 84602, USA
_______________________________________________________________
Tue Feb 18 22:33:34 2003
Input format: Paup matrix file
** Log Likelihood scores **
+I +G +I+G
JC = 254.8371 254.8371 252.2073 252.2073
F81 = 248.8631 248.8631 245.8168 245.8168
K80 = 253.0150 253.0150 250.3367 250.3367
HKY = 247.5020 247.5020 244.4027 244.4027
TrNef = 251.4816 251.4816 249.2190 249.2190
TrN = 245.8147 245.8147 243.2256 243.2256
K81 = 253.0004 253.0004 250.3363 250.3363
K81uf = 247.4607 247.4607 244.2478 244.2478
TIMef = 251.4680 251.4680 249.2189 249.2189
TIM = 245.7733 245.7733 243.0745 243.0745
TVMef = 249.7966 249.7966 247.2471 247.2471
TVM = 246.5513 246.5513 243.6032 243.6032
SYM = 248.3923 248.3923 246.0338 246.0338
GTR = 245.0729 245.0729 242.4911 242.4911
** Hierarchical Likelihood Ratio Tests (hLRTs) **
Equal base frequencies
Null model = JC -lnL0 = 255.5512
Alternative model = F81 -lnL1 = 249.1400
2(lnL1-lnL0) = 12.8225 df = 3
P-value = 0.005037
Ti=Tv
Null model = F81 -lnL0 = 249.1400
Alternative model = HKY -lnL1 = 247.7302
2(lnL1-lnL0) = 2.8196 df = 1
P-value = 0.093121
Equal rates among sites
Null model = F81 -lnL0 = 249.1400
Alternative model = F81+G -lnL1 = 245.8168
2(lnL1-lnL0) = 6.6465 df = 1
Using mixed chi-square distribution
P-value = 0.004968
No Invariable sites
Null model = F81+G -lnL0 = 245.8168
Alternative model = F81+I+G -lnL1 = 245.8168
2(lnL1-lnL0) = 0.0000 df = 1
Using mixed chi-square distribution
P-value = >0.999999
Model selected: F81+G
-lnL = 245.8168
Base frequencies:
freqA = 0.4043
freqC = 0.1639
freqG = 0.1654
freqT = 0.2665
Substitution model:
All rates equal
Among-site rate variation
Proportion of invariable sites = 0
Variable sites (G)
Gamma distribution shape parameter = 1.6296
-----
上の尤度比検定は，JC対F81に始まり（F81採用），F81対HKY（F81採用），F81対F81+G（F81+G採用）と有意な尤度差が出るかぎり検定を進めます．続くF81+G対F81+I+Gの検定で有意差が出なかったため，与えられたデータからは，定常塩基頻度の同一性を仮定しないF81（Felsenstein 1981）モデルにガンマ分布による置換速度のバラツキ補正をしたモデルが最適であると判定されます．
この解析結果を踏まえた PAUP* での最尤法のオプション設定が下記の通り出力されるので，再び PAUP* に復帰することで最尤系統樹の計算に進むことができます．
-----
PAUP* Commands Block: If you want to implement the previous estimates as likelihod settings in PAUP*, attach the next block of commands after the data in your PAUP file:
[!Likelihood settings from best-fit model (F81+G) selected by hLRT in Modeltest Version 3.06]
BEGIN PAUP;
Lset Base=(0.4043 0.1639 0.1654) Nst=1 Rates=gamma Shape=1.6296 Pinvar=0;
END;
--------------
赤池情報量基準
--------------
赤池情報量基準（ＡＩＣ）に基づくモデル選択では，与えられたデータのもとでの最大尤度ではなく，データの分布に基づく期待最大尤度をモデルの判定基準とします．
データの分布はもともと未知であるため，本来ならば期待最大尤度は認識論的には知ることができないのですが，尤度関数のテーラー展開による線形化と漸近的なχ二乗近似による期待値計算を適用することにより，期待最大尤度は認識論的に「最大尤度−パラメーター数」の期待値として求められます．ＡＩＣはこの期待最大尤度の−２倍として定義されます．
（ＡＩＣに関する数理統計学の論議は坂元慶行他 1986『情報量統計学』共立出版に詳しいですが，死亡率はきっと高いでしょう．ナミダなしのＡＩＣ？――これなら生き残れるぞ！ Forster & E. Sober 1994. Brit. J. Phil. Sci., 45: 1-35 の中の「Akaike without tears（pp.2-11）」を拝むべし．）
さて，ModelTest の出力リストの後半はＡＩＣに基づくモデル選択の結果が示されています：
-----
** Akaike Information Criterion (AIC) **
Model selected: TrN+G
-lnL = 243.2256
AIC = 498.4513
Base frequencies:
freqA = 0.3808
freqC = 0.1784
freqG = 0.1539
freqT = 0.2869
Substitution model:
Rate matrix
R(a) [A-C] = 1.0000
R(b) [A-G] = 0.8460
R(c) [A-T] = 1.0000
R(d) [C-G] = 1.0000
R(e) [C-T] = 0.1625
R(f) [G-T] = 1.0000
Among-site rate variation
Proportion of invariable sites = 0
Variable sites (G)
Gamma distribution shape parameter = 1.6368
-----
ＡＩＣの基準のもとでは，尤度比検定の結果とは異なり，ガンマ補正した田村-根井モデル（TrN+G）が最適の塩基置換モデルであると判定されます．
この結果もまた PAUP* の最尤法オプション設定として出力されます．
-----
[!Likelihood settings from best-fit model (TrN+G) selected by AIC in Modeltest Version 3.06
]
BEGIN PAUP;
Lset Base=(0.3808 0.1784 0.1539) Nst=6 Rmat=(1.0000 0.8460 1.0000 1.0000 0.1625) Rates=gamma Shape=1.6368 Pinvar=0;
END;
-----
尤度比検定では比較されるモデルの間に複雑性に関する階層性があることを仮定しますが，ＡＩＣではそのような仮定は必要ありません．また，尤度比検定であっても，χ二乗近似ではなく，パラメトリック・ブーツストラップ法を用いて検定統計量δの帰無分布を構築するならば階層性の仮定は不要になります．
------------------------------
尤度比検定と情報量基準の温度差
------------------------------
上の例のように，尤度基準と情報量基準では得られる結論が異なることが当然あり得ます．ここで，どちらの基準がすぐれているのかを論じることはきっと無意味でしょう．分析の目的というか，合理的判断のスタンスが両者では根本的に別次元にあると思われるからです．
尤度比検定は，あくまでも現時点でのデータとモデルとの適合度（fitness）を判断基準として，より適合度の高いモデルを判定しようという意図があります．一方，情報量基準は，いまのデータに対するあてはまりの程度を重視するのではなく，データが変動した場合のモデルの平均的な適合度の大きさを評価することで，モデルの善し悪しを判断します．いまあるデータにたとえ正確に適合しているモデル（尤度基準ではいい評価を得る）であっても，それがデータのバラツキに対してまずい予測的ふるまい（over-fittingみたいな性質）が見られるとすると，情報量基準のもとでは相対的に低い評価を受けることになるでしょう．
したがって，尤度比検定と情報量基準のちがいは，より深く考えるならば，分子系統推定における背景モデルの設定論拠を１）モデルの現時点でのデータとの適合性に置くのか，それとも２）モデルの潜在的予測能力に置くのか，という問題とからんでくるでしょう．両者は必ずしも矛盾なく一致した結論を導くとはかぎりません．
赤池の情報量基準は「最節約基準」の発現であるとみなす見解があります：
Forster, M. and E. Sober 1994. How to tell when simpler, more unified, or less ad hoc theories will provide more accurate predictions. The British Journal for the Philosophy of Science, 45: 1-35.
Sober, E. 1996. Parsimony and predictive equivalence. Erkenntnis, 44: 167-197.
Sober, E. 2001. What is the problem of simplicity? Pp.13-32 in: A. Zeller, H. Keuzenkamp and M. McAleer (eds.), Simplicity, Inference, and Modelling: Keeping It Sophistically Simple. Cambridge University Press, Cambridge.
この立場に立つならば，系統推定のための背景仮定やモデル導入を最少にとどめようという姿勢（三中・鈴木 2002 で論議）はある支持を得られるのかもしれません．
三中信宏・鈴木邦雄 2002. 生物体系学におけるポパー哲学の比較受容. 所収：日本ポパー哲学研究会（編），『批判的合理主義・第２巻：応用的諸問題』, pp.71-124. 未來社，東京.
今回の講義で用いた ModelTest は，最尤系統推定や分子進化研究に必要となる塩基置換モデルを選択するためのツールとして普及しているようですが，同時に系統樹構築のための背景モデルの意味を考える上でもまたとない機会をユーザーに提供しているといえるでしょう．
--- Additional Note ---
PAUP*4.0beta10とModelTest3.06を併用する場合の注意として，PAUP*の最終テスト版（beta10）は，最尤法スコアのファイル出力にバグがあることが指摘されています：
http://paup.csit.fsu.edu/problems.html#Anchor-General-28928
具体的には，最尤法によって計算した rates をファイル出力しないため，ModelTest 3.06 のように，その出力結果を要求するソフトと併用すると「入力データが不完全である」というエラーが生じるというバグです．
上記 PAUP* サイトに対処が書かれているので，下記のModelTestのPAUPブロックにそれを組み込みました．塩基配列データファイル（たとえば model3.nex）を読みこんだ後で，オープンする ModelTest のコマンドファイル（たとえば model3block.nex）の下記部分に PAUP* のスコア出力を制御する１行を挿入すればＯＫです：
--- model3block-debugged.nex ---
#NEXUS
[! ***** MODELFIT BLOCK -- MODELTEST 3.0 *****]
[The following command will calculate a NJ tree using the
JC69 model of evolution]
BEGIN PAUP;
log file= modelfit.log replace;
DSet distance=JC objective=ME base=equal rates=equal pinv=0
subst=all negbrlen=setzero;
NJ showtree=no breakties=random;
End;
[!
***** BEGIN TESTING 56 MODELS OF EVOLUTION ***** ]
BEGIN PAUP;
Set criterion=like;
default lscores longfmt=yes; ←※この１行を挿入
[!
** Model 1 of 56 * Calculating JC **]
lscores 1/ nst=1 base=equal rates=equal pinv=0
scorefile=model.scores replace;
《以下略（model3block.nexと同一）》
これによって，正しくスコアがファイル出力され，ModelTest 3.06 による分子進化モデルの尤度比検定とＡＩＣ計算ができます．
生態学会第50回大会（つくば）でポスター発表を見ていたら，上記バグのために PAUP* による最尤系統樹の計算を断念して Tree PUZZLE によろめきましたというハナシを耳にしましたので，上記のバグレポートと対処法をメモしておきました．潜在的被害者のみなさんがこれで救われればいいのですが．