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

重回帰分析について。

【目次】

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

単回帰分析については次の記事。

cochineal19.hatenablog.com

計算式等

重回帰分析は、多変量 { $x_{1},\ x_{2},\ ...,\ x_{k},\ y$ } の関係を直線で求める分析方法。
モデルは次のとおり。

■ 重回帰モデル
$\quad y=\beta_{0}+\beta_{1}x_{1}+\beta_{2}x_{2}+...+\beta_{k}x_{k} + \varepsilon\ ,\quad\sim\varepsilon: N\left(0,\sigma^{2}\right)$

■ 重回帰式
$\quad \widehat{y}_{i}=b_{0}+b_{1}x_{1i}+b_{2}x_{2i}+...+b_{k}x_{ki}$

$\quad b_{0}=\overline{y} - \left( b_{1}\overline{x_{1}}+b_{2}\overline{x_{2}}+...+b_{k}\overline{x_{k}} \right)$
$\quad \begin{bmatrix} b_{1} \\ b_{2} \\ \ \vdots \\ b_{k} \end{bmatrix} = \begin{bmatrix} S_{x_{1}x_{1}} \quad S_{x_{1}x_{2}} \quad \ldots \quad S_{x_{1}x_{k}} \\ S_{x_{2}x_{1}} \quad S_{x_{2}x_{2}} \quad \ldots \quad S_{x_{2}x_{k}} \\ \ \vdots \qquad \ \vdots \qquad \ \ddots \qquad \ \vdots \\ S_{x_{k}x_{1}} \quad S_{x_{k}x_{2}} \quad \ldots \quad S_{x_{k}x_{k}} \end{bmatrix}^{-1} \begin{bmatrix} S_{x_{1}y} \\ S_{x_{2}y} \\ \ \vdots \\ S_{x_{k}y} \\ \end{bmatrix}$

$b_{0}$ は切片（intercept）、 $b_{1},\ b_{2},\ ...,\ b_{k}$ は回帰係数（regression coefficient）であり、回帰係数は x が１増えた時に y がどのくらい増えるかを表す。
回帰係数の式は、「回帰係数」＝「説明変数間の分散共分散行列の逆行列」×「説明変数と目的変数の共分散行列」。

なお、正規方程式から次のように回帰係数を求めることもできる。
$\quad \left( X^{t}X\right) \widehat{\beta }=X^{t}Y \verb| （正規方程式）であるから、|$
$\quad \widehat{\beta }=\left( X^{t}X\right) ^{-1}\left( X^{t}Y\right)$
　※X：説明変数の行列、Y：目的変数の行列、β：回帰係数の行列、t：転置

また、説明変数で表現できないものを残差 $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）

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

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

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

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

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

AIC（Akaike's Information Criterion, 赤池情報量規準）
$\quad AIC = -2 \times \left( MLL \right) + 2 \times \left(k \right)$
　　※ MLL : 最大対数尤度（Maximum Log Likelihood）, k : 説明変数の数
　相対的な指標であり、値が小さいほど良い。2×(k) が説明変数の数に対する罰則項である。

VIF（Variance Inflation Factor, 分散拡大係数）
$\quad VIF_{i} = \dfrac{1}{1−r_{i}^{2}}$
　　※ r : 相関係数
　多重共線性の評価に用いる。目安として10以上で多重共線性があると考える。

ノート

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

　x1 = {1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3}
　x2 = {2, 1, 1, 2, 1, 3, 4, 2, 4, 3, 3, 6, 6, 4, 5}
　x3 = {5, 6, 5, 2, 3, 5, 7, 5, 4, 6, 5, 9, 6, 5, 5}
　x4 = {2, 1.5, 2, 2, 1.5, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4}
　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:20210101231218p:plain — 重回帰分析

モデルの作成

分散共分散行列と逆行列は次のとおり。

f:id:cochineal19:20210101231437p:plain

$\quad \begin{bmatrix} b_{1} \\ b_{2} \\ b_{3} \\ b_{4} \end{bmatrix} = \begin{bmatrix} 0.714 \quad 1.214 \quad 0.643 \quad 0.786 \\ 1.214 \quad 2.838 \quad 1.471 \quad 1.367 \\ 0.643 \quad 1.471 \quad 2.600 \quad 0.693 \\ 0.786 \quad 1.367 \quad 0.693 \quad 0.888 \end{bmatrix}^{-1} \begin{bmatrix} 0.921 \\ 0.989 \\ 0.669 \\ 1.031 \\ \end{bmatrix} = \begin{bmatrix} \ \ \ 53.431 \quad \ \ 0.058 \quad -0.786 \quad -46.747 \\ \ \ \ \ 0.058 \quad \ \ 1.526 \quad -0.301 \quad \ -2.165 \\ \ -0.786 \quad -0.301 \quad \ \ 0.556 \quad \ \ \ \ \ 0.725 \ \ \\ -46.747 \quad -2.165 \quad \ \ 0.725 \quad \ \ \ 45.249 \end{bmatrix}\begin{bmatrix} 0.921 \\ 0.989 \\ 0.669 \\ 1.031 \\ \end{bmatrix} = \begin{bmatrix} 0.582 \\ -0.871 \\ 0.097 \\ 1.910 \\ \end{bmatrix}$

