Quel est l'équivalent lme4 :: lmer d'une ANOVA à trois mesures répétées?

11

Ma question est basée sur cette réponse qui a montré quel lme4::lmermodèle correspond à une mesure répétée bidirectionnelle ANOVA:

require(lme4)
set.seed(1234)
d <- data.frame(
    y = rnorm(96),
    subject = factor(rep(1:12, 4)),
    a = factor(rep(1:2, each=24)),
    b = factor(rep(rep(1:2, each=12))),
    c = factor(rep(rep(1:2, each=48))))

# standard two-way repeated measures ANOVA:
summary(aov(y~a*b+Error(subject/(a*b)), d[d$c == "1",]))

# corresponding lmer call:
anova(lmer(y ~ a*b+(1|subject) + (1|a:subject) + (1|b:subject), d[d$c == "1",]))

Ma question est maintenant de savoir comment étendre cela au cas d'une ANOVA à trois voies:

summary(aov(y~a*b*c+Error(subject/(a*b*c)), d))
## [...]
## Error: subject:a:b:c
##           Df Sum Sq Mean Sq F value Pr(>F)
## a:b:c      1  0.101  0.1014   0.115  0.741
## Residuals 11  9.705  0.8822

L'extension naturelle ainsi que ses versions ne correspondent pas aux résultats de l'ANOVA:

anova(lmer(y ~ a*b*c +(1|subject) + (1|a:subject) + (1|b:subject) + (1|c:subject), d))
## [...]
## a:b:c  1 0.1014  0.1014  0.1500

anova(lmer(y ~ a*b*c +(1|subject) + (1|a:subject) + (1|b:subject) + (1|c:subject) + 
               (1|a:b:subject) + (1|a:c:subject) + (1|b:c:subject), d))
## [...]
## a:b:c  1 0.1014  0.1014  0.1539

Notez que a été posé une question similaire avant . Cependant, il manquait des exemples de données (qui sont fournis ici).

r anova mixed-model repeated-measures lme4-nlme Henrik
la source
Êtes-vous sûr de ne pas vouloir que votre modèle à deux niveaux soit y ~ a*b + (1 + a*b|subject), d[d$c == "1",]? Ou peut-être que je manque quelque chose?
Rasmus Bååth
@ RasmusBååth Go Ahead et essayez de l'adapter, lmerse plaindra car les effets aléatoires ne sont plus identifiés. Au départ, je pensais aussi que c'était le modèle que je voulais, mais ce n'est pas le cas. Si vous comparez le modèle lmer que je propose pour le cas à 2 voies avec l'ANOVA standard, vous verrez que les valeurs F correspondent exactement . Comme indiqué dans la réponse, j'ai lié.
Henrik
3
Pour le problème à trois voies, le premier lmermodèle que vous avez écrit (qui exclut les interactions bidirectionnelles aléatoires) ne devrait pas être équivalent à une RM-ANOVA à 3 voies, mais le second que vous avez écrit (qui inclut le hasard les interactions bidirectionnelles) devraient l'être. Quant à savoir pourquoi il y a un écart, même avec ce modèle, j'ai une idée de ce qu'est le problème.Je vais prendre le dîner puis j'examinerai un peu plus l'ensemble de données sur les jouets.
Jake Westfall

Réponses:

18

La réponse directe à votre question est que le dernier modèle que vous avez écrit,

anova(lmer(y ~ a*b*c +(1|subject) + (1|a:subject) + (1|b:subject) + (1|c:subject) + 
           (1|a:b:subject) + (1|a:c:subject) + (1|b:c:subject), d))

Je pense que c'est "en principe" correct, bien qu'il s'agisse d'un paramétrage étrange qui ne semble pas toujours bien fonctionner dans la pratique.

