Équivalence entre un modèle anova à mesures répétées et un modèle mixte: lmer vs lme, et symétrie composée

J'ai du mal à obtenir des résultats équivalents entre un aovmodèle de mesures répétées inter-intra et un lmermodèle mixte.

Mes données et mon script se présentent comme suit

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

Fondamentalement, 3 x 6 sujets ( id) ont été soumis à trois FITNESSprogrammes d'entraînement différents chacun et leur a PULSEété mesuré après avoir effectué trois types différents d'endurance TEST.

J'ai ensuite monté le aovmodèle suivant :

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

Le résultat que j'obtiens en utilisant

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

est identique à cela.

Un modèle mixte exécuté en utilisant nlmedonne un résultat directement équivalent, par exemple en utilisant 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

Ou en utilisant gls:

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

Cependant, le résultat que j'obtenir l' aide lme4« s lmerest différent:

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

Quelqu'un pense-t-il à ce que je fais mal avec le lmermodèle? Ou d'où vient la différence? Pourrait-il avoir à faire quoi que ce soit lmersans autoriser les corellations intraclasses négatives ou quelque chose comme ça? Étant donné que nlmel » glset lmene retournent le résultat correct, cependant, je me demande comment cela est différent glset lme? Est-ce que l'option les correlation=corCompSymm(form=~1|id)amène à estimer directement la corrélation intraclasse, qui peut être positive ou négative, alors qu'elle lmerestime une composante de variance, qui ne peut pas être négative (et finit par être estimée nulle dans ce cas)?

Comme Ben Bolker l'a déjà mentionné dans les commentaires, le problème est comme vous le pensez: le lmer()modèle est déclenché parce qu'il tente d'estimer un modèle de composante de variance, les estimations de la composante de variance étant contraintes d'être non négatives. Ce que j'essaierai de faire est de donner une compréhension quelque peu intuitive de ce qui concerne votre ensemble de données qui mène à cela, et pourquoi cela pose un problème pour les modèles de composants de variance.

Voici un tracé de votre jeu de données. Les points blancs sont les observations réelles et les points noirs sont les moyens sujets.

Pour rendre les choses plus simples, mais sans changer l'esprit du problème, je vais soustraire les effets fixes (c'est-à-dire les effets FITNESSet TEST, ainsi que la moyenne générale) et traiter les données résiduelles comme un problème d'effets aléatoires à sens unique . Voici donc à quoi ressemble le nouvel ensemble de données:

Examinez attentivement les motifs de cette intrigue. Réfléchissez à la façon dont les observations prises sur le même sujet diffèrent des observations prises sur différents sujets. Plus précisément, notez le schéma suivant: comme l'une des observations pour un sujet est supérieure (ou inférieure) au-dessus (ou en dessous) de la moyenne du sujet, les autres observations de ce sujet ont tendance à être du côté opposé de la moyenne du sujet. Et plus cette observation est éloignée de la moyenne du sujet, plus les autres observations tendent à être éloignées de la moyenne du sujet de l'autre côté. Cela indique une corrélation intra-classe négative. En fait, deux observations tirées du même sujet ont tendance à être moins similaires, en moyenne, que deux observations tirées uniquement au hasard de l'ensemble de données.

Une autre façon de penser à ce modèle est en termes d'amplitudes relatives de la variance inter-sujet et intra-sujet. Il semble qu'il y ait une variance intra-sujet beaucoup plus importante que la variance inter-sujet. Bien sûr, nous nous attendons à ce que cela se produise dans une certaine mesure. Après tout, la variance intra-sujet est basée sur la variation des points de données individuels, alors que l'écart entre-sujet est basée sur la variation du moyen des points de données individuels (les moyens en question), et nous savons que la variance une moyenne aura tendance à diminuer à mesure que le nombre de choses en moyenne augmente. Mais dans cet ensemble de données, la différence est assez frappante: il y a moyenplus de variation intra-sujet qu'entre sujets. En fait, cette différence est exactement la raison pour laquelle une corrélation intra-classe négative émerge.

D'accord, voici donc le problème. Le modèle de composante de variance suppose que chaque point de données est la somme d'un effet sujet et d'une erreur:$y_{ij}=u_j+e_{ij}$, où $u_j$ est l'effet de la $j$e sujet. Réfléchissons donc à ce qui se passerait s'il y avait vraiment 0 variance dans les effets du sujet - en d'autres termes, si la véritable composante de la variance entre les sujets était 0. Étant donné un ensemble de données réel généré sous ce modèle, si nous devions calculer des moyennes d'échantillon pour les données observées de chaque sujet, ces moyennes d'échantillon auraient toujours une variance non nulle, mais elles ne reflèteraient que la variance d'erreur, et non aucune «vraie» variance du sujet (parce que nous avons supposé qu'il n'y en avait pas).

Alors, dans quelle mesure devrions-nous nous attendre à ce que ces moyens sujets soient? Eh bien, fondamentalement, chaque effet sujet estimé est une moyenne, et nous connaissons la formule de la variance d'une moyenne:$\text{var}(\bar{X})=\text{var}(X_i)/n$, où $n$est le nombre de choses en moyenne. Appliquons maintenant cette formule à votre ensemble de données et voyons combien de variance nous nous attendrions à voir dans les effets de sujet estimés si la véritable composante de variance entre sujets était exactement 0.

La variance intra-sujet se révèle être $348$et chaque effet sujet est calculé comme la moyenne de 3 observations. Ainsi, l'écart type attendu dans le sujet signifie - en supposant que la vraie variance entre les sujets est 0 - se révèle être d'environ$10.8$. Maintenant, comparez ceci à l'écart-type dans le sujet signifie que nous avons effectivement observé:$4.3$! La variation observée est nettement inférieure à la variation attendue lorsque nous avons supposé 0 variance inter-sujet. Pour un modèle à composante de variance, la seule façon dont on pourrait s'attendre à ce que la variation observée soit aussi faible que ce que nous avons réellement observé est si la vraie variance inter-sujet était en quelque sorte négative . Et c'est là que réside le problème. Les données impliquent qu'il existe en quelque sorte une composante de variance négative, mais le logiciel (raisonnablement) ne permettra pas d'estimations négatives des composantes de variance, car une variance ne peut en fait jamais être négative. Les autres modèles que vous ajustez évitent ce problème en estimant directement la corrélation intra-classe plutôt qu'en supposant un modèle de composante de variance simple.

Si vous voulez voir comment vous pourriez réellement obtenir l'estimation de la composante de variance négative impliquée par votre ensemble de données, vous pouvez utiliser la procédure que j'illustre (avec le Rcode d' accompagnement ) dans cette autre réponse récente . Cette procédure n'est pas totalement triviale, mais pas trop difficile non plus (pour une conception équilibrée comme celle-ci).