Линейная модель матрицы, где каждый уровень категориального контрастируется со средним - PullRequest
0 голосов
/ 21 сентября 2018

У меня есть данные xy, где y - непрерывный ответ, а x - категориальная переменная:

set.seed(1)
df <- data.frame(y = rnorm(27), group = c(rep("A",9),rep("B",9),rep("C",9)), stringsAsFactors = F)

Я хотел бы подогнать линейную модель: y ~ group к ней, в которой каждый изуровни в df$group контрастируют со средним.

Я думал, что использование Deviation Coding делает это:

lm(y ~ group,contrasts = "contr.sum",data=df)

Но оно пропускает контрастирование группы A со средним:

> summary(lm(y ~ group,contrasts = "contr.sum",data=df))

Call:
lm(formula = y ~ group, data = df, contrasts = "contr.sum")

Residuals:
    Min      1Q  Median      3Q     Max 
-1.6445 -0.6946 -0.1304  0.6593  1.9165 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  -0.2651     0.3457  -0.767    0.451
groupB        0.2057     0.4888   0.421    0.678
groupC        0.3985     0.4888   0.815    0.423

Residual standard error: 1.037 on 24 degrees of freedom
Multiple R-squared:  0.02695,   Adjusted R-squared:  -0.05414 
F-statistic: 0.3324 on 2 and 24 DF,  p-value: 0.7205

Есть ли какая-нибудь функция, которая строит model matrix, чтобы сопоставить каждый из уровней df$group со средним значением в сводке?

Все, что я могу придумать, этовручную добавив «средний» уровень к df$group и установив его в качестве базовой линии с помощью Dummy Coding :

df <- df %>% rbind(data.frame(y = mean(df$y), group ="mean"))
df$group <- factor(df$group, levels = c("mean","A","B","C"))
summary(lm(y ~ group,contrasts = "contr.treatment",data=df))

Call:
lm(formula = y ~ group, data = df, contrasts = "contr.treatment")

Residuals:
     Min       1Q   Median       3Q      Max 
-2.30003 -0.34864  0.07575  0.56896  1.42645 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.14832    0.95210   0.156    0.878
groupA       0.03250    1.00360   0.032    0.974
groupB      -0.06300    1.00360  -0.063    0.950
groupC       0.03049    1.00360   0.030    0.976

Residual standard error: 0.9521 on 24 degrees of freedom
Multiple R-squared:  0.002457,  Adjusted R-squared:  -0.1222 
F-statistic: 0.01971 on 3 and 24 DF,  p-value: 0.9961

Аналогично, предположим, у меня есть данные с двумя категориальными переменными:

set.seed(1)
df <- data.frame(y = rnorm(18),
                 group = c(rep("A",9),rep("B",9)),
                 class = as.character(rep(c(rep(1,3),rep(2,3),rep(3,3)),2)))

и я хотел бы оценить эффект взаимодействия для каждого уровня: (т.е. class1:groupB, class2:groupB и class3:groupB для:

lm(y ~ class*group,contrasts = c("contr.sum","contr.treatment"),data=df)

Как бы я его получил

1 Ответ

0 голосов
/ 21 сентября 2018

Используйте +0 в формуле lm, чтобы пропустить перехват, тогда вы должны получить ожидаемое контрастное кодирование:

summary(lm(y ~ 0 + group, contrasts = "contr.sum", data=df))

Результат:

Call:
lm(formula = y ~ 0 + group, data = df, contrasts = "contr.sum")

Residuals:
    Min      1Q  Median      3Q     Max 
-2.3000 -0.3627  0.1487  0.5804  1.4264 

Coefficients:
       Estimate Std. Error t value Pr(>|t|)
groupA  0.18082    0.31737   0.570    0.574
groupB  0.08533    0.31737   0.269    0.790
groupC  0.17882    0.31737   0.563    0.578

Residual standard error: 0.9521 on 24 degrees of freedom
Multiple R-squared:  0.02891,   Adjusted R-squared:  -0.09248 
F-statistic: 0.2381 on 3 and 24 DF,  p-value: 0.8689

Если вы хотитечтобы сделать это для взаимодействия, вот один из способов:

lm(y ~ 0 + class:group,
    contrasts = c("contr.sum","contr.treatment"), 
    data=df)
...