stat-tests-permutationstest.qmd

```{r echo = FALSE}
pacman::p_load(tidyverse, readxl, knitr, kableExtra)
set.seed(20230526)
```

# Permutationstest & Bootstraping {#sec-permut}

*Letzte Änderung am `r format(fs::file_info("stat-tests-permutationstest.qmd")$modification_time, '%d. %B %Y um %H:%M:%S')`*

> *"Wer zu spät kommt, den bestraft das Leben." --- Gorbatschow zu Honecker beim Staatsbesuch in der DDR am 7. Oktober 1989*

------------------------------------------------------------------------

![](images/caution.png){fig-align="center" width="50%"}

**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 @fay2018biologist zu sagen:

> *"\[...\], if statistics were being invented today, bootstrapping, as well as related re-sampling methods, would be the standard go-to approach." --- @fay2018biologist*

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.

::: callout-tip
## 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](https://mac-theobio.github.io/QMEE/lectures/permutation_examples.notes.html). Eine der ersten Anlaufstellen, wenn man sich mit dem Thema Permutatuonstests in R beschäftigen will.
-   Das [R Paket `{infer}`](https://infer.netlify.app/) 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](http://www.wormbook.org/chapters/www_statisticalanalysis/statisticalanalysis.html#sec6-7) aus @fay2018biologist 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}`](https://rsample.tidymodels.org/index.html) 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.
:::

## Genutzte R Pakete

Wir wollen folgende R Pakete in diesem Kapitel nutzen.

```{r echo = TRUE}
#| message: false
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.

## 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.

```{r}
#| message: false

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.

```{r}
duckweeds_tbl <- read_excel("data/duckweeds_density.xlsx")
```

In der @tbl-duckweeds 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.

```{r}
#| echo: false
#| message: false
#| warning: false
#| label: tbl-duckweeds
#| tbl-cap: "Auszug aus Wasserlinsendatensatz."

rbind(head(duckweeds_tbl, n = 3),
      rep("...", times = ncol(duckweeds_tbl)),
      tail(duckweeds_tbl, n = 3)) |> 
  kable(align = "c", "pipe")
```

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 @fig-gam-duckweeds-01 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.

```{r}
#| echo: true
#| warning: false
#| message: false
#| label: fig-gam-duckweeds-01
#| fig-align: center
#| fig-height: 5
#| fig-width: 6
#| fig-cap: "Visualisierung der Sensorwerte nach Wasserlinsendichte. Pro Dichtewert liegen drei Sensormessungen vor."

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

## 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.

### ...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.

```{r}
model_tbl |>
  group_by(animal) |> 
  summarise(mean_jump = mean(jump_length)) |> 
  pull(mean_jump) |> 
  diff()
```

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.

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

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.

```{r}
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 @fig-permut-01 sehen wir die Verteilung aller Mittelwertsdifferenzen, die aus unserem permutierten Datensätzen herausgekommem sind.

```{r}
#| echo: false
#| message: false
#| label: fig-permut-01
#| fig-align: center
#| fig-height: 5
#| fig-width: 5
#| fig-cap: "Histogramm der permutierten Mittelwertsdifferenzen"

ggplot(mean_perm_tbl , aes(mean_diff)) +
  theme_minimal() +
  geom_histogram() +
  geom_vline(xintercept = 2.6119, color = "red") +
  labs(x = "Permutierte Mittelwertsdifferenz",
       y = "Anzahl")
```

```{r}
sum(mean_perm_tbl$mean_diff >= 2.6119)/n_sim
```

Ist das Gleiche als wenn wir dann den Mittelwert berechnen.

```{r}
mean(mean_perm_tbl$mean_diff >= 2.6119) 
```

::: column-margin
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.

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

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.

### ... mit `{infer}`

Hier wollen wir noch das [R Paket `{infer}`](https://infer.netlify.app/) vorstellen.

## 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.

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

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.

```{r}
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 @fig-permut-02 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.

```{r}
#| echo: false
#| message: false
#| label: fig-permut-02
#| fig-align: center
#| fig-height: 5
#| fig-width: 5
#| fig-cap: "Histogramm der permutierten Bestimmtheitsmaße $R^2$"

ggplot(r2_perm_tbl, aes(r2)) +
  theme_minimal() +
  geom_histogram(binwidth = 0.01) +
  xlim(0, 1) +
  ylim(0, 200) +
  geom_vline(xintercept = 0.291164, color = "red") +
  labs(x = expression(Permutiertes~Bestimmtheitsmaß~R^2),
       y = "Anzahl")
```

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

```{r}
sum(r2_perm_tbl$r2 >= 0.291164)/n_sim 
```

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

```{r}
mean(r2_perm_tbl$r2 >= 0.291164) 
```

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.

## 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.

### 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.

```{r}
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.

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

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$.

### Prädiktionsintervall

Idee vom Botstraping und Resample [Working with resampling sets](https://rsample.tidymodels.org/articles/Working_with_rsets.html)

[Bootstrap confidence intervals für nicht-lineare Daten](https://rsample.tidymodels.org/articles/Applications/Intervals.html) für unseren Datensatz @tbl-duckweeds

[Das R Paket `{resample}`](https://rsample.tidymodels.org/index.html)

[R Paket `{workboots}`](https://markjrieke.github.io/workboots/)

[Understanding Prediction Intervals](https://www.bryanshalloway.com/2021/03/18/intuition-on-uncertainty-of-predictions-introduction-to-prediction-intervals/)

[Simulating Prediction Intervals](https://www.bryanshalloway.com/2021/04/05/simulating-prediction-intervals/)

[Quantile Regression Forests for Prediction Intervals](https://www.bryanshalloway.com/2021/04/21/quantile-regression-forests-for-prediction-intervals/)

[P-values for prediction intervals](https://robjhyndman.com/hyndsight/forecasting-pvalues.html) machen keinen Sinn

[R Paket `{spin}`](https://github.com/brshallo/spin)

## Referenzen {.unnumbered}