Was ist das lme4 :: lmer-Äquivalent einer Drei-Wege-ANOVA mit wiederholten Messungen?

11

Meine Frage basiert auf dieser Antwort, die zeigte, welches lme4::lmerModell einer Zwei-Wege-ANOVA mit wiederholten Messungen entspricht:

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",]))

Meine Frage ist nun, wie ich dies auf den Fall einer Drei-Wege-ANOVA ausweiten kann:

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

Die natürliche Erweiterung sowie deren Versionen stimmen nicht mit den ANOVA-Ergebnissen überein:

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

Beachten Sie, dass eine sehr ähnliche Frage gestellt wurde , bevor . Es fehlten jedoch Beispieldaten (die hier bereitgestellt werden).

r anova mixed-model repeated-measures lme4-nlme Henrik
quelle
Sind Sie sicher, dass Ihr zweistufiges Modell nicht sein soll y ~ a*b + (1 + a*b|subject), d[d$c == "1",]? Oder fehlt mir vielleicht etwas?
Rasmus Bååth
@ RasmusBååth Gehen Sie voran und versuchen Sie es zu passen, lmerwird sich beschweren, da die zufälligen Effekte nicht mehr identifiziert werden. Anfangs dachte ich auch, dass dies das Modell ist, das ich will, aber es ist nicht so. Wenn Sie das von mir für den 2-Wege-Fall vorgeschlagene Modell mit der Standard-ANOVA vergleichen, werden Sie feststellen, dass die F-Werte genau übereinstimmen. Wie in der Antwort gesagt, habe ich verlinkt.
Henrik
3
Für das Drei-Wege-Problem lmerwird nicht erwartet, dass das erste Modell, das Sie geschrieben haben (das die zufälligen Zwei-Wege-Interaktionen ausschließt) , einer 3-Wege-RM-ANOVA entspricht, sondern das zweite, das Sie geschrieben haben (das den Zufall enthält) wechselseitige Wechselwirkungen) sollte sein. Was die Gründe angeht, warum es selbst bei diesem Modell eine Diskrepanz gibt, habe ich eine Ahnung, was das Problem ist. Wenn ich zum Abendessen gehe, schaue ich mir den Spielzeugdatensatz noch einmal an.
Jake Westfall

Antworten:

18

Die direkte Antwort auf Ihre Frage lautet: Das letzte Modell, das Sie geschrieben haben:

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))

Ich glaube, es ist "im Prinzip" richtig, obwohl es eine seltsame Parametrisierung ist, die in der Praxis nicht immer gut zu funktionieren scheint.

aov()Ich denke, es gibt zwei Gründe, warum die Ausgabe, die Sie von diesem Modell erhalten, nicht mit der Ausgabe übereinstimmt.

  1. Ihr einfacher simulierter Datensatz ist insofern pathologisch, als das am besten passende Modell negative Varianzkomponenten impliziert, nach denen gemischte Modelle passen lmer() (und den meisten anderen gemischten Modellprogrammen) nicht zugelassen werden.
  2. Selbst mit einem nicht pathologischen Datensatz scheint die Art und Weise, wie Sie das Modell eingerichtet haben, wie oben erwähnt, in der Praxis nicht immer gut zu funktionieren, obwohl ich zugeben muss, dass ich nicht wirklich verstehe, warum. Es ist meiner Meinung nach auch nur allgemein seltsam, aber das ist eine andere Geschichte.

Lassen Sie mich zunächst die Parametrisierung demonstrieren, die ich in Ihrem ersten Zwei-Wege-ANOVA-Beispiel bevorzuge. Angenommen, Ihr Datensatz dist geladen. Ihr Modell (beachten Sie, dass ich von Dummy- zu Kontrastcodes gewechselt habe) war:

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

was hier gut funktionierte, da es mit der aov()Ausgabe übereinstimmte . Das Modell, das ich bevorzuge, beinhaltet zwei Änderungen: Manuelles Kontrastcodieren der Faktoren, damit wir nicht mit R-Faktor-Objekten arbeiten (was ich in 100% der Fälle empfehle) und unterschiedliche Angabe der Zufallseffekte:

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)

Die beiden Ansätze sind im einfachen Zwei-Wege-Problem völlig gleichwertig. Jetzt gehen wir zu einem 3-Wege-Problem über. Ich habe bereits erwähnt, dass der von Ihnen angegebene Beispieldatensatz pathologisch war. Bevor ich mich also mit Ihrem Beispieldatensatz befasse, möchte ich zunächst einen Datensatz aus einem tatsächlichen Varianzkomponentenmodell generieren (dh wenn Varianzkomponenten ungleich Null in das wahre Modell integriert sind). Zuerst werde ich zeigen, wie meine bevorzugte Parametrisierung besser zu funktionieren scheint als die von Ihnen vorgeschlagene. Dann werde ich einen anderen Weg zur Schätzung der Varianzkomponenten demonstrieren, der nicht impliziert, dass sie nicht negativ sein dürfen. Dann sind wir in der Lage, das Problem mit dem ursprünglichen Beispieldatensatz zu sehen.

