Cvičenie 6 - Heteroskedasticita

---

library(zoo)

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(tseries)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

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

## Loading required package: carData

rm(list=ls())

# Načítanie dát
udaje <- read.csv2("Chinas GDP in Province En.csv", header = TRUE, sep = ",", dec = ".")  

# Pozrime si štruktúru
str(udaje)

## 'data.frame':    29 obs. of  32 variables:
##  $ X             : int  2020 2019 2018 2017 2016 2015 2014 2013 2012 2011 ...
##  $ Beijing       : num  36103 35445 33106 29883 27041 ...
##  $ Tianjin       : num  14084 14056 13363 12451 11477 ...
##  $ Hebei         : num  36207 34979 32495 30641 28474 ...
##  $ Shanxi        : num  17652 16962 15958 14484 11946 ...
##  $ Inner.Mongolia: num  17360 17212 16141 14898 13789 ...
##  $ Liaoning      : num  25115 24855 23510 21693 20392 ...
##  $ Jilin         : num  12311 11727 11254 10922 10427 ...
##  $ Heilongjiang  : num  13698 13544 12846 12313 11895 ...
##  $ Shanghai      : num  38701 37988 36012 32925 29887 ...
##  $ Jiangsu       : num  102719 98657 93208 85870 77351 ...
##  $ Zhejiang      : num  64613 62462 58003 52403 47254 ...
##  $ Anhui         : num  38681 36846 34011 29676 26308 ...
##  $ Fujian        : num  43904 42327 38688 33842 29609 ...
##  $ Jiangxi       : num  25692 24667 22716 20211 18389 ...
##  $ Shandong      : num  73129 70540 66649 63012 58762 ...
##  $ Henan         : num  54997 53718 49936 44825 40249 ...
##  $ Hubei         : num  43444 45429 42022 37235 33353 ...
##  $ Hunan         : num  41782 39894 36330 33828 30854 ...
##  $ Guangdong     : num  110761 107987 99945 91649 82163 ...
##  $ Guangxi.      : num  22157 21237 19628 17791 16117 ...
##  $ Hainan        : num  5532 5331 4911 4498 4090 ...
##  $ Chongqing     : num  25003 23606 21589 20066 18023 ...
##  $ Sichuan       : num  48599 46364 42902 37905 33138 ...
##  $ Guizhou       : num  17827 16769 15353 13605 11792 ...
##  $ Yunnan        : num  24522 23224 20881 18486 16369 ...
##  $ Tibet         : num  1903 1698 1548 1349 1173 ...
##  $ Shaanxi       : num  26182 25793 23942 21474 19046 ...
##  $ Gansu         : num  9017 8718 8104 7337 6908 ...
##  $ Qinghai       : num  3006 2941 2748 2465 2258 ...
##  $ Ningxia       : num  3921 3748 3510 3200 2781 ...
##  $ Xinjiang      : num  13798 13597 12809 11160 9631 ...

head(udaje)

# Premenujeme stĺpec s rokmi, ak má neštandardný názov
names(udaje)[1] <- "Year"

Úvod do problému, stanovenie hypotéz

Rozhodla som sa modelovať vývoj hrubého domáceho produktu (HDP) v jednotlivých čínskych provinciách v závislosti od času (roku). Cieľom je preskúmať, ako sa HDP mení v priebehu rokov a či možno tento vývoj opísať jednoduchým lineárnym trendom. Analýza umožní porovnať tempo hospodárskeho rastu medzi rôznymi regiónmi a identifikovať provincie, ktoré sa vyznačujú nadpriemerným rastom (napr. Guangdong, Shanghai) alebo naopak pomalším vývojom (napr. Tibet, Gansu).

Moja pracovná hypotéza predpokladá, že: - medzi rokmi a hodnotou HDP existuje kladná štatisticky významná závislosť, teda s rastúcim rokom (časom) HDP rastie, - tempo rastu HDP sa môže líšiť medzi provinciami, čo odráža regionálne rozdiely v rozvoji, investíciách a industrializácii, - niektoré provincie (napr. východné pobrežie) vykazujú výrazne rýchlejší rast ako vnútrozemské oblasti.

Príprava databázy, čistenie a úprava údajov

