Duplicating meta-regression results from PROC MIXED with lmer

Wed, May 6, 2009 8:00 PM

Ken,

Thanks very much for the great advice.  You essentially suggest the  
following:

1) Add interaction term b/t covariates and treatment arm variable
2) Center continuous covariates
3) Use a wide dataset (improves computation time dramatically over a  
long dataset)

The estimates obtained with your method would seem at first glance to  
compare favorably with the published results.  However, on closer  
analysis I wonder if the interpretation of your model is the same as  
Van Houwelingen's.




YOUR MODEL, WHERE YEAR AND LATITUDE ARE ACTUALLY CENTERED VARIABLES

Generalized linear mixed model fit by the Laplace approximation
Formula: y ~ bcg * latitude + bcg * year + (bcg | study)
    Data: newbcg
    AIC   BIC logLik deviance
  105.6 116.9  -43.8     87.6
Random effects:
  Groups Name        Variance   Std.Dev. Corr
  study  (Intercept) 0.60620909 0.778594
         bcgyes      0.00045554 0.021343 1.000
Number of obs: 26, groups: study, 13

Fixed effects:
                  Estimate Std. Error z value Pr(>|z|)
(Intercept)     -4.116305   0.221649 -18.571  < 2e-16 ***
bcgyes          -0.733330   0.049815 -14.721  < 2e-16 ***
latitude         0.024016   0.021007   1.143  0.25294
year            -0.099948   0.027955  -3.575  0.00035 ***
bcgyes:latitude -0.034332   0.003991  -8.601  < 2e-16 ***
bcgyes:year     -0.001489   0.005811  -0.256  0.79781
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.?  
0.1 ? ? 1

Correlation of Fixed Effects:
             (Intr) bcgyes latitd year   bcgys:l
bcgyes       0.046
latitude    -0.004  0.006
year        -0.015  0.014  0.658
bcgyes:lttd  0.002  0.156  0.088  0.059
bcgyes:year  0.020 -0.234  0.049  0.078  0.706


VAN HOUWELINGEN's MODEL RESULTS (I BELIEVE THAT HE DID NOT CENTER THE  
VARIABLES):

5.2.4. Regression on latitude and year. When both covariates latitude  
and year are put into
the model, the residual between-study variance becomes only 0.002,  
corresponding with an
explained variance of 99.3 per cent, only slightly more than by  
latitude alone. The regression
coe??cients for the intercept, latitude and year are, respectively,  
0.494 (standard error = 0:529),
?0:034 (standard error = 0:004) and ?0:001 (standard error = 0:006).



After taking into account centering, do you believe that your  
regression equation is interpreted in the same way as his?

Brant

On May 6, 2009, at 8:02 AM, Ken Beath wrote:

On 06/05/2009, at 8:43 PM, Brant Inman wrote:

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

Duplicating meta-regression results from PROC MIXED with lmer

Thread (8 messages)