Econometrics in R - cvičenie 5

Základná diagnostika modelu: ceny bývania na Slovensku

Autor

Riano Vajas
(podľa predlohy V. Gazdu)

Publikované

9. novembra 2026

Pokračujem s databázou Eurostatu House price index (prc_hpi_a), index cien bývania, 2015 = 100.

# Spustite raz v konzole, iba ak tieto baliky chybaju.
install.packages(c("broom", "knitr", "kableExtra", "ggplot2", "dplyr",
                   "tidyr", "tseries", "car"))
library(broom)
library(knitr)
library(kableExtra)
library(ggplot2)
library(dplyr)
library(tidyr)

Príprava databázy

Rovnako ako v minulom cvičení načítam súbor udaje/ceny_byvania_eu.csv (ak chýba, stiahne sa z Eurostatu) a vytvorím tabuľku, kde jeden riadok je jeden rok.

cesta <- "udaje/ceny_byvania_eu.csv"
if (!file.exists(cesta)) {
  url <- paste0(
    "https://ec.europa.eu/eurostat/api/dissemination/sdmx/2.1/data/prc_hpi_a/",
    "A.TOTAL+DW_NEW+DW_EXST.I15_A_AVG.SK+CZ+HU+PL+AT+DE+EU27_2020/",
    "?format=SDMX-CSV&startPeriod=2010"
  )
  dir.create("udaje", showWarnings = FALSE)
  write.table(read.csv(url), cesta, sep = ";", dec = ".", row.names = FALSE)
}

udaje <- read.csv2(cesta, sep = ";", dec = ".",
                   na.strings = c("", "NA"), stringsAsFactors = FALSE)

udaje <- udaje[udaje$purchase == "TOTAL", c("geo", "TIME_PERIOD", "OBS_VALUE")]
udaje.roky <- as.data.frame(pivot_wider(udaje, names_from = geo,
                                        values_from = OBS_VALUE))
udaje.roky <- udaje.roky[order(udaje.roky$TIME_PERIOD), ]
names(udaje.roky)[names(udaje.roky) == "TIME_PERIOD"] <- "Rok"
names(udaje.roky)[names(udaje.roky) == "EU27_2020"] <- "EU"

udaje.roky <- na.omit(udaje.roky[, c("Rok", "SK", "CZ", "EU")])
rownames(udaje.roky) <- udaje.roky$Rok   # v grafoch sa zobrazia roky
head(udaje.roky)
      Rok     SK    CZ     EU
2010 2010  97.06  95.2 100.26
2011 2011  95.60  95.3 100.88
2012 2012  92.94  93.9  98.76
2013 2013  93.69  93.9  97.06
2014 2014  94.96  96.2  97.75
2015 2015 100.00 100.0 100.00

Uloženie workspace, čistenie workspace

  • Aby som nemusel v každom cvičení znova načítavať a upravovať databázu, uložím si celý obsah pamäte (workspace) do súboru:
save.image("pamat.Rdata")
  • Mazanie niektorých objektov z pamäte:
pokus1 <- 6
pokus2 <- "text"           # tieto dva objekty chcem vymazat

rm(list = c("pokus1", "pokus2"))

ls()                       # co zostalo v pamati
[1] "cesta"      "ciarka"     "pformat"    "udaje"      "udaje.roky"
  • Čistenie celej pamäte:
rm(list = ls())
ls()
character(0)
  • Opätovné načítanie uloženého workspace:
load("pamat.Rdata")
ls()
[1] "cesta"      "ciarka"     "pformat"    "udaje"      "udaje.roky"

Na začiatku ďalších cvičení teda stačí load("pamat.Rdata") a na konci dokumentu save.image("pamat.Rdata"). Nepotrebné objekty budem priebežne mazať príkazom rm().

Lineárna regresia

Odhadujem model z minulého cvičenia:

\[ SK_t = \beta_0 + \beta_1 CZ_t + \beta_2 EU_t + u_t \]

kde \(SK_t\), \(CZ_t\) a \(EU_t\) sú indexy cien bývania na Slovensku, v Česku a priemer EÚ v roku \(t\).

model <- lm(SK ~ CZ + EU, data = udaje.roky)
tidy(model) |>
  mutate(
    Significance = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01  ~ "**",
      p.value < 0.05  ~ "*",
      p.value < 0.10  ~ ".",
      TRUE ~ ""
    )
  ) |>
  kable(
    digits = 4,
    col.names = c("Variable", "Estimate", "Std. Error", "t value",
                  "p-value", "Significance"),
    caption = "Regresné koeficienty"
  ) |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Regresné koeficienty
