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