Úvod

Táto štúdia skúma, či a ako je adopcia elektromobility v krajinách V4 spojená s kvalitou životného prostredia (EPI) a do akej miery vzťah súvisí s ekonomickou úrovňou (HDP/obyv.). Vychádzam z prečisteného panelového datasetu (krajina × rok) s hlavnými premennými EPI, BEV, BEV+PHEV a HDP/obyv.

Hlavnú asociáciu odhadujem pomocou OLS so semi-log špecifikáciou (log1p na počtoch vozidiel) a robustnými smerodajnými chybami (HC1); diagnostiku dopĺňam o RESET, testy heteroskedasticity a kontrolu vplyvných bodov. Keďže BEV a BEV+PHEV sú silno korelované, vyhodnocujem ich v oddelených modeloch a problém multikolinearity ilustrujem cez VIF, condition number a doplnkovo PCA. Nakoniec pomocou hierarchického zhlukovania (Ward.D2) profilujem krajiny podľa trojice EPI – (P)HEV – HDP. Limitáciou je absencia populácie, preto pracujem s absolútnymi počtami (zmiernené log-transformáciou a štandardizáciou).

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.

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.

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

1) Deskriptívna štatistika

Premenna	n	miss	mean	sd	q25	median	q75	min	max
Deskriptívna štatistika (číselné premenné)
Rok	50	0	2.02K	2.90	2.02K	2.02K	2.02K	2.01K	2.02K
BEV	50	0	15.93K	30.19K	911.00	4.16K	14.08K	113.00	163.25K
EPI	50	0	67.88	9.18	63.50	68.30	73.28	46.00	85.42
HDP	50	0	19.47K	6.66K	14.52K	17.32K	21.79K	11.36K	33.30K
BEV_PHEV	50	0	27.71K	54.79K	1.58K	7.20K	28.48K	201.00	290.22K
co2	50	0	734.50M	1.32B	54.59M	117.72M	129.41M	28.34M	3.57B
res_develop	50	0	1.52	0.48	1.07	1.45	1.93	0.79	2.31
r_sources	50	0	628.47	529.85	317.81	437.32	524.84	129.45	2.01K
unvest_man	50	0	0.15	0.20	0.05	0.08	0.11	0.03	1.00
invest_sources	50	0	0.15	0.07	0.09	0.14	0.19	0.04	0.35
invest_transport	50	0	0.27	0.13	0.17	0.26	0.32	0.09	0.72

Agregovaný stĺpcový graf ukazuje prudké zrýchlenie adopcie (P)HEV po roku 2018 a najmä po 2021, keď celkové počty vo V4 rastú takmer exponenciálne. Čiarový graf podľa krajín naznačuje, že Poľsko ťahá regionálny trend – po 2021 má jednoznačne najstrmšiu trajektóriu. Česko zrýchľuje od roku 2021, Slovensko rastie plynule z nižšej základne a Maďarsko vykazuje miernejší rast (s možným jednorazovým výkyvom okolo 2019).

Dôležité je, že ide o absolútne počty: poradie krajín preto odráža aj ich veľkosť a trh; pre „férové“ porovnanie by bolo vhodné doplniť ukazovatele na obyvateľa.

Celkovo grafy potvrdzujú, že V4 vstúpila do fázy rýchlej adopcie (P)HEV až v posledných rokoch, čo je dôležitý kontext pri interpretácii vzťahu k EPI a HDP v ďalších častiach analýzy.

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

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.

2) 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ášť:

Term	Estimate	Std_Error	t_value	p_value	CI_low	CI_high
Model A: EPI ~ log1p(BEV) + log1p(HDP)
robustné SE (HC1), 95% CI
(Intercept)	−16.736	37.345	−0.448	0.656	−89.933	56.461
log1p(BEV)	−3.628	0.762	−4.764	0.000	−5.120	−2.135	***
log1p(HDP)	11.659	4.112	2.835	0.007	3.599	19.719	**
Pozn.: * p<0.001, p<0.01, * p<0.05, · p<0.10.

Term	Estimate	Std_Error	t_value	p_value	CI_low	CI_high
Model B: EPI ~ log1p(BEV_PHEV) + log1p(HDP)
robustné SE (HC1), 95% CI
(Intercept)	−24.786	37.461	−0.662	0.511	−98.210	48.637
log1p(BEV_PHEV)	−3.815	0.779	−4.896	0.000	−5.342	−2.288	***
log1p(HDP)	12.828	4.159	3.084	0.003	4.676	20.980	**
Pozn.: * p<0.001, p<0.01, * p<0.05, · p<0.10.

Porovnanie kvality modelov (nižší AIC je lepší)
Model	df	AIC
m_bev	4	346.38
m_phev	4	343.50

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.

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

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

