¿Cómo se debe calcular el Tipo II SS en un modelo mixto?

Question

Tengo un conjunto de datos (y el correspondiente modelo mixto) que obtiene valores p muy diferentes para una de las interacciones bidireccionales cuando se prueba usando el Tipo I (secuencial, teniendo cuidado de que sea el último), y el Tipo II (usando lmerTest y car paquetes).

La secuencial y car los paquetes "están de acuerdo" y el lmerTest El paquete es diferente.

¿Qué hacen estos métodos, por qué son diferentes, y cuál es "correcto"?

(Nótese que hay muchas preguntas sobre el Tipo I/II/III aquí, pero ninguna que yo vea que responda a esta pregunta específica sobre el Tipo II. Además, mientras estoy usando R, y estoy interesado en las implementaciones particulares usadas en R, la pregunta más grande es estadística, no sobre la codificación).

library(lme4)
library(lmerTest)
library(car)

d <- data.frame(
  Y = c(6, 4, 5, 7, 6, 8, 9, 0, 10, 9, 10, 8, 7, 7, 6, 5, 6, 7, 8, 
        7, 9, 10, 10, 6, 5, 8, 10, 4, 6, 7, 7, 10, 6, 10, 10, 7, 6, 3, 
        6, 7, 5, 7, 8, 9, 11, 10, 7, 8, 5, 6, 6, 7, 8, 5, 4, 8, 8, 8, 
        7, 9, 5, 4, 4, 6, 5, 7, 8, 3, 9, 8, 8, 7, 5, 7, 4, 5, 5, 4, 6, 
        8, 8, 6, 7, 8, 4, 2, 4, 4, 3, 6, 7, 8, 9, 8, 7, 6, 4, 2, 3, 3, 
        3, 0, 1, 4, 6, 4, 2, 0, 9, 9, 10, 8, 10, 9, 7, 9, 6, 5, 9, 7, 
        6, 2, 2, 3, 5, 7, 8, 9, 9, 8, 7, 8, 1, 0, 3, 2, 2, 5, 3, 5, 8, 
        6, 6, 6, 8, 6, 2, 3, 4, 5, 7, 6, 7, 8, 3, 3),
  A = rep(c(2, 5, 4, 4, 2, 3, 3, 5, 9, 3, 7, 2, 11), each=12),
  B = factor(rep(1:2, each=6)),
  C = factor(1:6),
  ID = factor(rep(1:13, each=12))
)
options(contrasts=c("contr.sum","contr.poly")) # for Type III car::Anova

Aquí está mi modelo, estoy notando las diferencias para la interacción B:C.

m1 <- lmer(Y ~ A * B * C + (1|ID), data=d)

Anova tipo I (de lmerTest ), con la interacción B:C en último lugar excepto para A:B:C, da p=0.04 para B:C. (Sin embargo, el orden no parece importar.)

anova(m1, type=1, ddf="Kenward-Roger")
#> Analysis of Variance Table of type I  with  Kenward-Roger 
#> approximation for degrees of freedom
#>        Sum Sq Mean Sq NumDF DenDF F.value    Pr(>F)    
#> A       8.794   8.794     1    11   2.908   0.11621    
#> B     136.641 136.641     1   121  45.178 6.269e-10 ***
#> C      14.128   2.826     5   121   0.934   0.46141    
#> A:B     0.381   0.381     1   121   0.126   0.72330    
#> A:C    45.874   9.175     5   121   3.034   0.01291 *  
#> B:C    34.821   6.964     5   121   2.303   0.04882 *  
#> A:B:C  21.860   4.372     5   121   1.446   0.21285    

Anova de tipo II de la car El paquete da los mismos resultados.

Anova(m1, type=2, test="F")
#> Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
#> 
#> Response: Y
#>             F Df Df.res    Pr(>F)    
#> A      2.9075  1     11   0.11621    
#> B     45.1784  1    121 6.269e-10 ***
#> C      0.9343  5    121   0.46141    
#> A:B    0.1259  1    121   0.72330    
#> A:C    3.0335  5    121   0.01291 *  
#> B:C    2.3026  5    121   0.04882 *  
#> A:B:C  1.4456  5    121   0.21285   

Anova de tipo II de lmerTest da p=0,88 para B:C, y es idéntico con el Tipo III en ambos lmerTest y car .

anova(m1, type=2, ddf="Kenward-Roger")
#> Analysis of Variance Table of type II  with  Kenward-Roger 
#> approximation for degrees of freedom
#>       Sum Sq Mean Sq NumDF DenDF F.value    Pr(>F)    
#> A      8.794   8.794     1    11  2.9075 0.1162113    
#> B     41.500  41.500     1   121 13.7214 0.0003209 ***
#> C     50.879  10.176     5   121  3.3645 0.0070072 ** 
#> A:B    0.381   0.381     1   121  0.1259 0.7233027    
#> A:C   45.874   9.175     5   121  3.0335 0.0129149 *  
#> B:C    5.174   1.035     5   121  0.3422 0.8864025    
#> A:B:C 21.860   4.372     5   121  1.4456 0.2128521   

