Nastavíme prostredie pre analýzu dát, zabezpečíme reprodukovateľnosť výsledkov a vyčistíme pracovná pamäť.

Naimportujeme potrebné knižnice pre štatistickú analýzu, regresné modely, prácu s dátami, tvorbu tabuliek a vizualizácií.

rm(list=ls())
library(lmtest)   #  podpora regresie
library(outliers) # analyza odlahlych hodnot (outliers)
library(gptstudio)
library(kableExtra)
library(knitr)
library(dplyr)
library(broom)
library(corrplot)
library(sandwich)
library(lmtest)

1 Príklad

Naimportujeme dataset Life expectancy a vyberieme si krajinu -> Albánsko. Vyberieme údaje pre Albánsko a vytvoríme premenné pre následnú analýzu vzťahu medzi dĺžkou života a konzumáciou alkoholu.

# import the dataset and create a data.frame udaje
#
udaje_svet <- read.csv("udaje/Life-Expectancy-Data-Updated.csv",header=TRUE,sep=",",dec=".",check.names = TRUE)
head(udaje_svet)
##      Country                        Region Year Infant_deaths Under_five_deaths
## 1    Turkiye                   Middle East 2015          11.1              13.0
## 2      Spain                European Union 2015           2.7               3.3
## 3      India                          Asia 2007          51.5              67.9
## 4     Guyana                 South America 2006          32.8              40.5
## 5     Israel                   Middle East 2012           3.4               4.3
## 6 Costa Rica Central America and Caribbean 2006           9.8              11.2
##   Adult_mortality Alcohol_consumption Hepatitis_B Measles  BMI Polio Diphtheria
## 1        105.8240                1.32          97      65 27.8    97         97
## 2         57.9025               10.35          97      94 26.0    97         97
## 3        201.0765                1.57          60      35 21.2    67         64
## 4        222.1965                5.68          93      74 25.3    92         93
## 5         57.9510                2.89          97      89 27.0    94         94
## 6         95.2200                4.19          88      86 26.4    89         89
##   Incidents_HIV GDP_per_capita Population_mln Thinness_ten_nineteen_years
## 1          0.08          11006          78.53                         4.9
## 2          0.09          25742          46.44                         0.6
## 3          0.13           1076        1183.21                        27.1
## 4          0.79           4146           0.75                         5.7
## 5          0.08          33995           7.91                         1.2
## 6          0.16           9110           4.35                         2.0
##   Thinness_five_nine_years Schooling Economy_status_Developed
## 1                      4.8       7.8                        0
## 2                      0.5       9.7                        1
## 3                     28.0       5.0                        0
## 4                      5.5       7.9                        0
## 5                      1.1      12.8                        1
## 6                      1.9       7.9                        0
##   Economy_status_Developing Life_expectancy
## 1                         1            76.5
## 2                         0            82.8
## 3                         1            65.4
## 4                         1            67.0
## 5                         0            81.7
## 6                         1            78.2
udaje_albania <- udaje_svet[udaje_svet$Country == "Albania",]
udaje_albania
##      Country         Region Year Infant_deaths Under_five_deaths
## 114  Albania Rest of Europe 2007          15.3              17.1
## 372  Albania Rest of Europe 2005          17.8              20.0
## 544  Albania Rest of Europe 2011          10.8              12.1
## 581  Albania Rest of Europe 2002          21.7              24.4
## 627  Albania Rest of Europe 2015           8.5               9.6
## 854  Albania Rest of Europe 2006          16.6              18.5
## 992  Albania Rest of Europe 2000          24.1              27.2
## 1003 Albania Rest of Europe 2013           9.3              10.4
## 1162 Albania Rest of Europe 2001          22.9              25.8
## 1615 Albania Rest of Europe 2014           8.8               9.9
## 2032 Albania Rest of Europe 2003          20.4              23.0
## 2038 Albania Rest of Europe 2008          14.1              15.8
## 2140 Albania Rest of Europe 2004          19.1              21.5
## 2164 Albania Rest of Europe 2009          12.9              14.5
## 2222 Albania Rest of Europe 2012          10.0              11.2
## 2358 Albania Rest of Europe 2010          11.8              13.3
##      Adult_mortality Alcohol_consumption Hepatitis_B Measles  BMI Polio
## 114          82.6975                5.46          98      95 25.8    99
## 372          84.4915                4.98          98      97 25.6    97
## 544          80.3490                5.03          99      99 26.2    99
## 581          87.1820                3.96          96      93 25.3    98
## 627          75.2050                4.33          99      98 26.6    99
## 854          83.5945                5.17          98      94 25.7    97
## 992          94.6955                3.92          96      90 25.2    97
## 1003         78.2430                4.28          99      99 26.4    99
## 1162         90.9390                4.51          96      90 25.2    97
## 1615         76.7240                4.10          98      98 26.5    98
## 2032         86.2855                4.33          97      93 25.4    97
## 2038         82.1105                5.56          99      98 25.9    99
## 2140         85.3885                4.42          99      96 25.5    98
## 2164         81.5235                5.79          98      98 26.0    98
## 2222         79.7620                4.43          99      99 26.3    99
## 2358         80.9365                4.88          99      98 26.1    99
##      Diphtheria Incidents_HIV GDP_per_capita Population_mln
## 114          98          0.03           3045           2.97
## 372          98          0.03           2676           3.01
## 544          99          0.03           3678           2.91
## 581          98          0.02           2247           3.05
## 627          99          0.03           3953           2.88
## 854          97          0.03           2851           2.99
## 992          97          0.01           1961           3.09
## 1003         99          0.03           3781           2.90
## 1162         97          0.01           2144           3.06
## 1615         98          0.03           3856           2.89
## 2032         97          0.02           2381           3.04
## 2038         99          0.03           3298           2.95
## 2140         97          0.02           2522           3.03
## 2164         98          0.03           3432           2.93
## 2222         99          0.03           3736           2.90
## 2358         99          0.03           3577           2.91
##      Thinness_ten_nineteen_years Thinness_five_nine_years Schooling
## 114                          1.6                      1.7       9.2
## 372                          1.8                      1.8       9.1
## 544                          1.4                      1.5       9.3
## 581                          2.0                      2.1       9.0
## 627                          1.2                      1.3       9.7
## 854                          1.7                      1.8       9.2
## 992                          2.1                      2.2       8.8
## 1003                         1.3                      1.4       9.7
## 1162                         2.1                      2.1       8.7
## 1615                         1.2                      1.3       9.7
## 2032                         1.9                      2.0       9.0
## 2038                         1.6                      1.6       9.2
## 2140                         1.8                      1.9       9.1
## 2164                         1.5                      1.6       9.3
## 2222                         1.3                      1.4       9.6
## 2358                         1.4                      1.5       9.3
##      Economy_status_Developed Economy_status_Developing Life_expectancy
## 114                         0                         1            75.6
## 372                         0                         1            75.2
## 544                         0                         1            76.9
## 581                         0                         1            74.6
## 627                         0                         1            78.0
## 854                         0                         1            75.4
## 992                         0                         1            74.0
## 1003                        0                         1            77.6
## 1162                        0                         1            74.3
## 1615                        0                         1            77.8
## 2032                        0                         1            74.8
## 2038                        0                         1            75.9
## 2140                        0                         1            75.0
## 2164                        0                         1            76.2
## 2222                        0                         1            77.3
## 2358                        0                         1            76.6
udaje_albania$y <- udaje_albania$Life_expectancy
udaje_albania$x <- udaje_albania$Alcohol_consumption

