必要なサンプルサイズの大きさは（実験編）

前回の記事で、調査を行い母集団における比率を求める際に、誤差を任意の範囲内に収めるための計算について書いた。比率については、

　　 $\displaystyle n = \left( \frac{1.96}{E} \right)^2 p(1-p)$

の式の、Eの部分に収めたい範囲の誤差を代入すると求まることが分かった。

式の上では理論的に求まったものの本当にこれで誤差が所定の範囲内に収まるのか試してみたくなるのが世の常である。ということで、毎度おなじみRを用いたシミュレーションである。

前回の記事で、95%水準で誤差を5%、3%、1%に収めたいときに必要なサンプルサイズは、上の式をそれぞれについて解いて385、1068、9604だと書いた。本当にこの通りになっているか確認してみよう。

まず、二項分布で乱数を発生させる。二項分布の乱数発生の関数はrbinom(n, size, p)である。引数には順番に「観察の回数」「試行の回数」「試行における成功確率」を指定する。今やりたいことに当てはめて考えると、nにはシミュレーションの回数を、sizeにはサンプルサイズを、pには母集団における確率を指定する。シミュレーションの回数はまぁ大体10000回もやればとりあえず良いでしょう。pには計算に用いた0.5を指定しておく。

dat <- rbinom(10000,385,0.5)
head(dat)
# [1] 207 190 198 193 203 185 
length(dat)
# [1] 10000

head関数で中身を除いた結果から分かるように、成功率0.5で385回の試行における成功数についてのデータがそれぞれ発生する。データの個数は当然10000である。

このままだと、比率のデータではないが、サンプルサイズで割るとこれが比率になる。で、ここで得られた標本比率が母集団の比率である50%より5%以上ずれているかを調べたい訳なので、45%より小さいあるいは55%より大きい標本比率の数を数えることにする。その数をシミュレーションの回数で割れば一体何%ぐらいが5%よりも大きな誤差であったかが分かるはずである。

dat <- dat/385
over <- length(dat[dat > 0.55]) + length(dat[dat < 0.45])
over/10000
# [1] 0.0509

結果を見ると分かるように、標本比率の内5%ほどのデータが5%より大きい誤差であることが分かる。逆に言えば、信頼係数である95%の標本は誤差5%以内に収まっているということである。

発生させた比率のデータのヒストグラムは以下のようになる。

f:id:nekomosyakushimo:20171001122622p:plain

正規分布っぽい形をしたこのヒストグラムの両裾が全体の内の5%になる点が、0.45と0.55であるということになろう。

誤差の水準およびサンプルサイズを変えて同様の実験を繰り返し、ヒストグラムを出力したものが次のもの。

f:id:nekomosyakushimo:20171001123034p:plain

f:id:nekomosyakushimo:20171001123045p:plain

x軸の値に注目していただければ分かるが、サンプルサイズが大きくなるほど、誤差が小さくなっていくのが分かる。それぞれ何パーセントぐらいの結果が、指定した誤差の水準を上回っているかを確認すると、

dat <- rbinom(10000,1068,0.5)
dat <- dat/1068
over <- length(dat[dat > 0.53]) + length(dat[dat < 0.47])
over/10000
# [1] 0.0474

dat <- rbinom(10000,9604,0.5)
dat <- dat/9604
over <- length(dat[dat > 0.51]) + length(dat[dat < 0.49])
over/10000
# [1] 0.0515

やはり5%ぐらいが指定した誤差の水準をオーバーしている。そういう風に値を設定したのだから当然と言えば当然だろうか。微妙に5%とはずれているが、シミュレーションの数を限りなく大きくしていけばこの値は5%に限りなく近づいていくだろう（正確には5%より少しだけ小さい値。サンプルサイズの計算のとき小数点以下を切り上げたから）。

そういえばサンプルサイズの決め方に焦点を当てた本というのがいくつか出ていて、読もう読もうと思って放置していたのだけれど今度どれかを読んでみようかしら。