Analýza dát o elektromobilite

cat(’

’)

Import údajov z mojej bakalárskej práce

Načítam si csv súbor s dátami o elektromobilite vo V4. Sú to dáta, ktoré som používala na bakalársku prácu, je tam počet batériových elektrických vozidiel, hybridných vozidiel, EPI index, ktorý rozpráva o ekologickej vyspelosti krajín, emisie CO2 atď.Dáta už boli prečistené, tak v tomto zadaní nebolo treba ďalšie dočisťovanie.

data<-read.csv("data_r_comma_utf8.csv")

Úvod do problému, stanovenie hypotéz

Rozhodla som sa skúmať, či rozšírenie elektromobility súvisí s kvalitou životného prostredia meranou indexom environmentálnej výkonnosti (EPI). Budem pracovať s panelom krajín a rokov v mojom datasete (elektromobilita a súvisiace ukazovatele). Zaujíma ma, či vyšší počet registrovaných batériových elektromobilov (BEV) a BEV_PHEV (súčet BEV a plug-in hybridov) prispieva k lepšiemu skóre EPI, po zohľadnení hospodárskej úrovne a ďalších faktorov.Sledujeme krajiny V4. Teoreticky očakávame, že prechod na vozidlá s nižšími emisiami výfukových plynov, spolu s investíciami do infraštruktúry a zdrojov, zlepší kvalitu ovzdušia a environmentálne výsledky – čo by sa malo prejaviť vyšším EPI. Uvedomujeme si však, že EPI je zložený index (klíma, ovzdušie, voda, biodiverzita, odpadová politika…), preto do modelu zaradíme aj kontrolné premenné.

Závislá premenná:

EPI – index environmentálnej výkonnosti.

Kľúčové vysvetľujúce premenné:

BEV – počet batériových elektromobilov.

BEV_PHEV – počet BEV + plug-in hybridov.

Kontrolné premenné (podľa dostupnosti v dátach):

HDP (HDP na obyvateľa),

co2 (emise CO₂),

investičné/rozvojové ukazovatele: res_develop, r_sources, unvest_man, invest_sources, invest_transport

Hypotézy

Nulová hypotéza (H0): Počet BEV ani BEV_PHEV nemá štatisticky významný vplyv na EPI. Hlavná hypotéza (H1): Vyšší počet BEV aj BEV_PHEV je spojený s vyšším EPI (pozitívny vplyv).

if (!exists("udaje")) {
udaje <- read.csv("data_r_comma_utf8.csv", header = TRUE, check.names = FALSE)
}
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))

Boxploty ukazujú, že hodnota indexu EPI je v mojom súbore sústredená okolo 70 bodov a celkový rozptyl je skôr menší; objavuje sa však jeden nižší odľahlý prípad (≈50), na ktorý ďalej upozorňujem v diagnostike modelu. Po log-transformácii počtových premenných BEV a BEV_PHEV (log1p) sa rozdelenia výrazne stabilizovali – boxy sú symetrickejšie, bez extrémnych chvostov a s menšou vzdialenosťou od kvartilov, čo je priaznivé pre lineárne modelovanie. Premenná HDP má po logu veľmi úzky rozsah, takže v modeli bude pôsobiť najmä ako kontrola úrovne rozvoja, bez výraznych odľahlostí.

Lineárna regresia

# (1) Načítam data
udaje <- read.csv("data_r_comma_utf8.csv", header = TRUE, check.names = FALSE)

# (2) Lineárna regresia – dve alternatívy bez kolinearity (log1p)
# 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)

# (3) Súhrn výsledkov (obe)
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

# Voliteľne: robustné SE (HC1) pre oba modely
library(lmtest); library(sandwich)
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

# (4) Model pre diagnostické grafy 
# napr. ak chceš diagnostikovať PHEV model:
model <- m_phev

# (5) Diagnostické grafy pre vybraný model
par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))

NA
NA

Residuals vs Fitted (reziduá vs. predikované): body sú okolo nuly, ale krivka je jemne zakrivená. → Model nie je úplne priamka; môže mu chýbať nejaká premena alebo jej log/štvorec. Rozptyl sa mierne mení pri vyšších hodnotách, čiže môže ísť o slabú neustálu varianciu (heteroskedasticitu).

