diff --git a/FINAL_ADVERSARIAL_AUDIT.md b/FINAL_ADVERSARIAL_AUDIT.md new file mode 100644 index 0000000..794719e --- /dev/null +++ b/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. diff --git a/R/build_final_adversarial_summary.R b/R/build_final_adversarial_summary.R new file mode 100644 index 0000000..a612d3d --- /dev/null +++ b/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() diff --git a/R/check_audit_dependencies.R b/R/check_audit_dependencies.R index 40eed10..2505528 100644 --- a/R/check_audit_dependencies.R +++ b/R/check_audit_dependencies.R @@ -9,12 +9,12 @@ REQUIRED_PKGS <- c( "jsonlite", "gmp", "dplyr", "tidyr", "purrr", "readr", "stringr", - "tibble", "ggplot2", "scales", "viridisLite" + "tibble", "ggplot2", "scales", "viridisLite", "viridis" ) OPTIONAL_PKGS <- c( # reporting / figures - "patchwork", "ggridges", "ggrepel", "viridis", "gt", "ragg", "svglite", + "patchwork", "ggridges", "ggrepel", "gt", "ragg", "svglite", "quarto", # statistical tooling "boot", "broom", "rsample", "future", "future.apply", "progressr" diff --git a/R/global_multiple_testing.R b/R/global_multiple_testing.R index 168a1a2..08e08e0 100644 --- a/R/global_multiple_testing.R +++ b/R/global_multiple_testing.R @@ -1,64 +1,87 @@ #!/usr/bin/env Rscript # global_multiple_testing.R # --------------------------------------------------------------------------- -# Multiplicity correction across the DECLARED analysis family (read from the -# registry). Three p-value-based procedures plus a family-wise max-statistic -# procedure: -# A. Benjamini-Hochberg FDR -# B. Holm -# C. Bonferroni -# D. Family-wise maximum-statistic +# Multiplicity correction across the FULL DECLARED searched family: +# A. Bonferroni FWER (valid without independence) +# B. Holm FWER (valid without independence) +# C. Benjamini-Hochberg FDR (standard FDR reference; dependence-sensitive) +# D. Benjamini-Yekutieli FDR (valid under arbitrary dependence) +# E. Family-wise maximum-statistic calibration # -# The family-max procedure simulates, for each family-level realisation, a null -# scan of every searched analysis and records the LARGEST statistic across the -# family (over species, ions, bins, sources, baseline settings, ...). Each -# observed statistic is then compared to that family-max null distribution. +# The family includes the complete predefined Fe II preprocessing/specification +# multiverse (bin counts, baseline sigma, polynomial degree, source field) plus +# neighbouring ion-II scans. Duplicate specifications are removed before the +# correction is computed. # -# DEPENDENCE CAVEAT (written to output): the analyses share the same underlying -# NIST line list and are therefore NOT independent. Bonferroni/Holm treat them -# as independent and are conservative; the family-max procedure here combines -# per-analysis parametric nulls and is an approximation of the fully-joint null. -# Both facts are reported rather than hidden. -# -# Writes: -# tables_r/statistical_audit/multiple_testing.csv -# tables_r/statistical_audit/family_max_null.csv -# figures_r/statistical_audit/multiplicity_comparison.png (+ fig12) +# DEPENDENCE CAVEAT: many Fe analyses reuse the same transition list. Holm and +# Bonferroni do not require independence; BY is the dependence-robust FDR +# control. The family-max simulation below combines marginal per-analysis +# parametric null draws and is therefore an approximation to a fully joint +# generative null. That approximation is reported rather than hidden. # --------------------------------------------------------------------------- .this <- sub("^--file=", "", commandArgs(FALSE)[grep("^--file=", commandArgs(FALSE))][1]) if (!exists("emp_p", mode = "function")) source(file.path(if (length(.this) == 0L || is.na(.this)) "R" else dirname(.this), "audit_utils.R")) -# assemble_family(): observed per-analysis global statistics for the searched -# family. Each row: species, ion, source, bins, sigma, degree, observed deltaD. assemble_family <- function(cfg) { root <- audit_repo_root() raw_cache <- new.env() get_lines <- function(sp, src) { - key <- paste(sp, src) + key <- paste(sp, src, sep = "::") if (is.null(raw_cache[[key]])) { path <- file.path(root, sprintf("data/%s_lines.csv", sp)) raw_cache[[key]] <- clean_lines_source(read_nist_csv(path), sp, 2L, src)$lines } raw_cache[[key]] } + fam <- list() - add <- function(id, sp, src, bins, sigma, degree) - fam[[length(fam) + 1L]] <<- data.frame(analysis_id = id, species = sp, source = src, - bins = bins, sigma = sigma, degree = degree, - stringsAsFactors = FALSE) - # bin grid (Fe wn) - for (b in cfg$bins_grid) add(sprintf("fe_ion2_wn_bingrid%d", b), "Fe", "wavenumber", b, cfg$baseline_sigma, cfg$degree) - # sources (Fe) - for (s in c("observed", "ritz")) add(sprintf("fe_ion2_%s_bin%d", s, cfg$bins_primary), "Fe", s, cfg$bins_primary, cfg$baseline_sigma, cfg$degree) - # neighbours + add <- function(sp, src, bins, sigma, degree, family_component) { + fam[[length(fam) + 1L]] <<- data.frame( + species = sp, source = src, bins = as.integer(bins), sigma = as.numeric(sigma), + degree = as.integer(degree), family_component = family_component, + stringsAsFactors = FALSE) + } + + # Predeclared Fe bin grid. + for (b in cfg$bins_grid) + add("Fe", "wavenumber", b, cfg$baseline_sigma, cfg$degree, "bin_grid") + + # Full predeclared preprocessing/model multiverse used by run_model_sensitivity.R. + for (b in c(120L, 160L, 200L)) + for (s in cfg$sigma_grid) + for (d in cfg$degree_grid) + add("Fe", "wavenumber", b, s, d, "specification_multiverse") + + # Alternate line-position representations at the primary canonical settings. + for (src in c("observed", "ritz")) + add("Fe", src, cfg$bins_primary, cfg$baseline_sigma, cfg$degree, "source_replication") + + # Neighbouring ion-II controls at canonical settings. for (sp in c("Cr", "Mn", "Co", "Ni", "Ti")) { lines <- get_lines(sp, "wavenumber") if (nrow(lines) >= cfg$min_lines) - add(sprintf("%s_ion2_wn_bin%d", tolower(sp), cfg$bins_primary), sp, "wavenumber", cfg$bins_primary, cfg$baseline_sigma, cfg$degree) + add(sp, "wavenumber", cfg$bins_primary, cfg$baseline_sigma, cfg$degree, "neighbour_control") } - famdf <- do.call(rbind, fam) - # compute observed deltaD + mu0 for each + + famdf <- unique(do.call(rbind, fam)[, c("species", "source", "bins", "sigma", "degree", "family_component")]) + # If the same numerical specification entered through two declared components, + # collapse it to one hypothesis and retain all component labels. + key <- with(famdf, paste(species, source, bins, sigma, degree, sep = "|")) + groups <- split(seq_len(nrow(famdf)), key) + famdf <- do.call(rbind, lapply(groups, function(ii) { + x <- famdf[ii[1], , drop = FALSE] + x$family_component <- paste(sort(unique(famdf$family_component[ii])), collapse = "+") + x + })) + rownames(famdf) <- NULL + + src_tag <- function(x) ifelse(x == "wavenumber", "wn", x) + famdf$analysis_id <- sprintf( + "%s_ion2_%s_bin%d_sig%s_deg%d", + tolower(famdf$species), src_tag(famdf$source), famdf$bins, + format(famdf$sigma, trim = TRUE, scientific = FALSE), famdf$degree) + k_grid <- audit_k_grid(cfg) obs <- lapply(seq_len(nrow(famdf)), function(i) { g <- famdf[i, ] @@ -71,19 +94,15 @@ assemble_family <- function(cfg) { list(famdf = famdf, obs = obs, k_grid = k_grid) } -# family_max_null(): per-analysis parametric null-max distributions, combined -# into the family-max distribution by independent-draw maximum across analyses. family_max_null <- function(fam, cfg, family_n, parallel = FALSE) { - k_grid <- fam$k_grid per_analysis <- lapply(seq_len(nrow(fam$famdf)), function(i) { o <- fam$obs[[i]]; g <- fam$famdf[i, ] - nd <- null_distribution(o$ell, o$y, o$baseline, o$mu0, k_grid, g$degree, + nd <- null_distribution(o$ell, o$y, o$baseline, o$mu0, fam$k_grid, g$degree, family_n, parallel = parallel) nd$max_vals }) - M <- do.call(cbind, per_analysis) # family_n x n_analyses - fam_max <- apply(M, 1L, max) - list(per_analysis = per_analysis, family_max = fam_max) + M <- do.call(cbind, per_analysis) + list(per_analysis = per_analysis, family_max = apply(M, 1L, max)) } main <- function(argv = commandArgs(TRUE)) { @@ -98,13 +117,11 @@ main <- function(argv = commandArgs(TRUE)) { fam <- assemble_family(cfg) fmn <- family_max_null(fam, cfg, family_n, parallel) + famdf <- fam$famdf fam_max <- fmn$family_max - famdf <- fam$famdf - # per-analysis scan-global p (own null) and family-adjusted p (family-max null) - famdf$scan_global_p <- vapply(seq_len(nrow(famdf)), function(i) { - emp_p(sum(fmn$per_analysis[[i]] >= famdf$observed_stat[i]), family_n) - }, numeric(1)) + famdf$scan_global_p <- vapply(seq_len(nrow(famdf)), function(i) + emp_p(sum(fmn$per_analysis[[i]] >= famdf$observed_stat[i]), family_n), numeric(1)) famdf$family_max_p <- vapply(famdf$observed_stat, function(s) emp_p(sum(fam_max >= s), family_n), numeric(1)) @@ -112,38 +129,47 @@ main <- function(argv = commandArgs(TRUE)) { famdf$bonferroni_p <- pmin(1, p * length(p)) famdf$holm_p <- stats::p.adjust(p, method = "holm") famdf$bh_fdr <- stats::p.adjust(p, method = "BH") + famdf$by_fdr <- stats::p.adjust(p, method = "BY") famdf$family_size <- length(p) - famdf$dependence_note <- "shared NIST line list; corrections approximate, see header" + famdf$dependence_note <- paste( + "shared transition lists create dependence; Holm/Bonferroni valid without independence;", + "BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling") + + # Put the canonical primary Fe II row first for stable downstream selection. + is_primary <- with(famdf, species == "Fe" & source == "wavenumber" & bins == cfg$bins_primary & + sigma == cfg$baseline_sigma & degree == cfg$degree) + famdf <- famdf[order(!is_primary, famdf$species, famdf$source, famdf$bins, famdf$sigma, famdf$degree), ] + rownames(famdf) <- NULL write_table(famdf, file.path(root, "tables_r/statistical_audit/multiple_testing.csv")) write_table(data.frame(family_max_deltaD = fam_max), file.path(root, "tables_r/statistical_audit/family_max_null.csv")) - # comparison figure: raw vs scan-global vs family-adjusted for top analyses long <- rbind( data.frame(analysis_id = famdf$analysis_id, type = "scan-global", p = famdf$scan_global_p), data.frame(analysis_id = famdf$analysis_id, type = "BH-FDR", p = famdf$bh_fdr), + data.frame(analysis_id = famdf$analysis_id, type = "BY-FDR", p = famdf$by_fdr), data.frame(analysis_id = famdf$analysis_id, type = "Holm", p = famdf$holm_p), + data.frame(analysis_id = famdf$analysis_id, type = "Bonferroni", p = famdf$bonferroni_p), data.frame(analysis_id = famdf$analysis_id, type = "family-max", p = famdf$family_max_p)) ord <- famdf$analysis_id[order(famdf$observed_stat)] long$analysis_id <- factor(long$analysis_id, levels = ord) pl <- ggplot2::ggplot(long, ggplot2::aes(p, analysis_id, colour = type, shape = type)) + - ggplot2::geom_point(size = 2.5, alpha = 0.85) + + ggplot2::geom_point(size = 2.1, alpha = 0.8) + ggplot2::geom_vline(xintercept = 0.05, linetype = 3, colour = "grey40") + viridis::scale_colour_viridis(discrete = TRUE, end = 0.85, name = NULL) + - ggplot2::scale_shape_manual(values = c(16, 17, 15, 18), name = NULL) + ggplot2::labs( - title = "Multiplicity-corrected p-values across the declared family", - subtitle = sprintf("family size = %d; family-max from %d realisations. Dotted line = 0.05.", + title = "Multiplicity correction across the full declared analysis family", + subtitle = sprintf("family size = %d; family-max from %d realisations; dotted line = 0.05", length(p), family_n), - x = "p-value", y = NULL, - caption = "Analyses share one line list; corrections are approximate (see table header).") + + x = "adjusted/global p-value", y = NULL, + caption = "Holm/Bonferroni need no independence; BY is dependence-robust FDR; family-max coupling is approximate.") + theme_audit() - save_fig(pl, file.path(root, "figures_r/statistical_audit/multiplicity_comparison.png")) - save_fig(pl, file.path(root, "figures_r/statistical_audit/fig12_multiplicity_comparison.png")) + save_fig(pl, file.path(root, "figures_r/statistical_audit/multiplicity_comparison.png"), height = 10) + save_fig(pl, file.path(root, "figures_r/statistical_audit/fig12_multiplicity_comparison.png"), height = 10) - cat(sprintf("[multiple_testing] family=%d; min BH-FDR=%.4g; min family-max p=%.4g\n", - length(p), min(famdf$bh_fdr), min(famdf$family_max_p))) + cat(sprintf("[multiple_testing] family=%d; primary family-max p=%.4g; primary Holm=%.4g; primary BY=%.4g\n", + length(p), famdf$family_max_p[1], famdf$holm_p[1], famdf$by_fdr[1])) invisible(famdf) } diff --git a/R/model_comparison.R b/R/model_comparison.R index 7eb9f02..7ce9099 100644 --- a/R/model_comparison.R +++ b/R/model_comparison.R @@ -15,6 +15,11 @@ # multiplicity correction. M0 is a smooth statistical baseline, NOT a complete # physical atomic model. # +# Held-out prediction uses the polynomial centering/scaling learned on the +# training block and reuses that exact transform on the test block. Recomputing +# the transform on the test block would put the training coefficients in a +# different coordinate system and invalidate the out-of-sample score. +# # Writes: # tables_r/statistical_audit/model_comparison.csv # figures_r/statistical_audit/model_comparison.png (+ fig11) @@ -32,6 +37,20 @@ aic_bic <- function(loglik, npar, n) { c(aic = aic, aicc = aicc, bic = bic) } +# Explicit polynomial transform for honest train -> test prediction. +poly_transform <- function(ell) { + center <- mean(ell) + z <- ell - center + scale <- sqrt(mean(z * z)) + if (!is.finite(scale) || scale <= 0) scale <- 1 + list(center = center, scale = scale) +} + +design_poly_from_transform <- function(ell, degree, transform) { + z <- (ell - transform$center) / transform$scale + do.call(cbind, lapply(0:degree, function(d) z^d)) +} + # fit_models(): fit M0 and M1 at a given k on (ell,y,baseline). Returns metrics. fit_models <- function(ell, y, baseline, degree, k) { X0 <- design_poly(ell, degree) @@ -49,32 +68,38 @@ fit_models <- function(ell, y, baseline, degree, k) { } # heldout_predictive(): blocked split (lower half train / upper half test). -# k is estimated on the training block and LOCKED, then both models are refit -# on training and scored by predictive Poisson log-likelihood on the test block. +# k is estimated on the training block and LOCKED, then both models are fit on +# training and scored by predictive Poisson log-likelihood on the test block. +# The polynomial transform is estimated on training only and frozen. heldout_predictive <- function(ell, y, baseline, degree, k_grid) { n <- length(y) mid <- floor(n / 2) tr <- seq_len(mid); te <- seq.int(mid + 1L, n) + # lock k on training sk <- scan_k(ell[tr], y[tr], baseline[tr], k_grid, degree) k_lock <- sk$best$k_best - X0tr <- design_poly(ell[tr], degree) + + tf <- poly_transform(ell[tr]) + X0tr <- design_poly_from_transform(ell[tr], degree, tf) X1tr <- cbind(X0tr, cos(k_lock * ell[tr]), sin(k_lock * ell[tr])) f0 <- fit_poisson_loglinear(y[tr], baseline[tr], X0tr) f1 <- fit_poisson_loglinear(y[tr], baseline[tr], X1tr) - # predict on test using training betas (design built on test ell, centered by - # the training transform is approximated by recomputing design on test ell) - predict_mu <- function(beta, ell_te, base_te, degree, k = NULL) { - X0 <- design_poly(ell_te, degree) + + predict_mu <- function(beta, ell_te, base_te, degree, transform, k = NULL) { + X0 <- design_poly_from_transform(ell_te, degree, transform) X <- if (is.null(k)) X0 else cbind(X0, cos(k * ell_te), sin(k * ell_te)) eta <- pmin(pmax(as.numeric(X %*% beta), -10), 10) pmax(base_te * exp(eta), EPS) } - mu0_te <- predict_mu(f0$beta, ell[te], baseline[te], degree) - mu1_te <- predict_mu(f1$beta, ell[te], baseline[te], degree, k_lock) + mu0_te <- predict_mu(f0$beta, ell[te], baseline[te], degree, tf) + mu1_te <- predict_mu(f1$beta, ell[te], baseline[te], degree, tf, k_lock) + list(k_lock = k_lock, ll0_test = poisson_loglik(y[te], mu0_te), - ll1_test = poisson_loglik(y[te], mu1_te)) + ll1_test = poisson_loglik(y[te], mu1_te), + train_center = tf$center, + train_scale = tf$scale) } main <- function(argv = commandArgs(TRUE)) { @@ -103,8 +128,12 @@ main <- function(argv = commandArgs(TRUE)) { tab$deltaBIC <- tab$BIC - min(tab$BIC) tab$selected_k <- k tab$heldout_k_lock <- ho$k_lock + tab$heldout_train_center <- ho$train_center + tab$heldout_train_scale <- ho$train_scale tab$deviance_difference <- fm$deviance0 - fm$deviance1 - tab$note <- "k selected by scan; AIC/BIC do not fully account for look-elsewhere" + tab$note <- paste( + "k selected by scan; AIC/BIC do not fully account for look-elsewhere;", + "held-out polynomial transform frozen from training block") write_table(tab, file.path(root, "tables_r/statistical_audit/model_comparison.csv")) plotdf <- data.frame( diff --git a/R/run_alternative_nulls.R b/R/run_alternative_nulls.R new file mode 100644 index 0000000..95d6a4a --- /dev/null +++ b/R/run_alternative_nulls.R @@ -0,0 +1,214 @@ +#!/usr/bin/env Rscript +# run_alternative_nulls.R +# --------------------------------------------------------------------------- +# Adversarial scan-global null-model stress test for the primary Fe II scan. +# +# The canonical audit uses a fitted-smooth Poisson parametric bootstrap with the +# observed Gaussian baseline held fixed. This module deliberately broadens the +# null family and, where indicated, re-estimates the smoothing baseline inside +# every synthetic replicate. The goal is not to select the most favourable +# null; it is to report the WORST (largest) calibrated scan-global p-value across +# a declared set of reasonable alternatives. +# +# Declared nulls: +# 1. poisson_fixed_baseline -- canonical fitted-Poisson bootstrap +# 2. poisson_refit_baseline -- same generator, but refit smooth baseline +# 3. conditional_multinomial -- fixes total line count, refits baseline +# 4. negative_binomial_refit -- overdispersed count null, refits baseline +# 5. block_residual_refit -- block-permuted Pearson residual structure, +# converted to a rate perturbation, then +# Poisson sampled and baseline-refit +# +# The block-residual construction is an adversarial diagnostic rather than an +# exact generative model. It preserves short-range residual chunks while +# destroying their global ordering, so it tests sensitivity to local baseline +# misspecification/correlation without preserving a coherent long-range phase. +# +# Writes: +# tables_r/statistical_audit/alternative_null_results.csv +# tables_r/statistical_audit/alternative_null_maxima.csv +# figures_r/statistical_audit/fig17_alternative_nulls.png +# --------------------------------------------------------------------------- + +.this <- sub("^--file=", "", commandArgs(FALSE)[grep("^--file=", commandArgs(FALSE))][1]) +if (!exists("emp_p", mode = "function")) { + source(file.path(if (length(.this) == 0L || is.na(.this)) "R" else dirname(.this), + "audit_utils.R")) +} + +estimate_nb_size <- function(y, mu0) { + mu0 <- pmax(as.numeric(mu0), EPS) + y <- as.numeric(y) + excess <- sum((y - mu0)^2) - sum(mu0) + if (!is.finite(excess) || excess <= EPS) return(Inf) + size <- sum(mu0^2) / excess + if (!is.finite(size) || size <= EPS) return(Inf) + size +} + +block_residual_mean <- function(y, mu0, block_len = 8L) { + y <- as.numeric(y); mu0 <- pmax(as.numeric(mu0), EPS) + block_len <- max(1L, as.integer(block_len)) + r <- (y - mu0) / sqrt(mu0) + r <- r - mean(r, na.rm = TRUE) + blocks <- split(seq_along(r), ceiling(seq_along(r) / block_len)) + ord <- sample.int(length(blocks), length(blocks), replace = FALSE) + rp <- unlist(lapply(blocks[ord], function(ii) r[ii]), use.names = FALSE) + rp <- rp[seq_along(r)] + mu_star <- pmax(mu0 + rp * sqrt(mu0), EPS) + # Preserve the fitted expected total count so only residual structure changes. + if (sum(mu_star) > 0) mu_star <- mu_star * sum(mu0) / sum(mu_star) + mu_star +} + +simulate_counts <- function(model, y, mu0, nb_size, block_len) { + if (model %in% c("poisson_fixed_baseline", "poisson_refit_baseline")) { + return(stats::rpois(length(mu0), mu0)) + } + if (model == "conditional_multinomial") { + prob <- pmax(mu0, EPS); prob <- prob / sum(prob) + return(as.numeric(stats::rmultinom(1L, size = as.integer(round(sum(y))), prob = prob))) + } + if (model == "negative_binomial_refit") { + if (!is.finite(nb_size)) return(stats::rpois(length(mu0), mu0)) + return(stats::rnbinom(length(mu0), mu = mu0, size = nb_size)) + } + if (model == "block_residual_refit") { + mu_star <- block_residual_mean(y, mu0, block_len) + return(stats::rpois(length(mu_star), mu_star)) + } + stop("unknown null model: ", model) +} + +summarize_null_model <- function(model, maxima, observed_stat, B, nb_size = NA_real_, + block_len = NA_integer_) { + tail <- sum(maxima >= observed_stat, na.rm = TRUE) + ci <- binom_ci(tail, B) + data.frame( + null_model = model, + B = as.integer(B), + observed_deltaD = observed_stat, + tail_count = as.integer(tail), + scan_global_p = emp_p(tail, B), + resolution_floor = resolution_floor(B), + exceedance_rate_ci_lo = unname(ci["lower"]), + exceedance_rate_ci_hi = unname(ci["upper"]), + null_median = stats::median(maxima, na.rm = TRUE), + null_q95 = unname(stats::quantile(maxima, 0.95, na.rm = TRUE)), + null_q99 = unname(stats::quantile(maxima, 0.99, na.rm = TRUE)), + null_max = max(maxima, na.rm = TRUE), + nb_size = nb_size, + block_len = block_len, + stringsAsFactors = FALSE + ) +} + +run_one_model <- function(model, res, k_grid, B, seed, block_len = 8L, + parallel = FALSE) { + setup_rng(seed) + seeds <- sample.int(.Machine$integer.max, B) + nb_size <- estimate_nb_size(res$y, res$mu0) + refit <- model != "poisson_fixed_baseline" + one <- function(i) { + set.seed(seeds[i]) + y0 <- simulate_counts(model, res$y, res$mu0, nb_size, block_len) + base0 <- if (refit) pmax(gaussian_filter_nearest(y0, res$baseline_sigma), EPS) else res$baseline + sk <- scan_k(res$ell, y0, base0, k_grid, res$degree) + sk$best$deltaD + } + vals <- unlist(audit_lapply(seq_len(B), one, parallel = parallel), use.names = FALSE) + list( + summary = summarize_null_model(model, vals, res$best$deltaD, B, + nb_size = if (model == "negative_binomial_refit") nb_size else NA_real_, + block_len = if (model == "block_residual_refit") block_len else NA_integer_), + maxima = data.frame(null_model = model, replicate = seq_len(B), + max_deltaD = vals, stringsAsFactors = FALSE) + ) +} + +model_metadata <- function() { + data.frame( + null_model = c("poisson_fixed_baseline", "poisson_refit_baseline", + "conditional_multinomial", "negative_binomial_refit", + "block_residual_refit"), + description = c( + "Fitted smooth Poisson null; observed smoothing baseline held fixed (canonical audit null).", + "Fitted smooth Poisson null; Gaussian smoothing baseline re-estimated in every replicate.", + "Conditional multinomial null with total line count fixed; smoothing baseline re-estimated.", + "Negative-binomial overdispersed count null using method-of-moments size; smoothing baseline re-estimated.", + "Block-permuted Pearson-residual rate perturbation followed by Poisson sampling; smoothing baseline re-estimated." + ), + interpretation = c( + "reference", + "tests baseline-estimation uncertainty", + "conditions on observed total count", + "tests extra-Poisson variance", + "adversarial local-structure / baseline-misspecification diagnostic" + ), stringsAsFactors = FALSE) +} + +plot_alternative_nulls <- function(maxima, results) { + lab <- setNames(sprintf("p=%.4g", results$scan_global_p), results$null_model) + maxima$null_model <- factor(maxima$null_model, levels = rev(results$null_model)) + ggplot2::ggplot(maxima, ggplot2::aes(max_deltaD, null_model)) + + ggplot2::geom_boxplot(outlier.shape = NA, width = 0.65, fill = "grey88") + + ggplot2::geom_vline(xintercept = results$observed_deltaD[1], linetype = 2, + colour = species_colour("Fe")) + + ggplot2::stat_summary(fun = stats::median, geom = "point", size = 2) + + ggplot2::geom_text(data = results, + ggplot2::aes(x = Inf, y = factor(null_model, levels = rev(results$null_model)), + label = lab[null_model]), inherit.aes = FALSE, + hjust = 1.05, vjust = -0.5, size = 3) + + ggplot2::coord_cartesian(clip = "off") + + ggplot2::labs( + title = "Adversarial alternative-null scan-global calibration", + subtitle = "Dashed line = observed Fe II maximum statistic; p-values use (r+1)/(B+1).", + x = "null replicate maximum deltaD across full k scan", y = NULL, + caption = "Largest p across declared nulls is the conservative summary; zero exceedances remain resolution-limited.") + + theme_audit() +} + +main <- function(argv = commandArgs(TRUE)) { + B <- 5000L; block_len <- 8L; seed <- 20260517L + i <- which(argv == "--null-n"); if (length(i) == 1L) B <- as.integer(argv[i + 1L]) + i <- which(argv == "--block-len"); if (length(i) == 1L) block_len <- as.integer(argv[i + 1L]) + i <- which(argv == "--seed"); if (length(i) == 1L) seed <- as.integer(argv[i + 1L]) + parallel <- any(argv == "--parallel") + cfg <- default_audit_config(seed = seed, null_n = B, fast = is_fast(argv)) + if (isTRUE(cfg$fast)) B <- 20L + if (parallel) configure_parallel(TRUE, "auto") + setup_rng(cfg$seed) + + root <- audit_repo_root() + lines <- clean_lines_source(read_nist_csv(file.path(root, "data/Fe_lines.csv")), + "Fe", 2L, "wavenumber")$lines + k_grid <- audit_k_grid(cfg) + res <- run_scan_analysis(lines, cfg$bins_primary, k_grid, cfg$degree, cfg$baseline_sigma) + + models <- model_metadata()$null_model + model_seeds <- sample.int(.Machine$integer.max, length(models)) + runs <- lapply(seq_along(models), function(j) + run_one_model(models[j], res, k_grid, B, model_seeds[j], block_len, parallel)) + + tab <- do.call(rbind, lapply(runs, `[[`, "summary")) + tab <- merge(tab, model_metadata(), by = "null_model", sort = FALSE) + tab <- tab[match(models, tab$null_model), ] + tab$worst_case <- tab$scan_global_p == max(tab$scan_global_p, na.rm = TRUE) + maxima <- do.call(rbind, lapply(runs, `[[`, "maxima")) + + write_table(tab, file.path(root, "tables_r/statistical_audit/alternative_null_results.csv")) + write_table(maxima, file.path(root, "tables_r/statistical_audit/alternative_null_maxima.csv")) + save_fig(plot_alternative_nulls(maxima, tab), + file.path(root, "figures_r/statistical_audit/fig17_alternative_nulls.png"), + width = 10, height = 6.5) + + worst <- tab[which.max(tab$scan_global_p), ] + cat(sprintf("[alternative_nulls] %d models x %d replicates; worst-case p=%.4g (%s)\n", + nrow(tab), B, worst$scan_global_p, worst$null_model)) + invisible(tab) +} + +.invoked_file <- sub("^--file=", "", commandArgs(FALSE)[grep("^--file=", commandArgs(FALSE))]) +if (length(.invoked_file) > 0L && grepl("run_alternative_nulls\\.R$", .invoked_file)) { + main() +} diff --git a/R/run_final_adversarial_audit.R b/R/run_final_adversarial_audit.R new file mode 100644 index 0000000..aa61931 --- /dev/null +++ b/R/run_final_adversarial_audit.R @@ -0,0 +1,231 @@ +#!/usr/bin/env Rscript +# run_final_adversarial_audit.R +# --------------------------------------------------------------------------- +# Strict one-command final adversarial statistical audit. +# +# This entry point has NO silent Monte Carlo caps and NO partial-success mode. +# Any module failure aborts immediately so stale tables from a prior run cannot +# be summarized as current final evidence. +# +# Default final budgets: +# bootstrap_n 5000 +# scan/null_n 5000 +# source_bootstrap_n 2000 +# sensitivity_null_n 1000 +# holdout_null_n 5000 +# family_n 5000 +# calibration_n 10000 +# injection_n 2000 +# injection_null_n 5000 +# alternative_null_n 5000 per declared null model +# +# Example: +# Rscript R/run_final_adversarial_audit.R --parallel true --render-report true +# --------------------------------------------------------------------------- + +# Locate R/ robustly both when invoked directly and when sys.source()d from the +# test suite. commandArgs() describes the top-level process, not necessarily the +# currently sourced file, so it cannot be the sole locator. +locate_final_r_dir <- function() { + all_args <- commandArgs(FALSE) + file_args <- sub("^--file=", "", all_args[grep("^--file=", all_args)]) + invoked_dir <- if (length(file_args)) dirname(file_args[1]) else NA_character_ + helper_dir <- if (exists("AUDIT_PATH", inherits = TRUE)) dirname(get("AUDIT_PATH", inherits = TRUE)) else NA_character_ + candidates <- unique(c(invoked_dir, helper_dir, "R", file.path(getwd(), "R"), ".")) + candidates <- candidates[!is.na(candidates) & nzchar(candidates)] + for (d in candidates) { + if (file.exists(file.path(d, "audit_utils.R")) && + file.exists(file.path(d, "render_statistical_audit.R"))) { + return(normalizePath(d, mustWork = TRUE)) + } + } + stop("could not locate repository R/ directory for final adversarial audit") +} + +R_DIR <- locate_final_r_dir() +if (!exists("audit_repo_root", mode = "function")) source(file.path(R_DIR, "audit_utils.R")) + +parse_args <- function(argv) { + d <- list( + `bootstrap-n` = 5000L, + `null-n` = 5000L, + `source-bootstrap-n` = 2000L, + `sensitivity-null-n` = 1000L, + `holdout-null-n` = 5000L, + `family-n` = 5000L, + `calibration-n` = 10000L, + `injection-n` = 2000L, + `injection-null-n` = 5000L, + `alternative-null-n` = 5000L, + `block-len` = 8L, + parallel = FALSE, + `render-report` = TRUE + ) + i <- 1L + while (i <= length(argv)) { + key <- sub("^--", "", argv[i]) + if (!key %in% names(d)) { i <- i + 1L; next } + if (key %in% c("parallel", "render-report")) { + nxt <- if (i + 1L <= length(argv)) tolower(argv[i + 1L]) else NA_character_ + if (!is.na(nxt) && nxt %in% c("true", "false")) { + d[[key]] <- nxt == "true"; i <- i + 2L + } else { + d[[key]] <- TRUE; i <- i + 1L + } + } else { + if (i + 1L > length(argv)) stop("missing value for --", key) + d[[key]] <- as.integer(argv[i + 1L]); i <- i + 2L + } + } + nums <- setdiff(names(d), c("parallel", "render-report")) + bad <- nums[!vapply(d[nums], function(x) is.finite(x) && x > 0, logical(1))] + if (length(bad)) stop("all final audit budgets must be positive: ", paste(bad, collapse = ", ")) + d +} + +run_checked <- function(name, argv = character(0)) { + cat(sprintf("\n=== FINAL ADVERSARIAL [%s] ===\n", name)) + t0 <- Sys.time() + e <- new.env(parent = globalenv()) + ok <- FALSE; err <- "" + tryCatch({ + sys.source(file.path(R_DIR, name), envir = e) + if (!exists("main", envir = e, inherits = FALSE)) stop("module has no main(): ", name) + e$main(argv) + ok <- TRUE + }, error = function(x) { err <<- conditionMessage(x) }) + secs <- as.numeric(difftime(Sys.time(), t0, units = "secs")) + row <- data.frame(module = name, ok = ok, seconds = secs, error = err, + stringsAsFactors = FALSE) + if (!ok) stop(sprintf("final adversarial audit aborted at %s: %s", name, err)) + row +} + +require_outputs <- function(root, names) { + td <- file.path(root, "tables_r/statistical_audit") + miss <- names[!file.exists(file.path(td, names))] + if (length(miss)) stop("final audit required outputs missing: ", paste(miss, collapse = ", ")) + invisible(TRUE) +} + +main <- function(argv = commandArgs(TRUE)) { + args <- parse_args(argv) + root <- audit_repo_root() + td <- file.path(root, "tables_r/statistical_audit") + fd <- file.path(root, "figures_r/statistical_audit") + setup_rng(20260517L) + pflag <- if (isTRUE(args$parallel)) "--parallel" else character(0) + status <- list() + run <- function(name, a = character(0)) { + status[[length(status) + 1L]] <<- run_checked(name, a) + } + + run("build_analysis_registry.R") + run("build_dataset_flow.R") + + run("bootstrap_peak_uncertainty.R", + c("--bootstrap-n", as.character(args[["bootstrap-n"]]), pflag)) + run("peak_stability.R") + + run("build_effect_size_table.R", + c("--null-n", as.character(args[["null-n"]]), pflag)) + run("run_bin_grid.R", + c("--null-n", as.character(args[["null-n"]]), pflag)) + run("run_model_sensitivity.R", + c("--null-n", as.character(args[["sensitivity-null-n"]]), pflag)) + + run("model_comparison.R") + run("run_observed_ritz_replication.R", + c("--null-n", as.character(args[["null-n"]]), + "--bootstrap-n", as.character(args[["source-bootstrap-n"]]), pflag)) + run("run_holdout_replication.R", + c("--null-n", as.character(args[["holdout-null-n"]]), pflag)) + + # Full family now includes the complete declared sigma/degree/bin multiverse. + run("global_multiple_testing.R", + c("--family-n", as.character(args[["family-n"]]), pflag)) + + run("calibrate_false_positive_rate.R", + c("--calibration-n", as.character(args[["calibration-n"]]), + "--null-n", as.character(args[["null-n"]]), pflag)) + run("run_injection_recovery.R", + c("--injection-n", as.character(args[["injection-n"]]), + "--null-n", as.character(args[["injection-null-n"]]), pflag)) + + run("run_alternative_nulls.R", + c("--null-n", as.character(args[["alternative-null-n"]]), + "--block-len", as.character(args[["block-len"]]), + "--seed", "20260517", pflag)) + + # Reuse deterministic parity/dashboard functions, but do not invoke the old + # general orchestrator because it contains development-oriented caps. + renderer <- new.env(parent = globalenv()) + sys.source(file.path(R_DIR, "render_statistical_audit.R"), envir = renderer) + renderer$build_parity(file.path(root, "outputs"), file.path(root, "outputs_r"), td, fd) + tryCatch(renderer$build_dashboard(fd, td), error = function(e) message("[dashboard] ", conditionMessage(e))) + # Preserve the historical claim matrix too; the new one is stricter. + tryCatch(renderer$build_claim_matrix(td), error = function(e) message("[legacy claim matrix] ", conditionMessage(e))) + + run("build_final_adversarial_summary.R") + + required <- c( + "analysis_registry.csv", "dataset_flow_by_species.csv", "peak_confidence_intervals.csv", + "significance_results.csv", "multiple_testing.csv", "null_calibration.csv", + "holdout_results.csv", "model_comparison.csv", "injection_recovery.csv", + "specification_results.csv", "alternative_null_results.csv", + "calibrated_global_evidence.csv", "final_adversarial_claim_matrix.csv", "failed_claims.csv" + ) + require_outputs(root, required) + + # Write CURRENT provenance before rendering so the report cannot embed an old + # manifest from a previous run. + stat <- do.call(rbind, status) + expected_report <- file.path(root, "reports/rendered/nist_final_adversarial_audit.html") + manifest <- data.frame( + timestamp = audit_timestamp(), + git_commit = git_commit(), + canonical_seed = 20260517L, + bootstrap_n = args[["bootstrap-n"]], + null_n = args[["null-n"]], + source_bootstrap_n = args[["source-bootstrap-n"]], + sensitivity_null_n = args[["sensitivity-null-n"]], + holdout_null_n = args[["holdout-null-n"]], + family_n = args[["family-n"]], + calibration_n = args[["calibration-n"]], + injection_n = args[["injection-n"]], + injection_null_n = args[["injection-null-n"]], + alternative_null_n = args[["alternative-null-n"]], + alternative_block_len = args[["block-len"]], + parallel = args$parallel, + all_modules_ok = all(stat$ok), + report_requested = args[["render-report"]], + report_expected_path = expected_report, + stringsAsFactors = FALSE + ) + write_table(manifest, file.path(td, "final_adversarial_run_manifest.csv")) + write_table(stat, file.path(td, "final_adversarial_module_status.csv")) + + report_path <- NA_character_ + if (isTRUE(args[["render-report"]])) { + report_path <- renderer$render_report(file.path(root, "reports/nist_final_adversarial_audit.qmd"), + root, FALSE) + if (is.na(report_path)) message("[report] final tables are valid, but HTML was not rendered") + } + + ev <- utils::read.csv(file.path(td, "calibrated_global_evidence.csv"), stringsAsFactors = FALSE) + fc <- utils::read.csv(file.path(td, "failed_claims.csv"), stringsAsFactors = FALSE) + cat("\n=====================================================\n") + cat("NIST FINAL ADVERSARIAL STATISTICAL AUDIT\n") + cat("=====================================================\n") + cat(sprintf("largest reported global p : %.6g\n", ev$largest_reported_global_p[1])) + cat(" (diagnostic maximum across reported global/multiplicity/alternative-null p-values; not a combined p)\n") + cat(sprintf("non-passing claims : %d\n", nrow(fc))) + cat(sprintf("modules : %d (%d ok)\n", nrow(stat), sum(stat$ok))) + cat(sprintf("report : %s\n", + if (is.na(report_path)) "not rendered; use CSV outputs" else report_path)) + cat("\nNo failed or mixed claim is suppressed from failed_claims.csv.\n") + invisible(list(status = stat, manifest = manifest, evidence = ev, failed_claims = fc)) +} + +.invoked_file <- sub("^--file=", "", commandArgs(FALSE)[grep("^--file=", commandArgs(FALSE))]) +if (length(.invoked_file) > 0L && grepl("run_final_adversarial_audit\\.R$", .invoked_file)) main() diff --git a/R/run_holdout_replication.R b/R/run_holdout_replication.R index 3d5bf22..a6f1017 100644 --- a/R/run_holdout_replication.R +++ b/R/run_holdout_replication.R @@ -1,21 +1,15 @@ #!/usr/bin/env Rscript # run_holdout_replication.R # --------------------------------------------------------------------------- -# Blocked holdout validation in log-wavenumber space (NOT random per-line -# splits alone). Designs: -# 1. lower-range discovery / upper-range validation -# 2. upper-range discovery / lower-range validation -# 3. alternating contiguous blocks -# 4. repeated blocked K-fold +# Blocked holdout validation in log-wavenumber space. k is learned on TRAINING +# bins and LOCKED before the test block is evaluated. Confirmatory test blocks +# are never rescanned to choose k; exploratory rescans remain labelled. # -# CRITICAL confirmatory rule: k is estimated on TRAINING bins and LOCKED before -# the test block is touched. The test block is NOT rescanned in confirmatory -# mode. An exploratory rescan of the test block is reported separately and -# clearly labelled. -# -# Writes: -# tables_r/statistical_audit/holdout_results.csv -# figures_r/statistical_audit/holdout_validation.png (+ fig09) +# Two fixed-k null calibrations are reported: +# fixed_k_test_p : historical calibration with test baseline fixed +# fixed_k_test_p_refit : adversarial calibration that re-estimates the +# smoothing baseline inside every null replicate +# The conservative p-value is max(the two) and drives final holdout claims. # --------------------------------------------------------------------------- .this <- sub("^--file=", "", commandArgs(FALSE)[grep("^--file=", commandArgs(FALSE))][1]) @@ -23,7 +17,6 @@ if (!exists("emp_p", mode = "function")) source(file.path(if (length(.this) == 0 poisson_loglik <- function(y, mu) sum(stats::dpois(y, pmax(mu, EPS), log = TRUE)) -# fit_block(): fit M0 and M1 (at a given k) on a set of bins; return metrics. fit_block <- function(ell, y, baseline, degree, k) { X0 <- design_poly(ell, degree) f0 <- fit_poisson_loglinear(y, baseline, X0) @@ -36,38 +29,49 @@ fit_block <- function(ell, y, baseline, degree, k) { ll0 = poisson_loglik(y, f0$mu), ll1 = poisson_loglik(y, f1$mu)) } -# fixed_k_test_p(): parametric null p-value for the test block deltaD AT the -# locked k (no rescan). y0 ~ Poisson(mu0_test); deltaD at fixed k. -fixed_k_test_p <- function(ell, y, baseline, degree, k, B) { +# Parametric fixed-k test. The observed statistic always uses the observed +# test-block baseline. Under refit_baseline=TRUE, every synthetic null dataset +# gets its own freshly estimated smoothing baseline before M0/M1 are fitted. +fixed_k_test_p <- function(ell, y, baseline, degree, k, B, + sigma = NULL, refit_baseline = FALSE) { + if (isTRUE(refit_baseline) && (is.null(sigma) || !is.finite(sigma))) + stop("sigma is required when refit_baseline=TRUE") obs <- fit_block(ell, y, baseline, degree, k)$deltaD X0 <- design_poly(ell, degree) mu0 <- fit_poisson_loglinear(y, baseline, X0)$mu seeds <- sample.int(.Machine$integer.max, B) vals <- vapply(seq_len(B), function(i) { set.seed(seeds[i]); y0 <- rpois(length(mu0), mu0) - fit_block(ell, y0, baseline, degree, k)$deltaD + base0 <- if (isTRUE(refit_baseline)) + pmax(gaussian_filter_nearest(y0, sigma), EPS) else baseline + fit_block(ell, y0, base0, degree, k)$deltaD }, numeric(1)) - list(p = emp_p(sum(vals >= obs), B), obs = obs, B = B) + list(p = emp_p(sum(vals >= obs), B), obs = obs, B = B, + refit_baseline = isTRUE(refit_baseline)) } -# circular phase difference in [0, pi] -phase_diff <- function(a, b) { - d <- abs(((a - b + pi) %% (2 * pi)) - pi); d -} +phase_diff <- function(a, b) abs(((a - b + pi) %% (2 * pi)) - pi) -# one_holdout(): train on `tr` bins, lock k, validate on `te` bins. one_holdout <- function(design, ell, y, sigma, degree, k_grid, tr, te, B) { base_tr <- pmax(gaussian_filter_nearest(y[tr], sigma), EPS) base_te <- pmax(gaussian_filter_nearest(y[te], sigma), EPS) - # lock k on training + + # Training-only frequency selection. sk <- scan_k(ell[tr], y[tr], base_tr, k_grid, degree) k_lock <- sk$best$k_best ftr <- fit_block(ell[tr], y[tr], base_tr, degree, k_lock) - # confirmatory: locked-k on test + + # Confirmatory test at locked k only. fte <- fit_block(ell[te], y[te], base_te, degree, k_lock) - ptest <- fixed_k_test_p(ell[te], y[te], base_te, degree, k_lock, B) - # exploratory: rescan test + p_fixed <- fixed_k_test_p(ell[te], y[te], base_te, degree, k_lock, B, + sigma = sigma, refit_baseline = FALSE) + p_refit <- fixed_k_test_p(ell[te], y[te], base_te, degree, k_lock, B, + sigma = sigma, refit_baseline = TRUE) + p_cons <- max(p_fixed$p, p_refit$p) + + # Exploratory only: test-block rescan cannot rescue confirmatory failure. sk_te <- scan_k(ell[te], y[te], base_te, k_grid, degree) + data.frame( design = design, n_train = length(tr), n_test = length(te), k_lock = k_lock, @@ -78,7 +82,10 @@ one_holdout <- function(design, ell, y, sigma, degree, k_grid, tr, te, B) { phase_diff = phase_diff(fte$phase, ftr$phase), predictive_loglik_gain = fte$ll1 - fte$ll0, direction_consistent = (fte$deltaD > 0), - fixed_k_test_p = ptest$p, test_B = B, + fixed_k_test_p = p_fixed$p, + fixed_k_test_p_refit = p_refit$p, + conservative_fixed_k_test_p = p_cons, + test_B = B, exploratory_rescan_k = sk_te$best$k_best, exploratory_rescan_deltaD = sk_te$best$deltaD, stringsAsFactors = FALSE) @@ -95,17 +102,15 @@ run_holdout <- function(cfg, bins = cfg$bins_primary, B = 500L, kfold = 5L) { sigma <- cfg$baseline_sigma; degree <- cfg$degree rows <- list() - # 1. lower train / upper test rows[[1]] <- one_holdout("lower_train_upper_test", ell, y, sigma, degree, k_grid, seq_len(half), seq.int(half + 1L, n), B) - # 2. upper train / lower test rows[[2]] <- one_holdout("upper_train_lower_test", ell, y, sigma, degree, k_grid, seq.int(half + 1L, n), seq_len(half), B) - # 3. alternating contiguous blocks (4 blocks: train on 1,3 ; test on 2,4) + blk <- cut(seq_len(n), breaks = 4, labels = FALSE) tr3 <- which(blk %in% c(1, 3)); te3 <- which(blk %in% c(2, 4)) rows[[3]] <- one_holdout("alternating_blocks", ell, y, sigma, degree, k_grid, tr3, te3, B) - # 4. repeated blocked K-fold + fold <- cut(seq_len(n), breaks = kfold, labels = FALSE) for (f in seq_len(kfold)) { te <- which(fold == f); tr <- which(fold != f) @@ -120,16 +125,16 @@ plot_holdout <- function(ho) { ggplot2::ggplot(ho, ggplot2::aes(test_deltaD_lockedk, design, fill = direction_consistent)) + ggplot2::geom_col(width = 0.7) + - ggplot2::geom_text(ggplot2::aes(label = sprintf("p=%.3f", fixed_k_test_p)), + ggplot2::geom_text(ggplot2::aes(label = sprintf("p_cons=%.3f", conservative_fixed_k_test_p)), hjust = -0.1, size = 3) + ggplot2::scale_fill_manual(values = c(`TRUE` = species_colour("Fe"), `FALSE` = "#999999"), name = "direction consistent") + - ggplot2::expand_limits(x = max(ho$test_deltaD_lockedk, na.rm = TRUE) * 1.25) + + ggplot2::expand_limits(x = max(ho$test_deltaD_lockedk, na.rm = TRUE) * 1.28) + ggplot2::labs( title = "Blocked holdout validation (k locked from training)", - subtitle = "Test-block deviance improvement at the locked k; fixed-k empirical p shown.", + subtitle = "Conservative fixed-k p = max(fixed-baseline, baseline-refit null calibration).", x = "test-block deltaD at locked k", y = NULL, - caption = "Confirmatory: no test-block rescan. Exploratory rescans reported in the table only.") + + caption = "No confirmatory test-block rescan. Exploratory rescans are table diagnostics only.") + theme_audit() } @@ -143,12 +148,11 @@ main <- function(argv = commandArgs(TRUE)) { write_table(ho, file.path(root, "tables_r/statistical_audit/holdout_results.csv")) save_fig(plot_holdout(ho), file.path(root, "figures_r/statistical_audit/holdout_validation.png")) save_fig(plot_holdout(ho), file.path(root, "figures_r/statistical_audit/fig09_holdout_validation.png")) - cat(sprintf("[holdout] %d designs; %d with consistent direction; median fixed-k p=%.3f\n", - nrow(ho), sum(ho$direction_consistent), stats::median(ho$fixed_k_test_p))) + cat(sprintf("[holdout] %d designs; %d direction-consistent; median conservative p=%.3f\n", + nrow(ho), sum(ho$direction_consistent), + stats::median(ho$conservative_fixed_k_test_p))) invisible(ho) } .invoked_file <- sub("^--file=", "", commandArgs(FALSE)[grep("^--file=", commandArgs(FALSE))]) -if (length(.invoked_file) > 0L && grepl("run_holdout_replication\\.R$", .invoked_file)) { - main() -} +if (length(.invoked_file) > 0L && grepl("run_holdout_replication\\.R$", .invoked_file)) main() diff --git a/R/run_resolution_mode_diagnostics.R b/R/run_resolution_mode_diagnostics.R new file mode 100644 index 0000000..61efb28 --- /dev/null +++ b/R/run_resolution_mode_diagnostics.R @@ -0,0 +1,316 @@ +#!/usr/bin/env Rscript +# run_resolution_mode_diagnostics.R +# --------------------------------------------------------------------------- +# FOLLOW-UP diagnostic for the Fe II resolution-dependent mode transition. +# +# This module is deliberately separate from run_final_adversarial_audit.R. +# It MUST NOT overwrite, rescue, or reinterpret the frozen adversarial verdicts. +# In particular, the original 60:240 bin-grid stability result remains the +# primary robustness test and its pass/fail criterion is unchanged. +# +# Questions addressed here: +# 1. Does the observed k~9.6 -> k~31.3 transition persist when Gaussian +# smoothing is held approximately fixed in physical ell=log(wavenumber) +# width rather than fixed in bin units? +# 2. Are k=9.602325620315224 and the canonical Fe primary k simultaneously +# present as competing local modes even when only one is the global winner? +# 3. Does the previously frozen GWTC k=9.602325620315224 survive a direct +# fixed-frequency null test in Fe and neighbouring ion-II datasets? +# +# The GWTC frequency existed before this NIST follow-up, but noticing its match +# to the coarse-bin NIST branch is post-hoc. Results are therefore exploratory +# cross-domain follow-up evidence, not original preregistration. +# +# Writes: +# tables_r/statistical_audit/resolution_mode_dense.csv +# tables_r/statistical_audit/resolution_mode_top_peaks.csv +# tables_r/statistical_audit/external_gwtc_k_fixed_tests.csv +# figures_r/statistical_audit/resolution_mode_transition.png +# figures_r/statistical_audit/resolution_mode_strength.png +# --------------------------------------------------------------------------- + +.this <- sub("^--file=", "", commandArgs(FALSE)[grep("^--file=", commandArgs(FALSE))][1]) +if (!exists("emp_p", mode = "function")) { + source(file.path(if (length(.this) == 0L || is.na(.this)) "R" else dirname(.this), + "audit_utils.R")) +} + +GWTC_K_FROZEN <- 9.602325620315224 + +matched_sigma <- function(bins, reference_bins = 160L, reference_sigma = 6.0) { + as.numeric(reference_sigma) * as.numeric(bins) / as.numeric(reference_bins) +} + +local_peak_indices <- function(deltaD) { + x <- as.numeric(deltaD) + n <- length(x) + if (n == 0L) return(integer()) + if (n == 1L) return(1L) + x[!is.finite(x)] <- -Inf + idx <- integer() + if (x[1] > x[2]) idx <- c(idx, 1L) + if (n > 2L) { + mid <- which(x[2:(n - 1L)] >= x[1:(n - 2L)] & + x[2:(n - 1L)] > x[3:n]) + 1L + idx <- c(idx, mid) + } + if (x[n] >= x[n - 1L]) idx <- c(idx, n) + unique(idx[is.finite(x[idx])]) +} + +extract_top_peaks <- function(scan, n_top = 5L) { + idx <- local_peak_indices(scan$deltaD) + if (length(idx) == 0L) idx <- which.max(scan$deltaD) + out <- data.frame( + k = scan$k[idx], + deltaD = scan$deltaD[idx], + stringsAsFactors = FALSE + ) + out <- out[order(out$deltaD, decreasing = TRUE), , drop = FALSE] + out <- head(out, max(1L, as.integer(n_top))) + out$peak_rank <- seq_len(nrow(out)) + rownames(out) <- NULL + out +} + +classify_peak_branch <- function(k, k_external = GWTC_K_FROZEN, k_primary, + tol_relative = 0.02) { + if (!is.finite(k)) return("other") + if (abs(k - k_external) / k_external <= tol_relative) return("external_gwtc_9p602") + if (abs(k - k_primary) / k_primary <= tol_relative) return("primary_fe") + "other" +} + +scan_value_at_k <- function(res, k_target) { + idx <- nearest_k_index(res$scan$k, k_target) + c(k_grid = res$scan$k[idx], deltaD = res$scan$deltaD[idx]) +} + +resolution_scan_one <- function(lines, bins, sigma_bins, smoothing_mode, cfg, + k_primary, k_external = GWTC_K_FROZEN, + n_top = 5L) { + k_grid <- audit_k_grid(cfg) + res <- run_scan_analysis(lines, bins, k_grid, cfg$degree, sigma_bins) + pext <- scan_value_at_k(res, k_external) + ppri <- scan_value_at_k(res, k_primary) + bin_width <- stats::median(res$binned$edge_hi - res$binned$edge_lo) + best <- res$best$deltaD + summary <- data.frame( + smoothing_mode = smoothing_mode, + bins = as.integer(bins), + sigma_bins = as.numeric(sigma_bins), + sigma_ell_approx = as.numeric(sigma_bins) * bin_width, + n_lines = res$n_lines, + k_best = res$best$k_best, + deltaD_best = best, + best_branch = classify_peak_branch(res$best$k_best, k_external, k_primary, + cfg$tol_primary), + k_external_target = k_external, + k_external_grid = unname(pext["k_grid"]), + deltaD_external = unname(pext["deltaD"]), + external_to_best = unname(pext["deltaD"]) / best, + k_primary_target = k_primary, + k_primary_grid = unname(ppri["k_grid"]), + deltaD_primary = unname(ppri["deltaD"]), + primary_to_best = unname(ppri["deltaD"]) / best, + stringsAsFactors = FALSE + ) + peaks <- extract_top_peaks(res$scan, n_top) + peaks$smoothing_mode <- smoothing_mode + peaks$bins <- as.integer(bins) + peaks$sigma_bins <- as.numeric(sigma_bins) + peaks$sigma_ell_approx <- as.numeric(sigma_bins) * bin_width + peaks$branch <- vapply(peaks$k, classify_peak_branch, character(1), + k_external = k_external, k_primary = k_primary, + tol_relative = cfg$tol_primary) + peaks <- peaks[, c("smoothing_mode", "bins", "sigma_bins", "sigma_ell_approx", + "peak_rank", "k", "deltaD", "branch")] + list(summary = summary, peaks = peaks) +} + +run_dense_resolution_map <- function(lines, cfg, k_primary, + bins_dense = seq(40L, 300L, by = 5L), + n_top = 5L, parallel = FALSE) { + jobs <- do.call(rbind, lapply(as.integer(bins_dense), function(b) { + data.frame( + bins = c(b, b), + smoothing_mode = c("fixed_sigma_bins", "matched_sigma_ell"), + sigma_bins = c(cfg$baseline_sigma, + matched_sigma(b, cfg$bins_primary, cfg$baseline_sigma)), + stringsAsFactors = FALSE + ) + })) + one <- function(i) { + g <- jobs[i, ] + resolution_scan_one(lines, g$bins, g$sigma_bins, g$smoothing_mode, + cfg, k_primary, GWTC_K_FROZEN, n_top) + } + ans <- audit_lapply(seq_len(nrow(jobs)), one, parallel = parallel, seed = cfg$seed) + list( + summary = do.call(rbind, lapply(ans, `[[`, "summary")), + peaks = do.call(rbind, lapply(ans, `[[`, "peaks")) + ) +} + +fixed_frequency_null <- function(res, k_target, B, seed, + refit_baseline = FALSE, parallel = FALSE) { + if (B <= 0L) return(list(observed = NA_real_, p = NA_real_, tail = NA_integer_)) + observed <- scan_k(res$ell, res$y, res$baseline, c(k_target), res$degree)$best$deltaD + setup_rng(seed) + seeds <- sample.int(.Machine$integer.max, B) + one <- function(i) { + set.seed(seeds[i]) + y0 <- stats::rpois(length(res$mu0), res$mu0) + baseline0 <- if (isTRUE(refit_baseline)) { + pmax(gaussian_filter_nearest(y0, res$baseline_sigma), EPS) + } else { + res$baseline + } + scan_k(res$ell, y0, baseline0, c(k_target), res$degree)$best$deltaD + } + sims <- unlist(audit_lapply(seq_len(B), one, parallel = parallel, seed = seed), + use.names = FALSE) + tail <- sum(sims >= observed) + list(observed = observed, p = emp_p(tail, B), tail = as.integer(tail)) +} + +fixed_external_target_one <- function(raw, species, source, cfg, B, + parallel = FALSE) { + lines <- clean_lines_source(raw, species, 2L, source)$lines + k_grid <- audit_k_grid(cfg) + res <- run_scan_analysis(lines, cfg$bins_primary, k_grid, cfg$degree, + cfg$baseline_sigma) + fixed <- fixed_frequency_null(res, GWTC_K_FROZEN, B, cfg$seed, + refit_baseline = FALSE, parallel = parallel) + refit <- fixed_frequency_null(res, GWTC_K_FROZEN, B, cfg$seed + 101L, + refit_baseline = TRUE, parallel = parallel) + data.frame( + species = species, + ion = 2L, + source = source, + bins = cfg$bins_primary, + sigma = cfg$baseline_sigma, + degree = cfg$degree, + n_lines = nrow(lines), + external_target = "GWTC_frozen_k", + k_fixed = GWTC_K_FROZEN, + k_selected_full_scan = res$best$k_best, + selected_deltaD = res$best$deltaD, + fixed_deltaD = fixed$observed, + fixed_to_selected_ratio = fixed$observed / res$best$deltaD, + fixed_baseline_p = fixed$p, + fixed_baseline_tail = fixed$tail, + refit_baseline_p = refit$p, + refit_baseline_tail = refit$tail, + conservative_fixed_p = max(fixed$p, refit$p, na.rm = TRUE), + null_n = as.integer(B), + posthoc_cross_domain_followup = TRUE, + stringsAsFactors = FALSE + ) +} + +run_external_target_tests <- function(root, cfg, B = 5000L, parallel = FALSE) { + species <- c("Fe", "Cr", "Mn", "Co", "Ni", "Ti") + rows <- list() + j <- 1L + for (sp in species) { + path <- file.path(root, "data", sprintf("%s_lines.csv", sp)) + raw <- read_nist_csv(path) + sources <- if (sp == "Fe") c("wavenumber", "observed", "ritz") else "wavenumber" + for (src in sources) { + rows[[j]] <- fixed_external_target_one(raw, sp, src, cfg, B, parallel) + cat(sprintf("[external-k] %s/%s k=%.10f p_fixed=%.4g p_refit=%.4g p_cons=%.4g\n", + sp, src, GWTC_K_FROZEN, + rows[[j]]$fixed_baseline_p, + rows[[j]]$refit_baseline_p, + rows[[j]]$conservative_fixed_p)) + j <- j + 1L + } + } + do.call(rbind, rows) +} + +plot_resolution_transition <- function(dense, k_primary) { + ggplot2::ggplot(dense, ggplot2::aes(bins, k_best, linetype = smoothing_mode)) + + ggplot2::geom_line(linewidth = 0.7) + + ggplot2::geom_point(size = 1.5) + + ggplot2::geom_hline(yintercept = GWTC_K_FROZEN, linetype = 3) + + ggplot2::geom_hline(yintercept = k_primary, linetype = 2) + + ggplot2::labs( + title = "Fe II resolution-dependent selected mode", + subtitle = "Original fixed-sigma-bin pipeline versus approximately fixed ell-space smoothing width.", + x = "bin count", y = "selected k (rad / ln cm^-1)", + linetype = "smoothing") + + theme_audit() +} + +plot_mode_strength <- function(dense) { + long <- rbind( + data.frame(smoothing_mode = dense$smoothing_mode, bins = dense$bins, + target = "GWTC k=9.6023256", ratio = dense$external_to_best), + data.frame(smoothing_mode = dense$smoothing_mode, bins = dense$bins, + target = "Fe primary k", ratio = dense$primary_to_best) + ) + ggplot2::ggplot(long, ggplot2::aes(bins, ratio, linetype = target)) + + ggplot2::geom_line(linewidth = 0.7) + + ggplot2::facet_wrap(~smoothing_mode) + + ggplot2::geom_hline(yintercept = 1, linetype = 3) + + ggplot2::labs( + title = "Strength of the two tracked modes relative to the winning peak", + subtitle = "A ratio near 1 means the tracked mode is competitive with the global winner.", + x = "bin count", y = "deltaD(target) / deltaD(best)", linetype = "target") + + theme_audit() +} + +arg_int <- function(argv, flag, default) { + i <- which(argv == flag) + if (length(i) == 0L) return(as.integer(default)) + as.integer(argv[i[1] + 1L]) +} + +main <- function(argv = commandArgs(TRUE)) { + parallel <- any(argv == "--parallel") + fast <- is_fast(argv) + if (parallel) configure_parallel(TRUE, "auto") + cfg <- default_audit_config(fast = fast) + root <- audit_repo_root() + + min_bin <- arg_int(argv, "--min-bin", if (fast) 60L else 40L) + max_bin <- arg_int(argv, "--max-bin", if (fast) 220L else 300L) + step_bin <- arg_int(argv, "--step-bin", if (fast) 20L else 5L) + fixed_null_n <- arg_int(argv, "--fixed-null-n", if (fast) 25L else 5000L) + bins_dense <- seq(min_bin, max_bin, by = step_bin) + + raw_fe <- read_nist_csv(file.path(root, "data/Fe_lines.csv")) + lines_fe <- clean_lines_source(raw_fe, "Fe", 2L, "wavenumber")$lines + k_primary <- run_scan_analysis(lines_fe, cfg$bins_primary, audit_k_grid(cfg), + cfg$degree, cfg$baseline_sigma)$best$k_best + + cat(sprintf("[resolution] Fe primary k=%.10f; frozen external GWTC k=%.10f\n", + k_primary, GWTC_K_FROZEN)) + dense <- run_dense_resolution_map(lines_fe, cfg, k_primary, bins_dense, + n_top = 5L, parallel = parallel) + write_table(dense$summary, + file.path(root, "tables_r/statistical_audit/resolution_mode_dense.csv")) + write_table(dense$peaks, + file.path(root, "tables_r/statistical_audit/resolution_mode_top_peaks.csv")) + save_fig(plot_resolution_transition(dense$summary, k_primary), + file.path(root, "figures_r/statistical_audit/resolution_mode_transition.png")) + save_fig(plot_mode_strength(dense$summary), + file.path(root, "figures_r/statistical_audit/resolution_mode_strength.png")) + + ext <- run_external_target_tests(root, cfg, fixed_null_n, parallel) + write_table(ext, + file.path(root, "tables_r/statistical_audit/external_gwtc_k_fixed_tests.csv")) + + cat(sprintf("[resolution] wrote %d dense resolution rows and %d top-peak rows\n", + nrow(dense$summary), nrow(dense$peaks))) + cat(sprintf("[resolution] external fixed-k tests=%d, B=%d; results are exploratory follow-up\n", + nrow(ext), fixed_null_n)) + invisible(list(dense = dense, external = ext, k_primary = k_primary)) +} + +.invoked_file <- sub("^--file=", "", commandArgs(FALSE)[grep("^--file=", commandArgs(FALSE))]) +if (length(.invoked_file) > 0L && grepl("run_resolution_mode_diagnostics\\.R$", .invoked_file)) { + main() +} diff --git a/RESOLUTION_MODE_DIAGNOSTICS.md b/RESOLUTION_MODE_DIAGNOSTICS.md new file mode 100644 index 0000000..d77d84e --- /dev/null +++ b/RESOLUTION_MODE_DIAGNOSTICS.md @@ -0,0 +1,103 @@ +# Resolution-mode follow-up diagnostics + +This follow-up exists because the frozen adversarial bin-grid test can reveal a resolution-dependent winner without explaining why it changes. It is intentionally separate from the final audit and cannot overwrite or rescue a failed robustness claim. + +The original `R/run_bin_grid.R` result remains authoritative for the declared bin-stability verdict. If fewer than 80% of the predefined bin counts select the Fe II reference region, the final claim matrix reports that robustness claim as `fail` even if a later diagnostic explains the transition. + +## Why this follow-up is needed + +The canonical pipeline bins `ell = log(wavenumber)` and then applies a Gaussian baseline with `baseline_sigma = 6` measured in **bins**. Changing the number of histogram bins therefore changes both resolution and the approximate smoothing width in physical `ell` space. + +The follow-up separates those effects by comparing: + +- `fixed_sigma_bins`: the original convention, sigma = 6 at every bin count; +- `matched_sigma_ell`: sigma scaled as + +```text +sigma_bins(B) = 6 * B / 160 +``` + +so the approximate smoothing width in `ell` is held fixed relative to the 160-bin primary analysis. + +For reference, this gives sigma values 2.25, 3.00, 3.75, 4.50, 5.25, 6.00, 6.75, 7.50, 8.25, and 9.00 for 60 through 240 bins in steps of 20. + +## Two tracked frequencies + +The diagnostic records the full winning mode and also the strength of two explicit targets at every resolution: + +1. the canonical Fe II primary frequency selected by the 160-bin analysis; +2. `k = 9.602325620315224`, frozen previously in the independent GWTC repository. + +The GWTC value is taken from `rickyjreyes/GWTC/tables/gwtc_frozen_mode.json` and the corresponding frozen holdout result. The GWTC frequency predates this NIST resolution follow-up. However, noticing that the coarse-bin NIST branch lies near that value is **post-hoc**, so the NIST/GWTC comparison is reported only as exploratory cross-domain follow-up evidence. + +The dense map stores the top five distinct local maxima rather than only the global winner. This distinguishes a true disappearance of one mode from a winner swap between two simultaneously strong modes. + +## Fixed external-frequency test + +The script also tests the exact frozen GWTC frequency without re-optimizing `k` in NIST. At 160 bins it evaluates: + +- Fe II wavenumber source; +- Fe II observed-wavelength-derived source; +- Fe II Ritz-wavelength-derived source; +- Cr II, Mn II, Co II, Ni II, and Ti II wavenumber sources. + +For each dataset it reports the fixed-frequency deviance improvement and two empirical Poisson-null calibrations: + +- observed smoothing baseline held fixed; +- smoothing baseline re-estimated in every null replicate. + +The reported conservative fixed-frequency p-value is the larger of those two p-values. Every empirical p-value uses the plus-one correction `(r + 1)/(B + 1)`. + +Because this test was motivated after seeing the NIST resolution behavior, it is **not** promoted to the original confirmatory family and is not used to repair the final audit verdict. + +## Run + +Run the final audit first and preserve its outputs. Then pull the latest branch and execute: + +```powershell +Rscript R/run_resolution_mode_diagnostics.R --parallel true --fixed-null-n 5000 +``` + +The dense observed-data resolution map defaults to 40 through 300 bins in steps of 5. The bounds can be changed explicitly: + +```powershell +Rscript R/run_resolution_mode_diagnostics.R --parallel true --min-bin 40 --max-bin 300 --step-bin 5 --fixed-null-n 5000 +``` + +A development-only smoke run is available with: + +```powershell +Rscript R/run_resolution_mode_diagnostics.R --fast +``` + +`--fast` is not final evidence. + +## Outputs + +All outputs are written alongside the statistical-audit products: + +- `tables_r/statistical_audit/resolution_mode_dense.csv` + - winner at each bin count under both smoothing conventions; + - approximate physical smoothing width; + - deviance improvement at the frozen GWTC target and Fe primary target; + - each target's strength relative to the winning peak. +- `tables_r/statistical_audit/resolution_mode_top_peaks.csv` + - top five distinct local maxima at every resolution; + - branch labels for the frozen GWTC target region, Fe primary region, and other modes. +- `tables_r/statistical_audit/external_gwtc_k_fixed_tests.csv` + - exact fixed-frequency cross-ion/source tests; + - fixed-baseline, refit-baseline, and conservative empirical p-values. +- `figures_r/statistical_audit/resolution_mode_transition.png` + - selected frequency versus bin count under both smoothing conventions. +- `figures_r/statistical_audit/resolution_mode_strength.png` + - strength of the two tracked modes relative to the global winner. + +## Interpretation rules + +The follow-up is designed to discriminate among several possibilities without protecting a preferred result: + +- If the `k ~ 31.3` branch returns at coarse bin counts after matching physical smoothing width, the original coarse-bin failure is likely entangled with the smoothing convention. +- If `k ~ 9.6` remains dominant under matched smoothing, it is a genuine resolution-robust competing mode within this pipeline. +- If both tracked modes remain strong while the winner switches, the data support a multimodal spectrum more strongly than a single universally dominant frequency. +- If the frozen GWTC target fails under baseline refitting or across neighboring ions/sources, the cross-domain coincidence weakens. +- If the frozen GWTC target survives conservative fixed-frequency tests, that is noteworthy follow-up evidence but still does not establish WCT, a universal law, or independent experimental confirmation. diff --git a/reports/nist_final_adversarial_audit.qmd b/reports/nist_final_adversarial_audit.qmd new file mode 100644 index 0000000..241d512 --- /dev/null +++ b/reports/nist_final_adversarial_audit.qmd @@ -0,0 +1,294 @@ +--- +title: "NIST Atomic-Line Log-Cosine Scan — Final Adversarial Statistical Audit" +subtitle: "Global calibration, frozen holdouts, alternative nulls, and preserved failures" +date: today +format: + html: + toc: true + toc-depth: 3 + code-fold: true + embed-resources: true + theme: cosmo +execute: + echo: false + warning: false + message: false +--- + +```{r setup} +#| include: false +knitr::opts_chunk$set(fig.align = "center", out.width = "92%") +here_root <- normalizePath(file.path(dirname(knitr::current_input(dir = TRUE)), ".."), + mustWork = FALSE) +if (is.na(here_root) || !dir.exists(here_root)) here_root <- normalizePath("..", mustWork = FALSE) +TAB <- file.path(here_root, "tables_r/statistical_audit") +FIG <- file.path(here_root, "figures_r/statistical_audit") +rd <- function(f) { + p <- file.path(TAB, f) + if (file.exists(p)) utils::read.csv(p, stringsAsFactors = FALSE, check.names = FALSE) else NULL +} +fig <- function(f) { + p <- file.path(FIG, f) + if (file.exists(p)) p else NULL +} +show_tab <- function(df, n = NULL) { + if (is.null(df)) return(cat("_table not available_")) + if (!is.null(n)) df <- utils::head(df, n) + knitr::kable(df) +} +``` + +# 1. Executive verdict + +This report is the **final adversarial statistical audit** of the NIST atomic-line +log-cosine scan. It is deliberately stricter than the development audit: +full-resolution Monte Carlo budgets are passed to every module without hidden +500/1,000-trial caps; the confirmatory frequency is frozen before each holdout; +the complete declared preprocessing multiverse enters multiplicity correction; +alternative null generators are tested; and every failed, mixed, inconclusive, +or not-established claim is retained in a separate machine-readable table. + +The central global-evidence table is shown below. `largest_reported_global_p` is +the largest among the primary scan-global p-value, family-max/adjusted +multiplicity p-values, and the worst declared alternative-null p-value. This is a +**conservative diagnostic summary only**, not a newly combined formal p-value. + +```{r global-evidence} +show_tab(rd("calibrated_global_evidence.csv")) +``` + +# 2. Claims — including failures + +No result is promoted merely because another analysis passes. In particular, +failed held-out predictive transfer, weak bin/specification stability, an +alternative-null failure, or a non-significant frozen holdout remains failed or +mixed even if the in-sample peak is strong. + +```{r final-claims} +show_tab(rd("final_adversarial_claim_matrix.csv")) +``` + +## 2.1 Non-passing claims retained explicitly + +`failed_claims.csv` is generated mechanically from the claim matrix and contains +every verdict other than `pass`. + +```{r failed-claims} +show_tab(rd("failed_claims.csv")) +``` + +# 3. Exact run provenance + +The manifest is written **before** this report is rendered, so the report cannot +silently embed a manifest from a previous run. + +```{r provenance} +show_tab(rd("final_adversarial_run_manifest.csv")) +show_tab(rd("final_adversarial_module_status.csv")) +show_tab(rd("analysis_registry.csv"), 30) +show_tab(rd("dataset_flow_by_species.csv")) +``` + +# 4. Primary peak and bootstrap uncertainty + +Bootstrap uncertainty measures sampling and peak-selection stability within the +same line list. It is not independent experimental replication. + +```{r bootstrap} +show_tab(rd("peak_estimates.csv")) +show_tab(rd("peak_confidence_intervals.csv")) +show_tab(rd("peak_stability.csv")) +``` + +```{r bootstrap-figure} +p <- fig("fig05_bootstrap_peak_distribution.png") +if (!is.null(p)) knitr::include_graphics(p) +``` + +# 5. Full-scan look-elsewhere correction + +Pointwise and scan-global evidence remain separate. The scan-global statistic is +the maximum deviance improvement over the complete declared frequency grid. +Empirical p-values use `(r + 1)/(B + 1)`; zero exceedances are reported at the +Monte Carlo resolution floor rather than as `p=0`. + +```{r global-scan} +show_tab(rd("significance_results.csv")) +``` + +# 6. Full declared-family multiplicity correction + +The corrected family contains the predefined Fe II bin grid, the complete +sigma/degree/bin preprocessing multiverse, observed/Ritz representations, and +neighbouring ion-II controls, with duplicate numerical specifications removed. +The audit reports: + +- Bonferroni and Holm FWER control, neither of which requires independence; +- Benjamini-Hochberg FDR as a standard reference; +- Benjamini-Yekutieli FDR for arbitrary dependence; +- an approximate family maximum-statistic calibration. + +The family-max simulation couples marginal per-analysis null draws rather than +constructing a fully joint synthetic line-list generator, so its dependence +structure remains an explicit limitation. Holm/Bonferroni provide the +arbitrary-dependence FWER cross-check. + +```{r multiplicity} +show_tab(rd("multiple_testing.csv")) +``` + +```{r multiplicity-figure} +p <- fig("fig12_multiplicity_comparison.png") +if (!is.null(p)) knitr::include_graphics(p) +``` + +# 7. Adversarial alternative null models + +The canonical fixed-baseline fitted-Poisson null is supplemented with four +stress tests: + +1. fitted Poisson with Gaussian smoothing re-estimated in every replicate; +2. conditional multinomial counts with total line count fixed; +3. negative-binomial overdispersion with a method-of-moments size estimate; +4. block-permuted Pearson-residual rate perturbations followed by Poisson + sampling and baseline re-estimation. + +The conservative result is the largest scan-global p-value among these declared +nulls. These alternatives are **post-signal adversarial stress tests**, not a +claim of historical preregistration or of a unique physical null. The +block-residual construction is specifically a misspecification diagnostic. + +```{r alt-nulls} +show_tab(rd("alternative_null_results.csv")) +``` + +```{r alt-nulls-figure} +p <- fig("fig17_alternative_nulls.png") +if (!is.null(p)) knitr::include_graphics(p) +``` + +# 8. Synthetic-null false-positive calibration + +Calibration asks whether the complete scan-global decision rule achieves its +nominal type-I error under the canonical fitted-Poisson null. The finite Monte +Carlo interval is reported. Robustness to other generators is assessed +separately in the alternative-null section rather than conflated with this +calibration exercise. + +```{r calibration} +show_tab(rd("null_calibration.csv")) +``` + +```{r calibration-figure} +p <- fig("fig13_null_calibration.png") +if (!is.null(p)) knitr::include_graphics(p) +``` + +# 9. Frozen blocked holdouts + +For every confirmatory holdout, `k` is selected only on the training bins and +then locked. The test block is evaluated at that fixed frequency. Exploratory +test-block rescans are reported separately and cannot rescue a failed locked-k +result. The final claim matrix uses a strict rule: all declared holdout designs +must have positive direction and fixed-k `p <= 0.05` to receive `pass`; partial +success is `mixed`. + +```{r holdout} +show_tab(rd("holdout_results.csv")) +``` + +```{r holdout-figure} +p <- fig("fig09_holdout_validation.png") +if (!is.null(p)) knitr::include_graphics(p) +``` + +# 10. Preprocessing and specification sensitivity + +The audit-declared multiverse varies baseline smoothing, polynomial degree, bin +count, and line-position representation. An 80% reference-region fraction is +used as a descriptive robustness threshold in the final claim matrix; this is an +audit criterion, not a universal physical threshold and not represented as +historical preregistration. + +```{r sensitivity} +show_tab(rd("bin_stability_summary.csv")) +show_tab(rd("bin_grid_results.csv")) +show_tab(rd("specification_results.csv")) +show_tab(rd("observed_ritz_replication.csv")) +``` + +```{r sensitivity-figures} +p <- fig("fig10_specification_curve.png"); if (!is.null(p)) knitr::include_graphics(p) +p <- fig("specification_heatmap.png"); if (!is.null(p)) knitr::include_graphics(p) +``` + +# 11. In-sample versus held-out model comparison + +AIC/BIC improvement after selecting a frequency on the same data does not by +itself demonstrate predictive transfer. Held-out log-likelihood is reported +separately. A negative `M1 - M0` held-out gain is retained as a failed claim even +if the in-sample improvement is very large. + +```{r model-comparison} +show_tab(rd("model_comparison.csv")) +``` + +```{r model-comparison-figure} +p <- fig("fig11_model_comparison.png") +if (!is.null(p)) knitr::include_graphics(p) +``` + +# 12. Injection and recovery + +Injection/recovery reports both detection of any significant peak and recovery +inside the correct frequency region. A detector that fires at the wrong +frequency has not recovered the injected structure. The zero-amplitude row is +also retained as a direct type-I diagnostic. + +```{r injection} +show_tab(rd("injection_recovery.csv")) +``` + +```{r injection-figure} +p <- fig("fig14_injection_recovery_power.png") +if (!is.null(p)) knitr::include_graphics(p) +``` + +# 13. Computational parity + +```{r parity} +show_tab(rd("python_r_parity.csv")) +``` + +```{r parity-figure} +p <- fig("fig15_python_r_parity.png") +if (!is.null(p)) knitr::include_graphics(p) +``` + +# 14. Interpretation boundaries + +This audit can support or weaken a claim of statistical structure in the +declared NIST analysis. It cannot by itself establish: + +- independent experimental confirmation; +- independence of observed/Ritz representations or blocked subsets drawn from + the same transition list; +- a unique physical null model; +- a WCT mechanism, causal explanation, or universal atomic law; +- NIST endorsement of the analysis. + +Computational parity, same-dataset bootstrap stability, frozen holdouts, +alternative-null robustness, and independent physical replication are therefore +treated as different evidentiary categories. + +# 15. Reproduction command + +The strict final runner fails closed and never substitutes development budgets. + +```{r reproduce} +cat("Rscript R/run_final_adversarial_audit.R --parallel true --render-report true\n\n") +manifest <- rd("final_adversarial_run_manifest.csv") +if (!is.null(manifest)) show_tab(manifest) +cat("\nSession information for this render:\n") +sessionInfo() +``` diff --git a/tables_r/statistical_audit/alternative_null_results.csv b/tables_r/statistical_audit/alternative_null_results.csv new file mode 100644 index 0000000..d0b24cb --- /dev/null +++ b/tables_r/statistical_audit/alternative_null_results.csv @@ -0,0 +1,6 @@ +"null_model","B","observed_deltaD","tail_count","scan_global_p","resolution_floor","exceedance_rate_ci_lo","exceedance_rate_ci_hi","null_median","null_q95","null_q99","null_max","nb_size","block_len","description","interpretation","worst_case" +"poisson_fixed_baseline",5000,315.727660603729,0,0.0001999600079984,0.0001999600079984,0,0.000737503801108105,8.73190201477571,14.6822768991152,18.0892965989488,23.6672128792788,NA,NA,"Fitted smooth Poisson null; observed smoothing baseline held fixed (canonical audit null).","reference",FALSE +"poisson_refit_baseline",5000,315.727660603729,0,0.0001999600079984,0.0001999600079984,0,0.000737503801108105,49.9931468107456,59.4571087228756,63.9473648369696,72.9850319621235,NA,NA,"Fitted smooth Poisson null; Gaussian smoothing baseline re-estimated in every replicate.","tests baseline-estimation uncertainty",FALSE +"conditional_multinomial",5000,315.727660603729,0,0.0001999600079984,0.0001999600079984,0,0.000737503801108105,49.7635764597087,59.1592288582691,63.2882828798395,70.6641970688671,NA,NA,"Conditional multinomial null with total line count fixed; smoothing baseline re-estimated.","conditions on observed total count",FALSE +"negative_binomial_refit",5000,315.727660603729,0,0.0001999600079984,0.0001999600079984,0,0.000737503801108105,61.2357956825111,97.3841922782264,123.343342559042,187.552366843889,14.4388561467183,NA,"Negative-binomial overdispersed count null using method-of-moments size; smoothing baseline re-estimated.","tests extra-Poisson variance",FALSE +"block_residual_refit",5000,315.727660603729,2,0.000599880023995201,0.0001999600079984,4.84455271182235e-05,0.0014441827958669,111.129777725806,187.002073839735,234.873979019818,388.221228496508,NA,8,"Block-permuted Pearson-residual rate perturbation followed by Poisson sampling; smoothing baseline re-estimated.","adversarial local-structure / baseline-misspecification diagnostic",TRUE diff --git a/tables_r/statistical_audit/calibrated_global_evidence.csv b/tables_r/statistical_audit/calibrated_global_evidence.csv new file mode 100644 index 0000000..5f0ce19 --- /dev/null +++ b/tables_r/statistical_audit/calibrated_global_evidence.csv @@ -0,0 +1,2 @@ +"analysis_id","observed_deltaD","scan_global_p","scan_global_tail_count","scan_global_B","family_max_p","holm_p","bonferroni_p","bh_fdr","by_fdr","family_size","worst_fwer_adjusted_p","worst_alternative_null","worst_alternative_null_p","worst_alternative_null_B","largest_reported_global_p","largest_reported_global_p_note","calibration_fpr_alpha_0_05","calibration_ci_lo","calibration_ci_hi","calibration_compatible","holdout_designs","holdout_sig_direction_consistent","holdout_median_historical_fixed_p","holdout_median_refit_p","holdout_median_conservative_p","heldout_model_loglik_gain_M1_minus_M0","bin_reference_fraction","bootstrap_peak_reference_fraction","specification_reference_fraction","injection_null_detection_prob","injection_null_n" +"fe_ion2_wn_bin160_sigma6_degree1",315.727660603729,0.0001999600079984,0,5000,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,0.0117976404719056,"block_residual_refit",0.000599880023995201,5000,0.0117976404719056,"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",0.0528,0.048497623646849,0.057365099095839,TRUE,8,7,0.000999000999000999,0.000999000999000999,0.000999000999000999,-12.1944606016229,0.7,1,0.893617021276596,0.0455,2000 diff --git a/tables_r/statistical_audit/failed_claims.csv b/tables_r/statistical_audit/failed_claims.csv new file mode 100644 index 0000000..800663e --- /dev/null +++ b/tables_r/statistical_audit/failed_claims.csv @@ -0,0 +1,6 @@ +"claim","analysis","result","threshold","verdict","limitations","output_reference" +"Frozen blocked holdout replication","run_holdout_replication","7/8 designs have conservative p<=0.05 and positive direction; median conservative p=0.000999","strict: every declared blocked design conservative p<=0.05 with positive direction","mixed","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" +"Held-out predictive transfer of smooth-plus-periodic model","model_comparison","held-out log-likelihood gain M1-M0 = -12.1945","gain > 0","fail","explicitly preserves the prior failure if the large in-sample model gain does not transfer out of sample","model_comparison.csv" +"Declared Fe II bin stability","bin_stability_summary","70.0% of audit-declared bins select the reference region",">= 80% (audit-declared descriptive threshold)","fail","descriptive robustness criterion, not historical preregistration or a universal physical threshold","bin_stability_summary.csv" +"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" +"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" diff --git a/tables_r/statistical_audit/final_adversarial_claim_matrix.csv b/tables_r/statistical_audit/final_adversarial_claim_matrix.csv new file mode 100644 index 0000000..885a1aa --- /dev/null +++ b/tables_r/statistical_audit/final_adversarial_claim_matrix.csv @@ -0,0 +1,13 @@ +"claim","analysis","result","threshold","verdict","limitations","output_reference" +"Fe II full-scan look-elsewhere evidence","significance_results","scan-global p=0.0002 (0/5000 exceedances)","p <= 0.05","pass","empirical p is bounded by Monte Carlo resolution; zero exceedances are not an exact smaller p","significance_results.csv" +"Fe II full declared-family FWER correction","global_multiple_testing","family-max=0.0002; Holm=0.0118; Bonferroni=0.0118; BH=0.0002034; BY=0.0009485; family=59","all reported FWER-adjusted p-values <= 0.05","pass","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" +"Fe II robustness to declared alternative null models","alternative_null_results","worst-case p=0.0005999 under block_residual_refit","all declared null-model scan-global p <= 0.05","pass","post-signal adversarial stress tests; not historical preregistration and not a claim that one generator is the unique physical null","alternative_null_results.csv" +"Synthetic-null false-positive calibration","null_calibration","FPR@0.05=0.0528 [0.0485, 0.0574]","nominal 0.05 lies within calibration CI","pass","calibration tests the canonical fitted-Poisson synthetic null; alternative-null robustness is reported separately","null_calibration.csv" +"Frozen blocked holdout replication","run_holdout_replication","7/8 designs have conservative p<=0.05 and positive direction; median conservative p=0.000999","strict: every declared blocked design conservative p<=0.05 with positive direction","mixed","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" +"Held-out predictive transfer of smooth-plus-periodic model","model_comparison","held-out log-likelihood gain M1-M0 = -12.1945","gain > 0","fail","explicitly preserves the prior failure if the large in-sample model gain does not transfer out of sample","model_comparison.csv" +"Declared Fe II bin stability","bin_stability_summary","70.0% of audit-declared bins select the reference region",">= 80% (audit-declared descriptive threshold)","fail","descriptive robustness criterion, not historical preregistration or a universal physical threshold","bin_stability_summary.csv" +"Bootstrap peak-region stability","peak_stability","100.0% of primary bootstrap resamples select the 2% reference region",">= 80% (audit-declared descriptive threshold)","pass","resampling the same dataset does not constitute independent replication","peak_stability.csv" +"Specification multiverse stability","run_model_sensitivity","89.4% of declared specifications select the reference region",">= 80% (audit-declared descriptive threshold)","pass","the threshold is a robustness convention; the full specification table remains primary evidence","specification_results.csv" +"Injection null type-I behaviour","run_injection_recovery","A=0 detection probability=0.0455 (2000 simulations)","<= 0.05","pass","finite simulation uncertainty applies; null_calibration.csv is the dedicated type-I assessment","injection_recovery.csv" +"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" +"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" diff --git a/tables_r/statistical_audit/holdout_results.csv b/tables_r/statistical_audit/holdout_results.csv index 5c5b9ad..e7bd745 100644 --- a/tables_r/statistical_audit/holdout_results.csv +++ b/tables_r/statistical_audit/holdout_results.csv @@ -1,9 +1,9 @@ -"design","n_train","n_test","k_lock","test_deltaD_lockedk","test_amplitude_lockedk","test_phase_lockedk","train_phase","phase_diff","predictive_loglik_gain","direction_consistent","fixed_k_test_p","test_B","exploratory_rescan_k","exploratory_rescan_deltaD" -"lower_train_upper_test",80,80,31.8037214885954,159.927665124833,0.228925607199974,0.49997989858866,1.2217643859818,0.721784487393137,79.9638325624163,TRUE,0.0384615384615385,25,23.7551020408163,184.703270475327 -"upper_train_lower_test",80,80,23.7551020408163,12.1468995373982,0.0855839522697742,1.59485001169588,-1.37363197920463,2.96848199090051,6.07344976869916,TRUE,0.0384615384615385,25,31.8037214885954,180.195668047178 -"alternating_blocks",80,80,30.9765906362545,188.015570219377,0.246199108508106,3.03606653089951,2.26751525753209,0.768551273367418,94.0077851096889,TRUE,0.0384615384615385,25,31.1674669867947,190.681194494396 -"kfold_1_of_5",128,32,31.1992797118848,25.0259850268126,0.182524886288597,0.0868142237215331,0.590442657954158,0.503628434232625,12.5129925134063,TRUE,0.0384615384615385,25,26.5546218487395,31.8077700540645 -"kfold_2_of_5",128,32,31.421968787515,125.93092915412,0.447538611852894,-1.46534859221941,-1.74127464786816,0.275926055648743,62.9654645770601,TRUE,0.0384615384615385,25,30.0222088835534,130.004085984641 -"kfold_3_of_5",128,32,31.3265306122449,84.8422055799866,0.377720434126854,-0.537173445299685,-0.732750338381396,0.19557689308171,42.4211027899933,TRUE,0.0384615384615385,25,27.5726290516206,103.641403115694 -"kfold_4_of_5",128,32,31.3265306122449,0.608418496619265,0.0382006562011213,-1.71472503014343,-0.637411312162968,1.07731371798047,0.304209248309633,TRUE,0.615384615384615,25,20.0648259303721,93.6675111183221 -"kfold_5_of_5",128,32,30.9129651860744,69.3641082305783,0.184046256680928,-2.60903888785698,3.122150969199,0.551995450123602,34.6820541152892,TRUE,0.0384615384615385,25,17.1062424969988,133.675314654169 +"design","n_train","n_test","k_lock","test_deltaD_lockedk","test_amplitude_lockedk","test_phase_lockedk","train_phase","phase_diff","predictive_loglik_gain","direction_consistent","fixed_k_test_p","fixed_k_test_p_refit","conservative_fixed_k_test_p","test_B","exploratory_rescan_k","exploratory_rescan_deltaD" +"lower_train_upper_test",80,80,31.8037214885954,159.927665124833,0.228925607199974,0.49997989858866,1.2217643859818,0.721784487393137,79.9638325624163,TRUE,0.000999000999000999,0.000999000999000999,0.000999000999000999,1000,23.7551020408163,184.703270475327 +"upper_train_lower_test",80,80,23.7551020408163,12.1468995373982,0.0855839522697742,1.59485001169588,-1.37363197920463,2.96848199090051,6.07344976869916,TRUE,0.001998001998002,0.000999000999000999,0.001998001998002,1000,31.8037214885954,180.195668047178 +"alternating_blocks",80,80,30.9765906362545,188.015570219377,0.246199108508106,3.03606653089951,2.26751525753209,0.768551273367418,94.0077851096889,TRUE,0.000999000999000999,0.000999000999000999,0.000999000999000999,1000,31.1674669867947,190.681194494396 +"kfold_1_of_5",128,32,31.1992797118848,25.0259850268126,0.182524886288597,0.0868142237215331,0.590442657954158,0.503628434232625,12.5129925134063,TRUE,0.000999000999000999,0.000999000999000999,0.000999000999000999,1000,26.5546218487395,31.8077700540645 +"kfold_2_of_5",128,32,31.421968787515,125.93092915412,0.447538611852894,-1.46534859221941,-1.74127464786816,0.275926055648743,62.9654645770601,TRUE,0.000999000999000999,0.000999000999000999,0.000999000999000999,1000,30.0222088835534,130.004085984641 +"kfold_3_of_5",128,32,31.3265306122449,84.8422055799866,0.377720434126854,-0.537173445299685,-0.732750338381396,0.19557689308171,42.4211027899933,TRUE,0.000999000999000999,0.000999000999000999,0.000999000999000999,1000,27.5726290516206,103.641403115694 +"kfold_4_of_5",128,32,31.3265306122449,0.608418496619265,0.0382006562011213,-1.71472503014343,-0.637411312162968,1.07731371798047,0.304209248309633,TRUE,0.729270729270729,0.751248751248751,0.751248751248751,1000,20.0648259303721,93.6675111183221 +"kfold_5_of_5",128,32,30.9129651860744,69.3641082305783,0.184046256680928,-2.60903888785698,3.122150969199,0.551995450123602,34.6820541152892,TRUE,0.000999000999000999,0.000999000999000999,0.000999000999000999,1000,17.1062424969988,133.675314654169 diff --git a/tables_r/statistical_audit/injection_recovery.csv b/tables_r/statistical_audit/injection_recovery.csv index b7aa20a..2633907 100644 --- a/tables_r/statistical_audit/injection_recovery.csv +++ b/tables_r/statistical_audit/injection_recovery.csv @@ -1,17 +1,33 @@ "freq_name","k_injected","amplitude","n_sims","detection_prob","det_ci_lo","det_ci_hi","correct_region_prob","sig_correct_prob","sigcor_ci_lo","sigcor_ci_hi","false_localisation_prob","median_recovered_k","freq_bias","amplitude_bias","alpha","null_n" -"fe_reference",31.3265306122449,0,8,0,0,0.369416647552819,0,0,0,0.369416647552819,0,52.2752100840336,20.9486794717887,0.0436049366503891,0.05,25 -"lower_control",12,0,8,0,0,0.369416647552819,0,0,0,0.369416647552819,0,27.6680672268908,15.6680672268908,0.043829520486473,0.05,25 -"higher_control",55,0,8,0,0,0.369416647552819,0.125,0,0,0.369416647552819,0,45.1650660264106,-9.83493397358944,0.0404788629380577,0.05,25 -"off_grid",31.7337755102041,0,8,0,0,0.369416647552819,0.125,0,0,0.369416647552819,0,32.7740096038415,1.04023409363746,0.0368152947303226,0.05,25 -"fe_reference",31.3265306122449,0.03,8,0.125,0.00315972353125191,0.526509670875207,0,0,0,0.369416647552819,0.125,36.5438175270108,5.2172869147659,0.015117831425137,0.05,25 -"lower_control",12,0.03,8,0.125,0.00315972353125191,0.526509670875207,0.125,0,0,0.369416647552819,0.125,29.1155462184874,17.1155462184874,0.0221649903866946,0.05,25 -"higher_control",55,0.03,8,0.125,0.00315972353125191,0.526509670875207,0.125,0,0,0.369416647552819,0.125,42.6836734693878,-12.3163265306122,0.016557197981419,0.05,25 -"off_grid",31.7337755102041,0.03,8,0.5,0.157012770487058,0.842987229512942,0.375,0.25,0.0318540262499442,0.650855794412824,0.25,32.2490996398559,0.515324129651862,0.0253871695052827,0.05,25 -"fe_reference",31.3265306122449,0.05,8,0.375,0.0852334141372536,0.755136783633448,0.625,0.375,0.0852334141372536,0.755136783633448,0,31.1674669867947,-0.159063625450184,0.00301977708662526,0.05,25 -"lower_control",12,0.05,8,0.5,0.157012770487058,0.842987229512942,0.125,0.125,0.00315972353125191,0.526509670875207,0.375,12.2866146458583,0.286614645858343,0.00820057325814869,0.05,25 -"higher_control",55,0.05,8,0.5,0.157012770487058,0.842987229512942,0.875,0.5,0.157012770487058,0.842987229512942,0,55.0747298919568,0.0747298919567783,0.00381707788314227,0.05,25 -"off_grid",31.7337755102041,0.05,8,0.5,0.157012770487058,0.842987229512942,1,0.5,0.157012770487058,0.842987229512942,0,31.9627851140456,0.229009603841536,0.00781105933072553,0.05,25 -"fe_reference",31.3265306122449,0.1,8,1,0.630583352447181,1,1,1,0.630583352447181,1,0,31.1992797118848,-0.127250900360146,-0.00247420313202751,0.05,25 -"lower_control",12,0.1,8,1,0.630583352447181,1,0.625,0.625,0.244863216366552,0.914766585862746,0.375,12.0639255702281,0.0639255702280899,-0.00437200529793286,0.05,25 -"higher_control",55,0.1,8,1,0.630583352447181,1,1,1,0.630583352447181,1,0,54.9474789915966,-0.0525210084033603,0.00050375557064572,0.05,25 -"off_grid",31.7337755102041,0.1,8,1,0.630583352447181,1,1,1,0.630583352447181,1,0,31.8196278511405,0.0858523409363734,0.0104140043824759,0.05,25 +"fe_reference",31.3265306122449,0,2000,0.0455,0.0367892780988175,0.0555721910065792,0.013,5e-04,1.26588238685579e-05,0.00278263983465895,0.045,40.2022809123649,8.87575030012005,0.0438334109675285,0.05,1000 +"lower_control",12,0,2000,0.046,0.037240902034369,0.0561188305434364,0.006,0,0,0.00184273979340594,0.046,40.9975990396158,28.9975990396158,0.0437636538244677,0.05,1000 +"higher_control",55,0,2000,0.0525,0.0431376180846189,0.063199902483287,0.0285,0.002,0.000545193125210999,0.00511280860534621,0.0505,42.5564225690276,-12.4435774309724,0.0436873075318457,0.05,1000 +"off_grid",31.7337755102041,0,2000,0.0525,0.0431376180846189,0.063199902483287,0.0155,0.001,0.00012112759055683,0.00360762856982853,0.0515,40.0909363745498,8.35716086434574,0.0439056577701072,0.05,1000 +"fe_reference",31.3265306122449,0.01,2000,0.0635,0.0532088422314592,0.0750922552651426,0.0265,0.0025,0.00081222597519391,0.00582444533451852,0.061,40.7271908763505,9.40066026410564,0.0342200706642423,0.05,1000 +"lower_control",12,0.01,2000,0.052,0.0426824397973374,0.0626567679337087,0.0155,5e-04,1.26588238685579e-05,0.00278263983465895,0.0515,40.2022809123649,28.2022809123649,0.0341855015983926,0.05,1000 +"higher_control",55,0.01,2000,0.0595,0.0495346672333621,0.0707795372688371,0.0445,0.002,0.000545193125210999,0.00511280860534621,0.0575,42.1587635054022,-12.8412364945978,0.03422363464716,0.05,1000 +"off_grid",31.7337755102041,0.01,2000,0.061,0.0509109937256597,0.0723982793153182,0.036,0.0045,0.00205968852849036,0.00852514124099513,0.0565,39.5660264105642,7.83225090036014,0.0341637178684732,0.05,1000 +"fe_reference",31.3265306122449,0.02,2000,0.071,0.0601297870904785,0.0831470663609534,0.074,0.0095,0.00572903552296345,0.0147958321722885,0.0615,36.734693877551,5.40816326530613,0.0249109518101419,0.05,1000 +"lower_control",12,0.02,2000,0.0855,0.0736096478366675,0.0986212441907105,0.039,0.004,0.00172844957347615,0.00786630382853737,0.0815,34.030612244898,22.030612244898,0.0250575486142362,0.05,1000 +"higher_control",55,0.02,2000,0.082,0.0703452001669976,0.0948966967162819,0.1145,0.012,0.0077033433702185,0.017802731356322,0.07,47.280612244898,-7.71938775510204,0.0248761047120981,0.05,1000 +"off_grid",31.7337755102041,0.02,2000,0.085,0.0731429139559128,0.0980895471881371,0.088,0.0155,0.0105551589199576,0.0219293891896312,0.0695,36.3688475390156,4.63507202881152,0.0247135293272626,0.05,1000 +"fe_reference",31.3265306122449,0.03,2000,0.1415,0.126508942339429,0.157550228613897,0.1905,0.063,0.0527488858587171,0.0745538425895165,0.0785,31.9468787515006,0.6203481392557,0.0167004956201055,0.05,1000 +"lower_control",12,0.03,2000,0.159,0.143232615616254,0.175773867888546,0.0955,0.0275,0.0207824551794987,0.0356459912315418,0.1315,20.5897358943577,8.58973589435774,0.0174505364111435,0.05,1000 +"higher_control",55,0.03,2000,0.143,0.127939456900426,0.159115185914567,0.2695,0.0765,0.0652286499566994,0.089030648346499,0.0665,54.1680672268908,-0.831932773109244,0.0170132337854053,0.05,1000 +"off_grid",31.7337755102041,0.03,2000,0.127,0.112711881889353,0.142391201459244,0.2055,0.0535,0.0440487130191427,0.0642854427090288,0.0735,32.3922569027611,0.658481392557025,0.0169307557817861,0.05,1000 +"fe_reference",31.3265306122449,0.04,2000,0.2825,0.262849396372472,0.302790165114762,0.377,0.1765,0.160025106270065,0.193928948484415,0.106,31.4855942376951,0.15906362545018,0.0103275659728751,0.05,1000 +"lower_control",12,0.04,2000,0.2855,0.265782073918126,0.305848633341781,0.1815,0.0965,0.0839078667396506,0.110288726483738,0.189,12.6047418967587,0.604741896758703,0.0108024807268139,0.05,1000 +"higher_control",55,0.04,2000,0.289,0.269204882171254,0.309415497251195,0.453,0.211,0.193302025510116,0.229549212967665,0.078,54.8361344537815,-0.163865546218489,0.0104128994547745,0.05,1000 +"off_grid",31.7337755102041,0.04,2000,0.267,0.247714671970102,0.286970664778765,0.3985,0.1795,0.162910106917439,0.197034981399766,0.0875,31.8673469387755,0.133571428571429,0.00997767227396348,0.05,1000 +"fe_reference",31.3265306122449,0.05,2000,0.5105,0.488339967855429,0.532629206779004,0.586,0.3695,0.348300784811261,0.391082531564708,0.141,31.390156062425,0.0636254501800693,0.00606688327995275,0.05,1000 +"lower_control",12,0.05,2000,0.5275,0.505343313939927,0.54957595142096,0.313,0.22,0.202016098561306,0.23880840862764,0.3075,12.1752701080432,0.175270108043216,0.00709604720519226,0.05,1000 +"higher_control",55,0.05,2000,0.4995,0.477351299108812,0.521650168760985,0.692,0.429,0.407176562246782,0.451031906658574,0.0705,54.9633853541417,-0.036614645858343,0.00577468563847403,0.05,1000 +"off_grid",31.7337755102041,0.05,2000,0.532,0.509848448653203,0.55405760484072,0.6185,0.4125,0.390816749868454,0.434440186580811,0.1195,31.8037214885954,0.0699459783913561,0.00661143004207381,0.05,1000 +"fe_reference",31.3265306122449,0.075,2000,0.947,0.936257207001227,0.956406956407424,0.929,0.891,0.876514699801219,0.90432726315567,0.056,31.3265306122449,0,0.00298846889005866,0.05,1000 +"lower_control",12,0.075,2000,0.9555,0.945521988642517,0.964113055495469,0.56,0.548,0.525881173806467,0.569977901072997,0.4075,12.016206482593,0.0162064825930361,0.00368742899470667,0.05,1000 +"higher_control",55,0.075,2000,0.9425,0.931381664528725,0.952297522619749,0.9715,0.9285,0.916317290730765,0.939407458848455,0.014,54.9951980792317,-0.00480192076830832,0.00274857562313875,0.05,1000 +"off_grid",31.7337755102041,0.075,2000,0.942,0.930841048250829,0.951839794677264,0.9285,0.891,0.876514699801219,0.90432726315567,0.051,31.7400960384154,0.00632052821128681,0.00231007205384336,0.05,1000 +"fe_reference",31.3265306122449,0.1,2000,0.9995,0.997217360165341,0.999987341176131,0.9915,0.991,0.985813285528485,0.994657556976605,0.0085,31.3265306122449,0,0.00188386528953283,0.05,1000 +"lower_control",12,0.1,2000,0.9995,0.997217360165341,0.999987341176131,0.711,0.711,0.690584502748805,0.730795117828746,0.2885,12.016206482593,0.0162064825930361,0.00260688558144129,0.05,1000 +"higher_control",55,0.1,2000,0.9995,0.997217360165341,0.999987341176131,0.999,0.9985,0.99562267989547,0.999690557050688,0.001,54.9951980792317,-0.00480192076830832,0.00175964407518051,0.05,1000 +"off_grid",31.7337755102041,0.1,2000,0.999,0.996392371430171,0.999878872409443,0.9945,0.9935,0.988910450728344,0.996534618791332,0.0055,31.7400960384154,0.00632052821128681,0.00280297615432337,0.05,1000 diff --git a/tables_r/statistical_audit/model_comparison.csv b/tables_r/statistical_audit/model_comparison.csv index 1167025..426211b 100644 --- a/tables_r/statistical_audit/model_comparison.csv +++ b/tables_r/statistical_audit/model_comparison.csv @@ -1,3 +1,3 @@ -"model","loglik","n_params","n_bins","deviance","AIC","AICc","BIC","heldout_loglik_test","deltaAIC","deltaBIC","selected_k","heldout_k_lock","deviance_difference","note" -"M0_smooth_null",-953.216372386261,2,160,1011.9344557793,1910.43274477252,1910.50917789354,1916.58309240299,-523.841554349368,311.727660603729,305.577312973261,31.3265306122449,31.9627851140456,315.727660603729,"k selected by scan; AIC/BIC do not fully account for look-elsewhere" -"M1_smooth_plus_logperiodic",-795.352542084396,4,160,696.206795175567,1598.70508416879,1598.96314868492,1611.00577942973,-560.677192339407,0,0,31.3265306122449,31.9627851140456,315.727660603729,"k selected by scan; AIC/BIC do not fully account for look-elsewhere" +"model","loglik","n_params","n_bins","deviance","AIC","AICc","BIC","heldout_loglik_test","deltaAIC","deltaBIC","selected_k","heldout_k_lock","heldout_train_center","heldout_train_scale","deviance_difference","note" +"M0_smooth_null",-953.216372386261,2,160,1011.9344557793,1910.43274477252,1910.50917789354,1916.58309240299,-605.656386645527,311.727660603729,305.577312973261,31.3265306122449,31.9627851140456,9.19388764423775,0.312842672339254,315.727660603729,"k selected by scan; AIC/BIC do not fully account for look-elsewhere; held-out polynomial transform frozen from training block" +"M1_smooth_plus_logperiodic",-795.352542084396,4,160,696.206795175567,1598.70508416879,1598.96314868492,1611.00577942973,-617.85084724715,0,0,31.3265306122449,31.9627851140456,9.19388764423775,0.312842672339254,315.727660603729,"k selected by scan; AIC/BIC do not fully account for look-elsewhere; held-out polynomial transform frozen from training block" diff --git a/tables_r/statistical_audit/multiple_testing.csv b/tables_r/statistical_audit/multiple_testing.csv index f703ae0..41a04f6 100644 --- a/tables_r/statistical_audit/multiple_testing.csv +++ b/tables_r/statistical_audit/multiple_testing.csv @@ -1,18 +1,60 @@ -"analysis_id","species","source","bins","sigma","degree","observed_stat","scan_global_p","family_max_p","bonferroni_p","holm_p","bh_fdr","family_size","dependence_note" -"fe_ion2_wn_bingrid60","Fe","wavenumber",60,6,1,964.939313927201,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"fe_ion2_wn_bingrid80","Fe","wavenumber",80,6,1,637.189834656421,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"fe_ion2_wn_bingrid100","Fe","wavenumber",100,6,1,406.356274810153,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"fe_ion2_wn_bingrid120","Fe","wavenumber",120,6,1,355.744301050131,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"fe_ion2_wn_bingrid140","Fe","wavenumber",140,6,1,327.831699186689,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"fe_ion2_wn_bingrid160","Fe","wavenumber",160,6,1,315.727660603729,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"fe_ion2_wn_bingrid180","Fe","wavenumber",180,6,1,287.502616304374,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"fe_ion2_wn_bingrid200","Fe","wavenumber",200,6,1,259.178756760527,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"fe_ion2_wn_bingrid220","Fe","wavenumber",220,6,1,225.35158661837,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"fe_ion2_wn_bingrid240","Fe","wavenumber",240,6,1,194.310229119281,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"fe_ion2_observed_bin160","Fe","observed",160,6,1,315.376115676902,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"fe_ion2_ritz_bin160","Fe","ritz",160,6,1,346.338782567562,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"cr_ion2_wn_bin160","Cr","wavenumber",160,6,1,228.486019036192,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"mn_ion2_wn_bin160","Mn","wavenumber",160,6,1,89.4593944441067,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"co_ion2_wn_bin160","Co","wavenumber",160,6,1,327.153172835048,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"ni_ion2_wn_bin160","Ni","wavenumber",160,6,1,24.149190957751,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" -"ti_ion2_wn_bin160","Ti","wavenumber",160,6,1,46.9882628355725,0.0625,0.0625,1,1,0.0625,17,"shared NIST line list; corrections approximate, see header" +"species","source","bins","sigma","degree","family_component","analysis_id","observed_stat","scan_global_p","family_max_p","bonferroni_p","holm_p","bh_fdr","by_fdr","family_size","dependence_note" +"Fe","wavenumber",160,6,1,"bin_grid+specification_multiverse","fe_ion2_wn_bin160_sig6_deg1",315.727660603729,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Co","wavenumber",160,6,1,"neighbour_control","co_ion2_wn_bin160_sig6_deg1",327.153172835048,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Cr","wavenumber",160,6,1,"neighbour_control","cr_ion2_wn_bin160_sig6_deg1",228.486019036192,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","observed",160,6,1,"source_replication","fe_ion2_observed_bin160_sig6_deg1",315.376115676902,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","ritz",160,6,1,"source_replication","fe_ion2_ritz_bin160_sig6_deg1",346.338782567562,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",60,6,1,"bin_grid","fe_ion2_wn_bin60_sig6_deg1",964.939313927201,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",80,6,1,"bin_grid","fe_ion2_wn_bin80_sig6_deg1",637.189834656421,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",100,6,1,"bin_grid","fe_ion2_wn_bin100_sig6_deg1",406.356274810153,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,4,0,"specification_multiverse","fe_ion2_wn_bin120_sig4_deg0",286.874999417522,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,4,1,"specification_multiverse","fe_ion2_wn_bin120_sig4_deg1",287.623695108705,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,4,2,"specification_multiverse","fe_ion2_wn_bin120_sig4_deg2",288.954398265615,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,5,0,"specification_multiverse","fe_ion2_wn_bin120_sig5_deg0",333.946669109476,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,5,1,"specification_multiverse","fe_ion2_wn_bin120_sig5_deg1",335.383912820845,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,5,2,"specification_multiverse","fe_ion2_wn_bin120_sig5_deg2",338.1520352881,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,6,0,"specification_multiverse","fe_ion2_wn_bin120_sig6_deg0",353.318474822445,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,6,1,"bin_grid+specification_multiverse","fe_ion2_wn_bin120_sig6_deg1",355.744301050131,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,6,2,"specification_multiverse","fe_ion2_wn_bin120_sig6_deg2",360.843774114349,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,7,0,"specification_multiverse","fe_ion2_wn_bin120_sig7_deg0",381.832625700486,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,7,1,"specification_multiverse","fe_ion2_wn_bin120_sig7_deg1",377.549280874837,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,7,2,"specification_multiverse","fe_ion2_wn_bin120_sig7_deg2",373.965760764014,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,8,0,"specification_multiverse","fe_ion2_wn_bin120_sig8_deg0",517.388433563542,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,8,1,"specification_multiverse","fe_ion2_wn_bin120_sig8_deg1",508.250627700088,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",120,8,2,"specification_multiverse","fe_ion2_wn_bin120_sig8_deg2",497.11296063639,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",140,6,1,"bin_grid","fe_ion2_wn_bin140_sig6_deg1",327.831699186689,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",160,4,0,"specification_multiverse","fe_ion2_wn_bin160_sig4_deg0",191.523785068614,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",160,4,1,"specification_multiverse","fe_ion2_wn_bin160_sig4_deg1",191.874169361293,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",160,4,2,"specification_multiverse","fe_ion2_wn_bin160_sig4_deg2",192.328355519125,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",160,5,0,"specification_multiverse","fe_ion2_wn_bin160_sig5_deg0",267.199946889302,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",160,5,1,"specification_multiverse","fe_ion2_wn_bin160_sig5_deg1",267.804479299904,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",160,5,2,"specification_multiverse","fe_ion2_wn_bin160_sig5_deg2",268.719880745123,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",160,6,0,"specification_multiverse","fe_ion2_wn_bin160_sig6_deg0",314.670991439908,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",160,6,2,"specification_multiverse","fe_ion2_wn_bin160_sig6_deg2",317.42904914475,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",160,7,0,"specification_multiverse","fe_ion2_wn_bin160_sig7_deg0",339.977141567488,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",160,7,1,"specification_multiverse","fe_ion2_wn_bin160_sig7_deg1",341.623637677083,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",160,7,2,"specification_multiverse","fe_ion2_wn_bin160_sig7_deg2",344.640594959049,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",160,8,0,"specification_multiverse","fe_ion2_wn_bin160_sig8_deg0",352.926421741162,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",160,8,1,"specification_multiverse","fe_ion2_wn_bin160_sig8_deg1",355.416151088066,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",160,8,2,"specification_multiverse","fe_ion2_wn_bin160_sig8_deg2",360.139096405204,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",180,6,1,"bin_grid","fe_ion2_wn_bin180_sig6_deg1",287.502616304374,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,4,0,"specification_multiverse","fe_ion2_wn_bin200_sig4_deg0",120.022205448983,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,4,1,"specification_multiverse","fe_ion2_wn_bin200_sig4_deg1",120.276693357547,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,4,2,"specification_multiverse","fe_ion2_wn_bin200_sig4_deg2",120.317384883968,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,5,0,"specification_multiverse","fe_ion2_wn_bin200_sig5_deg0",194.465004004777,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,5,1,"specification_multiverse","fe_ion2_wn_bin200_sig5_deg1",194.861685671457,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,5,2,"specification_multiverse","fe_ion2_wn_bin200_sig5_deg2",195.016893114809,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,6,0,"specification_multiverse","fe_ion2_wn_bin200_sig6_deg0",258.556795908811,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,6,1,"bin_grid+specification_multiverse","fe_ion2_wn_bin200_sig6_deg1",259.178756760527,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,6,2,"specification_multiverse","fe_ion2_wn_bin200_sig6_deg2",259.548360871709,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,7,0,"specification_multiverse","fe_ion2_wn_bin200_sig7_deg0",304.248319746702,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,7,1,"specification_multiverse","fe_ion2_wn_bin200_sig7_deg1",305.158000043472,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,7,2,"specification_multiverse","fe_ion2_wn_bin200_sig7_deg2",306.019884666066,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,8,0,"specification_multiverse","fe_ion2_wn_bin200_sig8_deg0",332.939018597347,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,8,1,"specification_multiverse","fe_ion2_wn_bin200_sig8_deg1",334.301739521917,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",200,8,2,"specification_multiverse","fe_ion2_wn_bin200_sig8_deg2",335.857373193256,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",220,6,1,"bin_grid","fe_ion2_wn_bin220_sig6_deg1",225.35158661837,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Fe","wavenumber",240,6,1,"bin_grid","fe_ion2_wn_bin240_sig6_deg1",194.310229119281,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Mn","wavenumber",160,6,1,"neighbour_control","mn_ion2_wn_bin160_sig6_deg1",89.4593944441067,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Ni","wavenumber",160,6,1,"neighbour_control","ni_ion2_wn_bin160_sig6_deg1",24.149190957751,0.000799840031993601,0.0351929614077185,0.0471905618876225,0.0117976404719056,0.000799840031993601,0.00372981703362133,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" +"Ti","wavenumber",160,6,1,"neighbour_control","ti_ion2_wn_bin160_sig6_deg1",46.9882628355725,0.0001999600079984,0.0001999600079984,0.0117976404719056,0.0117976404719056,0.0002034075943432,0.000948531055964046,59,"shared transition lists create dependence; Holm/Bonferroni valid without independence; BY controls FDR under arbitrary dependence; family-max uses approximate marginal-null coupling" diff --git a/tables_r/statistical_audit/null_calibration.csv b/tables_r/statistical_audit/null_calibration.csv index 9ad878b..b35ca60 100644 --- a/tables_r/statistical_audit/null_calibration.csv +++ b/tables_r/statistical_audit/null_calibration.csv @@ -1,4 +1,4 @@ "nominal_alpha","n_outer_sims","inner_null_n","observed_count","expected_count","observed_fpr","ci_lo","ci_hi","calibration_ratio","mc_uncertainty","compatible" -0.1,25,25,0,2.5,0,0,0.137185171530713,0,0.06,TRUE -0.05,25,25,0,1.25,0,0,0.137185171530713,0,0.0435889894354067,TRUE -0.01,25,25,0,0.25,0,0,0.137185171530713,0,0.0198997487421324,TRUE +0.1,10000,5000,1034,1000,0.1034,0.0974984115194039,0.109533672417909,1.034,0.003,TRUE +0.05,10000,5000,528,500,0.0528,0.048497623646849,0.057365099095839,1.056,0.00217944947177034,TRUE +0.01,10000,5000,112,100,0.0112,0.0092307351214203,0.013461140389643,1.12,0.00099498743710662,TRUE diff --git a/tables_r/statistical_audit/observed_ritz_replication.csv b/tables_r/statistical_audit/observed_ritz_replication.csv index 598c169..db95b80 100644 --- a/tables_r/statistical_audit/observed_ritz_replication.csv +++ b/tables_r/statistical_audit/observed_ritz_replication.csv @@ -1,4 +1,4 @@ "source","n_lines","ell_min","ell_max","k_best","delta_log_x","scale_ratio","amplitude","deltaD","global_p","null_n","ci_lo","ci_hi","same_reference_region","abs_diff_from_wn","rel_diff_from_wn","overlap_with_wavenumber" -"wavenumber",9447,8.65875967500558,10.8128190954591,31.3265306122449,0.200570736190098,1.22210005588453,0.261305834349196,315.727660603729,0.0384615384615385,25,31.1658763505402,31.4855942376951,TRUE,0,0,9447 -"observed",8994,8.65903289132403,10.8131428632729,31.3265306122449,0.200570736190098,1.22210005588453,0.267872314470838,315.376115676902,0.0384615384615385,25,31.1674669867947,31.5174069627851,TRUE,0,0,1 -"ritz",9858,8.65903557021301,10.8131420591133,31.3265306122449,0.200570736190098,1.22210005588453,0.268053804171389,346.338782567562,0.0384615384615385,25,31.2302971188475,31.4863895558223,TRUE,0,0,0 +"wavenumber",9447,8.65875967500558,10.8128190954591,31.3265306122449,0.200570736190098,1.22210005588453,0.261305834349196,315.727660603729,0.0001999600079984,5000,31.1674669867947,31.4855942376951,TRUE,0,0,9447 +"observed",8994,8.65903289132403,10.8131428632729,31.3265306122449,0.200570736190098,1.22210005588453,0.267872314470838,315.376115676902,0.0001999600079984,5000,31.1356542617047,31.4855942376951,TRUE,0,0,1 +"ritz",9858,8.65903557021301,10.8131420591133,31.3265306122449,0.200570736190098,1.22210005588453,0.268053804171389,346.338782567562,0.0001999600079984,5000,31.1674669867947,31.4855942376951,TRUE,0,0,0 diff --git a/tables_r/statistical_audit/peak_stability.csv b/tables_r/statistical_audit/peak_stability.csv index 0c14424..0d93a8f 100644 --- a/tables_r/statistical_audit/peak_stability.csv +++ b/tables_r/statistical_audit/peak_stability.csv @@ -1,5 +1,5 @@ "tolerance_type","tolerance","is_primary","k_reference","n_valid","pct_in_reference","pct_largest_competitor","competitor_centre","median_k","iqr_k","ci_lo","ci_hi","stability_class" -"relative",0.01,FALSE,31.3265306122449,80,1,0,NA,31.3265306122449,0.103391356542627,31.1356542617047,31.4855942376951,"high" -"relative",0.02,TRUE,31.3265306122449,80,1,0,NA,31.3265306122449,0.103391356542627,31.1356542617047,31.4855942376951,"high" -"relative",0.05,FALSE,31.3265306122449,80,1,0,NA,31.3265306122449,0.103391356542627,31.1356542617047,31.4855942376951,"high" -"absolute",0.15906362545018,FALSE,31.3265306122449,80,0.8625,NA,NA,31.3265306122449,0.103391356542627,31.1356542617047,31.4855942376951,"high" +"relative",0.01,FALSE,31.3265306122449,5000,0.9996,4e-04,31,31.3265306122449,0.127250900360199,31.1674669867947,31.4855942376951,"high" +"relative",0.02,TRUE,31.3265306122449,5000,1,0,NA,31.3265306122449,0.127250900360199,31.1674669867947,31.4855942376951,"high" +"relative",0.05,FALSE,31.3265306122449,5000,1,0,NA,31.3265306122449,0.127250900360199,31.1674669867947,31.4855942376951,"high" +"absolute",0.15906362545018,FALSE,31.3265306122449,5000,0.9188,NA,NA,31.3265306122449,0.127250900360199,31.1674669867947,31.4855942376951,"high" diff --git a/tables_r/statistical_audit/significance_results.csv b/tables_r/statistical_audit/significance_results.csv index fc90c70..89f9779 100644 --- a/tables_r/statistical_audit/significance_results.csv +++ b/tables_r/statistical_audit/significance_results.csv @@ -1,9 +1,9 @@ "analysis_id","pointwise_statistic","pointwise_tail_count","pointwise_B","pointwise_p","global_statistic","global_tail_count","global_B","global_p","empirical_resolution","zero_exceedance","scan_k_min","scan_k_max","n_scanned_frequencies" -"fe_ion2_wn_bin120",355.744301050131,0,25,0.0384615384615385,355.744301050131,0,25,0.0384615384615385,0.0384615384615385,TRUE,0.5,80,2500 -"fe_ion2_wn_bin160",315.727660603729,0,25,0.0384615384615385,315.727660603729,0,25,0.0384615384615385,0.0384615384615385,TRUE,0.5,80,2500 -"fe_ion2_wn_bin200",259.178756760527,0,25,0.0384615384615385,259.178756760527,0,25,0.0384615384615385,0.0384615384615385,TRUE,0.5,80,2500 -"cr_ion2_wn_bin160",228.486019036192,0,25,0.0384615384615385,228.486019036192,0,25,0.0384615384615385,0.0384615384615385,TRUE,0.5,80,2500 -"mn_ion2_wn_bin160",89.4593944441067,0,25,0.0384615384615385,89.4593944441067,0,25,0.0384615384615385,0.0384615384615385,TRUE,0.5,80,2500 -"co_ion2_wn_bin160",327.153172835048,0,25,0.0384615384615385,327.153172835048,0,25,0.0384615384615385,0.0384615384615385,TRUE,0.5,80,2500 -"ni_ion2_wn_bin160",24.149190957751,0,25,0.0384615384615385,24.149190957751,0,25,0.0384615384615385,0.0384615384615385,TRUE,0.5,80,2500 -"ti_ion2_wn_bin160",46.9882628355725,0,25,0.0384615384615385,46.9882628355725,0,25,0.0384615384615385,0.0384615384615385,TRUE,0.5,80,2500 +"fe_ion2_wn_bin120",355.744301050131,0,5000,0.0001999600079984,355.744301050131,0,5000,0.0001999600079984,0.0001999600079984,TRUE,0.5,80,2500 +"fe_ion2_wn_bin160",315.727660603729,0,5000,0.0001999600079984,315.727660603729,0,5000,0.0001999600079984,0.0001999600079984,TRUE,0.5,80,2500 +"fe_ion2_wn_bin200",259.178756760527,0,5000,0.0001999600079984,259.178756760527,0,5000,0.0001999600079984,0.0001999600079984,TRUE,0.5,80,2500 +"cr_ion2_wn_bin160",228.486019036192,0,5000,0.0001999600079984,228.486019036192,0,5000,0.0001999600079984,0.0001999600079984,TRUE,0.5,80,2500 +"mn_ion2_wn_bin160",89.4593944441067,0,5000,0.0001999600079984,89.4593944441067,0,5000,0.0001999600079984,0.0001999600079984,TRUE,0.5,80,2500 +"co_ion2_wn_bin160",327.153172835048,0,5000,0.0001999600079984,327.153172835048,0,5000,0.0001999600079984,0.0001999600079984,TRUE,0.5,80,2500 +"ni_ion2_wn_bin160",24.149190957751,0,5000,0.0001999600079984,24.149190957751,2,5000,0.000599880023995201,0.0001999600079984,FALSE,0.5,80,2500 +"ti_ion2_wn_bin160",46.9882628355725,0,5000,0.0001999600079984,46.9882628355725,0,5000,0.0001999600079984,0.0001999600079984,TRUE,0.5,80,2500 diff --git a/tables_r/statistical_audit/specification_results.csv b/tables_r/statistical_audit/specification_results.csv index 354fd66..a7d147f 100644 --- a/tables_r/statistical_audit/specification_results.csv +++ b/tables_r/statistical_audit/specification_results.csv @@ -1,21 +1,48 @@ "sigma","degree","bins","source","k_best","deltaD","amplitude","global_p","in_reference_region","converged","n_lines","spec_null_n","spec_id" -7,1,120,"wavenumber",9.72569027611044,377.549280874837,0.314799960165075,0.0769230769230769,FALSE,TRUE,9447,12,"sig7_deg1_bin120_wavenumber" -5,1,200,"wavenumber",31.3265306122449,194.861685671457,0.205367807242324,0.0769230769230769,TRUE,TRUE,9447,12,"sig5_deg1_bin200_wavenumber" -5,2,200,"wavenumber",31.3265306122449,195.016893114809,0.205715675790948,0.0769230769230769,TRUE,TRUE,9447,12,"sig5_deg2_bin200_wavenumber" -6,1,200,"wavenumber",31.3265306122449,259.178756760527,0.236764872134194,0.0769230769230769,TRUE,TRUE,9447,12,"sig6_deg1_bin200_wavenumber" -6,2,200,"wavenumber",31.3265306122449,259.548360871709,0.237211782538629,0.0769230769230769,TRUE,TRUE,9447,12,"sig6_deg2_bin200_wavenumber" -5,1,160,"wavenumber",31.3265306122449,267.804479299904,0.2406364935376,0.0769230769230769,TRUE,TRUE,9447,12,"sig5_deg1_bin160_wavenumber" -5,2,160,"wavenumber",31.3265306122449,268.719880745123,0.241345459579706,0.0769230769230769,TRUE,TRUE,9447,12,"sig5_deg2_bin160_wavenumber" -7,1,200,"wavenumber",31.2947178871549,305.158000043472,0.256894848867292,0.0769230769230769,TRUE,TRUE,9447,12,"sig7_deg1_bin200_wavenumber" -7,2,200,"wavenumber",31.2947178871549,306.019884666066,0.257602568321779,0.0769230769230769,TRUE,TRUE,9447,12,"sig7_deg2_bin200_wavenumber" -6,1,160,"observed",31.3265306122449,315.376115676902,0.267872314470838,0.0769230769230769,TRUE,TRUE,8994,12,"sig6_deg1_bin160_observed" -6,1,160,"wavenumber",31.3265306122449,315.727660603729,0.261305834349196,0.0769230769230769,TRUE,TRUE,9447,12,"sig6_deg1_bin160_wavenumber" -6,2,160,"wavenumber",31.2947178871549,317.42904914475,0.262380813619649,0.0769230769230769,TRUE,TRUE,9447,12,"sig6_deg2_bin160_wavenumber" -5,1,120,"wavenumber",31.3265306122449,335.383912820845,0.269453735775903,0.0769230769230769,TRUE,TRUE,9447,12,"sig5_deg1_bin120_wavenumber" -5,2,120,"wavenumber",31.2947178871549,338.1520352881,0.271009489799855,0.0769230769230769,TRUE,TRUE,9447,12,"sig5_deg2_bin120_wavenumber" -7,1,160,"wavenumber",31.2947178871549,341.623637677083,0.271862261275608,0.0769230769230769,TRUE,TRUE,9447,12,"sig7_deg1_bin160_wavenumber" -7,2,160,"wavenumber",31.2947178871549,344.640594959049,0.273524902123414,0.0769230769230769,TRUE,TRUE,9447,12,"sig7_deg2_bin160_wavenumber" -6,1,160,"ritz",31.3265306122449,346.338782567562,0.268053804171389,0.0769230769230769,TRUE,TRUE,9858,12,"sig6_deg1_bin160_ritz" -6,1,120,"wavenumber",31.3265306122449,355.744301050131,0.27761532337233,0.0769230769230769,TRUE,TRUE,9447,12,"sig6_deg1_bin120_wavenumber" -6,2,120,"wavenumber",31.2947178871549,360.843774114349,0.28016270726631,0.0769230769230769,TRUE,TRUE,9447,12,"sig6_deg2_bin120_wavenumber" -7,2,120,"wavenumber",31.2629051620648,373.965760764014,0.285442245038734,0.0769230769230769,TRUE,TRUE,9447,12,"sig7_deg2_bin120_wavenumber" +7,1,120,"wavenumber",9.72569027611044,377.549280874837,0.314799960165075,0.00199600798403194,FALSE,TRUE,9447,500,"sig7_deg1_bin120_wavenumber" +7,0,120,"wavenumber",9.75750300120048,381.832625700486,0.311679481156615,0.00199600798403194,FALSE,TRUE,9447,500,"sig7_deg0_bin120_wavenumber" +8,2,120,"wavenumber",9.63025210084034,497.11296063639,0.363102113637611,0.00199600798403194,FALSE,TRUE,9447,500,"sig8_deg2_bin120_wavenumber" +8,1,120,"wavenumber",9.69387755102041,508.250627700088,0.364117658509001,0.00199600798403194,FALSE,TRUE,9447,500,"sig8_deg1_bin120_wavenumber" +8,0,120,"wavenumber",9.72569027611044,517.388433563542,0.361552396396774,0.00199600798403194,FALSE,TRUE,9447,500,"sig8_deg0_bin120_wavenumber" +4,0,200,"wavenumber",31.3583433373349,120.022205448983,0.16116250830519,0.00199600798403194,TRUE,TRUE,9447,500,"sig4_deg0_bin200_wavenumber" +4,1,200,"wavenumber",31.3583433373349,120.276693357547,0.161470076000754,0.00199600798403194,TRUE,TRUE,9447,500,"sig4_deg1_bin200_wavenumber" +4,2,200,"wavenumber",31.3583433373349,120.317384883968,0.161733798706019,0.00199600798403194,TRUE,TRUE,9447,500,"sig4_deg2_bin200_wavenumber" +4,0,160,"wavenumber",31.3583433373349,191.523785068614,0.203451054945692,0.00199600798403194,TRUE,TRUE,9447,500,"sig4_deg0_bin160_wavenumber" +4,1,160,"wavenumber",31.3583433373349,191.874169361293,0.20376319815596,0.00199600798403194,TRUE,TRUE,9447,500,"sig4_deg1_bin160_wavenumber" +4,2,160,"wavenumber",31.3583433373349,192.328355519125,0.204255547193865,0.00199600798403194,TRUE,TRUE,9447,500,"sig4_deg2_bin160_wavenumber" +5,0,200,"wavenumber",31.3265306122449,194.465004004777,0.205031374017348,0.00199600798403194,TRUE,TRUE,9447,500,"sig5_deg0_bin200_wavenumber" +5,1,200,"wavenumber",31.3265306122449,194.861685671457,0.205367807242324,0.00199600798403194,TRUE,TRUE,9447,500,"sig5_deg1_bin200_wavenumber" +5,2,200,"wavenumber",31.3265306122449,195.016893114809,0.205715675790948,0.00199600798403194,TRUE,TRUE,9447,500,"sig5_deg2_bin200_wavenumber" +6,0,200,"wavenumber",31.3265306122449,258.556795908811,0.236352819321805,0.00199600798403194,TRUE,TRUE,9447,500,"sig6_deg0_bin200_wavenumber" +6,1,200,"wavenumber",31.3265306122449,259.178756760527,0.236764872134194,0.00199600798403194,TRUE,TRUE,9447,500,"sig6_deg1_bin200_wavenumber" +6,2,200,"wavenumber",31.3265306122449,259.548360871709,0.237211782538629,0.00199600798403194,TRUE,TRUE,9447,500,"sig6_deg2_bin200_wavenumber" +5,0,160,"wavenumber",31.3265306122449,267.199946889302,0.2402462412111,0.00199600798403194,TRUE,TRUE,9447,500,"sig5_deg0_bin160_wavenumber" +5,1,160,"wavenumber",31.3265306122449,267.804479299904,0.2406364935376,0.00199600798403194,TRUE,TRUE,9447,500,"sig5_deg1_bin160_wavenumber" +5,2,160,"wavenumber",31.3265306122449,268.719880745123,0.241345459579706,0.00199600798403194,TRUE,TRUE,9447,500,"sig5_deg2_bin160_wavenumber" +4,0,120,"wavenumber",31.3265306122449,286.874999417522,0.249023789247506,0.00199600798403194,TRUE,TRUE,9447,500,"sig4_deg0_bin120_wavenumber" +4,1,120,"wavenumber",31.3583433373349,287.623695108705,0.249468191257671,0.00199600798403194,TRUE,TRUE,9447,500,"sig4_deg1_bin120_wavenumber" +4,2,120,"wavenumber",31.3265306122449,288.954398265615,0.250370654098415,0.00199600798403194,TRUE,TRUE,9447,500,"sig4_deg2_bin120_wavenumber" +7,0,200,"wavenumber",31.2947178871549,304.248319746702,0.256388063283114,0.00199600798403194,TRUE,TRUE,9447,500,"sig7_deg0_bin200_wavenumber" +7,1,200,"wavenumber",31.2947178871549,305.158000043472,0.256894848867292,0.00199600798403194,TRUE,TRUE,9447,500,"sig7_deg1_bin200_wavenumber" +7,2,200,"wavenumber",31.2947178871549,306.019884666066,0.257602568321779,0.00199600798403194,TRUE,TRUE,9447,500,"sig7_deg2_bin200_wavenumber" +6,0,160,"wavenumber",31.3265306122449,314.670991439908,0.260735397624097,0.00199600798403194,TRUE,TRUE,9447,500,"sig6_deg0_bin160_wavenumber" +6,1,160,"observed",31.3265306122449,315.376115676902,0.267872314470838,0.00199600798403194,TRUE,TRUE,8994,500,"sig6_deg1_bin160_observed" +6,1,160,"wavenumber",31.3265306122449,315.727660603729,0.261305834349196,0.00199600798403194,TRUE,TRUE,9447,500,"sig6_deg1_bin160_wavenumber" +6,2,160,"wavenumber",31.2947178871549,317.42904914475,0.262380813619649,0.00199600798403194,TRUE,TRUE,9447,500,"sig6_deg2_bin160_wavenumber" +8,0,200,"wavenumber",31.2947178871549,332.939018597347,0.268235860533237,0.00199600798403194,TRUE,TRUE,9447,500,"sig8_deg0_bin200_wavenumber" +5,0,120,"wavenumber",31.3265306122449,333.946669109476,0.268735879080348,0.00199600798403194,TRUE,TRUE,9447,500,"sig5_deg0_bin120_wavenumber" +8,1,200,"wavenumber",31.2947178871549,334.301739521917,0.268922475749762,0.00199600798403194,TRUE,TRUE,9447,500,"sig8_deg1_bin200_wavenumber" +5,1,120,"wavenumber",31.3265306122449,335.383912820845,0.269453735775903,0.00199600798403194,TRUE,TRUE,9447,500,"sig5_deg1_bin120_wavenumber" +8,2,200,"wavenumber",31.2947178871549,335.857373193256,0.269946298958297,0.00199600798403194,TRUE,TRUE,9447,500,"sig8_deg2_bin200_wavenumber" +5,2,120,"wavenumber",31.2947178871549,338.1520352881,0.271009489799855,0.00199600798403194,TRUE,TRUE,9447,500,"sig5_deg2_bin120_wavenumber" +7,0,160,"wavenumber",31.2947178871549,339.977141567488,0.27106431914111,0.00199600798403194,TRUE,TRUE,9447,500,"sig7_deg0_bin160_wavenumber" +7,1,160,"wavenumber",31.2947178871549,341.623637677083,0.271862261275608,0.00199600798403194,TRUE,TRUE,9447,500,"sig7_deg1_bin160_wavenumber" +7,2,160,"wavenumber",31.2947178871549,344.640594959049,0.273524902123414,0.00199600798403194,TRUE,TRUE,9447,500,"sig7_deg2_bin160_wavenumber" +6,1,160,"ritz",31.3265306122449,346.338782567562,0.268053804171389,0.00199600798403194,TRUE,TRUE,9858,500,"sig6_deg1_bin160_ritz" +8,0,160,"wavenumber",31.2947178871549,352.926421741162,0.276225631297536,0.00199600798403194,TRUE,TRUE,9447,500,"sig8_deg0_bin160_wavenumber" +6,0,120,"wavenumber",31.2947178871549,353.318474822445,0.276486668007853,0.00199600798403194,TRUE,TRUE,9447,500,"sig6_deg0_bin120_wavenumber" +8,1,160,"wavenumber",31.2947178871549,355.416151088066,0.277366980803298,0.00199600798403194,TRUE,TRUE,9447,500,"sig8_deg1_bin160_wavenumber" +6,1,120,"wavenumber",31.3265306122449,355.744301050131,0.27761532337233,0.00199600798403194,TRUE,TRUE,9447,500,"sig6_deg1_bin120_wavenumber" +8,2,160,"wavenumber",31.2629051620648,360.139096405204,0.279782691043789,0.00199600798403194,TRUE,TRUE,9447,500,"sig8_deg2_bin160_wavenumber" +6,2,120,"wavenumber",31.2947178871549,360.843774114349,0.28016270726631,0.00199600798403194,TRUE,TRUE,9447,500,"sig6_deg2_bin120_wavenumber" +7,2,120,"wavenumber",31.2629051620648,373.965760764014,0.285442245038734,0.00199600798403194,TRUE,TRUE,9447,500,"sig7_deg2_bin120_wavenumber" diff --git a/tests/testthat/test-final-adversarial-audit.R b/tests/testthat/test-final-adversarial-audit.R new file mode 100644 index 0000000..8f13ee4 --- /dev/null +++ b/tests/testthat/test-final-adversarial-audit.R @@ -0,0 +1,154 @@ +# test-final-adversarial-audit.R +# Fast unit tests for the final adversarial audit. These tests exercise pure +# helpers and small synthetic vectors only; they do not run Monte Carlo scans. + +.src_final_module <- function(name) { + p <- file.path(dirname(AUDIT_PATH), name) + e <- new.env(parent = globalenv()) + sys.source(p, envir = e) + e +} + +test_that("final audit modules source and CLI runner parses without executing", { + # Library-style modules are safe to source directly under testthat. + for (name in c("global_multiple_testing.R", "run_alternative_nulls.R", + "run_holdout_replication.R", "build_final_adversarial_summary.R", + "run_resolution_mode_diagnostics.R", "model_comparison.R")) { + e <- .src_final_module(name) + expect_true(exists("main", envir = e, inherits = FALSE), info = name) + } + + # The top-level CLI runner intentionally resolves repository paths at source + # time. Under testthat the working directory differs from normal CLI use, so + # syntax/parse safety is the appropriate no-side-effect CI check here. + runner <- file.path(dirname(AUDIT_PATH), "run_final_adversarial_audit.R") + expect_silent(parse(file = runner)) +}) + +test_that("negative-binomial size reduces to Poisson when no overdispersion is estimated", { + alt <- .src_final_module("run_alternative_nulls.R") + mu <- rep(20, 50) + y <- rep(20, 50) + expect_true(is.infinite(alt$estimate_nb_size(y, mu))) +}) + +test_that("negative-binomial size is finite under clear extra-Poisson dispersion", { + alt <- .src_final_module("run_alternative_nulls.R") + mu <- rep(20, 100) + y <- rep(c(2, 38), 50) + size <- alt$estimate_nb_size(y, mu) + expect_true(is.finite(size)) + expect_gt(size, 0) +}) + +test_that("block-residual adversarial mean remains positive and preserves expected total", { + alt <- .src_final_module("run_alternative_nulls.R") + setup_rng(123) + mu <- seq(10, 30, length.out = 40) + y <- round(mu + rep(c(-3, 2, 4, -2), 10)) + star <- alt$block_residual_mean(y, mu, block_len = 5L) + expect_length(star, length(mu)) + expect_true(all(is.finite(star))) + expect_true(all(star > 0)) + expect_equal(sum(star), sum(mu), tolerance = 1e-10) +}) + +test_that("alternative-null summary uses plus-one empirical p and never reports zero", { + alt <- .src_final_module("run_alternative_nulls.R") + maxima <- rep(10, 99) + s <- alt$summarize_null_model("x", maxima, observed_stat = 100, B = 99) + expect_equal(s$tail_count, 0L) + expect_equal(s$scan_global_p, 1 / 100) + expect_equal(s$resolution_floor, 1 / 100) + expect_gt(s$scan_global_p, 0) +}) + +test_that("alternative-null model registry contains the canonical and adversarial generators", { + alt <- .src_final_module("run_alternative_nulls.R") + models <- alt$model_metadata()$null_model + expect_setequal(models, c( + "poisson_fixed_baseline", + "poisson_refit_baseline", + "conditional_multinomial", + "negative_binomial_refit", + "block_residual_refit" + )) +}) + +test_that("strict holdout verdict cannot hide a partial failure", { + s <- .src_final_module("build_final_adversarial_summary.R") + expect_equal(s$holdout_verdict(c(0.01, 0.02), c(TRUE, TRUE)), "pass") + expect_equal(s$holdout_verdict(c(0.01, 0.20), c(TRUE, TRUE)), "mixed") + expect_equal(s$holdout_verdict(c(0.20, 0.30), c(TRUE, TRUE)), "fail") + expect_equal(s$holdout_verdict(c(0.01, 0.02), c(TRUE, FALSE)), "mixed") +}) + +test_that("descriptive robustness thresholds are explicit pass/fail rather than softened", { + s <- .src_final_module("build_final_adversarial_summary.R") + expect_equal(s$threshold_verdict(0.80, 0.80, TRUE), "pass") + expect_equal(s$threshold_verdict(0.70, 0.80, TRUE), "fail") + expect_equal(s$threshold_verdict(0.049, 0.05, FALSE), "pass") + expect_equal(s$threshold_verdict(0.051, 0.05, FALSE), "fail") +}) + +test_that("primary multiplicity row recognises canonical metadata", { + s <- .src_final_module("build_final_adversarial_summary.R") + x <- data.frame( + analysis_id = c("fe_ion2_wn_bin160_sig5_deg1", "fe_ion2_wn_bin160_sig6_deg1"), + species = c("Fe", "Fe"), source = c("wavenumber", "wavenumber"), + bins = c(160, 160), sigma = c(5, 6), degree = c(1, 1), + family_max_p = c(0.2, 0.01), stringsAsFactors = FALSE) + row <- s$primary_row(x) + expect_equal(row$analysis_id, "fe_ion2_wn_bin160_sig6_deg1") + expect_equal(row$family_max_p, 0.01) +}) + +test_that("baseline-refit holdout test requires a smoothing sigma", { + ho <- .src_final_module("run_holdout_replication.R") + expect_error( + ho$fixed_k_test_p(c(1, 2, 3), c(4, 5, 6), c(4, 5, 6), + degree = 1, k = 2, B = 2, refit_baseline = TRUE), + "sigma is required" + ) +}) + +test_that("held-out polynomial prediction reuses the training transform", { + mc <- .src_final_module("model_comparison.R") + tr <- c(1, 2, 3, 4) + te <- c(10, 11) + tf <- mc$poly_transform(tr) + + # The explicit frozen transform must reproduce the canonical design on train. + expect_equal(mc$design_poly_from_transform(tr, 2, tf), design_poly(tr, 2), + tolerance = 1e-12) + + # Test points are mapped with TRAINING center/scale, not their own statistics. + expected_z <- (te - mean(tr)) / sqrt(mean((tr - mean(tr))^2)) + x_te <- mc$design_poly_from_transform(te, 1, tf) + expect_equal(x_te[, 2], expected_z, tolerance = 1e-12) + expect_false(isTRUE(all.equal(x_te[, 2], design_poly(te, 1)[, 2]))) +}) + +test_that("matched smoothing holds approximate ell-space width fixed", { + d <- .src_final_module("run_resolution_mode_diagnostics.R") + expect_equal(d$matched_sigma(c(60, 80, 100, 160, 240)), + c(2.25, 3.00, 3.75, 6.00, 9.00)) +}) + +test_that("resolution diagnostic extracts distinct ranked local peaks", { + d <- .src_final_module("run_resolution_mode_diagnostics.R") + scan <- data.frame(k = 1:7, deltaD = c(0, 3, 1, 4, 1, 2, 0)) + peaks <- d$extract_top_peaks(scan, 3L) + expect_equal(peaks$k, c(4, 2, 6)) + expect_equal(peaks$peak_rank, 1:3) +}) + +test_that("resolution branch labels distinguish frozen external and Fe primary modes", { + d <- .src_final_module("run_resolution_mode_diagnostics.R") + expect_equal(d$classify_peak_branch(9.63, d$GWTC_K_FROZEN, 31.3265, 0.02), + "external_gwtc_9p602") + expect_equal(d$classify_peak_branch(31.30, d$GWTC_K_FROZEN, 31.3265, 0.02), + "primary_fe") + expect_equal(d$classify_peak_branch(20, d$GWTC_K_FROZEN, 31.3265, 0.02), + "other") +})