統計数理（2002）
第 50 巻 第 1 号 33–44
c 2002 統計数理研究所


特集「DNA 配列の統計解析」［研究詳解］

ブートストラップ法によるクラスタ分析のバラツキ評価†下平 英寿（受付 2001 年 10 月 1 日；改訂 2002 年 1 月 25 日）

要

旨クラスタリングにおけるバラツキを確率値（p-value）として定量的に評価する方法を解説する．もし仮に母集団からデータを何回もサンプルできるとすると，それをクラスタ分析した結果は観測値毎に異なる可能性がある．つまりクラスタ分析の結果得られる樹状図やそれから導かれる群（クラスタ）はデータや特徴量のサンプリングによるバラツキの影響を受けている．そこで観測値から得られた結果がどれほど信頼できるのかを 0 から 1 の範囲の実数を値にとる確率値として表現する．これはクラスタ分析という手法の性能評価をしているのではなく，データが本来持っている情報の不確実性を定量的に評価している．この方法はデータが仮説を支持するかしないかを示す二値関数とブートストラップ法によるリサンプリングだけを使っている
ので，クラスタ分析に限らずかなり広いクラスの問題に適用可能である．仮説を表す母数空間の領域の近似的に不偏な検定から確率値は計算される．基礎となっているのは Efron（1985）と Efron and Tibshirani （1998）による符号付距離と曲率の理論である．これを実用的な手法にするためのアイデアが Shimodaira（2000, 2002）のマルチスケールブートストラップ法である．生物の DNA から進化を推定する分子系統樹の問題を例題として取り上げる．
キーワード：クラスタ分析，ブートストラップ法，マルチスケールブートストラップ法，近似的に不偏な検定，分子系統樹．

1.

はじめに

クラスタ分析（例えば竹内（1989）, p. 381）では分類対象の類似性の情報を用いてその個体をいくつかの階層的な群（クラスタ）に分けることが行われる．すなわち互いに似たものは同じ群に含まれ，さらにこれらの群を幾つかの群に分けるということを階層的に行い，分類の結果
は樹状図を用いて表される．応用研究ではしばしばクラスタ分析の結果得られたひとつの樹状図が示され，これをもとに分類対象間の関係が議論される．ところがもし仮に母集団からデー
タを何回もサンプルできるとすると，それをクラスタ分析した結果は観測値毎に異なる可能性がある．すなわちわれわれが観測したデータから得られた樹状図は，データのサンプリングに関するバラツキや分類に利用した特徴量の選択に関するバラツキに影響されている．本稿では†

統計数理研究所（現 東京工業大学 情報理工学研究科：〒152–8552 東京都目黒区大岡山 2–12–1; shimo@

is.titech.ac.jp）

34

統計数理 第 50 巻

第1号

2002

このバラツキの影響を定量的に評価し，得られた分類結果がどれほど信頼出来るのかを議論する．これはクラスタ分析という手法の性能評価をしているのではなく，データが本来持っている情報の不確実性を定量的に評価する試みである．
分類対象の個数を M とし，それぞれの個体が長さ N からなるデータを考える．これは M ×Nの大きさのデータ行列
X = (Xit ; i = 1, . . . , M, t = 1, . . . , N )として表される．Xit は個体 i の t 番目のデータを表す．一般的なクラスタ分析では，Xit は t
番目の特徴量であり，個体 i と個体 j の間の類似度を例えば

X
N(Xit − Xjt )2

t=1

