Práca s regresiou

Úvod do regresie

Práca s údajmi z Japonska za roky 2000-2015

V tejto práci sa budeme zaoberať prvými krokmi v odhade regresnej funkcie. Budeme využívať databázu [Life Expectancy (WHO) Fixed dataset] (https://www.kaggle.com/datasets/lashagoch/life-expectancy-who-updated), ktorá obsahuje ukazovatele na úrovni krajín, ako je priemerná dĺžka života, HDP, školské vzdelávanie a výdavky na zdravotníctvo. Vyberáme krajinu Japonsko.

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

rm(list=ls())
library(lmtest)   #  podpora regresie
library(outliers) # analyza odlahlych hodnot (outliers)
library(gptstudio)
library(kableExtra)
library(knitr)
library(dplyr)
library(broom)
library(corrplot)

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("udaje/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 už spomínané Japonsko:

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

Tabuľka uvedená nižšie nám poskytuje základné popisné štatistiky vybraných kvantitatívnych premenných. Vylúčená je Hepatitída B nakoľko údaj bol konštantný.

# niektoré štatistiky a ich prehľad v tabuľke KableExtra
library(kableExtra)
udaje %>%
  select(Adult_mortality,Alcohol_consumption,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	Measles	BMI	Polio	Diphtheria	GDP_per_capita	Population_mln	Life_expectancy	Schooling
	Min. :53.59	Min. :6.980	Min. :78.00	Min. :22.50	Min. :81.00	Min. :85.00	Min. :31416	Min. :126.8	Min. :81.10	Min. :10.70
	1st Qu.:61.20	1st Qu.:7.638	1st Qu.:86.00	1st Qu.:22.57	1st Qu.:95.75	1st Qu.:96.00	1st Qu.:31785	1st Qu.:127.4	1st Qu.:81.88	1st Qu.:11.07
	Median :64.70	Median :7.970	Median :88.00	Median :22.60	Median :97.50	Median :97.00	Median :33055	Median :127.7	Median :82.55	Median :11.35
	Mean :64.52	Mean :7.932	Mean :88.81	Mean :22.62	Mean :96.00	Mean :96.25	Mean :32972	Mean :127.6	Mean :82.46	Mean :11.47
	3rd Qu.:69.23	3rd Qu.:8.402	3rd Qu.:92.25	3rd Qu.:22.70	3rd Qu.:99.00	3rd Qu.:98.00	3rd Qu.:33666	3rd Qu.:127.9	3rd Qu.:82.95	3rd Qu.:11.85
	Max. :72.95	Max. :8.650	Max. :97.00	Max. :22.70	Max. :99.00	Max. :99.00	Max. :34961	Max. :128.1	Max. :83.80	Max. :12.50

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,Measles,BMI,Polio,Diphtheria,GDP_per_capita,Population_mln,Life_expectancy,Schooling), use = "complete.obs")
corrplot(cor_matrix, method = "color", type = "upper", tl.col =
"darkblue", 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 zobrazení je uvedený graf vývoja očakávanej dĺžky dožitia v Japonsku v rokoch 2000-2015. Vidíme, že očakávaná dĺžka života sa zvyšovala od roku 2000 z 81 na 83 rokov v roku 2015, čo je pozitívny trend.

#  graf vývoja očakávanej dĺžky dožitia v Japonsku 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 Japonsku (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)	-242.582	17.119	-14.170	0
Year	0.162	0.009	18.987	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.963	0.96	0.157	360.514	0	1	7.965	-9.93	-7.613	0.346	14	16

NA

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 Japonsku sa zvyšovala v priebehu rokov 2000-2015. Jej priemerný ročný nárast dosahoval 0,162. Hodnota R-squared hodnota nám hovorí, že model vysvetľuje 96 % 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 = "green") +
  labs(title = "Vývoj očakávanej dĺžky dožitia v Japonsku (2000-2015) s regresnou priamkou",
       x = "Rok",
       y = "Očakávaná dĺžka dožitia") +
  theme_minimal()

###Vývoj spotreby alkoholu Teraz sa pozrieme, ako sa vyvýja konzumácia alkoholu v čase od roku 2000 do roku 2015 v Japonsku. Nakoľko alkohol je faktorom negatívne vplývajúcim na dĺžku dožitia.

model_alkohol <- lm(Alcohol_consumption ~ Year, data = udaje)
# kvalita regresie
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)	-91.115	53.198	-1.713	0.109
Year	0.049	0.026	1.862	0.084

# a teraz model quality statistics R squared
r2<- summary(model_alkohol)$r.square
adj_r2 <-summary(model_alkohol)$adj.r.square
#printujeme koeficient determinancie a upravený koeficient determinancie
cat ("R-squared:",round (r2,3), "\n")

R-squared: 0.198 

cat ("Adjusted R-squared:",round (adj_r2,3), "\n")

Adjusted R-squared: 0.141 

