forest plot

R code of the ggplot2 example

library(ggplot2)
library(dplyr)
library(HistData)

## Rate ratio of vaccinated children to their controls, with its 95%
## interval on the log scale: 1 / a + 1 / b is the variance of log(ratio).
ratio <- function(a, n_a, b, n_b) {
  r <- (a / n_a) / (b / n_b)
  half <- qnorm(0.975) * sqrt(1 / a + 1 / b)
  c(ratio = r, low = r * exp(-half), high = r * exp(half))
}
trials <- PolioTrials |>
  filter(Group %in% c("Vaccinated", "Placebo", "Controls"))
rows <- bind_rows(lapply(c("Paralytic", "NonParalytic"), function(case) {
  bind_rows(lapply(c("RandomizedControl", "ObservedControl"), function(design) {
    t <- filter(trials, Experiment == design)
    v <- t[t$Group == "Vaccinated", ]
    k <- t[t$Group != "Vaccinated", ]
    tibble(case = case, design = design,
           !!!as.list(ratio(v[[case]], v$Population, k[[case]], k$Population)))
  }))
})) |>
  mutate(label = paste0(ifelse(case == "Paralytic", "Paralytic", "Non-paralytic"), ", ",
                        ifelse(design == "RandomizedControl", "placebo trial", "observed controls")),
         label = factor(label, levels = rev(label)))

p <- ggplot(rows, aes(x = ratio, y = label)) +
  geom_vline(xintercept = 1, colour = "#535353", linewidth = 0.4, linetype = "dashed") +
  geom_errorbar(aes(xmin = low, xmax = high), width = 0.2, orientation = "y", colour = "#9b2367", linewidth = 0.5) +
  geom_point(shape = 15, size = 2.4, colour = "#9b2367") +
  scale_x_log10(breaks = c(0.25, 0.5, 1, 2), labels = c("0.25", "0.5", "1", "2")) +
  labs(title = "The vaccine cut paralytic polio by 62% to 71%",
       subtitle = "Each row is a separate comparison (trial design × form of polio), 1954; rate ratio, vaccinated / controls, with 95% intervals",
       x = "Rate ratio (log scale)", y = NULL,
       caption = "Data: Francis et al. 1955, via the HistData R package") +
  theme_lex_viz(grid = "x", arrows = "")