5 Multiple Regression I: Grundlagen linearer Modelle

📢 Zielsetzung dieser Einheit

In dieser Einheit sollen die Grundlagen multipler Regression - insbesondere die Beurteilung der Modellgüte und die Überprüfung eines Modells anhand der prinzipiellen Modellannahmen linearer Regression - an einem Beispiel behandelt werden.

tl;dr: Her mit dem Code!

5.1 Ouvertüre

Zentrales Ziel einer multiplen linearen Regression ist das Ableiten einer abhängigen (= zu erklärenden) aus mehreren unabhängigen (= erklärenden) Variablen. Um diesen Prozess näher kennenzulernen, wollen wir versuchen, die regionale Variabilität der Corona-Schutzimpfungsquoten in Österreich zu erklären.

Neben ersten individuellen Motivstudien (Universität Wien, 2021) liegen zur Erklärung der regionalen Variabilität von kommunalen und regionalen Impfquoten zum Stand Oktober 2021 erste Vermutungen zu möglichen Einflußfaktoren vor (zB Momentum Institut, 2021):

In dieser Einheit wollen wir einige dieser vermuteten Einflußfaktoren nutzen, um die unterschiedlichen Impfquoten (Stand 24.10.21) in den politischen Bezirken Österreichs zu erklären. Dazu greifen wir auf folgende Daten zurück:

Diese Datensätze wurden in einer Excel-Datei gesammelt, mittels Aggregation auf die einheitliche Bezugsebene der politischen Bezirke gebracht und im Tabellenblatt “ex” miteinander verknüpft.

🚩 Die Excel-Datei kann hier heruntergeladen werden 🚩

Neben diesen numerischen Daten nutzen wir in dieser Einheit auch folgenden Geodatensatz zur Visualisierung der Impfquoten:

Den Inhalt dieses ZIP-Archivs extrahieren wir in unserem “data” Ordner in den Unterordner “bez”.

👉 Anmerkung: Wir gehen in dieser Einheit von folgender Verzeichnisstruktur aus:

**Projektfolder**
| skript_1.R
| ...
| skript_n.R    
+-- data
|     bez
|       | geodatensatz_1.xyz
|       | ...
|       | geodatensatz_n.xyz
|     | datensatz_1.xyz
|     | ...
|     | datensatz_n.xyz
+-- output

5.2 Daten importieren

Über das readxl-Package laden wir zunächst die Daten aus dem Tabellenblatt “ex”:

library(readxl)     # Excel-Dateien lesen
library(tidyverse)  # https://www.tidyverse.org/packages/

daten <- read_excel("data/corona_bez_regression_v1.xlsx", sheet = "ex")

Damit erhalten wir:

head(daten)
## # A tibble: 6 x 10
##   bez_id bez_txt anteil_immun tote_100k bev_anteil_65pl~ bev_avg_alter anteil_noaut
##    <dbl> <chr>          <dbl>     <dbl>            <dbl>         <dbl>        <dbl>
## 1    101 Eisens~         69.7      73.9             19.5          43.5        15.7 
## 2    102 Rust(S~         70.4       0               25.9          47.2         7.9 
## 3    103 Eisens~         71.7      59.3             21.9          45.3         8.82
## 4    104 Güssing         70.1     233.              25.6          47.8         7.92
## 5    105 Jenner~         67.0     152.              23.3          47.4         6.35
## 6    106 Matter~         68.7      89.3             21.0          44.6         8.84
## # ... with 3 more variables: bildung_anteil_hochschule <dbl>,
## #   bildung_anteil_pflicht <dbl>, nrw19_anteil_fpoe <dbl>

Also viele Daten. Einen schnellen Einblick in die Struktur dieser Daten erhalten wir hiermit:

  str(daten)
## tibble [117 x 10] (S3: tbl_df/tbl/data.frame)
##  $ bez_id                   : num [1:117] 101 102 103 104 105 106 107 108 109 201 ...
##  $ bez_txt                  : chr [1:117] "Eisenstadt(Stadt)" "Rust(Stadt)" "Eisenstadt-Umgebung" "Güssing" ...
##  $ anteil_immun             : num [1:117] 69.7 70.3 71.7 70.1 67 ...
##  $ tote_100k                : num [1:117] 73.9 0 59.3 232.8 152 ...
##  $ bev_anteil_65plus        : num [1:117] 19.5 25.9 21.9 25.6 23.3 ...
##  $ bev_avg_alter            : num [1:117] 43.5 47.2 45.3 47.8 47.4 44.6 44.7 46.8 45.5 44.2 ...
##  $ anteil_noaut             : num [1:117] 15.67 7.9 8.82 7.92 6.35 ...
##  $ bildung_anteil_hochschule: num [1:117] 23.26 11.17 13.32 9.07 7.29 ...
##  $ bildung_anteil_pflicht   : num [1:117] 12.6 13.8 13.3 19.5 18.6 ...
##  $ nrw19_anteil_fpoe        : num [1:117] 13.7 19 15.6 16.5 21.2 ...

5.3 Daten validieren

Wir werfen eine Blick auf die Daten, v.a. auf fehlende Werte. Dazu gibt es mehrere Vorgehensweisen. Beispielsweise können über die ::colSums::-Funktion Spalten aufsummiert werden:

colSums(is.na(daten)) %>%
  knitr::kable()
