4.1.1 Daten importieren

Weitere Informationen - wie das Codebook (d.h. die Erläuterungen zu den Variablen) - finden Sie auf dem Github-Repo des RKI.

hosp_file <- "https://raw.githubusercontent.com/robert-koch-institut/COVID-19-Hospitalisierungen_in_Deutschland/master/Archiv/2021-11-24_Deutschland_adjustierte-COVID-19-Hospitalisierungen.csv"

hosp_adj <- read_csv(hosp_file)  # adj wie "adjustiert"

4.1.2 EDA

hosp_adj %>% 
  ggplot() +
  aes(x = Datum, y = PS_adjustierte_7T_Hospitalisierung_Inzidenz) +
  geom_line() 

4.2.1 Daten importieren

hosp_latest_file <- "https://github.com/robert-koch-institut/COVID-19-Hospitalisierungen_in_Deutschland/raw/master/Archiv/2021-11-27_Deutschland_COVID-19-Hospitalisierungen.csv"

hosp_latest <- read_csv(hosp_latest_file)

4.2.2 EDA

hosp_latest %>% 
  filter(Bundesland == "Bundesgebiet") %>% 
  ggplot() +
  aes(x = Datum, y = `7T_Hospitalisierung_Inzidenz`) +
  geom_line(aes(color = Altersgruppe)) +
  scale_color_brewer(type = "qual")

hosp_latest %>% 
  filter(Bundesland != "Bundesgebiet") %>% 
  ggplot() +
  aes(x = Datum, y = `7T_Hospitalisierung_Inzidenz`) +
  geom_line(aes(color = Altersgruppe)) +
  scale_color_brewer(type = "qual") +
  facet_wrap(~ Bundesland)

hosp_latest %>% 
  filter(Bundesland != "Bundesgebiet") %>% 
  filter(Datum == max(Datum)) %>% 
  select(-Bundesland_Id, -`7T_Hospitalisierung_Faelle`) %>% 
  mutate(Bundesland = factor(Bundesland)) %>% 
  group_by(Bundesland) %>% 
  summarise(Inz_7T_Hosp = mean(`7T_Hospitalisierung_Inzidenz`)) %>% 
  ggplot() +
  aes(y = fct_reorder(Bundesland, Inz_7T_Hosp), x = Inz_7T_Hosp) +
  geom_col()

hosp_latest <-
hosp_latest %>% 
  mutate(Monat = month(Datum),
         Woche = week(Datum),
         Quartal = quarter(Datum),
         Jahr = year(Datum)) %>% 
  select(Datum, Monat, Woche, Quartal, Jahr, everything())

5.1.1 Daten laden

Von diesem Github-Repo kann man sich die Impfdaten beim RKI herunterladen.

Die Daten werden sogar pro Landkreis bereitgestellt. Wir begnügen uns hier aber mit den Daten pro Bundesland.

vacc_file_latest <- "https://raw.githubusercontent.com/robert-koch-institut/COVID-19-Impfungen_in_Deutschland/master/Archiv/2021-11-27_Deutschland_Impfquoten_COVID-19.csv"

vacc_latest <- read_csv(vacc_file_latest)

glimpse(vacc_latest)

