Using mixed models, is there a way to explicitly ignore one pair of crossed fixed factors in a partially-crossed experiment?

Question

The data and model

I am working with data from a partially-crossed repeated measures experimental design. I would like to use mixed models to evaluate how my dependent variable, plant biomass, is affected by a three way interaction (or a simpler combination) between fixed effects - one continuous variable and two factors:

• time ("yr", continuous measure of time since treatment began)
• fire regime (3 levels: "con", "ex", and "re", where "con" is the control)
• fertilizer (2 levels: "none" and "fert", where "none" is the control)

This data set has nested random effects based on the spatial location of the samples: One sample was collected randomly from within each stratified subplot area ("subplot_id") each year. 5-7 subplots were located within an individual "plot".

> my.df
# A tibble: 312 x 6
      yr  plot subplot_id fire_trt fert_trt biomass
   <dbl> <dbl> <chr>      <fct>    <fct>      <dbl>
 1     0    41 41-1       ex       none      1.36  
 2     0    41 41-2       ex       none      0.542 
 3     0    41 41-3       ex       none      1.09  
 4     0    41 41-4       ex       none      2.02  
 5     0    41 41-5       ex       none      0.691 
 6     0    42 42-2       ex       fert      0.864 
 7     0    42 42-4       ex       fert      0.598 
 8     0    42 42-5       ex       fert      0.253 
 9     0    42 42-6       ex       fert      0.0383
10     0    42 42-7       ex       fert      0.559 
# ... with 302 more rows

Sample data are at the bottom of this post. In reality I will be working with multiple data sets from different sites & types of biomass, each similarly structured (& similarly unbalanced) to the sample data.

Here's my current model (the 3-way interaction is indeed significant for this particular data set):

mod <- lmer(biomass ~ fire_trt*fert_trt*yr + (1|plot) + 
             (1|plot:subplot_id), data = my.df)

#fixed-effect model matrix is rank deficient so dropping 
#2 columns / coefficients

 summary(mod)$varcor
 Groups          Name        Std.Dev.
 plot:subplot_id (Intercept) 0.29966 
 plot            (Intercept) 0.11822 
 Residual                    0.51916 

 summary(mod)$coefficients
                              Estimate Std. Error        df    t value     Pr(>|t|)
(Intercept)                 1.74306791 0.11721606  15.12193 14.8705548 1.956683e-10
fire_trtex                 -1.00366818 0.20302418  15.12193 -4.9435894 1.725986e-04
fire_trtre                 -0.83597681 0.20302418  15.12193 -4.1176219 8.988053e-04
fert_trtfert               -0.22440504 0.21400628  15.12193 -1.0485909 3.108374e-01
yr                          0.09798873 0.02721118 227.71514  3.6010465 3.890172e-04
fire_trtex:fert_trtfert    -0.16962079 0.35810099  15.12193 -0.4736675 6.425031e-01
fire_trtex:yr              -0.16170996 0.04713115 227.71514 -3.4310636 7.138130e-04
fire_trtre:yr               0.03139258 0.04713115 227.71514  0.6660687 5.060414e-01
fert_trtfert:yr            -0.21785079 0.04968059 227.71514 -4.3850283 1.772429e-05
fire_trtex:fert_trtfert:yr  0.21953446 0.08313153 227.71514  2.6408086 8.843392e-03

The problem

All models that include an interaction effect between fert_trt and fire_trt factors are rank deficient (where fire_trt = "con", fert_trt = "none" by design, leading to colinearity in the model matrix). The rank deficiency is a concern for a few reasons:

• The coefficient estimate of particular interest, fire_trtre:fert_trtfert:yr is dropped, presumably due to the rank deficiency.
• I'd like to also change the covariance structure to account for temporal autocorrelation, but nlme won't run rank deficient models.
• Most importantly, we're not interested in this interaction of fertilizer on control plots in the first place.

