Equivalencia entre un modelo anova de medidas repetidas y un modelo mixto: lmer vs lme, y simetría compuesta

Tengo algunos problemas para obtener resultados equivalentes entre un aovmodelo de medidas repetidas entre dentro y un lmermodelo mixto.

Mis datos y script se ven de la siguiente manera

data=read.csv("https://www.dropbox.com/s/zgle45tpyv5t781/fitness.csv?dl=1")
data$id=factor(data$id)
data
   id  FITNESS      TEST PULSE
1   1  pilates   CYCLING    91
2   2  pilates   CYCLING    82
3   3  pilates   CYCLING    65
4   4  pilates   CYCLING    90
5   5  pilates   CYCLING    79
6   6  pilates   CYCLING    84
7   7 aerobics   CYCLING    84
8   8 aerobics   CYCLING    77
9   9 aerobics   CYCLING    71
10 10 aerobics   CYCLING    91
11 11 aerobics   CYCLING    72
12 12 aerobics   CYCLING    93
13 13    zumba   CYCLING    63
14 14    zumba   CYCLING    87
15 15    zumba   CYCLING    67
16 16    zumba   CYCLING    98
17 17    zumba   CYCLING    63
18 18    zumba   CYCLING    72
19  1  pilates   JOGGING   136
20  2  pilates   JOGGING   119
21  3  pilates   JOGGING   126
22  4  pilates   JOGGING   108
23  5  pilates   JOGGING   122
24  6  pilates   JOGGING   101
25  7 aerobics   JOGGING   116
26  8 aerobics   JOGGING   142
27  9 aerobics   JOGGING   137
28 10 aerobics   JOGGING   134
29 11 aerobics   JOGGING   131
30 12 aerobics   JOGGING   120
31 13    zumba   JOGGING    99
32 14    zumba   JOGGING    99
33 15    zumba   JOGGING    98
34 16    zumba   JOGGING    99
35 17    zumba   JOGGING    87
36 18    zumba   JOGGING    89
37  1  pilates SPRINTING   179
38  2  pilates SPRINTING   195
39  3  pilates SPRINTING   188
40  4  pilates SPRINTING   189
41  5  pilates SPRINTING   173
42  6  pilates SPRINTING   193
43  7 aerobics SPRINTING   184
44  8 aerobics SPRINTING   179
45  9 aerobics SPRINTING   179
46 10 aerobics SPRINTING   174
47 11 aerobics SPRINTING   164
48 12 aerobics SPRINTING   182
49 13    zumba SPRINTING   111
50 14    zumba SPRINTING   103
51 15    zumba SPRINTING   113
52 16    zumba SPRINTING   118
53 17    zumba SPRINTING   127
54 18    zumba SPRINTING   113

Básicamente, 3 x 6 sujetos ( id) fueron sometidos a tres FITNESSesquemas de entrenamiento diferentes cada uno y PULSEse midieron después de llevar a cabo tres tipos diferentes de resistencia TEST.

Luego ajusté el siguiente aovmodelo:

library(afex)
library(car)
set_sum_contrasts()
fit1 = aov(PULSE ~ FITNESS*TEST + Error(id/TEST),data=data)
summary(fit1)
Error: id
          Df Sum Sq Mean Sq F value   Pr(>F)    
FITNESS    2  14194    7097   115.1 7.92e-10 ***
Residuals 15    925      62                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: id:TEST
             Df Sum Sq Mean Sq F value   Pr(>F)    
TEST          2  57459   28729   253.7  < 2e-16 ***
FITNESS:TEST  4   8200    2050    18.1 1.16e-07 ***
Residuals    30   3397     113                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

El resultado que obtengo usando

set_sum_contrasts()
fit2=aov.car(PULSE ~ FITNESS*TEST+Error(id/TEST),data=data,type=3,return="Anova")
summary(fit2)

Es idéntico a esto.

Una ejecución de modelo mixto usando nlmeda un resultado directamente equivalente, por ejemplo, usando lme:

library(lmerTest)    
lme1=lme(PULSE ~ FITNESS*TEST, random=~1|id, correlation=corCompSymm(form=~1|id),data=data)
anova(lme1)
             numDF denDF   F-value p-value
