R - Как контрастировать коэффициенты кода и сохранять значимые метки в итоговой сводке

Question

Хорошо, раз и навсегда, как вы (акцент на вас, потому что я уверен, что существует более чем один способ достижения этого). (обработка, сумма, helmert и т. д.) и сохранить значимую метку фактора (чтобы вы могли делать осмысленные интерпретации эффектов) в функции 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 

я решил использовать «А-С» и «В-С» в качестве меток, но вы можете изменить это в коде, если вам нравится. Затем работает

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