nonlinear meta-regression

Question

Can someone point me to a basic explanation of theory and methods for fitting nonlinear curves (particularly quadratic functions) to meta-analytic data? I have a set of effect sizes that are clearly hump-shaped with a peak at the intermediate value of the explanatory variable in a meta-regression. I can only find examples of linear models for meta-regression, which are clearly not suited to the data.

Or, am I missing something and there's a reason that I can't find nonlinear meta-regressions?

Answer 1

Are you trying to fit a quadratic polynomial or a truly non-linear model? If you are trying to fit a quadratic polynomial, then this is easy with pretty much any meta-analysis software that allows you to specify a (mixed-effects) meta-regression model with multiple predictors. You just include the explanatory variable and its squared version as predictors in the model. Here is an example using the ''metafor'' package in R:

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

The results are:

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

Not surprisingly, the quadratic term is highly significant (p < .0001).

Let's see what this looks like in a plot:

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))