はじめに
積み本消化のために、Rで学ぶデータサイエンスシリーズのブートストラップ入門を読みましたので、まとめています。
近年、統計学の界隈ではp値の使用を禁止したり、乱用を減らしていこうという流れが大きくなってきています。
ブートストラップ法は、その流れの中で使用することが推奨されるようになってきていると感じます。
今回は、本題のブートストラップ法を仮説検定についてです。
仮説検定とブートストラップ法
統計的仮説検定では、データと帰無仮説$H_0$との乖離を表現した検定統計量$T$に基づいて行われます。
$T$の値が大きくなるほど、帰無仮説$H_0$に反する根拠が強くなることを意味しています。
検定統計量$T$の観測値を$t$とすると、$H_0$に反する強さを有意確率、または$p$値は次のように表現されます。
p=Pr(T\geq t |H_0)
そして、パラメトリックモデルが仮定できるときには、帰無仮説と対立仮説は次のように、分布のパラメータに関する式で表現します。
これらの仮説に基づき検定を進めていきます。
H_0:\theta=\theta_0 \\
H_1:\theta \ne \theta_0 \\
一方でノンパラメトリックモデルでは、分布の形が不明な場合が多く、適切な検定統計量$T$の設計が難しいです。
また、データが従う分布$F_0$が決まらないような複合帰無仮説を検討する際には、上記の$p$値の定義では問題がある場合があります。
そのため、データから$H_0$を満たす分布を$\hat F_0$で推定し、$p$値を近似する必要があります。
そこで、標本のサイズ$n$が増えていく時の検定統計量の極限分布で、$p$値の近似計算をおこないます。
この計算をリサンプリング法=ブートストラップ法により行います。
ブートストラップ法による仮説検定は、通常の算出の精度では満足できない場合や近似法が存在していないときに利用できる有力な手法であるとされています。
モンテカルロ検定
帰無仮説で規定される分布から独立に抽出したい$B$組の標本から計算される$T$の値$t^{\ast 1},\cdots , t^{\ast B}$を比較します。
$H_0$のもとで$B+1$個の値$t,t^{\ast 1},\cdots , t^{\ast B}$が同様に確からしい値となり、$T^{\ast 1},\cdots , T^{\ast B}$の大きさ順に並べた時の$b$番目に小さい値を$T^\ast_{(b)}$とすると次の関係が成り立ちます。
Pr(T<T^\ast_{(b)}|H_0)=\frac{b}{B+1}
したがって、モンテカルロ$p$値と呼ばれる値が得られます。
p=Pr(T \geq t|H_0) \approx \frac{k+1}{B+1}
ブートストラップ法を用いたモンテカルロ$p$値は次にように算出されます。
p_{mc} = \frac{1}{B+1} \left( 1+\sum_{b=1}^B I\{t^{\ast b } \geq t\} \right)
Rでモンテカルロ検定
ここで用いるのはbootパッケージに含まれるデータセットfirは、バルサムモミの苗木の数を記録したものです。
library(boot)
fir$count
# > [1] 0 1 2 3 4 3 4 2 2 1 0 2 0 2 4 2 3 3 4 2 1 1 1 1 4 1
# [27] 5 2 2 3 4 1 2 5 2 0 3 2 1 1 3 1 4 3 1 0 0 2 7 0
hist(fir$count)
この分布がポアソン分布に従うかどうかを検証します。
平均値と分散を確認すると近い値であることがわかります。
> mean(fir$count)
[1] 2.14
> var(fir$count)
[1] 2.408571
それではブートストラップ法によるモンテカルロ$p$値を算出し、ポアソン分布に従うかどうかを確認したいと思います。
set.seed(42)
B <- 1999
stat.s <- numeric(B)
s <- sum(fir$count)
stat <- function(d){49*var(d)/mean(d)}
y.s <- rmultinom(B, size = s, prob = rep(1/50,50))
for (b in 1:B) {
stat.s[b] <- stat(y.s[,b])
}
(1+sum(stat.s >= stat(fir$count)))/(B+1)
> [1] 0.2775
$p_{mc}=0.2775$とポアソン分布に従うとは言えないと結論づけがることができます。
並べ替え検定
基本的な２標本問題を考えます。
母集団$F$と$G$からの独立な無作為標本$x_1,\cdots,x_m \overset{iid}{\sim}F$と$y_1,\cdots,y_n \overset{iid}{\sim}G$に基づいて、FとGに差がないという次の帰無仮説を検定します。
H_0:F=G
ここで、もし帰無仮説が正しければ、２つの標本$x_i,y_i$はどちらの母集団から抽出されたものなのか区別がつかないはずです。
$x_i$と$y_i$を混合した$m+n$個のデータを混合します。
その中から$m$個の標本を非復元抽出して得られて$X^\ast_1, \cdots ,X^\ast_m$を一群から無作為標本とし、残りの$Y^\ast_1, \cdots ,Y^\ast_m$を第2群からの無作為標本とみなした２群のデータは、帰無仮説を反映した抽出であると言えます。
このように抽出した2標本から算出した統計量$\hat \theta^\ast$を求め、繰り返し計算したものを$\hat \theta^{\ast 1}, \cdots ,\hat \theta^{\ast B}$とします。
ここで、このように混合された標本において、$x_i$から選択された標本かどうかのラベルのベクトルを$g=(g_1,\cdots,g_N)$とします。
このベクトルの全パターンは、以下のように表現できます。
U =
\begin{pmatrix}
N \\
m \\
\end{pmatrix}
=\frac{N!}{m!n!}
並べ替え統計量の値も全部で$U$通り存在し、これらの全ての値に確率$\frac{1}{U}$を与える分布を$\hat \theta^\ast$の並べ替え分布と呼びます。
帰無仮説$H_0:F=G$の元では、並べ替え$p$値は$\hat \theta^\ast $が$\hat \theta$以上となる確率として次式のように定義されます。
p_{perm}=Pr_{perm} \{\hat \theta^\ast \geq \hat \theta\}
= \frac{ \# \{\hat \theta^\ast \geq \hat \theta\} }{\begin{pmatrix}
N \\
m \\
\end{pmatrix}}
Nが大きい際には計算量が多くなってしまうため、次式のようなモンテカルロ近似値によって算出されます。
\hat p_{perm}=\frac{1}{B} \sum_{b=1}^B I\{\hat \theta^\ast \geq \hat \theta\}
アルゴリズム
1.１群と2群のデータを混合した$z=(x,y) $から大きさ$m$の標本$x^{\ast b}$を無作為復元抽出し、抽出されなかった残りの標本を$y^{\ast b}$とする。
2.手順1を$B$回繰り返し行い、それぞれのブートストラップ標本において、ブートストラップ検定統計量を計算する。
\hat \theta^\ast = t(x^{\ast b}, y^{\ast b})
3.ブートストラップp値のモンテカルロ近似値を次式により計算する。
\hat p_{perm} = \frac{1}{B} \sum_{b=1}^B I \left\{\hat \theta^\ast \geq \hat \theta\ \right \}
4.与えられる有意水準$\alpha$と$\hat p_{perm}$から、$H_0$の棄却・採択を決定する。
Rで並べ替え検定
ここで用いるデータはbootstrapパッケージに含まれる、手術を行なったマウスデータと手術を行なわなかったマウスのデータを用います。
それぞれ平均生存時間を$\mu_T,\mu_C$として、帰無仮説$H_0:\mu_T=\mu_c$に対する検定を考えます。
library(bootstrap)
mouse.t #手術を受けたマウスのデータ
# [1] 94 197 16 38 99 141 23
mouse.c #手術を受けなかったマウスのデータ
# [1] 52 104 146 10 50 31 40 27 46
mean(mouse.t)-mean(mouse.c)
# [1] 30.63492
差分だけをみると大きな差であるように感じますが、p値を計算して確認をしてみます。
並べ替え統計量$\hat \theta^\ast$は、ブートストラップ標本２群の差とします。
# 並べ替え検定
N <- 2000
mean(mouse.t)-mean(mouse.c)
data_join <- c(mouse.t, mouse.c)
theta <- numeric(N)
for(b in 1:N){
idx <- sample(1:16, 7)
c <- (1:16)[-idx]
data1_boot <- data_join[idx]
data2_boot <- data_join[c]
theta[b] <- mean(data1_boot)-mean(data2_boot)
}
# 並べ替えp値の算出
diff_tc <- mean(mouse.t)-mean(mouse.c)
p_perm <- sum(theta >= diff_tc)/N
p_perm
# [1] 0.144
有意水準を$\alpha=0.05$と考えると、有意な差ではないと言えそうです。
なお、通常のt検定を行うと、次のような結果になっている。
> t.test(mouse.t, mouse.c, var.equal=T)
Two Sample t-test
data: mouse.t and mouse.c
t = 1.1214, df = 14, p-value = 0.281
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-27.95786 89.22770
sample estimates:
mean of x mean of y
86.85714 56.22222
ブートストラップ検定
分布$F,G$ともつ母集団からの無作為標本$x,y$に基づいて、$H_0:F=G$を検定する問題を考えます。
$x$と$y$を混合した標本を$z$と置き、$x,y,z$に対応する確率変数を$X,Y,Z$、検定統計量を$t(Z)=t(X,Y)$すると、$p$値は次のように求めることができる。
p=Pr\{t(Z^\ast) \geq t(z) | H_0 \}
ここで、$t(z)$は検定統計量$t(Z)$の観測値であり、$Z^\ast$は帰無仮説$H_0$により定まる分布$F_0$に従う確率変数である。
ブートストラップ検定では、次式のように$F_0$を$\hat F_0$で推定した場合の$p$値を近似計算することになります。
p_{boot}=Pr\{t(Z^\ast) \geq t(z) | \hat F_0 \}
ここで、$\hat F_0$は$z$に基づく経験分布関数とすればよく、ブートストラップ法を用いたブートストラップ$p$値は次にように算出されます。
\hat p_{perm} = \frac{1}{B} \sum_{b=1}^B I \left\{ t(z^{\ast b}) \geq t(z) \right \}
アルゴリズム
1.$z={x,y}$から大きさm,nの標本$x^{\ast b}, y^{\ast b}$を無作為復元抽出し、これを$z^{\ast b}=( x^{\ast b}, y^{\ast b} )$とおく。
2.手順1を$B$回繰り返し行い、ブートストラップ検定統計量を計算する。
\hat \theta^{\ast b} = t(z^{\ast b})=t(x^{\ast b}, y^{\ast b})
3.ブートストラップp値のモンテカルロ近似値を次式により計算する。
\hat p_{boot} = \frac{1}{b} \sum_{b=1}^B I \left\{\hat \theta^{\ast b} \geq \theta_0 \right \}
4.与えられる有意水準$\alpha$と$\hat p_{boot}$から、$H_0$の棄却・採択を決定する。
Rでブートストラップ検定
ここでも並べ替え検定で用いたマウスのデータを用います。
それぞれ平均生存時間を$\mu_T,\mu_C$として、帰無仮説$H_0:\mu_T=\mu_c$に対する検定を考えます。
統計量はブートストラップ法により得られた手術群と非手術群に対応する値の平均値の差とします。
# 平均値の差に基づくブートストラップ検定
library(bootstrap)
N <- 2000
data_join <- c(mouse.t, mouse.c)
theta <- numeric(N)
for(b in 1:N){
data1_boot <- sample(data_join, 7, replace = TRUE)
data2_boot <- sample(data_join, 9, replace = TRUE)
theta[b] <- mean(data1_boot)-mean(data2_boot)
}
# ブートストラップp値の算出
diff_tc <- mean(mouse.t)-mean(mouse.c)
p_boot <- sum(theta >= diff_tc)/N
p_boot
# [1] 0.1175
ブートストラップ統計量を可視化してみます。
パーセンタイル信頼区間に基づく検定
パラメータ$\theta$の推定値$\hat \theta$が正規分布にしたがうとすると、$\theta$の$100(1-2\alpha)%$信頼区間は次のようになります。
(\hat \theta_{lo}, \hat \theta_{up})=(\hat \theta-z_{1-\alpha}\sigma, \hat \theta-z_\alpha\sigma)
ここで、$\hat \theta^\ast \sim N(\hat \theta_{lo},\sigma^2)$が成立するとすると、次にような関係式が得られます。
Pr[\hat \theta^\ast \ge \hat \theta|\hat \theta_{lo}]=Pr\left[ \frac{\hat \theta^\ast-\hat \theta_{lo}}{\sigma} \ge \frac{\hat \theta-\hat \theta_{lo}}{\sigma} \right]=
Pr[Z \ge z_{1-\alpha}]=\alpha
ここで、$Z$は$N(0,1)$に従う確率変数です。$\theta_{up}$においても同様の関係式が得られます。
$\theta < \hat \theta_{lo}$の時、実際に観測された値$\hat \theta$以上の値を観測する確率は上記の式から次にように表現されます。
Pr[\hat \theta^\ast \ge \hat \theta|\hat \theta_{lo}-\epsilon]<\alpha
ここで、$\epsilon>0$です。この時、帰無仮説は棄却されるべきとなります。
同様に、$\theta > \hat \theta_{up}$の時、実際に観測された値$\hat \theta$以下の値を観測する確率は次にように表現されます。
Pr[\hat \theta^\ast \ge \hat \theta|\hat \theta_{lo}-\epsilon]>\alpha
この時、帰無仮説は棄却されるべきとなります。
B組のブートストラップ標本によって得られるパーセンタイル信頼区間は$(\hat \theta^\ast_{B\alpha}, \hat \theta^\ast_{B(1-\alpha)})$であり、$\hat \theta^\ast_{B\alpha}=\theta_0$となるようなパーセンタイル法による推定値は次のように表現できます。
\hat p_{boot} = \alpha_P=\frac{1}{B} \sum_{b=1}^B I \{\hat \theta^{\ast b} < \hat \theta_0 \}
アルゴリズム
1.${x,y}$からそれぞれ大きさm,nの標本$x^{\ast b}, y^{\ast b}$を無作為復元抽出する。
2.手順1を$B$回繰り返し行い、ブートストラップ検定統計量を計算する。
\hat \theta^{\ast b} = t(z^{\ast b})=t(x^{\ast b}, y^{\ast b})
3.ブートストラップp値のモンテカルロ近似値を次式により計算する。
\hat p_{boot} = \frac{1}{b} \sum_{b=1}^B I \left\{ \hat \theta^{\ast b} < \theta_0 \right\}
4.与えられる有意水準$\alpha$と$\hat p_{boot}$から、$H_0$の棄却・採択を決定する。
Rでブートストラップ検定
ここでも並べ替え検定で用いたマウスのデータを用います。
それぞれ平均生存時間を$\mu_T,\mu_C$として、帰無仮説$H_0:\mu_T=\mu_c$に対する検定を考えます。
統計量は$x^{\ast b}$,$y^{\ast b}$のそれぞれの平均値の差とします。
# パーセンタイル信頼区間によるp値
B <- 2000
theta1 <- numeric(B)
theta2 <- numeric(B)
for(b in 1:B){
x_boot <- sample(mouse.t, replace = TRUE)
y_boot <- sample(mouse.c, replace = TRUE)
theta1[b] <- mean(x_boot)-mean(y_boot)
}
boot_alpha <- sum(theta1<0)/B
boot_alpha
# [1] 0.1195
こちらの場合でも、有意な差はないと判断できるようです。
ブートストラップ統計量の分布も確認してみます。
終わりに
通常の仮説検定よりも結果が直感的でわかりやすいですね。
普段はあまり仮説検定をすることはないですが、行うときには積極的に使っていこうと思います。