A Relative Error-Based Evaluation Framework of Heterogeneous Treatment Effect Estimators

## Abstract

Abstract While significant progress has been made in heterogeneous treatment effect (HTE) estimation, the evaluation of HTE estimators remains underdeveloped. In this article, we propose a robust evaluation framework based on relative error, which quantifies performance differences between two HTE estimators. We first derive the key theoretical conditions on the nuisance parameters that are necessary to achieve a robust estimator of relative error. Building on these conditions, we introduce novel loss functions and design a neural network architecture to estimate nuisance parameters and obtain robust estimation of relative error, thereby achieving reliable evaluation of HTE estimators. We provide the large sample properties of the proposed relative error estimator. Furthermore, beyond evaluation, we propose a new learning algorithm for HTE that leverages both the previously HTE estimators and the nuisance parameters learned through our neural network architecture. Extensive experiments demonstrate that our evaluation framework supports reliable comparisons across HTE estimators, and the proposed learning algorithm for HTE exhibits desirable performance.

## 1 Introduction

The estimation of heterogeneous treatment effects (HTEs) has attracted substantial attention across a range of disciplines, including economics (Imbens & Rubin, 2015), marketing (Wager & Athey, 2018b), biology (Rosenbaum, 2020), and medicine (Hernán & Robins, 2020), due to its critical role in understanding individual-level treatment heterogeneity and supporting personalized, context-specific decision-making.
Various methods have been developed to estimate HTEs; see Kunzel et al. (2019); Caron et al. (2022) for comprehensive reviews. Despite their growing popularity, the evaluation and comparison of HTE estimators remain relatively underexplored (Gao, 2025). Assessing estimator performance is crucial in real-world applications, as a reliable evaluation framework can identify the most suitable methods (Curth & Van Der Schaar, 2023), directly impacting downstream tasks.

Evaluating HTEs is inherently challenging, as the ground truth is not available: only one potential outcome is observed for each individual, while HTEs are defined as the difference between two. To address this, researchers often rely on stringent model assumptions (Saito & YasuiAuthors, 2023; Mahajan et al., 2024) or preprocessing techniques (e.g., matching) (Rolling & Yang, 2014) to approximate the unobserved counterfactuals, and obtain an estimated treatment effect.
Our work is motivated by Gao (2025), who introduced relative error to quantify the performance difference between two estimators, thereby reducing the bias caused by using inaccurately estimated treatment effects as ground truth.

Despite the significant contributions of Gao (2025), a notable limitation remains unaddressed. Their estimator
requires that all nuisance parameter estimators (propensity score and outcome regression models) are consistent at a rate faster than $n^{-1/4}$ to achieve consistency and valid confidence intervals for the relative error, which may be too stringent for real-world applications.
In practice, the outcome regression models for potential outcomes heavily rely on model extrapolation. These models are
trained separately within the treated and control groups, yet their predictions are applied across the entire dataset.
When there exists a significant distributional difference between the treated and control groups (Jeong & Namkoong, 2020; Jing Qin & Huang, 2024), the extrapolated predictions from these models are prone to inaccuracy and bias, potentially leading to unreliable conclusions.
Therefore, it is desirable to develop methods that reduce reliance on such extrapolation to ensure more robust and trustworthy evaluations.

To address this limitation, we propose a reliable evaluation approach for HTE estimation that retains the desirable properties of the method in Gao (2025), while relaxing the requirement for consistent outcome regression models. We show that the proposed estimator of relative error is $\sqrt{n}$-consistent, asymptotically normal, and yields valid confidence intervals, provided that the propensity score model is consistent at a rate faster than $n^{-1/4}$, even if the outcome regression model is inconsistent.

This robustness is achieved by carefully exploring the relationships between nuisance parameter models. We first derive the key conditions necessary for robustness and then design a novel loss function for estimating outcome regression models. Moreover, since the proposed method still requires a consistent propensity score model, we introduce novel balance regularizers to mitigate this reliance by encouraging the learned propensity scores to satisfy the balance property (Imai & Ratkovic, 2014), i.e., ensuring that the expectation of measurable functions of covariates, weighted by the inverse propensity scores, are equal between treated and control groups.
Furthermore, by combining the novel loss function with balance regularizers, we design a new neural network architecture that more accurately estimates outcome regression and propensity score models, enabling more reliable relative error estimation and, in turn, more robust HTE evaluations.
The main contributions are summarized as follows.

- •
We reveal the limitations of existing methods and, through theoretical analysis, derive key conditions for estimating the relative error that mitigate these limitations.
- •
We propose a reliable HTE evaluation method by designing novel loss functions and introducing a new neural network, enabling more robust estimation of relative error.
- •
We conduct extensive experiments to demonstrate the effectiveness of the proposed method.

## 2 Preliminaries

### 2.1 Problem Setting

We introduce notations to formulate the problem of interest. For each individual $i$, let $A_{i}\in\mathcal{A}=\{0,1\}$ denote the binary treatment variable, where $A_{i}=1$ and $A_{i}=0$ denote treatment and control. Let $X_{i}\in\mathcal{X}\subset\mathbb{R}^{d}$ be the pre-treatment covariates, and $Y_{i}\in\mathbb{R}$ be the outcome. We adopt the potential outcome framework in causal inference (Rubin, 1974; Neyman, 1990), defining $Y_{i}(0)$ and $Y_{i}(1)$ as the potential outcomes under $A_{i}=0$ and $A_{i}=1$, respectively. Since each individual receives either the treatment or the control, the observed outcome $Y_{i}$ satisfies $Y_{i}=A_{i}Y_{i}(1)+(1-A_{i})Y_{i}(0)$.

The individual treatment effect (ITE) is defined as $Y_{i}(1)-Y_{i}(0)$, which represents the treatment effect for a specific individual $i$. However, since only one of $(Y_{i}(0),Y_{i}(1))$ is observable, ITE is not identifiable without imposing strong assumptions (Hernán & Robins, 2020; Pearl, 2009).
In practice, the conditional average treatment effect (CATE) is often used to characterize “individual" treatment effects, defined by

$$ $\tau(x)=\mathbb{E}[Y_{i}(1)-Y_{i}(0)|X_{i}=x],$ $$

which captures how treatment effects vary across individuals with different covariate values.

###### Assumption 1 (Strongly Ignorability, (Rosenbaum & Rubin, 1983 ) ) .

(i) $(Y_{i}(0),Y_{i}(1))\mid X_{i}$; (ii) $0<e(x)\triangleq\mathbb{P}(A_{i}=1\mid X_{i}=x)<1$ for all $x\in\mathcal{X}$, where
$e(x)$ is the propensity score.

Under the standard strong ignorability assumption, CATE is identified as $\mu_{1}(x)-\mu_{0}(x)$ where $\mu_{a}(x)=\mathbb{E}[Y_{i}\mid X_{i}=x,A_{i}=a]$ for $a=0,1$ are the outcome regression functions, and various methods have been developed for estimating CATE (Wager & Athey, 2018a; Shalit et al., 2017a).
Suppose we have a set of candidate CATE estimators trained on a training set, denoted by $\{\hat{\tau}_{1}(x),\cdots,\hat{\tau}_{K}(x)\}.$ We aim to select the estimator with the highest accuracy in a test dataset $\{(X_{i},A_{i},Y_{i}),i=1,\dots,n\}$, which is of size $n$ and sampled from the super-population $\mathbb{P}$, and is independent of the training set.

###### Assumption 1 (Strongly Ignorability, (Rosenbaum & Rubin, 1983 ) ) .

### 2.2 Evaluation Metrics: Absolute Error and Relative Error

For a given estimator $\hat{\tau}(x)$, its accuracy is typically evaluated using the MSE defined by

$$ $\phi(\hat{\tau})\triangleq\mathbb{E}[(\hat{\tau}(X)-\tau(X))^{2}].$ $$

For any two estimators $\hat{\tau}_{1}(x)$ and $\hat{\tau}_{2}(x)$, the difference in their MSE is

$$ $\displaystyle\delta(\hat{\tau}_{1},\hat{\tau}_{2})\triangleq{}$ $\displaystyle\phi(\hat{\tau}_{1})-\phi(\hat{\tau}_{2})=\mathbb{E}[\hat{\tau}_{1}^{2}(X)-\hat{\tau}_{2}^{2}(X)-2(\hat{\tau}_{1}(X)-\hat{\tau}_{2}(X))\tau(X)].$ $$

