はじめに

【統計検定2級・準1級】Pythonで回帰分析実習(1) - Qiita の続きです。
今回は、独立変数（説明変数）が2個以上ある重回帰分析をPythonを使って実際に試し、統計的な解釈を考えます。今回もnumpyやpandasなどのライブラリを使って計算部分は数式をもとに自分で書いていきます。

参考書

長谷川勝也『ゼロからはじめてよくわかる多変量解析』技術評論社 (2004) Chapter 2

理論的な部分の参考としては以下も参照しています。
東京大学教養学部統計学教室『統計学入門 (基礎統計学Ⅰ)』東京大学出版会 (1991) 第13章

題材

前回使ったJRの路線データでできればよかったのですが、独立変数を集めてくるのが大変そうだったので、今回は別のデータを使うことにしました。
統計局ホームページ/日本の統計 2019－第2章　人口・世帯

以下の都道府県ごとのデータを組み合わせ、CSV形式に加工して使用します。

2-2　都道府県別人口と人口増減率
2-13　都道府県別昼間人口と自宅外就業・通学者数
2-14　都道府県別転出入者数
2-16　都道府県別出生・死亡数と婚姻・離婚件数

(2020/01/07) 人口と人口性比の列のみ、行が1つずれていましたので修正しました。

都道府県,人口,人口性比,人口増減率,昼夜間人口比率,転入超過数,出生率,死亡率,自然増減率,婚姻率,離婚率
北海道,5320,89.1,-0.455,99.945,-6569,6.4,11.8,-5.4,4.5,1.92
青森,1278,88.6,-0.965,99.849,-6075,6.3,13.8,-7.5,4.0,1.64
:
:
沖縄,1626,88.5,0.582,99.968,-1112,11.3,8.4,3.0,5.7,2.44

人口（単位：千人）は平成29年の推計人口です。
人口増減率は ((平成27年の人口/平成22年の人口)^(1/5) - 1)×100 で計算した、1年当たりの平均増減率とします（元のデータの数値ではありません）。

これを前回と同様、文字コードをUTF-8として population.csv というファイル名で保存することとします。

実習

以下、Python 3.6.8で動作確認しています。

データ入力

import numpy as np
import pandas as pd
df = pd.read_csv("population.csv")
print(df) # 取り込んだデータの内容を確認

多変数の回帰分析 (2-3節)

このデータを使って、独立変数が2個以上の場合の回帰分析、つまり「重回帰分析」を考えてみましょう。
独立変数と従属変数の組み合わせは色々考えられるのですが、まずは手法を確認するという意味で、分かりやすい関係が出そうな

独立変数：出生率¹と死亡率¹
従属変数：人口増減率

で試してみましょう。

出生率を $x_1$、死亡率を $x_2$、人口増減率を $y$ として、$y = \beta_0 + \beta_1 x_1 + \beta_2 x_2$ という回帰モデルを考えます。1変数の場合と同じで、各データ点の残差（従属変数の観測値（実現値）と理論値（予測値）との差）の二乗和を最小とするように偏回帰係数 $\beta_0, \beta_1, \beta_2$ の値を定めます。
各データ点（観測値） $(x_1^{(1)}, x_2^{(1)}, y^{(1)}), ..., (x_1^{(n)}, x_2^{(n)}, y^{(n)})$ について、従属変数の理論値を $\hat{y}^{(i)} = \beta_0 + \beta_1 x_1^{(i)} + \beta_2 x_2^{(i)}$ とすると、残差の二乗和 $L$ は

\begin{align}
L &= \sum_{i=1}^n (y^{(i)} - \hat{y}^{(i)})^2 \\
&= \sum_{i=1}^n (y^{(i)} - \beta_0 - \beta_1 x_1^{(i)} - \beta_2 x_2^{(i)})^2
\end{align}

となります。$\beta_0, \beta_1, \beta_2$ についてそれぞれ偏微分係数を0とおいたときの解を推定値 $\hat\beta_0, \hat\beta_1, \hat\beta_2$ とすると

\begin{align}
\hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y}
\end{align}

が得られます（$^\top$は行列・ベクトルの転置を表すものとします）。ただし

