wol*_*oor 6 regression r linear-regression lm categorical-data
我不想要主效应,因为它与更精细的因子固定效果共线,所以有这些是烦人的NA.
在这个例子中:
lm(y ~ x * z)
我想要x(数字)和z(因素)的相互作用,但不是主要的效果z.
R文件?formula说:
'*'运算符表示因子交叉:'a*b'被解释为'a + b + a:b
因此,通过执行以下操作之一,听起来主要影响是直截了当的:
a + a:b ## main effect on `b` is dropped
b + a:b ## main effect on `a` is dropped
a:b ## main effects on both `a` and `b` are dropped
真的吗？不不不(太简单,太天真).实际上它取决于变量类a和b.
这种行为是由于一个叫做魔术的函数model.matrix.default,它从公式构造一个设计矩阵.数值变量仅包含在列中,但因子变量自动编码为多个虚拟列.正是这种虚拟重新编码是一种魔力.通常认为我们可以启用或禁用对比来控制它,但不是真的.即使在这个最简单的例子中,我们也失去了对比度的控制.问题是model.matrix.default在进行虚拟编码时有自己的规则,并且它对指定模型公式的方式非常敏感.正是由于这个原因,当存在两个因素之间的相互作用时,我们不能放弃主效应.
从您的问题来看,x是数字,z是一个因素.您可以指定具有交互的模型,但不能指定zby的主效应
y ~ x + x:z
既然x是数字,那就相当于do
y ~ x:z
这里唯一的区别是参数化(或model.matrix.default虚拟编码如何).考虑一个小例子:
set.seed(0)
y <- rnorm(10)
x <- rnorm(10)
z <- gl(2, 5, labels = letters[1:2])
fit1 <- lm(y ~ x + x:z)
#Coefficients:
#(Intercept) x x:zb
# 0.1989 -0.1627 -0.5456
fit2 <- lm(y ~ x:z)
#Coefficients:
#(Intercept) x:za x:zb
# 0.1989 -0.1627 -0.7082
从系数的名称我们看到,在第一个规范中,z对比所以它的第一级"a"不是虚拟编码,而在第二个规范中,z没有对比,并且两个级别"a"和"b"都是虚拟编码.鉴于两个规范都以三个系数结束,它们实际上是等价的(从数学上讲,两种情况下的设计矩阵具有相同的列空间),您可以通过比较它们的拟合值来验证这一点:
all.equal(fit1$fitted, fit2$fitted)
# [1] TRUE
那么为什么z在第一种情况下进行对比呢？因为否则我们有两个虚拟列x:z,并且这两列的总和只是x,与x公式中的现有模型术语别名.事实上,在这种情况下即使你要求你不想要对比,model.matrix.default也不会服从:
model.matrix.default(y ~ x + x:z,
contrast.arg = list(z = contr.treatment(nlevels(z), contrasts = FALSE)))
# (Intercept) x x:zb
#1 1 0.7635935 0.0000000
#2 1 -0.7990092 0.0000000
#3 1 -1.1476570 0.0000000
#4 1 -0.2894616 0.0000000
#5 1 -0.2992151 0.0000000
#6 1 -0.4115108 -0.4115108
#7 1 0.2522234 0.2522234
#8 1 -0.8919211 -0.8919211
#9 1 0.4356833 0.4356833
#10 1 -1.2375384 -1.2375384
那么为什么在第二种情况下z没有对比呢？因为如果是,我们在构建交互时会失去级别"a"的效果.即使你需要对比,model.matrix.default也会忽略你:
model.matrix.default(y ~ x:z,
contrast.arg = list(z = contr.treatment(nlevels(z), contrasts = TRUE)))
# (Intercept) x:za x:zb
#1 1 0.7635935 0.0000000
#2 1 -0.7990092 0.0000000
#3 1 -1.1476570 0.0000000
#4 1 -0.2894616 0.0000000
#5 1 -0.2992151 0.0000000
#6 1 0.0000000 -0.4115108
#7 1 0.0000000 0.2522234
#8 1 0.0000000 -0.8919211
#9 1 0.0000000 0.4356833
#10 1 0.0000000 -1.2375384
哦,很棒model.matrix.default.它能够做出正确的决定!
让我重申一下:当存在交互时,没有办法放弃主效应.
我不会在这里提供额外的例子,因为我有一个在为什么我得到NA系数以及如何lm降低参考水平以进行交互.请参阅那里的"交互对比"部分.简而言之,以下所有规格都给出了相同的模型(它们具有相同的拟合值):
~ year:treatment
~ year:treatment + 0
~ year + year:treatment
~ treatment + year:treatment
~ year + treatment + year:treatment
~ year * treatment
特别是,第1规格导致NA系数.
因此,一旦RHS ~包含一个year:treatment,你永远不会要求model.matrix.default放弃主要效果.
model.matrix.default有些人认为model.matrix.default很烦人,因为它似乎在伪编码中没有一致的方式.他们认为"一致的方式"是始终放弃第一个因素水平.好吧,没问题,您可以model.matrix.default通过手动执行虚拟编码来绕过,并将生成的虚拟矩阵作为变量提供给lm,等等.
但是,你仍然需要model.matrix.default的帮助可以轻松地为做虚拟编码一个(是的,只有一个)因子变量.例如,对于z前一个示例中的变量,其完整的虚拟编码如下,您可以保留其全部或部分列以进行回归.
Z <- model.matrix.default(~ z + 0) ## no contrasts (as there is no intercept)
# za zb
#1 1 0
#2 1 0
#3 1 0
#4 1 0
#5 1 0
#6 0 1
#7 0 1
#8 0 1
#9 0 1
#10 0 1
#attr(,"assign")
#[1] 1 1
#attr(,"contrasts")
#attr(,"contrasts")$z
#[1] "contr.treatment"
回到我们简单的例子,如果我们不想为对比z中y ~ x + x:z,我们可以做
Z2 <- Z[, 1:2] ## use "[" to remove attributes of `Z`
lm(y ~ x + x:Z2)
#Coefficients:
#(Intercept) x x:Z2za x:Z2zb
# 0.1989 -0.7082 0.5456 NA
毫不奇怪,我们看到了NA(因为colSums(Z2)有别名x).如果我们想强制执行对比y ~ x:z,我们可以执行以下任一操作:
Z1 <- Z[, 1]
lm(y ~ x:Z1)
#Coefficients:
#(Intercept) x:Z1
# 0.34728 -0.06571
Z1 <- Z[, 2]
lm(y ~ x:Z1)
#Coefficients:
#(Intercept) x:Z1
# 0.2318 -0.6860
但是,我并不真的推荐这种黑客攻击.当您将模型公式传递给lm等时,model.matrix.default正试图为您提供最明智的构造.而且,实际上我们想用拟合模型进行预测.如果你做了假编码自己,你会提供时也很难newdata到predict.