R-如何对比代码因子并在输出摘要中保留有意义的标签

Question

好的，各位，一劳永逸地，你们如何(强调你，因为我确信有不止一种方法可以实现这一点)对比代码(治疗，总和，帮助，等等)并在glm函数中保留一个有意义的因子标签(这样您就可以对效果进行有意义的解释)？

我知道我可以使用level()来理解哪个因子级别是引用，但当我开始涉及5个或10个级别的因子以及它们之间的相互作用时，这就很繁琐了。

下面是一个简单的两因素示例，说明我的意思

outcome <- c(1,0,0,1,1,0,0,0,1, 0, 0, 1)
firstvar <- c("A", "B", "C", "C", "B", "B", "A", "A", "C", "A", "C", "B")
secondvar <- c("D", "D", "E", "F", "F", "E", "D", "E", "F", "F", "D", "E")
df <- as.data.frame(cbind(outcome, firstvar, secondvar))

df$firstvar <- as.factor(df$firstvar)
df$secondvar <- as.factor(df$secondvar)

#not coded manually (and default appears to be dummy or treatment coding)
#gives meaningful factor labels in summary function
summary(glm(outcome ~ firstvar*secondvar, data=df, family="binomial"))

#effects coded
#does not give meaningful factor labels
contrasts(df$firstvar)=contr.sum(3)
contrasts(df$secondvar)=contr.sum(3)
summary(glm(outcome ~ firstvar*secondvar, data=df, family="binomial"))

#dummy coded
contrasts(df$firstvar)=contr.treatment(3); 
contrasts(df$secondvar)=contr.treatment(3); 
summary(glm(outcome ~ firstvar*secondvar, data=df, family="binomial"))

任何和所有的建议都将不胜感激。这个问题已经困扰了我一段时间，我相信有一个简单的解决方案。

Answer 1

嗯，简单的答案(至少对于contr.treatment )是，您应该将因子级别传递给函数，而不仅仅是总计数。在大多数情况下，这将正确设置级别名称。例如

contr.treatment(levels(df$firstvar))

#   B C
# A 0 0
# B 1 0
# C 0 1

然后R使用列名作为回归汇总中系数的标签/后缀。但是，即使在传递标签时，contr.sum也不喜欢设置列名。在这里，我们可以创建自己的包装器。

named.contr.sum<-function(x, ...) {
    if (is.factor(x)) {
        x <- levels(x)
    } else if (is.numeric(x) & length(x)==1L) {
        stop("cannot create names with integer value. Pass factor levels")
    }
    x<-contr.sum(x, ...)
    colnames(x) <- apply(x,2,function(x) 
         paste(names(x[x>0]), names(x[x<0]), sep="-")
    )
    x
}

在这里，我们基本上调用了contr.sum，只是将列名添加到结果中(加上一些错误检查)。您可以使用以下命令运行此命令

named.contr.sum(levels(df$firstvar))

#   A-C B-C
# A   1   0
# B   0   1
# C  -1  -1

我决定使用"A-C“和"B-C”作为标签，但如果您愿意，可以在代码中更改它。然后运行

contrasts(df$firstvar)=named.contr.sum(levels(df$firstvar))
contrasts(df$secondvar)=named.contr.sum(levels(df$secondvar))

summary(glm(outcome ~ firstvar*secondvar, data=df, family="binomial"))

会给你带来

呼叫：

glm(formula = outcome ~ firstvar * secondvar, family = "binomial", 
    data = df)

Coefficients:
                           Estimate Std. Error z value Pr(>|z|)
(Intercept)              -6.855e+00  5.023e+03  -0.001    0.999
firstvarA-C              -6.855e+00  6.965e+03  -0.001    0.999
firstvarB-C               6.855e+00  6.965e+03   0.001    0.999
secondvarD-F             -6.855e+00  6.965e+03  -0.001    0.999
secondvarE-F             -6.855e+00  6.965e+03  -0.001    0.999
firstvarA-C:secondvarD-F  2.057e+01  8.473e+03   0.002    0.998
firstvarB-C:secondvarD-F -1.371e+01  1.033e+04  -0.001    0.999
firstvarA-C:secondvarE-F  7.072e-10  1.033e+04   0.000    1.000
firstvarB-C:secondvarE-F  6.855e+00  8.473e+03   0.001    0.999

Answer 2

这是一个古老的话题，但我需要为稍后阅读的任何人澄清。

即使在传递标签时，contr.sum也不喜欢设置列名。

这种现象的发生有很好的理由，即和编码将n-1个级别与大均值进行比较，例如B均值、C均值，因此以其他方式标记是错误的。

有关各种编码R LIBRARY CONTRAST CODING SYSTEMS FOR CATEGORICAL VARIABLES的参考，请参阅UCLA手册。