1. Úvod

Po autokoreácii a heteroskedasticite rezíduí je multikolinearita tretím závažným porušením predpokladov použitia metódy najmenších štvorcov. Tu sa okrem iného predpokadá, že matica \(\mathbf X\) je tvorená lineárne nezávislými riadkami a tiež stĺpcami, čo zabezpečí regularitu matice \(\mathbf X^T\mathbf X\) a teda možnosť jej inverzie. Tá sa používa pri odhadoch regresných koeficientov. V praxi sa ale môže stať, že vzniká "takmer singulárna matica \(\mathbf X^T\mathbf X\) ", t.j. matica \(\mathbf X\) je tvorená "približne" lineárne závislými stĺpcami, t.j. existuje taká ich lineárna kombinácia, v ktorej

\[ x_{il} = \alpha_0 + \alpha_1 x_{i1} + \dots + \alpha_{l-1} x_{i,(l-1)} + \alpha_{l+1} x_{i,(l+1)} + \alpha_k x_{i,k} + \nu_i \]

Tu \(\nu_i\) sú rádovo menšie čísla , než regresory \(x_{.}\), t.j. \((\forall i)(\nu_i << x_{.,i})\). V tomto prípade je inverzná matica \((\mathbf X^T\mathbf X)^{-1}\) veľmi nestabilná a obsahuje na hlavnej diagonále veľmi veľké hodnoty. Táto matica sa používa pri výpočtoch \[\hat \beta = (\mathbf X^T\mathbf X)^{-1} \mathbf X^T \mathbf y\] a tiež \[\text{std}(\beta_i) = \sqrt{\sigma^2 (\mathbf X^T \mathbf X)^{-1}_{ii}}.\] To spôsobuje nestabilitu odhadovaných regresných koeficientov a ich nadhodnotené rozptyly.

Tento problém nazývame problémom multikolinearity.

2. Dôsledky multikolinearity

Multikolinearita patrí medzi najčastejšie problémy viacnásobnej lineárnej regresie.
Je dôležité jasne rozlišovať dva fakty:

  1. Nespôsobuje skreslené (biased) odhady koeficientov
  2. Nadhodnocuje odhady štandardných odchýlok regresných koeficientov a vedie potom k falošnému neprijímaniu alternatívnej hypotézy o štatistickej významnosti jednotlivých regresorov.
  3. Odhadované regresné koeficienty sú nestabilné - pri malej zmene údajov sa sa prudko menia koeficienty ako aj ich znamienka.
  4. Interpretácia regresného modelu je z dôvodu vyššie uvedených dôvodov nespoľahlivá.

3. Východiskový model a údaje

Budeme pracovať s regresným modelom z minulých cvičení, pričom

\[ lead\_time_i = \beta_0 + \beta_1 arrival\_date\_week\_number_i + \beta_2 arrival\_date\_day\_of\_month_i + \beta_3 previous\_cancellations_i + \beta_4 previous\_booking\__not\_canceled_i + u_i \]

Prierezové údaje máme z databázy hotel_booking a po ich nadčítaní sme ich uložili do data.frame udaje. Z dát sme vyselektovali údaje relevantné pre našu regresnú analýzu a uložili ich do premennej udaje.new.

udaje <- read.csv("Ekonometria/hotel_bookings.csv", header = TRUE, sep = ",", dec = ".")
udaje1 <- udaje[udaje$previous_cancellations > 0,]
udaje2 <- udaje1[udaje$arrival_date_year == 2016,]
udaje.new <- udaje2[udaje2$previous_bookings_not_canceled > 0, c("lead_time", "arrival_date_week_number", "arrival_date_day_of_month", "previous_cancellations", "previous_bookings_not_canceled")]

column_medians <- sapply(udaje.new, median, na.rm = TRUE)

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

udaje.new <- udaje_imputed
udaje <- udaje.new

4. Odhad základného regresného modelu

library(ggplot2)

