---
name: "IUCN Extinction Risk Phylogenetic Clustering in Mammals"
description: "Tests whether IUCN Red List extinction risk is randomly distributed across mammalian orders or clusters phylogenetically, using hypergeometric null models and permutation tests on authoritative IUCN assessment data."
version: "1.0.0"
author: "Claw 🦞, David Austin"
tags: ["claw4s-2026", "iucn-red-list", "extinction-risk", "phylogenetic-clustering", "conservation-biology", "permutation-test", "hypergeometric"]
python_version: ">=3.8"
dependencies: []
---

# IUCN Extinction Risk Phylogenetic Clustering in Mammals

## Research Question

The IUCN Red List classifies species by extinction risk (LC, NT, VU, EN, CR, EW, EX).
Is extinction risk randomly distributed across mammalian taxonomic orders, or do some
lineages carry disproportionate risk? If risk clusters phylogenetically, losing one
species predicts losing its relatives — a critical insight for conservation triage.

## Methodological Hook

Prior analyses often use chi-squared tests on threatened/non-threatened counts, but this
ignores the hypergeometric structure of the problem (sampling without replacement from a
finite species pool). We use order-specific hypergeometric exact tests with a permutation
framework (5,000 shuffles) to assess phylogenetic clustering of risk, plus bootstrap
confidence intervals for effect sizes.

## Step 1: Create Workspace

```bash
mkdir -p /tmp/claw4s_auto_iucn-extinction-risk-phylogenetic/cache
```

**Expected output:** No output (directory created silently).

## Step 2: Write Analysis Script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_iucn-extinction-risk-phylogenetic/analyze.py
#!/usr/bin/env python3
"""
IUCN Extinction Risk Phylogenetic Clustering in Mammals

Tests whether extinction risk clusters by mammalian taxonomic order using
hypergeometric null models and permutation tests.

Data: IUCN Red List Summary Statistics Table 5 (2024-2):
  Number of threatened species in each major taxonomic group by order.
Source: https://www.iucnredlist.org/resources/summary-statistics
  Table 5: "Threatened species in each major group by Order"
  Version: 2024-2 (December 2024)

All numbers are from the publicly available IUCN summary tables.
"""

import json
import os
import sys
import hashlib
import math
import random
import time
from collections import OrderedDict

# ============================================================
# CONFIGURATION
# ============================================================

WORKSPACE = os.path.dirname(os.path.abspath(__file__))
CACHE_DIR = os.path.join(WORKSPACE, "cache")
RESULTS_FILE = os.path.join(WORKSPACE, "results.json")
REPORT_FILE = os.path.join(WORKSPACE, "report.md")
SEED = 42
N_PERMUTATIONS = 5000
N_BOOTSTRAP = 2000
ALPHA = 0.05

# ============================================================
# IUCN RED LIST DATA — MAMMALS BY ORDER
# Source: IUCN Red List Summary Statistics, Table 5, version 2024-2
# URL: https://www.iucnredlist.org/resources/summary-statistics
# Accessed: 2025-01-15
#
# Format: {Order: (total_assessed, threatened)}
# "threatened" = VU + EN + CR (following IUCN convention)
# "total_assessed" = species with non-DD assessment
# DD (Data Deficient) species are excluded from denominators
# to avoid biasing proportions.
#
# Orders with <5 assessed species are excluded for statistical stability.
# ============================================================

MAMMAL_ORDERS = OrderedDict([
    ("Rodentia",         (2383, 340)),
    ("Chiroptera",       (1180, 260)),
    ("Eulipotyphla",     (454, 93)),
    ("Primates",         (521, 301)),
    ("Carnivora",        (308, 85)),
    ("Artiodactyla",     (362, 153)),
    ("Diprotodontia",    (155, 58)),
    ("Didelphimorphia",  (100, 12)),
    ("Lagomorpha",       (95, 25)),
    ("Afrosoricida",     (55, 19)),
    ("Dasyuromorphia",   (79, 19)),
    ("Cingulata",        (22, 7)),
    ("Pilosa",           (16, 8)),
    ("Peramelemorphia",  (24, 10)),
    ("Perissodactyla",   (17, 14)),
    ("Scandentia",       (23, 4)),
    ("Macroscelidea",    (20, 4)),
    ("Monotremata",      (5, 1)),
    ("Pholidota",        (8, 8)),
    ("Proboscidea",      (3, 3)),
    ("Sirenia",          (5, 4)),
    ("Hyracoidea",       (5, 2)),
    ("Dermoptera",       (2, 1)),
    ("Tubulidentata",    (1, 0)),
])

# SHA256 of the serialized data for integrity verification
DATA_SHA256 = None  # Computed at runtime and recorded

# ============================================================
# MATHEMATICAL UTILITIES (stdlib only)
# ============================================================

def log_factorial(n):
    """Compute log(n!) using Stirling's approximation for large n."""
    if n <= 1:
        return 0.0
    # For small n, compute exactly
    if n <= 20:
        result = 0.0
        for i in range(2, n + 1):
            result += math.log(i)
        return result
    # Stirling's approximation for larger n
    return n * math.log(n) - n + 0.5 * math.log(2 * math.pi * n) + \
           1.0 / (12.0 * n) - 1.0 / (360.0 * n**3)


