機械学習と情報技術

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

多変量解析

ある新薬の効果を臨床試験で評価するとしましょう。治療の効果は「収縮期血圧の変化」と「拡張期血圧の変化」の2つの指標で測定されます。この2つの指標について、プラセボ群と薬剤群の間に差があるかを検定したいとします。

一つのアプローチは、各指標ごとに個別にt検定を行うことです。しかし、この方法には2つの問題があります。第一に、2回の検定を行うことで多重検定の問題が生じます。各検定の有意水準を5%とすると、少なくとも一方で偽陽性が出る確率は $1 – (1-0.05)^2 \approx 9.75\%$ に膨らみます。エンドポイントが10個あれば、約40%にもなります。第二に、収縮期と拡張期の血圧は互いに強く相関しており、この相関構造を考慮しないことで検定の効率（検出力）が低下します。

ホテリングのT²検定（Hotelling’s T² test）は、$p$ 個の変数の平均ベクトルを同時に検定する多変量の検定手法です。1931年にハロルド・ホテリングが提案しました。この検定は1変数のt検定を多変量に自然に拡張したものであり、変数間の相関構造を共分散行列を通じて適切に考慮します。

ホテリングのT²検定を理解すると、以下のような幅広い応用が開けます。

  • 臨床試験: 複数のエンドポイント（主要・副次評価項目）を同時に検定し、全体として薬剤効果を評価する
  • 品質管理: 製品の複数の特性値（寸法・重量・強度）が仕様の目標値に合致しているかを同時に検定する（多変量管理図の理論的基盤）
  • 心理学・教育学: 複数の尺度スコア（IQ下位検査スコアなど）の群間比較
  • 環境科学: 複数の汚染物質濃度が基準値を超えているかの同時検定
  • 金融工学: 複数の資産リターンの平均が理論値と一致するかの検定

本記事では、1変数のt検定からの自然な拡張としてT²統計量を導出し、F分布との正確な関係を証明します。1標本・2標本・対応のある標本の3つのケースを解説し、信頼領域（信頼楕円）との関係も示します。Pythonで実装し、個別t検定との比較を通じてT²検定の利点を確認します。

本記事の内容

  • t検定からT²検定への自然な拡張
  • 1標本T²統計量の導出とF分布との関係の証明
  • マハラノビス距離との関係
  • 2標本T²検定の導出とプール共分散行列
  • 対応のあるT²検定
  • 信頼領域（信頼楕円）の導出
  • T²検定の検出力と標本サイズ設計
  • Pythonでの実装と可視化

前提知識

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

t検定からT²検定への拡張

t統計量の二乗形での理解

1変数のt検定を多変量に拡張するために、まずt統計量を「二乗形」で理解しましょう。

$n$ 個の独立なデータ $x_1, x_2, \dots, x_n$ が $N(\mu, \sigma^2)$ に従うとき、$H_0: \mu = \mu_0$ を検定するt統計量は:

$$ \begin{equation} t = \frac{\bar{x} – \mu_0}{s/\sqrt{n}} \end{equation} $$

ここで $\bar{x}$ は標本平均、$s^2 = \frac{1}{n-1}\sum_{i=1}^{n}(x_i – \bar{x})^2$ は不偏分散です。$H_0$ のもとで $t \sim t(n-1)$ に従います。

$t$ を二乗すると:

$$ \begin{equation} t^2 = n(\bar{x} – \mu_0) \cdot \frac{1}{s^2} \cdot (\bar{x} – \mu_0) = n\frac{(\bar{x} – \mu_0)^2}{s^2} \end{equation} $$

この式の構造に注目しましょう。$t^2$ は次の3つの要素から構成されています。

  1. サンプルサイズ $n$: データが多いほど統計量が大きくなる
  2. 平均の差 $(\bar{x} – \mu_0)$: 帰無仮説からのずれ
  3. 分散の逆数 $1/s^2$: データの散らばりで標準化

多変量への一般化

$t^2$ の各要素を多変量に拡張します。

