的R - 如何对比代码的因素,并保留在输出总之有意义的标签
问题描述:
好大家好,一劳永逸,你怎么了(重视你,因为我敢肯定有实现这一目标的方法不止一种)的对比度代码(治疗,总和,helmert等),并保留一个有意义的因子标签(所以你可以对glm函数做出有意义的解释)?的R - 如何对比代码的因素,并保留在输出总之有意义的标签
我明白我可以使用水平(),以了解哪些因素电平为参考,但得到的单调乏味,当我开始用的因素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)解决方案。
答
嗯,简单的答案(用于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
感谢的人!当你说你“决定用‘AC’和‘BC’作为标签,你在哪里调用的代码?我习惯使用relevel并重新运行我的对比参考电平进行分类。 – gh0strider18
我应该有在'paste(names(x [x> 0]),names(x [x MrFlick