「多重共線性があると推定が不安定になる」とは、係数推定量の分散が大きくなることで、入力データ次第で係数推定値が大きく（時に正負さえ）変わってしまうことを意味する
多重共線性があってもサンプルサイズを増やせば推定量分散が小さくなり、不安定性は軽減する
Ridge 回帰は推定量分散を小さくするが、代わりにバイアスが増大し係数推定値が真の値から離れていく。

多重共線性と推定量

多重共線性とは？

多重共線性とは何かを説明するにあたり、まずは「完全な多重共線性」と呼ばれるものについて解説する。

完全な多重共線性とは、説明変数同士が線型結合で書き表せてしまうことを意味する。
もう少しわかりやすく説明すると、例えば説明変数の組 $(X_1, X_2, X_3, X_4)$ があったとき $X_1$ が $(X_2, X_3, X_4)$ を用いて $X_1 = \alpha_2 X_2 + \alpha_3 X_3 + \alpha_4 X_4$ のような形で表せてしまうような状況を指す。

具体例として、 $X_1$ が円単位でのある商品の価格に対して、 $X_2$ が同じ商品のドル単位での価格であったとする。そうすると、 $X_1$ と $X_2$ は実質的に同一の意味を持った説明変数であり、2023/12/10時点で $X_1 = 145 X_2$ という関係が成り立っていて完全な多重共線性が認められる。

一方、説明変数が別の説明変数の線型結合で「ピッタリ」説明できてしまうまで行かなくても、「大体」説明できてしまうという場合も考えられる。
上記の商品価格の例で言うと、為替レートは日々変動するため $X_1 = 145 X_2$ という関係は常には成り立たないものの、2023年でいうと大体1ドル130~150円の間で推移しているので大まかに $X_1 \sim 140 X_2$ という関係が成り立つと見做せる。^[1]
理論的には特定の説明変数同士の相関が高かったり、ある説明変数について他の説明変数いくつかを用いて線形回帰を行ったみたところ決定係数が1に近い、といった状況が当てはまる。

このように完全な多重共線性ではないがそれに近い性質を「完全でない多重共線性」と呼ぶ。

本記事では主にこのような「完全でない多重共線性」を扱うため、以降「多重共線性」と言った時には特に断りのない限り完全でないケースを指すものとする。

多重共線性が線形回帰に及ぼす影響

では、多重共線性があると線形回帰において何が問題だろうか？

冒頭で述べたように、多重共線性があると回帰係数の推定が「不安定」になる。
それが何を意味するかというと、わずかな入力データの違いで推定結果が大きく変わってしまう。
このことを簡単な例で示してみる。

先ほどの例を参考に、説明変数 $X_1$ を円単位の月収、 $X_2$ をドル単位の月収とし、これらをもとに目的変数として家賃額 $Y$ を予測することを考える。

ID	$X_1$ （月収、万円）	$X_2$ （月収、ドル）	参考：為替レート（1ドルあたり円）	$Y$ （家賃、月当たり、万円）
1	30	2,142	140	11.0
2	20	1,418	141	7.9
3	15	1,079	139	6.5
4	45	3,260	138	15.5
5	25	1,760	142	9.4

上記のような架空のデータセットが与えられたとする。これをもとに、

Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \varepsilon

という線形回帰モデルを立てる（ $\varepsilon$ は誤差項）。
それにより、「円単位で見た月収が高い人ほど、どれだけ支払う家賃が変わるか？」という問いに答えたい。

上記のモデルをもとに最小二乗法で定数項・回帰係数の推定を行なってみると、それぞれ $(\hat{\beta}_0, \hat{\beta}_1, \hat{\beta}_2) = (1.99, 0.163, 1.89\times 10^{-3})$ であった。
これを素直に解釈すると、 $\hat{\beta}_1=0.163$ から「月収が10,000円高くなると家賃がその16.3%つまり1,630円高いところに住む傾向がある」ということになりそうである。

さて、ここまで予測を行ったところで家賃のデータに誤りがあることがわかり、正しいデータは以下のとおりであった：

