Neyman–Pearson Causal Discovery

• Neyman–Pearson causal discovery is a framework that uses controlled error rates to determine causal links through hypothesis testing.
• It integrates methods like asymmetric edge-error control, calibrated conditional independence testing, and paired directional testing to enhance support recovery and causal direction estimation.
• The approach emphasizes maintaining level-α validity in inference, balancing false positives and negatives while addressing computational challenges in complex models.

Searching arXiv for the cited papers and topic context. Neyman–Pearson causal discovery denotes a family of causal inference frameworks that import the logic of classical Neyman–Pearson testing into causal graph learning and causal direction estimation. Across recent work, the central idea is to treat causal discovery as an error-controlled decision problem in which one type of mistake is constrained while another is minimized, or in which competing causal explanations are compared through level-$\alpha$ tests, directional power, and calibrated rejection behavior. In contemporary formulations, this perspective appears in at least three distinct but related forms: asymmetric edge-error control for graph support recovery, calibrated conditional independence testing inside constraint-based discovery with incomplete data, and paired directional testing for bivariate causal direction estimation (Shaska et al., 29 Jul 2025, Robinson et al., 6 May 2026, Prakash et al., 13 May 2026).

1. Conceptual definition and scope

Neyman–Pearson causal discovery is most explicitly defined in work on unequal edge error tolerance, where the goal is to minimize false negatives subject to a user-specified upper bound on false positives. In that formulation, the causal discovery problem is cast as the analogue of classical Neyman–Pearson testing: among all graph estimators whose false positive rate is at most $\epsilon$, select one with the smallest false negative rate (Shaska et al., 29 Jul 2025).

A broader usage of the term includes methods that do not claim a classical Neyman–Pearson lemma or uniformly most powerful rule, but nonetheless formulate causal discovery through level-$\alpha$ hypothesis tests, directional rejection probabilities, and calibrated inferential modules. Two examples are especially relevant. First, PAIR-CI treats conditional independence testing as the inferential engine of PC-style causal discovery and emphasizes that a CI oracle must behave like a level-$\alpha$ hypothesis test under incomplete data; its contribution is described explicitly through a Neyman–Pearson lens (Robinson et al., 6 May 2026). Second, CDSP formulates bivariate causal direction detection as a pair of hypothesis tests, defines a causal analogue of power and effect size, and bases direction choice on which directional null is more detectable (Prakash et al., 13 May 2026).

This suggests that “Neyman–Pearson causal discovery” is not a single algorithmic doctrine. Rather, it is a technical orientation in which causal discovery is organized around controlled error rates, explicit null and alternative hypotheses, and the asymmetry between different mistakes. A plausible implication is that the unifying concern is not merely graph recovery accuracy, but whether the inferential decisions embedded in a discovery procedure preserve the meaning of a nominal significance threshold.

2. Asymmetric edge-error control in graph support recovery

The most direct formulation appears in “Causal Link Discovery with Unequal Edge Error Tolerance” (Shaska et al., 29 Jul 2025). The paper studies a DAG

$$\mathcal{G}=(\mathcal{V},\mathcal{E}), \qquad |\mathcal{V}|=d,$$

with weighted adjacency matrix

$\boldsymbol{A}\in\mathbb{R}^{d\times d},$

and defines the support matrix $\boldsymbol{\chi}$ by

$$\chi_{i,j}=1 \iff \boldsymbol{A}_{i,j}\neq 0 \ \text{or}\ \boldsymbol{A}_{j,i}\neq 0.$$

The target is therefore support / causal link discovery, not orientation at first (Shaska et al., 29 Jul 2025).

The aggregate edgewise error rates are defined as

$\mbox{false positive rate}: \qquad \epsilon^+ = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=1,\ \boldsymbol{\chi}_{i,j}=0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}=0\}},$

and

$\mbox{false negative rate}: \qquad \epsilon^- = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=0,\ \boldsymbol{\chi}_{i,j}\neq 0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}\neq 0\}}.$

These are explicitly described as aggregated expected edgewise error rates, normalized by the expected numbers of absent and present edges, rather than family-wise error probabilities or FDR/FDP (Shaska et al., 29 Jul 2025).

The central optimization problem is

$\epsilon$0

where $\epsilon$1 (Shaska et al., 29 Jul 2025). This is the paper’s formal definition of Neyman–Pearson causal discovery. The constrained quantity is the false positive rate, and the optimized quantity is the false negative rate. The motivation is application-dependent asymmetry: in some domains, declaring a nonexistent causal link is more harmful than omitting a real one.

