例題をやってみる その1の続き ブートストラップ法

MikuHatsune2011-11-25


p値の続きのこちらの続き。
 
ブートストラップ法を用いてp値を計算してみようと試みている。
信頼区間のときにやったのと同じようなことである。
 
本では、12人からの体温サンプルを500回リサンプリングしている。
このとき、わずか7つが平均が37.0を超えたらしい。
だから(片側)p値が0.014だと言っている。

L<- (1:10)*1000 # リサンプリングする回数
trials<- 100 # 1リサンプリングあたりのシミュレーション回数
data<- matrix(0,trials,length(L))
# 3分くらいかかります。
for(l in L){
	for(i in 1:trials){
		boot<- matrix(sample(tm,size=l*length(tm),replace=TRUE),nc=length(tm))
		pvalue<- sum(apply(boot,1,mean)>37)/l
		data[i,which(L==l)]<- pvalue
	}
}
matplot(t(data),type="p",pch=1,cex=0.2,col=1,axes=FALSE)
axis(1,1:length(L),L)
axis(2)
lines(1:length(L),apply(data,2,mean),col=4,lwd=3)

0.01に収束している感じ。

# 片側検定で、観測平均が仮の母平均より小さいと既に期待する。
t.test(tm,mu=37,alternative="less",paired=FALSE,conf.level=0.95)

	One Sample t-test

data:  tm 
t = -2.0169, df = 11, p-value = 0.03439
alternative hypothesis: true mean is less than 37 
95 percent confidence interval:
     -Inf 36.97443 
sample estimates:
mean of x 
 36.76667 

p値が0.034とかなり差があるけどいいのかな?