Open5ヶ月前にコメント追加2

多重共線性と線形回帰係数推定量の分布

完全な多重共線性がある場合の推定量の分布

$p$ 次元の説明変数ベクトル $\boldsymbol{x}_i$ と目的変数 $y_i$ が $n$ 組について、線形回帰式

y_i = \boldsymbol{x}_i^\top \boldsymbol{\beta} + \varepsilon_i

が与えられたとき、 $\boldsymbol{Y}=(y_1, y_2, ..., y_n)^\top$ 、 $\boldsymbol{X}=(\boldsymbol{x}_1, \boldsymbol{x}_2, ..., \boldsymbol{x}_n)^\top$ 、 $\boldsymbol{\varepsilon} = (\varepsilon_1, \varepsilon_2, ..., \varepsilon_n)^\top$ を用いて以下のように表すことができる：

\boldsymbol{Y} = \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}. \tag{1}

ここで、 $\{\varepsilon_i\}_{i=1}^n$ は分散 $\sigma^2$ の誤差項である。

このとき、最小二乗法による係数 $\boldsymbol{\beta}$ の推定量は、以下のように表すことができる：

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

$\boldsymbol{X}^+$ は $\boldsymbol{X}$ の Moore-Penrose 逆行列、 $\boldsymbol{I}_p$ は $p\times p$ の単位行列、 $\boldsymbol{w}$ は任意の $p$ 次元ベクトルである。

説明変数間に完全な多重共線性がある場合には係数推定値は一意に定まらないが、上記の式(2)は同じ予測値を与える同値な全ての係数推定値を表現できる。

今回は、誤差項が独立同時な正規分布に従うと仮定したうえで、そのような完全な多重共線性がある場合の推定量の分布について見ていく。

結果

誤差項 $\{\varepsilon_i\}_{i=1}^n$ は正規分布に従い、その分散を $\sigma^2$ とする（ $\varepsilon_i \sim N(0, \sigma^2)$ ）。

このとき、式(2)で表された最小二乗推定量 $\hat{\boldsymbol{\beta}}$ の期待値と共分散行列は以下のように求まる：

E[\hat{\boldsymbol{\beta}}] = \boldsymbol{X}^+ \boldsymbol{X}\boldsymbol{\beta} + (\boldsymbol{I}_p - \boldsymbol{X}^+ \boldsymbol{X})\boldsymbol{w}, \tag{3}

V[\hat{\boldsymbol{\beta}}] = \sigma^2 \boldsymbol{X}^+ (\boldsymbol{X}^+)^\top. \tag{4}

導出

式(2)に式(1)を代入し、

\hat{\boldsymbol{\beta}} = \boldsymbol{X}^+ \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{X}^+ \boldsymbol{\varepsilon} + (\boldsymbol{I}_p - \boldsymbol{X}^+ \boldsymbol{X})\boldsymbol{w}

従って、両辺に期待値を取れば式(3)が成り立つ。

次に、共分散行列については、

\begin{aligned} V[\hat{\boldsymbol{\beta}}]&=E[\boldsymbol{X}^+ \boldsymbol{\varepsilon}(\boldsymbol{X}^+ \boldsymbol{\varepsilon})^\top]\\ &= \boldsymbol{X}^+ E[\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^\top](\boldsymbol{X}^+)^\top \end{aligned}

であるので、 $E[\boldsymbol{\varepsilon}\boldsymbol{\varepsilon}^\top]=\sigma^2 \boldsymbol{I}_n$ であることから式(4)が成り立つ。

なお、 $\hat{\boldsymbol{\beta}}$ の分布は、 $(\boldsymbol{I}_p - \boldsymbol{X}^+\boldsymbol{X})\boldsymbol{w}$ と直交する方向に広がる。

証明

V[\hat{\boldsymbol{\beta}}] (\boldsymbol{I}_p - \boldsymbol{X}^+\boldsymbol{X})\boldsymbol{w} = \boldsymbol{0}

