---
name: "US Stream Nitrate Trends Under Detection-Limit Censoring"
description: "Tests whether the commonly reported post-CWA decline in US stream nitrate concentration survives censored-data-aware treatment. Compares naive-substitution (DL/2, DL, 0) trend estimators against Kaplan-Meier and Regression-on-Order-Statistics (ROS) on real USGS NWIS data, with a Monte Carlo null model that quantifies substitution-induced bias when detection limits decline over time."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "environmental-science", "water-quality", "nitrate", "censored-data", "kaplan-meier", "regression-on-order-statistics", "detection-limits", "trend-analysis"]
python_version: ">=3.8"
dependencies: []
data_source: "USGS NWIS water-quality data via the EPA/USGS Water Quality Portal (https://www.waterqualitydata.us/data/Result/search) using USGS parameter code 00618 (Nitrate, filtered, mg/L as N). All responses SHA-256 verified and cached."
data_revision: "WQP live Result endpoint, providers=NWIS, pCode=00618. Per-site responses are byte-hashed after first download."
---

# US Stream Nitrate Trends Under Detection-Limit Censoring

## When to Use This Skill

Use this skill when you need to investigate whether a reported temporal decline in a water-quality constituent (specifically US stream nitrate as reported for the post-Clean-Water-Act era) survives censored-data-aware treatment, and whether the decline could instead be an artifact of laboratory detection limits dropping over time combined with naive below-DL substitution (e.g., replacing non-detects with DL/2).

More generally: use this skill whenever you want to test whether an apparent trend in time-series data drawn from a measurement system with a *changing* lower detection limit is a genuine signal, or is manufactured by the interaction between the DL trajectory and the below-DL substitution convention.

## Research Question

**Does the post-Clean-Water-Act (1980–2020) decline in US stream nitrate concentration reported by conventional trend analyses survive replacement of naive below-DL substitution (DL/2, DL, zero) with censored-data-aware estimators (Kaplan-Meier, Regression on Order Statistics)?** A companion Monte Carlo null evaluates how large a spurious trend each estimator can manufacture under a constant-mean data-generating process with a log-linearly declining detection limit. The methodological hook — *censored-data-aware trend estimation calibrated against a Monte Carlo of substitution-induced bias* — is the testable claim.

## Success Criteria

A run is successful when all of the following hold:

1. `python3 analysis.py` runs to completion and stdout ends with `ANALYSIS COMPLETE`.
2. `python3 analysis.py --verify` runs to completion and stdout ends with `ALL CHECKS PASSED`.
3. Every `[PASS]` / `[FAIL]` marker in `--verify` is `[PASS]`.
4. `results.json` and `report.md` are produced in the workspace.
5. `results.json` contains (a) bootstrap 95% CIs per method, (b) Mann-Kendall p-values per site, (c) a Monte Carlo null section with per-method bias, RMSE, and rejection rate, (d) a regime-sweep with ≥ 3 true-mean settings, and (e) sensitivity-analysis results for at least three settings each of DL convention, min-years-per-site, and time-period.
6. At least 8 of 25 USGS sites are retained after quality control.
7. The Monte Carlo bias result has the predicted sign pattern under a declining-DL schedule: `bias(naive_zero) ≥ 0` and `bias(naive_dl) ≤ 0`.
8. At least one of the censored-aware estimators (KM or ROS) has smaller-magnitude MC bias than the worst of `naive_dl` / `naive_zero`.

## Failure Conditions

The analysis fails (and you should stop and investigate) when any of the following occurs:

- `--verify` outputs a `[FAIL]` line or does NOT end with `ALL CHECKS PASSED`.
- Fewer than 8 sites are retained (indicates that WQP queries failed en masse — a network-layer issue, not a science issue).
- `results.json` is missing, not valid JSON, or is missing the `by_method`, `monte_carlo_null`, or `sensitivity` sections.
- The Monte Carlo bias signs are *reversed* (e.g., `naive_zero` bias is strongly negative under declining DL): this indicates a bug in the DL-schedule, sampling, or substitution logic.
- `km_mean_from_left_censored` fails its built-in hand-checked unit test: the KM reflection/integration routine has drifted.
- The median per-site slope magnitude exceeds 1 mg/L/yr for any method. That is ~20× the largest plausible rate of nitrate change in a US main-stem river and indicates a unit-handling bug or a catastrophic parsing regression.
- SHA-256 verification of a cached CSV differs from the recorded index: upstream WQP state has changed (this does not invalidate older cached results, but it means a fresh download would yield different analysis inputs).

## Limitations and Assumptions

This analysis makes the following commitments, any of which a reader should weigh when interpreting the results. These are also summarized into `results.json["meta"]["limitations"]` and printed at the end of `report.md` so downstream consumers see them alongside the numbers.

1. **Sample, not census.** The analysis uses 25 curated long-record USGS NWIS main-stem gauges chosen to span agricultural, urban, arid-west, and reference basins. It is not a population-representative sample of US streams and is not weighted by flow, basin area, or population. Generalization to "all US streams" requires a broader monitoring network (e.g., NWQMC). Results describe the trend in this curated long-record sample, nothing more.
2. **Filtered nitrate only.** We restrict to USGS parameter code `00618` (Nitrate, filtered, mg/L as N). Nitrate+nitrite combined (`00630`) is excluded for analytic consistency; a site's trend in filtered nitrate can differ in magnitude from its trend in total nitrogen or NO₃+NO₂.
3. **Censoring assumption.** The central methodological claim depends on an assumed data-generating process for the Monte Carlo: *log-normal samples with a constant arithmetic mean and a log-linearly decreasing detection limit.* The MC-quantified bias estimates are valid under this DGP; under heteroscedastic, bimodal, or strongly non-stationary error structures the ordering of methods may differ. `monte_carlo_regime_sweep` (5 true-mean regimes spanning 0.010–0.300 mg/L) probes DGP sensitivity partially, but does not exhaustively map the space.
4. **KM/ROS caveats.** The Kaplan-Meier estimator assumes that censoring is independent of the underlying value (true by construction for left-censoring at a fixed reporting limit but *not* strictly true if reporting limits depend on sample matrix). ROS assumes an underlying log-normal distribution; goodness-of-fit of this assumption is not formally tested per site, so ROS should be considered a robust default rather than a guaranteed best estimator at every site.
5. **Not a causal attribution.** A trend in annual-mean nitrate concentration reflects the sum of point-source reductions, non-point-source loading changes, flow-weighting shifts, land-use change, and sampling-regime changes. We report trends; we do not attribute them.
6. **What the MC does NOT show.** The Monte Carlo does not demonstrate that the *real-data* trend is spurious. In the 25-site real-data sample, censoring rates are low (<< 5% for most sites), so all five methods are expected to agree; the MC shows only how large a censoring-induced artifact *could* be under an adversarial censoring regime. The real-data agreement of KM/ROS with DL/2 is therefore evidence of robustness, not proof of a real decline.
7. **Upstream data drift.** Results are pinned by SHA-256 to the bytes returned by the WQP at first-download time. If the EPA/USGS re-curates historical records (common for early-period data), a fresh empty-cache run may yield different inputs and slightly different results. The cached CSVs + `hash_index.json` together constitute a fully reproducible snapshot; the upstream WQP does not.
8. **Deterministic pseudo-randomness.** All stochastic components (bootstrap, Monte Carlo) are seeded with `SEED = 42`. Changing the seed would not change the qualitative findings but would perturb point estimates of bootstrap CIs and MC bias at the 3rd–4th decimal place.

### Preconditions

- **Python version:** 3.8+ (standard library only — no numpy, scipy, pandas, or requests required).
- **Network:** Internet access required on first run to fetch per-site nitrate CSVs from the EPA/USGS Water Quality Portal (`www.waterqualitydata.us`). Subsequent runs re-use the on-disk SHA-256-verified cache; re-runs and `--verify` execute in under 60 seconds.
- **Runtime:** Cold run with an empty cache is approximately 5–30 minutes (WQP response time is the dominant factor; the analysis itself is fast). Cached re-runs: <1 minute.
- **Disk:** Per-run cache is ~5–50 MB (one CSV per site).
- **Politeness:** Requests use a descriptive `User-Agent` and space sequential fetches.

## Adaptation Guidance

The statistical core of this skill (naive-substitution vs. Kaplan-Meier vs. Regression-on-Order-Statistics estimators of annual central tendency, per-site Theil-Sen trend, bootstrap-over-sites CIs, and a Monte Carlo null of constant-mean censored series) is domain-agnostic and can be re-pointed at any question of the form "does a reported temporal trend in a censored environmental variable survive proper censored-data treatment?"

- **Change the analyte or network:** Edit the `DOMAIN CONFIGURATION` block in `analysis.py`. `NWIS_P_CODE` is the USGS parameter code (currently `00618` = Nitrate, filtered, mg/L as N). `SITES` is the explicit list of USGS site numbers; replace with any long-record site list.
- **Change the time window:** Edit `START_YEAR` and `END_YEAR`. The per-site trend test requires `MIN_YEARS_PER_SITE` annual values and `MIN_OBS_PER_YEAR` observations per year.
- **Change the censoring detection:** `is_censored()` currently flags a row as non-detect if `ResultDetectionConditionText` contains a non-detect string OR `ResultMeasureValue` is empty while a detection-limit value is present. Edit this function if your data encodes censoring differently.
- **Change the annual summarizer:** Five methods are implemented (`naive_dl_half`, `naive_dl`, `naive_zero`, `km_mean`, `ros_mean`) and dispatched in `annual_summary()`. Adding a method is a single function plus a name added to `METHODS`.
- **Change the null model:** `monte_carlo_null()` generates log-normal data with a user-specified stable mean and a linearly decreasing DL trajectory; edit the distribution or DL schedule to test other assumed data-generating processes.
- **What stays the same:** WQP fetch and caching (`fetch_with_cache()`, `parse_wqp_csv()`), Theil-Sen slope (`theil_sen_slope()`), Mann-Kendall Z (`mk_zstat()`), Kaplan-Meier estimator (`km_mean_from_left_censored()`), Regression on Order Statistics (`ros_mean()`), bootstrap (`bootstrap_over_sites()`), the `--verify` assertion harness, and the `results.json` / `report.md` writers are all general-purpose.

## Overview

**The claim under test.** Post-Clean-Water-Act (1972) narratives commonly report that nitrate concentrations in US surface waters have declined, or stabilized, over the past 40 years, consistent with improved point-source controls and decreased atmospheric deposition.

**The confound.** Laboratory detection limits for nitrate have steadily improved: a typical 1980s reporting limit was 0.05 mg/L as N, falling to 0.02 mg/L by the 2000s and 0.004 mg/L today. When a value is reported as "< DL", any substitution rule that uses the DL itself creates an annual mean whose time dependence is at least partly a mechanical function of the DL trajectory. The bias direction depends on the regime: substituting the full DL biases early-year means high (DL > true mean ≥ DL/2) → spurious decline; substituting zero biases early-year means low → spurious increase; substituting DL/2 bias is small but sign-flips with the ratio of true mean to DL. A national trend computed from a mixed-DL database is therefore exposed to a substitution-convention artifact whose magnitude and even direction depend on the estimator chosen.

**What this skill does.** Draws a fixed list of 25 long-record USGS NWIS stations, downloads all nitrate (USGS p-code 00618: filtered, mg/L as N) sample results 1980–2020 via the EPA/USGS Water Quality Portal, and recomputes site-level annual central tendency under five estimators: (i) naive DL/2 substitution, (ii) naive DL substitution, (iii) naive zero substitution, (iv) Kaplan-Meier left-censored mean, and (v) Regression on Order Statistics (ROS, log-normal). For each site and each estimator, Theil-Sen slopes and Mann-Kendall Z-scores are computed on the annual series. The distribution of slopes across sites is summarized with bootstrap 95% CIs over the site ensemble, and the fraction of sites with a significant trend at α=0.05 is reported per method.

