Plotten zur Überprüfung der Homoskedastizitätsannahme für ANOVA mit wiederholten Messungen in R.

Ich habe mit der aov()Funktion eine ANOVA mit wiederholten Messungen innerhalb der Probanden durchgeführt . Meine abhängige Variable ist nicht normal verteilt, daher bin ich sehr daran interessiert, Annahme-Tests für meine Analyse durchzuführen. Es scheint, dass das bloße Aufrufen plot()der Ausgabe bei wiederholten Messungen nicht funktioniert. Daher habe ich die Residuen und die angepassten Werte für ein interessierendes Modell manuell genommen und gegeneinander aufgetragen. Ich gehe davon aus, dass ich auf diese Weise planen würde, die Annahme der Homoskedastizität zu testen.

Das Diagramm enthält 2 vertikale Bänder (siehe Abbildung unten). Es stellt sich heraus, dass die angepassten Werte alle um 2 Werte zentriert sind (obwohl ==sie nicht genau gleich sind), wobei einer der negative des anderen ist.

Ich habe 2 Fragen:

1) Ist dies der richtige Weg, um die Annahme der Homoskedastizität manuell zu testen? Wenn nicht, wie würde ich bei Designs mit wiederholten Messungen vorgehen (da nur das Aufrufen plot()nicht funktioniert)?

2) Wenn es richtig ist, was sagt mir diese Handlung? Warum sind die angepassten Werte so gruppiert? Was kann ich daraus schließen?

Vielen Dank für jede Eingabe hier. Wenn Sie bessere Möglichkeiten kennen, um Annahmen in rm-ANOVAs zu überprüfen (vorzugsweise zu zeichnen), wären dies ebenfalls nützliche Informationen.

Ich habe hier einige Scheindaten eingefügt, um das Szenario zu replizieren:

#Create mock data (there's probably a more efficient way to do this.. would also be nice to know! :) )
p <- sort(rep(1:20,8))
y <- rep(rep(1:2,4),20)
z <- rep(rep(c(1,1,2,2),2),20)
w <- rep(c(1,1,1,1,2,2,2,2),20)
x <- rnorm(160,10,2)

d <- data.frame(x,p=factor(p),y=factor(y),z=factor(z),w=factor(w))

#Run repeated-measures ANOVA
ex.aov <- aov(x ~ y*z*w + Error(p/(y*z*w)), d)

#Try to plot full object (doesn't work)
plot(ex.aov)

#Try to plot section of object (doesn't work)
plot(ex.aov[["p:y:z"]])

#Plot residuals against fitted (custom "skedasticity" plot - works)
plot(residuals(ex.aov[["p:y:z"]])~fitted(ex.aov[["p:y:z"]]))

Beginnen Sie mit der Bearbeitung

In Anbetracht der von @Stefan bereitgestellten Informationen habe ich nachfolgend einige zusätzliche Details hinzugefügt, wobei ich die von ihm vorgeschlagene verbesserte Datenstruktur verwendet habe:

# Set seed to make it reproducible
set.seed(12)

#New variable names and generation
subj <- sort(factor(rep(1:20,8)))
x1 <- rep(c('A','B'),80)
x2 <- rep(c('A','B'),20,each=2)
x3 <- rep(c('A','B'),10, each=4)
outcome <- rnorm(80,10,2)

d3 <- data.frame(outcome,subj,x1,x2,x3)

#Repeated measures ANOVA
ex.aov <- aov(outcome ~ x1*x2*x3 + Error(subj/(x1*x2*x3)), d3)

#proj function
ex.aov.proj <- proj(ex.aov)

# Check for normality by using last error stratum
qqnorm(ex.aov.proj[[9]][, "Residuals"])
# Check for heteroscedasticity by using last error stratum
plot(ex.aov.proj[[9]][, "Residuals"])

Die resultierenden Diagramme sind unten:

Kann jemand die Bilder oben interpretieren (besonders das letzte)? Es sieht so aus, als ob es Clustering und Musterstruktur gibt. Kann es verwendet werden, um auf das Vorhandensein von Heteroskedastizität zu schließen?

Ich gehe davon aus, dass ein Modell, das mit der Error()Funktion innerhalb angepasst wurde, aov()bei Verwendung in nicht funktioniert, plot()da Sie mehr als eine Fehlerschicht erhalten, aus der Sie auswählen können. Nach diesen Informationen hier sollte man nun die proj()Funktion verwenden, die Ihnen die Residuen für jede Fehlerschicht gibt, die dann für diagnostische Diagramme verwendet werden können.

Bearbeiten Sie 1 Start

Weitere Informationen zu Multistratum-Modellen und zur proj()Funktion finden Sie in Venables und Ripley, Seite 284 (jedoch ab Seite 281): Residuen in Multistratum-Analysen: Projektionen . Im zweiten Satz schreiben sie (ich habe fett hervorgehoben):

Somit sind fitted(oats.aov[[4]])und resid(oats.aov[[4]])sind Vektoren der Länge 54, die angepasste Werte und Residuen aus der letzten Schicht darstellen , basierend auf 54 orthonormalen linearen Funktionen des ursprünglichen Datenvektors. Es ist nicht möglich, sie eindeutig mit den Darstellungen des ursprünglichen Experiments zu verknüpfen. Die Funktion projnimmt ein angepasstes Modellobjekt und findet die Projektionen des ursprünglichen Datenvektors auf die Teilräume, die durch jede Zeile in der Analyse der Varianztabellen definiert sind (einschließlich für Multistratum-Objekte die unterdrückte Tabelle nur mit dem großen Mittelwert). Das Ergebnis ist eine Liste von Matrizen, eine für jede Schicht, wobei die Spaltennamen für jede die Komponentennamen aus der Analyse der Varianztabellen sind.

Für Ihr Beispiel bedeutet das:

ex.aov.proj <- proj(ex.aov)

# Check number of strata 
summary(ex.aov.proj)

# Check for normality by using last error stratum
qqnorm(ex.aov.proj[[9]][, "Residuals"])
# Check for heteroscedasticity by using last error stratum
plot(ex.aov.proj[[9]][, "Residuals"])

Dies wird jedoch auch zu Handlungen führen, die ich nicht vollständig interpretieren kann (insbesondere die zweite).

In ihrem Fall war die letzte Schicht die WithinSchicht. Da Ihr Modell dies nicht abschätzen kann (vermutlich aufgrund Ihres Fehlerausdrucks), bin ich mir nicht sicher, ob die einfache Verwendung Ihrer letzten Schicht gültig ist.

Hoffentlich kann jemand anderes das klären.

1 Ende bearbeiten

Bearbeiten 2 starten

Nach dieser Quelle sollten Residuen zur Beurteilung der Normalität und Heteroskedastizität ohne die Error()Funktion überprüft werden .

Um Annahmen zu überprüfen, müssen Sie den Fehlerbegriff nicht verwenden. Sie können den Begriff fehlerfrei hinzufügen, aber die F-Tests sind falsch. Die Annahmeprüfung ist jedoch in Ordnung.

Das scheint mir vernünftig, aber ich hoffe, jemand anderes könnte es klarstellen.

Bearbeiten Sie 2 Ende

Mein alternativer Vorschlag:

Zuerst habe ich Ihren Datensatz leicht geändert und einen Startwert festgelegt, um ihn reproduzierbar zu machen (möglicherweise hilfreich für einige Probleme, die Sie in Zukunft haben):

# Set seed to make it reproducible
set.seed(12)

# I changed the names of your variables to make them easier to remember
# I also deleted a few nested `rep()` commands. Have a look at the `each=` argument.
subj <- sort(factor(rep(1:20,8)))
x1 <- rep(c('A','B'),80)
x2 <- rep(c('A','B'),20,each=2)
x3 <- rep(c('A','B'),10, each=4)
outcome <- rnorm(80,10,2)