によって定める．個体間の類似度から階層的なクラスタを構成する方式（計算機ソフトウエア，
アルゴリズム）には様々なものが提案されているが，いずれにしても樹状図をひとつ出力する．
どのような類似度を用いどのような方式を使うかは個々の応用では大変重要な問題であるがこ
こではそれには立ち入らず，問題に応じた適切な類似度と方式を採用していると仮定する．
あらかじめ候補となる樹状図があり，それが「本来」の樹状図であるかどうかを検定したいという場合を考える．その樹状図を T で表す．もし複数の候補 T1 , T2 , . . . がある場合には，以下の議論をそれぞれの候補で繰り返し行う．候補 T は「仮説」を表している．一方，データ Xのクラスタ分析の結果得られた樹状図を T (X) で表す．バラツキの影響で T (X) = T の場合もあるだろうしそうでない場合もある．もしたまたま T (X) = T であったなら仮説 T がもっともらしいと判断し，T (X) = T ならば仮説 T が疑わしいと考えるのがごく自然である．本稿ではさらに進めて，仮説 T が真実であるかどうかを定量的に評価し，0 から 1 の範囲の値をとる確率値（p-value）を計算する方法を述べる．
仮説は樹状図として与えられるとは限らず，群として与えられるほうが一般的かもしれない．
樹状図 T に含まれる階層的な群をT = {G1 , G2 , . . . , Gg }と表す．ただし Gi ⊂ {1, . . . , M } はひとつの群であり，Gi ⊂ Gj または Gj ⊂ Gi またはGi ∩ Gj = ∅ を満たす必要がある．g は T に含まれる群の数であり，一般的な樹状図では
g = M − 1 である．仮説となるある群 G ⊂ {1, . . . , M } が，データから得られた樹状図 T (X) にたまたま含まれていて G ∈ T (X) ならば，我々は仮説 G をもっともらしいと考えるだろうし，
逆に G ∈ T (X) ならば仮説 G を疑うだろう．しかし単に G が T (X) に含まれているかどうかだけで判断するより，信頼性の程度を定量的に確率値で与えたほうがより望ましい．
仮説を樹状図 T とするか群 G とするかいずれにしても，我々はデータ X から得られた樹状
図 T (X) が仮説を支持したときにその仮説をもっともらしいと考え，支持しなかったときに疑
わしいと考える．そこでデータ X が仮説を支持するときに値 1 をとり，そうでないときに値0 をとるような関数 S(X) を考えることにする．つまり仮説が T の場合には，T (X) = T のとき S(X) = 1，T (X) = T のとき S(X) = 0 である．仮説が G の場合には，G ∈ T (X) のとき
S(X) = 1，G ∈ T (X) のとき S(X) = 0 である．仮説が具体的にどのような形式をとるにせよ，
この関数 S(X) のみを使って我々は確率値の計算法を与える．したがってこの方法はクラスタ
分析だけに限らず，実はかなり広いクラスの問題を扱える．
確率値の計算の基礎は Efron（1979）のブートストラップ法によってデータを複製することである．具体的には 3 節で説明するが簡単に言えば，データ X の複製を乱数を用いたリサンプ
リングによって多数生成し，それがどのくらいの頻度で仮説を支持するかを数えて確率値を計

ブートストラップ法によるクラスタ分析のバラツキ評価

35

算する方法である．このような確率値の計算法は Felsenstein（1985）によって提案されて以来広く用いられている．ところがこの素朴な方法だと確率値のバイアスが一次の精度（漸近的にO(N −1/2 ) のオーダ）しかない．我々が与える計算法は確率値のバイアスが三次の精度（漸近的に O(N −3/2 ) のオーダ）である．理論的な基礎は Efron（1985）と Efron and Tibshirani（1998）
であり，仮説を支持する母数空間の領域に関する符号付距離と境界の曲率というものが関わっている．これらの理論に基づき Efron et al.（1996）では二次の精度（漸近的に O(N −1 ) のオーダ）の計算法が与えられていたが，我々の方法はそれより実装が簡単でかつ精度が良い．この近似的に不偏な検定（approximately unbiased test; AU test）の確率値を容易に計算可能にしたアイデアの中心は 4 節で述べるマルチスケールブートストラップ法である．
2.

分子系統樹の推定

クラスタ分析の具体例として，DNA から生物の進化を推定する問題を取り上げる．データ
は Shimodaira and Hasegawa（1999）で用いた図 1 にあるような 6 種の哺乳類のアミノ酸シー
ケンスである．したがってデータ行列 X = (Xit ) の各要素はアミノ酸に対応した 20 種のアルファベットを値にとる．M = 6 種の哺乳類について長さ N = 3414 のシーケンスを用いたので，行列の大きさは 6 × 3414 である．
ここでは類似度からクラスタ分析する方法ではなく，進化の確率モデルに基づいた最尤法を用いてクラスタ分析を行う（Cavalli-Sforza and Edwards（1967）, Felsenstein（1981）, 長谷川・
岸野（1996））．これを簡単に述べると以下のようになる．最尤法では各樹状図 T 毎に対数尤度
L(T, X) を計算する．対数尤度というのは，もっともらしさをあらわす量だと思えばよい．M
個の分類対象の可能な樹状図の数を n とする．候補となるすべての樹状図 T1 , T2 , . . . , Tn に対して対数尤度 L(T1 , X), . . . , L(Tn , X) を計算し，その中で対数尤度を最大にする樹状図を採用す
る．手続きは一見複雑だがデータ X から樹状図 T (X) を一つ選ぶことには変わりない．いずれ

図 1.

