#初めに
初投稿です。お手柔らかに願います。

皆様、今日も元気にスーパーデータサイエンティスト目指していますでしょうか。

確認ですが、本記事におけるスーパーデータサイエンティストとは__「手法を理解した上で、適切な場面で適切な手法を用いて分析を行い、有益な示唆を得られる一流のデータサイエンティスト」__を指します。

決して「モデルにデータを突っ込むだけの三流データサイエンティスト」や「数学？なにそれ、おいしいの？なデータサイエンティスト(笑)」ではないので注意してください。

さて、初投稿で何について書こうか悩みましたが、「多変量解析の入り口といえば回帰分析でしょ」って偉い人が言っていたので、回帰分析について書くことにします。

#アプローチ

先に断っておきますが、本記事で扱う回帰分析は重回帰分析です。ご認識おき願います。
さて、まずは回帰分析とは何ぞやって話ですが、当然皆様ご存知だと思います。
まさか回帰分析が何かも知らないのに普段仕事で回帰分析していたりしないですよね？？
今ギクッとした方は「重回帰分析 wikipedia」でググっといてください。

ということで、早速数式を見ていきましょう。
目的変数$Y$と説明変数$X_1, X_2, \ldots, X_n$が存在する時、$Y$の推定値$\hat{Y}$を下式のように求める方法を回帰分析といいます。親の顔より見た式ですね。

\hat{Y} = b_0 + b_1X_1 + b_2X_2 + \ldots + b_nX_n

上式の$b_0, b_1, \ldots, b_n$は未知のパラメータで、これらを__いい感じ__に推定することが回帰分析の目的です。

では__いい感じ__とは何でしょうか？？

今、$N$個の観測値(目的変数と説明変数のセット)があり、それらを$(Y_1, X_{11}, X_{12}, \ldots, X_{1n}),~(Y_2, X_{21}, X_{22}, \ldots, X_{2n}),~\dots,~(Y_N, X_{N1}, X_{N2}, \ldots, X_{Nn})$としましょう。

下記の観測された実測値$Y$と推定値$\hat{Y}$が__近くなる__ならば__いい感じ__の推定といえるのではないでしょうか。


Y=Y_1, Y_2, \ldots, Y_N, ~~~~Y=\hat{Y}_1, \hat{Y}_2, \ldots, \hat{Y}_N

近くなる=__差が小さくなる__と考えれば、次式$Q$を最小化するような$b_0, b_1, \ldots, b_n$を求めることが目的となります。
これがかの有名な最小二乗法というやつです。

Q = E(Y-\hat{Y})^2 = \frac{1}{N}\sum_{i=1}^N(Y_i-\hat{Y}_i)^2

ここで、$E(・)$は期待値を表し、$e = Y-\hat{Y}$を__残差__と呼びます。また、2乗しているのは__差__を正数にするためです。

#回帰式の導出
さて、では$Q$を最小化するような$b_0, b_1, \ldots, b_n$を求めていきましょう。
$Q = \frac{1}{N}\sum_{i=1}^N(Y_i-\hat{Y}_i)^2$に一番最初の式を代入すれば、

Q = \frac{1}{N}\sum_{i=1}^N(Y_i-b_0-\sum_{j=1}^nb_jX_{ij})^2

となります。ここから$Q$の最小化を目指しますが、最小化といえば、そうです、微分です。

何故微分をするか疑問に思った皆さん。
気にしたら負けです。制約なし連続最適化において初手微分は基本です。

早速$\frac{\partial Q}{\partial b_0}=\frac{\partial Q}{\partial b_1}=\ldots=\frac{\partial Q}{\partial b_n}=0$となるような$b_0, b_1, \ldots, b_n$を求めていきます。

まず、$\frac{\partial Q}{\partial b_0}=0$を求めていきましょう。上式を$b_0$で偏微分すれば、

\begin{align}
\frac{\partial Q}{\partial b_0} &= \frac{\partial}{\partial b_0}\frac{1}{N}\sum_{i=1}^N(Y_i-b_0-\sum_{j=1}^nb_jX_{ij})^2\\
&= -2\times\frac{1}{N}\times\sum_{i=1}^N(Y_i-b_0-\sum_{j=1}^nb_jX_{ij})\\
&= -2\times(\frac{1}{N}\sum_{i=1}^NY_i-\frac{1}{N}\sum_{i=1}^Nb_0-\frac{1}{N}\sum_{i=1}^N\sum_{j=1}^nb_jX_{ij})\\
&= -2\times(\bar{Y}-b_0-\sum_{j=1}^nb_j\bar{X}_j)\\