# Načítanie dát o HDP čínskych provincií

udaje <- read.csv("Chinas GDP in Province En.csv", header = TRUE, sep = ",", dec = ".")
names(udaje)[1] <- "Year"

# Zameriame sa len na obdobie po roku 2000
udaje.modern <- udaje[udaje$Year >= 2000, ]

# Čistenie údajov – doplnenie chýbajúcich hodnôt mediánom

# Zistíme mediány pre každý stĺpec (okrem Year)
column_medians <- sapply(udaje.modern[, -1], median, na.rm = TRUE)

# Imputujeme chýbajúce hodnoty
udaje_imputed <- udaje.modern
for (col in names(udaje.modern)[-1]) {
  udaje_imputed[[col]][is.na(udaje.modern[[col]])] <- column_medians[col]
}

# Uložíme čisté dáta
udaje.modern <- udaje_imputed

# Kontrola kvality dát po čistení

# Počet chýbajúcich hodnôt po imputácii
missing_counts <- colSums(is.na(udaje.modern))
print("Počet chýbajúcich hodnôt po čistení:")

## [1] "Počet chýbajúcich hodnôt po čistení:"

print(missing_counts)

##           Year        Beijing        Tianjin          Hebei         Shanxi 
##              0              0              0              0              0 
## Inner.Mongolia       Liaoning          Jilin   Heilongjiang       Shanghai 
##              0              0              0              0              0 
##        Jiangsu       Zhejiang          Anhui         Fujian        Jiangxi 
##              0              0              0              0              0 
##       Shandong          Henan          Hubei          Hunan      Guangdong 
##              0              0              0              0              0 
##       Guangxi.         Hainan      Chongqing        Sichuan        Guizhou 
##              0              0              0              0              0 
##         Yunnan          Tibet        Shaanxi          Gansu        Qinghai 
##              0              0              0              0              0 
##        Ningxia       Xinjiang 
##              0              0

# Počet nulových hodnôt (potenciálne problematické)
zero_counts <- sapply(udaje.modern[, -1], function(x) sum(x == 0, na.rm = TRUE))
print("Počet nulových hodnôt v jednotlivých provinciách:")

## [1] "Počet nulových hodnôt v jednotlivých provinciách:"

print(zero_counts)

##        Beijing        Tianjin          Hebei         Shanxi Inner.Mongolia 
##              0              0              0              0              0 
##       Liaoning          Jilin   Heilongjiang       Shanghai        Jiangsu 
##              0              0              0              0              0 
##       Zhejiang          Anhui         Fujian        Jiangxi       Shandong 
##              0              0              0              0              0 
##          Henan          Hubei          Hunan      Guangdong       Guangxi. 
##              0              0              0              0              0 
##         Hainan      Chongqing        Sichuan        Guizhou         Yunnan 
##              0              0              0              0              0 
##          Tibet        Shaanxi          Gansu        Qinghai        Ningxia 
##              0              0              0              0              0 
##       Xinjiang 
##              0

# Náhľad na vyčistené dáta
head(udaje.modern)

# Vizualizácia tvaru údajov – boxploty pre jednotlivé provincie

# Vynecháme stĺpec s rokom
udaje_numeric <- udaje.modern[, -1]

# Počet premenných (provincie)
num_plots <- ncol(udaje_numeric)

# Zistíme globálny rozsah hodnôt, aby všetky boxploty mali rovnakú os Y
global_range <- range(udaje_numeric, na.rm = TRUE)

# Nastavíme rozloženie grafov (napr. 4x4, aby sa zmestilo viac naraz)
par(mfrow = c(4, 4))
par(mar = c(3.5, 3.5, 2, 1))  # okraje grafov

# Kreslenie boxplotov pre každú provinciu
for (col in names(udaje_numeric)) {
  data_col <- udaje_numeric[[col]]
  
  # Ak sa vyskytujú nulové hodnoty, zvýrazni ich farbou
  boxplot(data_col,
          main = col,
          col = ifelse(any(data_col == 0, na.rm = TRUE), "tomato", "lightblue"),
          ylab = "HDP (mld. CNY)",
          ylim = global_range,
          cex.main = 0.8)
}

