#まず、pの小さめなデータを作る
Niter<-1000
library(Rassoc)
st<-rep(0,Niter)
p<-rep(0,Niter)
for(i in 1:Niter){
af<-runif(1)*0.6+0.2
delta<-rnorm(1)
af1<-af+af*0.05*delta
af2<-af-af*0.05*delta
case<-sample(c(0,1,2),1000,c(af1^2,2*af1*(1-af1),(1-af1)^2),replace=TRUE)
cont<-sample(c(0,1,2),1000,c(af2^2,2*af2*(1-af2),(1-af2)^2),replace=TRUE)
t<-matrix(c(length(case[case==0]),length(case[case==1]),length(case[case==2]),
length(cont[cont==0]),length(cont[cont==1]),length(cont[cont==2])),nrow=2,byrow=TRUE)
cattout<-CATT(t)
st[i]<-(cattout$statistic)^2
p[i]<-cattout$p
}
#"p.adjust()"関数でFDR
pfdr<-p.adjust(p,"fdr")
peven<-ppoints(length(p),a=0)
a=0.05
plot(peven,sort(peven,decreasing=TRUE),xlim=c(0,1),ylim=c(0,1),type="l")
par(new=T)
#この線がFDR(BH法)の基準線。これより下なら採択。
plot(peven,a*((length(p)-(1:length(p))+1)/length(p)),xlim=c(0,1),ylim=c(0,1),type="l",col="red")
par(new=T)
plot(peven,sort(pfdr,decreasing=TRUE),xlim=c(0,1),ylim=c(0,1),type="l",col="blue")
abline(h=a)
#p.adjust()関数の返り値は基準直線としてaをいくつに設定したら交わるか、の値を返す。
#水平直線が0.05の高さになっていることがわかる
abline(v=peven[length(which(pfdr>a))])
