回帰分析はかく語りき Part2 重回帰分析

2024-06-27 10:54Ogo実務Tips_運用 Playbook

こんにちは、小郷です。サケやマスでは、川を下り海で成長した後、生まれた川に帰ってきて産卵して命を終える母川回帰と呼ばれる現象が知られています。(某色塗りゲームで有名な)サーモンランという名前でも知られています。
主に匂い物質を嗅ぎ分けて故郷に帰ることができると複数の研究から示唆されていますが、視覚や海流、地磁気など複数の説も唱えられています。これらの全て/一部の要素が複雑に絡み合っているのかもしれません。まるで重回帰モデルのようですね。

前回は説明変数が一つの場合の単回帰分析の基礎について解説しました。今回は二つ以上の説明変数を持つ回帰である、重回帰分析について解説します。本記事では基礎理論と、多重共線性についての解説に重点を置いてみます。
今回からLaTeXを書けるようにしてみました。

重回帰分析の基礎理論

まずは線形単回帰モデルの一般式をおさらいします。

$$ y_i = \beta _0 + \beta _{1}x_i + \epsilon_i $$

それぞれ以下を表します。

β₀ : 切片
β₁ : 傾き
xᵢ : 説明変数
yᵢ : 応答変数
ϵ : 誤差項

これに対し、線形重回帰モデルは以下の式で表すことができます。

$$ y_i = \beta _0 + \beta _{1}x_{1i} + \beta _{2}x_{2i} + … + \beta_{p}x_{pi} + \epsilon_i $$

家賃を予想するケースを考えてみましょう。家賃を左右する要素として「都道府県」「築年数」「駅からの距離」「部屋の大きさ」「主要駅からの距離」「告知事項の有無」などがありますが、これらは説明変数として(x₁,x₂,…,xₚ)に該当します。
対して(β₀,β₁,β₂,…,βₚ)は偏回帰係数と呼ばれます。ある説明変数(例えば、x₁)以外のを固定した場合に、その変数を1単位変化させた時にどのくらいの重み(β₁)で応答変数が変動するかを表す係数と解釈できます。例えば「都道府県」と「最寄駅からの距離」では、前者の方が家賃決定への寄与が大きいことは、感覚的になんとなくわかると思います。この寄与度を決める値が偏回帰係数です。
ただしこの値は、他の説明変数の影響を受けるなど、簡単に求められる値ではありません。詳しくみていきましょう。かなり丁寧に解説していると自負しています。

最小二乗法による導出(簡易版)

単回帰と同様に、最小二乗法で求めることができます。かなり長いので、数学的な話は興味ないよ！という方は「偏回帰係数の値は、説明変数と従属変数の分散と共分散で決まる」事を念頭に置いて、この項目はすっ飛ばしてしまってください。

今回は説明変数を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$

β₀,β₁,β₂それぞれについて、偏微分を行います。この偏微分を行った結果が0になる時、RSSも最小になります。

β₀についての偏微分

β₀は切片なので、微分しても1にしかなりません。

$$ -2\sum^n_{i=1} (y_i – (\beta_0 + \beta_1x_{1i} + \beta_2x_{2i}) ) = 0 $$

-2で両片を割って削除して綺麗にする他、β₀はn回足されるだけなので、∑の外に出せます。符号も整理しましょう。以降も同様に進めていきます。
$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_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_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 $$

連立方程式を解く

ここまでに得られた上記の偏微分の式(正規方程式)を、連立方程式として解いていきます。まずは最も簡単なので、β₀の値を求めてみましょう。β₀を左辺におき、全体を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} $$

β₁の正規方程式と、β₂の正規方程式はちょっと面倒です。
まずはxᵢ1とxᵢ2とyᵢを、以下のように変換します。

$$ x_1i = (x_1i – \bar{x_1}) + \bar{x_1} $$

$$ x_2i = (x_2i – \bar{x_2}) + \bar{x_2} $$

また、以下を適用します。n個説明変数の合計 = 説明変数の平均×n、という以上の意味はありません。

まずはβ₁の正規方程式を見ていきます。

$$ 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 \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 \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₁の分散、 $\sum^n_{i=1} (x_{1i} – \bar{x_1})(x_{2i} – \bar{x_2})$ をnで割ったものはx₁とx₂の共分散、 $\sum^n_{i=1} (x_{1i} – \bar{x_1})(y_i – \bar{y})$ をnで割ったものはx₁と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_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}$

β₁とβ₂を求める

かなり複雑な計算をしてきましたが、あともう少しです。
先に求めた二種類の分散と共分散と平均の式に以下のβ₀を代入すると、最終的に分散と共分散だけの式が作られ、これを連立方程式とします。