# your regression model
model <- lm(
  lead_time ~ arrival_date_week_number + arrival_date_day_of_month +
    previous_cancellations + previous_bookings_not_canceled,
  data = udaje.new,
  na.action = na.exclude
)
summary(model)
## 
## Call:
## lm(formula = lead_time ~ arrival_date_week_number + arrival_date_day_of_month + 
##     previous_cancellations + previous_bookings_not_canceled, 
##     data = udaje.new, na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -86.222  -0.032  -0.032  -0.032 294.449 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    -10.87357    0.42600 -25.525   <2e-16 ***
## arrival_date_week_number         0.66584    0.01366  48.747   <2e-16 ***
## arrival_date_day_of_month       -0.03851    0.02349  -1.639    0.101    
## previous_cancellations           8.89961    0.05898 150.879   <2e-16 ***
## previous_bookings_not_canceled  -0.78493    0.01350 -58.150   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.03 on 54391 degrees of freedom
## Multiple R-squared:  0.3003, Adjusted R-squared:  0.3002 
## F-statistic:  5836 on 4 and 54391 DF,  p-value: < 2.2e-16

Vo výsledkoch pozorujeme veľké rádové rozdiely v odhade regresných koeficientov.

5. Korelačná matica

Korelácia dokáže zachytiť párové vzťahy medzi premennými. Ak medzi niektorými vysvetľujúcimi premennými je vysoká korelácia (signalizujúca multikolinearitu), potom je najjednoduchšie ju zo zoznamu regresorov vylúčiť. Korelácie sa dajú aj testovať, alebo len vyčísliť a potom podľa intuitívneho pravidla vylúčiť jednu premennú, ktorá má koreláciu s inou premennou v absolútnej hodnote vyššiu ako 0.8, resp. 0.9.

xvars <- udaje[, c("arrival_date_week_number", "arrival_date_day_of_month", "previous_cancellations", "previous_bookings_not_canceled")]
round(cor(xvars), 3)
##                                arrival_date_week_number
## arrival_date_week_number                          1.000
## arrival_date_day_of_month                         0.201
## previous_cancellations                           -0.141
## previous_bookings_not_canceled                    0.160
##                                arrival_date_day_of_month previous_cancellations
## arrival_date_week_number                           0.201                 -0.141
## arrival_date_day_of_month                          1.000                 -0.108
## previous_cancellations                            -0.108                  1.000
## previous_bookings_not_canceled                     0.053                  0.367
##                                previous_bookings_not_canceled
## arrival_date_week_number                                0.160
## arrival_date_day_of_month                               0.053
## previous_cancellations                                  0.367
## previous_bookings_not_canceled                          1.000

Z výsledkov vieme povedať, že naše premenné nemajú výrazné korelácie. K najvyššej korelácií dochádza medzi previous_cancellations a previous_bookings_not_canceled s hodnotou 0.367, a teda dochádza k veľmi miernej kladnej korelácii.

Korelačný vzťah sa dá vytušiť aj z jednoduchých párových scatterplotov ako je to na nasledujúcom obrázku.

pairs(xvars,
      main = "Scatterplotová matica")

Opticky vidíme miernu koreláciu medzi previous_cancellations a previous_bookings_not_cancelled. Taktiež vidíme veľmi miernu zhodu v niektorých dátach arrival_date_week_number a previous_bookings_not_cancelled. Táto zhoda je však náhodná, a ani logicky nedáva korelácia medzi týmito premennými veľký zmysel.

6. VIF

Indikátorom multikolinearity u premennej, ktorá multikolinearitu zapríčiňuje, je Variance Inflation Factor (VIF). Pre premennú \(x_j\) je potom

\[ VIF_j = \frac{1}{1 - R_j^2} \]

kde \(R_j^2\) pochádza z regresie:

\[ X_j = \gamma_0 + \gamma_1 X_1 + \cdots + \gamma_{j-1} X_{j-1} + \gamma_{j+1} X_{j+1} + u. \]

library(car)
vif(model)
##       arrival_date_week_number      arrival_date_day_of_month 
##                       1.109678                       1.052978 
##         previous_cancellations previous_bookings_not_canceled 
##                       1.225388                       1.224777

Intuitívnym kritériom, ktoré signalizuje prítomnosť multikolinearity, je podmienka VIF > 5 (prísne kritérium), alebo VIF > 10 (menej prísne kritérium). V našom prípade to nespĺňa žiadna z vysvetľujúcich veličín.

7. Condition Number

