← 目次 ｜ ← 前: 09 ｜ 次: 11 →

10. 重回帰分析（行列演算）

テーマ: 賃金関数の推定（教育年数・経験年数・週労働時間）

学習点: NumPy の行列演算, 正規方程式, 分散共分散行列, VIF, F検定

依存: NumPy, SciPy ／ 難易度: 中級

実行方法

uv run 10_ols_multiple.py

スクリプト冒頭の PEP 723 メタデータ（# /// script）により、必要なライブラリは uv が自動的に仮想環境へ導入します。事前の pip install は不要です。

解説

何をするプログラムか

労働経済学の定番である賃金関数（Mincer 型）を重回帰で推定するスクリプトです。教育年数・経験年数・経験年数の 2 乗・週労働時間の 4 変数で賃金の対数 ln(wage) を説明します。データは真のモデル（教育 1 年あたり +8.5% など）に正規乱数の誤差を加えて 200 人分を人工的に生成しているため、推定値が真の値をどの程度復元できるかを答え合わせできます。

前章 09 が総和記号で書ける単回帰だったのに対し、本章は行列演算で任意の本数の説明変数を一気に扱います。係数の推定・標準誤差・t 検定・F 検定に加え、説明変数どうしの相関の強さを測る VIF（分散拡大要因）も計算する、計量経済学の実習に近い内容です。

コードの読みどころ

• X = np.column_stack([np.ones(n), educ, exper, exper ** 2, hours]) で、定数項の 1 の列と各説明変数を並べた 200×5 の計画行列を作ります。exper ** 2 のように配列全体のべき乗が 1 式で書けるのが NumPy です。
• ols() の中核は beta = XtX_inv @ X.T @ y の 1 行で、正規方程式の解 β̂ = (X’X)⁻¹X’y を @ 演算子（行列積）でそのまま表現しています。説明変数が何本に増えてもこのコードは変わりません。
• 標準誤差は np.sqrt(np.diag(sigma2 * XtX_inv)) で求めます。分散共分散行列 σ²(X’X)⁻¹ の対角成分が各係数の分散になる、という理論がコードに直結しています。
• rng = np.random.default_rng(42) で乱数シードを固定しているため、誰が実行しても同じ推定結果が再現されます。数値実験の基本作法です。
• vif() は各説明変数を残りの変数に回帰した R² から VIF = 1/(1−R²) を計算します。列の削除に np.delete(X, j, axis=1)、補助回帰の求解に np.linalg.lstsq を使っており、「VIF は補助回帰の当てはまりの良さ」という定義がコードから読み取れます。
• 2 次項を入れたモデルでは、peak = -beta[2] / (2 * beta[3]) で経験年数の限界効果がゼロになる頂点（放物線の頂点の公式）を計算しています。

理論的背景

重回帰の OLS 推定量は正規方程式の解 β̂ = (X’X)⁻¹X’y で、その分散は Var(β̂) = σ²(X’X)⁻¹、σ² は RSS/(n−k) で推定します（k は定数項込みのパラメータ数）。F 検定は「定数項以外の係数がすべてゼロ」という帰無仮説をモデル全体で検定します。被説明変数が対数のため、係数は近似的に「説明変数 1 単位あたりの変化率」と読めます（半対数モデル）。VIF_j = 1/(1−R²_j) は多重共線性の指標で、一般に 10 を超えると係数の標準誤差が大きく膨らんでいる可能性を警戒します。

実行結果の見方

教育年数の係数 0.08189 は真の値 0.085 にほぼ一致し、t = 8.69 で 1% 水準でも有意です（教育 1 年で賃金約 8.2% 上昇）。一方、経験年数（p = 0.0743）と経験年数²（p = 0.3193）は単独では有意性が弱く、これは VIF がそれぞれ 11.77・11.88 と高いこと、つまり 2 つの変数がほぼ同じ情報を持ち標準誤差が膨らんでいることと表裏です。ただし末尾の注記のとおり、2 次項による設計上の相関なので直ちに問題ではありません。モデル全体では F(4, 195) = 22.54、p = 2.554e-15 と強く有意で、R² = 0.3162 は個人データの賃金回帰としては標準的な水準です。経験年数の効果のピークは 28.8 年と、賃金プロファイルが逓増から逓減へ転じる年数も推定できています。

ソースコード

# /// script
# requires-python = ">=3.11"
# dependencies = [
#     "numpy",
#     "scipy",
# ]
# ///
"""10: 重回帰分析（行列演算） ---------------------------------------------
テーマ: 賃金関数の推定（教育年数・経験年数・週労働時間）
学習点: NumPy の行列演算, 正規方程式, 分散共分散行列, VIF, F検定
根拠: β̂ = (X'X)^{-1} X'y、 Var(β̂) = σ²(X'X)^{-1}、σ² = RSS/(n-k)。
      VIF_j = 1/(1-R²_j) は多重共線性の指標で、一般に 10 超で警戒。
"""
import numpy as np
from scipy import stats

rng = np.random.default_rng(42)   # 再現性のため乱数シードを固定


