Úvod do regresie

Práca s údajmi z Albánska 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 Albánsko.

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)

 udaje_svet <- udaje_svet[-992,]

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é Albánsko:

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

Tabuľka uvedená nižšie nám poskytuje 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. :75.20	Min. :3.960	Min. :96.00	Min. :90.00	Min. :25.20	Min. :97.0	Min. :97.00	Min. :2144	Min. :2.880	Min. :74.30	Min. :8.700
1st Qu.:80.06	1st Qu.:4.330	1st Qu.:98.00	1st Qu.:94.50	1st Qu.:25.55	1st Qu.:97.5	1st Qu.:97.50	1st Qu.:2599	1st Qu.:2.905	1st Qu.:75.10	1st Qu.:9.100
Median :82.11	Median :4.510	Median :98.00	Median :98.00	Median :25.90	Median :98.0	Median :98.00	Median :3298	Median :2.950	Median :75.90	Median :9.200
Mean :82.36	Mean :4.749	Mean :98.13	Mean :96.33	Mean :25.90	Mean :98.2	Mean :98.13	Mean :3145	Mean :2.961	Mean :76.08	Mean :9.273
3rd Qu.:84.94	3rd Qu.:5.100	3rd Qu.:99.00	3rd Qu.:98.00	3rd Qu.:26.25	3rd Qu.:99.0	3rd Qu.:99.00	3rd Qu.:3707	3rd Qu.:3.020	3rd Qu.:77.10	3rd Qu.:9.450
Max. :90.94	Max. :5.790	Max. :99.00	Max. :99.00	Max. :26.60	Max. :99.0	Max. :99.00	Max. :3953	Max. :3.060	Max. :78.00	Max. :9.700

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 =
"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 Albánsku v rokoch 2000-2015. Vidíme, že očakávaná dĺžka života sa zvyšovala od roku 2000 z 74 na 78 rokov v roku 2015, čo je pozitívny trend.

#  graf vývoja očakávanej dĺžky dožitia v Albánsku 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 Albánsku (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+Alcohol_consumption+Adult_mortality+Incidents_HIV,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)	-578.820	23.965	-24.152	0.000
Year	0.324	0.011	28.549	0.000
Alcohol_consumption	-0.215	0.024	-8.829	0.000
Adult_mortality	0.062	0.015	4.143	0.002
Incidents_HIV	4.009	3.985	1.006	0.338

# 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.035	4416.244	0	4	32.249	-52.499	-48.251	0.012	10	15

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 Albánsku sa zvyšovala v priebehu rokov 2000-2015. Jej priemerný ročný nárast dosahoval 0,085. Hodnota R-squared hodnota nám hovorí, že model vysvetľuje 99 % variability modelu. Podľa hodnoty p-value môžeme povedať, že model ako celok je štatisticky významný. Hodnoty BMI a konzumácia alkoholu vykazujú záporné hodnoty.

# 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 Albánsku (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 Teraz sa pozrieme, ako sa v Albánsku vyvýjala konzumácia alkoholu v čase od roku 2000 do roku 2015. Nakoľko alkohol je jedným z faktorov, ktpré sme zaradili do modelu, negatívne vplývajúcich 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)	9.984	70.238	0.142	0.889
Year	-0.003	0.035	-0.075	0.942

# 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

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

## Adjusted R-squared: -0.076

Výsledky regresie ukazujú, že koeficient pre premennú rok je pozitívny, čo naznačuje veľmi mierny rast spotreby alkoholu v Albánsku v období rokov 2000 – 2015. Odhadovaný koeficient 0,016 znamená, že spotreba alkoholu sa v priemere zvyšovala približne o 0,016 jednotky ročne. Hodnota p-value (0,0627) však naznačuje, že tento koeficient nie je štatisticky významný na hladine významnosti 5 %, preto nie je možné 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 = "blue", size = 1) +
  geom_line(aes(y = fitted_vals), color = "magenta", size = 1, linetype = "dashed") +
  labs(title = "Skutočné vs. Vyrovnané hodnoty spotreby alkoholu v Albánsku (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.

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,¹ 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)

Residuals vs Fitted

Test na normalitu a test na autokoreláciu

#Test na normalitu reziduí 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.2980539

#Test na normalitu rezuduí pomocou Jarque-Berra testu
library(tseries)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

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.5488288

