I have a vegetation data set that consists of 150 plots that were sampled 1-3 times over a three year period. Plots are my unit of observation and they are unbalanced (since plots were sampled either once, twice, or three times). I would like to use mixed-effects models in order to account for variation in both plots and sampling year and to keep my sample size large (instead of conducting my analysis within individual years).
My response variables are cover of vegetation functional groups and predictors include variables related to fire and treatment history. Additionally, I am not interested in how plots change over time per se, but rather, in aggregating sampling from all three years to increase my sample size and to account for the spatial/temporal correlation that arises from doing so. It is my assumption that treating plot as a random effect (intercept only) accounts for variation that arises from potential spatial autocorrelation, but my main question is how to account for the repeated measures and if I need to account for the grouping of cells within sampling years:
model <- lmer(response~covariates + (1|Plot) + (1|Year).
However, I know that is not appropriate to use a random effect with only three levels, year in this case. I'm hoping for recommendations on how to incorporate year as a random effect. Is including (Year | Plot) recommended? And if so, how might I interpret that effect, i.e., is it accounting for variation introduced by different sampling year or variation in plots over sampling year?