x
bez_id 0
bez_txt 0
anteil_immun 0
tote_100k 23
bev_anteil_65plus 0
bev_avg_alter 0
anteil_noaut 0
bildung_anteil_hochschule 0
bildung_anteil_pflicht 0
nrw19_anteil_fpoe 0

Alternativ ermöglicht das dplyr-Package es uns mittels ::across:: beliebige Funktionen auf (ausgewählte - hier exemplarisch nur numerische) Spalten anzuwenden:

daten %>%
  summarise(across(where(is.numeric), ~ sum(is.na(.)))) %>%
  t() %>%
  knitr::kable()
bez_id 0
anteil_immun 0
tote_100k 23
bev_anteil_65plus 0
bev_avg_alter 0
anteil_noaut 0
bildung_anteil_hochschule 0
bildung_anteil_pflicht 0
nrw19_anteil_fpoe 0

Egal wie ermittelt, wir sehen, dass für die 23 Wiener Stadtbezirke keine Informationen zur Corona-bedingten Sterblichkeiten vorliegen. Da diese Informationen jedoch für die Stadt Wien gesamt vorliegen, können wir die Wiener Stadtbezirke (bez_id von 901 bis 923) bei unseren weiteren Analyse ausschließen.

5.3.1 Ein Exkurs: Die räumliche Variabilität der Impfquoten

Bevor wir das Regressionsmodell bilden, wollen wir noch einen Blick auf die räumliche Verteilung der Impfquoten legen. Aber wie?

Mittels des Simple Features (sf) Packages können gängige Geodatenformate (zB SHP-Dateien) in R gelesen und beispielsweise über das Thematic Maps (tmap) Package visualisiert werden:

library(sf)
library(tmap)
bez <- read_sf("data/bez/STATISTIK_AUSTRIA_POLBEZ_20210101.shp")
bez$id <- as.integer(bez$id)
tm_shape(bez) +
  tm_polygons()

Die PerfektionistInnen haben es bereits erkannt: Wien wird wieder als Summe seiner Stadtbezirke dargestellt. Ein Blick in den Geodatensatz verrät uns …

bez %>%
  filter(id >= 900)
## Simple feature collection with 24 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 611820.2 ymin: 473150.1 xmax: 641316 ymax: 496035.8
## Projected CRS: MGI / Austria Lambert
## # A tibble: 24 x 3
##       id name                                                              geometry
##  * <int> <chr>                                                   <MULTIPOLYGON [m]>
##  1   900 Wien(Stadt)           (((633092.9 476663.5, 633095.2 476662.3, 633097.3 4~
##  2   901 Wien  1.,Innere Stadt (((625604.8 482342.7, 625602.5 482342.8, 625572.3 4~
##  3   902 Wien  2.,Leopoldstadt (((626800.3 483669.1, 626762.9 483645.8, 626760.2 4~
##  4   903 Wien  3.,Landstraße   (((629078 480766.2, 629075.5 480761.7, 629073.1 480~
##  5   904 Wien  4.,Wieden       (((625990.5 480617.3, 625990 480617, 625988.9 48061~
##  6   905 Wien  5.,Margareten   (((625375.6 480306.9, 625358.7 480296.9, 625342.9 4~
##  7   906 Wien  6.,Mariahilf    (((623857.6 480935.6, 623782.8 480925.1, 623761 480~
##  8   907 Wien  7.,Neubau       (((623411.1 481718.7, 623411.1 481718.7, 623408.8 4~
##  9   908 Wien  8.,Josefstadt   (((624396.9 482942.4, 624335.8 482927.2, 624320.8 4~
## 10   909 Wien  9.,Alsergrund   (((624209.9 483872.5, 624200.4 483873, 624085.1 483~
## # ... with 14 more rows

… dass für den politischen Bezirk Wien bereits ein Polygon vorhanden ist, jedoch von den Polygonen der Stadtbezirke überlagert wird. Um diese Überlagerung zu vermeiden, entfernen (= filtern) wir die Wiener Stadtbezirke aus dem Geodatensatz “bez”:

bez.sel <- bez %>%
  filter(id <= 900)
tm_shape(bez.sel) +
  tm_polygons()

Sieht doch gleich besser aus 😎

Jetzt müssen wir nur noch unsere Attributdaten - die Impfquoten aus dem daten-Tibble - an die Geometrien der Bezirke hängen. Wie aus der Geoinformatik bekannt, verwenden wir dazu einen Join (konkret: die left_join Funktion des dplyr Packages):

joined_bez.sel <- left_join(bez.sel, daten,
                            by = c("id" = "bez_id"))
