Práca s regresiou 7.tyzden
Lucia Faltýnková (s pomocou vybraných zdrojov)
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
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("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)| 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)| 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 |
NAVý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. Jejnárast dosiahol až takmer štyri roky. 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()
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)| 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") +
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")
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)
Test na normalitu a test na autokoreláciu
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.
# test na normalitu rezíduí pomocou Shapiro-Wilk testu
shapiro_test <- shapiro.test(residuals(model))
cat("Normalita: Shapiro-Wilk test p-value:", shapiro_test$p.value, "\n")Normalita: Shapiro-Wilk test p-value: 0.1405092
# test na normalitu rezíduí pomocou Jarque-Berra testu
library(tseries)
jarque_bera_test <- jarque.bera.test(residuals(model))
cat("Normalita: Jarque-Bera test p-value:", jarque_bera_test$p.value, "\n")Normalita: Jarque-Bera test p-value: 0.8529862 # test na autokoreláciu rezíduí pomocou Breusch-Godfrey
bg_test <- bgtest(model, order = 1)Error in bgtest(model, order = 1) : could not find function "bgtest"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:
| 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)
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)
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.
# test na heteroskedasticitu pomocou Breusch-Pagan testu
bp_test <- bptest(model)
cat("Heteroskedasticita: Breusch-Pagan test p-value:",bp_test$p.value, "\n")Heteroskedasticita: Breusch-Pagan test p-value: 0.1987857 Residuals vs Leverage
Graf Residuals vs Leverage pomáha identifikovať vplyvné pozorovania. Leverage to predstavuje vzdialenosť 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)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.
Test na odľahlé hodnoty (z knižnice “car”)
LOESS (LOcal regrESSion) krivka predstavuje lokálne vyhladený trend medzi premennými bez predpokladu konkrétneho funkčného tvaru.↩︎