\end{align}

となります。上式の$\bar{Y}, \bar{X}_j$はそれぞれ$Y, X_j$の平均値を表します。あ

現在の目的は$\frac{\partial Q}{\partial b_0}=0$となる$b_0$を求めることなので、以下のように変形することで、$b_0$について解きます。

\begin{align}
0&=-2\times(\bar{Y}-b_0-\sum_{j=1}^nb_j\bar{X}_j)\\
b_0 &= \bar{Y}-\sum_{j=1}^nb_j\bar{X}_j

\end{align}

さて、これで$b_0$を$b_1, b_2, \ldots b_n$を用いて表せたので、残りの$b_1, \ldots, b_n$を一気に求めていきましょう。

と、その前に、せっかく$b_0$を求めたので、$Q$の式に代入しておきます。

\begin{align}
Q &= \frac{1}{N}\sum_{i=1}^N(Y_i-(\bar{Y}-\sum_{j=1}^nb_j\bar{X}_j)-\sum_{j=1}^nb_jX_{ij})^2\\
&= \frac{1}{N}\sum_{i=1}^N\{(Y_i-\bar{Y})-\sum_{j=1}^nb_j(X_{ij}-\bar{X}_j)\}^2\\
\end{align}

ここから行列表記を用いますが、__"行列を見ると蕁麻疹が出る！"__という方は、ちょっと先に行列を用いない表記も書いておくので、そちらを参照してください。
⇒面倒くさいのでやめました。行列の勉強してください。

簡単の為、$y_i=Y_i-\bar{Y},~x_{ij}=X_{ij}-\bar{X}_j$と置き、ベクトル(/行列)$y,~X,~b$そして残差ベクトル$e$を以下のようにあらわします。するとQは以下のようにあらわすことができます。

e = \begin{bmatrix}
e_1\\
e_2\\
\dots\\
e_N
\end{bmatrix} = 
\begin{bmatrix}
Y_1-\hat{Y}_1\\
Y_2-\hat{Y}_2\\
\dots\\
Y_N-\hat{Y}_N
\end{bmatrix}

,~

y = \begin{bmatrix}
y_1\\
y_2\\
\dots\\
y_N
\end{bmatrix}

,~

X = \begin{bmatrix}
x_{11}~~x_{12}~~\dots~~x_{1n}\\
x_{21}~~x_{22}~~\dots~~x_{2n}\\
\dots\\
x_{N1}~~x_{N2}~~\dots~~x_{Nn}
\end{bmatrix}

,~

b = \begin{bmatrix}
b_1\\
b_2\\
\dots\\
b_N
\end{bmatrix}

\frac{1}{N}e^Te=\frac{1}{N}(y-Xb)^T(y-Xb)

記号「$^T$」は転置を表します。あとは$b$で偏微分したとき0になればよいので、

\begin{align}
\frac{\partial}{\partial b}(\frac{1}{N}e^Te) &= \frac{\partial}{\partial b}\{\frac{1}{N}(y-Xb)^T(y-Xb)\}\\
&= \frac{\partial}{\partial b}(\frac{1}{N}y^Ty-\frac{2}{N}b^TX^Ty+\frac{1}{N}b^TX^TXb)\\
&= -\frac{2}{N}X^Ty+\frac{2}{N}X^TXb
\end{align}

より、

\begin{align}
0 &= -\frac{2}{N}X^Ty+\frac{2}{N}X^TXb\\
X^TXb &= X^Ty\\
b &= (X^TX)^{-1}X^Ty

\end{align}

これで$b=b_1, b_2,\ldots, b_n$を推定することができました。あとは$b_0$の式に今求めた$b_1, b_2,\ldots, b_n$を代入してあげることで全ての未知パラメータを求めることができます。

#正規方程式と分散・共分散
先程出てきた$X^TXb = X^Ty$ですが、この方程式を__正規方程式__といいます。さてこの正規方程式ですが、分散・共分散と密接な関係があります。
まずは両辺をデータ数$N$で割って以下のように表記します。
$$\frac{1}{N}X^TXb = \frac{1}{N}X^Ty$$
ここで左辺の$\frac{1}{N}X^TX$は各データ項目間の分散と共分散を表した行列となります。