ID	$Y$ （家賃、月当たり、万円）
1	10.9
2	8.0
3	6.6
4	15.6
5	9.4

これをもとに計算し直すとどうだろう、推定値は $(\hat{\beta}_0, \hat{\beta}_1, \hat{\beta}_2) = (2.17, -0.0510, 4.82\times 10^{-3})$ のようになった。
$X_1$ の係数 $\hat{\beta}_1$ は、先ほどの0.163から打って変わってマイナスの値-0.0510になってしまった。まさか月収が高い人ほど家賃が安い家に住むというわけもあるまい。

そこで再度家賃データを見直して修正したところ、以下のようになった：

ID	$Y$ （家賃、月当たり、万円）
1	11.2
2	8.0
3	6.5
4	15.2
5	9.5

これをもとに再度計算し直すと、 $(\hat{\beta}_0, \hat{\beta}_1, \hat{\beta}_2) = (2.07, 0.576, -3.91\times 10^{-3})$ 、今度は $\hat{\beta}_1=0.567$ であり家賃が月収増加額に対して57.6%上昇するという随分大きな値が出てしまった。

このように、入力データが少し変わっただけで、線形回帰係数の推定値が大きく変化してしまう。
この原因は、説明変数として円単位の月収 $X_1$ とほぼ同じ意味を持ち相関も強い $X_2$ 、つまりドル単位の月収も含まれていることにある。
なお、もし $X_1$ だけ用いて予測を行っていれば、いずれの家賃データを入力として用いても回帰係数 $\hat{\beta}_1$ の推定値はおおよそ 0.30、つまり「月収が10,000円高ければ家賃がその約30%の3,000円程度高いところに住む傾向にある」という結果になる。

このように、多重共線性があるとほんの少しの入力データの違いで係数推定値が大きく変わってしまい、このことを「推定が不安定になる」と言い表す。

数値実験：　2変数線形回帰における多重共線性と推定量分布

前節で多重共線性があるとわずかな入力データの違いで係数推定値が大きく変化してしまうことを見たが、このことは理論的には推定量の分散が大きいことを意味する。

以下、多重共線性がある場合の推定量の分布の性質について、2変数の重回帰による数値実験例をもとに時折数式も交えながら見ていく。

実験設定

以下のような定数項あり・2変数の線形回帰を考える：

y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \varepsilon_i.

$\varepsilon_i$ は期待値0で分散 $\sigma^2$ の独立同時な正規分布 $N(0, \sigma^2)$ に従う誤差項である。

今回の実験では、真の定数項・回帰係数を $\beta_0=1, \beta_1=1, \beta_2=1$ 、誤差項分散を $\sigma^2=(1.5)^2$ と置いた。

このような設定のもと、サンプルサイズ $n$ のデータセット $\{(x_{1i}, x_{2i}, y_i)\}_{i=1}^n$ を生成して最小二乗法にて回帰係数推定量 $(\hat{\beta}_1, \hat{\beta}_2)$ を求める操作を考える。
この操作を異なる乱数で生成した誤差項 $\{ \varepsilon_i \}_{i=1}^n$ を用いて1000回繰り返した時に、推定量 $(\hat{\beta}_1, \hat{\beta}_2)$ がどのような分布に従うかをシミュレーションした。

相関係数 $\rho$ をパラメータとした多重共線性の強さのコントロール

上記に加えて、説明変数 $x_1, x_2$ が多重共線性を持つ場合を考えるために、以下のような設定を採り入れた：

条件付き期待値で見ると、 $x_{2i}$ は $x_{1i}$ の2倍（ $E[x_{2i} \vert x_{1i}] = 2 x_{1i}$ ）
$x_1$ と $x_2$ の相関係数 $\rho$ はあらかじめ決められた値で、パラメータとしてコントロールする（ $0<\rho \le 1$ ）^[2]

なおこの設定では、 $\rho$ が1に近いほど多重共線性が強く（完全に近く）なる。

上記の設定を実現するには、サンプル $\{(x_{1i}, x_{2i})\}_{i=1}^n$ を以下のように生成すればよい：