For each pair $\epsilon$2, the edgewise detection problem is written as

$\epsilon$3

$\epsilon$4

The paper then proves an analogue of the Neyman–Pearson lemma. The optimal detector is a weighted likelihood-ratio threshold rule involving $\epsilon$5, prior-induced weights $\epsilon$6, a global threshold $\epsilon$7, and boundary randomization by $\epsilon$8, with $\epsilon$9 selected so that $\alpha$0 (Shaska et al., 29 Jul 2025). The presence of randomization is explicitly noted as parallel to classical NP testing.

This is a genuine decision-theoretic formulation, not merely an interpretive analogy. It differs from most score-based or constraint-based graph learning procedures because the false-positive tolerance is fixed in advance and the detector is derived to satisfy that tolerance.

3. Information-theoretic characterization and the $\alpha$1-CUT algorithm

A major theoretical feature of the same framework is its characterization of performance limits through Rényi divergence (Shaska et al., 29 Jul 2025). For the optimal detector and any $\alpha$2, the paper gives an upper bound on $\alpha$3 in terms of $\alpha$4, the threshold $\alpha$5, and $\alpha$6. It also provides a converse lower bound for any detector satisfying $\alpha$7, again expressed through $\alpha$8, $\alpha$9, and $\alpha$0 (Shaska et al., 29 Jul 2025). The stated interpretation is that the tradeoff between constrained false positives and achievable false negatives is fundamentally governed by likelihood-ratio interpolation between the edge-absent and edge-present distributions.

Because the exact optimal detector depends on conditional edge distributions that are mixtures over graph structures and edge weights, it is computationally intractable in general. This motivates the algorithm $\alpha$1-CUT, designed for linear additive Gaussian noise models with equal-variance Gaussian noise: $\alpha$2 where $\alpha$3 are i.i.d. zero-mean Gaussian with covariance $\alpha$4, and $\alpha$5 is assumed known (Shaska et al., 29 Jul 2025).

The method exploits the fact that, after regressing a node on its true parents, the residual sum of squares has a chi-squared distribution. If $\alpha$6 is the least-squares estimate using the true parent set $\alpha$7, then

$\alpha$8

satisfies

$\alpha$9

conditioned on the parent observations, where $$\mathcal{G}=(\mathcal{V},\mathcal{E}), \qquad |\mathcal{V}|=d,$$0 (Shaska et al., 29 Jul 2025).

The resulting procedure tests every unordered pair $$\mathcal{G}=(\mathcal{V},\mathcal{E}), \qquad |\mathcal{V}|=d,$$1 by comparing residual energies under candidate parent sets and declares $$\mathcal{G}=(\mathcal{V},\mathcal{E}), \qquad |\mathcal{V}|=d,$$2 if any candidate pair of parent sets yields sufficiently close residual sums of squares. Otherwise it declares $$\mathcal{G}=(\mathcal{V},\mathcal{E}), \qquad |\mathcal{V}|=d,$$3 (Shaska et al., 29 Jul 2025). The threshold is selected via a chi-squared cdf equation so that the false-positive tolerance equals $$\mathcal{G}=(\mathcal{V},\mathcal{E}), \qquad |\mathcal{V}|=d,$$4.

The main finite-sample guarantee is:

• Theorem 3: for any sample size $$\mathcal{G}=(\mathcal{V},\mathcal{E}), \qquad |\mathcal{V}|=d,$$5 and any false positive tolerance $$\mathcal{G}=(\mathcal{V},\mathcal{E}), \qquad |\mathcal{V}|=d,$$6, the false positive rate of $$\mathcal{G}=(\mathcal{V},\mathcal{E}), \qquad |\mathcal{V}|=d,$$7-CUT satisfies

$$\mathcal{G}=(\mathcal{V},\mathcal{E}), \qquad |\mathcal{V}|=d,$$8

The paper further states a stronger per-edge result: if there is truly no edge between $$\mathcal{G}=(\mathcal{V},\mathcal{E}), \qquad |\mathcal{V}|=d,$$9 and $\boldsymbol{A}\in\mathbb{R}^{d\times d},$0, then

$\boldsymbol{A}\in\mathbb{R}^{d\times d},$1

