Multikolinearita v lineárnej regresii

Import súboru

Naimportujeme si dataset s údajmi všetkých krajín.

# import the dataset and create a data.frame udaje
udaje_svet <- read.csv("udaje/Life-Expectancy-Data-Updated.csv",header=TRUE,sep=",",dec=".",check.names = TRUE)
head(udaje_svet)

 udaje_svet <- udaje_svet[-992,]

Vyberieme si len pozorovania pre Albánsko.

# z databázy udaje_svet si vyberieme len tie pozorovania, ktoré sa týkajú Abánska 
udaje <- subset(udaje_svet, Country == "Albania")

Zvolíme si tri premenné ktoré chceme pozorovať, nainštalujeme potrebné knižnice a vytvoríme regresný model.

# vyrovnanie priebehu očakávanej dĺžky dožitia v čase
model <- lm(Life_expectancy ~ Alcohol_consumption+Adult_mortality+Incidents_HIV,data = udaje)
library(broom)
library(knitr)
library(kableExtra)

# koeficienty regresie
tidy(model) %>%
  kable(digits = 3, caption = "Odhadnuté koeficienty regresie") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

Odhadnuté koeficienty regresie

term	estimate	std.error	statistic	p.value
(Intercept)	105.254	3.678	28.614	0.000
Alcohol_consumption	-0.018	0.202	-0.091	0.929
Adult_mortality	-0.339	0.042	-8.006	0.000
Incidents_HIV	-44.297	31.242	-1.418	0.184

# kvalita vyrovnania
glance(model) %>%
  kable(digits = 3, caption = "Ukazovatele kvality vyrovnania") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

Ukazovatele kvality vyrovnania

r.squared	adj.r.squared	sigma	statistic	p.value	df	logLik	AIC	BIC	deviance	df.residual	nobs
0.953	0.941	0.299	74.885	0	3	-0.847	11.694	15.234	0.983	11	15

NA

Pre x vyberieme tri už zvolené premenné z dát, vypočítame medzi nimi korelačnú maticu s ignorovaním chýbajúcich hodnôt a výsledné korelácie zaokrúhlime na štyri desatinné miesta.

# výber premenných
X <- udaje[, c("Alcohol_consumption", "Adult_mortality", "Incidents_HIV")]

# výpočet korelačnej matice
cor_matrix <- cor(X, use = "complete.obs")

# zaokrúhlenie
round(cor_matrix, 4)

                    Alcohol_consumption Adult_mortality Incidents_HIV
Alcohol_consumption              1.0000          0.0106        0.4130
Adult_mortality                  0.0106          1.0000       -0.8066
Incidents_HIV                    0.4130         -0.8066        1.0000

Vidíme, že najsilnejší vzťah je medzi Adult_mortality a Incidents_HIV (korelácia −0.8066), čo znamená silnú negatívnu závislosť, keď jedna premenná rastie, druhá klesá.

Vzťah medzi Alcohol_consumption a Incidents_HIV je stredne silný pozitívny (0.4130) a medzi Alcohol_consumption a Adult_mortality je takmer nulový (0.0106), teda prakticky žiadny.

S pomocou knižnice knitr si môžeme vytvoriť aj vizuálne krajšiu a prehľadnejšiu tabuľku pre korelačnú maticu.

library(knitr)
round(cor_matrix, 4) %>%
  kable(caption = "Korelačná matica")

Korelačná matica

	Alcohol_consumption	Adult_mortality	Incidents_HIV
Alcohol_consumption	1.0000	0.0106	0.4130
Adult_mortality	0.0106	1.0000	-0.8066
Incidents_HIV	0.4130	-0.8066	1.0000

Následne spravíme cor testy kedy ždy kombinujeme každú premennú s každou a pozorujeme výsledky p- value.

cor.test (udaje$Adult_mortality, udaje$Alcohol_consumption)

    Pearson's product-moment correlation

data:  udaje$Adult_mortality and udaje$Alcohol_consumption
t = 0.038131, df = 13, p-value = 0.9702
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.5044203  0.5200209
sample estimates:
      cor 
0.0105751 

cor.test (udaje$Alcohol_consumption, udaje$Incidents_HIV)

    Pearson's product-moment correlation

data:  udaje$Alcohol_consumption and udaje$Incidents_HIV
t = 1.6351, df = 13, p-value = 0.126
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1258822  0.7636995
sample estimates:
      cor 
0.4130138 

cor.test (udaje$Adult_mortality, udaje$Incidents_HIV )

    Pearson's product-moment correlation

data:  udaje$Adult_mortality and udaje$Incidents_HIV
t = -4.9197, df = 13, p-value = 0.0002801
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.9332439 -0.5015452
sample estimates:
       cor 
-0.8065793 

Medzi premennými Adult_mortality a Alcohol_consumption je p-value 0.9702, čo je veľmi vysoká hodnota, takže vzťah nie je štatisticky významný. Medzi Alcohol_consumption a Incidents_HIV je p-value 0.126, čo je tiež viac ako 0.05, preto ani tento vzťah nie je štatisticky významný. Naopak medzi Adult_mortality a Incidents_HIV je p-value 0.0002801, čo je menej ako 0.05, takže tento vzťah je štatisticky významný.

Multicollinearity detection and testing

V prípade potreby nainštalujeme balík corrplot, následne ho načítame a zobrazíme korelačnú maticu ako graf s číselnými hodnotami v hornej časti.

if (!requireNamespace("corrplot", quietly = TRUE)) {
  install.packages("corrplot")}

library(corrplot)

corrplot(cor_matrix, method = "number", type = "upper")

Vidíme výrazný vzťah medzi Adult_mortality a Incidience_HIV s hodnotou -0.81, ktorá sa zobrazuje červenou farbou.

