【統計】単回帰分析 - こちにぃるの日記

単回帰分析について。

【目次】

計算式等
ノート
プログラムコード
回帰係数の詳細
参考

計算式等

単回帰分析は、２変量 {x,y} の関係を直線（一次式）で求める分析方法。
モデル、式は次のとおり。

■ 単回帰モデル

$\quad y=\beta_{0}+\beta_{1}x + \varepsilon\ ,\quad\sim\varepsilon: N\left(0,\sigma^{2}\right)$

■ 単回帰式

$\quad \widehat{y}_{i}=b_{0}+b_{1}x_{i}$

$\quad b_{0}=\overline{y}-b_{1}\overline{x}$

$\quad b_{1}=\dfrac{S_{xy}}{S_{xx}}=\dfrac{\sum \left( x_{i}-\overline{x}\right) \left( y_{i}-\overline{y}\right)}{\sum \left( x_{i}-\overline{x}\right) ^{2}}$

$b_{0}$ は切片（intercept）、 $b_{1}$ は回帰係数（regression coefficient）であり、回帰係数は x が１増えた時に y がどのくらい増えるかを表す。
また、説明変数で表現できないものを残差 $e_{i}$ と言い、この残差は独立に正規分布 $N\left( 0,\sigma ^{2}\right)$ に従う。

単回帰分析は次の分散分析表で表現することができる。
この分散分析表は一元配置分散分析で示したものと同じ。

■ 分散分析表

変動因	平方和	自由度	平均平方	F比
回帰R	$S_{R}=\sum \left( \widehat{y}_{i}-\overline{y}\right) ^{2}=\dfrac{S_{xy}^{2}}{S_{xx}}$	$f_{R}=k$	$V_{R}=\dfrac{S_{R}}{f_{R}}$	$F_{0}=\dfrac{V_{R}}{V_{e}}$
残差e	$S_{e}=\sum \left( y_{i}-\widehat{y}_{i}\right) ^{2}$	$f_{e}=n-k-1$	$V_{e}=\dfrac{S_{e}}{f_{e}}$
全体T	$S_{T}=\sum \left( y_{i}-\overline{y}\right) ^{2}=S_{yy}$	$f_{T}=n-1$

　k=説明変数の数（単回帰分析では１）

なお、平方和と自由度は $S_{T}=S_{R}+S_{e} ,\ f_{T}=f_{R}+f_{e}$ の関係が成り立つ。

モデルの評価指標としては次のようなものがある（一例）。

寄与率もしくは決定係数（R-square）

　 $R^{2}=1- \dfrac{S_{e}}{S_{T}}=\dfrac{S_{R}}{S_{T}}$
　0～1の値を取り、1に近いほど良い。

自由度調整済寄与率もしくは自由度調整済決定係数（adjusted R-square）

　 $adj\ R^{2}=1- \dfrac{S_{e}/f_{e}}{S_{T} / f_{T}}$
　0～1の値を取り、1に近いほど良い。

RMSE（Root Mean Square Error, 二乗平均平方根誤差）

　 $RMSE=\sqrt{\dfrac{1}{f_{e}}\sum \left( y_{i}-\widehat{y_{i}}\right) ^{2}}$
　相対的な指標であり、値が小さいほど良い。予測値と実測値の誤差に関する指標。

ノート

今回も架空データで簡単に計算方法を見てみる。

　x = {1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3}
　y = {1.1, 2.4, 3.3, 2.5, 1.3, 3.1, 2.7, 4.9, 2.4, 3.5, 5.8, 4.4, 3.1, 4.9, 5.3}

Excelでの計算：

f:id:cochineal19:20201230020235p:plain — 単回帰分析

モデルの作成

$\overline{x}=2,\ \overline{y}=3.38,\ S_{xx}=10,\ S_{xy}=12.9$ となるため、 $b_{0},\ b_{1}$ は次のようになる。

$\quad b_{1}=\dfrac{S_{xy}}{S_{xx}}=\dfrac{12.9}{10}=1.29$
$\quad b_{0}=\overline{y}-b_{1}\overline{x}=3.38 - 1.29 \times 2=0.8$

よって、今回の単回帰式は次のとおり。

$\quad \widehat{y}_{i}=0.8+1.29 x_{i}$

グラフにしたもの。

