multicolinearity

library(tidyverse)
library(data.table)
library(recipes)
library(broom)

다중공선성 의미

다중공선성은 설명변수 간에 높은 상관관계가 존재하는 경우를 의미한다. 보통 설명변수 간의 correlation plot을 참고해서 다중공선성 여부를 1차적으로 확인할 수 있다. 다만 다중공선성은 각 변수의 pair 관계의 상관관계가 아닌 다변수간의 상관관계를 의미하는 포괄적인 개념이므로, correlation plot과 분산팽창계수($VI{F}_k$), condition number 등을 이용해서 판단한다.

다중공선성이 있을 경우 발생하는 문제점

다중공선성이 있을 경우 회귀계수 추정량의 분산이 커지는 문제가 있다. 추정량의 분산이 크다는 의미는 데이터가 바뀌었을 때, 값의 변동이 크다는 의미이고, 이는 데이터를 통해 구한 회귀계수를 신뢰할 수 없다는 것을 의미한다.

Example

먼저 다중회귀모형에서 다중공선성이 있을 때, 회귀계수의 추정량의 분산이 어떻게 변화하는지를 알아보자.

간단하게, 설명변수 $X_1,X_2$, 절편이 존재하는 다중회귀모형의 $\beta$의 공분산은 다음과 같이 구할 수 있다.

$$\begin{array}{r}cov(\hat{\beta })=\left[\begin{array}{ccc}var(\widehat{{\beta }_0})& cov({\hat{\beta }}_0,{\hat{\beta }}_1)& cov({\hat{\beta }}_0,{\hat{\beta }}_2)\\ cov({\hat{\beta }}_1,{\hat{\beta }}_0)& var({\hat{\beta }}_1)& cov({\hat{\beta }}_1,{\hat{\beta }}_2)\\ cov({\hat{\beta }}_2,{\hat{\beta }}_0)& cov({\hat{\beta }}_2,{\hat{\beta }}_1)& var({\hat{\beta }}_2)\end{array}\right]={\sigma }^2(X^{t}X{)}^{-1}\end{array}$$

표준화된 설명변수 $X_1,X_2$에 대한 상관계수가 $0.99,0.95,0.85$로 높은 경우를 가정해보자.

xtx <- matrix(c(1, 0.99, 0.85, 0.99, 1, 0.95, 0.85, 0.95, 1), ncol = 3)
xtx

     [,1] [,2] [,3]
[1,] 1.00 0.99 0.85
[2,] 0.99 1.00 0.95
[3,] 0.85 0.95 1.00

solve(xtx)

       [,1]   [,2]    [,3]
[1,] -15.60  29.20 -14.480
[2,]  29.20 -44.40  17.360
[3,] -14.48  17.36  -3.184

$X$는 표준화된 설명변수라고 할 때, 상관계수가 $0.99,0.95,0.85$로 높을 경우 역행렬을 구해보면 역행렬의 각 element는 모두 큰 값으로 커지게 된다($\sigma =1$). 즉, 위의 식에 대입해서 보면 $\hat{\beta }$의 분산이 매우 커지는 것을 볼 수 있다.

분산팽창계수(variance inflation factor)

다중공선성을 측정하는 대표적인 측도로 분산팽창계수가 있다.

$$\begin{array}{r}VI{F}_k=\frac{1}{1-R_k^2}\end{array}$$

$R_k^2$는 $k$번째 설명변수를 반응변수로, 나머지 설명변수를 이용해서 회귀모형을 적합했을 때의 결정계수를 의미한다.

Example

set.seed(1)
X1 <- rnorm(10)
X2 <- X1*runif(10)
y <- rnorm(10)
dat <- data.frame(X1, X2, y)

cor(dat)

           X1         X2          y
X1  1.0000000  0.8207918 -0.2047086
X2  0.8207918  1.0000000 -0.5307451
y  -0.2047086 -0.5307451  1.0000000

$X_1$, $X_2$의 경우 상관계수 값이 $0.82$로 매우 높은 양의 상관관계가 존재한다고 볼 수 있다. 이 경우 회귀모형을 적합시켰을 때의 결과를 살펴보자.

fit <- lm(y~., dat)

library(car)
vif(fit)

      X1       X2 
3.064656 3.064656 

fit2 <- lm(X1 ~ .-y, dat)
result <- glance(fit2)

1/(1-result$r.squared)

[1] 3.064656

해석을 해보면 $X_k$가 다른 설명변수와 관계가 있을 경우 $R_k^2\approx 1$이 될 것이다. $R_k^2\approx 1$이 되면 $VI{F}_k$의 분모의 값은 매우 작아지기 때문에, 결과적으로 $VI{F}_k$는 커지게 된다. 예를 들어 $R_k^2=0.9$일 경우, $VI{F}_k=10$이 된다. 보통 $VI{F}_k>10$의 경우 다중공선성이 있다고 판단하지만, 변수 간의 관계, 의미를 보고 주관적으로 결정할 수 있다.