Quant à savoir pourquoi la sortie que vous obtenez de ce modèle est différente de la aov()sortie, je pense qu'il y a deux raisons.

  1. Votre jeu de données simulé simple est pathologique dans la mesure où le modèle le mieux adapté est celui qui implique des composantes de variance négatives, auxquelles les modèles mixtes s'ajustent lmer()(et la plupart des autres programmes de modèles mixtes) ne le permettront pas.
  2. Même avec un ensemble de données non pathologique, la façon dont vous avez configuré le modèle, comme mentionné ci-dessus, ne semble pas toujours bien fonctionner dans la pratique, bien que je dois admettre que je ne comprends pas vraiment pourquoi. C'est aussi généralement étrange à mon avis, mais c'est une autre histoire.

Permettez-moi d'abord de démontrer le paramétrage que je préfère sur votre exemple ANOVA bidirectionnel initial. Supposons que votre jeu de données dest chargé. Votre modèle (notez que j'ai changé de code factice en code de contraste) était:

options(contrasts=c("contr.sum","contr.poly"))
mod1 <- lmer(y ~ a*b+(1|subject) + (1|a:subject) + (1|b:subject),
         data = d[d$c == "1",])
anova(mod1)
# Analysis of Variance Table
#     Df  Sum Sq Mean Sq F value
# a    1 2.20496 2.20496  3.9592
# b    1 0.13979 0.13979  0.2510
# a:b  1 1.23501 1.23501  2.2176

qui a bien fonctionné ici en ce qu'il correspondait à la aov()sortie. Le modèle que je préfère implique deux changements: coder manuellement les facteurs de contraste afin que nous ne travaillions pas avec des objets de facteur R (ce que je recommande de faire dans 100% des cas), et spécifier les effets aléatoires différemment:

d <- within(d, {
  A <- 2*as.numeric(paste(a)) - 3
  B <- 2*as.numeric(paste(b)) - 3
  C <- 2*as.numeric(paste(c)) - 3
})
mod2 <- lmer(y ~ A*B + (1|subject)+(0+A|subject)+(0+B|subject),
             data = d[d$c == "1",])
anova(mod2)
# Analysis of Variance Table
# Df  Sum Sq Mean Sq F value
# A    1 2.20496 2.20496  3.9592
# B    1 0.13979 0.13979  0.2510
# A:B  1 1.23501 1.23501  2.2176

logLik(mod1)
# 'log Lik.' -63.53034 (df=8)
logLik(mod2)
# 'log Lik.' -63.53034 (df=8)

Les deux approches sont totalement équivalentes dans le problème simple à 2 voies. Nous allons maintenant passer à un problème à 3 voies. J'ai mentionné plus tôt que l'exemple de jeu de données que vous avez donné était pathologique. Donc, ce que je veux faire avant d'aborder votre exemple de jeu de données est de générer d'abord un jeu de données à partir d'un modèle de composants de variance réel (c'est-à-dire où les composants de variance non nuls sont intégrés dans le vrai modèle). Je vais d'abord montrer comment mon paramétrage préféré semble mieux fonctionner que celui que vous avez proposé. Je montrerai ensuite une autre façon d'estimer les composantes de la variance qui n'impose pas qu'elles doivent être non négatives. Ensuite, nous serons en mesure de voir le problème avec l'exemple de jeu de données d'origine.

Le nouvel ensemble de données sera de structure identique, sauf que nous aurons 50 sujets:

set.seed(9852903)
d2 <- expand.grid(A=c(-1,1), B=c(-1,1), C=c(-1,1), sub=seq(50))
d2 <- merge(d2, data.frame(sub=seq(50), int=rnorm(50), Ab=rnorm(50),
  Bb=rnorm(50), Cb=rnorm(50), ABb=rnorm(50), ACb=rnorm(50), BCb=rnorm(50)))
d2 <- within(d2, {
  y <- int + (1+Ab)*A + (1+Bb)*B + (1+Cb)*C + (1+ABb)*A*B +
    (1+ACb)*A*C + (1+BCb)*B*C + A*B*C + rnorm(50*2^3)
  a <- factor(A)
  b <- factor(B)
  c <- factor(C)
})

Les ratios F auxquels nous voulons correspondre sont:

