Úvod do problému a stanovenie hypotéz

Budeme modelovať vývoj HDP jednotlivých provincií Číny v závislosti od roku. Zameriame sa na skúmanie trendu rastu a rozdielov medzi regiónmi.

Pracovná hypotéza:

-> HDP všetkých provincií rastie v čase (očakávame kladný trend) -> Rýchlosť rastu sa však líši medzi regiónmi (napr. Guangdong vs. Tibet)

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

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

# Z údajov si vyberieme obdobie po roku 2000 (modernejšia ekonomika)
udaje.modern <- udaje[udaje$Year >= 2000, ]

# Základné čistenie - nahradíme chýbajúce hodnoty mediánom každej provincie
column_medians <- sapply(udaje.modern[,-1], median, na.rm = TRUE)
udaje_imputed <- udaje.modern

for (col in names(udaje.modern)[-1]) {
  udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]
}

udaje.modern <- udaje_imputed

Analýza a vizualizácia

# Načítame potrebné knižnice
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)

# Premeníme dáta do long formátu (rok, provincia, HDP)
udaje_long <- udaje.modern %>%
  pivot_longer(cols = -Year, names_to = "Province", values_to = "GDP")

# Pridáme nový ukazovateľ - medziročný rast HDP (%)
udaje_long <- udaje_long %>%
  group_by(Province) %>%
  arrange(Year) %>%
  mutate(GDP_growth = 100 * (GDP - lag(GDP)) / lag(GDP))

Ranking provincií podľa rastu HDP a vizualizácia

# Priemerný medziročný rast HDP za celé obdobie
avg_growth <- udaje_long %>%
  group_by(Province) %>%
  summarise(AverageGrowth = mean(GDP_growth, na.rm = TRUE)) %>%
  arrange(desc(AverageGrowth))

# Top 10 provincií s najrýchlejším rastom
top10 <- head(avg_growth, 10)

# Vizualizácia: Top 10 rastúcich provincií
ggplot(top10, aes(x = reorder(Province, AverageGrowth), y = AverageGrowth, fill = Province)) +
  geom_col(show.legend = FALSE) +
  coord_flip() +
  labs(title = "Top 10 čínskych provincií podľa priemerného rastu HDP (2000–2020)",
       x = "Provincia", y = "Priemerný medziročný rast HDP (%)") +
  theme_minimal(base_size = 13)

Modelovanie trendu pre vybranú provinciu

# Vyberme napríklad Guangdong (ekonomický líder)
guangdong <- udaje_long %>% filter(Province == "Guangdong")

# Jednoduchý lineárny model rastu HDP ~ rok
model <- lm(GDP ~ Year, data = guangdong)
summary(model)
## 
## Call:
## lm(formula = GDP ~ Year, data = guangdong)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -6759  -4303  -2994   4132  11842 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.050e+07  4.198e+05  -25.02 5.26e-16 ***
## Year         5.250e+03  2.088e+02   25.14 4.80e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5795 on 19 degrees of freedom
## Multiple R-squared:  0.9708, Adjusted R-squared:  0.9693 
## F-statistic: 632.1 on 1 and 19 DF,  p-value: 4.801e-16
# Vizualizácia trendu
ggplot(guangdong, aes(x = Year, y = GDP)) +
  geom_line(color = "#0072B2", size = 1.2) +
  geom_smooth(method = "lm", se = TRUE, color = "#E69F00") +
  labs(title = "Trend HDP provincie Guangdong (2000–2020)",
       subtitle = paste0("Odhadovaný ročný rast: ",
                         round(coef(model)[2], 1), " mld. CNY / rok"),
       x = "Rok", y = "HDP (mld. CNY)") +
  theme_minimal(base_size = 13)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `geom_smooth()` using formula = 'y ~ x'

Vizualizácia tvaru údajov a kontrola možných nezrovnalostí

Predpokladáme, že máme dátový rámec udaje.modern (rok + provincie) alebo jeho verziu po imputácii chýbajúcich údajov

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

# Zistíme globálny rozsah hodnôt HDP vo všetkých provinciách
global_range <- range(udaje_numeric, na.rm = TRUE)

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

# Nastavíme rozloženie grafov – 4x4, aby sme videli 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]]
  
  # Boxplot pre danú provinciu
  boxplot(data_col,
          main = col,
          col = ifelse(any(data_col == 0, na.rm = TRUE), "tomato", "lightblue"),
          ylab = "HDP (mld. CNY)",
          cex.main = 0.8,
          ylim = global_range)  # <<< rovnaký rozsah pre všetky
}

# Globálny titulok pre všetky grafy
mtext("Boxploty HDP jednotlivých provincií Číny (2000–2020)", 
      outer = TRUE, cex = 1.3, font = 2)

# Obnovíme pôvodné nastavenie grafov
par(mfrow = c(1, 1))

## Kontrola nulových hodnôt v dataset-e

# Počet nulových hodnôt v každej provincii
nulove <- sapply(udaje_numeric, function(x) sum(x == 0, na.rm = TRUE))

# Zobraz provincie s nulami (ak existujú)
nulove_df <- data.frame(Provincia = names(nulove),
                        PocetNul = nulove)

