Špecifikácia modelu

library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
library(ggplot2)
library(dplyr)
library(tidyr)
library(MASS)   # Box–Cox

rm(list = ls())

Úvod

Pokračujem v predchádzajúcej úlohe, kde na mojej databáze WDI – Education, Health & Employment
skúmam, ako sú výdavky na zdravotníctvo ovplyvnené výdavkami na vzdelávanie.

Závislá premenná:

• health_exp – Current health expenditure (% of GDP).

Vysvetľujúce premenné:

• 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.

1. Príprava údajov

# načítanie CSV súboru z priečinka data
udaje_raw <- read.csv(
  "data/wdi_data.csv",
  dec = ".",
  sep = ",",
  header = TRUE,
  check.names = FALSE   # necháme pôvodné názvy stĺpcov
)

# vyberiem len rok 2015
udaje_2015_all <- udaje_raw[udaje_raw$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")

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

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()

Moje zhrnutie boxplotov:
Vidno, že všetky premenné majú pomerne sústredené rozdelenie, ale v každej z nich sa objaví niekoľko krajín s nezvyčajne nízkymi alebo vysokými hodnotami (bodky mimo boxu). To sú potenciálne odľahlé krajiny, ktoré sa môžu neskôr prejaviť v diagnostických grafoch.

2. Základná regresia

Model špecifikacia:

\[ \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 ~ educ_prim + educ_sec + educ_ter,
  data = udaje.2015
)

summary(model)

## 
## Call:
## lm(formula = health_exp ~ 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

3. Ramsey RESET test

resettest(model)

## 
##  RESET test
## 
## data:  model
## RESET = 0.58969, df1 = 2, df2 = 13, p-value = 0.5687

Ako tomu rozumiem:
Ak je p-hodnota malá (pod 0,05), RESET test naznačuje, že model môže byť zle špecifikovaný – buď chýba nejaká premenná, alebo by sa hodilo zaviesť nelineárny člen (napr. druhú mocninu alebo logaritmus).

4. Diagnostické grafy

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

par(mfrow = c(1, 1))

Moja interpretácia diagnostických grafov

• Residuals vs Fitted
  Reziduá sa síce držia okolo nuly, ale červená LOESS krivka je mierne zakrivená. V strede ide trochu nad nulou, na koncoch pod ňou. Vnímam to ako náznak nelineárneho vzťahu, ktorý jednoduchá lineárna špecifikácia úplne nezachytí. Zároveň vidno niekoľko bodov s väčšími reziduami – to sú pravdepodobne tie odľahlé krajiny z boxplotov.

• Normal Q-Q
  Väčšina bodov leží pomerne blízko priamky, ale na začiatku a na konci sa už body odchyľujú. To znamená, že v chvostoch rozdelenia rezíduí sú trochu ťažšie chvosty a normalita nie je ideálna, ale nie je to katastrofa.

• Scale-Location
  Červená čiara nie je úplne vodorovná – pri niektorých fitted hodnotách je rozptyl rezíduí o niečo väčší. Rozptyl teda nie je úplne konštantný, ale ide skôr o miernu heteroskedasticitu, nie o dramatický problém.

• Residuals vs Leverage
  Väčšina pozorovaní má malý pákový efekt, zopár bodov má väčší leverage a zároveň väčšie reziduá. Žiadny bod však extrémne neprekračuje kontúry Cookovej vzdialenosti, takže nevidím vyslovene „nebezpečne“ vplyvné pozorovania, skôr pár krajín, ktoré si zaslúžia pozornosť.

5. Component + Residual Plots

crPlots(model)

Čo z týchto grafov vyčítam

• C+R graf pre educ_prim
  Ružová hladká krivka je výrazne zakrivená – najmä v strede hodnôt educ_prim ide najprv prudko nahor a potom klesá. To mi hovorí, že vzťah medzi educ_prim a health_exp je nelineárny a lineárny člen sám o sebe nestačí. Práve tu by dávala zmysel transformácia – napríklad druhá mocnina alebo logaritmus educ_prim.

• C+R graf pre educ_sec
  Tu je nelinearita slabšia, ale krivka tiež nie je úplne rovná – mierne sa ohýba. Vzťah so sekundárnym vzdelávaním je teda tiež trochu nelineárny, ale menej dramaticky ako pri educ_prim.

• C+R graf pre educ_ter
  Krivka je najbližšie k priamke, takže v prípade terciárneho vzdelávania by jednoduchý lineárny člen mohol stačiť. Nevidím tu silný dôvod na špeciálnu transformáciu.

6. Box–Coxov transformačný test

boxcox(model)

Interpretácia Box–Cox grafu

Na vodorovnej osi je parameter \(\lambda\), na zvislej hodnota log-vierohodnosti. Krivka má maximum približne pri \(\lambda \approx 1\) a prerušovaná horizontálna čiara označuje 95 % interval spoľahlivosti.

Pre moje údaje:

• maximum log-vierohodnosti je blízko \(\lambda = 1\),
• interval spoľahlivosti obsahuje hodnotu 1.

To znamená, že netreba robiť drastickú transformáciu závislej premennej health_exp. Ak by som veľmi chcela experimentovať, mohla by som skúsiť napríklad miernu mocninovú transformáciu (napr. \(\lambda = 1{,}2\)), ale očakávam, že výsledky by boli veľmi podobné ako pri pôvodnej premennej.

7. Krátke záverečné zhrnutie

• Základný lineárny model vysvetľuje vzťah medzi výdavkami na vzdelávanie a výdavkami na zdravotníctvo, ale podľa RESET testu a diagnostických grafov nie je úplne ideálne špecifikovaný.
  
• C+R grafy naznačujú, že najväčší problém je nelinearita pri educ_prim (a čiastočne aj pri educ_sec).
  
• Box–Cox test neukazuje potrebu silnej transformácie závislej premennej, takže úpravy by som skôr robila na strane vysvetľujúcich premenných (logaritmus alebo kvadrát educ_prim).

Ak by som chcela model ďalej vylepšovať, ďalším krokom by bola práve nelineárna špecifikácia (napríklad pridať I(log(educ_prim)) alebo I(educ_prim^2)) a porovnať, či sa zlepší upravené \(R^2\) a správanie diagnostických grafov.