Co to jest lme4 :: lmer odpowiednik ANOVA z powtarzanymi pomiarami w trzech kierunkach?

Moje pytanie opiera się na tej odpowiedzi, która wykazała, który lme4::lmermodel odpowiada dwukierunkowej ANOVA z powtarzanymi pomiarami:

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

Moje pytanie dotyczy teraz, jak rozszerzyć to na przypadek trójdrożnej ANOVA:

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

Naturalne rozszerzenie oraz jego wersje nie pasują do wyników 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

Należy zauważyć, że bardzo podobne pytanie zostało zadane wcześniej . Brakowało jednak przykładowych danych (które podano tutaj).

— Henrik
źródło

Czy na pewno nie chcesz, aby Twój model dwupoziomowy był y ~ a*b + (1 + a*b|subject), d[d$c == "1",]? A może coś mi brakuje?

— Rasmus Bååth

@ RasmusBååth Go Ahead i spróbuj go dopasować, lmernarzeka, ponieważ losowe efekty nie są już identyfikowane. Początkowo myślałem również, że to model, który chcę, ale tak nie jest. Jeśli porównasz lmer model, który proponuję dla skrzynki 2-drożnej ze standardową ANOVA, zobaczysz, że wartości F dokładnie pasują. Jak powiedziałem w odpowiedzi, połączyłem.

— Henrik

W przypadku problemu trójdrożnego nie oczekuje się , że pierwszy lmernapisany przez ciebie model (który wyklucza losowe interakcje dwukierunkowe) będzie równoważny 3-kierunkowej RM-ANOVA, ale drugi napisany przez ciebie model (który obejmuje losowy dwustronne interakcje) powinny być. Jeśli chodzi o to, dlaczego istnieje rozbieżność nawet w przypadku tego modelu, mam przeczucie, na czym polega problem, pójdę na obiad, a następnie przejrzę zestaw danych z zabawkami.

— Jake Westfall

Bezpośrednią odpowiedzią na twoje pytanie jest to, że ostatni model, który napisałeś,

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

Uważam, że jest „w zasadzie” poprawny, chociaż jest to dziwna parametryzacja, która nie zawsze wydaje się dobrze sprawdzać w praktyce.

Jeśli chodzi o to, dlaczego dane wyjściowe z tego modelu są niezgodne z danymi aov()wyjściowymi, myślę, że istnieją dwa powody.

Twój prosty symulowany zestaw danych jest patologiczny, ponieważ najlepiej dopasowany model to taki, który sugeruje komponenty ujemnej wariancji, na które lmer()nie pozwalają ze sobą modele mieszane (i większość innych programów modeli mieszanych).
Nawet w przypadku niepatologicznego zestawu danych sposób, w jaki skonfigurowano model, jak wspomniano powyżej, nie zawsze wydaje się dobrze działać w praktyce, chociaż muszę przyznać, że tak naprawdę nie rozumiem, dlaczego. Moim zdaniem jest to również po prostu dziwne, ale to już inna historia.

Pozwól mi najpierw zademonstrować preferowaną parametryzację na twoim początkowym dwukierunkowym przykładzie ANOVA. Załóżmy, że Twój zestaw danych djest załadowany. Twój model (zwróć uwagę, że zmieniłem kody zastępcze na kontrastowe) to:

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

które działało tutaj dobrze, ponieważ pasowało do aov()wyniku. Model, który preferuję, obejmuje dwie zmiany: ręczne kodowanie kontrastu czynników, abyśmy nie pracowali z obiektami współczynnika R (co zalecam robić w 100% przypadków) i inne określenie efektów losowych:

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)

Te dwa podejścia są całkowicie równoważne w prostym problemie dwukierunkowym. Teraz przejdziemy do problemu trójstronnego. Wspomniałem wcześniej, że podany przykładowy zestaw danych był patologiczny. Tak więc, co chcę zrobić przed zaadresowaniem twojego przykładowego zestawu danych, to najpierw wygenerować zestaw danych z rzeczywistego modelu komponentów wariancji (tj. Gdzie niezerowe komponenty wariancji są wbudowane w prawdziwy model). Najpierw pokażę, w jaki sposób moja preferowana parametryzacja wydaje się działać lepiej niż ta, którą zaproponowałeś. Wtedy będę wykazać w inny sposób szacowania składników wariancji, która ma nie narzucać, że muszą być nieujemne. Wtedy będziemy w stanie zobaczyć problem z oryginalnym przykładowym zestawem danych.

