如何在使用 R 的回归分析中为我的变量设置对比？
How to set contrasts for my variable in regression analysis with R?
在编码过程中，我需要更改分配给一个因子的虚拟值。但是，以下代码不起作用。有什么建议吗？
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 的系数仍然相同。到底是怎么回事？谢谢
关于您的问题，有几处令人困惑。 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可以用来查看嗯）：
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
方法三：设置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 个级别的特定情况下，这有效地执行了 "flipping".
方法 4：手动重新编码您的因子水平
我本来不想提这个的，因为它太琐碎了，但是既然我添加了"Method 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
正如 Zheyuan 所指出的，对比仅控制分类预测变量（x 值）的虚拟值分配，而不控制 glm 建模中的分类响应（y 值）。我已将此问题报告给 R 核心团队。
要为预测变量手动分配虚拟值，另一种方法可以是通过 vector/matrix 直接分配，例如
contrasts(test_mx$a) = c(1,0)
但是，这样做存在风险：如果您稍后在代码中尝试使用 test_mx$a 作为建模中的响应值，虚拟值赋值可能会造成混淆，因为那里的赋值将不匹配和
contrasts(test_mx$a)。
在编码过程中，我需要更改分配给一个因子的虚拟值。但是，以下代码不起作用。有什么建议吗？
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 的系数仍然相同。到底是怎么回事？谢谢
关于您的问题，有几处令人困惑。 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可以用来查看嗯）：
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
方法三：设置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 个级别的特定情况下，这有效地执行了 "flipping".
方法 4：手动重新编码您的因子水平
我本来不想提这个的，因为它太琐碎了，但是既然我添加了"Method 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
正如 Zheyuan 所指出的，对比仅控制分类预测变量（x 值）的虚拟值分配，而不控制 glm 建模中的分类响应（y 值）。我已将此问题报告给 R 核心团队。
要为预测变量手动分配虚拟值，另一种方法可以是通过 vector/matrix 直接分配，例如
contrasts(test_mx$a) = c(1,0)
但是，这样做存在风险：如果您稍后在代码中尝试使用 test_mx$a 作为建模中的响应值，虚拟值赋值可能会造成混淆，因为那里的赋值将不匹配和
contrasts(test_mx$a)。