重回帰分析と変数増減法

# 概要
モデル

```math
y_i=\beta_0 + \beta_1 x_{i1}+ \beta_2 x_{i2} +\epsilon
```
とモデル0

```math
y_i=\beta_0 + \beta_1 x_{i1}+ \beta_2 x_{i2}+ \beta_3 x_{i3}+\epsilon
```

２つのモデルが建てられたとき、どちらを採用するだろうか。

当てはまり具合は

```math
\begin{eqnarray}
S_e = \sum (y_i - \hat y)^2 : 元モデル\\
S_{e0} = \sum (y_i - \hat y)^2  : モデル0\\
\end{eqnarray}
```

で求められる。一方でモデル0のほうが変数が多く、$S_e \ge S_{e0}$と予想できる。

## 変数増減法
仮定として元のモデル

```math
y_i=\beta_0 + \beta_1 x_{i1}+ \beta_2 x_{i2}
```
が正しいとすると、

```math
\frac{(S_e - S_{e0})/自由度}{(S_e)/自由度} \sim F(a,b)
```
である。式自体の意味は残差平方和$S_e$が変数の追加でどれだけ減少したかということである。

```math
\frac{(S_e - S_{e0})/自由度}{(S_e)/自由度} \sim F(1,2) \ge 2
```
となる確率は

https://keisan.casio.jp/exec/system/1161228871

を元にすると、$29$%ほどである。


----
まとめると

変数$x_3$の追加が無意味であるならば、

```math
P(\frac{(S_e - S_{e0})/自由度}{(S_e)/自由度} \sim F(1,2) \le 2) = 0.71
```
であり、$\frac{(S_e - S_{e0})/自由度}{(S_e)/自由度} \le 2$ならば$x_3$を追加する必要はほぼない。
その程度の残差平方和$S_e$の減少はモデルにとって間違っていても71%で起こり得ることである。
$x_3$を追加する意味はほぼない。

という結論が得られる。

## 71%
今回は71%であったが、元のモデルの変数が多くなるほどこの確率は大きくなる。
実際に計算すると


|2変数|3変数|4変数|5変数|6変数|10変数|
|---|---|---|---|---|---|
|71%|74%|77%|78%|79%|81%|

と求まる。

パチンコをしてる身だと81%を一回目から外すことは体感余裕なので、$\frac{(S_e - S_{e0})/自由度}{(S_e)/自由度} \ge 2$は偶然とりえる感覚です。

そんな事を言ってもどうしようもないですが、2以上なら変数を採用、2以下なら不採用でまぁ悪くないかなと言う感じです。
ちなみに95%点を計算すると

|2変数|3変数|4変数|5変数|6変数|10変数|
|---|---|---|---|---|---|
|18.5|10.1|7.71|6.61|5.99|4.96|

が求まります。

## そもそも２以上とり得るのか？

```math
\frac{(S_e - S_{e0})/自由度}{(S_e)/自由度} 
```
は一見して


```math
1-\frac{( S_{e0})/自由度}{(S_e)/自由度} 
```
なので1以上取らなそうです。


元モデルが変数2の場合は

```math
\frac{(S_e - S_{e0})/自由度}{(S_e)/自由度}=\frac{(S_e - S_{e0})/1}{(S_e)/2} = 2(1-\frac{ S_{e0}}{S_e})
```
となるので、2以上取ることはないです。


# F分布に従うことの証明
ここが一番書きたかった内容です。
F分布に従うことを証明するのにかなり悩みました。
多分また引っかかる日が来るので書いておきます。



モデル

```math
\begin{eqnarray}
y=X\theta +\epsilon\\
y_i=\beta_0 + \beta_1 x_{i1}+ \beta_2 x_{i2}\\

X=
\begin{pmatrix}
1 & x_11 & x12\\
1 & x_21 & x22\\
\vdots \\
1 & x_n1 & xn2\\
\end{pmatrix}
\end{eqnarray}
```
とモデル0

```math
\begin{eqnarray}
y=X_0 \theta_0 +\epsilon\\
y_i=\beta_0 + \beta_1 x_{i1}+ \beta_2 x_{i2}+ \beta_3 x_{i3}\\

X=
\begin{pmatrix}
1 & x_11 & x12 & x13\\
1 & x_21 & x22 & x23\\
\vdots \\
1 & x_n1 & xn2 & xn3\\
\end{pmatrix}
\end{eqnarray}
```

の2つがあります。
現象としては元モデルが正しいとします。


$X$を特異値分解する


```math
X=U\begin{pmatrix}
\Sigma\\
0\\
\end{pmatrix}
V\\
```
残差平方和$S_e$は

```math
S_e = y^T(I - X(X^TX)^{-1}X^T)y\\
=y^T U
\begin{pmatrix}
I & 0\\
0& 0\\
\end{pmatrix}
U^Ty
```
になる。

ここで、$U$行列は$XX^T$の固有ベクトルを並べたユニタリ行列である。

先程

```math
X=
\begin{pmatrix}
1 & x_11 & x12\\
1 & x_21 & x22\\
\vdots \\
1 & x_n1 & xn2\\
\end{pmatrix}
```
としたが、これを

```math
\begin{eqnarray}
X=
\begin{pmatrix}
l_1 & x_1 & x_2\\
\end{pmatrix}\\

l_1 = 
\begin{pmatrix}
1\\
1\\
\vdots\\
1\\
\end{pmatrix}\\
x_1 = 
\begin{pmatrix}
x_{11}\\
x_{21}\\
\vdots\\
x_{n1}\\
\end{pmatrix}\\

\end{eqnarray}
```
としてみる。