Variable Estimate Std. Error t value p-value Significance
(Intercept) -42.7681 33.6774 -1.2699 0.2264
CZ 0.1638 0.2413 0.6788 0.5091
EU 1.2439 0.5690 2.1860 0.0477 *

Interpretácia výsledkov

Stĺpce tabuľky: Estimate je odhad koeficientu \(\hat\beta_j\), Std. Error jeho štandardná chyba (neistota odhadu), t value je testovacia štatistika \(t_j = \hat\beta_j / SE(\hat\beta_j)\) a p-value je p-hodnota t-testu.

t-test pre každý koeficient testuje

\[ H_0:\beta_j=0 \qquad H_1:\beta_j\neq 0. \]

Ak je p-hodnota menšia ako 0,05, zamietam \(H_0\) a koeficient považujem za štatisticky významný. Pri platnosti \(H_0\) má štatistika t-rozdelenie s \(n-k\) = 13 stupňami voľnosti.

Koeficient pri CZ: Ak je index v Česku vyšší o 1 bod (pri rovnakom indexe EÚ), index na Slovensku je v priemere vyšší o 0,164 bodu. Koeficient nie je štatisticky významný na 5 % hladine (p = 0,509).

Koeficient pri EU: Ak je index EÚ vyšší o 1 bod (pri rovnakom indexe v Česku), index na Slovensku je v priemere vyšší o 1,244 bodu. Koeficient je štatisticky významný na 5 % hladine (p = 0,048).

Absolútny člen (Intercept): hodnota indexu SK, ak by indexy CZ aj EÚ boli nulové. Má len technický význam (priesečník s osou).

Pozor: Všetky tri indexy rastú v čase, takže spolu silno súvisia (multikolinearita). Štandardné chyby sú preto väčšie a koeficienty môžu vyjsť nevýznamné alebo s nečakaným znamienkom, aj keď model ako celok funguje dobre.

Kvalita vyrovnania modelu

glance(model) |>
  select(-deviance) |>
  kable(
    digits = 4,
    col.names = c("R.squared", "Adj.R.squared", "Sigma", "F Stat.", "p.value",
                  "df", "logLik", "AIC", "BIC", "df.residual", "N"),
    caption = "Model summary"
  ) |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Model summary
R.squared Adj.R.squared Sigma F Stat. p.value df logLik AIC BIC df.residual N
0.9946 0.9938 2.9474 1207.462 0 2 -38.3368 84.6737 87.764 13 16

Koeficient determinácie \(R^2 = 1 - RSS/TSS\) = 0,995. Model vysvetľuje približne 99,5 % variability indexu cien na Slovensku.

Upravený \(\bar R^2\) = 0,994 zohľadňuje počet premenných. Je trochu nižší ako \(R^2\), lebo mám len 16 pozorovaní.

Sigma \(\hat\sigma = \sqrt{RSS/(n-k)}\) = 2,95: skutočný index sa od modelu typicky odchyľuje asi o 2,9 bodu indexu.

F-test overuje, či sú vysvetľujúce premenné spoločne významné:

\[ H_0:\beta_1=\beta_2=0 \qquad H_1: \text{aspoň jeden } \beta_j \neq 0 \]

\[ F=\frac{R^2/(k-1)}{(1-R^2)/(n-k)} = 1207,5, \qquad p < 0{,}001. \]

Keďže p < 0,05, zamietam nulovú hypotézu: model ako celok je štatisticky významný.

df = 2 (počet vysvetľujúcich premenných), df.residual = \(n-k\) = 13, N = 16 rokov.

AIC a BIC použijem neskôr pri porovnaní modelov (nižšia hodnota = lepší model). logLik samostatne neinterpretujem.

Overenie vzťahu medzi RSS (deviance) a sigmou:

deviance(model)
[1] 112.9353
sigma(model)
[1] 2.947428
sqrt(deviance(model) / df.residual(model))
[1] 2.947428

Grafická analýza kvality modelu

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

Diagnostické grafy regresného modelu
par(mfrow = c(1, 1))

