[R-meta] Calculating marginal/least squares means with two factors - rma.mv | metafor

[Viechtbauer, Wolfgang (SP)]

Dear Liam,

First, please see this thread where I just provided a similar response: https://stat.ethz.ch/pipermail/r-sig-meta-analysis/2019-January/001392.html

To elaborate on your case, there is no need for lsmeans/emmeans. If one could make these packages interface with metafor (I haven't tried), it might automate some things, but getting marginal means can all be done with existing functions in metafor.

Thank you for providing a reproducible example. Using this example:

### overall marginal mean
predict(mod, newmods = colMeans(model.matrix(mod))[-1])

### marginal means for factor A
predict(mod, newmods = c(0,0, colMeans(model.matrix(mod))[4:6])) # level B
predict(mod, newmods = c(1,0, colMeans(model.matrix(mod))[4:6])) # level C
predict(mod, newmods = c(0,1, colMeans(model.matrix(mod))[4:6])) # level D

### marginal means for factor B
predict(mod, newmods = c(colMeans(model.matrix(mod))[2:3], 0,0,0)) # level one
predict(mod, newmods = c(colMeans(model.matrix(mod))[2:3], 1,0,0)) # level two
predict(mod, newmods = c(colMeans(model.matrix(mod))[2:3], 0,1,0)) # level three
predict(mod, newmods = c(colMeans(model.matrix(mod))[2:3], 0,0,1)) # level four

Test marginal means against each other is the same thing as testing certain coefficients against each other based on the model. For example, in the output of 'mod', the coefficient for 'AC' is the test of level C vs B of factor A while accounting for factor B. For completeness sake, here are all three pairwise comparisons:

### testing marginal means against each other for factor A
anova(mod, L=c(0,  1,0, 0,0,0)) # C vs B (same as in 'mod' output)
anova(mod, L=c(0,  0,1, 0,0,0)) # D vs B (same as in 'mod' output)
anova(mod, L=c(0, -1,1, 0,0,0)) # D vs C

And for factor B:

### testing marginal means against each other for factor B
anova(mod, L=c(0, 0,0, 1,0,0))  # 2 vs 1
anova(mod, L=c(0, 0,0, 0,1,0))  # 3 vs 1
anova(mod, L=c(0, 0,0, 0,0,1))  # 4 vs 1
anova(mod, L=c(0, 0,0, 1,-1,0)) # 3 vs 2
...

Best,
Wolfgang

-----Original Message-----
From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On Behalf Of Liam Kendall
Sent: Monday, 21 January, 2019 23:43
To: r-sig-meta-analysis using r-project.org
Subject: [R-meta] Calculating marginal/least squares means with two factors - rma.mv | metafor

Dear all,

I am conducting a meta analysis with two factorial moderator/predictor
variables. I would like some help calculating marginal/least squares means
for the group combinations to see if they are significantly more or less
than zero (e.g. Factor A1 + Factor B1 > 0) as well as if they are
significantly different from one another(i.e. pairwise comparisons). I
realise I can re-level the factors to achieve my first aim but this will
not be enable me to look at pairwise comparisons.

I tried and failed to adapt the emmeans code for an rma.mv object following
this guide (
https://cran.r-project.org/web/packages/emmeans/vignettes/xtending.html) so
I was wondering if anyone has had success adapting lsmeans/emmeans for a
rma.mv / metafor model object?

Any help would be much appreciated. I have provided an example data frame
and model below as a start.

Many thanks,

Liam

#Data frame

example=structure(list(yi = c(-0.548565951748838, -0.164303051291276,
-0.0254520632036632, -0.0233556263771351, -0.0285885006937935,
0.0117786991926127, -0.000590493078864322, 0.0222231367847103,
-0.0641931576424949, -0.0164387263440607, -0.151696143283091,
-0.0826917158451134, -0.111225635110224, 0.0741079721537218,
-0.167054084663166, -0.220542769614152, -0.0689928714869514,
-0.154150679827258, -0.0100503358535015, -0.0304592074847086,
-0.139939388897495, -0.167587694187896, -0.0525401403821952,
-0.061295387688433, 0.0878002427858518, 0.0790449954796141,
-0.0463247786360805,
-0.0271157097728215, -0.0336594753512927, 0.0215062052209637),
    vi = c(3.65256775588963e-05, 2.06660766466872e-05,
1.56332665015693e-06,
    1.56417471740522e-06, 1.56206449445605e-06, 1.10168792085364e-06,
    6.23432607076493e-06, 3.96607651507311e-05, 0.00302789999884411,
    0.00206654950572419, 8.2358635771675e-07, 1.60767239099821e-07,
    0.00396607651507311, 0.00396607651507311, 7.34458613902428e-06,
    7.34458613902428e-06, 7.34458613902428e-06, 7.34458613902428e-06,
    2.20337584170728e-05, 2.20337584170728e-05, 2.62914175316163e-07,
    2.31792109687421e-07, 2.38346536711187e-06, 2.26344954289566e-06,
    1.85676190504622e-06, 1.73674608083001e-06, 1.23468002879564e-06,
    1.18221030259889e-06, 5.29878610033144e-07, 2.21377164546603e-07
    ), A = structure(c(1L, 3L, 2L, 2L, 2L, 2L, 1L, 3L, 2L, 2L,
    2L, 2L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,
    1L, 2L, 2L, 3L, 3L), .Label = c("B", "C", "D"), class = "factor"),
    B = structure(c(3L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 3L, 1L, 1L,
    1L, 4L, 4L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L,
    3L, 3L, 1L, 1L), .Label = c("one", "two", "three", "four"
    ), class = "factor")), row.names = c(1L, 2L, 4L, 5L, 6L,
7L, 9L, 10L, 12L, 13L, 16L, 17L, 22L, 23L, 32L, 33L, 34L, 35L,
36L, 37L, 38L, 39L, 42L, 43L, 44L, 45L, 48L, 49L, 56L, 58L), class =
"data.frame")

#Model

mod=rma.mv(yi,vi,mods=~A+B,data=example)