[R-sig-ME] Duplicating meta-regression results from PROC MIXED with lmer
Ken Beath
ken at kjbeath.com.au
Wed May 6 14:02:57 CEST 2009
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