残差回帰（もしくは Frisch-Waugh-Lovell 定理）

はじめに

重回帰モデルが

$$yi=β0+β1x1,i+β2x2,i+β3x3,i+εi y_i = \beta_0 + \beta_1 x_{1,i} + \beta_2 x_{2,i} + \beta_3 x_{3, i} + \varepsilon_i$$

のような形で与えられた時、係数 $β1,β2,β3\beta_1, \beta_2, \beta_3$ の解釈はしばしば混乱を生む。

よくある誤解として、例えば $β3\beta_3$ は

$$yi=α3x3,i+ϵi y_i = \alpha_3 x_{3,i} + \epsilon_i$$

のような単回帰の係数 $α3\alpha_3$ に一致するというものがあるが、これは一般には誤りである。

実際には、$β3\beta_3$ は $x3x_3$ から $x1,x2x_1, x_2$ の影響を除去した残りの成分で $yiy_i$ に回帰する際の係数であり、「$x1,x2x_1, x_2$ を固定して $x3x_3$ だけを動かした際の変化率」の意味合いに近い。^[1]

このことを説明する際によく用いられるのが、Frisch-Waugh-Lovell(FWL)定理としても知られる残差回帰である。

残差回帰による係数の求め方とその解釈

いきなり残差回帰の一般論に入らず、具体的に冒頭の3変数の重回帰モデルを例に解説を行う。

$$yi=β0+β1x1,i+β2x2,i+β3x3,i+εi,(i=1,2,...,n) y_i = \beta_0 + \beta_1 x_{1,i} + \beta_2 x_{2,i} + \beta_3 x_{3, i} + \varepsilon_i, \quad (i=1, 2, ..., n)$$

$εi\varepsilon_i$ は誤差項である。
この時の最小二乗推定量 $β^=(β^0,β^1,β^2,β^3)⊤\hat{\boldsymbol{\beta}} = (\hat{\beta}_0, \hat{\beta}_1, \hat{\beta}_2, \hat{\beta}_3)^\top$ は、ストレートには以下のように求められる：^[2]

$$(β^0β^1β^2β^3)=(n∑ix1,i∑ix2,i∑ix3,i∑ix1,i∑ix1,i2∑ix1,ix2,i∑ix1,ix3,i∑ix2,i∑ix2,ix1,i∑ix2,i2∑ix2,ix3,i∑ix3,i∑ix3,ix1,i∑ix3,ix2,i∑ix3,i2)−1(∑iyi∑ix1,iyi∑ix2,iyi∑ix3,iyi) \begin{aligned} \begin{pmatrix} \hat{\beta}_0\\ \hat{\beta}_1\\ \hat{\beta}_2\\ \hat{\beta}_3 \end{pmatrix} &= \begin{pmatrix} n & \sum_i x_{1,i} & \sum_i x_{2,i} & \sum_i x_{3,i}\\ \sum_i x_{1,i} & \sum_i x_{1,i}^2 & \sum_i x_{1,i} x_{2,i} & \sum_i x_{1,i} x_{3,i}\\ \sum_i x_{2,i} & \sum_i x_{2,i} x_{1,i} & \sum_i x_{2,i}^2 & \sum_i x_{2,i} x_{3,i}\\ \sum_i x_{3,i} & \sum_i x_{3,i} x_{1,i} & \sum_i x_{3,i} x_{2,i} & \sum_i x_{3,i}^2 \end{pmatrix}^{-1} \begin{pmatrix} \sum_i y_i\\ \sum_i x_{1,i} y_i\\ \sum_i x_{2,i} y_i\\ \sum_i x_{3,i} y_i \end{pmatrix} \end{aligned}$$

一方で、次のように求めることもできる。

まず、$x3x_3$ を $x1,x2x_1, x_2$ で回帰することを考える：

$$x3,i=δ0+δ1x1,i+δ2x2,i+ϵi x_{3, i} = \delta_0 + \delta_1 x_{1, i} + \delta_2 x_{2,i} +\epsilon_i$$

$ϵi\epsilon_i$ は、$εi\varepsilon_i$ とは異なる誤差項である。
この時の係数 $δ0,δ1,δ2\delta_0, \delta_1, \delta_2$ の最小二乗推定量を $δ^0,δ^1,δ^2\hat{\delta}_0, \hat{\delta}_1, \hat{\delta}_2$ とすると、この線形回帰モデルにおける $x3,ix_{3,i}$ の残差は、

