Bayesian Network Propensity Score (BNPS)

• BNPS is a method that uses Bayesian Networks to nonparametrically estimate propensity scores in observational studies.
• It efficiently models complex dependency structures among covariates, maintaining consistency, asymptotic normality, and efficiency.
• Paired with the Hajek estimator, BNPS enhances inference accuracy and outperforms traditional logistic regression and machine-learning methods.

Bayesian Network Propensity Score (BNPS) denotes a propensity-score methodology in which the treatment assignment mechanism is estimated through Bayesian Networks rather than through a rigid parametric model such as logistic regression. In this setting, a Bayesian Network is a directed acyclic graph encoding the joint distribution of treatment and covariates, and the propensity score is recovered as the conditional probability of treatment given the observed covariates. BNPS is motivated by observational studies with unknown (and likely unbalanced) designs and complex dependency structures among covariates, where conventional propensity-score models may be misspecified. Recent work presents BNPS as a flexible yet interpretable alternative that can be combined with inverse-probability estimators, especially the Hajek estimator, while preserving statistical properties such as consistency, asymptotic normality and asymptotic efficiency in the settings studied (Vicard et al., 2023, Serio et al., 3 Sep 2025).

1. Conceptual basis and motivation

BNPS arises from a familiar causal-inference problem: treatment assignment in observational data is rarely governed by a simple linear logit with additivity and low-order interactions. Standard propensity-score estimation by logistic regression can therefore be fragile in the presence of strong conditional dependencies, higher-order interactions, or hierarchical relationships among covariates. The BNPS literature frames this as a misspecification problem: rigid parametric assumptions can lead to biased estimates and inefficiencies when adjusting for confounding, thereby compromising causal inference (Serio et al., 3 Sep 2025).

The proposed remedy is to estimate the propensity score through a Bayesian Network (BN). In this usage, the BN is a graphical model over treatment and pre-treatment covariates. Its directed acyclic graph represents conditional independence structure, and its factorization induces a nonparametric estimate of the treatment assignment mechanism. This allows flexible modeling of conditional dependencies among covariates and between covariates and treatment, while retaining a transparent graphical representation of the assignment process (Vicard et al., 2023).

A central comparison in the literature is with contemporary machine-learning alternatives such as boosted CART, random forests, and Stable Balancing Weights. These methods are described as attractive in a predictive perspective, but the BNPS papers emphasize that they often lack theoretical guarantees for consistency, asymptotic normality, efficiency, and valid variance estimation. BNPS is positioned as an attempt to combine flexibility with those inferential guarantees (Serio et al., 3 Sep 2025).

2. Statistical construction

The propensity score remains the usual conditional probability

$e(\mathbf X)=P(T=1\mid \mathbf X).$

What changes in BNPS is the model used to estimate it. Instead of specifying a parametric regression for $T\mid \mathbf X$, BNPS models the joint distribution of treatment and covariates through a Bayesian Network: $$P(T,\mathbf X)=\prod_{j=1}^{p} P(X_j \mid \mathrm{Pa}(X_j)) \cdot P(T \mid \mathrm{Pa}(T)),$$ where $\mathrm{Pa}(\cdot)$ denotes parent nodes in the directed acyclic graph (Serio et al., 3 Sep 2025).

For the discrete setting developed explicitly in the testing paper, if $X_0=T$ and $(X_1,\ldots,X_L)$ are covariates, the fitted BN yields

$$\widehat{p}_k(\mathbf x) = \frac{\widehat p(k,x_1,\ldots,x_L)} {\widehat p(1,x_1,\ldots,x_L)+\widehat p(0,x_1,\ldots,x_L)}, \quad k\in\{0,1\},$$

which is the estimated propensity score for treatment level $k$ at covariate profile $\mathbf x$ (Vicard et al., 2023).

In the implementation described for BNPS, structure learning is performed from the observed treatment and covariates using Tabu search under the BIC criterion, as implemented in the R bnlearn package, and parameter learning proceeds by maximum-likelihood estimation of the conditional probability tables (Serio et al., 3 Sep 2025).

Component	BNPS specification	Output
Structure learning	Tabu search with BIC criterion	Learned DAG for treatment and covariates
Parameter learning	Maximum likelihood for conditional probability tables	Fitted joint distribution
Propensity score extraction	$P(T=1\mid \mathbf X)$ computed from the fitted BN	Unit-level propensity scores
Effect estimation	Inverse-probability estimators, especially Hajek	ATE or treatment-group means

Once $T\mid \mathbf X$0 is available, BNPS is paired with inverse-probability weighting. Two estimators recur in the literature. For binary outcomes and treatment-group mean $T\mid \mathbf X$1, the Horvitz–Thompson estimator is

$T\mid \mathbf X$2

and the Hajek estimator is

