多重共線性を出してみよう

観測数100ぐらいでも、誤差項の分散など次第では多重共線性が推定結果に影響を与える例を作ってみました。

多重共線性が無いケース

set.seed(2103)
n <- 100
x1 <- 0.3 * runif(n, min=0, max=10) + rnorm(n, sd=1)
x2 <- 0.7 * runif(n, min=0, max=10) + rnorm(n, sd=1)
cor(x1, x2) # 説明変数間は低い相関係数
e <- rnorm(n, sd=10)
y <- 1 + 2*x1 + 3*x2 + e
summary(lm(y ~ x1 + x2))

真のモデルにx2が無い場合に、x1とx2で推定してもそれらしい推定結果になります。

y <- 1 + 2*x1 + e
summary(lm(y ~ x1 + x2))

分散拡大係数（VIF）は1.02弱と小さいです。

1/(1-summary(lm(x1 ~ x2))$r.squared)

多重共線性があるケース

係数の大きさが変化して、有意性が無くなるケース。誤差項の分散は無いケースと同じ程度です。

set.seed(2103)
n <- 100
x <- runif(n, min=0, max=10)
x1 <- 0.3 * x + rnorm(n, sd=0.5)
x2 <- 0.7 * x + rnorm(n, sd=0.5)
cor(x1, x2) # 説明変数間の相関係数が高い
e <- rnorm(n, sd=10)

真のモデルにx2が無い場合に、x1とx2で推定するとx1の有意性も消えます。

y <- 1 + 2*x1 + e
summary(lm(y ~ x1))
summary(lm(y ~ x1 + x2))

VIFは2.81強と大きくは無いですが、観測数が小さいため影響しています。

1/(1-summary(lm(x1 ~ x2))$r.squared)

多重共線性が深刻なケース

さすがにレアなのですが、真のモデルが $y = 1 + 2 x_1 + \epsilon$ のとき、 $x_1$ と相関する変数 $x_2$ を加えて $y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \epsilon$ を推定すると、 $x_1$ の有意性が消えて係数が逆転し、 $x_2$ の係数が有意になるデータをつくります。誤ったモデルの方が、自由度調整済み重相関係数は改善することに注意してください。

set.seed(737)
n <- 100
x <- runif(n, min=0, max=10)
x1 <- 0.3 * x + rnorm(n, sd=0.5)
x2 <- 0.7 * x + rnorm(n, sd=0.5)
e <- rnorm(n, sd=10)
# 説明変数の高い相関，高いVIF
sprintf("cor:%f, VIF:%f", cor(x1, x2), 1/(1-summary(lm(x1 ~ x2))$r.squared))
# 真のモデル/データ生成モデル
y <- 1 + 2*x1 + e
# 真のモデルに近い係数が有意に得られる
summary(lm(y ~ x1))
# x1の係数が逆転し有意性が消え、x2の係数が有意になる
summary(lm(y ~ x1 + x2))

サンプルサイズが小さければ小さいほど、誤差項の分散が大きければ大きいほど、VIFが高ければ高いほど色々とおきます。例ではR²は低いですが、R²が0.6ぐらいあっても、VIFが極端に高いと多重共線性が問題になる可能性はあります*1。

LOOCV

多重共線性が深刻なケースでLOOCVをかけるとx2だけで回帰した方がよくなります。

library(boot)
df01 <- data.frame(y, x1, x2)
frml <- c("y ~ x1", "y ~ x2", "y ~ x1 + x2")
delta <- matrix(NA, length(frml), 2)
rownames(delta) <- frml
colnames(delta) <- c("raw", "adjusted")
for(i in 1:length(frml)){
delta[i, ] <- cv.glm(df01, glm(formula(frml[i]), data=df01))$delta
}
print(delta)