Der neue Datensatz wird in seiner Struktur identisch sein, außer dass wir 50 Probanden haben werden:

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)
})

Die F-Verhältnisse, die wir anpassen möchten, sind:

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

Hier sind unsere beiden Modelle:

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)

Wie wir sehen können, stimmt nur die zweite Methode mit der Ausgabe von überein aov() , obwohl sich die erste Methode zumindest im Baseballstadion befindet. Die zweite Methode erreicht auch eine höhere Log-Wahrscheinlichkeit. Ich bin mir nicht sicher, warum diese beiden Methoden unterschiedliche Ergebnisse liefern, da ich wiederum denke, dass sie "im Prinzip" äquivalent sind, aber vielleicht aus numerischen / rechnerischen Gründen. Oder vielleicht irre ich mich und sie sind selbst im Prinzip nicht gleichwertig.

Jetzt werde ich eine andere Möglichkeit zeigen, die Varianzkomponenten basierend auf traditionellen ANOVA-Ideen zu schätzen. Grundsätzlich nehmen wir die erwarteten mittleren Quadratgleichungen für Ihr Design, ersetzen die beobachteten Werte der mittleren Quadrate und lösen die Varianzkomponenten. Um die erwarteten mittleren Quadrate zu erhalten, verwenden wir eine R-Funktion, die ich vor einigen Jahren geschrieben habe EMS()und die HIER dokumentiert ist . Unten gehe ich davon aus, dass die Funktion bereits geladen ist.

# 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

Okay, jetzt kehren wir zum ursprünglichen Beispiel zurück. Die F-Verhältnisse, die wir zu erreichen versuchen, sind:

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

Hier sind unsere beiden Modelle:

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)

In diesem Fall liefern die beiden Modelle grundsätzlich die gleichen Ergebnisse, obwohl die zweite Methode eine sehr geringfügig höhere Log-Wahrscheinlichkeit aufweist. Keine der beiden Methoden stimmt überein aov(). Aber schauen wir uns an, was wir erhalten, wenn wir wie oben nach den Varianzkomponenten suchen. Verwenden Sie dazu das ANOVA-Verfahren, das die Varianzkomponenten nicht als nicht negativ einschränkt (sondern nur in ausgeglichenen Designs ohne kontinuierliche Prädiktoren und ohne verwendet werden kann) fehlende Daten; die klassischen ANOVA-Annahmen).

# 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

Jetzt können wir sehen, was am ursprünglichen Beispiel pathologisch ist. Das am besten passende Modell impliziert, dass mehrere der zufälligen Varianzkomponenten negativ sind. Aber lmer()(und die meisten anderen gemischten Modellprogramme) beschränken die Schätzungen der Varianzkomponenten auf nicht negativ. Dies wird allgemein als sinnvolle Einschränkung angesehen, da Abweichungen natürlich niemals wirklich negativ sein können. Eine Konsequenz dieser Einschränkung ist jedoch, dass gemischte Modelle Datensätze mit negativen Korrelationen innerhalb der Klasse nicht genau darstellen können, d. H. Datensätze, bei denen die Beobachtungen aus demselben Cluster geringer sind(eher als mehr) im Durchschnitt ähnlich wie Beobachtungen, die zufällig aus dem Datensatz gezogen wurden, und folglich, wenn die Varianz innerhalb des Clusters die Varianz zwischen den Clustern wesentlich übersteigt. Solche Datensätze sind durchaus vernünftige Datensätze, die man gelegentlich in der realen Welt findet (oder versehentlich simuliert!), Sie können jedoch von einem Varianzkomponentenmodell nicht sinnvoll beschrieben werden, da sie negative Varianzkomponenten implizieren. Sie können jedoch von solchen Modellen "unsinnig" beschrieben werden, wenn die Software dies zulässt. aov()erlaubt es. lmer()nicht.

