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?

  • 2
    $\begingroup$ A meta-regression is typically just a random effects model. Nothing stopping you from adding a quadratic, or going semiparametric. $\endgroup$ Oct 31, 2014 at 15:47
  • $\begingroup$ Hmm. No complications due to the additional between-study variance introduced by the meta-analytic nature of the data? $\endgroup$ Oct 31, 2014 at 18:19
  • $\begingroup$ 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). $\endgroup$ Oct 31, 2014 at 18:22

1 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:


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


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

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

figure of quadratic relationship

  • $\begingroup$ 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. $\endgroup$ Nov 1, 2014 at 19:02
  • $\begingroup$ Aha! My variances were entered as a factor for some reason. Thanks for your concise and clear answer. $\endgroup$ Nov 1, 2014 at 20:00
  • $\begingroup$ Glad to hear you found a solution. $\endgroup$
    – Wolfgang
    Nov 1, 2014 at 20:02
  • $\begingroup$ 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. $\endgroup$ Nov 9, 2014 at 20:29
  • $\begingroup$ 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()). $\endgroup$
    – Wolfgang
    Dec 17, 2015 at 9:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.