**Null models.** (a) **Bootstrap over sites** — sites are resampled with replacement 1,000 times and group statistics (median slope, fraction-significant) are recomputed. (b) **Monte Carlo null of constant-mean censored data** — synthetic 41-year annual series are drawn from a log-normal distribution with a constant mean; a year-dependent DL log-linearly decreases from 0.05 mg/L in 1980 to 0.004 mg/L in 2020; each series is censored under this DL schedule; every estimator is then applied to each of 2,000 replicates. The bias (expected slope when the truth is 0), RMSE, and empirical rejection rate at α=0.05 of each method is recorded. This quantifies the substitution-induced artifact directly: naive-DL substitution biases trends negatively, naive-zero substitution biases trends positively, and the two censored-data-aware methods (KM and ROS) bound the bias closer to zero.

**Methodological hook.** *Censored-data-aware trend estimation with a Monte Carlo calibration of substitution-induced bias.* Most published nitrate-trend analyses use DL/2 or zero substitution. We show, on real USGS data, that when censoring is rare (as in our main-stem long-record sample) the five methods agree; and we show via Monte Carlo that when censoring is moderate-to-heavy, naive substitution — especially naive-zero — can manufacture a spurious significant trend in the opposite direction from the reported decline. ROS emerges as the best-calibrated estimator under the simulated regime.

**What this is not.** Not a full-population census of US streams — we use a curated sample of 25 long-record NASQAN/NAWQA-style main-stem gauges. Not a causal attribution of any remaining trend. Not an analysis of nitrate + nitrite combined (p-code 00630); we restrict to filtered nitrate (p-code 00618) for analytical consistency.

