Econometrics in R - cvičenie 5

Úvod a cieľ

Tento dokument nadväzuje na logiku zadania 5 (vizualizácie → štatistiky → testovanie hypotéz → regresia → diagnostika), ale pracuje výlučne s lokálnymi dátami zo súboru Databaza.csv. Nevykonávame žiadne mapovanie názvov premenných na iné názvy – používame presne tie stĺpce, ktoré v súbore sú: Okres, Obec a mesačné stĺpce 2024M12 až 2024M01.

Cieľom je ilustrovať: 1. prípravu dát (wide → long) a rýchly prehľad, 2. vizualizácie rozdelení a vzťahov medzi mesiacmi, 3. základné štatistiky po mesiacoch, 4. testovanie hypotéz (porovnanie dvoch mesiacov a ANOVA naprieč mesiacmi), 5. lineárnu regresiu s decembrom 2024 ako závislou veličinou, 6. diagnostiku modelu a stručné interpretačné poznámky.

Balíčky

suppressPackageStartupMessages({
  library(dplyr)
  library(readr)
  library(tidyr)
  library(ggplot2)
  library(knitr)
  library(kableExtra)
  library(broom)
  library(stringr)
  library(lmtest)
  library(sandwich)
  library(tseries)
  library(car)
})

Dáta a kontrola štruktúry

Očakávame CSV s bodkočiarkovým oddeľovačom a hlavičkou: Okres, Obec, 2024M12, …, 2024M01.

udaje_raw <- readr::read_delim("Databaza.csv", delim = ";", show_col_types = FALSE)

stopifnot(all(c("Okres", "Obec") %in% names(udaje_raw)))

mesiac_cols <- grep("^2024M\\d{2}$", names(udaje_raw), value = TRUE)
if (length(mesiac_cols) == 0) stop("Nenašli sa žiadne stĺpce typu 2024Mxx.")
mesiac_cols <- sort(mesiac_cols, decreasing = TRUE)  # 2024M12 ... 2024M01
mesiac_cols

 [1] "2024M12" "2024M11" "2024M10" "2024M09" "2024M08" "2024M07" "2024M06"
 [8] "2024M05" "2024M04" "2024M03" "2024M02" "2024M01"

Konverzia do long formátu

Long formát uľahčí agregácie a grafy podľa mesiaca.

udaje_long <- udaje_raw |>
  tidyr::pivot_longer(cols = all_of(mesiac_cols),
                      names_to = "YM",
                      values_to = "Hodnota") |>
  mutate(
    Rok     = as.integer(stringr::str_sub(YM, 1, 4)),
    Mesiac  = as.integer(stringr::str_sub(YM, 6, 7)),
    Hodnota = suppressWarnings(as.numeric(Hodnota))
  )

dplyr::glimpse(udaje_long)

Rows: 2,400
Columns: 6
$ Okres   <chr> "Okres Bratislava I", "Okres Bratislava I", "Okres Brati…
$ Obec    <chr> "Bratislava - mestská časť Staré Mesto", "Bratislava - m…
$ YM      <chr> "2024M12", "2024M11", "2024M10", "2024M09", "2024M08", "…
$ Hodnota <dbl> 47634, 47632, 47565, 47564, 47525, 47488, 47503, 47484, …
$ Rok     <int> 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2024, 20…
$ Mesiac  <int> 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 12, 11, 10, 9, 8,…

Vizualizácie

Rozdelenie hodnôt podľa mesiacov (boxplot)

ggplot(udaje_long, aes(x = factor(Mesiac), y = Hodnota)) +
  geom_boxplot(fill = "lightblue", colour = "darkblue") +
  theme_minimal() +
  labs(title = "Rozdelenie hodnôt podľa mesiacov (2024)",
       x = "Mesiac (1–12)", y = "Hodnota")

Vzťah medzi mesiacmi (scatter 2024M11 vs 2024M12)

if (all(c("2024M11","2024M12") %in% names(udaje_raw))) {
  udaje_scatter <- udaje_raw |>
    dplyr::select(Okres, Obec, `2024M11`, `2024M12`) |>
    mutate(`2024M11` = as.numeric(`2024M11`),
           `2024M12` = as.numeric(`2024M12`))

  ggplot(udaje_scatter, aes(x = `2024M11`, y = `2024M12`, colour = Okres)) +
    geom_point(alpha = 0.85) +
    theme_minimal() +
    labs(title = "Vzťah 2024M11 vs 2024M12",
         x = "2024M11", y = "2024M12", colour = "Okres")
} else {
  plot.new(); title("Na scatter chýbajú 2024M11 a/alebo 2024M12")
}

Základné štatistiky

Súhrnné štatistiky po mesiacoch

stat_tbl <- udaje_long |>
  group_by(Mesiac) |>
  summarise(
    n      = n(),
    mean   = mean(Hodnota, na.rm = TRUE),
    sd     = sd(Hodnota, na.rm = TRUE),
    min    = min(Hodnota, na.rm = TRUE),
    q25    = quantile(Hodnota, 0.25, na.rm = TRUE),
    median = median(Hodnota, na.rm = TRUE),
    q75    = quantile(Hodnota, 0.75, na.rm = TRUE),
    max    = max(Hodnota, na.rm = TRUE),
    .groups = "drop"
  )