\begin{align}
\hat{\boldsymbol{\beta}} &= (\hat\beta_0, \hat\beta_1, \hat\beta_2)^\top \\
\boldsymbol{X} &= \begin{pmatrix}
1 & x_1^{(1)} & x_2^{(1)} \\
\vdots & \vdots & \vdots \\
1 & x_1^{(n)} & x_2^{(n)} \\
\end{pmatrix}
\\
\boldsymbol{y} &= (y^{(1)}, ..., y^{(n)})^\top
\end{align}

で、$\boldsymbol{X}^\top \boldsymbol{X}$は逆行列を持つものとします。
この結果を用いて、偏回帰係数を計算してみると

n = len(df)
X = np.c_[np.ones(n), df[["出生率", "死亡率"]].values]
y = df["人口増減率"]
b = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
print(b) # Output: [ 1.00978055  0.10175178 -0.18246347]

よって、人口増減率の理論値 $\hat{y}$ の回帰式として $\hat{y} = 1.010 + 0.102 x_1 - 0.182 x_2$ が得られます。

人口1000人当たり出生率が1増加すると、1年当たりの人口増減率は0.102ポイント増加する。
人口1000人当たり死亡率が1増加すると、1年当たりの人口増減率は0.182ポイント減少する。

1変数の場合と同様に決定係数を求めてみると

TSS = ((y - y.mean()) ** 2).sum()        # 全変動
ESS = ((X.dot(b) - y.mean()) ** 2).sum() # 回帰変動
RSS = ((y - X.dot(b)) ** 2).sum()        # 残差変動
R_2 = ESS / TSS                          # 決定係数
print(R_2) # 0.9014307808651625

となり、約0.901と求められました。
この値が大きいのか小さいのか。後で他の独立変数を取った場合も試してみますので、値を比較してみてください。

偏回帰係数の計算式の導出

前述の参考書では詳細がすっ飛ばされているのですが、偏回帰係数 $\beta_0, \beta_1, \beta_2$ の計算式の導出をしてみます。普通に偏微分して連立方程式として解いても求められるはずですが、独立変数1個でも面倒だった式変形がますます面倒になるので、「ベクトルの微分」を使って見通しよく整理しましょう。

「ベクトルで微分・行列で微分」公式まとめ - Qiita

ここでは、一般化して独立変数の数を $p$ とし、偏回帰係数 $\beta_0, \beta_1, ..., \beta_p$ の計算式を求めます。

\begin{align}
\boldsymbol{\beta} &= (\beta_0, \beta_1, ...,  \beta_p)^\top \\
\boldsymbol{x}^{(i)} &= (1, x_1^{(i)}, ..., x_p^{(i)})^\top \\
\boldsymbol{X} &= (\boldsymbol{x}^{(1)}, ..., \boldsymbol{x}^{(n)})^\top = \begin{pmatrix}
1 & x_1^{(1)} & \cdots & x_p^{(1)} \\
\vdots & \vdots & & \vdots \\
1 & x_1^{(n)} & \cdots & x_p^{(n)} \\
\end{pmatrix}
\\
\boldsymbol{y} &= (y^{(1)}, ..., y^{(n)})^\top \\
\epsilon^{(i)}
&= y^{(i)} - \hat{y}^{(i)} = y^{(i)} - (\beta_0 + \beta_1 x_1^{(i)} + ... + \beta_p x_p^{(i)} ) \\
&= y^{(i)} - \boldsymbol{\beta}^\top \boldsymbol{x}^{(i)} \\
\boldsymbol{\epsilon} &= (\epsilon^{(1)}, ..., \epsilon^{(n)})^\top \\ 
&= \boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta}\\
\end{align}

とします。式を簡単にするため、$\boldsymbol{x}^{(i)}$ に定数項に相当する成分を1次元加えていることに注意してください。
また零ベクトル（すべての成分が0のベクトル）を $\boldsymbol{0}$ と書くことにします。$\beta_0, \beta_1, ..., \beta_p$ についてそれぞれ偏微分係数を0とおくことは、$\frac{\partial L}{\partial \boldsymbol{\beta}} = \boldsymbol{0}$ とおくことと同じです。

上記ページの公式を適宜用いて$\boldsymbol{\beta}$を求めます。以下の内容の式変形を少し丁寧に書いたものになります。
統計学入門−第7章 | 我楽多頓陳館＞雑学の部屋＞雑学コーナー

