import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import seaborn as sns

# 设置中文显示
plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False    # 用来正常显示负号

# 数据表格
data = np.array([
    # 成分A低剂量
    [2.4, 2.7, 2.3, 2.5,    4.6, 4.2, 4.9, 4.7,    4.8, 4.5, 4.4, 4.6],
    # 成分A中剂量
    [5.8, 5.2, 5.5, 5.3,    8.9, 9.1, 8.7, 9.0,    9.1, 9.3, 8.7, 9.4],
    # 成分A高剂量
    [6.1, 5.7, 5.9, 6.2,    9.9, 10.5, 10.6, 10.1,  13.5, 13.0, 13.3, 13.2]
])

# 将数据重塑为适合分析的格式
times = []  # 缓解时间
factor_A = []  # A的剂量水平
factor_B = []  # B的剂量水平

# 剂量水平的文本表示
A_levels = ['低', '中', '高']
B_levels = ['低', '中', '高']

for i in range(3):  # A的三个水平
    for j in range(3):  # B的三个水平
        for k in range(4):  # 每组4个样本
            index = j * 4 + k
            times.append(data[i][index])
            factor_A.append(A_levels[i])
            factor_B.append(B_levels[j])

# 创建数据框
df = pd.DataFrame({
    '缓解时间': times,
    '成分A': factor_A,
    '成分B': factor_B
})

# 计算每个组合的均值
means = df.groupby(['成分A', '成分B'])['缓解时间'].mean().unstack()
print("各组均值:")
print(means)

# 将均值矩阵转为数组，方便后续处理
mean_array = means.values

print("\n" + "="*50)
print("第(1)题：计算每个水平组合(A_i, B_j)上的均值μ_ij的估计值y_ij，做出形如图3.2的图形")
print("="*50)

# 计算各水平组合的均值估计值
print("\n各水平组合(A_i, B_j)的均值估计值y_ij:")
for i in range(3):
    for j in range(3):
        cell_mean = mean_array[i, j]
        print(f"A_{i+1}B_{j+1} (A={A_levels[i]}, B={B_levels[j]}): {cell_mean:.2f}")

# 绘制交互作用图
plt.figure(figsize=(14, 6))

# 方法1: 以A为x轴，不同B为线条
plt.subplot(1, 2, 1)
for j, b_level in enumerate(B_levels):
    plt.plot(A_levels, mean_array[:, j], marker='o', label=f'B-{b_level}剂量')
plt.xlabel('成分A剂量')
plt.ylabel('平均缓解时间(h)')
plt.title('交互作用图 (A为x轴)')
plt.legend()
plt.grid(True)

# 方法2: 以B为x轴，不同A为线条
plt.subplot(1, 2, 2)
for i, a_level in enumerate(A_levels):
    plt.plot(B_levels, mean_array[i, :], marker='o', label=f'A-{a_level}剂量')
plt.xlabel('成分B剂量')
plt.ylabel('平均缓解时间(h)')
plt.title('交互作用图 (B为x轴)')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.savefig('交互作用图.png')
plt.show()

# 根据图形判断交互效应
print("\n根据交互作用图的观察:")
print("从图形可以看出，成分A和成分B的剂量水平线不平行，这表明可能存在交互效应。")
print("特别是在A为高剂量时，随着B剂量的增加，缓解时间增长更为明显，这暗示交互效应可能是显著的。")

print("\n" + "="*50)
print("第(2)题：假设所给数据服从方差分析模型，建立方差分析表，A与B的交互效应在显著性水平α=0.05下是否显著？")
print("="*50)

# 进行双因素方差分析
model = ols('缓解时间 ~ C(成分A) + C(成分B) + C(成分A):C(成分B)', data=df).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print("\n方差分析结果:")
print(anova_table)

# 判断交互效应是否显著 (α=0.05)
p_interaction = anova_table.loc['C(成分A):C(成分B)', 'PR(>F)']
print(f"\n交互效应的p值: {p_interaction:.4f}")
if p_interaction < 0.05:
    conclusion = "在显著性水平α=0.05下，交互效应显著"
else:
    conclusion = "在显著性水平α=0.05下，交互效应不显著"
