Cvičenie 5

S využitím databázy World Development Indicators (WDI) Svetovej banky, z ktorej čerpáme:

• Life expectancy at birth, total (years) – stredná dĺžka života
  
• GDP per capita (current US$) – HDP na obyvateľa
  
• Expected years of schooling – očakávané roky štúdia
  
• Prevalence of overweight (% of adults) – podiel dospelých s BMI > 25 (použijeme ako premennú typu BMI).

Pri ďalšej práci budeme používať knižnice:

V hlavnom menu môžeme Session nastaviť na Source File Location,
alebo to urobiť cez Console, napr.:

setwd("/Cloud/project/tyzdne/tyzden5")

Organizácia priečinkov a podpriečinkov sa môže líšiť podľa projektu – v tomto príklade dáta nesťahujeme z lokálneho CSV, ale priamo z WDI API Svetovej banky.

Úvod do problému, stanovenie hypotéz

Rozhodol som sa modelovať strednú dĺžku života Life.expectancy v závislosti od troch vysvetľujúcich premenných, a to:

• BMI (v našom prípade: prevalencia nadváhy v % dospelých),
• HDP na obyvateľa GDP,
• stredný počet rokov štúdia Schooling.

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

• u premenných GDP a Schooling očakávame pozitívny vplyv (kladné znamienko koeficienta),
• v prípade BMI očakávame negatívny vplyv (záporné znamienko koeficienta).

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

Dáta získame z WDI pre rok 2015 a následne doplníme chýbajúce hodnoty mediánmi.

# Definícia indikátorov z WDI
indicators <- c(
  "SP.DYN.LE00.IN",  # Life expectancy at birth, total (years)
  "NY.GDP.PCAP.CD",  # GDP per capita (current US$)
  "SE.SCH.LIFE",     # Expected years of schooling
  "SH.STA.OWAD.ZS"   # Prevalence of overweight (% of adults)
)

# Stiahnutie údajov pre všetky krajiny za rok 2015
wdi_raw <- WDI(
  country   = "all",
  indicator = indicators,
  start     = 2015,
  end       = 2015,
  extra     = FALSE
)

# Vyberieme len rok 2015 (pre istotu) a príslušné stĺpce
udaje_2015 <- subset(wdi_raw, year == 2015)

udaje_2015 <- udaje_2015[, c(
  "SP.DYN.LE00.IN",
  "SH.STA.OWAD.ZS",
  "NY.GDP.PCAP.CD",
  "SE.SCH.LIFE"
)]

# Premenujeme stĺpce na jednoduchšie názvy
names(udaje_2015) <- c("Life.expectancy", "BMI", "GDP", "Schooling")

# Imputácia chýbajúcich hodnôt mediánom danej premennej
column_medians <- sapply(udaje_2015, median, na.rm = TRUE)

udaje_imputed <- udaje_2015

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

udaje_2015 <- udaje_imputed

str(udaje_2015)
summary(udaje_2015)

Teraz chceme vidieť tvar údajov (či nie sú v nich nejaké nezrovnalosti – napríklad extrémne hodnoty).

# Počet premenných
num_plots <- length(names(udaje_2015))

# Nastavenie layoutu: 2 riadky × 2 stĺpce
par(mfrow = c(2, 2))
par(mar = c(4, 4, 2, 1))

# Vykreslenie boxplotu pre každú premennú
for (col in names(udaje_2015)) {
  boxplot(
    udaje_2015[[col]],
    main = col,
    xlab = "Hodnota",
    col = "lightblue"
  )
}

# Hlavný nadpis
mtext("Boxploty jednotlivých premenných", outer = TRUE, cex = 1.4, font = 2)

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

Lineárna regresia

Model odhadujeme príkazom lm():

model <- lm(
  Life.expectancy ~ 1 + BMI + GDP + Schooling,
  data = udaje_2015
)

summary(model)

Objekt triedy lm nám poskytuje niekoľko výsledkov:

1. vektor odhadnutých koeficientov model$coefficients,
   
2. vektor rezíduí model$residuals,
   
3. vektor vyrovnaných hodnôt vysvetľovanej veličiny model$fitted.values,
   
4. maticu X – model.matrix(model).

coef(model)
head(residuals(model))
head(fitted(model))
head(model.matrix(model))

Súhrn odhadovaného modelu nám poskytuje súbor odhadovaný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.

Na základe Q–Q grafu získavame dojem o možných problémoch porušenia normality rezíduí.

# Rozloženie 2 x 2
par(mfrow = c(2, 2))

# Štandardné diagnostické grafy
plot(model)

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

Residuals vs. fitted – interpretácia

• Centrovanie okolo nuly
  Reziduály kolíšu približne okolo 0, čo je dobré. Naznačuje to, že model nemá výrazné systematické skreslenie v predikciách.

• Nelinearita
  Červená hladká čiara môže byť mierne zakrivená, čo naznačuje možnú miernu nelinearitu – modelu môže chýbať nelineárny tvar niektorej premennej.

