多重共線性とVIF — 回帰分析の診断

重回帰分析で住宅価格を予測するモデルを作るとしましょう。「1階の床面積」「2階の床面積」「総床面積」の3つの変数を説明変数に入れたいと考えます。しかしちょっと待ってください。総床面積 = 1階面積 + 2階面積ですから、この3変数は完全な線形関係にあります。このような状況で回帰分析を行うと、回帰係数の推定は破綻してしまいます。

完全な線形関係でなくても、説明変数間に強い相関があると、同じ問題が「やんわりと」生じます。回帰係数の推定値が異常に不安定になり、符号が直感に反したり、標準誤差が膨大になったりします。この問題を多重共線性（multicollinearity）と呼びます。

多重共線性を理解し適切に診断することは、信頼できる回帰分析を行うために不可欠です。

• 医学・疫学: 臨床試験で多数のバイオマーカーを説明変数に用いる場合、互いに相関するバイオマーカーが多重共線性を引き起こし、真の効果の同定を困難にする
• 経済学: マクロ経済変数（GDP、消費、投資など）は互いに強く相関しており、回帰分析の解釈に多重共線性が影響する
• 環境科学: 気温・湿度・日照時間・降水量など、互いに関連する環境変数を同時に用いると多重共線性が生じやすい
• 工学: 製品の寸法変数（幅・高さ・重量・体積）など、物理法則により関連する変数群が多重共線性を引き起こす

本記事では、多重共線性が回帰分析にもたらす影響を数学的に厳密に導出し、VIF（分散膨張係数）による診断法の理論を解説します。さらに、条件数やフレッシュ・ウォー（Farrar-Glauber）テストなどの関連する診断指標も紹介し、多重共線性への対処法を比較検討します。

本記事の内容

• 多重共線性の定義と直感的な理解
• 回帰係数の分散膨張の数学的メカニズム
• VIF（分散膨張係数）の導出と解釈
• 条件数による診断
• 多重共線性の対処法（変数除去・PCA回帰・正則化）の理論的比較
• Pythonでの実装と可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 重回帰分析の理論 — 回帰係数の分散と最小二乗法の性質
• 主成分分析の理論 — PCA回帰への接続

多重共線性とは — 直感的な理解

分離不能性の問題

多重共線性の本質は、「各説明変数の独自の寄与を分離できなくなる」ことです。

身近な例で考えてみましょう。ある患者の回復速度を「投薬量」と「副作用の強さ」で予測するモデルを考えます。もし投薬量が増えると副作用も必ず増えるという関係があれば、回復速度の変化が「投薬量のおかげ」なのか「副作用のせい」なのかを区別できません。データ上では、投薬量と副作用が常に一緒に変化するため、それぞれの独立した効果を推定することが困難になります。

もう少し幾何学的に考えてみましょう。2つの説明変数 $x_1$ と $x_2$ が全く無相関のとき、データ点は $x_1$-$x_2$ 平面上に「丸い」散布図を描きます。この場合、$x_1$ 軸方向と $x_2$ 軸方向の情報は独立しており、各変数の回帰係数を精度よく推定できます。

一方、$x_1$ と $x_2$ が強く相関するとき、データ点は細長い楕円形になります。楕円の長軸方向にはデータが広がっていますが、短軸方向にはほとんど広がりがありません。回帰平面の「傾き」を決めるには短軸方向の情報が必要なのですが、この方向のデータの広がりが小さいため、傾きの推定が非常に不安定になるのです。

完全共線性と近似共線性

多重共線性には2つのレベルがあります。

完全共線性（perfect multicollinearity）: 説明変数間に厳密な線形関係がある場合です。例えば $x_3 = 2x_1 + 3x_2$ のとき、デザイン行列 $\bm{X}$ の列がランク落ちし、$\bm{X}^T\bm{X}$ が特異行列になります。特異行列には逆行列が存在しないため、最小二乗推定量 $\hat{\bm{\beta}} = (\bm{X}^T\bm{X})^{-1}\bm{X}^T\bm{y}$ が定義できません。