1変数 多変量
スカラーの差 $(\bar{x} – \mu_0)$ ベクトルの差 $(\bar{\bm{x}} – \bm{\mu}_0)$
スカラーの分散 $s^2$ 共分散行列 $\bm{S}$
スカラーの分散の逆数 $1/s^2$ 共分散行列の逆行列 $\bm{S}^{-1}$
スカラー積 $(差) \times (1/s^2) \times (差)$ 二次形式 $(差)^T\bm{S}^{-1}(差)$

この対応により、T²統計量が自然に導かれます。

$$ \begin{equation} T^2 = n(\bar{\bm{x}} – \bm{\mu}_0)^T\bm{S}^{-1}(\bar{\bm{x}} – \bm{\mu}_0) \end{equation} $$

ここで $\bar{\bm{x}} = \frac{1}{n}\sum_{i=1}^{n}\bm{x}_i$ は標本平均ベクトル、$\bm{S} = \frac{1}{n-1}\sum_{i=1}^{n}(\bm{x}_i – \bar{\bm{x}})(\bm{x}_i – \bar{\bm{x}})^T$ は不偏標本共分散行列です。

$p = 1$（1変数）のとき、$\bm{S}^{-1} = 1/s^2$ となるので $T^2 = t^2$ に一致します。したがってT²検定はt検定の真の多変量拡張です。

マハラノビス距離との関係

T²統計量はマハラノビス距離と密接に関係しています。$\bar{\bm{x}}$ と $\bm{\mu}_0$ の間のマハラノビス距離の二乗は:

$$ \begin{equation} d_M^2(\bar{\bm{x}}, \bm{\mu}_0) = (\bar{\bm{x}} – \bm{\mu}_0)^T\bm{S}^{-1}(\bar{\bm{x}} – \bm{\mu}_0) \end{equation} $$

したがって $T^2 = n \cdot d_M^2(\bar{\bm{x}}, \bm{\mu}_0)$ です。

マハラノビス距離は、共分散構造を考慮した「標準化された距離」であり、変数間の相関と各変数のスケールの違いを自然に補正します。ユークリッド距離では、分散が大きい変数の差が過大評価されますが、マハラノビス距離ではこの問題が解消されます。

直感的には、T²は「標本平均が帰無仮説の平均からマハラノビス距離でどれだけ離れているか」を $n$ 倍にスケーリングしたものです。この距離が大きいほど、帰無仮説が正しくない証拠が強くなります。

T²統計量の意味が理解できたので、次にその帰無分布を厳密に導出しましょう。

T²統計量の帰無分布

定理の述べ方

定理: $\bm{x}_1, \bm{x}_2, \dots, \bm{x}_n$ が独立に $N_p(\bm{\mu}_0, \bm{\Sigma})$ に従うとき（帰無仮説 $H_0: \bm{\mu} = \bm{\mu}_0$ のもと）:

$$ \begin{equation} F = \frac{n – p}{p(n – 1)}T^2 \sim F(p, n – p) \end{equation} $$

証明

帰無仮説のもとで、標本平均の分布は:

$$ \bar{\bm{x}} \sim N_p\left(\bm{\mu}_0, \frac{\bm{\Sigma}}{n}\right) $$

したがって:

$$ \bm{z} = \sqrt{n}(\bar{\bm{x}} – \bm{\mu}_0) \sim N_p(\bm{0}, \bm{\Sigma}) $$

$\bm{\Sigma} = \bm{L}\bm{L}^T$（コレスキー分解）とすると、$\bm{u} = \bm{L}^{-1}\bm{z} \sim N_p(\bm{0}, \bm{I})$ は標準多変量正規分布に従います。

一方、$(n-1)\bm{S}$ はウィシャート分布 $W_p(\bm{\Sigma}, n-1)$ に従い、$\bar{\bm{x}}$ と独立です（正規分布における標本平均と標本共分散の独立性）。

T²統計量を $\bm{\Sigma}$ を用いて書き直すと:

$$ T^2 = n(\bar{\bm{x}} – \bm{\mu}_0)^T\bm{S}^{-1}(\bar{\bm{x}} – \bm{\mu}_0) = \bm{z}^T\bm{S}^{-1}\bm{z} $$

$\bm{S} = \frac{1}{n-1}\sum_{i}(\bm{x}_i – \bar{\bm{x}})(\bm{x}_i – \bar{\bm{x}})^T$ であり、$(n-1)\bm{S} \sim W_p(\bm{\Sigma}, n-1)$ なので、$\bm{S}^{-1} = (n-1)[(n-1)\bm{S}]^{-1}$ です。

多変量正規ベクトル $\bm{z}$ と独立なウィシャート行列 $(n-1)\bm{S}$ の比として構成されるT²は、ホテリングのT²分布に従うことが示されます。さらに、ホテリングのT²分布とF分布の間には次の正確な関係があります。

$$ F = \frac{n-p}{p(n-1)}T^2 $$

これが自由度 $(p, n-p)$ のF分布に従うことは、$\bm{z}$ と $(n-1)\bm{S}$ を同時対角化し、F分布の定義（独立な2つのカイ二乗分布の比）に帰着させることで証明できます。

特殊ケースの確認

$p = 1$ の場合: $T^2 = t^2$ であり:

$$ F = \frac{n-1}{1 \cdot (n-1)}T^2 = T^2 = t^2 $$

$F(1, n-1)$ は $t^2(n-1)$ に等しいため、1変数のt検定と一致します。

$\bm{\Sigma}$ が既知の場合: $\bm{S}$ の代わりに $\bm{\Sigma}$ を使うと:

$$ n(\bar{\bm{x}} – \bm{\mu}_0)^T\bm{\Sigma}^{-1}(\bar{\bm{x}} – \bm{\mu}_0) \sim \chi^2(p) $$

カイ二乗分布に従います。$\bm{S}$ で推定する不確かさが加わるため、カイ二乗分布からF分布に変わるのです。

1標本の理論が確立されたので、2標本の場合に拡張しましょう。

2標本T²検定

独立2群の平均の比較

2つの独立な群からデータが得られたとします。

  • 群1: $\bm{x}_{11}, \bm{x}_{12}, \dots, \bm{x}_{1n_1} \sim N_p(\bm{\mu}_1, \bm{\Sigma})$
  • 群2: $\bm{x}_{21}, \bm{x}_{22}, \dots, \bm{x}_{2n_2} \sim N_p(\bm{\mu}_2, \bm{\Sigma})$

等共分散の仮定: 2つの群の母共分散行列が等しい（$\bm{\Sigma}_1 = \bm{\Sigma}_2 = \bm{\Sigma}$）と仮定します。この仮定は1変数の等分散性の仮定の多変量版です。

帰無仮説は $H_0: \bm{\mu}_1 = \bm{\mu}_2$（2群の母平均ベクトルが等しい）です。

プール共分散行列

等共分散の仮定のもとで、共通の $\bm{\Sigma}$ を推定するために、2群のデータをプール（統合）して共分散行列を推定します。

$$ \begin{equation} \bm{S}_p = \frac{(n_1 – 1)\bm{S}_1 + (n_2 – 1)\bm{S}_2}{n_1 + n_2 – 2} \end{equation} $$

ここで $\bm{S}_1$ と $\bm{S}_2$ はそれぞれの群の標本共分散行列です。プール共分散行列は $n_1 + n_2 – 2$ の自由度を持つ $\bm{\Sigma}$ の不偏推定量です。各群の共分散行列をサンプルサイズで重み付け平均しているので、サンプルサイズが大きい群の共分散行列により大きな重みが置かれます。

T²統計量とF分布

2標本のT²統計量は:

$$ \begin{equation} T^2 = \frac{n_1 n_2}{n_1 + n_2}(\bar{\bm{x}}_1 – \bar{\bm{x}}_2)^T\bm{S}_p^{-1}(\bar{\bm{x}}_1 – \bar{\bm{x}}_2) \end{equation} $$

