3
$\begingroup$

I am new to analyzing complex datasets using mixed models in R, so please forgive me if I am asking something very basic.

I'm currently trying to write a linear mixed model using the lmer function of package lme4. However, I cannot find much on the correct notation of nested fixed factors in the model. I only seem to find examples on how to include nested random effects into the model.

In my experiment, we analyzed the soil microbial diversity (Shannon) at three different distances from plant roots. I'd like to separately analyze what's happening at each distance, so I made three subsets of the data. I'll be using subset 1 as an example here. Within this subset, we had 3 different soils and sowed 5 plant individuals per soil type. The plants all have the same 5 root types (A, B, C, D, E), and we sampled 1 root of each root type per plant. We also measured the growth rate of each sampled root. We then divided each root into two sections (tip, base) and sampled the microbiome. So: Shannon is the response variable. Soil, root type, root section and growth rate are fixed factors. Plant is a random effect. I would like to examine:

  • whether the relationship between the root parameters (root type, root section, growth rate) and the response variable microbial diversity (Shannon) varies depending on the soil type (soil).
>   str(sub1)
'data.frame':   117 obs. of  6 variables:
 $ plant       : Factor w/ 14 levels "L1","L3","L4",..: 1 1 1 1 1 5 5 5 5 5 ...
 $ growth.rate : num  0.0141 0.0141 0.2425 3.161 0.7612 ...
 $ soil        : Factor w/ 3 levels "Loam","Sand",..: 1 1 1 1 1 2 2 2 2 2 ...
 $ root.type   : Factor w/ 5 levels "C","D","E","A",..: 4 4 5 1 3 4 5 5 1 1 ...
 $ root.section: Factor w/ 2 levels "base","tip": 1 2 2 1 2 2 1 2 1 2 ...
 $ Shannon     : num  4.85 4.67 4.5 4.8 3.72 ...

I attached an image here to show what I believe are multiple levels of nestedness: The random effect plant is nested in soil. Plus, if I understand it correctly, both fixed factors growth rate and root section are nested in root type.

study design

I've found online that it is possible to include nested fixed factors into an lmer by writing (nesting factor/nested factor). So I tried to fit the model:

>   M1 <- lme4::lmer(Shannon ~  soil:root.type + 
                                soil:(root.type/root.section) + 
                                soil:(root.type/growth.rate) + 
                                (1|soil/plant), 
                     data=sub1, REML=FALSE) 

>   Anova(M1, type=3)
Analysis of Deviance Table (Type III Wald chisquare tests)

Response: Shannon
                              Chisq Df Pr(>Chisq)    
(Intercept)                 834.625  1  < 2.2e-16 ***
soil:root.type               46.058 14  2.743e-05 ***
soil:root.type:root.section  83.052 14  7.653e-12 ***
soil:root.type:growth.rate   33.930 15   0.003483 **

However, the output of both summary(M1) and Anova(M1) is exactly identical to the output of M2, which I thought "only" modelled three-way interactions?:

>   M2 <- lme4::lmer(Shannon ~  soil:root.type + 
                                soil:root.type:root.section + 
                                soil:root.type:growth.rate + 
                                (1|soil/plant), 
                     data=sub1, REML=FALSE)

So I am wondering:

  1. Most importantly: Is M1 the correct notation of nested fixed factors in lmer's? (Why is the output identical to M2 and why does the Anova output look like a test for a three-way interaction? Am I fundamentally misunderstanding something?)
  2. How to continue with a posthoc test from here? So far, I've tried function emmeans() from the emmeans package but I am not sure whether it is applicable for interactions and/or nested fixed factors. Running e.g. > emmeans(M1, pairwise ~ soil : root.section) gives me an exhaustive list containing 406 p-values of all possible pairwise comparisons between soil, root section AND root type ... Is there a better way?

So sorry for the long post. Any ideas / thoughts on the above are highly appreciated!

Improve this question
$\endgroup$
2
  • 2
    $\begingroup$ Some quick thoughts. First, you don't have to specify nesting for fixed effects. See Ben Bolker's explanation. Second, I think you might be specifying soil as both a fixed and a random effect (the latter in (1|soil/plant)). As you have labeled all the plants uniquely, lmer() can figure out the nesting with (1|plant). Third, you seem to be trying to suppress intercepts by using : for interactions without "main effects," but that can lead to confusion in anything but the simplest model. Try those changes. $\endgroup$
    – EdM
    Dec 16, 2022 at 19:38
  • $\begingroup$ Thanks a lot, your quick thoughts already help me a great deal. Thanks for the link also, it is quite the relief not having to specify nesting in the model. My plants are indeed labeled uniquely, so I am happy to drop soil as a random effect. I might try lmer(Shannon ~ soil*root.type + soil*root.section + soil*growth.rate + (1|plant) now, to not exclude the main effects. $\endgroup$
    – Anis
    Dec 21, 2022 at 19:56

1 Answer 1

Reset to default
2
$\begingroup$

This is not exactly an answer, but it may help you get a bit further.

First, I seriously doubt that root.type and root.section are nested in soil. When B is nested in A, it means the the very meaning of B depends of which A we have. Here we have root.section (with levels base and tip) and soil (with levels sand, loam); are you really saying that "root base" with a sandy soil has a different meaning than "base" with a loamy soil? I doubt it; instead, I think root.section just describes which part of the root you are looking at, and that doesn't have anything to do with which soil.. These are just combinations of factors. You can have all combinations of root.section and soil; and, I suspect that's true with root.type. It doesn't even make sense that something about the root should depend on something about the soil.

I also don't like that you have no main effects at all in your model. And growth.rate is a covariate -- you measured that, right? It is not a nested factor, it's just an observation. And soil is definitely not a random effect. I think your model should be something like Shannon ~ growth.rate + soil*root.type*root.section + (1|plant) -- though if you have plants indexed so that "plant 1" occurs several times with plants that are actually different, the random term should be (1|plant:soil)

Then consider simplifying the model. Maybe you don't need the three-=way interaction. Maybe some two-factor interactions are not necessary. Doing this simplifying might make it possible to sensibly produce marginal means of one or two of the factors at a time.

The second caution I offer has to do with emmeans(). For some reason, using of pairwise ~ <everything> seems to come from some blogger's quick recipe for using emmeans(), and I wish they'd change it. Please, please, please leave out the pairwise. You don't need pairwise comparisons of a zillion cell means. Just get the means first:

EMM <- emmeans(model, ~ soil*root.type*root.section)
EMM   # display the means
emmip(EMM, root.type ~ soil|root.section)   # visualize them

Then get the contrasts you need, e.g., contrast(EMM, "pairwise", simple = "soil") or equivalently, pairs(EMM, simple = "soil")

Again, you are better off proceeding in small steps and avoid recipes that give you a torrent of unwanted information.

And if you could tell me where you got the advice to use pairwise ~, I'd be interested to know.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you a lot for your advice, especially on the correct use of the emmeans() function which makes so much more sense now! Sadly I cannot find the forum site where pairwise ~ everything was suggested anymore. Concerning the nesting: I didn't mean to nest root.type and root.section in soil; I was wondering whether root.section was nested in root.type (as a root tip from root type A is probably different from a root tip of root type C). But since EdM suggested that nesting fixed factors is overrated I guess I can stop worrying about that now? What do you think? $\endgroup$
    – Anis
    Dec 21, 2022 at 21:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.