アルゴリズム (線形多重回帰)

線形多重回帰（重回帰）モデル

線形多重回帰（重回帰）モデルは、単回帰モデルを拡張したもので、複数の独立変数が存在する場合に使用されます。これは、複数の独立変数 $x_1, x_2, \dots, x_k$ が従属変数yに及ぼす影響を解析するために用いられます。与えられたデータセット $(y, x_1, x_2, \dots, x_k)$ に対して、線形多重回帰は以下のモデルにフィットさせます：

$y_i=\beta _0+\beta _1x_{1_i}+\beta _2x_{2_i}+\ldots +\beta _kx_{k_i}+\varepsilon_i$	(1)

ここで、 $\beta _0\,\!$ はy切片で $\beta _1\,\!$ $\beta _2\,\!$ ,...はパラメータです。 $\beta _k\,\!$ は部分回帰係数と呼ばれます。次のように行列形式で書くこともできます。

$Y=XB+E\,\!$	(2)

ここで

$Y=\begin{bmatrix} y_1\\ y_2\\ \vdots \\ y_n \end{bmatrix}$ , $X=\begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1k}\\ 1 & x_{21} & x_{22} & \cdots & x_{2k}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_{n1} & x_{n2} & \cdots & x_{nk} \end{bmatrix}$ $B=\begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_k \end{bmatrix}$ , $E=\begin{bmatrix} \varepsilon _1\\ \varepsilon _2\\ \vdots \\ \varepsilon _n \end{bmatrix}$

ここで $\varepsilon_i$ は、 $\bar{E}=0$ かつ $Var[E]=\sigma^2$ の正規分布に従う独立同分布と仮定します。 $B$ に関して $\left \|E\right \|$ を最小化するために、次の関数を解きます。

