Duplicating meta-regression results from PROC MIXED with lmer

Rolf,

I actually don't believe that this is a SAS vs R issue since I have  
3 sources that report the same results.  I know that STATA, SAS and  
the mima function from R can all be used to give the correct  
results.  The question is related more to how I can get similar  
results with lmer.  Currently, the code I provided generates VERY  
different results, especially related to between study variance  
estimation (tau) and regression coefficient standard errors.  
Basically, I think I got it all wrong in my coding of the models  
with covariates (f3 and f4), but I can't understand why.

Brant

This is closer but maybe not quite equivalent. It also avoids using  
the huge data file.

Ken

bcg <- read.csv("bcg.csv")
bcg$study <- paste(bcg$author, bcg$year)

# this is a bit messy but requires almost no thought
bcgnames <- bcg[,c(1,2,7:10)]
bcgvacc <- bcg[,3:4]
names(bcgvacc) <- c("tb","n")
bcgnovacc <- bcg[,5:6]
names(bcgnovacc) <- c("tb","n")

newbcg <- cbind(rbind(bcgnames,bcgnames),rbind(bcgvacc,bcgnovacc))
newbcg$bcg <- factor(rep(c("yes","no"),each=13))
newbcg$notb <- newbcg$n-newbcg$tb

newbcg$latitude <- newbcg$latitude-mean(newbcg$latitude)
newbcg$year <- newbcg$year-mean(newbcg$year)

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(cbind(tb,notb) ~ bcg + (1 | study), family=binomial,  
data=newbcg)
	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(cbind(tb,notb) ~ bcg + (bcg | study), family=binomial,  
data=newbcg) # 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

# Random effects model, 1 covariate
f3 <- lmer(cbind(tb,notb) ~ latitude*bcg + (bcg | study),  
family=binomial, data=newbcg)
	summary(f3)
	exp(fixef(f3)) # OR

# Random effects model, 1 covariate
f4 <- lmer(cbind(tb,notb) ~ latitude*bcg + year*bcg + (bcg | study),  
family=binomial, data=newbcg)
	summary(f4)
	exp(fixef(f4)) # OR