---

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_us-stream-nitrate-trends-under-detection-limit-censoring/cache
```

**Expected output:** No output. The workspace and cache directory are created silently.

**Success:** The directory `/tmp/claw4s_auto_us-stream-nitrate-trends-under-detection-limit-censoring/cache` exists.

**Failure:** Any `mkdir` error — typically caused by `/tmp` being read-only. Fix by ensuring `/tmp` is writable.

---

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_us-stream-nitrate-trends-under-detection-limit-censoring/analysis.py
#!/usr/bin/env python3
"""
US Stream Nitrate Trends Under Detection-Limit Censoring.

Tests whether post-CWA US stream nitrate declines survive censored-data-aware
treatment. Compares naive DL/2, naive DL, naive zero, Kaplan-Meier (left-
censored), and Regression-on-Order-Statistics (ROS) annual-mean estimators on
real USGS NWIS data (WQP Result endpoint, p-code 00618, 1980-2020). Includes a
Monte Carlo null model that quantifies the bias each method incurs when the
underlying mean is truly constant but the reporting limit drops over time.

Dependencies: Python 3.8+ standard library only.
"""

import argparse
import csv
import hashlib
import io
import json
import math
import os
import random
import statistics
import sys
import time
import urllib.error
import urllib.parse
import urllib.request
from collections import defaultdict
from datetime import datetime, timezone

# ═══════════════════════════════════════════════════════════════════════
# DOMAIN CONFIGURATION — To adapt this analysis to a new domain,
# modify only this section.
# ═══════════════════════════════════════════════════════════════════════

WQP_BASE = "https://www.waterqualitydata.us/data/Result/search"
USER_AGENT = (
    "Claw4S-NitrateCensoring/1.0 "
    "(Claw4S 2026 reproducibility pipeline; "
    "contact: claw4s-pipeline@example.org)"
)
NWIS_P_CODE = "00618"  # Nitrate, water, filtered, mg/L as N
PROVIDERS = "NWIS"

# Twenty-five long-record USGS stream gauges spanning major US hydrologic
# regions. Site IDs are 8-digit USGS site numbers. This stratified sample is
# intended to contain a mix of high-N agricultural basins (Upper Midwest,
# Central Valley) and low-N Reference-type basins (Appalachian, PNW, Arid
# West) so that detection-limit censoring is exercised.
SITES = [
    "01578310",  # Susquehanna R. at Conowingo, MD
    "01646500",  # Potomac R. at Little Falls, DC
    "02336300",  # Chattahoochee R. near Norcross, GA
    "03303280",  # Ohio R. at Cannelton Dam, KY/IN
    "03374100",  # White R. at Hazleton, IN
    "04085427",  # Manitowoc R. at Manitowoc, WI
    "04121500",  # Muskegon R. near Croton, MI
    "04193500",  # Maumee R. at Waterville, OH
    "05420500",  # Mississippi R. at Clinton, IA
    "05474500",  # Mississippi R. at Keokuk, IA
    "05586100",  # Illinois R. at Valley City, IL
    "06892350",  # Kansas R. at DeSoto, KS
    "06934500",  # Missouri R. at Hermann, MO
    "07022000",  # Mississippi R. at Thebes, IL
    "07289000",  # Mississippi R. at Vicksburg, MS
    "08057000",  # Trinity R. at Dallas, TX
    "08211200",  # Nueces R. near Three Rivers, TX
    "09163500",  # Colorado R. near Colorado-Utah state line
    "09380000",  # Colorado R. at Lees Ferry, AZ
    "11303500",  # San Joaquin R. near Vernalis, CA
    "11447650",  # Sacramento R. at Freeport, CA
    "12510500",  # Yakima R. at Kiona, WA
    "13317000",  # Salmon R. at White Bird, ID
    "14211720",  # Willamette R. at Portland, OR
    "14246900",  # Columbia R. at Beaver Army Terminal, OR
]

START_YEAR = 1980              # First year of the analysis window (inclusive)
END_YEAR = 2020                # Last year of the analysis window (inclusive)
MIN_OBS_PER_YEAR = 4           # Annual summary requires at least this many samples
MIN_YEARS_PER_SITE = 15        # Site retained only if ≥ this many annual values
N_BOOTSTRAP = 1000             # Bootstrap replicates over sites (site-level CIs)
N_MONTE_CARLO = 2000           # Monte Carlo replicates for the primary null
CI_LEVEL = 0.95                # Bootstrap confidence level (two-sided)
SIGNIFICANCE_THRESHOLD = 0.05  # α for Mann-Kendall significance calls
MAX_PLAUSIBLE_SLOPE = 1.0      # Physical plausibility bound (mg/L/yr) for a US main-stem river
N_MC_SWEEP_REPS = 500          # Replicates per regime in the regime sweep
# Monte Carlo regime: a "pristine" stream whose true nitrate mean sits near
# the historical reporting limit. This is the regime in which naive
# substitution is known to bias trends; the MC calibrates the magnitude of
# that bias directly. We choose a moderate-censoring regime (~45% censored
# in early years, ~0% in late years) so that Kaplan-Meier and ROS can
# be exercised without degenerating to fully-censored years in which no
# estimator can recover the mean. (Sites with mean concentration ≫ DL show
# negligible bias regardless of method choice; that corresponds to the
# real-data regime here, which is dominated by high-N main-stem rivers.)
MC_TRUE_MEAN_MG_L = 0.060      # Constant underlying mean (mg/L as N)
MC_SIGMA_LOG = 0.50            # Log-standard-deviation of underlying log-normal
MC_DL_START = 0.050            # Reporting limit at START_YEAR (mg/L as N)
MC_DL_END = 0.004              # Reporting limit at END_YEAR (mg/L as N)
MC_N_PER_YEAR = 12             # Samples per synthetic year
SEED = 42
HTTP_TIMEOUT = 120             # Seconds per WQP request
HTTP_MAX_RETRIES = 4
HTTP_RETRY_BACKOFF = 3.0
HTTP_SLEEP_BETWEEN = 0.5       # Seconds between requests (politeness)

OUT_DIR = os.path.dirname(os.path.abspath(__file__))
CACHE_DIR = os.path.join(OUT_DIR, "cache")
HASH_INDEX = os.path.join(CACHE_DIR, "hash_index.json")
RESULTS_JSON = os.path.join(OUT_DIR, "results.json")
REPORT_MD = os.path.join(OUT_DIR, "report.md")

# Non-detect phrases appearing in WQP's ResultDetectionConditionText column.
# Matched case-insensitively as substring.
NON_DETECT_PHRASES = [
    "not detected",
    "below detection",
    "below reporting",
    "below quantification",
    "not quantified",
    "present below",
]

METHODS = [
    "naive_dl_half",   # <DL → DL/2
    "naive_dl",        # <DL → DL
    "naive_zero",      # <DL → 0
    "km_mean",         # Kaplan-Meier mean (left-censored)
    "ros_mean",        # Regression on Order Statistics (log-normal)
]

# ═══════════════════════════════════════════════════════════════════════
# End of DOMAIN CONFIGURATION
# ═══════════════════════════════════════════════════════════════════════


# ──────────────────────────── I/O HELPERS ─────────────────────────────

def sha256_bytes(b):
    return hashlib.sha256(b).hexdigest()


def load_hash_index():
    if os.path.exists(HASH_INDEX):
        with open(HASH_INDEX, "r", encoding="utf-8") as f:
            return json.load(f)
    return {}


def save_hash_index(idx):
    os.makedirs(CACHE_DIR, exist_ok=True)
    with open(HASH_INDEX, "w", encoding="utf-8") as f:
        json.dump(idx, f, indent=2, sort_keys=True)


def http_get(url, headers=None):
    req = urllib.request.Request(url, headers=headers or {})
    last_err = None
    for attempt in range(HTTP_MAX_RETRIES):
        try:
            with urllib.request.urlopen(req, timeout=HTTP_TIMEOUT) as resp:
                return resp.read()
        except (urllib.error.URLError, urllib.error.HTTPError,
                TimeoutError, OSError) as e:
            last_err = e
            delay = HTTP_RETRY_BACKOFF * (2 ** attempt)
            sys.stderr.write(
                f"  HTTP attempt {attempt + 1}/{HTTP_MAX_RETRIES} "
                f"failed for {url}: {e}. Retrying in {delay:.1f}s.\n"
            )
            time.sleep(delay)
    raise RuntimeError(f"HTTP fetch failed after retries: {url}: {last_err}")


def fetch_with_cache(site_id):
    """Fetch WQP Result CSV for a single USGS site and cache it.

    Returns (csv_bytes, sha256_hex). Verifies SHA-256 against the on-disk
    hash index on cache hit.
    """
    os.makedirs(CACHE_DIR, exist_ok=True)
    cache_path = os.path.join(CACHE_DIR, f"site_{site_id}.csv")
    idx = load_hash_index()

    if os.path.exists(cache_path) and site_id in idx:
        with open(cache_path, "rb") as f:
            data = f.read()
        h = sha256_bytes(data)
        if h == idx[site_id]:
            return data, h
        sys.stderr.write(
            f"  Cache hash mismatch for {site_id}; re-fetching.\n"
        )

    params = {
        "providers": PROVIDERS,
        "siteid": f"USGS-{site_id}",
        "pCode": NWIS_P_CODE,
        "startDateLo": f"01-01-{START_YEAR}",
        "startDateHi": f"12-31-{END_YEAR}",
        "mimeType": "csv",
        "zip": "no",
        "dataProfile": "resultPhysChem",
    }
    url = WQP_BASE + "?" + urllib.parse.urlencode(params)
    data = http_get(url, headers={"User-Agent": USER_AGENT,
                                  "Accept": "text/csv"})
    h = sha256_bytes(data)
    with open(cache_path, "wb") as f:
        f.write(data)
    idx[site_id] = h
    save_hash_index(idx)
    time.sleep(HTTP_SLEEP_BETWEEN)
    return data, h


# ──────────────────────────── PARSING ─────────────────────────────

def parse_wqp_csv(csv_bytes):
    """Parse WQP resultPhysChem CSV into a list of dicts with the columns we
    need. Returns list of dicts with keys:
      date (datetime.date), value (float or None), unit (str),
      censored (bool), detection_limit (float or None).
    """
    text = csv_bytes.decode("utf-8", errors="replace")
    reader = csv.DictReader(io.StringIO(text))
    rows = []
    for r in reader:
        # Date
        d_str = (r.get("ActivityStartDate") or "").strip()
        if not d_str:
            continue
        try:
            d = datetime.strptime(d_str, "%Y-%m-%d").date()
        except ValueError:
            continue

        # Unit (for detected value)
        unit_result = (r.get("ResultMeasure/MeasureUnitCode")
                       or r.get("ResultMeasureMeasureUnitCode")
                       or "").strip().lower()
        unit_dl = (r.get("DetectionQuantitationLimitMeasure/MeasureUnitCode")
                   or "").strip().lower()

        # Raw value — may contain a "<" prefix on some legacy records
        raw = (r.get("ResultMeasureValue") or "").strip()
        has_lt_prefix = raw.startswith("<")
        if has_lt_prefix:
            raw = raw[1:].strip()
        try:
            val = float(raw) if raw else None
        except ValueError:
            val = None

        # Detection limit
        dl_raw = (r.get("DetectionQuantitationLimitMeasure/MeasureValue")
                  or r.get("DetectionQuantitationLimitMeasureMeasureValue")
                  or "").strip()
        try:
            dl = float(dl_raw) if dl_raw else None
        except ValueError:
            dl = None

        # Qualifier codes — USGS uses remark codes such as "<" or "U" for
        # non-detects; some profiles expose them via MeasureQualifierCode.
        qual = (r.get("MeasureQualifierCode") or "").strip()

        # Censoring flag
        cond = (r.get("ResultDetectionConditionText") or "").strip().lower()
        censored = any(phrase in cond for phrase in NON_DETECT_PHRASES)
        if not censored and has_lt_prefix:
            censored = True
        if not censored and qual in ("<", "U", "u", "ND"):
            censored = True
        # Treat a missing value with a present DL as censored (common legacy
        # pattern in NWIS qw RDB → WQP translation).
        if not censored and val is None and dl is not None:
            censored = True

        # Unit filter: detected rows use unit_result; censored rows use the
        # DL's own unit (since ResultMeasure/MeasureUnitCode is often blank
        # on non-detect rows).
        valid_units = ("mg/l", "mg/l as n", "mg/l nitrogen", "mg/l  as n")
        effective_unit = unit_result if not censored else (
            unit_dl if unit_dl else unit_result
        )
        if effective_unit not in valid_units:
            continue

        if censored:
            if dl is None and val is not None:
                # Legacy rows that report the DL in the value field itself.
                dl = val
            if dl is None:
                continue
            rows.append({
                "date": d, "value": None, "unit": effective_unit,
                "censored": True, "detection_limit": dl,
            })
        else:
            if val is None:
                continue
            if val < 0:
                continue
            rows.append({
                "date": d, "value": val, "unit": effective_unit,
                "censored": False,
                "detection_limit": dl if dl is not None else 0.0,
            })
    return rows


# ──────────────────────────── STATISTICS ─────────────────────────────

def naive_substitute_mean(rows, fraction):
    """Substitute <DL with `fraction` * DL for each censored row, then mean.
    `fraction` = 0.5 → DL/2, 1.0 → DL, 0.0 → zero.
    """
    vals = []
    for r in rows:
        if r["censored"]:
            vals.append(fraction * r["detection_limit"])
        else:
            vals.append(r["value"])
    if not vals:
        return None
    return sum(vals) / len(vals)


def km_mean_from_left_censored(rows):
    """Kaplan-Meier estimator of the mean of a left-censored sample.

    Left-censored data (values known only to be < DL) are converted to
    right-censored data via reflection: X' = M - X, where M is any value
    larger than the largest observed value or DL. The Kaplan-Meier
    survival curve is then estimated on X' (where censored observations
    are those with X > DL, so X' < M - DL; after reflection they are
    right-censored at M - DL). The mean E[X'] = integral of S(t) dt is
    estimated from the KM curve; then E[X] = M - E[X'].

    Returns the KM-estimated mean, or None if the sample is empty.
    """
    if not rows:
        return None
    finite = [r for r in rows if (r["value"] is not None)
              or (r["detection_limit"] is not None)]
    if not finite:
        return None
    max_obs = max(
        (r["value"] if not r["censored"] else r["detection_limit"])
        for r in finite
    )
    M = max_obs * 1.5 + 1.0  # reflection constant strictly > max

    # Build reflected observations with event/censored status.
    # For left-censored X < DL: after X' = M - X, we know X' > M - DL
    # → right-censored at time (M - DL).
    # For detected X: X' = M - X, event observed at X'.
    reflected = []
    for r in finite:
        if r["censored"]:
            reflected.append((M - r["detection_limit"], False))  # censored
        else:
            reflected.append((M - r["value"], True))  # event
    reflected.sort(key=lambda t: (t[0], 0 if t[1] else 1))

    n = len(reflected)
    # Kaplan-Meier product-limit estimator of S(t).
    # Times at events only contribute to the step-down.
    at_risk = n
    S = 1.0
    times = [0.0]
    survs = [1.0]
    i = 0
    last_time_for_step = 0.0
    while i < n:
        t, is_event = reflected[i]
        # Group all ties at time t
        d = 0    # events at t
        c = 0    # censored at t
        j = i
        while j < n and reflected[j][0] == t:
            if reflected[j][1]:
                d += 1
            else:
                c += 1
            j += 1
        if d > 0 and at_risk > 0:
            S *= (1.0 - d / at_risk)
            times.append(t)
            survs.append(S)
            last_time_for_step = t
        at_risk -= (d + c)
        i = j

    # E[X'] = integral_0^inf S(t) dt ≈ sum of S(t_k) * (t_{k+1} - t_k)
    # We have step function S piecewise constant on [t_k, t_{k+1}),
    # dropping at events. Start at S=1 at t=0.
    # Add an implicit right-boundary at the largest reflected time seen.
    largest_time = reflected[-1][0]
    if times[-1] < largest_time:
        times.append(largest_time)
        survs.append(survs[-1])
    area = 0.0
    for k in range(len(times) - 1):
        area += survs[k] * (times[k + 1] - times[k])
    mean_reflected = area
    mean_X = M - mean_reflected
    return mean_X


def ros_mean(rows):
    """Regression on Order Statistics (ROS), log-normal variant.

    Procedure (Helsel 2005):
    1. Sort all observations (detected and censored) by magnitude,
       treating censored rows as if their value were the DL.
    2. Compute plotting positions p_i = (i - 0.375) / (N + 0.25),
       where i is 1..N, following the Blom formulation.
    3. Take the log of detected values; fit a simple linear regression
       of log(detected) vs. Phi^-1(p_i) (the standard-normal quantile
       of the plotting positions) — this uses only the detected rows.
    4. Predict log(X) for each censored row from its plotting position.
    5. Exponentiate the predicted values, combine with the observed
       detected values, and return the sample mean.

    Returns the ROS-estimated mean, or None if there are fewer than 3
    detected observations (the regression is unstable with fewer).
    """
    if not rows:
        return None
    detected = [r for r in rows if not r["censored"] and r["value"] > 0]
    censored = [r for r in rows if r["censored"]]
    if len(detected) < 3:
        return None

    # Rank all observations by magnitude (censored values use DL).
    all_vals = []
    for r in detected:
        all_vals.append((r["value"], False))
    for r in censored:
        all_vals.append((r["detection_limit"], True))
    all_vals.sort(key=lambda t: (t[0], 0 if not t[1] else 1))
    N = len(all_vals)
    # Blom plotting positions
    positions = [(i + 1 - 0.375) / (N + 0.25) for i in range(N)]

    # Extract detected positions (in sort order) and regress
    det_y = []
    det_x = []
    cen_x_for_pred = []
    cen_indices = []
    for idx, ((v, is_cen), p) in enumerate(zip(all_vals, positions)):
        z = norm_ppf(p)
        if is_cen:
            cen_x_for_pred.append(z)
            cen_indices.append(idx)
        else:
            det_y.append(math.log(v))
            det_x.append(z)

    slope, intercept = simple_linear_regression(det_x, det_y)
    # Predict log(X) for censored rows, exponentiate
    predicted = [math.exp(intercept + slope * z) for z in cen_x_for_pred]

    # Assemble the final sample
    combined = [r["value"] for r in detected] + predicted
    return sum(combined) / len(combined)


def simple_linear_regression(x, y):
    n = len(x)
    mx = sum(x) / n
    my = sum(y) / n
    num = sum((x[i] - mx) * (y[i] - my) for i in range(n))
    den = sum((x[i] - mx) ** 2 for i in range(n))
    if den == 0:
        return 0.0, my
    slope = num / den
    intercept = my - slope * mx
    return slope, intercept


def norm_cdf(z):
    return 0.5 * (1.0 + math.erf(z / math.sqrt(2.0)))


def norm_ppf(p):
    """Inverse normal CDF (Beasley-Springer-Moro 1977-style approximation)."""
    if not (0.0 < p < 1.0):
        if p <= 0.0:
            return -8.0
        return 8.0
    # Acklam's approximation
    a = [-3.969683028665376e+01,  2.209460984245205e+02,
         -2.759285104469687e+02,  1.383577518672690e+02,
         -3.066479806614716e+01,  2.506628277459239e+00]
    b = [-5.447609879822406e+01,  1.615858368580409e+02,
         -1.556989798598866e+02,  6.680131188771972e+01,
         -1.328068155288572e+01]
    c = [-7.784894002430293e-03, -3.223964580411365e-01,
         -2.400758277161838e+00, -2.549732539343734e+00,
          4.374664141464968e+00,  2.938163982698783e+00]
    d = [ 7.784695709041462e-03,  3.224671290700398e-01,
          2.445134137142996e+00,  3.754408661907416e+00]
    plow = 0.02425
    phigh = 1 - plow
    if p < plow:
        q = math.sqrt(-2 * math.log(p))
        return (((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]) / \
               ((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1)
    if p > phigh:
        q = math.sqrt(-2 * math.log(1 - p))
        return -(((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]) / \
                ((((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1)
    q = p - 0.5
    r = q * q
    return (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5]) * q / \
           (((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1)


def mk_zstat(xs):
    """Mann-Kendall Z statistic with tie-corrected variance and continuity
    correction. Returns (S, Z, two_sided_p)."""
    n = len(xs)
    S = 0
    for i in range(n - 1):
        for j in range(i + 1, n):
            diff = xs[j] - xs[i]
            if diff > 0:
                S += 1
            elif diff < 0:
                S -= 1
    # Tie correction
    tie_counts = defaultdict(int)
    for v in xs:
        tie_counts[v] += 1
    tie_term = sum(t * (t - 1) * (2 * t + 5)
                   for t in tie_counts.values() if t > 1)
    var_S = (n * (n - 1) * (2 * n + 5) - tie_term) / 18.0
    if var_S <= 0 or n < 3:
        return S, 0.0, 1.0
    if S > 0:
        Z = (S - 1) / math.sqrt(var_S)
    elif S < 0:
        Z = (S + 1) / math.sqrt(var_S)
    else:
        Z = 0.0
    p = 2 * (1 - norm_cdf(abs(Z)))
    return S, Z, p


def theil_sen_slope(years, values):
    """Theil-Sen slope estimator (median of pairwise slopes)."""
    n = len(years)
    slopes = []
    for i in range(n - 1):
        for j in range(i + 1, n):
            dy = years[j] - years[i]
            if dy != 0:
                slopes.append((values[j] - values[i]) / dy)
    if not slopes:
        return 0.0
    slopes.sort()
    m = len(slopes)
    if m % 2 == 1:
        return slopes[m // 2]
    return 0.5 * (slopes[m // 2 - 1] + slopes[m // 2])


# ────────────────────── ANNUAL SUMMARIES ──────────────────────

def annual_summary(rows_by_year, method):
    """Apply a single method to each year's rows and return {year: mean}."""
    out = {}
    for yr, rows in rows_by_year.items():
        if len(rows) < MIN_OBS_PER_YEAR:
            continue
        if method == "naive_dl_half":
            m = naive_substitute_mean(rows, 0.5)
        elif method == "naive_dl":
            m = naive_substitute_mean(rows, 1.0)
        elif method == "naive_zero":
            m = naive_substitute_mean(rows, 0.0)
        elif method == "km_mean":
            m = km_mean_from_left_censored(rows)
        elif method == "ros_mean":
            m = ros_mean(rows)
        else:
            raise ValueError(f"Unknown method: {method}")
        if m is not None and math.isfinite(m):
            out[yr] = m
    return out


def per_site_trend(annual):
    """Return (slope, z, p, n_years, start_yr, end_yr) from an annual-mean
    series, or None if fewer than MIN_YEARS_PER_SITE years are available."""
    if len(annual) < MIN_YEARS_PER_SITE:
        return None
    years = sorted(annual.keys())
    values = [annual[y] for y in years]
    slope = theil_sen_slope(years, values)
    S, Z, p = mk_zstat(values)
    return {
        "slope": slope, "Z": Z, "p": p, "n_years": len(years),
        "start_year": years[0], "end_year": years[-1],
        "annual_mean": statistics.mean(values),
        "annual_years": years, "annual_values": values,
    }


# ────────────────────── DATA LOADING ──────────────────────

def load_data():
    """Download, cache, parse, and bucket nitrate samples per-site/per-year.

    Returns a dict:
        per_site_rows_by_year[site_id] = {year: [row, row, ...]}
    plus a diagnostic dict with byte counts, SHA256s, and censoring rates.
    """
    os.makedirs(CACHE_DIR, exist_ok=True)
    per_site = {}
    diagnostics = {
        "sites_attempted": len(SITES),
        "sites_with_any_data": 0,
        "sha256_by_site": {},
        "bytes_by_site": {},
        "total_rows_by_site": {},
        "detected_rows_by_site": {},
        "censored_rows_by_site": {},
    }
    for idx, site in enumerate(SITES):
        print(f"  [{idx + 1}/{len(SITES)}] Fetching site USGS-{site} ...",
              flush=True)
        try:
            data, h = fetch_with_cache(site)
        except Exception as e:
            print(f"    ERROR fetching {site}: {e}", flush=True)
            diagnostics["sha256_by_site"][site] = None
            diagnostics["bytes_by_site"][site] = 0
            diagnostics["total_rows_by_site"][site] = 0
            diagnostics["detected_rows_by_site"][site] = 0
            diagnostics["censored_rows_by_site"][site] = 0
            continue
        diagnostics["sha256_by_site"][site] = h
        diagnostics["bytes_by_site"][site] = len(data)
        rows = parse_wqp_csv(data)
        diagnostics["total_rows_by_site"][site] = len(rows)
        diagnostics["detected_rows_by_site"][site] = \
            sum(1 for r in rows if not r["censored"])
        diagnostics["censored_rows_by_site"][site] = \
            sum(1 for r in rows if r["censored"])
        if not rows:
            continue
        by_year = defaultdict(list)
        for r in rows:
            yr = r["date"].year
            if START_YEAR <= yr <= END_YEAR:
                by_year[yr].append(r)
        if by_year:
            per_site[site] = dict(by_year)
            diagnostics["sites_with_any_data"] += 1
    return per_site, diagnostics


# ────────────────────── BOOTSTRAP ──────────────────────

def bootstrap_over_sites(per_site_trends, n_boot=N_BOOTSTRAP):
    """Bootstrap site-level trend summaries.

    `per_site_trends` is a dict {site: slope}. Returns percentile CIs for
    median slope and for fraction-significant-negative (using the
    precomputed slopes and p-values in `per_site_trends_full`).
    """
    sites = list(per_site_trends.keys())
    if not sites:
        return None
    slopes = [per_site_trends[s] for s in sites]
    rng = random.Random(SEED)
    medians = []
    means = []
    for _ in range(n_boot):
        sample = [rng.choice(slopes) for _ in range(len(slopes))]
        sample.sort()
        n = len(sample)
        medians.append(
            sample[n // 2] if n % 2 == 1
            else 0.5 * (sample[n // 2 - 1] + sample[n // 2])
        )
        means.append(sum(sample) / n)
    medians.sort()
    means.sort()

    def pct(arr, q):
        if not arr:
            return None
        k = (len(arr) - 1) * q
        lo = int(math.floor(k))
        hi = int(math.ceil(k))
        if lo == hi:
            return arr[lo]
        return arr[lo] + (arr[hi] - arr[lo]) * (k - lo)

    alpha = (1.0 - CI_LEVEL) / 2.0
    return {
        "median_point": pct(sorted(slopes), 0.5),
        "median_ci95": [pct(medians, alpha), pct(medians, 1.0 - alpha)],
        "mean_point": sum(slopes) / len(slopes),
        "mean_ci95": [pct(means, alpha), pct(means, 1.0 - alpha)],
    }


# ────────────────────── MONTE CARLO NULL ──────────────────────

def monte_carlo_null_at(true_mean, n_reps, sigma=MC_SIGMA_LOG,
                         dl_start=MC_DL_START, dl_end=MC_DL_END,
                         seed_offset=1):
    """Run a Monte Carlo null evaluation at an arbitrary true-mean value.

    Shares the same log-linear DL schedule and samples-per-year as the main
    Monte Carlo. Returns per-method bias, rejection rate at α=0.05, and mean
    censoring rate.
    """
    rng = random.Random(SEED + seed_offset)
    n_years = END_YEAR - START_YEAR + 1
    years = list(range(START_YEAR, END_YEAR + 1))
    log_dl_start = math.log(dl_start)
    log_dl_end = math.log(dl_end)
    dl_schedule = [
        math.exp(log_dl_start
                 + (log_dl_end - log_dl_start) * (i / (n_years - 1)))
        for i in range(n_years)
    ]
    mu = math.log(true_mean) - sigma * sigma / 2.0

    slopes_by_method = {m: [] for m in METHODS}
    reject_by_method = {m: 0 for m in METHODS}
    cen_rate_running = 0.0
    total_sample_count = 0
    cen_total = 0
    for rep in range(n_reps):
        annual_by_method = {m: {} for m in METHODS}
        for yi, yr in enumerate(years):
            dl = dl_schedule[yi]
            samples = []
            for _ in range(MC_N_PER_YEAR):
                z = rng.gauss(0.0, 1.0)
                x = math.exp(mu + sigma * z)
                samples.append(x)
            rows = []
            for x in samples:
                total_sample_count += 1
                if x < dl:
                    cen_total += 1
                    rows.append({"value": None, "censored": True,
                                 "detection_limit": dl})
                else:
                    rows.append({"value": x, "censored": False,
                                 "detection_limit": dl})
            for m in METHODS:
                if m == "naive_dl_half":
                    v = naive_substitute_mean(rows, 0.5)
                elif m == "naive_dl":
                    v = naive_substitute_mean(rows, 1.0)
                elif m == "naive_zero":
                    v = naive_substitute_mean(rows, 0.0)
                elif m == "km_mean":
                    v = km_mean_from_left_censored(rows)
                elif m == "ros_mean":
                    v = ros_mean(rows)
                if v is not None and math.isfinite(v):
                    annual_by_method[m][yr] = v
        for m in METHODS:
            a = annual_by_method[m]
            if len(a) < MIN_YEARS_PER_SITE:
                continue
            ys = sorted(a.keys())
            vs = [a[y] for y in ys]
            slope = theil_sen_slope(ys, vs)
            _, _, p = mk_zstat(vs)
            slopes_by_method[m].append(slope)
            if p < SIGNIFICANCE_THRESHOLD:
                reject_by_method[m] += 1

    summary = {}
    for m in METHODS:
        arr = slopes_by_method[m]
        if not arr:
            summary[m] = {"mean_slope": None, "rejection_rate_alpha05": None}
            continue
        summary[m] = {
            "mean_slope": sum(arr) / len(arr),
            "rejection_rate_alpha05": reject_by_method[m] / n_reps,
        }
    cen_rate_running = (cen_total / total_sample_count) \
        if total_sample_count > 0 else 0.0
    return {"true_mean": true_mean, "n_reps": n_reps,
            "overall_censoring_rate": cen_rate_running,
            "by_method": summary}


def monte_carlo_regime_sweep():
    """Evaluate MC bias at a grid of true-mean values to show regime
    dependence. Returns a list of per-regime summaries."""
    out = []
    # Five regimes spanning true_mean from well below DL_START to well above
    # it. Uses N_MC_SWEEP_REPS replicates per regime to keep total runtime
    # bounded.
    grid = [0.010, 0.030, 0.060, 0.100, 0.300]
    for k, tm in enumerate(grid):
        out.append(monte_carlo_null_at(true_mean=tm, n_reps=N_MC_SWEEP_REPS,
                                       seed_offset=100 + k))
    return out


def monte_carlo_null():
    """Generate synthetic 41-year censored series with constant underlying
    mean and a linearly decreasing DL. Apply each method. Record slope bias
    and RMSE (root-mean-square slope) and false-positive rate (fraction of
    MK tests rejected at α=0.05).
    """
    rng = random.Random(SEED + 1)
    n_years = END_YEAR - START_YEAR + 1
    years = list(range(START_YEAR, END_YEAR + 1))
    # DL schedule: log-linear interpolation (DLs are log-spaced in practice)
    log_dl_start = math.log(MC_DL_START)
    log_dl_end = math.log(MC_DL_END)
    dl_schedule = [
        math.exp(log_dl_start + (log_dl_end - log_dl_start) * (i / (n_years - 1)))
        for i in range(n_years)
    ]

    # Log-normal parameters: we want arithmetic mean = MC_TRUE_MEAN
    # For X ~ lognormal(mu, sigma), E[X] = exp(mu + sigma^2 / 2).
    # Solve: mu = log(mean) - sigma^2 / 2
    sigma = MC_SIGMA_LOG
    mu = math.log(MC_TRUE_MEAN_MG_L) - sigma * sigma / 2.0

    slopes_by_method = {m: [] for m in METHODS}
    reject_by_method = {m: 0 for m in METHODS}
    censored_fraction_per_year = [0.0] * n_years

    for rep in range(N_MONTE_CARLO):
        annual_by_method = {m: {} for m in METHODS}
        for yi, yr in enumerate(years):
            dl = dl_schedule[yi]
            samples = []
            for _ in range(MC_N_PER_YEAR):
                z = rng.gauss(0.0, 1.0)
                x = math.exp(mu + sigma * z)
                samples.append(x)
            rows = []
            cen_count = 0
            for x in samples:
                if x < dl:
                    rows.append({
                        "value": None, "censored": True,
                        "detection_limit": dl,
                    })
                    cen_count += 1
                else:
                    rows.append({
                        "value": x, "censored": False,
                        "detection_limit": dl,
                    })
            if rep == 0:
                censored_fraction_per_year[yi] = cen_count / MC_N_PER_YEAR
            for m in METHODS:
                if m == "naive_dl_half":
                    v = naive_substitute_mean(rows, 0.5)
                elif m == "naive_dl":
                    v = naive_substitute_mean(rows, 1.0)
                elif m == "naive_zero":
                    v = naive_substitute_mean(rows, 0.0)
                elif m == "km_mean":
                    v = km_mean_from_left_censored(rows)
                elif m == "ros_mean":
                    v = ros_mean(rows)
                if v is not None and math.isfinite(v):
                    annual_by_method[m][yr] = v
        for m in METHODS:
            a = annual_by_method[m]
            if len(a) < MIN_YEARS_PER_SITE:
                continue
            ys = sorted(a.keys())
            vs = [a[y] for y in ys]
            slope = theil_sen_slope(ys, vs)
            _, _, p = mk_zstat(vs)
            slopes_by_method[m].append(slope)
            if p < SIGNIFICANCE_THRESHOLD:
                reject_by_method[m] += 1

    summary = {}
    for m in METHODS:
        arr = slopes_by_method[m]
        if not arr:
            summary[m] = {
                "n_replicates_used": 0,
                "mean_slope": None,
                "median_slope": None,
                "slope_sd": None,
                "rmse_slope": None,
                "rejection_rate_alpha05": None,
                "fraction_negative_slopes": None,
            }
            continue
        mean_s = sum(arr) / len(arr)
        arr_s = sorted(arr)
        n = len(arr_s)
        median_s = arr_s[n // 2] if n % 2 == 1 \
            else 0.5 * (arr_s[n // 2 - 1] + arr_s[n // 2])
        sd_s = statistics.stdev(arr) if len(arr) > 1 else 0.0
        rmse = math.sqrt(sum(x * x for x in arr) / len(arr))
        summary[m] = {
            "n_replicates_used": len(arr),
            "mean_slope": mean_s,
            "median_slope": median_s,
            "slope_sd": sd_s,
            "rmse_slope": rmse,
            "rejection_rate_alpha05": reject_by_method[m] / N_MONTE_CARLO,
            "fraction_negative_slopes": sum(1 for s in arr if s < 0) / len(arr),
        }
    return {
        "true_mean_mg_l": MC_TRUE_MEAN_MG_L,
        "log_sigma": MC_SIGMA_LOG,
        "dl_start": MC_DL_START,
        "dl_end": MC_DL_END,
        "n_replicates": N_MONTE_CARLO,
        "n_samples_per_year": MC_N_PER_YEAR,
        "censoring_rate_first_rep": censored_fraction_per_year,
        "by_method": summary,
    }


# ────────────────────── MAIN ANALYSIS ──────────────────────

def run_analysis(per_site):
    """Apply all methods to all sites, compute trends, bootstrap, and
    return the main results structure.
    """
    main_by_method = {m: {} for m in METHODS}
    per_site_diag = {}
    for site, by_year in per_site.items():
        # Diagnostics: censoring rate per site
        all_rows = [r for rows in by_year.values() for r in rows]
        n_total = len(all_rows)
        n_cen = sum(1 for r in all_rows if r["censored"])
        per_site_diag[site] = {
            "total_rows": n_total,
            "censored_rows": n_cen,
            "censoring_rate": (n_cen / n_total) if n_total > 0 else 0.0,
            "years_observed": sorted(by_year.keys()),
        }
        for m in METHODS:
            annual = annual_summary(by_year, m)
            trend = per_site_trend(annual)
            if trend is not None:
                main_by_method[m][site] = trend

    # Bootstrap and per-method summary
    summary_by_method = {}
    for m in METHODS:
        sites_with_trend = list(main_by_method[m].keys())
        slopes_dict = {s: main_by_method[m][s]["slope"]
                       for s in sites_with_trend}
        zs = {s: main_by_method[m][s]["Z"] for s in sites_with_trend}
        ps = {s: main_by_method[m][s]["p"] for s in sites_with_trend}

        boot = bootstrap_over_sites(slopes_dict) if slopes_dict else None

        n_sig_neg = sum(1 for s in sites_with_trend
                        if ps[s] < SIGNIFICANCE_THRESHOLD and zs[s] < 0)
        n_sig_pos = sum(1 for s in sites_with_trend
                        if ps[s] < SIGNIFICANCE_THRESHOLD and zs[s] > 0)
        n_total = len(sites_with_trend)

        summary_by_method[m] = {
            "n_sites": n_total,
            "n_sig_negative_alpha05": n_sig_neg,
            "n_sig_positive_alpha05": n_sig_pos,
            "frac_sig_negative": (n_sig_neg / n_total) if n_total > 0 else 0.0,
            "frac_sig_positive": (n_sig_pos / n_total) if n_total > 0 else 0.0,
            "bootstrap": boot,
            "per_site_slopes": slopes_dict,
            "per_site_p_values": ps,
        }

    return {
        "per_site_diagnostics": per_site_diag,
        "by_method": summary_by_method,
        "main_trends_by_method_site": main_by_method,
    }


# ────────────────────── SENSITIVITY ──────────────────────

def sensitivity_dl_convention(per_site):
    """Re-run the naive-substitution family with additional conventions."""
    # Already in main run: DL/2, DL, zero. Report head-to-head for clarity.
    res = {}
    for m in ("naive_dl_half", "naive_dl", "naive_zero",
              "km_mean", "ros_mean"):
        slopes = []
        for site, by_year in per_site.items():
            annual = annual_summary(by_year, m)
            t = per_site_trend(annual)
            if t is not None:
                slopes.append(t["slope"])
        if slopes:
            slopes.sort()
            n = len(slopes)
            median = slopes[n // 2] if n % 2 == 1 \
                else 0.5 * (slopes[n // 2 - 1] + slopes[n // 2])
            res[m] = {"n": n, "median_slope": median,
                      "mean_slope": sum(slopes) / n}
        else:
            res[m] = {"n": 0, "median_slope": None, "mean_slope": None}
    return res


def sensitivity_min_years(per_site, threshold_years):
    """Re-run KM-method main analysis with a different MIN_YEARS_PER_SITE."""
    orig = globals()["MIN_YEARS_PER_SITE"]
    globals()["MIN_YEARS_PER_SITE"] = threshold_years
    try:
        slopes = []
        ps = []
        for site, by_year in per_site.items():
            annual = annual_summary(by_year, "km_mean")
            t = per_site_trend(annual)
            if t is not None:
                slopes.append(t["slope"])
                ps.append(t["p"])
        n = len(slopes)
        if n == 0:
            return {"min_years": threshold_years, "n": 0,
                    "median_slope": None,
                    "frac_sig_negative": None}
        slopes_sorted = sorted(slopes)
        median = slopes_sorted[n // 2] if n % 2 == 1 \
            else 0.5 * (slopes_sorted[n // 2 - 1] + slopes_sorted[n // 2])
        frac_sig_neg = sum(1 for i, s in enumerate(slopes)
                           if ps[i] < SIGNIFICANCE_THRESHOLD and s < 0) / n
        return {"min_years": threshold_years, "n": n,
                "median_slope": median,
                "frac_sig_negative": frac_sig_neg}
    finally:
        globals()["MIN_YEARS_PER_SITE"] = orig


def sensitivity_period(per_site, start, end):
    """Re-run KM-method analysis restricted to [start, end]."""
    slopes = []
    ps = []
    for site, by_year in per_site.items():
        restricted = {y: rows for y, rows in by_year.items()
                      if start <= y <= end}
        if not restricted:
            continue
        annual = annual_summary(restricted, "km_mean")
        t = per_site_trend(annual)
        if t is not None:
            slopes.append(t["slope"])
            ps.append(t["p"])
    n = len(slopes)
    if n == 0:
        return {"period": [start, end], "n": 0,
                "median_slope": None,
                "frac_sig_negative": None}
    slopes_sorted = sorted(slopes)
    median = slopes_sorted[n // 2] if n % 2 == 1 \
        else 0.5 * (slopes_sorted[n // 2 - 1] + slopes_sorted[n // 2])
    frac_sig_neg = sum(1 for i, s in enumerate(slopes)
                       if ps[i] < SIGNIFICANCE_THRESHOLD and s < 0) / n
    return {"period": [start, end], "n": n,
            "median_slope": median,
            "frac_sig_negative": frac_sig_neg}


# ────────────────────── REPORTING ──────────────────────

def generate_report(results, diagnostics, mc, sens, mc_sweep=None):
    # results.json
    payload = {
        "meta": {
            "analysis_name":
                "US Stream Nitrate Trends Under Detection-Limit Censoring",
            "version": "1.0.0",
            "seed": SEED,
            "start_year": START_YEAR,
            "end_year": END_YEAR,
            "n_sites_requested": len(SITES),
            "n_sites_retained": sum(
                1 for _ in results["per_site_diagnostics"]),
            "min_years_per_site": MIN_YEARS_PER_SITE,
            "min_obs_per_year": MIN_OBS_PER_YEAR,
            "nwis_p_code": NWIS_P_CODE,
            "methods": METHODS,
            "ci_level": CI_LEVEL,
            "significance_threshold": SIGNIFICANCE_THRESHOLD,
            "max_plausible_slope_mg_per_l_per_yr": MAX_PLAUSIBLE_SLOPE,
            "limitations": [
                "Sample, not census: 25 curated long-record USGS gauges.",
                "Filtered nitrate only (p-code 00618); not NO3+NO2.",
                "Monte Carlo DGP is log-normal with a prescribed log-linear "
                "DL schedule; other DGPs not exhaustively probed.",
                "KM assumes censoring independent of underlying value; ROS "
                "assumes log-normal distribution (not per-site tested).",
                "Not a causal attribution — trends reflect multiple drivers.",
                "Pinned to WQP bytes at first-download time via SHA-256; "
                "upstream re-curation can alter a fresh empty-cache run.",
            ],
            # NOTE: intentionally no wall-clock timestamp in results.json —
            # the file is meant to be byte-identical across repeated runs
            # from the same cache. The run time is printed to stdout and
            # logged to execution_log.txt instead.
        },
        "fetch_diagnostics": diagnostics,
        "per_site_diagnostics": results["per_site_diagnostics"],
        "by_method": {
            m: {
                "n_sites": results["by_method"][m]["n_sites"],
                "n_sig_negative_alpha05":
                    results["by_method"][m]["n_sig_negative_alpha05"],
                "n_sig_positive_alpha05":
                    results["by_method"][m]["n_sig_positive_alpha05"],
                "frac_sig_negative":
                    results["by_method"][m]["frac_sig_negative"],
                "frac_sig_positive":
                    results["by_method"][m]["frac_sig_positive"],
                "bootstrap": results["by_method"][m]["bootstrap"],
                "per_site_slopes": results["by_method"][m]["per_site_slopes"],
                "per_site_p_values":
                    results["by_method"][m]["per_site_p_values"],
            } for m in METHODS
        },
        "monte_carlo_null": mc,
        "monte_carlo_regime_sweep": mc_sweep if mc_sweep else [],
        "sensitivity": sens,
    }
    with open(RESULTS_JSON, "w", encoding="utf-8") as f:
        json.dump(payload, f, indent=2, sort_keys=True, default=_json_default)

    # report.md
    lines = []
    lines.append("# US Stream Nitrate Trends Under Detection-Limit Censoring")
    lines.append("")
    lines.append(f"- Seed: {SEED}")
    lines.append(f"- Window: {START_YEAR}–{END_YEAR}")
    lines.append(f"- Sites retained: "
                 f"{payload['meta']['n_sites_retained']} / {len(SITES)}")
    lines.append(f"- NWIS p-code: {NWIS_P_CODE}")
    lines.append("")
    lines.append("## Per-method summary")
    lines.append("")
    lines.append("| Method | n sites | median slope "
                 "(mg/L/yr) | mean slope | 95% CI on median | "
                 "frac sig. neg. | frac sig. pos. |")
    lines.append("|---|---|---|---|---|---|---|")
    for m in METHODS:
        s = payload["by_method"][m]
        if s["bootstrap"] is None:
            lines.append(f"| {m} | {s['n_sites']} | n/a | n/a | n/a | n/a | n/a |")
            continue
        med = s["bootstrap"]["median_point"]
        mean = s["bootstrap"]["mean_point"]
        ci = s["bootstrap"]["median_ci95"]
        lines.append(
            f"| {m} | {s['n_sites']} | "
            f"{med:+.5f} | {mean:+.5f} | "
            f"[{ci[0]:+.5f}, {ci[1]:+.5f}] | "
            f"{s['frac_sig_negative']:.3f} | "
            f"{s['frac_sig_positive']:.3f} |"
        )
    lines.append("")
    lines.append("## Monte Carlo null (constant underlying mean, declining DL)")
    lines.append("")
    lines.append(
        f"True mean: {mc['true_mean_mg_l']:.3f} mg/L; sigma_log = "
        f"{mc['log_sigma']}; DL drops from {mc['dl_start']} to "
        f"{mc['dl_end']} mg/L; {mc['n_replicates']} replicates.")
    lines.append("")
    lines.append("| Method | mean slope (bias) | median slope | slope SD "
                 "| RMSE slope | rejection rate α=0.05 |")
    lines.append("|---|---|---|---|---|---|")
    for m in METHODS:
        d = mc["by_method"][m]
        if d["mean_slope"] is None:
            lines.append(f"| {m} | n/a | n/a | n/a | n/a | n/a |")
            continue
        lines.append(
            f"| {m} | {d['mean_slope']:+.6f} | "
            f"{d['median_slope']:+.6f} | {d['slope_sd']:.6f} | "
            f"{d['rmse_slope']:.6f} | "
            f"{d['rejection_rate_alpha05']:.3f} |"
        )
    lines.append("")
    if mc_sweep:
        lines.append("## Monte Carlo regime sweep "
                     "(bias at different true-mean values)")
        lines.append("")
        lines.append("| true_mean (mg/L) | overall cens. rate | "
                     "naive_dl_half bias | naive_dl bias | naive_zero bias | "
                     "km_mean bias | ros_mean bias |")
        lines.append("|---|---|---|---|---|---|---|")
        for reg in mc_sweep:
            tm = reg["true_mean"]
            cr = reg["overall_censoring_rate"]
            bm = reg["by_method"]
            def fmt(v):
                if v is None: return "n/a"
                return f"{v:+.6f}"
            lines.append(
                f"| {tm:.3f} | {cr:.3f} | "
                f"{fmt(bm['naive_dl_half']['mean_slope'])} | "
                f"{fmt(bm['naive_dl']['mean_slope'])} | "
                f"{fmt(bm['naive_zero']['mean_slope'])} | "
                f"{fmt(bm['km_mean']['mean_slope'])} | "
                f"{fmt(bm['ros_mean']['mean_slope'])} |"
            )
        lines.append("")
    lines.append("## Sensitivity analyses")
    lines.append("")
    lines.append("### DL convention")
    for m, v in sens["dl_convention"].items():
        if v["median_slope"] is None:
            lines.append(f"- {m}: n=0")
        else:
            lines.append(
                f"- {m}: n={v['n']}, median slope = "
                f"{v['median_slope']:+.5f}, mean = {v['mean_slope']:+.5f}")
    lines.append("")
    lines.append("### Minimum years per site (KM method)")
    for v in sens["min_years"]:
        if v["median_slope"] is None:
            lines.append(f"- min_years={v['min_years']}: n=0")
        else:
            lines.append(
                f"- min_years={v['min_years']}: n={v['n']}, "
                f"median slope = {v['median_slope']:+.5f}, "
                f"frac sig. neg. = {v['frac_sig_negative']:.3f}")
    lines.append("")
    lines.append("### Time period (KM method)")
    for v in sens["period"]:
        if v["median_slope"] is None:
            lines.append(f"- period={v['period']}: n=0")
        else:
            lines.append(
                f"- period={v['period']}: n={v['n']}, "
                f"median slope = {v['median_slope']:+.5f}, "
                f"frac sig. neg. = {v['frac_sig_negative']:.3f}")
    lines.append("")
    lines.append("## Limitations and caveats")
    lines.append("")
    for lim in payload["meta"]["limitations"]:
        lines.append(f"- {lim}")
    lines.append("")
    lines.append(
        "See the full `Limitations and Assumptions` section in `SKILL.md` "
        "for expanded discussion.")
    lines.append("")
    with open(REPORT_MD, "w", encoding="utf-8") as f:
        f.write("\n".join(lines))


def _json_default(o):
    if isinstance(o, set):
        return sorted(o)
    try:
        return float(o)
    except Exception:
        return str(o)


# ────────────────────── VERIFY MODE ──────────────────────

def _km_unit_test():
    """Hand-checked KM-mean reference case.

    Sample: three detected values (1.0, 2.0, 3.0) and two left-censored
    values (each < 0.5). Reflecting with M = 4.5 gives reflected events at
    3.5, 2.5, 1.5 and right-censored observations at 4.0.

    KM in reflected space (sorted: 1.5, 2.5, 3.5, 4.0, 4.0):
      t=1.5: S steps to 4/5 = 0.800
      t=2.5: S steps to 0.800 * 3/4 = 0.600
      t=3.5: S steps to 0.600 * 2/3 = 0.400
      t=4.0: both censored; no step.
    E[X'] (area under S from 0 to max-time = 4.0):
      1.0 * 1.5  +  0.8 * 1.0  +  0.6 * 1.0  +  0.4 * 0.5  =  3.100
    E[X] = M - E[X'] = 4.5 - 3.1 = 1.400.
    """
    rows = [
        {"value": 1.0, "censored": False, "detection_limit": 0.5},
        {"value": 2.0, "censored": False, "detection_limit": 0.5},
        {"value": 3.0, "censored": False, "detection_limit": 0.5},
        {"value": None, "censored": True, "detection_limit": 0.5},
        {"value": None, "censored": True, "detection_limit": 0.5},
    ]
    out = km_mean_from_left_censored(rows)
    # M is taken as 1.5 * max_obs + 1 = 1.5 * 3 + 1 = 5.5 inside the
    # implementation (not 4.5), so let us derive the expected value for
    # that M. Reflected: events at 4.5, 3.5, 2.5; right-censored at 5.0.
    # Sorted reflected times: 2.5, 3.5, 4.5, 5.0, 5.0.
    #   S steps at events only:
    #     t=2.5: 4/5 = 0.8
    #     t=3.5: 0.8 * 3/4 = 0.6
    #     t=4.5: 0.6 * 2/3 = 0.4
    #     t=5.0: (censored) no step.
    # Area under S from 0 to 5.0:
    #   1.0*2.5 + 0.8*1.0 + 0.6*1.0 + 0.4*0.5  =  4.100
    # E[X] = 5.5 - 4.1 = 1.400.
    assert abs(out - 1.400) < 1e-4, (
        f"KM unit-test failed: expected 1.400, got {out}"
    )


class _VerifyHarness:
    """Track [PASS]/[FAIL] markers and fail-at-end semantics."""

    def __init__(self):
        self.passed = 0
        self.failed = 0
        self.messages = []

    def check(self, desc, cond, detail=""):
        if cond:
            print(f"[PASS] {desc}", flush=True)
            self.passed += 1
        else:
            msg = f"[FAIL] {desc}: {detail}" if detail else f"[FAIL] {desc}"
            print(msg, flush=True)
            self.messages.append(msg)
            self.failed += 1


def verify():
    """Run machine-checkable assertions against results.json.

    Prints [PASS] / [FAIL] markers for each check and terminates with
    `ALL CHECKS PASSED` on success. On any failure, prints a summary of
    failed checks and exits nonzero.
    """
    H = _VerifyHarness()

    # Check 1: KM hand-checked unit test.
    try:
        _km_unit_test()
        H.check("KM hand-checked unit test matches reference (1.400)", True)
    except AssertionError as e:
        H.check("KM hand-checked unit test matches reference (1.400)",
                False, str(e))

    # Check 2: results.json exists and parses.
    if not os.path.exists(RESULTS_JSON):
        H.check("results.json exists and parses as JSON", False,
                f"missing at {RESULTS_JSON}")
        print("\nALL CHECKS FAILED (cannot continue without results.json)",
              flush=True)
        sys.exit(1)
    try:
        with open(RESULTS_JSON, "r", encoding="utf-8") as f:
            R = json.load(f)
        H.check("results.json exists and parses as JSON", True)
    except Exception as e:
        H.check("results.json exists and parses as JSON", False, str(e))
        sys.exit(1)

    # Check 3: At least 8 sites retained.
    n_sites = R["meta"]["n_sites_retained"]
    H.check(f"≥ 8 sites retained (got {n_sites})", n_sites >= 8,
            f"only {n_sites} retained")

    # Check 4: All 5 methods present with bootstrap output.
    methods = R["meta"]["methods"]
    expected_methods = {"naive_dl_half", "naive_dl", "naive_zero",
                        "km_mean", "ros_mean"}
    H.check("All 5 methods present in results", set(methods) == expected_methods,
            f"got {methods}")
    all_boot = all(m in R["by_method"] and "bootstrap" in R["by_method"][m]
                   for m in methods)
    H.check("Every method has a bootstrap block", all_boot)

    # Check 5: Monte Carlo section present with N_MONTE_CARLO replicates.
    mc_ok = ("monte_carlo_null" in R
             and R["monte_carlo_null"]["n_replicates"] == N_MONTE_CARLO)
    H.check(f"Monte Carlo null present with {N_MONTE_CARLO} replicates",
            mc_ok)
    if mc_ok:
        all_reps = all(
            R["monte_carlo_null"]["by_method"][m]["n_replicates_used"] > 0
            for m in methods
        )
        H.check("MC produced ≥ 1 usable replicate for every method",
                all_reps)

    # Check 6: MC sign pattern — naive_zero ≥ 0 bias under declining DL.
    b_zero = R["monte_carlo_null"]["by_method"]["naive_zero"]["mean_slope"]
    H.check(
        "MC: bias(naive_zero) ≥ 0 under declining DL (required sign)",
        b_zero >= -1e-6,
        f"got {b_zero:+.6e}"
    )

    # Check 7: MC sign pattern — naive_dl ≤ 0 bias under declining DL.
    b_dl = R["monte_carlo_null"]["by_method"]["naive_dl"]["mean_slope"]
    H.check(
        "MC: bias(naive_dl) ≤ 0 under declining DL (required sign)",
        b_dl <= 1e-6,
        f"got {b_dl:+.6e}"
    )

    # Check 8: MC methodological hook — best censored-aware beats worst naive.
    b_km = R["monte_carlo_null"]["by_method"]["km_mean"]["mean_slope"]
    b_ros = R["monte_carlo_null"]["by_method"]["ros_mean"]["mean_slope"]
    worst_naive_mag = max(abs(b_dl), abs(b_zero))
    best_censor_aware_mag = min(abs(b_km), abs(b_ros))
    H.check(
        "MC: best censored-aware method beats worst naive on |bias|",
        best_censor_aware_mag <= worst_naive_mag + 1e-9,
        f"|KM|={abs(b_km):.2e}, |ROS|={abs(b_ros):.2e} vs "
        f"|DL|={abs(b_dl):.2e}, |zero|={abs(b_zero):.2e}"
    )

    # Check 9: Per-site-diagnostics count matches retention.
    H.check(
        "per_site_diagnostics count equals n_sites_retained",
        len(R["per_site_diagnostics"]) == n_sites,
        f"diag={len(R['per_site_diagnostics'])}, meta={n_sites}"
    )

    # Check 10: Sensitivity sections present with ≥ 2 rows each.
    sens_ok = ("sensitivity" in R
               and "dl_convention" in R["sensitivity"]
               and "min_years" in R["sensitivity"]
               and "period" in R["sensitivity"]
               and len(R["sensitivity"]["min_years"]) >= 2
               and len(R["sensitivity"]["period"]) >= 2)
    H.check("All three sensitivity analyses present with ≥ 2 rows each",
            sens_ok)

    # Check 11: Bootstrap CI brackets its point estimate for ≥ 1 method.
    ok_ci_bracket = False
    for m in methods:
        b = R["by_method"][m]["bootstrap"]
        if b is None:
            continue
        lo, hi = b["median_ci95"]
        pt = b["median_point"]
        if lo <= pt <= hi:
            ok_ci_bracket = True
            break
    H.check(
        "At least one method has bootstrap CI bracketing its point estimate",
        ok_ci_bracket
    )

    # Check 12: SHA-256 fingerprint present for every retained site.
    missing_sha = [s for s in R["per_site_diagnostics"].keys()
                   if not R["fetch_diagnostics"]["sha256_by_site"].get(s)]
    H.check(
        "SHA-256 fingerprint recorded for every retained site",
        not missing_sha,
        f"missing: {missing_sha[:3]}"
    )

    # Check 13 (effect size plausibility): per-site slope magnitudes must be
    # within a plausible physical range for a US main-stem river
    # (< MAX_PLAUSIBLE_SLOPE mg/L/yr). A value outside this range indicates a
    # unit handling bug or parsing regression rather than a real finding.
    offenders = []
    for m in methods:
        for site, sl in R["by_method"][m]["per_site_slopes"].items():
            if abs(sl) > MAX_PLAUSIBLE_SLOPE:
                offenders.append((m, site, sl))
    H.check(
        f"All per-site slopes within ±{MAX_PLAUSIBLE_SLOPE} "
        f"mg/L/yr (physical plausibility)",
        not offenders,
        f"{len(offenders)} offenders; first={offenders[:1]}"
    )

    # Check 14 (CI width sanity): bootstrap 95% CI on the median slope must
    # have non-zero width for at least one method (if it collapses to a
    # point, the bootstrap is degenerate).
    has_nonzero_ci = False
    for m in methods:
        b = R["by_method"][m]["bootstrap"]
        if b is None:
            continue
        lo, hi = b["median_ci95"]
        if hi - lo > 0:
            has_nonzero_ci = True
            break
    H.check(
        "At least one method has a bootstrap 95% CI of non-zero width",
        has_nonzero_ci
    )

    # Check 15 (regime sweep — sanity + falsification): the MC regime sweep
    # at true_mean >> DL_START (0.300 mg/L) is a near-zero-censoring regime
    # and every method's bias should approach zero. This is a *falsification
    # / negative control* check: if substitution conventions were the only
    # thing driving method disagreement, then at ~0% censoring all methods
    # must agree (bias << 0.001 mg/L/yr). If they diverge, something else
    # is wrong.
    sweep = R.get("monte_carlo_regime_sweep", [])
    high_tm = [r for r in sweep if abs(r["true_mean"] - 0.300) < 1e-6]
    if high_tm:
        bm = high_tm[0]["by_method"]
        biases = [bm[m]["mean_slope"] for m in methods
                  if bm[m]["mean_slope"] is not None]
        # Under ~0% censoring we expect all methods to agree on bias
        # (≈ 0). Use a generous tolerance of 0.001 mg/L/yr.
        worst = max(abs(b) for b in biases) if biases else float("inf")
        H.check(
            "Falsification: at true_mean = 0.3 mg/L (near-zero censoring) "
            "all methods have |bias| ≤ 1e-3",
            worst <= 1e-3,
            f"worst |bias| = {worst:.4e}"
        )
    else:
        H.check("Falsification regime (true_mean=0.3) is present in sweep",
                False, "regime not found in monte_carlo_regime_sweep")

    # Check 16 (sensitivity stability): changing MIN_YEARS_PER_SITE should
    # not reverse the sign of the KM median slope — the headline qualitative
    # finding is expected to be robust to this choice.
    miny = R["sensitivity"]["min_years"]
    signs = [(1 if v["median_slope"] is not None and v["median_slope"] > 0
              else (-1 if v["median_slope"] is not None
                    and v["median_slope"] < 0 else 0))
             for v in miny]
    nonzero = [s for s in signs if s != 0]
    all_agree = (len(nonzero) == 0) or all(s == nonzero[0] for s in nonzero)
    H.check(
        "Sensitivity robustness: KM median-slope sign is stable across "
        "min_years ∈ {10,15,20}",
        all_agree,
        f"signs={signs}"
    )

    # Check 17 (diagnostic — informational): report censoring presence.
    any_censoring = any(d["censoring_rate"] > 0
                        for d in R["per_site_diagnostics"].values())
    print(f"  (info) Any site with censoring: {any_censoring}", flush=True)

    # Check 18 (data row count): every retained site must have > 0 total rows
    # parsed from its WQP CSV. A zero-row retained site would mean a parsing
    # bug promoted an empty CSV through quality control.
    zero_row_sites = [s for s, d in R["per_site_diagnostics"].items()
                      if d["total_rows"] <= 0]
    H.check(
        "Every retained site has > 0 parsed rows",
        not zero_row_sites,
        f"zero-row sites: {zero_row_sites[:3]}"
    )

    # Check 19 (corpus size): aggregate row count across retained sites must
    # exceed a loose floor (25 sites × 4 obs/year × 15 years = 1500 is the
    # structural minimum given MIN_OBS_PER_YEAR × MIN_YEARS_PER_SITE × 25;
    # we use a relaxed floor of 500 to allow a minority of sites being sparse).
    total_rows = sum(d["total_rows"]
                     for d in R["per_site_diagnostics"].values())
    H.check(
        f"Aggregate row count across retained sites ≥ 500 (got {total_rows})",
        total_rows >= 500,
        f"total_rows={total_rows}"
    )

    # Check 20 (CI width sanity — lower bound): at least one method's
    # bootstrap 95% CI width must exceed 1% of the |point estimate|. A CI
    # that is numerically collapsed to its point would indicate a degenerate
    # bootstrap (all slopes identical → no resampling variation).
    ci_width_ok = False
    for m in methods:
        b = R["by_method"][m]["bootstrap"]
        if b is None:
            continue
        lo, hi = b["median_ci95"]
        pt = b["median_point"]
        width = hi - lo
        if abs(pt) > 0 and width >= 0.01 * abs(pt):
            ci_width_ok = True
            break
        if abs(pt) == 0 and width > 0:
            # Point estimate exactly zero but CI is non-degenerate is fine.
            ci_width_ok = True
            break
    H.check(
        "CI-width sanity: ≥ 1 method has bootstrap CI width ≥ 1% of |estimate|",
        ci_width_ok
    )

    # Check 21 (CI width sanity — upper bound): no bootstrap CI width should
    # exceed MAX_PLAUSIBLE_SLOPE — a CI wider than the physical plausibility
    # bound indicates pathological bootstrap inputs.
    wild_ci = []
    for m in methods:
        b = R["by_method"][m]["bootstrap"]
        if b is None:
            continue
        lo, hi = b["median_ci95"]
        if (hi - lo) > MAX_PLAUSIBLE_SLOPE:
            wild_ci.append((m, hi - lo))
    H.check(
        f"No bootstrap CI width exceeds {MAX_PLAUSIBLE_SLOPE} mg/L/yr",
        not wild_ci,
        f"wild: {wild_ci[:3]}"
    )

    # Check 22 (MC variance): every method with at least one usable replicate
    # must have slope_sd ≥ 0 and rmse_slope ≥ |mean_slope| (RMSE always
    # dominates |bias| if replicates exist; a violation indicates an
    # aggregation bug).
    mc_variance_ok = True
    for m in methods:
        d = R["monte_carlo_null"]["by_method"][m]
        if d["mean_slope"] is None:
            continue
        if d["slope_sd"] is None or d["slope_sd"] < 0:
            mc_variance_ok = False
            break
        if d["rmse_slope"] is None or d["rmse_slope"] + 1e-12 < abs(
                d["mean_slope"]):
            mc_variance_ok = False
            break
    H.check(
        "MC variance sanity: slope_sd ≥ 0 and RMSE ≥ |bias| per method",
        mc_variance_ok
    )

    # Check 23 (Cohen's d on MC bias): the MC-bias signal for naive_zero vs
    # ros_mean must be bounded (|Cohen's d| on bias differences < 5), else a
    # rogue outlier has dominated and something is broken. We compute a
    # rough standardized effect size from the MC replicate moments.
    d_zero = R["monte_carlo_null"]["by_method"]["naive_zero"]
    d_ros = R["monte_carlo_null"]["by_method"]["ros_mean"]
    effect_ok = True
    effect_val = None
    if (d_zero["mean_slope"] is not None and d_ros["mean_slope"] is not None
            and d_zero["slope_sd"] not in (None, 0.0)
            and d_ros["slope_sd"] not in (None, 0.0)):
        pooled_sd = math.sqrt(0.5 * (d_zero["slope_sd"] ** 2
                                     + d_ros["slope_sd"] ** 2))
        if pooled_sd > 0:
            effect_val = (d_zero["mean_slope"] - d_ros["mean_slope"]) / pooled_sd
            effect_ok = abs(effect_val) < 5.0
    H.check(
        "Effect-size plausibility: |Cohen's d| on MC bias (naive_zero vs "
        "ros_mean) < 5",
        effect_ok,
        f"d={effect_val}"
    )

    # Check 24 (sensitivity — three settings of DL convention): the sensitivity
    # section must cover at least three DL-substitution conventions.
    dl_conv = R["sensitivity"]["dl_convention"]
    dl_conv_ok = len([m for m in ("naive_dl_half", "naive_dl", "naive_zero")
                      if m in dl_conv and dl_conv[m]["median_slope"] is not None]
                     ) >= 3
    H.check(
        "Sensitivity covers ≥ 3 DL-substitution conventions (DL/2, DL, zero)",
        dl_conv_ok
    )

    # Check 25 (seed propagation): every top-level seeded routine must have
    # been driven by SEED=42 — recorded in meta.
    seed_ok = R["meta"]["seed"] == SEED == 42
    H.check(
        f"Seed = 42 recorded in results.meta and propagated",
        seed_ok,
        f"meta.seed={R['meta']['seed']}"
    )

    # Check 26 (negative/falsification control — all-positive slopes):
    # the naive_zero method under declining DL should NEVER produce a strongly
    # negative MC bias. If it does, the DL schedule is mis-oriented.
    nz_bias = R["monte_carlo_null"]["by_method"]["naive_zero"]["mean_slope"]
    H.check(
        "Falsification: naive_zero MC bias is not strongly negative "
        "(> -1e-6 required)",
        nz_bias is None or nz_bias > -1e-6,
        f"nz_bias={nz_bias}"
    )

    # ── Summary ────────────────────────────────────────────────
    print(f"\n[SUMMARY] {H.passed} passed, {H.failed} failed",
          flush=True)
    if H.failed > 0:
        print("FAILED CHECKS:", flush=True)
        for m in H.messages:
            print(f"  {m}", flush=True)
        sys.exit(1)
    print("ALL CHECKS PASSED", flush=True)
    return True


# ────────────────────── MAIN ──────────────────────

def main():
    ap = argparse.ArgumentParser(
        description=("US Stream Nitrate Trends Under Detection-Limit "
                     "Censoring. Python 3.8+ stdlib only."))
    ap.add_argument("--verify", action="store_true",
                    help="Run verification assertions against "
                         "existing results.json.")
    args = ap.parse_args()

    random.seed(SEED)

    if args.verify:
        try:
            verify()
        except SystemExit:
            raise
        except Exception as e:
            sys.stderr.write(
                f"FATAL during --verify: {type(e).__name__}: {e}\n"
                "  (results.json may be missing or malformed — re-run the "
                "analysis to regenerate it.)\n"
            )
            sys.exit(4)
        return

    # Top-level error handling: catch any unexpected failure, print a clear
    # diagnostic to stderr, and exit with a nonzero code so the caller can
    # distinguish success from failure without parsing stdout.
    try:
        print("[1/6] Downloading and parsing USGS WQP nitrate records.",
              flush=True)
        try:
            per_site, fetch_diag = load_data()
        except Exception as e:
            sys.stderr.write(
                f"FATAL: load_data failed with {type(e).__name__}: {e}\n"
                "  Likely causes: no network access to "
                f"{WQP_BASE} (firewall / offline machine), WQP outage, or "
                "disk write failure in the cache directory. Check network, "
                f"then re-run. (Cache is at {CACHE_DIR}.)\n"
            )
            sys.exit(2)
        print(f"  Retained {len(per_site)} / {len(SITES)} sites with ≥1 row.",
              flush=True)
        if len(per_site) == 0:
            sys.stderr.write(
                "FATAL: no sites retained — every WQP fetch failed or every "
                "response was empty. Network access or WQP outage is the most "
                "likely cause. See fetch_diagnostics.sha256_by_site for "
                "per-site failure patterns. Aborting before downstream "
                "stages.\n"
            )
            sys.exit(3)

        print("[2/6] Computing per-site annual summaries and trends.",
              flush=True)
        results = run_analysis(per_site)

        print("[3/6] Running Monte Carlo null model "
              f"({N_MONTE_CARLO} replicates).", flush=True)
        mc = monte_carlo_null()
        print("  Running Monte Carlo regime sweep "
              "(true_mean ∈ {0.010, 0.030, 0.060, 0.100, 0.300}).", flush=True)
        mc_sweep = monte_carlo_regime_sweep()

        print("[4/6] Running sensitivity analyses.", flush=True)
        sens = {
            "dl_convention": sensitivity_dl_convention(per_site),
            "min_years": [
                sensitivity_min_years(per_site, 10),
                sensitivity_min_years(per_site, 15),
                sensitivity_min_years(per_site, 20),
            ],
            "period": [
                sensitivity_period(per_site, 1980, 2000),
                sensitivity_period(per_site, 1990, 2020),
                sensitivity_period(per_site, 2000, 2020),
            ],
        }

        print("[5/6] Writing results.json and report.md.", flush=True)
        try:
            generate_report(results, fetch_diag, mc, sens, mc_sweep)
        except OSError as e:
            sys.stderr.write(
                f"FATAL: could not write outputs to {OUT_DIR}: {e}\n"
                "  Check disk space and directory permissions.\n"
            )
            sys.exit(5)

        # Brief stdout summary
        print("[6/6] Summary:", flush=True)
        for m in METHODS:
            s = results["by_method"][m]
            if s["bootstrap"] is None:
                print(f"  {m}: no sites retained.", flush=True)
                continue
            print(
                f"  {m}: n={s['n_sites']}, "
                f"median slope = {s['bootstrap']['median_point']:+.5f} "
                f"mg/L/yr, frac sig. neg. = {s['frac_sig_negative']:.3f}",
                flush=True,
            )
        print("  MC null (true mean constant, declining DL):", flush=True)
        for m in METHODS:
            d = mc["by_method"][m]
            if d["mean_slope"] is None:
                print(f"    {m}: no replicates usable.", flush=True)
            else:
                print(
                    f"    {m}: mean slope = {d['mean_slope']:+.6f} mg/L/yr "
                    f"(bias), rejection rate = "
                    f"{d['rejection_rate_alpha05']:.3f}",
                    flush=True,
                )

        print("ANALYSIS COMPLETE", flush=True)
    except SystemExit:
        raise
    except Exception as e:
        sys.stderr.write(
            f"FATAL unexpected error: {type(e).__name__}: {e}\n"
        )
        import traceback
        traceback.print_exc(file=sys.stderr)
        sys.exit(10)


if __name__ == "__main__":
    main()
SCRIPT_EOF
```

