Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
45,506 characters · 14 sections · 20 citation commands
Causal EpiNets: Precision-corrected Bounds on Individual Treatment Effects using Epistemic Neural Networks
Estimating causal effects at the level of individuals is central in domains such as healthcare, policy design, and legal decision-making. While modern machine learning methods have made substantial progress in estimating conditional average treatment effects (CATEs), average effects do not directly answer the question that often motivates individual-level decisions: what is the probability that treatment would help this specific individual?
Average effects mask heterogeneity in individual responses. An average effect of zero, for instance, is consistent with both universal null effects and a mixture of strong positive and negative responders. For many decisions, the relevant quantity is not an average over similar individuals, but the probability that treatment would change the outcome for this individual.
The Probability of Necessity and Sufficiency (PNS), introduced by tian2000probabilities, provides a principled characterization of this notion. In the binary treatment and outcome setting, PNS is the probability that an individual would experience the outcome under treatment but not under control. A key result of tian2000probabilities shows that PNS is not point identified, but can be sharply bounded by combining experimental and observational data. Experimental data identify interventional response probabilities, while observational data constrain the joint distribution of treatment and outcome. Together, these sources restrict the set of admissible structural causal models. A profound consequence of this mathematical structure is that it allows for the computation of bounds on the PNS without requiring strong, untestable assumptions such as unconfoundedness or positivity—assumptions that otherwise dominate the causal inference literature. These identification results, however, assume exact knowledge of the relevant probabilities. In practice, the bounds must be estimated from finite samples, often with high-dimensional covariates. A natural approach is to estimate each probability from data and substitute these estimates into the bound expressions. Such naive plug-in estimators, however, suffer from two fundamental pathologies. First, they yield probability estimates that are incompatible with any underlying structural causal model. Second, PNS bounds are defined through maxima and minima over multiple estimated terms. This intersection-bound structure induces systematic extremum bias, causing lower bounds to be overestimated and upper bounds to be underestimated chernozhukov2013intersection. Valid inference therefore requires procedures that explicitly account for the envelope structure of the bounds.
We present the first complete framework for finite-sample estimation and inference of Tian–Pearl bounds on the Probability of Necessity and Sufficiency (PNS) with high-dimensional covariates. Our approach integrates three key components. First, we introduce an anchored, constrained multi-head neural architecture that jointly models observational and interventional quantities while guaranteeing compatibility with probability axioms and a shared structural causal model by construction. Second, we show how PNS bounds fit naturally into the intersection-bounds framework and adapt precision-corrected inference to this setting, yielding confidence intervals that explicitly account for bias induced by max–min operators. Third, we develop a scalable uncertainty quantification procedure based on Epistemic Neural Networks (ENN) osband2023epistemic, which provides a practical surrogate for resampling-based inference and enables precision correction in modern over-parameterised models. Together, these components provide the first end-to-end approach to finite-sample PNS-bound estimation with modern machine learning models, and establish a principled baseline for individual-level causal inference beyond average effects.
The remainder of the paper is organized as follows. Section 2 introduces the problem setup and reviews population-level PNS bounds. Section 3 presents the proposed anchored architecture, uncertainty modelling, and precision-corrected inference procedure. Section 4 reports empirical results. Section 5 discusses related work, and Section 6 concludes the paper.
We consider a binary treatment $X \in \{0,1\}$, a binary outcome $Y \in \{0,1\}$, and pre-treatment covariates $Z \in \mathcal{Z}$. We work within a structural causal model in which the potential outcomes $Y_x$ and $Y_{x'}$ are well defined, corresponding to the outcomes that would be observed under the interventions $\mathrm{do}(X=1)$ and $\mathrm{do}(X=0)$, respectively.
For each covariate value $z$, we summarize the relevant conditional quantities as follows. The interventional success probabilities are
and the observational joint probabilities are
so that $p_{10}(z) + p_{11}(z) = P(X = 1 \mid Z = z)$ and $p_{01}(z) + p_{00}(z) = P(X = 0 \mid Z = z)$.
Following Tian and Pearl tian2000probabilities, the probability of necessity and sufficiency (PNS) of $X$ for $Y$ is defined as
This quantity represents the fraction of units for which the treatment is both necessary and sufficient for the outcome.
We define the covariate-conditional PNS as
and view the marginal PNS as $\mathrm{PNS} = \mathbb{E}_Z\!\left[\mathrm{PNS}(Z)\right]$.
We adopt the structural assumptions used in tian2000probabilities, formulated conditionally on $Z = z$.
For a fixed $z$, the conditional PNS is a linear functional of the latent probabilities. Tian and Pearl derive sharp bounds by solving the corresponding linear programs.
In our notation, the bounds are
These satisfy
and the interval is sharp, and the marginal PNS satisfies $\mathrm{PNS} \in \big[\beta_{\mathrm{L}}, \beta_{\mathrm{U}}\big]$, where $\beta_{\mathrm{L}} = \mathbb{E}_{Z}[\mathrm{PNS}_{\mathrm{L}}(z)]$ and $\beta_{\mathrm{U}} = \mathbb{E}_{Z}[\mathrm{PNS}_{\mathrm{U}}(z)]$.
The same framework extends to bounds on the probability of harm $P(Y_{x'} = 1, Y_x = 0 \mid Z = z)$. In the binary setting, $\mathrm{PNS}(z)$ coincides with the probability that the individual treatment effect is positive. Since this quantity is not point-identified, the Tian--Pearl bounds provide the tightest possible interval compatible with the observed data.
While exact in the population limit, Tian–Pearl bounds must in practice be estimated from finite samples. A naive plug-in approach—fitting separate models for experimental and observational distributions—suffers from two fundamental structural flaws. First, independently estimated models for $p_{xy}(z)$ and $\mu_x(z)$ routinely violate the cross-regime feasibility constraints (Eq. (ref)). Such estimates fail the compatibility requirements of a shared structural causal model, yielding internally inconsistent bounds.
Second, the max/min structure of PNS bounds induces systematic extremum bias even with consistent probability estimates. To illustrate (Fig. (ref)), let the true lower-bound terms be $\ell \in \{0.05, 0.40, 0.28, 0.30\}$, yielding $\mathrm{PNS}_L(z) = 0.40$. A noisy estimator yielding $\hat{\ell} \in \{0.02, 0.40, 0.68, 0.22\}$ produces a plug-in bound of $\hat{\mathrm{PNS}}_L(z) = \max \hat{\ell} = 0.68 > 0.40$. Because extremum operators naturally select the most severe estimation errors, this "winner's curse" systematically overestimates lower bounds and underestimates upper bounds, yielding spuriously tight intervals. Valid inference therefore demands a framework that guarantees structural causal compatibility, quantifies finite-sample uncertainty, and rigorously corrects for extremum bias. We achieve this by unifying a constrained multi-head architecture with scalable neural bootstrap surrogates osband2023epistemic and precision-corrected intersection bounds.
The structural constraints (Eq. (ref)) impose algebraic relationships between the observational joint probabilities $p_{xy}(z)$ and the interventional probabilities $\mu_x(z)$. Rather than estimating these six quantities independently, we design an architecture that enforces these relationships by construction. The key idea is to treat the joint probabilities $p_{xy}(z)$ as anchors: primary outputs trained directly on observational data. The interventional probabilities $\mu_x(z)$ are then parameterized as constrained functions of these anchors.
Our architecture consists of a shared representation network $\phi_\theta : \mathcal{Z} \to \mathbb{R}^d$ that maps covariates to a learned embedding, followed by three output heads. The first head outputs logits over the four joint outcomes $(X,Y) \in \{0,1\}^2$, which are passed through a softmax to produce the joint probabilities $p_{00}(z), p_{01}(z), p_{10}(z), p_{11}(z)$. The remaining two heads output unconstrained scalar values $\delta_0(z), \delta_1(z) \in \mathbb{R}$, which are transformed into the interventional probabilities $\mu_0(z), \mu_1(z)$ via a constraint-enforcing reparameterization layer.
We enforce the Tian--Pearl constraints by reparameterizing the interventional probabilities as convex combinations of their feasible bounds:
where $\sigma(\cdot)$ is the sigmoid function. For each fixed $z$, the map $\delta_x(z) \mapsto \mu_x(z)$ is smooth and strictly increasing, since $\mu_x(z)$ is an affine function of $\sigma\bigl(\delta_x(z)\bigr)$ and $\sigma$ is smooth and strictly increasing on $\mathbb{R}$. Writing the lower and upper feasible bounds for $\mu_x(z)$ as $\ell_x(z)$ and $u_x(z)$, the derivative $\partial \mu_x(z) / \partial \delta_x(z) = \bigl(u_x(z) - \ell_x(z)\bigr)\, \sigma\bigl(\delta_x(z)\bigr)\bigl(1 - \sigma\bigl(\delta_x(z)\bigr)\bigr)$ is strictly positive whenever the interval $[\ell_x(z), u_x(z)]$ is non-degenerate. Thus the reparameterization is a smooth, order-preserving bijection from $\mathbb{R}$ onto the open feasible interval for $\mu_x(z)$. Gradient-based optimization over $\delta_x(z)$ therefore explores the entire feasible region for $\mu_x(z)$ while preserving the constraints. Since $\sigma(\delta) \in (0,1)$, each $\mu_x(z)$ lies within its feasible interval for any network weights, and the Tian--Pearl consistency constraints hold automatically.
The joint-probability head is trained on observational data using a cross-entropy loss over the four-class outcome $(X,Y)$. The auxiliary heads are trained on experimental data using binary cross-entropy losses on the derived interventional probabilities $\mu_0(z)$ and $\mu_1(z)$. We keep the gradients from these losses separate: the experimental loss updates only the auxiliary heads, treating the joint probabilities as fixed when computing the feasible bounds. This preserves a clean anchoring of $\mu_x(z)$ to $p_{xy}(z)$ while still allowing the representation $\phi_\theta$ to be shared across tasks.
Our inference procedure builds on the intersection bounds methodology of chernozhukov2013intersection. Let $\psi^\star$ denote a scalar target parameter that is only partially identified. Suppose $\psi^\star$ is known to lie in the intersection
where $\psi_\ell(v)$ and $\psi_u(v)$ are lower and upper bounding functions indexed by $v \in \mathcal{V}$. Equivalently, the identified set for $\psi^\star$ is the interval $[\psi_\ell^\star, \psi_u^\star]$ with
Given estimators $\hat{\psi}_\ell(v)$ and $\hat{\psi}_u(v)$ with pointwise standard errors $s_\ell(v)$ and $s_u(v)$, a naive plug-in approach forms $\sup_{v} \hat{\psi}_\ell(v)$ and $\inf_{v} \hat{\psi}_u(v)$. These plug-in estimators are systematically biased. The supremum tends to pick out values of $v$ where estimation noise is positive ($\hat{\psi}_\ell(v) > \psi_\ell(v)$), producing an upward-biased estimate of the lower bound. Conversely, the infimum tends to pick out values where noise is negative, producing a downward-biased estimate of the upper bound. The resulting bounds are therefore spuriously tight relative to the true identified set.
The precision-corrected estimators of chernozhukov2013intersection adjust for this optimization bias by penalizing regions of high uncertainty. The uncertainty at each $v$ is summarized by $s_\ell(v)$ and $s_u(v)$, which measure the sampling variability of $\hat{\psi}_\ell(v)$ and $\hat{\psi}_u(v)$. The precision-corrected lower and upper bound estimators are defined as
where $\kappa_\alpha$ is a critical value chosen to ensure coverage at level $1-\alpha$. Intuitively, the penalty $\kappa_\alpha\,s_\ell(v)$ (resp.\ $\kappa_\alpha\,s_u(v)$) downweights regions of $\mathcal{V}$ where the estimator is imprecise, preventing the supremum (infimum) from concentrating on points with favourable but unreliable estimation error. The theoretical choice of $\kappa_\alpha$ is based on the standardized error process. Focusing on the upper bound (the lower bound is analogous), define
Under regularity conditions, chernozhukov2013intersection show that $T_n(v)$ can be strongly approximated by a Gaussian process $T_n^\star(v)$ satisfying $\sup_{v \in \mathcal{V}} \bigl|T_n(v) - T_n^\star(v)\bigr| \xrightarrow{p} 0$. This approximation allows $\kappa_\alpha$ to be chosen as the $(1-\alpha)$-quantile of $\sup_{v \in \mathcal{V}} T_n^\star(v)$. In practice, one estimates the covariance structure of $T_n^\star$ from the data, simulates $B$ draws from the approximating Gaussian process, computes the simulated suprema, and takes their empirical $(1-\alpha)$-quantile as $\kappa_\alpha$.
\paragraph{Coverage guarantee.} Let $\psi_u^\star = \inf_{v \in \mathcal{V}} \psi_u(v)$ denote the true upper bound. The precision-corrected estimator $\hat{\psi}_{u,\alpha}$ satisfies
By construction, $\kappa_\alpha$ is the $(1-\alpha)$ - quantile of $\sup_{v \in \mathcal{V}} T_n^\star(v)$, and the Gaussian approximation implies $P(\sup_{v} T_n(v) > \kappa_\alpha) \approx P(\sup_{v} T_n^\star(v) > \kappa_\alpha) = \alpha$. Hence the interval $[\hat{\psi}_{\ell,\alpha}, \hat{\psi}_{u,\alpha}]$ achieves asymptotically valid coverage for the identified set of $\psi^\star$ at level $1-\alpha$.
\paragraph{Application to PNS bounds.} The Tian--Pearl bounds on $\mathrm{PNS}(z)$ take the form of a maximum over four lower bound terms and a minimum over four upper bound terms:
where each $\ell_z(j)$ and $u_z(k)$ is a known function of the observational and interventional probabilities at covariate value $z$. This is a finite dimensional instance of the intersection bounds problem. The precision correction serves two roles in this application. As an estimator, it removes the systematic bias induced by the max and min operators selecting components with favourable noise; this debiasing holds for any estimator with pointwise standard errors. As an inference procedure, it provides simultaneous coverage of the bounding functionals which, combined with the algebraic validity of the Tian--Pearl inequalities at the true distribution, yields coverage of $\mathrm{PNS}(z)$. This second guarantee requires the regularity conditions R1--R3 established in Appendix (ref), where we show that M-estimation regularity on the heads of our anchored network is sufficient (Proposition (ref), Corollary (ref), Theorem (ref)). To our knowledge, this provides the first formal coverage guarantee for PNS bounds using finite-sample estimation.
To apply the precision corrected framework, we require estimates of the standard errors $s_{\ell_z}(j)$ and $s_{u_z}(m)$, as well as the critical values $\kappa_\alpha$ which depend on the joint distribution of the underlying influence-function process. The classical approach estimates standard errors using the asymptotic variance of the estimator. This variance is obtained via an influence-function representation that requires Hessian inversion. Critical values are then computed using a multiplier bootstrap over the resulting influence-function process. For a network with $p$ parameters and $n$ samples, this incurs a cost of $O(np^2 + p^3)$ per covariate value. In high-dimensional causal estimation where $p \gg 10^5$, this is computationally intractable. To make the precision correction practical in this regime, we require an alternate approach for estimating $s_{\ell_z}(j)$ and $s_{u_z}(m)$ that scales favourably with model size while capturing epistemic uncertainty arising from finite data. The next section presents such a method.
To avoid the prohibitive computational cost of Hessian-based influence-function methods, we adopt an alternative framework based on Epistemic Neural Networks (ENN) osband2023epistemic. Our use of ENNs is motivated by a frequentist perspective. Rather than interpreting epistemic variability as a Bayesian posterior, we use the ENN as an amortized mechanism for generating randomized perturbations of a fitted estimator, in the same spirit as classical multiplier bootstrap procedures. This perspective emphasizes the role of ENNs as a scalable approximation to resampling-based inference, rather than as a model of belief uncertainty. In the following sections, we describe how the ENN architecture can be structured to produce coherent perturbations of the empirical estimator, providing a practical substitute for explicit bootstrap reweighting in high-dimensional settings.
The ENN framework augments a base network with an epinet and an explicit epistemic index. This construction is agnostic to the internal architecture of the base network and allows us to sample from our constrained multi-head model without architectural changes. In our setting, ENNs are used to generate stochastic perturbations of a fitted estimator that approximate finite-sample variability osband2018randomized,osband2023epistemic. These properties make ENNs a natural choice for our precision-corrected PNS bounds.
A standard neural network defines a function $f_\theta \colon \mathcal{Z} \to \mathbb{R}^k$. An ENN extends this to a map $f_\theta \colon \mathcal{Z} \times \mathcal{E} \to \mathbb{R}^k$, where $\mathcal{E}$ is the space of epistemic indices and $P_\zeta$ is a reference distribution on $\mathcal{E}$. For fixed $z$, the distribution of $f_\theta(z,\zeta)$ under $\zeta \sim P_\zeta$ represents the model's epistemic uncertainty about the prediction at $z$. When we consider a sequence of inputs $z_1,\ldots,z_\tau$, the ENN defines a joint predictive distribution
Dependence of the predictions across $t$ through a shared $\zeta$ encodes reducible uncertainty that is coherent across inputs. In our setting the ENN outputs the six probability estimates that define the Tian--Pearl bounds. For each $(z,\zeta)$ we obtain \[ f_\theta(z,\zeta) = \bigl( p_{xy}(z;\theta,\zeta)_{x,y\in\{0,1\}}, \; \mu_0(z;\theta,\zeta), \; \mu_1(z;\theta,\zeta) \bigr). \] and the anchored architecture ensures that these probabilities satisfy the algebraic constraints for every $\zeta$. For a fixed $z$ we draw $\zeta_1,\ldots,\zeta_M \sim P_\zeta$. For each $m$ we compute the four lower bound terms $\ell_z^{(m)}(j)$ and four upper bound terms $u_z^{(m)}(k)$ from the corresponding probabilities. We approximate the mean and standard deviation of each bound term under the epistemic distribution by Monte Carlo averages:
The quantities $\bar{\ell}_z(j)$ and $\bar{u}_z(k)$ play the role of $\hat{\psi}_\ell(v)$ and $\hat{\psi}_u(v)$ in the generic intersection bounds notation. The quantities $s_{\ell_z}(j)$ and $s_{u_z}(k)$ are ENN based analogues of the standard errors $s_\ell(v)$ and $s_u(v)$. Note that compared to the $O(np^2 + p^3)$ Hessian operations required by influence function computation, ENNs require only $O(M)$ forward passes, where $M$ is the number of epistemic index samples at test time.
\paragraph{Critical value computation.} The precision corrected framework also requires a joint critical value that respects the dependence structure across the bound terms. In chernozhukov2013intersection this critical value is derived from a Gaussian approximation to the studentized estimation process. In our setting we cannot construct such a process directly, because we do not have an influence function representation. Instead we use the empirical distribution of ENN samples as a surrogate.
For each epistemic draw $m$ we compute standardized deviations {
} We then form the maxima and minima
{
}
The critical value $\kappa_\alpha^L$ is chosen as the empirical $(1-\alpha)$ quantile of the sample $\{\tilde{W}_z^{L,(m)}\}_{m=1}^M$. The critical value $\kappa_\alpha^U$ is defined analogously from $\{\tilde{W}_z^{U,(m)}\}_{m=1}^M$. This construction mirrors the Gaussian process simulation step in chernozhukov2013intersection, but uses the ENN induced empirical distribution instead of an analytic Gaussian coupling.
\paragraph{Precision-corrected PNS bounds.} The precision corrected PNS bounds at covariate value $z$ penalise terms with high epistemic variability. They are given by
The interval $[\hat{L}_\alpha(z), \hat{U}_\alpha(z)]$ provides a precision corrected estimate of the identified set for $\mathrm{PNS}(z)$ that adjusts for selection bias induced by the max and min operators. Complete pseudocode for the proposed procedure is provided in Appendix (ref) - Algorithm (ref)
\paragraph{Why ENNs?} ENNs provide a scalable mechanism to generate the jointly coherent, function-level perturbations required for multiplier-bootstrap-style calibration. In classical frequentist theory, sampling variability is typically approximated via the influence function $\psi(x) = H^{-1} \nabla_{\theta} \ell(x,\hat{\theta})$. Rather than explicitly linearizing the estimator or inverting the Hessian $H$, ENNs directly produce stochastic functional realizations that perturb the empirical risk minimizer. Constrained by a shared architecture, these realizations enforce strict functional coherence across covariates and estimands, yielding a practical approximation of finite-sample variability. This function-space perspective is grounded in Neural Tangent Kernel (NTK) theory, which establishes asymptotically linear, jointly Gaussian limits for over-parameterized networks jacot2018neural, and aligns with coherent representations of predictive uncertainty charpentier2020posterior. While not formally equivalent to the multiplier bootstrap, we posit that ENNs offer a computationally efficient surrogate for the joint variability needed to calibrate the precision-corrected critical value $k_{\alpha}$, a conjecture empirically validated in Section 4.
This section evaluates the finite-sample behavior of machine learning estimators for PNS bounds. The experiments are designed to answer four questions. First, do standard plug-in estimators produce probability estimates that satisfy the algebraic constraints of the Tian–Pearl system, and does constraint violation translate into degraded PNS bound quality? Second, how much of the coverage gain is attributable to precision correction alone, and how severe is the extremum bias without it? Third, does the EpiNets uncertainty estimator produce coverage comparable to the multiplier bootstrap with exact influence functions, and at what computational cost? Fourth, how does bound quality vary with the amount of available experimental data?
We evaluate our approach in two settings. In both, the structural causal model (SCM) is fully specified, but the learner observes only a subset of covariates. This induces hidden confounding while allowing exact computation of ground-truth targets. The first setting uses the synthetic SCM from li2022probabilities, featuring a binary treatment $X$, a binary outcome $Y$, and a 20-dimensional binary covariate vector $Z$. The learner observes 15 dimensions of $Z$; the remaining 5 are hidden confounders. The second setting uses 200 real-valued covariates from the ACIC2019 dataset, augmented with 5 synthetic dimensions. Treatments and potential outcomes are generated from an SCM over all 205 variables. The learner observes only the original 200 covariates, leaving the 5 synthetic dimensions hidden. Because the data-generating SCM is known in both settings, the ground-truth $\text{PNS}(z)$ and sharp Tian–Pearl bounds are exact. This permits direct evaluation of bound quality and coverage. All results are evaluated on disjoint test splits and averaged over $K$ independent replicates. Additional details regarding the data generating processes can be found in Appendix (ref), and (ref).
Metrics. For each test point $z_i$, the evaluated method outputs an interval $[\ell_i, u_i]$. Given the true $\mathrm{PNS}(z_i)$ and oracle identification region $[L(z_i), U(z_i)]$, we report: Constraint violation rate, the fraction of points where estimated probabilities violate the constraints given in Eq. (ref); Point coverage and Identified set coverage, the fractions of points where $\mathrm{PNS}(z_i) \in [\ell_i, u_i]$ and $[L(z_i), U(z_i)] \subseteq [\ell_i, u_i]$ respectively, both targeting a 95% nominal level; Bound width, the average length $u_i - \ell_i$; Interval score winkler1972decision, a proper scoring rule penalising mis-coverage and width. Full definitions are in Appendix (ref).
We compare four estimators structured as an ablation: each adds one component of the proposed method, allowing the contribution of constraint enforcement and precision correction to be isolated.
S-learner and T-learner plug-in. These are standard baselines that estimate the quantities entering the Tian--Pearl bounds independently, without any structural coordination. The S-learner fits a single model for $P(Y=1 \mid X, Z)$ and obtains $\hat{\mu}_1(z)$ and $\hat{\mu}_0(z)$ by evaluating it at $X=1$ and $X=0$ respectively. The T-learner fits separate models on the treated and control arms to estimate $\mu_1(z)$ and $\mu_0(z)$ directly. In both cases, the joint observational probabilities $p_{xy}(z) = P(X=x, Y=y \mid Z=z)$ are estimated by a separate multinomial model trained on the observational data. All estimates are substituted directly into the Tian--Pearl formulas with no constraint enforcement and no precision correction.
Anchored network (no CLR correction). The anchored architecture introduced in Section 3.1 jointly estimates all six atomic probabilities $(\mu_1, \mu_0, p_{11}, p_{10}, p_{01}, p_{00})$ and enforces the Tian--Pearl consistency constraints by construction via the reparameterisation in Eq. (10)--(11). Point estimates are plugged into the bound formulas to obtain interval endpoints. No precision correction is applied. This baseline isolates the contribution of constraint enforcement from that of the CLR correction: any improvement over the plug-in methods is due to the anchored architecture alone, and any remaining gap to the full method is attributable to the absence of precision correction.
Anchored + ENN + CLR (full method) MB Full. The complete proposed method as described in Section 3: the anchored architecture with an ENN for epistemic uncertainty quantification, combined with precision-corrected inference via the CLR construction.
Last-layer influence function bootstrap (MB Last Lyr). For the high-dimensional comparison in Section 4.7, we additionally include a frequentist uncertainty estimator that approximates the influence function by differentiating only the final linear layer of the anchored network, reducing the Hessian solve to a $p_{\mathrm{last}}$-dimensional problem. Standard errors derived from a multiplier bootstrap over these approximate influence functions are used in place of ENN-based epistemic uncertainty, with the same CLR precision correction applied. Full network influence function computation is not run in this setting due to the $O(np^2 + p^3)$ cost of Hessian inversion, which is prohibitive at the network sizes required for the ACIC covariate structure.
Additional details regarding the experimental setup can be found in Appendix (ref).
Table (ref) reports results for the low-dimensional Li et al. benchmark ($n_\mathrm{obs}$ = 100k, $n_\mathrm{exp}$ = 50k) and the high-dimensional ACIC setting ($n_\mathrm{obs}$ = 2M, $n_\mathrm{exp}$ = 1M). Bracketed values indicate 95% confidence intervals over 1000 independent replicates. Replicates for the low dimensional setting use fresh samples from the exact SCM, and the ACIC replicates are drawn by subsampling a fixed super-population of $3$M observational and $3$M experimental records. In the low dimensional setting, S- and T-learner plug-ins achieve point coverage of $0.745$ and $0.724$ (interval scores $0.730$, $0.987$); their narrow widths reflect infeasible probability configurations rather than sharpness. The anchored network eliminates constraint violations and improves coverage to 0.823, but its interval score (0.810) confirms constraint enforcement alone is insufficient. Causal EpiNet achieves the best performance: valid point coverage at $0.982$ $[0.980, 0.986]$, nominal identified set coverage at $0.956 \ [0.953, 0.958]$, and a score of 0.347. The full-network bootstrap achieves comparable point coverage (0.985) but produces conservative, overly wide intervals (width 0.370 vs. 0.330), resulting in a substantially worse score $(0.639)$. The last-layer bootstrap undercovers at $0.881 \ [0.877, 0.883]$, consistent with uncaptured representation uncertainty. In the high-dimensional ACIC setting, computing the full-network influence function is computationally intractable; it is therefore omitted. Causal EpiNet uniquely maintains valid point coverage $(0.973 \ [0.968, 0.977])$ with a score of $0.363$. Its identified set coverage drops to $0.860$, isolating residual high-dimensional bias in estimating $p_{xy}(z)$ and $\mu_x(z)$ rather than precision miscalibration. The uncorrected estimator undercovers at $0.858 [0.857, 0.891]$ score $(0.782)$. The last-layer bootstrap undercovers at $0.875$ $[0.855, 0.895]$ (score $0.720$); computing influence functions only through the final linear layer produces an incomplete estimate of the asymptotic variance of $\hat{p}_{xy}(z)$ and $\hat{\mu}_x(z)$, which systematically underestimates $s_\ell(j)$ and $s_u(k)$.
We evaluate Causal EpiNet across $(n_{\mathrm{obs}}, n_{\mathrm{exp}})$ regimes to assess behavior under the practically common condition of scarce experimental data. Figures (ref) and (ref) show that as $n_{\mathrm{exp}}$ increases from 20k to 100k with $n_{\mathrm{obs}} = 100\text{k}$ fixed, identified set coverage remains nominal throughout while bound width decreases monotonically. This confirms that ENN uncertainty estimates correctly track experimental data availability without sacrificing calibration. Full results, including balanced scaling up to 1M/1M, are reported in Appendix (ref). To assess computational efficiency, we record the training and inference wall-clock times for all evaluated methods. Because they require only a single forward pass rather than repeated model retraining or costly influence function evaluations, ENNs incur substantially lower computational overhead than the bootstrap baselines. A comprehensive runtime breakdown is deferred to the Appendix (ref).
Note that the point coverage intrinsically exceeds identified (ID) set coverage. Because $\mathrm{PNS}(z) \in [L(z), U(z)]$ deterministically, $[L(z_i), U(z_i)] \subseteq [\ell_i, u_i]$ implies $\mathrm{PNS}(z_i) \in [\ell_i, u_i]$, but not conversely; an interval can miss an oracle endpoint yet still cover the true parameter. The CLR correction specifically targets point coverage. Thus, nominal ID set coverage in the low-dimensional setting reflects precise first-stage estimation rather than inferential conservatism. In the ACIC benchmark, residual bias in $\hat{p}_{xy}(z)$ and $\hat{\mu}_x(z)$ degrades ID set coverage to $0.860$, but point coverage remains valid at $0.973$ because the true $\mathrm{PNS}(z)$ lies strictly interior to the estimated bounds. Sub-nominal ID set coverage therefore diagnoses first-stage estimation error, not inferential failure.
We have presented Causal EpiNets, a framework for finite-sample inference on individual-level causal effects via PNS bounds. By combining a constrained anchored architecture with ENN-based uncertainty quantification and CLR precision-corrected inference, the method achieves near-nominal coverage without requiring unconfoundedness or positivity on the observational regime — assumptions that are routinely violated in practice and that methods based on CATE estimation cannot dispense with. The only requirements are access to experimental and observational data from the same population and the stability assumption of Section 2. The framework extends directly to the probability of harm $P(Y_{x'} = 1, Y_x = 0 \mid Z = z)$ and, together with PNS, provides a complete characterisation of individual-level treatment effect heterogeneity in the binary setting. We view this as a strong baseline for a line of work that has lacked practical estimation methods despite a rich 20-year-old identification theory.
Limitations. The current framework is restricted to binary treatments and outcomes. While theoretical progress exists for multi-valued and continuous treatments kawakami24a, li2024probabilities_nonbinary, developing a practical estimation procedure remains an open problem. Furthermore, scalability to high-dimensional unstructured covariates, such as text or images, requires further investigation. In these settings, representation learning becomes the dominant challenge. Although our anchored architecture enforces consistency constraints over a learned embedding, the generalization of these embeddings and the behavior of ENN epistemic uncertainty over them require future study.
Future work. A primary obstacle to progress in this domain is the absence of large-scale benchmark datasets combining experimental and observational data. This data scarcity restricts empirical evaluation to synthetic and semi-synthetic environments. We encourage the development of public benchmarks to facilitate the rigorous evaluation of individual-level causal inference methods. Methodologically, natural next steps include extending identification and estimation theory beyond the binary setting and incorporating modern representation learning for multi-modal and document-level covariates into the anchored ENN framework.