哺乳類（ヒト，アザラシ，ウシ，ウサギ，マウス，オポッサム）のミトコンドリア DNA のア
ミノ酸シーケンス．解析に用いた長さ N = 3414 のデータのうち t = 20 から t = 99 の部分までを示した．

(a)
図 2.

(b)

哺乳類（ヒト=1，アザラシ=2，ウシ=3，ウサギ=4，マウス=5，オポッサム=6）の系統樹．
（a）最尤法で選ばれた系統樹．
（b）最新のデータが支持している系統樹．

統計数理 第 50 巻

36

第1号

2002

にしても最尤法を自動的に行う計算機プログラム（Adachi and Hasegawa（1996）, Yang（1997）,）が開発されているので，手続きが複雑かどうかということはさほど問題では
Swoﬀord（1998）
ない．
哺乳類のデータから計算した T (X) は図 2（a）に示した．DNA から推定した生物種の樹状
図は分子系統樹とも呼ばれる．T (X) は生物が進化の過程で分化していった順序を表している．
ところが，この T (X) をそのまま真実として受け入れるのは危険である．データ X は進化の確率モデルで定義される確率変数の実現値であり，サンプリングによるバラツキがある．した
がって X から計算される T (X) にもバラツキがある．実際，その後新たに得られたデータや生
物学的な知識を動員すると，どうやら図 （2 b）が真実ではないかと現在では考えられている（長谷川政美（私信）, Cao et al.（2000）, Madsen et al.（2001）, Murphy et al.（2001））．この最新デー
タを入手する以前に戻って考えると，図 1 のデータから推定された樹状図をそのまま信じていたら誤った結論に導かれていた可能性が高い．このような早まった結論を避けるためには，バ
ラツキを考慮して T (X) の信頼性を評価することが必要になる．
3.

ブートストラップ法

データ X から計算する樹状図 T (X) のバラツキを見るには，母集団からデータを何回もサンプルして樹状図がどのように分布するかを見ればよい．ところが実際には母集団から得られるのはひとつのデータ X だけである．そこで Efron（1979）はデータからのリサンプリングによって X の複製を何回でも生成できる一般的な方法を考え，ブートストラップ法と名づけた．
これを以下で説明する．
データ行列 X = (Xit ) を
X = (x1 , x2 , . . . , xN )と書く．ただし，xt = (Xit ; i = 1, . . . , M ) はデータ行列の t 列目をあらわす．X の複製 X ∗ は
X ∗ = (xt1 , xt2 , . . . , xtN )と書かれる．ここで t1 , . . . , tN は 1, . . . , N のどれかの値を重複を許してランダムに取ることにするので，X ∗ の 1, 2, . . . 列目は X の t1 , t2 , . . . 列目を取り出したものである（図 3）．
t1 , . . . , tN は次のように乱数を使って作られる．まず 1, . . . , N からランダムに数を選びそれを t1 とする．同様に 1, . . . , N からランダムに数を選びそれを t2 とする．ここで重複を許すので t1 と t2 がたまたま同じ数になることもある．これを N 回繰り返して t1 , . . . , tN を生成する．
X x x
1 2x3 x4 x5x6
X 1*

x3 x1 x2 x2 x6 x3

10000 copiesx 3 x 5 x 1 x 2 x 1x 4

S

output

S

count
S

*
X 10000

図 3.

ブートストラップ法．ここでは N = 6 のデータから複製を B = 10000 個生成している．複製が仮説を何回支持したかがカウントされる．

ブートストラップ法によるクラスタ分析のバラツキ評価
表 1.

上位 15 個の樹状図と確率値．

表 2.

37

上位９個の群と確率値．

つまり，N 通りの目があるサイコロを N 回振ってその出た目を記録するのと同じである．そして出た目の X の列を順に取り出して X ∗ が作られる．
ブートストラップ法では X の複製 X ∗ を作る手続きを B 回繰り返し，B 個の複製∗
X1∗ , X2∗ , . . . , XBを生成する．ただし B は十分に大きな数（例えば B = 10000）とする．この多数の複製のバラ
ツキは，母集団における X のバラツキを近似的に表していると考えられる．従って∗T (X1∗ ), T (X2∗ ), . . . , T (XB)のバラツキを調べることによって，T (X) がどれほど信頼できるかが評価できる．
ここで 1 節で説明したように問題を少し一般化して仮説の支持または不支持を表す関数 S(X)を用いる．そして∗S(X1∗ ), S(X2∗ ), . . . , S(XB)

のうち値が 1 になった回数を C とする．つまり，∗
C = S(X1∗ ) + · · · + S(XB)

