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)
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
#
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.
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.
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ší.
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í.
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.