使用虚拟变量交互项回归时的 NA 值答案

【问题标题】：NA values when regressing with dummy variable interaction term使用虚拟变量交互项回归时的 NA 值
【发布时间】：2017-12-26 08:56:00
【问题描述】：

我正在尝试估算影响纽约和芝加哥居民幸福度差异的因素。

数据如下所示。

  Happiness     City Gender Employment   Worktype      Holiday
1        60 New York      0        0     Unemployed   Unemployed
2        80  Chicago      1        1     Whitecolor 1 day a week
3        39  Chicago      0        0     Unemployed   Unemployed
4        40 New York      1        0     Unemployed   Unemployed
5        69  Chicago      1        1     Bluecolor  2 day a week
6        90  Chicago      1        1     Bluecolor  2 day a week
7       100 New York      0        1     Whitecolor 2 day a week
8        30 New York      1        1     Whitecolor 1 day a week

幸福程度是因变量，“城市”是人们居住的地方。 “性别”编码为 0 = 男人 1 = 女人。 “就业”是 0 = 失业和 1 = 就业。 “工作类型”是三个级别的因素：“失业”、“白色”、“蓝色”。 “假期”是一个人一周休息多少天。这里的“城市”、“性别”、“工作类型”和“假期”变量都是因素。 'Happiness' 和 'Employment' 变量类型是数值型的。

我要估计的模型是

lm(Happiness ~ City + Gender + Employment:(Worktype + Holiday))

我将“Employment”值保留为数值，因此如果“Employment”等于 0(Unemployment)，0:(Worktype + Holiday) = 0，则模型自动缩减为

lm(Happiness ~ City + Gender)

适用于失业人员。

但是，回归结果返回 NA 值。

Coefficients: (2 not defined because of singularities)
                               Estimate Std. Error t value Pr(>|t|)
(Intercept)                       56.75      23.56   2.408    0.138
CityNew York                     -14.50      27.21  -0.533    0.647
Gender1                           -2.25      35.99  -0.063    0.956
Employment:WorktypeBluecolor      25.00      43.02   0.581    0.620
Employment:WorktypeUnemployed        NA         NA      NA       NA
Employment:WorktypeWhitecolor     57.75      35.99   1.604    0.250
Employment:Holiday1 day a week   -50.00      54.42  -0.919    0.455
Employment:Holiday2 day a week       NA         NA      NA       NA

这似乎是由于“Worktype”和“Holiday”变量中的“Unemployment”值造成的。但是，我不确定为什么 R 不将就业：WorktypeUnemployed 显然是 0:Worktype = 0 视为零并且不将其从模型中删除。这是因为 R 将就业：假日失业设置为基线，并且两者都是完全多重共线性的吗？（我不得不为“工作类型”和“假期”输入“失业”值，因为我想看看“工作类型”和“假期”与“失业”人相比的效果。如果我删除“失业”值，NA 会消失，但是基线将是“Whitecolor”和“每周 1 天”，因此与“失业”相比，我看不到效果。）

如果是这样，为什么我的“Employement:Holiday2 day a week”的系数会得到 NA？似乎与“失业”价值无关。

我可以在仅删除 NA 系数的同时依赖此结果吗？

以下是可重现的代码。

Happiness <- c(60, 80, 39, 40, 69, 90, 100, 30)

City <- as.factor(c("New York", "Chicago", "Chicago", "New York", "Chicago",         
                  "Chicago", "New York", "New York"))
Gender <- as.factor(c(0, 1, 0, 1, 1, 1, 0, 1)) # 0 = man, 1 = woman.
Employment <- c(0,1, 0, 0, 1 ,1 , 1 , 1) # 0 = unemployed, 1 = employed.
Worktype <- as.factor(c("Unemployed", "Whitecolor", "Unemployed",     
          "Unemployed", "Bluecolor", "Bluecolor", "Whitecolor","Whitecolor"))
Holiday <- as.factor(c(0, 1, 0, 0, 2, 2, 2, 1))
levels(Holiday) <- c("Unemployed", "1 day a week", "2 day a week")

