目次
- binom.test()関数の基本レシピ
- binom.test()関数が出力するオブジェクトについて
- 複数の比率データに対する二項検定を一括で行い、結果オブジェクトを再利用する（mapplyとの併せ技）
binom.test()関数の基本レシピ
Rで二項検定を行うための関数が、binom.test()である。まずは基本的な使い方をおさえよう、よーしパパ頑張っちゃうゾ。
## 試行回数100回、うち成功20回のとき
> binom.test(c(20,80)) ##引数を長さ2のベクトルとして与える
Exact binomial test
data: c(20, 80)
number of successes = 20, number of trials = 100, p-value = 1.116e-09
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.1266556 0.2918427
sample estimates:
probability of success
0.2
## 試行回数 80回、うち成功20回のとき
> binom.test(20,80) ##2つの引数をそのままコンマで区切って与える
Exact binomial test
data: 20 and 80
number of successes = 20, number of trials = 80, p-value = 8.581e-06
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.1598796 0.3593635
sample estimates:
probability of success
0.25
このように人間にとっては可読性の高い文字列として、検定結果が出力される。注意しなければならないのは、引数2つをベクトル形式 c(数字, 数字) で与えると「成功した数, 失敗した数」になるが、実数値のままコンマで併記すると「成功した数, 試行の数」になり、結果が変化する。
binom.test()関数が出力するオブジェクトについて
では次に、この検定結果がどのような要素を含んでいるのか、それらをいかに取り出すかを検討してみよう。
> example <- binom.test(c(20,80)) ##検定結果をオブジェクト"example"に代入
> names(example) ##"example"に含まれている各要素の名前を吐かせる
[1] "statistic" "parameter" "p.value" "conf.int" "estimate"
[6] "null.value" "alternative" "method" "data.name"
##それぞれの中身を覗いてみる
> example$statistic
number of successes
20
> example$parameter
number of trials
100
> example$p.value
[1] 1.115909e-09
> example$conf.int
[1] 0.1266556 0.2918427
attr(,"conf.level")
[1] 0.95
> example$estimate
probability of success
0.2
> example$null.value
probability of success
0.5
> example$alternative
[1] "two.sided"
> example$method
[1] "Exact binomial test"
> example$data.name
[1] "c(20, 80)"
上から順に、「成功数」「試行数」「P値」「95%信頼区間」「確率の推定値」「帰無仮説」「代替仮説（デフォルトでは両側）」「検定法の名前」「最初に与えたデータ」である。
例として、信頼区間の値を取り出す方法を考えてみる。上で出力した "example$conf.int" は配列の形をしているので、信頼区間の上限もしくは下限の値を単独で出力したいならば、その配列内の要素を添え字で指定しよう。
> example$conf.int[1] ##95%信頼区間の下限の値を出力する
[1] 0.1266556
> example$conf.int[2] ##95%信頼区間の上限の値を出力する
[1] 0.2918427
この、names()関数で結果のオブジェクトに含まれる要素名を調べ、"オブジェクト名$要素名"でダイレクトに出力させるという方法は、ある計算結果を次の計算へ受け渡したい場合や、計算結果を用いてグラフを自動描画したい場合によく使われる。これが奴等のやり口なのである。
複数の比率データに対する二項検定を一括で行い、結果を出力する（mapplyとの併せ技）
ここまで見てきたようにbinom.test()関数は少なくとも2つの引数（成功数, 例数）を食べ、複雑な構造の結果オブジェクトを吐き出す（信頼区間の幅や検定の片側・両側等も指定できるので、場合によってはさらに複雑となる）。ここで、対象としたい比率データ自体がたくさんある事例を考える。
> (n.success <- c (48, 60, 5, 199) ) ##成功した数を与える（×4つ）
[1] 48 60 5 199
> (n.trial <- c (350, 102, 90, 233) ) ##試行の数を与える（×4つ）
[1] 350 102 90 233
>
> binom.test(n.success[1],n.trial[1]) ##引数を1つずつ与える場合は問題ない
Exact binomial test
data: n.success[1] and n.trial[1]
number of successes = 48, number of trials = 350, p-value < 2.2e-16
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.1028793 0.1776893
sample estimates:
probability of success
0.1371429
> binom.test(n.success,n.trial) ##一度に複数回の処理を行わせようとすると......
以下にエラー binom.test(n.success, n.trial) : 'x' の長さが不正確です
これはbinom.test()関数自体が、そもそも一括処理に対応していないせいである。強引に for (i in 1:4) { print( binom.test( n.success[i], n.trial[i] ) ) } などとしても良いのだが、Rで繰り返し文を濫用するのは計算高速化の面で不向きという意見もあるし、面倒くさい。
こんなときは、apply関数ファミリーの一種であるmapply()を用いるのが吉。
> result <- mapply( binom.test, n.success, n.trial ) ##( 適用する関数の名前、引数1、引数2、3つめ以降の引数が必要なら以下同様)
> result
[,1] [,2]
statistic 48 60
parameter 350 102
p.value 3.704736e-46 0.09183768
conf.int Numeric,2 Numeric,2
estimate 0.1371429 0.5882353
null.value 0.5 0.5
alternative "two.sided" "two.sided"
method "Exact binomial test" "Exact binomial test"
data.name "dots[[1L]][[1L]] and dots[[2L]][[1L]]" "dots[[1L]][[2L]] and dots[[2L]][[2L]]"
[,3] [,4]
statistic 5 199
parameter 90 233
p.value 7.532842e-20 1.450662e-29
conf.int Numeric,2 Numeric,2
estimate 0.05555556 0.8540773
null.value 0.5 0.5
alternative "two.sided" "two.sided"
method "Exact binomial test" "Exact binomial test"
data.name "dots[[1L]][[3L]] and dots[[2L]][[3L]]" "dots[[1L]][[4L]] and dots[[2L]][[4L]]"
ちょっと結果が判り難いが、正しく4回分の計算が実行されている。次いで、この結果オブジェクトから欲しい部分だけを抜粋する方法。
> dim(result) ##このオブジェクトは配列形式データなので、まず次元を調べる。
[1] 9 4
##たとえば上記3回目の計算結果について、信頼区間の情報を取り出す
> result[4 ,3]
$conf.int
[1] 0.01828255 0.12490267
attr(,"conf.level")
[1] 0.95
##3回目の計算における信頼区間の下限値だけを取り出す
> result[4,3][[1]][1] ##アドレス指定方法その1
[1] 0.01828255
> result[,3]$conf.int[1] ##アドレス指定方法その2
[1] 0.01828255
##それぞれの試行における上限値・下限値の順に、繋がった1本のベクトルとして出力する
> unlist(result[4,])
[1] 0.10287932 0.17768928 0.48643537 0.68477386 0.01828255 0.12490267 0.80210056 0.89677988
> c.kukan <- matrix ( unlist(result[4,]) , byrow=F, nrow=2 ) ##上のベクトルを行列にする
[,1] [,2] [,3] [,4]
[1,] 0.1028793 0.4864354 0.01828255 0.8021006
[2,] 0.1776893 0.6847739 0.12490267 0.8967799
これを c.kukan[1,] または c.kukan[2,] などとすることで、下限値ないし上限値にまとめてアクセスすることが可能となる。
ここまでの結果を利用して、たとえば多数の比率データから信頼区間を自動計算し、その値をエラーバーとして、既存の棒グラフ上へ一括プロットするような操作が簡単に実施できる。まあ実際にはfor文など使ったほうが速いかもしれませんけどね。