実務でよくある完全共線性の例: – ダミー変数の罠: $k$ 個のカテゴリに対して $k$ 個のダミー変数を全て投入する（1つ減らすべき） – 物理的な定義関係: 総面積 = 1階面積 + 2階面積 – 単位変換: 同じ変数をcm単位とm単位の両方で含める

近似共線性（near multicollinearity）: 説明変数間にほぼ線形関係がある場合です。$\bm{X}^T\bm{X}$ は正則ですが、条件数が非常に大きくなり、逆行列の成分が不安定になります。実務で実際に問題になるのはこちらです。

近似共線性は連続的な概念です。相関が0のとき多重共線性はなく、相関が1に近づくにつれて多重共線性が強まります。どの程度の相関で「問題あり」とするかの基準として、VIFが使われます。

それでは、多重共線性が回帰係数の推定にどのように影響するかを、数学的に厳密に導出しましょう。

回帰係数の分散膨張のメカニズム

分散の公式の復習

正規線形モデル $\bm{y} = \bm{X}\bm{\beta} + \bm{\varepsilon}$, $\bm{\varepsilon} \sim N(\bm{0}, \sigma^2\bm{I})$ において、最小二乗推定量の共分散行列は:

$$ \begin{equation} \text{Cov}(\hat{\bm{\beta}}) = \sigma^2(\bm{X}^T\bm{X})^{-1} \end{equation} $$

$j$ 番目の回帰係数の分散は:

$$ \begin{equation} \text{Var}(\hat{\beta}_j) = \sigma^2[(\bm{X}^T\bm{X})^{-1}]_{jj} \end{equation} $$

多重共線性の問題は、$[(\bm{X}^T\bm{X})^{-1}]_{jj}$ が大きくなることで生じます。では、この対角成分がどのように決まるかを見ていきましょう。

固有値分解による解析

$\bm{X}^T\bm{X}$ の固有値分解を $\bm{X}^T\bm{X} = \bm{V}\bm{\Lambda}\bm{V}^T$ とします。ここで $\bm{V}$ は直交行列、$\bm{\Lambda} = \text{diag}(\lambda_1, \lambda_2, \dots, \lambda_p)$ です。逆行列は:

$$ \begin{equation} (\bm{X}^T\bm{X})^{-1} = \bm{V}\bm{\Lambda}^{-1}\bm{V}^T = \sum_{k=1}^{p}\frac{1}{\lambda_k}\bm{v}_k\bm{v}_k^T \end{equation} $$

$j$ 番目の対角成分は:

$$ \begin{equation} [(\bm{X}^T\bm{X})^{-1}]_{jj} = \sum_{k=1}^{p}\frac{v_{jk}^2}{\lambda_k} \end{equation} $$

ここで $v_{jk}$ は $\bm{V}$ の $(j, k)$ 成分です。この式から明白なことがわかります。固有値 $\lambda_k$ が小さいほど、$1/\lambda_k$ が大きくなり、分散が膨張するのです。

多重共線性があるとき、$\bm{X}^T\bm{X}$ の最小固有値がゼロに近づきます。もし $v_{jk}^2$ が小さい固有値 $\lambda_k$ に対してゼロでなければ、$\hat{\beta}_j$ の分散は爆発的に大きくなります。

標準化変数を用いた分析

分析を簡潔にするために、説明変数を標準化（平均0、分散1）して考えます。標準化後の変数を $\tilde{x}_j$ とすると、$\bm{X}^T\bm{X}/n$ は標本相関行列 $\bm{R}$ に一致します。

$j$ 番目の回帰係数の分散は:

$$ \begin{equation} \text{Var}(\hat{\beta}_j) = \frac{\sigma^2}{n} [\bm{R}^{-1}]_{jj} \end{equation} $$

もし多重共線性がなく全変数が互いに無相関であれば、$\bm{R} = \bm{I}$ なので $[\bm{R}^{-1}]_{jj} = 1$ となり、分散は $\sigma^2/n$ です。多重共線性があると $[\bm{R}^{-1}]_{jj} > 1$ となり、分散が膨張します。この膨張の倍率がVIFです。