Thus the method provides true finite-sample false-positive control for support recovery in the specified linear Gaussian SEM regime (Shaska et al., 29 Jul 2025).

The framework’s scope is correspondingly narrow. The guarantees are for support rather than orientation; direction estimation is deferred to a post-processing stage in which an unrestricted maximum-likelihood DAG estimate is masked by the support estimate: $\boldsymbol{A}\in\mathbb{R}^{d\times d},$2 The paper also notes exponential complexity in $\boldsymbol{A}\in\mathbb{R}^{d\times d},$3, with at most

$\boldsymbol{A}\in\mathbb{R}^{d\times d},$4

linear regressions and overall complexity

$\boldsymbol{A}\in\mathbb{R}^{d\times d},$5

(Shaska et al., 29 Jul 2025).

4. Neyman–Pearson logic in conditional independence testing for causal graphs

A second line of work applies Neyman–Pearson reasoning not to direct edge detection, but to the local conditional independence tests that drive constraint-based causal discovery. “PAIR-CI: Calibrated Conditional Independence Testing for Causal Discovery with Incomplete Data” develops a CI test specifically for constraint-based causal discovery with incomplete data and frames its contribution through a Neyman–Pearson lens: the objective is to provide a CI oracle that behaves like a level-$\boldsymbol{A}\in\mathbb{R}^{d\times d},$6 test when standard incomplete-data workflows do not (Robinson et al., 6 May 2026).

The paper formalizes the scientific null as

$\boldsymbol{A}\in\mathbb{R}^{d\times d},$7

with alternative

$\boldsymbol{A}\in\mathbb{R}^{d\times d},$8

Its strongest motivation is Proposition 1 (Miscalibration of impute-then-test), which states that if the imputation model is asymptotically misspecified and imputation error induces spurious dependence under the imputed distribution, then any consistent CI test satisfies

$\boldsymbol{A}\in\mathbb{R}^{d\times d},$9

The point is not that the CI test is defective on complete data, but that the imputation step changes the effective null (Robinson et al., 6 May 2026). This is an explicitly Neyman–Pearson-style failure: the nominal significance level no longer refers to the intended scientific null.

PAIR-CI addresses this by integrating multiple imputation into the inferential construction. It generates

$\boldsymbol{\chi}$0

then, on each imputed dataset and fold, compares a full model for predicting $\boldsymbol{\chi}$1 using $\boldsymbol{\chi}$2 with a partial model using the same imputed conditioning set and a placebo variable $\boldsymbol{\chi}$3 produced by conditional permutation. The fold-level loss contrast is

$\boldsymbol{\chi}$4

and in empirical form

$\boldsymbol{\chi}$5

Aggregation gives

$\boldsymbol{\chi}$6

The conceptual heart is the paired design: both models use the same imputed $\boldsymbol{\chi}$7, so imputation error is shared and is intended to cancel in the loss difference (Robinson et al., 6 May 2026).

The paper distinguishes the internal null

$\boldsymbol{\chi}$8

from the scientific null

$\boldsymbol{\chi}$9

Conditional permutation gives calibration for the internal null, and imputation consistency extends this to the scientific null (Robinson et al., 6 May 2026).

For inference, the combined variance is

$$\chi_{i,j}=1 \iff \boldsymbol{A}_{i,j}\neq 0 \ \text{or}\ \boldsymbol{A}_{j,i}\neq 0.$$0

where

$$\chi_{i,j}=1 \iff \boldsymbol{A}_{i,j}\neq 0 \ \text{or}\ \boldsymbol{A}_{j,i}\neq 0.$$1

The studentized statistic is

$$\chi_{i,j}=1 \iff \boldsymbol{A}_{i,j}\neq 0 \ \text{or}\ \boldsymbol{A}_{j,i}\neq 0.$$2

compared against a $$\chi_{i,j}=1 \iff \boldsymbol{A}_{i,j}\neq 0 \ \text{or}\ \boldsymbol{A}_{j,i}\neq 0.$$3-distribution with Barnard–Rubin adjusted degrees of freedom $$\chi_{i,j}=1 \iff \boldsymbol{A}_{i,j}\neq 0 \ \text{or}\ \boldsymbol{A}_{j,i}\neq 0.$$4, using the one-sided rejection region

$$\chi_{i,j}=1 \iff \boldsymbol{A}_{i,j}\neq 0 \ \text{or}\ \boldsymbol{A}_{j,i}\neq 0.$$5