V grafoch sú pozorovania označené rokmi. Čo sledujem v jednotlivých grafoch:

  • Residuals vs Fitted: rezíduá by mali kolísať okolo nuly bez tvaru. Zakrivená červená čiara naznačuje nelinearitu. Pri časovom rade si všímam aj to, či rezíduá susedných rokov nemajú rovnaké znamienko (to by naznačovalo autokoreláciu).
  • Normal Q-Q: body by mali ležať blízko priamky. Odchýlky na koncoch znamenajú ťažšie chvosty alebo odľahlé hodnoty.
  • Scale-Location: ak je červená čiara približne vodorovná, rozptyl rezíduí je konštantný (homoskedasticita). Lievik naznačuje heteroskedasticitu.
  • Residuals vs Leverage: body za prerušovanými čiarami Cookovej vzdialenosti (0,5 a 1) by výrazne ovplyvňovali odhad.

Konkrétne pre môj model: najvyššiu leverage má rok 2020 (h = 0,387; priemer je p/n = 0,188) a najväčšiu Cookovu vzdialenosť rok 2025 (D = 0,863). Žiadne pozorovanie neprekračuje hranicu 1, takže žiadny rok neovplyvňuje odhad neprimerane.

Test normality rezíduí

Rezíduá \(\mathbf e = (\mathbf I - \mathbf H)\mathbf u\) sú lineárnou kombináciou náhodných zložiek, preto normalitu \(\mathbf u\) testujem nepriamo cez rezíduá.

Jarqueov–Berov test

  • \(H_0\): rezíduá sú normálne rozdelené.
  • \(H_1\): rezíduá nie sú normálne rozdelené.

\[ JB = \frac{n}{6}\left(S^2+\frac{(K-3)^2}{4}\right) \overset{a}{\sim} \chi^2(2) \]

jb_test <- tseries::jarque.bera.test(residuals(model))

tidy(jb_test) |>
  select(statistic, p.value, method) |>
  kable(digits = 4, col.names = c("JB statistic", "p-value", "Test"),
        caption = "Jarque-Bera test normality rezíduí") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Jarque-Bera test normality rezíduí
JB statistic p-value Test
0.501 0.7784 Jarque Bera Test

p-hodnota (0,778) je väčšia ako 0,05, preto hypotézu o normalite rezíduí nezamietam. Treba však pamätať, že JB test je asymptotický a pri iba 16 pozorovaniach má malú silu, preto ho dopĺňam Q-Q grafom.

Test odľahlých pozorovaní – outlierTest()

Test využíva externe študentizované rezíduá \(t_i = e_i / (\hat\sigma_{(i)}\sqrt{1-h_{ii}})\) a Bonferroniho korekciu \(p_{Bonf} = \min(1, n \cdot p_i)\).

  • \(H_0\): pozorovanie nie je odľahlé.
  • \(H_1\): pozorovanie je odľahlé.
outlier_test <- car::outlierTest(model)

outlier_table <- data.frame(
  Observation = names(outlier_test$rstudent),
  rstudent = unname(outlier_test$rstudent),
  unadjusted.p = unname(outlier_test$p),
  bonf.p = unname(outlier_test$bonf.p)
)

outlier_table |>
  kable(digits = 4,
        col.names = c("Rok", "Studentized residual",
                      "Unadjusted p-value", "Bonferroni p-value"),
        caption = "Test odľahlých pozorovaní (Bonferroni)") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Test odľahlých pozorovaní (Bonferroni)
Rok Studentized residual Unadjusted p-value Bonferroni p-value
2025 2.7792 0.0167 0.2668

Najextrémnejšie je pozorovanie z roku 2025. Bonferroniho p-hodnota je 0,267, teda väčšia ako 0,05: v modeli nemám štatisticky významne odľahlé pozorovanie.

Úprava modelu

Podobne ako v predlohe skúsim logaritmickú transformáciu vysvetľujúcich premenných:

\[ SK_t = \beta_0 + \beta_1 \ln(CZ_t) + \beta_2 \ln(EU_t) + u_t \]

Interpretácia: ak sa index v Česku zvýši o 1 %, index na Slovensku sa zmení približne o \(\beta_1/100\) bodu.

model2 <- lm(SK ~ log(CZ) + log(EU), data = udaje.roky)

tidy(model2) |>
  kable(digits = 4,
        col.names = c("Variable", "Estimate", "Std. Error", "t value", "Pr(>|t|)"),
        caption = "Model s logaritmami") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Model s logaritmami
