Lineárna regresia na energetických údajoch (OWID)
V tomto cvičení nadväzujeme na databázu z Cvičenia 3 (energetický dataset Our World in Data (OWID) Energy Dataset) a vyskúšame si, ako odhadnúť a interpretovať jednoduchú a viacnásobnú lineárnu regresiu.
Sústredíme sa na nasledujúce kroky: - Načítanie a kontrola pripravených údajov. - Grafické znázornenie závislosti medzi vybranými premennými. - Odhad jednoduchej regresnej priamky a jej interpretácia. - Rozšírenie na viacnásobnú lineárnu regresiu. - Diagnostika modelu a práca s rezíduami.
library(broom)
library(knitr)
library(kableExtra)
library(ggplot2)
library(dplyr)
Attaching package: 'dplyr'The following object is masked from 'package:kableExtra':
group_rowsThe following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, unionlibrary(dplyr)
cesta <- "owid_energy.csv"
if (!file.exists(cesta)) cesta <- "owid-energy.csv"
owid_data <- read.csv(cesta, sep = ",", dec = ".",
na.strings = c("", "NA"),
stringsAsFactors = FALSE)
vybrane_krajiny <- c("Slovakia", "Czechia", "Austria", "Poland", "Germany")
udaje_vyber <- owid_data %>%
filter(country %in% vybrane_krajiny, year >= 2000) %>%
select(country, year, population, gdp, electricity_generation_twh, renewables_electricity_twh)
library(mice)
Attaching package: 'mice'The following object is masked from 'package:stats':
filterThe following objects are masked from 'package:base':
cbind, rbindcifrene_data <- udaje_vyber %>% select(population, gdp, electricity_generation_twh, renewables_electricity_twh)
imp <- mice(cifrene_data, m = 1, maxit = 5, seed = 123, printFlag = FALSE)
cifrene_imputovane <- complete(imp, 1)
# 4. Finálny objekt pre analýzu a regresie
udaje_analysa <- bind_cols(
country = udaje_vyber$country,
year = udaje_vyber$year,
cifrene_imputovane
)head(udaje_analysa) |>
kable(format = "html", digits = 2, row.names = FALSE,
caption = "Prvých šesť riadkov analytického datasetu") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| country | year | population | gdp | electricity_generation_twh | renewables_electricity_twh |
|---|---|---|---|---|---|
| Austria | 2000 | 8012929 | 280410163017 | 59.64 | 43.44 |
| Austria | 2001 | 8042233 | 285335695766 | 60.98 | 42.24 |
| Austria | 2002 | 8081640 | 291231427735 | 60.63 | 41.90 |
| Austria | 2003 | 8122604 | 294489280951 | 57.97 | 35.24 |
| Austria | 2004 | 8173496 | 303692479841 | 61.85 | 39.62 |
| Austria | 2005 | 8229204 | 311441138652 | 64.37 | 40.88 |
Predpokladáme, že s rastúcim HDP krajiny rastieme aj v celkovej výrobe elektrickej energie. Pozrime sa na túto závislosť graficky.
graf_regresia <- ggplot(udaje_analysa, aes(x = gdp, y = electricity_generation_twh)) +
geom_point(aes(color = country), size = 2, alpha = 0.75) +
labs(x = "HDP (štátny produkt)",
y = "Celková výroba elektriny (TWh)",
color = "Krajina") +
theme_minimal(base_size = 13)
graf_regresia
Tento graf znázorňuje vzťah a porovnanie celkovej výroby elektriny (v TWh) pre vybrané stredoeurópske krajiny (Rakúsko, Česko, Nemecko, Poľsko, Slovensko) v čase:
Nemecko vykazuje suverénne najvyššiu produkciu elektriny, ktorá sa pohybuje v rozmedzí od 500 do vyše 600 TWh.
Poľsko sa drží v strednej časti s produkciou okolo 150 až 180 TWh.
Česko a Rakúsko produkujú približne 60 až 90 TWh.
Slovensko má najnižšiu celkovú produkciu v rámci porovnatej vzorky, pohybujúcu sa tesne pod hranicou 40 TWh.
Najprv použijeme jednoduchú regresiu na modelovanie závislosti celkovej výroby elektriny od HDP:\[Electricity_i = \beta_0 + \beta_1 GDP_i + u_i\]
# HDP si prepočítame na miliardy
udaje_regresia <- udaje_vyber %>%
mutate(gdp_miliardy = gdp / 1e9) # ak by to bolo príliš malé/veľké, vieme prispôsobiť
# 2. Jednoduchá lineárna regresia
model_jednoduchy <- lm(electricity_generation_twh ~ gdp_miliardy, data = udaje_regresia)
# Zobrazenie koeficientov v tabuľke
tidy(model_jednoduchy) |>
kable(format = "html", digits = 4, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta", "Št. odchýlka", "Štatistika", "P-hodnota"),
caption = "Jednoduchá regresia: výroba elektriny a HDP (v miliardách)") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta | Št. odchýlka | Štatistika | P-hodnota |
|---|---|---|---|---|
| (Intercept) | 19.7444 | 4.2548 | 4.6405 | 0 |
| gdp_miliardy | 0.1719 | 0.0027 | 63.6883 | 0 |
Odhadnutý model ukazuje, že konštanta dosahuje hodnotu 19,74 TWh, čo predstavuje základnú úroveň výroby elektriny. Regresný koeficient pri premennej gdp_miliardy je 0,1719, čo znamená, že s nárastom HDP krajiny o 1 miliardu sa predikovaná výroba elektrickej energie zvyšuje v priemere o 0,1719 TWh. Extrémne nízka p-hodnota (blížiaca sa nule) potvrdzuje, že tento vzťah je štatisticky vysoko významný.
graf_regresia +
geom_smooth(method = "lm", color = "firebrick", se = FALSE, linewidth = 1)`geom_smooth()` using formula = 'y ~ x'
Výrobu elektriny ovplyvňuje nielen výkon ekonomiky (HDP), ale aj veľkosť populácie. Rozšírim model o premennú population:\[Electricity_i = \beta_0 + \beta_1 GDP_i + \beta_2 Population_i + u_i\]
model_viacnasobny <- lm(electricity_generation_twh ~ gdp_miliardy + population, data = udaje_regresia)
# Zobrazenie viacnásobného modelu
tidy(model_viacnasobny) |>
kable(format = "html", digits = 6, row.names = FALSE,
col.names = c("Člen modelu", "Odhad koeficienta", "Št. odchýlka", "Štatistika", "P-hodnota"),
caption = "Viacnásobná regresia: HDP a populácia") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Člen modelu | Odhad koeficienta | Št. odchýlka | Štatistika | P-hodnota |
|---|---|---|---|---|
| (Intercept) | 4.679783 | 4.352216 | 1.075264 | 0.284567 |
| gdp_miliardy | 0.120170 | 0.008426 | 14.261149 | 0.000000 |
| population | 0.000002 | 0.000000 | 6.390426 | 0.000000 |
Po rozšírení modelu o premennú populácie vidíme, že parciálny koeficient pri gdp_miliardy dosahuje hodnotu 0,1202, čo znamená, že pri nezmenenej veľkosti populácie zvýšenie HDP o 1 miliardu zvyšuje produkciu elektriny o 0,1202 TWh. Koeficient pri premennej population je kladný a štatisticky významný (p-hodnota < 0,001), čo potvrdzuje, že krajiny s väčším počtom obyvateľov vyprodukujú v priemere viac elektrickej energie.
suhrn <- glance(model_viacnasobny)
suhrn[c("nobs", "r.squared", "sigma")] |>
kable(format = "html", digits = c(0, 4, 2), row.names = FALSE,
col.names = c("Počet pozorovaní", "R²", "Odhad SD chyby (TWh)"),
caption = "Základné charakteristiky viacnásobného modelu") |>
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Počet pozorovaní | R² | Odhad SD chyby (TWh) |
|---|---|---|
| 115 | 0.9801 | 30.66 |
Pre každé pozorovanie si porovnáme reálnu hodnotu výroby s modelovou predikciou a vypočítame rezíduum (\(e_i = y_i - \hat{y}_i\)).
# 1. Vyberieme len tie riadky, ktoré nemajú NA v použitých premenných
udaje_ciste <- udaje_regresia %>%
filter(!is.na(electricity_generation_twh) & !is.na(gdp_miliardy) & !is.na(population))
# 2. Vytvorektor porovnania s rovnakým počtom riadkov (115)
porovnanie <- data.frame(
Krajina = udaje_ciste$country,
Rok = udaje_ciste$year,
Skutocnost = udaje_ciste$electricity_generation_twh,
Predikcia = fitted(model_viacnasobny),
Reziduum = residuals(model_viacnasobny)
)
head(porovnanie, 8) |>
kable(format = "html", digits = 2, row.names = FALSE,
col.names = c("Krajina", "Rok", "Skutočnosť (TWh)", "Predikcia (TWh)", "Rezíduum"),
caption = "Porovnanie skutočných a predikovaných hodnôt") |>
kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)| Krajina | Rok | Skutočnosť (TWh) | Predikcia (TWh) | Rezíduum |
|---|---|---|---|---|
| Austria | 2000 | 59.64 | 56.65 | 2.99 |
| Austria | 2001 | 60.98 | 57.31 | 3.67 |
| Austria | 2002 | 60.63 | 58.11 | 2.52 |
| Austria | 2003 | 57.97 | 58.59 | -0.62 |
| Austria | 2004 | 61.85 | 59.82 | 2.03 |
| Austria | 2005 | 64.37 | 60.88 | 3.49 |
| Austria | 2006 | 62.21 | 62.39 | -0.18 |
| Austria | 2007 | 62.83 | 64.02 | -1.19 |
Tabuľka porovnáva reálnu výrobu elektriny (Skutočnosť) s hodnotami odhadnutými viacnásobným regresným modelom (Predikcia). Rezíduum predstavuje rozdiel medzi týmito hodnotami (\(e_i = y_i - \hat{y}_i\)). Napríklad v roku 2000 v Rakúsku bola skutočná výroba 59,64 TWh, zatiaľ čo model predikoval 56,65 TWh, čo znamená kladné rezíduum 2,99 TWh (model mierne podhodnotil skutočnosť). Nízke hodnoty rezíduí naznačujú, že viacnásobný model dokáže pomerne presne kopírovať reálny vývoj výroby elektriny.