JINSテックブログ
🐴

多重共線性の一般的な問題点と対処法

2024/12/18に公開

この投稿は、2024年JINSのアドベントカレンダー18日目の記事です。

「多重共線性」がデータ分析の敵になる理由

データ分析をしていると、こんな経験をしたことはありませんか?
モデルの精度は高いのに、個々の変数の結果がどうも納得いかない…。
「この変数が重要だと思ったのに、なぜ有意じゃないの?」
「似たような変数を入れた途端、結果が不安定になった!」

その原因、実は 多重共線性 によるものかもしれません。
多重共線性とは、モデルの説明変数(独立変数)の間に強い相関があることで起こる問題です。この問題が分析に及ぼす影響は大きく、例えば次のようなことが起こります:

各変数の影響を正確に測定できない
回帰係数の信頼性が低下する
モデルの予測力が過信されやすくなる
この記事では、シンプルなシミュレーションを通じて、
「多重共線性がモデルに与える影響」 を実験データで示しながら、
その回避方法について分かりやすく解説していきます。

実際のRコードも紹介しますので、あなたのデータ分析にすぐ役立てることができます。多重共線性を理解して、もっと強い分析スキルを手に入れましょう!

多重共線性のまとめ

説明変数の中に相関が強い変数が含まれると両者の識別が難しくなる

  • 係数の標準誤差が大きくなる
    • t値が小さくなる → 係数が有意になりにくくなる(※ 後で実験を行う)
    • 係数が理論から予想される値と大きく乖離することがある
  • 個別係数のt値が小さいにも関わらず、決定係数が大きくなる
  • データのわずかな変動や観測期間の変更などで係数が大きく変化する

シミュレーションデータを使った実験

実験1a: 多重共線性でない独立なデータを300レコード発生させ、回帰分析を行う

# オブザベーション数
N <- 300

# 多重共線性が起こっていないデータ
set.seed(123)
x1 <- rnorm(N)
x2 <- rnorm(N)
x3 <- rnorm(N)
e <- rnorm(N)
y  <- 1000 + 100*x1 + 80*x2 + 60*x3 + 100*e
df.1 <- data.frame(x1, x2, x3, e, y)

# 回帰分析
model_lm <- lm(y~x1+x2+x3, data=df.1)
summary(model_lm)

# 結果の可視化
par(mfrow=c(1,2))  # グラフを並べる

## 予測値と実測値プロット
plot(df.1$y, predict(model_lm), xlab="Actual", ylab="Predict", xlim=c(200,1800), ylim=c(200,1800))
abline(0, 1, col="red")  # 理想的な直線

## 残差プロット
plot(predict(model_lm), resid(model_lm), xlab="Predict", ylab="Residual")
abline(0, 0, col="red")  # 残差の基準線

シミュレーションの結果を確認すると、想定通りの結果になっていることが確認できますね。

> summary(model_lm)
Call:
lm(formula = y ~ x1 + x2 + x3, data = df.1)
Residuals:
     Min       1Q   Median       3Q      Max
-256.535  -71.289    0.933   56.108  311.957
Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 1002.441
x1           102.728
x2            91.891
x3            64.558
---
5.874  170.65
6.234   16.48
5.961   15.41
5.700   11.33
<2e-16 ***
<2e-16 ***
<2e-16 ***
<2e-16 ***
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 101.7 on 296 degrees of freedom Multiple R-squared: 0.6617, Adjusted R-squared: 0.6583 F-statistic: 193 on 3 and 296 DF, p-value: < 2.2e-16

予測値と実測値のプロットや残差プロットは、回帰分析や予測モデルの性能を評価する際に役立ちます。これらのプロットを併用することで、モデルの性能や改善の余地をより的確に評価できます。

予測値と実測値プロット
実測値(横軸)と予測値(縦軸)をプロットし、それらが一致していればモデルが正確に予測していることを示します。理想的には、データ点が対角線上に分布します。
良い結果の特徴

  • 点が対角線付近に密集している:モデルの予測値が実測値に近い。
  • 点が広がらない:偏差が小さいことを示している。

残差プロット
残差(実測値 - 予測値)を横軸に、予測値または実測値を縦軸にプロットします。モデルの誤差(残差)がどのように分布しているかを確認します。
良い結果の特徴

  • 残差がランダムに散らばる:残差がゼロを中心に均一に散らばり、特定のパターンが見られない場合、モデルが実測値を適切に予測している可能性が高い。
  • 残差の分散が一様:予測値や実測値の大小にかかわらず、残差の広がりが一定であることが望ましい(等分散性)。

    図1:予測値と実測値のプロットや残差プロット(多重共線性を含まないデータ)

