Я анализирую результаты эксперимента времени реакции в R.

Я провел повторные измерения ANOVA (1 внутри-субъектный фактор с 2 уровнями и 1 между субъектный фактор с 2 уровнями). Я запустил аналогичную линейную смешанную модель, и я хотел обобщить результаты lmer в виде таблицы ANOVA с использованием lmerTest::anova.

Не поймите меня неправильно: я не ожидал идентичных результатов, однако я не уверен в степени свободы в lmerTest::anovaрезультатах. Мне кажется, это скорее отражает ANOVA без агрегации на предметном уровне.

Мне известно о том, что вычисление степеней свободы в моделях со смешанными эффектами является сложным, но lmerTest::anovaупоминается как одно из возможных решений в обновленной ?pvaluesтеме ( lme4пакете).

Является ли этот расчет правильным? Есть ли результаты lmerTest::anovaправильно отражают указанную модель?

Обновление: я сделал индивидуальные различия больше. Степени свободы в lmerTest::anovaбольшей степени отличаются от простой ановой, но я все еще не уверен, почему они так велики для внутрисубъектного фактора / взаимодействия.

# mini example with ANT dataset from ez package
library(ez); library(lme4); library(lmerTest)

# repeated measures ANOVA with ez package
data(ANT)
ANT.2 <- subset(ANT, !error)
# update: make individual differences larger
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]

anova.ez <- ezANOVA(data = ANT.2, dv = .(rt), wid = .(subnum), 
  within = .(direction), between = .(group))
anova.ez

# similarly with lmer and lmerTest::anova
model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
lmerTest::anova(model)

# simple ANOVA on all available data
m <- lm(rt ~ group * direction, data = ANT.2)
anova(m)

Результаты кода выше [ обновлено ]:

anova.ez

$ ANOVA

           Effect DFn DFd         F          p p<.05          ges
2           group   1  18 2.6854464 0.11862957       0.1294475137
3       direction   1  18 0.9160571 0.35119193       0.0001690471
4 group:direction   1  18 4.9169156 0.03970473     * 0.0009066868

lmerTest :: ANOVA (модель)

Analysis of Variance Table of type 3  with  Satterthwaite 
approximation for degrees of freedom
                Df Sum Sq Mean Sq F value Denom Pr(>F)
group            1  13293   13293  2.6830    18 0.1188
direction        1   1946    1946  0.3935  5169 0.5305
group:direction  1  11563   11563  2.3321  5169 0.1268

ANOVA (м)

Analysis of Variance Table

Response: rt
                  Df   Sum Sq Mean Sq  F value Pr(>F)    
group              1  1791568 1791568 242.3094 <2e-16 ***
direction          1      728     728   0.0985 0.7537    
group:direction    1    12024   12024   1.6262 0.2023    
Residuals       5187 38351225    7394                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Я думаю, что lmerTestэто правильно и ezanovaнеправильно в этом случае.

• результаты lmerTestсогласны с моей интуицией / пониманием
• два разных вычисления в lmerTest(Satterthwaite и Kenward-Roger) согласны
• они также согласны с nlme::lme
• когда я запускаю его, ezanovaвыдает предупреждение, которое я не совсем понимаю, но которое не следует игнорировать ...

Повторный пример:

library(ez); library(lmerTest); library(nlme)
data(ANT)
ANT.2 <- subset(ANT, !error)
set.seed(101)  ## for reproducibility
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]

Выяснить экспериментальный дизайн

with(ANT.2,table(subnum,group,direction))

Таким образом, похоже, что индивидуумы ( subnum) помещены в контрольную или лечебную группы, и каждая из них проверяется для обоих направлений - то есть направление может быть проверено в пределах отдельных лиц (знаменатель df велико), но группа и группа: направление могут быть проверены только среди индивидуумы

(anova.ez <- ezANOVA(data = ANT.2, dv = .(rt), wid = .(subnum), 
    within = .(direction), between = .(group)))
## $ANOVA
##            Effect DFn DFd         F          p p<.05          ges
## 2           group   1  18 2.4290721 0.13651174       0.1183150147
## 3       direction   1  18 0.9160571 0.35119193       0.0002852171
## 4 group:direction   1  18 4.9169156 0.03970473     * 0.0015289914

Здесь я получаю, что Warning: collapsing data to cell means. *IF* the requested effects are a subset of the full design, you must use the "within_full" argument, else results may be inaccurate. знаменатель DF выглядит немного странно (все равны 18): я думаю, что они должны быть больше для направления и группы: направления, которое можно проверить независимо (но будет ли меньше, если вы добавите (direction|subnum)в модель)?

