Econometrics in R - cvičenie 5

Základná diagnostika modelu

Author

V. Gazda

Pokračujeme s využívaním WHO Life Expactancy Data databázy.

library(broom)
library(knitr)
library(kableExtra)
library(ggplot2)
library(dplyr)

Príprava databázy

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

Uloženie všetkých objektov môjho workspace, čistenie workspace

  • Ak nechcem do každého projektu dávať opätovné importovanie a čistenie databázy, odožím si celý obsah - všetky objekty - pamäte (workspace) do súboru a ten na začiatku zadania nabudúce nadčítam:
save.image("pamat.Rdata")
  • Mazanie niektorých objektov z pamäte:
aaa <- 6
bbb <- "text"     #    aaa a bbb chcem vymazat

rm(list = c("aaa", "bbb"))  # vynazem tieto dva objekty


ls()  # v pamati zostalo
[1] "col"            "column_medians" "udaje"          "udaje_imputed" 
[5] "udaje.2015"
  • Čistenie pamäte - celej
rm(list=ls())
ls()
character(0)
  • Opätovne si nadčítam uložený obsah pamäte - workspace
load("pamat.Rdata")

Samozrejme, load() budem v budúcnosti používať na začiatku dokumentu a save.image() na konci dokumentu - tak, aby som si prenášal celý obsah pamäte z jedného cvičenia na druhé.

To predpokladá, že si budete robiť priebežne poriadok v pamäti a nepotrebné objekty budete mazať príkazom rm(). - pozrite príklad uvedený vyššie.

Lineárna regresia

Odhadujeme regresný model:

\[ LifeExpectancy_i = \beta_0+\beta_1 BMI_i + \beta_2 GDP_i + \beta_3 Schooling_i+u_i \]

Model odhadujeme príkazom lm(). V našom prípade predpokladáme, že

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

Výsledky regresie - odhady regresných koeficientov - môžeme vidieť v tabuľke

tidy(model) |>
  mutate(
    Significance = case_when(
      p.value < 0.001 ~ "***",
      p.value < 0.01  ~ "**",
      p.value < 0.05  ~ "*",
      p.value < 0.10  ~ ".",
      TRUE ~ ""
    )
  ) |>
  kable(
    digits = 4,
    col.names = c(
      "Variable",
      "Estimate",
      "Std. Error",
      "t value",
      "p-value",
      "Significance"
    ),
    caption = "Regression coefficients",
    escape = TRUE
  ) |>
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = FALSE
  )
Regression coefficients
Variable Estimate Std. Error t value p-value Significance
(Intercept) 44.1936 1.8389 24.0332 0.0000 ***
BMI 0.0495 0.0214 2.3078 0.0222 *
GDP 0.0001 0.0000 1.8132 0.0715 .
Schooling 1.9210 0.1645 11.6762 0.0000 ***

Interpretácia výsledkov regresného modelu

Funkcia broom::tidy(model) poskytuje tabuľku odhadnutých regresných koeficientov a výsledkov testovania ich štatistickej významnosti.

Vysvetlenie jednotlivých stĺpcov

Stĺpec Vysvetlenie
Variable Názov premennej, ku ktorej sa vzťahuje regresný koeficient. (Intercept) predstavuje absolútny člen regresnej rovnice.
Estimate Bodový odhad regresného koeficientu \(\hat{\beta}_j\). Vyjadruje zmenu očakávanej hodnoty vysvetľovanej premennej pri jednotkovej zmene vysvetľujúcej premennej, za predpokladu nezmenených ostatných premenných.
Std. Error Štandardná chyba odhadu regresného koeficientu \(SE(\hat{\beta}_j)\). Vyjadruje neistotu odhadu. Čím je štandardná chyba menšia, tým je odhad presnejší.
t value Hodnota testovacej štatistiky Studentovho t-testu. Používa sa na testovanie hypotézy o nulovej hodnote regresného koeficientu.
p-value Pravdepodobnosť získania rovnako alebo viac extrémnej hodnoty testovacej štatistiky za predpokladu platnosti nulovej hypotézy. Malá p-hodnota predstavuje dôkaz proti nulovej hypotéze.
Significance Symbolické vyjadrenie štatistickej významnosti regresného koeficientu pomocou hviezdičiek.

