Agricultural Crop Yield ANOVA

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

np.random.seed(1313)

# Yield in tons per hectare
control = np.random.normal(3.2, 0.6, 80)
organic = np.random.normal(3.8, 0.55, 80)
synthetic = np.random.normal(4.5, 0.7, 80)
bio_enhanced = np.random.normal(4.2, 0.5, 80)

F_stat, p_value = stats.f_oneway(control, organic, synthetic, bio_enhanced)

fig, ax = plt.subplots(figsize=(12, 7), facecolor='#ffffff')
ax.set_facecolor='#ffffff'

colors = ['#9C2007', '#6CF527', '#276CF5', '#F5D327']

parts = ax.violinplot([control, organic, synthetic, bio_enhanced], 
                       positions=[1, 2, 3, 4], showmeans=True, showmedians=True, widths=0.7)

for i, pc in enumerate(parts['bodies']):
    pc.set_facecolor(colors[i])
    pc.set_alpha(0.65)
    pc.set_edgecolor(colors[i])
    pc.set_linewidth(2)

parts['cmeans'].set_color='#F5276C'
parts['cmeans'].set_linewidth(2.5)
parts['cmedians'].set_color('#1f2937')
for partname in ['cbars', 'cmins', 'cmaxes']:
    parts[partname].set_color('#9ca3af')

# Baseline reference
ax.axhline(y=3.0, color='#ef4444', linestyle=':', alpha=0.6, linewidth=1.5)
ax.text(4.55, 3.0, 'Baseline', fontsize=9, color='#ef4444', va='center')

labels = ['Control\n(No Fertilizer)', 'Organic', 'Synthetic\nNPK', 'Bio-Enhanced']
means = [control.mean(), organic.mean(), synthetic.mean(), bio_enhanced.mean()]

# Improvement percentages
baseline = control.mean()
for i, (mean, color) in enumerate(zip(means[1:], colors[1:])):
    improvement = ((mean - baseline) / baseline) * 100
    ax.text(i+2, mean + 0.4, f'+{improvement:.0f}%', ha='center', fontsize=10, 
            color=color, fontweight='bold')

# Cost annotations
costs = ['$0', '$180', '$120', '$250']
for i, (cost, color) in enumerate(zip(costs, colors)):
    ax.text(i+1, 1.5, f'Cost/ha: {cost}', ha='center', fontsize=8, color='#6b7280')

stats_text = f"YIELD ANALYSIS\n{'─'*16}\nF = {F_stat:.2f}\np = {p_value:.4f}\n\nHighest Yield:\nSynthetic NPK\n(μ = {synthetic.mean():.2f} t/ha)"
bbox = dict(boxstyle="round,pad=0.5", facecolor='#fefce8', edgecolor='#F5D327', lw=2)
ax.text(0.02, 0.98, stats_text, transform=ax.transAxes, fontsize=10, color='#1f2937',
        va='top', fontfamily='monospace', bbox=bbox)

ax.set_xticks([1, 2, 3, 4])
ax.set_xticklabels(labels, fontsize=10, color='#1f2937')
ax.set_ylabel('Yield (Tons/Hectare)', fontsize=12, color='#1f2937', fontweight='500')
ax.set_title('Wheat Yield by Fertilizer Treatment\n2-Year Field Trial Results', 
             fontsize=14, color='#1f2937', fontweight='bold', pad=15)

ax.tick_params(colors='#374151')
for spine in ax.spines.values():
    spine.set_color('#e5e7eb')
ax.yaxis.grid(True, color='#f3f4f6', linewidth=0.8)
ax.set_axisbelow(True)
ax.set_ylim(1.3, 6.5)

plt.tight_layout()
plt.show()