Random-effects-meta-analysis-simulation: zero-estimates for tau^2

Keywords: r statistics simulation random-effects

Question: 

I am working on simulating a random-effects-model for comparison of the DerSimonian-Laird-method vs. Hartung-Knapp-Sidik-Jonkman-method in R. To do so, I chose different combinations of mu (true treatment effect), tau^2 (between-study-variance) and p1 (baseline risk of group 1). For in total 8 combinations of those parameters I want to run 2000 sets of meta-analysis. The problem is: in 80-90% percent of the cases the estimator for tau^2 is equal to zero and I just cannot see why. I am not quite sure whether my implementation of the 2 methods using the meta-package is correct.

My code might be ineffective and not exactly elegant but I hope it is good to read and for my purposes it is fast enough. Also, I think that the code is not so much of a problem, since the problem above seems to originate right after the implementation of the metagen-function. Also, the rma-function from the metafor-package returns different results than metagen, but the solution (Differences between results from meta and metafor packages in R) does not help me. Any help would be much appreciated as I am completely stuck with this one.

# set szenario
k = c(5)  #number of studies
N = c(40, 40, 40, 40, 40) # total N for each study

#define dataframe for results
resultsHfin=data.frame()
resultsDfin=data.frame()

#alter mu (real risk difference), ttau (real between-study-variance) and 
p1 (baseline-risk of control group) 
for (i in 0:1){
mu = (i*2/3)
for (j in 1:2){
ttau = (j*0.05)
for (t in 1:2) {
  p1 = ((t^2)*0.05)

#repeat 1000 times
for (r in 1:1000) {

# simulate sample effect for each study
mui = rnorm(n=k, mean=mu, sd=sqrt(ttau))

##### Simulate Group 1 ##### 
# assume equal numbers in each treatment arm (so n1 is N/2)
n1 = floor(N/2)

# simulate a, number of successes in group 1

a = rbinom(n=k, size=n1, prob=p1)

##### Simulate Group 2 ##### 
# calculate n for this group
n2 = N - n1

# figure out p(success) in this group using log odds
p2 = p1/(p1-exp(mui)*p1+exp(mui))

# simulate c, number of successes in group 2
c = rbinom(n=k, size=n2, prob=p2)

##### Do Analysis #####
require(metafor)
require(meta)

# calculate descriptive statistics
desc = escalc(measure="OR", ai=a, bi=n1-a, ci=c, di=n2-c, to = "if0all", 
add = 1/2, drop00 = FALSE)

# get observed treatment effect and estimated within-study variances for 
each study
ote = desc$yi
ewsv = desc$vi

# compute DSL-method
DSL = rma(yi = ote, vi = ewsv, method = "DL", measure = "OR", to = 
"if0all", add = 1/2, drop00 = FALSE)
DSL1 = metagen(TE = ote, seTE = ewsv, comb.fixed = FALSE, comb.random = 
TRUE, 
method.tau = "DL", hakn = FALSE, prediction=TRUE, sm="OR")

# compute HKSJ-method
HKSJ = metagen(TE = ote, seTE = ewsv, comb.fixed = FALSE, comb.random = 
TRUE, method.tau = "SJ", hakn = TRUE, prediction=TRUE, sm="OR")

#extract tau, est. mu and CI from HKSJ
resultsHKSJ = data.frame(HKSJ$TE.random, HKSJ$seTE.random, 
HKSJ$lower.random, 
HKSJ$upper.random, (HKSJ$tau)^2, HKSJ$I2, HKSJ$pval.random)

#extract tau, est. mu and CI from DSL

resultsDSL = data.frame(DSL$b[1], DSL$se, DSL$ci.lb, DSL$ci.ub, DSL$tau2, 
DSL$I2, DSL$pval)

#add results

resultsHfin=rbind(resultsHfin, resultsHKSJ)

resultsDfin=rbind(resultsDfin, resultsDSL)

}

#then i addded some code to sum up the results and their means in a 
dataframe, as i did with resultsHfin and resultsDfin 
}}}

Answers: