# Help for possible nested mixed effect model

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.

model1 <- with(.data,