Import súboru

# 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)
 udaje_svet <- udaje_svet[-992,]
# z databázy udaje_svet si vyberieme len tie pozorovania, ktoré sa týkajú Moldavsko 
udaje <- subset(udaje_svet, Country == "Moldova")
# vyrovnanie priebehu očakávanej dĺžky dožitia v čase
model <- lm(Life_expectancy ~ Alcohol_consumption+Adult_mortality+Incidents_HIV,data = udaje)
library(broom)
library(knitr)
library(kableExtra)

# koeficienty regresie
tidy(model) %>%
  kable(digits = 3, caption = "Odhadnuté koeficienty regresie") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Odhadnuté koeficienty regresie
term estimate std.error statistic p.value
(Intercept) 80.406 0.629 127.796 0.000
Alcohol_consumption -0.009 0.039 -0.242 0.813
Adult_mortality -0.056 0.002 -27.784 0.000
Incidents_HIV 1.342 0.900 1.492 0.162

# kvalita vyrovnania
glance(model) %>%
  kable(digits = 3, caption = "Ukazovatele kvality vyrovnania") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Ukazovatele kvality vyrovnania
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.988 0.985 0.194 318.592 0 3 5.853 -1.706 2.157 0.451 12 16
NA
# výber premenných
X <- udaje[, c("Alcohol_consumption", "Adult_mortality", "Incidents_HIV")]

# výpočet korelačnej matice
cor_matrix <- cor(X, use = "complete.obs")

# zaokrúhlenie
round(cor_matrix, 4)
                    Alcohol_consumption Adult_mortality Incidents_HIV
Alcohol_consumption              1.0000         -0.0964        0.3940
Adult_mortality                 -0.0964          1.0000       -0.3885
Incidents_HIV                    0.3940         -0.3885        1.0000
library(knitr)
round(cor_matrix, 4) %>%
  kable(caption = "Korelačná matica")
Korelačná matica
Alcohol_consumption Adult_mortality Incidents_HIV
Alcohol_consumption 1.0000 -0.0964 0.3940
Adult_mortality -0.0964 1.0000 -0.3885
Incidents_HIV 0.3940 -0.3885 1.0000
cor.test (udaje$Adult_mortality, udaje$Alcohol_consumption)

    Pearson's product-moment correlation

data:  udaje$Adult_mortality and udaje$Alcohol_consumption
t = -0.36249, df = 14, p-value = 0.7224
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.5651204  0.4193211
sample estimates:
        cor 
-0.09642769 
cor.test (udaje$Alcohol_consumption, udaje$Incidents_HIV)

    Pearson's product-moment correlation

data:  udaje$Alcohol_consumption and udaje$Incidents_HIV
t = 1.6041, df = 14, p-value = 0.131
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1263452  0.7443520
sample estimates:
      cor 
0.3940387 
cor.test (udaje$Adult_mortality, udaje$Incidents_HIV )

    Pearson's product-moment correlation

data:  udaje$Adult_mortality and udaje$Incidents_HIV
t = -1.5777, df = 14, p-value = 0.137
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.7414342  0.1327472
sample estimates:
       cor 
-0.3885244 

Multicollinearity detection and testing

if (!requireNamespace("corrplot", quietly = TRUE)) {
  install.packages("corrplot")}

library(corrplot)

corrplot(cor_matrix, method = "number", type = "upper")

Variance Inflation Factor

Pre premennú \(x_j\) definujeme

\[ VIF_j = \frac{1}{1-R_j^2}, \]

kde \(R_j^2\) je koeficient determinácie z pomocnej regresie, v ktorej je \(x_j\) vysvetľovaná ostatnými regresormi.

Ak je \(R_j^2\) blízko jednej, potom je \(VIF_j\) veľký a premenná \(x_j\) je silno lineárne vysvetliteľná ostatnými premennými.

install.packages(“car”) # run once if not installed

library(car)
# Variance Inflation Factors
vif_values <- vif(model)

vif_values
Alcohol_consumption     Adult_mortality       Incidents_HIV 
           1.189129            1.183085            1.387520 

Interpretation:

Číslo podmienenosti