$\quad \overline{y}=3.38,\ \overline{x_{1}}=2,\ \overline{x_{2}}=3.13,\ \overline{x_{3}}=5.2,\ \overline{x_{4}}=2.93\quadであるため、$
$\quad b_{0}=3.380 - 0.582 \times 2.0 - \left(-0.871\right) \times 3.133 - 0.097 \times 5.2 - 1.910 \times 2.933 = -1.163$

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

$\quad \widehat{y}_{i}=-1.163 - 0.582 x_{1i} + 0.871 x_{2i} - 0.097 x_{3i} - 1.910 x_{4i}$

モデルの評価

今回の結果についてVIFを調べてみると x1=38.2, x2=4.3, x3=1.5, x4=40.2 であり、x1 と x4 に多重共線性が疑われる（VIF≧10を目安）。x1 と x4 は相関係数=0.987 であり、かなり強い相関関係がある様子。
したがって、x1 と x4 の一方を外したモデルを作って再評価してみる。また、x3 が非有意（p>0.05）だったため、x3 を外したモデルも作ってみる。
（※この記事では変数選択の手法には触れない。）

＜評価モデル＞
　res1 : y = x1 + x2 + x3 + x4
　res2 : y = x1 + x2 + x3
　res3 : y = x2 + x3 + x4
　res4 : y = x1 + x2
　res5 : y = x2 + x4

実行結果は次のとおり。

	res1	res2	res3	res4	res5
VIF≧10	あり	なし	なし	なし	なし
R²	0.853	0.812	0.849	0.808	0.839
adj R²	0.794	0.761	0.808	0.776	0.812
RMSE	0.643	0.692	0.620	0.670	0.613
AIC	35.256	36.873	33.573	35.192	32.558

x1 と x4 の一方を外すことで VIF≧10 の変数は無くなり、多重共線性の問題が解消されたと考えられる。
また、自由度調整済寄与率 $adj\ R^{2}$ の最大値、RMSEとAICの最小値は何れも res5 であり、今回は res5: y=x2 + x4 が最良モデルと考えられる。
なお、res5のモデルは統計的に有意であり（p<0.001）、各回帰係数も有意（x2: p<0.05, x4: p<0.001）。
※今回のモデル比較では、説明変数の数が異なるため、寄与率は $R^{2}$ は使わず、自由度調整済寄与率 $adj\ R^{2}$ を評価指標としている。

Rでの実行

> ads <- data.frame(x1  = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3)
+                   ,x2 = c(2, 1, 1, 2, 1, 3, 4, 2, 4, 3, 3, 6, 6, 4, 5)
+                   ,x3 = c(5, 6, 5, 2, 3, 5, 7, 5, 4, 6, 5, 9, 6, 5, 5)
+                   ,x4 = c(2, 1.5, 2, 2, 1.5, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4)
+                   ,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))
> res1 <- lm(y ~ x1 + x2 + x3 + x4, data=ads)
> summary(res1)
Call:
lm(formula = y ~ x1 + x2 + x3 + x4, data = ads)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.8829 -0.3648 -0.2275  0.4359  0.9455 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -1.16286    1.11986  -1.038  0.32354   
x1           0.58168    1.25690   0.463  0.65342   
x2          -0.87065    0.21243  -4.099  0.00215 **
x3           0.09705    0.12824   0.757  0.46665   
x4           1.91008    1.15667   1.651  0.12968   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6434 on 10 degrees of freedom
Multiple R-squared:  0.8525,	Adjusted R-squared:  0.7935 
F-statistic: 14.45 on 4 and 10 DF,  p-value: 0.0003674

> AIC(res1)
[1] 35.25556
> vif(res1)
       x1        x2        x3        x4 
38.165148  4.331555  1.446171 40.185695 

#他のモデル
# res2 <- lm(y ~ x1 + x2 + x3, data=ads)
# res3 <- lm(y ~ x2 + x3 + x4, data=ads)
# res4 <- lm(y ~ x1 + x2, data=ads)
# res5 <- lm(y ~ x2 + x4, data=ads)

SASでの実行

