必要なパッケージを読み込む。
pacman::p_load(tidyverse,
broom)多重検定の問題を確認するために、シミュレーションを実施しよう。
まず、標本サイズ\(N\)で、\(K\)個の互いに独立な説明変数\(x_1, x_2, \dots, x_{K}\) を作る。\(N=100\)、\(K=20\)で試してみる。
set.seed(2021-11-12)
N <- 100 # 標本サイズ
K <- 20 # 説明変数の数
X <- matrix(rnorm(N * K, mean = 0, sd = 2), ncol = K)
colnames(X) <- paste0("x", 1:K)
myd <- as_tibble(X)これらの説明変数 \(x_k\) (\(k = 1, 2, \dots, K\)) とは無関係に \(y\)を作り、データフレームに加える。
myd$y <- rnorm(N)このデータを使い、以下の回帰モデルを推定する。 \[ \begin{aligned} Y_i &\sim \mbox{Normal}(\mu_i, \sigma) \\ \mu_i &= \beta_0 + \beta_1 x_{i1} + \cdots + \beta_K x_{iK} \end{aligned} \]
ここで、\(\beta_k = 0\) という帰無仮説を、すべての\(k\) (\(k = 1, 2, \dots, K\)) について個別に検定するとする。そうすると、検定の対象となるファミリーは、 \[ \mathcal{F} = \{\beta_1 = 0, \dots, \beta_K = 0\} \] となる。
このシミュレーションでは、\(y\)をどの\(x_k\)とも無関係に作っているので、すべての帰無仮説が正しい。よって、このファミリーに含まれる帰無仮説は1つも棄却したくない。
実際に回帰分析を行い、有意水準0.05で1つひとつの仮説を検定しよう。y ~ .
とすると、y
をそのデータに含まれるy以外のすべての変数に回帰する。
fit <- lm(y ~ ., data = myd)
tidy(fit) %>%
knitr::kable()| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 0.0632212 | 0.1064784 | 0.5937468 | 0.5543778 |
| x1 | 0.0013358 | 0.0477982 | 0.0279475 | 0.9777745 |
| x2 | 0.0139943 | 0.0660608 | 0.2118397 | 0.8327777 |
| x3 | -0.0006238 | 0.0479122 | -0.0130198 | 0.9896448 |
| x4 | -0.0730475 | 0.0514828 | -1.4188723 | 0.1598696 |
| x5 | -0.0054763 | 0.0532819 | -0.1027800 | 0.9183980 |
| x6 | 0.0946855 | 0.0537225 | 1.7624910 | 0.0818542 |
| x7 | 0.0183856 | 0.0552615 | 0.3327021 | 0.7402404 |
| x8 | -0.1972879 | 0.0583854 | -3.3790612 | 0.0011320 |
| x9 | -0.0418925 | 0.0522673 | -0.8015058 | 0.4252431 |
| x10 | -0.0051098 | 0.0553812 | -0.0922659 | 0.9267204 |
| x11 | -0.0301027 | 0.0641550 | -0.4692188 | 0.6402057 |
| x12 | 0.0119326 | 0.0548061 | 0.2177242 | 0.8282054 |
| x13 | 0.0048011 | 0.0515746 | 0.0930897 | 0.9260679 |
| x14 | -0.0089601 | 0.0579230 | -0.1546894 | 0.8774607 |
| x15 | -0.0299819 | 0.0608123 | -0.4930237 | 0.6233635 |
| x16 | 0.1100211 | 0.0591603 | 1.8597128 | 0.0666487 |
| x17 | 0.0204935 | 0.0603651 | 0.3394919 | 0.7351401 |
| x18 | -0.0383977 | 0.0598645 | -0.6414093 | 0.5231134 |
| x19 | -0.0253316 | 0.0610053 | -0.4152359 | 0.6790949 |
| x20 | 0.0224067 | 0.0599897 | 0.3735087 | 0.7097703 |
x8の\(p\)値 (p.value) が
0.05未満なので、\(\beta_8 = 0\)という帰無仮説が誤って棄却されてしまう。
少なくとも1つの帰無仮説を誤って棄却するかどうかは、次のようにして調べることができる。\(p\)値（p.value）が有意水準を下回っていれば、帰無仮説が棄却される。
pv <- tidy(fit) %>%
filter(term != "(Intercept)") %>%
pull(p.value)
ifelse(sum(pv < 0.05) > 0, TRUE, FALSE)## [1] TRUE
この
TRUE
が、少なくとも1つの帰無仮説を誤って棄却してしまっていることを示す。FALSE
ならすべての帰無仮説が「正しく」保留されているということを意味する。
このシミュレーションを繰り返し実施するために、関数を作ろう。
Nは標本サイズ、Kは説明変数の数（つまり、ファミリーに含まれる帰無仮説の数）、alpha
は個々の帰無仮説の検証に使う有意水準である。alphaの既定値は0.05にしておく。
sim_mc <- function(N = 100, K = 1, alpha = 0.05) {
X <- matrix(rnorm(N * K, mean = 0, sd = 2), ncol = K)
y <- rnorm(N)
fit <- lm(y ~ X)
pv <- tidy(fit) %>%
filter(term != "(Intercept)") %>%
pull(p.value)
ifelse(sum(pv < alpha) > 0, TRUE, FALSE)
}この関数を1回だけ使ってみよう。
sim_mc(N = 100, K = 20)## [1] TRUE
TRUE
という答えが返ってきたので、少なくとも1つの帰無仮説が誤って棄却されたことがわかる。
この関数を使い、シミュレーションを実施する。繰り返し回数は1万回とする（各自のパソコンの性能に応じて回数は調整されたい）。 まずは、\(K=1\)の場合（つまり、多重比較ではない場合）について確認する。
res_k1 <- replicate(10000,
sim_mc(N = 100,
K = 1,
alpha = 0.05))少なくとも1つの帰無仮説が誤って棄却される割合を計算する。
mean(res_k1)## [1] 0.0478有意水準に設定した0.05に近い値が得られた。「帰無仮説が正しいのに誤って帰無仮説を棄却する確率（危険率）」が有意水準なので、この結果は当然である。
次に、\(K=20\)の場合について確認する。
res_k20 <- replicate(10000,
sim_mc(N = 100,
K = 20,
alpha = 0.05))少なくとも1つの帰無仮説が誤って棄却される割合を計算する。
mean(res_k20)## [1] 0.5954
ファミリーの有意水準が0.59になった。1つひとつの検定で利用するalphaの値を「めったにない」と言えそうな5%にしても、全体ではめったにないとは決して言えない59%に
なってしまっており、多重比較の調整が必要であることがわかる。
ボンフェローニの方法で有意水準を補正するとどうなるだろうか。
\(K=20\)
のときにファミリーの有意水準を0.05にするには、alpha = 0.05 / 20
にすればよいはずだ。やってみよう。
bonferroni_k20 <- replicate(10000,
sim_mc(N = 100,
K = 20,
alpha = 0.05 / 20))少なくとも1つの帰無仮説が誤って棄却される割合を計算する。
mean(bonferroni_k20)## [1] 0.0484約5%になっており、多重比較の問題は解消されていることがわかる。
上では、すべての帰無仮説が正しい場合について考えた。 では、いくつかの帰無仮説が正しくない場合にはどうなるだろうか。
20個の説明変数のうち、3つは\(y\)に関係があるという状況でシミュレーションを行う。 そのための関数を作る。
sim_mc2 <- function(N = 100, K = 4, alpha = 0.05, type = 1) {
X <- matrix(rnorm(N * K, mean = 0, sd = 2), ncol = K)
y <- 0.8 * X[, 1] - 0.7 * X[, 2] + 0.4 * X[, 3] + rnorm(N)
fit <- lm(y ~ X)
pv <- tidy(fit) %>%
filter(term != "(Intercept)") %>%
pull(p.value)
if (type == 1) ifelse(sum(pv[4:K] < alpha) > 0, TRUE, FALSE)
else if (type == 2) ifelse(sum(pv[1:3] < alpha) == 3, FALSE, TRUE)
else stop("type must be either 1 or 2.")
}この関数を1回だけ使ってみよう。
set.seed(2021-11-12)
sim_mc2(N = 50, K = 20)## [1] FALSEFALSEという答えが返ってきた。正しい帰無仮説はすべて保留された。
この関数を使い、シミュレーションを実施する。繰り返し回数は1万回とする（各自のパソコンの性能に応じて回数は調整されたい）。\(K=20\)の場合について確認する。
res_k20b <- replicate(10000,
sim_mc2(N = 50,
K = 20,
alpha = 0.05))正しい帰無仮説のうち少なくとも1つが誤って棄却される割合を計算する。
mean(res_k20b)## [1] 0.4908
ファミリーの有意水準が0.49になった。1つひとつの検定で利用するalphaの値を「めったにない」と言えそうな5%にしても、全体では「めったにない」とは決して言えない52%に
なってしまっており、多重比較の調整が必要であることがわかる。
そこで、ボンフェローニの方法で有意水準を補正する。
\(K=20\) なのでalpha = 0.05 / 20
にすればよいはずだ。
bonferroni_k20b <- replicate(10000,
sim_mc2(N = 50,
K = 20,
alpha = 0.05 / 20))正しい帰無仮説のうち少なくとも1つが誤って棄却される割合を計算する。
mean(bonferroni_k20b)## [1] 0.0372約3.7%になっており、正しい帰無仮説を誤って棄却するという問題は解消されていることがわかる。
では、正しくない帰無仮説は棄却できているだろうか。
第2種の誤りの割合を計算したいので、type = 2 を指定する。
bonferroni_k20c <- replicate(10000,
sim_mc2(N = 50,
K = 20,
alpha = 0.05 / 20,
type = 2))
mean(bonferroni_k20c)## [1] 0.1912約19%について、正しくない帰無仮説を少なくとも1つは棄却し損ねていることがわかる。
ボンフェローニの補正を行う前についても同じ割合を計算してみよう。
res_k20c <- replicate(10000,
sim_mc2(N = 50,
K = 20,
alpha = 0.05,
type = 2))
mean(res_k20c)## [1] 0.0244補正前には、この種の誤りは2%ほどしかなかったことがわかる。ボンフェローニの方法を用いると、正しい帰無仮説を誤って棄却するという間違い（第1種のエラー）が減る代わりに、正しくない帰無仮説を棄却することに失敗する（第2種のエラー）が増えるということだ。
ボンフェローニの方法は、帰無仮説を棄却する基準を厳しくし過ぎる結果、正しくない帰無仮説を棄却することにも失敗している。この点を改良した、ホルムの方法を試してみよう。
ボンフェローニの方法とホルムの方法比較するために、データを生成して回帰分析の結果を保存する関数を作る。
sim_mc3 <- function(N = 50, K = 4) {
X <- matrix(rnorm(N * K, mean = 0, sd = 2), ncol = K)
y <- 0.8 * X[, 1] - 0.7 * X[, 2] + 0.3 * X[, 3] + rnorm(N)
fit <- lm(y ~ X)
tidy(fit) %>%
filter(term != "(Intercept)") %>%
select(term, estimate, p.value)
}この関数を1回だけ使ってみよう。
set.seed(2021-11-14)
res1 <- sim_mc3(N = 50, K = 20)
res1## # A tibble: 20 × 3
## term estimate p.value
## <chr> <dbl> <dbl>
## 1 X1 1.05 1.47e-12
## 2 X2 -0.711 1.02e- 8
## 3 X3 0.276 5.00e- 3
## 4 X4 -0.0177 8.44e- 1
## 5 X5 -0.0987 2.79e- 1
## 6 X6 -0.0106 8.87e- 1
## 7 X7 0.0717 4.98e- 1
## 8 X8 0.00894 9.18e- 1
## 9 X9 0.101 2.49e- 1
## 10 X10 -0.119 3.13e- 1
## 11 X11 -0.0503 5.90e- 1
## 12 X12 0.0641 5.27e- 1
## 13 X13 0.0583 5.06e- 1
## 14 X14 -0.0653 5.15e- 1
## 15 X15 0.156 4.88e- 2
## 16 X16 0.212 3.19e- 2
## 17 X17 -0.102 2.03e- 1
## 18 X18 0.0725 4.28e- 1
## 19 X19 -0.0346 6.93e- 1
## 20 X20 -0.000597 9.94e- 1有意水準0.05 で個々の帰無仮説が棄却されるかどうか確かめる。
res1$p.value < 0.05## [1] TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE
最初の3つ（正しくない帰無仮説）は棄却されている。残りの17個の帰無仮説は正しいが、そのうち2つが誤って棄却されている（X15とX16の検定結果がTRUEすなわち棄却になっている）。
ボンフェローニの方法で補正するとどうなるだろうか。
p.adjust()
関数を使うと、多重比較の補正を行うことができる。
p.adjust(res1$p.value, method = "bonferroni")## [1] 2.944308e-11 2.041376e-07 9.994142e-02 1.000000e+00 1.000000e+00
## [6] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
## [11] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 9.754977e-01
## [16] 6.374704e-01 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00これらの値が、ボンフェローニの方法で検定に使うべき\(p\)値である。これを使って検定を行う。
p.adjust(res1$p.value, method = "bonferroni") < 0.05## [1] TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE補正なしだと誤って棄却してしまった正しい帰無仮説は保留されている。代わりに、補正無しでは棄却できていた正しくない帰無仮説 (\(\beta_3 = 0\)) が保留されている。
p.adjust() でホルムの方法も利用できる。
p.adjust(res1$p.value, method = "holm")## [1] 2.944308e-11 1.939307e-07 8.994728e-02 1.000000e+00 1.000000e+00
## [6] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
## [11] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 7.803981e-01
## [16] 5.418499e-01 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00これらの値が、ホルムの方法で検定に使うべき\(p\)値である。これを使って検定を行う。
p.adjust(res1$p.value, method = "holm") < 0.05## [1] TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE補正なしだと誤って棄却してしまった正しい帰無仮説は保留されている。ボンフェローニの方法と同様に、補正無しでは棄却できていた正しくない帰無仮説 (\(\beta_3 = 0\)) が保留されている。
2つの方法の\(p\)値を比較してみよう。
cbind(p.adjust(res1$p.value, method = "bonferroni"),
p.adjust(res1$p.value, method = "holm"))## [,1] [,2]
## [1,] 2.944308e-11 2.944308e-11
## [2,] 2.041376e-07 1.939307e-07
## [3,] 9.994142e-02 8.994728e-02
## [4,] 1.000000e+00 1.000000e+00
## [5,] 1.000000e+00 1.000000e+00
## [6,] 1.000000e+00 1.000000e+00
## [7,] 1.000000e+00 1.000000e+00
## [8,] 1.000000e+00 1.000000e+00
## [9,] 1.000000e+00 1.000000e+00
## [10,] 1.000000e+00 1.000000e+00
## [11,] 1.000000e+00 1.000000e+00
## [12,] 1.000000e+00 1.000000e+00
## [13,] 1.000000e+00 1.000000e+00
## [14,] 1.000000e+00 1.000000e+00
## [15,] 9.754977e-01 7.803981e-01
## [16,] 6.374704e-01 5.418499e-01
## [17,] 1.000000e+00 1.000000e+00
## [18,] 1.000000e+00 1.000000e+00
## [19,] 1.000000e+00 1.000000e+00
## [20,] 1.000000e+00 1.000000e+00第1列がボンフェローニの方法による\(p\)値、第2列がホルムの方法による\(p\)値である。 2つの方法を比較すると、ホルムの方法の\(p\)値のほうが小さいことがわかる。 つまり、ホルムの方法のほうが、帰無仮説を棄却しやすい。 帰無仮説を棄却しにくくなるというボンフェローニの方法の欠点が改良されていることがわかる。
このシミュレーションを繰り返し実行するための関数を作る。
帰無仮説が3つ未満しか棄却されなかったときにTRUE
を返すことにする。（厳密には、これは知りたいこととは異なる。たとえば、正しくない帰無仮説がすべて保留され、正しい帰無仮説が誤って3つ棄却されるという場合も考えられる。）
sim_b_h <- function(N = 50, K = 20, alpha = 0.05, method) {
res <- sim_mc3(N = N, K = K)
sum(p.adjust(res$p.value, method = method) < alpha) < 3
}この関数を1回だけ使ってみよう。
sim_b_h(N = 100, K = 20, method = "holm")## [1] FALSE
FALSE
が返されたので、正しくない帰無仮説はすべて棄却されたことがわかる。
この関数を使い、シミュレーションを実施する。繰り返し回数は1万回とする（各自のパソコンの性能に応じて回数は調整されたい）。まず、ボンフェローニの方法を試す。
res_b <- replicate(10000,
sim_b_h(N = 50,
K = 20,
method = "bonferroni"))正しくない帰無仮説のうち、少なくとも1つが保留されてしまう割合を計算する。
mean(res_b)## [1] 0.493149.8%について、正しくない帰無仮説を棄却することに失敗している。
次に、ホルムの方法を試す。
res_h <- replicate(10000,
sim_b_h(N = 50,
K = 20,
method = "holm"))少なくとも1つの帰無仮説が誤って棄却される割合を計算する。
mean(res_h)## [1] 0.482248.2%について、正しくない帰無仮説を棄却することに失敗している。 ボンフェローニの方法よりも少しだけマシなようだ。
上で説明したとおり、分析結果に対して多重比較の補正を行う際は、p.adjust()
を使う。 ランダムに生成したデータで実行してみよう。
まず、データを作る。
set.seed(2021-11-13)
N <- 50 # 標本サイズ
K <- 8 # 説明変数の数
X <- matrix(rnorm(N * K, mean = 0, sd = 2), ncol = K)
colnames(X) <- paste0("x", 1:K)
myd <- as_tibble(X) %>%
mutate(y = 0.5 * x1 + 0.4 * x2 + rnorm(n()))回帰分析を実行する。
ols <- lm(y ~ .,
data = myd)係数の推定値、補正前の\(p\)値、ボンフェローニの方法による\(p\)値、ホルムの方法による\(p\)値を並べて表示する。
ols %>%
tidy() %>%
filter(term != "(Intercept)") %>%
select(term, estimate, p.value) %>%
mutate(Bonferroni = p.adjust(p.value, method = "bonferroni"),
Holm = p.adjust(p.value, method = "holm")) %>%
rename(`説明変数` = term,
`推定値` = estimate,
`p値` = p.value)## # A tibble: 8 × 5
## 説明変数 推定値 p値 Bonferroni Holm
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 x1 0.342 0.00120 0.00960 0.00840
## 2 x2 0.494 0.000000844 0.00000675 0.00000675
## 3 x3 0.0113 0.906 1 1
## 4 x4 0.0728 0.415 1 1
## 5 x5 -0.0914 0.244 1 1
## 6 x6 -0.0247 0.815 1 1
## 7 x7 0.0939 0.227 1 1
## 8 x8 -0.0529 0.652 1 1それぞれの\(p\)値が有意水準を下回る場合について帰無仮説を棄却する。 この例では、補正の有無によって検定結果は変わらない。