aovMod1 <- aov(y ~ a*b*c + Error(factor(sub)/(a*b*c)), data = d2)
tab <- lapply(summary(aovMod1), function(x) x[[1]][1,2:4])
do.call(rbind, tab)
#                          Sum Sq Mean Sq F value
# Error: factor(sub)       439.48    8.97        
# Error: factor(sub):a     429.64  429.64  32.975
# Error: factor(sub):b     329.48  329.48  27.653
# Error: factor(sub):c     165.44  165.44  17.924
# Error: factor(sub):a:b   491.33  491.33  49.694
# Error: factor(sub):a:c   305.46  305.46  41.703
# Error: factor(sub):b:c   466.09  466.09  40.655
# Error: factor(sub):a:b:c 392.76  392.76 448.101

Voici nos deux modèles:

mod3 <- lmer(y ~ a*b*c + (1|sub)+(1|a:sub)+(1|b:sub)+(1|c:sub)+
  (1|a:b:sub)+(1|a:c:sub)+(1|b:c:sub), data = d2)
anova(mod3)
# Analysis of Variance Table
#       Df Sum Sq Mean Sq F value
# a      1  32.73   32.73  34.278
# b      1  21.68   21.68  22.704
# c      1  12.53   12.53  13.128
# a:b    1  60.93   60.93  63.814
# a:c    1  50.38   50.38  52.762
# b:c    1  57.30   57.30  60.009
# a:b:c  1 392.76  392.76 411.365

mod4 <- lmer(y ~ A*B*C + (1|sub)+(0+A|sub)+(0+B|sub)+(0+C|sub)+
  (0+A:B|sub)+(0+A:C|sub)+(0+B:C|sub), data = d2)
anova(mod4)
# Analysis of Variance Table
#       Df Sum Sq Mean Sq F value
# A      1  28.90   28.90  32.975
# B      1  24.24   24.24  27.653
# C      1  15.71   15.71  17.924
# A:B    1  43.56   43.56  49.694
# A:C    1  36.55   36.55  41.703
# B:C    1  35.63   35.63  40.655
# A:B:C  1 392.76  392.76 448.101

logLik(mod3)
# 'log Lik.' -984.4531 (df=16)
logLik(mod4)
# 'log Lik.' -973.4428 (df=16)

Comme nous pouvons le voir, seule la deuxième méthode correspond à la sortie de aov(), bien que la première méthode soit au moins dans le stade approximatif. La deuxième méthode permet également d'obtenir une log-vraisemblance plus élevée. Je ne sais pas pourquoi ces deux méthodes donnent des résultats différents, car là encore je pense qu'elles sont "en principe" équivalentes, mais c'est peut-être pour des raisons numériques / informatiques. Ou peut-être que je me trompe et qu'ils ne sont pas équivalents, même en principe.

Je vais maintenant montrer une autre façon d'estimer les composantes de la variance sur la base des idées ANOVA traditionnelles. Fondamentalement, nous prendrons les équations des carrés moyens attendus pour votre conception, remplacerons les valeurs observées des carrés moyens et résoudrons les composantes de la variance. Pour obtenir les carrés moyens attendus, nous utiliserons une fonction R que j'ai écrite il y a quelques années, appelée EMS(), qui est documentée ICI . Ci-dessous, je suppose que la fonction est déjà chargée.