$$x~3,i=x3,i−δ^0−δ^1x1,i−δ^2x2,i \tilde{x}_{3, i} = x_{3,i} - \hat{\delta}_0 - \hat{\delta}_1 x_{1, i} - \hat{\delta}_2 x_{2,i}$$

と表すことができる。

ここでさらに、残差 $x~3,i\tilde{x}_{3, i}$ で $yiy_i$ に回帰する以下のようなモデルを考える（$ηi\eta_i$ は誤差項）：

$$yi=γ3x~3,i+ηi. y_i = \gamma_3 \tilde{x}_{3,i} + \eta_i.$$

この時の係数 $γ3\gamma_3$ の最小二乗推定量は、

$$γ^3=∑ix~3,iyi∑ix~3,i2=∑i(x3,i−δ^0−δ^1x1,i−δ^2x2,i)yi∑i(x3,i−δ^0−δ^1x1,i−δ^2x2,i)2 \begin{aligned} \hat{\gamma}_3 &= \frac{\sum_i \tilde{x}_{3, i} y_i}{\sum_i \tilde{x}_{3, i}^2}\\ &= \frac{\sum_i (x_{3,i} - \hat{\delta}_0 - \hat{\delta}_1 x_{1, i} - \hat{\delta}_2 x_{2,i}) y_i}{\sum_i (x_{3,i} - \hat{\delta}_0 - \hat{\delta}_1 x_{1, i}- \hat{\delta}_2 x_{2,i})^2}\\ \end{aligned}$$

のように求まるが、これが重回帰モデルにおける係数 $β3\beta_3$ の推定量と一致するのである（$β^3=γ^3\hat{\beta}_3 = \hat{\gamma}_3$）。

したがって、係数 $β3\beta_3$ は、説明変数 $x3,ix_{3,i}$ から残りの説明変数 $x1,i,x2,ix_{1,i}, x_{2,i}$ の影響を差し引いた成分（残差 $x~3,i\tilde{x}_{3,i}$）で $yiy_i$ に回帰した際の係数と解釈することができる。

一般の場合の残差回帰

以下では一般的な形で残差回帰を解説する。

設定

変数が $pp$ 個の線形回帰

$$yi=β1x1,i+β2x2,i+...+βpxp,i+εi,i=1,2,...,nE[εi]=0,E[εiεj]={σ2(i=j)0(i≠j) y_i = \beta_1 x_{1,i} + \beta_2 x_{2,i} + ... + \beta_p x_{p, i} + \varepsilon_i, \quad i=1, 2, ..., n\\ E[\varepsilon_i]=0, \quad E[\varepsilon_i \varepsilon_j] = \begin{cases} \sigma^2 \quad (i=j)\\ 0 \quad (i\ne j) \end{cases}$$

を、計画行列 $X∈Rn×p\boldsymbol{X} \in \mathbb{R}^{n \times p}$ を用いて以下のように表す：

$$Y=Xβ+ε(1) \boldsymbol{Y} = \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon} \tag{1}$$

ただし、

$$Y=(y1y2⋮yn),X=(x1,1x2,1...xp,1x1,2x2,2...xp,2⋮⋮x1,nx2,n...xp,n), \boldsymbol{Y} = \begin{pmatrix} y_1\\ y_2\\ \vdots\\ y_n \end{pmatrix}, \quad \boldsymbol{X} = \begin{pmatrix} x_{1,1} & x_{2,1} & ... & x_{p, 1}\\ x_{1,2} & x_{2,2} & ... & x_{p, 2}\\ \vdots & & & \vdots\\ x_{1,n} & x_{2,n} & ... & x_{p, n} \end{pmatrix},$$

$$β=(β1,β2,...,βp)⊤,ε=(ε1,ε2,...,εn)⊤. \boldsymbol{\beta} = (\beta_1, \beta_2, ..., \beta_p)^\top, \quad \boldsymbol{\varepsilon} = (\varepsilon_1, \varepsilon_2, ..., \varepsilon_n)^\top.$$

このとき、 係数 $β\boldsymbol{\beta}$ の最小二乗推定量を $β^\hat{\boldsymbol{\beta}}$ は

$$β^=(X⊤X)−1X⊤Y \hat{\boldsymbol{\beta}} = \left( \boldsymbol{X}^\top \boldsymbol{X} \right)^{-1} \boldsymbol{X}^\top \boldsymbol{Y}$$

のように求まり^[2:1]、これを用いて $Y\boldsymbol{Y}$ を以下のように表すことができる：

