二元配置分散分析の理論と交互作用をわかりやすく解説

ある作物の収量を改善したいとき、肥料の種類と灌漑方法の2つの要因を同時に変えて実験することがあります。このとき、「肥料の効果」と「灌漑の効果」をそれぞれ独立に評価したいのはもちろんですが、もう一つ重要な問いがあります。それは「特定の肥料と特定の灌漑方法の組み合わせが、それぞれの効果の単純な足し算以上の効果を持つか」という問いです。

この「足し算以上(または以下)の効果」を交互作用(interaction)と呼びます。交互作用が存在すると、一方の要因の効果が他方の要因の水準によって変わります。たとえば、肥料Aは通常の灌漑では効果が小さいが、点滴灌漑と組み合わせると劇的な効果を発揮する、という状況です。

このような2つの要因の効果と交互作用を同時に分析する手法が二元配置分散分析(two-way ANOVA)です。

二元配置分散分析を理解すると、以下のような応用が可能です。

  • 実験計画法: 2つの因子を同時に変える要因実験(factorial experiment)の解析の基礎です
  • 品質工学: 製造条件(温度×圧力、材料×処理方法など)の最適化に不可欠です
  • 医学研究: 治療法と患者の特性(年齢群、性別など)の交互作用を評価します
  • 教育研究: 教授法と学習者の特性の組み合わせ効果を分析します

本記事の内容

  • 交互作用の概念と直感的な理解
  • 二元配置分散分析のモデルと平方和の分解
  • F統計量の導出と検定手順
  • 繰り返しなし・繰り返しありの場合の違い
  • Pythonによる実装と交互作用プロットの可視化

前提知識

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

交互作用の直感的な理解

交互作用がない場合 — 加法モデル

まず、交互作用がない単純なケースを考えましょう。

2種類の肥料(A, B)と2種類の灌漑(標準, 点滴)で作物を育てる実験を想像してください。交互作用がない場合、各セルの期待値は次のような加法的な構造を持ちます。

標準灌漑 点滴灌漑 灌漑の効果
肥料A 10 15 +5
肥料B 14 19 +5
肥料の効果 +4 +4

この例では、灌漑を点滴に変えると収量が $+5$ 増え、肥料をBに変えると $+4$ 増えます。重要なのは、灌漑の効果 $+5$ は肥料の種類に依存せず、肥料の効果 $+4$ は灌漑方法に依存しません。2つの要因の効果が独立に加算されています。

交互作用がある場合

次に、交互作用がある場合を考えましょう。

標準灌漑 点滴灌漑 灌漑の効果
肥料A 10 13 +3
肥料B 14 22 +8
肥料の効果 +4 +9

今度は、灌漑の効果が肥料の種類によって異なります。肥料Aでは灌漑を変えても $+3$ しか増えませんが、肥料Bでは $+8$ も増えます。同様に、肥料の効果も灌漑方法によって異なります。

この「一方の要因の効果が他方の要因の水準によって変わる」という現象が交互作用です。交互作用があると、要因の主効果だけでは結果を正しく説明できません。「肥料Bは全体として良い」というだけでなく、「肥料Bは点滴灌漑と組み合わせると特に効果的」という情報が重要になります。

交互作用プロットによる視覚的判断

交互作用の有無を視覚的に判断する最も直感的な方法は、交互作用プロット(interaction plot)です。一方の要因を横軸に、各セルの平均値を縦軸にとり、他方の要因の各水準を別の線で描きます。

  • 交互作用なし: 線が平行になる
  • 交互作用あり: 線が交差するか、非平行になる

線の非平行度合いが交互作用の大きさに対応します。

交互作用の概念を理解したところで、二元配置ANOVAの数学的モデルを定式化しましょう。

二元配置ANOVAのモデル

記号の定義

  • 因子Aは $a$ 個の水準を持つ($i = 1, 2, \dots, a$)
  • 因子Bは $b$ 個の水準を持つ($j = 1, 2, \dots, b$)
  • 各セル $(i, j)$ に $r$ 個の繰り返し観測がある($k = 1, 2, \dots, r$)
  • 観測値を $Y_{ijk}$ と書く
  • 総観測数は $N = abr$

