invisible header
ESS8-Datensatz einlesen.
# Working directory setzen (z.B. "c:/daten")
setwd("mein_laufwerk/mein_datenverzeichnis")
# Daten einlesen
library(haven)
datensatz <- read_dta("daten.dta")
Ein neues Übungsskript erstellen.
# Statistik 2: R Tutorat
# Übungsskript zum Thema Regressionsanalyse
# Datum: 13.04.2021
# AutorIn: XXX
Installation und Aktivierung der Packages tidyverse und ggplot2.
# install.packages("tidyverse")
library(tidyverse)
# install.packages("ggplot2")
library(ggplot2)
1. Wir untersuchen das Merkmal “Klimawandelverklärung” (ccgdbd) und fragen uns, ob dies in der Schweiz vom Alter (agea) abhängt. Selektiere dazu zunächst alle SchweizerInnen aus dem ESS. Benenne den neuen Teildatensatz “ess8.CH”.
ess8_CH <- filter(ess8, cntry == "CH")
2. Schränke den Datensatz auf die relevanten Variablen ein und benenne den erstellten Teildatensatz “ess8_CH_ss”. Recherchiere, was genau die Variable ccgdbd misst & wie sie kodiert ist.
ess8_CH_ss <- select(ess8_CH, cntry, idno, agea, ccgdbd)
3. Inspiziere und bereinige die Daten: Wie viele Merkmalsträger (Fälle) gibt es? Was sind die Verteilungseigenschaften bzw. Kennwerte der Variablen? Müssen eventuell Missings rekodiert werden?
dim(ess8_CH_ss)
## [1] 1525 4
Der Datensatz enthält nun 1525 Merkmalsträger und 4 Variablen.
attributes(ess8_CH_ss$agea)
## $label
## [1] "Age of respondent, calculated"
##
## $format.stata
## [1] "%13.0g"
##
## $labels
## Not available
## NA
##
## $class
## [1] "haven_labelled" "vctrs_vctr" "double"
attributes(ess8_CH_ss$ccgdbd)
## $label
## [1] "Climate change good or bad impact across world"
##
## $format.stata
## [1] "%14.0g"
##
## $labels
## Extremely bad 1 2 3 4
## 0 1 2 3 4
## 5 6 7 8 9
## 5 6 7 8 9
## Extremely good Not applicable Refusal Don't know No answer
## 10 NA NA NA NA
##
## $class
## [1] "haven_labelled" "vctrs_vctr" "double"
Der attributes()-Output ist ziemlich unübersichtlich, enthält aber ein paar sinnvolle Information, z.B. die Variablenlabels und die verschiedenen Ausprägungen. Die Variable zur Klimaverklärung geht von 0 (‘Die Auswirkungen des Klimawandels auf der ganzen Welt sind äusserst negativ’) bis 10 (‘Die Auswirkungen des Klimawandels auf der ganzen Welt sind äusserst positiv’).
summary(ess8_CH_ss$agea)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 15.00 32.00 48.00 47.83 62.00 94.00 6
sd(ess8_CH_ss$agea, na.rm = TRUE)
## [1] 18.78213
Die Altersverteilung im Datensatz reicht von 15 bis 94 Jahren, liegt im Mittel bei etwas unter 48 und enthält 6 ausgewiesene fehlende Werte. Eine Altersmessung weicht im Schnitt um etwa 19 Jahre vom Durchschnittsalter ab.
summary(ess8_CH_ss$ccgdbd)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.000 2.000 3.000 3.227 5.000 10.000 55
sd(ess8_CH_ss$ccgdbd, na.rm = TRUE)
## [1] 2.101591
Im Schnitt ordnen sich die Personen bei einem Wert von etwa 3.2 auf der Klimaverklärungsskala ein - überwiegend werden die Auswirkungen des Klimawandels also negativ beurteilt. Die Standardabweichung liegt bei etwa 2.1, insgesamt 55 fehlende Angaben sind auch als solche ausgewiesen. Dies ist noch im Rahmen des üblichen und akzeptablen, muss nicht weiter beachtet oder berichtet werden.
Rekodierungen müssen nicht vorgenommen werden, alle vorhandenen numerischen Ausprägungen verweisen auf gültige Messungen. Entsprechend können wir sofort mit dem na.omit()-Befehl fortfahren, wobei sich die Fallzahl der gültigen Fälle von 1525 auf 1465 reduziert.
ess8_noNA <- na.omit(ess8_CH_ss)
dim (ess8_noNA)
## [1] 1465 4
4. Starte in die Regressionsanalyse mit der Visualisierung des Zusammenhangs. Erstelle zunächst ein einfaches Streudiagramm, dann diskutiere dies kritisch und biete eine modifizierte Darstellung an.
scatterplot <- ggplot(ess8_noNA,
aes(x = agea, y = ccgdbd)) +
geom_point() +
scale_x_continuous(breaks = seq(15,94,5)) +
scale_y_continuous(breaks = seq(0,10,1)) +
labs(x = "Alter in Jahren",
y = "Klimaverklärung",
title = "Alter und Klimaverklärung in der Schweiz",
subtitle ="Hat der Klimawandel negative oder positive Folgen?",
caption = "Quelle: ESS 2016 Teilstichprobe CH, n = 1482 \nKodierung: 0-sehr negative Folgen; 10-sehr positive Folgen.") +
theme_classic()
scatterplot
Oh weh, Overplotting! Viele Punkte liegen aufeinander, der Zusammenhang kann so nicht gedeutet werden. Hier ein Paar Lösungsvorschläge:
alphaplot <- ggplot(ess8_noNA,
aes(x = agea, y = ccgdbd)) +
geom_point(alpha = 0.1, size = 3) +
scale_x_continuous(breaks = seq(15,94,5)) +
scale_y_continuous(breaks = seq(0,10,1), labels=c("sehr negative Folgen","1", "2", "3", "4", "5", "6", "7", "8", "9", "sehr positive Folgen")) +
scale_fill_gradient(low = "lightblue",high = "darkblue") +
labs(x = "Alter in Jahren",
y = "Klimaverklärung",
title = "Alter und Klimaverklärung in der Schweiz",
subtitle ="Hat der Klimawandel negative oder positive Folgen?",
caption = "Quelle: ESS 2016 Teilstichprobe CH, n = 1482") +
theme_classic()
alphaplot
jitterplot <- ggplot(ess8_noNA,
aes(x = agea, y = ccgdbd)) +
geom_jitter(width=0, alpha=0.2, size = 3) +
scale_x_continuous(breaks = seq(15,94,5)) +
scale_y_continuous(breaks = seq(0,10,1), labels=c("sehr negative Folgen","1", "2", "3", "4", "5", "6", "7", "8", "9", "sehr positive Folgen")) +
scale_fill_gradient(low = "lightblue",high = "darkblue") +
labs(x = "Alter in Jahren",
y = "Klimaverklärung",
title = "Alter und Klimaverklärung in der Schweiz",
subtitle ="Hat der Klimawandel negative oder positive Folgen?",
caption = "Quelle: ESS 2016 Teilstichprobe CH, n = 1482") +
theme_classic()
jitterplot
heatplot <- ggplot(ess8_noNA,
aes(x = agea, y = ccgdbd)) +
geom_bin2d(bins = 10) +
scale_x_continuous(breaks = seq(15,94,5)) +
scale_y_continuous(breaks = seq(0,10,1), labels=c("sehr negative Folgen","1", "2", "3", "4", "5", "6", "7", "8", "9", "sehr positive Folgen")) +
scale_fill_gradient(low = "lightblue",high = "darkblue") +
labs(x = "Alter in Jahren",
y = "Klimaverklärung",
title = "Alter und Klimaverklärung in der Schweiz",
subtitle ="Hat der Klimawandel negative oder positive Folgen?",
caption = "Quelle: ESS 2016 Teilstichprobe CH, n = 1482") +
theme_classic()
heatplot
Zur Beschriftung der Skalenpunkte werden per Default die Wertelabels der Variablen verwendet. Alternativ haben wir hier die Unterfunktion labels=c(“Zeichenvektor”) der scale_y_continuous()-Funktion verwendet, um die Endpunkte der abhängigen Variablen zu verbalisieren. Dabei ist zu beachten, dass die Länge des Zeichenvektors gleich der in breaks= festgelegten Sequenz sein muss.
Es wird deutlich, dass die Schwerpunkte der Verteilung der Klimaverklärung durchgehend in der unteren Hälfte der Skala liegen. Gleichwohl lässt sich eine ganz leichte Verlagerung des Schwerpunktes der Verteilung nach oben (zu verklärter Beurteilung) mit zunehmendem Alter erkennen. Allerdings scheint der Zusammenhang, wenn überhaupt, nur sehr schwach ausgeprägt zu sein.
5. Berechne die Parameter der Regressionsgerade.
fit <- lm(ccgdbd ~ agea, data = ess8_noNA)
summary(fit)
##
## Call:
## lm(formula = ccgdbd ~ agea, data = ess8_noNA)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7214 -1.3349 -0.2325 1.4264 6.9266
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.686852 0.149362 17.989 < 2e-16 ***
## agea 0.011369 0.002913 3.903 9.93e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.093 on 1463 degrees of freedom
## Multiple R-squared: 0.01031, Adjusted R-squared: 0.009629
## F-statistic: 15.23 on 1 and 1463 DF, p-value: 9.926e-05
Mithilfe des summary()-Befehls lassen sich alle relevanten Parameter anzeigen.
6. Nutzt die grafische Aufbereitung aus der Aufgabe I.4 und setzt in diese eine Regressionsgerade ein, um euch bei den folgenden Interpretationsversuchen eine visuelle Stütze zu schaffen.
alphaplot + geom_smooth(method = "lm", se = FALSE, color = "red")
## `geom_smooth()` using formula 'y ~ x'
jitterplot + geom_smooth(method = "lm", se = FALSE, color = "red")
## `geom_smooth()` using formula 'y ~ x'
heatplot + geom_smooth(method = "lm", se = FALSE, color = "red")
## `geom_smooth()` using formula 'y ~ x'
7. Bietet sinnvolle Interpretationen des Regressionskoeffizienten bzw. der Regressionsgleichung an.
[Basis-Interpretation] Die Regressionsgerade hat einen Steigungskoeffizienten von b=0.011: mit jedem zusätzlichen Lebensjahr steigt der Grad der Klimaverklärung im Schnitt um etwa 0.01 Skalenpunkte.
[Inverse Vorhersage] Um also einen ganzen Skalenpunkt auf der 11er-Skala zu überschreiten, muss eine Person etwa 100 Jahre älter werden (1/0.01 = 100) - auf den ersten Blick mag dies als Ausdruck eines geringen Einflusses des Alters gedeutet werden.
Allerdings [Vorhersagebasierte Interpretation]:
library(ggeffects)
## Warning: package 'ggeffects' was built under R version 3.6.3
ggpredict(fit, terms = "agea[20, 65]")
## # Predicted values of Climate change good or bad impact across world
##
## agea | Predicted | 95% CI
## -------------------------------
## 20 | 2.91 | [2.72, 3.11]
## 65 | 3.43 | [3.28, 3.57]
Eine 20-jährige Person trennt laut Regressionsvorhersage fast ein halber Skalenpunkt von einer 65-jährigen Person [zusätzlich Standardisierung der Bezugseinheit]: Dies entspricht immerhin einer viertel Standardabweichung der Klimaverklärung (\frac{Differenz}{sd(Klimaverklärung)} = \frac{0.5}{2.1} = 0.23).
Conforti, E., Dürr, R., Siefart, F., Strassmann Rocha, D., Giesselmann, M. (2021): “Regressionsanalysen mit R”
Unter Mitarbeit von Norma De Min und Sebastian Senn