set.seed(12326)

### number of studies
n <- 40

### simulate explanatory variable
xi <- round(runif(n, 1, 10), 1)

### simulate sampling variances
vi <- rgamma(n, 2, 2)/20

### simulate estimates (quadratic relationship + residual heterogeneity + sampling error)
yi <- 0.5 + 0.3 * xi - 0.03 * xi^2 + rnorm(n, 0, 0.1) + rnorm(n, 0, sqrt(vi))

library(metafor)

res <- rma(yi, vi, mods = ~ xi + I(xi^2))
res

Mixed-Effects Model (k = 40; tau^2 estimator: REML)

tau^2 (estimated amount of residual heterogeneity):     0.0234 (SE = 0.0135)
tau (square root of estimated tau^2 value):             0.1530
I^2 (residual heterogeneity / unaccounted variability): 43.17%
H^2 (unaccounted variability / sampling variability):   1.76
R^2 (amount of heterogeneity accounted for):            69.27%

Test for Residual Heterogeneity:
QE(df = 37) = 61.0438, p-val = 0.0077

Test of Moderators (coefficient(s) 2,3):
QM(df = 2) = 37.5033, p-val < .0001

Model Results:

         estimate      se     zval    pval    ci.lb    ci.ub
intrcpt    0.2655  0.2024   1.3119  0.1896  -0.1311   0.6621
xi         0.4073  0.0834   4.8824  <.0001   0.2438   0.5708  ***
I(xi^2)   -0.0402  0.0073  -5.4753  <.0001  -0.0546  -0.0258  ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

plot(xi, yi, pch=19, cex=.2/sqrt(vi), xlim=c(1,10), ylim=c(0,2))
xi.new <- seq(1,10,length=100)
pred <- predict(res, newmods = cbind(xi.new, xi.new^2))
lines(xi.new, pred$pred)
lines(xi.new, pred$ci.lb, lty="dashed", col="gray")
lines(xi.new, pred$ci.ub, lty="dashed", col="gray")
points(xi, yi, pch=19, cex=.2/sqrt(vi))