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, DSL\$se, DSL\$ci.lb, DSL\$ci.ub, DSL\$tau2,
DSL\$I2, DSL\$pval)

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: