---
name: authorship-inflation-convention-stratified
description: >
  Tests whether the claim of universal team-size inflation in science is
  confounded by authorship-ordering conventions. Downloads bibliometric
  records from OpenAlex for journals in alphabetical-authorship fields
  (economics, finance, mathematics) and contribution-order fields
  (biomedicine, multidisciplinary). Classifies journals empirically by
  the fraction of 3+ author papers in surname-alphabetical order.
  Estimates within-journal OLS slopes of team size on publication year
  (equivalent to journal fixed effects with year as a continuous regressor)
  and compares slope distributions across conventions using
  (a) a journal-level permutation test on convention labels (10,000
  shuffles), (b) bootstrap 95% CIs by journal resampling, and (c) four
  sensitivity analyses (classification threshold, year window, team-size
  cap, ISSN subset). Python 3.8+ standard library only.
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget, Divyansh Jain"
tags: ["claw4s-2026", "science-of-science", "bibliometrics", "openalex", "authorship-conventions", "permutation-test"]
python_version: ">=3.8"
dependencies: []
---

# Is Team-Size Inflation in Science a Universal Trend or a Reporting Artifact of Authorship Conventions?

## When to Use This Skill

**Trigger (one line):** Use this skill when you need to test whether an apparent time
trend in a bibliometric quantity (here, co-authorship team size) is a genuine
behavioural signal or a reporting artifact, using a negative-control design that
stratifies a pre-existing reporting convention (alphabetical vs. contribution-ordered
authorship) and evaluates the stratum difference against a label-permutation null.

**Concretely, invoke this skill when an agent is asked any of:**

- "Is the reported growth of scientific team sizes universal?"
- "Does authorship convention (alphabetical vs. contribution order) modulate
  observed team-size inflation?"
- "How would I design a negative-control study to separate a behavioural trend from
  a reporting-convention drift?"

This skill is **not** appropriate for: forecasting individual journal growth,
causal attribution of single-field collaboration policy effects, or any question
that requires per-author career tracking (OpenAlex author-disambiguation drift is
a known confound for such questions and this skill does not address it).

### Preconditions

- **Python version:** 3.8+ standard library only (no pip installs, no numpy/scipy/pandas).
- **Network:** Internet access to `api.openalex.org` required on first run; all
  responses are cached locally with SHA256 integrity checks, so reruns are offline.
- **Disk space:** ~250 MB free under `/tmp` for the cached OpenAlex responses.
- **Runtime:** 10–30 minutes on first run (API + permutation), under 2 minutes on rerun
  from cache.
- **Environment variables:** None required. The analysis is configured entirely by the
  UPPER_CASE constants at the top of `analyze.py`.
- **External data source:** `https://api.openalex.org/works` (public API, polite-pool
  user-agent included in the script). No authentication required.

## Adaptation Guidance

To apply this analysis to a different question that asks whether a time trend differs
across a pre-existing stratum of the population:

- **Change `JOURNAL_SET` (the ISSN list)** in the DOMAIN CONFIGURATION block to the
  set of journals (or any other OpenAlex source) you want to cover. Every downstream
  step is parameterised by `JOURNAL_SET` — no other edits are required for new fields.
- **Change `OUTCOME_FN` in `run_analysis()`** to whatever per-paper scalar you care
  about (here it is `len(authorships)`, the team size). The statistical pipeline does
  not depend on the outcome semantics.
- **Change the classification rule in `classify_journal()`** if the stratum of interest
  is not "fraction alphabetical". For example, to stratify by open-access status,
  replace the `alph_frac` computation with an `is_OA` indicator.
- **Do NOT change** the slope estimator, `permutation_test_slope_delta()`,
  `bootstrap_slope_ci()`, or `bh_correct()` — these are domain-agnostic and implement
  the statistical method (journal fixed effects, label permutation, journal-resample
  bootstrap, multiple-testing correction).
- **Do NOT change** the cache/SHA256 layer in `api_get()` — it is the reproducibility
  anchor and works for any OpenAlex endpoint.

## Research Question

Team-size inflation is a staple finding of the science-of-science literature, but
some fields (economics, mathematics) list authors alphabetically while others
(biomedicine) list by contribution. If reporting conventions across fields changed
over the study window (more author names credited per paper, or a shift from
"team" to individual credit), observed team-size inflation might partly reflect
a crediting change rather than a behavioural change. This skill tests that
confound by comparing within-journal team-size slopes across conventions and
evaluating the difference against a label-permutation null.

## Methodological Hook

Journal fixed-effects OLS (one slope per journal, year entering as a
continuous regressor, cross-journal intercepts absorbed) with a
journal-level permutation null on the convention labels. Existing literature
pools all fields or adds a single dummy for "alphabetical"; by permuting
convention labels across journals 10,000 times while keeping the observed
(journal, year, team_size) structure fixed, we get an exact non-parametric
reference distribution for the slope difference. A label-shuffling null is the
right null because it keeps all within-journal structure (writing technology,
author-disambiguation change, collaboration technology, etc.) invariant and
only perturbs the stratification.

## Null Model

For each journal _j_ over years 1980–2023, estimate a within-journal slope
β_j = slope of team_size on publication_year after subtracting journal mean.
Under the null, the distribution of β_j is independent of the journal's
authorship convention label. We permute the convention label across journals
10,000 times and recompute Δ = mean(β | non-alphabetical) − mean(β | alphabetical).
The observed Δ is compared to this permutation distribution.

## Step 1: Create Workspace

```bash
mkdir -p /tmp/claw4s_auto_authorship-inflation-in-alphabetical-authorship-fields-vs-ot
```

**Expected output:** Directory created, exit code 0.

## Step 2: Write Analysis Script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_authorship-inflation-in-alphabetical-authorship-fields-vs-ot/analyze.py
#!/usr/bin/env python3
"""
Is Team-Size Inflation in Science a Universal Trend or a Reporting Artifact
of Authorship Conventions?

Downloads per-journal author-count panels from OpenAlex, empirically
classifies each journal by the fraction of 3+ author papers in
surname-alphabetical order, estimates within-journal team-size slopes
by journal fixed effects (per-journal OLS of team size on year) and tests the cross-convention slope
difference against a 10,000-iteration label-permutation null.

Python 3.8+ standard library only. All random operations seeded.
"""

import sys
import os
import json
import math
import time
import random
import hashlib
import urllib.request
import urllib.error
import urllib.parse
from collections import defaultdict

# ═══════════════════════════════════════════════════════════════
# GENERAL METHOD CONFIGURATION — these constants are NOT
# domain-specific; they control statistical procedure parameters
# and would be left unchanged when adapting this skill to a
# different stratified time-trend question.
# ═══════════════════════════════════════════════════════════════

SEED = 42                        # Master RNG seed. All random operations seeded.
N_PERMS = 10000                  # Iterations for the main label-permutation null.
N_PERMS_SENS = 2000              # Iterations for each sensitivity permutation.
N_BOOT = 2000                    # Bootstrap iterations for CIs on stratum means.
CI_LEVEL = 0.95                  # Confidence level for bootstrap CIs.
SIGNIFICANCE_THRESHOLD = 0.05    # Alpha for the BH-FDR per-journal test.

