Implement final adversarial statistical audit

FINAL_ADVERSARIAL_AUDIT.md

@@ -0,0 +1,105 @@
+# Final adversarial statistical audit
+
+The repository already contains a broad statistical-audit suite in `R/`. Development runs remain useful for checking code paths quickly, but their small Monte Carlo budgets are not final evidence.
+
+`R/run_final_adversarial_audit.R` is the strict final entry point. It prevents development-resolution outputs, hidden per-module caps, or partial failed runs from being mistaken for the final audit.
+
+## Run
+
+```powershell
+Rscript R/run_final_adversarial_audit.R --parallel true --render-report true
+```
+
+Default final budgets are explicit:
+
+| Stage | Replicates |
+|---|---:|
+| dataset bootstrap | 5,000 |
+| primary scan-global null | 5,000 |
+| observed/Ritz bootstrap | 2,000 |
+| specification-grid null per specification | 1,000 |
+| frozen blocked-holdout null per design | 5,000 |
+| declared-family maximum-statistic null | 5,000 |
+| outer null-calibration datasets | 10,000 |
+| injection/recovery per amplitude/frequency | 2,000 |
+| injection decision null | 5,000 |
+| alternative null per model | 5,000 |
+
+The runner aborts on any module failure. Budgets can be overridden explicitly, but the final runner never silently reduces them.
+
+## Adversarial requirements
+
+The final runner executes and reports all of the following:
+
+- **Look-elsewhere correction:** maximum statistic over the complete declared `k` scan.
+- **Full-family multiplicity:** the family now includes the Fe bin grid, every declared sigma/degree/bin preprocessing specification, observed/Ritz representations, and neighbouring ion-II controls. It reports family maximum-statistic, Holm, Bonferroni, BH-FDR, and dependence-robust BY-FDR corrections.
+- **Bootstrap uncertainty:** peak frequency, amplitude, phase, period/scale-ratio uncertainty, and peak-region selection stability.
+- **Null calibration:** synthetic-null false-positive calibration with Monte Carlo intervals.
+- **Frozen holdouts:** `k` is selected only on the training block and locked before confirmatory test-block evaluation. Exploratory test rescans remain separate.
+- **Injection/recovery:** any-peak detection and correct-frequency recovery are reported separately.
+- **Preprocessing/specification sensitivity:** bin counts, Gaussian baseline widths, polynomial degrees, and wavenumber/observed/Ritz representations.
+- **Alternative nulls:** fixed-baseline Poisson, baseline-refit Poisson, conditional multinomial, overdispersed negative binomial, and block-residual baseline-misspecification stress tests.
+- **Failure preservation:** every verdict other than `pass` is copied to `failed_claims.csv`. A failed holdout or predictive comparison cannot be rescued by a strong in-sample result.
+
+The alternative-null family is a **post-signal adversarial stress-test family**, not historical preregistration. The 80% stability threshold used in the final claim matrix is likewise an audit-declared descriptive robustness criterion, not a universal physical threshold.
+
+## Dependence handling
+
+Many Fe analyses reuse the same transition list, so they are dependent. The audit therefore reports several complementary corrections:
+
+- **Holm and Bonferroni:** FWER control that does not require independence.
+- **Benjamini–Yekutieli:** FDR control valid under arbitrary dependence.
+- **Benjamini–Hochberg:** standard FDR reference, retained for comparison.
+- **Family maximum statistic:** useful but explicitly approximate here because it couples marginal per-analysis parametric-null draws rather than simulating one fully joint synthetic line-list experiment.
+
+This distinction is kept in the output rather than collapsing every correction into one number.
+
+## Primary outputs
+
+All tables live under `tables_r/statistical_audit/`.
+
+- `calibrated_global_evidence.csv` — compact primary evidence dashboard.
+- `alternative_null_results.csv` — scan-global p-value under each declared null.
+- `multiple_testing.csv` — full declared-family corrections.
+- `null_calibration.csv` — observed false-positive rate versus nominal alpha.
+- `holdout_results.csv` — locked-frequency blocked validation.
+- `injection_recovery.csv` — power and localisation recovery.
+- `specification_results.csv` — preprocessing/model multiverse.
+- `final_adversarial_claim_matrix.csv` — explicit claim-by-claim verdicts.
+- `failed_claims.csv` — all failed, mixed, inconclusive, and not-established claims.
+- `final_adversarial_run_manifest.csv` — exact budgets, commit, timestamp, and run settings.
+- `final_adversarial_module_status.csv` — per-module success and runtime.
+
+The rendered report is `reports/rendered/nist_final_adversarial_audit.html` when Quarto is available.
+
+## Interpreting global evidence
+
+`calibrated_global_evidence.csv` exposes distinct quantities for the primary Fe II analysis: the full-scan maximum-statistic p-value; family-max, Holm, Bonferroni, BH, and BY multiplicity adjustments; and the worst p-value among the declared alternative null models.
+
+It also reports `largest_reported_global_p`, defined as the maximum across those reported global/multiplicity/alternative-null values. This is deliberately conservative but is only a **diagnostic summary**. It is not a mathematically new combined p-value and must not be described as one.
+
+Every empirical Monte Carlo p-value uses
+
+```text
+p = (r + 1) / (B + 1)
+```
+
+so zero exceedances are reported at the resolution floor `1/(B+1)`, never as `p=0` or as an unsupported extrapolation beyond the simulation budget.
+
+## Claim boundaries
+
+This audit can evaluate the robustness of statistical structure in the declared NIST analysis. It does **not** establish independent experimental confirmation, a unique physical null, a WCT mechanism, a universal atomic law, or NIST endorsement.
+
+## Resolution-mode follow-up is separate
+
+A resolution-dependent winner discovered by the frozen bin-grid audit is investigated by the standalone `R/run_resolution_mode_diagnostics.R` module documented in `RESOLUTION_MODE_DIAGNOSTICS.md`.
+
+That follow-up deliberately does **not** run inside `run_final_adversarial_audit.R` and cannot overwrite its claim matrix. It compares the original fixed-sigma-bin pipeline with an approximately fixed `ell`-space smoothing width, tracks both the Fe primary mode and the previously frozen GWTC `k = 9.602325620315224`, stores the top five local peaks across a dense resolution grid, and performs conservative fixed-frequency null tests across Fe/neighboring-ion sources.
+
+Run it only after preserving the final audit outputs:
+
+```powershell
+Rscript R/run_resolution_mode_diagnostics.R --parallel true --fixed-null-n 5000
+```
+
+The NIST/GWTC frequency comparison is explicitly labeled post-hoc exploratory follow-up even though the GWTC target itself was frozen previously.