$\{x_{1i}\}_{i=1}^n$ は期待値0分散 $\sigma_1^2$ の正規分布に従うように生成
$\{x_{2i}\}_{i=1}^n$ は正規分布に従う誤差項 $\eta_i$ を用いて $x_{2i} = 2 x_{1i} + \eta_i$ のように生成
- ただし、 $\eta_i$ の分散は $\sigma_{\eta}^2 = 4\sigma_1^2(1/\rho^2 - 1)$

なぜこのように生成するか？

$x_2$ の生成方法から、 $E[x_2\vert x_1] = 2x_1$ が成り立つことはすぐわかる。
$x_2 = 2 x_1 + \eta$ に両辺 $E[\cdot \vert x_1]$ をとればよい。

以下、 $x_1, x_2$ の相関係数が $\rho$ になることを確認する。

$x_1, x_2$ の分散 $V[x_1], V[x_2]$ と共分散 ${\rm Cov}[x_1, x_2]$ をそれぞれ求めると以下のようになる。

V[x_1] = \sigma_1^2

\begin{aligned} V[x_2] &= E[(2x_1 + \eta - 0)^2]\\ &= 4E[x_1^2] + E[\eta^2]\\ &= 4\sigma_1^2 + \sigma_{\eta}^2\\ &=4\sigma_1^2/\rho^2 \end{aligned}

\begin{aligned} {\rm Cov}[x_1, x_2] &= E[x_1x_2] - E[x_1]E[x_2]\\ &= E[x_1 (2x_1 + \eta)] - 0\\ &= 2E[x_1^2]\\ &= 2\sigma_1^2 \end{aligned}

したがって、 $x_1$ と $x_2$ の相関係数は、

\frac{{\rm Cov}[x_1, x_2]}{\sqrt{V[x_1]}\sqrt{V[x_2]}} = \rho

となる。^[3]

このような設定のもと、相関係数 $\rho$ を変えて多重共線性の強さを変化させるとともに回帰係数推定量 $(\hat{\beta}_1, \hat{\beta}_2)$ の分布がどのように変わっていくかを調べた。

推定量分布の可視化

サンプルサイズを $n=10$ として、説明変数 $\{(x_{1i}, x_{2i})\}_{i=1}^n$ は固定で誤差項 $\{\varepsilon_i\}_{i=1}^n$ を変えて1,000回シミュレーションを行い得られた係数推定値 $(\hat{\beta}_1, \hat{\beta}_2)$ をプロットした。

これを相関係数を $\rho = 0.5, 0.75, 0.9, 0.95$ とした時の結果を並べたものが以下である。

青点が各試行で得られた回帰係数推定値 $(\hat{\beta}_1, \hat{\beta}_2)$ 、赤星マークが回帰係数の真の値 $(\beta_1, \beta_2) = (1,1)$ である。

灰色の楕円は推定量分布の95%領域であり、推定値の散らばり具合を表すために描いた。
なお、多重共線性があっても最小二乗法による推定量の不偏性は保たれるため、分布の中心は常に真の値 $(\beta_1, \beta_2) = (1,1)$ である。

推定量の分布と不偏性の数式を用いた理解

計画行列 $\boldsymbol{X}$ を

\begin{aligned} \boldsymbol{X} = \begin{pmatrix} 1 & x_{11} & x_{21}\\ 1 & x_{12} & x_{22}\\ \vdots & \vdots & \vdots\\ 1 & x_{1,n} & x_{2,n} \end{pmatrix} \end{aligned}

とおく。

すると、最小二乗法による推定量 $\hat{\boldsymbol{\beta}} = (\hat{\beta}_0, \hat{\beta}_1, \hat{\beta}_2)^\top$ の期待値・分散共分散行列は以下のように表すことができる：^[4]

\begin{aligned} E[\hat{\boldsymbol{\beta}}] &= \boldsymbol{\beta}\\ V[\hat{\boldsymbol{\beta}}] &= \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1}. \end{aligned}

このうちひとつ目の式は、最小二乗法による推定量 $\hat{\boldsymbol{\beta}}$ が不偏推定量であることを意味している。