Q–Q plot: väčšina bodov sedí na priamke, konce sa však odchyľujú. → Reziduá nie sú dokonale normálne — v „chvostoch“ je viac extrémov, než by ideálne bolo. To väčšinou nevadí dramaticky, ale je dobré o tom vedieť.

Scale–Location: čiara nie je úplne vodorovná. → Potvrdzuje to, že rozptyl chýb nie je úplne rovnaký pre všetky predikované hodnoty (slabá heteroskedasticita).

Residuals vs Leverage: je tam jeden bod s veľkým vplyvom (označený napr. 41). → Tento záznam môže ťahať výsledky modelu. Oplatí sa skontrolovať, či je to chyba v dátach alebo legitímne extrémny prípad.

Krajšie typy grafov pre lineárnu regresiu s použitím farby

library(broom)
library(dplyr)
library(ggplot2)
library(patchwork)

clr <- "#0E6B4D"

# diagnostický dataframe + vlastné ID
diag <- augment(model) %>%
  mutate(id = dplyr::row_number())   # naše popisky

# prah pre Cookovu vzdialenosť
cook_cut <- 4 / nrow(diag)
lab_df <- filter(diag, .cooksd > cook_cut)

# 1) Residuals vs Fitted
p1 <- ggplot(diag, 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"))

# 2) Normal Q–Q
p2 <- ggplot(diag, 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"))

# 3) Scale–Location
p3 <- ggplot(diag, 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"))

# 4) Residuals vs Leverage + popisky vplyvných bodov
p4 <- ggplot(diag, 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"))

# 2×2 panel
(p1 + p2) / (p3 + p4)

Inferenciu uvádzam s robustnými štandardnými chybami (HC1).

library(lmtest); library(sandwich)
coeftest(model, vcov = vcovHC(model, 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

cat("\nPozn.: V zátvorkách sú robustné štandardné chyby (HC1). * p<0,10; ** p<0,05; *** p<0,01.\n")

Pozn.: V zátvorkách sú robustné štandardné chyby (HC1). * p<0,10; ** p<0,05; *** p<0,01.

Použila som robustné štandardné chyby (HC1), lebo grafy ukázali, že chyby nemajú všade rovnaký rozptyl a pár bodov má väčší vplyv. Bežné chyby by v takom prípade podhodnotili neistotu a p-hodnoty by boli príliš optimistické. Robustné chyby preto dávajú spoľahlivejšie p-hodnoty a intervaly pri tých istých koeficientoch.

Q-Q plot

Q–Q plot porovnáva kvantily štandardizovaných rezíduí s teoretickými normálnymi kvantilmi. Väčšina bodov leží pri priamke, no na koncoch sú mierne odchýlky, čo naznačuje ľahkú odchýlku od normality v chvostoch. Predpoklad normality je približne splnený; pre istotu uvádzam aj robustné smerodajné chyby.

# Q–Q plot rezíduí (samostatne)

qqnorm(rstandard(model)); qqline(rstandard(model))

NA
NA

Na Q–Q grafe vidno, že väčšina bodov leží blízko referenčnej priamky, takže rozdelenie rezíduí je približne normálne v strede. Na oboch koncoch sa body od priamky odchýlia – vľavo dole sú o niečo nižšie a vpravo hore vyššie než by boli pri dokonalej normalite. To znamená, že v chvostoch je pár extrémnejších hodnôt (mierne ťažšie chvosty), ale celkový tvar je inak veľmi blízky normalite.

plot(model, which = 2)  # which=2 = Q–Q plot

Na Q–Q grafe vidno, že väčšina boodov leží blízko priamky, takže moje reziduá sa správajú približne normálne .Menšie odchýlky sú na oboch koncoch grafu – vľavo dole a vpravo hore – čo znamená, že v chvostoch mám pár extrémnejších hodnôt.

Scale-Location plot

Body sú rozmiestnené pomerne rovnomerne naprieč osou X, bez výrazneho lievika.Červená LOESS krivka je len jemne naklonená – variancia rezíduí sa mierne mení, preto reportujem robustné štandardné chyby

Beriem to tak, že predpoklad konštantnej variability je približne splnený; pre istotu však uvádzam aj výsledky s robustnými smerodajnými chybami. Niekoľko bodov je mierne nad hodnotou 1,5 na osi Y, ale nejde o extrémne anomálie.

if (!exists("model")) {
if (!exists("udaje")) {
udaje <- read.csv("data_r_comma_utf8.csv", header = TRUE, check.names = FALSE)
}
model <- lm(EPI ~ log1p(BEV) + log1p(BEV_PHEV) + HDP, data = udaje)
}

# Scale–Location plot (which = 3)

plot(model, which = 3)

Graf Residuals vs Leverage ukazuje, aký veľký vplyv môžu mať jednotlivé pozorovania na priamku. Väčšina bodov je pri nízkom leverage a má malé reziduá, takže väčšina záznamov model neťahá.Vidieť však 1–2 body bližšie k bodkovaným krivkám Cookovej vzdialenosti – to sú potenciálne vplyvné pozorovania. Tie si kontrolujem testami a porovnám výsledky modelu s/bez nich.Záver: dáta sú vo všeobecnosti v poriadku; kvôli týmto bodom uvádzam aj robustné smerodajné chyby a robím doplnkovú kontrolu vplyvných bodov

# Residuals vs Leverage (s Cookovou vzdialenosťou)

plot(model, which = 5)

äčšina bodov má nízky leverage; zopár pozorovaní (20, 31, 41, 48) má vyšší leverage a/vyššie Cookovo D, teda môžu mať väčší vplyv na odhady. Hodnoty Cookovej vzdialenosti sú zvýšené, no nie extrémne.

library(car)
library(lmtest); library(sandwich)

# Kombinácia leverage × rezíduum × Cook’s D (označí najpodozrivejšie body)

influencePlot(model, main = "Influence plot")

# Najsilnejší outlier (Bonferroni p-value)

outlierTest(model)

No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:

# Robustné smerodajné chyby pre istotu

coeftest(model, vcov = vcovHC(model, 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

Bonferroni test neidentifikoval žiadny štatisticky významný outlier (p ≥ 0,05), takže problémom nie je jeden „chybný“ bod s extrémnym reziduom, skôr kombinácia leverage a rezíduí pri pár pozorovaniach. Pri log-špecifikácii vychádza HDP pozitívny a štatisticky významný, čo je v súlade s očakávaním (vyššia ekonomická úroveň súvisí s lepším EPI). Premenná log1p(BEV) nie je štatisticky významná – nulovú hypotézu o nulovom efekte nedokážeme zamietnuť. V tomto modeli teda nemáme dôkaz, že počet BEV súvisí s EPI.

# Model A: EPI ~ log(BEV) + HDP
m_bev  <- lm(EPI ~ log1p(BEV) + HDP, data = udaje)
summary(m_bev); coeftest(m_bev, vcov = vcovHC(m_bev, 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

# Model B: EPI ~ log(BEV_PHEV) + HDP
m_phev <- lm(EPI ~ log1p(BEV_PHEV) + HDP, data = udaje)
summary(m_phev); coeftest(m_phev, vcov = vcovHC(m_phev, 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

# porovnanie kvality
AIC(m_bev, m_phev)

NA

Vysvetlenie kontroly: Aby som sa vyhla kolinearite, odhadla som samostatne model s BEV a model s BEV_PHEV. Porovnanie AIC ukázalo, že lepšie sedí model s BEV_PHEV. V ňom je HDP pozitívne a významné, kým BEV_PHEV vychádza negatívne a vysoko významne. To je v rozpore s pôvodným očakávaním pozitívneho efektu elektromobility na EPI, a v diskusii uvádzam možné dôvody (napr. načasovanie – rýchly nárast PHEV bez okamžitého zlepšenia EPI, reálne emisie PHEV, iné environmentálne domény v EPI)

Načasovanie: oneskorený vplyv elektromobility (lag o 1 rok)

library(dplyr)

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

Štandardizované koeficienty = porovnateľná sila efektov

Uvidíme, ktorý prediktor má najväčší „zásah“ do EPI.

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

Najväčší „zásah“ má v tomto štandardizovanom modeli log1p(BEV_PHEV).

Po štandardizácii vychádza, že 1 SD nárast log(BEV_PHEV) súvisí s poklesom EPI o ~0.79 SD (p < 0.001), zatiaľ čo 1 SD nárast HDP súvisí s nárastom EPI o ~0.44 SD (p = 0.002); model vysvetľuje ~42 % variability EPI

library(dplyr)
library(kableExtra)

## 1) Slovné zhrnutie diagnostiky a výsledkov

t1 <- tibble(
`Q–Q plot` = "Väčšina bodov pri priamke; mierne odchýlky v chvostoch → normalita približne OK (reportujem robustné SE).",
`Scale-Location` = "LOESS jemne naklonený → slabá heteroskedasticita (riešená robustnými SE).",
`Residuals vs Leverage` = "Vplyvné body (20, 31, 35, 48); Bonferroni neoznačil extrémny outlier.",
`Výsledky regresie` = "Lepší je model s log1p(BEV_PHEV)+HDP (AIC ≈ 343 < 346). HDP: pozitívny a významný; log1p(BEV_PHEV): negatívny a vysoko významný; BEV nevýznamné."
)

emerald_light <- "#ECFDF5"  # pozadie buniek
emerald_head  <- "#D1FAE5"  # hlavička
emerald_text  <- "#065F46"  # tmavozelený text

kbl(t1, 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	Výsledky regresie
Väčšina bodov pri priamke; mierne odchýlky v chvostoch → normalita približne OK (reportujem robustné SE).	LOESS jemne naklonený → slabá heteroskedasticita (riešená robustnými SE).	Vplyvné body (20, 31, 35, 48); Bonferroni neoznačil extrémny outlier.	Lepší je model s log1p(BEV_PHEV)+HDP (AIC ≈ 343 < 346). HDP: pozitívny a významný; log1p(BEV_PHEV): negatívny a vysoko významný; BEV nevýznamné.

## 2) Rozhodnutia o hypotézach

t2 <- tibble::tribble(
  ~Premenná,  ~`H1 (očakávanie)`,            ~`Rozhodnutie o H0`,            ~`Rozhodnutie o H1`,                  ~Poznámka,
  "BEV",      "β_BE V > 0 (pozitívny)",       "ZAMIETAM (významné)",          "ZAMIETAM (efekt je záporný)",       "V modeli EPI ~ log1p(BEV)+log1p(HDP) je log1p(BEV) významne negatívny.",
  "BEV_PHEV", "β_PHEV > 0 (pozitívny)",       "ZAMIETAM (významné)",          "ZAMIETAM (efekt je záporný)",       "log1p(BEV_PHEV) významne negatívny; model s PHEV má nižší AIC."
)

kbl(t2, caption = "Hypotézy a rozhodnutia (robustné SE, log-špecifikácie)", 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(t2), background = emerald_light, color = emerald_text)

Hypotézy a rozhodnutia (robustné SE, log-špecifikácie)

Premenná	H1 (očakávanie)	Rozhodnutie o H0	Rozhodnutie o H1	Poznámka
BEV	β_BE V > 0 (pozitívny)	ZAMIETAM (významné)	ZAMIETAM (efekt je záporný)	V modeli EPI ~ log1p(BEV)+log1p(HDP) je log1p(BEV) významne negatívny.
BEV_PHEV	β_PHEV > 0 (pozitívny)	ZAMIETAM (významné)	ZAMIETAM (efekt je záporný)	log1p(BEV_PHEV) významne negatívny; model s PHEV má nižší AIC.

NA

Záver

Výsledky naznačujú, že HDP má na EPI pozitívny a štatisticky významný vplyv. V modeli s BEV (log-transformácia) je koeficient log1p(BEV) záporný a významný; v alternatívnom modeli s BEV_PHEV vychádza koeficient takisto záporný a významný. Podľa AIC dátam lepšie sedí verzia s BEV_PHEV. Štandardizované koeficienty ukazujú, že +1 SD log(BEV_PHEV) súvisí s ≈ −0,79 SD EPI (p < 0,001), zatiaľ čo +1 SD HDP s ≈ +0,44 SD EPI (p ≈ 0,002). Model vysvetľuje približne 40–42 % variability EPI. Diagnostika rezíduí ukázala len mierne odchýlky v chvostoch a slabú heteroskedasticitu; preto uvádzam robustné štandardné chyby (HC1). Niekoľko vplyvných bodov kvalitatívne závery nemení. Celkovo výsledky podporujú význam ekonomickej úrovne pre environmentálnu výkonnosť, no priame priaznivé prepojenie elektromobility s EPI sa v tomto súbore nepreukázalo; skôr vychádza negatívna asociácia. Možné dôvody: načasovanie efektov (prejav neskôr), reálne emisie PHEV a fakt, že EPI zachytáva viac domén než len dopravu.