Elegante R-Programmierung mit purrr::map und genisteten Datensätzen

Hex-Sticker des R-Pakets purrr

2016 machte Hadley Wickham eine Idee populär, von der er zunächst selbst nicht sicher war, ob sie gut ist: genistete Datensätze (nested data frames). Das Prinzip ist einfach: Eine Spalte eines Datensatzes kann selbst ein Datensatz sein. Was zunächst umständlich oder verwirrend klingt, kann zum mächtigen Werkzeug werden – vor allem, wenn man viele gleich aufgebaute Modelle für Untergruppen erstellen möchte.

Hier ein Beitrag von Hadley Wickham dazu von der Plotcon 2016.

tidyr, ein mächtiges Werkzeug zur Datenanalyse

Das R-Paket tidyr war mir hauptsächlich als Werkzeug zur Datenaufbereitung für Diagramme mit ggplot2 vertraut. Die Funktion gather() bildet Schlüssel-Wert-Paare, um beispielsweise aus Variablen wie Kategorie 1, Kategorie 2, Kategorie 3 (unterschiedliche Spalten im Datensatz) eine Variable Kategorie sowie eine Variable Wert zu bilden, sodass anschließend Farben, Formen, Linientypen, Legenden etc. auf die Kategorie Bezug nehmen können. Gegenstück zu gather(): spread(). Neuere Entwicklungen: pivot_longer() und pivot_wider, siehe hier.

In unserem Zusammenhang geht es um die Funktion nest(), die sogenannte list columns anlegt: Variablen, die selbst einen Datensatz enthalten.

Die Daten: Gapminder – in Erinnerung an Hans Rosling

Ich kann der Versuchung nicht widerstehen, den gleichen Datensatz zu verwenden wie Hadley Wickham in seinem Beispiel: Die Gapminder-Daten, die Hans Rosling populär gemacht hat. Wer das folgende Video noch nicht kennt, dem seien die weniger als fünf Minuten wärmsten empfohlen. Hier demonstriert Rosling, wie man in so kurzer Zeit 120.000 Datenpunkte präsentieren und dabei eine spannende Geschichte erzählen kann. Leider ist dieser inspirierende Mann bereits verstorben.

Ziel: Lineare Regressionsmodelle für 142 Länder

Die gapminder-Daten aus dem gleichnamigen R-Paket enthalten Lebenserwartung (lifeExp), Jahre (year, in Fünferschritten von 1952 – 2007), Bevölkerung (pop) sowie den Wohlstand (gdpPercap in US-Dollar, inflationsbereinigt / BIP pro Kopf) für 142 Länder (country) aus fünf Kontinenten.

Forschungsfrage: Wie verändert sich die Lebenserwartung über die Zeit je Land?

Methode: Lineare Regressionsmodelle für alle 142 Länder, um anhand einer einfachen Kennzahl (R-Quadrat) zu vergleichen, wie gut sich die Lebenserwartung nur mit dem Erhebungsjahr vorhersagen bzw. erklären lässt.

Erste Idee: Einzelnes Regressionsmodell aufstellen

library(tidyverse)
library(gapminder)
data(gapminder)
str(gapminder)

d <- gapminder %>% 
   filter(country == "Germany")
d_model <- lm(lifeExp ~ year, data = d)
d_model_summary <- summary(d_model)
d_model_summary$r.squared

Hier habe ich das Land Deutschland ausgewählt, ein Regressionsmodell aufgestellt (d_model), das summary aufgerufen und daraus das R² extrahiert. Ohne elegantere Technik müsste ich das nun für alle Länder wiederholen. Das widerspricht dem Programmierprinzip DRY: Don’t Repeat Yourself – wiederhole Dich nicht. Wie geht es also eleganter? Mit genisteten Datensätzen.

Genistete Datensätze am Beispiel gapminder

Wie sieht das aus?

by_country <- gapminder %>% 
   group_by(country) %>% 
   nest()

head(by_country)

# A tibble: 6 x 2
 country       data             
 <fct>         <list>           
 1 Afghanistan <tibble [12 x 5]>
 2 Albania     <tibble [12 x 5]> 
 3 Algeria     <tibble [12 x 5]> 
 4 Angola      <tibble [12 x 5]> 
 5 Argentina   <tibble [12 x 5]> 
 6 Australia   <tibble [12 x 5]> 