また、誤差項 $\varepsilon_i$ が正規分布に従う場合は、推定量 $\hat{\boldsymbol{\beta}}$ は以下のように正規分布する：

\hat{\boldsymbol{\beta}} \sim N(\boldsymbol{\beta}, \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1}).

図の灰色の楕円はこの分布に基づいて描いたものである。

以下、上記の関係式を示す。

計画行列 $\boldsymbol{X}$ を用いると、残差二乗和 $||\boldsymbol{Y} - \boldsymbol{X}\boldsymbol{\beta}||^2$ を $\boldsymbol{\beta}$ について最小化することで最小二乗推定量 $\hat{\boldsymbol{\beta}}$ は以下のように表すことができる：^[5]

\hat{\boldsymbol{\beta}} = (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{Y}.

ただし、 $\boldsymbol{Y}$ は目的変数ベクトル $\boldsymbol{Y}=(y_1, y_2, ..., y_n)^\top$ である。

ここで、誤差項ベクトル $\boldsymbol{\varepsilon} = (\varepsilon_1, \varepsilon_2, ..., \varepsilon_n)^\top$ を定義すると、重回帰モデルの式は以下のように表される：

\boldsymbol{Y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}

これを先ほどの最小二乗推定量の式に代入すると、

\hat{\boldsymbol{\beta}} = \boldsymbol{\beta} + (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\varepsilon}

のように表すことができる。

これについて両辺期待値をとると、 $E[\boldsymbol{\varepsilon}]=\boldsymbol{0}$ であることから

E[\hat{\boldsymbol{\beta}}] = \boldsymbol{\beta}

となり、最小二乗推定量 $\hat{\boldsymbol{\beta}}$ が不偏推定量であることがわかる。

分散共分散行列については、 $E[\boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^\top] = \sigma^2 \boldsymbol{I}_n$ に着目して以下のように求める：

\begin{aligned} V[\hat{\boldsymbol{\beta}}] &= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top E[\boldsymbol{\varepsilon} \boldsymbol{\varepsilon}^\top] \left\{ (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \right\}^\top\\ &= \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{X}(\boldsymbol{X}^\top \boldsymbol{X})^{-1}\\ &=\sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \end{aligned}

斜めに引いた灰色の破線は、記事末の補足にて後述するものの、現時点では補助線程度に見ていただきたい。

これを見ると、 $\rho=0.5$ で多重共線性が弱いときは $\hat{\beta}_1,\hat{\beta}_2$ は真の値 $(1,1)$ を中心に比較的等方的に広がっているものの、 $\rho$ の値が増大し多重共線性の強さが増すにつれて灰色点線方向に沿って分布が伸び広がっていくことがわかる。

このように、多重共線性の度合いが強く完全な多重共線性に近づくほど、異なる入力データに対して係数推定値のばらつきが大きくなることが視覚的に理解できる。
このような状況が「推定が不安定になる」ことの正体である。

多重共線性への対処法についてのあれこれ

このように、多重共線性があると係数推定値のばらつきが多くなった。

一方、多重共線性があっても安定的に推定を行うための対処方法として、以下の2つがよく知られている：

サンプルサイズを増やす
Ridge 回帰など正則化を含んだ手法を使う

それぞれ、どのような意味があるのだろうか？

以降では、前節の数値実験における $\rho=0.95$ の場合を用いて再度実験してみる。

サンプルサイズを増やした場合

サンプルサイズを増やすと、推定量の分散は小さくなっていき推定の安定性が増す。

以下は、先ほどの多重共線性が強かった $\rho=0.95$ の場合について、サンプルサイズを10から順に増やしていった場合の推定量分布を見たものである。

サンプルサイズが増えるにつれ、真の値 $(\beta_1, \beta_2)=(1,1)$ 周りの推定量 $(\hat{\beta}_1, \hat{\beta}_2)$ の散らばりが小さくなっていることが見てとれる。

したがって、多重共線性がある場合でもサンプルサイズを増やすことで安定した推定を行えるようになる。

理論的な理解

以下、サンプルサイズが増えた際に推定量の分散も小さくなることを数式を用いて確認する。

推定量 $\hat{\beta}_1, \hat{\beta}_2$ の分散はそれぞれ以下のように表すことができる：

V[\hat{\beta}_1] = \frac{1}{1 - \hat{\rho}^2} \frac{\sigma^2}{\sum_{i=1}^n (x_{1i} - \bar{x}_1)^2}, \quad V[\hat{\beta}_2] = \frac{1}{1 - \hat{\rho}^2} \frac{\sigma^2}{\sum_{i=1}^n (x_{2i} - \bar{x}_2)^2}.

ただし、 $\hat{\rho}$ は　 $x_1$ と $x_2$ の標本相関係数で以下のように定義した：

\hat{\rho} = \frac{\sum_{i=1}^n (x_{1i} - \bar{x}_1)(x_{2i} - \bar{x}_2)}{\sqrt{\sum_{i=1}^n (x_{1i} - \bar{x}_1)^2}\sqrt{\sum_{i=1}^n (x_{2i} - \bar{x}_2)^2}}.

ここで、 $\hat{\rho}$ は $n\to\infty$ と共に $\rho$ に確率収束する一方で、 $\sum_{i=1}^n (x_{1i} - \bar{x}_1)^2, \, \sum_{i=1}^n (x_{2i} - \bar{x}_2)^2$ はともにサンプルサイズ $n$ が増えると増大する量なので、 $n\to \infty$ に対して $V[\hat{\beta}_1] \to 0, \, V[\hat{\beta}_2] \to 0$ が言える。

したがって、サンプルサイズ $n$ が大きくなるほど推定量分散は小さくなる。

以下、推定量分散 $V[\hat{\beta}_1], \, V[\hat{\beta}_2]$ が上述の形で表せることを示す。

これまで通り、 $\boldsymbol{X}$ を計画行列、 $\boldsymbol{Y}$ を目的変数ベクトルとする。

ここで、 $x_1, x_2$ のサンプル平均をそれぞれ $\bar{x}_1, \bar{x}_2$ とする：

\bar{x}_1 = \frac{1}{n} \sum_{i=1}^n x_{1i}, \qquad \bar{x}_2 = \frac{1}{n} \sum_{i=1}^n x_{2i}.

これを用いて、中心化した計画行列 $\tilde{\boldsymbol{X}}$ を以下のように定義する：

\tilde{\boldsymbol{X}} = \begin{pmatrix} x_{11} - \bar{x}_1 & x_{21} - \bar{x}_2\\ x_{12} - \bar{x}_1& x_{22} - \bar{x}_2\\ \vdots & \vdots\\ x_{1,n} - \bar{x}_1& x_{2,n} - \bar{x}_2 \end{pmatrix}.

これを用いることで、回帰推定量 $(\hat{\beta}_1, \hat{\beta}_2)^\top$ は以下のように求められる：

\begin{aligned} V\left[ \begin{pmatrix} \hat{\beta}_1\\ \hat{\beta}_2 \end{pmatrix} \right] &= \sigma^2 (\tilde{\boldsymbol{X}}^\top \tilde{\boldsymbol{X}})^{-1}\\ &= \sigma^2 \begin{pmatrix} \sum_{i=1}^n (x_{1i} - \bar{x}_1)^2 & \sum_{i=1}^n (x_{1i} - \bar{x}_1)(x_{2i} - \bar{x}_2)\\ \sum_{i=1}^n (x_{1i} - \bar{x}_1)(x_{2i} - \bar{x}_2) & \sum_{i=1}^n (x_{2i} - \bar{x}_2)^2 \end{pmatrix}^{-1}\\ &=\frac{\sigma^2}{\sum_{i=1}^n (x_{1i} - \bar{x}_1)^2\sum_{i=1}^n (x_{2i} - \bar{x}_2)^2 (1 - \hat{\rho}^2)} \begin{pmatrix} \sum_{i=1}^n (x_{2i} - \bar{x}_2)^2 & -\sum_{i=1}^n (x_{1i} - \bar{x}_1)(x_{2i} - \bar{x}_2)\\ -\sum_{i=1}^n (x_{1i} - \bar{x}_1)(x_{2i} - \bar{x}_2) & \sum_{i=1}^n (x_{1i} - \bar{x}_1)^2 \end{pmatrix}. \end{aligned}

これを整理することで、上述の $V[\hat{\beta}_1], \, V[\hat{\beta}_1]$ が求められる。

なお、 $\hat{\beta}_1, \hat{\beta}_2$ の共分散は以下のように表すことができる：

{\rm Cov}[\hat{\beta}_1, \hat{\beta}_2] = - \frac{\hat{\rho}^2}{1 - \hat{\rho}^2} \frac{\sigma^2}{\sum_{i=1}^n (x_{1i} - \bar{x}_1)(x_{2i} - \bar{x}_2)}.

Ridge 回帰を用いた場合

Ridge 回帰では、正則化係数 $\alpha$ を大きくするほど推定量分散が小さくなり安定性は増す一方、バイアスが大きくなり推定値が真の値から外れていってしまう。
以下、このことをシミュレーション結果と共に示す。

Ridge 回帰では、推定量を求める際に残差二乗和に回帰係数 $\boldsymbol{\beta}$ の二乗ノルムに比例する項 $\alpha || \boldsymbol{\beta} ||^2$ を足したものを最小化する。
これにより、推定量 $\hat{\boldsymbol{\beta}}^{(\alpha)}$ が大きくなりすぎず安定して推定を行うことができるとされている。

理論的な解説

これまでと同様に計画行列 $\boldsymbol{X}$ 、目的変数ベクトル $\boldsymbol{Y}$ を用いて議論する。

Ridge 回帰では、以下のような損失関数を最小化する：

L(\boldsymbol{\beta}) = ||\boldsymbol{Y} - \boldsymbol{X} \boldsymbol{\beta}||^2 + \alpha ||\boldsymbol{\beta}||^2 \quad (\alpha>0: {\rm const.}).

$\partial L(\boldsymbol{\beta}) / \partial \boldsymbol{\beta} = \boldsymbol{0}$ の解として得られる推定量は、以下のようになる：^[6]

\hat{\boldsymbol{\beta}}^{(\alpha)} = (\boldsymbol{X}^\top \boldsymbol{X} + \alpha \boldsymbol{I}_3)^{-1} \boldsymbol{X}^\top \boldsymbol{Y}.

ここで、 $\boldsymbol{I}_3$ は $3 \times 3$ の単位行列である。

では、正則化項を加えることによって回帰係数推定量の分布に実際にどのような変化がもたらされるのだろうか？
これを見るために、正則化項の係数 $\alpha$ を変えた時の推定量分布をこれまでと同様の手順で可視化した。

その結果が以下である。

上記の図では、Ridge 回帰係数推定量 $(\hat{\beta}^{(\alpha)}_1, \hat{\beta}^{(\alpha)}_2)$ を青点で表している。

赤のx印は理論的に求めたRidge回帰係数推定量の期待値で、灰色の影は推定量分布の95%領域である。

理論式とその導出

Ridge回帰推定量 $\hat{\boldsymbol{\beta}}^{(\alpha)}$ の分布は、以下の期待値・分散共分散行列を持つ正規分布になる：

\begin{aligned} E\left[\hat{\boldsymbol{\beta}}^{(\alpha)}\right] &= (\boldsymbol{X}^\top \boldsymbol{X} + \alpha \boldsymbol{I}_3)^{-1} \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta}\\ V\left[\hat{\boldsymbol{\beta}}^{(\alpha)}\right] &= \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X} + \alpha \boldsymbol{I}_3)^{-1} \boldsymbol{X}^\top \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X} + \alpha \boldsymbol{I}_3)^{-1}. \end{aligned}

