1. Fragestellung

Der weltweite Plastikverbrauch wird zunehmend ein Problem (Geyer, Jambeck, and Law 2017), so dass auch die Politik handelt, siehe z.B. die Europäische Kommission - Pressemitteilung vom 21. Mai 2019.

Wie kann diese Zeitreihe modelliert werden?

2. Datenbeschaffung

Die Internetseite Our World in Data bietet dazu einen Artikel an (Ritchie and Roser 2019). Die Zeitreihe der Annual global polymer resin and fiber production (plastic production), measured in metric tonnes per year der Jahre 1950 bis 2015 lässt sich von der Seite https://ourworldindata.org/plastic-pollution#global-plastic-production als csv Datei herunterladen – und in R einlesen.

Plastic <- read.csv("global-plastics-production.csv")

Zur Vereinfachung wird die Variable Global.plastics.production..million.tonnes...tonnes. umbenannt und nur die relevanten Variablen selektiert, und in die Einheit Millionen Tonnen umgerechnet.

# Ggfs. Paket mosaic installieren
# install.packages(mosaic)

# Paket laden
library("mosaic")


Plastic <- Plastic %>%
  rename(Production = Global.plastics.production..million.tonnes...tonnes.) %>%
  mutate(Production = Production / 1e06) %>%
  select(Year, Production)

Eine erste Übersicht liefert ein Liniendiagramm:

gf_line(Production ~ Year, data = Plastic)

Man erkennt den Einbruch in der Wirtschaftskrise 2008 sowie evtl. Auswirkungen der Ölpreiskrise Mitte und Ende der 1970er Jahre.

3. Datenmodellierung

Eine lineare Modellierung über die Zeit scheint nicht sinnvoll:

gf_point(Production ~ Year, data = Plastic) %>%
  gf_lm(Production ~ Year, data = Plastic)

Hinweis: Ein lineares Modell (\(y_t \approx\beta_0 + \beta_1 \cdot t\)) würde beinhalten, dass die Änderung des Mittelwertes von \(y_t\) jedes Jahr konstant ist. Diese Annahme ist z.B. bei vielen Veränderungsprozessen nicht plausibel.¹

Betrachte dazu die absoluten Änderungen \(\delta_{y_t}=y_t-y_{t-1}\):

Plastic %>%
  mutate(delta.Production = Production - lag(Production)) %>%
  gf_point(delta.Production ~ Year) %>%
  gf_smooth()

Aber auch die relativen Änderungen \(\delta_{y_t}=\frac{y_t-y_{t-1}}{y_{t-1}}\) sind nicht konstant über die Zeit:²

Plastic %>%
  mutate(delta.Production.rel = (Production - lag(Production)) / lag(Production)) %>%
  gf_point(delta.Production.rel ~ Year) %>%
  gf_smooth()

Aufgrund der konvexen Form kann z.B. eine quadratische Funktion, d.h. \(y_t \approx \beta_0 + \beta_1 \cdot t + \beta_2 \cdot t^2\) angepasst werden:

Plastic.model <- lm(Production ~ I(Year^2) + Year, data = Plastic) 
plotModel(Plastic.model)

summary(Plastic.model)

## 
## Call:
## lm(formula = Production ~ I(Year^2) + Year, data = Plastic)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.8633  -4.3627   0.0199   3.9697  14.9621 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.001e+05  1.000e+04   40.00   <2e-16 ***
## I(Year^2)    1.046e-01  2.545e-03   41.09   <2e-16 ***
## Year        -4.091e+02  1.009e+01  -40.54   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.709 on 63 degrees of freedom
## Multiple R-squared:  0.9966, Adjusted R-squared:  0.9965 
## F-statistic:  9127 on 2 and 63 DF,  p-value: < 2.2e-16

Das Modell lautet also

\[\widehat{Production}_{Year}=4.0005517\times 10^{5} + (0.1045809)\cdot Year^2 + (-409.0843245)\cdot Year + e_{Year}\]

Auch wenn die Modellanpassung mit einem \(R^2=0.9966\) sehr gut ist, so ist zu beachten, dass insbesondere Polynomiale Modelle zur Überanpassung neigen und nicht unmittelbar zur Extrapolation geeignet sind. Auch sind in der Regel die Residuen (auto)korreliert, und somit nicht unabhängig (eine der Annahmen innerhalb der Inferenz):

gf_point(resid(Plastic.model) ~ Year, data = Plastic) %>%
  gf_hline(yintercept=0)

Eine kurzfristige Vorhersage für die Jahre 2016-2020 ergibt folgende Werte (in Millionen Tonnen):

Zeitraum <- data.frame(Year = c(2016, 2017, 2018, 2019, 2020))

predict(Plastic.model, newdata = Zeitraum, interval = "prediction")

##        fit      lwr      upr
## 1 384.8407 370.4957 399.1857
## 2 397.5313 383.0717 411.9909
## 3 410.4310 395.8458 425.0162
## 4 423.5399 408.8177 438.2621
## 5 436.8579 421.9868 451.7291

Anhang: Versionshinweise

• Datum erstellt: 2019-06-15
• R Version: 3.6.0

Verwendete Pakte:

• mosaic Version: 1.5.0
• ggformula Version: 0.9.1