[R-meta] Bivariate generalized linear mixed model with {metafor}

[Viechtbauer, Wolfgang (SP)]

Ok, so here is what I came up with:

library(metafor)
library(lme4)

dat <- structure(list(study.name = c("Cutler1995", "Dahlof1991", "Dowson2002", "Ensink1991", "Geraud2000", "Goadsby1991", "Goadsby2000", "Havanka2000", "Kaniecki2006", "Mathew2003", "Myllyla1998", "Nappi1994", "Patten1991", "Pfaffenrath1998", "Sandrini2002", "Sargent1995", "Sheftell2005", "Tfelt-Hansen1995", "Tfelt-Hansen1998", "Visser1996"), n1 = c(66L, 275L, 194L, 131L, 504L, 89L, 129L, 98L, 131L, 831L, 42L, 142L, 142L, 277L, 170L, 46L, 902L, 122L, 388L, 72L), n2 = c(65L, 182L, 99L, 78L, 56L, 93L, 142L, 91L, 127L, 419L, 41L, 81L, 101L, 91L, 84L, 47L, 892L, 126L, 160L, 85L), r1 = c(37L, 180L, 123L, 60L, 304L, 45L, 63L, 59L, 64L, 471L, 33L, 73L, 95L, 175L, 85L, 26L, 649L, 63L, 239L, 33L), r2 = c(17L, 48L, 42L, 14L, 24L, 9L, 30L, 28L, 47L, 105L, 12L, 25L, 22L, 27L, 25L, 8L, 375L, 30L, 64L, 15L)), row.names = c(NA, 20L), class = "data.frame")

# double-check that we can reproduce OR = 3.50 (95% CI, 2.94 to 4.16) for
# the two-stage random-effects model

res <- rma(measure="OR", ai=r1, n1i=n1, ci=r2, n2i=n2, data=dat, method="DL")
predict(res, transf=exp, digits=2)

# the model being fitted is essentially this one

res <- rma.glmm(measure="OR", ai=r1, n1i=n1, ci=r2, n2i=n2, data=dat, model="UM.RS", cor=TRUE)
summary(res)
predict(res, transf=exp, digits=2)

# we can also fit this model directly with glmer()

dat <- to.long(measure="OR", ai=r1, n1i=n1, ci=r2, n2i=n2, data=dat)

res <- glmer(cbind(out1,out2) ~ group + (group | study), data=dat, family=binomial)
summary(res)
round(exp(c(fixef(res)[2], confint(res, parm="group1", method="Wald"))), 2)

# note the exact same pooled OR and 95% CI

# however, the authors report the marginal OR, so for this we fit the model as

res <- glmer(cbind(out1,out2) ~ 0 + group + (0 + group | study), data=dat, family=binomial)
summary(res)

# and then we can estimate the marginal OR as described in B.2

cval <- 16*sqrt(3) / (15*pi)
p1hat <- 1 / (1 + exp(-fixef(res)[1] / sqrt(1 + cval^2 * VarCorr(res)[[1]][1,1])))
p2hat <- 1 / (1 + exp(-fixef(res)[2] / sqrt(1 + cval^2 * VarCorr(res)[[1]][2,2])))
round((p2hat / (1 - p2hat)) / (p1hat / (1 - p1hat)), 2)

# there is the 3.48 that the authors report; getting the CI for this marginal OR
# requires a bit more work using the delta method

