Return to Question

model <- nlme::lme(soil_property ~ Grazing*Tree_cover*Season,
               random = ~1|Random_site,
               data = df,
               na.action = na.omit,
               method = "REML")

model <- nlme::lme(soil_property ~ Grazing*Tree_cover*Season,
               random = ~1|Random_site,
               weights = varIdent(form = ~1|Tree_cover),
               data = df,
               na.action = na.omit,
               method = "REML")

model <- nlme::lme(soil_property ~ Grazing*Season*Tree_cover,
               random = ~1|Random_site,
               data = df,
               na.action = na.omit,
               method = "REML")

emmeans (model, pairwise ~ Grazing | Tree_cover)

model <- nlme::lme(soil_property ~ Grazing*Tree_cover*Season,
               random = ~1|Random_site,
               data = df,
               na.action = na.omit,
               method = "REML")

emmeans::emmeans (model, pairwise ~ Grazing | Tree_cover)