Regression V: nichtlineare Regression - exponentielles Wachstum

Welche der folgenden Aussagen ist richtig?

• Die Gleichung \(10^x = 100\) können wir bspw. mit dem 10er-Logarithmus \(log_{10}(100)\) berechnen. Die Antwort ist 2.
• Die Gleichung \(10^x = 100\) können wir bspw. mit dem 10er-Logarithmus \(log_{10}(100)\) berechnen. Die Antwort ist 3.
• Logarithmen können NICHT in einander umgerechnet werden.
• Es gilt: \(log(a+b)=log(a)log(b)\).
• Es gilt: \(log(ab)=log(a)+log(b)\).
• \(log_a(x)\) ist die Umkehrfunktion zu \(a^x\). Dies bedeutet, dass \(log_a(a^x)=x.\)
• Es gilt: \(a^{c+d}=a^c+a^d\).
• Es gilt: \(a^{c+d}=a^ca^d\).
• Es gilt immer \(a^x < b^x\) für a < b und alle x.
• Es gilt immer \(a^x < b^x\) für 1 < a < b und alle x.
• Es gilt immer \(a^x < b^x\) für 0 < a < b und alle x > 0.
• Es gilt immer \(a^x > b^x\) für 0 < a < b < 1 und alle x.
• Es gilt immer \(a^x < b^x\) für 0 < a < b < 1 und alle x.
• Es gilt immer \(a^x > b^x\) für 0 < a < b und alle x < 0.

In fett, kursiv und mit einem werden Ihnen die richtigen Antworten angezeigt:

• Die Gleichung \(10^x = 100\) können wir bspw. mit dem 10er-Logarithmus \(log_{10}(100)\) berechnen. Die Antwort ist 2.
• Die Gleichung \(10^x = 100\) können wir bspw. mit dem 10er-Logarithmus \(log_{10}(100)\) berechnen. Die Antwort ist 3.
• Logarithmen können NICHT in einander umgerechnet werden.
• Es gilt: \(log(a+b)=log(a)log(b)\).
• Es gilt: \(log(ab)=log(a)+log(b)\).
• \(log_a(x)\) ist die Umkehrfunktion zu \(a^x\). Dies bedeutet, dass \(log_a(a^x)=x.\)
• Es gilt: \(a^{c+d}=a^c+a^d\).
• Es gilt: \(a^{c+d}=a^ca^d\).
• Es gilt immer \(a^x < b^x\) für a < b und alle x.
• Es gilt immer \(a^x < b^x\) für 1 < a < b und alle x.
• Es gilt immer \(a^x < b^x\) für 0 < a < b und alle x > 0.
• Es gilt immer \(a^x > b^x\) für 0 < a < b < 1 und alle x.
• Es gilt immer \(a^x < b^x\) für 0 < a < b < 1 und alle x.
• Es gilt immer \(a^x > b^x\) für 0 < a < b und alle x < 0.

Vielleicht ist in den Daten anstatt eines exponentiellen Trends 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 <- lm(Population ~ poly(Year,2), data = WorldPopulation) # quadratischer Verlauf
summary(m_q)

## 
## Call:
## lm(formula = Population ~ poly(Year, 2), data = WorldPopulation)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -646432518 -267430205   89039661  288698332  475103733 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    2.617e+09  2.120e+07  123.47   <2e-16 ***
## poly(Year, 2)1 2.437e+10  3.152e+08   77.32   <2e-16 ***
## poly(Year, 2)2 1.238e+10  3.152e+08   39.30   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 315200000 on 218 degrees of freedom
## Multiple R-squared:  0.9718,	Adjusted R-squared:  0.9716 
## F-statistic:  3761 on 2 and 218 DF,  p-value: < 2.2e-16

Das Varianzinkrement im Vergleich zum linearen Verlauf beläuft sich auf:

summary(m_q)$r.squared - summary(m_l)$r.squared  # Inkrement 

## [1] 0.1995076

Da der Parameter des quadratischen Anteils signifikant ist, gehen wir davon aus, dass das Modell durch den quadratischen Verlauf verbessert wird. Wir sichern dies noch einmal mit einem Modellvergleich ab:

anova(m_l, m_q)