$$Y=Xβ^+ε^(2) \boldsymbol{Y} = \boldsymbol{X} \hat{\boldsymbol{\beta}} + \hat{\boldsymbol{\varepsilon}} \tag{2}$$

ただし、残差を $ε^=Y−Xβ^\hat{\boldsymbol{\varepsilon}} = \boldsymbol{Y} - \boldsymbol{X} \hat{\boldsymbol{\beta}}$ とした。

なお、以降の議論では $p≤np\le n$ かつ $rank X=p{\rm rank}\, \boldsymbol{X} = p$ で $X\boldsymbol{X}$ は full-rank であり、逆行列 $(X⊤X)−1(\boldsymbol{X}^\top \boldsymbol{X})^{-1}$ が存在するものとする。^[3]

Frisch-Waugh-Lovell (FWL) の定理

上記の設定のもと、$pp$ 個の説明変数（定数含む）を $p1p_1$ 個と $p−p1p - p_1$ 個の２グループに分割し、計画行列を $X=(X1,X2)\boldsymbol{X} = \left(\boldsymbol{X}_1, \boldsymbol{X}_2 \right)$ のように表す。
つまり、

$$X1=(x1,1x2,1...xp1,1x1,2x2,2...xp1,2⋮⋮x1,nx2,n...xp1,n),X2=(xp1+1,1xp1+2,1...xp,1xp1+1,2xp1+2,2...xp,2⋮⋮xp1+1,nxp1+2,n...xp,n). \boldsymbol{X}_1 = \begin{pmatrix} x_{1,1} & x_{2,1} & ... & x_{p_1, 1}\\ x_{1,2} & x_{2,2} & ... & x_{p_1, 2}\\ \vdots & & & \vdots\\ x_{1,n} & x_{2,n} & ... & x_{p_1, n} \end{pmatrix}, \quad \boldsymbol{X}_2 = \begin{pmatrix} x_{p_1 + 1,1} & x_{p_1 + 2,1} & ... & x_{p, 1}\\ x_{p_1 + 1,2} & x_{p_1 + 2,2} & ... & x_{p, 2}\\ \vdots & & & \vdots\\ x_{p_1 + 1,n} & x_{p_1 + 2,n} & ... & x_{p, n} \end{pmatrix}.$$

こうすることで、式(2) を以下のように書き換えることができる：

$$Y=X1β^1+X2β^2+ε^.(3) \boldsymbol{Y} = \boldsymbol{X}_1 \hat{\boldsymbol{\beta}}_1 + \boldsymbol{X}_2 \hat{\boldsymbol{\beta}}_2 + \hat{\boldsymbol{\varepsilon}} \tag{3}.$$

ただし、

$$β^=(β^1β^2),β^1,=(β^1β^2⋮β^p1),β^2=(β^p1+1β^p1+2⋮β^p) \hat{\boldsymbol{\beta}} = \begin{pmatrix} \hat{\boldsymbol{\beta}}_1\\ \hat{\boldsymbol{\beta}}_2 \end{pmatrix}, \quad \hat{\boldsymbol{\beta}}_1, = \begin{pmatrix} \hat{\beta}_1\\ \hat{\beta}_2\\ \vdots\\ \hat{\beta}_{p_1}\\ \end{pmatrix}, \quad \hat{\boldsymbol{\beta}}_2 = \begin{pmatrix} \hat{\beta}_{p_1 + 1}\\ \hat{\beta}_{p_1 + 2}\\ \vdots\\ \hat{\beta}_p\\ \end{pmatrix}$$

このとき、

$$M1=In−X1(X1⊤X1)−1X1⊤X~2=M1X2 \begin{aligned} \boldsymbol{M}_1 &= \boldsymbol{I}_n - \boldsymbol{X}_1 \left( \boldsymbol{X}_1^\top \boldsymbol{X}_1 \right)^{-1} \boldsymbol{X}_1^\top\\ \tilde{\boldsymbol{X}}_2 &= \boldsymbol{M}_1 \boldsymbol{X}_2 \end{aligned}$$

と置く。

すると、後半 $p−p1p - p_1$ 個の説明変数グループの係数推定量 $β^2\hat{\boldsymbol{\beta}}_2$ は、以下のように求めることができる：

$$β^2=(X~2⊤X~2)−1X~2⊤Y.(4) \hat{\boldsymbol{\beta}}_2 = (\tilde{\boldsymbol{X}}_2^\top \tilde{\boldsymbol{X}}_2)^{-1} \tilde{\boldsymbol{X}}_2^\top \boldsymbol{Y}. \tag{4}$$

