Econometrics in R - cvičenie 5

S využitím databázy Databaza.csv (so štruktúrou zhodnou s hárokom v_om7101mr_00_00_00_sk) s mesačnými údajmi za rok 2024 na úrovni obcí a okresov na Slovensku. Stĺpce 2024M01 až 2024M12 obsahujú mesačné hodnoty miery nezamestnanosti (%). V analýze pracujeme s rokom 2024, pričom pozorovaním je kombinácia obec × mesiac.

Pri práci používame knižnice:

library(readr)
library(dplyr)
library(tidyr)
library(stringr)
library(ggplot2)
library(lmtest)
library(sandwich)
library(car)
library(tseries)
rm(list = ls())

V hlavnom menu mám Session nastavené na Source File Location. Alternatívne možno použiť v Console príkaz, napr.:

setwd("/Cloud/project/tyzdne/tyzden5")

Údaje sú v súbore udaje/Databaza.csv (odporúčam mať súbor v podpriečinku udaje).

Úvod do problému a hypotézy

Chceme skúmať, ako sa mesačná miera nezamestnanosti v roku 2024 líši naprieč mesiacmi a okresmi/obcami.

• Budeme modelovať Mieru nezamestnanosti (v %) v závislosti od:
  • časového trendu v rámci roka (index mesiaca 1–12),
  • fixných efektov okresov (zachytia priestorové rozdiely),
  • a alternatívne aj sezónnych efektov mesiacov (faktor mesiaca).

Hypotézy: 1. Existujú významné rozdiely medzi okresmi (fixné efekty budú štatisticky významné). 2. V priebehu roka pozorujeme mierny trend/sezónnosť (v niektorých mesiacoch je nezamestnanosť vyššia/nižšia). 3. Diagnostika môže ukázať heteroskedasticitu (rozptyl rezíduí sa líši medzi obcami/okresmi), preto overíme aj robustné (HC) smerodajné chyby a transformáciu.

Načítanie, čistenie a príprava údajov

# ===== Načítanie, čistenie a príprava údajov (robustná verzia) =====

library(readr)
library(dplyr)
library(tidyr)
library(stringr)

# zvoľ cestu: ak existuje udaje/Databaza.csv, použi ju; inak hľadaj v koreňovom priečinku
FILE <- if (file.exists("udaje/Databaza.csv")) "udaje/Databaza.csv" else "Databaza.csv"

# 1) načítanie CSV; ak by bolo ; a desatinná čiarka, prepneme na read_csv2
raw <- suppressMessages(readr::read_csv(FILE, show_col_types = FALSE, locale = locale(encoding = "UTF-8")))
if (ncol(raw) == 1) {
  # pravdepodobne ; ako oddeľovač
  raw <- suppressMessages(readr::read_csv2(FILE, show_col_types = FALSE,
                                           locale = locale(encoding = "UTF-8", decimal_mark = ",", grouping_mark = " ")))
}

# 2) očistenie názvov stĺpcov (trim, na malé, nahradiť medzery podčiarkovníkom)
names(raw) <- gsub("\\s+", "_", tolower(trimws(names(raw))), perl = TRUE)

# 3) nájdi stĺpce s mesiacmi 2024: akceptuj 2024M01, 2024-01, 202401, 2024_01, 2024 01
month_cols <- names(raw)[grepl(
  "^2024(m)?(0[1-9]|1[0-2])$|^2024[-_ ]?(0[1-9]|1[0-2])$|^2024(0[1-9]|1[0-2])$",
  names(raw), perl = TRUE
)]

stopifnot(length(month_cols) > 0)

# 4) nájdi stĺpce okres/obec podľa názvov (fallback: prvé dva nemesačné stĺpce)
cand_okres <- names(raw)[grepl("^(okres|district)", names(raw), perl = TRUE)]
cand_obec  <- names(raw)[grepl("^(obec|municip|nazov_obce|nazov_obce|obec_nazov)", names(raw), perl = TRUE)]

non_month <- setdiff(names(raw), month_cols)

if (length(cand_okres) == 0) cand_okres <- non_month[1]
if (length(cand_obec)  == 0) cand_obec  <- non_month[2]

stopifnot(!is.na(cand_okres), !is.na(cand_obec))

# 5) premenuj na štandardné mená Okres/Obec
raw <- raw |>
  rename(Okres = all_of(cand_okres),
         Obec  = all_of(cand_obec))

# 6) wide -> long
df_long <- raw |>
  tidyr::pivot_longer(all_of(month_cols), names_to = "Mesacnik", values_to = "Miera") |>
  mutate(
    # Mesacnik ako "2024M01" / "2024-01" / "202401" -> vyťaž mesiac
    Year  = 2024L,
    Month = as.integer(str_extract(Mesacnik, "(0[1-9]|1[0-2])")),
    MonthIndex = Month
  ) |>
  select(Okres, Obec, Year, Month, MonthIndex, Miera)

