# Empirical power considering Sterne et al approach 

library(metafor)

M<-45 # M the number of meta-analysis
nPlus<-3 # nPlus the mean number of trials with the characteristic of interest
nMinus<-5 # nMinus the mean number of trials without the characteristic of interest
beta<-log(0.8) # overall impact term
phi<-0.04 # phi^2 the mean between meta-analysis variance
kappa<-0.16 # kappa^2 the mean between-trial variance
sigma<-0.5 # sigma^2 mean within-study variance of the treatment effect estimates

decision=matrix(NA,nrow=2000,ncol=1)




nsim<-0

while (nsim <2000){
  nsim<-nsim+1
  index<-0
  K<-0
  # generate meta-epidemiological datasets
base <- data.frame(matrix(nrow=0, ncol=5)) 
names(base)<-c("ma", "trial", "x", "y.nm", "sigma.nm")

  for (m in 1:M){
    check<-0
    # condition to have more than 2 trials
    # at least 1 with the characteristic 
    # at least 1 without the characteristic
    while (check==0){ 
      # number of trials with the characteristic, median=3
      n.plus<-round(rlnorm(1, meanlog = log(nPlus)-1.056/2, sdlog = sqrt(1.056)),0)
      # number of trials without the characteristic, median=5 
      n.minus<-round(rlnorm(1, meanlog = log(nMinus)-1.056/2, sdlog = sqrt(1.056)),0)

     n.m<-n.plus+n.minus # number of trials in meta-analysis m
      if (n.m>2 & n.plus!=0 & n.minus!=0) check=1
    }

    # between-trial variance of the effect estimate
    tau2<-rlnorm(1,-2.56,sqrt(1.74))
    phi2<-phi*rnorm(1,0,1)  

   for (n in 1:n.m){
      index<-index+1
      temp1<-m
      temp2<-n
      temp3<-1*(temp2>n.minus) # n.minus trials without characteristics
      temp5<-sqrt(runif(1,min=0.25*sigma^2,max=1.75*sigma^2))
      temp4<-(beta+phi2+kappa*rnorm(1,0,1) )*temp3+sqrt(tau2+temp5^2)*rnorm(1,0,1)
      base[index,]<-c(temp1,temp2,temp3,temp4,temp5)
    }

    K<-K+n.m
  }  

  # two-stage model
  betas<-rep(NA,M)  
  sebetas<-rep(NA,M)

  # first stage: estimation of beta.m and its variance, for each meta-analyses m
  for (m in 1:M){
    res<-rma(yi=y.nm, sei=sigma.nm, mods = ~ x, data=base[base$ma==m,],method="DL")
betas[m]<-res$b[2]
    sebetas[m]<-res$se[2]
}

  dat<-cbind(betas,sebetas)

  # second stage: estimation of beta and its variance    
  res2<-rma(yi=betas, sei=sebetas, data=dat,method="DL")
  pval<-res2$pval
  decision[nsim]<-ifelse(pval<0.05,1,0)
}


colMeans(decision)