# Estimating Long-term Heterogeneous Dose-response Curve: Generalization Bound
Leveraging Optimal Transport Weights

Zeqin Yang1 Weilin Chen1∗ Ruichu Cai1,2 Yuguang Yan1 Zhifeng Hao3
Zhipeng Yu4 Zhichao Zou4 Peng Zhen4 Jiecheng Guo4
1 School of Computer Science, Guangdong University of Technology
2 Pazhou Laboratory (Huangpu)
3 College of Science, Shantou University
4 Didi Chuxing
<EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS>Equal
contributionsCorresponding authors.

###### Abstract

Long-term causal effect estimation is a significant but challenging problem in
many applications. Existing methods rely on ideal assumptions to estimate
long-term average effects, e.g., no unobserved confounders or a binary
treatment, while in numerous real-world applications, these assumptions could
be violated and average effects are unable to provide individual-level
suggestions. In this paper, we address a more general problem of estimating
the long-term heterogeneous dose-response curve (HDRC) while accounting for
unobserved confounders. Specifically, to remove unobserved confounding in
observational data, we introduce an optimal transport weighting framework to
align the observational data to the experimental data with theoretical
guarantees. Furthermore, to accurately predict the heterogeneous effects of
continuous treatment, we establish a generalization bound on counterfactual
prediction error by leveraging the reweighted distribution induced by optimal
transport. Finally, we develop an HDRC estimator building upon the above
theoretical foundations. Extensive experimental studies conducted on multiple
synthetic and semi-synthetic datasets demonstrate the effectiveness of our
proposed method.

##  1 Introduction

Long-term causal inference is practically significant across many domains [27,
19, 13]. Randomized control experiments are the gold standard for estimating
causal effects. However, it is practically difficult or even impossible to
collect the long-term outcome of interest through experiments due to cost
considerations, e.g., taking years or even decades of follow-up to fully
reveal mortality in acquired immunodeficiency syndrome drug clinical trials.
In contrast, long-term observational data are often easier or cheaper to
acquire, and thus the long-term observational data has been taken as the
complementary of short-term experimental data. A typical causal graph of long-
term effect estimation is shown in Fig. 1(a), where short-term experimental
and long-term observational data are accessible.

There has been considerable interest in using short-term outcomes as
replacement or supplement for the long-term outcome in different ways via data
combination. With all observed confounders, [32] focuses on a binary treatment
setting, showing the benefit of short-term outcomes in estimating average
long-term effects. Further, [60] extends it to the continuous treatment
setting. As a counterpart, allowing the existence of unobserved confounders,
[3, 28, 30] assume that short-term outcomes, to some extent, capture the
information of unobserved confounders to deal with the confounding problem,
nevertheless, they still all focus on estimating average effect of binary
treatment.

(a) Experimental and Observational data

Works |  Continuous Treatment |  Unobserved Counfounder |  Heterogeneous Effect
---|---|---|---
Kallus and Mao [32] | ✗ | ✗ | ✗
Zeng et al. [60] | ✔ | ✗ | ✗
Athey et al.[3] | ✗ | ✔ | ✗
Hu et al. [28] | ✗ | ✔ | ✗
Imbens et al. [30] | ✗ | ✔ | ✗
Ours | ✔ | ✔ | ✔

(b) Setting Comparison

Figure 1: Fig. 1(a) are typical causal graphs of experimental and
observational data in long-term effect estimation scenarios. Fig. 1(a) shows a
ride-hailing scenario where treatment $A$ is continuous and unobserved
confounders $\mathbf{U}$ exist. Table 1(b) shows the differences between
existing works and ours.

Previous works, however, lack practical applicability, since they have largely
focused on estimating average treatment effects under ideal assumptions, e.g.,
no unobserved confounders in observational data or binary treatment, as shown
in Table 1(b). In many long-run applications, such assumptions may be
violated. Decision-makers usually leverage long-term individual-level effects
estimated with unobserved confounders and continuous treatment to develop
policies for individuals. For the example shown in Fig. 1(a), in a ride-
hailing platform, in order to increase drivers’ long-term retention (long-term
outcome), the platform has to evaluate the effect of various incomes
(continuous treatment) of different drivers (heterogeneity). For effect
estimation, the platform can conduct short-term experiments to obtain drivers’
short-term behaviors (short-term outcomes) and collect historical long-term
observational data. In the observational data, even though short/long-term
outcomes can all be observed, confounders are partially observed due to
information privacy, e.g., drivers’ family backgrounds (unobserved
confounders). In such an application, with unobserved confounders, it is still
understudied how to accurately estimate the long-term effect of various
incomes on the long-term retention of specific drivers, i.e., heterogeneous
dose-response curve (HDRC), via data combination.

Motivated by the example above, in this paper, we focus on estimating long-
term HDRC with unobserved confounders via data combination as shown in Table
1(b). Long-term heterogeneous effect estimation poses significant challenges,
primarily due to the unobserved confounders in the observational dataset and
the lack of counterfactual outcomes in estimating the heterogeneous curves of
continuous treatment. Specifically, the first and biggest challenge is that
the long-term outcome only exists in the observational dataset with unobserved
confounders and is missing in the experimental dataset, causing the
identification problem. The second challenge is that since treatment is
continuous, the short and long-term counterfactual outcomes are potentially
infinite, causing a large counterfactual prediction error in estimating the
heterogeneous curves with only the actual outcome. Our solution and
contribution can be summarized as follows:

  * •

We address the problem of estimating long-term heterogeneous dose-response
curves with unobserved confounders via data combination. To the best of our
knowledge, this is the first work to address this problem as shown in Table
1(b).

  * •

We theoretically show that the unobserved confounding can be removed by
aligning the conditional distributions of short-term outcomes between
observational and experimental groups. To achieve such alignment, under a
reweighting framework of optimal transport (OT), we establish a mini-batch and
joint distribution-based upper bounds on the discrepancy of conditional
distributions

  * •

Based on the OT-induced weighted distribution, we derive the generalization
bound on the counterfactual prediction error in estimating long-term HDRC,
which inspires our final model, called Long-term hEterogeneous dose-response
curve estimAtor with Reweighting and represeNtation learning (LEARN).

##  2 Related Work

Optimal Transport Recently, optimal transport has shown powerful abilities in
different kinds of applications, such as computer vision [26, 21, 38, 48],
domain adaption [14, 15, 17, 10], generative models [29, 22, 37, 1] and so on.
There are also some works trying to introduce optimal transport into causal
inference. [25] employs unbalanced optimal transport for matching. [39]
proposes to infer counterfactual outcomes via transporting the factual
distribution to the counterfactual one. [18, 58] apply optimal transport to
learn sample weights to achieve distribution balancing. [55] solves the
problems of mini-batch sampling and unobserved confounders in CFR [46] through
optimal transport. To the best of our knowledge, this is the first work
exploiting optimal transport to solve the unobserved confounder problem via
data combination under the long-term estimation scenario.

Long-term Causal Inference To date, a common approach to estimating long-term
causal effects is based on short-term outcomes. By making an additional
surrogacy assumption that short-term outcomes can mediate the full effect of
treatment, [2] combines multiple short-term outcomes to robustly estimate the
average long-term effects, [12] adapt RNN to learn the relationship between
the short-term and long-term outcomes. However, the surrogacy assumption can
be easily violated in practice. In this setting, [32] investigates the role of
short-term outcomes for long-term causal inference and [60] serves as a
counterpart for [32] in the continuous treatment setting. [9] considers a more
realistic scenario where some valid surrogates are unobserved and employ an
identifiable variational autoencoder to recover them. Further, To solve the
unobserved confounder problem, [3] proposes a relaxed latent unconfoundedness
assumption, and proposes different estimators based on the imputation,
weighting, and control function methods. [23] proposes other feasible
assumptions to replace the assumption used in [3]. [30] exploits the
sequential structure of short-term outcomes and develops a proximal learning-
based method. [50] performs instrumental variable regression to estimate
average long-term effects. Although many different settings of long-term
effect estimation have been explored, to the best of our knowledge, we are the
first to focus on the long-term HDRC estimation problem with unobserved
confounders.

##  3 Preliminary

###  3.1 Optimal Transport Weighting

Optimal transport (OT) [52] is widely used to quantify distribution
discrepancy as the minimum cost of transporting one distribution to another.
Among the rich theory of OT, we focus on the Kantorovich Problem [33] for
discrete distribution in this paper. Recently, thanks to its many advantages,
such as geometric sensitivity and stability, OT has been explored for learning
weights to align distribution [18, 58, 26] and demonstrated excellent
performance. Specifically, suppose we want to learn weights
$\mathbf{w}_{\alpha}$ for samples in distribution $\alpha$, such that the
reweighted distribution is aligned with another distribution $\beta$. To
achieve this goal, we view the to-be-learned weights $\mathbf{w}_{\alpha}$ as
the probability measure of distribution $\alpha$, and represent the
distribution $\beta$ with a uniform probability measure $\bm{\mu}$. By doing
so, the learning of $\mathbf{w}_{\alpha}$ can be formulated as the process of
minimizing the OT distance between distributions $\alpha$ and $\beta$:

$\displaystyle\min_{\mathbf{w}_{\alpha}}OT(\alpha,\beta)=\min_{\mathbf{w}_{\alpha}}\min_{\Gamma\in\Pi(\alpha,\beta)}\langle\Gamma,C\rangle,$
$\displaystyle s.t.~{}~{}\Pi(\alpha,\beta)=\\{\Gamma\in{\mathbb{R}}^{n\times
m}|\Gamma\mathbf{1}_{m}=\mathbf{w}_{\alpha},\Gamma^{T}\mathbf{1}_{n}=\bm{\mu},\Gamma_{ij}\in\left[0,1\right]\\}$
(1)

where $\langle\cdot,\cdot\rangle$ is the Frobenius inner product, $n$ and $m$
are the sample numbers of distribution $\alpha$ and $\beta$, the transport
cost matrix $C$ is the unit-wise distance between $\alpha$ and $\beta$, the
transport probability matrix $\Gamma$ satisfying $\Pi(\alpha,\beta)$ is
learned by minimizing $OT(\alpha,\beta)$, which reflects how to transport
samples from the distribution $\alpha$ to $\beta$.

###  3.2 Long-term Heterogeneous Dose-response Curve

Let $A\in{\mathcal{A}}$ be the 1-dimensional continuous treatment,
$\mathbf{X}\in{\mathcal{X}}$ be the observed covariates,
$\mathbf{U}\in{\mathcal{U}}$ be the unobserved covariates,
$\mathbf{S}(a)\in{\mathcal{S}}$ be the potential short-term outcomes measured
at timesteps $1,2,...,t_{0}$, and $Y(a)\in{\mathcal{Y}}$ be the potential
long-term outcome measured at timestep $T$. Let lowercase letters (e.g.,
$a,\mathbf{x},\mathbf{u},\mathbf{s}(a),y(a)$) denote the value of the above
random variables. Following potential outcome framework [44], the observed
short and long-term outcomes $\mathbf{S},Y$ are the potential outcomes
$\mathbf{S}(a),Y(a)$ corresponding to the actually received treatment $a$.

Following [3, 11, 23, 30], we consider a data combination setting where
$G=\\{o,e\\}$ distinguishes between two types of datasets: a large
observational dataset
$O=\\{a_{i},\mathbf{x}_{i},\mathbf{s}_{i},y_{i},G_{i}=o\\}_{i=1}^{n_{o}}$ and
a small experimental dataset
$E=\\{a_{i},\mathbf{x}_{i},\mathbf{s}_{i},G_{i}=e\\}_{i=1}^{n_{e}}$ ($n_{e}\ll
n_{o}$). That is, the treatment $A$, covariates $\mathbf{X}$ and short-term
outcomes $\mathbf{S}$ are available in both datasets, while long-term outcome
$Y$ is only available in observational dataset. Our goal is to estimate the
long-term heterogeneous dose-response curve (HDRC), which is defined as:

$\displaystyle\mu(\mathbf{x},a)={\mathbb{E}}[Y(a)|\mathbf{X}=\mathbf{x}].$ (2)

Notably, the long-term HDRC cannot be estimated solely from the experimental
dataset due to the missingness of long-term outcome $Y$, and also cannot be
estimated solely from the observational dataset due to the possible unobserved
confounders $\mathbf{U}$. To overcome these challenges, we learn long-term
HDRC via data combination, requiring the following assumptions:

######  Assumption 1 (Consistency, Positivity).

If $A=a$, we have $Y=Y(a),\mathbf{S}=\mathbf{S}(a)$. And $\forall a,x$, we
have
$0<P(A=a|\mathbf{X}=\mathbf{x})<1,~{}0<P(G=o|A=a,\mathbf{X}=\mathbf{x})<1$.

######  Assumption 2 (Weak internal validity of observational data).

$A\perp\\!\\!\\!\perp\\{Y(a),\mathbf{S}(a)\\}|\mathbf{X},\mathbf{U},G=o$.

######  Assumption 3 (Internal validity of experimental data).

