Как выбрать структуру со случайными и фиксированными эффектами в линейных смешанных моделях?

Рассмотрим следующие данные из двустороннего дизайна предметов:

df <- "http://personality-project.org/r/datasets/R.appendix4.data"
df <- read.table(df,header=T)
head(df)

Observation Subject Task Valence Recall
1           1     Jim Free     Neg      8
2           2     Jim Free     Neu      9
3           3     Jim Free     Pos      5
4           4     Jim Cued     Neg      7
5           5     Jim Cued     Neu      9
6           6     Jim Cued     Pos     10

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

# different fixed effects with random-intercept
a0 <- lmer(Recall~1 + (1|Subject), REML=F,df)
a1 <- lmer(Recall~Task + (1|Subject), REML=F,df)
a2 <- lmer(Recall~Valence + (1|Subject), REML=F,df)
a3 <- lmer(Recall~Task+Valence + (1|Subject), REML=F,df)
a4 <- lmer(Recall~Task*Valence + (1|Subject), REML=F,df)

# different fixed effects with random-intercept-random-slope
b0 <- lmer(Recall~1 + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b1 <- lmer(Recall~Task + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b2 <- lmer(Recall~Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b3 <- lmer(Recall~Task+Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)
b4 <- lmer(Recall~Task*Valence + (1|Subject) + (0+Task|Subject) + (0+Valence|Subject), REML=F,df)

# different fixed effects with random-intercept-random-slope including variance-covariance matrix
c0 <- lmer(Recall~1 + (1 + Valence + Task|Subject), REML=F,df)
c1 <- lmer(Recall~Task + (1 + Valence + Task|Subject), REML=F,df)
c2 <- lmer(Recall~Valence + (1 + Valence + Task|Subject), REML=F,df)
c3 <- lmer(Recall~Task+Valence + (1 + Valence + Task|Subject), REML=F,df)
c4 <- lmer(Recall~Task*Valence + (1 + Valence + Task|Subject), REML=F,df)

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

2. Рекомендуется ли сначала найти подходящую структуру с фиксированными эффектами, а затем подходящую структуру случайных эффектов или наоборот (я нашел ссылки на оба варианта ...)?

3. Каков рекомендуемый способ сообщения результатов? Сообщение значения p из теста логарифмического отношения правдоподобия, сравнивающего полную смешанную модель (с рассматриваемым эффектом) с уменьшенной моделью (без рассматриваемого эффекта). Или лучше использовать тест логарифмического отношения правдоподобия, чтобы найти наилучшую подходящую модель, а затем использовать lmerTest для сообщения значений p от эффектов в наилучшей подходящей модели?

Я не уверен, что есть действительно канонический ответ на это, но я дам ему шанс.

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

Это зависит от ваших целей.

• В общем, вы должны быть очень , очень осторожны в выборе модели (см., Например, этот ответ , или этот пост , или просто Google "Harrell stepwise" ...).
• Если вы заинтересованы в том, чтобы ваши значения p были значимыми (т. Е. Вы проводите проверку подтверждающих гипотез), вам не следует выбирать модель. Тем не менее : мне не очень понятно, насколько плохи процедуры выбора модели, если вы делаете выбор модели на нефокальных участках модели , например, выбираете модель по случайным эффектам, если ваш основной интерес - вывод о фиксированных эффектах.
• Нет такой вещи, как «помещение всех моделей в один тест отношения правдоподобия» - тестирование отношения правдоподобия является попарной процедурой. Если вы хотите сделать выбор модели (например) на случайных эффектах, я бы, вероятно, рекомендовал подход «все сразу» с использованием информационных критериев, как в этом примере, - который по крайней мере избегает некоторых проблем поэтапных подходов (но не выбор модели в целом).
• Барр и соавт. В журнале памяти и языка «Keep it maximal» 2013 года (doi: 10.1016 / j.jml.2012.11.001) рекомендуется использовать максимальную модель (только).
• Шраван Васишт не согласен , утверждая, что такие модели будут недостаточно мощными и, следовательно, проблематичными, если набор данных не очень большой (а отношение сигнал / шум высокое)
• Другой разумно оправданный подход состоит в том, чтобы подогнать большую, но разумную модель, а затем, если подгонка является единственной, удалить термины, пока она больше не будет
• С некоторыми оговорками (перечисленными в FAQ по GLMM ) вы можете использовать информационные критерии для сравнения не вложенных моделей с различными случайными эффектами (хотя Брайан Рипли не согласен: см. Нижнюю часть стр. 6 здесь )

Рекомендуется ли сначала найти подходящую структуру с фиксированными эффектами, а затем подходящую структуру случайных эффектов или наоборот (я нашел ссылки на оба варианта ...)?

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

Каков рекомендуемый способ сообщения результатов? Сообщение значения p из теста логарифмического отношения правдоподобия, сравнивающего полную смешанную модель (с рассматриваемым эффектом) с уменьшенной моделью (без рассматриваемого эффекта). Или лучше использовать тест логарифмического отношения правдоподобия, чтобы найти наилучшую подходящую модель, а затем использовать lmerTest для сообщения значений p от эффектов в наилучшей подходящей модели?

Это (увы) еще один сложный вопрос. Если вы сообщаете о предельных эффектах, о которых сообщалось lmerTest, вам придется беспокоиться о маргинальности (например, являются ли оценки основных эффектов Aи Bимеют ли смысл смысл при наличии в модели Aвзаимного Bвзаимодействия); Это огромная банка червей, но она несколько смягчается, если вы используете ее contrasts="sum"в соответствии с рекомендациями afex::mixed(). Сбалансированные конструкции тоже немного помогают. Если вы действительно хотите покрыть все эти трещины, думаю, я бы порекомендовал afex::mixed, который дает вам похожий результат lmerTest, но пытается решить эти проблемы.

Обновление, май 2017 : Как оказалось, большая часть того, что я здесь написал, немного ошибочна . Некоторые обновления сделаны по всему посту.

Я полностью согласен с тем, что уже сказал Бен Болкер (спасибо за привет afex::mixed()), но позвольте мне добавить несколько более общих и конкретных мыслей по этому вопросу.

Сосредоточьтесь на фиксированных и случайных эффектах и ​​на том, как сообщать о результатах

Для типа экспериментального исследования, которое представлено в примере набора данных от Джонатана Барона, вы используете важный вопрос, как правило, имеет ли управляемый фактор общий эффект. Например, мы находим общий основной эффект или взаимодействие Task? Важным моментом является то, что в этих наборах данных обычно все факторы находятся под полным экспериментальным контролем и распределяются случайным образом. Следовательно, в центре внимания обычно находятся фиксированные эффекты.
Напротив, компоненты случайных эффектов можно рассматривать как «неприятные» параметры, которые улавливают систематическую дисперсию (то есть межиндивидуальные различия в размере эффекта), которые не обязательно важны для основного вопроса. С этой точки зрения предложение об использовании структуры максимальных случайных эффектов, высказанное Barr et al. следует несколько естественно. Легко представить, что экспериментальные манипуляции не влияют на всех людей одинаково, и мы хотим контролировать это. С другой стороны, число факторов или уровней, как правило, не слишком велико, так что опасность переоснащения кажется сравнительно небольшой.

Следовательно, я бы последовал предложению Barr et al. и укажите в качестве моих основных результатов максимальную структуру случайных эффектов и опубликуйте тесты с фиксированными эффектами. Чтобы проверить фиксированные эффекты, я бы также предложил использовать в afex::mixed()качестве отчета о тестах эффектов или факторов (вместо тестирования параметров) и вычислять эти тесты несколько разумным способом (например, использовать одну и ту же структуру случайных эффектов для всех моделей, в которых удаляется один эффект, используются контрасты от суммы до нуля, предлагаются различные методы для вычисления p- значений, ...).

Как насчет данных примера

Проблема с данными примера, которые вы привели, состоит в том, что для этого набора данных структура максимальных случайных эффектов приводит к перенасыщенной модели, поскольку в каждой ячейке проекта имеется только одна точка данных:

> with(df, table(Valence, Subject, Task))
, , Task = Cued

       Subject
Valence Faye Jason Jim Ron Victor
    Neg    1     1   1   1      1
    Neu    1     1   1   1      1
    Pos    1     1   1   1      1

, , Task = Free

       Subject
Valence Faye Jason Jim Ron Victor
    Neg    1     1   1   1      1
    Neu    1     1   1   1      1
    Pos    1     1   1   1      1

Следовательно, lmerдроссели на максимальной структуре случайных эффектов:

> lmer(Recall~Task*Valence + (Valence*Task|Subject), df)
Error: number of observations (=30) <= number of random effects (=30) for term
(Valence * Task | Subject); the random-effects parameters and the residual variance
(or scale parameter) are probably unidentifiable

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

1. Первое решение может состоять в том, чтобы удалить наибольший случайный наклон и протестировать эффекты для этой модели:

   require(afex)
   mixed(Recall~Task*Valence + (Valence+Task|Subject), df)
           Effect    F ndf  ddf F.scaling p.value
   1         Task 6.56   1 4.00      1.00     .06
   2      Valence 0.80   2 3.00      0.75     .53
   3 Task:Valence 0.42   2 8.00      1.00     .67

   Тем не менее, это решение нерегулярно и не слишком мотивировано.

   Обновление май 2017: этот подход я сейчас поддерживаю. См. Этот пост в блоге и черновик главы, в которой я являюсь соавтором , раздел «Структуры случайных эффектов для традиционных конструкций ANOVA».

2. Альтернативное решение (и то, которое можно было бы рассматривать в качестве аргумента в обсуждении Барра и др.) Может заключаться в том, чтобы всегда удалять случайные наклоны для наименьшего эффекта. Это имеет две проблемы: (1) Какую структуру случайных эффектов мы используем, чтобы выяснить, каков наименьший эффект, и (2) R не хочет удалять эффект более низкого порядка, такой как основной эффект, если эффекты более высокого порядка, такие как взаимодействие этого эффекта присутствует (см. здесь ). Как следствие, необходимо вручную настроить эту структуру случайных эффектов и передать построенную таким образом матрицу модели вызову lmer.

3. Третьим решением может быть использование альтернативной параметризации части случайных эффектов, а именно той, которая соответствует модели RM-ANOVA для этих данных. К сожалению (?), lmerНе допускает «отрицательных отклонений», поэтому эта параметризация не совсем соответствует RM-ANOVA для всех наборов данных , см. Обсуждение здесь и в других местах (например, здесь и здесь ). «Lmer-ANOVA» для этих данных будет:

   > mixed(Recall~Task*Valence + (1|Subject) + (1|Task:Subject) + (1|Valence:Subject), df)
           Effect    F ndf  ddf F.scaling p.value
   1         Task 7.35   1 4.00      1.00     .05
   2      Valence 1.46   2 8.00      1.00     .29
   3 Task:Valence 0.29   2 8.00      1.00     .76

Учитывая все эти проблемы, я просто не стал бы использовать lmerдля подбора наборов данных, для которых имеется только одна точка данных на ячейку проекта, если не доступно более согласованное решение проблемы структуры максимальных случайных эффектов.

4. Вместо этого я бы также мог использовать классический ANOVA. Использование одной из оболочек car::Anova()в afexрезультатах будет:

   > aov4(Recall~Task*Valence + (Valence*Task|Subject), df)
           Effect         df  MSE      F  ges   p
   1      Valence 1.44, 5.75 4.67   1.46  .02 .29
   2         Task       1, 4 4.08 7.35 +  .07 .05
   3 Valence:Task 1.63, 6.52 2.96   0.29 .003 .71

   Обратите внимание, что afexтеперь также разрешается возвращать модель, с aovкоторой она может быть передана lsmeansдля специальных испытаний (но для проверки воздействий те, о которых сообщается, car::Anovaвсе же более разумны):

   > require(lsmeans)
   > m <- aov4(Recall~Task*Valence + (Valence*Task|Subject), df, return = "aov")
   > lsmeans(m, ~Task+Valence)
    Task Valence lsmean       SE   df lower.CL upper.CL
    Cued Neg       11.8 1.852026 5.52  7.17157 16.42843
    Free Neg       10.2 1.852026 5.52  5.57157 14.82843
    Cued Neu       13.0 1.852026 5.52  8.37157 17.62843
    Free Neu       11.2 1.852026 5.52  6.57157 15.82843
    Cued Pos       13.6 1.852026 5.52  8.97157 18.22843
    Free Pos       11.0 1.852026 5.52  6.37157 15.62843
   
   Confidence level used: 0.95