このような係数推定量の求め方を残差回帰と呼ぶほか、上記の関係は Frisch-Waugh-Lovell (FWL) の定理としても知られる。

意味・解釈

$X~2\tilde{\boldsymbol{X}}_2$ は $X2\boldsymbol{X}_2$ を $X1\boldsymbol{X}_1$ で回帰を行った際の残差を意味する。
したがって、式(4) から $β^2\hat{\boldsymbol{\beta}}_2$ は $Y\boldsymbol{Y}$ を残差 $X~2\tilde{\boldsymbol{X}}_2$ で回帰した際の係数として解釈できる。

$X~2\tilde{\boldsymbol{X}}_2$ が $X1\boldsymbol{X}_1$ による回帰の残差であることは以下のようにして分かる。
回帰モデル

$$X2=X1δ+ϵ \boldsymbol{X}_2 = \boldsymbol{X}_1 \boldsymbol{\delta} + \boldsymbol{\epsilon}$$

において、係数 $δ\boldsymbol{\delta}$ の最小二乗推定量 $δ^\hat{\boldsymbol{\delta}}$ が

$$δ^=(X1⊤X1)−1X1⊤X2 \hat{\boldsymbol{\delta}} = \left(\boldsymbol{X}_1^\top \boldsymbol{X}_1\right)^{-1} \boldsymbol{X}_1^\top \boldsymbol{X}_2$$

のように求められるため、その際の残差は、

$$X2−X1δ^=X2−X1(X1⊤X1)−1X1⊤X2={In−X1(X1⊤X1)−1X1⊤}X2=M1X2=X~2 \begin{aligned} \boldsymbol{X}_2 - \boldsymbol{X}_1\hat{\boldsymbol{\delta}} &= \boldsymbol{X}_2 - \boldsymbol{X}_1\left(\boldsymbol{X}_1^\top \boldsymbol{X}_1\right)^{-1} \boldsymbol{X}_1^\top \boldsymbol{X}_2\\ &= \left\{ \boldsymbol{I}_n - \boldsymbol{X}_1 \left( \boldsymbol{X}_1^\top \boldsymbol{X}_1 \right)^{-1} \boldsymbol{X}_1^\top \right\} \boldsymbol{X}_2\\ &= \boldsymbol{M}_1 \boldsymbol{X}_2 = \tilde{\boldsymbol{X}}_2 \end{aligned}$$

となり、確かに $X~2\tilde{\boldsymbol{X}}_2$ と一致することがわかる。

なお、冒頭の例は、

$$X1=(1x1,1x2,11x1,2x2,2⋮⋮1x1,nx2,n),X2=(x3,1x3,2⋮x3,n) \boldsymbol{X}_1 = \begin{pmatrix} 1 & x_{1, 1} & x_{2, 1}\\ 1 & x_{1, 2} & x_{2, 2}\\ \vdots & \vdots\\ 1 & x_{1, n} & x_{2, n} \end{pmatrix}, \qquad \boldsymbol{X}_2= \begin{pmatrix} x_{3,1}\\ x_{3,2}\\ \vdots\\ x_{3,n} \end{pmatrix}$$

とした場合に相当する。

証明

以下では、式(4) が成り立つことを2通りの方法で証明する。

証明1: 正攻法

$Y\boldsymbol{Y}$ を $X1\boldsymbol{X}_1$ で回帰した際の係数 $δ^\hat{\boldsymbol{\delta}}$ は

$$δ^=(X1⊤X1)−1X1⊤Y \hat{\boldsymbol{\delta}} = \left(\boldsymbol{X}_1^\top \boldsymbol{X}_1 \right)^{-1} \boldsymbol{X}_1^\top \boldsymbol{Y}$$

のように表され、さらにその際の残差 $Y~=Y−X1δ^\tilde{\boldsymbol{Y}} = \boldsymbol{Y} - \boldsymbol{X}_1 \hat{\boldsymbol{\delta}}$ は、以下のように表すことができる：

$$Y~={In−X1(X1⊤X1)−1X1⊤}Y=M1Y \begin{aligned} \tilde{\boldsymbol{Y}} &= \left\{ \boldsymbol{I}_n - \boldsymbol{X}_1 \left( \boldsymbol{X}_1^\top \boldsymbol{X}_1 \right)^{-1} \boldsymbol{X}_1^\top \right\} \boldsymbol{Y}\\ &= \boldsymbol{M}_1 \boldsymbol{Y} \end{aligned}$$