# Zobrazíme len provincie, kde sa nuly vyskytli
nulove_df <- subset(nulove_df, PocetNul > 0)

if(nrow(nulove_df) > 0) {
  print("Provinicie 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 – model rastu HDP v čase

# Vyberieme si konkrétnu provinciu, napr. Guangdong (ekonomický líder)
province_name <- "Guangdong"

# Vytvoríme podmnožinu dát pre danú provinciu
data_prov <- data.frame(
  Year = udaje.modern$Year,
  GDP = udaje.modern[[province_name]]
)

# Skontrolujeme, či tam nie sú chýbajúce hodnoty
sum(is.na(data_prov$GDP))
## [1] 0
# Odhadneme lineárny model HDP ~ rok
model <- lm(GDP ~ Year, data = data_prov)

# Výpis základných výsledkov
cat("✅ Odhadnuté koeficienty modelu pre provinciu", province_name, ":\n")
## ✅ Odhadnuté koeficienty modelu pre provinciu Guangdong :
print(model$coefficients)
##   (Intercept)          Year 
## -10502068.979      5250.519
cat("\nPriemerný medziročný rast HDP:", round(model$coefficients[2], 2), "mld. CNY\n")
## 
## Priemerný medziročný rast HDP: 5250.52 mld. CNY
# Zhrnutie modelu (t-statistiky, R² atď.)
summary(model)
## 
## Call:
## lm(formula = GDP ~ Year, data = data_prov)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -6759  -4303  -2994   4132  11842 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.050e+07  4.198e+05  -25.02 5.26e-16 ***
## Year         5.250e+03  2.088e+02   25.14 4.80e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5795 on 19 degrees of freedom
## Multiple R-squared:  0.9708, Adjusted R-squared:  0.9693 
## F-statistic: 632.1 on 1 and 19 DF,  p-value: 4.801e-16

VIZUALIZÁCIA REGRESNÉHO MODELU

# Nastavíme grafické parametre
par(mfrow = c(1, 1))
plot(data_prov$Year, data_prov$GDP,
     main = paste("Lineárna regresia rastu HDP -", province_name),
     xlab = "Rok", ylab = "HDP (mld. CNY)",
     pch = 16, col = "skyblue")

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

# Doplníme text s R²
r2 <- summary(model)$r.squared
legend("topleft", legend = paste("R² =", round(r2, 3)),
       bty = "n", text.col = "red")

DIAGNOSTICKÉ GRAFY MODELU

par(mfrow = c(2, 2))
plot(model)
mtext(paste("Diagnostické grafy pre model provincie", province_name),
      outer = TRUE, cex = 1.2, font = 2)

par(mfrow = c(1, 1))

Predikcia HDP do roku 2025

# Vytvoríme nový dataset s budúcimi rokmi
future_years <- data.frame(Year = seq(2021, 2025, by = 1))

# Predikcia pomocou modelu
pred <- predict(model, newdata = future_years, interval = "prediction")

# Spojíme s rokmi
pred_df <- data.frame(future_years, pred)

cat("\n📈 Odhadovaný HDP provincie", province_name, "do roku 2025:\n")
## 
## 📈 Odhadovaný HDP provincie Guangdong do roku 2025:
print(pred_df)
##   Year      fit       lwr      upr
## 1 2021 109229.6  95916.15 122543.0
## 2 2022 114480.1 101002.64 127957.6
## 3 2023 119730.6 106077.10 133384.1
## 4 2024 124981.1 111139.99 138822.3
## 5 2025 130231.7 116191.77 144271.5
# Vizualizácia predikcie
plot(data_prov$Year, data_prov$GDP,
     main = paste("Predikcia HDP do roku 2025 –", province_name),
     xlab = "Rok", ylab = "HDP (mld. CNY)",
     pch = 16, col = "skyblue", xlim = c(2000, 2025))

abline(model, col = "red", lwd = 2)
lines(pred_df$Year, pred_df$fit, col = "darkgreen", lwd = 2, lty = 2)
lines(pred_df$Year, pred_df$lwr, col = "gray50", lty = 3)
lines(pred_df$Year, pred_df$upr, col = "gray50", lty = 3)

legend("topleft", legend = c("Skutočné údaje", "Regresná priamka", "Predikcia 2021–2025"),
       col = c("skyblue", "red", "darkgreen"), pch = c(16, NA, NA), lty = c(NA, 1, 2),
       bty = "n")

Residuals vs. Fitted

# Vykreslenie základného diagnostického grafu
plot(model$fitted.values, model$residuals,
     main = paste("Rezíduá vs. prispôsobené hodnoty –", province_name),
     xlab = "Predikované hodnoty HDP (mld. CNY)",
     ylab = "Rezíduá (skutočné - predikované)",
     pch = 19, col = "steelblue")
abline(h = 0, col = "red", lwd = 2)

# Pridáme hladkú čiaru LOESS, ktorá pomáha odhaliť nelinearity
lines(lowess(model$fitted.values, model$residuals), col = "orange", lwd = 2)

Interpretácia

Rezíduá oscilujú okolo nulovej osi → model nie je systematicky skreslený. Oranžová LOESS krivka je relatívne plochá → lineárny model je postačujúci. Vertikálny rozptyl rezíduí je približne rovnaký → nepozorujeme heteroskedasticitu. Ak by sa body rozširovali do tvaru lievika, znamenalo by to problém s variabilitou chýb. Niekoľko bodov výrazne mimo centrálneho pásma môže byť potenciálne odľahlé pozorovanie.

Q-Q plot

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

# Test normality rezíduí
shapiro_result <- shapiro.test(model$residuals)

Interpretácia Q–Q grafu

Väčšina bodov leží blízko červenej čiary → rezíduá sú približne normálne rozložené. Mierne odchýlky na okrajoch (vľavo dole a vpravo hore) naznačujú len malé odchýlky od normálneho tvaru. Predpoklad normality je teda približne splnený.

SCALE–LOCATION PLOT (variabilita rezíduí)

# Vypočítame štandardizované rezíduá
std_resid <- rstandard(model)

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

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 = "skyblue")
lines(lowess(model$fitted.values, sqrt_abs_std_resid), col = "red", lwd = 2)

Interpretácia grafu Scale–Location

Body sú približne rovnomerne rozptýlené po osi X – variancia rezíduí sa javí ako konštantná (homoskedasticita). Červená LOESS čiara je relatívne plochá – čo potvrdzuje stabilný rozptyl chýb. Niekoľko bodov s mierne vyššími hodnotami (nad 1,5) môže indikovať roky s atypickým rastom HDP, ale bez extrémnych odchýlok.

RESIDUALS VS LEVERAGE (vplyv pozorovaní)

plot(model, which = 5, main = paste("Rezíduá vs. Pákový efekt –", province_name))

Interpretácia grafu Residuals vs. Leverage

Väčšina rokov má nízky pákový efekt (<0.05), čo znamená, že jednotlivé roky výrazne neovplyvňujú tvar regresie. Niekoľko bodov vpravo s vyššou hodnotou leverage (napr. po roku 2015) môže mať väčší vplyv na sklon regresnej priamky – ide o roky s neobvykle vysokým HDP. Bodkované kontúry Cookovej vzdialenosti ukazujú, že žiadne pozorovanie výrazne neprekračuje prah 0.5 → žiadny extrémny vplyv.

TESTY NORMALITY A ODĽAHLÝCH HODNÔT

library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
residuals <- residuals(model)

# Jarque–Bera test normality
jb_test <- jarque.bera.test(residuals)
cat("\n📊 Jarque–Bera test normality:\n")
## 
## 📊 Jarque–Bera test normality:
print(jb_test)
## 
##  Jarque Bera Test
## 
## data:  residuals
## X-squared = 2.4675, df = 2, p-value = 0.2912
# Outlier test
outlier_test <- outlierTest(model)
cat("\n📊 Test odľahlých pozorovaní (Bonferroni):\n")
## 
## 📊 Test odľahlých pozorovaní (Bonferroni):
if (!is.null(outlier_test)) print(outlier_test) else cat("✅ Žiadne významné odľahlé hodnoty.\n")
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##    rstudent unadjusted p-value Bonferroni p
## 21 2.561709           0.019613      0.41188

NOVÝ MODEL S LOGARITMICKOU TRANSFORMÁCIOU HDP

Ak sa ukáže, že rezíduá nie sú normálne rozdelené alebo majú heteroskedasticitu, transformujeme HDP pomocou logaritmu. Tým dostaneme model, ktorý lepšie vystihuje percentuálny rast HDP.

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.203248 -0.081222 -0.005315  0.096480  0.150612 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.327e+02  7.602e+00  -30.62   <2e-16 ***
## Year         1.211e-01  3.782e-03   32.01   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1049 on 19 degrees of freedom
## Multiple R-squared:  0.9818, Adjusted R-squared:  0.9808 
## F-statistic:  1025 on 1 and 19 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
jb_test2 <- jarque.bera.test(residuals(model2))
cat("\n📈 Jarque–Bera test pre log-transformovaný model:\n")
## 
## 📈 Jarque–Bera test pre log-transformovaný model:
print(jb_test2)
## 
##  Jarque Bera Test
## 
## data:  residuals(model2)
## X-squared = 0.98665, df = 2, p-value = 0.6106

Interpretácia výsledkov

Interpretácia modelu

Lineárny model HDP ~ Year vystihuje základný trend rastu v čase. Rezíduá oscilujú okolo nuly a nepreukazujú výraznú heteroskedasticitu. Predpoklad normality je približne splnený (ak p > 0.05 v JB teste). Odľahlé roky (napr. prudké skoky HDP v 2008 alebo 2020) môžu mať menší vplyv, ale neovplyvňujú výrazne celkový trend.

Interpretácia modelu 2 (log(HDP))

Po logaritmickej transformácii sa rozdelenie rezíduí približuje normalite. Koeficient pri Year teraz vyjadruje priemerný percentuálny rast HDP (elasticitu). Model vysvetľuje časový trend ešte lepšie, čo vidno aj na zlepšení R² a stabilnejších rezíduách.