機械学習と情報技術

重回帰分析の理論と最小二乗法の導出

多変量解析

住宅の価格は何によって決まるのでしょうか。面積が広ければ高い、駅に近ければ高い、築年数が浅ければ高い — 直感的にはいくつかの要因が思い浮かびます。しかし、面積が10平米増えると価格はいくら上がるのか、築年数が1年増えると価格はいくら下がるのか、これらの影響を同時に定量化するにはどうすればよいでしょうか。

重回帰分析（Multiple Regression Analysis）は、複数の説明変数から1つの目的変数を予測する統計手法です。単回帰が1つの説明変数しか扱えないのに対し、重回帰は複数の要因の個別の効果を同時に推定できます。「他の条件が等しいとき、面積が1平米増えると価格はいくら変わるか」といった偏回帰係数の解釈は、重回帰分析ならではの強みです。

重回帰分析を理解すると、以下のような応用が開けます。

  • 経済学: 消費支出を所得・金利・物価で説明する計量経済モデル
  • 疫学: 疾患リスクを年齢・喫煙・運動量で説明し、各リスク要因の独立した効果を推定
  • 工学: 製品の品質を製造条件（温度・圧力・速度）で予測する実験計画
  • マーケティング: 売上を広告費・価格・季節要因で説明する

本記事では、重回帰モデルを行列表現で定式化し、最小二乗推定量を正規方程式から導出します。推定量の統計的性質、決定係数、偏回帰係数の解釈、仮説検定までを体系的に解説し、Pythonで実装します。

本記事の内容

  • 重回帰モデルの行列表現
  • 正規方程式と最小二乗推定量の導出
  • ガウス=マルコフの定理（BLUE性）
  • 決定係数と自由度調整済み決定係数
  • 偏回帰係数の解釈と偏相関との関係
  • Pythonでの実装と可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

重回帰モデルの定式化

単回帰から重回帰へ

単回帰モデル $y = \beta_0 + \beta_1 x + \epsilon$ は1つの説明変数 $x$ で目的変数 $y$ を予測します。しかし現実の多くの問題では、結果に影響する要因は複数あります。住宅価格は面積だけでなく、築年数、駅からの距離、階数など、多くの要因に依存します。

重回帰モデルは、$p$ 個の説明変数 $x_1, x_2, \dots, x_p$ を使って目的変数 $y$ を予測します。

$$ \begin{equation} y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_p x_p + \epsilon \end{equation} $$

$\beta_0$ は切片、$\beta_j$（$j = 1, 2, \dots, p$）は偏回帰係数（partial regression coefficient）です。偏回帰係数 $\beta_j$ は「他の説明変数を一定に保ったとき、$x_j$ が1単位増加すると $y$ が $\beta_j$ だけ変化する」ことを意味します。

行列表現

$n$ 個の観測データ $(\bm{x}_i, y_i)$（$i = 1, 2, \dots, n$）に対して、モデルを行列で書きます。

$$ \begin{equation} \bm{y} = \bm{X}\bm{\beta} + \bm{\epsilon} \end{equation} $$

各行列・ベクトルの定義は以下の通りです。

$$ \bm{y} = \begin{pmatrix}y_1\\y_2\\\vdots\\y_n\end{pmatrix}, \quad \bm{X} = \begin{pmatrix}1 & x_{11} & x_{12} & \cdots & x_{1p}\\1 & x_{21} & x_{22} & \cdots & x_{2p}\\\vdots & \vdots & \vdots & \ddots & \vdots\\1 & x_{n1} & x_{n2} & \cdots & x_{np}\end{pmatrix}, \quad \bm{\beta} = \begin{pmatrix}\beta_0\\\beta_1\\\vdots\\\beta_p\end{pmatrix}, \quad \bm{\epsilon} = \begin{pmatrix}\epsilon_1\\\epsilon_2\\\vdots\\\epsilon_n\end{pmatrix} $$

$\bm{X}$ は $n \times (p+1)$ の計画行列（design matrix）で、第1列はすべて1（切片に対応）です。

誤差項の仮定

最小二乗法の導出には誤差項について次の仮定を置きます。

$$ \begin{equation} E[\bm{\epsilon}] = \bm{0}, \quad \text{Cov}(\bm{\epsilon}) = \sigma^2\bm{I}_n \end{equation} $$