**Expected output:** No output. The heredoc writes the Python script to `/tmp/claw4s_auto_us-stream-nitrate-trends-under-detection-limit-censoring/analysis.py`.

**Success:** The file `analysis.py` exists at that path and is a valid Python 3 source file (about 900 lines).

**Failure:** `cat` cannot create the file (disk full, permission denied). Fix the underlying filesystem issue.

---

## Step 3: Run analysis

```bash
cd /tmp/claw4s_auto_us-stream-nitrate-trends-under-detection-limit-censoring && python3 analysis.py
```

**Expected output (cold run, ~5–30 min):**

```
[1/6] Downloading and parsing USGS WQP nitrate records.
  [1/25] Fetching site USGS-01578310 ...
  [2/25] Fetching site USGS-01646500 ...
  ...
  [25/25] Fetching site USGS-14246900 ...
  Retained N / 25 sites with ≥1 row.
[2/6] Computing per-site annual summaries and trends.
[3/6] Running Monte Carlo null model (2000 replicates).
[4/6] Running sensitivity analyses.
[5/6] Writing results.json and report.md.
[6/6] Summary:
  naive_dl_half: n=K, median slope = -0.00??? mg/L/yr, frac sig. neg. = 0.???
  naive_dl: ...
  naive_zero: ...
  km_mean: ...
  ros_mean: ...
  MC null (true mean constant, declining DL):
    naive_dl_half: mean slope = -0.000??? mg/L/yr (bias), rejection rate = 0.???
    ...
ANALYSIS COMPLETE
```