# similarly with lmer and lmerTest::anova
model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
lmerTest::anova(model)
##                 Df  Sum Sq Mean Sq F value Denom Pr(>F)
## group            1 12065.7 12065.7  2.4310    18 0.1364
## direction        1  1952.2  1952.2  0.3948  5169 0.5298
## group:direction  1 11552.2 11552.2  2.3299  5169 0.1270

Dfколонка здесь относится к числитель ЦФ, Denom(второй до последнего) дает оценочную знаменатель ФР; они согласны с классической интуицией. Что более важно, мы также получаем разные ответы для значений F ...

Мы также можем перепроверить с Кенвордом-Роджером ( очень медленно, потому что это требует переоснащения модели несколько раз)

lmerTest::anova(model,ddf="Kenward-Roger")

Результаты идентичны.

Для этого примера lme(из nlmeпакета) действительно отлично работает угадывание соответствующего знаменателя df (значения F и p очень немного отличаются):

model3 <- lme(rt ~ group * direction, random=~1|subnum, data = ANT.2)
anova(model3)[-1,]
##                 numDF denDF   F-value p-value
## group               1    18 2.4334314  0.1362
## direction           1  5169 0.3937316  0.5304
## group:direction     1  5169 2.3298847  0.1270

Если я соответствую взаимодействие между directionи subnumФР для directionи group:directionгораздо меньше (я бы подумал , что было бы 18, но может быть , я получаю что - то неправильно):

model2 <- lmer(rt ~ group * direction + (direction | subnum), data = ANT.2)
lmerTest::anova(model2)
##                 Df  Sum Sq Mean Sq F value   Denom Pr(>F)
## group            1 20334.7 20334.7  2.4302  17.995 0.1364
## direction        1  1804.3  1804.3  0.3649 124.784 0.5469
## group:direction  1 10616.6 10616.6  2.1418 124.784 0.1459

Я в целом согласен с анализом Бена, но позвольте мне добавить пару замечаний и немного интуиции.

Во-первых, общие результаты:

1. Результаты lmerTest с использованием метода Satterthwaite верны
2. Метод Кенварда-Роджера также верен и согласуется с Satterthwaite

Бен обрисовывает в общих чертах проект, в который subnumвложен в groupто время как direction и group:directionпересечен subnum. Это означает, что естественный термин ошибки (т. Е. Так называемый "уровень ошибок включения") для groupявляется, в subnumто время как уровень ошибок включения для других терминов (включая subnum) является остатками.

Эта структура может быть представлена ​​в так называемой диаграмме фактор-структуры:

names <- c(expression("[I]"[5169]^{5191}),
           expression("[subnum]"[18]^{20}), expression(grp:dir[1]^{4}),
           expression(dir[1]^{2}), expression(grp[1]^{2}), expression(0[1]^{1}))
x <- c(2, 4, 4, 6, 6, 8)
y <- c(5, 7, 5, 3, 7, 5)
plot(NA, NA, xlim=c(2, 8), ylim=c(2, 8), type="n", axes=F, xlab="", ylab="")
text(x, y, names) # Add text according to ’names’ vector
# Define coordinates for start (x0, y0) and end (x1, y1) of arrows:
x0 <- c(1.8, 1.8, 4.2, 4.2, 4.2, 6, 6) + .5
y0 <- c(5, 5, 7, 5, 5, 3, 7)
x1 <- c(2.7, 2.7, 5, 5, 5, 7.2, 7.2) + .5
y1 <- c(5, 7, 7, 3, 7, 5, 5)
arrows(x0, y0, x1, y1, length=0.1)

Здесь случайные термины заключены в квадратные скобки, 0представляют общее среднее значение (или точку пересечения), [I]представляют термин ошибки, числа суперскриптов представляют собой количество уровней, а числа подскриптов представляют собой количество степеней свободы, предполагающих сбалансированный дизайн. Диаграмма показывает, что естественный член ошибки (включающий в себя уровень ошибок) для groupis subnumи что числитель df для subnum, который равен знаменателю df group, равен 18: 20 минус 1 df для groupи 1 df для общего среднего. Более полное введение в диаграммы факторной структуры доступно в главе 2 здесь: https://02429.compute.dtu.dk/eBook .

Если бы данные были точно сбалансированы, мы могли бы построить F-тесты из SSQ-декомпозиции, как указано в anova.lm. Поскольку набор данных очень близко сбалансирован, мы можем получить приблизительные F-тесты следующим образом:

