Práca s regresiou

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("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 Rakúsko (Slovensko v databáze nemá zastúpenie):

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

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. :64.69	Min. :11.60	Min. :33.00	Min. :34.00	Min. :25.00	Min. :71.00	Min. :81.00	Min. :38843	Min. :8.010	Min. :78.10	Min. : 9.00
	1st Qu.:71.15	1st Qu.:12.07	1st Qu.:83.00	1st Qu.:46.75	1st Qu.:25.18	1st Qu.:83.00	1st Qu.:83.00	1st Qu.:40442	1st Qu.:8.158	1st Qu.:79.08	1st Qu.: 9.75
	Median :77.31	Median :12.12	Median :84.00	Median :61.50	Median :25.30	Median :84.50	Median :84.50	Median :42995	Median :8.310	Median :80.25	Median :11.55
	Mean :78.31	Mean :12.23	Mean :81.06	Mean :60.94	Mean :25.29	Mean :86.00	Mean :86.75	Mean :42376	Mean :8.296	Mean :79.97	Mean :10.83
	3rd Qu.:85.53	3rd Qu.:12.36	3rd Qu.:89.75	3rd Qu.:74.25	3rd Qu.:25.43	3rd Qu.:89.75	3rd Qu.:89.75	3rd Qu.:44258	3rd Qu.:8.400	3rd Qu.:80.92	3rd Qu.:11.80
	Max. :94.88	Max. :13.35	Max. :98.00	Max. :88.00	Max. :25.50	Max. :98.00	Max. :98.00	Max. :44550	Max. :8.640	Max. :81.50	Max. :12.10

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 Rakúsku 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 Rakúsku 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 Rakúsku (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)	-368.170	25.474	-14.453	0
Year	0.223	0.013	17.592	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.957	0.954	0.234	309.495	0	1	1.606	2.789	5.106	0.766	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 Rakúsku 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 Rakúsku (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)

Odhadnuté koeficienty regresie pre spotrebu alkoholu

term	estimate	std.error	statistic	p.value
(Intercept)	107.022	35.895	2.982	0.010
Year	-0.047	0.018	-2.641	0.019

# 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.333 

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

Adjusted R-squared: 0.285 

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 Rakúsku (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")