Econometrics in R - cvičenie 5

S využitím databázy WHO Life Expactancy Data database.

Pri ďalšej práci budeme používať knižnice

library(zoo)
library(tseries)
library(lmtest)
library(sandwich)
library(car)
rm(list=ls())

V hlavnom menu som Session nastavil na Source File Location. Môžeme to urobiť aj interaktívnym spôsobom - napísať dole do Console (alebo to zahrnúť priamo do skriptu) ako

setwd(“/Cloud/project/tyzdne/tyzden5”)

Organizácia priečinkov a podpriečinkov sa však môže líšiť v závislosti od projektu, preto som ich nezaradil do Chunk-u.

Údaje o očakávanej dĺžke života sú usporiadané v súbore csv, stĺpce sú oddelené znakom “,” a používajú desatinnú čiarku. V pracovnom priečinku som vytvoril podpriečinok s názvom udaje, aby som oddelil údaje od zvyšku projektu. Môj priečinok sa tiež nazýva udaje.

Nie všetky údaje budú použité, preto som vybral len niektoré stĺpce pre neskoršie použitie.

Úvod do problému, stanovenie hypotéz

Rozhodol som sa modelovať strednú dĺžku života Life.expectancy v závislosti od troch vysvetľujúcich premenných a to BMI, HDP na obyvateľa GDP a stredný počet rokov štúdia Schooling.

Naša pracovná hypotéza hovorí o štatisticky významnom vplyve všetkých troch vysvetľujúcich premenných, pričom u premenných GDP a Schooling by malo ísť o pozitívny vplyv (očakávame kladné znamienko odhadovaného regresného koeficienta) a v prípade BMI by malo ísť of negatívny vplyv (so záporným znamienkom)

Príprava databázy, čistenie a úprava údajov

Keďže niektoré údaje chýbajú, doplnil som ich mediánovými hodnotami premennej, ktorú zvažujem. lm

udaje <- read.csv("udaje/Life_Expectancy_Data.csv",dec=".",sep=",",header = TRUE)
# select just the record from 2015
udaje.2015 <- udaje[udaje$Year==2015,c("Life.expectancy","BMI","GDP","Schooling")]

# data imputation

# Compute column medians
#column_medians <- sapply(udaje.2015, median, na.rm = TRUE)

# Impute missing values with column medians
# Compute column medians
column_medians <- sapply(udaje.2015, median, na.rm = TRUE)

# Impute missing values with column medians
udaje_imputed <- udaje.2015
for (col in names(udaje.2015)) {
  udaje_imputed[[col]][is.na(udaje_imputed[[col]])] <- column_medians[col]
}

udaje.2015 <- udaje_imputed

Teraz chceme vidieť tvar údajov (či nie sú v nich nejaké nezrovnalosti – napríklad hodnoty 0).

# Suppose udaje.2015 is your data frame

# Determine number of plots
num_plots <- length(names(udaje.2015))

# Set the layout: 2 rows × 2 columns
par(mfrow = c(2, 2))
par(mar = c(4, 4, 2, 1))  # Adjust margins

# Loop through columns and plot each boxplot with variable name as title
for (col in names(udaje.2015)) {
  boxplot(udaje.2015[[col]], 
          main = col,           # variable name as title
          xlab = "Value", 
          col = "lightblue")
}

# Add a global title
mtext("Boxploty jednotlivých premenných", outer = TRUE, cex = 1.4, font = 2)

Lineárna regresia

Model odhadujeme príkazom lm()

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

Objekt triedy lm() nám poskytuje niekoľko výsledkov:

1. Vector odhadnutých koeficientov model$coefficients
2. Vektor rezíduí model$ residuals
3. Vektor vyrovnaných hodnôt vysvetľovanej veličiny model$fitted.values
4. Maticu X model$x

#print("Odhadnuté koeficienty sú: ")
#      print(model$coefficients)
#print("Odhadnuté rezíduá: ")
#print(model$residuals)
#print("Vyrovnané hodnoty vysvetľovanej premennej sú: ")
#print(model$fitted.values)
#print("matica model$xlevels: ")
#print(model.matrix(model))
#X <- model.matrix(model)
#diag(X %*% solve(t(X) %*% X) %*% t(X))

