Úvod

Cieľom analýzy je preskúmať, či rozšírenie elektromobility v krajinách V4 súvisí s kvalitou životného prostredia meranou indexom environmentálnej výkonnosti (EPI). Pracujem s prečisteným panelovým datasetom z mojej bakalárskej práce (krajiny V4 × roky).

Charakteristika dát

Zdroj a pokrytie: Vlastný dataset zostavený pre BP (krajiny V4 × roky).
Hlavné premenné: EPI (zložený index kvality životného prostredia), BEV (počet batériových EV), BEV_PHEV (BEV + plug‑in hybridy), HDP na obyvateľa.
Stav dát: Dáta boli prečistené pred týmto cvičením; v tomto dokumente už neriešime imputačné či čistiace kroky.

Premenné

Závislá premenná: EPI – zložený index environmentálnej výkonnosti krajín.

Kľúčové vysvetľujúce premenné: BEV – počet batériových elektromobilov; BEV_PHEV – súčet BEV a plug-in hybridov (PHEV).

Kontrolné premenné: HDP na obyvateľa (ekonomická úroveň); podľa dostupnosti aj ďalšie sprievodné ukazovatele (CO₂, investičné/rozvojové indikátory).

Výskumná otázka a hypotézy

Súvisí vyšší počet (P)HEV s vyšším EPI, po zohľadnení ekonomickej úrovne krajín?

H0: Počet BEV/BEV_PHEV nemá významný vplyv na EPI.
H1: Vyšší BEV/BEV_PHEV je spojený s vyšším EPI.

Metodika v skratke

Používam lineárne modely OLS s log‑transformáciou počtových premenných (log1p), aby sa stabilizovala variancia a tlmili sa extrémy. Aby som sa vyhla kolinearite, odhadujem model s BEV a model s BEV_PHEV oddelene; HDP figuruje ako kontrola úrovne rozvoja. Diagnostiku robím cez štandardné grafy (Q–Q, Scale–Location, Leverage) a formálne testy (Breusch–Pagan). Inferenciu reportujem s robustnými SE (HC1), prípadne klastrovanými podľa krajiny pri panelových špecifikáciách.

Import údajov z mojej bakalárskej práce

Načítam csv súbor s dátami o elektromobilite vo V4 (počet BEV, BEV_PHEV, EPI, HDP, emisie CO₂, investičné ukazovatele…). Dáta sú už prečistené.

Interaktívna tabuľka s dátami

Deskriptívna štatistika

library(dplyr)
library(tidyr)

# 1) vyber číselné premenné

num_df <- udaje %>% dplyr::select(where(is.numeric))

# 2) deskriptíva

desc <- num_df %>%
summarise(across(
everything(),
list(
n      = ~sum(!is.na(.)),
miss   = ~sum(is.na(.)),
mean   = ~mean(., na.rm = TRUE),
sd     = ~sd(., na.rm = TRUE),
q25    = ~quantile(., 0.25, na.rm = TRUE),
median = ~median(., na.rm = TRUE),
q75    = ~quantile(., 0.75, na.rm = TRUE),
min    = ~min(., na.rm = TRUE),
max    = ~max(., na.rm = TRUE)
),
.names = "{.col}__{.fn}"
)) %>%
pivot_longer(everything(),
names_to = c("Premenna", "Statistika"),
names_sep = "__") %>%
pivot_wider(names_from = "Statistika", values_from = value) %>%
select(Premenna, n, miss, mean, sd, q25, median, q75, min, max)

# 3) prezentácia cez gt (ak je dostupné), inak jednoduchá kable

if (requireNamespace("gt", quietly = TRUE)) {
library(gt)
library(scales)

num_cols <- c("mean","sd","q25","median","q75","min","max")

gt(desc) |>
fmt_number(columns = all_of(num_cols), decimals = 2) |>
data_color(
columns = all_of(num_cols),
colors = col_numeric(
palette = c("#4575b4", "#f7f7f7", "#d73027"),  # modrá–biela–červená
domain = range(as.matrix(desc[, num_cols]), na.rm = TRUE)
)
) |>
tab_header(
title = "Deskriptívna štatistika (číselné premenné)"
)
} else {

# fallback bez farieb
  fig.width=6
  fig.height=4
  fig.align='center'
  dpi=96
plot(hc)

library(kableExtra)
desc %>%
mutate(across(where(is.numeric), ~round(.x, 2))) %>%
kbl(caption = "Deskriptívna štatistika (číselné premenné)") %>%
kable_classic(full_width = TRUE)
}

Premenna	n	miss	mean	sd	q25	median	q75	min	max
Deskriptívna štatistika (číselné premenné)
Rok	50	0	2,018.50	2.90	2,016.00	2,018.50	2,021.00	2,014.00	2,023.00
BEV	50	0	15,927.08	30,185.85	911.00	4,155.43	14,079.75	113.00	163,245.52
EPI	50	0	67.88	9.18	63.50	68.30	73.28	46.00	85.42
HDP	50	0	19,471.60	6,655.54	14,522.50	17,325.00	21,787.50	11,360.00	33,300.00
BEV_PHEV	50	0	27,711.88	54,787.41	1,575.50	7,198.50	28,483.25	201.00	290,224.40
co2	50	0	734,496,946.31	1,319,764,399.20	54,592,713.42	117,718,215.94	129,411,626.75	28,339,113.57	3,565,184,295.47
res_develop	50	0	1.52	0.48	1.07	1.45	1.93	0.79	2.31
r_sources	50	0	628.47	529.85	317.81	437.32	524.84	129.45	2,010.71
unvest_man	50	0	0.15	0.20	0.05	0.08	0.11	0.03	1.00
invest_sources	50	0	0.15	0.07	0.09	0.14	0.19	0.04	0.35
invest_transport	50	0	0.27	0.13	0.17	0.26	0.32	0.09	0.72

Prieskumné grafy - BOXPLOT

Prečo: pred modelovaním overujem, či mierky a rozdelenia dávajú zmysel a či nie sú viditeľné extrémy. Počtové premenné (BEV, BEV_PHEV) vizualizujem v log‑mierke (log1p).

par(mfrow = c(2, 2))
boxplot(udaje$EPI,                main = "EPI (bez log)",    col = "darkgreen",  horizontal = TRUE)
boxplot(log1p(udaje$BEV),         main = "log1p(BEV)",       col = "lightblue", horizontal = TRUE)
boxplot(log1p(udaje$BEV_PHEV),    main = "log1p(BEV_PHEV)",  col = "darkgreen",  horizontal = TRUE)
boxplot(log1p(udaje$HDP),         main = "log1p(HDP)",       col = "lightblue", horizontal = TRUE)

par(mfrow = c(1, 1))

Boxplot zobrazuje rozdelenie premennej – hrubá čiara v strede je medián, „krabička“ (box) je IQR (od 1. po 3. kvartil), „fúzy“ ukazujú typický rozsah a bodky mimo sú odľahlé hodnoty. V našich grafoch je EPI sústredené okolo ~70 s jedným nižším outlierom; po log-transformácii sú log1p(BEV) a log1p(BEV_PHEV) pekne stabilné a symetrické (bez dlhých chvostov) a log1p(HDP) má veľmi úzky rozsah – bude slúžiť najmä ako kontrola úrovne rozvoja.

Lineárna regresia (dve alternatívy)

Motivácia výberu špecifikácie: EPI je zložený index a elektromobilita môže byť korelovaná s celkovou vyspelosťou krajiny. Preto zahrniem HDP ako kontrolu. Keďže BEV a BEV_PHEV spolu silno súvisia, odhadujem ich oddelené modely, aby som predišla kolinearite. Koeficienty interpretujem ako semi‑elasticity: zmena v log1p(x) približne zodpovedá percentuálnej zmene počtu vozidiel.

Metóda A (Model A: EPI ~ log1p(BEV) + log1p(HDP)) Overuje, či samotné batériové elektromobily (BEV) súvisia s EPI po zohľadnení ekonomickej úrovne (HDP). Teda: keď dve krajiny majú rovnaký HDP, má vyšší počet BEV (v % zmenách – vďaka log1p) spojený vyšší/nižší EPI?

Metóda B (Model B: EPI ~ log1p(BEV_PHEV) + log1p(HDP)) Overuje vzťah medzi EPI a celkovým rozšírením (P)HEV (BEV + plug-in hybridy), opäť pri rovnakej úrovni HDP. Tento model hovorí, čo spraví spoločný“ (P)HEV indikátor s EPI.

Aby som sa vyhla kolinearite, odhadujem dve špecifikácie zvlášť:

# Model A: EPI ~ log1p(BEV) + log1p(HDP)
m_bev  <- lm(EPI ~ log1p(BEV)      + log1p(HDP), data = udaje)

# Model B: EPI ~ log1p(BEV_PHEV) + log1p(HDP)
m_phev <- lm(EPI ~ log1p(BEV_PHEV) + log1p(HDP), data = udaje)

