37 Permutationstest & Bootstraping

Letzte Änderung am 02. April 2024 um 09:52:43

“Wer zu spät kommt, den bestraft das Leben.” — Gorbatschow zu Honecker beim Staatsbesuch in der DDR am 7. Oktober 1989

Dieses Kapitel ist teilweise eine Baustelle und wird es vermutlich auch über das Sommersemester 2024 bleiben. Daher können Teile des Codes und des Textes kurzfristig keinen Sinn ergeben oder nicht funktional sein. Es ist geplant eine fertige Version im Sommer 2024 vorliegen zu haben.

In diesem Kapitel wollen wir ein mächtiges Verfahren kennen lernen, welches zu den altvorderen Zeiten von Fisher und Co. nicht zu Verfügung stand: der Permutationstest oder aber das Bootstraping. Dabei war die Theorie schon bekannt, aber die technische Durchführung war unmöglich. Erst seit der Entwicklung von Computern, die hunderte bis tausende Rechenschritte in wenigen Sekunden durchführen können, sind Permutationstests erst denkbar. Damit erlauben uns die heutigen technischen Gegebenheiten mehr als es vor über hundert Jahren der Fall war, oder um es mit Fay und Gerow (2018) zu sagen:

“[…], if statistics were being invented today, bootstrapping, as well as related re-sampling methods, would be the standard go-to approach.” — Fay und Gerow (2018)

Was macht also Permutationstests oder allgemeiner gesprochen resampling Methoden so wirkmächtig? Permutationstests beruhen nicht auf Annahmen über die Verteilung unseres Outcomes \(y\), wie es bei einigen anderen Tests, wie zum Beispiel dem T-test, der Fall ist. Permutationstests funktionieren, indem die beobachteten Daten mehrmals neu gezogen werden, um einen \(p\)-Wert für den Test zu bestimmen. Wir haben es also mit einem iterativen Prozess zu tun. Oder anders gesagt, wir simulieren hier Daten und damit die Nullhypothese. Erinnern wir uns daran, dass der \(p\)-Wert als die Wahrscheinlichkeit definiert ist, Daten zu erhalten, die genauso extrem sind wie die beobachteten Daten, wenn die Nullhypothese wahr ist. Wenn die Daten gemäß der Nullhypothese mehrmals gemischt werden, kann die Anzahl der Fälle mit Daten, die genauso extrem sind wie die beobachteten Daten, gezählt und ein \(p\)-Wert berechnet werden.

Die Vorteile von Permutationstests, neben der relative Einfachheit bei Durchführung und Interpretation, sind:

das Fehlen von Annahmen über die Verteilung der zugrunde liegenden Daten. Wir müssen eigentlich uns keine Gedanken über die Verteilung des Outcomes \(y\) machen.
ihre Flexibilität bei der Verarbeitung verschiedener Arten von Daten wie nominal, ordinal oder intervall-/ratioskaliert. Daher geht prinzipiell alles mit einem Permutationstest.

Die Nachteile von Permutationstests sind:

die begrenzte Komplexität der Designs, die sie verarbeiten können. Das heißt, es geht viel, aber wenn unser experimentelles Design genested ist oder aber eine komplexe faktorielle Struktur hat, geht es leider nicht mehr.
und die Unbekanntheit vieler Anwender mit dieser Methode. Der Permutationstest ist nicht sehr bekannt. Daher wird der Permutationstest auch immer etwas zwiespältig gesehen. Der Anwender baut ja seinen Test selber und nutzt nicht einen bekannten Algorithmus.

In diesem Kapitel gebe ich eine kurze Übersicht über den Permutationstest und andere resampling Methoden. Es gibt wie immer die ein oder andere Literatur und auch Kritik sowie Diskussionen. Ich finde den Permutationstest als eine gute Lösung, wenn i) genug Fallzahl vorliegt und damit mehr als 10 Beobachtungen pro Gruppe sowie ii) es keinen passenden Test gibt. Warum muss es genug Beobachtungen in den Gruppen bzw. dann in dem Datensatz geben? Wir müssen ja genug zufällige Daten permutieren. Wenn wir zu wenig Beobachtungen haben, dann wiederholen sich unsere zufälligen Datensätze sehr schnell. Wir laufen dann in das gleiche Problem wie bei den nicht-parametrischen Tests, die auch eine Unterklasse von Permutationsverfahren mit einer geschlossenen mathematischen Formel sind.

