Tento notebook vysvetľuje tri nadväzujúce oblasti diagnostiky a nápravy lineárneho regresného modelu:
Pedagogická logika notebooku je zámerne nasledovná:
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:
Preto existujú dve základné stratégie:
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 \]
To 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.
\[ H_0: \mathrm{Var}(u_i) = \sigma^2 \]
\[ H_1: \mathrm{Var}(u_i) \neq \sigma^2 \]
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} + \alpha_2 z_{2i} + \cdots + \alpha_m z_{mi} + v_i \]
V praxi sa často volí:
\[ z_{ji} = x_{ji} \]
Breusch–Paganov test je založený na LM štatistike:
\[ LM = nR^2 \]
kde \(R^2\) pochádza z pomocnej regresie. Za platnosti nulovej hypotézy platí:
\[ LM \sim \chi^2(m) \]
kde \(m\) je počet regresorov v pomocnej regresii bez interceptu.
Uvažujme model:
\[ y_t = \beta_0 + \beta_1 x_{1t} + \cdots + \beta_k x_{kt} + u_t \]
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 \]
\[ H_0: \rho_1 = \rho_2 = \cdots = \rho_p = 0 \]
\[ H_1: \text{aspoň jedno } \rho_j \neq 0 \]
Najprv odhadneme pôvodný model a získame rezíduá \(e_t\). Potom odhadneme pomocnú regresiu:
\[ e_t = \gamma_0 + \gamma_1 x_{1t} + \cdots + \gamma_k x_{kt} + \rho_1 e_{t-1} + \cdots + \rho_p e_{t-p} + v_t \] kde \(v_t\) je náhodná zložka pomocnej regresie.
Aj tu používame LM princíp:
\[ LM = nR^2 \]
kde \(n\) je počet pozorovaní v predchádzajúcej testovanej rovnici (pričom tu dochádza k redukcii o \(p\)), pričom za nulovej hypotézy platí:
\[ LM \sim \chi^2(p) \]
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.
Pripomínam, že pri zanedbaní narušenia vyššie uvedených 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.
Whiteove odhady sú robustné voči heteroskedasticite. Koeficienty regresie sa nemenia, menia sa teda len:
Ak chceme korigovať štandardné chyby súčasne na heteroskedasticitu aj autokoreláciu, používame HAC odhady, napríklad Newey–West.
Uvažujme lineárny regresný model v maticovom tvare:
\[ y = X\beta + u \]
Pri klasických predpokladoch metódy najmenších štvorcov (OLS) platí:
\[ \mathrm{Var}(u) = \sigma^2 I \]
To znamená, že:
Ak je však tento predpoklad porušený, teda:
\[ \mathrm{Var}(u) = \Omega \]
kde \(\Omega \neq \sigma^2 I\), potom:
Cieľom GLS je transformovať pôvodný model tak, aby náhodná zložka v transformovanom modeli spĺňala klasické predpoklady.
Uvažujme maticu \(P\), pre ktorú platí:
\[ P'P = \Omega^{-1} \]
Následne vynásobíme pôvodný model zľava touto maticou:
\[ Py = PX\beta + Pu \]
Označme:
\[ y^* = Py, \quad X^* = PX, \quad u^* = Pu \]
Dostávame transformovaný model:
\[ y^* = X^*\beta + u^* \]
Pre náhodnú zložku transformovaného modelu platí:
\[ \mathrm{Var}(u^*) = P \Omega P' \]
Keďže platí \(P'P = \Omega^{-1}\), dostávame:
\[ \mathrm{Var}(u^*) = I \]
Teda:
a preto spĺňajú klasické predpoklady metódy najmenších štvorcov.
Na transformovaný model môžeme aplikovať OLS. Výsledný GLS odhad má tvar:
\[ \hat{\beta}_{GLS} = (X'\Omega^{-1}X)^{-1}X'\Omega^{-1}y \]
V praxi matica \(\Omega\) zvyčajne nie je známa, preto sa používa:
gls() z balíka
nlmeInterpretácia :
GLS nemení iba štandardné chyby ako robustné metódy, ale mení samotný odhad modelu.
Využíva informáciu o štruktúre chýb, čím dosahuje efektívnejšie odhady než OLS. ## Intuícia GLS
GLS možno chápať ako transformáciu modelu tak, aby v transformovanom modeli mala chyba „sférický“ tvar. Inými slovami, problémová kovariančná štruktúra sa zohľadní už pri odhade, nie až dodatočne pri výpočte štandardných chýb.
library(lmtest)
library(sandwich)
library(nlme)
Vytvoríme dáta s heteroskedasticitou, kde rozptyl chyby rastie s premennou \(x\).
set.seed(123)
n <- 200
x <- runif(n, 1, 10)
u <- rnorm(n, mean = 0, sd = 0.5 * x)
y <- 2 + 3 * x + u
data_bp <- data.frame(y, x)
m_bp <- lm(y ~ x, data = data_bp)
summary(m_bp)
##
## Call:
## lm(formula = y ~ x, data = data_bp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.3828 -1.5748 -0.1864 1.3994 9.6110
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.16353 0.46233 4.68 5.32e-06 ***
## x 2.96535 0.07607 38.98 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.646 on 198 degrees of freedom
## Multiple R-squared: 0.8847, Adjusted R-squared: 0.8841
## F-statistic: 1520 on 1 and 198 DF, p-value: < 2.2e-16
bp_res <- bptest(m_bp)
bp_res
##
## studentized Breusch-Pagan test
##
## data: m_bp
## BP = 24.741, df = 1, p-value = 6.557e-07
Ak je p-hodnota testu menšia než zvolená hladina významnosti, napríklad \(0.05\), zamietame nulovú hypotézu homoskedasticity. To znamená, že rozptyl rezíduí pravdepodobne závisí od vysvetľujúcej premennej a v modeli je prítomná heteroskedasticita.
ehat2 <- residuals(m_bp)^2
aux_bp <- lm(ehat2 ~ x, data = data_bp)
summary(aux_bp)
##
## Call:
## lm(formula = ehat2 ~ x, data = data_bp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.098 -6.205 -2.066 2.126 78.775
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.615 1.975 -1.324 0.187
## x 1.718 0.325 5.287 3.27e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.3 on 198 degrees of freedom
## Multiple R-squared: 0.1237, Adjusted R-squared: 0.1193
## F-statistic: 27.95 on 1 and 198 DF, p-value: 3.273e-07
n * summary(aux_bp)$r.squared
## [1] 24.7411
Ak pomocná regresia vysvetľuje podstatnú časť variability \(e_i^2\), znamená to, že rozptyl chýb nie je konštantný. Vysvetľujúca premenná teda pomáha predpovedať veľkosť chýb.
white_hc0 <- coeftest(m_bp, vcov = vcovHC(m_bp, type = "HC0"))
white_hc1 <- coeftest(m_bp, vcov = vcovHC(m_bp, type = "HC1"))
white_hc0
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.163533 0.373375 5.7945 2.664e-08 ***
## x 2.965347 0.080514 36.8301 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
white_hc1
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.16353 0.37526 5.7655 3.088e-08 ***
## x 2.96535 0.08092 36.6455 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
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.
V balíku nlme môžeme modelovať heteroskedasticitu
napríklad pomocou štruktúry varPower.
gls_het <- gls(
y ~ x,
data = data_bp,
weights = varPower(form = ~ x)
)
summary(gls_het)
## Generalized least squares fit by REML
## Model: y ~ x
## Data: data_bp
## AIC BIC logLik
## 918.1181 931.2711 -455.059
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~x
## Parameter estimates:
## power
## 0.750207
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 2.142065 0.26107494 8.20479 0
## x 2.969510 0.06265433 47.39513 0
##
## Correlation:
## (Intr)
## x -0.85
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.0121760 -0.6537245 -0.1084361 0.6083592 3.1140187
##
## Residual standard error: 0.7054642
## Degrees of freedom: 200 total; 198 residual
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 premennou \(x\), a preto môžu byť odhady efektívnejšie.
Teraz vytvoríme dáta s autokorelovanými chybami prvého rádu.
set.seed(123)
n <- 200
x <- rnorm(n)
u <- numeric(n)
eps <- rnorm(n, sd = 1)
rho <- 0.7
u[1] <- eps[1]
for (t in 2:n) {
u[t] <- rho * u[t - 1] + eps[t]
}
y <- 1 + 2 * x + u
data_bg <- data.frame(
y = y,
x = x,
t = 1:n
)
m_bg <- lm(y ~ x, data = data_bg)
summary(m_bg)
##
## Call:
## lm(formula = y ~ x, data = data_bg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7937 -0.9739 -0.0469 0.8434 3.3399
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.15262 0.09073 12.70 <2e-16 ***
## x 2.06205 0.09644 21.38 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.283 on 198 degrees of freedom
## Multiple R-squared: 0.6978, Adjusted R-squared: 0.6963
## F-statistic: 457.2 on 1 and 198 DF, p-value: < 2.2e-16
bg_res_1 <- bgtest(m_bg, order = 1)
bg_res_1
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: m_bg
## LM test = 79.796, df = 1, p-value < 2.2e-16
bg_res_2 <- bgtest(m_bg, order = 2)
bg_res_2
##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: m_bg
## LM test = 80.315, df = 2, p-value < 2.2e-16
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é.
ehat <- residuals(m_bg)
aux_bg_data <- data.frame(
ehat = ehat[-1],
x = data_bg$x[-1],
ehat_lag1 = ehat[-length(ehat)]
)
aux_bg <- lm(ehat ~ x + ehat_lag1, data = aux_bg_data)
summary(aux_bg)
##
## Call:
## lm(formula = ehat ~ x + ehat_lag1, data = aux_bg_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.51584 -0.69443 -0.03918 0.66062 2.39618
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.01436 0.07011 -0.205 0.838
## x -0.07870 0.07477 -1.053 0.294
## ehat_lag1 0.63547 0.05515 11.522 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.989 on 196 degrees of freedom
## Multiple R-squared: 0.4038, Adjusted R-squared: 0.3978
## F-statistic: 66.39 on 2 and 196 DF, p-value: < 2.2e-16
nrow(aux_bg_data) * summary(aux_bg)$r.squared
## [1] 80.36491
Ak je koeficient pri oneskorenom rezíduu významný, znamená to, že minulé chyby pomáhajú vysvetliť súčasné chyby. To je typický znak autokorelácie.
hac_nw <- coeftest(m_bg, vcov = NeweyWest(m_bg))
hac_nw
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.15262 0.19178 6.0101 8.768e-09 ***
## x 2.06205 0.10274 20.0713 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Aj tu zostávajú koeficienty 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.
Pri autokorelácii môžeme v gls() špecifikovať napríklad
korelačnú štruktúru AR(1).
gls_ar1 <- gls(
y ~ x,
data = data_bg,
correlation = corAR1(form = ~ t)
)
summary(gls_ar1)
## Generalized least squares fit by REML
## Model: y ~ x
## Data: data_bg
## AIC BIC logLik
## 574.9673 588.1203 -283.4836
##
## Correlation Structure: AR(1)
## Formula: ~t
## Parameter estimate(s):
## Phi
## 0.6476351
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1.159708 0.19671994 5.89523 0
## x 1.981823 0.06075463 32.62011 0
##
## Correlation:
## (Intr)
## x 0.007
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.17738142 -0.71702998 -0.03610559 0.66145653 2.58582938
##
## Residual standard error: 1.298291
## Degrees of freedom: 200 total; 198 residual
GLS teraz priamo modeluje to, že chyby v čase spolu súvisia. Preto sa môžu zmeniť nielen štandardné chyby, ale aj samotné odhadnuté koeficienty.
coef_ols_bp <- coef(m_bp)
coef_gls_bp <- coef(gls_het)
coef_ols_bg <- coef(m_bg)
coef_gls_bg <- coef(gls_ar1)
coef_ols_bp
## (Intercept) x
## 2.163533 2.965347
coef_gls_bp
## (Intercept) x
## 2.142065 2.969510
coef_ols_bg
## (Intercept) x
## 1.152621 2.062053
coef_gls_bg
## (Intercept) x
## 1.159708 1.981823
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.
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í:
Takto majú študenti k dispozícii ucelenú logickú cestu od diagnostiky až po nápravné metódy.
V prítomnosti heteroskedasticity je klasický OLS odhad rozptylu nekonzistentný. Na riešenie tohto problému používame heteroskedasticity-konzistentné (HC) odhady kovariančnej matice.
Uvažujme lineárny regresný model:
\[ y = X\beta + u \]
Whiteov (1980) robustný odhad kovariančnej matice má tvar:
\[ \widehat{\mathrm{Var}}(\hat{\beta}) = (X'X)^{-1} \left( \sum_{i=1}^n x_i x_i' e_i^2 \right) (X'X)^{-1} \]
Tento odhad sa nazýva HC0.
Aj keď je HC0 asymptoticky konzistentný, v malých výberoch má tendenciu:
Na nápravu sa používajú alternatívne odhady, ktoré upravujú štvorce rezíduí \(e_i^2\).
\[ e_i^2 \]
\[ \frac{n}{n - k} e_i^2 \]
\[ \frac{e_i^2}{1 - h_i} \]
kde \(h_i\) je leverage (diagonálny prvok tzv. hat matice).
\[ \frac{e_i^2}{(1 - h_i)^2} \]
| Odhad | Úprava | Typické použitie |
|---|---|---|
| HC0 | Žiadna | Veľké výbery |
| HC1 | Stupne voľnosti | Bežná prax |
| HC2 | Korekcia leverage | Stredné výbery |
| HC3 | Silná korekcia leverage | Malé výbery (odporúčaný) |
| HC4+ | Extrémne leverage | Špeciálne prípady |
Interpretácia :
Whiteove robustné štandardné chyby nemenia odhady koeficientov, ale upravujú ich štandardné chyby.
Ak sa po tejto úprave zmení štatistická významnosť, pôvodné OLS štandardné chyby boli nespoľahlivé.
Uvažujme lineárny regresný model:
\[ y_t = \beta_0 + \beta_1 x_{1t} + \cdots + \beta_k x_{kt} + u_t \]
Klasický predpoklad OLS je, že chyby sú nekorelované:
\[ \mathrm{Cov}(u_t, u_s) = 0 \quad \text{pre } t \neq s \]
Ak tento predpoklad neplatí, hovoríme o autokorelácii rezíduí.
Najčastejšie sa predpokladá autoregresná štruktúra chýb, napríklad AR(1):
\[ u_t = \rho u_{t-1} + \varepsilon_t, \quad |\rho| < 1 \]
kde:
Pre AR(1) proces platí:
\[ \mathrm{Cov}(u_t, u_{t-1}) = \rho \sigma^2 \]
a všeobecnejšie:
\[ \mathrm{Cov}(u_t, u_{t-s}) = \rho^s \sigma^2 \]
To znamená, že:
Pre model s autokoreláciou má kovariančná matica tvar:
\[ \Omega = \begin{pmatrix} \sigma^2 & \rho\sigma^2 & \rho^2\sigma^2 & \cdots \\ \rho\sigma^2 & \sigma^2 & \rho\sigma^2 & \cdots \\ \rho^2\sigma^2 & \rho\sigma^2 & \sigma^2 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix} \]
Táto matica:
Pri autokorelácii platí:
Existujú dva hlavné prístupy:
Použijeme HAC (Newey–West) štandardné chyby:
coeftest(m_bg, vcov = NeweyWest(m_bg))
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.15262 0.19178 6.0101 8.768e-09 ***
## x 2.06205 0.10274 20.0713 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
alebo GLS
library(nlme)
gls_ar1 <- gls(
y ~ x,
data = data_bg,
correlation = corAR1(form = ~ t)
)
summary(gls_ar1)
## Generalized least squares fit by REML
## Model: y ~ x
## Data: data_bg
## AIC BIC logLik
## 574.9673 588.1203 -283.4836
##
## Correlation Structure: AR(1)
## Formula: ~t
## Parameter estimate(s):
## Phi
## 0.6476351
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1.159708 0.19671994 5.89523 0
## x 1.981823 0.06075463 32.62011 0
##
## Correlation:
## (Intr)
## x 0.007
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.17738142 -0.71702998 -0.03610559 0.66145653 2.58582938
##
## Residual standard error: 1.298291
## Degrees of freedom: 200 total; 198 residual