\begin{align}
L &= \sum_{i=1}^n (\epsilon^{(i)})^2  = \boldsymbol{\epsilon}^\top \boldsymbol{\epsilon} \\
&= (\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta})^\top (\boldsymbol{y} - \boldsymbol{X} \boldsymbol{\beta}) \\
&= \boldsymbol{y}^\top \boldsymbol{y} - \boldsymbol{y}^\top (\boldsymbol{X} \boldsymbol{\beta}) - (\boldsymbol{X} \boldsymbol{\beta})^\top \boldsymbol{y} + \boldsymbol{\beta}^\top \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta} \\
&= \boldsymbol{y}^\top \boldsymbol{y} - 2 \boldsymbol{y}^\top \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\beta}^\top \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta} \\
\frac{\partial L}{\partial \boldsymbol{\beta}} &= \frac{\partial}{\partial \boldsymbol{\beta}} \boldsymbol{y}^\top \boldsymbol{y} - 2 \frac{\partial}{\partial \boldsymbol{\beta}}(\boldsymbol{y}^\top \boldsymbol{X}) \boldsymbol{\beta} + \frac{\partial}{\partial \boldsymbol{\beta}}\boldsymbol{\beta}^\top (\boldsymbol{X}^\top \boldsymbol{X}) \boldsymbol{\beta} \\
&= -2 (\boldsymbol{y}^\top \boldsymbol{X})^\top + (\boldsymbol{X}^\top \boldsymbol{X} + (\boldsymbol{X}^\top \boldsymbol{X})^\top) \boldsymbol{\beta} \\
&= -2 \boldsymbol{X}^\top \boldsymbol{y} + 2\boldsymbol{X}^\top \boldsymbol{X}\boldsymbol{\beta}
\end{align}

よって、$\boldsymbol{X}^\top \boldsymbol{X}$に逆行列が存在するとき、$\frac{\partial L}{\partial \boldsymbol{\beta}} = \boldsymbol{0}$ の解を推定値 $\hat{\boldsymbol{\beta}}$ とすると

\begin{align}
\hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y}
\end{align}

が得られます。

回帰分析の統計的理論（2-4節）

回帰式を得ることはできましたが、次は「得られた回帰式がどれだけ信頼できるか」を統計的に考えましょう。
ここまで統計学というより線形代数学の話ばかりでしたので。

今までは、従属変数の観測値から最適な回帰式を1つ求めるという考え方をとっていました。
ここで、従属変数をある確率分布に従う確率変数と考えるとどうなるでしょう。つまり、独立変数$\boldsymbol{X}$の値が同じであっても、何度も観測を行えば毎回異なる観測値$\boldsymbol{y}$が得られると考えます。このとき、前述の式で得られる偏回帰係数 $\hat{\boldsymbol{\beta}}$ も観測ごとに変わる確率変数となります。真の偏回帰係数の値を我々は知ることができませんが、偏回帰係数の従う確率分布を考えることはできます。このあたりの考え方は統計検定2級範囲の母平均の区間推定（真の母平均を知ることはできず、母平均の確率分布を求める）などと同じです。

例えば、実際に回帰式を求めたときに、回帰式で重要な役割を果たしている独立変数と、ほとんど影響しない独立変数があるかもしれません。仮に $\beta_j = 0$ であるならば、対応する独立変数 $x_j$ は従属変数の値に何の影響もありません。そこで、観測値から回帰式を得たときに、帰無仮説 $\beta_j = 0$、対立仮説 $\beta_j \neq 0$ を立て、帰無仮説が棄却されれば、その独立変数はその回帰式で確かに従属変数を説明している、と考えることができます。そのためには、観測値から得られた $\hat{\beta_j}$ の確率分布を考える必要が出てきます。

ここでは、各データの誤差 $\epsilon^{(i)}$ が互いに独立に平均0, 分散$\sigma^2$の正規分布に従うと仮定します。すると

（確率変数としての）$y^{(i)}$ は互いに独立に平均$\hat{y}^{(i)}$, 分散$\sigma^2$の正規分布に従う
$\hat{\boldsymbol{\beta}} = (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y}$ より、$\hat\beta_j$ は　$y^{(1)}, ..., y^{(n)}$ の線形和になり、やはり正規分布に従う