We are interested in the interaction of fertilizer treatment with fire treatments "re" and "ex". But we want to keep control plot data in the model because it will allow us to see the "background" levels of change in the response variable over time and assess whether the slopes, beginnings, and end points of treated plots are significantly different from control (through emmeans).

Some possible solutions?

I have found a few related-but-unhelpful posts and perhaps some potential solutions. I'd really appreciate some wisdom on what's most appropriate here. Some options:

1. Stick with the current covariance structure & rank deficient model, and just acknowledge this in our results.

2. Somehow specify or pre-define which columns in the model.matrix to drop/ignore so we only compare the effects of interest? Similar to the question here, is it okay to explicitly remove the colinear lower-level fixed effects from my model? For example, the following model is rank deficient but drops the two columns I'm NOT interested in: "fire_trtcon:fert_trtfert" and "fire_trtcon:fert_trtfert:yr".

mod2 <- lmer(biomass ~ fire_trt*fert_trt*yr - fert_trt - fert_trt*yr +
(1|plot) + (1|plot:subplot_id), data = my.df)

> summary(mod2)$coefficients
                               Estimate Std. Error        df    t value     Pr(>|t|)
(Intercept)                 1.743067912 0.11721606  15.12193 14.8705548 1.956683e-10
fire_trtex                 -1.003668176 0.20302418  15.12193 -4.9435894 1.725986e-04
fire_trtre                 -0.835976809 0.20302418  15.12193 -4.1176219 8.988053e-04
fire_trtex:fert_trtfert    -0.394025829 0.28711955  15.12193 -1.3723407 1.899716e-01
fire_trtre:fert_trtfert    -0.224405037 0.21400628  15.12193 -1.0485909 3.108374e-01
fire_trtcon:yr              0.097988725 0.02721118 227.71514  3.6010465 3.890172e-04
fire_trtex:yr              -0.063721235 0.03848242 227.71514 -1.6558531 9.912893e-02
fire_trtre:yr               0.129381305 0.03848242 227.71514  3.3620886 9.072242e-04
fire_trtex:fert_trtfert:yr  0.001683661 0.06665351 227.71514  0.0252599 9.798698e-01
fire_trtre:fert_trtfert:yr -0.217850795 0.04968059 227.71514 -4.3850283 1.772429e-05

3. Restructure my data by creating a single factor representing the interaction, such that the two factors fire_trt and fert_trt are treated as 1 factor with five different levels. But won't this change the lower order interactions and main effects, as well my ability to directly compare AIC values to simpler versions (say, biomass ~ fire_trt*time + fert_trt ...)?

4. Restructure my data so some of the factors are hierarchical/nested? Maybe something like where a 2-level factor "trtmt" (with levels "con" vs "trtd") contains another 2-level factor "fire_trt" (with levels "re" & "ex") within the "trtd" factor level. And then just looking at the interaction between fertilizer treatment and the two levels of "treated" plots. But I'm not sure if this is a stretch, conceptually, to nest the data in this way. And I'm not sure this would be better. My attempt to run a model with the data modified as described also produced rank deficient results.

my.df2 <- my.df %>% 
     mutate(trtmt = if_else(fire_trt == "con", "con", "trtd"))

lmer(biomass ~ trtmt:fire_trt*fert_trt*yr + (1|plot/subplot_id), 
     data = my.df2)

Any guidance would be appreciated!

Sample data

