Power inflation index

powerinflationindex

## true effect size:
D<-seq(1,250,by=1)
## SE from a study:
se<-46
## sample size
n<-37
## typical SD:
stddev<-se*sqrt(n)

## rejection value (absolute) under null:
qnull<-abs(qnorm(0.025,mean=0,sd=se))

## sample repeatedly to get an estimated d:
nsim<-10000

pii<-true_power<-rep(NA,length(D))
pii_CI<-matrix(rep(NA,length(D)*2),ncol=2)

## for each D:
for(j in 1:length(D)){
  currentD<-D[j]
  ## repeatedly sample to compute obs power:
  drep<-rep(NA,nsim)
  for(i in 1:nsim){
     drep[i]<-mean(rnorm(n,mean=currentD,sd=stddev))
  }

  ## actual power for currentD under repeated sampling:
  true_power[j]<-pow<-mean(ifelse(abs(drep/se)>2,1,0))
  ## observed power from *each* repeated study:
  obspower<-pnorm(qnull,mean=drep,sd=se,lower.tail=FALSE)+
  pnorm(-qnull,mean=drep,sd=se,lower.tail=TRUE)
  ## power inflation index:
  pii[j]<-mean(obspower/pow)
  ## uncertainty bounds of PII:
  pii_CI[j,]<-quantile(obspower/pow,probs=c(0.025,0.975))
}

plot(true_power,pii,type="p",ylim=c(0,14),ylab="power inflation index",
     xlab="true power")
arrows(x0=true_power,x1=true_power,y0=pii_CI[,1],y1=pii_CI[,2],code=3,angle=90,length=0.01,col="gray")