How to analyse an on-farm strip trial in Python

Problem statement

A strip trial is the only way to measure how a field responds to a rate change, and the statistics are harder than they look. Yield varies across a field by more than most treatment effects, so a treated strip that happens to lie on the better land shows a response that is the land, not the treatment.

Two things fix it. Pairing each treated strip with its immediate neighbours removes most of the spatial trend, because adjacent strips share their soil. And a test that accounts for spatial autocorrelation gives an honest standard error, because yield points a few metres apart are not independent observations — treating 20,000 of them as independent makes almost any difference significant.

Quick answer

Compare each treated strip with its neighbours, then test on strip means:

import numpy as np, pandas as pd
from scipy import stats

season = strips[strips.year == 1999]
centre = season.set_index("strip_id").geometry.centroid
treated_strips = season.loc[season.treatment == "N5", "strip_id"]
controls = season.loc[season.treatment == "N0", "strip_id"]
neighbours = {s: list(centre[controls].distance(centre[s]).nsmallest(2).index) for s in treated_strips}

paired = []
for s in treated_strips:
    t = yields[yields.strip_id == s]["yield_t_ha"].mean()
    c = yields[yields.strip_id.isin(neighbours[s])]["yield_t_ha"].mean()
    paired.append({"strip": s, "treated": t, "control": c, "diff": t - c})

p = pd.DataFrame(paired)
res = stats.ttest_rel(p["treated"], p["control"])
print(f"{len(p)} pairs; mean difference {p['diff'].mean():+.3f} t/ha "
      f"(sd {p['diff'].std():.3f}), p = {res.pvalue:.3f}")

3 pairs; mean difference +0.948 t/ha (sd 0.239), p = 0.021

The degrees of freedom are the number of strip pairs, not the number of yield points. Three pairs is three observations however many thousand data points the combine recorded — here the 1999 season of the Las Rosas nitrogen trial, with the highest rate (125–131 kg N/ha) against none.

Step-by-step solution

1. Lay the trial out so it can be analysed

Strips the width of the harvester (or a multiple), running across the dominant gradient, alternating or replicated, with at least six to eight replicates. A design with two strips cannot be analysed however good the statistics.

2. Clean the yield data first

Lag correction, start and stop trimming, partial-width removal. The sensor lag smears yield along the pass by tens of metres, so the start and end of every strip carry the previous treatment — see Yield monitor data explained and the cleaning guide.

3. Trim the strip ends — and the sides if the strips are wide

The lag and the rate transition contaminate the first and last tens of metres along each strip; pass and spreader overlap contaminate its sides. Trim both where the geometry allows. Buffering inwards only works on strips several passes wide: the Las Rosas strips are one harvester pass, 9.8 m, and a 5 m inward buffer would delete them entirely, so there the trim is along the strip only.

4. Aggregate to strip means before testing

Each strip is one experimental unit. Testing on individual yield points treats hundreds or thousands of correlated observations as independent and produces a p-value that is meaningless:

from scipy import stats

for year in (1999, 2001):
    season = yields[yields.year == year]
    raw = stats.ttest_ind(season.loc[season.treatment == "N5", "yield_t_ha"],
                          season.loc[season.treatment == "N0", "yield_t_ha"])
    means = season.groupby(["rep", "treatment"])["yield_t_ha"].mean().unstack()
    paired = stats.ttest_rel(means["N5"], means["N0"])
    print(f"{year}: {len(season[season.treatment.isin(['N5', 'N0'])])} points, p = {raw.pvalue:.1e} "
          f"| 3 strip pairs, p = {paired.pvalue:.3f}")

1999: 582 points, p = 1.6e-35 | 3 strip pairs, p = 0.030
2001: 566 points, p = 5.2e-03 | 3 strip pairs, p = 0.019

In 1999 the points claim a certainty of 1 in 10³⁵; the three strip pairs support a p-value of 0.03. Both years' responses are real, but only the second number is an honest measure of how sure the trial can be.

5. Pair with neighbours rather than pooling

The mean of all treated strips against the mean of all control strips is vulnerable to any gradient across the field. Comparing each treated strip with the mean of its immediate neighbours removes most of it.

6. Check the residual spatial structure

If the differences still show a spatial pattern, the pairing has not removed the trend and a spatial model — a linear mixed model with a spatial correlation structure, or simply including position as a covariate — is needed.

7. Report the effect size and the interval, not just the p-value