#

1.1 Klasický odhad

Týmto príkazom sa vytvoríme klasický lineárny regresný model na analýzu vzťahu medzi dĺžkou života a konzumáciou alkoholu.

fit_klasicky <- lm(y ~ x, data = udaje_albania)
summary(fit_klasicky)
## 
## Call:
## lm(formula = y ~ x, data = udaje_albania)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8701 -0.9692 -0.2836  1.0312  2.0877 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  75.4669     2.8031  26.922 1.85e-13 ***
## x             0.1029     0.5925   0.174    0.865    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.338 on 14 degrees of freedom
## Multiple R-squared:  0.002148,   Adjusted R-squared:  -0.06913 
## F-statistic: 0.03013 on 1 and 14 DF,  p-value: 0.8647

Výsledky klasického lineárneho regresného modelu ukázali, že medzi konzumáciou alkoholu a dĺžkou života v Albánsku nebol zistený štatisticky významný vzťah (p = 0,865). Koeficient determinácie (R^2 = 0,002) naznačuje, že model vysvetľuje len približne 0,2 % variability dĺžky života.

1.1.1 Test na autokoreláciu rezíduí (Durbin-Watson test)

dwtest(fit_klasicky)
## 
##  Durbin-Watson test
## 
## data:  fit_klasicky
## DW = 3.1171, p-value = 0.9949
## alternative hypothesis: true autocorrelation is greater than 0

Výsledok Durbinovho-Watsonovho testu (DW = 3,1171; p = 0,9949) naznačuje, že v reziduách regresného modelu nebola potvrdená pozitívna autokorelácia.

1.2 Odhad Koyckovej transformácie

udaje_albania$y_lag <- c(NA, head(udaje_albania$y, -1))   # posunie to o jedno obdobie dozadu
udaje2 <- na.omit(udaje_albania)                  # vynecha sa prve pozorovanie (NA)

