多重共線性の仕組みとシミュレーション - 統計学といくつかのよしなしごと

多重共線性の確認

　マーケティングリサーチで重回帰分析を用いる際、各説明変数の係数に関心が行くことは多いと思う。「これこれの説明変数は他に比べて影響度が強いですね、コイツを上げる（下げる）施策を考えましょう*1」などと言って施策の方向性を定めていくという具合である。
　そこで必ず問題になるのがこの多重共線性で、一般に説明変数間に強い相関がある際に偏回帰係数の推定が不安定になるという状況を指す。この数理については『統計的方法のしくみ』の解説が抜群にわかりやすい。そこで多重共線性の仕組みについて、シミュレーションで確かめつつ記事に残しておくことにする。

統計的方法のしくみ―正しく理解するための30の急所

作者: 永田靖
出版社/メーカー: 日科技連出版社
発売日: 1996/10/01
メディア: 単行本
購入: 27人クリック: 226回
この商品を含むブログ (7件) を見る

推定の流れ

　この記事では2つの説明変数で重回帰分析を行うことを考える。回帰モデルは

${\displaystyle y_i = β_0 + β_1x_{1i} + β_2x_{2i} + ε_i,　ε_i ～ N(0, σ^2) }$

で表現される。最小二乗法では残差平方和を最小にするようパラメータを求めるので、残差を

${\displaystyle e_i = y_i - \hat{y_i} }$

としたとき、残差平方和

${\displaystyle S_e = \sum_{i=1}^n{e_i}^{2} = \sum_{i=1}^n{\{y_i - (\hat{β_0} + \hat{β_1}x_{1i} + \hat{β_2}x_{2i})}\}^2 }$

を最小にする $\hat{β_0},\, \hat{β_1},\, \hat{β_2}$ を求める。各パラメータで偏微分して=0と置いて求めれば良いのでそのまま云々とやっていくと、行列を用いて以下のように表現できる。

${\displaystyle \begin{pmatrix} S_{11} & S_{12} \\ S_{12} & S_{22} \\ \end{pmatrix} \begin{pmatrix} \hat{β_1} \\ \hat{β_2} \\ \end{pmatrix} = \begin{pmatrix} S_{1y} \\ S_{2y} \\ \end{pmatrix} }$

ここで $S_{11}, S_{22}$ はそれぞれ $x_1, x_2$ の偏差平方和、 $S_{12}$ は $x_1$ と $x_2$ の偏差積和、 $S_{1y}, S_{2y}$ はそれぞれ $x_1$ と $y$ 、 $x_2$ と $y$ の偏差積和。あとは逆行列を左から掛けてやれば良いので、

${\displaystyle \begin{eqnarray} \begin{pmatrix} \hat{β_1} \\ \hat{β_2} \\ \end{pmatrix} &=& \begin{pmatrix} S_{11} & S_{12} \\ S_{12} & S_{22} \\ \end{pmatrix} ^{-1} \begin{pmatrix} S_{1y} \\ S_{2y} \\ \end{pmatrix} \\ &=& \dfrac{1}{S_{11}S_{22} - S_{12}^{2}} \begin{pmatrix} S_{22} & {-S}_{12} \\ {-S}_{12} & S_{11} \\ \end{pmatrix} \begin{pmatrix} S_{1y} \\ S_{2y} \\ \end{pmatrix} \\ &=& \dfrac{1}{S_{11}S_{22} - S_{12}^{2}} \begin{pmatrix} S_{22}S_{1y}-S_{12}S_{2y} \\ {-S_{12}S_{1y}+S_{11}S_{2y}} \\ \end{pmatrix} \end{eqnarray} }$

ふぅ*2。これで偏回帰係数が求められた。

多重共線性のメカニズムと相関との関係

　先の式を見ると行列に $\displaystyle{1/(S_{11}S_{22}-S_{12}^{2})}$ が掛かっており、これが0のとき（すなわち偏差平方和・積和行列に逆行列が存在しないとき）は解が無数に存在する（不定）か解が存在しない（不能）ということになる。この状況を多重共線性が存在するという。
　では、この $\displaystyle{1/(S_{11}S_{22}-S_{12}^{2})}=0$ という状況と説明変数間の相関の関係がどうなっているかというと、

$\displaystyle{ \begin{eqnarray} S_{11}S_{22}-S_{12}^{2}=0 &\Longleftrightarrow& \dfrac{S_{12}^{2}}{S_{11}S_{22}}=1 \\ &\Longleftrightarrow& r_{x_{1}x_{2}}^{2} = \left\{\dfrac{S_{12}}{\sqrt{S_{11}S_{22}}}\right\}^{2} = 1 \\ &\Longleftrightarrow& r_{x_{1}x_{2}} = ±1 \end{eqnarray} }$

これにより、説明変数間の相関が±1という状況で多重共線性が発生するということがわかった*3。また $\displaystyle{ r_{x_{1}x_{2}}}$ が±1に極めて近ければ $\displaystyle{ S_{11}S_{22}-S_{12}}$ が0に極めて近いということになるので、解の変化幅が大きくなり安定しないということになる*4。

