model.matrix如何为交互术语选择级别

Question

如果交互术语中包含低阶项，则model.matrix返回的级别较少。如果双因素变量分别具有na和nb水平。在包含交互项的完整model.matrix中，

model.matrix(~ A + B + A:B)，我不应该有(na-1) + (nb-1) + (na*nb-1)吗？

在下面的示例中，a和b都有三个级别。加起来，它们有九个级别。

data(mtcars)
a <- as.factor(mtcars$gear)
b <- as.factor(mtcars$cyl)
table (a,b)

   b
a    4  6  8
  3  1  2 12
  4  8  4  0
  5  2  1  2

对于只有相互作用项的模型矩阵，它具有所有九个级别。

mod.I <- model.matrix(~ a:b)
colnames(mod.I)
[1] "(Intercept)" "a3:b4"       "a4:b4"       "a5:b4"       "a3:b6"      
[6] "a4:b6"       "a5:b6"       "a3:b8"       "a4:b8"       "a5:b8"  

然而，对于只有一个低阶项的model.matrix，它也会从其他变量中降低水平。在这种情况下，b没有b= 4的项。

mod.a <- model.matrix(~ a + a:b)
colnames(mod.a)
[1] "(Intercept)" "a4"          "a5"          "a3:b6"       "a4:b6"      
[6] "a5:b6"       "a3:b8"       "a4:b8"       "a5:b8" 

这相当于完整的model.matrix。

mod.ab <- model.matrix(~ a + b + a:b)
colnames(mod.ab)
[1] "(Intercept)" "a4"          "a5"          "b6"          "b8"         
[6] "a4:b6"       "a5:b6"       "a4:b8"       "a5:b8"

我读到它与对比有关，但是，对比不会在交互项上独立运行吗？另外，如果我想知道a4:b4的系数，我该怎么做呢？

Answer 1

您可以查看model.matrix的输出，以了解它在任何特定情况下到底做了什么，但是在任何情况下，要计算完整的系数集，可以使用dummy.coef，也可以选择使用use.na=TRUE参数。请参阅?dummy.coef

fm <- lm(mpg ~ a + a:b, mtcars)
dummy.coef(fm)

给予：

Full coefficients are 

(Intercept):       21.5                                                                
a:                    3       4       5                                                
                  0.000   5.425   6.700                                                
a:b:                3:4     4:4     5:4     3:6     4:6     5:6     3:8     4:8     5:8
                  0.000   0.000   0.000  -1.750  -7.175  -8.500  -6.450   0.000 -12.800