## Rows: 18
## Columns: 22
## $ Datum                   <date> 2021-11-26, 2021-11-26, 2021-11-26, 2021-11-2…
## $ Bundesland              <chr> "Deutschland", "Schleswig-Holstein", "Hamburg"…
## $ BundeslandId_Impfort    <chr> "00", "01", "02", "03", "04", "05", "06", "07"…
## $ Impfungen_gesamt        <dbl> 121165530, 4496736, 2866180, 11883793, 1147170…
## $ Impfungen_gesamt_min1   <dbl> 59121827, 2189818, 1412757, 5868583, 563007, 1…
## $ Impfungen_gesamt_voll   <dbl> 56862343, 2115139, 1373514, 5621902, 544534, 1…
## $ Impfungen_gesamt_boost  <dbl> 8611480, 323441, 176822, 812783, 84998, 200130…
## $ Impfquote_gesamt_min1   <dbl> 71.1, 75.2, 76.3, 73.3, 82.8, 75.1, 71.0, 72.3…
## $ Impfquote_12bis17_min1  <dbl> 51.9, 65.3, 53.2, 59.1, 55.2, 58.8, 52.0, 53.8…
## $ Impfquote_18plus_min1   <dbl> 81.8, 87.2, 88.5, 84.0, 95.6, 86.3, 81.8, 82.9…
## $ Impfquote_18bis59_min1  <dbl> 74.3, 81.0, 83.5, 74.5, 92.2, 78.7, 75.1, 74.2…
## $ Impfquote_60plus_min1   <dbl> 87.2, 90.2, 88.7, 90.3, 94.8, 90.7, 87.3, 89.9…
## $ Impfquote_gesamt_voll   <dbl> 68.4, 72.7, 74.1, 70.2, 80.1, 71.7, 67.5, 68.1…
## $ Impfquote_12bis17_voll  <dbl> 45.9, 58.6, 48.8, 52.6, 50.9, 53.2, 45.6, 46.6…
## $ Impfquote_18plus_voll   <dbl> 79.0, 84.6, 86.2, 80.8, 92.7, 82.6, 78.1, 78.4…
## $ Impfquote_18bis59_voll  <dbl> 75.3, 81.5, 85.8, 76.6, 92.3, 79.4, 74.8, 73.9…
## $ Impfquote_60plus_voll   <dbl> 86.0, 90.0, 87.4, 88.5, 93.3, 88.8, 84.7, 86.5…
## $ Impfquote_gesamt_boost  <dbl> 10.4, 11.1, 9.5, 10.2, 12.5, 11.2, 9.1, 10.0, …
## $ Impfquote_12bis17_boost <dbl> 1.0, 0.4, 1.2, 1.1, 0.5, 1.1, 1.0, 1.2, 0.9, 1…
## $ Impfquote_18plus_boost  <dbl> 12.3, 13.3, 11.4, 12.1, 14.9, 13.4, 10.9, 11.8…
## $ Impfquote_18bis59_boost <dbl> 6.9, 6.5, 6.9, 6.8, 7.6, 7.6, 6.0, 6.3, 5.9, 8…
## $ Impfquote_60plus_boost  <dbl> 22.6, 25.2, 22.8, 22.0, 29.8, 24.4, 20.6, 21.9…

Die Datei gibt die Impfquote für einen Tag (hier 27.11.21) wieder.

5.1.2 EDA

vacc_latest %>% 
  #filter(Bundesland != "Deutschland") %>% 
  mutate(Bundesland = factor(Bundesland)) %>% 
  ggplot() +
  aes(y = fct_reorder(Bundesland, Impfquote_gesamt_voll), 
      x = Impfquote_gesamt_voll) +
  geom_col() +
  labs(title = "Impfquote am 27.11.2021",
       caption = "Datenquelle: RKI")

Und hier als Vektor der Impfquoten (genauer gesagt: Impfquote_gesamt_voll), absteigend sortiert:

vacc_latest_df <-
  vacc_latest %>% 
  filter(!Bundesland %in% c("Bundesressorts", "Deutschland")) %>% 
  mutate(Bundesland = factor(Bundesland)) %>% 
  select(Bundesland, Impfquote_gesamt_voll) %>% 
  arrange(-Impfquote_gesamt_voll)

vacc_latest_df %>%   
  gt() %>% 
  tab_options(
              heading.align = "left"
              )

Und wir speichern die Reihenfolge der Bundesländer nach Impfquote für spätere Verwendung in einem Vektor:

bundeslaender_nach_impfquote <- 
  vacc_latest_df %>% 
  pull(Bundesland) %>% 
  fct_inorder()  # Faktorstufen ordnen entsprechend ihrer empirischen Reihenfolge