data ads;
input x1 x2 x3 x4 y @@;
cards;
1 2 5 2 1.1 1 1 6 1.5 2.4 1 1 5 2 3.3 1 2 2 2 2.5 1 1 3 1.5 1.3
2 3 5 3 3.1 2 4 7 3 2.7 2 2 5 3 4.9 2 4 4 3 2.4 2 3 6 3 3.5
3 3 5 4 5.8 3 6 9 4 4.4 3 6 6 4 3.1 3 4 5 4 4.9 3 5 5 4 5.3
;run;
proc reg data=ads;
 model y = x1 x2 x3 x4 / SELECTION=ADJR RSQUARE AIC RMSE VIF;
run;

f:id:cochineal19:20210102001834p:plain

プログラムコード

■ Rのコード

library(car)
res <- lm(y ~ x1 + x2, data=InputDS)
summary(res)
vif(res) #-- VIFを算出
AIC(res) #-- AICを算出

■ SASのコード

ods listing close;
ods output SubsetSelSummary=OutputSummary; /* 全変数パターンの評価指標値を取得 */
proc reg data=[InputDS];
 model y = x1 x2 / selection=ADJR RSQUARE AIC RMSE VIF; /* VIFはmodelの式のみ */
 output out=[OutputDS] p=pred r=residuals; /* modelに書いた式での予測値・残差取得 */
run;quit;
ods listing;

■ Pythonのコード

整備中

偏回帰係数の詳細

$b_{0}, \ b_{1}, \ b_{2}, \ ... \ , b_{k}$ の求め方についてもう少し突っ込んでみる。

精度の良いモデルとは、残差 $e_{i}$ が小さいモデルであるため、残差平方和が最小の時の $b_{0}, \ b_{1}, \ b_{2}, \ ... \ , b_{k}$ を用いれば良いと考えられる（最小二乗法）。

$\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}+b_{2}x_{i}+...+b_{k}x_{k\ i}\right) \right) ^{2}$

したがって、 $b_{0}, \ b_{1}, \ b_{2}, \ ... \ , b_{k}$ について偏微分した式を連立方程式にして最良値を求める。
　元： $\quad S_{e}=\sum \left( y_{i}-\left( b_{0}+b_{1}x_{1i}+b_{2}x_{2i}+...+b_{k}x_{ki}\right) \right) ^{2}$

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

ここから -2 を除去して変形すると次のようになり、

$\quad \begin{cases}\sum \left( y_{i}\right) =\sum \left( b_{0}+b_{1}x_{1i}+b_{2}x_{2i}+...+b_{k}x_{ki}\right) 　・・・①\\ \sum \left( x_{1i}y_{i}\right) =\sum x_{1i}\left( b_{0}+b_{1}x_{1i}+b_{2}x_{2i}+...+b_{k}x_{ki}\right) 　・・・②\\ \sum \left( x_{2i}y_{i}\right) =\sum x_{2i}\left( b_{0}+b_{1}x_{1i}+b_{2}x_{2i}+...+b_{k}x_{ki}\right) 　・・・③\\ \ \ :\\ \sum \left( x_{ki}y_{i}\right) =\sum x_{ki}\left( b_{0}+b_{1}x_{1i}+b_{2}x_{2i}+...+b_{k}x_{ki}\right) 　・・・④\end{cases}$

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

$\quad \sum \left( y_{i}\right) =\sum \left( b_{0}+b_{1}x_{1i}+b_{2}x_{2i}+...+b_{k}x_{ki}\right)$
$\qquad \Leftrightarrow \sum \left( y_{i}\right)=b_{0}\sum \left( 1\right)+b_{1} \sum \left( x_{1i}\right)+b_{2} \sum \left( x_{2i}\right)+...+b_{k} \sum \left( x_{ki}\right)　←nで割る。$
$\qquad \Leftrightarrow \overline{y}=b_{0}+b_{1} \overline{x_{1}}+b_{2} \overline{x_{2}}+...+b_{k} \overline{x_{k}}$
$\qquad \Leftrightarrow b_{0}=\overline{y} - \left( b_{1}\overline{x_{1}}+b_{2}\overline{x_{2}}+...+b_{k}\overline{x_{k}} \right)　・・・⑤$

この⑤を②に代入して、次の関係式を求める。

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

$\quad \sum \left( x_{1i}y_{i}\right) =\sum x_{1i}\left( b_{0}+b_{1}x_{1i}+b_{2}x_{2i}+...+b_{k}x_{ki}\right)$
$\qquad \Leftrightarrow \sum \left( x_{1i}y_{i}\right) =\sum x_{1i}[ \{ \overline{y} - \left( b_{1}\overline{x_{1}}+b_{2}\overline{x_{2}}+...+b_{k}\overline{x_{k}} \right) \} +b_{1}x_{1i}+b_{2}x_{2i}+...+b_{k}x_{ki}]$
$\qquad \Leftrightarrow \left( \sum x_{1i}y_{i} - \overline{y}\sum x_{1i} \right) = b_{1} \left( \sum x_{1i} x_{1i} - \overline{x_{1}} \sum x_{1i} \right) + b_{2} \left( \sum x_{1i} x_{2i} - \overline{x_{2}} \sum x_{1i} \right) + ... + b_{k} \left( \sum x_{1i} x_{ki} - \overline{x_{k}} \sum x_{1i} \right)$
$\qquad \Leftrightarrow S_{x_{1}y}=b_{1}S_{x_{1}x_{1}}+b_{2}S_{x_{1}x_{2}}+...+b_{k}S_{x_{1}x_{k}}　・・・⑥$