f:id:cochineal19:20201231134321p:plain

この計算、ExcelのLINEST関数でも実行できるのでメモ（個人の興味）。
貼付したExcelではP4:Q8が該当箇所。

 =LINEST(目的変数の範囲, 説明変数の範囲, 切片項有無, 補正項有無)
 [引数]
　第１、２引数：特に迷わないと思います。単純なデータ指定です。
　第３引数：切片を推定したい場合TRUEとします。FALSEにすると切片＝０です。 
　第４引数：FALSEにすると回帰係数と切片のみを返し、TRUEにすると寄与率、F値、平方和等の情報も計算してくれます。
　※今回の例では第３、４引数は両方ともTRUEにしました。

LINEST関数は５行×２列（列数は切片＋説明変数の数）を範囲指定し、Ctrl + Shift + Enter で実行する。実行後、各セルに次の結果が得られる。

出力結果＜再掲＞

各セルの詳細

回帰係数（b₁）	切片項（b₀）
b₁の標準誤差	b₀の標準誤差
寄与率（R²）	誤差項
F値	残差自由度
回帰平方和（S_R）	残差平方和（S_e）

出力結果を見てみると $b_{0},\ b_{1}$ は手計算と同じ結果が得られている。
また、F＝18.9384（回帰自由度＝1、残差自由度=13）なので、p≒0.0007843となる（セルP12=F.DIST.RT(F値のセル, 1, 残差自由度のセル) ）。

なお、回帰係数÷標準誤差で t値を求め、p値≒0.0007843 を得ることができる（セルP10:11）。

モデルの評価

LINEST関数の結果で得られる寄与率は0.5930。

「寄与率＝回帰平方和÷全体平方和」の式から分かるように、寄与率は本モデルでデータ全体のどのくらいを説明できるかを示す指標（１に近いほど良い）。
今回のモデルでは６割近くを説明できているという解釈になる。

なお、寄与率は「実測値と予測値の相関係数の二乗」に一致する。

$\quad R^{2}= \left( r _{y \widehat{y}} \right) ^{2}$

※寄与率は説明変数を増やすだけで大きくなってしまうため、その影響を考慮した指標として自由度調整済寄与率がある。

Rでの実行

> ads <- data.frame(x=c(1,1,1,1,1,2,2,2,2,2,3,3,3,3,3)
+   　　　　　　　　,y=c(1.1, 2.4, 3.3, 2.5, 1.3, 3.1, 2.7, 4.9, 2.4, 3.5, 5.8, 4.4, 3.1, 4.9, 5.3))
> res <- lm(y ~ x, data=ads)
> summary(res)
Call:
lm(formula = y ~ x, data = ads)

Residuals:
   Min     1Q Median     3Q    Max 
-1.570 -0.735  0.120  0.520  1.520 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.8000     0.6404   1.249 0.233580    
x             1.2900     0.2964   4.352 0.000784 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9374 on 13 degrees of freedom
Multiple R-squared:  0.593,	Adjusted R-squared:  0.5617 
F-statistic: 18.94 on 1 and 13 DF,  p-value: 0.0007843

> #-- グラフはこちら
> ggplot(ads, aes(x=x, y=y)) + geom_point() + geom_smooth(method="lm", se=F)

SASでの実行

data ads;
input x y @@;
cards;
1 1.1 1 2.4 1 3.3 1 2.5 1 1.3 2 3.1 2 2.7 2 4.9 2 2.4 2 3.5 3 5.8 3 4.4 3 3.1 3 4.9 3 5.3
;
run;
proc reg data=ads;
 model y = x;
 output out=OutputDS p=pred r=resid;
run;

f:id:cochineal19:20201230015746p:plain

プログラムコード

■ Rのコード

res <- lm(y ~ x, data=InputDS)
summary(res)

■ SASのコード

proc reg data=[InputDS];
 model y = x;
 output out=[OutputDS] p=pred r=residuals;
run;

■ Pythonのコード

整備中

回帰係数の詳細

$b_{0}, \ b_{1}$ の求め方についてもう少し突っ込んでみる。
精度の良いモデルとは、実測値 $y_{i}$ と予測値 $\widehat{y}_{i}$ との差（残差 $e_{i}$ ）が小さいものであるため、残差平方和が最小の時の $b_{0}, b_{1}$ を用いれば良いと考えられる（最小二乗法）。

