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")

```