構造モデル

二元配置ANOVAの構造モデルは次のように表されます。

$$ \begin{equation} Y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \varepsilon_{ijk} \end{equation} $$

各項の意味は以下の通りです。

  • $\mu$: 全体平均(grand mean)
  • $\alpha_i$: 因子Aの第 $i$ 水準の主効果($\sum_{i=1}^{a}\alpha_i = 0$)
  • $\beta_j$: 因子Bの第 $j$ 水準の主効果($\sum_{j=1}^{b}\beta_j = 0$)
  • $(\alpha\beta)_{ij}$: 交互作用効果($\sum_i (\alpha\beta)_{ij} = 0$, $\sum_j (\alpha\beta)_{ij} = 0$)
  • $\varepsilon_{ijk}$: 誤差項($\varepsilon_{ijk} \overset{\text{i.i.d.}}{\sim} N(0, \sigma^2)$)

制約条件 $\sum \alpha_i = 0$, $\sum \beta_j = 0$, $\sum_i (\alpha\beta)_{ij} = 0$, $\sum_j (\alpha\beta)_{ij} = 0$ はパラメータの一意性を保証するためのものです。

3つの帰無仮説

二元配置ANOVAでは、3つの帰無仮説を同時に検定します。

  1. 因子Aの主効果: $H_{0A}: \alpha_1 = \alpha_2 = \cdots = \alpha_a = 0$(因子Aのすべての水準の効果が等しい)
  2. 因子Bの主効果: $H_{0B}: \beta_1 = \beta_2 = \cdots = \beta_b = 0$(因子Bのすべての水準の効果が等しい)
  3. 交互作用: $H_{0AB}: (\alpha\beta)_{ij} = 0$ for all $i, j$(交互作用がない)

通常、まず交互作用の検定を行い、交互作用が有意であれば主効果の解釈には注意が必要です。

各効果を検定するために、全変動を各成分に分解しましょう。

平方和の分解

各平均値の定義

以下の平均値を定義します。

$$ \bar{Y}_{i\cdot\cdot} = \frac{1}{br}\sum_{j=1}^{b}\sum_{k=1}^{r} Y_{ijk} \quad \text{(因子Aの第}i\text{水準の平均)} $$

$$ \bar{Y}_{\cdot j \cdot} = \frac{1}{ar}\sum_{i=1}^{a}\sum_{k=1}^{r} Y_{ijk} \quad \text{(因子Bの第}j\text{水準の平均)} $$

$$ \bar{Y}_{ij\cdot} = \frac{1}{r}\sum_{k=1}^{r} Y_{ijk} \quad \text{(セル}(i,j)\text{の平均)} $$

$$ \bar{Y}_{\cdots} = \frac{1}{N}\sum_{i,j,k} Y_{ijk} \quad \text{(全体平均)} $$

全変動の分解

全変動(Total Sum of Squares)を各成分に分解します。

$$ \begin{equation} SS_T = SS_A + SS_B + SS_{AB} + SS_E \end{equation} $$

各平方和は次のように定義されます。

全変動:

$$ SS_T = \sum_{i=1}^{a}\sum_{j=1}^{b}\sum_{k=1}^{r}(Y_{ijk} – \bar{Y}_{\cdots})^2 $$

因子Aの主効果による変動:

$$ SS_A = br\sum_{i=1}^{a}(\bar{Y}_{i\cdot\cdot} – \bar{Y}_{\cdots})^2 $$

因子Bの主効果による変動:

$$ SS_B = ar\sum_{j=1}^{b}(\bar{Y}_{\cdot j\cdot} – \bar{Y}_{\cdots})^2 $$

交互作用による変動:

$$ SS_{AB} = r\sum_{i=1}^{a}\sum_{j=1}^{b}(\bar{Y}_{ij\cdot} – \bar{Y}_{i\cdot\cdot} – \bar{Y}_{\cdot j\cdot} + \bar{Y}_{\cdots})^2 $$

残差(誤差)変動:

$$ SS_E = \sum_{i=1}^{a}\sum_{j=1}^{b}\sum_{k=1}^{r}(Y_{ijk} – \bar{Y}_{ij\cdot})^2 $$

