🔖 キーワード索引

この用語ページの主要トピックを一覧から飛べます。

🎨 直感で掴む

単回帰は「散布図に 直線 を引く」絵で理解できます。 では説明変数が 2 個になったら？ y, x₁, x₂ の 3 次元空間で 平面 を引くイメージです。 3 個以上だと視覚化はできませんが、 数式上は 超平面 を引いていることになります。

🎬 ストーリー：47 都道府県の死亡率を 4 変数で説明する

SSDSE-B-2026 で死亡率（A4200 系の代理指標）を 4 変数で説明したい：

• x₁＝高齢化率（A1303 ÷ A1101 × 100）：高齢化が進めば死亡率は上がる（+ 期待）
• x₂＝人口10万人あたり病院数（I510120 由来）：医療体制が手厚いほど下がる？（− 期待）
• x₃＝新規求職申込件数（F3101）：都市部ほど件数が多く高齢化が低い → 見かけ上の負相関？（− 期待）
• x₄＝年平均気温（B4101）：寒冷地ほど高い？（− 期待）

単純に「死亡率 vs 求人倍率」の単相関を取ると、 たとえば r = +0.31 と出るかもしれません。 でも本当に「不景気 → 死亡率上昇」なのでしょうか？ 実は東北など雪国では同時に高齢化が進んでいる。 つまり「高齢化率」という第 3 変数が両方を動かしている 交絡（confounding） の典型です。

そこで重回帰の出番。 x₁（高齢化率）を式に同時投入すると、 求人倍率の偏回帰係数 β₃ は「高齢化を一定にしたうえでの求人倍率の追加効果」になり、 多くの場合大きく縮みます。 偏回帰係数こそ、 政策議論で本当に欲しい数値。 これが重回帰の威力です。

🎨 視覚的比喩：「点雲を斜めに切るハサミ」

47 都道府県を 4 次元空間に浮かべた 点雲とイメージしてください。 重回帰は、 この雲を 最も低い残差二乗和で貫く超平面 をハサミで探す作業。 ハサミで切った断面が「y の予測値 ŷ」になります。 各説明変数の β は その軸方向への平面の傾き に対応します。

🎬 もう一つの例：ラーメン屋の売上を 3 変数で予測

理論を覚える前に、 もっと身近な例を 1 つ：

• y = 日次売上（円）
• x₁ = 客数（人）
• x₂ = 平均気温（℃）
• x₃ = 平日／休日ダミー

「単純に客数だけで回帰したら β = 950 円」だったとしても、 重回帰で気温・曜日を入れると β = 1,050 円かもしれない。 なぜなら、 暑い休日には客数だけでなく 客単価も上がる 傾向があり、 単回帰は「客数だけ」の効果ではなく「客数＋それに付随するもの」を混ぜていたから。 重回帰で他要因を一定にして、 ようやく 客 1 人あたりの真の売上貢献が見えてくる。

📐 平面の式を絵にする — p=2 の具体例

$\hat{y} = 2 + 3x_1 - 1.5x_2$ を 3 次元空間で描くと、 切片 2 で、 x₁ 方向に 3 ずつ上がり、 x₂ 方向に 1.5 ずつ下がる「斜めの板」になります。 残差は各点から板までの「縦方向の距離」。 これらを二乗して合計した値が最小になる板を選ぶのが OLS です。 視覚化ツール（plotly や matplotlib の 3D 散布図）で一度自分の手で描いてみると、 偏回帰係数の「軸方向の傾き」という意味が腑に落ちます。

🎬 第三の例：マンション価格モデル（不動産業界の主軸）

「マンション価格を予測する」のは典型的な重回帰タスク。 説明変数として築年数・面積・駅距離・階数・周辺所得・最寄り路線などを使う。 例えば

$\hat{\text{価格}} = 800 + 30 \cdot \text{面積}_{m^2} - 5 \cdot \text{築年数} - 2 \cdot \text{駅距離}_{分} + 10 \cdot \text{階数} + \dots$

これだと「同じ面積・駅距離・階数・所得帯のマンションで比べたとき、 築 1 年増えるごとに価格が 5 万円下がる」という解釈ができる。 査定エンジンの古典的バックボーン。 SUUMO や HOMES の自動査定 API はこの拡張版を使っている。

🎨 R282 補強: 重回帰の視覚的理解

重回帰分析を「式」だけでなく「形」で掴むために、 4 つの観点から視覚化を整理する。 ここでは SSDSE-B-2026（47 都道府県）で実際に確認できる現象を中心に扱う。

1. 単相関 vs 偏相関の対比

単相関は「他の変数を一切無視した素朴な関係」、 偏相関は「他の変数を統計的に固定したうえでの純粋な関係」。 例: 死亡率 y と平均気温 x₄ の単相関は r = −0.44（寒い県ほど死亡率が高い）。 しかし高齢化率 x₁ を固定すると偏相関は ≈ 0 に縮む → 気温の効果は実は高齢化率の影づくり。

2. 残差プロットで仮定診断

残差 eᵢ = yᵢ − ŷᵢ を予測値 ŷᵢ の関数としてプロットしたとき、 ランダムに散らばっていれば線形性・等分散性が概ね成立。 曲線パターンが見えたら非線形項追加、 ラッパ状なら HC3 ロバスト SE か WLS。

3. VIF と多重共線性

VIFⱼ = 1 / (1 − Rⱼ²) は「変数 xⱼ を残りの説明変数で回帰したときの説明されにくさの逆数」。 5 を超えたら警戒、 10 で深刻。 多重共線性が深刻な場合、 個別の β̂ は不安定だが予測そのものは依然有効、 という点を見落とさない。

4. 標準化係数で重要度比較

単位が異なる x（人口、 温度、 倍率）を直接比較できないため、 z 標準化した X で回帰し直すか、 β*ⱼ = β̂ⱼ × (sₓⱼ / s_y) を計算する。 SSDSE 例では x₁: +1.02 が圧倒的、 x₂-x₄ は ≤ 0.07 で実質寄与しない。

これら 4 つの視点を毎回確認すれば、 重回帰の出力を「数字の羅列」から「物語」に変えられる。

5. 🐍 補強コード（R282）

このコードでやること： SSDSE-B-2026 を読み込み、 4 説明変数で重回帰し、 偏相関・残差プロット・VIF・標準化係数を一気に出す。

📥 入力: data/raw/SSDSE-B-2026.csv（cp932、 2023 年抽出）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

import pandas as pd, numpy as np
from statsmodels.regression.linear_model import OLS
from statsmodels.stats.outliers_influence import variance_inflation_factor as vif

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023].copy()
y = df['A4200'] / df['A1101'] * 1000           # 死亡率 (千分率)
X = pd.DataFrame({
    'age':  df['A1303'] / df['A1101'] * 100,    # 高齢化率 %
    'hosp': df['I510120'] / df['A1101'] * 1e5,  # 10万人あたり病院数
    'job':  df['F3101'],                         # 新規求職申込件数
    'temp': df['B4101'],                         # 平均気温
})
X = pd.concat([pd.Series(1, index=X.index, name='const'), X], axis=1)
res = OLS(y, X).fit()
print(res.summary())
print('VIF:', {c: round(vif(X.values, i), 2) for i, c in enumerate(X.columns) if c != 'const'})

📤 実行例:

💬 高齢化率だけが圧倒的に有意。 他 3 変数は単相関では見えても重回帰では消える → 疑似相関。

このコードでやること： 偏相関係数を全変数ペアで計算し、 単相関 r と比較する。

📥 入力: 同じ X, y

1
2
3
4
5
6
7
8

from pingouin import partial_corr

for v in ['age', 'hosp', 'job', 'temp']:
    covars = [c for c in ['age','hosp','job','temp'] if c != v]
    r_simple = X[v].corr(y)
    pc = partial_corr(data=pd.concat([y.rename('y'), X.drop(columns='const')], axis=1),
                       x=v, y='y', covar=covars)
    print(f'{v}: r={r_simple:+.3f}  partial r={pc["r"].iloc[0]:+.3f}')

📤 実行例:

💬 偏相関により「他変数を一定にした純粋な効果」が見える。 気温の効果は消える。

このコードでやること： 残差プロットと QQ プロットで仮定診断。

📥 入力: res（前段の OLS 結果オブジェクト）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

import matplotlib.pyplot as plt
import scipy.stats as stats

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
ax[0].scatter(res.fittedvalues, res.resid)
ax[0].axhline(0, color='red', ls='--')
ax[0].set_xlabel('Fitted'); ax[0].set_ylabel('Residual')
ax[0].set_title('Residual vs Fitted')

stats.probplot(res.resid, dist='norm', plot=ax[1])
ax[1].set_title('Normal QQ Plot')

plt.tight_layout(); plt.savefig('multireg_diagnostics.png', dpi=130)

📤 実行例: 残差は方向性なく散布、 QQ プロットはほぼ直線 → 仮定 OK

💬 残差プロットで曲線や広がりが見えなければ、 線形性・等分散性は許容範囲。

このコードでやること： 標準化係数で重要度を可視化。

📥 入力: 同じ X, y

1
2
3
4
5
6
7

from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LinearRegression

Xs = StandardScaler().fit_transform(X.drop(columns='const'))
ys = (y - y.mean()) / y.std()
beta_std = LinearRegression().fit(Xs, ys).coef_
print(dict(zip(['age','hosp','job','temp'], beta_std.round(3))))

📤 実行例:

💬 標準化 β は age が圧倒的（1.02）、 他は ≤0.07 → 実質「高齢化率 1 変数モデル」と等価。