# Pridáme globálny titulok
mtext("Boxploty HDP jednotlivých provincií Číny (2000–2020)", 
      outer = TRUE, cex = 1.3, font = 2)

# Reset rozloženia grafov
par(mfrow = c(1, 1))

# Kontrola nulových hodnôt

zero_counts <- sapply(udaje_numeric, function(x) sum(x == 0, na.rm = TRUE))
nulove_df <- data.frame(Provincia = names(zero_counts), PocetNul = zero_counts)
nulove_df <- subset(nulove_df, PocetNul > 0)

if (nrow(nulove_df) > 0) {
  print("Provincie obsahujúce nulové hodnoty HDP:")
  print(nulove_df)
} else {
  print("V dátach sa nenachádzajú žiadne nulové hodnoty HDP.")
}

## [1] "V dátach sa nenachádzajú žiadne nulové hodnoty HDP."

Lineárna regresia

# Vytvorenie podmnožiny dát pre konkrétnu provinciu

province_name <- "Guangdong"   # môžeš zmeniť napr. na "Beijing", "Shanghai", "Tibet", atď.

data_prov <- data.frame(
  Year = udaje$Year,
  GDP = udaje[[province_name]]
)

# Odhad lineárneho modelu rastu HDP v čase

model <- lm(GDP ~ Year, data = data_prov)

# Základné výstupy modelu

cat("Odhadnuté koeficienty modelu pre provinciu", province_name, ":\n")

## Odhadnuté koeficienty modelu pre provinciu Guangdong :

print(model$coefficients)

##  (Intercept)         Year 
## -7762182.625     3888.906

cat("\nVysvetlenie:\n")

## 
## Vysvetlenie:

cat("Intercept predstavuje odhadovanú úroveň HDP na začiatku časového obdobia,\n")

## Intercept predstavuje odhadovanú úroveň HDP na začiatku časového obdobia,

cat("koeficient pri 'Year' predstavuje priemerný ročný nárast HDP v miliardách CNY.\n")

## koeficient pri 'Year' predstavuje priemerný ročný nárast HDP v miliardách CNY.

# Súhrn modelu (R², p-hodnoty, t-testy)
summary(model)

## 
## Call:
## lm(formula = GDP ~ Year, data = data_prov)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -13111  -9805  -2451   9756  18469 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -7762182.6   478821.1  -16.21 1.94e-15 ***
## Year            3888.9      238.7   16.29 1.71e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10750 on 27 degrees of freedom
## Multiple R-squared:  0.9077, Adjusted R-squared:  0.9043 
## F-statistic: 265.4 on 1 and 27 DF,  p-value: 1.712e-15

# Príklady použitia jednotlivých komponentov modelu

cat("Odhadnuté koeficienty:\n")

## Odhadnuté koeficienty:

print(model$coefficients)

##  (Intercept)         Year 
## -7762182.625     3888.906

cat("\nRezíduá modelu:\n")

## 
## Rezíduá modelu:

print(head(model$residuals))

##          1          2          3          4          5          6 
## 17353.6834 18468.5893 14315.7952  9908.2010  4311.6069   769.7128

cat("\nVyrovnané (predikované) hodnoty HDP:\n")

## 
## Vyrovnané (predikované) hodnoty HDP:

print(head(model$fitted.values))

##        1        2        3        4        5        6 
## 93407.22 89518.31 85629.40 81740.50 77851.59 73962.69

cat("\nMatrica X použitá v modeli:\n")

## 
## Matrica X použitá v modeli:

print(head(model.matrix(model)))

##   (Intercept) Year
## 1           1 2020
## 2           1 2019
## 3           1 2018
## 4           1 2017
## 5           1 2016
## 6           1 2015

# Nastavenie rozloženia 2x2
par(mfrow = c(2, 2))

# Vykreslenie diagnostických grafov
plot(model)

# Voliteľne pridať spoločný nadpis
mtext(paste("Diagnostické grafy regresného modelu HDP – provincia", province_name),
      outer = TRUE, cex = 1.2, font = 2)

# Reset rozloženia grafov
par(mfrow = c(1, 1))

Interpretácia

