Three ways a careful pipeline misleads itself

Companion methods note to the cloud-organization study. Each failure below happened in that study's lab book, was caught, corrected and left on the record; none of them changes the paper's result. They are general, so the mechanisms and checks are written for any domain where several instruments, comparators and frozen protocols meet.

Technical note · cloud-organization project · 3 November 2026 · site edition 2026-11-03 · minimal reproduction included below

1. A referee that agrees with both sides

Setup

Two passive cloud masks, the image mask (I) and the MODIS mask (M), disagree on part of each scene, and part of that disagreement repeats across the two overpasses three hours apart. To ask whether the repeating part is one instrument's error, both were compared with an active lidar (CALIPSO, L) along its track, and the scene's repeating leftover was correlated with each instrument's departure from the lidar, controlling the track cloud amounts of L, I and M (4,947 box-days, rank partial correlations, year bootstrap).

Symptom

The leftover correlated with both departures (+0.065 for the image mask, +0.069 for MODIS; reading BOTH), and in the two regions where the effect was clearest the same held. It was read as a definitional difference shared by both passive instruments.

Mechanism

|I − M| ≤ |I − L| + |M − L|

On any track pixel, wherever I and M disagree, at least one of them departs from L, whichever is right. The leftover itself correlates with the track's I–M disagreement (Spearman +0.155), so the departures must rise with it; with errors of similar size on both sides, both rise. The comparison could not separate which instrument is wrong from how much they disagree.

Correction

With the track's own I–M disagreement added to the controls, the association went to +0.008 and -0.022 (reading NEITHER). The Aqua-pass leftover went from IMAGE SIDE to NEITHER, and every region read NEITHER. The earlier reading was withdrawn in the lab book (entry 36); the repeating leftover's source is open.

Minimal reproduction

Synthetic truth, two instruments with independent errors whose size varies by scene, a near-perfect referee, and a leftover that partly tracks how much the instruments disagree. Nothing is shared between I and M, yet the leftover correlates with both departures; controlling |I−M| removes it (first line of the output below).

2. A comparator that manufactures persistence

Setup

On each pass, the image mask and the MODIS mask disagree on part of the scene. "Persistence" was the correlation, across Terra and Aqua, of each scene's leftover disagreement after a comparator (cloud amounts and image statistics, the structural kernel's 29 coordinates and a glint term) was regressed out on each pass separately. A persistent leftover was read as a property of the scene that both passes share.

Mechanism

If the comparator misses the same nonlinear structure on both passes, both residuals contain the same misfit, and they correlate even when the two passes' errors are independent. Persistence then measures the comparator, not the scene.

Check

Fit a flexible map from the comparator's inputs to each pass's outcome, simulate outcomes as that map plus independent noise of the observed size, and recompute persistence (20 draws). Persistence that misfit alone can produce:

In samples B and C misfit explains little of the observed value; in sample A the observed value is small and misfit could account for most of it. A second check concerns what persistence means: across resolution, threshold and background settings of the image mask, the correlation stayed between 0.31 and 0.46 while the shared covariance varied about 32-fold. A stable correlation does not mean an unchanged shared amount, and "about 40% shared" reads a correlation as a variance share, which needs an equal-error model.

Minimal reproduction

One scene property drives both passes nonlinearly; errors are independent. A linear comparator leaves residuals that correlate strongly across passes; a flexible one leaves none (second line of the output).

3. A dry run that checked the wrong thing

Practice

Every frozen protocol in the study was run once with each outcome replaced by random noise; any named reading other than the null one counted as a pipeline fault to fix before freezing.

Failure

In the last protocol, a control added cubic splines of morning cloud cover with knots fixed on [0, 1]. Cover never falls below about 0.2, so several basis columns are zero almost everywhere. Inside leave-one-year-out ridge with standardisation, a column that varies on a few training rows is divided by a tiny standard deviation, and predictions for the held-out year explode. On noise, every reading still came out null, so the check passed; but the intervals of the arms containing splines were up to 150,696 wide. In the frozen run the same arms gave intervals up to 48,776,111 wide in two samples.

Correction

A declared amendment moved the knots to outcome-free percentiles and dropped near-constant columns, was checked on noise first (widest interval 0.039), and re-ran only the two readings that contained the splines; the other readings stand from the frozen run. All three records (frozen run, amendment, re-run) are in the lab book (entry 39).

Minimal reproduction

Two rare low-cover days in one year are enough: on a pure-noise outcome the out-of-fold R² with fixed knots falls far below zero and single predictions run to tens of standard deviations, while percentile knots stay at zero (third line of the output).

4. Checklist

1. When a third instrument referees two others, control the two instruments' own disagreement on the same footprint, and ask what the comparison gives when errors are independent.
2. Before reading a shared residual as a shared property, simulate outcomes from the fitted comparator plus independent noise and measure how much "persistence" misfit alone produces.
3. Report shared covariance with any persistence correlation; a stable correlation can hide a large change in the shared amount.
4. In a noise dry run, check interval widths and the size of predictions, not only the named readings. On noise, an R² interval much wider than about 0.1 is a fault.
5. Fix basis-function knots from the data's own range (outcome-free), and drop near-constant columns before standardising.
6. When a frozen run fails mechanically, declare the amendment, test it on noise, re-run only what it touches, and keep the failed run on the record.

Reproduce

The script below uses synthetic data only (numpy, scipy, scikit-learn; seed fixed). Output:

• 1. referee: leftover vs |I-L|, |M-L| given track amounts +0.403 +0.382; adding |I-M| +0.055 +0.041
• 2. misfit: residual correlation across passes, linear comparator +0.671; flexible comparator -0.011 (errors independent)
• 3. dry run on noise: fixed knots R2 -0.344, largest |prediction| 30.6 (smallest column sd 5.4e-07); percentile knots R2 -0.0001, largest |prediction| 0.52