実験1b: 多重共線性のあるデータを300レコード発生させ、回帰分析を行う

# データ生成
set.seed(123)
N <- 300

x1 <- rnorm(N)  # 独立変数1
x2 <- x1 + 0.4 * rnorm(N)  # 独立変数2(x1と高い相関を持つ)
x3 <- rnorm(N)  # 独立変数3
e <- rnorm(N)   # 誤差項

y <- 1000 + 100 * x1 + 80 * x2 + 60 * x3 + 100 * e  # 従属変数

df <- data.frame(x1, x2, x3, e, y)

# 回帰分析
model <- lm(y ~ x1 + x2 + x3, data = df)
summary(model)

# 多重共線性の確認
library(car)
vif(model)

# 結果の可視化
par(mfrow=c(1,2))  # グラフを並べる
plot(df$y, predict(model), xlab="Actual", ylab="Predicted", 
     xlim=c(min(df$y), max(df$y)), ylim=c(min(df$y), max(df$y)))
abline(0, 1, col="red")  # 理想的な直線

plot(predict(model), resid(model), xlab="Predicted", ylab="Residuals")
abline(h=0, col="red")  # 残差の基準線

回帰係数は、シミュレーション結果と異なるが、予測結果は問題ない。多重共線性が起こっているx1, x2において、係数の標準誤差が大きくなっている(t値が小さくな る)ことが確認できる。

> # 回帰分析
> model <- lm(y ~ x1 + x2 + x3, data = df)
> summary(model)

