The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
45,506 characters
Causal EpiNets: Precision-corrected Bounds on Individual Treatment Effects using Epistemic Neural Networks
\maketitle
\begin{abstract}
Individual treatment effects are not point-identified from data.
The Probability of Necessity and Sufficiency (PNS) circumvents
this limitation by characterizing individual-level causality through
intersection bounds derived from combined experimental and observational
data. In finite samples, however, standard plug-in estimators systematically
fail: they violate structural probability constraints and suffer from extremum
bias induced by max–min operators, yielding spuriously narrow intervals.
We propose a neural framework for finite-sample PNS estimation that resolves
both pathologies. We introduce an anchored neural architecture that guarantees
structural constraint satisfaction by construction. To correct extremum bias,
we employ precision-corrected intersection-bound inference, leveraging Epistemic
Neural Networks for scalable, high-dimensional uncertainty quantification.
Empirical evaluations confirm that this approach maintains nominal coverage
and exact constraint validity in high-dimensional regimes where standard
estimators systematically undercover.
\end{abstract}
\section{Introduction}
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:
\emph{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 \emph{Probability of Necessity and Sufficiency} (PNS),
introduced by \citet{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 \citet{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
\citep{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) \citep{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.
\section{Problem Setup}
\label{sec:pns-setup}
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 \emph{interventional success probabilities}
are
\begin{equation}
\begin{aligned}
\mu_1(z) &:= P\!\left(Y = 1 \mid \mathrm{do}(X = 1), Z = z\right), \\
\mu_0(z) &:= P\!\left(Y = 1 \mid \mathrm{do}(X = 0), Z = z\right),
\end{aligned}
\end{equation}
and the \emph{observational joint probabilities} are
\begin{equation}
p_{xy}(z) := P(X = x, Y = y \mid Z = z), \ \ x,y \in \{0,1\},
\end{equation}
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)$.
\subsection{Probability of necessity and sufficiency}
Following Tian and Pearl~\cite{tian2000probabilities}, the
\emph{probability of necessity and sufficiency} (PNS) of $X$ for $Y$
is defined as
\begin{equation}
\mathrm{PNS}
:= P\!\left(Y_x = 1,\; Y_{x'} = 0\right).
\label{eq:pns-def}
\end{equation}
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
\begin{equation}
\mathrm{PNS}(z)
:= P\!\left(Y_x = 1,\; Y_{x'} = 0 \mid Z = z\right),
\label{eq:pns-z-def}
\end{equation}
and view the marginal PNS as
$\mathrm{PNS} = \mathbb{E}_Z\!\left[\mathrm{PNS}(Z)\right]$.
\subsubsection{Structural assumptions and constraint system}
\label{subsec:pns-assumptions}
We adopt the structural assumptions used in
\cite{tian2000probabilities}, formulated conditionally on $Z = z$.
\begin{assumption}[{\bf Experimental identification}]
In the experimental regime, treatment assignment is randomized, possibly
conditional on $Z$, so that $X$ is independent of the potential outcomes
given $Z$.
Under this protocol, the interventional probabilities are identified by
\begin{align}
\nonumber \mu_x(z)
&:= P(Y = 1 \mid \mathrm{do}(X = x), Z = z)\\
&= P^{\mathrm{exp}}(Y = 1 \mid X = x, Z = z),
\qquad x \in \{0,1\}.
\end{align}
\end{assumption}
\begin{assumption}[{\bf Stability across regimes}]
Experimental and observational data share the same structural equation
for $Y$ and the same distribution of exogenous variables.
The treatment assignment mechanism for $X$ may differ across regimes:
in the experimental regime, $X$ is assigned at random, while in the
observational regime $X$ may depend on $Z$ and unobserved factors.
Under this shared structural causal model, the observational joint
distribution $(X,Y)\mid Z$ restricts the set of admissible interventional
probabilities $(\mu_0(z), \mu_1(z))$ through the compatibility constraints
\begin{equation}
\begin{aligned}
p_{11}(z) &\le \mu_1(z) \le 1 - p_{10}(z), \\
p_{01}(z) &\le \mu_0(z) \le 1 - p_{00}(z),
\end{aligned}
\label{eq:observational-constraints}
\end{equation}
which must hold for any values of $(\mu_0(z), \mu_1(z))$ consistent with
the observational distribution.
\end{assumption}
\subsubsection{Conditional and marginal PNS bounds}
\label{subsec:pns-bounds}
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
\begin{equation}
\begin{aligned}
\mathrm{PNS}_{\mathrm{L}}(z)
:= \max\Big\{&
0,\;
\mu_1(z) - \mu_0(z), \\
& p_{11}(z) + p_{01}(z) - \mu_0(z), \\
& \mu_1(z) - p_{11}(z) - p_{01}(z)
\Big\},
\end{aligned}
\label{eq:pns-lower-1}
\end{equation}
\begin{equation}
\begin{aligned}
\mathrm{PNS}_{\mathrm{U}}(z)
:= \min\Big\{&
\mu_1(z),\;
1 - \mu_0(z), \\
& p_{11}(z) + p_{00}(z), \\
& \mu_1(z) - \mu_0(z)
+ p_{10}(z) + p_{01}(z)
\Big\}.
\end{aligned}
\label{eq:pns-upper-1}
\end{equation}
These satisfy
\begin{equation}
\mathrm{PNS}_{\mathrm{L}}(z)
\le \mathrm{PNS}(z)
\le \mathrm{PNS}_{\mathrm{U}}(z),
\end{equation}
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.
\section{Proposed Approach/Contribution}
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{eq:observational-constraints}). Such estimates fail the compatibility
requirements of a shared structural causal model, yielding internally inconsistent bounds.
\begin{figure}[ht]
\centering
\includegraphics[width=0.8\columnwidth]{images/selection_bias_diagram.png}
\caption{Overestimation bias in the lower bound.}
\label{fig:bias_illustration}
\vspace{-1em}
\end{figure}
Second, the max/min structure of PNS bounds induces systematic
extremum bias even with consistent probability estimates.
To illustrate (Fig.~\ref{fig:bias_illustration}), 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 \citep{osband2023epistemic}
and precision-corrected intersection bounds.
\subsection{Ensuring estimates satisfy probability axioms}
The structural constraints (Eq.~\ref{eq:observational-constraints})
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 \emph{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:
\begin{align}
\mu_0(z)
&= p_{01}(z)
+ \bigl(1 - p_{00}(z) - p_{01}(z)\bigr)\,\sigma\bigl(\delta_0(z)\bigr),
\\
\mu_1(z)
&= p_{11}(z)
+ \bigl(1 - p_{10}(z) - p_{11}(z)\bigr)\,\sigma\bigl(\delta_1(z)\bigr),
\end{align}
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.
\subsection{Bias correction and uncertainty quantification}
Our inference procedure builds on the intersection bounds methodology of
\cite{chernozhukov2013intersection}.
Let $\psi^\star$ denote a scalar target parameter that is only partially
identified.
Suppose $\psi^\star$ is known to lie in the intersection
\begin{align}
\psi^\star
\in \bigcap_{v \in \mathcal{V}}
\bigl[\,\psi_\ell(v),\,\psi_u(v)\,\bigr],
\end{align}
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
\begin{align}
\psi_\ell^\star
= \sup_{v \in \mathcal{V}} \psi_\ell(v),
\qquad
\psi_u^\star
= \inf_{v \in \mathcal{V}} \psi_u(v).
\end{align}
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 \cite{chernozhukov2013intersection}
adjust for this \textbf{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
\begin{minipage}{0.48\linewidth}
\begin{equation}
\hat{\psi}_{\ell,\alpha}
= \sup_{v \in \mathcal{V}}
\bigl[\hat{\psi}_\ell(v) - \kappa_\alpha\,s_\ell(v)\bigr]
\end{equation}
\end{minipage}\hfill
\begin{minipage}{0.48\linewidth}
\begin{equation}
\hat{\psi}_{u,\alpha}
= \inf_{v \in \mathcal{V}}
\bigl[\hat{\psi}_u(v) + \kappa_\alpha\,s_u(v)\bigr]
\end{equation}
\end{minipage}
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
\begin{equation}
T_n(v)
= \frac{\hat{\psi}_u(v) - \psi_u(v)}{s_u(v)}.
\end{equation}
Under regularity conditions, \citeauthor{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
\begin{equation}
P\bigl(\psi_u^\star \leq \hat{\psi}_{u,\alpha}\bigr)
\geq 1 - \alpha + o(1).
\end{equation}
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:
\begin{minipage}{0.48\linewidth}
\begin{equation}
L(z)
= \max_{j \in \{1,\ldots,4\}} \ell_z(j)
\end{equation}
\end{minipage}\hfill
\begin{minipage}{0.48\linewidth}
\begin{equation}
U(z)
= \min_{k \in \{1,\ldots,4\}} u_z(k),
\end{equation}
\end{minipage}
\noindent
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{sec:pns-theory}, where we show that
M-estimation regularity on the heads of our anchored network is
sufficient (Proposition~\ref{prop:uniform-mest},
Corollary~\ref{cor:pns-parametric},
Theorem~\ref{thm:uniform-pns}).
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.
\vspace{-1em}
\subsubsection{Epistemic Neural Networks for Precision-Corrected Inference}
To avoid the prohibitive computational cost of Hessian-based influence-function
methods, we adopt an alternative framework based on Epistemic Neural Networks
(ENN) \citep{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 \citep{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
\begin{equation}
\hat{P}_{1:\tau}(y_{1:\tau})
= \int_{\mathcal{E}}
\prod_{t=1}^\tau P\bigl(y_t \mid f_\theta(z_t,\zeta)\bigr)
\, dP_\zeta(\zeta).
\end{equation}
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:
\begin{equation}
\resizebox{\columnwidth}{!}{
$\begin{aligned}
\bar{\ell}_z(j)
&= \frac{1}{M}
\sum_{m=1}^M \ell_z^{(m)}(j),
&
s_{\ell_z}(j)
&= \sqrt{
\frac{1}{M-1}
\sum_{m=1}^M
\bigl(\ell_z^{(m)}(j) - \bar{\ell}_z(j)\bigr)^2
},
\\
\bar{u}_z(k)
&= \frac{1}{M}
\sum_{m=1}^M u_z^{(m)}(k),
&
s_{u_z}(k)
&= \sqrt{
\frac{1}{M-1}
\sum_{m=1}^M
\bigl(u_z^{(m)}(k) - \bar{u}_z(k)\bigr)^2
}.
\end{aligned}$
}
\end{equation}
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 \cite{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
{\small
\begin{equation}
\tilde{W}_z^{(m)}(j) = \frac{\ell_z^{(m)}(j) - \bar{\ell}_z(j)}{s_{\ell_z}(j)},
\quad
\tilde{V}_z^{(m)}(m) = \frac{u_z^{(m)}(m) - \bar{u}_z(m)}{s_{u_z}(m)}.
\end{equation}
}
We then form the maxima and minima
{
\small
\begin{equation}
\tilde{W}_z^{L,(m)} = \max_{j \in \{1,\ldots,4\}} \tilde{W}_z^{(m)}(j),
\quad
\tilde{W}_z^{U,(m)} = \max_{m \in \{1,\ldots,4\}} \tilde{V}_z^{(m)}(m).
\end{equation}
}
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
\cite{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
\begin{align}
\hat{L}_{\alpha, z}
&= \max_{j \in \{1,\ldots,4\}}
\bigl[
\bar{\ell}_z(j)
- \kappa_\alpha^L\,s_{\ell_z}(j)
\bigr]\\
\hat{U}_{\alpha, z}
&= \min_{m \in \{1,\ldots,4\}}
\bigl[
\bar{u}_z(m)
+ \kappa_\alpha^U\,s_{u_z}(m)
\bigr].
\end{align}
\noindent
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{app:exp-details} - Algorithm \ref{alg:pns_bounds}
\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 \citep{jacot2018neural},
and aligns with coherent representations of predictive uncertainty
\citep{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.
\vspace{-0.5em}
\section{Experimental Evaluation}\label{sec:experiments}
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?
\subsection{Data-Generating Processes and Protocol}
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 \citep{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 \citet{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{app:dgp-low-dim}, and \ref{app:dgp-high-dim}.
\textbf{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: \textit{Constraint violation rate}, the fraction of points where estimated probabilities
violate the constraints given in Eq.~\ref{eq:observational-constraints}; \textit{Point coverage}
and \textit{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; \textit{Bound width}, the average length $u_i - \ell_i$;
\textit{Interval score} ~\cite{winkler1972decision}, a proper scoring rule penalising mis-coverage and width. Full definitions are in Appendix~\ref{app:metrics}.
\vspace{-1em}
\subsection{Methods}
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.
\begin{table*}[htbp]
\caption{Results on low-dimensional (top) and high-dimensional (bottom) datasets.}
\label{tab:results}
\vskip 0.1in
\begin{center}
\begin{small}
\begin{sc}
\begin{tabular}{lccccc}
\toprule
Method & \% Valid (\faThumbsUp\,\faArrowUp) & Pt Cov.\ (\faThumbsUp\,\faArrowUp) & ID cov (\faThumbsUp\,\faArrowUp) & Width (\faThumbsDown\,\faArrowDown) & Int.\ Score (\faThumbsDown\,\faArrowDown) \\
\midrule
\multicolumn{6}{c}{\textit{Low-dimensional synthetic}} \\
\midrule
S-learner & 76.29 {\tiny [75.81, 76.78] } & 74.52 {\tiny [73.90, 75.13] } & 22.49 {\tiny [22.21, 22.77] } & 0.1730 {\tiny [0.1716, 0.1745] } & 0.7295 {\tiny [0.7132, 0.7458] } \\
T-learner & 74.63 {\tiny [74.16, 75.11] } & 72.39 {\tiny [71.82, 72.95] } & 27.46 {\tiny [27.22, 27.70] } & 0.1782 {\tiny [0.1772, 0.1793] } & 0.9874 {\tiny [0.9695, 1.0054] } \\
Anchored & 100 & 82.34 {\tiny [82.14, 82.54] } & 36.36 {\tiny [36.11, 36.61] } & 0.2041 {\tiny [0.2001, 0.2081] } & 0.8102 {\tiny [0.8062, 0.8142] } \\
MB Full & 100 & 98.50 {\tiny [98.10, 98.90] } & 97.40 {\tiny [97.0, 97.8] } & 0.3695 {\tiny [0.3665, 0.3725] } & 0.6385 {\tiny [0.6365, 0.6405] } \\
MB Last Lyr. & 100 & 88.10 {\tiny [87.70, 88.30] } & 73.41 {\tiny [73.38, 73.44] } & 0.2045 {\tiny [0.2015, 0.2075] } & 0.5112 {\tiny [0.5072, 0.5152] } \\
Anch.\ EpiNet. & 100 & \textbf{98.19} {\tiny [97.99, 98.59] } & \textbf{95.59} {\tiny [95.30, 95.80] } & \textbf{0.3303} {\tiny [0.3263, 0.3323] } & \textbf{0.3471} {\tiny [0.3451, 0.3501] } \\
\midrule
\multicolumn{6}{c}{\textit{High-dimensional ACIC}} \\
\midrule
Anchored & 100 & 85.80 {\tiny [85.7, 69.1] } & 21.50 {\tiny [20.4, 22.6] } & 0.1630 {\tiny [0.160, 0.165] } & 0.7815 {\tiny [0.7820, 0.7810] } \\
MB Last Lyr. & 100 & 87.50 {\tiny [85.5, 89.5] } & 35.40 {\tiny [33.4, 37.4] } & 0.2010 {\tiny [0.180, 0.220] } & 0.7200 {\tiny [0.7000, 0.7400] } \\
Anch.\ EpiNet. & 100 & \textbf{97.30} {\tiny [96.80, 97.70] } & \textbf{86.00} {\tiny [85.3, 86.7] } & \textbf{0.3060} {\tiny [0.301, 0.310] } & \textbf{0.3630} {\tiny [0.3520, 0.3740] } \\
\bottomrule
\end{tabular}
\end{sc}
\end{small}
\end{center}
\vskip -0.1in
\end{table*}
\textbf{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.
\textbf{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.
\textbf{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.
\textbf{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{app:exp-details}.
\vspace{-1em}
\subsection{Results}
\begin{figure}[htbp]
\centering
\includegraphics[width=\columnwidth]{images/interval_vary_nexp.pdf}
\caption{Bound width and interval score vs.\ experimental sample size
($n_\mathrm{obs} = 100\mathrm{k}$ fixed). Both decrease monotonically as
$n_\mathrm{exp}$ grows, with the gap between them remaining stable.}
\label{fig:sens_int_nexp}
\end{figure}
\begin{figure}[htbp]
\centering
\includegraphics[width=\columnwidth]{images/coverage_vary_nexp.png}
\caption{Coverage vs.\ joint sample size ($n_\mathrm{obs} = n_\mathrm{exp}$).
ID set coverage remains near nominal across all balanced regimes.}
\label{fig:sens_cov_balanced}
\end{figure}
\vspace{-1em}
Table \ref{tab:results} 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{fig:sens_cov_balanced} and~\ref{fig:sens_int_nexp} 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{app:sensitivity}. 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{app:compute}.
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.
\vspace{-0.5em}
\section{Discussion and Conclusion}
\vspace{-0.5em}
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.
\textbf{Limitations.} The current framework is restricted to binary treatments
and outcomes. While theoretical progress exists for multi-valued and continuous
treatments~\citep{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.
\textbf{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.
\newpage
\bibliography{uai2026-template}
\newpage