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

11

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.

Tal Galili
źródło

Odpowiedzi:

5

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.

DWin
źródło
Dziękuję Dwin. Zastanawiam się, która z różnych pozostałości powinna zostać zbadana (oprócz Wewnątrz). Pozdrawiam, Tal
Tal Galili
npk.aovE to lista trzech elementów. Pierwsze dwa są szacunkami parametrów i nie zakłada się dla nich normalności, więc testowanie czegokolwiek poza $ Within nie wydaje się właściwe. names(npk.aovE)zwraca „[1]” (przechwycenie) „” blok ”„ Wewnątrz ”„
DW
@Dwin, czy mógłbyś sprawdzić ostatnie podejście, które napisał trev (ostatnia odpowiedź)? Wykorzystuje rodzaj projekcji do obliczenia reszt. Jest to dla mnie najłatwiejsze podejście do wykreślenia wykresu z jednorodnością różnic między grupami.
toto_tico,
@Dwin, również qqplot wydaje się akceptować tylko lm, a nie aov.
toto_tico,
2

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

Nico
źródło
1

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)
George Dontas
źródło
Dzięki gd047 - pytanie brzmi: co robimy, gdy mamy bardziej złożony scenariusz aov (wydajność ~ N P K + Błąd (blok / (N + K)), npk) czy test, który zaproponujesz, zadziała?
Tal Galili
Czy byłbyś na tyle uprzejmy, aby wyjaśnić różnicę między scenariuszami? Niestety, nie jestem wystarczająco zaznajomiony z użyciem terminu Błąd w modelu (a propos, czy możesz zasugerować dobrą książkę na ten temat?). Właśnie zaproponowałem SPSS sposób sprawdzenia założenia normalności, tak jak się tego nauczyłem. Zobacz to jako przykład goo.gl/p45Bx
George Dontas
Cześć gd047. Dziękuję za link. Wszystkie rzeczy, które wiem o tych modelach, są tutaj powiązane: r-statistics.com/2010/04/… Jeśli poznasz inne zasoby - chciałbym o nich wiedzieć. Pozdrawiam, Tal
Tal Galili
1

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.

trev
źródło
zgodnie z tym powinien to być [[3]], a nie [[4]]. Przetestowałem to i działa tylko dla [[3]]
toto_tico,
1

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)

wprowadź opis zdjęcia tutaj

toto_tico
źródło