$A\perp\\!\\!\\!\perp\\{Y(a),\mathbf{S}(a)\\}|\mathbf{X},G=e$.

######  Assumption 4 (External validity of experimental data).

$G\perp\\!\\!\\!\perp\\{Y(a),\mathbf{S}(a)\\}|\mathbf{X}$.

######  Assumption 5 (Latent Unconfoundedness).

$A\perp\\!\\!\\!\perp Y(a)|\mathbf{X},\mathbf{S}(a),G=o$.

The assumptions above are mild and widely used in existing literature, e.g.,
[3, 23]. Assumption 1 is a standard and basis assumption. Importantly,
Assumption 2 allows unobserved confounders to exist in the observational
dataset, and Assumption 3 guarantees that the treatment in the experimental
dataset is unconfounded after controlling the covariates. Assumption 4 allows
us to generalize the conditional distribution of potential outcomes between
different groups $G$. Assumption 5 is proposed by [3] to achieve average long-
term effect identification under a binary-treatment setting. Different from
them, we reuse this assumption to identify HDRC by a novel reweighting method.

##  4 Proposed Method

To estimate long-term HDRC, we present our method exploiting reweighting to
handle the unobserved confounders and representation learning to handle the
observed confounders, respectively. In section 4.1, we first propose a novel
weighting technique that can deal with the unobserved confounders in the
observational dataset, making HDRC identifiable. We then provide a practical
and computationally efficient OT method with theoretical guarantees to learn
such weights via data combination. In section 4.2, we derive the
generalization bound on counterfactual regression error to address the bias
caused by the observed confounders based on the reweighted observational
distribution.

###  4.1 Learning Optimal Transport Weights for Unobserved Confounders via
Data Combination

####  4.1.1 Identifiable Long-term HDRC via Reweighting

We first study how to identify the long-term HDRC. Given that $Y$ is only
observed in the observational dataset, identifying HDRC requires the
unconfoundedness assumption in the observational dataset, i.e.,
$A\perp\\!\\!\\!\perp\\{\mathbf{S}(a),Y(a)\\}|\mathbf{X},G=o$, which is
unsatisfied as stated in Assumption 2. Fortunately,
$A\perp\\!\\!\\!\perp\\{\mathbf{S}(a),Y(a)\\}|\mathbf{X},G=e$ holds, though
$Y$ is missing in the experimental dataset. It raises a question: can we
remove the unobserved confounding in the observational dataset to achieve
identification, by utilizing the experimental dataset which includes only
short-term outcomes? We find that, even though $\mathbf{U}$ is unobserved, its
influence can be captured by the (in)dependence between $\mathbf{S}$ and $G$
conditional on $\mathbf{X}$ and $A$, as illustrated in the following
proposition.

######  Proposition 1.

Under Assumptions 1, 2,3 4, and 5, given a set of weights
$\mathbf{w}=\\{\mathbf{w}_{o},~{}\bm{\mu}\\}$ consisting of the learnable
weights $\mathbf{w}_{o}$ for observational units and uniform weights
$\bm{\mu}$ for experimental units, which makes
$\mathbf{w}P(\mathbf{S},G|\mathbf{X},A)=P^{\mathbf{w}}(\mathbf{S},G|\mathbf{X},A)=P^{\mathbf{w}}(\mathbf{S}|\mathbf{X},A)P^{\mathbf{w}}(G|\mathbf{X},A)$,
i.e., $\mathbf{S}\perp\\!\\!\\!\perp G|\mathbf{X},A$, then
$\mathbb{E}[\mathbf{S}(a)|\mathbf{X},A=a,G=o]=\mathbb{E}[\mathbf{S}(a)|\mathbf{X},G=o]$
holds on the weighted observational distribution.

The proof is given in Appendix A.1. Proposition 1 is intuitive. Since
unconfoundedness holds in $G=e$, then
$P(\mathbf{S}(a)|\mathbf{X},G=e)=P(\mathbf{S}|\mathbf{X},A=a,G=e)$, but it
does not hold when $G=o$. If we can learn weight $\mathbf{w}$ aligning short-
term outcomes between observational data and experimental data, i.e., making
$P(\mathbf{S}|\mathbf{X},A,G=e)=P(\mathbf{S}|\mathbf{X},A,G=o)$, the
unconfoundedness condition will also hold in $G=o$ based on Assumption 4.
Proposition 1 provides an effective way to remove the unobserved confounding
in observational dataset $G=o$. Based on Proposition 1, we can further show
HDRC is identifiable as follows.

######  Theorem 1.

Suppose Assumptions in Proposition 1 hold, we have the unconfoundedness
$\mathbb{E}[Y(a)|\mathbf{X},A=a,G=o]=\mathbb{E}[Y(a)|\mathbf{X},G=o]$, then
the long-term HDRC can be identified.

The proof is given in Appendix A.2. Theorem 1 is built on Proposition 1,
indicating that if we could learn suitable weights making
$\mathbf{S}\perp\\!\\!\\!\perp G|\mathbf{X},A$, the unobserved confounding
will be removed, and then long-term HDRC can be identified. Therefore, our
goal now is to learn a set of weights $\mathbf{w}_{o}$ for observational units
such that the conditional independence $\mathbf{S}\perp\\!\\!\\!\perp
G|\mathbf{X},A$ holds based on the reweighted distribution, i.e.,
$P^{\mathbf{w}_{o}}(\mathbf{S}|\mathbf{X},A,G=o)=P^{\bm{\mu}}(\mathbf{S}|\mathbf{X},A,G=e)$.

####  4.1.2 Learning Optimal Transport Weights

Based on the above theoretical analyses, we need learn weights
$\mathbf{w}_{o}$ for observational units to minimize the distances between
$P(\mathbf{S}|\mathbf{X},A,G=o)$ and $P(\mathbf{S}|\mathbf{X},A,G=e)$.
However, most existing works study how to measure the distances between
(joint) distributions, instead of conditional distributions. Surprisingly, we
find that under the optimal transport framework, the distance between
conditional distributions can be bounded by that of joint distributions. As a
result, we propose a mini-batch joint distribution-based OT method that not
only helps align the above conditional distributions but also can be easily
embedded into the training of deep learning with theoretical guarantees.

To begin with, following the objective function in Eq. (3.1), we model our
goal as a problem of minimizing the conditional OT distance, i.e.,
$OT(P(\mathbf{S}|\mathbf{X},A,G=o),P(\mathbf{S}|\mathbf{X},A,G=e))$. The
intuitive approach to solve this problem is to minimize every sub-problem
$OT(P(\mathbf{S}|\mathbf{X}=\mathbf{x},A=a,G=o),P(\mathbf{S}|\mathbf{X}=\mathbf{x},A=a,G=e))$
for every realization $\mathbf{X}=\mathbf{x},A=a$ as follows:

$\displaystyle
OT^{con}_{\mathbf{x},a}=\sum_{\mathbf{x},a}OT(P(\mathbf{S}|\mathbf{X}=\mathbf{x},A=a,G=o),P(\mathbf{S}|\mathbf{X}=\mathbf{x},A=a,G=e)).$
(3)

The First Challenge However, solving the above conditional OT problem is
infeasible. The reason is that the conditional set $(\mathbf{X},A)$ in our
conditional OT distance are multi-dimensional continuous variables, leading to
the existence of infinitely sub-problems in Eq. (3). To overcome this
challenge, we establish a connection between the OT distance regarding
conditional distributions and that regarding joint distributions:

######  Theorem 2.

Assuming the cost matrix of the joint distribution
$P(\mathbf{S},\mathbf{X},A)$ is separable, i.e.,
$C(\mathbf{s},\mathbf{x},a;\tilde{\mathbf{s}},\tilde{\mathbf{x}},\tilde{a})=C(\mathbf{s};\tilde{\mathbf{s}})+C(\mathbf{x};\tilde{\mathbf{x}})+C(a;\tilde{a})$,
we have:

$\displaystyle OT^{con}_{\mathbf{x},a}\leq
OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e)).$ (4)

In addition, under the positivity assumption, we have
$OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))\to 0$ as
$n_{o}\to\infty$, leading to $OT^{con}_{\mathbf{x},a}\to 0$.

The proof can be found in Appendix A.3. Theorem 2 shows that the conditional
OT distance can be bounded by the OT distance between the corresponding joint
distributions, which motivates us to learn weights $\mathbf{w}_{o}$ by
minimizing the upper bound
$OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))$ instead of
the infeasible $OT^{con}_{\mathbf{x},a}$. And when $n_{o}\to\infty$, the
convergence property $OT^{con}_{\mathbf{x},a}\to 0$ means that the probability
measure $\mathbf{w}_{o}$ of the observational distribution
$P^{\mathbf{w}_{o}}(\mathbf{S}|\mathbf{X},A,G=o)$ will converge to that of
experimental distribution $P(\mathbf{S}|\mathbf{X},A,G=e)$, ensuring HDRC
identifiable based on Proposition 1.

The Second Challenge While we solve the problem of conditional OT, another
challenge emerges. Specifically, solving
$OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))$ will cause
a high computational cost since the whole optimized transport matrix
$\Gamma\in{\mathbb{R}}^{n_{o}\times n_{e}}$ at each iteration can be large
with the large size of observational dataset $n_{o}$. Additionally, optimizing
the whole matrix at each iteration is also not suitable for a deep learning
framework that is based on mini-batch training. To tackle this challenge,
inspired by [59], we propose a mini-batch OT based on a randomly sub-sampled
observational dataset and provide its theoretical analysis as follows:

######  Definition 1.

Consider two empirical distributions $\alpha$ and $\beta$ with $n$ and $m$
units, we assume the batch size $b$ of the first distribution $\alpha$ satisfy
$b\mid n$ and let $k=n/b$ be the number of batches. Let $\mathcal{B}_{i}$ be
index set of the $i$-th batch in $\alpha$ and the corresponding empirical
distribution is $\alpha_{\mathcal{B}_{i}}$, then the mini-batch OT problem is
defined as:

$\displaystyle
m\mbox{-}OT(\alpha,\beta)=\frac{1}{k}\sum_{i=1}^{k}OT(\alpha_{\mathcal{B}_{i}},\beta).$
(5)

######  Theorem 3.

Let $\gamma_{i}$ be the optimal transport probability matrix of the $i$-th
batch OT problem of
$m\mbox{-}OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))$,
i.e.,
$OT(P_{\mathcal{B}_{i}}(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))$.
With extending $\gamma_{i}$ to a $n_{o}\times n_{e}$ matrix $\Gamma_{i}$ that
pads zero entries to the row whose index does not belong to $\mathcal{B}_{i}$,
we have:

$\displaystyle\frac{1}{k}\sum_{i=1}^{k}\Gamma_{i}\in\Pi(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e)),$
(6)

and

$\displaystyle
OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))\leq
m\mbox{-}OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e)).$
(7)

The proof is in Appendix A.4. Theorem 3 implies that the OT distance of the
joint distribution $P(\mathbf{S},\mathbf{X},A)$ between observational and
experimental datasets is upper-bounded by its corresponding $m\mbox{-}OT$
problem. Solving the $m\mbox{-}OT$ problem not only significantly reduces the
high optimization cost, but also allows us to embed the OT distance into the
deep learning framework, where only batch observational units and full
experimental units are considered at each training iteration.

Conclusion We have overcome two challenges about the conditional OT distance
and high computational cost via theorem 2 and theorem 3 respectively.
Combining these two theorems, we could learn the weights $\mathbf{w}_{o}$ by
minimizing
$m\mbox{-}OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))$,
which is upper bound of the original $OT^{con}_{\mathbf{x},a}$. As a result,
our estimated weights $\mathbf{w}_{o}$ could achieve Proposition 1
approximately, which leads to the identification of HDRC according to Theorem
1.

###  4.2 Generalization Bound on Long-term HDCR Estimation Error

Ideally, learning the OT-induced weights $\mathbf{w}_{o}$ removes the
unobserved confounding and makes HDRC identifiable. However, in estimating
long-term HDRC, the counterfactual prediction error is still large. The reason
is that the regression can only be trained with one of the counterfactual
outcomes (factual outcomes) while the counterfactual outcomes are infinite due
to the continuous treatment and the confounding bias exists since
$\mathbf{X}\not\perp\\!\\!\\!\perp A$. To address it, in this section, we
derive the generalization bound on the counterfactual prediction error, based
on the OT-based reweighted distribution of the observational dataset.