The paper proves:

• Proposition 2 (Calibration):

$$\chi_{i,j}=1 \iff \boldsymbol{A}_{i,j}\neq 0 \ \text{or}\ \boldsymbol{A}_{j,i}\neq 0.$$6

• Proposition 3 (Consistency):

$$\chi_{i,j}=1 \iff \boldsymbol{A}_{i,j}\neq 0 \ \text{or}\ \boldsymbol{A}_{j,i}\neq 0.$$7

• Theorem 4 (Unified inference under cross-validation and multiple imputation):

$$\chi_{i,j}=1 \iff \boldsymbol{A}_{i,j}\neq 0 \ \text{or}\ \boldsymbol{A}_{j,i}\neq 0.$$8

• Corollary 5 (PC consistency):

$$\chi_{i,j}=1 \iff \boldsymbol{A}_{i,j}\neq 0 \ \text{or}\ \boldsymbol{A}_{j,i}\neq 0.$$9

The paper’s empirical emphasis is on null calibration. In standalone CI experiments with about 30% missingness, the calibration table reports average false positive rates

$\mbox{false positive rate}: \qquad \epsilon^+ = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=1,\ \boldsymbol{\chi}_{i,j}=0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}=0\}},$0

compared with MNAR false positive rates of $\mbox{false positive rate}: \qquad \epsilon^+ = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=1,\ \boldsymbol{\chi}_{i,j}=0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}=0\}},$1 for FZ-single, $\mbox{false positive rate}: \qquad \epsilon^+ = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=1,\ \boldsymbol{\chi}_{i,j}=0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}=0\}},$2 for FZ-Rubin, $\mbox{false positive rate}: \qquad \epsilon^+ = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=1,\ \boldsymbol{\chi}_{i,j}=0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}=0\}},$3 for imputed GCM, and $\mbox{false positive rate}: \qquad \epsilon^+ = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=1,\ \boldsymbol{\chi}_{i,j}=0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}=0\}},$4 for imputed KCI (Robinson et al., 6 May 2026). The abstract summarizes this as existing imputation-based CI tests exhibiting false positive rates of 28–45% under MNAR, whereas PAIR-CI averages below the nominal 5% level. When used inside PC, the method reduces structural Hamming distance by 8% on 10-variable nonlinear graphs, 15% on 30-variable nonlinear graphs, and up to 44% on the 56-variable HAILFINDER network (Robinson et al., 6 May 2026).

This line of work therefore locates Neyman–Pearson causal discovery at the level of the CI oracle: what matters is not only predictive performance, but whether the local test actually retains its level-$\mbox{false positive rate}: \qquad \epsilon^+ = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=1,\ \boldsymbol{\chi}_{i,j}=0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}=0\}},$5 meaning under the scientific null.

5. Directional power and effect-size asymmetry in bivariate discovery

A third formulation appears in “Causal Discovery via Statistical Power (CDSP)”, which addresses bivariate observational causal direction inference (Prakash et al., 13 May 2026). The paper does not claim classical NP optimality, but explicitly says that it formulates causal direction detection as a pair of hypothesis tests and introduces directional analogues of power, effect size, and directional evidence.

Under the standing assumption that the true direction is $\mbox{false positive rate}: \qquad \epsilon^+ = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=1,\ \boldsymbol{\chi}_{i,j}=0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}=0\}},$6, the data-generating mechanism is

$\mbox{false positive rate}: \qquad \epsilon^+ = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=1,\ \boldsymbol{\chi}_{i,j}=0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}=0\}},$7

with $\mbox{false positive rate}: \qquad \epsilon^+ = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=1,\ \boldsymbol{\chi}_{i,j}=0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}=0\}},$8. The reverse-direction working model is

$\mbox{false positive rate}: \qquad \epsilon^+ = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=1,\ \boldsymbol{\chi}_{i,j}=0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}=0\}},$9

where $\mbox{false negative rate}: \qquad \epsilon^- = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=0,\ \boldsymbol{\chi}_{i,j}\neq 0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}\neq 0\}}.$0 is the risk minimizer over a model class $\mbox{false negative rate}: \qquad \epsilon^- = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=0,\ \boldsymbol{\chi}_{i,j}\neq 0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}\neq 0\}}.$1, and $\mbox{false negative rate}: \qquad \epsilon^- = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=0,\ \boldsymbol{\chi}_{i,j}\neq 0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}\neq 0\}}.$2 (Prakash et al., 13 May 2026).

