Testujesz założenie normalności dla powtarzanych pomiarów anova? (w R)

Zakładając, że istnieje sens w testowaniu założenia normalności dla anova (patrz 1 i 2 )

Jak można to przetestować w R?

Spodziewałbym się zrobić coś takiego:

## From Venables and Ripley (2002) p.165.
utils::data(npk, package="MASS")
npk.aovE <- aov(yield ~  N*P*K + Error(block), npk)
residuals(npk.aovE)
qqnorm(residuals(npk.aov))

Co nie działa, ponieważ „reszty” nie mają metody (ani nie przewidują, jeśli o to chodzi) w przypadku anova powtarzanych pomiarów.

Co więc należy zrobić w tym przypadku?

Czy reszty można po prostu wyodrębnić z tego samego modelu dopasowania bez terminu Błąd? Nie znam wystarczająco literatury, aby wiedzieć, czy jest to ważne, czy nie, z góry dziękuję za wszelkie sugestie.

Możesz nie otrzymać prostej odpowiedzi, residuals(npk.aovE)ale to nie znaczy, że w tym obiekcie nie ma żadnych pozostałości. Zrób stri zobacz, że na poziomach nadal są pozostałości. Wyobrażam sobie, że najbardziej interesował Cię poziom „Wewnątrz”

> residuals(npk.aovE$Within)
          7           8           9          10          11          12 
 4.68058815  2.84725482  1.56432584 -5.46900749 -1.16900749 -3.90234083 
         13          14          15          16          17          18 
 5.08903669  1.28903669  0.35570336 -3.27762998 -4.19422371  1.80577629 
         19          20          21          22          23          24 
-3.12755705  0.03910962  2.60396981  1.13730314  2.77063648  4.63730314 

Moje własne szkolenie i praktyka nie polegały na stosowaniu testów normalności, zamiast na stosowaniu wykresów QQ i testowaniu równoległym za pomocą niezawodnych metod.

Inną opcją byłoby użycie lmefunkcji nlmepakietu (a następnie przekazanie uzyskanego modelu do anova). Możesz użyć residualsna jego wyjściu.

Myślę, że założenie normalności można ocenić dla każdej z powtarzanych miar przed wykonaniem analizy. Przekształciłem ramkę danych, aby każda kolumna odpowiadała powtarzanej mierze, a następnie wykonam test shapiro dla każdej z tych kolumn.

apply(cast(melt(npk,measure.vars="yield"), ...~N+P+K)[-c(1:2)],2,function(x) shapiro.test(x)$p.value)

Venables i Ripley wyjaśniają, w jaki sposób przeprowadzić diagnostykę resztkową dla projektu z powtarzanymi pomiarami w dalszej części swojej książki (s. 284), w części dotyczącej efektów losowych i mieszanych.

Funkcja reszt (lub reszta) jest implementowana dla wyników aov dla każdej warstwy:

z ich przykładu: oats.aov <- aov(Y ~ N + V + Error(B/V), data=oats, qr=T)

Aby uzyskać dopasowane wartości lub wartości resztkowe:

„Tak fitted(oats.aov[[4]])i resid(oats.aov[[4]])wektory długości 54 reprezentujące dopasowanych wartości i pozostałości z ostatniej warstwy”.

Co ważne, dodają:

„Nie można jednoznacznie powiązać ich z fabułami oryginalnego eksperymentu”.

Do diagnostyki używają projekcji:

plot(fitted(oats.aov[[4]]), studres(oats.aov[[4]]))
abline(h=0, lty=2)
oats.pr <- proj(oats.aov)
qqnorm(oats.pr[[4]][, "Residuals"], ylab = "Stratum 4 residuals")
qqline(oats.pr[[4]][, "Residuals"])

Pokazują również, że model można wykonać za pomocą lme, jak napisał inny użytkownik.

Oto dwie opcje, z aov i z lme (myślę, że druga jest preferowana):

require(MASS) ## for oats data set
require(nlme) ## for lme()
require(multcomp) ## for multiple comparison stuff

Aov.mod <- aov(Y ~ N * V + Error(B/V), data = oats)
the_residuals <- aov.out.pr[[3]][, "Residuals"]

Lme.mod <- lme(Y ~ N * V, random = ~1 | B/V, data = oats)
the_residuals <- residuals(Lme.mod)

Oryginalny przykład przyszedł bez interakcji ( Lme.mod <- lme(Y ~ N * V, random = ~1 | B/V, data = oats)), ale wydaje się, że z nim współpracuje (i daje różne wyniki, więc coś robi).

I to wszystko...

ale dla kompletności:

1 - Podsumowania modelu

summary(Aov.mod)
anova(Lme.mod)

2 - Test Tukeya z powtarzanymi pomiarami anova (3 godziny szukania tego !!).

summary(Lme.mod)
summary(glht(Lme.mod, linfct=mcp(V="Tukey")))

3 - Wykresy normalności i homoscedastyczności

par(mfrow=c(1,2)) #add room for the rotated labels
aov.out.pr <- proj(aov.mod)                                            
#oats$resi <- aov.out.pr[[3]][, "Residuals"]
oats$resi <- residuals(Lme.mod)
qqnorm(oats$resi, main="Normal Q-Q") # A quantile normal plot - good for checking normality
qqline(oats$resi)
boxplot(resi ~ interaction(N,V), main="Homoscedasticity", 
        xlab = "Code Categories", ylab = "Residuals", border = "white", 
        data=oats)
points(resi ~ interaction(N,V), pch = 1, 
       main="Homoscedasticity",  data=oats)