Bootstrap による検定のシミュレーション

bootstrap で検定ってどうやるんでしょうか？

— R言語徹底解説8章熟読中 (@hoxo_m) 2016年3月4日

@hoxo_m なるほど…各標本点に第一または第二グループという「ラベルがついてる」状態ですよね？全データ混ぜて無作為抽出で二グループ作り平均値の差を記録→この操作をくりかえして「差」の分布を生成→この95%区間内に「観測された平均値の差」が入るかどうか調べる…のはどうでしょう？

— 久保拓弥 (@KuboBook) 2016年3月4日

サンプル混ぜすに双方からリサンプリングして差を出すを繰り返して信頼区間がゼロまたぐかを出すとかも考えられるかも

— R言語徹底解説8章熟読中 (@hoxo_m) 2016年3月4日

というわけで調べてみた。

kubo_boot <- function(x, y, B=1000) {
  diff <- mean(x) - mean(y)
  replicate(B, {
    mix <- c(x,y)
    data_list <- split(mix, sample(1:2, length(mix), replace=TRUE))
    xb <- data_list[[1]]
    yb <- data_list[[2]]
    mean(xb) - mean(yb)
  }) -> boot_diff
  sum(boot_diff <= -abs(diff) | abs(diff) <= boot_diff) / B
}

hoxo_boot <- function(x, y, B=1000) {
  diff <- mean(x) - mean(y)
  replicate(B, {
    xb <- sample(x, replace=TRUE)
    yb <- sample(y, replace=TRUE)
    mean(xb) - mean(yb)
  }) -> boot_diff
  min(sum(boot_diff <= 0), sum(boot_diff >= 0)) / B
}

久保先生の方法と私の方法を調べてみる。

検定として正しければ、p値は一様分布するはずである。

replicate(1000, {
  x <- rnorm(100)
  y <- rnorm(100)
  kubo_boot(x, y)
}) -> pvalues_normal_kubo
hist(pvalues_normal_kubo)

mean(pvalues_normal_kubo <= 0.05)

## [1] 0.047

replicate(1000, {
  x <- rnorm(100)
  y <- rnorm(100)
  hoxo_boot(x, y)
}) -> pvalues_normal_hoxo
hist(pvalues_normal_hoxo)

mean(pvalues_normal_hoxo <= 0.05)

## [1] 0.083

第一種過誤率は久保手法では 0.05 くらいだが、私の手法だとずれている。

ゆがんだ分布に対しても調べてみる。

replicate(1000, {
  x <- rlnorm(100)
  y <- rlnorm(100)
  kubo_boot(x, y)
}) -> pvalues_log_normal_kubo
hist(pvalues_log_normal_kubo)

mean(pvalues_log_normal_kubo <= 0.05)

## [1] 0.052

replicate(1000, {
  x <- rlnorm(100)
  y <- rlnorm(100)
  hoxo_boot(x, y)
}) -> pvalues_log_normal_hoxo
hist(pvalues_log_normal_hoxo)

mean(pvalues_log_normal_hoxo <= 0.05)

## [1] 0.102

こちらも第一種過誤率は久保手法はよくて私の方はずれてる。

検出力 0.8 が出るか調べる。

delta <- power.t.test(100, sig.level=0.05, power=0.8)$delta

replicate(1000, {
  x <- rnorm(100)
  y <- rnorm(100, delta)
  kubo_boot(x, y)
}) -> pvalues_power_kubo
mean(pvalues_power_kubo <= 0.05)

## [1] 0.798

replicate(1000, {
  x <- rnorm(100)
  y <- rnorm(100, delta)
  hoxo_boot(x, y)
}) -> pvalues_power_hoxo
mean(pvalues_power_hoxo <= 0.05)

## [1] 0.869

久保先生は 0.8 くらいになるが私の方はずれてる。

結論：久保先生の手法が良い。