Diagnostika modelu, robustné štandardné chyby a GLS odhad
Livia
1 Úvod
Tento notebook vysvetľuje tri nadväzujúce oblasti diagnostiky a nápravy lineárneho regresného modelu:
- Breusch–Paganov test pre heteroskedasticitu,
- Breusch–Godfreyho test pre autokoreláciu,
- robustné štandardné chyby a GLS odhad ako nápravné metódy.
Pedagogická logika notebooku je zámerne nasledovná:
- uvedieme základný teoretický koncept,
- odhadneme model metódou OLS,
- overíme predpoklady modelu,
- ukážeme, ako opraviť inferenciu pomocou robustných štandardných chýb,
- napokon ukážeme, ako zmeniť samotný odhad modelu pomocou GLS.
2 Teoretické východisko
Uvažujme lineárny regresný model v maticovom tvare:
\[ y = X\beta + u \]
Pri klasických predpokladoch platí:
\[ E(u) = 0, \qquad \mathrm{Var}(u) = \sigma^2 I \]
Za týchto predpokladov je OLS odhad daný vzťahom:
\[ \hat{\beta}_{OLS} = (X'X)^{-1}X'y \]
a jeho kovariančná matica je:
\[ \mathrm{Var}(\hat{\beta}_{OLS}) = \sigma^2 (X'X)^{-1} \]
Ak však neplatí \(\mathrm{Var}(u) = \sigma^2 I\), ale všeobecnejšie \(\mathrm{Var}(u) = \Omega\), kde \(\Omega\) nie je násobkom jednotkovej matice, potom:
- OLS odhady môžu zostať nestranné,
- ale už nemusia byť efektívne,
- a klasické štandardné chyby nemusia byť správne.
Preto existujú dve základné stratégie:
- ponechať OLS koeficienty a opraviť iba štandardné chyby,
- zmeniť samotný spôsob odhadu modelu pomocou GLS.
3 Breusch–Paganov test (o heteroskedasticite)
3.1 Myšlienka testu
Uvažujme lineárny regresný model:
\[ y_i = \beta_0 + \beta_1 x_{1i} + \cdots + \beta_k x_{ki} + u_i \]
Za klasického predpokladu homoskedasticity platí \(\mathrm{Var}(u_i) = \sigma^2\), čo znamená, že rozptyl náhodnej zložky je konštantný pre všetky pozorovania. Breusch–Paganov test skúma, či sa tento rozptyl nemení v závislosti od vysvetľujúcich premenných.
3.2 Hypotézy
\[ H_0: \mathrm{Var}(u_i) = \sigma^2 \]
\[ H_1: \mathrm{Var}(u_i) \neq \sigma^2 \]
3.3 Postup testu
Najprv odhadneme pôvodný model metódou OLS a získame rezíduá \(e_i = y_i - \hat{y}_i\). Následne odhadneme pomocnú regresiu:
\[ e_i^2 = \alpha_0 + \alpha_1 z_{1i} + \cdots + \alpha_m z_{mi} + v_i \]
3.4 Testovacia štatistika
\[ LM = nR^2 \sim \chi^2(m) \]
kde \(R^2\) pochádza z pomocnej regresie a \(m\) je počet regresorov bez interceptu.
3.5 Rozhodovanie
- veľká p-hodnota \(\Rightarrow\) nulovú hypotézu (o homoskedasticite) nezamietame,
- malá p-hodnota \(\Rightarrow\) zamietame nulovú hypotézu a usudzujeme na heteroskedasticitu.
4 Breusch–Godfreyho test (o autokorelácii)
4.1 Myšlienka testu
Predpokladajme, že chyba môže mať autoregresnú štruktúru rádu \(p\):
\[ u_t = \rho_1 u_{t-1} + \rho_2 u_{t-2} + \cdots + \rho_p u_{t-p} + \varepsilon_t \]
4.2 Hypotézy
\[ H_0: \rho_1 = \rho_2 = \cdots = \rho_p = 0 \]
\[ H_1: \text{aspoň jedno } \rho_j \neq 0 \]
4.3 Testovacia štatistika
\[ LM = nR^2 \sim \chi^2(p) \]
4.4 Rozhodovanie
- veľká p-hodnota \(\Rightarrow\) nulovú hypotézu o neprítomnosti autokorelácie nezamietame,
- malá p-hodnota \(\Rightarrow\) zamietame nulovú hypotézu a usudzujeme na autokoreláciu.
5 Robustné štandardné chyby
5.1 Základná idea
Ak sú OLS koeficienty stále konzistentné, ale problém je v nesprávnej kovariančnej matici, môžeme ponechať odhad \(\hat{\beta}_{OLS}\) a korigovať iba štandardné chyby. Všeobecný robustný tvar kovariančnej matice je:
\[ \widehat{\mathrm{Var}}(\hat{\beta}) = (X'X)^{-1} X' \hat{\Omega} X (X'X)^{-1} \]
kde \(\hat{\Omega}\) zachytáva heteroskedasticitu alebo aj autokoreláciu.
Pri zanedbaní narušenia predpokladov by sme túto maticu odhadovali vzorcom \(\widehat{\mathrm{Var}}(\hat{\beta}) = \sigma^2(X'X)^{-1}\) a dopustili by sme sa chyby v odhade štandardných chýb, čo by mohlo viesť k nesprávnym záverom o štatistickej významnosti koeficientov.
5.2 Whiteove robustné štandardné chyby
Whiteove odhady sú robustné voči heteroskedasticite. Koeficienty regresie sa nemenia, menia sa iba štandardné chyby, t-štatistiky a p-hodnoty.
5.3 HAC odhady
Ak chceme korigovať štandardné chyby súčasne na heteroskedasticitu aj autokoreláciu, používame HAC odhady, napríklad Newey–West.
6 GLS odhad
6.1 Teoretický základ
Cieľom GLS je transformovať pôvodný model tak, aby náhodná zložka v transformovanom modeli spĺňala klasické predpoklady. Pre maticu \(P\), kde \(P'P = \Omega^{-1}\), dostávame transformovaný model \(y^* = X^*\beta + u^*\), pre ktorý platí \(\mathrm{Var}(u^*) = I\). Výsledný GLS odhad má tvar:
\[ \hat{\beta}_{GLS} = (X'\Omega^{-1}X)^{-1}X'\Omega^{-1}y \]
GLS „váži” jednotlivé pozorovania podľa variability chýb —
pozorovania s vyšším rozptylom dostávajú menšiu váhu. V praxi sa matica
\(\Omega\) odhaduje (FGLS), napríklad
pomocou gls() z balíka nlme.
6.2 OLS s robustnými SE verzus GLS
- OLS + White/HAC SE: koeficienty sa nemenia, opravuje sa iba inferencia.
- GLS: menia sa aj samotné odhady koeficientov, pretože modelujeme štruktúru chýb.
7 Použité balíky
library(lmtest)
library(sandwich)
library(nlme)8 Načítanie dát a odhad modelu
Uvažujeme lineárny regresný model, v ktorom vysvetľujeme trhovú cenu hráča (v mil. €) pomocou jeho veku a výšky:
\[ price_i = \beta_0 + \beta_1 \cdot age_i + \beta_2 \cdot height_i + \varepsilon_i \]
data <- read.csv("top5_leagues_player.csv")
data <- data[complete.cases(data[, c("price", "age", "height")]), ]
data$t <- 1:nrow(data)
m_ols <- lm(price ~ age + height, data = data)
summary(m_ols)9 Príklad 1: Heteroskedasticita
9.1 Breusch–Paganov test
bp_res <- bptest(m_ols)
bp_res##
## studentized Breusch-Pagan test
##
## data: m_ols
## BP = 21.833, df = 2, p-value = 1.816e-059.1.1 Interpretácia
Ak je p-hodnota testu menšia než zvolená hladina významnosti (0,05), zamietame nulovú hypotézu homoskedasticity. To znamená, že rozptyl rezíduí pravdepodobne závisí od vysvetľujúcich premenných a v modeli je prítomná heteroskedasticita.
9.2 Pomocná regresia
##
## Call:
## lm(formula = ehat2 ~ age + height, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -449.4 -232.3 -177.7 -92.0 27822.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1189.669 580.673 2.049 0.0406 *
## age -21.002 4.568 -4.597 4.49e-06 ***
## height -213.851 313.563 -0.682 0.4953
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1073 on 2588 degrees of freedom
## Multiple R-squared: 0.008426, Adjusted R-squared: 0.00766
## F-statistic: 11 on 2 and 2588 DF, p-value: 1.756e-05n * summary(aux_bp)$r.squared## [1] 21.832659.2.1 Interpretácia
Ak pomocná regresia vysvetľuje podstatnú časť variability \(e_i^2\), znamená to, že rozptyl chýb nie je konštantný. Vysvetľujúce premenné (vek a výška) teda pomáhajú predpovedať veľkosť chýb.
9.3 Vizualizácia heteroskedasticity
Ak sa rozptyl štvorcov rezíduí mení so zmenou vyrovnaných hodnôt (napr. rastie pri drahších hráčoch), ide o indikáciu heteroskedasticity. Červená LOESS krivka pomáha identifikovať systematický vzťah medzi predikovanou cenou a variabilitou chýb.
9.4 Whiteove robustné štandardné chyby
white_hc0 <- coeftest(m_ols, vcov = vcovHC(m_ols, type = "HC0"))
white_hc1 <- coeftest(m_ols, vcov = vcovHC(m_ols, type = "HC1"))
white_hc0##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.31313 8.49704 5.0974 3.691e-07 ***
## age -0.42759 0.05673 -7.5373 6.593e-14 ***
## height -11.66483 4.58475 -2.5443 0.01101 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1white_hc1##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.313134 8.501964 5.0945 3.749e-07 ***
## age -0.427590 0.056763 -7.5330 6.813e-14 ***
## height -11.664826 4.587405 -2.5428 0.01105 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 19.4.1 Interpretácia
Koeficienty regresie zostávajú rovnaké ako pri OLS, ale štandardné chyby sa upravia tak, aby boli robustné voči heteroskedasticite. Ak koeficient zostáva významný aj po tejto úprave, jeho významnosť je robustnejšia.
9.5 GLS pri heteroskedasticite
gls_het <- gls(
price ~ age + height,
data = data,
weights = varConstPower(form = ~ age)
)
summary(gls_het)## Generalized least squares fit by REML
## Model: price ~ age + height
## Data: data
## AIC BIC logLik
## 21276.07 21311.23 -10632.04
##
## Variance function:
## Structure: Constant plus power of variance covariate
## Formula: ~age
## Parameter estimates:
## const power
## 6.812175e-11 -1.913112e+00
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 47.63640 7.044516 6.762197 0.0000
## age -0.80946 0.054828 -14.763501 0.0000
## height -8.15670 3.811179 -2.140205 0.0324
##
## Correlation:
## (Intr) age
## age -0.147
## height -0.973 -0.083
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -0.8426302 -0.5837932 -0.3753471 0.1436377 9.7667987
##
## Residual standard error: 7441.558
## Degrees of freedom: 2591 total; 2588 residual9.5.1 Interpretácia
Na rozdiel od Whiteových štandardných chýb teraz nemeníme iba inferenciu, ale aj samotný odhad modelu. GLS priamo zohľadňuje, že rozptyl chýb sa mení s predikovanou hodnotou trhovej ceny hráčov, a preto môžu byť odhady efektívnejšie.
10 Príklad 2: Autokorelácia
10.1 Breusch–Godfreyho test
bg_res_1 <- bgtest(m_ols, order = 1)
bg_res_1##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: m_ols
## LM test = 343.97, df = 1, p-value < 2.2e-16bg_res_2 <- bgtest(m_ols, order = 2)
bg_res_2##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: m_ols
## LM test = 465.6, df = 2, p-value < 2.2e-1610.1.1 Interpretácia
Ak je p-hodnota malá, zamietame nulovú hypotézu neprítomnosti autokorelácie. To znamená, že súčasné rezíduum súvisí s minulými rezíduami a klasické OLS štandardné chyby nemusia byť spoľahlivé. Keďže dáta nepredstavujú časový rad, prípadná autokorelácia môže byť spôsobená systematickým usporiadaním pozorovaní v datasete (napr. podľa klubu alebo ligy).
10.2 Pomocná regresia pre BG test
nrow(aux_bg_data) * summary(aux_bg)$r.squared## [1] 344.499710.2.1 Interpretácia
Ak je koeficient pri oneskorenom rezíduu (ehat_lag1)
štatisticky významný, znamená to, že minulé chyby pomáhajú vysvetliť
súčasné chyby. To je typický znak autokorelácie.
10.3 HAC štandardné chyby
hac_nw <- coeftest(m_ols, vcov = NeweyWest(m_ols))
hac_nw##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.313134 9.585999 4.5184 6.511e-06 ***
## age -0.427590 0.074287 -5.7559 9.631e-09 ***
## height -11.664826 4.861362 -2.3995 0.01649 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 110.3.1 Interpretácia
Koeficienty zostávajú rovnaké ako pri OLS, ale štandardné chyby sú upravené tak, aby boli robustné voči heteroskedasticite aj autokorelácii. Výsledná inferencia je preto dôveryhodnejšia.
10.4 GLS pri autokorelácii
gls_ar1 <- gls(
price ~ age + height,
data = data,
correlation = corAR1(form = ~ t)
)
summary(gls_ar1)## Generalized least squares fit by REML
## Model: price ~ age + height
## Data: data
## AIC BIC logLik
## 21234.64 21263.94 -10612.32
##
## Correlation Structure: AR(1)
## Formula: ~t
## Parameter estimate(s):
## Phi
## 0.3721605
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 19.759992 8.484911 2.328839 0.0199
## age -0.395535 0.058705 -6.737685 0.0000
## height 0.742287 4.551958 0.163070 0.8705
##
## Correlation:
## (Intr) age
## age -0.180
## height -0.982 -0.002
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -0.9110579 -0.5204650 -0.3469207 0.1495574 10.7420537
##
## Residual standard error: 15.67778
## Degrees of freedom: 2591 total; 2588 residual10.4.1 Interpretácia
GLS teraz priamo modeluje to, že chyby v datasete spolu súvisia (napr. hráči z rovnakého klubu). Preto sa môžu zmeniť nielen štandardné chyby, ale aj samotné odhadnuté koeficienty.
11 Porovnanie prístupov
11.1 OLS verzus robustné štandardné chyby verzus GLS
cat("=== OLS koeficienty ===\n")## === OLS koeficienty ===coef(m_ols)## (Intercept) age height
## 43.31313 -0.42759 -11.66483cat("\n=== GLS koeficienty (heteroskedasticita) ===\n")##
## === GLS koeficienty (heteroskedasticita) ===coef(gls_het)## (Intercept) age height
## 47.636403 -0.809458 -8.156704cat("\n=== GLS koeficienty (autokorelácia AR1) ===\n")##
## === GLS koeficienty (autokorelácia AR1) ===coef(gls_ar1)## (Intercept) age height
## 19.7599925 -0.3955345 0.742287511.1.1 Interpretácia
Pri OLS a OLS s robustnými štandardnými chybami ostávajú odhadnuté koeficienty rovnaké. Pri GLS sa však môžu zmeniť aj samotné koeficienty, pretože GLS využíva dodatočné informácie o štruktúre chýb.
12 Zhrnutie
12.1 Kedy použiť jednotlivé prístupy?
- OLS — základný odhad modelu.
- Breusch–Paganov test — test heteroskedasticity.
- Breusch–Godfreyho test — test autokorelácie.
- Whiteove štandardné chyby — OLS koeficienty zachované, opravená inferencia pri heteroskedasticite.
- HAC štandardné chyby — OLS koeficienty zachované, opravená inferencia pri heteroskedasticite aj autokorelácii.
- GLS — priame modelovanie štruktúry chýb, efektívnejší odhad modelu.
13 Záver
Diagnostika regresného modelu nekončí odhadom koeficientov. Po odhade je potrebné overiť, či sú splnené predpoklady o štruktúre chýb, a ak nie sú, zvoliť vhodnú nápravu. V tomto notebooku sme ukázali tri úrovne práce s problémami rezíduí: (1) diagnostika pomocou BP a BG testu, (2) oprava inferencie pomocou Whiteových a HAC štandardných chýb, (3) zmena samotného odhadu pomocou GLS.
V kontexte modelu trhovej ceny futbalových hráčov je porušenie homoskedasticity prirodzené — pri lacnejších hráčoch sú ceny zoskupené tesnejšie, zatiaľ čo pri drahších hráčoch sa variabilita výrazne zvyšuje. Robustné metódy a GLS umožňujú získať spoľahlivejšie závery napriek tomuto porušeniu predpokladov. ```