ここで、式(3) を代入すると、

$$Y~=M1(X1β^1+X2β^2+ε^)=M1X2β^2+ε^∵M1X1=0, M1ε^=ε^ \begin{aligned} \tilde{\boldsymbol{Y}} &= \boldsymbol{M}_1 \left( \boldsymbol{X}_1 \hat{\boldsymbol{\beta}}_1 + \boldsymbol{X}_2 \hat{\boldsymbol{\beta}}_2 + \hat{\boldsymbol{\varepsilon}}\right)\\ &= \boldsymbol{M}_1 \boldsymbol{X}_2 \hat{\boldsymbol{\beta}}_2 + \hat{\boldsymbol{\varepsilon}} \qquad \because \boldsymbol{M}_1 \boldsymbol{X}_1 = \boldsymbol{0}, \, \boldsymbol{M}_1 \hat{\boldsymbol{\varepsilon}} = \hat{\boldsymbol{\varepsilon}} \end{aligned}$$

となり^[4][5]、さらに $M1X2=X~2\boldsymbol{M}_1 \boldsymbol{X}_2 = \tilde{\boldsymbol{X}}_2$ であったので、$Y~=X~2β^2+ε^\tilde{\boldsymbol{Y}} = \tilde{\boldsymbol{X}}_2 \hat{\boldsymbol{\beta}}_2 + \hat{\boldsymbol{\varepsilon}}$ が成り立っていることがわかる。

この両辺に左から $X~2⊤\tilde{\boldsymbol{X}}_2^\top$ をかけて $X~2⊤Y~=X~2⊤X~2β^2\tilde{\boldsymbol{X}}_2^\top\tilde{\boldsymbol{Y}} = \tilde{\boldsymbol{X}}_2^\top\tilde{\boldsymbol{X}}_2 \hat{\boldsymbol{\beta}}_2$ とし^[6]、さらに $(X~2⊤X~2)−1\left( \tilde{\boldsymbol{X}}_2^\top\tilde{\boldsymbol{X}}_2 \right)^{-1}$ を左からかけることで、

$$β^2=(X~2⊤X~2)−1X~2⊤Y~ \hat{\boldsymbol{\beta}}_2 = (\tilde{\boldsymbol{X}}_2^\top \tilde{\boldsymbol{X}}_2)^{-1} \tilde{\boldsymbol{X}}_2^\top \tilde{\boldsymbol{Y}}$$

が成立することがわかる。

ここで、行列 $M1\boldsymbol{M}_1$ は対称 $M1⊤=M1\boldsymbol{M}_1^\top = \boldsymbol{M}_1$ かつ冪等 $M12=M1\boldsymbol{M}_1^2 = \boldsymbol{M}_1$ なので、これを利用すると、

$$X~2⊤Y~=(M1X2)⊤M1Y=X2M1⊤M1Y=X2M1Y=(M1X2)⊤Y=X~2⊤Y \begin{aligned} \tilde{\boldsymbol{X}}_2^\top \tilde{\boldsymbol{Y}} &= \left( \boldsymbol{M}_1 \boldsymbol{X}_2 \right)^\top \boldsymbol{M}_1 \boldsymbol{Y}\\ &= \boldsymbol{X}_2 \boldsymbol{M}_1^\top \boldsymbol{M}_1 \boldsymbol{Y}\\ &= \boldsymbol{X}_2 \boldsymbol{M}_1 \boldsymbol{Y}\\ &= \left( \boldsymbol{M}_1 \boldsymbol{X}_2 \right)^\top \boldsymbol{Y}\\ &= \tilde{\boldsymbol{X}}_2^\top \boldsymbol{Y} \end{aligned}$$

が成り立ち、したがって

$$β^2=(X~2⊤X~2)−1X~2⊤Y \hat{\boldsymbol{\beta}}_2 = (\tilde{\boldsymbol{X}}_2^\top \tilde{\boldsymbol{X}}_2)^{-1} \tilde{\boldsymbol{X}}_2^\top \boldsymbol{Y}$$

が示される。

証明2: ブロック行列を用いた方法

推定量 $β^\hat{\boldsymbol{\beta}}$ は $β^=(X⊤X)−1X⊤Y\hat{\boldsymbol{\beta}} = \left( \boldsymbol{X}^\top \boldsymbol{X}\right)^{-1} \boldsymbol{X}^\top \boldsymbol{Y}$ のように得られたので、以下が成り立つ：