We first define some notations used in the following discussion. Let
$\phi:{\mathcal{X}}\rightarrow{\mathcal{Z}}$ be a representation function with
inverse $\psi:{\mathcal{Z}}\rightarrow{\mathcal{X}}$, where ${\mathcal{Z}}$ is
the representation space. Let
$g:{\mathcal{Z}}\times{\mathcal{A}}\rightarrow{\mathcal{S}}$ and
$h:{\mathcal{Z}}\times{\mathcal{A}}\times{\mathcal{S}}\rightarrow{\mathcal{Y}}$
be the hypotheses that predict the short and long-term outcomes, respectively.
Note that the HRDC estimator aims to predict potential long-term outcome
$Y(a)$ for all possible treatments $a$ on each unit $\mathbf{x}$, which
requires that we could already predict the potential short-term outcomes
$\mathbf{S}(a)$ accurately. As a result, we define the combined loss for
learning these two hypotheses $g,h$ as
$\ell_{\phi,g,h}(\mathbf{x},a)=\ell^{s}_{\phi,g}(\mathbf{x},a)+\ell^{y}_{\phi,h}(\mathbf{x},a)$,
where
$\ell_{\phi,g}^{s}(\mathbf{x},a)=\int_{{\mathcal{S}}}L(\mathbf{s}(a),g(\phi(\mathbf{x}),a))P(\mathbf{s}(a)|\mathbf{x})d\mathbf{s}(a)$
and
$\ell_{\phi,h}^{y}(\mathbf{x},a)=\int_{{\mathcal{S}}\times{\mathcal{Y}}}L(y(a),h(\phi(\mathbf{x}),a,\mathbf{s}(a)))P(y(a),\mathbf{s}(a)|\mathbf{x})d\mathbf{s}(a)dy(a)$
($L(\cdot,\cdot)$ is the error function). Therefore, the ideal prediction
error over all treatment $a$ for unit $\mathbf{x}$ can be defined as
${\mathcal{E}}(\mathbf{x})={\mathbb{E}}_{a\sim
P(A)}[\ell_{\phi,g,h}(\mathbf{x},a)]$, and then the target of counterfactual
prediction is minimizing ${\mathcal{E}}_{cf}={\mathbb{E}}_{\mathbf{x}\sim
P(\mathbf{X})}[{\mathcal{E}}(\mathbf{x})]$.

However, ${\mathcal{E}}_{cf}$ is practically incomputable, because the actual
datasets we access lack the counterfactual outcomes. As a result, we can only
obtain the prediction error on the reweighted factual distribution as
${\mathcal{E}}_{f}^{\mathbf{w}_{o}}={\mathbb{E}}_{\mathbf{x},a\sim
P^{\mathbf{w}_{o}}(\mathbf{X},A)}[\ell_{\phi,g,h}(\mathbf{x},a)]$. Next,
following [46, 6], we derive a generalization bound to bridge a connection
between ${\mathcal{E}}_{cf}$ and ${\mathcal{E}}_{f}^{\mathbf{w}_{o}}$ in the
long-term HRDC scenario:

######  Theorem 4.

Assuming a family ${\mathcal{M}}$ of function
$m:{\mathcal{Z}}\times{\mathcal{A}}\rightarrow\mathbb{R}$, and there exists a
constant $B_{\phi}>0$ such that
$\frac{1}{B_{\Phi}}\cdot\ell_{\phi,g,h}(\mathbf{x},a)\in{\mathcal{M}}$, then
we have:

$\displaystyle{\mathcal{E}}_{cf}\leq{\mathcal{E}}_{f}^{\mathbf{w}_{o}}+IPM_{\mathcal{M}}(P(\mathbf{Z})P(A),~{}\mathbf{w}_{o}\cdot
P(\mathbf{Z},A)),$ (8)

where
$IPM_{{\mathcal{M}}}(p,q)=\sup_{m\in{\mathcal{M}}}\left|\int_{\mathcal{X}}g(x)(p(x)-q(x))dx\right|$
is integral probability metric for a chosen ${\mathcal{M}}$, $P(\mathbf{Z})$
and $P(\mathbf{Z},A)$ are distributions induced by the map $\phi$ from
$P(\mathbf{X})$ and $P(\mathbf{X},A)$.

The proof is given in Appendix A.5. Theorem 4 states that the counterfactual
error ${\mathcal{E}}_{cf}$ is upper-bounded by the factual error
${\mathcal{E}}_{f}^{\mathbf{w}_{o}}$ plus an IPM term that measures the
dependence between treatment $A$ and covariates $\mathbf{X}$. The theorem
inspires a classical balanced representation method based on deep learning
[46], which handles the observed confounding bias by learning the balanced
representation $\mathbf{Z}=\phi(\mathbf{X})$ that is independent of treatment
$A$, i.e., $\mathbf{Z}\perp\\!\\!\\!\perp A$.

Figure 2: Model architecture of the proposed LEARN.

##  5 Model Architecture

We have presented the theoretical results of our method, which handle the
unobserved and observed confounders by reweighting and representation
learning, respectively. In this section, we summarize the model architecture
of LEARN inspired by these theoretical results. As shown in Fig. 2, LEARN can
be divided into three modules: OT weight module, balanced representation
module and long-term estimation module.

OT Weight Module This module is designed based on Theorems 2 and 3, which aims
to estimate weights $\mathbf{w}_{o}$ for removing the unobserved confounding
bias. In order to embed the mini-batch OT problem $m\mbox{-}OT(\alpha,\beta)$
into the deep learning framework, following [59], we approximate it by solving
$OT(\alpha_{\mathcal{B}},\beta)$ over each mini-batch. As a result, we learn
the weights $\mathbf{w}_{o,\mathcal{B}}=\\{w_{o,1},w_{o,2},...,w_{o,b}\\}$ by
minimizing
$OT(P_{\mathcal{B}}(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))$
between the batch observational dataset and full experimental dataset at each
iteration. However, optimizing the OT problem directly usually induces a
sparse solution, which means only a limited number of observational units are
transported, suffering from low data efficiency [8, 54]. Motivated by [58], we
apply a negative entropy regularization on the marginal distributions
$\mathbf{w}_{o,\mathcal{B}}$, i.e.,
$\Omega(\mathbf{w}_{o,\mathcal{B}})=\sum_{j=1}^{b}w_{o,j}(\log w_{o,j}-1)$, to
encourage more observational samples to be transported, and also avoid the
heavy computation of solving the linear programming problem [16]. Finally, we
learning weights $\mathbf{w}_{o,\mathcal{B}}$ as follows:

$\displaystyle\min_{\mathbf{w}_{o,\mathcal{B}}}\min_{\gamma\in\Pi(\alpha_{\mathcal{B}},\beta)}\langle\gamma,C\rangle+\lambda_{e}\Omega(\gamma),$
$\displaystyle s.t.~{}$
$\displaystyle\Pi(\alpha_{\mathcal{B}},\beta)=\\{\gamma\in{\mathbb{R}}^{b\times
n_{e}}|\gamma\mathbf{1}_{n_{e}}=\mathbf{w}_{o,\mathcal{B}},\gamma^{T}\mathbf{1}_{b}=\bm{\mu},\gamma_{ij}\in\left[0,1\right]\\},$
(9)

where $\lambda_{e}$ is the hyperparameter for the entropic constraint. We
construct the cost matrix $\mathbf{C}$ reflecting the transport cost between
$P_{\mathcal{B}}(\mathbf{S},\mathbf{X},A|G=o)$ and
$P(\mathbf{S},\mathbf{X},A|G=e)$ as
$C_{ij}=\beta_{r}d(\bar{\mathbf{r}}_{i},\bar{\mathbf{r}}_{j})+\beta_{x}d(\mathbf{x}_{i},\mathbf{x}_{j})+\beta_{a}d(a_{i},a_{j})$,
where $\beta_{r},\beta_{x},\beta_{a}$ are the hyperparameters,
$d(\cdot,\cdot)$ is the euclidean distance, and the distance of short-term
outcomes is calculated from the mean embedding of RNN in the long-term
estimation module as
$\bar{\mathbf{r}}=\frac{1}{t_{0}}\sum_{t=1}^{t_{0}}\mathbf{r}_{t}$. As a
result, we need to pre-train the model firstly to ensure that the RNN’s
embedding is meaningful.

Following [58], we develop a projected mirror descent [40, 43] to solve the
problem (5), which firstly performs proximal gradient descent with the Bregman
divergence [5], and then obtains a feasible solution in the set
$\Pi(\alpha_{\mathcal{B}},\beta)$ by projection. The details can be found in
Appendix B.

Balanced Representation Module This module is designed based on Theorem 4,
which aims to correct the observed confounding bias at each iteration through
learning the balanced representation $\mathbf{Z}$ via minimizing
$IPM(P(\mathbf{Z})P(A),\mathbf{w}_{o,\mathcal{B}}P(\mathbf{Z},A))$. In
practice, we use the Wasserstein distance [53] as the implementation of the
IPM term. To calculate it, we simulate samples standing for the product of
marginal distributions $P(\mathbf{Z})P(A)$ by randomly permuting the observed
treatment $A$, and the original samples are drawn from the joint distribution
$P(\mathbf{Z},A)$.

Long-Term Estimation Module This module aims to estimate the potential long-
term outcome $Y(a)$ based on potential short-term outcomes $\mathbf{S}(a)$ for
different units, and we follow LTEE [12] to model the relationship between the
short and long-term outcomes. Specifically, the balanced representations
$\mathbf{Z}$ are firstly fed into the RNN $q(\cdot)$ to obtain the short-term
representations $\mathbf{R}_{t}$ at each timestep $t$, which is then used to
predict the corresponding short-term outcome $\hat{S}_{t}$ with a shared MLP
$g(\cdot)$. After that, we leverage the attention mechanism $f(\cdot)$ [4] to
construct the long-term representation $\mathbf{R}_{T}$ from
$\\{\mathbf{R}_{1},\mathbf{R}_{2},...,\mathbf{R}_{t_{0}}\\}$, and then predict
the long-term outcome $\hat{Y}_{T}$ through an MLP $h(\cdot)$.

Different from the original LTEE which only works with binary treatment, we
follow the varying coefficient networks [41] to extend this module to a
continuous treatment scenario. In particular, we use splines to model the
parameters $\theta(a)$ of the above four network structures
$q(\cdot),g(\cdot),f(\cdot),h(\cdot)$ as
$\theta(a)=\sum_{l=1}^{L}\alpha_{l}\phi_{l}(a)$, where $\phi_{l}$ is the
polynomial basis function and $\alpha_{l}$ is the coefficient. As a result,
the influence of the treatment information is not lost in high-dimensional
representations, which in practice has been shown to lead to better
performance.

The training loss at each iteration of the balanced representation module and
long-term estimation module is:

$\displaystyle\mathcal{L}_{\theta}=$
$\displaystyle\frac{1}{b}\sum_{i=1}^{b}w_{o,i}(Y_{i}-\hat{Y}_{i})^{2}+\frac{\lambda_{o}}{b}\sum_{i=1}^{b}\sum_{t=1}^{t_{0}}w_{o,i}(S_{i,t}-\hat{S}_{i,t})^{2}+\frac{1-\lambda_{o}}{n_{e}}\sum_{i=1}^{n_{e}}\sum_{t=1}^{t_{0}}(S_{i,t}-\hat{S}_{i,t})^{2}$
$\displaystyle+\lambda_{b}IPM(\\{\phi(\mathbf{X}_{i}),\tilde{A}_{i}\\}_{i=0}^{b},\\{w_{o,i}\phi(\mathbf{X}_{i}),w_{o,i}A_{i}\\}_{i=0}^{b}),$
(10)

where
$\theta=\\{\theta_{\phi},\theta_{q}(a),\theta_{g}(a),\theta_{f}(a),\theta_{h}(a)\\}$,
$\lambda_{b}$ means the strength of IPM, and $\tilde{A}_{i}$ is the samples
obtained by randomly permuting original treatment samples $A_{i}$. To make
full use of the datasets, in addition to the observational dataset, we also
train the model based on the prediction error of the short-term outcomes on
the experimental dataset, and use $\lambda_{o}$ to control its strength. We
also provide more detail and the pseudocode of LEARN in Appendix D.

##  6 Experiments

###  6.1 Experimental Setup

Dataset Generation Since the true HDRC is unavailable for real-world datasets,
following existing literature [45, 7, 41, 56, 34], we use five synthetic
datasets with varying $\beta_{U}\in\\{1,1.25,1.5,1.75,2\\}$ that denotes
different strengths of unobserved confounding bias, and two semi-synthetic
datasets, News [45] and TCGA [57], to demonstrate the effectiveness of our
model. The details of the dataset generation can be found in Appendix E.1.

Baselines and Metrics We first compare our model LEARN with the following
baselines designed for continuous treatment, including DRNet [45], VCNet [41],
SCIGAN [7], ADMIT [56] and ACFR [34]. Besides, since there is a lack of work
on estimating long-term HDRC, we modify LTEE [12] into a baseline suitable for
continuous treatment by introducing the varying coefficient structure, denoted
as LTEE(VCNet). Moreover, to verify the effectiveness of the OT weight module
and IPM term, we implement two variants of our model: LEARN(None) is the
variant without any module to correct confounding bias, and LEARN(IPM) is the
variant with only the IPM term to correct observed confounding bias without
considering the unobserved one. For metrics, we first split the observational
dataset with a $6/2/2$ ratio into training, validation and test sets, and we
report the mean integrated square error (MISE) and the deviations (std) of
HDRC estimation on the test set by performing 10 replications:
$\mbox{MISE}=\frac{1}{s}\sum_{i=1}^{s}\int_{a_{1}}^{a_{2}}(Y(a)-\hat{Y}(a))^{2}da$,
where $s$ is the test sample size and $[a_{1},a_{2}]$ is the sampling interval
of treatment values.

