# nonlinear meta-regression

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?

• A meta-regression is typically just a random effects model. Nothing stopping you from adding a quadratic, or going semiparametric. Commented Oct 31, 2014 at 15:47
• Hmm. No complications due to the additional between-study variance introduced by the meta-analytic nature of the data? Commented Oct 31, 2014 at 18:19
• And, does anyone know of a meta-analysis specific software that can include quadratic terms in random effects models? I'm working in Meta-Win right now (I know, I know... dark ages). Commented Oct 31, 2014 at 18:22

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


• I get this error when defining the model. ind.mod = rma(d, Var.d., mods = ~ Prod.g.m2 + I(Prod.g.m2^2)) Error in rma(d, Var.d., mods = ~Prod.g.m2 + I(Prod.g.m2^2)) : Processing terminated since k = 0 In addition: Warning messages: 1: In Ops.factor(Prod.g.m2, 2) : ^ not meaningful for factors 2: In rma(d, Var.d., mods = ~Prod.g.m2 + I(Prod.g.m2^2)) : Studies with NAs omitted from model fitting. Without "I" it runs, but gives a crazy # of predictors. Thrilled you wrote metafor! Much thanks for any help. Commented Nov 1, 2014 at 19:02
• Aha! My variances were entered as a factor for some reason. Thanks for your concise and clear answer. Commented Nov 1, 2014 at 20:00
• Glad to hear you found a solution. Commented Nov 1, 2014 at 20:02
• Could you clarify what the function predict is doing in the plot code you provided? I read the predict.rma help and see that it returns "the estimated (average) outcomes for values of the moderator." I suspect these are not 95% conf intervals, though, since so few of the yi values are within the dashed lines. Commented Nov 9, 2014 at 20:29
• The dashed lines are pointwise 95% CIs for the estimated average outcome as a function of the moderator(s). This doesn't imply that the individual observations (yi values) are expected to fall within those bounds most of the time. You seem to be thinking of prediction/credibility intervals here, which are a different thing (those are also provided by predict()). Commented Dec 17, 2015 at 9:00