Multikolinearita v regresných modeloch

1. Východiskový model a údaje

Budeme pracovať s regresným modelom, ktorý sme definovali v predchádzajúcich cvičeniach, pričom skúmame faktory ovplyvňujúce zdravie:

\[\text{Healthy life expectancy at birth}_i = \beta_0 + \beta_1 \left(\text{Social support}\right)_i + \beta_2 \left(\text{Log GDP per capita}\right)_i + \beta_3 \left(\text{Freedom to make life choices}\right)_i + u_i\] Prierezové údaje pochádzajú z World Happiness Report 2005-2021 a sú filtrované pre rok 2015.

library(readr)  
library(dplyr)  
library(tidyr) 

udaje_all <- read_csv("World Happiness Report 2005-2021.csv")

udaje.2015 <- udaje_all %>%
  filter(Year == 2015) %>%

  dplyr::select(
    `Healthy life expectancy at birth`, 
    `Log GDP per capita`,            
    `Social support`,                 
    `Freedom to make life choices`    ) %>% 
  
  mutate(across(
    .cols = everything(), 
    .fns = ~tidyr::replace_na(., median(., na.rm = TRUE)) ))

udaje <- udaje.2015

print(head(udaje))

print(summary(udaje))

 Healthy life expectancy at birth Log GDP per capita Social support  
 Min.   :48.10                    Min.   : 6.742     Min.   :0.4344  
 1st Qu.:58.90                    1st Qu.: 8.507     1st Qu.:0.7292  
 Median :65.10                    Median : 9.444     Median :0.8255  
 Mean   :63.44                    Mean   : 9.390     Mean   :0.7984  
 3rd Qu.:68.30                    3rd Qu.:10.338     3rd Qu.:0.8996  
 Max.   :73.60                    Max.   :11.637     Max.   :0.9873  
 Freedom to make life choices
 Min.   :0.3889              
 1st Qu.:0.6577              
 Median :0.7759              
 Mean   :0.7493              
 3rd Qu.:0.8498              
 Max.   :0.9799              

Vďaka úspešnému spusteniu funkcií v R kóde sú dáta teraz pripravené na štatistickú analýzu (napr. regresiu). Všetky chýbajúce hodnoty (NA) v týchto stĺpcoch boli nahradené mediánmi stĺpcov vypočítanými z roku 2015.

---

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

model <- lm(`Healthy life expectancy at birth` ~ 
              `Log GDP per capita` + 
              `Social support` + 
              `Freedom to make life choices`,
            data = udaje)

summary(model)

Call:
lm(formula = `Healthy life expectancy at birth` ~ `Log GDP per capita` + 
    `Social support` + `Freedom to make life choices`, data = udaje)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.7871  -1.3770   0.2173   2.4253   6.6278 

Coefficients:
                               Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     18.1878     2.4714   7.359 1.47e-11 ***
`Log GDP per capita`             4.2031     0.3327  12.632  < 2e-16 ***
`Social support`                 6.1845     3.1527   1.962   0.0518 .  
`Freedom to make life choices`   1.1288     2.2757   0.496   0.6207    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.376 on 139 degrees of freedom
Multiple R-squared:  0.7306,    Adjusted R-squared:  0.7248 
F-statistic: 125.6 on 3 and 139 DF,  p-value: < 2.2e-16

Model vykazuje veľmi silnú vysvetľovaciu silu, keďže spoločne vybrané nezávislé premenné vysvetľujú približne 72.48% variability v zdravej dĺžke života (Healthy life expectancy at birth), čo naznačuje upravený koeficient determinácie (\(R^2 = 0.7248\)). Model je tiež štatisticky vysoko významný ako celok (nízka \(p\)-hodnota \(F\)-štatistiky).

Pri pohľade na individuálne prediktory zistíme, že:

• Ekonomický faktor (Log GDP per capita) je najsilnejším a vysoko štatisticky významným prediktorom. So zvýšením Log HDP na obyvateľa o jednu jednotku sa zdravá dĺžka života pri narodení v priemere zvýši o 4.20 roka (pri fixných ostatných faktoroch), čo je koeficient s najvyššou signifikanciou (\(p < 2\times 10^{-16}\)).

• Sociálna podpora (Social support) má silný, ale len hranične štatisticky významný pozitívny vplyv. Zvýšenie sociálnej podpory o jednu jednotku koreluje so zvýšením zdravej dĺžky života o 6.18 roka, avšak tento vplyv je len na hranici prijatej 5% hladiny významnosti (\(p = 0.0518\)).

