¿Son correctos los grados de libertad en lmerTest::anova? Son muy diferentes a los de RM-ANOVA

Question

Estoy analizando los resultados de un experimento de tiempo de reacción en R.

He realizado un ANOVA de medidas repetidas (1 factor intra-sujeto con 2 niveles y 1 factor inter-sujeto con 2 niveles). Corrí un modelo lineal mixto similar y quería resumir los resultados del lmer en forma de tabla de ANOVA usando lmerTest::anova .

No me malinterpreten: no esperaba resultados idénticos, sin embargo no estoy seguro de los grados de libertad en lmerTest::anova resultados. Me parece que refleja más bien un ANOVA sin agregación a nivel de sujeto.

Soy consciente de que el cálculo de los grados de libertad en los modelos de efectos mixtos es complicado, pero lmerTest::anova se menciona como una posible solución en la actualización ?pvalues tema ( lme4 paquete).

¿Es correcto este cálculo? ¿Los resultados de lmerTest::anova reflejan correctamente el modelo especificado?

Actualización: Hice las diferencias individuales más grandes. Los grados de libertad en lmerTest::anova son más diferentes del anova simple, pero todavía no estoy seguro de por qué son tan grandes para el factor/interacción dentro del sujeto.

# 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 del código anterior [ actualizado ]:

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

Answer 1

I piense en que lmerTest es acertar y ezanova se equivoca en este caso.

• los resultados de lmerTest coinciden con mi intuición/comprensión
• dos cálculos diferentes en lmerTest (Satterthwaite y Kenward-Roger) están de acuerdo
• también están de acuerdo con nlme::lme
• cuando lo ejecute, ezanova hace una advertencia, que no comprendo del todo, pero que no debe ser ignorada...

Repetición del ejemplo:

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

Determinar el diseño experimental

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

Así que parece que los individuos ( subnum ) se colocan en grupos de control o de tratamiento, y cada uno de ellos se somete a pruebas en ambas direcciones, es decir, la dirección puede probarse dentro de los individuos (el denominador df es grande), pero el grupo y el grupo:la dirección sólo pueden probarse entre los individuos