$n_1 n_2/(n_1 + n_2)$ という係数は、1変数の2標本t検定における $1/(1/n_1 + 1/n_2)$ に対応するものです。

F分布への変換は:

$$ \begin{equation} F = \frac{n_1 + n_2 – p – 1}{p(n_1 + n_2 – 2)}T^2 \sim F(p, n_1 + n_2 – p – 1) \end{equation} $$

第1自由度は $p$（変数の数）、第2自由度は $n_1 + n_2 – p – 1$ です。

等共分散の仮定が成り立たない場合

等共分散の仮定 $\bm{\Sigma}_1 = \bm{\Sigma}_2$ が成り立たない場合（多変量のBehrens-Fisher問題）には、正確な帰無分布を導出するのが困難になります。実用的な対処法として:

  1. BoxのM検定: $H_0: \bm{\Sigma}_1 = \bm{\Sigma}_2$ を検定できますが、正規性の仮定に敏感でサンプルサイズが大きいと過度に有意になる傾向があります
  2. 修正T²検定: WelchのtをT²に拡張したBehrens-Fisher近似があります。自由度を修正したF近似を使います
  3. ノンパラメトリック法: 分布の仮定を置かない方法（置換検定など）を使います

次に、対応のある標本（対応のある2群）の場合を見ましょう。

対応のあるT²検定

定式化

対応のある2群（例えば治療前後の同一被験者のデータ）では、差のベクトル $\bm{d}_i = \bm{x}_{1i} – \bm{x}_{2i}$ を計算し、$\bm{d}_i$ に対して1標本T²検定を適用します。

帰無仮説: $H_0: \bm{\mu}_d = \bm{0}$（差の平均がゼロ）

$$ \begin{equation} T^2 = n\bar{\bm{d}}^T\bm{S}_d^{-1}\bar{\bm{d}} \end{equation} $$

ここで $\bar{\bm{d}} = \frac{1}{n}\sum_{i=1}^{n}\bm{d}_i$ は差の標本平均、$\bm{S}_d$ は差の標本共分散行列、$n$ は対の数です。

F変換は1標本の場合と同じです。

$$ F = \frac{n-p}{p(n-1)}T^2 \sim F(p, n-p) $$

対応のある検定は、個体差（各被験者の基本レベルの違い）を除去できるため、独立2標本検定よりも検出力が高い場合が多いです。

信頼領域（信頼楕円）

同時信頼領域の導出

T²検定の非棄却域は、母平均ベクトル $\bm{\mu}$ の同時信頼領域に対応します。信頼水準 $1 – \alpha$ の信頼領域は:

$$ \begin{equation} \left\{\bm{\mu} : n(\bar{\bm{x}} – \bm{\mu})^T\bm{S}^{-1}(\bar{\bm{x}} – \bm{\mu}) \leq \frac{p(n-1)}{n-p}F_\alpha(p, n-p)\right\} \end{equation} $$

この領域は $p$ 次元空間における楕円体です。$p = 2$ のとき、楕円の形状と方向は共分散行列 $\bm{S}$ によって決まります。

楕円体の主軸方向は $\bm{S}$ の固有ベクトル、各主軸の半径は対応する固有値の平方根に比例します。したがって、変数間の相関が強いほど楕円は細長くなり、相関がない場合は軸に平行な楕円になります。

同時信頼区間

信頼楕円から、各変数の同時信頼区間（Bonferroni区間よりも広い）を導出できます。任意の方向 $\bm{a}$ に対する $\bm{a}^T\bm{\mu}$ の同時信頼区間は:

$$ \bm{a}^T\bar{\bm{x}} \pm \sqrt{\frac{p(n-1)}{n(n-p)}F_\alpha(p, n-p)} \cdot \sqrt{\bm{a}^T\bm{S}\bm{a}} $$

$\bm{a} = \bm{e}_j$（第 $j$ 単位ベクトル）とすれば、$\mu_j$ の同時信頼区間が得られます。

これらの理論をPythonで実装してみましょう。

Pythonによる実装と可視化

1標本T²検定の実装

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

np.random.seed(42)

