6 Multiple Regression II: Qualitative Variablen einbinden

📢 Zielsetzung dieser Einheit

In dieser Einheit wird die Nutzung qualitativer - also nominal oder ordinal skalierter - Variablen in Regressionsmodellen anhand des Beispiel zur Erklärung der bezirksweisen Impfquoten aus Einheit 5 behandelt.

tl;dr: Her mit dem Code!

6.1 Ouvertüre

Wie in Einheit 5 nutzen wir erneut den Datensatz zu den bezirksweiten Impfquoten (Stand 24.10.21). In diesem Datensatz wurden mehrere Teildatensätze miteinander verknüpft. Eine genaue Beschreibung dieser Teildatensätze findet sich in Kapitel 5.1. Diese Teildatensä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 🚩

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

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

6.2 Gedankliche Modellbildung

In der vorigen Einheit konnten wir in Kapitel 5.3.1 einen ersten Überblick auf die räumliche Verteilung der bezirksweisen Impfquoten gewinnen:

Anteil der Immunisierten in den österreichischen Bezirken (Stand: Oktober 2021)

Dabei sehen wir, dass vor allem die Bundesländer Salzburg, Oberösterreich und Kärnten mit niedrigen und das Burgenland sowie Teile Niederösterreichs und der Steiermark mit hohe Impfquoten aufweisen. Wir wollen daher überprüfen, ob die Hinzunahme der “Bundeslandzugehörigkeit” unser Regressionsmodell (vgl. Kapitel 5)

a) hinsichtlich der Varianzaufklärung verbessert;

b) und zu Veränderungen in den Regressionsparametern führt.

Folgender Workflow bietet sich dazu an

Wir importieren die Ausgangsdaten;
ermitteln die Bundeslandzugehörigkeit der Bezirke;
ergänzen unser Regressionsmodell um diese Zugehörigkeit;
und prüfen, ob wir damit die Annahmen linearer Modelle erfüllen.

6.3 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 den bekannten Datensatz zu den bezirksweisen Impfquoten zum Stand Oktober 2021:

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>

6.4 Die Bundeslandzugehörigkeit der Bezirke ermitteln

Zu welchem Bundesland ein Bezirk gehört, können wir in zwei Schritten ermitteln:

Zunächst extrahieren wir aus der Bezirkskennzahl “bez_id” das jeweilige Bundesland. Wir erinnern uns: Die erste Stelle der Bezirkskennzahl steht für das jeweilige Bundesland.

daten <- daten %>%
  mutate(bld = floor(bez_id/100))

Um das nominale Skalenniveau dieser Bundeslandzugehörigkeit auch noch korrekt abzubilden, überführen wir die Variable “bld” in einen Faktor. Um dabei zu wissen, welche Zahl für welches Bundesland steht: Die Verwaltungseinheiten werden dazu in Österreich immer alphabetisch aufsteigend sortiert und nummeriert.

daten <- daten %>%
  mutate(bld = factor(bld, labels = c("Bgld.", "Ktn.", "NÖ", "OÖ",
                                      "Sbg.", "Stmk.", "T", "Vbg.", "W")))
head(daten$bld)

## [1] Bgld. Bgld. Bgld. Bgld. Bgld. Bgld.
## Levels: Bgld. Ktn. NÖ OÖ Sbg. Stmk. T Vbg. W

Damit bleibt uns nur mehr, den Datensatz für die weitere Analyse etwas auszudünnen - also nicht benötigte Variablen und Records (vgl. Kapitel 5.7.1) zu entfernen:

sel_daten <- daten %>%
  filter(bez_id <= 900 & bez_id != 709) %>%
  select(bez_id, bez_txt, bld, anteil_immun, tote_100k, bev_anteil_65plus,
         anteil_noaut, bildung_anteil_hochschule, nrw19_anteil_fpoe)

6.5 Qualitative Variablen in Regressionsmodellen nutzen

Prinzipiell können qualitative - also nominal und ordinal skalierte - Variablen in Regressionsmodellen berücksichtigt werden. Man sollte aber folgende zwei Punkte nicht aus den Augen verlieren:

Regressionsanalyse misst den Einfluss einer oder mehrerer quantitativer unabhängiger auf eine quantitative abhängige Variable (vgl. Kapitel 1 in @Backhaus2018). Mit ihr kann beispielsweise der Einfluss der Körpergröße auf das Gewicht einer Person untersucht werden.
Qualitative Variablen können als sgn. Dummy-Variablen berücksichtigt werden, sollten jedoch nicht den überwiegenden Teil der unabhängigen Variablen stellen.
Varianzanalyse misst den Einfluss qualitativer unabhängiger Variablen auf eine quantitative abhängige Variable (vgl. Kapitel 3 in @Backhaus2018). Mit ihr kann beispielsweise der Einfluss des Geschlechts auf Einkommen einer Person untersucht werden. Analog zur Regressionsanalyse können dabei auch quantitative Variablen (als sgn. “Kovariaten”) berücksichtigt werden; sie sollten aber nicht den überwiegenden Teil der unabhängigen Variablen stellen.

Wir entscheiden uns, das bestehende Regressionsmodell zu erweitern und müssen dafür die Bundeslandzugehörigkeit in sgn. Dummy-Codes überführen. Dazu wird jede Merkmalsausprägung D einer qualitativen Variable in einer logischen Variable (binär) abgebildet:

\[ D=\begin{cases} 1 & \text{Merkmal vorhanden} \\ 0 & \text{Merkmal nicht vorhanden} \end{cases} \]

Glücklicherweise überführt R Faktoren in Regressionsmodellen automatisch in Dummy-Codes.

📚 Exkurs: Dummy-Codes from the inside

Um diese Dummy-Codes etwas besser zu verstehen, werfen wir zunächst einen Blick auf die Levels des Faktors Bundesland:

levels(sel_daten$bld)

## [1] "Bgld." "Ktn."  "NÖ"    "OÖ"    "Sbg."  "Stmk." "T"     "Vbg."  "W"

Diese neun Merkmalsausprägungen werden von R automatisch in n-1 (sprich: Enminuseeeeins mit n = Anzahl der Faktor-Levels), also 8 Dummy-Variablen überführt:

model.matrix(~ bld, sel_daten) %>%  # Dummy-Coding liefert ...
  .[,-1] %>%        # ... eine Matrix ...
  as_tibble() %>%   # ... die ein einen Tibble überführt wird ...
  bind_cols(sel_daten[c("bez_txt", "bld")], .) %>%  # ... der noch einleitende Spalten bekommt
  head()

## # A tibble: 6 x 10
##   bez_txt            bld   bldKtn. bldNÖ bldOÖ bldSbg. bldStmk.  bldT bldVbg.  bldW
##   <chr>              <fct>   <dbl> <dbl> <dbl>   <dbl>    <dbl> <dbl>   <dbl> <dbl>
## 1 Eisenstadt(Stadt)  Bgld.       0     0     0       0        0     0       0     0
## 2 Rust(Stadt)        Bgld.       0     0     0       0        0     0       0     0
## 3 Eisenstadt-Umgebu~ Bgld.       0     0     0       0        0     0       0     0
## 4 Güssing            Bgld.       0     0     0       0        0     0       0     0
## 5 Jennersdorf        Bgld.       0     0     0       0        0     0       0     0
## 6 Mattersburg        Bgld.       0     0     0       0        0     0       0     0

Am Bezirk Eisenstadt (Stadt) sehen wir, dass alle acht vergebenen Dummy-Codes 0 entsprechen.

🤔 Kann das stimmen?

Jep, da wir ja nur n-1 Dummy-Codes nutzen um die Bundeslandzugehörigkeiten der Bezirke abzubilden. Im Fall von Eisenstadt (Stadt) heißt das, dass sich R dazu entschieden hat, das erste Level des Faktors “bld” (= “Bgld.”) als Linearkombination der restlichen acht Levels dieses Faktors zu interpretieren. Die Logik dahinter: Wenn ein Bezirk nicht zu den acht Bundesländer Kärnten bis Wien gehört, muss er zwangsweise ein burgenländischer Bezirk sein.
Inhaltlich spielt es keine Rolle, welches Level eines Faktors als Linearkombination interpretiert wird. Jedoch: Bei der Interpretation der Regressionsparameter ist es wichtig zu wissen, welches Level eines Faktors als Linearkombination abgebildet wurde. Mehr dazu gleich …

6.6 Der Einfluß des Bundeslandes auf die Impfquote

Kommen wir nun zur eigentlichen Modellbildung. Dazu werden wir folgende Schritte durchlaufen:

Wir standardisieren die metrischen Variablen in unserem Datensatz, um die Vergleichbarkeit der Regressionsparameter sicherzustellen;
Wir reproduzieren das in Einheit 5 optimierte Regressionsmodell (lm0), um Vergleichswerte zur Varianzaufklärung und der Regressionsparameter zu erhalten ;
Wir bilden ein neues Regressionsmodell (lm1), das die Bundeslandzugehörigkeit der Bezirke berücksichtigt;
optimieren dieses (lm2 & lm3);
und vergleichen dieses neue Modell mit dem im Schritt 2 ermittelten Referenzmodell (aus Einheit 5).

6.6.1 Standardisieren der metrischen Variablen

Analog zu Kapitel 5.6.1 nutzen wir hierzu eine Z-Transformation:

sel_daten_trans <- sel_daten %>%
  mutate(across(c("anteil_immun", "tote_100k", "bev_anteil_65plus",
         "anteil_noaut", "bildung_anteil_hochschule", "nrw19_anteil_fpoe"),scale))

6.6.2 Reproduktion unseres Referenzmodells aus Einheit 5

Dazu greifen wir auf das in Kapitel 5.7.1 formulierte Modell zurück:

lm0 <- lm(anteil_immun ~ tote_100k + bev_anteil_65plus + anteil_noaut 
          + bildung_anteil_hochschule, data = sel_daten_trans)
summary(lm0)

## 
## 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.8033 -0.5695 -0.0117  0.5479  1.5982 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -5.745e-16  8.457e-02   0.000 1.000000    
## tote_100k                 -2.590e-01  9.070e-02  -2.855 0.005361 ** 
## bev_anteil_65plus          4.832e-01  9.671e-02   4.996 2.95e-06 ***
## anteil_noaut              -2.784e-01  1.178e-01  -2.364 0.020279 *  
## bildung_anteil_hochschule  4.436e-01  1.122e-01   3.952 0.000156 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8155 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

6.6.3 Erweiterung des Modells um die Bundeslandzugehörigkeit der Bezirke

Da R die Dummy-Codierung des Faktors bld automatisch vornimmt, können wir diesen einfach zur Modellgleichung hinzufügen:

lm1 <- lm(anteil_immun ~ tote_100k + bev_anteil_65plus + anteil_noaut 
          + bildung_anteil_hochschule + bld, data = sel_daten_trans)
summary(lm1)

## 
## Call:
## lm(formula = anteil_immun ~ tote_100k + bev_anteil_65plus + anteil_noaut + 
##     bildung_anteil_hochschule + bld, data = sel_daten_trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.61748 -0.26840 -0.02744  0.29904  1.44567 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                1.25425    0.22666   5.534 3.85e-07 ***
## tote_100k                 -0.15211    0.08420  -1.807 0.074580 .  
## bev_anteil_65plus          0.26523    0.09252   2.867 0.005298 ** 
## anteil_noaut              -0.18189    0.10007  -1.818 0.072859 .  
## bildung_anteil_hochschule  0.34643    0.09362   3.700 0.000394 ***
## bldKtn.                   -2.03304    0.29124  -6.981 7.73e-10 ***
## bldNÖ                     -0.91696    0.24846  -3.691 0.000407 ***
## bldOÖ                     -1.92249    0.29944  -6.420 9.00e-09 ***
## bldSbg.                   -1.65277    0.36251  -4.559 1.82e-05 ***
## bldStmk.                  -0.93820    0.29806  -3.148 0.002312 ** 
## bldT                      -1.56069    0.34817  -4.483 2.43e-05 ***
## bldVbg.                   -0.97117    0.43766  -2.219 0.029318 *  
## bldW                      -1.21916    0.71975  -1.694 0.094182 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6085 on 80 degrees of freedom
## Multiple R-squared:  0.678,  Adjusted R-squared:  0.6297 
## F-statistic: 14.04 on 12 and 80 DF,  p-value: 3.822e-15

🤔 Was sehen wir hier?

Vergleichen wir zunächst einmal die Varianzaufklärung unserer beiden Modelle lm0 und lm1:

cbind(c("lm0", "lm1"), c(summary(lm0)$adj.r.squared, summary(lm1)$adj.r.squared))

##      [,1]  [,2]               
## [1,] "lm0" "0.334897082123657"
## [2,] "lm1" "0.629715621183297"

Die Hinzunahme der Bundeslandzugehörigkeit hat also zu einer Verdoppelung der Varianzaufklärung geführt. Mit knapp 63% Varianzaufklärung liegt unser Modell nicht schlecht (klar über 50%), wenn auch noch deutlich Luft nach oben gegeben ist.

Kommen wir zu den Regressionskoeffizienten:

coeff_lm1 <- as_tibble(summary(lm1)$coefficients, rownames = "erklaerende") %>%
  janitor::clean_names() %>%    # Spaltennamen bereinigen
  mutate(erklaerende = factor(erklaerende),   # in Faktor überführen für Diagramm
         erklaerende = forcats::fct_reorder(erklaerende, estimate),  #Sortieren für Diagramm
         vis_sig = cut(pr_t, breaks = c(0, 0.05, 1), 
                       labels = c("signifikant (p<=0,05)", "nicht signifikant (p>0,05)")))  # Signifikanz als Faktor für Diagramm