$\quad S_{e}=\sum e_{i}^{2}=\sum \left( y_{i}-\widehat{y}_{i}\right) ^{2}=\sum \left( y_{i}-\left( b_{0}+b_{1}x_{i}\right) \right) ^{2}=\sum \left( y_{i}- b_{0}-b_{1}x_{i} \right) ^{2}$

f:id:cochineal19:20201230233701p:plain — 最小値のイメージ

よって、 $b_{0}, \ b_{1}$ について偏微分した式を連立方程式にして最良値を求める。
上のイメージのように、残差平方和は二乗のため放物線の形をとり、傾き＝０が最小値となる。
　元： $S_{e}=\sum \left( y_{i}- b_{0}-b_{1}x_{i} \right) ^{2}$

$\quad \begin{cases} \dfrac{\partial S_{e}}{\partial b_{0}}=\sum 2\left( y_{i}-b_{0}-b_{1}x_{i}\right) \left(-1\right)=-2\sum \left( y_{i}-b_{0}-b_{1}x_{i}\right)=0 \\ \dfrac{\partial S_{e}}{\partial b_{1}}=\sum 2 \left( y_{i}-b_{0}-b_{1}x_{i}\right)\left(-x_{i}\right)=-2\sum x_{i} \left( y_{i}-b_{0}-b_{1}x_{i}\right)=0 \end{cases}$
　※合成関数の微分により、[n]×[関数^(n-1)]×[関数の中身の偏微分]。傾き最小を求めるので右辺は０。

ここから -2 を除去すると、具体的な名称として「残差の合計」と「残差×説明変数の合計」が０になる解を求めることが分かる。

$\quad \begin{cases}\sum \left( y_{i}-b_{0}-b_{1}x_{i}\right) =\sum e_{i}=0\\ \sum x_{i}\left( yi-b_{0} -b_{1} x_{i}\right) =\sum x_{i}e_{i}=0\end{cases}$

この連立方程式を変形すると次のようになり、
$\quad \begin{cases}\left( \sum 1\right) b_{0}+\left( \sum x_{i}\right) b_{1}=\sum y_{i}　・・・①\\ \left( \sum x_{i}\right) b_{0}+\left( \sum x_{i}^{2}\right) b_{1}=\sum x_{i}y_{i}　・・・②\end{cases}$

①の式を用いて、 $b_{0}, \ b_{1}$ を次の関係式で表すことができる。
$\quad \sum 1=n,\ \dfrac{\sum x_{i}}{n}=\overline{x}, \ \dfrac{\sum y_{i}}{n}=\overline{y}\ であるため、$

$\quad \left( \sum 1\right) b_{0}+\left( \sum x_{i}\right) b_{1}=\sum y_{i}$
$\qquad \Leftrightarrow nb_{0}+\left( \sum x_{i}\right) b_{1}=\sum y_{i}$
$\qquad \Leftrightarrow b_{0}+b_{1}\overline{x}=\overline{y}$
$\qquad \Leftrightarrow b_{0}=\overline{y} - b_{1}\overline{x}　・・・③$

この③を②に代入すると $b_{1}$ が次のように求まる。

$\quad \left( \sum x_{i}^{2}-\overline{x}\sum x_{i}\right) = S_{xx}, \ \left(\sum x_{i}y_{i}-\overline{y}\sum x_{i} \right) = S_{xy} \ であるため、$

$\quad \left( \sum x_{i}\right) b_{0}+\left( \sum x_{i}^{2}\right) b_{1}=\sum x_{i}y_{i}$
$\qquad \Leftrightarrow \left( \sum x_{i}\right) \left( \overline{y}-b_{1}\overline{x}\right) +\left( \sum x_{i}^{2}\right) b_{1}=\sum x_{i} y_{i}$
$\qquad \Leftrightarrow \left( \sum x_{i}^{2}-\overline{x}\sum x_{i}\right) b_{1}=\sum x_{i}y_{i}-\overline{y}\sum x_{i}$
$\qquad \Leftrightarrow S_{xx}b_{1}=S_{xy}$
$\qquad \Leftrightarrow b_{1} = \dfrac {S_{xy}}{S_{xx}}　・・・④$