Výsledky regresie ukazujú, že koeficient pre premennú rok je pozitívny, čo naznačuje mierny rast spotreby alkoholu v Japonsku v období rokov 2000 – 2015. Odhadovaný koeficient 0,049 znamená, že spotreba alkoholu sa v priemere zvyšovala približne o 0,049 jednotky ročne. Hodnota p-value (0,084) však naznačuje, že tento koeficient nie je štatisticky významný na hladine významnosti 5 %, preto nie je možné s dostatočnou istotou potvrdiť existenciu štatisticky významného trendu rastu spotreby alkoholu v sledovanom období.

#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 = "darkblue", size = 1) +
  geom_line(aes(y = fitted_vals), color = "orange", size = 1, linetype = "dashed") +
  labs(title = "Skutočné vs. Vyrovnané hodnoty spotreby alkoholu v Japonsku (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")

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,¹ 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)

Reziduá sú rozložené približne okolo nulovej hodnoty, avšak pri najmenších aj najväčších odhadnutých hodnotách model systematicky nadhodnocuje, resp. podhodnocuje skutočné hodnoty. Do určitej miery to môže byť ovplyvnené odľahlými pozorovaniami, najmä pozorovaním č. 1411 a 530. Počet pozorovaní je však relatívne malý, preto túto skutočnosť zatiaľ nepovažujeme za zásadný problém a podrobnejšie sa k nej môžeme vrátiť pri analýze ďalších diagnostický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)

Upravený popis podľa grafu:

---

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ý. Centrálna časť rozdelenia rezíduí sa teda správa podobne ako normálne rozdelenie.

Na ľavom konci grafu dochádza k výraznejšiemu odchýleniu bodov pod referenčnú priamku. Najmä pozorovanie č. 530 predstavuje extrémne záporné reziduum, čo naznačuje ťažší ľavý chvost rozdelenia oproti normálnemu rozdeleniu. Menšia odchýlka je viditeľná aj pri pozorovaní č. 566.

Na pravom konci grafu sú body relatívne blízko referenčnej priamky. Pozorovanie č. 1411 sa nachádza mierne pod priamkou, čo poukazuje len na malú odchýlku od normality a nepredstavuje výrazný problém.

Celkovo možno konštatovať, že normalita rezíduí je približne splnená, avšak s miernymi odchýlkami v chvostoch rozdelenia, najmä v ľavom chvoste spôsobenými odľahlými pozorovaniami. Z praktického hľadiska ide o mierne porušenie predpokladu normality, ktoré však vzhľadom na dobré správanie väčšiny rezíduí nemusí predstavovať zásadný problém pre ďalšiu analýzu.

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)

Tu je upravený popis presnejšie podľa grafu:

Body sú rozptýlené pomerne rovnomerne v celom rozsahu vyrovnaných hodnôt, bez výrazného systematického vzoru. Červená vyhladzovacia krivka mierne rastie pri nižších hodnotách, následne sa stabilizuje a ku koncu len veľmi mierne kolíše, bez výrazného rastúceho trendu pri vyšších hodnotách.

Niektoré body (napr. pozorovania označené 530, 566 a 1411) sa od ostatných mierne odchyľujú, môžu teda predstavovať potenciálne odľahlé alebo vplyvné pozorovania, avšak ich vplyv nie je extrémny.

Nevzniká typický „lievikovitý“ tvar, ktorý by naznačoval výrazné porušenie predpokladu konštantného rozptylu. Variabilita rezíduí sa javí ako relatívne stabilná naprieč hodnotami, čo naznačuje, že predpoklad homoskedasticity je splnený, prípadne len veľmi mierne narušený.

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)

Tu je upravený popis podľa zobrazeného grafu:

V našom prípade má väčšina pozorovaní nízke až stredné hodnoty leverage a je rozptýlená najmä v intervale približne od 0.05 do 0.20, čo naznačuje, že žiadne pozorovanie nemá extrémne vysokú pákovú hodnotu. Dáta teda nepôsobia výrazne nevyvážene z hľadiska vysvetľujúcich premenných.

Z hľadiska rezíduí väčšina bodov nevykazuje extrémne hodnoty, avšak niektoré pozorovania (najmä 530) majú výraznejšie záporné rezíduum. Okrem toho body ako 1197 a 566 majú vyššie leverage, čo znamená, že môžu mať väčší vplyv na odhad modelu.

Kontúry Cookovej vzdialenosti ukazujú mieru vplyvu jednotlivých pozorovaní. Žiadny bod zjavne neprekračuje hranicu hodnoty 1, avšak pozorovanie 530 sa nachádza v oblasti vyšších hodnôt Cookovej vzdialenosti (medzi 0.5 a 1), čo naznačuje, že môže byť potenciálne vplyvné a zaslúži si bližšie preskúmanie.

Celkovo možno konštatovať, že model neobsahuje extrémne vplyvné pozorovania, ktoré by zásadne narúšali jeho stabilitu, avšak existujú jednotlivé body (najmä 530, prípadne 1197 a 566), ktoré môžu mať zvýšený vplyv na výsledky a je vhodné im venovať pozornosť.

---

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