ANOVA Split-Plot: testy porównawcze modeli w R.

Jak mogę testować efekty w analizie ANOVA typu split-plot, używając odpowiednich porównań modeli do użycia z argumentami Xi w R? Znam i Dalgaard (2007) [1]. Niestety szczotkuje tylko projekty Split-Plot. Robi to w całkowicie losowy sposób z dwoma czynnikami wewnątrz badanych:Manova.mlm()?anova.mlm

N  <- 20  # 20 subjects total
P  <- 3   # levels within-factor 1
Q  <- 3   # levels within-factor 2
DV <- matrix(rnorm(N* P*Q), ncol=P*Q)           # random data in wide format
id <- expand.grid(IVw1=gl(P, 1), IVw2=gl(Q, 1)) # intra-subjects layout of data matrix

library(car)        # for Anova()
fitA <- lm(DV ~ 1)  # between-subjects design: here no between factor
resA <- Anova(fitA, idata=id, idesign=~IVw1*IVw2)
summary(resA, multivariate=FALSE, univariate=TRUE)  # all tests ...

Poniższe porównania modeli prowadzą do tych samych wyników. Model ograniczony nie obejmuje danego efektu, ale wszystkie inne efekty tego samego lub niższego rzędu, pełny model dodaje dany efekt.

anova(fitA, idata=id, M=~IVw1 + IVw2, X=~IVw2, test="Spherical") # IVw1
anova(fitA, idata=id, M=~IVw1 + IVw2, X=~IVw1, test="Spherical") # IVw2
anova(fitA, idata=id, M=~IVw1 + IVw2 + IVw1:IVw2,
                      X=~IVw1 + IVw2, test="Spherical")          # IVw1:IVw2

Projekt Split-Splot z jednym czynnikiem wewnątrz i jednym między podmiotami:

idB  <- subset(id, IVw2==1, select="IVw1")          # use only first within factor
IVb  <- gl(2, 10, labels=c("A", "B"))               # between-subjects factor
fitB <- lm(DV[ , 1:P] ~ IVb)                        # between-subjects design
resB <- Anova(fitB, idata=idB, idesign=~IVw1)
summary(resB, multivariate=FALSE, univariate=TRUE)  # all tests ...

Są to anova()polecenia do replikacji testów, ale nie wiem, dlaczego działają. Dlaczego testy następujących porównań modeli prowadzą do takich samych wyników?

anova(fitB, idata=idB, X=~1, test="Spherical") # IVw1, IVw1:IVb
anova(fitB, idata=idB, M=~1, test="Spherical") # IVb

Dwa czynniki wewnętrzne i jeden czynnik wewnętrzny:

fitC <- lm(DV ~ IVb)  # between-subjects design
resC <- Anova(fitC, idata=id, idesign=~IVw1*IVw2)
summary(resC, multivariate=FALSE, univariate=TRUE)  # all tests ...

Jak zreplikować wyniki podane powyżej z odpowiednimi porównaniami modeli do użycia z argumentami Xi ? Jaka jest logika tych porównań modeli?Manova.mlm()

EDYCJA: suncoolsu wskazało, że dla wszystkich praktycznych celów dane z tych projektów powinny być analizowane przy użyciu modeli mieszanych. Jednak nadal chciałbym zrozumieć, jak powielić wyniki summary(Anova())z anova.mlm(..., X=?, M=?).

[1]: Dalgaard, P. 2007. Nowe funkcje analizy wielowymiarowej. R News, 7 (2), 2-7.

XI Mw zasadzie określić dwa modele, które chcesz porównać, ale tylko pod względem wewnątrzlaboratoryjnych temat skutków; następnie pokazuje wyniki interakcji wszystkich efektów między podmiotami (w tym przechwytywania) z efektami wewnątrz podmiotu, które zmieniły się między Xi M.

Twoje przykłady fitBsą łatwiejsze do zrozumienia, jeśli dodamy wartości domyślne Xi M:

anova(fitB, idata=idB, M=~1, X=~0, test="Spherical") # IVb
anova(fitB, idata=idB, M=diag(3), X=~1, test="Spherical") # IVw1, IVw1:IVb

Pierwszym modelem jest zmiana z braku efektów wewnątrz podmiotu (wszystkie mają tę samą średnią) na inną średnią dla każdego, dlatego dodaliśmy idefekt losowy, który jest właściwy do przetestowania ogólnego przechwytywania i ogólnego efektu między podmiotami na.

Drugi model reklamuje id:IVw1interakcję, którą należy przetestować, IVw1a także IVw1:IVbwarunki. Ponieważ istnieje tylko jeden efekt wewnątrz podmiotu (z trzema poziomami), uwzględniany jest domyślny diag(3)w drugim modelu; byłoby to równoważne z uruchomieniem

anova(fitB, idata=idB, M=~IVw1, X=~1, test="Spherical") # IVw1, IVw1:IVb

Wydaje mi się fitC, że te polecenia odtworzą Anovapodsumowanie.

anova(fitC, idata=id, M=~1, X=~0, test="Spherical") #IVb
anova(fitC, idata=id, M=~IVw1 + IVw2, X=~IVw2, test="Spherical") # IVw1
anova(fitC, idata=id, M=~IVw1 + IVw2, X=~IVw1, test="Spherical") # IVw2
anova(fitC, idata=id, M=~IVw1 + IVw2 + IVw1:IVw2,
                  X=~IVw1 + IVw2, test="Spherical")          # IVw1:IVw2

Teraz, jak odkryłeś, te polecenia są naprawdę trudne. Na szczęście nie ma już powodu, aby z nich korzystać. Jeśli chcesz założyć sferyczność, powinieneś po prostu użyć aov, lub dla jeszcze łatwiejszej składni, po prostu użyć lmi obliczyć odpowiednie testy F. Jeśli nie chcesz zakładać sferyczności, używanie lmejest naprawdę dobrym rozwiązaniem, ponieważ masz znacznie większą elastyczność niż w przypadku poprawek GG i HF.

Na przykład, oto kod aovi lmdla twojego fitA. Najpierw musisz mieć dane w długim formacie; Oto jeden ze sposobów, aby to zrobić:

library(reshape)
d0 <- data.frame(id=1:nrow(DV), DV)
d0$IVb <- IVb
d0 <- melt(d0, id.vars=c(1,11), measure.vars=2:10)
id0 <- id
id0$variable <- factor(levels(d0$variable), levels=levels(d0$variable))
d <- merge(d0, id0)
d$id <- factor(d$id)

A oto lm andkod aov:

anova(lm(value ~ IVw1*IVw2*id, data=d))
summary(aov(value ~ IVw1*IVw2 + Error(id/(IVw1*IVw2)), data=d))

Projekty podzielonej działki powstały w rolnictwie, stąd nazwa. Ale często się zdarzają i powiedziałbym - koń większości badań klinicznych. Główny wykres traktuje się poziomem jednego czynnika, podczas gdy poziomy niektórych innych czynników mogą się różnić w zależności od wykresów podrzędnych. Projekt powstaje w wyniku ograniczeń pełnej randomizacji. Na przykład: pole można podzielić na cztery wykresy cząstkowe. Możliwe jest sadzenie różnych odmian w podplotach, ale na całym polu można zastosować tylko jeden rodzaj nawadniania. Nie rozróżnienie między podziałami i blokami. Bloki są cechami jednostek eksperymentalnych, z których mamy możliwość skorzystać w projekcie eksperymentalnym, ponieważ wiemy, że tam są. Z drugiej strony podziały nakładają ograniczenia na możliwe przypisania czynników. Nakładają na projekt wymagania, które uniemożliwiają całkowitą randomizację.

