Multikolinearita v lineárnej regresii

1. Multikolinearita – definícia

Uvažujme lineárny regresný model

\[ y_i = \beta_0 + \beta_1 x_{1i} + \cdots + \beta_k x_{ki} + u_i, \]

alebo maticovo

\[ y = X\beta + u. \]

Multikolinearita znamená, že niektoré vysvetľujúce premenné v matici \(X\) sú silno lineárne závislé. Dokonalá multikolinearita nastáva vtedy, keď existuje nenulový vektor \(a\), pre ktorý

\[ Xa = 0. \]

Vtedy matica \(X'X\) nie je regulárna a OLS odhad

\[ \hat\beta = (X'X)^{-1}X'y \]

neexistuje.

V praxi častejšie pozorujeme nedokonalú, ale silnú multikolinearitu. Potom \(X'X\) je síce invertibilná, ale zle podmienená. Odhady existujú, no môžu byť veľmi nepresné.

Simulované dáta

Vytvoríme hlavný dátový súbor data_main, ktorý budeme používať v ďalších častiach.

Premenné:

• y je závislá premenná,
• x1, x2, x3 sú vysvetľujúce premenné,
• x2 je konštruovaná tak, aby bola silno korelovaná s x1,
• x3 je približne nezávislá od x1 a x2.

n <- 200

x1 <- rnorm(n)
x2 <- 0.95 * x1 + rnorm(n, sd = 0.15)
x3 <- rnorm(n)

u <- rnorm(n)

y <- 1 + 2*x1 - 1.5*x2 + 0.8*x3 + u

data_main <- data.frame(y = y, x1 = x1, x2 = x2, x3 = x3)

head(data_main)

2.Dôsledky multikolinearity pre OLS odhad

Pri multikolinearite OLS odhad zostáva nestranný za štandardných exogenitných predpokladov,

\[ E(u)=0, \]

ale jeho rozptyl môže byť veľký:

\[ Var(\hat\beta)=\sigma^2(X'X)^{-1}. \]

Ak sú stĺpce matice \(X\) takmer lineárne závislé, niektoré prvky matice \((X'X)^{-1}\) sú veľké.

X <- cbind(1, x1, x2, x3)
colnames(X) <- c("intercept", "x1", "x2", "x3")

XtX <- t(X) %*% X

XtX_inv <- solve(XtX)

XtX_inv

              intercept           x1           x2            x3
intercept  0.0050154322  0.001439290 -0.001476833 -0.0001822219
x1         0.0014392904  0.207599270 -0.213616763 -0.0015627618
x2        -0.0014768332 -0.213616763  0.225963071  0.0018298489
x3        -0.0001822219 -0.001562762  0.001829849  0.0054181885

To vedie k veľkým štandardným chybám, nízkym t-štatistikám a nestabilným odhadom \(\beta\) koeficientov. Pripomíname

\[var(\hat\beta_i) = \sigma^2 (X^TX)^{-1}_{ii}\]

model_ols <- lm(y ~ x1 + x2 + x3, data = data_main)
summary(model_ols)

Call:
lm(formula = y ~ x1 + x2 + x3, data = data_main)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.53250 -0.67022  0.07388  0.70545  2.68627 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.97304    0.07277  13.372  < 2e-16 ***
x1           1.42700    0.46815   3.048  0.00262 ** 
x2          -0.94615    0.48842  -1.937  0.05416 .  
x3           0.83698    0.07563  11.067  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.027 on 196 degrees of freedom
Multiple R-squared:  0.4671,    Adjusted R-squared:  0.459 
F-statistic: 57.27 on 3 and 196 DF,  p-value: < 2.2e-16

3. Multicollinearity detection and testing

3.1 Korelačná analýza vysvetľujúcich premenných

Najjednoduchší diagnostický krok je korelačná matica medzi vysvetľujúcimi premennými.

X <- data_main[, c("x1", "x2", "x3")]

cor_matrix <- cor(X)
round(cor_matrix, 4)

        x1      x2      x3
x1  1.0000  0.9863 -0.0302
x2  0.9863  1.0000 -0.0384
x3 -0.0302 -0.0384  1.0000

Vysoká absolútna hodnota korelácie, napríklad \(|r_{ij}| > 0.8\), môže indikovať problém multikolinearity. Sama o sebe však nestačí, pretože multikolinearita môže byť aj viacrozmerná: premenná môže byť lineárnou kombináciou viacerých ostatných premenných.

cor.test(data_main$x1, data_main$x2)

    Pearson's product-moment correlation

data:  data_main$x1 and data_main$x2
t = 84.039, df = 198, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9818871 0.9895986
sample estimates:
      cor 
0.9862705 

cor.test(data_main$x1, data_main$x3)

    Pearson's product-moment correlation

data:  data_main$x1 and data_main$x3
t = -0.4256, df = 198, p-value = 0.6709
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1682678  0.1089655
sample estimates:
        cor 
-0.03023254 

cor.test(data_main$x2, data_main$x3)

    Pearson's product-moment correlation

data:  data_main$x2 and data_main$x3
t = -0.54143, df = 198, p-value = 0.5888
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1762504  0.1008293
sample estimates:
       cor 
-0.0384496 

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

library(corrplot)

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

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

install.packages(“car”) # run once if not installed

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

vif_values

      x1       x2       x3 
36.74934 36.77011  1.00366 

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\)).

3.3 Čí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}}}. \]

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

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

[1] 1.98865618 0.99764855 0.01369527

condition_number

[1] 12.05021

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

3.4 Farrar–Glauberov test

Farrar–Glauberov test skúma multikolinearitu pomocou korelačnej matice regresorov. Nech \(R\) je korelačná matica vysvetľujúcich premenných, \(n\) je počet pozorovaní a \(k\) počet regresorov.

Celkový \(\chi^2\) test

Testujeme nulovú hypotézu, že regresory nie sú navzájom korelované:

\[ H_0: R = I_k. \]

Testovacia štatistika je

\[ \chi^2 = -\left(n-1-\frac{2k+5}{6}\right)\ln |R|. \]

Stupne voľnosti sú

\[ df = \frac{k(k-1)}{2}. \]

fg_chi_square <- function(X) {
  X <- as.data.frame(X)
  R <- cor(X)
  n <- nrow(X)
  k <- ncol(X)

  chi_stat <- -(n - 1 - (2*k + 5)/6) * log(det(R))
  df <- k * (k - 1) / 2
  p_value <- 1 - pchisq(chi_stat, df = df)

  data.frame(
    chi_square = chi_stat,
    df = df,
    p_value = p_value,
    det_R = det(R)
  )
}

fg_chi <- fg_chi_square(X)
round(fg_chi, 4)

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.

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

4.1 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ť.

model_reduced <- lm(y ~ x1 + x3, data = data_main)

summary(model_ols)

Call:
lm(formula = y ~ x1 + x2 + x3, data = data_main)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.53250 -0.67022  0.07388  0.70545  2.68627 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.97304    0.07277  13.372  < 2e-16 ***
x1           1.42700    0.46815   3.048  0.00262 ** 
x2          -0.94615    0.48842  -1.937  0.05416 .  
x3           0.83698    0.07563  11.067  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.027 on 196 degrees of freedom
Multiple R-squared:  0.4671,    Adjusted R-squared:  0.459 
F-statistic: 57.27 on 3 and 196 DF,  p-value: < 2.2e-16

summary(model_reduced)

Call:
lm(formula = y ~ x1 + x3, data = data_main)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.75173 -0.64762  0.09045  0.69569  2.69396 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.96686    0.07320  13.208  < 2e-16 ***
x1           0.53254    0.07780   6.845 9.46e-11 ***
x3           0.84464    0.07605  11.106  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.035 on 197 degrees of freedom
Multiple R-squared:  0.4569,    Adjusted R-squared:  0.4514 
F-statistic: 82.87 on 2 and 197 DF,  p-value: < 2.2e-16

Pozor: vynechanie relevantnej premennej môže viesť k vychýleniu odhadov. Preto sa premenné nemajú vyhadzovať mechanicky iba podľa korelácie.

4.2 Ridge regression

Ridge regresia rieši problém nestability tým, že minimalizuje penalizovanú sumu štvorcov:

\[ \hat\beta^{ridge} = \arg\min_{\beta} \left[ \sum_{i=1}^n (y_i - x_i'\beta)^2 + \lambda \sum_{j=1}^k \beta_j^2 \right]. \]

Parameter \(\lambda \geq 0\) určuje silu penalizácie. Ak \(\lambda=0\), dostaneme OLS. Ak \(\lambda\) rastie, koeficienty sa zmenšujú smerom k nule.

Maticovo, pri centrovaných a škálovaných regresoroch,

\[ \hat\beta^{ridge} = (X'X+\lambda I)^{-1}X'y. \]

Ridge znižuje rozptyl odhadov za cenu určitého vychýlenia. Preto je užitočná najmä pri predikcii alebo pri veľmi nestabilných OLS koeficientoch.

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

library(glmnet)

X_mat <- model.matrix(y ~ x1 + x2 + x3, data = data_main)[, -1]
y_vec <- data_main$y

ridge_cv <- cv.glmnet(
  x = X_mat,
  y = y_vec,
  alpha = 0,
  standardize = TRUE
)

ridge_cv$lambda.min

[1] 0.07977274

ridge_cv$lambda.1se

[1] 0.8164987

coef(ridge_cv, s = "lambda.min")

4 x 1 sparse Matrix of class "dgCMatrix"
            lambda.min
(Intercept) 0.96790458
x1          0.47731517
x2          0.02760042
x3          0.79830075

plot(ridge_cv)

Porovnajme OLS a ridge koeficienty.

ols_coef <- coef(model_ols)

ridge_coef <- as.matrix(coef(ridge_cv, s = "lambda.min"))

comparison <- data.frame(
  coefficient = rownames(ridge_coef),
  OLS = as.numeric(ols_coef[rownames(ridge_coef)]),
  Ridge = as.numeric(ridge_coef[, 1])
)

round_numeric(comparison, 4)

Ridge regresia typicky nemení cieľ z inferencie na jednotlivé nepenalizované parametre, ale pomáha stabilizovať odhady a zlepšiť predikčnú výkonnosť v prítomnosti silnej multikolinearity.

4.3 Predikčné porovnanie OLS a ridge

set.seed(321)

train_id <- sample(1:nrow(data_main), size = 0.7*nrow(data_main))

train <- data_main[train_id, ]
test <- data_main[-train_id, ]

ols_train <- lm(y ~ x1 + x2 + x3, data = train)

X_train <- model.matrix(y ~ x1 + x2 + x3, data = train)[, -1]
y_train <- train$y

X_test <- model.matrix(y ~ x1 + x2 + x3, data = test)[, -1]
y_test <- test$y

ridge_train_cv <- cv.glmnet(
  x = X_train,
  y = y_train,
  alpha = 0,
  standardize = TRUE
)

pred_ols <- predict(ols_train, newdata = test)
pred_ridge <- predict(ridge_train_cv, newx = X_test, s = "lambda.min")

mse_ols <- mean((y_test - pred_ols)^2)
mse_ridge <- mean((y_test - pred_ridge)^2)

data.frame(
  model = c("OLS", "Ridge"),
  test_MSE = c(mse_ols, as.numeric(mse_ridge))
)

Zhrnutie

Multikolinearita znamená silnú lineárnu závislosť medzi vysvetľujúcimi premennými. Pri dokonalej multikolinearite OLS odhad neexistuje. Pri silnej, ale nedokonalej multikolinearite OLS existuje, no štandardné chyby môžu byť veľké a odhady nestabilné.

V praxi je vhodné použiť viacero diagnostík: korelačnú maticu, VIF, číslo podmienenosti a Farrar–Glauberov test. Pri riešení problému treba zohľadniť ekonomickú teóriu. Ridge regresia je vhodná najmä vtedy, keď je cieľom stabilná predikcia.