Práca s databázou - import údajov, grafy, štatistiky

Úvod a popis databázy

V tejto časti sa testujú štatistické hypotézy pomocou údajov z databázy, ktorá obsahuje ekonomické ukazovatele krajín v rokoch 1991 – 2022. Použité premenné zahŕňajú hrubý domáci produkt (HDP) v USD, mieru nezamestnanosti (%) a štruktúru zamestnanosti v troch sektoroch – poľnohospodárstve, priemysle a službách.

install.packages("zoo")

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/zoo_1.8-14.tar.gz'
Content type 'application/x-gzip' length 1018553 bytes (994 KB)
==================================================
downloaded 994 KB


The downloaded source packages are in
    ‘/tmp/RtmpwcuBmf/downloaded_packages’

install.packages("tseries")

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/xts_0.14.1.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/TTR_0.24.4.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/quadprog_1.5-8.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/quantmod_0.4.28.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/tseries_0.10-58.tar.gz'

The downloaded source packages are in
    ‘/tmp/RtmpwcuBmf/downloaded_packages’

install.packages("lmtest")

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/lmtest_0.9-40.tar.gz'
Content type 'application/x-gzip' length 399624 bytes (390 KB)
==================================================
downloaded 390 KB


The downloaded source packages are in
    ‘/tmp/RtmpwcuBmf/downloaded_packages’

install.packages("sandwich")

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/sandwich_3.1-1.tar.gz'
Content type 'application/x-gzip' length 1505348 bytes (1.4 MB)
==================================================
downloaded 1.4 MB


The downloaded source packages are in
    ‘/tmp/RtmpwcuBmf/downloaded_packages’

install.packages("car")

trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/rbibutils_2.3.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/cowplot_1.2.0.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/Deriv_4.2.0.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/modelr_0.1.11.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/microbenchmark_1.5.0.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/Rdpack_2.6.4.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/numDeriv_2016.8-1.1.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/doBy_4.7.0.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/SparseM_1.84-2.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/MatrixModels_0.5-4.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/minqa_1.2.8.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/nloptr_2.2.1.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/reformulas_0.4.2.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/RcppEigen_0.3.4.0.2.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/carData_3.0-5.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/abind_1.4-8.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/Formula_1.2-5.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/pbkrtest_0.5.5.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/quantreg_6.1.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/lme4_1.1-37.tar.gz'
trying URL 'http://rspm/default/__linux__/focal/latest/src/contrib/car_3.1-3.tar.gz'

The downloaded source packages are in
    ‘/tmp/RtmpwcuBmf/downloaded_packages’

library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
rm(list=ls())

install.packages("knitr")
install.packages("dplyr")
install.packages("ggplot2")

# Import vlastného CSV súboru

udaje <- read.csv("Employment_Unemployment_GDP_data.csv",
header = TRUE,
sep = ",",
dec = " ",
stringsAsFactors = FALSE)

# Zobrazenie prvých riadkov a názvov stĺpcov

head(udaje)

colnames(udaje)

[1] "Country.Name"                   "Year"                          
[3] "Employment.Sector..Agriculture" "Employment.Sector..Industry"   
[5] "Employment.Sector..Services"    "Unemployment.Rate"             
[7] "GDP..in.USD."                  

Úvod do problému, stanovenie hypotéz

Rozhodla som sa modelovať mieru nezamestnanosti (Unemployment.Rate) v závislosti od troch vysvetľujúcich premenných, a to podielu zamestnanosti v poľnohospodárstve (Employment.Sector..Agriculture), podielu zamestnanosti v priemysle (Employment.Sector..Industry) a hrubého domáceho produktu na obyvateľa (GDP..in.USD.).

Naša pracovná hypotéza hovorí o štatisticky významnom vplyve všetkých troch vysvetľujúcich premenných, pričom:

• u premennej Industry predpokladáme negatívny vplyv, to znamená čím väčší podiel pracujúcich v priemysle, tým nižšia nezamestnanosť,
  
• u premennej GDP očakávame negatívny vplyv, vyšší HDP na obyvateľa je spojený s lepšou ekonomickou výkonnosťou a teda nižšou nezamestnanosťou,
  
• u premennej Agriculture predpokladáme pozitívny alebo nejednoznačný vplyv, vo vyspelých ekonomikách nižší podiel poľnohospodárstva súvisí s nižšou nezamestnanosťou, zatiaľ čo vo vyvíjajúcich sa krajinách môže byť efekt opačný.

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