## 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 (robustné SE HC1)
Term	Estimate	Std. Error (HC1)	t value	p-hodnota	Signif
(Intercept)	86.8957	5.4334	15.9928	0.0000	***
log1p(BEV)	-3.5940	0.7567	-4.7497	0.0000	***
HDP	0.0005	0.0002	2.8144	0.0071	**

Model B — EPI ~ log1p(BEV_PHEV) + HDP (robustné SE HC1)
Term	Estimate	Std. Error (HC1)	t value	p-hodnota	Signif
(Intercept)	89.2739	5.8010	15.3893	0.0000	***
log1p(BEV_PHEV)	-3.8006	0.7855	-4.8384	0.0000	***
HDP	0.0006	0.0002	3.0789	0.0035	**

Porovnanie kvality modelov
Model	R²	Adj. R²	AIC	RMSE (resid. se)	n
A: EPI ~ log1p(BEV) + HDP	0.3823	0.3560	346.5049	7.3673	50
B: EPI ~ log1p(BEV_PHEV) + HDP	0.4190	0.3942	343.4470	7.1454	50

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.

Lagovaný model: EPI ~ log1p(lag_BEV_PHEV) + HDP + FE(krajina, rok); robustné SE klastrované podľa krajiny
Term	Estimate	Std. Error (CL-HC1)	t value	p-hodnota	Signif
log1p(lag_BEV_PHEV) (t-1)	-2.0426	1.5286	-1.3363	0.1915
HDP	0.0021	0.0038	0.5386	0.5942

Súhrnné metriky lagovaného FE modelu
R²	Adj. R²	AIC	RMSE (resid. se)	n	FE: krajina	FE: rok
0.8533	0.7849	275.1799	4.4185	45	Áno	Áno

## <p style='color:#065F46; font-size: 90%'><em>Pozn.:</em> Štandardné chyby sú clusterované podľa krajiny (HC1). 
## Fixné efekty pre krajinu a rok sú v modeli zahrnuté, ale ich koeficienty sa nezobrazujú.</p>

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

Štandardizovaný OLS: EPI ~ z(log1p(BEV_PHEV)) + z(HDP) (robustné SE, HC1)
Term	Estimate	Std. Error	t value	p value	95% CI low	95% CI high	Stars
Intercept	0.000	0.110	0.000	1.000	-0.216	0.216
log1p(BEV_PHEV) (SD)	-0.792	0.164	-4.838	0.000	-1.113	-0.471	***
HDP (SD)	0.441	0.143	3.079	0.003	0.160	0.722	**
Pozn.:
Koeficienty sú štandardizované (beta). Interpretácia: zmena v EPI v SD jednotkách pri zmene vysvetľujúcej premennej o 1 SD. Hviezdičky: * p<0.001, p<0.01, * p<0.05, · p<0.10.

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.

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

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.

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

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.

4) 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).

4.1) Ramsey RESET – kontrola špecifikácie

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

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

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.

4.3) Kvadratické členy – jemné zakrivenie

Porovnanie kvality modelov (AIC – nižšie je lepšie)
Model	df	AIC
m_lin	4	343.50
m_quad	6	345.03

ANOVA porovnanie modelov (m_quad vs. m_lin)
Model	Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
1	47	2402.115	NA	NA	NA	NA
2	45	2286.639	2	115.4755	1.136253	0.330056

Lineárny model: EPI ~ log1p(BEV_PHEV) + log1p(HDP) (robustné SE, HC1)
Term	Estimate	Std_Error	t_value	p_value	CI_low	CI_high	Stars
(Intercept)	-24.786	37.461	-0.662	0.511	-98.210	48.637
log1p(BEV_PHEV)	-3.815	0.779	-4.896	0.000	-5.342	-2.288	***
log1p(HDP)	12.828	4.159	3.084	0.003	4.676	20.980	**
Pozn.:
Hviezdičky: * p<0.001, p<0.01, * p<0.05, · p<0.10.

Kvadratický model: + [log1p(BEV_PHEV)]² a [log1p(HDP)]² (robustné SE, HC1)
Term	Estimate	Std_Error	t_value	p_value	CI_low	CI_high
(Intercept)	1175.792	1220.852	0.963	0.341	-1217.077	3568.662
log1p(BEV_PHEV)	3.883	6.579	0.590	0.558	-9.011	16.778
I(log1p(BEV_PHEV)^2)	-0.442	0.382	-1.157	0.253	-1.190	0.307
log1p(HDP)	-236.257	249.855	-0.946	0.349	-725.972	253.459
I(log1p(HDP)^2)	12.571	12.626	0.996	0.325	-12.175	37.317
Pozn.:
Hviezdičky: * p<0.001, p<0.01, * p<0.05, · p<0.10.

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.4) Zlom v sklone – iný efekt (P)HEV pri vyšších hodnotách

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

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

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)

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.

5) Zhluková analýza (cluster analysis)

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

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

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

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

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

5.1)Konfigurácia

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

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

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

