Econometrics in R - cvičenie 5

Základná diagnostika modelu

Author

Aleksei Solberg

Pokračujeme s využívaním World Bank World Development Indicators databázy.

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

Príprava databázy

cesta <- "world_bank_development_indicators.csv"
if (!file.exists(cesta)) cesta <- file.path("udaje", "world_bank_development_indicators.csv")

udaje <- read.csv(cesta, dec=".", sep=",", header = TRUE, na.strings = c("", "NA"), check.names = FALSE, stringsAsFactors = FALSE)
names(udaje) <- make.names(trimws(names(udaje)), unique = TRUE)

# select just the record from 2013
rok <- substr(as.character(udaje$date), 1, 4)
udaje.2013 <- udaje[!is.na(rok) & rok == "2013", c("life_expectancy_at_birth", "individuals_using_internet.", "government_health_expenditure.", "GDP_current_US")]
names(udaje.2013) <- c("life_expectancy", "internet", "health_exp", "gdp")
udaje.2013 <- udaje.2013[!is.na(udaje.2013$life_expectancy), ]

# data imputation

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

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

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

udaje.2013 <- udaje_imputed
udaje.2013$gdp_mld <- udaje.2013$gdp / 1e9

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

Odožím si celý obsah - všetky objekty - pamäte (workspace) do súboru

save.image("pamat.Rdata")

Lineárna regresia

Odhadujeme regresný model:

\[ LifeExpectancy_i = \beta_0+\beta_1 Internet_i + \beta_2 HealthExp_i + \beta_3 GDP\_mld_i+u_i \]

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

model <- lm(life_expectancy ~ +1 + internet + health_exp + gdp_mld, data=udaje.2013)

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)	61.0602	0.5368	113.7585	0.0000	***
internet	0.2233	0.0132	16.9074	0.0000	***
health_exp	0.3591	0.1766	2.0330	0.0431	*
gdp_mld	0.0000	0.0000	0.7568	0.4499

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.6915	0.6878	4.4681	189.7503	0	3	-750.2873	1510.575	1528.339	254	258

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] 258

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=\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] 5070.846

sigma(model)

[1] 4.468105

sqrt(deviance(model) / df.residual(model))

[1] 4.468105

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:

Vysvetľovacia schopnosť modelu: R.squared, Adj.R.squared.
Štatistická významnosť a presnosť modelu: Sigma, F Stat., p.value, df, df.residual.
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í. Červená LOWESS (Locally Weighted Scatterplot Smooting) krivka vzniká vyhladzovaním približne 2/3 najbližších pozorovaní pomocou váženej metódy najmenších štvorcov (Weighted Least Squares), kde najvyššie váhy majú najbližšie pozorovania k referenčnému pozorovaniu.

# 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 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á, to naznačuje možnú miernu nelinearitu – modelu môže chýbať nelineárny tvar nejakej premennej (napr. HDP).
Rozptyl rezíduí - Vertikálny rozptyl (variancia rezíduí) sa javí ako približne konštantný v rámci prispôsobených hodnôt, prípadne s mierne vyšším rozptylom pri nižších hodnotách.
Ak by sa rezíduá výrazne rozprestierali (tvorili kónus), naznačovalo by to heteroskedasticitu.
Odľahlé hodnoty - Niekoľko bodov 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

Interpretácia grafu

Celkový tvar

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

Krajné hodnoty (extrémy)

Body na koncoch (najmä vľavo dole) sa 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 overujeme hypotézu o normálnom rozložení modelových chýb, ako aj prítomnosť extrémnych hodnôt.

Scale location plot

Interpretácia grafu

Horizontálne rozptýlenie

Body sú rozptýlené po osi x bez výrazného vytvorenia lievika. To naznačuje približne konštantnú varianciu.
Červená hladká čiara je pomerne plochá – znak, že pri zvyšovaní vyrovnaných hodnôt nedochádza k dramatickej systematickej zmene variancie.
Odľahlé hodnoty - niekoľko bodov je mierne nad 1,5, čo zodpovedá krajinám s väčšou odchýlkou.

residuals vs leverage

Interpretácia grafu

Rozloženie vplyvu

Väčšina pozorovaní má nízky vplyv (páku) — typické pre veľké vzorky.

Niekoľko bodov (krajiny s extrémne vysokým celkovým HDP, ako USA alebo Čína) vyniká vysokou pákou, čo znamená, že ich hodnoty \(X\) sú ďaleko od väčšiny údajov.

