library(metafor)
dat <- data.frame(
    study = paste0("S", 1:12),
    yi = c(
      0.42, 0.31, 0.35, 0.18,
      0.27, 0.14, 0.22, 0.05,
      0.15, -0.02, 0.08, -0.10
    ),
    
    sei = c(
      0.12, 0.10, 0.11, 0.09,
      0.10, 0.08, 0.09, 0.08,
      0.09, 0.07, 0.08, 0.07
    ),
    
    age = c(
      25, 28, 30, 32,
      35, 38, 40, 43,
      45, 48, 50, 55
    ),
    
    design = c(
      "Non-RCT", "RCT", "Non-RCT", "RCT",
      "Non-RCT", "RCT", "Non-RCT", "RCT",
      "Non-RCT", "RCT", "Non-RCT", "RCT"
    )
)
dat$vi <- dat$sei^2
dat$age_c <- dat$age - mean(dat$age)
dat$design <- factor(
  dat$design,
  levels = c("Non-RCT", "RCT")
)

model <- rma(
  yi = yi,
  vi = vi,
  mods = ~ age_c * design,
  data = dat,
  method = "REML"
)
summary(model)
model$QE
model$QM

model0 <- rma(
  yi = yi,
  vi = vi,
  data = dat,
  method = "REML"
)
model0

library(ggplot2)

newdat <- expand.grid(
  age_c = seq(
    min(dat$age_c),
    max(dat$age_c),
    length.out = 100
  ),
  design = levels(dat$design)
)

pred <- predict(
  model,
  newmods = model.matrix(
    ~ age_c * design,
    newdat
  )[,-1]
)

newdat$pred <- pred$pred
newdat$ci.lb <- pred$ci.lb
newdat$ci.ub <- pred$ci.ub


ggplot(dat,
       aes(x = age_c,
           y = yi,
           size = 1/sei,
           shape = design)) +
  
  geom_point(alpha = 0.8) +
  
  geom_ribbon(
    data = newdat,
    aes(
      x = age_c,
      ymin = ci.lb,
      ymax = ci.ub,
      fill = design
    ),
    alpha = 0.15,
    inherit.aes = FALSE
  ) +
  
  geom_line(
    data = newdat,
    aes(
      x = age_c,
      y = pred,
      linetype = design
    ),
    linewidth = 1,
    inherit.aes = FALSE
  ) +
  
  labs(
    x = "Centered mean age",
    y = "Effect size",
    size = "Precision",
    shape = "Study design",
    linetype = "Study design",
    fill = "Study design",
    title = "Advanced Meta-Regression"
  ) +
  
  theme_minimal()