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

\[ Life.expectancy_i = \beta_0 + \beta_1 BMI_i + \beta_2 GDP_i + \beta_3 Schooling_i + u_i \]

Prierezové údaje máme z databázy WHO a po ich nadčítaní sme ich uložili do data.frame udaje.

udaje <- read.csv("udaje/Life_Expectancy_Data.csv",dec=".",sep=",",header = TRUE)
# select just the record from 2015
udaje.2015 <- udaje[udaje$Year==2015,c("Life.expectancy","BMI","GDP","Schooling")]

# data imputation

# Compute column medians
#column_medians <- sapply(udaje.2015, median, na.rm = TRUE)

# Impute missing values with column medians
# Compute column medians
column_medians <- sapply(udaje.2015, median, na.rm = TRUE)

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

udaje.2015 <- udaje_imputed
udaje <- udaje.2015

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

model <- lm(Life.expectancy ~ BMI + GDP + Schooling,
            data = udaje)
summary(model)
## 
## Call:
## lm(formula = Life.expectancy ~ BMI + GDP + Schooling, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.6471  -2.4024   0.4563   3.2355  12.8591 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.419e+01  1.839e+00  24.033   <2e-16 ***
## BMI         4.948e-02  2.144e-02   2.308   0.0222 *  
## GDP         6.972e-05  3.845e-05   1.813   0.0715 .  
## Schooling   1.921e+00  1.645e-01  11.676   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.033 on 179 degrees of freedom
## Multiple R-squared:  0.6224, Adjusted R-squared:  0.6161 
## F-statistic: 98.36 on 3 and 179 DF,  p-value: < 2.2e-16

Vo výsledkoch pozorujeme veľké rádové rozdiely v odhade regresných koeficientov - porovnajte si rády koeficientov pri GDP a pri Schooling.

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("BMI", "GDP", "Schooling")]
round(cor(xvars), 3)
##             BMI   GDP Schooling
## BMI       1.000 0.311     0.526
## GDP       0.311 1.000     0.435
## Schooling 0.526 0.435     1.000

V našom prípade teda nevidíme dve premenné ktoré by boli silne korelované.

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 – premenné BMI, GDP, Schooling")

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)
##       BMI       GDP Schooling 
##  1.398736  1.247834  1.558912

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] 2692.719

Keďže u nás tento indikátor presahuje významne hodnotu 100, signalizuje prítomnosť závažnej multikolinearity.

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.

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

model_noBMI <- lm(Life.expectancy ~ GDP + Schooling, data = udaje)
summary(model_noBMI)
## 
## Call:
## lm(formula = Life.expectancy ~ GDP + Schooling, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.3746  -2.4939   0.5427   3.3179  10.8549 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4.401e+01  1.859e+00  23.672   <2e-16 ***
## GDP         7.920e-05  3.869e-05   2.047   0.0421 *  
## Schooling   2.094e+00  1.481e-01  14.142   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.094 on 180 degrees of freedom
## Multiple R-squared:  0.6112, Adjusted R-squared:  0.6069 
## F-statistic: 141.5 on 2 and 180 DF,  p-value: < 2.2e-16
model_noSchooling <- lm(Life.expectancy ~ GDP + BMI, data = udaje)
summary(model_noSchooling)
## 
## Call:
## lm(formula = Life.expectancy ~ GDP + BMI, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.2929  -4.1549   0.4721   4.0853  23.2326 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.317e+01  1.140e+00  55.429  < 2e-16 ***
## GDP         2.206e-04  4.793e-05   4.603 7.85e-06 ***
## BMI         1.638e-01  2.524e-02   6.490 8.06e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.662 on 180 degrees of freedom
## Multiple R-squared:  0.3349, Adjusted R-squared:  0.3275 
## F-statistic: 45.31 on 2 and 180 DF,  p-value: < 2.2e-16

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

Š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$BMI_c       <- scale(udaje$BMI, center=TRUE, scale=TRUE)
udaje$GDP_c       <- scale(udaje$GDP, center=TRUE, scale=TRUE)
udaje$Schooling_c <- scale(udaje$Schooling, center=TRUE, scale=TRUE)

