Ekonometria Master dokument
Filip Jurkáček
Poznámka: Pred kompiláciou majte v pracovnom priečinku súbor Databaza.csv (rovnaký ako v zadaniach).
1 Cvičenie 4: Práca s databázou - import údajov, grafy, štatistiky
Tento report načíta súbor Databaza.csv (stĺpce: Okres, Obec, 2024M12, …, 2024M01), prevedie údaje do “long” formátu a vykoná základné grafy, štatistiky, testy a regresiu.
2 Balíčky a dáta
suppressPackageStartupMessages({
library(dplyr)
library(readr)
library(tidyr)
library(ggplot2)
library(knitr)
library(kableExtra)
library(broom)
library(stringr)
})# CSV je bodkočiarkové
udaje_raw <- read_delim("Databaza.csv", delim = ";", show_col_types = FALSE)
# validácie
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šiel som žiadne stĺpce typu 2024Mxx (napr. 2024M12, 2024M01).")
}
# Long formát
udaje_long <- udaje_raw |>
pivot_longer(
cols = all_of(mesiac_cols),
names_to = "YM",
values_to = "Hodnota"
) |>
mutate(
Rok = as.integer(str_sub(YM, 1, 4)),
Mesiac = as.integer(str_sub(YM, 6, 7)),
Hodnota = suppressWarnings(as.numeric(Hodnota))
)
# Wide podmnožina pre jednoduché modely
udaje_dec <- udaje_raw |>
dplyr::select(Okres, Obec, any_of(c("2024M12","2024M11","2024M06","2024M01"))) |>
mutate(across(-c(Okres, Obec), ~ suppressWarnings(as.numeric(.))))3 Vizualizácie
3.1 Scatter plot: 2024M11 vs 2024M12
if (all(c("2024M12","2024M11") %in% names(udaje_dec))) {
ggplot(udaje_dec, aes(x = `2024M11`, y = `2024M12`, colour = Okres)) +
geom_point(alpha = 0.8) +
theme_minimal() +
labs(
title = "Vzťah medzi 2024M11 a 2024M12 naprieč obcami",
x = "Hodnota v 2024M11",
y = "Hodnota v 2024M12",
colour = "Okres"
)
} else {
plot.new(); title("Na scatter plot chýbajú 2024M11 a/alebo 2024M12")
}
3.2 Boxplot: rozdelenie hodnôt podľa mesiacov
ggplot(udaje_long, aes(x = factor(Mesiac), y = Hodnota)) +
geom_boxplot(fill = "lightblue", colour = "darkblue") +
labs(
title = "Rozdelenie hodnôt podľa mesiacov (rok 2024)",
x = "Mesiac",
y = "Hodnota"
) +
theme_minimal()
4 Základné štatistiky
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 tabuľka
kable(stat_tbl, digits = 2, caption = "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 |
# kableExtra tabuľka
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))| 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 |
5 Testovanie hypotéz
5.1 t-test: porovnanie priemerov (2024M12 vs 2024M06)
if (all(c("2024M12","2024M06") %in% names(udaje_raw))) {
t.test.result <- t.test(
udaje_long$Hodnota[udaje_long$YM == "2024M12"],
udaje_long$Hodnota[udaje_long$YM == "2024M06"]
)
t.test.result
} else {
"Na t-test chýba 2024M12 a/alebo 2024M06."
}##
## Welch Two Sample t-test
##
## data: udaje_long$Hodnota[udaje_long$YM == "2024M12"] and udaje_long$Hodnota[udaje_long$YM == "2024M06"]
## 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.615.2 ANOVA: Hodnota ~ Mesiac
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 1404398586 Lineárna regresia
Predikujeme 2024M12 pomocou 2024M11 a 2024M01 (ak sú dostupné).
if (all(c("2024M12","2024M11","2024M01") %in% names(udaje_dec))) {
model <- lm(`2024M12` ~ `2024M11` + `2024M01`, data = udaje_dec)
summary(model)
} else {
model <- NULL
"Na regresiu chýbajú niektoré z 2024M12, 2024M11, 2024M01."
}##
## Call:
## lm(formula = `2024M12` ~ `2024M11` + `2024M01`, data = udaje_dec)
##
## Residuals:
## Min 1Q Median 3Q Max
## -51.294 -2.214 0.940 3.381 31.350
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.778536 0.633074 -2.809 0.00546 **
## `2024M11` 1.066640 0.006654 160.291 < 2e-16 ***
## `2024M01` -0.066194 0.006661 -9.938 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.195 on 197 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 2.085e+08 on 2 and 197 DF, p-value: < 2.2e-166.1 Koeficienty a intervaly spoľahlivosti
if (!is.null(model)) {
coef.tbl <- broom::tidy(model, conf.int = TRUE) |>
dplyr::mutate(
term = dplyr::recode(term,
"(Intercept)" = "Intercept",
"`2024M11`" = "Hodnota v 2024M11",
"`2024M01`" = "Hodnota v 2024M01"),
stars = dplyr::case_when(
p.value < 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "*",
p.value < 0.1 ~ "·",
TRUE ~ ""
)
) |>
dplyr::transmute(
Term = 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 = "OLS koeficienty: 2024M12 ~ 2024M11 + 2024M01") |>
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 nefitol).")
}| Term | Estimate | Std. Error | t value | p value | 95% CI | Sig |
|---|---|---|---|---|---|---|
| Intercept | -1.779 | 0.633 | -2.809 | 0.005 | [-3.027, -0.53] | ** |
| Hodnota v 2024M11 | 1.067 | 0.007 | 160.291 | 0.000 | [1.054, 1.08] | *** |
| Hodnota v 2024M01 | -0.066 | 0.007 | -9.938 | 0.000 | [-0.079, -0.053] | *** |
| Note: | ||||||
| Signif. codes: *** p<0.001, ** p<0.01, * p<0.05, · p<0.1. |
6.2 Štatistiky kvality modelu
if (!is.null(model)) {
fit.tbl <- broom::glance(model) |>
dplyr::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 nefitol).")
}| R-squared | Adj. R-squared | F-statistic | F p-value | AIC | BIC | Num. obs. |
|---|---|---|---|---|---|---|
| 1 | 1 | 208498789 | 0 | 1413.939 | 1427.133 | 200 |
7 Cvičenie 6 (v podkladoch označené ako Econometrics in R - cvičenie 5)
7.1 Cieľ
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.
7.2 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)
})7.3 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"7.3.1 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 Bratislava …
## $ Obec <chr> "Bratislava - mestská časť Staré Mesto", "Bratislava - mestská…
## $ YM <chr> "2024M12", "2024M11", "2024M10", "2024M09", "2024M08", "2024M0…
## $ Hodnota <dbl> 47634, 47632, 47565, 47564, 47525, 47488, 47503, 47484, 47411,…
## $ Rok <int> 2024, 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, 7, 6,…7.4 Vizualizácie
7.4.1 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")
7.4.2 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")
}
7.5 Základné štatistiky
7.5.1 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)")| 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))| 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 |
7.6 Testovanie hypotéz
7.6.1 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.617.6.2 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 1404398587.7 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-167.7.1 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)."
}
## rstudent unadjusted p-value Bonferroni p
## 9 -6.805608 1.2082e-10 2.4164e-08
## 16 4.992611 1.3175e-06 2.6349e-04
## 3 4.939881 1.6768e-06 3.3536e-04
## 1 -4.767936 3.6350e-06 7.2701e-04
## 5 4.168518 4.6075e-05 9.2150e-03
## 10 -3.933493 1.1633e-04 2.3266e-027.7.2 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).")
}| 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. |
7.7.3 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).")
}| R-squared | Adj. R-squared | F-statistic | F p-value | AIC | BIC | Num. obs. |
|---|---|---|---|---|---|---|
| 1 | 1 | 141101301 | 0 | 1411.92 | 1428.411 | 200 |
7.8 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.
8 Cvičenie 7: Špecifikácia modelu
library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
# rm(list = ls()) # vypnuté v spoločnom reportePokračujeme v predchádzajúcich úlohách, kde s pomocou regresie skúmame, nakoľko je počet obyvateľov v obci na konci roka 2024 (stav v decembri 2024) ovplyvnený počtom obyvateľov:
- na začiatku roka (január 2024),
- v polovici roka (jún 2024),
- na začiatku jesene (september 2024).
Pracujeme teda s prierezovými dátami za obce Slovenska, pričom za každú obec máme stav obyvateľov v jednotlivých mesiacoch roka 2024. Výsledok základnej regresie je uvedený nižšie.
#######################################################################
# PRÍPRAVA ÚDAJOV
#######################################################################
# načítame databázu (tvoj súbor Databaza.csv)
# check.names = FALSE zabezpečí, že stĺpce zostanú pomenované 2024M12, 2024M01 atď.
udaje_raw <- read.csv("Databaza.csv", dec = ".", sep = ";", header = TRUE,
check.names = FALSE)
# skontrolujeme štruktúru
str(udaje_raw)## 'data.frame': 200 obs. of 14 variables:
## $ Okres : chr "Okres Bratislava I" "Okres Bratislava II" "Okres Bratislava II" "Okres Bratislava II" ...
## $ Obec : chr "Bratislava - mestská časť Staré Mesto" "Bratislava - m. č. Podunajské Biskupice" "Bratislava - mestská časť Ružinov" "Bratislava - mestská časť Vrakuňa" ...
## $ 2024M12: int 47634 23095 83365 20172 45539 26481 6050 17065 35538 34815 ...
## $ 2024M11: int 47632 23110 83251 20161 45480 26466 6052 17071 35575 34831 ...
## $ 2024M10: int 47565 23130 83192 20159 45444 26429 6045 17040 35579 34816 ...
## $ 2024M09: int 47564 23132 83163 20162 45455 26432 6049 17036 35636 34788 ...
## $ 2024M08: int 47525 23137 83137 20161 45442 26430 6040 17027 35649 34748 ...
## $ 2024M07: int 47488 23147 83056 20190 45442 26394 6036 17021 35627 34782 ...
## $ 2024M06: int 47503 23162 82867 20216 45411 26376 6030 16999 35633 34812 ...
## $ 2024M05: int 47484 23178 82789 20211 45376 26346 6024 17005 35624 34843 ...
## $ 2024M04: int 47411 23179 82641 20227 45414 26291 6031 17027 35600 34854 ...
## $ 2024M03: int 47408 23201 82612 20209 45439 26245 6034 17007 35575 34899 ...
## $ 2024M02: int 47391 23226 82553 20196 45410 26231 6037 17030 35570 34924 ...
## $ 2024M01: int 47375 23276 82483 20221 45342 26214 6038 17067 35572 34942 ...head(udaje_raw)# Vytvoríme si premenné pre jednotlivé mesiace, s ktorými budeme pracovať v modeli:
# obyv_12 – počet obyvateľov v decembri 2024 (vysvetľovaná premenná)
# obyv_01 – počet obyvateľov v januári 2024
# obyv_06 – počet obyvateľov v júni 2024
# obyv_09 – počet obyvateľov v septembri 2024
udaje <- within(udaje_raw, {
obyv_12 <- `2024M12`
obyv_01 <- `2024M01`
obyv_06 <- `2024M06`
obyv_09 <- `2024M09`
})
# vyberieme len premenné, ktoré idú do modelu
udaje_model <- udaje[, c("obyv_12", "obyv_01", "obyv_06", "obyv_09")]
# (prípadná) imputácia chýbajúcich hodnôt – v našej databáze zrejme chýbajúce údaje nie sú,
# ale pre ukážku ponecháme postup s doplnením mediánov
column_medians <- sapply(udaje_model, median, na.rm = TRUE)
udaje_imputed <- udaje_model
for (col in names(udaje_model)) {
udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]
}
udaje <- udaje_imputed
str(udaje)## 'data.frame': 200 obs. of 4 variables:
## $ obyv_12: num 47634 23095 83365 20172 45539 ...
## $ obyv_01: num 47375 23276 82483 20221 45342 ...
## $ obyv_06: num 47503 23162 82867 20216 45411 ...
## $ obyv_09: num 47564 23132 83163 20162 45455 ...8.0.1 Základná regresia
V základnom modeli predpokladáme, že počet obyvateľov v obci v decembri 2024 je lineárne vysvetlený počtom obyvateľov v januári, júni a septembri:
\[ \text{obyv\_12}_i = \beta_0 + \beta_1 \text{obyv\_01}_i + \beta_2 \text{obyv\_06}_i + \beta_3 \text{obyv\_09}_i + u_i \]
kde index \(i\) označuje obec.
#######################################################################
# ZÁKLADNÁ REGRESIA
#######################################################################
attach(udaje)
model <- lm(obyv_12 ~ 1 + obyv_01 + obyv_06 + obyv_09, data = udaje)
summary(model)##
## Call:
## lm(formula = obyv_12 ~ 1 + obyv_01 + obyv_06 + obyv_09, data = udaje)
##
## Residuals:
## Min 1Q Median 3Q Max
## -118.937 -4.211 0.368 5.820 41.782
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.34570 1.10739 -0.312 0.755242
## obyv_01 -0.11970 0.03708 -3.228 0.001459 **
## obyv_06 -0.34533 0.09994 -3.455 0.000673 ***
## obyv_09 1.46539 0.06681 21.934 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.26 on 196 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 4.588e+07 on 3 and 196 DF, p-value: < 2.2e-16Interpretácia: koeficienty pri premenných obyv_01,
obyv_06 a obyv_09 hovoria, o koľko sa v
priemere zmení počet obyvateľov v decembri, ak sa daná vysvetľujúca
premenná zvýši o jednu osobu (pri ostatných premenných nezmenených).
8.1 1. Test RESET (test špecifikácie funkčnej formy)
Chceme otestovať, či je lineárna špecifikácia vhodná, t. j. či medzi premennými nie je výraznejší nelineárny vzťah, ktorý by sme mali zachytiť pomocou transformácií (napr. kvadratické členy, logaritmy).
Základný model:
\[ y_i = \beta_0 + \beta_1 x_{1i} + \dots + \beta_k x_{ki} + u_i \]
Ak je model správne špecifikovaný, pridaním mocnín vyrovnaných hodnôt (napr. \(\hat{y}_i^2\), \(\hat{y}_i^3\)) by sa model nemal štatisticky významne zlepšiť. RESET test formálne testuje:
- \(H_0\): model je správne špecifikovaný,
- \(H_1\): model je nesprávne špecifikovaný (chybí nelinearita alebo ďalšie premenné).
#######################################################################
# RESET TEST
#######################################################################
library(lmtest)
resettest(model)##
## RESET test
##
## data: model
## RESET = 8.0852, df1 = 2, df2 = 194, p-value = 0.000424Interpretácia:
Ak p-hodnota < 0,05, zamietame \(H_0\) a usudzujeme, že lineárny model môže
byť nesprávne špecifikovaný – buď chýbajú niektoré vysvetľujúce
premenné, alebo je potrebná nelineárna transformácia premenných.
Ak p-hodnota ≥ 0,05, nemáme dôvod tvrdiť, že je špecifikácia modelu
nevhodná.
8.2 2. Grafická analýza
8.2.1 Graf Residuals vs. Fitted
Graf rezíduí voči vyrovnaným hodnotám nám pomáha odhaliť nelinearity, heteroskedasticitu či iné systematické vzory v rezíduách.
plot(model, which = 1)
Ak rezíduá vykazujú nenáhodný vzor (napr. zakrivenie), lineárna špecifikácia môže byť nevhodná a má zmysel uvažovať o transformácii niektorej premennej alebo o zaradení ďalšej premennej.
8.2.2 Grafy C+R (Component + Residual Plots)
Tieto grafy nám pomáhajú odhadnúť, pri ktorej vysvetľujúcej premennej by mohol byť vzťah s vysvetľovanou premennou nelineárny.
car::crPlots(model)
Vertikálna os zobrazuje hodnoty typu \(\hat{\beta}_i x_{i} + e_i\), horizontálna os zobrazuje hodnotu príslušnej vysvetľujúcej premennej.
Ak sa hladká krivka (lokálna aproximácia) výrazne odchyľuje od
priamky, zvážime transformáciu danej premennej (napr. kvadratický člen,
logaritmus a pod.). V našom prípade nás môžu zaujímať najmä premené
obyv_06 (počet obyvateľov v júni) a obyv_09
(počet obyvateľov v septembri), ak sa pri nich objavuje výraznejšie
zakrivenie.
8.3 3. Nelineárna špecifikácia
Často môžeme nelineárne vzťahy približovať pomocou polynómov. V prípade kvadratických členov má model tvar:
\[ y_i = \beta_0 + \beta_1 x_{1i} + \dots + \beta_k x_{ki} + \gamma_i x_{ki}^2 + u_i \]
V našom príklade môžeme zaradiť kvadrát niektorých mesačných stavov
obyvateľov, napríklad obyv_06 (jún) a obyv_09
(september).
8.4 4. Porovnanie základného a modifikovaného modelu
Predpokladajme, že podľa C+R grafov sa rozhodneme pridať kvadrát
premenných obyv_06 a obyv_09. Náš rozšírený
model bude:
#######################################################################
# MODELY S KVADRATICKÝMI ČLENMI
#######################################################################
model_linear <- lm(obyv_12 ~ 1 + obyv_01 + obyv_06 + obyv_09)
model_kvadr <- lm(obyv_12 ~ 1 + obyv_01 + obyv_06 + obyv_09 +
I(obyv_06^2) + I(obyv_09^2))
summary(model_kvadr)##
## Call:
## lm(formula = obyv_12 ~ 1 + obyv_01 + obyv_06 + obyv_09 + I(obyv_06^2) +
## I(obyv_09^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -108.896 -4.359 -0.382 5.112 57.656
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.075e+00 1.189e+00 0.904 0.366992
## obyv_01 -1.471e-01 3.627e-02 -4.056 7.21e-05 ***
## obyv_06 -4.099e-01 1.127e-01 -3.637 0.000353 ***
## obyv_09 1.557e+00 8.553e-02 18.201 < 2e-16 ***
## I(obyv_06^2) 9.685e-07 4.695e-07 2.063 0.040454 *
## I(obyv_09^2) -9.601e-07 4.699e-07 -2.043 0.042392 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.62 on 194 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 3.017e+07 on 5 and 194 DF, p-value: < 2.2e-16anova(model_linear, model_kvadr)resettest(model_kvadr)##
## RESET test
##
## data: model_kvadr
## RESET = 12.357, df1 = 2, df2 = 192, p-value = 8.949e-06Porovnávame:
- upravený koeficient determinácie \(R^2_{adj}\),
- ANOVA test dvoch vnoreneých modelov (lineárny vs. s kvadrátmi),
- výsledok RESET testu pre rozšírený model.
Ak rozšírený model vykazuje vyšší \(R^2_{adj}\), ANOVA test je významný a RESET test poukazuje na lepšiu špecifikáciu, môžeme kvadratické členy považovať za vhodné.
Ak by napríklad kvadrát jednej premennej vyšiel štatisticky nevýznamný, môžeme model zjednodušiť a ponechať len tie kvadratické členy, ktoré zlepšujú štatistické vlastnosti modelu.
# Príklad zjednodušenia – ponechajme len kvadrát obyv_06
model_kvadr_simple <- lm(obyv_12 ~ 1 + obyv_01 + obyv_06 + obyv_09 +
I(obyv_06^2))
summary(model_kvadr_simple)##
## Call:
## lm(formula = obyv_12 ~ 1 + obyv_01 + obyv_06 + obyv_09 + I(obyv_06^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -108.106 -4.552 -0.013 5.115 57.819
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.604e+00 1.170e+00 1.370 0.17217
## obyv_01 -1.508e-01 3.652e-02 -4.128 5.42e-05 ***
## obyv_06 -2.906e-01 9.719e-02 -2.990 0.00315 **
## obyv_09 1.441e+00 6.462e-02 22.300 < 2e-16 ***
## I(obyv_06^2) 9.292e-09 2.298e-09 4.044 7.58e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.74 on 195 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 3.711e+07 on 4 and 195 DF, p-value: < 2.2e-168.5 5. Rozšírený model s kvadratickými členmi
Môžeme ešte skúsiť systematickejšiu rozšírenú špecifikáciu, kde pridáme kvadratické členy ku všetkým vysvetľujúcim premenným:
#######################################################################
# ROZŠÍRENÝ MODEL S VIACERÝMI KVADRATICKÝMI ČLENMI
#######################################################################
model_rozsireny <- lm(obyv_12 ~ obyv_01 + obyv_06 + obyv_09 +
I(obyv_01^2) + I(obyv_06^2) + I(obyv_09^2))
summary(model_rozsireny)##
## Call:
## lm(formula = obyv_12 ~ obyv_01 + obyv_06 + obyv_09 + I(obyv_01^2) +
## I(obyv_06^2) + I(obyv_09^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -93.720 -4.850 0.383 5.847 56.068
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.625e-02 1.183e+00 0.014 0.989049
## obyv_01 -3.015e-01 5.350e-02 -5.636 6.09e-08 ***
## obyv_06 -3.479e-02 1.466e-01 -0.237 0.812710
## obyv_09 1.336e+00 1.008e-01 13.258 < 2e-16 ***
## I(obyv_01^2) 5.639e-06 1.476e-06 3.821 0.000179 ***
## I(obyv_06^2) -1.306e-05 3.699e-06 -3.530 0.000519 ***
## I(obyv_09^2) 7.418e-06 2.239e-06 3.313 0.001102 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.17 on 193 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 2.69e+07 on 6 and 193 DF, p-value: < 2.2e-16anova(model_linear, model_rozsireny)resettest(model_rozsireny)##
## RESET test
##
## data: model_rozsireny
## RESET = 15.208, df1 = 2, df2 = 191, p-value = 7.436e-07Tu sledujeme, či kvadratické členy vychádzajú štatisticky významné a či ANOVA test ukáže štatisticky významné zlepšenie oproti pôvodnému lineárnemu modelu. Výsledok RESET testu nám zároveň napovie, či aj po týchto úpravách zostáva problém nesprávnej špecifikácie.
8.6 6. Transformácia pomocou dummy premennej a lineárnej lomenej funkcie
Predpokladajme, že podľa C+R grafu pri premennej obyv_06
(počet obyvateľov v júni) vidíme, že pre menšie obce a
väčšie obce môže byť vzťah odlišný – teda pre určité
„kritické“ veľkosti obcí sa mení sklon.
Zaveďme dummy premennú:
\[ DUM_i = \begin{cases} 1, & \text{ak obec má v júni viac ako 10\,000 obyvateľov}, \\ 0, & \text{inak}. \end{cases} \]
Pomocou tejto premennej môžeme modelovať:
- Zlom v autonómnom člene – posun regresnej nadroviny medzi menšími a väčšími obcami.
- Zlom v sklone – odlišný vplyv júnového počtu obyvateľov na december podľa veľkosti obce.
#######################################################################
# DUMMY PREMENNÁ A LOMENÁ FUNKCIA
#######################################################################
# definujeme dummy podľa júnového stavu – hranica 10 000 obyvateľov je ilustračná
DUM <- ifelse(obyv_06 > 10000, 1, 0)
udaje$DUM <- DUM
# Model so zlomom v autonómnom člene
modelD_auto <- lm(obyv_12 ~ 1 + DUM + obyv_01 + obyv_06 + obyv_09, data = udaje)
summary(modelD_auto)##
## Call:
## lm(formula = obyv_12 ~ 1 + DUM + obyv_01 + obyv_06 + obyv_09,
## data = udaje)
##
## Residuals:
## Min 1Q Median 3Q Max
## -115.108 -4.408 0.171 5.410 46.743
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.24123 1.08877 -0.222 0.824885
## DUM -13.85798 4.90100 -2.828 0.005180 **
## obyv_01 -0.12770 0.03654 -3.495 0.000588 ***
## obyv_06 -0.31432 0.09881 -3.181 0.001708 **
## obyv_09 1.44264 0.06614 21.812 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.02 on 195 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 3.564e+07 on 4 and 195 DF, p-value: < 2.2e-16# Model so zlomom v sklone (interakcia DUM * obyv_06)
modelD_sklon <- lm(obyv_12 ~ 1 + obyv_01 + obyv_06 + I(DUM * obyv_06) + obyv_09,
data = udaje)
summary(modelD_sklon)##
## Call:
## lm(formula = obyv_12 ~ 1 + obyv_01 + obyv_06 + I(DUM * obyv_06) +
## obyv_09, data = udaje)
##
## Residuals:
## Min 1Q Median 3Q Max
## -119.941 -3.962 0.428 5.884 41.831
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6221606 1.5304139 0.407 0.684798
## obyv_01 -0.1257484 0.0376748 -3.338 0.001012 **
## obyv_06 -0.3363435 0.1004595 -3.348 0.000977 ***
## I(DUM * obyv_06) 0.0005266 0.0005745 0.917 0.360462
## obyv_09 1.4619095 0.0669453 21.837 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.27 on 195 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 3.438e+07 on 4 and 195 DF, p-value: < 2.2e-16anova(model_linear, modelD_sklon)resettest(modelD_sklon)##
## RESET test
##
## data: modelD_sklon
## RESET = 7.9324, df1 = 2, df2 = 193, p-value = 0.000489Interpretácia:
- v modeli
modelD_autoskúmame, či sa líši úroveň počtu obyvateľov v decembri (autonómny člen) medzi menšími a väčšími obcami, - v modeli
modelD_sklonskúmame, či sa líši sklon pri premennejobyv_06(t. j. či nárast počtu obyvateľov v júni o jednu osobu zvyšuje počet obyvateľov v decembri rovnakým spôsobom v malých aj veľkých obciach).
Ak ANOVA test medzi model_linear a
modelD_sklon vyjde štatisticky významný a interakčný člen
je významný, môžeme usúdiť, že zavedenie zlomu v sklone zlepšilo
model.
8.7 7. Box–Coxov transformačný test
Box–Coxov test nám pomáha rozhodnúť, či je vhodné transformovať
vysvetľovanú premennú (u nás obyv_12),
napríklad pomocou logaritmu, odmocniny alebo reciproka.
#######################################################################
# BOX–COXOVA TRANSFORMÁCIA
#######################################################################
library(MASS)
boxcox(model_linear)
Interpretácia:
- ak \(\lambda \approx 1\) – transformácia nie je potrebná,
- ak \(\lambda \approx 0\) – vhodná je logaritmická transformácia \(\log(y)\),
- ak \(\lambda \approx 0{,}5\) – vhodná je odmocnina \(\sqrt{y}\),
- ak \(\lambda \approx -1\) – vhodná je recipročná transformácia \(1/y\).
Predpokladajme, že Box–Coxov test naznačí nejakú hodnotu \(\lambda\); ilustračné nastavenie:
# Ilustračná transformácia s lambda = 0.5 (odmocnina), ak by nám to Box–Cox odporučil
model_lambda <- lm(I((obyv_12^0.5 - 1) / 0.5) ~ 1 + obyv_01 + obyv_06 + obyv_09,
data = udaje)
summary(model_lambda)##
## Call:
## lm(formula = I((obyv_12^0.5 - 1)/0.5) ~ 1 + obyv_01 + obyv_06 +
## obyv_09, data = udaje)
##
## Residuals:
## Min 1Q Median 3Q Max
## -175.16 -21.13 -4.00 16.55 98.26
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 69.71777 2.51086 27.766 < 2e-16 ***
## obyv_01 -0.26908 0.08407 -3.201 0.00160 **
## obyv_06 0.60113 0.22660 2.653 0.00864 **
## obyv_09 -0.32490 0.15148 -2.145 0.03320 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 32.34 on 196 degrees of freedom
## Multiple R-squared: 0.8732, Adjusted R-squared: 0.8712
## F-statistic: 449.7 on 3 and 196 DF, p-value: < 2.2e-16resettest(model_lambda)##
## RESET test
##
## data: model_lambda
## RESET = 857.29, df1 = 2, df2 = 194, p-value < 2.2e-16Pri porovnávaní s pôvodným modelom sledujeme zmenu \(R^2_{adj}\) a výsledok RESET testu. Ak síce dôjde k malej zmene \(R^2_{adj}\), ale výrazne sa skomplikuje interpretácia (už nepracujeme priamo s počtom obyvateľov, ale s transformáciou), v praktických aplikáciách sa často preferuje jednoduchší model, pokiaľ štatisticky „stačí“.
9 Cvičenie 8: Zhluková analýza (Cluster Analysis)
9.1 Úvod
V tejto práci predstavíme zhlukovú analýzu pri analýze údajov o počte obyvateľov okresov Slovenskej republiky z hľadiska ich mesačného vývoja v roku 2024. Pracujeme s databázou, ktorá obsahuje okres, obec a počet obyvateľov v jednotlivých mesiacoch roku 2024 (od januára po december). Údaje najskôr agregujeme na úroveň okresov a následne náhodne vyberieme 10 okresov, na ktorých demonštrujeme celý postup hierarchického zhlukovania. V Tab. 1 uvádzame celú nami používanú databázu 10 okresov.
library(knitr)
library(kableExtra)
library(dplyr)# rm(list = ls()) # vypnuté v spoločnom reporte
# Načítanie údajov z databázy
# Predpoklad: súbor "Databaza.csv" je v rovnakom priečinku ako tento Rmd súbor.
# Používame separátor ";" keďže ide o CSV v európskom formáte.
udaje_raw <- read.csv("Databaza.csv",
sep = ";",
stringsAsFactors = FALSE)
# R automaticky dopĺňa X pred názvy stĺpcov začínajúcich číslom,
# preto ho odstránime, aby sme mali späť názvy 2024M01, 2024M02, ...
names(udaje_raw) <- sub("^X", "", names(udaje_raw))
# Premenovanie mesiacov na bežné názvy
mesiace_map <- c(
"2024M01" = "Január",
"2024M02" = "Február",
"2024M03" = "Marec",
"2024M04" = "Apríl",
"2024M05" = "Máj",
"2024M06" = "Jún",
"2024M07" = "Júl",
"2024M08" = "August",
"2024M09" = "September",
"2024M10" = "Október",
"2024M11" = "November",
"2024M12" = "December"
)
# premenujeme len stĺpce s mesiacmi
names(udaje_raw)[match(names(mesiace_map), names(udaje_raw))] <- mesiace_map
month_names <- unname(mesiace_map) # v poradí január–december
# Agregácia na úroveň okresov (súčet obyvateľov za všetky obce)
udaje_okresy <- udaje_raw %>%
group_by(Okres) %>%
summarise(across(all_of(month_names), sum, na.rm = TRUE))
# Náhodný výber 10 okresov
set.seed(123) # pre reprodukovateľnosť
okresy_vyber <- sample(unique(udaje_okresy$Okres), 10)
udaje_vyber <- subset(udaje_okresy, Okres %in% okresy_vyber)
# Usporiadanie podľa názvu okresu (nie je nutné, len pre prehľadnosť)
udaje_vyber <- udaje_vyber[order(udaje_vyber$Okres), ]Tab. 1. Vybrané okresy a ich počet obyvateľov v jednotlivých mesiacoch roku 2024
kable(udaje_vyber, caption = "Vybrané okresy a ich počet obyvateľov (Január–December 2024)") %>%
kable_styling(full_width = FALSE)| Okres | Január | Február | Marec | Apríl | Máj | Jún | Júl | August | September | Október | November | December |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Okres Bratislava I | 47375 | 47391 | 47408 | 47411 | 47484 | 47503 | 47488 | 47525 | 47564 | 47565 | 47632 | 47634 |
| Okres Bratislava II | 125980 | 125975 | 126022 | 126047 | 126178 | 126245 | 126393 | 126435 | 126457 | 126481 | 126522 | 126632 |
| Okres Bratislava III | 77594 | 77678 | 77718 | 77736 | 77746 | 77817 | 77872 | 77912 | 77936 | 77918 | 77998 | 78070 |
| Okres Bratislava IV | 105043 | 105011 | 104992 | 105029 | 105039 | 105064 | 105066 | 105076 | 105140 | 105120 | 105165 | 105118 |
| Okres Bratislava V | 122048 | 122009 | 121976 | 121993 | 121936 | 121934 | 121918 | 121871 | 121869 | 121864 | 121836 | 121885 |
| Okres Dunajská Streda | 127025 | 127173 | 127217 | 127264 | 127307 | 127397 | 127460 | 127512 | 127576 | 127692 | 127772 | 127803 |
| Okres Hlohovec | 27236 | 27182 | 27133 | 27125 | 27084 | 27056 | 27027 | 27004 | 26961 | 26921 | 26904 | 26876 |
| Okres Malacky | 79689 | 79717 | 79761 | 79782 | 79849 | 79899 | 79923 | 79930 | 79917 | 79935 | 79961 | 79994 |
| Okres Pezinok | 69953 | 69936 | 69926 | 69945 | 69972 | 70006 | 70021 | 70070 | 70055 | 70067 | 70099 | 70105 |
| Okres Senec | 105075 | 105281 | 105448 | 105640 | 105848 | 106062 | 106188 | 106382 | 106488 | 106649 | 106781 | 106872 |
Hierarchická zhluková analýza pracuje s mierami vzdialenosti medzi pozorovaniami. Aby boli tieto vzdialenosti porovnateľné, je potrebné, aby všetky premenné boli definované na rovnakej škále. V našom prípade sú premennými mesačné počty obyvateľov (Január až December). Použijeme tzv. z-škálovanie, pričom transformované \(z\) hodnoty (skóre) vypočítame nasledovne
\[z = \frac{x-\mu}{\sigma}\]
kde \(\mu\) je stredná hodnota a \(\sigma\) je štandardná odchýlka pozorovaní \(x\). Predpokladáme pritom, že súbor údajov už neobsahuje NA hodnoty, ktoré boli ošetrené v predchádzajúcich krokoch.
Touto operáciou získame škálované pozorovania, pričom ich rozloženie je znázornené nasledovne.
# =======================================================
## 1) Príprava údajov a data.frame so škálovanými údajmi
## ======================================================
# Vyberieme len numerické stĺpce s mesačnými údajmi (Január–December)
udaje10 <- udaje_vyber[, month_names]
# Názvy riadkov nastavíme na názvy okresov
rownames(udaje10) <- udaje_vyber$Okres
# Odstránime prípadné chýbajúce hodnoty
udaje_complete <- na.omit(udaje10)
# Škálovanie (z-škóre) jednotlivých premenných
udaje_scaled <- scale(udaje_complete)Obr. 1. Boxploty škálovaných numerických premenných (mesiace roku 2024)
num_vars <- as.data.frame(udaje_scaled)
num_plots <- ncol(num_vars)
# 12 mesiacov = 3x4 mriežka
par(mfrow = c(3, 4),
mar = c(4, 4, 2, 1),
oma = c(0, 0, 3, 0))
for (col in names(num_vars)) {
boxplot(num_vars[[col]],
main = col,
col = "lightblue",
horizontal = TRUE)
}
mtext("Boxploty numerických premenných (mesiace roku 2024)",
outer = TRUE, cex = 1.3, font = 2)
Tentokrát odľahlé hodnoty nevylučujeme, nakoľko definujú konkrétny okres. Napríklad okres s výrazne vyšším počtom obyvateľov je z pohľadu analýzy zaujímavý a chceme ho v dátach zachovať.
Pri zhlukovej analýze je dôležitá korelačná matica premenných. Vysoká korelácia zvýhodňuje pri zhlukovej analýze korelované premenné. Preto pri korelácii nad 0,8 alebo 0,9 často vylúčime jednu z korelovaných premenných. V Tab. 2 ukážeme korelačnú maticu medzi jednotlivými mesiacmi roku 2024. Ak by sa medzi mesiacmi vyskytovali veľmi vysoké korelácie, bolo by možné uvažovať napríklad o znížení dimenzie pomocou analýzy hlavných komponentov.
V prípade, ak máme väčší počet významne korelovaných premenných, sa odporúča transformácia pomocou Analýzy hlavných komponentov (Principal Component Analysis).
Tab. 2. Korelačná matica medzi mesiacmi roku 2024
cor_mat <- cor(udaje_scaled, use = "pairwise.complete.obs")
cor_mat <- round(cor_mat, 2)
print(cor_mat)## Január Február Marec Apríl Máj Jún Júl August September Október
## Január 1 1 1 1 1 1 1 1 1 1
## Február 1 1 1 1 1 1 1 1 1 1
## Marec 1 1 1 1 1 1 1 1 1 1
## Apríl 1 1 1 1 1 1 1 1 1 1
## Máj 1 1 1 1 1 1 1 1 1 1
## Jún 1 1 1 1 1 1 1 1 1 1
## Júl 1 1 1 1 1 1 1 1 1 1
## August 1 1 1 1 1 1 1 1 1 1
## September 1 1 1 1 1 1 1 1 1 1
## Október 1 1 1 1 1 1 1 1 1 1
## November 1 1 1 1 1 1 1 1 1 1
## December 1 1 1 1 1 1 1 1 1 1
## November December
## Január 1 1
## Február 1 1
## Marec 1 1
## Apríl 1 1
## Máj 1 1
## Jún 1 1
## Júl 1 1
## August 1 1
## September 1 1
## Október 1 1
## November 1 1
## December 1 1Každému okresu zodpovedá jeden riadok pozorovaní. Vzdialenosť medzi okresmi \(i\) a \(j\) je:
\[ d^{ij} = \sqrt{\sum_k (x^i_k - x^j_k)^2} \]
kde \(x^i_k\) je \(k\)-tá premenná vstupujúca do výpočtu (v našom prípade ide o mesačné počty obyvateľov: Január, Február, …, December) okresu \(i\). Tento typ vzdialenosti nazývame Euklidovská vzdialenosť. Vzdialenosti medzi jednotlivými okresmi sa súhrnne vyjadrujú aj v matici vzdialenosti, čo je v našom prípade uvedené v Tab. 3. Analýzou tejto tabuľky môžeme identifikovať dvojice okresov, ktoré sú si podľa priebehu počtu obyvateľov v roku 2024 veľmi podobné (malá vzdialenosť), resp. veľmi odlišné (veľká vzdialenosť).
Tab. 3. Euklidovská matica vzdialeností medzi okresmi
## ============================
## 3) Distance matrix
## ============================
# ponecháme pôvodné názvy okresov ako rownames
dist_mat <- round(dist(udaje_scaled, method = "euclidean"), 2)
dist_mat## 1 2 3 4 5 6 7 8 9
## 2 7.94
## 3 3.06 4.88
## 4 5.80 2.14 2.75
## 5 7.50 0.44 4.44 1.70
## 6 8.06 0.12 5.00 2.25 0.56
## 7 2.06 10.00 5.12 7.86 9.56 10.12
## 8 3.26 4.68 0.20 2.54 4.24 4.79 5.32
## 9 2.27 5.67 0.79 3.53 5.23 5.79 4.33 0.99
## 10 5.90 2.04 2.85 0.11 1.60 2.15 7.96 2.64 3.639.2 Princíp hierarchického zhlukovania (Wardova metóda)
Zhlukovanie v prípade Wardovej metódy prebieha zdola smerom nahor, t. j. začíname s jednočlennými klastrami, ktoré postupne zlučujeme. Táto metóda patrí medzi aglomeratívne hierarchické metódy. Minimalizuje nárast vnútornej variability pri spojení dvoch klastrov, pričom využíva nasledovné výpočty:
Wardova metóda minimalizuje sumu štvorcov chýb (Error Sum of Squares – ESS)
\[ESS(C) = \sum_{i \in C} \lVert x_i - \bar{x}_C \rVert^2,\]
kde \(C\) je zvažovaný klaster (zhluk). V každom kroku zlučovania dvoch klastrov Wardova metóda hľadá minimálny prírastok sumy štvorcov chýb (\(\Delta ESS\)), pričom
\[\Delta ESS = ESS(A \cup B) - ESS(A) - ESS(B).\]
Dvojica zhlukov, ktorá tejto podmienke o minimalizácii vyhovuje, je následne zlúčená a prechádza sa k ďalšiemu kroku. To spravidla vedie k vytváraniu relatívne homogénnych zhlukov, pričom nedochádza k odtrhávaniu odľahlých hodnôt tak, ako pri niektorých iných zhlukovacích metódach.
Obr. 2. Hierarchické zhlukovanie – dendrogram. Červená čiara určuje rez definujúci tri klastre.
## ============================
## 4) Hierarchical clustering
## ============================
hc <- hclust(dist_mat, method = "ward.D2")
plot(hc, labels = rownames(udaje_scaled),
main = "Hierarchical clustering of districts (Ward.D2)",
xlab = "", sub = "")
# počet klastrov
k <- 3
h_cut <- hc$height[length(hc$height) - (k - 1)]
abline(h = h_cut, col = "red", lwd = 2, lty = 2)
klaster_membership <- cutree(hc, k = k)
udaje_klasters <- data.frame(
Okres = rownames(udaje_complete),
udaje_complete,
klaster = factor(klaster_membership)
)Tab. 4. Príslušnosť okresov do klastrov
data_prac <- data.frame(
Okres = udaje_klasters$Okres,
klaster = udaje_klasters$klaster
)
data_prac9.3 Deskriptívne štatistiky výsledkov
Na základe Tab. 5 môžeme posúdiť, ako dobre
jednotlivé mesiace (premenné) prispievajú k rozlišovaniu medzi
klastrami. Vnútroklastrová variabilita (WSS) vyjadruje, nakoľko sú
okresy v rámci jedného klastra podobné, zatiaľ čo medzi-klastrová
variabilita (BSS) vyjadruje rozdiely medzi klastrami. Podiel
medzi-klastrovej variability na celkovej variabilite (stĺpec
Prop_Between) nám ukazuje, do akej miery daný mesiac
prispieva k separácii klastrov.
Tab. 5. Vysvetlenie vnútroklastrovej variability z hľadiska jednotlivých premenných (mesiacov)
## ============================
## 5) Variability measures
## ============================
ssq <- function(x, m) sum((x - m)^2)
var_names <- colnames(udaje_scaled)
TSS <- sapply(var_names, function(v) ssq(udaje_scaled[, v], mean(udaje_scaled[, v])))
WSS <- sapply(var_names, function(v) {
x <- udaje_scaled[, v]
tapply(x, klaster_membership, function(z) ssq(z, mean(z))) |> sum()
})
BSS <- TSS - WSS
ss_table <- data.frame(
Variable = var_names,
TSS = TSS,
WSS = WSS,
BSS = BSS,
Prop_Between = BSS / TSS
)
ss_tableVyššie hodnoty podielu Prop_Between znamenajú, že daný
mesiac výraznejšie odlišuje jednotlivé klastre. Mesiace s nižšou
hodnotou tohto podielu predstavujú obdobia, v ktorých sú okresy medzi
sebou relatívne podobné a daná premenná nie je takým silným
„separátorom“.
Na záver sa pozrieme na centroidy, t. j. priemerné hodnoty sledovaných premenných v jednotlivých klastroch. Tie nám poskytujú prehľadný opis typického okresu v každom klastri.
# Vytvoríme dátový rámec s pôvodnými (neškálovanými) hodnotami a priradeným klastrom
udaje_complete_klaster <- data.frame(
Okres = rownames(udaje_complete),
udaje_complete,
klaster = factor(klaster_membership)
)Tab. 6. Centroidy – priemerné počty obyvateľov v jednotlivých mesiacoch podľa klastrov
descriptives <- udaje_complete_klaster %>%
group_by(klaster) %>%
summarise(
across(
.cols = where(is.numeric),
.fns = list(mean = ~ mean(.x, na.rm = TRUE)),
.names = "{.col}_{.fn}"
)
)
descriptives9.4 Záver
Predložená analýza sa zaoberá počtom obyvateľov vybraných okresov Slovenskej republiky v závislosti od ich mesačného vývoja v roku 2024. Na základe hierarchickej zhlukovej analýzy sme okresy rozdelili do troch klastrov. Klastre sa líšia najmä úrovňou priemerného počtu obyvateľov a podobnosťou priebehu v čase.
Takáto analýza môže slúžiť ako podklad pre:
- plánovanie verejných služieb (školy, zdravotníctво, doprava),
- priestorové plánovanie a rozvoj územia,
- porovnávanie typicky „malých“, „stredných“ a „väčších“ okresov podľa zvolených ukazovateľov.
Zhluková analýza teda umožňuje identifikovať skupiny podobných okresov a následne cieleným spôsobom analyzovať a plánovať verejné politiky či rozvojové projekty podľa ich príslušnosti ku klastrom.
10 Cvičenie 9: Econometrics in R - cvičenie 9
S využitím databázy Databaza.csv, ktorá obsahuje počty obyvateľov slovenských obcí podľa okresu a mesiaca v roku 2024.
Pri ďalšej práci budeme používať knižnice
library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
library(dplyr)
library(tidyr)
# rm(list = ls()) # vypnuté v spoločnom reporte10.1 1. Úvod a údaje
Údaje o počte obyvateľov sú uložené v súbore
Databaza.csv.
Každý riadok predstavuje jednu obec, stĺpce:
Okres– názov okresu,Obec– názov obce,2024M01…2024M12– počet obyvateľov v danom mesiaci roku 2024.
Súbor predpokladáme uložený v podpriečinku
udaje v pracovnom priečinku (aby boli
údaje oddelené od zvyšku projektu).
10.1.1 1.1 Úvod do problému, stanovenie hypotéz
Rozhodli sme sa sledovať vývoj počtu obyvateľov v
čase vo vybranej obci.
Konkrétne budeme modelovať mesačný počet obyvateľov v
obci
„Bratislava - mestská časť Staré Mesto“ v závislosti od
času (index mesiaca).
Označme:
- \(Y_t\): počet obyvateľov v mesiaci \(t\),
- \(t\): časový index (1 = január 2024, 2 = február 2024, …, 12 = december 2024).
Naša pracovná hypotéza:
- predpokladáme, že počet obyvateľov sa v priebehu roka mierne mení (napr. vplyvom migrácie),
- očakávame hladký vývoj v čase, teda silnú autokoreláciu, t. j. počet obyvateľov v mesiaci \(t\) je veľmi podobný počtu v mesiaci \(t-1\),
- očakávame kladný trendový koeficient (\(\beta_1 > 0\)), ak počet obyvateľov má mierne rastúcu tendenciu.
10.1.2 1.2 Príprava databázy, úprava údajov
Databáza je v tzv. „wide“ tvare (mesiace v stĺpcoch). Pre časovú analýzu jednej obce je výhodnejšie mať údaje v „long“ tvare: každý riadok = jedna obec v jednom mesiaci.
Najskôr načítame údaje, potom vyberieme jednu obec a preklopíme mesačné stĺpce do časového radu.
# načítanie databázy
udaje <- read.csv("Databaza.csv",
sep = ";",
dec = ",",
header = TRUE,
check.names = FALSE)
# vyberieme jednu konkrétnu obec (možeš zmeniť na inú)
udajeMesto <- udaje[udaje$Obec == "Bratislava - mestská časť Staré Mesto", ]
# preklopenie z wide na long (mesiace do riadkov)
udajeMesto_long <- udajeMesto %>%
pivot_longer(
cols = dplyr::starts_with("2024M"),
names_to = "RokMesiac",
values_to = "Pocet_obyvatelov"
) %>%
arrange(RokMesiac) %>%
mutate(
Rok = as.numeric(substr(RokMesiac, 1, 4)),
Mes = as.numeric(substr(RokMesiac, 6, 7)),
Time = dplyr::row_number() # 1 = január, 12 = december
)
udajeMesto_long10.2 2. Lineárna regresia v základnom tvare
Ide o odhad rovnice
\[ \text{Pocet\_obyvatelov}_t = \beta_0 + \beta_1 \text{Time}_t + e_t, \]
kde:
- \(\text{Pocet\_obyvatelov}_t\) – počet obyvateľov vo vybranom mesiaci,
- \(\text{Time}_t\) – index mesiaca (1–12),
- \(e_t\) – náhodná zložka (rezíduum).
library(ggplot2)
# lineárny trendový model
model <- lm(Pocet_obyvatelov ~ Time, data = udajeMesto_long)
summary(model)##
## Call:
## lm(formula = Pocet_obyvatelov ~ Time, data = udajeMesto_long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.232 -12.214 1.946 12.539 23.684
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 47339.47 10.96 4318.14 < 2e-16 ***
## Time 24.44 1.49 16.41 1.47e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 17.81 on 10 degrees of freedom
## Multiple R-squared: 0.9642, Adjusted R-squared: 0.9606
## F-statistic: 269.2 on 1 and 10 DF, p-value: 1.472e-0810.3 3. Autokorelácia rezíduí
V tejto časti sa pozrieme na ďalší dôležitý predpoklad klasického lineárneho regresného modelu – nezávislosť rezíduí. V časových radoch (ako je náš mesačný počet obyvateľov) sa často stáva, že chyba v čase \(t\) je systematicky spätá s chybou v čase \(t-1\), čo nazývame autokoreláciou rezíduí.
10.3.1 3.1 Čo je autokorelácia rezíduí?
Autokorelácia rezíduí je situácia, keď platí:
\[ \operatorname{Corr}(e_t, e_{t-k}) \neq 0 \quad \text{pre niektoré } k \neq 0. \]
Najčastejšie sa skúma autokorelácia prvého rádu:
\[ e_t = \rho e_{t-1} + \nu_t,\quad |\rho| < 1. \]
10.3.2 3.2 Dôsledky autokorelácie rezíduí
- odhady koeficientov \(\hat{\beta}\) sú nestranné,
- ale neefektívne,
- štandardné chyby sú skreslené (často podhodnotené), teda p-hodnoty sa javia menšie,
- t- a F-testy sú potom skreslené.
10.3.3 3.3 Detekcia autokorelácie
10.3.3.1 Grafická informácia
# pridáme do dát fitted values z modelu
udajeMesto_long$fitted <- fitted(model)
# scatterplot + vyrovnaná trendová čiara
ggplot(udajeMesto_long, aes(x = Time, y = Pocet_obyvatelov)) +
geom_point(color = "steelblue", size = 2) +
geom_line(aes(y = fitted), color = "red", size = 1) +
labs(
title = "Počet obyvateľov: empirické údaje (modrá) vs. odhadnutý trend (červená)",
x = "Mesiac (2024)",
y = "Počet obyvateľov"
) +
theme_minimal()
Analýzou obrázka vidíme, že mesačné hodnoty počtu obyvateľov na seba plynulo nadväzujú – v susedných mesiacoch sú si veľmi podobné. Ak by sme si zobrazili rezíduá (rozdiel medzi empirickými a vyrovnanými hodnotami), je prirodzené očakávať, že počas niekoľkých po sebe idúcich mesiacov budú mať podobné znamienko. To je typický vizuálny signál možnej autokorelácie.
# uložíme si reziduá z pôvodného modelu
res <- residuals(model)10.3.3.2 ACF graf (Autocorrelation Function)
Táto funkcia priraďuje odhad korelácie, ktorá je medzi jednotlivými rezíduami v aktuálnom období a období posunutom (Lag) o \(k\) období späť.
acf(res, lag.max = 4, main = "Autokorelačná funkcia rezíduí")
Na tomto grafe je modrou prerušovanou čiarou vyjadrený aj 95 % interval spoľahlivosti pre hodnotu autokorelačného koeficientu s príslušným posunom. Ak sa niektorý stĺpec nachádza mimo tohto pásma, naznačuje to štatisticky významnú autokoreláciu pri danom posune.
10.3.3.3 Durbin–Watsonov test
Durbin–Watsonov koeficient (DW) je vypočítaný z rezíduí podľa vzorca
\[ DW = \frac{\sum_{t=2}^{n} (e_t - e_{t-1})^{2}}{\sum_{t=1}^{n} e_t^{2}} \]
kde \(n\) je počet pozorovaní. Medzi koeficientom autokorelácie dvoch susedných rezíduí a DW platí približný vzťah
\[ \hat{\rho} \approx 1 - \frac{DW}{2}. \]
Hodnoty:
- blízke nule → silná pozitívna autokorelácia,
- blízke 4 → silná negatívna autokorelácia,
- okolo 2 → žiadny výrazný problém s autokoreláciou.
V praxi sa často používa intuitívne pravidlo: ak sa DW nachádza v intervale približne 1,8 až 2,2, problém autokorelácie väčšinou nepovažujeme za vážny.
library(lmtest)
dwtest(model)##
## Durbin-Watson test
##
## data: model
## DW = 2.4468, p-value = 0.6766
## alternative hypothesis: true autocorrelation is greater than 0DW test má isté obmedzenia – regresory nesmú byť časovo posunuté a nesmú obsahovať oneskorené pozorovania vysvetľovanej veličiny ako regresory. Tieto obmedzenia neplatia pre Breusch–Godfreyov test.
10.3.3.4 Breusch–Godfreyov test (BG test)
10.3.3.4.1 Čo test testuje
BG test je formálnym testom autokorelácie (sériovej korelácie) rezíduí:
\[ u_t = \rho_1 u_{t-1} + \rho_2 u_{t-2} + \cdots + \rho_p u_{t-p} + \varepsilon_t. \]
Testuje, či sú rezíduá korelované v čase až do rádu \(p\), čo je typické pre časové rady.
10.3.3.4.2 Hypotézy
Pre test autokorelácie až do rádu \(p\):
Nulová hypotéza \(H_0\): žiadna sériová korelácia
\[ \rho_1 = \rho_2 = \cdots = \rho_p = 0 \]Alternatívna hypotéza \(H_1\): séria je autokorelovaná
Aspoň jedna z \(\rho_j\) je odlišná od nuly.
10.3.3.4.3 Ako funguje BG test
- Odhadneme pôvodnú regresiu (u nás trendový model) a získame rezíduá \(e_t\).
- Spustíme pomocnú regresiu: \[ e_t = \alpha_0 + \alpha_1 x_{1t} + \cdots + \alpha_k x_{kt} + \rho_1 e_{t-1} + \cdots + \rho_p e_{t-p} + v_t, \] kde \(x_{jt}\) sú pôvodné regresory.
- Z tejto regresie vypočítame: \[ \text{BG} = (n - p) R^2_{\text{aux}} \]
- Pod \(H_0\) platí približne: \[ \text{BG} \sim \chi^2_p. \]
10.3.3.4.4 Interpretácia
- Veľká hodnota testovej štatistiky BG (malá p-hodnota) → zamietame \(H_0\) → autokorelácia prítomná.
- Malá hodnota BG (veľká p-hodnota) → nezamietneme \(H_0\) → nie je dôkaz autokorelácie.
10.3.3.4.5 Praktický výpočet v R
bgtest(model, order = 1)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model
## LM test = 0.70926, df = 1, p-value = 0.3997V našich údajoch (mesačný počet obyvateľov vo vybranej obci v roku 2024) sa ukáže, či BG test potvrdí alebo vyvráti prítomnosť autokorelácie prvého rádu. Aj keď pri krátkom časovom rade (len 12 mesiacov) nebýva test príliš silný, pre demonštračné účely použijeme aj postup odstránenia dôsledkov autokorelácie pomocou dynamizácie modelu.
10.4 Ako riešiť autokoreláciu
10.4.1 Koyckova transformácia a Koyckova rovnica
V predchádzajúcej časti sme sa venovali autokorelácii rezíduí. Teraz rozšírime model o dynamickú štruktúru založenú na tzv. geometricky klesajúcom rozložení oneskorených efektov. Takýto model sa nazýva Koyckov model a jeho ústredným prvkom je Koyckova transformácia.
10.4.1.0.1 Východiskový model s distribuovaným oneskorením
Uvažujme model:
\[ y_t = \alpha + \beta_0 x_t + \beta_1 x_{t-1} + \beta_2 x_{t-2} + \cdots + u_t. \]
10.4.1.0.2 Koyckova štruktúra koeficientov
Koyck navrhol, že oneskorené efekty majú geometricky klesajúcu podobu:
\[ \beta_k = \lambda^k \beta_0, \qquad 0 < \lambda < 1. \]
Tým sa distribučné oneskorenie zjednoduší tak, že namiesto nekonečného počtu parametrov odhadujeme len:
- \(\beta_0\) – okamžitý efekt,
- \(\lambda\) – rýchlosť tlmenia efektov oneskorenia.
10.4.1.0.3 Koyckova transformácia
Pôvodný model:
\[ y_t = \alpha + \beta_0 x_t + \lambda \beta_0 x_{t-1} + \lambda^2 \beta_0 x_{t-2} + \cdots + u_t. \]
Vynásobíme model konštantou \(\lambda\) a posunieme časový index o 1 dozadu:
\[ \lambda y_{t-1} = \lambda\alpha + \beta_0 x_{t-1} + \lambda \beta_0 x_{t-2} + \cdots + \lambda u_{t-1}. \]
Odčítaním oboch výrazov dostaneme:
\[ y_t - \lambda y_{t-1} = \alpha(1-\lambda) + \beta_0 (x_t - \lambda x_{t-1}) + (u_t - \lambda u_{t-1}). \]
To je Koyckova transformácia.
10.4.1.0.4 Koyckova rovnica (autoregresívny model)
Po algebraickej úprave:
\[ y_t = \alpha(1-\lambda) + \beta_0 x_t + \lambda y_{t-1} + v_t, \]
kde:
\[ v_t = u_t - \lambda u_{t-1}. \]
Toto je Koyckova rovnica – dynamický autoregresívny model so závislosťou na minulosti.
V našom kontexte:
- \(y_t\) = počet obyvateľov v mesiaci \(t\),
- \(x_t\) – napríklad časový trend (Time),
- \(y_{t-1}\) – počet obyvateľov v predchádzajúcom mesiaci.
10.4.2 Odstraňovanie problému autokorelácie rezíduí
10.4.2.1 Odhad Koyckovho modelu v R
Najjednoduchšia implementácia pre naše údaje:
udajeMesto_long <- udajeMesto_long %>%
arrange(Time) %>%
mutate(
Pocet_obyvatelov_lag1 = dplyr::lag(Pocet_obyvatelov)
)
model_koyck <- lm(Pocet_obyvatelov ~ Time + Pocet_obyvatelov_lag1,
data = udajeMesto_long)
summary(model_koyck)##
## Call:
## lm(formula = Pocet_obyvatelov ~ Time + Pocet_obyvatelov_lag1,
## data = udajeMesto_long)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.591 -16.219 2.003 13.554 23.918
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 58792.7569 15755.3637 3.732 0.00577 **
## Time 30.9479 8.3103 3.724 0.00584 **
## Pocet_obyvatelov_lag1 -0.2422 0.3330 -0.727 0.48779
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.75 on 8 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.9609, Adjusted R-squared: 0.9512
## F-statistic: 98.4 on 2 and 8 DF, p-value: 2.329e-06dwtest(model_koyck)##
## Durbin-Watson test
##
## data: model_koyck
## DW = 2.2885, p-value = 0.5587
## alternative hypothesis: true autocorrelation is greater than 0V dynamizovanom modeli vystupuje aj oneskorená hodnota vysvetľovanej
premennej (Pocet_obyvatelov_lag1). Ak je jej koeficient
kladný a menší ako 1, interpretujeme to tak, že
časť úrovne z minulého mesiaca sa prenáša do aktuálneho
mesiaca (zotrvačnosť). Pri hodnotení kvality modelu je vhodné
porovnať napríklad Adjusted R-squared pôvodného a
dynamického modelu, prípadne znova overiť autokoreláciu rezíduí (DW, BG
test).
10.4.2.2 Newey–West robustné štandardné chyby
Ešte jednou možnosťou, ako reagovať na prítomnosť autokorelácie (a heteroskedasticity), je použiť robustné štandardné chyby podľa Newey–Westa. Tie nemenia samotné odhady koeficientov, ale korigujú odhady smerodajných odchýlok.
library(sandwich)
library(lmtest)
coeftest(model, vcov = NeweyWest(model))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.7339e+04 6.5101e+00 7271.681 < 2.2e-16 ***
## Time 2.4441e+01 7.7275e-01 31.628 2.347e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 110.5 4. Záver
Na jednoduchom príklade mesačného počtu obyvateľov vo vybranej slovenskej obci v roku 2024 sme ukázali:
- ako odhadnúť trendový model v čase,
- ako graficky a štatisticky testovať autokoreláciu rezíduí (ACF, Durbin–Watson, BG test),
- ako možno model dynamizovať pomocou Koyckovej rovnice,
- a ako využiť Newey–Westove robustné štandardné chyby.
Autokorelácia rezíduí a z nej vyplývajúca potreba dynamizácie modelov má v ekonometrii veľký význam. Tu sme uviedli len základné prístupy – k zložitejším patrí napríklad Almonov model distribuovaného oneskorenia, metóda Cochran–Orcutt alebo všeobecnejšie ARIMA modely, ktoré sa používajú pri dlhších časových radoch.
11 Cvičenie 10: Cvičenie 10 – Multikolinearita v regresných modeloch
12 1. Úvod
Po autokorelácii a heteroskedasticite rezíduí je multikolinearita tretím závažným porušením predpokladov použitia metódy najmenších štvorcov. Tu sa okrem iného predpokladá, že matica \(\mathbf X\) je tvorená lineárne nezávislými riadkami a tiež stĺpcami, čo zabezpečí regularitu matice \(\mathbf X^T\mathbf X\) a teda možnosť jej inverzie. Tá sa používa pri odhadoch regresných koeficientov. V praxi sa ale môže stať, že vzniká „takmer singulárna“ matica \(\mathbf X^T\mathbf X\), t. j. matica \(\mathbf X\) je tvorená „približne“ lineárne závislými stĺpcami, t. j. existuje taká ich lineárna kombinácia, v ktorej
\[ x_{il} = \alpha_0 + \alpha_1 x_{i1} + \dots + \alpha_{l-1} x_{i,(l-1)} + \alpha_{l+1} x_{i,(l+1)} + \alpha_k x_{i,k} + \nu_i \]
Tu \(\nu_i\) sú rádovo menšie čísla než regresory \(x_{.}\), t. j. \((\forall i)(\nu_i << x_{.,i})\). V tomto prípade je inverzná matica \((\mathbf X^T\mathbf X)^{-1}\) veľmi nestabilná a obsahuje na hlavnej diagonále veľmi veľké hodnoty. Táto matica sa používa pri výpočtoch \[ \hat \beta = (\mathbf X^T\mathbf X)^{-1} \mathbf X^T \mathbf y \] a tiež \[ \text{std}(\beta_i) = \sqrt{\sigma^2 (\mathbf X^T \mathbf X)^{-1}_{ii}}. \] To spôsobuje nestabilitu odhadovaných regresných koeficientov a ich nadhodnotené rozptyly.
Tento problém nazývame problémom multikolinearity.
13 2. Dôsledky multikolinearity
Multikolinearita patrí medzi najčastejšie problémy viacnásobnej
lineárnej regresie.
Je dôležité jasne rozlišovať dva fakty:
- Nespôsobuje skreslené (biased) odhady koeficientov.
- Nadhodnocuje odhady štandardných odchýlok regresných koeficientov a vedie potom k falošnému neprijímaniu alternatívnej hypotézy o štatistickej významnosti jednotlivých regresorov.
- Odhadované regresné koeficienty sú nestabilné – pri malej zmene údajov sa prudko menia koeficienty ako aj ich znamienka.
- Interpretácia regresného modelu je z dôvodu vyššie uvedených dôvodov nespoľahlivá.
14 3. Východiskový model a údaje
Budeme pracovať s prierezovým regresným modelom, kde pre každú obec uvažujeme počet obyvateľov v decembri 2024 ako závislú premennú a počty obyvateľov v predchádzajúcich mesiacoch toho istého roku ako vysvetľujúce premenné:
\[ \text{Pop12}_i = \beta_0 + \beta_1 \text{Pop11}_i + \beta_2 \text{Pop10}_i + \beta_3 \text{Pop09}_i + u_i, \]
kde:
- \(\text{Pop12}_i\) – počet
obyvateľov obce \(i\) v decembri
2024,
- \(\text{Pop11}_i\) – počet
obyvateľov v novembri 2024,
- \(\text{Pop10}_i\) – počet
obyvateľov v októbri 2024,
- \(\text{Pop09}_i\) – počet obyvateľov v septembri 2024.
Údaje máme z databázy počtu obyvateľov slovenských obcí podľa mesiacov (rok 2024). Po načítaní ich uložíme do data.frame udaje a zároveň premenujeme stĺpce na kratšie názvy.
udaje <- read.csv("Databaza.csv",
sep = ";",
dec = ".",
header = TRUE,
stringsAsFactors = FALSE)
names(udaje)## [1] "Okres" "Obec" "X2024M12" "X2024M11" "X2024M10" "X2024M09"
## [7] "X2024M08" "X2024M07" "X2024M06" "X2024M05" "X2024M04" "X2024M03"
## [13] "X2024M02" "X2024M01"# vyberieme si len stĺpce, ktoré potrebujeme
udaje <- udaje[, c("Okres", "Obec", "X2024M12", "X2024M11", "X2024M10", "X2024M09")]
# premenujeme stĺpce na jednoduchšie názvy
colnames(udaje) <- c("Okres", "Obec", "Pop12", "Pop11", "Pop10", "Pop09")
head(udaje)V našej databáze nemáme chýbajúce hodnoty, takže nie je potrebná žiadna imputácia.
15 4. Odhad základného regresného modelu
model <- lm(Pop12 ~ Pop11 + Pop10 + Pop09,
data = udaje)
summary(model)##
## Call:
## lm(formula = Pop12 ~ Pop11 + Pop10 + Pop09, data = udaje)
##
## Residuals:
## Min 1Q Median 3Q Max
## -44.442 -2.185 0.749 3.659 50.908
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.70441 0.68551 -2.486 0.0137 *
## Pop11 1.24863 0.06192 20.164 < 2e-16 ***
## Pop10 0.09468 0.12167 0.778 0.4374
## Pop09 -0.34283 0.08212 -4.175 4.48e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.856 on 196 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 1.19e+08 on 3 and 196 DF, p-value: < 2.2e-16Vo výsledkoch zvyčajne vidíme veľmi vysoký koeficient determinácie a koeficienty pri Pop11, Pop10 a Pop09 sú si veľmi podobné. To je intuitívne, pretože počty obyvateľov v jednotlivých mesiacoch sa menia len veľmi málo a sú takmer perfektné lineárne kombinácie.
16 5. Korelačná matica
Korelácia dokáže zachytiť párové vzťahy medzi premennými. Ak medzi niektorými vysvetľujúcimi premennými je vysoká korelácia (signalizujúca multikolinearitu), potom je najjednoduchšie ju zo zoznamu regresorov vylúčiť. Korelácie sa dajú aj testovať, alebo len vyčísliť a potom podľa intuitívneho pravidla vylúčiť jednu premennú, ktorá má koreláciu s inou premennou v absolútnej hodnote vyššiu ako 0.8, resp. 0.9.
xvars <- udaje[, c("Pop11", "Pop10", "Pop09")]
round(cor(xvars), 3)## Pop11 Pop10 Pop09
## Pop11 1 1 1
## Pop10 1 1 1
## Pop09 1 1 1V našom prípade uvidíme, že korelácie medzi mesačnými počtami obyvateľov sú prakticky rovné 1. To je extrémny príklad multikolinearity – všetky vysvetľujúce premenné obsahujú takmer tú istú informáciu.
Korelačný vzťah sa dá vytušiť aj z jednoduchých párových scatterplotov ako je to na nasledujúcom obrázku – body ležia takmer na priamke.
pairs(xvars,
main = "Scatterplotová matica – premenné Pop11, Pop10, Pop09")
17 6. VIF
Indikátorom multikolinearity u premennej, ktorá multikolinearitu zapríčiňuje, je Variance Inflation Factor (VIF). Pre premennú \(x_j\) je potom
\[ VIF_j = \frac{1}{1 - R_j^2} \]
kde \(R_j^2\) pochádza z regresie:
\[ X_j = \gamma_0 + \gamma_1 X_1 + \cdots + \gamma_{j-1} X_{j-1} + \gamma_{j+1} X_{j+1} + u. \]
library(car)
vif(model)## Pop11 Pop10 Pop09
## 1367712 5278985 2404844Intuitívnym kritériom, ktoré signalizuje prítomnosť multikolinearity, je podmienka VIF > 5 (prísne kritérium), alebo VIF > 10 (menej prísne kritérium). V našich dátach sú VIF-y typicky veľmi vysoké, čo zodpovedá takmer dokonalej korelácii medzi mesačnými hodnotami počtu obyvateľov.
18 7. Condition Number
Pri existencii multikolinearity sa model prejavuje tak, že koeficient determinácie je síce vysoký a zdá sa, že model je veľmi dobrý, ale regresné koeficienty nemusia byť stabilné a ich štandardné odchýlky môžu byť nafúknuté. Uvedomíme si to, ak sa pozrieme, ako sa počítajú – t. j.
\[ \text{std}(\beta_i) = \sqrt{\sigma^2 (\mathbf X^T \mathbf X)^{-1}_{ii}}, \]
kde rozhodujúci je \(i\)-ty prvok hlavnej diagonály matice \((\mathbf X^T \mathbf X)^{-1}\). Tie prvky sú v prípade podobnosti vysvetľujúcich premenných mimoriadne veľké. Túto situáciu zachytáva nasledovný ukazovateľ.
Pri výpočte Condition number \(\kappa\) sa používa vzorec
\[ \kappa = \frac{\theta_{\text{max}}}{\theta_{\text{min}}} \]
kde \(\theta_.\) sú vlastné čísla matice. Condition number nie je test, je to len indikátor, ktorý posudzuje mieru multikolinearity medzi premennými. Používame intuitívne pravidlo. Ak Condition number je
- < 10 → nízka multikolinearita,
- 10–30 → mierna,
- 30–100 → silná,
100 → veľmi vážna.
V našom prípade to vypočítame nasledovne:
X <- model.matrix(model)[, -1]
XtX <- t(X) %*% X
eig <- eigen(XtX)
condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number## [1] 5315.386Keďže u nás tento indikátor typicky výrazne presahuje hodnotu 100, signalizuje prítomnosť závažnej multikolinearity medzi mesačnými údajmi.
Vlastné číslo štvorcovej matice \(\mathbf X^T \mathbf X\) je číslo \(\theta_j\), pre ktoré platí \((\mathbf X^T \mathbf X)\mathbf h^j = \theta_j\mathbf h^j\). \(\mathbf h^j\) je tzv. vlastný vektor tejto matice. Máme toľko vlastných čísel (teda aj vlastných vektorov), koľko obsahuje matica \(\mathbf X^T \mathbf X\) riadkov (stĺpcov).
Môže sa stať, že VIF faktor nesignalizuje multikolinearitu u žiadnej „jednej“ vysvetľujúcej veličiny, ale sú navzájom prepojené cyklickými lineárnymi závislosťami všetky premenné. To zachytáva práve Condition Number.
19 8. Riešenia multikolinearity
19.1 Vynechanie premennej
Pokúsme sa vynechať postupne dve premenné, ktoré majú najvyšší VIF, a porovnajme následne upravené koeficienty determinácie oboch nových modelov:
model_noPop10 <- lm(Pop12 ~ Pop11 + Pop09, data = udaje)
summary(model_noPop10)##
## Call:
## lm(formula = Pop12 ~ Pop11 + Pop09, data = udaje)
##
## Residuals:
## Min 1Q Median 3Q Max
## -43.171 -2.310 0.847 3.530 50.940
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.68757 0.68448 -2.465 0.0145 *
## Pop11 1.28666 0.03798 33.876 < 2e-16 ***
## Pop09 -0.28619 0.03799 -7.534 1.74e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.847 on 197 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 1.789e+08 on 2 and 197 DF, p-value: < 2.2e-16model_noPop11 <- lm(Pop12 ~ Pop10 + Pop09, data = udaje)
summary(model_noPop11)##
## Call:
## lm(formula = Pop12 ~ Pop10 + Pop09, data = udaje)
##
## Residuals:
## Min 1Q Median 3Q Max
## -53.303 -4.884 -0.385 4.526 85.897
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.8557 1.1967 -0.715 0.475
## Pop10 2.0312 0.1307 15.546 < 2e-16 ***
## Pop09 -1.0305 0.1307 -7.888 2.07e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.49 on 197 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 5.836e+07 on 2 and 197 DF, p-value: < 2.2e-16V našej databáze zvyčajne uvidíme, že vynechanie jednej z veľmi podobných mesačných premenných zníži upravený koeficient determinácie len minimálne. To je typické správanie pri silnej multikolinearite – jedna premenná dokáže takmer dokonale „zastúpiť“ inú.
19.2 Škálovanie premenných
Škálovanie môže byť veľmi efektívne z hľadiska numerickej stability (najmä pri veľmi rozdielnych rádoch premenných), znižuje však interpretovateľnosť modelu. Ide o úpravu premenných podľa nasledovného vzorca:
\[ x^{scale} = \frac{x-M}{\sqrt{D}} \]
kde \(M\) je stredná hodnota (priemer) a \(D\) je rozptyl premennej.
V našej databáze sú všetky premenné v počte obyvateľov, teda už sú v porovnateľných jednotkách.
udaje$Pop11_c <- scale(udaje$Pop11, center = TRUE, scale = TRUE)
udaje$Pop10_c <- scale(udaje$Pop10, center = TRUE, scale = TRUE)
udaje$Pop09_c <- scale(udaje$Pop09, center = TRUE, scale = TRUE)
model_centered <- lm(Pop12 ~ Pop11_c + Pop10_c + Pop09_c,
data = udaje)
summary(model_centered)##
## Call:
## lm(formula = Pop12 ~ Pop11_c + Pop10_c + Pop09_c, data = udaje)
##
## Residuals:
## Min 1Q Median 3Q Max
## -44.442 -2.185 0.749 3.659 50.908
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4930.9000 0.6262 7874.461 < 2e-16 ***
## Pop11_c 14803.9118 734.1607 20.164 < 2e-16 ***
## Pop10_c 1122.4034 1442.3446 0.778 0.437
## Pop09_c -4064.1462 973.5030 -4.175 4.48e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.856 on 196 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 1.19e+08 on 3 and 196 DF, p-value: < 2.2e-16vif(model_centered)## Pop11_c Pop10_c Pop09_c
## 1367712 5278985 2404844Condition Number je:
X <- model.matrix(model_centered)[, -1]
XtX <- t(X) %*% X
eig <- eigen(XtX)
condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number## [1] 4909.792Z výsledkov vidíme, že škálovanie môže zlepšiť numerické vlastnosti matice (Condition Number), ale multikolinearitu medzi mesačnými hodnotami ako takú neodstraňuje – dáta stále obsahujú takmer tú istú informáciu.
19.3 Iná úprava premennej, ktorá zachová interpretovateľnosť
Ak sa chceme vyhnúť strate interpretovateľnosti, môžeme niekedy zmeniť jednotky premennej, ktorá je v úplne iných rádoch než ostatné. V pôvodnom príklade to bola premenná GDP, ktorá bola v dolároch. V našej databáze sú všetky premenné v počte obyvateľov, takže tento problém nemáme. Pre ilustráciu si však ukážeme prevod jednej premennej na „tisíce obyvateľov“:
udaje$Pop11k <- udaje$Pop11 / 1000
head(udaje[, c("Obec", "Pop12", "Pop11", "Pop11k")])Potom lineárny model dosiahne výsledky:
model_Pop11k <- lm(Pop12 ~ Pop11k + Pop10 + Pop09,
data = udaje)
summary(model_Pop11k)##
## Call:
## lm(formula = Pop12 ~ Pop11k + Pop10 + Pop09, data = udaje)
##
## Residuals:
## Min 1Q Median 3Q Max
## -44.442 -2.185 0.749 3.659 50.908
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.70441 0.68551 -2.486 0.0137 *
## Pop11k 1248.62877 61.92243 20.164 < 2e-16 ***
## Pop10 0.09468 0.12167 0.778 0.4374
## Pop09 -0.34283 0.08212 -4.175 4.48e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.856 on 196 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 1.19e+08 on 3 and 196 DF, p-value: < 2.2e-16vif(model_Pop11k)## Pop11k Pop10 Pop09
## 1367712 5278985 2404844V tomto prípade sa koeficient pri Pop11k zmení len v jednotkách (na tisíce obyvateľov), ale samotná multikolinearita nezmizne, pretože Pop11, Pop10 a Pop09 sú stále veľmi podobné.
Poznámka: Dummy premenné (0/1) neškálujeme.
### očistenie databázy od „pracovných“ stĺpcov
library(dplyr)
udaje <- udaje %>%
dplyr::select(Okres, Obec, Pop12, Pop11, Pop10, Pop09)19.4 PCA (voliteľné)
X_pca <- scale(udaje[, c("Pop11", "Pop10", "Pop09")],
center = TRUE, scale = TRUE)
pca_res <- prcomp(X_pca)
summary(pca_res)## Importance of components:
## PC1 PC2 PC3
## Standard deviation 1.732 0.000992 0.0003528
## Proportion of Variance 1.000 0.000000 0.0000000
## Cumulative Proportion 1.000 1.000000 1.0000000udaje$PC1 <- pca_res$x[, 1]
model_pca <- lm(Pop12 ~ PC1, data = udaje)
summary(model_pca)##
## Call:
## lm(formula = Pop12 ~ PC1, data = udaje)
##
## Residuals:
## Min 1Q Median 3Q Max
## -80.322 -3.752 -0.201 3.493 112.824
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4930.9000 1.1468 4300 <2e-16 ***
## PC1 -6848.6258 0.6638 -10318 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.22 on 198 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 1.065e+08 on 1 and 198 DF, p-value: < 2.2e-16# VIF sa nedá počítať, lebo model má len jeden vysvetľujúci regresor
# vif(model_pca)19.5 Hrebeňová regresia (Ridge Regression – voliteľné)
Hrebeňová regresia je modifikácia regresie, ktorá zavádza perturbácie do matice \(\mathbf X^T \mathbf X\) tak, aby znížila dôsledky multikolinearity. Treba ale upozorniť, že odhadované regresné koeficienty sú skreslené. Perturbácia vyzerá nasledovne:
\[ (\mathbf X^T \mathbf X + \lambda \mathbf I) \]
V tomto kroku si vypíšeme výsledky odhadov regresných koeficientov s meniacimi sa parametrami \(\lambda\), aby sme získali určitú predstavu o číselnom ráde prvkov:
library(MASS)
X <- as.matrix(udaje[, c("Pop11", "Pop10", "Pop09")])
y <- udaje$Pop12
ridge_fit <- lm.ridge(y ~ X, lambda = seq(0, 10, 1))
ridge_fit## XPop11 XPop10 XPop09
## 0 -1.704409 1.2486288 0.09468307 -0.3428288
## 1 8.529458 0.3331203 0.33297224 0.3328418
## 2 16.706357 0.3324811 0.33243251 0.3323613
## 3 24.856002 0.3319011 0.33188566 0.3318343
## 4 32.978628 0.3313373 0.33133840 0.3312969
## 5 41.074389 0.3307811 0.33079207 0.3307565
## 6 49.143424 0.3302296 0.33024709 0.3302155
## 7 57.185869 0.3296815 0.32970366 0.3296750
## 8 65.201856 0.3291362 0.32916186 0.3291353
## 9 73.191517 0.3285934 0.32862173 0.3285968
## 10 81.154981 0.3280529 0.32808330 0.3280597Nastavovanie rôznych hodnôt \(\lambda\) je predmetom hlbšej analýzy, tento prehľad je len jej časťou.
20 10. Zhrnutie
- Multikolinearita nezavádza bias, ale zvyšuje štandardné odchýlky
odhadovaných regresných koeficientov a robí ich nestabilnými.
- Testy (resp. diagnostické ukazovatele) ako VIF a Condition Number
umožňujú diagnostiku multikolinearity.
- Riešenia zahŕňajú: vynechanie premennej, centrovanie/škálovanie, zmenu jednotiek, prípadne použitie PCA alebo hrebeňovej regresie. PCA a hrebeňová regresia vyžadujú hlbšie vedomosti, preto ich tu uvádzame len formálne a nevyžadujeme ich detailný výpočet.