【統計】単回帰分析

本記事では、単回帰分析について記載します。

【目次】

  

計算式等


単回帰分析は、2変量 {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 が1増えた時に 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=説明変数の数(単回帰分析では1)

なお、平方和と自由度は  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(目的変数の範囲, 説明変数の範囲, 切片項有無, 補正項有無)
 [引数]
 第1、2引数:特に迷わないと思います。単純なデータ指定です。  第3引数:切片を推定したい場合TRUEとします。FALSEにすると切片=0です。   第4引数:FALSEにすると回帰係数と切片のみを返し、TRUEにすると寄与率、F値、平方和等の情報も計算してくれます。  ※今回の例では第3、4引数は両方ともTRUEにしました。

LINEST関数は5行×2列(列数は切片+説明変数の数)を範囲指定し、Ctrl + Shift + Enter で実行します。

すると、各セルに次の結果が得られます。

出力結果<再掲> 各セルの詳細
f:id:cochineal19:20201230171243p:plain
回帰係数(b1 切片項(b0
bの標準誤差 bの標準誤差
寄与率(R2 誤差項
F値 残差自由度
回帰平方和(SR 残差平方和(Se

出力結果を見てみると  b_{0},\ b_{1} は手計算と同じ結果が得られています。
また、F=18.9384(回帰自由度=1、残差自由度=13)なので、p≒0.0007843となります(セルP12=F.DIST.RT(F値のセル, 1, 残差自由度のセル) )。これによって、説明変数 x は目的変数 y に影響を与える(回帰直線が傾き=0ではない)と判断できそうです。

なお、回帰係数÷標準誤差で t値を求めることができ、こちらを使っても同じく p値≒0.0007843 を得ることができます(セルP10:11)。

 

モデルの評価

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

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

なお、寄与率は「実測値と予測値の相関係数の二乗」に一致します。
このことからも、実測値と予測値の直線関係を評価できることが分かります。

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

※寄与率は説明変数を増やすだけで大きくなってしまうため、その影響を考慮した指標として自由度調整済寄与率があります。重回帰分析や変数選択によるモデル比較等においては、特にこちらが指標になります。自由度調整済寄与率とRMSEは別記事で扱おうと思います。

 

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} について偏微分した式を連立方程式にして最良値を求めます。上のイメージのように、残差平方和は二乗のため放物線の形をとり、傾き=0が最小値となります。
 元: 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)]×[関数の中身の偏微分]。傾き最小を求めるので右辺は0。

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

\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}} ・・・④

よって、③と④から  b_{0}, \ b_{1} は次となります。

\quad b_{1}=\dfrac{S_{xy}}{S_{xx}}
\quad b_{0}=\overline{y} - b_{1}\overline{x}

 

参考


統計学(出版:東京図書), 日本統計学会編

多変量解析実務講座テキスト, 実務教育研究所

市川伸一他, SASによるデータ解析入門 [第3版](出版:東京大学出版会

本ブログは個人メモです。 本ブログの内容によって生じた損害等の一切の責任を負いかねますのでご了承ください。