library(learnr) library(ez) library(emmeans) library(gradethis) knitr::opts_chunk$set(exercise.checker = gradethis::grade_learnr) source('startup.R') data(conspiracy, package = 'PsyBSc7')
In der letzten Sitzung haben wir die einfaktorielle Varianzanalyse behandelt. Die spezifische Benennung als einfaktoriell verdeutlicht schon, dass wir hier ansetzen und Erweiterungen vornehmen können. In dieser Sitzung geht es vor allem um die zweifaktorielle Varianzanalyse. Ziel dieser Analyse ist es gleichzeitig Gruppenunterschiede auf mehreren Variablen zu untersuchen und dabei zu überprüfen, ob Kombinationen von Gruppen besondere Auswirkungen haben.
Wir arbeiten wieder mit dem conspiracy
Datensatz, der im Paket integriert ist:
data(conspiracy, package = 'PsyBSc7')
Wenn Sie diesen Datensatz runtergeladen haben und aus der Datei laden wollen, können Sie den üblichen load
-Befehl nutzen:
load('conspiracy.rda')
Eine kurze Übersicht über den Datensatz zeigt:
dim(conspiracy) head(conspiracy)
Der Datensatz enthält die Werte von 2451 Personen auf 9 Variablen. Er stammt aus einer Untersuchung zum Thema verschwörungstheoretische Überzegungen. Die ersten vier Variablen enthalten Informationen über den demographischen Hintergrund der Personen: höchster Bildungsabschluss, Typ des Wohnortes, Geschlecht und Alter. Die fünf restlichen Variablen sind Skalenwerte bezüglich verschiedener subdimensionen verschwörungstheoretischer Überzeugungen: GM (goverment malfeasance), MG (malevolent global conspiracies), ET (extraterrestrial cover-up), PW (personal well-being) und CI (control of information).
In der letzten Sitzung zeigte sich, dass die Überzeugung, dass die Existenz von Außerirdischen durch eine globale Verschwörung verdeckt wird (ET
), von dem höchsten Bildungsabschluss (edu
) und von der Art des Wohngebiets (urban
) abhängig ist. Zur Berechnung der einfaktoriellen ANOVAs wurde das ez
-Paket verwendet. Dieses Paket brauchen wir weiterhin:
library(ez)
Wie schon in der letzten Sitzung, ist es zunächst erforderlich eine Personen-ID zu erzeugen. In diesem Fall kann einfach die Zeilennummer einer Person genutzt werden:
conspiracy$id <- as.factor(1:nrow(conspiracy))
Führen Sie noch einmal die einfaktorielle ANOVA bezüglich des Wohnortes durch, um sich zu vergegenwärtigen, wie der ezANOVA
-Befehl funktioniert!
conspiracy$id <- as.factor(1:nrow(conspiracy))
# ezANOVA für einfaktorielle ANOVAs braucht vier Argumente: ezANOVA(data = ..., wid = ..., dv = ..., between = ...)
# wid kennzeichnet den Within-Identifier, eine Personen ID Variable
# dv ist die dependent (abhängige) Variable
# between ist die unabhängige Variable, die zwischen Personen unterscheidet # (eine Gruppierungsvariable)
conspiracy$id <- as.factor(1:nrow(conspiracy)) solu1 <- ezANOVA(conspiracy, wid = id, dv = ET, between = urban) solu2 <- ezANOVA(conspiracy, wid = id, dv = ET, between = urban, detailed = TRUE) solu3 <- ezANOVA(conspiracy, wid = id, dv = ET, between = edu) solu4 <- ezANOVA(conspiracy, wid = id, dv = ET, between = edu, detailed = TRUE) grade_result( pass_if(~ identical(.result, solu1), 'Durch den detailed = TRUE können Sie zusätzlich die Quadratsummen erhalten.'), pass_if(~ identical(.result, solu2)), fail_if(~ identical(.result, solu3), 'Sie haben die falsche Gruppierungsvariable genutzt!'), fail_if(~ identical(.result, solu4), 'Sie haben die falsche Gruppierungsvariable genutzt!'), fail_if(~ !(identical(.result, solu1) | identical(.result, solu2))), correct = 'Sehr gut! So wird die einfaktorielle ANOVA mit dem ez-Paket umgesetzt.', incorrect = 'Leider falsch.', glue_correct = '{.correct} {.message}', glue_incorrect = '{.incorrect} {.message}')
Die Ergebnisse zeigen, dass die Annahme der Homoskedastizität nicht verworfen werden muss und dass es bedeutsame Unterschiede zwischen verschiedenen Wohnorten hinsichtlich der verschwörungstheoretischen Überzeugung gibt, dass die Existenz Außerirdischer absichtlich verschleiert wird.
Die Ergebnisse aus den Übungsaufgaben ergaben bezüglich des Bildungabschlusses:
ezANOVA(conspiracy, wid = id, dv = ET, between = edu, white.adjust = TRUE)
Hier musste die Annahme der Homoskedastizität verworfen werden, sodass eine Adjustierung der Inferenzstatistik durchgeführt werden sollte. In beiden Fällen erhalten wir das generalisierte $\eta^2$, also einen Schätzer der Effektstärke. Um dies in der Skala der ursprünglichen Variablen zu interpretieren, können wir uns rein deskriptiv die gruppenspezifischen Mittelwerte angucken. Dazu kann mit einer Vielzahl von Funktionen gearbeitet werden. Eine gängige Variante ist der tapply
-Befehl:
tapply(conspiracy$ET, list(conspiracy$urban), mean) tapply(conspiracy$ET, list(conspiracy$edu), mean)
Probieren Sie gerne aus, ob Ihnen noch Alternativen einfallen! In den Tipps finden Sie zwei weitere Möglichkeiten:
# Mithilfe des aggregate-Befehls aggregate(conspiracy$ET, list(conspiracy$urban), mean) aggregate(conspiracy$ET, list(conspiracy$edu), mean)
# Mithilfe des describeBy-Befehls aus dem psych-Paket library(psych) describeBy(conspiracy$ET, conspiracy$urban) describeBy(conspiracy$ET, conspiracy$edu)
In der mehrfaktoriellen ANOVA steht nicht nur der Vergleich von Gruppen anhand einer unabhängigen Variable im Mittelpunkt, sondern der Fokus liegt auf der Kombination von Gruppierungen anhand mehrerer unabhängiger Variablen. Deksriptiv können die Mittelwerte aus Gruppenkombinationen ebenfalls mit der tapply
-Funktion bestimmt werden:
# Gruppierungskombinationen erstellen kombi <- conspiracy[, c('urban', 'edu')] # Kombinationsspezifische Mittelwertetabelle tapply(conspiracy$ET, kombi, mean)
Im ez
-Paket sind neben den Funktionen zur direkten Berechnung von Varianzanalysen auch einige zusätzliche Hilfefunktionen integriert. Dazu gehört auch die ezStats
-Funktion, die die Darstellung von Gruppengrößen, Mittelwerten und Standardabweichungen innerhalb der einzelnen Gruppenkombinationen erlaubt. Die Argumente, die diese Funktion erwartet sind analog zu denen in der ezANOVA
-Funktion:
data =
: der genutzte Datensatzwid =
: eine Personen ID-Variabledv =
: die abhängige Variable (dependent variable)between =
: eine Gruppierungsvariable (die zwischen Personen unterscheidet)Um mehrere Variablen als unabhängige Variablen zu deklarieren, kann mit dem c()
ein Vektor eröffnet werden, der an das Argument between
weitergegeben wird.
ezStats(conspiracy, dv = ET, wid = id, between = c(urban, edu))
Neben $N$, $\bar{x}$ und $\sigma$ wird in der Ausgabe auch Fishers Least Significant Difference ausgebeben. Diese kennzeichnet den minimalen Mittelwertsunterschied, der im direkten Vergleich zweier Gruppen signifkant wäre. Schon an dieser Stelle werden wir von ez
darauf hingewiesen, dass die Gruppen ungleich groß sind und dies in der ANOVA zu Problemen führen könnte.
Für eine grafische Darstellung der Mittelwerte, kann ezPlot
benutzt werden. Der Befehl nimmt die gleichen Argumente entgegen wie ezStats
, benötigt aber zusätzlich eine Aussage darüber, welche Variable auf der x-Achse abgetragen werden soll (x =
) und welche Variable farblich unterschieden werden soll (split =
). Mit dieser Funktion wird dann ein ggplot
erstellt. Diesen könnten Sie mit dem, was wir in der 2. Sitzung besprochen haben auch händisch erstellen! ezPlot
nimmt Ihnen hier aber ein wenig Arbeit ab:
ezPlot(conspiracy, dv = ET, wid = id, between = c(urban, edu), x = urban, split = edu)
Die FLSD wird hier in Form von Error-Bars dargestellt - durch diese kann also abgeschätzt werden, welche Mittelwerte sich statistisch bedeutsam unterscheiden.
Mithilfe der zweifaktoriellen Varianzanalyse können drei zentralen Fragen beantwortet werden (Eid, Gollwitzer & Schmitt, 2015, S. 432):
Im Beispiel wären die Fragen also:
ET
) auf die Art des Wohnorts (urban
) zurückführen?ET
) auf das Bildungsniveau (edu
) zurückführen?urban
) zwischen den Bildungsniveaus (edu
)?Deskriptiv lässt sich ein Hinweis auf eine Antwort zur 3. Frage in der Abbildung der Mittelwerte darin erkennen, dass innerhalb der Gruppe mit ländlichem Wohnort (rural
) Personen ohne Highschool Abschluss eine höhere verschwörungstheoretische Überzeugung aufweisen als Personen mit höheren Bildungsabschlüssen, wohingegen sie in den beiden anderen Wohnort-Gruppen einen niedrigeren ET-Wert aufweisen als Personen mit Highschool Abschluss.
Etwas technischer Ausgedrückt, lassen sich die drei Fragen in Hypothesenpaaren formulieren (Eid, Gollwitzer & Schmitt, 2015, S. 442):
Die drei Nullhypothesen werden in der zweifaktoriellen ANOVA geprüft. Die für ezStats
genutzten Argumente können auch für ezANOVA
benutzt werden. Um eine etwas detailliertere Ausgabe zu erhalten, kann zudem detailed = TRUE
gesetzt werden.
ezANOVA(conspiracy, dv = ET, wid = id, between = c(urban, edu), detailed = TRUE)
Der Levene Test fällt in diesem Fall statistisch bedeutsam aus, sodass die Homoskedastizitätsannahme (in diesem Fall: die Varianz ist in allen 9 Gruppen identisch) verworfen werden muss. ezANOVA
liefert eine eingebaute Korrekturmöglichkeit (HC3 von MacKinnon J. G. & White H., 1985), die mithilfe white.adjust = TRUE
angefordert werden kann. Probieren Sie das unten aus! Wie verändert das die Ergebnisse, die Ihnen ausgegeben werden?
ezANOVA(conspiracy, dv = ET, wid = id, between = c(urban, edu), detailed = TRUE)
In diesem Fall werden beide Haupteffekte statistisch bedeutsam, die Interaktion allerdings nicht. Inhaltlich heißt das, dass sowohl die Art des Wohnorts als auch das Bildungsniveau einen Einfluss auf die verschwörungstheoretische Überzeugung haben. Über die jeweiligen Effekte hinaus, ist die spezifische Kombination aus Wohnort und Bildungsniveau für diese Überzeugung irrelevant.
Mit Tukeys Honest Significant Difference können, wie auch letzte Sitzung, alle möglichen Gruppenkombinationen verglichen werden.
TukeyHSD(aov(ET ~ urban*edu, conspiracy))
Leider ist das Ergebnis etwas unübersichtlich, weil sich in diesem Fall 36 Vergleiche ergeben. Mit ein paar Funktionen aus dem emmeans
-Paket können wir versuchen, das optisch etwas aufzubereiten. Dafür zunächst das Paket laden:
library(emmeans)
In diesem Paket gibt es, die wenig überraschend benannte emmeans
-Funktion, mit der wir alle weiteren Analysen vorbereiten müssen:
emm <- emmeans(aov(ET ~ urban*edu, conspiracy), ~ urban * edu)
Diese Funktion nimmt als erstes Argument das gleiche aov
-Objekt entgegen wie TukeyHSD
. Als Zweites müssen wir definieren, welche unabhängigen Variablen uns in der Post-Hoc Analyse interessieren. Weil wir hier alle Gruppen betrachten möchten, können wir einfach die gleiche Struktur der unabhängigen Variablen wiederholen: ~ urban * edu
.
Wenn wir uns das entstandene Objekt angucken, sehen wir eine Tabelle mit 7 Spalten:
emm
Die ersten beiden Spalten (urban
, edu
) geben an, welche Ausprägungen unsere beiden unabhängigen Variablen haben. Die erste Zeile bezieht sich also auf Personen aus einem ländlichen Gebiet, die keinen Highschool-Abschluss haben. Die dritte Spalte (emmean
) ist der Gruppenmittelwert, die vierte der dazugehörige Standardfehler (SE
), dann die fünfte die Freiheitsgrade (df
) und in die letzten beiden Spalten geben das Konfidenzintervall des Mittelwerts an (lower.CL
, upper.CL
).
Mittelwerte und Konfidenzintervalle können wir uns sehr einfach direkt plotten lassen:
plot(emm)
Diese Abbildung können wir um eine Aussage über die direkten Vergleiche erweitern:
plot(emm, comparisons = TRUE)
Die neu hinzugekommenen roten Pfeile geben uns einen Hinweis dazu, welche Gruppen sich unterscheiden. Wenn zwei rote Pfeile überlappen, gibt es keinen statistisch bedeutsamen Unterschied. Wenn Sie das nicht tun, unterscheiden sich die beiden Gruppenmittelwerte auf dem festgeleten $\alpha$-Fehlerniveau (per Voreinstellung 5%) statistisch bedeutsam.
Eine zweite Möglichkeit, die Ergebnisse ein wenig übersichtlicher zu gestalten sind pairwise $p$-value plots. Im emmeans
-Paket werden diese über pwpp
angefordert:
pwpp(emm)
In dieser Abbildung ist auf der x-Achse der $p$-Wert des Mittelwertvergleichs dargestellt. Auf der y-Achse werden die Gruppen anhand ihrer deskriptiven Mittelwerte sortiert und abgetragen (dabei sind alle Abstände zwischen zwei Gruppen gleich groß, egal wie groß der Mittelwertsunterschied auf der abhängigen Variable tatsächlich ist).
In Situationen mit vielen Gruppen ist es außerordentlich ineffizient alle Vergleiche durchzuführen. Durch die Grundgedanken des Nullhypothesentestens muss das $\alpha$-Fehlerniveau auf alle durchgeführten Tests korrigiert werden. Nach klassischer Bonferroni-Korrektur wäre das korrigierte $\alpha$-Fehlerniveau in diesem Fall also $\frac{.05}{36} = r round(.05/36, 4)
$ um eine echte Irrtumswahrscheinlichkeit von 5% aufrecht zu erhalten. An dieser Stelle können daher geplante Kontraste genutzt werden, um a-priori definierte, theoretisch relevante Vergleiche durchzuführen. So kann die Anzahl der Tests auf die theoretisch notwendigen reduziert werden und mit einer weniger restriktiven Korrektur gerechnet werden.
Um Kontraste definieren zu können, müssen wir zunächst in Erfahrung bringen, in welcher Reihenfolge die Gruppenkombinationen intern repräsentiert werden. Diese Reihenfolge haben wir bereits im emm
-Objekt gesehen:
emm
Mithilfe eines 9 Elemente langen Vektors können Kontraste festgelegt werden. Um z.B. die Gruppe "rural, not highschool" mit der Gruppe "suburban, not highschool" zu vergleichen, kann folgender Vektor angelegt werden:
cont1 <- c(1, -1, 0, 0, 0, 0, 0, 0, 0)
Die Nullhypothese, die durch diesen Vektor geprüft wird, lässt sich mithilfe der Reihenfolge der Gruppen leicht zusammenstellen. Wenn $j$ die drei Stufen von urban
indiziert (1 = rural, 2 = suburban, 3 = urban) und $k$ die drei Stufen von edu
(1 = not highschool, 2 = highschool, 3 = college), ist die durch cont1
festgelegte Nullhypothese:
$H_0: 1 \cdot \mu_{11} - 1 \cdot \mu_{21} + 0 \cdot \mu_{31} + 0 \cdot \mu_{12} + 0 \cdot \mu_{22} + 0 \cdot \mu_{32} + 0 \cdot \mu_{13} + 0 \cdot \mu_{23} + 0 \cdot \mu_{33} = 0$
Oder gekürzt:
$H_0: \mu_{11} - \mu_{21} = 0$
Mit dem contrast
-Befehl kann der festgelegte Kontrast geprüft werden:
contrast(emm, list(cont1))
Dieser Kontrast entspricht dem ersten Vergleich des oben durchgeführten TukeyHSD
, unterscheidet sich jedoch im $p$-Wert. Der hier bestimmte $p$-Wert ist nicht korrigiert (weil nur ein Kontrast geprüft wurde), der oben aufgeführte ist hingegen auf 36 Tests Tukey-korrigiert.
Probieren Sie es einmal selbst aus, indem Sie einen Kontrast aufstellen, der den Unterschied zwischen Personen mit Highschool-Abschluss und Personen ohne Highschool-Abschluss, die in Städten wohnen, prüft!
emm <- emmeans(aov(ET ~ urban*edu, conspiracy), ~ urban * edu) cont1 <- c(1, -1, 0, 0, 0, 0, 0, 0, 0)
# In emm sehen sie die Reihenfolge der Gruppen
emm
# Relevant sind die Gruppen urban - not highschool und urban - highschool
# Stellen Sie sicher, dass ihr Kontrastvektor in der Summe 0 ergibt # Für das vorherige Beispiel: sum(cont1) == 0
Mithilfe der Kontrast-Vektoren können auch komplexe Hypothesen geprüft werden. Beispielsweise könnten wir vergleichen, inwiefern sich Personen aus städtischer Umgebung ($j = 3$) mit mindestens High School Abschluss von Personen ohne High School Abschluss unterscheiden:
cont2 <- c(0, 0, 1, 0, 0, -.5, 0, 0, -.5)
oder in Hypothesenform: $H_0: \mu_{31} - .5 \cdot \mu_{32} - .5 \cdot \mu_{33} = 0$ bzw. $H_0: \mu_{31} - \frac{\mu_{32} + \mu_{33}}{2} = 0$. Eine generelle Daumenregel besagt, dass die Summe des Kontrastvektors 0 sein sollte.
Weil sowohl cont1
als auch cont2
durchgeführt werden, muss für das multiple Testen der beiden korrigiert werden. Das kann dadurch erreicht werden, dass im contrast
-Befehl alle Kontraste gleichzeitig eingeschlossen werden und mit adjust = 'bonferroni'
z.B. die Bonferroni-Korrektur ausgewählt wird:
contrast(emm, list(cont1, cont2), adjust = 'bonferroni')
Nutzen Sie für die Aufgaben den Datensatz nature
.
data(nature, package = 'PsyBSc7') str(nature)
Dieser Datensatz enthält Beobachtungen von 743 Personen auf 8 Variablen und ist ein Teil eines größeren Datensatzes, der auf https://openpsychometrics.org/ zur Verfügung gestellt wird.
Die sechs Variablen Q1A bis Q6A beziehen sich auf die Skala zur Erhebung der Naturverbundenheit (z.B. "My relationship to nature is an important part of who I am") von Nisbet & Zelenski (2013). Das Antwortformat ist 5-stufig, wobei höhere Werte stets eine stärkere Naturverbundenheit zum Ausdruck bringen. Darüber hinaus enthält der Datensatz zwei nominalskalierte Variablen: urban
- der Typ des Ortes, an dem eine Person aufgewachsen ist - und continent
- der Kontinent auf dem die teilnehmende Person aufgewachsen ist.
NV
(NaturVerbundenheit) an. Erstellen Sie außerdem - zur Vorbereitung der ez
-Funktionen eine ID-Variable mit dem Namen ID
.urban
und continent
.Bei mehrfaktoriellen ANOVAs können die Quadratsummen auf unterschiedliche Arten berechnet werden. Verbreitet sind dabei 3 Typen, zwischen denen man sich anhand der inhaltlichen Hypothesen entscheiden sollte.
Typ I berücksichtigt in der Berechnung der Quadratsummen nur die vorherigen unabhängigen Variablen. Dies entspricht konzeptuell der sequentiellen Aufnahme von Prädiktoren in der Regression.
# QS-Typ 1, Reihenfolge 1 ezANOVA(conspiracy, dv = ET, wid = id, between = c(urban, edu), type = 1) # QS-Typ 1, Reihenfolge 2 ezANOVA(conspiracy, dv = ET, wid = id, between = c(edu, urban), type = 1)
Typ II berücksichtigt in der Berechnung alle anderen unabhängigen Variablen. In der Berechnung der einzelnen Quadratsummen wird allerdings angenommen, dass alle Interaktionen, an denen dieser Term beteiligt ist, 0 sind. Typ II ist in ezANOVA
voreingestellt.
# QS-Typ 2 ezANOVA(conspiracy, dv = ET, wid = id, between = c(urban, edu), type = 2)
Typ III unterscheidet sich von Typ II nur darin, dass bei der Berechnung nicht angenommen wird, dass die Interaktionen 0 sind. Typ III ist z.B. in SPSS voreingestellt.
# QS-Typ 3 ezANOVA(conspiracy, dv = ET, wid = id, between = c(urban, edu), type = 3)
Generell ist Typ II besser geeignet um die Quadratsummen von Haupteffekten zu bestimmen, wenn Interaktionen empirisch nicht von 0 verschieden sind. Wenn Interaktionen von 0 verschieden sind, wird (unabhängig vom QS-Typ) davon abgeraten die Haupteffekte zu interpretieren, sodass deren Bestimmung in diesem Fall wenig Relevanz hat. Die Terme höchster Ordnung (hier die Interaktion) sind zwischen Typ II und Typ III identisch, sodass die Interpretation der Interaktion durch die Wahl nicht beeinflusst wird. Von Typ I wird generell abgeraten (wie Ihnen ezANOVA
auch direkt mitteilt).
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.