Pri existencii multikolinearity sa model prejavuje tak, že koeficient determinácie je síce vysoký a zdá sa, že model je veľmi dobrý, ale regresné koeficienty nie sú štatisticky významné - t.j. štandardné odchýlky regresných koeficientov sú veľmi veľké. Uvedomíme si to, ak sa pozrieme, ako sa počítajú - t.j. \(\text{std}(\beta_i) = \sqrt{\sigma^2 (\mathbf X^T \mathbf X)^{-1}_{ii}}\), kde rozhodujúci je \(i\)ty prvok hlavnej diagonály matice \((\mathbf X^T \mathbf X)^{-1}\). Tie prvky sú ale v prípade podobnosti vysvetľujúcich premenných mimoriadne veľké. Túto situáciu zachytáva nasledovný ukazovateľ.

Pri výpočte Condition number \(\kappa\) sa používa vzorec

\[\kappa = \frac{\theta_{\text{max}}}{\theta_{\text{min}}}\]

kde \(\theta_.\) sú vlastné čísla matice (vysvetlené nižšie). Conditional number nie je test, je to len indikátor, ktorý posudzuje mieru multikolinearity medzi premennými. Používame intuitívne pravidlo. Ak Conditional number je

V našom prípade to vypočítame nasledovne

X <- model.matrix(model)[, -1]
XtX <- t(X) %*% X
eig <- eigen(XtX)

condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number
## [1] 93.79626

Naše condition_number spadá do intervalu 30-100, a teda signalizuje silnú multikolinearitu.

Vlastné číslo štvorcovej matice \(\mathbf X^T \mathbf X\) je číslo \(\theta_j\), pre ktoré platí \((\mathbf X^T \mathbf X)\mathbf h^j = \theta_i\mathbf h^j\). \(\mathbf h^j\) je tzv vlastný vektor tejto matice. Máme toľko vlastných čísel (teda aj vlastných vektorov), koľko obsahuje matica \(\mathbf X^T \mathbf X\) riadkov (stĺpcov).

Môže sa stať, že VIF faktor nesignalizuje multikolinearitu u žiadnej z vysvetľujúcich veličín, ale sú navzájom prepojené cyklickými lineárnymi závislosťami všetky premenné. To zachytáva práve Condition Number.

Vyzerá to, že práve toto sa deje aj v našom príklade. Poďme sa na to teda pozrieť.

8. Riešenia multikolinearity

Vynechanie premennej

Pokúsme sa vynechať postupne dve premenné, ktoré majú najvyšší VIF a porovnajme následne upravené koeficienty determinácie oboch nových modelov.

Premenné s najvyšším VIF sú v našom modeli previous_cancellations a previous_bookings_not_cancelled so skoro identickou hodnotou VIF.

model_no.previous_cancellations <- lm(
  lead_time ~ arrival_date_week_number + arrival_date_day_of_month + previous_bookings_not_canceled,
  data = udaje.new,
  na.action = na.exclude
)
summary(model_no.previous_cancellations)
## 
## Call:
## lm(formula = lead_time ~ arrival_date_week_number + arrival_date_day_of_month + 
##     previous_bookings_not_canceled, data = udaje.new, na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.758  -0.138  -0.138  -0.138 307.117 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                     4.90492    0.49185   9.972   <2e-16 ***
## arrival_date_week_number        0.26165    0.01595  16.403   <2e-16 ***
## arrival_date_day_of_month      -0.38874    0.02785 -13.961   <2e-16 ***
## previous_bookings_not_canceled  0.03423    0.01472   2.325   0.0201 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.8 on 54392 degrees of freedom
## Multiple R-squared:  0.007442,   Adjusted R-squared:  0.007387 
## F-statistic: 135.9 on 3 and 54392 DF,  p-value: < 2.2e-16
model_no.previous_bookings_not_canceled <- lm(
  lead_time ~ arrival_date_week_number + arrival_date_day_of_month +
    previous_cancellations,
  data = udaje.new,
  na.action = na.exclude
)
summary(model_no.previous_bookings_not_canceled)
## 
## Call:
## lm(formula = lead_time ~ arrival_date_week_number + arrival_date_day_of_month + 
##     previous_cancellations, data = udaje.new, na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -72.164  -0.023  -0.023  -0.023 290.794 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -11.58417    0.43886 -26.396  < 2e-16 ***
## arrival_date_week_number    0.49453    0.01375  35.977  < 2e-16 ***
## arrival_date_day_of_month  -0.11948    0.02417  -4.943  7.7e-07 ***
## previous_cancellations      7.52004    0.05566 135.116  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.154 on 54392 degrees of freedom
## Multiple R-squared:  0.2568, Adjusted R-squared:  0.2568 
## F-statistic:  6265 on 3 and 54392 DF,  p-value: < 2.2e-16

