06.05

2014

Złożoność modelu a dokładność prognoz

Autor: adam

Wybór optymalnego modelu dla danych to podstawowe zadanie, z którym musi zmierzyć się każdy analityk. Dopasowując model powinniśmy pamiętać przede wszystkim o ryzyku przeparametryzowania/przeuczenia modelu (ang. overfitting). Choć bardziej skomplikowany model często daje lepsze dopasowanie, nadmierna liczba parametrów prowadzi do mniej dokładnych oszacowań (estymatorów parametrów), a co gorsza, może również skutkować złą jakością prognoz.

Pomimo dość powszechnej świadomości ryzyka przeparametryzowania modelu, nierzadko nie zdajemy sobie sprawy z tego, jak poważne mogą być tego konsekwencje. Postaram się przybliżyć ten problem na prostym przykładzie szeregu czasowego zawierającego informacje o liczbie ludności Polski na przestrzeni ostatnich kilkunastu lat.

Na bazie obserwacji zarejestrowanych w latach 1995-2012 spróbujemy dopasować odpowiedni model dla trendu (model trendu wielomianowego). Następnie, wykorzystując dopasowane modele wyznaczymy prognozowaną wielkość populacji dla następnych 5-ciu lat. Oczywiście przeprowadzona analiza zawiera sporo uproszczeń i nie jest nawet próbą optymalnego prognozowania liczby ludności w Polsce. Powinna być traktowana raczej jako przykład ilustrujący wybrane kwestie metodologiczne.

1 Wprowadzenie

  • Podobnie jak dla innych danych, również w przypadku szeregów czasowych złożoność (stopień skomplikowania) modelu może mieć istotny wpływ na dokładność wyznaczonych prognoz.
  • Dotyczy to zarówno modeli dla funkcji trendu (tendencji długoterminowej), jak i modelowania korelacji/zależności czasowej danych.
  • Przeanalizujmy ten problem na przykładzie danych dotyczących populacji Polski w ostatnich kilkunastu latach.

2 Dane

  • Dane zawierają informacje o liczbie ludności Polski, zarejestrowanej w latach 1995-2012.
  • Analizowany szereg zawiera $n=19$ obserwacji. Są to dane roczne, a każda obserwacja przedstawia stan na dzień 31 XII.
  • Źródło: Dane pobrano ze strony GUS: STAN LUDNOŚCI I PROGNOZY.
  • Rysunek 1 przedstawia wykres analizowanego szeregu.
Rysunek 1: Populacja Polski w latach 1995-2012.

Rysunek 1: Populacja Polski w latach 1995-2012.

3 Modelowanie trendu

  • Spróbujmy dopasować dla tych danych prosty model postaci: $X(t) = f (t) + Z(t),$ gdzie:
    • $f(t)$ – deterministyczny trend,
    • $Z(t)$ – zakłócenie losowe typu biały szum.
  • Dla uproszczenia nie będziemy szczegółowo analizowali czy założenie o białoszumowości zakłócenia losowego jest w tym przypadku spełnione.
  • Skupimy się na porównaniu dokładności prognoz uzyskanych dla różnych (mniej lub bardziej złożonych) modeli trendu.
  • Dopasujmy dla naszego szeregu trend wielomianowy postaci: $f(t)=\beta_0+\beta_1 t+\beta_2 t^2+…+\beta_k t^k.$

3.1 Dopasowanie trendów wielomianowych

  • Zaczynamy od dopasowania trendów wielomianowych różnych stopni. Na potrzeby przykładu dopasujemy trendy stopnia: k=2 (trend kwadratowy), k=4 oraz k=6.
czas <- as.vector(time(populacja))
 
# dopasujmy trend wielomianowy roznych stopni
pop.lm2 <- lm(populacja ~ poly(czas, 2))
pop.lm4 <- lm(populacja ~ poly(czas, 4))
pop.lm6 <- lm(populacja ~ poly(czas, 6))

3.2 Porównanie dopasowania (wykres)

  • Rysunek 2 przedstawia porównanie dopasowanych funkcji trendu na tle rzeczywistych wartości szeregu.
Rysunek 2: Porównanie dopasowania trendów wielomianowych.