summary(m_bev)

## 
## Call:
## lm(formula = EPI ~ log1p(BEV) + log1p(HDP), data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.1351  -4.4467  -0.8368   4.9316  15.2734 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -16.7362    36.6646  -0.456  0.65016    
## log1p(BEV)   -3.6278     0.6707  -5.409 2.09e-06 ***
## log1p(HDP)   11.6591     4.0185   2.901  0.00564 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.358 on 47 degrees of freedom
## Multiple R-squared:  0.3839, Adjusted R-squared:  0.3577 
## F-statistic: 14.64 on 2 and 47 DF,  p-value: 1.14e-05

summary(m_phev)

## 
## Call:
## lm(formula = EPI ~ log1p(BEV_PHEV) + log1p(HDP), data = udaje)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.005  -4.389  -1.012   4.244  14.395 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -24.7863    35.9496  -0.689  0.49392    
## log1p(BEV_PHEV)  -3.8147     0.6564  -5.811 5.19e-07 ***
## log1p(HDP)       12.8282     3.9667   3.234  0.00224 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.149 on 47 degrees of freedom
## Multiple R-squared:  0.4184, Adjusted R-squared:  0.3936 
## F-statistic:  16.9 on 2 and 47 DF,  p-value: 2.947e-06

# Robustné SE (HC1)
coeftest(m_bev,  vcov = vcovHC(m_bev,  type = "HC1"))

## 
## t test of coefficients:
## 
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -16.73620   37.34531 -0.4481 0.656106    
## log1p(BEV)   -3.62776    0.76151 -4.7639 1.86e-05 ***
## log1p(HDP)   11.65906    4.11228  2.8352 0.006733 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(m_phev, vcov = vcovHC(m_phev, type = "HC1"))

## 
## t test of coefficients:
## 
##                  Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)     -24.78627   37.46093 -0.6617  0.511422    
## log1p(BEV_PHEV)  -3.81473    0.77914 -4.8961 1.194e-05 ***
## log1p(HDP)       12.82820    4.15931  3.0842  0.003413 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Porovnanie kvality
aic_tbl <- AIC(m_bev, m_phev)
aic_tbl

# Zvolený model na diagnostiku
model <- m_phev

Z lineárnych regresií vychádza, že HDP na obyvateľa má na EPI pozitívny a štatisticky významný vplyv, kým počet (P)HEV (a podobne aj samotných BEV) je záporný a významný po zohľadnení HDP; model s BEV_PHEV sedí mierne lepšie (nižšie AIC) a závery sa nemenia ani pri robustných smerodajných chybách. Modely vysvetľujú zhruba 38–42 % variability EPI.

# ak m_lin ešte neexistuje, použijem lineárny model s (P)HEV
if (!exists("m_lin")) m_lin <- m_phev
# 1) Skontrolujeme a prípadne opravíme názov stĺpca EPI
if (!"EPI" %in% names(udaje)) {
  col_epi <- grep("^EPI\\b", names(udaje), value = TRUE)[1]
  if (!is.na(col_epi)) names(udaje)[names(udaje) == col_epi] <- "EPI"
}

# 2) Vyberieme model, z ktorého chceme predikovať
mod <- m_lin   # <- prípadne zmeň na: m_phev, m_slope, m_quad, ...

# 3) Predikčná krivka: EPI ~ BEV_PHEV pri mediáne HDP
pdat <- data.frame(
  BEV_PHEV = seq(min(udaje$BEV_PHEV, na.rm = TRUE),
                 max(udaje$BEV_PHEV, na.rm = TRUE),
                 length.out = 200),
  HDP = median(udaje$HDP, na.rm = TRUE)
)

pr <- predict(mod, newdata = pdat, se.fit = TRUE)
pdat$yhat <- pr$fit
pdat$lo   <- pr$fit - 1.96 * pr$se.fit
pdat$hi   <- pr$fit + 1.96 * pr$se.fit

library(ggplot2)
ggplot(pdat, aes(x = BEV_PHEV, y = yhat)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.2) +
  geom_line() +
  labs(x = "BEV_PHEV", y = "Predikované EPI",
       title = "EPI vs. (P)HEV (pri mediáne HDP)") +
  theme_minimal()

