# Simulation of random-effects meta-analysis yields biased tau^2

I need to realistically simulate study effect sizes and within-study variances for a random-effects meta-analysis in which the outcome is a relative risk. My question: why does this simulation lead to severe underestimation of $\tau^{2}$, and how can I fix it?

## Simulation approach

• Simulate total sample sizes for each of $k$ studies. Assume equal numbers in control and treatment groups for all studies.
• Draw a population effect, $\mu_{i}$, for each study from a Normal with a specified mean and variance $V$. (Keep your eye on the variance -- it is the heart of my woes.) This is on the log-RR scale.
• Draw the number of successes in the control group from a binomial with success probability $p_{0}$ (fixed across all studies).
• Draw the number of success in the treatment group from a binomial with success probability $p_{1i} = \exp\{ \mu_{i} \} \cdot p_{0}$.
• Compute observed effect sizes, $y_{i}$, for each study from the above data.
• Fit a random-effects meta-analysis using Viechtbauer's metaforpackage.
• Become upset that $\tau^{2} \ne V$.

## Explanations I've considered

Note that metafor by default estimates $\tau^{2}$ via REML. Using other options, such as Dersimonian-Laird, does not help. Also, the estimated SE of $\tau^{2}$ is also small, so it's not a precision problem. Finally, fiddling with the parameters doesn't help either.

## Reproducible example

# set parameters
.k = c(1000)  # huge number of studies to eliminate concerns about asymptotics
.Mt = log(1.5)  # mean of true effects (log-RR)
.V = 0.3  # variance of true effects
p0=0.08  # P(success) in control group
seed = 131457

# simulate total N for each study
# right-skewed with minimum 20
N = round( rchisq(.k, df=2) * 30 + 20 )

# simulate population effect for each study
Mi = rnorm( n=.k, mean=.Mt, sd=sqrt(.V) )

##### Simulate Control Groups #####
# assume equal numbers in each treatment arm (so denom is N/2)
n0 = floor(N/2)

# simulate d0, number of successes in group 0
d0 = rbinom( n=.k, size=n0, prob=p0 )

##### Simulate Treatment Groups #####
# calculate n for this group
n1 = N - n0

# figure out p(success) in this group using population RR = exp(Mi)
# use pmin to ensure probability isn't above 1
# I used a small p0 and smallish .Mt to ensure that there is very little truncation
p1 = pmin( exp(Mi) * p0, 1 )

# simulate d1, number of successes in group 1
d1 = rbinom( n=.k, size=n1, prob=p1 )

# calculate true ES using metafor
require(metafor)
temp = escalc( measure="RR", ai=d1, bi=n1-d1, ci=d0, di=n0-d0)
#head( log( (d1/n1) / (d0/n0) ) ); head(temp)  # yup, matches manual approach

# get observed effect sizes and within-study variances for each study
th.t = temp$yi vyi = temp$vi

# see if true ES are hitting correct tau^2
RE = rma.uni( yi=th.t, vi=vyi, measure="RR")

# these should be the same (var of true ES):