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

Question

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):
print(RE$tau2); print(.V)
# LOOKS COMPLETELY HORRIBLE. WHY IS TAU^2 SO SMALL??
# ITS SE IS SMALL, TOO, SO THERE IS NO EXCUSE?

Edit: A sort-of answer

The culprit is the distribution of the sample sizes (which affect the within-study variances). The above code produces a highly right-skewed distribution of $n$. By generating large samples sizes from a more symmetric distribution, estimation of $\tau^{2}$ is unbiased.

Generating huge sample sizes from a Normal works:

N = round( rnorm(.k, mean=10000, sd=800)   )

Or even from a uniform:

N = round( runif(.k, 1000, 3000) )

But not generating small sample sizes from a uniform:

N = round( runif(.k, 20, 200) )

I'd still very much appreciate a theoretical explanation for this behavior. As far as I know, the random-effects model does not make any assumptions on the within-study variances, so it's disturbing that this seems to matter so much.

Answer 1

The problem is that you are simulating data where the risks/probabilities and the group sizes are both relatively low. As a result, the sampling distributions of the log relative risk are not very well approximated by normal distributions and the sampling variances are poorly estimated. Take a look at:

sum(d0 == 0 | d1 == 0)

You will find that 0 events in at least one of the two groups are relatively common (in more than 15% of the cases). Things start to break down under the 'normal-normal' model under such conditions.

Try making the risks not quite so low and the studies a bit larger:

p0=0.18  # P(success) in control group
N = round( rchisq(.k, df=2) * 30 + 100 )

Things start to look much better then.

If you find yourself in a situation where both risks and sample sizes are low, you may have to switch to a 'binomial-normal' model. For relative risks, this is a bit more difficult, since you would have to fit a mixed-effects logistic regression model with a log link, which is a bit more tricky than the usual logit link (which is used to model odds ratios). For 'binomial-normal' models with a logit link, take a look at the rma.glmm() function. For example, if you are patient, you could try:

rma.glmm(measure="OR", ai=d1, bi=n1-d1, ci=d0, di=n0-d0, model="CM.EL", verbose=TRUE)

Note that you won't be estimating $\tau^2 = .30$ then, since that applies to the (log) relative risks, not the (log) odds ratios. But when risks are low, then relative risks and odds ratios are similar, so you will get something not too far off.