Pri ďalšej práci budeme používať knižnice

library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(lmtest)
library(sandwich)
library(car)
## Loading required package: carData
rm(list=ls())

1. Úvod a údaje

Údaje o očakávanej dĺžke života sú usporiadané v súbore csv, stĺpce sú oddelené znakom “;” a používajú bodku na oddelenie desatinných čisel.

1.1 Úvod do problému, stanovenie hypotéz

V tomto cvičení sa rozhodli sme modelovať počet denných prípadov Covid-19 na 100 tisíc obyvateľov cases_per_100k v závislosti od jednej vysvetľujúcej premennej, ktorou je časový trend označený ako t. Tento trend predstavuje poradové číslo dňa v časovom rade a zachytáva dlhodobý vývoj epidémie v čase.

Naša pracovná hypotéza predpokladá, že vývoj počtu prípadov má rastúci trend, teda že odhadnutý koeficient pri premennej t bude kladný. Inými slovami, očakávame, že:

čas má štatisticky významný pozitívny vplyv na počet prípadov cases_per_100k,

Okrem toho predpokladáme, že keďže pracujeme s časovým radom, reziduá môžu vykazovať autokoreláciu, čo bude predmetom ďalších testov v nasledujúcich častiach.

1.2 Príprava databázy, úprava údajov

covid <- read.csv("Covid.csv", sep = ";", header = TRUE)

# výber krajiny Slovensko a relevantných premenných

udajeSlovakia <- covid[covid$country == "Slovakia",
c("date", "cases", "population")]

# zoradenie podľa dátumu

udajeSlovakia$date <- as.Date(udajeSlovakia$date)
udajeSlovakia <- udajeSlovakia[order(udajeSlovakia$date), ]

# vytvorenie novej premennej cases_per_100k a časového trendu t

udajeSlovakia$cases_per_100k <- udajeSlovakia$cases /
udajeSlovakia$population * 100000

udajeSlovakia$t <- 1:nrow(udajeSlovakia)

# kontrolný výstup

udajeSlovakia

2. Lineárna regresia v základnom tvare

Ide o odhad rovnice

\[cases_per100k,t=β0​+β1​t+e\]

library(ggplot2)

# your regression model
model <- lm(cases_per_100k ~ t, data = udajeSlovakia)
summary(model)
## 
## Call:
## lm(formula = cases_per_100k ~ t, data = udajeSlovakia)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -9122.1 -3383.2  -250.6  4423.3  8544.0 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -8583.2990   300.5483  -28.56   <2e-16 ***
## t              39.2973     0.4759   82.57   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4965 on 1091 degrees of freedom
## Multiple R-squared:  0.862,  Adjusted R-squared:  0.8619 
## F-statistic:  6817 on 1 and 1091 DF,  p-value: < 2.2e-16

3 Autokorelácia rezíduí

V tejto časti sa pozrieme na ďalší dôležitý predpoklad klasického lineárneho regresného modelu – nezávislosť rezíduí. V časových radoch sa často stáva, že chyba v čase \(t\) je systematicky spätá s chybou v čase \(t-1\), čo nazývame autokoreláciou rezíduí.

3.1 Čo je autokorelácia rezíduí?

Autokorelácia rezíduí je situácia, keď platí:

\[ \operatorname{Corr}(e_t, e_{t-k}) \neq 0 \quad \text{pre niektoré } k \neq 0. \]

Najčastejšie sa skúma autokorelácia prvého rádu:

\[ e_t = \rho e_{t-1} + \nu_t,\quad |\rho| < 1. \]

3.2 Dôsledky autokorelácie rezíduí

  • odhady koeficientov \(\hat{\beta}\) sú nestranné,
  • ale neefektívne,
  • štandardné chyby sú podhodnotené teda p hodnoty sa javia menšie, než by mali byť - dochádza k falošnému zdaniu o štatistickej významnosti.
  • teda t- a F-testy sú skreslené.

3.3 Detekcia autokorelácie

Grafická informácia

library(ggplot2)

# pridáme vyrovnané hodnoty do databázy

udajeSlovakia$fitted <- fitted(model)

# scatterplot + fitted line

ggplot(udajeSlovakia, aes(x = date, y = cases_per_100k)) +
geom_point(color = "steelblue", size = 2) +

# regression fitted line

geom_line(aes(y = fitted), color = "red", size = 1) +

labs(
title = "Covid-19 Slovakia: Empirical Data (blue) vs. Fitted Data (red)",
x = "Dátum",
y = "Prípady na 100 tis. obyvateľov"
) +
theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Analýzou obrázka vidíme, že empirické hodnoty sa počas viacerých období držia nad alebo pod vyrovnanou hodnotou — napríklad počas nárastov či poklesov v jednotlivých vlnách pandémie. Keď sa dajú identifikovať kompaktné úseky s rovnakým znamienkom odchýlok (všetky nad alebo pod trendovou čiarou), je to typický signál možnej autokorelácie náhodných zložiek.

# uložíme si reziduá z pôvodného modelu - nazvaného model
res <- residuals(model)

ACF graf (Autocorrelation Function)

Táto funkcia priradzuje odhad korelácie, ktorá je medzi jednotlivými rezíduami v aktuálnom období a období posunutom (Lag) o \(k\) období späť.

acf(res, lag.max = 30, main = "Autokorelačná funkcia rezíduí")

Na grafe je modrou prerušovanou čiarou vyjadrený aj 95 % interval spoľahlivosti pre autokorelačný koeficient s príslušným posunom.