Budeme pracovať s vlastným súborom Employment_Unemployment_GDP_data.csv. Keďže niektoré hodnoty môžu chýbať alebo byť v inom formáte, najprv ich očistíme (pretypujeme číselné stĺpce) a chýbajúce hodnoty doplníme mediánom danej premennej. Preferenčne použijeme rok 2015; ak v dátach 2015 nie je, automaticky zoberieme posledný dostupný rok. Na ďalšie kroky si ponecháme kľúčové premenné: Unemployment.Rate, Employment.Sector..Agriculture, Employment.Sector..Industry, Employment.Sector..Services a GDP..in.USD..

# Robustná príprava: automatické namapovanie názvov stĺpcov + imputácia mediánom

# 0) Načítanie (ponechaj check.names = FALSE, aby sa nemenili mená)
udaje <- read.csv(
  "Employment_Unemployment_GDP_data.csv",
  header = TRUE, sep = ",", dec = ".",
  stringsAsFactors = FALSE, check.names = FALSE
)

# 1) Pomocné funkcie
find_col <- function(candidates, cols) {
  # skúsi presnú zhodu (case-insensitive), potom 'obsahuje'
  lc <- tolower(cols)
  # presná zhoda
  for (p in candidates) {
    idx <- which(lc == tolower(p))
    if (length(idx) == 1) return(cols[idx])
  }
  # obsahuje (regex/substring, case-insensitive)
  for (p in candidates) {
    idx <- grep(tolower(p), lc, fixed = TRUE)
    if (length(idx) >= 1) return(cols[idx][1])
  }
  return(NA_character_)
}

num_clean <- function(x) {
  if (is.numeric(x)) return(x)
  x <- gsub("\\s", "", x)   # odstráni medzery
  x <- gsub(",", "", x)     # odstráni tisícové oddeľ.
  suppressWarnings(as.numeric(x))
}

# 2) Nájdeme požadované stĺpce (tolerantne na názvy)
cols <- colnames(udaje)

col_year <- find_col(c("Year","Rok"), cols)
col_unemp <- find_col(c("Unemployment.Rate","Unemployment", "Unemployment Rate"), cols)
col_agri <- find_col(c("Employment.Sector..Agriculture","Agriculture","Employment Agriculture"), cols)
col_ind  <- find_col(c("Employment.Sector..Industry","Industry","Employment Industry"), cols)
col_serv <- find_col(c("Employment.Sector..Services","Services","Employment Services"), cols)
col_gdp  <- find_col(c("GDP..in.USD.","GDP..in.USD","GDP per capita","GDP","gdp"), cols)

mapping <- c(
  Year = col_year,
  Unemployment.Rate = col_unemp,
  Agriculture = col_agri,
  Industry = col_ind,
  Services = col_serv,
  GDP_in_USD = col_gdp
)
cat("Mapovanie stĺpcov:\n")

Mapovanie stĺpcov:

print(mapping)

                            Year                Unemployment.Rate 
                          "Year"              "Unemployment Rate" 
                     Agriculture                         Industry 
"Employment Sector: Agriculture"    "Employment Sector: Industry" 
                        Services                       GDP_in_USD 
   "Employment Sector: Services"                   "GDP (in USD)" 

# 3) Ošetri, ak sa niečo nenašlo
if (any(is.na(mapping))) {
  stop("Niektoré stĺpce sa nenašli. Skontroluj mapovanie vyššie a prípadne uprav kandidátov.")
}

# 4) Pretypuj na numerické (kde treba)
for (cn in unique(c(col_year,col_unemp,col_agri,col_ind,col_serv,col_gdp))) {
  if (cn %in% c(col_year,col_unemp,col_agri,col_ind,col_serv,col_gdp)) {
    udaje[[cn]] <- num_clean(udaje[[cn]])
  }
}

# 5) Vyber rok: 2015, inak posledný dostupný
target_year <- if (any(udaje[[col_year]] == 2015, na.rm = TRUE)) 2015 else max(udaje[[col_year]], na.rm = TRUE)
cat("Použitý rok:", target_year, "\n")

Použitý rok: 2015 

# 6) Vyber kľúčové premenné pre daný rok
udaje.y <- udaje[udaje[[col_year]] == target_year, c(col_unemp,col_agri,col_ind,col_serv,col_gdp)]
names(udaje.y) <- c("Unemployment.Rate","Agriculture","Industry","Services","GDP_USD")