Číslo podmienenosti je založené na vlastných číslach matice \(X'X\). Ak sú vlastné čísla veľmi rozdielne, matica je zle podmienená.

\[ \kappa = \sqrt{\frac{\lambda_{\max}}{\lambda_{\min}}}. \]

X_scaled <- scale(X)
eigen_values <- eigen(cor(X_scaled))$values

condition_number <- sqrt(max(eigen_values) / min(eigen_values))
eigen_values
[1] 1.6036759 0.9035825 0.4927416
condition_number
[1] 1.804051

Interpretácia čísla podmienenosti: - \(\kappa \approx 1\): žiadna multikolinearita, - \(\kappa > 10\): mierna multikolinearita, - \(\kappa > 30\): silná multikolinearita.

Malá p-hodnota znamená, že zamietame hypotézu \(R=I_k\), teda medzi vysvetľujúcimi premennými existuje štatisticky významná korelačná štruktúra.

Eliminácia problému multikolinearity

Multikolinearita nie je porušením exogenity. Nie je teda automaticky dôvodom na zamietnutie OLS modelu. Problém je hlavne inferenčný: veľké štandardné chyby a nestabilné individuálne koeficienty.

Možné riešenia:

  1. získať viac dát,
  2. odstrániť alebo spojiť silno korelované premenné,
  3. transformovať premenné,
  4. použiť teoretické obmedzenia,
  5. použiť regularizačné metódy, napríklad ridge regresiu.

Odstránenie jednej z korelovaných premenných

Ak sú dve premenné takmer rovnaké a ekonomická teória nevyžaduje obe, môžeme jednu z nich vynechať.

model_reduced <- lm(Life_expectancy ~ Alcohol_consumption + Incidents_HIV, data = udaje)

summary(model)

Call:
lm(formula = Life_expectancy ~ Alcohol_consumption + Adult_mortality + 
    Incidents_HIV, data = udaje)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.30738 -0.07191  0.00501  0.06394  0.33087 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)         80.406412   0.629177 127.796  < 2e-16 ***
Alcohol_consumption -0.009374   0.038803  -0.242    0.813    
Adult_mortality     -0.056383   0.002029 -27.784 2.92e-12 ***
Incidents_HIV        1.342110   0.899627   1.492    0.162    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1938 on 12 degrees of freedom
Multiple R-squared:  0.9876,    Adjusted R-squared:  0.9845 
F-statistic: 318.6 on 3 and 12 DF,  p-value: 1.06e-11
summary(model_reduced)