Vyššie uvedené testy nám poskytujú štatistické hodnoty p-value, ktoré nám pomáhajú posúdiť, či reziduály modelu spĺňajú predpoklad normality. V našom prípade sú p-hodnoty (Shapiro-Wilkov test: p = 0.3405562; Jarque-Bera test: p = 0.5548655) nižšie než zvolená hladina významnosti (napr. 0,05), čo vedie k zamietnutiu nulovej hypotézy o normalite rezíduí. Napriek tomu, vzhľadom na to, že ide o relatívne malé odchýlky od hranice významnosti, budeme tento predpoklad v ďalšej analýze prehliadať a považovať ho za dostatočne splnený.

#Test na autokoreláciu reziduí pomocou Breusch-Godfrey
bg_test <- bgtest(model, order = 1)
cat("Autokorelace: Breusch-Godfrey test p-value:", bg_test$p.value, "\n")

## Autokorelace: Breusch-Godfrey test p-value: 0.3019012

Test na autokoreláciu nám poskytuje p-value, ktoré nám pomáhajú posúdiť, či reziduály modelu vykazujú autokoreláciu. Keďže je p-value väčia než zvolená hladina významnosti (0,1), nemôžeme prijať alternatívnu hypotézu o prítomnosti autokorelacie rezíduí.

Reziduá sú rozložené približne okolo nulovej hodnoty, avšak pri najmenších aj najväčších odhadnutých hodnotách model systematicky nadhodnocuje, resp. podhodnocuje skutočné hodnoty. Do určitej miery to môže byť ovplyvnené odľahlými pozorovaniami, najmä pozorovaním č. 114, 1003 a 2164. Počet pozorovaní je však relatívne malý, preto túto skutočnosť zatiaľ nepovažujeme za zásadný problém a podrobnejšie sa k nej môžeme vrátiť pri analýze ďalších diagnostických grafov.

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:

Interpretácia Q–Q grafu
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)

Normal Q-Q plot

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ý. Centrálna časť rozdelenia rezíduí sa teda správa podobne ako normálne rozdelenie.

Na pravom konci grafu dochádza k miernemu odchýleniu bodov pod referenčnú priamku. Najmä pozorovanie č. 1003 predstavuje záporné reziduum, čo naznačuje ťažší pravý chvost rozdelenia oproti normálnemu rozdeleniu. Menšia odchýlka je viditeľná aj pri pozorovaní č. 2164kroeé leží takmer na priamke.

Na ľavom konci grafu je bod relatívne blízko referenčnej priamky. Pozorovanie č. 1162 sa nachádza mierne nad priamkou, čo poukazuje len na malú odchýlku od normality a nepredstavuje výrazný problém.

Celkovo možno konštatovať, že normalita rezíduí je približne splnená, avšak s miernymi odchýlkami v chvostoch rozdelenia. Z praktického hľadiska ide o mierne porušenie predpokladu normality, ktoré však vzhľadom na dobré správanie väčšiny rezíduí nemusí predstavovať zásadný problém pre ďalšiu analýzu.

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)

Scale-Location plot

Body sú rozptýlené pomerne rovnomerne v celom rozsahu vyrovnaných hodnôt, bez výrazného systematického vzoru. Červená vyhladzovacia krivka mierne klesá pri nižších hodnotách, následne sa stabilizuje a ku koncu len veľmi mierne kolíše, bez výrazného rastúceho trendu pri vyšších hodnotách.

Niektoré body (napr. pozorovania označené 1162, 2164 a 1003) sa od ostatných mierne odchyľujú, môžu teda predstavovať potenciálne odľahlé alebo vplyvné pozorovania, avšak ich vplyv nie je extrémny.

Nevzniká typický „lievikovitý“ tvar, ktorý by naznačoval výrazné porušenie predpokladu konštantného rozptylu. Variabilita rezíduí sa javí ako relatívne stabilná naprieč hodnotami, čo naznačuje, že predpoklad homoskedasticity je splnený, prípadne len veľmi mierne narušený.

Heteroskedaticita sa nepotvrdila, je porušená.

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

Keďže Breusch-Pagan test nám poskytuje p-value, ktoré je väčšie než zvolená hladina významnosti (0,1), nemôžeme zamietnuť nulovú hypotézu a konštatovať, že rozptyl rezíduí nie je konštantný.

Residuals vs Leverage

Graf Residuals vs Leverage pomáha identifikovať vplyvné pozorovania. Leverage to predstavuje vzialenosť 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)

Residuals vs Leverage

Väčšina pozorovaní má nízke až stredné hodnoty leverage (približne do 0,4) a je rozptýlená okolo nulových štandardizovaných rezíduí, prevažne v intervale od −1 do 1. To naznačuje, že prevažná časť dát je modelom zachytená pomerne dobre a nevykazuje výrazné odchýlky.

