Úvod

V tomto cvičení pracujem s údajmi z databázy
WHO Life Expectancy Data.

Dataset som si stiahla z Kaggle a následne som ho nahrala do svojho projektu
v Posit Cloud, kde sa nachádza v priečinku:

udaje/Life Expectancy Data.csv

Údaje obsahujú informácie o očakávanej dĺžke života a o rôznych
socio-ekonomických a zdravotných ukazovateľoch pre jednotlivé krajiny.

Analyzujem rok 2015 a premenné:

Moja hypotéza je, že vyšší GDP a viac Schooling sú spojené
s dlhšou dĺžkou života, zatiaľ čo BMI môže mať nejednoznačný vplyv.

Balíčky

library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(lmtest)
library(sandwich)
library(car)
## Loading required package: carData
rm(list = ls())

Načítanie a úprava údajov

udaje <- read.csv("udaje/Life Expectancy Data.csv",
                  dec = ".", sep = ",", header = TRUE)

names(udaje) <- make.names(names(udaje))

udaje.2015 <- udaje[udaje$Year == 2015,
                    c("Life.expectancy", "BMI", "GDP", "Schooling")]

udaje.2015$GDP[udaje.2015$GDP <= 0] <- NA

column_medians <- sapply(udaje.2015, median, na.rm = TRUE)

udaje_imputed <- udaje.2015
for (col in names(udaje_imputed)) {
  udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]
}

udaje.2015 <- udaje_imputed
summary(udaje.2015)
##  Life.expectancy      BMI             GDP             Schooling    
##  Min.   :51.00   Min.   : 2.50   Min.   :   33.68   Min.   : 4.90  
##  1st Qu.:65.75   1st Qu.:24.30   1st Qu.: 1152.31   1st Qu.:11.10  
##  Median :73.90   Median :48.60   Median : 2916.23   Median :13.10  
##  Mean   :71.62   Mean   :42.82   Mean   : 6508.81   Mean   :12.94  
##  3rd Qu.:76.95   3rd Qu.:61.35   3rd Qu.: 5881.99   3rd Qu.:14.85  
##  Max.   :88.00   Max.   :77.60   Max.   :66346.52   Max.   :20.40

Boxploty premenných

par(mfrow = c(2, 2))
par(mar = c(4, 4, 2, 1))

for (col in names(udaje.2015)) {
  boxplot(udaje.2015[[col]],
          main = col,
          col  = "lightblue",
          xlab = "Hodnota")
}

mtext("Boxploty jednotlivých premenných (2015)",
      outer = TRUE, cex = 1.2, font = 2)

par(mfrow = c(1, 1))

Boxploty ukazujú, že najviac odľahlých hodnôt má premenná GDP,
čo môže ovplyvniť vlastnosti modelu pri klasickej regresii.

Lineárny model 1

Life.expectancy ~ BMI + GDP + Schooling

model <- lm(Life.expectancy ~ BMI + GDP + Schooling,
            data = udaje.2015)

summary(model)
## 
## Call:
## lm(formula = Life.expectancy ~ BMI + GDP + Schooling, data = udaje.2015)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.6471  -2.4024   0.4563   3.2355  12.8591 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.419e+01  1.839e+00  24.033   <2e-16 ***
## BMI         4.948e-02  2.144e-02   2.308   0.0222 *  
## GDP         6.972e-05  3.845e-05   1.813   0.0715 .  
## Schooling   1.921e+00  1.645e-01  11.676   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.033 on 179 degrees of freedom
## Multiple R-squared:  0.6224, Adjusted R-squared:  0.6161 
## F-statistic: 98.36 on 3 and 179 DF,  p-value: < 2.2e-16

V prvom modeli vychádza premenná Schooling ako významná,
GDP má pozitívny vplyv a BMI má málo jasný efekt.

Diagnostické grafy modelu 1

par(mfrow = c(2, 2))
plot(model)

par(mfrow = c(1, 1))

Diagnostické grafy naznačujú miernu nenormalitu rezíduí
a možné odľahlé hodnoty najmä pri vysokom HDP.

Test normality a outlierov (model 1)

residuals_m1 <- residuals(model)
jarque.bera.test(residuals_m1)
## 
##  Jarque Bera Test
## 
## data:  residuals_m1
## X-squared = 21.391, df = 2, p-value = 2.265e-05
outlierTest(model)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##      rstudent unadjusted p-value Bonferroni p
## 433 -3.642954         0.00035336     0.064664

Jarque–Bera test zvyčajne zamieta normalitu rezíduí.

Lineárny model 2

Life.expectancy ~ log(GDP) + Schooling

Model upravujem:
- BMI vynechávam,
- GDP transformujem logaritmom, aby som znížila vplyv extrémov.

model2 <- lm(Life.expectancy ~ I(log(GDP)) + Schooling,
             data = udaje.2015)

summary(model2)
## 
## Call:
## lm(formula = Life.expectancy ~ I(log(GDP)) + Schooling, data = udaje.2015)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.7351  -2.4761   0.5209   3.4551  10.6277 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  40.0434     2.1992  18.208   <2e-16 ***
## I(log(GDP))   0.6338     0.3002   2.112   0.0361 *  
## Schooling     2.0561     0.1557  13.203   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.09 on 180 degrees of freedom
## Multiple R-squared:  0.6118, Adjusted R-squared:  0.6075 
## F-statistic: 141.8 on 2 and 180 DF,  p-value: < 2.2e-16

Premenné log(GDP) aj Schooling sú štatisticky významné
a smerujú k lepšiemu správaniu modelu ako v prvej verzii.

Diagnostické grafy modelu 2

par(mfrow = c(2, 2))
plot(model2)

par(mfrow = c(1, 1))

Reziduá modelu 2 sú lepšie rozložené
a grafy ukazujú menšie problémy s odľahlými hodnotami.

Test normality a outlierov (model 2)

residuals_m2 <- residuals(model2)
jarque.bera.test(residuals_m2)
## 
##  Jarque Bera Test
## 
## data:  residuals_m2
## X-squared = 30.107, df = 2, p-value = 2.899e-07
outlierTest(model2)
##      rstudent unadjusted p-value Bonferroni p
## 433 -3.828567         0.00017805     0.032583

Log-transformácia HDP znižuje extrémne odchýlky v reziduách.

Záver

Na základe analýzy môžem povedať, že:

Model sa javí ako vhodný na opis vzťahu medzi dĺžkou života
a socio-ekonomickými faktormi krajín v roku 2015.

Ďakujem za dnešné zadanie, teším sa na ďalšie. Vanessa