print(f"结论: {conclusion}")

print("\n" + "="*50)
print("第(3)题：若A与B的交互效应显著，分别就A的各水平给出在B的各水平上的均值的置信度为95%的置信区间以及两两之差的置信度不小于95%的Bonferroni同时置信区间")
print("="*50)

# 判断是否需要进行第(3)题的分析
if p_interaction < 0.05:
    # 对B的每个水平，分析A的简单主效应
    print("\n1. 固定B，比较不同水平的A (简单主效应分析):")
    
    for j, b_level in enumerate(B_levels):
        print(f"\n在B={b_level}剂量水平下:")
        
        # 提取当前B水平下的数据
        subset = df[df['成分B'] == b_level]
        
        # 进行单因素方差分析
        model_simple = ols('缓解时间 ~ C(成分A)', data=subset).fit()
        anova_simple = sm.stats.anova_lm(model_simple, typ=1)
        
        print(f"A的简单主效应方差分析:")
        print(anova_simple)
        
        # 计算A在此B水平下各水平的均值和95%置信区间
        print(f"\nA在B={b_level}水平下各水平的均值和95%置信区间:")
        for i, a_level in enumerate(A_levels):
            group = subset[subset['成分A'] == a_level]['缓解时间']
            n = len(group)
            mean = group.mean()
            std = group.std()
            sem = std / np.sqrt(n)
            df_t = n - 1
            t_crit = stats.t.ppf(0.975, df_t)  # 95% CI
            margin = t_crit * sem
            
            print(f"A={a_level}: 均值={mean:.3f}, 95% CI=[{mean-margin:.3f}, {mean+margin:.3f}]")
        
        # 如果A的简单主效应显著，进行多重比较
        if anova_simple.loc['C(成分A)', 'PR(>F)'] < 0.05:
            print(f"\nA的效应在B={b_level}水平下显著，进行多重比较:")
            
            # Tukey HSD多重比较
            tukey = pairwise_tukeyhsd(endog=subset['缓解时间'], groups=subset['成分A'], alpha=0.05)
            print(tukey)
            
            # Bonferroni校正的置信区间
            print("\nBonferroni同时置信区间:")
            a_levels = sorted(subset['成分A'].unique())
            for i in range(len(a_levels)):
                for j in range(i+1, len(a_levels)):
                    group1 = subset[subset['成分A'] == a_levels[i]]['缓解时间']
                    group2 = subset[subset['成分A'] == a_levels[j]]['缓解时间']
                    
                    # 计算差异的95%置信区间
                    n1, n2 = len(group1), len(group2)
                    mean1, mean2 = group1.mean(), group2.mean()
                    var1, var2 = group1.var(), group2.var()
                    
                    # 合并方差
                    pooled_var = ((n1-1)*var1 + (n2-1)*var2) / (n1+n2-2)
                    
                    # 标准误差
                    se = np.sqrt(pooled_var * (1/n1 + 1/n2))
                    
                    # Bonferroni校正的临界值 (k=3 for 3 pairs)
                    df_t = n1 + n2 - 2
                    alpha_bonf = 0.05 / 3  # 校正后的α值
                    t_crit = stats.t.ppf(1 - alpha_bonf/2, df_t)
                    
                    # 计算置信区间
                    margin = t_crit * se
                    ci_lower = (mean1 - mean2) - margin
                    ci_upper = (mean1 - mean2) + margin
                    
                    print(f"A={a_levels[i]}与A={a_levels[j]}的差异: {mean1-mean2:.3f}, 95% Bonferroni CI: [{ci_lower:.3f}, {ci_upper:.3f}]")
    
    # 对A的每个水平，分析B的简单主效应和置信区间
    print("\n2. 固定A，比较不同水平的B (简单主效应分析):")
    
    for i, a_level in enumerate(A_levels):
        print(f"\n在A={a_level}剂量水平下:")
        
        # 提取当前A水平下的数据
        subset = df[df['成分A'] == a_level]
        
        # 进行单因素方差分析
        model_simple = ols('缓解时间 ~ C(成分B)', data=subset).fit()
        anova_simple = sm.stats.anova_lm(model_simple, typ=1)
        
        print(f"B的简单主效应方差分析:")
        print(anova_simple)
        
        # 计算B在此A水平下各水平的均值和95%置信区间
        print(f"\nB在A={a_level}水平下各水平的均值和95%置信区间:")
        for j, b_level in enumerate(B_levels):
            group = subset[subset['成分B'] == b_level]['缓解时间']
            n = len(group)
            mean = group.mean()
            std = group.std()
            sem = std / np.sqrt(n)
            df_t = n - 1
            t_crit = stats.t.ppf(0.975, df_t)  # 95% CI
            margin = t_crit * sem
            
            print(f"B={b_level}: 均值={mean:.3f}, 95% CI=[{mean-margin:.3f}, {mean+margin:.3f}]")
        
        # 如果B的简单主效应显著，进行多重比较
        if anova_simple.loc['C(成分B)', 'PR(>F)'] < 0.05:
            print(f"\nB的效应在A={a_level}水平下显著，进行多重比较:")
            
            # Tukey HSD多重比较
            tukey = pairwise_tukeyhsd(endog=subset['缓解时间'], groups=subset['成分B'], alpha=0.05)
            print(tukey)
            
            # Bonferroni校正的置信区间
            print("\nBonferroni同时置信区间:")
            b_levels = sorted(subset['成分B'].unique())
            for i in range(len(b_levels)):
                for j in range(i+1, len(b_levels)):
                    group1 = subset[subset['成分B'] == b_levels[i]]['缓解时间']
                    group2 = subset[subset['成分B'] == b_levels[j]]['缓解时间']
                    
                    # 计算差异的95%置信区间
                    n1, n2 = len(group1), len(group2)
                    mean1, mean2 = group1.mean(), group2.mean()
                    var1, var2 = group1.var(), group2.var()
                    
                    # 合并方差
                    pooled_var = ((n1-1)*var1 + (n2-1)*var2) / (n1+n2-2)
                    
                    # 标准误差
                    se = np.sqrt(pooled_var * (1/n1 + 1/n2))
                    
                    # Bonferroni校正的临界值 (k=3 for 3 pairs)
                    df_t = n1 + n2 - 2
                    alpha_bonf = 0.05 / 3  # 校正后的α值
                    t_crit = stats.t.ppf(1 - alpha_bonf/2, df_t)
                    
                    # 计算置信区间
                    margin = t_crit * se
                    ci_lower = (mean1 - mean2) - margin
                    ci_upper = (mean1 - mean2) + margin
                    
                    print(f"B={b_levels[i]}与B={b_levels[j]}的差异: {mean1-mean2:.3f}, 95% Bonferroni CI: [{ci_lower:.3f}, {ci_upper:.3f}]")
    
    # 第三小题结束部分，选出最佳的剂量组合
    print("\n根据分析结果推荐最佳剂量组合:")
    best_mean = df.groupby(['成分A', '成分B'])['缓解时间'].mean().max()
    best_combo = df.groupby(['成分A', '成分B'])['缓解时间'].mean().idxmax()
    print(f"最佳组合为: 成分A={best_combo[0]}, 成分B={best_combo[1]}, 平均缓解时间: {best_mean:.2f}h")
    
    # 寻找最佳水平组合
    print("\n寻找最佳剂量组合（假设缓解时间越长越好）:")
    for a_level in A_levels:
        subset_a = df[df['成分A'] == a_level]
        for b_level in B_levels:
            subset_ab = subset_a[subset_a['成分B'] == b_level]
            mean_ab = subset_ab['缓解时间'].mean()
            print(f"A={a_level}, B={b_level}: 平均缓解时间 = {mean_ab:.2f}h")
else:
    print("\n由于交互效应不显著，不需要进行简单主效应分析。")
    print("可以直接分析A和B的主效应。")
    
    # 分析A的主效应
    print("\nA的主效应:")
    a_means = df.groupby('成分A')['缓解时间'].mean()
    print(a_means)
    
    # 分析B的主效应
    print("\nB的主效应:")
    b_means = df.groupby('成分B')['缓解时间'].mean()
    print(b_means)