4 Theoretisch ja, praktisch auch! Tests in R
Üblicherweise geht es nach den deskriptiven Statistiken in den meisten Fällen ans Eingemachte - die statistischen Tests. Nehmen wir an, in zwei Gruppen beobachten wir unterschiedliche Mittelwerte. Eine mögliche Erklärung könnte sein, dass der bloße Zufall hier Unterschiede erscheinen lässt, wo in Wirklichkeit keine sind. Eine andere Möglichkeit wäre dann, dass in der Tat Unterschiede zwischen den Gruppen bestehen, und diese nicht alleine dem Zufall zu schulden sind.
Die gesamten theoretischen Grundlagen hier zu wiederholen würde sicher den Rahmen sprengen, insofern beschränken wir uns im Folgenden mit der Durchführung verschiedener statistischer Tests in R. Welcher Test wann der geeignete ist, ist manchmal nicht leicht zu erkennen. Für eine erste grobe Einschätzung kann aber die folgende Tabelle hilfreich sein.
Mögliche statistische Tests in Abhängigkeit der Variablen. (Tabelle ursprünglich aus Wikipedia unter Creative Commons Lizenz)
4.1 Lernziele
- Daten auf Normalverteilung prüfen (Shapiro-Test)
- Einfache statistische Tests rechnen (T-Test, Wilcoxon-Test)
- Effektstärke berechnen und eigene Funktionen schreiben
- ANOVA durchführen
4.2 Normalverteilung prüfen
Wie aus obiger Tabelle ersichtlich ist, kann es mitunter zur Auswahl des korrekten statistischen Tests wichtig sein zunächst zu prüfen, ob die Werte einer Variablen normalverteilt sind. Natürlich können wir versuchen visuell abzuschätzen, ob dies der Fall ist. Für ein fiktives Beispiel wären wir zudem interessiert, die Normalverteilung für einen Wert innerhalb verschiedener Gruppen zu untersuchen. Eine schnelle deskriptive Statistik bzw. Visualisierung könnte wie folgt aussehen.
Wir nutzen hier die skew()-Funktion aus dem psych Paket, um die Schiefe der Verteilung zu erhalten (bei einer Normalverteilung nahe 0, bei positiven Werten ist der Ausläufer nach rechts länger, bei negativen Werten der Ausläufer nach links).
library(psych)
daten %>%
group_by(untersuchung) %>%
summarise(schiefe_dlp = skew(dlp))## # A tibble: 6 × 2
## untersuchung schiefe_dlp
## <chr> <dbl>
## 1 Abdomen 2.10
## 2 Becken 3.39
## 3 Gesicht 2.57
## 4 Kopf 0.223
## 5 Thorax 3.58
## 6 Traumaspirale 0.0290daten %>%
ggplot(aes(x = dlp, fill = untersuchung)) +
geom_density() +
facet_wrap(vars(untersuchung)) +
labs(x = "DLP [mGy*cm]",
y = "",
fill = "Untersuchungsbezeichnung") +
theme_bw()
Wir sehen hier bereits, dass die DLP-Werte für Untersuchungen des Kopfes und bei Traumaspiralen vermutlich normalverteilt sind. Doch um wirklich sicher zu gehen, sollte ein entsprechender statistischer Test verwendet werden. Für diese Frage bietet sich der Shapiro-Wilk-Test an.
Der Shapiro-Wilk-Test nimmt als sog. Nullhypothese an, dass die Daten normalverteilt sind. Ein kleiner p-Wert zeigt in diesem Fall dann entsprechen an, dass die Alternativhypothese anzunehmen ist, also dass die Daten nicht normalverteilt sind. Erfreulicherweise findet sich für die meisten statistischen Tests in R eine entsprechend benannte Funktion. Um also bspw. die DLP-Werte im gesamten Datensatz auf Normalverteilung zu untersuchen, können wir einfach die Funktion shapiro.test() nutzen.
shapiro.test(daten$dlp)##
## Shapiro-Wilk normality test
##
## data: daten$dlp
## W = 0.8234, p-value < 2.2e-16In diesem Fall zeigt also der p-Wert < 2,2 x 10-16 an, dass die Alternativhypothese angenommen werden sollte - die DLP-Werte im Gesamtdatensatz also nicht normalverteilt sind.
4.2.1 Ein kurzer Ausflug in Ergebnisobjekte
In vielen Fällen sind die Ausgaben solcher statistischen Tests etwas unintuitiv formatiert und die einzelnen Werte weiterzuverwenden erscheint zunächst schwierig. Es lohnt daher an dieser Stelle ein kurzer Ausflug, denn wir können durchaus auf Einzelteile der Ausgabe direkt zugreifen. Um einen etwas tieferen Einblick zu erhalten, speichern wir zunächst das Ergebnis der shapiro.test() Funktion in einer Variablen und schauen uns die Struktur der Variablen an.
shapiro.ergebnis <- shapiro.test(daten$dlp)
str(shapiro.ergebnis)## List of 4
## $ statistic: Named num 0.823
## ..- attr(*, "names")= chr "W"
## $ p.value : num 8.48e-46
## $ method : chr "Shapiro-Wilk normality test"
## $ data.name: chr "daten$dlp"
## - attr(*, "class")= chr "htest"Wie man erkennen kann, handelt es sich um eine Liste von Werten, die einfach—wenn sie nicht in einer Variablen gespeichert sind—gemeinsam auf der Konsole ausgegeben werden. Das Dollarzeichen $ in der Ausgabe der str() Funktion zeigt hier aber bereits an, dass wir auch einzelne Teile des Ergebnisses direkt referenzieren können, ähnlich den Spalten eines Dataframes.
shapiro.ergebnis$p.value## [1] 8.48211e-46shapiro.ergebnis$p.value < 0.05## [1] TRUEWären wir also bspw. in obigem fiktiven Beispiel daran interessiert die Normalverteilung innerhalb der Subgruppen zu prüfen, könnten wir in Tidyverse-Syntax eine Verkettung von Befehlen schreiben, an deren Ende eine ansprechend formatierte Tabelle stünde, die in einer Spalte auch gleich prüft, ob der p-Wert unterhalb des Signifikanzniveau von 0.05 liegt.
daten %>%
group_by(untersuchung) %>%
summarise(shapiro_p.wert = shapiro.test(dlp)$p.value) %>%
mutate(signifikant_0.05 = shapiro_p.wert < 0.05)## # A tibble: 6 × 3
## untersuchung shapiro_p.wert signifikant_0.05
## <chr> <dbl> <lgl>
## 1 Abdomen 1.48e-30 TRUE
## 2 Becken 4.46e-10 TRUE
## 3 Gesicht 1.01e-12 TRUE
## 4 Kopf 3.06e-20 TRUE
## 5 Thorax 4.29e-41 TRUE
## 6 Traumaspirale 5.19e- 5 TRUEWir sehen also, dass—entgegen der ursprünglichen Annahme—auch in den Gruppen, die rein visuell und deskriptiv normalverteilt erschienen (Kopf und Traumaspirale), in Wirklichkeit keine Normalverteilung vorliegt.
4.3 T-Test und Wilcoxon-Test
Nachdem nun die Daten auf Normalverteilung geprüft wurden, können wir mit dem entsprechenden statistischen Test untersuchen, ob zwischen einzelnen Gruppen ein signifikanter Unterschied in der Verteilung der Werte für eine Variable existiert. In unserem Falle wäre, da die Daten nicht normalverteilt sind, der Wilcoxon-Mann-Whitney-Test die richtige Wahl. Ähnlich wie für den Shapiro-Wilk-Test shapiro.test() existiert in diesem Falle die Funktion wilcoxon.test(). Die Funktion erwartet als Parameter entweder zwei Vektoren mit Zahlenwerten (oder ein Dataframe und eine sogenannte Formel, aber dazu später). Für ein einfaches Beispiel wollen wir den Vergleich der DLP-Werte von Thorax-Untersuchungen mit denen von Abdomen-Untersuchungen untersuchen. Hierzu ziehen wir uns zunächst die entsprechenden Werte aus dem Dataframe heraus und speichern sie in separaten Variablen, anschließend können diese Vektoren der Funktion wilcoxon.test() übergeben werden. Im Übrigen geht R bei T-Test und Wilcoxon-Test davon aus, dass es sich um unverbundene Stichproben handelt. Sicherheitshalber sollte man sich aber angewöhnen, dies auch explizit mithilfe des Parameters paired = FALSE anzugeben.
# erst die Funktion pull() extrahiert aus dem Dataframe einen einfachen Vektor
dlp.abdomen <- daten %>% filter(untersuchung == "Abdomen") %>% pull(dlp)
dlp.thorax <- daten %>% filter(untersuchung == "Thorax") %>% pull(dlp)
wilcox.test(dlp.abdomen, dlp.thorax, paired = FALSE)##
## Wilcoxon rank sum test with continuity correction
##
## data: dlp.abdomen and dlp.thorax
## W = 445818, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0Wie erwartet, zeigt hier der p-Wert < 2,2 x 10-16 an, dass die DLP-Werte von Thorax-Untersuchungen und Abdomen-Untersuchungen signifikant verschieden sind.
Genau so einfach ließe sich entsprechend auch ein T-Test t.test() rechnen, der zwar in diesem Falle vielleicht nicht optimal wäre, weil wie bereits besprochen keine Normalverteilung der Werte vorliegt, aber das ignorieren wir einmal.
t.test(dlp.abdomen, dlp.thorax, paired = FALSE)##
## Welch Two Sample t-test
##
## data: dlp.abdomen and dlp.thorax
## t = 14.719, df = 1032.2, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 554.7071 725.3559
## sample estimates:
## mean of x mean of y
## 1037.5276 397.49614.3.1 Ein kurzer Ausflug in Formeln
Mit Formeln werden in R nicht übliche Formeln in der Art f(x) = mx + b bezeichnet, sondern bei der Definition von statistischen Tests oder Modellen ein Ausdruck, der beschreibt welche Zielvariable von welchen Einflussvariablen abhängen soll. Eine Formel in R besteht immer aus einer Tilde ~ und einem oder mehreren Ausdrücken, die jeweils rechts (RHS = right hand side) bzw. links (LHS = left hand side) der Tilde ~ stehen. Die Zielvariable steht dabei links (LHS) und die möglichen Einflussvariablen rechts (RHS).
Wollten wir also beispielsweise ausdrücken, dass wir das DLP in Abhängigkeit der Untersuchungsbeschreibung auf statistische Signifikanz untersuchen wollen würden, könnten wir das folgendermaßen ausdrücken: dlp ~ untersuchung. Wollten wir also in unserem Beispiel die beiden Gruppen Thorax und Abdomen vergleichen, könnten wir den Ausdruck daten %>% filter(untersuchung %in% c("Abdomen", "Thorax")) nutzen, um ein Dataframe zu erzeugen, das nur diese beiden Gruppen enthält. Bei dem Aufruf der Funktion t.test() müssten wir dieses Dataframe als Parameter data = übergeben und als ersten Parameter die oben genannte Formel dlp ~ untersuchung.
# auch hier hilft eine Schreibweise mit Zeilenumbrüchen nicht den Überblick zu verlieren
t.test(
dlp ~ untersuchung,
data = daten %>% filter(untersuchung %in% c("Abdomen", "Thorax")),
paired = FALSE
)##
## Welch Two Sample t-test
##
## data: dlp by untersuchung
## t = 14.719, df = 1032.2, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Abdomen and group Thorax is not equal to 0
## 95 percent confidence interval:
## 554.7071 725.3559
## sample estimates:
## mean in group Abdomen mean in group Thorax
## 1037.5276 397.4961Das mag zunächst etwas kompliziert wirken. Aber solche Schreibweisen haben den Vorteil, dass nicht für jeden Test neue Variablen definiert werden müssen, was bei langen und komplizierten Auswertungen schonmal dazu führen kann, dass man den Überblick verliert. Für einfache statistische Tests ist die Formelschreibweise nicht unbedingt nötig, und beide Herangehensweisen liefern identische Ergebnisse. Im späteren Kapitel zu Machine Learning werden wir allerdings um die Verwendung dieser Formelschreibweise nicht herum kommen.
4.3.2 Ein kurzer Ausflug in Funktionen
Im Live-Webinar hatten wir darüber gesprochen, dass bspw. beim T-Test neben dem p-Wert auch die sog. Effektstärke relevant sein kann. Zwar existiert für R beispielsweise das Paket effsize, das Funktionen enthält um die Effektstärke zu berechnen, aber gelegentlich kommt vielleicht man an einen Punkt, wo auf die Schnelle kein Paket zu finden ist, dass genau die Funktion enthält, die man sucht. Für solche Fälle bietet R die Möglichkeit, eigene Funktionen zu definieren, die der Vollständigkeit hier genannt werden soll, obwohl diese spezielle Funktionalität für den Anfang vielleicht etwas komplex ist.
Für ein kurzes Beispiel wollen wir einen Spezialfall betrachten, in dem für unverbundene Stichproben gleicher Gruppengröße die Effektstärke anhand von T-Statistik und df-Wert des T-Testes brechnet werden kann. Wir wollen dabei das Ergebnisobjekt der Funktion t.test() nehmen und direkt die Werte für T-Statistik und df extrahiert und sodann die Effektstärke berechnet werden kann.
Eine Funktion kann einfach mithilfe der Funktion function() definiert werden, und in einer Variable gespeichert werden, die genutzt werden kann, um die Funktion aufzurufen. Innerhalb der runden Klammern () von function() können Parameter definiert werden, die dann innerhalb der Funktion verfügbar sind. Beim Aufruf kann auch an selbst definierte Formeln mithilfe des Verkettungsoperators %>% ein Parameter übergeben werden. Die Möglichkeiten, die sich hierdurch ergeben sind schier endlos…
meine_cohens_d_funktion <- function(ergebnis_t.test) {
# Die Funktion unname() ist hier für eine schönere Ausgabe am Ende
# der Funktion nötig. Kann aber auch ignoriert werden.
t_statistik <- ergebnis_t.test$statistic %>% unname()
df_wert <- ergebnis_t.test$parameter %>% unname()
# Formel für Cohens d
# siehe auch https://www.uccs.edu/lbecker/
2 * t_statistik / sqrt(df_wert)
}
t.test(
dlp ~ untersuchung,
data = daten %>% filter(untersuchung %in% c("Abdomen", "Thorax")),
paired = FALSE
) %>% meine_cohens_d_funktion()## [1] 0.91628414.4 ANOVA
Was aber wenn wir nun nicht nur zwei Gruppen auf signifikante Unterschiede hin untersuchen wollen, sondern gleich mehrere? Der Gedanke einfach bspw. mehrere T-Tests hintereinander zu machen, läge natürlich nahe. Aber dabei ergäbe sich ein zentrales Problem: Wenn wir für unsere Test ein Signifikanzniveau von 0.05 ansetzen, die Wahrscheinlichkeit keinen Typ I Fehler zu machen also bei 95% liegt, würde sich bei einem Vergleich von nur drei Gruppen schon eine Wahrscheinlichkeit von knapp 14% ergeben mindestens eine Nullhypothese falsch abzulehnen.
Für genau solche Fälle eignen sich sogenannte Varianzanalysen (im Englischen auch als ANOVA = analysis of variance bezeichnet). Diese Verfahren testen im Grunde genommen die Nullhypothese, dass die Mittelwerte aller Gruppen gleich sind. In R kann eine solche Varianzanalyse sehr einfach mit der Funktion aov() durchgeführt werden, wobei der p-Wert hier nur ausgegeben wird, wenn man sich das Ergebnisobjekt genauer ansieht.
aov(dlp ~ untersuchung, data = daten) %>% summary()## Df Sum Sq Mean Sq F value Pr(>F)
## untersuchung 5 1.454e+09 290784001 408.5 <2e-16 ***
## Residuals 2435 1.733e+09 711883
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1In diesem Fall sehen wir als, dass wir unsere Nullhypothese verwerfen können, d.h. für das DLP gibt es in Abhängigkeit von der Untersuchungsart signifikante Unterschiede der Mittelwerte zwischen den Gruppen. Dieses Ergebnis ist für sich genommen natürlich wenig hilfreich, denn das wirklich interessante wäre nun zu wissen, welche Gruppen sich auf welche Art unterscheiden. Hierfür können sogenannte post-hoc Tests durchgeführt werden, die die Ergebnisse der Varianzanalyse weiter untersuchen und automatisch für das multiple Testen korrigieren. Ein Beispiel für einen solchen post-hoc Test ist der Tukey-Test, für den R die Funktion TukeyHSD() enthält.
aov(dlp ~ untersuchung, data = daten) %>% TukeyHSD()## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = dlp ~ untersuchung, data = daten)
##
## $untersuchung
## diff lwr upr p adj
## Becken-Abdomen -548.90104 -904.7938 -193.0082 0.0001650
## Gesicht-Abdomen -596.16468 -857.2855 -335.0439 0.0000000
## Kopf-Abdomen 644.24473 504.3822 784.1072 0.0000000
## Thorax-Abdomen -640.03148 -764.8412 -515.2218 0.0000000
## Traumaspirale-Abdomen 1787.49991 1615.1627 1959.8371 0.0000000
## Gesicht-Becken -47.26364 -469.0022 374.4749 0.9995580
## Kopf-Becken 1193.14576 833.6449 1552.6466 0.0000000
## Thorax-Becken -91.13045 -445.0468 262.7859 0.9776902
## Traumaspirale-Becken 2336.40094 1963.0654 2709.7365 0.0000000
## Kopf-Gesicht 1240.40941 974.3921 1506.4267 0.0000000
## Thorax-Gesicht -43.86680 -302.2873 214.5537 0.9967232
## Traumaspirale-Gesicht 2383.66459 2099.2287 2668.1005 0.0000000
## Thorax-Kopf -1284.27621 -1419.0301 -1149.5223 0.0000000
## Traumaspirale-Kopf 1143.25518 963.5854 1322.9250 0.0000000
## Traumaspirale-Thorax 2427.53139 2259.3137 2595.7491 0.0000000Das Ergebnis zeigt uns nun, dass außer für die Vergleiche Gesicht/Becken, Thorax/Becken und Thorax/Gesicht, jeweils signifikante Unterschiede der DLP-Mittelwerte bestehen (die p-Werte sind in der letzten Spalte zu finden). Zudem wird die mittlere Differenz sowie das 95%-Konfidenzintervall der Unterschiede mit ausgegeben.
4.4.1 Ein zweiter kurzer Ausfulg in Ergebnisobjekte
Gerade solche Ausgaben, wie die der TukeyHSD() Funktion sind manchmal etwas unübersichtlich, obwohl sie ja eigentlich sogar eine Tabelle mit ausgeben. Diese Tabelle aber wiederum weiter zu manipulieren, ist mit R Basisfunktionen kaum möglich. Erfreulicherweise gibt es auch hier aus dem Tidyverse einige praktische Helferfunktionen, die in dem Paket broom enthalten sind. Insbesondere die Funktion tidy() hilft Ergebnistabellen statistischer Tests in saubere Dataframes zu überführen, mit denen dann wie gewohnt weitergearbeitet werden kann.
library(broom)
# zunächst speichern wir das Ergebnis des Tukey-Tests, wobei natürlich auch eine
# direkte Verkettung zur tidy()-Funktion, usw. möglich wäre
tukey.ergebnis <- aov(dlp ~ untersuchung, data = daten) %>% TukeyHSD()
# Tukey Ergebnisse im tidy-Format
# (sieht in der Ausgabe hier online nicht so sehr anders aus,
# in RStudio ist der Unterschied deutlicher)
tukey.ergebnis %>%
tidy()## # A tibble: 15 × 7
## term contrast null.value estimate conf.low conf.high adj.p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 untersuchung Becken-Abdom… 0 -549. -905. -193. 1.65e-4
## 2 untersuchung Gesicht-Abdo… 0 -596. -857. -335. 1.23e-9
## 3 untersuchung Kopf-Abdomen 0 644. 504. 784. 0
## 4 untersuchung Thorax-Abdom… 0 -640. -765. -515. 0
## 5 untersuchung Traumaspiral… 0 1787. 1615. 1960. 0
## 6 untersuchung Gesicht-Beck… 0 -47.3 -469. 374. 1.00e+0
## 7 untersuchung Kopf-Becken 0 1193. 834. 1553. 0
## 8 untersuchung Thorax-Becken 0 -91.1 -445. 263. 9.78e-1
## 9 untersuchung Traumaspiral… 0 2336. 1963. 2710. 0
## 10 untersuchung Kopf-Gesicht 0 1240. 974. 1506. 0
## 11 untersuchung Thorax-Gesic… 0 -43.9 -302. 215. 9.97e-1
## 12 untersuchung Traumaspiral… 0 2384. 2099. 2668. 0
## 13 untersuchung Thorax-Kopf 0 -1284. -1419. -1150. 0
## 14 untersuchung Traumaspiral… 0 1143. 964. 1323. 0
## 15 untersuchung Traumaspiral… 0 2428. 2259. 2596. 0# sind die Tukey Ergebnisse einmal in ein sauberes Dataframe überführt, können wir
# alle Tidyverse-Funktionen in gewohnter Weise nutzen
# bspw. nur die signifkanten Ergebnisse anzeigen, nach Größe der mittleren Differenz
# sortiert, und das ganze dann graphisch aufbereitet
tukey.ergebnis %>%
tidy() %>%
filter(adj.p.value < 0.05) %>%
select(contrast, estimate, adj.p.value) %>%
ggplot(aes(x = fct_reorder(contrast, estimate), y = estimate, fill = adj.p.value)) +
geom_col() +
coord_flip() +
labs(x = "Vergleich",
y = "Mittlere Differenz DLP",
fill = "p-Wert") +
theme_bw()