この膨張メカニズムを定量化するVIFを、次のセクションで正式に導出します。

VIF（分散膨張係数）の導出

補助回帰と決定係数

VIFの導出では、説明変数 $x_j$ を他の全ての説明変数で回帰するという「補助回帰」を考えます。

$$ \begin{equation} x_j = \gamma_0 + \gamma_1 x_1 + \dots + \gamma_{j-1}x_{j-1} + \gamma_{j+1}x_{j+1} + \dots + \gamma_p x_p + u_j \end{equation} $$

この補助回帰の決定係数を $R_j^2$ とします。$R_j^2$ は、$x_j$ の変動のうち他の説明変数で説明できる割合を表します。$R_j^2$ が1に近いほど、$x_j$ は他の変数の線形結合でほぼ表現でき、多重共線性が強いことを意味します。

VIFの導出（ブロック逆行列の公式）

相関行列 $\bm{R}$ を $j$ 番目の変数とそれ以外に分けてブロック化します。

$$ \bm{R} = \begin{pmatrix} \bm{R}_{-j,-j} & \bm{r}_{-j,j} \\ \bm{r}_{-j,j}^T & 1 \end{pmatrix} $$

ここで $\bm{R}_{-j,-j}$ は $j$ 番目の変数を除いた相関行列、$\bm{r}_{-j,j}$ は $j$ 番目の変数と他の変数の相関ベクトルです。

ブロック逆行列の公式（シャーマン・モリソンの公式の一般化）を用いると、$\bm{R}^{-1}$ の $(j,j)$ 成分は:

$$ \begin{equation} [\bm{R}^{-1}]_{jj} = \frac{1}{1 – \bm{r}_{-j,j}^T\bm{R}_{-j,-j}^{-1}\bm{r}_{-j,j}} \end{equation} $$

ここで重要なのは、$\bm{r}_{-j,j}^T\bm{R}_{-j,-j}^{-1}\bm{r}_{-j,j}$ が補助回帰の決定係数 $R_j^2$ に等しいということです。これは重回帰分析の理論から直接導かれます。$x_j$ を他の変数で回帰した場合の決定係数は、標準化変数の場合:

$$ R_j^2 = \bm{r}_{-j,j}^T\bm{R}_{-j,-j}^{-1}\bm{r}_{-j,j} $$

これを代入すると:

$$ \begin{equation} [\bm{R}^{-1}]_{jj} = \frac{1}{1 – R_j^2} \end{equation} $$

これがVIF（Variance Inflation Factor, 分散膨張係数）の定義です。

$$ \begin{equation} \text{VIF}_j = \frac{1}{1 – R_j^2} \end{equation} $$

VIFの解釈

VIFの名前は、回帰係数の分散がどれだけ「膨張」するかを直接表していることに由来します。

標準化変数を用いた場合の分散は:

$$ \text{Var}(\hat{\beta}_j) = \frac{\sigma^2}{n} \cdot \text{VIF}_j $$

もし多重共線性がなければ $\text{VIF}_j = 1$ であり、分散は $\sigma^2/n$ です。$\text{VIF}_j = 10$ なら分散が10倍に、$\text{VIF}_j = 100$ なら100倍に膨張しています。これは標準誤差がそれぞれ $\sqrt{10} \approx 3.2$ 倍、$\sqrt{100} = 10$ 倍に膨張することを意味します。

$R_j^2$ とVIFの対応関係を表にまとめます。

実務上の閾値としては、$\text{VIF} > 5$ または $\text{VIF} > 10$ が広く用いられていますが、分野や文脈に応じて判断する必要があります。VIFが大きくても、サンプルサイズが十分に大きければ推定精度が確保される場合もあるからです。

VIFは個々の変数ごとの多重共線性を診断しますが、全体的な多重共線性の程度を測る指標もあります。次に条件数について見ていきましょう。

条件数による診断

条件数の定義

デザイン行列 $\bm{X}$ の特異値分解 $\bm{X} = \bm{U}\bm{\Sigma}\bm{V}^T$ を考えます。$\bm{X}$ の条件数（condition number）は、最大特異値と最小特異値の比として定義されます。