$$(β^1β^2)={(X1⊤X2⊤)(X1X2)}−1(X1⊤X2⊤)Y=(X1⊤X1X1⊤X2X2⊤X1X2⊤X2)−1(X1⊤YX2⊤Y) \begin{aligned} \begin{pmatrix} \hat{\boldsymbol{\beta}}_1\\ \hat{\boldsymbol{\beta}}_2 \end{pmatrix} &= \left\{ \begin{pmatrix} \boldsymbol{X}_1^\top\\ \boldsymbol{X}_2^\top \end{pmatrix} \begin{pmatrix} \boldsymbol{X}_1 & \boldsymbol{X}_2 \end{pmatrix} \right\}^{-1} \begin{pmatrix} \boldsymbol{X}_1^\top\\ \boldsymbol{X}_2^\top \end{pmatrix} \boldsymbol{Y}\\ &= \begin{pmatrix} \boldsymbol{X}_1^\top \boldsymbol{X}_1 & \boldsymbol{X}_1^\top \boldsymbol{X}_2\\ \boldsymbol{X}_2^\top \boldsymbol{X}_1 & \boldsymbol{X}_2^\top \boldsymbol{X}_2 \end{pmatrix}^{-1} \begin{pmatrix} \boldsymbol{X}_1^\top \boldsymbol{Y}\\ \boldsymbol{X}_2^\top \boldsymbol{Y} \end{pmatrix} \end{aligned}$$

ここで、ブロック化した行列の逆行列について、以下の関係式が成り立つ：

$$(TUVW)−1=(T−1+T−1UQ−1VT−1−T−1UQ−1−Q−1VT−1Q−1) \begin{pmatrix} \boldsymbol{T} & \boldsymbol{U}\\ \boldsymbol{V} & \boldsymbol{W} \end{pmatrix}^{-1} = \begin{pmatrix} \boldsymbol{T}^{-1} + \boldsymbol{T}^{-1} \boldsymbol{U} \boldsymbol{Q}^{-1} \boldsymbol{V} \boldsymbol{T}^{-1} & - \boldsymbol{T}^{-1}\boldsymbol{U}\boldsymbol{Q}^{-1}\\ - \boldsymbol{Q}^{-1} \boldsymbol{V} \boldsymbol{T}^{-1} & \boldsymbol{Q}^{-1} \end{pmatrix}$$

ただし、$T,W\boldsymbol{T}, \boldsymbol{W}$ は正方行列、$Q=W−VT−1U\boldsymbol{Q} = \boldsymbol{W} - \boldsymbol{V} \boldsymbol{T}^{-1} \boldsymbol{U}$ と定義され、さらに $T\boldsymbol{T}$ と $Q\boldsymbol{Q}$ は逆行列を持つものとする。

https://zenn.dev/link/comments/4f90e72bf8b9b4

これを用いると、

$$(X1⊤X1X1⊤X2X2⊤X1X2⊤X2)−1=((X1⊤X1)−1+(X1⊤X1)−1X1⊤X2Q−1X2⊤X1(X1⊤X1)−1−(X1⊤X1)−1X1⊤X2Q−1−Q−1X2⊤X1(X1⊤X1)−1Q−1) \begin{pmatrix} \boldsymbol{X}_1^\top \boldsymbol{X}_1 & \boldsymbol{X}_1^\top \boldsymbol{X}_2\\ \boldsymbol{X}_2^\top \boldsymbol{X}_1 & \boldsymbol{X}_2^\top \boldsymbol{X}_2 \end{pmatrix}^{-1} = \begin{pmatrix} \left(\boldsymbol{X}_1^\top \boldsymbol{X}_1\right)^{-1} + \left(\boldsymbol{X}_1^\top \boldsymbol{X}_1\right)^{-1} \boldsymbol{X}_1^\top \boldsymbol{X}_2 \boldsymbol{Q}^{-1} \boldsymbol{X}_2^\top \boldsymbol{X}_1 \left(\boldsymbol{X}_1^\top \boldsymbol{X}_1\right)^{-1} & - \left(\boldsymbol{X}_1^\top \boldsymbol{X}_1\right)^{-1} \boldsymbol{X}_1^\top \boldsymbol{X}_2 \boldsymbol{Q}^{-1} \\ - \boldsymbol{Q}^{-1} \boldsymbol{X}_2^\top \boldsymbol{X}_1\left(\boldsymbol{X}_1^\top \boldsymbol{X}_1\right)^{-1} & \boldsymbol{Q}^{-1} \end{pmatrix}$$

