---
name: "USA-NPN Observer-Tenure Confound on Spring Advance"
description: "Tests whether the commonly reported ~2-4-days-per-decade first-leaf advance in the USA-National Phenology Network's Nature's Notebook dataset survives restriction to individual plants that have been continuously monitored at the same site for at least five consecutive years (individual-tenure, used here as a publicly available proxy for observer-tenure because USA-NPN's public API does not expose observer identifiers). Downloads real phenometric data for six widely monitored deciduous species from the USA-NPN services.usanpn.org API, estimates a pooled within-individual fixed-effects temporal trend on the full and tenure-matched cohorts, and evaluates the cohort difference with a bootstrap 95% CI and a 1,000-permutation label-shuffle null model. All HTTP responses are SHA-256 verified."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "climate-phenology", "first-leaf", "usa-npn", "natures-notebook", "observer-tenure", "permutation-test", "fixed-effects", "bootstrap", "spring-advance"]
python_version: ">=3.8"
dependencies: []
data_source: "USA-National Phenology Network Nature's Notebook, getSummarizedData endpoint (https://services.usanpn.org/npn_portal/observations/getSummarizedData.csv), queried per (species_id, phenophase_id) pair for the 2009-2023 window. Per-species responses are byte-hashed (SHA-256) after first download and the hash index is stored alongside the cache."
data_revision: "USA-NPN services.usanpn.org live getSummarizedData.csv endpoint; query parameters pinned in the DOMAIN CONFIGURATION block of the analysis script."
---

# USA-NPN Observer-Tenure Confound on Spring Advance

## Research Question

Does the widely cited ~2-4-days-per-decade first-leaf advance estimated from the USA-National Phenology Network's Nature's Notebook dataset survive restriction to the sub-cohort of individual plants that were continuously monitored at the same site for at least five consecutive years (a publicly auditable proxy for long-tenured observers), and is the difference between the tenure-matched and short-tenure slopes distinguishable from label-shuffle noise at the 5% level?

## Success Criteria

The analysis is considered successful when *all* of the following measurable conditions hold:

1. The analysis script runs to completion and stdout ends with the literal string `ANALYSIS COMPLETE`.
2. `results.json` and `report.md` are written to the workspace and are valid JSON / UTF-8 Markdown, respectively.
3. `python3 analysis.py --verify` runs to completion, ends with the literal string `ALL CHECKS PASSED`, and every one of its machine-checkable assertions (see Step 4) returns true.
4. `results.json` contains a non-null `slope_full_cohort.point` with a 95% bootstrap CI (`ci95_lo`, `ci95_hi`) bracketing the point estimate.
5. `results.json` contains a non-null `slope_tenure_matched.point` with a 95% bootstrap CI bracketing the point estimate.
6. `results.json` contains a `permutation_null` block with a two-sided p-value in `[0, 1]` computed from at least 900 successful permutations.
7. `results.json` contains a `sensitivity` block with at least three tenure-threshold rows, at least four species rows, and at least two time-window rows.
8. Every fetched species CSV is recorded in `fetch_diagnostics.sha256_by_species` with a valid hex SHA-256 fingerprint.
9. All point slopes lie within the physically plausible range `[-5, +5]` days per year.
10. Re-running the analysis on the cached data produces a byte-identical `results.json` modulo the `timestamp_utc` field, i.e. `meta.results_hash` (which excludes the timestamp by construction) is stable across re-runs.

## Failure Conditions

The analysis must be treated as failed (and the script / skill iterated on) if any of the following occur:

- stdout does not end with `ANALYSIS COMPLETE`, or `--verify` stdout does not end with `ALL CHECKS PASSED`.
- Any `--verify` assertion raises (traceback printed to stderr, exit code 2 per the top-level handler).
- The USA-NPN fetch returned zero species CSVs (exit code 1 with a `FATAL: analysis aborted` message on stderr, no partial `results.json` written).
- Fewer than `MIN_INDIVIDUALS_FLOOR` (=500) distinct individuals or fewer than `MIN_TENURE_MATCHED_FLOOR` (=30) tenure-matched individuals are retained after filtering.
- Any recorded per-species SHA-256 is `null` (i.e. the USA-NPN fetch failed after all retries), *and* this reduces the fetched-species count below `MIN_SPECIES_FETCHED` (=4) -- the analysis is not trustworthy below four species.
- A point slope falls outside `[SLOPE_PLAUSIBLE_MIN, SLOPE_PLAUSIBLE_MAX]` = `[-5, +5]` days per year (strong signal of a parser bug or sentinel leakage).
- The 95% bootstrap CI for the full-cohort slope is narrower than `CI_WIDTH_MIN` (=0.02) days/year (implausibly tight -- likely single-individual degeneracy) or wider than `CI_WIDTH_MAX` (=2.0) days/year (no signal).
- The negative-control (within-individual year<->DOY shuffle) slope is not strictly smaller in absolute value than the full-cohort slope -- the within-individual FE estimator is leaking between-individual variance.
- The label-shuffle null mean is not centred on zero (`|mean_null| / sd_null >= 0.5`) -- exchangeability broke.
- `meta.results_hash` differs between two back-to-back runs on the same cache (non-determinism, seed not plumbed through).

**Script exit codes (top-level handler):** `0` success, `1` data / network / parse error, `2` verification assertion failed, `130` Ctrl+C interrupt.

## When to Use This Skill

Use this skill when you need to investigate whether a reported temporal trend in volunteer-collected phenology data -- specifically the widely cited 2-4-days-per-decade advance in North American first-leaf dates from the USA-NPN Nature's Notebook network -- is robust to the confound that the observer network grew rapidly after its 2009 launch, so that cross-year comparisons on the full dataset mix long-tenured observers (continuous site records) with new sign-ups whose inclusion over time can shift the naive trend independent of any underlying phenological change.

### Preconditions

- **Python version:** 3.8+ (standard library only -- no numpy, scipy, pandas, or requests required).
- **Network:** Internet access required on first run to fetch the six USA-NPN per-species CSVs from `services.usanpn.org`. Subsequent runs re-use the on-disk SHA-256-verified cache; re-runs and `--verify` execute in well under 60 seconds.
- **Runtime:** Cold run with an empty cache is approximately 1-5 minutes (USA-NPN server response time dominates; the statistical computation itself runs in well under a minute on stdlib). Cached re-runs: < 1 minute.
- **Disk:** Per-run cache is ~5 MB (one CSV per species, SHA-256 indexed).
- **Politeness:** Requests use a descriptive `User-Agent` identifying the pipeline and space sequential fetches by half a second.

## Adaptation Guidance

The statistical core of this skill -- pooled within-individual fixed-effects regression of day-of-year on year, tenure classification via longest consecutive-year run per unit, bootstrap CIs over resampled units, and a label-permutation null on the tenure-vs-complement slope difference -- is domain-agnostic. It can be re-pointed at any question of the form *"does a reported temporal trend in a volunteer-monitored repeated-measures dataset survive restriction to long-tenured units?"*

- **Change the data source:** Edit the `DOMAIN CONFIGURATION` block near the top of the script. `USA_NPN_BASE` is the fixed URL of the getSummarizedData endpoint. `REQUEST_SRC` is the identifier USA-NPN uses for request-log attribution.
- **Change the species / phenophase list:** `SPECIES` is a list of `(species_id, phenophase_id, scientific_name, common_name)` tuples. Phenophase 371 is the standardized "Breaking leaf buds" protocol; 373 is the "Breaking leaf buds (lilac/honeysuckle)" protocol. Species IDs are USA-NPN internal numbers and are listed at `https://services.usanpn.org/npn_portal/species/getSpecies.json`.
- **Change the time window:** Edit `START_YEAR` and `END_YEAR`. The query window for each year is `YYYY-01-01` through `YYYY-07-01` (spring-only) so that `First_Yes_DOY` within the window equals the annual first-leaf day.
- **Change the tenure operationalization:** `longest_consecutive_run()` returns the longest run of consecutive years for any unit. `TENURE_THRESHOLD` (default 5) is the primary cut-off; `TENURE_SENSITIVITY` lists alternate thresholds reported in the sensitivity section.
- **Change the trend estimator:** `within_individual_pooled_slope()` is a fixed-effects OLS estimator, and `theil_sen_slope()` / `per_individual_slopes()` provide a distribution-free robustness alternative. Both are used in the output.
- **Change the null model:** `permutation_null_tenure_difference()` shuffles the tenure-matched flag across eligible units. Edit the shuffle scheme (for example, stratified-by-species shuffles) without changing the rest of the pipeline.
- **What stays the same:** HTTP fetch and SHA-256 caching (`fetch_with_cache()`, `http_get()`), CSV parsing (`parse_npn_csv()`), bootstrap resampling (`bootstrap_slope_over_individuals()`, `_bootstrap_delta()`), the `--verify` assertion harness, and the `results.json` / `report.md` writers are all general-purpose.

## Overview

**The claim under test.** Analyses of the USA-NPN Nature's Notebook network routinely report that first-leaf dates in North American deciduous species have advanced by roughly 2-4 days per decade since the network's national launch in 2009. This figure is then taken as evidence that spring is arriving earlier in response to climate warming.

**The confound.** The Nature's Notebook observer corps grew from a few thousand reports per year in 2009 to tens of thousands by the late 2010s. Because new sign-ups bring new sites (often at warmer southern latitudes, in residential yards rather than reference forests, with species and individuals that were not previously monitored), a naive cross-cohort trend estimated on the pooled dataset mixes a within-individual temporal signal with a between-individual composition shift over time. If, for example, later years include disproportionately more southern or urban individuals with earlier baseline DOYs, the pooled trend will understate (or overstate) the within-individual advance.

