In einem vorherigen Post habe ich bereits die einfaktorielle Varianzanalyse in R erklärt. Der nächste logische Schritt ist die zweifaktorielle Varianzanalyse. Während wir durch die einfaktorielle Varianzanalyse berechnen konnten, ob Gruppenunterschiede zwischen Gruppen unwahrscheinlich hoch sind, können wir anhand der zweifaktoriellen Varianzanalyse berechnen, ob Gruppenunterschiede nicht nur bei einem Faktor, sondern bei zwei Faktoren überzufällig sind. Meistens wird die zweifaktorielle Varianzanalyse angewandt, um ein experimentelles Design mit zwei Gruppen zu testen. Nehmen wir zum Beispiel an, du möchtest heraus finden, ob Rauchen und das Ausmaß an sportlicher Tätigkeit einen Effekt auf die Zeit hat, einen Halbmarathon zu laufen. Dein experimentelles Design umfasst zwei Gruppen: Die erste Gruppe unterscheidet sich in Raucher und Nichtraucher. Die zweite Gruppe unterscheidet sich im Ausmaß der sportlichen Betätigung:

Raucher (1) Nichtraucher (2)
Kein Sport(0)
Wenig Sport (1)
Viel Sport (2)

Die erste Gruppe soll ein paar Monate vor dem Halbmarathon regelmäßig rauchen, die andere Gruppe darf gar nicht rauchen (die ethischen Fragen für ein solches Design lassen wir einmal außen vor). Diese beiden Gruppen unterteilst du zusätzlich folgendermaßen: Ein Drittel pro Gruppe darf gar kein Sport treiben, ein weiteres Drittel macht nur ein bisschen Sport und eine dritte Gruppe macht viel Sport.

Nach 4 Monaten Training und Rauchen laufen alle 6 Gruppen einen Halbmarathon. Du stoppst am Ende die Zeit, die sie für den Halbmarathon gebraucht haben. Versuchen wir als nächstes heraus zu finden, ob Rauchen und die sportliche Betätigung einen Effekt auf die Zeit hat, einen Halbmarathon zu laufen. Der Einfachheit halber werden wir zunächst die zweifaktorielle Varianzanalyse ohne eine Interaktion berechnen.

Zunächst laden wir uns in R die nötigen Pakete:

library(tidyverse)
library(jmv)

Hier ist der Datensatz:

data <- tibble(
  endurance = c(120, 130, 150, 110, 99, 140, 100, 90, 120,
              240, 200, 180, 150, 170, 220, 119, 130, 115),
  smoker  = factor(c(rep(0, 9), rep(1, 9))),
  sports  = factor(c(rep(0, 3), rep(1, 3), rep(2, 3)) %>% rep(., 2))
)
  • endurace: Die Zeit in Minuten, die die Probanden benötigen, um den Halbmarathon zu Ende zu laufen
  • smoker: 0 = Kein Raucher, 1 = Raucher
  • sports: 0 = kein Sport, 1 = ein wenig Sport, 2 = viel Sport

Deskriptive Analyse

Schauen wir uns zunächst die deskriptive Statistik an:

jmv::descriptives(data, vars = c("endurance", "smoker",
                                "sports"),
                  sd = TRUE, missing = TRUE,
                  freq = TRUE)
DESCRIPTIVES

 Descriptives                                            
 ─────────────────────────────────────────────────────── 
                         endurance    smoker    sports   
 ─────────────────────────────────────────────────────── 
   N                            18        18        18   
   Missing                       0         0         0   
   Mean                        144                       
   Median                      130                       
   Standard deviation         42.9                       
   Minimum                    90.0                       
   Maximum                     240                       
 ─────────────────────────────────────────────────────── 


 FREQUENCIES

 Frequencies of smoker                              
 ────────────────────────────────────────────────── 
   Levels    Counts    % of Total    Cumulative %   
 ────────────────────────────────────────────────── 
   0              9          50.0            50.0   
   1              9          50.0           100.0   
 ────────────────────────────────────────────────── 


 Frequencies of sports                              
 ────────────────────────────────────────────────── 
   Levels    Counts    % of Total    Cumulative %   
 ────────────────────────────────────────────────── 
   0              6          33.3            33.3   
   1              6          33.3            66.7   
   2              6          33.3           100.0   
 ────────────────────────────────────────────────── 

Es scheint so, als ob wir keine fehlenden Werte haben. Im Mittel brauchen die Probanden zudem 144 Minuten, um den Halbmarathon abzuschließen. Bevor wir allerdings die Effekte berechnen, sollten wir die Daten zunächst auf Ausreißer prüfen: Zudem sollten wir sicherstellen, dass die Voraussetzungen für eine zweifaktorielle Varianzanalyse gelten. Da dies ein Testdatensatz ist, verzichten wir allerdings auf diese Annahmen. Mehr Informationen findest du hier.

Lass uns dennoch kurz prüfen, ob es Ausreißer gibt:

ggplot(data, aes(smoker, endurance)) +
  geom_boxplot(aes(fill = smoker), show.legend = FALSE) +
  geom_point() +
  facet_wrap(~ sports)

In Anbetracht unseres minimalen Datensatzes gibt es keine Ausreißer.

Zweifaktorielle Varianzanalyse

Mit Hilfe des Jamovi-Pakets in R können wir relativ problemlos, die zweifaktorielle Varianzanalyse berechnen:

model <- jmv::anova(data = data, 
           dep = "endurance", 
           factors = c("smoker", "sports"),
           modelTerms = list(
             "smoker", 
             "sports"),
           effectSize = "partEta",
           emMeans = list(
             c("smoker", "sports")))
model$main
ANOVA                                                                            
 ──────────────────────────────────────────────────────────────────────────────── 
                Sum of Squares    df    Mean Square    F        p         η²p     
 ──────────────────────────────────────────────────────────────────────────────── 
   smoker                12013     1          12013    18.60    < .001    0.571   
   sports                10172     2           5086     7.88     0.005    0.529   
   Residuals              9042    14            646                               
 ────────────────────────────────────────────────────────────────────────────────

Die p-Werte zeigen uns für beide Faktoren einen signifikanten Wert (p). Zusätzlich haben wir einen sehr großen Effekt für beide Faktoren (η²p). Wir können uns ebenfalls die Marginal Means ausgeben lassen, um zu sehen, wie sich die Gruppenmittelwerte voneinander unterscheiden:

model$emm

Die Frage ist nun, wie diese Werte eigentlich berechnet werden? Es ist das eine, einen inferenzstatistisches Verfahren zu berechnen, es ist etwas anderes nachvollziehen zu können, was im Hintergrund eigentlich passiert.

Zweifaktorielle Varianzanalyse händisch

Die Varianzanalyse ist im Grunde nichts anderes als das Verhältnis der Varianz, die wir manipulieren (systematisch), durch die unsystematische Varianz:

Varianzanalyse=systematischeVarianzunsystematischeVarianzVarianzanalyse = \frac{systematische Varianz}{unsystematische Varianz}

Die Summe zwischen der systematischen Varianz und der unsystematischen Varianz ist die Gesamtvarianz (SS-Total) in den Daten. Die Gesamtvarianz ist die quadrierte Abweichung der einzelnen Punkten vom Mittelwert der abhängigen Variable:

data <- rowid_to_column(data, "id")

ggplot(data, aes(id, endurance)) +
  geom_point() +
  geom_hline(yintercept = mean(data$endurance)) +
  geom_rect(aes(xmin = id, 
                xmax = id + (abs(endurance - mean(data$endurance))),   
                ymin = endurance, 
                ymax = mean(data$endurance), 
                alpha = .1), 
            fill = "#9999ff",
            data = data) +
  xlab("Personen ID") +
  ylab("IQ") +
  # theme_bw() +
  coord_fixed() +
  guides(alpha = FALSE, fill = FALSE)

Diese Varianz lässt sich unterteilen in die systematische Varianz und die unsystematische Varianz. Die systematische Varianz ist Summe der Varianz zwischen den Gruppen des Faktors A und des Faktors B. Wir berechnen zum Beispiel die quadrierte Abweichung des Faktors A vom Mittelwert der abhängigen Variablen diesen Faktors für jede Ausprägung dieses Faktors. Die unsystematische Varianz ist die Varianz, die wir nicht erklären können. Nicht jeder Wert in jeder unserer sechs Gruppen ist gleich. Obwohl Menschen nicht Rauchen und gar kein Sport treiben variieren sie trotzdem in der Zeit, die sie für einen Halbmarathon benötigen.

Im nächsten Schritte versuchen wir den Output unserer Analyse Stück für Stück zu berechnen.

Sum of Squares

Zunächst müssen wir die Sum of Squares berechnen:

ANOVA                                                                            
 ───────────────────────────────
                Sum of Squares       
 ───────────────────────────────
   smoker                12013    
   sports                10172       
   Residuals              9042                  
 ───────────────────────────────

Dazu müssen wir die Einzelwerte jeder Gruppe von dem Mittelwert der Gruppe abziehen und quadrieren und diese Werte am Ende aufsummieren:

SSbetween=j=1pnj(xjˉxˉ)2SS_{between} = \sum_{j=1}^p n_j (\bar{x_j} - \bar{x})^2
mean_smoker <- data %>%
  group_by(smoker) %>%
  summarise(mean_smoker = mean(endurance))

mean_sports <- data %>%
  group_by(sports) %>%
  summarise(mean_sports = mean(endurance))
  
data <- data %>%
  left_join(mean_smoker, by = "smoker") %>%
  left_join(mean_sports, by = "sports")

ss_between_smoker   <- sum((data$mean_smoker - mean(data$mean_smoker))**2)
ss_between_sports <- sum((data$mean_sports - mean(data$mean_sports))**2)

glue("SS_Between_Smoker ist {round(ss_between_smoker, 2)}")
glue("SS_Between_Sports ist {round(ss_between_sports, 2)}")
SS_Between_Smoker ist 12012.5
SS_Between_Sports ist 10172.33