つまり、各誤差の期待値はゼロ、分散は共通の $\sigma^2$（等分散性）、異なるデータ点の誤差は無相関です。正規分布の仮定 $\epsilon_i \sim N(0, \sigma^2)$ を追加すると、仮説検定が可能になります。

この定式化をもとに、パラメータ $\bm{\beta}$ を推定しましょう。

最小二乗推定量の導出

残差二乗和の最小化

最小二乗法は、残差二乗和（RSS: Residual Sum of Squares）を最小化するパラメータ $\hat{\bm{\beta}}$ を求めます。

$$ \begin{equation} \text{RSS}(\bm{\beta}) = \sum_{i=1}^{n}(y_i – \hat{y}_i)^2 = (\bm{y} – \bm{X}\bm{\beta})^\top(\bm{y} – \bm{X}\bm{\beta}) \end{equation} $$

行列を展開します。

$$ \text{RSS}(\bm{\beta}) = \bm{y}^\top\bm{y} – 2\bm{\beta}^\top\bm{X}^\top\bm{y} + \bm{\beta}^\top\bm{X}^\top\bm{X}\bm{\beta} $$

正規方程式の導出

$\bm{\beta}$ で微分してゼロとおきます。ベクトルの微分公式を使うと次のようになります。

$$ \frac{\partial \text{RSS}}{\partial \bm{\beta}} = -2\bm{X}^\top\bm{y} + 2\bm{X}^\top\bm{X}\bm{\beta} = \bm{0} $$

整理すると正規方程式（normal equations）が得られます。

$$ \begin{equation} \bm{X}^\top\bm{X}\hat{\bm{\beta}} = \bm{X}^\top\bm{y} \end{equation} $$

$\bm{X}^\top\bm{X}$ が正則（逆行列が存在する）であれば、最小二乗推定量が一意に求まります。

$$ \begin{equation} \hat{\bm{\beta}} = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\bm{y} \end{equation} $$

この式は「観測データから回帰係数を推定する」ための基本公式であり、統計学における最も重要な結果の一つです。

予測値とハット行列

予測値 $\hat{\bm{y}}$ は次のように書けます。