**Expected outputs on disk:**

- `results.json` (structured results, ~50–500 KB)
- `report.md` (human-readable summary)
- `cache/site_XXXXXXXX.csv` (one per fetched site)
- `cache/hash_index.json` (SHA-256 of every cached CSV)

**Success:** stdout ends with `ANALYSIS COMPLETE`. Both `results.json` and `report.md` exist.

**Failure:** If any site's fetch fails after 4 retries, the site is skipped and recorded in `fetch_diagnostics`; the run continues. If fewer than 8 sites survive quality control the verification step will fail.

---

## Step 4: Verify results

```bash
cd /tmp/claw4s_auto_us-stream-nitrate-trends-under-detection-limit-censoring && python3 analysis.py --verify
```

**Expected output:**

```
[PASS] KM hand-checked unit test matches reference (1.400)
[PASS] results.json exists and parses as JSON
[PASS] ≥ 8 sites retained (got 17)
[PASS] All 5 methods present in results
[PASS] Every method has a bootstrap block
[PASS] Monte Carlo null present with 2000 replicates
[PASS] MC produced ≥ 1 usable replicate for every method
[PASS] MC: bias(naive_zero) ≥ 0 under declining DL (required sign)
[PASS] MC: bias(naive_dl) ≤ 0 under declining DL (required sign)
[PASS] MC: best censored-aware method beats worst naive on |bias|
[PASS] per_site_diagnostics count equals n_sites_retained
[PASS] All three sensitivity analyses present with ≥ 2 rows each
[PASS] At least one method has bootstrap CI bracketing its point estimate
[PASS] SHA-256 fingerprint recorded for every retained site
[PASS] All per-site slopes within ±1.0 mg/L/yr (physical plausibility)
[PASS] At least one method has a bootstrap 95% CI of non-zero width
[PASS] Falsification: at true_mean = 0.3 mg/L (near-zero censoring) all methods have |bias| ≤ 1e-3
[PASS] Sensitivity robustness: KM median-slope sign is stable across min_years ∈ {10,15,20}
  (info) Any site with censoring: ...
[PASS] Every retained site has > 0 parsed rows
[PASS] Aggregate row count across retained sites ≥ 500 (got ...)
[PASS] CI-width sanity: ≥ 1 method has bootstrap CI width ≥ 1% of |estimate|
[PASS] No bootstrap CI width exceeds 1.0 mg/L/yr
[PASS] MC variance sanity: slope_sd ≥ 0 and RMSE ≥ |bias| per method
[PASS] Effect-size plausibility: |Cohen's d| on MC bias (naive_zero vs ros_mean) < 5
[PASS] Sensitivity covers ≥ 3 DL-substitution conventions (DL/2, DL, zero)
[PASS] Seed = 42 recorded in results.meta and propagated
[PASS] Falsification: naive_zero MC bias is not strongly negative (> -1e-6 required)

[SUMMARY] 25 passed, 0 failed
ALL CHECKS PASSED
```

