2014-07-01 19 views
6

Ok, ein für allemal, wie geht es dir (Betonung auf dich, weil ich sicher bin, dass es mehr als einen Weg gibt, dies zu erreichen) Kontrast-Code (Behandlung, Summe, Helmert usw.) und behalten Sie eine aussagekräftige Faktorbezeichnung (damit Sie sinnvolle Interpretationen von Effekten vornehmen können) in der glm-Funktion?R - Wie man Code-Faktoren kontrastiert und aussagekräftige Beschriftungen in der Ausgabezusammenfassung erhält

Ich verstehe, dass ich Level() verwenden kann, um zu verstehen, welche Faktorstufe die Referenz ist, aber das wird langweilig, wenn ich anfange, Faktoren mit 5 oder 10 Ebenen und deren Interaktionen einzubeziehen.

Hier ist ein kurzer Zwei-Faktor-Beispiel dafür, was ich meine

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

Jede und alle Vorschläge wird geschätzt. Dieses Problem hat mich eine Weile gestört, und ich bin mir sicher, dass es eine einfache (ish) Lösung gibt.

Antwort

4

Nun, die einfache Antwort (für contr.treatment mindestens), ist, dass Sie die Faktorstufen an die Funktion übergeben sollten, anstatt nur die Gesamtzahl. In den meisten Fällen werden dadurch die Ebenennamen korrekt festgelegt. Beispiel:

und R verwendet dann die Spaltennamen als Bezeichnungen/Suffixe für die Koeffizienten in der Regressionsübersicht. Auch wenn Etiketten übergeben werden, möchte contr.sum keine Spaltennamen festlegen. Hier können wir jedoch einen eigenen Wrapper erstellen.

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 
} 

Hier fordern wir grundsätzlich contr.sum aufrufen und das Hinzufügen von nur Spaltennamen zu dem Ergebnis (plus einige Fehlerprüfung). Sie können diese laufen mit

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

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

entschied ich mich „A-C“ und „B-C“ als Etiketten zu verwenden, aber man konnte, dass in den Code ändern, wenn Sie möchten. Dann

läuft
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")) 

geben Ihnen

Call:

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

Thanks man! Wenn Sie sagen, dass Sie sich entschieden haben, "AC" und "BC" als Beschriftungen zu verwenden, wo haben Sie das im Code genannt? Ich bin es gewohnt, revel zu verwenden und meine Kontraste erneut zu verwenden, um den Referenzpegel zu klassifizieren. – gh0strider18

+0

sollte ich haben war klarer, es ist in der 'einfügen (Namen (x [x> 0]), Namen (x [x <0]), sep =" - ")" Linie. Ich verwende den rowname des Wertes mit "1" minus der rowname mit dem Wert "-1". Die Paste setzt das "-" zwischen diesen Werten. – MrFlick