knitr::kable(stat_tbl, digits = 2, caption = "Základné štatistiky podľa mesiacov (2024)")

Základné štatistiky podľa mesiacov (2024)

Mesiac	n	mean	sd	min	q25	median	q75	max
1	200	4912.38	11845.18	102	880.50	1701.5	3250.50	112794
2	200	4913.83	11845.32	103	878.75	1702.5	3260.50	112740
3	200	4915.00	11844.50	103	881.75	1703.5	3257.00	112686
4	200	4916.87	11844.14	103	879.75	1707.5	3264.00	112661
5	200	4919.06	11846.72	103	882.00	1707.0	3266.00	112584
6	200	4921.61	11848.82	103	886.50	1710.0	3261.50	112553
7	200	4923.44	11853.09	101	887.75	1713.0	3267.50	112528
8	200	4925.03	11853.64	101	890.00	1717.0	3269.25	112469
9	200	4926.24	11854.74	100	889.50	1719.0	3268.25	112449
10	200	4927.26	11854.32	100	886.75	1717.5	3273.75	112436
11	200	4929.35	11856.14	101	892.00	1719.0	3276.50	112400
12	200	4930.90	11862.18	101	891.50	1720.0	3279.00	112447

stat_tbl |>
  kable(digits = 2, caption = "Základné štatistiky podľa mesiacov (2024)") |>
  kable_styling(full_width = FALSE, bootstrap_options = c("striped","hover","condensed")) |>
  column_spec(1, bold = TRUE) |>
  row_spec(0, bold = TRUE, background = "#f2f2f2") |>
  add_header_above(c(" " = 2, "Štatistiky hodnôt" = 7))

Základné štatistiky podľa mesiacov (2024)

		Štatistiky hodnôt
Mesiac	n	mean	sd	min	q25	median	q75	max
1	200	4912.38	11845.18	102	880.50	1701.5	3250.50	112794
2	200	4913.83	11845.32	103	878.75	1702.5	3260.50	112740
3	200	4915.00	11844.50	103	881.75	1703.5	3257.00	112686
4	200	4916.87	11844.14	103	879.75	1707.5	3264.00	112661
5	200	4919.06	11846.72	103	882.00	1707.0	3266.00	112584
6	200	4921.61	11848.82	103	886.50	1710.0	3261.50	112553
7	200	4923.44	11853.09	101	887.75	1713.0	3267.50	112528
8	200	4925.03	11853.64	101	890.00	1717.0	3269.25	112469
9	200	4926.24	11854.74	100	889.50	1719.0	3268.25	112449
10	200	4927.26	11854.32	100	886.75	1717.5	3273.75	112436
11	200	4929.35	11856.14	101	892.00	1719.0	3276.50	112400
12	200	4930.90	11862.18	101	891.50	1720.0	3279.00	112447

Testovanie hypotéz

t-test: porovnanie priemerov (2024M12 vs 2024M06)

if (all(c("2024M12","2024M06") %in% names(udaje_raw))) {
  t12 <- udaje_long$Hodnota[udaje_long$YM == "2024M12"]
  t06 <- udaje_long$Hodnota[udaje_long$YM == "2024M06"]
  t.test(t12, t06)
} else {
  "Na t-test chýbajú stĺpce 2024M12 a/alebo 2024M06."
}

    Welch Two Sample t-test

data:  t12 and t06
t = 0.007836, df = 398, p-value = 0.9938
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2321.433  2340.013
sample estimates:
mean of x mean of y 
  4930.90   4921.61 

ANOVA: rozdiely medzi mesiacmi

anova.result <- aov(Hodnota ~ factor(Mesiac), data = udaje_long)
summary(anova.result)

                 Df    Sum Sq   Mean Sq F value Pr(>F)
factor(Mesiac)   11 8.663e+04      7875       0      1
Residuals      2388 3.354e+11 140439858               

Lineárna regresia (štýl zadania 5)

Budeme modelovať december 2024 (2024M12) ako funkciu vybraných mesiacov. Zvolíme 2024M11, 2024M06 a 2024M01 (iba ak existujú).

target_y <- "2024M12"
x_candidates <- c("2024M11","2024M06","2024M01")
x_cols <- intersect(x_candidates, names(udaje_raw))

udaje_wide <- udaje_raw |>
  dplyr::select(Okres, Obec, any_of(c(target_y, x_cols))) |>
  mutate(across(-c(Okres, Obec), ~ suppressWarnings(as.numeric(.))))

if (target_y %in% names(udaje_wide) && length(x_cols) > 0) {
  form <- as.formula(paste0("`", target_y, "` ~ ", paste(sprintf("`%s`", x_cols), collapse = " + ")))
  model <- lm(form, data = udaje_wide |> tidyr::drop_na(all_of(c(target_y, x_cols))))
  summary(model)
} else {
  model <- NULL
  "Regresiu nie je možné odhadnúť (chýba závislá premenná alebo vysvetľujúce stĺpce)."
}