Table 1: Results of HDRC estimation on five synthetic and two semi-synthetic
datasets. The MISE is reported as $\mbox{Mean}\pm\mbox{Std}$. The best-
performing method is bolded.

Dataset | Synthetic | News | TCGA
---|---|---|---
$\beta_{U}=1$ | $\beta_{U}=1.25$ | $\beta_{U}=1.5$ | $\beta_{U}=1.75$ | $\beta_{U}=2$
DRNet | $66.77\pm 13.26$ | $69.23\pm 13.64$ | $72.92\pm 13.73$ | $78.43\pm 14.68$ | $86.67\pm 18.61$ | $26.47\pm 10.07$ | $0.788\pm 0.781$
SCIGAN | $75.92\pm 14.17$ | $78.29\pm 12.85$ | $82.03\pm 14.13$ | $87.82\pm 14.88$ | $94.64\pm 15.61$ | $26.61\pm 10.09$ | $1.747\pm 1.823$
VCNet | $63.31\pm 12.77$ | $65.78\pm 12.64$ | $69.55\pm 12.39$ | $74.19\pm 12.62$ | $79.86\pm 13.00$ | $9.42\pm 3.13$ | $0.187\pm 0.146$
ADMIT | $62.91\pm 12.13$ | $66.65\pm 13.69$ | $70.24\pm 13.54$ | $74.59\pm 12.92$ | $81.15\pm 14.18$ | $16.58\pm 8.72$ | $0.143\pm 0.081$
ACRF | $64.64\pm 8.33$ | $67.38\pm 8.14$ | $71.30\pm 8.28$ | $76.17\pm 9.52$ | $82.97\pm 11.29$ | $22.81\pm 9.39$ | $0.465\pm 0.379$
LTEE(VCNet) | $61.79\pm 12.48$ | $63.98\pm 12.90$ | $67.55\pm 12.03$ | $73.34\pm 12.67$ | $78.79\pm 13.33$ | $10.11\pm 2.76$ | $0.137\pm 0.097$
LEARN(None) | $62.49\pm 12.24$ | $64.27\pm 12.69$ | $68.08\pm 12.09$ | $74.08\pm 12.95$ | $79.99\pm 12.82$ | $8.86\pm 3.18$ | $0.145\pm 0.111$
LEARN(IPM) | $61.75\pm 12.64$ | $63.90\pm 12.13$ | $67.41\pm 12.01$ | $73.05\pm 12.82$ | $78.37\pm 13.00$ | $8.76\pm 3.14$ | $0.141\pm 0.106$
LEARN | $60.57\pm 11.98$ | $62.74\pm 12.24$ |  $66.02\pm 11.59$ | $71.00\pm 12.17$ | $76.76\pm 12.90$ | $8.31\pm 2.90$ | $0.131\pm 0.101$

###  6.2 Result and Analysis

Overall Performance In Table 1, we compare our model with baselines on both
synthetic and semi-synthetic datasets, and the main observations are as
follows. Firstly, compared with the baselines only designed for continuous
treatment, LEARN obtains the highest accuracy in estimating the HRDC in all
datasets. The accuracy improvement comes from the designed long-term
estimation module, employing the attention mechanism to capture the
relationship between the short-term and long-term outcomes, and this is also
the reason why LTEE(VCNet) performs relatively well. Secondly, we observe that
LEARN(IPM) consistently outperforms LEARN(None), indicating that learning
balanced representations through minimizing the IPM term indeed reduces the
observed confounding bias, resulting in smaller MISE. Finally, LEARN achieves
superior performance compared with LEARN(IPM), and the improvement increases
with the increase of $\beta_{U}$, which reflects the necessity of considering
the unobserved confounders. The improvement verifies the effectiveness of the
reweighting based on our designed OT weight module.

Performance of Confounding Bias Correction In Table 2, we evaluate the
effectiveness of our design modules in correcting confounding bias using the
synthetic dataset ($\beta_{U}=1$), including the observed and unobserved ones.
Regarding observed confounding bias, it arises from the correlation between
$\mathbf{X}$ and $A$, and we verify the effectiveness of a balanced
representation module in correcting such a bias. We use the Hilbert-Schmidt
Independence Criterion (HSIC[24]) to measure the independence between the
learned $\mathbf{Z}=\phi(\mathbf{X})$ and $A$. The result reveals an $86.11\%$
reduction in the correlation between $\mathbf{X}$ and $A$, indicating the
remarkable decrease in the observed confounding bias, which is also consistent
with the comparison results between LEARN(None) and LEARN(IPM) in Table 1.

Table 2: Results of the model performance in reducing observed and unobserved
confounding bias, using HSIC as a metric of correlation.

| Observed | Unobserved
---|---|---
| $\mbox{HSIC}(\mathbf{X},A)$ | $\mbox{HSCONIC}(\mathbf{S}(a),A|\mathbf{X})$ | $\mbox{HSCONIC}(Y(a),A|\mathbf{X})$
Reduce ratio | $86.11\%$ | $61.55\%$ | $66.17\%$

Regarding unobserved confounders, our designed OT weight module handles them
by reweighting. The effectiveness of this module is evaluated using the
Hilbert–Schmidt Conditional Independence Criterion (HSCONIC, [20, 47]), which
measures the degree to which reweighting makes the unconfoundedness in the
observational dataset satisfied, i.e.,
$A\perp\\!\\!\\!\perp\\{Y(a),\mathbf{S}(a)\\}|\mathbf{X},G=o$. The result
indicates a substantial decrease in the correlations between $\mathbf{S}(a)$
and $A$, as well as between $Y(a)$ and $A$ given $\mathbf{X}$ , by $61.55\%$
and $66.17\%$ respectively, suggesting an enhanced satisfaction of the
unconfoundedness assumption after reweighting, consistent with the comparison
results between LEARN(IPM) and LEARN in Table 1.

Impact of Sample Size on Experimental Dataset Since our proposed LEARN relies
on the data combination with an experimental dataset, we exploit the impact of
the experimental dataset’s sample size on MISE using the synthetic dataset
($\beta_{U}=1$) and report the results in Table 3. As expected, the MISE
decreases as the sample size increases, attributable to the larger
experimental dataset providing a more reliable distribution for aligning the
observational dataset, which allows us to learn more effective weights for
eliminating unobserved confounding bias.

Table 3: The impact of experimental dataset’s sample size on the MISE of HDRC
estimation.

Number | 100 | 250 | 500 | 1000 | 2000
---|---|---|---|---|---
MISE | $61.24\pm 11.74$ | $60.95\pm 11.53$ | $60.57\pm 11.98$ | $59.46\pm 11.29$ | $59.02\pm 11.39$

Hyperparameter Sensitivity Study We investigate the role of three
hyperparameters in our model: $\lambda_{o},\lambda_{b},\lambda_{e}$. The
details can be found in Appendix E.2.

##  7 Conclusion

In this paper, we provide a practical solution to estimate the long-term HDRC
with unobserved confounders on observational data via data combination, which
is not well-studied in existing work. Specifically, to remove unobserved
confounders and make long-term HDRC identifiable, we first propose a novel
weighting schema aligning the conditional distribution of short-term outcomes
between observational and experimental groups. We solve the reweighting
problem under an optimal transport framework, showing that the conditional
distribution discrepancy can be bounded by mini-batch joint distribution
discrepancy, which is computationally efficient. Further, to handle the
observed confounders, we derive a generalization bound on the counterfactual
outcome regression error based on the OT-induced reweighted distribution.
Building upon the above theoretical results, we develop our model LEARN to
estimate long-term HDRC accurately. Extensive experimental results verify the
correctness of our theory and the effectiveness of LEARN.

## Acknowledgements

This research was supported in part by National Key R&D Program of China
(2021ZD0111501), National Science Fund for Excellent Young Scholars
(62122022), Natural Science Foundation of China (62206064, 62206061),
Guangdong Basic and Applied Basic Research Foundation (2024A1515011901),
Guangzhou Basic and Applied Basic Research Foundation (2023A04J1700), and CCF-
DiDi GAIA Collaborative Research Funds (CCF-DiDi GAIA 202311).

