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)

Error in file(file, "rt") : cannot open the connection

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

Lineárna regresia

Model odhadujeme príkazom lm()

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

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

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

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.

a teraz model so zlogaritmizovanou premennou GDP.

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

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

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.