## Analysis of Variance Table
## 
## Model 1: Population ~ Year
## Model 2: Population ~ poly(Year, 2)
##   Res.Df        RSS Df  Sum of Sq      F    Pr(>F)    
## 1    219 1.7504e+20                                   
## 2    218 2.1652e+19  1 1.5338e+20 1544.3 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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!

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 
     geom_smooth(method="lm", formula = "y~poly(x,2)", col = "darkblue")  # plotte quadratischen Verlauf

Durch den quadratischen Verlauf lassen sich 97.18% der Variation der Bevölkerungsdichte erklären, was einem signifikantem Varianzinkrement von 19.95% 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.

Den quadratischen Verlauf können wir, wie auch den linearen, nicht inferenzstatistisch gegen den exponentiellen Verlauf testen. Dies liegt daran, dass die Modelle nicht auseinander hervorgehen (sie sind nicht geschachtelt). Allerdings war es so, dass das lineare Modell für die logarithmierte Bevölkerungsdichte zwar gut zu den Daten zu passen schien. Insgesamt können 92.98 % der Variation in den Daten durch den Zeitverlauf erklärt werden — das war allerdings etwas weniger als durch das quadratische Wachstum, durch welches 97.18 % der Variation in den Daten durch den Zeitverlauf erklärt werden konnten.

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.

Die Vorhersage mit dem diskontinuierlichen Modell exponentiellen Wachstums ist immer noch nicht perfekt. Insbesondere der vorhergesagte sprunghafte Anstieg von 1950 zu 1951 scheint unrealistisch. Inspizieren wir noch einmal den Residuenplot aus dem exponentiellen Modell:

residualPlot(m_log, col = "red") # Residualplot

Dem Residualplot ist zu entnehmen, dass ggf. sogar ein nicht-linearer Verlauf angemessen wäre (wobei wir uns hier auf die logarithmische Skala beziehen, die Ursprungsdaten werden natürlich bereits nicht-linear, nämlich exponentiell, modelliert). Wir gehen dieser Vermutung nach, indem wir dem Log-Plot einen quadratischen Verlauf hinzufügen. Auch dem geom_smooth Befehl kann einfach das modifizierte Regressionsmodell übergeben werden. Dem Argument formula müssen wir dazu jediglich die Formel für die Regression übergeben, wobei y das Kriterium und poly(x,2) (dies hatten wir in der vergangen Sitzung zu quadratischer oder moderierter Regression kennengelernt) der Prädiktor mit linearem und quadratischen Anteil beschreibt. Mit col = "gold3" wird die neue Linie auch extra eingefärbt.

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 in das logarithmierte Modell einen quadratischen Trend der Zeit ein.

Dies wollen wir schätzen, indem wir im Regressionsmodell der logarithmierten Populationsgröße einen quadratischen Trend (poly(...,2)) für die Zeit definieren. Das Modell nennen wir, mein unbegrenztes Einfallsreichtum beweisend, m_log_quad:

m_log_quad <- lm(log_Population ~ poly(Year, 2), data = WorldPopulation) # lineares Modell mit log(y) als AV (logarithmische Skala)
summary(m_log_quad)

## 
## Call:
## lm(formula = log_Population ~ poly(Year, 2), data = WorldPopulation)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.097348 -0.037302  0.008814  0.039090  0.070984 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    21.478533   0.003111 6903.34   <2e-16 ***
## poly(Year, 2)1  8.787459   0.046253  189.99   <2e-16 ***
## poly(Year, 2)2  2.315495   0.046253   50.06   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04625 on 218 degrees of freedom
## Multiple R-squared:  0.9944,	Adjusted R-squared:  0.9943 
## F-statistic: 1.93e+04 on 2 and 218 DF,  p-value: < 2.2e-16

Wir können dieses Modell nun gegen m_log testen und so die Modellverbesserung des quadratisch exponentiellen Verlaufs auf Signifikanz prüfen (und zwar mit dem anova-Befehl).

anova(m_log, m_log_quad)

## Analysis of Variance Table
## 
## Model 1: log_Population ~ Year
## Model 2: log_Population ~ poly(Year, 2)
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    219 5.8279                                  
## 2    218 0.4664  1    5.3615 2506.1 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Durch den quadratisch-exponentiellen Verlauf (bzw. den quadratischen Verlauf in den logarithmierten Daten) lassen sich 99.44% der Variation der Bevölkerungsdichte erklären, was einem signifikantem Varianzinkrement von 6.46% 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 Modellen ist in Appendix C einzusehen.

