Случайные эффекты-мета-анализ-моделирование: нулевые оценки для тау ^ 2

Question

Я работаю над моделированием модели случайных эффектов для сравнения метода DerSimonian-Laird с методом Hartung-Knapp-Sidik-Jonkman в R. Для этого я выбрал разные комбинации mu (истинное лечение эффект), тау ^ 2 (дисперсия между исследованиями) и р1 (базовый риск для группы 1). Всего для 8 комбинаций этих параметров я хочу запустить 2000 наборов мета-анализа. Проблема в том, что в 80-90% случаев оценка для tau ^ 2 равна нулю, и я просто не понимаю, почему. Я не совсем уверен, правильна ли моя реализация двух методов с использованием метапакета.

Мой код может быть неэффективным и не совсем элегантным, но я надеюсь, что он хорош для чтения, и для моих целей он достаточно быстр. Кроме того, я думаю, что код не является большой проблемой, так как вышеупомянутая проблема, кажется, возникает сразу после реализации функции metagen. Кроме того, функция rma из пакета metafor возвращает результаты, отличные от metagen, но решение ( Различия между результатами из пакетов meta и metafor в R ) мне не помогает. Любая помощь будет высоко ценится, так как я полностью застрял с этим.

# 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()
}}}