library(tidyverse)
## 図のなかで日本語を使えるようにする
## フォントの設定はお好みで
## （Unix/Linux ではIPAexフォントのインストールが必要かも）
if (.Platform$OS.type == "windows") { # Windows
library(fontregisterer)
my_font <- "Yu Gothic"
} else if (capabilities("aqua")) { # macOS
my_font <- "HiraginoSans-W3"
} else { # Unix/Linux
my_font <- "IPAexGothic"
}
theme_set(theme_gray(base_size = 9,
base_family = my_font))11 統計的検定と平均値の比較
今回の目標
- 2つの集団の平均値が等しいかどうか統計的に判断する方法を理解しよう！
11.1 準備
必要なパッケージを読み込み、図の日本語が正しく表示されるようにする。
11.2 統計的検定
11.2.1 有意水準の意味を理解する：帰無仮説が正しい場合の仮説検定シミュレーション
標準正規分布に従う変数 \(X\) について考える。つまり、
\[ X \sim \mbox{Normal}(\mu = 0, \sigma = 1) \] と表記される確率変数\(X\)を考える。
私たちが母数（パラメタ）を知らないと仮定して、この集団からランダムに標本を抽出し、以下の仮説を検証する。
- 帰無仮説：\(\mu = 0\)
- 対立仮説：\(\mu \neq 0\)
有意水準を5% に設定し、検定を行う。実際には帰無仮説のほうが正しいので、帰無仮説を棄却しないことが望ましい。母数を知らずに検定を行った場合、どのような結果が得られるだろうか。
11.2.1.1 1つの標本で検定を行う
まず、標本を抽出する。標本サイズ \(N = 20\) に設定しよう。
N <- 20
smp_1 <- rnorm(N, mean = 0, sd = 1) # 標準正規分布からの乱数生成この標本の標本平均は、
mean(smp_1)[1] 0.303926である。 この標本平均は0ではない。つまり \(\bar{x} \neq 0\) であるが、ここから、\(\mu \neq 0\)と言えるだろうか？
標本平均の分布は自由度 \(N - 1\) の \(t\) 分布に従うので、\(t\) 分布を利用して検定を行う（このような検定を、\(t\) 検定と呼ぶ）。
ある標本を利用して、「母平均が0」という帰無仮説を検証したいときは、t.test() という関数を使う。
対立仮説が「母平均は0ではない」のときは、次のようにする。
(res_1 <- t.test(smp_1))
One Sample t-test
data: smp_1
t = 1.4158, df = 19, p-value = 0.173
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
-0.1453650 0.7532171
sample estimates:
mean of x
0.303926 この結果の読み方を説明しよう。まず、data: smp_1 というところに、検定に利用した標本が示される。次の行を見ると、検定統計量 \(T\) の値が t = で 示された後、自由度 (df) と\(p\) 値 (p-value) が表示される（\(p\)値については「計量経済学」の講義で説明する [説明しないかもしれないが、\(p\) 値は使わない方針なので問題ない]）。検定統計量の絶対値 \(|T|\) が、自ら決めた有意水準における臨界値の絶対値 |\(c\)| より大きいとき、帰無仮説を棄却する。検定統計量自体は、次のようにすれば取り出せる。
res_1$statistic t
1.41584 また、自由度 \(N - 1\) の \(t\) 分布で、有意水準が5%のときに必要な臨界値は、
qt(p = 0.05 / 2, df = N - 1)[1] -2.093024である。この2つの値の絶対値同士を比べると、\[|T| < |c|\] だから、帰無仮説は棄却されない。つまり、\(\mu = 0\) を否定する証拠はない。
実際、\(\mu = 0\)ということを私たちは知っているので、これは妥当な結論である。標本抽出と仮説検定を繰り返したら、どんなことが起きるだろうか？
Important
\(p\)値というのは統計学の中でもとりわけ難しい概念である。大学教員のなかにものその意味を正しく理解していない人がいるくらいである（皆が統計学を専門としているわけではないのでしかたない）。よって、かなり真剣に統計学を勉強する気がある者以外は、\(p\)値には関わらないほうがよいだろう。意味もわからずに卒業論文に\(p\)値を掲載すると、質疑応答で困った事態になることが予想される。
「計量経済学」を受講してくれれば説明はする予定である（気が変わるかもしれないが）。
11.2.1.2 仮説検定のシミュレーション
まず、標本抽出と仮説検定を1回行い、帰無仮説 (null hypothesis) を棄却しないときは null を、帰無仮説を棄却して対立仮説 (alternative hypothesis) を採用するときは alt を返す関数を作る。 有意水準 (level) と標本サイズは自分で設定できるようにしておく。
smp_test <- function(N = 10, level = 0.05) {
## 標本を抽出し、「母平均 = 0」という仮説を検定する関数
## 引数：N = 標本サイズ（既定値は10）
## level = 有意水準（既定値は0.05）
smp <- rnorm(N, mean = 0, sd = 1) # 標準正規分布から乱数を生成する
test <- t.test(smp)
judge <- abs(test$statistic) > abs(qt(p = level / 2, df = N - 1))
res <- ifelse(judge, "alt", "null")
return(res)
}
Note
abs() は絶対値 (absolute values) を計算する関数である。
試しにこの関数を使ってみると、
smp_test(N = 20, level = 0.05) t
"alt" というように、1つの標本からどちらの仮説を採用したかがわかる。
標本サイズ \(N = 20\)、有意水準5% で、この作業を10,000回繰り返してみよう。関数を繰り返し実行する際は、replicate() を使うのが便利である。
n_sims <- 1e4
sim_1 <- replicate(n_sims, smp_test(N = 20, level = 0.05))sim_1 の中に、10,000個の検定結果が保存されたはずである。最初の5個だけ確認してみよう。
sim_1[1 : 5] t t t t t
"null" "null" "null" "null" "null" 結果が保存されていることが確認できる。
では、それぞれの仮説はそれぞれ何回ずつ採用されただろうか？表 (table) にしてみよう。
table(sim_1)sim_1
alt null
477 9523 割合を表示するために、シミュレーション回数で割ると、
table(sim_1) / n_simssim_1
alt null
0.0477 0.9523 となる。 このように、対立仮説が約5%、帰無仮説が95%の検定で採用されている。
私たちは、上で \(\mu=0\) と設定しているので、帰無仮説こそが正しい仮説である。しかし、帰無仮説が正しくても、5%分の検定では、誤って帰無仮説を棄却し、対立仮説を採用してしまう。この誤りのことを、「第1種の誤り (Type I error)」と呼ぶ。第1種の誤りの確率は、有意水準に一致する。私たちは有意水準5%を選んだので、約5%の検定で第1種の誤りをおかしてしまう。
有意水準を1%にして、本当にそうなるか確かめてみよう。
sim_2
alt null
0.0096 0.9904 第1種の誤り（null が正しいのに alt を採用）の割合は、約1%である。
有意水準を14%にすると、
sim_3
alt null
0.1414 0.8586 第1種の誤り（null が正しいのに alt を採用）の割合は、約14%になる。
実習課題： 第1種の誤りは、標本サイズの影響を受けるだろうか? 有意水準を固定し、標本サイズ \(N\) を10, 20, 100, 200 と変えて、第1種の誤りの割合をシミュレーションしてみよう！
11.2.2 危険率、検出力、標本サイズ：対立仮説が正しい場合の仮説検定シミュレーション
平均3、標準偏差2の正規分布に従う変数 \(Y\) について考える。すなわち、
\[ Y \sim \mbox{Normal}(\mu = 3, \sigma = 2) \] となる変数 \(Y\) について考える。
私たちが母数を知らないと仮定して、この集団からランダムに標本を抽出し、以下の仮説を検証する。
- 帰無仮説：\(\mu = 0\)
- 対立仮説：\(\mu \neq 0\)
有意水準を5% に設定し、検定を行う。実際には対立仮説のほうが正しいので、帰無仮説を棄却したい。母数を知らずに検定を行った場合、どのような結果が得られるだろうか。
11.2.2.1 1つの標本で検定を行う
まず、標本を抽出する。標本サイズ \(N = 20\) に設定しよう。
N <- 20
smp_2 <- rnorm(N, mean = 3, sd = 2) # N(3, 2) からの乱数生成この標本の標本平均は、
mean(smp_2)[1] 2.338013である。 この標本平均は0ではない。つまり \(\bar{y} \neq 0\) であるが、ここから、\(\mu \neq 0\)と言えるだろうか？
標本平均の分布は自由度 \(N - 1\) の \(t\) 分布に従うので、\(t\) 分布を利用して検定を行う。
先ほど同様、t.test() を使う。仮説は先ほどと同じなので、次のようにする。
(res_2 <- t.test(smp_2))
One Sample t-test
data: smp_2
t = 4.8381, df = 19, p-value = 0.0001143
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
1.326553 3.349472
sample estimates:
mean of x
2.338013 検定統計量 \(T\) の値は、
res_2$statistic t
4.838076 である。 また、自由度 \(N - 1\) の \(t\) 分布で、有意水準が5%のときに必要な臨界値は、
qt(p = 0.05 / 2, df = N - 1)[1] -2.093024である。この2つの値の絶対値同士を比べると、\[|T| > |c|\] だから、帰無仮説は棄却される。
実際、\(\mu = 3\)ということを私たちは知っているので、これは妥当な結論であるが、標本抽出と仮説検定を繰り返したら、どんなことが起きるだろうか？
11.2.2.2 仮説検定のシミュレーション
まず、標本抽出と仮説検定を1回行い、帰無仮説 (null hypothesis) を棄却しないときは null を、帰無仮説を棄却して対立仮説 (alternative hypothesis) を採用するときは alt を返す関数を作る。 先ほど作った関数を改良し、母平均と母標準偏差も自分で設定できるようにする。
smp_test_2 <- function(N = 10, level = 0.05, mu = 0, sigma = 1) {
## 標本を抽出し、「母平均 = 0」という仮説を検定する関数
## 引数：N = 標本サイズ（既定値は10）
## level = 有意水準（既定値は0.05）
## mu = 母平均（既定値は0）
## sigma = 母標準偏差（既定値は1）
smp <- rnorm(N, mean = mu, sd = sigma) # 正規分布から乱数を生成する
test <- t.test(smp)
judge <- abs(test$statistic) > abs(qt(p = level / 2, df = N - 1))
res <- ifelse(judge, "alt", "null")
return(res)
}試しにこの関数を使ってみると、
smp_test_2(N = 20, level = 0.05, mu = 2, sigma = 3) t
"alt" というように、1つの標本からどちらの仮説を採用したかがわかる。
標本サイズ \(N = 20\)、有意水準5% で、この作業を10,000回繰り返してみよう。
n_sims <- 1e4
sim_11 <- replicate(n_sims, smp_test_2(N = 20, level = 0.05, mu = 2, sigma = 3))sim_11 の中に、10,000個の検定結果が保存されたはずである。最初の5個だけ確認してみよう。
sim_11[1:5] t t t t t
"alt" "alt" "alt" "alt" "alt" 結果が保存されていることが確認できる。
では、それぞれの仮説はそれぞれ何割ずつ採用されただろうか？表 (table) にしてみよう。
table(sim_11) / n_simssim_11
alt null
0.8104 0.1896 このように、対立仮説が約81%、帰無仮説が約19%の検定で採用されている。
私たちは、上で \(\mu=2\) と設定しているので、対立仮説こそが正しい仮説である。しかし、対立仮説が正しくても、19%分の検定では、帰無仮説を棄却せず、対立仮説を採用できない。この誤りのことを、「第2種の誤り (Type II error)」と呼ぶ。第2種の誤りの確率は、どうやって決まるのだろうか？
有意水準を1%にして、第2種の誤りの割合を計算してみよう。
sim_12 <- replicate(n_sims, smp_test_2(N = 20, level = 0.01, mu = 2, sigma = 3))
table(sim_12) / n_simssim_12
alt null
0.5577 0.4423 第2種の誤り（null が正しいのに alt を採用 alt が正しいのに null を採用）の割合は、約44%である。有意水準を小さくしたところ、第2種の誤りの割合が大きくなった。
今度は、有意水準を10%にしてみよう。
sim_13 <- replicate(n_sims, smp_test_2(N = 20, level = 0.1, mu = 2, sigma = 3))
table(sim_13) / n_simssim_13
alt null
0.8902 0.1098 第2種の誤り（null が正しいのに alt を採用）の割合は、約11%である。有意水準を大きくしたところ、第2種の誤りの割合が小さくなった。
このように、第2種の誤りの確率は、有意水準と負の相関関係をもつ。第2種の誤りをおかす確率を小さくするためには、有意水準を大きくすればよい。しかし、有意水準は、それ自体が第1種の誤りの確率を表しているので、有意水準を大きくするということは、第1種の誤りの確率が大きくなるということである。
有意水準を変えずに、第2種の誤りの確率を小さくすることはできないだろうか。 第2種の誤りを小さくするには、「0との違い」を見つけ出す力を強めればよい。つまり、小さな違いでも違いとして見つけ出すことができれば、帰無仮説を棄却することができるようになり、第2種の誤りの確率は小さくなる。このような力のことを検出力 (power) と呼ぶ。検出力は標本サイズに依存する。
有意水準を5%にして、標本サイズを 10, 20, 50, 100, 200 と変えて、第2種の誤りの割合を計算してみよう。
\(N = 10\) のとき。
sim_14 <- replicate(n_sims, smp_test_2(N = 10, level = 0.05, mu = 2, sigma = 3))
table(sim_14) / n_simssim_14
alt null
0.4667 0.5333 \(N = 20\) のとき。
sim_15 <- replicate(n_sims, smp_test_2(N = 20, level = 0.05, mu = 2, sigma = 3))
table(sim_15) / n_simssim_15
alt null
0.8059 0.1941 \(N = 50\) のとき。
sim_16 <- replicate(n_sims, smp_test_2(N = 50, level = 0.05, mu = 2, sigma = 3))
table(sim_16) / n_simssim_16
alt null
0.9956 0.0044 \(N = 100\) のとき。
sim_17 <- replicate(n_sims, smp_test_2(N = 100, level = 0.05, mu = 2, sigma = 3))
table(sim_17) / n_simssim_17
alt
1 \(N = 200\) のとき。
sim_18 <- replicate(n_sims, smp_test_2(N = 200, level = 0.05, mu = 2, sigma = 3))
table(sim_18) / n_simssim_18
alt
1 このように、標本サイズを大きくすると、第2種の誤りの確率は下がる。したがって、統計的検定における誤りの確率を減らすためには、有意水準を小さくし、標本サイズを大きくすることが重要になる（しかし、このことは、標本サイズさえ大きくすれば、些細な違いも検出されてしまうということも意味しているので、実用上は注意が必要である）。
11.3 平均値の差の検定
2つのグループの平均値に差があるかどうか、統計的検定で判断しよう。これまで同様、t.test() を使って\(t\)検定を行うが、問題によって指定しなけばいけない引数が異なるので注意が必要である。
11.3.1 対応のないデータの平均値の差の検定
例題： 2つのコーヒーチェーンで、どちらのコーヒーが美味しいか調べるため、無作為に選んだ10人にD店のコーヒーを、無作為に選んだ別の10人にS店のコーヒーを飲んでもらった。10点満点で点数をつけてもらったところ、次のようなデータが得られた。
d とs がそれぞれのコーヒーに対する10人の評価である。D店のコーヒーを飲んだ10人と、S店のコーヒーを飲んだ10人は異なるので、これは対応のないデータであると考えられる。
2つの店のコーヒーの味に対する評価の平均値は、
であり、単純に比較すれば、D店のコーヒーの方が美味しいということになりそうである。しかし、この結果は、1つの標本から得られた結果に過ぎない。言い換えると、今回抽出された人たちは、たまたまD店の味を好んだだけかもしれない。そこで、統計的検定の方法を使って、2店の美味しさに違いがあるかどうか判断しよう。
各コーヒー店のコーヒーの味の評価の母数（パラメタ、真実の値）をそれぞれ \(D\) と \(S\) とすると、今回検証する仮説は、
- 帰無仮説：\(D = S\)（つまり、\(D - S = 0\)）
- 対立仮説：\(D \neq S\)（つまり、\(D - S \neq 0\)）
ということになる。これを、有意水準5% で検定する。また、2つのグループの分散は異なると仮定する（2つのグループの分散が等しいと仮定できるときは、var.equal = TRUE を指定するが、通常は母分散を知らないので等しいと仮定しない）。
以下のコマンドで Welch の \(t\) 検定を実行する。
(eg_1 <- t.test(d, s, var.equal = FALSE))
Welch Two Sample t-test
data: d and s
t = 1.6953, df = 14.848, p-value = 0.1109
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.4392736 3.8392736
sample estimates:
mean of x mean of y
7.4 5.7 検定統計量は、1.6953 である。また、この検定に必要な \(t\) 分布の自由度 (df) は、14.848 である。この自由度を使って検定に使う臨界値を求めると、
qt(p = 0.05 / 2, df = 14.848)[1] -2.133352である。\(|T| = 1.6953 < 2.133352 = |c|\) だから、帰無仮説は棄却されない。したがって、D店のコーヒーのほうが美味しいとは言えない。
ちなみに、t.test() の結果として表示される95%信頼区間から、検定結果を出すこともできる。 この例では、DとSの差の95%信頼区間が表示されており、それが \([-0.439, 3.839]\) である。この区間が0（つまり、2つの美味しさに差がない）を含んでいるので、有意水準5%で帰無仮説を棄却できない。
11.3.2 対応のあるデータの平均値の差の検定
例題： 2つのコーヒーチェーンで、どちらのコーヒーが美味しいか調べるため、無作為に選んだ10人にD店のコーヒーとS店のコーヒーを1杯ずつ飲んでもらった。それぞれのコーヒーに10点満点で点数をつけてもらったところ、次のようなデータが得られた。
survey <- tibble(
id = 1:10,
d = c(8, 7, 8, 6, 4, 8, 9, 10, 7, 7),
s = c(6, 10, 3, 10, 4, 4, 5, 7, 2, 6)
)
survey# A tibble: 10 × 3
id d s
<int> <dbl> <dbl>
1 1 8 6
2 2 7 10
3 3 8 3
4 4 6 10
5 5 4 4
6 6 8 4
7 7 9 5
8 8 10 7
9 9 7 2
10 10 7 6id は、回答者番号、d はD店のコーヒーに対する評価、s はS店のコーヒーに対する評価である。
同じ10人がD店のコーヒーとS店のコーヒーの両方を飲んで評価しているので、これは対応のあるデータであると考えられる。
2つの店のコーヒーの味に対する評価の平均値は、
であり、単純に比較すれば、D店のコーヒーのほうが美味しいということになりそうである。しかし、この結果は、1つの標本から得られた結果に過ぎない。言い換えると、今回抽出された人たちは、たまたまD店の味を好んだだけかもしれない。そこで、統計的検定の方法を使って、2店の美味しさに違いがあるかどうか判断しよう。
各コーヒー店のコーヒーの味の評価の母数（パラメタ、真実の値）をそれぞれ \(D\) と \(S\) とすると、今回検証する仮説は、
- 帰無仮説：\(D = S\)（つまり、\(D - S = 0\)）
- 対立仮説：\(D \neq S\)（つまり、\(D - S \neq 0\)）
ということになる。これを、有意水準5% で検証する。対応のあるデータであることをRに伝えるため、paired = TRUE を指定する）。
(eg_2 <- t.test(survey$d, survey$s, paired = TRUE))
Paired t-test
data: survey$d and survey$s
t = 1.6805, df = 9, p-value = 0.1272
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
-0.5883968 3.9883968
sample estimates:
mean difference
1.7 検定統計量は、1.6805 である。また、この検定に必要な \(t\) 分布の自由度 (df) は、9 である。この自由度を使って検定に使う臨界値を求めると、
qt(p = 0.05 / 2, df = 9)[1] -2.262157である。\(|T| = 1.6953 < 2.262157 = |c|\) だから、帰無仮説は棄却されない。したがって、D店のコーヒーのほうが美味しいとは言えない。
2つの美味しさの差の95%信頼区間を見てみると、 \([-0.588, 3.988]\) である。この区間が0（美味しさに差がない） を含んでいるので、有意水準5%で帰無仮説は棄却されず、保留される。