$\frac{\partial E'E}{\partial B}=0$	(3)

結果 $\hat B$ はベクトルBの最小二乗推定値であり、次の線形方程式を解いたものです。

$\hat B=\begin{bmatrix} \hat \beta_0 \\ \hat \beta_1 \\ \vdots \\ \hat \beta_k \end{bmatrix}=(X'X)^{-1}X^{\prime }Y$	(4)

ここで、X'はXの転置行列です。与えられたXに対するYの予測値は次のようになります。

$\hat{Y}=X\hat{B}$	(5)

式(4)の $P$ を代入することで行列 $P$ を定義することができます。

$\hat{Y}=[X(X'X)^{-1}X']Y=PY$	(6)

残差は次のように定義されます。

$res_i=Y-\hat{Y}$	(7)

そして残差平方和は次のように表されます。

$RSS=\left \\|E \right \\|^2={Y}'Y-\hat{B}'X'X\hat{B}$	(8)

フィット制御

誤差を重みとする

フィット時に各 $y_i$ に重みを与えることができます。yEr±誤差列 $\sigma_i$ がある場合、この列は各 $y_i$ に対する重み $w_i$ として扱われます。もしyEr±列が存在しない場合、すべての $i$ で $w_i$ は 1 と見なされます。

重み付きフィットにおける解 $\hat{B}$ は、次のように表されます。

$\hat{B}=(X'WX)^{-1}X'WY$	(9)

ここで

$W=\begin{bmatrix} w_1& 0 & \dots &0 \\ 0 & w_2 & \dots &0 \\ \vdots& \vdots &\ \ddots &\vdots \\ 0& 0 &\dots & w_n \end{bmatrix}$

重み付けなし

エラーバーは計算において重みとして扱われません。

直接重み付け

$w_i=\sigma_i$	(10)

機械的

$w_i=\frac 1{\sigma_i^2}$	(11)

切片固定

切片 $\beta_0$ を固定値に設定すると、その分の自由度が減少し、自由度n*=n-1になります。

二乗値（自由度あたりカイ二乗値）でのスケールエラー

重み付きでフィッティングする場合、二乗値（自由度あたりカイ二乗値）でのスケールエラーが利用できます。このチェックボックスはパラメータの誤差推定値にのみ影響し、フィット計算やデータ自体には影響を与えません。デフォルトでオンになっており、 $E$ の分散 $\sigma^2$ を考慮してパラメータ誤差が計算されます。オフにすると、 $E$ は考慮されません。共分散行列を例にとると、

二乗値（自由度あたりカイ二乗値）でのスケールエラー

$\sigma^2=\frac{RSS}{n^{*}-k}$
$Cov(\beta _i,\beta _j)=\sigma^2 (X^{\prime }X)^{-1}$	(12)

二乗値（自由度あたりカイ二乗値）でのスケールエラーでない

$Cov(\beta _i,\beta _j)=(X'X)^{-1}\,\!$	(13)

重み付けフィットでは、 $(X'X)^{-1}\,\!$ の代わりに $(X'WX)^{-1}\,\!$ が使用されます。

フィット結果

フィットパラメータ

フィット値

式 (4)

パラメータの標準誤差

各パラメータにおいて、標準誤差は以下のように得られます。

$s_{\hat \beta _j}=s_\varepsilon \sqrt{C_{jj}}$	(14)

ここで、 $C_{jj}$ は $(X'X)^{-1}$ のj番目の対角要素です（重み付きフィットでは $(X'WX)^{-1}$ を使用します）。残差の標準偏差 $s_\varepsilon$ （“標準誤差”や “推定の標準誤差”、“MSEの平方根”とも呼ばれます）は以下で求めます。

$s_\varepsilon =\sqrt{\frac{RSS}{df_{Error}}}=\sqrt{\frac{RSS}{n^{*}-k}}$	(15)

この $s_\varepsilon^2$ は、 $E$ の分散 $\sigma ^2$ の推定値です。

Note:自由度(df)の詳細については、ANOVA表のdfErrorを参照してください。

t値と信頼水準

回帰の前提が成り立つと仮定すれば、回帰係数に対して以下の帰無仮説と対立仮説に基づく t検定を行うことができます。

$H_0:\beta _j=0\,\!$

$H_\alpha :\beta _j\neq 0$

t値は次式で求められます。

$t=\frac{\hat \beta _j-0}{s_{\hat \beta _j}}$	(16)

計算されたt値を使って、対応する帰無仮説を棄却するかどうかを決定できます。通常与えられた有意水準 $\alpha\,\!$ に対して、 $|t|>t_{\frac \alpha 2}$ のときに $H_0\,\!$ を棄却します。さらに、p値は $アルファ値$ より小さくなります。

Prob>|t|

この値は、t検定における帰無仮説 $H_0$ が正しい確率（p値）を示します。

$prob=2(1-tcdf(\|t\|,df_{Error}))\,\!$	(17)

ここで、 $tcdf(|t|,df_{Error})$ は、値|t|におけるスチューデントt分布の累積分布関数を、誤差の自由度 $(df_{Error})$ で計算します。

LCLおよびUCL

t値から各パラメータの信頼区間 $(1-\alpha )\times 100\%$ を次の式で計算できます。

$\hat \beta _j-t_{(\frac \alpha 2,n^{}-k)}\varepsilon _{\hat \beta _j}\leq \hat \beta _j\leq \hat \beta _j+t_{(\frac \alpha 2,n^{}-k)}\varepsilon _{\hat \beta _j}$	(18)

ここで $UCL$ と $LCL$ はそれぞれ上限信頼区間と下限信頼区間の略です。

CI半幅

信頼区間の半幅は以下の通りです。

$CI=\frac{UCL-LCL}2$	(19)

VIF

分散拡大係数は以下のように計算されます。

$VIF_i=\frac{1}{1-R^2_i}$	(20)

ここで、 $R^2_i$ はi番目の独立変数を他のすべての独立変数で回帰した際の決定係数（補正なし）です。

フィット統計

以下に、主要なフィット統計の数式をまとめます。

自由度

誤差変動の自由度。詳細はANOVA表を参照してください。

自由度あたりカイ二乗

$\sigma^2=\frac{RSS}{n^{*}-k}$	(21)

残差平方和

残差平方和については、式 (8) を参照してください。

R二乗(COD)

フィットの良さ（適合度）は、決定係数(COD) $R^2$ によって評価できます。以下の式で定義されます。

$R^2=\frac{Explained\, variation}{Total \, variation}=1-\frac{RSS}{TSS}$	(22)

補正R二乗

補正 $R^2$ は自由度を考慮して $R^2$ を調整した指標です。以下の式で求められます。

$\bar R^2=1-\frac{RSS/df_{Error}}{TSS/df_{Total}}$	(23)

R値

R値は単純に $R^2$ の平方根です。

$R=\sqrt{R^2}$	(24)

Root-MSE(SD)

Root-MSE（残差標準偏差）は以下の式で求められます

$RootMSE=\sqrt{\frac{RSS}{n^*-k}}$	(25)

残差ノルム

残差平方和の平方根です。

$Norm \,of \,Residuals=\sqrt{RSS}$	(26)

ANOVA表

線形フィットのANOVA表は次のとおりです。

	df	平方和	平均平方	F値	Prob > F
モデル	k	$SS_{reg} = TSS-RSS$	$MS_{reg} = SS_{reg} / k$	$MS_{reg} / MSE$	P値
誤差	n - k*	$RSS$	$MSE = RSS / (n^ - k)$*
合計	n*	$TSS$

Note:モデルに切片（定数項）を含む場合、自由度n*=n-1。含まない場合は、n *= nであり、で、TSS（全体平方和）は「補正なし」のものになります。

ここで、総平方和（TSS）は次のようになります。

$TSS=\sum_{i=1}^n w_iy_i^2$ (補正なし)
$TSS =\sum_{i=1}^nw_i(y_i -\frac{\sum_{i=1}^n w_i y_i} {\sum_{i=1}^n w_i})^2$ (補正あり)	(27)

ここでのF値は、フィッティングモデルが定数モデル（y = 定数）と有意に異なるかどうかを検定するための指標です。

また、p値または有意水準がF検定とともに出力されます。p 値があらかじめ設定された有意水準 $\alpha\,\!$ 未満であれば、帰無仮説（y = 定数）を棄却し、フィットモデルが有意であると判断できます。

切片を固定した場合（たとえば0に固定など）、F検定による p 値は有意な意味を持たず、切片を自由にフィットする通常の重回帰とは異なる解釈となります。

適合度検定表

適合度検定を実行するには、同じX値に対して複数の観測が存在する繰り返しデータ（replicate data）が必要です。つまり、データセット内、または「連結フィット」モードを選択している場合は複数のデータセットにおいて、少なくとも1つのX値が繰り返されている必要があります。

複製データでフィットに使われている表記：

$y_{ij}$ はi番目のX値に対して行われたj番目の測定値
$\bar{y}_{i}$ はi番目のX値における全てのy値の平均
$\hat{y}_{ij}$ はi番目のX値、j番目の測定に対するモデルによる予測値

平方和は、次の通りです。

$RSS=\sum_{i}\sum_{j}(y_{ij}-\hat{y}_{ij})^2$
$LFSS=\sum_{i}\sum_{j}(\bar{y}_{i}-\hat{y}_{ij})^2$
$PESS=\sum_{i}\sum_{j}(y_{ij}-\bar{y}_{i})^2$

線形フィッティングの適合度検定表：

	DF	平方和	平均平方	F値	Prob > F
不適合度	c-k-1	LFSS	MSLF = LFSS / (c - k - 1)	MSLF / MSPE	p値
純誤差	n - c	PESS	MSPE = PESS / (n - c)
誤差	n-k*	RSS

Note:

モデルに切片が含まれている場合、n *= n-1です。それ以外の場合は、n *= nで全体平方和は補正されません。傾きが固定されている場合、 $df_{Model}$ = 0になります。

cはXの異なる値の個数（水準数）を表します。切片を固定している場合、不適合度の自由度はc-kとなります。

共分散行列と相関行列

線形多重回帰（重回帰）の共分散行列は、次のように計算できます。

$Cov(\beta _i,\beta _j)=\sigma ^2(X^{\prime }X)^{-1}$	(28)

任意の2つのパラメータ間の相関係数（ピアソンの相関）は次の式で求められます。

$\rho (\beta _i,\beta _j)=\frac{Cov(\beta _i,\beta _j)}{\sqrt{Cov(\beta _i,\beta _i)}\sqrt{Cov(\beta _j,\beta _j)}}$	(29)

残差解析

$r_i$ は通常の残差 $res_i$ 。

標準化

$r_i^{\prime }=\frac{r_i}s_\varepsilon$	(30)

スチューデント化

内部スチューデント化残差とも呼ばれます。

$r_i^{\prime }=\frac{r_i}{s_\varepsilon\sqrt{1-h_i}}$	(31)

スチューデント化削除

外部スチューデント化残差とも呼ばれる。

$r_i^{\prime }=\frac{r_i}{s_{\varepsilon-i}\sqrt{1-h_i}}$	(32)

スチューデント化、スチューデント化削除での $h_i$ は、行列 $P$ の i 番目の対角成分です。

$P=X(X'X)^{-1}X^{\prime }$	(33)

また、 $s_{\varepsilon-i}$ はi 番目のデータ点を除いて算出された誤差の分散です。

プロット

部分レバレッジプロット

線形多重回帰（重回帰）において、部分レバレッジプロットは、ある独立変数と従属変数との関係を検討するために使用されます。このプロットでは、Yの部分残差を Xまたは切片の部分残差に対してプロットします。ここでの「部分残差」とは、その変数をモデルから除いた状態での回帰残差です。

たとえば、モデル $y=\beta _0+\beta _1x_1+\beta _2x_2\,\!$ を考えます。このとき $x_1\,\!$ に対する部分レバレッジプロットは、 $y=\beta _0+\beta _2x_2\,\!$ の残差回帰を $x_1=\beta _0+\beta _2x_2\,\!$ の残差に対してプロットすることで作成されます。