見出し画像

【Python入門:主成分分析編#6】PCAを回帰分析に活かす!主成分回帰(PCR)で多重共線性を攻略しよう

「変数同士が似すぎていて、回帰分析がうまくいかない……」

データ分析をやっていると、こんな壁にぶつかることがあります。これが多重共線性という問題です。

実はこの厄介な問題、前回までに学んできた**PCA(主成分分析)**で解決できるんです。

最終回の今回は、PCAを回帰分析に応用する主成分回帰(PCR: Principal Component Regression)という手法を、Pythonで手を動かしながら学んでいきましょう!



📦 ライブラリの準備

今回は回帰分析も行うので、いつものPCA関連に加えてscikit-learnの回帰モジュールも読み込みます。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
import seaborn as sns

from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_squared_error, r2_score

# VIF計算用
from statsmodels.stats.outliers_influence import variance_inflation_factor

🏠 多重共線性って何? ~マンション価格予測で考える~

身近な例で理解する

あなたが不動産会社のデータサイエンティストだとします。マンションの価格を予測するモデルを作りたい。

説明変数として「延床面積」「部屋数」「築年数」「駅からの距離」などを使おうと考えます。

ここで問題が起きます。

「延床面積」が大きいマンションは「部屋数」も多いですよね。つまり、この2つの変数は強く相関しています。

このように、説明変数同士が強く相関している状態多重共線性(multicollinearity)といいます。

多重共線性があると、こんな困ったことが起きます。

📉 回帰係数がフラフラ不安定になる(ちょっとデータが変わるだけで係数が激変)
📉 係数の標準誤差が大きくなる(統計的検定の信頼性ダウン)
📉 係数の符号が常識と逆になる(「面積が大きいほど価格が下がる」とか…)

なぜこんなことが起きるのか?(数式で確認)

重回帰分析では、目的変数 $${\mathbf{y}}$$ を説明変数の行列 $${X}$$ を使って次のように表します。

$${\mathbf{y} = X\boldsymbol{\beta} + \boldsymbol{\varepsilon}}$$

最小二乗法による回帰係数の推定値は:

$${\hat{\boldsymbol{\beta}} = (X^\top X)^{-1} X^\top \mathbf{y}}$$

ポイントはこの式の中の $${X^\top X}$$ です。

説明変数同士が強く相関していると、$${X}$$ の列ベクトルが似た方向を向いてしまいます。すると $${X^\top X}$$ の行列式がほぼ0(=ほぼ特異行列)になり、逆行列の計算が不安定になるんです。

これが、係数がフラフラする根本原因です。


🔧 多重共線性のあるデータを作って体験してみよう

実際にデータを作って問題を体感してみましょう。

2つの潜在変数から、互いに相関の高い5つの説明変数を生成します。

以下、コードです。

# サンプル数
n = 100

# 基礎となる潜在変数(2つ)
z1 = np.random.normal(0, 1, n)
z2 = np.random.normal(0, 1, n)

# 5つの説明変数(互いに相関がある)
x1 = z1 + np.random.normal(0, 0.1, n)
x2 = z1 + np.random.normal(0, 0.15, n)  # x1と強く相関
x3 = z1 + np.random.normal(0, 0.2, n)   # x1, x2と強く相関
x4 = z2 + np.random.normal(0, 0.1, n)
x5 = z2 + np.random.normal(0, 0.15, n)  # x4と強く相関

# 説明変数をまとめる
X = np.column_stack([x1, x2, x3, x4, x5])
変数名リスト = ['x1', 'x2', 'x3', 'x4', 'x5']

# 目的変数(真の関係: y = 2*z1 + 3*z2 + noise)
y = 2 * z1 + 3 * z2 + np.random.normal(0, 0.5, n)

print(f"説明変数の形状: {X.shape}")
print(f"目的変数の形状: {y.shape}")

以下、実行結果です。

説明変数の形状: (100, 5)
目的変数の形状: (100,)

x1・x2・x3が互いに強く相関し、x4・x5も互いに強く相関するように設計しています。

相関行列をヒートマップで確認してみましょう。

