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"))

任何和所有建议将不胜感激.这个问题困扰了我一段时间,我确信有一个简单的(ish)解决方案.

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

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

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