# Spustite raz v konzole, iba ak tieto balíky chýbajú.
install.packages(c("broom", "knitr", "kableExtra", "ggplot2", "dplyr",
"tseries", "car"))Econometrics in R - cvičenie 5
Základná diagnostika modelu: Kanada, 2014–2024
Pokračujeme s ročnými makroekonomickými údajmi Kanady za roky 2014–2024 z predchádzajúceho cvičenia.
library(broom)
library(knitr)
library(kableExtra)
library(ggplot2)
library(dplyr)
# Desatinná čiarka vo výstupoch.
options(OutDec = ",")
# Pomocné funkcie na formátovanie čísel a textov.
fmt <- function(x, d = 2) formatC(x, format = "f", digits = d, decimal.mark = ",")
vyznam <- function(x) ifelse(x < 0.05, "štatisticky významný", "štatisticky nevýznamný")Príprava databázy
Údaje zadáme priamo v kóde. Jeden riadok zodpovedá jednému roku. Tabuľka je úplná, takže nie je potrebné dopĺňať chýbajúce hodnoty.
udaje <- data.frame(
Year = 2014:2024,
GDP = c(1805.7, 1556.5, 1528, 1649.3, 1725.3, 1743.7,
1655.7, 2022.4, 2190.4, 2173.3, 2243.6),
GDP.pc = c(50961, 43594, 42314, 45130, 46539, 46353,
43538, 52887, 56257, 54220, 54340),
Inflation = c(1.91, 1.13, 1.43, 1.6, 2.27, 1.95,
0.72, 3.4, 6.8, 3.88, 2.38),
Unemployment = c(7.02, 6.95, 7.04, 6.43, 5.84, 5.69,
9.66, 7.53, 5.28, 5.42, 6.35),
Trade.balance = c(-0.91, -2.46, -2.36, -2.19, -1.95, -1.47,
-2.21, 0.1, 0.57, 0.07, -0.26)
)
sum(is.na(udaje))[1] 0Uloženie všetkých objektov môjho workspace, čistenie workspace
- Ak nechcem do každého projektu dávať opätovné importovanie a čistenie databázy, uložím si celý obsah - všetky objekty - pamäte (workspace) do súboru a ten na začiatku zadania nabudúce nadčítam:
save.image("pamat.Rdata")- Mazanie niektorých objektov z pamäte:
aaa <- 6
bbb <- "text" # aaa a bbb chcem vymazať
rm(list = c("aaa", "bbb")) # vymažem tieto dva objekty
ls() # v pamäti zostalo[1] "fmt" "udaje" "vyznam"- Čistenie pamäte - celej
rm(list = ls())
ls()character(0)- Opätovne si nadčítam uložený obsah pamäte - workspace
load("pamat.Rdata")Lineárna regresia
Odhadujeme regresný model:
\[ GDPpc_t = \beta_0 + \beta_1 Unemployment_t + \beta_2 Inflation_t + \beta_3 TradeBalance_t + u_t \]
kde \(t\) označuje rok. Ukazovateľ GDP (v mld. USD) do modelu nedávame, lebo je to takmer to isté ako HDP na obyvateľa.
model <- lm(GDP.pc ~ +1 + Unemployment + Inflation + Trade.balance,
data = udaje)Výsledky regresie - odhady regresných koeficientov - vidíme v tabuľke:
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 = "Regression coefficients",
escape = TRUE
) |>
kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)| Variable | Estimate | Std. Error | t value | p-value | Significance |
|---|---|---|---|---|---|
| (Intercept) | 58052,4170 | 3499,1310 | 16,5905 | 0,0000 | *** |
| Unemployment | -437,2028 | 366,5512 | -1,1927 | 0,2718 | |
| Inflation | -363,6264 | 445,2711 | -0,8166 | 0,4410 | |
| Trade.balance | 4624,8778 | 593,2856 | 7,7954 | 0,0001 | *** |
Interpretácia výsledkov regresného modelu
Označenie štatistickej významnosti v poslednom stĺpci:
| Symbol | p-hodnota |
|---|---|
*** |
menšia ako 0,001 |
** |
menšia ako 0,01 |
* |
menšia ako 0,05 |
. |
menšia ako 0,10 |
| bez symbolu | 0,10 a viac (štatisticky nevýznamný) |
b <- coef(model)
p <- tidy(model)$p.value
names(p) <- tidy(model)$termKoeficient pri Unemployment: Ak miera nezamestnanosti vzrastie o 1 percentuálny bod, očakávané HDP na obyvateľa sa zmení približne o -437 USD, pri nezmenenej inflácii a obchodnej bilancii (ceteris paribus). Koeficient je na hladine 5 % štatisticky nevýznamný (p-hodnota 0,272).
Koeficient pri Inflation: Ak inflácia vzrastie o 1 percentuálny bod, očakávané HDP na obyvateľa sa zmení približne o -364 USD, ceteris paribus. Koeficient je štatisticky nevýznamný (p-hodnota 0,441).
Koeficient pri Trade.balance: Ak obchodná bilancia (v % HDP) vzrastie o 1 percentuálny bod, očakávané HDP na obyvateľa sa zmení približne o 4625 USD, ceteris paribus. Koeficient je štatisticky významný (p-hodnota 0,0001).
Absolútny člen (Intercept): Predstavuje predikované HDP na obyvateľa pri nulovej nezamestnanosti, inflácii a obchodnej bilancii. Má skôr technický význam, lebo nulová nezamestnanosť leží mimo pozorovaného rozsahu.
Poznámka
Máme len 11 pozorovaní, preto sú odhady citlivé na jednotlivé roky (napríklad na rok 2020 alebo 2022). Znamienka koeficientov pri inflácii a obchodnej bilancii nemusia zodpovedať učebnicovým očakávaniam a nedokazujú kauzalitu.
Kvalita vyrovnania modelu
Kvalitu vyrovnania posúdime pomocou \(R^2\), ale aj výsledku F-testu (pozri p.value vedľa F Stat.).
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
)| R.squared | Adj.R.squared | Sigma | F Stat. | p.value | df | logLik | AIC | BIC | df.residual | N |
|---|---|---|---|---|---|---|---|---|---|---|
| 0,9642 | 0,9489 | 1149,392 | 62,8738 | 0 | 3 | -90,6393 | 191,2786 | 193,268 | 7 | 11 |
g <- glance(model)Interpretácia súhrnných charakteristík
- R.squared: Model vysvetľuje približne 96,4 % variability HDP na obyvateľa. Zvyšných 3,6 % zostáva nevysvetlených.
- Adj.R.squared: Upravený koeficient determinácie je 0,949. Je o niečo nižší ako \(R^2\), lebo zohľadňuje malý počet pozorovaní voči počtu premenných.
- Sigma: Typická veľkosť odchýlky skutočného HDP na obyvateľa od modelu je približne 1149 USD.
- F Stat. a p.value: Celkový F-test overuje, či má aspoň jedna vysvetľujúca premenná nenulový koeficient. P-hodnota je 0,00002, takže model je ako celok na hladine 5 % štatisticky významný.
- df, df.residual, N: počet vysvetľujúcich premenných (3), reziduálne stupne voľnosti (7) a počet pozorovaní (11).
- logLik, AIC, BIC: Používajú sa pri porovnávaní modelov odhadnutých na rovnakých údajoch s rovnakou vysvetľovanou premennou (nižšie AIC a BIC znamená lepší model). Pri jednom modeli ich neinterpretujeme.
Charakteristiku deviance (reziduálny súčet štvorcov) sme z tabuľky vynechali. V R si môžeme overiť jej súvis so Sigma:
deviance(model)[1] 9247713sigma(model)[1] 1149,392sqrt(deviance(model) / df.residual(model))[1] 1149,392Posledné dva príkazy poskytujú rovnaký výsledok.
Grafická analýza kvality modelu
Pozrime sa na štyri diagnostické grafy modelu.
# Nastaviť rozloženie 2 x 2
par(mfrow = c(2, 2))
# Vykresliť všetky 4 diagnostické grafy modelu
plot(model)
# Resetovať layout
par(mfrow = c(1, 1))Na čo sa pri jednotlivých grafoch pozrieť:
- Residuals vs. fitted: Rezíduá by mali kolísať okolo nuly bez zjavného tvaru. Pozrite, či sa červená čiara výrazne zakrivuje a ktoré roky ležia ďaleko od ostatných.
- Q-Q plot: Body by mali ležať približne na prerušovanej priamke. Pozrite, ktoré roky sa od nej odchyľujú na koncoch.
- Scale-Location: Body by mali byť rozptýlené rovnomerne bez lievikového tvaru.
- Residuals vs. leverage: Pozrite, či niektorý rok nevybočuje vpravo (vysoký vplyv) alebo neprekračuje čiarkované hranice Cookovej vzdialenosti.
cook <- cooks.distance(model)
paka <- hatvalues(model)
rok_cook <- udaje$Year[which.max(cook)]V odhadnutom modeli má najväčšiu Cookovu vzdialenosť rok 2022 (hodnota 1,48, pákový efekt 0,88). Hodnota Cookovej vzdialenosti nad 1 naznačuje, že rok výrazne ovplyvňuje odhad koeficientov. Pri 11 pozorovaniach je to pochopiteľné.
Testy normality rezíduí
Jarqueov–Berov test
Testujeme hypotézu, že rezíduá majú normálne rozdelenie (\(H_0\)), proti alternatíve, že nie sú normálne rozdelené (\(H_1\)). Ak je p-hodnota menšia ako 0,05, normalitu zamietame.
res <- residuals(model)
jb_test <- tseries::jarque.bera.test(res)
broom::tidy(jb_test) |>
dplyr::select(statistic, p.value, method) |>
knitr::kable(
digits = 4,
col.names = c("JB statistic", "p-value", "Test"),
caption = "Jarque-Bera test of residual normality"
) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)| JB statistic | p-value | Test |
|---|---|---|
| 1,1482 | 0,5632 | Jarque Bera Test |
P-hodnota testu je 0,563, takže hypotézu o normalite rezíduí na hladine 5 % nezamietame.
Poznámka: Nezamietnutie nulovej hypotézy neznamená, že normalitu rezíduí sme dokázali. Pri 11 pozorovaniach má test malú silu, preto výsledok doplňte pohľadom na Q-Q graf.
Test odľahlých pozorovaní – outlierTest()
Funkcia car::outlierTest() identifikuje pozorovania, ktorých hodnota HDP na obyvateľa sa výrazne odlišuje od hodnoty predikovanej modelom. Rozhodujeme sa podľa Bonferroniho p-hodnoty: ak je menšia ako 0,05, rok považujeme za štatisticky významne odľahlý.
outlier_test <- car::outlierTest(model)
if (is.null(outlier_test)) {
cat("No observations returned by outlierTest.")
} else {
outlier_table <- data.frame(
Observation = names(outlier_test$rstud),
rstudent = unname(outlier_test$rstud),
unadjusted.p = unname(outlier_test$p),
bonf.p = unname(outlier_test$bonf.p)
)
outlier_table |>
knitr::kable(
digits = 4,
col.names = c(
"Observation",
"Studentized residual",
"Unadjusted p-value",
"Bonferroni p-value"
),
caption = "Outlier test with Bonferroni correction"
) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)
}| Observation | Studentized residual | Unadjusted p-value | Bonferroni p-value |
|---|---|---|---|
| 6 | -2,0763 | 0,0832 | 0,9148 |
Číslo v stĺpci Observation je poradie riadku v tabuľke udaje, teda 1 zodpovedá roku 2014, 6 roku 2019 atď.
Úprava modelu
Zlogaritmujeme infláciu (je vždy kladná, takže logaritmus existuje):
\[ GDPpc_t = \beta_0 + \beta_1 Unemployment_t + \beta_2 \ln(Inflation_t) + \beta_3 TradeBalance_t + u_t \]
Logaritmická transformácia mení interpretáciu koeficientu pri inflácii: ak sa inflácia zmení o 1 %, očakávané HDP na obyvateľa sa zmení približne o \(\beta_2/100\) USD. Obchodnú bilanciu nelogaritmujeme, lebo obsahuje záporné hodnoty.
model2 <- lm(
GDP.pc ~ Unemployment + log(Inflation) + Trade.balance,
data = udaje
)
tidy(model2) |>
kable(
digits = 4,
col.names = c(
"Variable",
"Estimate",
"Std. Error",
"t value",
"Pr(>|t|)"
),
caption = "Regression results – GDP per capita"
) |>
kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)| Variable | Estimate | Std. Error | t value | Pr(>|t|) |
|---|---|---|---|---|
| (Intercept) | 61226,6957 | 5819,6388 | 10,5207 | 0,0000 |
| Unemployment | -749,0418 | 551,1315 | -1,3591 | 0,2163 |
| log(Inflation) | -2089,8243 | 2078,6491 | -1,0054 | 0,3482 |
| Trade.balance | 5015,6448 | 849,7968 | 5,9022 | 0,0006 |
Koeficient pri log(Inflation) je približne -2090. Zvýšenie inflácie o 1 % sa teda spája so zmenou očakávaného HDP na obyvateľa približne o -20,9 USD, ceteris paribus.
Grafická diagnostika upraveného modelu
par(mfrow = c(2, 2))
plot(model2)
par(mfrow = c(1, 1))Test normality rezíduí upraveného modelu
res2 <- residuals(model2)
jb_test2 <- tseries::jarque.bera.test(res2)
broom::tidy(jb_test2) |>
dplyr::select(statistic, p.value, method) |>
knitr::kable(
digits = 4,
col.names = c("JB statistic", "p-value", "Test"),
caption = "Jarque-Bera test of residual normality – upravený model"
) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)| JB statistic | p-value | Test |
|---|---|---|
| 0,6824 | 0,7109 | Jarque Bera Test |
P-hodnota je 0,711, takže hypotézu o normalite rezíduí upraveného modelu na hladine 5 % nezamietame.
Test odľahlých hodnôt upraveného modelu
outlier_test2 <- car::outlierTest(model2)
if (is.null(outlier_test2)) {
cat("No observations returned by outlierTest.")
} else {
outlier_table2 <- data.frame(
Observation = names(outlier_test2$rstud),
rstudent = unname(outlier_test2$rstud),
unadjusted.p = unname(outlier_test2$p),
bonf.p = unname(outlier_test2$bonf.p)
)
outlier_table2 |>
knitr::kable(
digits = 4,
col.names = c(
"Observation",
"Studentized residual",
"Unadjusted p-value",
"Bonferroni p-value"
),
caption = "Outlier test with Bonferroni correction – upravený model"
) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)
}| Observation | Studentized residual | Unadjusted p-value | Bonferroni p-value |
|---|---|---|---|
| 6 | -2,6781 | 0,0366 | 0,4029 |
Malé úlohy
1. Porovnanie modelov pomocou AIC
Obidva modely majú rovnakú vysvetľovanú premennú (GDP.pc) a sú odhadnuté na rovnakých pozorovaniach, preto ich môžeme porovnať pomocou AIC. Lepší je model s nižšou hodnotou.
AIC(model, model2) |>
kable(
digits = 2,
col.names = c("Stupne voľnosti", "AIC"),
caption = "Porovnanie modelov pomocou AIC"
) |>
kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)| Stupne voľnosti | AIC | |
|---|---|---|
| model | 5 | 191,28 |
| model2 | 5 | 190,80 |
aic1 <- AIC(model)
aic2 <- AIC(model2)AIC pôvodného modelu je 191,28 a AIC upraveného modelu je 190,80. Podľa AIC je lepší upravený model s log(Inflation). Rozdiel je však len 0,48, čo pri 11 pozorovaniach neznamená výraznú výhodu.
2. Iná modifikácia modelu
Inflácia vyšla v pôvodnom modeli štatisticky nevýznamná. Skúsime ju z modelu vynechať a porovnať všetky tri modely:
model3 <- lm(GDP.pc ~ Unemployment + Trade.balance, data = udaje)
tidy(model3) |>
kable(
digits = 4,
col.names = c("Variable", "Estimate", "Std. Error", "t value", "Pr(>|t|)"),
caption = "Regression results – model bez inflácie"
) |>
kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)| Variable | Estimate | Std. Error | t value | Pr(>|t|) |
|---|---|---|---|---|
| (Intercept) | 55730,2657 | 1996,2588 | 27,9174 | 0,0000 |
| Unemployment | -294,5532 | 315,4730 | -0,9337 | 0,3778 |
| Trade.balance | 4233,7935 | 342,8433 | 12,3491 | 0,0000 |
AIC(model, model2, model3) |>
kable(
digits = 2,
col.names = c("Stupne voľnosti", "AIC"),
caption = "Porovnanie troch modelov pomocou AIC"
) |>
kable_styling(
bootstrap_options = c("striped", "hover"),
full_width = FALSE
)| Stupne voľnosti | AIC | |
|---|---|---|
| model | 5 | 191,28 |
| model2 | 5 | 190,80 |
| model3 | 4 | 190,28 |
aic_vsetky <- c("pôvodný model" = AIC(model),
"model s log(Inflation)" = AIC(model2),
"model bez inflácie" = AIC(model3))Najnižšie AIC (190,28) má model bez inflácie. Rozdiely medzi modelmi sú malé, takže výber nepovažujeme za jednoznačný.
Záver
V tejto úlohe sme odhadli viacnásobnú lineárnu regresiu pre Kanadu v rokoch 2014–2024 a predviedli základnú diagnostiku modelu. Model dobre vyrovnáva údaje (R² približne 0,96), ale štatisticky významná je predovšetkým obchodná bilancia, nie nezamestnanosť či inflácia. Pri 11 pozorovaniach sú testy málo silné a jednotlivé roky (najmä 2022) môžu odhad výrazne ovplyvniť. Logaritmická transformácia inflácie zlepšila AIC len minimálne, takže výsledky treba brať ako ukážku postupu, nie ako spoľahlivý ekonomický záver.
save.image("pamat.Rdata")
rm(list = ls())