Graf ukazuje predikovaný EPI v závislosti od počtu (P)HEV pri fixovanom (mediánovom) HDP z nášho semi-log modelu. Krivka má zreteľne klesajúci a splošťujúci sa tvar: pri nízkych hodnotách (P)HEV je pokles EPI strmší, no s rastúcim (P)HEV sa účinok zmierňuje (diminishing returns), čo je typické pre špecifikáciu s log⁡(1+(P)HEV. Šedé pásmo vyznačuje 95 % interval predikcie a smerom k vyšším hodnotám sa rozširuje, teda neistota rastie (v tých oblastiach máme menej pozorovaní alebo väčší rozptyl). Celkovo graf vizuálne potvrdzuje negatívnu asociáciu medzi rozšírením (P)HEV a EPI v našich dátach po zohľadnení HDP; ide o asociáciu, nie dôkaz kauzality.

Diagnostika (krajšie grafy)

Residuals vs Fitted

– ukazuje, či model nechýba tvar (nelinearita) a či sú chyby rovnomerne okolo nuly. Hľadáme „náhodný mrak“ bez vzoru.

Normal Q–Q

– porovnáva rozdelenie rezíduí s ideálnou normálnou krivkou; body pri priamke = približná normalita, odchýlky v chvostoch signalizujú extrémy.

Scale–Location (√|štandardizované rezíduá| vs. fitted)

– test rovnakého rozptylu chýb (homoskedasticita). Rovná LOESS krivka ≈ konštantný rozptyl; stúpajúci/lievik = heteroskedasticita.

Residuals vs Leverage

– identifikuje vplyvné pozorovania: kombinácia veľkého „leverage“ (hat hodnoty) a veľkých rezíduí. Pomáha rozhodnúť, či niekoľko bodov neťahá regresnú priamku (podľa Cookovej vzdialenosti).

clr <- "#0E6B4D"

# diagnostický dataframe 
diag_df <- augment(model) %>% mutate(id = dplyr::row_number())

cook_cut <- 4 / nrow(diag_df)
lab_df   <- dplyr::filter(diag_df, .cooksd > cook_cut)

p1 <- ggplot(diag_df, aes(.fitted, .resid)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey60") +
  geom_point(alpha = .75, color = clr) +
  geom_smooth(method = "loess", se = FALSE, color = clr) +
  labs(title = "Residuals vs Fitted", x = "Fitted values", y = "Residuals") +
  theme_minimal(base_size = 12) + theme(plot.title = element_text(face = "bold"))

p2 <- ggplot(diag_df, aes(sample = .std.resid)) +
  stat_qq(color = clr, alpha = .9) +
  stat_qq_line(color = "grey45") +
  labs(title = "Normal Q–Q", x = "Theoretical Quantiles", y = "Standardized residuals") +
  theme_minimal(base_size = 12) + theme(plot.title = element_text(face = "bold"))

p3 <- ggplot(diag_df, aes(.fitted, sqrt(abs(.std.resid)))) +
  geom_point(alpha = .75, color = clr) +
  geom_smooth(method = "loess", se = FALSE, color = clr) +
  labs(title = "Scale–Location", x = "Fitted values",
       y = expression(sqrt("|Standardized residuals|"))) +
  theme_minimal(base_size = 12) + theme(plot.title = element_text(face = "bold"))

p4 <- ggplot(diag_df, aes(.hat, .std.resid)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey60") +
  geom_point(alpha = .75, color = clr) +
  geom_smooth(method = "loess", se = FALSE, color = clr) +
  geom_text(data = lab_df, aes(label = id), vjust = -0.5, size = 3) +
  labs(title = "Residuals vs Leverage", x = "Hat values", y = "Standardized residuals") +
  theme_minimal(base_size = 12) + theme(plot.title = element_text(face = "bold"))

(p1 + p2) / (p3 + p4)

Diagnostika – čítanie grafov:

Residuals vs Fitted: jemná nelinearita a mierne kolísanie rozptylu.
Q–Q: drobné odchýlky v chvostoch → približná normalita.
Scale–Location: slabá heteroskedasticita → použijem robustné SE.
Residuals vs Leverage: zopár vplyvných bodov (označené ID), ale bez extrémov.

Keď zohľadníme HDP, viac (P)HEV sa v našich dátach spája s nižším EPI. Kontrolné grafy neukázali vážne problémy a drobné odchýlky sme pokryli robustnými chybami, takže výsledok berieme ako spoľahlivý.

Vplyvné pozorovania

influencePlot(model, main = "Influence plot")

outlierTest(model)  # Bonferroni p-value

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##     rstudent unadjusted p-value Bonferroni p
## 20 -2.300723           0.025991           NA

Krátke čítanie grafu a tabuľky:

Influence plot (bubble chart)

-X-os (Hat-values) = „páka“ (leverage): ak je bod viac vpravo, má neštandardnú kombináciu vysvetľujúcich premenných a vie viac strhnúť priamku.

-Y-os (Studentized residuals) = veľkosť chyby: nad +2 alebo pod −2 je podozrivé.

-Veľkosť/odtieň bubliny = Cookovo D: čím väčšia a tmavšia bublina, tým väčší vplyv jednotky na odhad celého modelu.

Graf vplyvu ukazuje niekoľko pozorovaní s vyššou „pákou“ a vplyvom (najmä body 31, 41 a 20), no zostávajú v rozumných medziach; Cookovo D nie je extrémne a hodnoty študentizovaných rezíduí sa pohybujú okolo hraníc ±2. Formálny outlierTest pritom nenašiel žiadny štatisticky významný odľahlý bod po Bonferroni korekcii (najväčší |rstudent| má bod 20 ≈ −2.30, neopravené p ~ 0.026, po korekcii nevýznamné). Inými slovami: v dátach sú jednotky, ktoré majú citeľnejší vplyv na odhad, ale nemáme dôkaz o skutočných outlieroch # Dodatočné porovnanie modelov

# Model A: EPI ~ log1p(BEV) + HDP (bez log na HDP – ako kontrola úrovne)
m_bev2  <- lm(EPI ~ log1p(BEV)      + HDP, data = udaje)
# Model B: EPI ~ log1p(BEV_PHEV) + HDP
m_phev2 <- lm(EPI ~ log1p(BEV_PHEV) + HDP, data = udaje)
summary(m_bev2); coeftest(m_bev2, vcov = vcovHC(m_bev2, type = "HC1"))

## 
## Call:
## lm(formula = EPI ~ log1p(BEV) + HDP, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.0909  -4.1711  -0.7445   5.2535  14.7114 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 86.8956784  4.7398188  18.333  < 2e-16 ***
## log1p(BEV)  -3.5939951  0.6665698  -5.392 2.21e-06 ***
## HDP          0.0005466  0.0001900   2.877  0.00603 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.367 on 47 degrees of freedom
## Multiple R-squared:  0.3823, Adjusted R-squared:  0.356 
## F-statistic: 14.54 on 2 and 47 DF,  p-value: 1.211e-05

## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept) 86.89567837  5.43343513 15.9928 < 2.2e-16 ***
## log1p(BEV)  -3.59399513  0.75667089 -4.7497  1.95e-05 ***
## HDP          0.00054656  0.00019420  2.8144  0.007115 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(m_phev2); coeftest(m_phev2, vcov = vcovHC(m_phev2, type = "HC1"))

## 
## Call:
## lm(formula = EPI ~ log1p(BEV_PHEV) + HDP, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.8667  -4.2839  -0.6805   4.3720  13.7936 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     89.2739287  4.7887883  18.642  < 2e-16 ***
## log1p(BEV_PHEV) -3.8005552  0.6530506  -5.820 5.04e-07 ***
## HDP              0.0006089  0.0001877   3.243  0.00218 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.145 on 47 degrees of freedom
## Multiple R-squared:  0.419,  Adjusted R-squared:  0.3942 
## F-statistic: 16.94 on 2 and 47 DF,  p-value: 2.877e-06

## 
## t test of coefficients:
## 
##                    Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)     89.27392868  5.80102043 15.3893 < 2.2e-16 ***
## log1p(BEV_PHEV) -3.80055516  0.78549751 -4.8384  1.45e-05 ***
## HDP              0.00060888  0.00019776  3.0789  0.003464 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC(m_bev2, m_phev2)

Po zohľadnení HDP vychádza negatívna asociácia EPI s (P)HEV aj s BEV; mierne lepšie sedí špecifikácia s BEV_PHEV.

Oneskorený vplyv (lag 1 rok) s fixnými efektmi

Prečo lag: zmeny v zložení vozového parku a budovaní infraštruktúry sa nemusia okamžite premietnuť do EPI. Jednoročný lag BEV_PHEV s fixnými efektmi krajín a rokov zachytáva nepozorovanú heterogenitu a spoločné šoky.

udaje_lag <- udaje %>%
  arrange(Krajina, Rok) %>%
  group_by(Krajina) %>%
  mutate(lag_BEV_PHEV = dplyr::lag(BEV_PHEV, 1)) %>%
  ungroup()

m_lag <- lm(EPI ~ log1p(lag_BEV_PHEV) + HDP + factor(Krajina) + factor(Rok), data = udaje_lag)
coeftest(m_lag, vcov = vcovCL(m_lag, cluster = ~ Krajina, type = "HC1"))

## 
## t test of coefficients:
## 
##                           Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)             50.8959137  65.0340670  0.7826 0.4399939    
## log1p(lag_BEV_PHEV)     -2.0425584   1.5285689 -1.3363 0.1915136    
## HDP                      0.0020646   0.0038335  0.5386 0.5941504    
## factor(Krajina)EU-avrg -20.6979705  37.5633150 -0.5510 0.5857037    
## factor(Krajina)MR        9.5261676  25.1708980  0.3785 0.7077543    
## factor(Krajina)PL        7.3166762  26.7654445  0.2734 0.7864475    
## factor(Krajina)SK        4.1908401  12.0927337  0.3466 0.7313421    
## factor(Rok)2016          9.5399727   2.4016859  3.9722 0.0004121 ***
## factor(Rok)2017         -1.4413973   3.0279243 -0.4760 0.6374981    
## factor(Rok)2018         -8.1105733   5.3325433 -1.5210 0.1387416    
## factor(Rok)2019        -10.1512429   6.6783228 -1.5200 0.1389742    
## factor(Rok)2020         -7.7457174   3.0950422 -2.5026 0.0180056 *  
## factor(Rok)2021        -15.1063011   6.3794152 -2.3680 0.0245284 *  
## factor(Rok)2022        -18.9950331   7.2001178 -2.6382 0.0130861 *  
## factor(Rok)2023         -6.9375471   8.2719932 -0.8387 0.4082826    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Aj keď sa pozrieme na oneskorený (1-ročný) efekt a zafixujeme rozdiely medzi krajinami a rokmi, nevidíme štatistický dôkaz, že viac (P)HEV zlepší EPI—koeficient je síce negatívny, ale nevýznamný. Pre našu prácu to znamená, že hlavný záver zostáva rovnaký: po zohľadnení HDP (a pri FE) sa v dostupných dátach V4 nepotvrdzuje krátkodobý pozitívny vplyv rozšírenia (P)HEV na EPI. Model teda skôr hovorí, že EPI formujú aj iné, širšie faktory než samotná elektromobilita

Štandardizované koeficienty

Zmysel: po štandardizácii viem porovnať veľkosť efektov naprieč rozdielnymi mierkami (SD jednotky).

m_std <- lm(scale(EPI) ~ scale(log1p(BEV_PHEV)) + scale(HDP), data = udaje)
summary(m_std)

## 
## Call:
## lm(formula = scale(EPI) ~ scale(log1p(BEV_PHEV)) + scale(HDP), 
##     data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.61936 -0.46663 -0.07412  0.47622  1.50247 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             4.914e-16  1.101e-01   0.000  1.00000    
## scale(log1p(BEV_PHEV)) -7.921e-01  1.361e-01  -5.820 5.04e-07 ***
## scale(HDP)              4.414e-01  1.361e-01   3.243  0.00218 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7783 on 47 degrees of freedom
## Multiple R-squared:  0.419,  Adjusted R-squared:  0.3942 
## F-statistic: 16.94 on 2 and 47 DF,  p-value: 2.877e-06

Po štandardizácii vychádza, že log1p(BEV_PHEV) má silnejší vplyv na EPI a je negatívny (~−0.79 SD za +1 SD v (P)HEV), kým HDP je pozitívne (~+0.44 SD za +1 SD v HDP); oba efekty sú štatisticky významné. Intercept je ~0 (lebo premenné sú centrované) a model vysvetľuje ~42 % variability EPI. Inými slovami: v našich dátach sa zmeny v (P)HEV spájajú s väčším (negatívnym) posunom EPI než rovnako veľké zmeny v HDP s pozitívnym smerom.

Heteroskedasticita – testy a grafy

Heteroskedasticita = rozptyl chýb nie je všade rovnaký. Ak by bola výrazná, bežné p-hodnoty môžu byť skreslené. Nižšie ju rýchlo overím (grafy + testy) a zohľadním cez robustné štandardné chyby (HC1).

Vizuálna diagnostika

library(broom); library(ggplot2); library(patchwork); clr <- "#0E6B4D"

h <- augment(model, data = udaje) |>
  dplyr::mutate(
    absres = sqrt(abs(.std.resid)),
    res2   = .resid^2,
    x_bev  = log1p(.data$BEV_PHEV),  # ak používame BEV: log1p(.data$BEV)
    x_hdp  = .data$HDP
  )

p_SL <- ggplot(h, aes(.fitted, absres)) +
  geom_point(alpha=.65) +
  geom_smooth(method="loess", se=FALSE, linewidth=1, color=clr) +
  labs(x="Predikované hodnoty", y=expression(sqrt("|rezíduá|")), title="Scale–Location") +
  theme_minimal()

p_bev <- ggplot(h, aes(x_bev, res2)) +
  geom_point(alpha=.65) +
  geom_smooth(method="loess", se=FALSE, linewidth=1, color=clr) +
  labs(x="log1p(BEV_PHEV)", y="rezíduá^2", title="Rezíduá^2 vs. log1p(BEV_PHEV)") +
  theme_minimal()

p_hdp <- ggplot(h, aes(x_hdp, res2)) +
  geom_point(alpha=.65) +
  geom_smooth(method="loess", se=FALSE, linewidth=1, color=clr) +
  labs(x="HDP", y="rezíduá^2", title="Rezíduá^2 vs. HDP") +
  theme_minimal()

p_SL / (p_bev | p_hdp)

LOESS krivka v Scale–Location má jemné „U“ – rozptyl rezíduí je najmenší pri stredných predikciách a väčší na okrajoch, takže ide o miernu heteroskedasticitu. Vo vzťahu rezíduá² vs. log1p(BEV_PHEV) je U-tvar len plytký, čiže s (P)HEV súvisí rozptyl skôr slabo. Naopak, graf rezíduá² vs. HDP ukazuje výraznejší U-tvar: variancia je najnižšia pri stredných hodnotách HDP a rastie pri nízkych aj vysokých, čo naznačuje, že prípadná heteroskedasticita je viazaná najmä na HDP. Celkovo nejde o silné porušenie, ale je rozumné používať robustné štandardné chyby (HC1); nižšie to ešte overíme formálnymi testami.

library(ggplot2)
library(patchwork)   

# model, ktorý chceme diagnostikovať

mod   <- m_phev     # EPI ~ log1p(BEV_PHEV) + log1p(HDP)
dta   <- udaje

# vypočítam rezíduá a pripravím log premennú pre BEV_PHEV (aby sedela s modelom)

dta$resid2   <- resid(mod)^2
dta$logBEVP  <- log1p(dta$BEV_PHEV)

# 1) Rezíduá^2 vs. log1p(BEV_PHEV)

p_bev <- ggplot(dta, aes(x = logBEVP, y = resid2)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "loess", se = FALSE, color = "firebrick") +
labs(x = "log1p(BEV_PHEV)", y = "Štvorce rezíduí",
title = "Rezíduá^2 vs. log1p(BEV_PHEV)") +
theme_minimal()

# 2) Rezíduá^2 vs. HDP  (pozor: v modeli je log(HDP) – tu úmyselne dávame aj „raw“ HDP,

# aby bolo jasne vidno, či problém súvisí s úrovňou bohatstva)

p_hdp <- ggplot(dta, aes(x = HDP, y = resid2)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "loess", se = FALSE, color = "firebrick") +
labs(x = "HDP na obyvateľa", y = "Štvorce rezíduí",
title = "Rezíduá^2 vs. HDP") +
theme_minimal()

p_bev | p_hdp

Z týchto dvoch grafov vidno, že veľkosť chýb (štvorce rezíduí) nie je všade rovnaká. Červená LOESS krivka má v oboch prípadoch tvar jemného „U“ – rozptyl rezíduí je najmenší v strede a väčší pri veľmi nízkych aj veľmi vysokých hodnotách. Efekt je zreteľnejší pri HDP, čo naznačuje, že nerovnaký rozptyl chýb súvisí skôr s úrovňou bohatstva krajín než s počtom (P)HEV. Nie je to extrémny „lievik“, ale ide o miernu heteroskedasticitu.

Heteroskedasticita - prehľad testov

library(lmtest); library(car); library(dplyr); library(knitr)

bp_fitted <- bptest(model)
bp_white  <- bptest(model, ~ fitted(model) + I(fitted(model)^2))
ncv       <- ncvTest(model)

tribble(
  ~Test,               ~Statistic,                   ~p_value,
  "Breusch–Pagan",     unname(bp_fitted$statistic),  bp_fitted$p.value,
  "BP (White štýl)",   unname(bp_white$statistic),   bp_white$p.value,
  "NCV (car)",         unname(ncv$ChiSquare),        ncv$p
) |>
  mutate(across(c(Statistic, p_value), ~round(.x, 4))) |>
  kable(caption = "Heteroskedasticita – prehľad testov")

Heteroskedasticita – prehľad testov
Test	Statistic	p_value
Breusch–Pagan	0.4876	0.7836
BP (White štýl)	5.6585	0.0591
NCV (car)	0.0032	0.9545

Formálne testy

Breusch–Pagan: p = 0,784 (nevýznamné).
White-štýl: p = 0,059 (hraničný náznak U-tvaru vo variancii).
NCV (car): p = 0,955 (nevýznamné).

Rozhodnutie: Testy nedávajú silný dôkaz heteroskedasticity; vzhľadom na jemný vzor v grafoch je primerané reportovať robustné štandardné chyby (HC1). Kvalitatívne závery modelu sa tým nemenia.

Nelineárne špecifikácie a testy funkčnej formy

Prečo toto robíme? Doteraz sme predpokladali, že vzťah medzi EPI a vysvetľujúcimi premennými (log1p(BEV_PHEV), log1p(HDP)) je lineárny. Nižšie skontrolujeme, či obyčajná „priamka“ naozaj stačí – a ak nie, pridáme jemné zakrivenie (kvadráty) alebo zlom v sklone (iný efekt pri vyšších hodnotách (P)HEV).

Zjednodušene: RESET test: rýchla kontrola, či priamka stačí. Malá p-hodnota ⇒ pridať nelineárne prvky. C+R grafy: ukážu, kde sa krivka ohýba – ktorú premennú transformovať. Kvadráty log-premenných: dovolia jemné zakrivenie. Zlom v sklone: dovolí iný vplyv (P)HEV pri nižších vs. vyšších hodnotách. Výber modelu: porovnáme AIC a (ak sú modely v tej istej mierke) aj ANOVA; koeficienty vždy s robustnými SE (HC1).

1) Ramsey RESET – kontrola špecifikácie

library(lmtest)

# Základný (lineárny) model – používame ten, s ktorým pracuješ v práci:

m_lin <- lm(EPI ~ log1p(BEV_PHEV) + log1p(HDP), data = udaje)

# Ramsey RESET test

resettest(m_lin)

## 
##  RESET test
## 
## data:  m_lin
## RESET = 4.9787, df1 = 2, df2 = 45, p-value = 0.01114

RESET vyšiel p = 0.011 → základná lineárna špecifikácia pravdepodobne chýba (nelinearita/nezahrnutý tvar). Preto ďalej cielene hľadáme kde by sa mohla krivka ohýbať a či sa oplatí pridať kvadráty alebo zlom.

2) Component + Residual (C+R) grafy – kde sa krivka ohýba

library(car)
car::crPlots(m_lin)   # pozrieme zakrivenie pri log1p(BEV_PHEV) a log1p(HDP)

C+R pre log1p(BEV_PHEV) ukazuje citeľné prehnutie okolo ~7–9 → kandidát na nelineárny tvar. C+R pre log1p(HDP) je prakticky priamka (len jemná krivka), takže výraznu nelinearitu pri HDP nečakáme. Z toho dôvodu testujeme hlavne úpravy pri (P)HEV.

3) Kvadratické členy – jemné zakrivenie