Rysunek 2: Porównanie dopasowania trendów wielomianowych.

  • Łatwo zauważamy, że dopasowany (najprostszy) trend kwadratowy nie w pełni odzwierciedla tendencję obserwowaną w analizowanym szeregu.
  • Dokładność dopasowania modelu trendu do danych poprawia się wraz ze wzrostem stopnia wielomianu.

3.3 Porównanie dopasowania (testy statystyczne)

  • Aby dokładniej porównać dopasowane modele trendu i wybrać najlepszy, porównamy szczegółowe informacje dotyczące jakości ich dopasowania.
  • Na początek spójrzmy na podstawowe informacje nt. modelu wykorzystując funkcję summary.
summary(pop.lm2)
## 
## Call:
## lm(formula = populacja ~ poly(czas, 2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -168244  -70539   -3756   26343  224726 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    38345618      25476 1505.18  < 2e-16 ***
## poly(czas, 2)1  -292908     108085   -2.71    0.016 *  
## poly(czas, 2)2   715210     108085    6.62  8.2e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 108000 on 15 degrees of freedom
## Multiple R-squared:  0.773,  Adjusted R-squared:  0.743 
## F-statistic: 25.6 on 2 and 15 DF,  p-value: 1.47e-05
summary(pop.lm4)
## 
## Call:
## lm(formula = populacja ~ poly(czas, 4))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -171433  -47414    6252   55582  148524 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    38345618      23321 1644.27  < 2e-16 ***
## poly(czas, 4)1  -292908      98942   -2.96    0.011 *  
## poly(czas, 4)2   715210      98942    7.23  6.7e-06 ***
## poly(czas, 4)3   182176      98942    1.84    0.089 .  
## poly(czas, 4)4  -121585      98942   -1.23    0.241    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 98900 on 13 degrees of freedom
## Multiple R-squared:  0.835,  Adjusted R-squared:  0.785 
## F-statistic: 16.5 on 4 and 13 DF,  p-value: 5.21e-05
summary(pop.lm6)
## 
## Call:
## lm(formula = populacja ~ poly(czas, 6))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -103845  -22356    1901   17009  165153 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    38345618      19009 2017.19  < 2e-16 ***
## poly(czas, 6)1  -292908      80650   -3.63   0.0039 ** 
## poly(czas, 6)2   715210      80650    8.87  2.4e-06 ***
## poly(czas, 6)3   182176      80650    2.26   0.0452 *  
## poly(czas, 6)4  -121585      80650   -1.51   0.1598    
## poly(czas, 6)5    38802      80650    0.48   0.6399    
## poly(czas, 6)6  -232828      80650   -2.89   0.0148 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 80700 on 11 degrees of freedom
## Multiple R-squared:  0.907,  Adjusted R-squared:  0.857 
## F-statistic:   18 on 6 and 11 DF,  p-value: 4.29e-05
  • Sprawdźmy teraz, wykorzystując test ANOVA, czy pomiędzy dopasowanymi modelami występują statystycznie istotne różnice.
anova(pop.lm2, pop.lm4, pop.lm6)
## Analysis of Variance Table
## 
## Model 1: populacja ~ poly(czas, 2)
## Model 2: populacja ~ poly(czas, 4)
## Model 3: populacja ~ poly(czas, 6)
##   Res.Df      RSS Df Sum of Sq    F Pr(>F)  
## 1     15 1.75e+11                           
## 2     13 1.27e+11  2  4.80e+10 3.69  0.059 .
## 3     11 7.15e+10  2  5.57e+10 4.28  0.042 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wnioski

  1. Kryteria oceniające jakość dopasowania (m.in. błąd standardowy reszt (RSE) i współczynnik dopasowania $R^2$ (adjusted R-squared)), wskazują na poprawę jakości dopasowania wraz ze wzrostem złożoności modelu.
  2. Najlepsze wartości kryteriów dopasowania otrzymujemy dla trendu wielomianowego 6-ego stopnia (k=6).
  3. Testy statystyczne (test t istotności współczynników oraz test ANOVA) nie pozwalają uznać model trendu dla k=4 za istotnie lepszy od trendu kwadratowego (k=2).
  4. Inaczej jest jednak w przypadku trendu wielomianowego stopnia k=6, dla którego jakość dopasowania jest istotnie lepsza w porównaniu do mniej złożonych modeli:
    • skorygowany współczynnik dopasowania (adjusted R2), który uwzględnia złożoność modelu wskazuje na model dla $k=6$ jako optymalny (85.7%),
    • testy istotności współczynników wskazują, że w tym modelu, współczynnik przy najwyższej potędze ($\beta_6$) jest statystycznie istotny (p-wartość=0.0148),
    • wynik testu ANOVA (p-wartość=0.042) również pozwala uznać model dla k=6 za lepiej dopasowany.