以下、コードです。

# 相関行列を計算
df_X = pd.DataFrame(X, columns=変数名リスト)
相関行列 = df_X.corr()

# ヒートマップで可視化
plt.figure(figsize=(8, 6))
sns.heatmap(
    相関行列, 
    annot=True, fmt='.2f', cmap='coolwarm',
    center=0, linewidths=1, 
    vmin=-1, vmax=1
)
plt.title('説明変数間の相関行列')
plt.show()

x1-x2-x3のグループとx4-x5のグループ内で相関が0.9以上!これが多重共線性の典型パターンです。


📊 VIF(分散膨張係数)で多重共線性を検出する

多重共線性の程度を数値で測る指標が**VIF(Variance Inflation Factor)**です。

VIFの目安:
1 → 多重共線性なし
🟡 1〜5 → 軽度
🟠 5〜10 → 中程度(要注意)
🔴 10以上 → 深刻(対処が必要)

以下、コードです。

def calculate_vif(X, 変数名リスト):
    vif_data = pd.DataFrame()
    vif_data['変数'] = 変数名リスト
    vif_data['VIF'] = [variance_inflation_factor(X, i) for i in range(X.shape[1])]
    return vif_data

vif_result = calculate_vif(X, 変数名リスト)
print("=== VIF(分散膨張係数) ===")
print(vif_result.round(2))

以下、実行結果です。

=== VIF(分散膨張係数) ===
   変数    VIF
0  x1  47.38
1  x2  32.68
2  x3  20.95
3  x4  29.48
4  x5  29.44

全変数のVIFが10を大きく超えている!多重共線性が深刻であることがはっきりわかります。


⚠️ 通常の重回帰分析を試してみると……

多重共線性のあるデータに通常の重回帰分析を適用するとどうなるか見てみましょう。

以下、コードです。

# データを標準化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# 通常の重回帰分析
ols = LinearRegression()
ols.fit(X_scaled, y)

# 結果を表示
print("=== 通常の重回帰分析(OLS)の結果 ===")
係数df = pd.DataFrame({
    '変数': 変数名リスト,
    '回帰係数': ols.coef_
})
print(係数df.round(3))
print(f"\nR²スコア: {r2_score(y, ols.predict(X_scaled)):.3f}")

以下、実行結果です。

   変数   回帰係数
0  x1  1.300
1  x2  0.082
2  x3  0.526
3  x4  1.575
4  x5  1.430

R²スコア: 0.971

R²は高いけれど、x1・x2・x3は同じ潜在変数から生成されたのに係数がバラバラ。x2の係数がほぼ0になっているのも不自然です。

ブートストラップで係数の不安定さを確認

データをちょっと変えるだけで係数がどれだけ変わるか確認してみましょう。

以下、コードです。

係数リスト = []

for i in range(100):
    idx = np.random.choice(n, n, replace=True)
    X_boot = X_scaled[idx]
    y_boot = y[idx]
    
    model = LinearRegression()
    model.fit(X_boot, y_boot)
    係数リスト.append(model.coef_)

係数配列 = np.array(係数リスト)

係数安定性 = pd.DataFrame({
    '変数': 変数名リスト,
    '係数の平均': 係数配列.mean(axis=0),
    '係数の標準偏差': 係数配列.std(axis=0)
})
print("=== ブートストラップによる係数の安定性評価 ===")
print(係数安定性.round(3))

以下、実行結果です。

   変数  係数の平均  係数の標準偏差
0  x1  1.315    0.419
1  x2  0.081    0.365
2  x3  0.525    0.291
3  x4  1.616    0.351
4  x5  1.390    0.340

x2は標準偏差が平均を大きく上回っている状態。つまり、データが少し変わるだけで係数が正にも負にもなってしまう。これは信頼できる結果とは言えませんね。


🎯 主成分回帰(PCR)で解決しよう!

ここで登場するのが**主成分回帰(PCR)**です。

アイデアはとてもシンプル。

① 説明変数にPCAを適用して主成分を抽出する ② 元の説明変数の代わりに主成分を使って回帰を行う

