Estou analisando os resultados de um experimento de tempo de reação em R.

Executei uma ANOVA de medidas repetidas (1 fator dentro do sujeito com 2 níveis e 1 fator entre sujeitos com 2 níveis). Eu executei um modelo misto linear semelhante e queria resumir os resultados anteriores na forma de tabela ANOVA usando lmerTest::anova.

Não me interpretem mal: não esperava resultados idênticos, mas não tenho certeza sobre os graus de liberdade nos lmerTest::anovaresultados. Parece-me que reflete antes uma ANOVA sem agregação no nível de assunto.

Estou ciente do fato de que calcular graus de liberdade em modelos de efeito misto é complicado, mas lmerTest::anovaé mencionado como uma solução possível no ?pvaluestópico atualizado ( lme4pacote).

Este cálculo está correto? Os resultados lmerTest::anovarefletem corretamente o modelo especificado?

Atualização: Eu fiz as diferenças individuais maiores. Os graus de liberdade lmerTest::anovasão mais diferentes da anova simples, mas ainda não tenho certeza, por que eles são tão grandes para o fator / interação dentro do sujeito.

# 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)

Resultados do código acima [ atualizado ]:

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 (modelo)

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 (m)

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

Eu acho que lmerTestestá acertando e ezanovaerrado nesse caso.

• os resultados de lmerTestconcordar com a minha intuição / compreensão
• dois cálculos diferentes em lmerTest(Satterthwaite e Kenward-Roger) concordam
• eles também concordam com nlme::lme
• quando o executo, ezanovaemite um aviso que não entendo completamente, mas que não deve ser desconsiderado ...

Exemplo de re-execução:

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)]

Descobrir o desenho experimental

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

Parece que indivíduos ( subnum) são colocados em grupos de controle ou de tratamento, e cada um é testado em ambas as direções - ou seja, a direção pode ser testada em indivíduos (o denominador df é grande), mas em grupo e grupo: a direção só pode ser testada entre indivíduos

(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

Aqui recebo 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. O denominador DF parece um pouco descolado (todos iguais a 18): acho que devem ser maiores para direção e grupo: direção, que pode ser testada independentemente (mas seria menor se você adicionasse (direction|subnum)ao modelo)?

# 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

a Dfcoluna aqui refere-se ao numerador df, Denom(penúltimo) fornece o denominador estimado df; eles concordam com a intuição clássica. Mais importante, também temos respostas diferentes para os valores de F ...

Também podemos verificar com Kenward-Roger ( muito lento porque envolve a reinstalação do modelo várias vezes)

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

Os resultados são idênticos.

Neste exemplo lme(do nlmepacote), na verdade, faz um bom trabalho adivinhando o denominador apropriado df (os valores F e p são muito ligeiramente diferentes):

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

Se eu ajustar uma interação entre directione subnumo df para directione group:directionfor muito menor (eu pensaria que eles teriam 18 anos, mas talvez eu esteja entendendo algo errado):

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

Eu geralmente concordo com a análise de Ben, mas deixe-me acrescentar algumas observações e um pouco de intuição.

Primeiro, os resultados gerais:

1. Os resultados do lmerTest usando o método Satterthwaite estão corretos
2. O método Kenward-Roger também está correto e concorda com Satterthwaite

Ben descreve o design no qual subnumestá aninhado por groupenquanto direction e group:directioné cruzado subnum. Isso significa que o termo de erro natural (isto é, o chamado "estrato de erro anexo") groupé para subnumenquanto o estrato de erro anexo para os outros termos (incluindo subnum) são os resíduos.

Essa estrutura pode ser representada no chamado diagrama de estrutura fatorial:

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)

Aqui, os termos aleatórios são colocados entre colchetes, 0representam a média geral (ou interceptação), [I]representam o termo do erro, os números do super script são o número de níveis e os números do sub script são o número de graus de liberdade, assumindo um design equilibrado. O diagrama indica que o termo de erro natural (que inclui o estrato de erro) para groupé subnume que o numerador df para subnum, que é igual ao denominador df para group, é 18: 20 menos 1 df para groupe 1 df para a média geral. Uma introdução mais abrangente aos diagramas de estrutura fatorial está disponível no capítulo 2 aqui: https://02429.compute.dtu.dk/eBook .

Se os dados fossem exatamente balanceados, poderíamos construir os testes F a partir de uma decomposição SSQ, conforme fornecido por anova.lm. Como o conjunto de dados é muito equilibrado, podemos obter testes F aproximados da seguinte maneira:

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

