class: center, middle, inverse, title-slide .title[ # Wissenschaftliches Programmieren mit R ] .subtitle[ ## Elementare Statistik ] .author[ ### Andreas Blaette ] .date[ ### 1. Juni 2023 ] --- # Vorbereitungen Die Folien nutzen den Datensatz `corona_by_county`, der ab v0.0.3 im [learningR]()-Paket enthalten ist. Die Tabelle enthält zur Corona-Pandemie etliche Variablen auf Kreisebene (Landkreise, kreisfreie Städte). ```r remotes::install_github("ablaette/learningR") # Installation der aktuellen Version ``` ```r library(learningR) dim(corona_by_county) ``` ``` ## [1] 400 13 ``` Die HTML-Variante dieser Folien ist im Paket als "Vignette" enthalten. Sie können Sie wie folgt aufrufen: ```r browseVignettes(package = "learningR") ``` Um den Code selbst nachzuvollziehen, öffnen Sie die R Markdown-Datei, auf welcher der Foliensatz beruht (Auswahl von "f_statistics.Rmd"): ```r learningR::open_rmd_file() ``` --- # Über den Corona-Regionaldatensatz Der Corona-Regionaldatensatz `corona_by_county`, der im Paket enthalten ist, umfasst die folgenden Variablen: - *RS*: Regionalschlüssel - *name*: Name des Landkreises, der kreisfreien Stadt - *region*: Name des Bundeslandes - *type*: Landkreis / kreisfreie Stadt - *population_total*: Bevölkerung insgesamt - *cases7_per_100k*: 7-Tage-Inzidenz (Stand: Mai 2021) - *death_rate*: Todesrate - *foreign_pop_total*: Bevölkerung ohne deutsche Staatsangehörigkeit - *afd*: Prozentualer Anteil der AfD am Zweitstimmenergebnis - *income*: Durchschnittliches Haushaltseinkommen - *per_100km2*: Bevölkerungsdichte - *agequot*: Altersquotient - Anteil der Bevölkerung der mindestens 65 ist Der Datensatz wurde mit diesem R-Skript durch die Verknüpfung verschiedener Datenquellen automatisch erstellt. Die konsistente Nutzung des Regionalschlüssels (RS) in offiziellen Daten auf Kreisebene erleichtert diese Verknüpfung von Daten ungemein. --- # Statistische Daten: 'dplyr' als Standard - Bei der Arbeit mit statistischen Daten etabliert sich das Vokabular des [dplyr](https://dplyr.tidyverse.org)-Pakets als de facto-Standard | Funktion | Kurzbeschreibung | | ---------|------------------| | mutate() | Definition neuer Variablen, die von existierenden Variablen abgeleitet sind | | select() | Auswahl von Variablen | | filter() | Auswahl von Fällen entsprechend ihrer Werte | | summarise() | Erstellen von "summary statistics" | | arrange() | Umgruppierung von Spalten (Reihenfolge von Variablen) | - Das dplyr-Paket ist eingebettet in das "tidyverse" als eigenständiges Ökosystem von R-Paketen: - Ausgiebige Nutzung des "Pipe"-Operators des [magrittr](https://magrittr.tidyverse.org)-Pakets - Der klassische R-`data.frame` wird überwölbt vom `tibble` als Datenstruktur (u.a. bessere `show()`-Methode) - Das ["dplyr"-Cheatsheet](https://www.rstudio.com/wp-content/uploads/2015/02/data-wrangling-cheatsheet.pdf) bietet eine großartige Handreichung! --- # Deskriptive Statistik ```r income <- corona_by_county[["income"]] ``` ```r mean(income, na.rm = TRUE) ``` ``` ## [1] 22651.04 ``` ```r sd(income, na.rm = TRUE) ``` ``` ## [1] 2649.815 ``` ```r c(min = min(income, na.rm = TRUE), max = max(income, na.rm = TRUE)) ``` ``` ## min max ## 16450 36883 ``` ... doch welcher Landkreis hat das höchste durchschnittliche Haushaltseinkommen (den höchsten Ausländeranteil, die niedrigste Inzidenz etc.)? --- # Verteilungen ```r hist(corona_by_county[["income"]]) ``` <!-- --> --- # Deskriptive Statistik II (for-Schleife) ```r by_county_data <- subset(corona_by_county, !is.na(region)) ``` ```r regions <- unique(by_county_data[["region"]]) ``` ```r county_data_aggr_li <- list() for (buland in regions){ corona_by_county_min <- subset(corona_by_county, region == buland) county_data_aggr_li[[buland]] <- data.frame( mean = mean(corona_by_county_min[["agequot"]], na.rm = TRUE), sd = sd(corona_by_county_min[["agequot"]], na.rm = TRUE) ) } county_data_aggr <- do.call(rbind, county_data_aggr_li) head(county_data_aggr, 3) ``` ``` ## mean sd ## SH 23.27333 2.226935 ## HH 18.20000 NA ## NI 22.68889 2.562334 ``` --- # Deskriptive Statistik II (lapply-Variante) ```r regions <- unique(by_county_data[["region"]]) ``` ```r county_data_aggr_li <- lapply( regions, function(buland){ corona_by_county_min <- subset(corona_by_county, region == buland) data.frame( region = buland, mean = mean(corona_by_county_min[["agequot"]], na.rm = TRUE), sd = sd(corona_by_county_min[["agequot"]], na.rm = TRUE) ) } ) county_data_aggr <- do.call(rbind, county_data_aggr_li) head(county_data_aggr, 3) ``` ``` ## region mean sd ## 1 SH 23.27333 2.226935 ## 2 HH 18.20000 NA ## 3 NI 22.68889 2.562334 ``` --- # Deskriptive Statistik II (dplyr-Variante) ```r library(dplyr, quietly = TRUE, warn.conflicts = FALSE) corona_by_county %>% filter(!is.na(region)) %>% group_by(region) %>% summarise(mean = mean(agequot), sd = sd(agequot)) %>% head(6) ``` ``` ## # A tibble: 6 × 3 ## region mean sd ## <chr> <dbl> <dbl> ## 1 BB 25.6 2.71 ## 2 BW 20.5 1.76 ## 3 BY 21.2 1.99 ## 4 HB 21.5 0.707 ## 5 HE 21.5 2.57 ## 6 HH 18.2 NA ``` --- # Deskriptive Statistik II (data.table) ```r library(data.table, warn.conflicts = FALSE) dt <- data.table(corona_by_county) dt[!is.na(region), list(mean = mean(.SD$agequot), sd = sd(.SD$agequot)), by = "region"] ``` ``` ## region mean sd ## 1: SH 23.27333 2.2269347 ## 2: HH 18.20000 NA ## 3: NI 22.68889 2.5623340 ## 4: HB 21.50000 0.7071068 ## 5: NW 21.42075 1.5197872 ## 6: HE 21.48462 2.5671295 ## 7: RP 22.37500 2.0793371 ## 8: BW 20.51591 1.7590885 ## 9: BY 21.24687 1.9945163 ## 10: SL 24.31667 0.9368387 ## 11: BB 25.56667 2.7062020 ## 12: MV 25.22500 1.0819955 ## 13: SN 27.33846 3.0357697 ## 14: ST 27.22143 2.3664900 ## 15: TH 26.63043 2.6929934 ``` --- # Formeln in R: Boxplot ```r boxplot(agequot ~ region , data = by_county_data, las = 2) ``` <!-- --> --- # Formeln in R: lattice ```r lattice::histogram(~ afd | region, data = corona_by_county) ``` <!-- --> --- # Formeln in R: lattice ```r library(lattice) vars <- c("cases7_per_100k", "foreign_pop_share", "afd", "income", "per_km2", "agequot") lattice::splom(~ corona_by_county[, vars], data = corona_by_county) ``` <!-- --> --- # Lineare Regression Fragestellungen / Hypothesen: - In der Diskussion um die Triebkräft der Pandemie wurde im Mai 2021 die Beobachtung diskutiert, Personen mit Migrationsbiographie wären besonders stark von COVD-19 betroffen. - Verhaltensvorschriften, welche die Ausbreitung der Pandemie dämpfen sollen, finden bei der AfD wenig Akzeptanz. Zum Teil werden diese offen abgelehnt. Führt eine mangelnde Akzeptanz von Verhaltensregeln (Kontaktreduktion, Tragen von Masken, Abstandsregeln), welche sich regional in einem hohen Anteil der AfD am Zweitstimmenergebnis ausdrücken könnte, zu einer höheren Inzidenz? Bei Berechnung einer Regression mit der `lm()`-Funktion kommt die bereits bekannte Formelnotation zum Einsatz: - Die "Formel" übersetzt die Regressionsgleichung. - Die abhängige Variable (AV) steht links vom "Tilde"-Zeichen (~), die unabhängigen Variablen (UV) rechts von diesem. ```r r <- lm(cases7_per_100k ~ afd, data = by_county_data) ``` --- # Lineare Regression - summary() ```r summary(r) ``` ``` ## ## Call: ## lm(formula = cases7_per_100k ~ afd, data = by_county_data) ## ## Residuals: ## Min 1Q Median 3Q Max ## -62.035 -20.107 -4.908 15.067 171.005 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 38.7553 4.0659 9.532 < 2e-16 *** ## afd 1.8520 0.2821 6.565 1.63e-10 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 30.07 on 398 degrees of freedom ## Multiple R-squared: 0.09771, Adjusted R-squared: 0.09544 ## F-statistic: 43.1 on 1 and 398 DF, p-value: 1.626e-10 ``` --- # Das lm-Objekt - str() ```r str(summary(r)) ``` ``` ## List of 11 ## $ call : language lm(formula = cases7_per_100k ~ afd, data = by_county_data) ## $ terms :Classes 'terms', 'formula' language cases7_per_100k ~ afd ## .. ..- attr(*, "variables")= language list(cases7_per_100k, afd) ## .. ..- attr(*, "factors")= int [1:2, 1] 0 1 ## .. .. ..- attr(*, "dimnames")=List of 2 ## .. .. .. ..$ : chr [1:2] "cases7_per_100k" "afd" ## .. .. .. ..$ : chr "afd" ## .. ..- attr(*, "term.labels")= chr "afd" ## .. ..- attr(*, "order")= int 1 ## .. ..- attr(*, "intercept")= int 1 ## .. ..- attr(*, "response")= int 1 ## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> ## .. ..- attr(*, "predvars")= language list(cases7_per_100k, afd) ## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric" ## .. .. ..- attr(*, "names")= chr [1:2] "cases7_per_100k" "afd" ## $ residuals : Named num [1:400] -47.1 -7.77 -11.83 -14.19 -45.67 ... ## ..- attr(*, "names")= chr [1:400] "1" "2" "3" "4" ... ## $ coefficients : num [1:2, 1:4] 38.755 1.852 4.066 0.282 9.532 ... ## ..- attr(*, "dimnames")=List of 2 ## .. ..$ : chr [1:2] "(Intercept)" "afd" ## .. ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)" ## $ aliased : Named logi [1:2] FALSE FALSE ## ..- attr(*, "names")= chr [1:2] "(Intercept)" "afd" ## $ sigma : num 30.1 ## $ df : int [1:3] 2 398 2 ## $ r.squared : num 0.0977 ## $ adj.r.squared: num 0.0954 ## $ fstatistic : Named num [1:3] 43.1 1 398 ## ..- attr(*, "names")= chr [1:3] "value" "numdf" "dendf" ## $ cov.unscaled : num [1:2, 1:2] 0.018287 -0.001179 -0.001179 0.000088 ## ..- attr(*, "dimnames")=List of 2 ## .. ..$ : chr [1:2] "(Intercept)" "afd" ## .. ..$ : chr [1:2] "(Intercept)" "afd" ## - attr(*, "class")= chr "summary.lm" ``` --- # Lineare Regression im Scatterplot ```r plot(by_county_data$afd, by_county_data$cases7_per_100k) abline(r) ``` <!-- --> --- # Multivariate Regression: Vorbereitungen Vor der Berechnung einer Regression sollten die Variablen normalisiert werden. ```r normalize <- function(x){ ((x - max(x, na.rm = TRUE)) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE))) + 1 } by_county_data_norm <- mutate( .data = by_county_data, cases7_per_100k = normalize(cases7_per_100k), foreign_pop_share = normalize(foreign_pop_share), afd = normalize(afd), income = normalize(income), per_km2 = normalize(per_km2), agequot = normalize(agequot) ) ``` Alternativ: ```r vars <- c("cases7_per_100k", "foreign_pop_share", "afd", "income", "per_km2", "agequot") normalize_col <- function(x){ x <- by_county_data[[x]] (max(x, na.rm = TRUE) - x) / (max(x, na.rm = TRUE) - min(x, na.rm = TRUE)) } df_norm <- data.frame(lapply(setNames(vars, vars), normalize_col)) ``` --- # Modelle im Vergleich Berechnung alternativer Modelle ```r m1 <- lm(cases7_per_100k ~ foreign_pop_share, data = by_county_data_norm) m2 <- lm(cases7_per_100k ~ foreign_pop_share + income, data = by_county_data_norm) m3 <- lm(cases7_per_100k ~ foreign_pop_share + income + afd, data = by_county_data_norm) m4 <- lm( cases7_per_100k ~ foreign_pop_share + income + afd + agequot + per_km2, data = by_county_data_norm ) ``` Für den Vergleich der Modelle nutzen wir das stargazer-Paket: ```r stargazer::stargazer( m1, m2, m3, m4, dep.var.labels.include = FALSE, type = "html", align = TRUE, single.row = TRUE, column.sep.width = "1pt", omit.stat = c("n", "rsq", "f", "ser") ) ``` --- # Modelle im Vergleich <table style="text-align:center"><caption><strong>Results</strong></caption> <tr><td colspan="5" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td colspan="4"><em>Dependent variable:</em></td></tr> <tr><td></td><td colspan="4" style="border-bottom: 1px solid black"></td></tr> <tr><td style="text-align:left"></td><td>(1)</td><td>(2)</td><td>(3)</td><td>(4)</td></tr> <tr><td colspan="5" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">foreign_pop_share</td><td>0.179^*** (0.043)</td><td>0.200^*** (0.045)</td><td>0.331^*** (0.043)</td><td>0.424^*** (0.069)</td></tr> <tr><td style="text-align:left">income</td><td></td><td>-0.071 (0.054)</td><td>0.064 (0.051)</td><td>0.024 (0.055)</td></tr> <tr><td style="text-align:left">afd</td><td></td><td></td><td>0.377^*** (0.039)</td><td>0.374^*** (0.043)</td></tr> <tr><td style="text-align:left">agequot</td><td></td><td></td><td></td><td>0.007 (0.052)</td></tr> <tr><td style="text-align:left">per_km2</td><td></td><td></td><td></td><td>-0.123^* (0.063)</td></tr> <tr><td style="text-align:left">Constant</td><td>0.204^*** (0.013)</td><td>0.221^*** (0.018)</td><td>0.042^* (0.025)</td><td>0.042 (0.035)</td></tr> <tr><td colspan="5" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Adjusted R²</td><td>0.040</td><td>0.043</td><td>0.224</td><td>0.227</td></tr> <tr><td colspan="5" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td colspan="4" style="text-align:right">^*p<0.1; ^**p<0.05; ^***p<0.01</td></tr> </table> --- # Diskussion - Die Korrelation von Einwanderungssitutation und Infektionsgeschehen wird bestätigt. Hinweise, dass die soziale Lage (beengte Wohnsituation etc., hier näherungsweise gemessen durch Haushaltseinkommen) die eigentliche Quelle der Infektionsdynamik ist können nicht bestätigt werden. - Der Regionaldatensatz, der keine Aussagen auf individueller Ebene zulässt, stößt hier an die Grenzen seiner Aussagekraft! - Auch bei dem "AfD-Effekt" (eigentlich: Ablehnung von Verhaltensvorschriften => Infektionsdynamik) ist noch Vorsicht geboten: Kreise mit einem starken Abschneiden der AfD sind oft grenznah. Durch eine weitere Kontrollvariable wäre zu prüfen, dass die Nähe zu Hochinzidenz-Gebieten im Ausland die eigentlich Ursache der Infektionen ist. - Im Datensatz sind weitere Annahmen gemacht, z.B. Stabilität des Wahlverhaltens 2017-2021 ... - Und natürlich haben wir die statistischen Möglichkeiten nur angerissen, es gibt noch viel zu entdecken (Sauer, 2019) --- # Logistische Regression: Daten Wenn die abhängige Variable binär codiert ist, kommt die logistische Regression zum Einsatz. Ein typisches Szenario hierfür ist in der Politikwissenschaft die Wahl einer Partei. Im folgenden Szenario interessieren wir uns für die möglichen Gründe einer Wahlentscheidung zugunsten der AfD bei der Bundestagswahl 2017. Als Datengrundlage verwenden wir die GLES-Nachwahlbefragung 2017. Für die Vorbereitung des Datensatzes nutzen wir Funktionen des dplyr-Pakets. ```r library(dplyr) library(gles) ``` ``` ## ## Attaching package: 'gles' ``` ``` ## The following object is masked from 'package:learningR': ## ## btw2017 ``` ```r bt_min <- filter(bt2017nw, !is.na(q19ba)) %>% mutate(AfD = as.character(haven::as_factor(q19ba)) == "AfD") %>% mutate(income = as.integer(haven::as_factor(q192))) %>% mutate(dissatisfaction = as.integer(haven::as_factor(q33))) %>% mutate(citysize = as.integer(haven::as_factor(q197))) ``` Beachte: Die `glm()`-Funktion könnte auch Faktoren "verdauen"! --- # Logistische Regression: Modelle Die abhängige Variable unseres ist die binär codierte Variable "AfD" mit den logischen Werten `TRUE` und `FALSE`. Wenn eine Befragte keine Angabe zur Parteiwahl gemacht hat (`NA`-Werte), haben wir den Fall aus der Analyse ausgeschlossen. Wir verfolgen - sehr grobschlächtig! - drei Hypothesen, welche die Wahl der AfD erklären könnten: - Die ökonomisch "Abgehängten" wählen AfD: Je höher das Einkommen, desto weniger wahrscheinlich ist die Wahl der Afd. - Die mit der Demokratie Unzufriedenen wählen die AfD: Je höher die Demokratiezufriedenheit, desto unwahrscheinlicher ist die Wahl der AfD. - Die AfD hat in den kosmopolitisch geprägten Städten weniger Chancen: Je größer die Stadt, desto weniger wahrscheinlich ist die Wahl der AfD. Wir berechnen hierzu zwei Modelle. Die Notation verläuft parallel zur linearen Regression ... ```r m1 <- glm(AfD ~ income, data = bt_min, family = binomial()) m2 <- glm(AfD ~ income + dissatisfaction + citysize, data = bt_min, family = binomial()) ``` --- # Literatur Field, A., J. Miles, and Z. Field (2012). _Discovering Statistics Using R_. Los Angeles et al.: SAGE. Sauer, S. (2019). _Moderne Datenanalyse mit R: Daten einlesen, aufbereiten, visualisieren, modellieren und kommunizieren_. Wiesbaden: Springer.