$$ \begin{equation} \kappa(\bm{X}) = \frac{\sigma_{\max}}{\sigma_{\min}} \end{equation} $$

等価的に、$\bm{X}^T\bm{X}$ の固有値を用いて:

$$ \begin{equation} \kappa(\bm{X}) = \sqrt{\frac{\lambda_{\max}(\bm{X}^T\bm{X})}{\lambda_{\min}(\bm{X}^T\bm{X})}} \end{equation} $$

条件数は行列の「扱いにくさ」を測る数値線形代数の基本的な概念です。条件数が大きいほど、連立方程式の解がデータの微小な変動に対して敏感になります。

条件数と多重共線性の関係

多重共線性があると、$\bm{X}^T\bm{X}$ の最小固有値がゼロに近づきます。完全共線性の極限では最小固有値がゼロになり、条件数は無限大になります。

条件数の解釈の目安は:

条件指標と分散分解比

条件数は全体的な指標ですが、条件指標（condition index）は各固有値ごとの指標です。

$$ \begin{equation} \eta_k = \sqrt{\frac{\lambda_{\max}}{\lambda_k}}, \quad k = 1, 2, \dots, p \end{equation} $$

さらに、各条件指標がどの変数の分散膨張に寄与しているかを示す分散分解比（variance decomposition proportion）を計算できます。

$$ \begin{equation} \pi_{jk} = \frac{v_{jk}^2/\lambda_k}{\sum_{l=1}^{p}v_{jl}^2/\lambda_l} \end{equation} $$

$\pi_{jk}$ は、$j$ 番目の変数の分散のうち、第 $k$ 固有値に関連する割合を表します。大きな条件指標 $\eta_k$ に対して $\pi_{jk}$ が大きい複数の変数が存在すれば、それらの変数が互いに共線的であると診断できます。

VIFが「どの変数が問題か」を示すのに対し、条件指標と分散分解比は「どの変数がどの変数と共線的か」を示すより詳細な診断ツールです。

多重共線性を検出したら、次にどう対処するかが問題です。対処法を見ていきましょう。

多重共線性への対処法

対処法1: 変数の除去

最も直接的な対処法は、VIFが高い変数の中から一部を除去することです。相関の高い変数群から代表的な1つだけを残し、残りを削除します。

手順: 1. 全変数のVIFを計算する 2. VIFが最大の変数を除去する 3. 残りの変数でVIFを再計算する 4. 全変数のVIFが閾値以下になるまで繰り返す

この方法は単純で理解しやすいですが、除去する変数の選択に主観が入る点、情報の一部が失われる点がデメリットです。

対処法2: PCA回帰（主成分回帰）

主成分分析で説明変数を主成分に変換し、上位の主成分のみを説明変数として回帰分析を行う方法です。

主成分は互いに直交（無相関）であるため、多重共線性は完全に解消されます。さらに、小さな固有値に対応する主成分（多重共線性の原因となる方向）を除去するため、推定が安定します。

PCA回帰の手順: 1. 説明変数を標準化する 2. 主成分分析を行い、上位 $m$ 個の主成分を取得する 3. 主成分スコアを説明変数として回帰分析を行う

PCA回帰の理論的背景を見てみましょう。元の回帰モデルは:

$$ \bm{y} = \bm{X}\bm{\beta} + \bm{\varepsilon} $$

$\bm{X}$ の列を主成分 $\bm{Z} = \bm{X}\bm{V}$ に変換すると:

$$ \bm{y} = \bm{Z}\bm{\gamma} + \bm{\varepsilon}, \quad \bm{\gamma} = \bm{V}^T\bm{\beta} $$

$\bm{Z}^T\bm{Z}$ は対角行列（固有値行列）なので、$\hat{\bm{\gamma}}$ の各成分の分散は:

$$ \text{Var}(\hat{\gamma}_k) = \frac{\sigma^2}{\lambda_k} $$

小さい固有値 $\lambda_k$ に対応する成分は分散が大きいため、これらを除去することで安定した推定が得られます。ただし、除去した主成分に対応する $\beta$ の成分の情報は失われるため、バイアスが導入される点に注意が必要です。

