Duplicating meta-regression results from PROC MIXED with lmer

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