model_centered <- lm(Life.expectancy ~ BMI_c + GDP_c + Schooling_c,
                     data = udaje)
summary(model_centered)
## 
## Call:
## lm(formula = Life.expectancy ~ BMI_c + GDP_c + Schooling_c, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.6471  -2.4024   0.4563   3.2355  12.8591 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  71.6169     0.3721 192.478   <2e-16 ***
## BMI_c         1.0183     0.4413   2.308   0.0222 *  
## GDP_c         0.7557     0.4168   1.813   0.0715 .  
## Schooling_c   5.4392     0.4658  11.676   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.033 on 179 degrees of freedom
## Multiple R-squared:  0.6224, Adjusted R-squared:  0.6161 
## F-statistic: 98.36 on 3 and 179 DF,  p-value: < 2.2e-16
vif(model_centered)
##       BMI_c       GDP_c Schooling_c 
##    1.398736    1.247834    1.558912

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] 2.0353

Z výsledkov vidíme, že ukazovateľ Condition number sa podstatne zlepšil a nesignalizuje už žiadnu multikolinearitu. Ukazovatele VIF aj naďalej ostali veľmi nízke, čo je dobrou vlastnosťou preškálovania premenných

Iná úprava premennej, ktorá zachová interpretovateľnosť

Ak sa chceme vynhúť strate interpretovateľnosti, môžeme sa pozrieť, v akých rádoch sa pohybuje vývoj jednotlivých vysvetľujúcich veličín a upraviť jednoduchým prevodom na iné jednotky tú premennú, ktorá sa odlišuje od ostatných. V našom prípade ide o premennú GDP, ktorá meria HDP na obyvateľa v dolároch. My to predelením 1000 prevedieme na HDP na obyvateľa vyjadrený v tisícoch dolárov. Tým dostaneme premenné vyjadrené v rádovo porovnateľných jednotkách tak, ako je to uvedené v nasledovnej tabuľke.

udaje$GDP1000 = udaje$GDP/1000
head(udaje)
##    Life.expectancy  BMI        GDP Schooling      BMI_c      GDP_c Schooling_c
## 1             65.0 19.1   584.2592      10.1 -1.1524374 -0.5466216  -1.0018056
## 17            77.8 58.0  3954.2278      14.2  0.7376778 -0.2356953   0.4461905
## 33            75.6 59.5  4132.7629      14.4  0.8105614 -0.2192230   0.5168244
## 49            52.4 23.3  3695.7937      11.4 -0.9483632 -0.2595394  -0.5426849
## 65            76.4 47.7 13566.9541      13.9  0.2372103  0.6512118   0.3402395
## 81            76.3 62.8 13467.1236      17.3  0.9709054  0.6420010   1.5410168
##       GDP1000
## 1   0.5842592
## 17  3.9542278
## 33  4.1327629
## 49  3.6957937
## 65 13.5669541
## 81 13.4671236

Potom lineárny model dosiahne výsledky

model_GDP_1000 <- lm(Life.expectancy ~ BMI + GDP1000 + Schooling,
                     data = udaje)
summary(model_GDP_1000)
## 
## Call:
## lm(formula = Life.expectancy ~ BMI + GDP1000 + Schooling, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.6471  -2.4024   0.4563   3.2355  12.8591 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 44.19364    1.83886  24.033   <2e-16 ***
## BMI          0.04948    0.02144   2.308   0.0222 *  
## GDP1000      0.06972    0.03845   1.813   0.0715 .  
## Schooling    1.92097    0.16452  11.676   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.033 on 179 degrees of freedom
## Multiple R-squared:  0.6224, Adjusted R-squared:  0.6161 
## F-statistic: 98.36 on 3 and 179 DF,  p-value: < 2.2e-16
vif(model_GDP_1000)
##       BMI   GDP1000 Schooling 
##  1.398736  1.247834  1.558912

kde všetky regresné koeficienty majú porovnateľné rády a tiež VIF je akceptovateľný vo všetkých prípadoch. Conditional number je potom

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

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

