Help for possible nested mixed effect model

Question

I'm super new to mixed effect models and I wanted to make sure I was interpreting R code correctly. I'm using the "lmer" function in the "lme4" R package to do my analyses.

I'm interested in physiological impacts of nitrogen and sulfur addition in tree species that have different mycorrhizal associations across three different sites. The issue that I'm struggling with figuring out is how to deal with species and mycorrhizal status, since a particular tree species usually only has one mycorrhizal type. Here's what I have for my code thus far. We've placed nitrogen treatments, sulfur treatments, mycorrhizal status, and species as fixed effects and site as the sole random effect:

model1 <- with(.data, 
               lmer(phys_param ~ nitrogen_treat * sulfur_treat * 
                      myc_status * species_code + (1 | site)))

Here is what my dataset looks like (it's been randomized to maintain anonymity):

  site nitrogen_treat sulfur_treat species_code myc_status   phys_param
1  BH1       nitrogen       sulfur           AB         EM         9.48
2  BH1       nitrogen       sulfur           SM         AM         3.70
3  BH1       nitrogen       sulfur           SM         AM         7.39
4  BH1         no_nit       sulfur           SM         AM         1.08
5  BH1       nitrogen       sulfur           SM         AM         8.74
6  BH1         no_nit    no_sulfur           SM         AM         9.95

As each species code has a single mycorrhizal status, I'm wondering if I need to nest myc_status within species_code and if so, how this could be represented. I'm also open to any suggestions to improve the general model.

Answer 1

I recently wrote a thread on Twitter about nested fixed effects in the context of lm() models which may come in handy to you: https://mobile.twitter.com/isabellaghement/status/1169752747375116289. See also Ben Bolker's answer on this thread: How to model nested fixed-factor with GLMM, which raises the legitmate point that you would need a lot of data to detect the presence of complex interactions.

Assuming you have sufficient data, if nitrogen_treat, sulfur_treat and species_code are fully crossed factors and myc_status is nested within species_code, then it seems to me that you should be able to fit this model:

model1 <- with(.data, 
               lmer(phys_param ~ nitrogen_treat * sulfur_treat * (species_code + myc_status) + (1 | site)))

This model is basically a collection of sub-models, one sub-model for each combination of levels of the four factors which is present in your data. Because the collection is smaller than what you would have gotten with a fully crossed design for the four factors, you'll be limited in the type of inferences you can make based on your model. To understand what kind of inferences you can make, it would help to clarify what specific combinations of levels of the four factors are actually present in your data.