• Sloboda (Freedom to make life choices) sa v tomto modeli nepreukázala ako štatisticky významný prediktor zdravej dĺžky života (\(p = 0.6207\)), čo naznačuje, že jej vplyv je už pravdepodobne pokrytý silnejšími ekonomickými a sociálnymi faktormi.

---

3. Korelačná matica

xvars_whr <- udaje %>%
  dplyr::select(`Log GDP per capita`, 
                `Social support`, 
                `Freedom to make life choices`)

cor_matrix <- cor(xvars_whr)
round(cor_matrix, 3)

                             Log GDP per capita Social support
Log GDP per capita                        1.000          0.678
Social support                            0.678          1.000
Freedom to make life choices              0.323          0.376
                             Freedom to make life choices
Log GDP per capita                                  0.323
Social support                                      0.376
Freedom to make life choices                        1.000

Najsilnejšia korelácia existuje medzi Log GDP per capita a Social support (\(r = 0.678\)). Táto hodnota je pomerne vysoká a naznačuje potenciálny problém multikolinearity v regresnom modeli. To znamená, že tieto dve premenné merajú do istej miery podobnú informáciu alebo sú silne spojené.

Ostatné vzťahy (Log GDP/Freedom a Social support/Freedom) sú na strednej alebo slabej úrovni (\(r \approx 0.3-0.4\)), čo obvykle nespôsobuje vážne problémy s multikolinearitou.

---

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_whr,
      main = "Scatterplotová matica – vplyvové premenné (WHR 2015)")

Scatterplotová matica vizuálne potvrdzuje vzájomné vzťahy medzi nezávislými premennými, ktoré boli predtým kvantifikované pomocou korelačnej matice.

• Vzťah medzi Log GDP per capita a Social support:

Silná pozitívna korelácia: Body tvoria jasný, smerom nahor stúpajúci oblak, čo vizuálne potvrdzuje silnú pozitívnu koreláciu (\(r=0.678\)). V krajinách s vyššou úrovňou HDP na obyvateľa je takmer vždy vyššia aj miera sociálnej podpory. Táto vizualizácia najzreteľnejšie poukazuje na potenciálnu multikolinearitu.

• Vzťah medzi Log GDP per capita a Freedom to make life choices (Stredný ľavý panel):

Slabá/Stredná korelácia: Body sú rozptýlenejšie, hoci je viditeľný mierny pozitívny sklon. Vplyv HDP na mieru slobody je menej silný a nie je tak konzistentný ako v prípade sociálnej podpory.

• Vzťah medzi Social support a Freedom to make life choices (Spodný stredný panel):

Stredná korelácia: Je badateľný pozitívny vzťah, ale rozptyl dát je pomerne veľký. Vyššia sociálna podpora je mierne spojená s vyššou mierou slobody, ale vzťah nie je dostatočne silný na to, aby predstavoval vážny problém multikolinearity.

---

4. VIF

library(car)
vif_result <- car::vif(model)

print(vif_result)

          `Log GDP per capita`               `Social support` 
                      1.867339                       1.948592 
`Freedom to make life choices` 
                      1.176299 

Napriek silnej korelácii \(r = 0.678\) medzi Log GDP per capita a Social support sa žiadna z premenných nepreukázala ako vážne multikolineárna.

Všetky hodnoty VIF sú výrazne pod kritickou hranicou 5. To naznačuje, že aj keď sú tieto premenné korelované, ich vzájomné prekrývanie nie je dostatočne silné na to, aby spôsobovalo nestabilitu odhadov regresných koeficientov alebo výrazné zvýšenie ich štandardných chýb. Z tohto pohľadu nie je nutné podnikať kroky na odstránenie multikolinearity, no za účelom demonštrácie ich vykonáme.

---

5. Condition Number

X <- model.matrix(model)[, -1]
XtX <- t(X) %*% X
eig <- eigen(XtX)
condition_number <- sqrt(max(eig$values) / min(eig$values))
condition_number

[1] 109.1818

Číslo podmienenosti s hodnotou \(109.1818\) je extrémne vysoké (nad kritickou hranicou 30). Táto vysoká hodnota signalizuje, že matica je zle podmienená a odhady koeficientov sú nestabilné a mimoriadne citlivé na malé zmeny v dátach, čo potvrdzuje existenciu závažnej multikolinearity v dátovej štruktúre, ktorú VIF nedokázal plne zachytiť.

---

6. Riešenia multikolinearity

Vynechanie premennej