# --- データ生成 ---
n = 50
p = 3
mu_true = np.array([1.5, -0.8, 0.3])
Sigma = np.array([[1.0, 0.3, 0.1],
                   [0.3, 1.5, 0.2],
                   [0.1, 0.2, 0.8]])
X = np.random.multivariate_normal(mu_true, Sigma, n)

# --- 帰無仮説: μ = 0 ---
mu0 = np.array([0, 0, 0])

# --- T²統計量の計算 ---
x_bar = X.mean(axis=0)
S = np.cov(X.T, ddof=1)
diff = x_bar - mu0
T2 = n * diff @ np.linalg.inv(S) @ diff

# --- F変換 ---
F_stat = (n - p) / (p * (n - 1)) * T2
p_value = 1 - stats.f.cdf(F_stat, p, n - p)

print("=== 1標本ホテリングのT²検定 ===")
print(f"帰無仮説: μ = {mu0}")
print(f"標本平均: [{x_bar[0]:.4f}, {x_bar[1]:.4f}, {x_bar[2]:.4f}]")
print(f"\nT² = {T2:.4f}")
print(f"F({p}, {n-p}) = {F_stat:.4f}")
print(f"p値 = {p_value:.2e}")
print(f"結論: {'H0を棄却' if p_value < 0.05 else 'H0を棄却しない'} (α=0.05)")

# --- 個別t検定との比較 ---
print("\n--- 各変数の個別t検定 ---")
for j in range(p):
    t_stat, p_val = stats.ttest_1samp(X[:, j], mu0[j])
    sig = '*' if p_val < 0.05 else 'n.s.'
    print(f"  x{j+1}: t = {t_stat:7.4f}, p = {p_val:.6f} {sig}")

# Bonferroni補正
alpha_bonf = 0.05 / p
print(f"\n--- Bonferroni補正後 (α/p = {alpha_bonf:.4f}) ---")
for j in range(p):
    _, p_val = stats.ttest_1samp(X[:, j], mu0[j])
    sig = '*' if p_val < alpha_bonf else 'n.s.'
    print(f"  x{j+1}: p = {p_val:.6f} {sig}")

この結果の比較から、T²検定の重要な特徴が確認できます。T²検定は全変数を同時に考慮するため、多重検定の問題を回避しつつ、変数間の相関を利用して検出力を維持します。個別t検定で有意でない変数があっても、T²検定では全体として有意になることがあります。これは、各変数単独では小さな差でも、多変量の相関構造を考慮すると統計的に有意な差となりうることを示しています。

2標本T²検定と信頼楕円の可視化

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse

np.random.seed(42)

# --- 2つの群のデータ生成 ---
n1, n2 = 40, 45
p = 2
mu1 = np.array([2.0, 3.0])
mu2 = np.array([3.5, 4.5])
Sigma = np.array([[1.5, 0.6], [0.6, 1.0]])

X1 = np.random.multivariate_normal(mu1, Sigma, n1)
X2 = np.random.multivariate_normal(mu2, Sigma, n2)

# --- 2標本T²検定 ---
x1_bar = X1.mean(axis=0)
x2_bar = X2.mean(axis=0)
S1 = np.cov(X1.T, ddof=1)
S2 = np.cov(X2.T, ddof=1)
S_pooled = ((n1-1)*S1 + (n2-1)*S2) / (n1+n2-2)

diff = x1_bar - x2_bar
T2 = (n1*n2)/(n1+n2) * diff @ np.linalg.inv(S_pooled) @ diff
F_stat = (n1+n2-p-1) / (p*(n1+n2-2)) * T2
df1, df2 = p, n1+n2-p-1
p_value = 1 - stats.f.cdf(F_stat, df1, df2)

print("=== 2標本ホテリングのT²検定 ===")
print(f"群1 (n={n1}): 平均 = {x1_bar.round(3)}")
print(f"群2 (n={n2}): 平均 = {x2_bar.round(3)}")
print(f"差 = {diff.round(3)}")
print(f"\nT² = {T2:.4f}")
print(f"F({df1}, {df2}) = {F_stat:.4f}")
print(f"p値 = {p_value:.2e}")
print(f"結論: {'H0を棄却' if p_value < 0.05 else 'H0を棄却しない'}")

