対応のあるt検定をシミュレーションする

MikuHatsune2011-09-18

遺伝学コースの、ただRだけをひたすらやる集中講義
やりたい検定、したい解析を見つける。
そのためにデータセットをシミュレーションして作る。
ということで、対応のあるt検定をやろう。
 
とりあえず、t検定にたどり着くまでの
とてもちいさな、とてもおおきな、とてもたいせつな、あいとゆうきのおとぎばなし
 
データの確認。
1群作って、それに対する治療効果をベクトルに収める。
 
シミュレーションデータの作成。
変化を適当に作成する。
 
という訳で中身はできた。あとは、対応のあるt検定をしてくれる関数を探そう。
その1
その2
 
帰無仮説としては、変化の平均は0、ということになろうか。
 
その2のリンクには
 
t 検定が意味を持ち、最も威力を発揮する状況は、「2群の分布が共に正規分布に近い」「分散がほぼ等しく、分布の違いが中心位置の違いに帰着できる」という状況である。t 検定は外れ値があるようなデータには弱いため、そのような場合は外れ値に強い(ロバストな)検定であるウィルコクソンの順位和検定を用いるのが良い。
 
とある。ウィルコクソンの順位和検定というのも知ってしまったので使ってみる。
マン・ホイットニーのU検定というのもあったが今回はやらない。
 
あとは、このシミュレーションを複数回行い、p値が一様に分布しているか確かめよう。

n1<- 25

Niter<- 3000
result_p_t<- rep(0,Niter)
result_p_w<- rep(0,Niter)

for(i in 1:Niter){
g1<- sample(150:200,size=n1,replace=TRUE)
delta<- sample(-50:50,size=n1,replace=TRUE)
g2<- g1 + delta
result_p_t[i]<- t.test(g1,g2,paired=TRUE)$p.value
result_p_w[i]<- wilcox.test(g1,g2,paired=TRUE)$p.value
}
par(mfcol=c(1,2))
plot(sort(result_p_t))
plot(sort(result_p_w))

par(mfcol=c(1,2))
plot(result_p_t,result_p_w,frame=FALSE,cex=0.5)
lim<- c(0,0.1)
plot(result_p_t,result_p_w,frame=FALSE,cex=0.5,xlim=lim,ylim=lim)

coplotして2つの手法のずれを描いた。