## References

  * Arjovsky et al. [2017] M. Arjovsky, S. Chintala, and L. Bottou.  Wasserstein generative adversarial networks.  In _International conference on machine learning_ , pages 214–223. PMLR, 2017.
  * Athey et al. [2019] S. Athey, R. Chetty, G. W. Imbens, and H. Kang.  The surrogate index: Combining short-term proxies to estimate long-term treatment effects more rapidly and precisely.  Technical report, National Bureau of Economic Research, 2019.
  * Athey et al. [2020] S. Athey, R. Chetty, and G. Imbens.  Combining experimental and observational data to estimate treatment effects on long term outcomes.  _arXiv preprint arXiv:2006.09676_ , 2020.
  * Bahdanau et al. [2014] D. Bahdanau, K. Cho, and Y. Bengio.  Neural machine translation by jointly learning to align and translate.  _arXiv preprint arXiv:1409.0473_ , 2014.
  * Banerjee et al. [2005] A. Banerjee, S. Merugu, I. S. Dhillon, J. Ghosh, and J. Lafferty.  Clustering with bregman divergences.  _Journal of machine learning research_ , 6(10), 2005.
  * Bellot et al. [2022] A. Bellot, A. Dhir, and G. Prando.  Generalization bounds and algorithms for estimating conditional average treatment effect of dosage.  _arXiv preprint arXiv:2205.14692_ , 2022.
  * Bica et al. [2020] I. Bica, J. Jordon, and M. van der Schaar.  Estimating the effects of continuous-valued interventions using generative adversarial networks.  _Advances in Neural Information Processing Systems_ , 33:16434–16445, 2020.
  * Blondel et al. [2018] M. Blondel, V. Seguy, and A. Rolet.  Smooth and sparse optimal transport.  In _International conference on artificial intelligence and statistics_ , pages 880–889. PMLR, 2018.
  * Cai et al. [2024] R. Cai, W. Chen, Z. Yang, S. Wan, C. Zheng, X. Yang, and J. Guo.  Long-term causal effects estimation via latent surrogates representation learning.  _Neural Networks_ , page 106336, 2024.
  * Chang et al. [2022] W. Chang, Y. Shi, H. Tuan, and J. Wang.  Unified optimal transport framework for universal domain adaptation.  _Advances in Neural Information Processing Systems_ , 35:29512–29524, 2022.
  * Chen and Ritzwoller [2023] J. Chen and D. M. Ritzwoller.  Semiparametric estimation of long-term treatment effects.  _Journal of Econometrics_ , 237(2):105545, 2023\.
  * Cheng et al. [2021] L. Cheng, R. Guo, and H. Liu.  Long-term effect estimation with surrogate representation.  In _Proceedings of the 14th ACM International Conference on Web Search and Data Mining_ , pages 274–282, 2021.
  * Chetty et al. [2011] R. Chetty, J. N. Friedman, N. Hilger, E. Saez, D. W. Schanzenbach, and D. Yagan.  How does your kindergarten classroom affect your earnings? evidence from project star.  _The Quarterly journal of economics_ , 126(4):1593–1660, 2011.
  * Courty et al. [2016] N. Courty, R. Flamary, D. Tuia, and A. Rakotomamonjy.  Optimal transport for domain adaptation.  _IEEE transactions on pattern analysis and machine intelligence_ , 39(9):1853–1865, 2016.
  * Courty et al. [2017] N. Courty, R. Flamary, A. Habrard, and A. Rakotomamonjy.  Joint distribution optimal transportation for domain adaptation.  _Advances in neural information processing systems_ , 30, 2017.
  * Cuturi [2013] M. Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  In _Annual Conference on Neural Information Processing Systems_ , pages 2292–2300, 2013.
  * Damodaran et al. [2018] B. B. Damodaran, B. Kellenberger, R. Flamary, D. Tuia, and N. Courty.  Deepjdot: Deep joint distribution optimal transport for unsupervised domain adaptation.  In _Proceedings of the European conference on computer vision (ECCV)_ , pages 447–463, 2018.
  * Dunipace [2021] E. Dunipace.  Optimal transport weights for causal inference.  _arXiv preprint arXiv:2109.01991_ , 2021.
  * Fleming et al. [1994] T. R. Fleming, R. L. Prentice, M. S. Pepe, and D. Glidden.  Surrogate and auxiliary endpoints in clinical trials, with potential applications in cancer and aids research.  _Statistics in medicine_ , 13(9):955–968, 1994\.
  * Fukumizu et al. [2007] K. Fukumizu, A. Gretton, X. Sun, and B. Schölkopf.  Kernel measures of conditional dependence.  _Advances in neural information processing systems_ , 20, 2007.
  * Gao et al. [2024] J. Gao, H. Zhao, Z. Li, and D. Guo.  Enhancing minority classes by mixing: An adaptative optimal transport approach for long-tailed classification.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Genevay et al. [2018] A. Genevay, G. Peyré, and M. Cuturi.  Learning generative models with sinkhorn divergences.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1608–1617. PMLR, 2018.
  * Ghassami et al. [2022] A. Ghassami, A. Yang, D. Richardson, I. Shpitser, and E. T. Tchetgen.  Combining experimental and observational data for identification and estimation of long-term causal effects.  _arXiv preprint arXiv:2201.10743_ , 2022.
  * Gretton et al. [2005] A. Gretton, O. Bousquet, A. Smola, and B. Schölkopf.  Measuring statistical dependence with hilbert-schmidt norms.  In _International conference on algorithmic learning theory_ , pages 63–77. Springer, 2005.
  * Gunsilius and Xu [2021] F. Gunsilius and Y. Xu.  Matching for causal effects via multimarginal unbalanced optimal transport.  _arXiv preprint arXiv:2112.04398_ , 2021.
  * Guo et al. [2022] D. Guo, Z. Li, H. Zhao, M. Zhou, H. Zha, et al.  Learning to re-weight examples with optimal transport for imbalanced classification.  _Advances in Neural Information Processing Systems_ , 35:25517–25530, 2022.
  * Hohnhold et al. [2015] H. Hohnhold, D. O’Brien, and D. Tang.  Focusing on the long-term: It’s good for users and business.  In _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 1849–1858, 2015.
  * [28] W. Hu, X. Zhou, and P. Wu.  Identification and estimation of treatment effects on long-term outcomes in clinical trials with external observational data.  _Statistica Sinica_.  doi: 10.5705/ss.202023.0006.
  * Huynh and Phung [2021] V. Huynh and D. Phung.  Optimal transport for deep generative models: State of the art and research challenges.  In _International Joint Conference on Artificial Intelligence 2021_ , pages 4450–4457. Association for the Advancement of Artificial Intelligence (AAAI), 2021.
  * Imbens et al. [2022] G. Imbens, N. Kallus, X. Mao, and Y. Wang.  Long-term causal inference under persistent confounding via data combination.  _arXiv preprint arXiv:2202.07234_ , 2022.
  * Imbens [2000] G. W. Imbens.  The role of the propensity score in estimating dose-response functions.  _Biometrika_ , 87(3):706–710, 2000.
  * Kallus and Mao [2020] N. Kallus and X. Mao.  On the role of surrogates in the efficient estimation of treatment effects with limited outcome data.  _arXiv preprint arXiv:2003.12408_ , 2020.
  * Kantorovich [2006] L. V. Kantorovich.  On the translocation of masses.  _Journal of mathematical sciences_ , 133(4):1381–1382, 2006.
  * Kazemi and Ester [2024] A. Kazemi and M. Ester.  Adversarially balanced representation for continuous treatment effect estimation.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 38, pages 13085–13093, 2024.
  * Kennedy et al. [2017] E. H. Kennedy, Z. Ma, M. D. McHugh, and D. S. Small.  Non-parametric methods for doubly robust estimation of continuous treatment effects.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 79(4):1229–1245, 2017.
  * Kingma and Ba [2014] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Korotin et al. [2019] A. Korotin, V. Egiazarian, A. Asadulaev, A. Safin, and E. Burnaev.  Wasserstein-2 generative networks.  _arXiv preprint arXiv:1909.13082_ , 2019.
  * Li et al. [2023] M. Li, D. Wang, X. Liu, Z. Zeng, R. Lu, B. Chen, and M. Zhou.  Patchct: Aligning patch set and label set with conditional transport for multi-label image classification.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 15348–15358, 2023.
  * Li et al. [2021] Q. Li, Z. Wang, S. Liu, G. Li, and G. Xu.  Causal optimal transport for treatment effect estimation.  _IEEE transactions on neural networks and learning systems_ , 34(8):4083–4095, 2021.
  * Nemirovskij and Yudin [1983] A. S. Nemirovskij and D. B. Yudin.  Problem complexity and method efficiency in optimization.  1983\.
  * Nie et al. [2021] L. Nie, M. Ye, qiang liu, and D. Nicolae.  Varying coefficient neural network with functional targeted regularization for estimating continuous treatment effects.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=RmB-88r9dL.
  * Owen [2013] A. B. Owen.  Monte carlo theory, methods and examples, 2013.
  * Raskutti and Mukherjee [2015] G. Raskutti and S. Mukherjee.  The information geometry of mirror descent.  _IEEE Transactions on Information Theory_ , 61(3):1451–1457, 2015.
  * Rubin [1974] D. B. Rubin.  Estimating causal effects of treatments in randomized and nonrandomized studies.  _Journal of educational Psychology_ , 66(5):688, 1974.
  * Schwab et al. [2020] P. Schwab, L. Linhardt, S. Bauer, J. M. Buhmann, and W. Karlen.  Learning counterfactual representations for estimating individual dose-response curves.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pages 5612–5619, 2020.
  * Shalit et al. [2017] U. Shalit, F. D. Johansson, and D. Sontag.  Estimating individual treatment effect: generalization bounds and algorithms.  In _International conference on machine learning_ , pages 3076–3085. PMLR, 2017.
  * Sun et al. [2007] X. Sun, D. Janzing, B. Schölkopf, and K. Fukumizu.  A kernel-based causal learning algorithm.  In _Proceedings of the 24th international conference on Machine learning_ , pages 855–862, 2007.
  * Tian et al. [2023] L. Tian, J. Feng, X. Chai, W. Chen, L. Wang, X. Liu, and B. Chen.  Prototypes-oriented transductive few-shot learning with conditional transport.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 16317–16326, 2023.
  * Tübbicke [2022] S. Tübbicke.  Entropy balancing for continuous treatments.  _Journal of Econometric Methods_ , 11(1):71–89, 2022.
  * Van Goffrier et al. [2023] G. Van Goffrier, L. Maystre, and C. M. Gilligan-Lee.  Estimating long-term causal effects from short-term experiments and long-term observational data with unobserved confounding.  In _Conference on Causal Learning and Reasoning_ , pages 791–813. PMLR, 2023.
  * Vegetabile et al. [2021] B. G. Vegetabile, B. A. Griffin, D. L. Coffman, M. Cefalu, M. W. Robbins, and D. F. McCaffrey.  Nonparametric estimation of population average dose-response curves using entropy balancing weights for continuous exposures.  _Health Services and Outcomes Research Methodology_ , 21:69–110, 2021.
  * Villani [2021] C. Villani.  _Topics in optimal transportation_ , volume 58.  American Mathematical Soc., 2021.
  * Villani et al. [2009] C. Villani et al.  _Optimal transport: old and new_ , volume 338.  Springer, 2009.
  * Vincent-Cuaz et al. [2022] C. Vincent-Cuaz, R. Flamary, M. Corneli, T. Vayer, and N. Courty.  Semi-relaxed gromov wasserstein divergence with applications on graphs.  In _International Conference on Learning Representations_ , 2022.
  * Wang et al. [2024] H. Wang, J. Fan, Z. Chen, H. Li, W. Liu, T. Liu, Q. Dai, Y. Wang, Z. Dong, and R. Tang.  Optimal transport for treatment effect estimation.  _Advances in Neural Information Processing Systems_ , 36, 2024.
  * Wang et al. [2022] X. Wang, S. Lyu, X. Wu, T. Wu, and H. Chen.  Generalization bounds for estimating causal effects of continuous treatments.  _Advances in Neural Information Processing Systems_ , 35:8605–8617, 2022.
  * Weinstein et al. [2013] J. N. Weinstein, E. A. Collisson, G. B. Mills, K. R. Shaw, B. A. Ozenberger, K. Ellrott, I. Shmulevich, C. Sander, and J. M. Stuart.  The cancer genome atlas pan-cancer analysis project.  _Nature genetics_ , 45(10):1113–1120, 2013.
  * Yan et al. [2024] Y. Yan, Z. Yang, W. Chen, R. Cai, Z. Hao, and M. K.-P. Ng.  Exploiting geometry for treatment effect estimation via optimal transport.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 38, pages 16290–16298, 2024.
  * Yang et al. [2023] Y. Yang, X. Gu, and J. Sun.  Prototypical partial optimal transport for universal domain adaptation.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 37, pages 10852–10860, 2023.
  * Zeng et al. [2024] Z. Zeng, D. Arbour, A. Feller, R. Addanki, R. Rossi, R. Sinha, and E. H. Kennedy.  Continuous treatment effects with surrogate outcomes.  _arXiv preprint arXiv:2402.00168_ , 2024.
  * Zhu et al. [2015] Y. Zhu, D. L. Coffman, and D. Ghosh.  A boosting algorithm for estimating generalized propensity scores with continuous treatments.  _Journal of causal inference_ , 3(1):25–40, 2015\.

##  Appendix A Proof

###  A.1 Proof of Proposition 1

######  Proposition 1.

Under Assumptions 1, 2,3 4, and 5, given a set of weights
$\mathbf{w}=\\{\mathbf{w}_{o},~{}\bm{\mu}\\}$ consisting of the learnable
weights $\mathbf{w}_{o}$ for observational units and uniform weights
$\bm{\mu}$ for experimental units, which makes
$\mathbf{w}P(\mathbf{S},G|\mathbf{X},A)=P^{\mathbf{w}}(\mathbf{S},G|\mathbf{X},A)=P^{\mathbf{w}}(\mathbf{S}|\mathbf{X},A)P^{\mathbf{w}}(G|\mathbf{X},A)$,
i.e., $\mathbf{S}\perp\\!\\!\\!\perp G|\mathbf{X},A$, then
$\mathbb{E}[\mathbf{S}(a)|\mathbf{X},A=a,G=o]=\mathbb{E}[\mathbf{S}(a)|\mathbf{X},G=o]$
holds on the weighted observational population.

###### Proof.

First, based on Assumption 1, we rewrite
$P^{\mathbf{w}}(\mathbf{S},G|\mathbf{X},A)=P^{\mathbf{w}}(\mathbf{S}|\mathbf{X},A)P^{\mathbf{w}}(G|X,A)$
as
$P^{\mathbf{w}}(\mathbf{S}(a),G|\mathbf{X},A=a)=P^{\mathbf{w}}(\mathbf{S}(a)|\mathbf{X},A=a)P^{\mathbf{w}}(G|\mathbf{X},A=a)$.
Based on chain rule, we have

$P^{\mathbf{w}}(\mathbf{S}(a),G|\mathbf{X},A=a)=P^{\mathbf{w}}(\mathbf{S}(a)|G,\mathbf{X},A=a)P^{\mathbf{w}}(G|\mathbf{X},A=a),$
(1)

then we have

$P^{\mathbf{w}}(\mathbf{S}(a)|\mathbf{X},A=a)=P^{\mathbf{w}}(\mathbf{S}(a)|G,\mathbf{X},A=a).$
(2)

Note that, based on Eq. (2), we have

$P^{\mathbf{w}}(\mathbf{S}(a)|\mathbf{X},A=a)=P^{\mathbf{w}}(\mathbf{S}(a)|G,\mathbf{X},A=a)=P^{\mathbf{w}}(\mathbf{S}(a)|G=e,\mathbf{X},A=a),$
(3)

and thus we have

$P^{\mathbf{w}}(\mathbf{S}(a),G|\mathbf{X},A=a)=P^{\mathbf{w}}(\mathbf{S}(a)|G=e,\mathbf{X},A=a)P^{\mathbf{w}}(G|\mathbf{X},A=a).$
(4)

Based on Assumption 3, we have

$P^{\mathbf{w}}(\mathbf{S}(a),G|\mathbf{X},A=a)=P^{\mathbf{w}}(\mathbf{S}(a)|G=e,\mathbf{X})P^{\mathbf{w}}(G|\mathbf{X},A=a).$
(5)

Based on Assumption 4, we have

$P^{\mathbf{w}}(\mathbf{S}(a),G|\mathbf{X},A=a)=P^{\mathbf{w}}(\mathbf{S}(a)|G,\mathbf{X})P^{\mathbf{w}}(G|\mathbf{X},A=a).$
(6)

Hence, based on Eq. (1) and (6), we conclude

$P^{\mathbf{w}}(\mathbf{S}(a)|A=a,G,\mathbf{X})=P^{\mathbf{w}}(\mathbf{S}(a)|G,\mathbf{X}),$
(7)

and thus
$\mathbb{E}[\mathbf{S}(a)|A=a,G=o,\mathbf{X}]=\mathbb{E}[\mathbf{S}(a)|G=o,\mathbf{X}]$
hold.

∎

###  A.2 Proof of Theorem 1

######  Theorem 1.

Suppose Assumptions in Proposition 1 hold, we have the unconfoundedness
$\mathbb{E}[Y(a)|\mathbf{X},A=a,G=o]=\mathbb{E}[Y(a)|\mathbf{X},G=o]$, then
the long-term HDRC can be identified.

###### Proof.

