Step2. 中級編
7. 重回帰分析

7-1. 重回帰分析1

複数の説明変数を含むデータがある場合、これらを同時に扱う方法として重回帰モデルを考えます。重回帰モデルでは、目的変数 $y$ が複数（ここでは $p$ 個）の説明変数 $x_{k}$ の線形結合によって表されると仮定します。

まず、モデル全体を次のように表します。

$\displaystyle y=\beta_{0}+\beta_{1}x_{1}+\cdots+\beta_{p}x_{p}+\varepsilon$

ここで、 $\beta_{k}$ は重回帰モデルにおける回帰係数（偏回帰係数）を表します。他の説明変数を一定に保ったときに、説明変数 $x_k$ が1単位変化した時の目的変数 $y$ の平均的な変化量を表します。 $\varepsilon$ は誤差項であり、説明変数では説明しきれないばらつきを表します。誤差項は平均0、分散 $\sigma^2$ の正規分布に従うと仮定します。

$\displaystyle \varepsilon \sim N(0, \sigma^{2})$

■重回帰モデル

$p$ 個の説明変数からなる $n$ 個の観測値がある場合、重回帰モデルは各観測値ごとに次のように表すことができます。

$\displaystyle y_i=\beta_{0}+\beta_{1}x_{1i}+\cdots+\beta_{p}x_{pi}+\varepsilon_{i} ~~~~~ (i=1, 2, \cdots, n)$

ここで、

$y_i$ は $i$ 番目の観測値における目的変数
$x_{ki}$ は $i$ 番目の観測値における $k$ 番目の説明変数
$\varepsilon_{i}$ は $i$ 番目の観測値における誤差項

を表します。誤差項 $\varepsilon_{i}$ は互いに独立で、正規分布 $N(0, \sigma^{2})$ に従うものとします。

$\displaystyle \varepsilon_{i} \sim N(0, \sigma^{2})$

この重回帰モデルにおける回帰係数 $\beta_{k} (k=1, 2, \cdots, p)$ を求めるための方法の1つが最小二乗法です。

■最小二乗法

最小二乗法とは観測値と予測値の誤差の二乗和が最小になるように、予測値の未知のパラメーターを求める方法のことです。すなわち、重回帰モデルの回帰係数 $\beta_k$ は直接観測することができないため、観測値にできるだけ近い予測値を与えるように回帰係数を推定する必要があります。回帰係数の推定値は一般に $\hat{\beta}_k$ と表します。

最小二乗法では、 $i$ 番目の観測値 $y_i$ に対する予測値 $\hat{y}_i$ を次のように表すことができます。

$\displaystyle \hat{y}_i=\hat{\beta}_{0}+\hat{\beta}_{1}x_{1i}+\cdots+\hat{\beta}_{p}x_{pi}$

このとき、実際の観測値 $y_i$ と予測値 $\hat{y}_i$ との差を残差と呼び、次のように定義します。

$\displaystyle e_i = y_{i}-\hat{y}_i=y_{i}-(\hat{\beta}_{0}+\hat{\beta}_{1}x_{1i}+\cdots+\hat{\beta}_{p}x_{pi})$

最小二乗法を用いて、残差の平方和（残差平方和） $S_e$ を最小にする $\beta_{k} (k=1, 2, \cdots, p)$ を計算します。

$\displaystyle S_e = \sum_{i=1}^{n} e_{i}^2 = \sum_{i=1}^{n} \{y_{i}-\hat{y}_i \}^2 = \{y_{i}-(\hat{\beta}_{0}+\hat{\beta}_{1}x_{1i}+\cdots+\hat{\beta}_{p}x_{pi}) \}^2$

具体的には $S_{e}$ を $\beta_{k} (k=0, 2, \cdots, p)$ のそれぞれで偏微分して0とおいた式を使います。

$\displaystyle \frac{\partial S_e}{\partial \beta_k}=0 ~~~~~ (k=1, 2, \cdots, p)$

実際に計算すると、次のような連立方程式が得られます。この連立方程式を正規方程式と呼びます。この正規方程式を解くと、各回帰係数の推定値 $\hat{\beta}_k$ が算出できます。このとき得られる回帰係数の推定値 $\hat{\beta}_k$ を最小二乗推定量（least squares estimator）と呼びます。

$\displaystyle \frac{\partial S_e}{\partial \beta_0}=-2\sum_{i=1}^{n} \{y_{i}-(\hat{\beta}_{0}+\hat{\beta}_{1}x_{1i}+\cdots+\hat{\beta}_{p}x_{pi}) \}=0 ~~~~~ (k=0)$

$\displaystyle \frac{\partial S_e}{\partial \beta_k}=-2\sum_{i=1}^{n} x_{ki} \{y_{i}-(\hat{\beta}_{0}+\hat{\beta}_{1}x_{1i}+\cdots+\hat{\beta}_{p}x_{pi}) \}=0 ~~~~~ (k=1, 2, \cdots, p)$

■行列表現による重回帰モデル

説明変数が複数ある場合、行列を使うと式が見やすくなります。まず各記号を次のように定義します。

目的変数

$\mathbf{y}= \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix}$

回帰係数

$\boldsymbol{\beta}= \begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_p \end{pmatrix}$

デザイン行列