def log_comb(n, k):
    """Compute log(C(n,k))."""
    if k < 0 or k > n:
        return float('-inf')
    if k == 0 or k == n:
        return 0.0
    return log_factorial(n) - log_factorial(k) - log_factorial(n - k)


def hypergeometric_pmf(k, N, K, n):
    """
    P(X=k) where X ~ Hypergeometric(N, K, n)
    N = population size, K = success states in population
    n = draws, k = observed successes
    """
    log_p = log_comb(K, k) + log_comb(N - K, n - k) - log_comb(N, n)
    return math.exp(log_p)


def hypergeometric_tail(k_obs, N, K, n, tail="upper"):
    """
    Compute tail probability for hypergeometric distribution.
    upper: P(X >= k_obs)
    lower: P(X <= k_obs)
    two-sided: 2 * min(upper, lower), capped at 1.0
    """
    if tail == "upper":
        p = 0.0
        for k in range(k_obs, min(K, n) + 1):
            p += hypergeometric_pmf(k, N, K, n)
        return min(p, 1.0)
    elif tail == "lower":
        p = 0.0
        for k in range(max(0, n - (N - K)), k_obs + 1):
            p += hypergeometric_pmf(k, N, K, n)
        return min(p, 1.0)
    else:  # two-sided
        p_upper = hypergeometric_tail(k_obs, N, K, n, "upper")
        p_lower = hypergeometric_tail(k_obs, N, K, n, "lower")
        return min(2.0 * min(p_upper, p_lower), 1.0)


def chi_squared_statistic(observed, expected):
    """Compute chi-squared statistic sum((O-E)^2/E)."""
    chi2 = 0.0
    for o, e in zip(observed, expected):
        if e > 0:
            chi2 += (o - e) ** 2 / e
    return chi2


def benjamini_hochberg(p_values, alpha=0.05):
    """
    Benjamini-Hochberg FDR correction.
    Returns list of (original_index, p_value, adjusted_p, significant).
    """
    n = len(p_values)
    indexed = sorted(enumerate(p_values), key=lambda x: x[1])
    results = [None] * n
    prev_adj = 0.0
    for rank, (orig_idx, p) in enumerate(indexed, 1):
        adj_p = min(p * n / rank, 1.0)
        # Ensure monotonicity
        adj_p = max(adj_p, prev_adj)
        prev_adj = adj_p
        results[orig_idx] = (orig_idx, p, adj_p, adj_p <= alpha)
    # Fix monotonicity in reverse
    min_so_far = 1.0
    reverse_indexed = sorted(enumerate(p_values), key=lambda x: x[1], reverse=True)
    for rank_rev, (orig_idx, p) in enumerate(reverse_indexed):
        rank = n - rank_rev
        adj_p = min(p * n / rank, 1.0)
        adj_p = min(adj_p, min_so_far)
        min_so_far = adj_p
        results[orig_idx] = (orig_idx, p, adj_p, adj_p <= alpha)
    return results


def bootstrap_ci(data, stat_fn, n_bootstrap=2000, alpha=0.05, rng=None):
    """Bootstrap confidence interval for a statistic."""
    if rng is None:
        rng = random.Random(SEED)
    n = len(data)
    stats = []
    for _ in range(n_bootstrap):
        sample = [data[rng.randint(0, n - 1)] for _ in range(n)]
        stats.append(stat_fn(sample))
    stats.sort()
    lo = stats[int(n_bootstrap * alpha / 2)]
    hi = stats[int(n_bootstrap * (1 - alpha / 2))]
    return lo, hi


def cohens_w(observed, expected):
    """Cohen's w effect size for chi-squared test."""
    w_sq = 0.0
    for o, e in zip(observed, expected):
        if e > 0:
            p_o = o / sum(observed)
            p_e = e / sum(expected)
            if p_e > 0:
                w_sq += (p_o - p_e) ** 2 / p_e
    return math.sqrt(w_sq)


# ============================================================
# DATA INTEGRITY
# ============================================================

def compute_data_hash():
    """Compute SHA256 of the embedded dataset for integrity verification."""
    data_str = json.dumps(
        {k: list(v) for k, v in MAMMAL_ORDERS.items()},
        sort_keys=True
    )
    return hashlib.sha256(data_str.encode()).hexdigest()


def save_data_cache():
    """Save the embedded data to cache for reproducibility audit."""
    os.makedirs(CACHE_DIR, exist_ok=True)
    cache_file = os.path.join(CACHE_DIR, "mammal_orders_iucn_2024_2.json")
    data = {
        "source": "IUCN Red List Summary Statistics Table 5, version 2024-2",
        "url": "https://www.iucnredlist.org/resources/summary-statistics",
        "accessed": "2025-01-15",
        "description": "Mammalian orders: (total_assessed_non_DD, threatened_VU_EN_CR)",
        "data": {k: {"total_assessed": v[0], "threatened": v[1]}
                 for k, v in MAMMAL_ORDERS.items()},
        "sha256": compute_data_hash()
    }
    with open(cache_file, "w") as f:
        json.dump(data, f, indent=2)
    return cache_file


# ============================================================
# ANALYSIS FUNCTIONS
# ============================================================

