# re: binary random-eff meta-analysis: Why is my heterogeneity estimate so much higher using GLMM compared to other methods? (metabin function)

I am trying to conduct a random effects meta-analysis for binary outcome data using the GLMM method in the meta package (metabin function / meta package).

meta package citation:
Balduzzi S, Rücker G, Schwarzer G (2019). “How to perform a meta-analysis with R: a practical tutorial.” Evidence-Based Mental Health, 153–160.https://doi.org/10.1136/ebmental-2019-300117

metafor package citation (rma.glmm called internally): Viechtbauer, W. (2010). Conducting meta-analyses in R with the metafor package.Journal of Statistical Software, 36(3), 1-48. https://doi.org/10.18637/jss.v036.i03

I am perplexed about why the heterogeneity estimate is so much higher for GLMM than what I get using other methods (Mantel-Haenszel or Inverse variance). Does anyone know why this might be? I am not sure if it is a coding error on my part or whether there's another explanation for this. I created a dummy dataset that replicates the pattern of findings from my data ("Data1" below). I also attached my output.

The OR estimates were similar across models, but the heterogeneity for GLMM is much larger (albeit, with a very large confidence interval)

OR estimates:

• GLMM: 1.13
• MH: 1.12
• Inverse-variance: 1.12

I^2 estimates:

• GLMM: 53.6% [95% CI: 5.3%; 77.3%]
• MH: 0.0% [95% CI: 0.0%; 62.4%]
• Inverse-variance: 0.0% [95% CI: 0.0%; 62.4%]

dummy data:

ID <- c(1,2,3,4,5,6,7,8,9,10)
Ne <- c(30,20,15,65,20,10,20,20,55,46)
Ee <- c(2,2,2,21,2,2,3,6,22,10)
Nc <- c(30,9,14,68,10,10,10,15,50,25)
Ec <- c(3,3,1,15,2,0,2,3,19,6)
Data1 <- data.frame(ID,Ne,Ee,Nc,Ec)


GLMM code and output:

m.bin.glmm <- metabin(event.e = Ee,
n.e = Ne,
event.c = Ec,
n.c = Nc,
studlab = ID,
data = Data1,
sm = "OR",
method = "GLMM",
fixed = FALSE,
random = TRUE,
method.tau = "ML",
method.incr = "only0",
model.glmm = "UM.RS",
method.random.ci = TRUE,
title = "test")
summary(m.bin.glmm)

Warning message:
Not possible to fit RE/ME model='UM.RS' with nAGQ > 1 with glmer(). nAGQ automatically set to 1.
> summary(m.bin.glmm)
Review:     test

OR             95%-CI
1  0.6429 [0.0995;   4.1530]
2  0.2222 [0.0297;   1.6646]
3  2.0000 [0.1608;  24.8709]
4  1.6864 [0.7778;   3.6561]
5  0.4444 [0.0528;   3.7383]
6  6.1765 [0.2599; 146.7771]
7  0.7059 [0.0978;   5.0957]
8  1.7143 [0.3510;   8.3725]
9  1.0877 [0.4958;   2.3861]
10 0.8796 [0.2772;   2.7911]

Number of studies combined: k = 10
Number of observations: o = 542
Number of events: e = 126

OR           95%-CI    t p-value
Random effects model 1.1340 [0.7041; 1.8262] 0.60  0.5653

Quantifying heterogeneity:
tau^2 = 0; tau = 0; I^2 = 53.6% [5.3%; 77.3%]; H = 1.47 [1.03; 2.10]

Test of heterogeneity:
Q d.f. p-value             Test
19.41    9  0.0219        Wald-type
10.16    9  0.3381 Likelihood-Ratio

Details on meta-analytical method:
- Mixed-effects logistic regression model (random study effects)
- Maximum-likelihood estimator for tau^2
- Random effects confidence interval based on t-distribution (df = 9)
- Continuity correction of 0.5 in studies with zero cell frequencies
(only used to calculate individual study results)


Mantel-Haenszel code and output:

m.bin_MH <- metabin(event.e = Ee,
n.e = Ne,
event.c = Ec,
n.c = Nc,
studlab = ID,
data = Data1,
sm = "OR",
method = "MH",
fixed = FALSE,
random = TRUE,
method.tau = "DL",
method.incr = "only0",
Q.Cochrane = TRUE,
#  model.glmm = "UM.RS",
method.random.ci = TRUE,
title = "test")
summary(m.bin_MH)

Review:     test