bundeslaender_nach_impfquote

##  [1] Bremen                 Saarland               Hamburg               
##  [4] Schleswig-Holstein     Nordrhein-Westfalen    Niedersachsen         
##  [7] Berlin                 Rheinland-Pfalz        Hessen                
## [10] Mecklenburg-Vorpommern Bayern                 Baden-Württemberg     
## [13] Sachsen-Anhalt         Thüringen              Brandenburg           
## [16] Sachsen               
## 16 Levels: Bremen Saarland Hamburg Schleswig-Holstein ... Sachsen

5.2.1 Daten laden 1

vacc_time_file <- "https://raw.githubusercontent.com/robert-koch-institut/COVID-19-Impfungen_in_Deutschland/master/Archiv/2021-11-24_Deutschland_Bundeslaender_COVID-19-Impfungen.csv"

vacc_time <- read_csv(vacc_time_file)

glimpse(vacc_time)

## Rows: 33,773
## Columns: 5
## $ Impfdatum            <date> 2020-12-27, 2020-12-27, 2020-12-27, 2020-12-27, …
## $ BundeslandId_Impfort <chr> "01", "02", "03", "04", "04", "05", "05", "06", "…
## $ Impfstoff            <chr> "Comirnaty", "Comirnaty", "Comirnaty", "Comirnaty…
## $ Impfserie            <dbl> 1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1…
## $ Anzahl               <dbl> 808, 443, 373, 511, 3, 8721, 9, 2712, 2, 1252, 26…

Diese Datei gibt uns nicht direkt die Impfquoten, sondern die Anzahl der Impfungen pro Bundesland pro Tag. Einfacher ist es vielleicht, die Dateien zu suchen, die uns die Impfquoten schon komfortabel zur Verfügung stellen.

5.2.2 Daten laden 2

vacc_dir <- "/Users/sebastiansaueruser/github-repos/COVID-19-Impfungen_in_Deutschland/Archiv"

search_term <- "^\\d{4}-\\d{2}-\\d{2}_Deutschland_Impfquoten_COVID-19.csv$"

vacc_time_files <- 
  list.files(path = vacc_dir,
             pattern = search_term)

head(vacc_time_files)

## [1] "2021-09-17_Deutschland_Impfquoten_COVID-19.csv"
## [2] "2021-09-18_Deutschland_Impfquoten_COVID-19.csv"
## [3] "2021-09-20_Deutschland_Impfquoten_COVID-19.csv"
## [4] "2021-09-21_Deutschland_Impfquoten_COVID-19.csv"
## [5] "2021-09-22_Deutschland_Impfquoten_COVID-19.csv"
## [6] "2021-09-23_Deutschland_Impfquoten_COVID-19.csv"

length(vacc_time_files)

## [1] 58

vacc_path <- "https://raw.githubusercontent.com/robert-koch-institut/COVID-19-Impfungen_in_Deutschland/master/Archiv"

vacc_time_files_w_path <-
  paste0(vacc_path, "/", vacc_time_files)

head(vacc_time_files_w_path)

## [1] "https://raw.githubusercontent.com/robert-koch-institut/COVID-19-Impfungen_in_Deutschland/master/Archiv/2021-09-17_Deutschland_Impfquoten_COVID-19.csv"
## [2] "https://raw.githubusercontent.com/robert-koch-institut/COVID-19-Impfungen_in_Deutschland/master/Archiv/2021-09-18_Deutschland_Impfquoten_COVID-19.csv"
## [3] "https://raw.githubusercontent.com/robert-koch-institut/COVID-19-Impfungen_in_Deutschland/master/Archiv/2021-09-20_Deutschland_Impfquoten_COVID-19.csv"
## [4] "https://raw.githubusercontent.com/robert-koch-institut/COVID-19-Impfungen_in_Deutschland/master/Archiv/2021-09-21_Deutschland_Impfquoten_COVID-19.csv"
## [5] "https://raw.githubusercontent.com/robert-koch-institut/COVID-19-Impfungen_in_Deutschland/master/Archiv/2021-09-22_Deutschland_Impfquoten_COVID-19.csv"
## [6] "https://raw.githubusercontent.com/robert-koch-institut/COVID-19-Impfungen_in_Deutschland/master/Archiv/2021-09-23_Deutschland_Impfquoten_COVID-19.csv"