A difference with its 95% interval is a result even when the interval is wide, and the difference the trial could have detected tells the farmer whether it was worth running. With three pairs per season, the Las Rosas trial could only have detected a response of 0.5–0.7 t/ha (Example 2). That is more useful than "significant" or "not significant".

Code examples

The data used on this page

The Las Rosas trial: yield-monitor data from a maize field near Río Cuarto, Córdoba, Argentina, harvested in 1999 and 2001, with six nitrogen rates (N0, none, to N5, 125–131 kg/ha) applied in strips and three replicates — 3,443 points, averaged along each strip to 9.8 m spacing, which is also the strip width. It was published by Bongiovanni and Lowenberg-DeBoer, distributed by the GeoDa Center and packaged in the R package agridat; the Rdatasets project mirrors it as CSV. Yield is in quintals per hectare, divided by 10 here for t/ha. The trial has no strip polygons, so the setup draws each strip as a 9.8 m wide band along the line through its points.

The design has a flaw worth knowing: the nitrogen order was randomised once, N5, N3, N4, N0, N2, N1, and repeated in all three replicates, so every treatment always has the same neighbours.

import os

import geopandas as gpd
import numpy as np
import pandas as pd
from shapely.geometry import LineString

if not os.path.exists("lasrosas.csv"):
    pd.read_csv("https://vincentarelbundock.github.io/Rdatasets/csv/agridat/lasrosas.corn.csv").to_csv(
        "lasrosas.csv", index=False)
raw = pd.read_csv("lasrosas.csv")
yields = gpd.GeoDataFrame(
    {"year": raw["year"], "rep": raw["rep"], "treatment": raw["nf"], "nitro": raw["nitro"],
     "yield_t_ha": raw["yield"] / 10, "strip_id": raw["year"].astype(str) + "-" + raw["rep"] + "-" + raw["nf"]},
    geometry=gpd.points_from_xy(raw["long"], raw["lat"]), crs=4326).to_crs(32720)

rows = []
for sid, pts in yields.groupby("strip_id"):
    xy = np.c_[pts.geometry.x, pts.geometry.y]
    axis = np.linalg.svd(xy - xy.mean(axis=0))[2][0]              # the strip's long axis
    along = (xy - xy.mean(axis=0)) @ axis
    line = LineString([xy.mean(axis=0) + along.min() * axis, xy.mean(axis=0) + along.max() * axis])
    rows.append({"strip_id": sid, "year": pts["year"].iloc[0], "rep": pts["rep"].iloc[0],
                 "treatment": pts["treatment"].iloc[0], "line": line,
                 "geometry": line.buffer(4.9, cap_style="flat")})
strips = gpd.GeoDataFrame(rows, crs=yields.crs)
print(f"{len(yields):,} yield points in {len(strips)} strips; strip length "
      f"{strips.line.map(lambda l: l.length).median():.0f} m; "
      f"order in 1999 rep R1: {list(strips[(strips.year == 1999) & (strips.rep == 'R1')].sort_values('strip_id', key=lambda s: strips.loc[s.index].geometry.centroid.y).treatment)}")

3,443 yield points in 36 strips; strip length 654 m; order in 1999 rep R1: ['N5', 'N3', 'N4', 'N0', 'N2', 'N1']

Example 1 — extract strip yields with the ends trimmed

import numpy as np, geopandas as gpd, pandas as pd

def strip_yields(yields, strips, end_trim_m=30, side_trim_m=0, min_points=50):
    """Mean yield per strip, dropping points near the ends (lag, rate change) and sides."""
    kept = []
    for strip in strips.itertuples():
        pts = yields[yields.strip_id == strip.strip_id]
        along = pts.geometry.map(strip.line.project)
        ok = (along >= end_trim_m) & (along <= strip.line.length - end_trim_m)
        if side_trim_m:
            ok &= pts.within(strip.geometry.buffer(-side_trim_m))
        kept.append(pts[ok].assign(treatment=strip.treatment))
    j = pd.concat(kept)
    g = j.groupby(["strip_id", "treatment"])["yield_t_ha"].agg(
        ["count", "mean", "std"]).reset_index()
    dropped = g[g["count"] < min_points]
    if len(dropped):
        print(f"dropping {len(dropped)} strips with fewer than {min_points} points")
    g = g[g["count"] >= min_points]
    print(g.groupby("treatment")[["count", "mean"]].agg(
        {"count": "sum", "mean": "mean"}).round(3).to_string())
    return g