Ein kurzer Check mit dem Output der Jamovi Funktion anova zeigt uns, dass die Werte stimmen. Ein Wert fehlt uns allerdings noch: Residuals. Die Residuen sind unsere unsystematische Varianz. Diese können wir berechnen, indem wir Varianzen zwischen den Gruppen von der totalen Varianz der Daten (siehe die Grafik weiter oben) abziehen. Die Gesamtvarianz ergibt sich aus der Summe der Sum of Squares der anova-Funktion:

SStotal=12013+10172+9042=31227SS_{total} = 12013 + 10172 + 9042 = 31227

Prüfen wir das zunächst:

ss_total <- sum((data$endurance - mean(data$endurance))**2)

glue("SS_Total ist {round(ss_total, 2)}")
SS_Total ist 31226.5

In der Tat, die Gesamtvarianz ist identisch. Nun können wir die unsystematische Varianz (Residuals) berechnen:

Residuals=SSTotalSSBetweensmokerSSBetweensportsResiduals = SS_{Total} - SSBetween_{smoker} - SSBetween_{sports}
sse <- ss_total - ss_between_smoker - ss_between_sports

glue("SSE ist {round(sse, 2)}")
SSE ist 9041.67

df

 ANOVA                                                                            
 ───────────────────────────────────── 
                Sum of Squares    df      
 ───────────────────────────────────── 
   smoker                12013     1     
   sports                10172     2     
   Residuals              9042    14                
 ──────────────────────────────────────

Als nächstes berechnen wir die Freiheitsgrade. Der Freiheitsgrad eines jeden Faktors wird durch die Anzahl der Ausprägungen des Faktors minus 1 berechnet. Da wir zwei Ausprägungen des Faktors smoker haben, haben wir einen Freiheitsgrad von 1 für diesen Faktor. Wir haben drei Ausprägungen des Faktors sports und daher zwei Freiheitsgrade für diesen Faktor. Der Faktor für die Residuen berechnet sich folgendermaßen: Anzahl der Probanden minus den Freiheitsgraden der beiden Faktoren minus 1:

dfresiduals=18121=14df_{residuals} = 18 - 1 - 2 - 1 = 14

Mean Square

 ANOVA                                                                            
 ───────────────────────────────────────────────────
                Sum of Squares    df    Mean Square    
 ────────────────────────────────────────────────────
   smoker                12013     1          12013    
   sports                10172     2           5086   
   Residuals              9042    14            646                               
 ────────────────────────────────────────────────────

Die Sum of Squares allein würden uns nicht helfen, den Quotienten zwischen der systematischen und der unsystematischen Varianz zu berechnen, da wir diese nicht anhand der Anzahl der Ausprägungen standardisiert haben. Daher teilen wir die Sum of Squares durch die Freiheitsgrade um die Mean Squares zu erhalten:

MeanSquaresmoker=120131=12013MeanSquare_{smoker} = \frac{12013}{1} = 12013 MeanSquaresports=101722=5086MeanSquare_{sports} = \frac{10172}{2} = 5086 MeanSquareresiduals=904214=646MeanSquare_{residuals} = \frac{9042}{14} = 646

F und p

ANOVA                                                                            
 ─────────────────────────────────────────────────────────────────────
                Sum of Squares    df    Mean Square    F        p           
 ─────────────────────────────────────────────────────────────────────
   smoker                12013     1          12013    18.60    < .001      
   sports                10172     2           5086     7.88     0.005  
   Residuals              9042    14            646                               
 ─────────────────────────────────────────────────────────────────────

Der F-Wert ist berechnet sich nun genau aus der Formel, die wir vorhin erklärt haben.

Varianzanalyse=systematischeVarianzunsystematischeVarianzVarianzanalyse = \frac{systematische Varianz}{unsystematische Varianz}

Die systematische Varianz ist die mittlere Abweichung der Abweichungen der Einzelwerte der Gruppen vom Gesamtmittelwert und die unsystematische Varianz ist die Varianz innerhalb der Gruppen, welche wir nicht erklären können. Je größer der Quotient der beiden Varianzen ist, desto unwahrscheinlicher ist es, dass unsere Gruppen aus der gleichen Population kommen.

F_smoker <- (ss_between_smoker / 1) / (sse / (nrow(data) - 1 - 2 - 1))

glue("F - Smoker ist {round(F_smoker, 2)}")
glue("p =  {round(1 - pf(F_smoker, 1, 72), 4)}")
F - Smoker ist 18.6
p =  1e-04
F_sports <- (ss_between_sports / 2) / (sse / (nrow(data) - 1 - 2 - 1))

glue("F - Sports ist {round(F_sports, 2)}")
glue("p =  {round(1 - pf(F_sports, 2, 72), 4)}")
F - Sports ist 7.88
p =  8e-04

Voilà, wir erhalten die gleichen Ergebnisse wie unsere Jamovi-Funktion.