pacman::p_load(tidyverse, magrittr, broom,
readxl, rstatix, coin,
effectsize, PMCMRplus, rcompanion)33 Der Friedman Test
Letzte Änderung am 02. April 2024 um 09:52:43
“It’s easy to lie with statistics. It’s hard to tell the truth without statistics.” — Andrejs Dunkels
Was macht der Friedman Test?
Der Friedman Test vergleicht die Mediane mehrerer beliebiger Verteilungen miteinander und berücksichtigt dabei auch einen Block oder Cluster. Im Prinzip die zweifaktorielle ANOVA für nicht-normalverteilte Daten.
Wann nutzen wir den Friedman Test? Wir nutzen den Friedman Test, wenn wir zwei Faktoren vorliegen haben und keine Normalverteilung annehmen können. Wichtig hierbei ist, dass wir zwei Faktoren vorliegen haben, die wir in unserem Modell berücksichtigen wollen. Du kannst den Friedman Test als das Äquivalent zu der zweifaktoriellen ANOVA sehen. Mit der Einschränkung, dass wir leider keine Interaktion berücksichtigen können.
In unserem Fall rechnen wir den Friedman-Test nicht per Hand. Es gibt dafür keinen Grund. Wir haben ja im vorherigen Kapitel uns mit der Berechnung des Kruskal-Wallis-Test beschäftigt, der Friedman-Test unterscheidet sich nur in Nuancen bei der Berechnung, so dass ich hier keinen Lerngewinn mehr sehe. Der Friedman-Test somit ein klarer Anwendungsfall. Wir nutzen den Friedman Test in unserer Abschlussarbeit, wenn wir zum Beispiel Boniturnoten auf der Likertskala auswerten wollen.
33.1 Genutzte R Pakete
Wir wollen folgende R Pakete in diesem Kapitel nutzen.
An der Seite des Kapitels findest du den Link Quellcode anzeigen, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst.
33.2 Daten für den Friedman Test
Für das Datenbeispiel bauen wir uns den Datensatz schnell selber zusammen. Wir wollen uns fünf Weizensorten \(A\) bis \(E\) anschauen und bonitieren, wie die Qualität der Pflanzen nach einer Dürreperiode aussieht. Daher vergeben wir einen Score auf der Likertskala von \(1\) bis \(9\), wobei \(9\) bedeutet, dass die Weizenpflanze vollkommen intakt ist. Wir haben unseren Versuch in vier Blöcken angelegt.
grade_tbl <- tibble(block = 1:4,
A = c(2,3,4,3),
B = c(7,9,8,9),
C = c(6,5,4,7),
D = c(2,3,1,2),
E = c(4,5,7,6)) |>
gather(key = variety, value = grade, A:E) |>
mutate(block = as_factor(block))Schauen wir uns nochmal unseren Datensatz an. Wir haben den Datensatz jetzt im Long-Format vorliegen und können dann gleich mit dem Datensatz weiterarbeiten.
grade_tbl# A tibble: 20 x 3
block variety grade
<fct> <chr> <dbl>
1 1 A 2
2 2 A 3
3 3 A 4
4 4 A 3
5 1 B 7
6 2 B 9
7 3 B 8
8 4 B 9
9 1 C 6
10 2 C 5
11 3 C 4
12 4 C 7
13 1 D 2
14 2 D 3
15 3 D 1
16 4 D 2
17 1 E 4
18 2 E 5
19 3 E 7
20 4 E 6In Abbildung 33.1 sehen wir nochmal die Visualisierung der Boniturnoten über die fünf Weizensorten und den vier Blöcken. Bei einer so kleinen Fallzahl von nur vier Blöcken pro Sorte entscheiden wir uns dann für einen Dotplot um uns die Daten einmal anzuschauen.
ggplot(grade_tbl, aes(variety, grade, fill = block)) +
theme_minimal() +
geom_dotplot(binaxis = "y", stackdir='center',
position=position_dodge(0.3)) +
scale_y_continuous(breaks = 1:9, limits = c(1,9))
Im Folgenden schauen wir usn nun einmal die Hypothesen des Friedman Tests an und rechen dann den Friedman Test in R einmal durch.
33.3 Hypothesen für den Friedman Test
Der Friedman Test betrachtet wie auch der Kruskal-Wallis-Test die Mediane \(\widetilde{y}\) und Ränge um einen Unterschied nachzuweisen. Daher haben wir in der Nullhypothese als Gleichheitshypothese. In unserem Beispiel lautet die Nullhypothese, dass die Mediane jedes Levels des Faktors variety gleich sind.
\[ H_0: \; \widetilde{y}_{A} = \widetilde{y}_{B} = \widetilde{y}_{C} = \widetilde{y}_{D} = \widetilde{y}_{E} \]
Die Alternative lautet, dass sich mindestens ein paarweiser Vergleich in den Medianen unterschiedet. Hierbei ist das mindestens ein Vergleich wichtig. Es können sich alle Mediane unterschieden oder eben nur ein Paar. Wenn ein Friedman Test die \(H_0\) ablehnt, also ein signifikantes Ergebnis liefert, dann wissen wir nicht, welche Mediane sich unterscheiden. Bein unseren fünf Weizensorten kommt da eine ganze Menge an Vergleichen zusammen. In der Folge nur ein Ausschnitt aller Alternativehypothesen.
\[ \begin{aligned} H_A: &\; \widetilde{y}_{A} \ne \widetilde{y}_{B}\\ \phantom{H_A:} &\; \widetilde{y}_{A} \ne \widetilde{y}_{C}\\ \phantom{H_A:} &\; ...\\ \phantom{H_A:} &\; \widetilde{y}_{D} \ne \widetilde{y}_{E}\\ \phantom{H_A:} &\; \mbox{für mindestens ein Paar} \end{aligned} \] Damit kommen uns die Hypothesen nicht so unbekannt vor. Die Hypothesenpaare sind die gleichen wie in einer ANOVA bzw. dem Kruskal-Wallis-Test.
33.4 Friedman Test in R
Der Friedman Test ist in R in verschiedenen Paketen implementiert. Wir nehmen die Standardfunktion friedman.test(). Wichtig ist wie wir in der Funktion das Modell definieren. Wir nutzen das Symbol | um die Behandlung von dem Block zu trennen. Vor dem | Symbol steht die Behandlung, hinter dem | steht der Block. Damit sieht das Modell etwas anders aus, aber im Prinzip ist die Definition des Modells einfach.
Das Tutorial von Salvatore S. Mangiafico zum Friedman Test liefert eine sehr ausführliche Anwendung über mehrere R Pakete hinweg.
friedman.test(grade ~ variety | block, data = grade_tbl)
Friedman rank sum test
data: grade and variety and block
Friedman chi-squared = 14.684, df = 4, p-value = 0.005403Nachdem wir die Funktion aufgerufen haben, erhalten wir auch gleich den \(p\)-Wert von \(0.005\) wieder. Da der \(p\)-Wert kleienr ist als das Signifikanzniveau \(\alpha\) von 5% können wir die Nullhypothese ablehnen. Es gibt mindestens einen paarweisen Unterschied. Das war auch anhand des Doplots zu erwarten. Nun stellt sich noch die Frage, wie groß der Effekt ist. Das können wir mit Kendalls \(W\) bestimmen. Auch die Funktionalität ist in der R Paket {effectsize} implementiert. Wir müssen nur wieder das ganze Modell in die Funktion kendalls_w() stecken.
kendalls_w(grade ~ variety | block, data = grade_tbl)Kendall's W | 95% CI
--------------------------
0.92 | [0.89, 1.00]
- One-sided CIs: upper bound fixed at [1.00].Nun kriegen wir ein Kendalls \(W\) von 0.92 raus. Wir immer ist die Frage nach der Interpretation. Dankenswerterweise gibt es auch die Funktion interpret_kendalls_w(), die uns mit Quelle ausgibt, wie stark der Effekt ist. Einfach nochmal für die Referenz in die Hilfeseite von der Funktion schauen.
interpret_kendalls_w(0.92)[1] "almost perfect agreement"
(Rules: landis1977)Wir haben also herausgefunden, dass wir einen starken Effekt haben und sich die Mediane der Gruppen unterscheiden.
33.5 Posthoc Test
Nachdem wir jetzt festgestellt haben, dass sich mindestens ein Gruppenunterschied zwischen den Weizensorten finden lassen muss, wollen wir noch feststellen wo dieser Unterschied liegt. Wir nutzen dafür den Siegel- und Castellan-Vergleichstests für alle Paare durch. Der Test ist mit der Funktion frdAllPairsSiegelTest() in dem R Paket {PMCMRplus} implementiert und einfach nutzbar. In dem Paket sind noch eine weitere Reihe an statistischen Test für paarweise nicht-parametrische Test enthalten. Leider sind das Paket und die Funktion schon älter, so dass wir nicht einmal die Formelschreibweise zu Verfügung haben. Wir müssen alle Spalten aus unserem Datensatz mit dem $ einzeln selektieren und den Optionen zuordnen.
siegel_test_res <- frdAllPairsSiegelTest(y = grade_tbl$grade,
groups = grade_tbl$variety,
blocks = grade_tbl$block,
p.adjust.method = "none")
siegel_test_res A B C D
B 0.0052 - - -
C 0.1461 0.1797 - -
D 0.5762 0.0008 0.0442 -
E 0.1797 0.1461 0.9110 0.0573Ich empfehle ja immer eine Adjustierung für multiple Vergleiche, aber du kannst das selber entscheiden. Wenn du für die multiplen Vergleiche adjustieren willst, dann nutze gerne die Option p.adjust.method = "bonferroni" oder eben statt bonferroni die Adjustierungsmethode fdr.
Damit wir auch das compact letter display aus unseren \(p\)-Werte berechnen können, müssen wir unser Ergebnisobjekt nochmal in eine Tabelle umwandeln. Dafür gibt es dann auch eine passende Funktion mit PMCMRTable. Wir sehen, es ist alles etwas altbacken.
siegel_test_tab <- PMCMRTable(siegel_test_res )
siegel_test_tab Comparison p.value
1 B - A = 0 0.00519
2 C - A = 0 0.146
3 D - A = 0 0.576
4 E - A = 0 0.18
5 C - B = 0 0.18
6 D - B = 0 0.000796
7 E - B = 0 0.146
8 D - C = 0 0.0442
9 E - C = 0 0.911
10 E - D = 0 0.0573Nachdem wir die Tabellenschreibweise der \(p\)-Werte vorliegen haben, können wir dann das compact letter display uns über die Funktion cldList() anzeigen lassen. Leider gibt es keine Möglichkeit die 95% Konfidenzintervalle zu erhalten. Wir auch beim Kruskal-Wallis-Test erhalten wir bei dem Friedman Test keine 95% Konfidenzintervalle und müssen im Zweifel mit dieser Einschränkung leben.
cldList(p.value ~ Comparison, data = siegel_test_tab) Group Letter MonoLetter
1 B a a
2 C ab ab
3 D c c
4 E abc abc
5 A bc bcAm compact letter display sehen wir im Prinzip das gleiche Muster wie in den \(p\)-Werten. Nur nochmal etwas anders dargestellt und mit dem Fokus auf die nicht Unterschiede. Wenn du mehr über das compact letter display wissen willst dann lese gerne nochmal im Kapitel 34.7 nach.