str(joined_bez.sel)
## sf [94 x 12] (S3: sf/tbl_df/tbl/data.frame)
##  $ id                       : num [1:94] 101 102 103 104 105 106 107 108 109 201 ...
##  $ name                     : chr [1:94] "Eisenstadt(Stadt)" "Rust(Stadt)" "Eisenstadt-Umgebung" "Güssing" ...
##  $ geometry                 :sfc_MULTIPOLYGON of length 94; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:1127, 1:2] 640982 640979 640912 640861 640842 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  $ bez_txt                  : chr [1:94] "Eisenstadt(Stadt)" "Rust(Stadt)" "Eisenstadt-Umgebung" "Güssing" ...
##  $ anteil_immun             : num [1:94] 69.7 70.3 71.7 70.1 67 ...
##  $ tote_100k                : num [1:94] 73.9 0 59.3 232.8 152 ...
##  $ bev_anteil_65plus        : num [1:94] 19.5 25.9 21.9 25.6 23.3 ...
##  $ bev_avg_alter            : num [1:94] 43.5 47.2 45.3 47.8 47.4 44.6 44.7 46.8 45.5 44.2 ...
##  $ anteil_noaut             : num [1:94] 15.67 7.9 8.82 7.92 6.35 ...
##  $ bildung_anteil_hochschule: num [1:94] 23.26 11.17 13.32 9.07 7.29 ...
##  $ bildung_anteil_pflicht   : num [1:94] 12.6 13.8 13.3 19.5 18.6 ...
##  $ nrw19_anteil_fpoe        : num [1:94] 13.7 19 15.6 16.5 21.2 ...
##  - 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:11] "id" "name" "bez_txt" "anteil_immun" ...

Damit können wir nun eine einfache Choroplethenkarte zur Impfquote erstellen:

tmap::qtm(joined_bez.sel, fill = "anteil_immun")

Und wer das gerne noch mit etwas mehr 🚀 🎉 möchte:

mymap <- tm_shape(joined_bez.sel) +
  tm_polygons("anteil_immun",
              title = "Anteil \nImmunisierte [%]",
              palette = "YlOrRd",
              legend.hist = TRUE) +
  tm_scale_bar(position = c("left", "bottom")) +
  tm_legend(outside = TRUE,
            legend.outside.size = 0.15,
            hist.width = 1,
            outer.margins = 0)
mymap

Was zeigen uns diese Karten?

  1. Man kann in R auch thematische Karten erstellen 😉
  2. Die Bezirke im Burgenland und dem nördlichen Niederösterreich weisen die höchsten Impfquoten auf. Weiters sehen wir, dass einige Bezirke in der Mur-Mürz-Furche sowie der Bezirk Schwaz in Tirol hervorstechen.

Zum Abschluß dieses Exkurses: Wie kann ich solche Karten speichern (zB in den Ordner “output”)?

tmap_save(mymap, filename = "output/impfquoten_bez_2021.png",
          units = "px", dpi = 300,
          width = 2000)

Moving on:

5.4 Die mathematischen Grundlagen linearer Modelle

Damit kommen wir zum eigentlich Kern: Wie können wir die Impfquoten (anteil_immun) aus den restlichen Variablen dieses Datensatzes ableiten?

Dazu werden wir ein lineares Modell nach diesem Vorbild schaffen:

\[\hat{Y} = b_0 + b_1x_1 + b_2x_2 + ... +b_Jx_J\]

\[\begin{aligned} & \hat{Y}~ ...~abhängige~Variable\\ & b_{0...j}~...~Koeffizienten\\ & x_{1...j}~...~unabhängige~Variablen\\ & J~...~Zahl~der~unabhängigen~Variablen \end{aligned}\]

Für dieses Modell wird die Summe der Differenzen (“Residuen”) zwischen die modellierten (“Fitted values”) und beobachteten Werten minimiert:

\[ \sum_{k=1}^{K}e_{k}^2 = \sum_{k=1}^{K}[y_k-(b_0 + b_1x_{1k} + b_2x_{2k} + ... +b_jx_{Jk})]^2 \rightarrow min \]

\[\begin{aligned} & e_k~...~Residuum~(k = 1, 2, ... K)\\ & y_k~...~Werte~abhängige~Variable~(k = 1, 2, ... K)\\ & b_0~...~Konstante~('Intercept')\\ & b_j~...~Koeffizienten~(j=1,2,...,J)\\ & x_{jk}~...~Werte~unabhängige~Variablen (j=1,2,...,J; k=1,2,...,K)\\ & J~...~Zahl~der~unabhängigen~Variablen\\ & K~...~Zahl~der~Beobachtungen \end{aligned}\]

Bevor es aber so weit ist:

5.5 Die gedankliche Modellbildung

Wir sollten uns zunächst darüber klar werden, mit welchen Variablen wir den Impfquoten der Bezirke erklären wollen. Dazu ein Blick auf die verfügbaren Daten:

cbind(colnames(daten))
##       [,1]                       
##  [1,] "bez_id"                   
##  [2,] "bez_txt"                  
##  [3,] "anteil_immun"             
##  [4,] "tote_100k"                
##  [5,] "bev_anteil_65plus"        
##  [6,] "bev_avg_alter"            
##  [7,] "anteil_noaut"             
##  [8,] "bildung_anteil_hochschule"
##  [9,] "bildung_anteil_pflicht"   
## [10,] "nrw19_anteil_fpoe"

In einer ersten Runde wollen wir zunächst folgende Variablen zur Erklärung nutzen:

  • Die Verstorbenen je 100.00 Einwohner (tote_100k)
  • Der Anteil der über 65-Jährigen (bev_anteil_65plus)
  • Der Anteil der Einwohner mit nicht-österreichischer Staatsbürgerschaft (anteil_noaut)
  • Der Anteil der HochschulabsolventInnen an der Gesamtbevölkerung (bildung_anteil_hochschule)
  • Der FPÖ-Stimmenanteil bei der Nationalratswahl 2019 (nrw19_anteil_fpoe)

