How to deal with nestedness of fixed factors in a linear mixed effects model (lmer in lme4)? 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.

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:

*

*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?)

*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!
 A: 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.