Z grafu však vidno niekoľko výnimiek. Pozorovanie 2164 má relatívne vysoké kladné rezíduum (nad 2) pri strednej hodnote leverage, čo môže naznačovať potenciálny problém z hľadiska prispôsobenia modelu. Naopak, pozorovanie 1162 má veľmi vysokú hodnotu leverage (okolo 0,8) a zároveň výrazne záporné rezíduum, čo z neho robí bod s potenciálne silným vplyvom na model.

Ďalšie body, ako napríklad 581, majú mierne zvýšené rezíduá aj leverage, no ich vplyv je menej výrazný.

Kontúry Cookovej vzdialenosti ukazujú, že žiadne pozorovanie neprekračuje hranicu hodnoty 1, takže sa v dátach nenachádzajú extrémne vplyvné body. Napriek tomu sa pozorovania 1162 a čiastočne aj 2164 nachádzajú v oblastiach vyšších hodnôt Cookovej vzdialenosti, čo naznačuje, že môžu mať citeľný vplyv na odhady modelu.

Celkovo je model relatívne stabilný, avšak niektoré jednotlivé pozorovania (najmä 1162 a 2164) by mali byť bližšie preskúmané, keďže môžu ovplyvňovať výsledky regresie.

Test na odľahlé hodnoty (z knižnice car)

## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##      rstudent unadjusted p-value Bonferroni p
## 2164 2.740913           0.022813       0.3422

Test na odľahlé hodnoty nám poskytuje informácie o pozorovaniach, ktoré môžu byť potenciálne odľahlé. V tomto prípade test identifikoval pozorovanie č. 2164 ako odľahlé (rstudent = 2,74) s neadjustovanou p-hodnotou 0,022813 a p-hodnotou po Bonferroniho korekcii 0,3422. Keďže po Bonferroniho korekcii je p-hodnota väčšia než zvolená hladina významnosti (napr. 0,05),presahuje hladinu významnosti. To znamená, že pri zohľadnení viacnásobného testovania nemožno toto pozorovanie považovať za štatisticky významnú odľahlú hodnotu.

Napriek tomu ide o pozorovanie s relatívne vysokou hodnotou studentizovaného rezídua, a preto ho možno označiť za potenciálne problematické. V kombinácii s diagnostickým grafom (Residuals vs Leverage), kde sa toto pozorovanie nachádza v oblasti vyšších rezíduí, je vhodné venovať mu bližšiu pozornosť a zvážiť jeho vplyv na model.

library(lmtest)
library(sandwich)

Breusch–Pagan

bptest(model)

## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 2.6119, df = 4, p-value = 0.6247

Breusch–Godfrey

bgtest(model)

## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model
## LM test = 1.0658, df = 1, p-value = 0.3019

Robustné štandardné chyby

White robustné štandardné chyby

coeftest(model, vcov = vcovHC(model, type = "HC1"))

## 
## t test of coefficients:
## 
##                        Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept)         -5.7882e+02  1.8161e+01 -31.8716 2.175e-11 ***
## Year                 3.2408e-01  8.6497e-03  37.4669 4.370e-12 ***
## Alcohol_consumption -2.1476e-01  3.0672e-02  -7.0019 3.707e-05 ***
## Adult_mortality      6.1531e-02  1.0455e-02   5.8853 0.0001541 ***
## Incidents_HIV        4.0095e+00  3.1092e+00   1.2896 0.2262278    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Newey-West robustné štandardné chyby

coeftest(model, vcov = NeweyWest(model))

## 
## t test of coefficients:
## 
##                        Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept)         -5.7882e+02  9.0008e+00 -64.3074 2.012e-14 ***
## Year                 3.2408e-01  4.3150e-03  75.1050 4.275e-15 ***
## Alcohol_consumption -2.1476e-01  2.3141e-02  -9.2803 3.138e-06 ***
## Adult_mortality      6.1531e-02  5.0098e-03  12.2821 2.347e-07 ***
## Incidents_HIV        4.0095e+00  1.9074e+00   2.1021   0.06186 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

vif(model)

##                Year Alcohol_consumption     Adult_mortality       Incidents_HIV 
##           30.275680            2.212030           44.615286            7.105055

Chyba špecifikácie

Jednoduchý príklad

model_1 <- lm(Life_expectancy ~ GDP_per_capita, data = udaje)
summary(model_1)