6. 📐 数式に値を入れて手で計算する: β = (X'X)^-1X'y の展開

重回帰の核心は 正規方程式 β = (X'X)^-1X'y。 これを SSDSE-B-2026 の 47 都道府県データで具体的に値を代入して計算する。 ここでは簡略化のため、 説明変数を 1 つだけ (高齢化率) に絞り、 切片を含む 2 次元線形回帰として手計算する。 47 件の実数値を全て展開すると紙面が膨大になるため、 最初の 5 県のみ 使用した縮小版で計算過程を示す (本来は 47 件で計算)。

Step 1: データ行列 X と目的変数 y の準備

最初の 5 都道府県の (高齢化率 [%], 死亡率 [千分率]) を取り出す:

これを X (5×2 行列、 切片列 + 高齢化率) と y (5×1 ベクトル) に並べる:

Step 2: X'X の計算 (2×2 対称行列)

X' は X の転置 (2×5)。 X' × X は (2×5) × (5×2) = (2×2):

Step 3: (X'X)^-1 の計算 (2×2 行列の逆行列公式)

2×2 行列 [[a,b],[c,d]] の逆行列は (1/det) × [[d,-b],[-c,a]]、 det = ad - bc。

Step 4: X'y の計算 (2×1 ベクトル)

X' × y = (2×5) × (5×1) = (2×1):

Step 5: β = (X'X)^-1 X'y を計算

解釈: 高齢化率が 1% 上がると、 死亡率は 約 0.69 千分率増加。 切片 -7.37 は (age=0 の理論値)、 外挿のため解釈に使わない。 これが手計算で得られた最初の 5 県だけの結果で、 47 県全体ではより精度高い推定 (本文の β₁ ≈ +0.649) になる。

7. 🐍 同じ計算を Python で再現

このコードでやること: 上記 Step 1〜5 を numpy で機械的に実行し、 手計算の結果と完全に一致することを確認する。

📥 入力: 5 都道府県の age と death の数値

 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
26

import numpy as np

# Step 1: データ準備 (5 都道府県の縮小版)
age   = np.array([33.0, 35.2, 35.0, 29.2, 39.1])
death = np.array([14.8, 17.6, 16.9, 12.7, 19.2])

# 切片列を追加した X 行列を作る
X = np.column_stack([np.ones(5), age])
y = death

# Step 2: X'X
XtX = X.T @ X
print("X'X =\n", XtX)

# Step 3: (X'X)^(-1)
XtX_inv = np.linalg.inv(XtX)
print("(X'X)^-1 =\n", XtX_inv)

# Step 4: X'y
Xty = X.T @ y
print("X'y =", Xty)

# Step 5: β = (X'X)^-1 X'y
beta = XtX_inv @ Xty
print("β (切片, 高齢化率) =", beta)
print(f"回帰式: death = {beta[0]:+.2f} + {beta[1]:+.2f} × age")

📤 実行例 (出力):

💬 手計算 (上の Step 5) と完全に一致。 numpy の np.linalg.inv は内部で LU 分解を使うため、 5×5 以下の小行列なら手計算との誤差は無視できる。 ただし n が大きいと数値安定性のため np.linalg.solve(XtX, Xty) や QR 分解 (np.linalg.lstsq(X, y, rcond=None)) を推奨。

8. 47 都道府県 全件を使った場合との対比

上の手計算は 5 県の縮小版。 47 県すべてで同じ計算をすると、 サンプルサイズが増えることで標準誤差が縮み、 推定値が安定する。 詳しい結果は本ページ前半 (5. 補強コード R282) の OLS 回帰結果を参照。 5 県 vs 47 県の β₁ を比較:

差 0.04 は 小標本誤差。 47 件で計算すれば、 標準誤差 0.046 を満たす中央推定 0.649 が得られる (推定値は ±2 SE で 0.56〜0.74 の範囲)。

📐 数式または定義

重回帰モデルの基本形：

行列形式（設計行列 X は n×(p+1)、 先頭列は切片用の 1）：

OLS（最小二乗法）の解析解：

推定量の分散共分散行列：

決定係数と調整済み R²：

全体検定（F検定）：帰無仮説 H₀: β₁ = β₂ = … = βₚ = 0

個別係数の t検定：

VIF（分散拡大因子）：x_j を他の説明変数で回帰したときの $R^2_j$ を用いて

標準化偏回帰係数 β*（変数を z 標準化してから推定）：

単位の違う変数間で「影響度」を比較するには β* を使う。

Ridge / Lasso 推定量（参考、 正則化版）：

Cook's distance（影響度指標）：

🔬 数式を言葉で読み解く

🔬 数式・定義を「言葉」で読み解く

数式 $y_i = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip} + \varepsilon_i$ は暗号ではなく、 1 文字 1 文字に明確な意味があります。 ここでは「重回帰の式を、 統計を知らない人にも説明できるレベルまで」分解します（4 段構え）。

1. 🎯 何を表す式か（一行で）

「i 番目の県の死亡率 yᵢ は、 切片 β₀ に、 高齢化率の寄与 β₁x₁ と、 病院数の寄与 β₂x₂ と…… を足したものに、 説明しきれないズレ ε を加えたもの」。 つまり 「目的変数 = 直線的な足し算 + 誤差」 というレシピを書いた式です。 「直線的」とは「説明変数を 2 倍したら y への寄与も 2 倍」という関係。 これが線形モデルの定義。

2. 📥 各記号の意味（記号 → 日本語）

• yᵢ：i 番目の観測ユニット（北海道 i=1、 青森 i=2、 …、 沖縄 i=47）の 目的変数の実測値。
• β₀：切片（intercept）。 すべての説明変数が 0 のときの y の値。 多くは「絶対値より、 他の係数の解釈の足場」として使う。 説明変数を z 標準化していれば β₀ は y̅（y の平均）と一致する。
• β₁, β₂, …, βₚ：偏回帰係数（partial regression coefficient）。 「他の x を全部固定したまま、 xⱼ を 1 単位増やしたとき y がどれだけ動くか」。 重回帰の主役。 単回帰係数とは違う値になる（交絡があるとき）。
• xᵢⱼ：i 番目のユニットの j 番目の 説明変数の値。 「青森県の高齢化率＝34.1%」など。 連続でもダミー（0/1）でも OK。
• εᵢ：誤差項（error / disturbance）。 モデルで説明しきれない部分。 観測誤差・未投入変数・偶然のゆらぎ、 すべてここに集約される。 平均 0、 分散 σ² の正規分布が仮定される。
• ŷᵢ（ハット記号付き）：実測 yᵢ に対するモデルの 予測値。 残差 eᵢ = yᵢ − ŷᵢ。
• X（大文字、 行列）：設計行列（design matrix）。 n 行 ×(p+1) 列。 先頭列は切片用の 1。

3. 🧠 言葉での全体読み — 物語として

この式は「料理のレシピ」だと思うと分かりやすい。 完成品（y、 たとえば死亡率）は、 ベースの量 β₀ から始めて、 素材ごとの分量 × 単価 β₁x₁、 β₂x₂、 … を足していくと組み上がる。 ただし、 どんなレシピも「実際に作るとちょっとズレる」。 そのズレが ε（誤差）です。

単回帰との 決定的な違いは、 β₁ が「単独で x₁ と y を見たときの傾き」ではないこと。 重回帰の β₁ は 「他の説明変数 x₂, x₃, …, xₚ を全部一定にした上で」の x₁ の効果。 これを partial（部分・偏）と呼ぶ理由がここにあります。 「高齢化率を一定にしたうえでの、 求人倍率の追加効果」— これは単純な相関よりも はるかに政策的に意味のある情報です。

直感的に 「補正後の効果」「条件付きの効果」 と読むと、 統計を知らない人にも伝わります。 「同じ高齢化率の県同士で比べたとき、 求人倍率の高い県の方が死亡率が 0.18 ポイント低い」のように。 これは Frisch-Waugh-Lovell の定理として正式に証明されており、 「重回帰の β₁ は、 (a) x₁ を他の x で回帰した残差と (b) y を他の x で回帰した残差で、 単回帰した時の傾きと一致する」という美しい性質を持つ。

4. 🔍 推定の「気持ち」 — OLS は何を最小化しているか

公式 $\hat{\boldsymbol{\beta}} = (X^\top X)^{-1} X^\top \mathbf{y}$ は天から降ってきたものではありません。 「残差二乗和 RSS = Σ(y − ŷ)² を最小化する β を解析的に解いた結果」です。 RSS を β で微分してゼロと置くと、 正規方程式（normal equations） $X^\top X \hat{\boldsymbol{\beta}} = X^\top \mathbf{y}$ が出てくる。 これを左から $(X^\top X)^{-1}$ で掛けると公式が得られる。

なぜ「二乗和」を最小化するのか？ — 大きく 3 つの理由があります：

1. 正規分布のもとで最尤と一致：ε が独立同分布の正規 N(0, σ²) に従うとき、 OLS 推定量は最尤推定量（MLE）と完全に一致する。
2. BLUE 性質：ガウス・マルコフ仮定（線形・独立・等分散）のもとで、 線形不偏推定量の中で 分散最小。 これが「OLS は最良（best）」と呼ばれる根拠。
3. 計算が解析的に閉じる：絶対値の和を最小化する LAD（最小絶対偏差回帰）と違い、 OLS は 1 ステップで解ける（逆行列計算）。 数値最適化が要らない。

5. 💬 仮定が崩れるとき — 何が起こるか

重回帰の「美しい性質」（BLUE）は、 ガウス・マルコフ仮定 が成り立つ前提です。 仮定とは：

• (A1) 線形性：y と x の関係が線形。
• (A2) 独立性：εᵢ どうしが独立。 時系列・パネルでは要注意。
• (A3) 等分散性：Var(εᵢ) = σ² で一定。 ヘテロスケダスティシティ（heteroscedasticity）の反対。
• (A4) 誤差の正規性：εᵢ ∼ N(0, σ²)。 t検定・F検定で必要。 推定だけなら不要。
• (A5) 多重共線性なし：説明変数同士が完全に線形従属でない。 強く相関すると数値的に不安定。
• (A6) 外生性：E[ε|X] = 0。 内生性（説明変数と誤差の相関）があると β̂ にバイアス。

仮定が破れた時の症状と対処：

つまり「数式を読む」とは、 単に記号を翻訳することではなく、 仮定が破れた時に何が壊れるか／どう直すかまで含めて理解すること。 ここまで読めて初めて、 重回帰を「使える」と言えます。

6. 📚 因果推論との橋渡し — βⱼ は「効果」か？

最後に重要な注意：βⱼ は 「効果」ではなく「条件付き連関」に過ぎません。 観察データから因果を読みたければ、 (i) 処置の無作為割り付け（RCT）、 または (ii) 適切な操作変数（IV）、 (iii) 差の差（DiD）、 (iv) 不連続回帰（RDD） のいずれかが必要。 「重回帰の係数」を「政策効果」と短絡しないのが、 統計リテラシーの第一歩です。

それでも観察研究で重回帰を使う価値はあります： (a) 予測モデルとして 高精度なら十分実用的、 (b) 仮説生成として どの変数が「効きそう」かを絞り込める、 (c) 政策シミュレーションの第一段階として 大まかな見当をつける、 などの用途で必要不可欠です。

🧮 実値で計算してみる — SSDSE-B-2026

SSDSE-B-2026 から 47 都道府県の死亡率系指標を 4 変数で重回帰した実例で、 数値を追います。 ここでは目的変数として 死亡数（A4200）を 総人口（A1101）で割った千分率を「死亡率」とし、 説明変数は高齢化率・病院数・新規求職申込件数・年平均気温の 4 つを取ります。

📥 使うデータ列（コード対応表）

🧮 1. 記述統計（n=47）

🧮 2. 単相関（y との r）

x₁: r = +0.972（強い正）／x₂: r = +0.52／x₃: r = −0.59／x₄: r = −0.44。 「高齢化率 1 変数で死亡率はほぼ説明できる」が見えます。 では他の 3 変数は「ついで」なのか？ 重回帰で検証します。

🧮 3. OLS 推定結果（4 変数）

R² = 0.947、 調整済み R² = 0.943、 F(4, 42) = 189.5、 p < 0.001。 モデル全体は強く有意。

🧮 4. 沖縄県を代入して予測

沖縄県：x₁=23.8, x₂=5.2, x₃=43877, x₄=23.8 →

$\hat{y} = -6.80 + 0.649 \times 23.8 - 0.010 \times 5.2 + 2.9\text{e-}6 \times 43877 + 0.017 \times 23.8$

= -6.80 + 15.45 − 0.05 + 0.13 + 0.40 = 9.13（‰）。 実測値 ≈ 10.29。 残差 = +1.16。 やや過小予測。

🧮 5. 秋田県を代入して予測

秋田県：x₁=39.1, x₂=5.3, x₃=30175, x₄=13.7 →

$\hat{y} = -6.80 + 0.649 \times 39.1 - 0.010 \times 5.3 + 2.9\text{e-}6 \times 30175 + 0.017 \times 13.7$

= -6.80 + 25.38 − 0.05 + 0.09 + 0.23 = 18.85（‰）。 実測値 ≈ 19.17。 残差 = +0.32。 高精度。

🧮 6. F検定の手計算

R² = 0.947, p = 4, n = 47 →

$F = \frac{R^2 / p}{(1-R^2) / (n-p-1)} = \frac{0.947/4}{0.053/42} = \frac{0.2369}{0.00125} \approx 189.5$

自由度 (4, 42) で F の 1% 臨界値は約 3.83、 0.1% 臨界値は約 5.85。 190 >> 5.85 → p < 0.001、 モデル全体は極めて強く有意。

🧮 7. VIF の手計算（β₁ の場合）

x₁（高齢化率）を x₂, x₃, x₄ で回帰 → R²₁ ≈ 0.77。 よって VIF₁ = 1/(1−0.77) ≈ 4.31。 5 を下回るので多重共線性は許容範囲。

🧮 8. 標準化 β の意味

β₁* = +1.02 は「高齢化率を 1 標準偏差（3.3 ポイント）動かすと、 死亡率が 1.02 標準偏差（≈ 2.14 ‰）動く」。 他の 3 変数の標準化 β は |0.01|, |0.07|, |0.02| なので、 高齢化率の 影響度は他を 15 倍以上上回る。 重回帰の標準化 β は、 政策議論で「どの要因が一番効くか」を語る最重要数値です。

🧮 9. 残差の正規性検定（Shapiro-Wilk）

残差の正規性を Shapiro-Wilk 検定で確認すると W = 0.946、 p = 0.030（n=47）。 0.05 をわずかに下回り、 帰無仮説「正規分布である」は 5% 水準で辛うじて棄却される（残差正規性はやや弱い）。 → t検定・F検定は大標本では頑健だが、 厳密には HC3 ロバスト SE の併用が無難。

🧮 10. 等分散性検定（Breusch-Pagan）

BP 統計量 = 1.82、 自由度 4 のカイ二乗で p = 0.77。 等分散性の帰無仮説は棄却されず、 通常の OLS SE をそのまま信頼できる。

🧮 数式に値を入れて手で計算する: 重回帰係数

合成 5 サンプル 2 説明変数で重回帰係数を計算する。

Step 1: データ

Step 2: 正規方程式 β = (XᵀX)⁻¹Xᵀy

🐍 Python で再現

1
2
3
4
5
6
7

import numpy as np
from sklearn.linear_model import LinearRegression
X = np.array([[1,2],[2,1],[3,3],[4,2],[5,4]])
y = np.array([5,4,9,9,14])
m = LinearRegression().fit(X, y)
print(f"切片: {m.intercept_:.2f}")
print(f"係数: {m.coef_}")

📤 実行結果

💬 手計算 (Step 2) [-0.40, 1.356, 1.889] と Python 出力が完全一致。

🐍 Python 実装

① statsmodels — 統計学的に正統な分析

 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
26
27
28
29
30

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.diagnostic import het_breuschpagan

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026']==2023].dropna(subset=['A4200','A1101','A1303','I510120','F3101','B4101']).copy()

# 派生変数
df['death_rate']    = df['A4200']    / df['A1101'] * 1000   # 千分率
df['aging_rate']    = df['A1303']    / df['A1101'] * 100    # %
df['hosp_per_100k'] = df['I510120']  / df['A1101'] * 100000

formula = 'death_rate ~ aging_rate + hosp_per_100k + F3101 + B4101'
model = smf.ols(formula, data=df).fit()
print(model.summary())

# VIF
X = sm.add_constant(df[['aging_rate','hosp_per_100k','F3101','B4101']])
vif = {c: variance_inflation_factor(X.values, i) for i, c in enumerate(X.columns)}
print('VIF:', {k: round(v, 2) for k, v in vif.items()})

# 等分散検定（Breusch-Pagan）
bp_stat, bp_p, _, _ = het_breuschpagan(model.resid, model.model.exog)
print(f'Breusch-Pagan p = {bp_p:.3f}  ({"等分散OK" if bp_p>0.05 else "ヘテロ"})')

# ロバスト標準誤差（HC3）
robust = smf.ols(formula, data=df).fit(cov_type='HC3')
print(robust.summary())

② scikit-learn — 予測重視 + Pipeline + 交差検証

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

from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.model_selection import cross_val_score, KFold
import numpy as np

X = df[['aging_rate','hosp_per_100k','F3101','B4101']].values
y = df['death_rate'].values
cv = KFold(5, shuffle=True, random_state=42)

models = {
    'OLS'   : LinearRegression(),
    'Ridge' : Ridge(alpha=1.0),
    'Lasso' : Lasso(alpha=0.1, max_iter=10000),
    'EN'    : ElasticNet(alpha=0.1, l1_ratio=0.5, max_iter=10000),
}
for name, est in models.items():
    pipe = Pipeline([('sc', StandardScaler()), ('m', est)])
    scores = cross_val_score(pipe, X, y, cv=cv, scoring='r2')
    pipe.fit(X, y)
    coefs = pipe['m'].coef_
    print(f'{name:6s}  CV R² = {scores.mean():.3f} ± {scores.std():.3f}   '
          f'coef = {dict(zip(["aging","hosp","jobs","temp"], coefs.round(3)))}')

③ numpy / scipy — 行列で「手計算」（理解の最終段階）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

import numpy as np
from scipy import stats, linalg

X_design = np.column_stack([np.ones(len(df)), X])           # (n, p+1)
beta = linalg.solve(X_design.T @ X_design, X_design.T @ y)  # β̂
residuals = y - X_design @ beta
n, p = X_design.shape
sigma2 = (residuals @ residuals) / (n - p)
cov_beta = sigma2 * linalg.inv(X_design.T @ X_design)
se = np.sqrt(np.diag(cov_beta))
t_stats = beta / se
p_values = 2 * (1 - stats.t.cdf(np.abs(t_stats), df=n-p))
t_crit = stats.t.ppf(0.975, df=n-p)
ci_low  = beta - t_crit * se
ci_high = beta + t_crit * se

names = ['intercept', 'aging', 'hosp', 'jobs', 'temp']
for nm, b, s, t_, p_, lo, hi in zip(names, beta, se, t_stats, p_values, ci_low, ci_high):
    print(f'{nm:10s}  β={b:8.3f}  SE={s:.3f}  t={t_:6.2f}  p={p_:.4f}  95%CI=[{lo:.2f}, {hi:.2f}]')

④ matplotlib — 残差プロットによる仮定診断

 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
26

import matplotlib.pyplot as plt
from scipy import stats as sst
import numpy as np

fig, ax = plt.subplots(2, 2, figsize=(12, 10))

# 残差 vs 予測値
ax[0,0].scatter(model.fittedvalues, model.resid)
ax[0,0].axhline(0, color='red', ls='--')
ax[0,0].set(xlabel='Fitted', ylabel='Residual', title='Residuals vs Fitted')

# QQ プロット
sst.probplot(model.resid, dist='norm', plot=ax[0,1])

# スケール - ロケーション
ax[1,0].scatter(model.fittedvalues, np.sqrt(np.abs(model.get_influence().resid_studentized_internal)))
ax[1,0].set(xlabel='Fitted', ylabel='√|Std. Residual|', title='Scale-Location')

# Cook's distance
infl = model.get_influence()
cooks = infl.cooks_distance[0]
ax[1,1].stem(range(len(cooks)), cooks)
ax[1,1].axhline(4/len(df), color='red', ls='--', label='4/n')
ax[1,1].set(xlabel='Observation', ylabel="Cook's D", title="Cook's distance")
ax[1,1].legend()
plt.tight_layout(); plt.savefig('mreg_diagnostics.png', dpi=120)

⑤ scikit-learn — Lasso のパス（係数の縮約過程）

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

from sklearn.linear_model import lasso_path
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler

X_std = StandardScaler().fit_transform(X)
alphas, coefs, _ = lasso_path(X_std, y, alphas=[0.001, 0.01, 0.05, 0.1, 0.5, 1.0, 5.0])
print('coefs shape:', coefs.shape)  # (4, len(alphas))

for i, name in enumerate(['aging','hosp','jobs','temp']):
    plt.plot(-np.log10(alphas), coefs[i], label=name)
plt.legend(); plt.xlabel('-log10(α)'); plt.ylabel('β'); plt.title('Lasso path')
plt.savefig('lasso_path.png', dpi=120)

⑥ statsmodels — 部分回帰プロット（CCPR）

1
2
3

import statsmodels.api as sm
fig = sm.graphics.plot_partregress_grid(model)
fig.tight_layout(); fig.savefig('partregress.png', dpi=120)

⑦ statsmodels — モデル比較（AIC / BIC）と階層的 F検定

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

m1 = smf.ols('death_rate ~ aging_rate', data=df).fit()
m2 = smf.ols('death_rate ~ aging_rate + hosp_per_100k', data=df).fit()
m3 = smf.ols('death_rate ~ aging_rate + hosp_per_100k + F3101 + B4101', data=df).fit()

print(f'{"Model":<6}{"AIC":>8}{"BIC":>8}{"adjR2":>8}')
for name, m in [('M1', m1), ('M2', m2), ('M3', m3)]:
    print(f'{name:<6}{m.aic:8.2f}{m.bic:8.2f}{m.rsquared_adj:8.3f}')

from statsmodels.stats.anova import anova_lm
print(anova_lm(m2, m3))

⑧ scipy / sklearn — Huber 損失でロバスト推定

1
2
3
4

from sklearn.linear_model import HuberRegressor
huber = HuberRegressor(epsilon=1.35)
huber.fit(StandardScaler().fit_transform(X), y)
print('Huber β:', dict(zip(['aging','hosp','jobs','temp'], huber.coef_.round(3))))

⑨ patsy で柔軟な式表記

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

# ── この抜粋で使うデータを用意します ──
import patsy
import statsmodels.api as sm
import pandas as pd

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023]          # 2023 年の 47 都道府県
df = df.copy()
df['death_rate']    = df['A4200'] / df['A1101'] * 1000     # 人口千人あたり死亡数
df['aging_rate']    = df['A1303'] / df['A1101'] * 100      # 高齢化率 %
df['hosp_per_100k'] = df['I510120'] / df['A1101'] * 100000 # 人口 10 万人あたり一般病院数
_E = {'北海道': '北', '青森県': '北', '岩手県': '北', '宮城県': '北', '秋田県': '北',
      '山形県': '北', '福島県': '北'}
df['region'] = df['Prefecture'].map(_E).fillna('その他')   # 2 水準の地域区分

y, X = patsy.dmatrices('death_rate ~ aging_rate + I(aging_rate**2) + aging_rate:hosp_per_100k + C(region)', data=df, return_type='dataframe')
model = sm.OLS(y, X).fit()
print(model.summary())

⑩ Bootstrap で頑健な信頼区間

 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
26
27
28
29
30

# ── この抜粋で使うデータを用意します ──
import numpy as np
import pandas as pd
from sklearn.utils import resample

df = pd.read_csv('data/raw/SSDSE-B-2026.csv', encoding='cp932', skiprows=[1])
df = df[df['SSDSE-B-2026'] == 2023]          # 2023 年の 47 都道府県
# const を除く 4 変数（下の name 一覧と数を合わせる）
X = np.column_stack([
    df['A1303'] / df['A1101'] * 100,          # aging: 高齢化率
    df['I510120'] / df['A1101'] * 100000,     # hosp : 人口10万人あたり一般病院数
    df['F3103'].astype(float) / 1000,         # jobs : 月間有効求人数（千件）
    df['B4101'].astype(float)])               # temp : 年平均気温
y = (df['A4200'] / df['A1101'] * 1000).to_numpy(dtype=float)   # 死亡率


n_boot = 1000
betas = []
for _ in range(n_boot):
    Xb, yb = resample(X, y, random_state=None)
    Xd = np.column_stack([np.ones(len(Xb)), Xb])
    beta_b = np.linalg.lstsq(Xd, yb, rcond=None)[0]
    betas.append(beta_b)
betas = np.array(betas)

ci_low  = np.percentile(betas, 2.5, axis=0)
ci_high = np.percentile(betas, 97.5, axis=0)
print('Bootstrap 95%CI:')
for i, name in enumerate(['const','aging','hosp','jobs','temp']):
    print(f'  {name}: [{ci_low[i]:.3f}, {ci_high[i]:.3f}]')

⚠️ 落とし穴 — よくある失敗 5 件

⚠️ 落とし穴 深掘り 10 項目

基本の落とし穴 5 つは押さえたとして、 さらに見落としがちな「重回帰の罠」を 10 個追加します。

1. シンプソンのパラドックス

サブグループでは正の効果なのに、 全体で集計すると負になる現象。 重回帰でも「グループダミーを入れ忘れる」と、 サブグループの符号と全体の符号が逆転することがある。 必ずグループ別の散布図を確認すること。

2. 範囲外への外挿

学習データの x の範囲外に予測を伸ばすと、 線形の仮定が崩れる可能性大。 「築 50 年のマンション」のデータで学習して「築 0 年」を予測すると、 切片次第で異常値が出る。 必ず学習データの分布範囲を確認。

3. ベルクソンバイアス

「両方の条件を満たすサンプル」のみ集めて回帰すると、 元集団には存在しない相関が出る。 病院に来た患者だけで「症状 A と病気 B の関係」を見ると、 母集団とは異なる結論になることがある。

4. 選択バイアス

アンケート回答者だけのデータで回帰すると、 「回答しなかった人」の特性が抜ける。 ヘックマンの 2 段階推定で補正可能だが、 手間がかかる。

5. 測定誤差バイアス（attenuation bias）

説明変数 x に測定誤差があると、 β̂ は 0 方向に縮む（attenuation）。 教育年数の自己申告など、 報告誤差のある変数では β を過小評価しがち。 構造方程式モデル (SEM) や IV で補正可能。

6. 逆因果（reverse causality）

y → x の因果関係が逆方向にも存在すると、 β̂ は混合された値になる。 「病院数→死亡率」の効果と「死亡率→病院数」（高齢化地域に病院誘致）の効果が同時に動く。 時間ラグを取る、 操作変数を使う、 で対処。

7. 欠損データの非ランダム性

欠損が「ランダムに起きる」(MCAR) ではなく、 「ある条件で起きる」(MAR / MNAR) と、 単純な listwise deletion でバイアス。 多重代入や逆確率重み付けで補正。

8. 多重比較問題

20 個の変数を入れて 1 個でも p<.05 になれば、 偶然の可能性が高い（family-wise error rate)。 Bonferroni、 Holm、 FDR で補正する。 重回帰では「F検定が有意か」を最初にチェック。

9. データのリーク

説明変数に「未来の情報」「目的変数の派生」が紛れると、 R² が見かけ上跳ね上がる。 たとえば「翌月の売上を当月の利益で予測」など。 時系列ではデータ分割の時点を厳密に管理。

10. 因果のはしごの誤読

Judea Pearl の「因果のはしご」では、 (1) 連関 → (2) 介入 → (3) 反実仮想、 の 3 段階。 重回帰は (1) のレベル。 介入効果を語るには (2)、 「もし〜だったら」を語るには (3) が必要。 重回帰 1 本で「効果」「介入」「反実仮想」を語ると論理破綻する。

🌐 関連手法・派生

• OLS／単回帰：説明変数 1 個の特殊形。 重回帰の基本構造を理解する出発点。 単回帰 / OLS
• Ridge 回帰：多重共線性のときに $\|\beta\|_2^2$ にペナルティ。 Ridge
• Lasso 回帰：$\|\beta\|_1$ にペナルティで自動変数選択。 Lasso
• ロジスティック回帰：y が 0/1 の二値版。 ロジスティック
• 一般化加法モデル（GAM）：非線形効果を許容。 GAM
• パネルデータ分析：時系列×断面の 2 次元データ。 パネル / 固定効果
• ステップワイズ：変数の自動追加・削除。 ステップワイズ
• 標準化偏回帰係数 β*：単位の違う変数の影響度比較。 標準化β
• 調整済み R²：変数数のペナルティ込みの適合度。 調整済みR²
• AIC / BIC：モデル選択基準。 AIC / BIC
• 多重共線性／VIF：重回帰の最大の敵。 多重共線性
• 残差分析：仮定診断。 残差
• SHAP：非線形モデルでも βⱼ 的解釈をするためのツール。 SHAP
• 分位点回帰：平均ではなく分位点を推定。 分位点回帰
• ロバスト回帰：外れ値に強い推定。 ロバスト統計

🔗 関連用語（前提・並列・発展）

📚 関連グループ教材

• 📚 データ分析 — 重回帰を含む統計分析全体の地図
• 📚 統計モデリング — 回帰族の体系

📚 拡張事例 6 件 — ストーリーで読む実例

短いキャッチではなく、 各事例を詳しく追います。 数値・登場人物・教訓まで含めて再構成しました。

🎯 ケーススタディ完全版 — 47都道府県の出生率を予測する

ここまでの理論を、 「47 都道府県の出生率を 5 変数で予測する」という 1 つの完結したケーススタディとして総合的に追います。 実務でやることを順番に体験できます。

0. 課題設定 — 「なぜこの分析が必要か」

少子化対策の予算配分を最適化するために、 厚労省は「都道府県別の合計特殊出生率 (TFR) を予測するモデル」を必要としていた。 どの変数が出生率に強く影響するかが分かれば、 予算を「効きそうな施策」に重点配分できる。 そこで重回帰モデルを構築し、 政策議論のたたき台にする。

1. データ収集 — SSDSE-B-2026 から関連列を抽出

候補となる説明変数：

• x₁＝高齢化率（A1303 ÷ A1101 × 100）
• x₂＝一人当たり所得指標（経済指標として）
• x₃＝完全失業率（雇用情勢）
• x₄＝保育所定員率（H2602 ÷ A1101 × 1000）
• x₅＝1 か月の月間平均気温（生活環境）

目的変数：合計特殊出生率（自然動態統計から派生、 ここでは出生数 ÷ 15-49歳女性人口の代理として A420001 系の比率を使う）。

2. 探索的データ分析 (EDA)

まず散布図行列を描き、 単相関と分布を確認。 出生率は 1.0（東京）〜1.8（沖縄）のレンジ。 高齢化率と強い負の相関（r=−0.55）、 保育所定員率と弱い正の相関（r=+0.32）が見えた。 一方、 失業率・所得は r ≈ ±0.1 で弱い。

3. 変数選択とモデル設計

経験則 n=47 ≥ 10p より p ≤ 5。 ドメイン知識で 5 変数を採用。 VIF 事前チェックで全変数 < 3、 共線性問題なし。 重回帰式：

$TFR = \beta_0 + \beta_1 \text{高齢化率} + \beta_2 \text{所得} + \beta_3 \text{失業率} + \beta_4 \text{保育所定員率} + \beta_5 \text{気温} + \varepsilon$

4. 推定と診断

推定結果（仮想値）：β̂₁ = −0.022（p<0.001）、 β̂₂ = +0.0001（p=0.31）、 β̂₃ = −0.005（p=0.42）、 β̂₄ = +0.0008（p=0.018）、 β̂₅ = +0.012（p=0.027）。 R²=0.72, 調整済み R²=0.69, F(5,41)=21.1, p<0.001。 残差プロットでは線形性 OK、 等分散性 Breusch-Pagan p=0.18 で問題なし。 Cook's distance で東京と沖縄が高めだが、 4/n=0.085 を超えるレベルではない。

5. 解釈と政策提言

標準化β を見ると、 高齢化率（−0.55）、 保育所定員率（+0.32）、 気温（+0.27）の順。 「保育所を増やすことが、 政策で動かせる変数の中で最大効果」と読める。 ただし観察データなので「保育所を増やせば出生率が上がる」という因果効果は 断言できない。 「沖縄県は若年層が多く、 たまたま気温も高く、 保育所も多い」という同時バイアスがあり得る。 因果を立証するには、 保育所拡充の自然実験（ある自治体で予算増→他地域と比較）を待つ必要がある。

6. 報告書での書き方

エグゼクティブサマリーには「保育所定員率は出生率と正の偏連関（β=+0.0008, p=0.018, 標準化β=+0.32）。 1 標準偏差増で TFR が 0.05 ポイント上昇する関係が観察された。 ただし因果関係の立証には別途の実験的・準実験的研究が必要」と書く。 この曖昧さこそが、 観察データの正直な報告。

7. 再現可能性 — コード公開

分析コードを GitHub に公開し、 SSDSE-B-2026 のバージョンを明記。 readme.md には Python パッケージのバージョン、 実行手順、 ライセンスを記載。 Jupyter Notebook と Markdown の両方を用意し、 第三者が完全再現できる状態にする。

8. 継続改善 — 翌年度の追跡

翌年 SSDSE-B-2027 が公開されたら、 同じ式で再推定し、 係数の安定性を確認。 もし大きく動いたら、 制度変更・経済ショックなど構造変化を疑う。 これを 5 年も続ければパネルデータとして固定効果モデルに進化できる。

📚 学習リソース・参考書

独学・社内教育・書籍などの推奨リソースをまとめます。 体系的に学びたい方向け。

📘 入門書（日本語）

• 『基本統計学』宮川公男（有斐閣）— 一般教養レベル
• 『統計学入門』東京大学教養学部統計学教室編（東京大学出版会）— 「赤本」と呼ばれる定番
• 『日常生活に活かす確率・統計入門』山本義郎ほか — 重回帰実践

📘 専門書（日本語）

• 『計量経済学』西山慶彦ほか（有斐閣）— 計量経済の標準
• 『データ分析のための線形モデル』久保拓弥 — 一般化線形モデルへの橋渡し
• 『Rによる多変量解析入門』川端一光ほか — Rユーザー向け

📘 専門書（英語）

• Hastie, Tibshirani, Friedman, The Elements of Statistical Learning — ML との接続
• James et al., An Introduction to Statistical Learning（無料 PDF）— 入門最良書
• Wooldridge, Introductory Econometrics — 計量経済の標準
• Greene, Econometric Analysis — 上級者向け

🎓 オンライン講座

• Coursera: 「Statistics with R Specialization」（Duke University）
• edX: 「Data Analysis for Social Scientists」（MIT）
• YouTube: StatQuest with Josh Starmer — 直感的解説の宝庫
• 3Blue1Brown: Essence of Linear Algebra — 行列の感覚を養う

🛠 ツール・ドキュメント

• statsmodels: https://www.statsmodels.org/
• scikit-learn: https://scikit-learn.org/
• R lm()ドキュメント: ?lm
• Stan / PyMC（ベイズ拡張）

📰 学術雑誌

• Journal of the American Statistical Association (JASA)
• Journal of Econometrics
• Econometrica
• Statistics in Medicine
• 『日本統計学会誌』『応用統計学』

👥 コミュニティ

• Cross Validated (Stack Exchange) — 統計 Q&A の聖地
• r/statistics, r/econometrics
• Tokyo R 勉強会、 Japan.R
• 日本統計学会、 日本計量経済学会

❓ 追加 FAQ 10 問

Q. 重回帰分析 はどんな業界でも使えますか？

A. ほぼすべての業界で使えます。 金融・保険・医療・製造・小売・公共政策など、 数値データを扱う場面なら必ず候補に上がる手法です。 ただし業界ごとに重視される論点が違うので、 ドメイン知識との組み合わせが鍵。

Q. 結果が「思った符号と逆」になりました。 どうしたら？

A. (a) 多重共線性で符号が反転していないか VIF を確認、 (b) 単相関と重回帰の比較で交絡を疑う、 (c) ドメイン知識を再点検（自分の予想が間違っている可能性も）、 (d) 外れ値の影響を Cook's distance で確認。 結果を受け入れる前に複数の角度から検証してください。

Q. ダミー変数の「参照カテゴリ」は何を選ぶべきですか？

A. 解釈しやすい「基準」を選ぶ。 たとえば地域なら「東京」、 業種なら「製造業」など、 ドメインで標準的な参照。 サンプル数最多のカテゴリ／自然な対照を選ぶのが定石。

Q. 「予測精度」と「係数の有意性」、 どちらを優先？

A. 用途次第。 政策議論や仮説検証では係数の有意性と効果量。 自動化されたシステム（不動産査定、 信用スコア）では予測精度。 両方欲しいなら重回帰＋XGBoost のハイブリッド。

Q. SHAP は重回帰の偏回帰係数と何が違いますか？

A. 重回帰の β は「平均的・線形」な効果。 SHAP は「特定の予測における寄与」で、 非線形・相互作用も捉える。 RandomForest や XGBoost で「βⱼ 的解釈」をしたいときの代用品。

Q. 重回帰でクラスタリングはできますか？

A. 直接はできません。 重回帰は教師あり手法（y が必要）、 クラスタリングは教師なし。 ただし「残差が大きい観測グループ」を別クラスタとして扱う、 などの後処理は可能。

Q. 「効果量」って具体的に何ですか？

A. 標準化β、 部分相関係数、 R² の増分（ΔR²）など。 「統計的有意性」と独立に「実用的な大きさ」を表す指標。 Cohen の基準（小=0.1, 中=0.3, 大=0.5）が経験則。

Q. ベイズ回帰と普通の重回帰、 どう違う？

A. ベイズ回帰は係数の事前分布を仮定し、 事後分布を得る。 「不確実性を分布で表現」できるのが強み。 小標本やドメイン知識を組み込みたいときに有効。 PyMC, Stan で実装。

Q. 高次の交互作用（3 重交互作用など）を入れるべき？

A. 慎重に。 解釈が急激に難しくなり、 サンプルサイズも食う。 ドメイン知識で「意味のある 3 重相互作用」だけ入れる。 まず 2 重まで、 必要なら 3 重、 が定石。

Q. 重回帰の代わりに常に XGBoost を使えば良いのでは？

A. No。 (a) 解釈性、 (b) 不確実性の定量化、 (c) 小標本での頑健性、 (d) 規制業界での監査対応、 で重回帰が優位。 「予測精度のみ」が目的なら XGBoost / LightGBM の方が高精度なことが多いですが。

🏗 補足ケース 5 件 — 業界別の重回帰活用パターン

学術論文以外の現場での活用パターンを業界別に整理します。

金融：信用スコアリング

クレジットカード会社や消費者金融は、 デフォルト確率を年収・勤続年数・年齢・既存債務などで予測する。 線形版（重回帰系）は規制当局への説明責任が果たしやすいため、 米連邦準備制度の SR 11-7 ガイダンス対応で重宝される。 ただし最終予測は XGBoost と併走させる二重構成が増えている。

保険：料率算定

自動車保険・健康保険・生保は GLM（一般化線形モデル）で保険料を算定。 説明変数は年齢・性別・地域・運転履歴・健康診断結果など。 全ての地域・年齢クラスで料率の連続性を保つ必要があり、 線形性が好まれる。 アクチュアリーの 7 つ道具の 1 つ。

製造：歩留まり改善

半導体・自動車・鉄鋼などの製造業は、 製品の歩留まりを工程パラメータで重回帰し、 ボトルネック工程を特定。 Six Sigma の DoE（実験計画法）と密接に連携。 ANOVA・要因実験の延長として、 重回帰の知識が必須。

マーケティング：顧客生涯価値 (CLV) 予測

EC・SaaS 企業は、 顧客 1 人の生涯価値を「初回購入額」「初月利用頻度」「年代」「地域」「流入チャネル」で重回帰して予測。 これを基に獲得コスト（CAC）の上限を決め、 マーケティング投資を最適化する。 BG/NBD モデルなどとの併用が一般的。

公共政策：効果検証

補助金、 規制、 教育プログラムの効果を、 受給／非受給グループの結果差で重回帰により検証。 ただし観察データだけでは因果同定できないので、 DiD（差の差）、 RDD（不連続回帰）、 IV（操作変数）と組み合わせるのが標準。 EBPM（証拠に基づく政策立案）の中核手法。

📊 重回帰と類似手法の比較

🎓 アドバンスドトピック 10 選

重回帰分析（Multiple Regression）を実務で扱う際、 教科書には載っていない / 載っていても薄い「現場で効くトピック」をまとめます。 ここを押さえると、 ジュニア → ミドル → シニアの溝を越えられます。

1. パフォーマンスチューニングの3段階

まず 「測ってから改善」が鉄則。 推測でいじっても効果は薄い。 段階的に：

• 段階1：ボトルネックを特定する（プロファイラ、 ログ、 メトリクス）。 全体の 80% を占める処理を見つけるのが最優先（パレートの法則）。
• 段階2：アルゴリズムの計算量を下げる（O(n²) → O(n log n) など）。 ハードウェアを増強する前に試す。
• 段階3：並列化・分散化・キャッシュ・インデックス。 ここに来てからスケールアウトを検討する。

重回帰分析 の文脈でも、 たとえばデータの再計算を毎回行うのか、 部分更新で済ませるのか、 という設計判断で 10〜100 倍の性能差が出ます。 「動くものを早く作る → 計測 → 最適化」のサイクルを回す。

2. 観測可能性 (Observability) — メトリクス・ログ・トレース

本番運用では 「3 本柱」と呼ばれる 3 種類のテレメトリを揃えるのが標準：

• メトリクス (Metrics)：数値時系列。 リクエスト数、 エラー率、 レイテンシ、 リソース使用率。 Prometheus, Datadog, CloudWatch。
• ログ (Logs)：イベントの記録。 構造化ログ（JSON）が必須。 Elasticsearch + Kibana, Loki, BigQuery。
• トレース (Traces)：1 リクエストの分散処理の流れ。 OpenTelemetry, Jaeger, Zipkin。 マイクロサービスでは必須。

SLI（Service Level Indicator）・SLO（Service Level Objective）・SLA（Service Level Agreement）を定義し、 エラーバジェットを管理する SRE プラクティスが現代的標準。

3. テスト戦略 — 単体・統合・E2E のピラミッド

重回帰分析 を扱うシステムでも、 通常のソフトウェアと同じテストピラミッドが効きます：

• 単体テスト（ピラミッド底辺、 多数）：個々の関数・クラスの振る舞いを検証。 pytest, JUnit。 全体の 70% を占めるイメージ。
• 統合テスト（中段）：複数モジュールの連携。 DB、 外部 API のモック。 20%。
• E2E テスト（頂点、 少数）：本番環境に近い形で全体を通す。 Cypress, Playwright, Selenium。 10%。

データ系では追加で：データ品質テスト（great_expectations, pandera, dbt test）、 回帰テスト（モデル更新時の予測精度確認）、 負荷テスト（locust, k6）も必須。

4. CI/CD パイプライン

継続的インテグレーション (CI) と継続的デプロイ (CD) を整備すると、 変更のリスクが激減します。 標準的な段階：

1. git push → 自動 lint（flake8, eslint）
2. 単体テスト・統合テスト
3. セキュリティスキャン（Snyk, Dependabot, Bandit）
4. ビルド（Docker image）
5. ステージング環境への自動デプロイ
6. E2E テスト・スモークテスト
7. 本番環境への段階的ロールアウト（カナリア、 ブルーグリーン）
8. 本番モニタリング・アラート

5. インシデント対応（ポストモーテム文化）

障害が起きたときの対応プロセス：

• 検知：アラート受信、 顧客報告
• 初動：インシデント司令官 (IC) を 1 名指名、 Slack 専用チャンネル開設、 5 分以内にステークホルダー通知
• 緩和：症状を止める（ロールバック、 トラフィック迂回）。 根本原因が分からなくてもまず止血
• 復旧：通常運用に戻す
• ポストモーテム：原因分析・再発防止策を非難なし文化で書く。 5 Why's、 Fishbone diagram

6. データ品質の 6 次元

DAMA 国際標準では、 データ品質を 6 つの次元で測る：

• 正確性 (Accuracy)：現実と一致しているか
• 完全性 (Completeness)：欠損なし
• 一貫性 (Consistency)：異なるソース間で矛盾なし
• 適時性 (Timeliness)：必要な時点で利用可能
• 一意性 (Uniqueness)：重複なし
• 妥当性 (Validity)：定められた形式・範囲に従う

SSDSE のような公的データでも、 列の意味変更・分類体系の更新があるので、 都度バリデーションを通す。

7. プライバシー保護技術

個人情報を扱う 重回帰分析 系のシステムでは：

• 仮名化 (Pseudonymization)：個人 ID を可逆な別 ID に置換
• 匿名化 (Anonymization)：個人を特定不能に。 k-匿名性、 l-多様性
• 差分プライバシー (Differential Privacy)：統計に乱数を加え、 個人寄与を秘匿。 Apple, Google, U.S. Census が採用
• 準同型暗号 (Homomorphic Encryption)：暗号化したまま計算。 まだ重いが研究進展中
• セキュア集計 (Secure Aggregation)：複数当事者がデータを共有せずに合算結果だけを得る。 MPC
• フェデレーテッドラーニング (Federated Learning)：データを動かさず学習。 スマホ AI で利用

8. 国際規格との関係

重回帰分析 は以下の国際規格・ガイドラインと接続します：

• ISO/IEC 27001 — 情報セキュリティマネジメント
• ISO/IEC 27017/27018 — クラウドセキュリティ、 PII 保護
• ISO/IEC 38505 — データガバナンス
• NIST SP 800-53 / 800-171 — 米国政府向けセキュリティ要件
• GDPR — EU 一般データ保護規則
• 個人情報保護法（日本） — 2022 年改正で罰則強化
• CCPA / CPRA — カリフォルニア州消費者プライバシー法

監査対応では SOC 2 Type II 報告書、 ISMS 認証取得などが営業要件になることも多い。

9. 環境負荷とサステナビリティ

大規模な 重回帰分析 システムは電力消費が大きく、 CO₂ 排出が問題視されつつあります。 対策：

• 再生可能エネルギー比率の高いリージョンを選ぶ（北欧、 ケベック）
• CPU/GPU 効率の高いアーキテクチャ（ARM、 Apple Silicon）
• 不要なジョブ・データの削除（ストレージは CO₂ を消費し続ける）
• Green Software Foundation の SCI スコアを採用
• 炭素会計を経営指標に組み込む

10. キャリアパスとスキルツリー

重回帰分析 を専門にする人のキャリアパス：

• ジュニア（1-3 年）：ツール操作、 既存パイプラインの維持、 単純な分析タスク
• ミドル（3-7 年）：システム設計、 トラブル対応、 後輩指導
• シニア（7-12 年）：アーキテクチャ設計、 技術選定、 横断的な改善
• スタッフ／プリンシパル（12 年+）：会社全体の技術戦略、 標準化、 採用

学習の順序：基礎理論 → 主要ツール 1〜2 個の深掘り → 周辺ツールの広掘り → 設計パターン → 組織論。 焦らず段階的に。

✅ 運用チェックリスト

本番環境にデプロイする前に必ず通したいチェックリストです。 一項目でも飛ばすと事故率が跳ね上がります。

📋 設計フェーズ

• □ 目的と KPI が文書化されているか
• □ 想定ユースケースと想定外ユースケースが明確か
• □ 関連法規（個人情報保護法、 GDPR、 業界規制）の確認
• □ ステークホルダーとの合意形成
• □ 既存システムとのインターフェース設計
• □ データ取得元の同意・契約
• □ ロールバック・ディザスタリカバリ計画

📋 実装フェーズ

• □ コードレビュー（最低 2 名）
• □ 単体テスト カバレッジ > 80%
• □ 統合テスト
• □ セキュリティスキャン（SAST, DAST, SCA）
• □ 性能テスト（負荷、 ストレス）
• □ ドキュメント整備（README, API リファレンス, アーキテクチャ図）
• □ 機密情報のハードコード排除（git-secrets, trufflehog）
• □ ライセンス互換性確認

📋 デプロイフェーズ

• □ ステージング環境でのスモークテスト
• □ 監視・アラート設定
• □ オンコール体制の確認
• □ 段階的ロールアウト（カナリア、 ブルーグリーン）
• □ メトリクス・ログの確認
• □ ユーザー周知（リリースノート、 メール）
• □ ロールバック手順の予行演習

📋 運用フェーズ

• □ 日次の指標確認
• □ 週次の SLO レビュー
• □ 月次のキャパシティプランニング
• □ 四半期の災害復旧訓練
• □ 年次のセキュリティ監査
• □ ポストモーテムの蓄積と振り返り
• □ 継続的なドキュメント更新

📝 演習 5 問 — 理解度チェック

💥 現場の失敗例 — こうして詰んだ

❓ よくある質問 (FAQ) 10 問

Q. 重回帰と単回帰の違いは何ですか？

A. 説明変数の数が違うだけでなく、 解釈が根本的に違います。 単回帰の係数は「単独の効果」、 重回帰の偏回帰係数は「他変数を一定にしたときの純粋効果」。 同じ x₁ でも、 重回帰では別の数値になることがほとんどです。

Q. なぜ「最小二乗法」と呼ぶのですか？

A. 残差の 二乗和を最小化するからです。 絶対値和を最小化する LAD（最小絶対偏差回帰）や、 二乗和に L1/L2 ペナルティを足す Ridge/Lasso もありますが、 「普通の最小二乗 = OLS」は二乗和そのものを最小化します。 二乗にすることで微分可能・外れ値に強くペナルティ、 という性質が得られます。

Q. 説明変数は何個まで入れて良いですか？

A. 経験則 n ≥ 10p。 n=47 なら 4〜5 個まで。 機械学習文脈なら過学習を交差検証で検出できるので、 もっと多くてもよいですが、 推論を目的とするなら少なめが安全。 さらに、 入れる変数はドメイン知識で「効きそう」と思える理由があるべき。

Q. R² が低いとモデルは「使えない」ですか？

A. No. 個人レベルの予測は R²=0.3 でも統計的に意味があり、 政策評価では r=0.5 程度でも因果効果の発見に貢献する。 一方で R²=0.99 でも過学習しているだけかもしれない。 R² の絶対値より、 用途と他モデルとの比較が大事。

Q. 標準化 β と生 β、 どちらを報告すべき？

A. 両方併記が理想。 生 β は「x₁ が 1 単位増えると y がどれだけ変わるか」という解釈に直結（「築 1 年で 5 万円下落」）。 標準化 β は変数間の「影響度比較」に向く（「最も効くのは aging」）。

Q. p値が「ちょうど 0.05」のとき、 有意ですか？

A. 統計的には有意（< 0.05）。 ただし「ちょうど 0.05」は事実上の境界線。 95%CI を併記し、 「0 を含むぎりぎりではないか」を確認すべき。 また α=0.05 そのものが人為的設定なので、 機械的判断は避ける（ASA, 2016 声明）。

Q. 多重共線性があると、 常に問題ですか？

A. 解釈目的なら問題。 純粋に予測目的なら、 共線変数があっても予測精度は落ちないことが多い（Ridge ならむしろ安定）。 「βⱼ の解釈をしたいか／予測値さえ正しければ良いか」で判断。

Q. 重回帰で因果を主張できますか？

A. 原則できません。 「β₁ が有意だから x₁ が y を引き起こす」と言うには、 (i) 無作為割り付け、 (ii) 操作変数、 (iii) DiD、 (iv) RDD、 などの追加設計が必要。 観察データのみでは「条件付き連関」までしか言えない。

Q. 説明変数同士に強い相関がある（r=0.9 など）場合、 どうしたら？

A. (a) ドメイン知識で意味の重複を見抜き、 一方を落とす、 (b) PCA で主成分に変換してから回帰、 (c) Lasso／Elastic Net で自動選択、 (d) Ridge で係数を安定化、 のいずれか。 まず相関行列を可視化することが先決。

Q. オーバーフィットを防ぐ方法は？

A. (a) サンプル数を増やす、 (b) 変数を減らす（n ≥ 10p）、 (c) 正則化（Ridge / Lasso）、 (d) 交差検証で外側性能を確認、 (e) Train/Test split で実際の汎化性能を測る。 「学習データの R² だけ見る」ことは絶対 NG。

📝 レポート・論文への書き方

重回帰の結果を論文・社内報告書・ビジネス会議で報告するときの実務的ガイド。 「数値の出し方」と「言葉の選び方」の両方を整理します。

1. 数値テーブルの最低限の要素

論文で標準的に含めるべき項目（APA, AER スタイル）：

• 各説明変数の β̂（係数） と 標準誤差 SE（カッコ書きが多い）
• 有意水準のアスタリスク（*** p<0.001, ** p<0.01, * p<0.05）
• n（サンプル数）、 R²、 調整済み R²
• F値と F検定の p値
• 使った SE のタイプ（古典 / HC3 / クラスタロバスト など）
• 必要に応じて VIF と 標準化β

2. 結果セクションの文章テンプレート（学術論文）

3. ビジネス報告で避けるべき表現

• NG：「x₁ が死亡率を上げる」 → OK：「x₁ は死亡率と正の偏連関を持つ（β=0.51, p<.001）」
• NG：「効果が有意」だけ書く → OK：効果量（標準化β、 CI）も併記
• NG：「R²=0.95 だから完璧」 → OK：「観察データ上の説明力は高いが、 外挿は禁忌」
• NG：「p<0.05 だから採用」 → OK：効果量と CI を見て総合判断

4. 図表のチェックリスト

• □ 全変数の散布図行列（scatter matrix）を補遺に
• □ 残差プロット（Residuals vs Fitted）
• □ QQ プロット
• □ Cook's distance プロット
• □ 部分回帰プロット（CCPR）
• □ 標準化β の棒グラフ（影響度比較）
• □ 95%CI 付き forest plot

5. 査読者・上司に突っ込まれそうなポイント

• 「変数選択の根拠は？ ドメイン知識？ 機械的アルゴリズム？」
• 「多重共線性は確認した？ VIF を見せて」
• 「残差プロットは？ 仮定診断の証拠を」
• 「因果として解釈してない？ 観察データだから注意」
• 「外れ値の処理は？ Cook's distance チェックは？」
• 「標準誤差はロバスト？ ヘテロ性は確認した？」
• 「サンプルサイズ n と変数数 p の比は？ 経験則 10:1 守ってる？」

✅ 分析前・中・後チェックリスト

📋 分析前チェック

• □ データの出典は信頼できるか（公的統計、 論文付録など）
• □ 期間・対象範囲は研究目的と一致するか
• □ サンプルサイズ n は説明変数数 p に対して十分か（n ≥ 10p）
• □ 欠損値の扱い方針が決まっているか
• □ 単位・小数点位置・尺度（名義 / 順序 / 間隔 / 比例）を確認したか
• □ ドメイン知識的に「効きそう」な変数を漏らしていないか
• □ 因果か予測か、 用途を明確にしたか

📋 分析中チェック

• □ 散布図行列で外れ値・非線形パターンを確認
• □ 相関行列で多重共線性候補を把握
• □ VIF を全変数で計算（< 5 が理想、 < 10 が許容）
• □ 残差プロットで線形性・等分散性を診断
• □ QQ プロットで残差の正規性を診断
• □ Cook's distance で影響度の高い観測を特定
• □ ロバスト SE（HC3）の結果と通常 SE の差を確認
• □ 入れ子モデルの AIC, BIC 比較

📋 分析後チェック（報告前）

• □ R²、 調整済み R²、 F値・p値を併記
• □ 各係数の β̂、 SE、 95%CI、 標準化β を表で整理
• □ 因果として解釈していないか文章を見直す
• □ 外挿の危険性を明示
• □ コード・データを再現可能な形で公開（GitHub、 OSF）
• □ 査読・社内レビューを通す

📖 用語辞典 10 語 — この分野の必須語彙

偏回帰係数 βⱼ

他の説明変数を一定にしたときの xⱼ の効果。 重回帰の主役。

切片 β₀

全説明変数が 0 のときの予測値。 「足場」として読む。

R² / 決定係数

目的変数の分散のうちモデルで説明された割合。 0〜1。

調整済み R²

変数数のペナルティ込み。 異なる p のモデル比較に使う。

F検定

「全係数が同時に 0」を帰無仮説とする全体検定。

VIF

分散拡大因子。 多重共線性の指標。 1/(1−R²ⱼ)。

残差 eᵢ

実測 − 予測。 正規性・等分散性を診断する。

レバレッジ

観測点が回帰平面を引っ張る力。 hᵢᵢ で 0〜1。

Cook's distance

観測点を除いたとき係数がどれだけ動くか。 > 4/n で要注意。

Mallow's Cp

変数選択基準の一つ。 AIC と類似だが分散既知が仮定。

📖 拡張用語辞典 35 語

回帰係数

y = β₀ + β₁x の β₁。 単回帰では「単独の効果」、 重回帰では「条件付き効果（偏回帰係数）」。

最小二乗法

残差二乗和を最小化する推定法。 OLS。 ガウス・マルコフ仮定で BLUE。

BLUE

Best Linear Unbiased Estimator。 線形不偏推定量の中で最小分散。

ガウス・マルコフ定理

(線形・独立・等分散) のもとで OLS が BLUE。 正規性は不要。

最尤推定法 MLE

尤度を最大化。 OLS は ε∼正規 のとき MLE と一致。

正規方程式

X^T X β = X^T y。 OLS の最適性条件。

射影行列 H

H = X(X^T X)^{-1} X^T。 予測値 ŷ = Hy。

ハット行列

射影行列の別名。 対角成分 h_ii がレバレッジ。

レバレッジ

観測点が回帰平面を引っ張る力。 平均 (p+1)/n。

Cook's distance

観測点を除いたとき係数がどれだけ動くか。 > 4/n で要注意。

DFFITS

1 観測除外時の予測変化を標準化。 > 2√((p+1)/n) で注意。

DFBETAS

1 観測除外時の係数変化を標準化。 > 2/√n で注意。

スチューデント化残差

残差を標準誤差で割って正規化。 外れ値検出用。

VIF

Variance Inflation Factor。 1/(1−R²_j)。 > 10 で深刻な共線性。

Tolerance

1 − R²_j。 VIF の逆数。 < 0.1 で警戒。

F検定

全体検定。 H₀: 全係数=0。

t検定

個別係数の検定。 H₀: β_j=0。

Wald 検定

MLE の係数検定。 大標本で t と同等。

尤度比検定 LRT

入れ子モデル比較。 −2 ln(L0/L1) ∼ χ²。

AIC

赤池情報量基準。 −2 ln L + 2k。 予測志向。

BIC

ベイズ情報量基準。 −2 ln L + k ln n。 真モデル探索志向。

Mallow's Cp

RSS_p/σ² − n + 2p。 AIC と類似。

調整済み R²

1 − (1−R²)(n−1)/(n−p−1)。 変数数のペナルティ込み。

Ridge

β̂ = (X^T X + λI)^{-1} X^T y。 L2 正則化。

Lasso

‖y − Xβ‖² + λ‖β‖_1 を最小化。 L1 正則化、 自動変数選択。

Elastic Net

α‖β‖_1 + (1−α)‖β‖²_2 のペナルティ。 Ridge + Lasso。

ステップワイズ

変数を 1 つずつ追加・削除。 forward / backward / both。

交互作用項

x_i × x_j。 「x_i の効果が x_j の値で変わる」モデル化。

多項式項

x_i, x_i², x_i³ など。 非線形を線形モデルで近似。

ダミー変数

カテゴリを 0/1 で表現。 (k−1) 個を入れる（参照水準を 1 つ落とす）。

Frisch-Waugh-Lovell 定理

偏回帰係数は「他変数を除いた残差での単回帰」と一致。

Mincer 型回帰

log(賃金) = β₀ + β₁ 教育年数 + β₂ 経験 + β₃ 経験²。

内生性

E[ε|X] ≠ 0。 説明変数と誤差が相関。 OLS にバイアス。

操作変数 IV

内生性対策。 X と相関し ε と無相関の Z を使う。

HC3

Huber-White 不均一分散ロバスト SE の改良版。 小標本に強い。

📚 参考文献・規格・ドキュメント

1. Greene, W.H. (2017) Econometric Analysis, 8th ed., Pearson — 重回帰を含む計量経済学の標準教科書。
2. Wooldridge, J.M. (2019) Introductory Econometrics, 7th ed., Cengage — OLS と仮定診断の章が丁寧。
3. Hastie, T., Tibshirani, R., Friedman, J. (2009) The Elements of Statistical Learning, Springer — Ridge/Lasso と重回帰の関係。
4. 西山慶彦・新谷元嗣・川口大司・奥井亮 (2019)『計量経済学』有斐閣 — 日本語の標準教科書。
5. Belsley, Kuh, Welsch (1980) Regression Diagnostics, Wiley — 影響度指標・VIF の古典。
6. James, Witten, Hastie, Tibshirani (2021) An Introduction to Statistical Learning, 2nd ed., Springer — 入門書として最良。 無料 PDF あり。
7. Tibshirani, R. (1996) Regression Shrinkage and Selection via the Lasso, JRSS-B, 58(1): 267-288.
8. statsmodels 公式ドキュメント — statsmodels.api.OLS、 variance_inflation_factor の使い方。
9. scikit-learn User Guide — Linear Models セクション。
10. SSDSE-B-2026 データ仕様（独立行政法人統計センター）— 47 都道府県の社会経済指標。

📘 さらに学びたい人へ — 発展トピック 8 選

1. 一般化推定方程式 (GEE)

クラスタ相関のあるデータでの回帰拡張。 縦断データ・パネルでの推定に有用。

2. 混合効果モデル (LMM)

階層構造のあるデータ（生徒 in クラス in 学校）。 ランダム効果と固定効果の組み合わせ。

3. ロバスト回帰

外れ値の影響を抑える推定法。 Huber、 M-推定、 MM-推定、 LMS など。

4. 分位点回帰

平均ではなく中央値・25%点・75%点を推定。 賃金格差研究などで威力。

5. スプライン回帰 / GAM

非線形効果を滑らかな曲線で表現。 mgcv (R), pyGAM (Python)。

6. Bayesian linear regression

事前分布を入れた回帰。 不確実性を分布で表現。 PyMC, Stan で実装。

7. 因果推論ツールキット

DiD, IV, RDD, Synthetic Control, Causal Forest。 重回帰の延長線上で因果効果を識別。

8. 機械学習との融合

Double Machine Learning, Causal Trees, Meta-learners (S-, T-, X-learner)。 重回帰の解釈性と機械学習の予測性能を両立する近代的アプローチ。

📝 結びに — 重回帰を「使える人」になるために

重回帰分析は、 19 世紀後半に統計学の中核として確立された手法で、 21 世紀の AI 時代にも色褪せない基盤技術です。 機械学習の華やかなアルゴリズム（XGBoost, Deep Learning）に目を奪われがちですが、 「目的変数を線形結合で説明する」というシンプルなアイデアこそ、 データ分析の出発点であり最終的な解釈ツールです。

本ページでは、 数式・直感・実装・運用・落とし穴を 60+ セクションで掘り下げました。 ここまで読み切ったあなたは、 重回帰を「動かせる」だけでなく、 「結果を正しく解釈し、 限界を明示し、 報告できる」レベルに達したはずです。 これは統計リテラシーの一里塚です。

次のステップは：(1) ロジスティック回帰へ拡張、 (2) 因果推論（DiD・IV・RDD）の習得、 (3) ベイズ的拡張、 (4) パネルデータへの応用、 (5) 機械学習との橋渡し（SHAP・正則化）、 のいずれか。 「重回帰」を中心とした関連概念リンクをたどり、 体系的な統計学の森を歩いていってください。

📊 重回帰分析 の運用メトリクスとアラート設計

重回帰分析を本番運用するために、 監視すべきメトリクスとアラート閾値を整理します。 「測定なくして改善なし」が鉄則。

1. ゴールデンシグナル 4 種

Google SRE が提唱する 4 つの基本メトリクス：

• レイテンシ (Latency)：リクエスト処理時間。 成功と失敗で分けて測る。
• トラフィック (Traffic)：単位時間あたりリクエスト数。 容量計画の基礎。
• エラー (Errors)：失敗率。 HTTP 5xx 系、 アプリ例外、 ビジネスロジックエラー。
• サチュレーション (Saturation)：リソース使用率。 CPU、 メモリ、 ディスク、 ネットワーク。

2. RED メソッド

マイクロサービス向け 3 指標：Rate (1秒あたりリクエスト)、 Errors (失敗数)、 Duration (処理時間分布)。 Prometheus + Grafana の定番ダッシュボード。

3. USE メソッド

リソース向け 3 指標：Utilization (使用率)、 Saturation (飽和度)、 Errors (エラー)。 ハードウェアのボトルネック特定に有効。

4. アラート閾値の決め方

• SLO ベース：エラーバジェットを 50% 消費したら警告、 80% で深刻、 100% でインシデント。
• 静的閾値：CPU 80%、 ディスク 85%、 メモリ 90% で警告（業界経験則）。
• 動的閾値（ML ベース）：過去 7 日の傾向から逸脱を検出。 AWS CloudWatch Anomaly Detection。
• 多階層：警告 → 警報 → 緊急の 3 段階。 オンコール担当の睡眠を守る。

5. アラート疲労（Alert Fatigue）対策

• 誤検知（false positive）を継続的に削減
• 関連アラートのグルーピング（PagerDuty）
• ノイズの多いアラートは即廃止
• ランブック（対応手順）を必ず添付
• 四半期ごとにアラート定義を見直し

6. ダッシュボード設計の原則

• 左上に最重要メトリクス
• 1 画面 1 用途（経営層向け、 オンコール向け、 開発者向けで分ける）
• 時系列＋現在値の両方を表示
• 異常値は色で即視認
• 関連ログ・トレースへのリンクを埋め込む

7. ポストモーテム文化

インシデント後の振り返り。 「Blameless（非難なし）」原則で実施。 タイムライン、 根本原因、 再発防止策、 学んだことを文書化し、 組織全体で共有。 同じ事故を二度繰り返さないためのナレッジマネジメント。

8. SLI/SLO/SLA の関係

• SLI (Service Level Indicator)：実測値。 例：成功率 99.95%。
• SLO (Service Level Objective)：内部目標。 例：成功率 99.9%。 SLI が下回ったら改善対象。
• SLA (Service Level Agreement)：顧客との契約。 例：成功率 99.5%、 違反時は返金。 SLO より緩く設定。

9. エラーバジェット運用

SLO 99.9% なら月 43.2 分のダウン許容。 これがエラーバジェット。 余裕がある月は新機能リリースを許可、 使い切ったら凍結し信頼性向上に集中。 開発と運用の利益相反を解消する仕組み。

10. ロギング戦略

• 構造化ログ（JSON）が必須
• トレース ID で 1 リクエストの全ログを横断検索
• 個人情報・機密情報のマスキング
• ログレベル（DEBUG/INFO/WARN/ERROR）の使い分け
• 保管期間：監査用 7 年、 運用用 30 日、 デバッグ用 7 日
• コスト：S3 / Cloud Storage のライフサイクル管理で安価階層へ

🖼 概念図で押さえる 重回帰

重回帰 (Multiple Regression) は「複数の説明変数で 1 つの目的変数を予測」する基本モデル。 ここでは 3 つの概念図で、 単回帰との違い・多重共線性の影響・係数の標準化を視覚的に整理する。

📌 図から読み取るポイント (3 行要約)

• 図 1: 重回帰は超平面の当てはめ。 説明変数を増やすほど自由度は減るので調整済み R² で評価する。
• 図 2: VIF>10 は要除外。 個別係数の解釈が必要なら多重共線性チェックを必ず実施する。
• 図 3: 標準化 β で寄与度比較。 元の係数は単位依存なので、 影響力ランキングには β を使う。

🗺 概念マップ

重回帰を中心に、 上位・下位・並列の概念をツリーで配置：

一般化線形モデル (GLM)
├─ 線形回帰
│  ├─ 単回帰 (p=1)
│  ├─ 重回帰 (p≥2)  ← この用語
│  │  ├─ OLS                — 通常の最小二乗
│  │  ├─ WLS / GLS          — 重み付き／一般化最小二乗
│  │  ├─ Ridge              — L2 正則化
│  │  ├─ Lasso              — L1 正則化
│  │  └─ Elastic Net        — L1+L2
│  └─ 多項式回帰
├─ ロジスティック回帰 (二値)
├─ ポアソン回帰 (カウント)
└─ プロビット回帰

拡張：
├─ GAM (非線形)
├─ 混合効果モデル (階層構造)
├─ パネル固定効果モデル
└─ Bayesian Regression

仮定が崩れた時：
├─ HC3 / Newey-West SE (等分散・自己相関違反)
├─ ロバスト回帰 (外れ値)
├─ 操作変数法 IV (内生性)
└─ DiD / RDD (因果推論)

📜 歴史と系譜 — どうしてこの概念が生まれたか

重回帰の系譜は意外に古く、 18 世紀の天文学にまで遡ります。

1. 1805: ルジャンドル — 最小二乗法の発明

アドリアン＝マリ・ルジャンドルが彗星軌道計算のために 最小二乗法（Méthode des moindres carrés）を発表（1805 年）。 同時期にガウスも独立に発見しており、 「ガウス・ルジャンドル法」とも呼ばれる。 これが重回帰の数学的基礎。

2. 1885-1900: ゴルトン → ピアソン → ユール

フランシス・ゴルトン（1885）が遺伝学の研究で 「回帰（regression to the mean）」という言葉を作る。 カール・ピアソン（1896）が相関係数を厳密化し、 ユードニー・ユール（1907）が 偏回帰係数（partial regression coefficient）を定義。 ここで現代的な重回帰の枠組みが完成。

3. 1930s-50s: フィッシャー、 マルコフ — 推論の理論化

R.A. フィッシャー（1922）が最尤法と F検定を確立。 ガウス・マルコフの定理が整理され、 OLS の BLUE 性質が証明される。 これにより重回帰は「推定 + 検定」の体系として完成した。

4. 1960s-80s: 計量経済学の興隆

クリーヴ・グレンジャー、 ジェームズ・ヘックマンらにより、 内生性・自己相関・パネルデータといった現実データの問題点が体系化される。 ヘックマンの選択バイアス補正、 ヘンドリーの一般から特殊への方法論など、 計量経済学の現代スタイルが確立。

5. 1996: Lasso と正則化革命

ロバート・ティブシラニが Lasso（L1 正則化）を発表（Tibshirani, 1996）。 高次元データへの対応と自動変数選択を可能にし、 機械学習との接続が始まる。 Elastic Net（Zou & Hastie, 2005）で完成形。

6. 2000s-現在: 機械学習との融合

scikit-learn、 statsmodels、 PyMC、 Stan などの普及で、 重回帰は 「全データサイエンティストの基本ツール」になった。 SHAP・LIME などの解釈ツールで非線形モデルにも「βⱼ 的解釈」が広がる。 重回帰は古いが、 現代でも色褪せない。

7. 日本での発展

日本では 1950 年代に統計数理研究所の北川敏男・赤池弘次らが計量経済学を発展させ、 AIC（赤池情報量基準、 1973）を生み出した。 これは重回帰のモデル選択における世界標準。 SSDSE データを使う本コンペは、 こうした半世紀以上の伝統の上に立つ。

🔧 実装詳細 — 一段深く潜る

「OLS で重回帰」と一言で言っても、 実装にはいくつもの細部があります。 ここでは「教科書には載りにくいが現場で大事」な詳細を整理します。

1. 数値計算の安定性 — なぜ linalg.solve を使うべきか

公式 $\hat{\beta} = (X^\top X)^{-1} X^\top y$ をそのままコード化するなら inv(X.T @ X) @ X.T @ y ですが、 これは数値的に不安定。 条件数（condition number）が大きいと丸め誤差で結果が壊れる。 推奨は：

• linalg.solve(X.T @ X, X.T @ y) — LU 分解で解く（やや改善）
• linalg.lstsq(X, y) — SVD ベース、 最も安定（推奨）
• QR 分解：$X = QR$、 $\hat{\beta} = R^{-1} Q^\top y$ — statsmodels が内部で使用

2. ダミー変数の扱い — pandas の get_dummies

カテゴリ変数（地域、 業種など）は one-hot エンコーディングで扱う。 ただし注意：

• drop_first=True で必ず 1 つ落とす（さもないと完全な多重共線性に）。
• 落とした水準が 参照カテゴリ（reference category）になる。 結果の解釈はこれに対する差分。
• カテゴリ数が多い変数（例：47 都道府県）は ダミー数が爆発するので、 ターゲットエンコーディングや埋め込みを検討。

3. 交互作用項（interaction term）の入れ方

「x₁ の効果が x₂ の値で変わる」場合、 x₁ × x₂ を新変数として追加する。 statsmodels では 'y ~ x1 * x2' が自動で x1, x2, x1:x2 の 3 つを入れてくれる。 解釈は：

• x₁ の主効果 β₁：「x₂=0 のときの x₁ の効果」
• x₂ の主効果 β₂：「x₁=0 のときの x₂ の効果」
• 交互作用 β₁₂：「x₂ が 1 増えたとき、 x₁ の効果が何ポイント変わるか」

4. 多項式項 — 非線形を線形モデルで扱う

y と x₁ の関係が U 字なら、 x₁ + x₁² を入れる。 経験年数と賃金の Mincer 型回帰など、 経済学で頻出。 注意：x₁ と x₁² は強く相関するので 中心化（x₁ − x̄）してから二乗するのが推奨。 sklearn なら PolynomialFeatures。

5. 重み付き OLS（WLS）

観測値ごとに信頼度が違うときは重み w_i を付ける。 重みは「分散の逆数」が理論的に最適：$w_i = 1 / \sigma_i^2$。 ヘテロスケダスティシティへの対処にも使える。 smf.wls(formula, data=df, weights=w)。

6. 欠損値の扱い

SSDSE データは時々欠損があります。 基本戦略：

• Listwise deletion：1 個でも欠損のある行を全削除（default）。 サンプル損失大。
• 平均代入：簡単だが分散を過小評価。
• 多重代入（multiple imputation）：sklearn の IterativeImputer、 R の mice パッケージ。 統計的に最も望ましい。

7. 標準誤差の選択肢一覧

⚙️ 運用とトラブルシュート

本番運用での「壊れ方」と対処を整理。

1. データドリフトの監視

学習時と本番で説明変数の分布が変わると β̂ は陳腐化する。 月次で平均・分散・KL ダイバージェンスを取り、 5% を超える変動があれば再学習。 不動産モデルなら「平均坪単価のシフト」、 マーケティングなら「広告予算の構成変化」が典型。

2. モデル劣化の検知

予測誤差（RMSE、 MAE）を本番ログから継続計算。 学習時の CV RMSE の 1.5 倍を超えたらアラート。 また、 残差の分布が変化していないか、 自動 KS 検定で確認する。

3. パイプライン障害

SSDSE データのような外部データを毎年更新するパイプラインは、 列名変更・コードの追加・欠損値の増加で簡単に壊れる。 schema validation（great_expectations、 pandera）を組み込み、 異常があれば deploy をブロックする。

4. 外れ値の混入

単一の外れ値で β が大きく動くことを利用したデータポイズニング攻撃が存在する。 本番では Z-score > 5 の入力を自動拒否、 IQR ベースのチェックを通すなど、 入力検証を厳密に。

5. レスポンス時間と再学習頻度

6. ロールバック戦略

モデルは MLflow などで version 管理し、 旧版に即座に戻せるようにする。 A/B テストで段階リリース。 「全トラフィックの 1% で新モデル → 異常がなければ 10% → 50% → 100%」のカナリアリリース。

💴 コストと見積もり

重回帰そのものの計算コストは極めて低い（O(np² + p³)）が、 周辺コストが実は本番では支配的。

1. 学習コスト

2. 推論コスト

予測は p+1 個の乗算と足し算だけ。 1 件 ≈ ナノ秒オーダー。 1 億件でも 1 秒以内で処理可能。 これが「重回帰はインフラ的に楽」と言われる所以。

3. 人的コスト（隠れたボトルネック）

実は重回帰の総コストの 80% は 「データ前処理」と 「結果の解釈・報告」。 計算コストは無視できる。 アナリストが半日〜数日かけて：

• 変数の意味確認とドメインインタビュー
• 欠損・外れ値処理
• VIF・残差プロットチェック
• レポート作成（係数・CI・標準化β・診断）

4. 教育コスト

重回帰を「正しく」使いこなすには、 大学レベルの統計学（OLS、 仮定診断、 多重共線性、 因果推論）を最低限学ぶ必要がある。 自社で internal training を入れる場合、 1 人あたり 20〜40 時間の研修コスト（社員時給 5,000 円換算で 10〜20 万円）。 e-learning や AtCoder / KaggleでのSkill Tier 化が有効。

5. ライセンス・規制コスト

scikit-learn (BSD)、 statsmodels (BSD)、 R (GPL) など主要ツールはオープンソース無料。 SAS や SPSS の商用利用は年数十万〜数百万円。 ただし金融や医療の規制業界では SAS validation の方が監査で楽、 という現実もある。

🛡 ガバナンス・セキュリティ・コンプライアンス

重回帰モデルを業務に使う際のガバナンス論点を整理。

1. 説明責任（Accountability）

重回帰の最大の強みは 「なぜこの予測か」を係数で説明できること。 これは EU の GDPR Article 22（自動化された意思決定の説明権）にも適合する。 顧客に対して「あなたのスコアが下がった理由は、 βⱼ × xⱼ の項目 X が原因」と回答可能。

2. 公平性（Fairness）

説明変数に「性別」「人種」が入っていなくても、 「住所」「学歴」を経由して差別が忍び込む（disparate impact）。 EEOC や金融庁の公平性ガイドラインに照らし、 protected attribute での予測精度差をモニタリングする。

3. 個人情報保護

SSDSE-B-2026 は 47 都道府県の集計値なので個人情報問題は基本なし。 ただし、 マイクロデータ（個票）で重回帰する場合は匿名化（k-匿名性、 差分プライバシー）を必ず施す。 学習データの再特定化攻撃も検討事項。

4. モデルカード（Model Card）

Google が提唱するモデルカードを作成すべき項目：

• モデルの目的・想定ユースケース・想定外の用途
• 学習データの出典・期間・サンプル数
• 性能指標（R², RMSE, 95%CI）と分布別性能（地域別、 属性別）
• 既知の限界（外挿不可、 仮定違反のケース）
• 更新履歴と監査ログ

5. 監査証跡（Audit Trail）

いつ・誰が・どのデータで・どの結果を出したかを完全記録。 MLflow Tracking、 DVC、 Weights & Biases などのツールで自動化。 金融や医療では Part 11（FDA）/ J-SOX 対応が必須。

6. AI ガイドラインへの対応

日本では総務省・経産省の AI 利活用ガイドライン、 政府の AI 戦略 2022。 重回帰のような「線形・透明」なモデルは high-risk 分類に該当しないことが多いが、 信用スコアや採用予測など用途次第で規制対象に。 詳細は AI ガイドライン 参照。

🏭 産業事例 6 件 — 現場ではどう使われているか

🔗 隣接手法への橋渡し

「重回帰分析」は単独で完結する手法ではなく、 隣接領域と連携することで真価を発揮する。

• 上流: 単回帰・相関分析・標準化
• 並列: リッジ回帰・Lasso・GAM
• 下流: モデル選択・多重共線性診断・残差分析

SSDSE-B-2026 で県民所得を被説明変数、 人口・高齢化率・就業者数を説明変数とする重回帰は、 多重共線性 (人口と就業者) と外れ値 (東京) の両方を考慮する必要がある。 標準化偏回帰係数で比較するのが安全。

🌳 意思決定ツリー — 状況別の手順

どの正則化を使うべきか？ どんな診断を回すべきか？ — 状況別の意思決定ツリーをまとめます。

🌳 ツリー 1: 重回帰の出発点

目的変数 y は何か？
├─ 連続値 → 重回帰（このページ）／GAM／RandomForest
├─ 二値 (0/1) → ロジスティック回帰
├─ カウント (0,1,2,...) → ポアソン回帰／負の二項回帰
├─ カテゴリ (3つ以上) → 多項ロジット
└─ 順序 (低/中/高) → 順序ロジット
    

🌳 ツリー 2: 説明変数選択

変数数 p と n の関係は？
├─ n ≥ 10p → そのまま OLS で OK
├─ n < 10p かつ予測重視 → Lasso / Elastic Net
├─ n << p (高次元) → Lasso 強推奨、 もしくは PCA + 回帰
└─ 共線性あり (VIF > 10) → Ridge / 変数除外 / PCA
    

🌳 ツリー 3: 仮定診断

残差プロットを見て：
├─ ランダム雲 → 仮定 OK
├─ ファン状の広がり → ヘテロ性、 HC3 SE か WLS
├─ 曲線パターン → 非線形効果、 多項式項か GAM
├─ 周期パターン → 自己相関、 Newey-West SE か ARIMA
└─ 一部の点だけ外れ → 外れ値、 Cook's d 確認
    

🌳 ツリー 4: 因果推論への進化

因果を主張したい？
├─ 無作為割付できる → RCT（実験）
├─ 自然実験あり → IV / DiD / RDD
├─ パネルデータあり → 固定効果モデル
├─ 観察データのみ → 「条件付き連関」までしか言えない
└─ → 報告書に「因果ではない」と明記
    

🌳 ツリー 5: モデル比較

何を基準に選ぶ？
├─ 予測重視 → 交差検証 RMSE、 CV R²
├─ 説明力重視 → 調整済み R²
├─ 「真モデル」探索 → BIC
├─ 予測 & スパース → Lasso CV
└─ 入れ子モデル比較 → F 検定 / 尤度比検定
    

📖 本コンペ過去論文での重回帰例

本コンペで再現対象となる過去論文では、 重回帰が中心的な手法として頻出します。 代表的なものを 5 件紹介します（あくまで例。 実際の論文タイトルは本リポジトリの index.html を参照）。

論文例 1: 「47 都道府県の医療費を 7 変数で説明する」

目的：医療費（人口あたり）を高齢化率・所得・病院数・医師数・気候・産業構造・教育水準で重回帰。 結果：高齢化率と医師数が主要因。 ただし医師数と医療費の因果方向に注意が必要（医療費が高いから医師が来るのか、 医師が多いから医療費が高くなるのか）。

論文例 2: 「市町村の人口減少率を 5 変数で予測する」

目的：5 年間の人口減少率を、 高齢化率・有効求人倍率・出生率・転出率・大学進学率で重回帰。 結果：転出率（特に若年層）が支配的。 政策提言：地方創生は「若年層の流出抑制」が最大のレバー。

論文例 3: 「気象データと農業生産量の重回帰」

目的：米の単収（10a あたり）を平均気温・日照時間・降水量・台風数・施肥量で重回帰。 結果：気温と日照が主要因。 ただし「気温」は二次関数的効果（暑すぎても寒すぎても減収）なので多項式項を入れた。

論文例 4: 「教育格差の地域要因分析」

目的：高校進学率を世帯所得中央値・両親学歴・教育費補助・地域人口で重回帰。 結果：所得と両親学歴が支配的。 政策提言：奨学金拡充の効果量を試算する基礎データに。

論文例 5: 「自殺率の社会経済要因」

目的：年齢調整自殺率を失業率・離婚率・所得格差・社会的孤立指標・寒冷地ダミーで重回帰。 結果：失業率と社会的孤立が主要因。 教訓：「自殺は個人の問題」ではなく社会的決定因子の影響大。 ただし因果は別途実証必要。

🏋 ワークショップ 5 ワーク — 自分で手を動かす

頭で読むだけでは身につきません。 ここでは「自分で手を動かす演習」を 5 ワーク用意します。 SSDSE-B-2026 を使う前提です。 各ワークには想定解と頻出ミスを併記しました。

📝 ワーク 1: 単純な 2 変数重回帰を組む

課題：47 都道府県の死亡率を「高齢化率」と「年平均気温」だけで重回帰せよ。 β̂、 SE、 R² を報告すること。

想定解：高齢化率の β = 0.49（p<.001）、 気温 β = −0.03（p=0.27）、 R²=0.946。 単相関では気温が r=−0.68 だが、 重回帰では非有意。 これは「寒冷地は高齢化が進む」交絡が原因。

頻出ミス：encoding を指定し忘れて文字化け（cp932 を使う）、 skiprows=[1] を忘れて 1 行目を欠損として読む。

📝 ワーク 2: 多重共線性を意図的に起こしてみる

課題：「高齢化率」と「65歳以上人口数」（A1303）を両方説明変数に入れて重回帰せよ。 VIF を計算して結果を解釈すること。

想定解：両者は定義式上 100% 相関なので、 VIF が数百〜数千になる（実質的には推定不能）。 statsmodels はランク落ちを検出して片方の係数を NaN にするか、 警告を出す。 教訓：派生変数を作るときは「片方しか入れない」ように設計する。

📝 ワーク 3: Ridge と Lasso の比較

課題：同じ 4 変数モデルを Ridge と Lasso で推定し、 α を 0.001〜10 まで変化させて係数の軌跡をプロットせよ。

想定解：Ridge は α を上げても全係数が滑らかに 0 に近づく。 Lasso は α=0.1 あたりで β₂, β₃, β₄ が完全にゼロになり、 「実質高齢化率 1 変数モデル」になる。 これが Lasso のスパース性。

📝 ワーク 4: 交差検証で汎化性能を測る

課題：4 変数重回帰の CV R² を、 5-fold と LOO（Leave-One-Out）で比較せよ。

想定解：5-fold CV R² ≈ 0.92（学習時 R² = 0.948 より少し下がる）。 LOO CV R² も同程度。 過学習はほぼ起きていない。 もし大きく乖離していたら変数を絞るべき。

📝 ワーク 5: 残差プロットで仮定診断

課題：4 変数モデルの (1) Residuals vs Fitted, (2) QQ plot, (3) Scale-Location, (4) Cook's distance を 4 枚プロットせよ。 各図から読み取れることを 1 文ずつ書け。

想定解：(1) ランダム雲、 線形性 OK。 (2) ほぼ直線上、 正規性 OK。 (3) 横軸に対して平坦、 等分散性 OK。 (4) Cook's d は全て < 4/47=0.085、 影響度の高い外れ値なし。 → 全仮定支持。