Ú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

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 (klaster 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, BEV_PHEV (alebo alternatívne BEV) 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 „lakeť testom“ a silhouette)

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

Metodika

Vstupné dáta: prečistený súbor data_r_comma_utf8.csv (z tvojej práce).

Použité premenné: EPI, BEV_PHEV (príp. BEV), HDP, plus identifikátory Krajina, Rok.

Log-transformácia počtových premenných: log1p(BEV_PHEV) (alebo log1p(BEV)) – obmedzí vplyv extrémov.

Štandardizácia (z-skóre) všetkých vstupných premenných pred výpočtom vzdialeností.

Metrika: euklidovská vzdialenosť; algoritmus: hierarchické zhlukovanie Ward.D2.

Výber počtu klastrov: graf „lakeť“ (WSS), silhouette a vizuálny dendrogram.

knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

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

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

# Dátový súbor z mojej práce:

udaje <- read.csv("data_r_comma_utf8.csv", header = TRUE, check.names = FALSE)

1)Konfigurácia

# Rok, v ktorom klastrujeme krajiny:

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")   # stĺpce s identifikáciou

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)

# Distance matrix

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) 
library(cluster); library(purrr); library(tibble); library(ggplot2)

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

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

sil_tbl <- map_dfr(k_grid, function(k){
  grp <- cutree(hc, k = k)
  sil <- silhouette(grp, D)
  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. 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

# Balíčky 
library(dplyr)
library(ggplot2)
library(purrr)
library(tibble)
library(cluster)

#  0) Príprava dát pre klastrovanie (aj pri knit v čistom prostredí) ---
# musíme mať načítaný data.frame 'udaje' (už ho máš z vyšších častí)
stopifnot(exists("udaje"))

# a) ak existuje X_scaled z predošlých krokov, použijeme ho
if (exists("X_scaled")) {
  udaje_scaled <- X_scaled
} else if (!exists("udaje_scaled")) {
  # b) inak si urobíme z-škálovanie z číselných stĺpcov
  num_df <- udaje %>% select(where(is.numeric)) %>% na.omit()
  stopifnot(nrow(num_df) > 2)
  udaje_scaled <- scale(num_df)
  rownames(udaje_scaled) <- rownames(num_df)
}
if (is.null(rownames(udaje_scaled))) rownames(udaje_scaled) <- seq_len(nrow(udaje_scaled))

# --- 1) Vzdialenosti a hierarchia Ward.D2 ---
D  <- dist(udaje_scaled, method = "euclidean")
hc <- hclust(D, method = "ward.D2")

# --- 2) Silhouette pre výber k (2 .. min(10, n-1)) ---
n      <- nrow(udaje_scaled)
k_grid <- 2:min(10, n - 1)

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

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

# Report + graf silhouette
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()

# --- 3) Priradenie do klastrov pre k_best ---
cl_best <- cutree(hc, k = k_best)
print(table(cl_best))

## cl_best
## 1 2 
## 4 1

# --- 4) Dendrogram (base R) + rámčeky okolo klastrov ---
plot(hc,
     labels = rownames(udaje_scaled),
     main = sprintf("Ward.D2 – rez pri k = %s", k_best),
     xlab = "", sub = "", cex = 0.8)
rect.hclust(hc, k = k_best, border = "#e41a1c")

8)Profilovanie a prezentácia klastrov

# ===== predpoklady =====
stopifnot(exists("udaje_scaled"), exists("cl_best"))

# 0) korektné "unscale" ak nemáme uložené pôvodné (neškálované) údaje
if (!exists("udaje_complete")) {
  center <- attr(udaje_scaled, "scaled:center")
  scalev <- attr(udaje_scaled, "scaled:scale")
  stopifnot(!is.null(center), !is.null(scalev))
  udaje_complete <- as.data.frame(
    sweep(sweep(udaje_scaled, 2, scalev, `*`), 2, center, `+`)
  )
  # zachováme názvy riadkov/ stĺpcov
  colnames(udaje_complete) <- colnames(udaje_scaled)
  rownames(udaje_complete) <- rownames(udaje_scaled)
}

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 
members <- tibble(
  Krajina = if (!is.null(rownames(udaje_scaled))) rownames(udaje_scaled) else seq_len(nrow(udaje_scaled)),
  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
1	1
2	1
3	1
4	1
5	2

# 3) Centroidy v z-skóre (ľahko porovnateľné)
centroids_z <- as.data.frame(udaje_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") %>% kable_classic()

Centroidy v z-skóre
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
centroids_raw <- udaje_complete %>%
  mutate(Klaster = factor(cl_best)) %>%
  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") %>%
  kable_classic()

Centroidy v pôvodných mierkach
Klaster	EPI	logBEVP	HDP
1	64.94	10.66	18105
2	70.30	12.58	33280

# 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))

members %>%
  dplyr::group_by(Klaster) %>%
  dplyr::summarise(Krajiny = paste(Krajina, collapse = ", ")) %>%
  kableExtra::kbl(caption = "Krajiny v klastroch") %>%
  kableExtra::kable_classic()

Krajiny v klastroch
Klaster	Krajiny
1	1, 2, 3, 4
2	5

Metodicky sme postupovali konzistentne: údaje sme štandardizovali (z-skóre), vzdialenosti sme rátali euklidovsky a hierarchické zhlukovanie sme robili metódou Ward.D2. Počet klastrov sme zvolili k = 2 na základe silhouette (max. priem. silueta = 0.386). Dendrogram s rezom tento výber potvrdzuje.

Rozdelenie do klastrov nie je extrémne nevyvážené: počty krajín podľa klastrov sú 4/1 (od klastru 1 po 2).

Profil klastrov (z-skóre, heatmapa):

Klaster 6 (≈4 krajiny) Výrazne nadpriemerné BEV a (P)HEV, vyšší HDP a R&D/inovačné ukazovatele, no EPI pod priemerom. Naznačuje to, že rýchla technologická a mobilitná expanzia sa do EPI krátkodobo nepremieta (oneskorenia efektov, dopravné externality).

Klaster 5 (≈10 krajín) EPI nad priemerom pri nižšom HDP a vyššom podiele OZE. „Enviro-efektívny“ profil: dobré výsledky v EPI bez nadpriemernej adopcie EV.

Klaster 4 (≈5 krajín) Nižší HDP, podpriemerná EV adopcia, zvýšené investície do dopravy a EPI mierne nad priemerom. Skôr fáza infraštruktúrneho dobiehania.

Klaster 3 (≈5 krajín) Mierne vyššie EV, nižší HDP aj EPI, vysoké investície v priemysle. Tranzitný profil s environmentálnymi výzvami.

Klaster 2 (≈10 krajín) Stredné hodnoty bez extrémov – „mainstream“ zhluk.

Klaster 1 (≈16 krajín) Najpočetnejší „baseline“ – väčšina premenných okolo priemeru.

Hlavné zistenie: Prierezovo nevidíme jednoznačne lineárne pozitívny vzťah medzi rozšírením (P)HEV a EPI. Objavujú sa dva režimy:

1)bohatšie/technologicky silné krajiny s rýchlou adopciou EV, kde sa prínosy do EPI krátkodobo nemusia ukázať,

2)krajiny s vyšším podielom OZE a nižším HDP, ktoré dosahujú lepšie EPI aj bez vysokej penetrácie EV.

Implikácie: pre zlepšovanie EPI nestačí len podpora EV – dôležité sú energetický mix, riadenie dopravného dopytu, mestská mobilita a časové oneskorenia (benefity EV sa do EPI premietajú postupne). Politikám dáva zhlukovanie rámec: pre každý klaster cielená kombinácia nástrojov (OZE vs. infra dobudovanie vs. manažment dopravy).

Limity: prierezový výber (jedno obdobie v časti analýzy), agregácia na úroveň krajín a obmedzený počet premenných. Pre robustné kauzálne tvrdenia by bolo vhodné pokračovať panelovým modelom s oneskoreniami a širšími kontrolami.

Stručná tabuľka výsledkov a záver

library(tibble)
library(kableExtra)

## 1) stručná diagnostika

sum_tbl <- tibble::tibble(
`Q–Q plot`                  = "Väčšina bodov na priamke; len mierne odchýlky v chvostoch.",
`Scale–Location`            = "Jemný U-tvar → slabá heteroskedasticita (riešime robustnými SE).",
`Residuals vs Leverage`     = "Pár vplyvných bodov, bez extrémov.",
`Heteroskedasticita – testy`= "Breusch–Pagan a NCV nevýznamné; White štýl hraničný → skôr slabý náznak.",
`RESET (špecifikácia)`      = "p ≈ 0.011 ⇒ lineárna špecifikácia má známky nesprávnej formy (skúšali sme jemné nelinearity)."
)

emerald_light <- "#ECFDF5"; emerald_head <- "#D1FAE5"; emerald_text <- "#065F46"

kbl(
sum_tbl,
caption = "Diagnostika – stručné zhrnutie",
booktabs = TRUE
) |>
kable_classic(full_width = TRUE, html_font = "Helvetica") |>
row_spec(0, background = emerald_head, color = emerald_text, bold = TRUE) |>
column_spec(1:ncol(sum_tbl), background = emerald_light, color = emerald_text, border_left = TRUE, border_right = TRUE)

Diagnostika – stručné zhrnutie
Q–Q plot	Scale–Location	Residuals vs Leverage	Heteroskedasticita – testy	RESET (špecifikácia)
Väčšina bodov na priamke; len mierne odchýlky v chvostoch.	Jemný U-tvar → slabá heteroskedasticita (riešime robustnými SE).	Pár vplyvných bodov, bez extrémov.	Breusch–Pagan a NCV nevýznamné; White štýl hraničný → skôr slabý náznak.	p ≈ 0.011 ⇒ lineárna špecifikácia má známky nesprávnej formy (skúšali sme jemné nelinearity).

## 2) porovnanie špecifikácií (AIC)

# predpoklad: objekty m_lin, m_quad, m_slope sú už odhadnuté v dokumente

aic_tbl <- tryCatch(
{
aa <- AIC(m_lin, m_quad, m_slope)
dplyr::rename(aa, Model = df, AIC = AIC) |>
dplyr::mutate(Model = rownames(aa))
},
error = function(e) tibble::tibble(
Model = c("Lineárny", "Kvadratický", "Zlom v sklone"),
AIC   = c(343.50, 345.03, 346.36)   # hodnoty z tvojich výstupov
)
)

kbl(aic_tbl, caption = "Porovnanie kvality modelov (nižší AIC = lepší)", booktabs = TRUE) |>
kable_classic(full_width = FALSE, html_font = "Helvetica") |>
row_spec(0, background = emerald_head, color = emerald_text, bold = TRUE) |>
column_spec(1:ncol(aic_tbl), background = emerald_light, color = emerald_text)

Porovnanie kvality modelov (nižší AIC = lepší)
	Model	AIC
m_lin	m_lin	343.4979
m_quad	m_quad	345.0346
m_slope	m_slope	346.3613

## 3) hypotézy

hyp_tbl <- tibble::tribble(
~Premenná,  ~`H1 (očakávanie)`,   ~`Rozhodnutie o H0`,               ~`Rozhodnutie o H1`,              ~Poznámka,
"BEV",      "β > 0",              "NEZAMIETAM / neisté",            "NEPODPORUJEME",                  "Efekt nevýznamný alebo záporný.",
"BEV_PHEV", "β > 0",              "ZAMIETAM (významné, záporné)",   "ZAMIETAM (pozitívny sa nepotvrdil)","log1p(BEV_PHEV) stabilne negatívny s HC1.",
"HDP",      "β > 0",              "ZAMIETAM (významné, pozitívne)", "PODPORUJEME",                    "HDP vychádza kladne a významne."
)

kbl(hyp_tbl, caption = "Hypotézy a rozhodnutia (robustné SE)", booktabs = TRUE) |>
kable_classic(full_width = FALSE, html_font = "Helvetica") |>
row_spec(0, background = emerald_head, color = emerald_text, bold = TRUE) |>
column_spec(1:ncol(hyp_tbl), background = emerald_light, color = emerald_text)

Hypotézy a rozhodnutia (robustné SE)
Premenná	H1 (očakávanie)	Rozhodnutie o H0	Rozhodnutie o H1	Poznámka
BEV	β > 0	NEZAMIETAM / neisté	NEPODPORUJEME	Efekt nevýznamný alebo záporný.
BEV_PHEV	β > 0	ZAMIETAM (významné, záporné)	ZAMIETAM (pozitívny sa nepotvrdil)	log1p(BEV_PHEV) stabilne negatívny s HC1.
HDP	β > 0	ZAMIETAM (významné, pozitívne)	PODPORUJEME	HDP vychádza kladne a významne.

Diskusia, limity a záver

Čo ukázala diagnostika a testy: Rozptyl chýb je najmä slabo nehomogénny, preto v celej práci reportujeme robustné smerodajné chyby (HC1). RESET test (p ≈ 0.011) signalizoval, že samotná priamka môže byť pre jednoduchý lineárny model tesná, preto sme skúsili jemné nelinearity (kvadrát log1p(BEV_PHEV), kvadrát log1p(HDP)) a aj zlom v sklone podľa úrovne (P)HEV.

Porovnanie špecifikácií: Napriek signálu z RESET-u sa kvalita modelov nezlepšila – AIC je najnižší pri lineárnom modeli a ANOVA porovnania kvadratickej/zlomovej špecifikácie voči lineárnej nevyšli významne. Predikčné krivky pre všetky tri verzie sú tvarovo podobné (mierne klesajúce), rozdiely sú malé a v rámci intervalov spoľahlivosti.

Interpretácia účinkov (robustné SE):

log1p(BEV_PHEV) je negatívny a štatisticky významný – v našom súbore sa s vyšším rozšírením (P)HEV spája nižší EPI.
HDP je pozitívny a významný – krajiny/roky s vyšším HDP dosahujú vyšší EPI.
Samotné BEV nevyšlo robustne pozitívne.

Zjednodušene: viac (P)HEV ≠ vyšší EPI v krátkom horizonte; EPI zrejme ovplyvňujú širšie faktory (ekonomická úroveň, štruktúra energií, zloženie dopravy atď.).

Limity: malá vzorka (V4 × roky), agregačná úroveň, možná meracia chyba v registroch vozidiel, chýbajúce kontroly (energetický mix, dopravné výkony, priemysel). Krátke obdobie môže skrývať oneskorené účinky.

Analýza dát o elektromobilite

Simona Vančová