Die Lebenserwartung der Menschen variiert, wie jeder weiß. Die Variation wird u.a. erklärt durch das Land, in dem jemand lebt, aber auch durch die Zeit – in den letzten Jahrzehnten ist die Lebenserwartung zum Glück deutlich gestiegen. Allerdings nicht in jedem Jahr und in jedem Land in gleichem Maße. In dieser Fallstudie soll der Einfluss von Land, Region und Zeit auf die Lebenserwartung untersucht werden. Konkret soll die Lebenserwartung im Alter von 10 Jahren modelliert werden.
Die Daten stammen von Our World in Data, genauer gesagt aus dieser Analyse. Weitere Quellen sind hier dokumentiert.
Es handelt sich um eine retrospektive Beobachtungsstudie. Die Daten wurden ex post erhoben.
Zuerst mosaic
laden; nicht vergessen! Falls eines der Pakete auf Ihrer Maschine nicht installiert ist, bitte mit install.packages("name_des_pakets")
installieren.
library(mosaic)
library(countrycode)
library(rnaturalearth)
library(tidyverse)
So können die Daten komfortabel in R importieren:
exp_raw <- read.csv("https://raw.githubusercontent.com/sebastiansauer/Statistiklehre/main/data/life-expectancy-at-age-10.csv")
head(exp_raw)
Ich habe die Tabelle exp_raw
genannt, um zu verdeutlichen, dass es sich um die “rohen”, unbearbeiteten Daten handelt. Häufig muss man die Daten noch modifizieren, da bietet es sich an, klar anzuzeigen, was die unbehandelten Daten und was die behandelten Daten sind.
Die Spalte e10..years.
fasst die Lebenserwartung zum Alter von 10 Jahren. Allerdings ist der Name der Spalte nicht so schön. Benennen wir die Spalte um:
exp <- exp_raw %>%
rename(e10 = e10..years.) # neu = alt
Es wäre doch schön, für jedes Land zu wissen, zu welchem Kontinent (oder Weltregion) es gehört. Natürlich kann man das händisch hinzufügen (z.B. in Excel). Aber komfortabler ist es, wenn man solche schnöde Arbeit den Computer erledigen lässt. In R gibt es ein Paket namens countrycode
(und eine Funktion mit gleichem Namen), das die Arbeit für uns übernimmt. Hier bilden wir den Landesnamen auf das zugehörige Kontinent ab; alternativ hätten wir auch die Spalte Code
verwenden können. Allerdings ist nicht ganz klar, um welchen Code es sich handelt (vielleicht ISO). Das zeigt, das ein gutes Codebuch zu jedem Datensatz gehört!
exp2 <- exp %>%
mutate(continent = countrycode(sourcevar = Entity, # Spalte, in der das Land steht
origin = "country.name", # Eingage: Name (nicht z.B. ISO-Code) des Landes
destination = "continent")) # gewünschte Ausgabe: Name des Kontinents
## Warning: Problem with `mutate()` input `continent`.
## ℹ Some values were not matched unambiguously: Channel Islands, Micronesia (country), Timor, World
##
## ℹ Input `continent` is `countrycode(sourcevar = Entity, origin = "country.name", destination = "continent")`.
## Warning in countrycode(sourcevar = Entity, origin = "country.name", destination = "continent"): Some values were not matched unambiguously: Channel Islands, Micronesia (country), Timor, World
Dazu haben wir eine Spalte mit mutate()
angelegt, die den Namen des Kontinents für jedes Land fasst. R informiert uns, dass einige Länder nicht einem Kontinent zugeordnet werden konnten. Schauen wir uns das genauer an. Wie viele Fälle liefern fehlende Werte zurück?
exp2 %>%
summarise(is_na_sum = is.na(continent) %>% sum())
Schauen wir uns ein Beispiel näher an:
exp2 %>%
filter(Entity == "Timor")
Vielleicht klappt es mit der Übersetzung der Spalte Code
besser? Gehen wir mal davon aus, dass es sich um ISO-3661-Codes handelt.
exp2 <- exp2 %>%
mutate(continent2 = countrycode(sourcevar = Code, # Spalte, in der das Land steht
origin = "iso3c", # Eingage: Name (nicht z.B. ISO-Code) des Landes
destination = "continent")) # gewünschte Ausgabe: Name des Kontinents
## Warning: Problem with `mutate()` input `continent2`.
## ℹ Some values were not matched unambiguously: OWID_CIS, OWID_WRL
##
## ℹ Input `continent2` is `countrycode(sourcevar = Code, origin = "iso3c", destination = "continent")`.
## Warning in countrycode(sourcevar = Code, origin = "iso3c", destination = "continent"): Some values were not matched unambiguously: OWID_CIS, OWID_WRL
exp2 %>%
summarise(is_na_sum = is.na(continent2) %>% sum())
Besser; 60 fehlende Werte. Aber welche Fälle sind übrig? Schauen wir uns die mal an:
exp2 %>%
filter(is.na(continent2)) %>%
distinct(Entity)
Auf Deutsch übersetzt heißt die Syntax oben:
Hey R, nimm die Tabelle exp2 UND DANN
filtere die Fälle, die keine Werte haben für continent2 UND DANN
zeige alle verschiedenen Werte für Entity
Laut dieser Quelle gehören die Channel Islands zu UK (GB). Es ist vielleicht pragmatisch, dies einfach zu übernehmen.
exp3 <- exp2 %>%
mutate(continent2 = as.character(continent2)) %>%
mutate(continent2 = case_when(
Entity == "Channel Islands" ~ "Europe",
Entity == "World" ~ "World",
TRUE ~ continent2))
Prüfen wir, ob jetzt alle fehlenden Werte bei continent2
beseitigt sind:
exp3 %>%
filter(is.na(continent2))
Gut.
Für welche Jahre liegen Daten vor?
exp3 %>%
distinct(Year) %>%
pull()
## [1] 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 2010 2015 2020
## [16] 2025 2030 2035 2040 2045 2050 2055 2060 2065 2070 2075 2080 2085 2090 2095
## [31] 1750 1755 1760 1765 1770 1775 1780 1785 1790 1795 1800 1805 1810 1815 1820
## [46] 1825 1830 1835 1840 1845 1850 1855 1860 1865 1870 1875 1880 1885 1890 1895
## [61] 1900 1905 1910 1915 1920 1925 1930 1935 1940 1945
Möchte man die Daten nicht als Tabelle haben (braucht weniger Platz im Output), sondern aus der Tabelle eine Spalte herausziehen, kann man dafür den Befehl pull()
verwenden.
Man sieht, dass die Daten offenbar nur in Fünf-Jahres-Schritten erhoben wurden. Wählen wir das Jahr 2015:
exp_2015 <- filter(exp3, Year == 2015)
Berechnen wir die typischen Deskriptivstatistiken pro Kontinent:
favstats(e10 ~ continent2, data = exp_2015)
Welche Länder gehören eigentlich zu “Oceania”?
filter(exp_2015, continent == "Oceania") %>%
pull(Entity)
## [1] "Australia" "Fiji" "French Polynesia" "Guam"
## [5] "Kiribati" "New Caledonia" "New Zealand" "Papua New Guinea"
## [9] "Samoa" "Solomon Islands" "Tonga" "Vanuatu"
Vermutlich würde es Sinn machen, Nord- und Südamerika getrennt auszuweisen.
gf_boxplot(e10 ~ continent2,
data = exp_2015) %>%
gf_point(stat = "summary",
color = "red",
size = 5) %>%
gf_jitter(width = .1,
alpha = .3) %>%
gf_labs(x= "Kontinent", y = "Lebenserwartung im Alter von 10 Jahren",
caption = "Daten von 'Our World in Data'",
title = "Lebenserwartung nach Kontinenten im Jahr 2015")
## No summary function supplied, defaulting to `mean_se()`
Da World
kein Kontinent ist, nehmen wir es heraus, aber zeichnen eine horizontale Linie für den Wert der Welt insgesamt.
exp_2015 %>%
filter(continent2 != "World") %>%
gf_boxplot(e10 ~ continent2) %>%
gf_point(stat = "summary",
color = "red",
size = 5) %>%
gf_jitter(width = .1,
alpha = .3) %>%
gf_labs(x= "Kontinent", y = "Lebenserwartung im Alter von 10 Jahren",
caption = "Daten von 'Our World in Data'. Die horizontale Linie gibt den weltweiten Wert wieder.",
title = "Lebenserwartung nach Kontinenten im Jahr 2015.") %>%
gf_hline(yintercept = ~65.3, linetype = "dashed")
## No summary function supplied, defaulting to `mean_se()`
Welche drei Länder weisen die geringste/höchste Lebenserwartung im Jahr 2015 auf?
exp_2015 %>%
arrange(e10) %>% # sortiert von klein nach groß bzgl. e10
slice(1:3)
slice()
“schneidet” eine “Scheibe” an Zeilen heraus, hier 1 bis 3.
In gleicher Manier:
exp_2015 %>%
arrange(-e10) %>%
slice
Das Minuszeichen dreht die Sortierreihenfolge um, d.h. von groß zu klein.
cor(e10 ~ Year, data = exp3)
## [1] NA
Oh, es gibt fehlende Werte im Datensatz, daher streckt der Befehl cor()
alle Viere von sich. Löschen wir mal alle Zeilen mit fehlenden Werten und hoffen, das wir nicht viele Daten verlieren:
exp_ohne_na <- exp3 %>%
na.omit()
Hm, ca. 500 Zeilen. Vielleicht doch lieber so:
cor(e10 ~ Year, data = exp3 %>% na.omit())
## [1] 0.7531835
Inferenzstatistik macht hier wenig Sinn, da es unklar klar ist, auf welche Grundgesamtheit verallgemeinert werden soll. Alle Länder der Welt sind ja schon enthalten. Zu argumentieren, wir verallgemeinern auf alle Zeiten (Jahre) ist gewagt, denn unsere Stichprobe an Jahren ist sicher alles andere als eine Zufallsstichprobe aus der Menge aller Jahre der Welt … Verzichten wir also auf eine Inferenzstatistik.
gf_point(e10 ~ Year, data = exp3, alpha = .2) %>%
gf_smooth() %>%
gf_facet_wrap(~ continent2)
## `geom_smooth()` using method = 'gam'
## Warning: Removed 419 rows containing non-finite values (stat_smooth).
## Warning: Removed 419 rows containing missing values (geom_point).
gf_facet_wrap()
macht ein “Teil-Bildchen” (Facette) pro Wert von Kontinent. Allerdings liegen die Daten für die meisten Länder erst ab ca. 1950 vor. Und die Prognosen in die Zukunft sehen wir mal kritisch. Vorhersagen sind bekanntlich schwierig. Vor allem, wenn sie die Zukunft betreffen, heißt es … gf_smooth()
legt eine Kurve in “die Mitte” des Punkteschwarms: Für jeden X-Werte wird der mittlere Y-Wert berechnet (nahe X-Werte fließen auch noch etwas ein), und dann werden die Punkte mit einer Linie verbunden.
exp3 %>%
filter(Year > 2016 & Year > 1949) %>%
gf_jitter(e10 ~ Year, alpha = .2) %>%
gf_smooth() %>%
gf_facet_wrap(~ continent2)
## `geom_smooth()` using method = 'loess'
Wir sehen durchweg einen Anstieg der Lebenserwartung - sehr erfreulich! Freilich fußt der Anstieg auf unterschiedlichen Sockeln, also Ausgangsniveaus. Ist die Höhe des Anstiegs unterschiedlich je Kontinent? Für diese Frage modellieren wir die Lebenserwartung als Funktion des Kontinents und des Jahres.
Bleiben wir beim realistischen Teil der Daten:
exp_1950_2015 <- exp3 %>%
filter(Year > 2016 & Year > 1949)
Berechnen wir dazu einige Modelle und vergleichen diese dann:
Year
als PrädikatorYear
und continent2
, aber ohne InteraktionseffektYear*continent2
, also mit Interaktionseffektlm0 <- lm(e10 ~ 1, data = exp_1950_2015)
lm1 <- lm(e10 ~ Year, data = exp_1950_2015)
lm2 <- lm(e10 ~ Year + continent2, data = exp_1950_2015)
lm3 <- lm(e10 ~ Year*continent2, data = exp_1950_2015)
Betrachten wir die \(R^2\)-Werte jedes Modells:
rsquared(lm0) %>% str()
## num 0
r2s_vector <- c(rsquared(lm0), rsquared(lm1), rsquared(lm2), rsquared(lm3))
r2s_vector %>% str()
## num [1:4] 0 0.219 0.531 0.532
Um die Werte in einer kleiner Grafik zu zeigen, erstellen wir zuerst eine Tabelle (tibble()
), weil gf_XXX()
nur Spaß an Tabellen hat.
r2s <- tibble(Modellname = c("lm0", "lm1", "lm2", "lm3"),
r2s = r2s_vector)
r2s
Und jetzt das Diagramm:
gf_point(r2s ~ Modellname, data = r2s, size = 4) %>%
gf_line(group = ~1) # Es gibt nur eine Gruppe, also alle Punkte sollen verbunden werden
Wie man sieht, ist der Zuwachs an erklärter Varianz von lm2
auf lm3
nicht groß. Der Interaktionseffekt scheint also nicht stark zu sein. Besser wir verzichten auf ihn und resümieren, dass die Daten eher die Hypothese stützen, dass der Zuwachs in Lebenserwartung gleich oder ähnlich groß ist zwischen den Kontinenten.
Wir können die Güte der Modelle bzw. die Unterschiede der Güte auf Signifikanz testen. Genau genommen werden die Likelihoods der Modelle verglichen und anhand der Freiheitsgrade normiert. Aber lassen wir die technischen Details an dieser Stelle. Diesen Modellvergleich können anhand der Funktion anova()
durchführen, der uns die Quadratsummen zurückliefert und einen p-Wert (letzte Spalte).
anova(lm0, lm1, lm2, lm3)
Wir sehen, dass das letzte Modell (Model 4
) nicht signifikant besser ist als das vorherige Modell. Das bestärkt unseren Schluss, die \(H_0\) der Gleichheit der Modellgüte der beiden Modelle, nicht zu verwerfen.
Analysiert man Länder (oder allgemeiner: Gegenden) der Erde, so bietet sich eine Geo-Visualisierung - wie eine Weltkarte - an. Das lässt sich recht einfach mit R bewerkstelligen.
Im Paket rnaturalearth
sind die Weltkarten gespeichert. Ziehen wir uns die Länder der “natural earth” (ne) zunächst heraus:
world <- ne_countries(scale = "medium", returnclass = "sf")
Die Tabelle ist sehr umfangreich; viele interessante Daten sind vorhanden. Die Geo-Daten verstecken sich in der Spalte geometry
. Die technischen Details interessieren uns hier nicht.
str(world)
## Classes 'sf' and 'data.frame': 241 obs. of 64 variables:
## $ scalerank : int 3 1 1 1 1 3 3 1 1 1 ...
## $ featurecla: chr "Admin-0 country" "Admin-0 country" "Admin-0 country" "Admin-0 country" ...
## $ labelrank : num 5 3 3 6 6 6 6 4 2 6 ...
## $ sovereignt: chr "Netherlands" "Afghanistan" "Angola" "United Kingdom" ...
## $ sov_a3 : chr "NL1" "AFG" "AGO" "GB1" ...
## $ adm0_dif : num 1 0 0 1 0 1 0 0 0 0 ...
## $ level : num 2 2 2 2 2 2 2 2 2 2 ...
## $ type : chr "Country" "Sovereign country" "Sovereign country" "Dependency" ...
## $ admin : chr "Aruba" "Afghanistan" "Angola" "Anguilla" ...
## $ adm0_a3 : chr "ABW" "AFG" "AGO" "AIA" ...
## $ geou_dif : num 0 0 0 0 0 0 0 0 0 0 ...
## $ geounit : chr "Aruba" "Afghanistan" "Angola" "Anguilla" ...
## $ gu_a3 : chr "ABW" "AFG" "AGO" "AIA" ...
## $ su_dif : num 0 0 0 0 0 0 0 0 0 0 ...
## $ subunit : chr "Aruba" "Afghanistan" "Angola" "Anguilla" ...
## $ su_a3 : chr "ABW" "AFG" "AGO" "AIA" ...
## $ brk_diff : num 0 0 0 0 0 0 0 0 0 0 ...
## $ name : chr "Aruba" "Afghanistan" "Angola" "Anguilla" ...
## $ name_long : chr "Aruba" "Afghanistan" "Angola" "Anguilla" ...
## $ brk_a3 : chr "ABW" "AFG" "AGO" "AIA" ...
## $ brk_name : chr "Aruba" "Afghanistan" "Angola" "Anguilla" ...
## $ brk_group : chr NA NA NA NA ...
## $ abbrev : chr "Aruba" "Afg." "Ang." "Ang." ...
## $ postal : chr "AW" "AF" "AO" "AI" ...
## $ formal_en : chr "Aruba" "Islamic State of Afghanistan" "People's Republic of Angola" NA ...
## $ formal_fr : chr NA NA NA NA ...
## $ note_adm0 : chr "Neth." NA NA "U.K." ...
## $ note_brk : chr NA NA NA NA ...
## $ name_sort : chr "Aruba" "Afghanistan" "Angola" "Anguilla" ...
## $ name_alt : chr NA NA NA NA ...
## $ mapcolor7 : num 4 5 3 6 1 4 1 2 3 3 ...
## $ mapcolor8 : num 2 6 2 6 4 1 4 1 1 1 ...
## $ mapcolor9 : num 2 8 6 6 1 4 1 3 3 2 ...
## $ mapcolor13: num 9 7 1 3 6 6 8 3 13 10 ...
## $ pop_est : num 103065 28400000 12799293 14436 3639453 ...
## $ gdp_md_est: num 2258 22270 110300 109 21810 ...
## $ pop_year : num NA NA NA NA NA NA NA NA NA NA ...
## $ lastcensus: num 2010 1979 1970 NA 2001 ...
## $ gdp_year : num NA NA NA NA NA NA NA NA NA NA ...
## $ economy : chr "6. Developing region" "7. Least developed region" "7. Least developed region" "6. Developing region" ...
## $ income_grp: chr "2. High income: nonOECD" "5. Low income" "3. Upper middle income" "3. Upper middle income" ...
## $ wikipedia : num NA NA NA NA NA NA NA NA NA NA ...
## $ fips_10 : chr NA NA NA NA ...
## $ iso_a2 : chr "AW" "AF" "AO" "AI" ...
## $ iso_a3 : chr "ABW" "AFG" "AGO" "AIA" ...
## $ iso_n3 : chr "533" "004" "024" "660" ...
## $ un_a3 : chr "533" "004" "024" "660" ...
## $ wb_a2 : chr "AW" "AF" "AO" NA ...
## $ wb_a3 : chr "ABW" "AFG" "AGO" NA ...
## $ woe_id : num NA NA NA NA NA NA NA NA NA NA ...
## $ adm0_a3_is: chr "ABW" "AFG" "AGO" "AIA" ...
## $ adm0_a3_us: chr "ABW" "AFG" "AGO" "AIA" ...
## $ adm0_a3_un: num NA NA NA NA NA NA NA NA NA NA ...
## $ adm0_a3_wb: num NA NA NA NA NA NA NA NA NA NA ...
## $ continent : chr "North America" "Asia" "Africa" "North America" ...
## $ region_un : chr "Americas" "Asia" "Africa" "Americas" ...
## $ subregion : chr "Caribbean" "Southern Asia" "Middle Africa" "Caribbean" ...
## $ region_wb : chr "Latin America & Caribbean" "South Asia" "Sub-Saharan Africa" "Latin America & Caribbean" ...
## $ name_len : num 5 11 6 8 7 5 7 20 9 7 ...
## $ long_len : num 5 11 6 8 7 13 7 20 9 7 ...
## $ abbrev_len: num 5 4 4 4 4 5 4 6 4 4 ...
## $ tiny : num 4 NA NA NA NA 5 5 NA NA NA ...
## $ homepart : num NA 1 1 NA 1 NA 1 1 1 1 ...
## $ geometry :sfc_MULTIPOLYGON of length 241; first list element: List of 1
## ..$ :List of 1
## .. ..$ : num [1:10, 1:2] -69.9 -69.9 -69.9 -70 -70.1 ...
## ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:63] "scalerank" "featurecla" "labelrank" "sovereignt" ...
Wir haben den Befehl angewiesen, dass wir keine super genauen Grenzen benötigen, sondern dass uns eine mittlere Genauigkeit ausreicht. Das Datenformat nennt sich hier sf
(simple feature), was zur Zeit ein Standardformat ist, um geografische Daten (GPS-Koordinaten für Grenzen und Punkte etc.) zu speichern. Aber wie gesagt, die technischen Details sind hier nicht interessant für uns. Mit dem Befehl gf_sf
lassen sich dann Geo-Daten plotten:
gf_sf(data = world)
Damit wir die Länder entsprechend ihrer Lebenserwartung einfärben können, müssen wir die Tabelle mit den Geo-Daten und die Tabelle mit den Lebensdaten zusammenführen. Das besorgt der Befehl full_join()
. Damit der Befehl weiß, welche Zeilen zusammen gehören, erklären wir ihm “Füge gleiche Länder zusammen, dazu schaue in der Spalte Entity
bzw. in der Spalte sovereignt
(Tabelle world
)”.
exp_2015_joined <- world %>%
left_join(exp_2015, by = c("sovereignt" = "Entity"))
Das eigentliche Plotten ist schnell erledigt:
gf_sf(fill = ~e10, data = exp_2015_joined)
Praktischerweise findet gg_sf()
die Geo-Daten selbständig in der Tabelle, wir müssen nicht extra erklären, welche Spalte gemeint ist. Sehr komfortabel.
Das Farbschema könnte noch schöner sein:
gf_sf(fill = ~e10, data = exp_2015_joined) %>%
gf_refine(scale_fill_viridis_c()) %>%
gf_labs(fill = "Lebenserwartung",
title = "Lebenserwatung im Alter von 10 Jahren")
Viel besser. Mit gf_refine()
“verfeinern` wir das Diagramm. In diesem Fall besteht das Verfeinern im Ändern des Farbschemas (wie nehmen Viridis), um die Füllfarbe zu ändern. Da es sich um eine kontinuierliche Variable handelt (d.h. metrisch), soll ein kontinuierliches Farbschema (mit fließenden, weichen Übergangen) verwendet werden (”_c" wie continuous).
Halt! Für die USA gibt es keine Werte! Kann das sein? Das ist ein gutes Beispiel dafür, dass eine Datenanalyse ein iterativer Prozess ist, d.h. einzelne Schritte - wie Datenvorverarbeitung - müssen immer wieder ausgeführt werden.
Versuchen wir, die USA in den Daten zu finden:
exp3 %>%
filter(str_detect(Entity, "USA|United States"))
die USA existieren also doch. Übrigens: str_detect(spalte, suchterm)
liefert für jeden Wert von spalte
zurück, ob sich der Suchterm (suchterm
) darin befindet. Man bekommt also einen Vektor mit TRUE, FALSE, TRUE, ..." und so weiter.
filter()` erlaubt, so einen logischen Vektor als Grundlage des Filterns herzunehmen.
Schauen wir, wie viele Länder zu unseren Suchterm passen (also TRUE
zurückliefern):
exp3 %>%
filter(str_detect(Entity, "USA|United States")) %>%
pull(Entity) %>%
unique() # zeigt nur unterschiedliche (unique) Werte an
## [1] "United States" "United States Virgin Islands"
Aha, “United States” und die “Virgin Islands”.
Schauen wir mal nach, ob es die USA auch in der Tabelle world
gibt:
world %>%
select(sovereignt) %>%
filter(str_detect(sovereignt, "USA|United States")) %>%
pull(sovereignt) %>%
unique()
## [1] "United States of America"
Ah - die USA heißen hier “United States of America”. Das ist nicht exakt zu “United States” aus exp3
. Daher hat das Verheiraten (der full_join()
) oben nicht funktioniert. Besser wir nehmen den ISO-Code zum Vereinigen. Codes sind (hoffentlich) “bruchsicher”, so dass das Vereinigen klappen sollte. Der relevant Code scheint in Spalte adm0_a3
zu Hause zu sein:
exp_2015_joined <- world %>%
left_join(exp_2015, by = c("adm0_a3" = "Code"))
Und Plotten:
gf_sf(fill = ~e10, data = exp_2015_joined) %>%
gf_refine(scale_fill_viridis_c())
Hat funktioniert!