【徹底解説】多重共線性とは

本記事は機械学習の徹底解説シリーズに含まれます。

初学者の分かりやすさを優先するため，多少正確でない表現が混在することがあります。もし致命的な間違いがあればご指摘いただけると助かります。

はじめに

多重共線性は，重回帰分析における概念です。そこで，本記事では最初に回帰分析について復習してから，再び多重共線性に戻ってくるという構成をとりたいと思います。

回帰式

回帰式とは，説明変数$x$を用いて表される目的変数$y$との関係$y=f(x)$のことを指します。$y$は連続変数であり，特に$x$が一次元の場合は単回帰式，$y$が二次元以上の場合は重回帰式とよばれます。単回帰分析は，単回帰式を用いて分析を行う手法のことを指します。

回帰分析を実現する最も基本的な手法に，線形回帰があります。以下，行列を用いて定式化を行います。線形回帰では，目的変数の予測値$\hat{\mY}\in\bbR^{n}$を$\vb\in\bbR^{d}$で重み付けした$\mX\in\bbR^{n\times d}$の線形和で表します。

しかし，回帰式に基づく予測値$\hat{\mY}$と，実際の観測値$\mY\in\bbR^{n}$には誤差$\boldsymbol{\varepsilon}\in\bbR^{n}$が生じてしまいます。

上式を各要素ごとに書き下せば，以下のように表されます。

線形回帰では，この誤差$\varepsilon$の二乗和$Q$を最小化するパラメータを推定します。これは最小二乗法と呼ばれます。

二乗和$Q$は凸な二次関数であるため，$Q$の微分を$0$と置くことで各パラメータの推定値を求めることができます。ただし，後の結果が綺麗な形になるように，偏微分の係数を$1/2$とすることに注意して下さい。

最終的に得られた条件式は，正規方程式と呼ばれています。正規方程式は，$\mX^T \mX$が正則な場合は以下のような解を与えます。

この解は偏回帰係数と呼ばれています。特に，観測データ$\mX$を標準化した場合に得られる偏回帰係数は標準偏回帰係数と呼ばれ，単回帰（$d=1$）の場合は相関係数と等しくなります。このことを証明していきましょう。

まず，単回帰の場合の正規方程式を考えます。計算しやすいように，$\vb$について整理する前の正規方程式を採用します。ただし，$\sum x_i/N=\bar{x}$，$\vx = (x_1, \ldots, x_n)$，$\vy = (y_1, \ldots, y_n)$とします。

具体的に計算を進めていきます。まずは$b_{0}$について整理します。

続いて，$b_{1}$について整理します。

ただし，$S_{\alpha \beta}$は以下のように定義しました。

ここで，相関係数の定義を確認しましょう。

さて，観測データ$\vx$および$\vy$が平均$0$，分散$1$に標準化されているとすれば，$S_{xx}=S_{yy}$となります。したがって，$b_1$は以下のように計算されます。

したがって，観測データが標準化された場合に限り，単回帰の偏回帰係数$b_1$は相関係数$r$と等しくなることが示されました。以上の議論から，偏回帰係数は説明変数の目的変数に対する影響力の大きさを表すことが分かります。ただし，重回帰の場合は偏回帰係数の解釈には注意が必要であり，後述する多重共線性が生じるために注意が必要です。

補足

最小二乗法は，誤差$\boldsymbol{\varepsilon}$の分布に平均$0$かつ等分散であるガウス分布を仮定した場合の最尤推定に相当します。最尤推定の目的関数は，以下の対数尤度です。

最尤推定では$-Q$の最大化を考えるため，結局$Q$の最小化と等価になり，誤差の二乗和を最小化する問題に相当することが示されました。

多重共線性

ここまでの議論は$\mX^T \mX$が正則な場合を考えました。しかし，$\mX^T \mX$が正則でない場合には多重回帰には問題が生じます。$\mX^T \mX$が正則でないとき，$\mX^T \mX$はランク落ちの可能性があります。これはつまり，説明変数$\mX$がランク落ちであるという状況（説明変数同士に相関がある/説明変数が線形従属である）に相当します。このとき，無理矢理$\mX^T \mX$の逆行列を計算しようとすると，本来逆行列を持たないはずの値を計算しているために，$\vb$の値が不安定に（大きく）なってしまいます。この問題を多重共線性と呼びます。

多重共線性を解決する手法の一つに，正則化付き最小二乗法があります。具体的には，以下の目的関数$\mathcal{L}$を最小化します。

