Skip to content
Prev 3927 / 5636 Next

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

Wolfgang Viechtbauer 
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

Thread (9 messages)

Arthur Albuquerque Bivariate generalized linear mixed model with {metafor} Mar 7 Wolfgang Viechtbauer Bivariate generalized linear mixed model with {metafor} Mar 8 Arthur Albuquerque Bivariate generalized linear mixed model with {metafor} Mar 9 Wolfgang Viechtbauer Bivariate generalized linear mixed model with {metafor} Mar 9 Arthur Albuquerque Bivariate generalized linear mixed model with {metafor} Mar 9 Wolfgang Viechtbauer Bivariate generalized linear mixed model with {metafor} Mar 9 Arthur Albuquerque Bivariate generalized linear mixed model with {metafor} Mar 9 Wolfgang Viechtbauer Bivariate generalized linear mixed model with {metafor} Mar 9 Arthur Albuquerque Bivariate generalized linear mixed model with {metafor} Mar 9