このことは、 $\hat{\boldsymbol{\beta}}^{(\alpha)} = (\boldsymbol{X}^\top \boldsymbol{X} + \alpha \boldsymbol{I}_3)^{-1} \boldsymbol{X}^\top \boldsymbol{Y}$ に $\boldsymbol{Y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}$ を代入することで

\hat{\boldsymbol{\beta}}^{(\alpha)} =(\boldsymbol{X}^\top \boldsymbol{X} + \alpha \boldsymbol{I}_3)^{-1} \boldsymbol{X}^\top \boldsymbol{X}\boldsymbol{\beta} + (\boldsymbol{X}^\top \boldsymbol{X} + \alpha \boldsymbol{I}_3)^{-1} \boldsymbol{X}^\top \boldsymbol{\varepsilon}

となることからわかる。

これを見ると、正則化係数 $\alpha$ が大きくなるにつれて推定量 $(\hat{\beta}^{(\alpha)}_1, \hat{\beta}^{(\alpha)}_2)$ の分布の広がりが小さくなっていることがわかる。
このことから、確かに正則化項を入れることで推定量のばらつきが小さくなり「推定が安定化」することが見てとれる。

しかし一方で、正則化係数 $\alpha$ が大きくなると共に分布の中心（赤x印）が赤★印で表した係数の真値 $(\beta_1, \beta_2)=(1,1)$ からずれていっていることに注意しなくてはならない。