# --- 可視化 ---
fig, axes = plt.subplots(1, 3, figsize=(17, 5.5))

# (a) 散布図と信頼楕円
ax = axes[0]
ax.scatter(X1[:, 0], X1[:, 1], c='#377eb8', alpha=0.4, s=20, label='Group 1')
ax.scatter(X2[:, 0], X2[:, 1], c='#e41a1c', alpha=0.4, s=20, label='Group 2')
ax.plot(*x1_bar, 'b+', markersize=15, markeredgewidth=2.5)
ax.plot(*x2_bar, 'r+', markersize=15, markeredgewidth=2.5)

# 95%信頼楕円
eigvals, eigvecs = np.linalg.eigh(S_pooled)
angle = np.degrees(np.arctan2(eigvecs[1, 1], eigvecs[0, 1]))
chi2_val = stats.chi2.ppf(0.95, p)

for mu_k, color, n_k in [(x1_bar, '#377eb8', n1), (x2_bar, '#e41a1c', n2)]:
    w = 2 * np.sqrt(chi2_val * eigvals[1] / n_k)
    h = 2 * np.sqrt(chi2_val * eigvals[0] / n_k)
    ell = Ellipse(xy=mu_k, width=w, height=h, angle=angle,
                  edgecolor=color, facecolor='none', linewidth=1.5, linestyle='--')
    ax.add_patch(ell)

