【コラム】多重共線性

URLをコピーする
URLをコピーしました!
zuka

こんにちは。
zuka(@beginaid)です。

この記事はコラムシリーズの1つになります。今回は大学の講義で多重共線性について学んだので,そちらの概要をメモ程度に記している内容になります。コラム記事のまとめは以下のページからご覧ください。なお,本記事の内容に関して正確性の保証はできません。

zuka

心優しき有識者の方々でもし間違いを見つけた方がいらっしゃれば,ご指摘いただけますと非常に助かります。

目次

多重共線性とは

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

zuka

多重共線性はできる限り起こってほしくない現象の一つだよ。

回帰式

回帰式とは,説明変数$x$を用いて表される目的変数$y$との関係$y=f(x)$のことを指します。$y$は連続変数であり,特に$x$が1次元の場合は単回帰式,$y$が2次元以上の場合は重回帰式,$y$が離散値の場合のタスクは分類と呼ばれます。単回帰分析は,単回帰式を用いて分析を行う手法のことを指します。

回帰分析を実現する最も基本的な手法に,線形回帰があります。以下,行列を用いて定式化を行います。線形回帰では,目的変数の予測値$\hat{\mY}$を$\mX$の線形和で表します。

\begin{align}
\hat{\mY} &= \mX \vb
\end{align}

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

\begin{align}
\boldsymbol{\varepsilon} &= \mY – \hat{\mY} \\
&= \mY – \mX \vb
\end{align}

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

\begin{align}
\left[\begin{array}{c}
\varepsilon_{1} \\
\varepsilon_{2} \\
\vdots \\
\varepsilon_{n}
\end{array}\right]
&=
\left[\begin{array}{c}
y_{1} \\
y_{2} \\
\vdots \\
y_{n}
\end{array}\right]
\left[\begin{array}{ccccc}
1 & x_{11} & x_{12} & \ldots & x_{1 d} \\
1 & x_{21} & x_{22} & \ldots & x_{2 d} \\
\vdots & \vdots & \vdots & & \vdots \\
1 & x_{n 1} & x_{n 2} & \ldots & x_{n d}
\end{array}\right]\left[\begin{array}{c}
b_{0} \\
b_{1} \\
\vdots \\
b_{d}
\end{array}\right]
\end{align}

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

\begin{align}
Q &= \boldsymbol{\varepsilon}^T \boldsymbol{\varepsilon} \\
&= (\mY – \mX \vb)^T (\mY – \mX \vb)\\
&= \mY^T \mY -2 \mY^T \mX \vb + \vb^T \mX^T \mX \vb
\end{align}

二乗和$Q$は凸な二次関数であるため,$Q$の微分を0と置くことで各パラメータを推定することができます。ただし,結果が正規方程式となるように偏微分の係数を1/2とすることで辻褄を合わせます。

\begin{align}
\frac{1}{2}\frac{\partial Q}{\partial \vb}
&= \mX^T \mX \vb – \mX^T \mY \\
&= 0
\end{align}

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

\begin{align}
\vb &= (\mX^T \mX)^{-1} \mX^T \mY
\end{align}

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

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

\begin{align}
\left[\begin{array}{cc}
N & N \bar{x} \\
N \bar{x} & \vx^{T} \vx
\end{array}\right]\left[\begin{array}{l}
b_{0} \\
b_{1}
\end{array}\right]=\left[\begin{array}{c}
N \bar{y} \\
\vx^{T} \vy
\end{array}\right]
\end{align}

具体的に計算を進めていきます。

\begin{align}
b_{0} &= \bar{y} -b_1\bar{x}\\
b_{1}&=\frac{\vx^{T} \vy-N \bar{x} \bar{y}}{\vx^{T} \vx-N \bar{x}^{2}}\\
&=\frac{(\vx-\bar{x} \mathbf{I})^{T}(\boldsymbol{y}-\bar{y} \mathbf{I})}{(\boldsymbol{x}-\bar{x} \mathbf{I})^{T}(\boldsymbol{x}-\bar{x} \mathbf{I})}\\
&=\frac{S_{x y}}{S_{x x}}
\end{align}

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

\begin{align}
S_{\alpha \beta} &= \sum_{i=1}^N (\alpha_n – \bar{\alpha}) (\beta_n – \bar{\beta}) \\
\end{align}

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

\begin{align}
r &= \frac{\mathrm{cov}[\vx, \vy]}{\sqrt{V[\vx]} \sqrt{V[\vy]}}\\
&= \frac{S_{xy}/N}{\sqrt{S_{xx}/N} \sqrt{S_{yy}/N}}\\
&= \frac{S_{xy}}{\sqrt{S_{xx}} \sqrt{S_{yy}}}
\end{align}

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

\begin{align}
b_1 &= \frac{S_{xy}}{S_{xx}} \\
&= \frac{S_{xy}}{\sqrt{S_{xx}}\sqrt{S_{xx}}}\\
&= \frac{S_{xy}}{\sqrt{S_{xx}}\sqrt{S_{yy}}}\\
&= r
\end{align}

したがって,観測データが標準化された場合に限り,単回帰の偏回帰係数$b_1$は相関係数$r$と等しくなります。

