こんにちは、小郷です。サケやマスでは、川を下り海で成長した後、生まれた川に帰ってきて産卵して命を終える母川回帰と呼ばれる現象が知られています。(某色塗りゲームで有名な)サーモンランという名前でも知られています。
主に匂い物質を嗅ぎ分けて故郷に帰ることができると複数の研究から示唆されていますが、視覚や海流、地磁気など複数の説も唱えられています。これらの全て/一部の要素が複雑に絡み合っているのかもしれません。まるで重回帰モデルのようですね。
前回は説明変数が一つの場合の単回帰分析の基礎について解説しました。今回は二つ以上の説明変数を持つ回帰である、重回帰分析について解説します。
今回からLaTeXを書けるようにしてみました。
重回帰分析の基礎理論
まずは線形単回帰モデルの一般式をおさらいします。
$$ y_i = \beta _0 + \beta _{1}x_i + \epsilon_i $$
それぞれ以下を表します。
\( \beta _0\) : 切片
\( \beta _1 \) : 傾き
\( x_i \) : 説明変数
\( y_i \) : 応答変数
\( \epsilon \) : 誤差項
これに対し、線形重回帰モデルは以下の式で表すことができます。
$$ y_i = \beta _0 + \beta _{1}x_{1i} + \beta _{2}x_{2i} + … + \beta_{p}x_{pi} + \epsilon_i $$
家賃を予想するケースを考えてみましょう。家賃を左右する要素として「都道府県」「築年数」「駅からの距離」「部屋の大きさ」「主要駅からの距離」「告知事項の有無」などがありますが、これらは説明変数として(\(x_1, x_2, …, x_p\))に該当します。
対して(\(\beta _0, \beta _1, \beta_2, …, \beta _p \))は偏回帰係数と呼ばれます。ある説明変数(例えば、\(x_1\))以外のを固定した場合に、その変数を1単位変化させた時にどのくらいの重み(\(\beta _1\))で応答変数が変動するかを表す係数と解釈できます。例えば「都道府県」と「最寄駅からの距離」では、前者の方が家賃決定への寄与が大きいことは、感覚的になんとなくわかると思います。この寄与度を決める値が偏回帰係数です。
ただしこの値は、他の説明変数の影響を受けるなど、簡単に求められる値ではありません。詳しくみていきましょう。かなり丁寧に解説していると自負しています。
最小二乗法による導出(簡易版)
単回帰と同様に、最小二乗法で求めることができます。かなり長いので、数学的な話は興味ないよ！という方は「偏回帰係数の値は、説明変数と従属変数の分散と共分散で決まる」事を念頭に置いて、この項目はすっ飛ばしてしまってください。
今回は説明変数を2つにして、シンプルに考えてみます。以下の残差平方和が最小になる\(\beta _0, \beta _1, \beta_2\)を求めます。
$$ RSS = \sum^n_{i=1} (y – y_i)^2 = \sum^n_{i=1} \{y_i – (\beta_0 + \beta_1x_{1i} + \beta_2x_{2i}) \}^2 $$
\(\beta _0, \beta _1, \beta_2\)それぞれについて、偏微分を行います。この偏微分を行った結果が0になる時、RSSも最小になります。
\(\beta _0\)についての偏微分
\(\beta_0\)は切片なので、微分しても1にしかなりません。
$$ \frac{\partial RSS}{\partial \beta_0} = \frac{\partial \sum^n_{i=1} \{y_i – (\beta_0 + \beta_1x_{1i} + \beta_2x_{2i}) \}^2}{\partial \beta_0} = 0 $$
$$ -2\sum^n_{i=1} (y_i – (\beta_0 + \beta_1x_{1i} + \beta_2x_{2i}) ) = 0 $$
-2で両片を割って削除して綺麗にする他、\(\beta_0\)はn回足されるだけなので、\(\sum\)の外に出せます。符号も整理しましょう。以降も同様に進めていきます。
$$ n\beta_0 + \beta_1 \sum^n_{i=1} x_{1i} + \beta_2 \sum^n_{i=1} x_{2i} = \sum^n_{i=1} y_i $$
\(\beta _1\)についての偏微分
$$ \frac{\partial RSS}{\partial \beta_1} = \frac{\partial \sum^n_{i=1} \{y_i – (\beta_0 + \beta_1x_{1i} + \beta_2x_{2i}) \}^2}{\partial \beta_1} = 0 $$
$$ \beta_0 \sum^n_{i=1} x_{1i} + \beta_1 \sum^n_{i=1} x_{1i}^2 + \beta_2 \sum^n_{i=1} x_{1i}x_{2i} = \sum^n_{i=1} x_{1i}y_i $$
\(\beta _2\)についての偏微分
$$ \frac{\partial RSS}{\partial \beta_2} = \frac{\partial \sum^n_{i=1} \{y_i – (\beta_0 + \beta_1x_{1i} + \beta_2x_{2i}) \}^2}{\partial \beta_2} = 0 $$
$$ \beta_0 \sum^n_{i=1} x_{2i} + \beta_1 \sum^n_{i=1} x_{1i}x_{2i} + \beta_2 \sum^n_{i=1} x_{2i}^2 = \sum^n_{i=1} x_{2i}y_i $$
連立方程式を解く
ここまでに得られた上記の偏微分の式(正規方程式)を、連立方程式として解いていきます。まずは最も簡単なので、\(\beta_0\)の値を求めてみましょう。\(\beta_0\)を左辺におき、全体をnで割ると、\(\sum^n_{i=1} y_i\)、\(\sum^n_{i=1} x_{1i}\)、\(\sum^n_{i=1} x_{2i}\)はそれぞれ、観測値の合計であることから、nで割ることで平均値が求められます。
$$ \beta_0 = \bar{y} – \beta_1\bar{x_1} – \beta_2\bar{x_2} $$
\(\beta_1\)の正規方程式と、\(\beta_2\)の正規方程式はちょっと面倒です。
まずは\(x_i1\)と\(x_i2\)と\(y_i\)を、以下のように変換します。
$$ x_1i = (x_1i – \bar{x_1}) + \bar{x_1} $$
$$ x_2i = (x_2i – \bar{x_2}) + \bar{x_2} $$
$$ y_i = (y_i – \bar{y}) + \bar{y} $$
また、以下を適用します。n個説明変数の合計 = 説明変数の平均×n、という以上の意味はありません。
$$ \sum^n_{i=1} x_1i = n\bar{x_1} $$
$$ \sum^n_{i=2} x_2i = n\bar{x_2} $$
まずは\(\beta_1\)の正規方程式を見ていきます。
$$ n\beta_0\bar{x_1} + \beta_1 \sum^n_{i=1} ((x_1i – \bar{x_1}) + \bar{x_1})^2 + \beta_2 \sum^n_{i=1} (x_1i – \bar{x_1} + \bar{x_1})(x_2i – \bar{x_2} + \bar{x_2}) = \sum^n_{i=1} (x_1i – \bar{x_1} + \bar{x_1})(y_i – \bar{y} + \bar{y}) $$
部分的に計算します。\(\beta_1\)の項については、以下のように展開できます。
$$ \beta_1 \sum^n_{i=1}((x_1i – \bar{x})^2 +2(x_1i – \bar{x})\bar{x} + \bar{x}^2 )$$
$$ \beta_1 (\sum^n_{i=1}(x_1i – \bar{x})^2 +\sum^n_{i=1}2(x_1i – \bar{x})\bar{x} + \sum^n_{i=1}\bar{x}^2 )$$
ここで、偏差の合計は0になることから、真ん中の式は削除することができます。
$$ \beta_1\sum^n_{i=1}(x_1i – \bar{x})^2 + \beta_1\sum^n_{i=1}\bar{x}^2 $$
\(\beta_2\)の項については、以下のように展開できます。
$$ \beta_2 \left( \sum^n_{i=1} (x_{1i} – \bar{x_1})(x_{2i} – \bar{x_2}) + \bar{x_2} \sum^n_{i=1} (x_{1i} – \bar{x_1}) + \bar{x_1} \sum^n_{i=1} (x_{2i} – \bar{x_2}) + n \bar{x_1} \bar{x_2} \right)$$
ここでも偏差の合計を消すことができるので、以下のようにシンプルになります。
$$ \beta_2\sum^n_{i=1} (x_{1i} – \bar{x_1})(x_{2i} – \bar{x_2}) + n\beta_2\bar{x_1} \bar{x_2} $$
右辺も同様に展開します。
$$ \sum^n_{i=1} (x_{1i} – \bar{x_1})(y_i – \bar{y}) + \bar{y} \sum^n_{i=1} (x_{1i} – \bar{x_1}) + \bar{x_1} \sum^n_{i=1} (y_i – \bar{y}) + n \bar{x_1} \bar{y} $$
$$ \sum^n_{i=1} (x_{1i} – \bar{x_1})(y_i – \bar{y}) + n \bar{x_1} \bar{y} $$
最終的な式は以下のようになります。
$$ n\beta_0 \bar{x_1} + \beta_1 \left( \sum^n_{i=1} (x_{1i} – \bar{x_1})^2 + n \bar{x_1}^2 \right) + \beta_2 \left( \sum^n_{i=1} (x_{1i} – \bar{x_1})(x_{2i} – \bar{x_2}) + n \bar{x_1} \bar{x_2} \right)
= \sum^n_{i=1} (x_{1i} – \bar{x_1})(y_i – \bar{y}) + n \bar{x_1} \bar{y} $$
なお、\(\sum^n_{i=1} (x_{1i} – \bar{x_1})^2\)をnで割ったものは\(x_1\)の分散、\(\sum^n_{i=1} (x_{1i} – \bar{x_1})(x_{2i} – \bar{x_2})\)をnで割ったものは\(x_1\)と\(x_2\)の共分散、 \(\sum^n_{i=1} (x_{1i} – \bar{x_1})(y_i – \bar{y})\)をnで割ったものは\(x_1\)とyの共分散になるので、両辺をnで割ることで、以下のように省略形で記述することもできます。
$$ \beta_0 \bar{x_1} + \beta_1 \left( Var(x_1) + \bar{x_1}^2 \right) + \beta_2 \left( Cov(x_1, x_2) + \bar{x_1} \bar{x_2} \right) = Cov(x_1, y) + \bar{x_1} \bar{y} $$
\(\beta_2\)の正規方程式については展開は省きますが、最終的に以下のようになります。
$$ \beta_0 \bar{x_2} + \beta_1 \left( Cov(x_1, x_2) + \bar{x_1}\bar{x_2} \right) + \beta_2 \left( Var(x_2) + \bar{x_2}^2 \right) = Cov(x_2, y) + \bar{x_2}\bar{y} $$
\(\beta_1\)と\(\beta_2\)を求める
かなり複雑な計算をしてきましたが、あともう少しです。
先に求めた二種類の分散と共分散と平均の式に以下の\(\beta_0\)を代入すると、最終的に分散と共分散だけの式が作られ、これを連立方程式とします。
$$ \beta_0 = \bar{y} – \beta_1\bar{x_1} – \beta_2\bar{x_2} $$
$$ \begin{cases} \beta_1 \text{Var}(x_1) + \beta_2 \text{Cov}(x_1, x_2) = \text{Cov}(x_1, y) \\ \beta_1 \text{Cov}(x_1, x_2) + \beta_2 \text{Var}(x_2) = \text{Cov}(x_2, y) \\ \end{cases} $$
こちらを簡単に解いて、\(\beta_1\)と\(\beta_2\)を得るためには、行列に変換することが有効です。
$$ \begin{pmatrix} \text{Var}(x_1) & \text{Cov}(x_1, x_2) \\ \text{Cov}(x_1, x_2) & \text{Var}(x_2) \\ \end{pmatrix} \begin{pmatrix} \beta_1 \\ \beta_2 \\ \end{pmatrix} = \begin{pmatrix} \text{Cov}(x_1, y) \\ \text{Cov}(x_2, y) \end{pmatrix} $$
左辺を\(\beta_1\)と\(\beta_2\)のみにするために、右辺に逆行列を掛けます。
$$ \begin{pmatrix} \beta_1 \\ \beta_2 \end{pmatrix} = \frac{1}{\text{Var}(x_1) \text{Var}(x_2) – \text{Cov}(x_1, x_2)^2} \begin{pmatrix} \text{Var}(x_2) & – \text{Cov}(x_1, x_2) \\ -\text{Cov}(x_1, x_2) & \text{Var}(x_1) \end{pmatrix} \begin{pmatrix} \text{Cov}(x_1, y) \\ \text{Cov}(x_2, y) \end{pmatrix} $$
最終的に以下の式が求められます。
$$ \beta_1 = \frac{\text{Var}(x_2) \text{Cov}(x_1, y) – \text{Cov}(x_1, x_2) \text{Cov}(x_2, y)}{\text{Var}(x_1) \text{Var}(x_2) – \text{Cov}(x_1, x_2)^2} $$
$$ \beta_2 = \frac{\text{Var}(x_1) \text{Cov}(x_2, y) – \text{Cov}(x_1, x_2) \text{Cov}(x_1, y)}{\text{Var}(x_1) \text{Var}(x_2) – \text{Cov}(x_1, x_2)^2} $$
多重共線性と共分散
ここで共分散の性質について改めて説明します。繰り返しになりますが、数式では以下のように表されます。
$$ \text{Cov}(x_1, x_2) = \frac{1}{n} \sum_{i=1}^{n} (x_{1i} – \bar{x}_1)(x_{2i} – \bar{x}_2) $$
この共分散の値は、その大きさによって以下の性質を持ちます。
正の値で大きい場合
\(x_1\)と\(x_2\)の間に、強い正の相関があることを示します。\(x_{1i}\)が大きい時に\(x_{2i}\)も大きくなる、あるいは\(x_{1i}\)が小さい時に\(x_{2i}\)も小さくなるならば、\((x_{1i} – \bar{x}_1)(x_{2i} – \bar{x}_2)\)の値は正の値を取りやすくなるので、結果として共分散の値は正の値で0から遠ざかります。
負の値で大きい場合
\(x_1\)と\(x_2\)の間に、強い負の相関があることを示します。\(x_{1i}\)が大きい時に、\(x_{2i}\)が小さくなる、あるいは\(x_{1i}\)が小さい時に、\(x_{2i}\)が大きくなるならば、\((x_{1i} – \bar{x}_1)(x_{2i} – \bar{x}_2)\)の値は負の値を取りやすくなるので、結果として共分散の値は負の値で0から遠ざかります。
0に近い場合
\(x_1\)と\(x_2\)の間の相関が弱い事を表します。\(x_{1i}\)が大きいとしても\(x_{2i}\)の値が大きくなるとは限らない、あるいは\(x_{1i}\)が小さいとしても\(x_{2i}\)の値が小さくなるとは限らないので、\((x_{1i} – \bar{x}_1)(x_{2i} – \bar{x}_2)\)の値は正にも負にもランダムになりやすく、結果としてその和は0に近づくと考えられます。
複数の説明変数間に強い相関があるために、それぞれの説明変数が従属変数に与える影響が評価しにくくなることを、多重共線性と言います。
多重共線性が及ぼす影響の数学的解説
\(\beta_1\)を導出する数式を例に解説します。
$$ \beta_1 = \frac{\text{Var}(x_2) \text{Cov}(x_1, y) – \text{Cov}(x_1, x_2) \text{Cov}(x_2, y)}{\text{Var}(x_1) \text{Var}(x_2) – \text{Cov}(x_1, x_2)^2} $$
\(\text{Cov}(x_1, x_2)\)が0から遠ざかる場合、\(\text{Cov}(x_1, x_2)^2\)の値が大きくなるため、\(\beta_1\)の分母が小さくなります。よって、少しの値の変動で\(\beta_1\)の値が変わってしまうことになり、\(\beta_1\)は非常に不安定になります。
\(x_1\)と\(x_2\)のいずれかが説明変数として正しくないと判断できるため、どちらかをモデルに含めない(次元削減)など、説明変数の選択を再度実施する必要があります。VIFの話は一旦ここでは省略しますが、VIFもとい\(R^2\)は分散と共分散から得られる値です。
その他の影響
多重共線性の話からは少し外れますが、\(\text{Cov}(x_1, y)\)が何を意味するかを簡単に解説します。
先述の\(\beta_1\)の値の分子に着目をすると、\(\text{Var}(x_2)\)\(\text{Cov}(x_1, y)\)が大きく、\(\text{Cov}(x_2, y)\)\(\text{Cov}(x_1, x_2)^2\)の値が小さいほど、\(\beta_1\)の値は大きくなります。これらの値は、以下を意味します。
- \(\text{Var}(x_2)\): \(\beta_1\)に、\(x_1\)以外の説明変数の影響を考慮している。
- \(\text{Cov}(x_1, y)\): \(x_1\)と\(y\)の相関の大きさを表す。
- \(\text{Cov}(x_2, y)\): \(x_2\)と\(y\)の相関の大きさを表す。
- \(\text{Cov}(x_1, x_2)^2\): \(x_1\)と\(x_2\)の相関の大きさを表す。
つまり、\(\text{Cov}(x_1, y)\)が小さければ、説明変数\(x_1\)は従属変数\(y\)の説明変数としてウェイトが低いと言えるため、\(\beta_1\)の値は小さくなります。
また、\(\text{Cov}(x_2, y)\)の値が大きい場合には\(\beta_1\)が小さくなるため、説明変数\(x_1\)は相対的に、説明変数\(x_2\)よりも従属変数\(y\)の説明変数としてウェイトが低くなります。
終わりに
ここまでで、重回帰モデルという強力なツールについて解説をしてきました。基本的にこんな複雑な計算は、コンピュータに任せるケースが大半だと思います。ただ、中身を知っているとモデル構築時のトラブルシューティングにも役立つのではないかと思うので、どこかで役立てていただけると、書いた奴が喜びます。
ところで、疲れた時に海に行きたくなったり、暖かくて狭くて暗い場所で小さくなったりしたくなるのは、胎内回帰願望というのだそうです。個人的には魂のルフランを想起します。有から無に帰るような極端なことはせずとも、たまには全てを忘れてリラックスしたいものです。
生と死・有と無のような二値に対する推定は、ロジスティック回帰と呼ばれるモデルで行うことができます。次回PART3では、このモデルについて解説したいと思います。
以降はおまけです。
最小二乗法による導出(一般式)
先に簡易版の解説をした偏回帰係数の導出を一般化しています。
以下の残差平方和を最小にする\(\beta _0, \beta _1, \beta_2, …, \beta _p\)を求めます。
$$ RSS = \sum^n_{i=1} \{y_i – (\hat{\beta_0} + \hat{\beta_1}x_{1i} + … + \hat{\beta_p}x_{pi}) \}^2 $$
最初のステップとして、上記の\(\beta _0, \beta _1, \beta_2, …, \beta _p\)のそれぞれで偏微分します。
\(j = 1, 2, … , p\)として、この結果が0になるとすると(正規方程式と呼びます)、以下のように置くことができます。
$$ \frac{\partial RSS}{\partial \beta_j} = \frac{ \partial \sum^n_{i=1} \{y_i – (\hat{\beta_0} + \hat{\beta_1}x_{1i} + … + \hat{\beta_p}x_{pi}) \}^2}{\partial \beta_j} = 0 $$
\(y_i – (\hat{\beta_0} + \hat{\beta_1}x_{1i} + … + \hat{\beta_p}x_{pi})\)の部分を\(\beta_j\)で微分すると、他の項は全て消えるので、\(-x_{ji}\)となります。つまり、以下の式が求められます。
$$ \frac{\partial RSS}{\partial \beta_j} = \sum^n_{i=1} -2x_{ji}\{y_i – (\hat{\beta_0} + \hat{\beta_1}x_{1i} + … + \hat{\beta_p}x_{pi}) \} = 0 $$
よりシンプルな式に変形します。。
$$ \sum^n_{i=1} x_{ji}\{y_i – (\hat{\beta_0} + \hat{\beta_1}x_{1i} + … + \hat{\beta_p}x_{pi}) \} = 0 $$
\(\beta _0, \beta _1, \beta_2, …, \beta _p\)について偏微分すると、以下のp元連立方程式が得られます。
$$ \left\{\begin{align*} &\sum^n_{i=1} \{y_i – (\hat{\beta_0} + \hat{\beta_1}x_{1i} + … + \hat{\beta_p}x_{pi}) \} = 0 \\ &\sum^n_{i=1} x_{1i}\{y_i – (\hat{\beta_0} + \hat{\beta_1}x_{1i} + … + \hat{\beta_p}x_{pi}) \} = 0 \\ &\sum^n_{i=1} x_{2i}\{y_i – (\hat{\beta_0} + \hat{\beta_1}x_{1i} + … + \hat{\beta_p}x_{pi}) \} = 0 \\& ・\\& ・\\& ・\\ &\sum^n_{i=1} x_{pi}\{y_i – (\hat{\beta_0} + \hat{\beta_1}x_{1i} + … + \hat{\beta_p}x_{pi}) \} = 0 \end{align*}\right.$$
証明は省きますが、最終的に偏回帰係数は以下の式で求められます。
$$ \hat{\beta}_j = \frac{\text{Var}(x_k) \text{Cov}(x_j, y) – \sum_{k \neq j} \text{Cov}(x_j, x_k) \text{Cov}(x_k, y)}{\prod_{k \neq j} \text{Var}(x_k) – \sum_{k \neq j} \text{Cov}(x_j, x_k)^2} $$
まぁこんなの面倒くさすぎるので、現代に生きる我々はコンピュータにやらせりゃいいですね。