# prepare coefficient matrix
r <- 1 # number of replicates
s <- 50 # number of subjects
a <- 2 # number of levels of A
b <- 2 # number of levels of B
c <- 2 # number of levels of C
CT <- EMS(r ~ a*b*c*s, random="s")
expr <- strsplit(CT[CT != ""], split="")
expr <- unlist(lapply(expr, paste, collapse="*"))
expr <- sapply(expr, function(x) eval(parse(text=x)))
CT[CT != ""] <- expr
CT[CT == ""] <- 0
mode(CT) <- "numeric"
# residual variance and A*B*C*S variance are confounded in
# this design, so remove the A*B*C*S variance component
CT <- CT[-15,-2]
CT
#        VarianceComponent
# Effect  e b:c:s a:c:s a:b:s a:b:c c:s b:s a:s b:c a:c a:b s   c   b   a
#   a     1     0     0     0     0   0   0   4   0   0   0 0   0   0 200
#   b     1     0     0     0     0   0   4   0   0   0   0 0   0 200   0
#   c     1     0     0     0     0   4   0   0   0   0   0 0 200   0   0
#   s     1     0     0     0     0   0   0   0   0   0   0 8   0   0   0
#   a:b   1     0     0     2     0   0   0   0   0   0 100 0   0   0   0
#   a:c   1     0     2     0     0   0   0   0   0 100   0 0   0   0   0
#   b:c   1     2     0     0     0   0   0   0 100   0   0 0   0   0   0
#   a:s   1     0     0     0     0   0   0   4   0   0   0 0   0   0   0
#   b:s   1     0     0     0     0   0   4   0   0   0   0 0   0   0   0
#   c:s   1     0     0     0     0   4   0   0   0   0   0 0   0   0   0
#   a:b:c 1     0     0     0    50   0   0   0   0   0   0 0   0   0   0
#   a:b:s 1     0     0     2     0   0   0   0   0   0   0 0   0   0   0
#   a:c:s 1     0     2     0     0   0   0   0   0   0   0 0   0   0   0
#   b:c:s 1     2     0     0     0   0   0   0   0   0   0 0   0   0   0
#   e     1     0     0     0     0   0   0   0   0   0   0 0   0   0   0

# get mean squares
(MSmod <- summary(aov(y ~ a*b*c*factor(sub), data=d2)))
#                   Df Sum Sq Mean Sq
# a                  1  429.6   429.6
# b                  1  329.5   329.5
# c                  1  165.4   165.4
# factor(sub)       49  439.5     9.0
# a:b                1  491.3   491.3
# a:c                1  305.5   305.5
# b:c                1  466.1   466.1
# a:factor(sub)     49  638.4    13.0
# b:factor(sub)     49  583.8    11.9
# c:factor(sub)     49  452.2     9.2
# a:b:c              1  392.8   392.8
# a:b:factor(sub)   49  484.5     9.9
# a:c:factor(sub)   49  358.9     7.3
# b:c:factor(sub)   49  561.8    11.5
# a:b:c:factor(sub) 49   42.9     0.9
MS <- MSmod[[1]][,"Mean Sq"]

# solve
ans <- solve(CT, MS)
cbind(rev(ans[c(grep("e",names(ans)),grep("s",names(ans)))])/
        c(1,2,2,2,4,4,4,1))
# s     1.0115549
# a:s   1.5191114
# b:s   1.3797937
# c:s   1.0441351
# a:b:s 1.1263331
# a:c:s 0.8060402
# b:c:s 1.3235126
# e     0.8765093
summary(mod4)
# Random effects:
#  Groups   Name        Variance Std.Dev.
#  sub      (Intercept) 1.0116   1.0058  
#  sub.1    A           1.5191   1.2325  
#  sub.2    B           1.3798   1.1746  
#  sub.3    C           1.0441   1.0218  
#  sub.4    A:B         1.1263   1.0613  
#  sub.5    A:C         0.8060   0.8978  
#  sub.6    B:C         1.3235   1.1504  
#  Residual             0.8765   0.9362  
# Number of obs: 400, groups:  sub, 50

D'accord, nous allons maintenant revenir à l'exemple d'origine. Les ratios F que nous essayons de faire correspondre sont les suivants:

aovMod2 <- aov(y~a*b*c+Error(subject/(a*b*c)), data = d)
tab <- lapply(summary(aovMod2), function(x) x[[1]][1,2:4])
do.call(rbind, tab)
#                       Sum Sq Mean Sq F value
# Error: subject       13.4747  1.2250        
# Error: subject:a      1.4085  1.4085  1.2218
# Error: subject:b      3.1180  3.1180  5.5487
# Error: subject:c      6.3809  6.3809  5.2430
# Error: subject:a:b    1.5706  1.5706  2.6638
# Error: subject:a:c    1.0907  1.0907  1.5687
# Error: subject:b:c    1.4128  1.4128  2.3504
# Error: subject:a:b:c  0.1014  0.1014  0.1149

