# =============================================================
# ANCOVA — ANÁLISIS DE COVARIANZA
# Efecto del CaCl2 sobre la firmeza poscosecha del tomate de árbol
# (Solanum betaceum Cav.) — día 20
# Covariable: peso del fruto (g)
# =============================================================

import numpy as np
import pandas as pd
from scipy import stats
from scipy.stats import studentized_range
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm
from itertools import combinations
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

# ─────────────────────────────────────────────────────────────
# 1. DATOS
# ─────────────────────────────────────────────────────────────
firmeza = {
    'Control (0 mM)': [32.2352,29.3778,32.9146,36.8536,28.9463,28.9464,37.1065,33.4535,27.8874,32.4415,27.9146,27.9042,31.0888,21.3902,22.2379],
    '570 mM CaCl2': [47.3010,45.1384,51.5084,45.6415,43.2209,57.0351,48.9163,50.3241,43.1612,47.3870,50.5324,44.4752,51.8034,47.1169,48.5999],
    '862 mM CaCl2': [49.2321,60.5205,51.9379,47.1345,55.7837,46.3841,52.9608,42.9855,45.8903,52.9056,55.3969,52.7883,51.4680,50.6149,45.1988]
}

peso = {
    'Control (0 mM)': [90.00,85.86,86.47,86.75,85.21,86.95,85.81,82.46,89.33,87.78,85.71,86.71,87.76,87.42,87.34],
    '570 mM CaCl2': [81.20,85.92,84.09,85.24,81.61,86.68,84.52,83.13,89.44,84.29,80.76,83.90,80.00,86.74,83.31],
    '862 mM CaCl2': [82.50,84.93,80.13,85.60,80.00,83.07,80.97,88.22,88.48,82.90,85.13,83.25,85.07,82.29,80.93]
}

trts = list(firmeza.keys())

records = []
for t in trts:
    for f, p in zip(firmeza[t], peso[t]):
        records.append({'tratamiento': t, 'firmeza': f, 'peso': p})
df = pd.DataFrame(records)

df['trt'] = pd.Categorical(df['tratamiento'], categories=trts, ordered=False)
REF = 'C(trt, Treatment(\'Control (0 mM)\'))'

print('='*65)
print('ANCOVA — DATOS COMPLETOS (N = 45)')
print('='*65)
print(df[['tratamiento','firmeza','peso']].to_string(index=False))

print('\n' + '='*65)
print('ESTADÍSTICAS DESCRIPTIVAS')
print('='*65)

for lbl, col in [('Firmeza (% peso)','firmeza'),('Peso (g)','peso')]:
    print(f'\n{lbl}:')
    d = df.groupby('tratamiento')[col].agg(n='count',media='mean',DE=lambda x:x.std(ddof=1),minimo='min',maximo='max').round(4)
    print(d.loc[trts].to_string())

r_glob, p_glob = stats.pearsonr(df['firmeza'], df['peso'])
print(f'\nCorrelación Pearson firmeza-peso: r={r_glob:.4f}, p={p_glob:.2e}')

m_add = ols(f'firmeza ~ {REF} + peso', data=df).fit()
aov = anova_lm(m_add)
params = m_add.params

print('\n' + '='*65)
print('TABLA ANCOVA')
print('='*65)
print(aov)

print(f'Pendiente beta = {params["peso"]:.4f}')

peso_bar = df['peso'].mean()
intercept = params['Intercept']

print('\nMedias ajustadas:')
for t in trts:
    print(f'{t}: {intercept + params["peso"]*peso_bar:.4f}')

plt.figure(figsize=(6,4))
plt.scatter(df['peso'], df['firmeza'])
plt.xlabel('Peso')
plt.ylabel('Firmeza')
plt.title('Relación peso-firmeza')
plt.grid(True)
plt.show()