čo sa nachádza niekde na hranici nášho intuitívneho kritéria medzi Nízka multikolinearita a Mierna multikolinearita. Regresné koeficienty si ale zachovali dobrú interpretovateľnosť.

Poznámka: Dummy premenné neškálujeme

###   ocistenie databazy od nadbytocnych - pracovnych stlpcov
library(dplyr)

udaje <- udaje %>%
  dplyr::select(-Schooling_c, -BMI_c, -GDP_c, GDP1000)

PCA (voliteľné)

PCA nám vytvára pracovný vektor, ktorý je definovaný ako vážený súťet dvoch nami špecifikovaných korelovaných premenných. Tento vektor potom vystupuje ako vysvetľujúca veličina, zatiaľ čo začlenené dve vysvetľujúce premenné vylúčime zo zoznamu regresorov. Tým sa redukuje počet vysvetľujúcich veličín a odstraňuje sa aj problém multikolinearity.

X_pca <- scale(udaje[, c("GDP", "Schooling")], center=TRUE, scale=TRUE)
pca_res <- prcomp(X_pca)

summary(pca_res)
## Importance of components:
##                           PC1    PC2
## Standard deviation     1.1980 0.7516
## Proportion of Variance 0.7176 0.2824
## Cumulative Proportion  0.7176 1.0000
udaje$PC1 <- pca_res$x[,1]

model_pca <- lm(Life.expectancy ~ BMI + PC1, data = udaje)
summary(model_pca)
## 
## Call:
## lm(formula = Life.expectancy ~ BMI + PC1, data = udaje)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.3086  -3.3559   0.7001   3.9098  16.0834 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 68.11003    1.07188  63.542  < 2e-16 ***
## BMI          0.08190    0.02311   3.545 0.000501 ***
## PC1          4.10533    0.39696  10.342  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.578 on 180 degrees of freedom
## Multiple R-squared:  0.5337, Adjusted R-squared:  0.5285 
## F-statistic:   103 on 2 and 180 DF,  p-value: < 2.2e-16
vif(model_pca)
##      BMI      PC1 
## 1.322726 1.322726

Hrebeňová regresia (Ridge Regression - voliteľné)

Hrebeňová regresia je modifikácia regresie, ktorá zavádza perturbácie do matice \(\mathbf X^T \mathbf X\) tak, aby znížila dôsledky multikolinearity. Treba ale upozorniť, že odhadované regresné koeficienty sú skreslené. Perturbácia vyzerá nasledovne

\[(\mathbf X^T \mathbf X + \lambda \mathbf I)\]

V tomto kroku si vypíšeme výsledky odhadov regresných koeficientov s meniacimi sa parametrami \(lambda\) tak, aby sme získlali určitú predstavu o číselnom ráde prvkov, ktoré obsahuje. Na základe toho sa potom rozhodneme pre nejakú voľbu parametra \(\lambda\).

Nastavovanie rôznych hodnôt \(\lambda\) je predmetom hlbšej analýzy, tento prehľad je len jej časťou.

ridge_fit <- lm.ridge(y ~ X, lambda = seq(0, 10, 1))
ridge_fit
##                   XBMI         XGDP XSchooling
##  0 44.19364 0.04947912 6.972179e-05   1.920974
##  1 44.34557 0.05009188 7.058986e-05   1.906765
##  2 44.49495 0.05068217 7.142956e-05   1.892841
##  3 44.64187 0.05125088 7.224193e-05   1.879194
##  4 44.78640 0.05179884 7.302793e-05   1.865813
##  5 44.92861 0.05232686 7.378848e-05   1.852689
##  6 45.06857 0.05283570 7.452448e-05   1.839815
##  7 45.20636 0.05332607 7.523678e-05   1.827183
##  8 45.34203 0.05379867 7.592618e-05   1.814785
##  9 45.47564 0.05425417 7.659347e-05   1.802613
## 10 45.60726 0.05469318 7.723938e-05   1.790661

10. Zhrnutie