なぜt検定を繰り返してはいけないのか?
なぜt検定を繰り返してはいけないのか
投稿日2017.2.24
修正日2017.7.9
sample は「石村卓夫;分散分析のはなし,東京図書,1992,p137」
A<-c(12,14,16)
B<-c(13,15,17)
C<-c(15,17,19)
D<-c(17,19,21)
x<-c(rep("A",3),rep("B",3),rep("C",3),rep("D",3))
y<-c(A,B,C,D)
stripchart(y~x,vertical=T,pch=1,cex=2,cex.axis=1.5)
視覚的にはAD間、BD間には差があるように見えます
以下、重要な前提です
Aの母平均をμa
Bの母平均をμb
Cの母平均をμc
Dの母平均をμd とします
「μa=μb=μc=μd」を検定する場合にt検定を繰り返してはいけないのです.
一元配置分散分析
oneway.test(y~x,var.equal = T)
#F = 3.6875, num df = 3, denom df = 8, p-value = 0.06214
#有意水準5%で差が「ない」という結果になりました
#全ての群の母分散は等しいと仮定してt検定を繰り返してみます
AB<- t.test ( A , B , paired=FALSE , var.equal=T , conf.level=0.95 )
AC<- t.test ( A , C , paired=FALSE , var.equal=T , conf.level=0.95 )
AD<- t.test ( A , D , paired=FALSE , var.equal=T , conf.level=0.95 )
BC<- t.test ( B , C , paired=FALSE , var.equal=T , conf.level=0.95 )
BD<- t.test ( B , D , paired=FALSE , var.equal=T , conf.level=0.95 )
CD<- t.test ( C , D , paired=FALSE , var.equal=T , conf.level=0.95 )
#それぞれの検定結果のP値のみを取り出してみます
str ( AB [ [3] ] ) ; str ( AC [ [3] ] ) ; str ( AD [ [3] ] ) ; str ( BC [ [3] ] ) ; str ( BD [ [3] ] ) ; str ( CD [ [3] ] )
AB 0.573
AC 0.14
AD 0.0376
BC 0.288
BD 0.0705
CD 0.288
やはりAD間には有意水準5%で「差がある」という結果になりました.しかし分散分析では差がないという結果でした.
分散分析の帰無仮説
有意水準0.5において「μa=μb=μc=μd」である.
t検定の帰無仮説
有意水準0.5において「μa=μb」かつ「μa =μc」かつ「μa =μd」かつ「μb=μc」かつ「μb=μd」かつ「μc=μd」である.言いかえれば、有意水準0.5で「6組全てに差が見られない」となり、対立仮説は「少なくとも一組には差がある」となります.
「μa=μb」かつ「μa =μc」かつ「μa =μd」かつ「μb=μc」かつ「μb=μd」かつ「μc=μd」となる確率は ( 1 - 0.05 ) ^6 です.したがって、有意水準0.5と設定した検定で少なくとも一組には差がある確率が 1 - ( 1 - 0.05 ) ^6 = 0.2649081 となり、有意水準26.5%の検定となってしまいます.つまりαエラーの危険性が高くなるのです.このような矛盾を解消するための手法が多重比較です.
多重比較(Tukey)を行ってみます
TukeyHSD ( aov (y ~ x) , ordered = TRUE )
Tukey multiple comparisons of means
95% family-wise confidence level
factor levels have been ordered
Fit: aov(formula = y ~ x)
$x
diff lwr upr p adj
B-A 1 -4.2294189 6.229419 0.9252929
C-A 3 -2.2294189 8.229419 0.3245304
D-A 5 -0.2294189 10.229419 0.0609499
C-B 2 -3.2294189 7.229419 0.6297636
D-B 4 -1.2294189 9.229419 0.1441838
D-C 2 -3.2294189 7.229419 0.6297636
有意水準0.05ではどのペアにも差はありませんでした.
T検定でよい?・・・研究目的が重要です
A、B、Cを比較する場合は、単純に多重比較と決めつけてしまうのではなく、求めようとすることを明確にすることが大切です。例えば、「A>BかつA>B」を証明するための研究計画であれば、AvsBとAvsBのt検定を行い、両方ともに有意差が言えればよいということになります。
参考 http://www.gen-info.osaka-u.ac.jp/MEPHAS/ave.html
2種類の既存薬AとBを組み合わせた配合薬Cの配合効果を評価する場合.この場合、既存薬A,Bのそれぞれの効果と配合薬Cの効果を比較します.ここで いいたいのはCが既存薬A、Bの両方よりも効果があるということです.「CがAよ りも優れている、かつCがBよりも優れている」ということを 示します.つまり帰無仮説はC=AかつC=Bとなり、A=B=Cとは異なります.このような場合には2標本t検定を繰り返して用いるほうが適切だと考えられます.
Rstudioでのデータ保存、読み込み
よく使うので書き留めておきます
事前準備
Rstudio 四つの画面から編成されていますので、勝手に名前をつけます.
①書き込み画面、②Consol画面、③環境画面、④ファイル画面
まずRスクリプト上に書いているベクトルを保存してみます
①にデータを書き込みます
t3.4<-c(150,132,144,139,118,135,123,133,152,136)
出典)出典)柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010, 表3.12
保存するために①に
save ( 保存するデータ , file = " 保存後のファイル名 .RData" )
save ( t3.12 , file = " 表3.12.RData " )
④にデータファイルができてます
完成したファイルをクリックすると③にデータを読み込みます.
最後に①にt3.12と入力するだけでベクトルになります.
データファイルの保管場所は、Rprojectのフォルダになります
もちろんデータフレームも同じ要領で保存できます.
x<-c(LETTERS[1:10])
y<-t3.12
z<-data.frame(x,y)
save(z,file = "表3.12名前入り.RData")
Zの部分をクリックすると①に次の表が出力されます
Zという一つのアルファベットの中に上記のデータセットが格納されているのです.
便利です!
シミュレーションの準備
まずはforの使い方(なかなか慣れません)
変数xに1から5を足す
0+1+2+3+4+5=15をfor関数で実行
x<-0
for(i in 1:5) x<-x+i
ベクトルに繰り返し数を付ける
x<-c()
i<-c(1,2,3,4,5)
x<-c(x,i)
for関数で実行
x<-c()
for(i in 1:5) x<-c(x,i)
関数として定義すると
myfunc01<-function(){
x<-c()
for(i in 1:5) { # for(i in 1:5){}で{}を5回繰り返す
x<-c(x,i)
}
return(x)
}
myfunc01()
次にif関数の基本です(else ifで条件を増やせます)
coin<-function(){
x<-runif(1) #一様乱数を5個準備
if (x<=1/2) shoubu<-1 #勝負!表=1
else shoubu<-0 #勝負!裏=0
return(shoubu)
}
coin() #表が出たら1、裏が出たら0
上記のcoin関数結果から表回数を集計したらシミュレーションができそうです
count<-0
if(shoubu==1) count<-count+1
coin関数を10回繰り返して、表の出た回数を合計してみます
myfunc02<-function(n){
count<-0
for(i in 1:n){ # for(i in 1:n){}で{}をn回繰り返すことになる
x<-runif(1) #一様乱数を5個準備
if (x<=1/2) shoubu<-1 #勝負!表=1
else shoubu<-0 #勝負!裏=0
if(shoubu==1) count<-count+1 #勝負!表=1ならカウントする
}
return(count)
}
myfunc02(10)
先にcoin関数を定義しておくと
coin<-function(){
x<-runif(1) #一様乱数を準備
if (x<=1/2) shoubu<-1 #勝負!表=1
else shoubu<-0 #勝負!裏=0
return(shoubu)
}
myfunc02を最適化できます
myfunc02<-function(n){
count<-0
for(i in 1:n){ # for(i in 1:n){}で{}をn回繰り返すことになる
x<-runif(1) #一様乱数を準備
count<-count+coin() #勝負!表=1ならカウントする(coin関数)
}
return(count)
}
myfunc02(10)
次回はモンテカルロシミュレーションで近似解を得る関数に挑戦したいと思います!