\begin{align}
\frac{1}{N}X^TX &=
\frac{1}{N}
\begin{bmatrix}
x_{11}~~x_{21}~~\dots~~x_{N1}\\
x_{12}~~x_{22}~~\dots~~x_{N2}\\
\dots\\
x_{1n}~~x_{2n}~~\dots~~x_{nN}
\end{bmatrix}

\begin{bmatrix}
x_{11}~~x_{12}~~\dots~~x_{1n}\\
x_{21}~~x_{22}~~\dots~~x_{2n}\\
\dots\\
x_{N1}~~x_{N2}~~\dots~~x_{Nn}
\end{bmatrix}\\

&= \begin{bmatrix}
\frac{1}{N}\sum_{i=1}^{N}x_{i1}^2~~\frac{1}{N}\sum_{i=1}^{N}x_{i1}x_{i2}~~\dots~~\frac{1}{N}\sum_{i=1}^{N}x_{i1}x_{in}\\
\frac{1}{N}\sum_{i=1}^{N}x_{i2}x_{i1}~~\frac{1}{N}\sum_{i=1}^{N}x_{i2}^2~~\dots~~\frac{1}{N}\sum_{i=1}^{N}x_{i2}x_{in}\\
\dots\\
\frac{1}{N}\sum_{i=1}^{N}x_{in}x_{i1}~~\frac{1}{N}\sum_{i=1}^{N}x_{in}x_{i2}~~\dots~~\frac{1}{N}\sum_{i=1}^{N}x_{in}^2\\
\end{bmatrix}

\end{align}

$x_{ij}=X_{ij}-\bar{X}_j$であることを思い出せば、行列の各要素は

\begin{align}
\frac{1}{N}\sum_{i=1}^{N}x_{ij}^2 &= \frac{1}{N}\sum_{i=1}^N(X_{ij}-\bar{X}_j)^2~~(対角成分)\\
\frac{1}{N}\sum_{i=1}^{N}x_{ij}x_{ik} &=\frac{1}{N}\sum_{i=1}^N(X_{ij}-\bar{X}_j)(X_{ik}-\bar{X}_k)~~(対角成分以外)

\end{align}

となり、親の顔より見た分散・共分散の式が現れます。
このように対角成分が同項目の分散、それ以外が異項目間の共分散を表す行列を__分散共分散行列__と呼びます。分散共分散行列は$jk$成分=$kj$成分となる対称行列です。

同様に$\frac{1}{N}X^Ty$についても考えてみます。

\begin{align}


\frac{1}{N}X^Ty &=
\frac{1}{N}
\begin{bmatrix}
x_{11}~~x_{21}~~\dots~~x_{N1}\\
x_{12}~~x_{22}~~\dots~~x_{N2}\\
\dots\\
x_{1n}~~x_{2n}~~\dots~~x_{nN}
\end{bmatrix}

\begin{bmatrix}
y_1\\
y_2\\
\dots\\
y_N
\end{bmatrix}\\

&= 

\begin{bmatrix}
\frac{1}{N}\sum_{i=1}^Nx_{i1}y_i\\
\frac{1}{N}\sum_{i=1}^Nx_{i2}y_i\\
\dots\\
\frac{1}{N}\sum_{i=1}^Nx_{i2}y_i\\
\end{bmatrix}

\end{align}

$y_i=Y_i-\bar{Y},~x_{ij}=X_{ij}-\bar{X}_j$なので、ベクトルの各要素は


\frac{1}{N}\sum_{i=1}^Nx_{ij}y_i = \frac{1}{N}\sum_{i=1}^N(X_{ij}-\bar{X}_j)(Y_i-\bar{Y})

となり、項目$X_j$と$y$の共分散になることがわかります。

最後に上で求めた$X$の分散共分散行列を$\Sigma$, $X$と$y$の共分散ベクトルを$\sigma_y$と置くことで、正規方程式の解をさらにスマートに書きましょう。

\begin{align}
\frac{1}{N}X^TXb &= \frac{1}{N}X^Ty\\
\Sigma b &= \sigma_y\\
b &= \Sigma^{-1}\sigma_y
\end{align}

回帰分析における未知パラメータ$b=b_1, b_2, \ldots, b_n$が各データ項目の__分散共分散行列の逆行列__から導出できるということが重要です。後述の多重共線性はこの分散共分散行列が原因で発生します。

#偏回帰係数と決定係数