と書いても良い．Felsenstein（1985）はブートストラップ確率をCp̃ =Bと定義し，これが 1 に近いほど仮説はもっともらしく，0 に近いほど仮説は疑わしいと考えた．
哺乳類のデータにこのブートストラップ確率を計算した結果を表 1 と表 2 に示した．6 種の哺乳類のラベルは図 2 で示したものである．表 1 では 105 通りの樹状図をそれぞれ仮説として確率値を計算し，表 2 では 25 通りの群をそれぞれ仮説としている．ここではオポッサムは常に一番外側に置いて群 {1,2,3,4,5} が常に正しいと仮定して分析してある．少しテクニカルな話になるがオポッサムはこの場合アウトグループと呼ばれ，reversible なマルコフ過程を進化のモデルに採用した最尤法では「無根系統樹」しか推定できないという制約上必要な処置である．

38

統計数理 第 50 巻

第1号

2002

図（
2 a）の樹状図は表 1 の 1 行目に対応し，図 （2 b）は 5 行目に対応する．ブートストラップ
確率 p̃ を見ると，1 行目では p̃ ≥ 0.05 であり有意水準 α = 0.05 では棄却されない．ところが 5行目では p̃ < 0.05 であり，この樹状図は棄却される．2 節で述べたように 5 行目の樹状図が最新のデータで支持されており，解析に用いた古いデータでは残念ながらそれと矛盾する結論を導いてしまったことになる．これは表 2 にも反映されていて，8 行目の群 {1,4,5} は p̃ < 0.05
で棄却されている．この群は図 2（b）に含まれるので，やはり最新データと矛盾した結論になる．すなわちバラツキを評価していても，必ずしも正しい結論に導かれるわけではない．例えば Graur et al.（1996）では同じような哺乳類のデータを分析して，ウサギとマウスからなる群{4,5} を有意に棄却したが，これも後になってみるとおかしいと考えられている．どのような方式でバラツキを評価してもこのような事態は起こりえるが，問題なのは素朴なブートストラッ
プ確率にバイアスがあり，誤った結論に導かれる可能性が必要以上に高いということである．
そこで次の節で述べるような改良が必要になる．
4.

マルチスケールブートストラップ法

ブートストラップ法では X からランダムに N 個の列を取り出して複製 X ∗ を作った．もし取り出す個数（つまり複製の長さ）を変えて N  とすると複製はX ∗ = (xt1 , xt2 , . . . , xtN  )となる．普通のブートストラップ法では N  = N であるが，もし N  = N とすると複製のバラ√ツキの程度（標準偏差）が変化する．例えば N  = 2N とすればバラツキの程度は 1/ 2 倍にな√り，逆に N  = N/2 とすればバラツキの程度は 2 倍になる．一般に N  = rN とすればバラツ√√キの程度は 1/ r 倍になる．1/ r を複製のスケールと呼び，それはデータ長の比 r によって制御できる．ブートストラップ確率 p̃ も N  ，つまり r に依存して変化する．実際 5 節で説明するように，N  = rN の時のブートストラップ確率の「理論値」は
√
√
(4.1)
π(r; d, c) = 1 − Φ(d r + c / r)で与えられる．ただし Φ(·) は標準正規分布関数，d は符号付距離，c は境界の曲率に関係した量であり，詳細は後ほど説明される．Shimodaira（2000, 2002）は p̃ の変化からより精度の高い AU test の確率値を求める方法を提案し，これをマルチスケールブートストラップ法と名づけた．この手続きは以下のようになる（図 4）．
ステップ 1．K 組のブートストラップを考える．データ長の比 r1 , r2 , . . . , rK ，複製の個数
B1 , B2 , . . . , BK を定める．以下の数値例では，K = 10, r1 = 0.5, r2 = 0.6, . . . , r10 = 1.4, B1 =· · · = B10 = 10, 000 を用いた．
ステップ 2．各 k = 1, . . . , K について，Bk 個の複製を N  = rk N を使って生成する．これを∗X1∗ (rk ), X2∗ (rk ), . . . , XB
(rk )kと書く．そして複製が仮説を支持するかしないかを∗S(X1∗ (rk )), S(X2∗ (rk )), . . . , S(XB(rk ))
k

によって調べる．仮説が支持された回数は∗C(rk ) = S(X1∗ (rk )) + S(X2∗ (rk )) + · · · + S(XB(rk ))kブートストラップ法によるクラスタ分析のバラツキ評価
X x x
1 2x3 x4 x5x6

S

39

outputX 1*
x5 x2 x1x3