model_noSupport <- lm(`Healthy life expectancy at birth` ~ 
                        `Log GDP per capita` + 
                        `Freedom to make life choices`, 
                      data = udaje)
summary(model_noSupport)

Call:
lm(formula = `Healthy life expectancy at birth` ~ `Log GDP per capita` + 
    `Freedom to make life choices`, data = udaje)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.8797  -1.5891   0.4567   2.3869   6.1292 

Coefficients:
                               Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     18.4814     2.4918   7.417 1.05e-11 ***
`Log GDP per capita`             4.6171     0.2598  17.769  < 2e-16 ***
`Freedom to make life choices`   2.1386     2.2392   0.955    0.341    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.411 on 140 degrees of freedom
Multiple R-squared:  0.7231,    Adjusted R-squared:  0.7192 
F-statistic: 182.8 on 2 and 140 DF,  p-value: < 2.2e-16

model_noFreedom <- lm(`Healthy life expectancy at birth` ~ 
                        `Log GDP per capita` + 
                        `Social support`, 
                      data = udaje)
summary(model_noFreedom)

Call:
lm(formula = `Healthy life expectancy at birth` ~ `Log GDP per capita` + 
    `Social support`, data = udaje)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.6750  -1.5098   0.2003   2.3566   6.6329 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)           18.5971     2.3233   8.005 4.14e-13 ***
`Log GDP per capita`   4.2195     0.3302  12.779  < 2e-16 ***
`Social support`       6.5383     3.0627   2.135   0.0345 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.367 on 140 degrees of freedom
Multiple R-squared:  0.7301,    Adjusted R-squared:  0.7262 
F-statistic: 189.4 on 2 and 140 DF,  p-value: < 2.2e-16

Model je bez premennej Sloboda (model_noFreedom) najlepší a najekonomickejší, pretože dosahuje najvyššie upravené \(R^2\) a všetky zahrnuté prediktory (Log GDP per capita a Social support) sú štatisticky významné. Premenná Freedom to make life choices neprispieva k vysvetleniu variability a je optimálne ju z modelu vylúčiť.

Škálovanie premenných

library(car)
library(dplyr) 

udaje <- udaje %>%
  mutate(
    Support_c  = scale(`Social support`, center = TRUE, scale = TRUE),
    GDP_c      = scale(`Log GDP per capita`, center = TRUE, scale = TRUE),
    Freedom_c  = scale(`Freedom to make life choices`, center = TRUE, scale = TRUE))

model_centered <- lm(`Healthy life expectancy at birth` ~ 
                       Support_c + 
                       GDP_c + 
                       Freedom_c,
                     data = udaje)

summary(model_centered)

Call:
lm(formula = `Healthy life expectancy at birth` ~ Support_c + 
    GDP_c + Freedom_c, data = udaje)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.7871  -1.3770   0.2173   2.4253   6.6278 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  63.4406     0.2823 224.688   <2e-16 ***
Support_c     0.7759     0.3955   1.962   0.0518 .  
GDP_c         4.8911     0.3872  12.632   <2e-16 ***
Freedom_c     0.1524     0.3073   0.496   0.6207    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.376 on 139 degrees of freedom
Multiple R-squared:  0.7306,    Adjusted R-squared:  0.7248 
F-statistic: 125.6 on 3 and 139 DF,  p-value: < 2.2e-16

cat("\n--- VIF pre centrovaný model ---\n")

--- VIF pre centrovaný model ---

vif(model_centered)

Support_c     GDP_c Freedom_c 
 1.948592  1.867339  1.176299 

Štandardizácia premenných neovplyvnila celkovú silu modelu, keďže upravené \(R^2\) zostalo nezmenené na hodnote 0.7248.

Aj signifikancia prediktorov zostala rovnaká: Log HDP je vysoko významný, Sociálna podpora je hranične významná a Sloboda je nevýznamná.

Vďaka centrálnej zmene je teraz Intercept ľahko interpretovateľný ako priemerná zdravá dĺžka života (63.44 roka) v prípade, že všetky nezávislé premenné sú na svojej priemernej hodnote.

Aj napriek štandardizácii zostali hodnoty VIF rovnaké (max. 1.949), čo potvrdzuje, že centrácia nerieši mieru korelácie medzi nezávislými premennými.

Numerická stabilita sa však vyhodnotí na základe zmeny čísla podmienenosti v ďalšom kroku:

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.462653

Hodnota čísla podmienenosti \(2.46\) spadá do rozsahu slabého stupňa multikolinearity (pod prahom 10).