ということがいえます。
この $\hat{\boldsymbol{\beta}}$ の平均は、真の偏回帰係数 $\boldsymbol{\beta}$ に一致します（証明は後述）。
また、$S_{ij}$ を $(\boldsymbol{X}^\top \boldsymbol{X})^{-1}$ の$(i,j)$成分と定める（ただし、$\beta_j$の番号付けに合わせて左上を(0,0)成分とします）と、各 $\hat{\beta_j}$ の分散は $\sigma^2 S_{jj}$ となります（証明は後述）。
よって、$\hat{\beta_j}$ を平均0、分散1に標準化した

Z_j = \frac{\hat{\beta_j} - \beta_j}{\sqrt{\sigma^2 S_{jj}}}

は標準正規分布に従うことになりますが、例によって誤差分散 $\sigma^2$ の値を知ることはできません。実際には観測値から推定するしかありませんが、不偏推定量にするには分母をデータの数から（独立変数+1）を引いた値とします。独立変数の個数を$p$とすると

s^2 = \frac{\sum_{i=1}^n (y^{(i)} -  \hat{y}^{(i)})^2}{N-p-1}

が $\sigma^2$ の不偏推定量となり、これを誤差分散の推定値とします。$s^2$ で $\sigma^2$ を置き換えた

t_j = \frac{\hat{\beta_j} - \beta_j}{\sqrt{s^2 S_{jj}}}

は、自由度 $n-p-1$ の $t$分布に従います。
これが$t$分布に従う理由は、以下の「t統計量 (10.4)」とほぼ同じです。
【統計検定2級】正規母集団に対する仮説検定に出てくる確率分布まとめ - Qiita

よって、$\beta_j = 0$ といえるかどうかを検定するときには、

t_j = \frac{\hat{\beta_j} - 0}{\sqrt{s^2 S_{jj}}} = \frac{\hat{\beta_j}}{\sqrt{(y^{(i)} - \sum_{i=1}^n y^{(i)})^2 S_{jj} / (N-p-1)}}

の値を計算して、求めた値を自由度$n-p-1$の$t$分布表から探して$P$値を求めます。$t$値の絶対値が大きいほど$P$値は小さくなり、その独立変数が意味を持っているということになります。

また、先程求めた誤差分散の推定値 $s^2$ は、回帰式で表現しきれない部分を表すわけですので、これが（相対的に）小さいほうが（大雑把には）良い回帰式といえるでしょう。

前置きが長くなりましたが、いよいよ実践。
先程のデータで、2つの独立変数に対応する偏回帰係数について、帰無仮説 $\beta_j = 0$ に対する$t$統計量と$P$値を計算してみましょう。試験で$P$値を求めるときは$t$分布表を用いますが、ここではscipy.statsモジュールで計算してしまいます。

from scipy import stats
p = X.shape[1] - 1                               # 独立変数の個数
S = np.linalg.inv(X.T.dot(X))                    # (X^T X) の逆行列
s_2 = ((y - X.dot(b)) ** 2).sum() / (n - p - 1)  # 誤差分散の推定値
t = b / np.sqrt(s_2 * S.diagonal())              # t統計量
P = (1 - stats.t.cdf(np.abs(t), n - p - 1)) * 2  # P値
print(s_2) # 0.015396963025934803
print(t)   # [  3.273055     3.94856867 -13.93114362]
print(P)   # [0.00207535 0.00028009 0.        ]

統計検定2級ではRの出力結果を読み取らせる問題が必ず1問は出てきますが、その出力はこういう計算で得られているのですね。

この結果を見ると、有意水準5%で $\beta_1 = 0, \beta_2 = 0$ はいずれも棄却されます。そして、出生率より死亡率のほうがよく効いている、と読み取れます。
一方、$s^2 = 0.0154$ より人口増減率 [%] の標準誤差は $s = 0.124$ と推定されるのですが、この結果がよいのかどうか、比較対象がないのでなんともいえませんね。

偏回帰係数の平均の導出

以下、ベクトルの期待値は各成分の期待値を並べたベクトルと思って読んでください。
各データの誤差 $\epsilon^{(i)}$ の平均を0とすると

\begin{align}
{\rm E}[\hat{\boldsymbol{\beta}}] &= {\rm E}[(\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y}] \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top {\rm E}[\boldsymbol{y}] \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top {\rm E}[\boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{\epsilon}] \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta} \\
&= \boldsymbol{\beta}
\end{align}

