I have the following data: Pine Forest
Biomass ~ Age | Plot:
Each black curve represents whole-plot biomass for each individual plot
I want to formally examine this trend using mixed effects model (using
lme4 function in R).
Thus far, I created a mixed model with random intercepts and slopes:
lme4(Biomass ~ I(Age - 51) + (1 + I(Age - 51) | Plot), data = dat, REML = F)
I added additional predictors to the model to better explain some of the trend:
lmer(Biomass ~ I(Age - 51) * LossRate + I(Age - 51) * I(scale(PineStems)) + (1 + I(Age - 51)|Plot),data = dat, REML = F)
LossRateis biomass lost per year due to mortality and
PineStemsare the number of pine stems in a given plot in a given year.
This model seems to perform just fine (based on significant 95% bootstrapped CIs for coefficient estimates and for CIs for the whole model). However, it doesn't account for the downward curve you can see in some plots.
I was thinking I could account for this curve by incorporating a quadratic term for age (i.e.,
I(scale((Age - 51)^2)) [we'll call it "
However, I'm not sure the best way to do this....
I'm pretty sure that I need to add a random effect for
Age^2 since only some plots have the curve, but do I also need to add
Age^2 as a fixed effect?
Adding it as both a fixed and random effect improves the AIC of the model a lot (vs. no
Age^2term), but adding just the random
Age^2as a fixed effect) improves the model even more.
In general, does adding a squared "time" predictor as a random effect but not adding it as a fixed effect make sense?? (read: "Can I do that??")
Further, I expect the largest impact of
PineStems on the model as being the main driver causing the downward slope in biomass in mature plots. In other words, I mostly expect
PineStems to drive the
Age^2 term (so that as
PineStems decreases to a certain point, the trend line turns downward). How can I additionally add that to my model? Just simply by creating an interaction (
Age^2 * PineStems)? Or does the model make this connection on its own?