なぜこれで多重共線性が解決するの?

PCAで得られる主成分は互いに**直交(=無相関)**しているからです。無相関な変数の間に多重共線性は存在しません。

📐 PCRの数式

標準化された説明変数の行列 $${X}$$ にPCAを適用すると、主成分得点 $${Z}$$ が得られます。

$${Z = XV}$$

$${V}$$ は固有ベクトル行列です。上位 $${k}$$ 個の主成分だけを取り出した $${Z_k}$$ を使って、PCRでは次の回帰モデルを推定します。

$${\mathbf{y} = Z_k \boldsymbol{\gamma} + \boldsymbol{\varepsilon}}$$

回帰係数の推定値:

$${\hat{\boldsymbol{\gamma}} = (Z_k^\top Z_k)^{-1} Z_k^\top \mathbf{y}}$$

ここで大事なのは、$${Z_k^\top Z_k}$$ が対角行列になること。対角行列の逆行列は常に安定して計算できるので、多重共線性の問題が解消されるんです。


🛠️ PCRをステップバイステップで実装する

Step 1:PCAで主成分を抽出

以下、コードです。

pca = PCA()
Z = pca.fit_transform(X_scaled)

print("=== PCAの結果 ===")
for i, ratio in enumerate(pca.explained_variance_ratio_):
    累積 = sum(pca.explained_variance_ratio_[:i+1])
    print(f"PC{i+1}: 寄与率 {ratio*100:.1f}%, 累積 {累積*100:.1f}%")

以下、実行結果です。

=== PCAの結果 ===
PC1: 寄与率 59.0%, 累積 59.0%
PC2: 寄与率 39.6%, 累積 98.6%
PC3: 寄与率 0.8%, 累積 99.4%
PC4: 寄与率 0.3%, 累積 99.7%
PC5: 寄与率 0.3%, 累積 100.0%

2つの潜在変数から生成されたデータなので、PC1とPC2だけで98.6%を説明できています。2成分で十分ですね。

Step 2:主成分で回帰を行う

以下、コードです。

成分数 = 2
Z_reduced = Z[:, :成分数]

pcr = LinearRegression()
pcr.fit(Z_reduced, y)

print(f"=== 主成分回帰(PCR, {成分数}成分)の結果 ===")
print("\n主成分空間での回帰係数:")
for i, coef in enumerate(pcr.coef_):
    print(f"  PC{i+1}: {coef:.3f}")
print(f"\nR²スコア: {r2_score(y, pcr.predict(Z_reduced)):.3f}")

以下、実行結果です。

=== 主成分回帰(PCR, 2成分)の結果 ===

主成分空間での回帰係数:
  PC1: 1.188
  PC2: 2.091

R²スコア: 0.970

主成分は互いに無相関なので、係数の推定が安定しています。R²もOLSとほぼ同等の0.970!

Step 3:元の変数空間の係数に変換する

主成分空間の係数 $${\hat{\boldsymbol{\gamma}}}$$ を、元の変数空間の係数 $${\hat{\boldsymbol{\beta}}}$$ に戻すには:

$${\hat{\boldsymbol{\beta}} = V_k \hat{\boldsymbol{\gamma}}}$$

固有ベクトル行列をかけるだけです。

以下、コードです。

V = pca.components_[:成分数].T 
beta_original = V @ pcr.coef_

print("=== 元の変数空間での係数(PCR vs OLS)===")
比較df = pd.DataFrame({
    '変数': 変数名リスト,
    'PCR係数': beta_original,
    'OLS係数': ols.coef_
})
print(比較df.round(3))

以下、実行結果です。

   変数  PCR係数  OLS係数
0  x1  0.664  1.300
1  x2  0.623  0.082
2  x3  0.627  0.526
3  x4  1.523  1.575
4  x5  1.498  1.430

PCRの係数は、x1・x2・x3がほぼ同じ値(0.62〜0.66)で安定しています。OLSでは0.08〜1.3とバラバラだったのと大違いですね。x4・x5もほぼ同じ値で、こちらも安定しています。


⚡ Pipelineでスッキリ実装