Gao (2025) refers to $\phi(\hat{\tau})$ and $\delta(\hat{\tau}_{1},\hat{\tau}_{2})$ as the absolute error and relative error, respectively. In practice, absolute error is used much more frequently than relative error. However, Gao (2025) demonstrated that using relative error as the evaluation metric is superior to using absolute error, both theoretically and experimentally, see Section [3](https://arxiv.org/html/2510.16419v1#S3) for more details. Intuitively, one can see that the key advantage of using relative error over absolute error lies in that it only relies on the first-order term of the unobserved $\tau$, which reduces the impact of estimation error in $\tau$.

Several studies (Gutierrez & Gérardy, 2017; Powers et al., 2017) have used $\mathbb{E}[(Y(1)-Y(0)-\hat{\tau}(X))^{2}]$
to evaluate the estimator $\hat{\tau}(x)$. However, its estimator requires knowing the values of $(Y(0),Y(1))$ that are not observable in real-world applications.
We note that

$$ $\mathbb{E}[(Y(1)-Y(0)-\hat{\tau}(X))^{2}=\mathbb{E}[(\hat{\tau}(X)-\tau(X))^{2}+\mathbb{E}[\text{Var}(Y(1)-Y(0)|X)],$ $$

where the second term on the right-hand side is independent of $\hat{\tau}(x)$. Thus, this metric is essentially equivalent to the absolute error $\phi(\hat{\tau})$, and we will not discuss it further. For clarity, we provide a notation summary tab, but due to limited space, we present it in Appendix [A](https://arxiv.org/html/2510.16419v1#A1).

## 3 Motivation

In this section, we briefly discuss the advantages of using relative error over absolute error,
and then analyze the limitations of the method in Gao (2025), which motivate this work.

The key theoretical advantage of relative error over absolute error is demonstrated through its semiparametric efficient estimators.
A semiparametric efficient estimator is considered optimal (or gold standard) in the sense that it has the smallest asymptotic variance under regularity conditions  (Newey, 1990; van der Vaart, 1998) given the observed test data. Let $\{\tilde{e}(x),\tilde{\mu}_{1}(x),\tilde{\mu}_{0}(x)\}$
be the estimators of $\{e(x),\mu_{1}(x),\mu_{0}(x)\}$, which are the nuisance functions to construct semiparametric efficient estimators of absolute error and relative error. Denote $\tilde{\tau}(x)=\tilde{\mu}_{1}(x)-\tilde{\mu}_{0}(x)$.

Absolute Error. Given $\hat{\tau}(x)$, an estimator of $\phi(\hat{\tau})$ is constructed as

$$ $\displaystyle\hat{\phi}(\hat{\tau})={}$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\{\tilde{\tau}(X_{i})-\hat{\tau}(X_{i})\}^{2}+2(\tilde{\tau}(X_{i})-\hat{\tau}(X_{i}))\left(\frac{A_{i}(Y_{i}-\tilde{\mu}_{1}(X_{i}))}{\tilde{e}(X_{i})}-\frac{(1-A_{i})(Y_{i}-\tilde{\mu}_{0}(X_{i}))}{1-\tilde{e}(X_{i})}\right).$ $$

Under Assumption [1](https://arxiv.org/html/2510.16419v1#Thmassumption1), $\hat{\phi}(\hat{\tau})$ is $\sqrt{n}$-consistent, asymptotically normal, and semiparametric efficient, provided that the estimated nuisance parameter satisfy the key Condition [1](https://arxiv.org/html/2510.16419v1#Thmcondition1).

###### Condition 1 .

$\mathbb{E}[(\tilde{e}(X)-e(X))^{2}]=o_{\mathbb{P}}(n^{-1/2})$, $\mathbb{E}[(\tilde{\mu}_{a}(X)-\mu_{a}(X))^{2}]=o_{\mathbb{P}}(n^{-1/2})$ for $a=0,1$.

Relative Error. Likewise, we can construct the estimator of $\delta(\hat{\tau}_{1},\hat{\tau}_{2})$ given as

$$ $\displaystyle\hat{\delta}(\hat{\tau}_{1},\hat{\tau}_{2})={}\frac{1}{n}\sum_{i=1}^{n}\varphi(Z_{i};\tilde{\mu}_{0},\tilde{\mu}_{1},\tilde{e}),\;\text{where}\;Z_{i}\triangleq(A_{i},X_{i},Y_{i}),$ $$

$$ $\displaystyle\varphi(Z_{i};\tilde{\mu}_{0},\tilde{\mu}_{1},\tilde{e})\triangleq{}$ $\displaystyle\{\hat{\tau}_{1}^{2}(X_{i})-\hat{\tau}_{2}^{2}(X_{i})\}-2(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\cdot$ $\displaystyle\left(\frac{A_{i}(Y_{i}-\tilde{\mu}_{1}(X_{i}))}{\tilde{e}(X_{i})}+\tilde{\mu}_{1}(X_{i})-\frac{(1-A_{i})(Y_{i}-\tilde{\mu}_{0}(X_{i}))}{1-\tilde{e}(X_{i})}-\tilde{\mu}_{0}(X_{i})\right).$ $$

Under Assumption [1](https://arxiv.org/html/2510.16419v1#Thmassumption1), the $\sqrt{n}$-consistency, asymptotic normality, and semiparametric efficiency of $\hat{\delta}(\hat{\tau}_{1},\tau_{2})$ rely on the key Condition [2](https://arxiv.org/html/2510.16419v1#Thmcondition2) below.

###### Condition 2 .

$\mathbb{E}[|\tilde{\mu}_{a}(X)-\mu_{a}(X)||\tilde{e}(X)-e(X)|]=o_{\mathbb{P}}(n^{-1/2})$.

Condition [2](https://arxiv.org/html/2510.16419v1#Thmcondition2) is strictly weaker than Condition [1](https://arxiv.org/html/2510.16419v1#Thmcondition1). Moreover, the estimator $\hat{\delta}(\hat{\tau}_{1},\hat{\tau}_{2})$ offers several additional advantages over $\hat{\phi}(\hat{\tau})$, see Appendix [B](https://arxiv.org/html/2510.16419v1#A2) for more details.

Motivation. Although $\hat{\delta}(\hat{\tau}_{1},\hat{\tau}_{2})$ has several desirable properties, a notable limitation is that Condition [2](https://arxiv.org/html/2510.16419v1#Thmcondition2) requires all nuisance parameter estimators to be consistent (as $\tilde{e}(x)$ and $\tilde{\mu}_{a}(x)$ generally converge at most at the rate of $n^{-1/2}$), which may be too stringent for real-world applications.
In practice, the outcome regression model $\tilde{\mu}_{a}(x)$ is learned from the data with $A=a$ and then applied to the entire data. It heavily relies on model extrapolation, as there is often a significant distributional difference between the data with $A=a$ and $A=1-a$ (Jeong & Namkoong, 2020; Jing Qin & Huang, 2024).
As a result, $\tilde{\mu}_{a}(x)$ is likely to be inaccurate and biased, violating Assumption [2](https://arxiv.org/html/2510.16419v1#Thmcondition2). Therefore, it is beneficial and practical to develop methods that rely less on model extrapolation. In contrast, the estimation of the propensity score does not depend on extrapolation, making it less susceptible to this issue.

A natural and practical question arises: Can we develop a method for estimating $\delta(\hat{\tau}_{1},\hat{\tau}_{2})$, that retains all the desirable properties of $\hat{\delta}(\hat{\tau}_{1},\hat{\tau}_{2})$, while allowing for bias in $\tilde{\mu}_{a}(x)$ (relaxing Condition [2](https://arxiv.org/html/2510.16419v1#Thmcondition2))? In this article, we show that this is achievable by carefully exploiting the connection between the propensity score and outcome regression models, and by designing appropriate loss functions.

###### Condition 1 .

###### Condition 2 .

## 4 Proposed Method

In this section, we propose a novel method for estimating the relative error $\delta(\hat{\tau}_{1},\hat{\tau}_{2})$ that retains the desirable properties of $\hat{\delta}(\hat{\tau}_{1},\hat{\tau}_{2})$ while simultaneously being robust to bias in $\tilde{\mu}_{a}(x)$ for $a=0,1$.
We consider the following working models for the propensity score and outcome regression functions,

$$ $\displaystyle e(X)=\mathbb{P}(A=1\mid X)={}$ $\displaystyle e(\Phi(X),\gamma)=\frac{\exp(\Phi(X)^{\intercal}\gamma)}{1+\exp(\Phi(X)^{\intercal}\gamma)},$ (1) $\displaystyle\mu_{a}(X)=\mathbb{E}(Y\mid X,A=a)={}$ $\displaystyle\mu_{a}(\Phi(X),\beta_{a})=\Phi(X)^{\intercal}\beta_{a},\quad a=0,1,$ (2) $$

where $\Phi(X)$ is the a representation of $X$.

To quantify the bias of $\tilde{\mu}_{a}(x)$, it is crucial to distinguish between the working model and the true model. We say a working model is misspecified if the true model does not belong to the working model class, and it is correctly specified if the true model is within the working model class. Example [1](https://arxiv.org/html/2510.16419v1#Thmexample1) provides a misspecified example.

###### Example 1 (A misspecified model) .

Consider $X$ as a scalar and assume that the true model is $\mu_{a}(X)=\mathbb{E}(Y|X,A=a)=X^{2}\beta_{a}^{*}$, which represents the true data-generating mechanism of $Y$ given $(X,A=a)$. However, if we learn $\mathbb{E}(Y|X,A=a)$ using a linear model, i.e. $\mu_{a}(X,\beta_{a}):=X\beta_{a},\beta_{a}\in\mathbb{R}$, we introduce an inductive bias, meaning we can never reach the true value of $\beta^{*}$, even though the estimator may converge.
Specifically, denote $\hat{\beta}_{a}$ as the least-square estimator of $\beta_{a}$.
By the property of least-square estimator, it converges to $\bar{\beta}_{a}:=\mathbb{E}[XX]^{-1}\mathbb{E}[XY]$, regardless of whether $\mu_{a}(X,\beta_{a})$ is correctly specified or not.
Since $\mu_{a}(X,\beta_{a})$ is misspecified, $\bar{\beta}_{a}\neq\beta_{a}^{*}$.

For models ([1](https://arxiv.org/html/2510.16419v1#S4.E1)) and ([2](https://arxiv.org/html/2510.16419v1#S4.E2)), let $\check{\gamma}$ and $\check{\beta}_{a}$ denote the estimators of $\gamma$ and $\beta_{a}$, respectively. Define $\bar{\gamma}$ and $\bar{\beta}_{a}$ as the probability limits of $\check{\gamma}$ and $\check{\beta}_{a}$, and denote $\bar{e}(X)=e(\Phi(X),\bar{\gamma})$ and $\bar{\mu}_{a}(X)=\mu_{a}(\Phi(X),\bar{\beta}_{a})$.
If model ([1](https://arxiv.org/html/2510.16419v1#S4.E1)) is specified correctly, $e(X)=\bar{e}(X)$; otherwise, $e(X)\neq\bar{e}(X)$ and their difference represents the systematic bias induced by model misspecification.
Similarly, if model ([2](https://arxiv.org/html/2510.16419v1#S4.E2)) is correctly specified, $\bar{\mu}_{a}(X)=\mu_{a}(X)$; otherwise, $\bar{\mu}_{a}(X)\neq\mu_{a}(X)$.
*It is important to note that $(\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})$ always converges to $(\bar{\gamma},\bar{\beta}_{0},\bar{\beta}_{1})$, regardless of whether models ([1](https://arxiv.org/html/2510.16419v1#S4.E1)) and ([2](https://arxiv.org/html/2510.16419v1#S4.E2)) are correctly specified.*

###### Example 1 (A misspecified model) .

### 4.1 Basic Idea

Before delving into the details, we outline its basic idea to provide an intuitive understanding.

First, to retain the semiparametric efficiency, the proposed estimator $\delta(\hat{\tau}_{1},\hat{\tau}_{2})$ preserves the same form of $\hat{\delta}(\hat{\tau}_{1},\tau_{2})$, which is given as

$$ $\displaystyle\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})={}\frac{1}{n}\sum_{i=1}^{n}\varphi(Z_{i};\check{\mu}_{0},\check{\mu}_{1},\check{e}),$ $$

where $\check{e}(X)=e(\Phi(X),\check{\gamma})$ and $\check{\mu}_{a}(X)=\mu_{a}(\Phi(X),\check{\beta}_{a})$ for $a=0,1$.
Although $\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})$ and $\delta(\hat{\tau}_{1},\hat{\tau}_{2})$ share the same form, they differ significantly in how they estimate nuisance parameters, which results in the robustness to biases in $\check{\mu}_{a}(X)$.

Second, we analyze the key conditions necessary to achieve robustness to biases in $\check{\mu}_{a}(X)$.
By a Taylor expansion of $\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})$ with respect to $(\bar{\gamma},\bar{\beta}_{0},\bar{\beta}_{1})$, we have that

$$ $\displaystyle\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})$ $\displaystyle-\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\bar{\gamma},\bar{\beta}_{0},\bar{\beta}_{1})=\Delta_{\gamma}^{\intercal}(\check{\gamma}-\bar{\gamma})+\Delta_{\beta_{0}}^{\intercal}(\check{\beta}_{0}-\bar{\beta}_{0})+\Delta_{\beta_{1}}^{\intercal}(\check{\beta}_{1}-\bar{\beta}_{1})$ $\displaystyle+O_{\mathbb{P}}((\check{\gamma}-\bar{\gamma})^{2}+(\check{\gamma}-\bar{\gamma})(\check{\beta}_{1}-\bar{\beta}_{1})+(\check{\gamma}-\bar{\gamma})(\check{\beta}_{0}-\bar{\beta}_{0})),$ $$

where

$$ $\displaystyle\Delta_{\gamma}={}$ $\displaystyle-\frac{1}{n}\sum_{i=1}^{n}2(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\Big(\frac{A_{i}(1-\bar{e}(X_{i}))(Y_{i}-\bar{\mu}_{1}(X_{i}))}{\bar{e}(X_{i})}+\frac{(1-A_{i})\bar{e}(X_{i})(Y_{i}-\bar{\mu}_{0}(X_{i}))}{1-\bar{e}(X_{i})}\Big)\Phi(X_{i}),$ $\displaystyle\Delta_{\beta_{0}}={}$ $\displaystyle-\frac{1}{n}\sum_{i=1}^{n}2(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\left(1-\frac{1-A_{i}}{1-\bar{e}(X_{i})}\right)\Phi(X_{i}),$ $\displaystyle\Delta_{\beta_{1}}={}$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}2(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\left(1-\frac{A_{i}}{\bar{e}(X_{i})}\right)\Phi(X_{i}).$ $$

Under mild conditions (see Theorem [1](https://arxiv.org/html/2510.16419v1#Thmtheorem1)),
the last term of above Taylor expansion is $o_{\mathbb{P}}(n^{-1/2})$. We note that $\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\bar{\gamma},\bar{\beta}_{0},\bar{\beta}_{1})$ is a $\sqrt{n}$-consistent and asymptotically normal estimator of $\delta(\hat{\tau}_{1},\hat{\tau}_{2})$ if either $\check{e}(x)$ is correctly specified, or $(\check{\mu}_{0}(x),\check{\mu}_{1}(x))$ is correctly specified. Thus, it is robust to biases in $\check{\mu}_{a}(x)$ for $a=0,1$ and is the ideal estimator we aim to obtain.
To ensure that $\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})$ has the same asymptotic properties as $\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\bar{\gamma},\bar{\beta}_{0},\bar{\beta}_{1})$, we require that

$$ $\Delta_{\gamma}^{\intercal}(\check{\gamma}-\bar{\gamma})+\Delta_{\beta_{0}}^{\intercal}(\check{\beta}_{0}-\bar{\beta}_{0})+\Delta_{\beta_{1}}^{\intercal}(\check{\beta}_{1}-\bar{\beta}_{1})=o_{\mathbb{P}}(n^{-1/2}),$ (3) $$

*even when $(\check{\mu}_{0}(x),\check{\mu}_{1}(x))$ is misspecified.*
Note that $(\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})$ always converges to $(\bar{\gamma},\bar{\beta}_{0},\bar{\beta}_{1})$.
To satisfy Eq. ([3](https://arxiv.org/html/2510.16419v1#S4.E3)), it suffices for $\Delta_{\gamma}$, $\Delta_{\beta_{0}}$, and $\Delta_{\beta_{1}}$ to converge to zero at a certain rate. By the central limit theorem, $\Delta_{\gamma}-\mathbb{E}[\Delta_{\gamma}]=O_{\mathbb{P}}(n^{-1/2})$, $\Delta_{\beta_{0}}-\mathbb{E}[\Delta_{\beta_{0}}]=O_{\mathbb{P}}(n^{-1/2})$, and $\Delta_{\beta_{1}}-\mathbb{E}[\Delta_{\beta_{1}}]=O_{\mathbb{P}}(n^{-1/2})$. Thus, Eq. ([3](https://arxiv.org/html/2510.16419v1#S4.E3)) holds provided that

$$ $\mathbb{E}[\Delta_{\gamma}]=0,~\mathbb{E}[\Delta_{\beta_{0}}]=0,~\mathbb{E}[\Delta_{\beta_{1}}]=0,$ $$

which is equivalent to the following equations:

$$ $\displaystyle\begin{cases}&\mathbb{E}\left[(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\left(\frac{A_{i}(1-\bar{e}(X_{i}))(Y_{i}-\bar{\mu}_{1}(X_{i}))}{\bar{e}(X_{i})}+\frac{(1-A_{i})\bar{e}(X_{i})(Y_{i}-\bar{\mu}_{0}(X_{i}))}{1-\bar{e}(X_{i})}\right)\Phi(X_{i})\right]=0,\\ &\mathbb{E}\left[(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\left(1-\frac{A_{i}}{\bar{e}(X_{i})}\right)\Phi(X_{i})\right]=0,\\ &\mathbb{E}\left[(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\left(1-\frac{1-A_{i}}{1-\bar{e}(X_{i})}\right)\Phi(X_{i})\right]=0.\end{cases}$ (4) $$

### 4.2 Novel Loss for Nuisance Parameter Estimation

To ensure that the first term in Eq. ([4](https://arxiv.org/html/2510.16419v1#S4.E4)) holds, we design the weighted least square loss function for $(\beta_{0},\beta_{1}$) as follows:

$$ $\displaystyle\mathcal{L}_{\text{wls}}(\beta_{0},\beta_{1};\check{\gamma})={}\frac{1}{n}\sum_{i=1}^{n}(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\left[\frac{(1-A_{i})\check{e}(X_{i})\{Y_{i}-\Phi(X)^{\intercal}\beta_{0}\}^{2}}{1-\check{e}(X_{i})}+\frac{A_{i}(1-\hat{e}(X_{i}))\{Y_{i}-\Phi(X)^{\intercal}\beta_{1}\}^{2}}{\hat{e}(X_{i})}\right].$ $$

These loss functions imply that $(\bar{\beta}_{0},\bar{\beta}_{1})\triangleq\arg\min_{\beta_{a}}\mathbb{E}[\mathcal{L}_{\text{wls}}(\beta_{0},\beta_{1};\bar{\gamma})]$.
By setting $\partial\mathbb{E}[\mathcal{L}_{\text{wls}}(\beta_{0},\beta_{1};\bar{\gamma})]/\partial\beta_{0}\big|_{\beta_{0}=\bar{\beta}_{0}}=0$ and $\partial\mathbb{E}[\mathcal{L}_{\text{wls}}(\beta_{0},\beta_{1};\bar{\gamma})]/\partial\beta_{1}\big|_{\beta_{1}=\bar{\beta}_{1}}=0$, one can see that the first term in Eq. ([4](https://arxiv.org/html/2510.16419v1#S4.E4)) holds even if $(\check{\mu}_{0}(x),\check{\mu}_{1}(x))$ is misspecified.

For learning $\gamma$, note that Eq. ([4](https://arxiv.org/html/2510.16419v1#S4.E4)) imposes $2d$ linear constraints, while $\gamma\in\mathbb{R}^{d}$ has only $d$ degrees of freedom. This makes the system over-constrained. To address this, following the soft-margin formulation of support vector machines (Murphy, 2022), we introduce slack variables $\xi,\eta\in\mathbb{R}^{d}$ to allow controlled constraint violations, and penalize their magnitudes in the objective.
Formally, we solve:

$$ $\displaystyle\min_{\gamma,\xi,\eta}\quad$ $\displaystyle-\frac{1}{n}\sum_{i=1}^{n}\left[A_{i}\log(e(X_{i}))+(1-A_{i})\log(1-e(X_{i}))\right]+c\sum_{j=1}^{d}(\xi_{j}+\eta_{j})$ s.t. $\displaystyle e(X_{i})=\frac{\exp(\Phi(X_{i})^{\top}\gamma)}{1+\exp(\Phi(X_{i})^{\top}\gamma)},\quad i=1,\dots,n,$ $\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\left(1-\frac{A_{i}}{e(X_{i})}\right)\Phi_{j}(X_{i})\right|\leq\xi_{j},\quad\forall j,$ $\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\left(1-\frac{1-A_{i}}{1-e(X_{i})}\right)\Phi_{j}(X_{i})\right|\leq\eta_{j},\quad\forall j,$ $\displaystyle\xi_{j},\eta_{j}\geq 0,\quad j=1,\dots,d.$ $$

where $c$ is a given hyperparameter.
In practice, we convert the above constrained optimization into two unconstrained loss terms:

$$ $\mathcal{L}_{\text{ce}}=-\frac{1}{n}\sum_{i=1}^{n}\left[A_{i}\log(e(X_{i}))+(1-A_{i})\log(1-e(X_{i}))\right],$ $$

$$ $\displaystyle\mathcal{L}_{\text{const}}=c\sum_{j=1}^{d}(\xi_{j}+\eta_{j})+\rho\cdot\left\|\begin{bmatrix}\max\left\{\left|\frac{1}{n}\sum_{i=1}^{n}(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\left(1-\frac{A_{i}}{e(X_{i})}\right)\Phi(X_{i})\right|-\xi,\;0\right\}\\ \max\left\{\left|\frac{1}{n}\sum_{i=1}^{n}(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\left(1-\frac{1-A_{i}}{1-e(X_{i})}\right)\Phi(X_{i})\right|-\eta,\;0\right\}\\ \max(-\xi,\;0)\\ \max(-\eta,\;0)\end{bmatrix}\right\|_{2},$ $$

where $\rho>0$ is a penalty parameter encouraging constraint satisfaction.

### 4.3 Constructing Neural Network

Building on the novel constraint loss introduced in Section [4.2](https://arxiv.org/html/2510.16419v1#S4.SS2), we propose a new neural network architecture, inspired by the Dragonnet structure (Shi et al., 2019a).
The proposed network takes input features $x\in\mathbb{R}^{d}$, and first passes them through multiple fully connected layers to produce the shared representation $\Phi(x)\in\mathbb{R}^{m}$. This representation is then fed into three separate heads: a control outcome head $\mu_{0}(x)$, predicting the potential outcome under control; a treated outcome head $\mu_{1}(x)$, predicting the potential outcome under treatment; a treatment head $e(x)$, estimating the propensity score via a sigmoid activation.

The control outcome head and the treated outcome head contribute to the weighted least square loss $\mathcal{L}_{\text{wls}}$, while $\mathcal{L}_{\text{ce}}$ and $\mathcal{L}_{\text{const}}$ are computed by the treatment head and the shared representation. During training, we minimize the total training loss given by:

$$ $\mathcal{L}=\mathcal{L}_{\mathrm{wls}}+\lambda_{1}\mathcal{L}_{\mathrm{ce}}+\lambda_{2}\mathcal{L}_{\mathrm{const}}.$ $$

This formulation encourages the propensity model $e(X)$ and the outcome model $\mu_{a}(X)$ to satisfy Eq. ([4](https://arxiv.org/html/2510.16419v1#S4.E4)), providing a reliable estimation that can be used in computing the estimated relative error $\hat{\delta}$ mentioned in Section [3](https://arxiv.org/html/2510.16419v1#S3). For clarity, we provide a schematic illustration of the network architecture in Appendix [C](https://arxiv.org/html/2510.16419v1#A3).

### 4.4 Theoretical Analysis

We analyze the large sample properties of the proposed estimator $\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})$.

###### Theorem 1 .

If the propensity score model is correctly specified, and $\check{\gamma}$, $\check{\beta}_{0}$ as well as $\check{\beta}_{1}$ converge to their probability limits at a rate faster than $n^{-1/4}$, then we have

$$ $\displaystyle\sqrt{n}\{\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})-\delta(\hat{\tau}_{1},\hat{\tau}_{2})\}\xrightarrow{d}\mathcal{N}(0,\sigma^{2}),$ $$

where $\sigma^{2}=\text{Var}\{\varphi(Z;\bar{u}_{0},\bar{u}_{1},\bar{e})\}$ and $\xrightarrow{d}$ means convergence in distribution.

Theorem [1](https://arxiv.org/html/2510.16419v1#Thmtheorem1) shows that the proposed estimator
is $\sqrt{n}$-consistent and asymptotically normal.
These properties hold even when the outcome regression model is misspecified, as long as $\check{\gamma}$, $\check{\beta}_{0}$, and $\check{\beta}_{1}$ converge to their respective probability limits at a rate faster than $n^{-1/4}$. This condition is readily satisfied, as $(\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})$ always converge to their probability limits $(\bar{\gamma},\bar{\beta}_{0},\bar{\beta}_{1})$, and a variety of flexible machine learning methods can achieve the required convergence rates (Chernozhukov et al., 2018; Semenova & Chernozhukov, 2021).

Based on Theorem [1](https://arxiv.org/html/2510.16419v1#Thmtheorem1), we can obtain a valid asymptotic $(1-\eta)$ confidence interval of $\delta(\hat{\tau}_{1},\hat{\tau}_{2})$.

###### Proposition 2 .

Under the conditions in Theorem [1](https://arxiv.org/html/2510.16419v1#Thmtheorem1),
a consistent estimator of $\sigma^{2}$ is

$$ $\displaystyle\hat{\sigma}^{2}=\frac{1}{n}\sum^{n}_{i=1}\left\{\varphi(Z_{i};\check{u}_{0},\check{u}_{1},\check{e})-\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})\right\}^{2},$ $$

an asymptotic $(1-\eta)$ confidence interval for $\delta(\hat{\tau}_{1},\hat{\tau}_{2})$ is $\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})\pm z_{\eta/2}\sqrt{\hat{\sigma}^{2}/n}$, where $z_{\eta/2}$ is the $(1-\eta/2)$ quantile of the standard normal distribution.

Proposition [2](https://arxiv.org/html/2510.16419v1#Thmtheorem2) shows that a valid asymptotic confidence interval for $\delta(\hat{\tau}_{1},\hat{\tau}_{2})$ is achievable even with a misspecified outcome model, unlike previous methods that require correct specification.
This further indicates the robustness of the proposed method.

###### Theorem 1 .

###### Proposition 2 .

## 5 Enhanced Estimation of Heterogeneous Treatment Effects

In this section, building on the evaluation framework proposed in Section [4](https://arxiv.org/html/2510.16419v1#S4), we extend the idea to develop a learning method for CATE.
In general, a reliable evaluation method can naturally serve as a basis for developing a learning method.
In our proposed approach, for any given pair of CATE estimators $\hat{\tau}_{k}(x)$ and $\hat{\tau}_{k^{\prime}}(x)$, the proposed neural network architecture introduced in Section [4.3](https://arxiv.org/html/2510.16419v1#S4.SS3) can output the corresponding estimates of outcome regression functions. We denote them as $\check{\mu}_{0}(x;\hat{\tau}_{k},\hat{\tau}_{k^{\prime}})$ and $\check{\mu}_{1}(x;\hat{\tau}_{k},\hat{\tau}_{k^{\prime}})$, emphasizing their dependence on $\hat{\tau}_{k}(x)$ and $\hat{\tau}_{k^{\prime}}(x)$. This leads to a new CATE estimator, defined as

$$ $\check{\tau}(x;\hat{\tau}_{k},\hat{\tau}_{k^{\prime}})=\check{\mu}_{1}(x;\hat{\tau}_{k},\hat{\tau}_{k^{\prime}})-\check{\mu}_{0}(x;\hat{\tau}_{k},\hat{\tau}_{k^{\prime}}).$ $$

Clearly, the performance of the estimator heavily depends on the choice of CATE estimators $\check{\tau}(x;\hat{\tau}_{k},\hat{\tau}_{k^{\prime}})$. However, due to the fundamental challenge in evaluating CATE (i.e., the absence of ground truth), it is difficult to develop a direct strategy for selecting them.
To mitigate this issue, we propose the following aggregation strategy for estimating CATE,

$$ $\check{\tau}(x)=\frac{2}{|\mathcal{K}|(|\mathcal{K}|-1)}\sum_{k,k^{\prime}\in\mathcal{K}}\check{\mu}_{1}(x;\hat{\tau}_{k},\hat{\tau}_{k^{\prime}})-\check{\mu}_{0}(x;\hat{\tau}_{k},\hat{\tau}_{k^{\prime}})$ $$

where $\mathcal{K}=\{1,2,\ldots,K\}$ is the index set for the candidate CATE estimators. This aggregated estimator aims to stabilize and improve the estimation of CATE by averaging over all pairs of candidate estimators. When $K$ is large, averaging over all pairs can be computationally burdensome. In such cases, one can randomly select a subset of pairs $\check{\tau}(x;\hat{\tau}_{k},\hat{\tau}_{k^{\prime}})$ and compute their average instead.
Surprisingly, our experiments show that this estimator performs exceptionally well, even surpassing the performance of any single candidate estimator.

## 6 Experiments

### 6.1 Experimental Setup

Datasets and Processing. Following previous studies (Yoon et al., 2018; Yao et al., 2018; Louizos et al., 2017), we choose one semi-synthetic dataset IHDP, and two real datasets, Twins and Jobs, to conduct our experiments. The Twins dataset is constructed from all twin births in the United States between 1989 and 1991 (Almond et al., 2005), owning 5271 samples with 28 different covariates.
The IHDP dataset is used to estimate the effect of specialist home visits on infants’ future cognitive test scores, containing 747 samples (139 treated and 608 control), each with 25 pre-treatment covariates, while the Jobs dataset focuses on estimating the impact of job training programs on individuals’ employment status, including 297 treated units, 425 control units from the experimental sample, and 2490 control units from the observational sample. We provide more dataset details in the Appendix [E.1](https://arxiv.org/html/2510.16419v1#A5.SS1). We randomly split each dataset into training and test sets in a 2:1 ratio, and repeat the experiments 50 times for the Twins, 100 times for the IHDP, and 20 times for the Jobs.

Evaluation Metrics.
We consider two classes of evaluation metrics below.

- •
We assess the proposed relative error estimator using two key metrics:
(i) the coverage probability of its confidence interval (named coverage rate), and
(ii) the probability of correctly identifying the better estimator (i.e., selecting the true winner, named selection accuracy).
In practice, we only pick the winner when the confidence interval for the relative error does not contain zero, otherwise, no selection will be made. We calculate the coverage rate of the targeted 90% confidence intervals and selection accuracy.
- •
For evaluating the performance of CATE estimation of our novel network, following previous studies (Shalit et al., 2017a; Shi et al., 2019b; Louizos et al., 2017), we compute the Precision in Estimation of Heterogeneous Effect (PEHE) (and, 2011), where $\sqrt{\epsilon_{\text{PEHE}}}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}\left(\hat{\tau}(x_{i})-\tau(x_{i})\right)^{2}}$, and
the absolute error on the ATE, $\epsilon_{\text{ATE}}=|\text{ATE}-\widehat{\text{ATE}}|$ ,
where
$\text{ATE}=\frac{1}{n}\sum_{i=1}^{n}(y_{i}^{1}-y_{j}^{0}),$
in which $y_{i}^{1}$ and $y_{j}^{0}$ are the true potential outcomes.

Figure: Figure 1: Coverage rate on IHDP and Twins.
Refer to caption: coverage.png

**Table 1: CATE estimation performance on the IHDP and Twins datasets (in-sample and out-of-sample). The best results are bolded.**
|  | IHDP | Twins |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Method | $\sqrt{\epsilon_{\text{PEHE}}^{\text{in}}}$ | $\epsilon_{\text{ATE}}^{\text{in}}$ | $\sqrt{\epsilon_{\text{PEHE}}^{\text{out}}}$ | $\epsilon_{\text{ATE}}^{\text{out}}$ | $\sqrt{\epsilon_{\text{PEHE}}^{\text{in}}}$ | $\epsilon_{\text{ATE}}^{\text{in}}$ | $\sqrt{\epsilon_{\text{PEHE}}^{\text{out}}}$ | $\epsilon_{\text{ATE}}^{\text{out}}$ |
| LinDML | 1.053 $\pm$ 0.134 | 0.580 $\pm$ 0.152 | 1.085 $\pm$ 0.187 | 0.574 $\pm$ 0.176 | 0.295 $\pm$ 0.005 | 0.013 $\pm$ 0.009 | 0.296 $\pm$ 0.008 | 0.013 $\pm$ 0.010 |
| SpaDML | 0.832 $\pm$ 0.119 | 0.252 $\pm$ 0.185 | 0.866 $\pm$ 0.112 | 0.280 $\pm$ 0.183 | 0.300 $\pm$ 0.008 | 0.046 $\pm$ 0.030 | 0.303 $\pm$ 0.010 | 0.046 $\pm$ 0.033 |
| CForest | 0.891 $\pm$ 0.121 | 0.419 $\pm$ 0.182 | 0.903 $\pm$ 0.127 | 0.403 $\pm$ 0.185 | 0.297 $\pm$ 0.005 | 0.012 $\pm$ 0.008 | 0.306 $\pm$ 0.008 | 0.013 $\pm$ 0.011 |
| X-Learner | 0.971 $\pm$ 0.178 | 0.196 $\pm$ 0.137 | 0.987 $\pm$ 0.196 | 0.207 $\pm$ 0.141 | 0.293 $\pm$ 0.005 | 0.022 $\pm$ 0.014 | 0.294 $\pm$ 0.008 | 0.024 $\pm$ 0.016 |
| S-Learner | 0.920 $\pm$ 0.102 | 0.212 $\pm$ 0.100 | 0.950 $\pm$ 0.111 | 0.205 $\pm$ 0.117 | 0.298 $\pm$ 0.011 | 0.057 $\pm$ 0.042 | 0.299 $\pm$ 0.010 | 0.059 $\pm$ 0.042 |
| TARNet | 0.896 $\pm$ 0.054 | 0.279 $\pm$ 0.084 | 0.920 $\pm$ 0.070 | 0.266 $\pm$ 0.117 | 0.292 $\pm$ 0.011 | 0.090 $\pm$ 0.047 | 0.294 $\pm$ 0.019 | 0.091 $\pm$ 0.045 |
| Dragonnet | 0.840 $\pm$ 0.046 | 0.124 $\pm$ 0.089 | 0.867 $\pm$ 0.087 | 0.134 $\pm$ 0.092 | 0.292 $\pm$ 0.004 | 0.080 $\pm$ 0.008 | 0.290 $\pm$ 0.007 | 0.092 $\pm$ 0.011 |
| DRCFR | 0.741 $\pm$ 0.068 | 0.186 $\pm$ 0.138 | 0.760 $\pm$ 0.090 | 0.185 $\pm$ 0.135 | 0.290 $\pm$ 0.004 | 0.075 $\pm$ 0.007 | 0.288 $\pm$ 0.007 | 0.076 $\pm$ 0.010 |
| SCIGAN | 0.898 $\pm$ 0.374 | 0.358 $\pm$ 0.509 | 0.919 $\pm$ 0.369 | 0.358 $\pm$ 0.502 | 0.296 $\pm$ 0.037 | 0.041 $\pm$ 0.044 | 0.293 $\pm$ 0.039 | 0.040 $\pm$ 0.047 |
| DESCN | 0.793 $\pm$ 0.187 | 0.133 $\pm$ 0.106 | 0.835 $\pm$ 0.197 | 0.140 $\pm$ 0.112 | 0.296 $\pm$ 0.060 | 0.059 $\pm$ 0.043 | 0.293 $\pm$ 0.063 | 0.058 $\pm$ 0.042 |
| ESCFR | 0.802 $\pm$ 0.041 | 0.111 $\pm$ 0.070 | 0.841 $\pm$ 0.074 | 0.135 $\pm$ 0.076 | 0.290 $\pm$ 0.004 | 0.075 $\pm$ 0.007 | 0.288 $\pm$ 0.007 | 0.076 $\pm$ 0.010 |
| Ours | 0.638 $\pm$ 0.138 | 0.090 $\pm$ 0.087 | 0.670 $\pm$ 0.150 | 0.105 $\pm$ 0.099 | 0.284 $\pm$ 0.005 | 0.009 $\pm$ 0.005 | 0.286 $\pm$ 0.007 | 0.009 $\pm$ 0.006 |

Baselines and Experimental Details.
To evaluate the performance of relative error estimation, we select three representative estimators from different methodological families: Causal Forest (tree-based) (Athey & Wager, 2019), X-Learner (meta-learner) (Künzel et al., 2019), and TARNet (representation learning) (Shalit et al., 2017a). We estimate their pairwise relative errors and evaluate the estimation performances. Although Gao’s work does not propose a concrete learning method, we follow their choice of nuisance estimators (Linear Regression, Boosting) to compute relative errors for reference (see Appendix [E.2](https://arxiv.org/html/2510.16419v1#A5.SS2)).

For CATE estimation, the baselines include Causal Forest (Athey & Wager, 2019), meta-learners (X-Learner, S-Learner) (Künzel et al., 2019), double machine learning (Linear DML, Sparse Linear DML) (Chernozhukov et al., 2024), TARNet (Shalit et al., 2017a), Dragonnet (Shi et al., 2019a), DR-CFR (Hassanpour & Greiner, 2020), SCIGAN (Bica et al., 2020), DESCN (Zhong et al., 2022) and ESCFR (Wang et al., 2023). In addition, see Appendix [E.5](https://arxiv.org/html/2510.16419v1#A5.SS5) for training details of hyperparameter tuning range.

### 6.2 Experimental Results

Quality of Relative Error Estimation. We first evaluate the performance of relative error estimation, comparing different pairs of HTE estimators.
In Figures [2](https://arxiv.org/html/2510.16419v1#S6.F2) and [2](https://arxiv.org/html/2510.16419v1#S6.F2), we present the coverage of the 90% confidence intervals
and the accuracy of selecting the better HTE estimator on the test sets, respectively, where TN stands for TARNet, CF stands for Causal Forest, and X stands for X-Learner, and the red dashed line marks the target level of 90%.
From these two figures, our method successfully achieves the target coverage, and
provide trustworthy advice on the selection across different pairs of HTE estimators.
These results demonstrate the validity of our uncertainty quantification and estimator selection.

Accuracy of the CATE Estimation. We then evaluate the performance of CATE estimation learned by our novel network and compare it with competing baselines.
We average over 100 realizations of our networks in IHDP and 50 realizations in Twins. The results are presented in Table [1](https://arxiv.org/html/2510.16419v1#S6.T1). Our proposed method achieves the best performance across all metrics, with the lowest $\sqrt{\epsilon_{\text{PEHE}}}$ and $\epsilon_{\text{ATE}}$ on both datasets. This demonstrates its ability to accurately estimate CATE. In addition, we report results on the Jobs dataset in Appendix [E.3](https://arxiv.org/html/2510.16419v1#A5.SS3) due to limited space.

**Table 2: Sensitivity analysis on the hyperparameter $\lambda_{2}$ (weight of constraint loss) for IHDP and Twins datasets. The best hyperparameter values and results are in bold.**
| IHDP | Twins |  |  |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Value | $\sqrt{\epsilon_{\text{PEHE}}^{\text{in}}}$ | $\epsilon_{\text{ATE}}^{\text{in}}$ | $\sqrt{\epsilon_{\text{PEHE}}^{\text{out}}}$ | $\epsilon_{\text{ATE}}^{\text{out}}$ | Coverage | Selection | Value | $\sqrt{\epsilon_{\text{PEHE}}^{\text{in}}}$ | $\epsilon_{\text{ATE}}^{\text{in}}$ | $\sqrt{\epsilon_{\text{PEHE}}^{\text{out}}}$ | $\epsilon_{\text{ATE}}^{\text{out}}$ | Coverage | Selection |
| 0.01 | 0.860 | 0.216 | 0.902 | 0.238 | 0.85 | 0.50 | 0.005 | 0.319 | 0.029 | 0.331 | 0.027 | 0.82 | 0.38 |
| 0.1 | 0.800 | 0.142 | 0.837 | 0.158 | 0.91 | 0.61 | 0.05 | 0.289 | 0.016 | 0.292 | 0.015 | 0.82 | 0.84 |
| 0.5 | 0.714 | 0.099 | 0.747 | 0.118 | 0.95 | 0.78 | 0.25 | 0.297 | 0.018 | 0.297 | 0.020 | 0.86 | 0.42 |
| 1 | 0.638 | 0.090 | 0.670 | 0.105 | 0.96 | 0.80 | 0.5 | 0.284 | 0.009 | 0.286 | 0.009 | 0.94 | 0.94 |
| 5 | 0.715 | 0.099 | 0.748 | 0.116 | 0.94 | 0.77 | 2.5 | 0.285 | 0.011 | 0.287 | 0.012 | 0.94 | 0.92 |
| 10 | 0.795 | 0.157 | 0.830 | 0.172 | 0.90 | 0.60 | 5 | 0.289 | 0.028 | 0.290 | 0.026 | 0.80 | 0.86 |
| 100 | 0.801 | 0.156 | 0.836 | 0.170 | 0.90 | 0.60 | 50 | 0.287 | 0.024 | 0.288 | 0.023 | 0.84 | 0.88 |

**Table 3: Ablation study results on the IHDP and Twins datasets.**
|  | IHDP | Twins |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Training Loss | $\sqrt{\epsilon_{\text{PEHE}}^{\text{in}}}$ | $\epsilon_{\text{ATE}}^{\text{in}}$ | $\sqrt{\epsilon_{\text{PEHE}}^{\text{out}}}$ | $\epsilon_{\text{ATE}}^{\text{out}}$ | Coverage | Selection | $\sqrt{\epsilon_{\text{PEHE}}^{\text{in}}}$ | $\epsilon_{\text{ATE}}^{\text{in}}$ | $\sqrt{\epsilon_{\text{PEHE}}^{\text{out}}}$ | $\epsilon_{\text{ATE}}^{\text{out}}$ | Coverage | Selection |
| $\mathcal{L}_{\text{wls}}$ & $\mathcal{L}_{\text{const}}$ | 0.725 | 0.101 | 0.758 | 0.122 | 0.92 | 0.71 | 0.284 | 0.013 | 0.287 | 0.013 | 0.94 | 0.92 |
| $\mathcal{L}_{\text{wls}}$ & $\mathcal{L}_{\text{ce}}$ | 3.495 | 2.879 | 3.531 | 2.900 | 0.88 | 0.14 | 0.319 | 0.028 | 0.328 | 0.026 | 0.82 | 0.14 |
| Full (Ours) | 0.638 | 0.090 | 0.670 | 0.105 | 0.96 | 0.80 | 0.284 | 0.009 | 0.286 | 0.009 | 0.94 | 0.94 |

Sensitive Analysis. The hyperparameter $\lambda_{2}$ before the constraint loss $\mathcal{L}_{\text{const}}$,
$\lambda_{1}$ before the cross entropy loss $\mathcal{L}_{\text{ce}}$ and the penalty weight $\rho$ in the constraint loss
play important roles in our training. In order to explore under which parameters our method has the best performance, we conduct sensitivity analysis experiments. We present
the results of $\lambda_{2}$ in Table [2](https://arxiv.org/html/2510.16419v1#S6.T2).
We observe that both the performance of CATE estimation and the relative error estimation remain relatively stable across a range of $\lambda_{2}$ values from 0.5 to 5, indicating robustness to this hyperparameter.
However, when $\lambda_{2}$ is extremely small (e.g., $\lambda_{2}=0.01$), the performance of the proposed method degrades significantly, indicating the importance of the constraint loss $\mathcal{L}_{\text{const}}$.
Also, we perform sensitive analysis for $\lambda_{1}$ and $\rho$, the associated results are provided in the Appendix [E.4](https://arxiv.org/html/2510.16419v1#A5.SS4).

Ablation Study. As shown in Section [4.3](https://arxiv.org/html/2510.16419v1#S4.SS3), the proposed method involves three loss functions: $\mathcal{L}_{\mathrm{wls}},\mathcal{L}_{\mathrm{ce}}$, and $\mathcal{L}_{\mathrm{const}}$.
We conduct an ablation study to assess the impact of $\mathcal{L}_{\mathrm{ce}}$ and $\mathcal{L}_{\mathrm{const}}$ on overall performance. The corresponding results are reported in Table [3](https://arxiv.org/html/2510.16419v1#S6.T3).
Specifically, removing $\mathcal{L}_{\text{const}}$ results in a notable drop in the accuracy of both outcome and relative error estimation, whereas removing $\mathcal{L}_{\text{ce}}$ only causes a moderate decline. These findings highlight the importance of the proposed novel loss $\mathcal{L}_{\text{const}}$, which not only improves HTE estimation accuracy but also facilitates the construction of narrower and more precise confidence intervals for relative error.

## 7 Conclusion

In this work, we addressed a key challenge in evaluating HTE estimators with less reliance on modeling assumptions for nuisance parameters. Building upon the relative error framework, we introduced a novel loss function and balance regularizers that encourage more stable and accurate learning of nuisance parameters. These components were integrated into a new neural network architecture tailored to enhance the reliability of HTE evaluation.
The proposed evaluation approach retains several desirable statistical properties while relaxing the stringent requirement for consistent outcome regression models, thereby facilitating more reliable comparisons and selection of estimators in real-world applications. A limitation of this work lies in the use of the simple averaging scheme over all estimator pairs for CATE estimation. While this approach improves stability, it may not fully exploit the varying strengths of individual estimators, potentially limiting overall efficiency and precision. Future research is warranted to further address this challenge.

## References

- A. Smith & E. Todd (2005)
Jeffrey A. Smith and Petra E. Todd.
Does matching overcome lalonde’s critique of nonexperimental
estimators?
*Journal of Econometrics*, 125(1):305–353,
2005.
ISSN 0304-4076.
doi: https://doi.org/10.1016/j.jeconom.2004.04.011.
URL
[https://www.sciencedirect.com/science/article/pii/S030440760400082X](https://www.sciencedirect.com/science/article/pii/S030440760400082X).
Experimental and non-experimental evaluation of economic policy and
models.
- Almond et al. (2005)
Douglas Almond, Kenneth Y. Chay, and David S. Lee.
The costs of low birth weight*.
*The Quarterly Journal of Economics*, 120(3):1031–1083, 08 2005.
ISSN 0033-5533.
doi: 10.1093/qje/120.3.1031.
URL [https://doi.org/10.1093/qje/120.3.1031](https://doi.org/10.1093/qje/120.3.1031).
- and (2011)
Jennifer L. Hill and.
Bayesian nonparametric modeling for causal inference.
*Journal of Computational and Graphical Statistics*, 20(1):217–240, 2011.
doi: 10.1198/jcgs.2010.08162.
URL [https://doi.org/10.1198/jcgs.2010.08162](https://doi.org/10.1198/jcgs.2010.08162).
- Athey & Wager (2019)
Susan Athey and Stefan Wager.
Estimating treatment effects with causal forests: An application,
2019.
URL [https://arxiv.org/abs/1902.07409](https://arxiv.org/abs/1902.07409).
- Bica et al. (2020)
Ioana Bica, James Jordon, and Mihaela van der Schaar.
Estimating the effects of continuous-valued interventions using
generative adversarial networks.
*CoRR*, abs/2002.12326, 2020.
URL [https://arxiv.org/abs/2002.12326](https://arxiv.org/abs/2002.12326).
- Caron et al. (2022)
Alberto Caron, Gianluca Baio, and Ioanna Manolopoulou.
Estimating individual treatment effects using non-parametric
regression models: A review.
*Journal of the Royal Statistical Society: Series A (Statistics
in Society)*, 185:1115–1149, 2022.
- Chernozhukov et al. (2018)
V. Chernozhukov, D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, and
J. Robins.
Double/debiased machine learning for treatment and structural
parameters.
*The Econometrics Journal*, 21:1–68, 2018.
- Chernozhukov et al. (2024)
Victor Chernozhukov, Denis Chetverikov, Mert Demirer, Esther Duflo, Christian
Hansen, Whitney Newey, and James Robins.
Double/debiased machine learning for treatment and causal parameters,
2024.
URL [https://arxiv.org/abs/1608.00060](https://arxiv.org/abs/1608.00060).
- Curth & Van Der Schaar (2023)
Alicia Curth and Mihaela Van Der Schaar.
In search of insights, not magic bullets: towards demystification of
the model selection dilemma in heterogeneous treatment effect estimation.
In *Proceedings of the 40th International Conference on Machine
Learning*, ICML’23. JMLR, 2023.
- Dehejia & Wahba (2002)
Rajeev H. Dehejia and Sadek Wahba.
Propensity score-matching methods for nonexperimental causal studies.
*The Review of Economics and Statistics*, 84(1):151–161, 02 2002.
ISSN 0034-6535.
doi: 10.1162/003465302317331982.
URL [https://doi.org/10.1162/003465302317331982](https://doi.org/10.1162/003465302317331982).
- Dorie (2016)
Vincent Dorie.
vdorie/npci, 2016.
URL [https://github.com/vdorie/npci](https://github.com/vdorie/npci).
GitHub repository.
- Gao (2025)
Zijun Gao.
Trustworthy assessment of heterogeneous treatment effect estimator
via analysis of relative error.
In *The 28th International Conference on Artificial Intelligence
and Statistics*, 2025.
URL [https://openreview.net/forum?id=kOTUgBknsK](https://openreview.net/forum?id=kOTUgBknsK).
- Gutierrez & Gérardy (2017)
Pierre Gutierrez and Jean-Yves Gérardy.
Causal inference and uplift modelling: A review of the literature.
In Claire Hardgrove, Louis Dorard, Keiran Thompson, and Florian
Douetteau (eds.), *Proceedings of The 3rd International Conference on
Predictive Applications and APIs*, volume 67 of *Proceedings of Machine
Learning Research*, pp. 1–13. PMLR, 11–12 Oct 2017.
URL [https://proceedings.mlr.press/v67/gutierrez17a.html](https://proceedings.mlr.press/v67/gutierrez17a.html).
- Hassanpour & Greiner (2020)
Negar Hassanpour and Russell Greiner.
Learning disentangled representations for counterfactual regression.
In *International Conference on Learning Representations*, 2020.
URL [https://openreview.net/forum?id=HkxBJT4YvB](https://openreview.net/forum?id=HkxBJT4YvB).
- Hernán & Robins (2020)
M.A. Hernán and J. M. Robins.
*Causal Inference: What If*.
Boca Raton: Chapman and Hall/CRC, 2020.
- Imai & Ratkovic (2014)
Kosuke Imai and Marc Ratkovic.
Covariate balancing propensity score.
*Journal of the Royal Statistical Society (Series B)*,
76(1):243–263, 2014.
- Imbens & Rubin (2015)
G. W. Imbens and D. B. Rubin.
*Causal Inference For Statistics Social and Biomedical Science*.
Cambridge University Press, 2015.
- Jeong & Namkoong (2020)
Sookyo Jeong and Hongseok Namkoong.
Robust causal inference under covariate shift via worst-case
subpopulation treatment effects.
In Jacob Abernethy and Shivani Agarwal (eds.), *Proceedings of
Thirty Third Conference on Learning Theory*, volume 125 of *Proceedings
of Machine Learning Research*, pp. 2079–2084. PMLR, 09–12 Jul 2020.
URL [https://proceedings.mlr.press/v125/jeong20a.html](https://proceedings.mlr.press/v125/jeong20a.html).
- Jing Qin & Huang (2024)
Moming Li Jing Qin, Yukun Liu and Chiung-Yu Huang.
Distribution-free prediction intervals under covariate shift, with an
application to causal inference.
*Journal of the American Statistical Association*, 0(0):1–2, 2024.
doi: 10.1080/01621459.2024.2356886.
URL [https://doi.org/10.1080/01621459.2024.2356886](https://doi.org/10.1080/01621459.2024.2356886).
- Kunzel et al. (2019)
Soren R Kunzel, Jasjeet S Sekhon, Peter J Bickel, and Bin Y.
Metalearners for estimating heterogeneous treatment effects using
machine learning.
*PNAS*, 116:4156–4165, 2019.
- Künzel et al. (2019)
Sören R. Künzel, Jasjeet S. Sekhon, Peter J. Bickel, and Bin Yu.
Metalearners for estimating heterogeneous treatment effects using
machine learning.
*Proceedings of the National Academy of Sciences*, 116(10):4156–4165, 2019.
doi: 10.1073/pnas.1804597116.
URL [https://www.pnas.org/doi/abs/10.1073/pnas.1804597116](https://www.pnas.org/doi/abs/10.1073/pnas.1804597116).
- LaLonde (1986)
Robert J. LaLonde.
Evaluating the econometric evaluations of training programs with
experimental data.
*The American Economic Review*, 76(4):604–620, 1986.
ISSN 00028282.
URL [http://www.jstor.org/stable/1806062](http://www.jstor.org/stable/1806062).
- Louizos et al. (2017)
Christos Louizos, Uri Shalit, Joris M Mooij, David Sontag, Richard Zemel, and
Max Welling.
Causal effect inference with deep latent-variable models.
*Advances in neural information processing systems*, 30, 2017.
- Mahajan et al. (2024)
Divyat Mahajan, Ioannis Mitliagkas, Brady Neal, and Vasilis Syrgkanis.
Empirical analysis of model selection for heterogeneous causal effect
estimation.
*arXiv preprint arXiv:2211.01939*, 2024.
- Murphy (2022)
Kevin P. Murphy.
*Probabilistic Machine Learning: An introduction*.
MIT Press, 2022.
URL [http://probml.github.io/book1](http://probml.github.io/book1).
- Newey (1990)
Whitney K. Newey.
Semiparametric efficiency bounds.
*Journal of Applied Econometrics*, 5:99–135, 1990.
- Neyman (1990)
Jerzy Splawa Neyman.
On the application of probability theory to agricultural experiments.
essay on principles. section 9.
*Statistical Science*, 5:465–472, 1990.
- Pearl (2009)
Judea Pearl.
*Causality*.
Cambridge university press, 2009.
- Powers et al. (2017)
Scott Powers, Junyang Qian, Kenneth Jung, Alejandro Schuler, Nigam H. Shah,
Trevor Hastie, and Robert Tibshirani.
Some methods for heterogeneous treatment effect estimation in
high-dimensions, 2017.
URL [https://arxiv.org/abs/1707.00102](https://arxiv.org/abs/1707.00102).
- Rolling & Yang (2014)
Craig A. Rolling and Yuhong Yang.
Model selection for estimating treatment effects.
*Journal of the Royal Statistical Society Series B: Statistical
Methodology*, 76(4):749–769, 2014.
- Rosenbaum (2020)
Paul R. Rosenbaum.
*Design of Observational Studies*.
Springer Nature Switzerland AG, second edition, 2020.
- Rosenbaum & Rubin (1983)
Paul R Rosenbaum and Donald B Rubin.
The central role of the propensity score in observational studies for
causal effects.
*Biometrika*, 70(1):41–55, 1983.
- Rubin (1974)
D. B. Rubin.
Estimating causal effects of treatments in randomized and
nonrandomized studies.
*Journal of educational psychology*, 66:688–701,
1974.
- Saito & YasuiAuthors (2023)
Yuta Saito and Shota YasuiAuthors.
Counterfactual cross-validation: stable model selection procedure for
causal inference models.
In *Proceedings of the 37th International Conference on Machine
Learning*, ICML’20. JMLR, 2023.
- Semenova & Chernozhukov (2021)
Vira Semenova and Victor Chernozhukov.
Debiased machine learning of conditional average treatment effects
and and other causal functions.
*The Econometrics Journal*, 24:264–289, 2021.
- Shalit et al. (2017a)
Uri Shalit, Fredrik D. Johansson, and David Sontag.
Estimating individual treatment effect: generalization bounds and
algorithms.
In Doina Precup and Yee Whye Teh (eds.), *Proceedings of the
34th International Conference on Machine Learning*, volume 70 of
*Proceedings of Machine Learning Research*, pp. 3076–3085. PMLR,
06–11 Aug 2017a.
URL [https://proceedings.mlr.press/v70/shalit17a.html](https://proceedings.mlr.press/v70/shalit17a.html).
- Shalit et al. (2017b)
Uri Shalit, Fredrik D. Johansson, and David Sontag.
Estimating individual treatment effect: generalization bounds and
algorithms, 2017b.
URL [https://arxiv.org/abs/1606.03976](https://arxiv.org/abs/1606.03976).
- Shi et al. (2019a)
Claudia Shi, David Blei, and Victor Veitch.
Adapting neural networks for the estimation of treatment effects.
In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (eds.), *Advances in Neural
Information Processing Systems*, volume 32. Curran Associates, Inc.,
2019a.
URL
[https://proceedings.neurips.cc/paper_files/paper/2019/file/8fb5f8be2aa9d6c64a04e3ab9f63feee-Paper.pdf](https://proceedings.neurips.cc/paper_files/paper/2019/file/8fb5f8be2aa9d6c64a04e3ab9f63feee-Paper.pdf).
- Shi et al. (2019b)
Claudia Shi, David Blei, and Victor Veitch.
Adapting neural networks for the estimation of treatment effects.
*Advances in neural information processing systems*, 32,
2019b.
- van der Vaart (1998)
Aad W. van der Vaart.
*Asymptotic statistics*.
Cambridge University Press, 1998.
- Wager & Athey (2018a)
Stefan Wager and Susan Athey.
Estimation and inference of heterogeneous treatment effects using
random forests.
*Journal of the American Statistical Association*, 113(523):1228–1242, 2018a.
doi: 10.1080/01621459.2017.1319839.
URL [https://doi.org/10.1080/01621459.2017.1319839](https://doi.org/10.1080/01621459.2017.1319839).
- Wager & Athey (2018b)
Stefan Wager and Susan Athey.
Estimation and inference of heterogeneous treatment effects using
random forests.
*Journal of the American Statistical Association*, 113:1228–1242, 2018b.
- Wang et al. (2023)
Hao Wang, Zhichao Chen, Jiajun Fan, Haoxuan Li, Tianqiao Liu, Weiming Liu,
Quanyu Dai, Yichao Wang, Zhenhua Dong, and Ruiming Tang.
Optimal transport for treatment effect estimation, 2023.
URL [https://arxiv.org/abs/2310.18286](https://arxiv.org/abs/2310.18286).
- Wu et al. (2022)
Anpeng Wu, Junkun Yuan, Kun Kuang, Bo Li, Runze Wu, Qiang Zhu, Yueting Zhuang,
and Fei Wu.
Learning decomposed representations for treatment effect estimation.
*IEEE Transactions on Knowledge and Data Engineering*,
35(5):4989–5001, 2022.
- Yao et al. (2018)
Liuyi Yao, Sheng Li, Yaliang Li, Mengdi Huai, Jing Gao, and Aidong Zhang.
Representation learning for treatment effect estimation from
observational data.
In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi,
and R. Garnett (eds.), *Advances in Neural Information Processing
Systems*, volume 31. Curran Associates, Inc., 2018.
URL
[https://proceedings.neurips.cc/paper_files/paper/2018/file/a50abba8132a77191791390c3eb19fe7-Paper.pdf](https://proceedings.neurips.cc/paper_files/paper/2018/file/a50abba8132a77191791390c3eb19fe7-Paper.pdf).
- Yoon et al. (2018)
Jinsung Yoon, James Jordon, and Mihaela Van Der Schaar.
Ganite: Estimation of individualized treatment effects using
generative adversarial nets.
In *International Conference on Learning Representations*, 2018.
- Zhong et al. (2022)
Kailiang Zhong, Fengtong Xiao, Yan Ren, Yaorong Liang, Wenqing Yao, Xiaofeng
Yang, and Ling Cen.
Descn: Deep entire space cross networks for individual treatment
effect estimation.
In *Proceedings of the 28th ACM SIGKDD Conference on Knowledge
Discovery and Data Mining*, KDD ’22, pp. 4612–4620. ACM, August 2022.
doi: 10.1145/3534678.3539198.
URL [http://dx.doi.org/10.1145/3534678.3539198](http://dx.doi.org/10.1145/3534678.3539198).

## Appendix A Notation Summary

**Table 4: Notation and their meanings.**
| Symbol | Meaning |
| --- | --- |
| $A$ | Binary treatment variable |
| $X$ | Pre-treatment covariates |
| $Y$ | Outcome |
| $\tau(x)$ | Individual treatment effect |
| $e(x)$ | Propensity score |
| $\mu_{a}(x)$ | Outcome regression function, i.e., $\mu_{a}(x)=\mathbb{E}[Y\mid X=x,A=a]$ for $a=0,1$ |
| $\delta(\hat{\tau}_{1},\hat{\tau}_{2})$ | Relative error between estimator $\hat{\tau}_{1}$ and $\hat{\tau}_{2}$ |
| $\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2})$ | Estimated relative error between estimator $\hat{\tau}_{1}$ and $\hat{\tau}_{2}$ |
| $\check{e},\check{\mu}_{a},\check{\gamma},\check{\beta}_{a}$ | Nuisance estimators for propensity score, conditional outcomes and their coefficients |
| $\bar{e},\bar{\mu}_{a},\bar{\gamma},\bar{\beta}_{a}$ | Probability limits of propensity score, conditional outcomes and their coefficients |
| $\Phi(X)$ | The shared representation of $X$, defined in Eq. (1) & (2) |

## Appendix B Merits of Relative Error

There are several advantages of using relative error over absolute error.

- •
(1) Weaker condition. Condition [2](https://arxiv.org/html/2510.16419v1#Thmcondition2) is strictly weaker than Condition [1](https://arxiv.org/html/2510.16419v1#Thmcondition1). Condition [1](https://arxiv.org/html/2510.16419v1#Thmcondition1) requires that *all nuisance parameter estimators* converge to their true values at a rate faster than $n^{-1/4}$. In contrast,
Condition [2](https://arxiv.org/html/2510.16419v1#Thmcondition2) imposes a weaker requirement—only that the product of the bias, $(\tilde{\mu}_{a}(x)-\mu_{a}(x))(\tilde{e}(x)-e(x))$, converges at a rate of order $o_{\mathbb{P}}(n^{-1/2})$ as well as the nuisance function estimators being consistent. This allows for cases where $\tilde{e}(x)$ converges at a rate of $o_{\mathbb{P}}(n^{-1/5})$ and $\tilde{\mu}_{a}(x)$ converges at a rate of $o_{\mathbb{P}}(n^{-1/3})$.
- •
(2) Easier to compare multiple estimators. When comparing two estimators $\hat{\tau}_{1}(x)$ and $\hat{\tau}_{2}(x)$ in terms of absolute error,
although both $\hat{\phi}(\hat{\tau}_{1})$ and $\hat{\phi}(\hat{\tau}_{2})$ are asymptotically normal, we cannot directly construct a confidence interval for $\hat{\phi}(\hat{\tau}_{1})-\hat{\phi}(\hat{\tau}_{2})$
due to their dependency (as they use the same test data and share the same nuisance parameter estimates). In contrast, $\hat{\delta}(\hat{\tau}_{1},\hat{\tau}_{2})$ does not suffer such a problem.
- •
(3) Double robustness. When we replace $o_{\mathbb{P}}(n^{-1/2})$ in Conditions [1](https://arxiv.org/html/2510.16419v1#Thmcondition1) and [2](https://arxiv.org/html/2510.16419v1#Thmcondition2) with $o_{\mathbb{P}}(1)$, both $\hat{\phi}(\hat{\tau}_{1})$ and $\hat{\delta}(\hat{\tau}_{1},\hat{\tau}_{2})$ are consistent (asymptotically unbiased) under their respective conditions. Thus, from Condition [2](https://arxiv.org/html/2510.16419v1#Thmcondition2), $\hat{\delta}(\hat{\tau}_{1},\hat{\tau}_{2})$ exhibits the property of double robustness, meaning it is a consistent estimator if either $\tilde{e}(x)$ is consistent or $\tilde{\mu}_{a}(x)$ for $a=0,1$ are consistent. However, $\hat{\phi}(\hat{\tau}_{1})$ dose not possess this property.

## Appendix C Illustration of Neural Network Structure

Figure [3](https://arxiv.org/html/2510.16419v1#A3.F3) shows the schematic structure of our proposed network.
The input covariates $X\in\mathbb{R}^{p}$ are passed through fully connected hidden layers to obtain a shared representation $\Phi(X)\in\mathbb{R}^{d}$.
This representation is fed into three heads: the control outcome head $\mu_{0}(X)$, the treated outcome head $\mu_{1}(X)$, and the treatment head $e(X)$.
The outcome heads contribute to the weighted least square loss $\mathcal{L}_{\text{wls}}$, the treatment head contributes to the cross entropy loss $\mathcal{L}_{\text{ce}}$, and the shared representation is regularized by the constraint loss $\mathcal{L}_{\text{const}}$.
The total objective is given by

$$ $\mathcal{L}=\mathcal{L}_{\mathrm{wls}}+\lambda_{1}\mathcal{L}_{\mathrm{ce}}+\lambda_{2}\mathcal{L}_{\mathrm{const}}.$ $$

Figure: Figure 3: Neural Network Structure
Refer to caption: network.png

## Appendix D Proof of Theorem 1

Theorem 1. If the propensity score model is correctly specified, and $\check{\gamma}$, $\check{\beta}_{0}$ as well as $\check{\beta}_{1}$ converge to their probability limits at a rate faster than $n^{-1/4}$, then we have

$$ $\displaystyle\sqrt{n}\{\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})-\delta(\hat{\tau}_{1},\hat{\tau}_{2})\}\xrightarrow{d}\mathcal{N}(0,\sigma^{2}),$ $$

where $\sigma^{2}=\text{Var}\{\varphi(Z;\bar{u}_{0},\bar{u}_{1},\bar{e})\}$ and $\xrightarrow{d}$ means convergence in distribution.

*Proof of Theorem 1*. As discussed in Section [4.1](https://arxiv.org/html/2510.16419v1#S4.SS1), we first show that

$$ $\displaystyle\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})-\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\bar{\gamma},\bar{\beta}_{0},\bar{\beta}_{1})=o_{\mathbb{P}}(n^{-1/2}).$ (A.1) $$

By a Taylor expansion of $\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})$ around $(\bar{\gamma},\bar{\beta}_{0},\bar{\beta}_{1})$, we obtain

$$ $\displaystyle\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})-\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\bar{\gamma},\bar{\beta}_{0},\bar{\beta}_{1})$ $\displaystyle=$ $\displaystyle\Delta_{\gamma}^{\intercal}(\check{\gamma}-\bar{\gamma})+\Delta_{\beta_{0}}^{\intercal}(\check{\beta}_{0}-\bar{\beta}_{0})+\Delta_{\beta_{1}}^{\intercal}(\check{\beta}_{1}-\bar{\beta}_{1})+O_{\mathbb{P}}((\check{\gamma}-\bar{\gamma})^{2}+(\check{\gamma}-\bar{\gamma})(\check{\beta}_{1}-\bar{\beta}_{1})+(\check{\gamma}-\bar{\gamma})(\check{\beta}_{0}-\bar{\beta}_{0})),$ $$

where

$$ $\displaystyle\Delta_{\gamma}={}$ $\displaystyle-\frac{1}{n}\sum_{i=1}^{n}2(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\Big(\frac{A_{i}(1-\bar{e}(X_{i}))(Y_{i}-\bar{\mu}_{1}(X_{i}))}{\bar{e}(X_{i})}+\frac{(1-A_{i})\bar{e}(X_{i})(Y_{i}-\bar{\mu}_{0}(X_{i}))}{1-\bar{e}(X_{i})}\Big)\Phi_{1}(X_{i}),$ $\displaystyle\Delta_{\beta_{0}}={}$ $\displaystyle-\frac{1}{n}\sum_{i=1}^{n}2(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\left(1-\frac{1-A_{i}}{1-\bar{e}(X_{i})}\right)\Phi_{2}(X_{i}),$ $\displaystyle\Delta_{\beta_{1}}={}$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}2(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\left(1-\frac{A_{i}}{\bar{e}(X_{i})}\right)\Phi_{2}(X_{i}).$ $$

By the condition that $\check{\gamma}=\bar{\gamma}+o_{\mathrm{\mathbb{P}}}(n^{-1/4})$, $\check{\beta}_{0}=\bar{\beta}_{0}+o_{\mathrm{\mathbb{P}}}(n^{-1/4})$ and $\check{\beta}_{0}=\bar{\beta}_{0}+o_{\mathrm{\mathbb{P}}}(n^{-1/4})$, we obtain
$O_{\mathbb{P}}((\check{\gamma}-\bar{\gamma})^{2}+(\check{\gamma}-\bar{\gamma})(\check{\beta}_{1}-\bar{\beta}_{1})+(\check{\gamma}-\bar{\gamma})(\check{\beta}_{0}-\bar{\beta}_{0}))=o_{\mathbb{P}}(n^{-1/2})$. Thus, we only need to deal with $\Delta_{\gamma}^{\intercal}(\check{\gamma}-\bar{\gamma})$, $\Delta_{\beta_{0}}^{\intercal}(\check{\beta}_{0}-\bar{\beta}_{0})$ and $\Delta_{\beta_{1}}^{\intercal}(\check{\beta}_{1}-\bar{\beta}_{1})$.

Since $\bar{\gamma}$, $\bar{\beta}_{0}$ and $\bar{\beta}_{1}$ are the probability limits of $\check{\gamma}$, $\check{\beta}_{0}$ and $\check{\beta}_{1}$, respectively, we obtain $\check{\gamma}-\bar{\gamma}=o_{\mathbb{P}}(1)$, $\check{\beta}_{0}-\bar{\beta}_{0}=o_{\mathbb{P}}(1)$ and $\check{\beta}_{1}-\bar{\beta}_{1}=o_{\mathbb{P}}(1)$.

Then, it sufficies to show that $\Delta_{\gamma}=O_{\mathbb{P}}(n^{-1/2})$, $\Delta_{\beta_{0}}=O_{\mathbb{P}}(n^{-1/2})$ and $\Delta_{\beta_{1}}=O_{\mathbb{P}}(n^{-1/2})$. By CLT, $\Delta_{\gamma}-\mathbb{E}(\Delta_{\gamma})=O_{\mathbb{P}}(n^{-1/2})$, $\Delta_{\beta_{0}}-\mathbb{E}(\Delta_{\beta_{0}})=O_{\mathbb{P}}(n^{-1/2})$ and $\Delta_{\beta_{1}}-\mathbb{E}(\Delta_{\beta_{1}})=O_{\mathbb{P}}(n^{-1/2})$, then, we only need to show that $\mathbb{E}(\Delta_{\gamma})=\mathbb{E}(\Delta_{\beta_{0}})=\mathbb{E}(\Delta_{\beta_{1}})=0$.

We first deal with $\mathbb{E}(\Delta_{\gamma})$.

$$ $\displaystyle\mathbb{E}(\Delta_{\gamma})=$ $\displaystyle\mathbb{E}\left(-\frac{1}{n}\sum_{i=1}^{n}2(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\Big(\frac{A_{i}(1-\bar{e}(X_{i}))(Y_{i}-\bar{\mu}_{1}(X_{i}))}{\bar{e}(X_{i})}+\frac{(1-A_{i})\bar{e}(X_{i})(Y_{i}-\bar{\mu}_{0}(X_{i}))}{1-\bar{e}(X_{i})}\Big)\Phi_{1}(X_{i})\right)$ $\displaystyle=$ $\displaystyle 2\mathbb{E}\left((\hat{\tau}_{1}(X)-\hat{\tau}_{2}(X))\Big(\frac{A(1-\bar{e}(X))(Y-\bar{\mu}_{1}(X))}{\bar{e}(X)}+\frac{(1-A)\bar{e}(X)(Y-\bar{\mu}_{0}(X))}{1-\bar{e}(X)}\Big)\Phi_{1}(X)\right)$ $\displaystyle=$ $\displaystyle 2\mathbb{E}\left((\hat{\tau}_{1}(X)-\hat{\tau}_{2}(X))\frac{A(1-\bar{e}(X))(Y-\bar{\mu}_{1}(X))}{\bar{e}(X)}\Phi_{1}(X)\right)$ $\displaystyle+2\mathbb{E}\left((\hat{\tau}_{1}(X)-\hat{\tau}_{2}(X))\frac{(1-A)\bar{e}(X)(Y-\bar{\mu}_{0}(X))}{1-\bar{e}(X)}\Phi_{1}(X)\right)$ $\displaystyle=$ $\displaystyle 0.$ $$

The last equation holds by the definition of $\bar{\beta}_{0}$ and $\bar{\beta}_{1}$ and the fact that $\Phi_{1}(X)$ is a sub-vector of $\Phi_{2}(X)$.

We then deal with $\mathbb{E}(\Delta_{\beta_{0}})$.

$$ $\displaystyle\mathbb{E}(\Delta_{\beta_{0}})={}$ $\displaystyle\mathbb{E}\left(-\frac{1}{n}\sum_{i=1}^{n}2(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\left(1-\frac{1-A_{i}}{1-\bar{e}(X_{i})}\right)\Phi_{2}(X_{i})\right)$ $\displaystyle=$ $\displaystyle-2\mathbb{E}\left((\hat{\tau}_{1}(X)-\hat{\tau}_{2}(X))\left(1-\frac{1-A}{1-\bar{e}(X)}\right)\Phi_{2}(X)\right)$ $\displaystyle=$ $\displaystyle-2\mathbb{E}_{X}\left(\left.\mathbb{E}\left((\hat{\tau}_{1}(X)-\hat{\tau}_{2}(X))\left(1-\frac{1-A}{1-\bar{e}(X)}\right)\Phi_{2}(X)\right|X\right)\right)$ $\displaystyle=$ $\displaystyle 0.$ $$

The last equation holds since the PS model is correct. Finally, we handle $\mathbb{E}(\Delta_{\beta_{1}})$.

$$ $\displaystyle\mathbb{E}(\Delta_{\beta_{1}})={}$ $\displaystyle\mathbb{E}\left(\frac{1}{n}\sum_{i=1}^{n}2(\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i}))\left(1-\frac{A_{i}}{\bar{e}(X_{i})}\right)\Phi_{2}(X_{i})\right)$ $\displaystyle=$ $\displaystyle 2\mathbb{E}\left((\hat{\tau}_{1}(X)-\hat{\tau}_{2}(X))\left(1-\frac{A}{\bar{e}(X)}\right)\Phi_{2}(X)\right)$ $\displaystyle=$ $\displaystyle 2\mathbb{E}_{X}\left(\mathbb{E}\left(\left.(\hat{\tau}_{1}(X)-\hat{\tau}_{2}(X))\left(1-\frac{A}{\bar{e}(X)}\right)\Phi_{2}(X)\right|X\right)\right)$ $\displaystyle=$ $\displaystyle 0.$ $$

The last equation holds since the PS model is correct. Therefore, equation ([A.1](https://arxiv.org/html/2510.16419v1#A4.E1)) holds.

We then want to show that

$$ $\displaystyle\sqrt{n}(\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\bar{\gamma},\bar{\beta}_{0},\bar{\beta}_{1})-\delta(\hat{\tau}_{1},\hat{\tau}_{2}))\rightarrow\mathcal{N}(0,\sigma^{2}).$ (A.2) $$

By definiton,

$$ $\displaystyle\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\bar{\gamma},\bar{\beta}_{0},\bar{\beta}_{1})$ $\displaystyle=$ $\displaystyle{}\frac{1}{n}\sum_{i=1}^{n}\{\hat{\tau}_{1}^{2}(X_{i})-\hat{\tau}_{2}^{2}(X_{i})\}$ $\displaystyle-{}$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left(2\{\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i})\}\cdot\left[\frac{A_{i}\{Y_{i}-\bar{\mu}_{1}(X_{i})\}}{\bar{e}(X_{i})}+\bar{\mu}_{1}(X_{i})-\frac{(1-A_{i})\{Y_{i}-\bar{\mu}_{0}(X_{i})\}}{1-\bar{e}(X_{i})}-\bar{\mu}_{0}(X_{i})\right]\right)$ $$

If model ([1](https://arxiv.org/html/2510.16419v1#S4.E1)) is correct,

$$ $\displaystyle\mathbb{E}\left[\{\hat{\tau}_{1}^{2}(X_{i})-\hat{\tau}_{2}^{2}(X_{i})\}\right]$ $\displaystyle-{}$ $\displaystyle\mathbb{E}\left\{\left(2\{\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i})\}\cdot\left[\frac{A_{i}\{Y_{i}-\bar{\mu}_{1}(X_{i})\}}{\bar{e}(X_{i})}+\bar{\mu}_{1}(X_{i})-\frac{(1-A_{i})\{Y_{i}-\bar{\mu}_{0}(X_{i})\}}{1-\bar{e}(X_{i})}-\bar{\mu}_{0}(X_{i})\right]\right)\right\}$ $\displaystyle=$ $\displaystyle{}\mathbb{E}\left[\{\hat{\tau}_{1}^{2}(X)-\hat{\tau}_{2}^{2}(X)\}\right]$ $\displaystyle-{}$ $\displaystyle\left.\mathbb{E}_{X_{i}}\mathbb{E}\left[\left\{\left(2\{\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i})\}\cdot\left[\frac{A_{i}\{Y_{i}-\bar{\mu}_{1}(X_{i})\}}{\bar{e}(X_{i})}+\bar{\mu}_{1}(X_{i})-\frac{(1-A_{i})\{Y_{i}-\bar{\mu}_{0}(X_{i})\}}{1-\bar{e}(X_{i})}-\bar{\mu}_{0}(X_{i})\right]\right)\right\}\right|X_{i}\right]$ $\displaystyle=$ $\displaystyle{}\mathbb{E}\left[\{\hat{\tau}_{1}^{2}(X)-\hat{\tau}_{2}^{2}(X)\}\right]$ $\displaystyle-{}$ $\displaystyle\mathbb{E}_{X_{i}}\mathbb{E}\left[\left\{\left(2\{\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i})\}\cdot\left[\frac{A_{i}\{Y_{i}(1)-\bar{\mu}_{1}(X_{i})\}}{\bar{e}(X_{i})}+\bar{\mu}_{1}(X_{i})\right.\right.\right.\right.$ $\displaystyle\left.\left.\left.\left.\left.-\frac{(1-A_{i})\{Y_{i}(0)-\bar{\mu}_{0}(X_{i})\}}{1-\bar{e}(X_{i})}-\bar{\mu}_{0}(X_{i})\right]\right)\right\}\right|X_{i}\right]$ $\displaystyle=$ $\displaystyle{}\mathbb{E}\left[\{\hat{\tau}_{1}^{2}(X)-\hat{\tau}_{2}^{2}(X)\}\right]$ $\displaystyle-{}$ $\displaystyle\mathbb{E}_{X_{i}}\mathbb{E}\left[\left\{\left(2\left\{\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i})\}\cdot\left[\{Y_{i}(1)-\bar{\mu}_{1}(X_{i})\}+\bar{\mu}_{1}(X_{i})-\{Y_{i}(0)-\bar{\mu}_{0}(X_{i})\}-\bar{\mu}_{0}(X_{i})\right]\right)\right\}\right|X_{i}\right]$ $\displaystyle=$ $\displaystyle{}\mathbb{E}\left[\{\hat{\tau}_{1}^{2}(X)-\hat{\tau}_{2}^{2}(X)\}\right]-\mathbb{E}_{X_{i}}\mathbb{E}[\{(2\{\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i})\}\tau(X_{i})|X_{i}]$ $\displaystyle=$ $\displaystyle{}\mathbb{E}\left[\{\hat{\tau}_{1}^{2}(X)-\hat{\tau}_{2}^{2}(X)\}-2\{\hat{\tau}_{1}(X)-\hat{\tau}_{2}(X)\}\tau(X)\right]$ $\displaystyle=$ $\displaystyle\delta(\hat{\tau}_{1},\hat{\tau}_{2}).$ $$

If model ([2](https://arxiv.org/html/2510.16419v1#S4.E2)) is correct,

$$ $\displaystyle\mathbb{E}\left[\{\hat{\tau}_{1}^{2}(X_{i})-\hat{\tau}_{2}^{2}(X_{i})\}\right]$ $\displaystyle-{}$ $\displaystyle\mathbb{E}\left\{\left(2\{\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i})\}\cdot\left[\frac{A_{i}\{Y_{i}-\bar{\mu}_{1}(X_{i})\}}{\bar{e}(X_{i})}+\bar{\mu}_{1}(X_{i})-\frac{(1-A_{i})\{Y_{i}-\bar{\mu}_{0}(X_{i})\}}{1-\bar{e}(X_{i})}-\bar{\mu}_{0}(X_{i})\right]\right)\right\}$ $\displaystyle=$ $\displaystyle{}\mathbb{E}\left[\{\hat{\tau}_{1}^{2}(X)-\hat{\tau}_{2}^{2}(X)\}\right]$ $\displaystyle-{}$ $\displaystyle\left.\mathbb{E}_{X_{i}}\mathbb{E}\left[\left\{\left(2\{\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i})\}\cdot\left[\frac{A_{i}\{Y_{i}-\bar{\mu}_{1}(X_{i})\}}{\bar{e}(X_{i})}+\bar{\mu}_{1}(X_{i})-\frac{(1-A_{i})\{Y_{i}-\bar{\mu}_{0}(X_{i})\}}{1-\bar{e}(X_{i})}-\bar{\mu}_{0}(X_{i})\right]\right)\right\}\right|X_{i}\right]$ $\displaystyle=$ $\displaystyle{}\mathbb{E}\left[\{\hat{\tau}_{1}^{2}(X)-\hat{\tau}_{2}^{2}(X)\}\right]$ $\displaystyle-{}$ $\displaystyle\mathbb{E}_{X_{i}}\mathbb{E}\left[\left\{\left(2\{\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i})\}\cdot\left[\frac{A_{i}\{Y_{i}(1)-\bar{\mu}_{1}(X_{i})\}}{\bar{e}(X_{i})}+\bar{\mu}_{1}(X_{i})\right.\right.\right.\right.$ $\displaystyle\left.\left.\left.\left.\left.-\frac{(1-A_{i})\{Y_{i}(0)-\bar{\mu}_{0}(X_{i})\}}{1-\bar{e}(X_{i})}-\bar{\mu}_{0}(X_{i})\right]\right)\right\}\right|X_{i}\right]$ $\displaystyle=$ $\displaystyle{}\mathbb{E}\left[\{\hat{\tau}_{1}^{2}(X)-\hat{\tau}_{2}^{2}(X)\}\right]$ $\displaystyle-{}$ $\displaystyle\mathbb{E}_{X_{i}}\mathbb{E}\left[\left\{\left(2\{\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i})\}\cdot\left[\frac{A_{i}\{\bar{\mu}_{1}(X_{i})-\bar{\mu}_{1}(X_{i})\}}{\bar{e}(X_{i})}+\bar{\mu}_{1}(X_{i})\right.\right.\right.\right.$ $\displaystyle\left.\left.\left.\left.\left.-\frac{(1-A_{i})\{\bar{\mu}_{0}(X_{i})-\bar{\mu}_{0}(X_{i})\}}{1-\bar{e}(X_{i})}-\bar{\mu}_{0}(X_{i})\right]\right)\right\}\right|X_{i}\right]$ $\displaystyle=$ $\displaystyle{}\mathbb{E}\left[\{\hat{\tau}_{1}^{2}(X)-\hat{\tau}_{2}^{2}(X)\}\right]-\mathbb{E}_{X_{i}}\mathbb{E}[\{(2\{\hat{\tau}_{1}(X_{i})-\hat{\tau}_{2}(X_{i})\}\tau(X_{i})|X_{i}]$ $\displaystyle=$ $\displaystyle{}\mathbb{E}\left[\{\hat{\tau}_{1}^{2}(X)-\hat{\tau}_{2}^{2}(X)\}-2\{\hat{\tau}_{1}(X)-\hat{\tau}_{2}(X)\}\tau(X)\right]$ $\displaystyle=$ $\displaystyle\delta(\hat{\tau}_{1},\hat{\tau}_{2}).$ $$

Therefore, when at least one of models ([1](https://arxiv.org/html/2510.16419v1#S4.E1)) or  ([2](https://arxiv.org/html/2510.16419v1#S4.E2)) is correct, $\bar{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\bar{\gamma},\bar{\beta}_{0},\bar{\beta}_{1})-\delta(\hat{\tau}_{1},\hat{\tau}_{2})$ is the average of i.i.d. observations with mean 0 and variance $\text{Var}\{\varphi(Z;\bar{u}_{0},\bar{u}_{1},\bar{e})\}$. By CLT, [A.2](https://arxiv.org/html/2510.16419v1#A4.E2) holds with $\sigma^{2}=\text{Var}\{\varphi(Z;\bar{u}_{0},\bar{u}_{1},\bar{e})\}$.

$\Box$

### D.1 Proof of Proposition 2

Proposition 2.
Under the conditions in Theorem [1](https://arxiv.org/html/2510.16419v1#Thmtheorem1),
a consistent estimator of $\sigma^{2}$ is

$$ $\displaystyle\hat{\sigma}^{2}=\frac{1}{n}\sum^{n}_{i=1}\left\{\varphi(Z_{i};\check{u}_{0},\check{u}_{1},\check{e})-\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})\right\}^{2},$ $$

an asymptotic $(1-\eta)$ confidence interval for $\delta(\hat{\tau}_{1},\hat{\tau}_{2})$ is $\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})\pm z_{\eta/2}\sqrt{\hat{\sigma}^{2}/n}$, where $z_{\eta/2}$ is the $(1-\eta/2)$ quantile of the standard normal distribution.

###### Proof of Proposition 2.

$$ $\displaystyle\hat{\sigma}^{2}-\sigma^{2}=$ $\displaystyle\frac{1}{n}\sum^{n}_{i=1}\left\{\varphi(Z_{i};\check{u}_{0},\check{u}_{1},\check{e})-\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})\right\}^{2}-\mathbb{E}\left\{\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})-\mathbb{E}\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})\right\}^{2}$ $\displaystyle=$ $\displaystyle\underbrace{\frac{1}{n}\sum^{n}_{i=1}\left\{\varphi(Z_{i};\check{u}_{0},\check{u}_{1},\check{e})-\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})\right\}^{2}-\frac{1}{n}\sum^{n}_{i=1}\left\{\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})-\mathbb{E}\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})\right\}^{2}}_{①}$ $\displaystyle\underbrace{+\frac{1}{n}\sum^{n}_{i=1}\left\{\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})-\mathbb{E}\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})\right\}^{2}-\mathbb{E}\left\{\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})-\mathbb{E}\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})\right\}^{2}}_{②}.$ $$

By the law of large numbers (LLN), $②\overset{p}{\rightarrow}0$, we only need to deal with $①$:

$$ $\displaystyle\frac{1}{n}\sum^{n}_{i=1}\left\{\varphi(Z_{i};\check{u}_{0},\check{u}_{1},\check{e})-\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})\right\}^{2}$ $\displaystyle=$ $\displaystyle\frac{1}{n}\sum^{n}_{i=1}\left\{\varphi(Z_{i};\check{u}_{0},\check{u}_{1},\check{e})-\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})+\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})-\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})\right\}^{2}$ $\displaystyle=$ $\displaystyle\underbrace{\frac{1}{n}\sum^{n}_{i=1}\left\{\varphi(Z_{i};\check{u}_{0},\check{u}_{1},\check{e})-\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})\right\}^{2}}_{ⓐ}+\underbrace{\frac{1}{n}\sum^{n}_{i=1}\left\{\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})-\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})\right\}^{2}}_{ⓑ}$ $\displaystyle+\underbrace{\frac{2}{n}\sum^{n}_{i=1}\left\{\varphi(Z_{i};\check{u}_{0},\check{u}_{1},\check{e})-\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})\right\}\left\{\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})-\check{\delta}(\hat{\tau}_{1},\hat{\tau}_{2};\check{\gamma},\check{\beta}_{0},\check{\beta}_{1})\right\}}_{ⓒ}.$ $$