# 7) typy a chýbajúce hodnoty
df_long <- df_long |>
  mutate(
    Okres = as.factor(Okres),
    Obec  = as.factor(Obec),
    Miera = suppressWarnings(as.numeric(Miera))
  )

# 8) imputácia chýbajúcich – medián obce; ak obec celý rok chýba, medián okresu; ak aj to chýba, celkový medián
overall_med <- median(df_long$Miera, na.rm = TRUE)
df_long <- df_long |>
  group_by(Okres) |>
  mutate(okres_med = median(Miera, na.rm = TRUE)) |>
  ungroup() |>
  group_by(Obec) |>
  mutate(obec_med = median(Miera, na.rm = TRUE)) |>
  ungroup() |>
  mutate(
    Miera = ifelse(is.na(Miera),
                   ifelse(!is.na(obec_med), obec_med,
                          ifelse(!is.na(okres_med), okres_med, overall_med)),
                   Miera)
  ) |>
  select(-okres_med, -obec_med)

# 9) rýchla kontrola
summary(df_long$Miera)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
    100     883    1716    4922    3264  112794 

Teraz si vizualizujeme rozloženie a extrémy:

# krabicové grafy pre maticu (vybrané premenné)
par(mfrow=c(1,3), mar=c(4,4,2,1))
boxplot(df_long$Miera, main="Miera nezamestnanosti (%)", xlab="", col="lightblue")
boxplot(df_long$MonthIndex, main="Index mesiaca (1–12)", xlab="", col="lightblue")

boxplot(tapply(df_long$Miera, df_long$Okres, median), main="Medián podľa okresu", col="lightblue")
mtext("Boxploty vybraných premenných", outer = TRUE, cex = 1.2, font = 2)

par(mfrow=c(1,1))

Lineárna regresia (panel 2024)

Zvolíme dva prirodzené modely:

• Model A (trend + okresy):
  Miera ~ 1 + MonthIndex + Okres
  — meria priemerný trend v roku 2024 s rozdielmi medzi okresmi.

• Model B (sezónnosť + okresy):
  Miera ~ 1 + factor(Month) + Okres
  — namiesto lineárneho trendu zachytí plnú sezónnosť mesiacov.

# robustne nájdi mesačné stĺpce (2024 + mesiac v rôznych formátoch)
month_cols <- names(raw)[grepl(
  "(^|[^0-9])2024\\D*(0[1-9]|1[0-2])($|\\D)", 
  names(raw), perl = TRUE
)]

if (length(month_cols) < 2) {
  message("Pozor: našiel som iba ", length(month_cols), " mesačných stĺpcov.")
  message("Hlavička dát:"); print(names(raw))
  stop("Potrebujem aspoň 2 mesačné stĺpce (napr. 2024M01, 2024M02, ...).")
}

library(stringr)

df_long <- raw |>
  tidyr::pivot_longer(all_of(month_cols), names_to = "Mesacnik", values_to = "Miera") |>
  mutate(
    Year  = 2024L,
    Month = as.integer(str_match(Mesacnik, "(0[1-9]|1[0-2])")[,1]),
    MonthIndex = Month
  ) |>
  select(Okres, Obec, Year, Month, MonthIndex, Miera)

df_long <- droplevels(df_long)

if (dplyr::n_distinct(df_long$Okres) < 2) {
  message("Pozor: Okres má len jednu úroveň. Model s Okres ako faktor nedáva zmysel.")
}

if (dplyr::n_distinct(df_long$Month) < 2) {
  message("Pozor: Month má len jednu úroveň – Model B (sezónnosť) sa preskočí.")
}

# Model A
modelA <- lm(Miera ~ 1 + MonthIndex + Okres, data = df_long)
summary(modelA)

Call:
lm(formula = Miera ~ 1 + MonthIndex + Okres, data = df_long)

Residuals:
   Min     1Q Median     3Q    Max 
-28712  -1740  -1034    135  82312 