tic()
vacc_time_df <- 
  vacc_time_files_w_path %>% 
  map_df(read_csv)
toc()

## 7.404 sec elapsed

write_csv(vacc_time_df,
          file = "/Users/sebastiansaueruser/datasets/Covid/vacc_time_df.csv")

write_csv(vacc_time_df,
          file = "static/datasets/vacc_time_df.csv")

vacc_time_df <-
  read_csv("https://raw.githubusercontent.com/sebastiansauer/sesa-blog/main/static/datasets/vacc_time_df.csv")

5.3.1 Impfquote Deutschland im Zeitverlauf

vacc_time_df %>% 
  filter(Bundesland == "Deutschland") %>% 
  ggplot() +
  aes(x = Datum, y = Impfquote_gesamt_voll) +
  geom_line()

5.3.2 Impfquote der Bundesländer im Zeitverlauf

vacc_time_df %>% 
  filter(!Bundesland %in% c("Deutschland", "Bundesressorts")) %>% 
  ggplot() +
  aes(x = Datum, y = Impfquote_gesamt_voll) +
  geom_line() +
  facet_wrap(~ Bundesland, ncol = 4)

vacc_time_df %>% 
  filter(Bundesland != "Deutschland") %>% 
  ggplot() +
  aes(x = Datum, y = Impfquote_gesamt_voll, color = Bundesland) +
  geom_line() +
  scale_color_viridis_d()

Die Daten reichen von … bis … :

vacc_time_df %>% 
  filter(Datum == min(Datum)) %>% 
  pull(Datum)

##  [1] "2021-09-16" "2021-09-16" "2021-09-16" "2021-09-16" "2021-09-16"
##  [6] "2021-09-16" "2021-09-16" "2021-09-16" "2021-09-16" "2021-09-16"
## [11] "2021-09-16" "2021-09-16" "2021-09-16" "2021-09-16" "2021-09-16"
## [16] "2021-09-16" "2021-09-16" "2021-09-16"

vacc_time_df %>% 
  filter(Datum == max(Datum)) %>% 
  pull(Datum)

##  [1] "2021-11-23" "2021-11-23" "2021-11-23" "2021-11-23" "2021-11-23"
##  [6] "2021-11-23" "2021-11-23" "2021-11-23" "2021-11-23" "2021-11-23"
## [11] "2021-11-23" "2021-11-23" "2021-11-23" "2021-11-23" "2021-11-23"
## [16] "2021-11-23" "2021-11-23" "2021-11-23"

Noch ein ganz nettes Detail: Sortieren wir die Legende nach der Höhe der Impfquote.

Zuerst brauchen wir eine nach Impfquote geordnete Liste. Die haben wir oben schon mal erstellt:

bundeslaender_nach_impfquote

##  [1] Bremen                 Saarland               Hamburg               
##  [4] Schleswig-Holstein     Nordrhein-Westfalen    Niedersachsen         
##  [7] Berlin                 Rheinland-Pfalz        Hessen                
## [10] Mecklenburg-Vorpommern Bayern                 Baden-Württemberg     
## [13] Sachsen-Anhalt         Thüringen              Brandenburg           
## [16] Sachsen               
## 16 Levels: Bremen Saarland Hamburg Schleswig-Holstein ... Sachsen

vacc_time_df %>% 
  filter(Bundesland != "Deutschland") %>% 
  ggplot() +
  aes(x = Datum, y = Impfquote_gesamt_voll, color = Bundesland) +
  geom_line() +
  scale_color_viridis_d(breaks = bundeslaender_nach_impfquote)