The directional nulls are

$\mbox{false negative rate}: \qquad \epsilon^- = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=0,\ \boldsymbol{\chi}_{i,j}\neq 0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}\neq 0\}}.$3

Thus each direction is evaluated by whether a valid structural model in that direction is compatible with the data (Prakash et al., 13 May 2026).

The paper defines Directional Power as

$\mbox{false negative rate}: \qquad \epsilon^- = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=0,\ \boldsymbol{\chi}_{i,j}\neq 0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}\neq 0\}}.$4

which is the joint success probability of retaining the correct-direction null while rejecting the reverse-direction null (Prakash et al., 13 May 2026).

Its central quantitative objects are the directional effect sizes

$\mbox{false negative rate}: \qquad \epsilon^- = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=0,\ \boldsymbol{\chi}_{i,j}\neq 0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}\neq 0\}}.$5

and the standardized critical values

$\mbox{false negative rate}: \qquad \epsilon^- = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=0,\ \boldsymbol{\chi}_{i,j}\neq 0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}\neq 0\}}.$6

The key assumption is the effect-size asymmetry assumption

$\mbox{false negative rate}: \qquad \epsilon^- = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=0,\ \boldsymbol{\chi}_{i,j}\neq 0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}\neq 0\}}.$7

Operationally, CDSP uses the directional detectability indices

$\mbox{false negative rate}: \qquad \epsilon^- = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=0,\ \boldsymbol{\chi}_{i,j}\neq 0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}\neq 0\}}.$8

and infers $\mbox{false negative rate}: \qquad \epsilon^- = \frac{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\hat{\boldsymbol{\chi}}_{i,j}=0,\ \boldsymbol{\chi}_{i,j}\neq 0\}}{\mathbb{E}\sum_{i,j}\mathbbm{1}\{\boldsymbol{\chi}_{i,j}\neq 0\}}.$9 when $\epsilon$00 (Prakash et al., 13 May 2026).

The main theorem states that under a specified misspecification regime,

$\epsilon$01

if and only if

$\epsilon$02

equivalently,

$\epsilon$03

The paper further proves directional consistency: $\epsilon$04 when $\epsilon$05 and the plug-in quantities are consistent (Prakash et al., 13 May 2026).

Uncertainty quantification is not given as a p-value or confidence set for the final direction. Instead the paper defines a bootstrap directional support probability

$\epsilon$06

for the direction $\epsilon$07, with an analogous expression for $\epsilon$08. The paper interprets this as the probability that the CDSP-favored direction remains favored under resampling (Prakash et al., 13 May 2026).

Empirically, the paper studies model misspecification in

$\epsilon$09

with $\epsilon$10. Reported simulation accuracies are $\epsilon$11 for $\epsilon$12, $\epsilon$13 for $\epsilon$14, and $\epsilon$15 for $\epsilon$16, while LiNGAM drops to $\epsilon$17 by $\epsilon$18 in the listed settings (Prakash et al., 13 May 2026). On the Tübingen Cause–Effect Pairs benchmark, CDSP achieves

$\epsilon$19

whereas LiNGAM yields $\epsilon$20 TDR and $\epsilon$21 FDR, corresponding to an approximately 18 percentage point reduction in FDR (Prakash et al., 13 May 2026).

CDSP therefore illustrates a distinct Neyman–Pearson-style idea: causal direction can be chosen by comparing the rejection behavior of two composite directional nulls, not by a single likelihood-ratio test between causal directions.

6. Related precursors and adjacent testing-based approaches

Earlier and adjacent work shows that the Neyman–Pearson perspective did not emerge in isolation. “Causal Discovery from Changes” studies learning from localized mechanism changes across a sequence of environments (Tian et al., 2013). It is not framed as a Neyman–Pearson method, but its pipeline depends on repeated binary decisions of the form “change” versus “no-change,” based on tests of

$\epsilon$22

with explicit discussion of significance level $\epsilon$23, Type I error, and Type II error (Tian et al., 2013). The method then translates those outcomes into descendant and nondescendant constraints under influentiality and local mechanism-change assumptions. This suggests an important precursor pattern: local hypothesis tests can be endowed with causal semantics and assembled into structural decisions, even if no explicit NP optimality result is proved.