summary(model)

Call:
lm(formula = Life.expectancy ~ +1 + 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

Súhrn odhadovaného modelu nám poskytuje súbor odhadovaných regresných koeficientov, ktorých znamienka budú rozoberané neskôr. Ak hovoríme o vlastnostiach modelu ako celku, pozrime sa najskôr na nasledujúce obrázky. Na základe Q-Q grafu získavame dojem o možných problémoch porušenia normality rezíduí.

# Nastaviť rozloženie 2 x 2
par(mfrow = c(2, 2))

# Vykresliť všetky 4 diagnostické grafy modelu
plot(model)

# (Voliteľné) pridať spoločný nadpis
#mtext("Diagnostické grafy regresného modelu", outer = TRUE, cex = 1.2, font = 2)

# Resetovať layout
par(mfrow = c(1, 1))

Residuals vs. fitted

Interpretácia vášho konkrétneho grafu

• Centrovanie okolo nuly

• Reziduály kolíšu približne okolo 0 – to je dobré.

• Naznačuje to, že model nemá výrazné skreslenie v predikciách.

• Červená hladká čiara je mierne zakrivená (vpravo sa ohýba nadol a uprostred nahor), to naznačuje možnú miernu nelinearitu – modelu môže chýbať nelineárny tvar nejakej premennej.

• Rozptyl rezíduí - Vertikálny rozptyl (variancia rezíduí) sa javí ako približne konštantný v rámci prispôsobených hodnôt čo je dobrý dôkaz homoscedasticity.

• Ak by sa rezíduá rozprestierali (tvorili kónus), naznačovalo by to heteroskedasticitu.

• Odľahlé hodnoty - Niekoľko bodov (napr. v blízkosti −20) leží ďaleko od ostatných – ide o potenciálne odľahlé hodnoty alebo vplyvné pozorovania.

Môžeme ich preskúmať pomocou outlierTest(model) (z balíka car).

Q-Q plot

Čo ukazuje

• Os X: Teoretické kvantily – čo by sme očakávali, ak by rezíduá boli dokonale normálne rozložené.
• Os Y: Štandardizované rezíduá – skutočné kvantily z vašej vzorky.
• 45° prerušovaná čiara: Ideálny prípad – ak sú rezíduá normálne rozložené, body by mali ležať tesne pozdĺž tejto čiary.

Interpretácia vášho konkrétneho grafu

Celkový tvar

Väčšina bodov leží blízko priamky — to naznačuje, že rezíduá sú približne normálne.

To je dobré: predpoklad normality sa zdá byť vo veľkej miere splnený.

Krajné hodnoty (extrémy)

Body na oboch koncoch (vľavo dole a vpravo hore) sa mierne odchyľujú od priamky.

To naznačuje miernu nenormálnosť v koncoch – možno niekoľko odľahlých hodnôt alebo ťažšie konce, ako je normálne (trochu špicatosť).

Stredná oblasť

Stredná časť grafu (−1 až +1 kvantily) sa veľmi dobre zhoduje – čo znamená, že väčšina rezíduí pekne zapadá do normálneho rozdelenia.

Zároveň vykonávame následné testy, v ktorých zamietame hypotézu o normálnom rozložení modelových chýb, ako aj prítomnosť extrémnych hodnôt.

Scale location plot

Čo to znázorňuje

Os X: Vyrovnan0 hodnoty Os Y: Druhá odmocnina absolútnych štandardizovaných rezíduí Červená čiara: LOESS (vyhladený) trend cez body

Interpretácia nášho konkrétneho grafu

Horizontálne rozptýlenie

• Body sú rovnomerne rozptýlené po osi x bez vytvorenia lievika alebo krivky. To naznačuje približne konštantnú varianciu – rezíduá sú homoscedastické.
• Červená hladká čiara je takmer rovná – ďalší znak, že pri zvyšovaní vyrovnaných hodnôt nedochádza k žiadnej systematickej zmene variancie.
• Odľahlé hodnoty - niekoľko bodov je mierne nad 1,5, ale žiadny z nich nie je extrémny – takže nedochádza k žiadnym závažným anomáliám variancie.

residuals vs leverage

Čo znázorňuje graf

Os X: Pákový efekt — meria, ako ďaleko je prediktorový vektor bodu od stredu všetkých 𝑥 x.

Os Y: Štandardizované rezíduá — ako ďaleko je pozorovaná hodnota od vyrovnanej hodnoty v jednotkách štandardnej odchýlky.

Bodkované krivky: Kontúry Cookovej vzdialenosti, ktoré udávajú, do akej miery pozorovanie ovplyvňuje regresnú priamku.

Červená čiara: Hladká LOESS čiara prechádzajúca rezíduami.lm

Interpretácia vášho konkrétneho grafu

Rozloženie vplyvu

Väčšina pozorovaní má nízky vplyv (pod 0,05) — typické pre veľké vzorky alebo dobre vyvážené údaje.

Jeden alebo dva body (napr. okolo 0,2) vynikajú – ide o pozorovania s vysokou pákou, čo znamená, že ich hodnoty sú ďaleko od väčšiny údajov.

Veľkosť rezíduí

Štandardizované rezíduá väčšinou medzi −2 a +2 – to je dobré (žiadne závažné výnimky v 𝑦 y).

Pozorovanie označené číslom 113 má zdanlivo stredný vplyv a relatívne veľkú rezíduálnu hodnotu – potenciálne vplyvný prípad.

Kontúry Cookovej vzdialenosti

Žiaden z bodov jasne neprekračuje vonkajšie línie Cookovej vzdialenosti (≈0,5 alebo 1,0).

Preto sa nezdá, že by niektoré pozorovanie neprimerane ovplyvňovalo regresné koeficienty.

{r}lm # normality tests residuals <- residuals(model) jb_test <- jarque.bera.test(residuals) jb_test # outlier test (see p-value for Bonferroni correction) outlier_test <- outlierTest(model) outlier_test

Keďže sa nepreukázala normalita rezíduí, pokúsme sa eliminovať odľahlé hodnoty v prípade GDP - dokážeme to logaritmickou transformáciou tejto premennej a vylúčením BMI, ktoré sa ukázalo ako neinterpretovateľné. Nová regresia bude mať tvar

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

Call:
lm(formula = Life.expectancy ~ +1 + 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

# Nastaviť rozloženie 2 x 2
par(mfrow = c(2, 2))

# Vykresliť všetky 4 diagnostické grafy modelu
plot(model2)

# (Voliteľné) pridať spoločný nadpis
#mtext("Diagnostické grafy regresného modelu", outer = TRUE, cex = 1.2, font = 2)

# Resetovať layout
par(mfrow = c(1, 1))

# normality tests
residuals <- residuals(model)
jb_test <- jarque.bera.test(residuals)
jb_test

    Jarque Bera Test

data:  residuals
X-squared = 21.391, df = 2, p-value = 2.265e-05

# outlier test (see p-value for Bonferroni correction)
outlier_test <- outlierTest(model)
outlier_test

No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:

Záver o analýze odľahlých hodnôt

Premenné GDP, a Schooling predlžujú štatisicky významne strednú dĺžku života. Na druhej strane BMI nám dávalo neinterpretovateľné výsledky. Rezíduá nevykazujú normálne rozdelenie, keďže však máme veľké množstvo pozorovaní, aj naďalej budeme pracovať s týmito údajmi. V modeli sa nepreukazujú žiadne významné nelinearity.

Heteroskedasticita

Prítomnosť heteroskedasticity (nekonštantného rozptylu náhodnej zložky) spôsobuje zlé vyhodnocovanie t-testov významnosti jednotlivých regresných koeficientov. Preto je nutné, aby sme heteroskedasticitu - detekovali (vizuálne a s pomocou testov) - a v prípade prítomnosti heteroskedasticity aby sme ju odstránili.

Aj v našom prípade by sme sa mohli pokúsiť o vizuálne vyhodnotenie nasledovných grafov (aj keď jeden graf sme už skúmali - bol to tzv. Scale-Location grafy uvedené vyššie).

Tentokrát sa pokúsime o vizuálne znázornenie závislosti štvorcov rezíduí a vysvetľujúcej premennej, u ktorej máme podozrenie, že môže heteroskedasticitu spôsobovať. Budeme posudzovať dva modely - a to model nazvaný model alebo model nazvaný model2. model2 má zlogaritmizovanú premennú GDP, čo sme robili z dôvodu odstránenia vplyvu odľahlých premenných v predchádzajúcich krokoch, model je pôvodným modelom.

library(ggplot2)
library(patchwork)  # install.packages("patchwork")

p1 <- ggplot(udaje.2015, aes(x = GDP, y = resid(model)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "GDP", 
       y = "Squared Residuals",
       title = "Sqiared Residuals vs GDP") +
  theme_minimal()

p2 <- ggplot(udaje.2015, aes(x = Schooling, y = resid(model)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Schooling", 
       y = "Squared Residuals",
       title = "Squared Residuals vs Schooling") +
  theme_minimal()

# Combine side by side
p1 + p2

a teraz model so zlogaritmizovanou premennou GDP.

library(ggplot2)
library(patchwork)  # install.packages("patchwork")

p1 <- ggplot(udaje.2015, aes(x = log(GDP), y = resid(model2)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "log(GDP)", 
       y = "Squared Residuals",
       title = "Residuals vs log(GDP)") +
  theme_minimal()

p2 <- ggplot(udaje.2015, aes(x = Schooling, y = resid(model2)^2)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, color = "red") +
  labs(x = "Schooling", 
       y = "Squared Residuals",
       title = "Residuals vs Schooling") +
  theme_minimal()

# Combine side by side
p1 + p2

Na tomto obrázku podľa vyhladených hodnôť štvorcov rezíduí (červená krivka) môžeme konštatovať, že neobsahuje žiaden významný vývoj s vysvetľujúcou premennou (či už \(log(GDP)\), alebo \(Schooling\)). Kvôli demonštrácii ale ukážme, že bez predchádzajúcej logaritmickej transformácie by to dopadlo inak, t.j. vychádzajme z pôvodného modelu označeného ako model nasledovne:

Testovanie prítomnosti heteroskedasticity

# Install (if not yet installed)
# install.packages("lmtest")

# Load the package
library(lmtest)

# Run the Breusch–Pagan test
bptest(model)

    studentized Breusch-Pagan test

data:  model
BP = 15.211, df = 3, p-value = 0.001645

# Install (if not yet installed)
# install.packages("lmtest")

# Load the package
library(lmtest)

# Run the Breusch–Pagan test
bptest(model2)

    studentized Breusch-Pagan test

data:  model2
BP = 3.9542, df = 2, p-value = 0.1385

Na základe výsledkov regresie môžeme povedať, že heteroskedasticita rezíduí nie je v model2 prítomná, ale v prípade model prétomná je. Ak by ale prítomná bola, a logaritmizácia premenných, alebo odstránenie odľahlých premenných by nám nepomohli, môžeme to riešiť s pomocou takzvanej White heteroskedasticity Consistent matrix, kde v t testoch významnosti regresných koeficientov sa používajú “hrubšie” odhady rozptylov regresných koeficientov. Urobíme to nasledovne

#install.packages("sandwich")
#install.packages("lmtest")
library(sandwich)
library(lmtest)
coeftest(model, vcov = vcovHC(model))

t test of coefficients:

              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 4.4194e+01 1.8560e+00 23.8108  < 2e-16 ***
BMI         4.9479e-02 2.4034e-02  2.0587  0.04097 *  
GDP         6.9722e-05 3.1291e-05  2.2282  0.02711 *  
Schooling   1.9210e+00 1.7450e-01 11.0087  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Všimnime si, že tentokrát je už každá premenná štatisticky významná. Na druhej strane treba podotknúť, že použitie tejto metódy vyžaduje veľký počet pozorovaní (> 100)

Menej používanou je možnosť, kedy predelíme všetky premenné v dátovom sete premennou, ktorá heteroskedasticitu spôsobuje, čo ale často vedie k zhoršeniu interpretačnej schopnosti modelu. V prípade napr. GDP, množstva opyvateľov a rozlohy krajiny to ale napríklad možné je. Ak napríklad (podľa grafov) zistíme, že premenná množstvo obuvateľov spôsobuje heteroskedasticitu, vypočítame GDP na obyvateľa, resp. počet km štvorcových na jedného obyvateľa a môžeme dostať nový model s odstránenou heteroskedasticitou.