To v praxi znamená, že centrácia premenných bola úspešná. Pôvodný vážny problém so zlou podmienenosťou matice \(X^TX\) (pôvodná hodnota bola 109.18) bol plne eliminovaný.Tento centrovaný regresný model je teraz numericky stabilný.

library(dplyr)

udaje <- udaje %>%
  dplyr::select(!ends_with("_c"))
  
head(udaje)

print(colnames(udaje))

[1] "Healthy life expectancy at birth" "Log GDP per capita"              
[3] "Social support"                   "Freedom to make life choices"    
[5] "Log GDP per capita_1000"         

Dátová sada udaje bola úspešne vrátená do “čistého stavu”. Obsahuje len pôvodné prierezové dáta z World Happiness Report pre rok 2015 a je pripravená na ďalšie modely, napríklad na analýzu hlavných komponentov (PCA), ktorá bude nasledovať.

PCA

library(car) 
library(dplyr) 

X_pca <- udaje %>%
  dplyr::select(`Log GDP per capita`, `Freedom to make life choices`) %>%
  scale(center = TRUE, scale = TRUE) 

pca_res <- prcomp(X_pca)
summary(pca_res)

Importance of components:
                          PC1    PC2
Standard deviation     1.1500 0.8231
Proportion of Variance 0.6613 0.3387
Cumulative Proportion  0.6613 1.0000

udaje$PC1 <- pca_res$x[, 1]
model_pca <- lm(`Healthy life expectancy at birth` ~ 
                  `Social support` + 
                  PC1, 
                data = udaje)

summary(model_pca)

Call:
lm(formula = `Healthy life expectancy at birth` ~ `Social support` + 
    PC1, data = udaje)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.9209  -2.5779   0.2318   2.6854   8.9936 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)        50.667      3.007  16.849  < 2e-16 ***
`Social support`   15.999      3.740   4.278 3.48e-05 ***
PC1                 2.873      0.408   7.041 7.86e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.259 on 140 degrees of freedom
Multiple R-squared:  0.5682,    Adjusted R-squared:  0.562 
F-statistic: 92.11 on 2 and 140 DF,  p-value: < 2.2e-16

cat("\n--- VIF pre PCA model ---\n")

--- VIF pre PCA model ---

vif_result_pca <- car::vif(model_pca)
print(vif_result_pca)

`Social support`              PC1 
        1.723204         1.723204 

Analýza hlavných komponentov (PCA) bola vykonaná na premenných Log GDP per capita a Freedom to make life choices, s cieľom nahradiť ich jedným ortogonálnym faktorom (PC1) v modeli a zlepšiť tak multikolinearitu.

Prvý hlavný komponent (PC1) vysvetľuje 66.13% celkovej variability pôvodných dvoch premenných - je dostatočne silný na to, aby bol použitý ako efektívny zástupca oboch pôvodných, vzájomne korelovaných premenných.

Model (model_pca) predpovedá zdravú dĺžku života pomocou Social support a PC1.

• Upravené \(R^2\) je \(0.562\). Oproti pôvodnému modelu (\(R^2=0.7248\)) klesla vysvetľujúca sila modelu o približne 16 percentuálnych bodov. K poklesu došlo, pretože PC1 neuchováva celú informáciu pôvodných dvoch premenných (Log GDP a Freedom), ale len 66.13% z ich spoločnej variability.

• Oba prediktory sú vysoko štatisticky významné (\(p < 0.001\)). Zvýšenie sociálnej podpory o jednotku vedie k zvýšeniu zdravej dĺžky života o 15.999 roka, čo je koeficient s výrazne vyššou silou ako v pôvodnom modeli.

Hodnoty VIF sú pre obe premenné 1.7232.

• Hodnoty sú nízke, a čo je dôležité, sú pre obe premenné identické, pretože PC1 je ortogonálny k pôvodným premenným. Použitie PCA výrazne znížilo multikolinearitu v rámci nezávislých premenných.

---

Hrebeňová regresia (Ridge Regression)

library(MASS)
library(dplyr) 

X_vars <- udaje %>%
  dplyr::select(`Social support`, `Log GDP per capita`, `Freedom to make life choices`)
ridge_fit <- MASS::lm.ridge(
  `Healthy life expectancy at birth` ~ `Social support` + `Log GDP per capita` + `Freedom to make life choices`,
  data = udaje,
  lambda = seq(0, 10, 1) )