A different adjacent approach is “Causal generalized linear models via Pearson risk invariance”, which identifies the causal parent set of a target variable $\epsilon$24 using a hybrid of invariance testing and model selection (Polinelli et al., 2024). The target obeys

$\epsilon$25

and the causal model is characterized by a Pearson risk identity

$\epsilon$26

together with expected-likelihood maximization on the true parent set (Polinelli et al., 2024). Finite-sample implementation uses the Pearson statistic

$\epsilon$27

with approximate null distribution

$\epsilon$28

and a two-sided acceptance region (Polinelli et al., 2024). This is not presented as Neyman–Pearson causal discovery, but it shows a closely related architecture: candidate causal models are screened by a calibrated null, then selected by BIC among the accepted models.

These adjacent cases clarify a common misconception. Neyman–Pearson causal discovery is not restricted to one specific graph-learning algorithm, nor does it require likelihood-ratio optimality in the classical simple-vs-simple sense. The recurring feature is the elevation of statistical testing, calibration, and asymmetric decision costs to first-class causal discovery principles.

7. Limitations, controversies, and technical boundaries

Several limitations recur across the literature.

The framework in (Shaska et al., 29 Jul 2025) is the most classically Neyman–Pearson in form, but its tractable algorithm $\epsilon$29-CUT is specialized to linear equal-variance Gaussian SEMs with known $\epsilon$30, and the finite-sample guarantee concerns support / skeleton recovery, not orientation. The optimal detector itself is generally computationally infeasible because the conditional edge distributions are mixtures over graphs and edge weights. The practical method also has exponential complexity in the number of variables (Shaska et al., 29 Jul 2025).

PAIR-CI (Robinson et al., 6 May 2026) addresses a different problem: preserving level-$\epsilon$31 behavior of CI tests under incomplete data. Its scientific-null calibration theorem formally relies on imputation consistency (A1), even though the paper argues that the paired construction is empirically robust beyond strict (A1). The authors are explicit that robustness under MNAR is practical rather than fully covered by the theorem. The method is also computationally heavier than parametric CI tests; on HAILFINDER, about 2400 tests take around 100 minutes per replicate. The paper does not provide a multiplicity treatment for the many sequential CI tests inside PC (Robinson et al., 6 May 2026).

CDSP (Prakash et al., 13 May 2026) likewise stops short of classical NP optimality. It provides no uniformly most powerful result, no global finite-sample Type I error control for the final direction decision, and no p-value for the selected causal direction. Its uncertainty measure is a bootstrap support probability rather than a confidence set. The effect-size asymmetry assumption is central but cannot be checked directly from observed data. Severe misspecification can reverse the asymmetry and induce confidently wrong decisions (Prakash et al., 13 May 2026).

A broader controversy concerns what exactly should be controlled in causal discovery. The explicit NP framework in (Shaska et al., 29 Jul 2025) constrains aggregated edgewise false positives. PAIR-CI emphasizes local test calibration because, in PC-like procedures, false CI rejections propagate through orientation rules and can be more damaging than false non-rejections (Robinson et al., 6 May 2026). CDSP instead compares directional success probabilities for paired nulls (Prakash et al., 13 May 2026). These are different inferential targets. This suggests that “error control” in causal discovery is not monolithic: support recovery, CI oracle behavior, and directional selection each instantiate distinct NP-style objects.

A second misconception is that better predictive performance is sufficient for reliable causal discovery. PAIR-CI explicitly argues that calibration matters at least as much as raw predictive performance, because miscalibrated CI tests create spurious edges that then propagate through orientation rules (Robinson et al., 6 May 2026). CDSP similarly shifts attention from raw dependence asymmetry to threshold-adjusted detectability indices (Prakash et al., 13 May 2026). The common implication is that causal discovery can fail even when its predictive subroutines are individually strong, if their rejection behavior is misaligned with the scientific nulls.

In summary, Neyman–Pearson causal discovery is best understood as a family of testing-centered formulations of causal inference in which causal claims are governed by constrained false-positive tolerance, level-$\epsilon$32 local null calibration, or directional power comparisons. Its most direct expression is the constrained optimization

$\epsilon$33

for support recovery (Shaska et al., 29 Jul 2025). Its broader significance lies in showing that causal discovery is often limited less by lack of predictive flexibility than by the lack of inferential procedures whose nulls, alternatives, and error rates remain meaningful under the data conditions at hand (Robinson et al., 6 May 2026, Prakash et al., 13 May 2026).