二元配置分散分析の理論と実装 — 交互作用の解釈

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

np.random.seed(42)

# 二元配置のデータ生成
# 要因A: 肥料（3水準）、要因B: 日照（2水準）
a_levels = ['Fertilizer A', 'Fertilizer B', 'Fertilizer C']
b_levels = ['Low light', 'High light']
n_rep = 10  # 各セルの繰り返し数

# セル平均（交互作用あり）
cell_means = {
    ('Fertilizer A', 'Low light'): 85,
    ('Fertilizer A', 'High light'): 95,
    ('Fertilizer B', 'Low light'): 100,
    ('Fertilizer B', 'High light'): 105,
    ('Fertilizer C', 'Low light'): 90,
    ('Fertilizer C', 'High light'): 110,  # Cは高日照で特に効果的
}
sigma = 6

# データ生成
data = {}
all_values = []
all_a = []
all_b = []

for a in a_levels:
    for b in b_levels:
        values = np.random.normal(cell_means[(a, b)], sigma, n_rep)
        data[(a, b)] = values
        all_values.extend(values)
        all_a.extend([a] * n_rep)
        all_b.extend([b] * n_rep)

all_values = np.array(all_values)

# 手動で二元配置ANOVA
grand_mean = np.mean(all_values)
N = len(all_values)
a, b, n = len(a_levels), len(b_levels), n_rep

# 各平均
a_means = {level: np.mean([data[(level, bl)] for bl in b_levels]) for level in a_levels}
b_means = {level: np.mean([data[(al, level)] for al in a_levels]) for level in b_levels}
cell_mean_obs = {(al, bl): np.mean(data[(al, bl)]) for al in a_levels for bl in b_levels}

# 平方和
SS_A = b * n * sum((a_means[al] - grand_mean)**2 for al in a_levels)
SS_B = a * n * sum((b_means[bl] - grand_mean)**2 for bl in b_levels)
SS_AB = n * sum((cell_mean_obs[(al,bl)] - a_means[al] - b_means[bl] + grand_mean)**2
                for al in a_levels for bl in b_levels)
SS_E = sum(np.sum((data[(al,bl)] - cell_mean_obs[(al,bl)])**2)
           for al in a_levels for bl in b_levels)
SST = np.sum((all_values - grand_mean)**2)

# 自由度
df_A = a - 1
df_B = b - 1
df_AB = (a-1)*(b-1)
df_E = a*b*(n-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_A = MS_A / MS_E
F_B = MS_B / MS_E
F_AB = MS_AB / MS_E

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)

print("Two-way ANOVA Table:")
print(f"{'Source':<15s} {'SS':>10s} {'df':>5s} {'MS':>10s} {'F':>8s} {'p':>8s}")
print("-" * 60)
print(f"{'Factor A':<15s} {SS_A:10.2f} {df_A:5d} {MS_A:10.2f} {F_A:8.3f} {p_A:8.4f}")
print(f"{'Factor B':<15s} {SS_B:10.2f} {df_B:5d} {MS_B:10.2f} {F_B:8.3f} {p_B:8.4f}")
print(f"{'A x B':<15s} {SS_AB:10.2f} {df_AB:5d} {MS_AB:10.2f} {F_AB:8.3f} {p_AB:8.4f}")
print(f"{'Error':<15s} {SS_E:10.2f} {df_E:5d} {MS_E:10.2f}")
print(f"{'Total':<15s} {SST:10.2f} {N-1:5d}")

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

# (a) 交互作用プロット
ax = axes[0]
x_pos = np.arange(len(a_levels))
for j, bl in enumerate(b_levels):
    means = [cell_mean_obs[(al, bl)] for al in a_levels]
    sems = [np.std(data[(al, bl)], ddof=1) / np.sqrt(n) for al in a_levels]
    color = ['#2196F3', '#FF9800'][j]
    ax.errorbar(x_pos, means, yerr=sems, marker='o', linewidth=2,
                markersize=8, capsize=5, color=color, label=bl)

ax.set_xticks(x_pos)
ax.set_xticklabels(a_levels, fontsize=10)
ax.set_ylabel('Yield (kg)', fontsize=12)
ax.set_title('Interaction Plot', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

# (b) 変動の分解
ax = axes[1]
ss_values = [SS_A, SS_B, SS_AB, SS_E]
ss_labels = ['$SS_A$\n(Fertilizer)', '$SS_B$\n(Light)', '$SS_{AB}$\n(Interaction)', '$SS_E$\n(Error)']
colors_bar = ['#4CAF50', '#2196F3', '#FF5722', '#9E9E9E']

bars = ax.bar(ss_labels, ss_values, color=colors_bar, alpha=0.7, edgecolor='gray')
for bar, val, p in zip(bars, ss_values, [p_A, p_B, p_AB, None]):
    sig = ''
    if p is not None:
        sig = ' ***' if p < 0.001 else ' **' if p < 0.01 else ' *' if p < 0.05 else ' n.s.'
    ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 30,
            f'{val:.0f}{sig}', ha='center', fontsize=9)

ax.set_ylabel('Sum of Squares', fontsize=12)
ax.set_title('Decomposition of Variance', fontsize=13)
ax.grid(True, alpha=0.3, axis='y')

# (c) 交互作用ありvs なしの比較
ax = axes[2]
# 交互作用なしのデータ
cell_means_no_int = {
    ('A', 'Low'): 85, ('A', 'High'): 95,
    ('B', 'Low'): 95, ('B', 'High'): 105,
    ('C', 'Low'): 90, ('C', 'High'): 100,
}
# 交互作用ありのデータ（元のまま）
cell_means_int = {
    ('A', 'Low'): 85, ('A', 'High'): 95,
    ('B', 'Low'): 100, ('B', 'High'): 105,
    ('C', 'Low'): 90, ('C', 'High'): 110,
}

for data_dict, title, style in [(cell_means_no_int, 'No interaction', '--'),
                                  (cell_means_int, 'With interaction', '-')]:
    for j, (bl, color) in enumerate(zip(['Low', 'High'], ['#2196F3', '#FF9800'])):
        means = [data_dict[(al, bl)] for al in ['A', 'B', 'C']]
        ax.plot([0, 1, 2], means, style, color=color, linewidth=2, alpha=0.7,
                marker='o', markersize=6,
                label=f'{bl} light ({title})' if j == 0 or style == '-' else None)

ax.set_xticks([0, 1, 2])
ax.set_xticklabels(['A', 'B', 'C'])
ax.set_xlabel('Fertilizer', fontsize=12)
ax.set_ylabel('Mean Yield', fontsize=12)
ax.set_title('Parallel (no int.) vs Non-parallel (int.)', fontsize=12)
ax.legend(fontsize=7, loc='upper left')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('anova_two_way.png', dpi=150, bbox_inches='tight')
plt.show()