x 1x 6 x 4 x 5 x 3x 5x 6x 6S

10000 copies

S

10000 copiesx6 x5x4 x1

Sx 2 x 5 x 4 x 2 x 6x 6 x 5 x 3

S*
X 10000
count

count

図 4.

マルチスケールブートストラップ法．ここでは N = 6 のデータから， N  = 4 と N  = 8 の複製を B1 = B2 = 10000 個ずつ生成している．各々のブートストラップ法で，複製が仮説を何回支持したかがカウントされる．

である．これからブートストラップ確率をp̃(rk ) =C(rk )
Bkと計算する．
ステップ 3．計算された p̃(rk ) をその理論値 π(rk ; d, c) の曲線に当てはめ，回帰係数 d と c を推定する．具体的には重みつき最小二乗法（WLS）を使って

X
K

RSS(d, c) =(Φ−1 (π(rk ; d, c)) − Φ−1 (p̃(rk )))2 /vk ,k=1を最小にするような d と c を計算する．ここで Φ−1 (·) は Φ(·) の逆関数であり，分散 vk は
vk = p̃(rk )(1 − p̃(rk ))/(φ(Φ−1 (p̃(rk )))2 Bk )で与えられる．φ(·) は標準正規密度関数である．
ステップ 4．推定した d と c を使い，補正した確率値をp̂ = 1 − Φ(d − c)で計算する．
ただしステップ 3 における WLS をやめて，Bk p̃(rk ) が母数 π(rk ; d, c) の二項分布に従うことを利用した最尤法（MLE）で d と c を推定してもよい．これには WLS で推定した d と c を初期値として，ニュートン法などで数値的にX
K

L(d, c) =Bk {p̃(rk ) log π(rk ; d, c) + (1 − p̃(rk )) log(1 − π(rk ; d, c))} ,k=1を最大化して d と c を計算する．マルチスケールブートストラップ法は WLS と MLE の両方とも計算機ソフトウエア CONSEL（Shimodaira and Hasegawa（2001））に実装されている．
CONSEL は分子系統樹のソフトウエアとの連携を意識して作られているが，その他の問題にも使える．ユーザは p̃(rk ) を自分の問題にあわせて計算すればよい．
補正した確率値 p̂ を哺乳類のデータで計算した結果は表 1 と表 2 に示されている．同じデー
タセットを用いた同様の計算は Shimodaira（2002）で行った．全般的に p̂ は p̃ より大きめの値

統計数理 第 50 巻

40

第1号

2002

になり各仮説は棄却されにくくなる．これは c ≥ 0 となっていることから理解できるのだが，じつは次節で述べる仮説の領域 H が凸であり曲率が正であることが関係している．表 1 の 5
行目を見ると p̂ ≥ 0.05 となって，この樹状図はもう棄却されなくなる．結果として最新データと矛盾しない結論を導いている．このことは表 2 の 8 行目にも反映されており，群 {1,4,5} はもう棄却されない．データのバラツキによって誤った仮説（表 1 の 1 行目，表 2 の 1, 2 行目)が最も高く支持されていたが，最新データによって最も高く支持されるようになった仮説（表1 の 5 行目，表 2 の 7, 8 行目）も否定されていなかったわけである．
5.

近似的に不偏な検定

このようにして補正した確率値 p̂ は，素朴なブートストラップ確率 p̃ より AU test としては一般的にずっと精度が良く，すなわち検定のバイアスが小さい．このことを以下で説明する．
まず X の適当な関数
Y = f (X)を考える．ただし Y はベクトルでその次元を m とする．そしてY ∼ Nm (µ, Im )(5.1)

のように未知の平均ベクトル µ，共分散が単位行列の m 次元多変量正規分布に従っていると仮定する．逆にいうと，少なくとも近似的に（5.1）式が適当な m で成り立つような関数 f (X) の存在を仮定する．そして m 次元空間の領域 H を考え，f (X) ∈ H なら S(X) = 1，f (X) ∈ H
なら S(X) = 0 とする．つまり，Y ∈ H ならデータは仮説を支持し，Y ∈ H なら支持しない（図 5）．領域 H の境界を ∂H と書き，境界上で Y へ最も近い点を µ̂ と書くことにする．そして境界 ∂H は滑らかであると仮定する．
さて確率値 p が領域 H の検定に関して不偏であるとは，任意の有意水準 0 < α < 1 に対して
Pr{p < α | µ} ≤ α ,µ∈H

Pr{p < α | µ} ≥ α ,µ ∈ H

Pr{p < α | µ} = α ,µ ∈ ∂Hが成り立つことを言う．従って(5.2)