By LLN, $ⓐ\overset{p}{\rightarrow}0$, and $ⓒ\overset{p}{\rightarrow}0$. Similarly, we can obtain

$$ $\displaystyle ⓑ=\frac{1}{n}\sum^{n}_{i=1}\left\{\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})-\mathbb{E}\varphi(Z_{i};\bar{u}_{0},\bar{u}_{1},\bar{e})\right\}^{2}+o_{p}(1).$ $$

Therefore $①\overset{p}{\rightarrow}0$, which leads to $\hat{\sigma}^{2}\overset{p}{\rightarrow}\sigma^{2}$.

The asymptotic $1-\eta$ confidence interval is constructed by the standard theory.

∎

###### Proof of Proposition 2.

## Appendix E Experimental Details

### E.1 Dataset Details

IHDP. The IHDP dataset is based on a randomized controlled trial conducted as part of the Infant Health and Development Program. The goal is to assess the impact of specialist home visits on children’s future cognitive outcomes. Following Hill (and, 2011), a subset of treated units is removed to introduce selection bias, creating a semi-synthetic evaluation setting. The dataset contains 747 samples (139 treated and 608 control), each with 25 pre-treatment covariates. The simulated outcome is the same as that in Shalit et al.(2017) (Shalit et al., 2017a), by setting “A” in the NPCI package (Dorie, 2016).