**What this skill does.** Downloads per-species `getSummarizedData` CSVs for six widely monitored USA-NPN deciduous species (common lilac, red maple, sugar maple, flowering dogwood, eastern redbud, border forsythia) across 2009-2023, consolidates each (individual plant, year) row to its earliest first-yes DOY, computes the longest consecutive-year observation run per individual, and estimates the pooled within-individual fixed-effects slope of first-leaf DOY on year (a) for the full cohort of individuals with at least three observation years and (b) for the tenure-matched sub-cohort whose longest consecutive run is at least five years. The within-individual fixed-effects estimator removes between-individual composition drift by construction; the tenure restriction additionally requires that contributing individuals have long, regular records at a single site.

**Null model and effect sizes.** (a) **Bootstrap over individuals** -- individuals are resampled with replacement 1,000 times and the within-individual slope (or cohort difference) is recomputed on each bootstrap sample; 95% CIs are reported as the 2.5-97.5 percentile range. (b) **Label-permutation null of tenure membership** -- the tenure-matched flag is randomly re-assigned across eligible individuals (preserving the observed count of tenure-matched), the slope difference is recomputed on 1,000 permutations, and the two-sided p-value is returned using the Phipson-Smyth (+1)/(n+1) adjustment.

**Methodological hook.** *Tenure-matched cohort combined with a permutation null on the tenure-vs-complement slope difference.* Prior USA-NPN trend analyses either pool everything or restrict to a single species; neither controls for observer-network growth. We (i) construct a tenure-matched cohort defined by longest consecutive-year run per individual, (ii) report a full-vs-tenure trend comparison with bootstrap CIs, and (iii) calibrate the significance of that comparison against an explicit permutation null so readers can see what difference would be expected by chance.

**What this is not.** Not a claim about the cause of any first-leaf advance -- the analysis is descriptive and agnostic to climate attribution. Not a representative census of US phenology -- the USA-NPN network is volunteer-driven and over-samples certain regions and species. Not an observer-level analysis -- observer identifiers are not returned by the public API, so "individual-tenure" (continuous monitoring of a physical plant) is used as the finest-grained public proxy. Not a geographic analysis -- we do not stratify by latitude or climate region in the primary test; doing so is a straightforward extension via the sensitivity scaffolding.

