ちょっと前に、回帰分析における多重共線性に関する解説記事を2本ほど書いたわけですが。
多重共線性そのものの問題点はこれでもかと論じている割に、その対処法についてはあまり触れていなかったなと気付いたのでした。ということで、今回の記事では遅ればせながら多重共線性への対処法をある程度網羅的に挙げていこうと思います。
- データセット
- VIFで多重共線性に寄与する変数を特定した上で削除orマージ（第一選択）
- PCAで変数を削除orマージ（要件次第）
- L1正則化で変数選択（要件次第）
- L2正則化で回帰係数同士のバランスをとる（非推奨）
- ベイズ回帰で事前分布を設定することで多重共線性によるバイアスを軽減させる（事前分布の蓋然性次第）
- コメントなど
データセット
まず、x1-5の5つの説明変数から成るデータセットを用意します。真の回帰係数はそれぞれ1, 2, -2, 5, -3.5とし、サンプルサイズは100とします。面倒だったので、要件だけ用意した上でGemini 2.5 Proに作らせました（笑）。良い時代になったものです。
# -------------------------------------------------------------------------- # ベースとなるデータの生成 # -------------------------------------------------------------------------- # 再現性のための乱数シード設定 set.seed(123) # サンプルサイズ n <- 100 # 説明変数の数 p <- 5 # 真の回帰係数 beta <- c(1.0, 2.0, -2.0, 5.0, -3.5) # 説明変数の生成 (行列として) X <- matrix(rnorm(n * p), nrow = n, ncol = p) # 誤差項の生成 epsilon <- rnorm(n, mean = 0, sd = 2) # 目的変数の生成 y <- X %*% beta + epsilon # データをデータフレームにまとめる sim_data <- as.data.frame(X) colnames(sim_data) <- paste0("x", 1:p) sim_data$y <- y # 線形回帰分析を行い、結果を出力する summary(lm(y ~., sim_data))
#R> Call: #R> lm(formula = y ~ ., data = sim_data) #R> #R> Residuals: #R> Min 1Q Median 3Q Max #R> -5.6518 -1.0858 0.0935 1.1618 5.1065 #R> #R> Coefficients: #R> Estimate Std. Error t value Pr(>|t|) #R> (Intercept) -0.1065 0.1909 -0.558 0.578 #R> x1 0.9825 0.2111 4.654 1.07e-05 *** #R> x2 2.2849 0.1954 11.696 < 2e-16 *** #R> x3 -1.9723 0.1986 -9.929 2.60e-16 *** #R> x4 4.8341 0.1802 26.822 < 2e-16 *** #R> x5 -3.0773 0.1945 -15.825 < 2e-16 *** #R> --- #R> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #R> #R> Residual standard error: 1.856 on 94 degrees of freedom #R> Multiple R-squared: 0.9364, Adjusted R-squared: 0.9331 #R> F-statistic: 277 on 5 and 94 DF, p-value: < 2.2e-16
まず、本来のデータセットに対して線形回帰した結果がこちら。加えた誤差項に依存する形で若干真の値からは推定値がずれていますが、とりあえずこれで良しとしましょう。次に、x2とx5にわずかなノイズを加えただけの新たな変数として、x6-7を追加します。
# -------------------------------------------------------------------------- # 多重共線性を持つ変数 (x6, x7) の生成 # -------------------------------------------------------------------------- # x6 の生成 (x2 に基づく) # シードを設定し、標準偏差の小さい正規乱数を加える set.seed(456) noise1 <- rnorm(n, mean = 0, sd = 0.1) # sd=0.1 の小さいノイズ x6 <- sim_data$x2 + noise1 # x7 の生成 (x5 に基づく) # 異なるシードを設定し、標準偏差の小さい正規乱数を加える set.seed(789) noise2 <- rnorm(n, mean = 0, sd = 0.1) # sd=0.1 の小さいノイズ x7 <- sim_data$x5 + noise2 # 元のデータに新しい変数を追加 multicol_data <- sim_data multicol_data$x6 <- x6 multicol_data$x7 <- x7 # 生成されたデータの最初の6行を表示 print(head(multicol_data)) # -------------------------------------------------------------------------- # 多重共線性の影響を確認 # -------------------------------------------------------------------------- # 新しいデータで回帰モデルを構築 model_multicol <- lm(y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7, data = multicol_data) # 結果の要約を表示 print(summary(model_multicol)) # VIFを計算 print(car::vif(model_multicol))
#R> Call: #R> lm(formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7, data = multicol_data) #R> #R> Residuals: #R> Min 1Q Median 3Q Max #R> -5.4663 -0.9905 0.0794 1.2047 5.2527 #R> #R> Coefficients: #R> Estimate Std. Error t value Pr(>|t|) #R> (Intercept) -0.1346 0.1919 -0.701 0.4850 #R> x1 1.0012 0.2111 4.743 7.66e-06 *** #R> x2 -0.2544 1.9065 -0.133 0.8942 #R> x3 -1.9681 0.1982 -9.932 3.18e-16 *** #R> x4 4.7738 0.1843 25.899 < 2e-16 *** #R> x5 -5.1358 1.9861 -2.586 0.0113 * #R> x6 2.5495 1.9110 1.334 0.1855 #R> x7 2.0782 1.9752 1.052 0.2955 #R> --- #R> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #R> #R> Residual standard error: 1.852 on 92 degrees of freedom #R> Multiple R-squared: 0.9381, Adjusted R-squared: 0.9334 #R> F-statistic: 199.2 on 7 and 92 DF, p-value: < 2.2e-16
そして、多重共線性を生じさせたデータセットに対して線形回帰した結果がこちら。物の見事にx2, x5とも回帰係数の値にバイアスが乗ってずれてしまっており、x2に至っては回帰係数の正負すら逆になっています。これは、x2とx6、x5とx7が極端に高い相関を持つことで、線形回帰における逆行列演算の結果が歪められたことによって起きている現象です。なお、これは線形回帰に限らずGLMやベイズ回帰など、原理的に同様の計算を行う回帰分析全てにおいて同じように起き得る問題です*1。
以前の記事でも指摘したように、多重共線性は「予測」には影響しませんが「説明」には大きく影響を及ぼします。そこで、「説明」を目的とすることが多い種々の回帰分析においては、多重共線性に遭遇した際は何かしらの方法で対処する必要があります。以下に、その対処法を挙げていこうと思います。
VIFで多重共線性に寄与する変数を特定した上で削除orマージ（第一選択）
現実的には、最も堅実なのが「VIF (variance inflation factor)で多重共線性に寄与する変数を特定した上で削除orマージする」というアプローチです。例えば、今回のデータセットに対してVIFを計算すると以下のようになります。
> print(car::vif(model_multicol)) #R> x1 x2 x3 x4 x5 x6 x7 #R> 1.071904 98.135113 1.023030 1.058597 111.489219 98.316420 111.880168
大体VIF > 7で多重共線性があると判定されますが、7どころか100前後という値をx2, x5, x6, x7が示していることが分かります。我々はこのデータセットの生成過程を知っているので「x6-7を削除すれば良い」とすぐ判断できますが、知らない場合は削除ではなくマージ（大抵は足し合わせor平均）することが多いと思われます。なお、(x2 + x6) / 2及び(x5 + x7) / 2で対処した場合、以下のような結果になります。
merged_data <- data.frame(y = multicol_data$y, x1 = multicol_data$x1, x2 = (multicol_data$x2 + multicol_data$x6) / 2, x3 = multicol_data$x3, x4 = multicol_data$x4, x5 = (multicol_data$x5 + multicol_data$x7) / 2) summary(lm(y ~., merged_data)) #R> Call: #R> lm(formula = y ~ ., data = merged_data) #R> #R> Residuals: #R> Min 1Q Median 3Q Max #R> -5.6413 -0.9014 0.0370 1.1701 5.3864 #R> #R> Coefficients: #R> Estimate Std. Error t value Pr(>|t|) #R> (Intercept) -0.1262 0.1918 -0.658 0.512 #R> x1 0.9787 0.2124 4.608 1.28e-05 *** #R> x2 2.3102 0.1969 11.731 < 2e-16 *** #R> x3 -1.9710 0.1997 -9.869 3.49e-16 *** #R> x4 4.8473 0.1812 26.748 < 2e-16 *** #R> x5 -3.0455 0.1952 -15.603 < 2e-16 *** #R> --- #R> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 #R> #R> Residual standard error: 1.866 on 94 degrees of freedom #R> Multiple R-squared: 0.9357, Adjusted R-squared: 0.9323 #R> F-statistic: 273.8 on 5 and 94 DF, p-value: < 2.2e-16 car::vif(lm(y ~., merged_data)) #R> x1 x2 x3 x4 x5 #R> 1.068414 1.026376 1.022840 1.007182 1.065008
線形回帰してみると、多少真の値からはずれているもののx2の回帰係数は正に戻り、x5の回帰係数も絶対値が上振れしていたのが抑えられているのが分かります。要件次第ではありますが、実務的にはこのアプローチが最も堅実かと思われます。
PCAで変数を削除orマージ（要件次第）
次に挙げるアプローチは「相関の高い変数同士のグループを定量的に見つけ出して削除orマージ」というものです。手法は色々考えられますが、PCAが最も無難でしょう。
biplot(princomp(multicol_data[, -6]))
こんな感じでx2 & x6、x5 & x7のペアが怪しいことが一目で分かります。ただ、今回生成したデータセットだと他にも微妙に相関の高い変数のグループが混じっており、ドメイン知識（データセットの生成過程の情報）がないと対処を決められないと思われます。その意味では、VIFベースのアプローチ同様にドメイン知識が必要だと言えそうです。
L1正則化で変数選択（要件次第）
これまでは個々の変数を一つひとつ吟味するアプローチを紹介してきましたが、ここからはアルゴリズミックなアプローチを挙げていきます。まず挙げるのがL1正則化（Lasso回帰）。L1ノルムをパラメータ推定に正則化項として加えることで、モデルの予測精度に寄与しない変数を削除する手法です。交差検証と組み合わせれば、ほぼ全自動で変数選択できます。
multicol_data.cv.glmnet1 <- glmnet::cv.glmnet(as.matrix(multicol_data[, -6]), as.matrix(multicol_data[, 6]), family = "gaussian", alpha = 1) coef(multicol_data.cv.glmnet1, s = multicol_data.cv.glmnet1$lambda.min) #R> 8 x 1 sparse Matrix of class "dgCMatrix" #R> s=0.0214674 #R> (Intercept) -0.141012413 #R> x1 0.967067128 #R> x2 0.007042437 #R> x3 -1.949233677 #R> x4 4.791074827 #R> x5 -3.039734696 #R> x6 2.276163921 #R> x7 .
glmnet::cv.glmnetで実際にやってみた結果がこちらなんですが……x7だけが削除された一方でx6が残ってしまっており、結果としてx2を初めとしてどの回帰係数にもバイアスが乗ってしまっています。基本的に正則化＋交差検証はアルゴリズミックなアプローチであり、ドメイン知識は一切加味されないという点に注意が必要です。よって「アルゴリズミックな対処で構わない」ケースに限るべきだというのが、個人的な意見です。
L2正則化で回帰係数同士のバランスをとる（非推奨）
これは佐和本でも紹介されている、L2ノルムをパラメータ推定に正則化項として加えることで、パラメータ同士のバランスをとる手法です（L2正則化：Ridge回帰）。これまた交差検証と組み合わせることで、ほぼ全自動で対処できます。
multicol_data.cv.glmnet2 <- glmnet::cv.glmnet(as.matrix(multicol_data[, -6]), as.matrix(multicol_data[, 6]), family = "gaussian", alpha = 0) coef(multicol_data.cv.glmnet2, s = multicol_data.cv.glmnet2$lambda.min) #R> 8 x 1 sparse Matrix of class "dgCMatrix" #R> s=0.5195386 #R> (Intercept) -0.1630073 #R> x1 0.9284621 #R> x2 1.0483262 #R> x3 -1.8539534 #R> x4 4.5201612 #R> x5 -1.7137011 #R> x6 1.2030075 #R> x7 -1.2417390
……なのですが、x6-7を残しながら全変数の回帰係数同士のバランスを取るように推定しているので、当然のように全ての回帰係数にしっかりバイアスが乗ってしまっています。これでは「説明」の用を為さないでしょう。ということで、色々な意見があろうかとは思いますが、個人的にはL2正則化で多重共線性に対処するのは「推奨しない」という立場です。
ベイズ回帰で事前分布を設定することで多重共線性によるバイアスを軽減させる（事前分布の蓋然性次第）
これは、アンコレさんのブログ記事で知りました。
要は「真の回帰係数の情報を事前分布として与えることで多重共線性によるバイアスの影響を抑える」という考え方ですね。念のためベイズ統計学の専門家にも聞いてみたところでは「アリ」とのことでしたので、試してみる価値はありそうです。
で、早速やってみようと思ったんですが、RStanだと計算時間が長くなりそうだなと思ったので、アンコレさんに倣ってMCMCpack::MCMCregressを使うことにしました。事前分布の与え方については、下記コード中にもありますが後述します。
# -------------------------------------------------------------------------- # ベイズ回帰：事前分布のパラメータを設定 # -------------------------------------------------------------------------- # モデル式: y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7 # 係数の順序: (Intercept), x1, x2, x3, x4, x5, x6, x7 # 事前分布の平均ベクトル (b0) # x2の係数(3番目)に 2.0, x5の係数(6番目)に -3.5 を設定 b0 <- c(0, 0, 2.0, 0, 0, -3.5, 0, 0) # 事前分布の精度行列 (B0) # 8x8の対角行列を作成し、対角成分を0.5とする（分散2に相当） # これにより、やや強い事前分布を与える B0 <- diag(8) * 0.5 # -------------------------------------------------------------------------- # MCMCregressによるベイズ線形回帰分析の実行 # -------------------------------------------------------------------------- # 再現性のためのシード設定 set.seed(111) # MCMCの実行 # burnin: 初期サンプルのうち、捨てる数 # mcmc: 最終的に保持するサンプルの数 model_bayes <- MCMCpack::MCMCregress( formula = y ~ x1 + x2 + x3 + x4 + x5 + x6 + x7, data = multicol_data, b0 = b0, # 事前分布の平均 B0 = B0, # 事前分布の精度 burnin = 2000, mcmc = 20000, verbose = 1000 # 1000回ごとに進捗を表示 ) # -------------------------------------------------------------------------- # 結果の確認 # -------------------------------------------------------------------------- # 結果の要約を表示 print(summary(model_bayes)) # 結果をプロットして確認 # 左側: 各係数の事後分布の密度プロット # 右側: MCMCサンプルのトレースプロット（サンプリングがうまく収束しているかを確認） plot(model_bayes)
#R> Iterations = 2001:22000 #R> Thinning interval = 1 #R> Number of chains = 1 #R> Sample size per chain = 20000 #R> #R> 1. Empirical mean and standard deviation for each variable, #R> plus standard error of the mean: #R> #R> Mean SD Naive SE Time-series SE #R> (Intercept) -0.1143 0.1916 0.001355 0.001355 #R> x1 0.9664 0.2109 0.001491 0.001491 #R> x2 1.6527 0.8815 0.006233 0.006233 #R> x3 -1.9418 0.1991 0.001408 0.001408 #R> x4 4.7402 0.1810 0.001280 0.001280 #R> x5 -3.6325 0.9035 0.006388 0.006340 #R> x6 0.6412 0.8828 0.006243 0.006243 #R> x7 0.5578 0.8986 0.006354 0.006350 #R> sigma2 3.4947 0.5211 0.003685 0.003932 #R> #R> 2. Quantiles for each variable: #R> #R> 2.5% 25% 50% 75% 97.5% #R> (Intercept) -0.49147 -0.24186 -0.1144 0.01164 0.2646 #R> x1 0.55185 0.82325 0.9669 1.10723 1.3837 #R> x2 -0.06377 1.05149 1.6547 2.24029 3.3989 #R> x3 -2.32957 -2.07666 -1.9424 -1.80940 -1.5418 #R> x4 4.38774 4.61784 4.7426 4.86155 5.0934 #R> x5 -5.40923 -4.24630 -3.6349 -3.02671 -1.8752 #R> x6 -1.10183 0.05615 0.6474 1.24342 2.3673 #R> x7 -1.19396 -0.05645 0.5596 1.17043 2.3258 #R> sigma2 2.62100 3.12232 3.4417 3.80502 4.6572
ということで、「x6-7を削除する」には敵いませんが、このアプローチでは残りの全ての他のアプローチと比べてもかなり真の結果に近い推定結果を返していることが分かるかと思います。
ただ、今回は分散2とかなり「強い」即ち裾の狭く尖った事前分布を与えた結果であるという点に注意が必要です。実際、分散100と弱い事前分布を与えると、x2の回帰係数が負のままになることが分かっており、事前分布への依存性が高い*2データセットであった可能性も否定できません。
その意味で言えば、「ベイズ回帰にして事前分布を設定することで多重共線性を軽減させる」というアプローチで最も重要なのは「事前分布の蓋然性」だということになるのでしょう。それは例えば介入実験などで予め真に近い値が得られているというような、ごくごく限られたシチュエーションでのみ達成されるものであろうと考えられます。
特に近年マーケティング分野で盛んなMMMではベイズ的手法が注目を集めており、それはイコール事前分布を援用できるということでもあるわけですが、そのためには予めきちんと条件統制されたマーケティング実験でCPAなどの「真の値」を「信頼できる形で」得ておくことが、必須とされそうです。
コメントなど
一通り試してみて思ったのが「意外とVIFベースで変数をマージするアプローチでもバイアスが残るんだ」ということと、「ベイズ回帰で事前分布を使うと思った以上にバイアスを補正できるんだ」ということ。これは個人的には望外の収穫であり、今後の回帰分析に活かしていこうと思った次第です。
……ところで。今回初めてこのブログの記事を書くに当たってGemini 2.5 Proにvibe codingで大半のRコードを書かせてみたんですが、この程度のレベル感なら間違わずに書いてくれるので大変に楽ですね〜。Geminiのおかげでブログの更新がめちゃくちゃ捗りそうです。
