はじめに
半年前ぐらいにselective inferenceについてまとめました。
その時、クラスタリングにおけるselective inferenceについて触れようと思いましたが、ブートストラップ法を用いたクラスタリングやマルチスケールブートストラップについての知識がなく、挫折しました。
ブートストラップ法を用いたクラスタリングにおける仮説検定は、生物学の分野ではよく用いられています。
特に最近では、新型コロナウィルスとSARSの遺伝子の違いを分析するのに用いられていたようです。
そこで今回ブートストラップ法を用いたクラスタリングについて調べたのでまとめます。
基本的には、以下の下平先生の解説を参考に(写経)しています。
以下の解説を読んでもらう方が早いかもしれません。
クタスタリングにおける問題
クラスタリングは、各サンプルの類似度からクラスターを作成する手法です。
形成されるクラスターは、得られているデータ内の関係性で変化します。
仮に母集団からデータを何回もサンプリングできるとすると,それぞれのサンプルでのクラスタ分析した結果は、ばらついてしまう可能性があります。
そのため、分析の結果によって得られたクラスターは、偶然同じクラスタであると判断されただけかもしれません。
しかし、ほとんどの場合は一回のクラスタリングにより得られた結果を用いて、生成したクラスターの有用性を検討している場合がほとんどです。
ブートストラップ法
ブートストラップ法は、得られているデータセットからリサンプリングを複数回行うことのより、統計量の推定値とそのばらつきの情報を得る、計算機統計学における方法です。
母集団の平均値を推定することを考えます。
はじめに標本集団から、重複を許したN個(元のサンプルと同じ数)の標本を無作為に抽出します。
これをブートストラップ標本と呼びます。
ブートストラップ標本毎に平均値を計算し、ブーストストラップ標本平均を計算します。
この操作をB回繰り返し、平均値をB個作成することで、統計量の分布を得ることができるようになります。
そして、平均値を求めることにより、母集団の平均値の推定値が得られます。
この時の分布は、繰り返し回数を大きくすると中心極限定理により正規分布に収束します。
これは、一般にはできない母集団からのサンプリングを擬似的に行うことを意味しています。
ブートストラップ法を用いた分析には統計量の推定や信頼区間の推定、検定、回帰分析がありますが、ここら辺を見てもらえればと思います。
クラスタリングにおけるブートストラップ検定
5つの対象を階層クラスタリングする場合を検討します。
それぞれの対象は個の変数により構成されています。
5つの対象により形成される可能性のある樹形図(デンドログラム)は、以下の通り15通りあります。
この15通りの樹形図の中のどれかが、真の樹形図です。
それぞれの樹形図とから構成される樹形図の類似度を計算することで、真の樹形図を推定することができます。
この時の類似度には、樹形図に対するデータの尤度を用います。
具体的には候補となるすべての樹状図 に対して対数尤度 を計算し,その中で対数尤度を最大となる樹形図を最もデータに当てはまりが良い、推定された真の樹形図であると言えます。
( 樹状図における最尤法については、こちらで解説されています「バイオインフォマティクス第6回 藤 博幸」)
しかし、得られているデータは母集団からサンプリングされた確率的な値であるため、揺らぎがあります。
そのデータから判断された樹形図が真の樹形図であるかにも、揺らぎが生じてしまいます。
そのため、ブートストラップ法によりが得られる確率の分布が推定します。
5つの対象の観測値のB回のリサンプリングにより得られる、ブートストラップ標本とします。
仮説の支持または不支持を表す関数 を用います。 上記の例では、観測値の樹形図とブートストラップ標本の樹形図が等しいかを返す関数となります。
これらのうち値が 1 になった回数Cを算出します。
ここからブートストラップ確率は次のように算出できます。
この値が有意水準を満たしていれば、得られている樹形図は有意に真の構造であると言えます。
しかし、この時に推定されるブートストラップ確率にはバイアスがあることが知られている。そのため、改良した次のマルチスケールブートストラップ法を適用する。
マルチスケールブートストラップ法
通常のブートストラップ法では、1回のリンサンプリングで回サンプリングを行います。
リサンプリングの数をからにすると、リサンプリングにより得られた統計量のばらつき(標準偏差)が変化します。
とすれば標準偏差の大きさは 倍になります。この時をリサンプリングのスケールと呼びます。
の変化に合わせてブートストラップ確率も変化します。実際には、次のようになります。
ここで、は標準正規分布関数,は符号付距離,は境界の曲率に関係した量です。
マルチスケールブートストラップ法による検定、AU検定のアルゴリズムは次のようになります。
(1) K組のそれぞれのブートストラップ標本に対応するスケールを考える。
(2) 各 について, 個のブートストラップ標本を に従って取得する.この時のぞれぞれの標本を次のように表します。
ブートストラップ標本が仮説を支持するかしないかを、次式を用いて調べる。
仮説が支持された回数は次式により計算し、ブートストラップ確率を算出する。
(3) 計算された をその理論値 の曲線に当てはめ,回帰係数とを推定する.具体的には、重みつき最小二乗法(WLS)を使って次式を最小にするようなとを計算する.
ここで、は標準正規分布関数の逆関数,は分散であり、次のように与えられる。
ここで、は標準正規密度関数である。
この時、WLSではなくが母数の二項分布に従うことを利用した最尤法でとを推定しても良い。
(4) 推定した とを使い,補正した確率を計算する。
この通常のブートストラップ法を用いた確率ではなく、修正した確率に基づくバイアスの少ない検定は、AU test(Approximately Unbiased test)と呼ばれています
マルチスケールブートストラップ法による推定量の普遍性
マルチスケールブートストラップ 法による推定量は3次の精度があります。一般に次の精度の確率値とは,サンプルサイズの増加にしたがい、検定のバイアスが [ tex: n ^{ \frac{-k}{2} } ]に比例して小さくなることを意味します.
ここで、の適当な関数があるとします。例えば、を用いて統計量を求める関数です。 ただし、は次元のベクトルであり、未知の平均ベクトル、分散が単位行列となる次のような次元多変量正規分に従っているとする。
下の図(a)のようなm次元空間の中に、ある仮説を意味する領域を考えます。 なら支持を指示して、なら仮説を支持しません。
また、領域の滑らかな境界をして、境界上でへ最も近い点をと書きます。
確率値が領域の検定に関して不変であるとは、任意の有意水準に対して次の関係が成り立つことを意味しています。
これは、不偏な検定では未知パラメータがちょうど仮説境界上にあるとき、仮説を棄却する確率がになることを意味します。
そして、がの外側にあり、離れていくほど棄却確率はよりも大きくなり、逆にがの内側に入っていくほど棄却確率はより小さくります。
この時の、補正した三次の精度を持つ確率値は次のように定義されます。
ここで、と定義され、符号はのとき負、のとき正のします。 はの曲率に関係した値で、実際にはである。
ここで、は の曲率を表す 行列の固有値で から、次のように計算できます。
しかし、とを解析的に求めるのは困難です。 そのため、マルチスケールブートストラップ法を用いてとを数値的に求めます。
ここでブートストラップ標本の数が十分に大きかった場合には、次のようにかけます。
との違いは、の符号だけです
とが一致するときは、の時であり仮説境界が平坦であると言えます。
逆に言えば、平面が曲がってくるととは、乖離してしまい推定の精度が悪くなってしまいます。
マルチスケールブートストラップ法では、とすることによりも変化します。
具体的には、Yのリサンプリングの標準偏差が倍になってしまう。 そのため、を算出する際にをとすることで、影響がなります。
これは、m次元空間の図を倍することを意味します。
これに伴い、はになり、はとなります。
つまり、空間全体が拡大されたことにより仮説領域も拡大され、境界面の曲率が小さくなります。
よって、の状態に近づけることができ、の精度の良い近似値が得られるようになる!
この曲率を小さくするために空間全体のスケールを変化させるという考え方は、最高にクールだと思いました。
全体的にさらっとした内容になりましたが、詳しい理論についてはこちらを参考にしてください。 私もあとで読みます。
Rで実施に分析してみる
マルチスケールブートストラップ法による階層クラスタリングにおける検定は、pvclustパッケージで実行ができます。
用いるデータはpvclustパッケージに含まれるlungデータです。
このデータは肺腫瘍のDNAマイクロアレイデータであり、67肺腫瘍を含む73肺組織ごとに916の遺伝子が観察されたDNAマイクロアレイデータです。
library(tidyverse) library(pvclust) library(parallel) #並列化のおまじない cl <- makeCluster(detectCores()) data(lung)
AU検定の実行はpvclust関数を用います。 引数のnbootはリサンプリングを何回行うかを指定できます。
また、method.distは距離関数 dist の引数(相関係数のような類似度関数)、method.hclustは関数 hclust の引数 method に対応するクラスタの形成方法を指定できます。
今回はデフォルトの設定で実行します。
N = 2000 sa <- 9^seq(-1,1,length=13) # クラスター分析とp値の計算 res_pvclust <- pvclust(lung, r=1/sa, nboot=N, parallel=cl)
手元のMBP(15-inch, 2018)で5分程度で終わりました。 結果を可視化してみます。
# 推定された樹形図を描く plot(res_pvclust, cex=0.5, cex.pv=0.5, offset=c(0.6,0.1,0.1,0.1)) # alphaが0.95以上のクラスターを線で囲む pvrect(res_pvclust, alpha = 0.95, pv = "au")
各分岐点のに記載されている数字は、赤はMultiscale Boostrap法によるAU testの結果、緑はBootstrap法による検定結果です。
この値が0.95以上だと分岐するクラスター同士に有意な差がない=同一のクラスタに属すると考えられます。
左から5番目の赤で囲われている二つの肺腫瘍のDNAは、BP確率では93%と有意ではありませんでしたが、AU testを用いることで、有意に同一のクラスタに属すると判別されています。
クラスタリングにおけるselective inference
ここで行っているクラスタリングにおける検定では、観測値に対してクラスタを形成し、そのクラスタに差があるかどうか検定を行っています。
この時のクラスタは、元々類似度が高いor低いで形成されるため、クラスタ間の検定を行った際に有意になりやすくなります。
これを仮説選択バイアスと呼ばれています。
このバイアスを除いた推定を行うのに必要なのがselective inferenceです。
で、最初の話題にやっと戻ってきました。
ここ話は、まだ理解できていないので、またの機会に...