#!/usr/bin/env python3 """ Do dam structural dimensions predict upstream river miles reconnected by removal? An analysis of the American Rivers Dam Removal Database. We test whether two structural dimensions of removed dams -- height and length -- are associated with the length of upstream river reconnected by the removal (Miles_Restored_US). The analysis is confirmatory in form but was specified after an initial profiling pass of the data (see narrative.md, Pre-registration). All tests, corrections, effect sizes, and robustness checks below were fixed before the reported run and are executed exactly as listed; no test selection depends on this run's outcome. Reproducibility: the script first attempts the canonical figshare download and verifies its md5; on any failure it falls back to a local pinned copy, also md5-gated (hard-fail on mismatch). All random seeds are fixed. Re-running reproduces every number in results.json. Source dataset: American Rivers Dam Removal Database, DOI 10.6084/m9.figshare.5234068.v13, licence CC BY 4.0. """ import hashlib import json import os import sys import urllib.request from pathlib import Path import numpy as np import pandas as pd from scipy import stats # ---------------------------------------------------------------------------- # Configuration and provenance # ---------------------------------------------------------------------------- SEED = 20260101 N_BOOT = 10_000 WINSOR_Q = 0.99 # robustness R2: winsorize US miles at 99th pct CI_LEVEL = 0.95 DATA_FILENAME = "ARDamRemovalList_Figshare_Mar2026.csv" CANONICAL_URL = "https://ndownloader.figshare.com/files/62712199" EXPECTED_MD5 = "0abfc0db243ae3bc5563098a1d8dd505" LOCAL_FALLBACK = "/Users/markhahnel/datasetpapers-handfeed/american-rivers/" + DATA_FILENAME HERE = Path(__file__).resolve().parent FIG_DIR = HERE / "figures" TBL_DIR = HERE / "tables" for d in (FIG_DIR, TBL_DIR): d.mkdir(exist_ok=True) rng = np.random.default_rng(SEED) # ---------------------------------------------------------------------------- # Data acquisition: canonical download, md5-gated local fallback # ---------------------------------------------------------------------------- def _md5(path): h = hashlib.md5() with open(path, "rb") as fh: for chunk in iter(lambda: fh.read(1 << 20), b""): h.update(chunk) return h.hexdigest() def acquire_data(): """Return (local_path, provenance_dict). Canonical download preferred; md5-gated local fallback. Hard-fail if neither yields the pinned md5.""" dest = HERE / DATA_FILENAME prov = {"expected_md5": EXPECTED_MD5, "canonical_url": CANONICAL_URL} # 1) canonical download try: print(f"[data] attempting canonical download: {CANONICAL_URL}", file=sys.stderr) req = urllib.request.Request(CANONICAL_URL) with urllib.request.urlopen(req, timeout=60) as resp: payload = resp.read() dest.write_bytes(payload) got = _md5(dest) if got == EXPECTED_MD5: prov.update(source="figshare_canonical_download", md5_ok=True, md5=got) print("[data] canonical download verified", file=sys.stderr) return dest, prov raise ValueError(f"md5 mismatch on download: {got} != {EXPECTED_MD5}") except Exception as exc: # noqa: BLE001 print(f"[data] canonical download unavailable ({exc}); using local fallback", file=sys.stderr) # 2) md5-gated local fallback if os.path.exists(LOCAL_FALLBACK): got = _md5(LOCAL_FALLBACK) if got != EXPECTED_MD5: raise SystemExit(f"FATAL: local fallback md5 mismatch: {got} != {EXPECTED_MD5}") prov.update(source="local_pinned_fallback", md5_ok=True, md5=got, local_path=LOCAL_FALLBACK) print("[data] local fallback verified", file=sys.stderr) return Path(LOCAL_FALLBACK), prov raise SystemExit("FATAL: no data source available (download failed, no local copy)") # ---------------------------------------------------------------------------- # Numeric coercion helpers # ---------------------------------------------------------------------------- def to_num(series): """Strip commas/units, extract the first signed decimal, coerce to float.""" s = series.astype(str).str.replace(",", "", regex=False) extracted = s.str.extract(r"([-+]?\d*\.?\d+)")[0] return pd.to_numeric(extracted, errors="coerce") def parse_year(series): """Extract a 4-digit year; textual values like 'pre-1959' -> NaN (ambiguous).""" s = series.astype(str).str.strip() # only accept a bare/clean 4-digit year token yr = s.str.extract(r"^(\d{4})$")[0] return pd.to_numeric(yr, errors="coerce") # ---------------------------------------------------------------------------- # Statistics # ---------------------------------------------------------------------------- def spearman(x, y): r = stats.spearmanr(x, y) return float(r.statistic), float(r.pvalue) def bootstrap_spearman_ci(x, y, n_boot=N_BOOT, level=CI_LEVEL, seed=SEED): """Percentile bootstrap CI for Spearman rho (paired resampling).""" x = np.asarray(x, float) y = np.asarray(y, float) n = len(x) local = np.random.default_rng(seed) boots = np.empty(n_boot) for i in range(n_boot): idx = local.integers(0, n, n) boots[i] = stats.spearmanr(x[idx], y[idx]).statistic lo = float(np.nanpercentile(boots, 100 * (1 - level) / 2)) hi = float(np.nanpercentile(boots, 100 * (1 + level) / 2)) return lo, hi def partial_spearman(x, y, z): """Rank-based partial correlation of x,y controlling for z. Rank all three, residualize ranks of x and y on ranks of z (OLS), then Pearson-correlate the residuals. Equivalent to partial Spearman.""" rx = stats.rankdata(x) ry = stats.rankdata(y) rz = stats.rankdata(z) Z = np.column_stack([np.ones_like(rz), rz]) bx, *_ = np.linalg.lstsq(Z, rx, rcond=None) by, *_ = np.linalg.lstsq(Z, ry, rcond=None) ex = rx - Z @ bx ey = ry - Z @ by r, p = stats.pearsonr(ex, ey) return float(r), float(p) def holm_bonferroni(pvals): """Return Holm-corrected p-values in original order.""" p = np.asarray(pvals, float) order = np.argsort(p) m = len(p) corrected = np.empty(m) running = 0.0 for rank, idx in enumerate(order): val = (m - rank) * p[idx] running = max(running, val) corrected[idx] = min(running, 1.0) return corrected.tolist() # ---------------------------------------------------------------------------- # Main # ---------------------------------------------------------------------------- def main(): data_path, provenance = acquire_data() df = pd.read_csv(data_path) n_total = len(df) height = to_num(df["Dam_Height_ft"]) length = to_num(df["Dam_Length_ft"]) us = to_num(df["Miles_Restored_US"]) year = parse_year(df["Year_Removed"]) # Analysis population: positive dam dimension and positive miles restored. # Miles_Restored_US == 0 means "no upstream reconnection reported" and a small # number of dimension==0 entries are placeholders; we require > 0 so the # size->reconnection question is well posed. Recorded transparently below. def complete(a, b): m = a.notna() & b.notna() & (a > 0) & (b > 0) return m m_h = complete(height, us) m_l = complete(length, us) results = { "meta": { "seed": SEED, "n_bootstrap": N_BOOT, "winsor_quantile": WINSOR_Q, "ci_level": CI_LEVEL, "n_records_total": int(n_total), "provenance": provenance, "question": ("Do dam structural dimensions (height, length) predict " "upstream river miles reconnected (Miles_Restored_US)?"), "inclusion": "height>0 (resp. length>0) AND Miles_Restored_US>0", }, "variable_summary": {}, "primary": {}, "robustness": {}, } # --- variable summary (tbl-3) --- for name, ser in [("Dam_Height_ft", height), ("Dam_Length_ft", length), ("Miles_Restored_US", us), ("Year_Removed", year)]: v = ser.dropna() results["variable_summary"][name] = { "n_nonmissing": int(ser.notna().sum()), "pct_missing": round(float(ser.isna().mean() * 100), 2), "median": float(np.median(v)) if len(v) else None, "q25": float(np.percentile(v, 25)) if len(v) else None, "q75": float(np.percentile(v, 75)) if len(v) else None, "min": float(v.min()) if len(v) else None, "max": float(v.max()) if len(v) else None, } # --- PRIMARY test 1: height ~ US --- hx, hy = height[m_h].to_numpy(), us[m_h].to_numpy() rho_h, p_h = spearman(hx, hy) ci_h = bootstrap_spearman_ci(hx, hy, seed=SEED) # --- PRIMARY test 2: length ~ US --- lx, ly = length[m_l].to_numpy(), us[m_l].to_numpy() rho_l, p_l = spearman(lx, ly) ci_l = bootstrap_spearman_ci(lx, ly, seed=SEED + 1) corr_p = holm_bonferroni([p_h, p_l]) results["primary"] = { "height_vs_us": { "test": "Spearman rank correlation", "predictor": "Dam_Height_ft", "outcome": "Miles_Restored_US", "rho": rho_h, "p_raw": p_h, "p_holm": corr_p[0], "ci95": [ci_h[0], ci_h[1]], "n": int(m_h.sum()), }, "length_vs_us": { "test": "Spearman rank correlation", "predictor": "Dam_Length_ft", "outcome": "Miles_Restored_US", "rho": rho_l, "p_raw": p_l, "p_holm": corr_p[1], "ci95": [ci_l[0], ci_l[1]], "n": int(m_l.sum()), }, "multiple_comparison": "Holm-Bonferroni across the two primary tests", } # --- ROBUSTNESS R1: partial Spearman on common subset --- m_both = complete(height, us) & complete(length, us) hh = height[m_both].to_numpy() ll = length[m_both].to_numpy() uu = us[m_both].to_numpy() pr_h, pp_h = partial_spearman(hh, uu, ll) # height~US | length pr_l, pp_l = partial_spearman(ll, uu, hh) # length~US | height results["robustness"]["partial_spearman"] = { "n": int(m_both.sum()), "height_vs_us_given_length": {"partial_rho": pr_h, "p": pp_h}, "length_vs_us_given_height": {"partial_rho": pr_l, "p": pp_l}, "note": "rank-residualization (rank all; OLS-residualize x,y on control; Pearson)", } # --- ROBUSTNESS R2: winsorize US miles at 99th pct, recompute primaries --- cap = float(np.percentile(us.dropna(), WINSOR_Q * 100)) us_w = us.clip(upper=cap) rho_hw, p_hw = spearman(height[m_h].to_numpy(), us_w[m_h].to_numpy()) rho_lw, p_lw = spearman(length[m_l].to_numpy(), us_w[m_l].to_numpy()) results["robustness"]["winsorized_us"] = { "winsor_cap_miles": cap, "height_vs_us": {"rho": rho_hw, "p": p_hw, "n": int(m_h.sum())}, "length_vs_us": {"rho": rho_lw, "p": p_lw, "n": int(m_l.sum())}, } # ------------------------------------------------------------------ # Tables # ------------------------------------------------------------------ tbl1 = pd.DataFrame([ {"predictor": "Dam_Height_ft", "outcome": "Miles_Restored_US", "test": "Spearman", "rho": rho_h, "p_raw": p_h, "p_holm": corr_p[0], "ci95_low": ci_h[0], "ci95_high": ci_h[1], "n": int(m_h.sum())}, {"predictor": "Dam_Length_ft", "outcome": "Miles_Restored_US", "test": "Spearman", "rho": rho_l, "p_raw": p_l, "p_holm": corr_p[1], "ci95_low": ci_l[0], "ci95_high": ci_l[1], "n": int(m_l.sum())}, ]) tbl1.to_csv(TBL_DIR / "tbl-1-primary-correlations.csv", index=False) tbl2 = pd.DataFrame([ {"analysis": "partial_spearman", "predictor": "Dam_Height_ft", "control": "Dam_Length_ft", "outcome": "Miles_Restored_US", "rho": pr_h, "p": pp_h, "n": int(m_both.sum())}, {"analysis": "partial_spearman", "predictor": "Dam_Length_ft", "control": "Dam_Height_ft", "outcome": "Miles_Restored_US", "rho": pr_l, "p": pp_l, "n": int(m_both.sum())}, {"analysis": "winsorized_us_p99", "predictor": "Dam_Height_ft", "control": "", "outcome": "Miles_Restored_US", "rho": rho_hw, "p": p_hw, "n": int(m_h.sum())}, {"analysis": "winsorized_us_p99", "predictor": "Dam_Length_ft", "control": "", "outcome": "Miles_Restored_US", "rho": rho_lw, "p": p_lw, "n": int(m_l.sum())}, ]) tbl2.to_csv(TBL_DIR / "tbl-2-robustness.csv", index=False) vs = results["variable_summary"] tbl3 = pd.DataFrame([ {"variable": k, **v} for k, v in vs.items() ]) tbl3.to_csv(TBL_DIR / "tbl-3-variable-summary.csv", index=False) write_datapackage() # ------------------------------------------------------------------ # Figures # ------------------------------------------------------------------ make_figures(hx, hy, lx, ly, rho_h, ci_h, rho_l, ci_l, pr_h, pr_l, rho_hw, rho_lw) # ------------------------------------------------------------------ # results.json # ------------------------------------------------------------------ with open(HERE / "results.json", "w") as fh: json.dump(results, fh, indent=2) print(json.dumps(results["primary"], indent=2)) print(json.dumps(results["robustness"], indent=2)) return results def write_datapackage(): """Emit a Frictionless tabular-data-package describing the three CSV tables.""" def fld(name, ftype, desc): return {"name": name, "type": ftype, "description": desc} dp = { "profile": "tabular-data-package", "name": "dam-size-vs-upstream-miles-reconnected", "title": "Dam size vs upstream river miles reconnected: derived result tables", "description": ("Derived statistical tables from an analysis of the American " "Rivers Dam Removal Database testing whether dam structural " "dimensions predict upstream river miles reconnected by removal."), "licenses": [{"name": "CC-BY-4.0", "title": "Creative Commons Attribution 4.0", "path": "https://creativecommons.org/licenses/by/4.0/"}], "sources": [{"title": "American Rivers Dam Removal Database", "path": "https://doi.org/10.6084/m9.figshare.5234068.v13"}], "resources": [ {"name": "tbl-1-primary-correlations", "path": "tbl-1-primary-correlations.csv", "format": "csv", "mediatype": "text/csv", "profile": "tabular-data-resource", "title": "Primary Spearman correlations", "schema": {"fields": [ fld("predictor", "string", "Structural dimension tested"), fld("outcome", "string", "Outcome variable (upstream miles reconnected)"), fld("test", "string", "Statistical test"), fld("rho", "number", "Spearman rank correlation coefficient"), fld("p_raw", "number", "Uncorrected two-sided p-value"), fld("p_holm", "number", "Holm-Bonferroni corrected p-value across the two primary tests"), fld("ci95_low", "number", "Lower bound, 95% bootstrap CI for rho (10,000 resamples)"), fld("ci95_high", "number", "Upper bound, 95% bootstrap CI for rho"), fld("n", "integer", "Complete-case sample size"), ]}}, {"name": "tbl-2-robustness", "path": "tbl-2-robustness.csv", "format": "csv", "mediatype": "text/csv", "profile": "tabular-data-resource", "title": "Robustness analyses (partial correlation and winsorized outcome)", "schema": {"fields": [ fld("analysis", "string", "Robustness analysis identifier"), fld("predictor", "string", "Structural dimension tested"), fld("control", "string", "Variable controlled for (partial analysis only)"), fld("outcome", "string", "Outcome variable"), fld("rho", "number", "(Partial) rank correlation coefficient"), fld("p", "number", "Two-sided p-value"), fld("n", "integer", "Complete-case sample size"), ]}}, {"name": "tbl-3-variable-summary", "path": "tbl-3-variable-summary.csv", "format": "csv", "mediatype": "text/csv", "profile": "tabular-data-resource", "title": "Analysis-variable summary and missingness", "schema": {"fields": [ fld("variable", "string", "Variable name"), fld("n_nonmissing", "integer", "Count of non-missing values in full dataset"), fld("pct_missing", "number", "Percent missing in full dataset"), fld("median", "number", "Median of non-missing values"), fld("q25", "number", "25th percentile"), fld("q75", "number", "75th percentile"), fld("min", "number", "Minimum"), fld("max", "number", "Maximum"), ]}}, ], } with open(TBL_DIR / "datapackage.json", "w") as fh: json.dump(dp, fh, indent=2) def make_figures(hx, hy, lx, ly, rho_h, ci_h, rho_l, ci_l, pr_h, pr_l, rho_hw, rho_lw): import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt # Colourblind-safe (Okabe-Ito) BLUE = "#0072B2" ORANGE = "#E69F00" GREEN = "#009E73" GREY = "#555555" def scatter_panel(x, y, rho, ci, xlabel, title, color, fname): fig, ax = plt.subplots(figsize=(6, 4.5)) ax.scatter(x, y, s=22, alpha=0.55, color=color, edgecolor="white", linewidth=0.4) ax.set_yscale("log") ax.set_xlabel(xlabel, fontsize=11) ax.set_ylabel("Upstream miles reconnected (log scale)", fontsize=11) ax.set_title(title, fontsize=11.5, weight="bold") txt = (f"Spearman \u03c1 = {rho:.3f}\n95% CI [{ci[0]:.3f}, {ci[1]:.3f}]\n" f"n = {len(x)}") ax.text(0.97, 0.03, txt, transform=ax.transAxes, ha="right", va="bottom", fontsize=9.5, bbox=dict(boxstyle="round", fc="white", ec=GREY, alpha=0.9)) ax.tick_params(direction="out") for s in ("top", "right"): ax.spines[s].set_visible(False) fig.tight_layout() fig.savefig(FIG_DIR / fname, dpi=200) plt.close(fig) scatter_panel(hx, hy, rho_h, ci_h, "Dam height (ft)", "Dam height vs upstream miles reconnected", BLUE, "fig-1-height-vs-miles.png") scatter_panel(lx, ly, rho_l, ci_l, "Dam length (ft)", "Dam length vs upstream miles reconnected", ORANGE, "fig-2-length-vs-miles.png") # fig-3: forest plot of rho across analyses (with CIs where available) labels = ["Height ~ US\n(primary)", "Length ~ US\n(primary)", "Height ~ US | Length\n(partial)", "Length ~ US | Height\n(partial)", "Height ~ US\n(US winsorized p99)", "Length ~ US\n(US winsorized p99)"] vals = [rho_h, rho_l, pr_h, pr_l, rho_hw, rho_lw] cis = [ci_h, ci_l, None, None, None, None] colors = [BLUE, ORANGE, BLUE, ORANGE, BLUE, ORANGE] fig, ax = plt.subplots(figsize=(7, 5)) ypos = np.arange(len(vals))[::-1] for yp, v, ci, c in zip(ypos, vals, cis, colors): if ci is not None: ax.plot([ci[0], ci[1]], [yp, yp], color=c, lw=2, alpha=0.8) ax.plot(v, yp, "o", color=c, ms=9, markeredgecolor="white", markeredgewidth=0.6) ax.axvline(0, color=GREY, lw=1, ls="--") ax.set_yticks(ypos) ax.set_yticklabels(labels, fontsize=9) ax.set_xlabel("Correlation coefficient (Spearman \u03c1 / partial \u03c1)", fontsize=11) ax.set_title("Association of dam size with upstream miles reconnected\n" "across primary and robustness analyses", fontsize=11.5, weight="bold") ax.set_xlim(-0.25, 0.35) for s in ("top", "right"): ax.spines[s].set_visible(False) fig.tight_layout() fig.savefig(FIG_DIR / "fig-3-robustness.png", dpi=200) plt.close(fig) if __name__ == "__main__": main()