Voici nos deux modèles:

mod5 <- lmer(y ~ a*b*c + (1|subject)+(1|a:subject)+(1|b:subject)+
  (1|c:subject)+(1|a:b:subject)+(1|a:c:subject)+(1|b:c:subject),
  data = d)
anova(mod5)
# Analysis of Variance Table
#       Df Sum Sq Mean Sq F value
# a      1 0.8830  0.8830  1.3405
# b      1 3.1180  3.1180  4.7334
# c      1 3.8062  3.8062  5.7781
# a:b    1 1.5706  1.5706  2.3844
# a:c    1 0.9620  0.9620  1.4604
# b:c    1 1.4128  1.4128  2.1447
# a:b:c  1 0.1014  0.1014  0.1539

mod6 <- lmer(y ~ A*B*C + (1|subject)+(0+A|subject)+(0+B|subject)+
  (0+C|subject)+(0+A:B|subject)+(0+A:C|subject)+
  (0+B:C|subject), data = d)
anova(mod6)
# Analysis of Variance Table
#       Df Sum Sq Mean Sq F value
# a      1 0.8830  0.8830  1.3405
# b      1 3.1180  3.1180  4.7334
# c      1 3.8062  3.8062  5.7781
# a:b    1 1.5706  1.5706  2.3844
# a:c    1 0.9620  0.9620  1.4604
# b:c    1 1.4128  1.4128  2.1447
# a:b:c  1 0.1014  0.1014  0.1539

logLik(mod5)
# 'log Lik.' -135.0351 (df=16)
logLik(mod6)
# 'log Lik.' -134.9191 (df=16)

Dans ce cas, les deux modèles donnent essentiellement les mêmes résultats, bien que la deuxième méthode ait une probabilité logarithmique très légèrement plus élevée. Aucune méthode ne correspond aov(). Mais regardons ce que nous obtenons lorsque nous résolvons les composants de variance comme nous l'avons fait ci-dessus, en utilisant la procédure ANOVA qui ne contraint pas les composants de variance à être non négatifs (mais qui ne peuvent être utilisés que dans des conceptions équilibrées sans prédicteurs continus et sans données manquantes; hypothèses classiques de l'ANOVA).

# prepare coefficient matrix
r <- 1 # number of replicates
s <- 12 # number of subjects
a <- 2 # number of levels of A
b <- 2 # number of levels of B
c <- 2 # number of levels of C
CT <- EMS(r ~ a*b*c*s, random="s")
expr <- strsplit(CT[CT != ""], split="")
expr <- unlist(lapply(expr, paste, collapse="*"))
expr <- sapply(expr, function(x) eval(parse(text=x)))
CT[CT != ""] <- expr
CT[CT == ""] <- 0
mode(CT) <- "numeric"
# residual variance and A*B*C*S variance are confounded in
# this design, so remove the A*B*C*S variance component
CT <- CT[-15,-2]

# get mean squares
MSmod <- summary(aov(y ~ a*b*c*subject, data=d))
MS <- MSmod[[1]][,"Mean Sq"]

# solve
ans <- solve(CT, MS)
cbind(rev(ans[c(grep("e",names(ans)),grep("s",names(ans)))])/
        c(1,2,2,2,4,4,4,1))
# s      0.04284033
# a:s    0.03381648
# b:s   -0.04004005
# c:s    0.04184887
# a:b:s -0.03657940
# a:c:s -0.02337501
# b:c:s -0.03514457
# e      0.88224787
summary(mod6)
# Random effects:
#  Groups    Name        Variance  Std.Dev. 
#  subject   (Intercept) 7.078e-02 2.660e-01
#  subject.1 A           6.176e-02 2.485e-01
#  subject.2 B           0.000e+00 0.000e+00
#  subject.3 C           6.979e-02 2.642e-01
#  subject.4 A:B         1.549e-16 1.245e-08
#  subject.5 A:C         4.566e-03 6.757e-02
#  subject.6 B:C         0.000e+00 0.000e+00
#  Residual              6.587e-01 8.116e-01
# Number of obs: 96, groups:  subject, 12

