Ekonometria – Cvičenie 5

Lineárna regresia a diagnostika modelu – udržateľná energetika

1 Úvod

V tomto cvičení nadväzujem na predchádzajúcu analýzu udržateľnej energetiky. Používam databázu Global Data on Sustainable Energy, ktorá obsahuje ekonomické a energetické ukazovatele krajín sveta.

Cieľom je vytvoriť lineárny regresný model, interpretovať jeho výsledky a vykonať základnú diagnostiku.

Výskumná otázka: Ako súvisia ekonomická úroveň krajiny a spotreba energie s podielom nízkouhlíkovej elektriny?

Analýza sa zameriava na rok 2019. Ide o prierezovú analýzu, kde jedno pozorovanie predstavuje jednu krajinu.

2 1. Príprava údajov

2.1 Načítanie databázy

Používam základné funkcie programu R bez potreby inštalácie ďalších balíkov.

subor <- "global-data-on-sustainable-energy (1).csv"

moznosti <- c(
  subor,
  file.path("udaje", subor),
  file.path("data", subor),
  "global-data-on-sustainable-energy.csv"
)

existujuce <- moznosti[file.exists(moznosti)]

if (length(existujuce) == 0) {
  stop(
    "CSV súbor sa nenašiel. Ulož ho do rovnakého priečinka ako Cvicenie5.qmd."
  )
}

udaje <- read.csv(
  existujuce[1],
  check.names = FALSE,
  stringsAsFactors = FALSE
)

cat("Počet riadkov:", nrow(udaje), "\n")

Počet riadkov: 3649

cat("Počet premenných:", ncol(udaje), "\n")

Počet premenných: 21

2.2 Výber premenných

Vysvetľovanou premennou je podiel nízkouhlíkovej elektriny. Medzi vysvetľujúce premenné patria HDP na obyvateľa, spotreba primárnej energie na obyvateľa a energetická intenzita.

d <- udaje[udaje$Year == 2019, ]

reg_data <- data.frame(
  krajina = d[["Entity"]],
  low_carbon = as.numeric(
    d[["Low-carbon electricity (% electricity)"]]
  ),
  gdp_pc = as.numeric(
    d[["gdp_per_capita"]]
  ),
  energy_pc = as.numeric(
    d[["Primary energy consumption per capita (kWh/person)"]]
  ),
  intensity = as.numeric(
    d[["Energy intensity level of primary energy (MJ/$2017 PPP GDP)"]]
  )
)

pocet_pred <- nrow(reg_data)

reg_data <- reg_data[
  complete.cases(reg_data) &
    reg_data$gdp_pc > 0 &
    reg_data$energy_pc > 0 &
    reg_data$intensity > 0 &
    reg_data$low_carbon >= 0 &
    reg_data$low_carbon <= 100,
]

reg_data$gdp_tisic <- reg_data$gdp_pc / 1000
reg_data$energy_tisic <- reg_data$energy_pc / 1000
reg_data$log_gdp <- log(reg_data$gdp_pc)
reg_data$log_energy <- log(reg_data$energy_pc)

rownames(reg_data) <- make.unique(reg_data$krajina)

cat("Počet krajín pred čistením:", pocet_pred, "\n")

Počet krajín pred čistením: 175

cat("Počet krajín po čistení:", nrow(reg_data), "\n")

Počet krajín po čistení: 158

cat("Vyradené pozorovania:", pocet_pred - nrow(reg_data), "\n")

Vyradené pozorovania: 17

knitr::kable(
  head(reg_data, 10),
  digits = 2,
  row.names = FALSE,
  caption = "Ukážka pripravených údajov"
)

Ukážka pripravených údajov
krajina	low_carbon	gdp_pc	energy_pc	intensity	gdp_tisic	energy_tisic	log_gdp	log_energy
Afghanistan	83.18	494.18	802.61	2.41	0.49	0.80	6.20	6.69
Albania	100.00	5395.66	11585.56	2.50	5.40	11.59	8.59	9.36
Algeria	1.10	3989.67	16259.05	5.32	3.99	16.26	8.29	9.70
Angola	71.56	2809.63	3667.42	2.94	2.81	3.67	7.94	8.21
Antigua and Barbuda	2.94	17376.65	33027.52	3.45	17.38	33.03	9.76	10.41
Argentina	31.99	10056.64	20725.17	3.35	10.06	20.73	9.22	9.94
Armenia	60.82	4604.65	16179.10	3.53	4.60	16.18	8.43	9.69
Australia	21.38	54875.29	66263.54	4.30	54.88	66.26	10.91	11.10
Austria	77.11	50114.40	48314.26	2.82	50.11	48.31	10.82	10.79
Azerbaijan	6.87	4805.75	17546.78	4.61	4.81	17.55	8.48	9.77

