lme4 девиантное / tratment контрастное кодирование с взаимодействиями в R - уровни отсутствуют

У меня есть модель смешанных эффектов (с lme4) с термином двустороннего взаимодействия, каждый член имеет несколько уровней (каждый по 4), и я хотел бы исследовать их эффекты относительно их общего среднего. Я представляю этот пример здесь из набора данных автомобиля и опускаю термин ошибки, так как он не нужен для этого примера:

## shorten data frame for simplicity
df=Cars93[c(1:15),]
df=Cars93[is.element(Cars93$Make,c('Acura Integra', 'Audi 90','BMW 535i','Subaru Legacy')),]
df$Make=drop.levels(df$Make)
df$Model=drop.levels(df$Model)

## define contrasts (every factor has 4 levels)
contrasts(df$Make) = contr.treatment(4)
contrasts(df$Model) = contr.treatment(4)

## model
m1 <- lm(Price ~ Model*Make,data=df)
summary(m1)

как видите, в термине взаимодействия опущены первые уровни. И я хотел бы, чтобы все 4 уровня на выходе были связаны с большим средним (часто называемым девиантным кодированием). Я просмотрел следующие источники: https://marissabarlaz.github.io/portfolio/contrastcoding/#coding-schemes и Как изменить контрастность для сравнения со средним значением всех уровней, а не с контрольным уровнем (R, lmer)?. Однако последняя ссылка не сообщает о взаимодействиях.


r lme4 contrast
person turquoisesquid    schedule 06.07.2020    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Ответ прост: то, чего вы хотите, напрямую невозможно. Вы должны использовать немного другой подход.

В модели с взаимодействиями вы хотите использовать контрасты, в которых среднее значение равно нулю, а не конкретному уровню. В противном случае эффекты более низкого порядка (т. Е. Основные эффекты) являются не основными эффектами, а простыми эффектами (оцениваются, когда уровень другого фактора находится на его эталонном уровне). Более подробно это объясняется в моей главе о смешанных моделях: http://singmann.org/download/publications/singmann_kellen-introduction-mixed-models.pdf

Чтобы получить то, что вы хотите, вы должны подобрать модель разумным образом, а затем передать ее emmeans для сравнения с перехватом (т.е. невзвешенным большим средним). Это также работает для взаимодействий, как показано ниже (поскольку ваш код не работал, я использую warpbreaks).

afex::set_sum_contrasts() ## uses contr.sum globally
library("emmeans")

## model
m1 <- lm(breaks ~ wool * tension,data=warpbreaks)
car::Anova(m1, type = 3)

coef(m1)[1]
# (Intercept) 
#    28.14815 

## both CIs include grand mean:
emmeans(m1, "wool") 
#  wool emmean   SE df lower.CL upper.CL
#  A      31.0 2.11 48     26.8     35.3
#  B      25.3 2.11 48     21.0     29.5
# 
# Results are averaged over the levels of: tension 
# Confidence level used: 0.95 

## same using test
emmeans(m1, "wool", null = coef(m1)[1], infer = TRUE) 
#   wool emmean   SE df lower.CL upper.CL null t.ratio p.value
#  A      31.0 2.11 48     26.8     35.3 28.1  1.372  0.1764 
#  B      25.3 2.11 48     21.0     29.5 28.1 -1.372  0.1764 
# 
# Results are averaged over the levels of: tension 
# Confidence level used: 0.95 

emmeans(m1, "tension", null = coef(m1)[1], infer = TRUE) 
#  tension emmean   SE df lower.CL upper.CL null t.ratio p.value
#  L         36.4 2.58 48     31.2     41.6 28.1  3.196  0.0025 
#  M         26.4 2.58 48     21.2     31.6 28.1 -0.682  0.4984 
#  H         21.7 2.58 48     16.5     26.9 28.1 -2.514  0.0154 
# 
# Results are averaged over the levels of: wool 
# Confidence level used: 0.95 

emmeans(m1, c("tension", "wool"), null = coef(m1)[1], infer = TRUE) 
#  tension wool emmean   SE df lower.CL upper.CL null t.ratio p.value
#  L       A      44.6 3.65 48     37.2     51.9 28.1  4.499  <.0001 
#  M       A      24.0 3.65 48     16.7     31.3 28.1 -1.137  0.2610 
#  H       A      24.6 3.65 48     17.2     31.9 28.1 -0.985  0.3295 
#  L       B      28.2 3.65 48     20.9     35.6 28.1  0.020  0.9839 
#  M       B      28.8 3.65 48     21.4     36.1 28.1  0.173  0.8636 
#  H       B      18.8 3.65 48     11.4     26.1 28.1 -2.570  0.0133 
# 
# Confidence level used: 0.95

Обратите внимание, что для coef() вы, вероятно, захотите использовать fixef() для lme4 моделей.

person Henrik    schedule 07.07.2020