はじめに

半年前ぐらいにselective inferenceについてまとめました。

qiita.com

その時、クラスタリングにおけるselective inferenceについて触れようと思いましたが、ブートストラップ法を用いたクラスタリングやマルチスケールブートストラップについての知識がなく、挫折しました。

ブートストラップ法を用いたクラスタリングにおける仮説検定は、生物学の分野ではよく用いられています。

特に最近では、新型コロナウィルスとSARSの遺伝子の違いを分析するのに用いられていたようです。

そこで今回ブートストラップ法を用いたクラスタリングについて調べたのでまとめます。

基本的には、以下の下平先生の解説を参考に(写経)しています。

以下の解説を読んでもらう方が早いかもしれません。

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

クタスタリングにおける問題

クラスタリングは、各サンプルの類似度からクラスターを作成する手法です。

形成されるクラスターは、得られているデータ内の関係性で変化します。

仮に母集団からデータを何回もサンプリングできるとすると，それぞれのサンプルでのクラスタ分析した結果は、ばらついてしまう可能性があります。

そのため、分析の結果によって得られたクラスターは、偶然同じクラスタであると判断されただけかもしれません。

しかし、ほとんどの場合は一回のクラスタリングにより得られた結果を用いて、生成したクラスターの有用性を検討している場合がほとんどです。

ブートストラップ法

ブートストラップ法は、得られているデータセットからリサンプリングを複数回行うことのより、統計量の推定値とそのばらつきの情報を得る、計算機統計学における方法です。

母集団の平均値を推定することを考えます。

はじめに標本集団 $X= \{ x_1, x_2,...,x_n \}$ から、重複を許したN個(元のサンプルと同じ数)の標本を無作為に抽出します。

これをブートストラップ標本 $X^{\ast b}$ と呼びます。

ブートストラップ標本毎に平均値を計算し、ブーストストラップ標本平均 $\bar X ^{\ast b}$ を計算します。

この操作をB回繰り返し、平均値 $\bar X ^{\ast b}$ をB個作成することで、統計量の分布を得ることができるようになります。

そして、平均値を求めることにより、母集団の平均値 $\theta$ の推定値 $\hat \theta$ が得られます。

$\displaystyle{ \hat \theta = \frac{1}{B} \sum_{b=1}^{B} \bar X ^{\ast b} }$

この時の分布は、繰り返し回数を大きくすると中心極限定理により正規分布に収束します。

これは、一般にはできない母集団からのサンプリングを擬似的に行うことを意味しています。

ブートストラップ法を用いた分析には統計量の推定や信頼区間の推定、検定、回帰分析がありますが、ここら辺を見てもらえればと思います。

qiita.com

クラスタリングにおけるブートストラップ検定

5つの対象 $X$ を階層クラスタリングする場合を検討します。

それぞれの対象は $N$ 個の変数により構成されています。

5つの対象により形成される可能性のある樹形図（デンドログラム） $T(X)$ は、以下の通り15通りあります。

f:id:saltcooky:20200308181752g:plain:w450

(出典)

この15通りの樹形図の中のどれかが、真の樹形図です。

それぞれの樹形図と $X$ から構成される樹形図の類似度を計算することで、真の樹形図を推定することができます。

この時の類似度には、樹形図に対するデータの尤度を用います。

具体的には候補となるすべての樹状図 $T_1, T_2,...,T_n$ に対して対数尤度 $L(T_1, X),...,L(T_n, X)$ を計算し，その中で対数尤度を最大となる樹形図を最もデータに当てはまりが良い、推定された真の樹形図であると言えます。

( 樹状図における最尤法については、こちらで解説されています「バイオインフォマティクス第6回藤博幸」)

しかし、得られているデータ $X$ は母集団からサンプリングされた確率的な値であるため、揺らぎがあります。

そのデータから判断された樹形図 $T_i(X)$ が真の樹形図であるかにも、揺らぎが生じてしまいます。

そのため、ブートストラップ法により $T_i(X)$ が得られる確率の分布が推定します。

5つの対象 $X$ の観測値のB回のリサンプリングにより得られる、ブートストラップ標本 $X^{\ast b}$ とします。

仮説の支持または不支持を表す関数 $S(X)$ を用います。上記の例では、観測値の樹形図とブートストラップ標本の樹形図が等しいかを返す関数となります。

$\displaystyle{ S(X_1^\ast ), S(X_2^\ast ),...,S(X_B^\ast) }$

これらのうち値が 1 になった回数Cを算出します。

$\displaystyle{ C = S(X_1^\ast ) + ··· + S(X_B^\ast) }$

ここからブートストラップ確率は次のように算出できます。

$\displaystyle{ \tilde p = \frac {C}{B} }$

