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 (optional)

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

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

# Reset layout to default (1 plot per figure)
par(mfrow = c(1, 1))

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|:

##Conclusion##

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.