3種類の肥料（A, B, C）を使って植物を育て、その効果を「茎の長さ」と「葉の数」の2つの指標で同時に評価したいとします。直感的には、各指標ごとにANOVA（分散分析）を行えばよさそうです。しかし、このアプローチには2つの問題があります。
1つ目は多重検定の問題です。2つの検定を独立に行うと、少なくとも一方で第1種の過誤が生じる確率は $1 – (1-0.05)^2 = 0.0975$ となり、名目水準の5%を大幅に超えてしまいます。指標が10個あれば、その確率は $1 – 0.95^{10} \approx 0.40$ にもなります。
2つ目は変数間の相関の無視です。「茎の長さ」と「葉の数」は互いに相関しているかもしれません。各変数を個別に検定すると、この相関構造が持つ情報が失われます。場合によっては、各変数単独では群間差が有意でなくても、2つの変数を組み合わせると有意な差が検出できることがあります。
多変量分散分析（MANOVA: Multivariate Analysis of Variance）は、$K$ 個のグループ間で複数の目的変数の平均ベクトルに差があるかを同時に検定する手法です。ホテリングのT²検定が2群の比較であったのに対し、MANOVAは3群以上に一般化した手法です。
MANOVAを理解すると、以下のような幅広い応用が開けます。
- 農学・生態学: 複数の生育指標（収量・茎長・根長・葉面積）に対する処理効果の同時検定
- 心理学: 複数の心理尺度（不安・抑うつ・ストレス）に対する介入効果の評価
- 医学: 複数の臨床指標（血圧・コレステロール・血糖値）に対する治療効果の比較
- 製造業: 複数の品質指標（強度・硬度・表面粗さ）に対する工程条件の効果の評価
- 教育学: 複数のテストスコア（読解・数学・科学）に対する教授法の効果の比較
本記事では、1元配置ANOVAからMANOVAへの拡張を丁寧に導出し、クラス間散布行列とクラス内散布行列の概念を明確にします。ウィルクスのラムダをはじめとする4つの検定統計量の理論と性質を比較し、帰無分布のF近似を解説します。さらにPythonで1元配置MANOVAを実装し、結果を可視化します。
本記事の内容
- 1元配置ANOVAの行列表現とMANOVAへの拡張
- クラス間散布行列 $\bm{B}$ とクラス内散布行列 $\bm{W}$ の定義と意味
- ウィルクスのラムダの導出と帰無分布
- ピライのトレース・ローリー=ホテリングのトレース・ロイの最大根
- 4つの検定統計量の比較 — どの統計量をいつ使うべきか
- MANOVAの前提条件と注意点
- Pythonでの実装と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- ホテリングのT²検定 — 2群の多変量平均の検定。MANOVAの特殊ケースです
- 多変量正規分布 — MANOVAの確率モデルの基盤
- 判別分析 — クラス間・クラス内散布行列の定義と幾何学的意味
ANOVAの行列表現とMANOVAへの拡張
1元配置ANOVAの復習
まず1変数の場合を復習し、それを多変数に拡張するというアプローチで進めます。
1変数の1元配置ANOVA（分散分析）では、$K$ 個のグループから得られたデータの全変動を群間変動（SSB: Sum of Squares Between）と群内変動（SSW: Sum of Squares Within）に分解します。
$$ \begin{equation} \underbrace{\sum_{k=1}^{K}\sum_{i=1}^{n_k}(x_{ki} – \bar{x})^2}_{\text{SST}} = \underbrace{\sum_{k=1}^{K}n_k(\bar{x}_k – \bar{x})^2}_{\text{SSB}} + \underbrace{\sum_{k=1}^{K}\sum_{i=1}^{n_k}(x_{ki} – \bar{x}_k)^2}_{\text{SSW}} \end{equation} $$
ここで $x_{ki}$ は第 $k$ グループの第 $i$ 番目のデータ、$\bar{x}_k$ は第 $k$ グループの平均、$\bar{x}$ は全体平均、$n_k$ は第 $k$ グループのサンプル数です。
帰無仮説 $H_0: \mu_1 = \mu_2 = \dots = \mu_K$（全グループの母平均が等しい）のもとでは、SSBは「偶然による変動」のみを反映し、対立仮説のもとではSSBは「処理効果による変動」を追加的に含みます。F統計量:
$$ F = \frac{SSB/(K-1)}{SSW/(n-K)} $$
が帰無仮説のもとで $F(K-1, n-K)$ 分布に従うことを利用して検定を行います。
スカラーから行列への拡張
MANOVAでは、各データ点がスカラー $x_{ki}$ ではなく $p$ 次元のベクトル $\bm{x}_{ki} \in \mathbb{R}^p$ になります。ANOVAの「二乗和」が「行列」に置き換わります。
ANOVAでの二乗和の分解:
$$ SST = SSB + SSW \quad \text{（スカラーの等式）} $$
MANOVAでの散布行列の分解:
$$ \begin{equation} \bm{T} = \bm{B} + \bm{W} \quad \text{（行列の等式）} \end{equation} $$
ここで $\bm{T}, \bm{B}, \bm{W}$ はそれぞれ $p \times p$ の行列であり、各方向のスカラーの分解を全ての方向について同時にまとめたものです。
この対応関係を理解するために、各行列を正式に定義しましょう。
クラス間散布行列とクラス内散布行列
全散布行列 $\bm{T}$
全散布行列（Total scatter matrix）は、全データの全体平均からの散らばりを表す行列です。
$$ \begin{equation} \bm{T} = \sum_{k=1}^{K}\sum_{i=1}^{n_k}(\bm{x}_{ki} – \bar{\bm{x}})(\bm{x}_{ki} – \bar{\bm{x}})^T \end{equation} $$
ここで $\bar{\bm{x}} = \frac{1}{n}\sum_{k=1}^{K}\sum_{i=1}^{n_k}\bm{x}_{ki}$ は全体平均ベクトル、$n = \sum_{k=1}^{K}n_k$ は全サンプル数です。
$\bm{T}$ の $(j, l)$ 成分は $\sum_{k,i}(x_{ki,j} – \bar{x}_j)(x_{ki,l} – \bar{x}_l)$ であり、第 $j$ 変数と第 $l$ 変数の間の全体の共変動を表しています。対角成分は各変数のSST（全変動）に対応します。
クラス内散布行列 $\bm{W}$
クラス内散布行列（Within-group scatter matrix）は、各グループ内での散らばりを全グループにわたって合計した行列です。
$$ \begin{equation} \bm{W} = \sum_{k=1}^{K}\sum_{i=1}^{n_k}(\bm{x}_{ki} – \bar{\bm{x}}_k)(\bm{x}_{ki} – \bar{\bm{x}}_k)^T = \sum_{k=1}^{K}\bm{S}_k \end{equation} $$
ここで $\bar{\bm{x}}_k = \frac{1}{n_k}\sum_{i=1}^{n_k}\bm{x}_{ki}$ は第 $k$ グループの平均ベクトル、$\bm{S}_k$ は第 $k$ グループの散布行列です。
$\bm{W}$ は処理効果を含まない「ノイズ」の変動を反映しています。帰無仮説のもとでも対立仮説のもとでも、$\bm{W}/(n-K)$ はプール共分散行列の推定量です。ANOVAのSSWに対応します。
クラス間散布行列 $\bm{B}$
クラス間散布行列（Between-group scatter matrix）は、各グループの平均が全体平均からどれだけ散らばっているかを表す行列です。
$$ \begin{equation} \bm{B} = \sum_{k=1}^{K}n_k(\bar{\bm{x}}_k – \bar{\bm{x}})(\bar{\bm{x}}_k – \bar{\bm{x}})^T \end{equation} $$
$\bm{B}$ は処理効果（群間の差）を反映しています。帰無仮説のもとでは $\bm{B}/(K-1)$ と $\bm{W}/(n-K)$ はどちらも共分散行列 $\bm{\Sigma}$ の推定量になります。対立仮説のもとでは $\bm{B}$ が追加の変動を含み、$\bm{B}/(K-1)$ は $\bm{\Sigma}$ よりも「大きく」なります。
$\bm{B}$ のランクは $\min(p, K-1)$ 以下です。$K$ 個のグループ平均から全体平均を引いた $K$ 個のベクトルは、$K-1$ 次元の部分空間に含まれるため（1つは他の線形結合で表せる）、$\bm{B}$ のランクは最大でも $K-1$ です。
分解の証明
$\bm{T} = \bm{B} + \bm{W}$ を証明しましょう。任意のデータ点 $\bm{x}_{ki}$ について:
$$ \bm{x}_{ki} – \bar{\bm{x}} = (\bm{x}_{ki} – \bar{\bm{x}}_k) + (\bar{\bm{x}}_k – \bar{\bm{x}}) $$
両辺の外積をとり、全データにわたって和をとります。
$$ \bm{T} = \sum_{k,i}(\bm{x}_{ki} – \bar{\bm{x}})(\bm{x}_{ki} – \bar{\bm{x}})^T $$
展開すると:
$$ \bm{T} = \sum_{k,i}\left[(\bm{x}_{ki} – \bar{\bm{x}}_k)(\bm{x}_{ki} – \bar{\bm{x}}_k)^T + (\bm{x}_{ki} – \bar{\bm{x}}_k)(\bar{\bm{x}}_k – \bar{\bm{x}})^T + (\bar{\bm{x}}_k – \bar{\bm{x}})(\bm{x}_{ki} – \bar{\bm{x}}_k)^T + (\bar{\bm{x}}_k – \bar{\bm{x}})(\bar{\bm{x}}_k – \bar{\bm{x}})^T\right] $$
交差項について、第 $k$ グループ内で $i$ について和をとると:
$$ \sum_{i=1}^{n_k}(\bm{x}_{ki} – \bar{\bm{x}}_k) = \bm{0} $$
ですから、交差項は全てゼロになります。残りの項をまとめると:
$$ \bm{T} = \underbrace{\sum_{k,i}(\bm{x}_{ki} – \bar{\bm{x}}_k)(\bm{x}_{ki} – \bar{\bm{x}}_k)^T}_{\bm{W}} + \underbrace{\sum_{k=1}^{K}n_k(\bar{\bm{x}}_k – \bar{\bm{x}})(\bar{\bm{x}}_k – \bar{\bm{x}})^T}_{\bm{B}} $$
これで $\bm{T} = \bm{B} + \bm{W}$ が示されました。
散布行列の分解が確立されたので、次にこれらの行列から検定統計量を構成する方法を見ていきましょう。
MANOVAの帰無仮説と検定の考え方
帰無仮説
MANOVAの帰無仮説は:
$$ \begin{equation} H_0: \bm{\mu}_1 = \bm{\mu}_2 = \dots = \bm{\mu}_K \end{equation} $$
全てのグループの母平均ベクトルが等しいという仮説です。対立仮説は「少なくとも1組のグループ間で平均ベクトルが異なる」です。
行列の「比」の問題
ANOVAでは $F = (SSB/\text{df}_B)/(SSW/\text{df}_W)$ というスカラーの比で検定しました。MANOVAでは $\bm{B}$ と $\bm{W}$ が行列であるため、単純な「比」は定義できません。
行列の「相対的な大きさ」を1つのスカラーに要約する方法として、$\bm{W}^{-1}\bm{B}$ の固有値を利用する方法が採られます。$\bm{W}^{-1}\bm{B}$ の固有値 $\lambda_1 \geq \lambda_2 \geq \dots \geq \lambda_s$（$s = \min(p, K-1)$）は、$\bm{W}$ で標準化した $\bm{B}$ の「大きさ」を各方向ごとに表しています。
帰無仮説のもとでは全ての固有値がゼロに近くなり、対立仮説のもとでは少なくとも一部の固有値が大きくなります。これらの固有値をどのようにスカラーに要約するかによって、異なる検定統計量が生まれます。
それでは、4つの主要な検定統計量を見ていきましょう。
ウィルクスのラムダ
定義と導出
最も広く使われるMANOVAの検定統計量がウィルクスのラムダ（Wilks’ Lambda）です。1932年にサミュエル・ウィルクスが提案しました。
$$ \begin{equation} \Lambda = \frac{|\bm{W}|}{|\bm{W} + \bm{B}|} = \frac{|\bm{W}|}{|\bm{T}|} \end{equation} $$
ここで $|\cdot|$ は行列式を表します。
この統計量の意味を理解するために、行列式の幾何学的解釈を使いましょう。$p$ 次元空間における散布行列の行列式は、データの「超楕円体の体積」の二乗に比例します。$|\bm{W}|$ はグループ内の散らばりの体積、$|\bm{T}|$ は全データの散らばりの体積です。$\Lambda$ はこの比であり、「全体の散らばりのうちグループ内の散らばりが占める割合」を測っています。
$\Lambda$ は0から1の値を取ります。
- $\Lambda = 1$: $\bm{B} = \bm{0}$、つまり全グループの平均が等しい。群間差なし
- $\Lambda \to 0$: $\bm{B}$ が非常に大きい。群間差が顕著
$\Lambda$ が十分に小さければ帰無仮説を棄却します。
固有値による表現
$\bm{W}^{-1}\bm{B}$ の固有値 $\lambda_1, \lambda_2, \dots, \lambda_s$ を用いると、ウィルクスのラムダは次のように表されます。
$$ \begin{equation} \Lambda = \prod_{i=1}^{s}\frac{1}{1 + \lambda_i} \end{equation} $$
この表現は、ウィルクスのラムダが全ての固有値の情報を「積」として統合していることを示しています。
導出を見てみましょう。$\bm{W}^{-1}\bm{B}$ の固有値を $\lambda_i$ とすると、$\bm{W}^{-1}(\bm{W} + \bm{B}) = \bm{I} + \bm{W}^{-1}\bm{B}$ の固有値は $1 + \lambda_i$ です。したがって:
$$ \frac{|\bm{W}|}{|\bm{W} + \bm{B}|} = \frac{|\bm{W}|}{|\bm{W}||\bm{I} + \bm{W}^{-1}\bm{B}|} = \frac{1}{\prod_i(1+\lambda_i)} = \prod_i\frac{1}{1+\lambda_i} $$
ここで $|\bm{W}+\bm{B}| = |\bm{W}||\bm{W}^{-1}(\bm{W}+\bm{B})| = |\bm{W}||\bm{I}+\bm{W}^{-1}\bm{B}|$ を用いました。
帰無分布のF近似
ウィルクスのラムダの正確な帰無分布は複雑ですが、いくつかの特殊ケースでは正確なF分布に変換できます。
$K = 2$（2群）の場合: $\Lambda = 1/(1+\lambda_1)$ であり、$T^2 = (n-2)\lambda_1$ はホテリングのT²統計量に一致します。
$p = 1$（1変数）の場合: $\Lambda = SSW/SST = 1 – SSB/SST$ であり、$F = (1-\Lambda)\cdot(n-K)/(\Lambda\cdot(K-1))$ は通常のF統計量に一致します。
一般の場合には、Raoのf近似が広く用いられます。
$$ F \approx \frac{1 – \Lambda^{1/t}}{\Lambda^{1/t}} \cdot \frac{df_2}{df_1} $$
ここで:
$$ t = \begin{cases} 1 & \text{if } \min(p, K-1) = 1 \\ \sqrt{\frac{p^2(K-1)^2 – 4}{p^2 + (K-1)^2 – 5}} & \text{otherwise} \end{cases} $$
$$ df_1 = p(K-1), \quad df_2 = \left(n – 1 – \frac{p+K}{2}\right)t – \frac{p(K-1)}{2} + 1 $$
この $F$ は近似的に $F(df_1, df_2)$ 分布に従います。
ウィルクスのラムダの他にも、固有値を異なる方法で要約する検定統計量があります。次にそれらを見ていきましょう。
その他の検定統計量
ピライのトレース（Pillai’s Trace）
$$ \begin{equation} V = \text{tr}(\bm{B}(\bm{W} + \bm{B})^{-1}) = \sum_{i=1}^{s}\frac{\lambda_i}{1 + \lambda_i} \end{equation} $$
ピライのトレースは固有値を「和」として統合します。$\lambda_i/(1+\lambda_i)$ は0から1の値を取るため、$V$ は0から $s$ の値を取ります。大きいほど群間差が大きいことを示します。
特徴: 4つの統計量の中で最も頑健（共分散の等質性の仮定や正規性の仮定からのずれに対して強い）です。仮定に疑いがある場合はピライのトレースが推奨されます。
ローリー=ホテリングのトレース（Lawley-Hotelling Trace）
$$ \begin{equation} U = \text{tr}(\bm{W}^{-1}\bm{B}) = \sum_{i=1}^{s}\lambda_i \end{equation} $$
固有値の単純な和です。$K = 2$ のとき、$U$ はホテリングのT²統計量を $(n-2)$ で割ったものに一致します。
特徴: ピライのトレースとウィルクスのラムダの中間的な性質を持ちます。
ロイの最大根（Roy’s Largest Root）
$$ \begin{equation} \Theta = \frac{\lambda_1}{1 + \lambda_1} \end{equation} $$
最大固有値のみを使う統計量です。
特徴: 群間の差が1つの方向に集中している場合（例えば処理効果が1つの判別関数の方向に沿っている場合）に最も検出力が高くなります。しかし、仮定からのずれに最も敏感であり、正確な帰無分布の計算も困難です。
4つの統計量の比較
| 統計量 | 固有値の要約方法 | 頑健性 | 検出力（1方向の差） | 検出力（多方向の差） |
|---|---|---|---|---|
| ウィルクスのラムダ | 積 | 中程度 | 中程度 | 中程度 |
| ピライのトレース | 和（正規化） | 最も高い | やや低い | 高い |
| ローリー=ホテリング | 和 | 中程度 | 中程度 | 中程度 |
| ロイの最大根 | 最大値のみ | 最も低い | 最も高い | 低い |
実用上の推奨:
- デフォルト: ウィルクスのラムダが最も標準的であり、多くのソフトウェアのデフォルトです
- 仮定に疑いがある場合: ピライのトレースを使う。共分散行列の等質性やサンプルサイズのバランスに問題がある場合に特に推奨
- 差が1方向に集中している確信がある場合: ロイの最大根が最も検出力が高い
4つの統計量が全て同じ結論（棄却/非棄却）を出す場合は安心できますが、結論が異なる場合はデータの性質を慎重に検討する必要があります。
検定統計量の理論的性質を理解した上で、MANOVAの前提条件と注意点を確認しましょう。
MANOVAの前提条件と注意点
前提条件
MANOVAは以下の前提条件を必要とします。
前提1: 多変量正規性
各グループのデータが多変量正規分布に従うこと。各グループのサンプルサイズが十分に大きければ（目安: 各グループ20以上）、中心極限定理により多少の逸脱は許容されます。
前提2: 共分散行列の等質性（等分散性）
$$ \bm{\Sigma}_1 = \bm{\Sigma}_2 = \dots = \bm{\Sigma}_K $$
全グループの母共分散行列が等しいこと。BoxのM検定で検証できますが、この検定はサンプルサイズに対して過度に感度が高いため、結果の解釈には注意が必要です。共分散行列の等質性が疑わしい場合は、ピライのトレースが最も頑健な選択肢です。
前提3: 独立性
データ点が互いに独立であること。実験計画によって保証される前提です。
前提4: 線形性
変数間の関係が線形であること。MANOVAは線形の群間差を検出する手法であり、非線形の関係は検出力が低下する可能性があります。
ANOVAとの比較
| 側面 | ANOVA | MANOVA |
|---|---|---|
| 目的変数 | 1つ（スカラー） | 複数（ベクトル） |
| 検定統計量 | F統計量（スカラーの比） | $\bm{W}^{-1}\bm{B}$ の固有値に基づく |
| 多重検定 | 変数ごとに独立に検定（問題あり） | 全変数を同時に検定 |
| 変数間の相関 | 考慮しない | 考慮する |
| 検出力 | 変数間の相関情報を使わない | 相関情報を活用して検出力が向上しうる |
MANOVAが有意であった場合、どの変数がどのグループ間で差があるかを特定するために、事後分析（post-hoc analysis）を行います。典型的には、判別分析を行って判別関数を解釈するか、各変数ごとにANOVAを行い（Bonferroni補正付き）、差の構造を探ります。
それでは、Pythonで1元配置MANOVAを実装して理論を確認しましょう。
Pythonによる実装と可視化
データ生成と散布行列の計算
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
np.random.seed(42)
# --- 3グループ x 3変数のデータ生成 ---
K = 3 # グループ数
p = 3 # 変数の数
n_per_group = 50 # 各グループのサンプル数
# 各グループの母平均ベクトル
mu_list = [
np.array([2.0, 3.0, 4.0]), # Group 1
np.array([4.0, 3.5, 3.0]), # Group 2
np.array([3.0, 5.0, 5.0]), # Group 3
]
# 共通の共分散行列
Sigma = np.array([
[1.0, 0.3, 0.2],
[0.3, 0.8, 0.1],
[0.2, 0.1, 1.2]
])
# データ生成
groups = []
labels = []
for k in range(K):
X_k = np.random.multivariate_normal(mu_list[k], Sigma, n_per_group)
groups.append(X_k)
labels.extend([k] * n_per_group)
X_all = np.vstack(groups)
labels = np.array(labels)
n = len(labels)
# --- 散布行列の計算 ---
x_bar = X_all.mean(axis=0) # 全体平均
W = np.zeros((p, p))
B = np.zeros((p, p))
for k in range(K):
X_k = groups[k]
x_bar_k = X_k.mean(axis=0)
n_k = len(X_k)
# クラス内
diff_within = X_k - x_bar_k
W += diff_within.T @ diff_within
# クラス間
diff_between = (x_bar_k - x_bar).reshape(-1, 1)
B += n_k * diff_between @ diff_between.T
T = W + B
# --- 検証: T = B + W ---
print("=== 散布行列の分解の検証 ===")
print(f"||T - (B + W)||_F = {np.linalg.norm(T - B - W):.2e}")
print(f"\nクラス内散布行列 W:")
print(W.round(2))
print(f"\nクラス間散布行列 B:")
print(B.round(2))
$\|\bm{T} – (\bm{B} + \bm{W})\|$ がほぼゼロになることが確認でき、散布行列の分解 $\bm{T} = \bm{B} + \bm{W}$ が数値的にも正確であることがわかります。クラス内散布行列 $\bm{W}$ は対角成分が大きく、各変数のグループ内変動を反映しています。クラス間散布行列 $\bm{B}$ は、グループ平均の差に基づく変動を表しています。
検定統計量の計算
import numpy as np
from scipy import stats
np.random.seed(42)
K, p, n_per_group = 3, 3, 50
mu_list = [np.array([2,3,4]), np.array([4,3.5,3]), np.array([3,5,5])]
Sigma = np.array([[1,0.3,0.2],[0.3,0.8,0.1],[0.2,0.1,1.2]])
groups = [np.random.multivariate_normal(mu_list[k], Sigma, n_per_group) for k in range(K)]
X_all = np.vstack(groups)
n = len(X_all)
x_bar = X_all.mean(axis=0)
W, B = np.zeros((p,p)), np.zeros((p,p))
for k in range(K):
X_k = groups[k]
x_bar_k = X_k.mean(axis=0)
W += (X_k - x_bar_k).T @ (X_k - x_bar_k)
d = (x_bar_k - x_bar).reshape(-1,1)
B += n_per_group * d @ d.T
T = W + B
# --- 固有値の計算 ---
eigvals = np.real(np.linalg.eigvals(np.linalg.inv(W) @ B))
eigvals = np.sort(eigvals)[::-1]
s = min(p, K - 1)
eigvals_pos = eigvals[:s]
# --- 4つの検定統計量 ---
wilks_lambda = np.prod(1 / (1 + eigvals_pos))
pillai_trace = np.sum(eigvals_pos / (1 + eigvals_pos))
lh_trace = np.sum(eigvals_pos)
roy_root = eigvals_pos[0] / (1 + eigvals_pos[0])
# --- ウィルクスのラムダのF近似（Rao近似） ---
df_h = p * (K - 1)
df_e = n - K
if min(p, K-1) == 1:
t = 1.0
else:
num = p**2 * (K-1)**2 - 4
den = p**2 + (K-1)**2 - 5
t = np.sqrt(num / den) if den > 0 else 1.0
df1 = p * (K - 1)
df2 = (n - 1 - (p + K) / 2) * t - df1 / 2 + 1
F_approx = ((1 - wilks_lambda**(1/t)) / wilks_lambda**(1/t)) * (df2 / df1)
p_value = 1 - stats.f.cdf(F_approx, df1, df2)
print("=== MANOVA検定統計量 ===")
print(f"W^(-1)B の固有値: {eigvals_pos.round(4)}")
print(f"\nウィルクスのΛ = {wilks_lambda:.4f}")
print(f"ピライのV = {pillai_trace:.4f}")
print(f"ローリー=ホテリングのU = {lh_trace:.4f}")
print(f"ロイのΘ = {roy_root:.4f}")
print(f"\n--- ウィルクスのΛのF近似 ---")
print(f"F = {F_approx:.4f}")
print(f"df1 = {df1:.0f}, df2 = {df2:.1f}")
print(f"p値 = {p_value:.6f}")
print(f"結論: {'H0を棄却（群間に有意な差あり）' if p_value < 0.05 else 'H0を棄却しない'}")
検定統計量の計算結果から、4つの統計量が全て帰無仮説の棄却を支持していることが確認できます。ウィルクスのラムダは0に近い値を取り、F近似のp値は非常に小さく、3つのグループの平均ベクトルに統計的に有意な差があることを示しています。
可視化
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
np.random.seed(42)
K, p, n_per_group = 3, 3, 50
mu_list = [np.array([2,3,4]), np.array([4,3.5,3]), np.array([3,5,5])]
Sigma = np.array([[1,0.3,0.2],[0.3,0.8,0.1],[0.2,0.1,1.2]])
groups = [np.random.multivariate_normal(mu_list[k], Sigma, n_per_group) for k in range(K)]
X_all = np.vstack(groups)
n = len(X_all)
x_bar = X_all.mean(axis=0)
W, B = np.zeros((p,p)), np.zeros((p,p))
for k in range(K):
X_k = groups[k]
x_bar_k = X_k.mean(axis=0)
W += (X_k - x_bar_k).T @ (X_k - x_bar_k)
d = (x_bar_k - x_bar).reshape(-1,1)
B += n_per_group * d @ d.T
eigvals = np.sort(np.real(np.linalg.eigvals(np.linalg.inv(W) @ B)))[::-1]
s = min(p, K-1)
eigvals_pos = eigvals[:s]
eigvecs = np.real(np.linalg.eig(np.linalg.inv(W) @ B)[1])
# 判別得点（W^{-1}Bの固有ベクトルへの射影）
_, evecs = np.linalg.eig(np.linalg.inv(W) @ B)
evecs = np.real(evecs)
# 固有値の大きい順にソート
sort_idx = np.argsort(np.real(np.linalg.eigvals(np.linalg.inv(W) @ B)))[::-1]
evecs = evecs[:, sort_idx]
disc_scores = X_all @ evecs[:, :2]
fig, axes = plt.subplots(2, 2, figsize=(14, 11))
colors = ['#e41a1c', '#377eb8', '#4daf4a']
group_names = ['Group 1', 'Group 2', 'Group 3']
# (a) 散布図行列（y1 vs y2）
ax = axes[0, 0]
for k in range(K):
ax.scatter(groups[k][:, 0], groups[k][:, 1], c=colors[k],
alpha=0.5, s=25, label=group_names[k])
ax.plot(mu_list[k][0], mu_list[k][1], '+', color=colors[k],
markersize=15, markeredgewidth=2.5)
ax.plot(x_bar[0], x_bar[1], 'k*', markersize=12, label='Grand mean')
ax.set_xlabel('$y_1$', fontsize=11)
ax.set_ylabel('$y_2$', fontsize=11)
ax.set_title('(a) $y_1$ vs $y_2$', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# (b) 散布図行列（y1 vs y3）
ax = axes[0, 1]
for k in range(K):
ax.scatter(groups[k][:, 0], groups[k][:, 2], c=colors[k],
alpha=0.5, s=25, label=group_names[k])
ax.plot(mu_list[k][0], mu_list[k][2], '+', color=colors[k],
markersize=15, markeredgewidth=2.5)
ax.set_xlabel('$y_1$', fontsize=11)
ax.set_ylabel('$y_3$', fontsize=11)
ax.set_title('(b) $y_1$ vs $y_3$', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# (c) 判別得点空間
ax = axes[1, 0]
for k in range(K):
mask = np.arange(k*n_per_group, (k+1)*n_per_group)
ax.scatter(disc_scores[mask, 0], disc_scores[mask, 1], c=colors[k],
alpha=0.5, s=25, label=group_names[k])
ax.set_xlabel('Discriminant function 1', fontsize=11)
ax.set_ylabel('Discriminant function 2', fontsize=11)
ax.set_title('(c) Discriminant scores', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
# (d) 固有値
ax = axes[1, 1]
ax.bar(range(1, s+1), eigvals_pos, color='steelblue', alpha=0.7,
edgecolor='black', width=0.6)
for i in range(s):
ax.text(i+1, eigvals_pos[i]+0.2, f'{eigvals_pos[i]:.3f}',
ha='center', fontsize=11)
pct = eigvals_pos[i] / eigvals_pos.sum() * 100
ax.text(i+1, eigvals_pos[i]/2, f'{pct:.1f}%',
ha='center', fontsize=10, color='white', fontweight='bold')
ax.set_xlabel('Discriminant function', fontsize=11)
ax.set_ylabel('Eigenvalue of $W^{-1}B$', fontsize=11)
ax.set_title('(d) Eigenvalues of $W^{-1}B$', fontsize=12)
ax.set_xticks(range(1, s+1))
ax.grid(True, alpha=0.3, axis='y')
plt.tight_layout()
plt.show()
4つのプロットから、MANOVAの結果を多角的に理解できます。
-
散布図 $y_1$ vs $y_2$（図a）: 3つのグループの平均（十字マーク）が異なる位置にあることが視覚的に確認できます。グループ間の重なりはありますが、平均の位置は明確に異なっています。全体平均（黒い星）は3つのグループ平均の重心に位置しています
-
散布図 $y_1$ vs $y_3$（図b）: 別の変数の組み合わせで見ても、グループ間の分離が確認できます。各変数ペアの散布図を確認することで、どの変数の組み合わせがグループ分離に最も寄与しているかの手がかりが得られます
-
判別得点空間（図c）: $\bm{W}^{-1}\bm{B}$ の固有ベクトルで定義される判別関数への射影です。この空間ではグループ間の分離が最大化されており、元の3次元空間よりもグループの違いが明確に見えます。第1判別関数（横軸）がグループ分離の大部分を担っていることが、グラフの広がり方からわかります
-
固有値（図d）: $\bm{W}^{-1}\bm{B}$ の固有値は、各判別方向での群間分散と群内分散の比を表しています。第1固有値が大きく、群間差の多くが第1判別関数の方向に集中していることがわかります。各固有値がどの程度の割合を占めるかを棒グラフ内のパーセンテージで表示しています
個別ANOVAとの比較
MANOVAの利点を確認するために、各変数ごとの個別ANOVAの結果と比較します。
import numpy as np
from scipy import stats
np.random.seed(42)
K, p, n_per_group = 3, 3, 50
mu_list = [np.array([2,3,4]), np.array([4,3.5,3]), np.array([3,5,5])]
Sigma = np.array([[1,0.3,0.2],[0.3,0.8,0.1],[0.2,0.1,1.2]])
groups = [np.random.multivariate_normal(mu_list[k], Sigma, n_per_group) for k in range(K)]
var_names = ['$y_1$', '$y_2$', '$y_3$']
print("=== 個別ANOVA vs MANOVA ===\n")
print("--- 各変数ごとのANOVA ---")
for j in range(p):
data_by_group = [groups[k][:, j] for k in range(K)]
F_val, p_val = stats.f_oneway(*data_by_group)
sig = '***' if p_val < 0.001 else ('**' if p_val < 0.01 else ('*' if p_val < 0.05 else 'n.s.'))
print(f" {var_names[j]}: F = {F_val:7.2f}, p = {p_val:.6f} {sig}")
# Bonferroni補正
alpha_bonf = 0.05 / p
print(f"\n--- Bonferroni補正後の有意水準: {alpha_bonf:.4f} ---")
for j in range(p):
data_by_group = [groups[k][:, j] for k in range(K)]
F_val, p_val = stats.f_oneway(*data_by_group)
sig = '有意' if p_val < alpha_bonf else '非有意'
print(f" {var_names[j]}: p = {p_val:.6f} -> {sig}")
print(f"\n--- MANOVA ---")
print(f" ウィルクスのΛ: 全変数を同時に考慮して1つの検定")
print(f" 多重検定の問題なし、変数間の相関を活用")
個別ANOVAとMANOVAの比較から、MANOVAの利点が確認できます。個別ANOVAでは各変数ごとに独立にF検定を行うため、多重検定の問題が生じます。Bonferroni補正を適用すると有意水準が厳しくなり、個別には有意でなくなる場合もあります。MANOVAは全変数を同時に考慮するため、多重検定の問題を回避しつつ、変数間の相関情報を活用して検出力を維持します。
まとめ
本記事では、多変量分散分析（MANOVA）の理論を体系的に解説しました。
- MANOVAの目的: $K$ 個のグループ間で、$p$ 個の目的変数の平均ベクトルに同時に差があるかを検定する
- 散布行列の分解: $\bm{T} = \bm{B} + \bm{W}$。全変動をクラス間変動とクラス内変動に分解する。ANOVAのスカラー分解の行列版
- 検定統計量: $\bm{W}^{-1}\bm{B}$ の固有値 $\lambda_1, \dots, \lambda_s$ を異なる方法で要約して構成される
- ウィルクスのラムダ $\Lambda = \prod(1+\lambda_i)^{-1}$（最も標準的）
- ピライのトレース $V = \sum\lambda_i/(1+\lambda_i)$（最も頑健）
- ローリー=ホテリングのトレース $U = \sum\lambda_i$
- ロイの最大根 $\Theta = \lambda_1/(1+\lambda_1)$（1方向の差に最強の検出力）
- ANOVAとの関係: $p=1$ でANOVA、$K=2$ でホテリングのT²検定に特殊化される
- 前提条件: 多変量正規性、共分散行列の等質性、独立性、線形性
MANOVAが有意であった場合は、事後分析として判別分析やBonferroni補正付き個別ANOVAを行い、差の構造を詳しく調べることが推奨されます。
次のステップとして、以下の記事も参考にしてください。
- ホテリングのT²検定 — MANOVAの2群版
- 判別分析 — MANOVAと表裏一体の分類手法
- 多変量正規分布 — MANOVAの確率モデルの基盤