• Rozptyl rezíduí (homoskedasticita)
  Vertikálny rozptyl sa javí ako približne konštantný v rámci rozsahu vyrovnaných hodnôt – to je dobrý dôkaz homoskedasticity.
  Ak by sa rezíduá rozširovali (tvorili “lievik”), naznačovalo by to heteroskedasticitu.

• Odľahlé hodnoty
  Niektoré body môžu ležať výrazne ďalej od ostatných – ide o potenciálne odľahlé hodnoty alebo vplyvné pozorovania.
  Môžeme ich preskúmať pomocou outlierTest(model) (z balíka car).

Q–Q plot – interpretácia

• Os X: teoretické kvantily – čo by sme očakávali, ak by rezíduá boli dokonale normálne rozložené.
  
• Os Y: štandardizované rezíduá – skutočné kvantily z našej vzorky.
  
• 45° prerušovaná čiara: ideálny prípad – ak sú rezíduá normálne rozložené, body ležia tesne na tejto čiare.

V našom prípade:

• Väčšina bodov leží blízko priamky – rezíduá sú približne normálne.
  
• V koncoch (extrémy) sú mierne odchýlky od priamky – môžu naznačovať odľahlé hodnoty alebo mierne ťažšie chvosty rozdelenia.
  
• Stredná časť grafu sa dobre zhoduje – väčšina rezíduí pekne zapadá do normálneho rozdelenia.

Scale–Location plot – interpretácia

• Os X: vyrovnané hodnoty (fitted values).
  
• Os Y: druhá odmocnina absolútnych štandardizovaných rezíduí.
  
• Červená čiara: LOESS (vyhladený trend).

V našom prípade:

• Body sú približne rovnomerne rozptýlené po osi X, bez typického lievika – naznačuje to konštantnú varianciu (homoskedasticitu).
  
• Červená čiara je relatívne rovná – pri zvyšovaní vyrovnaných hodnôt sa variancia zásadne nemení.
  
• Niekoľko bodov môže byť vyššie, ale nejde o extrémne hodnoty.

Residuals vs. leverage – interpretácia

• Os X: pákový efekt (leverage) – ako ďaleko je vektor vysvetľujúcich premenných daného bodu od stredu všetkých \(x\).
  
• Os Y: štandardizované rezíduá – ako ďaleko je pozorovanie od vyrovnanej hodnoty v jednotkách štandardnej odchýlky.
  
• Bodkované krivky: kontúry Cookovej vzdialenosti – mieru vplyvu pozorovania na odhadnuté koeficienty.

V našom prípade:

• Väčšina pozorovaní má nízky pákový efekt – typické pre väčšie vyvážené vzorky.
  
• Jednotlivé body s vyšším leverage sú potenciálne vplyvné, ale pokiaľ neprekračujú vysoké hodnoty Cookovej vzdialenosti, väčšinou neovplyvňujú výrazne regresné koeficienty.
  
• Žiadne pozorovanie zjavne neprekračuje vonkajšie krivky Cookovej vzdialenosti, preto sa nezdá, že by niektorý bod neprimerane ovplyvňoval model.

Test normality a odľahlých hodnôt – model

# Normality test (Jarque-Bera)
residuals_model <- residuals(model)
jb_test <- jarque.bera.test(residuals_model)
jb_test

# Outlier test (Bonferroni p-value)
outlier_test <- outlierTest(model)
outlier_test

Keďže sa normalita rezíduí nepreukázala ideálne, pokúsime sa eliminovať vplyv odľahlých hodnôt v prípade GDP.
Urobíme logaritmickú transformáciu tejto premennej a vylúčime BMI, ktoré sa ukázalo ako ťažšie interpretovateľné.

Nová regresia bude mať tvar:

\[ Life.expectancy = \beta_0 + \beta_1 \log(GDP) + \beta_2 Schooling + u \]

model2 <- lm(
  Life.expectancy ~ 1 + I(log(GDP)) + Schooling,
  data = udaje_2015
)

summary(model2)

# Rozloženie 2 x 2
par(mfrow = c(2, 2))

# Diagnostické grafy pre model2
plot(model2)

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

# Normality test pre model2
residuals_model2 <- residuals(model2)
jb_test2 <- jarque.bera.test(residuals_model2)
jb_test2

# Outlier test pre model2
outlier_test2 <- outlierTest(model2)
outlier_test2

Záver

Premenné GDP a Schooling štatisticky významne predlžujú strednú dĺžku života.
Na druhej strane premenná BMI (prevalencia nadváhy v populácii) v pôvodnom modeli poskytovala nejednoznačné a ťažko interpretovateľné výsledky, preto sme ju z finálneho modelu vylúčili.

Rezíduá modelov nevykazujú dokonalé normálne rozdelenie.
Keďže však máme relatívne veľký počet pozorovaní, porušenie normality nie je kritické a naďalej môžeme s týmito modelmi pracovať (spoliehame sa na asymptotické vlastnosti odhadov).

V modeli sa nepreukazujú výrazné nelinearity ani závažné problémy s heteroskedasticitou.
Transformácia GDP pomocou logaritmu prispela k lepšiemu správaniu rezíduí a interpretovateľnosti koeficientov.