Chýbajúce hodnoty boli z pracovnej vzorky odstránené. Nenahrádzam ich nulami, pretože by to mohlo skresliť regresnú analýzu.

2.3 Opisná štatistika

opis <- data.frame(
  Ukazovatel = c(
    "Počet krajín",
    "Priemerný podiel nízkouhlíkovej elektriny (%)",
    "Medián podielu nízkouhlíkovej elektriny (%)",
    "Priemerné HDP na obyvateľa (USD)",
    "Priemerná spotreba energie (kWh/os.)"
  ),
  Hodnota = c(
    nrow(reg_data),
    mean(reg_data$low_carbon),
    median(reg_data$low_carbon),
    mean(reg_data$gdp_pc),
    mean(reg_data$energy_pc)
  )
)

knitr::kable(
  opis,
  digits = 2,
  caption = "Základné opisné štatistiky"
)

Základné opisné štatistiky
Ukazovatel	Hodnota
Počet krajín	158.00
Priemerný podiel nízkouhlíkovej elektriny (%)	40.67
Medián podielu nízkouhlíkovej elektriny (%)	35.82
Priemerné HDP na obyvateľa (USD)	15496.82
Priemerná spotreba energie (kWh/os.)	26329.69

3 2. Ukladanie a obnovenie workspace

Workspace predstavuje pracovné prostredie programu R, v ktorom sú uložené všetky vytvorené objekty.

Funkcia save.image() uloží pracovné prostredie do súboru. Funkcia load() ho dokáže obnoviť.

save.image("pamat_cvicenie5.RData")

pomocny_objekt <- 123

rm(pomocny_objekt)

load("pamat_cvicenie5.RData")

cat(
  "Údaje boli úspešne obnovené:",
  exists("reg_data"),
  "\n"
)

Údaje boli úspešne obnovené: TRUE

Na odstránenie všetkých objektov možno použiť rm(list = ls()). V dokumente tento príkaz nespúšťam, pretože by odstránil objekty potrebné na ďalšie výpočty.

4 3. Prvý regresný model

Odhadujem model:

\[ LowCarbon_i = \beta_0 + \beta_1 GDPpc_i + \beta_2 EnergyPC_i + u_i \]

HDP na obyvateľa a spotrebu energie používam v tisícoch pôvodných jednotiek.

model1 <- lm(
  low_carbon ~ gdp_tisic + energy_tisic,
  data = reg_data
)

summary(model1)


