『統計学実践ワークブック』第16章の勉強メモです。
重回帰分析
ガウス・マルコフモデル
ガウス・マルコフモデルは以下で表される線形回帰モデル：
y=\beta_0+\beta_1 x_1+\cdots +\beta_d x_d +\epsilon
ここで、$y$は目的変数、$x_1,\cdots,x_d$は説明変数、$\beta_1,\cdots,\beta_d\in\mathbb{R}$は回帰係数、$\beta_0\in\mathbb{R}$は切片項、$\epsilon\sim N(0,\sigma^2)$は誤差項。
ガウス・マルコフモデルにおいて、$d=1$のとき単回帰、$d>1$のとき、重回帰という。
最小二乗推定量
ガウス・マルコフモデルの重回帰分析において最も基本的な推定量は最小二乗推定量となる。
推定量$\beta=(\beta_0,\beta_0,\cdots,\beta_d)^\top\in\mathbb{R}^{d+1}$を推定するためにサイズ$n$のデータあると
\begin{array}{cc}
y_i=\beta_0+\beta_1 x_{i,1}+\cdots +\beta_d x_{i,d} +\epsilon_i &(i=1,2,\cdots,n)
\end{array}
とできる。ここで、
\begin{multline}
\begin{split}
{Y}&=(y_1,\cdots,y_n)^\top\in\mathbb{B}^n \\
{X}&=
\begin{pmatrix}
1 & x_{1,1} & x_{1,2} & \cdots & x_{1,d}\\
1 & x_{2,1} & x_{2,2} & \cdots & x_{2,d}\\
\vdots & \vdots & \vdots & \vdots &\vdots \\
1 & x_{n,1} & x_{n,2} & \cdots & x_{n,d}\\
\end{pmatrix}
\in\mathbb{R}^{n\times(d+1)}
\end{split}
\end{multline}
とおく。
推定量は
\hat{\beta}\in\arg_{\beta\in\mathbb{R}^{d+1}}\min \|{Y}-{X}\beta\|
と定義する。
簡単のため$X^\top X$が可逆であるとすると、
\begin{multline}
\begin{split}
\|Y-X\beta\|^2&=(Y-X\beta)^\top(Y-X\beta)\\
&= \|Y\|^2-Y^\top X\beta-(X\beta)^\top Y+(X\beta)^\top X\beta\\
&=\|Y\|^2 -Y^\top X(X^\top X)^{-1}X^\top Y+Y^\top X(X^\top X)^{-1}X^\top Y-Y^\top X\beta-(X\beta)^\top Y+(X\beta)^\top X\beta\\
&=\|Y\|^2 -Y^\top X(X^\top X)^{-1}X^\top Y+(X^\top Y)^\top(X^\top X)^{-1}X^\top Y-(X^\top Y)^\top\beta -\beta^\top(X^\top Y)+\beta^\top X^\top X\beta\\
&=\|Y\|^2 -Y^\top X(X^\top X)^{-1}X^\top Y+(X^\top Y)^\top \{(X^\top X)^{-1}X^\top Y-\beta\}-\beta^\top \{X^\top Y - X^\top X\beta\}\\
&=\|Y\|^2 -Y^\top X(X^\top X)^{-1}X^\top Y+(X^\top Y)^\top \{(X^\top X)^{-1}X^\top Y-\beta\}-\beta^\top (X^\top X) \{(X^\top X)^{-1}X^\top Y - \beta\}\\
&=\|Y\|^2 -Y^\top X(X^\top X)^{-1}X^\top Y +\{(X^\top Y)^\top-\beta^\top (X^\top X)\}\{(X^\top X)^{-1}X^\top Y-\beta\}\\
&=\|Y\|^2 -Y^\top X(X^\top X)^{-1}X^\top Y +\{(X^\top Y)^\top(X^\top X)^{-1}-\beta^\top\} (X^\top X)\{(X^\top X)^{-1}X^\top Y-\beta\}\\
&=\|Y\|^2 -Y^\top X(X^\top X)^{-1}X^\top Y +\{(X^\top X)^{-1}X^\top Y-\beta\}^\top (X^\top X)\{(X^\top X)^{-1}X^\top Y-\beta\}\\
\end{split}
\end{multline}
となる。また、$X^\top X$が正定値対象行列であるので、最小二乗推定量は
\hat{\beta}=(X^\top X)^{-1}X^\top Y
とできる。
回帰分析では、（全変動）＝（回帰変動）＋（残差変動）となることが知られている。
つまり、$\sum_{i=1}^n(y_i-\bar{y_i})^2=\sum_{i=1}^n((x_i-\bar{x})^\top\hat{\beta_{1:d}})^2+\sum_{i=1}^ne_i^2$となる。
ここで、$\bar{y}=\frac{1}{n}\sum_{i=1}^ny_i$、$\bar{x}=\frac{1}{n}\sum_{i=1}^nx_i$、$\hat{\beta_{1:d}}=(\hat{\beta_j})_{j=1}^d$となる。
回帰変動は回帰変数$\hat{\beta}$を用いた予測値の分散に相当し、残差変動は節目見変数では説明できない目的変数の分散に相当する。目的変数の変動を説明変数で説明できる部分とそれ以外の部分に分けることでモデルの説明力を調べることがで、決定係数$R^2$を以下のように定義する。
R^2=\frac{\sum_{i=1}^n((x_i-\bar{x})^\top\hat{\beta_{1:d}})^2}{\sum_{i=1}^n(y_i-\bar{y_i})^2}=1-\frac{\sum_{i=1}^n(y_i-(1,x_i^\top)\hat{\beta})^2}{\sum_{i=1}^n(y_i-\bar{y_i})^2}
決定変数は$0\leq R^2\leq 1$となり、1に近ければデータへの当てはまりが良いことを示す。
決定変数は変数を増やせば増やすほど増大するので、変数の数について調整した自由度調整済み決定変数$R^{*2}$は
R^{*2}=1-\frac{(\sum_{i=1}^n(y_i-(1,x_i^\top)\hat{\beta})^2)/(n-d-1)}{(\sum_{i=1}^n(y_i-\bar{y})^2)/(n-1)}
で定義される。
最小二乗推定量の性質
ガウス・マルコフモデルが真の分布を表現しており、$y$がある${\beta}^* \in \mathbb{R}^{d+1}$を用いて$y={\beta}_0^* +\beta_1^* x_1+\cdots +\beta_d^* x_d+\epsilon$に従っているとする。
- 不偏性
最小二乗推定量は不偏推定量である。
つまり、推定量の期待値が
\begin{multline}
\begin{split}
E[\hat{\beta}]&=E[(X^\top X)^{-1}X^\top Y]\\
&=(X^\top X)^{-1}X^\top E[Y]\\
&=(X^\top X)^{-1}X^\top E[X\beta^* +\epsilon]\\
&=(X^\top X)^{-1}X^\top X\beta^*+(X^\top X)^{-1}X^\top E[\epsilon]\\
&=\beta^*
\end{split}
\end{multline}
真の値となる。
- 最小分散
最小二乗推定量の分散が不偏推定量の中で最小である。
$\hat{\beta}$の分散共分散行列は
\begin{multline}
\begin{split}
Cov[\hat{\beta}]&=Cov[(X^\top X)^{-1}X^\top Y]\\
&=(X^\top X)^{-1}X^\top((X^\top X)^{-1}X^\top)^\top Cov[Y]\\
&=(X^\top X)^{-1}X^\top((X^\top X)^{-1}X^\top)^\top Cov[X\beta^* +\epsilon]\\
&=(X^\top X)^{-1}X^\top((X^\top X)^{-1}X^\top)^\top X\beta^* Cov[1]+(X^\top X)^{-1}X^\top((X^\top X)^{-1}X^\top)^\top Cov[\epsilon]\\
&=(X^\top X)^{-1}X^\top(X(X^\top X)^{-1}) \sigma^2\\
&= \sigma^2(X^\top X)^{-1}
\end{split}
\end{multline}
となる。
任意の不偏推定量$\tilde{\beta}$は$C\in\mathbb{R}^{(d+1)\times n}$を用いて、$\tilde{\beta}=CY$とできる。
不偏推定量であるので、
\begin{multline}
\begin{split}
E[\tilde{\beta}]&=E[CY]=C E[X\beta^* +\epsilon]\\
&= CX\beta^* = \beta^*
\end{split}
\end{multline}
となるので、$CX=I$が成り立つ。ここで、分散共分散行列は、
\begin{multline}
\begin{split}
Cov[\tilde{\beta}]&=E[(\tilde{\beta}-\beta^* )(\tilde{\beta}-\beta^* )^\top]\\
&=E[(CY-\beta^* )(CY-\beta^* )^\top]=E[(C(X\beta^* +\epsilon)-\beta^* )(C(X\beta^* +\epsilon)-\beta^* )^\top]\\
&=E[C\epsilon(C\epsilon)^\top]\\
&=C E[\epsilon\epsilon^\top] C^\top\\
&=\sigma^2 CC^\top
\end{split}
\end{multline}
となる。ここで、$\hat{C}=(X^\top X)^{-1}X^\top$
\begin{multline}
\begin{split}
(C-\hat{C})\hat{C}^\top&=(C-\hat{C})((X^\top X)^{-1}X^\top)^\top=(C-\hat{C})X(X^\top X)^{-1}\\
&=(CX-\hat{C}X)(X^\top X)^{-1}\\
&=(I-(X^\top X)^{-1}X^\top X)(X^\top X)^{-1} = 0
\end{split}
\end{multline}
となるので、
\begin{multline}
\begin{split}
CC^\top&=(C-\hat{C}+\hat{C})(C-\hat{C}+\hat{C})^\top\\
&=(C-\hat{C})(C-\hat{C})^\top +(C-\hat{C})\hat{C}^\top + \hat{C}(C-\hat{C})^\top + \hat{C}\hat{C}^\top\\
&=(C-\hat{C})(C-\hat{C})^\top+ \hat{C}\hat{C}^\top
\end{split}
\end{multline}
となり、任意の$C\in\mathbb{R}^{(d+1)\times n}$について$CC^\top\succeq\hat{C}\hat{C}^\top$が成り立つ。
$\sigma^2\hat{C}\hat{C}^\top=\sigma^2(X^\top X)^{-1}=Cov[\hat{\beta}]$となるので、この分散共分散行列は任意の不偏推定量の中で最小となる。
- 分散の不偏推定量
$\frac{|e|^2}{n-d-1}$は$\sigma^2$の不偏推定量となる。
\begin{multline}
\begin{split}
E[\|e\|^2]&=E[\|Y-\hat{Y}\|^2]=E[\|(I-P_X)Y\|^2]\\
&=E[\|(I-P_X)(X\beta^* +\epsilon)\|^2]=E[\|(I-P_X)\epsilon\|^2]\\
&=\sigma^2 Tr[(I-P_X)^2]=sigma^2 Tr[I-P_X]\\
&=(n-d-1)\sigma^2
\end{split}
\end{multline}
正規方程式
最小二乗推定量は正規方程式からも求めることができる。目的関数$||Y-X\beta||^2$を$\beta$に関して微分し、それが$0$になるとき、最小二乗推定量$\hat\beta$となる。
\|Y-X\beta\|^2= \|Y\|^2-Y^\top X\beta-(X\beta)^\top Y+(X\beta)^\top X\beta\\
ここで、$Y^\top X\beta$はスカラーとなるので、$Y^\top X\beta=(Y^\top X\beta)^\top$となる。
\begin{multline}
\begin{split}
\nabla_{\beta}\|Y-X\beta\|^2&=-2X^\top Y+2X^\top X\beta=0\\
X^\top X\beta&=X^\top Y
\end{split}
\end{multline}
この方程式$X^\top X\beta=X^\top Y$を正規方程式という。
これは$X^\top X$が可逆でなくても成り立つ。この時、一般化逆行列$(X^\top X)^+$を考える。
正規方程式の解は
\hat\beta = (X^\top X)^+X^\top Y
で与えられる。
正規行列による最小二乗推定量の導出は逆行列$(X^\top X)^{-1}$を求める必要があり、計算コストが大きい。そのため、小規模・中規模のデータセットでは正規方程式を使い、大規模なデータセットに対しては最急降下法などによるパラメータ推定を求める必要がある。
非線形回帰
ある非線形基底関数$\phi_k$$(k=1,\cdots,d)$について、
y=\beta_0 + \beta_1\phi_1(x)+\cdots+\beta_d\phi_d(x)+\epsilon
というモデルに対して重回帰分析を適用できる。
非線形な基底としては、$\phi_k(x)=x^k$として多項式回帰や$\phi_k(x)=\cos(kx)$とした三角関数を用いた回帰などがある。
重回帰分析の検定
変数$\beta_k$がゼロか非ゼロかを検定することにより、変数$x_k$が$y$に有意に影響しているかどうかを判断できる。
変数の有意性検定を一般化すると、
$q<d+1$をみたすある$q$に対して、行列$A\in\mathbb{R}^{q\times(d+1)}$はランクが$rank(A)=q$かつ像が$Im(A^\top)\subset Im(X^\top)$を満たすとする。この$A$を用いて
\begin{matrix}
帰無仮説H_0：A\beta = 0 & vs. & 対立仮説H_1：A\beta \neq 0
\end{matrix}
を考える。このとき、$\Theta_0={\beta\in\mathbb{R}^{d+1}|A\beta =0}$として、
R_0^2=\min_{\beta\in\Theta_0}\|Y-X\beta\|^2
とおき、
R_1^2=\min_{\beta\in\mathbb{R}^{d+1}}\|Y-X\beta\|^2
とおく。もし$R_0^2$と$R_1^2$が大きく変わる場合、帰無仮説は間違っていると棄却する。
検定統計量として
T=\frac{(R_0^2-R_1^2)/q}{R_1^2/(n-d-1)}
を考える。帰無仮説のもと
T\sim F(q,n-d-1)
が成り立つことが知られている。ここで、$F(a,b)$は自由度$(a,b)$のF分布を表す。この$T$をF統計量という。
これを用いて変数の仮説検定を構成することができる。
なお、
$$S^{-1}=\left(\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})(x_i-\bar{x})^\top\right)^{-1}$$
の$(k,k')$要素を$S^{k,k'}$と表し、$\sigma^2$の不偏推定量を
\hat{\sigma}^2=\|e\|^2/(n-d-1)=R_1^2/(n-d-1)
とする。
- $\beta_k=0\quad (k\in{1,\cdots,d})$の検定
$A=(0,\cdots,0,1,0,\cdots,0)$($k+1$番目の成分のみが1)とすると
T=\frac{\hat{\beta}_k^2}{\hat{\sigma}^2S^{k,k}/n}\sim F(1,n-d-1)
を得る。$F_\alpha(1,n-d-1)$を有意水準$\alpha\in(0,1)$に対するF分布の上側$100\alpha$%点とすると、
T\geq F_\alpha(1,n-d-1)
のときに$H_0$を棄却すればよい。
- $\beta_0=0$の検定
$A=(1,0,\cdots,0)$とすると、
T=\frac{\hat{\beta}_0^2}{\hat{\sigma}^2(1+\sum_{k,k'=1}^d\bar{x_k}\bar{x_{k'}}S^{k,k'})/n}\sim F(1,n-d-1)
が成り立つ。
- $\beta_1=\beta_2=\cdots =\beta_d=0$の検定
$A=(0|I_d)$とすると
T=\frac{\|X\hat{\beta}-1\bar{y}\|^2/d}{\hat{\sigma}^2}\sim F(d,n-d-1)
が成り立つ。
正則化
$X^\top X$の条件数が悪い場合でも安定した推定を行うために正則化が有用である。代表的な正則化法として、$L_2$-正則化と$L_1$-正則化があり、$L_2$-正則化を用いた重回帰をリッジ回帰、$L_1$-正則化を用いた重回帰をLasso回帰という。
リッジ回帰推定量は
\hat{\beta}_R=\arg \min_{\beta\in\mathbb{R}^{d+1}}(\| Y-X\beta \|^2+\lambda\| \beta \|_2^2)
で与えられる。ここで、$|\beta|2 :=\sqrt{\sum_{j=0}^d\beta_j}$、$\lambda>0$は正則化パラメータ。$\lambda$は交差検証法やMallows'Cp規準によって選択される。
Lass推定量は
\hat{\beta}_L = \arg\min_{\beta\in\mathbb{R}^{d+1}}(\|Y-X\beta\|^2+\lambda\|\beta\|_1)
で与えられる。ここで、$|\beta|1:=\sum_{j=0}^d|\beta_j|$。
正則化項$||\beta||$を加えることで推定量の分散を抑えることができる。そのため、推定に用いるデータに合わせすぎて予測誤差が悪くなる過適合を抑え予測精度を向上させることができる。一方で、正則化が強すぎると推定量が原点に向けて縮小されすぎてしまい予測精度はかえって悪くなる過小適合が起こる。
リッジ正則化において、正則化パラメータを選ぶ規準としてMallows'Cp規準がある。
統計的自由度
M(\lambda)=Tr [X(X^\top X + \lambda I_d)^{-1}X^\top]
を用いて、Mallows'Cp規準は
Cp(\lambda)=\|Y-X\hat{\beta}_R\|^2+2\sigma^2M(\lambda)
で与えられる。Mallows'Cp規準を用いた$\lambda$の選択は$Cp(\lambda)$を最小化することで得られる。
Lasso回帰は$L_1$-正則化を用いることで、$\beta$はスパースになり、凸最適化によって変数選択が可能になり、高次元問題において有効となる。
$L_1$-正則化と$L_2$-正則化を混ぜたElastic-Netと呼ばれる手法があり、正則化項として、
\lambda\left(\alpha\|\beta\|_1+\frac{1-\alpha}{2}\|\beta\|_2^2\right)
を用いる。ただし、$\lambda>0$、$0\leq\alpha\leq 1$。
Lasso回帰の亜種にFused lassoがある。時系列解析で用いられ、$\beta_k$を時刻$k$での信号とし、$\beta_k$と$\beta_{k+1}$がなるべく等しくなるようにする正則化を用いる。その目的関数は
\| Y-X\beta\|^2+\lambda\sum_{k=1}^{d-1}|\beta_k-\beta_{k+1}|
で与えられる。
一般化逆行列
線形代数学における逆行列の一般化。疑似逆行列ともいう。
$m\times n$行列$A$に対し、$A$の随伴行列（複素共役かつ転置行列）を$A^*$とするとき、以下の条件を満たす$n\times m$行列$A^+$はただ一つに定まる。
-
$A$と$A^+$は広義可逆元である
\begin{multline} \begin{split} AA^+A&=A \\ A^+AA^+&=A^+ \end{split} \end{multline} -
$AA^+$および$A^+A$はエルミーと行列である
\begin{multline} \begin{split} (AA^+)^*&=AA^+ \\ (A^+A)^*&=A^+A \end{split} \end{multline}
この行列$A^+$を疑似逆行列という。$A$が正則ならば、逆行列$A^{-1}$はこれらの条件を満たす。
射影
線形空間上の線形写像$P:V\to V$が冪等律：$P^2=P$を満たすとき、写像$P$を射影作用素という。また、$P$が自己共役であるとき、写像$P$は直交射影作用素となる。すなわち、$V$が実ベクトル空間の時は、写像$P$が対象行列：$P^\top =P$の時となる。
最小二乗推定量について
$P_X =X(X^\top X)^{-1}X^\top$とすると、これは直交射影行列となる。
冪等律については
\begin{multline}
\begin{split}
P_X^2 &= X(X^\top X)^{-1}X^\top X(X^\top X)^{-1}X^\top\\
&=X(X^\top X)^{-1}X^\top\\
&=P_X
\end{split}
\end{multline}
となり、成り立つ。対称性については、
\begin{multline}
\begin{split}
P_X^\top &= (X(X^\top X)^{-1}X^\top)^\top \\
&=(X^\top)^\top((X^\top X)^{-1})^\top X^\top\\
&=X((X^\top X)^\top)^{-1} X^\top\\
&=P_X
\end{split}
\end{multline}
となり、成り立つ。
ここで、$\hat Y = X\hat\beta$とおくと、$\hat{Y}=P_X Y$となる。最小二乗推定量$\hat\beta$による$Y$の予測値$\hat Y$は$Y$の$P_X$による射影になっており、実績値と予測値の残差を
e=Y-X\hat\beta =Y-\hat Y=(I-P_X)Y
とする。$e$は$Im(X)$の直交補空間に含まれ、$X^\top e=0$が成り立つ。
\begin{multline}
\begin{split}
X^\top e &= X^\top (I-P_X)Y \\
&=(X^\top - X^\top X(X^\top X)^{-1}X^\top)Y \\
&=(X^\top -X^\top)Y \\
&=0
\end{split}
\end{multline}
また、$X_i=(1,x_{i,1},\cdots,x_{i,d})$より、$\mathbb{1}=(1,1,\cdots,1)\subset X$となるので、$\mathbb{1}^\top e=0$も成り立つ。
これにより、
\mathbb{1}^\top (Y-X\hat\beta)=n(\bar{y}-(1,\bar{x}^\top)\hat\beta)=0
となる。ただし、$\bar{y}:=\frac{1}{n}\sum_{i=1}^n y_i$、$\bar{x}=\frac{1}{n}\sum_{i=1}^n x_i$。よって、$\bar{y}=(1,\bar{x}^\top)\hat\beta$が成り立つ。これより、
Y-\mathbb{1}\bar{y}=X\hat\beta -\mathbb{1}(1,\bar{x}^\top)\hat\beta +e=
\left(
\begin{matrix}
0 & (x_1-\bar{x})^\top \\
\vdots & \vdots \\
0 & (x_n-\bar{x})^\top
\end{matrix}
\right)\hat\beta +e
となる。ここで、$X\hat\beta -\mathbb{1}(1,\bar{x}^\top)\hat\beta \in Im(X)$となるので、$||Y-\mathbb{1}\bar{y}||$は$Im(X)$とその直交補空間に分解できるため、
\|Y-\mathbb{1}\bar{y}\|^2=\|X\hat\beta -\mathbb{1}(1,\bar{x}^\top)\hat\beta\|^2+\|e\|^2
とできる。よって、総変動$\sum_{i=1}^n(y_i-\bar{y})^2$、回帰変動$\sum_{i=1}^n((x_i-\bar{x})^\top\hat\beta_{1:d})^2$、残差変動$\sum_{i=1}^n e_i^2$について、以下のようにできる。
\sum_{i=1}^n(y_i-\bar{y})^2 = \sum_{i=1}^n((x_i-\bar{x})^\top\hat\beta_{1:d})^2 + \sum_{i=1}^n e_i^2
例題
問１６．１
[1]
結果の自由度調整済み決定変数（Adjusted R-squared）が一番大きいモデルを選べばよいので、モデル３となる。
[2]
有意水準$\alpha=0.05$よりP-値が小さいものを選べばよいので、GEN、AMTとなる。
問１６．２
[1]
最小二乗推定量はEstimateなので、0.27388となる。
t検定統計量は$t=\frac{\sqrt{n}\bar{X}}{s}$で求められる。ここでは、最小二乗推定量/標準偏差とすればよいので、$0.27388/0.22967\approx 1.1925$となる。自由度（111-5-1=105）なので、P-値は約0.12となり、有意水準$alpha=0.1$とすると有意とは言えない。
[2]
予測誤差の観点からAICを比べるのが一番良い。よって、モデル１が良い。
問１６．３
[1]
予測誤差の観点から交差検証スコアを見る必要がある。よって、$\lambda=1$となる。
[2]
(a):$\lambda=e^{-2}$、$\alpha=0$：スパース性が弱いため。
(b):$\lambda=0$：分散が一番大きいため。
(c):$\lambda=e^{-2}$、$\alpha=1$：スパース性が一番強いため。
(d):$\lambda=e^{-2}$、$\alpha=0.5$
[3]
(a):$\alpha=0$：スパース性が弱いため
(b):$\alpha=1$：スパース性が強いため
(c):$\alpha=0.5$