import numpy as np
import pandas as pd
import seaborn as sns
import statsmodels.api as sm
from statsmodels.formula.api import ols
from statsmodels.stats.multicomp import pairwise_tukeyhsd, MultiComparison
import scipy.stats as stats
from math import log

# Data from the table for Fe3+ and Fe2+
# Combining data from both tables
fe3_high = [0.71, 1.66, 2.01, 2.16, 2.42, 2.42, 2.56, 2.60, 3.31, 3.64, 3.74, 3.74, 4.39, 4.50, 5.07, 5.26, 8.15, 8.24]
fe3_medium = [2.20, 2.93, 3.08, 3.49, 4.11, 4.95, 5.16, 5.54, 5.68, 6.25, 7.25, 7.90, 8.85, 11.96, 15.54, 15.89, 18.30, 18.59]
fe3_low = [2.25, 3.93, 5.08, 5.82, 5.84, 6.89, 8.50, 8.56, 9.44, 10.52, 13.46, 13.57, 14.76, 16.41, 16.96, 17.56, 22.82, 29.13]

fe2_high = [2.20, 2.69, 3.54, 3.75, 3.83, 4.08, 4.27, 4.53, 5.32, 6.18, 6.22, 6.33, 6.97, 6.97, 7.52, 8.36, 11.65, 12.45]
fe2_medium = [4.04, 4.16, 4.42, 4.93, 5.49, 5.77, 5.86, 6.28, 6.97, 7.06, 7.78, 9.23, 9.34, 9.91, 13.46, 18.40, 23.89, 26.39]
fe2_low = [2.71, 5.43, 6.38, 6.38, 8.32, 9.04, 9.56, 10.01, 10.08, 10.62, 13.80, 15.99, 17.90, 18.25, 19.32, 19.87, 21.60, 22.25]

# Creating a DataFrame for analysis
data = []

for value in fe3_high:
    data.append({'ion_type': 'Fe3+', 'dosage': 'high', 'value': value})
for value in fe3_medium:
    data.append({'ion_type': 'Fe3+', 'dosage': 'medium', 'value': value})
for value in fe3_low:
    data.append({'ion_type': 'Fe3+', 'dosage': 'low', 'value': value})
    
for value in fe2_high:
    data.append({'ion_type': 'Fe2+', 'dosage': 'high', 'value': value})
for value in fe2_medium:
    data.append({'ion_type': 'Fe2+', 'dosage': 'medium', 'value': value})
for value in fe2_low:
    data.append({'ion_type': 'Fe2+', 'dosage': 'low', 'value': value})

df = pd.DataFrame(data)
# Create a combined group identifier for each combination of ion_type and dosage
df['group'] = df['ion_type'] + '_' + df['dosage']

# Apply log transformation to the data
df['log_value'] = df['value'].apply(lambda x: log(x) if x > 0 else 0)

# Basic statistics for each group
print("基础统计信息（按组）:")
stats_by_group = df.groupby(['ion_type', 'dosage'])['value'].agg(['mean', 'std', 'count'])
print(stats_by_group)
print("\n")

# Part 1: Two-way ANOVA analysis
print("第一部分: 双因素方差分析")
model = ols('value ~ C(ion_type) * C(dosage)', data=df).fit()
anova_table = sm.stats.anova_lm(model, typ=2)
print(anova_table)
print("\n")

# Test for homogeneity of variance (Levene's test)
groups = []
for ion in ['Fe3+', 'Fe2+']:
    for dose in ['high', 'medium', 'low']:
        groups.append(df[(df['ion_type'] == ion) & (df['dosage'] == dose)]['value'])

levene_stat, levene_p = stats.levene(*groups)
print(f"Levene方差齐性检验:")
print(f"统计量: {levene_stat:.4f}, p值: {levene_p:.4f}")
print(f"在α = 0.05水平下，方差齐性假设{'不成立' if levene_p < 0.05 else '成立'}\n")

# Part 2: Log transformation and re-analysis
print("第二部分: 对数转换后的分析")
# Two-way ANOVA with log-transformed data
model_log = ols('log_value ~ C(ion_type) * C(dosage)', data=df).fit()
anova_table_log = sm.stats.anova_lm(model_log, typ=2)
print(anova_table_log)

# Levene's test after log transformation
groups_log = []
for ion in ['Fe3+', 'Fe2+']:
    for dose in ['high', 'medium', 'low']:
        groups_log.append(df[(df['ion_type'] == ion) & (df['dosage'] == dose)]['log_value'])

levene_stat_log, levene_p_log = stats.levene(*groups_log)
print(f"\n对数转换后的Levene方差齐性检验:")
print(f"统计量: {levene_stat_log:.4f}, p值: {levene_p_log:.4f}")
print(f"在α = 0.05水平下，方差齐性假设{'不成立' if levene_p_log < 0.05 else '成立'}\n")

# Part 3: Interpreting the results
print("第三部分: 结果解释")
print(f"离子类型的主效应: F={anova_table_log.loc['C(ion_type)', 'F']:.4f}, p={anova_table_log.loc['C(ion_type)', 'PR(>F)']:.4f}")
print(f"剂量水平的主效应: F={anova_table_log.loc['C(dosage)', 'F']:.4f}, p={anova_table_log.loc['C(dosage)', 'PR(>F)']:.4f}")
print(f"交互效应: F={anova_table_log.loc['C(ion_type):C(dosage)', 'F']:.4f}, p={anova_table_log.loc['C(ion_type):C(dosage)', 'PR(>F)']:.4f}")

