←Back to the bookChapter 60 · Crop Yield by Fertilizer (F / ANOVA)

Chapter 60 · Case Study

Crop Yield by Fertilizer: F & ANOVA 🌾

A t-test compares two means. To compare FOUR fertilizers at once we need ANOVA, which splits the total variance into between-group and within-group pieces and forms the F-ratio. Then Tukey tells us exactly which fertilizers differ.

Statistics, Data Science and AI: A Visual Handbook · John Fisher · 2026

⚙️ Setup & data

import numpy as np, pandas as pd, matplotlib.pyplot as plt, itertools
from scipy import stats
AMBER="#d97706"; TEAL="#0d9488"; INK="#1a2138"; GRID="#e6e9f2"; PINK="#db2777"
plt.rcParams.update({"figure.facecolor":"white","axes.facecolor":"white","figure.dpi":110,"font.size":11,
   "axes.edgecolor":GRID,"axes.grid":True,"grid.color":GRID,"axes.axisbelow":True,"axes.spines.top":False,
   "axes.spines.right":False,"axes.titlesize":12,"axes.titleweight":"bold","legend.frameon":False})
BASE="https://raw.githubusercontent.com/johnfisher-ai/Statistics-Data-Science-AI-Visual-Book/main/data/"
try:    df = pd.read_csv("../../data/agricultural_yield_optimization.csv")
except FileNotFoundError: df = pd.read_csv(BASE+"agricultural_yield_optimization.csv")
print("loaded:", df.shape)
df.head()

loaded: (1000, 4)

	plot_id	soil_ph	fertilizer_type	crop_yield_bushels
0	PLOT_6000	6.69	Fertilizer_B	49.09
1	PLOT_6001	6.95	Fertilizer_D	50.27
2	PLOT_6002	7.32	Fertilizer_D	49.43
3	PLOT_6003	6.83	Fertilizer_A	51.77
4	PLOT_6004	6.42	Fertilizer_B	49.49

BEAT 1 · FOUR GROUPS, FOUR MEANS

Look at the yields before testing

Each of 1,000 plots got one of four fertilizers and reported its crop_yield_bushels. We group by fertilizer and look at the means and spread, the question is whether the gaps between groups are bigger than the noise within them.

g = "fertilizer_type"; y = "crop_yield_bushels"
summary = df.groupby(g)[y].agg(["count","mean","std"]).round(2)
print(summary)
print(f"\ngrand mean = {df[y].mean():.2f} bushels")

                 count   mean   std
fertilizer_type                    
Fertilizer_A       257  49.62  5.03
Fertilizer_B       250  52.50  5.44
Fertilizer_C       266  48.34  4.72
Fertilizer_D       227  54.86  5.35

grand mean = 51.19 bushels

order = sorted(df[g].unique())
data = [df[df[g]==t][y].values for t in order]
fig,ax=plt.subplots(figsize=(7,3.4))
bp=ax.boxplot(data,tick_labels=[t.replace("Fertilizer_","") for t in order],patch_artist=True,
              boxprops=dict(facecolor="#fde68a"),medianprops=dict(color=PINK,linewidth=2))
ax.axhline(df[y].mean(),color=TEAL,ls="--",lw=1.5,label=f"grand mean {df[y].mean():.1f}")
ax.set_xlabel("fertilizer"); ax.set_ylabel("crop yield (bushels)"); ax.set_title("Yield by fertilizer: D highest, C lowest"); ax.legend()
plt.tight_layout(); plt.show()

The group means range from 48.3 bushels (Fertilizer C) to 54.9 (Fertilizer D), a spread of more than 6 bushels, while each group's own scatter is about 5 bushels. The between-group gaps look larger than the within-group noise, which is exactly what ANOVA quantifies.

BEAT 2 · PARTITION THE VARIANCE

Between-group vs within-group

ANOVA splits the total variation into SSB (how far group means sit from the grand mean) and SSW (scatter inside each group). The F-ratio is MSB/MSW, big when groups differ more than chance allows.

# CHECK ANOVA ASSUMPTIONS before trusting F (independence by design: each plot is separate).
resid = np.concatenate([v - v.mean() for v in data])   # pooled within-group residuals
lev = stats.levene(*data); sh = stats.shapiro(resid); kw = stats.kruskal(*data)
print(f"Levene equal-variance   p = {lev.pvalue:.3f}  -> "
      + ("variances similar; ANOVA valid" if lev.pvalue>0.05 else "variances differ; use Welch ANOVA"))
print(f"Shapiro residual-normal p = {sh.pvalue:.4f}  -> "
      + ("residuals normal" if sh.pvalue>0.05 else "NON-normal (see note); large n + Kruskal rescue it"))
print(f"Kruskal-Wallis cross-check p = {kw.pvalue:.2e} (distribution-free, should agree)")
fig, ax = plt.subplots(figsize=(4.8,3.8)); stats.probplot(resid, dist="norm", plot=ax)
ax.set_title("Residual QQ plot"); plt.tight_layout(); plt.show()

Levene equal-variance   p = 0.134  -> variances similar; ANOVA valid
Shapiro residual-normal p = 0.0053  -> NON-normal (see note); large n + Kruskal rescue it
Kruskal-Wallis cross-check p = 3.65e-42 (distribution-free, should agree)

