V tejto práci sa budeme zaoberať prvými krokmi v odhade regresnej funkcie. Budeme využívať databázu Life Expectancy (WHO) Fixed dataset, ktorá obsahuje ukazovatele na úrovni krajín, ako je priemerná dĺžka života, HDP, školské vzdelávanie a výdavky na zdravotníctvo.

Využívanie niektorých knižníc

rm(list=ls())
library(lmtest)   #  podpora regresie
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(outliers) # analyza odlahlych hodnot (outliers)
library(gptstudio)
library(kableExtra)
library(knitr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:kableExtra':
## 
##     group_rows
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(broom)
library(corrplot)
## corrplot 0.95 loaded

Príprava údajov - import z csv súboru

Súbor Life_Expectancy_Data obsahuje databázu determinantov očakávanej dĺžky života. Import údajov urobíme nasledovne

# import the dataset and create a data.frame udaje
#
udaje_svet <- read.csv("Life-Expectancy-Data-Updated.csv",header=TRUE,sep=",",dec=".",check.names = TRUE)
head(udaje_svet)
#

Databáza obsahuje údaje o 2938 pozorovaniach a 22 premenných. V tejto práci sa budeme zaoberať len časťou z nich, konkrétne tými, ktoré súvisia s dĺžkou dožitia. Na zažiatku si vyberieme krajinu, ktorej zdravotný stav chceme analyzovať. V tomto prípade ide o Portugalsko (Slovensko v databáze nemá zastúpenie):

# z databázy udaje_svet si vyberieme len tie pozorovania, ktoré sa týkajú Portugalska 
udaje <- subset(udaje_svet, Country == "Portugal")

Tabuľka uvedená nižšie nám poskuytuje základné popisné štatistiky vybraných kvantitatívnych premenných.

# niektoré štatistiky a ich prehľad v tabuľke KableExtra
library(kableExtra)
udaje %>%
  select(Adult_mortality,Alcohol_consumption, Hepatitis_B, Measles, BMI,Polio,Diphtheria,GDP_per_capita,Population_mln,Life_expectancy,Schooling) %>%
  summary() %>%
  kable() %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Adult_mortality Alcohol_consumption Hepatitis_B Measles BMI Polio Diphtheria GDP_per_capita Population_mln Life_expectancy Schooling
Min. : 73.28 Min. : 9.52 Min. :58.00 Min. :48.00 Min. :25.50 Min. :93.00 Min. :93.00 Min. :18585 Min. :10.29 Min. :76.30 Min. :6.800
1st Qu.: 83.07 1st Qu.:10.73 1st Qu.:94.00 1st Qu.:77.25 1st Qu.:25.70 1st Qu.:96.00 1st Qu.:96.75 1st Qu.:18830 1st Qu.:10.41 1st Qu.:77.58 1st Qu.:7.350
Median : 90.46 Median :11.52 Median :97.00 Median :95.00 Median :25.75 Median :96.50 Median :97.00 Median :19167 Median :10.49 Median :78.45 Median :7.750
Mean : 91.73 Mean :11.41 Mean :91.56 Mean :85.00 Mean :25.73 Mean :96.50 Mean :97.00 Mean :19199 Mean :10.47 Mean :78.74 Mean :7.831
3rd Qu.:100.28 3rd Qu.:12.23 3rd Qu.:97.25 3rd Qu.:95.25 3rd Qu.:25.80 3rd Qu.:97.25 3rd Qu.:98.00 3rd Qu.:19399 3rd Qu.:10.54 3rd Qu.:80.42 3rd Qu.:8.350
Max. :109.63 Max. :12.77 Max. :98.00 Max. :96.00 Max. :25.80 Max. :98.00 Max. :99.00 Max. :19986 Max. :10.57 Max. :81.10 Max. :9.100

Vyššie uvedená tabuľka nám poskytuje prehľad o základných štatistických charakteristikách vybraných premenných, ako sú priemerné hodnoty, rozptyl, minimum a maximum. Tieto informácie nám pomáhajú lepšie pochopiť rozdelenie a rozsah hodnôt v našich dátach. Na druhej strane je zaujímavá aj informácia o vzájomných vzťahoch medzi týmito premennými, čo môžeme merať pomocou korelačnej matice.

# grafický prehľad o korelačných vzťahoch vyjadruje nasledovný obrázok

cor_matrix <- cor(udaje %>% select(Adult_mortality,Alcohol_consumption
,Hepatitis_B, Measles, BMI,Polio,Diphtheria,GDP_per_capita,Population_mln,Life_expectancy,Schooling), use = "complete.obs")
corrplot(cor_matrix, method = "color", type = "upper", tl.col =
"black", tl.srt = 45, title = "Korelačná matica vybraných premenných", mar = c(0,0,1,0))

Uvedený graf nám poskytuje vizuálny prehľad o korelačných vzťahoch medzi vybranými premennými. Farby a intenzita farieb nám umožňujú rýchlo identifikovať silné pozitívne alebo negatívne korelácie. Upozorňujeme, že korelácia neznamená kauzalitu.

V nasledovnom je uvedený graf vývoja očakávanej dĺžky dožitia v Portugalsku v rokoch 2000-2015. Vidíme, že očakávaná dĺžka života sa zvyšovala, čo je pozitívny trend.

#  graf vývoja očakávanej dĺžky dožitia v Portugalsku v rokoch 2000-2015.
library(ggplot2)
ggplot(udaje, aes(x = Year, y = Life_expectancy)) +
  geom_line() +
  geom_point() +
  labs(title = "Vývoj očakávanej dĺžky dožitia v Portugalsku (2000-2015)",
       x = "Rok",
       y = "Očakávaná dĺžka dožitia") +
  theme_minimal()

Na začiatku sa pokúsme o vyrovnanie priebehu tejto premennej v čase pomocou lineárnej regresie, kde nezávislou premennou bude rok a závislou premennou bude očakávaná dĺžka dožitia. Odhadneme koeficienty tejto regresie a posúdime kvalitu vyrovnania pomocou ukazovateľov, ako je R-squared a p-value.

# vyrovnanie priebehu očakávanej dĺžky dožitia v čase
model <- lm(Life_expectancy ~ Year, data = udaje)

library(broom)
library(knitr)
library(kableExtra)

# koeficienty regresie
tidy(model) %>%
  kable(digits = 3, caption = "Odhadnuté koeficienty regresie") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Odhadnuté koeficienty regresie
term estimate std.error statistic p.value
(Intercept) -575.170 35.209 -16.336 0
Year 0.326 0.018 18.572 0 
# kvalita vyrovnania
glance(model) %>%
  kable(digits = 3, caption = "Ukazovatele kvality vyrovnania") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Ukazovatele kvality vyrovnania
r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance df.residual nobs
0.961 0.958 0.323 344.936 0 1 -3.573 13.145 15.463 1.464 14 16

Výsledky regresie nám ukazujú, že koeficient pre rok je pozitívny a štatisticky významný, čo naznačuje, že očakávaná dĺžka dožitia v Portugalsku sa zvyšovala v priebehu rokov 2000-2015. Jej priemerný ročný nárast dosahoval až takmer štvrť roka. Hodnota R-squared hodnota nám hovorí, že model vysvetľuje 95 % variability modelu. Podľa hodnoty p-value môžeme povedať, že model ako celok je štatisticky významný.

# teraz vyššie uvedený obrázok doplníme o regresnú priamku
ggplot(udaje, aes(x = Year, y = Life_expectancy)) +
  geom_line() +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(title = "Vývoj očakávanej dĺžky dožitia v Portugalsku (2000-2015) s regresnou priamkou",
       x = "Rok",
       y = "Očakávaná dĺžka dožitia") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

Vývoj spotreby alkoholu

Alkohol je predpokladaný negatívz faktor….

model_alkohol <-lm(Alcohol_consumption ~ Year, data = udaje)
tidy(model_alkohol) %>%
  kable(digits = 3, caption = "Odhadnuté koeficienty regresie pre spotrebu alkoholu") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
Odhadnuté koeficienty regresie pre spotrebu alkoholu
term estimate std.error statistic p.value
(Intercept) 416.781 38.734 10.760 0
Year -0.202 0.019 -10.465 0 
# and now the model quality statistics - just R sqared
r2 <- summary(model_alkohol)$r.squared
adj_r2 <- summary(model_alkohol)$adj.r.squared
# printujeme koeficient determinácie a upravený koeficient determinácie
cat("R-squared:", round(r2, 3), "\n")
## R-squared: 0.887
cat("Adjusted R-squared:", round(adj_r2, 3), "\n")
## Adjusted R-squared: 0.879

Obrázok o skutočných vyrovnaných hodnotách spotreby alkoholu

fitted_vals <- fitted(model_alkohol)
# vykreslenie skutočných a vyrovnaných hodnôt - Alcohol_consumption a fitted_values
ggplot(udaje, aes(x = Year)) +
  geom_line(aes(y = Alcohol_consumption), color = "blue", size = 1) +
  geom_line(aes(y = fitted_vals), color = "red", size = 1, linetype = "dashed") +
  theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

labs(title = "Skutočné vs. Vyrovnané hodnoty spotreby alkoholu v Portugalsku (2000-2015)",
       x = "Rok",
       y = "Spotreba alkoholu") +
  theme_minimal() +
  scale_y_continuous(limits = c(0, max(udaje$Alcohol_consumption, fitted_vals) * 1.1)) +
  theme(legend.position = "none")
## NULL

Diagnostické grafy

Residuals vs Fitted

Graf Residuals vs Fitted slúži na posúdenie linearity modelu a približnej konštantnosti rozptylu rezíduí. Reziduály by mali byť rozložené náhodne okolo nulovej hodnoty bez výrazného systematického vzoru. Ak sa objaví vizuálne významné zakrivenie červenej vyhladzovacej čiary,1 ktorá slúži na vizuálnu identifikáciu systematických odchýlok od náhodného rozloženia rezíduí, signalizuje to určité systematické problémy. Ak je LOESS približne horizontálna, model je špecifikovaný vhodne; ak sa zakrivuje, môže to naznačovať nelinearitu, chybnú špecifikáciu (teda výber vysvetľujúcich premenných) alebo iný problém v modeli. Ak sa body rozširujú alebo zužujú v tvare lievika, môže ísť o heteroskedasticitu.

plot(model, which = 1)
Residuals vs Fitted

Residuals vs Fitted

Reziduály sú rozložené približne okolo nulovej hodnoty, avšak pri najmenších, ako aj najväčších vyrovnaných hodnotách model systematicky nadhodnocuje vyrovnanú veličinu. Do určitej miery to môže ovplyvňovať odľahlé pozorovanie - a to jednak pozorovanie č. 1306, resp. 467. Naopak, pozorovanie 2782 z prostriedku vyrovanných hodnôt podhodnocuje významne vyrované hodnoty. Naše pozorovania sú ale málo početné a preto to zatiaľ špeciálne neriešime a podrobne sa k tomu môžeme vrítiť pri analýza ďalších grafov.

Q-Q plot

Q-Q graf porovnáva rozdelenie štandardizovaných rezíduí s teoretickým normálnym rozdelením. Ak body ležia približne na priamke, predpoklad normality rezíduí je približne splnený. Výraznejšie odchýlky na krajoch naznačujú možné problémy s extrémnymi hodnotami alebo s ťažšími koncami rozdelenia. Určitou pomôckou tu môže byť nasledovná tabuľka:

Interpretácia Q–Q grafu
Prípad Tvar.Q.Q.grafu Ľavý.chvost Pravý.chvost Interpretácia
Normálne rozdelenie Priamka Na priamke Na priamke Dáta majú normálne rozdelenie
Ľavostranná šikmosť (negatívna) Prehnutý nadol (∩ tvar) Nad priamkou Pod priamkou Dlhý ľavý chvost (extrémne nízke hodnoty)
Pravostranná šikmosť (pozitívna) Prehnutý nahor (∪ tvar) Pod priamkou Nad priamkou Dlhý pravý chvost (extrémne vysoké hodnoty)
Symetrické, vysoká špicatosť (> 3) S-tvar Pod priamkou Nad priamkou Ťažké chvosty (viac extrémnych hodnôt)
Symetrické, nízka špicatosť (< 3) Obrátený S-tvar Nad priamkou Pod priamkou Ľahké chvosty (menej extrémnych hodnôt)

Niekedy hovoríme o pomôcke:

  • Priamka → normálne rozdelenie
  • ∪ → pravostranná šikmosť
  • ∩ → ľavostranná šikmosť
  • S → vysoká špicatosť (ťažké chvosty)
  • obrátené S → nízka špicatosť
plot(model, which = 2)
Normal Q-Q plot

Normal Q-Q plot

Q-Q graf porovnáva empirické rozdelenie štandardizovaných rezíduí s teoretickým normálnym rozdelením. V prípade splnenia predpokladu normality by mali body ležať približne na diagonálnej priamke.

V našom prípade väčšina bodov v strednej časti grafu leží pomerne blízko priamky, čo naznačuje, že pre väčšinu pozorovaní je predpoklad normality približne splnený. To je priaznivý výsledok, keďže centrálna časť rozdelenia rezíduí sa správa podobne ako normálne rozdelenie.

Na druhej strane, na pravom konci grafu (pri vysokých kvantiloch) dochádza k výraznému odchýleniu bodov nad referenčnú priamku. To znamená, že pravý chvost rozdelenia je „ťažší“, než by sme očakávali pri normálnom rozdelení, teda v dátach sa nachádzajú extrémnejšie kladné reziduály. Niektoré konkrétne pozorovania (napr. označené body) sa odchyľujú výraznejšie, čo naznačuje prítomnosť odľahlých hodnôt.

Aj na ľavej strane sú odchýlky pod referenčnou čiarou, čo môže signalizovaťpravostranné zošikmenie rezíduí. Celkovo teda možno konštatovať, že normalita rezíduí nie je úplne splnená, najmä kvôli extrémnym hodnotám v pravom chvoste rozdelenia.

Z praktického hľadiska to znamená, že inferenčné štatistické testy založené na normalite môžu byť čiastočne ovplyvnené, avšak vzhľadom na relatívne dobré správanie v strednej časti rozdelenia nemusí ísť o zásadný problém, najmä pri väčších vzorkách.

Scale-Location plot

Scale-Location graf sa používa na posúdenie homoskedasticity, teda konštantnosti rozptylu rezíduí. Ak sú body rozložené približne rovnomerne a červená čiara je relatívne vodorovná, ide o priaznivý výsledok. Systematický rast alebo pokles naznačuje, že rozptyl rezíduí sa mení s úrovňou predikovaných hodnôt.

plot(model, which = 3)
Scale-Location plot

Scale-Location plot

V našom prípade sú body rozptýlené relatívne rovnomerne, avšak červená vyhladzovacia krivka má má pri vyšších vyrovnaných hodnotách vyrovnanej premennej výrazne rastúci trend. Môže to byť spojené s odľahlou hodnotou č. 1306.

Nevzniká typický „lievikovitý“ tvar, ktorý by signalizoval silné porušenie predpokladu konštantného rozptylu. Väčšina bodov sa nachádza v podobnom rozsahu hodnôt, čo naznačuje, že predpoklad homoskedasticity je porušený len mierne a to zrejme jednou odľahlou hodnotou.

Residuals vs Leverage

Graf Residuals vs Leverage pomáha identifikovať vplyvné pozorovania. Leverage to predstavuje vzialenosť daného pozorovania všetkých vysvetľujúcich premenných (teda vektora) od ostatných, resp. od stredného vektora všetkych pozorovaní. Títo vzdialenosť je normovaná na interval (0,1) a väčšina pozorovaní má veľmi malý leverage. Ak sa vyskytujú pozorovania s vysokým leverage, identifikujeme ich na vodorovnej osi a môžu znamenať potenciálne nebezpečenstvo v skreslení parametrov vyrovnávajúcej nadroviny (teda odhadovaných \(\beta\) koeficientov). Pozorovania s vysokou hodnotou leverage a súčasne veľkými rezíduami môžu môžu spôsobovať tento problém. Ak sa niektoré body nachádzajú blízko alebo za krivkami Cookovej vzdialenosti, je vhodné ich podrobnejšie preskúmať.

plot(model, which = 5)
Residuals vs Leverage

Residuals vs Leverage

V našom prípade väčšina pozorovaní má veľmi nízku hodnotu leverage a je sústredená blízko ľavej časti grafu, čo je typické a naznačuje, že väčšina dát reprezentujúca vysvetľujúce veličiny je dobre vyvážená. Tieto pozorovania zároveň nevykazujú extrémne rezíduá.

Kontúry Cookovej vzdialenosti naznačujú mieru vplyvu jednotlivých pozorovaní. Žiadny bod zjavne neprekračuje najvyššie hranice (napr. hodnotu 1), čo znamená, že model nie je dramaticky ovplyvnený jedným extrémnym pozorovaním. Odľahlé pozorovanie č.1306 sa nachádza v hraničnej oblasti s Coockovou vzdialenosťou medzi 0.5 a 1.

Celkovo možno konštatovať, že model neobsahuje výrazne problematické vplyvné pozorovania, avšak existuje jeden potenciálne rizikový bod, ktorému je vhodné venovať dodatočnú pozornosť. Hodnoty Leverage ako aj Cookove vzdialenosti súvisia úzko s tzv Hat maticou.

  1. LOESS (LOcal regrESSion) krivka predstavuje lokálne vyhladený trend medzi premennými bez predpokladu konkrétneho funkčného tvaru.↩︎