So, to summarize: The model they describe *is* 'Model 6'. You can also fit this with rma.glmm() as shown above (with the 'cor=TRUE' option, which I hadn't implemented yet when we wrote that 2018 paper). I don't know why they focus on the marginal OR, but I didn't read the article fully.

Best,
Wolfgang

>-----Original Message-----
>From: Arthur Albuquerque [mailto:arthurcsirio using gmail.com]
>Sent: Wednesday, 09 March, 2022 20:37
>To: r-sig-meta-analysis using r-project.org; Viechtbauer, Wolfgang (SP)
>Subject: RE: [R-meta] Bivariate generalized linear mixed model with {metafor}
>
>Dear Wolfgang,
>
>Thank you in advance. I contacted Reference [1] authors, and they sent me this
>data:
>
>structure(list(study.name = c("Cutler1995", "Dahlof1991", "Dowson2002",
>"Ensink1991", "Geraud2000", "Goadsby1991", "Goadsby2000", "Havanka2000",
>"Kaniecki2006", "Mathew2003", "Myllyla1998", "Nappi1994", "Patten1991",
>"Pfaffenrath1998", "Sandrini2002", "Sargent1995", "Sheftell2005", "Tfelt-
>Hansen1995", "Tfelt-Hansen1998", "Visser1996"), n1 = c(66L, 275L, 194L, 131L,
>504L, 89L, 129L, 98L, 131L, 831L, 42L, 142L, 142L, 277L, 170L, 46L, 902L, 122L,
>388L, 72L), n2 = c(65L, 182L, 99L, 78L, 56L, 93L, 142L, 91L, 127L, 419L, 41L,
>81L, 101L, 91L, 84L, 47L, 892L, 126L, 160L, 85L), r1 = c(37L, 180L, 123L, 60L,
>304L, 45L, 63L, 59L, 64L, 471L, 33L, 73L, 95L, 175L, 85L, 26L, 649L, 63L, 239L,
>33L), r2 = c(17L, 48L, 42L, 14L, 24L, 9L, 30L, 28L, 47L, 105L, 12L, 25L, 22L,
>27L, 25L, 8L, 375L, 30L, 64L, 15L)), row.names = c(NA, 20L), class =
>"data.frame")
>
>Best,
>
>Arthur M. Albuquerque
>
>Medical student
>Universidade Federal do Rio de Janeiro, Brazil
>
>On Mar 8, 2022, 7:07 PM -0300, Viechtbauer, Wolfgang (SP)
><wolfgang.viechtbauer using maastrichtuniversity.nl>, wrote:
>
>Hi Arthur,
>
>I don't want to start re-reading all those articles right now, but who cites who
>isn't the best indication of whether people are talking about the same model or
>not.
>
>In [1], the authors use data from a Cochrane review to illustrate the model. If
>you provide these data here using dput(), so that I can immediately reproduce the
>exact dataset, then I could see whether I can reproduce their results. It's time-
>consuming extracting data manually from pdfs, so I'll leave this up to you
>whether you want to do this.
>
>Best,
>Wolfgang
>
>-----Original Message-----
>From: Arthur Albuquerque [mailto:arthurcsirio using gmail.com]
>Sent: Tuesday, 08 March, 2022 2:07
>To: r-sig-meta-analysis using r-project.org; Viechtbauer, Wolfgang (SP)
>Subject: RE: [R-meta] Bivariate generalized linear mixed model with {metafor}
>
>Hi Wolfgang,
>
>It’s me again about this bivariate model. I am having a hard time trying to
>figure out if I understood it correctly.
>
>To recap, I wanted to fit a bivariate meta-analysis model (hereafter, mod1)
>described in Reference [1] below. You replied suggesting it was the "Model 6: the
>Van Houwelingen bivariate” (mod2) in your article with Jackson et al (Reference
>[2]).
>
>However, I am now re-reading all these articles and I believe mod1 and mod2 are
>different. Reference [1] cites Thompson et al. (Reference [3]), and does not cite
>van Houwelingen. You cited Van Houwelingen et al (Reference [4]). To my
>knowledge, they seem different models indeed.
>
>In fact, Van Houwelingen in Reference [5] directly cites Thompson suggesting
>these models are distinct:
>
>"The mix of many 􏱪fixed and a few random e􏱨ffects as proposed by Thompson et al.
>… are more in the spirit of the functional approach. These methods are meant to
>impose no conditions on the distribution of the true baseline risks… In a letter
>to the editor by Van Houwelingen and Senn following the article of Thompson et
>al. , Van Houwelingen and Senn argue that putting Bayesian priors on all nuisance
>parameters, as done by Thompson et al., does not help solving the inconsistency
>problem."
>
>Are they indeed different model?
>
>Please ignore this email if this question is out of the scope of your mailing
>list. Sorry in advance.
>
>Kind regards,
>
>Arthur M. Albuquerque
>
>Medical student
>Universidade Federal do Rio de Janeiro, Brazil
>
>References:
>
>[1] Xiao, Mengli, Yong Chen, Stephen R Cole, Richard F MacLehose, David B
>Richardson, and Haitao Chu. ‘Controversy and Debate: Questionable Utility of the
>Relative Risk in Clinical Research: Paper 2: Is the Odds Ratio “Portable” in
>Meta-Analysis? Time to Consider Bivariate Generalized Linear Mixed Model’.
>Journal of Clinical Epidemiology 142 (February 2022): 280–87.
>https://doi.org/10.1016/j.jclinepi.2021.08.004
>
>[2] Jackson, Dan, Martin Law, Theo Stijnen, Wolfgang Viechtbauer, and Ian R.
>White. ‘A Comparison of Seven Random-Effects Models for Meta-Analyses That
>Estimate the Summary Odds Ratio’. Statistics in Medicine 37, no. 7 (30 March
>2018): 1059–85. https://doi.org/10.1002/sim.7588
>
>[3] Thompson, Simon G., Teresa C. Smith, and Stephen J. Sharp. ‘Investigating
>Underlying Risk as a Source of Heterogeneity in Meta-Analysis’. Statistics in
>Medicine 16, no. 23 (15 December 1997): 2741–58.
>https://doi.org/10.1002/(SICI)1097-0258(19971215)16:23<2741::AID-SIM703>3.0.CO;2-
>0
>
>[4] Van Houwelingen, Hans C., Koos H. Zwinderman, and Theo Stijnen. ‘A Bivariate
>Approach to Meta-Analysis’. Statistics in Medicine 12, no. 24 (30 December 1993):
>2273–84. https://doi.org/10.1002/sim.4780122405
>
>[5] Houwelingen, Hans C. van, Lidia R. Arends, and Theo Stijnen. ‘Advanced
>Methods in Meta-Analysis: Multivariate Approach and Meta-Regression’. Statistics
>in Medicine 21, no. 4 (28 February 2002): 589–624.
>https://doi.org/10.1002/sim.1040