このように、偏回帰係数の平均は真の偏回帰係数 $\boldsymbol{\beta}$ と一致します（不偏推定量）。
偏回帰係数を推定する方法は最小二乗法以外にも色々考えられますが、不偏推定量が得られることは最小二乗法の大きなメリットといえます。

偏回帰係数の分散の導出

このページを参考に、表記を合わせて途中計算を補完しました。
２.４節　偏回帰係数とその検定 | データ解析/基礎と応用

以下、${\rm V}[\hat{\boldsymbol{\beta}}]$ と書けば、分散共分散行列を表すこととします。つまり、$(i, j)$成分が、$\hat{\beta_i}, \hat{\beta_j}$の共分散（$(i, j)$成分は$\hat{\beta_i}$の分散。ただし、$\beta_j$の番号付けに合わせて左上を(0,0)成分とします）であるような行列です。

\begin{align}
\hat{\boldsymbol{\beta}} &= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{y} \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top (\boldsymbol{\epsilon} + \boldsymbol{X} \boldsymbol{\beta}) \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\epsilon} + (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{X} \boldsymbol{\beta} \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\epsilon} + \boldsymbol{\beta}
\end{align}

となることを用いて

\begin{align}
{\rm V}[\hat{\boldsymbol{\beta}}] &= {\rm E}[(\hat{\boldsymbol{\beta}} - {\rm E}[\hat{\boldsymbol{\beta}}])(\hat{\boldsymbol{\beta}} - {\rm E}[\hat{\boldsymbol{\beta}}])^\top] \\
&= {\rm E}[(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})(\hat{\boldsymbol{\beta}} - \boldsymbol{\beta})^\top] \\
&= {\rm E}[((\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\epsilon})((\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\epsilon})^\top] \\
&= {\rm E}[(\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\epsilon} \boldsymbol{\epsilon}^\top \boldsymbol{X} ((\boldsymbol{X}^\top \boldsymbol{X})^{-1})^\top] \\
&= {\rm E}[(\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top \boldsymbol{\epsilon} \boldsymbol{\epsilon}^\top \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1}] \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top {\rm E}[\boldsymbol{\epsilon} \boldsymbol{\epsilon}^\top] \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \\
&= (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \boldsymbol{X}^\top (\sigma^2 \boldsymbol{I}) \boldsymbol{X} (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \\
&= \sigma^2 (\boldsymbol{X}^\top \boldsymbol{X})^{-1} \\

\end{align}

と整理できます（$\boldsymbol{I}$は単位行列）。ここで $S_{ij}$ を $(\boldsymbol{X}^\top \boldsymbol{X})^{-1}$ の$(i,j)$成分と定義すれば、各 $\hat{\beta_j}$ の分散は $\sigma^2 S_{jj}$ とわかります。

回帰式の検討 (2-6節)

さて、今まで人口増減率を推定するために、独立変数として出生率と死亡率の2つを用いました。ここでは、この独立変数で本当によいのか、他の独立変数の組み合わせのほうがよいということはないか、を検討します。もしかすると、出生率・死亡率以外にも人口増減率に影響を及ぼす要素があるかもしれませんよね。

仮に人口増減率（と、名義尺度である都道府県名）以外の変数から独立変数を任意に2個選ぶとした場合、どのような組み合わせがよいのでしょうか？

import itertools
y = df["人口増減率"]
vars = [x for x in df.columns if x not in ["都道府県", "人口増減率"]]
for exp_vars in itertools.combinations(vars, 2):     # 2つの独立変数を選ぶ
    X = np.c_[np.ones(n), df[list(exp_vars)].values]
    b = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)    # 偏回帰係数を計算する
    TSS = ((y - y.mean()) ** 2).sum()                # 全変動
    ESS = ((X.dot(b) - y.mean()) ** 2).sum()         # 回帰変動
    R_2 = ESS / TSS                                  # 決定係数
    p = X.shape[1] - 1                               # 独立変数の個数
    S = np.linalg.inv(X.T.dot(X))                    # (X^T X) の逆行列
    s_2 = ((y - X.dot(b)) ** 2).sum() / (n - p - 1)  # 誤差分散の推定値
    t = b / np.sqrt(s_2 * S.diagonal())              # t統計量
    P = (1 - stats.t.cdf(np.abs(t), n - p - 1)) * 2  # P値
    print("独立変数:        ", exp_vars)
    print("偏回帰係数:      ", b)
    print("決定係数:        ", R_2)
    print("t統計量:         ", t)
    print("P値:             ", P)
    print("誤差分散の推定値:", s_2)
    print()

すべての結果を出力すると長すぎるので、決定係数が0.9以上の組み合わせに対する出力を抜粋すると以下のようになります。
(2020/01/07) データ修正に伴い結果を差し替えました。

独立変数:         ('人口', '自然増減率')
偏回帰係数:       [1.35396495e-01 2.92466698e-05 1.38354513e-01]
決定係数:         0.9202813093155153
t統計量:          [ 2.50498414  4.06232118 16.12219163]
P値:              [0.01602233 0.00019698 0.        ]
誤差分散の推定値: 0.012452424232620031

独立変数:         ('転入超過数', '自然増減率')
偏回帰係数:       [2.36387964e-01 6.30652813e-06 1.43339245e-01]
決定係数:         0.9235317925311585
t統計量:          [ 6.24974791  4.36740999 18.56543006]
P値:              [1.44848286e-07 7.53394644e-05 0.00000000e+00]
誤差分散の推定値: 0.011944683881961054

独立変数:         ('出生率', '死亡率')
偏回帰係数:       [ 1.00978055  0.10175178 -0.18246347]
決定係数:         0.9014307808651792
t統計量:          [  3.273055     3.94856867 -13.93114362]
P値:              [0.00207535 0.00028009 0.        ]
誤差分散の推定値: 0.015396963025934803

独立変数:         ('出生率', '自然増減率')
偏回帰係数:       [ 1.02915847 -0.08316346  0.18254551]
決定係数:         0.9030150527461578
t統計量:          [ 3.35453348 -2.39398938 14.07003206]
P値:              [0.00164462 0.02099475 0.        ]
誤差分散の推定値: 0.01514949250939055

独立変数:         ('死亡率', '自然増減率')
偏回帰係数:       [ 1.00487249 -0.08024897  0.10178654]
決定係数:         0.902445538157104
t統計量:          [ 3.30081126 -2.33256484  4.02629601]
P値:              [0.00191782 0.02430887 0.0002203 ]
誤差分散の推定値: 0.01523845329397937

独立変数:         ('自然増減率', '婚姻率')
偏回帰係数:       [-0.56327258  0.12603311  0.16092148]
決定係数:         0.9001397075218234
t統計量:          [-1.3456451   7.29460315  2.0734612 ]
P値:              [1.85310560e-01 4.23937485e-09 4.40154160e-02]
誤差分散の推定値: 0.015598634589388555

この中だと、転入超過数と自然増減率の組み合わせが、決定係数が大きく、誤差分散の推定値も小さいのでよいですかね？しかし、転入超過数は割合ではなく実数なので、元の人口によって自然増減率への影響は変わるはずです。ここは、代わりに人口で割ったものを転入超過率とでも呼んで独立変数とすべきでしょう。

df["転入超過率"] = df["転入超過数"] / (df["人口"] * 1000) * 100
n = len(df)
exp_vars = ("転入超過率", "自然増減率")
y = df["人口増減率"]
X = np.c_[np.ones(n), df[exp_vars].values]
b = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
TSS = ((y - y.mean()) ** 2).sum()                # 全変動
ESS = ((X.dot(b) - y.mean()) ** 2).sum()         # 回帰変動
R_2 = ESS / TSS                                  # 決定係数
p = X.shape[1] - 1                               # 独立変数の個数
S = np.linalg.inv(X.T.dot(X))                    # (X^T X) の逆行列
s_2 = ((y - X.dot(b)) ** 2).sum() / (n - p - 1)  # 誤差分散の推定値
t = b / np.sqrt(s_2 * S.diagonal())              # t統計量
P = (1 - stats.t.cdf(np.abs(t), n - p - 1)) * 2  # P値
print("独立変数:        ", exp_vars)
print("偏回帰係数:      ", b)
print("決定係数:        ", R_2)
print("t統計量:         ", t)
print("P値:             ", P)
print("誤差分散の推定値:", s_2)
print()

独立変数:         ('転入超過率', '自然増減率')
偏回帰係数:       [0.20976966 0.75632976 0.10990394]
決定係数:         0.9580723972439072
t統計量:          [ 7.50186534  8.42827532 14.31096232]
P値:              [2.11598050e-09 9.87150361e-11 0.00000000e+00]
誤差分散の推定値: 0.006549283387531298

転入超過数をそのまま使う場合と比較すると、決定係数も誤差分散の推定値も良くなっています。
転入超過率と自然増減率は、ともに人口を増やす効果があると考えられるので、対応する偏回帰係数が確かに正の値になっていることに注意しましょう。

独立変数を3個使うとどうでしょうか？

for exp_vars in itertools.combinations(vars, 3):     # 3つの独立変数を選ぶ
    X = np.c_[np.ones(n), df[list(exp_vars)].values]
    b = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)    # 偏回帰係数を計算する
    TSS = ((y - y.mean()) ** 2).sum()                # 全変動
    ESS = ((X.dot(b) - y.mean()) ** 2).sum()         # 回帰変動
    R_2 = ESS / TSS                                  # 決定係数
    p = X.shape[1] - 1                               # 独立変数の個数
    S = np.linalg.inv(X.T.dot(X))                    # (X^T X) の逆行列
    s_2 = ((y - X.dot(b)) ** 2).sum() / (n - p - 1)  # 誤差分散の推定値
    t = b / np.sqrt(s_2 * S.diagonal())              # t統計量
    P = (1 - stats.t.cdf(np.abs(t), n - p - 1)) * 2  # P値
    print("独立変数:        ", exp_vars)
    print("偏回帰係数:      ", b)
    print("決定係数:        ", R_2)
    print("t統計量:         ", t)
    print("P値:             ", P)
    print("誤差分散の推定値:", s_2)
    print()

