臨床統計もおもしろいですよ、その2

これで、
# 2x2のクロス表からstanを使ってMCMCサンプリングしてFisherテストしてp.valueの配列からその分布を表示させる

library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

contingency2stan2pf <- function(mat,alpha=0.05,print=T){
N=apply(mat,1,sum)
data <- list(mat=mat,N=N)
stanModel=readRDS('contingency.rds')
fit=sampling(stanModel,data=data, chain=1, iter=25000,warmup=5000)
pars=c('y1','n1_y1','y0','n0_y0','d','RR','OR','NNT')
if(print) print(fit, pars=pars,digits=4,prob=c(0.025,0.5,0.975))
ms=rstan::extract(fit)
p.sim=NULL
for(i in 1:length(p11)){
MAT=round(with(ms,matrix(c(y1[i],n1_y1[i],y0[i],n0_y0[i]),ncol=2,byrow=TRUE)))
p.sim=append(p.sim,fisher.test(MAT)$p.value)
}
BEST::plotPost(p.sim,compVal = 0.05, xlab='p.value')
cat('Pr(p.value < ',alpha,') = ', mean(p.sim<alpha),'\n')
cat('Original p.value = ',fisher.test(mat)$p.value,'\n')
print(summary(p.sim))
invisible(p.sim)
}

サンプル数を増やしてｐを小さくして信用させるという手法を見破れるかな。