多重共線性(マルチコ)について / 重回帰分析

　多重共線性(multicollinearity)とは、説明変数同士の相関が強いことをいう。VIFが10以上となると多重共線性があるとされる。これに従うと、多重共線性の正確な定義が可能だ。説明変数X1について多重共線性が発生しているとは、以下の「削ぎ落としの回帰」での決定係数が0.9以上になるとき、多重共線性が発生していると定義できる。

$$重回帰モデル　Y=\beta_0 +\beta_1 X_1 +\beta_2 X_2 +\cdots + \beta_kX_k+U$$

$$削ぎ落としの回帰　x_{1i}=\widehat{\delta_0} +\widehat{\delta_2} x_{2i} +\cdots + \widehat{\delta_k}x_{ki}+\widehat{u_{1i}}$$

$$VIF_1=\frac{1}{1-R^2_1}　R^2_1は削ぎ落としの回帰の決定係数$$

　多重共線性の本質的な問題は、回帰係数の標準誤差が大きなり、推定結果が安定しなくなることだ。多重共線性の実務的な問題は、統計的に有意な結果が現れにくくなることだ。シミュレーションしてみよう。次の母集団に従う標本を乱数発生させて最小二乗推定値を計算するのを、5000回繰り返して、最小二乗推定値の分布を求めた。その際は、多重共線性なし（VIF=2）、多重共線性あり（VIF=10)、深刻な多重共線性あり（VIF=50）で比較する。その結果が、図1である。図１からは、多重共線性が標準誤差を大きくしている一方で、不偏性は失っていないことがわかる。

$$母集団モデル：Y=2 + X_1- 2X_2 +U$$

　多重共線性は、統計的に有意な結果が出ているのなら、気にしなくてよい。統計的に有意な結果が出ているのならば、十分に標準誤差が小さいわけであり、上記の問題はそもそも発生していない。なお、「最小二乗推定量の不偏性・一致性・漸近正規性について」にある通り、多重共線性があっても、最小二乗推定量に不偏性や一致性は失われない。

【追記】

　削ぎ落としの回帰の決定係数が1となることを、完全な共線関係（perfect collinearity）と呼ぶ。これは多重共線性(multicollinearity)と異なり、最小二乗法による推定に致命的な影響を与える。推定量の分散が無限大に発散して、推定不可能となるのだ。以下の式が標準誤差である。（詳細は「重回帰分析の標準誤差について」）

$$SE \left( \widehat{\beta_1} \right)=\sqrt{ \frac{\widehat{Var(U)}}{ var(x_1) (1-R_1^2)n }} $$

　ダミー変数をすべて入れると、完全な共線関係が起きる。例えば、男ダミーと女ダミーを入れたとしよう。片方が決まれば、片方が決まるので、決定係数1となる。

　描画とシミュレーションのRコードは以下の通り。

#パッケージの呼び出し
library(ggplot2)
　
#準備
sample_size0 <- 300
trial_number <- 5000
vif1 <- rep(NA, trial_number)　
vif2 <- rep(NA, trial_number)　
vif3 <- rep(NA, trial_number)
beta_1_x2 <- rep(NA, trial_number)　
beta_1_x3 <- rep(NA, trial_number)　
beta_1_x4 <- rep(NA, trial_number)
　
#vif調整
sample_size0 <- 300
sd1 <- 0.63
sd2 <- 1.5
sd3 <- 5.2
for(j in 1:trial_number){ 
  x2 <- x1+rnorm(sample_size0, mean =0, sd=sd1)
  x3 <- x1+rnorm(sample_size0, mean =0, sd=sd2)
  x4 <- x1+rnorm(sample_size0, mean =0, sd=sd3)
  vif1[j] <- 1/(1-cor(x1,x2)) #x1とx2のVIF
  vif2[j] <- 1/(1-cor(x1,x3)) #x1とx3のVIF
  vif3[j] <- 1/(1-cor(x1,x4)) #x1とx4のVIF
}
mean(vif1)
mean(vif2)
mean(vif3)
　
#VIF=50
for(j in 1:trial_number){ 
  x1<- rnorm(sample_size0, 15,3)
  x2 <- x1+rnorm(sample_size0, mean =0, sd=sd1)
  u <- rnorm(sample_size0, mean =0, sd=3)
  y <- 2 + 1*x1- 2*x2 +u
  answer <- lm(y~x1+x2)
  beta_1_x2[j] <- summary(answer)$coef["x1","Estimate"]
}
　
#VIF=10
for(j in 1:trial_number){ 
  x1<- rnorm(sample_size0, 15,3)
  x3 <- x1+rnorm(sample_size0, mean =0, sd=sd2)
  u <- rnorm(sample_size0, mean =0, sd=3)
  y <- 2 + 1*x1- 2*x3 +u
  answer <- lm(y~x1+x3)
  beta_1_x3[j] <- summary(answer)$coef["x1","Estimate"]
}
　
#VIF=2
for(j in 1:trial_number){ 
  x1<- rnorm(sample_size0, 15,3)
  x4 <- x1+rnorm(sample_size0, mean =0, sd=sd3)
  u <- rnorm(sample_size0, mean =0, sd=3)
  y <- 2 + 1*x1- 2*x4 +u
  answer <- lm(y~x1+x4)
  beta_1_x4[j] <- summary(answer)$coef["x1","Estimate"]
}
　
#推定量の分布の描画
library(tidyverse)
tibble(VIF50=beta_1_x2, VIF10=beta_1_x3, VIF02=beta_1_x4)%>% 
  gather(key=n,value=value,VIF50,VIF10,VIF02)%>% 
  ggplot() +
  geom_histogram(aes(x = value, y =..density..),bins=100,fill="skyblue", color= "black") +
  stat_function(fun = dnorm, args=list(mean=mean(beta_1_x2), sd=sqrt(var(beta_1_x2))))+
  stat_function(fun = dnorm, args=list(mean=mean(beta_1_x3), sd=sqrt(var(beta_1_x3))))+
  stat_function(fun = dnorm, args=list(mean=mean(beta_1_x4), sd=sqrt(var(beta_1_x4))))+
  geom_vline(xintercept = 1, color = "red") +
  facet_grid(n~.)+
  xlab("最小二乗推定量の分布　黒の曲線は正規分布") + 
  ylab ("確率密度") +
  theme_grey(base_family = "HiraKakuPro-W3")