例えばこういう組み合わせが見つかります。（これが一番よいというわけではありません。結果の解釈がしやすい一例です）

独立変数:         ('自然増減率', '婚姻率', '離婚率')
偏回帰係数:       [-0.30054533  0.12984192  0.19858489 -0.25130329]
決定係数:         0.9109243128862279
t統計量:          [-0.72219064  7.82613187  2.61426129 -2.28169067]
P値:              [4.74086687e-01 8.36975600e-10 1.22779282e-02 2.75178580e-02]
誤差分散の推定値: 0.014237611977607154

自然増減率と婚姻率が高いほど人口は増え、離婚率が高いほど人口は減ります。納得できそうな結果ですね。偏回帰係数は（定数項を除き）3つとも5%有意です。
一方で、あまり良くない例も見ていきましょう。

独立変数:         ('出生率', '死亡率', '自然増減率')
偏回帰係数:       [ 1.03883742 -0.2559058   0.17133451  0.35353742]
決定係数:         0.9034996234109349
t統計量:          [ 3.3482056  -0.68534239  0.4646739   0.96013752]
P値:              [0.00169918 0.49680614 0.6445094  0.34235404]
誤差分散の推定値: 0.015424353851135318

出生率が上がると人口が減り、死亡率が上がると人口が増えるという奇妙な結果になっています。
これは、「自然増減率 = 出生率 - 死亡率」という独立変数間の関係があるからでして、多重共線性というやっかいな問題を引き起こします。
多重共線性とは何？ Weblio辞書

変数選択を行わない場合には，独立変数相互間に相関の高いものは含めないほうがよい。
　もし，それらの中に独立でないものが含まれていると（例えば変数 A，B とその合計値 C = A + B が共に含まれていると）分析は失敗する。
　場合によっては，各独立変数と従属変数との相関係数の符号と，偏回帰係数の符号が一致しない場合が生ずる。（中略）このようなことが起きるのは，独立変数間に相関の高いものが混ざっていることが原因である（ある変数で予測しすぎた部分を別の変数で打消しているような場合がある）。しかし，このようなことは因果関係を考える上では不都合なので，符号が一致しない独立変数を除いた重回帰式を探索するとよいであろう。

多重共線性が起きているかどうかは、ある程度は検出しようがあるでしょうが、「出生率が高いほど人口が増えないとおかしい」といった話を機械はなかなか分かりません。「符号が一致しない独立変数を除いた重回帰式」って実は結構難しいような。最終的な結果の解釈は、人間がちゃんと考える必要があるところです。

ということで

長かった回帰分析、とりあえず終わり。
次回はChapter 3「主成分分析」の予定です。

人口1,000人あたりの出生数・死亡数 ↩