Cvičenie 10 – Multikolinearita v regresných modeloch

1. Úvod

Po autokorelá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 predpokladá, ž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 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 prierezovým regresným modelom, kde pre každú obec uvažujeme počet obyvateľov v decembri 2024 ako závislú premennú a počty obyvateľov v predchádzajúcich mesiacoch toho istého roku ako vysvetľujúce premenné:

\[ \text{Pop12}_i = \beta_0 + \beta_1 \text{Pop11}_i + \beta_2 \text{Pop10}_i + \beta_3 \text{Pop09}_i + u_i, \]

kde:

• \(\text{Pop12}_i\) – počet obyvateľov obce \(i\) v decembri 2024,
  
• \(\text{Pop11}_i\) – počet obyvateľov v novembri 2024,
  
• \(\text{Pop10}_i\) – počet obyvateľov v októbri 2024,
  
• \(\text{Pop09}_i\) – počet obyvateľov v septembri 2024.

Údaje máme z databázy počtu obyvateľov slovenských obcí podľa mesiacov (rok 2024). Po načítaní ich uložíme do data.frame udaje a zároveň premenujeme stĺpce na kratšie názvy.

udaje <- read.csv("Databaza.csv",
                  sep = ";",
                  dec = ".",
                  header = TRUE,
                  stringsAsFactors = FALSE)

names(udaje)

 [1] "Okres"    "Obec"     "X2024M12" "X2024M11" "X2024M10" "X2024M09" "X2024M08"
 [8] "X2024M07" "X2024M06" "X2024M05" "X2024M04" "X2024M03" "X2024M02" "X2024M01"

# vyberieme si len stĺpce, ktoré potrebujeme
udaje <- udaje[, c("Okres", "Obec", "X2024M12", "X2024M11", "X2024M10", "X2024M09")]

# premenujeme stĺpce na jednoduchšie názvy
colnames(udaje) <- c("Okres", "Obec", "Pop12", "Pop11", "Pop10", "Pop09")

head(udaje)

V našej databáze nemáme chýbajúce hodnoty, takže nie je potrebná žiadna imputácia.

---

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

model <- lm(Pop12 ~ Pop11 + Pop10 + Pop09,
            data = udaje)
summary(model)

Call:
lm(formula = Pop12 ~ Pop11 + Pop10 + Pop09, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-44.442  -2.185   0.749   3.659  50.908 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.70441    0.68551  -2.486   0.0137 *  
Pop11        1.24863    0.06192  20.164  < 2e-16 ***
Pop10        0.09468    0.12167   0.778   0.4374    
Pop09       -0.34283    0.08212  -4.175 4.48e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.856 on 196 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 1.19e+08 on 3 and 196 DF,  p-value: < 2.2e-16

Vo výsledkoch zvyčajne vidíme veľmi vysoký koeficient determinácie a koeficienty pri Pop11, Pop10 a Pop09 sú si veľmi podobné. To je intuitívne, pretože počty obyvateľov v jednotlivých mesiacoch sa menia len veľmi málo a sú takmer perfektné lineárne kombinácie.

---

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("Pop11", "Pop10", "Pop09")]
round(cor(xvars), 3)

      Pop11 Pop10 Pop09
Pop11     1     1     1
Pop10     1     1     1
Pop09     1     1     1

V našom prípade uvidíme, že korelácie medzi mesačnými počtami obyvateľov sú prakticky rovné 1. To je extrémny príklad multikolinearity – všetky vysvetľujúce premenné obsahujú takmer tú istú informáciu.

---

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

pairs(xvars,
      main = "Scatterplotová matica – premenné Pop11, Pop10, Pop09")

---

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)

  Pop11   Pop10   Pop09 
1367712 5278985 2404844 

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šich dátach sú VIF-y typicky veľmi vysoké, čo zodpovedá takmer dokonalej korelácii medzi mesačnými hodnotami počtu obyvateľov.

---

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 nemusia byť stabilné a ich štandardné odchýlky môžu byť nafúknuté. 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ú 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. Condition number nie je test, je to len indikátor, ktorý posudzuje mieru multikolinearity medzi premennými. Používame intuitívne pravidlo. Ak Condition number je

• < 10 → nízka multikolinearita,
• 10–30 → mierna,
• 30–100 → silná,

• 100 → veľmi vážna.

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