Firstly, we show that
$\mathbb{E}[Y(a)|\mathbf{X},A=a,G=o]=\mathbb{E}[Y(a)|\mathbf{X},G=o]$ holds
based on Proposition 1. Specifically,

$\displaystyle P^{\mathbf{w}}(Y(a)|A=a,G=o,\mathbf{X})$ $\displaystyle=\int
P^{\mathbf{w}}(Y(a),\mathbf{S}(a)|A=a,G=o,\mathbf{X})d\mathbf{S}(a)$ (8)
$\displaystyle=\int
P^{\mathbf{w}}(Y(a)|\mathbf{S}(a),A=a,G=o,\mathbf{X})P^{\mathbf{w}}(\mathbf{S}(a)|A=a,G=o,\mathbf{X})d\mathbf{S}(a)$
$\displaystyle=\int
P^{\mathbf{w}}(Y(a)|\mathbf{S}(a),G=o,\mathbf{X})P^{\mathbf{w}}(\mathbf{S}(a)|G=o,\mathbf{X})d\mathbf{S}(a)$
$\displaystyle=\int
P^{\mathbf{w}}(Y(a),\mathbf{S}(a)|G=o,\mathbf{X})d\mathbf{S}(a)$
$\displaystyle=P^{\mathbf{w}}(Y(a)|G=o,\mathbf{X}),$

where the second equality is based on the chain rule, the third equality is
based on Assumption 5 and Eq. (7). Based on the above equality, we can
directly conclude

$\displaystyle\mathbb{E}[Y(a)|A=a,G=o,\mathbf{X}]=\mathbb{E}[Y(a)|G=o,\mathbf{X}].$
(9)

Then HDRC is identified:

$\displaystyle\mathbb{E}[Y(a)|X=x]$ $\displaystyle=\mathbb{E}[Y(a)|X=x,G=o]$
(10) $\displaystyle=\mathbb{E}[Y(a)|X=x,A=a,G=o]$
$\displaystyle=\mathbb{E}[Y|X=x,A=a,G=o],$

where the first equality is based on Assumption 4, the second equality is
based on the result of Eq. (9), and the third equality is based on Assumption
1. ∎

###  A.3 Proof of Theorem 2

We firstly define the empirical observational distribution
$P(\mathbf{S},\mathbf{X},A|G=o)$ as
$\alpha_{\mathbf{w}_{o}}=\sum_{i=1}^{n_{o}}w_{o,i}\delta(\mathbf{s}_{i},\mathbf{x}_{i},a_{i})$,
where $\mathbf{w}_{o}$ is the to-be-learned weights, and we define the
empirical experimental distribution $P(\mathbf{S},\mathbf{X},A|G=e)$ as
$\beta=\sum_{i=1}^{n_{e}}\frac{1}{n_{e}}\delta(\mathbf{s}_{i},\mathbf{x}_{i},a_{i})$.
We also have the following lemma, which states the convergence of the
importance sampling weights:

######  Lemma 1.

Under the positivity assumption, let the importance sampling weights as
$\breve{w}^{\star}_{o,i}=\frac{dP(\mathbf{s}_{i},\mathbf{x}_{i},a_{i}|G=e)}{dP(\mathbf{s}_{i},\mathbf{x}_{i},a_{i}|G=o)}$
and the corresponding self-normalized weights as
$w^{\star}_{o,i}=\frac{\breve{w}^{\star}_{o,i}}{\sum_{i}\breve{w}^{\star}_{o,i}}$.
Then we have $P(\lim_{n_{o}\rightarrow\infty}\alpha_{w^{\star}_{o}}=\beta)=1$.

Based on Theorem 9.2 in [42], Lemma 1 is easy to establish by viewing
$\delta(\mathbf{s}_{i},\mathbf{x}_{i},a_{i})$ as the function of interest
$f(\cdot)$ in the importance sampling.

Now we are ready to prove Theorem 2.

######  Theorem 2.

Assuming the cost matrix of the joint distribution
$P(\mathbf{S},\mathbf{X},A)$ is separable, i.e.,
$C(\mathbf{s},\mathbf{x},a;\tilde{\mathbf{s}},\tilde{\mathbf{x}},\tilde{a})=C(\mathbf{s};\tilde{\mathbf{s}})+C(\mathbf{x};\tilde{\mathbf{x}})+C(a;\tilde{a})$,
we have:

$\displaystyle OT^{con}_{\mathbf{x},a}\leq
OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e)),$

In addition, under the positivity assumption, we have
$OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))\to 0$ as
$n_{o}\to\infty$, leading to $OT^{con}_{\mathbf{x},a}\to 0$.

###### Proof.

With the separable assumption of the cost matrix, we abbreviate the cost
matrix about the joint distribution $P(\mathbf{S},\mathbf{X},A)$ as
$C_{\mathbf{s},\mathbf{x},a}=C_{\mathbf{s}}+C_{\mathbf{x}}+C_{a}$. Then for
any
$\Gamma_{\mathbf{s},\mathbf{x},a,\tilde{\mathbf{s}},\tilde{\mathbf{x}},\tilde{a}}\in\Pi(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))$,
we have:

$\displaystyle
OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))$
$\displaystyle=\min_{\Gamma}~{}\langle\Gamma_{\mathbf{s},\mathbf{x},a,\tilde{\mathbf{s}},\tilde{\mathbf{x}},\tilde{a}},~{}C_{\mathbf{s},\mathbf{x},a}\rangle$
$\displaystyle=\min_{\Gamma}~{}\langle\Gamma_{\mathbf{s},\mathbf{x},a,\tilde{\mathbf{s}},\tilde{\mathbf{x}},\tilde{a}},~{}C_{\mathbf{s}}+C_{\mathbf{x}}+C_{a}\rangle$
$\displaystyle\geqslant\min_{\Gamma}~{}\langle\Gamma_{\mathbf{s},\mathbf{x},a,\tilde{\mathbf{s}},\tilde{\mathbf{x}},\tilde{a}},~{}C_{\mathbf{s}}\rangle$
$\displaystyle\geqslant\min_{\Gamma}~{}\langle\Gamma_{\mathbf{s},\mathbf{x},a,\tilde{\mathbf{s}},\tilde{\mathbf{x}},\tilde{a}}\mathbb{I}(\mathbf{x}=\tilde{\mathbf{x}},a=\tilde{a}),~{}C_{\mathbf{s}}\rangle$
$\displaystyle=\min_{\Gamma}~{}\sum_{\mathbf{x},a,\tilde{\mathbf{x}},\tilde{a}}\mathbb{I}(\mathbf{x}=\tilde{\mathbf{x}},a=\tilde{a})\langle\Gamma_{\mathbf{s}|x,a,\tilde{\mathbf{s}}|\tilde{\mathbf{x}},\tilde{a}},~{}C_{\mathbf{s}}\rangle$
$\displaystyle=\min_{\Gamma}~{}\sum_{\mathbf{x},a}\langle\Gamma_{\mathbf{s}|\mathbf{x},a,\tilde{\mathbf{s}}|\mathbf{x},a},~{}C_{\mathbf{s}}\rangle$
$\displaystyle=OT^{con}_{\mathbf{x},a}.$

Next, we establish the convergence property. Recall that the definition of
$\hat{\mathbf{w}}_{o}=\min_{\mathbf{w}_{o}}OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))=\min_{\mathbf{w}_{o}}OT(\alpha_{\mathbf{w}_{o}},\beta)$,
we have $OT(\alpha_{\hat{\mathbf{w}}_{o}},\beta)\leq
OT(\alpha_{\mathbf{w}^{\star}_{o}},\beta)$. Based on Lemma 1, when
$n_{o}\rightarrow\infty$, we could obtain:

$\displaystyle\alpha_{\mathbf{w}^{\star}_{o}}\rightarrow\beta\Longrightarrow$
$\displaystyle~{}~{}~{}OT(\alpha_{\mathbf{w}^{\star}_{o}},\beta)\rightarrow 0$
$\displaystyle\Longrightarrow$
$\displaystyle~{}~{}~{}OT(\alpha_{\hat{\mathbf{w}}_{o}},\beta)\rightarrow 0,$

which means that with the estimated weights $\hat{\mathbf{w}}_{o}$, we have
$OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))\rightarrow
0$. Therefore, we also have $OT^{con}_{\mathbf{x},a}\rightarrow 0$, since it
is the lower bound of
$OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))$. ∎

###  A.4 Proof of Theorem 3

######  Theorem 3.

Let $\gamma_{i}$ be the optimal transport probability matrix of the $i$-th
batch OT problem of
$m\mbox{-}OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))$,
i.e.,
$OT(P_{\mathcal{B}_{i}}(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))$.
We extend $\gamma_{i}$ to a $n_{o}\times n_{e}$ matrix $\Gamma_{i}$ that pads
zero entries to the row whose index does not belong to $\mathcal{B}_{i}$, then
we have:

$\displaystyle\frac{1}{k}\sum_{i=1}^{k}\Gamma_{i}\in\Pi(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e)),$
(11)

and

$\displaystyle
OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))\leq
m\mbox{-}OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e)).$
(12)

###### Proof.

Assuming the probability measure of $P(\mathbf{S},\mathbf{X},A|G=o)$ and
$P(\mathbf{S},\mathbf{X},A|G=e)$ are $\alpha,\beta$, respectively. Then
according to the definition, the proof of Eq. (11) is equivalent to prove:

$\displaystyle\left(\frac{1}{k}\sum_{i=1}^{k}\Gamma_{i}\right)\mathbf{1}_{n_{e}}=\alpha,~{}~{}~{}\left(\frac{1}{k}\sum_{i=1}^{k}\Gamma_{i}\right)^{T}\mathbf{1}_{n_{o}}=\beta.$

Note that
$\gamma_{i}\in\Pi(P_{\mathcal{B}_{i}}(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))$
satisfies:

$\displaystyle\gamma_{i}\mathbf{1}_{n_{e}}=\alpha_{\mathcal{B}_{i}},~{}~{}~{}\gamma_{i}^{T}\mathbf{1}_{b}=\beta.$

Combining it with the definition of $\Gamma_{i}$, we have:

$\displaystyle\Gamma_{i}\mathbf{1}_{n_{e}}$
$\displaystyle=\bar{\alpha}_{\mathcal{B}_{i}}=\eta_{i}\odot\alpha,$ (13)
$\displaystyle\Gamma_{i}^{T}\mathbf{1}_{n_{o}}$ $\displaystyle=\beta,$ (14)

where $\bar{\alpha}_{\mathcal{B}_{i}}\in\mathbb{R}^{n_{o}}$ is the extension
of $\alpha_{\mathcal{B}_{i}}$ by padding zero entries to the dimension whose
index does not belong to $\mathcal{B}_{i}$, $\odot$corresponds to entry-wise
product and $\eta_{i}$ is a $n_{o}$ dimensional vector with element satisfying
that:

$\displaystyle\eta_{i}^{j}=\begin{cases}\frac{n_{o}}{b}=k,&if\
j\in\mathcal{B}_{i},\\\ 0,&otherwise.\end{cases}$ (15)

According Eq. (13) and (14), we have:

$\displaystyle(\frac{1}{k}\sum_{i=1}^{k}\Gamma_{i})\mathbf{1}_{n_{e}}$
$\displaystyle=\frac{1}{k}\sum_{i=1}^{k}(\Gamma_{i}\mathbf{1}_{n_{e}})=(\frac{1}{k}\sum_{i=1}^{k}\eta_{i})\odot\alpha,$
(16)
$\displaystyle(\frac{1}{k}\sum_{i=1}^{k}\Gamma_{i})^{T}\mathbf{1}_{n_{o}}$
$\displaystyle=\frac{1}{k}\sum_{i=1}^{k}(\Gamma_{i}^{T}\mathbf{1}_{n_{o}})=\frac{1}{k}\sum_{i=1}^{k}\beta=\beta,$

Combining Eq. (15) with conditions
$\mathcal{B}_{i}\cap\mathcal{B}_{j}=\emptyset$ and
$\cup_{i=1}^{k}\mathcal{B}_{i}=\\{1,2,...,n_{o}\\}$, we obtain that:

$\displaystyle(\frac{1}{k}\sum_{i=1}^{k}\eta_{i})^{j}=k\sum_{i=1}^{k}\mathbb{I}(j\in\mathcal{B}_{i})=k,$

which means that the term $\frac{1}{k}\sum_{i=1}^{k}\eta_{i}$ in Eq. (16) is
$\frac{1}{k}\sum_{i=1}^{k}\eta_{i}=\mathbf{1}_{n_{o}}$. As a result, we have
already proved Eq. (11) because the following definition holds:

$\displaystyle\left(\frac{1}{k}\sum_{i=1}^{k}\Gamma_{i}\right)\mathbf{1}_{n_{e}}=\alpha,~{}~{}~{}\left(\frac{1}{k}\sum_{i=1}^{k}\Gamma_{i}\right)^{T}\mathbf{1}_{n_{o}}=\beta.$

