Econometrics in R - cvičenie 5

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

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: EPI ~ BEV + BEV_PHEV + HDP
model <- lm(EPI ~ 1 + BEV + BEV_PHEV + HDP, data = udaje)

# (3) Súhrn výsledkov
summary(model)

Call:
lm(formula = EPI ~ 1 + BEV + BEV_PHEV + HDP, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-20.337  -4.515  -1.171   5.483  17.717 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  6.378e+01  4.389e+00  14.531   <2e-16 ***
BEV          4.213e-05  4.561e-04   0.092    0.927    
BEV_PHEV    -9.857e-05  2.553e-04  -0.386    0.701    
HDP          3.166e-04  2.304e-04   1.374    0.176    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.771 on 46 degrees of freedom
Multiple R-squared:  0.1432,    Adjusted R-squared:  0.0873 
F-statistic: 2.562 on 3 and 46 DF,  p-value: 0.06626

# (4) Diagnostické grafy 
par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))

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.

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

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.

Scare-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 trochu mení, ale nie dramaticky.

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)

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)  6.3776e+01  4.5684e+00 13.9603   <2e-16 ***
BEV          4.2129e-05  4.3809e-04  0.0962   0.9238    
BEV_PHEV    -9.8571e-05  2.4122e-04 -0.4086   0.6847    
HDP          3.1662e-04  2.1474e-04  1.4744   0.1472    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Pri log-špecifikácii vychádza HDP pozitívne a štatisticky významne, čo je v súlade s očakávaním (vyššia ekonomická úroveň súvisí s lepším EPI). Premenná log1p(BEV) nevyšla štatisticky významne, takže hlavnú hypotézu pre BEV nezamietam ani nepotvrdzujem (efekt sa nepreukázal).

# 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)",       "NEZAMIETAM (nevýznamné)",      "NEPODPORUJEM",                   "V modeli EPI ~ log1p(BEV)+HDP je BEV nevýznamné.",
"BEV_PHEV",  "β_PHEV > 0 (pozitívny)",       "ZAMIETAM (významné)",          "ZAMIETAM (efekt je záporný)",    "log1p(BEV_PHEV) významne negatívne; 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)	NEZAMIETAM (nevýznamné)	NEPODPORUJEM	V modeli EPI ~ log1p(BEV)+HDP je BEV nevýznamné.
BEV_PHEV	β_PHEV > 0 (pozitívny)	ZAMIETAM (významné)	ZAMIETAM (efekt je záporný)	log1p(BEV_PHEV) významne negatívne; 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 (po log-transformácii) jeho efekt nie je štatisticky významný – nulovú hypotézu o nulovom efekte nevieme zamietnuť. V alternatívnom modeli s BEV_PHEV vychádza koeficient významne negatívny. Porovnanie AIC ukazuje, že model s BEV_PHEV a HDP lepšie sedí dátam ako verzia s BEV.

Štandardizované koeficienty naznačujú, že 1 SD nárast log(BEV_PHEV) súvisí s ≈ 0,79 SD poklesom EPI, zatiaľ čo 1 SD nárast HDP súvisí s ≈ 0,44 SD nárastom EPI. Model vysvetľuje približne 40–42 % variability EPI.

Diagnostika rezíduí poukazuje na mierne odchýlky v chvostoch a slabú heteroskedasticitu; preto uvádzam robustné štandardné chyby. Zistené vplyvné pozorovania nemenia kvalitatívne závery.

Celkovo výsledky podporujú význam ekonomickej úrovne krajiny pre environmentálnu výkonnosť. Priamy pozitívny vplyv elektromobility na EPI sa v tomto súbore dát nepreukázal – pri BEV efekt nevieme preukázať a pri BEV_PHEV vychádza negatívna asociácia. Tento nesúlad s pôvodným očakávaním môže súvisieť s načasovaním (efekty sa prejavia neskôr), s reálnymi emisiami PHEV alebo s tým, že EPI zachytáva aj oblasti mimo dopravy.