## 
## Call:
## lm(formula = Life_expectancy ~ GDP_per_capita, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47292 -0.20920  0.07593  0.14890  0.40610 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    7.023e+01  3.418e-01  205.45  < 2e-16 ***
## GDP_per_capita 1.864e-03  1.088e-04   17.13  8.7e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2859 on 14 degrees of freedom
## Multiple R-squared:  0.9544, Adjusted R-squared:  0.9512 
## F-statistic: 293.4 on 1 and 14 DF,  p-value: 8.697e-11

resettest(model_1)

## 
##  RESET test
## 
## data:  model_1
## RESET = 44.312, df1 = 2, df2 = 12, p-value = 2.877e-06

Rozšírený model

model_2 <- lm(Life_expectancy ~ GDP_per_capita + Schooling + Adult_mortality, data = udaje)
summary(model_2)

## 
## Call:
## lm(formula = Life_expectancy ~ GDP_per_capita + Schooling + Adult_mortality, 
##     data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.25267 -0.10272 -0.00742  0.09824  0.35555 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     5.237e+01  7.032e+00   7.447 7.76e-06 ***
## GDP_per_capita  1.176e-03  2.551e-04   4.609 0.000601 ***
## Schooling       1.977e+00  5.252e-01   3.763 0.002705 ** 
## Adult_mortality 2.043e-02  3.734e-02   0.547 0.594406    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1998 on 12 degrees of freedom
## Multiple R-squared:  0.9809, Adjusted R-squared:  0.9762 
## F-statistic: 205.9 on 3 and 12 DF,  p-value: 1.395e-10

resettest(model_2)

## 
##  RESET test
## 
## data:  model_2
## RESET = 30.261, df1 = 2, df2 = 10, p-value = 5.733e-05

model_3 <- lm(Life_expectancy ~  Alcohol_consumption+ BMI+ Adult_mortality, data = udaje)
summary(model_3)

## 
## Call:
## lm(formula = Life_expectancy ~ Alcohol_consumption + BMI + Adult_mortality, 
##     data = udaje)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.128981 -0.049255 -0.004123  0.047679  0.174974 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         14.06603    5.48110   2.566  0.02471 *  
## Alcohol_consumption -0.17772    0.04208  -4.224  0.00118 ** 
## BMI                  2.51203    0.16260  15.449 2.77e-09 ***
## Adult_mortality     -0.02686    0.01522  -1.765  0.10304    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08956 on 12 degrees of freedom
## Multiple R-squared:  0.9962, Adjusted R-squared:  0.9952 
## F-statistic:  1040 on 3 and 12 DF,  p-value: 9.264e-15

resettest(model_3)

## 
##  RESET test
## 
## data:  model_3
## RESET = 28.053, df1 = 2, df2 = 10, p-value = 7.921e-05

model_many <- lm(Life_expectancy ~ GDP_per_capita +Alcohol_consumption+ BMI + Adult_mortality + Schooling + Year, data = udaje)
summary(model_many)

## 
## Call:
## lm(formula = Life_expectancy ~ GDP_per_capita + Alcohol_consumption + 
##     BMI + Adult_mortality + Schooling + Year, data = udaje)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.033643 -0.022957 -0.005425  0.025682  0.040025 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -1.078e+03  1.693e+02  -6.367  0.00013 ***
## GDP_per_capita       2.029e-04  1.527e-04   1.329  0.21664    
## Alcohol_consumption -2.105e-01  2.809e-02  -7.494 3.71e-05 ***
## BMI                 -3.348e+00  7.656e-01  -4.373  0.00179 ** 
## Adult_mortality      5.283e-02  1.410e-02   3.748  0.00457 ** 
## Schooling            1.873e-01  1.447e-01   1.294  0.22786    
## Year                 6.151e-01  9.318e-02   6.601 9.91e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03287 on 9 degrees of freedom
## Multiple R-squared:  0.9996, Adjusted R-squared:  0.9994 
## F-statistic:  3873 on 6 and 9 DF,  p-value: 7.896e-15

resettest(model_many)

## 
##  RESET test
## 
## data:  model_many
## RESET = 2.8365, df1 = 2, df2 = 7, p-value = 0.1252

Porovnanie modelov pomocou AIC

AIC(model_1, model_2, model_3, model_many)

Automatický výber podľa AIC

model_start <- lm(Life_expectancy ~ GDP_per_capita + Schooling + BMI +
                    Adult_mortality + Alcohol_consumption + Year, data = udaje)

model_aic <- stepAIC(model_start, direction = "both", trace = FALSE)
summary(model_aic)