d3 <- data.frame(outcome,subj,x1,x2,x3)

Zweitens habe ich stattdessen ein lineares Mischeffektmodell verwendet, da Sie wiederholte Messungen und damit einen zufälligen Begriff haben, den Sie verwenden können:

require(lme4)
# I specified `subj` as random term to account for the repeated measurements on subject.
m.lmer<-lmer(outcome ~ x1*x2*x3 + (1|subj), data = d3)
summary(m.lmer)

# Check for heteroscedasticity
plot(m.lmer)

# or
boxplot(residuals(m.lmer) ~ d3$x1 + d3$x2 + d3$x3)

# Check for normality
qqnorm(residuals(m.lmer))

Mit dem afexPaket können Sie auch die festen Effekte im ANOVA-Tabellenformat erhalten (Sie können die Anova()Funktion aus dem carPaket auch als weitere Option verwenden):

require(afex)
mixed(outcome ~ x1*x2*x3 + (1|subj), data = d3, method="LRT")

Fitting 8 (g)lmer() models:
[........]
    Effect df    Chisq p.value
1       x1  1     0.04     .84
2       x2  1     2.53     .11
3       x3  1  7.68 **    .006
4    x1:x2  1  8.34 **    .004
5    x1:x3  1 10.51 **    .001
6    x2:x3  1     0.31     .58
7 x1:x2:x3  1     0.12     .73

Überprüfen Sie ?mixeddie verschiedenen Optionen, die Sie auswählen können. Auch in Bezug auf gemischte Modelle gibt es hier viele Informationen zu Cross Validated.

Vollständiger Haftungsausschluss: Ich liebe es, R für viele verschiedene Analysen zu verwenden, aber ich mache keine ANOVAs in R.

Frage 1 : Im analytischen Kontext von ANOVAs bin ich eher damit vertraut, diese Annahme durch Tests der Homogenität von Varianzen zu bewerten, als Homo / Heteroskedastizität zu zeichnen und visuell zu bewerten. Obwohl es mehrere Tests zur Homogenität der Varianz gibt, sehe ich am meisten den Levene-Test. In R scheint es, dass Sie dies über das carPaket mit der leveneTestFunktion tun können.

Basierend auf Ihren Daten würde es so aussehen : leveneTest(x ~ y*z*w, d). Beachten Sie, dass ich nicht glaube, dass Sie die Fehlerstruktur für wiederholte Messungen in dieser Funktion angeben können, und ehrlich gesagt bin ich mir nicht sicher, ob / inwieweit dies für den Levene-Test von Bedeutung ist. Im Vergleich zu anderen statischen Analysesoftware scheint es eine gewisse Variabilität hinsichtlich der Durchführung des Levene-Tests bei ANOVA mit wiederholten Messungen zu geben. SPSS bietet beispielsweise separate Levene-Tests zwischen Gruppen für jede Ebene Ihrer wiederholten Messung, während die leveneTestFunktion einen umfassenden Test aller Ebenen aller Variablen bietet - andere Software bietet möglicherweise auch alternative Ansätze. Auf jeden Fall scheint der SPSS-Ansatz auch die Abhängigkeit der Daten zu ignorieren, indem nur die Homogenität der Varianz zwischen den Gruppen bewertet wird.

Frage 2 : Wenn Sie einen Test der Homogenität der Varianz verwenden möchten - Levene oder auf andere Weise -, wäre es wahrscheinlich informativer, einfache Balkendiagramme der Varianzen nach jeder Ebene Ihrer Variablen zu erstellen (denn das ist es, was Sie tun Homogenität des Varianztests wird explizit bewertet). Sie können dies leicht tun, indem Sie die Varianz Ihres Ergebnisses für jede Kombination der Ebenen Ihrer Variablen schätzen und sie dann in Basis R zeichnen oder das ggplot2Paket verwenden.