Variable Estimate Std. Error t value Pr(>&#124;t&#124;)
(Intercept) -832.9038 217.4374 -3.8305 0.0021
log(CZ) 0.2143 49.8285 0.0043 0.9966
log(EU) 201.6146 96.6306 2.0864 0.0572

Grafická diagnostika upraveného modelu

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

Diagnostické grafy upraveného modelu
par(mfrow = c(1, 1))

Test normality a odľahlých hodnôt upraveného modelu

jb_test2 <- tseries::jarque.bera.test(residuals(model2))
tidy(jb_test2) |>
  select(statistic, p.value, method) |>
  kable(digits = 4, col.names = c("JB statistic", "p-value", "Test")) |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
JB statistic p-value Test
4.463 0.1074 Jarque Bera Test 
car::outlierTest(model2)
     rstudent unadjusted p-value Bonferroni p
2025 5.105662         0.00025935    0.0041497

Malé úlohy

Iná modifikácia modelu

V minulom cvičení som videl, že rast cien sa v čase zrýchľoval (rezíduá trendu mali tvar písmena U). Skúsim preto kvadratický trend:

\[ SK_t = \beta_0 + \beta_1 Rok_t + \beta_2 Rok_t^2 + u_t \]

udaje.roky$t <- udaje.roky$Rok - 2010     # cas od roku 2010 (0, 1, 2, ...)
model3 <- lm(SK ~ t + I(t^2), data = udaje.roky)

tidy(model3) |>
  kable(digits = 4,
        col.names = c("Variable", "Estimate", "Std. Error", "t value", "Pr(>|t|)"),
        caption = "Model s kvadratickým trendom") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Model s kvadratickým trendom
Variable Estimate Std. Error t value Pr(>&#124;t&#124;)
(Intercept) 93.2335 3.2838 28.3922 0.000
t -0.7352 1.0158 -0.7237 0.482
I(t^2) 0.5481 0.0653 8.3913 0.000

Rok som posunul tak, aby začínal nulou (t = Rok - 2010). Výsledky sa tým nemenia, len sa čísla v R lepšie počítajú (rok 2025² je veľmi veľké číslo).

Porovnanie modelov pomocou AIC

Všetky tri modely majú rovnakú vysvetľovanú premennú (SK) aj rovnaké pozorovania, preto ich môžem porovnať pomocou AIC (a BIC).

porovnanie <- data.frame(
  Model = c("SK ~ CZ + EU", "SK ~ log(CZ) + log(EU)", "SK ~ t + t^2"),
  R2 = c(glance(model)$r.squared, glance(model2)$r.squared, glance(model3)$r.squared),
  Adj.R2 = c(glance(model)$adj.r.squared, glance(model2)$adj.r.squared,
             glance(model3)$adj.r.squared),
  AIC = c(AIC(model), AIC(model2), AIC(model3)),
  BIC = c(BIC(model), BIC(model2), BIC(model3)),
  JB.p = c(jb_test$p.value, jb_test2$p.value,
           tseries::jarque.bera.test(residuals(model3))$p.value)
)

porovnanie |>
  kable(digits = 3, caption = "Porovnanie modelov") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)
Porovnanie modelov
Model R2 Adj.R2 AIC BIC JB.p
SK ~ CZ + EU 0.995 0.994 84.674 87.764 0.778
SK ~ log(CZ) + log(EU) 0.990 0.988 95.436 98.527 0.107
SK ~ t + t^2 0.985 0.983 101.180 104.270 0.204

Podľa AIC je najlepší model SK ~ CZ + EU (AIC = 84,67). Pôvodný a logaritmický model sa líšia o 10,76 bodu AIC, teda logaritmická transformácia model nezlepšila.

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

Diagnostické grafy modelu s kvadratickým trendom
par(mfrow = c(1, 1))

Záver

Odhadol som viacnásobnú lineárnu regresiu indexu cien bývania na Slovensku a urobil základnú diagnostiku: t-testy, F-test, \(R^2\), grafickú analýzu rezíduí, Jarqueov–Berov test a test odľahlých hodnôt. Model vysvetľuje veľkú časť variability, ale treba si uvedomiť dve obmedzenia:

  1. Mám len 16 pozorovaní, takže testy majú malú silu.
  2. Ide o časové rady, ktoré spolu rastú. Vysoké \(R^2\) preto môže byť čiastočne spôsobené spoločným trendom. V ďalších cvičeniach by bolo vhodné otestovať aj autokoreláciu rezíduí (napr. Durbin-Watsonov test).

Na konci uložím workspace, aby som ho mohol nabudúce načítať príkazom load("pamat.Rdata").

save.image("pamat.Rdata")
rm(list = ls())