分解の証明のスケッチ

全変動の分解は、恒等式

$$ Y_{ijk} – \bar{Y}_{\cdots} = (\bar{Y}_{i\cdot\cdot} – \bar{Y}_{\cdots}) + (\bar{Y}_{\cdot j\cdot} – \bar{Y}_{\cdots}) + (\bar{Y}_{ij\cdot} – \bar{Y}_{i\cdot\cdot} – \bar{Y}_{\cdot j\cdot} + \bar{Y}_{\cdots}) + (Y_{ijk} – \bar{Y}_{ij\cdot}) $$

の両辺を2乗して全 $i, j, k$ について和を取ることで得られます。交差項がすべてゼロになることを示す必要がありますが、これは各項の直交性(例えば $\sum_i (\bar{Y}_{i\cdot\cdot} – \bar{Y}_{\cdots}) = 0$)から証明できます。

自由度

各平方和の自由度は次の通りです。

$$ \underbrace{N – 1}_{SS_T} = \underbrace{a – 1}_{SS_A} + \underbrace{b – 1}_{SS_B} + \underbrace{(a-1)(b-1)}_{SS_{AB}} + \underbrace{ab(r-1)}_{SS_E} $$

自由度もちょうど加法的に分解されています。

平方和を自由度で割ったものが平均平方(mean square)です。

$$ MS_A = \frac{SS_A}{a-1}, \quad MS_B = \frac{SS_B}{b-1}, \quad MS_{AB} = \frac{SS_{AB}}{(a-1)(b-1)}, \quad MS_E = \frac{SS_E}{ab(r-1)} $$

これらの平均平方からF統計量を構成しましょう。

F統計量と検定

各効果のF統計量

各帰無仮説に対するF統計量は次の通りです。

$$ F_A = \frac{MS_A}{MS_E}, \quad F_B = \frac{MS_B}{MS_E}, \quad F_{AB} = \frac{MS_{AB}}{MS_E} $$

帰無仮説のもとで、各F統計量はそれぞれ次のF分布に従います。

$$ F_A \sim F(a-1,\, ab(r-1)) $$

$$ F_B \sim F(b-1,\, ab(r-1)) $$

$$ F_{AB} \sim F((a-1)(b-1),\, ab(r-1)) $$

F統計量の直感

F統計量は「要因による変動 / 偶然の変動」の比です。もし因子Aに効果があれば、$MS_A$ は $MS_E$ よりも大きくなり、$F_A$ は1よりも大きくなります。$F_A$ が十分大きければ、因子Aの効果は「偶然では説明できない」と判断し、帰無仮説を棄却します。

ANOVA表

結果は通常、以下のANOVA表にまとめます。

要因 平方和 自由度 平均平方 F値 p値
因子A $SS_A$ $a-1$ $MS_A$ $F_A$
因子B $SS_B$ $b-1$ $MS_B$ $F_B$
A×B $SS_{AB}$ $(a-1)(b-1)$ $MS_{AB}$ $F_{AB}$
残差 $SS_E$ $ab(r-1)$ $MS_E$
全体 $SS_T$ $N-1$

検定の解釈の順序

二元配置ANOVAの結果を解釈する際は、次の順序で進めることが推奨されます。

  1. まず交互作用 $F_{AB}$ を確認する
  2. 交互作用が有意でない場合: 主効果 $F_A$, $F_B$ をそれぞれ解釈する
  3. 交互作用が有意な場合: 主効果の解釈には注意が必要。交互作用が存在すると、主効果の平均は実際の効果を正しく反映しない可能性がある。この場合は、単純主効果の分析(simple main effects analysis)を行い、一方の因子の各水準における他方の因子の効果を個別に検定する

繰り返しがない場合の特殊性について、次に見ていきましょう。

繰り返しなしの場合

$r = 1$ の問題

各セルに1つの観測値しかない場合($r = 1$)、$SS_E = 0$ となり、分母がゼロになって検定が実行できません。これは、セル内の変動(誤差変動)を推定するためのデータがないためです。

交互作用なしの仮定