Häufig ist es sinnvoll, Liniendiagramme mit mehreren Linien gleich mit dem passenden Label zu versehen. Probieren wir das:

vacc_time_df %>% 
  filter(Bundesland != "Deutschland") %>% 
  ggplot() +
  aes(x = Datum, y = Impfquote_gesamt_voll, color = Bundesland) +
  geom_line() +
  scale_color_viridis_d(breaks = NULL) +
  geom_label(data = vacc_time_df %>% filter(Datum == max(Datum)),
             aes(label = Bundesland),
             size = 2)

Hm. Die Position der Labels ist nicht gut. Vielleicht etwas dezenter?

vacc_time_df %>% 
  filter(Bundesland != "Deutschland") %>% 
  ggplot() +
  aes(x = Datum, y = Impfquote_gesamt_voll, color = Bundesland) +
  geom_line() +
  scale_color_viridis_d(breaks = NULL) +
  geom_text(data = vacc_time_df %>% filter(Datum == max(Datum),
                                           Bundesland != "Deutschland"),
             aes(label = Bundesland),
             size = 2,
             nudge_x = +5)

Naja, ok.

Oder vielleicht mit geom_label_repel(), das versucht, Überschneidungen der Labels zu verhindern:

vacc_time_df %>% 
  filter(Bundesland != "Deutschland") %>% 
  ggplot() +
  aes(x = Datum, y = Impfquote_gesamt_voll, color = Bundesland) +
  geom_line() +
  scale_color_viridis_d(breaks = NULL) +
  geom_label_repel(data = vacc_time_df %>% filter(Datum == max(Datum),
                                                  Bundesland != "Deutschland"),
             aes(label = Bundesland))

Hm, alles andere als perfekt. Aber lassen wir es damit gut sein und wenden uns wichtigeren Themen zu. Es ist allzu leicht, sich in (vielleicht ganz unterhaltsamen) Details zu verstricken.

Hier gibt es noch ein interessantes Paket, um Liniendiagramme mit Labels zu versehen.

Probieren wir das doch mal aus:

vacc_time_df %>% 
  filter(Bundesland != "Deutschland") %>% 
  ggplot() +
  aes(x = Datum, 
      y = Impfquote_gesamt_voll, 
      color = Bundesland) +
  geom_line() +
  scale_x_date(limits = c(min(vacc_time_df$Datum),
                                max(vacc_time_df$Datum)+30)) +
  scale_color_viridis_d(breaks = NULL) +
  geom_dl(aes(label = Bundesland),
          method = list(dl.trans(x = x + 0.2), "last.points", cex = 0.8)) 

7.1.1 Im Zeitverlauf

Hier stellt sich die Frage, auf welchen Zeitrahmen wir aggregieren, um eine Korrelation zu berechnen? Je kürzer der Zeitrahmen, desto instabiler die Schätzwerte. Je länger der Zeitrahmen, desto mehr werden Veränderungen im Zeitverlauf unter den Teppich gekehrt. In so einer Situation ist es vielleicht am besten, sich den Koeffizienten in hoher zeitlicher Auflösung anzuschauen, um dann zu entscheiden, wie sehr man aggregieren kann.

korr_vacc_hosp <- 
  hosp_vacc %>% 
  select(Jahr,
         Woche, Impfquote_gesamt_voll, 
         Inz_7T_Hosp = `7T_Hospitalisierung_Inzidenz`) %>% 
  filter(Jahr == 2021) %>% 
  group_by(Woche) %>% 
  summarise(Korr_Woche = cor(Impfquote_gesamt_voll,
                             Inz_7T_Hosp))

korr_vacc_hosp_joined <- 
  korr_vacc_hosp %>% 
  left_join(hosp_vacc %>% select(Datum, Woche, Monat, Quartal)) %>% 
  group_by(Woche) %>% 
  filter(row_number() == 1)