def ols(X: np.ndarray, y: np.ndarray):
    n, k = X.shape                       # k は定数項込みのパラメータ数
    XtX_inv = np.linalg.inv(X.T @ X)
    beta = XtX_inv @ X.T @ y
    resid = y - X @ beta
    rss = float(resid @ resid)
    sigma2 = rss / (n - k)
    se = np.sqrt(np.diag(sigma2 * XtX_inv))
    tss = float(((y - y.mean()) ** 2).sum())
    r2 = 1 - rss / tss
    adj = 1 - (1 - r2) * (n - 1) / (n - k)
    f = (tss - rss) / (k - 1) / sigma2
    return beta, se, r2, adj, f, n, k


def vif(X: np.ndarray, names) -> None:
    """各説明変数を他の説明変数に回帰して VIF を計算する。"""
    print("\n[多重共線性 VIF]")
    for j in range(1, X.shape[1]):           # 0 列目は定数項
        others = np.delete(X, j, axis=1)
        b = np.linalg.lstsq(others, X[:, j], rcond=None)[0]
        r = X[:, j] - others @ b
        r2 = 1 - (r @ r) / ((X[:, j] - X[:, j].mean()) ** 2).sum()
        print(f"  {names[j]:<10} VIF = {1 / (1 - r2):>6.2f}")


def main() -> None:
    n = 200
    educ = rng.normal(14, 2.2, n).clip(9, 22)
    exper = rng.normal(15, 7, n).clip(0, 40)
    hours = 40 + 0.3 * exper + rng.normal(0, 4, n)   # 経験と相関する変数
    # 真のモデル: ln(wage) = 1.5 + 0.085*educ + 0.02*exper - 0.0003*exper^2 + u
    lnwage = (1.5 + 0.085 * educ + 0.02 * exper - 0.0003 * exper ** 2
              + 0.004 * hours + rng.normal(0, 0.25, n))

    names = ["定数項", "教育年数", "経験年数", "経験年数^2", "週労働時間"]
    X = np.column_stack([np.ones(n), educ, exper, exper ** 2, hours])
    beta, se, r2, adj, f, n_, k = ols(X, lnwage)
    t = beta / se
    p = 2 * (1 - stats.t.cdf(np.abs(t), n_ - k))

    print("=== 重回帰: ln(wage) の推定 ===")
    print(f"{'変数':<12}{'係数':>11}{'標準誤差':>11}{'t値':>9}{'p値':>10}  ")
    print("-" * 56)
    for nm, b, s, tv, pv in zip(names, beta, se, t, p):
        star = "***" if pv < 0.01 else "**" if pv < 0.05 else "*" if pv < 0.1 else ""
        print(f"{nm:<12}{b:>11.5f}{s:>11.5f}{tv:>9.2f}{pv:>10.4f}  {star}")
    print("-" * 56)
    print(f"n = {n_},  R² = {r2:.4f},  調整済R² = {adj:.4f}")
    print(f"F({k-1}, {n_-k}) = {f:.2f},  p = {1 - stats.f.cdf(f, k-1, n_-k):.3e}")
    print(f"\n解釈: 教育年数1年増で賃金は約 {beta[1]*100:.1f}% 上昇"
          "（半対数モデルの係数は近似的に変化率）。")

    # 経験年数の限界効果がゼロになる年数（頂点）
    peak = -beta[2] / (2 * beta[3])
    print(f"経験年数の効果のピーク: {peak:.1f} 年（以降は逓減）")

    vif(X, names)
    print("  ※ 経験年数と経験年数^2 は定義上強く相関するため VIF は高くなるが、"
          "\n     これは設計上の相関であり推定の妥当性を直ちに損なうものではない。")


if __name__ == "__main__":
    main()

実行結果

=== 重回帰: ln(wage) の推定 ===
変数                   係数       標準誤差       t値        p値  
--------------------------------------------------------
定数項             1.34301    0.24548     5.47    0.0000  ***
教育年数            0.08189    0.00943     8.69    0.0000  ***
経験年数            0.01576    0.00878     1.79    0.0743  *
経験年数^2         -0.00027    0.00027    -1.00    0.3193  
週労働時間           0.00960    0.00450     2.14    0.0339  **
--------------------------------------------------------
n = 200,  R² = 0.3162,  調整済R² = 0.3022
F(4, 195) = 22.54,  p = 2.554e-15

解釈: 教育年数1年増で賃金は約 8.2% 上昇（半対数モデルの係数は近似的に変化率）。
経験年数の効果のピーク: 28.8 年（以降は逓減）

[多重共線性 VIF]
  教育年数       VIF =   1.02
  経験年数       VIF =  11.77
  経験年数^2     VIF =  11.88
  週労働時間      VIF =   1.33
  ※ 経験年数と経験年数^2 は定義上強く相関するため VIF は高くなるが、
     これは設計上の相関であり推定の妥当性を直ちに損なうものではない。

---

← 目次 ｜ ← 前: 09 ｜ 次: 11 →