Aqui todos os valores de F e p são calculados assumindo que todos os termos têm os resíduos como estrato de erro anexo, e isso é verdade para todos, exceto para 'grupo'. O teste F 'equilibrado-correto' para o grupo é:

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 

onde usamos o subnumMS em vez do ResidualsMS no denominador de valor F.

Observe que esses valores correspondem muito bem aos resultados de 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

As diferenças restantes devem-se ao fato de os dados não serem exatamente equilibrados.

O OP compara anova.lmcom anova.lmerModLmerTest, o que é bom, mas para comparar como com nós temos que usar os mesmos contrastes. Nesse caso, há uma diferença entre anova.lme anova.lmerModLmerTestuma vez que eles produzem testes tipo I e III por padrão, respectivamente, e para esse conjunto de dados há uma diferença (pequena) entre os contrastes do tipo I e 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

Se o conjunto de dados tivesse sido completamente equilibrado, os contrastes do tipo I teriam sido os mesmos dos contrastes do tipo III (que não são afetados pelo número observado de amostras).

Uma última observação é que a 'lentidão' do método Kenward-Roger não se deve ao reajuste do modelo, mas porque envolve cálculos com a matriz de variância-covariância marginal das observações / resíduos (5191x5191 neste caso) que não é o caso do método de Satterthwaite.

Em relação ao modelo2

Quanto ao modelo2, a situação se torna mais complexa e acho que é mais fácil iniciar a discussão com outro modelo em que incluí a interação 'clássica' entre subnume 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

Como a variação associada à interação é essencialmente zero (na presença do subnumefeito principal aleatório), o termo de interação não afeta o cálculo dos graus de liberdade do denominador, valores F e valores 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

No entanto, subnum:directioné o estrato de erro que o encerra, subnumportanto, se removermos subnumtodo o SSQ associado, voltaremos asubnum:direction

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

Agora, o termo de erro natural para group, directione group:directioné subnum:directione com nlevels(with(ANT.2, subnum:direction))= 40 e quatro parâmetros, os graus de liberdade do denominador para esses termos devem ser cerca de 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

Esses testes F também podem ser aproximados com os testes F 'com balanceamento correto' :

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

agora voltando para o model2:

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

Este modelo descreve uma estrutura de covariância de efeito aleatório bastante complicada com uma matriz de covariância de variância 2x2. A parametrização padrão não é fácil de lidar e é melhor refazer a parametrização do modelo:

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

Se compararmos model2a model4, eles têm igualmente muitos de efeitos aleatórios; 2 para cada um subnum, ou seja, 2 * 20 = 40 no total. Embora model4estipule um único parâmetro de variância para todos os 40 efeitos aleatórios, model2estipula que cada subnumpar de efeitos aleatórios tem uma distribuição normal bi-variável com uma matriz de variância-covariância 2x2 cujos parâmetros são dados por

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

Isso indica excesso de ajuste, mas vamos deixar isso para outro dia. O ponto importante aqui é que model4é um caso especial de model2 e que modelé também um caso especial de model2. Falar livremente (e intuitivamente) (direction | subnum)contém ou captura a variação associada ao efeito principal subnum , bem como à interação direction:subnum. Em termos dos efeitos aleatórios, podemos pensar nesses dois efeitos ou estruturas como capturando variações entre linhas e linhas por colunas, respectivamente:

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

Nesse caso, essas estimativas de efeito aleatório e as estimativas dos parâmetros de variação indicam que realmente temos apenas um efeito principal aleatório de subnum(variação entre linhas) aqui apresentado. O que tudo isso leva a é que os graus de liberdade do denominador Satterthwaite em

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

é um compromisso entre essas estruturas de efeito principal e interação: O grupo DenDF permanece em 18 (aninhado subnumpor design), mas o directione group:directionDenDF são compromissos entre 36 ( model4) e 5169 ( model).

Acho que nada aqui indica que a aproximação de Satterthwaite (ou sua implementação no lmerTest ) está com defeito.

A tabela equivalente com o método Kenward-Roger fornece

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

Não é de surpreender que KR e Satterthwaite possam diferir, mas para todos os efeitos práticos a diferença nos valores de p é pequena. Minha análise acima indica que o valor DenDFfor directione group:directionnão deve ser menor que ~ 36 e provavelmente maior que o dado que basicamente temos apenas o efeito principal aleatório do directionpresente, portanto, se alguma coisa eu acho que isso é uma indicação de que o método KR fica DenDFmuito baixo nesse caso. Mas lembre-se de que os dados não suportam realmente a (group | direction)estrutura, portanto a comparação é um pouco artificial - seria mais interessante se o modelo fosse realmente suportado.