library(tidyverse) library(metafor) delta <- .5 tau <- .2 diamond_height <- .3 studies <- tibble( study = LETTERS, ni = round(runif(n = length(study), min = 100, max = 200)), true_i = delta + rnorm(n = length(study), mean = 0, sd = tau), yi = rnorm(n = length(study), mean = true_i, sd = sqrt(psychmeta::var_error_d(d = true_i, n1 = ni, n2 = ni))), se = sqrt(psychmeta::var_error_d(d = yi, n1 = ni, n2 = ni)), ci.lb = yi - qnorm(.975) * se, ci.ub = yi + qnorm(.975) * se ) plot_dat <- bind_rows( Studies = studies, Aggregate = rma(yi ~ 1, sei = se, data = studies, method = "REML", test = "knha") |> predict() |> as.data.frame() |> rename(yi = pred) |> mutate(study = "Random effects"), .id = "type" ) |> mutate(type = as_factor(type), study = as_factor(study), ypos = datawizard::reverse(as.integer(study))) ggplot(plot_dat) + aes(y = forcats::fct_rev(study), x = yi, xmin = ci.lb, xmax = ci.ub) + geom_linerange(data = \(x) filter(x, type == "Studies")) + geom_point(data = \(x) filter(x, type == "Studies")) + geom_linerange( aes(xmin = pi.lb, xmax = pi.ub), data = \(x) filter(x, type == "Aggregate"), linetype = "twodash" ) + geom_polygon( aes(x = x.diamond, y = y.diamond, group = study), data = \(x) filter(x, type == "Aggregate") |> group_by(study) |> reframe( x.diamond = c(ci.ub, yi, ci.lb, yi), y.diamond = ypos + c(0, diamond_height, 0, -diamond_height) ), inherit.aes = FALSE ) + labs(y = NULL, x = "Standardized mean difference") + scale_y_discrete(drop = FALSE) + # this one is key! see::theme_modern() + theme(strip.background = element_blank(), strip.text = element_blank())