$r = 1$ の場合の標準的な対処法は、交互作用がないと仮定することです。

$$ Y_{ij} = \mu + \alpha_i + \beta_j + \varepsilon_{ij} $$

この仮定のもとでは、$SS_{AB}$ を誤差の推定量として使います。

$$ F_A = \frac{MS_A}{MS_{AB}}, \quad F_B = \frac{MS_B}{MS_{AB}} $$

ここで $MS_{AB} = \frac{SS_{AB}}{(a-1)(b-1)}$ は誤差分散の推定量として機能します。

ただし、実際に交互作用が存在する場合、この方法は適切ではありません。交互作用の存在を確認するにはテューキーの加法性検定(Tukey’s test for additivity)を使うことができます。

Pythonによる実装

二元配置ANOVAのスクラッチ実装

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

def two_way_anova(data, a, b, r):
    """
    二元配置分散分析(等繰り返し)のスクラッチ実装

    Parameters
    ----------
    data : ndarray, shape (a, b, r)  観測データ
    a, b : int  因子Aと因子Bの水準数
    r : int  繰り返し数

    Returns
    -------
    anova_table : dict  ANOVA表の各値
    """
    N = a * b * r
    grand_mean = np.mean(data)

    # 各水準の平均
    mean_A = np.mean(data, axis=(1, 2))     # shape (a,)
    mean_B = np.mean(data, axis=(0, 2))     # shape (b,)
    mean_AB = np.mean(data, axis=2)          # shape (a, b)

    # 平方和
    SS_A = b * r * np.sum((mean_A - grand_mean)**2)
    SS_B = a * r * np.sum((mean_B - grand_mean)**2)
    SS_AB = r * np.sum((mean_AB - mean_A[:, None] - mean_B[None, :] + grand_mean)**2)
    SS_T = np.sum((data - grand_mean)**2)
    SS_E = SS_T - SS_A - SS_B - SS_AB

    # 自由度
    df_A = a - 1
    df_B = b - 1
    df_AB = (a - 1) * (b - 1)
    df_E = a * b * (r - 1)

    # 平均平方
    MS_A = SS_A / df_A
    MS_B = SS_B / df_B
    MS_AB = SS_AB / df_AB
    MS_E = SS_E / df_E

    # F値
    F_A = MS_A / MS_E
    F_B = MS_B / MS_E
    F_AB = MS_AB / MS_E

    # p値
    p_A = 1 - stats.f.cdf(F_A, df_A, df_E)
    p_B = 1 - stats.f.cdf(F_B, df_B, df_E)
    p_AB = 1 - stats.f.cdf(F_AB, df_AB, df_E)

    return {
        "SS": [SS_A, SS_B, SS_AB, SS_E, SS_T],
        "df": [df_A, df_B, df_AB, df_E, N - 1],
        "MS": [MS_A, MS_B, MS_AB, MS_E, None],
        "F": [F_A, F_B, F_AB, None, None],
        "p": [p_A, p_B, p_AB, None, None],
        "means": {"A": mean_A, "B": mean_B, "AB": mean_AB, "grand": grand_mean}
    }

# --- テストデータの生成 ---
np.random.seed(42)
a, b, r = 3, 2, 5  # 3水準×2水準×5繰り返し

# 真のパラメータ
mu = 50
alpha = np.array([-5, 0, 5])       # 因子Aの効果
beta = np.array([-3, 3])            # 因子Bの効果
interaction = np.array([             # 交互作用
    [-2, 2],
    [0, 0],
    [2, -2]
])

# データ生成
data = np.zeros((a, b, r))
for i in range(a):
    for j in range(b):
        data[i, j, :] = (mu + alpha[i] + beta[j] + interaction[i, j]
                          + np.random.randn(r) * 3)

# ANOVA
result = two_way_anova(data, a, b, r)

