I am attempting to fit a linear model, using forward stepwise ANOVA and AIC selection criterion, to a fairly large unbalanced data set, with interaction terms. There are 13,072 observations. Here's how it's set up:
Response variable dayseclosion is continuous numeric
Explanatory variables are all categorical: host (4 levels), site (25 levels), year (5 levels), monoverwinter (4 levels).
> glimpse(dat)
Observations: 13,072
Variables: 5
$ site <chr> "10", "10", "10", "10", "10", "14", "14", "15", "15",…
$ year <chr> "2014", "2014", "2014", "2014", "2014", "2017", "2017…
$ host <chr> "3", "3", "4", "3", "3", "1", "1", "1", "1", "1", "1"…
$ monoverwinter <chr> "6", "6", "6", "6", "6", "5", "5", "5", "5", "5", "5"…
$ dayseclose <dbl> 11, 12, 17, 19, 19, 20, 20, 20, 20, 21, 21, 21, 21, 2…
m0 –> lm(dayseclose ~ 1, data = dat)
model.aic.forward –> step(m0, direction = "forward", trace = 1, scope = ~ host * monoverwinter * site * year)
Now, when I say the data are unbalanced, I mean that not all 4host were collected at all 25site in all 5year, and there was unequal representation of 4monoverwinter experimental treatments among other factor levels. However, there are still large numbers of observations (=hundreds) within each factor level.
Everything seems to run fine—no Warnings and no Errors. The following model is selected:
## Step: AIC=61561.31
## dayseclose ~ host + site + monoverwinter + year + host:site +
## host:monoverwinter + site:monoverwinter + site:year
##
## Df Sum of Sq RSS AIC
## <none> 1433157 61561
## + host:year 1 8.059 1433149 61563
## + host:monoverwinter:site 3 242.885 1432914 61565
The problem is when I inspect the summary() and anova() tables, which reveal that one of the factor levels associated with year is being mysteriously removed (year2020). NOTE, this is not the level being used to estimate the intercept (that is year2016). There are 3,581 observations for that year, yet in the summary(model.aic.forward) table, the coefficients are (only partial shown):
## year2017 6.36787 2.44775 2.602 0.009292 **
## year2018 -0.13757 1.85568 -0.074 0.940906
## year2019 -10.56667 3.45693 -3.057 0.002243 **
## year2020 NA NA NA NA
Also not shown here, but all interactions with year2020 specifically, also appear NA.
Oddly, based on the F-stat degrees of freedom it appears all observations, including year2020 are being used to fit the model (79 + 12992 = 13071):
## Residual standard error: 10.5 on 12992 degrees of freedom
## Multiple R-squared: 0.3105, Adjusted R-squared: 0.3063
## F-statistic: 74.06 on 79 and 12992 DF, p-value: < 2.2e-16
Finally (I know this is long), the df for year in the anova table is 3 but should be 4, given five factor levels in the data:
d =anova(model.aic.forward)
as_tibble(d, rownames = "Predictors")
## # A tibble: 9 x 6
## Predictors Df `Sum Sq` `Mean Sq` `F value` `Pr(>F)`
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 host 3 497707. 165902. 1504. 0.
## 2 site 24 43276. 1803. 16.3 3.77e-67
## 3 monoverwinter 3 34395. 11465. 104. 1.73e-66
## 4 year 3 2422. 807. 7.32 6.71e- 5
## 5 host:site 16 43445. 2715. 24.6 1.08e-72
## 6 host:monoverwinter 3 7319. 2440. 22.1 2.80e-14
## 7 site:monoverwinter 12 9229. 769. 6.97 9.13e-13
## 8 site:year 15 7646. 510. 4.62 6.30e- 9
## 9 Residuals 12992 1433157. 110. NA NA
Am I interpreting this wrong? What is happening to year2020? Could the unbalanced nature of the data be causing this? I can't imagine how to provide a minimal reproducible example, since it is likely the amount and complexity of the data that is causing the issue?
Thank you for reading and potentially help solve this puzzle.