Смешанные модели множественных сравнений для взаимодействия между непрерывным и категориальным предиктором

11

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

m <- lmer(weight ~ Time * Diet + (1 | Chick), data=ChickWeight, REML=F)

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

summary(glht(m, linfct=mcp(Diet = "Tukey")))

и, действительно, они незначительно отличаются, но как я могу сделать аналогичный тест на Time:Dietэффект? Простое введение термина взаимодействия mcpприводит к ошибке:

summary(glht(m, linfct=mcp('Time:Diet' = "Tukey")))
Error in summary(glht(m, linfct = mcp(`Time:Diet` = "Tukey"))) : 
  error in evaluating the argument 'object' in selecting a method for function
 'summary': Error in mcp2matrix(model, linfct = linfct) : 
Variable(s) ‘Time:Diet’ have been specified in ‘linfct’ but cannot be found in ‘model’!
r mixed-model multiple-comparisons Дэн М.
источник
Имеет Time*Diet, что является просто упрощением Time + Diet + Time:Diet. Использование anova(m)или summary(m)подтверждает, что термин взаимодействия находится в модели.
Дан М.

Ответы:

8

По умолчанию lmerобрабатывает опорный уровень категориального предиктора как базовый уровень и оценивает параметры для других уровней. Таким образом , вы получаете некоторые попарные сравнения в выходе по умолчанию , и вы можете получить с помощью других , relevelчтобы определить новый опорный уровень и вновь подгонку модели. Преимущество этого состоит в том, что вы можете использовать сравнение моделей или MCMC для получения значений p, но не исправляет множественные сравнения (хотя впоследствии вы можете применить свою собственную коррекцию).

Чтобы использовать multcomp, вам нужно определить контрастную матрицу. Каждая строка в контрастной матрице представляет веса для эффектов, которые вы получаете при выводе модели по умолчанию, начиная с Intercept. Так что, если вы хотите эффект, который уже включен в базовый вывод, вы просто помещаете «1» в позицию, соответствующую этому эффекту. Поскольку оценки параметров относятся к общему эталонному уровню, вы можете получить сравнение между любыми двумя другими уровнями, установив вес одного на «-1», а другого - «1». Вот как это будет работать для Time:Dietусловий в ChickWeightпримере:

contrast.matrix <- rbind("Time:Diet1 vs. Time:Diet2" =  c(0, 0, 0, 0, 0, 1, 0, 0),
                           "Time:Diet1 vs. Time:Diet3" =  c(0, 0, 0, 0, 0, 0, 1, 0),
                           "Time:Diet1 vs. Time:Diet4" =  c(0, 0, 0, 0, 0, 0, 0, 1),
                           "Time:Diet2 vs. Time:Diet3" =  c(0, 0, 0, 0, 0, -1, 1, 0),
                           "Time:Diet2 vs. Time:Diet4" =  c(0, 0, 0, 0, 0, -1, 0, 1),
                           "Time:Diet3 vs. Time:Diet4" =  c(0, 0, 0, 0, 0, 0, -1, 1))
summary(glht(m, contrast.matrix))

Предостережение: этот подход использует нормальную аппроксимацию для получения значений p, которая является несколько антиконсервативной, а затем применяет некоторую коррекцию для множественных сравнений. В результате этот метод дает вам простой доступ к столько парных оценок параметров и стандартных ошибок, сколько вы хотите, но p-значения могут или не могут быть тем, что вы хотите.

(Спасибо Скотту Джексону из r-ling-lang-L за помощь в этом)

Дэн М.
источник