ggplot(coeff_lm1, aes(x = erklaerende, y = estimate, fill = vis_sig)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  labs(x = "Erklärende Variablen\n",
       y = "\nRegressionskoeffizienten",
       fill = "Signifikanz")

Bei den metrischen Variablen sehen wir, dass der Anteil der HochschulabsolventInnen und der Anteil der über 65-Jährigen signifikant positiv auf die Impfquote wirken. Die Mortalität und der Anteil der Nicht-ÖsterreicherInnen würden negativ auf die Impfquote einwirken, wenn sie denn signifikant wären. Bei unserer nominalen Bundeslandzugehörigkeit sehen wir, dass nur der Intercept - also die Konstante der Regressionsgleichung - klar positiv auf die Impfquote einwirkt. Und hier knüpfen wir nochmals an die vorherigen Ausführungen aus Kapitel 6.5 an: Dieser Intercept steht inhaltlich für die Linearkombination der acht verwendeten Dummy-Codes, also für das Bundesland Burgenland. Die Lage eines Bezirkes im Burgenland trägt also signifikant zu einer Erhöhung der Impfquote. Alle restlichen Bundesländer weisen negative Regressionskoeffizienten auf und verringern damit die Impfquote. Am deutlichsten in den Bundesländern Tirol, Salzburg, Oberösterreich und Kärnten.

6.6.4 Modelloptimierung

Um unser Modell im Sinne von Ockhams Razor noch etwas schlanker zu gestalten, wollen wir noch die nicht signifikanten metrischen Variablen entfernen:

lm2 <- lm(anteil_immun ~ tote_100k + bev_anteil_65plus 
          + bildung_anteil_hochschule + bld, data = sel_daten_trans)
summary(lm2)

## 
## Call:
## lm(formula = anteil_immun ~ tote_100k + bev_anteil_65plus + bildung_anteil_hochschule + 
##     bld, data = sel_daten_trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.50221 -0.27361 -0.01373  0.29486  1.46838 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                1.30888    0.22783   5.745 1.55e-07 ***
## tote_100k                 -0.16023    0.08526  -1.879 0.063806 .  
## bev_anteil_65plus          0.27093    0.09377   2.889 0.004954 ** 
## bildung_anteil_hochschule  0.24021    0.07418   3.238 0.001745 ** 
## bldKtn.                   -2.04719    0.29525  -6.934 9.05e-10 ***
## bldNÖ                     -0.91416    0.25196  -3.628 0.000498 ***
## bldOÖ                     -2.00430    0.30022  -6.676 2.82e-09 ***
## bldSbg.                   -1.80298    0.35795  -5.037 2.82e-06 ***
## bldStmk.                  -0.93645    0.30226  -3.098 0.002677 ** 
## bldT                      -1.72894    0.34038  -5.079 2.38e-06 ***
## bldVbg.                   -1.21734    0.42205  -2.884 0.005023 ** 
## bldW                      -1.54396    0.70706  -2.184 0.031882 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6171 on 81 degrees of freedom
## Multiple R-squared:  0.6647, Adjusted R-squared:  0.6192 
## F-statistic:  14.6 on 11 and 81 DF,  p-value: 4.44e-15

Das Entfernen des Anteils der Nicht-ÖsterreicherInnen bewirkt, dass im Vergleich zu lm1 die Bundeslandzugehörigkeit zu Wien signifikant zu einer verringerten Impfquote beiträgt. Nennesswerte Verschiebungen in den Regressionskoeffizienten der Dummy-Variablen zur Bundeslandzugehörigkeit können wir aber nicht beobachten.

Bliebe noch die Mortalität zu entfernen:

lm3 <- lm(anteil_immun ~ bev_anteil_65plus + bildung_anteil_hochschule + bld, 
          data = sel_daten_trans)
summary(lm3)

## 
## Call:
## lm(formula = anteil_immun ~ bev_anteil_65plus + bildung_anteil_hochschule + 
##     bld, data = sel_daten_trans)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.62136 -0.29075 -0.02028  0.35966  1.69079 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                1.37810    0.22827   6.037 4.37e-08 ***
## bev_anteil_65plus          0.25726    0.09492   2.710 0.008187 ** 
## bildung_anteil_hochschule  0.27473    0.07297   3.765 0.000312 ***
## bldKtn.                   -2.19986    0.28820  -7.633 3.71e-11 ***
## bldNÖ                     -0.93197    0.25564  -3.646 0.000467 ***
## bldOÖ                     -2.06545    0.30302  -6.816 1.45e-09 ***
## bldSbg.                   -1.90508    0.35922  -5.303 9.44e-07 ***
## bldStmk.                  -1.18207    0.27672  -4.272 5.19e-05 ***
## bldT                      -1.67698    0.34445  -4.869 5.38e-06 ***
## bldVbg.                   -1.16605    0.42762  -2.727 0.007818 ** 
## bldW                      -1.74226    0.70985  -2.454 0.016229 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6266 on 82 degrees of freedom
## Multiple R-squared:  0.6501, Adjusted R-squared:  0.6074 
## F-statistic: 15.24 on 10 and 82 DF,  p-value: 5.616e-15

🤔 Was hat uns das Entfernen dieser beiden metrischen Variablen nun gebracht?

Einerseits habt sich die Varianzaufklärung unseres Modells dadurch graduell von knapp 63% auf 60% verschlechtert. Gleichzeitig konnten wir aber auch die Anzahl der erklärenden Variablen von fünf auf drei reduzieren. Ein Tradeoff den wir hier einmal akzeptieren können: Mit fast der Hälfte an Variablen nahezu die gleiche Varianzaufklärung zu liefern, ist kein schlechter Deal 😉

In Summe hat sich gezeigt, dass die Bundeslandzugehörigkeit unser ursprünglich rein metrisches Modell deutlich verbessert hat. Die im Laufe unserer Analyse entfernten Variablen wie der Anteil der FPÖ-WählerInnen, die Corona-bedingte Mortalität oder der Anteil der Nicht-ÖsterreicherInnen in den Bezirken erwiesen sich zur Erklärung der bezirksweisen Impfquoten als nicht aussagekräftig.

6.7 Abschließendes Prüfen der Modellannahmen

Analog zu Kapitel 5.8 wollen wir zuletzt noch überprüfen, ob unser Modell lm3 folgende Modellannahmen erfüllt:

Linearer Zusammenhang gegeben
Residuen = normalverteilt
Erwartungswert der Residuen = 0
Residuen haben konstante Varianz (Homoskedastiziät)
Residuen sind unkorreliert zu Beobachtungen (Autokorrelation)

6.7.1 Prüfung des linearen Zusammenhangs

Anhand des “Residuals vs. Fitted”-Plot sehen wir, dass

kein klares Muster in der Punktverteilung ersichtlich ist;
und die Regressionslinie naher der Null-Linie zu liegen kommt.

plot(lm3, 1, labels.id = sel_daten_trans$bez_txt)

Wir können also von einem linearen Zusammenhang ausgehen.

6.7.2 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”. In unserem Fall sehen wir eine hinreichend gute Annäherung an eine Normalverteilung - also an die strichlierte Hauptdiagonale:

plot(lm3, 2, labels.id = sel_daten_trans$bez_txt)

## Warning: not plotting observations with leverage one:
##   93

📚 Exkurs: … Leverage one

Die Adleraugen haben sicherlich den Fehler im obigen Output erkannt: Warning: not plotting observations with leverage one: 93. Diese Fehlermeldung lässt sich in unserem Fall auf die Datenstruktur zurückführen: Kategoriale Variablen bei denen ein Level nur mit einer Beobachtung besetzt sind, verursachen Probleme bei der Berechnung von Residuen. In unserem Fall trifft das auf den Bezirk Wien (Rekord Nr. 93), da er der einzige Bezirk im Bundesland Wien ist.
Kurz und gut: Aufgrund unserer speziellen Datenstruktur können wir diese Fehlermeldung ignorieren.

6.7.3 Prüfung des Erwartungswerts der Residuen

Der Erwartungswert (= Summe) der Residuen sollte bei 0 liegen. Um dies zu überprüfen:

sum(residuals(lm3))

## [1] 4.614364e-16

Jep, bei 15 Nullen nach dem Komma, können wir von einer hinreichenden Approximation von 0 ausgehen 😉 .

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

Im Scale-Location-Plot sehen wir …

plot(lm3, 3, labels.id = sel_daten_trans$bez_txt)

… einigermaßen gleichmäßig über die Fitted-Values verteilte Residuen (= annähernd horizontal verlaufende Regressionsgerade). Wir können somit von einer Homogenität in den Varianzen (“Homoskedastizität”) ausgehen.

6.7.5 Prüfung auf Autokorrelation

Eine weitere Annahme linearer Modelle ist die

Um final zu überprüfen, ob die Unabhängigkeit der Residuen von der Reihenfolge der Beobachtungen. gegeben ist:

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

Und wieder ist ein “Jep!” angebracht: Wir sehen in der Verteilung der Residuen über die Beobachtungen keine auffälligen Muster, was auch durch die Lagen der Regressionsgeraden entlang der Nullline bestätigt wird. Wir können also eine Autokorrelation ausschließen.

🏆 Nun wissen wir, …

wie wir direkt in R Datenmanipulationen wie die Transformation einer numerischen Variable in einen Faktor vornehmen können;
wie die Dummy-Kodierung von Faktoren in R abläuft;
wie wir Visualisierungen zur Interpretation von Modellparametern nutzen können;
dass die Bundeslandzugehörigkeit eines Bezirkseinen starken Beitrag zur Erklärung der Impfquote liefert.

Einfach faszinierend!