Weitere Tutorien für den Permutationstest & Bootstrap

Wir immer geht natürlich mehr als ich hier Vorstellen kann. Du findest im Folgenden Tutorien, die mich hier in dem Kapitel inspiriert haben. Ich habe mich ja in diesem Kapitel auf die Durchführbarkeit in R und die allgemeine Verständlichkeit konzentriert. Es geht aber natürlich wie immer auch mathematischer…

Ein umfangreiches Tutorium von Ben Bolker Simple permutation tests in R. Eine der ersten Anlaufstellen, wenn man sich mit dem Thema Permutatuonstests in R beschäftigen will.
Das R Paket {infer} ist ein guter Startpunkt sich einmal klar zu machen, wie die Permutationstest funktionieren. Das Oaket hat ein sehr sauberes Framework und ich stelle hier auch noch ein Beispiel vor.
Fear not the bootstrap aus Fay und Gerow (2018) ist ein kurzen Text über den Bootstrap und was die Idee dahinter ist. Das ganze Buch ist zu empfehlen, aber der Abschnitt aus dem Buch ist kurz und übersichtlich.
Das R Paket {resample} liefert dann die aktuellste Form in R dann alle möglichen resampling Methoden durchzuführen. Ein Beispiel werde ich dann wieder vorstellen, aber es lohnt sich die Seite einmal zu besuchen.

37.1 Genutzte R Pakete

Wir wollen folgende R Pakete in diesem Kapitel nutzen.

pacman::p_load(tidyverse, magrittr, conflicted, broom)
conflicts_prefer(dplyr::filter)

An der Seite des Kapitels findest du den Link Quellcode anzeigen, über den du Zugang zum gesamten R-Code dieses Kapitels erhältst.

37.2 Daten

Für unsere Demonstration des Permutationstest nutzen wir den Datensatz flea_dog_cat_length_weight.xlsx. Zum einen wollen wir einen Gruppenvergleich zwischen den Sprungweiten der Hunde- und Katzenflöhe rechnen. Also einen klassischen t-Test für einen Gruppenvergleich. Nur eben hier als einen Permutationstest. Als zweites wollen wir einen \(p\)-Wert für das Bestimmtheitsmaß \(R^2\) abschätzen. Per se gibt es keinen \(p\)-Wert für das Bestimmtheitsmaß \(R^2\), aber der Permutationstest liefert hier eine Lösung für das Problem. Daher schauen wir uns in einer simplen linearen Regression den Zusammenhang zwischen einem \(y\) und einem \(x_1\) an. Daher wählen wir aus dem Datensatz die beiden Spalten jump_length und weight. Wir wollen nun feststellen, ob es einen Zusammenhang zwischen der Sprungweite in [cm] und dem Flohgewicht in [mg] gibt. In dem Datensatz finden wir 400 Flöhe wir wählen aber nur zufällig 20 Tiere aus.

model_tbl <- read_csv2("data/flea_dog_cat_length_weight.csv") |>
  select(animal, jump_length, weight) |> 
  filter(animal %in% c("dog", "cat")) |> 
  slice_sample(n = 20)

Neben dem einfachen Datensatz zu den Hunde- und Katzenflöhen haben wir noch Daten zu dem Wachstum von Wasserlinsen. Wir haben einmal händisch die Dichte bestimmt duckweeds_density und einmal mit einem Sensor gemessen. Dabei sind die Einheiten der Sensorwerte erstmal egal, wir wollen aber später eben nur mit einem Sensor messen und dann auf den Wasserlinsengehalt zurückschließen. Wir haben hier eher eine Sätigungskurve vorliegen, denn die Dichte der Wasserlinsen ist ja von der Oberfläche begrenzt. Auch können sich die Wasserlinsen nicht beliebig teilen, es gibt ja nur eine begrenzte Anzahl an Ressourcen.

duckweeds_tbl <- read_excel("data/duckweeds_density.xlsx")