が成り立つ．一般に p < α のとき仮説は棄却される．不偏な検定では未知パラメタ µ がちょう
ど仮説の境界上 ∂H にあるとき，仮説を棄却する確率が α になる．そして µ が H の外側に出て

(a)
図 5.

(b)

データベクトル Y が仮説を支持する領域 H; Shimodaira (2002)．領域の境界 ∂H 上の点でY に最も近い点 µ̂ から Y までの符号付距離が d．
（a）境界 ∂H が滑らか．
（b）境界 ∂H が尖っている．

ブートストラップ法によるクラスタ分析のバラツキ評価

41

離れていくほど棄却確率は α より大きくなり，逆に µ が H の内側に入っていくほど棄却確率は α より小さくなる．結論から言うと，ブートストラップ確率 p̃ は（5.2）の誤差が O(N −1/2 )であるが，補正した確率値 p̂ は（5.2）の誤差が O(N −3/2 ) になる．適当に大きな数を N に代入するとわかるが，N −1/2 より N −3/2 のほうが小さな値になる．つまり p̃ より p̂ のほうが誤差が小さい．
Efron（1985）の「補題」もしくは Efron and Tibshirani（1998）の（2.16）式によれば，三次の精度をもつ補正した確率値はp̂ = 1 − Φ(d − c)(5.3)

と書ける．ただし，d = ±Y − µ̂ で定義し，符号は Y ∈ H のとき負，Y ∈ H のとき正とす
る．c は ∂H の曲率に関係した量である．実際 c = c1 − dc2 と書けて，c1 = λ1 + · · · + λM −1と c2 = λ21 + · · · + λ2M −1 は µ̂ における ∂H の曲率を表す (M − 1) × (M − 1) 行列の固有値λ1 , . . . , λM −1 から計算できる．
たしかに式（5.3）は精度の高い確率値を与えているが，実際の応用では d や c を解析的に与えることは非常に困難であるから，このままでは（5.3）は役に立たない．そこで 4 節のマル
チスケールブートストラップ法が開発され，d と c を現実の問題で数値的に計算することが可能になった．その仕組みを理解するために，ブートストラップ確率が十分大きな B でp̃ = 1 − Φ(d + c)(5.4)

と書けるという Efron and Tibshirani（1998）の（2.19）式を利用する．
（5.3）と（5.4）の違い
は c の符号だけである．したがってもし境界 ∂H が平坦で c = 0 ならば p̃ = p̂ となる．ところが ∂H が曲がってくると p̃ と p̂ は逆の方向へ変化してしまうので，p̃ は p̂ の推定値としては精度が悪くなる．
マルチスケールブートストラップでデータ長を N  = rN とすると（5.4）も変化する．Y の√複製のバラツキが 1/ r 倍となってしまうので，これを元に戻して（5.4）を利用するには，Y√√の代わりに r Y を考えればよい．つまり図 5 の絵全体を r 倍拡大するのと同じである．こ√√うすると同時に d が r d になり c が c/ r になってしまう．こうして（4.1）で与えたブート
ストラップ確率の理論値 π(r; d, c) が出てくる．これで 4 節の方法が理解できたことになる．
仮定したモデル（5.1）は制約が強すぎるように見えるかもしれない．例えば共分散行列が単位行列というのはとても強い制約である．ところが共分散行列が一般の場合でもそれを単位行列にするような Y の線形変換が存在する．従ってその線形変換も f (·) に含めてしまえば結局（5.1）に帰着できる．任意の滑らかな非線形変換を f (·) として使えるので，かなり広い範囲の問題が（5.1）に帰着できる．そして変換 f (·) におけるブートストラップ確率の不変性より，4
節の方法はそのまま使えることになる．
しかしまったく問題が無いわけではない．領域の境界はしばしば滑らかでなく，図 5（b）のように尖っている．滑らかな変換 f (·) ではこの特異性を消すことができず，結局図 5（b）を図5（a）で近似することになる．これが検定のバイアスにつながり，特異性が無視できない場合には補正した確率値は必ずしも高い精度をもつわけではない．
6.

おわりに

データが仮説を支持する (S(X) = 1) かしない (S(X) = 0) かという情報と，ブートストラッ
プ法によるリサンプリングだけを使って AU test として精度の良い確率値を計算する方法を解説した．問題設定のシンプルさから，この方法はクラスタ分析に限らず，かなり広いクラスの問

42

統計数理 第 50 巻

第1号

2002

