# Spustite raz v konzole, iba ak balíky chýbajú:
# install.packages(c("knitr", "broom", "kableExtra", "ggplot2"))
library(broom)
library(knitr)
library(kableExtra)
library(ggplot2)Lineárna regresia: úroková sadzba a vznik nových podnikov na Slovensku
1. Predmet výskumu
V tomto cvičení odhadneme a interpretujeme lineárnu regresiu na vlastných údajoch z Eurostatu. Postupujeme podľa predlohy, ale s jedným dôležitým rozdielom: naše údaje sú časový rad (jeden riadok je jeden štvrťrok), nie prierezové údaje.
Otázky:
- Súvisí vyššia 3-mesačná úroková sadzba s nižším počtom registrácií nových podnikov?
- Grafické znázornenie závislosti.
- Odhad regresnej priamky a jej interpretácia.
- Pridanie ďalších vysvetľujúcich premenných.
- Predikcia.
Pracovná hypotéza: vyššie úroky zvyšujú cenu financovania, a preto znižujú počet nových firiem. V tomto cvičení sa nevenujeme diagnostickým testom.
2. Príprava údajov
2.1 Balíky
2.2 Načítanie CSV
Súbory načítame z podpriečinka udaje, alebo ak tam nie sú, z priečinka vedľa dokumentu.
najdi <- function(subor) {
cesty <- c(file.path("udaje", subor), subor)
existujuce <- cesty[file.exists(cesty)]
if (length(existujuce) == 0) stop("Súbor ", subor, " nie je v priečinku udaje ani vedľa QMD.")
existujuce[1]
}
reg <- read.csv(najdi("sts_rb_q$defaultview_linear_2_0.csv"),
na.strings = c("", "NA"), check.names = FALSE,
stringsAsFactors = FALSE)
irt <- read.csv(najdi("irt_st_q$defaultview_linear_2_0.csv"),
na.strings = c("", "NA"), check.names = FALSE,
stringsAsFactors = FALSE)2.3 Výber premenných a spojenie
Z obidvoch súborov vyberieme len Slovensko (registrácie a bankroty) a eurozónu (3-mesačná sadzba), potom ich spojíme podľa štvrťroka.
registracie <- reg[reg$geo == "SK" & reg$indic_bt == "REG" & reg$unit == "I21",
c("TIME_PERIOD", "OBS_VALUE")]
names(registracie) <- c("obdobie", "registracie")
bankroty <- reg[reg$geo == "SK" & reg$indic_bt == "BKRT" & reg$unit == "I21",
c("TIME_PERIOD", "OBS_VALUE")]
names(bankroty) <- c("obdobie", "bankroty")
urok <- irt[irt$geo == "EA" & irt$int_rt == "IRT_M3",
c("TIME_PERIOD", "OBS_VALUE")]
names(urok) <- c("obdobie", "urok")
udaje <- merge(merge(registracie, urok, by = "obdobie"), bankroty, by = "obdobie")
udaje <- udaje[order(udaje$obdobie), ]
# časový index a oneskorená úroková sadzba (o 1 štvrťrok)
udaje$t <- seq_len(nrow(udaje))
udaje$urok_lag1 <- c(NA, head(udaje$urok, -1))
# dátum pre grafy
rok <- as.integer(substr(udaje$obdobie, 1, 4))
stvrtrok <- as.integer(substr(udaje$obdobie, 7, 7))
udaje$datum <- as.Date(sprintf("%d-%02d-01", rok, (stvrtrok - 1) * 3 + 1))
head(udaje) |>
kable(format = "html", digits = 2, row.names = FALSE,
caption = "Prvých šesť štvrťrokov vo vybraných údajoch") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| obdobie | registracie | urok | bankroty | t | urok_lag1 | datum |
|---|---|---|---|---|---|---|
| 2015-Q1 | 65.3 | 0.05 | 34.8 | 1 | NA | 2015-01-01 |
| 2015-Q2 | 71.7 | -0.01 | 24.5 | 2 | 0.05 | 2015-04-01 |
| 2015-Q3 | 73.2 | -0.03 | 21.8 | 3 | -0.01 | 2015-07-01 |
| 2015-Q4 | 69.1 | -0.09 | 32.5 | 4 | -0.03 | 2015-10-01 |
| 2016-Q1 | 72.3 | -0.19 | 21.4 | 5 | -0.09 | 2016-01-01 |
| 2016-Q2 | 76.6 | -0.26 | 26.5 | 6 | -0.19 | 2016-04-01 |
Jeden riadok predstavuje jeden štvrťrok. Takejto organizácii údajov hovoríme časový rad (time series). Na rozdiel od prierezových údajov v predlohe tu záleží na poradí pozorovaní.
| Premenná | Význam |
|---|---|
obdobie |
Štvrťrok (napr. 2015-Q1). |
registracie |
Index registrácií nových podnikov na Slovensku (2021 = 100, sezónne a kalendárne očistený); vysvetľovaná premenná. |
urok |
3-mesačná medzibanková sadzba eurozóny v % p. a.; hlavná vysvetľujúca premenná. |
bankroty |
Index vyhlásených bankrotov na Slovensku (2021 = 100). |
t |
Časový index 1, 2, 3, … (jednoduchý trend). |
urok_lag1 |
Úroková sadzba oneskorená o jeden štvrťrok. |
2.4 Chýbajúce hodnoty
chybajuce <- data.frame(
Premenna = c("registracie", "urok", "bankroty", "urok_lag1"),
Pocet_NA = colSums(is.na(udaje[c("registracie", "urok", "bankroty", "urok_lag1")]))
)
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 |
|---|---|
| registracie | 0 |
| urok | 0 |
| bankroty | 0 |
| urok_lag1 | 1 |
Základné premenné nemajú chýbajúce hodnoty. Chýbajúca hodnota je len v prvom štvrťroku pre urok_lag1, a to prirodzene (pred rokom 2015-Q1 nemáme v údajoch sadzbu). Nič nedopĺňame, tento riadok sa pri modeli s oneskorením jednoducho vynechá.
Poznámka ku kvalite dát: všetky slovenské hodnoty registrácií majú v súbore príznak i, teda ide o hodnoty doplnené Eurostatom alebo inou prijímajúcou inštitúciou. Výsledky preto treba interpretovať opatrne.
table(reg$OBS_FLAG[reg$geo == "SK" & reg$indic_bt == "REG" & reg$unit == "I21"],
useNA = "ifany")
i
47 2.5 Opisná štatistika
premenne <- c("registracie", "urok", "bankroty")
opis <- data.frame(
Premenna = premenne,
Minimum = sapply(udaje[premenne], min),
Median = sapply(udaje[premenne], median),
Priemer = sapply(udaje[premenne], mean),
Maximum = sapply(udaje[premenne], max)
)
opis |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Premenná", "Minimum", "Medián", "Priemer", "Maximum"),
caption = "Opisná štatistika (47 štvrťrokov)") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Premenná | Minimum | Medián | Priemer | Maximum |
|---|---|---|---|---|
| registracie | 51.60 | 96.9 | 92.53 | 125.20 |
| urok | -0.57 | -0.3 | 0.76 | 3.96 |
| bankroty | 18.00 | 84.7 | 88.69 | 192.40 |
Úroková sadzba bola dlho okolo nuly alebo záporná (obdobie nízkych sadzieb po roku 2015) a prudko vzrástla po roku 2022.
3. Grafy
3.1 Vývoj v čase
ggplot(udaje, aes(x = datum, y = registracie)) +
geom_line(color = "steelblue", linewidth = 1) +
geom_point(color = "steelblue", size = 1.8) +
labs(x = "Štvrťrok", y = "Index registrácií (2021 = 100)") +
theme_minimal(base_size = 13)
ggplot(udaje, aes(x = datum, y = urok)) +
geom_line(color = "darkorange", linewidth = 1) +
geom_point(color = "darkorange", size = 1.8) +
labs(x = "Štvrťrok", y = "Úroková sadzba (%)") +
theme_minimal(base_size = 13)
3.2 Bodový graf závislosti
Očakávame, že vyššia sadzba sa spája s nižšou hodnotou registrácií. Ide zatiaľ o pracovnú hypotézu.
graf <- ggplot(udaje, aes(x = urok, y = registracie)) +
geom_point(color = "steelblue", size = 2, alpha = 0.75) +
labs(x = "Úroková sadzba (%)", y = "Index registrácií (2021 = 100)") +
theme_minimal(base_size = 13)
graf
Otázky: Vidíte rastúcu alebo klesajúcu závislosť? Majú všetky štvrťroky s podobnou sadzbou aj rovnaký počet registrácií? Z grafu je vidno veľký rozptyl, závislosť nie je jasná.
4. Jednoduchá lineárna regresia
4.1 Zápis modelu
\[ registracie_t = \beta_0 + \beta_1 \, urok_t + u_t \]
- \(t\) označuje štvrťrok,
- \(\beta_0\) je konštanta (priesečník s osou Y),
- \(\beta_1\) je smernica: zmena predikovaného indexu registrácií pri zvýšení sadzby o 1 percentuálny bod,
- \(u_t\) zachytáva časť registrácií, ktorú model nevysvetľuje.
4.2 Odhad pomocou lm()
model1 <- lm(registracie ~ urok, data = udaje)
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: registrácie a úroková sadzba") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 90.997 |
| urok | 2.009 |
Interpretácia smernice: Zvýšenie sadzby o 1 percentuálny bod sa v tomto jednoduchom modeli spája s rozdielom približne 2.01 bodu v indexe registrácií.
Interpretácia konštanty: Pri sadzbe 0 % priamka predikuje index približne 91. Nula leží v pozorovanom rozsahu sadzieb (od približne -0,6 % do 4 %), takže konštanta má tu aj vecný zmysel.
PoznámkaPoznámka
Znamienko smernice nezodpovedá pôvodnej hypotéze (kladné namiesto záporného) a model vysvetľuje len malú časť variability (pozri R² v časti 6). Toto je dôležitý výsledok, nie chyba. Úroky aj registrácie sa v čase vyvíjali spolu s ďalšími faktormi (trend, pandémia, ceny energií), ktoré jednoduchý model nezohľadňuje. Súvislosť ukazuje len spoločný pohyb, nie príčinu.
4.3 Regresná priamka v grafe
graf +
geom_abline(intercept = coef(model1)["(Intercept)"],
slope = coef(model1)["urok"],
color = "firebrick", linewidth = 1)
5. Viacnásobná lineárna regresia
Úrokové sadzby mohli ovplyvňovať registrácie spolu s dlhodobým trendom a s oneskorením (firma sa nezaloží hneď po zmene sadzby). Rozšírime model o časový trend t a o sadzbu oneskorenú o 1 štvrťrok:
\[ registracie_t = \beta_0 + \beta_1 \, urok_t + \beta_2 \, urok\_lag1_t + \beta_3 \, t + u_t \]
model2 <- lm(registracie ~ urok + urok_lag1 + t, data = udaje)
vysledky2 <- tidy(model2)
vysledky2[c("term", "estimate")] |>
kable(format = "html", digits = 3, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Viacnásobná regresia: sadzba, oneskorená sadzba a trend") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 80.498 |
| urok | -17.629 |
| urok_lag1 | 16.286 |
| t | 0.594 |
Model s oneskorením používa 46 pozorovaní, lebo prvý štvrťrok nemá oneskorenú hodnotu.
Interpretácia:
- Koeficient pri
urokhovorí o rozdiele registrácií pri zvýšení aktuálnej sadzby o 1 p. b. pri rovnakej oneskorenej sadzbe a rovnakom trende. - Koeficient pri
urok_lag1hovorí to isté pre sadzbu spred jedného štvrťroka. - Koeficient pri
tje priemerný rast indexu registrácií za jeden štvrťrok pri nezmenených sadzbách.
Všimnite si: koeficient pri urok sa po pridaní ďalších premenných výrazne zmenil oproti jednoduchému modelu (zmenilo sa aj znamienko). Jednoduchý a viacnásobný model odpovedajú na odlišné otázky. Aktuálna a oneskorená sadzba sú navyše silno prepojené (sadzba sa mení pomaly), takže ich samostatné koeficienty sú nestabilné a nedajú sa spoľahlivo interpretovať izolovane. Na toto sa vrátime pri multikolinearite.
Pre porovnanie odhadneme aj jednoduchší model s trendom:
model3 <- lm(registracie ~ urok + t, data = udaje)
tidy(model3)[c("term", "estimate")] |>
kable(format = "html", digits = 3, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Regresia: sadzba a trend") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 78.402 |
| urok | -2.182 |
| t | 0.658 |
Po pridaní trendu vychádza koeficient pri urok záporný, v súlade s pôvodnou hypotézou. Trend t v jednoduchom modeli chýbal a skresľoval odhad. Stále však ide len o spoločný pohyb premenných, nie o dôkaz príčiny.
6. Súhrn modelov pomocou glance()
porovnanie_modelov <- rbind(
cbind(Model = "1: urok", glance(model1)[c("nobs", "r.squared", "sigma")]),
cbind(Model = "3: urok + t", glance(model3)[c("nobs", "r.squared", "sigma")]),
cbind(Model = "2: urok + urok_lag1 + t", glance(model2)[c("nobs", "r.squared", "sigma")])
)
porovnanie_modelov |>
kable(format = "html", digits = c(0, 0, 3, 3), row.names = FALSE,
col.names = c("Model", "Počet pozorovaní", "R²", "Odhad SD chyby (body indexu)"),
caption = "Základné charakteristiky modelov") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Model | Počet pozorovaní | R² | Odhad SD chyby (body indexu) |
|---|---|---|---|
| 1: urok | 47 | 0.048 | 14.226 |
| 3: urok + t | 47 | 0.229 | 12.947 |
| 2: urok + urok_lag1 + t | 46 | 0.326 | 11.895 |
nobsje počet pozorovaní použitých pri odhade.r.squaredje podiel výberovej variability registrácií zachytený modelom.sigmaje odhad štandardnej odchýlky chyby modelu v bodoch indexu.
R² v modeli 1 je veľmi nízke (približne 0,05): samotná sadzba takmer nič nevysvetľuje. Pridanie trendu a oneskorenia ho zvyšuje, ale ani najlepší z modelov nevysvetľuje väčšinu variability. R² neznamená „percento správnych predikcií“ a nedokazuje kauzalitu.
7. Skutočná hodnota, predikcia a rezíduum
Pre každý štvrťrok porovnáme skutočnú hodnotu s hodnotou z modelu 3 (sadzba + trend):
\[ \hat{u}_t = registracie_t - \widehat{registracie}_t \]
porovnanie <- data.frame(
Obdobie = udaje$obdobie,
Skutocnost = udaje$registracie,
Predikcia = fitted(model3),
Reziduum = residuals(model3)
)
head(porovnanie, 8) |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Štvrťrok", "Skutočnosť", "Predikcia", "Rezíduum"),
caption = "Skutočné a predikované hodnoty indexu registrácií") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Štvrťrok | Skutočnosť | Predikcia | Rezíduum |
|---|---|---|---|
| 2015-Q1 | 65.3 | 78.96 | -13.66 |
| 2015-Q2 | 71.7 | 79.73 | -8.03 |
| 2015-Q3 | 73.2 | 80.44 | -7.24 |
| 2015-Q4 | 69.1 | 81.23 | -12.13 |
| 2016-Q1 | 72.3 | 82.10 | -9.80 |
| 2016-Q2 | 76.6 | 82.91 | -6.31 |
| 2016-Q3 | 75.4 | 83.66 | -8.26 |
| 2016-Q4 | 77.7 | 84.35 | -6.65 |
Kladné rezíduum znamená, že skutočná hodnota je vyššia než predikcia modelu, záporné znamená opak.
Úloha: Vyberte jeden riadok a overte, že Skutocnost - Predikcia = Reziduum.
# štvrťrok s najvyšším kladným a najnižším záporným rezíduom
porovnanie[which.max(porovnanie$Reziduum), ] Obdobie Skutocnost Predikcia Reziduum
24 2020-Q4 125.2 95.33988 29.86012porovnanie[which.min(porovnanie$Reziduum), ] Obdobie Skutocnost Predikcia Reziduum
33 2023-Q1 51.6 94.37034 -42.770348. Predikcia pre zadanú hodnotu
Použijeme jednoduchý model 1 a vypočítame predikciu registrácií pre sadzby -0,5 %, 0 %, 2 % a 4 % (všetky ležia v pozorovanom rozsahu):
nove_udaje <- data.frame(urok = c(-0.5, 0, 2, 4))
nove_udaje$Predikcia <- predict(model1, newdata = nove_udaje)
nove_udaje |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Úroková sadzba (%)", "Predikcia indexu registrácií"),
caption = "Predikcie jednoduchého regresného modelu") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Úroková sadzba (%) | Predikcia indexu registrácií |
|---|---|
| -0.5 | 89.99 |
| 0.0 | 91.00 |
| 2.0 | 95.02 |
| 4.0 | 99.03 |
Ide o bodové predikcie, nie o isté výsledky. Skutočný štvrťrok môže mať inú hodnotu. O neistote predikcie sa budeme učiť neskôr.
9. Malé cvičenia
1. Ktorá premenná je vysvetľovaná a ktoré sú vysvetľujúce? Vysvetľovaná premenná je registracie. Vysvetľujúce sú urok, urok_lag1 a t (podľa modelu).
2. Odhadnite jednoduchý model registracie ~ bankroty, zobrazte koeficienty a nakreslite graf s priamkou.
model_b <- lm(registracie ~ bankroty, data = udaje)
tidy(model_b)[c("term", "estimate")] |>
kable(format = "html", digits = 3, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta"),
caption = "Jednoduchá regresia: registrácie a bankroty") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta |
|---|---|
| (Intercept) | 78.666 |
| bankroty | 0.156 |
ggplot(udaje, aes(x = bankroty, y = registracie)) +
geom_point(color = "steelblue", size = 2, alpha = 0.75) +
geom_abline(intercept = coef(model_b)["(Intercept)"],
slope = coef(model_b)["bankroty"],
color = "firebrick", linewidth = 1) +
labs(x = "Index bankrotov (2021 = 100)", y = "Index registrácií (2021 = 100)") +
theme_minimal(base_size = 13)
3. Vysvetlite koeficient urok v jednoduchom aj viacnásobnom modeli. Prečo pri viacnásobnom dopĺňame „pri rovnakých ostatných premenných“? V jednoduchom modeli porovnávame štvrťroky len podľa sadzby, pričom sa ostatné faktory (trend) môžu líšiť. Vo viacnásobnom modeli porovnávame štvrťroky s rovnakým trendom (a oneskorenou sadzbou), preto koeficient odpovedá na inú otázku.
4. Predikcia pre sadzbu 1 % z modelu 1 a overenie ručným dosadením:
predict(model1, newdata = data.frame(urok = 1)) 1
93.00629 unname(coef(model1)["(Intercept)"] + coef(model1)["urok"] * 1)[1] 93.006295. Jedna krajina… resp. jeden štvrťrok s kladným a jeden so záporným rezíduom: Pozri tabuľku v časti 7. Kladné rezíduum znamená, že registrácií bolo viac, než model predpovedal, záporné, že ich bolo menej.
6. Prečo z koeficienta nemožno automaticky vyvodiť kauzálny účinok? V časovom rade sa mnohé premenné menia spoločne s časom (trend, ekonomický cyklus, pandémia, energetická kríza). Koeficient zachytáva spoločný pohyb, nie príčinu. Navyše sú naše údaje o registráciách z veľkej časti doplnené (príznak i) a vzorka má len 47 pozorovaní.
10. Záver
- Jednoduchý model s úrokovou sadzbou nepotvrdil hypotézu o zápornom vzťahu a vysvetlil veľmi málo variability.
- Po zohľadnení trendu sa koeficient pri sadzbe stal záporným, čo je v súlade s hypotézou, ale závery sú predbežné a bez testov (ich témou budú nasledujúce hodiny).
- Obmedzenia: len 47 pozorovaní, hodnoty registrácií sú z veľkej časti doplnené Eurostatom, sadzba eurozóny nie je špecifická pre Slovensko a údaje sú časový rad, pri ktorom treba v ďalších krokoch skúmať stacionaritu a autokoreláciu rezíduí.