Plot slope and intercept of a model with log and percentage transformations

Question

This is the structure of my data

 > dput(test)
structure(list(MAT = c(4.9, 4.9, 15.5, 14.1, 14.1, 14.1, 11.5, 
11.5, 11.5, 17, 6.1, 2.7, 2.2, 2.2, 14.1, 14.1, 14.1, 9.5, 9.5, 
9.5, 9.5, 9.3, 8.3, 8.266666651, 8.266666651, 4.3, 4.3, 22.3, 
14.1, 14.1, 14.1, 8.5, 8.5, 8.5, 8.5, 21.5, 21.5, 3.8, 3.8, 6, 
6, 6, 6, 6), es = c(0.29603085763985, 0.421393627439682, 0.189653473156549, 
0.226685054608428, 0.291373762079697, 0.166533544378467, 0.250586529054368, 
0.146320008054403, 0.199565119644333, -0.0819047677231083, 0.15963948187092, 
-0.154628141843561, 0.201121044198443, 0.0867981239977565, 0.543870310978598, 
0.34547921143505, 0.37557241352574, -0.287318919407836, 0.207937483228907, 0.190143660810163, 0.276182673435993, 0.128596803172119, 0.454753165843559, 
0.399237234440439, 0.32075358541748, 0.362664873575803, -0.0865925288159671, 
0.51290512543514, 0.186308318839249, 0.147936083867325, 0.243792477087184, 
0.625169403695832, 0.110317782120045, 0.217836235313289, 0.171468156841181, 
0.50548821117127, 0.164418265301427, -0.00246305543239786, 0.325552346507191, 
0.381240606108843, 0.19337350462531, 0.0408803528990759, 0.321815078821239, 
0.307642815014319), var = c(0.00496277337027962, 0.0130962311273343, 
0.0180149624217804, 0.0134568083459063, 0.00139708925143695, 
0.000725862546533828, 0.00670831011660164, 0.0190783110089115, 
0.0641568910090007, 0.0121596544795352, 0.0653909966557582, 0.0514610437228611, 
0.0231592619167496, 0.0108989891148006, 0.0588577146414195, 0.0695760532112402, 
0.0744256820906048, 0.00997789089155498, 0.00928124381998638, 
0.0145009450673482, 0.00652956018299188, 0.0111886178917916, 
0.0265943757419349, 0.142676904340634, 0.110705177803624, 0.0576538348777718, 
0.0625171635976251, 0.0131652117394448, 0.00947904166717649, 
0.00813569411386797, 0.00444289889858652, 0.0673007030900184,0.00545169559098343, 0.240046081413733, 0.00561125010476281, 
    0.0185516235174018, 0.0179989506841957, 0.0496806959944248, 0.022478393723115, 
    0.0521209786580004, 0.282298667080106, 0.0151428845076692, 0.00992945920656693, 0.0145544965304081), MAP = c(810, 810, 1140, 1750, 1750, 1750, 
    1034, 1034, 1034, 720, 645, 645, 645, 645, 1000, 1000, 1000, 
    691, 691, 691, 691, 1134, 1750, 1326, 1326, 1140, 1140, 1310, 
    1750, 1750, 1750, 1003, 1003, 1003, 1003, 1750, 1750, 1750, 1750, 
    1750, 1750, 1750, 1750, 1750), CO2dif = c(162L, 162L, 190L, 165L, 
    165L, 165L, 200L, 200L, 200L, 150L, 335L, 335L, 335L, 335L, 348L, 
    348L, 348L, 200L, 200L, 200L, 200L, 220L, 350L, 350L, 350L, 350L, 
    350L, 350L, 180L, 180L, 180L, 130L, 130L, 130L, 130L, 320L, 320L, 
    360L, 360L, 345L, 345L, 350L, 348L, 348L)), row.names = c(NA, 
    -44L), class = "data.frame", .Names = c("MAT", "es", "var", "MAP", 
    "CO2dif"))

I run model selection using meta-analysis, and the best model to predict the effects size is:

library(metafor)

summary(rma(es, var, data=test ,control=list(stepadj=.5), mods= ~ 1 + log(MAT) + MAP + CO2dif + log(MAT):CO2dif, knha=TRUE))

Model Results:

                 estimate      se     tval    pval    ci.lb    ci.ub     
intrcpt            1.2556  0.3719   3.3758  0.0017   0.5033   2.0080   **
log(MAT)          -0.5740  0.1694  -3.3882  0.0016  -0.9167  -0.2313   **
MAP                0.0001  0.0001   2.5181  0.0160   0.0000   0.0003    *
CO2dif            -0.0042  0.0013  -3.2932  0.0021  -0.0067  -0.0016   **
log(MAT):CO2dif    0.0020  0.0005   3.7500  0.0006   0.0009   0.0031  ***

Now I want to plot es vs MAT, with an example with this model, assuming that MAP=1200 mm and CO2dif=350

MAPi <- 1200
CO2i <- 350
make_pct <- function(x) (exp(x) - 1) * 100

ggplot(test, aes(x = log(MAT), y = make_pct(es))) +
  geom_abline(aes(intercept = make_pct(1.2556 + 0.0001 * MAPi - 0.0042 * CO2i),
                  slope = make_pct(log(0.0020 * CO2i)) - make_pct(log(0.5740))) , 
                  color = "red", size=0.8) +
  geom_point() +
  theme_classic()

Effect size (es) is in log format, and I want percentage, so I transform it with the function make_pct. MAT, on the other hand, has to be log-transformed in the plot as indicated in the model output. Is the slope of the ggplot above correct with the log and percentage transformations? It seems to me that the slope is rather low. I am not very familiar with this type of plots and transformations, so any tips are welcome. Thanks

Answer 1

The relationship between exp(es)-1 and the explanatory variable log(MAT) is not linear.
For a given set of values of MAP and CO2dif, this relationship is of the form:
y = exp(es)-1 = k1*exp(k2*log(MAT)).
This function can be plotted as follows:

library(metafor)
library(ggplot2)

modfit <- rma(es, var, data=test ,control=list(stepadj=.5), 
   mods= ~ 1 + MAP + log(MAT)*CO2dif, knha=TRUE)
pars <- coef(modfit)

MAPi <- 1200
CO2i <- 350

make_pct <- function(x) (exp(x) - 1) * 100

mod_fun <- function(MAP, MAT, CO2dif, pars) {
    y <- pars[1]+pars[2]*MAP+pars[3]*log(MAT)+
         pars[4]*CO2dif+pars[5]*log(MAT)*CO2dif
    make_pct(y)
}

test$ESpct <- mod_fun(MAPi, test$MAT, CO2i, coef(modfit))

ggplot(test, aes(x = log(MAT), y = make_pct(es))) +
  geom_line(aes(y=ESpct), color = "red", size=0.8) +
  geom_point() + theme_classic()