必要なパッケージを読み込み、図の日本語が正しく表示されるようにする。
library(tidyverse)
if (.Platform$OS.type == "windows") {
if (require(fontregisterer)) {
my_font <- "Yu Gothic"
} else {
my_font <- "Japan1"
}
} else if (capabilities("aqua")) {
my_font <- "HiraginoSans-W3"
} else {
my_font <- "IPAexGothic"
}
theme_set(theme_gray(base_size = 9,
base_family = my_font))標準正規分布に従う変数 \(X\) について考える。これを \[ X \sim \mbox{Normal}(\mu = 0, \sigma = 1) \] と表記する。
私たちが母数（パラメタ）を知らないと仮定して、この集団からランダムに標本を抽出し、以下の仮説を検証する。
有意水準を5% に設定し、検定を行う。実際には帰無仮説のほうが正しいので、帰無仮説を棄却しないことが望ましい。母数を知らずに検定を行った場合、どのような結果が得られるだろうか。
まず、標本を抽出する。標本サイズ \(N = 20\) に設定しよう。
N <- 20
smp_1 <- rnorm(N, mean = 0, sd = 1) # 標準正規分布からの乱数生成この標本の標本平均は、
mean(smp_1)## [1] -0.1246344である。 この標本平均は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 = -0.43663, df = 19, p-value = 0.6673
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.7220753 0.4728066
## sample estimates:
## mean of x
## -0.1246344この結果の読み方を説明しよう。まず、data: smp_1 というところに、検定に利用した標本が示される。次の行を見ると、検定統計量 \(T\) の値が t = で 示された後、自由度 (df) と\(p\) 値 (p-value) が表示される（\(p\)値については「計量経済学」の講義で説明する）。検定統計量の絶対値 \(|T|\) が、自ら決めた有意水準における臨界値の絶対値 |\(c\)| より大きいとき、帰無仮説を棄却する。 検定統計量自体は、次のようにすれば取り出せる。
res_1$statistic## t
## -0.4366334また、自由度 \(N - 1\) の \(t\) 分布で、有意水準が5%のときに必要な臨界値は、
qt(p = 0.05 / 2, df = N - 1)## [1] -2.093024である。この2つの値の絶対値同士を比べると、\[|T| < |c|\] だから、帰無仮説は棄却されない。
また、同じことだが、\(p\) 値が自分で設定した有意水準より小さいとき、帰無仮説を棄却する。ここでは、有意水準は5% (0.05) で、\(p\) 値は、
res_1$p.value## [1] 0.6672988なので、帰無仮説は棄却されない。つまり、\(\mu = 0\) を否定する証拠はない。
実際、\(\mu = 0\)ということを私たちは知っているので、これは妥当な結論である。標本抽出と仮説検定を繰り返したら、どんなことが起きるだろうか？
まず、標本抽出と仮説検定を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)) # abs() は絶対値 (absolute values)
res <- ifelse(judge, "alt", "null")
return(res)
}試しにこの関数を使ってみると、
smp_test(N = 20, level = 0.05)## t
## "null"というように、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
## 479 9521割合を表示するために、シミュレーション回数で割ると、
table(sim_1) / n_sims## sim_1
## alt null
## 0.0479 0.9521となる。 このように、対立仮説が約5%、帰無仮説が95%の検定で採用されている。
私たちは、上で \(\mu=0\) と設定しているので、帰無仮説が正しい仮説である。しかし、帰無仮説が正しくても、5%分の検定では、誤って帰無仮説を棄却し、対立仮説を採用してしまう。この誤りのことを、「第1種の誤り (Type I error)」と呼び。第1種の誤りの確率は、有意水準に一致する。私たちは有意水準5%を選んだので、約5%の検定で、第1種の誤りをおかしてしまう。
有意水準を1%にして、本当にそうなるか確かめてみよう。
sim_2 <- replicate(n_sims, smp_test(N = 20, level = 0.01))
table(sim_2) / n_sims## sim_2
## alt null
## 0.0106 0.9894第1種の誤り（null が正しいのに alt を採用）の割合は、約1%である。
有意水準を14%にすると、
sim_3 <- replicate(n_sims, smp_test(N = 20, level = 0.14))
table(sim_3) / n_sims## sim_3
## alt null
## 0.1383 0.8617第1種の誤り（null が正しいのに alt を採用）の割合は、約14%になる。
実習課題：第1種の誤りは、標本サイズの影響を受けるだろうか? 有意水準を固定し、標本サイズ \(N\) を10, 20, 100, 200 と変えて、第1種の誤りの割合をシミュレーションしてみよう！
平均3、標準偏差2の正規分布に従う変数 \(Y\) について考える。すなわち、 \[ Y \sim \mbox{Normal}(\mu = 3, \sigma = 2) \] となる変数 \(Y\) について考える。
私たちが母数を知らないと仮定して、この集団からランダムに標本を抽出し、以下の仮説を検証する。
有意水準を5% に設定し、検定を行う。実際には対立仮説のほうが正しいので、帰無仮説を棄却したい。母数を知らずに検定を行った場合、どのような結果が得られるだろうか。
まず、標本を抽出する。標本サイズ \(N = 20\) に設定しよう。
N <- 20
smp_2 <- rnorm(N, mean = 3, sd = 2) # N(3, 2) からの乱数生成この標本の標本平均は、
mean(smp_2)## [1] 3.161174である。 この標本平均は0ではない。つまり \(\bar{x} \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 = 5.6942, df = 19, p-value = 1.729e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 1.999213 4.323134
## sample estimates:
## mean of x
## 3.161174検定統計量 \(T\) の値は、
res_2$statistic## t
## 5.694179である。 また、自由度 \(N - 1\) の \(t\) 分布で、有意水準が5%のときに必要な臨界値は、
qt(p = 0.05 / 2, df = N - 1)## [1] -2.093024である。この2つの値の絶対値同士を比べると、\[|T| > |c|\] だから、帰無仮説は棄却される。
また、同じことだが、\(p\) 値は、
res_2$p.value## [1] 1.728524e-05で、有意水準である0.05より小さいので、帰無仮説は棄却される。つまり、\(\mu = 0\) を否定する証拠があるので、\(\mu \neq 0\) と考える。
実際、\(\mu = 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)) # abs() は絶対値 (absolute values)
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_sims## sim_11
## alt null
## 0.8129 0.1871このように、対立仮説が約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_sims## sim_12
## alt null
## 0.5524 0.4476第2種の誤り（null が正しいのに alt を採用 alt が正しいのに null を採用）の割合は、約45%である。有意水準を小さくしたところ、第2種の誤りの割合が大きくなった。
今度は、有意水準を10%にしてみよう。
sim_13 <- replicate(n_sims, smp_test_2(N = 20, level = 0.1, mu = 2, sigma = 3))
table(sim_13) / n_sims## sim_13
## alt null
## 0.8896 0.1104第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_sims## sim_14
## alt null
## 0.474 0.526\(N = 20\) のとき。
sim_15 <- replicate(n_sims, smp_test_2(N = 20, level = 0.05, mu = 2, sigma = 3))
table(sim_15) / n_sims## sim_15
## alt null
## 0.8138 0.1862\(N = 50\) のとき。
sim_16 <- replicate(n_sims, smp_test_2(N = 50, level = 0.05, mu = 2, sigma = 3))
table(sim_16) / n_sims## sim_16
## alt null
## 0.996 0.004\(N = 100\) のとき。
sim_17 <- replicate(n_sims, smp_test_2(N = 100, level = 0.05, mu = 2, sigma = 3))
table(sim_17) / n_sims## sim_17
## alt
## 1$ = 200$ のとき。
sim_18 <- replicate(n_sims, smp_test_2(N = 200, level = 0.05, mu = 2, sigma = 3))
table(sim_18) / n_sims## sim_18
## alt
## 1このように、標本サイズを大きくすると、第2種の誤りの確率は下がる。したがって、統計的検定における誤りの確率を減らすためには、有意水準を小さくし、標本サイズを大きくすることが重要になる（しかし、このことは、標本サイズさえ大きくすれば、些細な違いも検出されてしまうということも意味しているので、実用上は注意が必要である）。
2つのグループの平均値に差があるかどうか、統計的検定で判断しよう。これまで同様、t.test() を使って\(t\)検定を行うが、問題によって指定しなけばいけない引数が異なるので注意が必要である。
例題： 2つのコーヒーチェーンで、どちらのコーヒーが美味しいか調べるため、無作為に選んだ10人にD店のコーヒーを、無作為に選んだ別の10人にS店のコーヒーを飲んでもらった。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)d とs がそれぞれのコーヒーに対する10人の評価である。D店のコーヒーを飲んだ10人と、S店のコーヒーを飲んだ10人は異なるので、これは対応のないデータであると考えられる。
2つの店のコーヒーの味に対する評価の平均値は、
mean(d)## [1] 7.4mean(s)## [1] 5.7であり、単純に比較すれば、D店のコーヒーの方が美味しいということになりそうである。しかし、この結果は、1つの標本から得られた結果に過ぎない。言い換えると、今回抽出された人たちは、たまたまD店の味を好んだだけかもしれない。そこで、統計的検定の方法を使って、2店の美味しさに違いがあるかどうか判断しよう。
各コーヒー店のコーヒーの味の評価の母数（パラメタ、真実の値）をそれぞれ \(S\) と \(D\) とすると、今回検証する仮説は、
ということになる。これを、有意水準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%で帰無仮説を棄却できない。
例題： 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 x 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つの店のコーヒーの味に対する評価の平均値は、
mean(survey$d)## [1] 7.4mean(survey$s)## [1] 5.7であり、単純に比較すれば、D店のコーヒーのほうが美味しいということになりそうである。しかし、この結果は、1つの標本から得られた結果に過ぎない。言い換えると、今回抽出された人たちは、たまたまD店の味を好んだだけかもしれない。そこで、統計的検定の方法を使って、2店の美味しさに違いがあるかどうか判断しよう。
各コーヒー店のコーヒーの味の評価の母数（パラメタ、真実の値）をそれぞれ \(S\) と \(D\) とすると、今回検証する仮説は、
ということになる。これを、有意水準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 difference in means is not equal to 0
## 95 percent confidence interval:
## -0.5883968 3.9883968
## sample estimates:
## mean of the differences
## 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%で帰無仮説は棄却されず、保留される。