Schauen wir uns den Residuenplot des quadratisch-exponentiellen Modells an (in logarithmierter Skala, in welcher das Modell auch geschätzt wurde!).

residualPlot(m_log_quad)

Es ist deutlich zu sehen, dass der Knick um 1950 auch im quadratisch-exponentiellen Modell noch zu sehen ist. Die Residuenplots zeigen außerdem, dass auch hier die Residuen nicht vollständig unsystematisch sind. Dennoch ist die Resiudalvarianz sehr klein. Dies erkennen wir an der y-Achse.

Stellen wir dies grafisch dar. Dazu müssen wir wieder die Vorhersage abspeichern. Diese können wir einem lm-Objekt mit predict entlocken. Wir nennen die neue Variable einfach pred_Pop_exp_quad.

WorldPopulation$pred_Pop_exp_quad <- exp(predict(m_log_quad)) # Abspeichern der retransformierten vorhergesagten Werten (wieder auf der Skala der Weltbevölkerung)
head(WorldPopulation)

##   Year Population log_Population pred_Pop_exp Post1950 pred_Pop_dis
## 1 1800  982454635       20.70556    768019076        0    940385460
## 2 1801  985895135       20.70906    775168245        0    946020970
## 3 1802  989555435       20.71277    782383962        0    951690251
## 4 1803  992954135       20.71620    789666848        0    957393507
## 5 1804  995485735       20.71874    797017526        0    963130942
## 6 1805 1000098535       20.72336    804436629        0    968902759
##   pred_Pop_exp_quad
## 1        1082904426
## 2        1082790682
## 3        1082769604
## 4        1082841189
## 5        1083005453
## 6        1083262439

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, exponentiellen, diskontinuierlichen und dem quadratisch-exponentiellen Verlauf dargestellt. Außerdem

ggplot(data = WorldPopulation, aes(x = Year, y = Population))+
     geom_point()+geom_smooth(method="lm", formula = "y~x")+         # plotte linearen 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) +
     geom_line(aes(x = Year, y = pred_Pop_dis), col = 'aquamarine3', lwd = 1.5)

# 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")

Was bedeuten nun die Parameter in unserem quadratisch-exponentiellen Modell? Der Regressionskoeffizient des linearen Trends von poly liegt bei 8.79 und der Koeffizient des quadratischen Trends bei 2.32. Dies spricht für ein beschleunigtes exponentielles Wachstum. Nach diesem Vorhersagemodell scheinen die Menschen sich schneller als exponentiell zu vermehren (zumindest von 1800 bis 2020). Eine Übersicht über die Modelle sehen Sie in Appendix D.

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 (siehe Appendix A). 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 7.88\(*10^{11}\) vorhersagen, was absolut nicht sinnvoll ist.

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_dis)$r.squared,
            summary(m_log_quad)$r.squared)
rownames(R2) <- c("linear", "quadratisch", "exponentiell (log. Modell)", "diskontinuierlich (log. Modell)", "quadratisch-exponentiell (quadratisches log. Modell)")
colnames(R2) <- "R^2"
round(R2, 4)

##                                                         R^2
## linear                                               0.7723
## quadratisch                                          0.9718
## exponentiell (log. Modell)                           0.9298
## diskontinuierlich (log. Modell)                      0.9961
## quadratisch-exponentiell (quadratisches log. Modell) 0.9944

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} e^{\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, diskontinuierliches Modell \(\text{ln}(Y_i) = \beta_0 + \beta_1t_i + \beta_2d_i + \beta_3 (t_i \cdot d_i) + \varepsilon_i\) In diesem Modell wird das logarithmische Modell um eine potentiell sprunghafte Veränderung (\(\beta_2\)) und eine veränderte Wachstumsrate (\(\beta_3\)) erweitert, um der Veränderung des Verlaufs ab 1950 Rechnung zu tragen. Die Umsetzung setzt voraus, dass im Datensatz ein dummy-kodierter (0/1) Prädiktor angelegt wurde, der die Zeit vor und nach dem Ereignis kodiert. Diese dummy-Variable (\(d_i\)) kann dann als regulärer Prädiktor mit Interaktionsterm aufgenommen werden.

• 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.