By denoting the cost matrix between $P(\mathbf{S},\mathbf{X},A|G=o)$ and
$P(\mathbf{S},\mathbf{X},A|G=e)$ as $C$, the above definition also leads to
the following inequality:

$\displaystyle\quad\;\;OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))$
$\displaystyle\leq\langle\frac{1}{k}\sum_{i=1}^{k}\Gamma_{i},C\rangle=\frac{1}{k}\sum_{i=1}^{k}\langle\Gamma_{i},C\rangle$
$\displaystyle=\frac{1}{k}\sum_{i=1}^{k}OT(P_{\mathcal{B}_{i}}(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e))$
(17)
$\displaystyle=m\mbox{-}OT(P(\mathbf{S},\mathbf{X},A|G=o),P(\mathbf{S},\mathbf{X},A|G=e)),$

where the Eq. (17) comes from the definition of $~{}\Gamma_{i}$. So far, we
have finished the proof of Eq. (12). ∎

###  A.5 Proof of Theorem 4

######  Theorem 4.

Assuming a family ${\mathcal{M}}$ of function
$m:{\mathcal{Z}}\times{\mathcal{A}}\rightarrow\mathbb{R}$, and there exists a
constant $B_{\phi}>0$ such that
$\frac{1}{B_{\Phi}}\cdot\ell_{\phi,g,h}(\mathbf{x},a)\in{\mathcal{M}}$, then
we have:

$\displaystyle{\mathcal{E}}_{cf}\leq{\mathcal{E}}_{f}^{\mathbf{w}_{o}}+IPM_{\mathcal{M}}(P(\mathbf{Z})P(A),~{}\mathbf{w}_{o}\cdot
P(\mathbf{Z},A)),$

where
$IPM_{{\mathcal{M}}}(p,q)=\sup_{m\in{\mathcal{M}}}\left|\int_{\mathcal{X}}g(x)(p(x)-q(x))dx\right|$
is integral probability metric for a chosen ${\mathcal{M}}$, $P(\mathbf{Z})$
and $P(\mathbf{Z},A)$ are distributions induced by the map $\phi$ from
$P(\mathbf{X})$ and $P(\mathbf{X},A)$.

###### Proof.

The inequality is equivalent to

$\displaystyle{\mathcal{E}}_{cf}-{\mathcal{E}}_{f}^{\mathbf{w}_{o}}$
$\displaystyle=\int_{\mathcal{X}}\int_{\mathcal{A}}P(\mathbf{x})P(a)\ell_{\phi,g,h}(\mathbf{x},a)\
d\mathbf{x}da-\int_{\mathcal{X}}\int_{\mathcal{A}}P^{\mathbf{w}_{o}}(\mathbf{x},a)\ell_{\phi,g,h}(\mathbf{x},a)\
d\mathbf{x}da$
$\displaystyle=\int_{\mathcal{X}}\int_{\mathcal{A}}\left(P(\mathbf{x})P(a)-P^{\mathbf{w}_{o}}(\mathbf{x},a)\right)\ell_{\phi,g,h}(\mathbf{x},a)\
d\mathbf{x}da$
$\displaystyle=\int_{\mathcal{Z}}\int_{\mathcal{A}}\left(P(\mathbf{z})P(a)-P^{\mathbf{w}_{o}}(\mathbf{z},a)\right)\ell_{\phi,g,h}(\psi(\mathbf{z}),a)\
d\mathbf{z}da$ (18) $\displaystyle\leqslant\sup_{m\in\mathcal{M}}\
\left|\int_{\mathcal{Z}}\int_{\mathcal{A}}\left(P(\mathbf{z})P(a)-P^{\mathbf{w}_{o}}(\mathbf{z},a)\right)m(\mathbf{z},a)\
d\mathbf{z}da\right|$ (19)
$\displaystyle=IPM_{\mathcal{M}}\left(P(\mathbf{Z})P(A),\
P^{\mathbf{w}_{o}}(\mathbf{Z},A)\right)$
$\displaystyle=IPM_{\mathcal{M}}(P(\mathbf{Z})P(A),\ \mathbf{w}_{o}\cdot
P(\mathbf{Z},A)),$

where Eq. (18) can be obtained by the standard change of variables formula,
using the determinant of the Jacobian of $\psi(\mathbf{z})$, Eq. (19) is
according to the definition of IPM with the assumption that
$\frac{1}{B_{\Phi}}\cdot\ell_{\phi,g,h}(\mathbf{x},a)\in{\mathcal{M}}$. ∎

##  Appendix B Derivation about Projected Mirror Descent

In this section, following [58], we develop a projected mirror descent [40,
43] based on the Kullback-Leibler (KL) divergence to solve the following
problem, which is non-trivial to address because of the equality constrints:

$\displaystyle\min_{\mathbf{w}_{o,\mathcal{B}}}\min_{\gamma\in\Pi(\alpha_{\mathcal{B}},\beta)}\langle\gamma,\mathbf{C}\rangle+\lambda_{e}\Omega(\gamma),$
$\displaystyle s.t.~{}$
$\displaystyle\Pi(\alpha_{\mathcal{B}},\beta)=\\{\gamma\in{\mathbb{R}}^{b\times
n_{e}}|\gamma\mathbf{1}_{n_{e}}=\mathbf{w}_{o,\mathcal{B}},\gamma^{T}\mathbf{1}_{b}=\bm{\mu},\gamma_{ij}\in\left[0,1\right]\\},$

Firstly, we rewrite the constrain about $\gamma$ as

$\displaystyle\Pi^{{}^{\prime}}(\alpha_{\mathcal{B}},\beta)=\\{\gamma\in{\mathbb{R}}^{b\times
n_{e}}|\gamma^{T}\mathbf{1}_{b}=\bm{\mu},\gamma_{ij}\in\left[0,1\right]\\},$
(20)

which does not consider the constraint
$\gamma\mathbf{1}_{n_{e}}=\mathbf{w}_{o,\mathcal{B}}$ since
$\mathbf{w}_{o,\mathcal{B}}$ are also parameters to be optimized. Based on
this, the OT problem with a negative entropy regularization is given as
follows:

$\displaystyle\min_{\gamma}\langle\gamma,\mathbf{C}\rangle+\lambda_{e}\Omega(\gamma),~{}~{}~{}s.t.~{}\gamma\in\Pi^{{}^{\prime}}(\alpha_{\mathcal{B}},\beta).$
(21)

For simplicity, we define the objective function in Problem (21) as

$\displaystyle
f(\gamma)=\langle\gamma,\mathbf{C}\rangle+\lambda_{e}\Omega(\gamma),$

and the $(i,j)$-th element of the gradient $\nabla f(\gamma)$ is denoted by
$\nabla_{ij}$,

$\displaystyle\nabla_{ij}=C_{ij}+\lambda_{e}\log
w_{i}=C_{ij}+\lambda_{e}\log\gamma_{i\cdot},$

where $\gamma_{i\cdot}$ is the sum of $i$-th row of $\gamma$. Then at each
iteration, we solve the following problem:

$\displaystyle\gamma^{k}=\min_{\gamma\in\Pi^{{}^{\prime}}(\alpha_{\mathcal{B}},\beta)}\eta\langle\nabla
f(\gamma^{k-1}),\gamma\rangle+\mathcal{D}(\gamma||\gamma^{k-1}),$ (22)

which firstly performs proximal gradient descent with the Bregman divergence
and the stepsize $\eta$, and then obtains a feasible solution in the set
$\Pi^{{}^{\prime}}(\alpha_{\mathcal{B}},\beta)$ by projection. Next, we
present the details of these two operations.

Proximal Gradient Descent Let $\upsilon^{k}$ be the solution to Problem (22)
without considering the constraint
$\gamma\in\Pi^{{}^{\prime}}(\alpha_{\mathcal{B}},\beta)$, i.e.,

$\displaystyle\upsilon^{k}=\min_{\gamma}~{}~{}\eta\langle\nabla
f(\gamma^{k-1}),\gamma\rangle+\mathcal{D}(\gamma||\gamma^{k-1}).$

We adopt the KL divergence between two distributions $\gamma$ and
$\gamma^{k-1}$ as the Bregman divergence $\mathcal{D}(\gamma||\gamma^{k})$,
then the closed-form solution to the above problem is given as:

$\displaystyle\upsilon^{k}=\gamma^{k-1}\odot\exp(-\eta\nabla
f(\gamma^{k-1})).$ (23)

Projection Operation To make sure $\gamma^{k}$ satisfies the constraints in
Eq. (20), we update $\gamma^{k}$ by finding
$\gamma\in\Pi^{{}^{\prime}}(\alpha_{\mathcal{B}},\beta)$ which is most close
to $\upsilon^{k}$ under the KL metric. This is achieved by solving the
following projection problem (ignore constraint
$\gamma_{ij}\in\left[0,1\right]$ for now):

$\displaystyle\min_{\gamma}\quad\mathcal{D}(\gamma||\upsilon^{k})=\sum_{i=1}^{b}\sum_{j=1}^{n_{e}}\gamma_{ij}\log(\frac{\gamma_{ij}}{\upsilon^{k}_{ij}})-\gamma_{ij}+\upsilon^{k}_{ij},~{}~{}~{}s.t.~{}~{}\gamma^{\top}\mathbf{1}_{b}=\bm{\mu}.$
(24)

By introducing the Lagrangian multipliers
$\bm{\lambda}=[\lambda_{1},\ldots,\lambda_{n_{t}}]^{\top}$ for the equality
constraint $\gamma^{\top}\mathbf{1}_{b}=\bm{\mu}$, we obtain the Lagrangian
$\mathcal{L}(\gamma,\bm{\lambda})$ as follows:

$\displaystyle\mathcal{L}(\gamma,\bm{\lambda})=\sum_{i=1}^{b}\sum_{j=1}^{n_{e}}\gamma_{ij}\log(\frac{\gamma_{ij}}{\upsilon^{k}_{ij}})-\gamma_{ij}+\upsilon^{k}_{ij}+\bm{\lambda}^{\top}(\gamma^{\top}\mathbf{1}_{b}-\bm{\mu}).$

By taking the partial derivative of $\mathcal{L}(\gamma,\bm{\lambda})$ with
respect to $\gamma_{ij}$ to zero, we obtain:

$\displaystyle\log\gamma_{ij}$
$\displaystyle=\log\upsilon^{k}_{ij}-\lambda_{j}$
$\displaystyle\Rightarrow\gamma_{ij}$
$\displaystyle=\upsilon^{k}_{ij}\exp(-\lambda_{j}).$ (25)

According to the equality constraint $\gamma^{\top}\mathbf{1}_{b}=\bm{\mu}$,
we have $\sum_{i=1}^{b}\gamma_{ij}=1/n_{e}$. Combining it with the above
result, we further obtain:

$\displaystyle\sum_{i=1}^{b}\gamma_{ij}=\sum_{i=1}^{b}\upsilon^{k}_{ij}\exp(-\lambda_{j})$
$\displaystyle=\exp(-\lambda_{j})\sum_{i=1}^{b}\upsilon^{k}_{ij}=\frac{1}{n_{e}}$
$\displaystyle\Rightarrow\exp(-\lambda_{j})$
$\displaystyle=\frac{1}{n_{e}\sum_{i=1}^{b}\upsilon^{k}_{ij}}.$ (26)

Combining Eq. (25) and (26), we obtain the closed-form solution:

$\displaystyle\gamma_{ij}=\upsilon^{k}_{ij}\Big{/}\Big{(}n_{e}\sum_{i=1}^{b}\upsilon^{k}_{ij}\Big{)}.$
(27)

For an initial value $\gamma_{ij}^{0}\geq 0$, given $\upsilon^{k}_{ij}>0$
which is guaranteed by the update rule in Eq. (20), it is obvious that the
solution obtained by Eq. (27) satisfies the box constraint
$\gamma_{ij}\in[0,1]$. Therefore, Problem (24) does not consider this
constraint explicitly.

For now, we have derived how to solve the Problem Eq. (21). That is, we repeat
steps Eq. (23) and (27) until converging.

##  Appendix C Additional Related Work

Since our work is also close to work that considers continuous treatment. We
provide additional related work to them as follows.

Dose-Response Curve Estimation Dose-response curve estimation is broadly
studied through statistical methods [51, 35, 61], and most of them are based
on the generalized propensity score (GPS) [31] or the entropy balancing for
continuous treatments (EBCT) [49]. Recently, deep learning-based methods also
attract the attention of the research community. DRNet [45] divides the
continuous treatment into several intervals, and trains one separate head for
each interval. In order to produce a continuous curve, VCNet [41] adapts a
varying coefficient model into neural networks to handle the continuous
treatment. SCIGAN [7] generates counterfactual outcomes for continuous
treatments based on the generative adversarial network (GAN) framework. ADMIT
[56] derives a counterfactual bound of estimating the average dose-response
curve with a discretized approximation of the IPM distance, and ACFR [34]
minimizes the KL-divergence using an adversarial game to extract balanced
representations for continuous treatment. Different from existing works, we
aim to estimate the HDRC under the long-term estimation scenario in this
paper.