R/build_final_adversarial_summary.R

@@ -0,0 +1,252 @@
+#!/usr/bin/env Rscript
+# build_final_adversarial_summary.R
+# ---------------------------------------------------------------------------
+# Cross-cutting final evidence summary for the adversarial NIST audit.
+# Produces a conservative global-evidence dashboard, an explicit claim matrix,
+# and a table containing EVERY non-passing claim. A failed result is never
+# converted to a neutral label merely because another analysis passes.
+# ---------------------------------------------------------------------------
+
+.this <- sub("^--file=", "", commandArgs(FALSE)[grep("^--file=", commandArgs(FALSE))][1])
+if (!exists("write_table", mode = "function")) {
+  source(file.path(if (length(.this) == 0L || is.na(.this)) "R" else dirname(.this),
+                   "audit_utils.R"))
+}
+
+is_true <- function(x) isTRUE(x) || identical(toupper(as.character(x)), "TRUE")
+
+read_audit_table <- function(root, name, required = TRUE) {
+  p <- file.path(root, "tables_r/statistical_audit", name)
+  if (!file.exists(p)) {
+    if (required) stop("required audit table missing: ", p)
+    return(NULL)
+  }
+  utils::read.csv(p, stringsAsFactors = FALSE, check.names = FALSE)
+}
+
+primary_row <- function(df) {
+  if (is.null(df) || nrow(df) == 0L) return(NULL)
+  need <- c("species", "source", "bins", "sigma", "degree")
+  if (all(need %in% names(df))) {
+    hit <- which(df$species == "Fe" & df$source == "wavenumber" & df$bins == 160 &
+                 abs(df$sigma - 6) < 1e-12 & df$degree == 1)
+    if (length(hit) > 0L) return(df[hit[1], , drop = FALSE])
+  }
+  if (all(c("species", "source", "bins") %in% names(df))) {
+    hit <- which(df$species == "Fe" & df$source == "wavenumber" & df$bins == 160)
+    if (length(hit) > 0L) return(df[hit[1], , drop = FALSE])
+  }
+  if ("analysis_id" %in% names(df)) {
+    hit <- which(df$analysis_id == "fe_ion2_wn_bin160")
+    if (length(hit) == 0L) hit <- grep("^fe_ion2_wn.*bin160", df$analysis_id)
+    if (length(hit) > 0L) return(df[hit[1], , drop = FALSE])
+  }
+  df[1, , drop = FALSE]
+}
+
+claim_row <- function(claim, analysis, result, threshold, verdict, limitations, ref) {
+  data.frame(claim = claim, analysis = analysis, result = result,
+             threshold = threshold, verdict = verdict,
+             limitations = limitations, output_reference = ref,
+             stringsAsFactors = FALSE)
+}
+
+holdout_verdict <- function(p, direction, alpha = 0.05) {
+  ok <- is.finite(p) & p <= alpha & direction
+  if (length(ok) == 0L || sum(ok) == 0L) return("fail")
+  if (all(ok)) return("pass")
+  "mixed"
+}
+
+threshold_verdict <- function(value, threshold, higher_is_better = TRUE) {
+  if (!is.finite(value)) return("inconclusive")
+  good <- if (higher_is_better) value >= threshold else value <= threshold
+  if (good) "pass" else "fail"
+}
+
+build_summary <- function(root) {
+  sig <- read_audit_table(root, "significance_results.csv")
+  mult <- read_audit_table(root, "multiple_testing.csv")
+  alt <- read_audit_table(root, "alternative_null_results.csv")
+  cal <- read_audit_table(root, "null_calibration.csv")
+  ho <- read_audit_table(root, "holdout_results.csv")
+  mc <- read_audit_table(root, "model_comparison.csv")
+  bins <- read_audit_table(root, "bin_stability_summary.csv")
+  peak <- read_audit_table(root, "peak_stability.csv")
+  spec <- read_audit_table(root, "specification_results.csv")
+  inj <- read_audit_table(root, "injection_recovery.csv")
+
+  s <- primary_row(sig)
+  m <- primary_row(mult)
+  a <- alt[which.max(alt$scan_global_p), , drop = FALSE]
+  c05 <- cal[which.min(abs(cal$nominal_alpha - 0.05)), , drop = FALSE]
+  mc0 <- mc[mc$model == "M0_smooth_null", , drop = FALSE]
+  mc1 <- mc[mc$model == "M1_smooth_plus_logperiodic", , drop = FALSE]
+  if (nrow(mc0) == 0L || nrow(mc1) == 0L) stop("model comparison rows M0/M1 missing")
+  heldout_gain <- mc1$heldout_loglik_test[1] - mc0$heldout_loglik_test[1]
+
+  primary_peak <- peak[peak$is_primary %in% c(TRUE, "TRUE"), , drop = FALSE]
+  if (nrow(primary_peak) == 0L) primary_peak <- peak[1, , drop = FALSE]
+  spec_frac <- mean(spec$in_reference_region %in% c(TRUE, "TRUE"), na.rm = TRUE)
+  inj0 <- inj[inj$freq_name == "fe_reference" & abs(inj$amplitude) < 1e-12, , drop = FALSE]
+  if (nrow(inj0) == 0L) inj0 <- inj[which.min(abs(inj$amplitude)), , drop = FALSE]
+
+  # Prefer the new adversarial holdout calibration.  Fallback keeps this summary
+  # readable against older development tables but the strict final runner will
+  # regenerate holdout_results.csv before calling us.
+  ho_p <- if ("conservative_fixed_k_test_p" %in% names(ho))
+    ho$conservative_fixed_k_test_p else ho$fixed_k_test_p
+  ho_refit_p <- if ("fixed_k_test_p_refit" %in% names(ho))
+    ho$fixed_k_test_p_refit else rep(NA_real_, nrow(ho))
+  ho_dir <- ho$direction_consistent %in% c(TRUE, "TRUE")
+  ho_ok <- is.finite(ho_p) & ho_p <= 0.05 & ho_dir
+
+  by <- if ("by_fdr" %in% names(m)) m$by_fdr[1] else NA_real_
+  multiplicity_fwer_worst <- max(c(m$family_max_p[1], m$holm_p[1], m$bonferroni_p[1]), na.rm = TRUE)
+  p_candidates <- c(s$global_p[1], m$family_max_p[1], m$holm_p[1],
+                    m$bonferroni_p[1], m$bh_fdr[1], by, a$scan_global_p[1])
+  largest_p <- max(p_candidates, na.rm = TRUE)
+
+  evidence <- data.frame(
+    analysis_id = "fe_ion2_wn_bin160_sigma6_degree1",
+    observed_deltaD = s$global_statistic[1],
+    scan_global_p = s$global_p[1],
+    scan_global_tail_count = s$global_tail_count[1],
+    scan_global_B = s$global_B[1],
+    family_max_p = m$family_max_p[1],
+    holm_p = m$holm_p[1],
+    bonferroni_p = m$bonferroni_p[1],
+    bh_fdr = m$bh_fdr[1],
+    by_fdr = by,
+    family_size = m$family_size[1],
+    worst_fwer_adjusted_p = multiplicity_fwer_worst,
+    worst_alternative_null = a$null_model[1],
+    worst_alternative_null_p = a$scan_global_p[1],
+    worst_alternative_null_B = a$B[1],
+    largest_reported_global_p = largest_p,
+    largest_reported_global_p_note = paste(
+      "maximum of scan-global, family-max, Holm, Bonferroni, BH, BY, and worst alternative-null p-values;",
+      "conservative diagnostic only, not a newly combined formal p-value"),
+    calibration_fpr_alpha_0_05 = c05$observed_fpr[1],
+    calibration_ci_lo = c05$ci_lo[1],
+    calibration_ci_hi = c05$ci_hi[1],
+    calibration_compatible = is_true(c05$compatible[1]),
+    holdout_designs = nrow(ho),
+    holdout_sig_direction_consistent = sum(ho_ok),
+    holdout_median_historical_fixed_p = stats::median(ho$fixed_k_test_p, na.rm = TRUE),
+    holdout_median_refit_p = stats::median(ho_refit_p, na.rm = TRUE),
+    holdout_median_conservative_p = stats::median(ho_p, na.rm = TRUE),
+    heldout_model_loglik_gain_M1_minus_M0 = heldout_gain,
+    bin_reference_fraction = bins$pct_in_reference[1],
+    bootstrap_peak_reference_fraction = primary_peak$pct_in_reference[1],
+    specification_reference_fraction = spec_frac,
+    injection_null_detection_prob = inj0$detection_prob[1],
+    injection_null_n = inj0$n_sims[1],
+    stringsAsFactors = FALSE
+  )
+
+  claims <- list()
+  add <- function(...) claims[[length(claims) + 1L]] <<- claim_row(...)
+
+  add("Fe II full-scan look-elsewhere evidence", "significance_results",
+      sprintf("scan-global p=%.4g (%d/%d exceedances)", s$global_p[1], s$global_tail_count[1], s$global_B[1]),
+      "p <= 0.05", threshold_verdict(s$global_p[1], 0.05, FALSE),
+      "empirical p is bounded by Monte Carlo resolution; zero exceedances are not an exact smaller p",
+      "significance_results.csv")
+
+  add("Fe II full declared-family FWER correction", "global_multiple_testing",
+      sprintf("family-max=%.4g; Holm=%.4g; Bonferroni=%.4g; BH=%.4g; BY=%s; family=%d",
+              m$family_max_p[1], m$holm_p[1], m$bonferroni_p[1], m$bh_fdr[1],
+              if (is.finite(by)) sprintf("%.4g", by) else "NA", m$family_size[1]),
+      "all reported FWER-adjusted p-values <= 0.05",
+      threshold_verdict(multiplicity_fwer_worst, 0.05, FALSE),
+      "Holm/Bonferroni are valid without independence; family-max uses an approximate marginal-null coupling; BY is also reported for arbitrary-dependence FDR",
+      "multiple_testing.csv")
+
+  bad_alt <- alt$null_model[alt$scan_global_p > 0.05]
+  add("Fe II robustness to declared alternative null models", "alternative_null_results",
+      sprintf("worst-case p=%.4g under %s%s", a$scan_global_p[1], a$null_model[1],
+              if (length(bad_alt)) paste0("; p>0.05 under: ", paste(bad_alt, collapse = ", ")) else ""),
+      "all declared null-model scan-global p <= 0.05",
+      if (length(bad_alt) == 0L) "pass" else "fail",
+      "post-signal adversarial stress tests; not historical preregistration and not a claim that one generator is the unique physical null",
+      "alternative_null_results.csv")
+
+  add("Synthetic-null false-positive calibration", "null_calibration",
+      sprintf("FPR@0.05=%.4f [%.4f, %.4f]", c05$observed_fpr[1], c05$ci_lo[1], c05$ci_hi[1]),
+      "nominal 0.05 lies within calibration CI",
+      if (is_true(c05$compatible[1])) "pass" else "fail",
+      "calibration tests the canonical fitted-Poisson synthetic null; alternative-null robustness is reported separately",
+      "null_calibration.csv")
+
+  hverd <- holdout_verdict(ho_p, ho_dir)
+  add("Frozen blocked holdout replication", "run_holdout_replication",
+      sprintf("%d/%d designs have conservative p<=0.05 and positive direction; median conservative p=%.4g",
+              sum(ho_ok), nrow(ho), stats::median(ho_p, na.rm = TRUE)),
+      "strict: every declared blocked design conservative p<=0.05 with positive direction",
+      hverd,
+      "k is frozen from training; conservative p=max(fixed-baseline, baseline-refit null calibration); blocks still come from one line list and are not independent experiments",
+      "holdout_results.csv")
+
+  add("Held-out predictive transfer of smooth-plus-periodic model", "model_comparison",
+      sprintf("held-out log-likelihood gain M1-M0 = %.4f", heldout_gain),
+      "gain > 0", if (heldout_gain > 0) "pass" else "fail",
+      "explicitly preserves the prior failure if the large in-sample model gain does not transfer out of sample",
+      "model_comparison.csv")
+
+  add("Declared Fe II bin stability", "bin_stability_summary",
+      sprintf("%.1f%% of audit-declared bins select the reference region", 100 * bins$pct_in_reference[1]),
+      ">= 80% (audit-declared descriptive threshold)",
+      threshold_verdict(bins$pct_in_reference[1], 0.80, TRUE),
+      "descriptive robustness criterion, not historical preregistration or a universal physical threshold",
+      "bin_stability_summary.csv")
+
+  add("Bootstrap peak-region stability", "peak_stability",
+      sprintf("%.1f%% of primary bootstrap resamples select the 2%% reference region", 100 * primary_peak$pct_in_reference[1]),
+      ">= 80% (audit-declared descriptive threshold)",
+      threshold_verdict(primary_peak$pct_in_reference[1], 0.80, TRUE),
+      "resampling the same dataset does not constitute independent replication",
+      "peak_stability.csv")
+
+  add("Specification multiverse stability", "run_model_sensitivity",
+      sprintf("%.1f%% of declared specifications select the reference region", 100 * spec_frac),
+      ">= 80% (audit-declared descriptive threshold)",
+      threshold_verdict(spec_frac, 0.80, TRUE),
+      "the threshold is a robustness convention; the full specification table remains primary evidence",
+      "specification_results.csv")
+
+  add("Injection null type-I behaviour", "run_injection_recovery",
+      sprintf("A=0 detection probability=%.4f (%d simulations)", inj0$detection_prob[1], inj0$n_sims[1]),
+      "<= 0.05", threshold_verdict(inj0$detection_prob[1], 0.05, FALSE),
+      "finite simulation uncertainty applies; null_calibration.csv is the dedicated type-I assessment",
+      "injection_recovery.csv")
+
+  add("Independent experimental confirmation", "none", "not attempted",
+      "independent dataset/experiment required", "not established",
+      "programming-language parity, alternate preprocessing, and same-line-list holdouts are not independent experiments",
+      "limitations")
+
+  add("WCT physical mechanism or universal atomic law", "none",
+      "not tested by this statistical audit", "physical mechanism evidence required", "not established",
+      "statistical structure alone cannot identify a causal WCT mechanism",
+      "limitations")
+
+  claim_tab <- do.call(rbind, claims)
+  nonpassing <- claim_tab[claim_tab$verdict != "pass", , drop = FALSE]
+  list(evidence = evidence, claims = claim_tab, nonpassing = nonpassing)
+}
+
+main <- function(argv = commandArgs(TRUE)) {
+  root <- audit_repo_root()
+  out <- build_summary(root)
+  td <- file.path(root, "tables_r/statistical_audit")
+  write_table(out$evidence, file.path(td, "calibrated_global_evidence.csv"))
+  write_table(out$claims, file.path(td, "final_adversarial_claim_matrix.csv"))
+  write_table(out$nonpassing, file.path(td, "failed_claims.csv"))
+  cat(sprintf("[final_summary] largest reported global p=%.4g; %d/%d claims non-passing\n",
+              out$evidence$largest_reported_global_p[1], nrow(out$nonpassing), nrow(out$claims)))
+  invisible(out)
+}
+
+.invoked_file <- sub("^--file=", "", commandArgs(FALSE)[grep("^--file=", commandArgs(FALSE))])
+if (length(.invoked_file) > 0L && grepl("build_final_adversarial_summary\\.R$", .invoked_file)) main()