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() a function/wrapper

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

Can further modify, add layers

# further example here...