Os graus de liberdade no lmerTest :: anova estão corretos? Eles são muito diferentes do RM-ANOVA

10

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
Jiri Lukavsky
fonte

Respostas:

13

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
Ben Bolker
fonte
Obrigado @Ben Bolker pela sua resposta. Vou pensar nos seus comentários e fazer mais algumas experiências. Entendo o ezAnovaaviso, pois você não deve executar o 2x2 anova se, de fato, seus dados forem do design 2x2x2.
Jiri Lukavsky
11
Possivelmente o aviso que vem com ezpode ser reformulado; na verdade, tem duas partes importantes: (1) que os dados estão sendo agregados e (2) itens sobre projetos parciais. O número 1 é mais pertinente à discrepância, pois explica que, para fazer uma anova tradicional de efeitos não mistos, é necessário agregar os dados em uma única observação por célula do design. Nesse caso, queremos uma observação por assunto por nível da variável "direção" (mantendo os rótulos de grupo para os assuntos). O ezANOVA calcula isso automaticamente.
Mike Lawrence
+1, mas não tenho certeza se o ezanova está errado. Eu corri summary(aov(rt ~ group*direction + Error(subnum/direction), data=ANT.2))e dá 16 (?) Dfs para groupe 18 para directione group:direction. O fato de haver ~ 125 observações por combinação de grupo / direção é praticamente irrelevante para o RM-ANOVA; veja, por exemplo, minha própria pergunta stats.stackexchange.com/questions/286280 : a direção é testada, por assim dizer, contra assuntos interação de direção.
Ameba
Ben, acompanhando meu comentário anterior: é realmente o que você quis dizer com "Eu pensaria que eles teriam 18 anos, mas talvez eu esteja entendendo algo errado"? Se sim, então estamos de acordo. Mais uma vez, 18 concorda com RM-ANOVA e discorda com as lmerTestestimativas de ~ 125 dfs.
Ameba
11
Atualize para o acima: lmerTest::anova(model2, ddf="Kenward-Roger")retorna 18.000 df para groupe 17.987df para os outros dois fatores, o que está em excelente concordância com o RM-ANOVA (conforme ezAnova). Minha conclusão é que a aproximação de Satterthwaite falha model2por algum motivo.
Ameba
10

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)

Diagrama da estrutura fatorial

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.

Rune H Christensen
fonte
+6, obrigado, muito interessante! Algumas perguntas. (1) Onde posso ler mais sobre "incluir estrato de erro"? Pesquisei esse termo no Google e essa resposta foi a única atingida. De um modo mais geral, que literatura você recomendaria para aprender sobre essas questões? (2a) Pelo que entendi, o RM-ANOVA clássico para esse design corresponde ao seu model3. No entanto, ele usa subnum:directioncomo o termo de erro para o teste direction. Enquanto aqui você pode forçar isso a acontecer apenas excluindo (1|subnum)como em model4. Por quê? (2b) Além disso, a RM-ANOVA produz df = 18 para direction, e não 36 quando você entra model4. Por quê?
Ameba
Para meus pontos (2a + 2b), consulte summary(aov(rt ~ group*direction + Error(subnum/direction), data=ANT.2)).
Ameba
11
(1) O tópico estratos de erro e quais termos estão incluídos nos estratos derivados das expressões Quadradas Médias Esperadas para um determinado modelo / design. Este é um material "padrão" de Design de Experiências (DoE), embora esses tópicos mais técnicos sejam frequentemente descartados em variantes fáceis ("aplicadas") de tais cursos. Veja, por exemplo, os capítulos 11 e 12 em users.stat.umn.edu/~gary/book/fcdae.pdf para obter uma introdução. Aprendi o assunto com o texto equivalente de DC Montgomery e os extensos materiais extras do (recente e lamentavelmente) Professor Henrik Spliid.
Rune H Christensen
11
... Para um tratamento mais completo, a Variance Components (1992 e 2006) de Searle et al é um clássico.
Rune H Christensen
Ahh, sim, eu deveria ter visto isso: se tivermos um modelo em que ambos subnume subnum:directionsejam diferentes de zero, então anova(lm(rt2 ~ group * direction + subnum + subnum:direction, data = ANT.2)) daremos 18 df para todos os três fatores e é isso que o método KR capta. Isso pode ser visto já com model3onde KR dá a 18 df baseada em design para todos os termos , mesmo quando a variância da interação é zero, enquanto Satterthwaite reconhece o termo variância desaparecendo e ajusta a df conformidade ....
Rune H Christensen