$T\mid \mathbf X$3

The average treatment effect is then estimated by

$T\mid \mathbf X$4

A more general Hajek ATE expression is also used in the 2025 BNPS paper: $T\mid \mathbf X$5 (Vicard et al., 2023, Serio et al., 3 Sep 2025).

3. Identification assumptions and asymptotic claims

BNPS relies on the standard causal assumptions used by propensity-score methods. The relevant papers state strong ignorability/unconfoundedness, meaning no unmeasured confounding given observed covariates, and positivity/common support, meaning every subject has nonzero probability of receiving each treatment. In the discrete-theory formulation, assumptions are presented as: covariates and potential outcomes are discrete/binary, treatment assignment is independent of potential outcomes given covariates, and propensity scores are bounded away from 0 and 1 (Vicard et al., 2023).

The theoretical claims attached to BNPS are unusually strong for a flexible propensity-score estimator. For discrete covariates, the BN maximum-likelihood estimator of the propensity score is stated to be consistent, asymptotically normal, and efficient, with asymptotic variance attaining the minimum possible variance among regular, consistent estimators. The 2025 BNPS paper generalizes the claim by stating that BNPS estimators are consistent and asymptotically normal when BN structure learning and parameter estimation are performed via maximum likelihood, and that BNPS achieves asymptotic efficiency when coupled with the Hajek estimator (Vicard et al., 2023, Serio et al., 3 Sep 2025).

The literature also distinguishes the estimators used after propensity-score estimation. Both Horvitz–Thompson and Hajek estimators are stated to be consistent and asymptotically normal under correct specification of the propensity score. Under misspecification, however, the Hajek estimator is described as less sensitive: its bias is smaller and is bounded by the maximal conditional outcome variation, whereas the Horvitz–Thompson estimator can exhibit larger bias. This difference is one reason the later BNPS work emphasizes BNPS + Hajek rather than BNPS + Horvitz–Thompson (Vicard et al., 2023).

A further claim concerns model specification. The 2025 paper states that BNPS mitigates the risk of misspecification by directly learning the dependency structure from data. This does not eliminate the need for the identification assumptions, but it shifts the burden from choosing a fixed functional form to learning a graph and conditional probability structure from observed data (Serio et al., 3 Sep 2025).

4. Position within Bayesian propensity-score research

BNPS belongs to a broader Bayesian propensity-score literature, but it addresses a different problem from many earlier Bayesian proposals. A foundational issue in that literature is the tension between the propensity score as a scalar balancing score and Bayes theorem as a requirement for a full generative model. Specifically, the literature shows that

$T\mid \mathbf X$6

so a naively joint Bayesian analysis based only on treatment, outcome, and the scalar propensity score can induce feedback, contaminate the propensity-score model with outcome information, and destroy the balancing interpretation (Zigler, 2013).

The response advocated in that literature is quasi-Bayesian, cut feedback, or modularized inference: estimate the propensity-score model using only treatment and covariates, then propagate uncertainty into the causal analysis without allowing the outcome to update the treatment model. This places much of Bayesian propensity-score methodology closer to a two-stage design-and-analysis workflow than to a single fully joint posterior (Zigler, 2013).

This perspective helps clarify BNPS. In its standard form, BNPS estimates the treatment assignment mechanism through a Bayesian Network over treatment and covariates and then uses inverse-probability estimators such as Hajek or Horvitz–Thompson for the outcome stage. This suggests that BNPS, as currently formulated, sits closer to the modular or design-based branch of Bayesian propensity-score analysis than to a full joint outcome-treatment model. Related Bayesian procedures pursue different strategies: posterior predictive treatment assignment randomizes inclusion into an overlap subset to address limited overlap, and Bayesian propensity score matching averages matches over posterior draws of the propensity score to reduce arbitrary dropping of units (Zigler et al., 2017, Alvarez et al., 2021).

5. Empirical evidence

The empirical case for BNPS rests on both simulation studies and an applied analysis in prostate cancer data. The 2025 paper reports extensive simulations across fifteen realistic scenarios and varying sample sizes, with $T\mid \mathbf X$7 and 1000 Monte Carlo datasets per scenario/sample size. The data-generating mechanism used 6 categorical $T\mid \mathbf X$8, treatment $T\mid \mathbf X$9, and binary outcome $$P(T,\mathbf X)=\prod_{j=1}^{p} P(X_j \mid \mathrm{Pa}(X_j)) \cdot P(T \mid \mathrm{Pa}(T)),$$0, with all variables simulated from a known BN structure. BNPS + Hajek was compared with logistic regression PS, Boosted CART PS, Random Forest PS, and SBW. Performance was assessed through Empirical Rejection Rate (ERR) and Empirical Coverage (EC) (Serio et al., 3 Sep 2025).

