Úvod

Tento notebook vysvetľuje tri nadväzujúce oblasti diagnostiky a nápravy lineárneho regresného modelu:

  1. Breusch–Paganov test pre heteroskedasticitu,
  2. Breusch–Godfreyho test pre autokoreláciu,
  3. robustné štandardné chyby a GLS odhad ako nápravné metódy.

Použité balíky

library(lmtest)
library(sandwich)
library(nlme)

Príklad 1: Heteroskedasticita

Použijeme naše dáta - company_data.csv

library(dplyr)
udaje_nove <- read.csv2("udaje/company_data.csv",header=TRUE,sep=",",dec=".")
model <- lm(Profit ~ Price + Order.Qty + Unit.Cost, data = udaje_nove)
summary(model)
## 
## Call:
## lm(formula = Profit ~ Price + Order.Qty + Unit.Cost, data = udaje_nove)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -17243   -365   -100    -23  42237 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -28.49945   20.64798   -1.38    0.168    
## Price        11.62553    0.06444  180.40   <2e-16 ***
## Order.Qty     9.48104    0.33894   27.97   <2e-16 ***
## Unit.Cost   -11.82672    0.13620  -86.83   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1799 on 14996 degrees of freedom
## Multiple R-squared:  0.7179, Adjusted R-squared:  0.7178 
## F-statistic: 1.272e+04 on 3 and 14996 DF,  p-value: < 2.2e-16

Breusch–Paganov test

bp_res <- bptest(model)
bp_res
## 
##  studentized Breusch-Pagan test
## 
## data:  model
## BP = 1106.1, df = 3, p-value < 2.2e-16

Interpretácia

Keďže p-hodnota testu (\(< 2.2 \times 10^{-16}\)) je výrazne menšia než zvolená hladina významnosti \(0.05\), zamietame nulovú hypotézu o homoskedasticite. To znamená, že rozptyl rezíduí pravdepodobne závisí od vysvetľujúcej premennej a v modeli je prítomná heteroskedasticita. ## Pomocná regresia

ehat2 <- residuals(model)^2
aux_bp <- lm(ehat2 ~ Price + Order.Qty + Unit.Cost, data = udaje_nove)
summary(aux_bp)
## 
## Call:
## lm(formula = ehat2 ~ Price + Order.Qty + Unit.Cost, data = udaje_nove)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
##  -86786841   -3300657      86733    1859752 1747952522 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3265306     340758  -9.582   <2e-16 ***
## Price          32552       1064  30.607   <2e-16 ***
## Order.Qty      78633       5594  14.058   <2e-16 ***
## Unit.Cost     -35902       2248 -15.972   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 29690000 on 14996 degrees of freedom
## Multiple R-squared:  0.07374,    Adjusted R-squared:  0.07356 
## F-statistic:   398 on 3 and 14996 DF,  p-value: < 2.2e-16
n <- nrow(udaje_nove)
n * summary(aux_bp)$r.squared
## [1] 1106.142

Interpretácia

Pomocná regresia vysvetľuje podstatnú časť variability \(e_i^2\), znamená to, že rozptyl chýb nie je konštantný. Vysvetľujúce premenné teda pomáhajú predpovedať veľkosť chýb.

Pomocná regresia nám v tomto prípade pomohla odhaliť, že práve tieto premenné sú zdrojom heteroskedasticity v našom modeli.

## Vizualizácia heteroskedasticity: štvorce rezíduí vs. vysvetľujúca premenná pre všetky 3 yvsvetlujúce premenné

# výpočet rezíduí a ich štvorcov
udaje_nove$uhat <- residuals(model)
udaje_nove$uhat2 <- udaje_nove$uhat^2

# graf 1 - vysvetlujúca premenná Price
plot(
  udaje_nove$Price, udaje_nove$uhat2,
  pch = 19,
  col = "darkblue",
  xlab = "Vysvetľujúca premenná Price",
  ylab = expression(hat(u)[i]^2),
  main = "Štvorce rezíduí vs. Premenná Price (indikácia heteroskedasticity)"
)