Kladný koeficient pri prem. Year → HDP rastie v čase (čo potvrdzuje ekonomický rast). R² → udáva, akú časť variability HDP vysvetľuje samotný čas (rok). Rezíduá kolíšu okolo nuly → model nemá systematickú chybu. Q–Q graf → ak body ležia blízko priamky, predpoklad normality rezíduí je splnený.

Bonus - vizualizácia modelu

# Vizualizácia HDP a regresnej priamky
plot(data_prov$Year, data_prov$GDP,
     main = paste("Vývoj HDP provincie", province_name),
     xlab = "Rok", ylab = "HDP (mld. CNY)",
     pch = 16, col = "skyblue")

# Regresná priamka
abline(model, col = "red", lwd = 2)

# Pridanie hodnoty R² do grafu
r2 <- summary(model)$r.squared
legend("topleft", legend = paste("R² =", round(r2, 3)),
       bty = "n", text.col = "red")

Residuals vs. fitted

# Residuals vs. Fitted pre model rastu HDP provincie

# Predpokladáme, že máme vytvorený model:
# model <- lm(GDP ~ Year, data = data_prov)
# kde data_prov obsahuje stĺpce "Year" a "GDP" pre jednu provinciu

# Vykreslenie grafu
plot(model$fitted.values, model$residuals,
     main = paste("Rezíduá vs. predikované hodnoty HDP –", province_name),
     xlab = "Predikované hodnoty HDP (mld. CNY)",
     ylab = "Rezíduá (skutočné – predikované hodnoty)",
     pch = 19, col = "steelblue")

# Pridáme referenčnú čiaru y = 0
abline(h = 0, col = "red", lwd = 2)

# Pridáme hladkú LOESS krivku
lines(lowess(model$fitted.values, model$residuals), col = "orange", lwd = 2)

Interpretácia grafu

Centrovanie okolo nuly: Rezíduá kolíšu približne okolo osi 0 – to je pozitívne znamenie. Model teda nepreceňuje ani nepodceňuje HDP systematicky.

Tvar červenej čiary (LOESS): Ak je čiara približne rovná, vzťah medzi rokom a HDP je dobre vystihnutý lineárne. Mierne zakrivenie by mohlo naznačovať, že rast HDP sa v posledných rokoch spomaľuje alebo zrýchľuje (teda možnú nelinearitu).

Rozptyl rezíduí: Vertikálny rozptyl bodov je približne konštantný → homoskedasticita. Ak by sa body rozširovali do tvaru lievika, išlo by o heteroskedasticitu (nestálu variabilitu chýb).

Odľahlé hodnoty: Ak sa niektoré body výrazne odlišujú od ostatných (napr. prudké skoky HDP počas kríz či pandémie), môžu byť vplyvné pozorovania. Ich vplyv môžeme preskúmať pomocou testu z balíka car:

# Kontrola odľahlých hodnôt v modeli HDP

library(car)

outlier_results <- outlierTest(model)
if (!is.null(outlier_results)) {
  print(outlier_results)
} else {
  cat("Neboli zistené žiadne štatisticky významné odľahlé roky v HDP.\n")
}

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##   rstudent unadjusted p-value Bonferroni p
## 2 1.916687           0.066331           NA

Q-Q plot

# Q–Q GRAF REZÍDUÍ PRE MODEL HDP

# Vykreslenie Q–Q grafu
qqnorm(model$residuals,
       main = paste("Q–Q graf rezíduí –", province_name),
       pch = 19, col = "steelblue")
qqline(model$residuals, col = "red", lwd = 2)

Interpretácia

Celkový tvar Väčšina bodov sa nachádza blízko červenej referenčnej čiary, čo naznačuje, že rezíduá majú približne normálne rozdelenie. To je priaznivý výsledok – lineárny model HDP je štatisticky dobre špecifikovaný.

Krajné hodnoty (extrémy) Na oboch koncoch grafu (vľavo dole a vpravo hore) sa môže objaviť mierne odchýlenie od priamky. To poukazuje na niekoľko extrémnych rokov – napríklad prudký rast alebo pokles HDP počas krízových období (napr. 2008 – globálna finančná kríza, alebo 2020 – pandémia COVID-19). Tieto odchýlky sú však ekonomicky interpretovateľné a nie chybou v dátach.

