2014-07-01 14 views
6

好大家好,一勞永逸,你怎麼了(重視你,因爲我敢肯定有實現這一目標的方法不止一種)的對比度代碼(治療,總和,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)解決方案。

回答

4

嗯,簡單的答案(用於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 
+1

感謝的人!當你說你「決定用‘AC’和‘BC’作爲標籤,你在哪裏調用的代碼?我習慣使用relevel並重新運行我的對比參考電平進行分類。 – gh0strider18

+0

我應該有在'paste(names(x [x> 0]),names(x [x <0]),sep =「 - 」)'line中,我使用值爲「1」的rowname,減去值爲「-1」的rowname。粘貼將「 - 」放在這些值之間。 – MrFlick