であり、 $\hat{\boldsymbol{\beta}}$ の共分散行列が作る行ベクトル空間（列ベクトル空間と一致する）が $(\boldsymbol{I}_p - \boldsymbol{X}^+\boldsymbol{X})\boldsymbol{w}$ と直交していることからわかる。
これを示すには、

(\boldsymbol{X}^+)^\top (\boldsymbol{I}_p - \boldsymbol{X}^+\boldsymbol{X}) = \boldsymbol{0}

が分かればよい。

\begin{aligned} (\boldsymbol{X}^+)^\top \boldsymbol{X}^+ \boldsymbol{X} &= (\boldsymbol{X}^+)^\top (\boldsymbol{X}^+ \boldsymbol{X})^\top \qquad \because (\boldsymbol{X}^+ \boldsymbol{X})^\top = \boldsymbol{X}^+ \boldsymbol{X}\\ &= \{(\boldsymbol{X}^+ \boldsymbol{X})\boldsymbol{X}^+\}^\top\\ &= (\boldsymbol{X}^+)^\top \qquad \because \boldsymbol{X}^+ \boldsymbol{X} \boldsymbol{X}^+ = \boldsymbol{X}^+ \end{aligned}

よって、

(\boldsymbol{X}^+)^\top (\boldsymbol{I}_p - \boldsymbol{X}^+\boldsymbol{X}) = (\boldsymbol{X}^+)^\top - (\boldsymbol{X}^+)^\top = \boldsymbol{0}.

具体例

定数項あり・2変数（ $p=3$ ）の線形回帰を考える：

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

ただし、 $x_{2,i} = 2x_{1,i}$ という関係が常に成り立っていて、完全な多重共線性をもつものとする。

ここで、係数を $\beta_0=1, \beta_1=1, \beta_2=1$ 、誤差項の標準偏差を $\sigma=1.5$ として目的変数 $\{y_i\}_{i=1}^n$ を生成し、それをもとに係数の最小二乗推定量 $\hat{\beta}_1, \hat{\beta}_2$ を求めた：

青点が $\boldsymbol{w}=0$ とした最小二乗推定量 $(\hat{\beta}^{(0)}_1, \hat{\beta}^{(0)}_2)$ で、青の実線が式(2)で表される推定量全体である。

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

すると、青点で表した推定量 $(\hat{\beta}^{(0)}_1, \hat{\beta}^{(0)}_2)$ は、直線状に散らばっていることが見て取れる。

きちんと計算すると、係数推定量 $\hat{\beta}_1, \hat{\beta}_2$ の期待値は

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

のように表され、これは $E[\hat{\beta}_1] + 2 E[\hat{\beta}_2] = 3$ で表される直線上にあることがわかる（図：青の点線）。
このうち $(E[\hat{\beta}^{(0)}_1], E[\hat{\beta}^{(0)}_2]) = (3/5, 6/5)$ である（赤の＋印でプロット）。
なお、係数の「真の」値 $(\beta_1, \beta_2)=(1,1)$ は赤の★印で表した。こちらも上記の直線上に乗っていることがわかる。

次に推定量の分散を見る。共分散行列は

\begin{aligned} V\left[ \begin{pmatrix} \hat{\beta}_1\\ \hat{\beta}_2 \end{pmatrix} \right]&= \frac{\sigma^2}{25 \sum_i (x_{1,i}- \bar{x}_1)^2} \begin{pmatrix} 1 & 2\\ 2 & 4 \end{pmatrix}\\ &=\frac{1}{\sqrt{5}} \begin{pmatrix} 1 & -2\\ 2 & 1 \end{pmatrix} \begin{pmatrix} \frac{\sigma^2}{5\sum_i(x_{1,i}-\bar{x}_1)^2} & 0\\ 0 & 0 \end{pmatrix} \frac{1}{\sqrt{5}} \begin{pmatrix} 1 & 2\\ -2 & 1 \end{pmatrix} \tag{6} \end{aligned}

