Pourquoi ces 2 approches d'application de modèles mixtes donnent-elles des résultats différents?

Je ré-analyse les données d'un collègue. Les données et le code R sont ici .

Il s'agit d'une conception 2x2x2x2x3 complètement intra-SS. L'une des variables prédictives cue, est une variable à deux niveaux qui, lorsqu'elle est réduite à un score de différence, reflète une valeur pertinente pour la théorie. Elle s'est auparavant effondrée cueà un score de différence dans chaque sujet et condition, puis a calculé une ANOVA, produisant un MSE qu'elle pourrait ensuite utiliser pour les comparaisons prévues du score de différence moyenne de chaque condition par rapport à zéro. Vous devrez me faire confiance qu'elle ne pêchait pas et avait en effet une bonne base théorique pour faire les 24 tests.

J'ai pensé voir s'il y avait une différence en utilisant des modèles à effets mixtes pour représenter les données. Comme indiqué dans le code, j'ai adopté deux approches:

Méthode 1 - Modéliser les données sous la forme d'un plan 2x2x2x2x3, obtenir des échantillons a posteriori de ce modèle, calculer le cuescore de différence pour chaque condition dans chaque échantillon, calculer l'intervalle de prédiction de 95% pour le score de différence de repère dans chaque condition.

Méthode 2 - Réduire cueà un score de différence dans chaque sujet et condition, modéliser les données sous la forme d'un plan 2x2x2x3, obtenir des échantillons a posteriori de ce modèle, calculer l'intervalle de prédiction de 95% pour le score de différence de repère dans chaque condition.

Il semble que la méthode 1 donne des intervalles de prédiction plus larges que la méthode 2, avec pour conséquence que si l'on utilise le chevauchement avec zéro comme critère de «signification», seuls 25% des scores de cuing sont «significatifs» selon la méthode 1 tandis que 75% des scores de cuing sont "significatifs" selon la méthode 2. De manière notable, les profils de signification obtenus par la méthode 2 sont plus proches des résultats originaux basés sur l'ANOVA que les modèles obtenus par la méthode 1.

Une idée de ce qui se passe ici?

Il n'est pas surprenant de voir une telle différence avec lmer ou lme. Un modèle simple avec une interception aléatoire (par exemple, (1 | id) dans votre cas) peut parfois ne pas capturer complètement les effets aléatoires. Pour voir pourquoi cela se produit, permettez-moi d'utiliser un ensemble de données beaucoup plus simple que le vôtre pour démontrer la différence subtile. Avec les données 'dat' du fil que je copie ici:

dat <- structure(list(sex = structure(c(1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L), .Label = c("f",
"m"), class = "factor"), prevalence = c(0, 0.375, 0.133333333333333,
0.176470588235294, 0.1875, 0, 0, 1, 1, 0.5, 0.6, 0.333333333333333,
0.5, 0, 0.333333333333333, 0, 0.5, 0, 0.625, 0.333333333333333,
0.5, 0, 0.333333333333333, 0.153846153846154, 0.222222222222222,
0.5, 1, 0.5, 0, 0.277777777777778, 0.125, 0, 0, 0.428571428571429,
0.451612903225806, 0.362068965517241), tripsite = structure(c(1L,
1L, 4L, 4L, 14L, 14L, 5L, 5L, 8L, 8L, 15L, 15L, 6L, 6L, 9L, 9L,
11L, 11L, 16L, 16L, 2L, 2L, 7L, 7L, 10L, 10L, 13L, 13L, 17L,
17L, 3L, 3L, 12L, 12L, 18L, 18L), .Label = c("1.2", "4.2", "5.2",
"1.3", "2.3", "3.3", "4.3", "2.4", "3.4", "4.4", "3.5", "5.5",
"4.6", "1.9", "2.9", "3.9", "4.9", "5.9"), class = "factor")), .Names =
c("sex","prevalence", "tripsite"), row.names = c(1L, 2L, 3L, 4L, 9L,
10L, 11L, 12L, 13L, 14L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L,
27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 38L, 39L, 40L,
41L, 42L, 43L, 45L, 46L), class = "data.frame")

un test t apparié (ou un cas particulier d'ANOVA unidirectionnel à l'intérieur d'un sujet / mesures répétées) serait comme votre méthode 2:

t0 <- with(dat,t.test(prevalence[sex=="f"],prevalence[sex=="m"],paired=TRUE,var.equal=TRUE))
(fstat0 <- t0$statistic^2)         #0.789627

Sa version lme correspondant à votre méthode 1 serait:

a1 <- anova(lme(prevalence~sex,random=~1|tripsite,data=dat,method="REML"))
(fstat1 <- a1[["F-value"]][2])   # 0.8056624

Même chose pour l'homologue lmer:

a2 <- anova(lmer(prevalence~sex+(1|tripsite), data=dat))
(fstat2 <- a2[["F value"]][2])  # 0.8056624

Bien que la différence avec cet exemple simple soit minime, elle montre que le test t apparié a une hypothèse beaucoup plus forte sur les deux niveaux ("f" et "m") du facteur ("sexe"), que les deux niveaux sont corrélés et cette hypothèse est absente dans le modèle lme / lmer ci-dessus. Une telle différence d'hypothèse existe également entre les deux méthodes dans votre cas.

Pour réconcilier la différence, on peut continuer à modéliser 'dat' avec une pente aléatoire (ou matrice symétrique ou même symétrie composée) en lme / lmer:

a3 <- anova(lme(prevalence~sex,random=~sex-1|tripsite,data=dat,method="REML"))
(fstat3 <- a3[["F-value"]][2]) # 0.789627

a31 <- anova(lme(prevalence~sex,random=list(tripsite=pdCompSymm(~sex-1)),data=dat,method="REML")))
(fstat31 <- a31[["F-value"]][2]) # 0.789627

a4 <- anova(lmer(prevalence~sex+(sex-1|tripsite), data=dat))
(fstat4 <- a4[["F value"]][2]) # 0.789627

Cependant, avec plusieurs facteurs dans votre cas, plusieurs pentes aléatoires (ou d'autres spécifications de structure à effets aléatoires) peuvent devenir difficiles à utiliser avec lme / lmer si ce n'est impossible.