library(sandwich)
library(lmtest)

m_quad <- lm(EPI ~ log1p(BEV_PHEV) + I(log1p(BEV_PHEV)^2) +
log1p(HDP)      + I(log1p(HDP)^2),
data = udaje)

# Porovnanie kvality (nižší AIC je lepší)

AIC(m_lin, m_quad)

# ANOVA:

anova(m_lin, m_quad)

# Koeficienty s robustnými SE (HC1)

coeftest(m_lin,  vcov = vcovHC(m_lin,  type = "HC1"))

## 
## t test of coefficients:
## 
##                  Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)     -24.78627   37.46093 -0.6617  0.511422    
## log1p(BEV_PHEV)  -3.81473    0.77914 -4.8961 1.194e-05 ***
## log1p(HDP)       12.82820    4.15931  3.0842  0.003413 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(m_quad, vcov = vcovHC(m_quad, type = "HC1"))

## 
## t test of coefficients:
## 
##                        Estimate Std. Error t value Pr(>|t|)
## (Intercept)          1175.79233 1220.85177  0.9631   0.3406
## log1p(BEV_PHEV)         3.88313    6.57887  0.5902   0.5580
## I(log1p(BEV_PHEV)^2)   -0.44168    0.38178 -1.1569   0.2534
## log1p(HDP)           -236.25676  249.85488 -0.9456   0.3494
## I(log1p(HDP)^2)        12.57104   12.62555  0.9957   0.3247