"""Minimal reproductions for the methods note (paper/methods_note.html): three ways a careful pipeline misleads itself.
Pure numpy/scipy/sklearn, synthetic data, seed fixed. Run: python methods_note_minimal.py"""
import numpy as np
from scipy.stats import spearmanr, rankdata
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler, SplineTransformer

rng = np.random.default_rng(0)


def partial_rank(x, y, Z):
    """Spearman partial correlation of x and y given columns Z (rank, regress out, correlate residuals)."""
    r = lambda v: rankdata(v)
    Zr = np.column_stack([r(z) for z in Z.T]) if Z.size else np.empty((len(x), 0))
    res = lambda v: r(v) - LinearRegression().fit(Zr, r(v)).predict(Zr) if Zr.shape[1] else r(v) - r(v).mean()
    return float(np.corrcoef(res(x), res(y))[0, 1])


def referee(n=5000):
    """1. A third instrument 'agrees with both sides'. Truth T; two instruments with independent errors whose size
    varies by scene; a near-perfect referee L. A scene-level leftover that tracks how much I and M disagree correlates
    with BOTH departures from L, although nothing is shared between I and M."""
    T = rng.uniform(0.1, 0.9, n); s = rng.uniform(0.01, 0.12, n)          # scene-dependent error size
    I = T + rng.normal(0, 1, n) * s; M = T + rng.normal(0, 1, n) * s; L = T + rng.normal(0, 0.01, n)
    D = np.abs(I - M); leftover = D + rng.normal(0, 0.05, n)             # leftover partly tracks disagreement
    Z = np.column_stack([L, I, M])                                        # controls: track amounts, as in the study
    before = (partial_rank(leftover, np.abs(I - L), Z), partial_rank(leftover, np.abs(M - L), Z))
    after = (partial_rank(leftover, np.abs(I - L), np.column_stack([Z, D])), partial_rank(leftover, np.abs(M - L), np.column_stack([Z, D])))
    return before, after


def misfit(n=4000):
    """2. A misfit comparator manufactures persistence. The same scene property x drives both passes nonlinearly;
    errors are independent. Fit a linear comparator to each pass; its residuals share the misfit."""
    x = rng.uniform(-2, 2, (n, 1)); f = np.sin(2 * x[:, 0]) + 0.5 * x[:, 0] ** 2
    yt = f + rng.normal(0, 0.6, n); ya = f + rng.normal(0, 0.6, n)
    rt = yt - LinearRegression().fit(x, yt).predict(x); ra = ya - LinearRegression().fit(x, ya).predict(x)
    ft = make_pipeline(SplineTransformer(n_knots=8), LinearRegression())
    rt2 = yt - ft.fit(x, yt).predict(x); ra2 = ya - ft.fit(x, ya).predict(x)
    return float(np.corrcoef(rt, ra)[0, 1]), float(np.corrcoef(rt2, ra2)[0, 1])


def dry_run(n=3000, years=8):
    """3. A noise dry run that only checks readings. Morning cover lives above about 0.2; cubic splines with knots fixed
    on [0, 1] leave basis columns that are almost always zero. Two rare low-cover days in one year are enough: in the
    folds that train on other years the column varies on a single row, standardisation divides by a tiny standard
    deviation, and the held-out year's predictions explode. On a pure-noise outcome the gain is still 'not bounded', so
    a readings-only check passes; the out-of-fold R2 and the interval width do not."""
    from sklearn.linear_model import RidgeCV
    r = np.random.default_rng(0)
    cover = np.clip(r.beta(5, 2, n) * 0.79 + 0.21, 0.21, 1.0); yr = r.integers(0, years, n)
    cover[np.nonzero(yr == 3)[0][:2]] = 0.19; cover[np.nonzero(yr == 5)[0][:1]] = 0.205
    y = r.normal(0, 1, n)
    def oof(X):
        p = np.empty(n)
        for t in range(years):
            tr, te = yr != t, yr == t
            p[te] = make_pipeline(StandardScaler(), RidgeCV(alphas=np.logspace(-3, 4, 30))).fit(X[tr], y[tr]).predict(X[te])
        return 1 - np.sum((y - p) ** 2) / np.sum((y - y.mean()) ** 2), float(np.abs(p).max())
    fixed = SplineTransformer(n_knots=6, degree=3, knots=np.linspace(0, 1, 6)[:, None], extrapolation='linear').fit_transform(cover[:, None])
    quant = SplineTransformer(n_knots=6, degree=3, knots=np.quantile(cover, np.linspace(.05, .95, 6))[:, None], extrapolation='linear').fit_transform(cover[:, None])
    quant = quant[:, quant.std(0) >= 1e-3]
    return oof(fixed), oof(quant), float(fixed.std(0).min())


if __name__ == '__main__':
    b, a = referee(); print(f'1. referee: leftover vs |I-L|, |M-L| given track amounts {b[0]:+.3f} {b[1]:+.3f}; adding |I-M| {a[0]:+.3f} {a[1]:+.3f}')
    m1, m2 = misfit(); print(f'2. misfit: residual correlation across passes, linear comparator {m1:+.3f}; flexible comparator {m2:+.3f} (errors independent)')
    (f, fm), (q, qm), sd = dry_run(); print(f'3. dry run on noise: fixed knots R2 {f:+.3f}, largest |prediction| {fm:.1f} (smallest column sd {sd:.1e}); percentile knots R2 {q:+.4f}, largest |prediction| {qm:.2f}')