r - 如何使用 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的coef是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 的 coef 仍然相同。到底是怎么回事？谢谢。

Answer 1

关于您的问题，有几个令人困惑的地方。您已经使用了 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 之间的翻转，这正是您对虚拟变量的期望。

<小时/>

方法二:直接给数据框设置“contrasts”属性

我们可以使用C或contrasts来设置“contrasts”属性(C仅用于设置，而contrasts code> 也可用于查看):

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