渋谷駅前で働くデータサイエンティストのブログ

元祖「六本木で働くデータサイエンティスト」です / 道玄坂→銀座→東京→六本木→渋谷駅前

VIF (Variance Inflation Factor)を計算する関数を書いてみた(備忘録)

多重共線性(multicolinearity)の代表的指標として頻繁に用いられるVIF (Variance Inflation Factor)というと、Rでは普通に{car}とか{usdm}とかのパッケージに実装された関数があるのでそれらを利用すれば良いのですが、ちょっと訳あって自分で実装してみることにしました。ということで、備忘録的に書き残しておきます。

定義


普通にWikipedia記事を読めばおしまいです。


 \displaystyle VIF_j = \frac{1}{1 - R_j^2} (ただし j番目の説明変数に対するVIFとする: R_j^2 j番目の説明変数に対して残りの説明変数全てで線形回帰モデルを推定した際に得られる自由度調整前の決定係数)

ちなみに他のどの資料を読んでも全く同じことが書いてあるので、これで十分だと思われます。


実装


定義の字面だけ見ると瞬殺に見えますが、Rだとlm関数が愚かにもformula式しかモデル定義を受け付けない(低レベルメソッドのlm.fitとかだと目的変数と説明変数を別々にmatrixで与えることができる)ので、手間はかかりますが入力したdata.frameから無理やりformula式を自動的に生成するようにしてやる必要があります。


そんなわけで、完成したクソコードがこちらです。

vif <- function(d, idx = NULL){
  d1 <- d
  if (!is.null(idx)){
    d1 <- d[, -idx]
  }
  dnam <- names(d1)
  m <- ncol(d1)
  
  fm <- c()
  for (i in 1:m){
    fm[[i]] <- as.formula(paste(paste(dnam[i], '~'), paste(dnam[-i], collapse = '+')))
  }
  
  vif_val <- data.frame(var = dnam, vif = 0)
  for (i in 1:m){
    vif_val$vif[i] <- 1 / ( 1- summary(lm(fm[[i]], d1))$r.square)
  }
  return(vif_val)
}

意味もなく目的変数の列をomitするための引数が入っていますが、普通にd[, -idx]とかやって渡せばいいのでぶっちゃけ完全にどうでも良いです。


動作確認


同じようにdata.frameを入力して動くusdm::vifと挙動を比較してみます。使ったデータは、適当に前処理して手元のローカル環境に置きっぱなしにしてあるBoston Housingです。

d <- read.csv('housing.csv')

vif(d, 14)
#R>        var      vif
#R> 1     CRIM 1.792192
#R> 2       ZN 2.298758
#R> 3    INDUS 3.991596
#R> 4     CHAS 1.073995
#R> 5      NOX 4.393720
#R> 6       RM 1.933744
#R> 7      AGE 3.100826
#R> 8      DIS 3.955945
#R> 9      RAD 7.484496
#R> 10     TAX 9.008554
#R> 11 PTRATIO 1.799084
#R> 12       B 1.348521
#R> 13   LSTAT 2.941491

usdm::vif(d[, -14])
#R>    Variables      VIF
#R> 1       CRIM 1.792192
#R> 2         ZN 2.298758
#R> 3      INDUS 3.991596
#R> 4       CHAS 1.073995
#R> 5        NOX 4.393720
#R> 6         RM 1.933744
#R> 7        AGE 3.100826
#R> 8        DIS 3.955945
#R> 9        RAD 7.484496
#R> 10       TAX 9.008554
#R> 11   PTRATIO 1.799084
#R> 12         B 1.348521
#R> 13     LSTAT 2.941491

完全に同じ結果を返すことが確認できました。


追記


usdm::vifのコードを見に行ったら同じようなやり方だったので、安堵しました(笑)。