Variance Inflation Factor

Pre premennú \(x_j\) definujeme

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

kde \(R_j^2\) je koeficient determinácie z pomocnej regresie, v ktorej je \(x_j\) vysvetľovaná ostatnými regresormi.

Ak je \(R_j^2\) blízko jednej, potom je \(VIF_j\) veľký a premenná \(x_j\) je silno lineárne vysvetliteľná ostatnými premennými.

V prípade potreby nainštalujeme balík car. # install.packages(“car”) # run once if not installed

Urobíme načítanie balíka car, následne vypočítame hodnoty VIF (Variance Inflation Factors) pre model model a zobrazíme ich, aby sme posúdili multikolinearitu medzi premennými.

library(car)
# Variance Inflation Factors
vif_values <- vif(model)

vif_values

Alcohol_consumption     Adult_mortality       Incidents_HIV 
           2.035214            4.830860            5.823734 

Interpretation:

• \(VIF_j = 1\): žiadna multikolinearita,
• \(1 < VIF_j < 5\): mierna multikolinearita
• \(VIF_j \geq 5\): silná multikolinearita (niektorí autori používajú prísnejší prah \(VIF_j \geq 10\)).

Číslo podmienenosti

Číslo podmienenosti je založené na vlastných číslach matice \(X'X\). Ak sú vlastné čísla veľmi rozdielne, matica je zle podmienená.

\[ \kappa = \sqrt{\frac{\lambda_{\max}}{\lambda_{\min}}}. \]

Urobíme štandardizáciu premenných v matici X, následne vypočítame vlastné čísla korelačnej matice a z nich určíme podmienené číslo (condition number), ktoré nám slúži na posúdenie multikolinearity a stability modelu, a potom zobrazíme vlastné čísla aj výsledné podmienené číslo.

X_scaled <- scale(X)
eigen_values <- eigen(cor(X_scaled))$values

condition_number <- sqrt(max(eigen_values) / min(eigen_values))
eigen_values
condition_number

Interpretácia čísla podmienenosti: - \(\kappa \approx 1\): žiadna multikolinearita, - \(\kappa > 10\): mierna multikolinearita, - \(\kappa > 30\): silná multikolinearita.

Malá p-hodnota znamená, že zamietame hypotézu \(R=I_k\), teda medzi vysvetľujúcimi premennými existuje štatisticky významná korelačná štruktúra.

Eliminácia problému multikolinearity

Multikolinearita nie je porušením exogenity. Nie je teda automaticky dôvodom na zamietnutie OLS modelu. Problém je hlavne inferenčný: veľké štandardné chyby a nestabilné individuálne koeficienty.

Možné riešenia:

1. získať viac dát,
2. odstrániť alebo spojiť silno korelované premenné,
3. transformovať premenné,
4. použiť teoretické obmedzenia,
5. použiť regularizačné metódy, napríklad ridge regresiu.

Odstránenie jednej z korelovaných premenných

Ak sú dve premenné takmer rovnaké a ekonomická teória nevyžaduje obe, môžeme jednu z nich vynechať. Odstránime premennú Adult_mortality a pozorujeme výsledky.

model_reduced <- lm(Life_expectancy ~ Alcohol_consumption + Incidents_HIV, data = udaje)

summary(model)

Call:
lm(formula = Life_expectancy ~ Alcohol_consumption + Adult_mortality + 
    Incidents_HIV, data = udaje)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.30738 -0.07191  0.00501  0.06394  0.33087 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)         80.406412   0.629177 127.796  < 2e-16 ***
Alcohol_consumption -0.009374   0.038803  -0.242    0.813    
Adult_mortality     -0.056383   0.002029 -27.784 2.92e-12 ***
Incidents_HIV        1.342110   0.899627   1.492    0.162    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1938 on 12 degrees of freedom
Multiple R-squared:  0.9876,    Adjusted R-squared:  0.9845 
F-statistic: 318.6 on 3 and 12 DF,  p-value: 1.06e-11

summary(model_reduced)

Call:
lm(formula = Life_expectancy ~ Alcohol_consumption + Incidents_HIV, 
    data = udaje)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4850 -1.0324 -0.5959  0.9277  2.5740 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         65.54601    2.57299  25.475 1.76e-12 ***
Alcohol_consumption -0.08151    0.30065  -0.271    0.791    
Incidents_HIV       10.91938    6.45281   1.692    0.114    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.505 on 13 degrees of freedom
Multiple R-squared:   0.19, Adjusted R-squared:  0.06535 
F-statistic: 1.524 on 2 and 13 DF,  p-value: 0.2542

V pôvodnom modeli, kde boli zahrnuté všetky tri premenné (Alcohol_consumption, Adult_mortality a Incidents_HIV), model vysvetľoval veľmi vysokú časť variability v očakávanej dĺžke života, keďže R² = 0.9876 a model bol celkovo štatisticky významný (p-value ≈ 1.06e-11). Z jednotlivých premenných sa však ako štatisticky významná ukázala iba Adult_mortality (p < 0.001), zatiaľ čo Alcohol_consumption aj Incidents_HIV neboli významné.

Po odstránení premennej Adult_mortality sa model výrazne zhoršil, keďže R² kleslo na 0.19 a upravené R² dokonca na 0.065, čo znamená, že model už vysvetľuje len malú časť variability. Celý model už nie je štatisticky významný (p-value = 0.2542) a ani jednotlivé premenné (Alcohol_consumption a Incidents_HIV) nie sú štatisticky významné.

Adult_mortality je teda kľúčová premenná v modeli, bez nej model stráca vysvetľovaciu silu a výrazne sa zhoršuje jeho kvalita.