strip_stats = strip_yields(yields, strips)
print(f"kept {int(strip_stats['count'].sum()):,} of {len(yields):,} points")

           count   mean
treatment              
N0           513  6.312
N1           517  6.694
N2           511  6.793
N3           515  6.891
N4           512  7.085
N5           515  7.104
kept 3,083 of 3,443 points

Example 2 — neighbour pairing and the paired test

import numpy as np, pandas as pd, geopandas as gpd
from scipy import stats

def neighbour_pairs(strip_stats, strips, treated="high", control="standard"):
    centroids = strips.set_index("strip_id").geometry.centroid
    stats_by_id = strip_stats.set_index("strip_id")

    rows = []
    for sid in strip_stats.query("treatment == @treated")["strip_id"]:
        here = centroids.loc[sid]
        others = strip_stats.query("treatment == @control")["strip_id"]
        d = centroids.loc[others].distance(here).sort_values()
        near = d.index[:2]
        if len(near) < 2:
            continue
        rows.append({"strip_id": sid,
                     "treated": stats_by_id.loc[sid, "mean"],
                     "control": float(stats_by_id.loc[near, "mean"].mean()),
                     "neighbours": list(near)})
    p = pd.DataFrame(rows)
    p["diff"] = p["treated"] - p["control"]
    return p

def paired_result(p, alpha=0.05):
    n = len(p)
    d = p["diff"].values
    t = stats.ttest_rel(p["treated"], p["control"])
    se = d.std(ddof=1) / np.sqrt(n)
    half = stats.t.ppf(1 - alpha / 2, n - 1) * se
    return {"pairs": n, "mean_diff": round(float(d.mean()), 3),
            "sd": round(float(d.std(ddof=1)), 3), "se": round(float(se), 3),
            "ci_low": round(float(d.mean() - half), 3),
            "ci_high": round(float(d.mean() + half), 3),
            "p_value": round(float(t.pvalue), 4),
            # (t for alpha/2 + t for 80% power) standard errors; 2.8 only when there are many pairs
            "detectable_at_80_power": round(float(
                (stats.t.ppf(1 - alpha / 2, n - 1) + stats.t.ppf(0.8, n - 1)) * se), 3)}

for year in (1999, 2001):
    season_ids = strips.loc[strips.year == year, "strip_id"]
    pairs = neighbour_pairs(strip_stats[strip_stats.strip_id.isin(season_ids)],
                            strips[strips.year == year], treated="N5", control="N0")
    print(year, paired_result(pairs))

1999 {'pairs': 3, 'mean_diff': 0.949, 'sd': 0.233, 'se': 0.134, 'ci_low': 0.372, 'ci_high': 1.527, 'p_value': 0.0194, 'detectable_at_80_power': 0.72}
2001 {'pairs': 3, 'mean_diff': 0.664, 'sd': 0.147, 'se': 0.085, 'ci_low': 0.298, 'ci_high': 1.029, 'p_value': 0.016, 'detectable_at_80_power': 0.456}

Both seasons show a response to the top rate — about 0.9 t/ha in 1999 and 0.7 in 2001 — with intervals that stay above zero, from only three pairs a season. The detectable difference is the other half of the story: a response below about 0.7 t/ha in 1999, or 0.5 in 2001, would have been invisible to this design.

detectable_at_80_power is the number worth quoting when the result is not significant: it is the difference this trial could have detected, and it tells the farmer whether the trial was worth running. The familiar rule of thumb, 2.8 standard errors, assumes many pairs. With three, the t quantiles on 2 degrees of freedom make it 5.4 standard errors — the first version of this function used 2.8 and reported half the true detectable difference.

Example 3 — check the residual spatial structure

import numpy as np, pandas as pd
from scipy.spatial import cKDTree

def residual_autocorrelation(p, strips, k=3):
    """Moran's I on the paired differences: is there structure the pairing missed?"""
    c = strips.set_index("strip_id").geometry.centroid
    xy = np.c_[c.loc[p.strip_id].x, c.loc[p.strip_id].y]
    d = p["diff"].values
    z = d - d.mean()

    tree = cKDTree(xy)
    _, idx = tree.query(xy, k=min(k + 1, len(xy)))
    W = np.zeros((len(xy), len(xy)))
    for i, row in enumerate(idx[:, 1:]):
        W[i, row] = 1
    W = (W + W.T) / 2

    n, s0 = len(z), W.sum()
    I = (n / s0) * (z @ W @ z) / (z @ z)
    expected = -1 / (n - 1)
    print(f"Moran's I on the differences: {I:+.3f} (expected {expected:+.3f})")
    if I > 0.3:
        print("! the pairing has not removed the spatial trend — use a spatial model")
    return float(I)