(Intercept)      1    30 12136.126  <.0001
FITNESS          2    15   115.127  <.0001
TEST             2    30   253.694  <.0001
FITNESS:TEST     4    30    18.103  <.0001


summary(lme1)
Linear mixed-effects model fit by REML
 Data: data 
       AIC      BIC    logLik
  371.5375 393.2175 -173.7688

Random effects:
 Formula: ~1 | id
        (Intercept) Residual
StdDev:    1.699959 9.651662

Correlation Structure: Compound symmetry
 Formula: ~1 | id 
 Parameter estimate(s):
       Rho 
-0.2156615 
Fixed effects: PULSE ~ FITNESS * TEST 
                                 Value Std.Error DF   t-value p-value
(Intercept)                   81.33333  4.000926 30 20.328628  0.0000
FITNESSpilates                 0.50000  5.658164 15  0.088368  0.9308
FITNESSzumba                  -6.33333  5.658164 15 -1.119327  0.2806
TESTJOGGING                   48.66667  6.143952 30  7.921069  0.0000
TESTSPRINTING                 95.66667  6.143952 30 15.570868  0.0000
FITNESSpilates:TESTJOGGING   -11.83333  8.688861 30 -1.361897  0.1834
FITNESSzumba:TESTJOGGING     -28.50000  8.688861 30 -3.280062  0.0026
FITNESSpilates:TESTSPRINTING   8.66667  8.688861 30  0.997446  0.3265
FITNESSzumba:TESTSPRINTING   -56.50000  8.688861 30 -6.502579  0.0000

O usando gls:

library(lmerTest)    
gls1=gls(PULSE ~ FITNESS*TEST, correlation=corCompSymm(form=~1|id),data=data)
anova(gls1)

Sin embargo, el resultado que obtengo usando lme4's lmeres diferente:

set_sum_contrasts()
fit3=lmer(PULSE ~ FITNESS*TEST+(1|id),data=data)
summary(fit3)
Linear mixed model fit by REML ['lmerMod']
Formula: PULSE ~ FITNESS * TEST + (1 | id)
   Data: data

REML criterion at convergence: 362.4

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept)  0.00    0.0     
 Residual             96.04    9.8     
...

Anova(fit3,test.statistic="F",type=3)
Analysis of Deviance Table (Type III Wald F tests with Kenward-Roger df)

Response: PULSE
                    F Df Df.res    Pr(>F)    
(Intercept)  7789.360  1     15 < 2.2e-16 ***
FITNESS        73.892  2     15 1.712e-08 ***
TEST          299.127  2     30 < 2.2e-16 ***
FITNESS:TEST   21.345  4     30 2.030e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

¿Alguien tiene alguna idea de lo que estoy haciendo mal con el lmermodelo? ¿O de dónde viene la diferencia? ¿Podría tener algo que ver con lmerno permitir correlaciones intraclase negativas o algo así? Sin embargo, dado que nlmees glsy lmedevuelve el resultado correcto, me pregunto cómo es esto diferente en glsy lme. ¿Es que la opción correlation=corCompSymm(form=~1|id)hace que estimen directamente la correlación intraclase, que puede ser positiva o negativa, mientras que lmerestima un componente de varianza, que no puede ser negativo (y termina siendo estimado como cero en este caso)?

Como Ben Bolker ya mencionó en los comentarios, el problema es como sospecha: el lmer()modelo se dispara porque intenta estimar un modelo de componente de varianza, con las estimaciones del componente de varianza restringidas para que no sean negativas. Lo que intentaré hacer es dar una comprensión algo intuitiva de lo que se trata de su conjunto de datos que conduce a esto, y por qué esto causa un problema para los modelos de componentes de varianza.

Aquí hay una gráfica de su conjunto de datos. Los puntos blancos son las observaciones reales y los puntos negros son los medios sujetos.

Para simplificar las cosas, pero sin cambiar el espíritu del problema, restaré los efectos fijos (es decir, los efectos FITNESSy TEST, así como la gran media) y trataré los datos residuales como un problema de efectos aleatorios unidireccionales. . Así es como se ve el nuevo conjunto de datos:

Mire detenidamente los patrones en esta trama. Piense en cómo las observaciones tomadas del mismo sujeto difieren de las observaciones tomadas de diferentes sujetos. Específicamente, observe el siguiente patrón: como una de las observaciones para un sujeto es más alta (o más baja) por encima (o debajo) de la media del sujeto, las otras observaciones de ese sujeto tienden a estar en el lado opuesto de la media del sujeto. Y cuanto más lejos esté la observación de la media del sujeto, más lejos estarán las otras observaciones de la media del sujeto en el lado opuesto. Esto indica una correlación negativa dentro de la clase. Dos observaciones tomadas del mismo sujeto en realidad tienden a ser menos similares, en promedio, que dos observaciones extraídas puramente al azar del conjunto de datos.

Otra forma de pensar sobre este patrón es en términos de las magnitudes relativas de la varianza entre sujetos y dentro del sujeto. Parece que hay una varianza sustancialmente mayor dentro del sujeto en comparación con la varianza entre sujetos. Por supuesto, esperamos que esto suceda hasta cierto punto. Después de todo, la varianza dentro del sujeto se basa en la variación en los puntos de datos individuales, mientras que la varianza entre sujetos se basa en la variación en los medios de los puntos de datos individuales (es decir, el sujeto significa), y sabemos que la varianza de una media tenderá a disminuir a medida que aumente el número de cosas promediadas. Pero en este conjunto de datos la diferencia es bastante sorprendente: hay formavariación más dentro del sujeto que entre sujetos. En realidad, esta diferencia es exactamente la razón por la cual surge una correlación negativa dentro de la clase.

Bien, entonces aquí está el problema. El modelo de componente de varianza supone que cada punto de datos es la suma de un efecto sujeto y un error: , donde es el efecto del sujeto . Entonces, pensemos en lo que sucedería si realmente hubiera 0 varianza en los efectos del sujeto; en otras palabras, si el componente de varianza entre sujetos real fuera 0. Dado un conjunto de datos real generado bajo este modelo, si tuviéramos que calcular medias de muestra para los datos observados de cada sujeto, esas medias muestrales todavía tendrían alguna varianza distinta de cero, pero reflejarían solo la varianza de error, y no cualquier varianza de sujeto "verdadera" (porque asumimos que no hay ninguna).$y_{ij}=u_j+e_{ij}$$u_j$$j$

Entonces, ¿qué tan variable esperaríamos que sea este sujeto? Bueno, básicamente cada efecto sujeto estimado es una media, y conocemos la fórmula para la varianza de una media: , donde es la cantidad de cosas que se promedian. Ahora apliquemos esta fórmula a su conjunto de datos y veamos cuánta varianza esperaríamos ver en los efectos de sujeto estimados si el componente de varianza entre sujetos verdadero fuera exactamente 0.$\text{var}(\bar{X})=\text{var}(X_i)/n$$n$

La varianza dentro del sujeto resulta ser , y cada efecto del sujeto se calcula como la media de 3 observaciones. Por lo tanto, la desviación estándar esperada en el sujeto significa, suponiendo que la verdadera varianza entre sujetos sea 0, resulta ser aproximadamente . Ahora compare esto con la desviación estándar en el tema significa que realmente observamos: ¡ ! La variación observada es sustancialmente menor que la variación esperada cuando asumimos 0 varianza entre sujetos. Para un modelo de componente de varianza, la única forma en que se podría esperar que la variación observada sea tan baja como la que realmente observamos es si la verdadera varianza entre sujetos fuera de alguna manera negativa$348$$10.8$$4.3$. Y ahí radica el problema. Los datos implican que de alguna manera hay un componente de varianza negativo, pero el software (sensiblemente) no permitirá estimaciones negativas de los componentes de varianza, ya que una varianza nunca puede ser negativa. Los otros modelos que se ajustan evitan este problema al estimar directamente la correlación intraclase en lugar de asumir un modelo de componente de varianza simple.

Si desea ver cómo podría realmente obtener la estimación del componente de varianza negativa implícita en su conjunto de datos, puede utilizar el procedimiento que ilustre (con el Rcódigo adjunto ) en esta otra respuesta reciente mía . Ese procedimiento no es totalmente trivial, pero tampoco es demasiado difícil (para un diseño equilibrado como este).