The reported findings are sharp. BNPS maintained correct type I error (ERR near 0.05 when ATE = 0), achieved high power (ERR $$P(T,\mathbf X)=\prod_{j=1}^{p} P(X_j \mid \mathrm{Pa}(X_j)) \cdot P(T \mid \mathrm{Pa}(T)),$$1 as $$P(T,\mathbf X)=\prod_{j=1}^{p} P(X_j \mid \mathrm{Pa}(X_j)) \cdot P(T \mid \mathrm{Pa}(T)),$$2 increases) when the ATE was nonzero, and produced coverage reliably near 95% across all scenarios and sample sizes. The same paper states that BNPS consistently outperformed benchmark methods in both empirical rejection rates and coverage accuracy. Machine-learning methods are reported to over-reject under no effect, especially as $$P(T,\mathbf X)=\prod_{j=1}^{p} P(X_j \mid \mathrm{Pa}(X_j)) \cdot P(T \mid \mathrm{Pa}(T)),$$3 increased, while SBW and wrongly specified logistic regressions showed increasing bias as sample size increased (Serio et al., 3 Sep 2025).

The earlier testing paper provides a more focused comparison between BN-based propensity scores estimated under AIC or BIC structure learning and two inverse-probability estimators. In that study, the Hajek estimator with BN-BIC PS estimation yielded coverage closest to nominal 95%, while the Hajek estimator was more robust than Horvitz–Thompson in small samples and under misspecification. The authors conclude that estimating the propensity score using BNs, especially with BIC structure, and using the Hajek estimator provided the best bias-variance tradeoff and the best inference properties in the simulated setting (Vicard et al., 2023).

The real-data application uses 7,162 prostate cancer patients from San Raffaele Hospital (Milan, Italy). The treatment is pelvic lymph node dissection during radical prostatectomy. The primary outcome is prolonged hospital stay (binary: >10 days) and the secondary outcome is 5-year biochemical recurrence (BCR). The learned BN is described as a complex web of dependencies among covariates and treatment, illustrating the graphical interpretability of the method (Serio et al., 3 Sep 2025).

Outcome	BNPS ATE	95% CI
Long hospital stay	0.13	[0.03, 0.23]
5-year biochemical recurrence	0.04	[-0.96, 1.03]

For hospital stay, the paper interprets the result as indicating that pelvic lymph node removal increases the risk of prolonged hospital stay, consistent with clinical knowledge. For biochemical recurrence, the unadjusted association was significant by Fisher test ($$P(T,\mathbf X)=\prod_{j=1}^{p} P(X_j \mid \mathrm{Pa}(X_j)) \cdot P(T \mid \mathrm{Pa}(T)),$$4), but after BNPS adjustment the paper reports no evidence for treatment effect after accounting for confounding (Serio et al., 3 Sep 2025).

6. Terminology, adjacent meanings, and limitations

The phrase Bayesian Network Propensity Score can be misread in two distinct ways. In the main BNPS literature, “Bayesian Network” refers to a directed acyclic graphical model used to estimate the treatment assignment mechanism. It does not refer to spillovers, interference, or treatment entanglement on social or biological networks. Those are separate research programs: one uses Bayesian generalized propensity scores for direct and spillover effects under neighborhood interference, and another models network evolution directly to handle entangled treatments (Forastiere et al., 2018, Toulis et al., 2018).

The acronym is also used differently in at least one paper on historical control borrowing. There, BNPS denotes a procedure that integrates propensity score adjustment, inverse probability weighting, Bayesian bootstrap, empirical Bayes, and a dynamic power prior for borrowing external controls. In that workflow, the propensity score is estimated by logistic regression for $$P(T,\mathbf X)=\prod_{j=1}^{p} P(X_j \mid \mathrm{Pa}(X_j)) \cdot P(T \mid \mathrm{Pa}(T)),$$5, not by a Bayesian Network DAG. The method is therefore distinct from BNPS in the Bayesian-Network sense, even though the acronym is shared (Wang et al., 2022).

The current BNPS evidence base also has clear scope conditions. The core theoretical development is for binary treatment, binary outcome, and discrete or categorical covariates. The real-data and simulation designs in the cited BNPS papers follow that template. The method still relies on unconfoundedness and positivity, and the general Bayesian propensity-score literature remains relevant because it highlights the broader danger of feedback-prone joint models and the importance of modularization (Vicard et al., 2023, Zigler, 2013).

The 2025 BNPS paper identifies continuous treatments, survival outcomes, and longitudinal data as promising avenues for future work. This suggests that the main unresolved question is not whether Bayesian Networks can define a propensity score, but how far the current asymptotic and finite-sample guarantees can be extended beyond discrete-treatment, discrete-covariate, inverse-probability-weighted settings (Serio et al., 3 Sep 2025).