重回帰分析の計算【多変量解析】

多次元の説明変数と1次元の被説明変数との間に線形関係があると仮定し，そのパラメータの値を推定することを，線形重回帰分析といいます．本稿では最小二乗法に基づく線形重回帰分析の計算法を詳述します．

【スマホでの数式表示について】

当サイトをスマートフォンなど画面幅が狭いデバイスで閲覧すると，数式が画面幅に収まりきらず，正確に表示されない場合があります．その際は画面を回転させ横長表示にするか，ブラウザの表示設定を「PCサイト」にした上でご利用ください．

重回帰分析とは何か

単回帰分析について思い出しておこう．線形単回帰分析(simple linear regression analysis)とは，ひとつの説明変数 $x$ とひとつの被説明変数 $y$ の間に

(1) $\begin{equation*} y = ax + b + \varepsilon \quad \left(\varepsilon \sim {\rm Norm}(0,\sigma^2) \right) \end{equation*}$

なる関係が存在すると仮定したとき，これらの変数 $(x,y)$ に関するデータセット

(2) $\begin{equation*} \{(x_i,y_i)|\; i=1,2,...,n \} \end{equation*}$

を用いて，パラメータ $a,b$ の推定値(estimate) $\hat a , \hat b$ を求めるというものであった．

これに対して，説明変数 $x$ を多次元の変数 $x_1,x_2,...,x_p$ に拡張し，これらと被説明変数 $y$ との間に線形関係

(3) $\begin{equation*} y = a_0 + a_1 x_1 + \cdots + a_j x_j + \cdots + a_p x_p + \varepsilon \quad \left(\varepsilon \sim {\rm Norm}(0,\sigma^2) \right) \end{equation*}$

があると仮定し，そのパラメータ $a_0,a_1,...,a_p$ の値を推定することを，線形重回帰分析(multiple linear regression analysis)という．

線形重回帰モデル

説明変数 $(x_1, ... , x_j, ... , x_p)$ と被説明変数 $y$ の間に，式(3)，すなわち

(4) $\begin{equation*} y = a_0 + \sum_{j=1}^{p} a_j x_j + \varepsilon \end{equation*}$

なる関係があると仮定する．ただし，