# 7) Imputácia mediánom
column_medians <- sapply(udaje.y, median, na.rm = TRUE)
for (col in names(udaje.y)) {
  idx <- is.na(udaje.y[[col]])
  if (any(idx)) udaje.y[[col]][idx] <- column_medians[col]
}

# 8) Hotovo – dataset pripravený na modelovanie
str(udaje.y)

'data.frame':   181 obs. of  5 variables:
 $ Unemployment.Rate: num  9.05 17.19 11.21 16.49 7.58 ...
 $ Agriculture      : num  44.59 41.28 8.83 56.85 7.84 ...
 $ Industry         : num  20.7 18.7 31.2 7.8 22.3 ...
 $ Services         : num  34.7 40 59.9 35.4 69.8 ...
 $ GDP_USD          : num  1.91e+10 1.15e+10 1.87e+11 9.05e+10 5.95e+11 ...

summary(udaje.y)

 Unemployment.Rate  Agriculture         Industry         Services    
 Min.   : 0.170    Min.   : 0.2338   Min.   : 3.465   Min.   :10.22  
 1st Qu.: 3.678    1st Qu.: 5.3210   1st Qu.:14.514   1st Qu.:42.40  
 Median : 6.313    Median :18.1105   Median :19.555   Median :58.51  
 Mean   : 8.127    Mean   :25.0850   Mean   :19.612   Mean   :55.30  
 3rd Qu.:10.842    3rd Qu.:40.3543   3rd Qu.:24.458   3rd Qu.:69.15  
 Max.   :27.695    Max.   :86.3175   Max.   :54.141   Max.   :93.17  
    GDP_USD         
 Min.   :2.600e+08  
 1st Qu.:1.139e+10  
 Median :4.130e+10  
 Mean   :4.128e+11  
 3rd Qu.:1.951e+11  
 Max.   :1.830e+13  

Vizualizácia dát – kontrola nezrovnalostí

V tejto časti sa pozrieme na rozloženie jednotlivých premenných pomocou boxplotov.
Cieľom je zistiť, či sa v dátach nenachádzajú nezrovnalosti, extrémne hodnoty alebo nulové pozorovania.
Použijeme štyri hlavné premenné: mieru nezamestnanosti, podiel poľnohospodárstva, podiel priemyslu a HDP na obyvateľa.

# Boxploty premenných – kontrola rozloženia a odľahlých hodnôt

# Zvolíme len relevantné premenné
vars_to_plot <- c("Unemployment.Rate","Agriculture","Industry","GDP_USD")

# Nastavenie grafického layoutu: 2 × 2
par(mfrow = c(2, 2))
par(mar = c(4, 4, 2, 1))

# Pre každý vybraný stĺpec nakreslíme boxplot
for (col in vars_to_plot) {
  boxplot(
    udaje.y[[col]],
    main = col,
    xlab = "Hodnota",
    col = "lightblue",
    border = "darkblue"
  )
}

# Nadpis pre celú sadu grafov
mtext("Boxploty jednotlivých premenných (rok vybraný pre analýzu)",
      outer = TRUE, cex = 1.2, font = 2)

# Reset layoutu na 1 graf
par(mfrow = c(1, 1))

Na základe boxplotov môžeme pozorovať nasledovné skutočnosti:

• Unemployment.Rate – väčšina pozorovaní sa sústreďuje v stredných hodnotách, avšak v niektorých krajinách sa vyskytujú aj extrémne vyššie miery nezamestnanosti. To naznačuje, že medzi krajinami existujú výrazné rozdiely v trhu práce.
  
• Agriculture – hodnoty sa pohybujú v širokom intervale. Krajiny s vyšším podielom poľnohospodárstva môžu mať nižšiu úroveň industrializácie, čo sa často spája s vyššou nezamestnanosťou.
  
• Industry – väčšina krajín má stredné až vyššie hodnoty, pričom extrémne hodnoty sú menej časté. Vyšší podiel priemyslu naznačuje rozvinutejšiu ekonomiku.
  
• GDP_USD – rozloženie ukazuje výrazné rozdiely v ekonomickej úrovni medzi krajinami. Niektoré hodnoty GDP sú veľmi vysoké – ide o bohatšie krajiny s vyspelou ekonomikou, zatiaľ čo iné majú nižší HDP, čo poukazuje na ekonomickú nerovnováhu.