# LOESS vyhladenie (trend)
lines(
  lowess(udaje_nove$Price, udaje_nove$uhat2),
  col = "red",
  lwd = 2
)

# graf 2 - vysvetlujúca premenná Order.Qty
plot(
  udaje_nove$Order.Qty, udaje_nove$uhat2,
  pch = 19,
  col = "darkblue",
  xlab = "Vysvetľujúca premenná Order.Qty",
  ylab = expression(hat(u)[i]^2),
  main = "Štvorce rezíduí vs. Premenná Order.Qty (indikácia heteroskedasticity)"
)

# LOESS vyhladenie (trend)
lines(
  lowess(udaje_nove$Order.Qty, udaje_nove$uhat2),
  col = "red",
  lwd = 2
)

# graf 3 - vysvetlujúca premenná Unit.Cost
plot(
  udaje_nove$Unit.Cost, udaje_nove$uhat2,
  pch = 19,
  col = "darkblue",
  xlab = "Vysvetľujúca premenná Unit.Cost",
  ylab = expression(hat(u)[i]^2),
  main = "Štvorce rezíduí vs. Premenná Unit.Cost (indikácia heteroskedasticity)"
)

# LOESS vyhladenie (trend)
lines(
  lowess(udaje_nove$Unit.Cost, udaje_nove$uhat2),
  col = "red",
  lwd = 2
)

Interpretácia >Rozptyl štvorcov rezíduí (\(e_i^2\)) sa mení so zmenou premennej Price, Order.Qty alebo Unit.Cost, ide o jasnú indikáciu heteroskedasticity. Červená LOESS krivka v každom z grafov pomáha identifikovať tento systematický vzťah medzi danou premennou a variabilitou chýb (rezíduí) modelu.“

Praktická implementácia v R

white_hc0 <- coeftest(model, vcov = vcovHC(model, type = "HC0"))
white_hc1 <- coeftest(model, vcov = vcovHC(model, type = "HC1"))

white_hc0
## 
## t test of coefficients:
## 
##              Estimate Std. Error  t value Pr(>|t|)    
## (Intercept) -28.49945   32.65044  -0.8729   0.3828    
## Price        11.62553    0.21498  54.0785   <2e-16 ***
## Order.Qty     9.48104    1.14709   8.2653   <2e-16 ***
## Unit.Cost   -11.82672    0.33244 -35.5752   <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) -28.49945   32.65480  -0.8727   0.3828    
## Price        11.62553    0.21500  54.0713   <2e-16 ***
## Order.Qty     9.48104    1.14724   8.2642   <2e-16 ***
## Unit.Cost   -11.82672    0.33249 -35.5705   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Koeficienty regresie zostávajú rovnaké ako pri pôvodnom OLS odhade, ale štandardné chyby boli upravené tak, aby boli robustné voči heteroskedasticite. Keďže všetky koeficienty premenných Price, Order.Qty a Unit.Cost zostávajú štatisticky významné aj po tejto úprave (p-hodnoty sú naďalej výrazne pod \(0.001\)), môžeme konštatovať, že ich významnosť je robustná voči prítomnosti heteroskedasticity. ## GLS pri heteroskedasticite

V balíku nlme môžeme modelovať heteroskedasticitu napríklad pomocou štruktúry varPower.

gls_het <- gls(
  Profit ~ Price + Order.Qty + Unit.Cost,
  data = udaje_nove,
  weights = varPower(form = ~ Price)
)

summary(gls_het)
## Generalized least squares fit by REML
##   Model: Profit ~ Price + Order.Qty + Unit.Cost 
##   Data: udaje_nove 
##        AIC      BIC    logLik
##   255340.9 255386.6 -127664.4
## 
## Variance function:
##  Structure: Power of variance covariate
##  Formula: ~Price 
##  Parameter estimates:
##     power 
## 0.4309234 
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) -43.30093  7.862672  -5.50715       0
## Price        12.18040  0.082751 147.19296       0
## Order.Qty     3.44687  0.058339  59.08300       0
## Unit.Cost   -12.29750  0.166349 -73.92611       0
## 
##  Correlation: 
##           (Intr) Price  Ordr.Q
## Price     -0.208              
## Order.Qty -0.366  0.056       
## Unit.Cost -0.066 -0.816  0.016
## 
## Standardized residuals:
##          Min           Q1          Med           Q3          Max 
## -21.05581522  -0.26577210  -0.05966600   0.01376072  25.03541618 
## 
## Residual standard error: 138.1981 
## Degrees of freedom: 15000 total; 14996 residual

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 premennou \(x\), a preto môžu byť odhady efektívnejšie.