Dieses reduzierte Datenset legen wir als eigenen Tibble ab und entfernen auch noch die Wiener Gemeindebezirke (bez_id > 900) daraus:

sel.daten <- daten %>%
  filter(bez_id <= 900) %>%
  select(bez_id, bez_txt, anteil_immun, tote_100k, bev_anteil_65plus,
         anteil_noaut, bildung_anteil_hochschule, nrw19_anteil_fpoe)

5.6 Ein Blick auf die gewählten Variablen

Nun wollen wir einen Blick auf die Wertverteilungen der ausgewählten Variablen werfen:

sel.daten %>%
  select(anteil_immun:nrw19_anteil_fpoe) %>%
  pivot_longer(cols = anteil_immun:nrw19_anteil_fpoe, names_to = "Variable", 
               values_to = "Messwert") %>%
  ggplot(., aes(x = Variable, y = Messwert)) +
  geom_boxplot() +
  geom_jitter(width=0.1,alpha=0.2, color="red") +
  theme(axis.text.x=element_text(angle = 45, hjust = 1))

Wir sehen, dass die Variablen durchaus unterschiedliche absolute Wertverteilungen aufweisen.

5.6.1 Standardisierung von Variablen

Um Variablen mit unterschiedlichen absoluten Wertverteilung besser vergleichbar zu machen, können diese standardisiert werden. Das einfachste Verfahren dazu ist die sgn. z-Transformation:

\[ z = \frac{x - \bar{x}}{s} \]

\[\begin{aligned} & z~...~\text{z-transformierter ('standardisierter') Wert}\\ & x~...~\text{beobachteter Wert}\\ & \bar{x}~...~\text{Mittelwert}\\ & s~...~\text{Standardabweichung} \end{aligned}\]

Man erhält dadurch Variablen mit einem Mittelwert von 0 und einer Standardabweichung von 1. Nutzt man diese standardisierten Variablen für die Regression, erhält man standardisierte Regressionskoeffizienten, die miteinander verglichen (genauer gesagt als partielle Korrelationskoeffizienten auslegt) werden können.

Die Standardisierung von Variablen können wir mit dem Befehl scale erzielen:

sel.daten.trans <- sel.daten %>%
  mutate(across(c("anteil_immun", "tote_100k", "bev_anteil_65plus",
         "anteil_noaut", "bildung_anteil_hochschule", "nrw19_anteil_fpoe"),scale))

Zur Kontrolle: Die Mittelwerte der so transformierten Variablen sollten bei 0 liegen:

# zur Kontrolle: Mittelwert = 0
sel.daten.trans %>%
  summarise(across(c("anteil_immun", "tote_100k", "bev_anteil_65plus",
         "anteil_noaut", "bildung_anteil_hochschule", "nrw19_anteil_fpoe"),mean)) %>%
  t() %>%
  knitr::kable()
anteil_immun 0
tote_100k 0
bev_anteil_65plus 0
anteil_noaut 0
bildung_anteil_hochschule 0
nrw19_anteil_fpoe 0

Ein kleiner Exkurs zur Veranschaulichung der Wirkung einer Standardisierung:

# Streuung der nicht-standardisierten Variablen abbilden
vis.data.1 <- sel.daten %>%
  select(anteil_immun:nrw19_anteil_fpoe) %>%
  pivot_longer(cols = anteil_immun:nrw19_anteil_fpoe, names_to = "Variable", 
               values_to = "Messwert")

p1 <- ggplot(vis.data.1, aes(x = Variable, y = Messwert)) +
  geom_boxplot() +
  labs(x = "Variablen\n", y = "\nbeobachtete Werte") +
  theme_gray(base_size = 16) +
  coord_flip()

# Streuung der standardisierten Variablen darstellen
vis.data.2 <- sel.daten.trans %>%
  select(anteil_immun:nrw19_anteil_fpoe) %>%
  pivot_longer(cols = anteil_immun:nrw19_anteil_fpoe, names_to = "Variable", 
               values_to = "Messwert")

p2 <- ggplot(vis.data.2, aes(x = Variable, y = Messwert)) +
  geom_boxplot() +
  labs(y = "\nz-transf. Werte") +
  theme_gray(base_size = 16) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank()) +
  coord_flip()

# beide Plots kombinieren damit der Vergleich leichter fällt mit dem grid_extra package
library(gridExtra)
grid.arrange(p1, p2, nrow = 1, widths = c(2, 1))

5.6.2 Beziehungen der Variablen zueinander

Bevor wir nun das Regressionsmodell befüllen, werfen wir noch einen Blick auf die Korrelationen zwischen den erklärenden Variablen.

Warum?
Eine Annahme linearer Regressionsmodelle ist, dass zwischen den erklärenden Variablen keine Abhängigkeiten bestehen sollen. Das können wir numerisch anhand der Korrelationskoeffizienten zwischen den erklärenden Variablen überprüfen. In der Old-School-Variante erreicht man dies beispielsweise anhand des rcorr Befehls:

# old school: numerisch
library(Hmisc)
rcorr(as.matrix(sel.daten[,4:8]))
##                           tote_100k bev_anteil_65plus anteil_noaut
## tote_100k                      1.00              0.28        -0.23
## bev_anteil_65plus              0.28              1.00        -0.44
## anteil_noaut                  -0.23             -0.44         1.00
## bildung_anteil_hochschule     -0.27             -0.28         0.64
## nrw19_anteil_fpoe              0.47              0.23        -0.38
##                           bildung_anteil_hochschule nrw19_anteil_fpoe
## tote_100k                                     -0.27              0.47
## bev_anteil_65plus                             -0.28              0.23
## anteil_noaut                                   0.64             -0.38
## bildung_anteil_hochschule                      1.00             -0.51
## nrw19_anteil_fpoe                             -0.51              1.00
## 
## n= 94 
## 
## 
## P
##                           tote_100k bev_anteil_65plus anteil_noaut
## tote_100k                           0.0064            0.0263      
## bev_anteil_65plus         0.0064                      0.0000      
## anteil_noaut              0.0263    0.0000                        
## bildung_anteil_hochschule 0.0086    0.0068            0.0000      
## nrw19_anteil_fpoe         0.0000    0.0231            0.0002      
##                           bildung_anteil_hochschule nrw19_anteil_fpoe
## tote_100k                 0.0086                    0.0000           
## bev_anteil_65plus         0.0068                    0.0231           
## anteil_noaut              0.0000                    0.0002           
## bildung_anteil_hochschule                           0.0000           
## nrw19_anteil_fpoe         0.0000

Etwas hübscher gelingt dies mittels des GGally-Packages:

# etwas hübscher: graphisch mittels GGally
library(GGally)
ggpairs(sel.daten, columns = 4:8)

👉 Eine Daumenregel zur Interpretation:
Korrelationen über 0,8 sind meist problematisch, da sie die Präzision der Koeffizienten negativ beeinflussen würden. In unserem Fall können wir also von keinen problematischen Korrelationen zwischen den erklärenden Variablen ausgehen.

Und damit kommen wir (endlich) zur:

5.7 Die Modellbildung

Wir bilden zunächst ein Modell mit den fünf ausgewählten erklärenden Variablen (vgl. “Die gedankliche Modellbildung”):

# options(scipen = 999)
# options(scipen = 0)

lm.v1.trans <- lm(anteil_immun ~ tote_100k + bev_anteil_65plus + anteil_noaut + bildung_anteil_hochschule + nrw19_anteil_fpoe,
                  data = sel.daten.trans)
summary(lm.v1.trans)
## 
## Call:
## lm(formula = anteil_immun ~ tote_100k + bev_anteil_65plus + anteil_noaut + 
##     bildung_anteil_hochschule + nrw19_anteil_fpoe, data = sel.daten.trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.73891 -0.58091  0.03349  0.50599  1.78176 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -4.931e-16  8.662e-02   0.000   1.0000    
## tote_100k                 -2.179e-01  1.006e-01  -2.166   0.0330 *  
## bev_anteil_65plus          4.550e-01  9.889e-02   4.602  1.4e-05 ***
## anteil_noaut              -2.656e-01  1.207e-01  -2.200   0.0304 *  
## bildung_anteil_hochschule  3.924e-01  1.218e-01   3.220   0.0018 ** 
## nrw19_anteil_fpoe         -6.361e-02  1.106e-01  -0.575   0.5667    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8398 on 88 degrees of freedom
## Multiple R-squared:  0.3327, Adjusted R-squared:  0.2948 
## F-statistic: 8.775 on 5 and 88 DF,  p-value: 8.746e-07

Wir sehen, dass unser erster Modellversuch ein adjustiertes R² von 0.2948 also rund 30% aufweist. Das ist zunächst nicht allzu berauschend, bedenkt man dass 70% der Varianz der Impfquoten durch dieses Modell nicht erklärt werden.

Ein Blick auf die Regressionskoeffizienten (Spalte “Estimates”) zeigt uns, dass

  • Der Anteil der über 65-Jährigen und der Anteil der HochschulabsolventInnen positiv,

  • und die Corona-bezogene Mortalität und der Anteil der Nicht-ÖsterreicherInnen negativ

auf die Impfquote einwirken.

Ein Blick auf die Signifikanz der erklärenden Variablen verrät uns darüber hinaus, dass der FPÖ-Stimmanteil nicht signifikant zur Erklärung der Impfquote beiträgt.

5.7.1 Ein alternatives Modell

Ausgehend von diesem ersten Modell wollen wir ein verbessertes Modell bilden. Dazu nehmen wir zwei Änderungen vor:

  1. Wir entfernen den nicht signifikanten FPÖ-Stimmanteil als erklärenden Variable:
lm.v1.trans <- lm(anteil_immun ~ tote_100k + bev_anteil_65plus + anteil_noaut + bildung_anteil_hochschule,
                  data = sel.daten.trans)
summary(lm.v1.trans)
## 
## Call:
## lm(formula = anteil_immun ~ tote_100k + bev_anteil_65plus + anteil_noaut + 
##     bildung_anteil_hochschule, data = sel.daten.trans)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8175 -0.5766 -0.0157  0.5337  1.8124 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -4.840e-16  8.629e-02   0.000 1.000000    
## tote_100k                 -2.404e-01  9.238e-02  -2.602 0.010854 *  
## bev_anteil_65plus          4.546e-01  9.852e-02   4.615 1.31e-05 ***
## anteil_noaut              -2.619e-01  1.201e-01  -2.181 0.031805 *  
## bildung_anteil_hochschule  4.162e-01  1.141e-01   3.647 0.000447 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8366 on 89 degrees of freedom
## Multiple R-squared:  0.3302, Adjusted R-squared:  0.3001 
## F-statistic: 10.97 on 4 and 89 DF,  p-value: 2.823e-07