ここまでで回帰式迄導出できるようになりましたが、スーパーデータサイエンティストになる為にはこれだけでは不十分です。
というよりここまではエクセルでも勝手にやってくれます。真のスーパーデータサイエンティストたるもの、導出した回帰式を適切に読みとり、有益な示唆を得られなければなりません。

そこで、回帰式の分析で最重要項目となる__"変回帰係数"と"決定係数"__についてみていきましょう。

###偏回帰係数
上記の方法で回帰分析を行い、以下の式が得られました。
$$\hat{Y} = b_0 + b_1X_1 + b_2X_2 + \ldots + b_nX_n$$
この時、推定した$b_1, b_2, \ldots, b_n$を__偏回帰係数__といいます。
突然ですがここで質問です。

Q.この偏回帰係数とは何を表しているのでしょうか？

「偏回帰係数は説明変数と目的変数の相関を表しているに決まっているじゃんバーカ」

と思っている方。ではもう一つ質問です。

Q2.それは単に説明変数と目的変数の相関係数を調べるのとどのように違うのですか？

回答できた方、以降は読まなくて大丈夫です。

さて回答できなかった方。結論を言ってしまうと、

偏回帰係数：対応する説明変数がその他の説明変数に与える影響を度外視した、説明変数-目的変数間の相関
相関係数　：説明変数がその他の説明変数に与える諸々の影響を考慮した、説明変数-目的変数間の相関

を表しています。つまり

$$相関係数 = 偏回帰係数^※＋説明変数がその他の変数に与えた影響によって目的変数に与えた影響(間接効果)\ ※正確には標準偏回帰係数です。$$

という関係が成り立ちます。

この相関係数と(標準)偏回帰係数を調べることで、いろんな示唆が得られたり得られなかったりするのですが、もう面倒くさくなってきたのでやる気が出れば後日別記事として書きたいと思います。気になる方は"疑似相関"とかでググってみてください。

###決定係数
エクセルで回帰分析を行った結果を見てみると"R^2"という項目があります。これは__$R^2$値(R-Squere値)__と呼ばれ、日本語では__決定係数__といいます。決定係数は$0$以上$1$以下の値をとり、回帰式全体がどれだけの説明力を持っているかを表しています。端的に言えば、回帰式がデータに対してどれだけフィットしているかです。

決定係数は__全変動の平方和__と__回帰変動の平方和__の比によって求まります。
ここで、全変動、回帰変動、残差変動について説明しておきます。全変動、回帰変動、残渣変動は式で書くと以下のようにあらわせます。

\begin{align}
(全変動) &= Y_i-\bar{Y} \Rightarrow 観測値と観測値の平均値の差\\
(回帰変動) &= \hat{Y}_i-\bar{Y} \Rightarrow 推定値と観測値の平均値の差\\
(残差変動) &= Y_i-\hat{Y}_i \Rightarrow 観測値と推定値の差
\end{align}

先述の通り、決定係数$R^2$は全変動の平方和と回帰変動の平方和の比によって以下のように導出されます。

R^2 = \frac{\sum_{i=1}^N(\hat{Y}_i-\bar{Y})^2}{\sum_{i=1}^N(Y_i-\bar{Y})^2}

また、$回帰変動 = 全変動 - 残差変動$であることを考えれば、

R^2 = \frac{\sum_{i=1}^N(\hat{Y}_i-\bar{Y})^2}{\sum_{i=1}^N(Y_i-\bar{Y})^2} = 1-\frac{\sum_{i=1}^N(Y_i-\hat{Y})^2}{\sum_{i=1}^N(Y_i-\bar{Y})^2}

と表すことをもできます。というよりこっちのほうが一般的によく使われています。

さてこの決定係数ですが、一つ重大な欠点があります。
それは「説明変数の数を増やすと、それがどんなに関係がない変数だとしても決定係数が上がってしまう」ということです。
つまり無限個の説明変数を用意すれば、理論的に決定係数1の回帰式が出来上がります。これが決定係数の罠です。
この問題を解決しているのが「自由度調整済み決定係数」です。自由度調整済み決定係数$R_f^2$は説明変数の数$N_X$を用いて以下のように計算されます。

R_f^2 = 1-\frac{\frac{\sum_{i=1}^N(Y_i-\hat{Y})^2}{N-N_X-1}}{\frac{\sum_{i=1}^N(Y_i-\bar{Y})^2}{N-1}}

説明変数の数が多い場合はこちらを用いるとよいでしょう。