Maintenant, nous pouvons voir ce qui est pathologique dans l'exemple original. Le modèle le mieux adapté est celui qui implique que plusieurs des composantes de la variance aléatoire sont négatives. Mais lmer()(et la plupart des autres programmes de modèles mixtes) contraint les estimations des composantes de la variance à être non négatives. Ceci est généralement considéré comme une contrainte sensible, car les variances ne peuvent bien sûr jamais vraiment être négatives. Cependant, une conséquence de cette contrainte est que les modèles mixtes sont incapables de représenter avec précision les ensembles de données qui présentent des corrélations intraclasses négatives, c'est-à-dire les ensembles de données où les observations du même cluster sont moins(plutôt que plus) similaires en moyenne aux observations tirées au hasard de l'ensemble de données, et par conséquent lorsque la variance intra-cluster dépasse considérablement la variance inter-cluster. Ces ensembles de données sont des ensembles de données parfaitement raisonnables que l'on rencontrera occasionnellement dans le monde réel (ou simuler accidentellement!), Mais ils ne peuvent pas être décrits de manière raisonnable par un modèle à composantes de variance, car ils impliquent des composantes de variance négatives. Ils peuvent cependant être décrits de manière «non sensible» par de tels modèles, si le logiciel le permet. aov()le permet. lmer()ne fait pas.

Jake Westfall
la source
+1. Re I am not sure why these two methods give different results, as again I think they are "in principle" equivalent, but maybe it is for some numerical/computational reasons- comprenez-vous peut-être mieux cela maintenant (deux ans plus tard)? J'ai essayé de comprendre quelle est la différence, mais je ne comprends pas non plus ...
Amoeba dit Réintégrer Monica
@amoeba Ma pensée actuelle est à peu près la même qu'alors: AFAIK, les deux modèles sont statistiquement équivalents (dans le sens où ils font les mêmes prédictions sur les données et impliquent les mêmes erreurs standard), même si les effets aléatoires sont paramétrés différemment. Je pense que les différences observées - qui semblent ne se produire que parfois - sont uniquement dues à des problèmes de calcul. En particulier, je soupçonne que vous pourriez jouer avec les paramètres d'optimisation (comme varier les points de départ ou utiliser des critères de convergence plus stricts) jusqu'à ce que les deux modèles retournent exactement la même réponse.
Jake Westfall
Merci pour votre réponse. Je ne suis pas assez convaincu: j'ai essayé de jouer avec les paramètres de l'optimiseur et je n'ai rien pu changer dans les résultats; mon impression est que les deux modèles sont bien convergés. Je pourrais peut-être poser une question distincte à un moment donné.
amibe dit Réintégrer Monica
Je continue à enquêter sur ce problème et j'ai réalisé ce qui suit: dès que les facteurs de mesures répétées ont plus de deux niveaux, ces deux méthodes deviennent entièrement différentes! Si Aa niveaux, alors estime seulement un paramètre de variance, alors qu'il estime paramètres de variance et corrélations entre eux. D'après ce que je comprends, l'ANOVA classique estime un seul paramètre de variance et devrait donc être équivalente à la première méthode, pas à la seconde. Droite? Mais pour deux méthodes estiment un paramètre, et je ne sais toujours pas pourquoi elles ne sont pas d'accord. k - 1 k ( k - 1 ) / 2 k = 2k(1|A:sub)(0+A|sub)k−1k(k−1)/2k=2
amibe dit Réintégrer Monica
Pour en revenir à ce problème ... J'ai remarqué que pour le cas à deux facteurs où deux lmerappels produisent une anova()sortie identique , les variances d'effets aléatoires sont néanmoins assez différentes: voir VarCorr(mod1)et VarCorr(mod2). Je ne comprends pas très bien pourquoi cela se produit; le faites vous? Pour mod3et mod4, on peut voir que quatre des sept variances pour mod3sont en fait égales à zéro (pour mod4les sept sont non nulles); cette "singularité" en mod3est probablement la raison pour laquelle les tables anova diffèrent. En dehors de cela, comment utiliseriez-vous votre "voie préférée" si aet bavait plus de deux niveaux?
amibe dit Réintégrer Monica le
1