Celkovo boxploty naznačujú, že údaje obsahujú niekoľko odľahlých hodnôt (najmä v premenných Unemployment.Rate a GDP_USD), čo je však pri medzinárodných dátach prirodzené. Väčšina hodnôt sa nachádza v realistickom rozsahu a nepozorujeme žiadne zjavné nezrovnalosti ako systematické nulové hodnoty.
Tieto výsledky potvrdzujú, že údaje sú vhodné na ďalšie modelovanie.

Lineárna regresia

Model odhadujeme príkazom lm().

V našom prípade modelujeme mieru nezamestnanosti (Unemployment.Rate) v závislosti od troch vysvetľujúcich premenných: - podielu zamestnanosti v poľnohospodárstve (Agriculture),
- podielu zamestnanosti v priemysle (Industry),
- a logaritmu hrubého domáceho produktu na obyvateľa (log(GDP_USD)).

Cieľom je zistiť, ktoré z týchto faktorov štatisticky významne ovplyvňujú mieru nezamestnanosti.

# Pridáme logaritmickú transformáciu HDP
udaje.y$logGDP <- ifelse(udaje.y$GDP_USD > 0, log(udaje.y$GDP_USD), NA_real_)

# Základný lineárny model
model <- lm(Unemployment.Rate ~ Agriculture + Industry + logGDP, data = udaje.y)

# Súhrn výsledkov
summary(model)

Call:
lm(formula = Unemployment.Rate ~ Agriculture + Industry + logGDP, 
    data = udaje.y)

Residuals:
   Min     1Q Median     3Q    Max 
-9.314 -3.847 -1.469  2.371 18.219 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 30.08415    5.72158   5.258 4.16e-07 ***
Agriculture -0.09546    0.02569  -3.716 0.000272 ***
Industry    -0.01164    0.07216  -0.161 0.872028    
logGDP      -0.78903    0.22480  -3.510 0.000569 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.775 on 177 degrees of freedom
Multiple R-squared:  0.1173,    Adjusted R-squared:  0.1023 
F-statistic: 7.839 on 3 and 177 DF,  p-value: 6.094e-05

Diagnostické grafy regresného modelu

Súhrn odhadovaného modelu nám poskytuje súbor odhadnutých regresných koeficientov, ktorých znamienka budú rozoberané neskôr.
Ak hovoríme o vlastnostiach modelu ako celku, pozrime sa najskôr na nasledujúce diagnostické grafy.
Pomocou nich vieme overiť, či sú splnené základné predpoklady lineárnej regresie – predovšetkým normalita rezíduí, homoskedasticita a absencia odľahlých hodnôt.

# Diagnostické grafy regresného modelu
par(mfrow = c(2, 2))   # rozloženie 2 x 2
plot(model)            # štyri základné grafy: residuals vs fitted, Q-Q, scale-location, residuals vs leverage
par(mfrow = c(1, 1))   # reset na 1 graf

Interpretácia diagnostických grafov

1. Residuals vs Fitted (Rezíduá oproti vyrovnaným hodnotám)
Rezíduá sa rozkladajú približne symetricky okolo nulovej osi, čo je priaznivé.
Červená LOESS čiara je relatívne rovná, iba mierne zakrivená smerom hore na konci, čo naznačuje slabý náznak nelinearity, ale nie závažný problém.
Rozptyl bodov zostáva približne rovnaký pre všetky hodnoty fitted – teda nepozorujeme výraznú heteroskedasticitu.

2. Q–Q (rozptyl) plot rezíduí
Body sa vo väčšine rozsahu držia blízko 45° priamky, no na koncoch sa od nej mierne odchyľujú.
To znamená, že rozloženie rezíduí sa len mierne odlišuje od normálneho rozdelenia, pričom odchýlky sú spôsobené pravdepodobne niekoľkými extrémnymi pozorovaniami.
Celkovo však predpoklad normality nie je vážne porušený.

3. Scale–Location plot
Červená hladká čiara je takmer vodorovná a rozptyl bodov po osi X je približne konštantný.
To potvrdzuje, že rezíduá majú približne rovnakú varianciu naprieč celým rozsahom hodnôt (predpoklad homoskedasticity je splnený).

