ある工場で製品の品質を3つの指標（強度・耐久性・外観スコア）で評価しているとします。これら3つの品質指標を、製造条件（温度・圧力・時間・材料配合比など）から同時に予測したい場合、各指標について個別に回帰分析を行うこともできますが、3つの指標は互いに相関しているかもしれません。相関構造を無視して個別に分析するのは情報のロスにつながります。
多変量回帰分析（Multivariate Regression）は、複数の説明変数から複数の目的変数を同時に予測するための手法です。目的変数が1つの通常の回帰分析（重回帰分析）の自然な拡張であり、目的変数間の相関構造も考慮できます。
多変量回帰分析を理解すると、以下のような応用が開けます。
- 品質工学: 複数の品質特性を製造条件から同時に予測・最適化する
- 経済学: GDP成長率・失業率・インフレ率を同時にモデル化する（VAR モデルの基礎）
- 生物学: 遺伝子発現データから複数の表現型を同時に予測する
- 心理学: 複数のテストスコアを背景変数で説明する
本記事では、多変量回帰モデルを行列表現で定式化し、最小二乗推定量を導出します。推定量の統計的性質と仮説検定についても解説し、Pythonで実装して動作を確認します。
本記事の内容
- 多変量回帰モデルの行列表現
- 最小二乗推定量の導出（行列微分を用いた厳密な導出）
- 推定量の統計的性質（不偏性、分散）
- 各目的変数の独立推定との関係
- Pythonでの実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
多変量回帰モデルの定式化
通常の回帰との違い
重回帰分析では、$n$ 個のデータに対して1つの目的変数 $\bm{y} \in \mathbb{R}^n$ を $p$ 個の説明変数で予測しました。多変量回帰では、目的変数が$q$ 個に増えます。
たとえば、3つの説明変数（温度 $x_1$、圧力 $x_2$、時間 $x_3$）から、2つの品質指標（強度 $y_1$、耐久性 $y_2$）を同時に予測する状況を考えてください。目的変数ごとに別々の回帰式を立てることもできますが、多変量回帰では1つのモデルとして定式化し、目的変数間の関係も同時に扱います。
行列表現
$n$ 個のデータに対して、$p$ 個の説明変数と $q$ 個の目的変数を持つ多変量回帰モデルは次のように書けます。
$$ \begin{equation} \bm{Y} = \bm{X}\bm{B} + \bm{E} \end{equation} $$
各行列のサイズと意味は以下の通りです。
- $\bm{Y} \in \mathbb{R}^{n \times q}$: 目的変数行列。第 $i$ 行は $i$ 番目のデータの $q$ 個の目的変数
- $\bm{X} \in \mathbb{R}^{n \times p}$: 計画行列（説明変数行列）。通常、第1列は切片のために $\bm{1}$ ベクトル
- $\bm{B} \in \mathbb{R}^{p \times q}$: 回帰係数行列。$B_{jk}$ は第 $j$ 説明変数が第 $k$ 目的変数に与える影響
- $\bm{E} \in \mathbb{R}^{n \times q}$: 誤差行列。第 $i$ 行 $\bm{\epsilon}_i^\top$ は $i$ 番目のデータの $q$ 個の誤差
成分で書くと、第 $i$ データの第 $k$ 目的変数は次のようになります。
$$ y_{ik} = \sum_{j=1}^{p}x_{ij}B_{jk} + \epsilon_{ik} $$
誤差項の仮定
誤差ベクトル $\bm{\epsilon}_i = (e_{i1}, e_{i2}, \dots, e_{iq})^\top$ に対して、以下を仮定します。
$$ \begin{equation} E[\bm{\epsilon}_i] = \bm{0}, \quad \text{Cov}(\bm{\epsilon}_i) = \bm{\Sigma} \end{equation} $$
$$ \begin{equation} \text{Cov}(\bm{\epsilon}_i, \bm{\epsilon}_j) = \bm{O} \quad (i \neq j) \end{equation} $$
ここで $\bm{\Sigma}$ は $q \times q$ の正定値行列で、目的変数間の誤差の共分散行列です。異なるデータ点の誤差は無相関ですが、同じデータ点の異なる目的変数の誤差は相関していてもかまいません。
この共分散構造を行列で表すと、$\text{vec}(\bm{E})$（$\bm{E}$ を列方向にベクトル化したもの）の共分散行列は $\bm{\Sigma} \otimes \bm{I}_n$ です（$\otimes$ はクロネッカー積）。
多変量回帰モデルの構造を理解したところで、パラメータ $\bm{B}$ の推定に進みましょう。
最小二乗推定量の導出
残差二乗和行列
最小二乗法の目標は、残差 $\bm{E} = \bm{Y} – \bm{X}\bm{B}$ をある意味で最小化することです。多変量の場合、残差二乗和はスカラーではなく行列になります。
残差二乗和行列を次のように定義します。
$$ \begin{equation} \bm{Q}(\bm{B}) = (\bm{Y} – \bm{X}\bm{B})^\top(\bm{Y} – \bm{X}\bm{B}) \end{equation} $$
$\bm{Q}(\bm{B})$ は $q \times q$ の行列で、その $(k, l)$ 成分は第 $k$ 目的変数と第 $l$ 目的変数の残差の内積です。対角成分 $Q_{kk}$ は第 $k$ 目的変数の残差二乗和です。
スカラーの目的関数としては、$\bm{Q}(\bm{B})$ のトレース（対角成分の和）を最小化するか、行列式を最小化するかの選択肢がありますが、最も基本的なのはトレースの最小化です。
$$ \min_{\bm{B}} \quad \text{tr}(\bm{Q}(\bm{B})) = \text{tr}[(\bm{Y} – \bm{X}\bm{B})^\top(\bm{Y} – \bm{X}\bm{B})] $$
トレースの最小化は、すべての目的変数の残差二乗和の合計を最小化することに対応します。
行列微分による導出
$\bm{Q}(\bm{B})$ のトレースを $\bm{B}$ で微分するために、まず展開します。
$$ \text{tr}(\bm{Q}) = \text{tr}[(\bm{Y} – \bm{X}\bm{B})^\top(\bm{Y} – \bm{X}\bm{B})] $$
展開すると次のようになります。
$$ \text{tr}(\bm{Q}) = \text{tr}(\bm{Y}^\top\bm{Y}) – 2\text{tr}(\bm{B}^\top\bm{X}^\top\bm{Y}) + \text{tr}(\bm{B}^\top\bm{X}^\top\bm{X}\bm{B}) $$
ここでトレースの性質 $\text{tr}(\bm{A}^\top\bm{B}) = \text{tr}(\bm{B}^\top\bm{A})$ を使いました。
$\bm{B}$ に関する行列微分の公式を適用します。
$$ \frac{\partial}{\partial \bm{B}}\text{tr}(\bm{B}^\top\bm{X}^\top\bm{Y}) = \bm{X}^\top\bm{Y} $$
$$ \frac{\partial}{\partial \bm{B}}\text{tr}(\bm{B}^\top\bm{X}^\top\bm{X}\bm{B}) = 2\bm{X}^\top\bm{X}\bm{B} $$
微分をゼロとおくと次の方程式が得られます。
$$ -2\bm{X}^\top\bm{Y} + 2\bm{X}^\top\bm{X}\hat{\bm{B}} = \bm{O} $$
これを $\hat{\bm{B}}$ について解くと、多変量最小二乗推定量が得られます。
$$ \begin{equation} \hat{\bm{B}} = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\bm{Y} \end{equation} $$
この結果は重回帰の場合 $\hat{\bm{\beta}} = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\bm{y}$ の直接的な拡張です。$\bm{y}$ が行列 $\bm{Y}$ に変わっただけで、形は全く同じです。
重要な含意 — 各目的変数の独立推定との同値性
推定量 $\hat{\bm{B}} = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\bm{Y}$ の第 $k$ 列は次のようになります。
$$ \hat{\bm{B}}_{:,k} = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\bm{Y}_{:,k} = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\bm{y}_k $$
これは、第 $k$ 目的変数 $\bm{y}_k$ だけに対して個別に重回帰分析を行った場合の推定量と完全に一致します。
つまり、多変量最小二乗推定量は、各目的変数について独立に最小二乗法を適用した結果を列として並べたものです。目的変数間の誤差の相関（共分散行列 $\bm{\Sigma}$）は、最小二乗推定量には影響を与えません。
これは一見驚くべき結果です。しかし、すべての目的変数が同じ説明変数 $\bm{X}$ を共有している場合、この結果は自然です。もし目的変数ごとに異なる説明変数を使いたい場合は、一見無関連回帰（SUR: Seemingly Unrelated Regressions）モデルを用い、一般化最小二乗法（GLS）で推定する必要があります。
ただし、$\bm{\Sigma}$ は推定量の精度（分散）や仮説検定には影響するため、$\bm{\Sigma}$ の推定も重要です。次にこの点を見ていきましょう。
推定量の統計的性質
不偏性
$\hat{\bm{B}}$ に真のモデル $\bm{Y} = \bm{X}\bm{B} + \bm{E}$ を代入すると次のようになります。
$$ \hat{\bm{B}} = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top(\bm{X}\bm{B} + \bm{E}) = \bm{B} + (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\bm{E} $$
期待値を取ると $E[\hat{\bm{B}}] = \bm{B} + (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top E[\bm{E}] = \bm{B}$ です。したがって $\hat{\bm{B}}$ は $\bm{B}$ の不偏推定量です。
誤差共分散行列の推定
誤差共分散行列 $\bm{\Sigma}$ の不偏推定量は、残差から次のように計算されます。
$$ \begin{equation} \hat{\bm{\Sigma}} = \frac{1}{n – p}(\bm{Y} – \bm{X}\hat{\bm{B}})^\top(\bm{Y} – \bm{X}\hat{\bm{B}}) = \frac{1}{n – p}\hat{\bm{E}}^\top\hat{\bm{E}} \end{equation} $$
分母の $n – p$ は自由度補正です（$p$ 個のパラメータを推定したため）。
回帰係数の分散
$\hat{\bm{B}}$ の第 $(j, k)$ 成分の分散は次のように書けます。
$$ \text{Var}(\hat{B}_{jk}) = \sigma_{kk}[(\bm{X}^\top\bm{X})^{-1}]_{jj} $$
ここで $\sigma_{kk}$ は $\bm{\Sigma}$ の $(k, k)$ 成分、$[(\bm{X}^\top\bm{X})^{-1}]_{jj}$ は $(\bm{X}^\top\bm{X})^{-1}$ の $(j, j)$ 成分です。
この結果も、各目的変数について個別に回帰分析を行った場合と同じです。
多変量仮説検定
多変量回帰では、「すべての目的変数に対して $x_j$ の効果がゼロ」という仮説を同時に検定することがあります。つまり $H_0: \bm{B}_{j,:} = \bm{0}$（$j$ 行目がすべてゼロ）を検定します。
この検定には以下の4つの多変量検定統計量が使われます。
ウィルクスのラムダ: 最も広く使われる統計量です。
$$ \begin{equation} \Lambda = \frac{|\bm{E}|}{|\bm{E} + \bm{H}|} \end{equation} $$
ここで $\bm{E}$ は残差平方和行列、$\bm{H}$ は仮説平方和行列です。$\Lambda$ は0から1の値を取り、小さいほど帰無仮説を棄却する証拠が強くなります。
ピライのトレース: $\text{tr}(\bm{H}(\bm{E} + \bm{H})^{-1})$
ホテリング=ローリーのトレース: $\text{tr}(\bm{H}\bm{E}^{-1})$
ロイの最大根: $\bm{H}\bm{E}^{-1}$ の最大固有値
これらの統計量は一般的な場合には異なる値を与えますが、仮説の自由度が1の場合にはすべて同じ結果になります。
SUR（一見無関連回帰）モデルとの関係
前節で、同じ $\bm{X}$ を共有する場合に多変量最小二乗推定量が各目的変数の独立推定と一致することを示しました。しかし、目的変数ごとに異なる説明変数を使いたい場合はどうでしょうか。
$$ \bm{y}_k = \bm{X}_k\bm{\beta}_k + \bm{\epsilon}_k, \quad k = 1, 2, \dots, q $$
このモデルは一見無関連回帰（SUR: Seemingly Unrelated Regressions）と呼ばれ、Zelner（1962）が提案しました。各方程式は見た目上は無関連ですが、誤差項が相関している可能性があります。
SURモデルでは、誤差の共分散構造 $\bm{\Sigma}$ を利用した一般化最小二乗法（GLS）で推定することで、各方程式を独立に推定するよりも効率的な推定量（より小さい分散）が得られます。ただし、すべての方程式が同じ説明変数を共有する場合は、SURのGLS推定量はOLS推定量に一致し、効率の改善はありません。
理論をPythonで実装して確認しましょう。
Pythonによる実装と可視化
多変量回帰の実装
合成データを使って多変量回帰を実装し、推定量の性質を確認します。
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# --- パラメータ設定 ---
n = 200 # サンプル数
p = 3 # 説明変数の数（切片含む）
q = 2 # 目的変数の数
# 真の回帰係数行列
B_true = np.array([
[1.5, -0.5], # 切片
[2.0, 1.0], # x1の係数
[-1.0, 3.0], # x2の係数
])
# 誤差の共分散行列（目的変数間に相関あり）
Sigma_true = np.array([[1.0, 0.6],
[0.6, 1.5]])
# --- データ生成 ---
X_raw = np.random.randn(n, p - 1) # 説明変数
X = np.column_stack([np.ones(n), X_raw]) # 切片を追加
E = np.random.multivariate_normal(np.zeros(q), Sigma_true, n)
Y = X @ B_true + E
# --- 最小二乗推定 ---
B_hat = np.linalg.inv(X.T @ X) @ X.T @ Y
# 残差と誤差共分散行列の推定
E_hat = Y - X @ B_hat
Sigma_hat = E_hat.T @ E_hat / (n - p)
# --- 各目的変数の独立推定（一致の確認） ---
B_indep = np.zeros((p, q))
for k in range(q):
B_indep[:, k] = np.linalg.inv(X.T @ X) @ X.T @ Y[:, k]
print("=== 多変量回帰の結果 ===")
print(f"\n真の回帰係数行列 B:")
print(np.array2string(B_true, precision=3))
print(f"\n推定値 B_hat (多変量):")
print(np.array2string(B_hat, precision=3))
print(f"\n推定値 B_hat (各目的変数独立):")
print(np.array2string(B_indep, precision=3))
print(f"\n差の最大値: {np.max(np.abs(B_hat - B_indep)):.2e}")
print(f"\n真の誤差共分散行列:")
print(np.array2string(Sigma_true, precision=3))
print(f"\n推定された誤差共分散行列:")
print(np.array2string(Sigma_hat, precision=3))
この結果から、以下の重要な点が確認できます。
-
推定値は真の値に近い: $n = 200$ のサンプルで、回帰係数の推定値が真の値によく一致しています。不偏推定量なので、十分なサンプルサイズがあれば真の値に収束します
-
多変量推定と独立推定が一致: 多変量回帰の推定量と、各目的変数について個別に推定した結果が数値誤差の範囲内で完全に一致しています。これは理論的に導いた結果を数値的に確認するものです
-
誤差共分散の推定: 推定された $\hat{\bm{\Sigma}}$ が真の $\bm{\Sigma}$ に近い値を示しています。目的変数間の正の相関（$\sigma_{12} = 0.6$）も適切に推定されています
予測と決定係数の可視化
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
n = 200; p = 3; q = 2
B_true = np.array([[1.5, -0.5], [2.0, 1.0], [-1.0, 3.0]])
Sigma_true = np.array([[1.0, 0.6], [0.6, 1.5]])
X_raw = np.random.randn(n, p-1)
X = np.column_stack([np.ones(n), X_raw])
E = np.random.multivariate_normal(np.zeros(q), Sigma_true, n)
Y = X @ B_true + E
B_hat = np.linalg.inv(X.T @ X) @ X.T @ Y
Y_hat = X @ B_hat
E_hat = Y - Y_hat
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# (a)(b) 実測値 vs 予測値
for k in range(q):
ax = axes[k]
ax.scatter(Y[:, k], Y_hat[:, k], alpha=0.4, s=15, color="steelblue")
lim = [min(Y[:, k].min(), Y_hat[:, k].min()) - 0.5,
max(Y[:, k].max(), Y_hat[:, k].max()) + 0.5]
ax.plot(lim, lim, "r--", linewidth=1.5, label="y = x")
SS_res = np.sum((Y[:, k] - Y_hat[:, k])**2)
SS_tot = np.sum((Y[:, k] - Y[:, k].mean())**2)
R2 = 1 - SS_res / SS_tot
ax.set_xlabel(f"Actual $y_{k+1}$", fontsize=11)
ax.set_ylabel(f"Predicted $\\hat{{y}}_{k+1}$", fontsize=11)
ax.set_title(f"({chr(97+k)}) $y_{k+1}$: $R^2={R2:.3f}$", fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_aspect("equal")
# (c) 残差の散布図（目的変数間の相関を確認）
ax = axes[2]
ax.scatter(E_hat[:, 0], E_hat[:, 1], alpha=0.4, s=15, color="gray")
corr_resid = np.corrcoef(E_hat[:, 0], E_hat[:, 1])[0, 1]
ax.set_xlabel("Residual of $y_1$", fontsize=11)
ax.set_ylabel("Residual of $y_2$", fontsize=11)
ax.set_title(f"(c) Residual correlation (r={corr_resid:.3f})", fontsize=12)
ax.grid(True, alpha=0.3)
ax.axhline(0, color="k", linewidth=0.5)
ax.axvline(0, color="k", linewidth=0.5)
plt.tight_layout()
plt.savefig("multivariate_regression.png", dpi=150, bbox_inches="tight")
plt.show()
この可視化から、以下の知見が得られます。
-
実測値 vs 予測値（図a, b）: 両方の目的変数について、予測値が実測値の周辺に集中しており、モデルがデータをよく説明していることがわかります。$R^2$ 値はそれぞれの目的変数に対する説明力を示しています
-
残差の相関（図c）: 残差の散布図に正の相関が見られます。相関係数は真の値 $\sigma_{12}/\sqrt{\sigma_{11}\sigma_{22}} = 0.6/\sqrt{1.0 \times 1.5} \approx 0.49$ に近い値です。この残差間相関は、多変量回帰の最小二乗推定量自体には影響しませんが、仮説検定やモデル比較では重要な役割を果たします
まとめ
本記事では、多変量回帰分析の理論を解説しました。
- モデル: $\bm{Y} = \bm{X}\bm{B} + \bm{E}$ で、$q$ 個の目的変数を同時にモデル化します
- 最小二乗推定量: $\hat{\bm{B}} = (\bm{X}^\top\bm{X})^{-1}\bm{X}^\top\bm{Y}$ は重回帰の公式と同じ形です
- 独立推定との同値性: 同じ説明変数を共有する場合、多変量推定は各目的変数の個別推定と一致します
- 誤差共分散: $\hat{\bm{\Sigma}} = \hat{\bm{E}}^\top\hat{\bm{E}}/(n-p)$ で目的変数間の誤差の相関構造を推定できます
多変量回帰は重回帰分析の自然な拡張であり、構造方程式モデリングやVARモデルなど、より高度な多変量手法の基礎となります。
次のステップとして、以下の記事も参考にしてください。