AIC sa zhoršil (345.0 vs. 343.5).
ANOVA: p = 0.33 → kvadráty nezlepšujú model.
Robustné SE (HC1): kvadratické koeficienty nevýznamné.

Záver: kvadráty nepriniesli prínos – nelinearitu takto nepotvrdzujeme.

4) Zlom v sklone – iný efekt (P)HEV pri vyšších hodnotách

# Prah: medián log1p(BEV_PHEV)

cut_bev <- median(log1p(udaje$BEV_PHEV), na.rm = TRUE)
udaje$D_highBEV <- as.integer(log1p(udaje$BEV_PHEV) >= cut_bev)

# (a) len posun (intercept shift)

m_shift <- lm(EPI ~ D_highBEV + log1p(BEV_PHEV) + log1p(HDP), data = udaje)

# (b) zlom v sklone (interakcia premenná × dummy)

m_slope <- lm(EPI ~ log1p(BEV_PHEV)*D_highBEV + log1p(HDP), data = udaje)

# Porovnanie AIC

AIC(m_lin, m_shift, m_slope)

# ANOVA (rovnaká vysvetľovaná premenná)

anova(m_lin, m_shift)

anova(m_lin, m_slope)

# Robustné SE

coeftest(m_slope, vcov = vcovHC(m_slope, type = "HC1"))

## 
## t test of coefficients:
## 
##                           Estimate Std. Error t value Pr(>|t|)   
## (Intercept)               -29.6808    39.5065 -0.7513 0.456389   
## log1p(BEV_PHEV)            -2.5330     1.5078 -1.6798 0.099918 . 
## D_highBEV                   7.4189    20.1690  0.3678 0.714720   
## log1p(HDP)                 12.4336     4.3666  2.8474 0.006619 **
## log1p(BEV_PHEV):D_highBEV  -1.1895     2.1857 -0.5442 0.588969   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Prah = medián log1p(BEV_PHEV). Skúsili sme (a) posun priamky a (b) zmenu sklonu (interakcia s dummy):

AIC je vyšší než v lineárnom modeli (horšie).

ANOVA: posun p = 0.41, zmena sklonu p = 0.52.

Robustné SE: interakcia nevýznamná.

Záver: dôkaz pre odlišný efekt (P)HEV pri „vyšších“ úrovniach sme nenašli.

5) Výber pracovného modelu a krátke zhrnutie

Keďže (i) RESET síce naznačil nesprávnu špecifikáciu, ale (ii) C+R ukázal len jemné ohyby a (iii) ani kvadráty, ani zlom štatisticky nezlepšili model, ako pracovný model ponechávame základnú lineárnu špecifikáciu s robustnými chybami (HC1). V nej vychádza log1p(BEV_PHEV) negatívne a významne, kým log1p(HDP) pozitívne a významne. Prakticky: v tomto súbore dát V4 nevychádza, že by mierne zakrivenia alebo prahové efekty menili hlavné závery; ak je nelinearita prítomná, je skôr slabá a naše testy ju nevedia presvedčivo uchopiť (limit vzorky).

Porovnávacia tabuľka modelov

# istota, že máme základný lineárny model
if (!exists("m_lin")) {
  m_lin <- lm(EPI ~ log1p(BEV_PHEV) + log1p(HDP), data = udaje)
}

mods <- list("Lineárny" = m_lin)
if (exists("m_quad"))   mods[["Kvadratický"]]   <- m_quad
if (exists("m_slope"))  mods[["Zlom v sklone"]] <- m_slope
if (exists("m_shift"))  mods[["Posun úrovne"]]  <- m_shift

library(purrr); library(lmtest); library(dplyr); library(kableExtra)

reset_p <- function(m) tryCatch(resettest(m)$p.value, error = function(e) NA_real_)
anova_vs_base <- function(m, base) {
  if (identical(formula(base)[[2]], formula(m)[[2]])) as.numeric(anova(base, m)$`Pr(>F)`[2]) else NA_real_
}

base_model <- mods[["Lineárny"]]
cmp_tbl <- tibble(
  model = names(mods),
  AIC   = map_dbl(mods, AIC),
  adjR2 = map_dbl(mods, ~ summary(.x)$adj.r.squared),
  RESET_p = map_dbl(mods, reset_p),
  ANOVA_vs_lin_p = map_dbl(mods, ~ anova_vs_base(.x, base_model))
) %>% arrange(AIC)

kbl(cmp_tbl, digits = 3, caption = "Porovnanie špecifikácií (nižšie AIC je lepšie)") %>%
  kable_classic(full_width = FALSE)

Porovnanie špecifikácií (nižšie AIC je lepšie)
model	AIC	adjR2	RESET_p	ANOVA_vs_lin_p
Lineárny	343.498	0.394	0.011	NA
Posun úrovne	344.750	0.390	0.025	0.409
Kvadratický	345.035	0.397	0.007	0.330
Zlom v sklone	346.361	0.381	0.032	0.600

Graf predikčných kriviek (Lineárny vs. Kvadratický vs. Zlom)

# Predikcie pri mediáne HDP
x_grid <- tibble(
  BEV_PHEV = seq(min(udaje$BEV_PHEV, na.rm = TRUE),
                 max(udaje$BEV_PHEV, na.rm = TRUE),
                 length.out = 200),
  HDP = median(udaje$HDP, na.rm = TRUE)
)

pred_one <- function(m, name) {
  # (ak model potrebuje dummy D_highBEV, dopočítame ju podľa prahu, ktorý používame v práci)
  newdat <- x_grid
  if ("D_highBEV" %in% names(model.frame(m))) {
    cut_bev <- median(log1p(udaje$BEV_PHEV), na.rm = TRUE)
    newdat$D_highBEV <- as.integer(log1p(newdat$BEV_PHEV) >= cut_bev)
  }
  pr <- predict(m, newdata = newdat, se.fit = TRUE)
  tibble(model = name,
         BEV_PHEV = newdat$BEV_PHEV,
         fit = pr$fit,
         lo  = pr$fit - 1.96*pr$se.fit,
         hi  = pr$fit + 1.96*pr$se.fit)
}

pred_all <- bind_rows(
  pred_one(m_lin,  "Lineárny"),
  pred_one(m_quad, "Kvadratický"),
  pred_one(m_slope,"Zlom v sklone")
)

ggplot(pred_all, aes(BEV_PHEV, fit, color = model, fill = model)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.10, color = NA) +
  geom_line(size = 1) +
  labs(x = "(P)HEV (počet)", y = "Predikované EPI",
       title = "Predikčné krivky: EPI vs. (P)HEV pri mediáne HDP",
       subtitle = "Porovnanie špecifikácií (lineárny, kvadratický, zlom v sklone)") +
  theme_minimal() +
  theme(legend.position = "bottom")

Graf porovnáva predikčné krivky EPI vs. (P)HEV pri mediánovom HDP pre tri špecifikácie (lineárna v loge, kvadratická, a „zlom v sklone“). Všetky tri modely ukazujú rovnaký príbeh: s rastúcim (P)HEV EPI klesá a efekt sa postupne splošťuje (najmä pri nízkych hodnotách je pokles strmší). Rozdiely medzi krivkami sú malé – kvadratická verzia je o niečo prudšia v chvoste, zatiaľ čo „zlom v sklone“ sa prakticky prekrýva s lineárnym modelom. Intervaly neistoty (priesvitné pásma) sa rozširujú v extrémoch, kde je menej dát. Spolu s výsledkami AIC/ANOVA vyššie to naznačuje, že jednoduchý lineárny model v loge postačuje; prípadné nelinearity nijako nemenia kvalitatívny záver o negatívnej asociácii (P)HEV a EPI po zohľadnení HDP.

Zhluková analýza (cluster analysis)

Cieľom je rozdeliť krajiny do homogénnych skupín (klastrov) podľa podobnosti v ukazovateľoch kvality prostredia a rozvoja.

Zaujíma nás: - ktoré krajiny sa správajú podobne z hľadiska EPI, rozšírenia elektromobility a HDP, - ako sa klastre líšia (priemery EPI, (P)HEV, HDP – tzv. centroidy), - či je rozumné mať 2, 3 alebo 4 klastre (pomôžeme si silhouette a dendrogramom).

Poznámka (limit dát):
Ideálne by bolo klastrovať ukazovatele na obyvateľa (napr. (P)HEV na milión obyvateľov), aby sme odstránili vplyv veľkosti krajiny. V mojom datasete však chýbajú údaje o populácii, takže takýto prepočet neviem urobiť. Ako kompromis preto používam: - log-transformáciu počtových premenných (log1p(BEV_PHEV) resp. log1p(BEV)),
- následnú štandardizáciu (z-skóre) pri výpočte vzdialeností.

Tým tlmím vplyv extrémov a rozdielnych mierok, aj keď úplne neodstraňujem fakt, že väčšie krajiny majú prirodzene viac vozidiel.

Budeme klastrovať krajiny v jednom vybranom roku (prierez).