condition number

다중공선성은 $X^{t}X$의 고유값을 통해 판단해볼 수 있다. 다중공선성이 존재할 경우 $X^{t}X$는 ill-conditioned matrix이므로, 다중공선성의 원인이 되는 열에 해당하는 고유값은 $0$에 가까워진다. condition number의 정의는 다음과 같다.

$n\times p$ 행렬 $X$를 SVD하여 만든 singular value를 $d_1\ge d_2\ge ,\cdots \ge d_p\ge 0$라 할 때, condition number $K(X)$는 다음과 같다. $X^{t}X$정방행렬이므로 singular value는 고유값으로 표현할 수 있다.

$$\begin{array}{r}K(X)=\frac{d_1}{d_p}=\sqrt{\frac{{\lambda }_1}{{\lambda }_p}}\end{array}$$

condition number는 R base 함수인 kappa()를 통해 구할 수 있다.

xtx <- t(model.matrix(fit))%*%model.matrix(fit)
kappa(xtx, exact = T)

[1] 33.00623

svd_d <- svd(xtx)
svd_d$d[1]/svd_d$d[3]

[1] 33.00623

보통 $K(X)\ge 10$일 경우 다중공선성이 존재한다고 판단한다.

generalized linear model에서의 다중공선성 문제

일반화선형모형의 경우 $var(\hat{\beta })$의 asymtotic result는 다음과 같다.

$$\begin{array}{r}var(\hat{\beta })=(X^{t}WX{)}^{-1}\end{array}$$

이전의 선형모형과 유사하게, 일반화선형모형의 $(X^{t}WX{)}^{-1}$를 svd 했을 때의 결과를 보면, 작은 singular value $d_i$는 선형모형과 동일하게, 일반화선형모형의 회귀계수 추정량의 분산에 영향을 주는 것을 볼 수 있다.

$$\begin{array}{rl}(X^{t}WX{)}^{-1}& =({\tilde{X}}^t\tilde{X}{)}^{-1},\tilde{X}=W^{\frac{1}{2}}X\\ & =(VD{U}^{t}UD{V}^t{)}^{-1}\\ & =(V{D}^2{V}^t{)}^{-1}\\ & =V{D}^{-2}{V}^t\\ & =\sum _{i=1}^p\frac{v_i{v}_i^t}{d_i^2}\\ & =\sum _{i=1}^p\frac{1}{d_i^2}{I}_{p\times p},V:orthonormal\end{array}$$

Example

set.seed(1)
X1 <- rnorm(10)
X2 <- X1*runif(10)
y <- rbinom(10, 1, 0.3)
dat <- data.frame(X1, X2, y)

fit <- glm(y~., family = binomial, dat)

$(X^{t}WX{)}^{-1}$는 R base 함수인 vcov()를 통해 구할 수 있다.

vcov(fit)

            (Intercept)         X1         X2
(Intercept)   0.9009862 -0.7664066  0.6210139
X1           -0.7664066  3.5509802 -7.4328529
X2            0.6210139 -7.4328529 28.2736799

$(X^{t}WX{)}^{-1}$를 직접 구해보면 다음과 같다. $w_{ii}={\hat{\pi }}_i(1-{\hat{\pi }}_i)$

p <- fit$fitted.values
w <- p*(1-p)
X <- model.matrix(fit)

solve(t(X)%*%diag(w)%*%X)

            (Intercept)         X1         X2
(Intercept)   0.9009921 -0.7664113  0.6210083
X1           -0.7664113  3.5509949 -7.4329015
X2            0.6210083 -7.4329015 28.2739839

vcov()의 $(1,1)$ element와 $var(\widehat{{\beta }_0})$의 결과가 같은 것을 볼 수 있다.

tidy(fit)$std.error[1]^2

[1] 0.9009862

따라서 GLM에서도 다중회귀모형과 마찬가지로, 모형 적합 시에 다중공선성도 체크해줘야 한다.

generalized variance inflation factor (GVIF)

일반적으로 vif는 범주형변수 혹은 polynomial 항에 대해서는 계산하지 않는다. 모델의 기본적인 구조에 의해 가상으로 생성되는 항이기 때문이다. 범주형변수 혹은 polynomial 항에 대해서 vif를 계산하고자 한다면 GVIF를 고려해볼 수 있다.

먼저 다음과 같은 회귀식을 고려해보자.

$$Y={\beta }_0+X_1{\beta }_1+X_2{\beta }_2+ϵ$$