Außer einer marginalen Verbesserung der adjustierten R²-Wertes können wir leichte Veränderungen der Koeffizienten beobachten.

  1. Als Experiment wollen wir auch noch den Bezirk Schwaz in Tirol (bez_id = 709) aus dem Sample entfernen.
    Warum?
    Dieser Bezirk wurde vor allen anderen Bezirken im März und April 2021 einer gezielten Impfkampagne unterzogen. Was uns zu folgender These verleiten kann: Durch diese Impfkampagne wurde die messbare Impfquote nach oben verschoben.
# Schwaz aus Modellbildung ausschließen
lmv2noSchwaz <- sel.daten.trans %>%
  filter(bez_id != 709)

lm.v2.trans <- lm(anteil_immun ~ tote_100k + bev_anteil_65plus + anteil_noaut + bildung_anteil_hochschule, 
                  data = lmv2noSchwaz)
summary(lm.v2.trans)
## 
## Call:
## lm(formula = anteil_immun ~ tote_100k + bev_anteil_65plus + anteil_noaut + 
##     bildung_anteil_hochschule, data = lmv2noSchwaz)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.80796 -0.57097 -0.01173  0.54936  1.60235 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -0.02029    0.08481  -0.239 0.811431    
## tote_100k                 -0.25885    0.09065  -2.855 0.005361 ** 
## bev_anteil_65plus          0.48591    0.09726   4.996 2.95e-06 ***
## anteil_noaut              -0.27791    0.11756  -2.364 0.020279 *  
## bildung_anteil_hochschule  0.44339    0.11219   3.952 0.000156 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8177 on 88 degrees of freedom
## Multiple R-squared:  0.3638, Adjusted R-squared:  0.3349 
## F-statistic: 12.58 on 4 and 88 DF,  p-value: 3.876e-08

Die Steigerung des adjustierten R²-Wertes auf rund 33% kann als ein erster Beleg für diese These gewertet werden.

5.8 Überprüfung der Modellannahmen

Nun gilt es für das gefundene Modell nachzuweisen, dass dieses die grundsätzlichen Modellannahmen linearer Regression nicht verletzt. Diese Annahmen sind:

  1. Modell ist ausreichend spezifiziert

    • Linearer Zusammenhang gegeben
    • Relevante erklärende Variablen eingebunden
    • Anzahl Koeffizienten < Anzahl Beobachtungen

  2. Erklärende Variablen = unkorreliert (Multikollinearität)

  3. Erwartungswert der Residuen = 0

  4. Residuen = normalverteilt

  5. Residuen haben konstante Varianz (Homoskedastiziät)

  6. Keine Korrelation zwischen erklärenden Variablen und Residuen

  7. Residuen sind unkorreliert zu Beobachtungen (Autokorrelation)

5.8.1 Prüfung des linearen Zusammenhangs

Der einfachste Weg hierzu ist ein graphischer Vergleich des Zusammenhangs zwischen der zu erklärenden und den erklärenden Variablen:

ggpairs(sel.daten, c(3:8))

Der ersten Spalte dieser Matrix entnehmen wird, dass keine der erklärenden Variablen einen klar erkennbaren exponentiellen oder einen sonstigen nicht-linearen Zusammenhang zur Impfquote aufweist.

Ein weitere oftmals genutzte Möglichkeit bietet der sgn. “Residuals vs. Fitted”-Plot:

plot(lm.v2.trans, 1, labels.id = sel.daten$bez_txt)

Bei der Gegenüberstellung der durch das Modell vorhergesagten Werte (“Fitted values”) und den dabei erzeugten Abweichungen von den beobachteten Werten (“Residuals”) sollte

  1. sich kein klares Muster in der Punktverteilung erkennen lassen und
  2. die rote Regressionslinie im Idealfall der Null-Linie entsprechen.

In unserem Fall sind beide Bedingungen noch hinreichend gut erfüllt. Ein linearer Zusammenhang kann damit als hwst. gegeben angesehen werden.

Sollte sich im “Residuals vs. Fitted”-Plot ein klar nicht-linearer Zusammenhang zeigen, kann eine Transformation der Ausgangsdaten (beispielsweise durch Logarithmieren) ggf. einen linearen Zusammenhang bewirken.

5.8.2 Prüfung der Unabhängigkeit der erklärenden Variablen (“Multikollinearität”)

Wie bereits zu Beginn erwähnt, sollen die erklärenden Variablen voneinander unabhängig sein. Dies kann am einfachsten über die paarweise Korrelationen abgeklärt werden:

ggpairs(sel.daten, columns = 4:8)

Anhand der bereits bekannten Daumenregel “Keine Korrelation größer als 0,8” können wir keine problematischen Korrelationen feststellen.

Generell kann das Auffüllen von Modellen mit vielen korrelierten Variablen zu einem “Overfitting” führen. Dabei kann das overfittete Modell zwar die beobachtete Stichprobe an Werten gut abbilden, bei einer Erweiterung der Stichprobe besteht aber die Gefahr einer sinkenden Erklärungskraft.