③④も同じ処理で次の関係式を求める。
$\quad S_{x_{2}y}=b_{1}S_{x_{2}x_{1}}+b_{2}S_{x_{2}x_{2}}+...+b_{k}S_{x_{2}x_{k}}　・・・⑦$
$\quad S_{x_{k}y}=b_{1}S_{x_{k}x_{1}}+b_{2}S_{x_{k}x_{2}}+...+b_{k}S_{x_{k}x_{k}}　・・・⑧$

これら⑥⑦⑧を行列表現すると次のようになり、
$\quad \begin{bmatrix} S_{x_{1}x_{1}} \quad S_{x_{1}x_{2}} \quad \ldots \quad S_{x_{1}x_{k}} \\ S_{x_{2}x_{1}} \quad S_{x_{2}x_{2}} \quad \ldots \quad S_{x_{2}x_{k}} \\ \ \vdots \qquad \ \vdots \qquad \ \ddots \qquad \ \vdots \\ S_{x_{k}x_{1}} \quad S_{x_{k}x_{2}} \quad \ldots \quad S_{x_{k}x_{k}} \end{bmatrix} \begin{bmatrix} b_{1} \\ b_{2} \\ \ \vdots \\ b_{k} \end{bmatrix} = \begin{bmatrix} S_{x_{1}y} \\ S_{x_{2}y} \\ \ \vdots \\ S_{x_{k}y} \\ \end{bmatrix}$

さらに変形して、各回帰係数の導出式を得る。
$\quad \begin{bmatrix} b_{1} \\ b_{2} \\ \ \vdots \\ b_{k} \end{bmatrix} = \begin{bmatrix} S_{x_{1}x_{1}} \quad S_{x_{1}x_{2}} \quad \ldots \quad S_{x_{1}x_{k}} \\ S_{x_{2}x_{1}} \quad S_{x_{2}x_{2}} \quad \ldots \quad S_{x_{2}x_{k}} \\ \ \vdots \qquad \ \vdots \qquad \ \ddots \qquad \ \vdots \\ S_{x_{k}x_{1}} \quad S_{x_{k}x_{2}} \quad \ldots \quad S_{x_{k}x_{k}} \end{bmatrix}^{-1} \begin{bmatrix} S_{x_{1}y} \\ S_{x_{2}y} \\ \ \vdots \\ S_{x_{k}y} \\ \end{bmatrix} 　・・・⑨$

つまり、「回帰係数」＝「説明変数間の分散共分散行列の逆行列」×「説明変数と目的変数の共分散行列」。

よって、⑤⑨から $b_{0}, \ b_{1}, \ b_{2}, \ ... \ , b_{k}$ は次となる。

なお、分散共分散行列を相関行列にすると標準偏回帰係数が得られる。

・・・その他、
①②③④を行列表現すると次のようになり、この方程式からも回帰係数を導出できる（この方程式を「正規方程式」と言う）。
＜再掲＞
$\quad \begin{cases}\sum \left( y_{i}\right) =\sum \left( b_{0}+b_{1}x_{1i}+b_{2}x_{2i}+...+b_{k}x_{ki}\right) 　・・・①\\ \sum \left( x_{1i}y_{i}\right) =\sum x_{1i}\left( b_{0}+b_{1}x_{1i}+b_{2}x_{2i}+...+b_{k}x_{ki}\right) 　・・・②\\ \sum \left( x_{2i}y_{i}\right) =\sum x_{2i}\left( b_{0}+b_{1}x_{1i}+b_{2}x_{2i}+...+b_{k}x_{ki}\right) 　・・・③\\ \ \ :\\ \sum \left( x_{ki}y_{i}\right) =\sum x_{ki}\left( b_{0}+b_{1}x_{1i}+b_{2}x_{2i}+...+b_{k}x_{ki}\right) 　・・・④\end{cases}$

＜正規方程式・回帰係数の解＞
$\quad \left( X^{t}X\right) \widehat{\beta }=X^{t}Y \quad\left(正規方程式\right)$
$\quad \widehat{\beta }=\left( X^{t}X\right) ^{-1}\left( X^{t}Y\right)$
　※X：説明変数の行列、Y：目的変数の行列、β：回帰係数の行列、t：転置