print(ridge_fit)

            `Social support` `Log GDP per capita`
 0 18.18784         6.184546             4.203129
 1 18.41699         6.426066             4.154874
 2 18.64172         6.654264             4.108264
 3 18.86223         6.870008             4.063202
 4 19.07874         7.074097             4.019600
 5 19.29142         7.267265             3.977380
 6 19.50046         7.450188             3.936465
 7 19.70600         7.623492             3.896786
 8 19.90821         7.787751             3.858280
 9 20.10722         7.943500             3.820887
10 20.30315         8.091233             3.784552
   `Freedom to make life choices`
 0                       1.128788
 1                       1.170350
 2                       1.211424
 3                       1.251981
 4                       1.291997
 5                       1.331451
 6                       1.370328
 7                       1.408613
 8                       1.446297
 9                       1.483371
10                       1.519831

–>

library(MASS)
library(dplyr) 

X_vars <- udaje %>%
  dplyr::select(`Social support`, `Log GDP per capita`, `Freedom to make life choices`)
ridge_fit <- MASS::lm.ridge(
  `Healthy life expectancy at birth` ~ `Social support` + `Log GDP per capita` + `Freedom to make life choices`,
  data = udaje,
  lambda = seq(0, 10, 1) )
print(ridge_fit)

Koeficienty pri \(\lambda = 0\)

• Keď je \(\lambda = 0\), Ridge Regresia je ekvivalentná klasickej OLS regresii (základný model).Koeficienty v tomto riadku (Intercepet: 18.18784, Log GDP: 4.203129, Social support: 6.184546) presne zodpovedajú koeficientom, ktoré sme získali v základnom modeli.

Vplyv Rastúceho \(\lambda\)

• S postupným rastom \(\lambda\) (od 1 do 10) sa uplatňuje penalizačný efekt Ridge Regresie, ktorý má za cieľ stabilizovať odhady koeficientov a riešiť multikolinearitu.

  • Log GDP per capita klesá (z 4.203 na 3.785)

  • Social support rastie (z 6.185 na 8.091)

  • Freedom to make life choices rastie (z 1.129 na 1.520)

  • Intercept rastie (z 18.188 na 20.303)

Zmeny koeficientov pre jednotlivé premenné s rastúcim \(\lambda\) ukazujú, že Ridge Regresia je úspešná v stabilizácii modelu, ktorý mal pôvodne “zlé” podmienené číslo (\(109.18\)). Ridge Regresia “ťahá” koeficienty bližšie k sebe, čím zmierňuje ich vzájomnú nestabilitu, ktorá bola spôsobená koreláciou medzi Log GDP per capita a Social support (\(r=0.678\)).

10. Zhrnutie

Analýza sa začala odhadom lineárneho modelu, ktorý preukázal veľmi silnú predikčnú silu s upraveným \(R^2\) na úrovni \(0.7248\). Z hľadiska individuálneho vplyvu bol Log GDP per capita vysoko významným prediktorom, zatiaľ čo Social support bol hranične významný a Freedom to make life choices sa ukázala ako nevýznamná premenná. Následná diagnostika odhalila problém s numerickou stabilitou: hoci hodnoty VIF boli nízke (max. \(1.949\)), kritické Číslo Podmienenosti \(109.18\) potvrdilo, že matica \(X^TX\) je zle podmienená a odhady koeficientov sú nestabilné.

Pre riešenie nestability a optimalizáciu modelu boli testované rôzne metódy. Najefektívnejším krokom bola centrácia a škálovanie premenných, ktorá drasticky znížila číslo podmienenosti na stabilnú hodnotu \(2.46\), čím sa úspešne eliminovala numerická nestabilita.

Z hľadiska optimalizácie a signifikancie sa za najlepší model ukázal ten, kde bola odstránená nevýznamná premenná Freedom to make life choices. Tento zjednodušený model dosiahol mierne vyššie \(R^2\) (\(0.7262\)) a obe zostávajúce premenné sa stali plne štatisticky významnými.

Ako alternatíva bola otestovaná Analýza Hlavných Komponentov (PCA), ktorá úspešne vyriešila multikolinearitu, ale viedla k výraznému poklesu predikčnej sily (\(R^2 = 0.562\)).

Hrebeňová regresia potvrdila prítomnosť multikolinearity tým, že s rastúcim penalizačným faktorom \(\lambda\) prerozdeľovala váhy medzi korelovanými premennými (Log GDP per capita a Social support), čím stabilizovala odhady.

Záverom, pre dosiahnutie maximálnej numerickej stability sa odporúča použiť model s centrovanými premennými, zatiaľ čo pre najlepší štatistický výkon a plnú signifikanciu sa odporúča zjednodušený model, ktorý neobsahuje premennú Freedom to make life choices.

---