[R-sig-ME] Duplicating meta-regression results from PROC MIXED with lmer

[Brant Inman]

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
>