Call:
lm(formula = y ~ x1 + x2 + x3, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-256.535  -71.289    0.933   56.108  311.957 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1002.441      5.874 170.650  < 2e-16 ***
x1            73.000     15.791   4.623 5.66e-06 ***
x2           109.728     14.904   7.363 1.80e-12 ***
x3            64.558      5.700  11.327  < 2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 101.7 on 296 degrees of freedom
Multiple R-squared:  0.7688,	Adjusted R-squared:  0.7665 
F-statistic: 328.1 on 3 and 296 DF,  p-value: < 2.2e-16


図2:予測値と実測値のプロットや残差プロット(多重共線性を含むデータ)

多重共線性を含まないデータ(図1)と多重共線性を含むデータ(図2)のプロットを比較すると、予測結果にほとんど違いがないことがわかります。このことから、多重共線性が存在しても、モデルの予測力そのものには必ずしも大きな悪影響が出るわけではありません。その理由は、説明変数間の強い相関が係数間の調整を引き起こし、結果として予測値を一定に保つ働きをするからです。さらに、多重共線性があっても、R²や予測誤差(RMSEなど)のようなモデル評価指標が極端に変化するとは限らないこともポイントです。

実験1c: 実験1bのデータを1000回発生させ、回帰分析を行う

# 必要なパッケージをロード
library(parallel)

# レコード数の設定
N <- 300

# シミュレーション回数の設定
num_simulations <- 1000

# 並列処理の設定
cl <- makeCluster(detectCores() - 1)  # コア数に応じたクラスター作成
clusterExport(cl, varlist=c("N"))  # クラスター内に必要な変数を共有

# シミュレーション関数
simulate_coefficients <- function(i) {
  tryCatch({
    # データ生成
    x1 <- rnorm(N)  # 独立変数1
    x2 <- x1 + 0.4 * rnorm(N)  # 独立変数2(x1に依存)
    x3 <- rnorm(N)  # 独立変数3
    e  <- rnorm(N)  # 誤差項
    y  <- 1000 + 100 * x1 + 80 * x2 + 60 * x3 + 100 * e  # 従属変数

    # データフレーム作成
    df <- data.frame(x1, x2, x3, e, y)

    # 線形回帰モデル
    model <- lm(y ~ x1 + x2 + x3, data=df)

    # 回帰係数を返す
    coef(model)
  }, error=function(e) {
    # エラーが発生した場合はNULLを返す
    warning(sprintf("Simulation %d failed: %s", i, e$message))
    return(NULL)
  })
}

# 並列処理でシミュレーションを実行
results <- parLapply(cl, 1:num_simulations, simulate_coefficients)

# クラスターを停止
stopCluster(cl)

# 結果を行列に変換(NULLを除外)
para <- do.call(rbind, results[!sapply(results, is.null)])

# 列名の設定
colnames(para) <- c("Intercept", "x1", "x2", "x3")

# 結果の要約
summary(para)

# 結果のプロット(例: ヒストグラム)
par(mfrow=c(2, 2))  # 4つのプロットを並べる設定
hist(para[, "Intercept"], main="Intercept", xlab="Value", col="lightblue", border="white")
hist(para[, "x1"], main="x1 Coefficient", xlab="Value", col="lightgreen", border="white")
hist(para[, "x2"], main="x2 Coefficient", xlab="Value", col="lightpink", border="white")
hist(para[, "x3"], main="x3 Coefficient", xlab="Value", col="lightyellow", border="white")

# x1とx2の回帰係数の散布図と相関
par(mfrow=c(1, 1))  # レイアウトを1つに戻す
plot(para[, "x1"], para[, "x2"], 
     main="Scatterplot of x1 and x2 Coefficients", 
     xlab="x1 Coefficient", ylab="x2 Coefficient", 
     pch=16, col=rgb(0.2, 0.4, 0.6, 0.6))

# 回帰係数間の相関係数
correlation <- cor(para[, "x1"], para[, "x2"])
text(min(para[, "x1"]), max(para[, "x2"]), 
     labels=sprintf("Correlation: %.3f", correlation), 
     pos=4, col="red", cex=1.2)



x1とx2の回帰係数を分析すると、平均値はほぼ正しい結果となっていることが確認できます。また、回帰係数を散布図で可視化したところ、x1とx2の回帰係数には強い負の相関が見られました。これは、元の変数x1とx2が正の相関を持つため、回帰モデル内で係数が調整された結果と考えられます。このように、強い多重共線性がある場合、回帰係数間の相関が逆転することがある点に注意が必要です。

実験1d: データ数300から100,000に増やし、実験3と同様のシミュレーションを行う

多重共線性の正しい対処法

  • データ量を増やす
  • 回帰係数に関する追加情報がある場合、それを使う(β1+β2=1)

多重共線性の誤った対処法

  • 説明変数を減らす
  • x1, x2の相関が高い場合、x1 - x2, x2を変数にする
  • 階差を取る
  • データ数を減らす
# レコード数の設定
# N <- 300
# データ数300から100,000に増やす
N <- 100000

データ数を300から100,000に増やした結果、切片やx1, x2, x3の分散が大幅に減少していることが確認できました。これは、サンプルサイズが少ない場合に多重共線性が原因で、回帰係数が真の値から大きく外れることがある一方、サンプルサイズを増やすことで推定結果が真の値の近辺に安定して分布するようになることを示しています。このように、データ数の増加は多重共線性の影響を軽減し、信頼性の高い推定結果を得るための有効な方法といえます。

おまけ:多重共線性を乗り越える!ムーア・ペンローズ逆行列とは?

ムーア・ペンローズ逆行列って何?

ムーア・ペンローズ逆行列は、逆行列が存在しない場合に「擬似的な逆行列」として使える便利な数学ツールです。これを使うと、多重共線性が原因で通常の逆行列が計算できない場合でも、回帰モデルを構築することができます。

どんな場面で役立つの?

多重共線性がある場合: 説明変数が互いに強い相関を持っていると、通常の方法では回帰係数を計算できませんが、ムーア・ペンローズ逆行列を使えば計算が可能になります。
ランク落ちしたデータ: データが不足している場合や変数が冗長な場合でも、この方法で安定した結果を得られます。

実際に使ってみよう(Rの例)

以下は、ムーア・ペンローズ逆行列を使って多重共線性の問題を解決する簡単なRコードです。

# 必要なパッケージ
library(MASS)

# データ生成(多重共線性あり)
set.seed(123)
N <- 300
x1 <- rnorm(N)
x2 <- x1 + 0.01 * rnorm(N)  # x1に強く依存
x3 <- rnorm(N)
y <- 1000 + 100 * x1 + 80 * x2 + 60 * x3 + rnorm(N)

# 設計行列の作成
X <- cbind(1, x1, x2, x3)

# ムーア・ペンローズ逆行列で回帰係数を計算
beta_mp <- ginv(t(X) %*% X) %*% t(X) %*% y
print(round(beta_mp, 3))

# 比較:lm関数を使った回帰
model <- lm(y ~ x1 + x2 + x3)
print(round(coef(model), 3))

注意点
ムーア・ペンローズ逆行列は便利ですが、多重共線性の根本的な解決策ではありません。
基本的には、多重共線性の正しい対処法で紹介した「データ量を増やす」ことを考えてみましょう。

ムーア・ペンローズ逆行列については奥が深く、まだ語り足りない部分も多いです。このテーマについては、また改めて詳しく解説する記事を投稿したいと思いますので、ぜひお楽しみに!

JINSテックブログ
JINSテックブログ

Discussion