Coefficients: (1 not defined because of singularities)
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                   47498       2599  18.276  < 2e-16 ***
MonthIndex                       NA         NA      NA       NA    
OkresOkres Bratislava II      -5405       3001  -1.801   0.0718 .  
OkresOkres Bratislava III    -21554       3001  -7.182 9.12e-13 ***
OkresOkres Bratislava IV     -29986       2807 -10.682  < 2e-16 ***
OkresOkres Bratislava V      -17016       2906  -5.856 5.39e-09 ***
OkresOkres Dunajská Streda   -45596       2618 -17.415  < 2e-16 ***
OkresOkres Galanta           -44850       2635 -17.022  < 2e-16 ***
OkresOkres Hlohovec          -44118       2757 -16.004  < 2e-16 ***
OkresOkres Malacky           -44427       2648 -16.774  < 2e-16 ***
OkresOkres Pezinok           -43380       2674 -16.221  < 2e-16 ***
OkresOkres Senec             -43841       2643 -16.585  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9003 on 2389 degrees of freedom
Multiple R-squared:  0.4226,    Adjusted R-squared:  0.4202 
F-statistic: 174.9 on 10 and 2389 DF,  p-value: < 2.2e-16

# Model B – iba ak je aspoň 2+ mesiacov
if (dplyr::n_distinct(df_long$Month) >= 2) {
  modelB <- lm(Miera ~ 1 + factor(Month) + Okres, data = df_long)
  summary(modelB)
} else {
  message("Model B (factor(Month)) preskakujem: v dátach je len jeden mesiac.")
}

names(raw) <- gsub("\\s+", "_", tolower(trimws(names(raw))), perl = TRUE)

Interpretácia:

• Koeficient MonthIndex (Model A) ukazuje, či má nezamestnanosť v roku 2024 rastúci/klesajúci trend.
• Faktor Month (Model B) umožní porovnať každý mesiac so základným (január), a tým zmerať sezónnosť.
• Koeficienty Okres zachytávajú trvalé rozdiely medzi okresmi (napr. štruktúra ekonomiky, veľkosť trhu práce a pod.).

Diagnostické grafy

par(mfrow=c(2,2))
plot(modelA)
par(mfrow=c(1,1))

Residuals vs Fitted

• Ak reziduá kmitajú okolo nuly bez kónusu, je rozptyl približne konštantný (homoskedasticita). Kónus –> heteroskedasticita.

Q–Q plot

• Odchýlky na koncoch znamenajú nenormalitu rezíduí (ťažšie chvosty).

Scale-Location

• Rovnomerná výška/šírka naznačuje konštantný rozptyl.

Residuals vs Leverage

• Sledujeme vplyvné pozorovania (Cookove vzdialenosti).

Formálne testy normality a odľahlostí

resA <- residuals(modelA)

# Jarque–Bera test normality
jb_A <- jarque.bera.test(resA)
jb_A

    Jarque Bera Test

data:  resA
X-squared = 139874, df = 2, p-value < 2.2e-16

# Outlier test (Bonferroni)
out_A <- outlierTest(modelA)
out_A

Ak normalita nevychádza (pri väčšom n je to bežné), môžeme skúsiť miernu transformáciu (napr. log(…)), ale pri percentách hrozí nula. Použijeme log(Miera + c), kde c je malé kladné číslo.

c_eps <- max(1e-3, 0.01*median(df_long$Miera, na.rm=TRUE))
df_long <- df_long |> mutate(Miera_log = log(Miera + c_eps))

modelA_log <- lm(Miera_log ~ 1 + MonthIndex + Okres, data = df_long)
summary(modelA_log)

Call:
lm(formula = Miera_log ~ 1 + MonthIndex + Okres, data = df_long)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.91782 -0.66110 -0.03617  0.52232  2.99933 

Coefficients: (1 not defined because of singularities)
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 10.7688     0.2766  38.940  < 2e-16 ***
MonthIndex                       NA         NA      NA       NA    
OkresOkres Bratislava II    -0.3385     0.3193  -1.060  0.28919    
OkresOkres Bratislava III   -0.8977     0.3193  -2.811  0.00498 ** 
OkresOkres Bratislava IV    -1.3931     0.2987  -4.664 3.28e-06 ***
OkresOkres Bratislava V     -1.8724     0.3092  -6.056 1.62e-09 ***
OkresOkres Dunajská Streda  -3.7277     0.2786 -13.380  < 2e-16 ***
OkresOkres Galanta          -3.3050     0.2804 -11.788  < 2e-16 ***
OkresOkres Hlohovec         -3.5681     0.2933 -12.164  < 2e-16 ***
OkresOkres Malacky          -3.1714     0.2818 -11.254  < 2e-16 ***
OkresOkres Pezinok          -3.0171     0.2846 -10.602  < 2e-16 ***
OkresOkres Senec            -3.0874     0.2813 -10.976  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.958 on 2389 degrees of freedom
Multiple R-squared:  0.3501,    Adjusted R-squared:  0.3474 
F-statistic: 128.7 on 10 and 2389 DF,  p-value: < 2.2e-16

par(mfrow=c(2,2)); plot(modelA_log); par(mfrow=c(1,1))

# normalita a odľahlé po transformácii
jarque.bera.test(residuals(modelA_log))

    Jarque Bera Test