ただし

$$Q=X2⊤X2−X2⊤X1(X1⊤X1)−1X1⊤X2=X2⊤{In−X1(X1⊤X1)−1X1⊤}X2=X2⊤M1X2=X2⊤M1⊤M1X2=X~2⊤X~2 \begin{aligned} \boldsymbol{Q} &= \boldsymbol{X}_2^\top \boldsymbol{X}_2 - \boldsymbol{X}_2^\top \boldsymbol{X}_1 \left( \boldsymbol{X}_1^\top \boldsymbol{X}_1 \right)^{-1} \boldsymbol{X}_1^\top \boldsymbol{X}_2\\ &= \boldsymbol{X}_2^\top \left\{ \boldsymbol{I}_n - \boldsymbol{X}_1 \left( \boldsymbol{X}_1^\top \boldsymbol{X}_1 \right)^{-1} \boldsymbol{X}_1^\top \right\} \boldsymbol{X}_2\\ &= \boldsymbol{X}_2^\top \boldsymbol{M}_1 \boldsymbol{X}_2\\ &= \boldsymbol{X}_2^\top \boldsymbol{M}_1^\top \boldsymbol{M}_1 \boldsymbol{X}_2\\ &= \tilde{\boldsymbol{X}}_2^\top\tilde{\boldsymbol{X}}_2 \end{aligned}$$

のようになる。

したがって、

$$β^2=−Q−1X2⊤X1(X1⊤X1)−1X1⊤Y+Q−1X2⊤Y=Q−1X2⊤{In−X1(X1⊤X1)−1X1⊤}Y=Q−1X2⊤M1Y=(X~2⊤X~2)−1X~2Y \begin{aligned} \hat{\boldsymbol{\beta}}_2 &= - \boldsymbol{Q}^{-1} \boldsymbol{X}_2^\top \boldsymbol{X}_1\left(\boldsymbol{X}_1^\top \boldsymbol{X}_1\right)^{-1}\boldsymbol{X}_1^\top \boldsymbol{Y} + \boldsymbol{Q}^{-1}\boldsymbol{X}_2^\top \boldsymbol{Y}\\ &= \boldsymbol{Q}^{-1} \boldsymbol{X}_2^\top \left\{ \boldsymbol{I}_n - \boldsymbol{X}_1\left(\boldsymbol{X}_1^\top \boldsymbol{X}_1\right)^{-1}\boldsymbol{X}_1^\top \right\} \boldsymbol{Y}\\ &= \boldsymbol{Q}^{-1} \boldsymbol{X}_2^\top \boldsymbol{M}_1 \boldsymbol{Y}\\ &= \left( \tilde{\boldsymbol{X}}_2^\top\tilde{\boldsymbol{X}}_2 \right)^{-1}\tilde{\boldsymbol{X}}_2 \boldsymbol{Y} \end{aligned}$$

が示される。

証明2の補足

なお、同様に $β^1\hat{\boldsymbol{\beta}}_1$ についても計算を行うと、

$$β^1=(X1⊤X1)−1X1⊤Y−(X1⊤X1)−1X1⊤X2(X~2⊤X~2)−1X~2Y=(X1⊤X1)−1X1⊤Y−(X1⊤X1)−1X1⊤X2β^2 \begin{aligned} \hat{\boldsymbol{\beta}}_1 &= \left(\boldsymbol{X}_1^\top \boldsymbol{X}_1\right)^{-1} \boldsymbol{X}_1^\top \boldsymbol{Y} - \left(\boldsymbol{X}_1^\top \boldsymbol{X}_1\right)^{-1} \boldsymbol{X}_1^\top \boldsymbol{X}_2 \left( \tilde{\boldsymbol{X}}_2^\top\tilde{\boldsymbol{X}}_2 \right)^{-1}\tilde{\boldsymbol{X}}_2 \boldsymbol{Y}\\ &= \left(\boldsymbol{X}_1^\top \boldsymbol{X}_1\right)^{-1} \boldsymbol{X}_1^\top \boldsymbol{Y} - \left(\boldsymbol{X}_1^\top \boldsymbol{X}_1\right)^{-1} \boldsymbol{X}_1^\top \boldsymbol{X}_2 \hat{\boldsymbol{\beta}}_2 \end{aligned}$$

が得られる。

特に $X1=(1,1,...,1)⊤\boldsymbol{X}_1 = (1, 1, ..., 1)^\top$ の場合、