In der Tabelle 55.2 siehst du dann einmal einen Auszug aus den Daten zu den Wasserlinsen. Es ist ein sehr einfacher Datensatz mit nur zwei Spalten. Wie du siehst, scheint sich hierbei um eine nicht lineare Regression zu handeln. Einen linearen Zusammenhang würde ich hier nicht annehmen.

Tabelle 37.1— Auszug aus Wasserlinsendatensatz.

duckweeds_density	sensor
4.8	0.4303
4.8	0.4763
4.8	0.4954
…	…
53.2	2.1187
53.2	2.1296
53.2	2.1246

Auch die Wasserlinsendaten wollen wir uns erstmal in einer Abbildung anschauen und dann sehen, ob wir eine Kurve durch die Punkte gelegt kriegen. Für die Visualisierung der Wasserlinsendaten in der Abbildung 55.10 verzichte ich einmal auf die logarithmische Darstellung. Auffällig ist erstmal, dass wir sehr viel weniger Beobachtungen und auch Dichtemesspunkte auf der \(x\)-Achse haben. Wir haben dann zu den jeweiligen Wasserlinsendichten dann drei Sensormessungen. Das könnte noch etwas herausfordernd bei der Modellierung werden.

ggplot(duckweeds_tbl, aes(duckweeds_density, sensor)) +
  geom_point() +
  theme_minimal() +
  labs(x = "Gemessene Dichte der Wasserlinsen", y = "Sensorwert")

Abbildung 37.1— Visualisierung der Sensorwerte nach Wasserlinsendichte. Pro Dichtewert liegen drei Sensormessungen vor.

37.3 Einfacher Mittelwertsvergleich

Jetzt wollen wir uns den Permutationstest einmal ganz einfach anschauen. Wir wollen einen t-Test für den Mittelwertsvergleich einmal händisch mit der Funktion sample() nachbauen. Dann schauen wir uns die Sachlage einmal in dem R Paket {infer} an.

37.3.1 …mit `sample()`

Wir wollen zuerst einmal mit einem einfachen Mittelwertsvergleich anfangen. Im Prinzip bauen wir hier kompliziert einen t-Test nach. Der t-Test testet, ob es einen signifikanten Mittelwertsunterschied gibt. Anstatt jetzt den t-Test zu rechnen, berechnen wir erstmal das \(\Delta\) und damit den Mittelwertsunterschied der Sprungweiten der Hunde- und Katzenflöhe.

model_tbl |>
  group_by(animal) |> 
  summarise(mean_jump = mean(jump_length)) |> 
  pull(mean_jump) |> 
  diff()

[1] 3.376154

Wir sehen, dass die Hunde- und Katzenflöhe im Mittel einen Unterschied in der Sprungweite von \(3.38cm\) haben. Das ist der Mittelwertsunterschied in unseren beobachteten Daten.

Jetzt wollen wir einmal einen Permutationstest rechnen. Die wichtigste Funktion hierfür ist die Funktion sample(). Die Funktion sample() durchmischt zufällig einen Vektor. Einmal ein Beispiel für die Zahlen 1 bis 10, die wir in die Funktion sample() pipen.

c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10) |> 
  sample()

 [1]  4  5  8  1  3  9  2  7  6 10

Das gleiche Durchmischen findet auch in der Funktion mutate() statt. Wir durchwirblen einmal die Zuordnung der Level des Faktors animal zu den jeweiligen Speungweiten. Dann berechnen wir die Mittelwertsdifferenz für diesen neuen Datensatz. Das machen wir dann \(n\_sim\) gleich 1000 Mal.

n_sim <- 1000

mean_perm_tbl <- map_dfr(1:n_sim, \(x) {
  mean_diff <- model_tbl |>
    ## Permutation Start
    mutate(animal = sample(animal)) |> 
    ## Permutation Ende
    group_by(animal) |> 
    summarise(mean_jump = mean(jump_length)) |> 
    pull(mean_jump) |> 
    diff()
  return(tibble(mean_diff))
})

In der Abbildung 37.2 sehen wir die Verteilung aller Mittelwertsdifferenzen, die aus unserem permutierten Datensätzen herausgekommem sind.

