[R-sig-ME] Duplicating meta-regression results from PROC MIXED with lmer
Brant Inman
brant.inman at me.com
Wed May 6 05:00:04 CEST 2009
R-experts:
In 2002, Hans Van Houwelingen et al. published a tutorial on how to do
meta-regression in Statistics in Medicine. They used the classic BCG
dataset of Colditz to demonstrate correct methodology and computed the
results using PROC MIXED in SAS. In trying to duplicate the results
presented in this paper, I have discovered that I can reproduce
certain items with lmer but not others. I was hoping that someone
might point out how I could correctly program R code to arrive at the
correct solution. I have placed the paper and the datasets used below
at the following website for easy access to any helpers: http://www.duke.edu/~bi6
I start by loading the data into R.
-----
bcg <- read.csv('bcg.csv')
bcg.long <- read.csv('bcg-long.csv')
bcg$study <- paste(bcg$author, bcg$year)
-----
I then perform standard meta-analysis using two different R functions:
(1) the metabin function (meta package) to pool odds ratios with
standard inverse variance techniques using the "wide" bcg dataset and
(2) the lmer function (lme4 package) to perform a multilevel meta-
analysis using the "long" dataset. The only reason that the lmer
function is possible here is because the outcome is binary
(disease .vs. no disease) and I could create a long dataset. A 3rd
option is the mima function, but that is not presented here since I am
interested in using lmer to extend to situations where there are study
level (level 2) and individual level (level 1) predictors, something
mima cannot currently handle.
-----
library(meta)
# Fixed and random effects models, no covariates
f0 <- metabin(bcg[,3], bcg[,4], bcg[,5], bcg[,6], sm='OR',
method='Inverse')
summary(f0)
library(lme4)
# Fixed effects model, no covariates
f1 <- lmer(tb ~ bcg + (1 | study), family=binomial, data=bcg.long)
summary(f1)
exp(fixef(f1)[2]) # OR
exp(f1 at fixef[2] - (1.96*sqrt(vcov(f1)[2,2]))) # lci
exp(f1 at fixef[2] + (1.96*sqrt(vcov(f1)[2,2]))) # uci
# Random effects model, no covariates.
f2 <- lmer(tb ~ bcg + (bcg | study), family=binomial,
data=bcg.long) # Random effects, no covariates
summary(f2)
exp(fixef(f2)[2]) # OR
exp(f2 at fixef[2] - (1.96*sqrt(vcov(f2)[2,2]))) # lci
exp(f2 at fixef[2] + (1.96*sqrt(vcov(f2)[2,2]))) # uci
as.numeric(summary(f2)@REmat[2,3]) # Tau
-----
So far these results look good and compare favorably to those obtained
by Van Houwelingen. It is also obvious that although metabin and lmer
give similar results, computational time is much longer with lmer than
meta since it must use the long version of the dataset. The problem
comes when two covariates are added to model f2, latitude and year of
publication.
-----
# Random effects model, 1 covariate
f3 <- lmer(tb ~ bcg + latitude + (bcg | study), family=binomial,
data=bcg.long)
summary(f3)
exp(fixef(f3)) # OR
# Random effects model, 1 covariate
f4 <- lmer(tb ~ bcg + year + (bcg | study), family=binomial,
data=bcg.long)
summary(f4)
exp(fixef(f4)) # OR
-----
I assumed, incorrectly, that models f3 and f4 would reproduce the
results of Van Houwelingen in sections 5.2.1 and 5.2.2. They do not
seem to do so. I would be very appreciative if someone pointed out my
error with models f3 and f4 and why they do not seem to be correct.
Incidentally, other sources (ex: Egger/Altman book on systematic
reviews) report results on this dataset similar to Van Houwelingen, so
I think that my code is definitely the problem.
Thanks,
Brant Inman
More information about the R-sig-mixed-models
mailing list