R 的 lm() 函数在不报告错误或警告的情况下删除因子水平(截距估计器除外)

R's lm() function removes factor levels (aside from intercept estimator) without reporting error or warning

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

响应变量 dayseclosion 为连续数值

解释变量都是分类的: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)

现在,当我说数据不平衡时,我的意思是并非所有 4host 都被收集,所有 25site 全部 5year,并且存在不平等的代表性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() table 时,它揭示了与year 被神秘删除 (year2020)。注意,这不是用于估计截距的级别(即 year2016)。那一年有 3,581 个观测值,但在 summary(model.aic.forward) table 中,系数为(仅显示部分):

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

最后(我知道这很长),方差分析 table 中 yeardf 是 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 发生了什么事?数据的不平衡性是否会导致这种情况?我无法想象如何提供一个最小的可重现示例,因为可能是导致问题的数据的数量和复杂性?

感谢您阅读并可能帮助解决这个难题。

是的,您的缺失可能导致了您的结果。

考虑以下可重现的示例。

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 因为缺失。