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 (my Code is similar to Simulation of random-effects meta-analysis yields biased tau^2 up to some point). Also, the rma-function from the metafor-package returns different results than metagen, but the solution (https://stackoverflow.com/questions/41962137/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 (to keep it simple it is just 30 times here)
for (r in 1:30) {
# 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)
#here you can see my problem
print(resultsDfin); print(resultsHfin)
}
#then i addded some code to sum up the results and their means in a
dataframe, as i did with resultsHfin and resultsDfin
#and finally got to clean results for the next iteration
resultsDfin=data.frame(); resultsHfin=data.frame()
}}}
EDIT: after cross-checking with RevMan I think the problem is caused by the composition. Could the problem be the combination of low study number, relatively low baseline risk in control group and a small true tau^2? So that the within study variance is large and therefore the small tau is "eclipsed"?