S využitím databázy WDI – Education, Health & Employment (2011–2021)
(súbor data/wdi_data.csv v podpriečinku data).

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

library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
library(ggplot2)
library(patchwork)   # na zloženie grafov vedľa seba
library(dplyr)
library(tidyr)

rm(list = ls())

Príprava databázy, čistenie a úprava údajov

Rozhodla som sa modelovať výdavky na zdravotníctvo

health_exp – Current health expenditure (% of GDP)

v závislosti od vysvetľujúcich premenných:

educ_prim – výdavky na primárne vzdelávanie
educ_sec – výdavky na sekundárne vzdelávanie
educ_ter – výdavky na terciárne vzdelávanie

Načítanie dát a výber roku 2015

# načítanie CSV súboru z priečinka data
udaje <- read.csv(
  "data/wdi_data.csv",
  dec = ".",
  sep = ",",
  header = TRUE,
  check.names = FALSE   # aby sa zachovali originálne názvy stĺpcov
)

# vyberiem len rok 2015
udaje_2015_all <- udaje[udaje$Time == 2015, ]

# názvy stĺpcov, ktoré potrebujem (podľa hlavičky WDI súboru)
health_col <- "Current health expenditure (% of GDP) [SH.XPD.CHEX.GD.ZS]"
prim_col   <- "Current education expenditure, primary (% of total expenditure in primary public institutions) [SE.XPD.CPRM.ZS]"
sec_col    <- "Current education expenditure, secondary (% of total expenditure in secondary public institutions) [SE.XPD.CSEC.ZS]"
ter_col    <- "Current education expenditure, tertiary (% of total expenditure in tertiary public institutions) [SE.XPD.CTER.ZS]"

# vytiahnem len tieto štyri stĺpce
udaje.2015 <- udaje_2015_all[, c(health_col, prim_col, sec_col, ter_col)]

# premenujem ich na kratšie názvy
names(udaje.2015) <- c("health_exp",
                       "educ_prim",
                       "educ_sec",
                       "educ_ter")

# rýchly popis
summary(udaje.2015)

##    health_exp       educ_prim        educ_sec        educ_ter    
##  Min.   : 3.466   Min.   :85.55   Min.   :86.52   Min.   :87.06  
##  1st Qu.: 9.724   1st Qu.:87.69   1st Qu.:89.63   1st Qu.:89.13  
##  Median :10.340   Median :94.13   Median :93.58   Median :90.35  
##  Mean   :10.093   Mean   :91.86   Mean   :93.13   Mean   :90.97  
##  3rd Qu.:10.788   3rd Qu.:94.45   3rd Qu.:96.77   3rd Qu.:91.69  
##  Max.   :16.491   Max.   :97.51   Max.   :97.59   Max.   :96.40  
##                   NA's   :6       NA's   :6       NA's   :5

Imputácia chýbajúcich hodnôt

# mediány jednotlivých premenných (bez NA)
column_medians <- sapply(udaje.2015, median, na.rm = TRUE)
column_medians

## health_exp  educ_prim   educ_sec   educ_ter 
##   10.34000   94.12962   93.58111   90.35175

# nahradím NA mediánom danej premennej
udaje_imputed <- udaje.2015
for (col in names(udaje.2015)) {
  udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]
}
udaje.2015 <- udaje_imputed

# kontrola po imputácii
summary(udaje.2015)

##    health_exp       educ_prim        educ_sec        educ_ter    
##  Min.   : 3.466   Min.   :85.55   Min.   :86.52   Min.   :87.06  
##  1st Qu.: 9.724   1st Qu.:90.26   1st Qu.:92.08   1st Qu.:90.09  
##  Median :10.340   Median :94.13   Median :93.58   Median :90.35  
##  Mean   :10.093   Mean   :92.58   Mean   :93.27   Mean   :90.81  
##  3rd Qu.:10.788   3rd Qu.:94.26   3rd Qu.:95.46   3rd Qu.:91.05  
##  Max.   :16.491   Max.   :97.51   Max.   :97.59   Max.   :96.40

Boxploty premenných

Teraz sa pozriem na rozdelenie jednotlivých premenných pomocou boxplotov.

box_data <- udaje.2015 %>%
  mutate(id = dplyr::row_number()) %>%
  pivot_longer(
    cols = c(health_exp, educ_prim, educ_sec, educ_ter),
    names_to  = "variable",
    values_to = "value"
  )

ggplot(box_data, aes(x = variable, y = value)) +
  geom_boxplot(fill = "lightblue") +
  facet_wrap(~ variable, scales = "free_y") +
  labs(
    title = "Boxploty jednotlivých premenných",
    x = "Premenná",
    y = "Hodnota"
  ) +
  theme_minimal()