(5) $\begin{equation*} \rho_{x_j x_{j'}}=\delta_{jj'}\; (^{\forall}j,j'\in I:=\{1,2,...,p\}) \end{equation*}$

(6) $\begin{equation*} \varepsilon \sim {\rm Norm}(0,\sigma^2) \end{equation*}$

であるとする． $\rho_{x_j x_{j'}}$ は $x_j$ と $x_{j'}$ の相関係数であり，また

(7) $\begin{equation*} \delta_{jj'} = \left\{ \begin{array}{cc} 1 &(j=j')\\ 0 &(j\not=j') \end{array} \right \end{equation*}$

はクロネッカーデルタである．なお，これら2つの条件式(5)，(6)は，それぞれ

式(5)：各説明変数が無相関であること
式(6)：説明変数に含まれない要因の影響は，平均 $0$ ，分散 $\sigma^2$ の正規分布に従う確率変数（揺らぎ）として記述可能であること

を意味している．このとき(3)式は，説明変数と被説明変数の関係を記述するモデルであり，これは線形重回帰モデル(multiple linear regression model)と呼ばれる．また，パラメータと説明変数をそれぞれ

(8) $\begin{equation*} {\bf a} = \left( \begin{array}{c} a_0\\ a_1\\ \vdots \\ a_j \\ \vdots \\ a_p \end{array} \right), \quad {\bf x} = \left( \begin{array}{c} 1\\ x_1\\ \vdots \\ x_j \\ \vdots \\ x_p \end{array} \right) \end{equation*}$

のようにベクトルで書くことにより，(4)式は

(9) $\begin{equation*} y = \;^t{\bf x} \cdot {\bf a} + \varepsilon \end{equation*}$

と表すこともできる．ただし， $\;^t{\bf x}$ は ${\bf x}$ の転置ベクトルである．

線形予測子

線形重回帰モデル(9)は，被説明変数 $y$ の値が，

説明変数の線形結合 $\;^t{\bf x} \cdot {\bf a}$ によって記述される決定論的(deterministic)な要因
$\varepsilon$ によって記述される確率論的(probabilistic)な要因

という2種類の要因の和として理解される，ということを意味する．

線形重回帰モデル(9)において，当初，パラメータ ${\bf a}$ の値は具体的には与えられていない．その後，特定の現象を説明するためのパラメータの値が，何らかの方法で得られたとしよう．これを ${\bf \hat a}$ と書くことにする．パラメータの値が ${\bf \hat a}$ によって具体的に決まったとすると，説明変数 ${\bf x}$ の値が与えられた際，それに対応する被説明変数 $y$ の値を予測する(predict)ことができる． $y$ の値を予測するためには，偶然に決まる確率的要因を考慮する必要はなく，説明変数から必然的に決まる決定論的要因のみ考慮すればよい．

すなわち， $y$ の予測値を $\hat y$ とすれば，これは

(10) $\begin{equation*} \hat y = \;^t{\bf x} \cdot {\bf \hat a} \end{equation*}$

のように書くことができる．つまり，パラメータの値が $\hat a_0,...,\hat a_p$ として決まっていれば，(10)式の説明変数 $x_1,...,x_p$ に数値を代入することにより，被説明変数 $y$ の数値を予測できることを意味する．このことから，(10)式は線形予測子(linear predictor)と呼ばれる．

線形重回帰モデルへの観測データの代入

さて，説明変数と被説明変数についての観測データ（具体的な数値）が，下の表(Table 1)のように得られたとしよう．エクセルシートやデータベース上のテーブルをイメージするとよい．

線形重回帰モデル(4)式あるいは(9)式において，最初の時点ではそのパラメータの値は未知であって，Table 1 のデータセットをよく説明するようなパラメータの値を推定することが，ここで行うべきことである．ただし，データの組の数 $n$ は，未知パラメータの数 $p+1$ よりも十分大きい ( $p+1\ll n$ )とする．

線形重回帰モデルがTable 1 に示される説明変数と被説明変数の間の関係を記述するのに適切だということは，Table 1 のデータを(9)式に代入したとき，(9)式の等式が成り立つことを意味する．(9)式にTable 1 で与えられる $n$ 組のデータを代入することにより，以下のような $n$ 個の式

(11) $\begin{equation*} \left\{ \begin{array}{ccl} y_1 &=& a_0 + x_{11}a_1 + \cdots + x_{1j}a_j + \cdots + x_{1p}a_p + \varepsilon_1 \\ \vdots&& \\ y_i &=& a_0 + x_{i1}a_1 + \cdots + x_{ij}a_j + \cdots + x_{ip}a_p + \varepsilon_i \\ \vdots&& \\ y_n &=& a_0 + x_{n1}a_1 + \cdots + x_{nj}a_j + \cdots + x_{np}a_p + \varepsilon_n \end{array} \right \end{equation*}$

を得る．また，(11)式の変数およびパラメータは，ベクトルや行列を用いて以下のように表記することができる：

(12) $\begin{equation*} {\bf y}= \left( \begin{array}{c} y_1 \\ \vdots \\ y_i \\ \vdots \\ y_n \end{array} \right), \quad {\bf x}_i= \left( \begin{array}{c} 1 \\ x_{i1} \\ \vdots \\ x_{ij} \\ \vdots \\ x_{ip} \end{array} \right), \quad {\bf a} = \left( \begin{array}{c} a_0\\ a_1\\ \vdots \\ a_j \\ \vdots \\ a_p \end{array} \right), \quad {\boldsymbol \epsilon} = \left( \begin{array}{c} \varepsilon_1 \\ \vdots \\ \varepsilon_i \\ \vdots \\ \varepsilon_n \end{array} \right), \end{equation*}$

(13) $\begin{equation*} X= \left( \begin{array}{cccccc} 1 & x_{11} & \cdots & x_{1j} & \cdots & x_{1p}\\ \vdots & \vdots && \vdots && \vdots\\ 1 & x_{i1} & \cdots & x_{ij} & \cdots & x_{ip}\\ \vdots & \vdots && \vdots && \vdots\\ 1 & x_{n1} & \cdots & x_{nj} & \cdots & x_{np}\\ \end{array} \right) = \left( \begin{array}{c} ^t{\bf x}_1\\ \vdots \\ ^t{\bf x}_i\\ \vdots\\ ^t{\bf x}_n\\ \end{array} \right). \end{equation*}$

なお，拡大行列 $[{\bf y},X]$ は，値1が並ぶ2列目（すなわち $X$ の1列目）を除いて，Table 1 に他ならない．これらの表記を用いれば，(11)式は

(14) $\begin{equation*} {\bf y}= \left( \begin{array}{c} y_1 \\ \vdots \\ y_i \\ \vdots \\ y_n \end{array} \right) = \left( \begin{array}{c} ^t{\bf x}_1 \cdot {\bf a} + \varepsilon_1\\ \vdots \\ ^t{\bf x}_i \cdot {\bf a} + \varepsilon_i\\ \vdots\\ ^t{\bf x}_n \cdot {\bf a} + \varepsilon_n\\ \end{array} \right) = \left( \begin{array}{c} ^t{\bf x}_1\\ \vdots \\ ^t{\bf x}_i\\ \vdots\\ ^t{\bf x}_n\\ \end{array} \right) {\bf a} + \left( \begin{array}{c} \varepsilon_1 \\ \vdots \\ \varepsilon_i \\ \vdots \\ \varepsilon_n \end{array} \right) =X{\bf a} + {\boldsymbol \epsilon} \end{equation*}$

より，

(15) $\begin{equation*} {\bf y}=X{\bf a} + {\boldsymbol \epsilon} \end{equation*}$

のように書くことができる．結局，線形重回帰モデル(9)式に観測データセットを代入することで
(25)式が得られる．

線形重回帰パラメータの最小二乗推定

観測データを代入された線形重回帰モデル(25)式において，ベクトル ${\bf a}$ およびベクトル ${\boldsymbol \epsilon}$ の値は未知である．我々の目的はパラメータ ${\bf a}$ の値を決めることであるが，観測データ $[{\bf y},X]$ やモデル(9)だけからこれを決定することはできない．そこで，

線形予測子(10)によって得られる被説明変数の予測値と実際の観測で得られている被説明変数の値の差ができるだけ小さくなるようにパラメータ ${\bf \hat a}$ の値を選ぶ

という方針を立て，これによって線形重回帰モデルのパラメータの値を決定することにする．これを実現するための手法として最小二乗法(least squares method)がある．

まず，被説明変数の予測値 $\hat y_i$ とその実際の観測値 $y_i$ の差は残差(residual)と呼ばれ，

(16) $\begin{equation*} r_i := y_i - \hat y_i \end{equation*}$

で定義される．また，観測データ全体で残差の絶対値を小さくするために残差平方和(residual sum of square; RSS)

(17) $\begin{equation*} R := \sum_{i=1}^{n} r_i^2 = \sum_{i=1}^{n} (y_i - \hat y_i)^2 \end{equation*}$

を定義し，これを極小化するパラメータ ${\bf \hat a}$ を求める．線形予測子

(18) $\begin{equation*} \hat y_i = \;^t{\bf x}_i\cdot {\bf \hat a} = \hat a_0 + x_{i1}\hat a_1 + \cdots + x_{ij}\hat a_j + \cdots + x_{ip}\hat a_p \end{equation*}$

を考慮すれば， $R$ は ${\bf \hat a}$ の関数であることが分かる．さらに， $R({\bf \hat a})$ は各 $\hat a_j$ の2次関数となり，その係数は $\hat a_0$ の係数が $n$ ， $j\not=0$ なる $\hat a_j$ の係数が $(\sum_{i=1}^{n} x_{ij}^2)$ であることから，特に $R({\bf \hat a})$ は各 $\hat a_j$ に関する下に凸の2次関数となることが分かる．すなわち， $R$ の極値条件

(19) $\begin{equation*} \frac{\partial R({\bf \hat a})}{\partial {\bf \hat a} }={\bf 0} \end{equation*}$

に基づいて端点 ${\bf \hat a}$ を求めれば，それが残差平方和 $R$ を極小化するような ${\bf a}$ の値である．(19)式左辺の計算は，微分演算子

(20) $\begin{equation*} \frac{\partial }{\partial {\bf \hat a} } =\left(\frac{\partial }{\partial \hat a_0 }, \frac{\partial }{\partial \hat a_1},..., \frac{\partial }{\partial \hat a_j},..., \frac{\partial }{\partial \hat a_p}\right) \end{equation*}$

および線形予測子(18)に注意しながら， $a$ の添字 $j$ について $j=0$ のときと $j\not=0$ のときとで場合分けをすればよい．

[ $j=0$ のとき]

(21) $\begin{eqnarray*} \frac{\partial R({\bf \hat a})}{\partial \hat a_0} &=& \frac{\partial}{\partial \hat a_0} \sum_{i=1}^{n} (y_i - \;^t{\bf x}_i\cdot {\bf \hat a})^2 \nonumber \\ &=& \frac{\partial}{\partial \hat a_0} \sum_{i=1}^{n} \left\{ y_i - (\hat a_0 + x_{i1}\hat a_1 + \cdots + x_{ij}\hat a_j + \cdots + x_{ip}\hat a_p )\right\}^2 \nonumber \\ &=& \sum_{i=1}^{n} \frac{\partial}{\partial \hat a_0} \left\{y_i - (\hat a_0 + x_{i1}\hat a_1 + \cdots + x_{ij}\hat a_j + \cdots + x_{ip}\hat a_p )\right\}^2 \nonumber \\ &=& \sum_{i=1}^{n} 2 \left\{y_i - (\hat a_0 + x_{i1}\hat a_1 + \cdots + x_{ij}\hat a_j + \cdots + x_{ip}\hat a_p )\right\}\cdot(-1) \nonumber \\ &=& 2 \left\{ - \sum_{i=1}^{n} y_i + \hat a_0 n + \hat a_1 \sum_{i=1}^{n} x_{i1} + \cdots + \hat a_j \sum_{i=1}^{n} x_{ij} + \cdots + \hat a_p \sum_{i=1}^{n} x_{ip} \right\} =0 \end{eqnarray*}$

したがって

(22) $\begin{equation*} \hat a_0 n+ \hat a_1 \sum_{i=1}^{n} x_{i1} + \cdots + \hat a_j \sum_{i=1}^{n} x_{ij} + \cdots + \hat a_p \sum_{i=1}^{n} x_{ip} = \sum_{i=1}^{n} y_i \end{equation*}$

を得る．

[ $j\not=0$ のとき]

(23) $\begin{eqnarray*} \frac{\partial R({\bf \hat a})}{\partial \hat a_j} &=& \frac{\partial}{\partial \hat a_j} \sum_{i=1}^{n} (y_i - \;^t{\bf x}_i\cdot {\bf \hat a})^2 \nonumber \\ &=& \frac{\partial}{\partial \hat a_j} \sum_{i=1}^{n} \left\{ y_i - (\hat a_0 + x_{i1}\hat a_1 + \cdots + x_{ij}\hat a_j + \cdots + x_{ip}\hat a_p )\right\}^2 \nonumber \\ &=& \sum_{i=1}^{n} \frac{\partial}{\partial \hat a_j} \left\{y_i - (\hat a_0 + x_{i1}\hat a_1 + \cdots + x_{ij}\hat a_j + \cdots + x_{ip}\hat a_p )\right\}^2 \nonumber \\ &=& \sum_{i=1}^{n} 2 \left\{y_i - (\hat a_0 + x_{i1}\hat a_1 + \cdots + x_{ij}\hat a_j + \cdots + x_{ip}\hat a_p )\right\}\cdot(-x_{ij}) \nonumber \\ &=& 2 \left\{ - \sum_{i=1}^{n} x_{ij} y_i + \hat a_0 \sum_{i=1}^{n} x_{ij} + \hat a_1 \sum_{i=1}^{n} x_{ij}x_{i1} + \cdots + \hat a_j \sum_{i=1}^{n} x_{ij}x_{ij} + \cdots + \hat a_p \sum_{i=1}^{n} x_{ij}x_{ip} \right\} =0 \nonumber \\ \end{eqnarray*}$

したがって

(24) $\begin{equation*} \hat a_0 \sum_{i=1}^{n} x_{ij} + \hat a_1 \sum_{i=1}^{n} x_{ij}x_{i1} + \cdots + \hat a_j \sum_{i=1}^{n} x_{ij}x_{ij} + \cdots + \hat a_p \sum_{i=1}^{n} x_{ij}x_{ip} = \sum_{i=1}^{n} x_{ij} y_i \end{equation*}$

を得る．

(22)式および(24)式より，求めるパラメータ ${\bf \hat a}=\;^t(\hat a_0, \hat a_1,...,\hat a_j,...,\hat a_p)$ の値は，次の $(p+1)$ 元1次連立方程式を解くことにより得られる：

(25) $\begin{eqnarray*} {}&& \hat a_0 n+ \hat a_1 \sum_{i=1}^{n} x_{i1} + \cdots + \hat a_j \sum_{i=1}^{n} x_{ij} + \cdots + \hat a_p \sum_{i=1}^{n} x_{ip} = \sum_{i=1}^{n} y_i \nonumber \\ {}&& \hat a_0 \sum_{i=1}^{n} x_{i1} + \hat a_1 \sum_{i=1}^{n} x_{i1}x_{i1} + \cdots + \hat a_j \sum_{i=1}^{n} x_{i1}x_{ij} + \cdots + \hat a_p \sum_{i=1}^{n} x_{i1}x_{ip} = \sum_{i=1}^{n} x_{i1} y_i \nonumber \\ {}&& \vdots \nonumber \\ {}&& \hat a_0 \sum_{i=1}^{n} x_{ij} + \hat a_1 \sum_{i=1}^{n} x_{ij}x_{i1} + \cdots + \hat a_j \sum_{i=1}^{n} x_{ij}x_{ij} + \cdots + \hat a_p \sum_{i=1}^{n} x_{ij}x_{ip} = \sum_{i=1}^{n} x_{ij} y_i \nonumber \\ {}&& \vdots \nonumber \\ {}&& \hat a_0 \sum_{i=1}^{n} x_{ip} + \hat a_1 \sum_{i=1}^{n} x_{ip}x_{i1} + \cdots + \hat a_j \sum_{i=1}^{n} x_{ip}x_{ij} + \cdots + \hat a_p \sum_{i=1}^{n} x_{ip}x_{ip} = \sum_{i=1}^{n} x_{ip} y_i. \end{eqnarray*}$

変数のベクトル・行列表記(12)および(13)を用いれば，連立方程式(25)は

(26) $\begin{equation*} \;^tXX{\bf \hat a} = \;^tX{\bf y} \end{equation*}$

となっていることがわかる．この左辺にある $(p+1)\times(p+1)$ 正方行列

(27) $\begin{equation*} \;^tXX= \left( \begin{array}{cccccc} n & \sum_i x_{i1} & \cdots & \sum_i x_{ij} & \cdots & \sum_i x_{ip}\\ \sum_i x_{i1} & \sum_i x_{i1}x_{i1} & \cdots & \sum_i x_{i1}x_{ij} & \cdots & \sum_i x_{i1}x_{ip}\\ \vdots&\vdots&&\vdots&&\vdots\\ \sum_i x_{ij} & \sum_i x_{ij}x_{i1} & \cdots & \sum_i x_{ij}x_{ij} & \cdots & \sum_i x_{ij}x_{ip}\\ \vdots&\vdots&&\vdots&&\vdots\\ \sum_i x_{ip} & \sum_i x_{ip}x_{i1} & \cdots & \sum_i x_{ip}x_{ij} & \cdots & \sum_i x_{ip}x_{ip}\\ \end{array} \right) \end{equation*}$

は対称行列である．

(26)式の両辺に左から $\left(\;^tXX\right)^{-1}$ を掛ければ

(28) $\begin{equation*} \left(\;^tXX\right)^{-1}\left(\;^tXX\right){\bf \hat a} = \left(\;^tXX\right)^{-1}\;^tX{\bf y}. \end{equation*}$

これにより，パラメータ ${\bf a}$ の最小二乗推定値 ${\bf \hat a}$ が

(29) $\begin{equation*} {\bf \hat a} = \left(\;^tXX\right)^{-1}\;^tX{\bf y} \end{equation*}$

のように求まる．

Remark: 残差について

統計学における一般的な術語の用法としては，統計的誤差(statistical error)または揺らぎ(disturbance)とは，ある観測量とその期待値との差のことであり，残差(residual)とは直接観測できない統計的誤差の推定値のことである．ここでは，重回帰モデルにおける揺らぎ $\varepsilon$ は直接観測できない量であり，残差 $r$ が揺らぎ $\varepsilon$ の推定値ということになる．

記事一覧

カテゴリー

重回帰分析の計算【多変量解析】