R: la equivalencia entre un aov entre-dentro de las medidas repetidas y un modelo lmer modelo mixto

Question

Tengo algunos problemas para obtener resultados equivalentes entre un aov entre-dentro de las medidas repetidas y un modelo lmer modelo mixto.

Mis datos y el script es el siguiente

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 tipos de FITNESS entrenamiento de los esquemas de cada uno y de sus PULSE fue medido después de la realización de tres tipos diferentes de resistencia TESTs.

Me ajustaron a continuación, el siguiente aov modelo de :

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 obtiene mediante

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

es idéntico a este.

Un modelo mixto de ejecución usando nlme da directamente un resultado equivalente, utilizando, por ejemplo, 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 el uso de gls :

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

Sin embargo, el resultado que obtiene mediante lme4's lmer es 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 alguna idea de qué estoy haciendo mal con el lmer modelo? O donde la diferencia viene de? No podía hacer nada con lmer , no permitiendo que los negativos intraclase corellations o algo por el estilo? Dado que nlme's gls y lme hacer devolver el resultado correcto, sin embargo, me pregunto cómo esto es diferente en gls y lme? Es que la opción correlation=corCompSymm(form=~1|id) hace directamente la estimación de la correlación intraclase, que pueden ser positivos o negativos, mientras que lmer las estimaciones de la varianza del componente, que puede ser negativa (y termina siendo estimado como cero en este caso)?

Answer 1

Como Ben Bolker ya se ha mencionado en los comentarios, el problema es como sospecha: La lmer() modelo obtiene tropezó porque se pretende estimar un modelo de componentes de varianza, con los componentes de varianza de las estimaciones obligado a ser no negativo. Lo que voy a intentar hacer es darle un poco de comprensión intuitiva de lo que es acerca de su conjunto de datos que lleva a esto, y por qué esto provoca un problema para los modelos de componentes de varianza.

Aquí está una parcela de su conjunto de datos. Los puntos blancos son las observaciones reales y los puntos negros son el tema de los medios.

Para hacer las cosas más simples, pero sin cambiar el espíritu de el problema, voy a restar los efectos fijos (es decir, el FITNESS y TEST efectos, así como la gran media) y de acuerdo con los datos residuales como una forma de efectos aleatorios problema. Así que aquí es lo que el nuevo conjunto de datos se parece a:

Cuesta mirar a los patrones en esta parcela. Piense acerca de cómo las observaciones tomadas a partir de la misma materia difieren de las observaciones tomadas a partir de diferentes materias. Específicamente, observe el siguiente patrón: Como una de las observaciones de un objeto es mayor (o menor) por encima (o por debajo) de la media de los sujetos, las otras observaciones de que los sujetos tienden a estar en el lado opuesto de la media de los sujetos. Y, además, que la observación es de la media de los sujetos, más que la de otras observaciones tienden a ser a partir de la media de los sujetos en el lado opuesto. Esto indica un negativo de la correlación intra clase. Dos observaciones tomadas a partir de la misma materia en realidad, tienden a ser menos similar, en promedio, de dos observaciones dibujado puramente al azar del conjunto de datos.

Otra forma de pensar acerca de este patrón es, en términos de las magnitudes relativas de la entre-sujetos y de sujeto de la varianza. Parece que no es sustancialmente mayor en materia de varianza en comparación con entre-sujetos de la varianza. Por supuesto, esperamos que esto ocurra en cierta medida. Después de todo, la de sujeto de la varianza se basa en la variación en los puntos de datos individuales, mientras que el entre-sujetos de la varianza se basa en la variación en medio de los puntos de datos individuales (es decir, el sujeto que significa), y sabemos que la varianza de la media tienden a disminuir a medida que el número de cosas que está haciendo el promedio de los aumentos. Pero en este conjunto de datos la diferencia es bastante notable: no Hay manera más dentro de los sujetos de entre-sujetos de la variación. En realidad, esta diferencia es exactamente la razón por la negativa correlación intra clase emerge.

Bueno, así que aquí está el problema. Los componentes de varianza modelo asume que cada punto de datos es la suma de un sujeto, efecto y un error: $y_{ij}=u_j+e_{ij}$ donde $u_j$ es el efecto de la $j$th tema. Así que vamos a pensar en qué pasaría si hubo verdaderamente 0 varianza en el tema de los efectos, en otras palabras, si la verdad entre-sujetos de componentes de varianza fueron 0. Dado un conjunto de datos reales generados bajo este modelo, si vamos a calcular la muestra significa para cada sujeto de los datos observados, los de la muestra significa que todavía tiene algunas no-cero, varianza, pero ellos sólo reflejan la varianza de error, y no cualquier "verdadero" sujeto de la varianza (debido a que hemos supuesto que no hay ninguno).

Entonces, ¿cómo variable sería esperamos que estos sujetos significa ser? Bueno, básicamente cada estima sujetos efecto es un decir, y sabemos que la fórmula para la varianza de una media: $\text{var}(\bar{X})=\text{var}(X_i)/n$ donde $n$ es el número de cosas que se promedien. Ahora vamos a aplicar esta fórmula a su conjunto de datos y ver qué parte de la varianza podemos esperar ver en la estimación de los sujetos en los efectos si la verdad entre-sujetos de componentes de varianza fueron exactamente 0.

La de sujeto de la varianza de las obras a ser $348$, y cada uno de los sujetos efecto se calcula como la media de las 3 observaciones. Así que a la espera de la desviación estándar en el sujeto significa -- suponiendo que la verdad entre los sujetos de la varianza es 0 -- obras a ser acerca de $10.8$. Ahora compara esto a la desviación estándar en el sujeto significa que en realidad observada: $4.3$! La variación observada es sustancialmente menor que la variación esperada cuando asumimos 0 entre los sujetos de la varianza. Para un modelo de componentes de varianza, la única manera de que la variación observada podría ser tan bajo como lo que en realidad se observa es que si la verdad entre los sujetos de la varianza de alguna manera fueron negativos. Y ahí radica el problema. Los datos indican que no es de alguna manera un negativo de componentes de varianza, pero el software (sensatez) no permitir la negativa de las estimaciones de los componentes de varianza, ya que una variación de hecho, puede nunca ser negativo. Los otros modelos que se ajustan evitar este problema directamente la estimación de la correlación intra clase en lugar de asumir un modelo de componentes de varianza simple.

Si quieres ver cómo se podría conseguir realmente la negativa de componentes de varianza de la estimación implícita por el conjunto de datos, puede utilizar el procedimiento que se ilustra (con el acompañamiento R código) en este otro reciente respuesta de la mina. Ese procedimiento no es totalmente trivial, pero no demasiado duro (para un diseño equilibrado como este).