rpredictionpredictmetaformeta-analysis

predict function in metafor for polynomials


I would like to create prediction intervals for a meta-analysis rma.mv object.

I can run this fine when my rma.mv object is linear.

However, when my rma.mv object is a cubic polynomial I get the error

Could not match up variable 'c_treattemp' uniquely to a variable in the model.

Could anyone help me?

Here is my code so far

### Linear model of effect of c_treattemp
meta_lm <- rma.mv(es, VCV,  mod= ~c_treattemp), random= list(~ 1|study_code,  ~1|obs), data= rdata, method= "REML")

### Cubic model of effect of c_treattemp
meta_3 <- rma.mv(es, VCV,  mod= ~poly(c_treattemp, degree=3, raw=TRUE), random= list(~ 1|study_code,  ~1|obs), data= rdata, method= "REML")

### Make sample of 100 datapoints between min/max of observed c_treattemp
sampledata <- as.matrix(data.frame(c_treattemp = seq(min(rdata$c_treattemp), max(rdata$c_treattemp), 
                                                     length.out = 100)))  

### Make predictions using the linear model.
preds <- data.frame(predict.rma(meta_lm, newmods = sampledata, digits = 2, addx=TRUE))  #<<< This works. 

### Make predictions using the cubic model.
preds <- data.frame(predict.rma(meta_3, newmods = sampledata, digits = 2, addx=TRUE))  #<<<  ERROR.  

I guess I need to create some sort of quadratic and cubic input in the sample data, but I am not sure how.


Solution

  • This tutorial on the metafor website shows how this can be done:

    https://www.metafor-project.org/doku.php/tips:non_linear_meta_regression

    In essence, you need to use:

    xs <- seq(min(rdata$c_treattemp), max(rdata$c_treattemp), length.out=100)
    preds <- data.frame(predict(meta_3, newmods=unname(poly(xs, degree=3, raw=TRUE)), addx=TRUE))
    

    Sidenote: The digits=2 part is superfluous, since this only affects printing, and if you turn the object into a data frame, then the values will always be unrounded anyway. Also, do not call predict.rma() directly - instead, call the generic function predict() and let R handle the dispatching to the appropriate method function (which is predict.rma()).