# Spustite raz v konzole, iba ak tieto balíky chýbajú.
install.packages(c("knitr", "rmarkdown", "broom", "kableExtra", "ggplot2", "tidyr"))Lineárna regresia
Ceny bývania na Slovensku v čase a v porovnaní s okolím
1 Predmet výskumu
Pokračujem s databázou z minulého zadania: index cien bývania Eurostatu (prc_hpi_a), kde 2015 = 100. Budem si klásť tieto otázky:
- Ako rýchlo rastú ceny bývania na Slovensku v čase?
- Grafické znázornenie tejto závislosti.
- Odhad regresnej priamky a jej interpretácia.
- Súvisia ceny na Slovensku s cenami v Česku a v priemere EÚ?
- Predikcia indexu na ďalšie roky.
Diagnostickým testom sa zatiaľ nevenujem. Cieľom je pochopiť základnú myšlienku regresie a prácu s funkciou
lm().
2 Príprava údajov
2.1 Balíky
library(broom)
library(knitr)
library(kableExtra)
library(ggplot2)
library(tidyr)2.2 Načítanie CSV
Používam súbor udaje/ceny_byvania_eu.csv, ktorý som si stiahol v minulom cvičení. Ak v priečinku chýba, stiahne sa znova z Eurostatu. Stĺpce sú oddelené bodkočiarkou a desatinné čísla majú bodku.
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)
# Ponecham iba index za vsetky nehnutelnosti (TOTAL)
# a krajiny dam do samostatnych stlpcov: jeden riadok = jeden rok.
udaje <- udaje[udaje$purchase == "TOTAL", c("geo", "TIME_PERIOD", "OBS_VALUE")]
udaje.roky <- pivot_wider(udaje, names_from = geo, values_from = OBS_VALUE)
udaje.roky <- as.data.frame(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 <- udaje.roky[, c("Rok", "SK", "CZ", "EU")]
udaje.roky |>
kable(format = "html", digits = 2, row.names = FALSE,
caption = "Index cien bývania (2015 = 100)") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Rok | SK | CZ | EU |
|---|---|---|---|
| 2010 | 97.06 | 95.2 | 100.26 |
| 2011 | 95.60 | 95.3 | 100.88 |
| 2012 | 92.94 | 93.9 | 98.76 |
| 2013 | 93.69 | 93.9 | 97.06 |
| 2014 | 94.96 | 96.2 | 97.75 |
| 2015 | 100.00 | 100.0 | 100.00 |
| 2016 | 106.67 | 107.2 | 104.31 |
| 2017 | 112.91 | 119.7 | 109.31 |
| 2018 | 121.32 | 130.0 | 114.81 |
| 2019 | 132.40 | 141.9 | 120.40 |
| 2020 | 145.02 | 153.9 | 127.15 |
| 2021 | 154.27 | 184.2 | 137.87 |
| 2022 | 175.41 | 215.3 | 148.86 |
| 2023 | 175.08 | 211.7 | 148.44 |
| 2024 | 181.67 | 222.3 | 153.42 |
| 2025 | 204.22 | 245.4 | 161.89 |
Jeden riadok predstavuje jeden rok. Na rozdiel od predlohy, kde boli krajiny v jednom roku (prierezové údaje), mám tu údaje zoradené v čase, teda časový rad (time series).
| Premenná | Význam |
|---|---|
| Rok | Kalendárny rok. |
| SK | Index cien bývania na Slovensku; vysvetľovaná premenná. |
| CZ | Index cien bývania v Česku. |
| EU | Index cien bývania, priemer EÚ (27 krajín). |
2.3 Chýbajúce hodnoty
chybajuce <- data.frame(
Premenna = names(udaje.roky),
Pocet_NA = colSums(is.na(udaje.roky))
)
chybajuce |>
kable(format = "html", col.names = c("Premenná", "Počet chýbajúcich hodnôt"),
row.names = FALSE, caption = "Chýbajúce hodnoty") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Premenná | Počet chýbajúcich hodnôt |
|---|---|
| Rok | 0 |
| SK | 0 |
| CZ | 0 |
| EU | 0 |
V predlohe sa chýbajúce hodnoty dopĺňali mediánom. Pri časovom rade by to skreslilo trend (doplnila by sa hodnota zo stredu obdobia), preto prípadné riadky s NA radšej vynechám.
udaje.roky <- na.omit(udaje.roky)
premenne <- c("SK", "CZ", "EU")
opis <- data.frame(
Premenna = premenne,
Minimum = sapply(udaje.roky[premenne], min),
Median = sapply(udaje.roky[premenne], median),
Priemer = sapply(udaje.roky[premenne], mean),
Maximum = sapply(udaje.roky[premenne], max)
)
opis |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Premenná", "Minimum", "Medián", "Priemer", "Maximum"),
caption = "Opisná štatistika indexov") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Premenná | Minimum | Medián | Priemer | Maximum |
|---|---|---|---|---|
| SK | 92.94 | 117.11 | 130.20 | 204.22 |
| CZ | 93.90 | 124.85 | 144.13 | 245.40 |
| EU | 97.06 | 112.06 | 120.07 | 161.89 |
3 Graf
Očakávam, že index cien bývania na Slovensku s časom rastie.
graf_cas <- ggplot(udaje.roky, aes(x = Rok, y = SK)) +
geom_point(color = "steelblue", size = 2.5) +
labs(x = "Rok", y = "Index cien bývania SK (2015 = 100)") +
theme_minimal(base_size = 13)
graf_casOtázky: Je rast rovnomerný? Nie celkom: približne do roku 2012 ceny mierne klesali, potom rástli stále rýchlejšie a v roku 2023 sa rast na chvíľu zastavil.
4 Jednoduchá lineárna regresia
4.1 Zápis modelu
\[ SK_t = \beta_0 + \beta_1 Rok_t + u_t. \]
- \(t\) označuje rok.
- \(\beta_0\) je konštanta (priesečník s osou Y).
- \(\beta_1\) je smernica: o koľko bodov sa v priemere zmení index za jeden rok.
- \(u_t\) je časť indexu, ktorú lineárny trend nevysvetlí.
4.2 Odhad pomocou lm()
model1 <- lm(SK ~ Rok, data = udaje.roky)
vysledky1 <- tidy(model1)
vysledky1[c("term", "estimate")] |>
kable(format = "html", digits = 3, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Jednoduchá regresia: index cien SK a rok") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | -14974.940 |
| Rok | 7.487 |
Odhadnutá rovnica:
\[ \widehat{SK}_t = -14974,9 + 7,487\,Rok_t. \]
Interpretácia smernice: Index cien bývania na Slovensku v rokoch 2010 – 2025 rástol v priemere o 7,49 bodu ročne. Keďže 2015 = 100, je to zhruba 7,5 % z úrovne roku 2015 každý rok.
Interpretácia konštanty: Konštanta je hodnota priamky pre rok 0, čo nemá ekonomický zmysel (je veľmi ďaleko od pozorovaných rokov). Chápem ju iba ako parameter polohy priamky.
Rok sám osebe ceny nespôsobuje. Za rastom stoja mzdy, úrokové sadzby, ponuka bytov a ďalšie faktory, ktoré sa s časom menia. Trend ich len spoločne zachytáva.
4.3 Regresná priamka v grafe
graf_cas +
geom_abline(intercept = coef(model1)["(Intercept)"],
slope = coef(model1)["Rok"],
color = "firebrick", linewidth = 1)5 Viacnásobná lineárna regresia
Ceny na Slovensku môžu súvisieť s vývojom v okolí. Ako vysvetľujúce premenné použijem index v Česku a priemer EÚ:
\[ SK_t = \beta_0 + \beta_1 CZ_t + \beta_2 EU_t + u_t. \]
model <- lm(SK ~ CZ + EU, data = udaje.roky)
vysledky <- tidy(model)
vysledky[c("term", "estimate")] |>
kable(format = "html", digits = 3, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Viacnásobná regresia: SK vysvetlené CZ a EÚ") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | -42.768 |
| CZ | 0.164 |
| EU | 1.244 |
\[ \widehat{SK}_t = -42,77 + 0,164\,CZ_t + 1,244\,EU_t. \]
Koeficient pri CZ interpretujem pri rovnakom indexe EÚ: ak je index v Česku vyšší o 1 bod, index na Slovensku je v priemere vyšší o 0,164 bodu.
Pozor: Všetky tri indexy rastú v čase, takže spolu silno súvisia. Koeficienty sa preto môžu výrazne meniť a niektorý môže mať aj nečakané znamienko. Model ukazuje súvislosť, nie to, že ceny v Česku spôsobujú ceny na Slovensku.
6 Stručný súhrn modelu pomocou glance()
suhrn <- glance(model)
suhrn[c("nobs", "r.squared", "sigma")] |>
kable(format = "html", digits = c(0, 3, 3), row.names = FALSE,
col.names = c("Počet rokov", "R²", "Odhad SD chyby (body indexu)"),
caption = "Základné charakteristiky viacnásobného modelu") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Počet rokov | R² | Odhad SD chyby (body indexu) |
|---|---|---|
| 16 | 0.995 | 2.947 |
Model zachytáva približne 99,5 % variability indexu na Slovensku. Pri časových radoch s rovnakým trendom býva \(R^2\) vysoké takmer vždy, preto to neznamená, že model je „dokonalý“.
7 Skutočná hodnota, predikcia a rezíduum
\[ e_t = y_t - \widehat y_t. \]
porovnanie <- data.frame(
Rok = udaje.roky$Rok,
Skutocnost = udaje.roky$SK,
Predikcia = fitted(model1),
Reziduum = residuals(model1)
)
porovnanie |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Rok", "Skutočnosť", "Predikcia", "Rezíduum"),
caption = "Skutočný index SK a hodnoty z lineárneho trendu") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Rok | Skutočnosť | Predikcia | Rezíduum |
|---|---|---|---|
| 2010 | 97.06 | 74.05 | 23.01 |
| 2011 | 95.60 | 81.54 | 14.06 |
| 2012 | 92.94 | 89.02 | 3.92 |
| 2013 | 93.69 | 96.51 | -2.82 |
| 2014 | 94.96 | 104.00 | -9.04 |
| 2015 | 100.00 | 111.48 | -11.48 |
| 2016 | 106.67 | 118.97 | -12.30 |
| 2017 | 112.91 | 126.46 | -13.55 |
| 2018 | 121.32 | 133.94 | -12.62 |
| 2019 | 132.40 | 141.43 | -9.03 |
| 2020 | 145.02 | 148.92 | -3.90 |
| 2021 | 154.27 | 156.41 | -2.14 |
| 2022 | 175.41 | 163.89 | 11.52 |
| 2023 | 175.08 | 171.38 | 3.70 |
| 2024 | 181.67 | 178.87 | 2.80 |
| 2025 | 204.22 | 186.35 | 17.87 |
Rezíduá nie sú náhodne rozhádzané: na začiatku a na konci obdobia sú skôr kladné, v strede záporné. Naznačuje to, že rast cien nie je úplne lineárny (zrýchľoval sa).
8 Predikcia pre zadanú hodnotu
Pomocou trendu (model1) vypočítam index pre roky 2026 – 2028:
nove_udaje <- data.frame(Rok = c(2026, 2027, 2028))
nove_udaje$Predikcia <- predict(model1, newdata = nove_udaje)
nove_udaje |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Rok", "Predikcia indexu SK"),
caption = "Predikcie z lineárneho trendu") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Rok | Predikcia indexu SK |
|---|---|
| 2026 | 193.84 |
| 2027 | 201.33 |
| 2028 | 208.82 |
Ide o bodové predikcie mimo pozorovaného obdobia (extrapolácia). Predpokladajú, že trend bude pokračovať rovnako, čo pri cenách nehnuteľností nemusí platiť.
9 Malé cvičenia
1. Ktorá premenná je vysvetľovaná a ktoré vysvetľujúce? Vysvetľovaná je index cien na Slovensku (SK). Vysvetľujúce sú Rok v jednoduchom modeli a CZ, EU vo viacnásobnom modeli.
2. Jednoduchý model SK ~ CZ s grafom a regresnou priamkou.
model2 <- lm(SK ~ CZ, data = udaje.roky)
tidy(model2)[c("term", "estimate")] |>
kable(format = "html", digits = 3, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Jednoduchá regresia: SK a CZ") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 30.70 |
| CZ | 0.69 |
ggplot(udaje.roky, aes(x = CZ, y = SK)) +
geom_point(color = "steelblue", size = 2.5) +
geom_abline(intercept = coef(model2)[1], slope = coef(model2)[2],
color = "firebrick", linewidth = 1) +
labs(x = "Index cien CZ", y = "Index cien SK") +
theme_minimal(base_size = 13)3. Prečo pri viacnásobnom modeli dopĺňame „pri rovnakých ostatných premenných“? Lebo koeficient pri CZ ukazuje súvislosť s indexom na Slovensku po zohľadnení indexu EÚ. V jednoduchom modeli SK ~ CZ by v koeficiente bol schovaný aj vplyv celoeurópskeho vývoja.
4. Predikcia pre rok 2026 a ručné overenie.
predict(model1, newdata = data.frame(Rok = 2026)) # pomocou predict() 1
193.8413
coef(model1)[1] + coef(model1)[2] * 2026 # rucne dosadenie(Intercept)
193.8413
5. Rok s kladným a so záporným rezíduom.
porovnanie[which.max(porovnanie$Reziduum), ] # najvacsie kladne reziduum Rok Skutocnost Predikcia Reziduum
1 2010 97.06 74.04831 23.01169
porovnanie[which.min(porovnanie$Reziduum), ] # najvacsie zaporne reziduum Rok Skutocnost Predikcia Reziduum
8 2017 112.91 126.4577 -13.54772
Kladné rezíduum znamená, že skutočný index bol v danom roku vyšší, ako predpovedá trend. Záporné znamená, že bol nižší.
6. Prečo z kladného koeficienta nemožno vyvodiť kauzalitu? Ceny v SK, CZ aj EÚ ovplyvňujú tie isté faktory (úrokové sadzby ECB, rast miezd, inflácia). Model zachytí, že sa hýbu spolu, ale neukáže, čo je príčinou.