Call:
lm(formula = form, data = tidyr::drop_na(udaje_wide, all_of(c(target_y, 
    x_cols))))

Residuals:
    Min      1Q  Median      3Q     Max 
-45.409  -2.011   0.838   3.181  31.438 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.70320    0.62947  -2.706  0.00742 ** 
`2024M11`    1.11619    0.02570  43.431  < 2e-16 ***
`2024M06`   -0.08987    0.04505  -1.995  0.04745 *  
`2024M01`   -0.02589    0.02126  -1.218  0.22470    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.133 on 196 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 1.411e+08 on 3 and 196 DF,  p-value: < 2.2e-16

Diagnostika modelu (grafy a testy)

if (!is.null(model)) {
  par(mfrow = c(2, 2))
  plot(model)
  par(mfrow = c(1, 1))

  # Normalita rezíduí: Jarque–Bera
  jb <- tseries::jarque.bera.test(residuals(model))
  jb

  # Potenciálne odľahlé pozorovania (Bonferroni)
  ot <- car::outlierTest(model)
  ot
} else {
  "Diagnostika nie je dostupná (model sa neodhadol)."
}

Koeficienty a intervaly spoľahlivosti

if (!is.null(model)) {
  coef_tbl <- broom::tidy(model, conf.int = TRUE) |>
    mutate(
      term = dplyr::recode(term, "(Intercept)" = "Intercept"),
      stars = dplyr::case_when(
        p.value < 0.001 ~ "***",
        p.value < 0.01  ~ "**",
        p.value < 0.05  ~ "*",
        p.value < 0.1   ~ "·",
        TRUE            ~ ""
      )
    ) |>
    transmute(
      Parameter    = term,
      Estimate     = estimate,
      `Std. Error` = std.error,
      `t value`    = statistic,
      `p value`    = p.value,
      `95% CI`     = stringr::str_c("[", round(conf.low, 3), ", ", round(conf.high, 3), "]"),
      Sig          = stars
    )

  coef_tbl |>
    kable(digits = 3, caption = paste0("OLS koeficienty: ", target_y, " ~ ", paste(x_cols, collapse = " + "))) |>
    kable_styling(full_width = FALSE, bootstrap_options = c("striped","hover","condensed")) |>
    column_spec(1, bold = TRUE) |>
    row_spec(0, bold = TRUE, background = "#f2f2f2") |>
    footnote(general = "Signif. codes: *** p<0.001, ** p<0.01, * p<0.05, · p<0.1.", threeparttable = TRUE)
} else {
  cat("Koeficienty nie sú k dispozícii (model sa neodhadol).")
}

OLS koeficienty: 2024M12 ~ 2024M11 + 2024M06 + 2024M01

Parameter	Estimate	Std. Error	t value	p value	95% CI	Sig
Intercept	-1.703	0.629	-2.706	0.007	[-2.945, -0.462]	**
`2024M11`	1.116	0.026	43.431	0.000	[1.066, 1.167]	***
`2024M06`	-0.090	0.045	-1.995	0.047	[-0.179, -0.001]	*
`2024M01`	-0.026	0.021	-1.218	0.225	[-0.068, 0.016]
Note:
Signif. codes: *** p<0.001, ** p<0.01, * p<0.05, · p<0.1.

Kvalita prispôsobenia modelu

if (!is.null(model)) {
  fit_tbl <- broom::glance(model) |>
    transmute(
      `R-squared`      = r.squared,
      `Adj. R-squared` = adj.r.squared,
      `F-statistic`    = statistic,
      `F p-value`      = p.value,
      `AIC`            = AIC,
      `BIC`            = BIC,
      `Num. obs.`      = nobs
    )

  fit_tbl |>
    kable(digits = 3, caption = "Model Fit Statistics") |>
    kable_styling(full_width = FALSE, bootstrap_options = c("condensed"))
} else {
  cat("Fit štatistiky nie sú k dispozícii (model sa neodhadol).")
}

Model Fit Statistics

R-squared	Adj. R-squared	F-statistic	F p-value	AIC	BIC	Num. obs.
1	1	141101301	0	1411.92	1428.411	200

Stručná interpretácia

• Boxploty ukazujú rozptyl hodnôt v jednotlivých mesiacoch; extrémy môžu indikovať špecifické obce/okresy.
• Scatter 2024M11 vs 2024M12 naznačí, či má november podobnú úroveň ako december (silná/ slabá lineárna väzba).
• t-test odpovie, či sa priemery dvoch mesiacov (napr. december vs. jún) štatisticky líšia.
• ANOVA testuje rozdiely naprieč všetkými mesiacmi.
• Regresia (december ~ iné mesiace) kvantifikuje, do akej miery vedia iné mesiace vysvetliť decembrovú hodnotu.
• Diagnostika pomáha overiť predpoklady OLS a identifikovať vplyvné/odľahlé pozorovania.