←戻る

二項検定の作業例

ここから統計学が始まったとされる出来事を、David S. Salsburg『The Lady Tasting Tea』/邦訳『統計学を拓いた異才たち』を参考に、以下の例で計算していく。

「ある婦人に対して、20杯のミルクティーについてミルクを先に入れたか紅茶を先に入れたかを当てさせたところ、16回が正答であった。」

n <- 20   # 問題数
x <- 16   # 正答数

以下の統計分析については佐藤信『官能検査入門』および『統計的官能検査法』によるところが大きいので、より詳しく知りたい向きはそちらをどうぞ。

この識別は有意か?

この婦人にミルクが先か紅茶が先かを言い当てる能力があるのだろうか、それとも今回は単なるまぐれだったのか、ということを調べる。ミルクと紅茶をでたらめ(無作為)に選んだのであれば、(p_0=½)である。この無作為な確率で選んだときに、まぐれで20回中16回以上正答が起こる確率(p)を計算する。

p0 <- 1/2   # 無作為の場合の確率
p <- pbinom(x-1, n, p0, lower.tail=FALSE)
print(p)
## [1] 0.005908966

(p)の値を見て、基準とする危険率(たとえば(\alpha=.05)や(\alpha=.01))より低ければ、「このように低い確率のことを実現するのはまぐれ当たりでは難しいので、婦人には識別能力があるのではないか」と考える。今回のケースでは(p=0.005909)と十分に低いので、婦人に識別能力があると判断する。

婦人の識別能力(p_1)の95%信頼区間を求める

今回の婦人の識別能力は(p_1=\frac{x}{n}=0.8)だったが、その両側信頼区間を計算する。

alpha = 1 - 0.95
f1 <- 2 * (x + 1)
f2 <- 2 * (n - x)
p.upper <- (f1 * qf(1-alpha/2, f1, f2)) / (f2 + f1 * qf(1-alpha/2, f1, f2))
f1 <- 2 * (n - x + 1)
f2 <- 2 * x
p.lower <- f2 / (f2 + f1 * qf(1-alpha/2, f1, f2))

を計算すると、95%信頼区間((p_\mathrm{lower}), (p_\mathrm{upper}))は(0.56, 0.94)となる。以下のようにbinom.testでも同様の計算が行える。

res <- binom.test(x, n, p=p0, alternative="two.sided", conf.level=0.95)
print(res$conf.int)
## [1] 0.563386 0.942666
## attr(,"conf.level")
## [1] 0.95

95%信頼区間は、今回同様の実験を100回繰り返したとき、そのうち95回において真の婦人の能力をとらえられる区間である。

第2種の誤りを求める

第1種の誤りは「婦人には識別能力がないのに“識別能力がある”と誤った判断をしてしまう」もので、上記で求めた(p=0.005909)がその確率になる。一方で、第2種の誤りは「婦人には識別能力があるにもかかわらず“識別能力はない”と誤った判断をしてしまう」もの。この確率は以下のように計算される。

p1 <- x / n   # 婦人の能力の推定値
beta <- pbinom(x-1, n, prob=p1, lower.tail=TRUE)

ここから、第2種の誤りをおこす確率は(\beta=0.37)となる。検定力は(1-\beta=0.63)である。

第1種の誤りと第2種の誤りを小さくしたい

識別能力も試行回数もそのままにした状態では、第2種の誤りを減らそうとすると第1種の誤りが増えてしまう。第1種・第2種の誤りを同時に小さくするためには識別能力を高くするか、試行回数を増やすしかない。

識別能力がこれまで通り0.80であるとして、第1種の誤りを0.05、第2種の誤りを0.10に抑えたいときに、最低必要となる試行回数(ミルクティーのカップ数)はいくつだろうか。まずは概算からはじめる。

alpha <- 0.05
beta <- 0.10
p0 <- 0.50
p1 <- 0.80
n <- (qnorm(1-alpha) * sqrt(p0*(1-p0)) + qnorm(1-beta)*sqrt(p1*(1-p1)))^2 / (p1 - p0)^2

(n=19.8039074)となるが、この数を整数にするときの切り上げ・切り捨てで、第1種の誤り・第2種の誤りの値が変化するため、試行回数(n)、正答回数(k)ともに適宜調整が必要になる。

概算の結果から、(n=20)として第1の誤りを0.05以下にするような(k)を探すと

pbinom(14-1, 20, prob=0.5, lower.tail=FALSE)
## [1] 0.05765915
pbinom(15-1, 20, prob=0.5, lower.tail=FALSE)
## [1] 0.02069473

なので、(k=15)にすれば良さそうに思う。ただし、第2の誤りについて計算すると

pbinom(15-1, 20, prob=0.8, lower.tail=TRUE)
## [1] 0.1957922

と0.10を越えてしまう。(n=21)に増加させて

pbinom(15-1, 21, prob=0.8, lower.tail=TRUE)
## [1] 0.1085125
pbinom(15-1, 22, prob=0.8, lower.tail=TRUE)
## [1] 0.0561446

と第2種の誤りを0.10以下に抑える(n)を探す。これを交互に行って(n)と(k)を追い込んでゆく。

pbinom(15-1, 22, prob=0.5, lower.tail=FALSE)
## [1] 0.06690025
pbinom(16-1, 22, prob=0.5, lower.tail=FALSE)
## [1] 0.0262394
pbinom(16-1, 22, prob=0.8, lower.tail=TRUE)
## [1] 0.1329508
pbinom(16-1, 23, prob=0.8, lower.tail=TRUE)
## [1] 0.07150584
pbinom(16-1, 23, prob=0.5, lower.tail=FALSE)
## [1] 0.04656982

最終的に、(n=23)、(k=16)のときに第1種・第2種の誤りの両条件を満たすことができた。23杯のミルクティーのうち、16杯以上について正解できた場合、第1種の誤りは0.05以下、第2種の誤りは0.10以下の条件を満たして“利きミルクティー”に成功したと言える。

ただし、最近のコンピュータのように計算速度が十分速い場合は、前半の概算を行わず後半の追い込み部分だけでも所望の(n)と(k)が計算できる。以下のプログラムで(n)と(k)を求められる。

alpha <- 0.05
beta <- 0.10
p0 <- 0.50
p1 <- 0.80

n <- 1
k <- 1
a <- pbinom(q=k-1, size=n, prob=p0, lower.tail=FALSE)
b <- pbinom(q=k-1, size=n, prob=p1, lower.tail=TRUE)
while (a > alpha || b > beta) {
  if (a > alpha) {
  k <- k+1
    a <- pbinom(q=k-1, size=n, prob=p0, lower.tail=FALSE)
    b <- pbinom(q=k-1, size=n, prob=p1, lower.tail=TRUE)
  }
  if (b > beta) {
    n <- n+1
    a <- pbinom(q=k-1, size=n, prob=p0, lower.tail=FALSE)
    b <- pbinom(q=k-1, size=n, prob=p1, lower.tail=TRUE)
  }
}

結果、(n=23)、(k=16)となった。



[←戻る](index.html)