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
species_code and if so, how this could be represented. I'm also open to any suggestions to improve the general model.