# Ein `{purrr}` Cookbook {#sec-purrr-furrr}
*Letzte Änderung am `r format(fs::file_info("programing-purrr-furrr.qmd")$modification_time, '%d. %B %Y um %H:%M:%S')`*
> *"Until I realized, it is the struggle itself that is most important. We must strive to be more than we are, Lal. It does not matter that we will never reach our ultimate goal. The effort yields its own rewards." --- Lieutenant Commander Data, Star Trek: The Next Generation, The Offspring*
{{< video https://youtu.be/gcKDbf4cL7M >}}
::: {.callout-caution collapse="true"}
## Was soll das hier? Was für ein Kochbuch?
{#fig-math-cartoon fig-align="center" width="100%"}
Dieses Kapitel dient dazu *fortgeschrittene* Programmierung in R zu präsentieren. Es ist eine Sammlung von Code und Rezepten, die ich immer mal wieder nutze. Deshalb hat dieses Kapitel auch den Titel *Cookbook*. Wenn dich also Programmieren interessiert, dann kannst du dir hier noch was anschauen. Das Kapitel ist eigentlich nie fertig, da ich sicherlich immer mal wieder was ergänzen werde. Ich habe mich übrigens dann auch für eine bessere Schriftart im RStudio entschieden. Ich nutze jetzt [Fira Code: free monospaced font with programming ligatures](https://github.com/tonsky/FiraCode) und es ist großartig.
:::
In diesem Kapitel geht es hauptsächlich um [Iterationen](https://r4ds.had.co.nz/iteration.html#iteration). Das heißt wir wollen immer das Gleiche auf verschiedene unterschiedliche Dinge anwenden. In unserem Fall ist das "Gleiche" eine Funktion `function(...){}` in R und das "Unterschiedliche" sind Einträge in einer Liste `lst()` oder einem Datensatz `tibble()`. Wir wollen also zum Beispiel auf verschiedenen Datensätzen mit fünf unterschiedlichen Outcomes immer wieder eine multiple lineare Regression rechnen. Anstatt also per Copy&Paste fünfmal den Code zu kopieren, wollen wir alle Datensätze in einer Liste speichern und die Liste dann in einem Schwung auswerten. Das geht natürlich nur eingeschränkt. Wir müssen ja noch beachten, dass jedes Outcome einer anderen Verteilung folgen könnte.
Wie wir die Daten in Gruppen zusammenfassen ist unterschiedlich möglich. Ich präsentiere hier zum einen die Funktion `split()`, die Daten in Listen aufteilt sowie die Funktion `nest()` die Daten in ein `tibble` zusammenfaltet. Beide Varianten haben so ihre Vor- und Nachteile. Wie immer kommt es dann auf die Anwendung an und was du machen willst.
Es lohnt sich hierbei die Listen bzw. die Datensätze, über die dann iteriert werden soll, so gleich wie irgendwie möglich zu bauen. Das heißt, dass wir die gleiche Anzahl an Spalten mit den gleichen Spaltennamen vorliegen haben wollen. Das erleichtert dann die spätere Anwendung mit `map()`.
## Genutzte R Pakete
Wir wollen folgende R Pakete in diesem Kapitel nutzen.
```{r echo = TRUE}
#| message: false
pacman::p_load(tidyverse, magrittr, rstatix,
janitor, purrr, furrr, see,
readxl, tictoc, multcompView,
parameters, scales,
conflicted)
conflicts_prefer(magrittr::extract)
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.
## Die Daten
Als Datenbeispiel schauen wir uns einmal ein einfaktoriellen Datensatz an und suchen uns auch nur acht Zeilen und drei Outcomes raus. Wir könnten die Analyse auch über den vollen Datensatz rechnen, aber dann wird hier alles sehr voll. Es geht ja hier mehr um die Demonstration.
```{r}
#| message: false
#| warning: false
soil_tbl <- read_excel("data/soil_1fac_data.xlsx") |>
mutate(variante = str_c(variante, "_", amount),
variante = as_factor(variante),
across(where(is.numeric), round, 2)) |>
select(-amount) |>
extract(1:8, 1:4) |>
pivot_longer(cols = fe:no3,
names_to = "outcome",
values_to = "rsp")
```
Als zweiten Datensatz nehmen wir noch eine zweifaktorilles Design mit einem Behandlungs- und einem Blockeffekt. Darüberhinaus haben wir dann noch verschiedene Outcomes und diese Outcomes dann auch an zwei Orten, einmal im Blatt und einmal im Stiel, gemessen. Das heißt, wir haben immer eine Outcome/Sample-Kombination vorliegen. Bei drei Outcomes und zwei Messorten macht das dann sechs Kombinationen auf denen wir dann immer unsere zweifaktoriellen Analysen rechnen wollen.
```{r}
spinach_tbl <- read_excel("data/spinach_metal_data.xlsx") |>
mutate(trt = as_factor(trt),
sample = as_factor(sample),
block = as_factor(block)) |>
pivot_longer(cols = fe:zn,
names_to = "outcome",
values_to = "rsp") |>
mutate(outcome = as_factor(outcome))
```
## Daten aufteilen...
In R haben wir zwei Möglichkeiten für `map()` die Daten aufzuteilen. Klar, wir können die Aufteilung sicherlich auch mit anderen Funktionen machen, aber diese beiden Funktionen sind sehr nützlich.
### ... mit `split()`
Mit der Funktion `split()` können wir einen Datensatz nach *einer* Faktorspalte in eine Liste aufspalten. Die Liste ist dann auch gleich so benannt wie das Level es Faktors für das wir die Aufteilung gemacht haben. Die Benamung der Liste ist dann praktisch, wenn wir später wieder einen Datensatz aus den Ergebnissen bauen.
```{r}
soil_lst <- soil_tbl |>
split(~outcome)
soil_lst
```
### ... mit `nest()`
Wenn wir mehr als eine Gruppierungsspalte haben, dann können wir die Funktion `nest()` nutzen. In unserem Beispiel haben wir die Spalte `outcome` und `sample`. Für jede Kombination der beiden Spalte wollen wir dann jeweils ein Modell rechnen. Hier meine ich mit Modell eine lineare Regression und dann eine ANOVA. Als erstes müssen wir unsere Daten gruppieren und dann können wir die Daten nesten. Mit `unnest()` lässt sich dann die genestete Struktur wieder in einen normalen Datensatz zurückführen.
```{r}
spinach_nest_tbl <- spinach_tbl |>
group_by(sample, outcome) |>
nest()
spinach_nest_tbl
```
## Mit `{purrr}` über Daten {#sec-purrr}
Das [R Paket `{purrr}`](https://purrr.tidyverse.org/) erlaubt es sehr effizient immer das Gleiche auf Listeneinträgen oder genereller auf Daten anzuwenden. Wir können uns dabei selber eine Funktion schreiben oder aber schon implementierte Funktionen anwenden. Gehen wir einmal alle Funktionen durch. Wir werden hier nicht alle zeigen, aber es ist gut einmal zu wissen, welche Funktionen es gibt. Schaue auch mal in das Cheat Sheet des R Paketes `{purrr}` rein: [Apply functions with purrr::cheat sheet](https://evoldyn.gitlab.io/evomics-2018/ref-sheets/R_purrr.pdf).
### Mit `map()` auf Listeneinträgen
Wir können die Funktion `map()` und deren verwandete Funktionen nutzen um zügig über Listen andere Funktionen anzuwenden.
- `map()` erlaubt über eine Liste von Datensätzen ein Funktion anzuwenden. Dabei können wir dann die einzelnen Listeneinträge über `.x` an die Funktionen weitergeben. Siehe hierzu auch [Basic map functions](https://dcl-prog.stanford.edu/purrr-basics.html).
- `map2()` erlaubt es über zwei gleichlange Vektoren zu laufen. Wir können hier zwei Optionen in der Form `.x`, `.y` an die Funktion weitergeben. Siehe hierzu auch [Map with multiple inputs](https://dcl-prog.stanford.edu/purrr-parallel.html).
- `pmap()` kann nun über eine Liste an Vektoren laufen und somit mehrere Inputoptionen verarbeiten. Damit ist `pmap()` die Generalisierung der `map()` Funktion. Siehe hierzu auch [Map with multiple inputs](https://dcl-prog.stanford.edu/purrr-parallel.html).
- `walk()` ist ein *silent* `map()`. Damit können wir Daten in eine Datei schreiben, ohne ein Output wieder zubekommen.
- `imap()` können wir nutzen, wenn wir den Index $i$ wieder haben wollen. Das heißt, wir wollen über einen Vektor laufen und brauchen dafür den Index. Hier hilft die `imap()` Familie.
- `modify()`können wir anwenden, wenn wir nur Spalten modifizieren oder mutieren wollen. Wir haben einen Datensatz und wollen alle `character` Spalten in einen Faktor umwandeln. Siehe hierzu auch [Modify elements selectively](https://purrr.tidyverse.org/reference/modify.html).
Schauen wir uns die Anwendung von der Funktion `map()` auf eine Liste an. In folgenden Code entfernen wir einmal in *jedem* Listeneintrag die Spalte `outcome`. Dann lassen wir uns von *jedem* Listeneintrag die erste Zeile wiedergeben.
```{r}
soil_lst |>
map(select, -outcome) |>
map(head, 1)
```
Wir können zum einen `map(head, 1)` schreiben oder aber den Listeneintrag als `.x` direkt in die Funktion weiterleiten. Dann schreiben wir `map(~head(.x, 1))` und müssen noch die Tilde `~` vor die Funktion setzen.
```{r}
soil_lst |>
map(~head(.x, 1))
```
Die richtige Stärke entwickelt dann `map()`, wenn wir mehrere Funktionen hineinander schalten. In unserem Fall rechnen wir einen Games-Howell Test und wollen uns dann das *compact letter display* wiedergeben lassen. Da wir am Ende dann einen Datensatz haben wollen, nutzen wir die Funktion `bind_rows()` um die Listeneinträge in einen Datensatz zusammen zukleben.
```{r}
soil_lst |>
map(~games_howell_test(rsp ~ variante, data = .x)) |>
map(~mutate(.x, contrast = str_c(.x$group1, "-", .x$group2))) |>
map(pull, p.adj, contrast) |>
map(~multcompLetters(.x)$Letters) |>
bind_rows(.id = "outcome")
```
Wir können auch auch auf genesteten Daten die Funktion `map()` anwenden. In diesem Fall belieben wir die ganze Zeit in einem `tibble`. Wir lagern unsere Ergebnisse sozusagen immer in einer Zelle und können auf diese Einträge dann immer wieder zugreifen. Einmal zu Demonstration rechnen wir sechs Mal ein lineares Modell mit der Funktion `lm()` und speichern das Ergebnis in der Spalte `model`. Wir haben jetzt dort jeweils `<lm>` stehen. Damit wissen wir auch, dass wir dort unser Modell drin haben. Wir können jetzt auf der Spalte `model` weiterechnen und uns neue Spalten mit Ergebnissen erschaffen.
```{r}
spinach_nest_tbl |>
mutate(model = map(data, ~lm(rsp ~ trt + block, data = .x)))
```
Im Folgenden rechnen wir ein lineares Modell, dann eine ANOVA und lassen uns das Ergebnis der ANOVA mit der Funktion `model_parameters()` aufhübschen. Wie du siehst, geben wie immer den Spaltennamen eine Funktion weiter. Dann wählen wir noch die Spalten, die wir dann unnesten wollen.
```{r}
spinach_nest_tbl <- spinach_nest_tbl |>
mutate(model = map(data, ~lm(rsp ~ trt + block, data = .x))) |>
mutate(anova = map(model, anova)) |>
mutate(parameter = map(anova, model_parameters)) |>
select(sample, outcome, parameter)
```
Zum Abschluss nutzen wir die Funktion `unnest()` um uns die aufgehübschten Ergebnisse der ANOVA wiedergeben zu lassen. Dann will ich noch, dass die Namen alle klein geschrieben sind und auch sonst sauber sind. Dafür nutze ich dann die Funktion `clean_names()`. Abschließend filtere ich und runde ich noch die Ergebnisse. Am Ende will ich dann nur die Kombinationen aus `sample` und `outcome` haben sowie den $p$-Wert aus der ANOVA.
```{r}
spinach_nest_tbl |>
unnest(parameter) |>
clean_names() |>
mutate(across(where(is.numeric), round, 2)) |>
filter(parameter != "Residuals") |>
select(sample, outcome, parameter, p)
```
Hier noch ein weiteres Beispiel für `split()`, `nest()` und `nest_by()` zum Ausprobieren und rumspielen. Wir wollen für hier einmal auf ganze vielen Behandlungen den Shapiro-Wilk-Tests für die Abweichung von der Normalverteilung rechnen. Dazu laden wir uns einmal die Daten `clove_germ_rate.xlsx`.
```{r}
clove_tbl <- read_excel("data/clove_germ_rate.xlsx") |>
mutate(clove_strain = as_factor(clove_strain),
germ_rate = as.numeric(germ_rate))
```
Wir rechnen einmal die Shapiro-Wilk-Tests über die Funktion `split()` und dann einer Liste.
```{r}
clove_tbl |>
split(~clove_strain) |>
map(~shapiro.test(.x$germ_rate)) |>
map(tidy) |>
bind_rows(.id = "test") |>
select(test, p.value) |>
mutate(decision = ifelse(p.value <= 0.05, "reject normal", "normal"),
p.value = pvalue(p.value, accuracy = 0.001))
```
Dann das selbe nochmal mit der Funktion `nest_by()`, die jetzt Vektoren generiert.
```{r}
clove_tbl |>
nest_by(clove_strain) |>
mutate(shapiro = map(data, ~shapiro.test(.x)),
clean = tidy(shapiro)) |>
reframe(clean)
```
Und nochmal mit der Pipe von `group_by()` zu `nest()`. Die Funktion `nest()` hate auch eine `.by =`-Option, so dass wir auch den Schritt mit `group_by()` weglassen könnten.
```{r}
clove_tbl |>
group_by(clove_strain) |>
nest() |>
mutate(shapiro = map(data, ~shapiro.test(.x$germ_rate)),
clean = map(shapiro, tidy)) |>
unnest(clean)
```
::: callout-tip
## Besseres `unnest()`
Es hilft hier auch sich einmal die Funktion [`unnest_wider()`](https://tidyr.tidyverse.org/reference/unnest_wider.html) sowie [`unnest_longer()`](https://tidyr.tidyverse.org/reference/unnest_longer.html) anzuschauen. Beide Funktionen liefern manchmal gleich das gewünschte Resultat. Besonders hilfreich sind beide Funktionen, wenn es um Listen geht. Aber da musst du dann einmal konkret schauen, was du brauchst.
:::
### Weitere Funktionen wie `keep`, `reduce` und Co.
Im Folgenden nochmal eine Sammlung von weiteren Funktionen, die ich im Rahmen meiner Data Science Analysen dann immer mal wieder nutze und gerne vergesse. Da aber die Funktionen so praktisch sind, habe ich mir mal alles hier mit aufgeschrieben. Tja, teilweise ist das Skript hier auch eine große Fundgrube für mich.
::: column-margin
Das [Data Wrangling with dplyr and tidyr - Cheat Sheet](https://www.rstudio.com/wp-content/uploads/2015/02/data-wrangling-cheatsheet.pdf) ist auch immer eine gute Hilfe. Auch ist das [Apply functions with purrr Cheat - Sheet](https://evoldyn.gitlab.io/evomics-2018/ref-sheets/R_purrr.pdf) ist auch immer ein Blick wert.
:::
Wenn ich einen Vektor habe mit Namen, der teilweise aus `map()` herauskommt, dann kann ich den Vektor über `enframe()` und der [Dokumentation zu `enframe()`](https://tibble.tidyverse.org/reference/enframe.html) in ein `tibble` umwandeln.
```{r}
named_vec <- c(dog = 2, cat = 4, dog = 3, fox = 5)
enframe(named_vec)
```
Das ist manchmal extrem praktisch, wenn eine Liste nur einen Eintrag hat. Dann können wir über `unlist()` einen benamten Vektor bauen und den Vektor dann in einen `tibble` umwandeln.
```{r}
lst(id237189 = 4,
id7629w0 = 5,
id790182 = 3) |>
unlist() |>
enframe()
```
Ich habe häufig Listen und möchte nach gewissen Merkmalen die Listeneinträge filtern - manchmal auch Datensätze. Hier hilft die Funktion `keep()` und der [Dokumentation zu `keep()`](https://purrr.tidyverse.org/reference/keep.html) mit der ich dann Listeneinträge nach einer Funktion, die `TRUE` oder `FALSE` ausgibt filtern kann. Das Gange geht dann mit `discard()` auch in der Verneinung von `keep()`. Ich will in der Folge nur die Listeneinträge behalten in denen das entsprechende `tibble` im Listeneintrag mehr als zwei Zeilen hat.
```{r}
data_lst <- lst(dog_tbl = tibble(id = c("id786", "id987", "id231", "id566"),
dog = rep("dog", 4)),
cat_tbl = tibble(id = c("id786", "id566"),
cat = rep("cat", 2)),
fox_tbl = tibble(id = c("id786", "id987", "id776", "id129", "id231", "id566"),
fox = rep("fox", 6)))
data_lst |>
keep(\(x) nrow(x) > 2)
```
Wenn ich in einer Liste Einträge finde, die keinen Werte haben, also `NULL` sind, dann kann ich über die Funktion `compact()` die leeren Listeneinträge einfach entfernen.
```{r}
data_null_lst <- lst(dog_tbl = tibble(animal = rep("dog", 4)),
cat_tbl = NULL,
fox_tbl = NULL)
data_null_lst |>
compact()
```
Die Funktion `reduce` und der [Dokumentation zu `reduce()`](https://purrr.tidyverse.org/reference/reduce.html) ermöglicht es mit Listeneinträge zu kombinieren. Wenn ich alle Listeneinträge nur untereinander packen will, dann nutze ich die Funktion `bind_rows()` oder gleich `map_dfr()`.
```{r}
data_lst |>
reduce(left_join, by = "id")
```
Wenn alle Listeneinträge nach einer Spalte zusammengführt werden sollen, dann Nutze ich die Funktion `left_join()`. Weiteres zum Zusammenführen (eng. *merge*) von Datensätzen in der [Dokumentation zu `join()`](https://dplyr.tidyverse.org/reference/mutate-joins.html). Beim *Mergen* von Datensätzen muss man recht viel Nachdenken und überlegen was man eigentlich will, deshalb kann ich hier keine vollständige Abhandlung liefern.
```{r}
sort_vec <- c("cat", "fox", "dog")
data_tbl <- tibble(animal = c("dog", "cat", "fox"),
jump_length = c(4.1, 5.8, 6.2))
data_tbl
```
Dann möchte ich meist noch einen Datensatz nach einer externen Spalte sortieren. Dafür nutze ich die Funktion `arrange()` wie folgt.
```{r}
data_tbl |>
arrange(factor(animal, levels = sort_vec))
```
Manchmal möchte man dann auch einem Long-Format wieder ein Wide-Format machen. Zum einen natürlich wieder der Verweis auf die Funktion `pivot_wider()` und der [Dokumentation zu `pivot_wider()`](https://tidyr.tidyverse.org/reference/pivot_wider.html). Wir nutzen hier aber auch das R Paket `{glue}` und der [Dokumentation zu `{glue}`](https://glue.tidyverse.org/) um uns die Spaltennamen zu bauen. In unserem Beispiel kriegt jede Kuh zwei Zeilen mit Werten. Wir wollen aber den die Spalte `type` auflösen und nur noch eine Zeile pro Kuh vorliegen haben.
```{r}
cow_tbl <- tibble(lom = c("276000355496007", "276000355496007"),
type = c("val", "lal"),
rz = c(128, 254),
lfd_nr = c(3, 4))
cow_tbl
```
Wir nehmen also die Namen aus der Spalte `type` und kleben (eng. *glue*) erst den Namen der Spalte und dann mit einem Unterstrich getrennt den Namen der beiden anderen Spalten `rz` und `lfd_nr` dran. Wir wollen ja auch die Werte von den beiden Spalten dann in das Wide-Format übertragen.
```{r}
cow_tbl |>
pivot_wider(names_from = type,
names_glue = "{type}_{.value}",
values_from = c(rz, lfd_nr))
```
## Mit `{furrr}` parallel über Daten {#sec-furrr}
Warum geht es den jetzt hier? Wenn du `purrr` und die Funktionen `map()` verstanden hast, dann geht natürlich alles auch in paralleler Berechnung. Die parallele Berechnung ist in dem [R Paket `{furrr}`](https://furrr.tidyverse.org/) implementiert. Das heißt wir müssen nur die Funktionsnamen ändern und schon rechnet sich alles in Parallel. Wir nutzen also nicht nur einen Kern von unseren Rechnern sondern eben alles was wir haben.
```{r}
no_cores <- availableCores() - 1
no_cores
```
Einmal das ganze in sequenzieller Programmierung. Also alles nacheinander gerechnet.
```{r}
plan(sequential)
tic()
nothingness <- future_map(c(2, 2, 2), ~Sys.sleep(.x))
toc()
```
Der folgende Code sollte ca. 2 Sekunden dauern, wenn der Code parallel läuft. Wir haben einen kleinen Overhead in `future_map()` durch das Senden von Daten an die einzelnen Kerne. Es gibt auch einmalige Zeitkosten für `plan(multisession)`, um die Kerne einzurichten.
```{r parallel-session}
#| cache: true
plan(multisession, workers = 3)
tic()
nothingness <- future_map(c(2, 2, 2), ~Sys.sleep(.x))
toc()
```
Wie du siehst, musst du nur `future_` vor die `map()` Funktion ergänzen und schon kannst du parallel rechnen.
## `{{}}`-Operator (curly-curly)
Auch hier nochmal für mich eine Erinnerung was der `{{}}`-Operator, ausgesprochen *curly-curly*-Operator, ist. Damit können dann Funktionen in R besser gebaut werden. Mehr dazu dann unter [The good and bad of tidy evaluation with `{rlang}`](https://www.tidyverse.org/blog/2019/06/rlang-0-4-0/) sowie der Hilfeseite mit [Programming with dplyr](https://dplyr.tidyverse.org/articles/programming.html) für eine umfangreiche Betrachtung.
## progressr: An Introduction
Die Funktion `map()` hat die Option `.progress = TRUE` mit der du dir auch einen Fortschritt anzeigen lassen kannst. Also wie lange noch die Funktion braucht um über alle Listeneinträge zu rechnen. Wenn du es noch schöner haben willst, dann schaue dir einmal das R Paket [`{progressr}`: An Introduction](https://cran.r-project.org/web/packages/progressr/vignettes/progressr-intro.html) an.