Anova(m1, type=3, test="F")
#> Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)
#> 
#> Response: Y
#>                   F Df Df.res    Pr(>F)    
#> (Intercept) 88.2434  1     11 1.376e-06 ***
#> A            2.9075  1     11 0.1162113    
#> B           13.7214  1    121 0.0003209 ***
#> C            3.3645  5    121 0.0070072 ** 
#> A:B          0.1259  1    121 0.7233027    
#> A:C          3.0335  5    121 0.0129149 *  
#> B:C          0.3422  5    121 0.8864025    
#> A:B:C        1.4456  5    121 0.2128521    

anova(m1, type=3, ddf="Kenward-Roger")
#> Analysis of Variance Table of type III  with  Kenward-Roger 
#> approximation for degrees of freedom
#>       Sum Sq Mean Sq NumDF DenDF F.value    Pr(>F)    
#> A      8.794   8.794     1    11  2.9075 0.1162113    
#> B     41.500  41.500     1   121 13.7214 0.0003209 ***
#> C     50.879  10.176     5   121  3.3645 0.0070072 ** 
#> A:B    0.381   0.381     1   121  0.1259 0.7233027    
#> A:C   45.874   9.175     5   121  3.0335 0.0129149 *  
#> B:C    5.174   1.035     5   121  0.3422 0.8864025    
#> A:B:C 21.860   4.372     5   121  1.4456 0.2128521

Answer 1

La diferencia entre las pruebas de Tipo II de car::Anova y el anova El método de lmerTest se debe a cómo se manejan las variables explicativas continuas. El primer paso en la sección de detalles de help(Anova) describe cómo Anova se encarga de ellos:

Las designaciones "tipo II" y "tipo III" se han tomado prestadas de SAS, pero las definiciones utilizadas aquí no se corresponden precisamente con las empleado por el SAS. Las pruebas de tipo II se calculan de acuerdo con la principio de marginalidad, probando cada término después de todos los demás, excepto ignorando los parientes de orden superior del término; las llamadas pruebas de tipo III violan la marginalidad, probando cada término del modelo después de todas las otros. Esta definición de pruebas de tipo II corresponde a las pruebas producido por SAS para los modelos de análisis de la varianza, donde todos los Los predictores son factores, pero no de manera más general (es decir, cuando hay predictores cuantitativos). Tengan mucho cuidado al formular el modelo para pruebas de tipo III, o las hipótesis probadas no tendrán sentido.

así que mientras lmerTest implementa la definición de SAS del Tipo II como se describe en este documento , car::Anova hace algo diferente. La definición de SAS se describe aquí lo que significa que una prueba para un término es marginal para todos los términos en los que no está contenido. La definición de contención se da en el documento del SAS como:

Teniendo en cuenta dos efectos F1 y F2 , F1 se dice que está contenida en F2 siempre que se cumplan las dos condiciones siguientes:

• Ambos efectos implican las mismas variables continuas (si las hay).
• F2 tiene más variables de CLASE [factores en R] que F1 lo hace, y si F1 tiene variables de CLASE, todas ellas aparecen en F2

Así que en tu modelo B:C no está contenida en A:B:C porque A es continua y la hipótesis de tipo II de B:C es por lo tanto marginal para A:B:C . car::Anova no distingue entre factores y covariables en este contexto, por lo que si A fue un factor car::Anova y lmerTest- anova está de acuerdo:

d$A2 <- cut(d$A, quantile(d$A), include.lowest = TRUE)
m2 <- lmer(Y ~ A2 * B * C + (1|ID), data=d)
anova(m2, type=2, ddf="Ken") 
Type II Analysis of Variance Table with Kenward-Roger's method
        Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
A2      17.210   5.737     3     9  1.7035   0.23533    
B      136.641 136.641     1    99 40.5756 5.987e-09 ***
C       14.128   2.826     5    99  0.8391   0.52512    
A2:B    11.137   3.712     3    99  1.1024   0.35193    
A2:C    38.816   2.588    15    99  0.7684   0.70879    
B:C     34.821   6.964     5    99  2.0680   0.07575 .  
A2:B:C  50.735   3.382    15    99  1.0044   0.45709    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Anova(m2, type=2, test="F")
Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)

Response: Y
             F Df Df.res    Pr(>F)    
A2      1.7035  3      9   0.23533    
B      40.5756  1     99 5.987e-09 ***
C       0.8391  5     99   0.52512    
A2:B    1.1024  3     99   0.35193    
A2:C    0.7684 15     99   0.70879    
B:C     2.0680  5     99   0.07575 .  
A2:B:C  1.0044 15     99   0.45709    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

El lmerTest::show_tests() puede utilizarse para ver exactamente qué función lineal de los coeficientes del modelo constituye la hipótesis probada. Por ejemplo, tenemos que