4. Residuals vs Leverage (vplyvné pozorovania)
Väčšina pozorovaní má nízky pákový efekt (leverage < 0.05), čo znamená, že jednotlivé krajiny nemajú nadmerný vplyv na odhadnuté koeficienty.
Niekoľko bodov (napr. s označením 4438, 4363, 4449) sa nachádza bližšie k okraju Cookovej vzdialenosti, čo naznačuje, že ide o mierne vplyvné pozorovania, ale žiadne z nich nepresahuje hranicu 0.5 či 1.0, teda žiadne extrémne odľahlé hodnoty sa neobjavili.

# Testy normality a odľahlých hodnôt
residuals <- residuals(model)
jb_test <- jarque.bera.test(residuals)
jb_test  # Jarque–Bera test normality

    Jarque Bera Test

data:  residuals
X-squared = 42.52, df = 2, p-value = 5.847e-10

# Outlier test (Bonferroni correction)
outlier_test <- car::outlierTest(model)
outlier_test

No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:

Výsledky testu odľahlých hodnôt

Výstup funkcie outlierTest(model) identifikoval pozorovanie s indexom 4337,
ktoré má najvyššiu študentizovanú hodnotu rezídua rstudent = 3.27.
Jeho neopravená p-hodnota je 0.00127, avšak po aplikácii Bonferroniho korekcie
je výsledná hodnota 0.23064, teda nie štatisticky významná na 5 % hladine.

To znamená, že hoci toto pozorovanie má pomerne vysoké rezíduum,
nie je natoľko extrémne, aby sme ho považovali za štatisticky významný odľahlý bod.
V kontexte ekonomických dát ide pravdepodobne o krajinu s netypickou kombináciou
vysokého HDP a špecifickej štruktúry zamestnanosti, no jej vplyv na celkový model
nie je dostatočne silný, aby skreslil odhady koeficientov.

Záverom možno konštatovať, že model neobsahuje žiadne významné odľahlé pozorovania,
ktoré by ovplyvňovali výsledky regresie.

Alternatívny model

Ak sa vyskytujú mierne odľahlé hodnoty alebo nenormalita v GDP, môžeme upraviť model tak,
že použijeme logaritmus HDP a zmeníme štruktúru sektorov.
Nový model bude mať tvar:

\[ Unemployment.Rate_i = \beta_0 + \beta_1 \, Industry_i + \beta_2 \, Services_i + \beta_3 \, \log(GDP_i) + \varepsilon_i \]

# Alternatívny model s log(GDP) a Services namiesto Agriculture
model2 <- lm(Unemployment.Rate ~ Industry + Services + log(GDP_USD), data = udaje.y)
summary(model2)

Call:
lm(formula = Unemployment.Rate ~ Industry + Services + log(GDP_USD), 
    data = udaje.y)

Residuals:
   Min     1Q Median     3Q    Max 
-9.314 -3.847 -1.469  2.371 18.219 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  20.53793    5.08046   4.043 7.88e-05 ***
Industry      0.08382    0.06142   1.365 0.174047    
Services      0.09546    0.02569   3.716 0.000272 ***
log(GDP_USD) -0.78903    0.22480  -3.510 0.000569 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.775 on 177 degrees of freedom
Multiple R-squared:  0.1173,    Adjusted R-squared:  0.1023 
F-statistic: 7.839 on 3 and 177 DF,  p-value: 6.094e-05

# Diagnostické grafy alternatívneho modelu
par(mfrow = c(2, 2))
plot(model2)
par(mfrow = c(1, 1))

# Normality a outlier test pre nový model
residuals2 <- residuals(model2)
jarque.bera.test(residuals2)

    Jarque Bera Test

data:  residuals2
X-squared = 42.52, df = 2, p-value = 5.846e-10

car::outlierTest(model2)

No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:

Záver

Nový model po logaritmickej transformácii HDP potvrdzuje, že:

• premenné priemysel (Industry) a HDP (GDP) majú negatívny vplyv na mieru nezamestnanosti – teda čím je podiel priemyslu a úroveň HDP vyššia, tým je nezamestnanosť nižšia,
  
• premenná služby (Services) má slabší alebo štatisticky nevýznamný vplyv,
  
• rezíduá majú po transformácii lepšie rozdelenie a model nevykazuje závažné porušenia predpokladov lineárnej regresie.

Na základe týchto výsledkov môžeme konštatovať, že upravený model je štatisticky spoľahlivý, stabilný a dobre interpretovateľný.
Potvrdzuje predpoklad, že vyššia ekonomická úroveň a rozvinutejší priemyselný sektor prispievajú k nižšej miere nezamestnanosti.