1. はじめに
農学系の実験では10個程度の小サンプルサイズの実験をすることが多いです。そこで、n=10の場合の結果について、統計的な解釈について考えてみます。
2. パッケージの準備
Rを使います。また、パッケージ"pwr"をインストールしてください。
(tiduverseはデフォルトです。）
library(tidyverse)
library(pwr)
3. 架空データの生成
温州ミカンの試験を仮定します。
栽培期間中にある処理を実施した結果、果実の糖度があがったかどうかを調査します。
対照区：A
処理区：B
試験区の規模：10果
調査内容：Brix(糖度)
＊対照区も処理区も果実糖度の分散は同じと仮定する。
(真実は処理区が0.5%高くなっているが分析者にはわからない）
目的は、この処理をすることで果実糖度に影響を与えるかどうかを知りたい。
set.seed(123456)
non_treat = round(rnorm(10,10,0.5),2)
treat = round(rnorm(10,10.5,0.5),2)
non_treat
treat
>[1] 10.42 9.86 9.82 10.04 11.13 10.42 10.66 11.25 10.58 9.79
>[1] 10.00 9.94 10.47 11.09 11.03 10.53 10.13 10.97 11.33 10.78
上段が対照区の糖度で下段が処理区の糖度になります。
それぞれ平均を計算します。
mean(non_treat)
mean(treat)
>[1] 10.397
>[1] 10.627
対照区が10.4%、処理区が10.6&になりました。これが統計的に有意な差であるか、t検定をしてみます。
t.test(non_treat,treat)
> Welch Two Sample t-test
>
>data: non_treat and treat
>t = -1.0125, df = 17.911, p-value = 0.3248
>alternative hypothesis: true difference in means is not equal to 0
>95 percent confidence interval:
> -0.7074091 0.2474091
>sample estimates:
>mean of x mean of y
> 10.397 10.627
p-value = 0.3248になったので、有意差はなかった。
よって、対照区と処理区の糖度に違いはなく、同じである...
（このデータの真実は違いがあるので同じではない）
このように「有意差がでない=同じである」と、つい解釈しがちですが、有意差がでた場合は、違いがあると言えますががでない場合の解釈は実は難しい。
この場合に統計的に問題となるのは、サンプルサイズが小さいことによる検出力が小さいことである。
4. t検定と検出力
検出力とは、「帰無仮説H0が正しくないときに、このH0を正しく棄却する確率」のこと。
直観的にわかりにくい表現ですが、要するに有意差がでない場合に、両者が同じである確率が、$\beta$ですので、その逆確率になります。
では、今回の調査の検出力はどうなのでしょうか？
検出力を計算するためには効果量の考え方が必要になります。
効果量とは、対照区と処理区との平均値の違いを標準化した量をいいます。
先ほどのミカンの糖度の平均値の違いは（真実の違い）
10.5 - 10 = 0.5
となりますが、これを標準化するとは、平均0、標準偏差1になるように調整します。
（こうすることによって、数字の大小の違いに影響されなくなる）
そのために、両者の違いを共通する標準偏差で割ります。
m1:コントロールの平均
m2:処理区の平均
s1:コントロールの標準偏差
s2:処理区の標準偏差
Cohens d:効果量(正規化された違いの量)
$$
Cohen's;d = \frac{(m_1-m_2)}{\sqrt{\frac{s_1^2+s_2^2}{2}}}
$$
では、上記の効果量を計算する関数を作成します。
cohens_d = function(m1,m2,s1,s2){
return(abs((m1-m2)/(sqrt((s1^2+s2^2)/2))))
}
母集団の標準偏差は両者ともに0.5％とする
10%と10.5%の違いによるCohens dは？
cohens_d(10,10.5,0.5,0.5)
>[1] 1
効果量が１の違いでした。この場合の検出力を計算します。
pwr.t.test(d=1,n=10,sig.level=0.05,type="two.sample",
alternative="two.sided")
> Two-sample t test power calculation
>
> n = 10
> d = 1
> sig.level = 0.05
> power = 0.5620066
> alternative = two.sided
>
>NOTE: n is number in *each* group
検出力(power)は0.56となりました。有意差がでない場合、帰無仮説を正しく棄却する確率は56%とのことです。
意味がまだわかりにくいので、上記条件で1000回調査を実施したとしたら、その内の何回で有意差がでるかシミュレーションで確かめてみます。
test=rep(NA,1000) #記録用の変数
m1=10
m2=10.5
n=10 #サンプルサイズ
for(i in 1:1000){
a = rnorm(n,m1,0.5) #母集団の標準偏差は0.5%
b = rnorm(n,m2,0.5) #母集団の標準偏差は0.5%
t=t.test(a,b) #t_testを実施
test[i]=t$p.value #p-値を記録する
}
1000回のt検定した結果のp値の分布を見てみましょう。
test = data.frame(value=test)
ggplot(test,aes(x=value)) + geom_histogram(aes(y =..density..)) + geom_vline(xintercept = 0.05,color="red") + theme_light()
赤い垂線が有意水準の閾値です。
1000回中、t検定をした結果、正しく有意差がでた回数は？
sum(test<= 0.05)
>[1] 567
1000回のサンプリング＆t検定を実施した結果、567回（57%）に有意差がでました。
検出力とほぼ一致しました。
検出力とは効果量（今回の場合d=1）の下、検定によって正しく有意差がでる確率とも言い換えられます。
効果量が小さくなるとどうなる？ 1から半分の0.5にしてみます。
pwr.t.test(d=0.5,n=10,sig.level=0.05,type="two.sample",
alternative="two.sided")
> Two-sample t test power calculation
>
> n = 10
> d = 0.5
> sig.level = 0.05
> power = 0.1850957
> alternative = two.sided
>
>NOTE: n is number in *each* group
検出力は0.185と小さくなりました。
この場合のミカンの糖度の差が効果量0.5とは実際はどれぐらい違う？
cohens_d(10,10.25,0.5,0.5)
>[1] 0.5
この場合、糖度0.25%の差が効果量0.5に相当します。
今度は効果量が1.5の場合の糖度は何％の違いになる？
cohens_dd(10,10.75,0.5,0.5)
[1] 1.5
効果量1.5とは0.75%の違いになります。検出力を計算します。
pwr.t.test(d=1.5,n=10,sig.level=0.05,type="two.sample",
alternative="two.sided")
> Two-sample t test power calculation
>
> n = 10
> d = 1.5
> sig.level = 0.05
> power = 0.8869702
> alternative = two.sided
>
>NOTE: n is number in *each* group
検出力は0.886になりました。このように想定する違いによって検出力が異なり、効果量が大きくなるほど検出力が大きくなります。
次に効果量同じ（d=1）でサンプル数を増やした場合はどうなるでしょう。
例えばn=20個にした場合は？
pwr.t.test(d=1,n=20,sig.level=0.05,type="two.sample",
alternative="two.sided")
>Two-sample t test power calculation
>
> n = 20
> d = 1
> sig.level = 0.05
> power = 0.868953
> alternative = two.sided
>
>NOTE: n is number in *each* group
検出力は0.56から0.87と上がりました。サンプル数を増やすと検出力が上がることがわかります。
しかし、実際は想定する効果量はわからないことが多いと思うので、サンプルサイズから効果量を推定してみます。
n=10と固定してdを変化させて、検出力が約0.8になるように効果量を調整すると1.33となった。
（違いがないことを言う場合、検出力は0.8とする場合が多い）
pwr.t.test(d=1.33,n=10,sig.level=0.05,type="two.sample",alternative="two.sided")
> Two-sample t test power calculation
>
> n = 10
> d = 1.33
> sig.level = 0.05
> power = 0.802969
> alternative = two.sided
>
>NOTE: n is number in *each* group
効果量1.3となる条件では糖度の違いは0.665%となります。
cohens_d(10,10.665,0.5,0.5)
>[1] 1.33
すなわち、この場合、10果実のサンプルサイズでは、有意差がでない場合において、0.665%の効果量を仮定して違いがない（同じである）との判断が統計的に言えます。
しかしこの場合、真実の糖度の違いは0.5％なので、0.67％の違いまで同じとすると、このサンプルサイズでは、真実の結果を見落とす可能性があり、不適切だと言えるでしょう。
もう少しサンプルサイズを増やす必要があります。
5. 終わりに
10個程度の小サンプルサイズで統計的に言えることは、効果量の仮定が重要で、効果量を1と仮定すると、その検出力は56%になるということでした。検出力が低くなるということは、有意差がたまたま出ればOKですが、出なければ解釈が曖昧なままになる可能性が高くなる（検出力が低いから実験で有意差が出にくい）ということになります。
以上で、間違っているところ等ありましたら、ご指摘をいただけたらありがたいです。