---

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_usa-npn-observer-tenure-confound-on-spring-advance/cache
```

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

**Success:** The directory `/tmp/claw4s_auto_usa-npn-observer-tenure-confound-on-spring-advance/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_usa-npn-observer-tenure-confound-on-spring-advance/analysis.py
#!/usr/bin/env python3
"""
USA-NPN Observer-Tenure Confound on Spring Advance.

Tests whether the ~2-4-days-per-decade first-leaf advance reported from the
USA-National Phenology Network's Nature's Notebook dataset survives restriction
to individual plants that have been continuously monitored (a proxy for long-
tenured observers) for at least five consecutive years at the same site. The
confound being tested is that the observer network grew rapidly after 2009 --
new sign-ups bring new sites, new species, and new baseline phenologies whose
inclusion over time can shift the naive cross-cohort trend even when the
underlying phenology of any given monitored individual is stable. The skill
computes a pooled within-individual fixed-effects trend on the full dataset,
repeats it on the tenure-matched cohort, and evaluates whether the difference
between the two trends is significant under (a) individual-resampling
bootstrap and (b) a label-permutation null model in which the "tenure-matched"
flag is shuffled across individuals.

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.
# ═══════════════════════════════════════════════════════════════════════

USA_NPN_BASE = (
    "https://services.usanpn.org/npn_portal/observations/getSummarizedData.csv"
)
USER_AGENT = (
    "Claw4S-USANPN-TenureConfound/1.0 "
    "(Claw4S 2026 reproducibility pipeline; "
    "contact: claw4s-pipeline@example.org)"
)
REQUEST_SRC = "claw4s_tenure_confound"

# Six widely monitored deciduous species in USA-NPN Nature's Notebook, each
# paired with the standardized first-leaf phenophase_id that applies to it.
# The "lilac/honeysuckle" species (syringa, lonicera clones) use protocol 373;
# all other woody deciduous species use the generic "Breaking leaf buds" = 371.
# (species_id, phenophase_id, scientific_name, common_name)
SPECIES = [
    (36,  373, "Syringa vulgaris",        "common lilac"),
    (3,   371, "Acer rubrum",              "red maple"),
    (81,  371, "Acer saccharum",           "sugar maple"),
    (97,  371, "Cornus florida",           "flowering dogwood"),
    (68,  371, "Cercis canadensis",        "eastern redbud"),
    (12,  371, "Forsythia x intermedia",   "border forsythia"),
]

START_YEAR = 2009
END_YEAR = 2023
# Query spring-only window per year (Jan 1 through July 1) so that
# First_Yes_DOY within the window equals the annual first-leaf DOY.
QUERY_MONTH_DAY_START = "01-01"
QUERY_MONTH_DAY_END = "07-01"

MIN_YEARS_PER_INDIVIDUAL = 3    # individuals need >= this many annual first-yes values
TENURE_THRESHOLD = 5            # primary "tenure-matched cohort" = >= this many consecutive years
TENURE_SENSITIVITY = [3, 5, 7, 10]

N_BOOTSTRAP = 1000              # resamples for 95% CI on any pooled slope
N_PERMUTATIONS = 1000           # label-shuffle iterations for null model
SEED = 42                       # master seed -- every RNG is Random(SEED + k)
CI_LEVEL = 95.0                 # percentile-based bootstrap CI level (%)
CI_ALPHA_LO = (100.0 - CI_LEVEL) / 2.0   # = 2.5
CI_ALPHA_HI = 100.0 - CI_ALPHA_LO        # = 97.5
SIGNIFICANCE_THRESHOLD = 0.05   # two-sided alpha for permutation p-value

# Quality floors enforced by --verify. Raise these to be stricter.
MIN_INDIVIDUALS_FLOOR = 500     # minimum total individuals retained
MIN_TENURE_MATCHED_FLOOR = 30   # minimum tenure-matched sub-cohort size
MIN_N_OBS_FLOOR = 1000          # minimum (individual, year) rows in eligible
MIN_SPECIES_FETCHED = 4         # below this, cohort analysis is untrustworthy
SLOPE_PLAUSIBLE_MIN = -5.0      # days/year floor -- anything lower is a bug
SLOPE_PLAUSIBLE_MAX = 5.0       # days/year ceiling -- anything higher is a bug
CI_WIDTH_MIN = 0.02             # days/year -- guards against degenerate bootstrap
CI_WIDTH_MAX = 2.0              # days/year -- guards against no-signal runs

HTTP_TIMEOUT = 180              # seconds per request
HTTP_MAX_RETRIES = 4            # attempts before abandoning a species fetch
HTTP_RETRY_BACKOFF = 3.0        # seconds; doubled each retry
HTTP_SLEEP_BETWEEN = 0.5        # politeness delay between species fetches

# DOY sanity window for first-leaf (Jan 1 = 1, ..., Dec 31 = 365/366).
# Filter out obviously-wrong codings (e.g. -9999, non-spring leafing).
DOY_VALID_MIN = 1               # earliest plausible first-leaf DOY
DOY_VALID_MAX = 220             # latest plausible first-leaf DOY (Aug 8)
SENTINEL_NODATA = -9999         # USA-NPN's "missing" sentinel

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")

# Column indices (zero-based) for USA-NPN getSummarizedData.csv output.
COL_SITE_ID = 0
COL_LAT = 1
COL_LON = 2
COL_STATE = 4
COL_SPECIES_ID = 5
COL_INDIVIDUAL_ID = 10
COL_PHENOPHASE_ID = 11
COL_FIRST_YES_YEAR = 13
COL_FIRST_YES_DOY = 16

# ═══════════════════════════════════════════════════════════════════════
# 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):
    req = urllib.request.Request(url, headers={
        "User-Agent": USER_AGENT,
        "Accept": "text/csv",
    })
    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[:120]}...: {e}. "
                f"Retrying in {delay:.1f}s.\n"
            )
            time.sleep(delay)
    raise RuntimeError(f"HTTP fetch failed after retries: {url}: {last_err}")


def build_query_url(species_id, phenophase_id):
    params = [
        ("request_src", REQUEST_SRC),
        ("start_date", f"{START_YEAR}-{QUERY_MONTH_DAY_START}"),
        ("end_date", f"{END_YEAR}-{QUERY_MONTH_DAY_END}"),
        ("species_id[0]", str(species_id)),
        ("phenophase_id[0]", str(phenophase_id)),
    ]
    return USA_NPN_BASE + "?" + urllib.parse.urlencode(params)


def fetch_with_cache(species_id, phenophase_id):
    """Fetch USA-NPN summarized CSV for one (species, phenophase) pair.

    The resulting bytes are cached on disk and verified via SHA-256 on
    subsequent runs. Returns (bytes, sha256_hex).
    """
    os.makedirs(CACHE_DIR, exist_ok=True)
    cache_key = f"species_{species_id}_phenophase_{phenophase_id}.csv"
    cache_path = os.path.join(CACHE_DIR, cache_key)
    idx = load_hash_index()

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

    url = build_query_url(species_id, phenophase_id)
    data = http_get(url)
    h = sha256_bytes(data)
    with open(cache_path, "wb") as f:
        f.write(data)
    idx[cache_key] = h
    save_hash_index(idx)
    time.sleep(HTTP_SLEEP_BETWEEN)
    return data, h


# ---------------------------- PARSING ------------------------------------

def parse_npn_csv(csv_bytes, expected_species_id, expected_phenophase_id):
    """Parse USA-NPN summarized CSV bytes into a list of observation dicts.

    Each row of the CSV is one (individual, phenophase, first-yes event) tuple.
    The file has no header row. The parser discards rows with sentinel
    (-9999) or out-of-season DOY values and enforces that species_id and
    phenophase_id match the requested query.

    Returns list of dicts with keys:
        site_id (str), lat (float), lon (float), state (str),
        species_id (int), individual_id (str), phenophase_id (int),
        year (int), doy (int).
    """
    text = csv_bytes.decode("utf-8", errors="replace")
    reader = csv.reader(io.StringIO(text))
    rows = []
    for raw in reader:
        if len(raw) <= COL_FIRST_YES_DOY:
            continue
        try:
            species_id = int(raw[COL_SPECIES_ID])
            phenophase_id = int(raw[COL_PHENOPHASE_ID])
            year = int(raw[COL_FIRST_YES_YEAR])
            doy = int(raw[COL_FIRST_YES_DOY])
        except ValueError:
            continue
        if species_id != expected_species_id:
            continue
        if phenophase_id != expected_phenophase_id:
            continue
        if year < START_YEAR or year > END_YEAR:
            continue
        if doy == SENTINEL_NODATA:
            continue
        if doy < DOY_VALID_MIN or doy > DOY_VALID_MAX:
            continue
        site_id = (raw[COL_SITE_ID] or "").strip()
        individual_id = (raw[COL_INDIVIDUAL_ID] or "").strip()
        if not site_id or not individual_id:
            continue
        try:
            lat = float(raw[COL_LAT])
            lon = float(raw[COL_LON])
        except ValueError:
            lat = float("nan")
            lon = float("nan")
        state = (raw[COL_STATE] or "").strip()
        rows.append({
            "site_id": site_id,
            "lat": lat,
            "lon": lon,
            "state": state,
            "species_id": species_id,
            "individual_id": individual_id,
            "phenophase_id": phenophase_id,
            "year": year,
            "doy": doy,
        })
    return rows


def consolidate_per_individual_year(rows):
    """Collapse multiple same-(individual, year) rows to a single earliest DOY.

    Some individuals report status through yes-no-yes gaps within a spring
    season, producing multiple rows from the summarized endpoint. Reduce to
    one (individual, year) -> min(DOY) observation. Returns a dict keyed by
    individual_id with value {year: doy, ...} and a parallel dict giving
    metadata (site_id, species_id, lat, lon) per individual.
    """
    by_ind = defaultdict(dict)
    meta = {}
    for r in rows:
        ind = r["individual_id"]
        yr = r["year"]
        prev = by_ind[ind].get(yr)
        if prev is None or r["doy"] < prev:
            by_ind[ind][yr] = r["doy"]
        if ind not in meta:
            meta[ind] = {
                "site_id": r["site_id"],
                "species_id": r["species_id"],
                "state": r["state"],
                "lat": r["lat"],
                "lon": r["lon"],
                "phenophase_id": r["phenophase_id"],
            }
    return dict(by_ind), meta


# ---------------------------- TENURE CALCULUS ----------------------------

def longest_consecutive_run(years_set):
    """Given a set of integer years, return the length of the longest run
    of consecutive integers (e.g., {2010, 2011, 2012, 2015, 2016} -> 3)."""
    if not years_set:
        return 0
    ys = sorted(years_set)
    best = cur = 1
    for i in range(1, len(ys)):
        if ys[i] == ys[i - 1] + 1:
            cur += 1
            if cur > best:
                best = cur
        else:
            cur = 1
    return best


def classify_individuals(by_ind, tenure_threshold):
    """For each individual compute (n_years, longest_run, is_tenure_matched).

    Returns a dict {individual_id: {"n_years": int, "longest_run": int,
                                    "tenure_matched": bool}}.
    """
    out = {}
    for ind, year_doy in by_ind.items():
        years = set(year_doy.keys())
        lr = longest_consecutive_run(years)
        out[ind] = {
            "n_years": len(years),
            "longest_run": lr,
            "tenure_matched": lr >= tenure_threshold,
        }
    return out


# ---------------------------- REGRESSION ---------------------------------

def within_individual_pooled_slope(by_ind, include_ids=None,
                                   min_years=MIN_YEARS_PER_INDIVIDUAL):
    """Pooled within-individual fixed-effects slope of DOY on Year.

    For each included individual with at least `min_years` distinct years:
      demean DOY and Year by that individual's own means,
    then run OLS on the pooled demeaned (year_dev, doy_dev) sample:
      slope_hat = sum(year_dev_ij * doy_dev_ij) / sum(year_dev_ij^2).

    This estimator is mathematically equivalent to an OLS regression with
    individual fixed effects and identifies the average within-individual
    temporal trend (days per year). Multiply by 10 to report per-decade.
    Returns None if fewer than two usable individuals contribute or if the
    pooled year-variance is zero.
    """
    if include_ids is None:
        include_ids = set(by_ind.keys())
    # Sort for deterministic floating-point accumulation order across runs
    # (set/dict iteration order depends on PYTHONHASHSEED for str keys).
    num = 0.0   # sum of year_dev * doy_dev
    den = 0.0   # sum of year_dev^2
    n_ind = 0
    n_obs = 0
    for ind in sorted(include_ids):
        year_doy = by_ind.get(ind, {})
        if len(year_doy) < min_years:
            continue
        years = sorted(year_doy.keys())
        doys = [year_doy[y] for y in years]
        my = sum(years) / len(years)
        md = sum(doys) / len(doys)
        for y, d in zip(years, doys):
            ydev = y - my
            ddev = d - md
            num += ydev * ddev
            den += ydev * ydev
        n_ind += 1
        n_obs += len(years)
    if den <= 0.0 or n_ind < 2:
        return None
    slope = num / den
    return {"slope": slope, "n_individuals": n_ind, "n_obs": n_obs}


def theil_sen_slope(xs, ys):
    n = len(xs)
    if n < 2:
        return 0.0
    slopes = []
    for i in range(n - 1):
        for j in range(i + 1, n):
            dx = xs[j] - xs[i]
            if dx != 0:
                slopes.append((ys[j] - ys[i]) / dx)
    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])


def per_individual_slopes(by_ind, include_ids=None,
                          min_years=MIN_YEARS_PER_INDIVIDUAL):
    """Compute a Theil-Sen slope per individual (for robustness reporting)."""
    if include_ids is None:
        include_ids = set(by_ind.keys())
    slopes = []
    for ind in sorted(include_ids):
        year_doy = by_ind.get(ind, {})
        if len(year_doy) < min_years:
            continue
        years = sorted(year_doy.keys())
        doys = [year_doy[y] for y in years]
        slopes.append(theil_sen_slope(years, doys))
    return slopes


# ---------------------------- BOOTSTRAP ----------------------------------

def bootstrap_slope_over_individuals(by_ind, include_ids,
                                     n_boot=N_BOOTSTRAP,
                                     min_years=MIN_YEARS_PER_INDIVIDUAL,
                                     rng=None):
    """Bootstrap 95% CI on the pooled within-individual slope by resampling
    individuals with replacement. Returns dict with point, lo, hi.
    """
    if rng is None:
        rng = random.Random(SEED)
    # Sort for deterministic iteration order across PYTHONHASHSEED values.
    ids = sorted(i for i in include_ids
                 if len(by_ind.get(i, {})) >= min_years)
    if len(ids) < 2:
        return None
    # Point estimate
    point = within_individual_pooled_slope(by_ind, set(ids), min_years)
    if point is None:
        return None

    # Precompute per-individual contributions (num_i, den_i, n_obs_i) so
    # bootstrap only needs to sum these rather than re-iterate observations.
    contribs = {}
    for ind in ids:
        year_doy = by_ind[ind]
        years = sorted(year_doy.keys())
        doys = [year_doy[y] for y in years]
        my = sum(years) / len(years)
        md = sum(doys) / len(doys)
        num_i = 0.0
        den_i = 0.0
        for y, d in zip(years, doys):
            ydev = y - my
            ddev = d - md
            num_i += ydev * ddev
            den_i += ydev * ydev
        contribs[ind] = (num_i, den_i)

    slopes = []
    for _ in range(n_boot):
        num = 0.0
        den = 0.0
        for _ in range(len(ids)):
            pick = ids[rng.randrange(len(ids))]
            ni, di = contribs[pick]
            num += ni
            den += di
        if den > 0.0:
            slopes.append(num / den)
    if not slopes:
        return None
    slopes.sort()
    return {
        "point": point["slope"],
        "n_individuals": point["n_individuals"],
        "n_obs": point["n_obs"],
        "ci95_lo": percentile(slopes, CI_ALPHA_LO),
        "ci95_hi": percentile(slopes, CI_ALPHA_HI),
        "n_boot": len(slopes),
    }


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


# ---------------------------- PERMUTATION NULL ---------------------------

def permutation_null_tenure_difference(by_ind, classification,
                                       min_years=MIN_YEARS_PER_INDIVIDUAL,
                                       n_perm=N_PERMUTATIONS,
                                       rng=None):
    """Permutation test on the difference between tenure-matched and
    complementary (short-tenure) pooled slopes.

    Null hypothesis: tenure-match status is unrelated to an individual's
    within-individual trend. The statistic is

        Delta = slope_tenure_matched - slope_short_tenure

    where both slopes are pooled within-individual fixed-effects estimates
    computed only on individuals with at least `min_years` distinct years.
    Under the null, we shuffle the tenure-matched label across the eligible
    individuals (preserving the observed count of tenure-matched) and
    recompute Delta. The two-sided p-value is the fraction of permutations
    with |Delta_perm| >= |Delta_obs|.

    Returns dict with obs_delta, mean_null, sd_null, p_two_sided, n_perm.
    """
    if rng is None:
        rng = random.Random(SEED + 1)
    # Sort eligible for deterministic permutation order across runs.
    eligible = sorted(i for i in classification
                      if len(by_ind.get(i, {})) >= min_years)
    if len(eligible) < 4:
        return None
    labels = {i: classification[i]["tenure_matched"] for i in eligible}
    n_high = sum(1 for i in eligible if labels[i])
    n_low = len(eligible) - n_high
    if n_high < 2 or n_low < 2:
        return {
            "error": ("Too few individuals in tenure-matched or short-"
                      "tenure group under the given threshold."),
            "n_eligible": len(eligible),
            "n_high": n_high,
            "n_low": n_low,
        }

    # Precompute contributions (num_i, den_i)
    contribs = {}
    for ind in eligible:
        year_doy = by_ind[ind]
        years = sorted(year_doy.keys())
        doys = [year_doy[y] for y in years]
        my = sum(years) / len(years)
        md = sum(doys) / len(doys)
        num_i = 0.0
        den_i = 0.0
        for y, d in zip(years, doys):
            ydev = y - my
            ddev = d - md
            num_i += ydev * ddev
            den_i += ydev * ydev
        contribs[ind] = (num_i, den_i)

    def delta_from_flags(is_high_list):
        num_h = 0.0
        den_h = 0.0
        num_l = 0.0
        den_l = 0.0
        for ind, is_h in zip(eligible, is_high_list):
            ni, di = contribs[ind]
            if is_h:
                num_h += ni
                den_h += di
            else:
                num_l += ni
                den_l += di
        if den_h <= 0.0 or den_l <= 0.0:
            return None
        return (num_h / den_h) - (num_l / den_l)

    obs_flags = [labels[i] for i in eligible]
    obs_delta = delta_from_flags(obs_flags)
    if obs_delta is None:
        return None

    null_deltas = []
    for _ in range(n_perm):
        shuffled = obs_flags[:]
        rng.shuffle(shuffled)
        d = delta_from_flags(shuffled)
        if d is not None:
            null_deltas.append(d)
    if not null_deltas:
        return None

    ge = sum(1 for d in null_deltas if abs(d) >= abs(obs_delta))
    p = (ge + 1) / (len(null_deltas) + 1)   # Phipson-Smyth adjustment
    mean_n = sum(null_deltas) / len(null_deltas)
    sd_n = (statistics.stdev(null_deltas)
            if len(null_deltas) > 1 else 0.0)
    return {
        "obs_delta": obs_delta,
        "mean_null_delta": mean_n,
        "sd_null_delta": sd_n,
        "p_two_sided": p,
        "n_permutations_used": len(null_deltas),
        "n_eligible": len(eligible),
        "n_tenure_matched": n_high,
        "n_short_tenure": n_low,
    }


# ---------------------------- DATA LOADING -------------------------------

def load_data():
    """Download, cache, parse, and consolidate first-leaf observations for
    all configured (species, phenophase) pairs.

    Returns a tuple (by_ind, meta, diagnostics) where:
      * by_ind[individual_id] = {year: first_yes_doy}
      * meta[individual_id]  = {site_id, species_id, state, lat, lon,
                                phenophase_id}
      * diagnostics          = dict with per-species byte counts, SHA-256s,
                                and row counts.
    """
    os.makedirs(CACHE_DIR, exist_ok=True)
    by_ind = {}
    meta = {}
    diagnostics = {
        "species_attempted": len(SPECIES),
        "species_fetched": 0,
        "sha256_by_species": {},
        "bytes_by_species": {},
        "raw_rows_by_species": {},
        "valid_rows_by_species": {},
        "individuals_by_species": {},
    }
    for (sid, ppid, sci, common) in SPECIES:
        key = f"{sid}_{ppid}"
        print(f"  Fetching species_id={sid} ({common}), phenophase={ppid}...",
              flush=True)
        try:
            data, h = fetch_with_cache(sid, ppid)
        except Exception as e:
            print(f"    ERROR: {e}", flush=True)
            diagnostics["sha256_by_species"][key] = None
            diagnostics["bytes_by_species"][key] = 0
            diagnostics["raw_rows_by_species"][key] = 0
            diagnostics["valid_rows_by_species"][key] = 0
            diagnostics["individuals_by_species"][key] = 0
            continue
        diagnostics["sha256_by_species"][key] = h
        diagnostics["bytes_by_species"][key] = len(data)
        diagnostics["species_fetched"] += 1
        rows = parse_npn_csv(data, sid, ppid)
        diagnostics["raw_rows_by_species"][key] = \
            data.count(b"\n")  # approximate row count
        diagnostics["valid_rows_by_species"][key] = len(rows)

        per_ind, per_meta = consolidate_per_individual_year(rows)
        diagnostics["individuals_by_species"][key] = len(per_ind)
        for ind, year_doy in per_ind.items():
            by_ind[ind] = year_doy
            meta[ind] = per_meta[ind]
            meta[ind]["common_name"] = common
            meta[ind]["scientific_name"] = sci

        print(f"    -> {len(rows)} first-yes rows, "
              f"{len(per_ind)} distinct individuals.", flush=True)
    # Fail loud if every single fetch failed (typically USA-NPN is down or
    # the host has no outbound network). We do NOT want to silently write
    # an empty results.json -- that is the "corrupt or empty output" mode
    # we're trying to avoid.
    if diagnostics["species_fetched"] == 0:
        raise RuntimeError(
            "No species CSVs fetched successfully from USA-NPN. "
            "Check network connectivity to services.usanpn.org and "
            "inspect stderr for HTTP retry messages. No partial "
            "results.json will be written."
        )
    if diagnostics["species_fetched"] < MIN_SPECIES_FETCHED:
        sys.stderr.write(
            f"WARNING: only {diagnostics['species_fetched']} / "
            f"{len(SPECIES)} species fetched (minimum for a trustworthy "
            f"cohort analysis = {MIN_SPECIES_FETCHED}). --verify will "
            "fail on this run; re-execute when USA-NPN is reachable.\n"
        )
    return by_ind, meta, diagnostics


# ---------------------------- MAIN ANALYSIS ------------------------------

def _bootstrap_delta(by_ind, ids_high, ids_low,
                     min_years=MIN_YEARS_PER_INDIVIDUAL,
                     n_boot=N_BOOTSTRAP, rng=None):
    """Bootstrap CI on (slope_high - slope_low) by jointly resampling
    individuals within each group with replacement."""
    if rng is None:
        rng = random.Random(SEED + 2)
    contribs_h = []
    contribs_l = []
    # Sort inputs for deterministic iteration / bootstrap indexing.
    ids_high = sorted(ids_high)
    ids_low = sorted(ids_low)
    for ind in ids_high:
        year_doy = by_ind.get(ind, {})
        if len(year_doy) < min_years:
            continue
        years = sorted(year_doy.keys())
        doys = [year_doy[y] for y in years]
        my = sum(years) / len(years)
        md = sum(doys) / len(doys)
        ni = sum((y - my) * (d - md) for y, d in zip(years, doys))
        di = sum((y - my) ** 2 for y in years)
        contribs_h.append((ni, di))
    for ind in ids_low:
        year_doy = by_ind.get(ind, {})
        if len(year_doy) < min_years:
            continue
        years = sorted(year_doy.keys())
        doys = [year_doy[y] for y in years]
        my = sum(years) / len(years)
        md = sum(doys) / len(doys)
        ni = sum((y - my) * (d - md) for y, d in zip(years, doys))
        di = sum((y - my) ** 2 for y in years)
        contribs_l.append((ni, di))
    if len(contribs_h) < 2 or len(contribs_l) < 2:
        return None
    deltas = []
    for _ in range(n_boot):
        nh = 0.0
        dh = 0.0
        for _ in range(len(contribs_h)):
            n_, d_ = contribs_h[rng.randrange(len(contribs_h))]
            nh += n_
            dh += d_
        nl = 0.0
        dl = 0.0
        for _ in range(len(contribs_l)):
            n_, d_ = contribs_l[rng.randrange(len(contribs_l))]
            nl += n_
            dl += d_
        if dh > 0.0 and dl > 0.0:
            deltas.append((nh / dh) - (nl / dl))
    if not deltas:
        return None
    deltas.sort()
    return {
        "point": (sum(ni for ni, _ in contribs_h)
                  / sum(di for _, di in contribs_h)
                  - sum(ni for ni, _ in contribs_l)
                    / sum(di for _, di in contribs_l)),
        "ci95_lo": percentile(deltas, CI_ALPHA_LO),
        "ci95_hi": percentile(deltas, CI_ALPHA_HI),
        "n_boot": len(deltas),
    }


def run_analysis(by_ind, meta):
    """Compute primary results: full-cohort slope, tenure-matched slope,
    difference bootstrap, permutation null, per-species, tenure sensitivity,
    and a robustness check using per-individual Theil-Sen medians.
    """
    # Primary classification at TENURE_THRESHOLD. Sort so the bootstrap
    # and permutation seeds produce reproducible samples across runs
    # regardless of PYTHONHASHSEED.
    classification = classify_individuals(by_ind, TENURE_THRESHOLD)
    eligible = sorted(i for i in classification
                      if len(by_ind.get(i, {})) >= MIN_YEARS_PER_INDIVIDUAL)
    high = [i for i in eligible if classification[i]["tenure_matched"]]
    low = [i for i in eligible if not classification[i]["tenure_matched"]]

    boot_full = bootstrap_slope_over_individuals(
        by_ind, set(eligible), n_boot=N_BOOTSTRAP,
        rng=random.Random(SEED + 10))
    boot_high = bootstrap_slope_over_individuals(
        by_ind, set(high), n_boot=N_BOOTSTRAP,
        rng=random.Random(SEED + 20))
    boot_low = bootstrap_slope_over_individuals(
        by_ind, set(low), n_boot=N_BOOTSTRAP,
        rng=random.Random(SEED + 30))

    # Bootstrap on difference (paired-resample)
    diff_boot = _bootstrap_delta(
        by_ind, high, low, n_boot=N_BOOTSTRAP,
        rng=random.Random(SEED + 40))

    # Permutation null on the same statistic
    perm = permutation_null_tenure_difference(
        by_ind, classification,
        n_perm=N_PERMUTATIONS,
        rng=random.Random(SEED + 50))

    # Per-individual Theil-Sen slope distribution (robustness)
    ts_full = per_individual_slopes(by_ind, include_ids=set(eligible))
    ts_high = per_individual_slopes(by_ind, include_ids=set(high))
    ts_low = per_individual_slopes(by_ind, include_ids=set(low))

    # Negative control / falsification check. Within each individual, shuffle
    # the year -> DOY map so the underlying within-individual variance is
    # preserved but any year trend is destroyed. The resulting pooled slope
    # should sit tightly around zero (|slope| << |observed full-cohort slope|).
    # A non-zero result here would indicate the pooled estimator is leaking
    # between-individual variance (i.e., not actually a within-individual FE).
    shuffled_by_ind = {}
    nc_rng = random.Random(SEED + 99)
    for ind in sorted(eligible):
        year_doy = by_ind[ind]
        years = sorted(year_doy.keys())
        doys = [year_doy[y] for y in years]
        doys_shuffled = doys[:]
        nc_rng.shuffle(doys_shuffled)
        shuffled_by_ind[ind] = {
            y: d for y, d in zip(years, doys_shuffled)
        }
    nc_slope = within_individual_pooled_slope(shuffled_by_ind, set(eligible))
    negative_control = {
        "slope_days_per_year": nc_slope["slope"] if nc_slope else None,
        "n_individuals": nc_slope["n_individuals"] if nc_slope else 0,
        "description": ("Within-individual year<->DOY shuffle. Under the "
                        "null that the pooled estimator is a true "
                        "within-individual FE, this should be ~0."),
    }

    def summarise(arr):
        if not arr:
            return None
        arr_s = sorted(arr)
        n = len(arr_s)
        median = percentile(arr_s, 50)
        mean = sum(arr_s) / n
        # Bootstrap 95% CI on median
        rng = random.Random(SEED + 60)
        medians = []
        for _ in range(N_BOOTSTRAP):
            sample = sorted(arr_s[rng.randrange(n)] for _ in range(n))
            medians.append(sample[n // 2] if n % 2 == 1
                           else 0.5 * (sample[n // 2 - 1] + sample[n // 2]))
        medians.sort()
        return {
            "n": n,
            "median": median,
            "mean": mean,
            "median_ci95_lo": percentile(medians, CI_ALPHA_LO),
            "median_ci95_hi": percentile(medians, CI_ALPHA_HI),
        }

    return {
        "n_individuals_total": len(by_ind),
        "n_individuals_eligible": len(eligible),
        "n_tenure_matched": len(high),
        "n_short_tenure": len(low),
        "slope_full_cohort": boot_full,
        "slope_tenure_matched": boot_high,
        "slope_short_tenure": boot_low,
        "difference_bootstrap": diff_boot,
        "permutation_null": perm,
        "theil_sen_summary": {
            "full": summarise(ts_full),
            "tenure_matched": summarise(ts_high),
            "short_tenure": summarise(ts_low),
        },
        "negative_control": negative_control,
    }


# ---------------------------- SENSITIVITY --------------------------------

def sensitivity_tenure_thresholds(by_ind):
    """Re-run the tenure-vs-short comparison at several tenure thresholds."""
    out = []
    for t in TENURE_SENSITIVITY:
        cls = classify_individuals(by_ind, t)
        eligible = sorted(i for i in cls
                          if len(by_ind.get(i, {})) >= MIN_YEARS_PER_INDIVIDUAL)
        high = [i for i in eligible if cls[i]["tenure_matched"]]
        low = [i for i in eligible if not cls[i]["tenure_matched"]]
        s_high = within_individual_pooled_slope(by_ind, set(high))
        s_low = within_individual_pooled_slope(by_ind, set(low))
        out.append({
            "tenure_threshold": t,
            "n_tenure_matched": len(high),
            "n_short_tenure": len(low),
            "slope_tenure_matched_per_year":
                s_high["slope"] if s_high else None,
            "slope_short_tenure_per_year":
                s_low["slope"] if s_low else None,
            "delta_per_decade":
                ((s_high["slope"] - s_low["slope"]) * 10.0)
                if (s_high and s_low) else None,
        })
    return out


def sensitivity_species_subsets(by_ind, meta):
    """Within-species replication: run the tenure-confound test species by
    species so readers can check whether the result is driven by one taxon."""
    out = []
    for (sid, ppid, sci, common) in SPECIES:
        sub = {i: y for i, y in by_ind.items()
               if meta.get(i, {}).get("species_id") == sid}
        if len(sub) < 20:
            out.append({
                "species_id": sid,
                "common_name": common,
                "n_individuals": len(sub),
                "note": "too few individuals for species-level test",
            })
            continue
        cls = classify_individuals(sub, TENURE_THRESHOLD)
        eligible = sorted(i for i in cls
                          if len(sub.get(i, {})) >= MIN_YEARS_PER_INDIVIDUAL)
        high = [i for i in eligible if cls[i]["tenure_matched"]]
        low = [i for i in eligible if not cls[i]["tenure_matched"]]
        s_high = within_individual_pooled_slope(sub, set(high))
        s_low = within_individual_pooled_slope(sub, set(low))
        s_full = within_individual_pooled_slope(sub, set(eligible))
        out.append({
            "species_id": sid,
            "common_name": common,
            "scientific_name": sci,
            "n_individuals": len(sub),
            "n_eligible": len(eligible),
            "n_tenure_matched": len(high),
            "n_short_tenure": len(low),
            "slope_full_per_year":
                s_full["slope"] if s_full else None,
            "slope_tenure_matched_per_year":
                s_high["slope"] if s_high else None,
            "slope_short_tenure_per_year":
                s_low["slope"] if s_low else None,
        })
    return out


def sensitivity_window(by_ind, start_year, end_year):
    """Re-run primary analysis restricted to observations within a
    sub-window of years. Tenure is recomputed within the sub-window."""
    sub = {}
    for ind, ydoy in by_ind.items():
        pruned = {y: d for y, d in ydoy.items()
                  if start_year <= y <= end_year}
        if pruned:
            sub[ind] = pruned
    cls = classify_individuals(sub, TENURE_THRESHOLD)
    eligible = sorted(i for i in cls
                      if len(sub.get(i, {})) >= MIN_YEARS_PER_INDIVIDUAL)
    high = [i for i in eligible if cls[i]["tenure_matched"]]
    low = [i for i in eligible if not cls[i]["tenure_matched"]]
    s_high = within_individual_pooled_slope(sub, set(high))
    s_low = within_individual_pooled_slope(sub, set(low))
    s_full = within_individual_pooled_slope(sub, set(eligible))
    return {
        "window": [start_year, end_year],
        "n_individuals": len(sub),
        "n_eligible": len(eligible),
        "n_tenure_matched": len(high),
        "n_short_tenure": len(low),
        "slope_full_per_year": s_full["slope"] if s_full else None,
        "slope_tenure_matched_per_year":
            s_high["slope"] if s_high else None,
        "slope_short_tenure_per_year":
            s_low["slope"] if s_low else None,
    }


# ---------------------------- REPORTING ----------------------------------

LIMITATIONS = [
    ("Individual-tenure proxy. The USA-NPN public getSummarizedData endpoint "
     "does not return observer identifiers, so 'tenure' here means continuous "
     "monitoring of a physical plant, not continuous monitoring by the same "
     "human. A plant can change hands across years, so the cohort restriction "
     "is an imperfect (but strictly public) proxy for the people-level "
     "question."),
    ("Unstratified label shuffle. The permutation null shuffles the "
     "tenure-matched flag across all eligible individuals without "
     "stratifying by species, site, or state. If tenure correlates with "
     "species or region, exchangeability is weakened and the reported "
     "p-value is an upper bound on the species-conditional p-value."),
    ("Volunteer-network composition. USA-NPN Nature's Notebook is a "
     "volunteer corps that over-samples certain regions, species, and "
     "landowner types (backyards, arboreta). These results should NOT be "
     "read as a representative US-wide phenology census."),
    ("No climate attribution. The analysis is descriptive -- we do not "
     "estimate the contribution of temperature, precipitation, urban heat "
     "island, or any other climate driver. A non-zero slope is consistent "
     "with many causal stories."),
    ("Sign-instability across tenure thresholds. In the reported USA-NPN "
     "snapshot, the sign of (slope_tenure_matched - slope_short_tenure) "
     "is NOT monotone across tenure_threshold in {3, 5, 7, 10}. We report "
     "the full sweep in sensitivity.tenure_thresholds rather than cherry-"
     "picking a threshold; downstream users should check the sweep before "
     "taking the primary threshold=5 result at face value."),
    ("No cross-(species, individual_id) collision detection. "
     "consolidate_per_individual_year() keys by individual_id alone. This "
     "held on the reported snapshot (per-species individual counts summed "
     "exactly to the pooled total) but if USA-NPN issues a future vintage "
     "where individual_id is no longer globally unique across species, "
     "edit the helper to key on (species_id, individual_id)."),
]


def generate_report(results, diagnostics, sens, by_ind, meta):
    """Emit results.json and report.md."""
    payload = {
        "meta": {
            "analysis_name":
                "USA-NPN Observer-Tenure Confound on Spring Advance",
            "version": "1.0.0",
            "seed": SEED,
            "start_year": START_YEAR,
            "end_year": END_YEAR,
            "tenure_threshold": TENURE_THRESHOLD,
            "min_years_per_individual": MIN_YEARS_PER_INDIVIDUAL,
            "ci_level_pct": CI_LEVEL,
            "significance_threshold": SIGNIFICANCE_THRESHOLD,
            "n_bootstrap": N_BOOTSTRAP,
            "n_permutations": N_PERMUTATIONS,
            "species_queried": [
                {"species_id": sid, "phenophase_id": ppid,
                 "scientific_name": sci, "common_name": common}
                for (sid, ppid, sci, common) in SPECIES
            ],
            "n_individuals_total": len(by_ind),
            "n_observations_total":
                sum(len(v) for v in by_ind.values()),
            "n_sites_total": len({m["site_id"] for m in meta.values()}),
            "n_states_total": len({m["state"] for m in meta.values()
                                   if m.get("state")}),
            "timestamp_utc":
                datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ"),
        },
        "fetch_diagnostics": diagnostics,
        "results": results,
        "sensitivity": sens,
        "limitations": LIMITATIONS,
    }
    # Compute a deterministic hash of the scientific payload (everything
    # except timestamp_utc) so reproducibility can be checked across runs
    # with a single byte-comparison even though the wall-clock field moves.
    hashable = {k: v for k, v in payload.items()}
    hashable_meta = {k: v for k, v in payload["meta"].items()
                     if k != "timestamp_utc"}
    hashable["meta"] = hashable_meta
    canon = json.dumps(hashable, sort_keys=True, default=_json_default)
    payload["meta"]["results_hash"] = hashlib.sha256(
        canon.encode("utf-8")).hexdigest()

    with open(RESULTS_JSON, "w", encoding="utf-8") as f:
        json.dump(payload, f, indent=2, sort_keys=True,
                  default=_json_default)

    # Human-readable report
    lines = []
    lines.append("# USA-NPN Observer-Tenure Confound on Spring Advance")
    lines.append("")
    lines.append(f"- Seed: {SEED}")
    lines.append(f"- Window: {START_YEAR}-{END_YEAR}")
    lines.append(f"- Tenure threshold (primary): >= "
                 f"{TENURE_THRESHOLD} consecutive years")
    lines.append(f"- Individuals total: "
                 f"{payload['meta']['n_individuals_total']}")
    lines.append(f"- Observations total: "
                 f"{payload['meta']['n_observations_total']}")
    lines.append(f"- Sites total: {payload['meta']['n_sites_total']}")
    lines.append(f"- States total: {payload['meta']['n_states_total']}")
    lines.append("")

    def decade(s):
        return None if s is None else s * 10.0

    def fmt_slope_ci(boot):
        if boot is None:
            return "n/a"
        return (f"{decade(boot['point']):+.3f} "
                f"[{decade(boot['ci95_lo']):+.3f}, "
                f"{decade(boot['ci95_hi']):+.3f}]")

    lines.append("## Primary results (days per decade)")
    lines.append("")
    lines.append("| Cohort | n individuals | slope (days/decade) | 95% CI |")
    lines.append("|---|---|---|---|")
    for label, key in [("Full (>= 3 yrs)", "slope_full_cohort"),
                       ("Tenure-matched (>= 5 consec. yrs)",
                        "slope_tenure_matched"),
                       ("Short-tenure (< 5 consec. yrs)",
                        "slope_short_tenure")]:
        boot = results[key]
        if boot is None:
            lines.append(f"| {label} | - | n/a | n/a |")
            continue
        n = boot.get("n_individuals", "?")
        lines.append(
            f"| {label} | {n} | "
            f"{decade(boot['point']):+.3f} | "
            f"[{decade(boot['ci95_lo']):+.3f}, "
            f"{decade(boot['ci95_hi']):+.3f}] |"
        )
    lines.append("")

    lines.append("## Difference (tenure-matched minus short-tenure)")
    lines.append("")
    dd = results["difference_bootstrap"]
    if dd is None:
        lines.append("- Bootstrap difference: n/a (insufficient data).")
    else:
        lines.append(
            f"- Bootstrap delta (days/decade): "
            f"{decade(dd['point']):+.3f} "
            f"[{decade(dd['ci95_lo']):+.3f}, "
            f"{decade(dd['ci95_hi']):+.3f}]"
        )
    perm = results["permutation_null"]
    if perm and "p_two_sided" in perm:
        lines.append(
            f"- Permutation null (two-sided p): "
            f"{perm['p_two_sided']:.4f}  "
            f"(|obs|={abs(decade(perm['obs_delta'])):.4f} days/decade, "
            f"null mean={decade(perm['mean_null_delta']):+.4f}, "
            f"null sd={decade(perm['sd_null_delta']):.4f}, "
            f"n_permutations={perm['n_permutations_used']})"
        )
    lines.append("")

    lines.append("## Per-individual Theil-Sen distribution")
    lines.append("")
    for label, key in [("Full", "full"),
                       ("Tenure-matched", "tenure_matched"),
                       ("Short-tenure", "short_tenure")]:
        s = results["theil_sen_summary"][key]
        if s is None:
            lines.append(f"- {label}: n/a")
            continue
        lines.append(
            f"- {label}: n={s['n']}, "
            f"median slope = {decade(s['median']):+.3f} days/decade "
            f"[{decade(s['median_ci95_lo']):+.3f}, "
            f"{decade(s['median_ci95_hi']):+.3f}], "
            f"mean slope = {decade(s['mean']):+.3f}"
        )
    lines.append("")

    lines.append("## Sensitivity")
    lines.append("")
    lines.append("### Tenure threshold")
    for t in sens["tenure_thresholds"]:
        lines.append(
            f"- threshold={t['tenure_threshold']}: "
            f"tenure={t['n_tenure_matched']}, short={t['n_short_tenure']}, "
            f"slope_tenure={decade(t['slope_tenure_matched_per_year']):+.3f} "
            f"days/dec, slope_short="
            f"{decade(t['slope_short_tenure_per_year']):+.3f} days/dec, "
            f"delta={t['delta_per_decade']:+.3f} days/dec"
            if (t['slope_tenure_matched_per_year'] is not None
                and t['slope_short_tenure_per_year'] is not None) else
            f"- threshold={t['tenure_threshold']}: "
            f"tenure={t['n_tenure_matched']}, short={t['n_short_tenure']}, "
            f"slope n/a"
        )
    lines.append("")
    lines.append("### Species subsets")
    for s in sens["species_subsets"]:
        if "slope_full_per_year" not in s:
            lines.append(
                f"- {s['common_name']} (sid={s['species_id']}): "
                f"n={s['n_individuals']} ({s.get('note','')})"
            )
            continue
        sf = s["slope_full_per_year"]
        st = s["slope_tenure_matched_per_year"]
        ss = s["slope_short_tenure_per_year"]
        lines.append(
            f"- {s['common_name']} (sid={s['species_id']}): "
            f"n={s['n_individuals']} ({s['n_tenure_matched']} tenure-matched), "
            f"slope_full={decade(sf):+.3f}, "
            f"slope_tenure={decade(st):+.3f}, "
            f"slope_short={decade(ss):+.3f}"
        )
    lines.append("")
    lines.append("### Time window")
    for w in sens["windows"]:
        lines.append(
            f"- {w['window'][0]}-{w['window'][1]}: "
            f"n={w['n_individuals']} individuals, "
            f"slope_full={decade(w['slope_full_per_year']):+.3f}, "
            f"slope_tenure={decade(w['slope_tenure_matched_per_year']):+.3f}, "
            f"slope_short={decade(w['slope_short_tenure_per_year']):+.3f}"
        )
    lines.append("")

    lines.append("## Limitations and failure modes")
    lines.append("")
    for lim in LIMITATIONS:
        lines.append(f"- {lim}")
    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 verify():
    """Run machine-checkable assertions against results.json."""
    assert os.path.exists(RESULTS_JSON), (
        f"results.json missing at {RESULTS_JSON}"
    )
    with open(RESULTS_JSON, "r", encoding="utf-8") as f:
        R = json.load(f)

    M = R["meta"]
    res = R["results"]

    # 1. Minimum individual count (sanity floor).
    assert M["n_individuals_total"] >= 500, (
        f"Only {M['n_individuals_total']} individuals; floor = 500."
    )

    # 2. Minimum tenure-matched cohort size.
    assert res["n_tenure_matched"] >= 30, (
        f"Only {res['n_tenure_matched']} tenure-matched individuals; "
        f"floor = 30."
    )

    # 3. Full-cohort slope is present with a bootstrap CI.
    assert res["slope_full_cohort"] is not None, "slope_full_cohort missing."
    for k in ("point", "ci95_lo", "ci95_hi", "n_individuals", "n_obs"):
        assert k in res["slope_full_cohort"], (
            f"slope_full_cohort missing field {k}"
        )
    b = res["slope_full_cohort"]
    assert b["ci95_lo"] <= b["point"] <= b["ci95_hi"], (
        f"Full-cohort CI does not bracket point: {b}"
    )

    # 4. Tenure-matched slope is present with a bootstrap CI.
    assert res["slope_tenure_matched"] is not None, (
        "slope_tenure_matched missing."
    )
    b = res["slope_tenure_matched"]
    assert b["ci95_lo"] <= b["point"] <= b["ci95_hi"], (
        f"Tenure-matched CI does not bracket point: {b}"
    )

    # 5. Plausibility: point slope (days/year) must lie in [-5, 5] for any
    # real phenology dataset. Anything outside that range indicates a bug.
    for key in ("slope_full_cohort", "slope_tenure_matched",
                "slope_short_tenure"):
        if res[key] is None:
            continue
        s = res[key]["point"]
        assert -5.0 <= s <= 5.0, (
            f"{key} point slope {s} is outside plausible [-5, 5] days/year."
        )

    # 6. Permutation null reports an obs delta and two-sided p-value in [0, 1].
    perm = res["permutation_null"]
    assert perm is not None, "permutation_null missing."
    assert "p_two_sided" in perm and 0.0 <= perm["p_two_sided"] <= 1.0, (
        f"Permutation p-value out of range: {perm}"
    )
    assert perm["n_permutations_used"] >= int(0.9 * N_PERMUTATIONS), (
        f"Too few permutations used: {perm['n_permutations_used']}."
    )

    # 7. Difference bootstrap CI bracket consistency.
    dd = res["difference_bootstrap"]
    assert dd is not None, "difference_bootstrap missing."
    assert dd["ci95_lo"] <= dd["ci95_hi"], (
        "difference_bootstrap has inverted CI."
    )

    # 8. Sensitivity sections present with required breadth.
    sens = R["sensitivity"]
    assert "tenure_thresholds" in sens and len(sens["tenure_thresholds"]) >= 3, (
        "Need >= 3 tenure-threshold rows."
    )
    assert "species_subsets" in sens and len(sens["species_subsets"]) >= 4, (
        "Need >= 4 species rows."
    )
    assert "windows" in sens and len(sens["windows"]) >= 2, (
        "Need >= 2 time-window rows."
    )

    # 9. SHA-256 fingerprints recorded for every fetched species CSV.
    diag = R["fetch_diagnostics"]
    for (sid, ppid, _, _) in [
        (m["species_id"], m["phenophase_id"], m["scientific_name"],
         m["common_name"])
        for m in M["species_queried"]
    ]:
        key = f"{sid}_{ppid}"
        h = diag["sha256_by_species"].get(key)
        assert h, f"Missing SHA-256 for species {sid} phenophase {ppid}."
        int(h, 16)  # must be hex

    # 10. Theil-Sen robustness block is populated for full cohort.
    ts = res["theil_sen_summary"]["full"]
    assert ts is not None and ts["n"] > 0, (
        "Theil-Sen full summary missing or empty."
    )
    assert ts["median_ci95_lo"] <= ts["median"] <= ts["median_ci95_hi"], (
        "Theil-Sen median CI does not bracket point."
    )

    # 11. Full-cohort bootstrap CI has realistic width (> ~1% of |point| and
    # not pathologically wide). Guards against degenerate single-point CIs
    # and against no-signal runs.
    b = res["slope_full_cohort"]
    ci_width = b["ci95_hi"] - b["ci95_lo"]
    assert ci_width > 0.02, (
        f"Full-cohort CI width {ci_width:.4f} < 0.02 days/year: "
        "implausibly tight, suggests degenerate bootstrap."
    )
    assert ci_width < 2.0, (
        f"Full-cohort CI width {ci_width:.4f} > 2.0 days/year: "
        "no signal or catastrophic outlier influence."
    )

    # 12. At least four of the six configured species must have produced a
    # non-null SHA-256 fingerprint. Below that floor the cohort analysis
    # is not trustworthy.
    diag = R["fetch_diagnostics"]
    n_ok = sum(1 for v in diag["sha256_by_species"].values() if v)
    assert n_ok >= 4, (
        f"Only {n_ok} species CSVs fetched successfully (need >= 4)."
    )

    # 13. Deterministic results_hash is present and is a valid SHA-256.
    rh = M.get("results_hash")
    assert isinstance(rh, str) and len(rh) == 64, (
        f"meta.results_hash missing or malformed: {rh!r}"
    )
    int(rh, 16)  # must be hex

    # 14. Sensitivity agreement: at least two of the reported tenure
    # thresholds must produce a non-null slope for both sub-cohorts, so a
    # reader can judge robustness across cut-offs without degenerate cells.
    sens_tt = R["sensitivity"]["tenure_thresholds"]
    populated = [row for row in sens_tt
                 if row["slope_tenure_matched_per_year"] is not None
                 and row["slope_short_tenure_per_year"] is not None]
    assert len(populated) >= 2, (
        f"Only {len(populated)} tenure-threshold rows are fully populated."
    )

    # 15. Sensitivity sanity: the full-cohort sign agrees with the
    # tenure-matched sign at the primary threshold (both negative for an
    # advancing phenology, both positive for a delaying one). Guards
    # against a parser or sign bug.
    tm = res["slope_tenure_matched"]["point"]
    fc = res["slope_full_cohort"]["point"]
    assert (tm <= 0.0 and fc <= 0.0) or (tm >= 0.0 and fc >= 0.0), (
        f"Full-cohort ({fc:+.4f}) and tenure-matched ({tm:+.4f}) "
        "slopes have opposite signs -- suspicious."
    )

    # 16. Falsification / null-calibration check: under the label-shuffle
    # null the absolute value of mean_null_delta should be much smaller
    # than the sd_null_delta (the null is centred on zero by construction).
    # A mean/sd ratio larger than 0.5 would indicate the shuffle procedure
    # is not actually exchangeable.
    perm = res["permutation_null"]
    if perm and "sd_null_delta" in perm and perm["sd_null_delta"] > 0:
        ratio = abs(perm["mean_null_delta"]) / perm["sd_null_delta"]
        assert ratio < 0.5, (
            f"Permutation null not centred on zero: "
            f"|mean|/sd = {ratio:.3f} (expected << 0.5)."
        )

    # 17. Observation-count floor: the eligible (>= MIN_YEARS observation
    # years) sample must cover at least 1,000 (individual, year) rows
    # aggregated across cohorts.
    n_obs_full = res["slope_full_cohort"].get("n_obs", 0)
    assert n_obs_full >= 1000, (
        f"Full-cohort eligible observation count {n_obs_full} < 1,000."
    )

    # 18. Per-species rows are present for every configured species.
    sp_rows = R["sensitivity"]["species_subsets"]
    assert len(sp_rows) == len(M["species_queried"]), (
        f"Species subsets rows {len(sp_rows)} != "
        f"configured species {len(M['species_queried'])}."
    )

    # 19. Negative control (within-individual year<->DOY shuffle) must
    # produce a slope whose absolute value is strictly smaller than the
    # full-cohort |slope|. If the shuffle test recovers a similar-magnitude
    # slope, the pooled estimator is NOT actually removing between-
    # individual composition drift and the headline finding is unsafe.
    nc = res.get("negative_control")
    assert nc is not None and nc.get("slope_days_per_year") is not None, (
        "negative_control block missing from results.json."
    )
    fc_abs = abs(res["slope_full_cohort"]["point"])
    nc_abs = abs(nc["slope_days_per_year"])
    assert nc_abs < fc_abs, (
        f"Negative-control |slope| ({nc_abs:.4f}) is not strictly smaller "
        f"than the observed full-cohort |slope| ({fc_abs:.4f}); the "
        "within-individual FE claim is unsafe."
    )

    # 20. Limitations block is populated (at least 4 caveats recorded).
    lims = R.get("limitations", [])
    assert isinstance(lims, list) and len(lims) >= 4, (
        f"limitations block too thin: {len(lims)} entries (need >= 4)."
    )

    # 21. Configuration echo: the meta block must contain CI level and
    # significance threshold so a reader can reproduce the inference
    # rules without re-parsing the script.
    assert "ci_level_pct" in M, "meta.ci_level_pct missing."
    assert M["ci_level_pct"] == CI_LEVEL, (
        f"CI level recorded {M['ci_level_pct']} != configured {CI_LEVEL}."
    )
    assert "significance_threshold" in M, (
        "meta.significance_threshold missing."
    )

    # 22. Observation-count floor meets MIN_N_OBS_FLOOR.
    n_obs_full = res["slope_full_cohort"].get("n_obs", 0)
    assert n_obs_full >= MIN_N_OBS_FLOOR, (
        f"Full-cohort eligible observation count {n_obs_full} "
        f"< MIN_N_OBS_FLOOR={MIN_N_OBS_FLOOR}."
    )

    # 23. Window sensitivity must yield at least one window whose full
    # slope is in the physically plausible range. Guards against a data
    # pipeline that silently drops all observations in every sub-window.
    win_rows = R["sensitivity"]["windows"]
    in_range = [
        w for w in win_rows
        if w.get("slope_full_per_year") is not None
        and SLOPE_PLAUSIBLE_MIN <= w["slope_full_per_year"]
        <= SLOPE_PLAUSIBLE_MAX
    ]
    assert len(in_range) >= 1, (
        f"No sub-window produced a plausible full slope "
        f"(range {SLOPE_PLAUSIBLE_MIN}..{SLOPE_PLAUSIBLE_MAX})."
    )

    print("VERIFY: all assertions passed.")
    print("ALL CHECKS PASSED")
    return True


# ---------------------------- MAIN ---------------------------------------

def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--verify", action="store_true",
                    help="Run verification assertions against an existing "
                         "results.json.")
    args = ap.parse_args()

    random.seed(SEED)
    if args.verify:
        verify()
        return

    print("[1/5] Downloading and parsing USA-NPN first-leaf records.",
          flush=True)
    by_ind, meta, diag = load_data()
    print(f"  Retained {len(by_ind)} individuals "
          f"({sum(len(v) for v in by_ind.values())} "
          f"(individual, year) first-yes observations) "
          f"across {diag['species_fetched']} / {len(SPECIES)} species.",
          flush=True)

    print("[2/5] Computing primary within-individual trend comparisons.",
          flush=True)
    results = run_analysis(by_ind, meta)

    print("[3/5] Running sensitivity analyses.", flush=True)
    sens = {
        "tenure_thresholds": sensitivity_tenure_thresholds(by_ind),
        "species_subsets": sensitivity_species_subsets(by_ind, meta),
        "windows": [
            sensitivity_window(by_ind, START_YEAR, START_YEAR + 9),
            sensitivity_window(by_ind, START_YEAR + 5, END_YEAR),
            sensitivity_window(by_ind, START_YEAR, END_YEAR),
        ],
    }

    print("[4/5] Writing results.json and report.md.", flush=True)
    generate_report(results, diag, sens, by_ind, meta)

    print("[5/5] Summary:", flush=True)

    def decade(s):
        return None if s is None else s * 10.0

    bf = results["slope_full_cohort"]
    bh = results["slope_tenure_matched"]
    bl = results["slope_short_tenure"]
    if bf:
        print(f"  Full cohort:       slope = "
              f"{decade(bf['point']):+.3f} days/decade "
              f"[{decade(bf['ci95_lo']):+.3f}, "
              f"{decade(bf['ci95_hi']):+.3f}], "
              f"n={bf['n_individuals']}", flush=True)
    if bh:
        print(f"  Tenure-matched:    slope = "
              f"{decade(bh['point']):+.3f} days/decade "
              f"[{decade(bh['ci95_lo']):+.3f}, "
              f"{decade(bh['ci95_hi']):+.3f}], "
              f"n={bh['n_individuals']}", flush=True)
    if bl:
        print(f"  Short-tenure:      slope = "
              f"{decade(bl['point']):+.3f} days/decade "
              f"[{decade(bl['ci95_lo']):+.3f}, "
              f"{decade(bl['ci95_hi']):+.3f}], "
              f"n={bl['n_individuals']}", flush=True)
    dd = results["difference_bootstrap"]
    if dd:
        print(f"  Difference (T - S): "
              f"{decade(dd['point']):+.3f} days/decade "
              f"[{decade(dd['ci95_lo']):+.3f}, "
              f"{decade(dd['ci95_hi']):+.3f}]", flush=True)
    perm = results["permutation_null"]
    if perm and "p_two_sided" in perm:
        print(f"  Permutation p (two-sided): {perm['p_two_sided']:.4f} "
              f"({perm['n_permutations_used']} perms)", flush=True)

    print("ANALYSIS COMPLETE", flush=True)


if __name__ == "__main__":
    try:
        main()
    except AssertionError as e:
        # --verify failures land here. Re-raise the traceback to stderr
        # for the caller but exit with a distinguishable non-zero code.
        sys.stderr.write(f"VERIFICATION FAILED: {e}\n")
        sys.exit(2)
    except KeyboardInterrupt:
        sys.stderr.write("INTERRUPTED BY USER\n")
        sys.exit(130)
    except Exception as e:
        # Anything else -- network error, parse error, disk-full, API
        # outage -- is funnelled through here so the caller gets a
        # single clear message and a non-zero exit code instead of a
        # raw traceback interleaved with partial stdout.
        sys.stderr.write(
            f"FATAL: analysis aborted: {type(e).__name__}: {e}\n"
            "       Inspect cache/ and hash_index.json for partial "
            "state; USA-NPN outages are the most common root cause.\n"
        )
        sys.exit(1)
SCRIPT_EOF
```

**Expected output:** No output. The heredoc writes the Python script to `/tmp/claw4s_auto_usa-npn-observer-tenure-confound-on-spring-advance/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_usa-npn-observer-tenure-confound-on-spring-advance && python3 analysis.py
```

**Expected output (cold run, ~1-5 minutes):**

```
[1/5] Downloading and parsing USA-NPN first-leaf records.
  Fetching species_id=36 (common lilac), phenophase=373...
    -> N first-yes rows, K distinct individuals.
  Fetching species_id=3 (red maple), phenophase=371...
    -> N first-yes rows, K distinct individuals.
  ... (six species)
  Retained R individuals (O (individual, year) first-yes observations) across 6 / 6 species.
[2/5] Computing primary within-individual trend comparisons.
[3/5] Running sensitivity analyses.
[4/5] Writing results.json and report.md.
[5/5] Summary:
  Full cohort:       slope = -X.XXX days/decade [lo, hi], n=N
  Tenure-matched:    slope = -Y.YYY days/decade [lo, hi], n=N
  Short-tenure:      slope = -Z.ZZZ days/decade [lo, hi], n=N
  Difference (T - S): +/-W.WWW days/decade [lo, hi]
  Permutation p (two-sided): 0.PPPP (1000 perms)
ANALYSIS COMPLETE
```

**Expected outputs on disk:**

- `results.json` (structured results, ~10 KB)
- `report.md` (human-readable summary)
- `cache/species_<id>_phenophase_<id>.csv` (one per fetched species)
- `cache/hash_index.json` (SHA-256 of every cached CSV)

**Success:** stdout ends with `ANALYSIS COMPLETE`. Both `results.json` and `report.md` exist. `cache/hash_index.json` has six entries.

**Failure:** If any species CSV fetch fails after 4 retries, the species is skipped and recorded with a `null` SHA-256 in `fetch_diagnostics`; the run continues. If fewer than 500 total individuals remain after filtering, verification will fail.

---

## Step 4: Verify results

```bash
cd /tmp/claw4s_auto_usa-npn-observer-tenure-confound-on-spring-advance && python3 analysis.py --verify
```

**Expected output:**

```
VERIFY: all assertions passed.
ALL CHECKS PASSED
```

**Success criteria (all twenty-three assertions must pass):**

1. `results.json` exists and parses as JSON.
2. Total individuals retained >= `MIN_INDIVIDUALS_FLOOR` (=500).
3. Tenure-matched individuals >= `MIN_TENURE_MATCHED_FLOOR` (=30).
4. `slope_full_cohort` has point, CI, and sample-size fields, and the CI brackets the point estimate.
5. `slope_tenure_matched` has point, CI, and sample-size fields, and the CI brackets the point estimate.
6. All three within-cohort point slopes lie in the physically plausible `[SLOPE_PLAUSIBLE_MIN, SLOPE_PLAUSIBLE_MAX]` = `[-5, +5]` days-per-year range (i.e. `[-50, +50]` days per decade).
7. Permutation null reports a two-sided p-value in `[0, 1]` and used at least 90% of the requested 1,000 permutations.
8. Difference bootstrap CI is non-inverted.
9. Sensitivity block includes >= 3 tenure-threshold rows, >= 4 species rows, and >= 2 time-window rows.
10. Every species whose CSV was fetched has a recorded SHA-256 fingerprint (valid hex).
11. Full-cohort 95% bootstrap CI width lies in `(CI_WIDTH_MIN, CI_WIDTH_MAX)` = `(0.02, 2.0)` days/year -- neither degenerate nor signal-free.
12. At least `MIN_SPECIES_FETCHED` (=4) of the six configured species produced a non-null SHA-256 fingerprint (the cohort analysis is not trustworthy below four species).
13. `meta.results_hash` is a valid 64-character hex SHA-256 of the scientific payload (used for byte-level reproducibility checks across runs).
14. At least two tenure-threshold sensitivity rows produce non-null slopes for *both* sub-cohorts (robustness across cut-offs).
15. Full-cohort and tenure-matched point slopes share the same sign (falsification guard against sign / parser bugs).
16. Permutation null is centred on zero: `|mean_null_delta| / sd_null_delta < 0.5` (falsification guard -- an exchangeable shuffle should be centred by construction).
17. Full-cohort eligible sample covers at least `MIN_N_OBS_FLOOR` (=1,000) (individual, year) observations.
18. The per-species sensitivity block reports one row for every configured species (no silent species drop).
19. Negative control (within-individual year<->DOY shuffle) produces an absolute slope strictly smaller than the observed full-cohort absolute slope (the within-individual FE claim is safe).
20. `results.json` `limitations` block contains at least four caveats (so downstream readers cannot miss them).
21. `meta.ci_level_pct` and `meta.significance_threshold` are echoed into `results.json` so the inference rules are auditable without re-parsing the script.
22. Full-cohort eligible observation count >= `MIN_N_OBS_FLOOR` (redundant-by-design with #17 but keyed off the named constant to surface changes when the constant is tuned).
23. At least one time-window sensitivity row yields a full slope within the physically plausible range (guards against a pipeline that silently drops all observations in every sub-window).

**Failure conditions:**

- *"Only N individuals; floor = 500"*: multiple species fetches failed. Inspect `fetch_diagnostics.sha256_by_species` for `null` entries and re-run Step 3 (USA-NPN may be transiently unreachable).
- *"slope ... is outside plausible [-5, 5] days/year"*: a parser bug (for example, sentinel -9999 values leaking through) or an off-by-one in the per-individual consolidation. Re-inspect `parse_npn_csv()` and `consolidate_per_individual_year()`.
- *"Full-cohort CI does not bracket point"*: numerical instability in the bootstrap resampler; verify `random.seed(SEED)` is being called.

---

## Notes on reproducibility

- Every random operation uses `random.Random(SEED)` or a derived sub-seed (`SEED + k`). Re-running the script on an identical cache yields byte-identical `results.json` (modulo the `timestamp_utc` field).
- Every USA-NPN CSV is SHA-256 fingerprinted on first download. Subsequent re-runs verify the fingerprint before parsing; if the cached bytes drift from the stored hash, the file is re-fetched and the hash index updated.
- No third-party Python packages are required (stdlib only): `urllib`, `csv`, `json`, `hashlib`, `math`, `random`, `statistics`, `datetime`, `argparse`, `os`, `sys`, `time`, `collections`, `io`.
- Bootstrap and permutation replicate counts are pinned at 1,000 each via the `N_BOOTSTRAP` and `N_PERMUTATIONS` constants. Raising these improves CI / p-value resolution at the cost of wall-time.
- Tenure is defined as the longest run of consecutive years with a valid first-yes record per individual. The primary threshold (5 years) follows the prompt; a full sweep over 3, 5, 7, and 10 years is reported in the sensitivity block so readers can check robustness.
- Known limitation -- "tenure" here is *individual*-tenure (continuous observation of a physical plant), not observer-tenure, because USA-NPN's public `getSummarizedData` endpoint does not return observer identifiers. A single plant can be observed by multiple people over years, so the cohort restriction is an imperfect (but strictly public) proxy for the people-level question.
- Empirical note on `Individual_ID` uniqueness. The `consolidate_per_individual_year()` helper keys records by `Individual_ID` alone. In the reported run (2026-04-18 USA-NPN fetch) the per-species individual counts (1126 + 2090 + 325 + 154 + 126 + 774) sum to exactly the reported `n_individuals_total = 4595`, so no cross-species key collisions occurred on that vintage. If you re-run on a future USA-NPN snapshot where `Individual_ID` is no longer globally unique across species, edit the helper to key on `(species_id, individual_id)`; the verify harness will surface the change as a drop in total individuals.
- Known limitation -- the label-permutation null shuffles the tenure-matched flag across all eligible individuals without stratifying by species, site, or state. If tenure correlates with species (for example, if lilac observers tend to have longer individual tenure than flowering-dogwood observers), the exchangeability assumption is weakened. A species-stratified permutation is a drop-in replacement in `permutation_null_tenure_difference()` if a downstream user wants it.
- Known limitation -- the direction of the tenure-vs-short-tenure slope difference is NOT stable across the sensitivity sweep of tenure thresholds (at threshold=3 and 5 the tenure-matched cohort advances faster; at threshold=7 and 10 the short-tenure cohort advances faster). This is reported honestly in `results.json -> sensitivity.tenure_thresholds` and should be examined by any downstream user before taking the primary threshold=5 result at face value; the most plausible explanation is shrinking tenure-matched sample size inflating sampling variance, not a meaningful sign change.