と。

Github: https://github.com/8-u8

多重共線性ってどれくらいのあたりから問題になりそうなのか？という話

01_統計モデリング 00_データ分析(全般) 07_R言語

今月まだ1個も書いてないので書きます。

こんなツイートをみました。

回帰分析ではよく耳にする問題「多重共線性」です。
「説明変数同士の相関が強いと、回帰係数の推定が不安定になってしまう」という問題ですね*1。
とはいえ、「ぶっちゃけどれくらい相関が強いとヤバいんだ？」みたいなものはあまりしられていません。
VIFなどの指標が10くらい、というのはありつつ、そもそも「多重共線性が起きている」としていい場合というのはあまり議論されることが多くないように思われます。
今回は説明変数同士の相関がどの程度あると多重共線性が起きたっぽい、というところの判別がつくのか、考えてみます。

ソース

後述のソースは以下でも参照できます。
ライブラリを一切使っていないので、誰でも走らせられると思われます。
github.com

多重共線性とは

そもそも明確な定義付けがなされている本が、きぬいとの手元では下記しかありませんでした……

計量経済学 New Liberal Arts Selection

計量経済学 New Liberal Arts Selection

ここでは多重共線性「ではない定義」を

「説明変数の間に線形関係がない」

という形で書いています。つまるところ説明変数で構成される行列が一次独立であれば、
最も理想的な形で、多重共線性「ではない」と言えます。
逆に、一次従属の場合、逆行列を求めることができなくなってしまうので、回帰分析のパラメータ推定法によっては、
推定値を得られないことがあります。

また、以下の本には「行列式」の講において、多重共線性の仕組みについて記述されています*2。

統計学のための数学入門30講 (科学のことばとしての数学)

統計学のための数学入門30講 (科学のことばとしての数学)

  • 作者:靖, 永田
  • 発売日: 2005/04/01
  • メディア: 単行本

そういえばこの辺のシンプルな事例は大昔にgitに上げていました。
多重共線性の話は何回かやってます＞＜；

github.com

一方、上記によれば完全な一次従属ではない場合でも、説明変数同士の相関が強ければ、
推定値が不安定になるとのこと。実務経験上もよくあります。
具体的には、推定値の標準誤差が大きくなり、統計的有意性を保証できなくなることがあります。
ただ、「どのくらい強い相関だとヤバいのか」はあまり知られていないように思います。
ちょっと知りたかったので適当にコードを書いてシミュレートしてみました。

## init
seed(1234)
x1 <- rnorm(100,0,1)

## multico_generator
multico_regressor <- function(x2_std = 100,
                              x1 = x1,
                              beta1 = 0.4, # x1の真のパラメータ
                              beta2 = 0.3, # x2の真のパラメータ
                              intercept = 1.0 # 切片の真のパラメータ
                              ){
  #### 多重共線性を起こすための説明変数 ####
  x2 <- x1 + rnorm(100, 0, x2_std)
  #### 真のパラメータでyを構成する ####
  y  <- beta1*x1 + beta2 * x2 + intercept + rnorm(100, 0, 0.5)
  print(paste0(cat("correlation of x1 and x2: \n"),round(cor(x1, x2), digits = 3)))
  cat(c("------------------------------------ \n"))
  usedata <- data.frame(y = y, 
                      x1 = x1,
                      x2 = x2)

  model <-  lm(y~., data = usedata)
  #### 出力 ####
  cat("true parameters are \n")
  print(paste0(c(intercept, beta1, beta2)))
  cat(c("------------------------------------\n", "coefficients are \n"))
  print(round(summary(model)$coefficients[,1], digits = 3))
  cat(c("------------------------------------\n", "standard error are \n"))
  print(round(summary(model)$coefficients[,2], digits = 3))
  cat(c("------------------------------------ \n"))
}

上記のコードでは、x2という変数をx1に正規分布に従う乱数を噛ませる形で表現しています。 正規分布に従う乱数なので、標準偏差を小さくとれば実現値の幅は狭いですし、
逆に大きくとれば実現値の幅が大きくなります。
これによって「x1とx2の相関関係を制御」しています。

今回は真の値として

  • x1の回帰係数: 0.4
  • x2の回帰係数: 0.3
  • 切片: 1.0

と設定しています。

相関が弱い場合

x2を作るときの標準偏差が100の場合、つまりx1とx2の相関が弱い場合の結果が以下です*3。
相関は0.117で、決して強いとは言えない数字です。

-*-*-*-*-* x2_std = 100 *-*-*-*-*-*- 
correlation of x1 and x2: 
[1] "0.117"
------------------------------------ 
true parameters are 
[1] "1"   "0.4" "0.3"
------------------------------------
 coefficients are 
(Intercept)          x1          x2 
      1.094       0.426       0.300 
------------------------------------
 standard error are 
(Intercept)          x1          x2 
      0.051       0.048       0.000 
------------------------------------ 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

パラメータ推定値はほぼ真の値になっており、標準誤差を見れば真の値の周りで推定できていることがわかります。

相関が強い場合

x2を作るときの標準偏差が1の場合、つまりx1とx2の相関がそれなりに強い場合の結果が以下です。
相関は0.788で、割と強めです。

-*-*-*-*-* x2_std = 1 *-*-*-*-*-*- 
correlation of x1 and x2: 
[1] "0.788"
------------------------------------ 
true parameters are 
[1] "1"   "0.4" "0.3"
------------------------------------
 coefficients are 
(Intercept)          x1          x2 
      1.011       0.418       0.275 
------------------------------------
 standard error are 
(Intercept)          x1          x2 
      0.050       0.077       0.054 
------------------------------------ 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

意外と大丈夫っぽそうです。意外と。 標準誤差は大きくなっているものの、おいても解釈を大きく変えるほどではないかと思われます。

相関が非常に強い場合

x2を作るときの標準偏差が1の場合、つまりx1とx2の相関が非常に強い場合の結果が以下です。
相関は0.995で、ほぼ1。つよい。

-*-*-*-*-* x2_std = 0.1 *-*-*-*-*-*- 
correlation of x1 and x2: 
[1] "0.995"
------------------------------------ 
true parameters are 
[1] "1"   "0.4" "0.3"
------------------------------------
 coefficients are 
(Intercept)          x1          x2 
      1.055       0.966      -0.260 
------------------------------------
 standard error are 
(Intercept)          x1          x2 
      0.054       0.493       0.492 
------------------------------------ 
-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

これはダメそうです。多重共線性によって、推定値の標準誤差が大きくなり、
推定値も真の値からは大きく外れてしまっています……
解釈をする場合にも、だいぶ裏にある構造とは異なる結果として解釈してしまいそうです。

まとめ

割と大丈夫。0.8くらいあっても平気なときは平気なようです。
ただ、多変量の場合はもうちょっと複雑な問題(3変数以上の連関)などもあるので、一概には言えないかもしれません。
交互作用項を入れる場合、確かに変数間の相関は強くなりがちなのですが、
今回の検証のようにべらぼうに相関関係が強くならない限りは、
推定値が不安定になっているという懸念はそこまで強くないのでは、と思っています。

深夜テンションでしたが今月もう1記事書かないと目標を達成できないので、がんばります。

*1:機械学習の領域でも、例えばFeature importanceの過小評価などでたまに問題になります

*2:あのブレインパッドでは、この本が新人の課題図書となっているようです。

*3:英語がガバガバなのは許してほしい