実は、Ridge 回帰を行うと推定量の不偏性は失われてしまうのである。

理論： Ridge 回帰推定量が不偏推定量でないこと

先ほどの解説より

\begin{aligned} E\left[\hat{\boldsymbol{\beta}}^{(\alpha)}\right] &= (\boldsymbol{X}^\top \boldsymbol{X} + \alpha \boldsymbol{I}_3)^{-1} \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta}\\ &= (\boldsymbol{X}^\top \boldsymbol{X} + \alpha \boldsymbol{I}_3)^{-1} (\boldsymbol{X}^\top \boldsymbol{X} + \alpha \boldsymbol{I}_3) \boldsymbol{\beta} - (\boldsymbol{X}^\top \boldsymbol{X} + \alpha \boldsymbol{I}_3)^{-1} \alpha \boldsymbol{I}_3 \boldsymbol{\beta}\\ &= \boldsymbol{\beta} - \alpha (\boldsymbol{X}^\top \boldsymbol{X} + \alpha \boldsymbol{I}_3)^{-1} \boldsymbol{\beta} \end{aligned}

であるので、 $E\left[\hat{\boldsymbol{\beta}}^{(\alpha)}\right] \ne \boldsymbol{\beta}$ であり不偏推定量とはならず、 $- \alpha (\boldsymbol{X}^\top \boldsymbol{X} + \alpha \boldsymbol{I}_3)^{-1} \boldsymbol{\beta}$ のバイアスがかかっている。