• $X_1$ : indicator variable or polynomial term

• $X_2$ : ${\beta }_0,X_1$을 제외한 나머지 수치형 변수

$$GVI{F}_1=\frac{det(R_{11})det(R_{22})}{det(R)}$$ - $R_{11}$ : $X_1$의 correlation matrix

• $R_{22}$ : $X_2$의 correlation matrix

• $R$ : ${\beta }_0$를 제외한 $X$ 전체의 correlation matrix

GVIF는 전체 변수의 상관행렬의 determinant에서 범주형변수의 상관행렬의 determinant가 차지하는 비율로 정의되는 것을 확인할 수 있다.

Example

먼저 iris 데이터로 GVIF를 계산결과를 확인해본다. iris 데이터에는 species 변수가 범주형 변수이므로, species 변수에 대해서 GVIF를 계산해본다.

iris %>% str()

'data.frame':   150 obs. of  5 variables:
 $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
 $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
 $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
 $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
 $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

iris_dat <- iris %>% 
    mutate(Species = as.factor(Species))

fit <- lm(Sepal.Length~., iris)
X <- model.matrix(fit)

먼저 design matrix를 구해서 X로 저장한다.

X1 <- X[, 5:6]
X2 <- X[, 2:4]
X_nonconst <- X[,-1]

(det(cor(X1))*det(cor(X2)))/det(cor(X_nonconst))

[1] 40.03918

car::vif(fit)

                  GVIF Df GVIF^(1/(2*Df))
Sepal.Width   2.227466  1        1.492470
Petal.Length 23.161648  1        4.812655
Petal.Width  21.021401  1        4.584910
Species      40.039177  2        2.515482

species에 대한 계산 결과는 car::vif()와 같은 것을 확인할 수 있다.

어떻게 처리해야 하는가?

1. 그냥 처리하지 않는다
   • 풀고자 하는 문제가 예측 task라면 갖고 있는 데이터가 대표성이 있다는 가정 하에 크게 문제가 되지 않는다.
2. 변수를 제거한다
   • 변수 간의 의미가 중복되는 경우 제거하는 것을 고려해볼 수 있다
3. 다항식 혹은 상호작용항의 경우 중심화를 통해 다중공선성을 보정해볼 수 있다
4. 주성분 분석을 통해 다중공선성이 높은 변수를 그룹화할 수 있다
5. LASSO regression을 통해 변수 선택을 고려해볼 수 있다

Example

3번 중심화의 영향을 예시로 알아보자

fit1 <- lm(Sepal.Length~. + Petal.Length*Sepal.Width, iris)
car::vif(fit1)

there are higher-order terms (interactions) in this model
consider setting type = 'predictor'; see ?vif

                               GVIF Df GVIF^(1/(2*Df))
Sepal.Width                7.338253  1        2.708921
Petal.Length             102.800499  1       10.139058
Petal.Width               22.012759  1        4.691776
Species                   41.330982  2        2.535531
Sepal.Width:Petal.Length  72.972392  1        8.542388

Sepal.Width:Petal.Length를 상호작용항으로 추가했을 때, vif는 8.54로 높은 것을 볼 수 있다.

library(recipes)

center_dat <- iris %>% 
    recipe(Sepal.Length~.) %>% 
    step_center(all_numeric_predictors()) %>% 
    step_interact(~Petal.Length:Sepal.Width) %>% 
    prep() %>% 
    juice()

fit2 <- lm(Sepal.Length~., center_dat)
car::vif(fit2)

                                GVIF Df GVIF^(1/(2*Df))
Sepal.Width                 2.396217  1        1.547972
Petal.Length               23.176426  1        4.814190
Petal.Width                22.012759  1        4.691776
Species                    41.330982  2        2.535531
Petal.Length_x_Sepal.Width  1.444820  1        1.202007

centering을 한 후에 결과를 보면 Sepal.Width:Petal.Length 항의 경우 vif=1.2로 크게 개선된 것을 확인할 수 있다.

참고자료

https://statisticaloddsandends.wordpress.com/2020/11/20/variance-of-coefficients-for-linear-models-and-generalized-linear-models/

colinearity in generalized linear models, p.32, chapter 2.5.1

http://web.vu.lt/mif/a.buteikis/wp-content/uploads/PE_Book/4-5-Multiple-collinearity.html

Citation

BibTeX citation:

@online{don2024,
  author = {Don, Don and Don, Don},
  title = {Multicolinearity},
  date = {2024-01-29},
  url = {https://dondonkim.netlify.app/posts/logistic_reg},
  langid = {en}
}

For attribution, please cite this work as:

Don, Don, and Don Don. 2024. “Multicolinearity.” January 29, 2024. https://dondonkim.netlify.app/posts/logistic_reg.