Nowy zestaw danych będzie identyczny pod względem struktury, z tym wyjątkiem, że będziemy mieli 50 tematów:

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

Współczynniki F, które chcemy dopasować, to:

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

Oto nasze dwa modele:

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)

Jak widzimy, tylko druga metoda pasuje do wyniku aov(), chociaż pierwsza metoda jest co najmniej na boisku. Druga metoda zapewnia również wyższe prawdopodobieństwo logarytmiczne. Nie jestem pewien, dlaczego te dwie metody dają różne wyniki, ponieważ znów uważam, że są one „w zasadzie” równoważne, ale być może dzieje się tak z kilku powodów liczbowych / obliczeniowych. A może się mylę i nawet w zasadzie nie są one równoważne.

Teraz pokażę inny sposób szacowania składników wariancji w oparciu o tradycyjne pomysły ANOVA. Zasadniczo weźmiemy oczekiwane równania średnich kwadratów dla twojego projektu, podstawimy obserwowane wartości średnich kwadratów i rozwiążemy składowe wariancji. Aby uzyskać oczekiwane średnie kwadraty, użyjemy funkcji R, którą napisałem kilka lat temu, zwanej EMS(), co jest udokumentowane TUTAJ . Poniżej zakładam, że funkcja jest już załadowana.

# 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

OK, teraz wrócimy do oryginalnego przykładu. Współczynniki F, które próbujemy dopasować, to:

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

Oto nasze dwa modele:

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)