Call:
lm(formula = Life_expectancy ~ Alcohol_consumption + Incidents_HIV, 
    data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4850 -1.0324 -0.5959  0.9277  2.5740 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         65.54601    2.57299  25.475 1.76e-12 ***
Alcohol_consumption -0.08151    0.30065  -0.271    0.791    
Incidents_HIV       10.91938    6.45281   1.692    0.114    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.505 on 13 degrees of freedom
Multiple R-squared:   0.19, Adjusted R-squared:  0.06535 
F-statistic: 1.524 on 2 and 13 DF,  p-value: 0.2542
LS0tCnRpdGxlOiAiTXVsdGlrb2xpbmVhcml0YSB2IGxpbmXDoXJuZWogcmVncmVzaWkiCmF1dGhvcjogIk5hdMOhbGlhIEt1bm92w6EiCm91dHB1dDoKICBodG1sX25vdGVib29rOgogICAgdG9jOiB0cnVlCiAgICB0b2NfZGVwdGg6IDMKLS0tCiMjICBJbXBvcnQgc8O6Ym9ydQpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFLCB3YXJuaW5nID0gRkFMU0UsIG1lc3NhZ2UgPSBGQUxTRSkKc2V0LnNlZWQoMTIzKQpgYGAKCmBgYHtyfQojIGltcG9ydCB0aGUgZGF0YXNldCBhbmQgY3JlYXRlIGEgZGF0YS5mcmFtZSB1ZGFqZQp1ZGFqZV9zdmV0IDwtIHJlYWQuY3N2KCJ1ZGFqZS9MaWZlLUV4cGVjdGFuY3ktRGF0YS1VcGRhdGVkLmNzdiIsaGVhZGVyPVRSVUUsc2VwPSIsIixkZWM9Ii4iLGNoZWNrLm5hbWVzID0gVFJVRSkKaGVhZCh1ZGFqZV9zdmV0KQogdWRhamVfc3ZldCA8LSB1ZGFqZV9zdmV0Wy05OTIsXQpgYGAKYGBge3J9CiMgeiBkYXRhYsOhenkgdWRhamVfc3ZldCBzaSB2eWJlcmllbWUgbGVuIHRpZSBwb3pvcm92YW5pYSwga3RvcsOpIHNhIHTDvWthasO6IE1vbGRhdnNrbyAKdWRhamUgPC0gc3Vic2V0KHVkYWplX3N2ZXQsIENvdW50cnkgPT0gIk1vbGRvdmEiKQpgYGAKCmBgYHtyfQojIHZ5cm92bmFuaWUgcHJpZWJlaHUgb8SNYWvDoXZhbmVqIGTEusW+a3kgZG/Fvml0aWEgdiDEjWFzZQptb2RlbCA8LSBsbShMaWZlX2V4cGVjdGFuY3kgfiBBbGNvaG9sX2NvbnN1bXB0aW9uK0FkdWx0X21vcnRhbGl0eStJbmNpZGVudHNfSElWLGRhdGEgPSB1ZGFqZSkKbGlicmFyeShicm9vbSkKbGlicmFyeShrbml0cikKbGlicmFyeShrYWJsZUV4dHJhKQoKIyBrb2VmaWNpZW50eSByZWdyZXNpZQp0aWR5KG1vZGVsKSAlPiUKICBrYWJsZShkaWdpdHMgPSAzLCBjYXB0aW9uID0gIk9kaGFkbnV0w6kga29lZmljaWVudHkgcmVncmVzaWUiKSAlPiUKICBrYWJsZV9zdHlsaW5nKGJvb3RzdHJhcF9vcHRpb25zID0gInN0cmlwZWQiLCBmdWxsX3dpZHRoID0gRkFMU0UpCgojIGt2YWxpdGEgdnlyb3ZuYW5pYQpnbGFuY2UobW9kZWwpICU+JQogIGthYmxlKGRpZ2l0cyA9IDMsIGNhcHRpb24gPSAiVWthem92YXRlbGUga3ZhbGl0eSB2eXJvdm5hbmlhIikgJT4lCiAga2FibGVfc3R5bGluZyhib290c3RyYXBfb3B0aW9ucyA9ICJzdHJpcGVkIiwgZnVsbF93aWR0aCA9IEZBTFNFKQoKYGBgCgpgYGB7cn0KIyB2w71iZXIgcHJlbWVubsO9Y2gKWCA8LSB1ZGFqZVssIGMoIkFsY29ob2xfY29uc3VtcHRpb24iLCAiQWR1bHRfbW9ydGFsaXR5IiwgIkluY2lkZW50c19ISVYiKV0KCiMgdsO9cG/EjWV0IGtvcmVsYcSNbmVqIG1hdGljZQpjb3JfbWF0cml4IDwtIGNvcihYLCB1c2UgPSAiY29tcGxldGUub2JzIikKCiMgemFva3LDumhsZW5pZQpyb3VuZChjb3JfbWF0cml4LCA0KQpgYGAKCmBgYHtyfQpsaWJyYXJ5KGtuaXRyKQpyb3VuZChjb3JfbWF0cml4LCA0KSAlPiUKICBrYWJsZShjYXB0aW9uID0gIktvcmVsYcSNbsOhIG1hdGljYSIpCmBgYAoKYGBge3J9CmNvci50ZXN0ICh1ZGFqZSRBZHVsdF9tb3J0YWxpdHksIHVkYWplJEFsY29ob2xfY29uc3VtcHRpb24pCmNvci50ZXN0ICh1ZGFqZSRBbGNvaG9sX2NvbnN1bXB0aW9uLCB1ZGFqZSRJbmNpZGVudHNfSElWKQpjb3IudGVzdCAodWRhamUkQWR1bHRfbW9ydGFsaXR5LCB1ZGFqZSRJbmNpZGVudHNfSElWICkKYGBgCgojIE11bHRpY29sbGluZWFyaXR5IGRldGVjdGlvbiBhbmQgdGVzdGluZwoKYGBge3IgY29ycmVsYXRpb24tcGxvdH0KaWYgKCFyZXF1aXJlTmFtZXNwYWNlKCJjb3JycGxvdCIsIHF1aWV0bHkgPSBUUlVFKSkgewogIGluc3RhbGwucGFja2FnZXMoImNvcnJwbG90Iil9CgpsaWJyYXJ5KGNvcnJwbG90KQoKY29ycnBsb3QoY29yX21hdHJpeCwgbWV0aG9kID0gIm51bWJlciIsIHR5cGUgPSAidXBwZXIiKQpgYGAKCiMjICBWYXJpYW5jZSBJbmZsYXRpb24gRmFjdG9yCgpQcmUgcHJlbWVubsO6IFwoeF9qXCkgZGVmaW51amVtZQoKXFsKVklGX2ogPSBcZnJhY3sxfXsxLVJfal4yfSwKXF0KCmtkZSBcKFJfal4yXCkgamUga29lZmljaWVudCBkZXRlcm1pbsOhY2llIHogcG9tb2NuZWogcmVncmVzaWUsIHYga3RvcmVqIGplIFwoeF9qXCkgdnlzdmV0xL5vdmFuw6Egb3N0YXRuw71taSByZWdyZXNvcm1pLgoKQWsgamUgXChSX2peMlwpIGJsw616a28gamVkbmVqLCBwb3RvbSBqZSBcKFZJRl9qXCkgdmXEvmvDvSBhIHByZW1lbm7DoSBcKHhfalwpIGplIHNpbG5vIGxpbmXDoXJuZSB2eXN2ZXRsaXRlxL5uw6Egb3N0YXRuw71taSBwcmVtZW5uw71taS4KCiMgaW5zdGFsbC5wYWNrYWdlcygiY2FyIikgICMgcnVuIG9uY2UgaWYgbm90IGluc3RhbGxlZAoKCmBgYHtyfQpsaWJyYXJ5KGNhcikKIyBWYXJpYW5jZSBJbmZsYXRpb24gRmFjdG9ycwp2aWZfdmFsdWVzIDwtIHZpZihtb2RlbCkKCnZpZl92YWx1ZXMKYGBgCgoqKkludGVycHJldGF0aW9uOioqCgotIFwoVklGX2ogPSAxXCk6IMW+aWFkbmEgbXVsdGlrb2xpbmVhcml0YSwKLSBcKDEgPCBWSUZfaiA8IDVcKTogbWllcm5hIG11bHRpa29saW5lYXJpdGEKLSBcKFZJRl9qIFxnZXEgNVwpOiBzaWxuw6EgbXVsdGlrb2xpbmVhcml0YSAobmlla3RvcsOtIGF1dG9yaSBwb3XFvsOtdmFqw7ogcHLDrXNuZWrFocOtIHByYWggXChWSUZfaiBcZ2VxIDEwXCkpLgoKIyMgIMSMw61zbG8gcG9kbWllbmVub3N0aQoKxIzDrXNsbyBwb2RtaWVuZW5vc3RpIGplIHphbG/FvmVuw6kgbmEgdmxhc3Ruw71jaCDEjcOtc2xhY2ggbWF0aWNlIFwoWCdYXCkuIEFrIHPDuiB2bGFzdG7DqSDEjcOtc2xhIHZlxL5taSByb3pkaWVsbmUsIG1hdGljYSBqZSB6bGUgcG9kbWllbmVuw6EuCgpcWwpca2FwcGEgPSBcc3FydHtcZnJhY3tcbGFtYmRhX3tcbWF4fX17XGxhbWJkYV97XG1pbn19fS4KXF0KCmBgYHtyIGNvbmRpdGlvbi1udW1iZXJ9Clhfc2NhbGVkIDwtIHNjYWxlKFgpCmVpZ2VuX3ZhbHVlcyA8LSBlaWdlbihjb3IoWF9zY2FsZWQpKSR2YWx1ZXMKCmNvbmRpdGlvbl9udW1iZXIgPC0gc3FydChtYXgoZWlnZW5fdmFsdWVzKSAvIG1pbihlaWdlbl92YWx1ZXMpKQplaWdlbl92YWx1ZXMKY29uZGl0aW9uX251bWJlcgpgYGAKICAKKipJbnRlcnByZXTDoWNpYSDEjcOtc2xhIHBvZG1pZW5lbm9zdGk6KiogCi0gXChca2FwcGEgXGFwcHJveCAxXCk6IMW+aWFkbmEgbXVsdGlrb2xpbmVhcml0YSwKLSBcKFxrYXBwYSA+IDEwXCk6IG1pZXJuYSBtdWx0aWtvbGluZWFyaXRhLAotIFwoXGthcHBhID4gMzBcKTogc2lsbsOhIG11bHRpa29saW5lYXJpdGEuCgogIApNYWzDoSBwLWhvZG5vdGEgem5hbWVuw6EsIMW+ZSB6YW1pZXRhbWUgaHlwb3TDqXp1IFwoUj1JX2tcKSwgdGVkYSBtZWR6aSB2eXN2ZXTEvnVqw7pjaW1pIHByZW1lbm7DvW1pIGV4aXN0dWplIMWhdGF0aXN0aWNreSB2w716bmFtbsOhIGtvcmVsYcSNbsOhIMWhdHJ1a3TDunJhLgoKCiMgIEVsaW1pbsOhY2lhIHByb2Jsw6ltdSBtdWx0aWtvbGluZWFyaXR5CgpNdWx0aWtvbGluZWFyaXRhIG5pZSBqZSBwb3J1xaFlbsOtbSBleG9nZW5pdHkuIE5pZSBqZSB0ZWRhIGF1dG9tYXRpY2t5IGTDtHZvZG9tIG5hIHphbWlldG51dGllIE9MUyBtb2RlbHUuIFByb2Jsw6ltIGplIGhsYXZuZSBpbmZlcmVuxI1uw706IHZlxL5rw6kgxaF0YW5kYXJkbsOpIGNoeWJ5IGEgbmVzdGFiaWxuw6kgaW5kaXZpZHXDoWxuZSBrb2VmaWNpZW50eS4KCk1vxb5uw6kgcmllxaFlbmlhOgoKMS4gesOtc2thxaUgdmlhYyBkw6F0LAoyLiBvZHN0csOhbmnFpSBhbGVibyBzcG9qacWlIHNpbG5vIGtvcmVsb3ZhbsOpIHByZW1lbm7DqSwKMy4gdHJhbnNmb3Jtb3ZhxaUgcHJlbWVubsOpLAo0LiBwb3XFvmnFpSB0ZW9yZXRpY2vDqSBvYm1lZHplbmlhLAo1LiBwb3XFvmnFpSByZWd1bGFyaXphxI1uw6kgbWV0w7NkeSwgbmFwcsOta2xhZCByaWRnZSByZWdyZXNpdS4KCiMjICBPZHN0csOhbmVuaWUgamVkbmVqIHoga29yZWxvdmFuw71jaCBwcmVtZW5uw71jaAoKQWsgc8O6IGR2ZSBwcmVtZW5uw6kgdGFrbWVyIHJvdm5ha8OpIGEgZWtvbm9taWNrw6EgdGXDs3JpYSBuZXZ5xb5hZHVqZSBvYmUsIG3DtMW+ZW1lIGplZG51IHogbmljaCB2eW5lY2hhxaUuCgpgYGB7ciByZWR1Y2VkLW1vZGVsfQptb2RlbF9yZWR1Y2VkIDwtIGxtKExpZmVfZXhwZWN0YW5jeSB+IEFsY29ob2xfY29uc3VtcHRpb24gKyBJbmNpZGVudHNfSElWLCBkYXRhID0gdWRhamUpCgpzdW1tYXJ5KG1vZGVsKQpzdW1tYXJ5KG1vZGVsX3JlZHVjZWQpCmBgYAo=