この $p$ 値が有意水準を満たしていれば、得られている樹形図は有意に真の構造であると言えます。

しかし、この時に推定されるブートストラップ確率にはバイアスがあることが知られている。そのため、改良した次のマルチスケールブートストラップ法を適用する。

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

通常のブートストラップ法では、１回のリンサンプリングで $N$ 回サンプリングを行います。

リサンプリングの数を $N$ から $N'$ にすると、リサンプリングにより得られた統計量のばらつき(標準偏差)が変化します。

$N = rN$ とすれば標準偏差の大きさは $1/ \sqrt{r}$ 倍になります。この時 $1/ \sqrt{r}$ をリサンプリングのスケールと呼びます。

$r$ の変化に合わせてブートストラップ確率 $\bar p$ も変化します。実際には、次のようになります。

$\displaystyle{ \pi(r; d, c)=1 − \Phi(d \sqrt{r} + c / \sqrt{r}) }$

ここで、 $\Phi (・)$ は標準正規分布関数， $d$ は符号付距離， $c$ は境界の曲率に関係した量です。

マルチスケールブートストラップ法による検定、AU検定のアルゴリズムは次のようになります。

(1) K組のそれぞれのブートストラップ標本に対応するスケール $r_1, r_2,...,r_K$ を考える。

(2) 各 $k = 1,...,K$ について， $B_k$ 個のブートストラップ標本を $N = r_kN$ に従って取得する．この時のぞれぞれの標本を次のように表します。

$\displaystyle{ X^{\ast}_1 (rk), X^{\ast}_2 (r_k),...,X^{\ast}_B (r_k) }$

ブートストラップ標本が仮説を支持するかしないかを、次式を用いて調べる。

$\displaystyle{ S(X_1^\ast (r_k)), S(X^{\ast}_2 (r_k)),...,S(X^{\ast}_{Bk} (r_k)) }$

仮説が支持された回数は次式により計算し、ブートストラップ確率を算出する。

$\displaystyle{ C(rk) = S(X^{\ast}_1 (r_k)) + S(X^{\ast}_2 (r_k)) + ··· + S(X^{\ast}_{B_k} (r_k)) }$

$\displaystyle{ \tilde p = \frac{C(r_k)}{B_k} }$

(3) 計算された $\tilde p(r_k)$ をその理論値 $\pi (r_k; d, c)$ の曲線に当てはめ，回帰係数 $d$ と $c$ を推定する．具体的には、重みつき最小二乗法(WLS)を使って次式を最小にするような $d$ と $c$ を計算する．

$\displaystyle{ RSS(d,c) = \sum ^K _{k=1} \frac{(\Phi ^{-1}(\pi (r_k;d,c))-\Phi^{-1}(\tilde p(r_k)))^2 }{v_k} }$

ここで、 $\Phi ^{-1}(・)$ は標準正規分布関数の逆関数， $v_k$ は分散であり、次のように与えられる。

$\displaystyle{ v_k = \frac{\tilde p (r_k)(1-\tilde p (r_k))}{(\phi (\Phi^{-1}(\tilde p (r_k)))^2 B_k)} }$

ここで、 $\phi$ は標準正規密度関数である。

この時、WLSではなく $B_k \tilde p(r_k)$ が母数 $\pi(r_k;d,c)$ の二項分布に従うことを利用した最尤法で $d$ と $c$ を推定しても良い。

(4) 推定した $d$ と $c$ を使い，補正した確率を計算する。

$\displaystyle{ \hat p = 1-\Phi (d-c) }$

この通常のブートストラップ法を用いた確率ではなく、修正した確率に基づくバイアスの少ない検定は、AU test(Approximately Unbiased test)と呼ばれています

マルチスケールブートストラップ法による推定量の普遍性

マルチスケールブートストラップ法による推定量は３次の精度があります。一般に $k$ 次の精度の確率値とは，サンプルサイズ $n$ の増加にしたがい、検定のバイアスが [ tex: n ^{ \frac{-k}{2} } ]に比例して小さくなることを意味します．

ここで、 $X$ の適当な関数 $Y=f(X)$ があるとします。例えば、 $X$ を用いて統計量を求める関数です。ただし、 $Y$ は $m$ 次元のベクトルであり、未知の平均ベクトル $\mu$ 、分散が単位行列となる次のような $m$ 次元多変量正規分に従っているとする。

$\displaystyle{ Y \sim N_m (\mu ,I_m)}$

下の図(a)のようなm次元空間の中に、ある仮説を意味する領域 $H$ を考えます。 $Y \in H$ なら支持を指示して、 $Y \notin H$ なら仮説を支持しません。

f:id:saltcooky:20200302230825p:plain:w450