Príklad 2: Autokorelácia

Breusch–Godfreyho test

bg_res_1 <- bgtest(model, order = 1)
bg_res_1
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model
## LM test = 1.505, df = 1, p-value = 0.2199
bg_res_2 <- bgtest(model, order = 2)
bg_res_2
## 
##  Breusch-Godfrey test for serial correlation of order up to 2
## 
## data:  model
## LM test = 4.5125, df = 2, p-value = 0.1047

Test prvého rádu nepotvrdil štatisticky významnú autokoreláciu. Avšak test druhého rádu preukázal štatisticky významnú autokoreláciu (p-hodnota < 0,05), teda 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é. ## Pomocná regresia pre BG test

ehat <- residuals(model)

aux_bg_data <- data.frame(
  ehat      = ehat[-1],               
  ehat_lag1 = ehat[-length(ehat)],               
  Price     = udaje_nove$Price[-1],   
  Order.Qty = udaje_nove$Order.Qty[-1],
  Unit.Cost = udaje_nove$Unit.Cost[-1]
)

aux_bg <- lm(ehat ~ Price + Order.Qty + Unit.Cost + ehat_lag1, data = aux_bg_data)
summary(aux_bg)
## 
## Call:
## lm(formula = ehat ~ Price + Order.Qty + Unit.Cost + ehat_lag1, 
##     data = aux_bg_data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -17242   -366   -101    -21  42236 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1687593 20.6494594  -0.008    0.993
## Price       -0.0003148  0.0644459  -0.005    0.996
## Order.Qty    0.0006486  0.3389411   0.002    0.998
## Unit.Cost    0.0023740  0.1362191   0.017    0.986
## ehat_lag1   -0.0100182  0.0081675  -1.227    0.220
## 
## Residual standard error: 1799 on 14994 degrees of freedom
## Multiple R-squared:  0.0001003,  Adjusted R-squared:  -0.0001664 
## F-statistic: 0.3761 on 4 and 14994 DF,  p-value: 0.8258
nrow(aux_bg_data) * summary(aux_bg)$r.squared
## [1] 1.504869

Interpretácia

Koeficient pri oneskorenom rezíduu je významný,to znamená, že minulé chyby pomáhajú vysvetliť súčasné chyby. To je typický znak autokorelácie.

HAC štandardné chyby

hac_nw <- coeftest(model, vcov = NeweyWest(model))
hac_nw
## 
## t test of coefficients:
## 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) -28.49945   33.14517  -0.8598    0.3899    
## Price        11.62553    0.21482  54.1179 < 2.2e-16 ***
## Order.Qty     9.48104    1.16989   8.1042 5.717e-16 ***
## Unit.Cost   -11.82672    0.33539 -35.2627 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretácia

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.

GLS pri autokorelácii

Pri autokorelácii môžeme v gls() špecifikovať napríklad korelačnú štruktúru AR(1).

