In martscht/PsyBSc7: Übungen für die Praktika im Sommersemester 2020
source('setup.R')
library(learnr) # package to generate interactive HTML
library(gradethis) # grade Code from input in HTML
library(shiny) # nice appearance in R
library(fontawesome) # nice fonts
library(car)
library(MASS)
library(ggplot2)
library(lm.beta) # erforderlich für standardiserte Gewichte
data("PISA2009", package = "PsyBSc7")
data("WorldPopulation", package = "PsyBSc7")
knitr::opts_chunk$set(exercise.checker = gradethis::grade_learnr)
Einleitung und Datensatz (PISA: wie in den Übungen zum letzten Termin)
$*$ zum Teil basierend auf Unterlagen von Prof. Johannes Hartig
In dieser Sitzung werden wir uns weitere nichtlineare Effekte in Regressionsmodellen ansehen. Dazu verwenden wir zunächst den Datensatz aus der Übung zur PsyBSc7::Sitzung_7, welcher in den Aufgaben verwendet wurde und dessen Analysen in PsyBSc7::Loesung_7 mit dem richtigen Passwort einsehbar waren.
Der Beispieldatensatz enthält Daten zu Lesekompetenz aus der deutschen Stichprobe der PISA-Erhebung in Deutschland 2009.
library(car)
library(MASS)
library(lm.beta) # erforderlich für standardiserte Gewichte
library(ggplot2)
# Datensatz laden
data("PISA2009", package = "PsyBSc7")
Quadratische Verläufe in der Vorhersage von Lesekompetenz mit individuellen Merkmalen der Schüler/innen
$**$ zum Teil basierend auf Unterlagen von Prof. Johannes Hartig
In der Übung zur letzten Sitzung hatten wir herausgefunden, dass der Sozialstatus (HISEI), der Bildungsabschluss der Mutter (MotherEdu) und die Zahl der Bücher zu Hause (Books) bedeutsame Prädiktoren für die Lesekompetenz der Schüler/innen sind, allerdings zeigten Analysen, dass nicht alle Voraussetzungen erfüllt waren:
# Berechnung des Modells und Ausgabe der Ergebnisse
m1 <- lm(Reading ~ HISEI + MotherEdu + Books, data = PISA2009)
summary(lm.beta(m1))
Die Residuenplots sowie die Testung auf quadratische Trends zeigen an, dass für den Bildungsabschluss der Mutter auch eine quadratische Beziehung mit der Lesekompetenz besteht:
# Residuenplots
residualPlots(m1, pch = 16)
Die Effekte von Sozialstatus und Büchern werden durch das lineare Modell gut wiedergegeben. Für den Bildungsabschluss der Mutter ist ein leicht nicht-linearer Zusammenhang zu erkennen, der quadratische Trend für die Residuen ist signifikant (signifikantes Ergebnis für den Bildungsabschluss der Mutter). Der Effekt ist dadurch charakterisiert, dass der Zuwachs der Lesekompetenz im unteren Bereich des mütterlichen Bildungsabschlusses stärker ist und im oberen Bereich abflacht.
Auch dem Histogramm war eine Schiefe zu entnehmen, welche durch nichtlineare Terme entstehen können (im niederen Bereich liegen mehr Werte; eine Linksschiefe/Rechtssteile ist zu erkennen).
res <- studres(m1) # Studentisierte Residuen als Objekt speichern
df_res <- data.frame(res) # als Data.Frame für ggplot
# Grafisch: Histogramm mit Normalverteilungskurve
library(ggplot2)
ggplot(data = df_res, aes(x = res)) +
geom_histogram(aes(y =..density..),
bins = 15, # Wie viele Balken sollen gezeichnet werden?
colour = "blue", # Welche Farbe sollen die Linien der Balken haben?
fill = "skyblue") + # Wie sollen die Balken gefüllt sein?
stat_function(fun = dnorm, args = list(mean = mean(res), sd = sd(res)), col = "darkblue") + # Füge die Normalverteilungsdiche "dnorm" hinzu und nutze den empirischen Mittelwert und die empirische Standardabweichung "args = list(mean = mean(res), sd = sd(res))", wähle dunkelblau als Linienfarbe
labs(title = "Histogramm der Residuen mit Normalverteilungsdichte", x = "Residuen") # Füge eigenen Titel und Achsenbeschriftung hinzu
# Test auf Abweichung von der Normalverteilung mit dem Shpiro Test
shapiro.test(res)
Die Frage ist nun, woher die Verstöße gegen die Normalverteilungsannahme kommen. Erste Indizen aus den Partialplots wiesen darauf hin, dass möglicherweise ein quadratischer Effekt des Bildungsabschlusses der Mutter besteht.
Aufnahme eines quadratischen Effekts
$***$ zum Teil basierend auf Unterlagen von Prof. Johannes Hartig
Wird für den Bildungsabschluss der Mutter mit der Funktion poly ein linearer und qudratischer Trend in das Regressionsmodell aufgenommen, wird der quadratische Trend signifikant und das Modell erklärt signifikant mehr Varianz als ohne den quadratischen Trend: Um diese Ergebnisse zu sehen müssen wir zunächst ein quadratisches Regressionsmodell schätzen. Wir interessieren uns anschließend für die standardisierten Ergebnisse (summary und lm.beta). Den quadratischen Verlauf erhalten wir, indem wir innerhalb des linearen Modells poly auf den Bildungsabschluss der Mutter anwenden. poly nimmt als zweites Argument die Potenz, für welche wir uns interessieren; hier 2:
# Modell mit linearem + quadratischem Trend für MotherEdu:
m1.b <-
...
# Verwenden Sie den Befehl
lm(...~...)
# Verwenden Sie den Befehl
lm(...~...) # sowie innerhalb von lm
poly(..., 2)
# Verwenden Sie den Befehl
lm(...~...) # sowie innerhalb von lm
poly(..., 2) # wenden Sie dies auf MotherEdu an
# Verwenden Sie den Befehl
lm(...~...) # sowie innerhalb von lm
poly(..., 2) # wenden Sie dies auf MotherEdu an
# wenden Sie anschließend
summary(lm.beta(...)) # an
# Finale Lösung
m1.b <- lm(Reading ~ HISEI + poly(MotherEdu, 2) + Books, data = PISA2009)
summary(lm.beta(m1.b))
# For further use in exercises
m1 <- lm(Reading ~ HISEI + MotherEdu + Books, data = PISA2009)
m1.b <- lm(Reading ~ HISEI + poly(MotherEdu, 2) + Books, data = PISA2009)
m1.b <- lm(Reading ~ HISEI + poly(MotherEdu, 2) + Books, data = PISA2009)
Mit dem folgenden Befehl können wir auf eine simple Weise das Inkrement bestimmen.
# Vergleich mit Modell ohne quadratischen Trend
summary(m1.b)$r.squared - summary(m1)$r.squared # Inkrement
Wir möchten dieses Inkrement auf Signifikanz prüfen. Dies geht mit dem anova Befehl.
# Intro into excercise
m1 <- lm(Reading ~ HISEI + MotherEdu + Books, data = PISA2009)
m1.b <- lm(Reading ~ HISEI + poly(MotherEdu, 2) + Books, data = PISA2009)
# Begin grading:
res <- anova(m1, m1.b)
sol2 <- anova(m1.b, m1)
grade_result(
fail_if(~ (!identical(.result, res) & !identical(.result, sol2)), 'Der Befehl anova muss in diesem Zusammenhang zwei Argumente entgegennehmen: die zwei Regressions-Objekte, die wir zuvor abgespeichert hatten. Schauen Sie in der Übung zuvor an, wie dieses Modell hieß, falls Sie sich nicht sicher sind.'),
fail_if(~ (identical(.result, sol2)), 'Das ist an sich RICHTIG, allerdings sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell (mit weniger Prädiktoren und Parametern, die zu schätzen sind) zuerst übergeben werden. Hier: m1, da sonst die df negativ sind und auch die Änderung in den Sum of Sq (Quadratsumme) negativ sind! `R` erkennt dies zwar und testet trotzdem die richtige Differenz auf Signifikanz, aber wir wollen uns besser vollständig korrekt aneignen!'),
pass_if(~ identical(.result, res)),
correct = 'Sehr gut! Dies ist der Modellvergleich des lineares und des quadratischen Modells. Hier sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell (mit weniger Prädiktoren und Parametern, die zu schätzen sind) zuerst übergeben werden. Hier: m1, da sonst die df negativ sind und auch die Änderung in den Sum of Sq (Quadratsumme) negativ sind! `R` erkennt dies zwar und testet trotzdem die richtige Differenz auf Signifikanz, aber wir wollen uns besser vollständig korrekt aneignen!',
incorrect = 'Leider falsch.',
glue_correct = '{.correct}',
glue_incorrect = '{.incorrect} {.message}')
Erzeugt man für das erweiterte Modell Residuenplots, ist der quadratische Trend beim Bildungsabschluss komplett verschwunden - er ist ja schon im Modell enthalten und bildet sich somit nicht mehr in den Residuen ab:
residualPlots(m1.b, pch = 16)
Was bedeutet nun dieser Effekt inhaltlich? Um dies genauer zu verstehen, stellen wir die um die anderen Variablen bereinigte Beziehung zwischen dem Bildungsabschluss der Mutter und der Leseleistung grafisch dar.
X <- scale(poly(PISA2009$MotherEdu, 2))
std_par_ME <- c(0.1588, -0.1436)
pred_effect_ME <- X %*% std_par_ME
std_ME <- X[,1]
data_ME <- data.frame(std_ME, pred_effect_ME)
ggplot(data = data_ME, aes(x = std_ME, y = pred_effect_ME)) + geom_point(pch = 16, col = "blue", cex = 4)+
labs(y = "std. Leseleistung | Others", x = "std. Bildungsabschluss der Mutter | Others",
title = "Standardisierte bedingte Beziehung zwischen\n Bildungsabschluss der Mutter und Leseleistung")
Für den Grafik-Code sowie weitere Informationen zu quadaratischen Effekten und Funktionen siehe Appendix A. Die Grafik zeigt die vorhergesagte Beziehung zwischen den standardisierten Werten des Bildungsabschlusses der Mutter sowie der Leseleistung. Hierbei steht erneut | für "gegeben" (wie etwa beim Partialplot mit avPlots aus der vergangenen Sitzung: PsyBSc7::Sitzung_7()). Wir sehen also den um die anderen Variablen im Modell bereinigten Effekt zwischen Bildungsabschluss und Leseleistung. Hierbei ist ein starker mittlerer Anstieg der Leseleistung (-1 bis ca. 0.1) für einen Anstieg des Bildungsabschlusses von deutlich unterdurchschnittlich bis durchschnittlich (von -2.5 bis 0) zu sehen. Danach ist die Beziehung zwischen Leseleistung und Bildungsabschluss fast horizontal (Veränderung geringer als 0.1), was dafür spricht, dass es für einen durchschnittlichen bis überdurchschnittlichen Bildungsabschluss der Mutter (von 0 bis 1.5) kaum eine Beziehung zwischen den Variablen gibt. Dies bedeutet, dass besonders im unterdurchschnittlichen Bereich der mütterlichen Bildung Unterschiede zwischen Müttern einen starken Zusammenhang mit der Leseleistung ihrer Kinder zeigen. Wenn das Bildungsniveau der Mutter jedoch durchschnittlich oder überdurchschnittlich ist, scheint der Zusammenhang beinahe zu verschwinden.
Die grobe Gestalt der Beziehung hätten wir auch aus dem Koeffizienten ablesen können. Der Koeffizient des quadratischen Teils war negativ, was für eine invers-u-förmige (konkave) Funktion steht. Das Einzeichnen hilft uns jedoch, das genaue Ausmaß zu verstehen (siehe auch Appendix A). Auch hatten wir gesehen, dass der lineare Teil des Bildungsabschlusses der Mutter keinen statistisch signifikanten Beitrag zur Vorhersage geleistet hat. Jedoch gehört zu einer quadratischen Funktion immer auch ihr linearer Anteil dazu. Aus diesem Grund können wir unsere Stichprobe nur angemessen beschreiben, wenn wir den linearen Trend des Bildungsabschlusses der Mutter im Regressionsmodell beibehalten. Um das genaue Ausmaß besser zu verstehen, manipulieren Sie doch einmal die Beziehung, die wir soeben grafisch gesehen haben im folgenden interaktiven Block. Hierbei können Sie den linearen und den quadratischen Effekt verändern und sich die Auswirkungen auf die Grafik (die Beziehung zwischen Bildungsabschluss der Mutter und Leseleistung) ansehen. Die Default-Einstellungen sind identisch zu der oberen Grafik. curve plottet eine Linie und nimmt x automatisch als Argument der Funktion, somit wird $f(x)=\text{linear}x+\text{quadratisch}x^2$ geplottet. Probieren Sie doch einmal aus, was passiert, wenn Sie den linearen Teil auf 0 setzen oder das Vorzeichen des quadratischen Anteils ändern!
linear <- .1588
quadratisch <- -.1436
curve(linear * x + quadratisch * x^2,
xlim = c(-2, 2))
Mit Hilfe von poly(X, p), lassen sich Polynome bis zum Grad $p$ (als $X, X^2,\dots,X^{p-1},X^p$) in die Regression mit aufnehmen, ohne, dass sich die Parameterschätzungen der anderen Potenzen von $X$ ändern. Wenn Sie noch mehr über die Funktion poly und ihre Vorteile erfahren möchten, dann schauen Sie sich doch mal den Appendix A an.
Mit Hilfe von I() lassen sich innerhalb des lm Befehls zu dem noch weitere Funktionale hinzufügen, ohne diese vorher erzeugen zu müssen. Beispielsweise ließe sich durch lm(Y ~ X + I(sin(X))) + I(exp(sqrt(X)) folgendes Regressionsmodell schätzen: $Y = \beta_0+\beta_1X + \beta_2\sin(X) + \beta_3e^{\sqrt{X}} + \varepsilon$. Allerdings lassen sich so nicht Wachstumsraten modellieren (z.B. exponentielles oder logarithmisches Wachstum) - hierzu müssten die Variablen tatsächlich transformiert werden.
Weitere nichtlineare Effekte
Im Folgenden werden zum Teil fortgeschrittene Inhalte verwendet. Bitte sehen Sie dies als Bonus zum Vorlesungsstoff an.
Ein sehr wichtiges Vorhersagemodell zu Zeiten der Corona-Pandemie oder in Anbetracht von starkem weltweitem Bevölkerungswachstum ist das exponentielle Wachstum. Um dieses genauer zu verstehen, werden einige Grundlagen im Umgang mit Exponenten benötigt. Zum Beispiel müssen wir uns überlegen, welche Rate und welche Basis ein exponentieller Verlauf hat. Beispielsweise wächst $10^x$ deutlich schneller als $2^x$ (für gleiches wachsendes $x$). Wie gut kennen Sie sich noch mit Exponenten und Logarithmen aus?
question("Welche Aussagen über Logarithmen bzw. Exponenten stimmen?",
answer("Die Gleichung $10^x = 100$ können wir bspw. mit dem 10er-Logarithmus $log_{10}(100)$ berechnen. Die Antwort ist 2.", correct = T),
answer("Die Gleichung $10^x = 100$ können wir bspw. mit dem 10er-Logarithmus $log_{10}(100)$ berechnen. Die Antwort ist 3."),
answer("Logarithmen können NICHT in einander umgerechnet werden."),
answer("Es gilt: $log(a+b)=log(a)log(b)$."),
answer("Es gilt: $log(ab)=log(a)+log(b)$.", correct=T),
answer("$log_a(x)$ ist die Umkehrfunktion zu $a^x$. Dies bedeutet, dass $log_a(a^x)=x.$", correct = T),
answer("Es gilt: $a^{c+d}=a^c+a^d$."),
answer("Es gilt: $a^{c+d}=a^ca^d$.", correct = T),
answer("Es gilt immer $a^x < b^x$ für a < b und alle x."),
answer("Es gilt immer $a^x < b^x$ für 1 < a < b und alle x.", correct = T),
answer("Es gilt immer $a^x > b^x$ für 0 < a < b < 1 und alle x.", correct = T),
answer("Es gilt immer $a^x < b^x$ für 0 < a < b < 1 und alle x.")
)
Im folgenden Exkurs haben Sie die Möglichkeit Ihr Wissen nochmal aufzufrischen und zu erweitern. Ihnen werden Quellen genannt, in welchen Sie weitere Inhalte nachlesen können.
Exkurs: Exponenten und Logarithmen
Um eine exponentielle Gleichung nach $x$ aufzulösen wird der Logarithmus verwendet, denn er ist die Umkehrfunktion zum Exponentiellen (siehe beispielsweise hier für eine Wiederholung zum Logarithmus und Logarithmusrechenregeln; Pierce (2018), Exponenten; Pierce (2019) oder Exponentenrechenregeln; Pierce (2018)). In der Schule wird zunächst der Logarithmus zur Basis 10 gelehrt. So kann man die Gleichung $10^x = 100$ leicht auflösen: Dazu müssen beide Seiten (zur Basis 10) logarithmiert werden: $\log_{10}(10^x) = \log_{10}(100)$, denn es gilt im Allgemeinen:
$$\log_a(a)=1$$
$$\log_a(bc)=\log_a(b)+\log_a(c)$$
$$\log_a(b^c)=c\log_a(b)$$
$$x=\log_a(a^x)=a^{\log_a(x)}$$
(hier ist $\log_a$ der Logarithmus zur Basis $a$ und $a, b$ und $c$ sind Zahlen, wobei der Einfachheit halber $a,b,c>0$ gelten soll)
Also rechnen wir: $\log_{10}(10^x) = x\log_{10}(10)=x*1=x$ auf der linken Seite und $\log_{10}(100)=2$ auf der rechten Seite. So erhalten wir $x=2$ als Lösung. Natürlich wussten wir schon vorher, dass $10^2=100$ ergibt; es ist aber eine schöne Veranschaulichung, wofür Logarithmen unter anderem gebraucht werden.
Es muss allerdings nicht immer der richtige Logarithmus zur jeweiligen Basis sein. Eine wichtige Logarithmusrechenregel besagt, dass der Logarithmus zur Basis $a$ leicht mit Hilfe des Logarithmus zur Basis $b$ berechnet werden kann:
$$\log_{a}(x) = \frac{\log_{b}(x)}{\log_{b}(a)},$$ siehe bspw. Logarithmus und Logarithmusrechenregeln; Pierce (2018). Da diese Gesetzmäßigkeit gilt, wird in der Mathematik und auch in der Statistik meist nur der natürliche Logarithmus (ln; logarithmus naturalis) verwendet. Die Basis ist hierbei $e$ ($\approx 2.718\dots$). Mit den obigen Rechenregeln sowie der Rechenregeln für Exponenten können wir jedes exponentielles Wachstum als Wachstum zur Basis e darstellen, denn
$$a^x = e^{\text{ln}(a^x)}=e^{x\text{ln}(a)},\quad\text{ z.B. ist }\quad 10^3=e^{\text{ln}(10^3)}=e^{3\text{ln}(10)}=1000$$
Da sich durch e und ln alle exponentiellen Verläufe darstellen lassen, wird in der Mathematik häufig log als das Symbol für den natürlich Logarithmus verwendet; so ist es auch in R: ohne weitere Einstellungen ist log der natürliche Logarithmus und exp ist die Exponentialfunktion. Mit Hilfe von log(..., base = 10) erhalten sie beispielsweise den 10er-Logarithmus. Probieren Sie die obige Gleichung selbst aus:
# gleiches Ergebnis:
10^3
exp(3*log(10))
# gleiches Ergebnis:
log(10^3, base = 10) # Logarithmus von 1000 zur Basis 10
log(10^3)/log(10) # mit ln
# gleiches Ergebnis:
log(9, base = 3) # Logarithmus von 9 zur Basis 3
log(9)/log(3) # mit ln
Das Modellieren von exponentiellem Wachstum
Wir betrachten nun eine allgemeine exponentielle Wachstumsfunktion $$f(x) = a\ b^{c\ x}$$ (hierbei ist $a$ ein Vorfaktor, der die Ausprägung an der Stelle $x=0$ beschreibt, $b$ ist die Basis des exponentiellen Wachstums und $c$ ist ein eigentlich redundanter Ratenparameter). Wegen der Beliebigkeit der Basis ist dies gleich $$f(x) = e^{\text{ln}(a\ b^{c\ x})}=e^{\text{ln}(a) + \text{ln}(b)cx}.$$ Nun sind $\text{ln}(a)$ und $\text{ln}(b)c$ zwei Konstanten, die wir einfach umbenennen dürfen. Wir können sie beispielsweise $\beta_0$ und $\beta_1$ nennen; also $\beta_0 := \text{ln}(a)$ und $\beta_1:=\text{ln}(b)c$. Folglich steht $$f(x) = e^{\beta_0 + \beta_1x}$$ für beliebiges exponentielle Wachstum (wir erhalten den Verlauf "$2^x$" indem wir $\beta_0 = \text{ln}(a) = 0$ und $\beta_1=\text{ln}(b)c = \text{ln}(2)$ wählen; also $a=1, b=2$ und $c=1$). Gleichzeit bedeutet dies, dass $$\text{ln}(f(x)) = \beta_0 + \beta_1x$$ eine lineare Funktion ist, welche wir ganz einfach mit einer Regressionanalyse, also lm in R, untersuchen können! Wir können somit sagen, dass wir durch Logarithmieren von $f(x)$ in der Lage sind, das exponentielle Wachstum zu Linearisieren, also das exponentielle Wachstum in eine lineare Funktion zu transformieren, welche mit Auswertungsinstrumenten für lineare Funktionen untersucht werden können. Dies bedeutet, dass wir nach Logarithmieren der abhängigen Variable in der Lage sind die gesamte Klasse der exponentiellen Funktionen/des exponentiellen Wachstums für die Modellierung unserer Daten zu verwenden. Die folgenden zwei Grafiken verdeutlichen dies (Sie können hier die Parameter beliebig verändern und sich den Effekt auf die Grafiken ansehen; die horizontale gestrichelte Linien repräsentiert jeweils den Wert $a$ bzw. $\text{ln}(a)$, an welchem die Kurve $f$ bzw. $\ln(f)$ die y-Achse schneiden):
##################
#### Einstellen der Koeffizienten und berechnen von f(x)
####
x <- seq(-1,2,0.1) # x = Variablen (als Zahlen zwischen -1 und 2)
a <- 2 # Vorfaktor, der die Ausprägung an der Stelle x=0 beschreibt
b <- 3 # Basis des exponentiellen Wachstums
c <- 1.5 # *eigentlich redundanter* Ratenparameter
f <- a*b^(c*x) # f(x), eine exponentiell-wachsende Funktion in x
##################
#### Plot von f(x) vs. x
####
plot(x = x, y = f, type = "l", col = "blue", lwd = 2, main = "Plot von f(x) vs. x") # plotte f(x) gegen x
abline(v = 0, lwd = 0.7) # y-Achse, v = 0 zeichnet eine vertikale Linie bei x = 0
abline(h = a, lty = 3) # im Punkt a schneidet f (das exponentielle Wachstum) die y-Achse (x=0), h = a zeichnet zu y = a eine horizontale Linie
##################
#### Plot von ln(f(x)) vs. x
####
plot(x = x, y = log(f), type = "l", col = "blue", lwd = 2, main = "Plot von ln(f(x)) vs. x") # plotte ln(f(x)) gegen x
abline(v = 0, lwd = 0.7) # y-Achse, v = 0 zeichnet eine vertikale Linie bei x = 0
abline(h = log(a), lty = 3) # im Punkt log(a) schneidet log(f) (das linearisierte exponentielle Wachstum) die y-Achse (x=0), h =llog(a) zeichnet zu y = log(a) eine horizontale Linie
Modellierung von exponentiellem Wachstum
Das Modellieren von exponentiellem Wachstum am Beispiel der Weltbevölkerung von 1800 bis 2020
Nun können wir beispielsweise die Entwicklung der Weltbevölkerung von 1800 bis 2020 modellieren. Dazu müssen wir zunächst die Daten laden:
Die Dokumentation des Datensatzes mit Datenquellen sind hier einzusehen: gapminder.org-Dokumentationen. Im Paket PsyBSc7 ist unter dem Namen WorldPopulation ein vorverarbeiter Datensatz enthalten. Diesen können wir ganz einfach wie folgt laden und ansehen:
data("WorldPopulation", package = "PsyBSc7") # Bereits zusammengefasster Datensatz von 1800-2020
head(WorldPopulation) # Datenstruktur
In der ersten Spalte steht das Jahr; in der 2. die Weltbevölkerungsgröße. Wir wollen uns dies grafisch ansehen. Verwenden Sie ggplot, um die Population (Population) gegen das Jahr (Year) abzutragen und verwenden Sie geom_point().
ggplot(...)
# Sie müssen den Datensatz sowie die Namen der Variablen verwenden
ggplot(data = ..., aes(x = ..., y = ...))
# Sie müssen den Datensatz sowie die Namen der Variablen verwenden
ggplot(data = ..., aes(x = ..., y = ...)) + geom_point()
ggplot(data = WorldPopulation, aes(x = Year, y = Population))+geom_point()
Das Diagramm lässt deutlich einen nichtlinearen Anstieg der Weltbevölkerung von 1800 bis 2020 vermuten. Auffällig ist auch der leichte Knick, der um 1950 zu vermuten ist und ab welchem die Bevölkerung, deskriptiv gesehen, noch stärker wächst. Dieser Knick ist zum Teil durch das Ende des Krieges, aber auch durch modernere Landwirtschaft und das Aufkommen von neuen Medikamenten (z.B. Penicilline) zu erklären.
Lineares Modell für das Bevölkerungswachstum
Wir wollen uns naiverweise ein lineares Regressionmodell, also einen linearen Verlauf, der Weltbevölkerung vorhergesagt durch das Jahr, ansehen.
ggplot(data = WorldPopulation, aes(x = Year, y = Population))+
geom_point()+geom_smooth(method="lm", formula = "y~x") # plotte linearen Verlauf
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf
summary(m_l)
#########
### Normalverteilung der Residuen?
##
res <- studres(m_l) # Studentisierte Residuen als Objekt speichern
df_res <- data.frame(res) # als Data.Frame für ggplot
# Grafisch: Histogramm mit Normalverteilungskurve
ggplot(data = df_res, aes(x = res)) +
geom_histogram(aes(y =..density..),
bins = 10, # Wie viele Balken sollen gezeichnet werden?
colour = "blue", # Welche Farbe sollen die Linien der Balken haben?
fill = "skyblue") + # Wie sollen die Balken gefüllt sein?
stat_function(fun = dnorm, args = list(mean = mean(res), sd = sd(res)), col = "darkblue") + # Füge die Normalverteilungsdiche "dnorm" hinzu und nutze den empirischen Mittelwert und die empirische Standardabweichung "args = list(mean = mean(res), sd = sd(res))", wähle dunkelblau als Linienfarbe
labs(title = "Histogramm der Residuen mit Normalverteilungsdichte\n für das lineare Modell", x = "Residuen") # Füge eigenen Titel und Achsenbeschriftung hinzu
Durch + geom_smooth(method="lm", formula = "y~x"), kann mit ggplot ein linearer Trend inklusive Konfidenzintervall hinzugefügt werden. Obwohl ein linearer Verlauf sehr unwahrscheinlich erscheint, können mit dem linearen Modell bereits r round(summary(m_l)$r.squared*100, 2)% der Variation der Bevölkerungsdichte durch die Jahreszahl eklärt werden. Wie der Grafik deutlich zu entnehmen ist, sind die Residuen in dieser Regressionanalyse stark abhängig von der Jahreszahl (negatives Residuum von a. 1860-1970 und positive Residuen sonst; Wiederholung: $\varepsilon_i=Y_i-\hat{Y}_i$, wobei $\hat{Y}_i$ der vorhergesagte Wert ist). Auch wenn wir uns das zugehörige Histogramm der Residuen ansehen, widerspricht dieses der Annahme auf Normalverteilung.
Quadratisches Modell für das Bevölkerungswachstum
Vielleicht ist auch hier ein quadratischer Effekt versteckt? Wir möchten dem auf den Grund gehen und nutzen wieder die Funktion poly, um ein Polynom 2. Grades (eine quadratische Funktion) der Jahreszahl in unsere Analysen mit aufzunehmen.
m_q <- ... # quadratischer Verlauf
summary(...)
# Finale Lösung
m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf
summary(m_q)
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf
m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf
m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf
Das Varianzinkrement beläuft sich auf:
summary(m_q)$r.squared - summary(m_l)$r.squared # Inkrement
Testen Sie dieses auf Signifikanz mit anova.
# Intro into excercise
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf
m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf
# Begin grading:
res <- anova(m_l, m_q)
sol2 <- anova(m_q, m_l)
grade_result(
fail_if(~ (!identical(.result, res) & !identical(.result, sol2)), 'Der Befehl anova muss in diesem Zusammenhang zwei Argumente entgegennehmen: die zwei Regressions-Objekte, die wir zuvor abgespeichert hatten. Schauen Sie in der Übung zuvor an, wie dieses Modell hieß, falls Sie sich nicht sicher sind.'),
fail_if(~ (identical(.result, sol2)), 'Das ist an sich RICHTIG, allerdings sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell (mit weniger Prädiktoren und Parametern, die zu schätzen sind) zuerst übergeben werden. Hier: m_l, da sonst die df negativ sind und auch die Änderung in den Sum of Sq (Quadratsumme) negativ sind! `R` erkennt dies zwar und testet trotzdem die richtige Differenz auf Signifikanz, aber wir wollen uns besser vollständig korrekt aneignen!'),
pass_if(~ identical(.result, res)),
correct = 'Sehr gut! Dies ist der Modellvergleich des linearen und des quadratischen Modells. Hier sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell (mit weniger Prädiktoren und Parametern, die zu schätzen sind) zuerst übergeben werden. Hier: m_l, da sonst die df negativ sind und auch die Änderung in den Sum of Sq (Quadratsumme) negativ sind! `R` erkennt dies zwar und testet trotzdem die richtige Differenz auf Signifikanz, aber wir wollen uns besser vollständig korrekt aneignen!',
incorrect = 'Leider falsch.',
glue_correct = '{.correct}',
glue_incorrect = '{.incorrect} {.message}')
Auch dem geom_smooth Befehl kann einfach das modifizierte Regressionsmodell übergeben werden. Stellen Sie die Weltbevölkerung sowie das lineare und das quadratische Modell dar und färben Sie die Kurve in dunkelblau (col = "darkblue").
ggplot(data = WorldPopulation, aes(x = Year, y = Population))+
geom_point()+geom_smooth(method="lm", formula = "y~x")+ # plotte linearen Verlauf (bis hier her ist es Code der Aufgabe zuvor!)
... # plotte quadratischen Verlauf
# Finale Lösung
ggplot(data = WorldPopulation, aes(x = Year, y = Population))+
geom_point()+geom_smooth(method="lm", formula = "y~x")+ # plotte linearen Verlauf
geom_smooth(method="lm", formula = "y~poly(x,2)", col = "darkblue") # plotte quadratischen Verlauf
Durch den quadratischen Verlauf lassen sich r round(summary(m_q)$r.squared*100, 2)% der Variation der Bevölkerungsdichte erklären, was einem signifikantem Varianzinkrement von r round(summary(m_q)$r.squared*100 - summary(m_l)$r.squared*100, 2)% entspricht (mit einer Irrtumswahrscheinlichkeit von 5% ist das Inkrement in der Population nicht null. Dies ist äquivialent zu folgdender Aussage: mit einer Irrtumswahrscheinlichkeit von 5% ist der Effektparameter (der Regressionskoeffizient) des quadratischen Verlaufs in der Population nicht null; dies spricht folglich für einen quadratischen im Gegensatz zu einem linearen Verlauf). Der Grafik ist deutlich zu entnehmen, dass der quadratische Verlauf nicht weit vom empirischen entfernt liegt.
Exponentielles Modell für das Bevölkerungswachstum
Nun wollen wir prüfen, ob ein exponentieller Verlauf die Daten noch besser beschreibt. Dazu müssen wir zunächst die Weltpopulation logarithmieren und können anschließend ein lineares Modell verwenden. Transformieren Sie hierzu die Variablen und speichern diese als weitere Spalte in unserem Datensatz ab. Ordnen Sie diese folgendem Namen zu: WorldPopulation$log_Population <- ... und geben Sie anschließend WorldPopulation$log_Population aus. Der Name ist entscheidend, da wir diesen später weiter verwenden wollen.
# Intro into excercise
# Begin grading:
res <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung
grade_result(
fail_if(~ (!identical(.result, res)), 'Sie müssen den Befehl `log` verwenden.'),
pass_if(~ identical(.result, res)),
correct = 'Sehr gut! Sie haben die logarithmierte Weltbevölkerung unserem Datensatz hinzugefügt.',
incorrect = 'Leider falsch.',
glue_correct = '{.correct}',
glue_incorrect = '{.incorrect} {.message}')
WorldPopulation$log_Population <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf
m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf
WorldPopulation$log_Population <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung
Wir wollen anschließend ein einfaches Regressionsmodell schätzen, in welchem die logarithmierte Bevölkerungszahl die abhängige Variable darstellt. Speichern Sie das lm Objekt unter dem Namen m_log ab. Geben Sie die summary für m_log aus.
m_log <- lm(... ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
m_log <- lm(log_Population ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
summary(m_log)
m_log <- lm(log_Population ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf
m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf
WorldPopulation$log_Population <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung
m_log <- lm(log_Population ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
Das lineare Modell für die logarithmierte Bevölkerunsdichte scheint gut zu den Daten zu passen. Insgesamt können r round(summary(m_log)$r.squared*100, 2) % der Variation in den Daten durch den Zeitverlauf erklärt werden; etwas weniger als durch das quadratische Wachstum, durch welches r round(summary(m_q)$r.squared*100, 2) % der Variation in den Daten durch den Zeitverlauf erklärt werden konnten. Die zugehörige Grafik sieht so aus:
ggplot(data = WorldPopulation, aes(x = Year, y = log_Population))+
geom_point()+geom_smooth(method="lm", formula = "y~x", col = "red")+
labs(title = "Logarithmierte Weltbevölkerung vs. Jahr")
Wieso erklärt nun das quadratische Modell mehr Variation? Wir wollen die Vorhersage auch retransformiert plotten, um sie mit dem quadratischen Verlauf grafisch vergleichen zu können. Dazu müssen wir die vorhergesagten Werte unseres logarithmischen Modells verwenden und mit mit der Umkehrfunktion des Logarithmus retransformieren. Diese können wir einem lm-Objekt mit predict entlocken.
WorldPopulation$pred_Pop_exp <- ... # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung)
WorldPopulation
# Verwenden Sie
predict(...)
# Verwenden Sie
predict(...) # und transformieren Sie mit
exp(...)# welche die Umkehrfunktion des Log ist
WorldPopulation$pred_Pop_exp <- exp(predict(m_log)) # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung)
WorldPopulation
WorldPopulation$pred_Pop_exp <- exp(predict(m_log)) # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung)
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf
m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf
WorldPopulation$log_Population <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung
m_log <- lm(log_Population ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
WorldPopulation$pred_Pop_exp <- exp(predict(m_log)) # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung)
Dem Datensatz haben wir nun eine neue Spalte hinzugefügt, welche die vorhergesagten Populationswerte enthält (retransformiert; nicht mehr in Log-Skala).
ggplot(data = WorldPopulation, aes(x = Year, y = Population))+
geom_point()+geom_smooth(method="lm", formula = "y~x")+ # plotte linearen Verlauf
geom_smooth(method="lm", formula = "y~poly(x,2)", col = "darkblue")+ # plotte quadratischen Verlauf
geom_line(aes(x = Year, y = ...), col = "red", lwd = 1.5)
ggplot(data = WorldPopulation, aes(x = Year, y = Population))+
geom_point()+geom_smooth(method="lm", formula = "y~x")+ # plotte linearen Verlauf
geom_smooth(method="lm", formula = "y~poly(x,2)", col = "darkblue")+ # plotte quadratischen Verlauf
geom_line(aes(x = Year, y = pred_Pop_exp), col = "red", lwd = 1.5)
Das Diagramm der retransformierten vorhergesagten Werten signalisiert, dass ein exponentielles Wachstumsmodell die Daten gut beschreibt, allerdings scheint der quadratische Trend vor allem ab ca. 1975 die Daten besser zu beschreiben. Wir schauen uns dazu die Residuen des logarithmischen Modells an.
residualPlot(..., col = "red") # Residualplot
# Finale Lösung
residualPlot(m_log, col = "red") # Residualplot
Dem Residualplot ist zu entnehmen, dass ein nicht-linearer Verlauf angemessen wäre. Wir gehen dieser Vermutung nach, indem wir dem Log-Plot einen quadratischen Verlauf hinzufügen.
ggplot(data = WorldPopulation, aes(x = Year, y = log_Population))+
geom_point()+geom_smooth(method="lm", formula = "y~x", col = "red")+
geom_smooth(method="lm", formula = "y~poly(x,2)", col = "gold3")+
labs(title = "Logarithmierte Weltbevölkerung vs. Jahr")
Vielleicht wächst die Bevölkerung sogar schneller als exponentiell? Die gelbe Linie im logarithmierten Plot lässt dies vermuten. Um dies genauer zu untersuchen fügen wir auch in das logarithmierte Modell einen quadratischen Trend der Zeit ein.
Quadratisch-exponentielles Modell für das Bevölkerungswachstum
Dies wollen wir schätzen, indem wir im Regressionsmodell der logarithmierten Populationsgröße einen quadratischen Trend (poly(...,2)) für die Zeit definieren. Geben Sie für dieses Regressionsmodell die summary aus.
m_log_quad <-
summary(...)
m_log_quad <- lm(log_Population ~ poly(Year, 2), data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
summary(m_log_quad)
m_log_quad <- lm(log_Population ~ poly(Year, 2), data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf
m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf
WorldPopulation$log_Population <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung
m_log <- lm(log_Population ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
WorldPopulation$pred_Pop_exp <- exp(predict(m_log)) # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung)
m_log_quad <- lm(log_Population ~ poly(Year, 2), data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
Testen Sie dieses Modell im Vergleich zu m_log auf Signifikanz mit anova .
# Intro into excercise
WorldPopulation$log_Population <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung
m_log <- lm(log_Population ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
m_log_quad <- lm(log_Population ~ poly(Year, 2), data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
# Begin grading:
res <- anova(m_log, m_log_quad)
sol2 <- anova(m_log_quad, m_log)
grade_result(
fail_if(~ (!identical(.result, res) & !identical(.result, sol2)), 'Der Befehl anova muss in diesem Zusammenhang zwei Argumente entgegennehmen: die zwei Regressions-Objekte, die wir zuvor abgespeichert hatten. Schauen Sie in der Übung zuvor an, wie dieses Modell hieß, falls Sie sich nicht sicher sind.'),
fail_if(~ (identical(.result, sol2)), 'Das ist an sich RICHTIG, allerdings sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell (mit weniger Prädiktoren und Parametern, die zu schätzen sind) zuerst übergeben werden. Hier: m_log, da sonst die df negativ sind und auch die Änderung in den Sum of Sq (Quadratsumme) negativ sind! `R` erkennt dies zwar und testet trotzdem die richtige Differenz auf Signifikanz, aber wir wollen uns besser vollständig korrekt aneignen!'),
pass_if(~ identical(.result, res)),
correct = 'Sehr gut! Dies ist der Modellvergleich des linearen und des quadratischen logarithmierten Modells. Hier sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell (mit weniger Prädiktoren und Parametern, die zu schätzen sind) zuerst übergeben werden. Hier: m_log, da sonst die df negativ sind und auch die Änderung in den Sum of Sq (Quadratsumme) negativ sind! `R` erkennt dies zwar und testet trotzdem die richtige Differenz auf Signifikanz, aber wir wollen uns besser vollständig korrekt aneignen!',
incorrect = 'Leider falsch.',
glue_correct = '{.correct}',
glue_incorrect = '{.incorrect} {.message}')
Durch den quadratisch-exponentiellen Verlauf (bzw. den quadratischen Verlauf in den logarithmierten Daten) lassen sich r round(summary(m_log_quad)$r.squared*100, 2)% der Variation der Bevölkerungsdichte erklären, was einem signifikantem Varianzinkrement von r round(summary(m_log_quad)$r.squared*100 - summary(m_log)$r.squared*100, 2)% im Vergleich zum reinen exponentiellen Verlaufsmodell entspricht. Interessant zu sehen ist, dass fast 100% der Variation im Datensatz erklärbar ist. Eine Übersicht über $R^2$ in den vier Modellen ist in Appendix B einzusehen.
Stellen wir dies grafisch dar. Dazu müssen wir wieder die Vorhersage abspeichern. Diese können wir einem lm-Objekt mit predict entlocken.
WorldPopulation$pred_Pop_exp_quad <- ... # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung)
WorldPopulation
# Verwenden Sie
predict(...)
# Verwenden Sie
predict(...) # und transformieren Sie mit
exp(...)# welche die Umkehrfunktion des Log ist
WorldPopulation$pred_Pop_exp_quad <- exp(predict(m_log_quad)) # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung)
WorldPopulation
WorldPopulation$pred_Pop_exp_quad <- exp(predict(m_log_quad))
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf
m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf
WorldPopulation$log_Population <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung
m_log <- lm(log_Population ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
WorldPopulation$pred_Pop_exp <- exp(predict(m_log)) # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung)
m_log_quad <- lm(log_Population ~ poly(Year, 2), data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
WorldPopulation$pred_Pop_exp_quad <- exp(predict(m_log_quad))
Dem Datensatz haben wir nun eine neue Spalte hinzugefügt, welche die vorhergesagten Populationswerte enthält, die durch das quadratisch-exponentielle Modell vorhergesagt werden (retransformiert; nicht mehr in Log-Skala). Im Folgenden ist die finale Grafik mit dem linearen, quadratischen, exponentiellen und dem quadratisch-exponentiellen Verlauf dargestellt. Außerdem schauen wir uns das Histogramm sowie den Residuenplot des quadratisch-exponentiellen Modells an (in logarithmierter Skala, in welcher das Modell auch geschätzt wurde!).
ggplot(data = WorldPopulation, aes(x = Year, y = Population))+
geom_point()+geom_smooth(method="lm", formula = "y~x")+ # plotte linearen Verlauf
geom_smooth(method="lm", formula = "y~poly(x,2)", col = "darkblue")+ # plotte quadratischen Verlauf
geom_line(aes(x = Year, y = pred_Pop_exp), col = "red", lwd = 1.5)+
geom_line(aes(x = Year, y = pred_Pop_exp_quad), col = "gold3", lwd = 2)
# nur quadratisch-exponentiell
ggplot(data = WorldPopulation, aes(x = Year, y = Population))+
geom_point()+
geom_line(aes(x = Year, y = pred_Pop_exp_quad), col = "gold3", lwd = 2)+
labs(title = "Beobachtetes und durch das quadratisch-exponentielle Modell\n vorhergesagtes Bevölkerungswachstum")
res <- studres(m_log_quad) # Studentisierte Residuen als Objekt speichern
df_res <- data.frame(res) # als Data.Frame für ggplot
# Grafisch: Histogramm mit Normalverteilungskurve
library(ggplot2)
ggplot(data = df_res, aes(x = res)) +
geom_histogram(aes(y =..density..),
bins = 10, # Wie viele Balken sollen gezeichnet werden?
colour = "blue", # Welche Farbe sollen die Linien der Balken haben?
fill = "skyblue") + # Wie sollen die Balken gefüllt sein?
stat_function(fun = dnorm, args = list(mean = mean(res), sd = sd(res)), col = "darkblue") + # Füge die Normalverteilungsdiche "dnorm" hinzu und nutze den empirischen Mittelwert und die empirische Standardabweichung "args = list(mean = mean(res), sd = sd(res))", wähle dunkelblau als Linienfarbe
labs(title = "Histogramm der Residuen mit Normalverteilungsdichte\n für das logarithmierte quadratische Modell", x = "Residuen") # Füge eigenen Titel und Achsenbeschriftung hinzu
residualPlot(m_log_quad)
Es ist deutlich zu sehen, dass der Knick um 1950 durch das quadratisch-exponentielle Wachstum am besten darzustellen/ zu modellieren ist.
Das Histogramm der Residuen des logarithmierten quadratischen Modells (welches das quadratisch-exponentielle Wachstum der Bevölkerung darstellt) ist etwas rechts-schief/links-steil. Die Residuenplots zeigen außerdem, dass auch hier die Residuen nicht vollständig unsystematisch sind. Dennoch ist die Resiudalvarianz sehr klein.
Da die Modelle des quadratischen und des quadratisch-exponentiellen Verlaufs nicht geschachtelt sind (sie gehen nicht auseinander hervor und sind nicht durch Restringieren, bzw. Nullsetzen, von Regressionsparametern aus einander überführbar), können sie nicht inferenzstatistisch miteinander verglichen werden. Da das quadratisch-exponentielle Modell mehr Variation aufklärt als das rein quadratische Modell und auch grafisch besser zum Verlauf der Daten passt (besonders von 1800 bis ca. 1850 sagt das quadratische Modell ein Verkleinern der Bevölkerung vorher), entscheiden wir uns final für dieses.
Parameterinterpretation des quadratisch-exponentielles Modell für das Bevölkerungswachstum
Was bedeuten nun die Parameter in unserem quadratisch-exponentiellen Modell? Der Regressionskoeffizient des linearen Trends von poly liegt bei r round(coef(m_log_quad)[2], 2) und der Koeffizient des quadratischen Trends bei r round(coef(m_log_quad)[3], 2). Dies spricht für ein beschleunigtes exponentielles Wachstum. Nach diesem Vorhersagemodel scheinen die Menschen sich schneller als exponentiell zu vermehren (zumindest von 1800 bis 2020). Eine Übersicht über die Modelle sehen Sie in Appendix C.
Wären beide Koeffizienten negativ, so würde dies für beschleunigten exponentiellen Zerfall/Abnahme sprechen. Ist der Koeffizient des quadratische Trends negativ, wird von exponentiellem Wachstum mit Dämpfung gesprochen.
Final ist zu sagen, dass die Menschheit am wahrscheinlichsten exponentiell wächst. Ein quadratisches Modell ist im Allgemeinen nicht sinnvoll. Da wir uns in dieser Übung nur einen kleine Zeitausschnitt angesehen haben, kam es zu diesen Abweichungen. Allerdings würde das geschätzte quadratische Modell für das Jahr 0 eine Bevölkerungszahl von r round(predict(object = m_q, newdata = data.frame(Year=c(0)))/10^11, 2)$*10^{11}$ vorhersagen, was absolut nicht sinnvoll ist.
Aufgaben
Die Aufgaben zu dieser Sitzung finden Sie im Appendix D.
Appendix A {#AppendixA}
Exkurs: Was genau macht poly?
X <- 1:10 # Variable X
X2 <- X^2 # Variable X hoch 2
X_poly <- poly(X, 2) # erzeuge Variable X und X hoch mit Hilfe der poly Funktion
colnames(X_poly) <- c("poly(X, 2)1", "poly(X, 2)2")
cbind(X, X2, X_poly)
Die Funktion poly erzeugt sogenannte orthogonale Polynome. Das bedeutet, dass zwar die $x$ und $x^2$ berechnet werden, diese Terme anschließend allerdings so transformiert werden, dass sie jeweils einen Mittelwert von 0 und die gleiche Varianz haben und zusätzlich noch unkorreliert sind:
round(apply(X = cbind(X, X2, X_poly), MARGIN = 2, FUN = mean), 2) # Mittelwerte über die Spalten hinweg berechnen
round(apply(X = cbind(X, X2, X_poly), MARGIN = 2, FUN = sd), 2) # Standardabweichung über die Spalten hinweg berechnen
round(cor(cbind(X, X2, X_poly)),2) # Korrelationen berechnen
Die Funktion apply führt an der Matrix, welche dem Argument X übergeben wird, entweder über die Zeilen MARGIN = 1 oder über die Spalten MARGIN = 2 (hier jeweils gewählt) die Funktion aus, welche im Argument FUN angegeben wird. So wird zunächst mit FUN = mean der Mittelwert und anschließend mit FUN = sd die Standardabweichung von $X, X^2$ sowie poly(X, 2) berechnet. Der Korrelationsmatrix ist zu entnehmen, dass $X$ und $X^2$ in diesem Beispiel sehr hoch miteinander korrelieren und somit gleiche lineare Informationen enthalten ($\hat{r}_{X,X^2}$ = cor(X, X2) = 0.97), während die linearen und die quadratischen Anteile in poly(X, 2) keinerlei lineare Gemeinsamkeiten haben - sie sind unkorreliert (cor(poly(X,2)1 , poly(X,2)2) = 0). Ein weiterer Vorteil ist deshalb, dass bei sukzessiver Aufnahme der Anteile von poly(X, 2) in ein Regressionmodell, sich die Parameterschätzungen des linearen Terms im Modell nicht (bzw. sehr wenig) ändern. Der Anteil erklärter Varianz bleibt jedoch in allen gleich - die Modelle sind äquivalent, egal auf welche Art und Weise quadratische Terme gebildet werden.
m1.b1 <- lm(Reading ~ HISEI + poly(MotherEdu, 1) + Books, data = PISA2009)
summary(lm.beta(m1.b1))
m1.b2 <- lm(Reading ~ HISEI + poly(MotherEdu, 2) + Books, data = PISA2009)
summary(lm.beta(m1.b2))
PISA2009$MotherEdu2 <- PISA2009$MotherEdu^2 # füge dem Datensatz den quadrierten Bildungsabschluss der Mutter hinzu
m1.c1 <- lm(Reading ~ HISEI + MotherEdu + Books, data = PISA2009)
summary(lm.beta(m1.c1))
m1.c2 <- lm(Reading ~ HISEI + MotherEdu + MotherEdu2 + Books, data = PISA2009)
summary(lm.beta(m1.c2))
rbind(coef(m1.b1), coef(m1.c1))
rbind(coef(m1.b2),coef(m1.c2))
Einordnung quadratischer Verläufe
Wie kommen wir nun auf die Interpretation der quadratischen Beziehung?
Eine allgemeine quadratische Funktion $f$ hat folgende Gestalt
$$f(x):=ax^2 + bx + c,$$
wobei $a\neq 0$, da es sich sonst nicht um eine quadratische Funktion handelt. Wäre $a=0$, würde es sich um eine lineare Funktion mit Achsenabschnitt $c$ und Steigung (Slope) $b$ handeln. Wäre zusätzlich $b=0$, so handelt es sich um eine horizontale Linie bei $y=f(x)=c$.
Für betraglich große $x$ fällt $x^2$ besonders ins Gewicht. Damit entscheidet das Vorzeichen von $a$, ob es sich um eine u-formige (falls $a>0$) oder eine umgekehrt-u-förmige (falls $a<0$) Beziehung handelt. Die betragliche Größe von $a$ entscheidet hierbei, wie gestaucht die u-förmige Beziehung (die Parabel) ist. Die reine quadratische Beziehung $f(x)=x^2$ sieht so aus:
a <- 1; b <- 0; c <- 0
x <- seq(-2,2,0.01)
f <- a*x^2 + b*x + c
data_X <- data.frame(x, f)
ggplot(data = data_X, aes(x = x, y = f)) +
geom_line(col = "black")+
labs(x = "x", y = "f(x)",
title = expression("f(x)="~x^2))
Wir werden diese Funktion immer als Referenz mit in die Grafiken einzeichnen.
a <- 0.5; b <- 0; c <- 0
x <- seq(-2,2,0.01)
f <- a*x^2 + b*x + c
data_X <- data.frame(x, f)
ggplot(data = data_X, aes(x = x, y = f)) +
geom_line(col = "blue", lwd = 2)+ geom_line(mapping = aes(x, x^2), lty = 3)+
labs(x = "x", y = "f(x)",
title = expression("f(x)="~0.5*x^2))
a <- 2; b <- 0; c <- 0
x <- seq(-2,2,0.01)
f <- a*x^2 + b*x + c
data_X <- data.frame(x, f)
ggplot(data = data_X, aes(x = x, y = f)) +
geom_line(col = "blue", lwd = 2)+ geom_line(mapping = aes(x, x^2), lty = 3)+
labs(x = "x", y = "f(x)",
title = expression("f(x)="~2*x^2))
a <- -1; b <- 0; c <- 0
x <- seq(-2,2,0.01)
f <- a*x^2 + b*x + c
data_X <- data.frame(x, f)
ggplot(data = data_X, aes(x = x, y = f)) +
geom_line(col = "blue", lwd = 2)+ geom_line(mapping = aes(x, x^2), lty = 3)+
labs(x = "x", y = "f(x)",
title = expression("f(x)="~-x^2))
Diese invers-u-förmige Beziehung ist eine konkave Funktion. Als Eselsbrücke für das Wort konkav, welches fast das englische Wort cave enthält, können wir uns merken: eine konkave Funktion stellt eine Art Höhleneingang dar.
$c$ bewirkt eine vertikale Verschiebung der Parabel:
a <- 1; b <- 0; c <- 1
x <- seq(-2,2,0.01)
f <- a*x^2 + b*x + c
data_X <- data.frame(x, f)
ggplot(data = data_X, aes(x = x, y = f)) +
geom_line(col = "blue", lwd = 2)+ geom_line(mapping = aes(x, x^2), lty = 3)+
labs(x = "x", y = "f(x)",
title = expression("f(x)="~x^2+1))
$b$ bewirkt eine horizontale und vertikale Verschiebung, die nicht mehr leicht vorhersehbar ist. Für $f(x)=x^2+x$ lässt sich beispielsweise durch Umformen leicht erkennen: $f(x)=x^2+x=x(x+1)$, dass diese Funktion zwei Nullstellen bei $0$ und $-1$ hat. Somit ist ersichtlich, dass die Funktion nach unten und links verschoben ist:
a <- 1; b <- 1; c <- 0
x <- seq(-2,2,0.01)
f <- a*x^2 + b*x + c
data_X <- data.frame(x, f)
ggplot(data = data_X, aes(x = x, y = f)) +
geom_line(col = "blue", lwd = 2)+ geom_line(mapping = aes(x, x^2), lty = 3)+
labs(x = "x", y = "f(x)",
title = expression("f(x)="~x^2+x))
Für die genaue Gestalt einer allgemeinen quadratischen Funktion $ax^2 + bx + c$ würden wir die Nullstellen durch Lösen der Gleichung $ax^2 + bx + c=0$ bestimmen (via p-q Formel oder a-b-c-Formel). Den Scheitelpunkt würden wir durch Ableiten und Nullsetzen der Gleichung bestimmen. Wir müssten also $2ax+b=0$ lösen und dies in die Gleichung einsetzen. Wir könnten auch die binomischen Formeln nutzen, um die Funktion in die Gestalt $f(x):=a'(x-b')^2+c'$ oder $f(x):=a'(x-b'_1)(x-b_2')+c'$ zu bekommen, falls die Nullstellen reell sind (also das Gleichungssystem lösbar ist), da wir so die Nullstellen ablesen können als $b'$ oder $b_1'$ und $b_2'$, falls $c=0$. Für die Interpretation der Ergebnisse reicht es zu wissen, dass $a$ eine Stauchung bewirkt und entscheind dafür ist, ob die Funktion u-förmig oder invers-u-förmig verläuft.
a <- -0.5; b <- 1; c <- 2
x <- seq(-2,2,0.01)
f <- a*x^2 + b*x + c
data_X <- data.frame(x, f)
ggplot(data = data_X, aes(x = x, y = f)) +
geom_line(col = "blue", lwd = 2)+ geom_line(mapping = aes(x, x^2), lty = 3)+
labs(x = "x", y = "f(x)",
title = expression("f(x)="~-0.5*x^2+x+2))
$\longrightarrow$ so ähnlich sieht die bedingte Beziehung (kontrolliert für die weiteren Prädiktoren im Modell) zwischen Bildungsabschluss der Mutter und Leseleistung aus.
Code für quadratische Verlaufsgrafik
Der Code, der die Grafik des standardisierten vorhergesagten bedingten Verlaufs des Bildungsabschlusses der Mutter erzeugt, sieht folgendermaßen aus:
X <- scale(poly(PISA2009$MotherEdu, 2))
std_par_ME <- c(0.1588, -0.1436)
pred_effect_ME <- X %*% std_par_ME
std_ME <- X[,1]
data_ME <- data.frame(std_ME, pred_effect_ME)
ggplot(data = data_ME, aes(x = std_ME, y = pred_effect_ME)) + geom_point(pch = 16, col = "blue", cex = 4)+
labs(y = "std. Leseleistung | Others", x = "std. Bildungsabschluss der Mutter | Others",
title = "Standardisierte bedingte Beziehung zwischen\n Bildungsabschluss der Mutter und Leseleistung")
Wir verwenden scale, um die linearen und quadratischen Anteile des Bildungsabschlusses der Mutter zu standardisieren und speichern sie in X. Anschließend ist das Interzept der quadratischen Funktion 0 ($c=0$, da wir standardisiert haben). Die zugehörigen standardisierten Koeffizienten sind $b=0.1588$ und $a=-0.1436$, die wir aus der standardisierten summary abgelesen haben. Somit wissen wir, dass es sich um eine invers-u-förmige Beziehung handelt (ohne die Grafik zu betrachten). Wir speichern die standardisierten Koeffizienten unter std_par_ME ab und verwenden anschließend das Matrixprodukt X %*% std_par_ME, um die vorhergesagten Werte via $y_{std,i}=0.1588 ME - 0.1436ME^2$ zu berechnen. Diese vorhergesagten Werte pred_effect_ME plotten wir nun gegen die standardisierten Werte des Bildungsabschlusses der Mutter std_ME, welche in der ersten Spalte von X stehen: X[, 1].
Appendix B {#AppendixB}
Übersicht über erklärte Varianzanteile
Hier sind nochmals die Anteile erklärter Varianz der Bevölkerungsdichte über die Zeit in den vier betrachteten Modellen dargestellt:
R2 <- rbind(summary(m_l)$r.squared,
summary(m_q)$r.squared,
summary(m_log)$r.squared,
summary(m_log_quad)$r.squared)
rownames(R2) <- c("linear", "quadratisch", "exponentiell (log. Modell)", "quadratisch-exponentiell (quadratisches log. Modell)")
colnames(R2) <- "R^2"
round(R2, 4)
Appendix C {#AppendixC}
Übersicht über die Modelle
Die angenommenen Modell pro Messzeitpunkt $i$ sind von der Konzeption deutlich verschieden. Insbesondere der Regressionfehler ist an einer anderen Stelle:
-
Lineares Modell: $Y_i = \beta_0 + \beta_1t_i + \varepsilon_i$
-
Quadratisches Modell: $Y_i = \beta_0 + \beta_1t_i^ + \beta_2t_i^{2} + \varepsilon_i$. Hier wurde mit poly(...,2) eine Transformation der Variable Zeit vorgenommen, was hier durch "$^*$" dargestellt werden soll.
-
Logarithmisches Modell: $\text{ln}(Y_i) = \beta_0 + \beta_1t_i + \varepsilon_i$ bzw. retransformiert $Y_i = e^{\beta_0 + \beta_1t_i + \varepsilon_i} = e^{\beta_0 + \beta_1t_i} \ e^{\varepsilon_i}$. Diesem Modell ist zu entnehmen, dass sich der Fehler in Abhängigkeit von der Ausprägung der unabhängigen Variablen (hier: $t_i$) unterschiedlich stark auf das Kriterium (die abhängige Variable, hier die Bevölkerungszahl) auswirkt - er hängt also mit der Ausprägung der unabhängigen Variable zusammen. Dies ist im logarithmierten Modell nicht mehr der Fall; es handelt sich hier bis auf die Transformation der abhängigen Variable um ein "normales" Regressionsmodell: Hier verliert sich diese Beziehung zwischen unabhängiger Variable und Fehler!
-
Logarithmisches Modell mit quadratischem Term: $\text{ln}(Y_i) = \beta_0 + \beta_1t_i^ + \beta_2t_i^{2} + \varepsilon_i$ bzw. retransformiert $Y_i = e^{\beta_0 + \beta_1t_i^ + \beta_2t_i^{2} + \varepsilon_i}$. Hier wurde mit poly(...,2) eine Transformation der Variable Zeit vorgenommen, was hier durch "$^*$" dargestellt werden soll.
Appendix D {#AppendixD}
Aufgaben
Corona-Datensatz (wie in der ggplot2-Sitzung; PsyBSc7::Sitzung_2())
Zur Vorbereitung laden Sie die Daten ein. Entweder via folgendem Code
library(ggplot2)
library(car)
library(MASS)
### Datensatz: Corona-Pandemie 2020
confirmed <- read.csv('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv')
deaths <- read.csv('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv')
### Long format ---
confirmed_long <- reshape(confirmed,
varying = names(confirmed)[-c(1:4)],
v.names = 'Confirmed',
timevar = 'Day',
idvar = names(confirmed)[1:4],
direction = 'long')
deaths_long <- reshape(deaths,
varying = names(deaths)[-c(1:4)],
v.names = 'Deaths',
timevar = 'Day',
idvar = names(deaths)[1:4],
direction = 'long')
### Merged data ----
long <- merge(confirmed_long, deaths_long,
by = c('Province.State', 'Country.Region', 'Lat', 'Long', 'Day'))
### Full data ----
covid <- aggregate(cbind(Confirmed, Deaths) ~ Country.Region + Day, data = long, FUN = 'sum')
### Only data until April 11th ----
covid_full <- covid
covid <- covid[covid$Day < 80, ]
### Subsets ----
covid_de <- covid[covid$Country.Region == 'Germany', ]
covid_sel <- covid[covid$Country.Region %in% c('France', 'Germany', 'Italy', 'Spain', 'United Kingdom'), ]
oder über PsyBSc7::Aufgaben_8(). Die tagesaktuellen Daten finden Sie unter covid_full. Wir wollen uns für die Aufgaben jedoch zeitlich einschränken bis zum 11. April. Betrachten Sie in den folgenden Aufgaben bitte immer nur die Daten bis zum 11. April in covid_sel und covid_de.
Aufgabe 1. Verschaffen Sie sich einen Überblick über die Daten:
a) Stellen Sie die Verläufe der bestätigten Fälle in den Ländern Frankreich, Deutschland, Italien, Spanien und Großbritannien gegen die Zeit dar (Subdatensatz covid_sel).
b) Stellen Sie nur den Verlauf der bestätigten Fälle von Deutschland grafisch dar (Subdatensatz covid_de).
c) Stellen Sie den logarithmierten Verlauf der bestätigten Fälle von Deutschland grafisch dar.
Tipp: der Logarithmus von 0 ist nicht definiert (log(0) = -Inf); diese Werte machen in einer Regressionsanalyse keinen Sinn, weswegen Sie diese durch fehlende Werte (NA) ersetzten müssen (Subdatensatz covid_de)
d) Löschen Sie die Tage aus dem Datensatz heraus, an welchen noch keine bestätigten Fälle in Deutschland aufgetreten sind.
Tipp: is.na(...) fragt ab, ob Einträge in einem Element fehlen. Entsprechend können mit !is.na(...) jene Fälle angesprochen werden, die vorhanden sind. (Subdatensatz covid_de)
Aufgabe 2. Schauen Sie sich den Verlauf der bestätigten Infektion in Deutschland genauer an, indem Sie die folgenden Regressionsmodelle berechnen.
a) Berechnen Sie eine lineare und quadratische Regression und vergleichen Sie die beiden Modelle. (Subdatensatz covid_de)
b) Führen Sie eine lineare Regression der logarithmierten "Confirmed" Cases durch. (Subdatensatz covid_de)
c) Passt das lineare logarithmierte Modell (bzw. das exponentielle Modell) zu den Daten? Analysieren Sie die Voraussetzungen. (Subdatensatz covid_de)
d) Stellen Sie den logarithmierten sowie den retransformierten Verlauf dar. (Subdatensatz covid_de)
Aufgabe 3. Prüfen Sie nun, ob Sie die Vorhersage durch ein quadratisch-exponentielles Modell verbessern können, indem Sie einen quadratischen Trend der Zeit aufnehmen und die Modelle vergleichen.
(Subdatensatz covid_de)
Bonusaufgabe.: Vergleichen Sie die Koeffizienten des exponentiellen Wachstums zwischen den Ländern.
(Subdatensatz covid_sel)
(Die Bonusaufgabe ist komplett freiwillig und muss nicht gelöst werden, um die Studienleistung zu erbringen.)
Literatur
Eid, M., Gollwitzer, M., & Schmitt, M. (2017). Statistik und Forschungsmethoden (5. Auflage, 1. Auflage: 2010). Weinheim: Beltz.
- Blau hinterlegte Autorenangaben führen Sie direkt zur universitätsinternen Ressource.
martscht/PsyBSc7 documentation built on Sept. 1, 2020, 10:50 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
source('setup.R') library(learnr) # package to generate interactive HTML library(gradethis) # grade Code from input in HTML library(shiny) # nice appearance in R library(fontawesome) # nice fonts library(car) library(MASS) library(ggplot2) library(lm.beta) # erforderlich für standardiserte Gewichte data("PISA2009", package = "PsyBSc7") data("WorldPopulation", package = "PsyBSc7") knitr::opts_chunk$set(exercise.checker = gradethis::grade_learnr)
Einleitung und Datensatz (PISA: wie in den Übungen zum letzten Termin)
$*$ zum Teil basierend auf Unterlagen von Prof. Johannes Hartig
In dieser Sitzung werden wir uns weitere nichtlineare Effekte in Regressionsmodellen ansehen. Dazu verwenden wir zunächst den Datensatz aus der Übung zur PsyBSc7::Sitzung_7, welcher in den Aufgaben verwendet wurde und dessen Analysen in PsyBSc7::Loesung_7 mit dem richtigen Passwort einsehbar waren.
Der Beispieldatensatz enthält Daten zu Lesekompetenz aus der deutschen Stichprobe der PISA-Erhebung in Deutschland 2009.
library(car) library(MASS) library(lm.beta) # erforderlich für standardiserte Gewichte library(ggplot2) # Datensatz laden data("PISA2009", package = "PsyBSc7")
Quadratische Verläufe in der Vorhersage von Lesekompetenz mit individuellen Merkmalen der Schüler/innen
$**$ zum Teil basierend auf Unterlagen von Prof. Johannes Hartig
In der Übung zur letzten Sitzung hatten wir herausgefunden, dass der Sozialstatus (HISEI), der Bildungsabschluss der Mutter (MotherEdu) und die Zahl der Bücher zu Hause (Books) bedeutsame Prädiktoren für die Lesekompetenz der Schüler/innen sind, allerdings zeigten Analysen, dass nicht alle Voraussetzungen erfüllt waren:
# Berechnung des Modells und Ausgabe der Ergebnisse m1 <- lm(Reading ~ HISEI + MotherEdu + Books, data = PISA2009) summary(lm.beta(m1))
Die Residuenplots sowie die Testung auf quadratische Trends zeigen an, dass für den Bildungsabschluss der Mutter auch eine quadratische Beziehung mit der Lesekompetenz besteht:
# Residuenplots residualPlots(m1, pch = 16)
Die Effekte von Sozialstatus und Büchern werden durch das lineare Modell gut wiedergegeben. Für den Bildungsabschluss der Mutter ist ein leicht nicht-linearer Zusammenhang zu erkennen, der quadratische Trend für die Residuen ist signifikant (signifikantes Ergebnis für den Bildungsabschluss der Mutter). Der Effekt ist dadurch charakterisiert, dass der Zuwachs der Lesekompetenz im unteren Bereich des mütterlichen Bildungsabschlusses stärker ist und im oberen Bereich abflacht.
Auch dem Histogramm war eine Schiefe zu entnehmen, welche durch nichtlineare Terme entstehen können (im niederen Bereich liegen mehr Werte; eine Linksschiefe/Rechtssteile ist zu erkennen).
res <- studres(m1) # Studentisierte Residuen als Objekt speichern df_res <- data.frame(res) # als Data.Frame für ggplot # Grafisch: Histogramm mit Normalverteilungskurve library(ggplot2) ggplot(data = df_res, aes(x = res)) + geom_histogram(aes(y =..density..), bins = 15, # Wie viele Balken sollen gezeichnet werden? colour = "blue", # Welche Farbe sollen die Linien der Balken haben? fill = "skyblue") + # Wie sollen die Balken gefüllt sein? stat_function(fun = dnorm, args = list(mean = mean(res), sd = sd(res)), col = "darkblue") + # Füge die Normalverteilungsdiche "dnorm" hinzu und nutze den empirischen Mittelwert und die empirische Standardabweichung "args = list(mean = mean(res), sd = sd(res))", wähle dunkelblau als Linienfarbe labs(title = "Histogramm der Residuen mit Normalverteilungsdichte", x = "Residuen") # Füge eigenen Titel und Achsenbeschriftung hinzu # Test auf Abweichung von der Normalverteilung mit dem Shpiro Test shapiro.test(res)
Die Frage ist nun, woher die Verstöße gegen die Normalverteilungsannahme kommen. Erste Indizen aus den Partialplots wiesen darauf hin, dass möglicherweise ein quadratischer Effekt des Bildungsabschlusses der Mutter besteht.
Aufnahme eines quadratischen Effekts
$***$ zum Teil basierend auf Unterlagen von Prof. Johannes Hartig
Wird für den Bildungsabschluss der Mutter mit der Funktion poly ein linearer und qudratischer Trend in das Regressionsmodell aufgenommen, wird der quadratische Trend signifikant und das Modell erklärt signifikant mehr Varianz als ohne den quadratischen Trend: Um diese Ergebnisse zu sehen müssen wir zunächst ein quadratisches Regressionsmodell schätzen. Wir interessieren uns anschließend für die standardisierten Ergebnisse (summary und lm.beta). Den quadratischen Verlauf erhalten wir, indem wir innerhalb des linearen Modells poly auf den Bildungsabschluss der Mutter anwenden. poly nimmt als zweites Argument die Potenz, für welche wir uns interessieren; hier 2:
# Modell mit linearem + quadratischem Trend für MotherEdu: m1.b <- ...
# Verwenden Sie den Befehl lm(...~...)
# Verwenden Sie den Befehl lm(...~...) # sowie innerhalb von lm poly(..., 2)
# Verwenden Sie den Befehl lm(...~...) # sowie innerhalb von lm poly(..., 2) # wenden Sie dies auf MotherEdu an
# Verwenden Sie den Befehl lm(...~...) # sowie innerhalb von lm poly(..., 2) # wenden Sie dies auf MotherEdu an # wenden Sie anschließend summary(lm.beta(...)) # an
# Finale Lösung m1.b <- lm(Reading ~ HISEI + poly(MotherEdu, 2) + Books, data = PISA2009) summary(lm.beta(m1.b))
# For further use in exercises m1 <- lm(Reading ~ HISEI + MotherEdu + Books, data = PISA2009) m1.b <- lm(Reading ~ HISEI + poly(MotherEdu, 2) + Books, data = PISA2009)
m1.b <- lm(Reading ~ HISEI + poly(MotherEdu, 2) + Books, data = PISA2009)
Mit dem folgenden Befehl können wir auf eine simple Weise das Inkrement bestimmen.
# Vergleich mit Modell ohne quadratischen Trend summary(m1.b)$r.squared - summary(m1)$r.squared # Inkrement
Wir möchten dieses Inkrement auf Signifikanz prüfen. Dies geht mit dem anova Befehl.
# Intro into excercise m1 <- lm(Reading ~ HISEI + MotherEdu + Books, data = PISA2009) m1.b <- lm(Reading ~ HISEI + poly(MotherEdu, 2) + Books, data = PISA2009) # Begin grading: res <- anova(m1, m1.b) sol2 <- anova(m1.b, m1) grade_result( fail_if(~ (!identical(.result, res) & !identical(.result, sol2)), 'Der Befehl anova muss in diesem Zusammenhang zwei Argumente entgegennehmen: die zwei Regressions-Objekte, die wir zuvor abgespeichert hatten. Schauen Sie in der Übung zuvor an, wie dieses Modell hieß, falls Sie sich nicht sicher sind.'), fail_if(~ (identical(.result, sol2)), 'Das ist an sich RICHTIG, allerdings sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell (mit weniger Prädiktoren und Parametern, die zu schätzen sind) zuerst übergeben werden. Hier: m1, da sonst die df negativ sind und auch die Änderung in den Sum of Sq (Quadratsumme) negativ sind! `R` erkennt dies zwar und testet trotzdem die richtige Differenz auf Signifikanz, aber wir wollen uns besser vollständig korrekt aneignen!'), pass_if(~ identical(.result, res)), correct = 'Sehr gut! Dies ist der Modellvergleich des lineares und des quadratischen Modells. Hier sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell (mit weniger Prädiktoren und Parametern, die zu schätzen sind) zuerst übergeben werden. Hier: m1, da sonst die df negativ sind und auch die Änderung in den Sum of Sq (Quadratsumme) negativ sind! `R` erkennt dies zwar und testet trotzdem die richtige Differenz auf Signifikanz, aber wir wollen uns besser vollständig korrekt aneignen!', incorrect = 'Leider falsch.', glue_correct = '{.correct}', glue_incorrect = '{.incorrect} {.message}')
Erzeugt man für das erweiterte Modell Residuenplots, ist der quadratische Trend beim Bildungsabschluss komplett verschwunden - er ist ja schon im Modell enthalten und bildet sich somit nicht mehr in den Residuen ab:
residualPlots(m1.b, pch = 16)
Was bedeutet nun dieser Effekt inhaltlich? Um dies genauer zu verstehen, stellen wir die um die anderen Variablen bereinigte Beziehung zwischen dem Bildungsabschluss der Mutter und der Leseleistung grafisch dar.
X <- scale(poly(PISA2009$MotherEdu, 2)) std_par_ME <- c(0.1588, -0.1436) pred_effect_ME <- X %*% std_par_ME std_ME <- X[,1] data_ME <- data.frame(std_ME, pred_effect_ME) ggplot(data = data_ME, aes(x = std_ME, y = pred_effect_ME)) + geom_point(pch = 16, col = "blue", cex = 4)+ labs(y = "std. Leseleistung | Others", x = "std. Bildungsabschluss der Mutter | Others", title = "Standardisierte bedingte Beziehung zwischen\n Bildungsabschluss der Mutter und Leseleistung")
Für den Grafik-Code sowie weitere Informationen zu quadaratischen Effekten und Funktionen siehe Appendix A. Die Grafik zeigt die vorhergesagte Beziehung zwischen den standardisierten Werten des Bildungsabschlusses der Mutter sowie der Leseleistung. Hierbei steht erneut | für "gegeben" (wie etwa beim Partialplot mit avPlots aus der vergangenen Sitzung: PsyBSc7::Sitzung_7()). Wir sehen also den um die anderen Variablen im Modell bereinigten Effekt zwischen Bildungsabschluss und Leseleistung. Hierbei ist ein starker mittlerer Anstieg der Leseleistung (-1 bis ca. 0.1) für einen Anstieg des Bildungsabschlusses von deutlich unterdurchschnittlich bis durchschnittlich (von -2.5 bis 0) zu sehen. Danach ist die Beziehung zwischen Leseleistung und Bildungsabschluss fast horizontal (Veränderung geringer als 0.1), was dafür spricht, dass es für einen durchschnittlichen bis überdurchschnittlichen Bildungsabschluss der Mutter (von 0 bis 1.5) kaum eine Beziehung zwischen den Variablen gibt. Dies bedeutet, dass besonders im unterdurchschnittlichen Bereich der mütterlichen Bildung Unterschiede zwischen Müttern einen starken Zusammenhang mit der Leseleistung ihrer Kinder zeigen. Wenn das Bildungsniveau der Mutter jedoch durchschnittlich oder überdurchschnittlich ist, scheint der Zusammenhang beinahe zu verschwinden.
Die grobe Gestalt der Beziehung hätten wir auch aus dem Koeffizienten ablesen können. Der Koeffizient des quadratischen Teils war negativ, was für eine invers-u-förmige (konkave) Funktion steht. Das Einzeichnen hilft uns jedoch, das genaue Ausmaß zu verstehen (siehe auch Appendix A). Auch hatten wir gesehen, dass der lineare Teil des Bildungsabschlusses der Mutter keinen statistisch signifikanten Beitrag zur Vorhersage geleistet hat. Jedoch gehört zu einer quadratischen Funktion immer auch ihr linearer Anteil dazu. Aus diesem Grund können wir unsere Stichprobe nur angemessen beschreiben, wenn wir den linearen Trend des Bildungsabschlusses der Mutter im Regressionsmodell beibehalten. Um das genaue Ausmaß besser zu verstehen, manipulieren Sie doch einmal die Beziehung, die wir soeben grafisch gesehen haben im folgenden interaktiven Block. Hierbei können Sie den linearen und den quadratischen Effekt verändern und sich die Auswirkungen auf die Grafik (die Beziehung zwischen Bildungsabschluss der Mutter und Leseleistung) ansehen. Die Default-Einstellungen sind identisch zu der oberen Grafik. curve plottet eine Linie und nimmt x automatisch als Argument der Funktion, somit wird $f(x)=\text{linear}x+\text{quadratisch}x^2$ geplottet. Probieren Sie doch einmal aus, was passiert, wenn Sie den linearen Teil auf 0 setzen oder das Vorzeichen des quadratischen Anteils ändern!
linear <- .1588 quadratisch <- -.1436 curve(linear * x + quadratisch * x^2, xlim = c(-2, 2))
Mit Hilfe von poly(X, p), lassen sich Polynome bis zum Grad $p$ (als $X, X^2,\dots,X^{p-1},X^p$) in die Regression mit aufnehmen, ohne, dass sich die Parameterschätzungen der anderen Potenzen von $X$ ändern. Wenn Sie noch mehr über die Funktion poly und ihre Vorteile erfahren möchten, dann schauen Sie sich doch mal den Appendix A an.
Mit Hilfe von I() lassen sich innerhalb des lm Befehls zu dem noch weitere Funktionale hinzufügen, ohne diese vorher erzeugen zu müssen. Beispielsweise ließe sich durch lm(Y ~ X + I(sin(X))) + I(exp(sqrt(X)) folgendes Regressionsmodell schätzen: $Y = \beta_0+\beta_1X + \beta_2\sin(X) + \beta_3e^{\sqrt{X}} + \varepsilon$. Allerdings lassen sich so nicht Wachstumsraten modellieren (z.B. exponentielles oder logarithmisches Wachstum) - hierzu müssten die Variablen tatsächlich transformiert werden.
Weitere nichtlineare Effekte
Im Folgenden werden zum Teil fortgeschrittene Inhalte verwendet. Bitte sehen Sie dies als Bonus zum Vorlesungsstoff an.
Ein sehr wichtiges Vorhersagemodell zu Zeiten der Corona-Pandemie oder in Anbetracht von starkem weltweitem Bevölkerungswachstum ist das exponentielle Wachstum. Um dieses genauer zu verstehen, werden einige Grundlagen im Umgang mit Exponenten benötigt. Zum Beispiel müssen wir uns überlegen, welche Rate und welche Basis ein exponentieller Verlauf hat. Beispielsweise wächst $10^x$ deutlich schneller als $2^x$ (für gleiches wachsendes $x$). Wie gut kennen Sie sich noch mit Exponenten und Logarithmen aus?
question("Welche Aussagen über Logarithmen bzw. Exponenten stimmen?", answer("Die Gleichung $10^x = 100$ können wir bspw. mit dem 10er-Logarithmus $log_{10}(100)$ berechnen. Die Antwort ist 2.", correct = T), answer("Die Gleichung $10^x = 100$ können wir bspw. mit dem 10er-Logarithmus $log_{10}(100)$ berechnen. Die Antwort ist 3."), answer("Logarithmen können NICHT in einander umgerechnet werden."), answer("Es gilt: $log(a+b)=log(a)log(b)$."), answer("Es gilt: $log(ab)=log(a)+log(b)$.", correct=T), answer("$log_a(x)$ ist die Umkehrfunktion zu $a^x$. Dies bedeutet, dass $log_a(a^x)=x.$", correct = T), answer("Es gilt: $a^{c+d}=a^c+a^d$."), answer("Es gilt: $a^{c+d}=a^ca^d$.", correct = T), answer("Es gilt immer $a^x < b^x$ für a < b und alle x."), answer("Es gilt immer $a^x < b^x$ für 1 < a < b und alle x.", correct = T), answer("Es gilt immer $a^x > b^x$ für 0 < a < b < 1 und alle x.", correct = T), answer("Es gilt immer $a^x < b^x$ für 0 < a < b < 1 und alle x.") )
Im folgenden Exkurs haben Sie die Möglichkeit Ihr Wissen nochmal aufzufrischen und zu erweitern. Ihnen werden Quellen genannt, in welchen Sie weitere Inhalte nachlesen können.
Exkurs: Exponenten und Logarithmen
Um eine exponentielle Gleichung nach $x$ aufzulösen wird der Logarithmus verwendet, denn er ist die Umkehrfunktion zum Exponentiellen (siehe beispielsweise hier für eine Wiederholung zum Logarithmus und Logarithmusrechenregeln; Pierce (2018), Exponenten; Pierce (2019) oder Exponentenrechenregeln; Pierce (2018)). In der Schule wird zunächst der Logarithmus zur Basis 10 gelehrt. So kann man die Gleichung $10^x = 100$ leicht auflösen: Dazu müssen beide Seiten (zur Basis 10) logarithmiert werden: $\log_{10}(10^x) = \log_{10}(100)$, denn es gilt im Allgemeinen:
$$\log_a(a)=1$$
$$\log_a(bc)=\log_a(b)+\log_a(c)$$ $$\log_a(b^c)=c\log_a(b)$$
$$x=\log_a(a^x)=a^{\log_a(x)}$$ (hier ist $\log_a$ der Logarithmus zur Basis $a$ und $a, b$ und $c$ sind Zahlen, wobei der Einfachheit halber $a,b,c>0$ gelten soll)
Also rechnen wir: $\log_{10}(10^x) = x\log_{10}(10)=x*1=x$ auf der linken Seite und $\log_{10}(100)=2$ auf der rechten Seite. So erhalten wir $x=2$ als Lösung. Natürlich wussten wir schon vorher, dass $10^2=100$ ergibt; es ist aber eine schöne Veranschaulichung, wofür Logarithmen unter anderem gebraucht werden.
Es muss allerdings nicht immer der richtige Logarithmus zur jeweiligen Basis sein. Eine wichtige Logarithmusrechenregel besagt, dass der Logarithmus zur Basis $a$ leicht mit Hilfe des Logarithmus zur Basis $b$ berechnet werden kann:
$$\log_{a}(x) = \frac{\log_{b}(x)}{\log_{b}(a)},$$ siehe bspw. Logarithmus und Logarithmusrechenregeln; Pierce (2018). Da diese Gesetzmäßigkeit gilt, wird in der Mathematik und auch in der Statistik meist nur der natürliche Logarithmus (ln; logarithmus naturalis) verwendet. Die Basis ist hierbei $e$ ($\approx 2.718\dots$). Mit den obigen Rechenregeln sowie der Rechenregeln für Exponenten können wir jedes exponentielles Wachstum als Wachstum zur Basis e darstellen, denn
$$a^x = e^{\text{ln}(a^x)}=e^{x\text{ln}(a)},\quad\text{ z.B. ist }\quad 10^3=e^{\text{ln}(10^3)}=e^{3\text{ln}(10)}=1000$$
Da sich durch e und ln alle exponentiellen Verläufe darstellen lassen, wird in der Mathematik häufig log als das Symbol für den natürlich Logarithmus verwendet; so ist es auch in R: ohne weitere Einstellungen ist log der natürliche Logarithmus und exp ist die Exponentialfunktion. Mit Hilfe von log(..., base = 10) erhalten sie beispielsweise den 10er-Logarithmus. Probieren Sie die obige Gleichung selbst aus:
# gleiches Ergebnis: 10^3 exp(3*log(10)) # gleiches Ergebnis: log(10^3, base = 10) # Logarithmus von 1000 zur Basis 10 log(10^3)/log(10) # mit ln # gleiches Ergebnis: log(9, base = 3) # Logarithmus von 9 zur Basis 3 log(9)/log(3) # mit ln
Das Modellieren von exponentiellem Wachstum
Wir betrachten nun eine allgemeine exponentielle Wachstumsfunktion $$f(x) = a\ b^{c\ x}$$ (hierbei ist $a$ ein Vorfaktor, der die Ausprägung an der Stelle $x=0$ beschreibt, $b$ ist die Basis des exponentiellen Wachstums und $c$ ist ein eigentlich redundanter Ratenparameter). Wegen der Beliebigkeit der Basis ist dies gleich $$f(x) = e^{\text{ln}(a\ b^{c\ x})}=e^{\text{ln}(a) + \text{ln}(b)cx}.$$ Nun sind $\text{ln}(a)$ und $\text{ln}(b)c$ zwei Konstanten, die wir einfach umbenennen dürfen. Wir können sie beispielsweise $\beta_0$ und $\beta_1$ nennen; also $\beta_0 := \text{ln}(a)$ und $\beta_1:=\text{ln}(b)c$. Folglich steht $$f(x) = e^{\beta_0 + \beta_1x}$$ für beliebiges exponentielle Wachstum (wir erhalten den Verlauf "$2^x$" indem wir $\beta_0 = \text{ln}(a) = 0$ und $\beta_1=\text{ln}(b)c = \text{ln}(2)$ wählen; also $a=1, b=2$ und $c=1$). Gleichzeit bedeutet dies, dass $$\text{ln}(f(x)) = \beta_0 + \beta_1x$$ eine lineare Funktion ist, welche wir ganz einfach mit einer Regressionanalyse, also lm in R, untersuchen können! Wir können somit sagen, dass wir durch Logarithmieren von $f(x)$ in der Lage sind, das exponentielle Wachstum zu Linearisieren, also das exponentielle Wachstum in eine lineare Funktion zu transformieren, welche mit Auswertungsinstrumenten für lineare Funktionen untersucht werden können. Dies bedeutet, dass wir nach Logarithmieren der abhängigen Variable in der Lage sind die gesamte Klasse der exponentiellen Funktionen/des exponentiellen Wachstums für die Modellierung unserer Daten zu verwenden. Die folgenden zwei Grafiken verdeutlichen dies (Sie können hier die Parameter beliebig verändern und sich den Effekt auf die Grafiken ansehen; die horizontale gestrichelte Linien repräsentiert jeweils den Wert $a$ bzw. $\text{ln}(a)$, an welchem die Kurve $f$ bzw. $\ln(f)$ die y-Achse schneiden):
################## #### Einstellen der Koeffizienten und berechnen von f(x) #### x <- seq(-1,2,0.1) # x = Variablen (als Zahlen zwischen -1 und 2) a <- 2 # Vorfaktor, der die Ausprägung an der Stelle x=0 beschreibt b <- 3 # Basis des exponentiellen Wachstums c <- 1.5 # *eigentlich redundanter* Ratenparameter f <- a*b^(c*x) # f(x), eine exponentiell-wachsende Funktion in x ################## #### Plot von f(x) vs. x #### plot(x = x, y = f, type = "l", col = "blue", lwd = 2, main = "Plot von f(x) vs. x") # plotte f(x) gegen x abline(v = 0, lwd = 0.7) # y-Achse, v = 0 zeichnet eine vertikale Linie bei x = 0 abline(h = a, lty = 3) # im Punkt a schneidet f (das exponentielle Wachstum) die y-Achse (x=0), h = a zeichnet zu y = a eine horizontale Linie ################## #### Plot von ln(f(x)) vs. x #### plot(x = x, y = log(f), type = "l", col = "blue", lwd = 2, main = "Plot von ln(f(x)) vs. x") # plotte ln(f(x)) gegen x abline(v = 0, lwd = 0.7) # y-Achse, v = 0 zeichnet eine vertikale Linie bei x = 0 abline(h = log(a), lty = 3) # im Punkt log(a) schneidet log(f) (das linearisierte exponentielle Wachstum) die y-Achse (x=0), h =llog(a) zeichnet zu y = log(a) eine horizontale Linie
Modellierung von exponentiellem Wachstum
Das Modellieren von exponentiellem Wachstum am Beispiel der Weltbevölkerung von 1800 bis 2020
Nun können wir beispielsweise die Entwicklung der Weltbevölkerung von 1800 bis 2020 modellieren. Dazu müssen wir zunächst die Daten laden:
Die Dokumentation des Datensatzes mit Datenquellen sind hier einzusehen: gapminder.org-Dokumentationen. Im Paket PsyBSc7 ist unter dem Namen WorldPopulation ein vorverarbeiter Datensatz enthalten. Diesen können wir ganz einfach wie folgt laden und ansehen:
data("WorldPopulation", package = "PsyBSc7") # Bereits zusammengefasster Datensatz von 1800-2020 head(WorldPopulation) # Datenstruktur
In der ersten Spalte steht das Jahr; in der 2. die Weltbevölkerungsgröße. Wir wollen uns dies grafisch ansehen. Verwenden Sie ggplot, um die Population (Population) gegen das Jahr (Year) abzutragen und verwenden Sie geom_point().
ggplot(...)
# Sie müssen den Datensatz sowie die Namen der Variablen verwenden ggplot(data = ..., aes(x = ..., y = ...))
# Sie müssen den Datensatz sowie die Namen der Variablen verwenden ggplot(data = ..., aes(x = ..., y = ...)) + geom_point()
ggplot(data = WorldPopulation, aes(x = Year, y = Population))+geom_point()
Das Diagramm lässt deutlich einen nichtlinearen Anstieg der Weltbevölkerung von 1800 bis 2020 vermuten. Auffällig ist auch der leichte Knick, der um 1950 zu vermuten ist und ab welchem die Bevölkerung, deskriptiv gesehen, noch stärker wächst. Dieser Knick ist zum Teil durch das Ende des Krieges, aber auch durch modernere Landwirtschaft und das Aufkommen von neuen Medikamenten (z.B. Penicilline) zu erklären.
Lineares Modell für das Bevölkerungswachstum
Wir wollen uns naiverweise ein lineares Regressionmodell, also einen linearen Verlauf, der Weltbevölkerung vorhergesagt durch das Jahr, ansehen.
ggplot(data = WorldPopulation, aes(x = Year, y = Population))+ geom_point()+geom_smooth(method="lm", formula = "y~x") # plotte linearen Verlauf m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf summary(m_l) ######### ### Normalverteilung der Residuen? ## res <- studres(m_l) # Studentisierte Residuen als Objekt speichern df_res <- data.frame(res) # als Data.Frame für ggplot # Grafisch: Histogramm mit Normalverteilungskurve ggplot(data = df_res, aes(x = res)) + geom_histogram(aes(y =..density..), bins = 10, # Wie viele Balken sollen gezeichnet werden? colour = "blue", # Welche Farbe sollen die Linien der Balken haben? fill = "skyblue") + # Wie sollen die Balken gefüllt sein? stat_function(fun = dnorm, args = list(mean = mean(res), sd = sd(res)), col = "darkblue") + # Füge die Normalverteilungsdiche "dnorm" hinzu und nutze den empirischen Mittelwert und die empirische Standardabweichung "args = list(mean = mean(res), sd = sd(res))", wähle dunkelblau als Linienfarbe labs(title = "Histogramm der Residuen mit Normalverteilungsdichte\n für das lineare Modell", x = "Residuen") # Füge eigenen Titel und Achsenbeschriftung hinzu
Durch + geom_smooth(method="lm", formula = "y~x"), kann mit ggplot ein linearer Trend inklusive Konfidenzintervall hinzugefügt werden. Obwohl ein linearer Verlauf sehr unwahrscheinlich erscheint, können mit dem linearen Modell bereits r round(summary(m_l)$r.squared*100, 2)% der Variation der Bevölkerungsdichte durch die Jahreszahl eklärt werden. Wie der Grafik deutlich zu entnehmen ist, sind die Residuen in dieser Regressionanalyse stark abhängig von der Jahreszahl (negatives Residuum von a. 1860-1970 und positive Residuen sonst; Wiederholung: $\varepsilon_i=Y_i-\hat{Y}_i$, wobei $\hat{Y}_i$ der vorhergesagte Wert ist). Auch wenn wir uns das zugehörige Histogramm der Residuen ansehen, widerspricht dieses der Annahme auf Normalverteilung.
Quadratisches Modell für das Bevölkerungswachstum
Vielleicht ist auch hier ein quadratischer Effekt versteckt? Wir möchten dem auf den Grund gehen und nutzen wieder die Funktion poly, um ein Polynom 2. Grades (eine quadratische Funktion) der Jahreszahl in unsere Analysen mit aufzunehmen.
m_q <- ... # quadratischer Verlauf summary(...)
# Finale Lösung m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf summary(m_q)
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf
Das Varianzinkrement beläuft sich auf:
summary(m_q)$r.squared - summary(m_l)$r.squared # Inkrement
Testen Sie dieses auf Signifikanz mit anova.
# Intro into excercise m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf # Begin grading: res <- anova(m_l, m_q) sol2 <- anova(m_q, m_l) grade_result( fail_if(~ (!identical(.result, res) & !identical(.result, sol2)), 'Der Befehl anova muss in diesem Zusammenhang zwei Argumente entgegennehmen: die zwei Regressions-Objekte, die wir zuvor abgespeichert hatten. Schauen Sie in der Übung zuvor an, wie dieses Modell hieß, falls Sie sich nicht sicher sind.'), fail_if(~ (identical(.result, sol2)), 'Das ist an sich RICHTIG, allerdings sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell (mit weniger Prädiktoren und Parametern, die zu schätzen sind) zuerst übergeben werden. Hier: m_l, da sonst die df negativ sind und auch die Änderung in den Sum of Sq (Quadratsumme) negativ sind! `R` erkennt dies zwar und testet trotzdem die richtige Differenz auf Signifikanz, aber wir wollen uns besser vollständig korrekt aneignen!'), pass_if(~ identical(.result, res)), correct = 'Sehr gut! Dies ist der Modellvergleich des linearen und des quadratischen Modells. Hier sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell (mit weniger Prädiktoren und Parametern, die zu schätzen sind) zuerst übergeben werden. Hier: m_l, da sonst die df negativ sind und auch die Änderung in den Sum of Sq (Quadratsumme) negativ sind! `R` erkennt dies zwar und testet trotzdem die richtige Differenz auf Signifikanz, aber wir wollen uns besser vollständig korrekt aneignen!', incorrect = 'Leider falsch.', glue_correct = '{.correct}', glue_incorrect = '{.incorrect} {.message}')
Auch dem geom_smooth Befehl kann einfach das modifizierte Regressionsmodell übergeben werden. Stellen Sie die Weltbevölkerung sowie das lineare und das quadratische Modell dar und färben Sie die Kurve in dunkelblau (col = "darkblue").
ggplot(data = WorldPopulation, aes(x = Year, y = Population))+ geom_point()+geom_smooth(method="lm", formula = "y~x")+ # plotte linearen Verlauf (bis hier her ist es Code der Aufgabe zuvor!) ... # plotte quadratischen Verlauf
# Finale Lösung ggplot(data = WorldPopulation, aes(x = Year, y = Population))+ geom_point()+geom_smooth(method="lm", formula = "y~x")+ # plotte linearen Verlauf geom_smooth(method="lm", formula = "y~poly(x,2)", col = "darkblue") # plotte quadratischen Verlauf
Durch den quadratischen Verlauf lassen sich r round(summary(m_q)$r.squared*100, 2)% der Variation der Bevölkerungsdichte erklären, was einem signifikantem Varianzinkrement von r round(summary(m_q)$r.squared*100 - summary(m_l)$r.squared*100, 2)% entspricht (mit einer Irrtumswahrscheinlichkeit von 5% ist das Inkrement in der Population nicht null. Dies ist äquivialent zu folgdender Aussage: mit einer Irrtumswahrscheinlichkeit von 5% ist der Effektparameter (der Regressionskoeffizient) des quadratischen Verlaufs in der Population nicht null; dies spricht folglich für einen quadratischen im Gegensatz zu einem linearen Verlauf). Der Grafik ist deutlich zu entnehmen, dass der quadratische Verlauf nicht weit vom empirischen entfernt liegt.
Exponentielles Modell für das Bevölkerungswachstum
Nun wollen wir prüfen, ob ein exponentieller Verlauf die Daten noch besser beschreibt. Dazu müssen wir zunächst die Weltpopulation logarithmieren und können anschließend ein lineares Modell verwenden. Transformieren Sie hierzu die Variablen und speichern diese als weitere Spalte in unserem Datensatz ab. Ordnen Sie diese folgendem Namen zu: WorldPopulation$log_Population <- ... und geben Sie anschließend WorldPopulation$log_Population aus. Der Name ist entscheidend, da wir diesen später weiter verwenden wollen.
# Intro into excercise # Begin grading: res <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung grade_result( fail_if(~ (!identical(.result, res)), 'Sie müssen den Befehl `log` verwenden.'), pass_if(~ identical(.result, res)), correct = 'Sehr gut! Sie haben die logarithmierte Weltbevölkerung unserem Datensatz hinzugefügt.', incorrect = 'Leider falsch.', glue_correct = '{.correct}', glue_incorrect = '{.incorrect} {.message}')
WorldPopulation$log_Population <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf WorldPopulation$log_Population <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung
Wir wollen anschließend ein einfaches Regressionsmodell schätzen, in welchem die logarithmierte Bevölkerungszahl die abhängige Variable darstellt. Speichern Sie das lm Objekt unter dem Namen m_log ab. Geben Sie die summary für m_log aus.
m_log <- lm(... ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
m_log <- lm(log_Population ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala) summary(m_log)
m_log <- lm(log_Population ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf WorldPopulation$log_Population <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung m_log <- lm(log_Population ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
Das lineare Modell für die logarithmierte Bevölkerunsdichte scheint gut zu den Daten zu passen. Insgesamt können r round(summary(m_log)$r.squared*100, 2) % der Variation in den Daten durch den Zeitverlauf erklärt werden; etwas weniger als durch das quadratische Wachstum, durch welches r round(summary(m_q)$r.squared*100, 2) % der Variation in den Daten durch den Zeitverlauf erklärt werden konnten. Die zugehörige Grafik sieht so aus:
ggplot(data = WorldPopulation, aes(x = Year, y = log_Population))+ geom_point()+geom_smooth(method="lm", formula = "y~x", col = "red")+ labs(title = "Logarithmierte Weltbevölkerung vs. Jahr")
Wieso erklärt nun das quadratische Modell mehr Variation? Wir wollen die Vorhersage auch retransformiert plotten, um sie mit dem quadratischen Verlauf grafisch vergleichen zu können. Dazu müssen wir die vorhergesagten Werte unseres logarithmischen Modells verwenden und mit mit der Umkehrfunktion des Logarithmus retransformieren. Diese können wir einem lm-Objekt mit predict entlocken.
WorldPopulation$pred_Pop_exp <- ... # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung) WorldPopulation
# Verwenden Sie predict(...)
# Verwenden Sie predict(...) # und transformieren Sie mit exp(...)# welche die Umkehrfunktion des Log ist
WorldPopulation$pred_Pop_exp <- exp(predict(m_log)) # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung) WorldPopulation
WorldPopulation$pred_Pop_exp <- exp(predict(m_log)) # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung)
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf WorldPopulation$log_Population <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung m_log <- lm(log_Population ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala) WorldPopulation$pred_Pop_exp <- exp(predict(m_log)) # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung)
Dem Datensatz haben wir nun eine neue Spalte hinzugefügt, welche die vorhergesagten Populationswerte enthält (retransformiert; nicht mehr in Log-Skala).
ggplot(data = WorldPopulation, aes(x = Year, y = Population))+ geom_point()+geom_smooth(method="lm", formula = "y~x")+ # plotte linearen Verlauf geom_smooth(method="lm", formula = "y~poly(x,2)", col = "darkblue")+ # plotte quadratischen Verlauf geom_line(aes(x = Year, y = ...), col = "red", lwd = 1.5)
ggplot(data = WorldPopulation, aes(x = Year, y = Population))+ geom_point()+geom_smooth(method="lm", formula = "y~x")+ # plotte linearen Verlauf geom_smooth(method="lm", formula = "y~poly(x,2)", col = "darkblue")+ # plotte quadratischen Verlauf geom_line(aes(x = Year, y = pred_Pop_exp), col = "red", lwd = 1.5)
Das Diagramm der retransformierten vorhergesagten Werten signalisiert, dass ein exponentielles Wachstumsmodell die Daten gut beschreibt, allerdings scheint der quadratische Trend vor allem ab ca. 1975 die Daten besser zu beschreiben. Wir schauen uns dazu die Residuen des logarithmischen Modells an.
residualPlot(..., col = "red") # Residualplot
# Finale Lösung residualPlot(m_log, col = "red") # Residualplot
Dem Residualplot ist zu entnehmen, dass ein nicht-linearer Verlauf angemessen wäre. Wir gehen dieser Vermutung nach, indem wir dem Log-Plot einen quadratischen Verlauf hinzufügen.
ggplot(data = WorldPopulation, aes(x = Year, y = log_Population))+ geom_point()+geom_smooth(method="lm", formula = "y~x", col = "red")+ geom_smooth(method="lm", formula = "y~poly(x,2)", col = "gold3")+ labs(title = "Logarithmierte Weltbevölkerung vs. Jahr")
Vielleicht wächst die Bevölkerung sogar schneller als exponentiell? Die gelbe Linie im logarithmierten Plot lässt dies vermuten. Um dies genauer zu untersuchen fügen wir auch in das logarithmierte Modell einen quadratischen Trend der Zeit ein.
Quadratisch-exponentielles Modell für das Bevölkerungswachstum
Dies wollen wir schätzen, indem wir im Regressionsmodell der logarithmierten Populationsgröße einen quadratischen Trend (poly(...,2)) für die Zeit definieren. Geben Sie für dieses Regressionsmodell die summary aus.
m_log_quad <- summary(...)
m_log_quad <- lm(log_Population ~ poly(Year, 2), data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala) summary(m_log_quad)
m_log_quad <- lm(log_Population ~ poly(Year, 2), data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf WorldPopulation$log_Population <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung m_log <- lm(log_Population ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala) WorldPopulation$pred_Pop_exp <- exp(predict(m_log)) # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung) m_log_quad <- lm(log_Population ~ poly(Year, 2), data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
Testen Sie dieses Modell im Vergleich zu m_log auf Signifikanz mit anova .
# Intro into excercise WorldPopulation$log_Population <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung m_log <- lm(log_Population ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala) m_log_quad <- lm(log_Population ~ poly(Year, 2), data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala) # Begin grading: res <- anova(m_log, m_log_quad) sol2 <- anova(m_log_quad, m_log) grade_result( fail_if(~ (!identical(.result, res) & !identical(.result, sol2)), 'Der Befehl anova muss in diesem Zusammenhang zwei Argumente entgegennehmen: die zwei Regressions-Objekte, die wir zuvor abgespeichert hatten. Schauen Sie in der Übung zuvor an, wie dieses Modell hieß, falls Sie sich nicht sicher sind.'), fail_if(~ (identical(.result, sol2)), 'Das ist an sich RICHTIG, allerdings sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell (mit weniger Prädiktoren und Parametern, die zu schätzen sind) zuerst übergeben werden. Hier: m_log, da sonst die df negativ sind und auch die Änderung in den Sum of Sq (Quadratsumme) negativ sind! `R` erkennt dies zwar und testet trotzdem die richtige Differenz auf Signifikanz, aber wir wollen uns besser vollständig korrekt aneignen!'), pass_if(~ identical(.result, res)), correct = 'Sehr gut! Dies ist der Modellvergleich des linearen und des quadratischen logarithmierten Modells. Hier sollte dem anova-Befehl immer das "kleinere" (restriktivere) Modell (mit weniger Prädiktoren und Parametern, die zu schätzen sind) zuerst übergeben werden. Hier: m_log, da sonst die df negativ sind und auch die Änderung in den Sum of Sq (Quadratsumme) negativ sind! `R` erkennt dies zwar und testet trotzdem die richtige Differenz auf Signifikanz, aber wir wollen uns besser vollständig korrekt aneignen!', incorrect = 'Leider falsch.', glue_correct = '{.correct}', glue_incorrect = '{.incorrect} {.message}')
Durch den quadratisch-exponentiellen Verlauf (bzw. den quadratischen Verlauf in den logarithmierten Daten) lassen sich r round(summary(m_log_quad)$r.squared*100, 2)% der Variation der Bevölkerungsdichte erklären, was einem signifikantem Varianzinkrement von r round(summary(m_log_quad)$r.squared*100 - summary(m_log)$r.squared*100, 2)% im Vergleich zum reinen exponentiellen Verlaufsmodell entspricht. Interessant zu sehen ist, dass fast 100% der Variation im Datensatz erklärbar ist. Eine Übersicht über $R^2$ in den vier Modellen ist in Appendix B einzusehen.
Stellen wir dies grafisch dar. Dazu müssen wir wieder die Vorhersage abspeichern. Diese können wir einem lm-Objekt mit predict entlocken.
WorldPopulation$pred_Pop_exp_quad <- ... # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung) WorldPopulation
# Verwenden Sie predict(...)
# Verwenden Sie predict(...) # und transformieren Sie mit exp(...)# welche die Umkehrfunktion des Log ist
WorldPopulation$pred_Pop_exp_quad <- exp(predict(m_log_quad)) # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung) WorldPopulation
WorldPopulation$pred_Pop_exp_quad <- exp(predict(m_log_quad))
m_l <- lm(Population ~ Year, data = WorldPopulation) # linearer Verlauf m_q <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf WorldPopulation$log_Population <- log(WorldPopulation$Population) # Logarithmus der Weltbevölkerung m_log <- lm(log_Population ~ Year, data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala) WorldPopulation$pred_Pop_exp <- exp(predict(m_log)) # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung) m_log_quad <- lm(log_Population ~ poly(Year, 2), data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala) WorldPopulation$pred_Pop_exp_quad <- exp(predict(m_log_quad))
Dem Datensatz haben wir nun eine neue Spalte hinzugefügt, welche die vorhergesagten Populationswerte enthält, die durch das quadratisch-exponentielle Modell vorhergesagt werden (retransformiert; nicht mehr in Log-Skala). Im Folgenden ist die finale Grafik mit dem linearen, quadratischen, exponentiellen und dem quadratisch-exponentiellen Verlauf dargestellt. Außerdem schauen wir uns das Histogramm sowie den Residuenplot des quadratisch-exponentiellen Modells an (in logarithmierter Skala, in welcher das Modell auch geschätzt wurde!).
ggplot(data = WorldPopulation, aes(x = Year, y = Population))+ geom_point()+geom_smooth(method="lm", formula = "y~x")+ # plotte linearen Verlauf geom_smooth(method="lm", formula = "y~poly(x,2)", col = "darkblue")+ # plotte quadratischen Verlauf geom_line(aes(x = Year, y = pred_Pop_exp), col = "red", lwd = 1.5)+ geom_line(aes(x = Year, y = pred_Pop_exp_quad), col = "gold3", lwd = 2) # nur quadratisch-exponentiell ggplot(data = WorldPopulation, aes(x = Year, y = Population))+ geom_point()+ geom_line(aes(x = Year, y = pred_Pop_exp_quad), col = "gold3", lwd = 2)+ labs(title = "Beobachtetes und durch das quadratisch-exponentielle Modell\n vorhergesagtes Bevölkerungswachstum") res <- studres(m_log_quad) # Studentisierte Residuen als Objekt speichern df_res <- data.frame(res) # als Data.Frame für ggplot # Grafisch: Histogramm mit Normalverteilungskurve library(ggplot2) ggplot(data = df_res, aes(x = res)) + geom_histogram(aes(y =..density..), bins = 10, # Wie viele Balken sollen gezeichnet werden? colour = "blue", # Welche Farbe sollen die Linien der Balken haben? fill = "skyblue") + # Wie sollen die Balken gefüllt sein? stat_function(fun = dnorm, args = list(mean = mean(res), sd = sd(res)), col = "darkblue") + # Füge die Normalverteilungsdiche "dnorm" hinzu und nutze den empirischen Mittelwert und die empirische Standardabweichung "args = list(mean = mean(res), sd = sd(res))", wähle dunkelblau als Linienfarbe labs(title = "Histogramm der Residuen mit Normalverteilungsdichte\n für das logarithmierte quadratische Modell", x = "Residuen") # Füge eigenen Titel und Achsenbeschriftung hinzu residualPlot(m_log_quad)
Es ist deutlich zu sehen, dass der Knick um 1950 durch das quadratisch-exponentielle Wachstum am besten darzustellen/ zu modellieren ist. Das Histogramm der Residuen des logarithmierten quadratischen Modells (welches das quadratisch-exponentielle Wachstum der Bevölkerung darstellt) ist etwas rechts-schief/links-steil. Die Residuenplots zeigen außerdem, dass auch hier die Residuen nicht vollständig unsystematisch sind. Dennoch ist die Resiudalvarianz sehr klein.
Da die Modelle des quadratischen und des quadratisch-exponentiellen Verlaufs nicht geschachtelt sind (sie gehen nicht auseinander hervor und sind nicht durch Restringieren, bzw. Nullsetzen, von Regressionsparametern aus einander überführbar), können sie nicht inferenzstatistisch miteinander verglichen werden. Da das quadratisch-exponentielle Modell mehr Variation aufklärt als das rein quadratische Modell und auch grafisch besser zum Verlauf der Daten passt (besonders von 1800 bis ca. 1850 sagt das quadratische Modell ein Verkleinern der Bevölkerung vorher), entscheiden wir uns final für dieses.
Parameterinterpretation des quadratisch-exponentielles Modell für das Bevölkerungswachstum
Was bedeuten nun die Parameter in unserem quadratisch-exponentiellen Modell? Der Regressionskoeffizient des linearen Trends von poly liegt bei r round(coef(m_log_quad)[2], 2) und der Koeffizient des quadratischen Trends bei r round(coef(m_log_quad)[3], 2). Dies spricht für ein beschleunigtes exponentielles Wachstum. Nach diesem Vorhersagemodel scheinen die Menschen sich schneller als exponentiell zu vermehren (zumindest von 1800 bis 2020). Eine Übersicht über die Modelle sehen Sie in Appendix C.
Wären beide Koeffizienten negativ, so würde dies für beschleunigten exponentiellen Zerfall/Abnahme sprechen. Ist der Koeffizient des quadratische Trends negativ, wird von exponentiellem Wachstum mit Dämpfung gesprochen.
Final ist zu sagen, dass die Menschheit am wahrscheinlichsten exponentiell wächst. Ein quadratisches Modell ist im Allgemeinen nicht sinnvoll. Da wir uns in dieser Übung nur einen kleine Zeitausschnitt angesehen haben, kam es zu diesen Abweichungen. Allerdings würde das geschätzte quadratische Modell für das Jahr 0 eine Bevölkerungszahl von r round(predict(object = m_q, newdata = data.frame(Year=c(0)))/10^11, 2)$*10^{11}$ vorhersagen, was absolut nicht sinnvoll ist.
Aufgaben
Die Aufgaben zu dieser Sitzung finden Sie im Appendix D.
Appendix A {#AppendixA}
Exkurs: Was genau macht poly?
X <- 1:10 # Variable X X2 <- X^2 # Variable X hoch 2 X_poly <- poly(X, 2) # erzeuge Variable X und X hoch mit Hilfe der poly Funktion colnames(X_poly) <- c("poly(X, 2)1", "poly(X, 2)2") cbind(X, X2, X_poly)
Die Funktion poly erzeugt sogenannte orthogonale Polynome. Das bedeutet, dass zwar die $x$ und $x^2$ berechnet werden, diese Terme anschließend allerdings so transformiert werden, dass sie jeweils einen Mittelwert von 0 und die gleiche Varianz haben und zusätzlich noch unkorreliert sind:
round(apply(X = cbind(X, X2, X_poly), MARGIN = 2, FUN = mean), 2) # Mittelwerte über die Spalten hinweg berechnen round(apply(X = cbind(X, X2, X_poly), MARGIN = 2, FUN = sd), 2) # Standardabweichung über die Spalten hinweg berechnen round(cor(cbind(X, X2, X_poly)),2) # Korrelationen berechnen
Die Funktion apply führt an der Matrix, welche dem Argument X übergeben wird, entweder über die Zeilen MARGIN = 1 oder über die Spalten MARGIN = 2 (hier jeweils gewählt) die Funktion aus, welche im Argument FUN angegeben wird. So wird zunächst mit FUN = mean der Mittelwert und anschließend mit FUN = sd die Standardabweichung von $X, X^2$ sowie poly(X, 2) berechnet. Der Korrelationsmatrix ist zu entnehmen, dass $X$ und $X^2$ in diesem Beispiel sehr hoch miteinander korrelieren und somit gleiche lineare Informationen enthalten ($\hat{r}_{X,X^2}$ = cor(X, X2) = 0.97), während die linearen und die quadratischen Anteile in poly(X, 2) keinerlei lineare Gemeinsamkeiten haben - sie sind unkorreliert (cor(poly(X,2)1 , poly(X,2)2) = 0). Ein weiterer Vorteil ist deshalb, dass bei sukzessiver Aufnahme der Anteile von poly(X, 2) in ein Regressionmodell, sich die Parameterschätzungen des linearen Terms im Modell nicht (bzw. sehr wenig) ändern. Der Anteil erklärter Varianz bleibt jedoch in allen gleich - die Modelle sind äquivalent, egal auf welche Art und Weise quadratische Terme gebildet werden.
m1.b1 <- lm(Reading ~ HISEI + poly(MotherEdu, 1) + Books, data = PISA2009) summary(lm.beta(m1.b1)) m1.b2 <- lm(Reading ~ HISEI + poly(MotherEdu, 2) + Books, data = PISA2009) summary(lm.beta(m1.b2)) PISA2009$MotherEdu2 <- PISA2009$MotherEdu^2 # füge dem Datensatz den quadrierten Bildungsabschluss der Mutter hinzu m1.c1 <- lm(Reading ~ HISEI + MotherEdu + Books, data = PISA2009) summary(lm.beta(m1.c1)) m1.c2 <- lm(Reading ~ HISEI + MotherEdu + MotherEdu2 + Books, data = PISA2009) summary(lm.beta(m1.c2)) rbind(coef(m1.b1), coef(m1.c1)) rbind(coef(m1.b2),coef(m1.c2))
Einordnung quadratischer Verläufe
Wie kommen wir nun auf die Interpretation der quadratischen Beziehung?
Eine allgemeine quadratische Funktion $f$ hat folgende Gestalt $$f(x):=ax^2 + bx + c,$$ wobei $a\neq 0$, da es sich sonst nicht um eine quadratische Funktion handelt. Wäre $a=0$, würde es sich um eine lineare Funktion mit Achsenabschnitt $c$ und Steigung (Slope) $b$ handeln. Wäre zusätzlich $b=0$, so handelt es sich um eine horizontale Linie bei $y=f(x)=c$. Für betraglich große $x$ fällt $x^2$ besonders ins Gewicht. Damit entscheidet das Vorzeichen von $a$, ob es sich um eine u-formige (falls $a>0$) oder eine umgekehrt-u-förmige (falls $a<0$) Beziehung handelt. Die betragliche Größe von $a$ entscheidet hierbei, wie gestaucht die u-förmige Beziehung (die Parabel) ist. Die reine quadratische Beziehung $f(x)=x^2$ sieht so aus:
a <- 1; b <- 0; c <- 0 x <- seq(-2,2,0.01) f <- a*x^2 + b*x + c data_X <- data.frame(x, f) ggplot(data = data_X, aes(x = x, y = f)) + geom_line(col = "black")+ labs(x = "x", y = "f(x)", title = expression("f(x)="~x^2))
Wir werden diese Funktion immer als Referenz mit in die Grafiken einzeichnen.
a <- 0.5; b <- 0; c <- 0 x <- seq(-2,2,0.01) f <- a*x^2 + b*x + c data_X <- data.frame(x, f) ggplot(data = data_X, aes(x = x, y = f)) + geom_line(col = "blue", lwd = 2)+ geom_line(mapping = aes(x, x^2), lty = 3)+ labs(x = "x", y = "f(x)", title = expression("f(x)="~0.5*x^2))
a <- 2; b <- 0; c <- 0 x <- seq(-2,2,0.01) f <- a*x^2 + b*x + c data_X <- data.frame(x, f) ggplot(data = data_X, aes(x = x, y = f)) + geom_line(col = "blue", lwd = 2)+ geom_line(mapping = aes(x, x^2), lty = 3)+ labs(x = "x", y = "f(x)", title = expression("f(x)="~2*x^2))
a <- -1; b <- 0; c <- 0 x <- seq(-2,2,0.01) f <- a*x^2 + b*x + c data_X <- data.frame(x, f) ggplot(data = data_X, aes(x = x, y = f)) + geom_line(col = "blue", lwd = 2)+ geom_line(mapping = aes(x, x^2), lty = 3)+ labs(x = "x", y = "f(x)", title = expression("f(x)="~-x^2))
Diese invers-u-förmige Beziehung ist eine konkave Funktion. Als Eselsbrücke für das Wort konkav, welches fast das englische Wort cave enthält, können wir uns merken: eine konkave Funktion stellt eine Art Höhleneingang dar.
$c$ bewirkt eine vertikale Verschiebung der Parabel:
a <- 1; b <- 0; c <- 1 x <- seq(-2,2,0.01) f <- a*x^2 + b*x + c data_X <- data.frame(x, f) ggplot(data = data_X, aes(x = x, y = f)) + geom_line(col = "blue", lwd = 2)+ geom_line(mapping = aes(x, x^2), lty = 3)+ labs(x = "x", y = "f(x)", title = expression("f(x)="~x^2+1))
$b$ bewirkt eine horizontale und vertikale Verschiebung, die nicht mehr leicht vorhersehbar ist. Für $f(x)=x^2+x$ lässt sich beispielsweise durch Umformen leicht erkennen: $f(x)=x^2+x=x(x+1)$, dass diese Funktion zwei Nullstellen bei $0$ und $-1$ hat. Somit ist ersichtlich, dass die Funktion nach unten und links verschoben ist:
a <- 1; b <- 1; c <- 0 x <- seq(-2,2,0.01) f <- a*x^2 + b*x + c data_X <- data.frame(x, f) ggplot(data = data_X, aes(x = x, y = f)) + geom_line(col = "blue", lwd = 2)+ geom_line(mapping = aes(x, x^2), lty = 3)+ labs(x = "x", y = "f(x)", title = expression("f(x)="~x^2+x))
Für die genaue Gestalt einer allgemeinen quadratischen Funktion $ax^2 + bx + c$ würden wir die Nullstellen durch Lösen der Gleichung $ax^2 + bx + c=0$ bestimmen (via p-q Formel oder a-b-c-Formel). Den Scheitelpunkt würden wir durch Ableiten und Nullsetzen der Gleichung bestimmen. Wir müssten also $2ax+b=0$ lösen und dies in die Gleichung einsetzen. Wir könnten auch die binomischen Formeln nutzen, um die Funktion in die Gestalt $f(x):=a'(x-b')^2+c'$ oder $f(x):=a'(x-b'_1)(x-b_2')+c'$ zu bekommen, falls die Nullstellen reell sind (also das Gleichungssystem lösbar ist), da wir so die Nullstellen ablesen können als $b'$ oder $b_1'$ und $b_2'$, falls $c=0$. Für die Interpretation der Ergebnisse reicht es zu wissen, dass $a$ eine Stauchung bewirkt und entscheind dafür ist, ob die Funktion u-förmig oder invers-u-förmig verläuft.
a <- -0.5; b <- 1; c <- 2 x <- seq(-2,2,0.01) f <- a*x^2 + b*x + c data_X <- data.frame(x, f) ggplot(data = data_X, aes(x = x, y = f)) + geom_line(col = "blue", lwd = 2)+ geom_line(mapping = aes(x, x^2), lty = 3)+ labs(x = "x", y = "f(x)", title = expression("f(x)="~-0.5*x^2+x+2))
$\longrightarrow$ so ähnlich sieht die bedingte Beziehung (kontrolliert für die weiteren Prädiktoren im Modell) zwischen Bildungsabschluss der Mutter und Leseleistung aus.
Code für quadratische Verlaufsgrafik
Der Code, der die Grafik des standardisierten vorhergesagten bedingten Verlaufs des Bildungsabschlusses der Mutter erzeugt, sieht folgendermaßen aus:
X <- scale(poly(PISA2009$MotherEdu, 2)) std_par_ME <- c(0.1588, -0.1436) pred_effect_ME <- X %*% std_par_ME std_ME <- X[,1] data_ME <- data.frame(std_ME, pred_effect_ME) ggplot(data = data_ME, aes(x = std_ME, y = pred_effect_ME)) + geom_point(pch = 16, col = "blue", cex = 4)+ labs(y = "std. Leseleistung | Others", x = "std. Bildungsabschluss der Mutter | Others", title = "Standardisierte bedingte Beziehung zwischen\n Bildungsabschluss der Mutter und Leseleistung")
Wir verwenden scale, um die linearen und quadratischen Anteile des Bildungsabschlusses der Mutter zu standardisieren und speichern sie in X. Anschließend ist das Interzept der quadratischen Funktion 0 ($c=0$, da wir standardisiert haben). Die zugehörigen standardisierten Koeffizienten sind $b=0.1588$ und $a=-0.1436$, die wir aus der standardisierten summary abgelesen haben. Somit wissen wir, dass es sich um eine invers-u-förmige Beziehung handelt (ohne die Grafik zu betrachten). Wir speichern die standardisierten Koeffizienten unter std_par_ME ab und verwenden anschließend das Matrixprodukt X %*% std_par_ME, um die vorhergesagten Werte via $y_{std,i}=0.1588 ME - 0.1436ME^2$ zu berechnen. Diese vorhergesagten Werte pred_effect_ME plotten wir nun gegen die standardisierten Werte des Bildungsabschlusses der Mutter std_ME, welche in der ersten Spalte von X stehen: X[, 1].
Appendix B {#AppendixB}
Übersicht über erklärte Varianzanteile
Hier sind nochmals die Anteile erklärter Varianz der Bevölkerungsdichte über die Zeit in den vier betrachteten Modellen dargestellt:
R2 <- rbind(summary(m_l)$r.squared, summary(m_q)$r.squared, summary(m_log)$r.squared, summary(m_log_quad)$r.squared) rownames(R2) <- c("linear", "quadratisch", "exponentiell (log. Modell)", "quadratisch-exponentiell (quadratisches log. Modell)") colnames(R2) <- "R^2" round(R2, 4)
Appendix C {#AppendixC}
Übersicht über die Modelle
Die angenommenen Modell pro Messzeitpunkt $i$ sind von der Konzeption deutlich verschieden. Insbesondere der Regressionfehler ist an einer anderen Stelle:
-
Lineares Modell: $Y_i = \beta_0 + \beta_1t_i + \varepsilon_i$
-
Quadratisches Modell: $Y_i = \beta_0 + \beta_1t_i^ + \beta_2t_i^{2} + \varepsilon_i$. Hier wurde mit
poly(...,2)eine Transformation der Variable Zeit vorgenommen, was hier durch "$^*$" dargestellt werden soll. -
Logarithmisches Modell: $\text{ln}(Y_i) = \beta_0 + \beta_1t_i + \varepsilon_i$ bzw. retransformiert $Y_i = e^{\beta_0 + \beta_1t_i + \varepsilon_i} = e^{\beta_0 + \beta_1t_i} \ e^{\varepsilon_i}$. Diesem Modell ist zu entnehmen, dass sich der Fehler in Abhängigkeit von der Ausprägung der unabhängigen Variablen (hier: $t_i$) unterschiedlich stark auf das Kriterium (die abhängige Variable, hier die Bevölkerungszahl) auswirkt - er hängt also mit der Ausprägung der unabhängigen Variable zusammen. Dies ist im logarithmierten Modell nicht mehr der Fall; es handelt sich hier bis auf die Transformation der abhängigen Variable um ein "normales" Regressionsmodell: Hier verliert sich diese Beziehung zwischen unabhängiger Variable und Fehler!
-
Logarithmisches Modell mit quadratischem Term: $\text{ln}(Y_i) = \beta_0 + \beta_1t_i^ + \beta_2t_i^{2} + \varepsilon_i$ bzw. retransformiert $Y_i = e^{\beta_0 + \beta_1t_i^ + \beta_2t_i^{2} + \varepsilon_i}$. Hier wurde mit
poly(...,2)eine Transformation der Variable Zeit vorgenommen, was hier durch "$^*$" dargestellt werden soll.
Appendix D {#AppendixD}
Aufgaben
Corona-Datensatz (wie in der ggplot2-Sitzung; PsyBSc7::Sitzung_2())
Zur Vorbereitung laden Sie die Daten ein. Entweder via folgendem Code
library(ggplot2) library(car) library(MASS) ### Datensatz: Corona-Pandemie 2020 confirmed <- read.csv('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv') deaths <- read.csv('https://raw.githubusercontent.com/CSSEGISandData/COVID-19/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv') ### Long format --- confirmed_long <- reshape(confirmed, varying = names(confirmed)[-c(1:4)], v.names = 'Confirmed', timevar = 'Day', idvar = names(confirmed)[1:4], direction = 'long') deaths_long <- reshape(deaths, varying = names(deaths)[-c(1:4)], v.names = 'Deaths', timevar = 'Day', idvar = names(deaths)[1:4], direction = 'long') ### Merged data ---- long <- merge(confirmed_long, deaths_long, by = c('Province.State', 'Country.Region', 'Lat', 'Long', 'Day')) ### Full data ---- covid <- aggregate(cbind(Confirmed, Deaths) ~ Country.Region + Day, data = long, FUN = 'sum') ### Only data until April 11th ---- covid_full <- covid covid <- covid[covid$Day < 80, ] ### Subsets ---- covid_de <- covid[covid$Country.Region == 'Germany', ] covid_sel <- covid[covid$Country.Region %in% c('France', 'Germany', 'Italy', 'Spain', 'United Kingdom'), ]
oder über PsyBSc7::Aufgaben_8(). Die tagesaktuellen Daten finden Sie unter covid_full. Wir wollen uns für die Aufgaben jedoch zeitlich einschränken bis zum 11. April. Betrachten Sie in den folgenden Aufgaben bitte immer nur die Daten bis zum 11. April in covid_sel und covid_de.
Aufgabe 1. Verschaffen Sie sich einen Überblick über die Daten:
a) Stellen Sie die Verläufe der bestätigten Fälle in den Ländern Frankreich, Deutschland, Italien, Spanien und Großbritannien gegen die Zeit dar (Subdatensatz covid_sel).
b) Stellen Sie nur den Verlauf der bestätigten Fälle von Deutschland grafisch dar (Subdatensatz covid_de).
c) Stellen Sie den logarithmierten Verlauf der bestätigten Fälle von Deutschland grafisch dar.
Tipp: der Logarithmus von 0 ist nicht definiert (log(0) = -Inf); diese Werte machen in einer Regressionsanalyse keinen Sinn, weswegen Sie diese durch fehlende Werte (NA) ersetzten müssen (Subdatensatz covid_de)
d) Löschen Sie die Tage aus dem Datensatz heraus, an welchen noch keine bestätigten Fälle in Deutschland aufgetreten sind.
Tipp: is.na(...) fragt ab, ob Einträge in einem Element fehlen. Entsprechend können mit !is.na(...) jene Fälle angesprochen werden, die vorhanden sind. (Subdatensatz covid_de)
Aufgabe 2. Schauen Sie sich den Verlauf der bestätigten Infektion in Deutschland genauer an, indem Sie die folgenden Regressionsmodelle berechnen.
a) Berechnen Sie eine lineare und quadratische Regression und vergleichen Sie die beiden Modelle. (Subdatensatz covid_de)
b) Führen Sie eine lineare Regression der logarithmierten "Confirmed" Cases durch. (Subdatensatz covid_de)
c) Passt das lineare logarithmierte Modell (bzw. das exponentielle Modell) zu den Daten? Analysieren Sie die Voraussetzungen. (Subdatensatz covid_de)
d) Stellen Sie den logarithmierten sowie den retransformierten Verlauf dar. (Subdatensatz covid_de)
Aufgabe 3. Prüfen Sie nun, ob Sie die Vorhersage durch ein quadratisch-exponentielles Modell verbessern können, indem Sie einen quadratischen Trend der Zeit aufnehmen und die Modelle vergleichen.
(Subdatensatz covid_de)
Bonusaufgabe.: Vergleichen Sie die Koeffizienten des exponentiellen Wachstums zwischen den Ländern.
(Subdatensatz covid_sel)
(Die Bonusaufgabe ist komplett freiwillig und muss nicht gelöst werden, um die Studienleistung zu erbringen.)
Literatur
Eid, M., Gollwitzer, M., & Schmitt, M. (2017). Statistik und Forschungsmethoden (5. Auflage, 1. Auflage: 2010). Weinheim: Beltz.
- Blau hinterlegte Autorenangaben führen Sie direkt zur universitätsinternen Ressource.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.