題に適用可能であろう．基礎となっているのは Efron（1985）と Efron and Tibshirani（1998）による符号付距離と曲率の理論である．これを実用的な手法にするためのアイデアが Shimodaira（2000, 2002）のマルチスケールブートストラップ法である．
マルチスケールブートストラップ法の構成要素として本稿では単純なブートストラップ法によるリサンプリングを用いたが，応用では問題ごとにサンプリング法の工夫が必要である．例
えば時系列データではブロックブートストラップ法を用いる．クラスタリングの特徴量を直接
サンプリングするような場合は，特徴量の影響度に応じた重み付けしたリサンプリングが意味を持つだろう．いずれにしても形式的にブートストラップ法を適用するのではなく，何に関するバラツキを評価したいのかを適切に考慮する必要がある．
問題によっては計算量を減らす工夫が必要である．例えば分子系統樹の解析では，樹状図の対数尤度 L(Ti , Xb∗ ) をすべての i = 1, . . . , n, b = 1, . . . , B に対して計算するのは非常にコストがかかる．そこで，一種の線形近似である RELL 法（Kishino et al.（1990））を用いて近似的に対数尤度を計算した．Hasegawa and Kishino（1994）の数値例や Shimodaira（2001）の補題 1 で示したように，N が大きい問題ではこの近似の精度は十分に良い．
分子系統樹の解析では対数尤度を最大にする樹状図を選択しているので，これは統計的モデ
ル選択の一例になる．下平（1993, 1999），Shimodaira（1998）はリサンプリングを利用して対数尤度の多重比較を行いモデル選択の信頼性を確率値として評価する方法を提案した．これはShimodaira-Hasegawa（SH）test（Shimodaira and Hasegawa（1999）, Goldman et al.（2000））として分子系統樹の分野で利用されはじめている．この分野では Kishino-Hasegawa（KH）test（Kishino and Hasegawa（1989））がブートストラップ確率に並ぶ標準手法として広く利用されているが，SH test は KH test で見落とされていた「選択バイアス」を多重比較法で補正したものである．本稿の AU test における補正した確率値 p̂ は，実は SH test と定性的には同じ役割を果たしている．つまり選択バイアスの大きさと境界の曲率とは定性的には同じものである．
ただし多重比較では図 5（b）の尖りの先端に µ があると仮定することによって「最悪ケース」を想定した補正を行っているのに対して，本稿の AU test では図 5（a）のように滑らかな境界を仮定したうえで µ̂ 周辺の「典型的ケース」を想定した補正を行っている．これらの方法から計算される確率値は想定したケースの違いを反映して定量的には異なってくる．SH test のほ
うが AU test より保守的な結果になる．
計算機ソフトウエア CONSEL（Shimodaira and Hasegawa（2001））の DOS バイナリと UNIX
ソースコードは著者より無償で入手可能である．これは AU test, KH test, SH test, ブートスト
ラップ確率などを同時に計算する．

参 考 文

献

Adachi, J. and Hasegawa, M.（1996）. MOLPHY version 2.3: Programs for molecular phylogenetics
based on maximum likelihood, Comput. Sci. Monographs, No. 28, The Institute of Statistical
Mathematics, Tokyo.
Cao, Y., Fujiwara, M., Nikaido, M., Okada, N. and Hasegawa, M. （2000）. Interordinal relationships
and timescale of eutherian evolution as inferred from mitochondrial genome data, Gene, 259,
149–158.
Cavalli-Sforza, L. L. and Edwards, A. W. F.（1967）. Phylogenetic analysis: Models and estimation
procedures, Evolution, 32, 550–570.
. Bootstrap methods: Another look at the jackknife, Ann. Statist., 7, 1–26.
Efron, B.（1979）. Bootstrap conﬁdence intervals for a class of parametric problems, Biometrika, 72,
Efron, B.（1985）

ブートストラップ法によるクラスタ分析のバラツキ評価

43