Abbildung 37.2— Histogramm der permutierten Mittelwertsdifferenzen

sum(mean_perm_tbl$mean_diff >= 2.6119)/n_sim

[1] 0.054

Ist das Gleiche als wenn wir dann den Mittelwert berechnen.

mean(mean_perm_tbl$mean_diff >= 2.6119)

[1] 0.054

Teilweise wird diskutiert, ob der \(p\)-Wert hier noch mal 2 genommen werden muss, da wir ja eigentlich zweiseitig Testen wollen, aber da gehen die Meinungen auseinander. Ich belasse es so wie hier.

Dann das ganze nochmal mit einem Student t-Test verglichen und wir sehen, dass wir dem \(p\)-Wert aus einem Student t-Test sehr nahe kommen. Wenn du jetzt noch die Anzahl an Simulationen erhöhen würdest, dann würde sich der \(p_{perm}\) dem \(p_{t-Test}\) immer weiter annähern.

t.test(jump_length ~ animal, data = model_tbl) |> 
  pluck("p.value") |> 
  round(3)

[1] 0.049

Am Ende bleibt dann die Frage, wie viele Permutationen sollen es denn sein? Auch hier sehen wir dann, dass der t-Test signifikant ist, aber der Permutationstest noch nicht. Vielleicht helfen da mehr Permutationen? Oder aber der Effekt ist dann doch zu gering. Hier musst du dann immer überlegen, ob du nicht zu sehr an dem Signifikanzniveau von 5% klebst. Am Ende muss dann der permutierte \(p\)-Wert zudammen mit dem Effekt im Kontext der Fragestellung diskutiert werden.

37.3.2 … mit `{infer}`

Hier wollen wir noch das R Paket {infer} vorstellen.

37.4 Bestimmtheitsmaß \(R^2\)

Das vorherige Beispiel mit dem Mittelwertsvergleich war im Prinzip nur eine Fingerübung für den Ablauf. Wir können auch einfach einen t-Test rechnen und dann ist gut. Anders sieht es aus, wenn wir keinen \(p\)-Wert geliefert bekommen und auch keinen \(H_0\) Testverteilung bekannt ist um einen \(p\)-Wert zu bestimmen. Das ist der Fall bei dem Bestimmtheitsmaß \(R^2\). Wir haben hier keine Teststatistik und dadurch einen \(p\)-Wert vorliegen. Dagegen können wir dann mit einem Permutationstest was tun. Bei dem 95% Konfidenzintervall wird es dann schwieriger, hier müssen wir dann etwas anders permutieren. Wir nutzen dann die Bootstrap Konfidenzintervalle im nächsten Abschnitt.

model_tbl %$%
  lm(jump_length ~ weight) |> 
  glance() |> 
  pull(r.squared)

[1] 0.3000498

Damit haben wir erstmal das Bestimmtheitsmaß aus unseren Daten berechnet. Jetzt stellt sich natürlich die Frage, wie wahrscheinlich ist es den dieses Bestimmtheitsmaß von 0.30 zu beobachten? Dafür lassen wir jetzt einen Permutationstest laufen in dem wir die Daten einmal durchmischen.

n_sim <- 1000

r2_perm_tbl <- map_dfr(1:n_sim, \(x) {
  r2 <- model_tbl |>
    ## Permutation Start
    mutate(weight = sample(weight)) %$%
    ## Permutation Ende    
    lm(jump_length ~ weight) |> 
    glance() |> 
    pull(r.squared)
  return(tibble(r2))
})

