Значение P для члена взаимодействия в моделях смешанных эффектов с использованием lme4

Я анализирую некоторые поведенческие данные, используя lme4in R, в основном следуя отличным учебникам Бодо Винтера , но не понимаю, правильно ли я обрабатываю взаимодействия. Хуже того, никто другой, участвующий в этом исследовании, не использует смешанные модели, поэтому я немного волнуюсь, чтобы убедиться, что все правильно.

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

• Во время написания я нашел этот вопрос , показывая, что nlmeболее прямо дают значения p для терминов взаимодействия, но я думаю, что все еще уместно задавать этот вопрос в отношении lme4.
• Livius'Ответ на этот вопрос предоставил ссылки на множество дополнительных материалов, которые я постараюсь пройти в ближайшие несколько дней, поэтому я прокомментирую любой прогресс, который принесет.

В моих данных у меня есть зависимая переменная dv, conditionманипуляция (0 = контроль, 1 = экспериментальное условие, которое должно привести к более высокому значению dv), а также предварительное условие, помеченное appropriate: испытания, закодированные 1для этого, должны показать эффект, но испытания, закодированные 0могут нет, потому что решающий фактор отсутствует.

Я также включил два случайных перехвата, для subjectи для targetотражения коррелированных dvзначений в каждом предмете и в каждой из 14 решенных задач (каждый участник решал контрольную и экспериментальную версию каждой проблемы).

library(lme4)
data = read.csv("data.csv")

null_model        = lmer(dv ~ (1 | subject) + (1 | target), data = data)
mainfx_model      = lmer(dv ~ condition + appropriate + (1 | subject) + (1 | target),
                         data = data)
interaction_model = lmer(dv ~ condition + appropriate + condition*appropriate +
                              (1 | subject) + (1 | target), data = data)
summary(interaction_model)

Вывод:

## Linear mixed model fit by REML ['lmerMod']
## ...excluded for brevity....
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 0.006594 0.0812  
##  target   (Intercept) 0.000557 0.0236  
##  Residual             0.210172 0.4584  
## Number of obs: 690, groups: subject, 38; target, 14
## 
## Fixed effects:
##                                Estimate Std. Error t value
## (Intercept)                    0.2518     0.0501    5.03
## conditioncontrol               0.0579     0.0588    0.98
## appropriate                   -0.0358     0.0595   -0.60
## conditioncontrol:appropriate  -0.1553     0.0740   -2.10
## 
## Correlation of Fixed Effects:
## ...excluded for brevity.

ANOVA затем показывает, interaction_modelчто он значительно лучше подходит, чем mainfx_model, из чего я заключаю, что присутствует значительное взаимодействие (р = 0,035).

anova(mainfx_model, interaction_model)

Вывод:

## ...excluded for brevity....
##                   Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
## mainfx_model       6 913 940   -450      901                          
## interaction_model  7 910 942   -448      896  4.44      1      0.035 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Оттуда я изолирую подмножество данных, для которых выполняется appropriateтребование (т. Е. appropriate = 1), И для него подходит нулевая модель, а модель, включающая conditionв качестве эффекта, снова сравниваю две модели с использованием ANOVA, и вот, я обнаружил, что conditionявляется значимым предиктором.

good_data = data[data$appropriate == 1, ]
good_null_model   = lmer(dv ~ (1 | subject) + (1 | target), data = good_data)
good_mainfx_model = lmer(dv ~ condition + (1 | subject) + (1 | target), data = good_data)

anova(good_null_model, good_mainfx_model)

Вывод:

## Data: good_data
## models:
## good_null_model: dv ~ (1 | subject) + (1 | target)
## good_mainfx_model: dv ~ condition + (1 | subject) + (1 | target)
##                   Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
## good_null_model    4 491 507   -241      483                          
## good_mainfx_model  5 487 507   -238      477  5.55      1      0.018 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Я не вижу слишком много, чтобы сказать здесь. Я думаю, что вы проделали хорошую работу.

Есть несколько способов, которые люди обсуждали, чтобы проверить эффекты и получить p-значения для сложных моделей смешанных эффектов. Существует хороший поверхностный обзор здесь . Лучше всего использовать вычислительно-интенсивные методы (самозагрузка или байесовские методы), но для большинства людей это более продвинутый метод. Второй лучший (и самый удобный) метод - это использовать критерий отношения правдоподобия. Это то, что anova()(технически ? Anova.merMod () ) делает. Важно использовать только критерий отношения правдоподобия моделей, которые соответствовали полному максимальному правдоподобию , а не ограниченному максимальному правдоподобию(REML). С другой стороны, для вашей окончательной модели и для интерпретации вы хотите использовать REML. Это сбивает с толку многих людей. В ваших выходных данных мы видим, что вы подгоняете свои модели к REML (это потому, что опция установлена TRUEпо умолчанию в lmer(). Это означает, что ваш тест недействителен, однако, поскольку это такая распространенная ошибка, anova.merMod()содержит refitаргумент, который По умолчанию установлено значение TRUE, и вы не изменили его. Таким образом, предвидение разработчиков пакетов спасло вас там.

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

Я сам новичок и следую советам Zuur et al .. Я использую lmeиз nlmeпакета вместо того, lme4когда мне нужно добавить иерархическую структуру ошибок в линейную модель. Мой ответ может быть далеко.

Два комментария:

(1) Я не уверен, что имеет смысл тестировать conditionв подмножестве, appropriate==1только когда . Если вы хотите получить p-значения для основных эффектов, вы можете использовать Anovaиз пакета 'car':

require(car)
Anova(M,type="III")# type III Sum of Squares. M was fitted using method=REML

Если вы хотите разрешить взаимодействие, вы можете запускать парные сравнения напрямую (?) Или делать то, что вы делали, но на обоих подмножествах (то есть также с подмножеством где appropriate==0).

(2) Возможно, вы захотите сначала выбрать структуру ошибок, а не предполагать, что (1 | subject) + (1 | target)это лучшая структура ошибок. Из того, что вы написали, я понимаю, что conditionэто внутрисубъектный фактор, а также фактор appropriateмежду субъектами или между объектами. Возможно, вы захотите добавить наклоны для факторов внутри объекта и / или внутри цели, например: добавьте dv ~ condition + appropriate + (1+condition | subject) + (1 | target)случайный наклон для фактора внутри объекта condition. Никаких наклонов не требуется для факторов между субъектами / целями.

ура