Úvod

Tento dokument vychádza z prečistených údajov z mojej bakalárskej práce o elektromobilite v krajinách V4 a o vzťahu k indexu environmentálnej výkonnosti (EPI). Cieľom je transparentne ukázať: (i) prípravu a pochopenie dát; (ii) odhad jednoduchých, dobre interpretovateľných lineárnych modelov; (iii) dôslednú diagnostiku predpokladov a prácu s robustnými smerodajnými chybami; (iv) krátku diskusiu a limity zistení.

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.

Výskumná otázka a hypotézy

Skúmam, či väčšie rozšírenie elektromobility súvisí s vyšším EPI po zohľadnení ekonomickej úrovne.

  • H0: Počet BEV/BEV_PHEV nemá významný vplyv na EPI.

  • H1: Vyšší BEV/BEV_PHEV je spojený s vyšším EPI.

  • Závislá premenná: EPI

  • Kľúčové vysvetľujúce premenné: BEV (batériové EV), BEV_PHEV (BEV + plug‑in hybridy)

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é.

# Balíčky a dáta (použité ďalej v práci)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.0     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(broom)
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(sandwich)
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some
library(patchwork)
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
udaje <- read.csv("data_r_comma_utf8.csv", header = TRUE, check.names = FALSE)

Prieskumné grafy

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.

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.

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)
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

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

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 znamená, že rozptyl chýb v modeli nie je všade rovnaký (niekde sú chyby malé, inde viac „lietajú“). Ak by bola výrazná, bežné p-hodnoty by mohli klamať. Preto to overujeme vizuálne (Scale–Location a grafy „rezíduá² vs. prediktory“) a formálne (Breusch–Pagan/White štýl).

Postup: (i) vizuálne diagnostiky; (ii) Breusch–Pagan / White‑štýl; (iii) reportovanie robustných SE. Ak by bola heteroskedasticita výrazná a transformácie nepomohli, zvážila by som aj vážené OLS (WLS).

# Vizualizácia: vzťah rozptylu s predikciami a prediktormi
h_diag <- augment(m_phev, data = udaje) %>%
  mutate(res2 = .resid^2, absres = sqrt(abs(.resid)), x_bev = log1p(BEV_PHEV), x_hdp = HDP)

p_SL <- ggplot(h_diag, 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_diag, 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_diag, 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)
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

# Formálne testy
bptest(m_phev)
## 
##  studentized Breusch-Pagan test
## 
## data:  m_phev
## BP = 0.48764, df = 2, p-value = 0.7836
bptest(m_phev, varformula = ~ x_bev + x_hdp, data = h_diag)
## 
##  studentized Breusch-Pagan test
## 
## data:  m_phev
## BP = 0.29849, df = 2, p-value = 0.8614
bptest(m_phev, varformula = ~ .fitted + I(.fitted^2), data = h_diag)  # White-štýl
## 
##  studentized Breusch-Pagan test
## 
## data:  m_phev
## BP = 5.6585, df = 2, p-value = 0.05906
ncvTest(m_phev)
## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.003248773, Df = 1, p = 0.95455

Heteroskedasticita je slabá až stredná, najmä viazaná na HDP. Preto v práci reportujeme robustné smerodajné chyby (HC1) a kvalitatívne závery modelu sa tým nemenia.

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

sum_tbl <- tibble::tibble(
  `Q–Q plot`            = "Väčšina bodov na priamke; mierne odchýlky v chvostoch.",
  `Scale–Location`      = "Jemný sklon/U-tvar → slabá heteroskedasticita (riešené robustnými SE).",
  `Residuals vs Leverage` = "Niekoľko vplyvných bodov; bez extrémnych outlierov.",
  `Modely`              = "Lepší je model s log1p(BEV_PHEV)+log1p(HDP) (nižší AIC); HDP +, log1p(BEV_PHEV) –."
)

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

kbl(sum_tbl, caption = "Diagnostika a výsledky – 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:4, background = emerald_light, color = emerald_text, border_left = TRUE, border_right = TRUE)
Diagnostika a výsledky – stručné zhrnutie
Q–Q plot Scale–Location Residuals vs Leverage Modely
Väčšina bodov na priamke; mierne odchýlky v chvostoch. Jemný sklon/U-tvar → slabá heteroskedasticita (riešené robustnými SE). Niekoľko vplyvných bodov; bez extrémnych outlierov. Lepší je model s log1p(BEV_PHEV)+log1p(HDP) (nižší AIC); HDP +, log1p(BEV_PHEV) –. 
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é)",         "ZAMIETAM (efekt je záporný)", "log1p(BEV_PHEV) významne negatívny; nižší AIC."
)

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é) ZAMIETAM (efekt je záporný) log1p(BEV_PHEV) významne negatívny; nižší AIC.

Diskusia, limity a záver

Interpretácia: HDP vychádza stabilne pozitívny a významný. log1p(BEV_PHEV) má v dátach V4 prekvapivo negatívnu asociáciu s EPI, kým samotné BEV nevyšlo robustne pozitívne. Vysvetlení je viac: (i) načasovanie – EPI reaguje oneskorene; (ii) PHEV nemusia v reálnej prevádzke prinášať očakávané emisie; (iii) EPI je široký index (voda, biodiverzita, odpad…) a doprava tvorí len časť.

Limity: veľkosť vzorky pre V4, agregačná úroveň, prípadná meracia chyba v registroch vozidiel, absencia ďalších kontrol (energetický mix, priemyselné štruktúry).

Záver: Zatiaľ nevidím dôkaz o krátkodobom pozitívnom vzťahu medzi rozšírením (P)HEV a EPI v krajinách V4. Výsledky sú však informatívne a metodicky konzistentné; po doplnení vysvetľujúcich premenných a dlhšej časovej osi sa môže obraz zmeniť.

HDP je pozitívny a významný determinant EPI. log1p(BEV_PHEV) vychádza záporný a významný; BEV samo o sebe sa neukazuje ako pozitívne významné. Diagnostika: mierna heteroskedasticita a pár vplyvných bodov → inferenciu uvádzam s robustnými SE (HC1; prípadne klastrované podľa krajiny pri paneli). Priamy priaznivý vzťah elektromobility a EPI sa na tomto súbore nepreukázal; možné dôvody: oneskorené efekty, reálne emisie PHEV, EPI zachytáva širšie domény než len dopravu.