Balíčky a načítanie dát

library(tidyverse)
library(cluster)      # silhouette
library(factoextra)   # fviz_* helpery (dendrogram, atď.)
library(stats)        # hclust, dist
library(kableExtra)   # pekné tabuľky

# ak by nebol objekt `udaje` dostupný z úvodu, načítam ho:
if (!exists("udaje")) {
  udaje <- read.csv("data_r_comma_utf8.csv", header = TRUE, check.names = FALSE)
}

1)Konfigurácia

# Rok, v ktorom klastrujeme krajiny – zoberiem najnovší rok v dátach
ROK <- max(udaje$Rok, na.rm = TRUE)

# Pracovné premenné:
use_bev_phev <- TRUE   # TRUE = BEV_PHEV, FALSE = BEV

if (use_bev_phev) {
  feature_vars <- c("EPI", "BEV_PHEV", "HDP")
} else {
  feature_vars <- c("EPI", "BEV", "HDP")
}

id_vars <- c("Krajina", "Rok")   # identifikátory

2)Filtrovanie na prierez (rok) a príprava vstupov

# Prierez: krajiny v zvolenom roku
df <- udaje %>%
  filter(Rok == ROK) %>% 
  select(all_of(c(id_vars, feature_vars))) %>%
  drop_na()

# Prehľad prvých riadkov
df %>% 
  head(10) %>% 
  kbl(caption = paste0("Náhľad dát pre rok ", ROK)) %>% 
  kable_classic()

Náhľad dát pre rok 2023
Krajina	Rok	EPI	BEV_PHEV	HDP
CZ	2023	60.70	33448.0	21660
SK	2023	69.05	15041.0	18750
MR	2023	63.70	66508.0	16060
PL	2023	66.30	98348.0	15950
EU-avrg	2023	70.30	290224.4	33280

3)Transformácie (log1p pre počty, štandardizácia)

# Log1p na počtovej prem. (P)HEV resp. BEV
if (use_bev_phev) {
  df <- df %>% mutate(logBEVP = log1p(BEV_PHEV))
} else {
  df <- df %>% mutate(logBEV = log1p(BEV))
}

# Vstupné stĺpce do klastrovania:
X <- if (use_bev_phev) {
  df %>% select(EPI, logBEVP, HDP)
} else {
  df %>% select(EPI, logBEV, HDP)
}

# Štandardizácia (z-skóre)
X_scaled <- scale(X)

# Euklidovská vzdialenosť
D <- dist(X_scaled, method = "euclidean")

4)Ward.D2 hierarchické zhlukovanie + vizualizácia

hc <- hclust(D, method = "ward.D2")

# Dendrogram
fviz_dend(
  hc,
  k = NULL,
  cex = 0.9,
  main = paste("Dendrogram – rok", ROK),
  color_labels_by_k = FALSE,
  rect = FALSE
)

5)Pomôcky na výber počtu klastrov (k)

# Krok 5: Silhouette + výber počtu klastrov (rovnaké dáta ako dendrogram)

n      <- nrow(X_scaled)
k_grid <- 2:min(10, n - 1)

sil_tbl <- purrr::map_dfr(k_grid, function(k) {
  grp <- cutree(hc, k = k)
  sil <- cluster::silhouette(grp, D)
  tibble::tibble(
    k        = k,
    mean_sil = mean(sil[, "sil_width"], na.rm = TRUE)
  )
})

k_best <- sil_tbl$k[which.max(sil_tbl$mean_sil)]

print(sil_tbl)

## # A tibble: 3 × 2
##       k mean_sil
##   <int>    <dbl>
## 1     2    0.386
## 2     3    0.217
## 3     4    0.200

message(sprintf("Vybrané k = %d (max. priemerná silhouette = %.3f)",
                k_best, max(sil_tbl$mean_sil, na.rm = TRUE)))

ggplot(sil_tbl, aes(k, mean_sil)) +
  geom_line() + 
  geom_point() +
  labs(
    x = "Počet klastrov (k)",
    y = "Priemerná silhouette",
    title = "Silhouette podľa k (Ward.D2)"
  ) +
  theme_minimal()

6,7)Priradenie do klastrov a dendrogram s rezom

# Predpokladáme, že už existujú:
# - X_scaled, D, hc  (z krokov 3–4)
# - k_best           (z kroku 5 – silhouette)

# 1) Priradenie krajín do klastrov pre zvolené k_best
cl_best <- cutree(hc, k = k_best)

# tabuľka počtu krajín v jednotlivých klastroch
table(cl_best)

## cl_best
## 1 2 
## 4 1

# 2) Dendrogram s rámčekmi okolo klastrov
plot(
  hc,
  labels = rownames(X_scaled),   # menovky = krajiny
  main   = sprintf("Ward.D2 – rez pri k = %s (rok %s)", k_best, ROK),
  xlab   = "", 
  sub    = "", 
  cex    = 0.8
)
rect.hclust(hc, k = k_best, border = "#e41a1c")

8)Profilovanie a prezentácia klastrov

library(dplyr)
library(tidyr)
library(ggplot2)
library(kableExtra)
library(tibble)

## 1) Veľkosť klastrov
sizes <- as.data.frame(table(Klaster = factor(cl_best)))

kbl(sizes, caption = "Veľkosť klastrov") %>%
  kable_classic()

Veľkosť klastrov
Klaster	Freq
1	4
2	1

ggplot(sizes, aes(Klaster, Freq)) +
  geom_col() +
  geom_text(aes(label = Freq), vjust = -0.3) +
  labs(
    x = "Klaster",
    y = "Počet krajín",
    title = "Veľkosť klastrov"
  ) +
  theme_minimal()

## 2) Priradenie krajín do klastrov
members <- tibble(
  Krajina = df$Krajina,          
  Klaster = factor(cl_best)
) %>%
  arrange(Klaster, Krajina)

kbl(members, caption = "Priradenie krajín do klastrov", booktabs = TRUE) %>%
  kable_classic(full_width = FALSE)

Priradenie krajín do klastrov
Krajina	Klaster
CZ	1
MR	1
PL	1
SK	1
EU-avrg	2

## 3) Centroidy v z-skóre (ľahko porovnateľné)
centroids_z <- as.data.frame(X_scaled) %>%
  mutate(Klaster = factor(cl_best)) %>%
  group_by(Klaster) %>%
  summarise(across(where(is.numeric), mean), .groups = "drop")

kbl(centroids_z, digits = 2,
    caption = "Centroidy v z-skóre (EPI, log(BEV/BEV_PHEV), HDP)") %>%
  kable_classic()

Centroidy v z-skóre (EPI, log(BEV/BEV_PHEV), HDP)
Klaster	EPI	logBEVP	HDP
1	-0.27	-0.34	-0.42
2	1.10	1.37	1.69

## 4) Centroidy v pôvodných mierkach
df_with_cl <- df %>%
  mutate(Klaster = factor(cl_best))

centroids_raw <- df_with_cl %>%
  group_by(Klaster) %>%
  summarise(
    across(where(is.numeric), ~ mean(.x, na.rm = TRUE)),
    .groups = "drop"
  )

kbl(centroids_raw, digits = 2,
    caption = "Centroidy v pôvodných mierkach (EPI, BEV/BEV_PHEV, HDP)") %>%
  kable_classic()

Centroidy v pôvodných mierkach (EPI, BEV/BEV_PHEV, HDP)
Klaster	Rok	EPI	BEV_PHEV	HDP	logBEVP
1	2023	64.94	53336.25	18105	10.66
2	2023	70.30	290224.40	33280	12.58

## 5) Heatmapa centroidov (z-skóre)
heat <- centroids_z %>%
  pivot_longer(-Klaster, names_to = "Premenná", values_to = "z_mean")