```math
XX^T = l_1 l_1^T + x_1 x_1^T + x_2 x_2^T
```
である。左辺に登場するベクトルを変形し、スペクトル分解と一致するようにする。
つまり、

```math
l_2 = x_1 -  ( x_1 \cdot l_1) l_1/n\\
```
とすれば、

```math
l_2 \cdot l_1 = 0
```
になる。
こんな感じで

```math
XX^T =\lambda_1 l_1 l_1^T + \lambda_2 l_2 l_2^T +\lambda_3 l_3 l_3^T
```
とできたとする。

$X_0$に対しても同様の考え方ができる。

```math
\begin{eqnarray}
X_0X_0^T = l_1 l_1^T + x_1 x_1^T + x_2 x_2^T + x_3 x_3^T\\
= \lambda_1 l_1 l_1^T + \lambda_2 l_2 l_2^T +\lambda_3 l_3 l_3^T + \lambda_4 l_4 l_4^T
\end{eqnarray}
```
ここで、$l_1 \sim l_3$のベクトルは$X$でも$X_0$でも同じベクトルである（導出過程よりわかる。自然になったというより、$l_1$から順に直交ベクトルを生成する流れからである）


以上より

```math
\begin{eqnarray}
X=U\begin{pmatrix}
\Sigma\\
0\\
\end{pmatrix}
V\\

X_0=U_0\begin{pmatrix}
\Sigma_0\\
0\\
\end{pmatrix}
V_0\\
\end{eqnarray}
```
それぞれの特異値分解に対して、

```math
\begin{eqnarray}
U=\begin{pmatrix}
l_1 & l_2 & l_3& \cdots\\
\end{pmatrix}\\
U_0=\begin{pmatrix}
l_1 & l_2 & l_3&l_4& \cdots\\
\end{pmatrix}\\

\end{eqnarray}
```
と言える。

残差平方和$S_e$は

```math
\begin{eqnarray}
S_e = y^T(I - X(X^TX)^{-1}X^T)y\\
=y^T U
\begin{pmatrix}
I & 0\\
0& 0\\
\end{pmatrix}
U^Ty\\
=\sum_{i=1}^{3} (l_i \cdot y)^2

\end{eqnarray}
```
モデル0の残差平方和は

```math
\begin{eqnarray}
S_{e0} = y^T(I - X_0 (X_0^T X_0)^{-1} X_0^T )y\\
=y^T U_0

\begin{pmatrix}
I & 0\\
0& 0\\
\end{pmatrix}

U_0^Ty\\

=\sum_{i=1}^{4} (l_i \cdot y)^2

\end{eqnarray}
```

となる。


```math
S_e - S_{e0} = (l_4 \cdot y)^2
```
となる。

## 残差平方和の分布
少し話を戻して、残差平方和の分布の話をしよう


```math
\begin{eqnarray}
e&=& y- \hat y\\
E[e] &=& E[y] - X \hat \beta\\
&=&X \beta - X \hat \beta\\
&=& 0\\

V[e] &=& V[Y-\hat  Y]\\
&=& V[Y - X(X^T X)^{-1} X^T Y]\\
&=& V[(I_n - X(X^T X)^{-1} X^T)Y]\\
&=& (I_n - X(X^T X)^{-1} X^T)Y Y^T(I_n - X(X^T X)^{-1} X^T)^T\\
&=& (I_n - X(X^T X)^{-1} X^T)\sigma^2 I_n (I_n - X(X^T X)^{-1} X^T)^T\\
&=& \sigma^2 (I_n - X(X^T X)^{-1} X^T)\\
&=&\sigma^2 (n-rank(X)) \sim \sigma^2 \chi(n-rank(X))
\end{eqnarray}
```
となる。

ここで、

```math
\begin{eqnarray}
V[e] = E[e^Te]\\
=E[S_e]\\
\end{eqnarray}
```
となる。よって、


```math
S_e \sim   \sigma^2 \chi(n-rank(X))
```
となる。


## もとの話に戻る

```math
\begin{eqnarray}

S_e = y^T(I - X(X^TX)^{-1}X^T)y\\
=y^T U
\begin{pmatrix}
I & 0\\
0& 0\\
\end{pmatrix}
U^Ty\\
=\sum_{i=1}^{3} (l_i \cdot y)^2 \sim \sigma^2 \chi^2 (n-rank(X))
\end{eqnarray}
```
となる。

```math
\begin{eqnarray}

Cov(l_4 \cdot y,l_3 \cdot y) =Cov(l_4^T y,l_3^T y) \\
=l_4^T l_3 Cov(y,y)\\
=0
\end{eqnarray}
```
となる。つまり、$S_e$と$S_{e0}$は**独立である変数の2乗である**と考えることができる。

```math
S_e - S_{e0} = (l_4 \cdot y)^2 \sim \sigma^2  \chi^2 (1)
```
となる。
よって、

```math
\frac{(S_e - S_{e0})/自由度}{(S_e)/自由度} \sim F(a,b)
```
が示せた。

# まとめ

1. F分布に従うことを証明できた

かなり雑に書いてしまった。
まぁでも、誰もあまり見ないだろうから許してほしい。
非常に悩んだ問題だったので解決できてよかった。