#多重共線性
最後に多重共線性(Multi-Colinearity)、通称マルチコについて説明していきます。
多重共線性とは説明変数間の相関が大きすぎるときに偏回帰係数や決定係数が正しく算出されない問題を言います。

回帰式の導出をした際に偏回帰係数$b$が分散共分散行列の逆行列$\Sigma^{-1}$を用いることで求められることを説明しました。
この分散共分散行列$\Sigma$ですが、説明変数間の相関が大きいと逆行列が存在しません。
結果、偏回帰係数(の絶対値)がぶっ飛んだ値となってしまいます。

さて、ここから多重共線性が起こる原因を数式で追っていきますが、その前に一つ正しておくべきことがあります。
多重共線性の発生条件は__2変数間の相関が高い__では不十分です。正確には説明変数が一次従属性を持っている時に多重共線性は発生します。説明変数が一次従属であるとは以下の式が$a_1 = a_2=\ldots = a_n = 0$以外で成立することを言います。

$$\sum_{k=1}^na_kX_k=0$$

偏別の表現をすると、ある説明変数が他の説明変数の定数倍の足し算で表せてしまう状況です。

それでは多重共線性の発生原因を数式で確かめていきましょう。
ここでの目的は、説明変数が一次従属性を持っているときに分散共分散行列$\Sigma$が逆行列$\Sigma^{-1}$を持たないことを証明することです。

前提として以下の知識が必要になります。

１．行列$A$が逆行列$A^{-1}$を持たない条件は、行列$A$の行列式$|A|$が0となることである。
２．行列$A$の固有値$\lambda$に0が存在する場合、行列式$|A|=0$となる。
３．行列の固有値がすべて0より大きい場合、その行列を正定値行列と呼ぶ。
４．行列$A$が正定値行列の時、0ではない任意のベクトル$v$に対して$v^TAv>0$となる。
５．行列の固有値がすべて0以上の時、その行列を半正定値行列と呼ぶ
６．分散共分散行列$\Sigma$は半正定値行列である。

以上の前提を用いて、以下の流れで分散共分散行列が逆行列を持たないことを証明していきます。

分散共分散行列が半正定値行列であり、かつ正定値行列ではない
$\Rightarrow$分散共分散行列は固有値に0を持つ
$\Rightarrow$分散共分散行列の行列式は0
$\Rightarrow$分散共分散行列は逆行列を持たない

以上より説明変数間に一次従属性があるとき、分散共分散行列$\Sigma$が正定値行列とならないことだけ示せばよさそうです。

さて、今説明変数間に一次従属性があるため、$a_1 = a_2 =\ldots = a_n = 0$以外で以下式を満たすような$a_1, a_2, \ldots, a_n$が存在します。
$$\sum_{j=1}^na_jX_j = 0$$

上式について、両辺期待値をとったものを用意し引き算すると

\begin{align}
\sum_{j=1}^na_jX_j - \sum_{j=1}^na_jE(X_j) &= 0\\
\sum_{j=1}^na_j(X_j-E(X_j)) &= 0\\
a^T\Sigma a &= 0 ~~(\because 両辺を2乗)
\end{align}

以上より、分散共分散行列$\Sigma$は0ベクトルではないベクトル$a$を用いて$a^T\Sigma a=0$であることが分かったため、正定値行列ではないことが示されました。

さて、今示したのは説明変数間に一次従属性があるときに分散共分散行列の逆行列が計算できず回帰式が導出できないことであり、いわば「2変数の相関係数が$1$または$-1$の時」の極端な話をしています。
先述の通り多重共線性は説明変数間の相関係数の絶対値が"大きい"時(正確には一次従属に近い時)に"偏回帰係数が非常に大きくなってしまう"現象なので、相関係数$0.9$の時にも発生します。

ということで、説明変数間で一次従属性がほぼ成り立つ時についても見ていきたいのですが、、
すみません。わかりませんでした。
(一次従属に近い関係の時に分散共分散行列の行列式が小さくなることを示せばよいと思うのですが、、)
わかる方いたら教えてください。

因みに多重共線性の見つけ方についても書こうと思ったんですが、例によって面倒くさいので、やる気が出れば疑似相関とかと一緒に記事書こうと思います。

#終わりに
さて、初めて記事を書いてみた感想ですが、、超面倒くさい!!
正直もう書きたくないです。間違っているところあったら教えてください。(やる気が出れば)直します。

あと、ベクトルを__ボールド__で書きたいんですが、数式内で太文字を書く方法を教えてくれる優しい方いたら教えてください。