Testovanie štatistickej významnosti koeficientov

Pre každý regresný koeficient testujeme hypotézy:

\[ H_0:\beta_j=0 \]

\[ H_1:\beta_j\neq 0 \]

Testovacia štatistika má tvar:

\[ t_j=\frac{\hat{\beta}_j}{SE(\hat{\beta}_j)} \]

Za platnosti klasických predpokladov vrátane normality náhodných chýb má testovacia štatistika pri platnosti nulovej hypotézy Studentovo t-rozdelenie s \(n-k\) stupňami voľnosti, kde \(n\) predstavuje počet pozorovaní a \(k\) počet odhadovaných regresných koeficientov vrátane absolútneho člena.

Rozhodovacie pravidlo: Ak je p-hodnota menšia ako zvolená hladina významnosti \(\alpha\), zamietame nulovú hypotézu.

Označenie štatistickej významnosti

Symbol p-hodnota Interpretácia
*** \(p < 0.001\) Štatisticky významný na hladine 0,1 %
** \(p < 0.01\) Štatisticky významný na hladine 1 %
* \(p < 0.05\) Štatisticky významný na hladine 5 %
. \(p < 0.10\) Štatisticky významný na hladine 10 %
Bez symbolu \(p \geq 0.10\) Štatisticky nevýznamný na hladine 10 %

Interpretácia koeficientov nášho modelu

Koeficient pri GDP:

Pri zvýšení HDP na obyvateľa o 1 dolár sa očakávaná dĺžka života zmení približne o \(\hat{\beta}_1/1000\) roka, za predpokladu nezmenených hodnôt ostatných premenných - ceteris paribus`. Dá sa to ale lepšie interpretovať aj ako ak HDP na obyvateľa narastie o 1000 dolárov, potom očakávaná dĺžka dožitia narastie o 1 rok.

Koeficient pri Schooling:

Ak je školská dochádzka meraná v rokoch, zvýšenie jej hodnoty o jeden rok je spojené so zmenou očakávanej dĺžky života takmer o \(2\) roky (ceteris paribus).

Absolútny člen (Intercept):

Koeficient \(\hat{\beta}_0\) predstavuje predikovanú hodnotu očakávanej dĺžky života pri nastavení všetkých vysvetľujúcich premenných na nulovú hodnotu. V našom modeli má skôr význam, pretože uvedená kombinácia hodnôt nemusí byť ekonomicky realistická. Ide o priesečník regresnej nadroviny so zvislou osou.

Koeficient pri BMI Ak premenná BMI narastie o 1 jednotku, potom sa očakávaná dĺžka života zvýši takmer o 5 stotín roka - t.j. približne o 15 dní. Je to ekonomicky neinterpretovateľné, preto sme túto premennú v ďaších analýzach vylúčili.

Kvalita vyrovnania modelu

Kvalita vyrovnania sa dá posúdiť pomocou \(R^2\), ale aj výsledku F-testu (pozri p.value vedľa F stat.)

# Model summary without deviance
glance(model) |>
  select(-deviance) |>
  kable(
    digits = 4,
    col.names = c(
      "R.squared",
      "Adj.R.squared",
      "Sigma",
      "F Stat.",
      "p.value",
      "df",
      "logLik",
      "AIC",
      "BIC",
      "df.residual",
      "N"
    ),
    caption = "Model summary"
  ) |>
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = FALSE
  )
Model summary
R.squared Adj.R.squared Sigma F Stat. p.value df logLik AIC BIC df.residual N
0.6224 0.6161 5.0334 98.3632 0 3 -553.3884 1116.777 1132.824 179 183

Interpretácia súhrnných charakteristík regresného modelu

Funkcia broom::glance(model) poskytuje súhrnné štatistické charakteristiky odhadnutého regresného modelu. Na rozdiel od funkcie tidy(), ktorá poskytuje informácie o jednotlivých regresných koeficientoch, funkcia glance() charakterizuje model ako celok.

Vysvetlenie jednotlivých stĺpcov

Stĺpec Vysvetlenie
R.squared Koeficient determinácie \(R^2\). Vyjadruje podiel variability vysvetľovanej premennej, ktorý vysvetľuje regresný model.
Adj.R.squared Upravený koeficient determinácie \(\bar R^2\). Zohľadňuje počet vysvetľujúcich premenných a počet pozorovaní. Penalizuje pridávanie premenných, ktoré dostatočne nezlepšujú vysvetľovaciu schopnosť modelu.
Sigma Reziduálna štandardná chyba \(\hat\sigma\). Predstavuje odhad smerodajnej odchýlky náhodnej zložky regresného modelu.
F Stat. Hodnota F-štatistiky testujúcej spoločnú štatistickú významnosť vysvetľujúcich premenných.
p.value P-hodnota celkového F-testu regresného modelu. Používa sa na rozhodnutie o zamietnutí nulovej hypotézy, že všetky smernicové regresné koeficienty sú nulové.
df Počet vysvetľujúcich premenných v modeli, bez absolútneho člena. Zodpovedá počtu stupňov voľnosti v čitateli celkového F-testu.
logLik Hodnota logaritmu vierohodnostnej funkcie modelu. Používa sa pri porovnávaní modelov a výpočte informačných kritérií.
AIC Akaikeho informačné kritérium. Hodnotí model na základe jeho vierohodnosti a počtu odhadovaných parametrov. Pri porovnávaní modelov uprednostňujeme nižšiu hodnotu AIC.
BIC Bayesovské informačné kritérium. Podobá sa AIC, ale spravidla prísnejšie penalizuje zložitosť modelu, najmä pri väčšom počte pozorovaní.
df.residual Reziduálne stupne voľnosti \(n-k\), kde \(n\) je počet pozorovaní a \(k\) počet odhadovaných regresných koeficientov vrátane absolútneho člena.
N Počet pozorovaní použitých pri odhade regresného modelu.
deviance Reziduálny súčet štvorcov (RSS). Táto charakteristika bola z výslednej tabuľky vynechaná, ale je dôležitá pri výpočte reziduálnej štandardnej chyby a koeficientu determinácie.

1. Koeficient determinácie – R.squared

Koeficient determinácie vyjadruje, akú časť variability vysvetľovanej premennej dokážeme vysvetliť pomocou regresného modelu.

\[ R^2=1-\frac{RSS}{TSS} \]

kde:

  • \(RSS=\sum_{i=1}^{n}(y_i-\hat y_i)^2\) je reziduálny súčet štvorcov,
  • \(TSS=\sum_{i=1}^{n}(y_i-\bar y)^2\) je celkový súčet štvorcov.

Pri OLS modeli s absolútnym členom platí:

\[ 0\leq R^2\leq 1 \]

Interpretácia: Ak \(R^2=0.62\), potom regresný model vysvetľuje 62 % variability vysvetľovanej premennej. Zvyšných 38 % variability zostáva nevysvetlených modelom.

Vysoká hodnota \(R^2\) sama osebe nezaručuje správnu špecifikáciu modelu ani kauzálny vzťah medzi premennými.

2. Upravený koeficient determinácie – Adj.R.squared

Upravený koeficient determinácie zohľadňuje počet vysvetľujúcich premenných.

\[ \bar R^2 = 1-(1-R^2)\frac{n-1}{n-k} \]

kde \(k\) predstavuje počet odhadovaných regresných koeficientov vrátane absolútneho člena.

Na rozdiel od obyčajného \(R^2\) nemusí po pridaní ďalšej vysvetľujúcej premennej vzrásť.

Interpretácia: Upravený koeficient determinácie je vhodný pri porovnávaní modelov s rôznym počtom vysvetľujúcich premenných a rovnakou vysvetľovanou premennou, odhadnutých na rovnakých pozorovaniach.

Jeho hodnota môže byť aj záporná. V našom prípade je táto hodnota blézka pôvodnému koeficientu determinácie, nakoľko máme veľké množstvo pozorovaní a málo vysvetľujúcich premenných.

3. Reziduálna štandardná chyba – Sigma

Reziduálna štandardná chyba predstavuje odhad smerodajnej odchýlky náhodnej zložky regresného modelu.

\[ \hat\sigma= \sqrt{\frac{\sum_{i=1}^{n}e_i^2}{n-k}} \]

alebo

\[ \hat\sigma=\sqrt{\frac{RSS}{n-k}} \]

Interpretácia: Vyjadruje typickú veľkosť reziduálnych odchýlok od regresnej funkcie. Má rovnakú jednotku merania ako vysvetľovaná premenná.

V našom modeli očakávanej dĺžky života je vyjadrená v rokoch a hovoríme, že chyba modelu je v priemere približne 5 rokov (teda odhadovaná smerodajná odchýlka náhodnej zložky predstavuje približne 5 rokov).

4. F-štatistika – F Stat.

F-test overuje spoločnú štatistickú významnosť vysvetľujúcich premenných.

Pre náš model:

\[ LifeExpectancy_i= \beta_0+\beta_1\ln(GDP_i)+\beta_2Schooling_i+u_i \]

testujeme hypotézy:

\[ H_0:\beta_1=\beta_2=0 \]

\[ H_1:\text{aspoň jeden z koeficientov }\beta_1,\beta_2\neq 0 \]

Testovacia štatistika má tvar:

\[ F= \frac{R^2/(k-1)} {(1-R^2)/(n-k)} \]

Pri platnosti nulovej hypotézy a klasických predpokladov vrátane normality chýb má F-rozdelenie s \(k-1\) a \(n-k\) stupňami voľnosti.

Interpretácia: Vysoká hodnota F-štatistiky predstavuje dôkaz proti nulovej hypotéze, že všetky vysvetľujúce premenné majú nulové smernicové koeficienty. Táto štatistika nemá špecifickú ekonomickú interpretáciu a má technický charakter.

5. P-hodnota celkového F-testu – p.value

P-hodnota vyjadruje pravdepodobnosť získania rovnako alebo viac extrémnej hodnoty F-štatistiky za predpokladu platnosti nulovej hypotézy

\[H_0: \beta_i = 0 \quad \forall (i>0)\] \[H_0: (\exists i>0) \quad \beta_i \neq 0 \]

Tradičné rozhodovacie pravidlo:

  • Ak \(p<0.05\), zamietame nulovú hypotézu na 5 % hladine významnosti.
  • Ak \(p\geq 0.05\), nulovú hypotézu na tejto hladine významnosti nezamietame.

Zamietnutie nulovej hypotézy znamená, že aspoň jedna vysvetľujúca premenná má štatisticky významný smernicový koeficient.

Pozor: Táto p-hodnota sa vzťahuje na celý regresný model, nie na jednotlivé regresné koeficienty. Ich významnosť testujeme pomocou t-testov.

6. Stupne voľnosti F-testu – df

V prípade štandardného OLS modelu s absolútnym členom predstavuje df počet vysvetľujúcich premenných.

\[ df=k-1 \]

V našom modeli máme tri vysvetľujúce premenné a preto

\[ df=3 \]

Táto hodnota predstavuje počet stupňov voľnosti v čitateli F-štatistiky.

7. Logaritmická vierohodnosť – logLik

Logaritmická vierohodnosť vyjadruje, ako dobre odhadnutý pravdepodobnostný model zodpovedá pozorovaným údajom.

\[ \ell(\theta)=\ln L(\theta) \]

kde \(L(\theta)\) je vierohodnostná funkcia a \(\theta\) označuje parametre modelu.

Pri lineárnom regresnom modeli sa štandardne vychádza z predpokladu normálneho rozdelenia náhodných chýb.

Interpretácia: Pri porovnávaní modelov odhadnutých na rovnakých údajoch a pre tú istú vysvetľovanú premennú znamená vyššia hodnota maximálnej logaritmickej vierohodnosti lepšie prispôsobenie modelu údajom.

Samotná hodnota logLik však nepenalizuje počet parametrov modelu.

8. Akaikeho informačné kritérium – AIC

Akaikeho informačné kritérium má tvar:

\[ AIC=-2\ell(\hat\theta)+2K \]

kde:

  • \(\ell(\hat\theta)\) je maximálna hodnota logaritmickej vierohodnosti,
  • \(K\) je počet odhadovaných parametrov pravdepodobnostného modelu.

Pri normálnom lineárnom regresnom modeli sa do \(K\) započítavajú regresné koeficienty aj parameter rozptylu náhodnej zložky.

Interpretácia: Pri porovnávaní vhodných modelov uprednostňujeme model s nižšou hodnotou AIC.

AIC zohľadňuje kompromis medzi kvalitou prispôsobenia modelu údajom a jeho zložitosťou.

9. Bayesovské informačné kritérium – BIC

Bayesovské informačné kritérium má tvar:

\[ BIC=-2\ell(\hat\theta)+K\ln(n) \]

kde:

  • \(K\) je počet odhadovaných parametrov,
  • \(n\) je počet pozorovaní.

Interpretácia: Rovnako ako pri AIC uprednostňujeme nižšiu hodnotu BIC.

BIC spravidla penalizuje zložité modely prísnejšie ako AIC, pretože penalizačný člen závisí od počtu pozorovaní.

Poznámka: Hodnoty AIC a BIC majú význam predovšetkým pri porovnávaní modelov odhadnutých na porovnateľných údajoch a s rovnakou vysvetľovanou premennou.

10. Reziduálne stupne voľnosti – df.residual

Reziduálne stupne voľnosti vypočítame:

\[ df_{\mathrm{residual}}=n-k \]

kde:

  • \(n\) je počet pozorovaní,
  • \(k\) je počet odhadovaných regresných koeficientov vrátane absolútneho člena.

Napríklad pri \(n=100\) pozorovaniach a štyroch regresných koeficientoch:

\[ df_{\mathrm{residual}}=183-4=179 \]

Reziduálne stupne voľnosti využívame pri odhade rozptylu náhodnej zložky, pri t-testoch aj F-testoch.

11. Počet pozorovaní – N

Charakteristika N predstavuje počet pozorovaní skutočne použitých pri odhade regresného modelu.

V R zodpovedá:

nobs(model)
[1] 183

Pozor: Počet použitých pozorovaní môže byť nižší ako počet riadkov pôvodného dátového súboru, napríklad ak niektoré pozorovania obsahujú chýbajúce hodnoty.

12. Reziduálny súčet štvorcov – deviance

Táto charakteristika bola z našej tabuľky vynechaná, ale patrí medzi základné charakteristiky OLS regresie.

Pre model odhadnutý pomocou lm() platí:

\[ deviance=RSS=\sum_{i=1}^{n}(y_i-\hat y_i)^2 \]

alebo

\[ RSS=\sum_{i=1}^{n}e_i^2 \]

Interpretácia: Reziduálny súčet štvorcov meria celkovú veľkosť štvorcov odchýlok pozorovaných hodnôt od hodnôt predikovaných regresným modelom.

Metóda najmenších štvorcov (OLS) hľadá také odhady regresných koeficientov, ktoré minimalizujú práve túto veličinu.

Z reziduálneho súčtu štvorcov následne vypočítame reziduálnu štandardnú chybu:

\[ \hat\sigma=\sqrt{\frac{deviance}{df_{\mathrm{residual}}}} \]

V R si tento vzťah môžeme overiť:

deviance(model)
[1] 4534.949
sigma(model)
[1] 5.03338
sqrt(deviance(model) / df.residual(model))
[1] 5.03338

Posledné dva príkazy poskytujú rovnaký výsledok, s výnimkou prípadných numerických rozdielov spôsobených zaokrúhľovaním.

Záverečné zhrnutie

Pri interpretácii regresného modelu rozlišujeme tri základné skupiny charakteristík:

  1. Vysvetľovacia schopnosť modelu: R.squared, Adj.R.squared.
  2. Štatistická významnosť a presnosť modelu: Sigma, F Stat., p.value, df, df.residual.
  3. Porovnávanie modelov a ďalšie informácie: logLik, AIC, BIC, N, deviance.

Funkcia broom::glance() teda dopĺňa funkciu broom::tidy() o informácie potrebné na posúdenie regresného modelu ako celku.

Grafická analýza kvality modelu

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)

Diagnostické grafy regresného modelu
# (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.

Testy normality rezíduí

Keďže rezíduá \(\mathbf e\) sú lineárnou kombináciou náhodných ložiek \(mathbf u\), teda

\[\mathbf e = \mathbf y - \hat{\mathbf y} = \mathbf X \beta + \mathbf u - \mathbf X \hat{\beta} = \mathbf X \beta + \mathbf u - \mathbf X (\mathbf X^T \mathbf X)^{-1} \mathbf X^T \mathbf y = \] \[ = (\mathbf I - \mathbf X (\mathbf X^T \mathbf X)^{-1} \mathbf X^T\mathbf X) \mathbf u \]

a vieme, 6e linerárna kombinácia lineárne rozdelených veličím (pozrime sa na výraz dole) je tiež lineárne rozdelená veličina, potom predpoklad normálne rozdelených (nepozorovateľných) veličín \(\mathbf u\) môžeme testovať sprostredkovane cez normalitu pozorovatených veličín - rezíduí \(\mathbf e\).

Jarqueov–Berov test normality rezíduí

Jarqueov–Berov test (JB test) overuje, či rozdelenie rezíduí regresného modelu zodpovedá normálnemu rozdeleniu. Test vychádza z dvoch charakteristík rozdelenia: šikmosti (skewness) a špicatosti (kurtosis).

Pri normálnom rozdelení je šikmosť rovná nule a špicatosť má hodnotu 3.

Testujeme hypotézy:

  • \(H_0\): Rozdelenie rezíduí je normálne.
  • \(H_1\): Rozdelenie rezíduí nie je normálne.

Testovacia štatistika:

\[ JB = \frac{n}{6}\left(S^2+\frac{(K-3)^2}{4}\right) \]

kde:

  • \(n\) – počet pozorovaní,
  • \(S\) – výberová šikmosť rezíduí,
  • \(K\) – výberová špicatosť rezíduí.

Pri platnosti nulovej hypotézy má testovacia štatistika asymptoticky chí-kvadrát rozdelenie s dvoma stupňami voľnosti:

\[ JB \overset{a}{\sim} \chi^2(2) \]

Rozhodovacie pravidlo:

Ak je \(p < 0.05\), zamietame nulovú hypotézu o normalite rezíduí na 5 % hladine významnosti.

Ak je \(p \geq 0.05\), nulovú hypotézu nezamietame. To však neznamená, že sme normalitu rezíduí dokázali.

# Jarque-Bera normality test
res <- residuals(model)

jb_test <- tseries::jarque.bera.test(res)

broom::tidy(jb_test) |>
  dplyr::select(statistic, p.value, method) |>
  knitr::kable(
    digits = 4,
    col.names = c("JB statistic", "p-value", "Test"),
    caption = "Jarque-Bera test of residual normality"
  ) |>
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = FALSE
  )
Jarque-Bera test of residual normality
JB statistic p-value Test
21.3907 0 Jarque Bera Test

Poznámka: Existujú aj ďalšie testy normality, napríklad Shapirov–Wilkov test, Andersonov–Darlingov test alebo Kolmogorovov–Smirnovov test. Jarqueov–Berov test volíme najmä pre jeho jednoduchú interpretáciu a časté využitie v ekonometrii. Umožňuje identifikovať odchýlky od normality prostredníctvom šikmosti a špicatosti, ktoré študenti už poznajú z opisnej štatistiky. Keďže test využíva asymptotické rozdelenie a nie je citlivý na všetky typy odchýlok od normality, je vhodné jeho výsledky doplniť napríklad Q-Q grafom.

Test odľahlých pozorovaní – outlierTest()

Test odľahlých pozorovaní slúži na identifikáciu pozorovaní, ktorých hodnoty vysvetľovanej premennej sa výrazne odlišujú od hodnôt predikovaných regresným modelom.

Funkcia car::outlierTest() využíva externe študentizované rezíduá, ktoré sa vypočítajú ako:

\[ t_i = \frac{e_i}{\hat{\sigma}_{(i)}\sqrt{1-h_{ii}}} \]

kde:

  • \(e_i = y_i-\hat{y}_i\) – rezíduum \(i\)-teho pozorovania,
  • \(h_{ii}\) – diagonálny prvok Hat matice (leverage),
  • \(\hat{\sigma}_{(i)}\) – odhad smerodajnej odchýlky náhodnej zložky získaný po vynechaní \(i\)-teho pozorovania.

Testované hypotézy:

  • \(H_0\): Pozorovanie nie je odľahlé vzhľadom na špecifikovaný regresný model.
  • \(H_1\): Pozorovanie je odľahlé vzhľadom na špecifikovaný regresný model.

Za predpokladu normálne rozdelených nezávislých chýb má externe študentizované rezíduum pri nulovej hypotéze Studentovo t-rozdelenie s \(n-k-1\) stupňami voľnosti, kde \(k\) je počet regresných koeficientov vrátane absolútneho člena.

Keďže testujeme viacero pozorovaní súčasne, funkcia používa Bonferroniho korekciu p-hodnôt:

\[ p_{i,\mathrm{Bonf}}=\min(1,n\cdot p_i) \]

kde \(p_i\) predstavuje neupravenú p-hodnotu individuálneho testu a \(n\) počet testovaných pozorovaní.

Rozhodovacie pravidlo:

Ak \(p_{i,\mathrm{Bonf}}<0.05\), príslušné pozorovanie identifikujeme ako štatisticky významne odľahlé na hladine významnosti 5 %.

# Outlier test with Bonferroni correction


outlier_test <- car::outlierTest(model)

if (is.null(outlier_test)) {
  cat("No observations returned by outlierTest.")
} else {

  outlier_table <- data.frame(
    Observation = names(outlier_test$rstud),
    rstudent = unname(outlier_test$rstud),
    unadjusted.p = unname(outlier_test$p),
    bonf.p = unname(outlier_test$bonf.p)
  )

  outlier_table |>
    knitr::kable(
      digits = 4,
      col.names = c(
        "Observation",
        "Studentized residual",
        "Unadjusted p-value",
        "Bonferroni p-value"
      ),
      caption = "Outlier test with Bonferroni correction"
    ) |>
    kableExtra::kable_styling(
      bootstrap_options = c("striped", "hover"),
      full_width = FALSE
    )
}
Outlier test with Bonferroni correction
Observation Studentized residual Unadjusted p-value Bonferroni p-value
433 -3.643 4e-04 0.0647

Poznámka: Odľahlé pozorovanie nemusí byť chybné. Môže obsahovať dôležitú ekonomickú informáciu, a preto by sme ho nemali automaticky odstrániť. Test identifikuje predovšetkým nezvyčajné hodnoty vysvetľovanej premennej vzhľadom na regresný model. Na posúdenie vplyvu pozorovania na odhad regresných koeficientov môžeme použiť napríklad Cookovu vzdialenosť.

3. Interpretácia jednotlivých stĺpcov

Stĺpec Význam
Observation Identifikátor pozorovania, pri ktorom sa skúma možná odľahlosť.
Studentized residual Externe študentizované rezíduum, pri ktorom sa rozptyl odhaduje bez príslušného pozorovania.
Unadjusted p-value P-hodnota individuálneho testu bez korekcie na viacnásobné testovanie.
Bonferroni p-value P-hodnota upravená Bonferroniho korekciou, ktorá zohľadňuje viacnásobné testovanie odľahlých pozorovaní.

Dôležité: Pri identifikácii odľahlých pozorovaní by sme sa mali rozhodovať podľa Bonferroni p-value. Ak je menšia ako 0,05, príslušné pozorovanie považujeme za štatisticky významne odľahlé na 5 % hladine významnosti.

Štatisticky významne odľahlé pozorovanie však nemusí byť chybné a nemali by sme ho automaticky odstrániť z dátového súboru.

Poznámka: Funkcia car::outlierTest() štandardne zobrazuje len najextrémnejšie pozorovania, nie všetky riadky dátového súboru. Ak nenájde Bonferroniho významné pozorovanie, môže napriek tomu zobraziť pozorovanie s najväčším absolútnym študentizovaným rezíduom.

V našom prípade je p-hodnota testu 433. pozorovania mierne nad 0.06, čo nepovažujeme za vážny problém. Napriek tomu sa pokúsime tento problém redukovať.

Úprava modelu

Problém nelinearity, alebo odľahlých hodnôt sa dá efektívne retukovať s logaritmickými transformáciami jednotlivých premenných.

\[ Life\_Expectancy_t = \beta_0 + \beta_1 \ln(GDP_t) + \beta_2 Schooling_t + u_t\] Logaritmicá transformácia zároveň znamená zmenu interpretácie zodpovedajúceho regresného koeficientu. V našom prípade by sa \(\beta_1\) interpretovalo nasledovne. Ak \(GDP_t\) sa zmení o 1 percento, potom očakávaná dĺžka života sa zmení o \(\beta_1\) rokov.

library(broom)
library(knitr)
library(kableExtra)

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

tidy(model2) |>
  kable(
    digits = 4,
    col.names = c(
      "Variable",
      "Estimate",
      "Std. Error",
      "t value",
      "Pr(>|t|)"
    ),
    caption = "Regression results – Life expectancy"
  ) |>
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = FALSE
  )
Regression results – Life expectancy
Variable Estimate Std. Error t value Pr(>&#124;t&#124;)
(Intercept) 40.0434 2.1992 18.2079 0.0000
log(GDP) 0.6338 0.3002 2.1115 0.0361
Schooling 2.0561 0.1557 13.2028 0.0000

Graphical diagnostic of the model

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

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

Diagnostické grafy regresného modelu
# (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))

Test nornality rezíduí upraveného modelu

# Jarque-Bera normality test
res <- residuals(model2)

jb_test <- tseries::jarque.bera.test(res)

broom::tidy(jb_test) |>
  dplyr::select(statistic, p.value, method) |>
  knitr::kable(
    digits = 4,
    col.names = c("JB statistic", "p-value", "Test"),
    caption = "Jarque-Bera test of residual normality"
  ) |>
  kableExtra::kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = FALSE
  )
Jarque-Bera test of residual normality
JB statistic p-value Test
30.1075 0 Jarque Bera Test

Logaritmická transformácia teda nevedie k priblíženiu rozdelenia rezíduí k normalite.

Test odľahlých hodnôt upraveného modelu

# Outlier test with Bonferroni correction


outlier_test <- car::outlierTest(model2)

if (is.null(outlier_test)) {
  cat("No observations returned by outlierTest.")
} else {

  outlier_table <- data.frame(
    Observation = names(outlier_test$rstud),
    rstudent = unname(outlier_test$rstud),
    unadjusted.p = unname(outlier_test$p),
    bonf.p = unname(outlier_test$bonf.p)
  )

  outlier_table |>
    knitr::kable(
      digits = 4,
      col.names = c(
        "Observation",
        "Studentized residual",
        "Unadjusted p-value",
        "Bonferroni p-value"
      ),
      caption = "Outlier test with Bonferroni correction"
    ) |>
    kableExtra::kable_styling(
      bootstrap_options = c("striped", "hover"),
      full_width = FALSE
    )
}
Outlier test with Bonferroni correction
Observation Studentized residual Unadjusted p-value Bonferroni p-value
433 -3.8286 2e-04 0.0326

Záver

V tejto úlohe sme kvantifikovali viacnásobnú lineárnu regresiu a demonštrovali sme základnú diagnostiku modelu. Logaritmická transformácia spôsobila vznik odľahlej hodnoty a teda nezlepšila tech nické charakteristiky modelu. Mali by sme sa vrátiť k pôvodnému lineárnemu modelu s vynechaním neinterpretovateľnej vysvetľujúcej premennej \(BMI\).

save.image("pamat.Rdata")

rm(list=ls())