both = pd.concat([neighbour_pairs(strip_stats[strip_stats.strip_id.isin(strips.loc[strips.year == y, "strip_id"])],
                                  strips[strips.year == y], treated="N5", control="N0") for y in (1999, 2001)])
residual_autocorrelation(both.reset_index(drop=True), strips, k=2)

Moran's I on the differences: -0.219 (expected -0.200)

With six pairs, two seasons laid over the same ground, the test has almost no power, so treat it as a warning rather than a decision; with eight to twelve pairs it is still weak. A strongly positive value means the differences themselves are spatially clustered, which the pairing was supposed to prevent.

Explanation

Why the pseudo-replication problem is so severe

A combine logs a point every second or so, and a strip 600 m long produces several thousand of them. They are not independent: neighbouring points share soil, weather and the same sensor calibration, and their correlation at short lags is very high. A t-test on the raw points has thousands of degrees of freedom and will report a p-value below 0.001 for a difference of a few kilograms. Aggregating to strip means restores the correct number of experimental units, which is the number of strips.

Why neighbour pairing beats a pooled comparison

Fields have gradients — soil type, drainage, headland effects, previous cropping. If the treated strips are on average slightly further up the gradient, a pooled comparison attributes the gradient to the treatment. Adjacent strips share the local conditions almost exactly, so their difference is close to the treatment effect. This is the same argument that makes randomised block designs standard in plot trials.

Why the edges have to be trimmed

Three effects converge on a strip boundary: the harvester's sensor lag displaces yield by tens of metres along the pass, the spreader or drill has a transition zone of several metres where the rate is changing, and adjacent passes overlap. The first few metres of every strip therefore contain a mixture of both treatments, and including them dilutes the measured effect towards zero.

Why to report the detectable difference

Most on-farm trials cannot detect the effects people hope to find. A trial with eight pairs and a standard error of 0.15 t/ha can detect about 0.49 t/ha at 80% power — 3.26 standard errors on 7 degrees of freedom, not the large-sample 2.8. The Las Rosas seasons, with three pairs each, could detect 0.5–0.7 t/ha. If the expected response is 0.2 t/ha, the trial was never going to show it, and saying so is far more useful than reporting a non-significant p-value that gets read as "no effect".

Edge cases or notes

• Six to eight replicates minimum. Fewer cannot support a conclusion.
• Strips must cross the gradient, not run along it.
• Strip width should be a multiple of the harvester width.
• Randomise the assignment in every replicate. Las Rosas repeated one random order in all three, so N5 always sits beside N3.
• Headlands are their own environment. Exclude or treat separately.
• A single year is a single year. Responses vary with the season.
• Record the as-applied data, not just the prescription.
• Report the detectable difference whenever the result is not significant.

Internal links

• How to clean yield monitor data in Python — the cleaning this depends on
• Yield monitor data explained: what the numbers really are — why the edges are contaminated
• Management zones explained — why zones need a trial to become rates
• How to write a variable-rate prescription map — applying the result
• Spatial autocorrelation explained — why points are not independent
• How to calculate Moran's I in Python — the residual check
• Sample design explained — replication and randomisation
• Spatial leakage explained — the same problem in modelling

FAQ

How do I analyse an on-farm strip trial?

Clean the yield data, trim the strip edges, aggregate to one mean per strip, pair each treated strip with its neighbours, and run a paired test on the strip means.

Why can I not test on the individual yield points?

They are not independent. A strip produces thousands of correlated points, and treating them as independent gives a p-value that is meaningless.

How many strips do I need?

Six to eight replicates of each treatment as a minimum. Fewer cannot support a conclusion whatever the analysis.

Why pair with neighbours instead of pooling?

Because fields have gradients. Adjacent strips share their soil, so their difference is close to the treatment effect; a pooled comparison can attribute the gradient to the treatment.

How much of the strip edge should I trim?

A few metres — enough to clear the application transition, the pass overlap and the harvester's lag smear.

What should I report if the result is not significant?

The effect size, the confidence interval and the difference the trial could have detected at reasonable power. "Not significant" on its own is read as "no effect", which is a different claim.