Boxploty jednotlivých premenných

Interpretácia boxplotov

health_exp – väčšina krajín míňa na zdravotníctvo okolo 10 % HDP. Na grafe sú však aj bodky nižšie a vyššie od boxu – tie predstavujú krajiny s nezvyčajne nízkymi alebo vysokými výdavkami na zdravotníctvo.
educ_prim – väčšina pozorovaní je sústredená okolo 90–95 % a len pár krajín má nižší alebo vyšší podiel výdavkov na primárne vzdelávanie. Variabilita je skôr malá, ale zopár odľahlých bodov existuje.
educ_sec – podobný obraz ako pri primárnom vzdelávaní, no vidno niekoľko krajín s výrazne nižšími výdavkami na sekundárne vzdelanie.
educ_ter – rozdelenie je relatívne úzke, ale opäť sa objavujú krajiny, ktoré dávajú na terciárne vzdelávanie podstatne menej alebo viac než väčšina.

Celkovo boxploty ukazujú, že údaje sú pomerne stabilné, ale existuje niekoľko odľahlých krajín, s ktorými treba pri regresii počítať.

Lineárna regresia – prvý model

Model odhadujem príkazom lm():

\[ \text{health\_exp}_i = \beta_0 + \beta_1 \text{educ\_prim}_i + \beta_2 \text{educ\_sec}_i + \beta_3 \text{educ\_ter}_i + u_i \]

model <- lm(health_exp ~ +1 + educ_prim + educ_sec + educ_ter,
            data = udaje.2015)

summary(model)

## 
## Call:
## lm(formula = health_exp ~ +1 + educ_prim + educ_sec + educ_ter, 
##     data = udaje.2015)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3302 -0.3082  0.1696  0.4648  6.4845 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  6.01846   27.38203   0.220    0.829
## educ_prim   -0.15982    0.40800  -0.392    0.701
## educ_sec     0.05971    0.46031   0.130    0.899
## educ_ter     0.14647    0.24278   0.603    0.555
## 
## Residual standard error: 2.576 on 15 degrees of freedom
## Multiple R-squared:  0.05594,    Adjusted R-squared:  -0.1329 
## F-statistic: 0.2963 on 3 and 15 DF,  p-value: 0.8275

Diagnostické grafy – prvý model

par(mfrow = c(2, 2))
plot(model)

Diagnostické grafy – prvý model

par(mfrow = c(1, 1))

Interpretácia diagnostických grafov – prvý model

Residuals vs Fitted – reziduá kolíšu okolo nuly, ale červená LOESS krivka je mierne zakrivená. To naznačuje miernu nelinearitu a prítomnosť niekoľkých odľahlých pozorovaní (napr. body s veľkými kladnými alebo zápornými reziduami).
Normal Q-Q – väčšina bodov leží pri priamke, no v dolnom aj hornom chvoste sú odchýlky. Rozdelenie rezíduí má ťažšie chvosty, takže normalita nie je úplne ideálna, ale nie je to katastrofa.
Scale-Location – krivka nie je úplne rovná, rozptyl rezíduí je pri niektorých fitted hodnotách väčší. To naznačuje miernu heteroskedasticitu.
Residuals vs Leverage – zopár bodov má vyšší leverage a zároveň väčšie reziduá, ale väčšinou neprekračujú vonkajšie čiary Cookovej vzdialenosti. Ide o potenciálne vplyvné krajiny, ale nie extrémne.

Testy normality a odľahlých hodnôt – prvý model

residuals1 <- residuals(model)

jb_test1      <- jarque.bera.test(residuals1)
outlier_test1 <- outlierTest(model)

list(
  jarque_bera_model1   = jb_test1,
  outlier_test_model1  = outlier_test1
)

## $jarque_bera_model1
## 
##  Jarque Bera Test
## 
## data:  residuals1
## X-squared = 10.358, df = 2, p-value = 0.005634
## 
## 
## $outlier_test_model1
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 95 3.440047          0.0039827     0.075671

Jarque–Bera test potvrdzuje, že reziduá nie sú dokonale normálne, a outlierTest označí niekoľko krajín ako štatisticky výrazné odľahlé pozorovania.

Druhý model – transformácia premenných

Podobne ako na cvičení upravím model:

logaritmujem výdavky na primárne vzdelávanie educ_prim,
vypustím educ_ter, ktorá môže byť kolineárna s ostatnými.

\[ \text{health\_exp}_i = \beta_0 + \beta_1 \log(\text{educ\_prim}_i) + \beta_2 \text{educ\_sec}_i + u_i \]

model2 <- lm(health_exp ~ +1 + I(log(educ_prim)) + educ_sec,
             data = udaje.2015)

