### loading packages
library(metafor)
library(ggplot2)
library(scales)
library(ggthemes) # for more layout presets
library(xkcd) # for even more layout presets

### run the meta-analysis and the prediction (as per www.metafor-project.org)
dat <- escalc(measure = "RR", ai = tpos, bi = tneg, ci = cpos, di = cneg, data = dat.bcg)

### fit mixed-effects model with absolute latitude as predictor
res <- rma(yi, vi, mods = ~ ablat, data = dat)

### calculate predicted risk ratios for 0 to 60 degrees absolute latitude
preds <- predict(res, newmods = c(0:60), transf = exp)
plot.preds <- as.data.frame(preds) # ggplot() can't read list.rma() objects
plot.preds$ablat <- c(0:60) # adding the newmods for plotting in a new row

### calculate point sizes by rescaling the standard errors
wi    <- 1/sqrt(dat$vi)
size  <- 0.5 + 3.0 * (wi - min(wi))/(max(wi) - min(wi))

### extract plotting data from escalc() model
plot.dat <- summary.escalc(dat)

### set up the plot
ggplot(plot.dat, aes(ablat, exp(yi))) +
  geom_point(pch = 19, cex = 1.5 * size) + # as, to my experience, ggplot() will render dots smaller than the plot() function
  theme_classic() + # insert your preferred style here
  geom_linerange(aes(ymin = exp(plot.dat$ci.lb), 
                     ymax = exp(plot.dat$ci.ub)), 
                 size = .25) + # gives straight lines w/o brackets at the ends. if „classical“ look is preferred, use geom_errorbar()
  scale_y_continuous(trans = log_trans(), 
                     limits = c(.2, 1.6),
                     breaks = c(.2, .4, .6, .8, 1, 1.2, 1.4, 1.6)) + # for adjustment of axes, see (as example) http://www.cookbook-r.com/Graphs/
  xlab("Absolute Latitude") + 
  ylab("risk ratio") +
  geom_hline(yintercept = 1, linetype = "dotted") + # draws the dotted line at y = 1 for equal risk ratio
  geom_line(data = plot.preds, aes(x = ablat, y = pred)) + # plots the regression line
  geom_line(data = plot.preds, aes(x = ablat, y = ci.lb), linetype = "dashed") + # plots the lower confidence interval boundary
  geom_line(data = plot.preds, aes(x = ablat, y = ci.ub), linetype = "dashed") # plots the upper confidence interval boundary