ax.set_xlabel('$x_1$', fontsize=11)
ax.set_ylabel('$x_2$', fontsize=11)
ax.set_title(f'(a) Data and 95% confidence ellipses', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# (b) T²統計量の帰無分布
ax = axes[1]
x_range = np.linspace(0, 30, 300)
scale = p*(n1+n2-2)/(n1+n2-p-1)
f_range = x_range / scale
pdf_vals = stats.f.pdf(f_range, df1, df2) / scale
ax.plot(x_range, pdf_vals, 'b-', linewidth=2, label='Null distribution')
ax.axvline(T2, color='red', linewidth=2, linestyle='--',
           label=f'$T^2 = {T2:.2f}$')
ax.fill_between(x_range, pdf_vals, where=x_range >= T2,
                alpha=0.3, color='red')
ax.set_xlabel('$T^2$', fontsize=11)
ax.set_ylabel('Density', fontsize=11)
ax.set_title(f'(b) Null distribution (p = {p_value:.2e})', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)

# (c) 差の信頼楕円
ax = axes[2]
# 差の共分散行列
S_diff = S_pooled * (1/n1 + 1/n2)
eigvals_d, eigvecs_d = np.linalg.eigh(S_diff)
angle_d = np.degrees(np.arctan2(eigvecs_d[1, 1], eigvecs_d[0, 1]))
F_crit = stats.f.ppf(0.95, p, n1+n2-p-1)
scale_factor = p*(n1+n2-2)/(n1+n2-p-1) * F_crit

w_d = 2 * np.sqrt(scale_factor * eigvals_d[1])
h_d = 2 * np.sqrt(scale_factor * eigvals_d[0])
ell_d = Ellipse(xy=diff, width=w_d, height=h_d, angle=angle_d,
                edgecolor='green', facecolor='lightgreen', alpha=0.3, linewidth=2)
ax.add_patch(ell_d)
ax.plot(*diff, 'g+', markersize=15, markeredgewidth=2.5, label='Observed difference')
ax.plot(0, 0, 'ko', markersize=8, label='H0: difference = 0')
ax.set_xlabel('$\\bar{x}_{1,1} - \\bar{x}_{2,1}$', fontsize=11)
ax.set_ylabel('$\\bar{x}_{1,2} - \\bar{x}_{2,2}$', fontsize=11)
ax.set_title('(c) 95% confidence region for difference', fontsize=12)
ax.legend(fontsize=9)
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

3つのプロットから、T²検定の理論と実践の対応が明確に確認できます。

  1. 散布図と信頼楕円（図a）: 2群のデータが2次元空間上で分布し、各群の標本平均の周りに95%信頼楕円が描かれています。2つの楕円が重ならないことは、2群の平均が統計的に有意に異なることを視覚的に示しています。楕円の形状（傾き）は共分散行列の固有構造を反映しており、$x_1$ と $x_2$ の正の相関が楕円を右上がりに傾けています

  2. 帰無分布（図b）: T²統計量の帰無分布（スケーリングされたF分布）と、観測されたT²値（赤い破線）が表示されています。赤く塗られた領域がp値に対応します。T²値が分布の右裾に位置していることは、帰無仮説のもとではこのような大きな値が観測される確率が非常に低いことを示しています

  3. 差の信頼楕円（図c）: 2群の平均の差 $\bar{\bm{x}}_1 – \bar{\bm{x}}_2$ に対する95%信頼楕円が描かれています。帰無仮説の値（原点: 差がゼロ）が楕円の外にあれば、帰無仮説は5%水準で棄却されます。このプロットから、信頼領域と仮説検定の双対性が視覚的に確認できます

検出力シミュレーション

T²検定と個別t検定（Bonferroni補正付き）の検出力を比較します。

import numpy as np
from scipy import stats

np.random.seed(42)

def hotelling_t2_test(X1, X2, alpha=0.05):
    """2標本T²検定"""
    n1, n2 = len(X1), len(X2)
    p = X1.shape[1]
    S1 = np.cov(X1.T, ddof=1)
    S2 = np.cov(X2.T, ddof=1)
    Sp = ((n1-1)*S1 + (n2-1)*S2) / (n1+n2-2)
    diff = X1.mean(0) - X2.mean(0)
    T2 = n1*n2/(n1+n2) * diff @ np.linalg.solve(Sp, diff)
    F = (n1+n2-p-1) / (p*(n1+n2-2)) * T2
    pval = 1 - stats.f.cdf(F, p, n1+n2-p-1)
    return pval < alpha

def bonferroni_ttest(X1, X2, alpha=0.05):
    """個別t検定（Bonferroni補正）"""
    p = X1.shape[1]
    for j in range(p):
        _, pval = stats.ttest_ind(X1[:, j], X2[:, j])
        if pval < alpha / p:
            return True
    return False

# パラメータ
p = 3
n_per_group = 30
Sigma = np.array([[1, 0.5, 0.3], [0.5, 1, 0.4], [0.3, 0.4, 1]])
n_sim = 2000

# 効果量を変えて検出力を計算
effect_sizes = np.linspace(0, 1.5, 12)
power_T2 = []
power_bonf = []

for delta in effect_sizes:
    mu_diff = np.array([delta, delta*0.7, delta*0.5])
    count_T2 = 0
    count_bonf = 0
    for _ in range(n_sim):
        X1 = np.random.multivariate_normal(np.zeros(p), Sigma, n_per_group)
        X2 = np.random.multivariate_normal(mu_diff, Sigma, n_per_group)
        count_T2 += hotelling_t2_test(X1, X2)
        count_bonf += bonferroni_ttest(X1, X2)
    power_T2.append(count_T2 / n_sim)
    power_bonf.append(count_bonf / n_sim)

fig, ax = plt.subplots(figsize=(9, 5.5))
ax.plot(effect_sizes, power_T2, 'bo-', linewidth=2, markersize=7,
        label="Hotelling's T²")
ax.plot(effect_sizes, power_bonf, 'rs--', linewidth=2, markersize=7,
        label="Individual t-tests (Bonferroni)")
ax.axhline(0.05, color='gray', linewidth=0.8, linestyle=':', label='α = 0.05')
ax.axhline(0.80, color='green', linewidth=0.8, linestyle=':', label='Power = 0.80')
ax.set_xlabel('Effect size $\\delta$', fontsize=12)
ax.set_ylabel('Power', fontsize=12)
ax.set_title('Power comparison: T² vs Bonferroni t-tests', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim(-0.02, 1.05)
plt.tight_layout()
plt.show()

検出力の比較から、T²検定の優位性が定量的に確認できます。

  1. 効果量ゼロ（$\delta = 0$）: 両手法とも棄却率が約5%（名目有意水準）であり、第1種の過誤率が適切に制御されています

  2. 小〜中程度の効果量: T²検定はBonferroni補正付きt検定よりも一貫して高い検出力を示しています。この差は、T²検定が変数間の相関を考慮して「情報を効率的に統合する」ことに起因します

  3. 大きな効果量: 両手法の検出力はともに1に近づきますが、T²検定の方が早く高い検出力に達します。T²検定は80%検出力を達成するために必要な効果量が小さいことがわかります

この検出力の優位性は、特に変数間の相関が強い場合に顕著になります。相関が強いほど、T²検定は個別検定と比べて効率的に情報を統合できるためです。

まとめ

本記事では、ホテリングのT²検定の理論を体系的に解説しました。

  • T²統計量: $T^2 = n(\bar{\bm{x}} – \bm{\mu}_0)^T\bm{S}^{-1}(\bar{\bm{x}} – \bm{\mu}_0)$ はt統計量の二乗の多変量拡張であり、マハラノビス距離の二乗に $n$ を掛けたもの
  • F分布との関係: $F = \frac{n-p}{p(n-1)}T^2 \sim F(p, n-p)$ が帰無仮説のもとで厳密に成り立つ
  • 2標本T²検定: プール共分散行列 $\bm{S}_p$ を用いて2群の平均ベクトルの差を検定する。$F = \frac{n_1+n_2-p-1}{p(n_1+n_2-2)}T^2 \sim F(p, n_1+n_2-p-1)$
  • 信頼領域: T²検定の非棄却域は母平均の信頼楕円体に対応する
  • 検出力: T²検定は個別t検定+Bonferroni補正よりも一般に高い検出力を持つ

T²検定と判別分析の関係

T²検定は判別分析とも密接に関連しています。2群の判別分析では、フィッシャーの線形判別関数 $\bm{a} = \bm{S}_p^{-1}(\bar{\bm{x}}_1 – \bar{\bm{x}}_2)$ を求めます。T²統計量はこの判別関数の方向にデータを射影したときのt統計量の二乗に他なりません。

$$ T^2 = \frac{n_1 n_2}{n_1 + n_2}\bm{a}^T\bm{S}_p\bm{a} \cdot \left(\frac{\bm{a}^T(\bar{\bm{x}}_1 – \bar{\bm{x}}_2)}{\bm{a}^T\bm{S}_p\bm{a}}\right)^2 \cdot \frac{1}{1/n_1 + 1/n_2} $$

つまり、T²検定は「2群の分離を最大化する方向に射影した上でt検定を行う」操作と等価です。この方向は判別分析が見つける最適な方向と一致します。

T²検定の限界と注意点

T²検定にはいくつかの限界があります。

サンプルサイズの要件: $\bm{S}^{-1}$ の計算に $n > p$ が必要です。$p$ が $n$ に近いと $\bm{S}$ の推定精度が低下し、検出力が大幅に落ちます。経験則として $n \geq 5p$ が推奨されます。

正規性への依存: T²検定は多変量正規性を仮定しています。裾の重い分布や歪んだ分布では、ロバスト版のT²検定やノンパラメトリック検定（Cramér検定や置換検定）が推奨されます。

等共分散の仮定（2標本）: 2群の共分散行列が等しいという仮定が必要です。この仮定が破られると、特にサンプルサイズが不均等な場合に第1種の過誤率が名目水準から乖離します。

T²検定はMANOVAの特殊ケース（$K = 2$）でもあります。3群以上の平均の比較にはMANOVAを用いる必要があります。

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

関連記事

前の記事

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

次の記事

対応分析の理論と解釈

関連タグ: F検定 ホテリングT2検定 仮説検定 多変量検定 多変量正規分布 多変量解析