korr_vacc_hosp_joined %>% 
  mutate(Quartal = factor(Quartal)) %>% 
  ggplot() +
  aes(y = Korr_Woche) +
  geom_rect(aes(
    x = NULL,
    y = NULL,
    xmin = Woche-1, 
    xmax = Woche, 
    ymin = -Inf,
    ymax = +Inf,
    fill = Quartal),
    alpha = .5) +
  scale_fill_manual(values = c("grey", "lightgreen", "darkgreen", "indianred4")) +
    geom_hline(yintercept = 0, color = "white", size = 2) +
  geom_point(alpha = .7, aes(x = Woche)) +
  geom_line(aes(x = Woche)) +
  labs(title = "Korrelation von Hospitalisierung und Impfquote pro Woche",
       caption = "Quartale des Jahres 2021 sind eingefärbt. Datenquelle: RKI",
       x = "Kalenderwoche 2021",
       y = "Korrelation (Wochenbasis)")

korr_vacc_hosp_monat <- 
  hosp_vacc %>% 
  select(Jahr,
         Monat, 
         Impfquote_gesamt_voll, 
         Inz_7T_Hosp = `7T_Hospitalisierung_Inzidenz`) %>% 
  filter(Jahr == 2021) %>% 
  group_by(Monat) %>% 
  summarise(Korr_Monat = cor(Impfquote_gesamt_voll,
                             Inz_7T_Hosp))

korr_vacc_hosp_monat %>% 
  gt() %>% 
  fmt_number(2, decimals = 2)

korr_vacc_hosp_monat_joined <- 
  korr_vacc_hosp_monat %>% 
  left_join(hosp_vacc %>% select(Datum, Monat, Quartal)) %>% 
  group_by(Monat) %>% 
  filter(row_number() == 1)

korr_vacc_hosp_monat_joined %>% 
  mutate(Quartal = factor(Quartal)) %>% 
  ggplot() +
  aes(y = Korr_Monat) +
  geom_rect(aes(
    x = NULL,
    y = NULL,
    xmin = Monat-1, 
    xmax = Monat, 
    ymin = -Inf,
    ymax = +Inf,
    fill = Quartal),
    alpha = .5) +
  scale_fill_manual(values = c("grey", "lightgreen", "darkgreen", "indianred4")) +
    geom_hline(yintercept = 0, color = "white", size = 2) +
  geom_point(alpha = .7, aes(x = Monat)) +
  geom_line(aes(x = Monat)) +
  labs(title = "Korrelation von Hospitalisierung und Impfquote pro Monat",
       caption = "Quartale des Jahres 2021 sind eingefärbt. Datenquelle: RKI",
       x = "Kalenderwoche 2021",
       y = "Korrelation (Wochenbasis)")

7.1.2 Korrelation nach Bundesland

hosp_vacc %>% 
  filter(Monat == 11, Bundesland != "Bundesgebiet") %>% 
  select(I = Impfquote_gesamt_voll,
         H = `7T_Hospitalisierung_Inzidenz`,
         B = Bundesland) %>% 
  group_by(B) %>% 
  summarise(I = mean(I),
            H = mean(H)) %>% 
  ggplot() +
  aes(x = I,
      y = H) +
  geom_point() +
  geom_smooth(method = "lm") +
  labs(x = "Impfquote (vollständig)",
       y = "7Tage-Hospitalisierungsquote",
       title = "Zusammenhang von Impfen und Hospitalisierung",
       subtitle = "Daten für November 2021, alle Altersgruppe",
       caption = "Datenquelle: RKI")

