Testen der Normalitätsannahme für wiederholte Messungen anova? (in R)

11

Unter der Annahme, dass es sinnvoll ist, die Normalitätsannahme für Anova zu testen (siehe 1 und 2 )

Wie kann es in R getestet werden?

Ich würde erwarten, etwas zu tun wie:

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

Was nicht funktioniert, da "Residuen" keine Methode (und auch keine Vorhersage) für den Fall wiederholter Anova-Messungen haben.

Was ist also in diesem Fall zu tun?

Können die Residuen einfach ohne den Fehlerterm aus demselben Anpassungsmodell extrahiert werden? Ich bin mit der Literatur nicht vertraut genug, um zu wissen, ob dies gültig ist oder nicht. Vielen Dank im Voraus für jeden Vorschlag.

Tal Galili
quelle

Antworten:

5

Möglicherweise erhalten Sie keine einfache Antwort auf residuals(npk.aovE), dies bedeutet jedoch nicht, dass dieses Objekt keine Residuen enthält. Tun Sie strund sehen Sie, dass es innerhalb der Ebenen noch Residuen gibt. Ich würde mir vorstellen, dass Sie am meisten an der Stufe "Innerhalb" interessiert waren

> 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 

Meine eigene Ausbildung und Praxis bestand nicht darin, Normalitätstests zu verwenden, sondern QQ-Diagramme und parallele Tests mit robusten Methoden zu verwenden.

DWin
quelle
Danke Dwin. Ich frage mich, welche der verschiedenen Residuen untersucht werden sollten (außer der Innerhalb). Prost, Tal
Tal Galili
npk.aovE ist eine Liste von drei Elementen. Die ersten beiden sind Parameterschätzungen, und für sie wird keine Normalität angenommen. Daher erscheint es nicht angebracht, etwas anderes als $ Within zu testen. names(npk.aovE)gibt `[1]" (Intercept) "" block "" Within "`
DWin
@Dwin, könntest du den letzten Ansatz überprüfen, den trev gepostet hat (letzte Antwort)? Es verwendet eine Art Projektion, um die Residuen zu berechnen. Für mich ist es der einfachste Ansatz, einen Graphen mit Homogenität der Varianzen zwischen den Gruppen zu zeichnen.
toto_tico
@Dwin, auch qqplot scheint nur lm zu akzeptieren und nicht aov.
toto_tico
2

Eine andere Möglichkeit wäre, die lmeFunktion des nlmePakets zu verwenden (und dann das erhaltene Modell an zu übergeben anova). Sie können residualsauf seiner Ausgabe verwenden.

nico
quelle
1

Ich denke, dass die Normalitätsannahme für jede der wiederholten Messungen bewertet werden kann, bevor die Analyse durchgeführt wird. Ich würde den Datenrahmen so umformen, dass jede Spalte einer wiederholten Messung entspricht, und dann für jede dieser Spalten einen Shapiro.-Test durchführen.

apply(cast(melt(npk,measure.vars="yield"), ...~N+P+K)[-c(1:2)],2,function(x) shapiro.test(x)$p.value)
George Dontas
quelle
Danke gd047 - die Frage ist, was wir tun, wenn wir ein komplexeres Szenario von aov haben (Ausbeute ~ N P K + Fehler (Block / (N + K)), npk). Würde der von Ihnen vorgeschlagene Test die Arbeit erledigen?
Tal Galili
Würden Sie so freundlich sein, den Unterschied zwischen den Szenarien zu erklären? Leider bin ich mit der Verwendung des Begriffs Fehler im Modell nicht vertraut genug (können Sie übrigens ein gutes Buch dazu vorschlagen?). Was ich gerade vorgeschlagen habe, ist die SPSS-Methode zur Überprüfung der Normalitätsannahme, wie ich sie gelernt habe. Sehen Sie dies als Beispiel goo.gl/p45Bx
George Dontas
Hallo gd047. Danke für den Link. Die Dinge, die ich über diese Modelle weiß, sind alle von hier aus verlinkt: r-statistics.com/2010/04/… Wenn Sie andere Ressourcen kennenlernen möchten, würde ich gerne davon erfahren. Prost, Tal
Tal Galili
1

Venables und Ripley erklären später in ihrem Buch (S. 284) im Abschnitt über zufällige und gemischte Effekte, wie eine Restdiagnose für ein Design mit wiederholten Messungen durchgeführt wird.

Die Residuenfunktion (oder Residuen) wird für die aov-Ergebnisse für jede Schicht implementiert:

aus ihrem Beispiel: oats.aov <- aov(Y ~ N + V + Error(B/V), data=oats, qr=T)

So erhalten Sie die angepassten Werte oder Residuen:

"Somit sind fitted(oats.aov[[4]])und resid(oats.aov[[4]])Vektoren der Länge 54 die angepassten Werte und Residuen der letzten Schicht."

Wichtig ist, dass sie hinzufügen:

"Es ist nicht möglich, sie eindeutig mit den Darstellungen des ursprünglichen Experiments zu verknüpfen."

Für die Diagnose verwenden sie eine Projektion:

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

Sie zeigen auch, dass das Modell mit lme erstellt werden kann, wie ein anderer Benutzer gepostet hat.

trev
quelle
nach dieser sollte es [[3]] und nicht [[4]] sein. Ich habe es getestet und es funktioniert nur für [[3]]
toto_tico
1

Hier sind zwei Optionen, mit aov und mit lme (ich denke, die zweite wird bevorzugt):

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)

Das ursprüngliche Beispiel kam ohne die Interaktion ( Lme.mod <- lme(Y ~ N * V, random = ~1 | B/V, data = oats)), aber es scheint damit zu arbeiten (und unterschiedliche Ergebnisse zu erzielen, also tut es etwas).

Und das ist es...

aber der Vollständigkeit halber:

1 - Die Zusammenfassungen des Modells

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

2 - Der Tukey-Test mit wiederholten Messungen anova (3 Stunden danach !!).

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

3 - Die Normalitäts- und Homoskedastizitätskurven

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)

Geben Sie hier die Bildbeschreibung ein

toto_tico
quelle