📖 もっと詳しく

重回帰分析（multiple regression）は、 複数の説明変数で目的変数を同時に予測する回帰モデル。 「死亡率に効く要因は何か？高齢化率、 病院数、 医療費、 求人倍率のうち、 真に効くのはどれ？」のような問いを扱えます。 本サイトの論文の多くが重回帰を中核に使います。

モデル式：$y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_k x_k + \varepsilon$

核心アイデア：各係数 $\beta_j$ は「他の変数を全て固定した上で、 $x_j$ を1単位増やすと y がどう変わるか」を表します。 つまり、 単純な相関では見えなかった「純粋な効果」を取り出せる。

例：交絡の制御

単変量相関では「有効求人倍率と死亡率に r = +0.308 の正相関」と出ます。 でもこれは疑似相関で、 真の犯人は「高齢化率」。 重回帰で高齢化率も同時に投入すると、 求人倍率の偏回帰係数は β ≈ +0.03（p > 0.4 で非有意）と急減します。 「高齢化を制御後は、 求人倍率は死亡率に独立した影響を持たない」ことが判明。 これが重回帰の威力です。

変数選択の指針：(i) 経済理論・先行研究に基づいて選ぶ、 (ii) AIC、 BIC、 LASSO で自動選択、 (iii) n ≥ 10×(変数の数) の経験則。 n=47 なら 4-5変数までが安全。

結果の解釈ステップ：

1. F検定：モデル全体が意味あるか（Prob < 0.05 か）
2. R²：説明力はどれくらいか
3. 各係数：符号は予想通りか、 大きさは経済的に意味あるか、 p値は？
4. VIF：多重共線性は大丈夫か（< 10 なら可、 < 5 が望ましい）
5. 残差プロット：仮定違反はないか

標準化偏回帰係数 β を使えば、 単位の違う変数同士の影響度を直接比較できます（「1標準偏差動かしたら、 y は何標準偏差動くか」）。 論文では生の係数と標準化係数の両方を併記するのが理想。

👁️ 直感 — 重回帰は「多変数で予測」

重回帰は単回帰の一般化。 複数の説明変数 x₁, ..., xₚ を使って目的変数 y を予測：

$$ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \cdots + \beta_p x_p + \varepsilon $$

各説明変数の影響を独立に定量化できます。 「他の変数を一定にしたときの効果」を測れるのが重要なポイント。

📐 偏回帰係数 — 重回帰の中核概念

重回帰での β_j は「他のすべての変数を一定にしたときに、 x_j が1単位増えると y がどれだけ変わるか」を表します。 単回帰の係数とは違うことに注意。

例：教育費 = β₀ + β₁ 食料費 + β₂ 住居費 + ε

• β₁ は「住居費を一定にしたときの食料費1単位増の効果」
• 単回帰で食料費だけで予測した時の係数とは違う値になることが多い
• 「住居費」の影響を取り除いた、 純粋な食料費の効果

📐 行列形式と OLS 解

$$ \mathbf{y} = X \boldsymbol{\beta} + \boldsymbol{\varepsilon} $$

X は n×(p+1) の設計行列（最初の列が全1で切片用）。 OLS の解：

$$ \hat{\boldsymbol{\beta}} = (X^T X)^{-1} X^T \mathbf{y} $$

⚠️ 多重共線性 — 重回帰の最大の敵

説明変数同士が強く相関していると、 係数の推定が不安定に。 兆候：

• 係数が異常に大きい/小さい
• 係数の符号が直感に反する
• F検定では有意なのに、 個別の t検定では全て非有意
• VIF（分散拡大因子）が 5〜10 を超える

対策

• 相関の高い変数の一方を除外
• Ridge 回帰（係数を縮小）
• LASSO（変数選択）
• PCA で次元削減

🎯 変数選択 — どの変数を入れるか

• ステップワイズ法：1変数ずつ追加/除外（前進・後退・両方向）
• 情報量規準：AIC、 BIC、 調整済みR² で比較
• LASSO：自動的に係数を0にする
• ドメイン知識：因果関係を考慮した変数選択

🐍 Python での重回帰

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['年度']==2023].dropna()

## R風の式（formula API）
model = smf.ols('D2101 ~ A4101 + I3101 + Q5101 + D1101', data=df).fit()
print(model.summary())

## VIF（多重共線性チェック）
from statsmodels.stats.outliers_influence import variance_inflation_factor
X = sm.add_constant(df[['A4101', 'I3101', 'Q5101', 'D1101']])
vifs = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
print(f'VIFs: {dict(zip(X.columns, vifs))}')

## ロバスト標準誤差（不均一分散に対応）
robust = smf.ols('D2101 ~ A4101 + I3101 + Q5101 + D1101', data=df).fit(cov_type='HC3')
print(robust.summary())

## 診断テスト
from statsmodels.stats.diagnostic import het_breuschpagan
bp = het_breuschpagan(model.resid, model.model.exog)
print(f'Breusch-Pagan p={bp[1]:.3f}（>0.05 で等分散）')