# 結果表示
print("=" * 65)
print("二元配置分散分析表")
print("=" * 65)
labels = ["因子A", "因子B", "A×B", "残差", "全体"]
print(f"{'要因':<8} {'SS':>10} {'df':>5} {'MS':>10} {'F':>10} {'p':>10}")
print("-" * 65)
for i, label in enumerate(labels):
    ss = f"{result['SS'][i]:.2f}"
    df = f"{result['df'][i]}"
    ms = f"{result['MS'][i]:.2f}" if result['MS'][i] is not None else ""
    f_val = f"{result['F'][i]:.3f}" if result['F'][i] is not None else ""
    p_val = f"{result['p'][i]:.6f}" if result['p'][i] is not None else ""
    sig = ""
    if result['p'][i] is not None:
        if result['p'][i] < 0.001: sig = " ***"
        elif result['p'][i] < 0.01: sig = " **"
        elif result['p'][i] < 0.05: sig = " *"
    print(f"{label:<8} {ss:>10} {df:>5} {ms:>10} {f_val:>10} {p_val:>10}{sig}")

このANOVA表から、因子A、因子B、交互作用のそれぞれについてp値が計算されています。F値が大きいほど(p値が小さいほど)、その要因の効果が統計的に有意です。

交互作用プロットの可視化

import numpy as np
import matplotlib.pyplot as plt

# 先ほどのANOVA結果を使用
mean_AB = result["means"]["AB"]
a_levels = [f"A{i+1}" for i in range(a)]
b_levels = [f"B{i+1}" for i in range(b)]

fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))

# 交互作用プロット(因子Aを横軸)
ax = axes[0]
x_pos = np.arange(a)
colors = ["#4C72B0", "#DD8452"]
for j in range(b):
    ax.plot(x_pos, mean_AB[:, j], "-o", color=colors[j],
            markersize=8, linewidth=2, label=b_levels[j])
    # 各セルの個別データも表示
    for i in range(a):
        jitter = np.random.normal(0, 0.05, r)
        ax.scatter(np.full(r, x_pos[i]) + jitter, data[i, j, :],
                   color=colors[j], alpha=0.3, s=20)

ax.set_xticks(x_pos)
ax.set_xticklabels(a_levels)
ax.set_xlabel("Factor A level", fontsize=12)
ax.set_ylabel("Response", fontsize=12)
ax.set_title("Interaction plot (Factor A on x-axis)", fontsize=13)
ax.legend(fontsize=11, title="Factor B")
ax.grid(True, alpha=0.3)

# 交互作用プロット(因子Bを横軸)
ax = axes[1]
x_pos = np.arange(b)
colors_a = ["#4C72B0", "#55A868", "#C44E52"]
for i in range(a):
    ax.plot(x_pos, mean_AB[i, :], "-o", color=colors_a[i],
            markersize=8, linewidth=2, label=a_levels[i])
    for j in range(b):
        jitter = np.random.normal(0, 0.05, r)
        ax.scatter(np.full(r, x_pos[j]) + jitter, data[i, j, :],
                   color=colors_a[i], alpha=0.3, s=20)

ax.set_xticks(x_pos)
ax.set_xticklabels(b_levels)
ax.set_xlabel("Factor B level", fontsize=12)
ax.set_ylabel("Response", fontsize=12)
ax.set_title("Interaction plot (Factor B on x-axis)", fontsize=13)
ax.legend(fontsize=11, title="Factor A")
ax.grid(True, alpha=0.3)

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

交互作用プロットから、以下のことが読み取れます。

  1. 左図(因子Aを横軸): 2本の線(B1, B2)が平行でなければ交互作用が存在します。線が交差していれば交互作用は特に強いです。この例では線がやや交差する傾向にあり、A1ではB2が高く、A3ではB1が高い(あるいはその逆)というパターンが見えれば、交互作用の存在を示唆します

  2. 右図(因子Bを横軸): 3本の線(A1, A2, A3)の傾きの違いが交互作用を表しています。すべての線が同じ方向に同じ傾きで変化していれば交互作用はなく、傾きが異なれば交互作用が存在します

  3. 個別データ点(薄い点)はセル内のばらつきを示している。このばらつきが誤差分散 $MS_E$ の推定に使われます。セル間の平均の違いがこのばらつきに比べて十分大きければ、効果は有意です

交互作用の有無による比較

交互作用がある場合とない場合のシミュレーションを比較します。

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

np.random.seed(42)
a, b, r = 3, 2, 8
sigma = 2.0