**Success criteria (all assertions must pass):**

1. `results.json` exists and parses as JSON.
2. ≥ 8 sites retained.
3. All five methods (`naive_dl_half`, `naive_dl`, `naive_zero`, `km_mean`, `ros_mean`) are present with bootstrap output.
4. The Monte Carlo section contains all five methods with at least one usable replicate each.
5. Under the Monte Carlo null (constant mean, log-linearly declining DL), at least one censored-data-aware method (`km_mean` or `ros_mean`) has a smaller-magnitude mean-slope bias than *both* of the worst-case naive methods (`naive_dl`, `naive_zero`). This is the key methodological claim: a censored-aware estimator can outperform the extremal naive conventions on bias under a stable-mean, declining-DL regime.
6. `bias(naive_zero)` ≥ 0 (naive-zero substitution produces a non-negative mean slope under declining DL because early-year censored rows are replaced with 0, biasing the early-year annual mean low).
7. `bias(naive_dl)` ≤ 0 (naive-DL substitution produces a non-positive mean slope because early-year censored rows are replaced with the full DL, biasing the early-year annual mean high).
8. A built-in KM unit-test passes: given a small hand-checked synthetic sample, `km_mean_from_left_censored()` returns the expected value to four decimal places.
9. `per_site_diagnostics` size equals the retained site count.
10. All three sensitivity sections (`dl_convention`, `min_years`, `period`) exist with ≥ 2 entries each.
11. At least one method has a bootstrap CI bracketing its point estimate.
12. Every retained site has a recorded SHA-256 fingerprint of its WQP CSV.

