5 Analysis of Variance (ANOVA)

5.1 Theorie in Kürze

  • Spezialfall der linearen Regression, bei der alle erklärenden Variablen kategoriell (also Faktoren) sind.
  • Aus historischen Gründen parallel zur linearen Regression immer noch sehr verbreitet.
  • Wir behandeln primär den “balancierten” Fall: Alle Gruppen haben gleich viele Beobachtungen.
  • “Ein-Weg ANOVA” wenn es nur eine erklärende Variable gibt; “Zwei-Weg ANOVA”, wenn es zwei erklärende Variablen (evtl. mit Wechselwirkung) gibt
  • Modell der Ein-Weg ANOVA: \[Y_{ij} = \mu + \alpha_i +\varepsilon_{ij}, \; \varepsilon_{ij} \; \textrm{i.i.d.} \sim N(0, \sigma^2)\] mit Nebenbedingung \(\sum_{i=1}^g \alpha_i = 0\) um die Parameter eindeutig identifizieren zu können (Vergleich zum Thema “identifizierbar”: Wenn man eine Gerade (mit zwei Parametern) mit nur einem Punkt fitten will, gibt es unendlich viele Möglichkeiten einer perfekt passenden Geraden und man sagt, dass die Parameter nicht eindeutig identifiziert werden können. Wenn eine Nebenbedingung hinzugefügt wird (z.B. Steigung muss gleich null sein), kann eine eindeutige Gerade gefunden werden, nämlich die horizontale Gerade, die durch den Punkt geht)
  • Nullhypothese: \(H_0: \alpha_1 =\alpha_2 = \ldots=\alpha_g=0\)
  • Teststatistik \(T = \frac{MS_B}{MS_W}\), wobei \(MS_B = SS_B/(g-1)\) und \(MS_W = SS_W/(g\cdot (p-1))\) bei \(g\) Gruppen mit je \(p\) Beobachtungen; man nennt die einzelnen Terme folgendermassen:
    • \(SS_B\) “between sum of squares”,
    • \(SS_W\) “within sum of squares”,
    • \(MS_B\) “between mean squares”,
    • \(MS_W\) “within mean squares”.
  • Falls \(H_0\) stimmt: \(T \sim F_{g-1, \, g \cdot (p-1)}\), also eine \(F\)-Verteilung mit \(g-1\) und \(g \cdot (p-1)\) “Degrees of freedom” (abgekürzt mit “df” oder “dof”)
  • Zwei-Weg ANOVA konzeptionell gleich, aber kompliziertere Formeln.
  • Häufig interessiert man sich neben obiger “globaler” Nullhypothese noch für spezifische Fragestellungen. Z.B. den paarweisen Vergleich aller Faktorstufen (z.B. mit Tukey HSD) oder sonstige Linearkombinationen der Faktorstufen (Kontraste); da hier potenziell viele Hypothesen getestet werden, sollte man für multiples Testen korrigieren.

  • Residuenanalyse wie bei linearer Regression.

  • Wichtige R Funktionen: aov, TukeyHSD und glht aus dem Paket multcomp.

5.2 Bsp: Ein-Weg ANOVA

In diesem Beispiel analysieren wir eine Ein-Weg ANOVA auf einem simulierten Datensatz.

5.2.1 Simulation

Zunächst simulieren wir Daten. Diesen Schritt würden wir in einer realen Auswertung natürlich nicht machen. Aber in der Erklärung der Methode hilft die Simulation uns zu sehen, ob wir die wahren Effekte finden konnten. Wir stellen uns eine kontinuierliche Zielgrösse vor (Rumpfkraft in einem standardisierten Test: Wie viele Sekunden kann man nach einer gewissen Behandlung in der Position “Planke” bleiben?). Zudem haben wir drei mögliche Behandlungen: Placebo, Therapie 1 und Therapie 2. Wir modellieren das mit einem Faktor x mit den drei Levels (Stufen): “P” (Placebo), “T1” (Therapie 1) und “T2” (Therapie 2). Pro Therapieform haben wir 10 Probanden gemessen (d.h., wir haben ein balanciertes Design mit 10 Beobachtungen pro Gruppe). Nach der Placebo Behandlung können Probanden im Mittel 30 Sekunden in der Planke bleiben. Nach Therapie 1 sind es 60 Sekunden und nach Therapie 2 sind es 65 Sekunden. Wir nehmen an, dass die gemessenen Werte pro Gruppe normalverteilt sind mit Standardabweichung 5 um den oben vorgegebenen Effekt.

