如何在使用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

关于您的问题有一些令人困惑的事情。您同时使用了a ~ b和b ~ a，那么您到底在看什么？

• 对比度仅适用于协变量/自变量，因为它与模型矩阵的构建有关；因此a ~ b，应当将对比应用于b，而b ~ a对比应当应用于a;
• 对比度仅适用于因子/逻辑变量，不适用于数值变量。因此，除非您有b其他考虑因素，否则您将无法与之形成对比。

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

方法1：使用和的contrasts参数glmlm

我们可以通过（的相同contrasts参数）控制对比处理：glmlm

## 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，这是您对虚拟变量的期望。

方法2：直接将“ contrasts”属性设置为数据帧

我们可以使用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

现在我们可以glm不用contrasts参数就适合：

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