如何在使用 R 的回归分析中为我的变量设置对比？

Question

在编码过程中，我需要更改分配给一个因子的虚拟值。但是，以下代码不起作用。有什么建议吗？

test_mx= data.frame(a= c(T,T,T,F,F,F), b= c(1,1,1,0,0,0))
test_mx
      a b
1  TRUE 1
2  TRUE 1
3  TRUE 1
4 FALSE 0
5 FALSE 0
6 FALSE 0

model= glm(b ~ a, data= test_mx, family= "binomial")
summary(model)

model= glm(a ~ b, data= test_mx, family= "binomial")
summary(model)

在这里我会得到 b 的系数是 47。现在如果我交换虚拟值，那么它应该是 -47。然而，这种情况并非如此。

test_mx2= test_mx
contrasts(test_mx2$a)
      TRUE
FALSE    0
TRUE     1
contrasts(test_mx2$a) = c(1,0)
contrasts(test_mx2$a)
      [,1]
FALSE    1
TRUE     0
model= glm(a ~ b, data= test_mx2, family= "binomial")
summary(model)

b 的系数仍然相同。到底是怎么回事？谢谢。

Answer 1

正如哲元所指出的，在 glm 建模中，对比度仅控制分类预测变量（x 值）的虚拟值分配，而不控制分类响应（y 值）。我已将此问题报告给 R 核心团队。

要手动为预测变量分配虚拟值，另一种方法是通过向量/矩阵直接分配，例如

contrasts(test_mx$a) = c(1,0)

但是，这样做存在风险：如果稍后在代码中您尝试使用 test_mx$a 作为建模中的响应值，则虚拟值分配可能会令人困惑，因为那里的分配不会与 contrasts(test_mx$a)。

Answer 2

关于您的问题，有几个令人困惑的地方。 a ~ b 和 b ~ a 你都用过，你到底在看什么？

• 对比仅适用于协变量/自变量，因为它与模型矩阵的构建有关；所以对于a ~ b，对比应该应用到b，而对于b ~ a，对比应该应用到a；
• 对比仅适用于因子/逻辑变量，而不是数值变量。因此，除非您将b 作为一个因素，否则您无法与之形成对比。

在不改变数据类型的情况下，很明显只有模型b ~ a 是合法的，可以进一步讨论。下面，我将展示如何为a设置对比度。

---

方法一：使用glm和lm的contrasts参数

我们可以通过glm 的contrasts 参数来控制对比处理（lm 也一样）：

## dropping the first factor level (default)
coef(glm(b ~ a, data = test_mx, family = binomial(),
     contrasts = list(a = contr.treatment(n = 2, base = 1))))
#(Intercept)          a2 
#  -24.56607    49.13214 

## dropping the second factor level
coef(glm(b ~ a, data = test_mx, family = binomial(),
     contrasts = list(a = contr.treatment(n = 2, base = 2))))
#(Intercept)          a1 
#   24.56607   -49.13214 

这里，contr.treatment 正在生成一个对比矩阵：

contr.treatment(n = 2, base = 1)
#  2
#1 0
#2 1

contr.treatment(n = 2, base = 2)
#  1
#1 1
#2 0

它们被传递给glm 以有效地改变model.matrix.default 的行为。让我们比较两种情况的模型矩阵：

model.matrix.default( ~ a, test_mx, contrasts.arg =
                     list(a = contr.treatment(n = 2, base = 1)))

#  (Intercept) a2
#1           1  1
#2           1  1
#3           1  1
#4           1  0
#5           1  0
#6           1  0

model.matrix.default( ~ a, test_mx, contrasts.arg =
                     list(a = contr.treatment(n = 2, base = 2)))

#  (Intercept) a1
#1           1  0
#2           1  0
#3           1  0
#4           1  1
#5           1  1
#6           1  1

a 的第二列只是 0 和 1 之间的翻转，这是您对虚拟变量的预期。

---

方法二：直接给数据框设置“对比度”属性

我们可以使用C或contrasts来设置“对比度”属性（C只用于设置，contrasts也可以用于查看）：

test_mx2 <- test_mx
contrasts(test_mx2$a) <- contr.treatment(n = 2, base = 1)
str(test_mx2)
#'data.frame':  6 obs. of  2 variables:
# $ a: Factor w/ 2 levels "FALSE","TRUE": 2 2 2 1 1 1
#  ..- attr(*, "contrasts")= num [1:2, 1] 0 1
#  .. ..- attr(*, "dimnames")=List of 2
#  .. .. ..$ : chr  "FALSE" "TRUE"
#  .. .. ..$ : chr "2"
# $ b: num  1 1 1 0 0 0

test_mx3 <- test_mx
contrasts(test_mx3$a) <- contr.treatment(n = 2, base = 2)
str(test_mx3)
#'data.frame':  6 obs. of  2 variables:
# $ a: Factor w/ 2 levels "FALSE","TRUE": 2 2 2 1 1 1
#  ..- attr(*, "contrasts")= num [1:2, 1] 1 0
#  .. ..- attr(*, "dimnames")=List of 2
#  .. .. ..$ : chr  "FALSE" "TRUE"
#  .. .. ..$ : chr "1"
# $ b: num  1 1 1 0 0 0

现在我们可以在不使用contrasts 参数的情况下拟合glm：

coef(glm(b ~ a, data = test_mx2, family = "binomial"))
#(Intercept)          a2 
#  -24.56607    49.13214 

coef(glm(b ~ a, data = test_mx3, family = "binomial"))
#(Intercept)          a1 
#   24.56607   -49.13214 

---

方法 3：设置options("contrasts") 进行全局更改

哈哈哈，@BenBolker 还提到了另一个选项，即通过设置 R 的全局选项。对于您的具体示例，因子仅涉及两个级别，我们可以使用 ?contr.SAS。

## using R default contrasts options
#$contrasts
#        unordered           ordered 
#"contr.treatment"      "contr.poly" 

coef(glm(b ~ a, data = test_mx, family = "binomial"))
#(Intercept)       aTRUE 
#  -24.56607    49.13214 

options(contrasts = c("contr.SAS", "contr.poly"))
coef(glm(b ~ a, data = test_mx, family = "binomial"))
#(Intercept)      aFALSE 
#   24.56607   -49.13214 

但我相信 Ben 只是提到这一点来完成图片；他不会在现实中采用这种方式，因为更改全局选项不利于获得可重现的 R 代码。

另一个问题是contr.SAS 只会将最后一个因子水平视为参考。在您只有 2 个级别的特定情况下，这有效地进行了“翻转”。

---

方法 4：手动重新编码因子水平

我本来不想提这个，因为它太琐碎了，但是我已经添加了“方法3”，所以我最好也添加这个。

test_mx4 <- test_mx
test_mx4$a <- factor(test_mx4$a, levels = c("TRUE", "FALSE"))
coef(glm(b ~ a, data = test_mx4, family = "binomial"))
#(Intercept)       aTRUE 
#  -24.56607    49.13214 

test_mx5 <- test_mx
test_mx5$a <- factor(test_mx5$a, levels = c("FALSE", "TRUE"))
coef(glm(b ~ a, data = test_mx5, family = "binomial"))
#(Intercept)      aFALSE 
#   24.56607   -49.13214