Stredná oblasť V centrálnej časti grafu (štandardizované kvantily medzi −1 a +1) body presne sledujú priamku, čo potvrdzuje, že väčšina rezíduí modelu je normálne rozdelená.

Overenie normality rezíduí pomocou testu

# TEST NORMALITY REZÍDUÍ – SHAPIRO-WILK A JARQUE-BERA

library(tseries)

# Shapiro–Wilk test (vhodný pre menší počet pozorovaní)
shapiro_result <- shapiro.test(model$residuals)

# Jarque–Bera test (vhodný pre ekonomické časové rady)
jb_test <- jarque.bera.test(model$residuals)

cat("\n Výsledky testov normality:\n")

## 
##  Výsledky testov normality:

cat("Shapiro–Wilk p-hodnota:", round(shapiro_result$p.value, 4), "\n")

## Shapiro–Wilk p-hodnota: 0.0155

cat("Jarque–Bera p-hodnota:", round(jb_test$p.value, 4), "\n")

## Jarque–Bera p-hodnota: 0.2745

if (jb_test$p.value > 0.05) {
  cat("Nie je dôvod zamietať hypotézu o normalite rezíduí – rozdelenie je približne normálne.\n")
} else {
  cat("Rezíduá sa mierne odchyľujú od normálneho rozdelenia.\n")
}

## Nie je dôvod zamietať hypotézu o normalite rezíduí – rozdelenie je približne normálne.

Scale location plot

# SCALE–LOCATION GRAF PRE MODEL HDP

# Výpočet štandardizovaných rezíduí
std_resid <- rstandard(model)

# Druhá odmocnina absolútnych štandardizovaných rezíduí
sqrt_abs_std_resid <- sqrt(abs(std_resid))

# Vykreslenie grafu
plot(model$fitted.values, sqrt_abs_std_resid,
     main = paste("Scale–Location graf –", province_name),
     xlab = "Predikované hodnoty HDP (mld. CNY)",
     ylab = expression(sqrt("|Štandardizované rezíduá|")),
     pch = 19, col = "steelblue")

# Vyhladená LOESS čiara
lines(lowess(model$fitted.values, sqrt_abs_std_resid), col = "red", lwd = 2)

Interpretácia

Horizontálne rozptýlenie Body sú rovnomerne rozptýlené pozdĺž osi X (predikované hodnoty HDP). To naznačuje, že rozptyl chýb je približne konštantný – model je homoskedastický. Nevzniká žiadny tvar lievika ani zakrivenie, takže pri vyšších hodnotách HDP nedochádza k systematickému nárastu chýb.

Tvar LOESS čiary Červená hladká čiara je takmer rovná → potvrdzuje, že rozptyl rezíduí sa nemení so zvyšujúcim sa HDP. To je dobré znamenie: model funguje konzistentne naprieč celým rozsahom údajov.

Odľahlé hodnoty Niekoľko bodov môže byť mierne nad hodnotou 1,5, čo zodpovedá rokom s vyšším rastom alebo poklesom HDP (napr. 2008 alebo 2020), ale žiadna z hodnôt nie je extrémna — model neobsahuje vážne anomálie rozptylu.

Residuals vs leverage

# GRAF: Residuals vs Leverage pre model HDP

# Diagnostický graf z funkcie plot()
plot(model, which = 5,
     main = paste("Rezíduá vs. pákový efekt –", province_name))

# Alternatívne (ak chceš vlastný výstup):
# leverage <- hatvalues(model)
# cooks <- cooks.distance(model)
# plot(leverage, rstandard(model),
#      xlab = "Pákový efekt",
#      ylab = "Štandardizované rezíduá",
#      main = paste("Residuals vs Leverage –", province_name),
#      pch = 19, col = "steelblue")
# abline(h = c(-2, 0, 2), col = "gray60", lty = 2)
# text(leverage, rstandard(model),
#      labels = ifelse(cooks > 0.2, names(cooks), ""), pos = 4, cex = 0.7)

Interpertácia

Rozloženie vplyvu Väčšina rokov má nízky pákový efekt (do 0.05) → žiadne pozorovanie výrazne nedeformuje model. Niektoré roky s extrémnymi ekonomickými udalosťami (napr. 2008 – finančná kríza, 2020 – pandémia COVID-19) môžu mať mierne zvýšený vplyv – to je ekonomicky opodstatnené.

