分散拡大係数（VIF）と相関係数Rの関係

一般的に，10以上で多重共線性（マルチコ）があると判断される。 この場合，回帰係数が不安定になり分析結果は信頼できない。 説明変数間の相関係数とVIFとの関係は次のグラフのとおり。

options(digits = 2)

VIF <- seq(1, 15, 0.01)

get.R <- function(vif) sqrt(1 - 1 / vif)

R <- get.R(VIF)

#png("vif_r.png")
par(cex = 1.5)

matplot(x = VIF, y = R, type = "l", 
        main = "分散拡大係数VIFと相関係数Rの関係")

abline(v = 10, col = "red")

grid()

#dev.off()

VIFが10以上では非常に高い相関（＞0.9）があることが分かる。

実験データ

N <- 100 # データサイズ

set.seed(1)

x1 <- 1:N                             # 変数1：初項1公差1の等差数列
x2 <- 1 + 2 * x1 + rnorm(N, sd = 4^2) # 変数2：変数1（x1）の線形結合＋ノイズ
x3 <- 1.1 ^ (1:N)                     # 変数3：指数系列

y <- x1 + x2 + x3 + rnorm(N)

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

library(DT)
datatable(round(d, 1))

回帰モデルのフィッティング

モデル式の「.」は全説明変数を利用する設定

fit <- lm(y ~ ., data = d)

summary(fit)
## 
## Call:
## lm(formula = y ~ ., data = d)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.888 -0.619 -0.131  0.577  2.282 
## 
## Coefficients:
##              Estimate Std. Error  t value Pr(>|t|)    
## (Intercept) -1.26e-01   2.20e-01    -0.58     0.57    
## x1           1.00e+00   1.46e-02    68.74   <2e-16 ***
## x2           1.00e+00   6.81e-03   146.83   <2e-16 ***
## x3           1.00e+00   4.67e-05 21402.06   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.97 on 96 degrees of freedom
## Multiple R-squared:     1,   Adjusted R-squared:     1 
## F-statistic: 3.18e+08 on 3 and 96 DF,  p-value: <2e-16

VIF（分散拡大係数：variance inflation factor）

library(car)

(v <- vif(fit))
## x1 x2 x3 
## 19 17  2
barplot(v, main = 'VIF')

変数1（x1）と変数2（x2）で高いVIFとなった。 どちらか一方の変数だけで十分である。

演習課題

車の燃費データを使って多重共線性の有無を分析せよ。
最後の列にある車名は説明変数から除くこと（d[, -ncol(d)]）。 なお，vif関数で出力される値はGVIF （1列目に入っている; v[, 1]）を参照せよ。 （説明変数にカテゴリ変数が含まれるとVIFは表形式で出力される）

データ

単位はオリジナル（出典参照）から馴染みのあるものに変換した。

出典：【UCI Machine Learning Repository】Auto MPG Data Set

番号 内容（単位）
1 燃費（km/L）
2 気筒数（本）
3 排気量（cc）
4 馬力（hp）
5 車体重量（kg）
6 加速性能：時速97km（60mile）に到達する時間（秒）
7 製造年（年）
8 アメ車，日本車，欧州車の3区分
9 車名（英語） 
library(DT)
d <- read.csv('https://stats.dip.jp/01_ds/data/car_mileage.csv')