Veľkosť rezíduí

Štandardizované rezíduá sú väčšinou medzi −2 a +2 – to je dobré.

Kontúry Cookovej vzdialenosti

Pozorovania s vysokou pákou môžu mať vyšší vplyv na sklon regresnej nadroviny pri premennej HDP, čo je ďalší dôvod na logaritmickú transformáciu HDP.

Jarqueov–Berov test normality rezíduí

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
20.9787	0	Jarque Bera Test

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

# 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
9010	-3.9697	1e-04	0.0242

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

model2 <- lm(
  life_expectancy ~ log(gdp_mld) + internet,
  data = udaje.2013
)

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(>\|t\|)
(Intercept)	61.0295	0.5925	102.9953	0.0000
log(gdp_mld)	0.1522	0.0991	1.5361	0.1258
internet	0.2377	0.0106	22.4610	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)

# (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 normality 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
13.7754	0.001	Jarque Bera Test

Na základe p-hodnoty Jarqueovho-Berovho testu posudzujeme, či logaritmická transformácia viedla k zlepšeniu normality rezíduí.

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
9010	-3.604	4e-04	0.0973

Malé úlohy

Porovnajte oba vyššie uvedené modely a ktorý možno považovať za lepší pri použití Akaike Information Criterion?

model3 <- lm(
  life_expectancy ~ log(gdp_mld) + internet + health_exp,
  data = udaje.2013
)

porovnanie_aic <- bind_rows(
  glance(model) |> mutate(Model = "model (lineárny)"),
  glance(model2) |> mutate(Model = "model2 (log GDP + internet)"),
  glance(model3) |> mutate(Model = "model3 (log GDP + internet + health_exp)")
) |>
  select(Model, r.squared, adj.r.squared, sigma, AIC, BIC)

porovnanie_aic |>
  kable(
    digits = 4,
    caption = "Porovnanie modelov pomocou informačných kritérií"
  ) |>
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = FALSE
  )

Porovnanie modelov pomocou informačných kritérií
Model	r.squared	adj.r.squared	sigma	AIC	BIC
model (lineárny)	0.6915	0.6878	4.4681	1510.575	1528.339
model2 (log GDP + internet)	0.6874	0.6849	4.4889	1511.985	1526.197
model3 (log GDP + internet + health_exp)	0.6930	0.6894	4.4570	1509.296	1527.060

Pri porovnaní prvých dvoch modelov má model hodnotu \(AIC = 1510.57\) a model2 hodnotu \(AIC = 1511.98\). Keďže pri Akaikeho informačnom kritériu uprednostňujeme model s nižšou hodnotou AIC, za lepší považujeme model s nižšou hodnotou v tabuľke vyššie.

Urobte inú modifikáciu modelu a výsledky porovnajte.

Ako ďalšiu modifikáciu sme odhadli model3, v ktorom sme ponechali logaritmus HDP (log(gdp_mld)), penetráciu internetu (internet) a zároveň sme vrátili vládne výdavky na zdravotníctvo (health_exp). Výsledky (vrátane \(R^2\), upraveného \(\bar R^2\), AIC a BIC) sú zobrazené v súhrnnej tabuľke vyššie.

Záver

V tejto úlohe sme kvantifikovali viacnásobnú lineárnu regresiu na údajoch Svetovej banky pre rok 2013 a demonštrovali sme základnú diagnostiku modelu. Overili sme štatistickú významnosť jednotlivých premenných, celkovú kvalitu modelu pomocou F-testu a koeficientu determinácie, ako aj predpoklady o normalite rezíduí a prítomnosti odľahlých pozorovaní. Následne sme porovnali pôvodnú špecifikáciu s logaritmicky transformovanými modelmi pomocou informačných kritérií AIC a BIC.

save.image("pamat.Rdata")

rm(list=ls())

Príprava databázy

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

Lineárna regresia

Kvalita vyrovnania modelu

11. Počet pozorovaní – N

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

Záverečné zhrnutie

Grafická analýza kvality modelu

Residuals vs. fitted

Interpretácia grafu

Q-Q plot

Interpretácia grafu

Scale location plot

Interpretácia grafu

residuals vs leverage

Interpretácia grafu

Jarqueov–Berov test normality rezíduí

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

Graphical diagnostic of the model

Test normality rezíduí upraveného modelu

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

Malé úlohy

Záver

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