Im Originaldatensatz umfasst jedes Land 12 Zeilen: je eine pro Jahr, von 1952 bis 2007 in Fünferschritten. Der genistete Datensatz enthält nur noch eine Zeile pro Land, mit den Daten (data) als Datensatz. Das R-Objekt heißt tibble – eine moderne Erweiterung des data frame.

Regressionsmodelle mit funktionaler Programmierung erstellen

R-Programmierung zur Datenanalyse ist im Wesentlichen funktionales Programmieren – getreu dem Motto von John Chambers, Entwickler von S (worauf R aufbaut) und Mitglied der R Core Group: „Everything that happens is a function call“ – Alles, was passiert, ist ein Funktionsaufruf.

Somit schreiben wir uns eine Funktion, die das R-Quadrat der Regression lifeExp ~ year zurückliefert:

get_rsq <- function(data) {
   model <- lm(lifeExp ~ year, data = data)
   model_summary <- summary(model)
   return(model_summary$r.squared)
 }

Das Erstellen von 142 Regressionsmodellen für die 142 Länder ist nun mit einem kurzen R-Aufruf möglich. Man könnte eine for- oder while-Schleife programmieren, die die Länder abarbeitet. Eleganter ist es, das Paket purrr (mit drei „r“ / englisch für Katzenschnurren) zu verwenden:

RQuadrate <- by_country %>% 
   mutate(rsq = map_dbl(data, get_rsq))

head(RQuadrate)

# A tibble: 6 x 3  
  country     data              rsq  
  <fct>       <list>            <dbl>  
1 Afghanistan <tibble [12 x 5]> 0.948  
2 Albania     <tibble [12 x 5]> 0.911  
3 Algeria     <tibble [12 x 5]> 0.985  
4 Angola      <tibble [12 x 5]> 0.888  
5 Argentina   <tibble [12 x 5]> 0.996  
6 Australia   <tibble [12 x 5]> 0.980   

Die map-Funktionen im purrr-Paket funktionieren ähnlich wie die apply-Funktionen in Base R, sind aber konsistenter in der Syntax und vor allem datentyp-stabil. Hier verwenden wir map_dbl, um sicherzustellen, dass wir numerische Werte (mit Dezimalstellen) zurückerhalten.

Ganze Modelle im Datensatz speichern statt nur R²

Alternativ kann man auch das komplette Modellergebnis im Datensatz ablegen:

Regmodell <- function(data) {
   lm(lifeExp ~ year, data = data)
}

RegModelle <- by_country %>% 
   mutate(Modell = map(data, Regmodell))

head(RegModelle)

A tibble: 6 x 3
 country       data              Modell
 <fct>         <list>             <list>                   
 1 Afghanistan  <tibble [12 x 5]> <lm>    
 2 Albania      <tibble [12 x 5]> <lm> 
 3 Algeria      <tibble [12 x 5]> <lm> 
 4 Angola       <tibble [12 x 5]> <lm> 
 5 Argentina    <tibble [12 x 5]> <lm> 
 6 Australia    <tibble [12 x 5]> <lm> 

Hier verwenden wir map statt map_dbl. map liefert immer eine Liste als Ergebnis zurück, ähnlich wie lapply.

Hier erkennt man die Modell-Ergebnisse nicht direkt – sie sind hinter dem Eintrag <lm> versteckt. Wie können wir sie elegant sichtbar und nutzbar für die Weiterverarbeitung machen? Mit dem R-Paket broom:

RegModelle_tidy <- RegModelle %>% 
   mutate(glance = map(Modell, broom::glance)) %>% 
   unnest(glance)

head(RegModelle_tidy)
broom::glance() wandelt Ergebnisse von Regressionsmodellen in Spalten eines Datensatzes um

That’s it! Mit diesen Werkzeugen kann man sich auf Interpretation und Visualisierung konzentrieren und verliert wenig Zeit mit Copy & Paste für gleichförmige Arbeitsschritte. Angenehme Nebeneffekte: Der Code wird wesentlich kürzer, besser lesbar und leichter zu pflegen und anzupassen.

Zur Vertiefung empfehle ich das Kapitel Many Models in R for Data Science (R4DS) von Hadley Wickham und Garrett Grolemund. Funktionales Programmieren ist Teil meines R-Kurses R-Programmierung für Fortgeschrittene.

Freue mich über Kommentare!