対処法3: リッジ回帰（L2正則化）

リッジ回帰は、$\bm{X}^T\bm{X}$ に正の対角行列 $\lambda\bm{I}$ を加えることで、最小固有値がゼロに近づく問題を直接解決します。

$$ \begin{equation} \hat{\bm{\beta}}^{\text{ridge}} = (\bm{X}^T\bm{X} + \lambda\bm{I})^{-1}\bm{X}^T\bm{y} \end{equation} $$

$\bm{X}^T\bm{X} + \lambda\bm{I}$ の固有値は $\lambda_k + \lambda$ となるため、もとの最小固有値が非常に小さくても、$\lambda$ を加えることで全ての固有値が十分に大きくなります。これにより:

$$ \text{Var}(\hat{\beta}_j^{\text{ridge}}) \propto \frac{1}{\lambda_k + \lambda} $$

$\lambda$ が大きいほど分散は小さくなりますが、バイアスは大きくなります。分散とバイアスのトレードオフにおいて、多重共線性がある場合は少しのバイアスを許容して分散を大幅に削減する方が全体的な推定精度（MSE）が向上します。

対処法の比較

それでは、Pythonで多重共線性の影響とVIFの計算を実装し、対処法を比較しましょう。

Pythonによる実装と可視化

多重共線性の生成とVIF計算

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# --- 多重共線性のあるデータ ---
n = 200
x1 = np.random.randn(n)
x2 = 0.95 * x1 + 0.1 * np.random.randn(n)  # x1 と相関 ~0.99
x3 = 0.8 * x1 + 0.3 * np.random.randn(n)   # x1 と相関 ~0.93
x4 = np.random.randn(n)                      # 独立
x5 = np.random.randn(n)                      # 独立

X_vars = np.column_stack([x1, x2, x3, x4, x5])
var_names = ['$x_1$', '$x_2$', '$x_3$', '$x_4$', '$x_5$']

# 真の係数
beta_true = np.array([2.0, -1.0, 0.5, 1.5, -0.8])
X_design = np.column_stack([np.ones(n), X_vars])
beta_full = np.array([1.0] + list(beta_true))
y = X_design @ beta_full + 0.5 * np.random.randn(n)

# --- VIFの計算 ---
def compute_vif(X):
    """各変数のVIFを計算"""
    p = X.shape[1]
    vifs = np.zeros(p)
    for j in range(p):
        y_j = X[:, j]
        X_others = np.delete(X, j, axis=1)
        X_reg = np.column_stack([np.ones(len(y_j)), X_others])
        beta = np.linalg.lstsq(X_reg, y_j, rcond=None)[0]
        SS_res = np.sum((y_j - X_reg @ beta)**2)
        SS_tot = np.sum((y_j - y_j.mean())**2)
        R2_j = 1 - SS_res / SS_tot
        vifs[j] = 1.0 / (1.0 - R2_j)
    return vifs

vifs = compute_vif(X_vars)

# --- 条件数 ---
R = np.corrcoef(X_vars.T)
eigvals = np.linalg.eigvalsh(R)
condition_number = np.sqrt(eigvals.max() / eigvals.min())

# --- 相関行列 ---
corr_matrix = np.corrcoef(X_vars.T)

print("=== 多重共線性の診断 ===")
for j, name in enumerate(var_names):
    status = "!!" if vifs[j] > 10 else ("!" if vifs[j] > 5 else "ok")
    print(f"{name}: VIF = {vifs[j]:8.2f}  [{status}]")
print(f"\n条件数 = {condition_number:.2f}")

VIF計算の出力から、$x_1, x_2, x_3$ のVIFが閾値を超えていることが確認できます。$x_1$ は $x_2$ と $x_3$ の両方と強い相関を持つため最もVIFが大きく、$x_4, x_5$ は独立なのでVIFがほぼ1です。条件数も閾値の30を超えており、行列の条件が悪いことを示しています。