Są one często stosowane w badaniach klinicznych, w których jeden czynnik można łatwo zmienić, podczas gdy inny czynnik wymaga znacznie więcej czasu. Jeśli eksperymentator musi wykonać wszystkie przebiegi dla każdego poziomu trudnego do zmiany współczynnika kolejno, powstaje projekt podziału wykresu, przy czym trudny do zmiany współczynnik reprezentuje cały współczynnik wykresu.

Oto przykład: W próbie rolniczej celem było określenie skutków dwóch odmian upraw i czterech różnych metod nawadniania. Osiem pól było dostępnych, ale dla każdego pola można zastosować tylko jeden rodzaj nawadniania. Pola mogą być podzielone na dwie części o różnych odmianach w każdej części. Cały współczynnik powierzchni stanowi nawadnianie, które należy losowo przypisać do pól. W obrębie każdego pola odmiana jest przypisana.

Oto jak to zrobić w R:

install.packages("faraway")
data(irrigation)
summary(irrigation)

library(lme4)

R> (lmer(yield ~ irrigation * variety + (1|field), data = irrigation))
Linear mixed model fit by REML 
Formula: yield ~ irrigation * variety + (1 | field) 
   Data: irrigation 
  AIC  BIC logLik deviance REMLdev
 65.4 73.1  -22.7     68.6    45.4
Random effects:
 Groups   Name        Variance Std.Dev.
 field    (Intercept) 16.20    4.02    
 Residual              2.11    1.45    
Number of obs: 16, groups: field, 8

Fixed effects:
                       Estimate Std. Error t value
(Intercept)               38.50       3.02   12.73
irrigationi2               1.20       4.28    0.28
irrigationi3               0.70       4.28    0.16
irrigationi4               3.50       4.28    0.82
varietyv2                  0.60       1.45    0.41
irrigationi2:varietyv2    -0.40       2.05   -0.19
irrigationi3:varietyv2    -0.20       2.05   -0.10
irrigationi4:varietyv2     1.20       2.05    0.58

Correlation of Fixed Effects:
            (Intr) irrgt2 irrgt3 irrgt4 vrtyv2 irr2:2 irr3:2
irrigation2 -0.707                                          
irrigation3 -0.707  0.500                                   
irrigation4 -0.707  0.500  0.500                            
varietyv2   -0.240  0.170  0.170  0.170                     
irrgtn2:vr2  0.170 -0.240 -0.120 -0.120 -0.707              
irrgtn3:vr2  0.170 -0.120 -0.240 -0.120 -0.707  0.500       
irrgtn4:vr2  0.170 -0.120 -0.120 -0.240 -0.707  0.500  0.500

Zasadniczo, co mówi ten model, nawadnianie i różnorodność są ustalonymi efektami, a różnorodność jest zagnieżdżona w nawadnianiu. Pola są efektami losowymi i obrazowo będą podobne

I_1 | I_2 | I_3 | I_4

V_1 V_2 | V_1 V_2 | V_1 V_2 | V_1 V_2

Był to jednak specjalny wariant ze stałym efektem całej fabuły i efektem podplotu. Mogą istnieć warianty, w których jeden lub więcej jest losowych. Mogą być bardziej skomplikowane projekty, takie jak split-split .. projekty fabuły. Zasadniczo możesz oszaleć i zwariować. Ale biorąc pod uwagę, że podstawowa struktura i rozkład (tj. Stały lub losowy, zagnieżdżony lub skrzyżowany, ...) jest jasno zrozumiany, model nie lmer-Ninjabędzie miał problemów z modelowaniem. Być może interpretacja będzie bałaganem.

Jeśli chodzi o porównania, powiedz, że masz lmer1i lmer2:

anova(lmer1, lmer2)

da ci odpowiedni test oparty na statystyce testu chi-sq ze stopniami swobody równymi różnicy parametrów.

por .: Faraway, J., Rozszerzanie modeli liniowych o R.

Casella, G., Projekt statystyczny