$$ \hat{\bm{y}} = \bm{X}\hat{\bm{\beta}} = \bm{X}(\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\bm{y} = \bm{H}\bm{y} $$

ここで $\bm{H} = \bm{X}(\bm{X}^\top\bm{X})^{-1}\bm{X}^\top$ をハット行列（hat matrix, projection matrix）と呼びます。$\bm{H}$ は $\bm{y}$ を $\hat{\bm{y}}$ に「帽子をかぶせる」（hat をつける）操作を行う射影行列で、$\bm{H}^2 = \bm{H}$（冪等性）、$\bm{H}^\top = \bm{H}$（対称性）を満たします。

幾何学的には、$\hat{\bm{y}}$ は $n$ 次元空間における $\bm{y}$ の、$\bm{X}$ の列ベクトルが張る部分空間への直交射影です。残差ベクトル $\hat{\bm{\epsilon}} = \bm{y} – \hat{\bm{y}} = (\bm{I} – \bm{H})\bm{y}$ はこの部分空間に直交しています。

推定量の最適性を保証する定理を見ていきましょう。

ガウス=マルコフの定理

BLUEとは

ガウス=マルコフの定理（Gauss-Markov theorem）は、線形モデルの仮定（$E[\bm{\epsilon}] = \bm{0}$、$\text{Cov}(\bm{\epsilon}) = \sigma^2\bm{I}$）のもとで、最小二乗推定量 $\hat{\bm{\beta}}$ がBLUE（Best Linear Unbiased Estimator）であることを主張します。

  • Linear（線形）: $\hat{\bm{\beta}}$ は $\bm{y}$ の線形関数（$\hat{\bm{\beta}} = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\bm{y}$）
  • Unbiased（不偏）: $E[\hat{\bm{\beta}}] = \bm{\beta}$
  • Best（最良）: すべての線形不偏推定量の中で分散が最小

つまり、$\bm{\beta}$ の任意の線形結合 $\bm{c}^\top\bm{\beta}$ に対して、$\bm{c}^\top\hat{\bm{\beta}}$ は $\bm{c}^\top\bm{\beta}$ の線形不偏推定量の中で分散が最小です。

不偏性の証明

$\hat{\bm{\beta}} = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\bm{y}$ に $\bm{y} = \bm{X}\bm{\beta} + \bm{\epsilon}$ を代入します。

$$ \hat{\bm{\beta}} = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top(\bm{X}\bm{\beta} + \bm{\epsilon}) = \bm{\beta} + (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\bm{\epsilon} $$

期待値を取ると次のようになります。

$$ E[\hat{\bm{\beta}}] = \bm{\beta} + (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top E[\bm{\epsilon}] = \bm{\beta} $$

$E[\bm{\epsilon}] = \bm{0}$ を使いました。したがって $\hat{\bm{\beta}}$ は不偏です。

推定量の分散

$\hat{\bm{\beta}} – \bm{\beta} = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\bm{\epsilon}$ より、共分散行列は次のようになります。

$$ \begin{equation} \text{Cov}(\hat{\bm{\beta}}) = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\text{Cov}(\bm{\epsilon})\bm{X}(\bm{X}^\top\bm{X})^{-1} = \sigma^2(\bm{X}^\top\bm{X})^{-1} \end{equation} $$

$\hat{\beta}_j$ の分散は $\sigma^2[(\bm{X}^\top\bm{X})^{-1}]_{jj}$ です。$\sigma^2$ の不偏推定量は $\hat{\sigma}^2 = \text{RSS}/(n – p – 1)$ です。

次に、モデルの全体的な適合度を測る指標を見ていきましょう。

決定係数

$R^2$ の定義

決定係数（coefficient of determination）$R^2$ は、モデルがデータの変動をどれだけ説明しているかを示す指標です。

全変動を次の3つに分解します。

$$ \begin{equation} \underbrace{\sum_{i=1}^n(y_i – \bar{y})^2}_{\text{SST}} = \underbrace{\sum_{i=1}^n(\hat{y}_i – \bar{y})^2}_{\text{SSR}} + \underbrace{\sum_{i=1}^n(y_i – \hat{y}_i)^2}_{\text{RSS}} \end{equation} $$

SST（Total Sum of Squares）は全変動、SSR（Regression Sum of Squares）は回帰による変動、RSS は残差変動です。

$$ \begin{equation} R^2 = \frac{\text{SSR}}{\text{SST}} = 1 – \frac{\text{RSS}}{\text{SST}} \end{equation} $$

$R^2$ は0以上1以下の値を取り、1に近いほどモデルの当てはまりがよいことを示します。$R^2 = 1$ は完全な適合（すべてのデータ点が回帰面上にある）、$R^2 = 0$ は説明変数が目的変数を全く説明していないことを意味します。

自由度調整済み $R^2$

$R^2$ には重大な欠点があります。説明変数を増やすと、たとえ無関係な変数を追加しても $R^2$ は決して減少しません。これは過学習のリスクを隠してしまいます。

自由度調整済み $R^2$（Adjusted $R^2$）はこの問題を修正します。

$$ \begin{equation} \bar{R}^2 = 1 – \frac{\text{RSS}/(n-p-1)}{\text{SST}/(n-1)} = 1 – \frac{n-1}{n-p-1}(1 – R^2) \end{equation} $$

分子と分母をそれぞれの自由度で割ることで、パラメータ数 $p$ の増加によるペナルティを導入しています。無関係な変数を追加すると $\bar{R}^2$ は減少し得るため、モデル選択の指標としてより適切です。

偏回帰係数と偏相関の関係

偏回帰係数 $\beta_j$ と偏相関係数 $r_{y \cdot x_j | \text{rest}}$ の間には密接な関係があります。偏回帰係数は「$x_j$ が1単位変わったときの $y$ の変化量」であり、偏相関係数は「他の変数を制御した後の $x_j$ と $y$ の線形関係の強さ」です。

両者の関係は次のように書けます。

$$ \begin{equation} t_j = \frac{\hat{\beta}_j}{\text{SE}(\hat{\beta}_j)} = r_{y \cdot x_j | \text{rest}} \sqrt{\frac{n – p – 1}{1 – r_{y \cdot x_j | \text{rest}}^2}} \end{equation} $$

つまり、偏回帰係数の $t$ 検定統計量と偏相関係数の $t$ 検定統計量は同じ値を取ります。「$\beta_j = 0$ かどうか」と「$x_j$ と $y$ の偏相関がゼロかどうか」は同じ問いなのです。

F検定とモデル全体の有意性

F検定は「すべての説明変数の効果が同時にゼロ」という帰無仮説 $H_0: \beta_1 = \beta_2 = \dots = \beta_p = 0$ を検定します。

$$ \begin{equation} F = \frac{\text{SSR}/p}{\text{RSS}/(n-p-1)} = \frac{R^2/p}{(1-R^2)/(n-p-1)} \end{equation} $$

この統計量は $H_0$ のもとで $F(p, n-p-1)$ 分布に従います。$R^2$ が大きいほど $F$ 値は大きくなり、帰無仮説を棄却する根拠が強くなります。

F検定は個別の $t$ 検定とは異なる問いに答えます。$t$ 検定は「この変数は必要か？」を変数ごとに検定しますが、F検定は「説明変数セット全体が目的変数の説明に役立っているか？」を問います。多重検定の問題を避けるため、まずF検定でモデル全体の有意性を確認し、次に個別の $t$ 検定で各変数の効果を調べるのが標準的な手順です。

Pythonでこれらをまとめて実装して確認しましょう。

Pythonによる実装と可視化

重回帰分析のスクラッチ実装

合成データに対して重回帰分析を実装し、推定量の性質を確認します。

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)

# --- 合成データの生成 ---
n = 200
p = 3  # 説明変数の数

# 真の回帰係数
beta_true = np.array([5.0, 2.0, -1.5, 0.8])  # [切片, β1, β2, β3]
sigma_true = 1.5

# 説明変数（互いにやや相関がある）
mean_x = [0, 0, 0]
cov_x = [[1.0, 0.3, 0.1],
          [0.3, 1.0, 0.2],
          [0.1, 0.2, 1.0]]
X_raw = np.random.multivariate_normal(mean_x, cov_x, n)
X = np.column_stack([np.ones(n), X_raw])  # 切片を追加

# 目的変数
y = X @ beta_true + sigma_true * np.random.randn(n)

# --- 最小二乗推定 ---
beta_hat = np.linalg.inv(X.T @ X) @ X.T @ y
y_hat = X @ beta_hat
residuals = y - y_hat

# --- 統計量の計算 ---
n_obs, p_total = X.shape
SST = np.sum((y - y.mean())**2)
SSR = np.sum((y_hat - y.mean())**2)
RSS = np.sum(residuals**2)
R2 = 1 - RSS / SST
R2_adj = 1 - (n_obs - 1) / (n_obs - p_total) * (1 - R2)
sigma2_hat = RSS / (n_obs - p_total)
se_beta = np.sqrt(sigma2_hat * np.diag(np.linalg.inv(X.T @ X)))
t_values = beta_hat / se_beta

print("=== 重回帰分析の結果 ===")
print(f"{'':>12s} {'真値':>8s} {'推定値':>8s} {'SE':>8s} {'t値':>8s}")
names = ["切片", "β1", "β2", "β3"]
for i in range(p_total):
    print(f"{names[i]:>12s} {beta_true[i]:8.3f} {beta_hat[i]:8.3f} "
          f"{se_beta[i]:8.3f} {t_values[i]:8.3f}")
print(f"\nR² = {R2:.4f}")
print(f"Adjusted R² = {R2_adj:.4f}")
print(f"σ² (推定) = {sigma2_hat:.4f} (真値 = {sigma_true**2:.4f})")

この結果から、以下が確認できます。

  1. 推定値は真の値に近い: 各回帰係数の推定値が真の値から大きく外れていないことがわかります。$n = 200$ は十分なサンプルサイズであり、不偏推定量の良好な性質が発揮されています

  2. t値: t値の絶対値が大きい係数は統計的に有意です。一般に $|t| > 2$ であれば5%水準で有意とみなせます

  3. $\sigma^2$ の推定: 推定された誤差分散が真の値 $\sigma^2 = 2.25$ に近いことが確認できます

可視化

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(42)
n = 200; p = 3
beta_true = np.array([5.0, 2.0, -1.5, 0.8])
sigma_true = 1.5
X_raw = np.random.multivariate_normal([0,0,0], [[1,0.3,0.1],[0.3,1,0.2],[0.1,0.2,1]], n)
X = np.column_stack([np.ones(n), X_raw])
y = X @ beta_true + sigma_true * np.random.randn(n)
beta_hat = np.linalg.inv(X.T @ X) @ X.T @ y
y_hat = X @ beta_hat
residuals = y - y_hat
SST = np.sum((y - y.mean())**2)
RSS = np.sum(residuals**2)
R2 = 1 - RSS / SST

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# (a) 実測値 vs 予測値
ax = axes[0]
ax.scatter(y, y_hat, alpha=0.4, s=15, color="steelblue")
lim = [min(y.min(), y_hat.min())-1, max(y.max(), y_hat.max())+1]
ax.plot(lim, lim, "r--", linewidth=1.5)
ax.set_xlabel("Actual $y$", fontsize=11)
ax.set_ylabel("Predicted $\\hat{y}$", fontsize=11)
ax.set_title(f"(a) Actual vs Predicted ($R^2={R2:.3f}$)", fontsize=12)
ax.grid(True, alpha=0.3)

# (b) 残差プロット
ax = axes[1]
ax.scatter(y_hat, residuals, alpha=0.4, s=15, color="gray")
ax.axhline(0, color="red", linewidth=1, linestyle="--")
ax.set_xlabel("Predicted $\\hat{y}$", fontsize=11)
ax.set_ylabel("Residual $e$", fontsize=11)
ax.set_title("(b) Residual plot", fontsize=12)
ax.grid(True, alpha=0.3)

# (c) 残差のQ-Qプロット
ax = axes[2]
from scipy import stats
res_sorted = np.sort(residuals)
theoretical_quantiles = stats.norm.ppf(np.arange(1, n+1) / (n+1))
ax.scatter(theoretical_quantiles, res_sorted, alpha=0.4, s=15, color="steelblue")
ax.plot([-3, 3], [-3*np.std(residuals), 3*np.std(residuals)], "r--", linewidth=1.5)
ax.set_xlabel("Theoretical quantiles", fontsize=11)
ax.set_ylabel("Sample quantiles", fontsize=11)
ax.set_title("(c) Q-Q plot of residuals", fontsize=12)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig("multiple_regression.png", dpi=150, bbox_inches="tight")
plt.show()

この可視化から、モデルの適合度と仮定の妥当性を診断できます。

  1. 実測値 vs 予測値（図a）: データ点が45度線の周辺に集中しており、モデルが良好な予測精度を持つことがわかります。散らばりの大きさは残差の標準偏差 $\hat{\sigma}$ に対応しています

  2. 残差プロット（図b）: 残差が予測値に対してランダムに散らばり、特定のパターン（曲線や漏斗形）が見られません。これは線形性と等分散性の仮定が妥当であることを示唆しています

  3. Q-Qプロット（図c）: 残差の分布が正規分布に近いかを視覚的に確認できます。データ点が直線上に並んでいれば正規性の仮定が支持されます

まとめ

本記事では、重回帰分析の理論を体系的に解説しました。

  • モデル: $\bm{y} = \bm{X}\bm{\beta} + \bm{\epsilon}$ で、$p$ 個の説明変数から目的変数を予測します
  • 最小二乗推定量: $\hat{\bm{\beta}} = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\bm{y}$ は正規方程式の解として得られます
  • ガウス=マルコフの定理: 標準的な仮定のもとで、最小二乗推定量はBLUE（最良線形不偏推定量）です
  • 決定係数: $R^2 = 1 – \text{RSS}/\text{SST}$ はモデルの説明力を測り、自由度調整済み $\bar{R}^2$ は変数数のペナルティを加えます
  • 幾何学的解釈: 予測値はデータベクトル $\bm{y}$ の列空間への直交射影であり、残差は列空間に直交します

重回帰分析は多変量統計手法の基礎であり、正則化回帰や一般化線形モデルなど、より高度な手法の出発点です。

次のステップとして、以下の記事も参考にしてください。

関連記事

前の記事

ロジスティック回帰の理論と最尤推定

次の記事

ホテリングのT²検定 — 多変量の平均検定

関連タグ: 多変量解析 最小二乗法 正規方程式 決定係数 統計学 重回帰分析