5.2.2 Auswertung

Zunächst stellt sich die Frage, ob es überhaupt irgendeinen Unterschied zwischen den drei Behandlungsmethoden gibt. Um diese Frage zu beantworten, machen wir eine Ein-Weg ANOVA:

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## x            2   6668    3334   140.2 5.49e-15 ***
## Residuals   27    642      24                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Die Nullhypothese “Alle drei Behandlungen haben gleichen Erwartungswert” kann deutlich verworfen werden (p-Wert ist nur \(5.49 \cdot 10^{-15}\)). Bleiben wir kurz bei der ANOVA Tabelle: Die degrees of freedom sind 2 (\(df_1=g-1=2\)) und 27 (\(df_2=g \cdot (p-1)=3 \cdot 9=27\)). Between Sum-of-Squares \(SS_B\) ist \(6668\). Within Sum-of-Squares \(SS_W\) ist \(642\). Daraus berechnen wir die mean-squares: \(MS_B = SS_B / df_1 = 3334\); \(MS_W = SS_W / df_2 = 24\). Der Quotient der mean squares ergibt den F value, also den Wert der Teststatistik im F-Test: \(F = 3334/24 = 140.2\). Mit der \(F\)-Verteilung mit \(2\) und \(27\) Freiheitsgraden (df’s) kann so der p-Wert berechnet werden (in R: ).

Wir wissen jetzt also, dass es Unterschiede zwischen den Behandlungen gibt. Wir wissen aber nicht wo diese Unterschiede genau sind, oder wie gross sie sind. Die einfachste Möglichkeit das herauszufinden besteht darin, alle möglichen Paare der Faktorstufen zu vergleichen und jedes mal z.B. einen Zwei-Stichproben t-Test zu machen. Sehr ähnlich aber technisch etwas eleganter geht es mit “Tukey Honest Significant Difference”. Hier wird die Differenz aller Behandlungspaare inkl. Vertrauensintervall berechnet. Dabei wird auch gleich für multiples Testen korrigiert (bei \(k\) Faktorstufen gibt es \(k(k-1)/2\) mögliche Paare; es gibt also recht schnell mal sehr viele Tests, sodass ein Korrigieren für multiples Testen wichtig wird; “95% family-wise confidence level” bedeutet, dass die Wa., dass alle angegebenen Vertrauensintervalle ihren wahren Parameterwert überdecken gleich 95% ist; ohne Korrektur wäre die Wahrscheinlichkeit pro Zeile 95%, d.h., wenn man viele Zeilen hat, sinkt die Wa., dass alle angegebenen VIs ihren wahren Parameterwert enthalten deutlich unter 95%).

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = y ~ x, data = df)
## 
## $x
##            diff       lwr       upr     p adj
## T1-P  30.669982 25.262704 36.077260 0.0000000
## T2-P  32.504077 27.096799 37.911355 0.0000000
## T2-T1  1.834096 -3.573182  7.241374 0.6812485

Wir sehen, dass sowohl \(T1\) als auch \(T2\) deutlich unterschiedlich von \(P\) sind. Allerdings ist der Unterschied zwischen \(T1\) und \(T2\) nicht signifikant. D.h., die Therapien funktionieren beide in etwa gleich gut und sind beide deutlich besser als Placebo. Z.B. sind die Probanden nach Therapie 1 im Mittel etwa \(30.7\) Sekunden besser als die Probanden der Placebo-Gruppe (95%-VI: \([25.3, 36.1]\), dabei wurde schon für multiples Testen korrigert). Der (für multiples Testen korrigierte) p-Wert (“p adj”) ist praktisch 0.

