Saubere Statistik mit rstatix

24.05.2023, von Christian Glahn

Beim Modernisieren meiner Arbeitsumgebung für meine Vorlesungen bin ich auf die rstatix-Bibliothek gestossen. Die Bibliothek ist zwar recht neu, stellt aber alle wichtigen Funktionen für die Basisstatistik mit tidy-Datensätzen bereit. Die Bibliothek rstatix schliesst eine Lücke für eine konsistente und durchgehende Methodik mit modernem R. Das ist deshalb wichtig, weil die tidyverse-Methodik die Datenaufbereitung und -Analyse mit R viel einfacher und verständlicher macht. In diesem Post stelle ich vor, warum mich diese Bibliothek begeistert.

Die Einführung in die praktische Statistik mit R vermittelt die wichtigsten Tests für die Beantwortung einfacher Forschungsfragen. Weil meine Einführung in die Datenverarbeitung mit R auf den tidyverse-Prinzipien aufbaut, versuche ich auch die Auswertung und Analyse nach den gleichen Prinzipien zu gestalten. Dazu verwende ich seit einiger Zeit die infer-Bibliothek aus dem tidymodels-Projekt. Leider stellt infer nur die beiden fundamentalen Tests t-Test und (\chi^2)-Test für verkettete Funktionsfolgen bereit. Für alle anderen Tests und Verfahren musste ich bisher auf eine Kombination aus stats-R und broom’s tidy()-Funktion zurückgreifen. Das ist zwar nicht besonders kompliziert, erfordert aber ein paar nicht-triviale Konzepte bei der Programmierung.

rstatix kommt ohne komplizierte Zusätze aus und vereinheitlicht die Datenanalyse stark. Die besondere Stärke von rstatix liegt darin, dass sie sich elegant in den übergeordneten empirischen Forschungsprozess einbettet. Aus sich der empirischen Forschung bilden Modelle die Grundlage für die Beantwortung von Forschungsfrage. Aus solchen Modellen lassen sich die vermuteten Beziehungen zwischen Variablen direkt erfassen. In R werden Modelle in der folgenden Form geschrieben:

abhaengigeVariable ~ unabhaengigeVariable

Der ~-Operator zeigt an, welche Einflüsse für eine Variable vermutet werden. Dabei ist der linke Operand immer der beeinflusste Wert und der rechte Operand der beeinflussende Wert. Mit diesem Operanden wird angezeigt, dass eine Beziehung vermutet wird, für die (noch) keine mathematische bzw. algorithmische Funktion bekannt ist. Diese vermuteten Beziehungen werden normalerweise vor der Durchführung einer Studie festgelegt und dokumentiert.

Weil schon einfache Studien oft recht komplizierte Beziehungsnetzwerke untersuchen verwende ich die DAG-Methode, um kausale Strukturen zu dokumentieren und zu visualisieren. Die R-Bibliothek ggdag stellt die notwendigen Funktionen für diesen Schritt zur Verfügung und basiert auf ggplot2 für die Visualisierung. Hierbei spielt die dagify() Funktion eine besondere Rolle, weil die Beziehungen zwischen Variablen direkt in der Modellschreibweise eingegeben werden. Hier ein ganz einfaches Beispiel.

library(tidyverse)
library(ggdag)

dagify(
   abhaengigeVariable ~ unabhaengigeVariable
) %>%
    ggdag()

Aus dem DAG und den festgelegten Skalenniveaus der Variablen ergeben sich sehr oft welche Tests für die Analyse eine Studie notwendig sind. Nach der Datenerhebung müssen die Daten oft so transformiert werden, damit die Variablennamen und Skalenniveaus in den Messungen und in der Dokumentation übereinstimmen. Nun kommt rstatix ins Spiel: Die Bibliothek verwendet für alle Tests immer die gleiche Schreibweise wie dagify(), um die zu untersuchenden Beziehungen zu beschreiben.

Nehmen wir an, dass die unabhaengigeVariable genau zwei Werte enthält, dann können über diese Variable zwei Gruppen gebildet werden. Nehmen wir ausserdem an, dass die abhaengigeVariable metrisch (bzw. kardinal)-skaliert ist. Dann kann die (ungerichtete) Beziehung zwischen den beiden Variablen mit einem t-Test untersucht werden. Die rstatix-Funktion t_test() erwartet die folgende Schreibweise:

library(rstatix)

messdaten %>% 
    t_test(abhaengigeVariable ~ unabhaengigeVariable)

Das ist im Vergleich zur Funktion t.test() aus der mitgelieferten stats Bibliothek noch nicht besonders spektakulär. Sehr viele Studien müssen aber mehrere Variablen mit dem gleichen Skalenniveau für die gleichen Gruppen untersuchen. rstatix macht diesen Schritt extrem einfach, indem sie gruppierte Stichproben korrekt behandelt. Dazu müssen wir die abhängigen Variablen erst mit pivot_longer() in die Langform gebracht und dann über den urspünglichen Vektornamen gruppiert werden:

messdaten %>% 
    pivotLonger( 
        starts_with("abhaengig_"), 
        values_to = "wert"
    ) %>%
    group_by(name) %>%
    t_test(wert ~ unabhaengigeVariable)

Das Ergebnis sind alle t-Tests für die mit starts_with("abhaengig_") ausgewählten Variablen über die durch die unabhaengigeVariable gebildeten Gruppen. Mit Base-R oder infer ist dieser Schritt komplizierter, weil immer ein Umweg über eingebette Listen und unnest() gegangen werden muss. Das führt zu vielen Fehlern, die nichts mit den statistischen Prinzipien zu tun haben. Mit rstatix entfallen solche Umwege, weshalb man sich ganz auf die Forschungsfragen und die Datenorganisation konzentrieren kann.

Beispiel

Veranschaulichen wir uns die Vorgehensweise an einem konkreten Beispiel:

mtcars %>% 
    as_tibble(rownames = "type") -> 
    messdaten

Dieser Schritt wandelt die mtcars-Matrix in ein tibble um. Der Parameter rownames sorgt dafür, dass die Zeilennamen in der Spalte type gespeichert werden.

Nun werten wir diese Daten für die Vektoren mpg, disp, drat, wt und qsec mit einem t-Test über die Variable am aus. Damit prüfen wir, ob die Getriebeart (Automatik: am == 1, Schaltung: am == 0) diese Variablen beeinflusst. Weil die Variablennamen keine gemeinsame Komponente haben, müssen sie für pivot_longer() einzeln angegeben werden.

messdaten %>% 
    pivot_longer(c(mpg, disp, drat, wt, qsec), values_to = "wert") %>% 
    group_by(name) %>% 
    t_test(wert ~ am)

Das Ergebnis zeigt uns für jede abhängige Variable die Kennwerte für den zugehörigen (ungerichteten) t-Test. Für den Hypothesentest ist der wichtigste Wert der p-Wert in der letzten Spalte. Die Spalte name enthält wie gewohnt die Gruppenwerte. In unserem Fall sind das die Namen der untersuchten Variablen. Die Spalten group1 und group2 entsprechen den beiden Werten in der Variable am.

Fazit

Mir gefällt an rstatix besonders die konsistente Modellierung für alle Tests und die Umsetzung der Tests (U-Test, Kruskal-Wallis-Test, Dunn-Test) für kategorische Daten. Dadurch wird die Basisstatistik mit R zum Spaziergang.