と表されることから、推定量は $\frac{1}{\sqrt{5}} (1, 2)^\top$ 方向に大きさ $\frac{\sigma^2}{5\sum_i(x_{1,i} - \bar{x}_1)^2}$ の分散を持っていることがわかる。
この分散をもとに計算した95%信頼区間を、図で影付きで表した。

計算方法

式(5)、(6)の計算方法について記述する。

$\boldsymbol{x}_1 = (x_{1,1}, x_{2,1}, ..., x_{n,1})^\top$ とすると、計画行列は $\boldsymbol{X} = (\boldsymbol{1}_n, \boldsymbol{x}_1, 2\boldsymbol{x}_1)$ のように表すことができる（ $\boldsymbol{1}_n$ は $n$ 次元で成分が全て1のベクトル）。

このとき残差回帰により、

\begin{aligned} \tilde{\boldsymbol{X}} &= \begin{pmatrix} x_{1,1} - \bar{x}_1 & x_{1,2} - \bar{x}_2\\ x_{2,1} - \bar{x}_1 & x_{2,2} - \bar{x}_2\\ \vdots & \vdots\\ x_{n,1} - \bar{x}_1 & x_{n,2} - \bar{x}_2\\ \end{pmatrix}\\ &= (\boldsymbol{x}_1 - \bar{x}_1 \boldsymbol{1}_n, 2\boldsymbol{x}_1 - 2\bar{x}_1 \boldsymbol{1}_n) \end{aligned}

を用いると、係数 $\beta_1, \beta_2$ の最小二乗推定量は、以下のように表すことができる：

\begin{pmatrix} \hat{\beta}_1\\ \hat{\beta}_2 \end{pmatrix} = \tilde{\boldsymbol{X}}^+\boldsymbol{Y} + (\boldsymbol{I}_2 - \tilde{\boldsymbol{X}}^+ \tilde{\boldsymbol{X}})\boldsymbol{w}.

従って、式(3)、(4)に従い、

E\left[ \begin{pmatrix} \hat{\beta}_1\\ \hat{\beta}_2 \end{pmatrix} \right] = \tilde{\boldsymbol{X}}^+ \tilde{\boldsymbol{X}} \begin{pmatrix} \beta_1\\ \beta_2 \end{pmatrix}+ (\boldsymbol{I}_2 - \tilde{\boldsymbol{X}}^+ \tilde{\boldsymbol{X}}) \begin{pmatrix} w_1\\ w_2 \end{pmatrix} \tag{7}

V\left[ \begin{pmatrix} \hat{\beta}_1\\ \hat{\beta}_2 \end{pmatrix} \right]=\sigma^2 \tilde{\boldsymbol{X}}^+(\tilde{\boldsymbol{X}}^+)^\top \tag{8}

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

Moore-Penrse 逆行列 $\tilde{\boldsymbol{X}}^+$

ここで、 $\tilde{\boldsymbol{X}}^+$ を求めるために $\tilde{\boldsymbol{X}}$ の特異値分解を考える。
${\rm rank}\,\tilde{\boldsymbol{X}}=1$ であるため、 $\tilde{\boldsymbol{X}}$ は1個の0でない特異値を持つ。
ここで、 $2\times 2$ 行列 $\tilde{\boldsymbol{X}}^\top\tilde{\boldsymbol{X}}$ について見ると、

\begin{aligned} \tilde{\boldsymbol{X}}^\top\tilde{\boldsymbol{X}}&= \sum_i (x_{1,i} - \bar{x}_1)^2 \begin{pmatrix} 1 & 2\\ 2 & 4 \end{pmatrix}\\ &=\frac{1}{\sqrt{5}} \begin{pmatrix} 1 & -2\\ 2 & 1 \end{pmatrix} \begin{pmatrix} 5\sum_i(x_{1,i}-\bar{x}_1)^2 & 0\\ 0 & 0 \end{pmatrix} \frac{1}{\sqrt{5}} \begin{pmatrix} 1 & 2\\ -2 & 1 \end{pmatrix} \end{aligned}