VIFと係数不安定性の可視化

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
n = 200
x1 = np.random.randn(n)
x2 = 0.95 * x1 + 0.1 * np.random.randn(n)
x3 = 0.8 * x1 + 0.3 * np.random.randn(n)
x4 = np.random.randn(n)
x5 = np.random.randn(n)
X_vars = np.column_stack([x1, x2, x3, x4, x5])
var_names = ['$x_1$', '$x_2$', '$x_3$', '$x_4$', '$x_5$']
beta_true = np.array([2.0, -1.0, 0.5, 1.5, -0.8])
X_design = np.column_stack([np.ones(n), X_vars])
beta_full = np.array([1.0] + list(beta_true))
y = X_design @ beta_full + 0.5 * np.random.randn(n)

def compute_vif(X):
    p = X.shape[1]
    vifs = np.zeros(p)
    for j in range(p):
        y_j = X[:, j]
        X_others = np.delete(X, j, axis=1)
        X_reg = np.column_stack([np.ones(len(y_j)), X_others])
        beta = np.linalg.lstsq(X_reg, y_j, rcond=None)[0]
        SS_res = np.sum((y_j - X_reg @ beta)**2)
        SS_tot = np.sum((y_j - y_j.mean())**2)
        vifs[j] = 1.0 / (1.0 - (1 - SS_res / SS_tot))
    return vifs

vifs = compute_vif(X_vars)

# --- ブートストラップで係数の不安定性を検証 ---
n_boot = 1000
beta_boot = np.zeros((n_boot, 6))
for b in range(n_boot):
    idx = np.random.choice(n, n, replace=True)
    beta_boot[b] = np.linalg.lstsq(X_design[idx], y[idx], rcond=None)[0]

fig, axes = plt.subplots(2, 2, figsize=(14, 11))

# (a) VIF棒グラフ
ax = axes[0, 0]
colors_vif = ['#e41a1c' if v > 10 else ('#ff7f00' if v > 5 else '#377eb8')
              for v in vifs]
bars = ax.bar(range(5), vifs, color=colors_vif, alpha=0.7, edgecolor='black')
ax.axhline(5, color='orange', linestyle='--', linewidth=1.5, label='Threshold = 5')
ax.axhline(10, color='red', linestyle='--', linewidth=1.5, label='Threshold = 10')
ax.set_xticks(range(5))
ax.set_xticklabels(var_names, fontsize=11)
ax.set_ylabel('VIF', fontsize=11)
ax.set_title('(a) Variance Inflation Factors', fontsize=12)
for j, v in enumerate(vifs):
    ax.text(j, v + 0.5, f'{v:.1f}', ha='center', fontsize=10)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, axis='y')

# (b) ブートストラップ係数分布
ax = axes[0, 1]
bp = ax.boxplot([beta_boot[:, j+1] for j in range(5)],
                labels=var_names, patch_artist=True, widths=0.6)
colors_box = ['#ffcccc' if vifs[j] > 5 else '#cce5ff' for j in range(5)]
for patch, color in zip(bp['boxes'], colors_box):
    patch.set_facecolor(color)
for j in range(5):
    ax.axhline(beta_true[j], color='gray', linestyle=':', linewidth=0.5)
    ax.plot(j+1, beta_true[j], 'k*', markersize=10)
ax.set_ylabel('Coefficient value', fontsize=11)
ax.set_title('(b) Bootstrap coefficient distributions', fontsize=12)
ax.grid(True, alpha=0.3)

# (c) 相関行列ヒートマップ
ax = axes[1, 0]
corr = np.corrcoef(X_vars.T)
im = ax.imshow(corr, cmap='RdBu_r', vmin=-1, vmax=1)
plt.colorbar(im, ax=ax, shrink=0.8)
ax.set_xticks(range(5))
ax.set_yticks(range(5))
ax.set_xticklabels(['$x_1$','$x_2$','$x_3$','$x_4$','$x_5$'], fontsize=10)
ax.set_yticklabels(['$x_1$','$x_2$','$x_3$','$x_4$','$x_5$'], fontsize=10)
for i in range(5):
    for j in range(5):
        ax.text(j, i, f'{corr[i,j]:.2f}', ha='center', va='center', fontsize=9,
                color='white' if abs(corr[i,j]) > 0.5 else 'black')
