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

Ejecuté un ANOVA de medidas repetidas (1 factor dentro del sujeto con 2 niveles y 1 factor entre sujetos con 2 niveles). Ejecuté un modelo lineal mixto similar y quería resumir los resultados de lmer en forma de tabla ANOVA usando lmerTest::anova.

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

Soy consciente del hecho de que calcular los grados de libertad en los modelos de efectos mixtos es complicado, pero lmerTest::anovase menciona como una posible solución en el ?pvaluestema actualizado ( lme4paquete).

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

Actualización: hice las diferencias individuales más grandes. Los grados de libertad en lmerTest::anovason más diferentes de los simples anova, 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

Creo que eso lo lmerTestestá haciendo bien y lo ezanovaestá haciendo mal en este caso.

• los resultados de lmerTestestar de acuerdo 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 ejecuto, me ezanovada una advertencia, que no entiendo del todo, pero que no debe pasarse por alto ...

Ejemplo de repetición:

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

Descubre el diseño experimental

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

Por lo tanto, parece que los individuos ( subnum) se colocan en grupos de control o de tratamiento, y cada uno se prueba para ambas direcciones, es decir, la dirección se puede probar dentro de los individuos (el denominador df es grande), pero el grupo y el grupo: la dirección solo se puede probar entre 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í obtengo 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 se ve 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 puede probar de forma independiente (pero sería más pequeño si se agrega (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

la Dfcolumna aquí se refiere al numerador df, Denom(penúltimo) da el denominador estimado df; Están de acuerdo con la intuición clásica. Más importante, también obtenemos diferentes respuestas para los valores de F ...

También podemos verificar con Kenward-Roger ( muy lento porque implica volver a colocar el modelo varias veces)

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

Los resultados son idénticos.

Para este ejemplo lme(del nlmepaquete) 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 ajusto una interacción entre directiony subnumel df para directiony group:directionson mucho más pequeños (hubiera pensado que serían 18, pero tal vez me está yendo algo mal):

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

Generalmente estoy de acuerdo con el análisis de Ben, pero permítanme agregar un par de comentarios y un poco de intuición.

Primero, 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 está de acuerdo con Satterthwaite

Ben describe el diseño en el que subnumestá anidado en groupwhile direction y con el que group:directionse cruza subnum. Esto significa que el término de error natural (es decir, el llamado "estrato de error de cierre") groupes subnummientras que el estrato de error de cierre para los otros términos (incluidos subnum) son los residuos.

Esta estructura se puede representar en un 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 están encerrados entre paréntesis, 0representa la media general (o intercepción), [I]representa el término de error, los números de superguión son el número de niveles y los números de 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 (que incluye el estrato de error) para groupes subnumy que el numerador df para subnum, que es igual al denominador df para group, es 18:20 menos 1 df para groupy 1 df para la media general. Una introducción más completa a los diagramas de estructura de factores está disponible en el capítulo 2 aquí: https://02429.compute.dtu.dk/eBook .

Si los datos estuvieran exactamente equilibrados, podríamos construir las pruebas F a partir de una descomposición SSQ según lo dispuesto 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 los valores de F y p se calculan suponiendo que todos los términos tienen los residuos como su estrato de error que los encierra, y eso es cierto para todos menos para 'grupo'. La prueba F 'equilibrada-correcta' para el grupo es:

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 usamos la subnumMS en lugar de la ResidualsMS en el denominador de valor F.

Tenga en cuenta 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 OP se compara anova.lmcon anova.lmerModLmerTest, lo cual está bien, pero para comparar like con like tenemos que usar los mismos contrastes. En este caso, hay una diferencia entre anova.lmy anova.lmerModLmerTestdado que producen pruebas de Tipo I y III por defecto, respectivamente, y para este conjunto de datos hay una diferencia (pequeña) 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 se hubiera equilibrado completamente, los contrastes de tipo I habrían sido los mismos que los contrastes de tipo III (que no se ven afectados por el número observado de muestras).

Una última observación es que la 'lentitud' del método de Kenward-Roger no se debe al reajuste del modelo, sino a que implica cálculos con la matriz de varianza-covarianza marginal de las observaciones / residuos (5191x5191 en este caso) que no es El caso del método de Satterthwaite.

En cuanto al modelo2

En cuanto al modelo 2, la situación se vuelve más compleja y creo que es más fácil comenzar la discusión con otro modelo donde he incluido la interacción 'clásica' entre subnumy 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

Debido a que la varianza asociada con la interacción es esencialmente cero (en presencia del subnumefecto principal aleatorio), el término de interacción no tiene efecto en el cálculo de los grados de libertad del denominador, los valores F y los 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

Sin embargo, subnum:directiones el estrato de error de cierre, por subnumlo que si eliminamos subnumtodos los SSQ asociados vuelven a caersubnum:direction

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

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

Estas pruebas F también se pueden aproximar con las pruebas F 'equilibradas correctas' :

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 volviendo al modelo2:

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

Este modelo describe una estructura de covarianza de efecto aleatorio bastante complicada con una matriz de varianza-covarianza de 2x2. La parametrización predeterminada no es fácil de manejar y estamos mejor con una nueva parametrización del modelo:

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

Si comparamos model2a model4, tienen igualmente muchos efectos aleatorios; 2 para cada uno subnum, es decir, 2 * 20 = 40 en total. Si bien model4estipula un único parámetro de varianza para los 40 efectos aleatorios, model2estipula que cada subnumpar de efectos aleatorios tiene una distribución normal bivariada con una matriz de varianza-covarianza 2x2 cuyos parámetros están dados por

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

Esto indica un ajuste excesivo, pero guardemos eso para otro día. El punto importante aquí es que model4es un caso especial de model2 y que tambiénmodel es un caso especial de . Hablar libremente (e intuitivamente) contiene o captura la variación asociada con el efecto principal , así como la interacción . En términos de los efectos aleatorios, podemos pensar en estos dos efectos o estructuras como capturar la variación entre filas y filas por columnas, respectivamente: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

En este caso, estas estimaciones de efectos aleatorios, así como las estimaciones de parámetros de varianza, indican que realmente solo tenemos un efecto principal aleatorio de subnum(variación entre filas) presente aquí. A lo que todo esto lleva es al denominador Satterthwaite grados de libertad 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 interacción y efecto principal: el grupo DenDF permanece en 18 (anidado subnumpor diseño) pero directiony group:directionDenDF son compromisos entre 36 ( model4) y 5169 ( model).

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

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 todos los efectos prácticos, la diferencia en los valores p es diminuta. Mi análisis anterior indica que el DenDFfor directiony group:directionno debería ser menor que ~ 36 y probablemente mayor que eso dado que básicamente solo tenemos el efecto principal aleatorio del directionpresente, por lo que, en todo caso, creo que esto es una indicación de que el método KR es DenDFdemasiado bajo en este caso. Pero tenga en cuenta que los datos realmente no son compatibles con la (group | direction)estructura, por lo que la comparación es un poco artificial: sería más interesante si el modelo fuera realmente compatible.