W tym przypadku dwa modele dają w zasadzie te same wyniki, chociaż druga metoda ma bardzo nieznacznie większe prawdopodobieństwo logarytmiczne. Żadna z metod nie pasuje aov(). Ale spójrzmy na to, co otrzymujemy, gdy rozwiązujemy dla składników wariancji, jak to zrobiliśmy powyżej, stosując procedurę ANOVA, która nie ogranicza składników wariancji, aby były nieujemne (ale które można stosować tylko w zrównoważonych projektach bez ciągłych predyktorów i bez brakujące dane; klasyczne założenia 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

Teraz możemy zobaczyć, co jest patologiczne w oryginalnym przykładzie. Model najlepiej dopasowany to taki, który sugeruje, że kilka losowych składników wariancji jest ujemnych. Ale lmer()(i większość innych programów modeli mieszanych) ogranicza szacunki składników wariancji do wartości nieujemnych. Jest to ogólnie uważane za rozsądne ograniczenie, ponieważ odchylenia oczywiście nigdy naprawdę nie mogą być ujemne. Jednak konsekwencją tego ograniczenia jest to, że modele mieszane nie są w stanie dokładnie przedstawić zestawów danych, które wykazują ujemne korelacje wewnątrzklasowe, to znaczy zestawów danych, w których obserwacje z tego samego skupienia są mniejsze(raczej niż bardziej) podobny średnio niż obserwacje losowe z zestawu danych, aw konsekwencji, gdy wariancja wewnątrz klastra znacznie przekracza wariancję między klastrami. Takie zestawy danych są całkowicie uzasadnionymi zestawami danych, które czasami można spotkać w świecie rzeczywistym (lub przypadkowo symulować!), Ale nie można ich rozsądnie opisać w modelu wariancji-składników, ponieważ implikują negatywne składniki wariancji. Można je jednak opisać „nie rozsądnie” przez takie modele, jeśli oprogramowanie na to pozwoli. aov()pozwala na to. lmer()nie.

— Jake Westfall
źródło

+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

- czy może to zrozumieć lepiej (dwa lata później)? Próbowałem dowiedzieć się, jaka jest różnica, ale też nie rozumiem ...

— ameba

@amoeba Moje obecne myślenie jest prawie takie samo jak wtedy: AFAIK, oba modele są statystycznie równoważne (w tym sensie, że przewidują te same dane i sugerują te same błędy standardowe), nawet jeśli losowe efekty są sparametryzowane różnie. Myślę, że zaobserwowane różnice - które zdają się występować tylko czasami - wynikają wyłącznie z problemów obliczeniowych. W szczególności podejrzewam, że można było bawić się przy ustawieniach optymalizatora (takich jak zmienianie punktów początkowych lub stosowanie bardziej rygorystycznych kryteriów konwergencji), dopóki oba modele nie zwrócą dokładnie tej samej odpowiedzi.

— Jake Westfall,

Dzięki za odpowiedź. Jestem raczej nieprzekonany: próbowałem majstrować przy ustawieniach optymalizatora i nie mogłem nic zmienić w wynikach; mam wrażenie, że oba modele są dobrze połączone. W pewnym momencie mogę zadać to osobne pytanie.

— ameba

Ciągle badam ten problem i zdaję sobie sprawę z tego, że: gdy tylko czynniki powtarzanego pomiaru mają więcej niż dwa poziomy, te dwie metody stają się zupełnie inne! Jeśli Ama poziomów następnie ocenia tylko jeden parametr wariancji, natomiast szacunki parametrów wariancji i korelacji między nimi. Rozumiem, że klasyczna ANOVA szacuje tylko jeden parametr wariancji, a zatem powinna być równoważna pierwszej metodzie, a nie drugiej. Dobrze? Ale dla obie metody szacują jeden parametr i nadal nie jestem pewien, dlaczego się nie zgadzają.

k

$k$ (1|A:sub)(0+A|sub)

k - 1

$k-1$

k (k - 1) / 2

$k(k-1)/2$

k = 2

$k=2$

— ameba

Wracając do tego problemu ... Zauważyłem, że dla przypadku dwuskładnikowego, w którym dwa lmerwywołania dają identyczne anova()wyjście, wariancje efektów losowych są jednak całkiem różne: patrz VarCorr(mod1)i VarCorr(mod2). Nie do końca rozumiem, dlaczego tak się dzieje; czy ty? Dla mod3i mod4widać, że cztery z siedmiu wariancji dla mod3są w rzeczywistości równe zeru (dla mod4wszystkich siedmiu są niezerowe); ta „osobliwość” mod3jest prawdopodobnie powodem, dla którego tabele anova różnią się. Poza tym, jak użyłbyś swojego „preferowanego sposobu”, gdyby ai bmiał więcej niż dwa poziomy?

— ameba

Są a, b, cstałe lub losowe efekty? Jeśli zostaną naprawione, twoja składnia będzie po prostu

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
źródło

Są to ustalone efekty. Jednak dopasowany model ANOVA nie jest modelem, który wydaje się być klasycznym modelem ANOVA z powtarzanymi pomiarami, patrz np . Tutaj . Zobacz warstwy błędów w twoim i moim przypadku.

— Henrik

W rzeczywistości sposób, w jaki to robią, jest nieprawidłowy. Jeśli masz w pełni przekreślony czynnikowy układ powtarzanych miar (lub losowy blokowy układ czynnikowy), powinieneś otrzymać tylko 1 warunek błędu, oprócz subject, dla wszystkich efektów (tj Within.). Zobacz Projekt eksperymentalny: Procedury dla nauk behawioralnych (2013) autorstwa Kirka, rozdział 10 (str. 458) lub mój post tutaj

— Masato Nakazawa

Pomińmy na chwilę to pytanie i załóżmy, że model, który zamontowałem, byłby właściwym modelem. Jak byś do tego pasował lmer? Mimo to dostanę moją kopię Kirka (tylko druga edycja) i zobaczę, co mówi.

— Henrik

Ciekawe, co myślisz o rozdziale Kirka. Myślę, że numer rozdziału w 2 wyd. jest inny. Tymczasem postaram się dopasować do różnych lmermodeli. Najlepszym sposobem sprawdzenia dopasowania modelu jest sprawdzenie ich dfs, lmerTestponieważ przybliżenie KR ma dać ci exactdfs, a zatem wartości p.

— Masato Nakazawa

Mam drugie wydanie Kirka, w którym moim zdaniem odpowiednia dyskusja znajduje się na str. 443–449, która omawia przykład dwukierunkowy (a nie trójstronny). Oczekiwane średnie kwadraty, przy założeniu addytywności A i B, lub nie, podano na p. 447. Zakładając, że A i B są stałe, a podmioty / bloki są losowe, z oczekiwanych średnich kwadratów wymienionych przez Kirka w „modelu nieaddytywnym” możemy zobaczyć, że testy A, B i AB obejmują różne terminy błędu, a mianowicie: odpowiednie interakcje z blokiem / tematem. Ta sama zasada rozciąga się na obecny przykład trójstronny.

— Jake Westfall