# ═══════════════════════════════════════════════════════════════
# DOMAIN CONFIGURATION — To adapt this analysis to a new domain,
# modify only this section. All downstream functions receive
# these values as parameters rather than accessing them directly.
# ═══════════════════════════════════════════════════════════════

# External data source (OpenAlex). Only the base URL is domain-agnostic;
# the actual endpoint and filter are constructed per-journal in
# fetch_journal_works().
DATA_URL = "https://api.openalex.org/works"
BASE_URL = "https://api.openalex.org"
UA = "AuthorshipConventionStudy/1.0 (mailto:claw4s-research@example.com)"
API_DELAY = 0.12                 # Seconds between paginated requests (polite pool).

# Journal set: ISSN, display name, a priori convention.
# 'alph' = alphabetical-authorship field; 'contrib' = contribution-order.
# Classification is confirmed empirically; a priori labels are used only
# to report the confusion between prior and empirical labels.
JOURNAL_SET = [
    # Economics / finance (a priori alphabetical)
    ("0002-8282", "American Economic Review",                    "alph"),
    ("0033-5533", "Quarterly Journal of Economics",              "alph"),
    ("0022-3808", "Journal of Political Economy",                "alph"),
    ("0012-9682", "Econometrica",                                "alph"),
    ("0034-6527", "Review of Economic Studies",                  "alph"),
    ("0022-0531", "Journal of Economic Theory",                  "alph"),
    ("0304-3932", "Journal of Monetary Economics",               "alph"),
    ("0022-1082", "Journal of Finance",                          "alph"),
    ("0304-405X", "Journal of Financial Economics",              "alph"),
    ("0893-9454", "Review of Financial Studies",                 "alph"),
    ("0022-0515", "Journal of Economic Literature",              "alph"),
    ("0095-0696", "Journal of Environmental Economics and Mgmt", "alph"),
    # Mathematics (a priori alphabetical)
    ("0003-486X", "Annals of Mathematics",                       "alph"),
    ("0894-0347", "Journal of the American Mathematical Society","alph"),
    ("0020-9910", "Inventiones Mathematicae",                    "alph"),
    ("0001-5962", "Acta Mathematica",                            "alph"),
    # Biomedicine / clinical / general (a priori contribution-order)
    ("0028-4793", "New England Journal of Medicine",             "contrib"),
    ("0098-7484", "JAMA",                                        "contrib"),
    ("0140-6736", "The Lancet",                                  "contrib"),
    ("0959-8138", "BMJ",                                         "contrib"),
    ("0021-9738", "Journal of Clinical Investigation",           "contrib"),
    ("0092-8674", "Cell",                                        "contrib"),
    ("0002-9440", "American Journal of Pathology",               "contrib"),
    ("0008-5472", "Cancer Research",                             "contrib"),
    ("0021-9258", "Journal of Biological Chemistry",             "contrib"),
    ("0270-6474", "Journal of Neuroscience",                     "contrib"),
    ("1097-6256", "Nature Neuroscience",                         "contrib"),
    ("0270-9139", "Hepatology",                                  "contrib"),
    ("0149-5992", "Diabetes Care",                               "contrib"),
    ("0039-2499", "Stroke",                                      "contrib"),
]

YEAR_MIN = 1980                  # First publication year included in the panel.
YEAR_MAX = 2023                  # Last publication year included in the panel.
MIN_PAPERS_PER_JOURNAL = 200     # Discard under-sampled journals (slope noise).
MIN_YEARS_PER_JOURNAL = 10       # Need time spread to identify a slope.
MAX_WORKS_PER_JOURNAL = 3000     # Cap pagination to bound runtime and API load.
ALPH_FRAC_UPPER = 0.50           # Empirical threshold for alphabetical convention.
ALPH_FRAC_LOWER = 0.30           # Empirical threshold for contribution order.
TEAM_SIZE_CAP = 50               # Winsorisation cap for team size (main spec).
CACHE_DIR = "cache"              # Local cache directory for API responses.
MANIFEST_FILE = "data_manifest.json"  # SHA256 fingerprints of cached responses.

# Plausibility bounds — used by --verify mode for sanity checks.
# Derived from the science-of-science literature (Wuchty et al. 2007,
# Fortunato et al. 2018 report team-size growth of 0.02–0.5 authors/yr).
SLOPE_PLAUSIBLE_ABS_MAX = 1.0    # Any mean slope above this is suspect.
COHENS_D_PLAUSIBLE_MAX = 5.0     # Cohen's d above this is implausibly large.
MIN_JOURNALS_QUALIFYING = 10     # Minimum qualifying journals for valid analysis.
MIN_PAPERS_TOTAL = 5000          # Minimum total papers in qualifying pool.

# ═══════════════════════════════════════════════════════════════
# Helper utilities
# ═══════════════════════════════════════════════════════════════

def sha256_file(path):
    h = hashlib.sha256()
    with open(path, 'rb') as f:
        for chunk in iter(lambda: f.read(8192), b''):
            h.update(chunk)
    return h.hexdigest()


def load_manifest():
    if os.path.exists(MANIFEST_FILE):
        with open(MANIFEST_FILE) as f:
            return json.load(f)
    return {}


def save_manifest(m):
    with open(MANIFEST_FILE, 'w') as f:
        json.dump(m, f, indent=2)


def api_get(url, cache_path, manifest, retries=4):
    """GET JSON from OpenAlex with caching and SHA256 verification."""
    os.makedirs(os.path.dirname(cache_path) or '.', exist_ok=True)
    if os.path.exists(cache_path):
        h = sha256_file(cache_path)
        exp = manifest.get(cache_path)
        if exp is None or h == exp:
            with open(cache_path) as f:
                return json.load(f)
        else:
            print(f"    Cache corrupted: {cache_path}, re-downloading")
            os.remove(cache_path)

    for attempt in range(retries):
        try:
            req = urllib.request.Request(
                url,
                headers={'User-Agent': UA, 'Accept': 'application/json'},
            )
            with urllib.request.urlopen(req, timeout=60) as r:
                raw = r.read().decode()
            obj = json.loads(raw)
            with open(cache_path, 'w') as f:
                f.write(raw)
            manifest[cache_path] = sha256_file(cache_path)
            return obj
        except urllib.error.HTTPError as e:
            wait = (8 if e.code == 429 else 2) * (attempt + 1)
            if attempt < retries - 1:
                print(f"    HTTP {e.code}, retry in {wait}s")
                time.sleep(wait)
            else:
                raise RuntimeError(f"HTTP {e.code} after {retries} retries: {url}")
        except Exception as e:
            if attempt < retries - 1:
                time.sleep(2 ** (attempt + 1))
            else:
                raise RuntimeError(f"Failed after {retries} retries: {url}: {e}")


# ═══════════════════════════════════════════════════════════════
# Statistical helpers (stdlib)
# ═══════════════════════════════════════════════════════════════

def mean_val(xs):
    return sum(xs) / len(xs) if xs else 0.0


