Show that FE model has nominal coverage for conditional inferences.

fe_coverage.r

rm(list=ls())

library(metafor)

iters <- 1000

### number of studies/effects
k <- 10

### SD of the true effects
tau <- 0.2

cov.w <- rep(NA, iters)
cov.u <- rep(NA, iters)

for (j in 1:iters) {

   ### simulate estimates with heterogeneous true effects
   thetai <- rnorm(k, 0, tau)
   vi <- runif(k, .001, 1)
   yi <- rnorm(k, thetai, sqrt(vi))

   ### estimate weighted mean of thetai values
   res <- rma(yi, vi, method="FE")

   ### check for coverage of the weighted mean
   wi <- 1/vi
   theta.mean <- weighted.mean(thetai, w=wi)
   cov.w[j] <- (res$ci.lb <= theta.mean) && (res$ci.ub >= theta.mean)

   ### estimate unweighted mean of thetai values
   res <- rma(yi, vi, method="FE", weighted=FALSE)

   ### check for coverage of the unweighted mean
   theta.mean <- mean(thetai)
   cov.u[j] <- (res$ci.lb <= theta.mean) && (res$ci.ub >= theta.mean)

}

### print coverage rates
print(mean(cov.w))
print(mean(cov.u))