V našom prípade vidíme, že viaceré hodnoty autokorelácie (najmä pri nízkych lagoch) výrazne presahujú hranice intervalov spoľahlivosti. To znamená, že autokorelácia rezíduí je štatisticky významná a predpoklad nezávislosti náhodných zložiek je porušený.

Durbin–Watsonov test

Durbin–Watsonov koeficient (DW) je vypočítaný z rezíduí podľa vzorca

\[ DW = \frac{\sum_{t=2}^{n} (e_t - e_{t-1})^{2}}{\sum_{t=1}^{n} e_t^{2}} \]

kde 𝑛je počet pozorovaní. Medzi koeficientom autokorelácie dvoch susedných rezíduí a DW existuje približný vzťah

\[ \hat{\rho} \approx 1 - \frac{DW}{2} \]

Durbin–Watsonov koeficient nadobúda hodnoty od 0 po 4:

hodnoty blízko 0 znamenajú silnú pozitívnu autokoreláciu,

hodnoty okolo 2 znamenajú žiadnu autokoreláciu,

hodnoty blízko 4 indikujú negatívnu autokoreláciu.

V praxi sa často používa jednoduché pravidlo: ak DW leží približne v intervale 1,8 až 2,2, problém autokorelácie nemusíme riešiť.

V našom prípade je hodnota DW veľmi nízka, čo signalizuje prítomnosť významnej pozitívnej autokorelácie rezíduí prvého rádu.

library(lmtest)
dwtest(model)
## 
##  Durbin-Watson test
## 
## data:  model
## DW = 0.00015457, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

DW test však má obmedzenia – regresory nesmú byť časovo posunuté a model nesmie obsahovať oneskorenú vysvetľovanú premennú. Tieto obmedzenia neplatí pre Breusch–Godfreyov test, ktorý použijeme v nasledujúcom kroku.

Praktický výpočet v R
bgtest(model, order = 1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model
## LM test = 1089.9, df = 1, p-value < 2.2e-16

V našich údajoch sme identifikovali týmto testom nízku úroveň autokorelácie s posunom o 1, teda závery DW testu a BG testu sú rozporuplné. Preto aspoň z demonštratívneho dôvodu pristúpime ku odstraňovaniu dôsledkov autokorelácie.

Odhad Koyckovho modelu v R

Najjednoduchšia implementácia:

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
udajeSlovakia <- udajeSlovakia %>%
arrange(date) %>%
mutate(
cases_per_100k_lag1 = lag(cases_per_100k)
)

model_koyck <- lm(cases_per_100k ~ t + cases_per_100k_lag1,
data = udajeSlovakia)

summary(model_koyck)
## 
## Call:
## lm(formula = cases_per_100k ~ t + cases_per_100k_lag1, data = udajeSlovakia)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -53.12 -33.23 -14.71   3.31 335.02 
## 
## Coefficients:
##                       Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)         -3.6964296  4.7702085   -0.775    0.439    
## t                    0.0972863  0.0153417    6.341 3.33e-10 ***
## cases_per_100k_lag1  0.9985682  0.0003624 2755.526  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 59.43 on 1089 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 2.757e+07 on 2 and 1089 DF,  p-value: < 2.2e-16
dwtest(model_koyck)
## 
##  Durbin-Watson test
## 
## data:  model_koyck
## DW = 0.13812, p-value < 2.2e-16
## alternative hypothesis: true autocorrelation is greater than 0

Dynamizovaný model identifikuje kladný koeficient pri oneskorenej premennej cases_per_100k_lag1, ktorý je zároveň nižší ako 1. Tento koeficient zachytáva zotrvačný vplyv minulých hodnôt počtu prípadov na súčasné hodnoty — teda určitá časť zmeny prípadov sa prenáša z jedného dňa do nasledujúceho.

V našich výsledkoch sa ukazuje, že dynamizácia síce znižuje autokoreláciu rezíduí, ale nerobí žiadny z regresorov výrazne štatisticky významným. Navyše, pri porovnaní kvality modelu cez Adjusted R-squared vidíme, že pôvodný (jednoduchší) model poskytuje lepšie výsledky než tento dynamický model.

Newey–West robustné štandardné chyby

library(sandwich)
library(lmtest)

coeftest(model_koyck, vcov = NeweyWest(model_koyck))
## 
## t test of coefficients:
## 
##                       Estimate Std. Error  t value Pr(>|t|)    
## (Intercept)         -3.6964296 12.7677052  -0.2895   0.7722    
## t                    0.0972863  0.0630267   1.5436   0.1230    
## cases_per_100k_lag1  0.9985682  0.0016563 602.8980   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4. Záver

Autokorelácia rezíduí je pri modelovaní časových radov veľmi častým problémom. V našom prípade jednoduchý lineárny model s časovým trendom ukázal výraznú pozitívnu autokoreláciu, čo znamená, že predpoklad nezávislosti náhodných zložiek nie je splnený.

Jednou z možností riešenia je dynamizácia modelu — v našom cvičení prostredníctvom Koyckovho modelu, ktorý zahŕňa oneskorenú hodnotu vysvetľovanej premennej. Dynamizovaný model síce autokoreláciu zmiernil, no úplne ju neodstránil a výkon modelu (meraný Adjusted R-squared) nebol lepší než pri pôvodnom modeli.

Robustné Newey–West štandardné chyby umožňujú vykonávať spoľahlivú štatistickú inferenciu aj v prípade autokorelácie a heteroskedasticity.

V praxi sa pri dlhších časových radoch alebo pri komplexnej dynamike často používajú pokročilejšie postupy, napríklad ARIMA modely, Almonova metóda distribuovaného oneskorenia alebo postup Cochran–Orcutt, ktoré sú vhodné najmä pri výraznej sériovej závislosti dát.