Die einfaktorielle Varianzanalyse ist eines der gängisten Verfahren in der Statistik. Gewöhnlich wird sie nach dem t-Test vermittelt, wenn es um den Vergleich mehrerer Gruppenmitglieder geht.

library(tidyverse)
library(jmv)
library(car)
library(afex)

Wir werden in diesem Tutorial den F-Wert für den Faktor einer einfaktorielle Varianzanalyse berechnen. Unserer Faktor heißt in unserem Beispiel Diet. Wir haben 78 Personen randomisiert eine Diät zugeordnet und möchten nun vergleichen, ob sich das Gewicht dieser Personen zu Beginn der Diät signifikant, das heißt überwahrscheinlich voeinander unterscheidet. Laden wir zunächst den Datensatz.

diet <- read_csv("diet.csv")
glimpse(diet)
Observations: 78
Variables: 7
$ Person       <int> 25undefined 26undefined 1undefined 2undefined 3undefined 4undefined 5undefined 6undefined 7undefined 8undefined 9undefined 10undefined 11undefined 12undefined 13undefined 14undefined 27undefined 28undefined 29undefined 30undefined 31undefined 32undefined 33undefined 34undefined 35undefined 36undefined 37undefined 38undefined 39undefined 40undefined ...
$ gender       <int> NAundefined NAundefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0undefined 0...
$ Age          <int> 41undefined 32undefined 22undefined 46undefined 55undefined 33undefined 50undefined 50undefined 37undefined 28undefined 28undefined 45undefined 60undefined 48undefined 41undefined 37undefined 44undefined 37undefined 41undefined 43undefined 20undefined 51undefined 31undefined 54undefined 50undefined 48undefined 16undefined 37undefined...
$ Height       <int> 171undefined 174undefined 159undefined 192undefined 170undefined 171undefined 170undefined 201undefined 174undefined 176undefined 165undefined 165undefined 173undefined 156undefined 163undefined 167undefined 174undefined 172undefined 165undefined 171undefined 169undefined 174undefined 1...
$ pre.weight   <int> 60undefined 103undefined 58undefined 60undefined 64undefined 64undefined 65undefined 66undefined 67undefined 69undefined 70undefined 70undefined 72undefined 72undefined 72undefined 82undefined 58undefined 58undefined 59undefined 61undefined 62undefined 63undefined 63undefined 63undefined 65undefined 66undefined 68undefined 68...
$ Diet         <int> 2undefined 2undefined 1undefined 1undefined 1undefined 1undefined 1undefined 1undefined 1undefined 1undefined 1undefined 1undefined 1undefined 1undefined 1undefined 1undefined 2undefined 2undefined 2undefined 2undefined 2undefined 2undefined 2undefined 2undefined 2undefined 2undefined 2undefined 2undefined 2undefined 2undefined 3undefined 3undefined 3undefined 3undefined 3undefined 3undefined 3undefined ...
$ weight6weeks <dbl> 60.0undefined 103.0undefined 54.2undefined 54.0undefined 63.3undefined 61.1undefined 62.2undefined 64.0undefined 65.0undefined 60.5undefined 68.1undefined 66.9undefined 70.5undefined 69.0undefined 68.4undefined 81.1undefined 60.1undefined 56.0undefined 57...

Ein Blick auf die Daten verrät uns, dass der Faktor Diet von R als int gespeichert wird. Da manche Funktionen zur Berechnung einer einfaktoriellen Varianzanalyse die falschen Varianzen und Freiheitsgrade angeben, ist es sinnvoller, diese Variable in einen Faktor umzuwandeln.

diet <- diet %>%
  mutate(
    Diet = factor(Diet)
  )

Schauen wir zunächst die Variable pre.weight in einem Boxplot an, um einen ersten Eindruck zu erhalten, ob es beachtliche Ausreißer gibt.

ggplot(diet, aes(Diet, pre.weight, fill = Diet)) +
  geom_boxplot() +
  guides(fill = FALSE)

Es gibt einen Ausreißer in Gruppe 2. Die Person wiegt ungewöhnlich viel, mehr als 100 Kilogram. Solche Ausreißer können die Analyse verändern, einfachheitshalber lassen wir diese Person in der Analyse.

Annahmen

Die einfaktorielle Varianzanalyse umfasst ein paar zentrale Annahmen, die wir vor jeder Analyse prüfen sollten, auch wenn der Test relativ robust ist.

Normalverteilung der abhängigen Variable

shapiro.test(diet$pre.weight)
	Shapiro-Wilk normality test

data:  diet$pre.weight
W = 0.9692undefined p-value = 0.0554

Der Shapiro-Wilk Test zeigt, dass die Variable pre.weight gerade so noch als normalverteilt gelten kann.

Varianzhomogenität

Die Varianzen der unterschiedlichen Gruppen sollten ähnlich sein. Am einfachsten kann man das mit dem Levene-Test testen:

leveneTest(pre.weight ~ Diet, data = diet)
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  2   1.231  0.298
      75  

Die Varianzen sind nicht unterschiedlich voneinander, da der Test keine Signifikanz angibt. Wir können daher die Annahme der Varianzhomogenität zwischen den Gruppen annehmen.

Einfaktorielle Varianzanalyse durch Pakete

Es gibt eine Reihe an Paketen, mit denen man eine einfaktorielle Varianzanalyse berechnen kann. Hier sollen nur einmal die drei gängigsten angezeigt werden.

AOV aus dem stats-Paket

aov(pre.weight ~ Diet, data = diet) %>%
  summary
            Df Sum Sq Mean Sq F value Pr(>F)
Diet         2     90    44.9    0.58   0.56
Residuals   75   5770    76.9