以上から、 Ridge 回帰を適用すると推定量のばらつきが小さくなるが、それと引き換えに推定の不偏性が失われることで期待値が真の値から外れていってしまうことがわかった。

まとめ

この記事では、2変数の重回帰による例を交えながら、多重共線性が強い（完全に近い）ほど回帰係数推定量の分散が大きくなり分布が広がることを見てきた。
これによりわずかな入力データの違いで大きく異なる推定値が得られてしまうことから、「多重共線性があると推定が不安定になる」と呼ばれる。

これに対してサンプルサイズを大きくすると推定量の分散が小さくなり「安定した推定」が行えるようになる一方で、Ridge 回帰では正則化係数を大きくするにつれて推定の安定性と引き換えにバイアスが大きくなり推定量が真の値から離れていってしまうことも見た。

多重共線性はデータ分析を行ううえでよく見かける現象であり、その発生メカニズムが明快な場合を除いて対処方法に苦慮するケースが多いと思うが、本記事の内容を念頭に分析結果にどのような影響を及ぼすかを想像しつつ、安易に説明変数を間引いたり Ridge 回帰を持ち出すなどをせず目的に応じて慎重に分析を進めていきたい。

補足：完全な多重共線性がある場合の推定量とその分布

本文のシミュレーションで相関係数が $\rho= \pm 1$ のときは、完全な多重共線性があるケースに相当する。

このような場合、「回帰係数が求まらない・推定を行えない」と言われることもしばしば見かけるが、実際は係数が一通りに定まらないだけで、推定量自体は一般化逆行列を用いて求めることができる。

