# Coding nested factors for 3-level in rma.mv function of metafor package

I'm working collaboratively on an arm-based, network meta-analysis in the agricultural sciences. After carefully reviewing the Konstantopoulos walkthrough and Viechtbauer's slides to develop the code, I am looking for feedback on the rma.mv structure below for a multi-level model in metafor.

The data are structured such that crop yield is the dependent variable, and each arm represents a specific treatment combination of fertilizer and water (i.e. Arm 1 = all the water, fertilizer, Arm 2 = some water, fertilizer, Arm 3 = no water, fertilizer …) . The arms are nested by study, which are nested by publication. One paper can have multiple studies (i.e. same plant studied over multiple years, same plant different sites, or different types of plants studied). Note, the arms are not a crossed factor because not all arms are reported in all studies.

Paper>Study>Treatments (Arms)>Crop Yield

There are several potential moderating variables, including the specific plant studied, as well as a categorical assignment of crop types (i.e. cereal, legume, fruit, etc.).

Consistent with the linear, mixed-effects model guidance for metafor (rma.mv), I've specified the random effects as study and paper, coding the model to include arm as the inner variable (and fixed effect).

Is this code accurate based on the described structure?

fit<- rma.mv(yi, vi, mods = ~ arm + crop.type - 1,
random = ~ arm | paper/study,
data = dat,
struct = "HCS"
)


I’m happy to provide any additional information to facilitate feedback.

EDIT: Providing data sample to facilitate discussion

paper study arm  yi     vi    Crop.Type
3   1   B   -0.596  0.348   Cereal
3   1   C   -1.809  0.470   Cereal
3   1   D   -1.948  0.491   Cereal
3   1   E   -5.697  1.686   Cereal
3   1   F   -4.805  1.295   Cereal
3   1   G   -4.899  1.333   Cereal
3   1   H   -3.417  0.820   Cereal
3   2   B   0.230   0.336   Cereal
3   2   C   -0.646  0.351   Cereal
3   2   D   0.471   0.343   Cereal
3   2   E   -1.704  0.454   Cereal
3   2   F   -0.909  0.368   Cereal
3   2   G   -1.053  0.380   Cereal
3   2   H   -1.386  0.413   Cereal
6   7   B   -9.371  7.984   Cereal
6   7   F   -29.326 72.336  Cereal
6   7   G   -24.988 52.700  Cereal
7   9   B   -5.647  3.324   Oilseed
7   9   D   5.448   3.140   Oilseed
7   9   F   -14.070 17.165  Oilseed
7   9   G   -10.364 9.618   Oilseed
7   9   H   -2.267  1.095   Oilseed
7   9   J   -24.315 49.935  Oilseed
7   9   K   -19.636 32.797  Oilseed
7   9   L   -12.677 14.058  Oilseed
7   10  B   -7.357  5.178   Oilseed
7   10  D   -2.218  1.077   Oilseed
7   10  F   -19.687 32.966  Oilseed
7   10  G   -14.070 17.165  Oilseed
7   10  H   -8.948  7.339   Oilseed
7   10  J   -25.770 56.006  Oilseed
7   10  K   -22.915 44.423  Oilseed
7   10  L   -19.994 33.980  Oilseed


The rma.mv() function does not currently handle terms like random = ~ var1 | var2/var3. I suggest you install the development version of metafor (see: https://wviechtb.github.io/metafor/#installation), which will flag this as an error. You can fit this however with:

dat$$study.in.paper <- paste0(dat$$paper, ".", dat\$study)
fit <- rma.mv(yi, vi, mods = ~ arm + crop.type - 1,
random = list(~ arm | paper, ~ arm | study.in.paper),
data = dat,
struct = "HCS"
)


With two terms of the form ~ inner | outer, struct should technically have two values (for the structure of the first and the second term). However, if you only specify a single value for struct, then this automatically gets expanded to struct=c("HCS", "HCS"), which I suppose is what you had in mind anyway.

Otherwise, looks fine to me. Not sure how much data you have, but this isn't a simple model, so I hope you have a decent number of papers and a decent number of papers that have multiple studies (I know that this is vague, but I cannot give you a hard cutoff here).

Note that if you think a variables such as crop.type is a moderator (i.e., it potentially influence the difference between arms), then you should not just add it as main effect, but also as an interaction with arm.

• Thank you for providing the code. I would like to double-check a few points. 1. You create the string variable paper.study. In my dataset, study ID is unique (new data sample in question edit). Given this, can I use study as-is? i.e. random = list( ~ (arm | paper), ~ (arm | study)) 2. Since the arm is nested within study and study nested within paper, I expected (study | paper), (arm | study). Can you describe why the hierarchical relationship is represented differently in the syntax you provided?
– Eli
Apr 20 '19 at 20:12
• 3. We may wish to specify the two covariances differently. If struct = c(A, B), does A refer to the covariance structure across arms within study and B to that within paper? 4. If we decide that the moderating variables influence the arms, how would this be implemented in a way that’s currently supported by metafor? 5. I will keep your thoughts about model complexity in mind: presently, the database has 20 papers with 65 studies, containing 579 rows.
– Eli
Apr 20 '19 at 20:12
• 1. Then you can use study as is. 2. You could use that too. I just went with random = ~ arm | paper/study, which (in lme() syntax) adds random effects for arm at the paper and study levels, so I only showed how to specify this for rma.mv(). 3. A refers to whatever is your first ~ inner | outer term, B to the second. 4. For example: mods = ~ arm*crop.type. Apr 22 '19 at 20:11
• @Wolfgang, just wondering if this model can be fit using rma.mv()? Nov 4 '19 at 1:28