$\mathbf{X}= \begin{pmatrix} 1 & x_{11} & x_{21} & \cdots & x_{p1} \\ 1 & x_{12} & x_{22} & \cdots & x_{p2} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{1n} & x_{2n} & \cdots & x_{pn} \end{pmatrix}$

誤差項

$\boldsymbol{\varepsilon}= \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\ \varepsilon_n \end{pmatrix}$

これらを用いると、

重回帰モデル

$\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}$

予測値ベクトル（ただし回帰係数の推定値を $\hat{\beta}$ とします）

$\hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol{\beta}}$

残差ベクトル

$\mathbf{e} = \mathbf{y} - \hat{\mathbf{y}} = \mathbf{y} - \mathbf{X}\hat{\boldsymbol{\beta}}$

と表すことができます。残差平方和は次のようになります。

$S_e = \sum_{i=1}^{n} e_i^2 = \mathbf{e}^\mathsf{T}\mathbf{e} = (\mathbf{y}-\mathbf{X}\boldsymbol{\beta})^\mathsf{T} (\mathbf{y}-\mathbf{X}\boldsymbol{\beta})$

ここで $S_e$ を $\beta$ で微分して0とおくと、正規方程式が得られます。

$\mathbf{X}^\mathsf{T}\mathbf{X}\boldsymbol{\beta}=\mathbf{X}^\mathsf{T}\mathbf{y}$

さらに、 $\mathbf{X}^\mathsf{T}\mathbf{X}$ が逆行列をもつとき，正規方程式を解くことで、重回帰モデルにおける回帰係数の最小二乗推定量である $\hat{\boldsymbol{\beta}}$ が算出できます。

$\hat{\boldsymbol{\beta}} = (\mathbf{X}^\mathsf{T}\mathbf{X})^{-1}\mathbf{X}^\mathsf{T}\mathbf{y}$

正規方程式 $\mathbf{X}^\mathsf{T}\mathbf{X}\boldsymbol{\beta}=\mathbf{X}^\mathsf{T}\mathbf{y}$ は、残差ベクトル $\mathbf{e} = \mathbf{y} - \hat{\mathbf{y}} = \mathbf{y} - \mathbf{X}\hat{\boldsymbol{\beta}}$ が $\mathbf{X}$ の各列（切片列を含む）と直交することを意味します。

■重回帰モデルにおける仮定

これまでに示した最小二乗解が統計的に意味をもつためには、重回帰モデルに対していくつかの仮定を置く必要があります。以下に基本的な仮定をまとめます。

誤差項 $\varepsilon_i$ の期待値は $E(\varepsilon_i)=0$
誤差項の分散は一定である（等分散性）
誤差項は互いに独立である
誤差項は正規分布 $N(0,\sigma^2)$ に従う
説明変数間に完全な線形従属が存在しない（ $\mathbf{X}^\mathsf{T}\mathbf{X}$ が逆行列をもつ）

上記の仮定のもとで，最小二乗推定量 $\hat{\beta}_k$ は不偏であり，分散が最小となります（Gauss–Markov 定理）。

■決定係数

目的変数 $y$ の変動である全変動の平方和 $S_T$ は次のように計算できます。 $\bar{y}$ は $y$ 全体の平均です。

$\displaystyle S_T = \sum_{i=1}^n (y_i-\bar{y})^2$

また、残差平方和 $S_e$ は次のように計算できます。

$\displaystyle S_e = \sum_{i=1}^n (y_i-\hat{y}_i)^2$

これらを用いると、決定係数 $R^2$ は次の式で定義されます。

$\displaystyle R^2 = 1 - \frac{S_e}{S_T}$

$R^2$ は目的変数の変動のうち，回帰モデルによって説明される割合を表します。説明変数の個数が増えると $R^2$ は必ず増加するため，モデルの複雑さを考慮した指標として自由度調整済み決定係数 $\bar{R}^2$ を用います。

$\displaystyle \bar{R}^2 = 1 - \frac{S_e/(n-p-1)}{S_T/(n-1)}$

ここまでの説明を行列で表してみます。行列で扱うために、平均からの偏差ベクトルを定義します。目的変数ベクトルを平均で中心化したものを

$\tilde{\mathbf{y}} = \mathbf{y} - \bar{y}\mathbf{1}$

とします。ここで $\mathbf{1}$ はすべての要素が1の $n$ 次元ベクトルです。すると、全変動の平方和は次のように計算できます。

$S_T = \tilde{\mathbf{y}}^\mathsf{T}\tilde{\mathbf{y}}$

最小二乗解による予測値 $\hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol{\beta}}$ を射影行列で書くと

$\hat{\mathbf{y}} = \mathbf{P}\mathbf{y}, \qquad \mathbf{P} = \mathbf{X}(\mathbf{X}^\mathsf{T}\mathbf{X})^{-1}\mathbf{X}^\mathsf{T}$

ここで $\mathbf{P}$ は $\mathbf{X}$ の列空間への直交射影行列です。残差ベクトルは $\mathbf{e} = \mathbf{y} - \hat{\mathbf{y}} = \mathbf{y}^\mathsf{T}(\mathbf{I}-\mathbf{P})\mathbf{y}$ であることから、残差の平方和は次のように計算できます。