重回帰モデル

---

定義

$p$ 個の説明変数 $x_1,x_2,\dots ,x_p$ で $y$ を推定する次のモデル：

$$y_i={\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots {\beta }_p{x}_{ip}+{ϵ}_i$$

を 重回帰モデル といいます。ここで、$i$ はサンプルのインデックス（$i=1,2,\dots ,n$ とする）を表し、${ϵ}_i$ は ${ϵ}_i\sim N(0,{\sigma }^2)$ をみたす独立同分布からの誤差項であるとします。これにより、所与として $i$ 番目のサンプル $x_{i1},x_{i2},\dots ,x_{ip}$ を与えた時の $y_i$ は、$y_i\sim N({\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots {\beta }_p{x}_{ip},{\sigma }^2)$ なる確率変数であると考えることができます。

---

行列表現

サンプルをすべてまとめて表現したい場合、次の行列表現を考えます。

$$\mathbf{y}:=(y_1,\dots ,y_n{)}^{\top },X:=\left(\begin{array}{cccc}1& x_{11}& \cdots & x_{1p}\\ 1& x_{21}& \cdots & x_{2p}\\ ⋮& ⋮& \ddots & ⋮\\ 1& x_{n1}& \cdots & x_{np}\end{array}\right)$$
$$\mathbf{\beta }:=({\beta }_0,{\beta }_1,\dots ,{\beta }_p{)}^{\top },\mathbf{ϵ}:=({ϵ}_1,\dots ,{ϵ}_n{)}^{\top }$$

すると、

$$\mathbf{y}=X\mathbf{\beta }+\mathbf{ϵ}$$

と表すことができます。

最小二乗法

---

重回帰モデルの 回帰係数 ${\beta }_0,{\beta }_1,\dots ,{\beta }_p$ の推定量 ${\hat{\beta }}_0,{\hat{\beta }}_1,\dots ,{\hat{\beta }}_p$ を求める方法として、一般的なものに 最小二乗法 があります。

残差・残差平方和

最小二乗法の要点は、もとのデータ $y_i$ と回帰係数の推定量 ${\hat{\beta }}_0,{\hat{\beta }}_1,\dots ,{\hat{\beta }}_p$ を使って得られた $y_i$ の推定値：

$${\hat{y}}_i={\hat{\beta }}_0+{\hat{\beta }}_1{x}_{i1}+\cdots +{\hat{\beta }}_p{x}_{ip}$$

の差分 $y_i-{\hat{y}}_i$ をできるだけ小さくする、という考え方にあります。この量 $y_i-{\hat{y}}_i$ を $i$ 番目の観測値の 残差 といいます。また、残差の2乗和：

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

を 残差平方和 といいます。

最小二乗推定量

残差平方和 $S_e$ を最小化するような推定量 ${\beta }_k$ として、最小二乗推定量 を次の条件のもと定めます：

$$\frac{\partial }{\partial {\beta }_k}\sum _{i=1}^n[y_i-({\beta }_0+{\beta }_1{x}_{i1}+\cdots +{\beta }_p{x}_{ip}){]}^2=0$$

---

行列表現

ある回帰係数 $\mathbf{\beta }$を用いて計算される残差 $\mathbf{y}-X\mathbf{\beta }$ に対し、残差平方和はその内積 $||\mathbf{y}-X\mathbf{\beta }|{|}^2$ と与えることができます。これを $\mathbf{\beta }$ で微分することにより、

$$\frac{\partial }{\partial \mathbf{\beta }}||\mathbf{y}-X\mathbf{\beta }|{|}^2=-2X^{\top }\mathbf{y}+2X^{\top }X\mathbf{\beta }=0$$

から、これをみたす $\mathbf{\beta }$ を $\hat{\mathbf{\beta }}$ として、

$$X^{\top }X\hat{\mathbf{\beta }}=X^{\top }\mathbf{y}⟺\hat{\mathbf{\beta }}=(X^{\top }X{)}^{-1}{X}^{\top }\mathbf{y}$$

を得ます。これを 正規方程式 といいます。

平方和の分解・決定係数

各平方和（変動）を次のように定義します：

	総平方和 $S_T$	回帰平方和 $S_R$	残差平方和 $S_e$
定義	$\sum _{i=1}^n(y_i-\bar{y}{)}^2$	$\sum _{i=1}^n({\hat{y}}_i-\bar{y}{)}^2$	$\sum _{i=1}^n(y_i-{\hat{y}}_i{)}^2$
自由度	$n-1$	$p$	$n-p-1$

各平方和の間には、次の分解：

$$S_T=S_R+S_e$$

が成り立ちます。これらを $n-1$ で除せば、$S_T$ は $y$ の分散、$S_R$ は $\hat{y}$ の分散に相当します。

これらを用いて、

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

という量を定義すると、これは $0\le R^2\le 1$ を満たし、これを 決定係数 といいます。上式から、決定係数は $1$ に近いほど $S_e$ が小さく、よい当てはまりを示すモデルであることを表します。重回帰モデルのように説明変数が多い場合には、

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

として各変動の自由度で除した 自由度調整済み決定係数 を用いる場合が多いです。

最小二乗推定量の性質

期待値と分散

$\hat{\mathbf{\beta }}$ は不偏性をもちます：

$$E[\hat{\mathbf{\beta }}]=\mathbf{\beta }$$

しかも、$\hat{\mathbf{\beta }}$ は最小分散であり、以下を満たします：

$$V[\hat{\mathbf{\beta }}]={\sigma }^2(X^{\top }X{)}^{-1}$$

ここで ${\sigma }^2$ は誤差 ${ϵ}_i$ が独立同一にしたがう正規分布の分散ですが、この推定量 ${\hat{\sigma }}^2$ としては、

$${\hat{\sigma }}^2=\frac{S_e}{n-p-1}$$

が不偏推定量となります（$E[{\hat{\sigma }}^2]={\sigma }^2$）。結局、$\hat{\mathbf{\beta }}$ がしたがう確率分布として、

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

なる $p+1$ 次元正規分布を得ます。

---

最尤推定量との関係

誤差 ${ϵ}_i$ が独立同一に正規分布 $N(0,{\sigma }^2)$ にしたがうという仮定のもとでは、

$$y_i\sim N({\beta }_0+{\beta }_1{x}_{i1}+\cdots {\beta }_p{x}_{ip},{\sigma }^2)$$

ですから、データ $\{y_i\}$ を与えたもとでの尤度関数 $L({\beta }_0,{\beta }_1,\dots ,{\beta }_p)$ は次のように与えることができます：

$$L({\beta }_0,{\beta }_1,\dots ,{\beta }_p)=\prod _{i=1}^n\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left[-\frac{[y_i-({\beta }_0+{\beta }_1{x}_{i1}+\cdots {\beta }_p{x}_{ip}){]}^2}{2{\sigma }^2}\right]$$

これより、対数尤度関数は

$$\ln L({\beta }_0,{\beta }_1,\dots ,{\beta }_p)=\sum _{i=1}^n\left(-\frac{1}{2}\ln (2\pi {\sigma }^2)-\frac{[y_i-({\beta }_0+{\beta }_1{x}_{i1}+\cdots {\beta }_p{x}_{ip}){]}^2}{2{\sigma }^2}\right)$$

となりますから、$\partial \ln L/\partial {\beta }_k=0$ をみたすようにして決められる最尤推定量 ${\hat{\beta }}_k$ は残差平方和を最小化することに相当し、最小二乗推定量に等しくなります。

回帰の妥当性の検討

---

回帰係数の有意性検定

回帰係数の推定量 $\hat{\mathbf{\beta }}$ の効果を検定することを考えます。いま、すべての回帰係数について、ゼロか非ゼロを検定するための仮説を次のように設定します：

$$H_0:{\beta }_1={\beta }_2=\cdots ={\beta }_p=0$$
$$H_1:{\beta }_k\ne 0\text{s.t.}{}^{\exists }k\in \{1,2,\dots ,p\}$$

このとき、次のような $F$ 値を検定統計量として用いることで、この仮説を検定することができます：

$$F=\frac{S_R/p}{S_e/(n-p-1)}\sim F(p,n-p-1)$$

レバレッジ

最小二乗推定量 $\hat{\mathbf{\beta }}=(X^{\top }X{)}^{-1}{X}^{\top }\mathbf{y}$ に対し、次の行列：

$$H:=X(X^{\top }X{)}^{-1}{X}^{\top }$$

は $\hat{\mathbf{y}}=X\hat{\mathbf{\beta }}=H\mathbf{y}$ をみたします。このような行列 $H$ において、$i$ 番目の対角要素 $h_{ii}$ を $i$ 番目の観測値における レバレッジ といい、観測値 $y_i$ が回帰係数 ${\beta }_1,{\beta }_2,\dots ,{\beta }_p$ へ及ぼす影響の大きさを表します。

レバレッジの値は $y_i$ を外れ値とみなすかどうかの指標として利用できます。

ペナルティ付き回帰モデル

---

以上のような重回帰モデルの推定は、訓練データ $X,\mathbf{y}$ を与えたもとでモデルを構成する（最適な回帰係数 $\mathbf{\beta }$ を決める）教師あり学習の手法と考えることができます。ここで、最小二乗推定の際に最小化する量：

$$E(\mathbf{\beta })=\frac{1}{n}||\mathbf{y}-X\mathbf{\beta }|{|}^2(=S_e/n)$$

を、機械学習の文脈で 損失関数^[1] とよびます。

重回帰モデルでは、一般に説明変数の数 $p$ が増えると、用いた訓練データへの適合が大きすぎるために回帰の精度が落ちるという 過学習 の問題が生じます。この問題に対処するために、損失関数にデータへの適合を押さえるための 正則化項（ペナルティ） を加えたモデルが存在します。

Ridge回帰

以下の関数：

$$\frac{1}{n}||\mathbf{y}-X\mathbf{\beta }|{|}^2+\lambda (||\mathbf{\beta }|{|}_2{)}^2$$

を最小化する回帰モデルは Ridge回帰 とよばれます。ここで $\lambda$ は正則化の強さを決めるハイパーパラメータです。Ridge回帰では正則化のために $L_2$-ノルム^[2] $||\mathbf{\beta }|{|}_2$ を用いており、このような形式を $L_2$-正則化 とよびます。

このとき、正則化項は $\lambda {\mathbf{\beta }}^{\top }\mathbf{\beta }$ と書きなおせますから、最小化の条件は

$$-2X^{\top }\mathbf{y}+2X^{\top }X\hat{\mathbf{\beta }}+2\lambda \hat{\mathbf{\beta }}=0⟺\hat{\mathbf{\beta }}=(X^{\top }X+\lambda I{)}^{-1}{X}^{\top }\mathbf{y}$$

となります（$I$ は単位行列）。Ridge回帰モデルは効果の大きい回帰係数にペナルティを課すことで過学習を防ぎます。

LASSO回帰

$L_1$-ノルムを正則化に用いた以下の関数：

$$\frac{1}{n}||\mathbf{y}-X\mathbf{\beta }|{|}^2+\lambda ||\mathbf{\beta }|{|}_1$$

を最小化する回帰モデルは LASSO回帰 とよばれ、このような正則化形式を $L_1$-正則化 といいます。$\hat{\mathbf{\beta }}$ を解析的に書くことは難しいですが、このモデルにはペナルティを受けた回帰係数の値をゼロにするという特徴があります。したがって、LASSO回帰モデルでは不要な説明変数を削ることで過学習を防ぐことができます。

LASSO回帰モデルの代表例として、時系列データの平滑化に用いられる Fused LASSO があります。時系列データ $\{y_t{\}}_{t=1,2,\dots ,T}$ に対応する信号を $\{{\beta }_t{\}}_{t=1,2,\dots ,T}$ とするとき、ペナルティとして

$$\lambda \sum _{t=1}^{T-1}|{\beta }_t-{\beta }_{t-1}|$$

を課すと、隣接する時刻で ${\hat{\beta }}_t$ は同じ値を取りやすくなり、時系列データのグラフは定数関数をつなぎ合わせたような形になります。

参考

---

書籍

• 『統計学実践ワークブック』（学術図書出版社, 2020）第16,17章
• 永田靖・棟近雅彦『多変量解析入門』（サイエンス社, 2001）第5章

脚注

---

1. 損失関数にはさまざまな種類があり、この形は最も基本的なもので 平均2乗誤差 とよばれるようです。 ↩︎

2. ベクトル $\mathbf{x}=(x_1,\dots ,x_n{)}^{\top }$ に対し、$||\mathbf{x}|{|}_d:=(|x_1{|}^d+\cdots +|x_n{|}^d{)}^{1/d}$ で定義される量を $\mathbf{x}$ の $L_d$-ノルムといいます。$L_2$-ノルムは通常の意味での距離（ユークリッド距離）、$L_1$-ノルムはマンハッタン距離に相当します。 ↩︎