summary(model2)

## 
## Call:
## lm(formula = health_exp ~ +1 + I(log(educ_prim)) + educ_sec, 
##     data = udaje.2015)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3663 -0.6028  0.1653  0.8434  6.4103 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)
## (Intercept)        76.73464  125.70854   0.610    0.550
## I(log(educ_prim)) -16.49006   35.86152  -0.460    0.652
## educ_sec            0.08591    0.44427   0.193    0.849
## 
## Residual standard error: 2.526 on 16 degrees of freedom
## Multiple R-squared:  0.03124,    Adjusted R-squared:  -0.08986 
## F-statistic: 0.258 on 2 and 16 DF,  p-value: 0.7758

Diagnostické grafy – druhý model

par(mfrow = c(2, 2))
plot(model2)

Diagnostické grafy – druhý model

par(mfrow = c(1, 1))

Interpretácia diagnostických grafov – druhý model

Residuals vs Fitted – LOESS krivka je oveľa rovnejšia a reziduá sú symetrickejšie okolo nuly. Nelinearita je slabšia než v prvom modeli.
Normal Q-Q – body sú bližšie k priamke, odchýlky v chvostoch sú menšie, takže reziduá sú bližšie k normálnemu rozdeleniu.
Scale-Location – červená krivka je pomerne plochá, rozptyl rezíduí je stabilnejší → problém heteroskedasticity sa zmiernil.
Residuals vs Leverage – pákové účinky sú podobné, ale reziduá pri týchto bodoch sú menšie, takže žiadne pozorovanie výrazne „neťahá“ regresiu.

Celkovo diagnostické grafy ukazujú, že druhý model sa správa lepšie než prvý.

Testy normality a odľahlých hodnôt – druhý model

residuals2 <- residuals(model2)

jb_test2      <- jarque.bera.test(residuals2)
outlier_test2 <- outlierTest(model2)

list(
  jarque_bera_model2   = jb_test2,
  outlier_test_model2  = outlier_test2
)

## $jarque_bera_model2
## 
##  Jarque Bera Test
## 
## data:  residuals2
## X-squared = 8.6256, df = 2, p-value = 0.0134
## 
## 
## $outlier_test_model2
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##     rstudent unadjusted p-value Bonferroni p
## 93 -3.412601          0.0038568     0.073278

Skúmanie heteroskedasticity (podľa profesora)

Teraz sa pozriem, či rozptyl rezíduí závisí od vysvetľujúcich premenných.

Štvorce rezíduí a vysvetľujúce premenné – prvý model

resid_sq1 <- residuals(model)^2

p1 <- ggplot(udaje.2015, aes(x = educ_prim, y = resid_sq1)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "educ_prim",
       y = "Štvorce rezíduí",
       title = "Squared residuals vs educ_prim") +
  theme_minimal()

p2 <- ggplot(udaje.2015, aes(x = educ_sec, y = resid_sq1)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "educ_sec",
       y = "Štvorce rezíduí",
       title = "Squared residuals vs educ_sec") +
  theme_minimal()

p1 + p2

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Skúmanie heteroskedasticity – prvý model

Interpretácia heteroskedasticity – prvý model

V grafe štvorce rezíduí vs. educ_prim má červená LOESS krivka výrazný oblúk – pri stredných hodnotách educ_prim sú štvorce rezíduí väčšie než pri nižších a vyšších hodnotách. Rozptyl chýb teda závisí od úrovne výdavkov na primárne vzdelávanie.
Podobne v grafe štvorce rezíduí vs. educ_sec sú najväčšie hodnoty rezíduí v strede rozsahu educ_sec.

Oba grafy naznačujú, že v prvom modeli je prítomná heteroskedasticita – variancia rezíduí nie je konštantná.

Štvorce rezíduí – druhý model (s log(educ_prim))

resid_sq2 <- residuals(model2)^2