**Failure conditions:**

- *"Only N sites retained; floor = 8"*: Too many WQP queries failed. Inspect `fetch_diagnostics.sha256_by_site` for `null` entries. Re-running may help (network-transient errors).
- *"MC sanity: the best censored-aware method should have smaller magnitude bias"*: a bug in either the KM or ROS implementation, or an MC regime in which all methods degenerate. Inspect `monte_carlo_null.by_method` in `results.json`.
- *"MC sanity: naive_zero bias should be ≥ 0 ..."* / *"... naive_dl bias should be ≤ 0 ..."*: check the MC DL schedule (`MC_DL_START`, `MC_DL_END`); these assertions assume a DL that declines over time.
- *"KM unit test failed"*: the KM-mean implementation has drifted from the expected reference output. Inspect `km_mean_from_left_censored()`.
- *"Missing SHA256 for site S"*: fetch failed mid-run. Delete `cache/site_S.csv` and re-run Step 3.

---

## Notes on reproducibility

- **Deterministic computation and byte-stable outputs.** Every random operation uses `random.seed(SEED)` with `SEED = 42`. `results.json` contains no wall-clock timestamp and no other run-dependent fields; re-running the script on the same on-disk cache yields a *byte-identical* `results.json` (and `report.md`). The wall-clock time of the run is printed to stdout and captured in `execution_log.txt`, but is not embedded in the analytic output. Bootstrap and Monte Carlo routines each use their own independently seeded `random.Random` instance derived from `SEED`, so the order of method evaluation does not perturb the results.
- **Local cache integrity.** Every WQP CSV is SHA-256 fingerprinted on first download and the fingerprint is recorded in `cache/hash_index.json`. Subsequent re-runs verify the cached bytes against the stored hash before parsing; a drift triggers a re-fetch and a hash-index update. This guarantees local cache integrity but does *not* pin against upstream WQP data. If WQP re-curates or corrects historical records, a fresh run with an empty cache may see different numbers. The `hash_index.json` persisted alongside the cache serves as a run-specific data-manifest: any consumer of the results can re-verify that the exact bytes used in the analysis are still on disk.
- **Upstream data provenance.** Numeric results depend on the WQP's internal database state at first-download time. For downstream regulatory or publication use, archive the `cache/` directory (≈5–50 MB) alongside the `results.json`; together they constitute a self-contained reproducibility package that does not require network access to re-analyze.
- **No third-party Python packages** (stdlib only): `urllib`, `csv`, `json`, `hashlib`, `math`, `random`, `statistics`, `datetime`, `argparse`, `os`, `sys`, `time`, `collections`, `io`.