In der Abbildung 37.3 sehen wir die Verteilung aller Bestimmtheitsmaße \(R^2\), die aus unserem permutierten Datensätzen herausgekommem sind. Wir erkennen sofort, dass es wenig zufällig bessere Datensätze gibt, die ein höheres Bestimmtheitsmaße \(R^2\) erzeugen.

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_bar()`).

Abbildung 37.3— Histogramm der permutierten Bestimmtheitsmaße \(R^2\)

Jetzt wollen wir einmal bestimmen wie viele Bestimmtheitsmaße \(R^2\) größer sind als unser Bestimmtheitsmaß \(R^2 = 0.29\) aus den Daten.

sum(r2_perm_tbl$r2 >= 0.291164)/n_sim

[1] 0.013

Die Berechnung ist das Gleiche, als wenn wir den Mittelwert aus der logischen Abfrage berechnen würden.

mean(r2_perm_tbl$r2 >= 0.291164)

[1] 0.013

Teilweise wird diskutiert, ob der \(p\)-Wert hier noch mal 2 genommen werden muss, da wir ja eigentlich zweiseitig Testen wollen, aber da gehen die Meinungen auseinander. Ich belasse es so wie hier. Damit haben wir unseren \(p\)-Wert für das Bestimmtheitsmaß \(R^2\) mit \(0.013\). Das ist was wir wollten und somit können wir dann auch sagen, dass wir einen signifikantes Bestimmtheitsmaß \(R^2\) vorliegen haben. Was noch fehlt ist das 95% Konfidenzintervall, was wir Mithilfe des Bootstrapverfahrens berechnen wollen.

37.5 Bootstrap

Bootstrap oder Bootstraping ist eine mächtige Methode um Vertrauensbereiche zu bestimmen. Wir wollen mehr oder minder wissen, in welchen Bereich unsere Beobachtungen fallen könnten, wenn wir unser Experiment wiederholen. Wir haben ja meistens nur einen sehr kleinen Datensatz, der nur einen Ausschnitt repräsentiert. Durch Bootstrap können wir jetzt auch ohne eine geschlossene mathematische Formel für die Konfidenzintervalle oder Prädiktionsintervalle die jeweiligen Intervalle berechnen. Dennoch musst du beachten, dass wir ein paar Beobachtungen brauchen. Wenn du nur fünf Beobachtungen pro Gruppe hast, wird es meistens sehr schnell sehr eng mit dem Bootstrapverfahren.

37.5.1 95% Konfidenzintervalle

Wir können die Methode des Bootstraping nutzen um uns die 95% Konfidenzintervalle über eine Simulation bzw. Permutation berechnen zu lassen. Haben wir in dem Permutatiosntest noch die Spalten permutiert so permutieren wir bei dem Bootstrap-Verfahren die Zeilen. Da wir aber keinen neuen Datensatz erhalten würden, wenn wir nur die Zeilen permutieren, ziehen wir einen kleineren Datensatz mit zurücklegen. Das heißt, dass wir auch Beobachtungen mehrfach in unseren gezogenen Bootstrapdatensatz haben können. Wir nutzen in R die Funktion slice_sample() in der wir dann auswählen, dass 80% der Beobachtungen mit zurücklegen gezogen werden sollen. Das Zurücklegen können wir mit der Option replace = TRUE einstellen. Wir führen das Bootstraping dann insgesamt 1000 mal durch.

n_boot <- 1000

r2_boot_tbl <- map_dfr(1:n_boot, \(x) {
  r2 <- model_tbl |>
    ## Bootstrap Start
    slice_sample(prop = 0.8, replace = TRUE) %$%
    ## Bootstrap Ende
    lm(jump_length ~ weight) |> 
    glance() |> 
    pull(r.squared)
  return(tibble(r2))
})

Nachdem wir nun unser Bootstrap gerechnet haben und eintausend Bestimmtheitsmaße bestimmt haben, können wir einfach das \(2.5\%\) und \(97.5\%\) Quantile bestimmen um unser 95% Konfidenzintervall zu bestimmen. Zwischen \(2.5\%\) und \(97.5\%\) liegen ja auch 95% der Werte der eintausend Bestimmtheitsmaße.

r2_boot_tbl %$% 
  quantile(r2, probs = c(0.025, 0.975)) |> 
  round(3)

 2.5% 97.5% 
0.016 0.620

Wir hatten ein Bestimmtheitsmaß \(R^2\) von \(0.30\) berechnet und können dann die untere und obere 95% Konfidenzgrenze ergänzen. Wir erhaltend dann \(0.300\, [0.016; 0.620]\), somit liegt unser beobachtetes Bestimmtheitsmaß \(R^2\) mit 95% Sicherheit zwischen \(0.016\) und \(0.620\).