Call:
lm(formula = low_carbon ~ gdp_tisic + energy_tisic, data = reg_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-68.444 -26.899  -5.776  22.903  78.644 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   41.0296     3.3132  12.384  < 2e-16 ***
gdp_tisic      0.5023     0.1676   2.997  0.00317 ** 
energy_tisic  -0.3094     0.1018  -3.039  0.00279 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 32.32 on 155 degrees of freedom
Multiple R-squared:  0.06516,   Adjusted R-squared:  0.0531 
F-statistic: 5.402 on 2 and 155 DF,  p-value: 0.005398

4.1 Tabuľka regresných koeficientov

koef1 <- as.data.frame(
  summary(model1)$coefficients
)

koef1$Premenna <- rownames(koef1)

koef1 <- koef1[
  ,
  c(5, 1, 2, 3, 4)
]

knitr::kable(
  koef1,
  digits = 4,
  row.names = FALSE,
  col.names = c(
    "Premenná",
    "Odhad",
    "Štandardná chyba",
    "t-štatistika",
    "p-hodnota"
  ),
  caption = "Koeficienty prvého regresného modelu"
)

Koeficienty prvého regresného modelu
Premenná	Odhad	Štandardná chyba	t-štatistika	p-hodnota
(Intercept)	41.0296	3.3132	12.3835	0.0000
gdp_tisic	0.5023	0.1676	2.9974	0.0032
energy_tisic	-0.3094	0.1018	-3.0389	0.0028

5 4. t-testy regresných koeficientov

t-testom overujem, či sú jednotlivé regresné koeficienty štatisticky významné.

Hypotézy:

\[ H_0: \beta_j = 0 \]

\[ H_1: \beta_j \neq 0 \]

Testovacia štatistika:

\[ t_j = \frac{\hat{\beta}_j}{SE(\hat{\beta}_j)} \]

Používam hladinu významnosti 5 %.

Ak je p-hodnota menšia ako 0,05, nulovú hypotézu zamietam.

tab <- summary(model1)$coefficients

for (i in 2:nrow(tab)) {

  nazov <- rownames(tab)[i]
  beta <- tab[i, 1]
  p <- tab[i, 4]

  cat("\n\n**Premenná:", nazov, "**\n\n")

  cat(
    "Odhad koeficientu:",
    round(beta, 4),
    ". P-hodnota:",
    signif(p, 4),
    ". "
  )

  if (p < 0.05) {
    cat(
      "Nulovú hypotézu zamietam. ",
      "Premenná je štatisticky významná."
    )
  } else {
    cat(
      "Nulovú hypotézu nezamietam. ",
      "Štatistická významnosť sa nepreukázala."
    )
  }

  cat(
    "\n\nZvýšenie tejto premennej o 1 000 pôvodných ",
    "jednotiek je spojené so zmenou podielu ",
    "nízkouhlíkovej elektriny o",
    round(beta, 4),
    "percentuálneho bodu pri nezmenenej druhej premennej."
  )
}

Premenná: gdp_tisic

Odhad koeficientu: 0.5023 . P-hodnota: 0.003173 . Nulovú hypotézu zamietam. Premenná je štatisticky významná.

Zvýšenie tejto premennej o 1 000 pôvodných jednotiek je spojené so zmenou podielu nízkouhlíkovej elektriny o 0.5023 percentuálneho bodu pri nezmenenej druhej premennej.

Premenná: energy_tisic

Odhad koeficientu: -0.3094 . P-hodnota: 0.002788 . Nulovú hypotézu zamietam. Premenná je štatisticky významná.

Zvýšenie tejto premennej o 1 000 pôvodných jednotiek je spojené so zmenou podielu nízkouhlíkovej elektriny o -0.3094 percentuálneho bodu pri nezmenenej druhej premennej.

Výsledky predstavujú štatistické súvislosti, nie dôkaz príčinného vzťahu.

6 5. Koeficient determinácie R²

Koeficient determinácie vyjadruje, akú časť variability vysvetľovanej premennej vysvetľuje regresný model.

\[ R^2 = 1 - \frac{RSS}{TSS} \]

RSS je reziduálny súčet štvorcov a TSS celkový súčet štvorcov.

r2 <- summary(model1)$r.squared

upr_r2 <- summary(model1)$adj.r.squared

cat(
  "Koeficient determinácie R² je **",
  round(r2, 4),
  "**.\n\n",
  sep = ""
)

Koeficient determinácie R² je 0.0652.

cat(
  "Model vysvetľuje približne **",
  round(r2 * 100, 2),
  " %** variability podielu nízkouhlíkovej elektriny.\n\n",
  sep = ""
)

Model vysvetľuje približne 6.52 % variability podielu nízkouhlíkovej elektriny.

cat(
  "Upravený koeficient determinácie je **",
  round(upr_r2, 4),
  "**.",
  sep = ""
)

Upravený koeficient determinácie je 0.0531.

7 6. F-test celkovej významnosti modelu

F-test overuje, či sú vysvetľujúce premenné ako celok štatisticky významné.

Hypotézy:

\[ H_0: \beta_1 = \beta_2 = 0 \]

\[ H_1: \text{aspoň jeden koeficient je nenulový} \]

Testovacia štatistika:

\[ F = \frac{R^2/(k-1)} {(1-R^2)/(n-k)} \]

Kde \(n\) predstavuje počet pozorovaní a \(k\) počet odhadovaných koeficientov vrátane konštanty.

f <- summary(model1)$fstatistic

f_hodnota <- unname(f[1])

p_f <- pf(
  unname(f[1]),
  unname(f[2]),
  unname(f[3]),
  lower.tail = FALSE
)

cat(
  "F-štatistika: **",
  round(f_hodnota, 4),
  "**\n\n",
  sep = ""
)

F-štatistika: 5.4017

cat(
  "P-hodnota: **",
  signif(p_f, 4),
  "**\n\n",
  sep = ""
)

P-hodnota: 0.005398

if (p_f < 0.05) {
  cat(
    "Na hladine 5 % zamietam nulovú hypotézu. ",
    "Regresný model je ako celok štatisticky významný."
  )
} else {
  cat(
    "Na hladine 5 % nulovú hypotézu nezamietam. ",
    "Spoločná štatistická významnosť sa nepreukázala."
  )
}

Na hladine 5 % zamietam nulovú hypotézu. Regresný model je ako celok štatisticky významný.

8 7. Grafická diagnostika regresného modelu

Na diagnostiku používam štyri štandardné grafy funkcie plot(model1).

par(mfrow = c(2, 2))

plot(model1)

par(mfrow = c(1, 1))

8.1 Residuals vs Fitted

Prvý graf zobrazuje rezíduá voči odhadnutým hodnotám.

Pri vhodne špecifikovanom modeli by mali byť rezíduá náhodne rozptýlené okolo nuly. Systematické zakrivenie môže signalizovať nelinearitu a lievikovitý tvar môže naznačovať heteroskedasticitu.

8.2 Normal Q-Q

Druhý graf slúži na vizuálne posúdenie normality rezíduí.

Ak sú rezíduá približne normálne rozdelené, body by mali sledovať referenčnú priamku. Výrazné odchýlky môžu naznačovať problém s normalitou.

8.3 Scale-Location

Tretí graf slúži na posúdenie konštantnosti rozptylu rezíduí.

Rovnomerné rozptýlenie bodov podporuje predpoklad homoskedasticity. Výrazný rastúci alebo klesajúci trend môže signalizovať heteroskedasticitu.

8.4 Residuals vs Leverage

Štvrtý graf pomáha identifikovať pozorovania, ktoré môžu výrazne ovplyvňovať regresný model.

Leverage vyjadruje neobvyklosť pozorovania vzhľadom na vysvetľujúce premenné. Cookova vzdialenosť pomáha hodnotiť jeho vplyv na odhadnutý model.

9 8. Jarqueov–Berov test normality

Jarqueov–Berov test overuje, či majú rezíduá normálne rozdelenie.

Hypotézy:

\[ H_0: \text{rezíduá majú normálne rozdelenie} \]

\[ H_1: \text{rezíduá nemajú normálne rozdelenie} \]

Testovacia štatistika:

\[ JB = \frac{n}{6} \left( S^2 + \frac{(K-3)^2}{4} \right) \]

Kde \(S\) predstavuje šikmosť a \(K\) špicatosť rozdelenia.

Nasledujúca funkcia vypočíta Jarqueov–Berov test pomocou základných funkcií R.

jb_test <- function(model) {

  e <- residuals(model)
  n <- length(e)

  priemer <- mean(e)
  odchylka <- sqrt(mean((e - priemer)^2))

  s <- mean((e - priemer)^3) / odchylka^3
  k <- mean((e - priemer)^4) / odchylka^4

  jb <- n / 6 * (s^2 + (k - 3)^2 / 4)

  p <- pchisq(
    jb,
    df = 2,
    lower.tail = FALSE
  )

  data.frame(
    JB_statistika = jb,
    P_hodnota = p
  )
}

jb1 <- jb_test(model1)

knitr::kable(
  jb1,
  digits = 4,
  caption = "Jarqueov–Berov test – model 1"
)

Jarqueov–Berov test – model 1
JB_statistika	P_hodnota
9.6552	0.008

9.1 Interpretácia normality

p <- jb1$P_hodnota

cat(
  "P-hodnota Jarqueovho–Berovho testu je **",
  signif(p, 4),
  "**.\n\n",
  sep = ""
)

P-hodnota Jarqueovho–Berovho testu je 0.008006.

if (p < 0.05) {
  cat(
    "Na hladine významnosti 5 % zamietam nulovú ",
    "hypotézu normality rezíduí. Výsledok naznačuje ",
    "odchýlku od normálneho rozdelenia."
  )
} else {
  cat(
    "Na hladine významnosti 5 % nulovú hypotézu ",
    "nezamietam. Test neposkytuje dostatočný dôkaz ",
    "proti normalite rezíduí."
  )
}

Na hladine významnosti 5 % zamietam nulovú hypotézu normality rezíduí. Výsledok naznačuje odchýlku od normálneho rozdelenia.

10 9. Identifikácia odľahlých pozorovaní

Na identifikáciu odľahlých pozorovaní používam externe študentizované rezíduá a Bonferroniho korekciu.

Táto metóda predstavuje rovnaký princíp, aký využíva funkcia car::outlierTest().

Hypotézy:

\[ H_0: \text{pozorovanie nie je odľahlé} \]

\[ H_1: \text{pozorovanie je odľahlé} \]

Pri Bonferroniho p-hodnote menšej ako 0,05 považujem pozorovanie za štatisticky významne odľahlé.

outlier_test <- function(model) {

  r <- rstudent(model)

  df <- df.residual(model) - 1

  p <- 2 * pt(
    abs(r),
    df = df,
    lower.tail = FALSE
  )

  p_bonf <- p.adjust(
    p,
    method = "bonferroni"
  )

  vysledok <- data.frame(
    Krajina = names(r),
    Studentizovane_reziduum = as.numeric(r),
    P_hodnota = as.numeric(p),
    Bonferroni_p = as.numeric(p_bonf)
  )

  vysledok <- vysledok[
    order(vysledok$Bonferroni_p),
  ]

  vysledok
}

out1 <- outlier_test(model1)

knitr::kable(
  head(out1, 10),
  digits = 4,
  row.names = FALSE,
  caption = "Pozorovania s najnižšími Bonferroniho p-hodnotami"
)

Pozorovania s najnižšími Bonferroniho p-hodnotami
Krajina	Studentizovane_reziduum	P_hodnota	Bonferroni_p
Afghanistan	1.3140	0.1908	1
Albania	1.8741	0.0628	1
Algeria	-1.1480	0.2527	1
Angola	0.9401	0.3486	1
Antigua and Barbuda	-1.1373	0.2572	1
Argentina	-0.2377	0.8125	1
Armenia	0.6974	0.4866	1
Australia	-0.8380	0.4033	1
Austria	0.8099	0.4193	1
Azerbaijan	-0.9674	0.3348	1

10.1 Interpretácia odľahlých pozorovaní

vyznamne <- out1[
  out1$Bonferroni_p < 0.05,
]

if (nrow(vyznamne) > 0) {

  cat(
    "Štatisticky významne odľahlé krajiny: **",
    paste(vyznamne$Krajina, collapse = ", "),
    "**.",
    sep = ""
  )

} else {

  cat(
    "Na hladine významnosti 5 % nebolo ",
    "identifikované žiadne štatisticky významne ",
    "odľahlé pozorovanie podľa Bonferroniho korekcie."
  )
}

Na hladine významnosti 5 % nebolo identifikované žiadne štatisticky významne odľahlé pozorovanie podľa Bonferroniho korekcie.

Odľahlé pozorovania automaticky neodstraňujem, pretože môžu predstavovať skutočné rozdiely medzi krajinami.

11 10. Leverage a Cookova vzdialenosť

Ako doplňujúcu diagnostiku vypočítam leverage a Cookovu vzdialenosť.

vplyv <- data.frame(
  Krajina = names(residuals(model1)),
  Leverage = as.numeric(hatvalues(model1)),
  Cookova_vzdialenost = as.numeric(
    cooks.distance(model1)
  )
)

vplyv <- vplyv[
  order(-vplyv$Cookova_vzdialenost),
]

knitr::kable(
  head(vplyv, 10),
  digits = 4,
  row.names = FALSE,
  caption = "Desať najvplyvnejších pozorovaní"
)

Desať najvplyvnejších pozorovaní
Krajina	Leverage	Cookova_vzdialenost
Iceland	0.1251	0.3227
Cayman Islands	0.1258	0.2461
Norway	0.0579	0.0494
Ireland	0.1052	0.0414
Canada	0.0416	0.0383
Switzerland	0.1190	0.0255
Sweden	0.0254	0.0217
Singapore	0.0988	0.0213
Israel	0.0233	0.0181
Bhutan	0.0106	0.0149

plot(
  cooks.distance(model1),
  type = "h",
  main = "Cookova vzdialenosť",
  xlab = "Poradie pozorovania",
  ylab = "Cookova vzdialenosť"
)

abline(
  h = 4 / nobs(model1),
  col = "red",
  lty = 2
)

Červená čiara predstavuje orientačnú hranicu \(4/n\). Jej prekročenie samo osebe neznamená, že pozorovanie treba odstrániť.

12 11. Druhý model – logaritmická transformácia

V druhom modeli používam logaritmickú transformáciu vysvetľujúcich premenných.

\[ LowCarbon_i = \beta_0 + \beta_1 \ln(GDPpc_i) + \beta_2 \ln(EnergyPC_i) + u_i \]

Logaritmická transformácia môže pomôcť zachytiť nelineárny vzťah medzi premennými.

model2 <- lm(
  low_carbon ~ log_gdp + log_energy,
  data = reg_data
)

summary(model2)


Call:
lm(formula = low_carbon ~ log_gdp + log_energy, data = reg_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-49.276 -27.861  -6.459  27.091  72.184 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   57.220     16.656   3.435 0.000759 ***
log_gdp        7.603      4.560   1.667 0.097448 .  
log_energy    -8.915      4.329  -2.059 0.041139 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 32.95 on 155 degrees of freedom
Multiple R-squared:  0.02845,   Adjusted R-squared:  0.01592 
F-statistic:  2.27 on 2 and 155 DF,  p-value: 0.1068

12.1 Regresné koeficienty

knitr::kable(
  as.data.frame(summary(model2)$coefficients),
  digits = 4,
  caption = "Regresné koeficienty druhého modelu"
)

Regresné koeficienty druhého modelu
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	57.2203	16.6561	3.4354	0.0008
log_gdp	7.6027	4.5595	1.6674	0.0974
log_energy	-8.9154	4.3293	-2.0593	0.0411

Pri zvýšení vysvetľujúcej premennej približne o 1 % sa očakávaný podiel nízkouhlíkovej elektriny zmení približne o \(\beta_j/100\) percentuálneho bodu pri nezmenenej druhej premennej.

12.2 Diagnostické grafy druhého modelu

par(mfrow = c(2, 2))

plot(model2)

par(mfrow = c(1, 1))

12.3 Test normality druhého modelu

jb2 <- jb_test(model2)

knitr::kable(
  jb2,
  digits = 4,
  caption = "Jarqueov–Berov test – model 2"
)

Jarqueov–Berov test – model 2
JB_statistika	P_hodnota
11.0154	0.0041

12.4 Test odľahlých pozorovaní druhého modelu

out2 <- outlier_test(model2)

knitr::kable(
  head(out2, 10),
  digits = 4,
  row.names = FALSE,
  caption = "Test odľahlých pozorovaní – model 2"
)

Test odľahlých pozorovaní – model 2
Krajina	Studentizovane_reziduum	P_hodnota	Bonferroni_p
Afghanistan	1.1838	0.2383	1
Albania	1.8691	0.0635	1
Algeria	-1.0024	0.3177	1
Angola	0.8271	0.4095	1
Antigua and Barbuda	-1.0907	0.2771	1
Argentina	-0.2030	0.8394	1
Armenia	0.7903	0.4306	1
Australia	-0.6077	0.5443	1
Austria	1.0380	0.3009	1
Azerbaijan	-0.8458	0.3990	1

13 12. Tretí model – rozšírenie o energetickú intenzitu

Ako ďalšiu modifikáciu pridávam energetickú intenzitu hospodárstva.

\[ LowCarbon_i = \beta_0 + \beta_1 \ln(GDPpc_i) + \beta_2 \ln(EnergyPC_i) + \beta_3 Intensity_i + u_i \]

Cieľom je overiť, či energetická intenzita prináša dodatočnú vysvetľovaciu schopnosť.

model3 <- lm(
  low_carbon ~ log_gdp + log_energy + intensity,
  data = reg_data
)

summary(model3)


Call:
lm(formula = low_carbon ~ log_gdp + log_energy + intensity, data = reg_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-49.284 -27.896  -6.477  27.056  72.192 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  57.3762    19.2331   2.983  0.00332 **
log_gdp       7.5625     5.1897   1.457  0.14709   
log_energy   -8.8853     4.7159  -1.884  0.06144 . 
intensity    -0.0184     1.1231  -0.016  0.98695   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 33.05 on 154 degrees of freedom
Multiple R-squared:  0.02845,   Adjusted R-squared:  0.009527 
F-statistic: 1.503 on 3 and 154 DF,  p-value: 0.2159

13.1 Regresné koeficienty

knitr::kable(
  as.data.frame(summary(model3)$coefficients),
  digits = 4,
  caption = "Regresné koeficienty tretieho modelu"
)

Regresné koeficienty tretieho modelu
	Estimate	Std. Error	t value	Pr(>\|t\|)
(Intercept)	57.3762	19.2331	2.9832	0.0033
log_gdp	7.5625	5.1897	1.4572	0.1471
log_energy	-8.8853	4.7159	-1.8841	0.0614
intensity	-0.0184	1.1231	-0.0164	0.9870

13.2 Diagnostické grafy tretieho modelu

par(mfrow = c(2, 2))

plot(model3)

par(mfrow = c(1, 1))

13.3 Test normality tretieho modelu

jb3 <- jb_test(model3)

knitr::kable(
  jb3,
  digits = 4,
  caption = "Jarqueov–Berov test – model 3"
)

Jarqueov–Berov test – model 3
JB_statistika	P_hodnota
11.0223	0.004

13.4 Test odľahlých pozorovaní tretieho modelu

out3 <- outlier_test(model3)

knitr::kable(
  head(out3, 10),
  digits = 4,
  row.names = FALSE,
  caption = "Test odľahlých pozorovaní – model 3"
)

Test odľahlých pozorovaní – model 3
Krajina	Studentizovane_reziduum	P_hodnota	Bonferroni_p
Afghanistan	1.1866	0.2372	1
Albania	1.8687	0.0636	1
Algeria	-1.0002	0.3188	1
Angola	0.8248	0.4108	1
Antigua and Barbuda	-1.0879	0.2783	1
Argentina	-0.2033	0.8392	1
Armenia	0.7890	0.4313	1
Australia	-0.6057	0.5456	1
Austria	1.0349	0.3024	1
Azerbaijan	-0.8448	0.3996	1

14 13. Porovnanie regresných modelov

Porovnávam tri modely pomocou koeficientu determinácie, upraveného koeficientu determinácie, AIC, BIC a testu normality rezíduí.

Akaikeho informačné kritérium (AIC) hodnotí kompromis medzi kvalitou prispôsobenia modelu a jeho zložitosťou.

Pri modeloch odhadnutých na rovnakých pozorovaniach s rovnakou vysvetľovanou premennou uprednostňujem nižšiu hodnotu AIC.

modely <- list(
  model1,
  model2,
  model3
)

nazvy <- c(
  "Model 1 – lineárny",
  "Model 2 – logaritmický",
  "Model 3 – rozšírený"
)

porovnanie <- data.frame(
  Model = nazvy,
  R2 = sapply(
    modely,
    function(m) summary(m)$r.squared
  ),
  Upr_R2 = sapply(
    modely,
    function(m) summary(m)$adj.r.squared
  ),
  AIC = sapply(modely, AIC),
  BIC = sapply(modely, BIC),
  JB_p = sapply(
    modely,
    function(m) jb_test(m)$P_hodnota
  ),
  N = sapply(modely, nobs)
)

knitr::kable(
  porovnanie,
  digits = 4,
  row.names = FALSE,
  caption = "Porovnanie troch regresných modelov"
)

Porovnanie troch regresných modelov
Model	R2	Upr_R2	AIC	BIC	JB_p	N
Model 1 – lineárny	0.0652	0.0531	1551.646	1563.896	0.0080	158
Model 2 – logaritmický	0.0285	0.0159	1557.730	1569.981	0.0041	158
Model 3 – rozšírený	0.0285	0.0095	1559.730	1575.043	0.0040	158

14.1 Výber modelu podľa AIC

naj <- which.min(porovnanie$AIC)

cat(
  "**Najnižšiu hodnotu AIC dosiahol ",
  porovnanie$Model[naj],
  " s hodnotou ",
  round(porovnanie$AIC[naj], 4),
  ".**\n\n",
  sep = ""
)

Najnižšiu hodnotu AIC dosiahol Model 1 – lineárny s hodnotou 1551.646.

cat(
  "Podľa informačného kritéria AIC je tento ",
  "model spomedzi porovnávaných špecifikácií ",
  "preferovaný. Pri konečnom hodnotení je však ",
  "potrebné zohľadniť aj diagnostiku rezíduí."
)

Podľa informačného kritéria AIC je tento model spomedzi porovnávaných špecifikácií preferovaný. Pri konečnom hodnotení je však potrebné zohľadniť aj diagnostiku rezíduí.

14.2 Porovnanie normality rezíduí

for (i in seq_len(nrow(porovnanie))) {

  cat(
    "\n\n**",
    porovnanie$Model[i],
    ":** ",
    sep = ""
  )

  if (porovnanie$JB_p[i] < 0.05) {
    cat(
      "Jarqueov–Berov test zamieta normalitu ",
      "rezíduí na hladine 5 %."
    )
  } else {
    cat(
      "Jarqueov–Berov test nezamieta normalitu ",
      "rezíduí na hladine 5 %."
    )
  }
}

Model 1 – lineárny: Jarqueov–Berov test zamieta normalitu rezíduí na hladine 5 %.

Model 2 – logaritmický: Jarqueov–Berov test zamieta normalitu rezíduí na hladine 5 %.

Model 3 – rozšírený: Jarqueov–Berov test zamieta normalitu rezíduí na hladine 5 %.

15 14. Záverečné zhodnotenie

V tomto cvičení som analyzoval súvislosť medzi ekonomickou úrovňou, spotrebou energie a podielom nízkouhlíkovej elektriny v krajinách sveta.

Pomocou t-testov som overil štatistickú významnosť jednotlivých regresných koeficientov.

F-testom som posúdil celkovú významnosť modelu a pomocou koeficientu determinácie jeho vysvetľovaciu schopnosť.

Následne som vykonal diagnostiku rezíduí pomocou štyroch štandardných grafov, Jarqueovho–Berovho testu a identifikácie odľahlých pozorovaní.

Vytvoril som tri regresné modely a porovnal ich pomocou informačného kritéria AIC.

cat(
  "Najnižšiu hodnotu AIC dosiahol **",
  porovnanie$Model[naj],
  "**.\n\n",
  sep = ""
)

Najnižšiu hodnotu AIC dosiahol Model 1 – lineárny.

cat(
  "Výsledky predstavujú štatistické súvislosti, ",
  "nie dôkaz príčinných vzťahov. ",
  "Energetická štruktúra krajín závisí aj od ",
  "energetickej politiky, prírodných podmienok ",
  "a historického vývoja, ktoré regresný model ",
  "nemusí dostatočne zachytávať."
)

Výsledky predstavujú štatistické súvislosti, nie dôkaz príčinných vzťahov. Energetická štruktúra krajín závisí aj od energetickej politiky, prírodných podmienok a historického vývoja, ktoré regresný model nemusí dostatočne zachytávať.

16 15. Uloženie výsledkov

Na konci cvičenia ukladám celé pracovné prostredie vrátane pripravených údajov a regresných modelov.

save.image("pamat_cvicenie5.RData")

cat("Pracovné prostredie bolo úspešne uložené.")

Pracovné prostredie bolo úspešne uložené.

17 Použité zdroje

Gazda, V. Econometrics in R – cvičenie 5: Základná diagnostika modelu. https://rpubs.com/gazdavladimir/1472929
Gazda, V. Hat matica, leverage a diagnostika vplyvných pozorovaní. https://rpubs.com/gazdavladimir/1474655
Global Data on Sustainable Energy. Kaggle. https://www.kaggle.com/datasets/anshtanwar/global-data-on-sustainable-energy
R Core Team. R: A Language and Environment for Statistical Computing.