Wenn man alle paarweisen Vergleiche anschaut, macht man potenziell sehr viele Tests, von denen einige evtl. vom Kontext her gar nicht interessant sind. Um beim Korrigieren für multiples Testen möglichst wenig Macht zu verlieren, sollte man nur die Hypothesen testen, die einen auch wirklich interessieren. Mit der Funktion glht (für General Linear Hypothesis Test) aus dem Paket multtest lassen sich beliebige Linearkombinationen der Faktorstufen (sog. Kontraste, also z.B. paarweise Vergleiche aber noch viel mehr) untersuchen. Im Folgenden vergleichen wir \(T1\) mit \(T2\). Zudem vergleichen wir \(P\) mit dem Mittelwert aus \(T1\) und \(T2\) (sind die Therapien im Mittel anders als Placebo?):

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = y ~ x, data = df)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)    
## T2-T1 == 0    1.834      2.181   0.841    0.645    
## T-P == 0     31.587      1.889  16.724   <1e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Auch hier werden wieder für multiples Testen kontrollierte p-Werte ausgegeben. Wir sehen (analog zu dem Ergebnis von vorhin), dass zwischen \(T1\) und \(T2\) kein nennenswerter Unterschied besteht. Allerdings ist der Unterschied zwischen \(T1\) und \(T2\) auf der einen Seite und \(P\) auf der anderen Seite sehr gross. Vertrauensintervalle für solche Kontraste können wie gewohnt mit dem Befehl berechnet werden:

## 
##   Simultaneous Confidence Intervals
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = y ~ x, data = df)
## 
## Quantile = 2.3635
## 95% family-wise confidence level
##  
## 
## Linear Hypotheses:
##            Estimate lwr     upr    
## T2-T1 == 0  1.8341  -3.3204  6.9886
## T-P == 0   31.5870  27.1231 36.0509

Z.B. ist der mittlere Unterschied zwischen \(T2\) und \(T1\) etwa \(1.8\) mit einem (für multiples Testen korrigierten) 95%-VI von \([-3.3, 7.0]\).

Wir hätten die Daten auch mit einer Linearen Regression auswerten können:

## 
## Call:
## lm(formula = y ~ x, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.8762  -3.2061  -0.2551   2.3306   8.3919 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   30.373      1.542   19.70  < 2e-16 ***
## xT1           30.670      2.181   14.06 6.07e-14 ***
## xT2           32.504      2.181   14.90 1.50e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.877 on 27 degrees of freedom
## Multiple R-squared:  0.9122, Adjusted R-squared:  0.9057 
## F-statistic: 140.2 on 2 and 27 DF,  p-value: 5.485e-15

Der F-Test in der letzten Zeile ist identisch mit dem Test der Ein-Weg ANOVA. Ansonsten lassen sich im summary output alle paarweisen Vergleiche mit dem Referenzlevel (hier Placebo) ablesen, allerdings wird nicht für multiples Testen korrigiert. Wir können den summary output der Linearen Regression auch mit gezielten Kontrasten mit der Funktion glht berechnen (der Unterschied bei den p-Werten kommt daher, dass glht für multiples Testen korrigiert):

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = y ~ x, data = df)
## 
## Linear Hypotheses:
##           Estimate Std. Error t value Pr(>|t|)    
## T1-P == 0   30.670      2.181   14.06   <1e-10 ***
## T2-P == 0   32.504      2.181   14.90   <1e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

