[R-sig-ME] Duplicating meta-regression results from PROC MIXED with lmer
Brant Inman
brant.inman at me.com
Thu May 7 05:00:38 CEST 2009
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
>
More information about the R-sig-mixed-models
mailing list