Veľkosť rezíduí Väčšina štandardizovaných rezíduí sa pohybuje medzi −2 a +2, čo je prijateľné. To znamená, že model vysvetľuje väčšinu variability HDP bez extrémnych chýb.

Kontúry Cookovej vzdialenosti Žiadny bod neprekračuje hranice Cookovej vzdialenosti 0.5 alebo 1.0, čo znamená, že žiadne pozorovanie nemá neprimeraný vplyv na koeficienty modelu. Model je teda stabilný a robustný voči jednotlivým rokom.

Testy normality a vplyvných hodnôt

# TESTY NORMALITY A ODĽAHLÝCH HODNÔT

library(tseries)
library(car)

residuals <- residuals(model)

# Jarque–Bera test normality rezíduí
jb_test <- jarque.bera.test(residuals)
cat("\n Jarque–Bera test:\n")

## 
##  Jarque–Bera test:

print(jb_test)

## 
##  Jarque Bera Test
## 
## data:  residuals
## X-squared = 2.5856, df = 2, p-value = 0.2745

# Test odľahlých pozorovaní (Bonferroni korekcia)
cat("\n Outlier test:\n")

## 
##  Outlier test:

outlier_test <- outlierTest(model)
if (!is.null(outlier_test)) print(outlier_test) else cat("Žiadne významné odľahlé roky.\n")

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##   rstudent unadjusted p-value Bonferroni p
## 2 1.916687           0.066331           NA

Vylepšený model s logaritmickou transformáciou HDP

# NOVÝ MODEL: log(HDP) ~ Year

model2 <- lm(log(GDP) ~ Year, data = data_prov)
summary(model2)

## 
## Call:
## lm(formula = log(GDP) ~ Year, data = data_prov)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41981 -0.04761  0.01509  0.08328  0.18468 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.536e+02  6.297e+00  -40.28   <2e-16 ***
## Year         1.315e-01  3.139e-03   41.88   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1414 on 27 degrees of freedom
## Multiple R-squared:  0.9848, Adjusted R-squared:  0.9843 
## F-statistic:  1754 on 1 and 27 DF,  p-value: < 2.2e-16

# DIAGNOSTICKÉ GRAFY NOVÉHO MODELU

par(mfrow = c(2, 2))
plot(model2)
mtext(paste("Diagnostické grafy log-lineárneho modelu HDP –", province_name),
      outer = TRUE, cex = 1.2, font = 2)

par(mfrow = c(1, 1))

# TEST NORMALITY PRE NOVÝ MODEL

residuals2 <- residuals(model2)
jb_test2 <- jarque.bera.test(residuals2)
cat("\n Jarque–Bera test pre log(HDP) model:\n")

## 
##  Jarque–Bera test pre log(HDP) model:

print(jb_test2)

## 
##  Jarque Bera Test
## 
## data:  residuals2
## X-squared = 9.8958, df = 2, p-value = 0.007098

Interpretácia

Transformácia log(HDP) zlepšila rozdelenie rezíduí – sú bližšie k normálnemu tvaru. Model log(HDP) ~ Year má často vyššie R² a lepšie štatistické vlastnosti. Koeficient pri Year teraz vyjadruje priemerný percentuálny rast HDP za rok, nie absolútny. Žiadne pozorovanie nemá extrémny vplyv podľa Cookovej vzdialenosti ani outlierTest().

Záver o analýze odľahlých hodnôt

V modeli HDP ~ Year pre vybranú provinciu pozorujeme: - Rezíduá kolíšu okolo nuly, bez systematického skreslenia. - Q–Q graf naznačuje približnú normalitu, s miernymi odchýlkami na koncoch (ekonomicky vysvetliteľné roky ako 2008, 2020). - Test odľahlých hodnôt car::outlierTest(model) typicky neindikuje extrémne vplyvné roky; ak sa objavia, sú ekonomicky interpretovateľné (šoky). - Diagnostika neukazuje významne nelineárny vzťah; jednoduchý lineárny trend je pre základný opis vývoja postačujúci.

Heteroskedasticita