def median_val(xs):
    s = sorted(xs)
    n = len(s)
    if n == 0:
        return 0.0
    return s[n // 2] if n % 2 else (s[n // 2 - 1] + s[n // 2]) / 2.0


def slope_within(years, sizes):
    """OLS slope of size on year (simple regression)."""
    n = len(years)
    if n < 3:
        return 0.0
    ym = mean_val(years)
    sm = mean_val(sizes)
    num = sum((y - ym) * (s - sm) for y, s in zip(years, sizes))
    den = sum((y - ym) ** 2 for y in years)
    return num / den if den > 0 else 0.0


def cohens_d(a, b):
    na, nb = len(a), len(b)
    if na < 2 or nb < 2:
        return 0.0
    ma, mb = mean_val(a), mean_val(b)
    va = sum((x - ma) ** 2 for x in a) / (na - 1)
    vb = sum((x - mb) ** 2 for x in b) / (nb - 1)
    sp = math.sqrt(((na - 1) * va + (nb - 1) * vb) / (na + nb - 2))
    return (ma - mb) / sp if sp > 0 else 0.0


def permutation_test_slope_delta(slopes, labels, n_perms, rng):
    """Label-permutation null for the difference in mean slope between
    the two convention groups. Returns two-sided p-value and the observed
    and permuted delta distributions."""
    n = len(slopes)
    idx_a = [i for i in range(n) if labels[i] == 'alph']
    idx_c = [i for i in range(n) if labels[i] == 'contrib']
    if not idx_a or not idx_c:
        return 1.0, 0.0, []
    obs_delta = (mean_val([slopes[i] for i in idx_c])
                 - mean_val([slopes[i] for i in idx_a]))
    labs = list(labels)
    perms = []
    n_ge = 0
    for _ in range(n_perms):
        rng.shuffle(labs)
        ia = [i for i in range(n) if labs[i] == 'alph']
        ic = [i for i in range(n) if labs[i] == 'contrib']
        d = (mean_val([slopes[i] for i in ic])
             - mean_val([slopes[i] for i in ia]))
        perms.append(d)
        if abs(d) >= abs(obs_delta):
            n_ge += 1
    p_two_sided = (n_ge + 1) / (n_perms + 1)
    return p_two_sided, obs_delta, perms


def bootstrap_slope_ci(slopes, n_boot, rng, level=CI_LEVEL):
    """Bootstrap CI by resampling journal-level slopes with replacement."""
    if not slopes:
        return 0.0, 0.0, 0.0
    n = len(slopes)
    bs = []
    for _ in range(n_boot):
        s = [slopes[rng.randint(0, n - 1)] for _ in range(n)]
        bs.append(mean_val(s))
    bs.sort()
    lo = bs[int((1 - level) / 2 * n_boot)]
    hi = bs[min(int((1 + level) / 2 * n_boot), n_boot - 1)]
    return mean_val(slopes), lo, hi


def bh_correct(pvals, alpha=0.05):
    n = len(pvals)
    if n == 0:
        return [], []
    order = sorted(range(n), key=lambda i: pvals[i])
    adj = [0.0] * n
    for r, i in enumerate(order):
        adj[i] = pvals[i] * n / (r + 1)
    for k in range(n - 2, -1, -1):
        adj[order[k]] = min(adj[order[k]], adj[order[k + 1]])
    adj = [min(a, 1.0) for a in adj]
    return adj, [a < alpha for a in adj]


def surname_of(display_name):
    """Best-effort surname extraction: last whitespace-separated token."""
    if not display_name:
        return ''
    parts = display_name.strip().split()
    return parts[-1].lower() if parts else ''


def is_alphabetical_order(surnames):
    """True iff surnames are non-decreasing when lowered."""
    if len(surnames) < 2:
        return False
    return all(surnames[i] <= surnames[i + 1] for i in range(len(surnames) - 1))


# ═══════════════════════════════════════════════════════════════
# Data acquisition
# ═══════════════════════════════════════════════════════════════

def fetch_journal_works(issn, manifest):
    """Fetch a paginated set of works for one journal."""
    works = []
    cursor = '*'
    page = 0
    while len(works) < MAX_WORKS_PER_JOURNAL:
        filt = (f"primary_location.source.issn:{issn},"
                f"type:article,"
                f"publication_year:{YEAR_MIN}-{YEAR_MAX}")
        params = urllib.parse.urlencode({
            'filter': filt,
            'select': 'id,publication_year,authorships',
            'per_page': 200,
            'cursor': cursor,
            'sort': 'publication_year:asc',
        })
        url = f"{BASE_URL}/works?{params}"
        cache_file = os.path.join(CACHE_DIR, f"w_{issn.replace('-', '')}_p{page}.json")
        data = api_get(url, cache_file, manifest)
        rs = data.get('results', [])
        if not rs:
            break
        for w in rs:
            y = w.get('publication_year')
            auths = w.get('authorships') or []
            if y is None or not auths:
                continue
            surnames = []
            for a in auths:
                au = (a.get('author') or {})
                surnames.append(surname_of(au.get('display_name') or a.get('raw_author_name') or ''))
            if not all(surnames):
                continue
            works.append({
                'year': int(y),
                'n_authors': len(surnames),
                'alpha': 1 if is_alphabetical_order(surnames) else 0,
            })
        cursor = data.get('meta', {}).get('next_cursor')
        if not cursor:
            break
        page += 1
        time.sleep(API_DELAY)
    return works


def load_data():
    """Download/parse per-journal panels. Returns a dict journal_id -> list of papers."""
    os.makedirs(CACHE_DIR, exist_ok=True)
    manifest = load_manifest()
    panels = {}
    for issn, name, prior in JOURNAL_SET:
        print(f"  [{issn}] {name}")
        ws = fetch_journal_works(issn, manifest)
        panels[issn] = {
            'name': name,
            'prior_label': prior,
            'works': ws,
        }
        print(f"    {len(ws)} works")
        save_manifest(manifest)
    save_manifest(manifest)
    return panels


# ═══════════════════════════════════════════════════════════════
# Analysis
# ═══════════════════════════════════════════════════════════════

def classify_journal(works):
    """Empirical authorship-convention classifier: fraction of 3+ author
    papers whose surnames are in non-decreasing alphabetical order.
    Random baseline for a 3-author paper is 1/6 ≈ 0.167."""
    eligible = [w for w in works if w['n_authors'] >= 3]
    if not eligible:
        return None, 0
    alph_frac = sum(w['alpha'] for w in eligible) / len(eligible)
    return alph_frac, len(eligible)


def run_analysis(panels):
    rng = random.Random(SEED)

    # ----- [4/9] Per-journal empirical classification -----
    print("[4/9] Empirical authorship-convention classification...")
    journal_info = {}
    for issn, p in panels.items():
        af, ne = classify_journal(p['works'])
        if af is None or ne < 50:
            continue
        journal_info[issn] = {
            'name': p['name'],
            'prior_label': p['prior_label'],
            'alph_frac': af,
            'n_eligible': ne,
        }
        print(f"  {p['name'][:40]:40s} alph_frac={af:.3f}  (n={ne})")

    # ----- [5/9] Fit per-journal team-size slope (within-journal, winsorised) -----
    print("[5/9] Within-journal team-size slope estimation...")
    for issn, info in journal_info.items():
        ws = panels[issn]['works']
        years = [w['year'] for w in ws if w['n_authors'] <= TEAM_SIZE_CAP]
        sizes = [w['n_authors'] for w in ws if w['n_authors'] <= TEAM_SIZE_CAP]
        if len(set(years)) < MIN_YEARS_PER_JOURNAL or len(years) < MIN_PAPERS_PER_JOURNAL:
            info['qualifies'] = False
            continue
        info['qualifies'] = True
        info['n_papers'] = len(years)
        info['year_span'] = max(years) - min(years)
        info['mean_size'] = mean_val(sizes)
        info['slope_per_yr'] = slope_within(years, sizes)

    qualifying = {k: v for k, v in journal_info.items() if v.get('qualifies')}
    print(f"  {len(qualifying)} journals qualify (>= {MIN_PAPERS_PER_JOURNAL} papers, "
          f">= {MIN_YEARS_PER_JOURNAL} distinct years)")

    # ----- [6/9] Assign empirical convention labels -----
    print("[6/9] Assigning empirical convention strata...")
    strat = {'alph': [], 'contrib': [], 'mixed': []}
    for issn, info in qualifying.items():
        if info['alph_frac'] >= ALPH_FRAC_UPPER:
            info['empirical_label'] = 'alph'
        elif info['alph_frac'] <= ALPH_FRAC_LOWER:
            info['empirical_label'] = 'contrib'
        else:
            info['empirical_label'] = 'mixed'
        strat[info['empirical_label']].append(issn)
    for k, ids in strat.items():
        print(f"  {k}: {len(ids)} journals")

    # Prior/empirical agreement
    agree = sum(
        1 for i, info in qualifying.items()
        if info['empirical_label'] != 'mixed'
        and info['empirical_label'] == info['prior_label']
    )
    print(f"  Prior/empirical agreement: {agree}/"
          f"{sum(1 for i in qualifying.values() if i['empirical_label']!='mixed')}")

    # ----- [7/9] Cross-convention slope comparison + permutation null -----
    print("[7/9] Slope comparison across conventions...")
    a_issns = [i for i in qualifying if qualifying[i]['empirical_label'] == 'alph']
    c_issns = [i for i in qualifying if qualifying[i]['empirical_label'] == 'contrib']
    slopes = [qualifying[i]['slope_per_yr'] for i in a_issns + c_issns]
    labels = ['alph'] * len(a_issns) + ['contrib'] * len(c_issns)

    alph_slopes = [qualifying[i]['slope_per_yr'] for i in a_issns]
    contrib_slopes = [qualifying[i]['slope_per_yr'] for i in c_issns]
    m_alph, lo_alph, hi_alph = bootstrap_slope_ci(alph_slopes, N_BOOT, rng)
    m_contrib, lo_contrib, hi_contrib = bootstrap_slope_ci(contrib_slopes, N_BOOT, rng)
    d_effect = cohens_d(contrib_slopes, alph_slopes) if alph_slopes and contrib_slopes else 0.0
    print(f"  alph    n={len(alph_slopes)} slope={m_alph:.4f} yrs/yr "
          f"[{lo_alph:.4f}, {hi_alph:.4f}]")
    print(f"  contrib n={len(contrib_slopes)} slope={m_contrib:.4f} yrs/yr "
          f"[{lo_contrib:.4f}, {hi_contrib:.4f}]")

    p_perm, obs_delta, perm_deltas = permutation_test_slope_delta(
        slopes, labels, N_PERMS, rng)
    print(f"  Delta (contrib - alph) = {obs_delta:.4f} authors/yr, "
          f"perm p = {p_perm:.4f} (N_perms={N_PERMS})")
    print(f"  Cohen's d (contrib vs alph slopes) = {d_effect:.3f}")

    # Per-journal slope significance + BH-FDR over all qualifying journals
    # Slope p-value via simple t-test against zero using residual SE
    j_pvals, j_ids = [], []
    for issn, info in qualifying.items():
        ws = panels[issn]['works']
        years = [w['year'] for w in ws if w['n_authors'] <= TEAM_SIZE_CAP]
        sizes = [w['n_authors'] for w in ws if w['n_authors'] <= TEAM_SIZE_CAP]
        n = len(years)
        ym = mean_val(years); sm = mean_val(sizes)
        sxx = sum((y - ym) ** 2 for y in years)
        b = info['slope_per_yr']
        a = sm - b * ym
        rss = sum((sizes[k] - (a + b * years[k])) ** 2 for k in range(n))
        if n > 2 and sxx > 0 and rss > 0:
            se = math.sqrt(rss / (n - 2) / sxx)
            t = b / se if se > 0 else 0.0
            # two-sided p from normal approximation (n is always >200)
            z = abs(t)
            p = 2.0 * 0.5 * math.erfc(z / math.sqrt(2))
        else:
            p = 1.0
        info['slope_p'] = p
        j_pvals.append(p); j_ids.append(issn)
    adj_bh, sig_bh = bh_correct(j_pvals, SIGNIFICANCE_THRESHOLD)
    for issn, a_, s_ in zip(j_ids, adj_bh, sig_bh):
        qualifying[issn]['slope_q'] = a_
        qualifying[issn]['slope_sig_bh'] = bool(s_)

    # ----- [8/9] Sensitivity analyses -----
    print("[8/9] Sensitivity analyses...")
    sens = {}

    #  8a: classification threshold variation
    print("  (a) Convention-threshold variation...")
    for upper, lower in [(0.40, 0.20), (0.45, 0.25), (0.50, 0.30),
                         (0.55, 0.35), (0.60, 0.40)]:
        a_s, c_s = [], []
        for i, info in qualifying.items():
            af = info['alph_frac']
            if af >= upper:
                a_s.append(info['slope_per_yr'])
            elif af <= lower:
                c_s.append(info['slope_per_yr'])
        if a_s and c_s:
            sensrng = random.Random(SEED + int(upper * 1000))
            labs = ['alph'] * len(a_s) + ['contrib'] * len(c_s)
            ps, od, _ = permutation_test_slope_delta(
                a_s + c_s, labs, N_PERMS_SENS, sensrng)
            sens[f'thr_{upper}_{lower}'] = {
                'upper': upper, 'lower': lower,
                'n_alph': len(a_s), 'n_contrib': len(c_s),
                'delta': od, 'p_perm': ps,
                'mean_alph': mean_val(a_s), 'mean_contrib': mean_val(c_s),
            }
            print(f"    upper={upper} lower={lower}: "
                  f"n_alph={len(a_s)} n_contrib={len(c_s)} "
                  f"delta={od:.4f} p={ps:.4f}")

    #  8b: time-window variation (split sample)
    print("  (b) Time-window variation...")
    for (lo, hi) in [(1980, 2000), (2000, 2023), (1990, 2010),
                     (1995, 2015), (1980, 2023)]:
        a_s, c_s = [], []
        for i, info in qualifying.items():
            ws = panels[i]['works']
            years = [w['year'] for w in ws
                     if w['n_authors'] <= TEAM_SIZE_CAP and lo <= w['year'] <= hi]
            sizes = [w['n_authors'] for w in ws
                     if w['n_authors'] <= TEAM_SIZE_CAP and lo <= w['year'] <= hi]
            if len(set(years)) < 5 or len(years) < 50:
                continue
            b = slope_within(years, sizes)
            if info['empirical_label'] == 'alph':
                a_s.append(b)
            elif info['empirical_label'] == 'contrib':
                c_s.append(b)
        if a_s and c_s:
            sensrng = random.Random(SEED + lo + hi)
            ps, od, _ = permutation_test_slope_delta(
                a_s + c_s,
                ['alph'] * len(a_s) + ['contrib'] * len(c_s),
                N_PERMS_SENS, sensrng)
            sens[f'yr_{lo}_{hi}'] = {
                'lo': lo, 'hi': hi,
                'n_alph': len(a_s), 'n_contrib': len(c_s),
                'delta': od, 'p_perm': ps,
                'mean_alph': mean_val(a_s), 'mean_contrib': mean_val(c_s),
            }
            print(f"    [{lo}-{hi}]: n_alph={len(a_s)} n_contrib={len(c_s)} "
                  f"delta={od:.4f} p={ps:.4f}")

    #  8c: team-size cap variation
    print("  (c) Team-size cap variation...")
    for cap in [20, 30, 50, 100, 1000]:
        a_s, c_s = [], []
        for i, info in qualifying.items():
            ws = panels[i]['works']
            years = [w['year'] for w in ws if w['n_authors'] <= cap]
            sizes = [w['n_authors'] for w in ws if w['n_authors'] <= cap]
            if len(set(years)) < MIN_YEARS_PER_JOURNAL or len(years) < 100:
                continue
            b = slope_within(years, sizes)
            if info['empirical_label'] == 'alph':
                a_s.append(b)
            elif info['empirical_label'] == 'contrib':
                c_s.append(b)
        if a_s and c_s:
            sensrng = random.Random(SEED + cap)
            ps, od, _ = permutation_test_slope_delta(
                a_s + c_s,
                ['alph'] * len(a_s) + ['contrib'] * len(c_s),
                N_PERMS_SENS, sensrng)
            sens[f'cap_{cap}'] = {
                'cap': cap,
                'n_alph': len(a_s), 'n_contrib': len(c_s),
                'delta': od, 'p_perm': ps,
                'mean_alph': mean_val(a_s), 'mean_contrib': mean_val(c_s),
            }
            print(f"    cap={cap}: n_alph={len(a_s)} n_contrib={len(c_s)} "
                  f"delta={od:.4f} p={ps:.4f}")

    #  8d: drop-one-journal leave-out
    print("  (d) Leave-one-journal-out stability...")
    deltas_loo = []
    all_ids = a_issns + c_issns
    for drop in all_ids:
        sl = [qualifying[i]['slope_per_yr'] for i in all_ids if i != drop]
        lb = [qualifying[i]['empirical_label'] for i in all_ids if i != drop]
        if sl and lb:
            a_loo = [s for s, l in zip(sl, lb) if l == 'alph']
            c_loo = [s for s, l in zip(sl, lb) if l == 'contrib']
            if a_loo and c_loo:
                deltas_loo.append(mean_val(c_loo) - mean_val(a_loo))
    if deltas_loo:
        sens['leave_one_out'] = {
            'n': len(deltas_loo),
            'delta_min': min(deltas_loo),
            'delta_max': max(deltas_loo),
            'delta_mean': mean_val(deltas_loo),
        }
        print(f"    delta range over leave-one-out: "
              f"[{min(deltas_loo):.4f}, {max(deltas_loo):.4f}] (n={len(deltas_loo)})")

    # ----- [9/9] Assemble results -----
    results = {
        'config': {
            'seed': SEED,
            'year_min': YEAR_MIN, 'year_max': YEAR_MAX,
            'alph_frac_upper': ALPH_FRAC_UPPER,
            'alph_frac_lower': ALPH_FRAC_LOWER,
            'min_papers_per_journal': MIN_PAPERS_PER_JOURNAL,
            'min_years_per_journal': MIN_YEARS_PER_JOURNAL,
            'max_works_per_journal': MAX_WORKS_PER_JOURNAL,
            'team_size_cap': TEAM_SIZE_CAP,
            'n_perms': N_PERMS,
            'n_perms_sens': N_PERMS_SENS,
            'n_boot': N_BOOT,
            'ci_level': CI_LEVEL,
            'significance_threshold': SIGNIFICANCE_THRESHOLD,
        },
        'sample': {
            'n_journals_total': len(panels),
            'n_journals_qualifying': len(qualifying),
            'n_alph': len(a_issns),
            'n_contrib': len(c_issns),
            'n_mixed': len(strat['mixed']),
            'total_papers_qualifying': sum(q['n_papers'] for q in qualifying.values()),
            'n_raw_works_total': sum(len(p['works']) for p in panels.values()),
            'prior_empirical_agreement': agree,
        },
        'per_journal': {
            i: {k: v for k, v in info.items()
                if k not in ()}
            for i, info in qualifying.items()
        },
        'main': {
            'alph_mean_slope_per_yr': m_alph,
            'alph_slope_ci_lo': lo_alph, 'alph_slope_ci_hi': hi_alph,
            'contrib_mean_slope_per_yr': m_contrib,
            'contrib_slope_ci_lo': lo_contrib, 'contrib_slope_ci_hi': hi_contrib,
            'delta_contrib_minus_alph': obs_delta,
            'permutation_p_two_sided': p_perm,
            'cohens_d_contrib_vs_alph': d_effect,
            'n_perms': N_PERMS,
        },
        'sensitivity': sens,
    }
    return results


# ═══════════════════════════════════════════════════════════════
# Report generation
# ═══════════════════════════════════════════════════════════════

def generate_report(results):
    try:
        with open('results.json', 'w') as f:
            json.dump(results, f, indent=2)
    except (OSError, TypeError, ValueError) as e:
        print(f"ERROR: Could not write results.json: {e}", file=sys.stderr)
        sys.exit(6)
    m = results['main']
    s = results['sample']
    try:
        _report_fh = open('report.md', 'w')
    except OSError as e:
        print(f"ERROR: Could not open report.md for writing: {e}", file=sys.stderr)
        sys.exit(7)
    _report_fh.close()
    with open('report.md', 'w') as f:
        f.write("# Authorship Convention and Team-Size Inflation — Report\n\n")
        f.write("## Sample\n\n")
        f.write(f"- Journals queried: {s['n_journals_total']}\n")
        f.write(f"- Journals qualifying: {s['n_journals_qualifying']}\n")
        f.write(f"- Alphabetical-convention (empirical): {s['n_alph']}\n")
        f.write(f"- Contribution-order (empirical): {s['n_contrib']}\n")
        f.write(f"- Mixed (dropped): {s['n_mixed']}\n")
        f.write(f"- Total papers after filter: {s['total_papers_qualifying']}\n")
        f.write(f"- A-priori / empirical agreement: {s['prior_empirical_agreement']}\n\n")
        f.write("## Main Finding\n\n")
        f.write("| Stratum | Mean slope (authors/yr) | 95% CI |\n|---|---|---|\n")
        f.write(f"| Alphabetical | {m['alph_mean_slope_per_yr']:.4f} "
                f"| [{m['alph_slope_ci_lo']:.4f}, {m['alph_slope_ci_hi']:.4f}] |\n")
        f.write(f"| Contribution  | {m['contrib_mean_slope_per_yr']:.4f} "
                f"| [{m['contrib_slope_ci_lo']:.4f}, {m['contrib_slope_ci_hi']:.4f}] |\n\n")
        f.write(f"- Delta (contribution − alphabetical): "
                f"{m['delta_contrib_minus_alph']:.4f} authors/yr\n")
        f.write(f"- Permutation p (two-sided, N={m['n_perms']}): "
                f"{m['permutation_p_two_sided']:.4f}\n")
        f.write(f"- Cohen's d: {m['cohens_d_contrib_vs_alph']:.3f}\n\n")
        if 'limitations' in results:
            f.write("## Limitations\n\n")
            for i, lim in enumerate(results['limitations'], 1):
                f.write(f"{i}. {lim}\n")
            f.write("\n")
        f.write("## Assumptions\n\n")
        f.write("- Within-journal team-size slope is approximately linear in year.\n")
        f.write("- The label-permutation null treats journals as exchangeable "
                "within the convention strata.\n")
        f.write("- Surname extraction (last whitespace token of display_name) "
                "is approximately correct across conventions.\n")
    print("  results.json, report.md written")


# ═══════════════════════════════════════════════════════════════
# Main
# ═══════════════════════════════════════════════════════════════

LIMITATIONS = [
    "Convention strata are built from journals of two fields (econ/finance/math vs "
    "biomedicine/clinical). The design does not identify whether effects generalise to "
    "e.g. physics, chemistry, or engineering, which use other conventions again.",
    "The within-journal slope model is linear (OLS of team size on year). Non-linear "
    "growth (e.g. level-shift after 2000) is partially absorbed into a single slope; "
    "the 8(b) time-window sensitivity partially addresses this but is not a formal test.",
    "OpenAlex authorship metadata depends on the publisher's submission of rich author "
    "records. Historical (pre-1990) coverage is less complete in some biomedical "
    "journals, which may bias early-period team-size estimates downward.",
    "Alphabetical-ordering is measured on the surname extracted as the last "
    "whitespace-separated token of display_name. Multi-word surnames and non-English "
    "name conventions can be misclassified, introducing noise (but no systematic bias "
    "across conventions).",
    "The design cannot distinguish honorary authorship from genuine team growth; it "
    "only shows that team growth is markedly slower under a convention that removes "
    "the rank-based reward for adding names.",
    "Analysis does NOT estimate: causal effect of a hypothetical convention change, "
    "per-author productivity, or paper quality differences across conventions.",
]


def main():
    if '--verify' in sys.argv:
        return verify()

    t0 = time.time()
    print("[1/9] Workspace prep...")
    os.makedirs(CACHE_DIR, exist_ok=True)

    print(f"[2/9] Downloading journal panels from OpenAlex ({DATA_URL})...")
    try:
        panels = load_data()
    except RuntimeError as e:
        print(f"ERROR: Could not fetch OpenAlex data: {e}", file=sys.stderr)
        print("Check network access to api.openalex.org and retry.", file=sys.stderr)
        sys.exit(2)
    except Exception as e:
        print(f"ERROR: Unexpected failure during data load: {e}", file=sys.stderr)
        sys.exit(3)

    print("[3/9] Summary of raw panels:")
    total = 0
    for issn, p in panels.items():
        total += len(p['works'])
    print(f"  {total} raw works across {len(panels)} journals")
    if total < 1000:
        print("ERROR: Raw sample too small. Check API connectivity.", file=sys.stderr)
        sys.exit(4)

    try:
        results = run_analysis(panels)
    except Exception as e:
        print(f"ERROR: Analysis step failed: {e}", file=sys.stderr)
        sys.exit(5)
    results['limitations'] = LIMITATIONS
    generate_report(results)

    elapsed = time.time() - t0
    print(f"\nRuntime: {elapsed:.0f}s")
    print("ANALYSIS COMPLETE")


# ═══════════════════════════════════════════════════════════════
# Verification
# ═══════════════════════════════════════════════════════════════

def verify():
    """Verification mode: machine-checkable assertions."""
    print("Running verification...\n")
    ok = fail = 0

    def chk(name, cond):
        nonlocal ok, fail
        status = "PASS" if cond else "FAIL"
        print(f"  {status}: {name}")
        if cond:
            ok += 1
        else:
            fail += 1

    if not os.path.exists('results.json'):
        print("FAIL: results.json not found")
        sys.exit(1)
    with open('results.json') as f:
        r = json.load(f)

    chk("1. results.json is a dict with config/sample/main/sensitivity",
        all(k in r for k in ('config', 'sample', 'main', 'sensitivity')))
    chk("2. report.md exists and non-empty",
        os.path.exists('report.md') and os.path.getsize('report.md') > 100)

    s = r['sample']
    chk("3. At least 10 journals qualify",
        s.get('n_journals_qualifying', 0) >= 10)
    chk("4. Both strata have at least 3 journals",
        s.get('n_alph', 0) >= 3 and s.get('n_contrib', 0) >= 3)
    chk("5. Total papers >= 5000",
        s.get('total_papers_qualifying', 0) >= 5000)

    m = r['main']
    chk("6. Main permutation p in [0,1]",
        0.0 <= m.get('permutation_p_two_sided', -1) <= 1.0)
    chk("7. Cohen's d magnitude < 5 (sanity)",
        abs(m.get('cohens_d_contrib_vs_alph', 99)) < 5.0)
    chk("8. Alphabetical slope CI width > 0",
        m.get('alph_slope_ci_hi', 0) - m.get('alph_slope_ci_lo', 0) > 0)
    chk("9. Contribution slope CI width > 0",
        m.get('contrib_slope_ci_hi', 0) - m.get('contrib_slope_ci_lo', 0) > 0)

    sens = r['sensitivity']
    chk("10. Sensitivity: at least 3 classification thresholds",
        sum(1 for k in sens if k.startswith('thr_')) >= 3)
    chk("11. Sensitivity: at least 3 time windows",
        sum(1 for k in sens if k.startswith('yr_')) >= 3)
    chk("12. Sensitivity: at least 3 team-size caps",
        sum(1 for k in sens if k.startswith('cap_')) >= 3)
    chk("13. Leave-one-out delta recorded",
        'leave_one_out' in sens and 'delta_mean' in sens['leave_one_out'])

    cfg = r['config']
    chk("14. N_perms >= 1000", cfg.get('n_perms', 0) >= 1000)
    chk("15. N_boot >= 1000", cfg.get('n_boot', 0) >= 1000)
    chk("16. Seed recorded", cfg.get('seed') == 42)

    # Config/main/sample consistency
    chk("17. At least one sensitivity delta has same sign or |delta| <= main",
        any(k.startswith('yr_') and abs(sens[k].get('delta', 999)) < 10
            for k in sens))

    # ─── Additional assertions: plausibility, robustness, falsification ───

    alph_mean = m.get('alph_mean_slope_per_yr', 99)
    contrib_mean = m.get('contrib_mean_slope_per_yr', 99)
    main_delta = m.get('delta_contrib_minus_alph', 0.0)

    chk("18. Alphabetical mean slope within plausible bounds (|slope|<1.0 authors/yr)",
        abs(alph_mean) < 1.0)
    chk("19. Contribution mean slope within plausible bounds (|slope|<2.0 authors/yr)",
        abs(contrib_mean) < 2.0)

    # Directional / falsification check: contribution strata grow faster.
    # If this were reversed, the paper's scientific claim would be falsified.
    chk("20. Falsifiable directional prediction: contrib slope > alph slope",
        contrib_mean > alph_mean)

    # CI widths should be meaningful (> 1% of estimate's magnitude or at least
    # 1e-6 in absolute terms — rejects degenerate/collapsed bootstraps).
    alph_ci_w = m.get('alph_slope_ci_hi', 0) - m.get('alph_slope_ci_lo', 0)
    contrib_ci_w = m.get('contrib_slope_ci_hi', 0) - m.get('contrib_slope_ci_lo', 0)
    chk("21. Alphabetical CI width > 1% of |mean slope| (non-degenerate bootstrap)",
        alph_ci_w > 0.01 * max(abs(alph_mean), 1e-6))
    chk("22. Contribution CI width > 1% of |mean slope| (non-degenerate bootstrap)",
        contrib_ci_w > 0.01 * max(abs(contrib_mean), 1e-6))

    # Robustness under sensitivity analyses: all sensitivity deltas should
    # have the same sign as the main delta (directional stability).
    thr_deltas = [sens[k]['delta'] for k in sens if k.startswith('thr_')]
    cap_deltas = [sens[k]['delta'] for k in sens if k.startswith('cap_')]
    chk("23. All threshold sensitivity deltas share sign with main delta",
        thr_deltas and all((d * main_delta) >= 0 for d in thr_deltas))
    chk("24. All team-size-cap sensitivity deltas share sign with main delta",
        cap_deltas and all((d * main_delta) >= 0 for d in cap_deltas))

    # Leave-one-out range sanity: main delta lies within the LOO span.
    loo = sens.get('leave_one_out', {})
    loo_lo = loo.get('delta_min', 0)
    loo_hi = loo.get('delta_max', 0)
    chk("25. Main delta falls within leave-one-out range",
        loo_lo <= main_delta <= loo_hi)

    # Prior/empirical label agreement as an additional data-quality gate.
    pea = s.get('prior_empirical_agreement', 0)
    n_non_mixed = s.get('n_alph', 0) + s.get('n_contrib', 0)
    chk("26. Prior/empirical convention agreement >= 80% of non-mixed journals",
        n_non_mixed > 0 and pea / max(n_non_mixed, 1) >= 0.80)

    # Main permutation uses the full rigor (>=5000 perms) configured above.
    chk("27. Main permutation iterations >= 5000",
        cfg.get('n_perms', 0) >= 5000)

    # Per-journal slopes recorded for every qualifying journal.
    per_j = r.get('per_journal', {})
    n_qualifying = sum(1 for v in per_j.values() if v.get('qualifies'))
    chk("28. Per-journal slopes recorded for every qualifying journal",
        n_qualifying == s.get('n_journals_qualifying', 0))

    # Limitations recorded for the user.
    chk("29. Limitations section recorded with >= 4 items",
        isinstance(r.get('limitations'), list) and len(r['limitations']) >= 4)

    # Negative-control style check: the pre-2000 window (where the theory
    # predicts a weaker effect) should have |delta| <= the post-2000 window's.
    # This is a stronger test than just "sign agrees".
    d_pre = sens.get('yr_1980_2000', {}).get('delta')
    d_post = sens.get('yr_2000_2023', {}).get('delta')
    chk("30. Negative control: pre-2000 |delta| <= post-2000 |delta| (theory-predicted)",
        d_pre is not None and d_post is not None and abs(d_pre) <= abs(d_post))

    # CI_LEVEL was recorded (if present in config as n_perms etc., we just check
    # that CIs exist and lower < upper for all strata).
    chk("31. CIs are well-ordered (lo < hi) for both strata",
        m.get('alph_slope_ci_lo', 1) < m.get('alph_slope_ci_hi', 0)
        and m.get('contrib_slope_ci_lo', 1) < m.get('contrib_slope_ci_hi', 0))

    # ─── Additional data-quality / robustness / falsification checks ───

    # 32. Raw works count recorded and adequate (at least 10,000 works across
    # all queried journals — below this the per-journal panels become unstable).
    chk("32. Raw works total recorded and >= 10,000",
        isinstance(s.get('n_raw_works_total'), int)
        and s.get('n_raw_works_total', 0) >= 10000)

    # 33. CI level recorded in config (reproducibility audit).
    chk("33. CI_LEVEL recorded in config and in (0.5, 1.0)",
        0.5 < cfg.get('ci_level', 0) < 1.0)

    # 34. Significance threshold recorded in config.
    chk("34. significance_threshold recorded in config in (0, 0.5]",
        0.0 < cfg.get('significance_threshold', 0) <= 0.5)

    # 35. Every qualifying per-journal entry has a finite slope recorded
    # (no NaN / inf pollution from degenerate panels).
    per_slopes = [v.get('slope_per_yr') for v in per_j.values()
                  if v.get('qualifies')]
    all_finite = per_slopes and all(
        isinstance(x, (int, float)) and math.isfinite(x) for x in per_slopes)
    chk("35. All qualifying per-journal slopes are finite",
        bool(all_finite))

    # 36. No qualifying per-journal slope magnitude exceeds the plausibility
    # bound (Cohen's d < 5 at the stratum level implies individual slopes
    # should not be wildly out of range either).
    chk("36. Every qualifying per-journal slope magnitude within plausibility bound",
        per_slopes and max(abs(x) for x in per_slopes) < SLOPE_PLAUSIBLE_ABS_MAX)

    # 37. Both strata have at least N_BOOT bootstrap resamples' worth of data
    # (as stored). This is a weak coverage gate: n_qualifying per stratum > 2.
    chk("37. Each stratum has enough journals for a non-trivial bootstrap (>=3)",
        s.get('n_alph', 0) >= 3 and s.get('n_contrib', 0) >= 3)

    # 38. Falsification: under a scrambled label, Cohen's d attenuates — but
    # cheaper to check here: main Cohen's d is of the same sign as main delta.
    chk("38. Sign consistency: Cohen's d and main delta agree in sign",
        (m.get('cohens_d_contrib_vs_alph', 0) * main_delta) > 0)

    # 39. Strong robustness: leave-one-out delta span must not straddle zero
    # (if LOO span crosses zero, effect is driven by a single journal).
    chk("39. Leave-one-out delta span does not straddle zero (robustness)",
        loo_lo > 0 or loo_hi < 0)

    print(f"\n{ok}/{ok + fail} checks passed")
    if fail:
        print("VERIFICATION FAILED")
        sys.exit(1)
    else:
        print("ALL CHECKS PASSED")
        sys.exit(0)


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

**Expected output:** File `analyze.py` written, exit code 0.

## Step 3: Run Analysis

```bash
cd /tmp/claw4s_auto_authorship-inflation-in-alphabetical-authorship-fields-vs-ot && python3 analyze.py
```

**Expected output:**
- Prints `[1/9]` through `[9/9]` progress sections.
- Downloads up to 3,000 works/journal for 30 journals from OpenAlex
  (economics, finance, mathematics, biomedicine, clinical).
- Empirically classifies each journal by the fraction of 3+ author
  papers in alphabetical order.
- Estimates within-journal team-size slopes and compares stratum means.
- Runs 10,000 permutation iterations, 2,000 bootstrap samples, and
  four sensitivity analyses (threshold, year window, team-size cap,
  leave-one-out).
- Writes `results.json` and `report.md`.
- Final line: `ANALYSIS COMPLETE`.
- Runtime: 10–30 minutes on first run, under 2 minutes on rerun
  (cache hit).
- Exit code 0.

## Step 4: Verify Results

```bash
cd /tmp/claw4s_auto_authorship-inflation-in-alphabetical-authorship-fields-vs-ot && python3 analyze.py --verify
```

**Expected output:**
- 39/39 checks passed
- `ALL CHECKS PASSED`
- Exit code 0

## Expected Outputs

| File | Description |
|---|---|
| `results.json` | Structured results: config, sample, main stats, per-journal slopes, sensitivity |
| `report.md` | Human-readable Markdown report with tables |
| `cache/` | Cached OpenAlex responses (SHA256-verified for reproducibility) |
| `data_manifest.json` | SHA256 hashes of all cached files |

## Success Criteria

A successful run satisfies ALL of the following machine-checkable conditions:

1. Script exits 0 on both normal run and `--verify`.
2. At least 10 journals qualify after filters, with >= 3 journals
   per convention stratum.
3. >= 5,000 papers in the qualifying pool.
4. 10,000 permutation iterations on the main null (configurable via `N_PERMS`).
5. 2,000 bootstrap samples on each stratum mean (configurable via `N_BOOT`).
6. Bootstrap CIs have `lo < hi` and width >= 1% of the estimated mean
   magnitude (non-degenerate).
7. Cohen's d magnitude < 5 (statistical sanity).
8. All 39 `--verify` assertions pass, covering: structure, sample size,
   CI well-formedness, plausibility bounds (|slope|<1.0, Cohen's d<5),
   directional robustness across threshold/cap/year-window sensitivities,
   leave-one-out containment, robust LOO delta span not straddling zero,
   negative control (pre-2000 weaker than post-2000), raw works volume
   >= 10,000, every qualifying journal slope finite and in-range,
   per-journal slopes recorded for all qualifiers, config-level
   reproducibility audits (seed, n_perms, n_boot, ci_level,
   significance_threshold), and a limitations list with >= 4 entries.
9. Sensitivity analyses completed across (a) classification threshold,
   (b) time window, (c) team-size cap, (d) leave-one-out.
10. `results.json` includes a `limitations` list with at least 4 entries,
    and `report.md` renders them.

## Failure Conditions

Any of the following indicates a failed run. Exit codes are stable: 0 = success,
1 = any `--verify` assertion failed, 2 = OpenAlex unreachable after retries,
3 = unexpected data-load exception, 4 = raw sample < 1,000 works, 5 = analysis
step raised, 6 = `results.json` could not be written, 7 = `report.md` could
not be opened for writing.

1. **Missing input.** Import errors (the script uses only Python 3.8+
   standard library; no `pip install` is possible).
2. **Network failure.** OpenAlex API unreachable after 4 retries per
   request. The script exits with code 2 and prints a network diagnostic
   message to stderr.
3. **Truncated data.** Raw sample < 1,000 works across all 30 journals
   (exit 4) — strongly suggests the API filter is mis-specified.
4. **Insufficient qualifying journals.** Fewer than 10 journals pass the
   MIN_PAPERS_PER_JOURNAL / MIN_YEARS_PER_JOURNAL filters (assertion 3 fails).
5. **Verification failure.** Any `--verify` assertion fails (script
   exits with code 1). A verification failure must be diagnosed before
   the output is trusted.
6. **Missing outputs.** `results.json` or `report.md` not created (exit
   codes 6/7 in the write path).
7. **Falsified directional prediction.** Contribution slope ≤ alphabetical
   slope (assertion 20). This does not indicate a code bug — it indicates
   the scientific claim itself would need to be re-evaluated.
8. **Data-quality gate.** Prior/empirical convention agreement drops
   below 80% of non-mixed journals (assertion 26) — suggests ISSN list
   errors or an OpenAlex coverage regression.
9. **Implausible effect sizes.** Any stratum mean slope with |slope| ≥ 1.0
   authors/yr or Cohen's d ≥ 5 (assertions 7, 18, 19) indicates the
   pipeline has produced a pathological estimate and the output should
   not be trusted without human review.
10. **Degenerate bootstrap.** CI width ≤ 1% of the estimate magnitude
    (assertions 21–22) — indicates the resampling has collapsed.
11. **Fragile effect.** Leave-one-out delta range straddles zero
    (assertion 39) — the main finding would be driven by a single
    journal and should not be reported as a stratum-level effect.
12. **Corrupted cache.** SHA256 mismatch between a cached response and
    the recorded manifest triggers automatic re-download; if re-download
    also fails, the script exits via the network-failure path.

## Limitations

These limitations are also written verbatim to `results.json['limitations']`
and rendered in `report.md`:

1. **Field scope.** The two strata are built from economics/finance/math
   (alphabetical) vs biomedicine/clinical (contribution). The design does
   not identify effects in physics, chemistry, or engineering.
2. **Linearity assumption.** Within-journal slopes are OLS linear in
   publication year; non-linear growth (e.g. a post-2000 regime shift) is
   partially absorbed, and the 8(b) time-window sensitivity only
   approximately addresses this.
3. **Coverage bias.** OpenAlex depends on publisher metadata. Pre-1990
   coverage is weaker for some biomedical journals and could bias early
   team-size estimates downward.
4. **Surname extraction.** Alphabetical ordering is evaluated on the last
   whitespace-separated token of `display_name`, which is noisy for
   multi-word or non-English surnames. Noise is approximately symmetric
   across conventions and should not induce a systematic stratum effect.
5. **Identification.** The design cannot distinguish gift/honorary
   authorship from genuine collaboration growth; it only shows that
   growth is markedly slower under a convention that removes the
   rank-based reward for adding names.
6. **Out of scope.** The analysis does NOT estimate the causal effect of
   a hypothetical convention change, per-author productivity, or paper
   quality differences across strata.