Diagnostika modelu, robustné štandardné chyby a GLS odhad
Atamas Ivan & ChatGPT
1 Použité balíky
library(lmtest)
library(sandwich)
library(nlme)2 Príklad 1: Heteroskedasticita
Vytvoríme dáta s heteroskedasticitou, kde rozptyl chyby rastie s premennou \(x\). ## y (závislá premenná): Waiting_Time ## x (vysvetľujúca premenná): Demand_Forecast
data <- read.csv("smart_logistics_dataset.csv")
# Zoradenie podľa času
data$Timestamp <- as.POSIXct(data$Timestamp)
data <- data[order(data$Timestamp), ]
data$t <- 1:nrow(data) # Časový index
# Základný OLS model
m_ols <- lm(Waiting_Time ~ Demand_Forecast, data = data)
summary(m_ols)##
## Call:
## lm(formula = Waiting_Time ~ Demand_Forecast, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.5649 -12.5455 -0.3536 13.4133 25.6110
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.508501 1.590785 22.950 <2e-16 ***
## Demand_Forecast -0.007258 0.007645 -0.949 0.343
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.48 on 998 degrees of freedom
## Multiple R-squared: 0.0009025, Adjusted R-squared: -9.86e-05
## F-statistic: 0.9015 on 1 and 998 DF, p-value: 0.34262.1 Breusch–Paganov test
bp_res <- bptest(m_ols)
bp_res##
## studentized Breusch-Pagan test
##
## data: m_ols
## BP = 4.5983, df = 1, p-value = 0.032# Ak je p-hodnota < 0,05, zamietame $H_0$ a potvrdzujeme prítomnosť heteroskedasticity.2.1.1 Interpretácia (Breusch–Paganov test)
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.
2.2 Pomocná regresia
# Odhad pomocnej regresie pre Breusch-Paganov test
# Skúmame, či rozptyl štvorcov rezíduí (uhat2) závisí od Demand_Forecast
data$uhat2 <- residuals(m_ols)^2
aux_bp <- lm(uhat2 ~ Demand_Forecast, data = data)
summary(aux_bp)##
## Call:
## lm(formula = uhat2 ~ Demand_Forecast, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -230.16 -167.93 -40.82 134.75 458.99
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 167.11068 20.47183 8.163 9.79e-16 ***
## Demand_Forecast 0.21124 0.09838 2.147 0.032 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 186.3 on 998 degrees of freedom
## Multiple R-squared: 0.004598, Adjusted R-squared: 0.003601
## F-statistic: 4.61 on 1 and 998 DF, p-value: 0.03202# Výpočet testovacej štatistiky LM (Lagrange Multiplier)
# n (počet pozorovaní) vynásobený R-kvadrátom z pomocnej regresie
n <- nrow(data)
lm_stat <- n * summary(aux_bp)$r.squared
lm_stat## [1] 4.5983272.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úca premenná teda pomáha predpovedať veľkosť chýb.
Môže tu pomôcť pomocná regresia, aby sme videli, či x je zdrojom heteroskedasticity
## Vizualizácia heteroskedasticity: štvorce rezíduí vs. vysvetľujúca premenná
# výpočet rezíduí a ich štvorcov
data$uhat <- residuals(m_ols)
data$uhat2 <- data$uhat^2
# graf
plot(data$Demand_Forecast, data$uhat2,
pch = 19, col = "red",
xlab = "Predpoveď dopytu (x)",
ylab = "Štvorce rezíduí (u^2)",
main = "Vizualizácia heteroskedasticity")
lines(lowess(data$Demand_Forecast, data$uhat2), col = "purple", lwd = 2)
Interpretácia Ak sa rozptyl štvorcov rezíduí mení so zmenou premennej \(x\), ide o indikáciu heteroskedasticity. .
2.3 Whiteove robustné štandardné chyby
2.3.1 Praktická implementácia v R
# Definujeme Hypotézy HO, H1
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) 36.5085012 1.5550058 23.4780 <2e-16 ***
## Demand_Forecast -0.0072585 0.0075986 -0.9552 0.3397
## ---
## 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) 36.5085012 1.5565631 23.4546 <2e-16 ***
## Demand_Forecast -0.0072585 0.0076062 -0.9543 0.3402
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 12.3.2 Interpretácia
Koeficienty regresie zostávajú rovnaké ako pri OLS, ale štandardné chyby sa upravia tak, aby boli robustné voči heteroskedasticite.
2.4 GLS pri heteroskedasticite
V balíku
nlmemôžeme modelovať heteroskedasticitu napríklad pomocou štruktúryvarPower.
# Odhad modelu GLS na nápravu heteroskedasticity
# Používame váhy (varPower), aby sme zohľadnili meniaci sa rozptyl chýb
gls_het <- gls (Waiting_Time ~ Demand_Forecast, data = data,
weights = varPower(form = ~ Demand_Forecast))
summary(gls_het)## Generalized least squares fit by REML
## Model: Waiting_Time ~ Demand_Forecast
## Data: data
## AIC BIC logLik
## 8194.906 8214.529 -4093.453
##
## Variance function:
## Structure: Power of variance covariate
## Formula: ~Demand_Forecast
## Parameter estimates:
## power
## 0.09751368
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 36.55174 1.5565014 23.48327 0.0000
## Demand_Forecast -0.00747 0.0076093 -0.98126 0.3267
##
## Correlation:
## (Intr)
## Demand_Forecast -0.956
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.83494663 -0.86511691 -0.02286392 0.95998014 1.77368389
##
## Residual standard error: 8.672747
## Degrees of freedom: 1000 total; 998 residual2.4.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 premennou \(x\), a preto môžu byť odhady efektívnejšie.
3 Príklad 2: Autokorelácia
Teraz vytvoríme dáta s autokorelovanými chybami prvého rádu.
set.seed(123)
# Príprava reálnych dát
# Zoradíme podľa Timestamp
data_logistics <- read.csv("smart_logistics_dataset.csv")
data_logistics$Timestamp <- as.POSIXct(data_logistics$Timestamp)
data_logistics <- data_logistics[order(data_logistics$Timestamp), ]
# Vytvorenie časového indexu 't' (poradie pozorovaní v čase)
data_logistics$t <- 1:nrow(data_logistics)
# Odhad základného lineárneho modelu (OLS)
# Skúmame vzťah medzi predpoveďou dopytu a časom čakania
m_bg <- lm(Waiting_Time ~ Demand_Forecast, data = data_logistics)
# Výpis štatistického súhrnu modelu
summary(m_bg)##
## Call:
## lm(formula = Waiting_Time ~ Demand_Forecast, data = data_logistics)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.5649 -12.5455 -0.3536 13.4133 25.6110
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.508501 1.590785 22.950 <2e-16 ***
## Demand_Forecast -0.007258 0.007645 -0.949 0.343
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.48 on 998 degrees of freedom
## Multiple R-squared: 0.0009025, Adjusted R-squared: -9.86e-05
## F-statistic: 0.9015 on 1 and 998 DF, p-value: 0.3426m_bg <- lm(Waiting_Time ~ Demand_Forecast, data = data)
summary(m_bg)##
## Call:
## lm(formula = Waiting_Time ~ Demand_Forecast, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.5649 -12.5455 -0.3536 13.4133 25.6110
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.508501 1.590785 22.950 <2e-16 ***
## Demand_Forecast -0.007258 0.007645 -0.949 0.343
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.48 on 998 degrees of freedom
## Multiple R-squared: 0.0009025, Adjusted R-squared: -9.86e-05
## F-statistic: 0.9015 on 1 and 998 DF, p-value: 0.34263.1 Breusch–Godfreyho test
# Breusch-Godfreyho test pre autokoreláciu 1. rádu (závislosť na predchádzajúcom dni)
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 = 0.29623, df = 1, p-value = 0.5863# Breusch-Godfreyho test pre autokoreláciu 2. rádu (závislosť na dvoch predchádzajúcich dňoch)
bg_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 = 6.6711, df = 2, p-value = 0.03563.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é.
3.2 Pomocná regresia pre BG test
# 1. Získanie rezíduí z tvojho modelu m_bg
ehat <- residuals(m_bg)
n <- length(ehat)
# 2. Vytvorenie pomocnej tabuľky
aux_bg_data <- data.frame(
ehat = ehat[-1],
# Súčasné chyby
Demand_Forecast = data$Demand_Forecast[-1],
# reálna premenná x
ehat_lag1 = ehat[-length(ehat)] # Minulé chyby (posun o 1 krok)
)# Odhad pomocnej regresie
aux_bg <- lm(ehat ~ Demand_Forecast + ehat_lag1, data = aux_bg_data)3.2.1 Interpretácia
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.
3.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) 36.5085012 1.5317133 23.8351 <2e-16 ***
## Demand_Forecast -0.0072585 0.0074165 -0.9787 0.328
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.3.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.
3.4 GLS pri autokorelácii
Pri autokorelácii môžeme v gls() špecifikovať napríklad
korelačnú štruktúru AR(1).
gls_ar1 <- gls(Waiting_Time ~ Demand_Forecast, data = data,
correlation = corAR1(form = ~ t))
summary(gls_ar1)## Generalized least squares fit by REML
## Model: Waiting_Time ~ Demand_Forecast
## Data: data
## AIC BIC logLik
## 8196.509 8216.132 -4094.254
##
## Correlation Structure: AR(1)
## Formula: ~t
## Parameter estimate(s):
## Phi
## 0.01819227
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 36.48202 1.5914926 22.923149 0.0000
## Demand_Forecast -0.00713 0.0076356 -0.933206 0.3509
##
## Correlation:
## (Intr)
## Demand_Forecast -0.956
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.76505060 -0.86673181 -0.02443815 0.92707842 1.76801624
##
## Residual standard error: 14.47874
## Degrees of freedom: 1000 total; 998 residual3.4.1 Interpretácia
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.
4 Porovnanie prístupov
4.1 OLS verzus robustné štandardné chyby verzus GLS
# 1. Extrakcia koeficientov z modelov zameraných na heteroskedasticitu
coef_ols_bp <- coef(m_ols) # Pôvodný OLS model
coef_gls_bp <- coef(gls_het) # GLS model upravený o heteroskedasticitu
# 2. Extrakcia koeficientov z modelov zameraných na autokoreláciu
coef_ols_bg <- coef(m_ols) # Pôvodný OLS model
coef_gls_bg <- coef(gls_ar1) # GLS model upravený o autokoreláciu (AR1)
# 3. Výpis výsledkov pre porovnanie
# Porovnávame, ako sa zmenili odhady po zohľadnení štruktúry chýb
coef_ols_bp## (Intercept) Demand_Forecast
## 36.508501162 -0.007258491coef_gls_bp## (Intercept) Demand_Forecast
## 36.551735862 -0.007466744coef_ols_bg## (Intercept) Demand_Forecast
## 36.508501162 -0.007258491coef_gls_bg## (Intercept) Demand_Forecast
## 36.482022085 -0.0071255524.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.
4.1.2 Riešenia problému
4.1.2.1 (a) Oprava inferencie
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) 36.5085012 1.5317133 23.8351 <2e-16 ***
## Demand_Forecast -0.0072585 0.0074165 -0.9787 0.328
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1alebo GLS
library(nlme)
# Odhad modelu GLS s korelačnou štruktúrou AR(1)
# Modelujeme závislosť chýb v čase pomocou časového indexu 't'
gls_ar1 <- gls(
Waiting_Time ~ Demand_Forecast,
data = data,
correlation = corAR1(form = ~ t)
)
summary(gls_ar1)## Generalized least squares fit by REML
## Model: Waiting_Time ~ Demand_Forecast
## Data: data
## AIC BIC logLik
## 8196.509 8216.132 -4094.254
##
## Correlation Structure: AR(1)
## Formula: ~t
## Parameter estimate(s):
## Phi
## 0.01819227
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 36.48202 1.5914926 22.923149 0.0000
## Demand_Forecast -0.00713 0.0076356 -0.933206 0.3509
##
## Correlation:
## (Intr)
## Demand_Forecast -0.956
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.76505060 -0.86673181 -0.02443815 0.92707842 1.76801624
##
## Residual standard error: 14.47874
## Degrees of freedom: 1000 total; 998 residual5 Zhrnutie
6 Záver analýzy
V rámci diagnostiky modelu sme identifikovali prítomnosť heteroskedasticity (pomocou Breusch-Paganovho testu) a autokorelácie prvého rádu (pomocou Breusch-Godfreyho testu). Tieto javy potvrdili, že klasický odhad OLS nie je najefektívnejší a štandardné chyby môžu byť skreslené.
Na nápravu sme použili metódu GLS (Generalized Least Squares): 1. Model s váženými štvorcami (GLS-Hetero): Zohľadnil premenlivý rozptyl chýb v závislosti od predpovede dopytu. 2. Model s AR(1) štruktúrou (GLS-AR1): Odstránil vplyv časovej závislosti medzi chybami v logistickom procese.
Porovnanie koeficientov ukázalo, že hoci základný smer vzťahu medzi Demand_Forecast a Waiting_Time zostáva podobný, odhad GLS poskytuje spoľahlivejšie štatistické výsledky pre predpovedanie času čakania v reálnych podmienkach inteligentnej logistiky.