show_tests(anova(m1, type=2))$`B:C`
      (Intercept) A B1 C1 C2 C3 C4 C5 A:B1 A:C1 A:C2 A:C3 A:C4 A:C5 B1:C1 B1:C2 B1:C3
B1:C1           0 0  0  0  0  0  0  0    0    0    0    0    0    0     1     0     0
B1:C2           0 0  0  0  0  0  0  0    0    0    0    0    0    0     0     1     0
B1:C3           0 0  0  0  0  0  0  0    0    0    0    0    0    0     0     0     1
B1:C4           0 0  0  0  0  0  0  0    0    0    0    0    0    0     0     0     0
B1:C5           0 0  0  0  0  0  0  0    0    0    0    0    0    0     0     0     0
      B1:C4 B1:C5 A:B1:C1 A:B1:C2 A:B1:C3 A:B1:C4 A:B1:C5
B1:C1     0     0       0       0       0       0       0
B1:C2     0     0       0       0       0       0       0
B1:C3     0     0       0       0       0       0       0
B1:C4     1     0       0       0       0       0       0
B1:C5     0     1       0       0       0       0       0

así que este contraste es marginal para todos los demás términos del modelo. Porque el car::Anova Prueba de tipo II para B:C es equivalente a la prueba de Tipo I, podemos ver que este contraste puede ser expresado como:

show_tests(anova(m1, type=1), fractions = TRUE)$`B:C`
      (Intercept) A     B1    C1    C2    C3    C4    C5    A:B1  A:C1  A:C2  A:C3 
B1:C1     0           0     0     0     0     0     0     0     0     0     0     0
B1:C2     0           0     0     0     0     0     0     0     0     0     0     0
B1:C3     0           0     0     0     0     0     0     0     0     0     0     0
B1:C4     0           0     0     0     0     0     0     0     0     0     0     0
B1:C5     0           0     0     0     0     0     0     0     0     0     0     0
      A:C4  A:C5  B1:C1 B1:C2 B1:C3 B1:C4 B1:C5 A:B1:C1 A:B1:C2 A:B1:C3 A:B1:C4 A:B1:C5
B1:C1     0     0     1   1/2   1/2   1/2   1/2 60/13   30/13   30/13   30/13   30/13  
B1:C2     0     0     0     1   1/3   1/3   1/3     0   60/13   20/13   20/13   20/13  
B1:C3     0     0     0     0     1   1/4   1/4     0       0   60/13   15/13   15/13  
B1:C4     0     0     0     0     0     1   1/5     0       0       0   60/13   12/13  
B1:C5     0     0     0     0     0     0     1     0       0       0       0   60/13  

que no es marginal para el término de interacción de tres vías.

Answer 2

La diferencia es que A no está centrada en el cero. Aún no estoy seguro de por qué esto importa, o si debería (no parece que deba...); otras respuestas que exploren esto más a fondo son muy bienvenidas.

d$Ac <- d$A - mean(d$A)
m2 <- lmer(Y ~ Ac*B*C + (1|ID), data=d)

anova(m2, type=1)
#> Type I Analysis of Variance Table with Satterthwaite's method
#>         Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
#> Ac       8.794   8.794     1    11  2.9075   0.11621    
#> B      136.641 136.641     1   121 45.1784 6.269e-10 ***
#> C       14.128   2.826     5   121  0.9343   0.46141    
#> Ac:B     0.381   0.381     1   121  0.1259   0.72330    
#> Ac:C    45.874   9.175     5   121  3.0335   0.01291 *  
#> B:C     34.821   6.964     5   121  2.3026   0.04882 *  
#> Ac:B:C  21.860   4.372     5   121  1.4456   0.21285    

anova(m2, type=2)
#> Type II Analysis of Variance Table with Satterthwaite's method
#>         Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
#> Ac       8.794   8.794     1    11  2.9075   0.11621    
#> B      136.641 136.641     1   121 45.1784 6.269e-10 ***
#> C       14.128   2.826     5   121  0.9343   0.46141    
#> Ac:B     0.381   0.381     1   121  0.1259   0.72330    
#> Ac:C    45.874   9.175     5   121  3.0335   0.01291 *  
#> B:C     34.821   6.964     5   121  2.3026   0.04882 *  
#> Ac:B:C  21.860   4.372     5   121  1.4456   0.21285    

anova(m2, type=3)
#> Type III Analysis of Variance Table with Satterthwaite's method
#>         Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
#> Ac       8.794   8.794     1    11  2.9075   0.11621    
#> B      136.641 136.641     1   121 45.1784 6.269e-10 ***
#> C       14.128   2.826     5   121  0.9343   0.46141    
#> Ac:B     0.381   0.381     1   121  0.1259   0.72330    
#> Ac:C    45.874   9.175     5   121  3.0335   0.01291 *  
#> B:C     34.821   6.964     5   121  2.3026   0.04882 *  
#> Ac:B:C  21.860   4.372     5   121  1.4456   0.21285