my.df <- structure(list(yr = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
5, 5, 5, 5, 5), plot = c(41, 41, 41, 41, 41, 42, 42, 42, 42, 
42, 43, 43, 43, 43, 43, 47, 47, 47, 47, 47, 46, 46, 46, 46, 46, 
48, 48, 48, 48, 48, 44, 44, 44, 44, 44, 45, 45, 45, 45, 45, 17, 
17, 17, 17, 17, 21, 21, 21, 21, 21, 25, 25, 25, 25, 25, 29, 29, 
29, 29, 29, 47, 47, 47, 47, 47, 47, 47, 25, 25, 25, 25, 25, 25, 
25, 44, 44, 44, 44, 44, 44, 44, 46, 46, 46, 46, 46, 46, 46, 43, 
43, 43, 43, 43, 43, 43, 45, 45, 45, 45, 45, 45, 45, 41, 41, 41, 
41, 41, 41, 41, 42, 42, 42, 42, 42, 42, 42, 48, 48, 48, 48, 48, 
48, 48, 17, 17, 17, 17, 17, 17, 17, 29, 29, 29, 29, 29, 29, 29, 
21, 21, 21, 21, 21, 21, 21, 17, 17, 17, 17, 17, 17, 17, 21, 21, 
21, 21, 21, 21, 21, 25, 25, 25, 25, 25, 25, 25, 29, 29, 29, 29, 
29, 29, 29, 41, 41, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 42, 
42, 43, 43, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 44, 44, 45, 
45, 45, 45, 45, 45, 45, 48, 48, 48, 48, 48, 48, 48, 47, 47, 47, 
47, 47, 47, 47, 46, 46, 46, 46, 46, 46, 46, 17, 17, 17, 17, 17, 
17, 17, 21, 21, 21, 21, 21, 21, 21, 25, 25, 25, 25, 25, 25, 25, 
29, 29, 29, 29, 29, 29, 29, 43, 43, 43, 43, 43, 43, 43, 44, 44, 
44, 44, 44, 44, 44, 41, 41, 41, 41, 41, 41, 41, 42, 42, 42, 42, 
42, 42, 42, 48, 48, 48, 48, 48, 48, 48, 45, 45, 45, 45, 45, 45, 
45, 47, 47, 47, 47, 47, 47, 47, 46, 46, 46, 46, 46, 46, 46), 
    subplot_id = c("41-1", "41-2", "41-3", "41-4", "41-5", "42-2", 
    "42-4", "42-5", "42-6", "42-7", "43-1", "43-2", "43-3", "43-6", 
    "43-7", "47-1", "47-2", "47-3", "47-4", "47-5", "46-1", "46-3", 
    "46-4", "46-6", "46-7", "48-1", "48-2", "48-4", "48-6", "48-7", 
    "44-2", "44-3", "44-4", "44-5", "44-7", "45-1", "45-2", "45-3", 
    "45-6", "45-7", "17-1", "17-2", "17-3", "17-4", "17-7", "21-1", 
    "21-2", "21-4", "21-5", "21-7", "25-1", "25-2", "25-3", "25-5", 
    "25-6", "29-1", "29-2", "29-3", "29-4", "29-6", "47-1", "47-2", 
    "47-3", "47-4", "47-5", "47-6", "47-7", "25-1", "25-2", "25-3", 
    "25-4", "25-5", "25-6", "25-7", "44-1", "44-2", "44-3", "44-4", 
    "44-5", "44-6", "44-7", "46-1", "46-2", "46-3", "46-4", "46-5", 
    "46-6", "46-7", "43-1", "43-2", "43-3", "43-4", "43-5", "43-6", 
    "43-7", "45-1", "45-2", "45-3", "45-4", "45-5", "45-6", "45-7", 
    "41-1", "41-2", "41-3", "41-4", "41-5", "41-6", "41-7", "42-1", 
    "42-2", "42-3", "42-4", "42-5", "42-6", "42-7", "48-1", "48-2", 
    "48-3", "48-4", "48-5", "48-6", "48-7", "17-1", "17-2", "17-3", 
    "17-4", "17-5", "17-6", "17-7", "29-1", "29-2", "29-3", "29-4", 
    "29-5", "29-6", "29-7", "21-1", "21-2", "21-3", "21-4", "21-5", 
    "21-6", "21-7", "17-1", "17-2", "17-3", "17-4", "17-5", "17-6", 
    "17-7", "21-1", "21-2", "21-3", "21-4", "21-5", "21-6", "21-7", 
    "25-1", "25-2", "25-3", "25-4", "25-5", "25-6", "25-7", "29-1", 
    "29-2", "29-3", "29-4", "29-5", "29-6", "29-7", "41-1", "41-2", 
    "41-3", "41-4", "41-5", "41-6", "41-7", "42-1", "42-2", "42-3", 
    "42-4", "42-5", "42-6", "42-7", "43-1", "43-2", "43-3", "43-4", 
    "43-5", "43-6", "43-7", "44-1", "44-2", "44-3", "44-4", "44-5", 
    "44-6", "44-7", "45-1", "45-2", "45-3", "45-4", "45-5", "45-6", 
    "45-7", "48-1", "48-2", "48-3", "48-4", "48-5", "48-6", "48-7", 
    "47-1", "47-2", "47-3", "47-4", "47-5", "47-6", "47-7", "46-1", 
    "46-2", "46-3", "46-4", "46-5", "46-6", "46-7", "17-1", "17-2", 
    "17-3", "17-4", "17-5", "17-6", "17-7", "21-1", "21-2", "21-3", 
    "21-4", "21-5", "21-6", "21-7", "25-1", "25-2", "25-3", "25-4", 
    "25-5", "25-6", "25-7", "29-1", "29-2", "29-3", "29-4", "29-5", 
    "29-6", "29-7", "43-1", "43-2", "43-3", "43-4", "43-5", "43-6", 
    "43-7", "44-1", "44-2", "44-3", "44-4", "44-5", "44-6", "44-7", 
    "41-1", "41-2", "41-3", "41-4", "41-5", "41-6", "41-7", "42-1", 
    "42-2", "42-3", "42-4", "42-5", "42-6", "42-7", "48-1", "48-2", 
    "48-3", "48-4", "48-5", "48-6", "48-7", "45-1", "45-2", "45-3", 
    "45-4", "45-5", "45-6", "45-7", "47-1", "47-2", "47-3", "47-4", 
    "47-5", "47-6", "47-7", "46-1", "46-2", "46-3", "46-4", "46-5", 
    "46-6", "46-7"), fire_trt = structure(c(2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("con", "ex", "re"
    ), class = "factor"), fert_trt = structure(c(1L, 1L, 1L, 
    1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("none", "fert"
    ), class = "factor"), biomass = c(1.356764544, 0.542272916, 
    1.089511468, 2.015284272, 0.691496644, 0.864402636, 0.597531532, 
    0.253247436, 0.038324524, 0.558570388, 0.368348332, 0.231984328, 
    1.099697388, 0.067991016, 0.128851888, 0.626052108, 1.431503732, 
    0.214668264, 0.29539168, 0.517699384, 1.171762772, 1.008533404, 
    0.153298096, 0.13241696, 0.908202092, 1.287627612, 1.216326172, 
    1.101861896, 1.608229444, 1.106190912, 3.415084328, 0.183473884, 
    0.423097652, 0.174561204, 1.319585936, 1.211233212, 0.096511592, 
    0.389356792, 0.203591076, 0.9612962, 2.019231316, 1.215434904, 
    1.9703389, 1.982180032, 2.268404384, 1.624272268, 3.799730132, 
    1.494656436, 2.263311424, 1.561119564, 2.492876596, 0.888848844, 
    0.894323776, 1.450093036, 0.807106836, 1.904512392, 0.712759752, 
    2.25809114, 1.185895736, 1.779607548, 0.40552694, 1.477340372, 
    0.594475756, 1.069139628, 1.768275712, 2.07920092, 0.807616132, 
    1.60300916, 1.329007912, 1.672528064, 3.053102196, 0.616375484, 
    2.01681216, 1.95633326, 0.347976492, 0.00152788799999999, 
    0.655845924, 0.561626164, 0.425389484, 0.634964788, 0.4647326, 
    1.026358764, 1.125798808, 1.151900228, 0.212885728, 0.512861072, 
    0.126432732, 0.283423224, 1.254396048, 0.069773552, 0.336262684, 
    0.360454244, 0.544564748, 0.060733548, 0.311179856, 0.261268848, 
    0.960277608, 0.906292232, 0.019862544, 0.2355494, 0.802395848, 
    0.769928228, 0.707157496, 0.260632228, 0.954802676, 0.334480148, 
    1.141332336, 0.649734372, 0.695571012, 0.00738479199999999, 
    0.112299768, 0.449581044, 0.061242844, 0.135090764, 0.035141424, 
    0.56913828, 1.376754412, 0.412784408, 0.276547728, 0.677872976, 
    0.582889272, 1.269929576, 0.659156348, 2.0626488, 1.922847048, 
    1.939908464, 1.5756345, 1.713017096, 1.752487536, 1.632802976, 
    2.283173968, 2.226005492, 2.932526368, 1.826080808, 1.745484716, 
    1.65075566, 1.087601608, 1.428829928, 2.702324576, 2.972760752, 
    2.147064612, 1.838303912, 1.308508748, 1.532598988, 2.980272868, 
    2.486255748, 2.471613488, 3.125549552, 1.38019216, 1.32034988, 
    1.241154352, 1.801379952, 1.758853736, 3.227790724, 3.64910584, 
    2.782538696, 1.931250432, 1.396871604, 1.987909612, 1.125416836, 
    1.00649622, 1.370770184, 1.746757956, 3.024963592, 2.12885728, 
    2.94564074, 2.688318936, 1.754779368, 2.94054778, 0.930611116, 
    1.433413592, 1.238989844, 1.109883308, 0.287370268, 1.924374936, 
    1.03514412, 1.011207208, 0.701045944, 1.063282724, 0.022790996, 
    0.000763944, 0.2578311, 0.38324524, 0.016806768, 0.02737466, 
    0.140565696, 0.39088468, 0.436339348, 0.13241696, 0.232111652, 
    1.740264432, 0.127069352, 0.113191036, 0.598677448, 0.434047516, 
    0.690860024, 0.035268748, 0.326331412, 0.188057548, 0.53730728, 
    0.076012428, 0.270308852, 0.27565646, 0.044945372, 0.127324, 
    0.47237204, 0.695952984, 1.749559084, 1.099570064, 0.035523396, 
    0.48892416, 1.687806944, 1.957988472, 1.783427268, 1.13063712, 
    1.293866488, 0.724346236, 0.420423848, 0.889740112, 3.57334806, 
    2.141462356, 0.34313818, 0.709067356, 0.110262584, 0.47046218, 
    0.149860348, 0.356252552, 0.22345362, 2.305073696, 1.465881212, 
    2.222695068, 2.4032405, 2.117270796, 1.533108284, 1.566212524, 
    2.769424324, 2.32875596, 2.087222332, 3.833598316, 4.686287144, 
    1.534763496, 0.773238652, 1.607720148, 1.186914328, 3.63637344, 
    2.807876172, 1.23758928, 2.99720696, 1.252486188, 2.263693396, 
    2.194047168, 2.792724616, 3.060741636, 1.92832198, 1.706650896, 
    1.525596168, 0.168449652, 0.015024232, 0.11586484, 0.275529136, 
    0.1241409, 0.17761698, 0.093837788, 1.162213472, 0.051311572, 
    0.417240748, 0.18971276, 0.270308852, 0.263178708, 0.13241696, 
    0.394322428, 0.63407352, 0.469698236, 0.456838512, 0.514516284, 
    0.512988396, 0.716324824, 0.193150508, 0.006748172, 0.095365676, 
    0.096893564, 0.040361708, 0.06429862, 0.031449028, 2.92717876, 
    0.993891144, 0.456583864, 1.401455268, 1.271075492, 2.122109108, 
    1.71378104, 0.151897532, 0.089508772, 0.018716628, 0.12668738, 
    0.023682264, 0.002419156, 0.195569664, 0.91927928, 1.742938236, 
    0.396486936, 0.933412244, 2.928324676, 1.470974172, 2.183097304, 
    0.59714956, 0.681056076, 0.741407652, 0.538453196, 0.15215218, 
    0.082887924, 0.333716204)), row.names = c(NA, -312L), class = c("tbl_df", 
"tbl", "data.frame"))

Answer 1

This is just to make the 5-level-factor approach more explicit. What you have is an empty cell:

> with(my.df, table(fert_trt, fire_trt))
        fire_trt
fert_trt con  ex  re
    none 104  52  52
    fert   0  26  78

Let's define a 5-level factor trt having levels "con", "ex.n", "re.n", "ex.f", and "re.f", and fit your model(s) with the fixed-effects part being ~ yr * trt. (Note: explicitly list the levels in factor() so they come out in the order shown rather than alphabetically.)

Later in doing post hoc analyses, use emmeans::emmeans(mod, ~ trt) to study marginal means, averaged over yr, and emmeans::emtrends(mod, ~ trt, var = "yr") to estimate marginal slopes. Each yields a set of 5 estimates. Then you can use emmeans::contrast() to study key contrasts; most importantly,

key.contrasts <- list(
  fire_eff = c(0, 1, -1, 1, -1) / 2,
  fert_eff = c(0, 1, 1, -1, -1) / 2,
  fire:fert_eff = c(0, 1, -1, -1, 1) / 2,
  con_eff = c(4, -1, -1, -1, -1) / 4)

The first three of these compare the column, row, and diagonal averages of the 2x2 table where we ignore the con column, and the last compares the (1,1) cell with the others.

I wouldn't worry about interpreting the regression coefficients (even though this is much simpler for this model than the 3-factor model). As long as you can estimate the effects you need, the coefficients don't need interpretation.

Answer 2

Your interactions of one 3-level categorical predictor with one 2-level categorical predictor (one combination of categories missing) and with a continuous predictor limits you to an intercept plus 9 coefficients. That has nothing to do with this being a mixed model. However you end up coding this, the results should all be identical in terms of outcome predictions at any set of predictor values when all of the interaction terms are taken into account. For example:

predict(mod,data.frame(fire_trt="re",fert_trt="fert",yr=3,subplot_id="41-1",plot="41"))
        1 
0.6235175 

predict(mod2,data.frame(fire_trt="re",fert_trt="fert",yr=3,subplot_id="41-1",plot="41"))
        1 
0.6235175

With respect to Option 3, creating a single 5-level categorical predictor, you say:

won't this change the lower order interactions and main effects...?

When there are interactions you should be wary of how you interpret regression coefficients for "main effects" or lower order interactions anyway. Those are the values only when all of the interacting predictors at higher levels are at 0 (continuous predictors) or at reference levels (categorical predictors). In your case, re-centering the yr values would change the intercept and the coefficients and some lower-level interactions (those that omit yr) of the other predictors. Try this to see:

mod3 <- lmer(biomass ~ fire_trt*fert_trt*I(yr-2.5) + (1|plot) (1|plot:subplot_id), data = my.df)

So don't be put off by that fact alone.

I would lean toward Option 3, a 5-level categorical predictor, with post hoc tests of particular coefficient combinations of interest. Russ Lenth illustrates that approach with your data in another answer.

I don't see a harm in using mod2 if you wish. Be very careful in how you do your model comparisons, however. Your mod2 replaces the yr coefficient in mod1 with a numerically identical fire_trtcon:yr interaction coefficient. It might be safest to specify each of the "main effects" and interaction terms as defined in mod2 explicitly. That way, if you try to compare models, you can specify what models to compare rather than depend on how R happens to remove columns to fix rank deficiency.

Finally, if you want to compare mixed models having different fixed effects via likelihood-based metrics (likelihood-ratio tests on nested models, or the AIC as you indicate in the question), you should not be using the default REML option for fitting the models you compare. See this page. As Russ Lenth notes in a comment, once model comparisons are done it would be best to fit the final model with REML for further analysis.