Keďže u nás tento indikátor typicky výrazne presahuje hodnotu 100, signalizuje prítomnosť závažnej multikolinearity medzi mesačnými údajmi.

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_j\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 „jednej“ vysvetľujúcej veličiny, 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 (u nás sú VIF-y vysoké pre všetky), a porovnajme následne upravené koeficienty determinácie oboch nových modelov:

model_noPop10 <- lm(Pop12 ~ Pop11 + Pop09, data = udaje)
summary(model_noPop10)

Call:
lm(formula = Pop12 ~ Pop11 + Pop09, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-43.171  -2.310   0.847   3.530  50.940 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.68757    0.68448  -2.465   0.0145 *  
Pop11        1.28666    0.03798  33.876  < 2e-16 ***
Pop09       -0.28619    0.03799  -7.534 1.74e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.847 on 197 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 1.789e+08 on 2 and 197 DF,  p-value: < 2.2e-16

model_noPop11 <- lm(Pop12 ~ Pop10 + Pop09, data = udaje)
summary(model_noPop11)

Call:
lm(formula = Pop12 ~ Pop10 + Pop09, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-53.303  -4.884  -0.385   4.526  85.897 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.8557     1.1967  -0.715    0.475    
Pop10         2.0312     0.1307  15.546  < 2e-16 ***
Pop09        -1.0305     0.1307  -7.888 2.07e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.49 on 197 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 5.836e+07 on 2 and 197 DF,  p-value: < 2.2e-16

V našej databáze zvyčajne uvidíme, že vynechanie jednej z veľmi podobných mesačných premenných zníži upravený koeficient determinácie len minimálne. To je typické správanie pri silnej multikolinearite – jedna premenná dokáže takmer dokonale „zastúpiť“ inú.

Škálovanie premenných

Škálovanie môže byť veľmi efektívne z hľadiska numerickej stability (najmä pri veľmi rozdielnych rádoch premenných), znižuje však 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.

V našej databáze sú všetky premenné v počte obyvateľov, teda už sú v porovnateľných jednotkách. Napriek tomu si škálovanie ukážeme, aby sme videli vplyv na multikolinearitu.

udaje$Pop11_c <- scale(udaje$Pop11, center = TRUE, scale = TRUE)
udaje$Pop10_c <- scale(udaje$Pop10, center = TRUE, scale = TRUE)
udaje$Pop09_c <- scale(udaje$Pop09, center = TRUE, scale = TRUE)

model_centered <- lm(Pop12 ~ Pop11_c + Pop10_c + Pop09_c,
                     data = udaje)
summary(model_centered)

Call:
lm(formula = Pop12 ~ Pop11_c + Pop10_c + Pop09_c, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-44.442  -2.185   0.749   3.659  50.908 

Coefficients:
              Estimate Std. Error  t value Pr(>|t|)    
(Intercept)  4930.9000     0.6262 7874.461  < 2e-16 ***
Pop11_c     14803.9118   734.1607   20.164  < 2e-16 ***
Pop10_c      1122.4034  1442.3446    0.778    0.437    
Pop09_c     -4064.1462   973.5030   -4.175 4.48e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.856 on 196 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 1.19e+08 on 3 and 196 DF,  p-value: < 2.2e-16

vif(model_centered)

Pop11_c Pop10_c Pop09_c 
1367712 5278985 2404844 

Condition 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] 4909.792

Z výsledkov vidíme, že škálovanie môže zlepšiť numerické vlastnosti matice (Condition Number), ale multikolinearitu medzi mesačnými hodnotami ako takú neodstraňuje – dáta stále obsahujú takmer tú istú informáciu.

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

Ak sa chceme vyhnúť strate interpretovateľnosti, môžeme niekedy zmeniť jednotky premennej, ktorá je v úplne iných rádoch než ostatné. V pôvodnom príklade to bola premenná GDP, ktorá bola v dolároch. V našej databáze sú všetky premenné v počte obyvateľov, takže tento problém nemáme. Pre ilustráciu si však ukážeme prevod jednej premennej na „tisíce obyvateľov“:

udaje$Pop11k <- udaje$Pop11 / 1000

head(udaje[, c("Obec", "Pop12", "Pop11", "Pop11k")])

Potom lineárny model dosiahne výsledky:

model_Pop11k <- lm(Pop12 ~ Pop11k + Pop10 + Pop09,
                   data = udaje)
summary(model_Pop11k)

Call:
lm(formula = Pop12 ~ Pop11k + Pop10 + Pop09, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-44.442  -2.185   0.749   3.659  50.908 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -1.70441    0.68551  -2.486   0.0137 *  
Pop11k      1248.62877   61.92243  20.164  < 2e-16 ***
Pop10          0.09468    0.12167   0.778   0.4374    
Pop09         -0.34283    0.08212  -4.175 4.48e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.856 on 196 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 1.19e+08 on 3 and 196 DF,  p-value: < 2.2e-16

vif(model_Pop11k)

 Pop11k   Pop10   Pop09 
1367712 5278985 2404844 

V tomto prípade sa koeficient pri Pop11k zmení len v jednotkách (na tisíce obyvateľov), ale samotná multikolinearita nezmizne, pretože Pop11, Pop10 a Pop09 sú stále veľmi podobné.

Poznámka: Dummy premenné (0/1) neškálujeme.

###   očistenie databázy od „pracovných“ stĺpcov
library(dplyr)

udaje <- udaje %>%
  dplyr::select(Okres, Obec, Pop12, Pop11, Pop10, Pop09)

PCA (voliteľné)

PCA nám vytvára nový pracovný vektor (komponent), ktorý je definovaný ako vážený súčet dvoch alebo viacerých korelovaných premenných. Tento komponent potom vystupuje ako vysvetľujúca veličina, zatiaľ čo pôvodné premenné vylúčime zo zoznamu regresorov. Tým sa redukuje počet vysvetľujúcich veličín a môže sa zmierniť problém multikolinearity.

V našom prípade vytvoríme prvú hlavnú komponentu z premenných Pop11, Pop10 a Pop09:

X_pca <- scale(udaje[, c("Pop11", "Pop10", "Pop09")],
               center = TRUE, scale = TRUE)
pca_res <- prcomp(X_pca)

summary(pca_res)

Importance of components:
                         PC1      PC2       PC3
Standard deviation     1.732 0.000992 0.0003528
Proportion of Variance 1.000 0.000000 0.0000000
Cumulative Proportion  1.000 1.000000 1.0000000

udaje$PC1 <- pca_res$x[, 1]

model_pca <- lm(Pop12 ~ PC1, data = udaje)
summary(model_pca)

Call:
lm(formula = Pop12 ~ PC1, data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-80.322  -3.752  -0.201   3.493 112.824 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4930.9000     1.1468    4300   <2e-16 ***
PC1         -6848.6258     0.6638  -10318   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 16.22 on 198 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 1.065e+08 on 1 and 198 DF,  p-value: < 2.2e-16

# VIF sa nedá počítať, lebo model má len jeden vysvetľujúci regresor
# vif(model_pca)

---

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\), aby sme získali určitú predstavu o číselnom ráde prvkov:

library(MASS)

X <- as.matrix(udaje[, c("Pop11", "Pop10", "Pop09")])
y <- udaje$Pop12

ridge_fit <- lm.ridge(y ~ X, lambda = seq(0, 10, 1))
ridge_fit

                XPop11     XPop10     XPop09
 0 -1.704409 1.2486288 0.09468307 -0.3428288
 1  8.529458 0.3331203 0.33297224  0.3328418
 2 16.706357 0.3324811 0.33243251  0.3323613
 3 24.856002 0.3319011 0.33188566  0.3318343
 4 32.978628 0.3313373 0.33133840  0.3312969
 5 41.074389 0.3307811 0.33079207  0.3307565
 6 49.143424 0.3302296 0.33024709  0.3302155
 7 57.185869 0.3296815 0.32970366  0.3296750
 8 65.201856 0.3291362 0.32916186  0.3291353
 9 73.191517 0.3285934 0.32862173  0.3285968
10 81.154981 0.3280529 0.32808330  0.3280597

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

---

10. Zhrnutie

• Multikolinearita nezavádza bias, ale zvyšuje štandardné odchýlky odhadovaných regresných koeficientov a robí ich nestabilnými.
  
• Testy (resp. diagnostické ukazovatele) ako VIF a Condition Number umožňujú diagnostiku multikolinearity.
  
• Riešenia zahŕňajú: vynechanie premennej, centrovanie/škálovanie, zmenu jednotiek, prípadne použitie PCA alebo hrebeňovej regresie. PCA a hrebeňová regresia vyžadujú hlbšie vedomosti, preto ich tu uvádzame len formálne a nevyžadujeme ich detailný výpočet.

---