Load metafor package
library("metafor")
## Loading required package: Matrix
## Loading 'metafor' package (version 1.9-7). For an overview
## and introduction to the package please type: help(metafor).
### load BCG vaccine data
data(dat.bcg)
### calculate log relative risks and corresponding sampling variances
dat <- escalc(measure="RR", ai=tpos, bi=tneg, ci=cpos, di=cneg, data=dat.bcg)
### random-effects model (method="REML" is default, so technically not needed)
resRMA = rma(yi, vi, data=dat)
resRMA
##
## Random-Effects Model (k = 13; tau^2 estimator: REML)
##
## tau^2 (estimated amount of total heterogeneity): 0.3132 (SE = 0.1664)
## tau (square root of estimated tau^2 value): 0.5597
## I^2 (total heterogeneity / total variability): 92.22%
## H^2 (total variability / sampling variability): 12.86
##
## Test for Heterogeneity:
## Q(df = 12) = 152.2330, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## -0.7145 0.1798 -3.9744 <.0001 -1.0669 -0.3622 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
forest(resRMA)
ggforest = function(x){
require("ggplot2")
# Function to convert REM results in `rma`-format into a data.frame
rma2df = function(x){
rbind(
data.frame(Study = "RE Model", LogFC = x$b, CILB=x$ci.lb, CIUB=x$ci.ub,
p = x$pval,
stringsAsFactors = FALSE),
data.frame(Study = x$slab, LogFC = x$yi,
CILB=x$yi - 2*sqrt(x$vi),
CIUB=x$yi + 2*sqrt(x$vi),
p = x$pval,
stringsAsFactors = FALSE)
)
}
remresdf = rma2df(x)
remresdf <- transform(remresdf, interval = CIUB - CILB)
remresdf <- transform(remresdf, RelConf = 1/interval)
p = ggplot(remresdf,
aes(LogFC, Study, xmax=CIUB, xmin=CILB)) +
coord_cartesian(xlim=c(-2, 2)) +
#scale_alpha_discrete(range = c(0.2, 1)) +
geom_vline(xintercept = 0.0, linetype=2, alpha=0.75) +
geom_errorbarh(alpha=0.5, color="black") +
geom_point(aes(size = RelConf)) +
geom_point(data = subset(remresdf, Study=="RE Model"), size=7) +
scale_size(range = c(2, 5), guide=FALSE) +
theme_bw() +
theme(text = element_text(size=16))
return(p)
}
pforest = ggforest(resRMA)
## Loading required package: ggplot2
pforest
# further example here...