p3 <- ggplot(udaje.2015, aes(x = log(educ_prim), y = resid_sq2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "log(educ_prim)",
       y = "Štvorce rezíduí",
       title = "Squared residuals vs log(educ_prim)") +
  theme_minimal()

p4 <- ggplot(udaje.2015, aes(x = educ_sec, y = resid_sq2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "educ_sec",
       y = "Štvorce rezíduí",
       title = "Squared residuals vs educ_sec (model2)") +
  theme_minimal()

p3 + p4

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Skúmanie heteroskedasticity – druhý model

Interpretácia heteroskedasticity – druhý model

V grafe štvorce rezíduí vs. log(educ_prim) je LOESS krivka podstatne rovnejšia než v prvom modeli – veľkosť rezíduí je podobná pri nízkych aj vysokých hodnotách logaritmu educ_prim.
V grafe štvorce rezíduí vs. educ_sec (model2) síce vidno mierny hrb, ale celkovo je krivka hladšia a výkyvy sú menšie ako pri prvom modeli.

To ukazuje, že po logaritmickej transformácii educ_prim je heteroskedasticita slabšia a druhý model lepšie spĺňa predpoklad konštantného rozptylu.

Breusch–Pagan test heteroskedasticity

bp1 <- bptest(model)
bp2 <- bptest(model2)

list(
  breusch_pagan_model1 = bp1,
  breusch_pagan_model2 = bp2
)

## $breusch_pagan_model1
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 1.4632, df = 3, p-value = 0.6908
## 
## 
## $breusch_pagan_model2
## 
##  studentized Breusch-Pagan test
## 
## data:  model2
## BP = 1.6778, df = 2, p-value = 0.4322

Ak má prvý model nízku p-hodnotu a druhý model vyššiu, podporuje to záver, že log-transformácia educ_prim pomohla problém heteroskedasticity znížiť.

White heteroskedasticity-consistent rozptyly (robustné štandardné chyby)

# robustné (White) štandardné chyby pre prvý model
coeftest(model,  vcov = vcovHC(model))

## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  6.01846   14.31319  0.4205   0.6801
## educ_prim   -0.15982    0.26712 -0.5983   0.5586
## educ_sec     0.05971    0.28826  0.2071   0.8387
## educ_ter     0.14647    0.16401  0.8931   0.3859

# robustné (White) štandardné chyby pre druhý model
coeftest(model2, vcov = vcovHC(model2))

## 
## t test of coefficients:
## 
##                     Estimate Std. Error t value Pr(>|t|)
## (Intercept)        76.734642  91.866855  0.8353   0.4159
## I(log(educ_prim)) -16.490062  26.815821 -0.6149   0.5472
## educ_sec            0.085911   0.321963  0.2668   0.7930

Zhrnutie a záver

V práci som použila údaje z databázy WDI pre rok 2015 a analyzovala som, ako rôzne typy výdavkov na vzdelávanie súvisia s výdavkami na zdravotníctvo:

v prvom modeli som zahrnula tri typy vzdelávacích výdavkov: educ_prim, educ_sec, educ_ter,
diagnostické grafy a testy naznačili miernu nelinearitu, nenormalitu rezíduí a prítomnosť heteroskedasticity,
v druhom modeli som logaritmovala educ_prim a vypustila educ_ter, čím som zmenšila vplyv extrémnych hodnôt a zjednodušila špecifikáciu modelu.

Boxploty ukázali, že údaje obsahujú niekoľko odľahlých krajín, ale žiadne úplne šialené hodnoty. Porovnanie heteroskedasticity (štvorcov rezíduí a Breusch–Pagan testu) aj diagnostické grafy ukazujú, že druhý model má:

stabilnejší rozptyl rezíduí,
reziduá bližšie k normálnemu rozdeleniu,
lepšie splnené predpoklady klasického lineárneho regresného modelu.

Pri interpretácii sa opieram o:

p-hodnoty koeficientov (štatistická významnosť vzťahov),
znamienka koeficientov (či vyššie výdavky na vzdelávanie súvisia s vyššími výdavkami na zdravotníctvo),
\(R^2\) a upravené \(R^2\) (ako dobre model vysvetľuje variabilitu výdavkov na zdravotníctvo).

Na základe týchto výsledkov považujem druhý model za vhodnejší opis vzťahu medzi výdavkami na vzdelávanie a výdavkami na zdravotníctvo v krajinách zahrnutých v databáze WDI pre rok 2015.

Econometrics in R - cvičenie 5 (moja databáza WDI)

Natália Soligova

Príprava databázy, čistenie a úprava údajov

Načítanie dát a výber roku 2015

Imputácia chýbajúcich hodnôt

Boxploty premenných

Interpretácia boxplotov

Lineárna regresia – prvý model

Diagnostické grafy – prvý model

Interpretácia diagnostických grafov – prvý model

Testy normality a odľahlých hodnôt – prvý model

Druhý model – transformácia premenných

Diagnostické grafy – druhý model

Interpretácia diagnostických grafov – druhý model

Testy normality a odľahlých hodnôt – druhý model

Skúmanie heteroskedasticity (podľa profesora)

Štvorce rezíduí a vysvetľujúce premenné – prvý model

Interpretácia heteroskedasticity – prvý model

Štvorce rezíduí – druhý model (s log(educ_prim))

Interpretácia heteroskedasticity – druhý model

Breusch–Pagan test heteroskedasticity

White heteroskedasticity-consistent rozptyly (robustné štandardné chyby)

Zhrnutie a záver