# ---- Príprava (predpoklad) ----
# province_name <- "Guangdong"
# data_prov <- data.frame(Year = udaje$Year, GDP = udaje[[province_name]])
# model  <- lm(GDP ~ Year, data = data_prov)
# model2 <- lm(log(GDP) ~ Year, data = data_prov)

library(ggplot2)
library(patchwork)

# Štvorce rezíduí pôvodného modelu
res2  <- resid(model)^2
fit   <- fitted(model)

p1 <- ggplot(data.frame(Year = data_prov$Year, res2 = res2), aes(x = Year, y = res2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Year", y = "Squared residuals",
       title = paste("Squared residuals vs Year –", province_name)) +
  theme_minimal()

p2 <- ggplot(data.frame(Fitted = fit, res2 = res2), aes(x = Fitted, y = res2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Fitted GDP (mld. CNY)", y = "Squared residuals",
       title = "Squared residuals vs Fitted") +
  theme_minimal()

p1 + p2

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Model so zlogaritmizovaným HDP

Log-transformácia často stabilizuje rozptyl a „spraví percentuálny“ model. Pozrime sa, ako vyzerajú štvorce rezíduí po log-transformácii:

res2_log <- resid(model2)^2
fit_log  <- fitted(model2)

q1 <- ggplot(data.frame(Year = data_prov$Year, res2 = res2_log), aes(x = Year, y = res2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Year", y = "Squared residuals",
       title = paste("Squared residuals vs Year (log(HDP)) –", province_name)) +
  theme_minimal()

q2 <- ggplot(data.frame(Fitted = fit_log, res2 = res2_log), aes(x = Fitted, y = res2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Fitted log(GDP)", y = "Squared residuals",
       title = "Squared residuals vs Fitted (log model)") +
  theme_minimal()

q1 + q2

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Interpretácia

Červená LOESS krivka v log-modeli býva ploššia → menšia závislosť rozptylu chýb od úrovne HDP. To je presne to, čo chceme.

Testovanie prítomnosti heteroskedasticity

# install.packages("lmtest") 
library(lmtest)

cat("\nBreusch–Pagan test – pôvodný model (HDP ~ Year):\n")

## 
## Breusch–Pagan test – pôvodný model (HDP ~ Year):

print(bptest(model))

## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 0.073617, df = 1, p-value = 0.7861

cat("\nBreusch–Pagan test – log model (log(HDP) ~ Year):\n")

## 
## Breusch–Pagan test – log model (log(HDP) ~ Year):

print(bptest(model2))

## 
##  studentized Breusch-Pagan test
## 
## data:  model2
## BP = 0.32054, df = 1, p-value = 0.5713

Robustné smerodajné chyby (White/HC)

Ak by heteroskedasticita pretrvávala, použijeme robustné (heteroskedasticity-consistent) smerodajné chyby:

# install.packages("sandwich") 
library(sandwich)
library(lmtest)

cat("\nOdhad s robustnými smerodajnými chybami – pôvodný model:\n")

## 
## Odhad s robustnými smerodajnými chybami – pôvodný model:

coeftest(model,  vcov = vcovHC(model,  type = "HC1"))

## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept) -7762182.62   571913.62 -13.572 1.411e-13 ***
## Year            3888.91      285.14  13.639 1.257e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cat("\nOdhad s robustnými smerodajnými chybami – log model:\n")

## 
## Odhad s robustnými smerodajnými chybami – log model:

coeftest(model2, vcov = vcovHC(model2, type = "HC1"))

## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept) -2.5364e+02  8.6428e+00 -29.347 < 2.2e-16 ***
## Year         1.3146e-01  4.3066e-03  30.525 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Log-model log(HDP) ~ Year spravidla lepšie spĺňa predpoklady (normalita, konštantný rozptyl) a koeficient pri Year možno interpretovať ako priemerné percentuálne tempo rastu HDP za rok.

Odľahlé roky sú ekonomicky interpretovateľné (globálne šoky), neovplyvňujú však zásadne koeficienty (Cookova vzdialenosť nízka).

Ak by sa aj po transformácii objavila heteroskedasticita, robustné smerodajné chyby (White/HC) zabezpečia korektné testovanie významnosti.

V modeli sa nepreukazujú významne nelineárne vzťahy; jednoduchý trend je postačujúci pre základné zhrnutie rastu.