詳しい計算はこちらのメモに委ねるが、 $\rho=1$ の場合の係数推定量 $(\hat{\beta}_1, \hat{\beta}_2)$ の期待値は、本記事におけるシミュレーションの設定では以下のように表すことができる：

E\left[ \begin{pmatrix} \hat{\beta}_1\\ \hat{\beta}_2 \end{pmatrix} \right] = \begin{pmatrix} 3/5\\ 6/5 \end{pmatrix} + \lambda \begin{pmatrix} -2\\ 1 \end{pmatrix}

$\lambda$ は任意の定数である。
ここから定数 $\lambda$ を消去することで $E[\hat{\beta}_1] + 2 E[\hat{\beta}_2] = 3$ という関係が成立するため、これを $(\hat{\beta}_1, \hat{\beta}_2)$ のグラフ上で表すと斜めの直線となることがわかる。
これが、これまでのグラフで「補助線」として引いていた灰色の点線の正体である。

以下、推定量 $(\hat{\beta}_1, \hat{\beta}_2)$ の分布の95%領域と共に上記の関係をグラフに表した：

多重共線性のあるデータに対して（正則化項を加えず）回帰を行った際の推定量の分布は、完全な多重共線性に近づく（今回の実験では $\rho \to 1$ ）につれ、この直線とその周りの帯状の領域に漸近していっているものという見方ができる。

理論詳細

詳細は前述の別記事に委ねて、ここでは結果だけ簡単に紹介する。

完全な多重共線性がある場合の最小二乗法による回帰係数推定量は、任意の定数ベクトル $\boldsymbol{w}$ を用いて以下のように表すことができる：

\hat{\boldsymbol{\beta}} = \boldsymbol{X}^+ \boldsymbol{Y} + (\boldsymbol{I}_3 - \boldsymbol{X}^+ \boldsymbol{X})\boldsymbol{w}.

ここで、 $\boldsymbol{X}^+$ は計画行列 $\boldsymbol{X}$ についての Moore-Penrose 逆行列である。

この時、係数推定量の期待値と分散共分散行列はそれぞれ以下のようになる：

\begin{aligned} E[\hat{\boldsymbol{\beta}}] &= \boldsymbol{X}^+ \boldsymbol{X}\boldsymbol{\beta} + (\boldsymbol{I}_3 - \boldsymbol{X}^+ \boldsymbol{X})\boldsymbol{w},\\ V[\hat{\boldsymbol{\beta}}] &= \sigma^2 \boldsymbol{X}^+ (\boldsymbol{X}^+)^\top. \end{aligned}

本文中に記述した $(\hat{\beta}_1, \hat{\beta}_2)$ の期待値や分布はこれらをもとに計算した。

実験・可視化コード

本記事で行ったシミュレーション・可視化のコードは、以下にアップしている：

参考文献

浅野皙・中村二朗、「計量経済学」第2版（有斐閣、2009）

https://amzn.asia/d/2xFQebO

佐和隆光、回帰分析(新装版) (統計ライブラリー、朝倉書店、2020)

https://amzn.asia/d/2vsEmZ6

脚注

$X_1 \sim 140 X_2$ を線形回帰式として解釈しようとした場合、誤差項の分布が独立同時な正規分布にならないが、ここではあくまでも例なので気にしないものとする。 ↩︎
今回の設定では $x_2$ に対する $x_1$ の回帰係数が正のため、相関係数は必ず正になる（後述）。 ↩︎
今回の場合 $x_1$ にかかる回帰係数が正の値のため $\rho>0$ として式展開を行ったが、回帰係数が負で $\rho <0$ となる場合も同様の関係式が成り立つ。実際、 $x_2 = \alpha x_1 + \eta$ とおくと、 ${\rm Cov}[x_1, x_2] /\sqrt{V[x_1] V[x_2]} =\frac{\alpha}{|\alpha|}|\rho|$ となり、 $\alpha$ の正負と相関係数の正負が一致することがわかる。 ↩︎
https://zenn.dev/link/comments/db9e98ef7a0301 を参照 ↩︎
前注と同様 ↩︎
詳しくは https://zenn.dev/link/comments/b36e340783fcad を参照。 ↩︎

この記事は？

サマリー