Twins. The Twins dataset is constructed from twin births in the U.S.. For each twin pair, the heavier twin is assigned as the treated unit ($t_{i}=1$), and the lighter twin as the control ($t_{i}=0$). We extract 28 covariates related to parental, pregnancy, and birth characteristics from the original data and generate an additional 10 covariates following (Wu et al., 2022). The outcome of interest is the one-year mortality of each child. We restrict the analysis to same-sex twins with birth weights below 2000g and without any missing features, yielding a final dataset with 5,271 samples.
The treatment assignment mechanism is defined as:
$t_{i}\mid x_{i}\sim\mathrm{Bern}\left(\sigma(w^{\top}X+n)\right),$
where $\sigma(\cdot)$ is the sigmoid function, $w^{\top}\sim\mathcal{U}((-0.1,0.1)^{38\times 1})$, and $n\sim\mathcal{N}(0,0.1)$.

Jobs. The Jobs dataset is a standard benchmark in causal inference, originally introduced by LaLonde (1986) (LaLonde, 1986). It evaluates the impact of job training on employment outcomes by combining data from a randomized study (National Supported Work program) with observational records (PSID), following the setup of Smith and Todd (2005) (A. Smith & E. Todd, 2005). The dataset includes 297 treated units, 425 control units from the experimental sample, and 2490 control units from the observational sample. Each record consists of 8 covariates, such as age, education, ethnicity, and pre-treatment earnings. The task is framed as a binary classification problem predicting unemployment status post-treatment, using features defined by Dehejia and Wahba (2002) (Dehejia & Wahba, 2002).