5.4)Ward.D2 hierarchické zhlukovanie + vizualizácia 5.5)Pomôcky na výber počtu klastrov (k)

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

5.6, 5.7)Priradenie do klastrov a dendrogram s rezom

## cl_best
## 1 2 
## 4 1

5.8)Profilovanie a prezentácia klastrov

Veľkosť klastrov
Klaster	Freq
1	4
2	1

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

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

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

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

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

6) Multikolinearita

Úvod

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

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

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

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

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

Východiskový model a pracovný dataset

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

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

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

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

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

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

Korelačná matica a scatterplotová matica

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

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

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

VIF

Premenna	VIF
Diagnostika multikolinearity (VIF)
Condition number κ ≈ 108.0 \| Pravidlo: VIF > 10 alebo κ > 30 = problém
logBEV	85.18
logBEV_PHEV	87.91
logHDP	1.56
Pozn.: VIF > 10 (prísnejšie > 5) a/alebo κ > 30 signalizuje silnú multikolinearitu.

## [1] 108.0268

VIF a Condition number – interpretácia

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

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

Riešenie: vynechanie premennej

Term	Estimate	Std_Error	t_value	p_value	CI_low	CI_high
Model: EPI ~ log1p(BEV_PHEV) + log1p(HDP)
robustné SE (HC1), 95% CI
(Intercept)	−24.786	37.461	−0.662	0.511	−98.210	48.637
logBEV_PHEV	−3.815	0.779	−4.896	0.000	−5.342	−2.288	***
logHDP	12.828	4.159	3.084	0.003	4.676	20.980	**
Pozn.: * p<0.001, p<0.01, * p<0.05, · p<0.10.

Premenná	VIF
VIF – model s (P)HEV
logBEV_PHEV	1.51
logHDP	1.51
Pravidlo: VIF > 10 (prísnejšie > 5) signalizuje problém.

Term	Estimate	Std_Error	t_value	p_value	CI_low	CI_high
Model: EPI ~ log1p(BEV) + log1p(HDP)
robustné SE (HC1), 95% CI
(Intercept)	−16.736	37.345	−0.448	0.656	−89.933	56.461
logBEV	−3.628	0.762	−4.764	0.000	−5.120	−2.135	***
logHDP	11.659	4.112	2.835	0.007	3.599	19.719	**
Pozn.: * p<0.001, p<0.01, * p<0.05, · p<0.10.

Premenná	VIF
VIF – model s BEV
logBEV	1.47
logHDP	1.47
Pravidlo: VIF > 10 (prísnejšie > 5) signalizuje problém.

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

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

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

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

PCA analýza

Komponenta	Vlastné_číslo	Podiel_variability	Kumulatívne
PCA (logBEV, logBEV_PHEV)
Vysvetlená variabilita – PC1 a PC2
PC1	1.994	99.7%	99.7%
PC2	0.006	0.3%	100.0%

Term	Estimate	Std_Error	t_value	p_value	CI_low	CI_high
Model: EPI ~ PC1_EV + log1p(HDP)
robustné SE (HC1), 95% CI
(Intercept)	−52.852	40.764	−1.297	0.201	−132.749	27.046
PC1_EV	5.034	1.043	4.827	0.000	2.990	7.078	***
logHDP	12.288	4.138	2.969	0.005	4.177	20.399	**
Pozn.: * p<0.001, p<0.01, * p<0.05, · p<0.10.

Premenná	VIF
VIF – model s PC1_EV
PC1_EV	1.49
logHDP	1.49
Pravidlo: VIF > 10 (prísnejšie > 5) signalizuje problém.

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

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

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

Záver práce

Vysvetlenie hypotéz

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

Analýza dát o elektromobilite

Simona Vančová

Úvod

Charakteristika dát

Premenné

Výskumná otázka a hypotézy

Import údajov z mojej bakalárskej práce

Interaktívna tabuľka s dátami

1) Deskriptívna štatistika

Prieskumné grafy - BOXPLOT

2) Lineárna regresia (dve alternatívy)

Diagnostika (krajšie grafy)

Residuals vs Fitted

Normal Q–Q

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

Residuals vs Leverage

Vplyvné pozorovania

Dodatočné porovnanie modelov

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

Štandardizované koeficienty

3) Heteroskedasticita – testy a grafy

Vizuálna diagnostika

Heteroskedasticita - prehľad testov

Formálne testy

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

4.1) Ramsey RESET – kontrola špecifikácie

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

4.3) Kvadratické členy – jemné zakrivenie

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

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

Porovnávacia tabuľka modelov

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

5) Zhluková analýza (cluster analysis)

6) Multikolinearita

Úvod

Východiskový model a pracovný dataset

Korelačná matica a scatterplotová matica

VIF

VIF a Condition number – interpretácia

Riešenie: vynechanie premennej

PCA analýza

Záver práce

Vysvetlenie hypotéz