data <- data.frame(Happiness, City, Gender, Employment, Worktype, Holiday)

head(data,8)
str(data)

reg <- lm(Happiness ~ City + Gender + Employment:(Worktype + Holiday))
summary(reg)

【问题讨论】：

标签： r interaction dummy-variable

【解决方案1】：

您不必担心Employment:WorktypeUnemployed 的 NA 值。 R 尝试自动计算所有交互作用，但该特定系数仍未确定，因为很明显，Employment=1 和 Worktype="Unemployed" 绝不是这种情况。它对其他系数的计算没有任何影响：您可以通过手动编码虚拟变量来验证：

> library(lme4) # for the convenient "dummy" function 
> data <- data.frame(data, 
+   dummy(Worktype, c("Bluecolor","Whitecolor")), 
+   h1=dummy(Holiday)[,1], 
+   h2=dummy(Holiday)[,2])
>   
> reg <- lm(Happiness ~ City + Gender + Employment:Bluecolor + Employment:Whitecolor  + Employment:h1 + Employment:h2 , data)
> summary(reg)

Call:
lm(formula = Happiness ~ City + Gender + Employment:Bluecolor + 
    Employment:Whitecolor + Employment:h1 + Employment:h2, data = data)

Residuals:
         1          2          3          4          5          6          7          8 
 1.775e+01  1.775e+01 -1.775e+01  8.882e-16 -1.050e+01  1.050e+01  4.441e-15 -1.775e+01 

Coefficients: (1 not defined because of singularities)
                      Estimate Std. Error t value Pr(>|t|)
(Intercept)              56.75      23.56   2.408    0.138
CityNew York            -14.50      27.21  -0.533    0.647
Gender1                  -2.25      35.99  -0.063    0.956
Employment:Bluecolor     25.00      43.02   0.581    0.620
Employment:Whitecolor    57.75      35.99   1.604    0.250
Employment:h1           -50.00      54.42  -0.919    0.455
Employment:h2               NA         NA      NA       NA

Residual standard error: 27.21 on 2 degrees of freedom
Multiple R-squared:  0.6798,    Adjusted R-squared:  -0.1208 
F-statistic: 0.8491 on 5 and 2 DF,  p-value: 0.619

即使Employment:WorktypeUnemployed 不再存在，估计的系数也是相同的。

但是，Employment:h2 的 NA 值仍然存在（等效于 Employment:Holiday2 day a week）。这似乎是因为在这个简化的数据集中，您最终会得到一个奇异的模型矩阵（即一列是其他列的线性组合）

> solve(crossprod(model.matrix(reg)))
Error in solve.default(crossprod(model.matrix(reg))) : 
  system is computationally singular: reciprocal condition number = 1.79897e-18

因此，较大的数据集可能不会出现此问题。最终，您可以尝试删除模型中的任何冗余（例如，是否有任何员工每周有 0 天的假期？如果没有，那么 1 天应该是基线，并且您将在代码中添加额外的列以表示假期 > 1）。您可以使用alias() 函数来检查是哪个词出现了问题。

【讨论】：

那么在Employment:worktypeunemployment中拥有NA基本上没有问题。我有大约的数据集。 9000 个样本。在实际分析中，我估计具有更多变量的更大模型。我仍然在几个变量中得到 NA。我很困惑，因为我为每个与工作相关的变量的每个最后一个交互项都获得了 NA。
对于相同观察的每个变量都有“失业”值是否没有问题？我的真实模型中有四个具有“失业”值的变量。
R 在存在奇异性问题时丢弃最后一个交互项，即当模型矩阵的一列是其他列的线性组合时。函数alias(reg) 可用于检查哪个术语引起了麻烦。为避免该问题，您需要调整编码以减少冗余。如果您使用reg <- lm(Happiness ~ City + Gender + Employment + Whitecolor + h2, data)，上面的示例可以正常工作，其中Employment 表示一周有1天假期的蓝领。