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("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 Indiu (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 == "India")

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. :177.9 Min. :0.930 Min. :29.00 Min. :27.00 Min. :20.70 Min. :57.00 Min. :58.0 Min. : 758.0 Min. :1057 Min. :62.50 Min. :4.400
1st Qu.:187.6 1st Qu.:1.177 1st Qu.:56.00 1st Qu.:34.00 1st Qu.:20.98 1st Qu.:58.00 1st Qu.:62.5 1st Qu.: 884.8 1st Qu.:1125 1st Qu.:64.00 1st Qu.:4.700
Median :199.5 Median :1.745 Median :67.00 Median :34.00 Median :21.20 Median :68.00 Median :67.5 Median :1084.5 Median :1192 Median :65.60 Median :5.100
Mean :199.5 Mean :1.978 Mean :61.56 Mean :38.81 Mean :21.20 Mean :69.62 Mean :71.0 Mean :1111.8 Mean :1189 Mean :65.59 Mean :5.169
3rd Qu.:211.0 3rd Qu.:2.962 3rd Qu.:70.75 3rd Qu.:36.50 3rd Qu.:21.43 3rd Qu.:79.00 3rd Qu.:82.0 3rd Qu.:1306.5 3rd Qu.:1254 3rd Qu.:67.20 3rd Qu.:5.450
Max. :221.6 Max. :3.040 Max. :87.00 Max. :69.00 Max. :21.70 Max. :86.00 Max. :87.0 Max. :1606.0 Max. :1310 Max. :68.60 Max. :6.300

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 Indii 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 Indii (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) -769.296 5.184 -148.412 0
Year 0.416 0.003 161.066 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.999 0.999 0.048 25942.12 0 1 27.08 -48.161 -45.843 0.032 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 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 Indii (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) -349.512 28.026 -12.471 0
Year 0.175 0.014 12.542 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.918
cat("Adjusted R-squared:", round(adj_r2, 3), "\n")
## Adjusted R-squared: 0.912

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 Indii (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")
## 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.