Tu vidíme, že vynechanie premennej previous_cancellations podstatne znižuje upravený koeficient determinácie, a preto v prípade vynechávania premennej by sme uprednostnili vynechanie premennej previous_bookings_not_canceled.

Škálovanie premenných

Škálovanie môže byť veľmi efektívne, znižuje ale interpretovateľnosť modelu. Ide o úpravu premenných podľa nasledovného vzorca:

\[x^{scale} = \frac{x-M}{\sqrt{D}}\] kde \(M\) je stredná hodnota (priemer) a \(D\) je rozptyl premennej.

udaje.new$arrival_date_week_number_c <- scale(udaje.new$arrival_date_week_number, center=TRUE, scale=TRUE)
udaje.new$arrival_date_day_of_month_c <- scale(udaje.new$arrival_date_day_of_month, center=TRUE, scale=TRUE)
udaje.new$previous_cancellations_c <- scale(udaje.new$previous_cancellations, center=TRUE, scale=TRUE)
udaje.new$previous_bookings_not_canceled_c <- scale(udaje.new$previous_bookings_not_canceled, center=TRUE, scale=TRUE)

model_centered <- lm(
  lead_time ~ arrival_date_week_number_c + arrival_date_day_of_month_c +
    previous_cancellations_c + previous_bookings_not_canceled_c,
  data = udaje.new,
  na.action = na.exclude
)
summary(model_centered)
## 
## Call:
## lm(formula = lead_time ~ arrival_date_week_number_c + arrival_date_day_of_month_c + 
##     previous_cancellations_c + previous_bookings_not_canceled_c, 
##     data = udaje.new, na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -86.222  -0.032  -0.032  -0.032 294.449 
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       5.14196    0.01728 297.567   <2e-16 ***
## arrival_date_week_number_c        0.88736    0.01820  48.747   <2e-16 ***
## arrival_date_day_of_month_c      -0.02906    0.01773  -1.639    0.101    
## previous_cancellations_c          2.88612    0.01913 150.879   <2e-16 ***
## previous_bookings_not_canceled_c -1.11205    0.01912 -58.150   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.03 on 54391 degrees of freedom
## Multiple R-squared:  0.3003, Adjusted R-squared:  0.3002 
## F-statistic:  5836 on 4 and 54391 DF,  p-value: < 2.2e-16
vif(model_centered)
##       arrival_date_week_number_c      arrival_date_day_of_month_c 
##                         1.109678                         1.052978 
##         previous_cancellations_c previous_bookings_not_canceled_c 
##                         1.225388                         1.224777

Conditional Number je

X <- model.matrix(model_centered)[, -1]
XtX <- t(X) %*% X
eig <- eigen(XtX)

condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number
## [1] 1.602874

Vidíme, že Condition number sa výrazne zlepšilo a momentálne vykazuje znak minimálnej až nulovej multikolinearity.

10. Záver

Náš model pôvodne vykazoval znaky multikolinearity. Model vieme vylepšiť vyhodením premennej previous_bookings_not_canceled, ktorá hovorí o počte nezrušených rezervácií daného zákazníka. Táto premenná by síce na prvý pohľad mala byť pre model rovnako relevantná, ako premenná previous_cancellations, ktorá hovorí o počte zrušených rezervácií toho istého zákazníka. Z našich výsledkov však vidíme, že zrušené rezervácie sú pre model oveľa výpovednejšie.

Náš preškálovaný model dosahuje veľmi dobré výsledky, ale jeho interpretovateľnosť sa znižuje. Overíme si teda na záver multikolinearitu nášho modelu bez premennej previous_bookings_not_canceled pomocou condition number.

X <- model.matrix(model_no.previous_bookings_not_canceled)[, -1]
XtX <- t(X) %*% X
eig <- eigen(XtX)

condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number
## [1] 80.17795

Vidíme, že hodnota condition number oproti pôvodnému modelu poklesla, ale stály hovorí o silnej multikolinearite. V rámci nášho modelu by sme si teda radšej vybrali preškálovaný model model_centered.