特に，$q=0$のときは素朴な最小二乗法（重回帰），$q=1$のときはLASSO，$q=2$のときはリッジ回帰と呼ばれます。正則化を施すことにより，回帰係数が不安定になる（大きくなりすぎる）ことを防ぐことができます。リッジ回帰の偏回帰係数は，素朴な最小二乗法のときと同様にして，以下のように求められます。

多重共線性を解決する他の手法としては，MP逆（一般逆行列）の利用も挙げられます。MPはMoore-Penroseの頭文字をとった略称です。Penrose（ペンローズ）がノーベル物理学賞を受賞したのは記憶に新しいですね。

リッジ回帰

ここからは，リッジ回帰を素朴な重回帰と比べることで，正則化の果たす役割を確認していきましょう。誤差関数に正則化項を付けたリッジ回帰では，偏回帰係数の大きさが抑制されることを具体的にみていきます。

データ作成

分析に用いるデータは，多重共線性を引き起こすデータが好ましいため，線形従属となるような説明変数・目的変数のペアを用意します。なお，今回は標準化後の観測データを作成するものとします。作成手順は以下の通りです。

このようにして得られた($\vx, y$)，($\vx^{\prime}, y^{\prime}$)には，多重共線性が生じると考えられます。そこで，リッジ回帰による正則化の働きを調べます。

パラメータ設定

本節で用いるパラメータは，以下のように設定します。

このパラメータに従って，$\vx$，$\vx^{\prime}$，$y$，$y^{\prime}$を作成します。

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
a = 2
c = 1/2

# xはN(0, 1)からサンプリングする
x = np.random.normal(0, 1, 10)

# xに線形従属であるx'を作成する
x_dash = 2*x + np.random.normal(0, 1, 10)

# y = 2x + epsilonに従って説明変数を作成する
y = sum(a*x) + np.random.normal(0, 1, 1)
y_dash = sum(a*x_dash) + np.random.normal(0, 1, 1)

# bには定数項も含めるためにxには1をappendしておく
x = np.insert(x, 0, 1)
x_dash = np.insert(x_dash, 0, 1)

# 2つの線形従属なデータ組みを行列にまとめる
X = np.array([x, x_dash])
Y = np.array([y, y_dash])

$\vx$は以下のようになりました。

$\vx^{\prime}$は以下のようになりました。

$y,y^{\prime}$は以下のようになりました。

以下，$\vb$を求めていきます。なぜなら，$\vb$に定数項を含めることで，$y = \vx^T\vb $，$y^{\prime} = (\vx^{\prime})^T \vb$と予測値$y$，$y^{\prime}$を求めることができるからです。

実験結果

$\lambda=0$の場合，つまり素朴な最小二乗法の場合は，変回帰係数の値と絶対値の平均は以下のようになります。

# 最小二乗法における正規方程式の解
b = np.linalg.inv(X.T @ X) @ X.T @ Y

$\vb$は以下のようになりました。

$|\vb|$の期待値は以下のように計算されます。

$y,y^{\prime}$は以下のようになりました。

説明変数間が線形従属であるため，$\mX^T \mX$がランク落ちして逆行列が計算できないような状況であるため，正しく$y$と$y^{\prime}$が推定できていないことが分かります。$\lambda=1$の場合，つまりリッジ回帰の場合は，変回帰係数の値と絶対値の平均は以下のようになります。

# リッジ回帰における正規方程式の解
lmd = 1
b = np.linalg.inv(X.T @ X + lmd * np.eye(X.shape[1])) @ X.T @ Y

$\vb$は以下のようになりました。

$|\vb|$の期待値は以下のように計算されます。

$y,y^{\prime}$は以下のようになりました。

この結果から，素朴な最小二乗法と比べて偏回帰係数の大きさがしっかりと抑制されていることが分かります。さらに，$y$と$y^{\prime}$もよく推定できていることも分かります。一方，$\lambda$を大きくすると，より正則化の効果が強まります。例えば，$\lambda=10^3$の場合は，変回帰係数の値と絶対値の平均は以下のようになります。

# リッジ回帰における正規方程式の解
lmd = 1e3
b = np.linalg.inv(X.T @ X + lmd * np.eye(X.shape[1])) @ X.T @ Y

$\vb$は以下のようになりました。

$|\vb|$の期待値は以下のように計算されます。

$y,y^{\prime}$は以下のようになりました。

これでは，回帰分析によって予測したい$\vy$に対してモデルを当てはめられたとは言えません。$y$と$y^{\prime}$もうまく推定できていません。そのため，$\lambda$には偏回帰係数の大きさを適切に抑制するのに妥当な値を利用する必要があります。