4 Prognozowanie

  • Sprawdźmy czy wnioski dotyczące porównania modeli dla trendu przełożą się na dokładność skonstruowanych prognoz.
  • Po dopasowaniu modelu dla funkcji trendu $f(t)$ możemy wyznaczyć prognozy dla kolejnych okresów ekstrapolując wartości trendu, tzn.: $\hat X(n+h)=\hat f(n+h)$ dla h=1,2,..
  • Spójrzmy zatem na prognozy dla 5-ciu kolejnych lat, wyznaczone na bazie dopasowanych modeli trendu.
  • Oprócz prognoz punktowych (point forecast) wyniki zawierają również przedziały predykcyjne, skonstruowane na poziomie ufności 80% i 95%.
library(forecast)
 
czas.h <- 2013:2017  # horyzont dla prognoz
 
# wyznaczamy prognozy na podstawie doapsowane trendu wielomianowego
prognozy.pop2 <- forecast(pop.lm2, newdata = data.frame(czas = czas.h))
prognozy.pop4 <- forecast(pop.lm4, newdata = data.frame(czas = czas.h))
prognozy.pop6 <- forecast(pop.lm6, newdata = data.frame(czas = czas.h))
 
prognozy.pop2
##   Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1       38664743 38479754 38849732 38370626 38958860
## 2       38792134 38589535 38994732 38470019 39114248
## 3       38933594 38709026 39158162 38576550 39290638
## 4       39089124 38838303 39339945 38690340 39487908
## 5       39258724 38977526 39539922 38811643 39705805
prognozy.pop4
##   Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1       38672901 38408164 38937638 38249303 39096498
## 2       38741980 38310041 39173920 38050847 39433113
## 3       38778820 38090555 39467086 37677547 39880093
## 4       38769944 37723537 39816351 37095620 40444269
## 5       38700433 37177320 40223547 36263345 41137522
prognozy.pop6
##   Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1       38091111 37663748 38518475 37401219 38781004
## 2       36829457 35657886 38001027 34938191 38720723
## 3       34171322 31532139 36810506 29910890 38431755
## 4       29308775 24089310 34528239 20882995 37734555
## 5       21157540 11717222 30597858  5918037 36397042
  • Porównajmy jeszcze na jednym wykresie prognozy wyznaczone na bazie trzech modeli trendu (patrz Rysunek 3).
Rysunek 3: Porównanie dokładności prognoz dla różnych modeli trendu.

Rysunek 3: Porównanie dokładności prognoz dla różnych modeli trendu.

Wnioski

  1. W przypadku trendu kwadratowego (k=2) i trendu wielomianowego 4-ego stopnia, prognozy krótkoterminowe (dla 2-3 najbliższych lat) są stosunkowe podobne. Oczywiście, zbadanie na ile dokładne są to prognozy wymagałaby przeprowadzenia dodatkowych analiz.
  2. Prognozy uzyskane na bazie — uznanego za najlepiej dopasowany — trendu stopnia k=6 są całkowicie nie do przyjęcia! Obserwujemy gwałtowny spadek wielkości populacji w kolejnych okresach.

5 Podsumowanie

  • Dopasowanie zbyt złożonego modelu do danych (ang. overfitting) może mieć bardzo niekorzystny wpływ na dokładność prognoz.
  • Mając do wyboru kilka modeli, o podobnej jakości dopasowania, na potrzeby konstrukcji prognoz powinniśmy wybrać model prostszy.
  • W ten sposób otrzymamy model łatwiejszy w interpretacji oraz zmniejszymy ryzyko uzyskania błędnych prognoz.
  • Aby ograniczyć ryzyko przeparametryzowania modelu możemy wykorzystać także ocenę dokładności prognoz na bazie (niezależnego) zbioru testowego.

Spróbuj ponownie