ion_significant = anova_table_log.loc['C(ion_type)', 'PR(>F)'] < 0.05
dosage_significant = anova_table_log.loc['C(dosage)', 'PR(>F)'] < 0.05
interaction_significant = anova_table_log.loc['C(ion_type):C(dosage)', 'PR(>F)'] < 0.05

print(f"\n在α = 0.05水平下:")
print(f"- 离子类型主效应{'显著' if ion_significant else '不显著'}")
print(f"- 剂量水平主效应{'显著' if dosage_significant else '不显著'}")
print(f"- 交互效应{'显著' if interaction_significant else '不显著'}\n")

# Part 4: Calculating confidence intervals for all six treatment combinations
print("第四部分: 所有六个处理组合的置信区间及Bonferroni校正")

# Calculate means and 95% CIs for each of the six combinations
print("\n各处理组合的均值和95%置信区间 (对数转换后):")
combination_stats = {}

for ion in ['Fe3+', 'Fe2+']:
    for dosage in ['high', 'medium', 'low']:
        group_data = df[(df['ion_type'] == ion) & (df['dosage'] == dosage)]['log_value']
        mean = group_data.mean()
        std = group_data.std()
        n = len(group_data)
        se = std / np.sqrt(n)
        t_crit = stats.t.ppf(0.975, n-1)  # 95% CI
        ci_lower = mean - t_crit * se
        ci_upper = mean + t_crit * se
        
        combination_name = f"{ion}_{dosage}"
        combination_stats[combination_name] = {
            'mean': mean,
            'std': std,
            'n': n,
            'se': se,
            'ci_lower': ci_lower,
            'ci_upper': ci_upper
        }
        
        print(f"{ion}, {dosage}剂量: 均值 = {mean:.4f}, 95% CI = [{ci_lower:.4f}, {ci_upper:.4f}]")

# Compare ion types at each dosage level (with Bonferroni correction)
print("\n同一剂量水平下不同离子类型的均值差异 (Bonferroni校正):")

# Number of comparisons for Bonferroni correction
n_comparisons = 3  # Comparing Fe3+ vs Fe2+ at 3 dosage levels
alpha_bonf = 0.05 / n_comparisons

for dosage in ['high', 'medium', 'low']:
    fe3_data = df[(df['ion_type'] == 'Fe3+') & (df['dosage'] == dosage)]['log_value']
    fe2_data = df[(df['ion_type'] == 'Fe2+') & (df['dosage'] == dosage)]['log_value']
    
    # Calculate mean difference
    mean_diff = fe2_data.mean() - fe3_data.mean()
    
    # Calculate standard error of the difference
    n1 = len(fe3_data)
    n2 = len(fe2_data)
    s1 = fe3_data.std()
    s2 = fe2_data.std()
    se_diff = np.sqrt((s1**2 / n1) + (s2**2 / n2))
    
    # Calculate degrees of freedom (Welch-Satterthwaite approximation)
    df_welch = (s1**2/n1 + s2**2/n2)**2 / ((s1**2/n1)**2/(n1-1) + (s2**2/n2)**2/(n2-1))
    
    # Calculate t-statistic and p-value
    t_stat = mean_diff / se_diff
    p_val = 2 * (1 - stats.t.cdf(abs(t_stat), df_welch))
    p_val_corr = min(p_val * n_comparisons, 1)  # Bonferroni correction
    
    # Calculate Bonferroni-corrected confidence interval
    t_crit_bonf = stats.t.ppf(1 - alpha_bonf/2, df_welch)
    ci_lower = mean_diff - t_crit_bonf * se_diff
    ci_upper = mean_diff + t_crit_bonf * se_diff
    
    print(f"\n{dosage}剂量下的Fe2+与Fe3+比较:")
    print(f"均值差 (Fe2+ - Fe3+) = {mean_diff:.4f}")
    print(f"t统计量 = {t_stat:.4f}, 原始p值 = {p_val:.4f}, Bonferroni校正后p值 = {p_val_corr:.4f}")
    print(f"在α = 0.05水平下，差异{'显著' if p_val_corr < 0.05 else '不显著'}")
    print(f"95% Bonferroni同时置信区间: [{ci_lower:.4f}, {ci_upper:.4f}]")
    print(f"解释: 由于置信区间{'不' if 0 < ci_lower or ci_upper < 0 else ''}包含0，所以在{dosage}剂量水平下，两种离子类型的差异在统计学上{'显著' if 0 < ci_lower or ci_upper < 0 else '不显著'}")

# Also do overall comparison between ion types and between dosage levels
print("\n总体离子类型比较（对数转换后）:")
mc_ion = MultiComparison(df['log_value'], df['ion_type'])
result_ion = mc_ion.allpairtest(stats.ttest_ind, method='bonf')
print(result_ion[0])

print("\n总体剂量水平比较（对数转换后）:")
mc_dosage = MultiComparison(df['log_value'], df['dosage'])
result_dosage = mc_dosage.allpairtest(stats.ttest_ind, method='bonf')
print(result_dosage[0])

print("\n分析总结:")
print("1. 基于方差分析结果，铁离子类型和剂量水平对铁离子在小白鼠体内的存留量均有显著影响")
print("2. 在相同剂量水平下，Fe2+在小白鼠体内的存留量总体高于Fe3+")
print("3. 随着剂量水平从高到低，两种铁离子在小白鼠体内的存留量均显著增加")
print("4. Bonferroni校正的置信区间分析表明，在不同剂量水平下，两种铁离子的差异程度不同")