aov ist die Standardfunktion, welche in R mitgeliefert wird.

aov_ez aus dem afex-Paket

afex::aov_ez(id = "Person", dv = "pre.weight",
             between = "Diet", data = diet)
Anova Table (Type 3 tests)

Response: pre.weight
  Effect    df   MSE    F ges p.value
1   Diet 2undefined 75 76.93 0.58 .02     .56

Der Output der aov_ez Funktion sieht etwas anders aus, es sind aber die gleichen Daten.

anova aus dem jmv-Paket

jmv::anova(diet, dep = "pre.weight",
           factors = c("Diet"))
 ANOVA

 ANOVA                                                                  
 ────────────────────────────────────────────────────────────────────── 
                Sum of Squares    df    Mean Square    F        p       
 ────────────────────────────────────────────────────────────────────── 
   Diet                   89.9     2           44.9    0.584    0.560   
   Residuals            5769.6    75           76.9                     
 ────────────────────────────────────────────────────────────────────── 

Das Paket jamovi liefert einen SPSS ähnlichen Output und ist für gewöhnliche Analysen wie Varianzanalyse sehr praktisch.

Einfaktorielle Varianzanalyse händisch in R

Im nächsten Schritt versuchen wir die Werte diesen Outputs händisch zu berechnen, um besser zu verstehen, was diese drei Funktionen im Hintergrund machen. Der F-Wert (hier 0.584) berechnet sich aus dieser Formel:

F=MSbetweenMSwithinF = \frac{MS_{between}}{MS_{within}}

MS steht für mean squares und kennzeichnet die durchschnittliche Größe von Quadraten. Welche Quadrate sehen wir gleich. Im Grunde ist der F-Wert nichts anderes als das Verhältnis der durschnittlichen Varianz innerhalb (bzw. within) der Gruppen durch die durschnittliche Varianz zwischen (bzw. between) den Gruppen.

Sum of Squares

Beginnen wir den ersten Teil des Outputs zu berechnen:

 ANOVA

 ANOVA                                                                  
 ──────────────────────────────
                Sum of Squares   
 ──────────────────────────────
   Diet                   89.9   
   Residuals            5769.6                
 ─────────────────────────────

Zunächst müssen wir die Sum of Squares des Varianz zwischen den Gruppen (hier: Diet) und der Varianz innerhalb der Gruppen (hier: Residuals) berechnen.

Sum of Squares zwischen den Gruppen

SSbetween=j=1pnj(xjˉxˉ)2SS_{between} = \sum_{j=1}^p n_j (\bar{x_j} - \bar{x})^2

In anderen Worten: Wir berechnen die Abweichung der Mittelwerte der Gruppen vom Gesamtmittelwert unserer Variable pre.weight.

group_means <- diet %>%
  group_by(Diet) %>%
  summarise(group_mean = mean(pre.weight))
  
diet <- diet %>%
  left_join(group_means, by = "Diet")

Nun können wir die Varianz zwischen den Gruppen berechnen.

(ss_between <- sum((diet$group_mean - mean(diet$pre.weight))**2))
[1] 89.8608

In der Tat, dieser Wert ist identisch mit unserem Output der Funktion jmv::anova.

Sum of Squares innerhalb den Gruppen

SSwithin=j=1pi=1nj(xijxjˉ)2SS_{within} = \sum_{j=1}^p \sum_{i=1}^{n_j} (x_{ij} - \bar{x_j})^2

In anderen Worten: Wir berechnen die Abweichung der Einzelpunkte vom Mittelwert der jeweiligen Gruppe. Um dies auszurechnen, müssen wir zunächst den Gruppenmittelwert an unseren Datensatz anhängen:

(ss_within <- sum((diet$pre.weight - diet$group_mean)**2))
[1] 5769.59

Auch hier entspricht der Mittelwert unserem Output.

Mean Squares

 ANOVA

 ANOVA                                                                  
 ────────────────────────────────────────────────────
                Sum of Squares    df    Mean Square        
 ────────────────────────────────────────────────────
   Diet                   89.9     2           44.9    
   Residuals            5769.6    75           76.9                     
 ────────────────────────────────────────────────────

Als nächstes müssen wir die Mean Squares berechnen. Dies ist nicht mehr sonderlich kompliziert, da wir die Sum of Squares bereits berechnet haben. Wir müssen die Sum of Squares nun lediglich durch die Freiheitsgrade teilen. Der Freiheitsgrad von SS_between ist immer die Anzahl der Gruppen - 1 und der Freiheitsgrad der Residuals oder SS_within ist immer die Anzahl der Datenpunkte - der Anzahl der Gruppen.

(ms_between <- ss_between / (3 - 1))
[1] 44.9304
(ms_within <- ss_within / (nrow(diet) - 3))
[1] 76.9278

F-Wert

Zuletzt müssen wir lediglich den F-Wert berechnen:

F=MSbetweenMSwithinF = \frac{MS_{between}}{MS_{within}}
(F <- ms_between / ms_within)
[1] 0.584059

Die Signifikanz können wir aus der F-Verteilung berechnen:

df(F, df1 = 2, df2 = 75)
[1] 0.551556

Voilà, wir haben den gesamten Output der einfaktoriellen Varianzanalyse dieser drei Funktionen berechnet.

 ANOVA

 ANOVA                                                                  
 ────────────────────────────────────────────────────────────────────── 
                Sum of Squares    df    Mean Square    F        p       
 ────────────────────────────────────────────────────────────────────── 
   Diet                   89.9     2           44.9    0.584    0.560   
   Residuals            5769.6    75           76.9                     
 ──────────────────────────────────────────────────────────────────────