state_names <- 
  tibble::tribble(
    ~Bundesland, ~Kurzform,                          ~state,
    "Deutschland",      "DE",                       "Germany",
    "Baden-Württemberg",      "BW",             "Baden-Württemberg",
    "Bayern",      "BY",          "Bavaria (FreeState)",
    "Berlin",      "BE",                        "Berlin",
    "Brandenburg",      "BB",                   "Brandenburg",
    "Bremen",      "HB",       "Bremen (Hanseatic City)",
    "Hamburg",      "HH",      "Hamburg (Hanseatic City)",
    "Hessen",      "HE",                         "Hesse",
    "Mecklenburg-Vorpommern",      "MV", "Mecklenburg-Western Pomerania",
    "Niedersachsen",      "NI",                  "Lower Saxony",
    "Nordrhein-Westfalen",      "NW",        "North Rhine-Westphalia",
    "Rheinland-Pfalz",      "RP",          "Rhineland-Palatinate",
    "Saarland",      "SL",                      "Saarland",
    "Sachsen",      "SN",           "Saxony (Free tate)",
    "Sachsen-Anhalt",      "ST",                 "Saxony-Anhalt",
    "Schleswig-Holstein",      "SH",            "Schleswig-Holstein",
    "Thüringen",      "TH",        "Thuringia (Free State)",
    "Bundesgebiet",                "De", "Federal area"
  )

state_names %>% 
  gt()

hosp_vacc2 <- 
  hosp_vacc %>% 
  left_join(state_names)

hosp_vacc_bundesland_nov <- 
hosp_vacc %>% 
  filter(Monat == 11, Bundesland != "Bundesgebiet") %>% 
  select(I = Impfquote_gesamt_voll,
         H = `7T_Hospitalisierung_Inzidenz`,
         Bundesland) %>% 
  group_by(Bundesland) %>% 
  summarise(I = mean(I),
            H = mean(H)) %>% 
  left_join(state_names %>% select(-state))


hosp_vacc_bundesland_nov %>% 
  ggplot() +
  aes(x = I,
      y = H) +
  geom_point(alpha = .5) +
  geom_smooth(method = "lm") +
  labs(x = "Impfquote (vollständig)",
       y = "7Tage-Hospitalisierungsquote",
       title = "Zusammenhang von Impfen und Hospitalisierung",
       subtitle = "Daten für November 2021, alle Altersgruppe",
       caption = "Datenquelle: RKI") +
  geom_text(aes(label = Kurzform))

hosp_vacc_bundesland_nov %>% 
  gt() %>% 
  fmt_number(3, decimals = 2)

7.1.3 Pro Quartal

Aggregieren wir den Datensatz pro Quartal und pro Bundesland;:

hosp_vacc_bundesland_pro_Quartal <- 
  hosp_vacc %>% 
  filter(Bundesland != "Bundesgebiet") %>% 
  select(I = Impfquote_gesamt_voll,
         H = `7T_Hospitalisierung_Inzidenz`,
         Bundesland,
         Quartal) %>% 
  group_by(Bundesland, Quartal) %>% 
  summarise(I = mean(I),
            H = mean(H)) %>% 
  left_join(state_names)

So sieht die zugehörige Tabelle aus (die ersten paar Zeilen)

hosp_vacc_bundesland_pro_Quartal %>% 
  ungroup() %>% 
  slice_head(n = 10) %>% 
  gt() %>% 
  fmt_number(3, decimals = 2)

Dann updaten wir das letzte Diagramm, nur eben dieses Mal erstellen pro Quartal ein Diagramm.

Zuvor legen wir noch schönere Bezeichnungen der Facetten an, nicht nur 1,2,3,4, sondern “1. Quartal,”2. Quartal” etc.:

Quartal_values <- c("1. Quartal", "2. Quartal",
                    "3. Quartal", "4. Quartal")
names(Quartal_values) <- c(1:4)

Ansonsten ist es das gleiche Vorgehen wie eben auch:

hosp_vacc_bundesland_pro_Quartal %>% 
  ggplot() +
  aes(x = I,
      y = H) +
  geom_point(alpha = .5) +
  geom_smooth(method = "lm") +
  labs(x = "Impfquote (vollständig)",
       y = "7Tage-Hospitalisierungsquote",
       title = "Zusammenhang von Impfen und Hospitalisierung",
       subtitle = "Aufgeteilt nach Quartale, aggregiert über alle Altersgruppe",
       caption = "Datenquelle: RKI") +
  geom_text(aes(label = Kurzform)) +
  facet_wrap(~ Quartal, labeller = labeller(Quartal = Quartal_values))

Wir sehen, dass es im 3. Quartal in der Inzidenz kaum Streuung gibt, was die sehr geringe Korrelation erklärt.

8.2.1 Kausalmodell 1

Kausalmodelle kann man mit Directed Acyclic Graphs (DAGs) darstellen.

Also, ein Modell, das mit unserem impliziten Kausalmodell d’accord geht, sagt, sinngemäß:

“Die Hospitalisierungsquote H wird von der Impfquote H beeinflusst. Natürlich spielen weitere, unbekannte Faktoren U eine Rolle; die haben aber nichts mit H zu tun.”

Als DAG:

dag1 <-
  dagitty(
    'dag {
    I [pos = "0,1"]
    H [pos = "1,0"]
    U [pos = "0, -1"]
    I -> H
    U -> H
    }
    '
  )
plot(dag1)

Nur wenn U nichts mit I zu tun hat (zumindest in einem bestimmten Sinn nichts zu tun hat), können wir mit der Korrelation von I und H den kausalen Effekt wahrnehmen, dem wir vermutlich auf der Fährte sind, zumindest ist das ein lohnendes Ziel.

8.2.2 Kausalmodell 2

Allerdings könnte die Wahrheit auch so aussehen:

dag2 <-
  dagitty(
    'dag {
    I [pos = "0,1"]
    H [pos = "1,0"]
    U1 [pos = "0, -1"]
    U2 [pos = "-1,0"]
    I -> H
    U1 -> H
    U2 -> U1
    U2 -> I
    }
    '
  )
plot(dag2)

Jetzt ist die Korrelation (allgemeiner: der statistische Zusammenhang) von I mit H nicht mehr so unbedarft. U2 hat seine Finger im Spiel. U1 ist nämlich über U2 (eine weitere unbekannte oder ungemessene Größe) mit I (und weiter nach H) verbunden. Man spricht von einer Konfundierung: In der Assoziation von I mit H mischt sich ein weiterer kausaler Pfad: Der von U1->U2->I->H. Damit ist die Assoziation von I mit H die Summe (oder besser: die Mischung) von zwei Pfaden und man kann nicht ohne Weiteres sagen, welcher Pfad wieviel der Assoziation versursacht. Schlechte Nachrichten. Wenn wir aber U1 und U2 messen können, können wir wiederum Hoffnung schöpfen, die kausalen Pfade klar zu ziehen.

8.2.3 Interpretation hängt vom Kausalmodell ab

Ob wir also diese (oder eine andere) Korrelation interpretieren können (und interpretieren, so glaube ich, ist oft oder in hohem Maße eine kausale Sache) hängt vom Kausalmodell ab, an das wir glauben.

Natürlich sind noch viele andere Kausalmodelle denkbar. Ja, die Sache ist kompliziert.

Das ist übrigens der Grund, warum die Menschen Experimente erfunden haben. Mit Experimenten kann man diese gordischen Kausalknoten durchschneiden, und dann vergleichsweise einfach die Kausalstrukturen klären. Leider sind Experimente nicht immer möglich.

Glücklicherweise im Fall der Frage “Was bringt Impfen gegen Corona?” schon. Da gibt es Experimente, die diese Frage untersuchen. Vermutlich sind diese Studien die besten Möglichkeiten, die kausale Frage nach der Wirksamkeit von Impfung zu klären.

Natürlich, wenn man an ein Kausalmodell glaubt - bzw. es gut begründen kann - das die Vergleichbarkeit der Gruppen nahelegt, dann ist wieder alles gut. Wir wären also zurück beim Kausalmodell 1, oben.