# --- ケース1: 交互作用なし ---
mu1 = 50
alpha1 = np.array([-4, 0, 4])
beta1 = np.array([-3, 3])
data_no_int = np.zeros((a, b, r))
for i in range(a):
    for j in range(b):
        data_no_int[i, j, :] = mu1 + alpha1[i] + beta1[j] + np.random.randn(r) * sigma

# --- ケース2: 強い交互作用あり ---
interaction2 = np.array([[-4, 4], [0, 0], [4, -4]])
data_with_int = np.zeros((a, b, r))
for i in range(a):
    for j in range(b):
        data_with_int[i, j, :] = (mu1 + alpha1[i] + beta1[j] +
                                   interaction2[i, j] +
                                   np.random.randn(r) * sigma)

fig, axes = plt.subplots(1, 2, figsize=(13, 5.5))

for ax, data_case, title in zip(axes,
                                  [data_no_int, data_with_int],
                                  ["No interaction", "Strong interaction"]):
    mean_AB = np.mean(data_case, axis=2)
    for j in range(b):
        ax.plot(range(a), mean_AB[:, j], "-o", markersize=8, linewidth=2,
                label=f"B{j+1}")
    ax.set_xticks(range(a))
    ax.set_xticklabels([f"A{i+1}" for i in range(a)])
    ax.set_xlabel("Factor A", fontsize=12)
    ax.set_ylabel("Response mean", fontsize=12)

    # ANOVA実行
    result_case = two_way_anova(data_case, a, b, r)
    p_int = result_case["p"][2]
    ax.set_title(f"{title}\n(Interaction p = {p_int:.4f})", fontsize=12)
    ax.legend(fontsize=10, title="Factor B")
    ax.grid(True, alpha=0.3)

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

この比較から、以下のことが視覚的に確認できます。

  1. 交互作用なし(左図): 2本の線がほぼ平行で、因子Bの効果が因子Aの水準によらず一定です。交互作用のp値は大きく、帰無仮説(交互作用なし)は棄却されません

  2. 強い交互作用(右図): 2本の線が交差しており、因子Bの効果が因子Aの水準によって逆転しています。A1ではB2の方が高いのに、A3ではB1の方が高くなっています。交互作用のp値は小さく、帰無仮説は棄却されます

二元配置ANOVAの仮定と診断

前提条件

二元配置ANOVAの妥当性は、以下の仮定に依存します。

  1. 正規性: 各セルのデータが正規分布に従う
  2. 等分散性: すべてのセルの分散が等しい($\sigma_{ij}^2 = \sigma^2$ for all $i, j$)
  3. 独立性: 観測値が互いに独立

正規性の仮定に対して頑健(ロバスト)であることが知られていますが、等分散性の仮定は重要であり、ルビーン検定やバートレット検定で確認すべきです。

効果量

二元配置ANOVAにおける効果量として、偏イータ二乗(partial eta-squared)が広く使われます。

$$ \eta^2_{\text{partial}, A} = \frac{SS_A}{SS_A + SS_E} $$

$$ \eta^2_{\text{partial}, AB} = \frac{SS_{AB}}{SS_{AB} + SS_E} $$

偏イータ二乗は0から1の範囲を取り、その要因が残差分散に対してどれだけの説明力を持つかを表します。

まとめ

本記事では、二元配置分散分析の理論を交互作用の概念から導出し、Pythonで実装しました。

  • 交互作用は「一方の要因の効果が他方の要因の水準によって変わる」現象であり、交互作用プロットで線が非平行になることとして視覚的に確認できます
  • 平方和の分解 $SS_T = SS_A + SS_B + SS_{AB} + SS_E$ により、全変動を各要因と交互作用の寄与に分離します
  • F統計量 $F = MS_{\text{要因}} / MS_E$ は、要因の変動が偶然の変動に比べて大きいかを評価します
  • 検定の解釈では、まず交互作用を確認し、交互作用が有意な場合は主効果だけでなく単純主効果の分析が必要です
  • 繰り返しなし ($r = 1$) の場合は交互作用の検定ができず、加法モデルを仮定する必要があります

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