## 
## Call:
## lm(formula = Life_expectancy ~ GDP_per_capita + Schooling + BMI + 
##     Adult_mortality + Alcohol_consumption + Year, data = udaje)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.033643 -0.022957 -0.005425  0.025682  0.040025 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -1.078e+03  1.693e+02  -6.367  0.00013 ***
## GDP_per_capita       2.029e-04  1.527e-04   1.329  0.21664    
## Schooling            1.873e-01  1.447e-01   1.294  0.22786    
## BMI                 -3.348e+00  7.656e-01  -4.373  0.00179 ** 
## Adult_mortality      5.283e-02  1.410e-02   3.748  0.00457 ** 
## Alcohol_consumption -2.105e-01  2.809e-02  -7.494 3.71e-05 ***
## Year                 6.151e-01  9.318e-02   6.601 9.91e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03287 on 9 degrees of freedom
## Multiple R-squared:  0.9996, Adjusted R-squared:  0.9994 
## F-statistic:  3873 on 6 and 9 DF,  p-value: 7.896e-15

resettest(model_aic)

## 
##  RESET test
## 
## data:  model_aic
## RESET = 2.8365, df1 = 2, df2 = 7, p-value = 0.1252

AIC(model_start, model_aic)

Porovnanie lineárneho a logaritmického tvaru

model_linear <- lm(Life_expectancy ~ GDP_per_capita, data = udaje)
model_log <- lm(log(Life_expectancy) ~ GDP_per_capita+ Schooling + BMI +
                    Adult_mortality + Alcohol_consumption, data = udaje)

summary(model_linear)

## 
## Call:
## lm(formula = Life_expectancy ~ GDP_per_capita, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.47292 -0.20920  0.07593  0.14890  0.40610 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    7.023e+01  3.418e-01  205.45  < 2e-16 ***
## GDP_per_capita 1.864e-03  1.088e-04   17.13  8.7e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2859 on 14 degrees of freedom
## Multiple R-squared:  0.9544, Adjusted R-squared:  0.9512 
## F-statistic: 293.4 on 1 and 14 DF,  p-value: 8.697e-11

summary(model_log)

## 
## Call:
## lm(formula = log(Life_expectancy) ~ GDP_per_capita + Schooling + 
##     BMI + Adult_mortality + Alcohol_consumption, data = udaje)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0011733 -0.0004857 -0.0001523  0.0005365  0.0018619 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          3.878e+00  1.430e-01  27.113 1.08e-10 ***
## GDP_per_capita       1.025e-05  3.850e-06   2.663 0.023775 *  
## Schooling           -5.866e-04  4.119e-03  -0.142 0.889577    
## BMI                  1.845e-02  6.343e-03   2.908 0.015618 *  
## Adult_mortality     -3.978e-04  1.973e-04  -2.016 0.071443 .  
## Alcohol_consumption -3.781e-03  7.418e-04  -5.097 0.000466 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0009722 on 10 degrees of freedom
## Multiple R-squared:  0.9978, Adjusted R-squared:  0.9967 
## F-statistic:   917 on 5 and 10 DF,  p-value: 5.71e-13

AIC(model_linear, model_log)

Aplikácia RESET testu

resettest(model_linear, power = 2, type = "fitted")

## 
##  RESET test
## 
## data:  model_linear
## RESET = 13.846, df1 = 1, df2 = 13, p-value = 0.002565

RESET test pri lepšom modeli

resettest(model_log, power = 2, type = "fitted")

## 
##  RESET test
## 
## data:  model_log
## RESET = 3.6796, df1 = 1, df2 = 9, p-value = 0.0873

LOESS (LOcal regrESSion) krivka predstavuje lokálne vyhladený trend medzi premennými bez predpokladu konkrétneho funkčného tvaru.↩︎

Albánsko - štatistika

Natália Kunová

Úvod do regresie

Využívanie niektorých knižníc

Príprava údajov - import z csv súboru

Diagnostické grafy

Residuals vs Fitted

Test na normalitu a test na autokoreláciu

Q-Q plot

Scale-Location plot

Residuals vs Leverage

Test na odľahlé hodnoty (z knižnice car)

Breusch–Pagan

Breusch–Godfrey

Robustné štandardné chyby

Chyba špecifikácie

Jednoduchý príklad

Rozšírený model

Porovnanie modelov pomocou AIC

Automatický výber podľa AIC

Porovnanie lineárneho a logaritmického tvaru

Aplikácia RESET testu

RESET test pri lepšom modeli