OR             95%-CI %W(random)
1  0.6429 [0.0995;   4.1530]        5.0
2  0.2222 [0.0297;   1.6646]        4.3
3  2.0000 [0.1608;  24.8709]        2.8
4  1.6864 [0.7778;   3.6561]       29.3
5  0.4444 [0.0528;   3.7383]        3.9
6  6.1765 [0.2599; 146.7771]        1.7
7  0.7059 [0.0978;   5.0957]        4.5
8  1.7143 [0.3510;   8.3725]        7.0
9  1.0877 [0.4958;   2.3861]       28.4
10 0.8796 [0.2772;   2.7911]       13.1

Number of studies combined: k = 10
Number of observations: o = 542
Number of events: e = 126

OR           95%-CI    t p-value
Random effects model 1.1210 [0.7413; 1.6953] 0.62  0.5477

Quantifying heterogeneity:
tau^2 = 0; tau = 0; I^2 = 0.0% [0.0%; 62.4%]; H = 1.00 [1.00; 1.63]

Test of heterogeneity:
Q d.f. p-value
6.60    9  0.6788

Details on meta-analytical method:
- Mantel-Haenszel method
- DerSimonian-Laird estimator for tau^2
- Mantel-Haenszel estimator used in calculation of Q and tau^2 (like RevMan 5)
- Hartung-Knapp (HK) adjustment for random effects model (df = 9)
- Continuity correction of 0.5 in studies with zero cell frequencies


Inverse-variance code and output:

m.bin_Inv <- metabin(event.e = Ee,
n.e = Ne,
event.c = Ec,
n.c = Nc,
studlab = ID,
data = Data1,
sm = "OR",
method = "Inverse",
fixed = FALSE,
random = TRUE,
method.tau = "DL",
method.incr = "only0",
#  model.glmm = "UM.RS",
method.random.ci = TRUE,
title = "test")
summary(m.bin_Inv)

Review:     test

OR             95%-CI %W(random)
1  0.6429 [0.0995;   4.1530]        5.0
2  0.2222 [0.0297;   1.6646]        4.3
3  2.0000 [0.1608;  24.8709]        2.8
4  1.6864 [0.7778;   3.6561]       29.3
5  0.4444 [0.0528;   3.7383]        3.9
6  6.1765 [0.2599; 146.7771]        1.7
7  0.7059 [0.0978;   5.0957]        4.5
8  1.7143 [0.3510;   8.3725]        7.0
9  1.0877 [0.4958;   2.3861]       28.4
10 0.8796 [0.2772;   2.7911]       13.1

Number of studies combined: k = 10
Number of observations: o = 542
Number of events: e = 126

OR           95%-CI    t p-value
Random effects model 1.1210 [0.7413; 1.6953] 0.62  0.5477

Quantifying heterogeneity:
tau^2 = 0 [0.0000; 1.1907]; tau = 0 [0.0000; 1.0912]
I^2 = 0.0% [0.0%; 62.4%]; H = 1.00 [1.00; 1.63]

Test of heterogeneity:
Q d.f. p-value
6.60    9  0.6792

Details on meta-analytical method:
- Inverse variance method
- DerSimonian-Laird estimator for tau^2
- Jackson method for confidence interval of tau^2 and tau
- Hartung-Knapp (HK) adjustment for random effects model (df = 9)
- Continuity correction of 0.5 in studies with zero cell frequencies

• metabin is not in the metafor package. Jul 30, 2023 at 13:52
• Yes - good point. Sorry for causing confusion - however, a function from the metafor package is used: "For GLMMs, the rma.glmm function from R package metafor (Viechtbauer, 2010) is called internally" (rdrr.io/cran/meta/man/metabin.html) Aug 1, 2023 at 4:01
• thanks for pointing this out - I edited the attributions to be more accurate. Aug 1, 2023 at 4:14

It appears that the heterogeneity diverges pretty strongly between the three methods in terms of $$I^2$$. However, when you are looking at the estimates of heterogeneity in terms of standard deviation or variance (tau = $$\tau$$ or tau^2 = $$\tau^2$$), all three models come to the same estimate (0). Tau aims to estimate the actual spread/dispersion of parameters in the statistical population that you assess, free of sampling error. $$I^2$$ on the other hand is a relative measure of $$I^2 = \frac{\tau^2}{\tau^2 + \nu} \cdot 100\%$$. Parameter $$\nu$$ describes the expected within-study random sampling variance. Since the $$\tau$$-estimates seem to be (nearly) identical across all three methods, it appears that the estimated $$\nu$$-values differ more strongly. Since $$\tau$$ is probably not exactly zero, but pretty close to it, if $$\nu$$ differs substantially large, huge differences in $$I^2$$ may not be surprising. In terms of actual spread of population parameters, the studies all seem to come to similar conclusions.