### E.2 Choosing Different Nuisance Estimators

Experimental Set-up. As for the nuisance estimators, we choose linear regression and gradient boosting as Gao (Gao, 2025) used in their paper. The evaluation metrics are the same as those in Section [6](https://arxiv.org/html/2510.16419v1#S6). We provide the results on the IHDP and the Twins.

Experimental Results. Table [5](https://arxiv.org/html/2510.16419v1#A5.T5) summarizes the results on the IHDP and Twins datasets. When plugging conventional nuisance estimators (linear regression and gradient boosting) into the relative error framework, the resulting procedures do achieve nominal coverage. Nevertheless, the corresponding variance is so large that the confidence intervals frequently include zero, making it essentially impossible to tell which candidate estimator is superior. These baselines therefore serve as valid but uninformative references. In contrast, our proposed method not only maintains well-calibrated coverage but also delivers much higher selection accuracy, producing confidence intervals that are substantially tighter and practically useful for identifying the winner.

**Table 5: Relative Error Estimation Performance with Different Nuisance Estimators on the IHDP and Twins datasets.**
|  | IHDP | Twins |  |  |
| --- | --- | --- | --- | --- |
| Nuisance Estimators | Coverage Rate | Selection Accuracy | Coverage Rate | Selection Accuracy |
| Linear Regression | 0.94 | 0.44 | 0.94 | 0.88 |
| Gradient Boosting | 0.95 | 0.48 | 0.94 | 0.86 |
| Ours | 0.96 | 0.80 | 0.94 | 0.94 |

### E.3 Results on Jobs

Evaluation Metrics.
For the Jobs datasets, as there are no counterfactual outcomes, we report the true Average Treatment Effect on the Treated (ATT) and the Policy Risk ($\mathcal{R}_{\text{pol}}$) recommended by Shalit et al.(Shalit et al., 2017b).
Specifically, the policy risk can be estimated using only the randomized subset of the Jobs dataset:

$$ $\hat{\mathcal{R}}{\text{pol}}=1-\left(\frac{1}{|A_{1}\cap T_{1}\cap E|}\sum_{x_{i}\in A_{1}\cap T_{1}\cap E}y_{1}^{(i)}\cdot\frac{|A_{1}\cap E|}{|E|}+\frac{1}{|A_{0}\cap T_{0}\cap E|}\sum_{x_{i}\in A_{0}\cap T_{0}\cap E}y_{0}^{(i)}\cdot\frac{|A_{0}\cap E|}{|E|}\right)$ $$

where E denotes units from the experimental group, $A_{1}=\{x_{i}:\hat{y}_{1}^{(i)}-\hat{y}_{0}^{(i)}>0\},A_{0}=\{x_{i}:\hat{y}_{1}^{(i)}-\hat{y}_{0}^{(i)}<0\}$, and $T_{1},T_{0}$ are the treated and control subsets, respectively. Since all treated units $T$ belong to the randomized subset $E$, the true Average Treatment Effect on the Treated (ATT) can be identified and computed as:

$$ $\text{ATT}=\frac{1}{|T|}\sum_{i\in T}y_{i}-\frac{1}{|C\cap E|}\sum_{i\in C\cap E}y_{i}$ $$

where C denotes the control group. We evaluate estimation accuracy using the ATT error:
$\epsilon_{\text{ATT}}=\left|\text{ATT}-\frac{1}{|T|}\sum_{i\in T}\left(f(x_{i},1)-f(x_{i},0)\right)\right|.$

Accuracy of the CATE Estimation. We evaluate the performance of CATE estimation by our network and compare it with baselines mentioned in Section [6](https://arxiv.org/html/2510.16419v1#S6). We average over 20 realizations of our network, and the results are presented in Table [6](https://arxiv.org/html/2510.16419v1#A5.T6). One can clearly see that our proposed method achieves the best performance across all metrics, having the lowest $\hat{\mathcal{R}}{\text{pol}}$ and $\epsilon_{\text{ATT}}$ in both training sets and test sets.

**Table 6: Performance on the Jobs dataset (in-sample and out-of-sample).**
| Method | $\mathcal{R}_{\text{pol}}^{\text{in}}$ | $\epsilon_{\text{ATT}}^{\text{in}}$ | $\mathcal{R}_{\text{pol}}^{\text{out}}$ | $\epsilon_{\text{ATT}}^{\text{out}}$ |
| --- | --- | --- | --- | --- |
| LinDML | 0.158 $\pm$ 0.015 | 0.019 $\pm$ 0.015 | 0.183 $\pm$ 0.040 | 0.053 $\pm$ 0.051 |
| SpaDML | 0.150 $\pm$ 0.024 | 0.131 $\pm$ 0.118 | 0.165 $\pm$ 0.046 | 0.144 $\pm$ 0.134 |
| CForest | 0.114 $\pm$ 0.016 | 0.025 $\pm$ 0.018 | 0.155 $\pm$ 0.028 | 0.058 $\pm$ 0.047 |
| X-Learner | 0.169 $\pm$ 0.037 | 0.026 $\pm$ 0.015 | 0.173 $\pm$ 0.034 | 0.053 $\pm$ 0.050 |
| S-Learner | 0.148 $\pm$ 0.026 | 0.095 $\pm$ 0.040 | 0.160 $\pm$ 0.027 | 0.115 $\pm$ 0.070 |
| TarNet | 0.141 $\pm$ 0.005 | 0.183 $\pm$ 0.047 | 0.145 $\pm$ 0.009 | 0.190 $\pm$ 0.074 |
| Dragonnet | 0.230 $\pm$ 0.011 | 0.021 $\pm$ 0.018 | 0.143 $\pm$ 0.009 | 0.172 $\pm$ 0.039 |
| DRCFR | 0.142 $\pm$ 0.005 | 0.122 $\pm$ 0.017 | 0.218 $\pm$ 0.021 | 0.048 $\pm$ 0.032 |
| SCIGAN | 0.144 $\pm$ 0.005 | 0.112 $\pm$ 0.025 | 0.220 $\pm$ 0.026 | 0.049 $\pm$ 0.034 |
| DESCN | 0.192 $\pm$ 0.029 | 0.098 $\pm$ 0.029 | 0.143 $\pm$ 0.011 | 0.065 $\pm$ 0.046 |
| ESCFR | 0.202 $\pm$ 0.023 | 0.086 $\pm$ 0.028 | 0.145 $\pm$ 0.011 | 0.076 $\pm$ 0.045 |
| Ours | 0.112 $\pm$ 0.019 | 0.018 $\pm$ 0.012 | 0.131 $\pm$ 0.030 | 0.053 $\pm$ 0.039 |

Sensitive Analysis and Ablation Study. We explore which value of $\lambda_{1}$, $\lambda_{2}$ and $\rho$ can achieve the best performance. The results are demonstrated in Table [7](https://arxiv.org/html/2510.16419v1#A5.T7). We can see that our model is not sensitive to the change of hyperparameters. That is, the performance of CATE estimation remains relatively stable across a range of hyperparameters. For the ablation study presented in Table [8](https://arxiv.org/html/2510.16419v1#A5.T8), as that in IHDP and Twins, taking $\mathcal{L}_{\text{ce}}$ off only causes a moderate decline, while removing $\mathcal{L}_{\text{cosnt}}$ brings a whole disaster.

**Table 7: Sensitivity analysis on the Jobs dataset with respect to $\lambda_{1}$, $\lambda_{2}$, and $\rho$.**
| $\lambda_{1}$ | $\lambda_{2}$ | $\rho$ |  |  |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Value | $\mathcal{R}_{\text{pol}}^{\text{in}}$ | $\epsilon_{\text{ATT}}^{\text{in}}$ | $\mathcal{R}_{\text{pol}}^{\text{out}}$ | $\epsilon_{\text{ATT}}^{\text{out}}$ | Value | $\mathcal{R}_{\text{pol}}^{\text{in}}$ | $\epsilon_{\text{ATT}}^{\text{in}}$ | $\mathcal{R}_{\text{pol}}^{\text{out}}$ | $\epsilon_{\text{ATT}}^{\text{out}}$ | Value | $\mathcal{R}_{\text{pol}}^{\text{in}}$ | $\epsilon_{\text{ATT}}^{\text{in}}$ | $\mathcal{R}_{\text{pol}}^{\text{out}}$ | $\epsilon_{\text{ATT}}^{\text{out}}$ |
| 0.1 | 0.113 | 0.018 | 0.132 | 0.054 | 0.1 | 0.109 | 0.021 | 0.129 | 0.056 | 10 | 0.124 | 0.020 | 0.141 | 0.051 |
| 0.5 | 0.113 | 0.022 | 0.130 | 0.058 | 0.5 | 0.109 | 0.020 | 0.128 | 0.054 | 50 | 0.112 | 0.019 | 0.131 | 0.052 |
| 1 | 0.112 | 0.018 | 0.131 | 0.053 | 1 | 0.112 | 0.018 | 0.131 | 0.053 | 100 | 0.112 | 0.018 | 0.131 | 0.053 |
| 2 | 0.115 | 0.020 | 0.135 | 0.054 | 2 | 0.117 | 0.019 | 0.135 | 0.050 | 200 | 0.115 | 0.020 | 0.132 | 0.053 |
| 10 | 0.121 | 0.027 | 0.140 | 0.060 | 10 | 0.123 | 0.027 | 0.144 | 0.060 | 1000 | 0.114 | 0.020 | 0.133 | 0.052 |

**Table 8: ablation studies Jobs**
| Training Loss | $\mathcal{R}_{\text{pol}}^{\text{within-s.}}$ | $\epsilon_{\text{ATT}}^{\text{within-s.}}$ | $\mathcal{R}_{\text{pol}}^{\text{out-of-s.}}$ | $\epsilon_{\text{ATT}}^{\text{out-of-s.}}$ |
| --- | --- | --- | --- | --- |
| $\mathcal{L}_{\text{wls}}$ & $\mathcal{L}_{\text{const}}$ | 0.114 | 0.023 | 0.134 | 0.053 |
| $\mathcal{L}_{\text{wls}}$ & $\mathcal{L}_{\text{ce}}$ | 0.121 | 0.029 | 0.141 | 0.055 |
| Full (Ours) | 0.112 | 0.018 | 0.131 | 0.053 |

### E.4 Extended Sensitive Analysis

In this section we present the results of the sensitive analysis of hyperparameter $\lambda_{1}$ and $\rho$ in the IHDP and Twins dataset. One can see from Table [9](https://arxiv.org/html/2510.16419v1#A5.T9) and Table [10](https://arxiv.org/html/2510.16419v1#A5.T10) that our model is robust to the change of $\lambda_{1}$ and $\rho$, remaining good performance in the CATE estimation as well as relative error prediction.

**Table 9: Sensitivity analysis of $\lambda_{1}$ on IHDP and Twins datasets.**
| IHDP | Twins |  |  |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Value | $\sqrt{\epsilon_{\text{PEHE}}^{\text{in}}}$ | $\epsilon_{\text{ATE}}^{\text{in}}$ | $\sqrt{\epsilon_{\text{PEHE}}^{\text{out}}}$ | $\epsilon_{\text{ATE}}^{\text{out}}$ | Coverage | Selection | Value | $\sqrt{\epsilon_{\text{PEHE}}^{\text{in}}}$ | $\epsilon_{\text{ATE}}^{\text{in}}$ | $\sqrt{\epsilon_{\text{PEHE}}^{\text{out}}}$ | $\epsilon_{\text{ATE}}^{\text{out}}$ | Coverage | Selection |
| 0.1 | 0.678 | 0.096 | 0.709 | 0.112 | 0.93 | 0.74 | 0.1 | 0.286 | 0.009 | 0.288 | 0.010 | 0.96 | 0.94 |
| 0.25 | 0.693 | 0.096 | 0.724 | 0.113 | 0.93 | 0.75 | 0.25 | 0.285 | 0.010 | 0.287 | 0.010 | 0.94 | 0.94 |
| 0.5 | 0.638 | 0.090 | 0.670 | 0.105 | 0.96 | 0.80 | 0.5 | 0.284 | 0.009 | 0.286 | 0.009 | 0.94 | 0.94 |
| 1 | 0.712 | 0.103 | 0.746 | 0.115 | 0.96 | 0.79 | 1 | 0.285 | 0.013 | 0.287 | 0.014 | 0.94 | 0.92 |
| 2.5 | 1.011 | 0.245 | 1.036 | 0.262 | 0.94 | 0.77 | 2.5 | 0.283 | 0.015 | 0.284 | 0.016 | 0.92 | 0.88 |

**Table 10: Sensitivity analysis of $\rho$ on IHDP and Twins datasets.**
| IHDP | Twins |  |  |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Value | $\sqrt{\epsilon_{\text{PEHE}}^{\text{in}}}$ | $\epsilon_{\text{ATE}}^{\text{in}}$ | $\sqrt{\epsilon_{\text{PEHE}}^{\text{out}}}$ | $\epsilon_{\text{ATE}}^{\text{out}}$ | Coverage | Selection | Value | $\sqrt{\epsilon_{\text{PEHE}}^{\text{in}}}$ | $\epsilon_{\text{ATE}}^{\text{in}}$ | $\sqrt{\epsilon_{\text{PEHE}}^{\text{out}}}$ | $\epsilon_{\text{ATE}}^{\text{out}}$ | Coverage | Selection |
| 10 | 0.698 | 0.108 | 0.735 | 0.123 | 0.96 | 0.78 | 10 | 0.299 | 0.015 | 0.306 | 0.015 | 0.92 | 0.62 |
| 50 | 0.711 | 0.098 | 0.745 | 0.116 | 0.95 | 0.79 | 50 | 0.289 | 0.011 | 0.291 | 0.012 | 0.90 | 0.88 |
| 100 | 0.638 | 0.090 | 0.670 | 0.105 | 0.96 | 0.80 | 100 | 0.284 | 0.009 | 0.286 | 0.009 | 0.94 | 0.94 |
| 200 | 0.737 | 0.103 | 0.772 | 0.123 | 0.94 | 0.76 | 200 | 0.286 | 0.012 | 0.288 | 0.013 | 0.94 | 0.92 |
| 1000 | 0.751 | 0.111 | 0.785 | 0.130 | 0.93 | 0.76 | 1000 | 0.284 | 0.010 | 0.285 | 0.011 | 0.92 | 0.94 |

### E.5 Model Implementation

We implement all models using PyTorch and optimize them with the Adam optimizer. The key hyperparameters include the size of each hidden layer, learning rate, the loss coefficients $\lambda_{1}$, $\lambda_{2}$, the penalty coefficient $\rho$, and the number of training epochs. These hyperparameters are manually tuned through empirical trials. The search ranges are as follows: hidden layer size in $\{30,40,50,60,70\}$, learning rate in $\{5\times 10^{-4},10^{-3},2\times 10^{-3},3\times 10^{-3}\}$, $\lambda_{1},\lambda_{2}$ in $\{0.1,0.25,0.5,1,2\}$, $\rho$ in $\{10,50,100,200\}$, and number of training epochs in $\{700,800,900,1000,1100\}$.

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. What does the robust evaluation exactly mean? After reading this paper, I guess the author wants to claim the relative error estimator is "doubly robust", so it is robust to nuisance parameter estimation. 

2. What's the difference between your loss and R-loss (R-loss is also doubly robust)?

3. The problem of robust evaluation of the CATE estimator has been studied in [1], where they also consider the worst-case performance of the CATE estimator selected by the proposed evaluation metric. Does this paper also provide any information on the worst-case performance? I guess Figure 1 provides such information, and I suggest emphasizing this point as it can strengthen the efficacy of "robust".

4. A brainstorm: I think the whole framework can be extended to policy/estimator adaptation, e.g., following the setting in [2]. Maybe we can design a new policy evaluation metric that is doubly robust to nuisance parameters, which is very useful when distribution shift presents, and it has not been discussed in previous literature. So this is also a good point of this paper, as it has the potential to be extended to other problems.

[1] Unveiling the Potential of Robustness in Selecting Conditional Average Treatment Effect Estimators

[2] Optimal Policy Adaptation under Covariate Shift
soundness: 3
strengths: 1. The idea of solving the fundamental and practical CATE evaluation problem is interesting. Evaluating CATE is not easy unless we impose additional structural assumptions. Instead of just using a standard doubly-robust estimator, the authors first derive the specific theoretical conditions required for their estimator to be robust to outcome model misspecification, and then they design a new loss function and a neural network architecture to force the nuisance parameter estimates to satisfy these conditions. This is, somehow, follows the idea of the design of Dragonnet but gives an interesting solution for CATE evaluation.

2. Theoretical results are solid. They establish a theoretical asymptotic property for their proposed estimator, under some assumptions (conditions).

3. The empirical evaluation is through. The experiments strongly support the paper's claims, showing that the proposed method achieves the target 90% confidence interval coverage and high "selection accuracy"
summary: This paper aims to address a very important and critical challenge of CATE estimator evaluation. The authors design a novel loss function and neural network architecture that produces a robust relative error estimate, which remains root-n-consistent.
weaknesses: 1. Can strengthen the connection with other studies (see below question 2&3).

2. The asymptotic property relies on some critical assumptions, for example, $\check{\gamma}, \check{\beta}_0, \check{\beta}_1$ should converge to the true one at a rate faster than $n^{-1/4}$. It is strong but also reasonable. It is reasonable because orthogonal ML literature always assumes the convergence rate of nuisance parameters. It is strong because we are treating them as "plug-in" quantities, instead of the nuisance parameters that the estimator should be doubly robust to. It might be useful to provide some explanations to justify that the condition is realistic.

Overall, I think these are minor weaknesses. The pros outweigh the cons.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: - Is there any theoretical or empirical evidence of what happens when the propensity score estimation model is misspecified?
- How do results depend on the number of candidate estimators?
soundness: 4
strengths: - This paper addresses a very important yet understudied problem, and an effective practical solution for the evaluation of HTE estimators is very valuable
 - The paper offers significant theoretical contributions offering theoretical guarantees on estimator consistency even with misspecified regression models
 - In addition to theoretical results on the evaluation of HTE estimators, the paper offers a learning algorithm and method for HTE estimation
 - The empirical evaluation on datasets available in the literature is comprehensive
summary: This paper extends the work of (Gao 2025) on the use of relative error for heterogeneous treatment effect estimation. Its main contribution is that it offers a relaxation of the requirement that the nuisance parameter estimators must be consistent. Instead, it is shown that robust estimation can be achieved if the a correctly specified propensity score estimator is provided. Moreover, the theoretical results lead to the development of an HTE method, while both theoretical and experimental results support the arguments of the paper.
weaknesses: - The need for a correctly specified propensity score model still remains a strong assumption, which may not hold or be guaranteed in practice in many cases.
 - One would expect a deeper evaluation on the comparison of the proposed method against the method by Gao 2025, as well as an empirical evaluation of what happens when the propensity score estimator is misspecified.
 - There is some confusion regarding the terminology, since there is mention of "balance regularizers" in the beginning but later this terminology is abandoned, where $\mathcal{L}_{\text{const}}$ is mentioned

---

confidence: 4
contribution: 2
initial_rating: 2
presentation: 1
questions: 1. Your asymptotic CI relies on the correct specification of propensities. What happens when this assumption is violated?
2. How is the ground truth $\tau(x)$ constructed for Twins and Jobs
3. Why does the method not require sample-splitting, a standard procedure in orthogonal estimators?
4. Theorem 1 is only correct when Eq.(4) is satisfied, right? But in your method, you cannot make Eq.4 satisfied. Is there any bound analysis?
5. Why is propensity score easier to estimate than the outcome? (line 172 - 175)
6. The complexity of the method seems high. What is the average time cost of achieving an enhanced estimator using your method? And how does it compare to the time cost of training a baseline CATE estimator?
soundness: 2
strengths: - The paper successfully identifies a problem, which is that the existing estimators require asymptotic rates of nuisances no less than n^-¼. 
- Through Taylor expansion, the paper shows the conditions under which higher error rate of \mu does not affect the asymptotic rate of n^-½ of the relative error estimator. This provides interesting theoretical insights.
- The paper designs a novel loss that reformulates the equation-based condition as a minimization problem which allows batch gradient descent.
- The paper shows that the method has practical significance in getting stronger CATE models.
summary: The paper proposes a way of making an estimator of the relative error of two given CATE estimators robust against the misspecification of the mean outcome nuisance \mu_a, relaxing the constraints of existing methods. Furthermore, the by-product of their algorithm could be used to craft a stronger CATE estimator.
weaknesses: 1. The paper has limited novelty compared with existing DR estimator Gao (2025). The fundamental form of estimator remains the same as Gao. And the paper seems to just reconstruct the nuisances to meet certain constraints. 
2. The parametrization of e(X) and \mu_a(X) in (1),(2) seem quite arbitrary and lack a decent justification on why these nuisances should share representation. 
3. There is NO solution to the moment-based conditions. To this, soft penalties are introduced for finite-sample optimization. Without the correct specification of propensity (which is never met in real world applications) or sensitivity analysis on the propensity score, I doubt the complete formulation might just go meaningless.
4. Misspecification can also come from the representation $\Phi(X)$.
5. Nuisance overfit could be a problem without sample splitting.
6. Typo: line 1017 should be L_const.
7. Experiment: 1) The paper does not perform sensitivity analysis (empirical) on the propensity score, which makes the practical utility doubtful. 2) The paper does not explain how they compute the \tau in real world datasets Twins and Jobs.3) The paper does not provide convergence analysis of the losses, especially L_const and L_ce. The whole method would only make sense if the paper could show, at least empirically, that under finite samples the two losses successfully approximate the original condition.

---

confidence: 5
contribution: 3
initial_rating: 8
presentation: 3
questions: - The paper focuses on HTE estimation but only provides formal definitions for ITE and CATE in the problem setup. For clarity, please provide a precise definition of the HTE you are targeting.

- There appears to be a typo in Assumption 1(i). It should likely be the standard unconfoundedness assumption: $(Y_i(0), Y_i(1)) \perp A_i \mid X_i$. Please clarify and correct this.

- Choice of Evaluation Objective: Line 109 states that the goal is to select the estimator with the highest accuracy on a given test dataset. Could the authors elaborate on why this objective was chosen over the more common goal of finding an estimator that accurately models the CATE function over the marginal distribution of covariates, $P(X)$?

- Example 1 illustrates model misspecification, which is a well-understood concept in statistical machine learning. To improve conciseness, this example could be simplified or moved to the appendix.

- The neural network architecture described in Section 4.3 appears to be the Dragonnet architecture. Could the authors clarify the novelty of their proposed network, or explicitly frame it as an adaptation of Dragonnet for their specific loss functions?

- Experimental Comparison: The paper's main theoretical contribution is relaxing the assumptions of Gao (2025). It would be valuable to include a baseline that represents the Gao (2025) approach or can you clarify how the proposed estimator outperform it?

- Minor Typos: Line 173: "violating Assumption 2" should likely refer to "violating Condition 2". Line 216: The text seems to have a typo and should probably read "...the proposed estimator of... ".
soundness: 4
strengths: This is a solid paper that addresses an important and interesting problem. The main strengths are:

- Valuable Problem: The paper focuses on the reliable evaluation of HTE estimators, which is a critical issue for real-world applications. The idea of using "relative error" to compare models is novel and meaningful.

- Strong Theoretical Contribution: The highlight of this paper is its theoretical work. It successfully relaxes the strict assumptions from prior work (Gao, 2025), notably by removing the requirement for the outcome model to be correctly specified. This makes the proposed framework much more robust and practical for real-world scenarios, which is a significant and elegant theoretical improvement.

- Clear and Solid Method: The paper is very clearly written and easy to follow. The proposed new loss function and network architecture are well-motivated and tightly connected to the theoretical derivations. The overall approach feels very solid.

- Thorough Experiments: The authors provide convincing results from experiments on several standard datasets. The experimental design is comprehensive, including ablation and sensitivity analyses, which effectively validate the proposed method's effectiveness.
summary: This paper introduces a new *evaluation framework* for HTE estimators, extending a recently proposed relative error–based framework that compares the performance of two estimators. The proposed approach relaxes a key assumption of the prior method, which required all nuisance parameter estimators to be consistent at a rate faster than n-1/4. In contrast, the new framework only requires the propensity score model to meet this condition, allowing the outcome regression and propensity model to be misspecified. Building on this evaluation framework, the authors further develop a new HTE estimation method and evaluate it on some very often used semi-synthetic datasets, including IHDP, Twins and Jobs, demonstrating improved performance compared to several baseline estimators. The paper also provides sensitivity and ablation analyses to support the empirical findings.
weaknesses: - The paper's primary weakness is the unclear practical motivation for using "relative error." While the metric is theoretically interesting for comparing two estimators, its real-world applicability is not well-established. The paper would be significantly stronger if it could provide clear use cases where practitioners would prefer this comparative measure over standard absolute performance metrics like PEHE. This is particularly important for the ICLR audience, which values practical impact.

- The paper is mathematically dense, and while this contributes to its theoretical rigor, it also appears to have led to several typos and minor inconsistencies. To ensure the paper's core contributions are communicated accurately, a thorough proofreading of all definitions, assumptions, and equations is highly recommended.


Decision: Accept (Poster)