udaje_nove <- read.csv2("udaje/company_data.csv",header=TRUE,sep=",",dec=".")
model <- lm(Profit ~ Price + Order.Qty + Unit.Cost, data = udaje_nove)
set.seed(123)
vzorka_index <- sample(1:nrow(udaje_nove), 2000)
data_vzorka <- udaje_nove[vzorka_index, ]
data_vzorka$t <- 1:nrow(data_vzorka)
gls_vzorka <- gls(
  Profit ~ Price + Order.Qty + Unit.Cost, 
  data = data_vzorka, 
  correlation = corAR1(form = ~ t)
)
summary(gls_vzorka)
## Generalized least squares fit by REML
##   Model: Profit ~ Price + Order.Qty + Unit.Cost 
##   Data: data_vzorka 
##       AIC      BIC    logLik
##   35422.3 35455.89 -17705.15
## 
## Correlation Structure: AR(1)
##  Formula: ~t 
##  Parameter estimate(s):
##          Phi 
## -0.008185088 
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) -80.68032  54.21904  -1.48804  0.1369
## Price        11.22447   0.16761  66.96902  0.0000
## Order.Qty    10.08750   0.84613  11.92186  0.0000
## Unit.Cost   -10.58010   0.36861 -28.70295  0.0000
## 
##  Correlation: 
##           (Intr) Price  Ordr.Q
## Price     -0.325              
## Order.Qty -0.327  0.049       
## Unit.Cost -0.177 -0.713  0.017
## 
## Standardized residuals:
##          Min           Q1          Med           Q3          Max 
## -5.001717111 -0.218940198 -0.058086599  0.005726058 14.134147404 
## 
## Residual standard error: 1697.083 
## Degrees of freedom: 2000 total; 1996 residual

Bohužiaľ, môj počítač nezvládol tento test spraviť na všetkých údajoch, tak som musela vybrať 2000 náhodných. Výsledky potvrdili, že odhady koeficientov sú štatisticky významné a ekonomicky interpretovateľné: rast cien a množstva objednávok pozitívne ovplyvňuje zisk, zatiaľ čo rast jednotkových nákladov ho znižuje. Tento odhad potvrdil robustnosť modelu aj po zohľadnení časovej korelácie v rezíduách. Na internete som našla aj podobný test, ktorý sa mi podarilo spustiť aj na všetkých mojich údajoch:

##HAC test

udaje_nove <- read.csv2("udaje/company_data.csv",header=TRUE,sep=",",dec=".")
model <- lm(Profit ~ Price + Order.Qty + Unit.Cost, data = udaje_nove)
hac_results <- coeftest(model, vcov = NeweyWest(model))
print(hac_results)
## 
## t test of coefficients:
## 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) -28.49945   33.14517  -0.8598    0.3899    
## Price        11.62553    0.21482  54.1179 < 2.2e-16 ***
## Order.Qty     9.48104    1.16989   8.1042 5.717e-16 ***
## Unit.Cost   -11.82672    0.33539 -35.2627 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretácia

Výsledky HAC testu potvrdili, že všetky premenné v modeli majú štatisticky významný vplyv na zisk, pričom aj po oprave štandardných chýb ostali tieto závery stabilné.To znamená, že autokorelácia v dátach nespôsobila výrazné skreslenie a pôvodný model je dostatočne spoľahlivý.

Porovnanie prístupov

OLS verzus robustné štandardné chyby verzus GLS

coef_ols_bp <- coef(model)
coef_gls_bp <- coef(gls_het)

coef_ols_bg <- coef(model)
coef_gls_bg <- coef(gls_vzorka)

coef_ols_bp
## (Intercept)       Price   Order.Qty   Unit.Cost 
##  -28.499447   11.625531    9.481036  -11.826720
coef_gls_bp
## (Intercept)       Price   Order.Qty   Unit.Cost 
##  -43.300927   12.180399    3.446866  -12.297500
coef_ols_bg
## (Intercept)       Price   Order.Qty   Unit.Cost 
##  -28.499447   11.625531    9.481036  -11.826720
coef_gls_bg
## (Intercept)       Price   Order.Qty   Unit.Cost 
##   -80.68032    11.22447    10.08750   -10.58010

Interpretácia

Výsledky porovnania koeficientov ukazujú, že odhady modelu sa medzi jednotlivými prístupmi prakticky nelíšia. V tabuľke vidíme, že hodnoty koeficientov pre Price (cca 11.63), Order.Qty (cca 9.48) a Unit.Cost (cca -11.83) zostávajú pri oboch spôsoboch odhadu totožné. Hoci sme v dátach diagnostikovali heteroskedasticitu a autokoreláciu, tieto javy nespôsobili významné skreslenie vplyvu našich vysvetľujúcich premenných na zisk.