$$ \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} $$

こちらを簡単に解いて、β₁とβ₂を得るためには、行列に変換することが有効です。

$$ \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} $$

左辺をβ₁とβ₂のみにするために、右辺に逆行列を掛けます。

$$ \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₁とx₂の間に、強い正の相関があることを示します。x₁ᵢが大きい時にx₂ᵢも大きくなる、あるいはx₁ᵢが小さい時にx₂ᵢも小さくなるならば、 $(x_{1i} – \bar{x}_1)(x_{2i} – \bar{x}_2)$ の値は正の値を取りやすくなるので、結果として共分散の値は正の値で0から遠ざかります。

負の値で大きい場合

x₁とx₂の間に、強い負の相関があることを示します。x₁ᵢが大きい時に、x₂ᵢが小さくなる、あるいはx₁ᵢが小さい時に、x₂ᵢが大きくなるならば、 $(x_{1i} – \bar{x}_1)(x_{2i} – \bar{x}_2)$ の値は負の値を取りやすくなるので、結果として共分散の値は負の値で0から遠ざかります。

0に近い場合

x₁とx₂の間の相関が弱い事を表します。x₁ᵢが大きいとしてもx₂ᵢの値が大きくなるとは限らない、あるいはx₁ᵢが小さいとしてもx₂ᵢの値が小さくなるとは限らないので、 $(x_{1i} – \bar{x}_1)(x_{2i} – \bar{x}_2)$ の値は正にも負にもランダムになりやすく、結果としてその和は0に近づくと考えられます。

複数の説明変数間に強い相関があるために、それぞれの説明変数が従属変数に与える影響が評価しにくくなることを、多重共線性と言います。

多重共線性が及ぼす影響の数学的解説

β₁を導出する数式を例に解説します。

$$ \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} $$

Cov(x₁,x₂)が0から遠ざかる場合、Cov(x₁,x₂)²の値が大きくなるため、β₁の分母が小さくなります。よって、少しの値の変動でβ₁の値が変わってしまうことになり、β₁は非常に不安定になります。
x₁とx₂のいずれかが説明変数として正しくないと判断できるため、どちらかをモデルに含めない(次元削減)など、説明変数の選択を再度実施する必要があります。VIFの話は一旦ここでは省略しますが、VIFもといR²は分散と共分散から得られる値です。

その他の影響

多重共線性の話からは少し外れますが、Cov(x₁,y)が何を意味するかを簡単に解説します。
先述のβ₁の値の分子に着目をすると、Var(x₂)Cov(x₁,y)が大きく、Cov(x₂,y)Cov(x₁,x₂)²の値が小さいほど、β₁の値は大きくなります。これらの値は、以下を意味します。

Var(x₂): β₁に、x₁以外の説明変数の影響を考慮している。
Cov(x₁,y): x₁とyの相関の大きさを表す。
Cov(x₂,y): x₂とyの相関の大きさを表す。
Cov(x₁,x₂)²: x₁とx₂の相関の大きさを表す。

つまり、Cov(x₁,y)が小さければ、説明変数x₁は従属変数yの説明変数としてウェイトが低いと言えるため、β₁の値は小さくなります。
また、Cov(x₂,y)の値が大きい場合にはβ₁が小さくなるため、説明変数x₁は相対的に、説明変数x₂よりも従属変数yの説明変数としてウェイトが低くなります。

終わりに

ここまでで、重回帰モデルという強力なツールについて解説をしてきました。基本的にこんな複雑な計算は、コンピュータに任せるケースが大半だと思います。ただ、中身を知っているとモデル構築時のトラブルシューティングにも役立つのではないかと思うので、どこかで役立てていただけると、書いた奴が喜びます。

ところで、疲れた時に海に行きたくなったり、暖かくて狭くて暗い場所で小さくなったりしたくなるのは、胎内回帰願望というのだそうです。個人的には魂のルフランを想起します。有から無に帰るような極端なことはせずとも、たまには全てを忘れてリラックスしたいものです。
生と死・有と無のような二値に対する推定は、ロジスティック回帰と呼ばれるモデルで行うことができます。次回PART3では、このモデルについて解説したいと思います。

以降はおまけです。

最小二乗法による導出(一般式)

先に簡易版の解説をした偏回帰係数の導出を一般化しています。

以下の残差平方和を最小にするβ₀,β₁,β₂,…,βₚを求めます。

$RSS = \sum^n_{i=1} \{y_i – (\hat{\beta_0} + \hat{\beta_1}x_{1i} + … + \hat{\beta_p}x_{pi}) \}^2$

最初のステップとして、上記のβ₀,β₁,β₂,…,βₚのそれぞれで偏微分します。
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$