Trotz aller Freude an multivariaten Verfahren sollte ein “lean model” - also ein Modell mit sparsamen Variableneinsatz - stets das Ziel der Analyse sein. Was solche “lean models” mit Rasiermessern zu tun haben, kann ::hier:: nachgelesen werden 😉

5.8.3 Prüfung der Normalverteilung der Residuen

Da die Residuen ja Abbild einer zufälligen Streuung sein sollten, müssten sie einer Normalverteilung folgen. Auch hier bietet R einen passenden Standard-Plot für lineare Modelle: Einen sgn. “Q-Q Plot”

plot(lm.v2.trans, 2, labels.id = sel.daten$bez_txt)

In diesem Plot werden die Residuen nach ihren Quantilszugehörigkeiten angeordnet. Die dabei entstehenden Punkte sollte bei vorliegender Normalverteilung auf der Gleichverteilungsgeraden zu liegen kommen. Abweichungen von dieser Geraden deuten auf eine Abweichung von der Normalverteilung hin.

In unserem Fall schlagen einige wenige Ausreißer nach oben und unten aus. In Summe entspricht die Verteilung jedoch hinreichend einer Normalverteilung. Dafür spricht auch ein hinreichend nahe Null liegender Mittelwert der Residuen:

mean(residuals(lm.v2.trans))
## [1] 2.143893e-17

5.8.4 Prüfung der Konstanz der Varianz der Residuen (“Homoskedastizität”)

Eine weitere Annahme linearer Regressionsmodelle besagt, dass die Residuen über die unabhängigen Variablen (und damit auch über die abhängige Variable) hinweg gleichverteilt sein sollen. Diese Annahme sichert die Homogenität in den Varianzen (“Homoskedastizität”) und damit letztlich, dass Messfehler bei einzelnen unabhängigen Variablen nicht Eingang in die Modellierung finden und zur Verzerrung der Regressionskoeffizienten beitragen.

Auch hierzu bietet R einen passenden Standard-Plot für lineare Modelle:

plot(lm.v2.trans, 3, labels.id = sel.daten$bez_txt)

Im Idealfall zeigt die Verteilung der Varianzen (Wurzel der “Standardized residuals”) über die prognostizierten Werte (“Fitted values”) eine annähernd horizontale Verteilungslinie mit gleichverteilten Punkten.

Aufgrund zweier Ausreißer im oberen Wertspektrum der fitted values erfüllt unser Modell dieses Kriterium nur knapp. Bei noch weiter streuenden Varianzen bieten sich üblicherweise zwei Möglichkeiten an:

  1. Über eine Transformation der Ausgangsdaten (beispielsweise ein Logarithmieren) kann eine homogenere Verteilung der Varianzen erzielt werden. In unserem Fall ist dies jedoch nicht (direkt) möglich, da auch negative Werte der zu erklärenden Variable möglich sind.
  2. Oftmals kann durch eine Bereinigung der Daten um Ausreißer oder eine Segmentierung der Daten in zwei oder mehr Subgruppen eine homogenere Verteilung der Varianzen erzielt werden.

Exkurs: Ausreißern auf der Spur

So Ausreißer in unserem Modell keine Rolle spielen, müsste in Modell mit einer getrimmten Wertverteilung der Ausgangsdaten (zB 95% oder 97,5%) zu ähnlichen Koeffizienten und Erklärungsgehalt führen. In unserem Fall trimmen wir die Variable mit den meisten Ausreißern (vgl. Ein Blick auf die gewählten Variablen) - den Anteil der HochschulabsolventInnen um die obersten 2,5% der Werte:

trimmed.data <- lmv2noSchwaz %>%
  filter(bildung_anteil_hochschule < quantile(bildung_anteil_hochschule, 0.975))

# Vergleich der Streuungen vor und nach der Trimmung
vis.min <- min(lmv2noSchwaz$bildung_anteil_hochschule)
vis.max <- max(lmv2noSchwaz$bildung_anteil_hochschule)

p3 <- ggplot(lmv2noSchwaz, aes(x = bildung_anteil_hochschule, y = "ungetrimmt")) +
  geom_boxplot() +
  labs(x = "Ausprägungen bildung_anteil_hochschule\n", y = "") +
  xlim(vis.min, vis.max) +
  theme_gray() +
  coord_flip()
p4 <- ggplot(trimmed.data, aes(x = bildung_anteil_hochschule, y = "getrimmt (97,5%)")) +
  geom_boxplot() +
  xlim(vis.min, vis.max) +
  labs(x = "", y = "") +
  coord_flip()
# beide Plots kombinieren damit der Vergleich leichter fällt mit dem grid_extra package
grid.arrange(p3, p4, nrow = 1)

# Regressionsmodell V2 mit getrimmten Ausgangswerten berechnen
lm.v2.trans.trimmed <- lm(anteil_immun ~ tote_100k + bev_anteil_65plus + anteil_noaut + bildung_anteil_hochschule, 
                  data = trimmed.data)
