どうも昨年末にあちこちで多重共線性についての議論がなされていたようなんですが、些事にかまけていた僕はすっかりそのウェーブに乗り損ねてしまっていたのでした。そこで、今年最初の記事では遅ればせながらそのウェーブに乗る形で、また今までに学んだり調べてきたりしてきたことの備忘録も兼ねて、多重共線性についてまとめてみようと思います。
多重共線性とは
このブログの読者の皆様なら多重共線性が何であるかはご存知かと思いますが、一応簡単な例だけ挙げておきます。ここではサンプルデータとしてBoston Housingを使います*1。
オールドタイプなコーディングで恐縮ですが、Rで適当に前処理した上で線形回帰するとこんな感じになります。
d <- read.csv('https://raw.githubusercontent.com/ozt-ca/tjo.hatenablog.samples/refs/heads/master/r_samples/public_lib/jp/R/housing.csv', sep = '\t') d$CHAS <- as.factor(d$CHAS) d$RAD <- as.integer(d$RAD) summary(d) #R> CRIM ZN INDUS CHAS NOX RM #R> Min. : 0.00632 Min. : 0.00 Min. : 0.46 0:471 Min. :0.3850 Min. :3.561 #R> 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19 1: 35 1st Qu.:0.4490 1st Qu.:5.886 #R> Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.5380 Median :6.208 #R> Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.5547 Mean :6.285 #R> 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.6240 3rd Qu.:6.623 #R> Max. :88.97620 Max. :100.00 Max. :27.74 Max. :0.8710 Max. :8.780 #R> AGE DIS RAD TAX PTRATIO #R> Min. : 2.90 Min. : 1.130 Min. : 1.000 Min. :187.0 Min. :12.60 #R> 1st Qu.: 45.02 1st Qu.: 2.100 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 #R> Median : 77.50 Median : 3.207 Median : 5.000 Median :330.0 Median :19.05 #R> Mean : 68.57 Mean : 3.795 Mean : 9.549 Mean :408.2 Mean :18.46 #R> 3rd Qu.: 94.08 3rd Qu.: 5.188 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 #R> Max. :100.00 Max. :12.127 Max. :24.000 Max. :711.0 Max. :22.00 #R> B LSTAT MEDV #R> Min. : 0.32 Min. : 1.73 Min. : 5.00 #R> 1st Qu.:375.38 1st Qu.: 6.95 1st Qu.:17.02 #R> Median :391.44 Median :11.36 Median :21.20 #R> Mean :356.67 Mean :12.65 Mean :22.53 #R> 3rd Qu.:396.23 3rd Qu.:16.95 3rd Qu.:25.00 #R> Max. :396.90 Max. :37.97 Max. :50.00 d.lm <- lm(MEDV ~ ., d) summary(d.lm) #R> #R> Call: #R> lm(formula = MEDV ~ ., data = d) #R> #R> Residuals: #R> Min 1Q Median 3Q Max #R> -15.595 -2.730 -0.518 1.777 26.199 #R> #R> Coefficients: #R> Estimate Std. Error t value Pr(>|t|) #R> (Intercept) 3.646e+01 5.103e+00 7.144 3.28e-12 *** #R> CRIM -1.080e-01 3.286e-02 -3.287 0.001087 ** #R> ZN 4.642e-02 1.373e-02 3.382 0.000778 *** #R> INDUS 2.056e-02 6.150e-02 0.334 0.738288 #R> CHAS1 2.687e+00 8.616e-01 3.118 0.001925 ** #R> NOX -1.777e+01 3.820e+00 -4.651 4.25e-06 *** #R> RM 3.810e+00 4.179e-01 9.116 < 2e-16 *** #R> AGE 6.922e-04 1.321e-02 0.052 0.958229 #R> DIS -1.476e+00 1.995e-01 -7.398 6.01e-13 *** #R> RAD 3.060e-01 6.635e-02 4.613 5.07e-06 *** #R> TAX -1.233e-02 3.760e-03 -3.280 0.001112 ** #R> PTRATIO -9.527e-01 1.308e-01 -7.283 1.31e-12 *** #R> B 9.312e-03 2.686e-03 3.467 0.000573 *** #R> LSTAT -5.248e-01 5.072e-02 -10.347 < 2e-16 *** #R> --- #R> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #R> #R> Residual standard error: 4.745 on 492 degrees of freedom #R> Multiple R-squared: 0.7406, Adjusted R-squared: 0.7338 #R> F-statistic: 108.1 on 13 and 492 DF, p-value: < 2.2e-16
連続値の説明変数が大半なので、試しに"CRIM"に適当な正規分布ノイズを加えたものを"X1", "NOX"にノイズを加えたものを"X2"として、説明変数に加えます。その上で同じことをやってみましょう。
d1 <- d d1$X1 <- d1$CRIM + rnorm(nrow(d1), 0, sd(d1$CRIM) / 10) d1$X2 <- d1$NOX + rnorm(nrow(d1), 0, sd(d1$NOX) / 10) summary(d1) # 省略 d1.lm <- lm(MEDV ~ ., d1) summary(d1.lm) #R> #R> Call: #R> lm(formula = MEDV ~ ., data = d1) #R> #R> Residuals: #R> Min 1Q Median 3Q Max #R> -15.6377 -2.7417 -0.4871 1.8472 26.1112 #R> #R> Coefficients: #R> Estimate Std. Error t value Pr(>|t|) #R> (Intercept) 3.626e+01 5.121e+00 7.081 4.98e-12 *** #R> CRIM -2.989e-01 2.404e-01 -1.243 0.214383 # ←ココ #R> ZN 4.632e-02 1.375e-02 3.368 0.000818 *** #R> INDUS 1.788e-02 6.167e-02 0.290 0.772004 #R> CHAS1 2.699e+00 8.669e-01 3.114 0.001957 ** #R> NOX -1.908e+01 1.895e+01 -1.007 0.314541 # ←ココ #R> RM 3.800e+00 4.190e-01 9.068 < 2e-16 *** #R> AGE 8.219e-04 1.324e-02 0.062 0.950514 #R> DIS -1.467e+00 2.004e-01 -7.320 1.02e-12 *** #R> RAD 3.037e-01 6.650e-02 4.567 6.27e-06 *** #R> TAX -1.219e-02 3.771e-03 -3.231 0.001315 ** #R> PTRATIO -9.503e-01 1.316e-01 -7.222 1.97e-12 *** #R> B 9.383e-03 2.692e-03 3.485 0.000536 *** #R> LSTAT -5.264e-01 5.083e-02 -10.355 < 2e-16 *** #R> X1 1.915e-01 2.388e-01 0.802 0.422984 #R> X2 1.590e+00 1.851e+01 0.086 0.931546 #R> --- #R> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #R> #R> Residual standard error: 4.752 on 490 degrees of freedom #R> Multiple R-squared: 0.741, Adjusted R-squared: 0.7331 #R> F-statistic: 93.45 on 15 and 490 DF, p-value: < 2.2e-16
物の見事に"CRIM"と"NOX"の回帰係数の推定結果がおかしなことになっています。どちらもp < 0.05であったはずがそうではなくなり、回帰係数そのものもバイアスがかかって異なる値になってしまいました。これこそが多重共線性(multicolinearity)です。
その定義については佐和本を参照するのが単純かつ確実ですが、要は「互いに相関が強い2つ以上の説明変数があると互いに識別するのが難しくなる」ということですね。本質的には回帰分析における逆行列演算のランク落ちに帰着させられますが、これは程度問題であって「ある程度以上互いに相関係数が高い説明変数同士の回帰係数の推定は不安定になる」と理解するのが妥当かと思われます。
なお、多重共線性を定量化する指標としてVIF (variance inflation factor)があり、Rではcar::vifで算出できます。そこで2つのモデルそれぞれのVIFを算出してみると、以下のような結果が得られます。
car::vif(d.lm) #R> CRIM ZN INDUS CHAS NOX RM AGE DIS RAD TAX #R> 1.792192 2.298758 3.991596 1.073995 4.393720 1.933744 3.100826 3.955945 7.484496 9.008554 #R> PTRATIO B LSTAT #R> 1.799084 1.348521 2.941491 car::vif(d1.lm) #R> CRIM ZN INDUS CHAS NOX RM AGE DIS #R> 95.646892 2.301523 4.003448 1.084371 107.809621 1.938677 3.104448 3.984259 #R> RAD TAX PTRATIO B LSTAT X1 X2 #R> 7.499034 9.034075 1.814976 1.351166 2.947177 94.988828 103.959498
オリジナルのモデルでは概ね穏当な値である一方、X1 & X2を加えたモデルでは"CRIM"と"NOX"のVIFが100近くと極めて異常な値に跳ね上がっており、X1 & X2が強い多重共線性をもたらしていることが見て取れます。
「予測」が目的なら多重共線性があっても問題ではない
そんな厄介者の多重共線性ですが、かつては「百計講じて排除すべき」とされてきました。しかし、同じ多変量解析でも回帰分析ではなくNNや樹木モデルなど機械学習系の手法が多用されるようになった現在では、必ずしもそうではないんですね。一般に、「予測」においては多重共線性は特に問題にならないとされます。
と言いますか、実はOLS線形回帰ですらも「予測」が目的なら多重共線性はあっても構わないとされます。これに関する議論はかなり以前からweb上を探せば見つかり、Cross Validatedにもそのものズバリの問答があります。
cv::cv(d.lm) # オリジナル R RNG seed set to 385600 #R> cross-validation criterion (mse) = 23.61569 cv::cv(d1.lm) # 多重共線性が存在するケース R RNG seed set to 458688 #R> cross-validation criterion (mse) = 23.93573
なお、前掲のBoston Housingの回帰分析の例でもcv::cvで10-fold CV MSEを比較してみると、ほとんど差がないことが分かります。多重共線性は「予測」には影響しないと見て良さそうです。
また、特に樹木モデルにおいて多重共線性が問題にならない理由についても良い問答がちらほらあります。要は線形モデルとは異なり樹木モデルはその構造上説明変数を逐次処理していくので、互いに相関の強い変数が複数あってもそれらから同時に影響を受けにくい、ということなんですね。
「説明」が目的なら多重共線性に対処する必要があるケース「も」ある
実は、前掲したHyndmanのブログからリンクされている論文"To Explain or To Predict?"の中に、以下のような端的なまとめが書かれています。僕個人はこれで全て言い尽くされていると思います。
Multicollinearity is not a problem unless either (i) the individual regression coefficients are of interest, or (ii) attempts are made to isolate the contribution of one explanatory variable to Y, without the influence of the other explanatory variables. Multicollinearity will not affect the ability of the model to predict.
（多重共線性は、(i) 個々の回帰係数が重要である場合、または (ii) 他の説明変数の影響なしに、1つの説明変数のYへの寄与を分離しようとする場合を除いて、問題にはならない。多重共線性は、モデルの予測能力に影響を与えない。）
なのですが、これは裏を返すと「個々の回帰係数が重要である場合」と「他の説明変数の影響なしに1つの説明変数の目的変数への寄与を分離しようとする場合」には、何かしらの多重共線性への対処が必要になるということです。
即ち前節との比較で言えば、「説明」が目的なら多重共線性に対処しなければならないということですね。ただ、現実には「関心ある特定の説明変数」に多重共線関係にある他の変数がなければ、放っておいても問題ないとされます。実際、冒頭のBoston Housingの回帰分析の多重共線性の例でも、例えば"RAD"や"TAX"などX1 & X2と相関の低い変数では回帰係数もp値もほとんど影響を受けていません。
ちなみに、先述したように樹木モデルは多重共線性の影響を受けにくいとされますが、SHAPやfeature importanceなど解釈性つまり「説明」のための指標を重視する場合には、やはり多重共線性への対処が必要なようです。
実務における多重共線性への対処法
ということで、「説明」が目的なら多重共線性への対処が必要だということが分かったわけですが、ビジネス実務におけるデータ分析でその対処が必要な回帰分析の代表例を挙げるとしたらやはりMMM (Media/Marketing Mix Models)でしょうか。理由は簡単で、MMMこそが「個々の回帰係数が重要」な分析の典型例だからです。
その対処法には様々な流儀がありますが、実務上は「VIFで多重共線関係にある説明変数を特定」した上で「ドメイン知識（変数の生成過程に関する情報）に基づいて削除orマージする変数を決定」することをお薦めしたいです。例えば、冒頭のBoston Housingの例で言えば、以下のVIFの算出結果から"CRIM", "NOX", X1, X2が多重共線関係にあることが分かります。
car::vif(d1.lm) #R> CRIM ZN INDUS CHAS NOX RM AGE DIS #R> 95.646892 2.301523 4.003448 1.084371 107.809621 1.938677 3.104448 3.984259 #R> RAD TAX PTRATIO B LSTAT X1 X2 #R> 7.499034 9.034075 1.814976 1.351166 2.947177 94.988828 103.959498
当然ですが、我々は既にデータ定義から"CRIM"（犯罪率）と"NOX"（大気汚染指標）は把握しているわけです。故に、X1とX2にそれを上回る重要性がなければ、合算平均するなどしてマージするか、さもなくば削除するというのが妥当でしょう。
また、実務上はステークホルダーの意向などの理由で「多重共線性が予め見込まれる」「残すべき変数が予め決まっている」ことも多いです。そのような場合も、VIFで多重共線関係にある変数を特定した上で、それらの実務上の事情や要件を勘案した上でマージor削除する変数を決定するのが無難でしょう。
特に広告マーケティング分野では、マス広告の出稿タイミングにシンクロさせて他カテゴリの広告も出稿させていることが多く、最初からがっつり多重共線性が乗っているというケースは珍しくありません。そのような場合は、分析にもとめられる要件と注意深く照らし合わせながら慎重にマージor削除する変数を決める必要があります。
なお、L1正則化(Lasso)で変数選択するという考え方もありますが、基本的にはお薦めしません。例えば、冒頭のBoston Housingの回帰分析における多重共線性の例にL1正則化を適用すると、X1 & X2を生成する際の乱数シード次第では以下のような結果になります。
d1.cv.glmnet <- glmnet::cv.glmnet(as.matrix(d1[, -13]), as.matrix(d1[, 13]), family = 'gaussian', alpha = 1) coef(d1.cv.glmnet, s = d1.cv.glmnet$lambda.min) #R> 16 x 1 sparse Matrix of class "dgCMatrix" #R> s1 #R> (Intercept) 36.006548447 #R> CRIM . #R> ZN 0.024192608 #R> INDUS 0.069390086 #R> CHAS 0.057447332 #R> NOX . #R> RM -2.306342512 #R> AGE 0.072678003 #R> DIS -0.347259995 #R> RAD 0.114193893 #R> TAX -0.003444762 #R> PTRATIO -0.213520854 #R> B -0.003507256 #R> MEDV -0.334299691 #R> X1 0.048139632 #R> X2 -1.283534711
このケースでは物の見事に、後から人為的に加えたX1 & X2が残った一方で、オリジナルの説明変数であるはずの"CRIM" & "NOX"が削られてしまっています。このように、L1正則化は多重共線性のある変数同士の中からランダムに削ってしまい、実際の個々の変数の生成過程や変数同士の関係は考慮してくれないので注意が必要です。AIC/BICなどによるstepwise法は、そもそもの変数選択の仕方に問題があるケースが目立つので、ここでは割愛します。
ちなみに、PCAなどの次元集約型の手法で変数をマージするアプローチもありますが、個人的には「説明変数が多過ぎる＆ドメイン知識に乏しい」場合に限った方が良いと考えています。可能な限り、個々の変数の生成過程や分析そのものの実務上の要件を勘案して対処するべきでしょう。
最後に
以上の議論は、以前に某所*2にまとめた内容を、改めて再構成したものです。個人的には良い勉強になったと感じている次第です。
*1:この有名なデータセットには微妙にpolitically incorrectな要素があり、近年はあまり用いられないようです
*2:LightGBMの出力結果に対して、SHAPを使った場合、SHAPで算出した寄与度は多重共線性の影響を受けるのでしょうか？ | mond
