From cc65ad54f5cecb20ed45cc9f1122a08e5fed202e Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 19:07:01 -0500 Subject: [PATCH 01/10] Reject saturated and degenerate local BIC models --- dagguard.py | 97 ++++++++++++++++++++++++++++++++++++----------------- 1 file changed, 67 insertions(+), 30 deletions(-) diff --git a/dagguard.py b/dagguard.py index 0afad54..dad4a6f 100644 --- a/dagguard.py +++ b/dagguard.py @@ -8,7 +8,7 @@ - ``method='exact'``: exact child-wise best-subset search using enumeration and branch-and-bound. -The public API enforces the regular full-column-rank condition used by the +The public API enforces the regular local-regression conditions used by the Gaussian-BIC theory. The numerical engine remains in ``local_bic_refinement.py`` for backward compatibility with the original NOTEARS-BP reproducibility commit. New code should import from this module. @@ -54,18 +54,18 @@ def _validate_candidate_full_rank(X, candidate_adjacency, rank_tolerance=None): - """Validate the regular full-rank condition for every candidate regression. + """Validate regular Gaussian-BIC conditions for every candidate regression. The conventional local BIC used by DAGGuard penalizes the nominal number of selected parents. To keep that score aligned with regular Gaussian BIC, the - public API requires each child's full candidate-parent design to have full - column rank after centering. Every deletion subset is then also full rank. - Rank-deficient candidates should first remove redundant predictors or use a - score explicitly designed for singular models. + public API requires each child's centered full candidate-parent design to + have full column rank, to satisfy ``q_j < n - 1``, and to have strictly + positive, numerically nondegenerate full-model residual variance. Every + deletion subset is then also nonsaturated and full rank. The rank check uses the same scale-stabilization principle as the score - engine: centered nonconstant predictor columns are normalized before the SVD. - This makes the numerical validation invariant to changes of measurement + engine: centered nonconstant predictor columns are normalized before the + SVD. This makes the numerical validation invariant to changes of measurement units. When ``rank_tolerance`` is supplied, it is interpreted as a relative singular-value cutoff, matching NumPy least-squares ``rcond`` semantics. """ @@ -73,6 +73,8 @@ def _validate_candidate_full_rank(X, candidate_adjacency, rank_tolerance=None): A = (np.asarray(candidate_adjacency) != 0).astype(np.int8) if X.ndim != 2: raise ValueError("X must be a two-dimensional numeric array") + if X.shape[0] < 2: + raise ValueError("X must contain at least two observations") if A.ndim != 2 or A.shape[0] != A.shape[1] or A.shape[0] != X.shape[1]: raise ValueError("candidate_adjacency must be square and match X") if not np.all(np.isfinite(X)): @@ -81,35 +83,69 @@ def _validate_candidate_full_rank(X, candidate_adjacency, rank_tolerance=None): raise ValueError("rank_tolerance must be nonnegative") Xc = X - X.mean(axis=0, keepdims=True) + n = Xc.shape[0] eps = np.finfo(float).eps + residual_ratio_floor = 100.0 * eps + for child in range(A.shape[0]): + y = Xc[:, child] + tss = float(y @ y) + if tss == 0.0: + raise ValueError( + "candidate child response is constant after centering for child " + f"{child}; conventional Gaussian BIC requires positive residual " + "variance" + ) + parents = np.flatnonzero(A[:, child]) q = int(len(parents)) - if q == 0: - continue - design = Xc[:, parents] - norms = np.linalg.norm(design, axis=0) - if np.any(norms <= eps): + if q >= n - 1: raise ValueError( - "candidate parent design is rank deficient for child " - f"{child}: at least one candidate parent is constant after " - "centering; remove redundant predictors before Gaussian-BIC " - "refinement" + "candidate local regression is saturated for child " + f"{child}: q_j={q} candidate parents with n={n}; require " + "q_j < n - 1 for conventional Gaussian BIC" ) - stable_design = design / norms - singular = np.linalg.svd(stable_design, compute_uv=False) - if singular.size == 0: - rank = 0 - elif rank_tolerance is None: - cutoff = max(stable_design.shape) * eps * singular[0] - rank = int(np.sum(singular > cutoff)) + + if q == 0: + full_rss = tss else: - rank = int(np.sum(singular > float(rank_tolerance) * singular[0])) - if rank < q: + design = Xc[:, parents] + norms = np.linalg.norm(design, axis=0) + if np.any(norms <= eps): + raise ValueError( + "candidate parent design is rank deficient for child " + f"{child}: at least one candidate parent is constant after " + "centering; remove redundant predictors before Gaussian-BIC " + "refinement" + ) + stable_design = design / norms + singular = np.linalg.svd(stable_design, compute_uv=False) + if singular.size == 0: + rank = 0 + elif rank_tolerance is None: + cutoff = max(stable_design.shape) * eps * singular[0] + rank = int(np.sum(singular > cutoff)) + else: + rank = int(np.sum(singular > float(rank_tolerance) * singular[0])) + if rank < q: + raise ValueError( + "candidate parent design is rank deficient for child " + f"{child}: rank {rank} < {q} candidate parents; remove " + "redundant predictors before Gaussian-BIC refinement" + ) + + coef, *_ = np.linalg.lstsq( + stable_design, y, rcond=rank_tolerance + ) + residual = y - stable_design @ coef + full_rss = float(residual @ residual) + + if full_rss <= 0.0 or full_rss / tss <= residual_ratio_floor: raise ValueError( - "candidate parent design is rank deficient for child " - f"{child}: rank {rank} < {q} candidate parents; remove " - "redundant predictors before Gaussian-BIC refinement" + "candidate local regression is degenerate for child " + f"{child}: full-model residual variance is numerically zero; " + "conventional Gaussian BIC requires strictly positive residual " + "variance" ) @@ -133,7 +169,8 @@ def refine_dag( Directed candidate adjacency with ``A[parent, child] = 1``. The graph must be acyclic. DAGGuard only deletes edges; it never adds or reverses an edge. For conventional Gaussian BIC, each child's centered candidate - parent design must have full column rank. + parent design must have full column rank, ``q_j < n - 1``, and positive + full-model residual variance. method : {"exact", "greedy"} Exact best-subset refinement or fast greedy deletion. enumeration_max_parents, branch_node_limit : int From 72febd385c80dd8adbfff76c8feb962a12ebc9ca Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 19:07:21 -0500 Subject: [PATCH 02/10] Add saturated and degenerate model validation tests --- tests/test_dagguard_api.py | 48 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/tests/test_dagguard_api.py b/tests/test_dagguard_api.py index cf575c7..76779eb 100644 --- a/tests/test_dagguard_api.py +++ b/tests/test_dagguard_api.py @@ -103,6 +103,54 @@ def test_public_pruning_pressure_is_scale_invariant(self): [(r["parent"], r["child"], r["below_cutoff"]) for r in changed_rows], ) + def test_saturated_candidate_is_rejected(self): + rng = np.random.default_rng(31) + n = 5 + parents = rng.normal(size=(n, n - 1)) + y = rng.normal(size=n) + X = np.column_stack([parents, y]) + child = n - 1 + candidate = np.zeros((n, n), dtype=int) + candidate[: n - 1, child] = 1 + with self.assertRaisesRegex(ValueError, "saturated"): + refine_dag(X, candidate, method="exact") + + def test_constant_response_is_rejected(self): + rng = np.random.default_rng(37) + x0 = rng.normal(size=100) + y = np.ones(100) + X = np.column_stack([x0, y]) + candidate = np.zeros((2, 2), dtype=int) + candidate[0, 1] = 1 + with self.assertRaisesRegex(ValueError, "constant after centering"): + refine_dag(X, candidate, method="exact") + + def test_exact_fit_is_rejected(self): + rng = np.random.default_rng(41) + n = 120 + x0 = rng.normal(size=n) + x1 = rng.normal(size=n) + y = 2.0 * x0 - 0.5 * x1 + X = np.column_stack([x0, x1, y]) + candidate = np.zeros((3, 3), dtype=int) + candidate[0, 2] = 1 + candidate[1, 2] = 1 + with self.assertRaisesRegex(ValueError, "residual variance is numerically zero"): + refine_dag(X, candidate, method="exact") + + def test_near_zero_full_model_rss_is_rejected(self): + rng = np.random.default_rng(43) + n = 150 + x0 = rng.normal(size=n) + x1 = rng.normal(size=n) + y = 1.5 * x0 - 0.3 * x1 + 1e-9 * rng.normal(size=n) + X = np.column_stack([x0, x1, y]) + candidate = np.zeros((3, 3), dtype=int) + candidate[0, 2] = 1 + candidate[1, 2] = 1 + with self.assertRaisesRegex(ValueError, "residual variance is numerically zero"): + refine_dag(X, candidate, method="exact") + @staticmethod def _scale_test_problem(): rng = np.random.default_rng(29) From 5f9faae8c7c6db744c02d20b4b92b0e673cbf492 Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:14:45 -0500 Subject: [PATCH 03/10] Document external comparator source audit --- benchmarks/seven_method/SOURCE_AUDIT.md | 33 +++++++++++++++++++++++++ 1 file changed, 33 insertions(+) create mode 100644 benchmarks/seven_method/SOURCE_AUDIT.md diff --git a/benchmarks/seven_method/SOURCE_AUDIT.md b/benchmarks/seven_method/SOURCE_AUDIT.md new file mode 100644 index 0000000..1d72544 --- /dev/null +++ b/benchmarks/seven_method/SOURCE_AUDIT.md @@ -0,0 +1,33 @@ +# External-method source audit + +This benchmark distinguishes direct source ports from published-parameter adaptations. The purpose of this file is to make the external-method comparisons auditable without implying that independently implemented or adapted methods are official author software. + +## NOTEARS + +The project linear NOTEARS implementation was checked against the public `xunzheng/notears` linear least-squares implementation. The objective, positive/negative variable split for the L1 penalty, matrix-exponential acyclicity function and gradient, augmented-Lagrangian updates, zero-diagonal bounds, `h_tol=1e-8`, and `rho_max=1e16` agree in substance. The reported primary settings `lambda1=0.1` and coefficient threshold `0.3` match the standard public example. The project adds a deterministic safeguard that drops threshold-passing edges only if a residual numerical cycle would otherwise remain; this has no effect when the thresholded solution is already acyclic. + +## Li and Wang (2009) PC-FDR + +`pc_fdr.py` is an independent implementation of the paper's Algorithm 3 for the skeleton. The source audit checked the following features against Algorithm 3 and its Gaussian appendix: ordered-pair testing with current neighborhoods, maximum conditional-independence p-value (`pmax`) updates, repeated FDR application after valid p-values have been accumulated, theorem-version retention of `pmax` values for removed edges, and Fisher-z testing with the `sqrt(n-|C|-3)` factor. The primary benchmark uses the paper's practical `H*=H` step-up choice at `q=0.05`; `q=0.10` and a Benjamini-Yekutieli sensitivity are reported separately. Finite-sample FDP in these simulations is treated as an empirical operating characteristic, not as a test of the paper's asymptotic guarantee. + +## Strobl, Spirtes, and Visweswaran (2019) PC-p + +`pcp_faithful.py` is a source-audited Python port, not official author software. The audit used the authors' public MATLAB repository `ericstrobl/PCp`, master tree SHA `42a179d7305641dbe6f8926e46a67ba622c66524`, and checked the workflow corresponding to `PC_with_pval.m`, `get_skeleton_stable.m`, `get_v_structures2.m`, `clamp_edges.m`, `orientation_rules.m`, `control_FDR.m`, `binary_search.m`, `get_BY_FDR.m`, and `gaussCItest.m`. The port preserves the official adaptive initial PC threshold, PC-stable neighborhood snapshots, p-value propagation, and Benjamini-Yekutieli FDR step. + +One rare source-code ambiguity is documented rather than hidden: the MATLAB conflict-handling code contains expressions that can be read literally as linear indexing although the surrounding logic indicates edge-coordinate indexing. The primary Python port uses the coordinate interpretation and exposes `literal_conflict_indexing=True` for audit purposes. No orientation conflicts occurred in any of the 240 primary simulation runs, so this ambiguity cannot affect the primary simulation table. + +Earlier exploratory PC-p results from an incomplete port were discarded and are not part of the manuscript or archived primary results. + +## Wang et al. (2026) hybrid structural pipeline + +No official public implementation was located. The benchmark therefore labels this method as a transparent published-parameter adaptation rather than an exact reproduction. `wang_full.py`/`wang_sensitivity.py` implement the published structural Steps 1-3: mutual-information skeleton screening, conditional-mutual-information collider identification and pruning, and local discrete-BIC parent pruning. The primary thresholds are the paper's setting A `(0.008, 0.005, 0.009)`; the other three settings tied for best BIC in the source application are included as sensitivity analyses. + +The source application is discrete and includes domain-specific root-node restrictions. For the generic continuous simulation benchmark, variables are discretized by empirical tertiles and no source-application root labels are transferred. The paper does not fully formalize the phrase describing exclusion of collider-related nodes in the Step-3 conditioning set, so the implementation uses a conservative documented interpretation. Step 4 is not benchmarked because it only orients the retained skeleton and the common endpoint is skeleton adjacency; its state-wise orientation rule is also not sufficiently specified for a generic continuous-variable adaptation. These choices are limitations of comparability, not claimed features of the authors' original implementation. + +## Ordinary PC + +The ordinary-PC baseline uses an original-style ordered-pair skeleton search with immediate graph updates and two-sided Gaussian Fisher-z tests at `alpha=0.05`. It is presented as a conventional PC skeleton baseline, not as a reproduction of a particular software package. + +## Common endpoint and data checks + +PC-family methods and the adapted Wang procedure need not produce the same type of fully oriented DAG as NOTEARS/DAGGuard. The primary seven-method comparison therefore uses skeleton adjacency for all methods. Simulation comparisons use the same 240 `(d, s, noise, rep, seed)` keys. On the authorized commercial data, the comparator methods were rerun from the pinned input bytes and their adjacency matrices matched the archived results entry-for-entry. From dc1260edb4cd0258b42eab4254cecb743056a155 Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:16:31 -0500 Subject: [PATCH 04/10] Validate all public refinement entry points --- dagguard.py | 51 +++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 43 insertions(+), 8 deletions(-) diff --git a/dagguard.py b/dagguard.py index dad4a6f..34283c1 100644 --- a/dagguard.py +++ b/dagguard.py @@ -25,10 +25,10 @@ candidate_indegree_summary, deletion_diagnostics, edge_jaccard, - exact_refine_dag, + exact_refine_dag as _exact_refine_dag_engine, gaussian_local_bic, graph_metrics, - greedy_refine_dag, + greedy_refine_dag as _greedy_refine_dag_engine, initial_pruning_pressure, is_acyclic, total_gaussian_bic, @@ -134,9 +134,7 @@ def _validate_candidate_full_rank(X, candidate_adjacency, rank_tolerance=None): "redundant predictors before Gaussian-BIC refinement" ) - coef, *_ = np.linalg.lstsq( - stable_design, y, rcond=rank_tolerance - ) + coef, *_ = np.linalg.lstsq(stable_design, y, rcond=rank_tolerance) residual = y - stable_design @ coef full_rss = float(residual @ residual) @@ -149,6 +147,44 @@ def _validate_candidate_full_rank(X, candidate_adjacency, rank_tolerance=None): ) +def exact_refine_dag( + X, + candidate_adjacency, + *, + enumeration_max_parents: int = 15, + branch_node_limit: int = 2_000_000, + score_tolerance: float = 1e-10, + rank_tolerance: float | None = None, +) -> RefinementResult: + """Validated public wrapper for exact fixed-candidate refinement.""" + _validate_candidate_full_rank(X, candidate_adjacency, rank_tolerance) + return _exact_refine_dag_engine( + X, + candidate_adjacency, + enumeration_max_parents=enumeration_max_parents, + branch_node_limit=branch_node_limit, + score_tolerance=score_tolerance, + rank_tolerance=rank_tolerance, + ) + + +def greedy_refine_dag( + X, + candidate_adjacency, + *, + score_tolerance: float = 1e-10, + rank_tolerance: float | None = None, +) -> RefinementResult: + """Validated public wrapper for greedy fixed-candidate refinement.""" + _validate_candidate_full_rank(X, candidate_adjacency, rank_tolerance) + return _greedy_refine_dag_engine( + X, + candidate_adjacency, + score_tolerance=score_tolerance, + rank_tolerance=rank_tolerance, + ) + + def refine_dag( X, candidate_adjacency, @@ -190,7 +226,6 @@ def refine_dag( Selected adjacency, score, runtime, search diagnostics, and exact-search status where applicable. """ - _validate_candidate_full_rank(X, candidate_adjacency, rank_tolerance) method = str(method).lower() if method == "exact": return exact_refine_dag( @@ -221,9 +256,9 @@ def pruning_pressure(X, candidate_adjacency, *, rank_tolerance: float | None = N def dagguard_exact(X, candidate_adjacency, **kwargs) -> RefinementResult: """Convenience wrapper for exact DAGGuard refinement.""" - return refine_dag(X, candidate_adjacency, method="exact", **kwargs) + return exact_refine_dag(X, candidate_adjacency, **kwargs) def dagguard_greedy(X, candidate_adjacency, **kwargs) -> RefinementResult: """Convenience wrapper for fast DAGGuard-Greedy refinement.""" - return refine_dag(X, candidate_adjacency, method="greedy", **kwargs) + return greedy_refine_dag(X, candidate_adjacency, **kwargs) From 4fb0145c1e2bbd6196065c73e76e2f823ae60f1b Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:16:50 -0500 Subject: [PATCH 05/10] Add benchmark source-audit regression tests --- tests/test_benchmark_source_audits.py | 57 +++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 tests/test_benchmark_source_audits.py diff --git a/tests/test_benchmark_source_audits.py b/tests/test_benchmark_source_audits.py new file mode 100644 index 0000000..dce0616 --- /dev/null +++ b/tests/test_benchmark_source_audits.py @@ -0,0 +1,57 @@ +import unittest + +import numpy as np +from scipy.stats import norm + +from benchmarks.seven_method.common import fisherz_p +from benchmarks.seven_method.pc_fdr import fdr_stepup +from benchmarks.seven_method.pcp_faithful import by_fdr +from benchmarks.seven_method.wang_full import local_disc_bic +from dagguard import exact_refine_dag, greedy_refine_dag + + +class BenchmarkSourceAuditTests(unittest.TestCase): + def test_fisher_z_matches_published_gaussian_formula(self): + correlation = np.array([[1.0, 0.30], [0.30, 1.0]]) + n = 100 + observed = fisherz_p(correlation, n, 0, 1, ()) + statistic = np.sqrt(n - 3) * np.arctanh(0.30) + expected = 2.0 * norm.sf(abs(statistic)) + self.assertAlmostEqual(observed, expected, places=14) + + def test_pc_fdr_stepup_matches_bh_rule(self): + pvalues = np.array([0.001, 0.010, 0.200]) + rejected = fdr_stepup(pvalues, 0.05, by=False) + np.testing.assert_array_equal(rejected, [True, True, False]) + + def test_pc_p_by_estimator_matches_official_formula(self): + pvalues = np.array([0.01, 0.02, 0.20]) + alpha = 0.05 + harmonic = 1.0 + 1.0 / 2.0 + 1.0 / 3.0 + expected = len(pvalues) * alpha * harmonic / 2.0 + self.assertAlmostEqual(by_fdr(pvalues, alpha), expected, places=14) + + def test_wang_discrete_bic_parameter_count(self): + # A binary parent perfectly predicts a binary child. The saturated + # conditional log-likelihood is zero and there are two free conditional + # Bernoulli parameters, giving BIC = -0.5 * 2 * log(4). + D = np.array([[0, 0], [0, 0], [1, 1], [1, 1]], dtype=int) + observed = local_disc_bic(D, child=1, parents=(0,)) + expected = -np.log(4.0) + self.assertAlmostEqual(observed, expected, places=14) + + def test_all_public_refinement_entrypoints_reject_saturation(self): + rng = np.random.default_rng(20260824) + n = 5 + parents = rng.normal(size=(n, n - 1)) + y = rng.normal(size=n) + X = np.column_stack([parents, y]) + candidate = np.zeros((n, n), dtype=int) + candidate[: n - 1, n - 1] = 1 + for refiner in (exact_refine_dag, greedy_refine_dag): + with self.assertRaisesRegex(ValueError, "saturated"): + refiner(X, candidate) + + +if __name__ == "__main__": + unittest.main() From d98323cea227bfec0a77e5c72be64cb504cb77f8 Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:17:13 -0500 Subject: [PATCH 06/10] Align README with final validation and comparator audit --- README.md | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/README.md b/README.md index 415fd9d..7cabd75 100644 --- a/README.md +++ b/README.md @@ -33,11 +33,11 @@ delete edge iff partial R^2 < 1 - n^(-1/n). DAGGuard is a score-based refinement procedure, not a finite-sample nominal FDR method. Its formal guarantees are conditional on a fixed or independently learned candidate. The current end-to-end empirical validation is specifically for NOTEARS-generated candidates. -For conventional Gaussian BIC, the public API requires each child's centered full candidate-parent design to have full column rank; rank-deficient candidate designs are rejected. Before numerical rank assessment, nonconstant candidate-parent columns are normalized by their Euclidean norms, so the validation is stable to changes of measurement units. `globally_optimal=True` means the exact search established the minimum score within the documented numerical tolerance and did not hit its search limit; it does not imply a unique representative when multiple subsets are numerically tied. +For conventional Gaussian BIC, every public refinement entry point requires each child's centered full candidate-parent design to have full column rank, `q_j < n - 1`, and strictly positive, numerically nondegenerate full-model residual variance. Rank-deficient, saturated, constant-response, and degenerate candidate regressions are rejected with informative exceptions. Before numerical rank assessment, nonconstant candidate-parent columns are normalized by their Euclidean norms, so validation is stable to changes of measurement units. `globally_optimal=True` means the exact search established the minimum score within the documented numerical tolerance and did not hit its search limit; it does not imply a unique representative when multiple subsets are numerically tied. ## Installation -Python 3.12 was used for the reported analyses. +Python 3.12 was used for the reported analyses. Core package versions are pinned in `requirements.txt`. ```bash python -m pip install -r requirements.txt @@ -57,14 +57,14 @@ print(exact.globally_optimal) print(exact.total_bic) ``` +The lower-level names `exact_refine_dag` and `greedy_refine_dag` exported by `dagguard.py` use the same validation policy. `local_bic_refinement.py` remains the tested numerical engine for backward compatibility with the earlier reproducibility commit. + Run the minimal example: ```bash python -m examples.dagguard_quickstart ``` -`local_bic_refinement.py` remains the tested numerical engine for backward compatibility with the earlier reproducibility commit. - ## Main empirical results ### Controlled fixed-candidate experiments @@ -89,7 +89,7 @@ The proprietary row-level data are not distributed. The real-data workflow recor ## Repository map -- `dagguard.py` - public DAGGuard API. +- `dagguard.py` - validated public DAGGuard API. - `local_bic_refinement.py` - backward-compatible numerical engine: local BIC, exact search, greedy search, pruning pressure, and graph metrics. - `candidate_contamination_simulations.py` - 12-setting fixed-candidate experiment. - `reproduce_simulations.py` - primary NOTEARS simulation workflow. @@ -97,10 +97,11 @@ The proprietary row-level data are not distributed. The real-data workflow recor - `notears_tuning_sensitivity.py` - NOTEARS penalty/threshold sensitivity analysis. - `realdata_postselection_diagnostics.py` - authorized swine-data analysis. - `synthetic_application_twin.py` - public 37-variable workflow without proprietary observations. -- `reproduce_submission.sh` - staged reproduction entrypoint. +- `reproduce_submission.sh` - staged reproduction entry point. - `benchmarks/seven_method/` - audited competitor implementations and real-data benchmark runner. +- `benchmarks/seven_method/SOURCE_AUDIT.md` - source-by-source comparator audit and documented adaptation choices. - `results/seven_method_benchmark/` - audited benchmark summary tables (no proprietary observations). -- `tests/` - deterministic unit and numerical-policy tests. +- `tests/` - deterministic numerical, public-API, and benchmark source-audit tests. - `REAL_DATA_SCHEMA.md` - construction of the 37 application variables. ## Reproduce the main DAGGuard analyses @@ -128,9 +129,9 @@ python realdata_postselection_diagnostics.py \ ## Seven-method benchmark provenance -The benchmark distinguishes the methods' inferential targets. PC-family procedures are compared by skeleton adjacency because they need not return a uniquely oriented DAG. The Wang et al. (2026) structural method is transparently adapted to continuous variables by empirical tertiles and is not represented as official author software. PC-p is a source-audited Python port of the authors' official MATLAB code because MATLAB/Octave was unavailable in the benchmark runtime. +The benchmark distinguishes the methods' inferential targets and implementation status. NOTEARS was checked against the public `xunzheng/notears` linear implementation. Li & Wang PC-FDR is an independent implementation of the published Algorithm 3. PC-p is a source-audited Python port of the authors' official MATLAB code because MATLAB/Octave was unavailable in the benchmark runtime. The Wang et al. (2026) structural method is a published-parameter adaptation to continuous variables by empirical tertiles and is not represented as official author software. PC-family procedures and the Wang adaptation are compared by skeleton adjacency because they need not return a uniquely oriented DAG comparable to NOTEARS/DAGGuard. -See `benchmarks/seven_method/README.md` and `results/seven_method_benchmark/method_implementation_provenance.csv`. +See `benchmarks/seven_method/SOURCE_AUDIT.md`, `benchmarks/seven_method/README.md`, `results/seven_method_benchmark/method_implementation_provenance.csv`, and `results/seven_method_benchmark/AUDIT_NOTES.md`. ## Tests @@ -140,7 +141,7 @@ python -m examples.dagguard_quickstart python synthetic_application_twin.py --out results/synthetic_application_twin ``` -The public API tests include scale-invariance checks for exact refinement, greedy refinement, and pruning pressure, together with duplicate-column, near-collinearity, and near-tie cases. +The public API tests cover scale invariance for exact refinement, greedy refinement, and pruning pressure; duplicate columns; near-collinearity; near ties; saturated local models; constant responses; and exact or numerically near-exact fits. Additional regression tests check the Gaussian Fisher-z formula, PC-FDR step-up rule, PC-p BY estimator, and the discrete-BIC parameter count used in the Wang adaptation. ## Data availability From 6bf1c2b14151859cfe9aabae93150a2f17fe2ba3 Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:17:24 -0500 Subject: [PATCH 07/10] Clarify comparator implementation status and source audit --- benchmarks/seven_method/README.md | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/benchmarks/seven_method/README.md b/benchmarks/seven_method/README.md index c7e1611..e190089 100644 --- a/benchmarks/seven_method/README.md +++ b/benchmarks/seven_method/README.md @@ -1,20 +1,22 @@ # Seven-method benchmark -This directory contains the comparator implementations used for the DAGGuard JDS benchmark. The common simulation endpoint is skeleton adjacency because PC, PC-FDR, and PC-p do not necessarily return a uniquely oriented DAG. +This directory contains the comparator implementations used for the DAGGuard JDS benchmark. The common simulation endpoint is skeleton adjacency because PC, PC-FDR, PC-p, and the adapted Wang procedure do not necessarily return a uniquely oriented DAG comparable to NOTEARS/DAGGuard. Primary methods: -1. NOTEARS (existing project implementation) +1. NOTEARS (project implementation checked against the public `xunzheng/notears` linear source) 2. NOTEARS + DAGGuard-Greedy 3. NOTEARS + DAGGuard-Exact 4. Wang et al. (2026) hybrid structural pipeline, transparently adapted to continuous variables by empirical-tertile discretization -5. Li & Wang (2009) PC-FDR -6. Strobl, Spirtes & Visweswaran (2019) PC-p +5. Li & Wang (2009) PC-FDR, independently implemented from published Algorithm 3 +6. Strobl, Spirtes & Visweswaran (2019) PC-p, source-audited Python port of the authors' MATLAB code 7. ordinary PC -`pcp_faithful.py` is a Python port audited against the authors' official MATLAB source (`ericstrobl/PCp`, master tree SHA `42a179d7305641dbe6f8926e46a67ba622c66524`). It is labeled as a port, not official author software. +`SOURCE_AUDIT.md` records the source-by-source audit, including which components are direct source ports, which are independent implementations, and which are adaptations. This distinction is intentional: none of the independently implemented or adapted comparators is represented as official author software. -The Wang et al. implementation reproduces structural Steps 1-3 with the published threshold setting A `(0.008, 0.005, 0.009)` and empirical-tertile discretization for continuous inputs. Step 4 only orients the resulting skeleton and was not used for the common skeleton benchmark because the publication does not fully specify a generic multilevel adaptation of that orientation rule. All four threshold settings tied as optimal in their application were checked as sensitivity analyses. +`pcp_faithful.py` was audited against the authors' official MATLAB repository (`ericstrobl/PCp`, master tree SHA `42a179d7305641dbe6f8926e46a67ba622c66524`). The port preserves the official adaptive PC threshold, PC-stable neighborhood snapshots, p-value propagation, and BY FDR procedure. A rare conflict-indexing ambiguity in the MATLAB source is documented in `SOURCE_AUDIT.md`; no orientation conflicts occurred in the 240 primary simulation runs, so it cannot affect the primary benchmark table. + +The Wang et al. implementation reproduces the published structural Steps 1-3 as a generic continuous-data adaptation using empirical-tertile discretization. The primary threshold setting is A `(0.008, 0.005, 0.009)` and all four source-paper settings tied as optimal in their application are checked as sensitivities. Source-application root-node labels are not transferred to the generic simulation benchmark. Step 4 only orients the retained skeleton and was not used for the common skeleton endpoint because the publication does not fully specify a generic multilevel adaptation of that orientation rule. For the proprietary swine data, run only with an authorized local copy: @@ -26,4 +28,4 @@ python -m benchmarks.seven_method.swine_benchmark \ The script writes only non-row-level adjacency matrices and summaries. The expected SHA256 for the analysis file used in the manuscript is `b933fd66f49fd381bb9698ee2b3f5835d0db8d01820a01d5e45c9ac3a7bf5156`. -See `results/seven_method_benchmark/` for the audited primary summary tables included with the repository. +See `SOURCE_AUDIT.md` and `results/seven_method_benchmark/` for the audit record and primary summary tables. From b1489bc531db955a36ce0d0dee404031eb76b08f Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:17:31 -0500 Subject: [PATCH 08/10] Refine comparator provenance wording --- .../method_implementation_provenance.csv | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/results/seven_method_benchmark/method_implementation_provenance.csv b/results/seven_method_benchmark/method_implementation_provenance.csv index 19dd562..136fdce 100644 --- a/results/seven_method_benchmark/method_implementation_provenance.csv +++ b/results/seven_method_benchmark/method_implementation_provenance.csv @@ -1,8 +1,8 @@ method,implementation,software_source,primary_parameters,output -NOTEARS,Pinned project implementation,historical pre-rename commit 509cb29c24967d12b99e2b53641349ce7bb470ed in ChongWangStat/DAGGuard,lambda1=0.1; edge threshold=0.3,DAG -NOTEARS + DAGGuard-Greedy,Pinned project implementation,same historical pinned commit; local_bic_refinement.py,Gaussian local BIC; repeated best single deletion,DAG deletion subgraph -NOTEARS + DAGGuard-Exact,Pinned project implementation,same historical pinned commit; exact enumeration/branch-and-bound,Gaussian local BIC; certified exact fixed-candidate optimum,DAG deletion subgraph -Wang et al. (2026) hybrid structural pipeline,Transparent published-parameter adaptation of structural Steps 1-3; no official public source code located. Step 4 orientation not benchmarked because it does not alter the primary skeleton endpoint and is under-specified for generic multilevel adaptation.,BMC Med Inform Decis Mak 26:257 (2026); paper states pgmpy and bnlearn used,"empirical tertile discretization; published A: eps_skeleton=.008, eps_collider=.005, eps_prune=.009; A-D sensitivity",final skeleton after structural Steps 1-3; orientation not compared -Li & Wang (2009) PC-FDR,Independent implementation of published Algorithm 3,JMLR 10:475-514 (2009),Fisher-z; q=.05 primary; practical H*=H; q=.10 and BY sensitivity,skeleton -"Strobl, Spirtes & Visweswaran (2019) PC-p",Faithful Python port audited line-by-line against official MATLAB source,ericstrobl/PCp master tree 42a179d7305641dbe6f8926e46a67ba622c66524; ACM TIST 10(5):46,Fisher-z; q=.05 primary; official adaptive alpha and BY FDR; q=.10 sensitivity,CPDAG adjacency/skeleton for common comparison -Ordinary PC,Tetrad-style ordered-pair skeleton search,PC skeleton logic aligned to Li & Wang Algorithm 1; Fisher-z implementation in benchmark code,alpha=.05,skeleton +NOTEARS,"Pinned project implementation checked against the public linear NOTEARS source","historical pre-rename commit 509cb29c24967d12b99e2b53641349ce7bb470ed in ChongWangStat/DAGGuard; source audit against xunzheng/notears linear.py","lambda1=0.1; edge threshold=0.3",DAG +NOTEARS + DAGGuard-Greedy,Pinned project implementation,"same historical pinned commit; local_bic_refinement.py","Gaussian local BIC; repeated best single deletion",DAG deletion subgraph +NOTEARS + DAGGuard-Exact,Pinned project implementation,"same historical pinned commit; exact enumeration/branch-and-bound","Gaussian local BIC; certified exact fixed-candidate optimum",DAG deletion subgraph +Wang et al. (2026) hybrid structural pipeline,"Transparent published-parameter adaptation of structural Steps 1-3; no official public source code located. Step 4 orientation not benchmarked because it does not alter the primary skeleton endpoint and is under-specified for generic multilevel adaptation.","BMC Med Inform Decis Mak 26:257 (2026); paper states pgmpy and bnlearn used","empirical tertile discretization; published A: eps_skeleton=.008, eps_collider=.005, eps_prune=.009; A-D sensitivity",final skeleton after structural Steps 1-3; orientation not compared +Li & Wang (2009) PC-FDR,Independent implementation of published Algorithm 3,"JMLR 10:475-514 (2009); audited against Algorithm 3 and Gaussian appendix","Fisher-z; q=.05 primary; practical H*=H; q=.10 and BY sensitivity",skeleton +"Strobl, Spirtes & Visweswaran (2019) PC-p","Source-audited Python port of official MATLAB source","ericstrobl/PCp master tree 42a179d7305641dbe6f8926e46a67ba622c66524; ACM TIST 10(5):46","Fisher-z; q=.05 primary; official adaptive alpha and BY FDR; q=.10 sensitivity",CPDAG adjacency/skeleton for common comparison +Ordinary PC,"Conventional ordered-pair skeleton search","PC skeleton logic aligned to Li & Wang Algorithm 1; Fisher-z implementation in benchmark code",alpha=.05,skeleton From 38a30b12f2844d7c15164a440f3c4ff11c3a6383 Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:17:42 -0500 Subject: [PATCH 09/10] Expand benchmark source and reproducibility audit --- results/seven_method_benchmark/AUDIT_NOTES.md | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/results/seven_method_benchmark/AUDIT_NOTES.md b/results/seven_method_benchmark/AUDIT_NOTES.md index 773dcbd..c70f86f 100644 --- a/results/seven_method_benchmark/AUDIT_NOTES.md +++ b/results/seven_method_benchmark/AUDIT_NOTES.md @@ -1,11 +1,13 @@ # Benchmark audit notes -- Uploaded `train2023cw_simple.csv` SHA256: `b933fd66f49fd381bb9698ee2b3f5835d0db8d01820a01d5e45c9ac3a7bf5156`. This exactly matches the pinned application input at NOTEARS-BP commit `509cb29c24967d12b99e2b53641349ce7bb470ed`. +- Uploaded `train2023cw_simple.csv` SHA256: `b933fd66f49fd381bb9698ee2b3f5835d0db8d01820a01d5e45c9ac3a7bf5156`. This exactly matches the pinned application input at historical commit `509cb29c24967d12b99e2b53641349ce7bb470ed`. - Independent preprocessing reproduces 2,592 source rows, 36 incomplete rows excluded, and the authoritative 2,556 x 37 analysis matrix. -- Ordinary PC, PC-FDR q=.05/q=.10/BY, PC-p q=.05/q=.10, and Wang published settings A-D were rerun on the newly uploaded data; every adjacency matrix matched the stored benchmark result entry-for-entry. -- PC, PC-FDR, PC-p and Wang simulation outputs use the same 240 `(d,s,noise,rep,seed)` keys. See the reproducibility archive for the common-key check. -- PC-FDR was re-audited against Li & Wang (2009) Algorithm 3 and Appendix B after its finite-sample FDP exceeded q. The implementation uses ordered pairs, current neighborhoods, pmax updates, theorem-version retention of removed-edge pmax values, and Fisher-z with `sqrt(n-|C|-3)*atanh(r)`. The observed FDP is therefore reported as a regime-specific finite-sample benchmark, not as a contradiction of the asymptotic theorem. -- PC-p uses a Python port audited against the authors' official MATLAB repository `ericstrobl/PCp`, master tree SHA `42a179d7305641dbe6f8926e46a67ba622c66524`. MATLAB/Octave was unavailable. Earlier exploratory PC-p port results were found to be too sparse because the port was incomplete; they were discarded and are not included in this package. -- Wang et al. (2026) has no official public source repository located. The benchmark implements published structural Steps 1-3, adapts continuous variables by empirical tertiles, uses published threshold setting A primarily, and includes all four tied application-optimal threshold settings as sensitivity analyses. Step 4 only orients the resulting skeleton and does not affect the common skeleton endpoint; it is not claimed as reproduced for generic multilevel variables because the paper does not fully specify that reduction. The thresholds were not re-tuned on the present datasets. +- Ordinary PC, PC-FDR q=.05/q=.10/BY, PC-p q=.05/q=.10, and Wang published settings A-D were rerun on the authorized data; every comparator adjacency matrix matched the stored benchmark result entry-for-entry. +- PC, PC-FDR, PC-p, and Wang simulation outputs use the same 240 `(d,s,noise,rep,seed)` keys. See the reproducibility archive for the common-key check. +- The project NOTEARS implementation was checked against the public `xunzheng/notears` linear least-squares source. The objective, positive/negative L1 parameterization, matrix-exponential acyclicity function and gradient, augmented-Lagrangian updates, and standard numerical tolerances agree in substance. The project additionally uses a deterministic safeguard against residual numerical cycles after coefficient thresholding; this is inert when the thresholded solution is already acyclic. +- PC-FDR was re-audited against Li & Wang (2009) Algorithm 3 and its Gaussian testing appendix after its finite-sample FDP exceeded q. The implementation uses ordered pairs, current neighborhoods, `pmax` updates, theorem-version retention of removed-edge `pmax` values, and Fisher-z with `sqrt(n-|C|-3)*atanh(r)`. The observed FDP is therefore reported as a regime-specific finite-sample benchmark, not as a contradiction of the paper's asymptotic theorem. +- PC-p uses a Python port audited against the authors' official MATLAB repository `ericstrobl/PCp`, master tree SHA `42a179d7305641dbe6f8926e46a67ba622c66524`. MATLAB/Octave was unavailable. Earlier exploratory PC-p port results were found to be too sparse because the port was incomplete; they were discarded and are not included in this package. A rare indexing ambiguity in the official orientation-conflict code is documented in `benchmarks/seven_method/SOURCE_AUDIT.md`; no orientation conflicts occurred in the 240 primary simulation runs. +- Wang et al. (2026) has no official public source repository located. The benchmark is explicitly a published-parameter adaptation of structural Steps 1-3, not official author software. It adapts continuous variables by empirical tertiles, uses published threshold setting A primarily, and includes all four source-paper threshold settings tied as optimal in their application as sensitivities. Domain-specific source-application root labels are not transferred to the generic simulation benchmark. Step 4 only orients the retained skeleton and does not affect the common skeleton endpoint; it is not claimed as reproduced for generic multilevel variables because the paper does not fully specify that reduction. +- A detailed source-by-source audit is maintained in `benchmarks/seven_method/SOURCE_AUDIT.md`. Deterministic tests check the Gaussian Fisher-z formula, the PC-FDR step-up rule, the PC-p BY estimator, and the discrete-BIC parameter count used in the Wang adaptation. - An independent real-data NOTEARS refit was attempted in the benchmark runtime but exceeded a 10-minute execution window without generating output. No incomplete result was used. The authoritative candidate/greedy/exact results are tied to the exact data bytes by the matching SHA256; the existing certified exact search took 1,915 seconds and need not be repeated. - No row-level proprietary data are included in this deliverable. From 742c5a6430cbfd0a11920838bfe50ec9109a3999 Mon Sep 17 00:00:00 2001 From: Chong Wang <87325417+ChongWangStat@users.noreply.github.com> Date: Mon, 24 Aug 2026 20:17:50 -0500 Subject: [PATCH 10/10] Update v1.0.0 release candidate notes --- RELEASE_NOTES.md | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index 9707b67..04e2c07 100644 --- a/RELEASE_NOTES.md +++ b/RELEASE_NOTES.md @@ -4,18 +4,21 @@ This repository state accompanies the JDS submission **“DAGGuard: Exact and Gr ## Reproducibility scope -- Public refinement API: `dagguard.py`. +- Validated public refinement API: `dagguard.py`. - Historical numerical engine and pinned application audit trail: `local_bic_refinement.py` and commit `509cb29c24967d12b99e2b53641349ce7bb470ed`. - Controlled candidate-contamination simulations: 1,200 replicates. - NOTEARS end-to-end Gaussian/non-Gaussian simulations: 240 datasets. - NOTEARS penalty/threshold sensitivity addressing dependence on upstream tuning. - Seven-method benchmark provenance and non-row-level summaries: `results/seven_method_benchmark/`. +- External-method source audit: `benchmarks/seven_method/SOURCE_AUDIT.md`. - Proprietary swine observations are not distributed; the authorized workflow verifies the pinned input SHA-256. - `synthetic_application_twin.py` provides a public 37-variable end-to-end example without confidential observations. -- `reproduce_submission.sh` is the staged reproduction entrypoint. +- `reproduce_submission.sh` is the staged reproduction entry point. ## Numerical policy -The public API requires each centered full candidate-parent design to have full column rank for conventional Gaussian BIC. Rank-deficient candidates are rejected. `globally_optimal=True` means that the exact search established the minimum objective value within the documented numerical tolerance and did not hit its search limit; it does not imply uniqueness among numerically tied subsets. +For conventional Gaussian BIC, every public refinement entry point requires each centered full candidate-parent design to have full column rank, `q_j < n - 1`, and strictly positive, numerically nondegenerate full-model residual variance. Rank-deficient, saturated, constant-response, and degenerate candidate regressions are rejected rather than silently scored. `globally_optimal=True` means that the exact search established the minimum objective value within the documented numerical tolerance and did not hit its search limit; it does not imply uniqueness among numerically tied subsets. -A GitHub tag/release should point to the final merged submission commit after the referee-revision branch is approved. +The deterministic test suite covers the public numerical policy, exact-versus-greedy invariants, branch-and-bound certification, and targeted source-audit checks for the comparator calculations. + +A GitHub tag/release should point to the final merged submission commit after the final acceptance-focused pull request passes CI and is merged.