(出典:p.40)

また、領域 $H$ の滑らかな境界を $\partial H$ して、境界上で $Y$ へ最も近い点を $\hat \mu$ と書きます。

確率値 $p$ が領域 $H$ の検定に関して不変であるとは、任意の有意水準 $\alpha$ に対して次の関係が成り立つことを意味しています。

$\displaystyle{ Pr \left \{ p < \alpha | \mu \right \} = \alpha , \mu \notin H }$

これは、不偏な検定では未知パラメータ $\mu$ がちょうど仮説境界上にあるとき、仮説を棄却する確率が $\alpha$ になることを意味します。

そして、 $\mu$ が $H$ の外側にあり、離れていくほど棄却確率は $\alpha$ よりも大きくなり、逆に $\mu$ が $H$ の内側に入っていくほど棄却確率は $\alpha$ より小さくります。

この時の、補正した三次の精度を持つ確率値は次のように定義されます。

$\displaystyle{ \hat p = 1-\Phi (d-c) }$

ここで、 $d=| Y-\hat \mu |$ と定義され、符号は $Y \in H$ のとき負、 $Y \notin H$ のとき正のします。 $c$ は $\partial H$ の曲率に関係した値で、実際には $c=c_1-dc_2$ である。

ここで、 $c_1, c_2$ は $\partial H$ の曲率を表す $(M − 1) × (M − 1)$ 行列の固有値で $\lambda_1 , \dots , \lambda _{M−1}$ から、次のように計算できます。

$\displaystyle{ c_1 = \lambda_1 + \dots + \lambda_{M−1} }$

$\displaystyle{ c_2 = \lambda ^2_1 + \dots + \lambda ^2_{M−1} }$

しかし、 $d$ と $c$ を解析的に求めるのは困難です。そのため、マルチスケールブートストラップ法を用いて $d$ と $c$ を数値的に求めます。

ここでブートストラップ標本の数 $B$ が十分に大きかった場合には、次のようにかけます。

$\displaystyle{ \tilde p = 1-\Phi (d+c) }$

$\hat p$ と $\tilde p$ の違いは、 $c$ の符号だけです

$\hat p$ と $\tilde p$ が一致するときは、 $c=0$ の時であり仮説境界 $\partial H$ が平坦であると言えます。

逆に言えば、平面 $\partial H$ が曲がってくると $\hat p$ と $\tilde p$ は、乖離してしまい推定の精度が悪くなってしまいます。

マルチスケールブートストラップ法では、 $N'=rN$ とすることにより $\tilde p$ も変化します。

具体的には、Yのリサンプリングの標準偏差が $1/\sqrt{r}$ 倍になってしまう。そのため、 $\tilde p$ を算出する際に $Y$ を $\sqrt{r}Y$ とすることで、影響がなります。

これは、m次元空間の図を $\sqrt{r}Y$ 倍することを意味します。

これに伴い、 $d$ は $\sqrt{r} d$ になり、 $c$ は $c/ \sqrt{r}$ となります。

つまり、空間全体が拡大されたことにより仮説領域も拡大され、境界面の曲率が小さくなります。

よって、 $c=0$ の状態に近づけることができ、 $\hat p$ の精度の良い近似値が得られるようになる！

この曲率を小さくするために空間全体のスケールを変化させるという考え方は、最高にクールだと思いました。

全体的にさらっとした内容になりましたが、詳しい理論についてはこちらを参考にしてください。私もあとで読みます。

マルチスケールブートストラップの漸近理論

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")

f:id:saltcooky:20200309224854p:plain

各分岐点のに記載されている数字は、赤はMultiscale Boostrap法によるAU testの結果、緑はBootstrap法による検定結果です。

この値が0.95以上だと分岐するクラスター同士に有意な差がない=同一のクラスタに属すると考えられます。

左から5番目の赤で囲われている二つの肺腫瘍のDNAは、BP確率では93%と有意ではありませんでしたが、AU testを用いることで、有意に同一のクラスタに属すると判別されています。

クラスタリングにおけるselective inference

ここで行っているクラスタリングにおける検定では、観測値に対してクラスタを形成し、そのクラスタに差があるかどうか検定を行っています。

この時のクラスタは、元々類似度が高いor低いで形成されるため、クラスタ間の検定を行った際に有意になりやすくなります。

これを仮説選択バイアスと呼ばれています。

このバイアスを除いた推定を行うのに必要なのがselective inferenceです。

で、最初の話題にやっと戻ってきました。

ここ話は、まだ理解できていないので、またの機会に...

⬜︎⬜︎⬜︎

データ分析とかの備忘録か, 趣味の話か, はたまた

マルチスケールブートストラップ法によるクラスタリングの有意差検定