ggplot(heat, aes(Premenná, Klaster, fill = z_mean)) +
  geom_tile() +
  geom_text(aes(label = round(z_mean, 2)), size = 3) +
  scale_fill_gradient2(
    low = "#4575b4",
    mid = "white",
    high = "#d73027",
    midpoint = 0
  ) +
  labs(
    title = "Centroidy (z-skóre) – heatmapa",
    x = NULL,
    y = NULL,
    fill = "z-priemer"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

## 6) Stručný prehľad krajín v klastroch
members %>%
  group_by(Klaster) %>%
  summarise(Krajiny = paste(Krajina, collapse = ", "), .groups = "drop") %>%
  kbl(caption = "Krajiny v jednotlivých klastroch") %>%
  kable_classic()

Krajiny v jednotlivých klastroch
Klaster	Krajiny
1	CZ, MR, PL, SK
2	EU-avrg

Zhlukovanie (Ward.D2) a silhouette graf ukázali, že najrozumnejšie je rozdeliť krajiny na dva klastre – pri k = 2 má priemerná silhouette najvyššiu hodnotu (≈ 0,39), pri 3 a 4 klastroch sa kvalita rozdelenia zhoršuje. Klaster 2 (1 krajina) má výrazne nadpriemerné EPI, HDP aj počet (P)HEV, ide teda o lídra s najvyššou environmentálnou výkonnosťou a zároveň najvyššou úrovňou elektromobility. Klaster 1 (4 krajiny) má naopak všetky tri ukazovatele pod priemerom vzorky, čo zodpovedá skupine menej bohatých krajín s nižším EPI aj nižšou adopciou (P)HEV. Zhluky sú preto nevyvážené (4 vs. 1 krajina) a výsledok do veľkej miery odráža veľkosť a ekonomickú silu krajín, keďže pracujeme s absolútnymi počtami vozidiel; pri detailnejšej politickej interpretácii by bolo vhodné použiť skôr normované ukazovatele (napr. (P)HEV na milión obyvateľov).

Multikolinearita

Úvod

Po autokorelácii a heteroskedasticite je multikolinearita tretím častým problémom lineárnej regresie.
V našom prípade sa týka najmä situácie, keď do modelu naraz zahrnieme:

log1p(BEV) – log počtu batériových elektromobilov,
log1p(BEV_PHEV) – log súčtu BEV a plug-in hybridov,
log1p(HDP) – log HDP na obyvateľa.

BEV a BEV_PHEV spolu prirodzene veľmi úzko súvisia (BEV je súčasťou BEV_PHEV), takže sú takmer lineárne závislé. To vedie k tomu, že matica \(\mathbf X^T \mathbf X\) je „takmer singulárna“ a jej inverzia je numericky nestabilná. Dôsledky:

Odhady \(\hat\beta\) sú nestabilné (malá zmena dát → veľká zmena koeficientov).
Štandardné chyby koeficientov sú nadhodnotené → regresory sa javia ako nevýznamné.
Interpretácia jednotlivých koeficientov (najmä BEV vs. BEV_PHEV) je veľmi problematická.

V ďalších krokoch si to ukážeme priamo na panelových dátach V4 (EPI, BEV, BEV_PHEV, HDP).

Východiskový model a pracovný dataset

Pre účely multikolinearity si zámerne zostavíme „problémový“ model, kde dáme naraz log1p(BEV) aj log1p(BEV_PHEV):

\[ EPI_{it} = \beta_0 + \beta_1 \log(1 + BEV_{it}) + \beta_2 \log(1 + BEV\_PHEV_{it}) + \beta_3 \log(1 + HDP_{it}) + u_{it} \]

# predpoklad: objekt `udaje` je už načítaný vyššie z data_r_comma_utf8.csv

library(dplyr)

udaje_mc <- udaje %>%
  dplyr::select(Krajina, Rok, EPI, BEV, BEV_PHEV, HDP) %>%
  mutate(
    logBEV       = log1p(BEV),
    logBEV_PHEV  = log1p(BEV_PHEV),
    logHDP       = log1p(HDP)
  ) %>%
  drop_na()

# základný "multikolineárny" model
mod_mc <- lm(EPI ~ logBEV + logBEV_PHEV + logHDP, data = udaje_mc)
summary(mod_mc)

## 
## Call:
## lm(formula = EPI ~ logBEV + logBEV_PHEV + logHDP, data = udaje_mc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.2169  -5.0787  -0.3093   4.7756  15.2159 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -31.059     35.722  -0.869   0.3891   
## logBEV         7.369      4.903   1.503   0.1397   
## logBEV_PHEV  -11.173      4.939  -2.262   0.0285 * 
## logHDP        13.830      3.971   3.483   0.0011 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.055 on 46 degrees of freedom
## Multiple R-squared:  0.4456, Adjusted R-squared:  0.4094 
## F-statistic: 12.32 on 3 and 46 DF,  p-value: 4.83e-06

V základnom modeli, kde súčasne používam log-transformované počty BEV aj (P)HEV, vychádza:

logHDP má očakávaný kladný a štatisticky významný vzťah k EPI,
logBEV_PHEV je záporný a významný,
samotný logBEV je kladný, ale štatisticky nevýznamný.

Model má relatívne slušné \(R^2 \approx 0{,}45\), ale dve veľmi podobné premenné (logBEV a logBEV_PHEV) majú odlišné znamienka a jedna z nich „spadne“ zo štatistickej významnosti. To je typický varovný signál multikolinearity – model nedokáže spoľahlivo rozlíšiť, ktorá z týchto dvoch úzko prepojených premenných nesie informáciu.

Korelačná matica a scatterplotová matica

xvars <- udaje_mc %>%
  dplyr::select(logBEV, logBEV_PHEV, logHDP)

round(cor(xvars), 3)

##             logBEV logBEV_PHEV logHDP
## logBEV       1.000       0.994  0.564
## logBEV_PHEV  0.994       1.000  0.582
## logHDP       0.564       0.582  1.000

pairs(xvars,
      main = "Scatterplotová matica – logBEV, logBEV_PHEV, logHDP")

Korelačná matica ukazuje, že medzi premennými logBEV a logBEV_PHEV je korelácia 0,994, teda prakticky dokonalá lineárna závislosť. To znamená, že počet BEV a (P)HEV spolu takmer dokonale spolubežia – informácia, ktorú nesú, je veľmi podobná. Korelácie s logHDP sú mierne až stredne vysoké (0,56–0,58), ale zďaleka nie také extrémne ako medzi oboma EV premennými.

Scatterplotová matica tento výsledok vizuálne potvrdzuje. V paneli logBEV vs. logBEV_PHEV body ležia takmer na jednej priamke – ide teda o takmer perfektný lineárny vzťah. Vzťahy logBEV – logHDP a logBEV_PHEV – logHDP sú zreteľne pozitívne, ale rozptýlené; medzi týmito dvojicami premenných nie je úplná redundancia. Z hľadiska multikolinearity je preto problémom najmä dvojica logBEV a logBEV_PHEV, ktoré do modelu prinášajú veľmi podobnú informáciu.

VIF

library(car)
vif(mod_mc)

##      logBEV logBEV_PHEV      logHDP 
##   85.177825   87.911765    1.556331

X   <- model.matrix(mod_mc)[, -1]   # bez interceptu
XtX <- t(X) %*% X
eig <- eigen(XtX)

condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number

## [1] 108.0268

VIF a Condition number – interpretácia

Hodnoty VIF jasne ukazujú extrémnu multikolinearitu medzi logBEV a logBEV_PHEV.
Obe premenné majú VIF približne 85–88, čo vysoko presahuje bežne používané orientačné prahy (5 alebo 10). To znamená, že variancia odhadov ich koeficientov je v porovnaní so situáciou bez multikolinearity nafúknutá približne 80-násobne. Naproti tomu logHDP má VIF ≈ 1,56, teda prakticky žiadny problém. Výsledok potvrdzuje, že problém multikolinearity v našom modeli spôsobuje najmä dvojica premenných logBEV a logBEV_PHEV, ktoré nesú takmer identickú informáciu o rozšírení elektromobility.

Tento záver dopĺňa aj condition number. Vypočítaná hodnota vychádza κ ≈ 108, čo podľa bežného pravidla (hodnoty nad 30–100 už považujeme za vážny problém) znamená silnú až veľmi výraznú multikolinearitu v matici \(X'X\). Matica je zle kondicionovaná – malé zmeny v dátach by viedli k veľkým zmenám v odhadoch regresných koeficientov. V kombinácii s takmer jednotkovou koreláciou medzi logBEV a logBEV_PHEV a extrémne vysokými VIF to potvrdzuje, že spoločné zaradenie týchto dvoch premenných do jedného modelu je štatisticky problematické a vedie k nestabilným a ťažko interpretovateľným odhadom.

Riešenie: vynechanie premennej

# model len s (P)HEV + HDP
mod_no_logBEV <- lm(EPI ~ logBEV_PHEV + logHDP, data = udaje_mc)
summary(mod_no_logBEV)

## 
## Call:
## lm(formula = EPI ~ logBEV_PHEV + logHDP, data = udaje_mc)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.005  -4.389  -1.012   4.244  14.395 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -24.7863    35.9496  -0.689  0.49392    
## logBEV_PHEV  -3.8147     0.6564  -5.811 5.19e-07 ***
## logHDP       12.8282     3.9667   3.234  0.00224 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.149 on 47 degrees of freedom
## Multiple R-squared:  0.4184, Adjusted R-squared:  0.3936 
## F-statistic:  16.9 on 2 and 47 DF,  p-value: 2.947e-06

library(car)
vif(mod_no_logBEV)

## logBEV_PHEV      logHDP 
##    1.512513    1.512513

# model len s BEV + HDP
mod_no_logBEV_PHEV <- lm(EPI ~ logBEV + logHDP, data = udaje_mc)
summary(mod_no_logBEV_PHEV)

## 
## Call:
## lm(formula = EPI ~ logBEV + logHDP, data = udaje_mc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.1351  -4.4467  -0.8368   4.9316  15.2734 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -16.7362    36.6646  -0.456  0.65016    
## logBEV       -3.6278     0.6707  -5.409 2.09e-06 ***
## logHDP       11.6591     4.0185   2.901  0.00564 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.358 on 47 degrees of freedom
## Multiple R-squared:  0.3839, Adjusted R-squared:  0.3577 
## F-statistic: 14.64 on 2 and 47 DF,  p-value: 1.14e-05

vif(mod_no_logBEV_PHEV)

##   logBEV   logHDP 
## 1.465476 1.465476

Po odstránení jednej z dvojice silno korelovaných premenných (logBEV, logBEV_PHEV) sa správanie modelu výrazne zlepšilo. V oboch zjednodušených špecifikáciách –

EPI ~ logBEV_PHEV + logHDP,
EPI ~ logBEV + logHDP –

vychádzajú oba regresory štatisticky významné: log počtu (P)HEV aj log počtu BEV má negatívny a logHDP pozitívny vplyv na EPI. VIF sa znižujú na hodnoty okolo 1,5, takže multikolinearita už nie je problémom. Upravený koeficient determinácie zostáva porovnateľný s pôvodným modelom s oboma EV premennými, pričom mierne lepšie (vyššie \(R^2\)) vychádza špecifikácia s logBEV_PHEV.

Z praktického hľadiska to znamená, že v dátach V4 môžem spoľahlivo odhadovať vzťah medzi EPI a celkovým rozšírením elektromobility (či už meraným BEV alebo (P)HEV) po zohľadnení HDP, ale nie je rozumné mať BEV a BEV_PHEV v jednom modeli súčasne. Preto v ďalšej analýze používam oddelené modely s BEV a s BEV_PHEV, čím sa vyhýbam silnej multikolinearite a získavam stabilnejšie a lepšie interpretovateľné odhady.

PCA analýza

# PCA len z logBEV a logBEV_PHEV (štandardizované)
X_pca   <- scale(udaje_mc[, c("logBEV", "logBEV_PHEV")],
                 center = TRUE, scale = TRUE)
pca_res <- prcomp(X_pca)

summary(pca_res)   # podiel variability vysvetlený PC1 a PC2

## Importance of components:
##                          PC1     PC2
## Standard deviation     1.412 0.07784
## Proportion of Variance 0.997 0.00303
## Cumulative Proportion  0.997 1.00000

# prvá hlavná zložka ako index elektromobility
udaje_mc$PC1_EV <- pca_res$x[, 1]

mod_pca <- lm(EPI ~ PC1_EV + logHDP, data = udaje_mc)
summary(mod_pca)

## 
## Call:
## lm(formula = EPI ~ PC1_EV + logHDP, data = udaje_mc)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.983  -4.378  -1.034   4.376  14.839 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -52.852     39.227  -1.347  0.18433    
## PC1_EV         5.034      0.895   5.624 9.92e-07 ***
## logHDP        12.288      3.991   3.079  0.00346 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.246 on 47 degrees of freedom
## Multiple R-squared:  0.4025, Adjusted R-squared:  0.3771 
## F-statistic: 15.83 on 2 and 47 DF,  p-value: 5.536e-06

library(car)
vif(mod_pca)

##   PC1_EV   logHDP 
## 1.490645 1.490645

Pri PCA nad dvojicou silno korelovaných premenných logBEV a logBEV_PHEV vychádza, že prvá hlavná zložka PC1 vysvetľuje približne 99,7 % variability pôvodných premenných, zatiaľ čo druhá zložka má zanedbateľný príspevok. Z toho dôvodu môžem PC1 (v dokumente ju označujem ako PC1_EV) interpretovať ako index celkovej úrovne elektromobility v krajine.

V modeli EPI ~ PC1_EV + logHDP má PC1_EV štatisticky významný vplyv (p < 0,001) a logHDP ostáva pozitívny a významný. Upravené \(R^2\) (~0,38–0,40) je veľmi podobné ako pri modeloch s pôvodnými premennými logBEV/logBEV_PHEV a hodnoty VIF sú nízke (≈1,5), takže multikolinearita už nie je problémom.

PCA teda potvrdzuje, že informáciu o počte BEV a (P)HEV je možné zhrnúť do jedného faktora a tým technicky vyriešiť multikolinearitu. Nevýhodou je však slabšia vecná interpretácia koeficientu pri PC1_EV (je to zmena v abstraktnom indexe, nie „percentá vozidiel“), preto v hlavnej analýze pracujem s jednoduchšími oddelenými modelmi s BEV a s BEV_PHEV a PCA uvádzam len ako doplnkový ilustratívny prístup.

library(factoextra)

# voliteľné: menovky podľa krajiny a roku 
rownames(pca_res$x) <- paste(udaje_mc$Krajina, udaje_mc$Rok, sep = " ")

fviz_pca_biplot(
  pca_res,
  geom.ind  = "point",      # pozorovania ako body
  col.ind   = "#075E54",    # smaragdovo zelená, ladí s témou
  pointshape = 19,
  pointsize  = 2,
  label     = "var",        # menovky len pre premenné (šípky)
  col.var   = "#D55E00",    # farba šípok (logBEV, logBEV_PHEV)
  repel     = TRUE          # odpudzovanie menoviek, aby sa neprekrývali
) +
  labs(
    title = "PCA biplot – logBEV a logBEV_PHEV",
    x = "PC1",
    y = "PC2"
  ) +
  theme_minimal()

Záver práce Vysvetlenie hypotéz

library(tibble)
library(kableExtra)

hyp_tbl <- tibble::tribble(
  ~Hypoteza, ~Formalne_znenie,                               ~Vysledok,                                                         ~Rozhodnutie,
  "H0",      "BEV/(P)HEV nemá významný vplyv na EPI.",       "Koeficienty pri log1p(BEV) aj log1p(BEV_PHEV) sú významné.",      "H0 zamietam.",
  "H1",      "Vyšší BEV/(P)HEV je spojený s vyšším EPI.",    "Odhadnuté koeficienty sú záporné → viac EV = nižší EPI.",          "H1 nepodporujem (pozitívny vplyv sa nepotvrdil).",
  "Doplnenie",
             "Vyšší HDP na obyvateľa zvyšuje EPI.",
             "Koeficient pri log1p(HDP) je stabilne kladný a štatisticky významný.",
             "Očakávanie potvrdzujem."
)

kbl(
  hyp_tbl,
  caption = "Zhrnutie hypotéz a výsledkov",
  col.names = c("Hypotéza", "Formálne znenie", "Výsledok v dátach V4", "Rozhodnutie")
) |>
  kable_classic(full_width = FALSE)

Zhrnutie hypotéz a výsledkov
Hypotéza	Formálne znenie	Výsledok v dátach V4	Rozhodnutie
H0	BEV/(P)HEV nemá významný vplyv na EPI.	Koeficienty pri log1p(BEV) aj log1p(BEV_PHEV) sú významné.	H0 zamietam.
H1	Vyšší BEV/(P)HEV je spojený s vyšším EPI.	Odhadnuté koeficienty sú záporné → viac EV = nižší EPI.	H1 nepodporujem (pozitívny vplyv sa nepotvrdil).
Doplnenie	Vyšší HDP na obyvateľa zvyšuje EPI.	Koeficient pri log1p(HDP) je stabilne kladný a štatisticky významný.	Očakávanie potvrdzujem.

Analýza dát o elektromobilite

Simona Vančová

Úvod

Charakteristika dát

Premenné

Výskumná otázka a hypotézy

Metodika v skratke

Import údajov z mojej bakalárskej práce

Interaktívna tabuľka s dátami

Deskriptívna štatistika

Prieskumné grafy - BOXPLOT

Lineárna regresia (dve alternatívy)

Diagnostika (krajšie grafy)

Residuals vs Fitted

Normal Q–Q

Scale–Location (√|štandardizované rezíduá| vs. fitted)

Residuals vs Leverage

Vplyvné pozorovania

Oneskorený vplyv (lag 1 rok) s fixnými efektmi

Štandardizované koeficienty

Heteroskedasticita – testy a grafy

Vizuálna diagnostika

Heteroskedasticita - prehľad testov

Formálne testy

Nelineárne špecifikácie a testy funkčnej formy

1) Ramsey RESET – kontrola špecifikácie

2) Component + Residual (C+R) grafy – kde sa krivka ohýba

3) Kvadratické členy – jemné zakrivenie

4) Zlom v sklone – iný efekt (P)HEV pri vyšších hodnotách

5) Výber pracovného modelu a krátke zhrnutie

Porovnávacia tabuľka modelov

Graf predikčných kriviek (Lineárny vs. Kvadratický vs. Zlom)

Zhluková analýza (cluster analysis)

Multikolinearita

Úvod

Východiskový model a pracovný dataset

Korelačná matica a scatterplotová matica

VIF

VIF a Condition number – interpretácia

Riešenie: vynechanie premennej

PCA analýza