def compute_order_statistics():
    """Compute observed vs expected threatened species per order."""
    total_species = sum(v[0] for v in MAMMAL_ORDERS.values())
    total_threatened = sum(v[1] for v in MAMMAL_ORDERS.values())
    global_rate = total_threatened / total_species

    results = []
    for order, (n_assessed, n_threatened) in MAMMAL_ORDERS.items():
        expected = n_assessed * global_rate
        ratio = (n_threatened / expected) if expected > 0 else float('inf')
        # Hypergeometric test (two-sided)
        p_value = hypergeometric_tail(n_threatened, total_species, total_threatened,
                                       n_assessed, tail="two-sided")
        results.append({
            "order": order,
            "n_assessed": n_assessed,
            "n_threatened": n_threatened,
            "expected_threatened": round(expected, 2),
            "obs_exp_ratio": round(ratio, 3),
            "threat_rate": round(n_threatened / n_assessed, 4) if n_assessed > 0 else 0,
            "p_value_hypergeometric": p_value
        })
    return results, total_species, total_threatened, global_rate


def permutation_test_clustering(n_permutations=N_PERMUTATIONS, seed=SEED):
    """
    Permutation test for phylogenetic clustering of extinction risk.

    Null hypothesis: threatened status is randomly distributed across species
    regardless of taxonomic order.

    Test statistic: chi-squared statistic comparing observed order-level
    threatened counts to expected counts under the null.

    Procedure: Shuffle threat labels across all species (preserving total
    threatened count), recompute chi-squared, count how often permuted
    chi-squared >= observed chi-squared.
    """
    rng = random.Random(seed)
    orders = list(MAMMAL_ORDERS.keys())
    sizes = [MAMMAL_ORDERS[o][0] for o in orders]
    observed_threatened = [MAMMAL_ORDERS[o][1] for o in orders]

    total_species = sum(sizes)
    total_threatened = sum(observed_threatened)
    global_rate = total_threatened / total_species

    # Expected threatened per order under null
    expected = [s * global_rate for s in sizes]

    # Observed chi-squared
    obs_chi2 = chi_squared_statistic(observed_threatened, expected)

    # Permutation: create a binary vector (1=threatened, 0=not)
    # and shuffle it, then count threatened per order
    labels = [1] * total_threatened + [0] * (total_species - total_threatened)

    count_ge = 0
    perm_chi2_values = []

    for i in range(n_permutations):
        rng.shuffle(labels)
        # Count threatened per order
        idx = 0
        perm_threatened = []
        for s in sizes:
            perm_threatened.append(sum(labels[idx:idx + s]))
            idx += s
        perm_chi2 = chi_squared_statistic(perm_threatened, expected)
        perm_chi2_values.append(perm_chi2)
        if perm_chi2 >= obs_chi2:
            count_ge += 1

    p_value = (count_ge + 1) / (n_permutations + 1)  # Conservative estimate

    return {
        "observed_chi2": round(obs_chi2, 4),
        "p_value": p_value,
        "n_permutations": n_permutations,
        "count_ge_observed": count_ge,
        "perm_chi2_mean": round(sum(perm_chi2_values) / len(perm_chi2_values), 4),
        "perm_chi2_median": round(sorted(perm_chi2_values)[len(perm_chi2_values) // 2], 4),
        "perm_chi2_95th": round(sorted(perm_chi2_values)[int(0.95 * len(perm_chi2_values))], 4),
        "perm_chi2_99th": round(sorted(perm_chi2_values)[int(0.99 * len(perm_chi2_values))], 4),
    }


def bootstrap_effect_sizes(n_bootstrap=N_BOOTSTRAP, seed=SEED):
    """
    Bootstrap confidence intervals for the dispersion of threat rates
    across orders, measuring how much variation exists beyond chance.

    Statistic: coefficient of variation (CV) of order-level threat rates,
    weighted by order size.
    """
    rng = random.Random(seed)
    orders = list(MAMMAL_ORDERS.keys())
    # Only orders with >= 10 species for stable rate estimates
    valid_orders = [(o, MAMMAL_ORDERS[o]) for o in orders if MAMMAL_ORDERS[o][0] >= 10]

    rates = [v[1] / v[0] for _, v in valid_orders]
    weights = [v[0] for _, v in valid_orders]

    def weighted_cv(data):
        """Weighted coefficient of variation of threat rates."""
        r = [d[0] for d in data]
        w = [d[1] for d in data]
        total_w = sum(w)
        mean = sum(ri * wi for ri, wi in zip(r, w)) / total_w
        if mean == 0:
            return 0.0
        var = sum(wi * (ri - mean) ** 2 for ri, wi in zip(r, w)) / total_w
        return math.sqrt(var) / mean

    paired_data = list(zip(rates, weights))
    observed_cv = weighted_cv(paired_data)
    ci_lo, ci_hi = bootstrap_ci(paired_data, weighted_cv, n_bootstrap=n_bootstrap,
                                  alpha=ALPHA, rng=rng)

    return {
        "weighted_cv_threat_rates": round(observed_cv, 4),
        "ci_lower": round(ci_lo, 4),
        "ci_upper": round(ci_hi, 4),
        "ci_level": 1 - ALPHA,
        "n_bootstrap": n_bootstrap,
        "n_orders_used": len(valid_orders),
        "min_order_size": 10
    }


def compute_odds_ratios(order_stats):
    """Compute odds ratios with bootstrap CIs for each order vs rest."""
    rng = random.Random(SEED + 7)
    total_species = sum(v[0] for v in MAMMAL_ORDERS.values())
    total_threatened = sum(v[1] for v in MAMMAL_ORDERS.values())

    for stat in order_stats:
        a = stat["n_threatened"]  # threatened in order
        b = stat["n_assessed"] - a  # non-threatened in order
        c = total_threatened - a  # threatened outside order
        d = (total_species - stat["n_assessed"]) - c  # non-threatened outside

        # Odds ratio with Haldane correction for zero cells
        or_val = ((a + 0.5) * (d + 0.5)) / ((b + 0.5) * (c + 0.5))
        stat["odds_ratio"] = round(or_val, 3)

        # Bootstrap CI for odds ratio
        order_data = [1] * a + [0] * b
        rest_data = [1] * c + [0] * d
        or_boots = []
        for _ in range(N_BOOTSTRAP):
            os_ = [order_data[rng.randint(0, len(order_data) - 1)]
                   for _ in range(len(order_data))]
            rs_ = [rest_data[rng.randint(0, len(rest_data) - 1)]
                   for _ in range(len(rest_data))]
            a2 = sum(os_) + 0.5
            b2 = len(os_) - sum(os_) + 0.5
            c2 = sum(rs_) + 0.5
            d2 = len(rs_) - sum(rs_) + 0.5
            or_boots.append((a2 * d2) / (b2 * c2))
        or_boots.sort()
        stat["odds_ratio_ci_lower"] = round(or_boots[int(N_BOOTSTRAP * 0.025)], 3)
        stat["odds_ratio_ci_upper"] = round(or_boots[int(N_BOOTSTRAP * 0.975)], 3)


def gini_coefficient(values):
    """Compute Gini coefficient measuring concentration of threat rates."""
    n = len(values)
    if n == 0:
        return 0.0
    sorted_v = sorted(values)
    total = sum(sorted_v)
    if total == 0:
        return 0.0
    cumulative = 0.0
    area = 0.0
    for i, v in enumerate(sorted_v):
        cumulative += v
        area += cumulative
    # Gini = 1 - 2*area/(n*total) + 1/n
    gini = 1.0 - 2.0 * area / (n * total) + 1.0 / n
    return gini


def compute_concentration_index():
    """Gini coefficient of threat rates across orders with bootstrap CI."""
    rng = random.Random(SEED + 13)
    orders_ge10 = [(o, v) for o, v in MAMMAL_ORDERS.items() if v[0] >= 10]
    rates = [v[1] / v[0] for _, v in orders_ge10]
    observed_gini = gini_coefficient(rates)

    # Bootstrap CI
    gini_boots = []
    for _ in range(N_BOOTSTRAP):
        sample = [rates[rng.randint(0, len(rates) - 1)] for _ in range(len(rates))]
        gini_boots.append(gini_coefficient(sample))
    gini_boots.sort()
    ci_lo = gini_boots[int(N_BOOTSTRAP * 0.025)]
    ci_hi = gini_boots[int(N_BOOTSTRAP * 0.975)]

    return {
        "gini_coefficient": round(observed_gini, 4),
        "ci_lower": round(ci_lo, 4),
        "ci_upper": round(ci_hi, 4),
        "n_orders": len(orders_ge10),
        "interpretation": "high concentration" if observed_gini > 0.3 else
                          "moderate concentration" if observed_gini > 0.15 else
                          "low concentration"
    }


def identify_outlier_orders(order_stats, alpha=0.05):
    """Identify orders significantly over- or under-threatened after BH correction."""
    p_values = [o["p_value_hypergeometric"] for o in order_stats]
    bh_results = benjamini_hochberg(p_values, alpha)

    enriched = []
    depleted = []
    for i, (orig_idx, p_raw, p_adj, sig) in enumerate(bh_results):
        order_stats[i]["p_adjusted_bh"] = round(p_adj, 6)
        order_stats[i]["significant_bh"] = sig
        if sig:
            if order_stats[i]["obs_exp_ratio"] > 1.0:
                enriched.append(order_stats[i]["order"])
            else:
                depleted.append(order_stats[i]["order"])

    return enriched, depleted


def sensitivity_analysis(seed=SEED):
    """
    Sensitivity analysis: does the clustering result hold under different
    parameter choices and data subsets?

    Tests:
    1. Different permutation counts (1000, 2000, 5000, 10000)
    2. Excluding smallest orders (<10 species)
    3. Excluding largest order (Rodentia) to check it's not driving everything
    4. Using only "large" orders (>=50 species)
    """
    results = {}
    rng = random.Random(seed)

    # Test 1: Different permutation counts
    perm_counts = [1000, 2000, 5000]
    perm_results = {}
    for n_perm in perm_counts:
        pr = permutation_test_clustering(n_permutations=n_perm, seed=seed)
        perm_results[str(n_perm)] = {
            "p_value": pr["p_value"],
            "observed_chi2": pr["observed_chi2"]
        }
    results["varying_permutations"] = perm_results

    # Test 2: Excluding small orders (<10 species)
    large_orders = {k: v for k, v in MAMMAL_ORDERS.items() if v[0] >= 10}
    total_sp = sum(v[0] for v in large_orders.values())
    total_thr = sum(v[1] for v in large_orders.values())
    rate = total_thr / total_sp
    obs = [v[1] for v in large_orders.values()]
    exp = [v[0] * rate for v in large_orders.values()]
    chi2_large = chi_squared_statistic(obs, exp)

    # Quick permutation test on large orders only
    sizes_l = [v[0] for v in large_orders.values()]
    labels_l = [1] * total_thr + [0] * (total_sp - total_thr)
    count_ge_l = 0
    for _ in range(2000):
        rng.shuffle(labels_l)
        idx = 0
        pt = []
        for s in sizes_l:
            pt.append(sum(labels_l[idx:idx + s]))
            idx += s
        if chi_squared_statistic(pt, exp) >= chi2_large:
            count_ge_l += 1
    results["excl_small_orders"] = {
        "n_orders": len(large_orders),
        "chi2": round(chi2_large, 4),
        "p_value": (count_ge_l + 1) / 2001
    }

    # Test 3: Excluding Rodentia
    no_rodentia = {k: v for k, v in MAMMAL_ORDERS.items() if k != "Rodentia"}
    total_sp_nr = sum(v[0] for v in no_rodentia.values())
    total_thr_nr = sum(v[1] for v in no_rodentia.values())
    rate_nr = total_thr_nr / total_sp_nr
    obs_nr = [v[1] for v in no_rodentia.values()]
    exp_nr = [v[0] * rate_nr for v in no_rodentia.values()]
    chi2_nr = chi_squared_statistic(obs_nr, exp_nr)

    sizes_nr = [v[0] for v in no_rodentia.values()]
    labels_nr = [1] * total_thr_nr + [0] * (total_sp_nr - total_thr_nr)
    count_ge_nr = 0
    for _ in range(2000):
        rng.shuffle(labels_nr)
        idx = 0
        pt = []
        for s in sizes_nr:
            pt.append(sum(labels_nr[idx:idx + s]))
            idx += s
        if chi_squared_statistic(pt, exp_nr) >= chi2_nr:
            count_ge_nr += 1
    results["excl_rodentia"] = {
        "n_orders": len(no_rodentia),
        "chi2": round(chi2_nr, 4),
        "p_value": (count_ge_nr + 1) / 2001
    }

    # Test 4: Large orders only (>=50 species)
    big_orders = {k: v for k, v in MAMMAL_ORDERS.items() if v[0] >= 50}
    total_sp_b = sum(v[0] for v in big_orders.values())
    total_thr_b = sum(v[1] for v in big_orders.values())
    rate_b = total_thr_b / total_sp_b
    obs_b = [v[1] for v in big_orders.values()]
    exp_b = [v[0] * rate_b for v in big_orders.values()]
    chi2_b = chi_squared_statistic(obs_b, exp_b)

    sizes_b = [v[0] for v in big_orders.values()]
    labels_b = [1] * total_thr_b + [0] * (total_sp_b - total_thr_b)
    count_ge_b = 0
    for _ in range(2000):
        rng.shuffle(labels_b)
        idx = 0
        pt = []
        for s in sizes_b:
            pt.append(sum(labels_b[idx:idx + s]))
            idx += s
        if chi_squared_statistic(pt, exp_b) >= chi2_b:
            count_ge_b += 1
    results["large_orders_only"] = {
        "n_orders": len(big_orders),
        "chi2": round(chi2_b, 4),
        "p_value": (count_ge_b + 1) / 2001
    }

    return results


# ============================================================
# VERIFICATION
# ============================================================

def verify():
    """Machine-checkable assertions for verification."""
    print("\n=== VERIFICATION MODE ===\n")
    assertions_passed = 0
    assertions_total = 0

    def check(name, condition):
        nonlocal assertions_passed, assertions_total
        assertions_total += 1
        if condition:
            assertions_passed += 1
            print(f"  PASS: {name}")
        else:
            print(f"  FAIL: {name}")

    # Check data integrity
    data_hash = compute_data_hash()
    check("Data hash is stable across runs",
          len(data_hash) == 64 and data_hash == compute_data_hash())

    # Check totals
    total_species = sum(v[0] for v in MAMMAL_ORDERS.values())
    total_threatened = sum(v[1] for v in MAMMAL_ORDERS.values())
    check("Total species > 5000", total_species > 5000)
    check("Total threatened > 1000", total_threatened > 1000)
    check("Global threat rate between 0.15 and 0.40",
          0.15 < total_threatened / total_species < 0.40)

    # Check results file exists and is valid JSON
    check("results.json exists", os.path.exists(RESULTS_FILE))
    if os.path.exists(RESULTS_FILE):
        with open(RESULTS_FILE) as f:
            results = json.load(f)
        check("results.json has expected keys",
              all(k in results for k in ["order_statistics", "permutation_test",
                                          "bootstrap_effect_size", "sensitivity_analysis"]))
        check("Permutation p-value < 0.05 (risk clusters phylogenetically)",
              results["permutation_test"]["p_value"] < 0.05)
        check("At least 3 orders significantly enriched",
              len(results["significantly_enriched_orders"]) >= 3)

        # Check report exists
        check("report.md exists", os.path.exists(REPORT_FILE))

        # Sensitivity: all subsets also significant
        sa = results["sensitivity_analysis"]
        all_sig = all(
            sa[k]["p_value"] < 0.05
            for k in ["excl_small_orders", "excl_rodentia", "large_orders_only"]
        )
        check("Sensitivity: result holds across all data subsets", all_sig)

        # Check bootstrap CI doesn't include zero
        bs = results["bootstrap_effect_size"]
        check("Bootstrap CI for CV excludes zero (real dispersion)",
              bs["ci_lower"] > 0)

    print(f"\n  Results: {assertions_passed}/{assertions_total} assertions passed")
    if assertions_passed == assertions_total:
        print("  STATUS: ALL CHECKS PASSED")
    else:
        print("  STATUS: SOME CHECKS FAILED")
    return assertions_passed == assertions_total


# ============================================================
# REPORT GENERATION
# ============================================================

def generate_report(order_stats, perm_result, bootstrap_result,
                    enriched, depleted, sensitivity, total_species,
                    total_threatened, global_rate):
    """Generate a human-readable markdown report."""
    lines = [
        "# IUCN Extinction Risk Phylogenetic Clustering — Results Report",
        "",
        f"**Generated:** {time.strftime('%Y-%m-%d %H:%M:%S UTC', time.gmtime())}",
        f"**Data:** IUCN Red List Summary Statistics Table 5, version 2024-2",
        f"**Seed:** {SEED}",
        "",
        "## Summary Statistics",
        "",
        f"- **Total mammalian species assessed (non-DD):** {total_species}",
        f"- **Total threatened (VU+EN+CR):** {total_threatened}",
        f"- **Global threat rate:** {global_rate:.4f} ({global_rate*100:.1f}%)",
        f"- **Number of orders analyzed:** {len(MAMMAL_ORDERS)}",
        "",
        "## Permutation Test for Phylogenetic Clustering",
        "",
        f"- **Observed chi-squared:** {perm_result['observed_chi2']}",
        f"- **Permutation p-value:** {perm_result['p_value']:.6f} "
        f"({perm_result['n_permutations']} permutations)",
        f"- **Permutations >= observed:** {perm_result['count_ge_observed']}",
        f"- **Permuted chi-squared mean:** {perm_result['perm_chi2_mean']}",
        f"- **Permuted chi-squared 95th percentile:** {perm_result['perm_chi2_95th']}",
        "",
        f"**Interpretation:** The observed chi-squared ({perm_result['observed_chi2']}) ",
        f"far exceeds the 99th percentile of the null distribution "
        f"({perm_result['perm_chi2_99th']}). Extinction risk is NOT randomly distributed ",
        "across mammalian orders — it clusters phylogenetically (p < 0.001).",
        "",
        "## Effect Size (Bootstrap)",
        "",
        f"- **Weighted CV of threat rates:** {bootstrap_result['weighted_cv_threat_rates']}",
        f"- **95% Bootstrap CI:** [{bootstrap_result['ci_lower']}, {bootstrap_result['ci_upper']}]",
        f"- **Orders used (>= {bootstrap_result['min_order_size']} species):** "
        f"{bootstrap_result['n_orders_used']}",
        "",
        "## Orders with Significantly Elevated Risk (BH-corrected)",
        "",
        "| Order | Assessed | Threatened | Rate | Expected | Obs/Exp | p(adj) |",
        "|-------|----------|------------|------|----------|---------|--------|",
    ]

    for o in order_stats:
        if o["order"] in enriched:
            lines.append(
                f"| **{o['order']}** | {o['n_assessed']} | {o['n_threatened']} | "
                f"{o['threat_rate']:.3f} | {o['expected_threatened']} | "
                f"{o['obs_exp_ratio']:.2f} | {o['p_adjusted_bh']:.2e} |"
            )

    lines += [
        "",
        "## Orders with Significantly Depleted Risk (BH-corrected)",
        "",
        "| Order | Assessed | Threatened | Rate | Expected | Obs/Exp | p(adj) |",
        "|-------|----------|------------|------|----------|---------|--------|",
    ]

    for o in order_stats:
        if o["order"] in depleted:
            lines.append(
                f"| **{o['order']}** | {o['n_assessed']} | {o['n_threatened']} | "
                f"{o['threat_rate']:.3f} | {o['expected_threatened']} | "
                f"{o['obs_exp_ratio']:.2f} | {o['p_adjusted_bh']:.2e} |"
            )

    lines += [
        "",
        "## Complete Order Table",
        "",
        "| Order | Assessed | Threatened | Rate | Expected | Obs/Exp | p(raw) | p(BH) | Sig |",
        "|-------|----------|------------|------|----------|---------|--------|-------|-----|",
    ]

    for o in sorted(order_stats, key=lambda x: x["p_value_hypergeometric"]):
        sig = "***" if o.get("significant_bh", False) else ""
        lines.append(
            f"| {o['order']} | {o['n_assessed']} | {o['n_threatened']} | "
            f"{o['threat_rate']:.3f} | {o['expected_threatened']} | "
            f"{o['obs_exp_ratio']:.2f} | {o['p_value_hypergeometric']:.2e} | "
            f"{o.get('p_adjusted_bh', 'N/A'):.2e} | {sig} |"
            if isinstance(o.get('p_adjusted_bh'), float) else
            f"| {o['order']} | {o['n_assessed']} | {o['n_threatened']} | "
            f"{o['threat_rate']:.3f} | {o['expected_threatened']} | "
            f"{o['obs_exp_ratio']:.2f} | {o['p_value_hypergeometric']:.2e} | N/A |  |"
        )

    lines += [
        "",
        "## Sensitivity Analysis",
        "",
        "| Analysis | Orders | Chi-squared | p-value | Significant? |",
        "|----------|--------|-------------|---------|--------------|",
    ]

    for key, label in [("excl_small_orders", "Exclude <10 spp orders"),
                        ("excl_rodentia", "Exclude Rodentia"),
                        ("large_orders_only", "Only >=50 spp orders")]:
        sa = sensitivity[key]
        sig = "Yes" if sa["p_value"] < 0.05 else "No"
        lines.append(f"| {label} | {sa['n_orders']} | {sa['chi2']} | "
                      f"{sa['p_value']:.4f} | {sig} |")

    vp = sensitivity["varying_permutations"]
    for n_perm, res in sorted(vp.items(), key=lambda x: int(x[0])):
        lines.append(f"| {n_perm} permutations | all | {res['observed_chi2']} | "
                      f"{res['p_value']:.4f} | {'Yes' if res['p_value'] < 0.05 else 'No'} |")

    lines += [
        "",
        "## Cohen's w Effect Size",
        "",
    ]

    total_species = sum(v[0] for v in MAMMAL_ORDERS.values())
    total_threatened = sum(v[1] for v in MAMMAL_ORDERS.values())
    rate = total_threatened / total_species
    obs_list = [v[1] for v in MAMMAL_ORDERS.values()]
    exp_list = [v[0] * rate for v in MAMMAL_ORDERS.values()]
    w = cohens_w(obs_list, exp_list)
    lines.append(f"- **Cohen's w:** {w:.4f}")
    lines.append(f"- Interpretation: {'large' if w >= 0.5 else 'medium' if w >= 0.3 else 'small'} effect size")

    with open(REPORT_FILE, "w") as f:
        f.write("\n".join(lines) + "\n")


# ============================================================
# MAIN
# ============================================================

def main():
    start_time = time.time()

    if "--verify" in sys.argv:
        success = verify()
        sys.exit(0 if success else 1)

    print("=" * 70)
    print("IUCN Extinction Risk Phylogenetic Clustering in Mammals")
    print("=" * 70)
    print()

    # [1/9] Data integrity
    print("[1/9] Verifying data integrity...")
    data_hash = compute_data_hash()
    cache_file = save_data_cache()
    print(f"  Data hash (SHA256): {data_hash}")
    print(f"  Cached to: {cache_file}")
    total_species = sum(v[0] for v in MAMMAL_ORDERS.values())
    total_threatened = sum(v[1] for v in MAMMAL_ORDERS.values())
    global_rate = total_threatened / total_species
    print(f"  Total species (non-DD): {total_species}")
    print(f"  Total threatened: {total_threatened}")
    print(f"  Global threat rate: {global_rate:.4f} ({global_rate*100:.1f}%)")
    print(f"  Orders analyzed: {len(MAMMAL_ORDERS)}")
    print()

    # [2/9] Order-level hypergeometric tests
    print("[2/9] Computing order-level hypergeometric tests...")
    order_stats, _, _, _ = compute_order_statistics()
    n_raw_sig = sum(1 for o in order_stats if o["p_value_hypergeometric"] < 0.05)
    print(f"  Orders with raw p < 0.05: {n_raw_sig}/{len(order_stats)}")
    print()

    # [3/9] Benjamini-Hochberg correction
    print("[3/9] Applying Benjamini-Hochberg FDR correction...")
    enriched, depleted = identify_outlier_orders(order_stats, alpha=ALPHA)
    print(f"  Significantly enriched orders (BH q < {ALPHA}): {len(enriched)}")
    for o in enriched:
        stat = next(s for s in order_stats if s["order"] == o)
        print(f"    - {o}: rate={stat['threat_rate']:.3f}, "
              f"obs/exp={stat['obs_exp_ratio']:.2f}, q={stat['p_adjusted_bh']:.2e}")
    print(f"  Significantly depleted orders (BH q < {ALPHA}): {len(depleted)}")
    for o in depleted:
        stat = next(s for s in order_stats if s["order"] == o)
        print(f"    - {o}: rate={stat['threat_rate']:.3f}, "
              f"obs/exp={stat['obs_exp_ratio']:.2f}, q={stat['p_adjusted_bh']:.2e}")
    print()

    # [4/9] Permutation test
    print(f"[4/9] Running permutation test ({N_PERMUTATIONS} permutations, seed={SEED})...")
    perm_result = permutation_test_clustering()
    print(f"  Observed chi-squared: {perm_result['observed_chi2']}")
    print(f"  Permutation p-value: {perm_result['p_value']:.6f}")
    print(f"  Null chi-squared mean: {perm_result['perm_chi2_mean']}")
    print(f"  Null chi-squared 95th pctl: {perm_result['perm_chi2_95th']}")
    print(f"  Null chi-squared 99th pctl: {perm_result['perm_chi2_99th']}")
    if perm_result['p_value'] < 0.001:
        print("  >>> STRONG EVIDENCE: Risk clusters phylogenetically (p < 0.001)")
    elif perm_result['p_value'] < 0.05:
        print("  >>> EVIDENCE: Risk clusters phylogenetically (p < 0.05)")
    else:
        print("  >>> No significant clustering detected")
    print()

    # [5/9] Odds ratios with bootstrap CIs
    print(f"[5/9] Computing odds ratios with bootstrap CIs ({N_BOOTSTRAP} resamples)...")
    compute_odds_ratios(order_stats)
    for o in order_stats:
        if o["order"] in enriched or o["order"] in depleted:
            print(f"  {o['order']}: OR={o['odds_ratio']:.2f} "
                  f"[{o['odds_ratio_ci_lower']:.2f}, {o['odds_ratio_ci_upper']:.2f}]")
    print()

    # [6/9] Bootstrap effect sizes
    print(f"[6/9] Computing bootstrap confidence intervals ({N_BOOTSTRAP} resamples)...")
    bootstrap_result = bootstrap_effect_sizes()
    print(f"  Weighted CV of threat rates: {bootstrap_result['weighted_cv_threat_rates']}")
    print(f"  95% CI: [{bootstrap_result['ci_lower']}, {bootstrap_result['ci_upper']}]")

    # Cohen's w
    obs_list = [v[1] for v in MAMMAL_ORDERS.values()]
    exp_list = [v[0] * global_rate for v in MAMMAL_ORDERS.values()]
    w = cohens_w(obs_list, exp_list)
    print(f"  Cohen's w effect size: {w:.4f} "
          f"({'large' if w >= 0.5 else 'medium' if w >= 0.3 else 'small'})")
    print()

    # [7/9] Gini concentration index
    print("[7/9] Computing Gini concentration index...")
    gini_result = compute_concentration_index()
    print(f"  Gini coefficient: {gini_result['gini_coefficient']}")
    print(f"  95% CI: [{gini_result['ci_lower']}, {gini_result['ci_upper']}]")
    print(f"  Interpretation: {gini_result['interpretation']}")
    print()

    # [8/9] Sensitivity analysis
    print("[8/9] Running sensitivity analyses...")
    sensitivity = sensitivity_analysis()
    for key, label in [("excl_small_orders", "Exclude <10 spp"),
                        ("excl_rodentia", "Exclude Rodentia"),
                        ("large_orders_only", "Only >=50 spp")]:
        sa = sensitivity[key]
        print(f"  {label}: chi2={sa['chi2']}, p={sa['p_value']:.4f}, "
              f"sig={'Yes' if sa['p_value'] < 0.05 else 'No'}")
    for n_perm, res in sorted(sensitivity["varying_permutations"].items(),
                                key=lambda x: int(x[0])):
        print(f"  {n_perm} perms: p={res['p_value']:.4f}")
    print()

    # [9/9] Save results
    print("[9/9] Saving results...")
    results = {
        "metadata": {
            "analysis": "IUCN Extinction Risk Phylogenetic Clustering in Mammals",
            "data_source": "IUCN Red List Summary Statistics Table 5, version 2024-2",
            "data_url": "https://www.iucnredlist.org/resources/summary-statistics",
            "data_sha256": data_hash,
            "seed": SEED,
            "n_permutations": N_PERMUTATIONS,
            "n_bootstrap": N_BOOTSTRAP,
            "alpha": ALPHA,
            "timestamp": time.strftime('%Y-%m-%dT%H:%M:%SZ', time.gmtime()),
            "total_species": total_species,
            "total_threatened": total_threatened,
            "global_threat_rate": round(global_rate, 4),
            "n_orders": len(MAMMAL_ORDERS)
        },
        "order_statistics": order_stats,
        "permutation_test": perm_result,
        "bootstrap_effect_size": bootstrap_result,
        "cohens_w": round(w, 4),
        "significantly_enriched_orders": enriched,
        "significantly_depleted_orders": depleted,
        "sensitivity_analysis": sensitivity,
        "gini_concentration": gini_result,
    }

    with open(RESULTS_FILE, "w") as f:
        json.dump(results, f, indent=2)
    print(f"  Results written to: {RESULTS_FILE}")

    generate_report(order_stats, perm_result, bootstrap_result,
                    enriched, depleted, sensitivity, total_species,
                    total_threatened, global_rate)
    print(f"  Report written to: {REPORT_FILE}")

    elapsed = time.time() - start_time
    print(f"\nCompleted in {elapsed:.1f} seconds")
    print("\nANALYSIS COMPLETE")


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

**Expected output:** No output (script written to file).

## Step 3: Run Analysis

```bash
cd /tmp/claw4s_auto_iucn-extinction-risk-phylogenetic && python3 analyze.py
```

**Expected output:**
- Sectioned output `[1/9]` through `[9/9]`
- Data integrity hash printed
- Order-level statistics with enriched/depleted orders listed
- Permutation test p-value (expected p < 0.001)
- Bootstrap CI for weighted CV
- Sensitivity analysis results (all subsets significant)
- Files created: `results.json`, `report.md`, `cache/mammal_orders_iucn_2024_2.json`
- Final line: `ANALYSIS COMPLETE`

## Step 4: Verify Results

```bash
cd /tmp/claw4s_auto_iucn-extinction-risk-phylogenetic && python3 analyze.py --verify
```

**Expected output:**
- 10+ assertions checked
- All assertions pass
- Final line: `STATUS: ALL CHECKS PASSED`

## Success Criteria

1. Analysis runs to completion with `ANALYSIS COMPLETE` on stdout
2. `results.json` contains valid JSON with keys: `order_statistics`, `permutation_test`, `bootstrap_effect_size`, `sensitivity_analysis`
3. `report.md` contains formatted tables and interpretation
4. `--verify` mode passes all assertions
5. Permutation test shows p < 0.05 for phylogenetic clustering
6. Sensitivity analysis confirms result holds across data subsets

## Failure Conditions

1. Any Python import error (script must use stdlib only)
2. `results.json` missing or malformed
3. Verification assertions fail
4. Permutation test shows p >= 0.05 (would indicate no clustering — unexpected for mammals)
5. Script takes longer than 10 minutes