であることから固有値 $\lambda_1 = 5\sum_i(x_{1,i}-\bar{x}_1)^2$ とそれに対応する固有ベクトル $\boldsymbol{v}_1 = \frac{1}{\sqrt{5}}(1,2)^\top$ を持つことがわかる。

さらに、 $n \times n$ 行列 $\tilde{\boldsymbol{X}}\tilde{\boldsymbol{X}}^\top$ も固有値 $\lambda_1$ をもち、それに対応する固有ベクトル $\boldsymbol{u}_1$ は以下のように表すことができる：

\begin{aligned} \boldsymbol{u}_1 &= \frac{1}{\sqrt{\lambda_1}} \tilde{\boldsymbol{X}} \boldsymbol{v}_1\\ &= \frac{1}{\sqrt{\sum_i(x_{1,i}-\bar{x}_1)^2}} (\boldsymbol{x}_1 - \bar{x}_1 \boldsymbol{1}_n). \end{aligned}

以上から、 $\tilde{\boldsymbol{X}}$ の特異値分解は

\begin{aligned} \tilde{\boldsymbol{X}} = \boldsymbol{u}_1 \sqrt{\lambda_1} \boldsymbol{v}_1^\top \end{aligned}

であり、従ってその Moore-Penrose 逆行列は

\begin{aligned} \tilde{\boldsymbol{X}}^+ &= \boldsymbol{v}_1 \frac{1}{\sqrt{\lambda_1}} \boldsymbol{u}_1^\top\\ &=\frac{1}{5\sum_i(x_{1,i}-\bar{x}_1)^2} \begin{pmatrix} \boldsymbol{x}_1^\top - \bar{x}_1 \boldsymbol{1}_n^\top\\ 2\boldsymbol{x}_1 ^\top- 2\bar{x}_1 \boldsymbol{1}_n^\top\\ \end{pmatrix} \end{aligned}

のように求まる。

推定量期待値

上記より

\tilde{\boldsymbol{X}}^+ \tilde{\boldsymbol{X}} = \frac{1}{5} \begin{pmatrix} 1 & 2\\ 2 & 4 \end{pmatrix}

であることから、式(7)にこれと $\beta_1=1, \beta_2=1$ を代入すると

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

となり、 $\alpha = (-2w_1 + w_2)/5$ と置き換えると式(5)が得られる。

推定量共分散行列

式(8)に先ほどの $\tilde{\boldsymbol{X}}^+$ を代入すればよい：

\begin{aligned} V\left[ \begin{pmatrix} \hat{\beta}_1\\ \hat{\beta}_2 \end{pmatrix} \right]&=\sigma^2 \tilde{\boldsymbol{X}}^+(\tilde{\boldsymbol{X}}^+)^\top\\ &=\sigma^2\boldsymbol{v}_1 \frac{1}{\sqrt{\lambda_1}} \boldsymbol{u}_1^\top\boldsymbol{u}_1 \frac{1}{\sqrt{\lambda_1}} \boldsymbol{v}_1^\top\\ &= \boldsymbol{v}_1 \frac{\sigma^2}{\lambda_1}\boldsymbol{v}_1^\top\\ &=\frac{1}{\sqrt{5}} \begin{pmatrix} 1 \\ 2 \end{pmatrix} \frac{\sigma^2}{5\sum_i(x_{1,i}-\bar{x}_1)^2} \frac{1}{\sqrt{5}} \begin{pmatrix} 1 & 2\\ \end{pmatrix}\\ &= \frac{\sigma^2}{25 \sum_i (x_{1,i}- \bar{x}_1)^2} \begin{pmatrix} 1 & 2\\ 2 & 4 \end{pmatrix}\\ \end{aligned}

これにより式(6)が求まった。

コード

シミュレーション・作図に用いた Python コードは以下にある：

畳屋民也 5ヶ月前

完全な多重共線性がある場合の Ridge 回帰推定量の分布

（とりあえずグラフだけ貼っておく）