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 <- 300

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

u <- rnorm(n)

y <- 3 + 4*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(5, 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  1.335728e-04  2.559411e-05 -5.521933e-05 -1.182672e-05
## x1         2.559411e-05  1.388586e-01 -1.437087e-01 -7.248936e-04
## x2        -5.521933e-05 -1.437087e-01  1.528524e-01  9.357597e-04
## x3        -1.182672e-05 -7.248936e-04  9.357597e-04  3.143605e-03

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.56535 -0.71289  0.00933  0.56108  3.11957 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.02441    0.05874  51.486  < 2e-16 ***
## x1           3.27416    0.37880   8.644 3.48e-16 ***
## x2          -0.70724    0.39743  -1.780   0.0762 .  
## x3           0.84558    0.05700  14.836  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.017 on 296 degrees of freedom
## Multiple R-squared:  0.8671, Adjusted R-squared:  0.8658 
## F-statistic: 643.9 on 3 and 296 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, 3)
##        x1     x2     x3
## x1  1.000  0.986 -0.045
## x2  0.986  1.000 -0.052
## x3 -0.045 -0.052  1.000

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 = 103.71, df = 298, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9829918 0.9891746
## sample estimates:
##       cor 
## 0.9864286
cor.test(data_main$x1, data_main$x3)
## 
##  Pearson's product-moment correlation
## 
## data:  data_main$x1 and data_main$x3
## t = -0.78061, df = 298, p-value = 0.4357
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.15760808  0.06841734
## sample estimates:
##         cor 
## -0.04517348
cor.test(data_main$x2, data_main$x3)
## 
##  Pearson's product-moment correlation
## 
## data:  data_main$x2 and data_main$x3
## t = -0.89129, df = 298, p-value = 0.3735
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.16384653  0.06204071
## sample estimates:
##         cor 
## -0.05156238
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.

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.99114930 0.99529999 0.01355071
condition_number
## [1] 12.12191

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 + 4)/5) * 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.56535 -0.71289  0.00933  0.56108  3.11957 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.02441    0.05874  51.486  < 2e-16 ***
## x1           3.27416    0.37880   8.644 3.48e-16 ***
## x2          -0.70724    0.39743  -1.780   0.0762 .  
## x3           0.84558    0.05700  14.836  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.017 on 296 degrees of freedom
## Multiple R-squared:  0.8671, Adjusted R-squared:  0.8658 
## F-statistic: 643.9 on 3 and 296 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.65975 -0.71767 -0.02682  0.61284  3.06464 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.02313    0.05895   51.28   <2e-16 ***
## x1           2.60923    0.06245   41.78   <2e-16 ***
## x3           0.84990    0.05715   14.87   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.02 on 297 degrees of freedom
## Multiple R-squared:  0.8657, Adjusted R-squared:  0.8648 
## F-statistic: 957.3 on 2 and 297 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.2424027
ridge_cv$lambda.1se
## [1] 0.5102345
coef(ridge_cv, s = "lambda.min")
## 4 x 1 sparse Matrix of class "dgCMatrix"
##             lambda.min
## (Intercept)   3.027200
## x1            1.507383
## x2            1.027079
## x3            0.782122
plot(ridge_cv)

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.