R 的 lm() 函数删除因子级别(除了截距估计器)而不报告错误或警告

Question

我正在尝试使用前向逐步方差分析和 AIC 选择标准，使用交互项将线性模型拟合到相当大的不平衡数据集。有 13,072 个观察值。设置方式如下:

响应变量 天数 是连续数字

解释变量都是分类的: host(4 个级别)、site(25 个级别)、year (5 级)，monoverwinter(4 级)。

> 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)

现在，当我说数据不平衡时，我的意思是并非所有 4主机 都在 5 年的所有 25站点 中收集>，并且其他因素水平中 4monoverwinter 实验处理的代表性不均等。然而，每个因子水平内仍然有大量观察结果(=数百个)。

一切似乎都运行良好 - 没有警告，也没有错误。选择以下型号:

## 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

问题是当我检查 summary() 和 anova() 表时，它们揭示了与year 被神秘删除 (year2020)。注意，这不是用于估计截距的水平(即year2016)。该年有 3,581 个观测值，但在 summary(model.aic.forward) 表中，系数为(仅显示部分):

## 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    

此处也未显示，但具体与 year2020 的所有交互也显示为 NA。

奇怪的是，根据 F-stat 自由度，似乎所有观测值(包括 year2020)都用于拟合模型 (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

最后(我知道这很长)，方差分析表中 year 的 df 是 3，但应该是 4，因为数据中有五个因子水平:

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

我的理解是错误的吗？ year2020 发生了什么？数据的不平衡性质是否会导致这种情况？我无法想象如何提供一个最小的可重现示例，因为可能是数据的数量和复杂性导致了问题？

感谢您的阅读，并且可能有助于解决这个难题。

Answer 1

是的，你的缺失可能导致了你的结果。

考虑以下可重现的示例。

x1 <- sample(1:5, 1000, replace=T)  
x2 <- sample(1:3, 1000, replace=T)
y <- 2*x1 + 3*x2 + rnorm(1000)
#no missings (everything works fine)
lm(y~as.factor(x1) + as.factor(x2))
lm(y~ -1 + as.factor(x1) + as.factor(x2)) # no intercept
#with missings
x1[x2==2]<-NA #create specific missingness
table(x1,useNA = "always")
table(x2,useNA = "always")
table(x1,x2,useNA = "always") #you see the missing pattern
lm(y~ as.factor(x1) + as.factor(x2))
lm(y~ -1 + as.factor(x1) + as.factor(x2))

如果运行代码，您将看到对于因子 x2，两个类别被省略。第一个是 ref.cat。第二个 (x2)2 因为缺失。