Before trusting the F-test, check its conditions. Independence holds (each plot is separate) and Levene's test (p = 0.13) confirms the four fertilizers have similar variances. But the residual QQ plot bends at the tails and Shapiro fails (p = 0.005), the residuals are not perfectly normal. In a small study that would be a real worry; here two things rescue us: with n = 1000 the Central Limit Theorem makes the F-test robust, and the distribution-free Kruskal-Wallis test agrees overwhelmingly (p about 4e-42). When non-normality is severe and the sample is small, report Kruskal-Wallis instead; unequal variances would instead call for Welch's ANOVA. This is the discipline: check first, and switch tools when a check fails rather than ignore it.

grand = df[y].mean(); N=len(df); k=df[g].nunique()
groups = {t:df[df[g]==t][y].values for t in order}
SSB = sum(len(v)*(v.mean()-grand)**2 for v in groups.values())
SSW = sum(((v-v.mean())**2).sum() for v in groups.values())
dfb, dfw = k-1, N-k
MSB, MSW = SSB/dfb, SSW/dfw
F_manual = MSB/MSW
print(f"SSB = {SSB:8.1f}   dfb = {dfb:3d}   MSB = {MSB:8.2f}")
print(f"SSW = {SSW:8.1f}   dfw = {dfw:3d}   MSW = {MSW:8.2f}")
print(f"F = MSB/MSW = {F_manual:.2f}")
# scipy one-liner agrees
F, p = stats.f_oneway(*data)
print(f"\nscipy f_oneway: F = {F:.2f},  p-value = {p:.3e}")
print(f"eta-squared (variance explained) = {SSB/(SSB+SSW):.3f}")

SSB =   6264.8   dfb =   3   MSB =  2088.28
SSW =  26201.1   dfw = 996   MSW =    26.31
F = MSB/MSW = 79.38

scipy f_oneway: F = 79.38,  p-value = 4.765e-46
eta-squared (variance explained) = 0.193

The F-statistic is 79.4 on (3, 996) degrees of freedom, with p ≈ 10⁻⁴⁶. Under H₀ (all fertilizers equal) F would hover near 1; a value of 79 is astronomically unlikely by chance, so at least one fertilizer differs. The η² ≈ 0.19 says fertilizer choice explains about 19% of the variation in yield, the rest is plot-to-plot noise.

BEAT 3 · WHICH ONES DIFFER? TUKEY HSD

ANOVA says 'some'; post-hoc says 'which'

A significant F only says the means are not all equal. Tukey's Honest Significant Difference compares every pair while controlling the family-wide error rate, so we can rank the fertilizers safely.

qcrit = stats.studentized_range.ppf(0.95, k, dfw)
print(f"Tukey q critical (k={k}, df={dfw}, alpha=0.05) = {qcrit:.3f}\n")
for a,b in itertools.combinations(order,2):
    na,nb = len(groups[a]),len(groups[b])
    se = np.sqrt(MSW/2*(1/na+1/nb))
    diff = groups[a].mean()-groups[b].mean()
    q = abs(diff)/se
    verdict = "SIGNIFICANT" if q>qcrit else "not sig"
    print(f"{a[-1]} vs {b[-1]}: mean diff = {diff:+5.2f}, q = {q:5.2f}  -> {verdict}")
# Tukey CI plot: each pairwise mean difference with its simultaneous 95% interval
labs=[]; diffs=[]; half=[]
for a,b in itertools.combinations(order,2):
    na,nb=len(groups[a]),len(groups[b]); se=np.sqrt(MSW/2*(1/na+1/nb))
    diffs.append(groups[a].mean()-groups[b].mean()); half.append(qcrit*se); labs.append(f"{a[-1]}-{b[-1]}")
diffs=np.array(diffs); half=np.array(half); ypos=np.arange(len(labs))
fig,ax=plt.subplots(figsize=(6,3.2))
ax.errorbar(diffs,ypos,xerr=half,fmt='o',color=AMBER,capsize=4,lw=1.6)
ax.axvline(0,color=PINK,ls='--'); ax.set_yticks(ypos); ax.set_yticklabels(labs)
ax.set_xlabel('mean yield difference (bushels)')
ax.set_title('Tukey 95% intervals: any interval clear of 0 is a real difference')
plt.tight_layout(); plt.show()

Tukey q critical (k=4, df=996, alpha=0.05) = 3.639

A vs B: mean diff = -2.88, q =  8.92  -> SIGNIFICANT
A vs C: mean diff = +1.28, q =  4.03  -> SIGNIFICANT
A vs D: mean diff = -5.23, q = 15.84  -> SIGNIFICANT
B vs C: mean diff = +4.15, q = 13.00  -> SIGNIFICANT
B vs D: mean diff = -2.36, q =  7.09  -> SIGNIFICANT
C vs D: mean diff = -6.51, q = 19.87  -> SIGNIFICANT

Every one of the six pairwise comparisons clears the Tukey threshold (q > 3.64), so the full ranking is real: D (54.9) > B (52.5) > A (49.6) > C (48.3). Fertilizer D is the clear winner, beating the worst performer C by more than 6 bushels a plot. Tukey keeps this honest: testing all six pairs with naive t-tests would inflate the false-positive rate, while the studentized-range adjustment holds the family-wide error at 5%.

---

✅ Case closed — and Part IX complete

ANOVA compares many means at once by partitioning variance: F = MSB/MSW = 79.4 (p ≈ 10⁻⁴⁶) proves the fertilizers are not all equal, and Tukey HSD ranks them D > B > A > C with every pair significant. This closes our tour of the distributions, from Bernoulli trials to the F-ratio. The same engine, the Central Limit Theorem, underlies the t, chi-square, and F tests alike, and carries us straight into statistical inference.

Statistics, Data Science and AI: A Visual Handbook · © 2026 John Fisher

←Back to the bookChapter 60 · Crop Yield by Fertilizer (F / ANOVA)