I am currently doing a meta-analysis looking at the effects of logging in tropical forests.
As part of this I have been testing hypotheses about whether the effects vary by region and method of logging used.
I am doing all of this using the metafor package in R.
My data looks like this:
structure(list(Method = structure(c(2L, 2L, 1L, 1L, 1L, 1L, 1L,
2L, 2L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L
), .Label = c("Conventional", "RIL"), class = "factor"), MU = c(192.96,
252.41, 235.6, 258, 258, 399, 313, 409.8, 420.4, 333.47, 327.54,
256, 228.1, 547.1, 453.3873094, 427.495, 346.8, 330.833333333333,
343.3, 221.5, 194.8, 51.1, 276), SSU = c(3, 3, 30, 3, 3, 2, 5,
17, 10, 4, 4, 4, 9, 15, 35, 10, 3, 3, 3, 3, 3, 3, 10), ML = c(157.03,
171.97, 219.5, 198, 148, 191, 204, 315.3647059, 386.22, 135.8,
211.78, 183.8, 159.9, 230.8, 97.00798294, 218.31, 279.933333333333,
261.4, 249.733333333333, 118.6, 42.9, 18.7, 128.4), SSL = c(3,
3, 10, 3, 3, 10, 5, 17, 10, 4, 4, 4, 9, 10, 131, 45, 3, 3, 3,
3, 3, 3, 10), Region = structure(c(3L, 3L, 2L, 2L, 2L, 3L, 2L,
2L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 3L, 3L, 3L, 3L
), .Label = c("Africa", "Americas", "Asia & Oceania"), class = "factor"),
SDU = c(7.69030558560582, 12.1243556529821, 74.4902678207026,
30, 30, 145, 107.33126291999, 64.9, 92.95, 40.73364703, 54.0371067,
53.6, 98.1, 193.8, 16.13693527, 109.3250955, 28.21329474,
30.91865671942, 32.220024829289, 37.065887281974, 96.4752299815865,
37.4122974434878, 91.706052144883), SDL = c(8.46972844901181,
7.81154914213564, 53.1262646908288, 18, 10, 324.8738217,
84.970583144992, 44.90907399, 109.0794186, 20.75198304, 18.6400617,
11.6, 88.2, 104.2, 4.008416039, 185.9464001, 29.85034897,
28.7292533839639, 15.297494348204, 37.7587076050015, 32.5625551822949,
7.44781847254617, 126.174878640718)), .Names = c("Method",
"MU", "SSU", "ML", "SSL", "Region", "SDU", "SDL"), row.names = c(NA,
23L), class = "data.frame")
I then used this to calculate effect sizes and associated SEs for each site I have data for, like this:
require("metafor")
ROM <- escalc(data=AGB, measure="ROM", m2i=MU, sd2i=SDU,
n2i=SSU, m1i=ML, sd1i=SDL, n1i=SSL, append=TRUE)
My problem is that I don't know how to interpret the treatment contrasts from a model with two categorical predictors.
My 'best' model (the one with lowest AIC) looks like this:
ROM.ma1 <- rma(yi,vi,mods=~Method+Region,method="ML",data=ROM)
using a random effects model.
summary(ROM.ma1)
gives us this:
Mixed-Effects Model (k = 23; tau^2 estimator: ML)
logLik deviance AIC BIC
-4.2852 65.8950 18.5705 24.2479
tau^2 (estimated amount of residual heterogeneity): 0.0634 (SE = 0.0241)
tau (square root of estimated tau^2 value): 0.2519
I^2 (residual heterogeneity / unaccounted variability): 90.58%
H^2 (unaccounted variability / sampling variability): 10.62
Test for Residual Heterogeneity:
QE(df = 19) = 616.2226, p-val < .0001
Test of Moderators (coefficient(s) 2,3,4):
QM(df = 3) = 17.0683, p-val = 0.0007
Model Results:
estimate se zval pval ci.lb ci.ub
intrcpt -0.4392 0.3150 -1.3944 0.1632 -1.0566 0.1781
MethodRIL 0.3544 0.1513 2.3420 0.0192 0.0578 0.6510 *
RegionAmericas 0.1027 0.3201 0.3208 0.7484 -0.5247 0.7301
RegionAsia & Oceania -0.3487 0.3068 -1.1365 0.2557 -0.9500 0.2526
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Now I understand that the intercept is a combination of the first levels of the factors Method and Region.
What I would like to be able to do is calculate coefficient estimates for each of these groups along with their confidence intervals so I can plot the results of this test.
Is there a way in which I could do this?
I have asked a number of colleagues and none of them have given me a useful response.
Thanks in advance.