Algorithm 1 Pseudocode of the proposed LEARN

1:Observational dataset $O$, experimental dataset $E$, initial model
parameters
$\theta=\\{\theta_{\phi},\theta_{q}(a),\theta_{g}(a),\theta_{f}(a),\theta_{h}(a)\\}$,
hyper-parameters $\\{\lambda_{o},\lambda_{b},\lambda_{e}\\}$, batch size $b$,
pre-training epoch $n_{pre}$, training epoch $n$, OT epoch $n_{ot}$, learning
rate $\alpha$

2:// Pre-training phase. Ensuring the short-term embeddings of RNN are
meaningful, which will be used in constructing the cost matrix $\mathbf{C}$ in
optimal transport at each training iteration

3:for $i=1,2,...,n_{pre}$ do

4: Sample a mini-batch $\mathcal{B}$ from the observational dataset $O$

5: Update
$\theta^{i+1}\leftarrow\theta^{i}-\alpha\nabla_{\theta}\mathcal{L}_{\theta}$
by setting
$\mathbf{w}_{\mathcal{B}}=\\{\frac{1}{b},\frac{1}{b},...,\frac{1}{b}\\}$.

6:end for

7:// Training phase. Removing the observed and unobserved confounding bias

8:for $i=1,2,...,n$ do

9: Sample a mini-batch $\mathcal{B}$ from the observational dataset $O$

10: // Learning weights $\mathbf{w}_{\mathcal{B}}$ for the batch observational
units via data combination

11: Compute the cost matrix $\mathbf{C}$ between $\mathcal{B}$ and $E$

12: for $k=1,2,...,n_{ot}$ do

13: Calculate $\upsilon^{k}$ according to Eq. (23)

14: Update $\gamma^{k}$ according to Eq. (27)

15: end for

16: $\hat{\mathbf{w}}_{\mathcal{B}}=\gamma^{k}\mathbf{1}_{n_{e}}$

17: // Weighted regression with IPM term

18: Update
$\theta^{i+1}\leftarrow\theta^{i}-\alpha\nabla_{\theta}\mathcal{L}_{\theta}$
by setting $\mathbf{w}_{\mathcal{B}}=\hat{\mathbf{w}}_{\mathcal{B}}$.

19:end for

20:return Learned model parameters $\theta$

##  Appendix D Model Details

Implementation Details Our code is implemented using the PyTorch 1.8.1
framework. Note that all MLP in our model is 2 fully connected layers with 50
hidden units and Elu activation function, then we build our model as follows:

For the balanced representation module, we use an MLP to build the
representation function $\phi$, and use the Wasserstein distance as the IPM
term. The code of the Wasserstein distance can be found in
https://github.com/rguo12/network-deconfounder-wsdm20/blob/master/utils.py.

For the OT weight module, at each iteration, we first construct the cost
matrix as
$C_{ij}=\beta_{r}d(\bar{\mathbf{r}}_{i},\bar{\mathbf{r}}_{j})+\beta_{x}d(\mathbf{x}_{i},\mathbf{x}_{j})+\beta_{a}d(a_{i},a_{j})$,
where the hyperparameters $\beta_{r}=10,\beta_{x}=0.1,\beta_{a}=0.1$,
$d(\cdot,\cdot)$ is the euclidean distance, and the distance of short-term
outcomes is calculated from the mean embedding of RNN as
$\bar{\mathbf{r}}=\frac{1}{t_{0}}\sum_{t=1}^{t_{0}}\mathbf{r}_{t}$. After
that, we obtain the optimal transport probability matrix $\gamma^{*}$ by
repeating steps Eq. (23) and (27) in Appendix B until converging. Finally, the
weights of observational samples could be obtained as
$w_{i}=\sum_{j=1}^{n_{e}}\gamma^{*}_{ij}$, where $w_{i}$ is the sum of the
$i$-th row of $\gamma^{*}$.

For the long-term estimation module, we use a shared MLP to build the function
$g(\cdot)$ predicting the short-term outcomes, use another MLP to build the
function $h(\cdot)$ predicting the long-term outcome, and use a bidirectional
GRU model to build the RNN $q(\cdot)$ extracting the short-term
representations. As for the attention mechanism $f(\cdot)$, we implement it as
follows:

$\displaystyle\tilde{\mathbf{R}}_{t}=\tanh{(\mathbf{W}\mathbf{R}_{t}+\mathbf{b})},~{}~{}~{}\alpha_{t}=\frac{\exp{\tilde{\mathbf{R}}_{t}^{\top}\mathbf{V}}}{\sum_{t=1}^{t_{0}}\exp{\tilde{\mathbf{R}}_{t}^{\top}\mathbf{V}}},~{}~{}~{}\mathbf{R}_{T}=\sum_{t=1}^{t_{0}}\alpha_{t}\mathbf{R}_{t},$

where $\mathbf{W},\mathbf{b},\mathbf{V}$ are the parameters of the attention
mechanism $f(\cdot)$.

Further, following [41], we model the parameters $\theta(a)$ of the above four
network structures $q(\cdot),g(\cdot),f(\cdot),h(\cdot)$ as the varying
coefficient structure, i.e., $\theta(a)=\sum_{l=1}^{L}\alpha_{l}\phi_{l}(a)$,
where $\phi_{l}$ is the basis function and $\alpha_{l}$ is the coefficient.
Specifically, we use truncated polynomial basis with degree 2 and two knots at
$\\{1/3,2/3\\}$ (altogether 5 basis) to build the varying coefficient
structure. The code of the varying coefficient structure is written based on
VCNet (https://github.com/lushleaf/varying-coefficient-net-with-functional-
tr/blob/main/models/dynamic_net.py).

Parameters Setting We set the hyperparameters in our model as follows: the
strength of the Wasserstein distance $\lambda_{b}\in\\{50,100,150,200\\}$, the
strength of negative entropy regularization $\lambda_{b}=100$, the learning
rate about the projected mirror descent as 0.001, and the iterations about the
projected mirror descent as $\\{50,100,500\\}$, the strength of short-term
outcome loss as $\lambda_{o}\in\\{0.25,0.5,0.75\\}$. Besides, we pre-train the
model for 100 epochs without weighting, and then the number of training epochs
is 400 with an early stop. We use Adam optimizer [36] with a learning rate of
$1e-3$ and weight decay of $5e-4$. All the experiments are run on the Tesla
P40 GPU.

Pseudocode We provide the pseudocode of our model in Algorithm 1.

##  Appendix E Experimental Details

###  E.1 Dataset Generation

Synthetic Data Generation We simulate synthetic data as follows. For each
observational unit $i\in\\{1,2,...,n_{o}\\}$ with $G=o$, we generate $p$
observed covariates and $q$ unobserved covariates from independent identical
distributions, i.e., $X^{o}\sim{\mathcal{N}}(\mathbf{0.1}_{p},\mathbf{I}_{p})$
and $U\sim{\mathcal{N}}(\mathbf{0.25}_{q},\mathbf{I}_{q})$, where $\mathbf{I}$
denotes the identity matrix. Similarly, we generate experimental units
$i\in\\{1,2,...,n_{e}\\}$ with $G=e$ as
$X^{e}\sim{\mathcal{N}}(\mathbf{0.5}_{p},\mathbf{I}_{p}),U\sim{\mathcal{N}}(\mathbf{0.25}_{q},\mathbf{I}_{q})$.
Then we generate treatment and outcomes for two types of datasets using the
following rules with
$W^{t}_{j}\sim\mathcal{U}(0,0.5),W^{y}_{j}\sim\mathcal{U}(0.5,1)$:

$\displaystyle
A=\sigma(\sum_{j=1}^{p/3}W_{j}^{t}\sin{X_{j}}+\sum_{j=p/3}^{2p/3}$
$\displaystyle W_{j}^{t}X_{j}^{2}+\beta_{U}\bar{U}\mathbb{I}(G=o)),$
$\displaystyle
Y_{t}=\sum_{j=p/3}^{2p/3}W_{j}^{y}X_{j}^{2}+(t+5)A(\sum_{j=2p/3}^{p}W_{j}^{y}X_{j})$
$\displaystyle+\frac{\beta_{U}tA}{5}(\sum_{j=1}^{q}W_{j}^{y}\cos{U_{j}})+0.25\bar{Y}_{1:t-1}+\mathcal{N}(0,0.5),$

where $n_{o}=10000,n_{e}=500,p=15,q=5$. We denote short-term outcomes as
$\mathbf{S}=\\{Y_{1},Y_{2},...,Y_{7}\\}$ and long-term outcome as
$Y_{T}=Y_{14}$. Besides, we design 5 simulation datasets with varying
$\beta_{U}\in\\{1,1.25,1.5,1.75,2\\}$ that controls different strengths of
unobserved confounding bias.

Semi-synthetic Data Generation In semi-synthetic datasets News and TCGA, we
reuse their covariates. Firstly, we randomly partition the dataset into
observational ($G=o$) and experimental ($G=e$) subsets at a 9:1 ratio,
followed by segregation of original covariates into observed ($X$) and
unobserved ($U$) ones at an 8:2 ratio. Then we follow [41, 56] to generate the
assigned treatments and their corresponding outcomes. After generating a set
of parameters $\bm{v}_{i,x/u}=\bm{u}_{i,x/u}/||\bm{u}_{i,x/u}||$ and
$i=1,2,3$, where $\bm{u}_{i,x/u}$ is sampled from a normal distribution
$\mathcal{N}(\mathbf{0},\mathbf{1})$, we have:

$\displaystyle A\sim
Beta(\gamma,\beta),\textit{where}~{}\gamma=2,\beta=\frac{\gamma-1}{d^{*}}$
$\displaystyle+2-\gamma~{}\textit{and}~{}d^{*}=\left|\frac{\bm{v}_{3,x}^{T}X}{2\bm{v}_{2,x}^{T}X}+\frac{\bm{v}_{3,u}^{T}U}{2\bm{v}_{2,u}^{T}U}\mathbb{I}(G=o)\right|,$
$\displaystyle
Y_{t}\sim\mathcal{N}(\mu,0.5),\textit{where}~{}\mu=4(A-0.5)^{2}\times$
$\displaystyle\sin{(\frac{\pi}{2}A)}\times
2(\max(-2,\exp(\frac{\bm{v}_{2,x}^{T}X}{\bm{v}_{3,x}^{T}X}+\frac{\bm{v}_{2,u}^{T}U}{\bm{v}_{3,u}^{T}U}-0.3))$
$\displaystyle+20A(\bm{v}_{1,x}^{T}X+\bm{v}_{1,u}^{T}U))$
$\displaystyle+0.5\bar{Y}_{1:t-1}.$

Same as the synthetic data generation, we collect
$\mathbf{S}=\\{Y_{1},Y_{2},...,Y_{7}\\}$ and $Y_{T}=Y_{14}$.

###  E.2 Hyperparameter Sensitivity Study

We investigate the role of three hyperparameters in our model:
$\lambda_{o},\lambda_{b},\lambda_{e}$, and the results are shown in Fig. 3.
Firstly, we observe that MISE increases when $\lambda_{b}$ becomes small or
large, indicating that there exists a trade-off between eliminating observed
confounding bias and estimating outcomes accurately. Secondly, the results of
varying $\lambda_{o}$ shows that considering the short-term outcomes of
observational and experimental datasets together can bring a certain degree of
improvement. Finally, applying the negative entropy regularization in OT can
bring a slight improvement, and the performance of our model with different
strengths is relatively stable.

(a) Impact of $\lambda_{b}$

(b) Impact of $\lambda_{o}$

(c) Impact of $\lambda_{e}$

Figure 3: Hyperparameter sensitivity. (a) Strength of IPM. (b) Strength of
short-term loss between observational and experimental datasets. (c) Strength
of negative entropy regularization.

##  Appendix F Broader Impact and Limitations

This paper studies the problem of estimating the long-term HDRC, which
utilizes the power of deep learning and optimal transport technologies to
assist better long-term decision-making in many domains. For example, in a
ride-hailing platform, the platform needs to evaluate the effect of various
incomes of different drivers on their retention after one year to keep the
long-term balance of demand and supply. In the medical field, a healthcare
institution might need to analyze the impact of a new medication on the five-
year survival rates of cancer patients to determine whether to adopt it as a
standard treatment protocol. However, this technology towards better
estimating and understanding long-term HDRC can be used negatively, where
someone wishing to cause harm can use the estimated outcomes to select the
worst outcome.

At the same time, although this work allows there exist unobserved confounders
in the observational dataset, some standard assumptions in causal inference
are still needed. For example, this work assumes the Positivity Assumption
holds, and this assumption may be violated in some real-world scenarios when
there exist units that only receive treatment $A=a$, e.g., in some cases,
critically ill patients have to be treated, leading to no samples of untreated
critically ill patients. Also, this work still relies on the existence of an
experimental dataset to remove unobserved confounding, and how to estimate
long-term HDRC only with a confounded observational dataset can be further
explored.