今までの記事と全く違う分野ですが…。
統計学を専門としている研究者が新しい統計解析手法を開発するような場合、想定した状況で正しい結果が求まることを示すためにシミュレーション研究を行います。しかしながらシミュレーションはそのような特殊な状況でのみ使うようなものではなく、他の分野でも様々に用いられますし、普段の解析の際にも解析方針を決めるにも役立ちます。この記事ではシミュレーションの基礎について、R言語を用いて簡単に説明したいと思います。
t検定の状況設定
統計学の教科書で、平均や分散などの記述統計の後、検定として最初に登場するのが2群の平均値の差があるかどうかをデータから判断するt検定ではないでしょうか。ここでは、t検定をシミュレーションで実行してみたいと思います。
身近な医学の題材を選び、2つの集団（A群とB群）の収縮期血圧に差があるかどうかを調べる状況を考えます。実際は真の分布は分からないのですが、ここでは2群の血圧が、
- A群は平均が120、標準偏差20の正規分布
- B群は平均が140、標準偏差20の正規分布
に従うと仮定します。現実にここまで分かっていれば、そもそもサンプルを取ってきて検定をする必要すらなく、2群の平均値の差は異なる、と結論付ければいいわけです。しかしながら、実際には真の分布は分からないために、手元のサンプルから判断するわけです。R言語の基礎については、機会があれば書きます（よってここでは省略します）。
1回のt検定を実施
R言語にて、正規分布からの乱数生成はrnorm関数で実行できます。それぞれから10サンプルを取ってきたとします。
dat_a <- rnorm(10, 120, 20) dat_b <- rnorm(10, 140, 20)
と実行すると、それぞれdat_aとdat_bに平均が120と140、標準偏差20の10人分のサンプルデータができます。実行ごとに結果は違うと思いますが、手元だとdat_aが
134.8699 115.3814 123.0973 138.3278 102.4526 108.6041 162.6671 160.2167 122.3463 133.6674
dat_bが
123.2808 134.9284 125.6528 162.6626 153.7979 149.5983 150.5219 131.6545 160.5555 220.7252
のようになりました。さて、各群10人分、合計20人分のデータのみが手元に得られたとして、この分布の平均は異なると結論付けることができるのか、さっそくt検定を実行してみましょう。
t.test(dat_a, dat_b, var.equal=TRUE)
と実行すると、
Two Sample t-test data: dat_a and dat_b t = -1.9362, df = 18, p-value = 0.06871 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -44.151005 1.801542 sample estimates: mean of x mean of y 130.1631 151.3378
という結果が得られます。（var.equal=TRUEの指定にはここでは深入りしませんが、分布を見てある程度ばらつきが異ならないようなら付けますし、かなり違う状況ではそもそもt検定を実施すること自体が危ぶまれると思われます）真の分布は異なるにも関わらず、手元のサンプルでは差があるとは言えないということです。データはばらつきますので、実際にはよくあることです。
では、この分布から同じ設定で10人を取り出して検定を繰り返すとどうなるでしょうか。普通研究は1回きりですが、仮に何回も行ったとするシミュレーションを実施してみます。t.test関数でt検定を実行すると、その結果がオブジェクトとして返されますが、その中でp値はp.valueに入っているので、以下のようにすると取り出すことができます。
t.test(dat_a, dat_b, var.equal=TRUE)$p.value
今回では0.06870732が得られました。
シミュレーションの開始
ここからが反復を行います。ここでは繰り返しを基本的なforループで書きたいと思います（そのうち機会があればapply系関数の使い方も記事にします）。ここの流れとしては、まず結果であるp値を格納する変数を作り、データ発生からt検定の実施、さらに結果の格納までをループ内に収めます。
result <- numeric(1000) #結果格納用変数
for(i in 1:1000){
dat_a <- rnorm(10, 120, 20)
dat_b <- rnorm(10, 140, 20)
result[i] <- t.test(dat_a, dat_b, var.equal=TRUE)$p.value #結果格納
}
これを実行すれば、resultに1000回同じことを実行した際のp値が格納されます。手元だと一瞬で終わりました。p値が0.05よりも小さく有意であった回数を数えてみます。あ、もちろん手でじゃなくて、以下を実行すればすぐに分かります。
sum(result < 0.05)/1000
手元では0.56となりました。これが検出力です。要は真の分布の平均は120と140で20も差があるんだけど、10人ずつのサンプルでは、個人のばらつきに比べた相対的な差があまり大きくなく、50%くらいしか検定では差があると結論付けられないことになります。
では分布を変えずに、どうすれば差があるとより結論付けられるようになるか、と言えば、サンプルサイズを増やすことです。これも、簡単に試すことができます。例えば、20人ずつに増やすのであれば、rnorm関数の最初の引数を10から20に変更するだけです。やってみますと、手元では0.854となり、80%を越える確率で有意な差があると結論付けられます。
このように分布（真の差やばらつきなど）を定め、サンプルサイズを変更すると、有意な結果が得られる確率、つまり検出力が変わってきます。逆にある一定以上の確率で有意な結果が得られるだけのサンプルサイズを求めることを、サンプルサイズ設計といいます。サンプルサイズ設計は、単純な状況では解析的に式で求まりますが（今回の状況も含む）、以上のようにシミュレーションで行うことも可能です。
このようにシミュレーションを行うことで、仮定の元でどのような結果が得られるか、ということを試みることができます。今回のように個々のデータの標準偏差が20で、両群の平均値の差が20くらいの状況でt検定にて有意な結果を得るのであれば、各群のデータは20くらいは欲しいということが分かります。
コメント