Übrigens lässt sich die Funktion glht genau gleich auch auf Objekte der Funktion lm anwenden. Somit können wir also auch bei der linearen Regression beliebige Levels vergleichen und Vertrauensintervalle für die Unterschiede berechnen. Z.B. der Vergleich zwischen \(T1\) und \(T2\) (der im summary output von lm ja nicht erscheint):

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: lm(formula = y ~ x, data = df)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)
## T1-T2 == 0   -1.834      2.181  -0.841    0.408
## (Adjusted p values reported -- single-step method)

Zum Vergleich: Die Berechnung mit dem Anova Objekt führt zum gleichen Ergebnis:

## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: User-defined Contrasts
## 
## 
## Fit: aov(formula = y ~ x, data = df)
## 
## Linear Hypotheses:
##            Estimate Std. Error t value Pr(>|t|)
## T1-T2 == 0   -1.834      2.181  -0.841    0.408
## (Adjusted p values reported -- single-step method)

5.3 Bsp: Zwei-Weg ANOVA

Die Auswertung einer Zwei-Weg ANOVA (also mit zwei Faktoren als erklärende Variablen) ist grundsätzlich sehr ähnlich wie bei der Ein-Weg ANOVA.

Allerdings gibt es einen sehr wichtigen Unterschied: Bei der Ein-Weg ANOVA spielt es keine Rolle, ob es pro Gruppe gleich viele Beobachtungen (“balanciertes Design”) oder unterschiedlich viele Beobachtungen (“unbalanciertes Design”) gibt. Das ist bei der Zwei-Weg ANOVA anders. Bei einem unbalancierten Design ändert sich der summary output, wenn man die Reihenfolge der erklärenden Variablen im aov-Aufruf ändert. Dazu ein kleines Beispiel: Wir haben fiktive Daten von einem Wettlauf. Die Zielgrösse ist die gerannte Zeit in Sekunden. Erklärende Variablen sind Geschlecht und ob jemand vorher Energy Drink A oder B zu sich genommen hat. Die Daten sind unbalanciert; es gibt wesentlich mehr Frauen und mehr Personen mit Energy Drink B. Zunächst zeigen wir die Auswertung unter Verwendung der Funktion summary (dabei werden die sog. “Typ-1 Sum-of-Squares” verwendet):

##              Df Sum Sq Mean Sq F value   Pr(>F)    
## gender        1 2024.0  2024.0 263.719  < 2e-16 ***
## drink         1  455.2   455.2  59.316 9.05e-11 ***
## gender:drink  1   29.1    29.1   3.791   0.0558 .  
## Residuals    66  506.5     7.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##              Df Sum Sq Mean Sq F value Pr(>F)    
## drink         1 1145.9  1145.9 149.299 <2e-16 ***
## gender        1 1333.4  1333.4 173.737 <2e-16 ***
## drink:gender  1   29.1    29.1   3.791 0.0558 .  
## Residuals    66  506.5     7.7                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Für eine unbalancierte 2-Weg ANOVA empfehlen wir anstelle der Funktion summary die Funktion drop1 zu verwenden (es werden dann die sog. “Typ-3 Sum-of-Squares” verwendet):

## Single term deletions
## 
## Model:
## y ~ gender * drink
##              Df Sum of Sq     RSS    AIC  F value    Pr(>F)    
## <none>                     506.54 146.54                       
## gender        1   1352.40 1858.95 235.55 176.2110 < 2.2e-16 ***
## drink         1    484.24  990.78 191.50  63.0935 3.347e-11 ***
## gender:drink  1     29.09  535.64 148.45   3.7908   0.05579 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Single term deletions
## 
## Model:
## y ~ drink * gender
##              Df Sum of Sq     RSS    AIC  F value    Pr(>F)    
## <none>                     506.54 146.54                       
## drink         1    484.24  990.78 191.50  63.0935 3.347e-11 ***
## gender        1   1352.40 1858.95 235.55 176.2110 < 2.2e-16 ***
## drink:gender  1     29.09  535.64 148.45   3.7908   0.05579 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Details dieser Problematik werden in der Wahlfach-Vorlesung “Applied Analysis of Variance and Experimental Design” (401-0625-01) behandelt.