$$β^1=1n∑i=1nyi−β^21n∑i=1nx2,i−β^31n∑i=1nx3,i−...−β^p1n∑i=1nxp,i=yˉ−β^2xˉ2−β^3xˉ3−...−β^pxˉp \begin{aligned} \hat{\beta}_1 &= \frac{1}{n}\sum_{i=1}^n y_i - \hat{\beta}_2 \frac{1}{n} \sum_{i=1}^n x_{2,i} - \hat{\beta}_3 \frac{1}{n} \sum_{i=1}^n x_{3,i} - ... - \hat{\beta}_p \frac{1}{n} \sum_{i=1}^n x_{p,i}\\ &= \bar{y} - \hat{\beta}_2 \bar{x}_2 - \hat{\beta}_3 \bar{x}_3 - ... - \hat{\beta}_p \bar{x}_p \end{aligned}$$

のように、よく見た定数項の求め方になる ($yˉ=1n∑iyi\bar{y}= \frac{1}{n}\sum_i y_i$, $xˉk=1n∑ixk,i\bar{x}_k = \frac{1}{n}\sum_i x_{k,i}$)。
なお、このとき残差 $X~2\tilde{\boldsymbol{X}}_2$ は

$$X~2=(x2,1−xˉ2x3,1−xˉ3...xp,1−xˉpx2,2−xˉ2x3,2−xˉ3...xp,2−xˉp⋮⋮x2,n−xˉ2x3,n−xˉ3...xp,n−xˉp) \tilde{\boldsymbol{X}}_2 = \begin{pmatrix} x_{2,1} - \bar{x}_2 & x_{3,1} - \bar{x}_3 & ... & x_{p,1} - \bar{x}_p\\ x_{2,2} - \bar{x}_2 & x_{3,2} - \bar{x}_3 & ... & x_{p,2} - \bar{x}_p\\ \vdots & & & \vdots\\ x_{2,n} - \bar{x}_2 & x_{3,n} - \bar{x}_3 & ... & x_{p,n} - \bar{x}_p\\ \end{pmatrix}$$

のように各説明変数の平均値を差し引いたものになる。

参考文献

• 高橋将宜、「統計的因果推論の理論と実装」（共立出版、2022）

https://amzn.asia/d/fmchigb

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

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

• 西山慶彦・新谷元嗣・川口大司・奥井亮、「計量経済学(New Liberal Arts Selection)」（有斐閣、2019）

https://amzn.asia/d/hbZoFo2

脚注

1. 参考文献に挙げた 高橋(2022) の Ch.6 では、統計的因果推論の文脈から重回帰について直感的にわかりやすい解説がなされている。

2. https://zenn.dev/link/comments/db9e98ef7a0301 参照

3. https://zenn.dev/tatamiya/articles/ca2ab1f8d3f069b78a82, https://zenn.dev/link/comments/116bc995370fff 参照

4. $M1X1=0\boldsymbol{M}_1 \boldsymbol{X}_1 = \boldsymbol{0}$ は、$M1=In−X1(X1⊤X1)−1X1⊤\boldsymbol{M}_1 = \boldsymbol{I}_n - \boldsymbol{X}_1 ( \boldsymbol{X}_1^\top \boldsymbol{X}_1)^{-1} \boldsymbol{X}_1^\top$ を代入すればすぐ示せる。

5. $M1ε^=ε^\boldsymbol{M}_1 \hat{\boldsymbol{\varepsilon}} = \hat{\boldsymbol{\varepsilon}}$ は、$X⊤ε^=0\boldsymbol{X}^\top \hat{\boldsymbol{\varepsilon}} = \boldsymbol{0}$ であることに注目して $X1⊤ε^=0\boldsymbol{X}_1^\top \hat{\boldsymbol{\varepsilon}} = \boldsymbol{0}$ を用いると示せる。

6. $X~2⊤ε^=0\tilde{\boldsymbol{X}}_2^\top\hat{\boldsymbol{\varepsilon}}=\boldsymbol{0}$ になる。なぜなら、$X~2⊤ε^=X2⊤M1ε^=X2⊤ε^=0\tilde{\boldsymbol{X}}_2^\top\hat{\boldsymbol{\varepsilon}}=\boldsymbol{X}_2^\top \boldsymbol{M}_1 \hat{\boldsymbol{\varepsilon}} = \boldsymbol{X}_2^\top \hat{\boldsymbol{\varepsilon}} = \boldsymbol{0}$