CODEWITHSUNDEEP
rmetafor
user115916
user115916

Reputation: 188

metafor provides 95%CI that are different from the original values

I am using metafor package for combining beta coefficients from the linear regression model. I used the following code. I supplied the reported se and beta values for the rma function. But, when I see the forest plot, the 95% confidence intervals are different from the ones reported in the studies. I also tried it using mtcars data set by running three models and combining the coefficients. Still, the 95%CI we see on the forest plot are different from the original models. The deviations are far from rounding errors. A reproducible example is below.

library(metafor)
library(dplyr)

lm1 <- lm(hp~mpg, data=mtcars[1:15,])
lm2 <- lm(hp~mpg, data=mtcars[1:32,])
lm3 <- lm(hp~mpg, data=mtcars[13:32,])

study <- c("study1", "study2", "study3")
beta_coef <- c(lm1$coefficients[2], 
               lm2$coefficients[2], 
               lm3$coefficients[2]) %>% as.numeric()
se <- c(1.856, 1.31,1.458) 

ci_lower <- c(confint(lm1)[2,1], 
              confint(lm2)[2,1], 
              confint(lm3)[2,1]) %>% as.numeric()

ci_upper <- c(confint(lm1)[2,2], 
              confint(lm2)[2,2], 
              confint(lm3)[2,2]) %>% as.numeric()

df <- cbind(study=study, 
            beta_coef=beta_coef, 
            se=se, 
            ci_lower=ci_lower, 
            ci_upper=ci_upper) %>% as.data.frame()

pooled <- rma(yi=beta_coef, vi=se, slab=study)

forest(pooled)

enter image description here

Compare the confidence intervals on the forest plot with the one on the data frame.

enter image description here

data frame

df <- cbind(study=study, 
            beta_coef=beta_coef, 
            se=se, 
            ci_lower=ci_lower, 
            ci_upper=ci_upper) %>% as.data.frame()

Upvotes: 2

Views: 443

Answers (1)

Wolfgang
Wolfgang

Reputation: 3395

Argument vi is for specifying the sampling variances, but you are passing the standard errors to the argument. So you should do:

pooled <- rma(yi=beta_coef, sei=se, slab=study)

But you will still find a discrepancy here, since the CIs in the forest plot are constructed based on a normal distribution, while the CIs you obtained from the regression model are based on t-distributions. If you want the exact same CIs in the forest plot, you could just pass the CI bounds to the function like this:

forest(beta_coef, ci.lb=ci_lower, ci.ub=ci_upper)

If you want to add a summary polygon from some meta-analysis to the forest plot, you can do this with addpoly(). So the complete code for this example would be:

forest(beta_coef, ci.lb=ci_lower, ci.ub=ci_upper, ylim=c(-1.5,6))
addpoly(pooled, row=-1)
abline(h=0)

Upvotes: 2

Related Questions