Statistik II
Varianzanalyse
AG Forschungsmethoden und Kognitive Psychologie, Institut für Psychologie, Universität Bremen
Willkommen bei Stats by Randolph. Hier geht es zur Hauptseite mit weiteren Informationen und Inhalten.
Autor:innen dieser Seite: An den Inhalten dieser Seite haben mitgearbeitet: Markus Janczyk und Valentin Koob. Der Inhalt dieser Seite wird in der Lehre in den Studiengängen Psychologie von der AG Forschungsmethoden und Kognitive Psychologie an der Universität Bremen verwendet, steht aber allen Interessierten zur Verfügung. Rückmeldungen/Fehler/Vorschläge können gesendet werden an randolph@uni-bremen.
Versionshistory:
- v1.0: erste online-gestellte Version (14.2.2024)
2 Einfaktorielle Varianzanalyse I
Bisher haben wir in Statistik I den Fall betrachtet, dass maximal zwei (un-)abhängige Stichproben erhoben und ausgewertet wurden (mit dem entsprechenden \(t\)-Test). Dies deckt bereits eine ganze Menge an möglichen Untersuchungen ab und manchmal können komplexere Designs auch auf diese Situation und einen \(t\)-Test heruntergebrochen werden. Allerdings ist es natürlich auch möglich, Studien mit mehr als zwei Gruppen/Bedingungen durchzuführen. Für diesen Fall stellt die Varianzanalyse (englisch: Analysis of Variance, abgekürzt daher oft als ANOVA) Möglichkeiten zur Auswertung bereit. Die Teile 2-5 von Statistik II befassen sich daher mit verschiedenen Varianten der Varianzanalyse.
Wir beginnen hier in Teil 2 mit dem einfachsten Fall der einfaktoriellen between-subject Varianzanalyse. Hier spielen bereits zwei Begriffe eine wichtige Rolle:
- einfaktoriell: Es gibt nur eine unabhängige Variable mit (i.d.R.) mehr als zwei Stufen. Da unabhhängige Variablen in der Varianzanalyse oft “Faktoren” genannt werden, reden wir von einer einfaktoriellen Varianzanalyse.
- between-subject: Jede Stufe oder Bedingung des Fakors wird an einer separaten Stichprobe erhoben und die Stichproben sind unabhängig voneinander gezogen (vgl. \(t\)-Test für zwei unabhängige Stichproben).
Streng genommen kommt noch hinzu, dass wir nur eine abhängige Variable betrachten und der Fall daher als univariat bezeichnet wird.
2.1 Grundlagen
2.1.1 Datensituation und Beispiel
Als ein praktisches Beispiel betrachten wir eine (fiktive) Studie, die untersucht, ob eine wiederholte Anwendung eines Sprachtrainings bei Kindern deren Fähigkeit in einem bestimmten sprachlichen Bereich mehr fördert, als die einmalige Anwendung. Dazu wurden drei Gruppen von jeweils \(n=10\) Kindern in die Studie mit einbezogen: eine Gruppe bekam das Training nur einen Tag, eine andere Gruppe zwei Tage und wiederum eine andere Gruppe drei Tage. Im Anschluss wurde bei allen Kindern ein Sprachtest durchgeführt, der einen Wert von 0-20 ergeben konnte (wobei höhere Werte für eine bessere Leistung stehen). Zusammengefasst haben wir also folgende Situation:
- abhängige Variable: Leistung im Sprachtest (0-25)
- Anzahl der Tage Training als unabhängige Variable/Faktor mit drei
Stufen:
- 1 Tag
- 2 Tage
- 3 Tage
- für jede der drei Stufen wurde eine unabhängige Stichprobe vom Umfang \(n=10\) gezogen
- zusätzlich wurde das Alter der Kinder erhoben (da wir es für Teil 4 benötigen werden)
Die resultierenden Daten inkl. der Gruppenmittelwerte und -varianzen sind in folgender Tabelle zusammengefasst (und können über den Link ganz am Anfang des Dokuments heruntergeladen werden):
Als erstes wollen wir die Daten visualisieren. Dazu laden wir die Daten
zunächst und berechnen die Mittelwerte für die drei Gruppen im
Anschluss:
# ./ als Beginn des Suchpfads bedeutet "beginne beim aktuellen Working Directory"
daten_training <- read.table("./Daten/daten_training.DAT",
sep = "\t", # Tab-getrennt (Default)
header = TRUE) # Variablennamen in Zeile 1
# hier werden die ersten 6 Versuchspersonen angezeigt
head(daten_training)# und wir berechnen die drei Mittelwerte
# in Klammern gesetzt: Wert wird direkt auch ausgegeben
( mittelwerte <- tapply(daten_training$Sprachtest,
daten_training$Trainingstage,
FUN = mean) )## 1 2 3
## 13.1 15.7 18.1Wir stellen diese Mittelwert nun grafisch dar und sehen, dass mit
zunehmender Trainingsdauer die Leistung im Sprachtest offenbar zunimmt:
Wir werden auf diese Daten später zurückkommen und entsprechende Varianzanalysen von Hand und mit R berechnen. Zuvor gehen wir noch einmal auf zentrale Begriffe und die Hypothesen der Varianzabalyse ein und klären die Frage, warum zum Testen der Nullhypothese nicht auch eine Reihe von \(t\)-Test durchgeführt werden sollte.
2.1.2 Zentrale Begriffe und Hypothesen
Im Kontext der Varianzanalyse werden manche Begriffe sehr speziell verwendet (wie bereits eingangs am Beispiel des “Faktors” erwähnt). Die folgenden zentralen Begriffe sind für unsere Belange die zunächst wichtigsten:
- Faktor: Bezeichnung für die unabhängige (oft
manipulierte) Variable in der Varianzanalyse
- im Beispiel: 1 Faktor “Tage Training” mit 3 Stufen
- allgemein hat ein Faktor \(J\) Stufen (im Beispielalso \(J=3\))
- Abhängige Variable: Variable, die gemessen wird
(üblicherweise als \(Y\) bezeichnet)
- im Beispiel: Leistung im Sprachtest
- Die Stichproben für jede der \(J\) Faktorstufen müssen nicht gleich groß sein, d.h. wir schreiben den Umfang jeder Stichprobe allgemein als \(n_j\).
Im allgemeinen Fall meinen wir im Folgenden mit \[y_{ij}\] den Wert der Person \(i\) in der Stufe/Stichprobe \(j\). Der allgemeine Fall mit \(J\) Stufen und jeweils \(n_j\) Versuchspersonen pro Stufe würde sich als Matrix also wie folgt darstellen:
\[\begin{equation*} \begin{array}{ccccc}\hline {\text{Stufe 1}} & \cdots & {\text{Stufe }j} & \cdots & {\text{Stufe }J} \\\hline y_{11} & \cdots & y_{1j} & \cdots&y_{1J}\\ \vdots& & \vdots& &\vdots\\ y_{i1}& \cdots &y_{ij}& \cdots &y_{iJ}\\ \vdots& & \vdots& &\vdots\\ y_{{n_1}1} & \cdots & y_{{n_j}j} &\cdots& y_{{n_J}J}\\\hline \end{array} \end{equation*}\]
Wie bereits beim \(t\)-Test auch, gehen wir von Anfang an von bestimmten Voraussetzungen bzgl. der Daten aus. Tatsächlich sind diese Voraussetzungen denen für den \(t\)-Test für zwei unabhängige Stichproben sehr ähnlich (siehe hier):
- Die \(J\) Stichproben sind unabhängig voneinander gezogen.
- Die abhängige Variable \(Y\) ist (mindestens) intervallskaliert.
- Die abhängige Variable ist in jeder Population \(j\) normalverteilt mit einer gemeinsamen Varianz \(\sigma^2\), d.h. “Varianzhomogenität” liegt vor, wobei die jeweiligen Erwartungswerte \(\mu_j\) sich unterscheiden können: \[Y_j \sim N(\mu_j, \sigma^2)\hspace{1cm}\forall j\in \{1,\ldots,J\}\]
Nun kommen wir auf die Hypothesen zu sprechen, die mit der Varianzanalyse getestet werden können (siehe hier für die Einführung aus Statistik I). Da die (between-subject) Varianzanalyse eine Verallgemeinerung des \(t\)-Tests für zwei unabhängige Stichproben auf mehr als zwei unabhängige Stichproben ist, erinnern wir uns die Hypothesen des \(t\)-Tests: Die Nullhypothese hatte im Wesentlichen die Gleichheit von Erwartungswerten eingeschlossen während die Alternativhypothese einen Unterschied postuliert hatte. Ganz ähnlich formulieren wir nun beide Hypothesen auch für den hier betrachteten Fall der Varianzanalyse mit \(J\) Stufen, indem wir die Nullhypothese verallgemeinern. Die Alternativhypothese ergibt sich dann als das logische Gegenteil der Nullhypothese: \[\begin{equation*} \begin{aligned} &H_0:\mu_1=\mu_2=\ldots=\mu_J\\ &H_1:\mu_k\neq\mu_m\text{ für mindestens ein Paar }k,m\in\{1,\ldots,J\} \text{ mit } k\neq m \end{aligned} \end{equation*}\]
Es bietet sich hier direkt an, sich klar zu machen, was diese Hypothesen eigentlich bedeuten. Die Nullhypothese \(H_0\) ist recht einfach. Sie besagt, dass alle \(\mu_j\) identisch sein müssen. Wichtig ist aber, dass die Alternativhypothese \(H_1\) nicht postuliert, dass alle \(\mu_j\) voneinander verschieden sein müssen. Vielmehr gilt die Nullhypothese bereits nicht mehr, wenn bereits zwei Erwartungswerte sich unterscheiden. Dies schließt natürlich nicht aus, dass sich alle Erwartungswerte unterscheiden; sie müssen es aber eben nicht nicht.
Nun betrachten wir wieder das Eingangsbeispiel mit \(J=3\) Gruppen bzw. Populationen. Die Nullhypothese würde nun also besagen \[ H_0:\mu_1=\mu_2=\mu_3 \] Diese Nullhypothese hat aber ihre Richtigkeit auch dann, wenn drei andere Nullhypothesen gelten würden, nämlich: \[ \mu_1=\mu_2\quad\text{und}\quad\mu_1=\mu_3\quad\text{und}\quad\mu_2=\mu_3 \] Warum können (oder besser: sollten) wir daher nicht einfach drei \(t\)-Tests durchführen und, wenn alle nicht-signifikant sind, darauf schließen, dass auch die \(H_0:\mu_1=\mu_2=\mu_3\) gilt? Dies bringt uns zum Konzept der \(\alpha\)-Kumulation und dem damit eng verwandten Konzept der \(\alpha\)-Adjustierung.
2.1.3 \(t\)-Test vs. Varianzanalyse: \(\alpha\)-Kumulation
Der Haken an der möglichen alternativen Vorgehensweise ist eine wahrscheinlichkeitstheoretische Einsicht: Die Wahrscheinlichkeit, dass ein Test (fälschlicherweise) signifikant wird, obwohl alle Nullhypothesen gelten, ist kleiner als die Wahrscheinlichkeit, dass mindestens einer der Tests (fälschlicherweise) signifikant wird, obwohl alle Nullhypothesen gelten.
Fragen wir uns als erstes: Wie groß ist die Wahrscheinlichkeit, dass ein Test fälschlicherweise signifikant wird, also ein Fehler 1. Art auftritt? Die Antwort auf diese Frage kennen wir bereits: \[\alpha\] Nun kann jeder einzelne Test als ein Bernoulli-Experiment aufgefasst werden, bei der ein “Erfolg” mit der Wahrscheinlichkeit \(\alpha\) eintritt. Betrachten wir die Serie an \(t\)-Tests als eine Bernoulli-Folge), dann können wir die Wahrscheinlichkeit dafür, dass \(m\) von \(k\) Tests signifikant werden, mit der Binomialverteilung berechnen.
Darüber hinaus ist die gesuchte Wahrscheinlichkeit \(P(m\geq 1)\), also dass mindestens einer der \(k\) Test signifikant wird, gleich der Gegenwahrscheinlichkeit, dass keiner der \(k\) Tests fälschlicherweise signifikant wird, also: \[ P(m\geq 1)=1-P(m=0) \] Nun haben wir alle beisammen, um die gesuchte Wahrscheinlichkeit zu berechnen: \(P(m=0)\) ergibt sich aus der Binomialverteilung wegen \({{k}\choose{0}}=1\) und \(\alpha^0=1\) folgt \[\begin{equation*} \begin{aligned} P(m=0)=&{{k}\choose{m}}\cdot \alpha^m\cdot (1-\alpha)^{k-m}\\ =&{{k}\choose{0}}\cdot \alpha^0\cdot (1-\alpha)^{k-0}\\ =&(1-\alpha)^k \end{aligned} \end{equation*}\] Damit ergibt sich insgesamt also: \[ \boxed{P(m\geq 1)=1-P(m=0)=1-(1-\alpha)^k} \] Nun kehren wir zum Eingangsbeispiel mit den \(k=3\) Tests zurück und fragen uns, wie groß die Wahrscheinlichkeit, dass mindestens einer der drei Tests signifikant wird ist, wenn wir für jeden einzelnen Test \(\alpha=0.05\) annehmen. Eingesetzt in die gerade hergeleitete Formel ergibt sich: \[ P(m\geq 1)=1-(1-0.05)^3=1-0.95^3=0.14 \] Dies bedeutet: Bei drei \(t\)-Tests, bei denen jeweils die Nullhypothese gilt, ist die Wahrscheinlichkeit mindestens einer Fehlentscheidung im Sinne eines Fehler 1. Art also bereits fast dreimal so groß wie \(\alpha\). Mit anderen Worten: Diese Vorgehensweise zu nutzen, um die \(H_0:\mu_1=\mu_2=\mu_3\) zu testen, erhöht die Gefahr einer Fehlentscheidung deutlich.
Um dem entgegenzuwirken nutzt man sog. \(\alpha\)-Adjustierungen. Das Grundprinzip hierbei ist, dass durch die Wahl eines kleineren \(\alpha'\) für jeden einzelnen Test gewährleistet wird, dass die Gesamtwahrscheinlichkeit mindestens einer Fehlentscheidung das eigentlich beabsichtige \(\alpha\) nicht überschreitet. Die vielleicht bekannteste Variante ist die sog. Bonferroni-Korrektur (elaboriertere Varianten der Korrektur finden sich ausführlich z.B. bei Baguley, 2012, ab S. 490): \[\alpha'=\frac{\alpha}{k}\] Für unser Beispiel mit den \(k=3\) Tests ergibt sich damit als neues \(\alpha'\): \[\alpha'=\frac{0.05}{3}=0.01\bar{6}\] Das heißt: Im Prinzip könnten wir tatsächlich jeden der drei \(t\)-Tests mit \(\alpha'\) durchführen. Die bessere Alternative ist allerdings eine Varianzanalyse, da diese nur einen Test für diese Hypothese benötigt. Darüber hinaus werden wir in den weiteren Teilen, dass die Varianzanalyse weiter ausgebaut werden kann und dann eine Menge weiterer Vorteile bietet.
2.1.4 Die konzeptuelle Idee der Varianzanalyse
Bevor wir uns an die Formeln einer Varianzanalyse und die
rechnerische Durchführung begeben, wollen wir nun die konzeptuelle Idee
einer Varianzanalyse verstehen. Die folgende Abbildung visualisiert in
beiden Diagrammen die drei Gruppenmittelwerte aus dem Eingangsbeispiel.
Die Fehlerbalken stellen (hypothetische) Varianzen dar, die als
Schätzung für die Populationsvarianz genutzt werden. Wir sehen also,
dass die Populationen im linken Diagramm eine größere Varianz haben als
im rechten Diagramm. Bei welcher Situation würde man nun intuitiv eher
an einen Unterschied der Populationsmittelwerte \(\mu_j\) glauben?
Vermutlich ist die intuitive Antwort, dass in der rechten Abbildung eher
von einem wirklichen Unterschied in den Populationsparametern auszugehen
ist: Da die Varianz sehr klein ist, liegen die einzelnen Werte also auch
eher nahe beieinander. Daher dürften auch die Mittelwerte eher nur in
einem kleineren Bereich liegen und die abgebildeten Gruppenmittelwerte
dürften daher die Populationsmittelwerte relativ gut repräsentieren. Bei
einer größeren Varianz wie im linken Teil der Abbildung liegen die
einzelnen Werte ja generell auch weiter auseinander. Dementsprechend
haben auch die Mittelwerte von Stichproben einen eher größeren Spielraum
was ihre konkreten Werte angeht.
Ausgehend von dieser Überlegung können wir zwei Charakteristika von Daten ermitteln, die wichtig für die Entwicklung einer geeigneten Prüfgröße sind:
- Relativ klar ist, dass große Unterschiede in den Gruppen-/Bedingungsmittelwerte auch eher für eine Unterschiedlichkeit Populationsmittelwerte sprechen.
- Die Gruppen-/Bedingungsmittelwerte sind in den beiden Abbildungen allerdings identisch, aber es kommt eben auch das gerade vorgebrachte Argument zur Geltung: Bei großer (Populations-)Varianz sind größere Unterschiede der Mittelwerte sowieso wahrscheinlich(er). Daher würden wir–bei gegebenen Mittelwerten–eher von Unterschieden in den Populationsmittelwerten ausgehen, wenn die Varianzen der Populationen eher gering ausfallen.
Um Punkt 2 noch besser zu illustrieren, schauen wir die Ergebnisse verschiedener Simulationen an. Zunächst betrachten wir dazu den Betrag von Mittelwertunterschieden zwischen zwei Gruppen (bei gültiger Nullhypothese) in Abhängigkeit verschiedener Populationsvarianzen (jeder Punkt repräsentiert den Mittelwert aus 5,000 Ziehungen):
n.repetitions <- 5000 # wie viele Wiederholungen pro Varianzwert
vars <- seq(500, 10000, by = 500) # welche Varianzen sollen betrachtet werden?
results_diff <- NULL # Container zum Speichern pro Varianzwert
set.seed(1) # seed setzen für exakte Replizierbarkeit
for (one_var in vars) { # für jede Varianz ...
abs_differences <- NULL # Container für die absoluten Differenzen pro Wiederholung
for (i in 1:n.repetitions) { # ... tue x-mal
one_sd <- sqrt(one_var)
m_sample1 <- mean(rnorm(20, 0, one_sd)) # Mittelwert der ersten Stichprobe
m_sample2 <- mean(rnorm(20, 0, one_sd)) # Mittelwert der zweiten Stichprobe
abs_diff <- abs(m_sample1 - m_sample2) # absolute Differenz berechnen
abs_differences <- c(abs_differences, abs_diff) # und zwischenspeichern pro Wiederholung
}
# letztlich den Mittelwert der absoluten Differenzen speichern
results_diff <- c(results_diff, mean(abs_differences))
}
par(mar = c(5, 5, 2, 1), mfrow = c(1, 1))
plot(vars, results_diff, type = "p", pch = 19, xlab = "Populationsvarianz",
ylim = c(0, 50), xlim = c(0, 10000), ylab = "Differenz der 2 Mittelwerte (Betrag)",
cex = 1.5, cex.lab = 1.5, axes = FALSE)
axis(1, las = 1, at = seq(0, 10000, by = 2000), cex.axis = 1.5) # x-Achse
axis(2, las = 1, at = seq(0, 50, by = 10), cex.axis = 1.5) # y-Achse
abline(h = 0)
Der wichtige Punkt an dieser Abbildung ist, dass die
Mittelwertdifferenzen tatsächlich ohnehin größer werden, wenn die
Populationsvarianz ansteigt. Bei mehr als zwei Gruppen/Bedingungen ist
eine Differenz allerdings nicht mehr zu berechnen. Allerdings bildet
auch die Varianz der Mittelwerte (1) die Unterschiedlichkeit der
Mittelwerte ab und (2) sie wäre zudem auch bei mehr als zwei
Gruppen/Bedingungen anwendbar. Die folgende Abbildung vergleicht
nocheinmal die Mittelwertdifferenz (links) mit der für die gleichen
Daten berechneten Varianz der beiden Mittelwerte (rechts). Offenbar
nimmt auch die Varianz von zwei Mittelwerten zu, wenn die zugrunde
liegende Populationsvarianz steigt:
Nun erweitern wir die Situation auf 3 Gruppen (rechter Teil der nächsten
Abbildung) und sehen ein ganz ähnliches Muster (im linken Teil ist zum
Vergleich noch einmal die Situation mit 2 Gruppen dargestellt):
n.repetitions <- 5000 # wie viele Wiederholungen pro Varianzwert
vars <- seq(500, 10000, by = 500) # welche Varianzen sollen betrachtet werden?
results_vars_12 <- NULL # NEU: Container zum Speichern für 2 Gruppen
results_vars <- NULL # NEU: Container zum Speichern für 3 Gruppen
set.seed(1) # seed setzen für exakte Replizierbarkeit
for (one_var in vars) { # für jede Varianz ...
vars_ms_12 <- NULL # NEU: Container der Varianz von 2 Mittelwerten
vars_ms <- NULL # NEU: Container der Varianz von 3 Mittelwerten
for (i in 1:n.repetitions) { # ... tue x-mal
one_sd <- sqrt(one_var)
m_sample1 <- mean(rnorm(20, 0, one_sd)) # Mittelwert der ersten Stichprobe
m_sample2 <- mean(rnorm(20, 0, one_sd)) # Mittelwert der zweiten Stichprobe
m_sample3 <- mean(rnorm(20, 0, one_sd)) # NEU: Mittelwert der dritten Stichprobe
var_m_12 <- var(c(m_sample1, m_sample2)) # NEU: Varianz der Mittelwerte für 2 Gruppen
var_m <- var(c(m_sample1, m_sample2, m_sample3)) # NEU: ... für 3 Gruppen
vars_ms_12 <- c(vars_ms_12, var_m_12) # NEU: zwischenspeichern pro Wiederholung
vars_ms <- c(vars_ms, var_m) # NEU: zwischenspeichern pro Wiederholung
}
# letztlich den Mittelwert der Varianzen speichern
results_vars_12 <- c(results_vars_12, mean(vars_ms_12))
results_vars <- c(results_vars, mean(vars_ms))
}
par(mar = c(5, 5, 2, 2), mfrow = c(1, 2))
plot(vars, results_vars_12, type = "p", pch = 19, xlab = "Populationsvarianz",
ylim = c(0, 500), xlim = c(0, 10000), ylab = "Varianz der 2 Mittelwerte",
cex = 1.5, cex.lab = 1.5, axes = FALSE)
axis(1, las = 1, at = seq(0, 10000, by = 2000), cex.axis = 1.5) # x-Achse
axis(2, las = 1, at = seq(0, 500, by = 100), cex.axis = 1.5) # y-Achse
abline(h = 0)
plot(vars, results_vars, type = "p", pch = 19, xlab = "Populationsvarianz",
ylim = c(0, 500), xlim = c(0, 10000), ylab = "Varianz der 3 Mittelwerte",
cex = 1.5, cex.lab = 1.5, axes = FALSE)
axis(1, las = 1, at = seq(0, 10000, by = 2000), cex.axis = 1.5) # x-Achse
axis(2, las = 1, at = seq(0, 500, by = 100), cex.axis = 1.5) # y-Achse
abline(h = 0)
Zusammengefasst bedeutet dies, dass für eine Unterschiedlichkeit der Populationsmittelwerte (d.h. für die Alternativ- und gegen die Nullhypothese) zwei Aspekte sprechen:
- Eine große Unterschiedlichkeit (d.h. eine große Varianz) der Gruppen- (bzw. Populations-)Mittelwerte
- Eine kleine Varianz innerhalb jeder Gruppe bzw. Population, da bei großer Varianz die Mittelwerte ohnehin eine größere Varianz, und damit Unterschiedlichkeit, aufweisen.
Aufbauend auf diesen Überlegungen konstruieren wir nun eine Prüfgröße (vgl. auch Teil 12 in Statistik I), die die beiden bekannten Eigenschaften erfüllt:
- Sie vereint die beiden Punkte von gerade eben so, dass sie umso extremere Werte annimmt, je mehr die Daten gegen die Nullhypothese sprechen.
- Aufgefasst als Zufallsvariable ist ihre Dichtefunktion bestimmbar.
Um den 1. Punkt zu erfüllen, eignet sich ein Bruch in besonderer Weise: \[ F=\frac{\text{Variabilität der Stichproben-/Populationsmittelwerte}}{\text{Variabilität in den Stichproben/Populationen}} \] Dieser Bruch \(F\) wird größer, wenn der Zähler groß wird und wenn der Nenner klein wird. Die Verteilung, der die entsprechende Zufallsvariable folgt wird die sog. \(F\)-Verteilung sein.
In etwas anderen Worten kann der \(F\)-Bruch auch beschrieben werden als \[ F=\frac{\text{Variabilität zwischen den Gruppen}}{\text{Variabilität innerhalb der Gruppen}} \] Was sind schließlich die Quellen dieser beiden Variabilitäten?
- Für die Variabilität innerhalb der Gruppen spielen (vor allem) der Messfehler eine Rolle und nicht-interessierende individuelle Unterschiede: Einerseits können wir den wahren Wert einer Versuchsperson nicht genau messen, sodass wir eine unsystematische Schwankung von Werten aufgrund von Messfehler haben. Andererseits haben nicht alle Personen, die aus der gleichen Population stammen, den gleichen Wert auf der abhängigen Variable. Das heißt, der wahre Wert mancher Personen liegt oberhalb des Populationsmittelwerts, während die Werte anderer Personen darunter liegen. Die Variabilität innerhalb der Gruppen stellt also eine unsystematische Schwankung von Werten innerhalb einer Stichprobe aufgrund von individuellen Unterschieden sowie einem Messfehler dar
- Für die Variabiliät zwischen den Gruppen spielt nach unseren obigen Überlegungen zum einen die Schwankung innerhalb der Gruppen eine Rolle; denn je stärker die Werte innerhalb einer Gruppe varriieren, desto größer wird auch die Differenz bzw. Varianz der Mittelwerte. Zum anderen aber unterscheiden sich die Gruppenmittelwerte aber natürlich besonders stark, wenn die wahren Populationsmittelwerte sich tatsächlich unterscheiden, der Faktor also einen Effekt auf die abhängige Variable hat.
Man kann also den \(F\)-Bruch auch verstehen als \[ F=\frac{\text{Variabilität zwischen den Gruppen}}{\text{Variabilität innerhalb der Gruppen}} = \frac{\text{Effekte+Messfehler+individuelle Unterschiede}}{\text{Messfehler+individuelle Unterschiede}}\] und aus dieser Schreibweise folgen zwei wichtige Aspekte:
- Wenn es keinen Effekt des Faktors gibt (also Effekt = 0), und damit keine Unterschiede zwischen den Populationsmittelwerten, dann nimmt \(F\) in etwa einen Wert von 1 an: \(F \approx 1\)
- Bei Anwesenheit von Effekten, also wenn \(F>1\) wird, wird der \(F\)-Bruch umso größer, je “stärker” die Effekte sind.
Wir werden nun zunächst eine Varaianzanalyse berechnen, indem wir zeigen, wie die beiden Variabilitäten quantifiziert werden, wie die entsprechende Zufallsvariable verteilt ist und wie wir zu einer Entscheidung zwischen den beiden Hypothesen kommen.
2.2 Rechnerische Durchführung
Die rechnerische Durchführung einer Varianzanalyse erarbeiten und zeigen wir anhand der Daten aus dem Eingangsbeispiel:
Im allgemeinen Fall gehen wir jedoch davon aus, dass es \(J\) Gruppen/Bedingungen mit jeweils \(n_j\) Versuchspersonen gibt. Wir werden die benötigten Formeln daher immer direkt für den allgemeinen Fall notieren, aber jeweils mit den Daten des Beispiels, also mit \(J=3\), rechnen.
Zunächst berechnen sich die Gruppenmittelwerte und -varianzen so wie üblich, wir müssen lediglich einige Indizes einführen, um klar zu machen, zu welcher Gruppe der Wert gehört:
- Gruppenmittelwert: \[M_j=\frac{1}{n_j}\sum_{i=1}^{n_j}y_{ij}\hspace{1cm}\forall j\in\{1,\ldots, J\}\]
- Gruppenvarianz: \[S_j^2=\frac{1}{n_j}\sum_{i=1}^{n_j}(y_{ij}-M_j)^2\hspace{1cm}\forall j\in\{1,\ldots, J\}\]
Schließlich werden wir noch den Gesamtmittelwert (engl.: grand mean) benötigen. Mit \(N\) meinen wir die Gesamtzahl der Beobachtungen/Versuchspersonen, also: \[N=\sum_{i=1}^Jn_j\] Dann kann der Gesamtmittelwert \(M\) berechnet werden als \[M=\frac{1}{N}\sum_{j=1}^J\sum_{i=1}^{n_j}y_{ij}\] bzw. als gewichtetes Mittel der Gruppenmittelwerte: \[\sum_{j=1}^J\frac{n_j}{N}M_j\]
Für die Beispieldaten ergibt sich: \[ \sum_{j=1}^J\frac{n_j}{N}M_j=\frac{10}{30}\cdot 13.1+\frac{10}{30}\cdot 15.7+\frac{10}{30}\cdot 18.1=15.63 \]
2.2.1 Quadratsummenzerlegung
Im Zuge der einfachen, linearen Regression (Teil 6 vo Statistik I hatten wir folgenden Zusammenhang hergeleitet: \[S_Y^2=S_{\hat{Y}}^2+S_E^2\] Dies hatten wir als ein Beispiel für sog. Varianzaufteilung bezeichnet: Die Varianz der \(Y\)-Werte setzt sich additiv zusammen aus der Varianz der vorhergesagten/erklärten \(\hat{Y}\)-Werte sowie der Varianz der Residuen \(E\) (dem Messfehler/dem nicht-erklärten Anteil).
Ganz ähnliches machen wir im Fall einer Varianzanalyse nun: Wir teilen die Gesamtvariabilität der Daten in einen systematischen/erklärten Anteil und einen Messfehleranteil auf. Diese Variabilitäten werden in der Regel Quadratsummen (QS) bzw. Sums of Squares (SS) genannt und wir berechnen hier zunächst drei solche Quadratsummen. Jede einzelne quantifiziert einen bestimmt Aneil der Variabilität der Daten:
- Die Gesamt-QS (Sums of Squares total) \(SS_{\text{tot}}\) quantifiziert die Variabilität aller Daten ohne Berücksichtigung der Gruppen/Bedingungen.
- Die QS zwischen den Gruppen (Sums of Squares between) \(SS_A\) quantifiziert die Variabiltät der Gruppen-/Bedingungsmittelwerte. Manchmal wird diese QS auch als \(SS_{\text{between}}\) bezeichnet. Da wir in Teil 4 auch den Fall mehrerer Faktoren betrachten, nutzen wir hier direkt die Notation \(SS_A\) (für die Quadratsumme zwischen den Gruppenmittelwerten des Faktors \(A\) und später dann entsprechend \(SS_B\) für den Faktor \(B\)).
- Die QS innerhalb der Gruppen (Sums of Squares within) \(SS_w\) quantfiziert schließlich die Variabilitäten innerhalb der Gruppen, wie also de einzelnen Werte um den Gruppenmittelwert herum variieren.
Im Folgenden führen wir immer die allgemeinen Formeln ein und zeigen dann anhand des Eingangsbeispiels eine konkrete Berechnung. Es ist hilfreic, diese Berechnungen selber komplett nachzuvollziehen.
Gesamt-QS (Sums of Squares total): \(SS_{\text{tot}}\) \[SS_{\text{tot}}=\sum_{j=1}^J\sum_{i=1}^{n_j}(y_{ij}-M)^2\] Bezogen auf das Beispiel ergibt sich: \[\begin{equation*} \begin{aligned} SS_{\text{tot}} &=\sum_{j=1}^3\sum_{i=1}^{n_j}(y_{ij}-15.63)^2\\ &=(17-15.63)^2+(12-15.63)^2+\ldots+(20-15.63)^2+(18-15.63)^2\\ &=264.97 \end{aligned} \end{equation*}\] Tatsächlich ist die \(SS_\text{tot}\) fast die (“normale”) Varianz aller Daten. Aufsummiert werden hier die quadrierten Abweichungen der Werte vom Gesamtmittelwert. Um die Varianz zu erhalten müsste schließlich noch die Division durch \(N\) erfolgen. Aus diesem Grund könnte man auch sagen, die \(SS_\text{tot}\) ist die \(N\)-fache Varianz der Daten.
QS zwischen den Gruppen (Sums of Squares between): \(SS_A\) \[SS_A=\sum_{j=1}^Jn_j(M_j-M)^2\] Bezogen auf das Beispiel ergibt sich:
\[\begin{equation*} \begin{aligned} SS_A&=\sum_{j=1}^3n_j(M_j-15.63)^2\\ &=10\cdot(13.1-15.63)^2+10\cdot(15.7-15.63)^2+10\cdot(18.1-15.63)^2\\ & = 125.07 \end{aligned} \end{equation*}\] Hier werden i.W. die quadrierten Abweichungen der Gruppen-/Bedingungsmittelwerte vom Gesamtmittlwert quadriert und aufsummiert. Da dieser Term für jede Person berücksichtigt wird, wird er innerhalb der Summe mit den jeweiligen Stichprobengrößen \(n_j\) gewichtet.
QS innerhalb der Gruppen (Sums of Squares within): \(SS_w\) \[ SS_w=\sum_{j=1}^J\sum_{i=1}^{n_j}(y_{ij}-M_j)^2=\sum_{j=1}^Jn_jS_j^2\] Bezogen auf das Beispiel ergibt sich:
\[\begin{equation*} SS_w=\sum_{j=1}^3n_jS_j^2=10\cdot 7.69+10\cdot 2.01 + 10\cdot 4.29=139.90 \end{equation*}\] Hier werden für jede Gruppe/Bedingung separat die quadrierten Abweichungen der Werte von den Gruppen-/Bedingungsmittelwerten aufsummiert. Würde für jede Gruppe/Bedingung noch durch \(n_j\) dividiert werden, würden die Varianzen \(S_j^2\) aufsummiert werden. Da die Division aber nicht stattfindet, werden also die mit \(n_j\) multiplizierten Varianzen aufsummiert.
Insgesamt ergibt sich mit diesen Formeln für die Quantifizierung der Variabilitäten die folgende wichtige Eigenschaft der “Quadratsummenzerlegung” (Dtreng genommen kann man zwischen einer Varianz- und einer Quadratsummenzerlegung unterscheiden. Das Prinzip ist allerdings das gleiche und wir unterscheiden hier nicht.): \[\boxed{SS_{\text{tot}} = SS_A+SS_w}\] Bezogen auf das Beispiel also \[264.97=125.07+139.90\]
2.2.2 Mittlere Quadratsummen
In die gerade berechneten Quadratsummen gehen jeweils unterschiedlich viele Werte ein. So gehen z.B. in die \(SS_A\) nur die jeweiligen Mittelwerte ein, während in die \(SS_w\) auch Werte einzelner Versuchspersonen eingehen. Die Quadratsummen sind daher nicht direkt miteinander vergleichbar. Allerdings unterscheiden sich die Freiheitsgrade der jeweiligen Quadratsummen und um sie besser vergleichbar zu machen, werden die Quadratsummen jeweils an ihren Freiheitsgraden relativiert und dann Mittlere Quadratsummen (MS) bzw. Mean Squares (MS) genannt:
\(SS_A\) hat \(J-1\) Freiheitsgrade: \[MS_A=\frac{SS_A}{J-1}\] Bezogen auf das Beispiel ergibt sich: \[\begin{equation*} MS_A=\frac{125.07}{3-1}=62.53 \end{equation*}\]
\(SS_w\) hat \(N-J\) Freiheitsgrade: \[MS_w=\frac{SS_w}{N-J}\] Bezogen auf das Beispiel ergibt sich: \[\begin{equation*} MS_w=\frac{139.90}{30-3}=5.18 \end{equation*}\] \(MS_w\) wird oft auch Fehlerterm oder Mean Square Error (MSE) genannt.
Was haben wir bisher mit der Berechnung von \(MS_A\) und \(MS_w\) nun gewonnen? Wir sollten für den nächsten Abschnitt festhalten, dass \(MS_A\) die Variabilität der Gruppenmittelwerte quantifiziert und dass \(MS_w\) die Variabilität innerhalb der Gruppen quantifiziert. Diese wiederum waren weiter oben als Zähler und Nenner des \(F\)-Bruchs erarbeitet worden. Dass wir nicht einfach Varianzen berechnen hat den Grund, dass wir dann Probleme bei der Bestimmung der Dichtefunktion der Zufallsvariablen des \(F\)-Bruchs bekommen würden.
2.2.3 der \(F\)-Bruch und die \(F\)-Verteilung
Der nächste Schritt ist nun, dass wir aus den Mittleren Quadratsummen den empirischen \(F\)-Bruch berechnen können: \[F=\frac{MS_A}{MS_w}\] Bezogen auf das Beispiel ergibt sich: \[\begin{equation*} F=\frac{62.53}{5.18}=12.07 \end{equation*}\] Ähnlich wie beim \(t\)-Bruch müssen wir uns nun noch fragen: Wie ist eine Zufallsvariable \(\mathbf{F}\) verteilt, die jeder Kombination von \(J\) Stichproben den empirischen \(F\)-Bruch zuweist? Hierzu schauen wir zunächst wieder eine Simulation an, die konzeptuell an das Beispiel anschließt:
- Die drei Stichproben werden hierbei aus der gleichen normalverteilten Population gezogen (\(\mu=15\), \(\sigma^2=25\) bzw. \(\sigma = 5\)), d.h. die Gültigkeit der Nullhypothese \(H_0\) wird angenommen.
- Diese Ziehung wiederholen wir 2,000 Mal, berechnen jeweils den \(F-\)Wert und stellen die Verteilung der empirischen \(F\)-Werte als Histogramm grafisch dar. (Die durchgezogene schwarze Linie stellt die tatsächliche theoretische Verteilung dar. Anmerkung: Die y-Achse repräsentiert keine richtigen relativen Häufigkeiten, sondern ist so skaliert, dass die Summe der Flächen aller Balken 1 ergibt; analog zur Eigenschaft, dass die Fläche unter einer Dichtefunktion 1 ergibt)
library(HistogramTools)
set.seed(1)
f_values <- NULL # Container der F-Werte
wie_oft <- 20000 # Wie oft soll die "Studie" durchgeführt werden
for (i in 1:2000) {
x1 <- rnorm(10, mean = 15, sd = 5) # Werte der Stichprobe 1 (Trainingstag 1)
x2 <- rnorm(10, mean = 15, sd = 5) # Werte der Stichprobe 2 (Trainingstag 2)
x3 <- rnorm(10, mean = 15, sd = 5) # Werte der Stichprobe 3 (Trainingstag 3)
av <- c(x1,x2,x3) # alle AV Werte
training <- rep(1:3, each = 10) # Variable
daten <- as.data.frame(cbind(training,av)) # Data.Frane der Werte
# ANOVA in R rechnen (kommt weiter unten im Text)
daten$training <- as.factor(daten$training)
aov_ergebnis <- aov(av ~ training,
data = daten)
summary_aov <- summary(aov_ergebnis)
f_value <- summary_aov[[1]][1,4] # F Wert extrahieren
f_values <- c(f_values, f_value) # und abspeichern
}
# Ergebnis plotten
par(mar = c(5,5,0,1), mfrow = c(1,1))
hist(f_values,
breaks = 20,
freq = FALSE,
col = "lightgray",
axes = FALSE,
ylim = c(0,0.8),
xlab = expression(italic(F)),
ylab = "relative Häufigkeit/Dichte",
main = "",
xlim = c(0, 5),
cex.lab = 1.5)
axis(1, # x-Achse
las = 1,
at = seq(0, 5, 1),
cex.axis = 1.2)
axis(2, # y-Achse
las = 1,
at = seq(0, 0.8, 0.1),
cex.axis = 1.2)
abline(h = 0)
x <- seq(0, 5, 0.1)
curve(df(x, 2, 27),
add = TRUE,
lwd = 2)
Zugegebenermaßen ist die Form der \(F\)-Verteilung bei 3 Gruppen etwas besonders (siehe dazu auch die Abbildung mit den echten \(F\)-Verteilungen weiter unten). Zur Illustration wiederholen wir die Simulation daher noch einmal, dieses Mal allerdings mit 5 Gruppen:
Wir können hier einige Merkmale der resultierenden Verteilung (bzw. der relativen Häufigkeiten) erkennen:
- Die Verteilung ist nicht symmetrisch um einen bestimmten Wert (wie z.B. die Normalverteilung oder eine zentrale \(t\)-Verteilung es sind).
- Negative \(F-\)Werte kommmen nicht vor.
- Auf den ersten Blick zeigt die Verteilung eine gewisse Ähnlichkeit zur \(\chi^2\)-Verteilung; sie ist aber keine.
Tatsächlich kann man zeigen, dass bei Gültigkeit der \(H_0\) der \(F\)-Bruch (zentral) \(F\)-verteilt ist. Während eine zentrale \(t\)-Verteilung nur durch einen Parameter – ihre Freiheitsgrade – gekennzeichnet ist, ist die zentrale \(F\)-Verteilung durch zwei Parameter charakterisiert:
- Die Zählerfreiheitsgrade entsprechenden den Freiheitsgraden von \(SS_A\), also \(J-1\).
- Die Nennerfreiheitsgrade entsprechen den Freiheitsgraden von \(SS_w\), also \(N-J\).
Insgesamt folgt, dass bei Gültigkeit der Nullhypothese \(\mathbf{F}\) \(F\)-verteilt mit \(J-1\) und \(N-J\) Freiheitsgraden ist: \[\mathbf{F}\overset{H_0}{\sim}F_{J-1,N-J}\]
Die folgende Abbildung visualisiert Beispiele für \(F\)-Verteilungen mit \(m\) Zähler- und \(n\) Nennerfreiheitsgraden. Auffällig ist dabei, dass die Form für \(m=2\) Zählerfreiheitsgrade anders ist, als für \(m>2\) Zählerfreiheitsgrade:
x <- seq(0, 10, length = 500) # x-Werte zum zeichnen
y1 <- df(x, 2, 2) # Dichtewerte für m = 2 und n = 2
y2 <- df(x, 5, 2) # Dichtewerte für m = 5 und n = 2
y3 <- df(x, 5, 10) # Dichtewerte für m = 5 und n = 10
y4 <- df(x, 10, 10) # Dichtewerte für m = 10 und n = 10
# Plotten aller Dichtefunktionen
par(mar = c(5, 5, 2, 1), mfrow = c(1, 1))
# Grundplot: mit einer Dichtefunktion
plot(x, y1, # 2,2
xlim = c(0, 3), ylim = c(0, 1),
type = "l", xlab = expression(italic(F)),
ylab = "Dichte", axes = FALSE, cex.lab = 1.5
)
abline(h = 0, col = "black")
axis(1, las = 1, at = c(0, 1, 2, 3), cex.axis = 1.3) # x-Achse
axis(2, las = 1, at = seq(0, 1, 0.2), cex.axis = 1.3) # y-Achse
# Hinzufügen der weiteren Dichtefunktionen
lines(x, y2, lwd = 2, col = "red") # 5,2
lines(x, y3, lwd = 2, col = "blue") # 5,10
lines(x, y4, lwd = 2, col = "orange") # 10,10
# Legende hinzufügen
legend("topright", bty = "n", lty = 1,
col = c("black", "red", "blue", "orange"), lwd = 2,
legend = c(expression(italic(F)(2, 2)),
expression(italic(F)(5, 2)),
expression(italic(F)(5, 10)),
expression(italic(F)(10, 10))),
title = expression(italic(F)(italic(m), italic(n))))
Ist eine Zufallsvariable \(\mathbf{F}\) \(F\)-verteilt mit \(m\) und \(n\) Freiheitsgraden, also \(\mathbf{F}\sim F_{m,n}\), dann sind ihr Erwartungswert und ihre Varianz: \[E(\mathbf{F})=\frac{n}{n-2}\hspace{1cm}\text{für }n>2\] \[Var(\mathbf{F})=\frac{2n^2(m+n-2)}{m(n-2)^2(n-4)}\hspace{1cm}\text{für }n>4\]
2.2.4 die Entscheidungsregel(n)
Beim \(t\)-Test haben wir zwei Arten von Entscheidungsregeln zwischen \(H_0\) und \(H_1\) kennengelernt:
- Entscheidungen mit kritischen \(t\)-Werten
- Entscheidungen auf Basis von \(p\)-Werten
Beide Entscheidungsregeln wenden wir nun ganz analog auch für die Varianzanalyse mit \(F\) als Prüfgröße an. Zur Erinnerung: Beim (gerichteten) \(t\)-Test haben wir einen Wert \(t_{krit}\) gesucht, rechts von dem noch 5% der Fläche unter der Dichteverteilung liegen (bei einem festgelegten \(\alpha=0.05\)):
Ganz ähnlich gilt hier nun: Beim \(F\)-Test der Varianzanalyse suchen wir einen Wert \(F_{krit}\), rechts von dem noch 5% der Fläche unter der Dichteverteilung liegen (bei einem festgelegten \(\alpha=0.05\)):
Gesucht ist also \(F_{krit}\), d.h. der Wert rechts von dem 5% und links von dem 95% der Fläche liegen, so dass gilt. \[ \int_{F_{krit}}^{+\infty}f(x)dx=0.05 \text{, wobei }f(x)\small{\text{ die Dichtefunktion der }F\text{-Verteilung ist}} \] Ganz analog zum Vorgehen beim \(t\)-Test und der \(t\)-Verteilung ist dieser Wert nun das \((1-\alpha)\)-Quantil der F-Verteilung mit \(m\) und \(n\) Freiheitsgraden: \[F_{krit}=F_{m,n;1-\alpha}\]
Das wiederum heißt:
- Wenn die \(H_0\) gilt, dann tritt \(F\geq F_{krit}\) mit einem relativen Anteil von \(\alpha=0.05\), also in 5% aller Fälle auf.
- In diesem Fall ist der F-Wert der einen Studie so unwahrscheinlich, dass wir an der Annahme der \(H_0\) zweifeln.
- In der Folge entscheiden wir uns stattdessen, von der Gültigkeit der \(H_1\) auszugehen: der F-Test/die Varianzanalyse liefern ein signifikantes Ergebnis.
Und auch ganz ähnlich wie beim \(t\)-Test (vgl. Teil 12) lautet dann die Entscheidungsregel 1 auf Basis kritischer Werte:
Entscheidungsregel 1 t-Test (gerichtete \(H_1:\mu_A>\mu_B\)): “Wenn \(t\geq t_{krit}\) ist, dann tritt der \(t\)-Wert (bzw. ein noch größerer) so selten auf, wenn die \(H_0\) gelten würde, dass wir an dieser Annahme Zweifel haben. Wir entscheiden uns daher für die \(H_1\).”
Entscheidungsregel 1 F-Test/Varianzanalyse: “Wenn \(F\geq F_{krit}\) ist, dann tritt der \(F\)-Wert (bzw. ein noch größerer) so selten auf, wenn die \(H_0\) gelten würde, dass wir an dieser Annahme Zweifel haben. Wir entscheiden uns daher für die \(H_1\).”
Den kritischen \(F\)-Wert können wir
– genau wie den kritischen \(t\)-Wert –
auf verschiedene Arten bestimmen. Prizipiell benötigen wir den Wert
\(F\), an dem die Verteilungsfunktion
den Wert \(F(F) = 0.95\) annimmt. Diese
Werte sind in Tabellen zu finden (auch bei Wikipedia)
oder können sehr leicht mit der Funktion qf() von R
bestimmt werden (hier für 3 Zähler- und 30 Nennerfreiheitsgrade):
# qf = Quantile der F-Verteilung
qf(p = 0.95, # welches Quantil?
df1 = 3, # Zählerfreiheitsgrade
df2 = 30) # Nennerfreiheitsgrade## [1] 2.922277Für die 2. Entscheidungsregel haben wir beim \(t\)-Test bereits den \(p\)-Wert berechnet. Dieser war die Wahrscheinlichkeit dafür, dass ein empirischer \(t\)-Wert oder ein größerer \(t\)-Wert auftritt, wenn die Nullhypothese gilt: \[p=\int^{+\infty}_{t_{\text{empirisch}}}f(x)dx\] wobei \(f(x)\) die Dichtefunktion der t-Verteilung ist.
Die gleiche Logik gilt nun auch beim \(F\)-Test/der Varianzanalyse: Der \(p\)-Wert ist die Wahrscheinlichkeit dafür, dass ein empirischer \(F\)-Wert oder ein größerer \(F\)-Wert auftritt, wenn die Nullhypothese gilt: \[p=\int^{+\infty}_{F_{\text{empirisch}}}f(x)dx\] wobei \(f(x)\) die Dichtefunktion der F-Verteilung ist.
Für einen bestimmten empirischen \(F\)-Wert können wir auch leicht eine
Berechnung des \(p\)-Werts mit der
Funktion pf() in R durchführen (für ein Beispiel \(F_{\text{empirisch}} = 3.20\) und 3 Zähler-
sowie 30 Nennerfreiheitsgraden):
# pf = Verteilungsfunktion der F-Verteilung
pf(q = 3.20, # berechnet die Fläche bis zu 3.20 unter der F-Verteilung
df1 = 3,
df2 = 30) ## [1] 0.96264231 - pf(q = 3.20, # berechnet die Fläche ab 3.20 unter der F-Verteilung,
df1 = 3, # also p
df2 = 30) ## [1] 0.03735765Durch den Vergleich von \(p\) mit dem vorab festgelegten \(\alpha\) können wir dann auch wieder ganz analog eine Entscheidungsregel 2 aufstellen (vgl. Teil 12:
Entscheidungsregel 2 t-Test: “Wenn \(p\leq\alpha\) ist, dann tritt der \(t\)-Wert bzw. ein noch größerer \(t\)-Wert so selten auf, wenn die \(H_0\) gelten würde, dass wir an dieser Annahme Zweifel haben. Wir entscheiden uns daher für die \(H_1\).”
Entscheidungsregel 2 F-Test/Varianzanalyse: “Wenn \(p\leq\alpha\) ist, dann tritt der \(F\)-Wert bzw. ein noch größerer \(F\)-Wert so selten auf, wenn die \(H_0\) gelten würde, dass wir an dieser Annahme Zweifel haben. Wir entscheiden uns daher für die \(H_1\).”
Für unser Eingangsbeispiel ergibt sich folgende Entscheidung (2 Zähler- und 27 Nennerfreiheitsgrade): \(F=12.07\) und \(\alpha=0.05\)
- Berechnung des kritischen \(F\)-Werts:
qf(p = 0.95, # F_krit für alpha = 0.05
df1 = 2,
df2 = 27) ## [1] 3.354131- Berechnung des \(p\)-Werts:
1 - pf(q = 12.07, # p
df1 = 2,
df2 = 27) ## [1] 0.0001799455In beiden Fällen fällt die Entscheidung also zugunsten von \(H_1\).
Die gleichen Anmerkungen wie schon bei den Entscheidungen im Falle eines \(t\)-Tests gelten auch hier:
- Beide Entscheidungsregeln kommen zur gleichen Entscheidung.
- Bei einer Entscheidung zugunsten der \(H_1\), redet man auch von einem signifikanten Ergebnis.
Und schließlich hat der p-Wert die gleiche Bedeutung wie schon beim \(t\)-Test:
- Er wurde bestimmt unter Annahme der Gültigkeit der \(H_0\): Der \(p\)-Wert ist daher eine bedingte Wahrscheinlichkeit über Daten, gegeben \(H_0\) gilt \[p=p(\text{Daten}|H_0)\]
- Der p-Wert sagt etwas über Wahrscheinlichkeiten von Daten aus….
- …aber nichts über die Wahrscheinlichkeit der Nullhypothese.
- Die Gültigkeit der Nullhypothese wird angenommen.
Bevor wir dann wieder etwas theoretischer werden, betrachten wir
abschließend noch, wie die Ergebnisse einer Varianzanalyse dargestellt
werden können. Zum einen werden Varianzanalysen manchmal vollständig in
Form einer Tabelle berichtet (Effektstärken führen wir dann in Teil
3 ein):
Sollen die Ergebnisse hingegen in Textform berichtet werden, dann wäre eine passende Formulierung z.B.: “Die mittlere Punktzahl im Sprachtest ist in Abb. XX dargestellt und sie nimmt mit der Anzahl der Trainingstage zu. Die Varianzanalyse ergab einen signifikanten Effekt der Gruppenzugehörigkeit/der Anzahl Trainingstage, \(F(2,27)=12.07\), \(p <.001\).”
2.2.5 Voraussetzungsverletzungen
Genau wie auch beim \(t\)-Test haben wir am Anfang der Varianzanalyse Voraussetzungen formuliert. Sind diese Voraussetzungen verletzt, so kann die Verteilung von \(\mathbf{F}\) nicht exakt bestimmt werden. Es lässt sich hier festhalten:
- Normalverteiltheitsannahme: Bei Verletzungen ist die Varianzanalyse dennoch recht robust und kommt zu weitgehend validen Entscheidungen, insbesondere bei \(n_j>25\); ansonsten helfen ggf.geeignete Transformationen der Daten oder die entsprechenden non-parametrischen Alternativenm hier der Kruskal-Wallis-Test.
- Wenn die Daten nur Ordinalskalenniveau aufweisen, ist der non-parametrische Kruskal-Wallis-Test das Verfahren der Wahl.
- Varianzhomogenität haben wir bereits beim \(t\)-Test getestet und ganz ähnlich können wir auch im Falle von \(J> 2\) den Levene-Test nutzen:
library(car)
daten_training$Trainingstage <- as.factor(daten_training$Trainingstage)
leveneTest(Sprachtest ~ Trainingstage,
data = daten_training,
center = "mean")Für die Beispieldaten ist der Levene-Test nicht signifikant, wir gehen daher weiterhin von Varianzhomogenität aus. Verletzungen der Varianzhomogenität sind insbesondere dann problematisch, wenn die Stichproben verschieden groß sind. Dies sollte daher idealerweise bereits bei Planung einer Studie berücksichtigt werden.
Im Fall des \(t\)-Tests für zwei unabhängige Stichproben hatten wir in Teil 12 von Statistik I bereits den Welch-Test kennen gelernt. Eine Erweiterung dieses Tests auf mehr als zwei unabhängige Gruppen gibt es auch als Möglichkeit des Umgangs mit Varianzheterogenität für die Varianzanalyse. Eine weitere Alternative stellt der sog. Brown-Forsythe-Test dar. Beide Varianten berechnen den \(F\)-Wert sowie die Freiheitsgrade auf eine andere Art und Weise und kompensieren so für die entstandene Liberalität des \(F\)-Tests (für Details siehe Eid et al., 2010).
Alternativ werden Daten manchmal transformiert, um Varianzen homogener zu bekommen. Dies ist insbesondere dann sinnvoll, wenn Mittelwerte und Varianzen proportional zueinander sind. Die genaue Art der Transformation hängt dabei von der Art der Daten ab (für Details, siehe z.B. Keppel et al., 2004):
- Häufigkeiten: \(Y'=\sqrt{Y+\frac{1}{2}}\)
- (Reaktions-)Zeiten: \(Y'=log(Y+1)\)
- Anteile (z.B. bei Fehlern): \(2\cdot \arcsin(\sqrt{Y})\)
2.2.6 Warum ist der \(F\)-Bruch \(F\)-verteilt?
Wir haben weiter oben einfach festgestellt, dass \(\mathbf{F}\) \(F\)-verteilt ist. Die komplette Begründung ist etwas aufwändiger, aber wir wollen hier kurz die grundlegende Idee wieder demonstieren.
Zunächst halten wir dazu, dass sich eine \(\chi^2\)-Verteilung aus der Standardnormalverteilung herleiten lässt: Die Summe \(n\)-vieler quadrierter standardnormalverteilter Zufallsvariablen ist \(\chi^2\)-verteilt mit \(n\) Freiheitsgraden. Seien also \(\mathbf{X_1},\ldots,\mathbf{X_n}\) unabhängige und standardnormalverteilte Zufallsvariablen, dann gilt: \[\sum_{i=1}^n\mathbf{X_i}^2 \sim \chi^2_n\] Weiter lässt sich eine F-Verteilung mit \(m\) Zählerfreiheitsgraden und \(n\) Nennenfreiheitsgraden formal definieren als: \[\mathbf{F}\equiv \frac{\frac{\mathbf{U}}{m}}{\frac{\mathbf{V}}{n}}\] wobei \(\mathbf{U}\sim\chi^2_m\) und \(\mathbf{V}\sim\chi^2_n\) seien. \(\mathbf{U}\) und \(\mathbf{V}\) seien ferner unabhängig voneinander.
Etwas komplizierter herzuleiten sind die beiden folgenden Behauptungen. Für die Verteilungen von \(SS_A\) und \(SS_w\) unter \(H_0\) kann man zeigen:
- \(\frac{SS_A}{\sigma^2}\) ist \(\chi^2\)-verteilt: \[\frac{SS_A}{\sigma^2}\sim\chi^2_{J-1}\]
- \(\frac{SS_w}{\sigma^2}\) ist \(\chi^2\)-verteilt: \[\frac{SS_w}{\sigma^2}\sim\chi^2_{N-J}\]
Da überdies \(SS_A\) und \(SS_w\) unabhängig sind, folgt die oben behauptete Verteilung für \(\mathbf{F}\): \[\frac{\frac{SS_A/\sigma^2}{J-1}}{\frac{SS_w/\sigma^2}{N-J}}=\frac{\frac{SS_A}{J-1}}{\frac{SS_w}{N-J}}=\frac{MS_A}{MS_w}= F\sim F_{J-1,N-J}\]
2.3 Berechnung mit R
Im letzten Abschnitt von Teil 2 zeigen wir, wie eine Varianzanalyse mit R berechnet wird. Hierbei verwenden wir wiederum die Daten des Eingangsbeispiels:
head(daten_training)Grundsätzlich gibt es in R viele Möglichkeiten zur Berechnung (einfaktorieller) Varianzanalysen:
- Die Funktion
aov()ist die R-eigene Funktion, mit der varianzanalytische Modelle mit der Formelsprache modelliert werden. - Die Funktion
oneway.test()bietet sich dazu für einfaktorielle Varianzanalysen an, da sie bei Verletzungen der Varianzhomogenität benutzt werden um \(F\)-Wert und Freiheitsgrade entsprechend anzupassen (vgl. Welch-Test als Alternative zum \(t\)-Test; Teil 12 in Statistik I). - Es gibt Pakete, die sog. Wrapperfunktionen bereitstellen, mit denen
Varianzanalysen einfach formuliert werden können:
- Wir werden das Paket
ez(Lawrence, 2016) und die darin enthaltene FunktionezANOVA()nutzen. - Eine Alternative ist das Paket
afex(Singmann et al., 2023).
- Wir werden das Paket
Wir beschränken uns hier auf die Funktion aov() sowie
die Funktionen des Paketes ez. Beispiele zur Verwendung von
oneway.text() sowie der Funktionen des Paketes
afex, wird es in Zukunft als Ergänzung geben.
2.3.1 Funktion
aov()
Die Funktion aov() benutzt die Modellsprache um die
abhängige Variable durch den Faktor zu modellieren. Eine ganz ähnliche
Variante hatten wir bereits beim \(t\)-Test kennengelernt. Das Ergebnis können
wir dann mit der Funktion summary() ausgeben (vergleichen
Sie die Ausgabe mit der händischen Berechnung weiter oben):
daten_training$Trainingstage <- as.factor(daten_training$Trainingstage) # faktorisieren
aov_ergebnis <- aov(Sprachtest ~ Trainingstage, # AV ~ Faktor
data = daten_training) # welche Daten?
summary(aov_ergebnis) # Ergebnis ausgeben## Df Sum Sq Mean Sq F value Pr(>F)
## Trainingstage 2 125.1 62.53 12.07 0.00018 ***
## Residuals 27 139.9 5.18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Mit dem aov()-Objekt können wir auch relativ leicht noch
deskriptive Darstellungen der Mittelwerte der Gruppen sowie des
Gesamtmittelwerts erhalten:
model.tables(aov_ergebnis,
type = "means") # Mittelwerte berechnen## Tables of means
## Grand mean
##
## 15.63333
##
## Trainingstage
## Trainingstage
## 1 2 3
## 13.1 15.7 18.1Interessieren uns stattdessen die Abweichungen der Mittelwerte vom Gesamtmittelwert (darauf werden wir in Teil 3 Effektgrößen aufbauen), so können wir diese ebenfalls leicht berechnen:
model.tables(aov_ergebnis,
type = "effects") # Effekte berechnen## Tables of effects
##
## Trainingstage
## Trainingstage
## 1 2 3
## -2.5333 0.0667 2.46672.3.2 Paket
ez
Die Funktion aov() ist relativ flexibel, erfordert aber
bei komplexeren Varianzanalysen ein gewisses Maß an Kenntnis zur
korrekten Verwendung. Um dies einfacher zu gestalten, stellt
ez mit der Funktion ezANOVA() eine sog.
“Wrapperfunktion” zur Verfügung, die wiederum eine einfache
Spezifizierung von Varianzanalysen erlaubt. Wir werden auch in den
folgenden Kapiteln mit dem Paket ez arbeiten. Die hier
betrachtete einfaktorielle (between-subject) Varianzanalyse kann wie
folgt berechnet werden:
suppressPackageStartupMessages(library(ez)) # library(ez) reicht auch...
daten_training$VP <- as.factor(daten_training$VP) # faktorisieren
daten_training$Trainingstage <- as.factor(daten_training$Trainingstage) # faktorisieren
ez_ergebnis <- ezANOVA(data = daten_training, # welche Daten?
wid = .(VP), # 'within-id' = Versuchsperson
dv = .(Sprachtest), # "dependent variable"
between = .(Trainingstage), # between-Faktor(en)
detailed = TRUE) # zusätzliche Info für schoRsch## Coefficient covariances computed by hccm()ez_ergebnis # Ergebnis ausgeben## $ANOVA
## Effect DFn DFd SSn SSd F p p<.05 ges
## 1 Trainingstage 2 27 125.0667 139.9 12.06862 0.0001800766 * 0.4720091
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 2 27 7.466667 53.5 1.884112 0.1714068Zum Aufruf und den Ergebnissen fügen wir drei Anmerkungen hinzu:
- Wird nur eine Variable z.B. für
between = ...angegeben, kann die.()-Konstruktion auch weggelassen werden (d.h.wid = VPoderbetween = Trainingstagewürde hier auch reichen). - Der Levene-Test entspricht dem Ergebnis der Funktion
leveneTest(), wenncenter = "median"verwendet wird. gessteht für “generalized \(\eta^2\)”, eine Effektstärke, die im Kontext der Varianzanalyse manchmal verwendet wird.
Eine kompaktere Formatierung erhalten wir mit der Funktion
anova_out() aus dem Paket schoRsch:
library(schoRsch)
anova_out(ez_ergebnis) # formatierte Ausgabe## $`--- ANOVA RESULTS ------------------------------------`
## Effect MSE df1 df2 F p petasq getasq
## 1 Trainingstage 5.181481 2 27 12.07 0.000 0.47 0.47
##
## $`--- SPHERICITY TESTS ------------------------------------`
## [1] "N/A"
##
## $`--- FORMATTED RESULTS ------------------------------------`
## Effect Text
## 1 Trainingstage F(2, 27) = 12.07, p < .001, np2 = .47
##
## $`NOTE:`
## [1] "Reporting unadjusted p-values."Auch hier Zwei Anmerkungen:
- “Sphericity Tests” spielt erst bei Varianzanalysen mit Messwiederholung (siehe Teil 5) eine Rolle.
np2steht für “partial \(\eta^2\)” und ist ebenfalls eine Effektstärke, die wir in Teil 3 kennenlernen werden.
Das Paket ez() stellt noch eine Reihe weiterer
nützlicher Funktionen zur Verfügung, die den gleichen generellen Aufbau
wie ezANOVA() haben. Mit ezStats() können wir
z.B. schnell deskriptive Statistiken wie den Mittelwert berechnen:
stats <- ezStats(data = daten_training, # welche Daten?
wid = .(VP), # 'within-id' = Versuchsperson
dv = .(Sprachtest), # dependent variable
between = .(Trainingstage)) # between-Faktor(en)## Coefficient covariances computed by hccm()statsstats ist dabei ein DataFrame, auf dessen Spalten auch
einzeln zugegriffen werden kann:
stats$Mean## [1] 13.1 15.7 18.1Dies ist z.B. sehr nützlich für die flexible Gestaltung von Abbildungen oder wenn die Mittelwerte für weitere Analysen benötigt werden.
Die Funktion ezPlot() erlaubt zudem die schnelle
grafische Darstellung von Daten:
ezPlot(data = daten_training, # welche Daten?
wid = .(VP), # 'within-id' = Versuchsperson
dv = .(Sprachtest), # dependent variable
between = .(Trainingstage), # between-Faktor(en)
x = .(Trainingstage)) # welche Variable auf die x-Achse?## Coefficient covariances computed by hccm()
Die hier generierte Abbildung ist zwar nicht sonderlich schön, aber
bietet einen schnellen Überblick über die Daten. Man kann aber auch die
Rückgabe von ezPlot() nutzen, um daraus mit
ggplot schönere Abbildungen zu generieren. Wir gehen hier
nicht weiter auf diese Möglichkeiten ein.
2.4 Literatur
Baguley, T. (2012). Serious stats. A guide to advanced statistics for the behavioral sciences. Palgrave Macmillan.
Eid, M., Gollwitzer, M. & Schmitt, M. (2010). Statistik und Forschungsmethoden. Beltz.
Keppel, G. & Wickens, T.D. (2004). Design and analysis. A researcher’s handbook. Pearson.
Lawrence MA (2016). ez: Easy analysis and visualization of dactorial experiments. R package version 4.4-0, https://CRAN.R-project.org/package=ez.
Singmann H, Bolker B, Westfall J, Aust F, & Ben-Shachar M (2023). afex: Analysis of factorial experiments. R package version 1.3-0, https://CRAN.R-project.org/package=afex.
3 Einfaktorielle Varianzanalyse II
In Teil 2 hatten wir die Grundlagen der einfaktoriellen Varianzanalyse erarbeitet, indem wir die Quadratsummenzerlegung, Mittlere Quadratsummen sowie den \(F\)-Bruch und die \(F\)-Verteilung kennengelernt haben. In Teil 3 bauen wir darauf auf und betrachten weitergehende Konzepte und wichtige Punkte im Kontext der Varianzanalyse. Den Abschluss bildet dann eine Betrachtung der sog. Modellgleichung, auf die wir dann in Teil 4 wieder zurückkommen werden.
3.1 \(t\)-Test vs. Varianzanalyse
3.1.1 gerichtete vs. ungerichtete Tests
Beim \(t\)-Test hatten wir zwischen ungerichteten und gerichteten Tests bzw. Alternativhypothesen unterschieden: \[H_1:\mu_A \neq \mu_B\hspace{1cm}\text{vs.}\hspace{1cm}H_1:\mu_A > \mu_B\] Ab \(J> 2\) Gruppen sind gerichteten Vorhersagen im Rahmen einer Varianzanalyse nicht mehr möglich, da eine eindeutige Richtung i.d.R. nicht bestimmbar ist (vgl. aber später: Kontraste). Das heißt, dass eine Varianzanalyse eigentlich immer ungerichtet testet und keine Vorhersagen über Richtungen von Unterschieden berücksichtigt.
3.1.2 Zusammenhang für \(J=2\)
Der \(t\)-Test für zwei unabhängige Stichproben eignet sich für zwei Stichproben. Die Varianzanalyse kann mehr als zwei Stichproben berücksichtigen, sie stellt also gewissermaßen eine Verallgemeinerung des \(t\)-Tests dar. Tatsächlich kann man den \(t\)-Test als “Spezialfall” der Varianzanalyse mit \(J=2\) auffassen. In diesem Fall gilt die einfache Beziehung: \[ F=t^2 \] Man kann also einen ungerichteten \(t\)-Test mit \(m\) Freiheitsgraden auch als Varianzanalyse mit \(J=2\) Gruppen auffassen. Der \(F\)-Test hat dann 1 Zählerfreiheitsgrad und \(m\) Nennerfreiheitsgrade und der \(t\)-Wert muss quadriert werden, um den \(F\)-Wert zu erhalten. Neben dem formalen Beweis, können wir dies auch anhand der Daten aus dem Eingangsbeispiel aus Teil 2 klarmachen:
daten_training <- read.table("./Daten/daten_training.DAT",
header = TRUE)
daten_training.2 <- subset(daten_training,
Trainingstage != 2) # nur 1 vs. 3 Tage Training
daten_training.2$Trainingstage <- as.factor(daten_training.2$Trainingstage)
str(daten_training.2$Trainingstage) # Faktor mit nur noch J = 2 Stufen## Factor w/ 2 levels "1","3": 1 1 1 1 1 1 1 1 1 1 ...# erst als t-Test
t.result <- t.test(Sprachtest ~ Trainingstage,
data = daten_training.2,
var.equal = TRUE)
t.result$statistic^2 # t-Wert quadriert = 18.7813## t
## 18.7813# dann als Varianzanalyse
aov.result <- aov(Sprachtest ~ Trainingstage,
data = daten_training.2)
summary(aov.result) # F = 18.78## Df Sum Sq Mean Sq F value Pr(>F)
## Trainingstage 1 125.0 125.00 18.78 4e-04 ***
## Residuals 18 119.8 6.66
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.2 Effektstärken, Power und Konfidenzintervalle
3.2.1 Effektstärken und ihre Schätzung
In Teil 14 von Statistik I haben wir \[\delta = \frac{\mu_A-\mu_B}{\sigma}\] als Effektstärke für Mittelwertunterschiede eingeführt und geschätzt durch \[d=\frac{M_A-M_B}{\hat{\sigma}}\] Bei mehr als zwei Gruppen, also \(J> 2\), gibt es aber keine einzelne Mittelwertdifferenz mehr, die alle Gruppen gleichermaßen erfasst. Mit \(\delta\) kommen wir also hier nicht weiter und wir müssen stattdessen einen anderen Ansatz einer Effekstärke wählen.
Wir beginnen dazu mit folgender Überlegung: Als Gesamterwartungswert \(\mu\) bezeichnen wir den (Populations-)Mittelwert aller \(J\) Populationen gemeinsam. Der Erwartungswert einer Zufallsvariablen \(\mathbf{M}\), die jeder Kombination von Stichproben den Gesamtmittelwert \(M\) zuweist, ist \[E(\mathbf{M})=\mu\] Da \(M\) auch als gewichtetes Mittel der Gruppenmittelwerte geschrieben werden kann (vgl. Teil 2), kann auch \(\mu\) entsprechend geschrieben werden: \[E(\mathbf{M})=\mu=\sum_{j=1}^J\frac{n_j}{N}\mu_j\]
Nun weicht jeder Populationsmittelwert \(\mu_j\) vom Gesamterwartungswert \(\mu\) ab. Diese Abweichung wird Effekt der Stufe \(j\) des Faktors genannt und mit \(\alpha_j\) bezeichnet (dieses \(\alpha\) hat nichts mit dem Signifikanzniveau \(\alpha\) zu tun): \[\boxed{\alpha_j=\mu_j-\mu\hspace{1cm}\forall j\in\{1,\ldots,J\}}\] Es gibt also genau so viele Effekte \(\alpha_j\) wie es Stufen des Faktors gibt. Die Effekte einfach addieren ist nicht zielführend, da ihre Summe Null ergibt: \[ \sum_{j=1}^J\alpha_j = 0 \] Daher müssen wir noch einen Schritt weitergehen und wir berechnen eine Effektvarianz: \[ \sigma^2_{\alpha}=\sum_{j=1}^J\frac{n_j}{N}(\underbrace{\mu_j-\mu}_{\alpha_j})^2 \] Für die Effektvarianz gilt: je näher \(\sigma_\alpha^2\) an Null liegt, desto weniger unterscheiden sich die Populationsmittelwerte \(\mu_j\) voneinander. Gibt es gar keine Effektvarianz, dann sind alle \(\mu_j\) und und auch \(\mu\) gleich. Dies ermöglicht auch eine alternative Formulierung der Hypothesen: \[H_0:\sigma_\alpha^2=0\hspace{1cm}\text{und}\hspace{1cm}H_1:\sigma_\alpha^2>0\] Die Effektstärke \(\delta\) hatten wir an der Standardabweichung relativiert. Ganz ähnlich gehen wir auch hier vor und die Effektvarianz wird noch an der Fehlervarianz relativiert. Dieses Maß wird Effektstärke \(\mathbf{f^2}\) genannt: \[f^2=\frac{\sigma_\alpha^2}{\sigma^2}\] Als weiteres Maß wird die Effektstärke \(\mathbf{\eta^2}\) benutzt: \[\eta^2=\frac{\sigma_\alpha^2}{\sigma_\alpha^2+\sigma^2}=\frac{\sigma_\alpha^2}{\sigma^2_{\text{tot}}}\] \(\eta^2\) ist also der Anteil der Gesamtvarianz der Daten, der durch die Effektvarianz erklärt wird. Da es sich um einen Anteil handelt, kann \(\eta^2\) nur Werte zwischen 0 und 1 annehmen.
Beide Maße können leicht ineinander umgerechnet werden: \[\eta^2=\frac{f^2}{1+f^2}\hspace{1cm}\text{bzw.}\hspace{1cm}f^2=\frac{\eta^2}{1-\eta^2}\] Schließlich haben wir für \(\delta\) (und auch \(r\)) bereits die Konventionen von Cohen (1988) kennengelernt, was als kleiner, mittlerer und großer Effekt gilt. Derartige Konventionen gibt es auch für die gerade eingeführten Effektstärken, wobei zu beachten ist, dass die Werte für \(f\), also der Wurzel aus \(f^2\), angegeben sind.
Genau wie bei \(\delta\) handelt es sich auch bei \(f^2\) und bei \(\eta^2\) um Populationsparameter. Für die Schätzung von \(\eta^2\) aus Daten gibt es zwei Varianten:
- \[\hat{\eta}^2=\frac{SS_A}{SS_{\text{tot}}}\]
Für das Beispiel aus Teil
2 ergibt sich damit: \(\hat{\eta}^2=\frac{125.07}{264.97}=.47\).
Dies ist der Wert, der von
ezANOVA()und auch vonanova_out()ausgegeben wird. (Im einfaktoriellen Fall unterscheiden sich \(\eta^2\), partial \(\eta^2\) (\(\eta^2_P\)) und generalized \(\eta^2\) nicht, daher entsprechen sich die Werte.) Allerdings überschätzt diese Variante den Populationseffekt systematisch: Auch wenn \(\eta^2=0\) ist, findet sich \(\hat{\eta}^2>0\). - \[\hat{\eta}^2=\frac{SS_A-(J-1)MS_w}{SS_{\text{tot}}+MS_w}\] Für das Beispiel aus Teil 2 ergibt sich damit: \(\hat{\eta}^2=\frac{125.07-2\cdot 5.18}{264.97+5.18}=.42\). Dieser Schätzer wird in der Literatur auch als \(\omega^2\) bezeichnet; er ist allerdings auch nicht erwartungstreu.
3.2.2 Power in der Varianzanalyse
3.2.2.1 Wiederholung: \(t\)-Test und Power
Weiter oben haben wir uns damit beschäftigt, wie der \(F\)-Bruch verteilt ist, wenn die Nullhypothese zutrifft, d.h. wenn alle \(\mu_j\) gleich sind. In Teil 14 von Statistik I haben wir uns aber auch damit beschäftigt, wie der \(t\)-Bruch verteilt ist, wenn eine bestimmte Alternativhypothese zutrifft, also ein bestimmter Effekt \(\delta \neq 0\) vorhanden ist. Angenommen wir hätten zwei normalverteilte Populationen \(A\) und \(B\), die sich um einen Effekt \(\delta=3.0\) unterscheiden:
Wenn dann zwei unabhängige Stichproben aus diesen beiden Populationen
mit verschiedenem Erwartungswert \(\mu\) stammen, also die Alternativhypothese
\(H_1: \mu_B=\mu_A+3.0\) gilt, dann ist
die Zufallsvariable \(\mathbf{t}\), die
jeder Kombination zweier Stichproben den \(t\)-Bruch zuordnet, nicht mehr zentral
\(t\)-verteilt, sondern
non-zentral \(t\)-verteilt. Diese non-zentrale \(t\)-Verteilung hat (auch) \(m\) Freiheitsgrade, aber zusätzlich einen
Nonzentralitätsparameter \(\Delta\),
der proportional zum Effekt \(\delta\)
in der Population ist. Die folgende Abbildung visualisiert noch einmal
verschiedene non-zentrale \(t\)-Verteilungen:
Die folgende Abbildung illustriert die Situation der Verteilung des \(t\)-Bruchs bei Gültigkeit der Nullhypothese (links) und bei Gültigkeit einer bestimmten Alternativhypothese (rechts) und fasst die Überlegungen zur Bestimmung im Fall des \(t\)-Tests nocheinmal zusammen:
- Zunächst wird der kritische \(t\)-Wert immer aufgrund des gewählten \(\alpha\) unter Annahme der Nullhypothese berechnet, d.h. mit Hilfe der zentralen \(t\)-Verteilung. Bei \(\alpha=0.05\) ist also die rote Fläche genau 5% der \(t\)-Verteilung.
- Die blaue Fläche unter der non-zentralen \(t\)-Verteilung bei einer bestimmten Alternativhypothese (die also einen bestimmten Effekt annimmt) ist die Wahrscheinlichkeit, bei Gültigkeit dieser Alternativhypothese ein signifikantes Ergebnis zu bekommen (weil \(t\geq t_{\text{krit}}\)). Diese Wahrscheinlichkeit wird Power bzw. Teststärke genannt und man schreibt für die Power auch \(1-\beta\).
- Die maximale Wahrscheinlichkeit eines Fehlers 2. Art ist dann \(\beta\), also die orangene Fläche unter der non-zentralen \(t\)-Verteilung bis zu \(t_{\text{krit}}\).
3.2.2.2 Power und die Varianzanalyse
Ganz ähnliche Überlegungen gibt es für die Varianzanalyse auch:
- Gilt die Nullhypothese, dann ist \(\mathbf{F}\) zentral \(F\)-verteilt. Diesen Fall haben wir bisher betrachtet.
- Gilt jedoch eine bestimmte Alternativhypothese – also eine, die einen bestimmten Effekt annimmt – dann ist \(\mathbf{F}\) non-zentral \(F\)-verteilt
Gena wie die non-zentrale \(t\)-Verteilung hat auch die non-zentrale \(F\)-Verteilung neben den (beiden) Freiheitsgraden einen zusätzlichen Nonzentralitätsparameter \(\Delta\), der proportional zur vorliegenden Effektstärke ist. Je größer dieser Nonzentralitätsparameter \(\Delta\) ist (bzw. die Effektstärke), desto weiter “nach rechts” wird die resultierende \(F\)-Verteilung verschoben:
Ganz analog wie oben für die \(t\)-Verteilung, illustriert die folgende Abbildung die Situation der Verteilung des \(F\)-Bruchs bei Gültigkeit der Nullhypothese (links) und bei Gültigkeit einer bestimmten Alternativhypothese (rechts) und fasst die Überlegungen zur Bestimmung der Power zusammen:
- Zunächst wird der kritische \(F\)-Wert immer aufgrund des gewählten \(\alpha\) unter Annahme der Nullhypothese berechnet, d.h. mit Hilfe der zentralen \(F\)-Verteilung. Bei \(\alpha=0.05\) ist also die rote Fläche genau 5% der \(F\)-Verteilung.
- Die blaue Fläche unter der non-zentralen \(F\)-Verteilung bei einer bestimmten Alternativhypothese (die also einen bestimmten Effekt annimmt) ist die Wahrscheinlichkeit, bei Gültigkeit dieser Alternativhypothese ein signifikantes Ergebnis zu bekommen (weil \(F\geq F_{\text{krit}}\)). Diese Wahrscheinlichkeit ist wiederum die Power bzw. Teststärke \(1-\beta\).
- Die maximale Wahrscheinlichkeit eines Fehlers 2. Art ist dann \(\beta\), also die orangene Fläche unter der non-zentralen \(F\)-Verteilung bis zu \(F_{\text{krit}}\).
Zusammengefasst spiegeln die vier Farben in der vorangegangenen Abbildung folgende Fälle wider:
- rot: \(\alpha\cdot 100\)% der Fläche ab kritischem \(F\)-Wert bei Gültigkeit der Nullhypothese \(\Rightarrow\) \(F\geq F_{\text{krit}}\) ist dann ein Fehler 1. Art
- grün: \((1-\alpha)\cdot 100\)% der Fläche bis zu kritischem \(F\)-Wert bei Gültigkeit der Nullhypothese \(\Rightarrow\) \(F< F_{\text{krit}}\) ist dann eine korrekte Entscheidung
- orange: Fläche \(\beta\) bis zu kritischem \(F\)-Wert bei Gültigkeit der Alternativhypothese \(\Rightarrow\) \(F< F_{\text{krit}}\) ist dann ein Fehler 2. Art
- blau: Fläche \(1-\beta\) ab kritischem \(F\)-Wert bei Gültigkeit der Alternativhypothese \(\Rightarrow\) \(F\geq F_{\text{krit}}\) ist dann eine korrekte Entscheidung
Eine ShinyApp zur interaktiven Visualisierung der Power-Berechnung finden Sie hier
3.2.2.3 Bestimmung der Power
Die Bestimmung der Power bei einer einfaktoriellen Varianzanalyse wie bisher betrachtet ist noch relativ einfach. Die “traditionelle Möglichkeit” hierzu bieten wieder die Tabellen im Buch von Cohen (1988). Moderner und einfacher ist es allerdings, die Power mit dem Computer zu bestimmen. Zum einen kann hierfür wieder das Programm G*Power benutzt werden.
Darüber hinaus gibt es zur Bestimmung der Power mit R z.B. die
Funktion pwr.anova.test() aus dem Paket pwr.
Hier werden ähnlich wie bei der Funktion power.t.test()
manche der Parameter übergeben, während eine der Größen als
NULL bezeichnet wird. Dies ist die Größe, die berechnet
werden soll; im folgenden Beispiel ist es die Anzahl der
Versuchspersonen n:
library(pwr)
pwr.anova.test(k = 3, # Anzahl Gruppen
n = NULL, # Anzahl VPen - wird bestimmt
f = 0.25, # (angenommener) Effekt (Achtung: nicht f^2!)
sig.level = 0.05, # Alpha
power = 0.80) # 1-Beta = Power##
## Balanced one-way analysis of variance power calculation
##
## k = 3
## n = 52.3966
## f = 0.25
## sig.level = 0.05
## power = 0.8
##
## NOTE: n is number in each groupIm Hauptskript zur Veranstaltung “Statistik 2” behandeln wir
lediglich die Power-Berechnung im einfaktoriellen Design. Die
Power-Berechnung in komplexeren Designs behandeln wir nicht explizit, da
sie deutlich komplizierter ausfällt und vom genauen Design abhängt. Ein
sehr mächtiges Paket ist in diesem Zusammenhang allerdings
Superpower, welches die Power-Berechnung für beliebige
ANOVA-Designs mit Hilfe von simulativen Methoden ermöglicht. Eine kurze
Einführung und weitere Informationen zum Paket finden Sie hier.
3.2.3 Konfidenzintervalle
In Teil 13 von Statistik I hatten wir Konfidenzintervalle betrachtet. Für den Fall zweier unabhängiger Stichproben hatten wir zwei Möglichkeiten für Konfidenzintervalle unterschieden:
Konfidenzintervalle für jeden einzelnen Mittelwert (äquivalent zum \(t\)-Test für eine Stichprobe auf einen bestimmten Populationswert \(\mu_0\)). Diese Möglichkeit besteht selbstverständlich auch bei \(J> 2\) und alle Berechnungen sind wie beschrieben in Teil 13 von Statistik I. Sofern wir allerdings von Varianzhomogenität ausgehen, ist \(MS_w\) ein Schätzer der Populationsvarianz, der auf den Daten aller Gruppen basiert und daher präziser ist. Wenn wir zusätzlich beachten, dass die Freiheitsgrade des Sicherheitsparameters dann auf \(MS_w\) beruhen und i.d.R. größer sind, können diese Konfidenzintervalle präziser wie folgt berechnet werden: \[\left[ M_j\pm t_{N-J;1-\frac{\alpha}{2}}\cdot \sqrt{\frac{MS_w}{n_j}}\right]\]
Das Konfidenzintervall für den Mittelwertunterschied auf den Wert 0 hatten wir als äquivalent zum \(t\)-Test für zwei unabhängige Stichproben identifiziert. Bei \(J> 2\) gibt es allerdings keinen einzelnen Mittelwertunterschied mehr, aber es gibt eine Korrespondenz zwischen Varianzanalyse und Konfidenzintervall für den Fall, dass die Stichprobengrößen in allen Gruppen gleich sind (\(n_j = n, \quad\forall j \in \{1,\ldots,J\}\)) (siehe dazu auch Loftus & Masson, 1994): zwei Mittelwerte unterscheiden sich signifikant, wenn gilt \[|M_k-M_m|>\sqrt{2}\cdot t_{N-J;1-\frac{\alpha}{2}}\cdot\sqrt{\frac{MS_w}{n}}\qquad k,m\in\{1,\ldots,J\}\]
3.3 Kontraste
Bisher haben wir im Rahmen der Varianzanalyse globale Nullhypothesen getestet; diese Form der Varianzanalyse wird auch Omnibus-Varianzanalyse genannt. Im Fall eines signifikanten Ergebnisses wissen wir allerdings lediglich, dass sich mindestens zwei der Erwartungswerte vermutlich unterscheiden. Welche dies sind, wissen wir nicht. Erinnern wir uns an das Ergebnis der Varianzanalyse aus Teil 2, wissen wir nur, dass mindestens zwei der Erwartungswerte der drei Populationen 1, 2 bzw. 3 Trainingstage sich vermutlich unterscheiden. Vielleicht unterscheiden sich aber auch nur die Leistungen im Sprachtest zwischen 1 und 3 Trainingstagen, vielleicht sind aber auch alle drei Stufen unterschiedlich. Wir wissen es bei einer signifikanten Varianzanalyse schlichtweg nicht. Hier kommen sog. Kontraste ins Spiel, mit denen präzisere Hypothesen getestet werden können. Kontraste sind dabei ein mächtiges Instrument, wir gehen hier aber nur auf die einfachsten Fälle ein. Mehr Informationen finden sich z.B. in Abschnitt 8.4 von Janczyk und Pfister (2020) oder bei Schad et al. (2020).
3.3.1 A priori vs. a posteriori Kontraste
Gezielte Unterschiede und noch weitergehende Fragen können also mit Kontrasten beantwortet werden. Konzeptuell werden zwei Arten von Kontrasten unterschieden, sog. a priori und a posteriori Kontraste.
- A posteriori Kontraste: Nach einer signifikanten Varianzanalyse (wie im Beispiel aus Teil 2) werden oft Kontraste berechnet, um die signifikanten Unterschiede zu lokalisieren. Dies sind dann sog. a posteriori (auch post-hoc genannte) Kontraste, die also im Nachgang einer signifikanten Varianzanalyse berechnet werden, ohne gezielte Hypothese gehabt zu haben. Sie können daher nur ungerichtet getestet werden und sollten eigentlich gar nicht als Hypothesentests aufgefasst werden, da keine Hypothesen über einzelne Erwartungswerte vorab formuliert wurden. Für solche Zwecke gibt es verschiedene Verfahren, die alle automatisch eine \(\alpha\)-Adjustierung beinhalten, sich aber darin unterscheiden, wie sie dies tun. Eine weitverbreitete Variante ist der sog. Scheffé-Test, der die Prüfgröße entsprechend so korrigiert, dass beliebig viele Vergleiche durchgeführt werden, ohne dass das nominelle \(\alpha\)-Niveau überschritten wird. Allerdings ist der Scheffé-Test sehr konservativ; eine weniger konservative Alternative ist bspw. Tukey’s “honestly significant difference (HSD)”-Methode. Beide Varianten werden weiter unten für die praktische Berechnung mit R demonstriert.
- A priori Kontraste: Wir können auch von vornherein (nur) theoretisch interessante Kontraste testen, deren Hypothesen auf theoretisch begründeten Vorhersagen beruhen. Dies sind sog. a priori Kontraste oder “geplante Vergleiche”. Sie beruhen zwar letztlich auf einem varianzanalytischem Design und nutzen dessen Vorteile wie z.B. die einer Varianzschätzung mit allen Daten durch \(MS_w\), testen aber gezielte Hypothesen über Unterschiede und werden im Idealfall sogar anstelle einer Omnibus-Varianzanalyse durchgeführt. Generell sind sie a posteriori Kontrasten vorzuziehen, da eben ihre Hypothesen auf inhaltlichen Überlegungen basieren.
3.3.2 Zwei einfache Beispiele für Kontraste
Um einen Überblick über verschiedene Möglichkeiten zur Formulierung von Kontrasten zu erlangen, betrachten wir wieder das Eingangsbeispiel aus Teil 2. Die Mittelwerte der drei Gruppen waren 13.1, 15.7 und 18.1 Punkte im Sprachtest, die Leistung in diesem Test nimmt also mit einer zunehmendem Anzahl an Trainingstagen zu und und die Varianzanalyse war signifikant. Wir beschreiben hier zwei einfache Szenarien, die theoretisch interessant sein könnten und mit Kontrasten erfassbar sind. Mit \(\mu_1\), \(\mu_2\) und \(\mu_3\) bezeichnen wir hierbei die Erwartungswerte der Populationen, aus denen die Stichproben der drei Faktorstufen 1, 2 bzw. 3 Trainingstage stammen:
- Beispiel 1: Ein einfacher Fall wäre die Frage, ob die Leistung bereits nach einer Trainingswiederholung (d.h., mit zwei statt nur einem Trainingstag) besser wird, also ob \(\mu_1<\mu_2\) ist.
- Beispiel 2: Vielleicht gäbe es aber auch Grund zur Annahme, die Leistung würde erst nach zweimaliger Wiederholung des Trainings steigen (d.h., am dritten Trainingstag) während nur einer oder zwei Trainingstage noch keinen Unterschied in der Leistung ausmachen. Insgesamt umfasst diese Vorhersage zwei Kontraste: Der zweite Teil besagt, dass \(\mu_1=\mu_2\) sein soll und der erste Teil besagt, dass \(\frac{\mu_1+\mu_2}{2}<\mu_3\) sein soll.
Was ist nun ein Kontrast? Ein Kontrast ist i.W. nichts anderes als eine bestimmte Kombination von Populationsparametern. Diese (Linear-)Kombination nennen wir \(\psi\) (ein kleines Psi), und sie ergibt sich daraus, dass alle \(\mu_j\) jeweils mit einem Koeffizienten \(c_j\) multipliziert und dann aufsummiert werden: \[ \psi=c_1\mu_1+\ldots+c_J\mu_J=\sum_{j=1}^Jc_j\mu_j\;. \] Wichtig ist allerdings, dass wir nur von einem Kontrast sprechen, wenn sich die \(c_j\) insgesamt zu 0 aufaddieren, also wenn: \[\sum_{j=1}^Jc_j=0\]
3.3.3 Hypothesen, Schätzung und Testen von Kontrasten
Wir betrachten zunächst das Hypothesenpaar einer Kontrastanalyse und kommen dann zu deren rechnerischer Durchführung. Bei allen folgenden Schritten gelten die gleichen Annahmen wie bei der Varianzanalyse.
Das zu testende Hypothesenpaar ist üblicherweise für alle Kontraste identisch: \(\psi\) wird auf den Wert 0 getestet. Im Fall einer ungerichteten Hypothese schreiben wir daher: \[ H_0:\psi=0\qquad\text{und}\qquad H_1:\psi\neq 0 \]
Im Falle einer gerichteten Hypothese schreiben wir hingegen \[ H_0:\psi\leq0\qquad\text{und}\qquad H_1:\psi> 0 \] bzw. \[ H_0:\psi\geq0\qquad\text{und}\qquad H_1:\psi< 0 \]
Die eigentliche Schwierigkeit bei Kontrasten besteht darin, die herkömmlichen Hypothesen in Kontrasthypothesen zu überführen. Dies geschieht über eine Spezifizierung der Koeffizienten \(c_j\). Wir betrachten dies an den beiden Beispielen:
Beispiel 1: Bei diesem Beispiel fragen wir uns, ob die Leistung bereits nach einmaliger Wiederholung besser wird (\(\mu_1<\mu_2\)). Das entsprechende Hypothesenpaar lautet hier also: \[ H_0:\mu_1\geq\mu_2\qquad\text{und}\qquad H_1:\mu_1<\mu_2 \] Wir können nun diese Hypothese umformulieren als Kontrasthypothese im Zuge der Varianzanalyse. Als ersten Schritt formulieren wir die Hypothese so um, dass die Parameter als Summen/Differenzen dargestellt werden: \[ H_0:\mu_1 - \mu_2\geq0 \qquad\text{und}\qquad H_1:\mu_1-\mu_2<0 \] Schließlich schreiben wir diese Hypothese so um, dass wir die Koeffizienten \(c_j\) explizit darstehen haben: \[ H_0:\underset{=\psi}{\underbrace{1\cdot\mu_1+(-1)\cdot\mu_2+0\cdot\mu_3}}\geq0 \qquad\text{und}\qquad H_1:\underset{=\psi}{\underbrace{1\cdot\mu_1+(-1)\cdot\mu_2+0\cdot\mu_3}}<0 \] Hieraus ergibt sich die Formel für den Kontrast \(\psi\), und die Koeffizienten sind \(c_1=1\), \(c_2=-1\) und \(c_3=0\).
Beispiel 2: Bei diesem Beispiel fragen wir uns, ob die Leistung nach zweimaliger Wiederholung des Trainings steigt (d.h., ab dem dritten Trainingstag). Damit diese Vermutung gilt, müssen zwei Aspekte gelten: \(\frac{\mu_1+\mu_2}{2}<\mu_3\) und \(\mu_1=\mu_2\). Die Kontrastkoeffizienten für den zweiten Teil sind analog zum Beispiel 1, nur, dass wir hier ungerichtetet testen. Für den interessanteren (ersten) Teil schreiben wir das folgende Hypothesenpaar: \[ H_0:\frac{\mu_1+\mu_2}{2}\geq\mu_3 \qquad\text{und}\qquad H_1:\frac{\mu_1+\mu_2}{2}<\mu_3 \] Nun formulieren wir die Hypothese so um, dass die Parameter als einzelne Summen/Differenzen dargestellt werden können: \[ H_0:\underset{=\psi}{\underbrace{0.5\cdot\mu_1+ 0.5\cdot\mu_2 - \mu_3}} \geq 0 \qquad\text{und}\qquad H_1:\underset{=\psi}{\underbrace{0.5\cdot\mu_1+ 0.5\cdot\mu_2 - \mu_3}} < 0 \] Als Koeffizienten ergeben sich dann also \(c_1=0.5\), \(c_2=0.5\) und \(c_3=-1\).
Die Kombination der Populationserwartungswerte \(\mu_j\) bzw. die Formulierung über \(\Psi\) gelten auf Populationsebene. Auf Stichprobenebene benötigen wir einen Schätzer für \(\psi\) wobei es direkt naheliegt, die \(\mu_j\) jeweils durch die korrespondierenden Mittelwerte \(M_j\) zu schätzen. Demnach wäre ein geeigneter Schätzer für \(\psi\): \[ \hat{\psi}=\sum_{j=1}^Jc_jM_j \] Betrachten wir nun eine Zufallsvariable \(\mathbf{\hat{\psi}}\), die jeder möglichen Stichprobenziehung aus der oder den betrachteten Population(en) \(\hat{\psi}\) zuordnet, dann ergeben sich als deren Erwartungswert und Varianz: \[ E(\mathbf{\hat{\psi}})=\psi\quad\text{und}\quad V(\mathbf{\hat{\psi}})=\sigma^2\sum_{j=1}^J\frac{c_j^2}{n_j} \] \(\hat{\psi}\) ist also ein erwartungstreuer Schätzer für \(\psi\). Da sich \(\mathbf{\hat{\psi}}\) zudem aus den normalverteilten Zufallsvariablen \(\mathbf{M_j}\) zusammensetzt, ist sie selbst auch wie folgt normalverteilt: \[ \mathbf{\hat{\psi}}\sim N\left(\psi,\sigma^2\sum_{j=1}^J\frac{c_j^2}{n_j}\right) \] Als erwartungstreue Schätzung für \(V(\mathbf{\hat{\psi}})\) benutzen wir schließlich \[ S_{\mathbf{\hat{\psi}}}^2=MS_w\sum_{j=1}^J\frac{c_j^2}{n_j} \] wobei \(MS_w\) die bekannte Mittlere Quadratsumme aus der Varianzanalyse ist.
Wie können wir nun mit diesen Informationen die Hypothesenpaare für Kontraste testen? Das Vorgehen ist wieder identisch mit dem Vorgehen bei \(t\)-Tests und der oben dargestellten Varianzanalyse: Wir suchen eine (Prüf-)Größe, in die \(\hat{\psi}\) eingeht und die zwei Bedingungen erfüllt:
- Ihr Wert wird besonders extrem, wenn die Daten gegen die Gültigkeit der \(H_0\) sprechen.
- Wir kennen die Verteilung einer entsprechenden Zufallsvariablen unter Annahme der Gültigkeit der \(H_0\).
In Teil 12 von Statistik I hatten wir bereits die allgemeine Form des \(t\)-Bruchs betrachtet. Angewendet auf Kontraste kann diese Form spezifiziert werden als: \[ t=\frac{T-\tau_0}{SE_T} =\frac{\hat{\psi}-0}{\sqrt{S^2_{\mathbf{\hat{\psi}}}}}= \frac{\hat{\psi}}{\sqrt{MS_w\sum_{j=1}^J\frac{c_j^2}{n_j}}} \]
Eine Zufallsvariable \(\mathbf{t}\), die jeder möglichen Kombination von Stichproben diesen \(t\)-Wert zuordnet, ist daher \(t\)-verteilt mit \(N-J\) Freiheitsgraden: \[t\overset{H_0}{\sim} t_{N-J},\] wobei \(N\) die Gesamtzahl der Versuchspersonen und \(J\) die Anzahl der Faktorstufen bezeichnet.
Die Entscheidungsregel im zweiseitigen Fall lautet: “Verwirf die \(H_0\), falls gilt \(|t|\geq t_{N-J;1-\frac{\alpha}{2}}\) bzw. falls \(p\leq\alpha\) ist.
Die Entscheidungsregel im einseitigen Fall für \(H_1: \Psi > 0\) lautet: “Verwirf die \(H_0\), falls gilt \(t\geq t_{N-J;1-\alpha}\) bzw. falls \(p\leq\alpha\) ist. Für den Fall \(H_1: \Psi < 0\) gilt analog:”Verwirf die \(H_0\), falls gilt \(t\leq t_{N-J;\alpha}\) bzw. falls \(p\leq\alpha\) ist.
Anmerkung: Insbesondere für das erste Beispiel könnte man i.Ü leicht einwenden, dass doch ein einfacher \(t\)-Test ausreichen würde um die Hypothesen zu testen. Tatsächlich wäre der Zähler in beiden Fällen identisch. Da im Nenner des \(t\)-Bruchs für Kontraste allerdings die gemeinsame Varianzschätzung mittels \(MS_w\) steht (also auf Basis aller Versuchspersonen), wird i.d.R. der Nenner kleiner und der \(t\)-Bruch entsprechend größer als im Fall eines einfachen \(t\)-Tests. Zudem wird eine \(t\)-Verteilung mit mehr Freiheitsgraden herangezogen, wodurch der kritische \(t\)-Wert kleiner wird. Zusammengefasst ist die Power von Kontrasten i.d.R. (aber nicht immer!) höher als die der entsprechenden \(t\)-Tests.
3.3.4 Praktische Berechnung von Hand und mit R
Wir betrachten nun die Umsetzung des Beispiels 1 von oben. Das Hypothesenpaar über die Kontrasthypothese lautet: \[ H_0:\underset{=\psi}{\underbrace{1\cdot\mu_1+(-1)\cdot\mu_2+0\cdot\mu_3}}\geq0 \qquad\text{und}\qquad H_1:\underset{=\psi}{\underbrace{1\cdot\mu_1+(-1)\cdot\mu_2+0\cdot\mu_3}}<0 \] Die Koeffizienten für den Kontrast lauten \(c_1=1\), \(c_2=-1\) und \(c_3=0\). Die Formel für den empirischen \(t\)-Bruch lautet, \[ t=\frac{T-\tau_0}{SE_T} =\frac{\hat{\psi}-0}{\sqrt{S^2_{\mathbf{\hat{\psi}}}}}= \frac{\hat{\psi}}{\sqrt{MS_w\sum_{j=1}^J\frac{c_j^2}{n_j}}} \quad , \] mit \[ \hat{\psi}=\sum_{j=1}^Jc_jM_j \quad. \] Die Koeffizienten haben wir bereits bestimmt. Was noch fehlt sind die Mittelwerte der Bedingungen, \(M_j\) sowie die mittlere Quadratsumme innerhalb der Gruppen, \(MS_w\). Diese wurden bereits in Teil 2 berechnet, aber wir berechnen diese hier nochmal erneut, um uns etwas Suchen zu sparen:
daten_training <- read.table("./Daten/daten_training.DAT", header = TRUE)
daten_training$Trainingstage <- as.factor(daten_training$Trainingstage) # faktorisieren
aov_ergebnis <- aov(Sprachtest ~ Trainingstage, # AV ~ Faktor
data = daten_training) # welche Daten?
summary(aov_ergebnis) # Ergebnis ausgeben## Df Sum Sq Mean Sq F value Pr(>F)
## Trainingstage 2 125.1 62.53 12.07 0.00018 ***
## Residuals 27 139.9 5.18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tapply(daten_training$Sprachtest,
daten_training$Trainingstage,
FUN = mean)## 1 2 3
## 13.1 15.7 18.1Die Mittelwerte lauten also \(13.1\), \(15.7\) und \(18.1\). Die gesuchte mittlere Quadratsumme \(MS_w = 5.18\). Damit haben wir alle Informationen, um \(t\) zu berechnen: \[ t=\frac{\hat{\psi}}{\sqrt{MS_w\sum_{j=1}^J\frac{c_j^2}{n_j}}} =\frac{1\cdot 13.1+(-1)\cdot 15.7}{\sqrt{5.18\cdot(\frac{1^2}{10}+\frac{-1^2}{10}+\frac{0^2}{10})}} =\frac{-2.6}{\sqrt{5.18\cdot 0.2}} =\frac{-2.6}{\sqrt{1.036}} =-2.554427 \] Analog mit R:
t.zaehler <- 1*13.1 + (-1)*15.7 + 0*18.1 # Zähler
t.nenner <- sqrt(5.18 * (1/10 + 1/10 + 0/10)) # Nenner
t <- t.zaehler / t.nenner # t-Wert berechnen
t # t-Wert ausgeben## [1] -2.554427Die Frage ist nun, ob dieser \(t\)-Wert auch statistisch signifikant wird. Hierfür berechnen wir entweder einen kritischen \(t\)-Wert oder einen \(p\)-Wert anhand der zentralen \(t\)-Verteilung mit \(N-J\) Freiheitsgraden (\(N\) ist hier 30 und \(J\) ist 3). Beachten müssen wir allerdings noch ob wir einseitig bzw. zweiseitig testen. Im Beispiel haben wir eine einseitige Hypothese formuliert mit \(H_1: \psi <0\). Kleine Werte für den Kontrast sprechen also für unsere \(H_1\) bzw. gegen die \(H_0\). Somit ergibt sich als kritischer \(t\)-Wert:
qt(0.05, df = 27)## [1] -1.703288bzw. für den \(p\)-Wert:
p <- pt(t, df = 27) # p-Wert berechnen (einseitig)
p # p-Wert ausgeben## [1] 0.008295685Der Test der Kontrasthypothese wäre also statistisch signifikant und wir würden davon ausgehen, dass das Training tatsächlich bereits nach der ersten Wiederholung zu einer besseren Sprachleistung führt.
Eine direkte Variante mit R bietet die Funktion glht()
aus dem Paket multcomp, welche direkt auf dem
aov()-Objekt aufbaut. Dazu berechnen wir die Varianzanalyse
zunächst mit aov() und formulieren dann den Kontrast für
das Beispiel 1 von oben:
suppressPackageStartupMessages(library(multcomp)) # für glht()
contrast <- glht(aov_ergebnis, # im R Chunk oben berechnet
linfct = mcp(Trainingstage = c(1,-1,0)))
summary(contrast)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: aov(formula = Sprachtest ~ Trainingstage, data = daten_training)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 == 0 -2.600 1.018 -2.554 0.0166 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)Für einen gerichteten Hypothesentest kann das Argument
alternative spezifiziert werden:
contrast <- glht(aov_ergebnis, # im R Chunk oben berechnet
linfct = mcp(Trainingstage = c(1,-1,0)),
alternative = "less") # oder "greater", je nach Richtung der H1
summary(contrast)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: aov(formula = Sprachtest ~ Trainingstage, data = daten_training)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(<t)
## 1 >= 0 -2.600 1.018 -2.554 0.0083 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)Mit glht() können auch direkt mehrere Kontraste getestet
werden. Dazu wird am besten zunächst eine Kontrastmatrix erstellt, die
dann dem Argument linftc übergeben wird. Standardmäßig
führt die Ausgabe mit summary() zur Verwendung einer
Korrektur bzgl. der \(\alpha\)-Kumulation; mit dem Argument
test = adjusted("none") wird dies hier verhindert, sodass
die \(p\)-Werte denen entsprechen, die
auch eine händische Berechnung ergeben würde:
contr.matrix <- rbind("1 vs 2" = c(1, -1, 0),
"1 vs 3" = c(1, 0, -1),
"2 vs 3" = c(0, 1, -1))
contrast2 <- glht(aov_ergebnis,
linfct = mcp(Trainingstage = contr.matrix))
summary(contrast2, test = adjusted("none"))##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: User-defined Contrasts
##
##
## Fit: aov(formula = Sprachtest ~ Trainingstage, data = daten_training)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 vs 2 == 0 -2.600 1.018 -2.554 0.0166 *
## 1 vs 3 == 0 -5.000 1.018 -4.912 3.86e-05 ***
## 2 vs 3 == 0 -2.400 1.018 -2.358 0.0259 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- none method)Bis hierher haben wir Kontraste immer so behandelt, als wären sie a
priori formuliert worden. Im Anschluss an eine signifikante
Varianzanalyse werden aber oft a posteriori Kontraste durchgeführt, um
zu untersuchen, welche \(\mu\)s sich
unterscheiden und für den signifikanten Effekt verantwortlich sind. Dazu
werden oft sämtliche Einzelvergleiche durchgeführt – so wie im letzten
Beispiel – allerdings muss dann eine Anpassung des \(\alpha\)-Niveaus durchgeführt werden. Die
Funktion PostHocTest() im Paket DescTools
bietet eine sehr einfache Variante zur Durchführung solcher Vergleiche.
Im nächsten Beispiel wird die Variante von Scheffé zur Anpassung
verwendet:
library(DescTools)
PostHocTest(aov_ergebnis, method = "scheffe")##
## Posthoc multiple comparisons of means: Scheffe Test
## 95% family-wise confidence level
##
## $Trainingstage
## diff lwr.ci upr.ci pval
## 2-1 2.6 -0.03661653 5.236617 0.05386 .
## 3-1 5.0 2.36338347 7.636617 0.00018 ***
## 3-2 2.4 -0.23661653 5.036617 0.07989 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Durch die Anpassung werden die Vergleiche natürlich konservativer durchgeführt und man sieht, dass nur der Vergleich 1 vs. 3 Trainingstagen signifikant wäre (während im Beispiel vorher alle Vergleiche signifikant waren). Eine etwas weniger konservative Methode ist Tukey’s HSD-Methode, bei der zumindest auch der Vergleich 1 vs.2 Trainingstage signifikant wird:
PostHocTest(aov_ergebnis, method = "hsd")##
## Posthoc multiple comparisons of means : Tukey HSD
## 95% family-wise confidence level
##
## $Trainingstage
## diff lwr.ci upr.ci pval
## 2-1 2.6 0.07598653 5.124013 0.04253 *
## 3-1 5.0 2.47598653 7.524013 0.00011 ***
## 3-2 2.4 -0.12401347 4.924013 0.06472 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.4 Modellgleichung der einfaktoriellen Varianzanalyse
Am Ende diesen Teils betrachten wir die Varianzanalyse noch aus einem etwas anderem Blickwinkel. Dies wird uns im nächsten Teil helfen, bestimmte Konzepte zu verstehen, die wir im nächsten Teil einführen werden.
Die nachfolgende Abbildung beruht wieder auf den Daten des Eingangsbeispiels. In ihr sind als kleine schwarze Punkte die Werte einzelner Versuchspersonen in Abhängigkeit von ihrer Gruppe dargestellt. Als horizontale Linien sind in schwarz die Gruppenmittelwerte (13.1, 15.7 und 18.1) und als rot der Gesamtmittelwert (15.6) eingezeichnet. Die vertikalen, gestrichelten Linien visualisieren Abweichungen:
- Die roten Linien sind die Effekte \(\alpha_j\), also die Abweichungen der Gruppenmittelwerte vom Gesamtmittelwert \(M\). Diese Abweichung wird oft die erklärte Abweichung genannt.
- Die schwarzen Linien sind die Abweichungen der einzelnen Versuchspersonen \(y_{ij}\) von ihrem Gruppenmittelwert \(M_j\).
Die Abweichungen der einzelnen Personen vom Gruppenmittelwert bezeichnen wir als Messfehler \(e_{ij}\): \[e_{ij}=y_{ij}-\mu_j\] Diese Abweichung wird dann oft nicht-erklärte Abweichung genannt, da sie nicht auf die Wirkung des untersuchten Faktors zurückgeht. Für den Messfehler \(e_{ij}\) gibt es eine interessante Verteilung: Zunächst gilt wegen der eingangs gemachten Voraussetzungen der Varianzanalyse für eine Zufallsvariable \(\mathbf{y_{ij}}\), dass auch diese normalverteilt ist: \[\mathbf{y_{ij}} \sim N(\mu_j,\sigma^2)\qquad\forall j\in\{1,\ldots,J\},\forall i\in\{1,\ldots,n_j\}\] Daraus folgt, dass – der Messfehler als Zufallsvariable aufgefasst – auch die Messfehler normalverteilt sind, und zwar mit einem Erwartungswert von Null: \[\mathbf{e_{ij}}\sim N(0,\sigma^2)\] In anderen Worten: Über alle (unendlich vielen) möglichen Realisierungen einer Person \(i\) in Gruppe \(j\) würde sich der Messfehler ausmitteln.
Insgesamt ergibt sich auf der Ebene der einzelnen Messwerte, dass sich eine einzelne Beobachtung \(y_{ij}\) additiv zusammensetzt aus (1) dem Gesamtmittelwert, (2) dem Effekt der Stufe \(j\) und (3) dem Messfehler \(e_{ij}\): \[y_{ij}=\mu+\alpha_j+e_{ij}\] Diese Gleichung bezeichnet man auch als Modellgleichung der einfaktoriellen Varianzanalyse auf Messwertebene. Betrachten wir auf beiden Seiten der Gleichung den Erwartungswert, dann erhalten wir die Modellgleichung auf Parameterebene: \[\mu_j=\mu+\alpha_j\] Auf diese Gleichung werden wir in Teil 4 zurückkommen.
3.5 Literatur
Cohen, J. (1988). Statistical power analysis for the behavioral sciences. (2. Auflage). Erlbaum.
Janczyk, M., & Pfister, R. (2020). Inferenzstatistik verstehen. Von A wie Signifikanztest bis Z wie Konfidenzintervall (3. überarbeitete und erweiterte Auflage). Springer.
Loftus, G. R., & Masson, M. E. J. (1994). Using confidence intervals in within-subject designs. Psychonomic Bulletin & Review, 1, 476–490.
Schad et al. (2020). How to capitalize on a priori contrasts in linear (mixed) models: A tutorial. Journal of Memory and Language, 110, 104038.
4 Mehrfaktorielle Varianzanalyse
In den Teilen 2 und 3 haben wir wichtige Konzepte der Varianzanalyse an ihrem einfachsten Fall – der einfaktoriellen Varianzanalyse – erarbeitet. Nun werden wir Konzepte wie die Quadratsummenzerlegung und den resultierenden \(F\)-Bruch erweitern auf komplexere Fälle. Diese komplexeren Fälle sind eigentlich auch diejenigen, bei denen die Varianzanalyse ihre besonderen Vorteile ausspielen kann.
Für einen solchen Fall betrachten wir wieder das Eingangsbeispiel aus Teil 2. Neben der Leistung im Sprachtest ist in den Daten das Alter der untersuchten Kinder erfasst worden. Wir haben bisher das Alter außer Acht gelassen und festgestellt, dass die Leistung im Sprachtest mit der Anzahl der Trainingstage zunahm. Nun könnte eine Theorie vorhersagen, dass dies besonders stark bei jüngeren Kindern bis einschließlich 7 Jahren der Fall ist, während ältere Kinder weniger stark von den Trainingstagen profitieren.
Schauen wir zunächst die Daten entsprechend in einer Abbildung an.
Dazu laden wir den Datensatz wieder und kodieren
Altersgruppe entsprechend als neue Variable mit den
Ausprägungen jünger und älter:
daten_training <- read.table("./Daten/daten_training.DAT",
header = TRUE)
# Kodieren nach Alter
daten_training$Altersgruppe <- "jünger"
daten_training$Altersgruppe[daten_training$Alter >= 8] <- "älter"
# jeweils 5 Fälle pro entstehender Gruppe Trainingstage x Altersgruppe
table(daten_training$Trainingstage, daten_training$Altersgruppe)##
## älter jünger
## 1 5 5
## 2 5 5
## 3 5 5head(daten_training,
n = 2) # die ersten zwei Fälle anzeigenDann können wir die Daten getrennt für beide Altersgruppen visualisieren und wir sehen, dass die formulierte Hypothese nicht abwegig ist. Dieses Szenario ist ein typischer Fall für eine sog. mehrfaktorielle Varianzanalyse, hier genauer: eine zweifaktorielle Varianzanalyse: Neben dem dreistufigen Faktor “Trainingsgruppe” gibt es nun zusätzlich den zweistufigen Faktor “Altersgruppe”
4.1 Grundlagen und Begriffe der zweifaktoriellen Varianzanalyse
Die Prinzipien erarbeiten wir uns an einer zweifaktoriellen Varianzanalyse; sie können aber danach leicht auf mehr Faktoren verallgemeinert werden.
4.1.1 Annahmen
Wir bezeichnen im Folgenden den ersten Faktor nun als Faktor \(A\) mit \(J\) Stufen und den zweiten Faktor als Faktor \(B\) mit \(K\) Stufen. Alle \(J\cdot K\) Kombinationen beider Faktoren werden dann an unabhängigen Stichproben realisiert. Bezogen auf das Beispiel ist “Trainingstage” der Faktor \(A\) mit 3 Stufen und “Altersgruppe” ist der Faktor \(B\) mit 2 Stufen. Es gibt also \(3\cdot 2=6\) Kombinationen bzw. Gruppen.
Ganz ähnlich wie in Teil 2 bereits nehmen wir an, dass die abhängige Variable in der jeweiligen Population normalverteilt ist, wobei wir für jede Population einen eigenen Erwartungswert zulassen, aber die gleiche Varianz annehmen: \[Y_{jk}\sim N(\mu_{jk},\sigma^2)\qquad \forall j\in\{1,\ldots,J\}, \forall k\in\{1,\ldots,K\}\] Darüber hinaus gehen wir auch von Intervallskalenniveau der abhängigen Variablen und Unabhängigkeit der einzelnen Stichproben aus. Wir machen hier weiter die Einschränkung, dass die einzelnen Gruppen alle gleich groß sind. Diese Einschränkung macht die Formeln für den hier betrachteten Ansatz etwas einfacher. Moderne Statistikprogramme benötigen diese Einschränkung nicht, weisen aber durchaus darauf hin, wenn die Gruppen unterschiedlich groß sind, weil dies durchaus Probleme mit sich bringen kann. Daher ist es bei der Planung einer Studie sinnvoll, von vornherein gleich große Gruppen zu planen.
4.1.2 Globale Erwartungswerte
Eine wichtige Rolle spielen im Folgenden die sog. globalen
Erwartungswerte \(\mu_{j\cdot}\) und \(\mu_{\cdot k}\), die für jede Stufe des
einen Faktors über alle Stufen des anderen Faktors hinweg berechnet
werden. Der Punkt ersetzt dabei den Index des Faktors, über dessen
Stufen hinweg gemittelt wird. Für den Faktor Trainingstage unseres
Beispiels lassen sich die globalen Erwartungswerte also berechnen als
\[\begin{equation*}
\mu_{1\cdot}=\frac{\mu_{11}+\mu_{12}}{2}
\quad\text{und}\quad \mu_{2\cdot}=\frac{\mu_{21}+\mu_{22}}{2}
\quad\text{und}\quad \mu_{3\cdot}=\frac{\mu_{31}+\mu_{32}}{2}\text{,}
\end{equation*}\] sowie für den Faktor Altersgruppe als \[\begin{equation*}
\mu_{\cdot 1}=\frac{\mu_{11}+\mu_{21}+\mu_{31}}{3} \quad\text{und}\quad
\mu_{\cdot 2}=\frac{\mu_{12}+\mu_{22}+\mu_{32}}{3} \text{.}
\end{equation*}\] Die folgende Tabelle fasst die (globalen)
Erwartungswerte für unser Beispiel zusammen:
4.1.3 Haupteffekte
Zunächst konzentrieren wir uns auf die sog. Haupteffekte. Im Beispiel könnte man zunächst eine einfaktorielle Varianzanalyse mit dem Faktor Trainingstage durchführen und dabei den Faktor Altersgruppe außer Acht lassen. Eine zweite Möglichkeit wäre eine solche Varianzanalyse mit dem Faktor Altersgruppe, wobei der Faktor Schlafentzug außer Acht gelassen wird. Ähnlich wird auch im Rahmen der Haupteffekte der zweifaktoriellen Varianzanalyse vorgegangen, allerdings steigt die Power durch die Berücksichtigung des jeweils anderen Faktors (dazu später mehr).
In der varianzanalytischen Terminologie wird mit Haupteffekten die Unterschiedlichkeit der globalen Erwartungswerte der einzelnen Stufen jedes Faktors getestet (in gewisser Weise sind die globalen Erwartungswerte die Gruppenerwartungswerte, wenn es den anderen Faktor nicht geben würde). Es gibt dabei immer so viele Haupteffekthypothesen, wie es Faktoren gibt. Ganz ähnlich zum Hypothesenpaar in der einfaktoriellen Varianzanalyse lauten sie in unserem Beispiel:
- Haupteffekt \(\boldsymbol{A}\) (Trainingstage): \[\begin{equation*} \begin{aligned} H_0^A\colon\quad& \mu_{1\cdot}=\mu_{2\cdot}=\mu_{3\cdot} \phantom{\text{sadasdasdasdahahahahaahaaajajaja}}&\phantom{at}\\ H_1^A\colon\quad& \mu_{m\cdot}\neq\mu_{n\cdot}\text{ für mindestens ein Paar }m,n\in \{1,2,3\}&\\ \end{aligned} \end{equation*}\]
- Haupteffekt \(\boldsymbol{B}\) (Altersgruppe): \[\begin{equation*} \begin{aligned} H_0^B\colon\quad& \mu_{\cdot 1}=\mu_{\cdot 2} &\phantom{\text{saadasdasdasdahahahahaahaaajajajajajajaja}}\\ H_1^B\colon\quad& \mu_{\cdot 1}\neq\mu_{\cdot 2}& \end{aligned} \end{equation*}\]
Wichtig ist: Die Haupteffekte sagen uns noch nichts darüber aus, ob – bezogen auf das Beispiel – der Effekt des Schlafentzugs unterschiedliche Wirkung bei jüngeren bzw. älteren Versuchspersonen zeigt. Vielmehr beschreiben Haupteffekte, wie sich Schlafentzug bzw. Altersgruppe generell auf die Leistung auswirken, würde es den anderen Faktor nicht geben.
Als nächsten Schritt führen wir – analog zur einfaktoriellen Varianzanalyse (vgl. Abschnitt 3.2.1) – Effekte ein. Die Effekte des Faktors \(A\) nennen wir wieder \(\alpha_j\) und die Effekte des Faktors \(B\) nennen wir \(\beta_k\): \[\begin{equation} \alpha_j=\mu_{j\cdot}-\mu\quad\forall j\in\{1,2,3\} \quad\text{und}\quad\beta_k=\mu_{\cdot k}-\mu\quad\forall k\in\{1,2\} \end{equation}\] Nehmen wir nun zunächst an, die empirischen Mittelwerte entsprächen den Populationsmittelwerten \(\mu_{jk}\). Diese Situation und die dann resultierenden globalen Erwartungswerte sind in der folgenden Tabelle dargestellt:
Nun können wir die jeweiligen Effekte berechnen, indem der Gesamterwartungswert von den globalen Erwartungswerten subtrahiert wird. Die jeweiligen Effekte addieren sich dabei (hier annähernd) zu 0 auf (dass sie dies hier nicht genau tun, liegt an Rundungsfehlern): \[ \quad\alpha_1 =13.10-15.63=-2.53\quad\alpha_2 =15.70-15.63=0.07 \quad\alpha_3 = 18.1-15.63=2.47 \] und \[\beta_1 =16.80-15.63=1.17\quad\beta_2 =14.47-15.63=-1.16\] Erinnern wir uns nun an 3.4. Dort hatten wir die Modellgleichung der einfaktoriellen Varianzanalyse auf Messwertebene als \[y_{ij}=\mu+\alpha_j+e_{ij}\] und auf Parameterebene als \[\mu_j=\mu+\alpha_j\] beschrieben und gesagt, der Gruppenerwartungswert setze sich additiv zusammen aus Gesamterwartungswert und Effekt der Stufe \(j\) bzw. der Messwert setze sich additiv zusammen aus Gesamterwartungswert, Effekt der Stufe \(j\) und dem Messfehler.
Auf Basis der Effekte \(\alpha_j\) und \(\beta_k\) formulieren wir nun eine ganz ähnliche vorläufige Modellgleichung. Die Gleichung dieses “Haupteffektmodells” lautet: \[\begin{equation} \mu_{jk}=\mu+\alpha_j+\beta_k \end{equation}\] Nun nehmen wir das Haupteffektmodell und berechnen die \(\mu_{jk}\) mit Hilfe des Gesamterwartungswertes und der Effekte. Für \(\mu_{11}\) ergibt sich also z.B. \[\mu_{11}=\mu+\alpha_1+\beta_1=15.63+(-2.53)+1.17=14.27\] Die vom Haupteffektmodell vorhergesagten Werte sehen dann wie folgt aus:
Vergleichen wir diese Werte mit den Werten der letzten Tabelle, wird klar, dass die empirischen und die vom Haupteffektmodell vorhergesagten Werte nicht übereinstimmen. Das Haupteffektmodell erklärt die Daten also nicht ausreichend. Daher werden wir es im folgenden Abschnitt um sog. Interaktionseffekte erweitern. Vorher schauen wir uns die Daten noch visualisiert an und stellen fest, dass beide Linien parallel zueinander verlaufen. Darauf werden wir noch zurückkommen:
4.1.4 Interaktionen
Die folgende Tabelle stellt die Abweichungen der vom Haupteffektmodell vorhergesagten Erwartungswerte von den tatsächlichen Erwartungswerten dar, also \[\mu_{jk}-(\mu+\alpha_j+\beta_k)\] dar. Klar ist, dass in jeder der sechs Zellen Abweichungen vorliegen:
Würden keine Abweichungen vorliegen, also jede Zelle den Wert 0 besitzen, dann würde das rein additive Haupteffektmodell bereits zur Erklärung der Daten ausreichen. Da es allerdings Abweichungen gibt, reicht diese reine Additivität nicht aus und diese Abweichungen sind das, was man Interaktion (hier der zwei Faktoren \(A\) und \(B\)) nennt. Wir bezeichnen die Interaktion für jede Zelle mit \((\alpha\beta)_{jk}\) und berechnen sie also als: \[(\alpha\beta)_{jk}=\mu_{jk}-(\mu+\alpha_j+\beta_k)=\mu_{jk}-\mu-\alpha_j-\beta_k\] Insgesamt sind sog. Interaktionseffekte nichts anderes als Abweichungen von der reinen Additivität der Haupteffekte. Berücksichtigen wir die Interaktionseffekte, können wir eine Modellgleichung formulieren, die immer gültig ist. Auf der Parameterebene lautet sie dann \[\mu_{jk}=\mu+\alpha_j+\beta_k+(\alpha\beta)_{jk}\] bzw. auf der Ebene einzelner Messwerte \[y_{ijk}=\mu+\alpha_j+\beta_k+(\alpha\beta)_{jk}+e_{ijk}\]
Zusätzlich zu den beiden Haupteffekten, können wir in der zweifaktoriellen Varianzanalyse also noch eine Interaktionshypothese testen:
- Interaktionseffekt \(\mathbf{AB}\) (Trainigstage \(\times\) Altersgruppe): \[\begin{equation*} \begin{aligned} H_0^{AB}\colon\quad& (\alpha\beta)_{jk}=0 \quad\forall j\in\{1,2,3\},\forall k\in\{1,2\}\phantom{atatatatatatatatat}\\ H_1^{AB}\colon\quad& \text{es gibt mindestens ein }(\alpha\beta)_{jk}\neq 0\quad j\in\{1,2,3\}, k\in\{1,2\} \end{aligned} \end{equation*}\]
In anderen Worten: Die Nullhypothese des Interaktionseffekts \(AB\) nimmt an, dass es keine Abweichung von reiner Additivität gibt, d.h. die Daten rein durch die Haupteffekte erklärt werden können. Nun werden wir nie einen realen Datensatz vorfinden, in dem die empirischen Interaktionen komplett verschwinden (also 0 werden), auch wenn dies in der Population tatsächlich zutreffen würde. Genau wie bei jedem anderen inferenzstatistischen Verfahren nehmen wir zum Test der Interaktionseffekthypothese die \(H_0\) an. Die Frage ist dann: Wie wahrscheinlich sind die beobachteten oder extremere Abweichungen (also die Interaktionen) unter dieser Annahme? Sind sie sehr unwahrscheinlich, zweifeln wir an der Annahme der \(H_0\) und gehen davon aus, dass eine Interaktion beider Faktoren vorliegt, also eine Abweichung von reiner Additivität.
Da die bisher durchgeführten Berechnungen auf den Beispieldaten basieren, liegt die Vermutung nahe, dass es eine Interaktion zwischen Schlafentzug und Altersgruppe gibt. Betrachtet man dann noch einmal die Visualisierung der Daten vom Anfang, fällt auf, dass die beiden Linien für jüngere und ältere Versuchspersonen nicht parallel zueinander verlaufen: Dies ist in der Tat ein Hinweis auf das Vorliegen einer Interaktion. Wären beide Linien parallel zueinander, würde man keine Interaktion vermuten, sondern sagen, beide Effekte verhalten sich additiv zueinander.
Inhaltlich lässt sich das Vorliegen einer Interaktion wie folgt interpretieren: Die Wirkung des einen Faktors unterscheidet sich je nachdem, welche Stufe des anderen Faktors man betrachtet. Diese Betrachtung kann in beide Richtungen gehen:
- Einerseits sieht man, dass die Wirkung der Trainingstage unterschiedlich stark für die beiden Altersgruppen ausfällt (stärker bei jungen Kindern als bei älteren Kindern)
- Andererseits sieht man, dass der Unterschied beider Altersgruppe mit zunehender Anzahl an Trainingstagen abnimmt: Während nach nur einem Trainingstag noch ein großer Unterschied besteht, sind die Mittelwerte nach drei Trainingstagen fast gleich.
Wären beide Linien parallel, ist die Wirkung des einen Faktors für alle Stufen des anderen Faktors gleich stark. Wichtig: Eine signifikante Zweifach-Interaktion bedeutet nicht, dass unter einer Stufe des einen Faktors ein signifikanter Effekt des anderen Faktors vorliegt, und unter anderen Stufen nicht. Es reicht, wenn die einzelnen Effekte unter den verschiedenen Stufen unterschiedlich stark sind. Dies impliziert übrigens auch: Um festzustellen, ob der Effekt eines Faktors unterschiedlich stark ist bei verschiedenen Stufen eines anderen Faktors, muss die Interaktion berechnet werden. Einzelne Vergleiche reichen nicht; selbst dann nicht, wenn ein Vergleich signifikant wird, ein anderer nicht.
4.2 Rechnerische Durchführung
Nachdem wir die Idee der zweifaktoriellen Varianzanalyse skizziert haben, kommen wir nun zu ihrer eigentlichen Berechnung. Das Vorgehen dabei ist konzeptuell identisch mit dem der einfaktoriellen Varianzanalyse: Zunächst erfassen wir die Gesamtvariation der Daten mit der totalen Quadratsumme \(SS_{\text{tot}}\) und zerlegen diese in einzelne Bestandteile. Diese Bestandteile überführen wir in die Mittleren Quadratsummen, und aus diesen berechnen wir schließlich \(F\)-Brüche, die wir zur Entscheidungsfindung heranziehen. Wir beginnen hier mit dem allgemeinen Fall und zeigen danach die einzelnen Schritte an den Daten des Beispiels, bevor die Auswertung mit R demonstriert wird.
4.2.1 Globale Erwartungswerte und Hypothesen
Im allgemeinen Fall hat der Faktor \(A\) \(J\)
Stufen und der Faktor \(B\) hat \(K\) Stufen. Damit ergibt sich ein
Erwartungswertschema wie nachfolgend dargestellt.
Die globalen Erwartungswerte berechnen sich demnach als \[\begin{equation*} \mu_{j\cdot}=\frac{1}{K}\sum_{k=1}^K\mu_{jk}\quad\forall j\in\{1,\ldots,J\}\quad\text{und}\quad\mu_{\cdot k}=\frac{1}{J}\sum_{j=1}^J\mu_{jk}\quad\forall k\in\{1,\ldots,K\}\text{,} \end{equation*}\] während sich der Gesamterwartungswert \(\mu\) berechnet als \[\begin{equation*} \mu=\frac{1}{JK}\sum_{j=1}^J\sum_{k=1}^K\mu_{jk} \quad\text{bzw.}\quad\mu=\frac{1}{J}\sum_{j=1}^J\mu_{j\cdot}=\frac{1}{K}\sum_{k=1}^K\mu_{\cdot k}\text{.} \end{equation*}\] Der Gesamterwartungswert ergibt sich also entweder als Mittelwert aller Erwartungswerte \(\mu_{jk}\) oder als Mittelwert der globalen Erwartungswerte \(\mu_{\cdot k}\) bzw. \(\mu_{j\cdot}\).
Die zu testenden Hypothesen im allgemeinen Fall lauten dann:
- Haupteffekt \(\boldsymbol{A}\):
- \(H_0^A\colon \quad \mu_{1\cdot}=\mu_{2\cdot}=\ldots=\mu_{J\cdot}\)
- \(H_1^A\colon \quad \mu_{m\cdot}\neq\mu_{n\cdot}\text{ für mindestens ein Paar }m,n\in \{1,\cdots,J\}\)
- Haupteffekt \(\boldsymbol{B}\):
- \(H_0^B\colon \quad \mu_{\cdot 1}=\mu_{\cdot 2}=\ldots=\mu_{\cdot K}\)
- \(H_1^B\colon \quad \mu_{\cdot m}\neq\mu_{\cdot n}\text{ für mindestens ein Paar }m,n\in \{1,\cdots,K\}\)
- Interaktionseffekt \(\boldsymbol{AB}\):
- \(H_0^{AB}\colon \quad (\alpha\beta)_{jk}=0\quad\forall j\in\{1,\ldots,J\},\forall k\in\{1,\ldots,K\}\)
- \(H_1^{AB}\colon \quad \text{es gibt mindestens ein }(\alpha\beta)_{jk}\neq 0\quad j\in\{1,\ldots,J\}, k\in\{1,\ldots,K\}\)
4.2.2 Quadratsummenzerlegung
In Teil 2 sind wir zum ersten Mal im Kontext der einfaktoriellen Varianzanalyse mit dem Begriff der Quadratsummen- bzw. Varianzzerlegung konfrontiert worden: Eine Art Gesamtvarianz, nämlich die \(SS_{\text{tot}}\), wurde hierbei als Summe zweier anderer varianzähnlicher Terme \(SS_A\) und \(SS_w\) ausgedrückt: \[SS_{\text{tot}}=SS_A+SS_w\].
Genau das Gleiche passiert auch im Fall der zweifaktoriellen
Varianzanalyse. Allerdings teilen wir \(SS_{\text{tot}}\) jetzt in insgesamt vier
Bestandteile auf: diejenigen Anteile, die auf die Haupteffekte der
Faktoren \(A\) und \(B\) (\(SS_A\) bzw. \(SS_B\)) sowie auf deren Interaktion \(AB\) (\(SS_{AB}\)) zurückgehen, und schließlich den
Anteil, der auf den Messfehler (\(SS_w\)) zurückgeht. Auch hier gilt die
gleiche Additivität: \[SS_{\text{tot}}=SS_A+SS_B+SS_{AB}+SS_w\]
Genau wie in der einfaktoriellen Varianzanalyse wird \(SS_w\) die Grundlage für den Fehlerterm
bilden. Die Varianzzerlegungen der ein- und der zweifaktoriellen
Varianzanalyse sind in der folgenden Abbildung einander gegenüber
gestellt:
Im Folgenden wird nun die konkrete Berechnung der Quadratsummen beschrieben. Dabei meinen wir mit \(y_{ijk}\) den Messwert von Versuchsperson \(i\) in der Gruppe, die sich durch die Kombination von Stufe \(j\) des Faktors \(A\) und Stufe \(k\) des Faktors \(B\) ergibt. Bei gleichem Stichprobenumfang \(n\) sind die Mittelwerte \(M_{jk}\) dann: \[\begin{gather*} M_{jk}=\frac{1}{n}\sum_{i=1}^ny_{ijk}\qquad \forall j\in\{1,\ldots,J\}, \forall k\in\{1,\ldots,K\}\text{.} \end{gather*}\] Als globale Mittelwerte der einzelnen Faktorstufen über die Stufen des jeweils anderen Faktors hinweg ergeben sich: \[\begin{gather*} \begin{aligned} M_{j\cdot}&=\frac{1}{K}\sum_{k=1}^KM_{jk}=\frac{1}{nK}\sum_{k=1}^K\sum_{i=1}^ny_{ijk} \qquad \forall j\in\{1,\ldots,J\}\\ \end{aligned} \end{gather*}\] und \[\begin{gather*} \begin{aligned} M_{\cdot k}&=\frac{1}{J}\sum_{j=1}^JM_{jk}=\frac{1}{nJ}\sum_{j=1}^J\sum_{i=1}^ny_{ijk} \qquad \forall k\in\{1,\ldots,K\}\text{.}\\ \end{aligned} \end{gather*}\] Schließlich kann der Gesamtmittelwert der Daten auf verschiedene Arten berechnet werden, z.B. als Mittelwert der (globalen) Mittelwerte \[\begin{gather*} M=\frac{1}{JK}\sum_{j=1}^J\sum_{k=1}^KM_{jk}\quad\text{bzw.}\quad M= \frac{1}{J}\sum_{j=1}^JM_{j\cdot}=\frac{1}{K}\sum_{k=1}^KM_{\cdot k} \end{gather*}\] oder direkt aus den Rohdaten: \[\begin{gather*} M=\frac{1}{nJK}\sum_{j=1}^J\sum_{k=1}^K\sum_{i=1}^n y_{ijk}\text{.} \end{gather*}\] Nach diesen Vorüberlegungen können wir nun die entsprechenden Quadratsummen berechnen. Diese Quadratsummen verfügen auch wiederum über bestimmte Freiheitsgrade, die hier direkt mit angegeben werden. Dabei meint \(N\) die Gesamtanzahl aller Datenpunkte, also \(N=nJK\): \[\begin{alignat*}{2} SS_{\text{tot}}&=\sum_{j=1}^J\sum_{k=1}^K\sum_{i=1}^n(y_{ijk}-M)^2\quad& &\text{und }\quad df_{\text{tot}}=N-1\\ SS_A&=nK\sum_{j=1}^J(M_{j\cdot}-M)^2\quad& &\text{und }\quad df_A=J-1\\ SS_B&=nJ\sum_{k=1}^K(M_{\cdot k}-M)^2\quad& &\text{und }\quad df_B=K-1\\ SS_{AB}&=n\sum_{j=1}^J\sum_{k=1}^K(M_{jk}-M_{j\cdot}-M_{\cdot k}+M)^2\quad& &\text{und }\quad df_{AB}=(J-1)(K-1)\\ SS_w&=\sum_{j=1}^J\sum_{k=1}^K\sum_{i=1}^n(y_{ijk}-M_{jk})^2\quad& &\\ &=\sum_{j=1}^J\sum_{k=1}^K n \cdot S_{jk}^2\quad & &\text{und }\quad df_w=JK(n-1)=N-JK\text{.} \end{alignat*}\] Wie oben bereits bemerkt, summieren sich die Quadratsummen der Effekte und die Fehlerquadratsumme \(SS_w\) zu \(SS_{\text{tot}}\) auf: \[\begin{gather*}\label{varianzzweifak} SS_{\text{tot}}=SS_A+SS_B+SS_{AB}+SS_w\text{.} \end{gather*}\]
4.2.3 Mittlere Quadratsummen
Aus den Quadratsummen berechnen wir nun im nächsten Schritt die entsprechenden Mittleren Quadratsummen, indem die Quadratsummen durch ihre jeweiligen Freiheitsgrade dividiert werden. Die Mittleren Quadratsummen der drei Effekte lauten daher \[\begin{gather*} MS_A=\frac{SS_A}{J-1}\text{,}\qquad MS_B=\frac{SS_B}{K-1}\qquad \text{und}\qquad MS_{AB}=\frac{SS_{AB}}{(J-1)(K-1)}\text{,} \end{gather*}\] und die Mittlere Quadratsumme des Fehlerterms ist: \[\begin{gather*} MS_w=\frac{SS_w}{N-JK}\text{.} \end{gather*}\]
4.2.4 \(F\)-Brüche
Im letzten Schritt werden – wiederum genau wie in der einfaktoriellen Varianzanalyse – aus den Mittleren Quadratsummen entsprechende \(F\)-Brüche gebildet. Zur Erinnerung: Bei einer zweifaktoriellen Varianzanalyse können wir drei Hypothesen testen: Die beiden Haupteffekthypothesen \(A\) und \(B\) sowie die Interaktions-hypothese \(AB\). Pro Hypothesenpaar gibt es nun jeweils einen \(F\)-Bruch, in dessen Zähler die Mittlere Quadratsumme des relevanten Effekts steht. Im Nenner steht dabei immer \(MS_w\); wir verwenden also immer den gleichen Fehlerterm.
Diese \(F\)-Brüche erfüllen wieder genau die zwei Eigenschaften, die jede Prüfgröße haben sollte:
- Sie nehmen besonders extreme Werte an, je stärker die Daten gegen die \(H_0\) sprechen.
- Ihre Verteilung unter der Annahme der Gültigkeit der entsprechenden \(H_0\) ist bekannt; alle drei \(F\)-Brüche sind \(F\)-verteilt, i.d.R. aber mit unterschiedlichen Zählerfreiheitsgraden (die Nennerfreiheitsgrade sind immer gleich, da in allen drei Fällen \(MS_w\) den Fehlerterm bildet): \[\begin{gather*} \begin{aligned} F^A&=\frac{MS_A}{MS_w}\overset{H_0^A}{\sim}F_{J-1,N-JK}\\ F^B&=\frac{MS_B}{MS_w}\overset{H_0^B}{\sim}F_{K-1,N-JK}\\ F^{AB}&=\frac{MS_{AB}}{MS_w}\overset{H_0^{AB}}{\sim}F_{(J-1)(K-1),N-JK}\text{.} \end{aligned} \end{gather*}\] Die Entscheidungsregeln für die drei Hypothesen sind ganz analog zur Entscheidungsregel bei der einfaktoriellen Varianzanalyse:
Verwirf die \(H_0\), falls \(F\geq F_{\text{krit}}\) bzw. falls \(p\leq\alpha\) ist.
4.2.5 Beispiel
4.2.5.1 Daten
Die Daten des Beispiels sind in der nachfolgenden Tabelle noch einmal
dargestellt, sowohl die Werte der einzelnen Personen als auch die
Gruppenmittelwerte und -varianzen:
Wir fassen dabei “Trainingstage” als den Faktor \(A\) mit \(J=3\) Stufen auf. Die “Altersgruppe” wird als Faktor \(B\) mit \(K=2\) Stufen aufgefasst.
4.2.5.2 Vorbereitende Schritte
Die Gruppenmittelwerte sind in der Tabelle bereits enthalten; wir berechnen daher nun den Gesamtmittelwert sowie die globalen Mittelwerte einzelner Faktorstufen. Zur Berechnung des Gesamtmittelwerts \(M\) können wir bspw. auf die Rohdaten zurückgreifen: \[ \begin{split} M&=\frac{1}{nJK}\sum_{j=1}^J\sum_{k=1}^K\sum_{i=1}^n y_{ijk} =\frac{1}{5\cdot 3\cdot 2}\cdot\sum_{j=1}^3\sum_{k=1}^2\sum_{i=1}^5 y_{ijk}\\ &=\frac{1}{30}\cdot (12+11+10+8+\ldots+19+18+14+18)=\frac{1}{30}\cdot 469\\ &= 15.633 \end{split} \] Nun folgen die globalen Mittelwerte. Wir beginnen dabei mit den globalen Mittelwerten für die drei Stufen des Faktors \(A\) (“Schlafentzug”) und unter Rückgriff auf die Rohdaten lautet die allgemeine Formel: \[ M_{j\cdot}=\frac{1}{nK}\sum_{k=1}^K\sum_{i=1}^n y_{ijk}\qquad \forall j\in \{1,\ldots,J\} \] Damit ergeben sich im konkreten Falle unseres Beispiels folgende Werte: \[\begin{align*} M_{1\cdot}&=\frac{1}{5\cdot 2}\sum_{k=1}^2\sum_{i=1}^5 y_{i1k}=\frac{1}{10}\cdot (12+11+10+8+13+17+16+15+13+16)=13.1\\ M_{2\cdot}&=\frac{1}{5\cdot 2}\sum_{k=1}^2\sum_{i=1}^5 y_{i2k}=\frac{1}{10}\cdot (15+15+14+16+13+18+17+16+16+17)=15.7\\ M_{3\cdot}&=\frac{1}{5\cdot 2}\sum_{k=1}^2\sum_{i=1}^5 y_{i3k}=\frac{1}{10}\cdot (17+18+16+19+20+22+19+18+14+18)=18.1 \end{align*}\] Für die globalen Mittelwerte des Faktors \(B\) (“Altersgruppe”) benutzen wir die allgemeine Formel \[ M_{\cdot k}=\frac{1}{nJ}\sum_{j=1}^J\sum_{i=1}^n y_{ijk}\qquad \forall k\in \{1,\ldots,K\}\text{ ,} \] und erhalten so: \[\begin{align*} M_{\cdot 1}&=\frac{1}{5\cdot 3}\cdot\sum_{j=1}^3\sum_{i=1}^5 y_{ij1} =\frac{1}{15}\cdot(12+11+10+\ldots+16+19+20)=14.467\\ M_{\cdot 2}&=\frac{1}{5\cdot 3}\cdot\sum_{j=1}^3\sum_{i=1}^5 y_{ij2} =\frac{1}{15}\cdot(17+16+15+\ldots+18+14+18)=16.8 \end{align*}\]
Zusammengefasst ergibt sich also folgendes Schema aus den Mittelwerten:
4.2.5.3 Quadratsummen
Nach diesen Vorbereitungen berechnen wir nun die Quadratsummen. Der Vollständigkeit halber beginnen wir dabei mit der Totalen Quadratsumme \(SS_{\text{tot}}\): \[ \begin{split} SS_{tot}&=\sum_{j=1}^J\sum_{k=1}^K\sum_{i=1}^n (y_{ijk}-M)^2= \sum_{j=1}^3\sum_{k=1}^2\sum_{i=1}^5 (y_{ijk}-15.633)^2\\ &=(12-15.633)^2+(11-15.633)^2+(10-15.633)^2+\ldots+\\ &\phantom{=|}(18-15.633)^2+(14-15.633)^2+(18-15.633)^2\\ &=264.9667 \end{split} \] Als Quadratsummen der Faktoren \(A\) und \(B\) ergeben sich \[ \begin{split} SS_A&=nK\sum_{j=1}^J(M_{j\cdot}-M)^2=5\cdot 2\cdot\sum_{j=1}^3(M_{j\cdot}-15.633)^2\\ &=10\cdot [(13.1-15.633)^2+(15.7-15.633)^2+(18.1-15.633)^2]\\ &=125.0667 \\ SS_B&=nJ\sum_{k=1}^K(M_{\cdot k}-M)^2=5\cdot 3\cdot\sum_{k=1}^2(M_{\cdot k}-15.633)^2\\ &=15\cdot [(16.8-15.633)^2+(14.467-15.633)^2]\\ &=40.82167 \end{split} \] und die Quadratsumme der Interaktion ergibt \[ \begin{split} SS_{AB}&=n\sum_{j=1}^J\sum_{k=1}^K(M_{jk}-M_{j\cdot}-M_{\cdot k}+M)^2\\ &=5\cdot [(M_{11}-M_{1\cdot}-M_{\cdot 1}+M)^2+(M_{12}-M_{1\cdot}-M_{\cdot 2}+M)^2+\\ &\phantom{=|}(M_{21}-M_{2\cdot}-M_{\cdot 1}+M)^2+(M_{22}-M_{2\cdot}-M_{\cdot 2}+M)^2+\\ &\phantom{=|}(M_{31}-M_{3\cdot}-M_{\cdot 1}+M)^2+(M_{32}-M_{3\cdot}-M_{\cdot 2}+M)^2]\\ &=5\cdot [(15.4-13.1-16.8+15.633)^2+(10.8-13.1-14.467+15.633)^2+\\ &\phantom{=|}(16.8-15.7-16.8+15.633)^2+(14.6-15.7-14.467+15.633)^2+\\ &\phantom{=|}((18.2-18.1-16.8+15.633)^2+(18.0-18.1-14.467+15.633)^2]\\ &=24.26668 \end{split} \] Schließlich benötigen wir noch die Fehlerquadratsumme: \[ \begin{split} SS_w&=\sum_{j=1}^J\sum_{k=1}^K\sum_{i=1}^n(y_{ijk}-M_{jk})^2= \sum_{j=1}^3\sum_{k=1}^2\sum_{i=1}^5(y_{ijk}-M_{jk})^2\\ &=(12-10.8)^2+(11-10.8)^2+\ldots+(14-18.2)^2+(18-18.2)^2\\ &=74.80 \end{split} \]
4.2.5.4 Mittlere Quadratsummen
Indem die Quadratsummen durch ihre Freiheitsgrade dividiert werden, erhalten wir in einem nächsten Schritt die Mittleren Quadratsummen: \[ \begin{split} MS_A&=\frac{SS_A}{J-1}=\frac{125.0667}{3-1}=62.53335\\ MS_B&=\frac{SS_B}{K-1}=\frac{40.82167}{2-1}=40.82167\\ MS_{AB}&=\frac{SS_{AB}}{(J-1)(K-1)}=\frac{24.26668}{(3-1)(2-1)}=12.13334\\ MS_w&=\frac{SS_w}{JK(n-1)}=\frac{74.8}{3\cdot 2\cdot (5-1)}=3.116667\\ \end{split} \]
4.2.5.5 Die \(F\)-Brüche
Der letzte Schritt besteht aus der Berechnung der drei \(F\)-Brüche. Dazu werden die Mittleren Quadratsumme der Faktoren durch die Mittlere Fehlerquadratsumme dividiert: \[ \begin{split} F^A&=\frac{MS_A}{MS_w}=\frac{62.53335}{3.116667}=20.06417\\ F^B&=\frac{MS_B}{MS_w}=\frac{40.82167}{3.116667}=13.09786\\ F^{AB}&=\frac{MS_{AB}}{MS_w}=\frac{12.13334}{3.116667}=3.89305\\ \end{split} \] Diese empirischen \(F\)-Brüche können dann zur Evaluierung der Hypothesen herangezogen werden. Dazu gelten die gleichen Entscheidungsregeln, die wir am Beispiel der einfaktoriellen Varianzanalyse eingeführt haben.
4.3 Weiterführende Anmerkungen
4.3.1 Effektstärken
Als ein wichtiges Maß der Effektstärke für Varianzanalysen hatten wir den Anteil der Effektvarianz an der Gesamtvarianz eingeführt: \[\begin{gather*} \eta^2=\frac{\sigma^2_{\alpha}}{\sigma^2_{\text{tot}}} \end{gather*}\] Im Fall mehrerer Faktoren setzt sich die Gesamtvarianz nun zusammen aus den Varianzen aller beteiligten Effekte sowie der Fehlervarianz. Daher wird hier oft auf ein anderes Maß zurückgegriffen, das als partielles \(\mathbf{\eta^2}\) bezeichnet wird. Hierbei wird die Effektvarianz nicht mehr an der Gesamtvarianz, sondern an der Summe aus relevanter Effektvarianz und Fehlervarianz relativiert: \[\begin{gather*} \eta^2_p=\frac{\sigma^2_{\text{Effekt}}}{\sigma^2_{\text{Effekt}}+\sigma^2} \end{gather*}\] Da der Nenner für \(\eta^2_p\) kleiner ist als im Fall von \(\eta^2\), gilt i.d.R. \(\eta^2_p>\eta^2\). Aus diesem Grund kann die Summe der \(\eta^2_p\) aller Effekte auch größer als 1 werden; für \(\eta^2\) ist dies nicht möglich. Zur Powerberechnung wird üblicherweise auch \(\eta^2_p\) herangezogen.
Ähnlich wie \(\eta^2\) kann auch \(\eta^2_p\) mit Hilfe der (Mittleren) Quadratsummen geschätzt werden. Die meisten Computerprogramme verwenden dabei folgende Formel: \[\begin{gather*} \hat{\eta}^2_p=\frac{SS_{\text{Effekt}}}{SS_{\text{Effekt}}+SS_w} \end{gather*}\] Alternativ kann auch wieder eine konservativere Schätzung genutzt werden: \[\begin{gather*} \hat{\eta}^2_p=\frac{SS_{\text{Effekt}}-df_{\text{Effekt}}\cdot MS_w} {SS_{\text{Effekt}}+(N-df_{\text{Effekt}})\cdot MS_w} \end{gather*}\]
4.3.2 Ergebnisdarstellung
Eine vollständige Darstellung der Ergebnisse einer zweifaktoriellen Varianzanalyse erfolgt wieder am besten in Form einer Tabelle:
4.3.3 Ein- vs. zweifaktorielle Varianzanalyse: Vorteile mehrfaktorieller Varianzanalysen
Ein erster Vorteil einer zweifaktoriellen Varianzanalyse (im Vergleich zu zwei einfaktoriellen Varianzanalysen) ist, dass zusätzlich zu den zwei Haupteffekthypothesen noch die Interaktion getestet werden kann – also das Zusammenspiel beider Faktoren und unterschiedliche Wirkungen eines Faktors in Abhängigkeit von der betrachteten Stufe des anderen Faktors.
Darüber hinaus gibt es noch einen anderen Vorteil der zwei
Haupteffekthypothesen gegenüber zwei einfaktoriellen Varianzanalysen.
Vergleichen Sie hierzu die gerade dargestellte Ergebnistabelle der
zweifaktoriellen Varianzanalyse mit der Ergebnistabelle der
einfaktoriellen Varianzanalyse mit dem Faktor Schlafentzug (vgl. auch
Teil 2):
Sowohl die Quadratsumme \(SS_A\) als
auch die Mittlere Quadratsumme \(MS_A\)
sind identisch. Der \(F\)-Bruch
allerdings ist deutlich größer, er würde also \(F_{\text{krit}}\) schneller überschreiten
und auch der \(p\)-Wert als die Fläche
rechts vom \(F\)-Wert unter der
relevanten \(F\)-Verteilung wäre
kleiner. Anders ausgedrückt: der Haupteffekt würde leichter signifikant
werden als bei der einfaktoriellen Varianzanalyse, er hätte also eine
größere Power.
Der größere \(F\)-Bruch rührt daher, dass der entsprechende Fehlerterm \(MS_w\) im Nenner des \(F\)-Bruchs im zweifaktoriellen Fall kleiner geworden ist. Da insgesamt also bei Berücksichtigung eines zweiten Faktors (also dem Rechnen einer zweifaktoriellen Varianzanalyse) die Wirkung des ersten Faktors nicht verändert wird, ist im einfaktoriellen Fall der Fehlerterm \(SS_w\) zusammengesetzt aus drei verschiedenen Termen der zweifaktoriellen Varianzanalyse, nämlich \(SS_B+SS_{AB}+SS_w\). Da \(SS_B\) und \(SS_{AB}\) aber immer größer Null sind (man sagt: “Varianz binden”), ist \(SS_w\) im zweifaktoriellen Fall kleiner als im einfaktoriellen Fall.
Durch die so verkleinerte Fehlerquadratsumme wird ein Effekt in einer mehrfaktoriellen Varianzanalyse i.d.R. also eher signifikant, da Teile der Gesamtvariabilität durch den zweiten Faktor bzw. die Interaktion beider Faktoren gebunden werden. Dies können Sie auch daran sehen, dass \(SS_w\) im einfaktoriellen Fall (139.90) genau der Summe \(SS_B+SS_{AB}+SS_w\) der zweifaktoriellen Varianzanalyse entspricht.
4.3.4 Varianzanalysen mit mehr als zwei Faktoren
Bisher haben wir in diesem Kapitel nur zwei Faktoren betrachtet, und eine entsprechende Varianzanalyse erlaubt bereits das Testen dreier Hypothesenpaare: Haupteffekt \(A\), Haupteffekt \(B\) sowie deren Interaktion \(AB\). Selbstverständlich können auch Varianzanalysen mit mehr als zwei Faktoren gerechnet werden, was die Varianzanalyse zu einem sehr flexiblen Werkzeug werden lässt.
Welche Hypothesenpaare können wir in einem dreifaktoriellen Design mit den Faktoren \(A\), \(B\) und \(C\) testen? Zunächst gibt es genau drei Haupteffekthypothesen \(A\), \(B\) und \(C\). Darüber hinaus kann jeder Faktor mit jedem anderen Faktor interagieren, also gibt es drei Interaktionshypothesen 1. Ordnung (d.h. über Zweifach-Interaktionen): \(AB\), \(AC\) und \(BC\). Schließlich – und das ist etwas Neues – gibt es auch eine Interaktion 2. Ordnung, nämlich die Dreifach-Interaktion \(ABC\). Wird diese Interaktion signifikant, bedeutet dies, dass die Interaktionen 1. Ordnung in Abhängigkeit von der Stufe des jeweils dritten Faktors unterschiedlich stark ausfallen. Ein solcher Fall ist in der folgenden Abbildung beispielhaft illustriert: Während unter der ersten Stufe des Faktors \(C\) (linker Teil der Abbildung) keine Interaktion der Faktoren \(A\) und \(B\) vorliegt (die Linien sind parallel!), interagieren diese beiden Faktoren unter der zweiten Stufe des Faktors \(C\) (rechter Teil der Abbildung; keine parallelen Linien). Wichtig ist hier (wir kommen in der Vorlesung Forschungsmethoden darauf zurück): Eine signifikante Dreifach-Interaktion muss nicht heißen, dass im einen Fall eine Zweifach-Interaktion anwesend ist, im anderen nicht – es reicht, wenn beide Zweifach-Interaktionen unterschiedlich stark sind.
Insgesamt können wir in einem dreifaktoriellen Design also bereits sieben Hypothesenpaare testen, und die entsprechende Quadratsummenzerlegung sieht wie folgt aus: \[\begin{gather*} SS_{\text{tot}}=SS_A+SS_B+SS_C+SS_{AB}+SS_{AC}+SS_{BC} +SS_{ABC}+SS_w\text{.} \end{gather*}\] Praktisch können wir aber auch festhalten: Varianzanalysen mit drei Faktoren lassen sich noch einigermaßen interpretieren; bei mehr als drei Faktoren kann die Interpretation zu einer großen Herausforderung werden.
4.4 Berechnung mit R
Die Berechnung einer zweifaktoriellen Varianzanalyse kann natürlich
sehr leicht mit R umgesetzt werden. Wir fokussieren hier auf die Nutzung
des Pakets ez und dessen Funktionen (mehr Informationen zur
Nutzung der Funktion aov() und des Paketes
afex finden sich hier). Die Daten wurden zu Beginn bereits
eingelesen und aufbereitet.
Die Funktion ezANOVA() erlaubt eine einfache
Spezifierung einer zweifaktoriellen Varianzanalyse. Die einzige Änderung
ist, dass nun zwei Faktoren dem Parameter between übergeben
werden:
daten_training$VP <- as.factor(daten_training$VP) # Faktorisieren
daten_training$Trainingstage <- as.factor(daten_training$Trainingstage) # Faktorisieren
daten_training$Altersgruppe <- as.factor(daten_training$Altersgruppe) # Faktorisieren
ez_ergebnis <- ezANOVA(data = daten_training, # welche Daten?
wid = .(VP), # Versuchsperson
dv = .(Sprachtest), # dependent variable
between = .(Trainingstage, Altersgruppe), # between-Faktor(en)
detailed = TRUE) # ...für schoRsch## Coefficient covariances computed by hccm()ez_ergebnis # Ergebnis ausgeben## $ANOVA
## Effect DFn DFd SSn SSd F p
## 1 Trainingstage 2 24 125.06667 74.8 20.064171 7.549788e-06
## 2 Altersgruppe 1 24 40.83333 74.8 13.101604 1.369112e-03
## 3 Trainingstage:Altersgruppe 2 24 24.26667 74.8 3.893048 3.433117e-02
## p<.05 ges
## 1 * 0.6257505
## 2 * 0.3531277
## 3 * 0.2449529
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 5 24 4.666667 34.8 0.6436782 0.6687841Eine schönere Formatierung erhalten Sie auch hier wieder mit der
Funktion anova_out() aus dem Paket
schoRsch:
library(schoRsch)
anova_out(ez_ergebnis) # formatierte Ausgabe## $`--- ANOVA RESULTS ------------------------------------`
## Effect MSE df1 df2 F p petasq getasq
## 1 Trainingstage 3.116667 2 24 20.06 0.000 0.63 0.63
## 2 Altersgruppe 3.116667 1 24 13.10 0.001 0.35 0.35
## 3 Trainingstage:Altersgruppe 3.116667 2 24 3.89 0.034 0.24 0.24
##
## $`--- SPHERICITY TESTS ------------------------------------`
## [1] "N/A"
##
## $`--- FORMATTED RESULTS ------------------------------------`
## Effect Text
## 1 Trainingstage F(2, 24) = 20.06, p < .001, np2 = .63
## 2 Altersgruppe F(1, 24) = 13.10, p = .001, np2 = .35
## 3 Trainingstage:Altersgruppe F(2, 24) = 3.89, p = .034, np2 = .24
##
## $`NOTE:`
## [1] "Reporting unadjusted p-values."Die Funktion ezStats() ist genauso wieder verwendbar, um
deskriptive Statistiken zu berechnen:
stats <- ezStats(data = daten_training, # welche Daten?
wid = .(VP), # Versuchsperson
dv = .(Sprachtest), # dependent variable
between = .(Trainingstage, Altersgruppe)) # between-Faktor(en)## Coefficient covariances computed by hccm()statsDie Funktion ezPlot() erlaubt auch wieder eine schnelle
grafische Darstellung:
ezPlot(data = daten_training, # welche Daten?
wid = .(VP), # Versuchsperson
dv = .(Sprachtest), # dependent variable
between = .(Altersgruppe, Trainingstage), # between-Faktor(en)
x = .(Trainingstage), # Faktor für die x-Achse?
split = .(Altersgruppe)) # Faktor für separate Linien?## Coefficient covariances computed by hccm()
5 Varianzanalyse mit Messwiederholung und weitere Designs
Bei der Einführung der einfaktoriellen Varianzanalyse in Teil 2 haben wir unterschiedliche Stichproben jeweils mit einer der drei Stufen des Faktors “Trainingstage” behandelt. Die einfaktorielle Varianzanalyse war letztlich die Verallgemeinerung des \(t\)-Tests für unabhängige Stichproben auf mehr als zwei Stichproben.
Wir haben in Teil 12 in Statistik I auch bereits den \(t\)-Test für zwei abhängige Stichproben kennengelernt, der dann Anwendung findet, wenn jede Versuchsperson Werte unter zwei Bedingungen liefert. Nun könnten wir die Untersuchung zur Wirkung der Anzahl Trainingstage auf die Leistung im Sprachtest auch so durchgeführt haben, dass jedes Kind nach dem ersten Trainingstag den Sprachtest bekommt, dann wieder einen Trainingstag und im Anschluss eine Parallelversion des Sprachtests und schließlich nach einem dritten Trainingstag abermals eine Parallelversion des Sprachtests. Jedes Kind hätte dann drei Datenpunkte geliefert, die allerdings voneinander abhängig sind: Ein Kind mit einem generell eher gutem sprachlichen Können wird in allen drei Erhebungen bessere Werte erhalten als ein Kind mit eher schlechten sprachlichen Leistungen. Zur Analyse solcher Daten benötigen wir die Varianzanalyse mit Messwiederholung (engl.: repeated measures oder within-subject Analysis of Variance) – also gewissermaßen die Verallgemeinerung des \(t\)-Tests für zwei abhängige Stichproben auf mehr als zwei abhängige Stichproben. Wie auch beim entsprechenden \(t\)-Test liegt hier der Fokus auf den bedingungsabhängigen Veränderungen innerhalb jeder Versuchsperson; weniger interessant ist das generelle Niveau der Versuchspersonen, also ob sie weniger oder mehr Punkte im Sprachtest erhalten.
Um das Prinzip der Varianzanalyse mit Messwiederholung zu verstehen, betrachten wir zunächst eine vereinfachte Methode ihrer Berechnung. Im Anschluss wenden wir uns der allgemeinen Vorgehensweise bei einfaktoriellen Varianzanalysen mit Messwiederholung zu, die große Ähnlichkeit zur zweifaktoriellen Varianzanalyse aus Teil 4 aufweist.
5.1 Ein einfacher Zugang: ipsative Werte
Wir betrachten nun einmal nur die Daten der älteren Kinder aus dem Eingangsbeispiel (siehe Teil 4) und stellen uns dabei vor, die Werte einer Zeile würden jeweils von einem Kind Versuchsperson stammen:
An den Mittelwerten der einzelnen Kinder können wir sehen, dass das allgemeine Leistungsniveau sich unterscheidet, d.h. dass große interindividuelle Unterschiede vorhanden sind. Trotzdem scheint es einen generellen Trend zu geben, von dem nur die Werte von Kind Nr. 4 abweichen. Dies wird auch im linken Teil der folgenden Abbildung klar, in welcher die grauen Linien die Daten einzelner Kinder darstellen und die dicke schwarze Linie die Mittelwerte der drei Bedingungen:
Da aber die interindividuellen Unterschiede gar nicht wirklich interessant sind, sondern nur das generelle Muster der Verläufe, können wir die Unterschiede auch eliminieren. Dazu subtrahieren wir von jedem Messwert einer Versuchsperson den Mittelwert dieser Versuchsperson:
Diese neuen Werte werden ipsative Werte genannt. Ihre zentrale Eigenschaft ist, dass der Mittelwert eines jeden Kindes nun 0 ist. Mit anderen Worten: Die Kinder unterscheiden sich nun nicht mehr in ihrem allgemeinen Leistungsniveau, denn die bestehenden interindividuellen Unterschiede wurden eliminiert. Das Muster der Unterschiede zwischen den drei Stufen innerhalb jeden Kinders und über alle Kinder hinweg bleibt jedoch erhalten. Dies wird auch aus dem rechten Teil der letzten Abbildung deutlich.
Zur weiteren Illustration wenden wir nun eine normale einfaktorielle
Varianzanalyse auf die ursprünglichen und auf die ipsativen Werte an
(wir benutzen dafür die Funktion oneway.test() die wir nur
kurz eingeführt haben):
daten_training$Trainingstage = as.factor(daten_training$Trainingstage)
# mit den Rohwerten
aov_ergebnis <- aov(Sprachtest ~ Trainingstage,
data = daten_training)
summary(aov_ergebnis)## Df Sum Sq Mean Sq F value Pr(>F)
## Trainingstage 2 19.6 9.800 2.625 0.113
## Residuals 12 44.8 3.733# mit den ipsativen Werten
aov_ipsativ <- aov(ipsativ ~ Trainingstage,
data = daten_training)
summary(aov_ipsativ)## Df Sum Sq Mean Sq F value Pr(>F)
## Trainingstage 2 19.6 9.800 11.31 0.00174 **
## Residuals 12 10.4 0.867
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Für die ursprünglichen Werte ist die leichte Leistungszunahme mit zunehmender Anzahl an Trainingstagen nicht signifikant, \(F(2,12)=2.63\), \(p=.113\), für die ipsativen Werte ergibt sich jedoch ein signifikantes Ergebnis, \(F(2,12)=11.31\), \(p=.002\). Was ist der Grund für diese unterschiedlichen Ergebnisse? Zunächst einmal ist festzuhalten, dass in beiden Fällen \(SS_A=19.60\) ist (Versuchen Sie einmal selber, dies nachzuvollziehen, indem Sie entweder die Quadratsummen von Hand berechnen oder mit R und einer entsprechenden Funktion, die die Quadratsummen mitausgibt.) Die Variabilität, die auf den Faktor Trainingstage zurückgeht, ist also identisch.
Während für die Ursprungsdaten aber \(SS_w=44.80\) ist, beträgt sie im Fall der ipsativen Werte nur noch \(SS_w=10.40\). Da der daraus resultierende Fehlerterm \(MS_w\) im Nenner des \(F\)-Bruchs steht, wird der \(F\)-Bruch größer. Die Eliminierung der interindividuellen Unterschiede hat demnach die Power des Tests erhöht, indem der Fehlerterm reduziert wurde. Tatsächlich sind interindividuelle Unterschiede eine große Quelle für den Fehlerterm der normalen Varianzanalyse und im Rahmen der Varianzanalyse mit Messwiederholung wird im Prinzip dafür gesorgt, dass diese interindividuellen Unterschiede keine Rolle spielen.
Ein Anmerkung dazu: Wir haben hier so getan, als könnten wir mit den ipsativen Werten eine normale einfaktorielle Varianzanalyse rechnen. Naheliegend ist es dann, zur Entscheidung die \(F\)-Verteilung mit \((J-1)\) und \((N-J)\) Freiheitsgraden zu verwenden. Allerdings sind die (ipsativen) Werte nicht unabhängig voneinander; dies ist aber eine Annahme der (normalen) einfaktoriellen Varianzanalyse. Streng genommen müssen daher die Nennerfreiheitsgrade reduziert werden auf \((J-1)(n-1)\), wobei \(n\) die Anzahl der Versuchspersonen ist. Hierdurch vergrößert sich zwar der kritische \(F\)-Wert; die Zunahme des empirischen \(F\)-Werts durch die Eliminierung der interindividuellen Unterschiede wiegt diesen Nachteil aber in den meisten Fällen auf.
5.2 Der Umgang mit interindividuellen Unterschieden
Zur Illustration des Prinzips einer Varianzanalyse mit Messwiederholung haben wir im vorangegangenen Abschnitt interindividuelle Unterschiede durch die Bildung ipsativer Werte “eliminiert”. Die Varianzanalyse mit Messwiederholung setzt ein ähnliches, aber flexibleres Verfahren ein: Hier werden die interindividuellen Unterschiede als separater Effekt aufgefasst, der wie die eigentlich interessierenden Effekte quantifiziert wird. So wird die dazugehörige Varianz aus den Daten genommen, auch wenn der Effekt der interindividuellen Unterschiede in den seltensten Fällen interessant ist und letztlich ignoriert wird.
Die einfaktorielle Varianzanalyse mit Messwiederholung wird also wie eine normale zweifaktorielle Varianzanalyse (siehe Teil 4) aufgefasst: 1. Der erste Faktor \(A\) ist dabei der interessierende Effekt mit \(J\) Stufen. 2. Den zweiten Faktor nennen wir Subjektfaktor \(S\) und er hat so viele Stufen, wie es Versuchspersonen gibt (also \(n\)). Im Prinzip hat das Design dann also \(J\cdot n\) Zellen und es folgt daraus, dass jede Zelle mit nur exakt einem Wert besetzt ist. Das heißt auch: Der Mittelwert pro Zelle entspricht dem einem Wert der Zelle und die Varianz jeder Zelle ist Null.
Nun fahren wir genauso fort, wie bei einer zweifaktoriellen Varianzanalyse und zerlegen die Gesamtvariabilität \(SS_{\text{tot}}\) in die Quadratsumme des Effektes A, also \(SS_A\), in die Quadratsumme des Subjektfaktors \(S\), also \(SS_S\), in die Quadratsumme der Interaktion \(AS\), also \(SS_{AS}\) sowie in die Quadratsumme innerhalb der Zellen, also \(SS_w\): \[ SS_{\text{tot}}=SS_A+SS_S+SS_{AS}+SS_w \]
Zur Illustration der Gemeinsamkeiten sind in der folgenden Abbildung die Quadratsummen- bzw. Varianzzerlegung der zweifaktoriellen Varianzanalyse und der einfaktoriellen Varianzanalyse mit Messwiederholung noch einmal gegenübergestellt:
Auch die zur Berechnung benötigten Formeln sind diejenigen, die wir in Teil 4 bereits eingeführt haben: Wir bestimmen \(SS_A\), \(SS_S\) sowie \(SS_{AS}\) und daraus die entsprechenden Mittleren Quadratsummen durch Division durch die jeweiligen Freiheitsgrade. Ignorieren wir nun die interindividuellen Unterschiede \(SS_S\) und fassen die außer \(SS_A\) verbleibenden Quadratsummen als Fehler auf, ist \(SS_{\text{Fehler}}=SS_{AS}+SS_w\).
Eine Besonderheit ist nun noch, dass jede Zelle ja mit nur einem Wert besetzt und die Varianz innerhalb jeder Zelle daher Null ist. Als Konsequenz gilt folglich auch \(SS_w=0\). Es bleibt also \(SS_{\text{Fehler}} = SS_{AS}\). Mit anderen Worten: Die Interaktion von Faktor \(A\) und Subjektfaktor \(S\) bildet die Grundlage für den Fehlerterm im \(F\)-Bruch der einfaktoriellen Varianzanalyse mit Messwiederholung: \[\begin{equation*} F=\frac{MS_A}{MS_{AS}} \end{equation*}\] Man kann auch hier zeigen, dass unter der Annahme der \(H_0\) identischer Populationsmittelwerte in allen Bedingungen dieser \(F\)-Bruch \(F\)-verteilt ist mit \(J-1\) Zählerfreiheitsgraden und \((J-1)(n-1)\) Nennerfreiheitsgraden. Die Entscheidungsregel und die Interpretation der Entscheidung entspricht dann wieder allen bisher behandelten Formen der Varianzanalyse.
5.3 Weiterführende Anmerkungen
5.3.1 Voraussetzungen und ihre Verletzungen
Für die in den Teilen 2-4 behandelten Varianzanalysen hatten wir die Unabhängigkeit der Stichproben vorausgesetzt. Diese Unabhängigkeit ist im Fall abhängiger Stichproben natürlich verletzt: Eine schnellere Versuchsperson wird in allen Bedingungen schneller sein als eine andere, langsamere Versuchsperson oder hat eine Versuchsperson hat bessere sprachliche Fähigkeiten als eine andere Versuchsperson. Die einzelnen Bedingungen sind eben nicht unabhängig voneinander. Um die Verteilung der entsprechenden Zufallsvariablen des \(F\)-Bruchs dennoch bestimmen zu können, wird stattdessen die folgende Voraussetzung gemacht:
- Die Varianzen aller paarweisen Differenzen (der Bedingungen) sind identisch.
Ist diese Voraussetzungen erfüllt, spricht man auch von Sphärizität (engl.: sphericity) Sphärizität ist immer gegeben, wenn sog. compound symmetry vorliegt, das heißt, wenn alle Varianzen und paarweisen Kovarianzen gleich sind (mehr Informationen dazu finden sich bei Lane, 2016). Daraus folgt, dass im Fall eines Faktors mit nur zwei Stufen Sphärizität immer gegeben ist. Ist die Sphärizität in einem Datensatz aber verletzt, verhält sich die Varianzanalyse mit Messwiederholung liberal, d.h., es kommt häufiger zu fehlerhaften Entscheidungen für die \(H_1\) – genau wie beim \(t\)-Test für unabhängige Stichproben oder der normalen Varianzanalyse, wenn Varianzhomogenität nicht gegeben ist und dann beim Welch-Test die Freiheitsgrade entsprechend angepasst werden. Ganz ähnliches wird dann auch im Fall der Verletzung von Sphärizität getan.
Ein gebräuchlicher Test der Sphärizitätsannahme ist Mauchly’s (1940) W-Test, der oft auch in den Ausgaben von Computerprogrammen auftaucht. Wichtig ist hierbei: Wird der Test signifikant, liegt keine Sphärizität vor. Dies ist ganz ähnlich wie beim Levene-Test auf Varianzhomogeniät. Im Fall einer Verletzung der Sphärizität wird üblicherweise eine Schätzung ihres Ausmaßes vorgenommen. Diese Schätzung wird mit \(\epsilon\) (ein kleines Epsilon) bezeichnet und kann Werte zwischen 0 und 1 annehmen. Schließlich werden dann die Freiheitsgrade des jeweiligen \(F\)-Bruchs mit \(\epsilon\) multipliziert; sie werden also kleiner, der kritische \(F\)-Wert steigt, und dies wirkt der Liberalität entgegen.
Es gibt hierbei verschiedene Methoden zur Schätzung von \(\epsilon\). Eine häufig empfohlene Schätzung ist die von Greenhouse und Geisser (1959), die auch von vielen Computerprogrammen ausgegeben wird. Eine alternative Variante ist die Huynh-Feldt-Korrektur (1976). Nonparametrische Alternativen zur Varianzanalyse mit Messwiederholung sind die Friedman-Rangvarianzanalyse oder der Kendall-W-Test.
5.3.2 Effektstärken
Für Varianzanalysen mit Messwiederholung werden i.d.R. die gleichen Effektstärken angegeben wie auch in den vorangegangenen Teilen. Einzig ist bei deren Schätzung zu beachten, dass der richtige Fehlerterm, also die Interaktion des Subjektfaktors mit dem interessierenden Faktor, verwendet wird.
5.3.3 Konfidenzintervalle
Das adäquate Vorgehen beim Bilden von Konfidenzintervallen im Fall abhängiger Stichproben ist in den letzten Jahren intensiv diskutiert worden. Ein sehr einflussreicher Vorschlag zur Berechnung von Konfidenzintervallen bei (mehr als zwei) abhängigen Stichproben geht auf Loftus und Masson (1994) zurück. Nach diesem Vorgehen wird das entsprechende Konfidenzintervall analog zu dem in der normalen Varianzanalyse (vgl. Teil 3) berechnet und basiert daher auf dem verwendeten Fehlerterm \(MS_{AS}\): \[\begin{gather*} \left[ M_j\pm t_{(J-1)(n-1);1-\frac{\alpha}{2}}\cdot\sqrt{\frac{MS_{AS}}{n}}\right]\;. \end{gather*}\] Ebenfalls gilt dann, dass sich zwei Mittelwerte~\(M_k\) und \(M_m\) signifikant voneinander unterscheiden, falls: \[\begin{gather*} |M_k-M_m| > \sqrt{2}\cdot t_{(J-1)(n-1);1-\frac{\alpha}{2}}\cdot\sqrt{\frac{MS_{AS}}{n}} \qquad k,m\in\{1,\ldots,J\}\;. \end{gather*}\]
5.4 Berechnung mit R
Wir beschreiben hier die Berechnung mit R direkt wieder mit den
Funktion aus dem Paket ez.
Wie zuvor verwenden wir dazu die Funktion ezANOVA(). Der
einzige Unterschied ist nun, dass wir keinen between-Faktor
haben, sondern stattdessen spezifizieren, dass es einen
within-Faktor gibt:
# faktorisieren
daten_training$VP <- as.factor(daten_training$VP)
daten_training$Trainingstage <- as.factor(daten_training$Trainingstage)
ez_ergebnis <- ezANOVA(daten_training,
wid = .(VP),
dv = .(Sprachtest),
within = .(Trainingstage), # Messwiederholungsfaktor
detailed = TRUE)
ez_ergebnis## $ANOVA
## Effect DFn DFd SSn SSd F p p<.05 ges
## 1 (Intercept) 1 4 4233.6 34.4 492.279070 2.442699e-05 * 0.9895288
## 2 Trainingstage 2 8 19.6 10.4 7.538462 1.444270e-02 * 0.3043478
##
## $`Mauchly's Test for Sphericity`
## Effect W p p<.05
## 2 Trainingstage 0.1873767 0.08110983
##
## $`Sphericity Corrections`
## Effect GGe p[GG] p[GG]<.05 HFe p[HF] p[HF]<.05
## 2 Trainingstage 0.5516866 0.04507414 * 0.6070624 0.03903435 *Im Gegensatz zu den bisherigen Ausgaben liefert
ezANOVA() nun allerdings drei Tabellen. Zunächst ist nur
die erste der drei Tabellen interessant, die die Ergebnisse der
Varianzanalyse im gewohnten Format angibt. Die angegebene Effektstärke
generalized \(\eta^2\)
(ges) ist für Varianzanalysen mit Messwiederholung nicht
identisch mit dem partiellen \(\eta^2\). Dieses müsste also noch per Hand
berechnet werden:
ez_ergebnis_kompakt <- ez_ergebnis$ANOVA[ez_ergebnis$ANOVA$Effect=="Trainingstage",]
ss_a <- ez_ergebnis_kompakt$SSn
ss_as <- ez_ergebnis_kompakt$SSd
eta2 <- ss_a / (ss_a + ss_as)
eta2## [1] 0.6533333Die zweite Tabelle (Mauchly's Test for Sphericity) gibt
an, ob die Sphärizitätsannahme verletzt ist. Dies ist in unserem
Beispiel nicht der Fall, \(p = .081\).
Wäre der Test signifikant geworden, hätten wir entsprechende
Korrekturverfahren aus der dritten Tabelle übernehmen müssen.
p[GG] steht hierbei für den korrigierten \(p\)-Wert nach Greenhouse- Geisser, der bei
verletzter Sphärizitätsannahme anstelle des üblichen \(p\)-Werts verwendet werden sollte).
GGe entspricht hierbei dem Korrekturfaktor \(\epsilon\) (hier also \(\epsilon=.552\)), mit dem die
ursprünglichen Freiheitsgrade multipliziert wurden um den korrigierten
\(p\)-Wert zu erhalten. In der Regel
werden jedoch unkorrigierte Freiheitsgrade, ggf. zusammen mit \(\epsilon\), berichtet, um eine leichtere
Interpretation zu ermöglichen.
Einfacher und übersichtlicher lässt sich die Ausgabe wiederum mit
anova_out() aus dem Paket schoRsch
formatieren, wobei direkt eine Schätzung des partiellen \(\eta^2\) berechnet wird. Außerdem liefert
anova_out() Informationen über Tests zur
Sphärizitätsannahme und korrigiert betroffene \(p\)-Werte automatisch
(anova_out() kann allerdings auch die Huynh-Feldt-Korrektur
oder korrigierte Freiheitsgrade ausgeben oder ein anderes \(\alpha\) für den Mauchly-Test verwenden;
siehe die Hilfefunktion ?anova_out):
anova_out(ez_ergebnis)## $`--- ANOVA RESULTS ------------------------------------`
## Effect MSE df1 df2 F p petasq getasq
## 1 (Intercept) 8.6 1 4 492.28 0.000 0.99 0.99
## 2 Trainingstage 1.3 2 8 7.54 0.014 0.65 0.30
##
## $`--- SPHERICITY TESTS ------------------------------------`
## Effect p_Mauchly GGEpsilon p_GG HFEpsilon p_HF
## 1 Trainingstage 0.081 0.552 0.045 0.607 0.039
##
## $`--- FORMATTED RESULTS ------------------------------------`
## Effect Text
## 1 (Intercept) F(1, 4) = 492.28, p < .001, np2 = .99
## 2 Trainingstage F(2, 8) = 7.54, p = .014, np2 = .65
##
## $`NOTE:`
## [1] "No adjustments necessary (all p_Mauchly > 0.05)."5.5 Weitere Designs
5.5.1 Mehrfaktorielle Varianzanalyse mit Messwiederholung
Wir haben nun die Varianzanalyse mit Messwiederholung eingeführt. Dabei hatten wir bisher nur einen (experimentellen) Faktor \(A\) betrachtet (Trainingstage) und den Fall dann wie eine zweifaktorielle Varianzanalyse behandelt, wobei der zweite Faktor der Subjektfaktor \(S\) war. Grundlage für den Fehlerterm war dann die Interaktion von \(A\) mit \(S\).
Selbstverständlich können wir auch mehrfaktorielle Varianzanalysen mit Messwiederholung durchführen. Dies kombiniert die Vorteile mehrfaktorieller Designs (z.B. Untersuchung von Interaktionen) und abhängiger Stichproben (z.B. höhere Power). Betrachten wir dazu das Beispiel einer zweifaktoriellen Varianzanalyse mit Messwiederholung auf den Faktoren \(A\) und \(B\). Ähnlich wie in der normalen Varianzanalyse können wir dann drei Hypothesen testen: die Haupteffekte \(A\) und \(B\) sowie deren Interaktion \(AB\).
Die zweifaktorielle Varianzanalyse mit Messwiederholung fassen wir nun als normale dreifaktorielle Varianzanalyse auf. Die drei Faktoren sind also \(A\), \(B\) und der Subjektfaktor \(S\). Der einzige Unterschied zur normalen zweifaktoriellen Varianzanalyse besteht in der Art der verwendeten Fehlerterme: Während in der normalen Varianzanalyse immer \(MS_w\) im Nenner der \(F\)-Brüche stand, wird im Fall der Messwiederholung immer die Interaktion des Subjektfaktors mit dem relevanten Faktor benutzt. Wir müssen also für jeden der drei \(F\)-Brüche einen anderen Fehlerterm verwenden und entsprechend zur Entscheidung auch jeweils andere \(F\)-Verteilungen heranziehen: \[\begin{gather*} F^A=\frac{MS_A}{MS_{AS}}\;,\qquad F^B=\frac{MS_B}{MS_{BS}}\qquad\text{und}\qquad F^{AB}=\frac{MS_{AB}}{MS_{ABS}} \end{gather*}\]
5.5.2 Gemischte Varianzanalyse
Betrachten wir noch einmal das Beispiel zur Messwiederholungs-Varianzanalyse und stellen uns vor, wir hätten doch auch die jüngeren Kinder mit einbezogen und alle Kinder den drei Trainingstagen mit jeweils anschließendem Sprachtest ausgesetzt. Während Trainigstage weiterhin ein Faktor mit Messwiederholung wäre, ist die Altersgruppe demgegenüber ein Faktor ohne Messwiederholung. Für den Fall (mindestens) eines Faktors, der mit abhängigen Stichproben realisiert wird, und (mindestens) eines weiteren Faktors, der mit unabhängigen Stichproben realisiert wird, werden sog. gemischte Varianzanalysen benötigt. Ein weiterer typischer Nutzungsfall sind Studien zur Wirksamkeit von Therapien: Hier gibt es dann eine behandelte Gruppe und eine Kontrollgruppe (d.h. einen Faktor mit zwei Stichproben). Beide Gruppen (bzw. alle Personen beider Gruppen) werden dann einmal vor und einmal nach der Therapie untersucht (d.h. dies ist ein Faktor mit Messwiederholung). Interessant ist dann vor allem oft die Interaktion beider Faktoren, da sie darüber Auskunft gibt, ob die Therapie in der behandelten Gruppe etwas (anderes) bewirkt hat, als in der Kontrollgruppe ohnehin passiert ist.
Auch bei gemischten Varianzanalysen ist das generelle Vorgehen identisch zu den bisher beschriebenen Varianzanalysen; der wesentliche Unterschied liegt in der Bestimmung der jeweiligen Fehlerterme für die einzelnen Effekte. Das Grundprinzip hierbei ist, dass Effekte des Gruppenfaktors anhand der Variabilität innerhalb der Gruppen getestet werden (ganz ähnlich wie bei der normalen Varianzanalyse), während Effekte von Faktoren mit Messwiederholung an entsprechenden Interaktion mit dem Subjektfaktor getestet werden. Bei letzteren Interaktionen ist allerdings zu bedenken, dass nicht alle Versuchspersonen in allen Stufen, sondern nur in jeweils einer Stufe des Gruppenfaktors, anzutreffen sind (man sagt auch, dass der Subjektfaktor “genested” ist). Als Konsequenz ist die Interaktion als Grundlage für den Fehlerterm nicht mehr mit den bekannten Formeln zu berechnen. Allerdings: Berechnet man die Quadratsumme der Interaktion des Messwiederholungsfaktors innerhalb jeder Gruppe einzeln und summiert dann auf, so erhält man die relevante Quadratsumme. Diese wird manchmal auch als \(SS_{\text{BS/A}}\) bezeichnet, wenn \(A\) der Gruppenfaktor und \(B\) der Messwiederholungsfaktor ist, und sie hat \(J(K-1)(n-1)\) Freiheitsgrade (wobei \(J\) bzw. \(K\) die jeweilige Anzahl der Stufen von \(A\) bzw. \(B\) sind).
Eine gemischte Varianzanalyse lässt sich mit R sehr einfach mit der
Funktion ezANOVA() berechnen. Bisher haben wir entweder das
Argument between = .() oder das Argument
within = .() verwendet; bei einer gemischten Varianzanalyse
würden wir einfach beide Argumente verwenden.
Zur Demonstration laden wir einen Datensatz, in welchem die Werte von 10 Versuchspersonen auf einer Wohlbefindens-Skala vorhanden sind. Jeweils 5 Versuchspersonen sind der Therapie- bzw. Kontrollgruppe zugeordnet und alle Versuchspersonen wurden zu zwei Zeitpunkten untersucht (vor und nachdem in der Therapiegruppe eine bestimmte Therapie durchgeführt wurde):
therapie_daten <- read.table("./Daten/therapie_daten.DAT",
header = TRUE)
head(therapie_daten) # ersten sechs Zeilen anzeigentherapie_daten$vp <- as.factor(therapie_daten$vp) # faktorisieren
therapie_daten$gruppe <- as.factor(therapie_daten$gruppe)
therapie_daten$messung <- factor(therapie_daten$messung,
levels = c("vorher","nachher"), # sorgt dafür, dass die....
ordered = TRUE) # ...Stufen nicht...
# ...alphabetisch sind. Zunächst verschaffen wir uns einen visuellen Überblick über die
Daten, wobei hier schon beide Argumente between und
within zum Einsatz kommen:
suppressPackageStartupMessages(library(ez))
library(schoRsch)
ezPlot(data = therapie_daten,
wid = .(vp),
dv = .(wohlbefinden),
between = .(gruppe),
within = .(messung),
x = .(messung),
split = .(gruppe))
Scheinbar gibt es in beiden Gruppen einen Zuwachs im Wohlbefinden über
den betrachteten Zeitraum hinweg; der Zuwachs fällt aber in der
Therapiegruppe deutlich größer aus (beide Linien sind nicht parallel!)
und wir könnten eine Interaktion zwischen Gruppe und Messzeitpunkt
erwarten:
ez.ergebnis <- ezANOVA(data = therapie_daten,
wid = .(vp),
dv = .(wohlbefinden),
between = .(gruppe),
within = .(messung),
detailed = TRUE)
anova_out(ez.ergebnis)## $`--- ANOVA RESULTS ------------------------------------`
## Effect MSE df1 df2 F p petasq getasq
## 1 (Intercept) 3.85 1 8 1077.19 0.000 0.99 0.99
## 2 gruppe 3.85 1 8 1.30 0.287 0.14 0.14
## 3 messung 0.15 1 8 85.33 0.000 0.91 0.29
## 4 gruppe:messung 0.15 1 8 33.33 0.000 0.81 0.14
##
## $`--- SPHERICITY TESTS ------------------------------------`
## [1] "N/A"
##
## $`--- FORMATTED RESULTS ------------------------------------`
## Effect Text
## 1 (Intercept) F(1, 8) = 1077.19, p < .001, np2 = .99
## 2 gruppe F(1, 8) = 1.30, p = .287, np2 = .14
## 3 messung F(1, 8) = 85.33, p < .001, np2 = .91
## 4 gruppe:messung F(1, 8) = 33.33, p < .001, np2 = .81
##
## $`NOTE:`
## [1] "Reporting unadjusted p-values."Diese Vermutungen spiegeln sich auch im Ergebnis der Varianzanalyse
wider: Der Haupteffekt messung entspricht den generell
höheren Werten in der “nachher”-Messung; die Interaktion zeigt an, dass
der Zuwachs in der Therapiegruppe stärker ausfällt.
5.6 Weitere Ergänzungen
5.6.1 Feste Faktoren (fixed factors) vs. Zufallsfaktoren (random factors)
Im bisherigen Verlauf haben wir implizit bereits sowohl feste Faktoren als auch Zufallsfaktoren kennengelernt, allerdings haben wir nie einen systematischen Unterschied zwischen beiden Typen diskutiert. Vorab ist festzuhalten, dass beide Typen idealisierte Konzepte darstellen. Von welchem Typ ein Faktor ist, ist aber mitunter auch Interpretationssache.
5.6.1.1 Feste Faktoren
Bei einem festen Faktor sind die Stufen des Faktors systematisch und absichtlich so gewählt worden, wie sie denn realisiert wurden. Sie stellen keine zufällig ausgewählte Stichprobe aus einer größeren Population dar und würden – wenn das Experiment wiederholt wird – genau so wieder zum Einsatz kommen. Im Allgemeinen gilt, dass experimentelle Manipulationen feste Effekte darstellen. So waren unsere Stufen 1, 2 oder 3 Trainingstage Beispiele für solche feste Faktoren.
5.6.1.2 Zufallsfaktoren
Bei einem Zufallsfaktor sind die Stufen unsystematisch und sind zufällig aus einer größeren Population gezogen worden. Würde das Experiment wiederholt werden, würden die Stufen anders ausfallen. Das typischste Beispiel hierfür sind Versuchspersonen. Im Rahmen der Varianzanalyse mit Messwiederholung haben wir den Subjektfaktor eingeführt, der soviele Stufen hat, wie es Versuchspersonen gibt: Hierbei handelt es sich um einen Zufallsfaktor, da bei Wiederholung des Experimentes die Stufen anders ausfallen würden, da wir andere Versuchspersonen ziehen würden.
Ein anderes typisches Beispiel stammt aus der Sprachpsychologie. Hier werden manchmal Stimuluswörter bzw. -phrasen verwendet, die für eine semantische Kategorie stehen, z.B. “gestern”, “früher”, “vor einiger Zeit”, … für “Vergangenheit”. Insofern diese Wörter auch zufällig aus allen möglichen Wörtern die für Vergangenheit stehen gezogen wurden, werden diese auch als Zufallsfaktor aufgefasst. Um dem Rechnung zu tragen, wird in der Sprachpsychologie oft eine sog. \(F_1\)-Varianzanalyse mit “Versuchsperson” als Zufallsfaktor und eine sog. \(F_2\)-Varianzanalyse mit “Stimulus” als Zufallsfaktor gerechnet und aus beiden \(F\)-Werten dann ein gemeinsamer Index gebildet und auf Signifikanz getestet. Alternativ werden sprachpsychologische Daten oft auch mit sog. gemischten linearen Modellen (engl.: linear mixed-effect models) ausgewertet, die eine simultane Berücksichtigung mehrerer Zufallsfaktoren erlauben.
5.6.2 Quadratsummen vom Typ I-III
Bei der Einführung der zweifaktoriellen Varianzanalyse sind wir davon ausgegangen, dass alle Zellen einen gleich großen Stichprobenumfang haben. Ist dies der Fall redet man auch von einem balancierten oder orthogonalen Design. Ein Merkmal diesen Falles ist, dass die Quadratsumme eines Effektes unabhängig davon ist, ob ein zweiter Faktor mit aufgenommen wird oder nicht. Unterscheiden sich die Stichproben in ihrem Umfang, sind die Effekte nicht mehr orthogonal, sondern voneinander abhängig: die Quadratsumme eines Effektes ändert sich dann auch durch die Hinzunahme eines zweiten Faktors (“in Abhängigkeit des zweiten Faktors”).
Das Problem kann verdeutlicht werden, wenn man bedenkt, dass bei
unseren bisherigen Betrachtungen Haupteffekthypothesen über globale
Erwartungswerte aufgestellt wurden, die ihrerseits mit globalen
Mittelwerten geschätzt werden. Als Beispiel nehmen wir einen Datensatz
bei dem von 12 Männern und 12 Frauen das Jahresgehalt in Einheiten von
1000€ ermittelt wurde. Zusätzlich wurde erfasst, ob die Personen einen
Universitätsabschluss haben oder nicht und wir kodieren eine zusätzliche
Variabe exclude, da wir im Folgenden dann auch Personen aus
dem Datensatz ausschließen wollen:
gehalt <- c(24,32,25,24,27,24,27,23,15,17,20,16,25,29,27,19,18,21,20,21,22,19,26,28)
geschlecht <- rep(c("m","w"), each = 12)
abschluss <- rep(c("ja","nein"), each = 6, times = 2)
exclude <- c(0,0,0,0,0,1,1,0,0,0,1,1,0,1,0,0,0,0,1,0,1,0,1,0)
daten <- data.frame(geschlecht, abschluss, gehalt, exclude)
head(daten)Interessieren wir nun eigentlich nur für einen möglichen Unterschied zwischen Geschlechtern, können wir uns die Mittelwerte anschauen:
aggregate(gehalt ~ geschlecht,
data = daten,
FUN = mean)Nun nehmen wir den zweiten Faktor dazu und erinnern uns daran, dass die globalen Mittelwerte als Mittelwerte der Zellmittelwerte berechnet werden:
# alle 4 Mittelwerte...
( means <- aggregate(gehalt ~ abschluss + geschlecht,
data = daten,
FUN = mean) )# ...und daraus die globalen Mittelwerte berechnen:
mean(means$gehalt[1:2]) # Männer## [1] 22.83333mean(means$gehalt[3:4]) # Frauen## [1] 22.91667Die globalen Mittelwerte entsprechen also den Mittelwerten, die wir oben ohne Berücksichtigung des zweiten Faktors erhalten haben und entsprechend würde sich die Quadratsumme auch nicht ändern. Dies ist genau dann der Fall, wenn alle Zellen gleich stark besetzt sind:
table(daten$geschlecht, daten$abschluss)##
## ja nein
## m 6 6
## w 6 6Nun schließen wir einen Teil der Personen aus dem Datensatz aus und erhalten dadurch ungleiche Zellbesetzungen:
daten <- subset(daten,
exclude == 0)
table(daten$geschlecht, daten$abschluss)##
## ja nein
## m 5 3
## w 5 3Nun können wir wieder die Mittelwerte für Männer und Frauen berechnen, ohne den zweiten Faktor zu berücksichtigen:
aggregate(gehalt ~ geschlecht,
data = daten,
FUN = mean)Schließlich berechnen wir nun die globalen Mittelwerte auf die gleiche Art wie oben und erhalten:
# alle 4 Mittelwerte...
( means <- aggregate(gehalt ~ abschluss + geschlecht,
data = daten,
FUN = mean) )# ...und daraus die globalen Mittelwerte berechnen:
mean(means$gehalt[1:2]) # Männer## [1] 22.36667mean(means$gehalt[3:4]) # Frauen## [1] 22.33333In diesem Fall ist nun der Unterschied deutlich kleiner, als er es ohne den zweiten Faktor war und dadurch ändert sich natürlich auch die Quadratsumme des respektiven Effektes. Dies kann genau bei ungleichen Zellbesetzungen passieren.
Um diesem Problem zu begegnen wurden verschiedene Varianten zur Berechnung der Quadratsummen vorgeschlagen, die oft als Typ I-III bezeichnet werden (die “inhaltlichen” Bezeichnungen variieren durchaus zwischen Autor:innen…) und je ihre eigenen Vor- und Nachteile haben. Die genaue Definition und dahinterstehende Logik ist leider sehr formal und wir belassen es daher hier bei verbalen Beschreibungen:
Typ I Quadratsummen werden auch als sequentiell bezeichnet. Hier wird zunächst möglichst viel der Variabilität dem ersten Effekt zugeschlagen; der nächste Effekt bindet dann noch einen entsprechenden Teil der verbleibenden Variabilität usw. Ein Vorteil dieser Variante ist, dass die additive Beziehung der Quadratsummen erhalten bleibt. Allerdings wird dies eingekauft durch den Nachteil, dass die Ergebnisse abhängig sind davon, in welcher Reihenfolge die Effekte in ein Modell gegeben wurden.
Typ II Quadratsummen werden auch als nicht-sequentiell hierarchisch bezeichnet. Zur Berechnung der Quadratsumme eines Effektes wird die Anwesenheit aller anderen Effekte berücksichtigt, die den aktuellen Effekt beinhalten.
Typ III Quadratsummen werden auch als nicht-sequentiell unique bezeichnet. Zur Berechnung der Quadratsumme eines Effektes wird die Anwesenheit aller anderen Effekte berücksichtigt. Das heißt, die Quadratsumme beinhaltet nur diejenige Variabilität, die eindeutig auf den Effekt zurückgeht.
Während bei gleichen Stichprobengrößen alle Typen zu gleichen
Ergebnissen kommen und sich additiv verhalten (d.h. sie sich zur totalen
Quadratsumme aufsummieren), trifft die Additivität bei den Typen II und
III im Falle ungleicher Stichprobengrößen nicht mehr zu. Die Funktion
ezANOVA() verwendet standardmäßig Typ II Quadratsummen
(kann aber mit type = 1 bzw. type = 3 andere
Quadratsummen berechnen); andere Statistiksoftware wie SPSS verwendet
standardmäßig Typ III Quadratsummen.
5.7 Literatur
Greenhouse, S., & Geisser, S. (1959). On methods in the analysis of profile data. Psychometrika, 24, 95–112.
Huynh, H. & Feldt, L. S. (1976). Estimation of the Box correction for degrees of freedom from sample data in randomized block and split-plot designs. Journal of Educational Statistics, 1, 69-82.
Lane, D. M. (2016). The assumption of sphericity in repeated-measures designs: What it means and what to do when it is violated. The Quantitative Methods for Psychology, 12, 114–122.
Loftus, G. R., & Masson, M. E. J. (1994). Using confidence intervals in within-subject designs. Psychonomic Bulletin & Review, 1, 476–490.
Mauchly, J. W. (1940). Significance test for sphericity of a normal n-variate distribution. The Annals of Mathematical Statistics, 11, 204–209.