/
stat-tests-friedman.qmd
158 lines (114 loc) · 9.25 KB
/
stat-tests-friedman.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
```{r echo = FALSE}
pacman::p_load(tidyverse, readxl, knitr, kableExtra, Hmisc)
```
# Der Friedman Test {#sec-friedman}
*Letzte Änderung am `r format(fs::file_info("stat-tests-friedman.qmd")$modification_time, '%d. %B %Y um %H:%M:%S')`*
> *"It's easy to lie with statistics. It's hard to tell the truth without statistics." --- Andrejs Dunkels*
::: callout-note
## 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.
## Genutzte R Pakete
Wir wollen folgende R Pakete in diesem Kapitel nutzen.
```{r echo = TRUE}
pacman::p_load(tidyverse, magrittr, broom,
readxl, rstatix, coin,
effectsize, PMCMRplus, rcompanion)
```
An der Seite des Kapitels findest du den Link *Quellcode anzeigen*, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst.
## 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.
```{r}
#| message: false
#| warning: false
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.
```{r}
grade_tbl
```
In @fig-fried-ggplot 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.
```{r}
#| echo: true
#| message: false
#| label: fig-fried-ggplot
#| fig-align: center
#| fig-height: 5
#| fig-width: 6
#| fig-cap: "Dotplot des Datenbeispiels für die Bonitur von fünf Weizensorten."
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.
## 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.
## 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.
::: column-margin
Das [Tutorial von Salvatore S. Mangiafico zum Friedman Test](https://rcompanion.org/handbook/F_10.html) liefert eine sehr ausführliche Anwendung über mehrere R Pakete hinweg.
:::
```{r}
#| message: false
#| warning: false
friedman.test(grade ~ variety | block, data = grade_tbl)
```
Nachdem 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.
```{r}
#| message: false
#| warning: false
kendalls_w(grade ~ variety | block, data = grade_tbl)
```
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.
```{r}
#| message: false
#| warning: false
interpret_kendalls_w(0.92)
```
Wir haben also herausgefunden, dass wir einen starken Effekt haben und sich die Mediane der Gruppen unterscheiden.
## 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.
```{r}
#| message: false
#| warning: false
siegel_test_res <- frdAllPairsSiegelTest(y = grade_tbl$grade,
groups = grade_tbl$variety,
blocks = grade_tbl$block,
p.adjust.method = "none")
siegel_test_res
```
Ich 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.
```{r}
siegel_test_tab <- PMCMRTable(siegel_test_res )
siegel_test_tab
```
Nachdem 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.
```{r}
cldList(p.value ~ Comparison, data = siegel_test_tab)
```
Am *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 @sec-compact-letter nach.