Sont a, b, cfixes ou effets aléatoires? S'ils sont corrigés, votre syntaxe sera simplement

summary(aov(y~a*b*c+Error(subject), d))
Error: subject
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 11  13.47   1.225               

Error: Within
          Df Sum Sq Mean Sq F value  Pr(>F)   
a          1   1.41   1.408   1.730 0.19235   
b          1   3.12   3.118   3.829 0.05399 . 
c          1   6.38   6.381   7.836 0.00647 **
a:b        1   1.57   1.571   1.929 0.16889   
a:c        1   1.09   1.091   1.339 0.25072   
b:c        1   1.41   1.413   1.735 0.19168   
a:b:c      1   0.10   0.101   0.124 0.72518   
Residuals 77  62.70   0.814  

library(lmerTest)
anova(lmer(y ~ a*b*c+(1|subject), data=d))
Analysis of Variance Table of type 3  with  Satterthwaite 
approximation for degrees of freedom
      Sum Sq Mean Sq NumDF  DenDF F.value   Pr(>F)   
a     1.4085  1.4085     1 76.991  1.7297 0.192349   
b     3.1180  3.1180     1 76.991  3.8291 0.053995 . 
c     6.3809  6.3809     1 76.991  7.8363 0.006469 **
a:b   1.5706  1.5706     1 76.991  1.9289 0.168888   
a:c   1.0907  1.0907     1 76.991  1.3394 0.250716   
b:c   1.4128  1.4128     1 76.991  1.7350 0.191680   
a:b:c 0.1014  0.1014     1 76.991  0.1245 0.725183
Masato Nakazawa
la source
Ce sont des effets fixes. Cependant, le modèle ANOVA que vous ajustez n'est pas le modèle qui semble être le modèle ANOVA à mesures répétées classiques, voir par exemple ici . Voir les strates d'erreur dans votre et mon cas.
Henrik
1
En fait, la façon dont ils le font est incorrecte. Si vous avez un plan de mesures répétitives factorielles entièrement croisé (ou plan factoriel à blocs randomisés), vous ne devriez obtenir qu'un seul terme d'erreur, à part subject, pour tous les effets (c.-à Within-d.). Voir Conception expérimentale: Procédures pour les sciences du comportement (2013) par Kirk, chapitre 10 (p.458) ou mon article ici
Masato Nakazawa
Laissons de côté cette question pour le moment et supposons que le modèle que j'ai adapté serait le bon modèle. Comment adapteriez-vous cela en utilisant lmer? J'obtiendrai néanmoins ma copie de Kirk (2e édition seulement) et verrai ce qu'elle dit.
Henrik
Je suis curieux de savoir ce que vous pensez du chapitre de Kirk. Je pense que le numéro de chapitre dans le 2e éd. est différent. En attendant, je vais essayer d'adapter différents lmermodèles. La meilleure façon de vérifier l'ajustement du modèle est de vérifier leur dfs en utilisant lmerTestparce que l'approximation KR est censée vous donner des exactdfs et donc des valeurs de p.
Masato Nakazawa
1
J'ai la 2e édition de Kirk, où je crois que la discussion pertinente se trouve aux pages 443 et 449, qui traite d'un exemple à deux voies (et non à trois voies). Les carrés moyens attendus, en supposant ou non l'additivité de A et B, sont donnés en p. 447. En supposant que A et B sont fixes et que les sujets / blocs sont aléatoires, nous pouvons voir d'après les carrés moyens attendus répertoriés par Kirk sous "modèle non additif" que les tests de A, B et AB impliquent chacun des termes d'erreur différents, à savoir, les interactions pertinentes avec le bloc / le sujet. Le même principe s'étend au présent exemple à trois voies.
Jake Westfall