以上の議論から,偏回帰係数は説明変数の目的変数に対する影響力の大きさを表すことが分かります。ただし,重回帰の場合は偏回帰係数の解釈には注意が必要であり,後述する多重共線性が生じるために注意が必要です。

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

\begin{align}
\log p(\boldsymbol{\varepsilon}) &=
\log \prod_{i=1}^n p(\varepsilon_i)\\
&= \sum_{i=1}^n \log p(\varepsilon_i) \\
&\propto \sum_{i=1}^n \log \exp\left( -\frac{\varepsilon_i^2}{2\sigma^2} \right) \\
&\propto – \sum_{i=1}^n \varepsilon_i^2 \\
&= – \boldsymbol{\varepsilon}^T \boldsymbol{\varepsilon}\\
&= -Q
\end{align}

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

Come back 多重共線性

ここまでの議論は$\mX^T \mX$が正則な場合を考えました。しかし,$\mX^T \mX$が正則でない場合には多重回帰には問題が生じます。

$\mX^T \mX$が正則でないとき,$\mX^T \mX$はランク落ちの可能性があります。これはつまり,説明変数$\mX$がランク落ちであるという状況(説明変数同士に相関がある/説明変数が線形従属である)に相当します。

このとき,無理矢理$\mX^T \mX$の逆行列を計算しようとすると,本来逆行列を持たないはずの値を計算しているために,$\vb$の値が不安定に(大きく)なってしまいます。この問題を多重共線性と呼びます。

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

\begin{align}
\mathcal{L} &= |\mY – \mX \vb|^2 + \lambda |\vb|^q
\end{align}

特に,$q=0$のときは素朴な最小二乗法(重回帰),$q=1$のときはLASSO,$q=2$のときはリッジ回帰と呼ばれます。正則化を施すことにより,回帰係数が不安定になる(大きくなりすぎる)ことを防ぐことができます。

リッジ回帰の偏回帰係数は,素朴な最小二乗法のときと同様にして,以下のように求められます。

\begin{align}
\vb &= (\mX^T \mX + \lambda \mI)^{-1} \mX^T \mY
\end{align}

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

リッジ回帰

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

データ作成

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

  • 真の回帰式を$y=\vx^T\va$とする
  • $\vx \sim \mathrm{N}(\bf{0},\bf{I})$として大元になる説明変数をサンプリングする
  • 大元になる観測データ$D=$($\vx$, $y$)を$y=\vx^T\va+\varepsilon$,$\varepsilon \sim \mathrm{N}(0,1)$としてサンプリングする
  • $\vx$と線形従属となる説明変数$\vx’$を$\vx’ = c\vx + \boldsymbol{\varepsilon}$としてサンプリングする
  • $y’=\vx’^T\va + \varepsilon$として$y’$を生成する
  • 上記操作を$l$回繰り返す(結果として$n=l+1$となる)

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

パラメータ設定

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

\begin{align}
(\va, \vc, n, d, l) &= ([2,\ldots,2], 1/2, 2, 10, 1)
\end{align}

このパラメータに従って,$\vx$,$\vx’$,$y$,$y’$を作成します。

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])

\begin{align}
\vx &= [ 1.00, 1.02, 2.27, -0.50 , 0.60 , 0.53, -0.20 , 0.06, 1.41, 0.40 , -0.72] \\
\vx’ &= [ 1.00, 0.66, 5.19, -0.28, 0.73, 2.36, -0.20, -1.24, 1.50, 1.93, -0.92] \\
y &= 11.21 \\
y’ &= 18.65
\end{align}

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

実験結果

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

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

\begin{align}
&\vb = [ 23.75, 0.60, -8.92, 84.07, 32.31, 1.60, -75.44, -13.45, 20.35, -9.30, -17.07]\\
E[|&\vb|] = 26.08 \\
y &= 33.8\\
y’ &= 41.1
\end{align}

説明変数間が線形従属であるため,$\mX^T \mX$がランク落ちして逆行列が計算できないような状況であるため,正しく$y$と$y’$が推定できていないことが分かります。

$\lambda=1$の場合,つまりリッジ回帰の場合は,変回帰係数の値と絶対値の平均は以下のようになります。

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

\begin{align}
&\vb = [ 0.92, 0.91, 2.25, -0.45, 0.55, 0.59, -0.19, -0.02, 1.30, 0.46, -0.67]\\
E[&\vb] = 0.51 \\
y &= 10.4\\
y’ &= 18.6
\end{align}

この結果から,素朴な最小二乗法と比べて偏回帰係数の大きさがしっかりと抑制されていることが分かります。さらに,$y$と$y’$もよく推定できていることも分かります。

一方,$\lambda$を大きくすると,より正則化の効果が強まります。例えば,$\lambda=10^3$の場合は,変回帰係数の値と絶対値の平均は以下のようになります。

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

\begin{align}
&\vb = [ 0.03, 0.02, 0.12, -0.01, 0.02, 0.05, -0.01, -0.02, 0.04, 0.04, -0.02] \\
E[&\vb] = 0.02 \\
y &= 0.44\\
y’ &= 0.96
\end{align}

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

シェアはこちらからお願いします!
URLをコピーする
URLをコピーしました!

コメント

コメントする

目次
閉じる