Rでシミュレーション

　ここまでで多重共線性の仕組みが分かったので、実際回帰式の推定にどのように影響するのかをシミュレーションで確かめてみる。シミュレーション内容はこんな感じ。

N=1,000で目的変数 $\displaystyle{y}$ と1つ目の説明変数 $\displaystyle{x_{1}}$ の正規乱数を作る（今回は $\displaystyle{r_{yx_{1}}=0.4}$ で固定しておく*5）

2つ目の説明変数 $\displaystyle{x_{2}}$ については $\displaystyle{r_{x_{1}x_{2}}}$ を0, 0.05, 0.1, ..., 1と0.05ずつ刻み、各々の値に対し100個の $\displaystyle{x_{2}}$ セットを作る

$\displaystyle{r_{x_{1}x_{2}}}$ 毎にパラメータがどのようにぶれるのか確認

#シミュレーション用データの作成
set.seed(221)

#目的変数と説明変数1の作成
y <- rnorm(1000, mean = 0, sd = 1) #目的変数
e <- rnorm(1000, mean = 0, sd = 1) #説明変数作成のために必要なデータ
p=0.4 #yとx1の相関係数を固定する
x1 <- p * y + sqrt(1 - p ^ 2) * e #yとの相関が0.4の説明変数x1を作成

#説明変数2を作成し、シミュレーションを実施
r = seq(0, 1, by = 0.05) #r_12 = 0, 0.05, ..., 1で設定(21段階)

b0 <- matrix(nrow = 100, ncol = 21) #β0を格納する箱。各相関係数に対し100個ずつ
b1 <- matrix(nrow = 100, ncol = 21) #β1を格納する箱。各相関係数に対し100個ずつ
b2 <- matrix(nrow = 100, ncol = 21) #β2を格納する箱。各相関係数に対し100個ずつ

for(i in 1:21){
  stock <- matrix(nrow = 1000, ncol = 100)
  for(j in 1:100){
    e2 = rnorm(1000, mean = 0, sd = 1)
    x2 <- r[i] * x1 + sqrt(1 - r[i] ^ 2) * e2
    stock[, j] <- x2
    
    dat <- as.data.frame(cbind(y, x1, x2 = stock[, j]))
    res <- lm(y ~ ., data = dat)
    b0[j, i] <- res$coefficients[1]
    b1[j, i] <- res$coefficients[2]
    b2[j, i] <- res$coefficients[3]
  }
}

これで $\displaystyle{\hatβ_{0}, \hatβ_{1}, \hatβ_{2}}$ についてそれぞれ、100（ある $\displaystyle{r_{x_{1}x_{2}}}$ 下でのデータセット数）×21（ $\displaystyle{r_{x_{1}x_{2}}}$ の段階数）の値が得られた。例えば $\displaystyle{\hatβ_{1}}$ の箱にはこんな感じで値が入っている。

> head(b1)
        r=0    r=0.05     r=0.1    r=0.15     r=0.2    r=0.25     r=0.3    r=0.35     r=0.4    r=0.45     r=0.5
1 0.4306679 0.4325078 0.4283519 0.4317270 0.4361706 0.4260468 0.4184288 0.4508224 0.4204306 0.4261973 0.4178097
2 0.4322390 0.4346633 0.4311487 0.4308968 0.4269444 0.4343129 0.4436880 0.4380596 0.4259308 0.4228499 0.4535279
3 0.4313652 0.4297881 0.4267221 0.4289261 0.4415001 0.4260835 0.4101387 0.4224972 0.4297717 0.4291471 0.4383937
4 0.4304963 0.4330273 0.4283483 0.4204201 0.4253410 0.4492216 0.4335296 0.4235167 0.4709619 0.4294702 0.4160259
5 0.4322628 0.4335394 0.4328325 0.4331517 0.4354849 0.4317349 0.4314967 0.4338696 0.4191118 0.4356418 0.3976430
6 0.4334116 0.4313089 0.4324587 0.4292935 0.4345609 0.4290567 0.4176773 0.4438420 0.4313572 0.4149576 0.4371784
     r=0.55     r=0.6    r=0.65     r=0.7    r=0.75     r=0.8    r=0.85     r=0.9    r=0.95       r=1
1 0.4111970 0.4093910 0.4339163 0.4184000 0.3734337 0.4756565 0.3922978 0.4544657 0.4121455 0.4322047
2 0.4391379 0.4516133 0.3873439 0.4027791 0.4100673 0.4379951 0.3401908 0.3585950 0.3741827 0.4322047
3 0.4076027 0.4311014 0.4450212 0.4278456 0.4481529 0.4148706 0.4109423 0.5119509 0.4214167 0.4322047
4 0.4324386 0.4308051 0.4209028 0.4485870 0.4614509 0.3932755 0.4713389 0.5392622 0.4512087 0.4322047
5 0.4360566 0.4087040 0.4727114 0.4316713 0.4444006 0.4488128 0.4205936 0.4277047 0.4467345 0.4322047
6 0.4354719 0.4178111 0.4086015 0.4483683 0.4363689 0.3716729 0.3623556 0.4071001 0.6429823 0.4322047