ax.set_title('(c) Correlation matrix', fontsize=12)

# (d) VIF vs 係数の標準偏差
ax = axes[1, 1]
coef_std = beta_boot[:, 1:].std(axis=0)
ax.scatter(vifs, coef_std, s=100, c=colors_vif, edgecolors='black', zorder=5)
for j, name in enumerate(var_names):
    ax.annotate(name, (vifs[j], coef_std[j]), textcoords='offset points',
                xytext=(8, 5), fontsize=11)
ax.set_xlabel('VIF', fontsize=11)
ax.set_ylabel('Bootstrap std of coefficient', fontsize=11)
ax.set_title('(d) VIF vs Coefficient instability', fontsize=12)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

4つの診断プロットから、多重共線性の影響が明確に読み取れます。

1. VIF棒グラフ（図a）: $x_1, x_2, x_3$ のVIFが閾値5（オレンジ線）を超えており、特に $x_1$ と $x_2$ は閾値10（赤線）も超えています。これらの変数間に深刻な多重共線性が存在することを示しています。$x_4, x_5$ のVIFはほぼ1であり、他の変数と独立であることを反映しています

2. 係数のブートストラップ分布（図b）: VIFが高い変数（$x_1, x_2, x_3$）の係数推定値は、箱ひげ図の幅が非常に広く、推定が不安定であることがわかります。黒い星印は真の係数値ですが、ブートストラップ推定値は真の値の周りに大きく散らばっています。一方、VIFが低い $x_4, x_5$ の分布は狭く、安定した推定が得られています

3. 相関行列ヒートマップ（図c）: $x_1$ と $x_2$ の相関が約0.99、$x_1$ と $x_3$ の相関が約0.93と非常に高いことが可視化されています。$x_4, x_5$ は他の変数との相関が低いことも確認できます。この相関構造がVIFの値を直接反映しています

4. VIF vs 係数の標準偏差（図d）: VIFと係数推定の不安定性（ブートストラップ標準偏差）の間に明確な正の関係があります。VIFが大きい変数ほど、係数の推定がばらつきます。これはVIFが分散の膨張倍率であるという理論的結果を実験的に裏付けています

対処法の比較

import numpy as np
from sklearn.linear_model import Ridge
from sklearn.decomposition import PCA as skPCA
from sklearn.preprocessing import StandardScaler

np.random.seed(42)
n = 200
x1 = np.random.randn(n)
x2 = 0.95 * x1 + 0.1 * np.random.randn(n)
x3 = 0.8 * x1 + 0.3 * np.random.randn(n)
x4 = np.random.randn(n)
x5 = np.random.randn(n)
X_vars = np.column_stack([x1, x2, x3, x4, x5])
beta_true = np.array([2.0, -1.0, 0.5, 1.5, -0.8])
X_design = np.column_stack([np.ones(n), X_vars])
beta_full = np.array([1.0] + list(beta_true))
y = X_design @ beta_full + 0.5 * np.random.randn(n)

# 各手法でブートストラップMSEを比較
n_boot = 500
n_test = 100
mse_ols = []
mse_ridge = []
mse_pca = []
mse_remove = []