fit <- lm(y ~ x + y_lag, data = udaje2)
summary(fit)
## 
## Call:
## lm(formula = y ~ x + y_lag, data = udaje2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.16144 -0.77411 -0.00892  0.77181  1.53252 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 120.0183    18.0369   6.654 2.35e-05 ***
## x             0.1524     0.5504   0.277   0.7866    
## y_lag        -0.5896     0.2345  -2.514   0.0272 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.164 on 12 degrees of freedom
## Multiple R-squared:  0.3489, Adjusted R-squared:  0.2403 
## F-statistic: 3.215 on 2 and 12 DF,  p-value: 0.07622

Výsledky regresného modelu ukázali, že premenná konzumácie alkoholu nebola štatisticky významná (p = 0,7866), zatiaľ čo oneskorená hodnota dĺžky života ((y_lag)) mala štatisticky významný vplyv na aktuálnu hodnotu dĺžky života (p = 0,0272). Model vysvetľoval približne 34,9 % variability závislej premennej ((R^2 = 0,3489)).

Odhadované parametre sú:

coef_fit <- coef(fit)

c_hat <- coef_fit["(Intercept)"]
beta_hat <- coef_fit["x"]
lambda_hat <- coef_fit["y_lag"]
alpha_hat <- c_hat / (1 - lambda_hat)
long_run_hat <- beta_hat / (1 - lambda_hat)

estimates <- data.frame(
  parameter = c("c", "beta", "lambda", "alpha = c/(1-lambda)", "dlhodobý efekt beta/(1-lambda)"),
  odhad = c(c_hat, beta_hat, lambda_hat, alpha_hat, long_run_hat)
)

estimates
##                        parameter        odhad
## 1                              c 120.01833989
## 2                           beta   0.15235880
## 3                         lambda  -0.58957754
## 4           alpha = c/(1-lambda)  75.50329388
## 5 dlhodobý efekt beta/(1-lambda)   0.09584861

Odhadnuté parametre modelu ukázali hodnotu interceptu (c = 120,02), koeficient pre premennú konzumácie alkoholu (= 0,152) a parameter oneskorenej premennej (= -0,590), hodnota naznačuje, že predchádzajúca hodnota dĺžky života mala v modeli určitý vplyv na súčasnú hodnotu. Vypočítaná dlhodobá rovnovážna hodnota () dosiahla približne 75,50 a dlhodobý efekt konzumácie alkoholu na dĺžku života bol približne 0,096.

Význam týchto výsledkov spočíva v tom, že model opisuje nielen okamžitý vzťah medzi konzumáciou alkoholu a dĺžkou života, ale aj vplyv predchádzajúcich hodnôt dĺžky života na aktuálny stav. Dlhodobý efekt (/(1-)) ukazuje, že pri zvýšení konzumácie alkoholu o jednu jednotku by sa dĺžka života v dlhodobom horizonte zmenila približne o 0,096 roka, pričom tento efekt je slabší.

2 Robustné štandardné chyby typu Newey-West

Pre dynamické modely časových radov je rozumné ukázať aj HAC štandardné chyby. Balíky sandwich a lmtest nie sú súčasťou úplne základnej inštalácie R, preto ich načítame len vtedy, ak sú dostupné.

  nw_vcov <- sandwich::NeweyWest(fit, lag = 4, prewhite = FALSE, adjust = TRUE)
  lmtest::coeftest(fit, vcov. = nw_vcov)
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 120.01834   12.69752  9.4521 6.56e-07 ***
## x             0.15236    0.51671  0.2949 0.773136    
## y_lag        -0.58958    0.17450 -3.3787 0.005481 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Výsledky regresného modelu s použitím Newey-Westových robustných štandardných chýb ukázali, že premenná konzumácie alkoholu nebola štatisticky významná (p = 0,773), zatiaľ čo oneskorená hodnota dĺžky života ((y_lag)) ostala štatisticky významná (p = 0,005). Použitie robustných štandardných chýb zabezpečilo spoľahlivejšie odhady pri možnej prítomnosti autokorelácie a heteroskedasticity reziduí.

3 Diagnostika rezíduí

Vkreslíme graf rezíduí transformovaného regresného modelu, ktorý slúži na diagnostiku modelu a kontrolu náhodnosti chýb v čase.

res <- residuals(fit)

plot(res, type = "l", main = "Rezíduá z transformovaného modelu",
     xlab = "čas", ylab = "rezíduum")
abline(h = 0, lty = 2)

dwtest(fit)
## 
##  Durbin-Watson test
## 
## data:  fit
## DW = 2.0191, p-value = 0.6309
## alternative hypothesis: true autocorrelation is greater than 0

Výsledok Durbinovho-Watsonovho testu (DW = 2,0191; p = 0,6309) naznačuje, že v reziduách modelu nebola zistená prítomnosť pozitívnej autokorelácie.