ANT.2 <- subset(ANT, !error)
set.seed(101)
baseline.shift <- rnorm(length(unique(ANT.2$subnum)), 0, 50)
ANT.2$rt <- ANT.2$rt + baseline.shift[as.numeric(ANT.2$subnum)]
fm <- lm(rt ~ group * direction + subnum, data=ANT.2)
(an <- anova(fm))
Analysis of Variance Table

Response: rt
                  Df   Sum Sq Mean Sq  F value Pr(>F)    
group              1   994365  994365 200.5461 <2e-16 ***
direction          1     1568    1568   0.3163 0.5739    
subnum            18  7576606  420923  84.8927 <2e-16 ***
group:direction    1    11561   11561   2.3316 0.1268    
Residuals       5169 25629383    4958                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

F_group <- an["group", "Mean Sq"] / an["subnum", "Mean Sq"]
c(Fvalue=F_group, pvalue=pf(F_group, 1, 18, lower.tail = FALSE))
   Fvalue    pvalue 
2.3623466 0.1416875 

где мы используем subnumMS вместо ResidualsMS в знаменателе F- значения.

Обратите внимание, что эти значения довольно хорошо совпадают с результатами Satterthwaite:

model <- lmer(rt ~ group * direction + (1 | subnum), data = ANT.2)
anova(model, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group           12065.3 12065.3     1    18  2.4334 0.1362
direction        1951.8  1951.8     1  5169  0.3936 0.5304
group:direction 11552.2 11552.2     1  5169  2.3299 0.1270

Остальные различия связаны с тем, что данные не являются точно сбалансированными.

ОП сравнивается anova.lmс anova.lmerModLmerTest, что нормально, но для сравнения с тем, как мы должны использовать те же контрасты. В этом случае есть разница между anova.lmи, anova.lmerModLmerTestпоскольку они производят тесты типа I и III по умолчанию соответственно, и для этого набора данных есть (небольшая) разница между контрастами типов I и III:

show_tests(anova(model, type=1))$group
               (Intercept) groupTreatment directionright groupTreatment:directionright
groupTreatment           0              1    0.005202759                     0.5013477

show_tests(anova(model, type=3))$group # type=3 is default
               (Intercept) groupTreatment directionright groupTreatment:directionright
groupTreatment           0              1              0                           0.5

Если бы набор данных был полностью сбалансирован, контрасты типа I были бы такими же, как контрасты типа III (на которые не влияет наблюдаемое количество образцов).

Последнее замечание заключается в том, что «медлительность» метода Кенварда-Роджера не связана с повторной подгонкой модели, а связана с тем, что она включает вычисления с предельной дисперсионно-ковариационной матрицей наблюдений / невязок (5191x5191 в данном случае), которая не является случай для метода Satterthwaite.

По поводу модели2

Что касается model2, ситуация становится более сложной, и я думаю, что легче начать обсуждение с другой моделью, где я включил «классическое» взаимодействие между subnumи direction:

model3 <- lmer(rt ~ group * direction + (1 | subnum) +
                 (1 | subnum:direction), data = ANT.2)
VarCorr(model3)
 Groups           Name        Std.Dev.  
 subnum:direction (Intercept) 1.7008e-06
 subnum           (Intercept) 4.0100e+01
 Residual                     7.0415e+01

Поскольку дисперсия, связанная с взаимодействием, по существу равна нулю (при наличии subnumслучайного основного эффекта) член взаимодействия не влияет на вычисление степеней свободы знаменателя, F- значений и p- значений:

anova(model3, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group           12065.3 12065.3     1    18  2.4334 0.1362
direction        1951.8  1951.8     1  5169  0.3936 0.5304
group:direction 11552.2 11552.2     1  5169  2.3299 0.1270

Тем не менее, subnum:directionявляется ли это ошибочной стратой для ошибки, subnumпоэтому, если мы удалим subnumвсе связанные SSQ, снова перейдем вsubnum:direction

model4 <- lmer(rt ~ group * direction +
                 (1 | subnum:direction), data = ANT.2)

Теперь естественным условием ошибки для group, directionа group:directionявляется subnum:directionи с nlevels(with(ANT.2, subnum:direction))= 40 и четырьмя параметрами знаменателя степени свободы для этих терминов должно быть около 36:

anova(model4, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
group           24004.5 24004.5     1 35.994  4.8325 0.03444 *
direction          50.6    50.6     1 35.994  0.0102 0.92020  
group:direction   273.4   273.4     1 35.994  0.0551 0.81583  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Эти F- тесты также могут быть аппроксимированы «сбалансированно-корректными» F- тестами:

an4 <- anova(lm(rt ~ group*direction + subnum:direction, data=ANT.2))
an4[1:3, "F value"] <- an4[1:3, "Mean Sq"] / an4[4, "Mean Sq"]
an4[1:3, "Pr(>F)"] <- pf(an4[1:3, "F value"], 1, 36, lower.tail = FALSE)
an4
Analysis of Variance Table

Response: rt
                   Df   Sum Sq Mean Sq F value Pr(>F)    
group               1   994365  994365  4.6976 0.0369 *  
direction           1     1568    1568  0.0074 0.9319    
group:direction     1    10795   10795  0.0510 0.8226    
direction:subnum   36  7620271  211674 42.6137 <2e-16 ***
Residuals        5151 25586484    4967                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Теперь переходим к model2:

model2 <- lmer(rt ~ group * direction + (direction | subnum), data = ANT.2)

Эта модель описывает довольно сложную ковариационную структуру со случайным эффектом с дисперсионно-ковариационной матрицей 2x2. С параметризацией по умолчанию нелегко иметь дело, и нам лучше перенастроить модель:

model2 <- lmer(rt ~ group * direction + (0 + direction | subnum), data = ANT.2)

Если мы сравним model2с model4, у них одинаково много случайных эффектов; 2 для каждого subnum, то есть 2 * 20 = 40 в общей сложности. В то время как model4оговаривается один параметр дисперсии для всех 40 случайных эффектов, model2оговаривается, что каждая subnumпара случайных эффектов имеет нормальное распределение с двумя вариациями с матрицей дисперсии-ковариации 2x2, параметры которой определяются

VarCorr(model2)
 Groups   Name           Std.Dev. Corr 
 subnum   directionleft  38.880        
          directionright 41.324   1.000
 Residual                70.405        

Это указывает на переоснащение, но давайте сохраним это на другой день. Важным моментом здесь является то , что model4это особый случай, model2 и что modelэто также особый случай model2. Свободно (и интуитивно) говоря, (direction | subnum)содержит или фиксирует изменения, связанные с основным эффектом, subnum а также взаимодействие direction:subnum. С точки зрения случайных эффектов мы можем думать об этих двух эффектах или структурах как о захвате изменений между строками и строками по столбцам соответственно:

head(ranef(model2)$subnum)
  directionleft directionright
1    -25.453576     -27.053697
2     16.446105      17.479977
3    -47.828568     -50.835277
4     -1.980433      -2.104932
5      5.647213       6.002221
6     41.493591      44.102056

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

anova(model2, type=1)
Type I Analysis of Variance Table with Satterthwaite's method
                 Sum Sq Mean Sq NumDF   DenDF F value Pr(>F)
group           12059.8 12059.8     1  17.998  2.4329 0.1362
direction        1803.6  1803.6     1 125.135  0.3638 0.5475
group:direction 10616.6 10616.6     1 125.136  2.1418 0.1458

является компромиссом между этими структурами основного эффекта и взаимодействия: группа DenDF остается в 18 (вложенная в subnumдизайн), но directionи group:directionDenDF являются компромиссами между 36 ( model4) и 5169 ( model).

Я не думаю, что здесь что-то указывает, что приближение Satterthwaite (или его реализация в lmerTest ) является ошибочным.

Эквивалентная таблица с методом Кенварда-Роджера дает

anova(model2, type=1, ddf="Ken")
Type I Analysis of Variance Table with Kenward-Roger's method
                 Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)
group           12059.8 12059.8     1 18.000  2.4329 0.1362
direction        1803.2  1803.2     1 17.987  0.3638 0.5539
group:direction 10614.7 10614.7     1 17.987  2.1414 0.1606

Неудивительно, что KR и Satterthwaite могут различаться, но для всех практических целей разница в p-значениях незначительна. Мой анализ выше показывает, что DenDFfor directionи group:directionне должны быть меньше, чем ~ 36 и, вероятно, больше, чем тот, что у нас есть только случайный основной эффект directionприсутствия, так что, если что-то мне кажется, это показатель того, что метод KR становится DenDFслишком низким в этом случае. Но имейте в виду, что данные на самом деле не поддерживают (group | direction)структуру, поэтому сравнение немного искусственное - было бы более интересно, если бы модель действительно поддерживалась.