理学療法士がまとめたノート

統計学備忘録（R言語のメモ）

since2016　ときどきTEXのメモ

符号検定

ﾉﾝﾊﾟﾗﾒﾄﾘｯｸ 確率

符号検定 ( sign test )　一標本問題
n個のうちｘ個以下の符号が他の符号と異なる確率を求める検定です．正と負が同じ確率で出るときに、正（または負）の値の数がxとなる確率が B ( n , 0.5 ) の2項分布に従うことを利用します．
例
y<- c (1.5, 1.1, 1.2, -1.2, 2.3)

yの中に含まれる負の値をもつ数は1個です．少なくとも4個は正になる確率を求めます． 5個中1個負になるか、または5個中全てが正になる確率 です．また例題では一つが負なのですが、一つが正になる場合も同じ確率になりるので両側確率を求めることになります．以下のグラフを見てみましょう．

x <- 0 : 5　                              #ｘ軸は成功回数
y <- dbinom ( x, 5, 0.5 )　    #ｙ軸は確率密度
plot ( x, y, type = "h", col = c(2, 2, 1, 1, 2, 2), lwd = 3 )
求める確率は以下の赤部分の合計

f:id:yoshida931:20170711111556p:plain

 

5個が全て正になる場合は  5C0 = choose( 5 , 0 ) = 1
5個が全て正になる確率は　( 1/2 )^5 = 1/32
5個中1個が負になる場合は  5C1 = choose( 5, 1 ) = 5 通り
5個中1個が負になる確率は　( 1/32 ) *5 = 5/32
5個中2個が負になる場合は  5C2 = choose( 5, 2 ) = 10 通り
5個中2個が負になる確率は　( 1/32 ) *10 = 10/32

上記のように考えてｐ値を算出します．
yは5個中1個が負の数になっています．したがって5個のうち1個以下の符号が他と異なる確率を求めることになります．つまり少なくとも4個が同じ符号である確率を求めます．

1個が正になる確率＋全てが負になる確率＝ 5/32 + 1/32
1個が負になる確率＋全てが正になる確率＝ 5/32 + 1/32
ｐ値 =  ( 5/32 + 1/32 ) * 2 =  0.375

Rの関数で確かめます

2*pbinom ( 1, 5, 0.5 )
=0.375

上式は以下のようになります

5個が全て正になる確率
= choose ( 5,0 ) * ( 0.5^0 ) * ( 0.5^5 )　　#成功回数0回
= dbinom ( 0, 5, 0.5 )
= 0.03125

5個中1個が負になる確率
= choose ( 5,1 ) * ( 0.5^1 ) * ( 0.5^4 )　　#成功回数1回
= dbinom ( 1 , 5 , 0.5 )
= 0.15625

両側確率なので×2
{ dbinom ( 0 , 5 , 0.5 ) +  dbinom  ( 1 , 5 , 0.5 )  } *2 
= 0.375

 

nが大きい場合
E[x] = np  ,  V[x] = np ( 1 – p ) より　
p = 1/2 を代入して正規近似にて検定を行います．

なぜｔ検定を繰り返してはいけないのか？

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)

f:id:yoshida931:20170709103144p:plain

視覚的にはAD間、BD間には差があるように見えます

以下、重要な前提です
Aの母平均をμa
Bの母平均をμb
Cの母平均をμc
Dの母平均をμd　とします
「μa=μb=μc=μd」を検定する場合にｔ検定を繰り返してはいけないのです．

 

一元配置分散分析

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」である．

ｔ検定の帰無仮説
有意水準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
２種類の既存薬ＡとＢを組み合わせた配合薬Ｃの配合効果を評価する場合．この場合、既存薬Ａ，Ｂのそれぞれの効果と配合薬Ｃの効果を比較します．ここで いいたいのはＣが既存薬Ａ、Ｂの両方よりも効果があるということです．「ＣがＡよ りも優れている、かつＣがＢよりも優れている」ということを 示します．つまり帰無仮説はC=AかつC=Bとなり、A=B=Cとは異なります．このような場合には２標本ｔ検定を繰り返して用いるほうが適切だと考えられます．

Rstudioでのデータ保存、読み込み

R言語 RStudio

よく使うので書き留めておきます

事前準備
Rstudio 四つの画面から編成されていますので、勝手に名前をつけます．
①書き込み画面、②Consol画面、③環境画面、④ファイル画面

f:id:yoshida931:20170620165251p:plain


まずRスクリプト上に書いているベクトルを保存してみます
①にデータを書き込みます
t3.4<-c(150,132,144,139,118,135,123,133,152,136)

f:id:yoshida931:20170620164100p:plain

　　出典）出典）柳川 堯 , 荒木 由布子; バイオ統計の基礎―医薬統計入門,近代科学社 ,2010, 表3.12

 

保存するために①に
save ( 保存するデータ , file = " 保存後のファイル名 .RData" )

save (  t3.12 , file = " 表3.12.RData " )

f:id:yoshida931:20170620164424p:plain

 

④にデータファイルができてます

f:id:yoshida931:20170620164601p:plain

 

完成したファイルをクリックすると③にデータを読み込みます．

f:id:yoshida931:20170620164731p:plain

 最後に①にt3.12と入力するだけでベクトルになります．
データファイルの保管場所は、Rprojectのフォルダになります
もちろんデータフレームも同じ要領で保存できます．

x<-c(LETTERS[1:10])
y<-t3.12
z<-data.frame(x,y)
save(z,file = "表3.12名前入り.RData")

f:id:yoshida931:20170620170456p:plain

Zの部分をクリックすると①に次の表が出力されます

f:id:yoshida931:20170620170607p:plain

Zという一つのアルファベットの中に上記のデータセットが格納されているのです．
便利です！