/
stat-tests-diagnostic.qmd
278 lines (197 loc) · 19.4 KB
/
stat-tests-diagnostic.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
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
```{r echo = FALSE}
pacman::p_load(tidyverse, readxl, knitr, kableExtra, Hmisc, see,
openxlsx)
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73",
"#F0E442", "#0072B2", "#D55E00", "#CC79A7")
```
# Der diagnostische Test {#sec-test-diag}
*Letzte Änderung am `r format(fs::file_info("stat-tests-diagnostic.qmd")$modification_time, '%d. %B %Y um %H:%M:%S')`*
> *"Im Zweifel für den Zweifel; Das Zaudern und den Zorn." --- Tocotronic, Im Zweifel für den Zweifel*
::: callout-note
## Was macht der diagnostische Test?
Ein diagnostischer Test gibt die Güte einer Klassifizierung wieder. Wir gut wurden Kranke als krank und Gesunde als gesund erkannt? Diese und ähnliche Fragen beantwortet das diagnostische Testen.
:::
In diesem Kapitel wollen wir uns mit dem diagnostischen Test beschäftigen. Eigentlich würde man meinen, dass Diagnostik nun eher in die Medizin gehört. Hauptsächlich wird das diagnostische Testen auch in der Labordiagnostik oder bei der Herstellung eines neuen medizinischen Tests verwendet. Deshalb gibt es hier auch erstmal (Stand Ende 2022) nur einen Ausschnitt aus dem diagnostischen Testen. Wir brauchen aber die Fachbegriffe, wenn wir uns später einmal mit dem maschinellen Lernen oder dem Klassifizieren beschäftigen wollen. Dann brauchen wir die hier verwendeten Fachbegriffe, wie Spezifität, Sensitivität, AUC und auch die ROC Abbildung.
## Genutzte R Pakete
Wir wollen folgende R Pakete in diesem Kapitel nutzen.
```{r echo = TRUE}
pacman::p_load(tidyverse, magrittr,
pROC, readxl)
```
An der Seite des Kapitels findest du den Link *Quellcode anzeigen*, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst.
## Die Daten für das diagnostische Testen
Die Datentabelle für das diagnostische Testen basiert wie der $\chi^2$-Test auf der 2x2 Kreuztabelle oder Verfeldertafel. Es ist dabei *unabdingbar*, dass oben links in der 2x2 Kreuztabelle immer die $T^+$ und $K^+$ Werte stehen. Sonst funktionieren alle Formeln in diesem Kapitel nicht. Wir schauen uns auch immer *das Schlechte an*. Daher wollen wir immer wissen, ist der Hund krank? Ist der Hund tot? Ist der Hund weggelaufen? Beide Voraussetzung sind wichtig, damit wir mit der 2x2 Kreuztabelle wie in @tbl-2x2-table-general gezeigt rechnen können.
| | | | | |
|:--------:|:---------:|:-----------------:|:-----------------:|:--------------:|
| | | **Krank** | | |
| | | $K^+$ (1) | $K^-$ (0) | |
| **Test** | $T^+$ (1) | $TP_{\;\Large a}$ | $FN_{\;\Large b}$ | $\mathbf{a+b}$ |
| | $T^-$ (0) | $FP_{\;\Large c}$ | $TN_{\;\Large d}$ | $\mathbf{c+d}$ |
| | | $\mathbf{a+c}$ | $\mathbf{b+d}$ | $\mathbf{n}$ |
: Eine 2x2 Tabelle oder Vierfeldertafel {#tbl-2x2-table-general}
Wir wollen die @tbl-2x2-table-general mit einem Beispiel von *Tollwut* an Hauskatzen in *ländlicher* Umgebung. Die Katzen haben also Auslauf und können sich auch mit Tollwut infizieren. Wir wollen einen neuen, nicht invasiven Labortesten *Tollda* darauf überprüfen, wie gut der diagnostische Test Tollwut bei Katzen im Frühstadium erkennt.
Wir haben jetzt folgende Informationen erhalten:
- Der diagnostische Test *TollDa* ist positiv $T^+$, wenn Tollwut vorliegt $K^+$ , in 80% der Fälle.
- Der diagnostische Test *TollDa* ist positiv $T^+$, wenn *keine* Tollwut vorliegt $K^-$, in 9.5% der Fälle.
- Abschließend haben noch 2% der Katzen in ländlicher Umgebung Tollwut. Wir haben eine Prävalenz der Erkrankung in der betrachteten Population von 2%.
Die Halterin einer Katze möchte nun wissen, wie groß ist die Wahrscheinlichkeit bei einem positiven Testergebnis, dass meine Katze Tollwut hat. Also die bedingte Wahrscheinlichkeit $Pr(K^+|T^+)$ oder die Wahrscheinlichkeit krank zu sein gegeben eines positiven Tests.
![Visualisierung der Informationen zur Tollwutdiagnostik in einem Doppelbaum. Gefragt ist nach der Wahrscheinlichkeit $Pr(K^+|T^+)$ oder die Wahrscheinlichkeit krank zu sein gegeben eines positiven Tests.](images/diagnostic/diag_testen_doppel.png){#fig-stat-diag-01 fig-align="center" width="80%"}
@fig-stat-diag-01 visualisiert unsere Frage in einem Doppelbaum. Wir haben $10000$ Katzen vorliegen. Ja wir waren fleißig und wir können damit besser rechnen. Du siehst aber auch, für die Diagnostik brauchen wir eine Menge Beobachtungen.
Von den $10000$ Katzen sind 2% also $200$ wirklich mit Tollwut infiziert, also haben den Status $K^+$. Damit sind $9800$ Katzen gesund oder nicht krank und haben den Status $K^-$. Wir wissen jetzt, dass 80% der $K^+$ Katzen als positiv vom Test erkannt werden. Damit werden $200 \cdot 0.8 = 160$ Katzen $T^+$. Im Umkehrschluss sind die anderen $40$ Katzen dann $T^-$. Von den $9800$ gesunden Katzen werden 9.5% falsch als krank erkannt, also $9800 \cdot 0.095 = 931$ Katzen. Wiederum im Umkehrschluss sind dann $9069$ richtig als gesunde Tiere erkannt.
Wir können nun den Doppelbaum nach unten ergänzen. Wir haben damit $1091 = 160 + 931$ positiv getestete Katzen $T^+$ sowie $9109 = 40 + 9069$ negativ getestete Katzen $T^-$.
Wie groß ist nun die Wahrscheinlichkeit $Pr(K^+|T^+)$ oder die Wahrscheinlichkeit krank zu sein gegeben eines positiven Tests? Wir können die Wahrscheinlichkeit $Pr(K^+|T^+)$ direkt im Baum ablesen. Wir haben $160$ kranke und positive Tier. Insgesamt sind $1091$ Tiere positiv getestet. Daher Wahrscheinlichkeit krank zu sein gegeben eines positiven Tests $\tfrac{160}{1091} = 0.147$ oder nur 14%.
Sammeln wir unsere Informationen um die @tbl-2x2-table-general zu füllen. Wir haben insgesamt $n = 10000$ Katzen untersucht. In @tbl-2x2-table-example sehen wir die Ergebnisse unseres Tests auf Tollwut zusammengefasst. Wir haben $160$ Katzen die Tollwut haben und vom Test als krank erkannt wurden. Dann haben wir $931$ Katzen, die keine Tollwut hatten und vom Test als positiv erkannt wurden. Darüber hinaus haben wir $40$ Katzen, die Tollwut haben, aber ein negatives Testergebnis. Sowie $8109$ gesunde Katzen mit einem negativen Testergebnis.
| | | | | |
|:--------:|:---------:|:--------------------:|:---------------------:|:---------------------:|
| | | **Krank** | | |
| | | $K^+$ (1) | $K^-$ (0) | |
| **Test** | $T^+$ (1) | 160 | 931 | $\mathbf{a+b} = 1091$ |
| | $T^-$ (0) | 40 | 9069 | $\mathbf{c+d} = 9109$ |
| | | $\mathbf{a+c} = 200$ | $\mathbf{b+d} = 9800$ | $\mathbf{n} = 10000$ |
: Eine 2x2 Tabelle oder Vierfeldertafel gefüllt mit dem Beispiel aus dem Doppelbaum. {#tbl-2x2-table-example}
Wie du sehen kannst ist die mittlere Reihe des Doppelbaums nichts anderes als die ausgefüllte 2x2 Kreuztabelle. In @fig-stat-diag-00 sehen wir die letzte Visualisierung des Zusammenhangs von Testscore auf der x-Achse und der Verteilung der tollwütigen Katzen $K^+$ und der gesunden Katzen $K^-$. Links und rechts von der Testenstscheidung werden Katzen falsch klassifiziert. Das heißt, die Katzen werden als krank oder gesund von dem Test erkannt, obwohl die Katzen diesen Status nicht haben. Ein idealer Test würde zwei Verteilungen der kranken und gesunden Tiere hervorbringen, die sich perfekt separieren lassen. Es gibt anhand des Testscores keine Überlappung der beiden Verteilungen
![Visualisierung des Zusammenhangs zwischen Testscore und den Verteilungen der kranken und gesunden Katzen sowie der Testenstscheidung. Links und rechts von der Testenstscheidung werden Katzen falsch klassifiziert.](images/diagnostic/diag_testen_dist.png){#fig-stat-diag-00 fig-align="center" width="80%"}
## Confusion matrix (deu. *Fehlermatrix*)
In der Statistik beziehungsweise der Epidemiologie gibt es eine Vielzahl an statistischen Maßzahlen für die 2x2 Kreuztabelle. In der Fachsprache wird die 2x2 Kreuztabelle auch [Confusion matrix](https://en.wikipedia.org/wiki/Confusion_matrix) genannt. Auf der Confusion Matrix können viele Maßzahlen berechnet werden wir konzentrieren uns hier erstmal auf die zwei wichtigsten Maßzahlen, nämlich der Spezifität und der Sensitivität.
In der wissenschaftlichen Fachsprache hat ein diagnostischer Test oder eine Methode, die erkrankte Personen sehr zuverlässig als krank ($K^+$) erkennt, eine hohe Sensitivität. Das heißt, sie übersieht kaum erkrankte Personen. Ein Test, der gesunde Personen zuverlässig als gesund ($K^-$) einstuft, hat eine hohe Spezifität. Wir können die Spezifität und die Sensitivität auf der 2x2 Kreuztabelle wir folgt berechnen.
$$
\mbox{Sensitivität} = \mbox{sens} = \cfrac{TP}{TP + FN} = \cfrac{160}{160 + 931} = 0.147
$$
$$
\mbox{Spezifität} = \mbox{spec} = \cfrac{TN}{TN + FP} = \cfrac{8869}{8869 + 40} = 0.995
$$
Beide statistischen Maßzahlen benötigen wir im Besonderen bei der Erstellung der Reciever Operator Curve (ROC) im folgenden Abschnitt.
## Receiver Operating Characteristic (ROC)
Die [Reciever Operator Curve (ROC)](https://en.wikipedia.org/wiki/Receiver_operating_characteristic) ist *die* Abbildung wenn es darum geht darzustellen wie gut eine Klassifikation funktioniert hat. Wir werden die Abbildung später im maschinellen Lernen und der Klssifikation wiedertreffen. Da die ROC Kurve aber ursprünglich zum diagnostischen Testen gehört, kommt die ROC Kurve hier auch rein.
### Daten für die Receiver Operating Characteristic (ROC)
Schauen wir uns erstmal ein Beispiel für die ROC Kurve an. In @tbl-data-roc-1 sehen wir Beispieldaten von vierzehn Hunden. Von den vierzehn Hunden haben sieben eine verdeckte Flohinfektion im Anfangsstadium und sieben Hunde sind flohfrei und gesund. Daher haben wir sieben kranke Hunde ($K^+$ oder Infektionsstatus ist $1$) und sieben gesunde Hunde ($K^-$ oder Infektionsstatus ist $0$). Wir *testen* jeden der Hunde mit dem Flohindikatormittel *FleaDa*. Das Flohindikatormittel gibt einen Farbwert als `test_score` zwischen $0$ und $1$ wieder. Wir wollen nun herausfinden, ob wir mit dem `test_score` die vierzehn Hunde korrekt nach ihrem Krankheitszustand klassifizieren können.
```{r}
#| message: false
#| echo: false
#| tbl-cap: Datensatz für vierzehn Hunde mit Flohinfektion und keiner Flohinfektion sowie dem Testscore des Flohindikatormittels *fleaDa*.
#| label: tbl-data-roc-1
roc_tbl <- tibble("1" = c(0.1, 0.15, 0.25, 0.3, 0.4, 0.5, 0.6),
"0" = c(0.35, 0.425, 0.47, 0.55, 0.7, 0.8, 0.82)) |>
gather(infected, test_score) |>
mutate(infected = as.numeric(infected))
roc_tbl |> kable(align = "c", "pipe")
```
Wir sehen die Daten visualisiert in @fig-labels-roc1. Auf der $y$-Achse ist der binäre Endpunkt $infected$ und auf der $x$-Achse der Testscore des Flohindikators *FleaDa*.
```{r}
#| warning: false
#| echo: false
#| message: false
#| label: fig-labels-roc1
#| fig-align: center
#| fig-height: 3
#| fig-width: 7
#| fig-cap: "Aufteilung der kranken und gesunden Hunde nach Infektionsstatus und dem Testscore."
roc_tbl <- tibble("1" = c(0.1, 0.15, 0.25, 0.3, 0.4, 0.5, 0.6),
"0" = c(0.35, 0.425, 0.47, 0.55, 0.7, 0.8, 0.82)) |>
gather(infected, test_score) |>
mutate(infected = as.numeric(infected)) |>
dplyr::arrange(test_score)
write.csv2(roc_tbl, "data/roc_data.csv", row.names = FALSE)
write.xlsx(roc_tbl, "data/roc_data.xlsx", rowNames = FALSE)
roc_tbl |>
mutate(infected = factor(infected, levels = c(0, 1))) |>
ggplot(aes(infected, test_score, shape = infected)) +
geom_point(size = 3) +
theme_minimal() +
coord_flip() +
labs(x = "Status der Infektion", y = "Testscore") +
theme(legend.position = "none",
axis.text.x = element_text(face="bold", size=10),
axis.text.y = element_text(face="bold", size=10))
```
```{r}
#| warning: false
#| echo: false
#| message: false
roc_res_tbl <- roc(infected ~ test_score, roc_tbl,
direction = ">")
kable_tbl <- roc_res_tbl |>
coords(ret = "all") |>
select(threshold:fp) |>
mutate_if(is.numeric, round, 2) |>
dplyr::arrange(threshold)
```
Die Idee der ROC Kurve ist nun für jeden möglichen Wert des Testscores die Spezifität und Sensitivität zu berechnen. Wir müssen das aber nicht für alle Werte des Testscores machen, sondern nur für die Wert bei denen sich der Status einer Beobachtung anhand des Testscores ändern würde. Klingt ein wenig schräg, schauen wir es uns einmal an. Zuerst brauchen wir jeweils die Grenzen an denen wir für den Testscore eine Klassifikation machen. In @fig-labels-roc3 sehen wir die Grenzen als gelbe Linie.
```{r}
#| warning: false
#| echo: false
#| message: false
#| label: fig-labels-roc3
#| fig-align: center
#| fig-height: 3
#| fig-width: 8
#| fig-cap: "Aufteilung der kranken und gesunden Hunde nach Infektionsstatus und dem Testscore."
roc_tbl |>
mutate(infected = factor(infected, levels = c(0, 1))) |>
ggplot(aes(infected, test_score, shape = infected)) +
geom_point(size = 3) +
geom_hline(yintercept = kable_tbl$threshold[-1], color = cbbPalette[2]) +
scale_y_continuous(breaks = kable_tbl$threshold[-1]) +
theme_minimal() +
coord_flip() +
labs(x = "Status der Infektion", y = "Testscore") +
theme(legend.position = "none",
axis.text.x = element_text(face="bold", size=10),
axis.text.y = element_text(face="bold", size=10))
```
Wir können jetzt für jede gelbe Linie als Threshold des Testscore die Spezifität und die Sensitivität berechen. berechnen @tbl-data-roc-2 sind die Werte von Spezifität und Sensitivität für jede gelbe Linie ausgegeben. Wenn wir als Entscheidungsgrenze einen Testscore von 0.12 nehmen würden, dann würden wir 1 Hund als krank und 13 als gesund klassifizieren. Wir hätten 7 TN, 1 TP, 6 FN und 0 FP. Aus diesen Zahlen, die eine 2x2 Kreuztabelle entsprechen, können wir dann die Spezifität und die Sensitivität berechnen. Wir berechnen dann die Spezifität und die Sensitivität für jeden Threshold.
```{r}
#| message: false
#| echo: false
#| tbl-cap: Die Spezifität und die Sensitivität berechnet für jeden Threshold.
#| label: tbl-data-roc-2
kable_tbl[c(-1, -15),] |>
select(-accuracy) |>
kable(align = "c", "pipe")
```
Wir können jetzt für jeden Threshold und damit jedes verbundene Spezifität und Sensitivität Paar einen Punkt in einen Plot einzeichnen. Wir lassen damit Stück für Stück die ROC Kurve wachsen.
In der @fig-roc-drawn-1 sehen wir die ROC Kurve mit dem Threshold von 0.32 vorliegen. Wir haben damit 7 TN, 4 TP, 3 FN und 0 FP klassifizierte Beobchtungen. Darauf ergibt sich eine Spezifität von 1 und eine Sensitivität von 0.57. In der ROC Kurve tragen wir auf die $x$-Achse die 1 - Spezifitätswerte auf. Damit haben wir eine schön ansteigende Kurve.
![Visualisierung des Anwachsen der ROC Kurve mit einem Threshold von 0.32 und damit 7 TN, 4 TP, 3 FN und 0 FP.](images/diagnostic/diag-roc-01.png){#fig-roc-drawn-1 fig-align="center" width="100%"}
In der @fig-roc-drawn-2 istr die ROC Kurve schon nach rechts gewachsen, da wir die ersten falsch positiv (FP) klassifizierten Beobachtungen bei einem Threshold von 0.48 erhalten.
![Visualisierung des Anwachsen der ROC Kurve mit einem Threshold von 0.48 und damit 4 TN, 5 TP, 2 FN und 3 FP.](images/diagnostic/diag-roc-02.png){#fig-roc-drawn-2 fig-align="center" width="100%"}
In @fig-roc-drawn-3 sind wir mit einem Threshold von 0.65 schon fast am rechten Rand der Verteilung des Testscores angekommen. Wir haben nur noch 3 richtig negative (TN) Beobachtungen, die jetzt mit steigenden Threshold alle zu falsch positiven (FP) klassifiziert werden.
![Visualisierung des Anwachsen der ROC Kurve mit einem Threshold von 0.65 und damit 3 TN, 7 TP, 0 FN und 4 FP.](images/diagnostic/diag-roc-03.png){#fig-roc-drawn-3 fig-align="center" width="100%"}
Die ROC Kurve endet dann oben rechts in der Abbildung und startet immer unten links. Wir können noch die Fläche unter der ROC Kurve berechnen. Wir nennen diese Fläche unter der Kurve AUC (eng. *area under the curce*) und wir brauchen diese Maßzahl, wenn wir verschiedene Testscores und die entsprechenden ROC Kurven miteinander vergleichen wollen. Eine AUC von 0.5 bedeutet, dass unser Testscore nicht in der Lage war die kranken von den gesunden Hunden zu trennen. Ein AUC von 1 bedeutet eine perfekte Trennung der beiden Gruppen durch den Testscore. Ein AUC von 0.5 bedeutet somit auch keine Trennung der beiden Gruppen durch den Testscore.
## Reciever Operator Curve (ROC) in R
Wir berechnen die ROC Kurve nicht per Hand. Das wäre eine sehr große Fleißarbeit für jeden Threshold zwischen den Beobachtungen die jeweilige Spezifität und Sensitivität Paare zu berechnen. Wir nutzen daher das R Paket `{pROC}`. Das Pakt erlaubt es uns später auch recht einfach ROC Kurven miteinander zu vergleichen und einen statistischen Test zu rechnen.
Die Daten von den obigen Beispiel sind in der Datei `roc_data.xlsx` gespeichert und können dann einfach verwendet werden. Die Funktion `roc()` akzeptiert ein binäres $y$ wo bei die $1$ das Schlechte bedeutet und die $0$ die Abwesenheit von dem Schlechten. Also eben krank oder gesund wie oben beschrieben. Wir wollen dann noch die AUC und die 95% Konfidenzintervalle für das AUC wiedergegeben haben.
```{r}
#| warning: false
#| echo: true
#| message: false
roc_tbl <- read_excel("data/roc_data.xlsx")
roc_obj <- roc(infected ~ test_score, roc_tbl, auc = TRUE, ci = TRUE)
roc_obj
```
Wir sehen, dass wir eine AUC von 0.837 haben, was auf eine gute Trennschärfe des Tests hindeutet. Wäre die AUC näher an 0.5 dann würden wir sagen, dass der Test die kranken Hunde nicht von den gesunden Hunden unterscheiden kann.
Im nächsten Schritt extrahieren wir noch alle wichtigen Informationen aus dem Objekt `roc_obj` damit wir die ROC Kurve in `{ggplot}` zeichnen können. Das Paket {pROC} hat auch eine eigne Plotfunktion. Es gibt eine Reihe von [Screenshots vom `{pROC}` Paket](https://web.expasy.org/pROC/screenshots.html), wo du noch andere Möglichkeiten siehst. Einfach mal ausprobieren.
```{r}
#| warning: false
#| echo: true
#| message: false
roc_res_tbl <- roc_obj |>
coords(ret = "all") |>
select(specificity, sensitivity)
```
Das Objekt `roc_obj` nutzen wir jetzt um die ROC Kurve einmal darzustellen. Achte drauf, dass auf der $x$-Achse die 1-Spezifität Werte stehen. Wir erhalten die ROC Kurve wie in @fig-labels-roc2 dargestellt. Wie zeichnen noch die Diagonale ein. Wenn die ROC nahe an der Diagonalen verläuft, dann ist der Testscore nicht in der Lage die Kranken von den Gesunden zu trennen. Eine perfekte ROC Kurve läuft senkrecht nach oben und dann waagerecht nach rechts und es ergibt sich eine AUC von 1.
```{r}
#| warning: false
#| echo: true
#| message: false
#| label: fig-labels-roc2
#| fig-align: center
#| fig-height: 5
#| fig-width: 5
#| fig-cap: "Abbildung der ROC Kurve."
ggplot(roc_res_tbl, aes(1-specificity, sensitivity)) +
geom_path() +
geom_abline(intercept = 0, slope = 1, alpha = 0.75, color = "red") +
theme_minimal() +
labs(x = "1 - Spezifität", y = "Sensitivität")
```
In späteren Kapiteln zur Klassifikation werden wir auch verschiedene ROC Kurven zusammen darstellen und miteinander Vergleichen. Dazu nutzen wir dann auch die ROC Kurve und die Möglichkeit auch einen `roc.test()` zu rechnen. Mehr dazu gibt es soweit erstmal auf der Hilfeseite des R Paketes `{pROC}` mit vielen [Screenshots vom `{pROC}` Paket](https://web.expasy.org/pROC/screenshots.html). Für dieses Kapitel soll die Darstellung einer ROC Kurve genügen.