45–58.
Efron, B. and Tibshirani, R.（1998）. The problem of regions, Ann. Statist., 26, 1687–1718.
. Bootstrap conﬁdence levels for phylogenetic trees,
Efron, B., Halloran, E. and Holmes, S.（1996）
Proc. Nat. Acad. Sci. U.S.A., 93, 13429–13434.
Felsenstein, J.（1981）. Evolutionary trees from DNA sequences: A maximum likelihood approach,
Journal of Molecular Evolution, 17, 368–376.
. Conﬁdence limits on phylogenies: An approach using the bootstrap, Evolution,
Felsenstein, J.（1985）
39, 783–791.
Goldman, N., Anderson, J. P. and Rodrigo, A. G.（2000）. Likelihood-based tests of topologies inphylogenetics, Systematic Biology, 49, 652–670.
. Phylogenetic position of the order Lagomorpha （rabbits,
Graur, D., Duret, L. and Gouy, M.（1996）
hares and allies）, Nature, 379, 333–335.
Hasegawa, M. and Kishino, H.（1994）. Accuracies of the simple methods for estimating the bootstrap
probability of a maximum likelihood tree, Molecular Biology and Evolution, 11, 142–145.
長谷川政美，岸野洋久（1996）. 『分子系統学』，岩波書店，東京.
Kishino, H. and Hasegawa, M.（1989）. Evaluation of the maximum likelihood estimate of the evolutionary tree topologies from DNA sequence data, and the branching order in Hominoidea,
Journal of Molecular Evolution, 29, 170–179.
. Maximum likelihood inference of protein phylogeny
Kishino, H., Miyata, T. and Hasegawa, M.（1990）
and the origin of chloroplasts, Journal of Molecular Evolution, 30, 151–160.
Madsen, O., Scally, M., Douady, C. J., Kao, D. J., DeBry, R. W., Adkins, R., Amrine, H. M., Stanhope,
M. J., de Jong, W. W. and Springer, M. S. （2001）. Parallel adaptive radiations in two major
clades of placental mammals, Nature, 409, 610–614.
Murphy, W. J., Eizirik, E., Johnson, W. E., Zhang, Y. P., Ryder, O. A. and O’Brien, S. J.（2001）.
Molecular phylogenetics and the origins of placental mammals, Nature, 409, 614–618.
下平英寿（1993）. モデルの信頼集合と地図によるモデル探索，統計数理，41, 131–147.
. An application of multiple comparison techniques to model selection, Ann.
Shimodaira, H.（1998）
Inst. Statist. Math., 50, 1–13.
下平英寿（1999）. モデル選択理論の新展開，統計数理，47, 3–27.
Shimodaira, H.（2000）. Another calculation of the p-value for the problem of regions using the scaled
bootstrap resamplings, Tech. Report, No. 2000-35, Stanford University, California.
Shimodaira, H.（2001）. Multiple comparisons of log-likelihoods and combining nonnested models withapplications to phylogenetic tree selection, Comm. Statist. Theory Methods, 30, 1751–1772.
. An approximately unbiased test of phylogenetic tree selection, Systematic
Shimodaira, H.（2002）
Biology, 51, 492–508.
Shimodaira, H. and Hasegawa, M.（1999）. Multiple comparisons of log-likelihoods with applications
to phylogenetic inference, Molecular Biology and Evolution, 16, 1114–1116.
. CONSEL: For assessing the conﬁdence of phylogenetic
Shimodaira, H. and Hasegawa, M.（2001）
tree selection, Bioinformatics, 17, 1246–1247.
Swoﬀord, D. L.（1998）. PAUP*. Phylogenetic Analysis Using Parsimony（*and Other Methods）,
Version 4, Sinauer Associates, Sunderland, Massachusetts.
竹内 啓 編（1989）. 『統計学辞典』，東洋経済新報社，東京.
. PAML: A program package for phylogenetic analysis by maximum likelihood,Yang, Z.（1997）
CABIOS, 13, 555–556.

44

Proceedings of the Institute of Statistical Mathematics Vol. 50, No. 1, 33–44 (2002)

Assessing the Uncertainty of the Cluster Analysis
Using the Bootstrap Resampling
Hidetoshi Shimodaira(The Institute of Statistical Mathematics)

This paper reviews the method of calculating the p-value for assessing the uncertainty
of cluster analysis. Considering that the dendrogram as well as the derived clusters
obtained by the cluster analysis is subject to change due to the fluctuation of the samplingof the data or that of the characters, the reliability of the result is represented as the pvalue, between 0 and 1. This method is applicable to a wide class of problems, and is
not limited to cluster analysis, since it uses only bootstrap resampling and the 0/1-valuefunction to indicate whether the data supports the hypothesis. The p-value is calculated
from the approximately unbiased test of the region in the parameter space representing
the hypothesis. The method is based on the theory of “signed distance” and “curvature”
by Efron (1985) and Efron and Tibshirani (1998). The key idea to convert the theory intoa practical algorithm is the multiscale bootstrap resampling of Shimodaira (2000, 2002).
The issue is illustrated by the phylogeny analysis to infer the history of evolution fromthe DNA sequences.

Key words: Cluster analysis, bootstrap, multiscale bootstrap, approximately unbiased test, molecularphylogeny.