> nrow(b1)
[1] 100

どのような結果になっているか箱ひげ図で確認してみると、 $\displaystyle{r_{x_{1}x_{2}}}$ が1に近づくほどブレ幅が大きくなることがわかる。

#b1の偏回帰係数の図示
b1 <- as.data.frame(b1)
colnames(b1) <- colnames_r

b1_long <- b1 %>% 
  tidyr::gather(key = r, value = b1)

ggplot(b1_long, aes(x = r, y = b1)) +
  geom_boxplot() +
  ggtitle("説明変数間の相関の強弱に伴うb1のバラツキ")

#b2の偏回帰係数の図示
b2 <- as.data.frame(b2)
colnames(b2) <- colnames_r

b2_long <- b2 %>% 
  tidyr::gather(key = r, value = b2)

ggplot(b2_long, aes(x = r, y = b2)) +
  geom_boxplot() +
  ggtitle("説明変数間の相関の強弱に伴うb2のバラツキ")

● $\displaystyle{\hatβ_{1}}$ の結果図示
f:id:yanbow221:20170805033034j:plain
● $\displaystyle{\hatβ_{2}}$ の結果図示
f:id:yanbow221:20170805033052j:plain

ちなみに $\displaystyle{\hatβ_{0}}$ について同じように図示するとこうなる。
● $\displaystyle{\hatβ_{0}}$ の結果図示
f:id:yanbow221:20170805033808j:plain

これは、 $\displaystyle{\hatβ_{0}}$ については
$\displaystyle{ \hat{β_0}=\bar{y}-\hat{β_1}\bar{x_1}-\hat{β_2}\bar{x_2} }$
を解くことで求まるが、 $\displaystyle{\bar{x_1},\bar{x_2}}$ についてはどのデータセットでも同じ値であること、そのうえで $\displaystyle{(\hat{β_1},\hat{β_2})}$ の組み合わせ（足し合わせ）が

> head(b1 + b2)
        r=0    r=0.05     r=0.1    r=0.15     r=0.2    r=0.25     r=0.3
1 0.4062166 0.4154052 0.3839790 0.4027693 0.4471804 0.4732325 0.4681899
2 0.4664639 0.4493094 0.4712963 0.4394332 0.4249584 0.4497007 0.4273785
3 0.4267857 0.4824163 0.4557214 0.4181341 0.4777092 0.4289483 0.4333954
4 0.4379343 0.4109607 0.4335442 0.4378974 0.4731251 0.4555382 0.4145050
5 0.4744823 0.4237386 0.4455729 0.4360844 0.4437070 0.4832751 0.4168798
6 0.4561966 0.4234361 0.4466039 0.3999932 0.4673920 0.4392135 0.4219600
     r=0.35     r=0.4    r=0.45     r=0.5    r=0.55     r=0.6    r=0.65
1 0.4214690 0.3901941 0.4546184 0.4290851 0.4504682 0.4186020 0.4276473
2 0.4738807 0.4074830 0.4196819 0.4178540 0.4350354 0.4318700 0.4474259
3 0.4628155 0.4458085 0.4428629 0.4416364 0.4580332 0.4516505 0.4254354
4 0.4202259 0.4383358 0.4085167 0.4249698 0.4341524 0.4211597 0.4240621
5 0.4645901 0.4214547 0.4467117 0.4362892 0.4431544 0.4008554 0.4301498
6 0.4132844 0.4117565 0.4340837 0.4367685 0.4224682 0.4045233 0.4415290
      r=0.7    r=0.75     r=0.8    r=0.85     r=0.9    r=0.95 r=1
1 0.4308617 0.4302392 0.4172367 0.4228712 0.4474329 0.4418713  NA
2 0.4246184 0.4306662 0.4290404 0.4310337 0.4285018 0.4379655  NA
3 0.4430312 0.4253368 0.4116594 0.4210978 0.4446660 0.4302211  NA
4 0.4359602 0.4359433 0.4402329 0.4316962 0.4229854 0.4388992  NA
5 0.4174786 0.4173395 0.4187569 0.4235048 0.4437448 0.4364285  NA
6 0.4327213 0.4251063 0.4441639 0.4309610 0.4294942 0.4313048  NA

とrによらずほぼ一定であることから理解できる。

終わりに

　ここまで多重共線性の仕組みとその影響を確認してきた。今回は値自体に意味はないため結果の実際的な解釈は行っていないが、多重共線性は係数の大きさや符号を不安定なものにする危険があることがよくわかった。オーソドックスな回帰でも多重共線性をはじめいくつか落とし穴や確認すべき点があるので、結果に基づいた意思決定を行う際は十分注意して分析していこうと思う次第です。

*1:説明変数と目的変数は因果性を示すものではないが、実務的にはそこに因果性を認めて設計することが多いと思う

*2:あまり慣れてないのでLaTeX書くの疲れる

*3:相関係数が±1ということは $\displaystyle{x_{1}とx_{2}}$ が一直線上に並んでいるということ。共通の直線上にあるから「共線」らしい