(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

Aquí tengo 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. El denominador DF parece un poco raro (todos iguales a 18): creo que deberían ser más grandes para la dirección y el grupo:dirección, que se pueden probar independientemente (pero serían más pequeños si se añadiera (direction|subnum) al 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

el Df se refiere al numerador df, Denom (penúltimo) da el denominador estimado df; coinciden con la intuición clásica. Y lo que es más importante, también obtenemos respuestas diferentes para los valores de F...

También podemos comprobarlo con Kenward-Roger ( muy lento porque implica reajustar el modelo varias veces)

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

Los resultados son idénticos.

Para este ejemplo lme (de la nlme ) en realidad hace un trabajo perfectamente bueno adivinando el denominador apropiado df (los valores F y p son muy ligeramente 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

---

Si encajo una interacción entre direction y subnum el df de direction y group:direction son mucho más pequeños (yo habría pensado que serían de 18, pero quizá me esté equivocando):

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

Answer 2

En general, estoy de acuerdo con el análisis de Ben, pero permítanme añadir un par de observaciones y una pequeña intuición.

En primer lugar, los resultados generales:

1. Los resultados de lmerTest con el método Satterthwaite son correctos
2. El método Kenward-Roger también es correcto y coincide con Satterthwaite

Ben describe el diseño en el que subnum está anidado en group mientras que direction y group:direction se cruzan con subnum . Esto significa que el término de error natural término de error (es decir, el llamado "estrato de error envolvente") para group es subnum mientras que el estrato de error adjunto para el otros términos (incluyendo subnum ) son los residuos.

Esta estructura puede representarse en el llamado diagrama de estructura factorial:

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)

Aquí los términos aleatorios van entre paréntesis, 0 representa la media global (o intercepción), [I] representa el término de error, los números del superguión son el número de niveles y los números del subguión son el número de grados de libertad asumiendo un diseño equilibrado. El diagrama indica que el término de error natural término de error (estrato de error adjunto) para group es subnum y que el numerador df para subnum que es igual al denominador df para group , es 18: 20 menos 1 df para group y 1 df para la media global. En el capítulo 2 se ofrece una introducción más completa a los diagramas de estructura factorial: https://02429.compute.dtu.dk/eBook .

Si los datos estuvieran exactamente equilibrados podríamos construir las pruebas F de una descomposición SSQ como la proporcionada por anova.lm . Dado que el conjunto de datos está muy equilibrado, podemos obtener pruebas F aproximadas de la siguiente manera:

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

Aquí todos F y p se calculan asumiendo que todos los términos tienen los residuos como su estrato de error de inclusión, y eso es cierto para todos menos para "grupo". El 'equilibrado-correcto' F -prueba para el grupo es en cambio:

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 

donde utilizamos el subnum MS en lugar del Residuals MS en el F -valor denominador.

Obsérvese que estos valores coinciden bastante bien con los 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

Las diferencias restantes se deben a que los datos no están exactamente equilibrados.

El PO compara anova.lm con anova.lmerModLmerTest , lo cual está bien, pero para comparar lo mismo con lo mismo tenemos que utilizar los mismos contrastes. En este caso hay una diferencia entre anova.lm y anova.lmerModLmerTest ya que producen pruebas de tipo I y III por defecto, respectivamente, y para este conjunto de datos hay una (pequeña) diferencia entre los contrastes de tipo I y 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

Si el conjunto de datos hubiera estado completamente equilibrado, los contrastes de tipo I habrían los contrastes de tipo III (que no se ven afectados por el número observado de número de muestras).

Una última observación es que la "lentitud" del método Kenward-Roger no se debe a de la matriz de varianza-covarianza de las observaciones. matriz de varianza-covarianza de las observaciones/residuos (5191x5191 en este caso) lo que no ocurre con el método de Satterthwaite.

Con respecto al modelo2

En cuanto al modelo2 la situación se vuelve más compleja y creo que es más fácil iniciar la discusión con otro modelo en el que he incluido la interacción "clásica" entre subnum y 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

Dado que la varianza asociada a la interacción es esencialmente cero (en la presencia de la subnum efecto principal aleatorio) la interacción no tiene ningún efecto en el cálculo de los grados de libertad del denominador, F -valores y p -valores:

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

Sin embargo, subnum:direction es el estrato de error que encierra a subnum por lo que si eliminamos subnum todo el SSQ asociado vuelve a caer en subnum:direction

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

Ahora el término de error natural para group , direction y group:direction es subnum:direction y con nlevels(with(ANT.2, subnum:direction)) = 40 y cuatro parámetros, los grados de libertad del denominador para esos términos deberían ser alrededor 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

Estos F -También se pueden aproximar las pruebas con el método "equilibrado-correcto". F -prueba:

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

Ahora pasamos al modelo 2:

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

Este modelo describe una estructura de covarianza de efectos aleatorios bastante complicada con una matriz de varianza-covarianza de 2x2. La parametrización por defecto no es fácil de tratar y es mejor que hagamos una re-parametrización del modelo:

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

Si comparamos model2 a model4 tienen igual número de efectos aleatorios; 2 para cada subnum es decir, 2*20=40 en total. Mientras que model4 estipula una única varianza para los 40 efectos aleatorios, model2 estipula que cada subnum -El par de efectos aleatorios tiene una distribución normal bivariada con una matriz de varianza-covarianza de 2x2 cuyos parámetros vienen dados por

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

Esto indica un exceso de ajuste, pero dejemos eso para otro día. Lo importante aquí es que model4 es un caso especial de model2 y que model es también un caso especial de model2 . En términos generales (e intuitivos) (direction | subnum) contiene o capta la variación asociada al efecto principal subnum así como la interacción direction:subnum . En cuanto a los efectos aleatorios, podemos pensar que estos dos efectos o estructuras capturan la variación entre filas y filas por columnas, 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

En este caso, estas estimaciones de efectos aleatorios, así como el parámetro de varianza indican que en realidad sólo tenemos un efecto principal aleatorio de subnum (variación entre filas) presente aquí. Todo esto nos lleva a que los grados de libertad del denominador de Satterthwaite en

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

es un compromiso entre estas estructuras de efecto principal y de interacción: El grupo DenDF se mantiene en 18 (anidado en subnum por diseño) pero el direction y group:direction DenDF son compromisos entre 36 ( model4 ) y 5169 ( model ).

No creo que nada aquí indique que la aproximación de Satterthwaite (o su implementación en lmerTest ) es defectuoso.

La tabla equivalente con el método Kenward-Roger da

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

No es sorprendente que KR y Satterthwaite puedan diferir, pero a efectos prácticos la diferencia en p -valores es diminuto. Mi análisis anterior indica que el DenDF para direction y group:direction no debería ser menor que ~36 y probablemente mayor que eso dado que básicamente sólo tenemos el efecto principal aleatorio de direction presente, por lo que, en todo caso, creo que esto es una indicación de que el método KR obtiene la DenDF demasiado bajo en este caso. Pero hay que tener en cuenta que los datos no apoyan realmente la (group | direction) estructura, por lo que la comparación es un poco artificial; sería más interesante si el modelo estuviera realmente respaldado.