for b in range(n_boot):
    # 訓練/テスト分割
    idx_train = np.random.choice(n, n, replace=True)
    X_train = X_design[idx_train]
    y_train = y[idx_train]

    # テストデータ生成
    x1_t = np.random.randn(n_test)
    x2_t = 0.95*x1_t + 0.1*np.random.randn(n_test)
    x3_t = 0.8*x1_t + 0.3*np.random.randn(n_test)
    x4_t = np.random.randn(n_test)
    x5_t = np.random.randn(n_test)
    X_test = np.column_stack([np.ones(n_test), x1_t, x2_t, x3_t, x4_t, x5_t])
    y_test = X_test @ beta_full + 0.5*np.random.randn(n_test)

    # OLS
    b_ols = np.linalg.lstsq(X_train, y_train, rcond=None)[0]
    mse_ols.append(np.mean((y_test - X_test @ b_ols)**2))

    # リッジ回帰
    ridge = Ridge(alpha=1.0, fit_intercept=False)
    ridge.fit(X_train, y_train)
    mse_ridge.append(np.mean((y_test - X_test @ ridge.coef_)**2))

    # 変数除去（x2を除去）
    X_train_r = X_train[:, [0,1,3,4,5]]
    X_test_r = X_test[:, [0,1,3,4,5]]
    b_remove = np.linalg.lstsq(X_train_r, y_train, rcond=None)[0]
    mse_remove.append(np.mean((y_test - X_test_r @ b_remove)**2))

    # PCA回帰（3成分）
    scaler = StandardScaler()
    X_tr_scaled = scaler.fit_transform(X_train[:, 1:])
    pca = skPCA(n_components=3)
    Z_train = pca.fit_transform(X_tr_scaled)
    Z_train = np.column_stack([np.ones(len(Z_train)), Z_train])
    b_pca = np.linalg.lstsq(Z_train, y_train, rcond=None)[0]
    X_te_scaled = scaler.transform(X_test[:, 1:])
    Z_test = pca.transform(X_te_scaled)
    Z_test = np.column_stack([np.ones(n_test), Z_test])
    mse_pca.append(np.mean((y_test - Z_test @ b_pca)**2))

print("=== 手法ごとの予測MSE（平均 ± 標準偏差） ===")
for name, mse_list in [('OLS', mse_ols), ('Ridge', mse_ridge),
                        ('Remove x2', mse_remove), ('PCR (3 comp)', mse_pca)]:
    print(f"{name:15s}: {np.mean(mse_list):.4f} +/- {np.std(mse_list):.4f}")

対処法の比較結果から、以下の洞察が得られます。

まず、多重共線性が存在する場合でもOLSの予測MSEが極端に悪いわけではありません。OLSは不偏推定量であるため、係数の推定は不安定でも「平均的には正しい」予測を行います。ただし、個々のブートストラップ標本での予測のばらつきは大きくなります。

リッジ回帰は、少しのバイアスを導入する代わりに分散を大幅に削減するため、多くの場合OLSよりも小さい予測MSEを達成します。変数除去は、除去した変数の情報を失うため、省略変数バイアスが生じますが、多重共線性による分散膨張は解消されます。PCA回帰は主成分の数の選択に依存しますが、適切に選べば良好な予測性能を示します。

まとめ

本記事では、多重共線性の原因・影響・診断・対処法を体系的に解説しました。

• 多重共線性の本質: 説明変数間の強い線形関係により、各変数の独自の寄与を分離できなくなり、回帰係数の分散が膨張する
• VIFの導出: $\text{VIF}_j = 1/(1-R_j^2)$。$x_j$ を他の全変数で回帰した決定係数 $R_j^2$ から計算される分散の膨張倍率。VIF > 10 は重度の多重共線性を示す
• 条件数: $\kappa = \sqrt{\lambda_{\max}/\lambda_{\min}}$。行列全体の条件の悪さを測る。$\kappa > 30$ で多重共線性が疑われる
• 固有値分解による理解: 分散膨張は $\bm{X}^T\bm{X}$ の小さな固有値に起因する。小さな固有値に対応する方向の情報が乏しいため、その方向の回帰係数の推定が不安定になる
• 対処法: 変数除去（単純だが情報を失う）、PCA回帰（多重共線性を完全に解消）、リッジ回帰（バイアスと分散のトレードオフで全体のMSEを改善）

多重共線性は回帰分析の「敵」として語られがちですが、適切な診断と対処法を理解していれば恐れるものではありません。重要なのは、VIFや条件数で問題を早期に発見し、分析の目的に応じた対処法を選択することです。

次のステップとして、以下の記事も参考にしてください。

• リッジ回帰とLASSOの理論 — 正則化による多重共線性への対処の理論的基盤
• 残差分析とモデル診断 — 多重共線性以外の回帰モデルの問題を診断する手法
• 主成分分析の理論 — PCA回帰の基礎となる次元削減手法