実務では、標準化→PCA→回帰をscikit-learnのPipelineにまとめると便利です。

以下、コードです。

pcr_pipeline = Pipeline([
    ('scaler', StandardScaler()),
    ('pca', PCA(n_components=2)),
    ('regressor', LinearRegression())
])

pcr_pipeline.fit(X, y)
y_pred = pcr_pipeline.predict(X)

print("=== Pipeline PCRの結果 ===")
print(f"R²スコア: {r2_score(y, y_pred):.3f}")

以下、実行結果です。

=== Pipeline PCRの結果 ===
R²スコア: 0.970

3行で前処理から予測まで完結!交差検証との組み合わせも簡単です。


🚀 発展編:PLS回帰との違い

PCRの発展形としてPLS回帰(Partial Least Squares Regression)があります。

両者の違いをひと言でいうと:

PCR(PCA) → $${X}$$ の分散を最大化する方向で成分を抽出。目的変数 $${y}$$ は考慮しない。
PLS → $${X}$$ と $${y}$$ の共分散を最大化する方向で成分を抽出。目的変数 $${y}$$ を考慮する。

数式で比較すると:

$${\text{PCA: } \max_{\mathbf{w}} \mathbf{w}^\top X^\top X \mathbf{w}}$$

$${\text{PLS: } \max_{\mathbf{w}} \mathbf{w}^\top X^\top \mathbf{y} \mathbf{y}^\top X \mathbf{w}}$$

PCAは「$${X}$$ のばらつきが大きい方向」を見つけますが、PLSは「$${X}$$ と $${y}$$ が一緒に変動する方向」を見つけます。

つまり、PLSは予測に直接役立つ成分を優先的に抽出するので、PCRより少ない成分数で高い精度を出せることが多いです。

以下、コードです。

from sklearn.cross_decomposition import PLSRegression

pls = PLSRegression(n_components=2)
pls.fit(X_scaled, y)

print("=== PCR vs PLS(2成分での比較) ===")
print(f"PCR R²: {r2_score(y, pcr.predict(Z_reduced)):.3f}")
print(f"PLS R²: {r2_score(y, pls.predict(X_scaled)):.3f}")

以下、実行結果です。

=== PCR vs PLS(2成分での比較) ===
PCR R²: 0.970
PLS R²: 0.970

今回のデータでは同じ精度ですが、実務データではPLSが有利になるケースがよくあります。


📚 参考にした情報源


🔜 次回予告

📝 まとめ

多重共線性 → 説明変数同士が強く相関している状態。回帰係数が不安定に。VIF≧10なら要対処。

主成分回帰(PCR) → PCAで無相関な主成分を抽出し、それを使って回帰。直交なので多重共線性が解消される。

Pipeline → 標準化→PCA→回帰を一括管理。実務ではこの形で使うのが便利。

PLS回帰 → PCRの発展形。目的変数との共分散を最大化する方向を抽出。予測が目的なら有利になることが多い。


🏁 シリーズ全体の振り返り

全6回にわたるシリーズ「Pythonで体感する次元削減入門」、これにて完走です!

第1回: なぜ「次元を減らす」のか? 情報の海で溺れないために
第2回: PCAのしくみを分散・共分散・固有値分解で理解する
第3回: 次元はいくつ残す? 寄与率とスクリープロットの読み方
第4回: scikit-learnでさくっとPCAを実施する
第5回: PCAと因子分析、何が違う?
第6回: 次元削減を回帰に活かす主成分回帰(PCR)← 今回

第1回で「データが広がる方向を見つける」というPCAの直感を学び、第2回で共分散行列と固有値分解という数学的基盤を理解しました。第3回で寄与率による次元数の決め方、第4回でscikit-learnでの実装、第5回でPCAと因子分析の違い、そして今回の第6回でPCAの回帰への応用を学びました。

次元削減は、データサイエンスの基礎でありながら奥が深いテーマです。ぜひ実際のデータ分析で活用してみてください!


最後まで読んでいただきありがとうございました! 「スキ」していただけると励みになります 🙌

いいなと思ったら応援しよう!