data:  residuals(modelA_log)
X-squared = 146.11, df = 2, p-value < 2.2e-16

outlierTest(modelA_log)

No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:

Heteroskedasticita

Otestujeme prítomnosť heteroskedasticity Breusch–Paganovým testom pre Model A aj transformovaný model:

bptest(modelA)

    studentized Breusch-Pagan test

data:  modelA
BP = 1025, df = 10, p-value < 2.2e-16

bptest(modelA_log)

    studentized Breusch-Pagan test

data:  modelA_log
BP = 100.3, df = 10, p-value < 2.2e-16

Ak je heteroskedasticita prítomná, použijeme White-HC robustné smerodajné chyby:

coeftest(modelA, vcov = vcovHC(modelA))      # robustné SE pre model A

t test of coefficients:

                             Estimate Std. Error   t value  Pr(>|t|)    
(Intercept)                 47498.333     27.059 1755.3841 < 2.2e-16 ***
OkresOkres Bratislava II    -5404.806   4954.148   -1.0910   0.27540    
OkresOkres Bratislava III  -21554.028   2757.620   -7.8162 8.096e-15 ***
OkresOkres Bratislava IV   -29986.347   1585.497  -18.9129 < 2.2e-16 ***
OkresOkres Bratislava V    -17016.271   6987.021   -2.4354   0.01495 *  
OkresOkres Dunajská Streda -45596.346    117.882 -386.7974 < 2.2e-16 ***
OkresOkres Galanta         -44850.377    161.300 -278.0559 < 2.2e-16 ***
OkresOkres Hlohovec        -44118.031    636.330  -69.3320 < 2.2e-16 ***
OkresOkres Malacky         -44426.676    223.010 -199.2142 < 2.2e-16 ***
OkresOkres Pezinok         -43379.926    394.388 -109.9930 < 2.2e-16 ***
OkresOkres Senec           -43841.109    221.210 -198.1875 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

coeftest(modelA_log, vcov = vcovHC(modelA_log))

t test of coefficients:

                              Estimate  Std. Error    t value  Pr(>|t|)    
(Intercept)                10.76880958  0.00056941 18912.3506 < 2.2e-16 ***
OkresOkres Bratislava II   -0.33853454  0.10899306    -3.1060  0.001919 ** 
OkresOkres Bratislava III  -0.89765407  0.14595436    -6.1502 9.042e-10 ***
OkresOkres Bratislava IV   -1.39311045  0.11792566   -11.8135 < 2.2e-16 ***
OkresOkres Bratislava V    -1.87238053  0.23760590    -7.8802 4.922e-15 ***
OkresOkres Dunajská Streda -3.72771556  0.03097039  -120.3639 < 2.2e-16 ***
OkresOkres Galanta         -3.30503736  0.04134388   -79.9402 < 2.2e-16 ***
OkresOkres Hlohovec        -3.56808417  0.11902083   -29.9787 < 2.2e-16 ***
OkresOkres Malacky         -3.17144336  0.05190112   -61.1055 < 2.2e-16 ***
OkresOkres Pezinok         -3.01705929  0.07296094   -41.3517 < 2.2e-16 ***
OkresOkres Senec           -3.08741108  0.05861326   -52.6743 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Vizualizácia heteroskedasticity (príklad)

# Reziduá^2 vs. fitted (Model A)
df_long$fitA <- fitted(modelA)
df_long$res2A <- residuals(modelA)^2

ggplot(df_long, aes(x = fitA, y = res2A)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "loess", se = FALSE) +
  labs(x = "Fitted values", y = "Squared residuals",
       title = "Heteroskedasticita: reziduá^2 vs. fitted (Model A)") +
  theme_minimal()

Zhrnutie výsledkov

• Okresy: Fixné efekty okresov sú spravidla štatisticky významné – regionálne rozdiely v miere nezamestnanosti v roku 2024 existujú.
• Čas:
  • V Modeli A koeficient MonthIndex ukáže, či v roku 2024 existuje trend (rast/klesanie).
    
  • V Modeli B porovnanie factor(Month) odhalí sezónny profil.
• Diagnostika: Normalita rezíduí sa pri veľkom n často formálne zamieta; dôležitejšia je robustnosť odhadov. BP test prípadnú heteroskedasticitu odhalí a v týchto prípadoch používame HC robustné chyby.
• Interpretácia je teraz priamo viazaná na tvoju databázu 2024M01–2024M12 na úrovni obec/okres, nie na WHO life expectancy.

Poznámky k rozšíreniam

• Vieme pridať efekty obcí, náhodné efekty (balík lme4), mapové vizualizácie podľa okresov, alebo porovnať prvú vs. druhú polovicu roka 2024.