Jake Westfall
quelle
+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- verstehen Sie vielleicht das besser jetzt (zwei Jahre später)? Ich habe versucht herauszufinden, was der Unterschied ist, aber verstehe es auch nicht ...
Amöbe sagt Reinstate Monica
@amoeba Mein derzeitiges Denken ist immer noch ziemlich dasselbe wie damals: AFAIK, die beiden Modelle sind statistisch äquivalent (in dem Sinne, dass sie dieselben Vorhersagen über die Daten treffen und dieselben Standardfehler implizieren), obwohl die zufälligen Effekte parametrisiert sind anders. Ich denke, die beobachteten Unterschiede - die nur manchmal auftreten - sind nur auf Rechenprobleme zurückzuführen. Insbesondere vermute ich, dass Sie mit Optimierungseinstellungen (wie dem Variieren von Startpunkten oder der Verwendung strengerer Konvergenzkriterien) herumspielen könnten, bis die beiden Modelle genau die gleiche Antwort zurückgegeben haben.
Jake Westfall
Danke für deine Antwort. Ich bin eher nicht überzeugt: Ich habe versucht, an den Optimierungseinstellungen herumzuspielen, und konnte an den Ergebnissen nichts ändern. Mein Eindruck ist, dass beide Modelle gut zusammenlaufen. Ich könnte fragen, dass dies irgendwann eine separate Frage ist.
Amöbe sagt Reinstate Monica
Ak(1|A:sub)(0+A|sub)k- -1k(k- -1)/.2k=2Beide Methoden schätzen einen Parameter, und ich bin mir immer noch nicht ganz sicher, warum sie nicht übereinstimmen.
Amöbe sagt Reinstate Monica
Zurück zu diesem Problem ... Ich habe festgestellt, dass für den Zwei-Faktor-Fall, in dem zwei lmerAufrufe identische anova()Ausgaben erzeugen , die zufälligen Effektvarianzen dennoch sehr unterschiedlich sind: siehe VarCorr(mod1)und VarCorr(mod2). Ich verstehe nicht ganz, warum das passiert; machst du? Für mod3und mod4kann man sehen, dass vier von sieben Varianzen für mod3tatsächlich gleich Null sind (für mod4alle sieben sind sie ungleich Null); Diese "Singularität" mod3ist wahrscheinlich der Grund, warum sich die Anova-Tabellen unterscheiden. Abgesehen davon, wie würden Sie Ihren „bevorzugten Weg“ verwenden , wenn aund bmehr als zwei Ebene haben?
Amöbe sagt Reinstate Monica
1

Ist a, b, cfest oder zufällige Effekte? Wenn sie festgelegt sind, lautet Ihre Syntax einfach

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
quelle
Sie sind feste Effekte. Das von Ihnen angepasste ANOVA-Modell ist jedoch nicht das Modell, das das klassische ANOVA-Modell mit wiederholten Messungen zu sein scheint, siehe z . B. hier . Siehe die Fehlerschichten in Ihrem und meinem Fall.
Henrik
1
Eigentlich ist es falsch, wie sie es machen. Wenn Sie ein vollständig gekreuztes faktorielles Design für wiederholte Messungen (oder ein randomisiertes Blockfaktor-Design) haben, sollten Sie abgesehen von subjectallen Effekten (dh Within) nur einen Fehlerterm erhalten . Siehe Experimentelles Design: Verfahren für Verhaltenswissenschaften (2013) von Kirk, Kapitel 10 (S.458) oder meinen Beitrag hier
Masato Nakazawa
Lassen Sie uns diese Frage für den Moment umgehen und davon ausgehen, dass das von mir angepasste Modell das richtige Modell ist. Wie würden Sie dies mit verwenden lmer? Ich werde trotzdem mein Exemplar von Kirk (nur 2. Auflage) bekommen und sehen, was darin steht.
Henrik
Ich bin gespannt, was Sie von Kirks Kapitel halten. Ich denke die Kapitelnummer in der 2. Aufl. ist anders. In der Zwischenzeit werde ich versuchen, verschiedene lmerModelle zu montieren . Der beste Weg, um die Modellanpassung zu überprüfen, besteht darin, ihre dfs mit zu überprüfen, lmerTestda die KR-Näherung Ihnen exactdfs und damit p-Werte geben soll.
Masato Nakazawa
1
Ich habe die 2. Ausgabe von Kirk, in der ich glaube, dass die relevante Diskussion auf den Seiten 443-449 stattfindet, in denen ein Zwei-Wege-Beispiel (kein Drei-Wege-Beispiel) erörtert wird. Die erwarteten mittleren Quadrate, entweder unter der Annahme einer Additivität von A und B oder nicht, sind auf p angegeben. 447. Unter der Annahme, dass A und B fest sind und Subjekte / Blöcke zufällig sind, können wir aus den erwarteten mittleren Quadraten, die Kirk unter "nichtadditives Modell" aufgelistet hat, sehen, dass Tests von A, B und AB jeweils unterschiedliche Fehlerterme beinhalten, nämlich die relevanten Interaktionen mit Block / Subjekt. Das gleiche Prinzip gilt auch für das vorliegende Drei-Wege-Beispiel.
Jake Westfall