summary(lm.v2.trans.trimmed)
## 
## Call:
## lm(formula = anteil_immun ~ tote_100k + bev_anteil_65plus + anteil_noaut + 
##     bildung_anteil_hochschule, data = trimmed.data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.79863 -0.55805  0.02059  0.56046  1.60745 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -0.01536    0.08823  -0.174 0.862210    
## tote_100k                 -0.26033    0.09338  -2.788 0.006542 ** 
## bev_anteil_65plus          0.48879    0.09919   4.928 4.06e-06 ***
## anteil_noaut              -0.28432    0.12152  -2.340 0.021644 *  
## bildung_anteil_hochschule  0.47051    0.13166   3.574 0.000583 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.828 on 85 degrees of freedom
## Multiple R-squared:  0.369,  Adjusted R-squared:  0.3393 
## F-statistic: 12.42 on 4 and 85 DF,  p-value: 5.306e-08
summary(lm.v2.trans)
## 
## Call:
## lm(formula = anteil_immun ~ tote_100k + bev_anteil_65plus + anteil_noaut + 
##     bildung_anteil_hochschule, data = lmv2noSchwaz)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.80796 -0.57097 -0.01173  0.54936  1.60235 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -0.02029    0.08481  -0.239 0.811431    
## tote_100k                 -0.25885    0.09065  -2.855 0.005361 ** 
## bev_anteil_65plus          0.48591    0.09726   4.996 2.95e-06 ***
## anteil_noaut              -0.27791    0.11756  -2.364 0.020279 *  
## bildung_anteil_hochschule  0.44339    0.11219   3.952 0.000156 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8177 on 88 degrees of freedom
## Multiple R-squared:  0.3638, Adjusted R-squared:  0.3349 
## F-statistic: 12.58 on 4 and 88 DF,  p-value: 3.876e-08

Im direkten Vergleich der getrimmten mit der ungetrimmten Modellvariante sehen wir nur geringe Unterschiede bei den Koeffizienten sowie beim Erklärungsgehalt (R²-Wert). Wir können daraus schließen, dass Ausreißer unsere Modellparameter nur gering verzerren.

Der Wirkung von Ausreißern - also extremen Beobachtungen - kann auch mittels eigener Kennzahlen wie der Cookschen Distanz oder der Leverage nachgegangen werden:

Anhand der Cookschen Distanz (“Cook’s distance”) kann auf die Veränderungen im Modellergebnis durch das Weglassen der n-ten Beobachtung geschlossen werden. Je größer dieser Wert ausfällt, umso mehr Impact hat die Beobachtung auf die Modellierung. Als Daumenregel sollte dieser Wert deutlich unter 1, idealerweise auch klar unter 0,5 liegen.

plot(lm.v2.trans, 4, labels.id = sel.daten$bez_txt)

In unserem Fall trifft dies auf alle Beobachtungen zu.

Die Abbildung “Residuals vs Leverage” vereinigt zwei interessante Indikatoren:

  • Die standardisierten Residuen: Berechnet als die durch den Standardfehler dividierten Residuen. Als Daumenregel sollten diese nicht größer als das 3 (James et al. 2004) ausfallen.

  • Die “Leverage”: Die auch als “hat values” bezeichneten Werte messen den Einfluss einzelner Werte auf die Regressionskoeffizienten. Als Daumenregel sollten diese Werte kleiner als 2(p + 1)/n ausfallen. Dabei steht

    • p für die Anzahl der erklärenden Variablen und
    • n für die Anzahl der Beobachtungen
plot(lm.v2.trans, 5, labels.id = sel.daten$bez_txt)

In unserem Fall liegt diese Schwelle also bei 2*(4 + 1)/93 = 0,11 welche nur durch wenige Ausreißer übertroffen wird. In der Zusammenschau aus standardisierten Residuen und der Leverage können jedoch starke Verzerrungen des Regressionsmodells ausgeschlossen werden.

5.8.5 Prüfung auf Autokorrelation

Eine weitere Annahme linearer Modelle ist die Unabhängigkeit der Residuen von der Reihenfolge der Beobachtungen. Dieser Effekt ist vor allem bei Zeitreihendaten zu beachten. Ist diese Unabhängigkeit nicht gegeben, spricht man von “Autokorrelation”, welche meist auf systematische Messfehler (z.B. wachsende Unaufmerksamkeit) zurückgeführt werden kann. Autokorrelation führt in linearen Modellen zu einer Verzerrung der Regressionskoeffizienten.

Zur Überprüfung werden die Residuen über die Beobachtungen aufgetragen:

ggplot(lmv2noSchwaz, aes(x=(1:length(bez_id)), y=residuals(lm.v2.trans))) +
  geom_point() +
  geom_smooth(method='lm', formula= y~x, se=F, color="red") +
  labs(x = "\nBeobachtungen", y = "Residuen\n") +
  theme_gray(base_size = 18)

Wie zu erwarten (wir nutzen keine Zeitreihendaten) sind die Residuen über die Beobachtungen gleichverteilt und die (rote) Regressionsgerade der Residuen entspricht der Null-Linie. Es liegen also keine Indizien für Autokorrelation vor.

🏆 Nun wissen wir, …

  • daß die Bildung von Regressionsmodellen meist ein iterativer Prozess ist.
  • daß man in diesem Prozess viel Informationen die kausale Abhängigkeiten zwischen Variablen überprüfen kann.
  • daß ein Regressionsmodell zwar in sich konsistent, aber aufgrund der gewählten Variablen nur einen geringen Erklärungswert besitzen kann.

Aber Moment 🤓, da war doch noch etwas …