Thong Pham<EMAIL_ADDRESS>
_Data Science and AI Innovation Research Promotion Center, Shiga University,
Japan_
_Center for Advanced Intelligence Project, RIKEN, Japan_
Shohei Shimizu<EMAIL_ADDRESS>
_Graduate School of Data Science, Shiga University, Japan_
_Center for Advanced Intelligence Project, RIKEN, Japan_
Hideitsu Hino<EMAIL_ADDRESS>
_The Institute of Statistical Mathematics_
_Center for Advanced Intelligence Project, RIKEN, Japan_
Tam Le<EMAIL_ADDRESS>
_The Institute of Statistical Mathematics_
_Center for Advanced Intelligence Project, RIKEN, Japan_

# Scalable Counterfactual Distribution Estimation
in Multivariate Causal Models

###### Abstract

We consider the problem of estimating the counterfactual joint distribution of
multiple quantities of interests (e.g., outcomes) in a multivariate causal
model extended from the classical difference-in-difference design. Existing
methods for this task either ignore the correlation structures among
dimensions of the multivariate outcome by considering univariate causal models
on each dimension separately and hence produce incorrect counterfactual
distributions, or poorly scale even for moderate-size datasets when directly
dealing with such multivariate causal model. We propose a method that
alleviates both issues simultaneously by leveraging a robust latent one-
dimensional subspace of the original high-dimension space and exploiting the
efficient estimation from the univariate causal model on such space. Since the
construction of the one-dimensional subspace uses information from all the
dimensions, our method can capture the correlation structures and produce good
estimates of the counterfactual distribution. We demonstrate the advantages of
our approach over existing methods on both synthetic and real-world data.

###### keywords:

multivariate counterfactual distribution, optimal transport, difference in
difference

##  1 Introduction

Causal inference has received explosive interest in the last decades, due to
the need to extract causal knowledge from data in various research fields,
such as statistics (Pearl, 2009), sociology (Gangl, 2010), biomedical
informatics (Kleinberg and Hripcsak, 2011), public health (Glass et al.,
2013), and machine learning (Schölkopf, 2022). One of the most popular causal
inference models in practice is the difference in difference (DiD) model,
which dates back to the works of John Snow in the 1850s (Snow, 1854, 1855;
Donald and Lang, 2007; Lechner, 2011; Roth et al., 2022). In this model, we
observe a quantity of interest (i.e., outcome) from two different groups,
i.e., the treatment group and the control group, at two different time steps,
i.e., before and after the intervention (i.e., treatment) event. More
precisely, the intervention is only applied on the treatment group while there
is no intervention applied on the control group. We assume that if there was
no intervention in the treatment group, its outcome variable would evolve the
same way as that of the control group. This is the so-called “ _parallel
trend_ ” assumption. The classical DiD model, however, requires that the
parallel trend is linear, i.e., the means of the outcome variable in the two
groups must evolve the same way when there is no intervention. This may limit
its application in practice. To extend the DiD model for non-linear settings,
several proposals have been developed in the literature (Abadie, 2005; Athey
and Imbens, 2006; Blundell and Costa Dias, 2009; Sofer et al., 2016), of which
notably is the Changes-in-Changes (CiC) model (Athey and Imbens, 2006).

The CiC model generalizes the DiD model to include non-linear parallel trends
that can act on the whole distribution of the outcome, i.e., in the absence of
intervention, the means of the outcome in control and treatment groups are
allowed to evolve in different ways, as long as the two outcome
_distributions_ evolve in the same way. This allows identifications of more
complex treatment effects that require information from the whole outcome
distributions, not just their first moments (Lechner, 2011). The standard CiC
model, however, is only designed for univariate outcomes. In order to extend
the CiC model for a multivariate outcome variable, a naive approach is to
tensorize univariate CiC models, i.e., considering independently a univariate
CiC model for each coordinate. Nevertheless, this naive approach fails to
capture correlations among coordinates of the outcome and thus is incapable of
modelling complex, multivariate parallel trends.

Recently, by leveraging the optimal transport (OT) theory (Villani, 2003,
2008), Torous et al. (2021) proved that the counterfactual outcome
distribution of the treatment group in the CiC model (i.e., the outcome
distribution of the treatment group _without_ receiving intervention) at the
post-intervention time stamp can be estimated by exploiting the _optimal
transport map_ which pushforwards the outcome distribution of the control
group at the pre-intervention time stamp to that at the post-intervention time
stamp. Consequently, it is natural to extend the CiC model for univariate
quantity of interest into that for multivariate one through the lens of OT,
since this would take into account the dependence structure of the dimensions
and be able to model complex parallel trends.

However, OT suffers a few drawbacks. It has a high computational complexity,
which is super cubic with respect to the number of supports of the input
distribution. A popular approach is to rely on the entropic regularization for
OT, a.k.a., Sinkhorn (Cuturi, 2013), to reduce its computational complexity to
quadratic. Yet, Sinkhorn yields a dense estimator for the optimal transport
plan, which is not a desirable property for counterfactual estimation in the
CiC model. Additionally, OT has a high sample complexity, i.e.,
${\mathcal{O}}(n^{-1/d})$ where $n$ is the number of samples and $d$ is the
dimension of samples in the probability measures.

In this work, in order to exploit the efficient computation of the CiC model
for the univariate quantity of interest, and alleviate the above-mentioned
challenges of the OT approach for the multivariate CiC model, we propose to
leverage the max-min robust OT approach (Paty and Cuturi, 2019). In
particular, we propose to lift the univariate CiC model to that for a
multivariate quantity of interest by seeking a _robust_ latent univariate
subspace. Unlike the naive tensorization approach, our approach can
incorporate the correlations of coordinates. Moreover, unlike the standard OT
approach as in (Torous et al., 2021), our estimator can preserve the efficient
computation as in the univariate CiC model since the optimal transport plan is
estimated on the robust latent one-dimensional subspace instead of its
original high-dimensional space.

Intuitively, our approach follows the _max-min robust OT_ approach (Paty and
Cuturi, 2019), which steams from the robust optimization (Ben-Tal et al.,
2009; Bertsimas et al., 2011) where there are uncertainty non-stochastic
parameters. The robust optimization has many roots in applied sciences, e.g.,
in robust control (Keel et al., 1988), machine learning (Morimoto and Doya,
2001; Xu et al., 2009; Panaganti et al., 2022). In the context of OT, several
advantages of the max-min (and its relaxation min-max) robust OT have been
reported. For example, (i) it makes the OT approach robust to noise (Paty and
Cuturi, 2019; Dhouib et al., 2020); and (ii) it also helps to reduce the
sample complexity (Paty and Cuturi, 2019; Deshpande et al., 2019). At a high
level, our contributions are two-fold as follows:

  * •

(i) We propose a max-min robust OT approach for the multivariate CiC model.
The proposed approach not only inherits properties of the OT approach for the
CiC model but also preserves the efficient computation as in the univariate
CiC model.

  * •

(ii) We evaluate our approach on both synthesized and real data to illustrate
the advantages of the proposed method.

The paper is organized as follows. After reviewing the multivariate CiC model
and existing methods for estimating the counterfactual distribution in this
model in Section 2, we discuss our proposed method in Section 3 and
demonstrate its benefits through synthetic data in Section 4. We apply it to
the classical dataset of Card and Krueger (1993) in Section 5 before giving
concluding remarks in Section 6.

Notations. We use the superscripts C and T to indicate the control group and
treatment group, respectively. We drop those superscripts when either the
context is clear or it is not necessary to distinguish these two groups.

##  2 The Causal Model

In this section, we describe the CiC causal model for multiple quantities of
interests (Torous et al., 2021), which is an extension based on OT theory from
the original, univariate model (Athey and Imbens, 2006).

We use a stochastic process to model the quantity of interests (i.e.,
outcomes) before the intervention, i.e., at the time stamp $t=0$, and after
the intervention, i.e., at the time stamp $t=1$. We denote them as
$\\{Y_{t}\\}_{t=0,1}$ where $Y_{t}$ is in the $\mathbb{R}^{d}$ space. For the
original CiC causal model, we have $d=1$ (Athey and Imbens, 2006). We let
$\mu_{t}$ be the distribution of $Y_{t}$ for $t=0,1$. Without loss of
generality, we assume that $Y_{t}$ is generated from a latent variable
$U_{t}\in\mathbb{R}^{d}$ which may change over time, but the distribution
$\nu$ of the latent variable $U_{t}$ is not changed over time, i.e., time-
invariant $U_{t}\sim\nu$ for $t=0,1$. Intuitively, $U_{t}$ may be regarded as
(unobserved) intrinsic features of a sample.

In the CiC model, we observe two groups:

  * •

(i) the control group: the stochastic process $\\{Y_{t}^{\texttt{C}}\\}$ is
solely affected by the _natural drift_. The evolution from
$\\{Y_{0}^{\texttt{C}}\\}$ at $t=0$ to $\\{Y_{1}^{\texttt{C}}\\}$ at $t=1$ is
independent of the treatment effects (of intervention).

  * •

(ii) the treatment group: the stochastic process $\\{Y_{t}^{\texttt{T}}\\}$ is
affected by both the _natural drift_ and the _treatment effects_.

For the CiC causal model, the goal of our causal inference is to deconvolve
the _natural drift_ and the _treatment effects_ in the treatment group. For
example, we would like to estimate the counterfactual distribution of the
control group at post-intervention under _only natural drift effect_ (i.e.,
without the treatment effects). By doing so, we can estimate the _treatment
effects_ of the intervention for the considered groups in application domains.

###  2.1 The Natural Drift Model

Natural drift is best explained in the stochastic process
$\\{Y_{t}^{\texttt{C}}\\}$ for $t=\\{0,1\\}$ since this process involves
solely the natural drift and is not affected by the treatment effect. The
change of $\\{Y_{t}^{\texttt{C}}\\}$ from the pre-intervention ($t=0$) to the
post-intervention ($t=1$) is modeled by assuming the existence of two
_production functions_ $h_{t}:\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$ with
$t\in\\{0,1\\}$ such that

$Y_{t}^{\texttt{C}}=h_{t}\left(U_{t}^{\texttt{C}}\right).$

Consequently, we have $\mu_{t}^{\texttt{C}}=(h_{t})_{\sharp}\nu^{\texttt{C}}$
where we introduce a new notation $\sharp$ as the _pushforward_ operator which
is defined as for any measurable set $A\subseteq\mathbb{R}^{d}$,
$\nu(A)=\mu(h_{t}^{-1}(A))$ (Peyré and Cuturi, 2019, Def. 2.1). In other
words, the distribution of the quantity of interests in the control group at
the time stamp $t$ (i.e., $\mu_{t}^{C}$) is the pushforward of the
distribution of the latent variable $U_{t}^{\texttt{C}}$ of the control group
at the time stamp $t$ (i.e., $\nu^{\texttt{C}}$) by the production function
$h_{t}$.

#### Natural drift map.

Assume that the production function $h_{0}$ is invertible, we have
$\nu^{\texttt{C}}=(h_{0}^{-1})_{\sharp}\mu_{0}^{C}$. Additionally, the
distribution $\nu^{\texttt{C}}$ of $U_{t}^{\texttt{C}}$ is time-invariant.
Thus, we have

$\mu_{1}^{\texttt{C}}=(h_{1}\circ h_{0}^{-1})_{\sharp}\mu_{0}^{\texttt{C}}.$
(1)

Equivalently, $\mu_{1}^{\texttt{C}}$ is the pushforward of
$\mu_{0}^{\texttt{C}}$ by the _natural drift map_

$\texttt{f}=h_{1}\circ h_{0}^{-1}.$ (2)

#### Natural drift in the treatment group.

We first introduce the concept of a _counterfactual distribution_ of the
outcome in this group, which is the outcome variable of the treatment group at
post-intervention $t=1$ under the purely hypothetical situation that the
treatment was never applied. Denote $\\{Y_{1}^{\texttt{T}^{*}}\\}$ as this
hypothetical stochastic process for the outcome and its distribution as
$\mu_{1}^{\texttt{T}^{*}}$. It is assumed that the change from
$\\{Y_{0}^{\texttt{T}}\\}$ to $\\{Y_{1}^{\texttt{T}^{*}}\\}$ is governed by
the same production functions $h_{0}$ and $h_{1}$ used in modeling natural
drift in the control group, namely $Y_{0}^{\texttt{T}}=h_{0}(U_{0}^{T})$ and
$Y_{1}^{\texttt{T}^{*}}=h_{1}(U_{1}^{T})$. This assumption generalizes the
parallel trend assumption in the classical DiD model. Then, similarly to Eq.
(1), we have

$\mu_{1}^{\texttt{T}^{*}}=(h_{1}\circ
h_{0}^{-1})_{\sharp}\mu_{0}^{\texttt{T}}=(\texttt{f})_{\sharp}\mu_{0}^{\texttt{T}},$
(3)

i.e., the counterfactual distribution $\mu_{1}^{\texttt{T}^{*}}$ is the
pushforward of $\mu_{0}^{\texttt{T}}$ by the natural drift map f.

Equations (1) and (3) suggest the following two-step method to estimate
$\mu_{1}^{\texttt{T}^{*}}$, the counterfactual distribution of outcomes in the
treatment group under only effects of the natural drift:

  * •

(i) We first estimate the natural drift map f from observed samples in the
stochastic processes $\\{Y_{0}^{\texttt{C}}\\},\\{Y_{1}^{\texttt{C}}\\}$ of
the distributions $\mu_{0}^{\texttt{C}},\mu_{1}^{\texttt{C}}$ respectively.
How to perform this estimation will be discussed in the next section.

  * •

(ii) We then use the estimated natural drift map f as the pushforward function
for the stochastic process $\\{Y_{0}^{\texttt{T}}\\}$ with distribution
$\mu_{0}^{\texttt{T}}$ to obtain an estimate for the distribution
$\mu_{1}^{\texttt{T}^{*}}$.

The schematic summary of the causal model is shown in Fig. 1.

Figure 1: Schematic of the multivariate CiC model. The latent distributions
$\nu^{\texttt{C}}$ and $\nu^{\texttt{T}}$ are time-invariant. At each time
$t$, the production function $h_{t}$ is applied to each latent distribution to
produce the corresponding outcome distribution in the absence of intervention.
Consequently, the natural drift map $\mathtt{f}=h_{1}\circ h_{0}^{-1}$
dictates the evolutions of the outcome distribution in both control and
treatment groups. This generalizes the linear parallel trend assumption in the
classical DiD model. For estimating treatment effects on the treatment group,
our goal is to estimate $\mu_{1}^{\texttt{T}^{*}}$, the counterfactual
distribution of outcomes in the treatment group at $t=1$. One can estimate
$\mathtt{f}$ from the observable empirical versions of $\mu_{0}^{\texttt{C}}$
and $\mu_{1}^{\texttt{C}}$ then use this estimation of $\mathtt{f}$ to push-
forward the observable empirical version of $\mu_{0}^{\texttt{T}}$ to obtain
an estimation of $\mu_{1}^{\texttt{T}^{*}}$.

If we further assume $U_{0}^{\texttt{C}}=U_{1}^{\texttt{C}}$, then we have
$Y_{1}^{\texttt{C}}=\texttt{f}\left(Y_{0}^{\texttt{C}}\right)$, i.e., the
natural drift map can be estimated by regress the control group at $t=1$ on
the control group at $t=0$. However, it requires coupled observations
$\left(Y_{0}^{\texttt{C}},Y_{1}^{\texttt{C}}\right)$, which may be not
available in practical applications, e.g., single-cell RNA-Seq data.

######  Remark 2.1.

For $x\in\mathbb{R}^{d}$, denote $x^{i}$ ($i=1,\cdots,d$) as the $i$-th
coordinate of $x$. The multivariate CiC model is equivalent to the naive
tensorization of univariate CiC models when (i) each production function
$h_{t}:\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$ ($t\in\\{0,1\\}$) can be
decomposed as $h_{t}(x)=[h_{t}^{1}(x^{1})\cdots h_{t}^{d}(x^{d})]^{T}$ for
univariate functions $h_{t}^{i}:\mathbb{R}\mapsto\mathbb{R}$, and (ii) each
coordinate of the latent variable (in both control and treatment groups) is
independent. Therefore, it is difficult for the tensorization of univariate
CiC models to express complex, multivariate natural drifts.

#### The CiC causal model and optimal transport.

For the uncoupled observations $Y_{0}^{\texttt{C}},Y_{1}^{\texttt{C}}$, given
$d=1$, the original CiC estimator (Athey and Imbens, 2006) assumes that
$h_{0}$ and $h_{1}$ are _monotone increasing_. Let $F_{t}^{\texttt{C}}$ be the
cumulative distribution function (cdf) of $Y_{t}^{\texttt{C}}$ for $t=0,1$.
Then, the unique monotone increasing _natural drift map_ is given by
$\left(F_{1}^{\texttt{C}}\right)^{-1}\circ F_{0}^{\texttt{C}}$ such that
$\mu_{1}^{\texttt{C}}=\left(\left(F_{1}^{\texttt{C}}\right)^{-1}\circ
F_{0}^{\texttt{C}}\right)_{\sharp}\mu_{0}^{\texttt{C}}$. Interestingly, the
natural drift map $\left(\left(F_{1}^{\texttt{C}}\right)^{-1}\circ
F_{0}^{\texttt{C}}\right)$ is also the optimal transport map between
$\mu_{0}^{\texttt{C}}$ and $\mu_{1}^{\texttt{C}}$. Therefore, the OT theory
provides a natural framework to extend the Changes-in-Changes causal model for
multivariate outcomes (Torous et al., 2021).

###  2.2 The OT approach for estimating the counterfactual distribution in
multivariate CiC

As briefly discussed in Section 2.1, for the univariate outcome, the OT map is
the natural drift map estimation for CiC causal model. Torous et al. (2021)
leveraged the OT theory to extend it for the multivariate quantity of
interests. We will discuss some mathematical background of OT for estimating f
in high dimensional setting for the CiC causal model.

Let $\mu,\nu$ be two probability measures supported on $\mathbb{R}^{d}$, the
OT problem between $\mu$ an $\nu$ with squared Euclidean ground cost is
defined as

$\text{OT}(\mu,\nu)=\min_{T}\int_{\mathbb{R}^{d}}\left\lVert
x-T(x)\right\rVert_{2}^{2}\,\text{d}\mu\quad\text{s.t.}\quad(T)_{\sharp}\mu=\nu,$
(4)

where $T$ is known as the transport plan.

#### For univariate model ($d=1$).

Given a probability measure $\mu$ supported on $\mathbb{R}$, we define its
cumulative distribution function (cdf) $F_{\mu}$ as

$F_{\mu}(x)=\mu((-\infty,x]).$ (5)

Note that the cdf is not always invertible, since it is not strictly
increasing. The pseudo-inverse of cdf $F^{-1}:[0,1]\mapsto\mathbb{R}$ is given
by

$F^{-1}(x)=\inf\left\\{t\in\mathbb{R}\mid F(t)\geq x\right\\}.$ (6)

When $d=1$, the OT admits a closed-form expression for the optimal map $T$ as
follows:

$T=\left(F_{\nu}\right)^{-1}\circ F_{\mu}.$ (7)

Note that $T$ is the unique increasing map such that
$(T)_{\sharp}{\mu}=\nu$.111The optimal condition for OT (i.e., existence of
the optimal map $T$) was described in (Gangbo, 1999) for $d=1$. Therefore, if
$\mu=\mu_{0}^{\texttt{C}}$, $\nu=\mu_{1}^{\texttt{C}}$, and f is increasing,
e.g., by choosing monotone $h_{0}$ and $h_{1}$, then f is the OT map between
$\mu_{0}^{\texttt{C}}$ and $\mu_{1}^{\texttt{C}}$.

#### For multivariate model ($d>1$).

The natural drift map f can be estimated via the OT map between
$\mu_{0}^{\texttt{C}}$ and $\mu_{1}^{\texttt{C}}$ (Torous et al., 2021).222The
optimal condition for OT (with $d>1$) was derived in (Brenier, 1991). However,
for high-dimensional space, OT suffers a few drawbacks: (i) computational
complexity, i.e., super cubic ${\mathcal{O}}(n^{3}\log(n))$ where $n$ is the
number of supports of input measures; (ii) high sample complexity, i.e.,
${\mathcal{O}}(n^{-1/d})$, which requires too many samples to precisely
estimate the OT between two continuous distributions (e.g.,
$\mu_{0}^{\texttt{C}},\mu_{1}^{\texttt{C}}$).

A naive approach to estimate the natural drift map in the multivariate CiC
causal model is to decompose the model for each dimension. Specifically, one
treats the multivariate CiC causal model with $d$-dimensional outcomes ($d>1$)
as $d$ independent univariate CiC causal models, a.k.a., the tensorization of
univariate CiC causal models. It is then efficient to estimate univariate OT
maps for these univariate CiC causal models via their closed-form expressions
of the corresponding univariate OT problems as in the original CiC causal
model (Athey and Imbens, 2006) instead of solving the high-complexity full OT
problem for measures supported in high-dimensional spaces (Torous et al.,
2021). One then tensorizes all the estimated univariate OT maps to create an
estimation of the multivariate OT map. However, this approach might fail to
capture the dependence structure among dimensions of the multivariate OT map,
i.e., the natural drift map f, and thus when one uses this tensorized map to
pushforward samples of $\mu_{0}^{\texttt{T}}$, one might produce a
counterfactual distribution with a wrong dependence structure (e.g., see Fig.
2).

In this work, we propose an efficient approach to leverage the advantages of
the univariate CiC causal model for the multivariate CiC by seeking a robust
latent $1$-dimensional space for OT estimation. More precisely, our approach
is inspired by the subspace robust OT (Paty and Cuturi, 2019) whose authors
proposed to estimate the OT map in low-dimensional subspace to reduce the
sample complexity for the OT problem, and further increase the robustness of
OT estimation with respect to noise.

##  3 The proposed method based on robust OT over latent $1$-dimensional
subspaces

In order to efficiently leverage the advantages of the univariate CiC causal
model and mitigate issues in the naive tensorization for the multivariate CiC,
in this work, we propose to seek a robust latent $1$-dimensional subspace as a
surrogate to estimate the univariate OT map to bridge the univariate and
multivariate CiC causal models. This approach is also known as _max-min robust
variant of OT_ (Paty and Cuturi, 2019).

Let $\mathcal{G}_{1}$ be the Grassmannian of the $1$-dimensional subspaces of
$\mathbb{R}^{d}$, defined as

$\mathcal{G}_{1}=\\{E\subset\mathbb{R}^{d}\mid\text{dim}(E)=1\\},$

where $\text{dim}(E)$ is the dimension of the space $E$. Given two measures
$\mu,\nu$ supported on the space $\mathbb{R}^{d}$, the max-min robust OT (Paty
and Cuturi, 2019) considers the maximal OT distance over all possible
$1$-dimensional projections of input probability measures. Denote $P_{E}$ as
the projector on the $1$-dimensional space $E$, the robust OT is defined as

$\widetilde{\text{ROT}}(\mu,\nu)=\sup_{E\in\mathcal{G}_{1}}\text{OT}({P_{E}}_{\sharp}\mu,{P_{E}}_{\sharp}\nu).$
(8)

One way to parameterize the projection on the $1$-dimensional space $P_{E}$
and the Grassmannian manifold $\mathcal{G}_{1}$ is to utilize a projected
direction vector $\omega$ on the sphere centering at the origin with a radius
$1$ as follows:

$\overline{\text{ROT}}(\mu,\nu)=\max_{\omega\in\mathbb{R}^{d}\,|\,\left\lVert\omega\right\rVert_{2}=1}\text{OT}({P_{\omega}}_{\sharp}\mu,{P_{\omega}}_{\sharp}\nu),$
(9)

where $P_{\omega}$ denotes the projector on the direction $\omega$. For a
support $x\in\mathbb{R}^{d}$, its projection by the projected direction
$\omega$ is computed by $\left<x,\omega\right>$. However, the problem
$\overline{\text{ROT}}$ is non-convex, which is usually approximated by the
first-order method in practical applications (Deshpande et al.,
2019).333$\overline{\text{ROT}}$ is also known as the max-sliced Wasserstein
(Deshpande et al., 2019).

Much as the sliced-Wasserstein (SW) (Rabin et al., 2011)444SW projects
supports into $1$-dimensional space and exploits the closed-form expression of
the univariate OT. which is usually approximated by averaging over a few
random directions in practical applications instead of integrating over all
possible directions on the sphere, and in order to optimize the robust OT
efficiently, we also seek a direction over a subset of projected directions
$\Omega$ for the robust OT, defined as follows:

$\text{ROT}(\mu,\nu)=\max_{\omega\in\Omega}\,\text{OT}({P_{\omega}}_{\sharp}\mu,{P_{\omega}}_{\sharp}\nu),$
(10)

where we construct $\Omega$ by randomly sampling $k$ directions similar as SW.
We observe that a small $k$, e.g., in the range of $10$ to $50$, provides a
good balance between computation speed and accuracy. For estimating the
natural drift map f, $\mu$ and $\nu$ are chosen to be the empirical versions
of $\mu_{0}^{\texttt{C}}$ and $\mu_{1}^{\texttt{C}}$ respectively.

While our approach inherits the fast computation of the CiC estimator for
univariate causal models and thus can scale up for large-scale causal
inference applications, it also gives good performances in estimating the
natural drift map f. This is in line with previous observations from various
applications of sliced-based OT (Paty and Cuturi, 2019; Deshpande et al.,
2019).

######  Remark 3.1.

A technical remark is that it may not exist the optimal transport map $T$ for
the OT problem in Eq. (4) (i.e., the Monge formulation of OT problem). In
practical applications, one can relax such the Monge problem of OT by the
Kantorovich formulation of OT problem (Peyré and Cuturi, 2019), in which the
optimal transport plan always exists.555We give a review about the Kantorovich
formulation of OT problem in Section C. Consequently, one can utilize such
optimal transport plan and leverage the barycentric map (Bonneel et al., 2016)
for the OT map in the Monge problem of OT.

##  4 Synthetic experiments

All experiments were carried out on commodity hardware and can be reproduced
by the code in the supplementary. From here, we will denote the naive
tensorization of the univariate CiC method for estimating multidimensional
counterfactual distributions simply as the CiC method.

###  4.1 Illustrative examples

We demonstrate our method through two 2D examples for ease of visualization.
For the latent distributions $\nu^{\texttt{C}}$ and $\nu^{\texttt{T}}$, we
choose bivariate Gamma distributions in the first example and 2D Gaussian
mixtures in the second example. These settings lead to an unimodal
counterfactual for the Gamma case, and a multimodal counterfactual
distribution for the Gaussian mixture case, as can be seen in the ground truth
panels of Fig. 2. More details on experiment settings can be found in Appendix
A.

Figure 2: Counterfactual distribution estimated by each method. While CiC
failed to capture the dependence structures between dimensions of the
counterfactual distributions in both examples, our method succeeded.

The proposed method produced a counterfactual distribution with the correct
correlation structure, while CiC failed to do so. This demonstrates the
1-dimensional subspace created by our method can indeed provide the type of
information that CiC cannot capture. To systematically evaluate the
performance of each method, we generate $10$ datasets and then measure the
running time as well as the OT distance between the estimated counterfactual
distribution of each method and the empirical version of the true
counterfactual distribution $\mu_{1}^{\texttt{T}^{*}}$. This empirical
distribution, denoted as $\widetilde{\mu_{1}^{\texttt{T}^{*}}}$, is regarded
as the ground truth in our experiments. The averaged values are reported in
Fig. 2.

Figure 3: Averaged running time and OT distance to ground truth. The OT
distance is measured between the estimation result of each method and the
ground truth $\widetilde{\mu_{1}^{\texttt{T}^{*}}}$, which is the empirical
version of the true counterfactual distribution. Values are calculated over
ten datasets. Our method outperforms the OT approach in terms of speed while
being significantly better than the CiC approach in terms of accuracy measured
by OT distance to ground truth.

The CiC estimator has the worst averaged OT distance to ground truth among the
three methods. This might be due to its failure to capture the correlation
structure between dimensions. The proposed method is about one order smaller
in OT distance to ground truth than the CiC while being about two orders
faster than the OT approach.

###  4.2 Varying the number of samples $n$

In this experiment, we check the findings in the illustrative examples by
running the same experiment setting, with $d$ fixed at $2$, for various values
of $n$. For each value of $n$, we generate $10$ datasets and report the
averaged running time as well as the averaged OT distance between the
estimated counterfactual distribution of each method and the ground truth
$\widetilde{\mu_{1}^{\texttt{T}^{*}}}$. We also add Sinkhorn as a baseline.
The results are shown in Fig. 4.

Figure 4: Effect of varying the number of samples $n$ for two types of latent
distributions. In both cases, the proposed method is much faster than both the
OT and Sinkhorn approaches, while being much more accurate than the CiC.

The running time of our method is close to that of the CiC and is faster than
OT and Sinkhorn for all values of $n$. This is in line with the theoretical
worst-case running time of each method. Regarding OT distance to ground truth,
while being worse than OT, our method outperforms CiC for all values of $n$.
This suggests that, while a large $n$ might help CiC in estimating the
marginals of the counterfactual distribution in each dimension, the advantages
of the 1-dimensional subspace constructed by our method do not diminish.

In comparison with Sinkhorn, our method is both faster and more accurate. We
caution that the performance of Sinkhorn depends heavily on the strength of
the entropic regularization term, and choosing a suitable value for the
entropic regularization hyperparameter is non-trivial. However, as stated
earlier, since the worst-case running time of Sinkhorn is quadratic, it is
reasonable to expect it to be generally slower than our method. Additional
results that include Sinkhorn with different hyperparameters are shown in
Appendix B.

###  4.3 Varying the dimension $d$

In order to preserve the computational efficiency of CiC, the robust subspace
in our method has to be one-dimensional. We investigate whether this subspace
can still capture meaningful information in high-dimensional cases, by varying
$d$ while keeping the number of samples $n$ fixed at $5000$. In this
experiment, the latent distribution is multivariate Gamma. For each $d$, we
generate $10$ datasets and report the averaged running time as well as the
averaged OT distance to ground truth. More details on experiment settings can
be found in Appendix A. The results are shown in Fig. 5. Additional results
that include Sinkhorn with different hyperparameters are shown in Appendix B.

Figure 5: Varying the dimension $d$ while fixing number of samples $n=5000$.
The latent distribution is multivariate Gamma. When $d$ is high, our method is
among the bests in both running time and accuracy. This suggests that the
robust 1-dimensional subspace of our method can capture meaningful information
even in high-dimensional situations.

It is interesting to observe that when $d$ is high, our method is the fastest
method, i.e., even faster than CiC. This is in line with the theoretical
worst-case running time of CiC being $\mathcal{O}(dn\log n)$ and ours being
$\mathcal{O}(n\log n+dn)$ when $k$ is fixed. As $d$ increases, the superiority
of OT over our method in terms of OT distance to the ground truth diminishes
and even reverses: our method is as accurate as, if not better than, OT when
$d=100$. Since the performance of CiC is still much worse, this reversal comes
from the robust one-dimensional subspace. This is consistent with previous
observations that found robustifications can alleviate the poor sample
complexity of OT (Paty and Cuturi, 2019).

###  4.4 Varying the number of projections $k$

In all previous experiments, the number of projections used in constructing
the robust 1-dimensional subspace has been fixed at $k=10$. Using the same
setting as the experiment above, for each $d$, we look at one dataset and
inspect how varying $k$ affects the quality of the method. One trade-off to
expect is that increasing $k$ might improve the accuracy of the estimation,
e.g., finding a better subspace or reducing the variance between each run, at
the cost of longer running time. We also investigate a closely related
approach to estimate the counterfactual distribution by using the
$\overline{\text{ROT}}$ objective function in Eq. (9) optimized by Adam
(Kingma and Ba, 2015) with different numbers of iterations. The results are
shown in Fig. 6.

Figure 6: Performances of our proposed method and an approach employed the
$\overline{\text{ROT}}$ function when varying $k$. Here $k$ has two different
meanings depending on which method is considered. In our proposed method that
uses the ROT function in Eq. (10), $k$ is the number of projections used to
construct the robust one-dimensional subspace. In the approach that employs
the $\overline{\text{ROT}}$ function in Eq. (9), $k$ is the number of
iterations of Adam first-order method. Our approach is as good as
$\overline{\text{ROT}}$ in terms of accuracy while being much faster.
Increasing the number of projections in our method to the hundreds marginally
improves the accuracy and variance of the results, while requiring more time.

We found that increasing $k$ to the hundreds indeed improves the variance and
accuracy of the result of our method, as can be seen from the mean and
variance of OT distance to ground truth when $d=100$. However, these
improvements are marginal and come at the cost of a linearly longer running
time. This is the reason we suggest choosing a small value for $k$, namely
around $10-50$, as a reasonable balance region between accuracy and running
time. We also observe that the $\overline{\text{ROT}}$ approach with Adam
offered no improvement in terms of variance and accuracy of the results, while
significantly running longer. This ineffectiveness of Adam might be due to the
non-convexity and non-smooth of $\overline{\text{ROT}}$.

##  5 A real dataset example

We demonstrate the working of our proposed method on the classical data of
Card and Krueger (1993) (CK). On April 1, 1992, New Jersey’s minimum wage rose
from $4.25 to $5.05 per hour, while Pennsylvania’s did not. This provided an
opportunity to estimate the causal impact of the rise on employment in fast-
food restaurants in New Jersey, by analyzing employment data of New Jersey,
i.e., the treatment group, and Pennsylvania, i.e., the control group, before
and after the rise. In CK and subsequent re-examinations (Lu and Rosenbaum,
2004; Card and Krueger, 2000), the number of full-time employees (FT) and the
number of part-time employees (PT) in each restaurant were converted to a
single number, the full-time equivalent employees (FTE), which is defined as
$FTE=FT+0.5\times PT$. The causal impact of the rise on FTE was then
investigated.

This conversion may lose fine details in the characteristics of restaurants.
For the same FTE, there might be a restaurant with high FT and low PT and
another one with low FT and high PT, depending on each restaurant’s
characteristics. These restaurants might respond differently to the increase
in minimum wages. Therefore, analyzing only the univariate FTE risks
confounding those different trends. Simultaneous analysis of FT and PT can be
expected to better capture different trends in the responses of restaurants by
dissecting the causal impacts of increasing minimum wage on FT and PT. The
estimation results of the $2$-dimensional counterfactual distribution of FT
and PT in New Jersey after the rise are shown in Fig. 7. More details on the
dataset can be found in Appendix A.

Figure 7: Counterfactual distribution of the numbers of full-time (FT) and
part-time (PT) employees in restaurants in New Jersey. The data is from Card
and Krueger (1993). The panel of our method shows the result of one typical
run. Our method captures well the correlation between FT and PT in the result
of OT, while CiC struggles in the region of high FT and PT. The OT distance
between the result of CiC and the result of OT is $72.26$. For our method,
this number averaged over $1000$ runs is $68.66\pm 1.42$.

We compare our method with CiC by measuring how close the estimation results
of our method and CiC to the estimation result of OT, which can be regarded as
the standard in terms of accuracy when one does not know the ground truth. A
quick visual inspection of Fig. 7 reveals that our result captures well the
relationship between PT and FT as found in the result of OT, while CiC
struggles in the region where both PT and FT are large. Since there is
randomness in our method, we measure the OT distance between our method and
the result of OT, averaging over $1000$ runs. The OT distance between the
result of CiC and the result of OT is $72.26$, while that of our method is
$68.66\pm 1.42$, where the confidence interval is two standard deviations.
These numbers reinforce the aforementioned visual impressions and offer
statistical evidence to support the conclusion that our method captured better
the relationship between FT and PT than CiC.

##  6 Concluding remarks

We proposed a method for estimating the counterfactual distribution in
multidimensional CiC models. Our method, like CiC, enjoys the computational
efficiency of one-dimensional optimal transports while utilizing correlation
information that is ignored under CiC. Through synthetic and real-dataset
experiments, our method is shown to consistently outperform CiC in terms of
accuracy, while running at a fraction of the time of the multidimensional OT
approach. In future works, we plan to explore the robustness of our proposed
method in the presence of outliers or noises, as well as in other causal
settings, such as the triple difference model (Gruber, 1994; Olden and Møen,
2022).

This work was partially supported by Shiga University Competitive Research
Fund, JST CREST JPMJCR22D2, JPMJCR2015, JST MIRAI program JPMJMI21G2, and ISM
joint-research 2023-ISMCRP-2010.

## References

  * Abadie (2005) Alberto Abadie.  Semiparametric difference-in-differences estimators.  _The Review of Economic Studies_ , 72(1):1–19, 2005.  ISSN 00346527, 1467937X.  URL http://www.jstor.org/stable/3700681.
  * Akbari et al. (2023) Sina Akbari, Luca Ganassali, and Negar Kiyavash.  Learning causal graphs via monotone triangular transport maps, 2023.
  * Angrist and Krueger (1999) Joshua D. Angrist and Alan B. Krueger.  Chapter 23 - empirical strategies in labor economics.  volume 3 of _Handbook of Labor Economics_ , pages 1277–1366. Elsevier, 1999.  https://doi.org/10.1016/S1573-4463(99)03004-7.  URL https://www.sciencedirect.com/science/article/pii/S1573446399030047.
  * Ashman et al. (2023) Matthew Ashman, Chao Ma, Agrin Hilmkil, Joel Jennings, and Cheng Zhang.  Causal reasoning in the presence of latent confounders via neural ADMG learning.  In _The Eleventh International Conference on Learning Representations_ , 2023.  URL https://openreview.net/forum?id=dcN0CaXQhT.
  * Athey and Imbens (2006) Susan Athey and Guido W. Imbens.  Identification and inference in nonlinear difference-in-differences models.  _Econometrica_ , 74(2):431–497, 2006.  ISSN 00129682, 14680262.  URL http://www.jstor.org/stable/3598807.
  * Ben-Tal et al. (2009) Aharon Ben-Tal, Laurent El Ghaoui, and Arkadi Nemirovski.  _Robust optimization_ , volume 28.  Princeton university press, 2009.
  * Bertsimas et al. (2011) Dimitris Bertsimas, David B Brown, and Constantine Caramanis.  Theory and applications of robust optimization.  _SIAM review_ , 53(3):464–501, 2011.
  * Blundell and Costa Dias (2009) Richard Blundell and Monica Costa Dias.  Alternative approaches to evaluation in empirical microeconomics.  _Journal of Human Resources_ , 44(3), 2009.  URL https://EconPapers.repec.org/RePEc:uwp:jhriss:v:44:y:2009:i3:p565-640.
  * Blundell and Macurdy (1999) Richard Blundell and Thomas Macurdy.  Chapter 27 - labor supply: A review of alternative approaches.  volume 3 of _Handbook of Labor Economics_ , pages 1559–1695. Elsevier, 1999.  https://doi.org/10.1016/S1573-4463(99)03008-4.  URL https://www.sciencedirect.com/science/article/pii/S1573446399030084.
  * Bonhomme and Sauder (2011) Stéphane Bonhomme and Ulrich Sauder.  Recovering distributions in difference-in-differences models: A comparison of selective and comprehensive schooling.  _The Review of Economics and Statistics_ , 93(2):479–494, 2011.  ISSN 00346535, 15309142.  URL http://www.jstor.org/stable/23015949.
  * Bonneel et al. (2016) Nicolas Bonneel, Gabriel Peyré, and Marco Cuturi.  Wasserstein barycentric coordinates: histogram regression using optimal transport.  _ACM Trans. Graph._ , 35(4):71–1, 2016.
  * Brenier (1991) Yann Brenier.  Polar factorization and monotone rearrangement of vector-valued functions.  _Communications on pure and applied mathematics_ , 44(4):375–417, 1991.
  * Callaway and Li (2019) Brantly Callaway and Tong Li.  Quantile treatment effects in difference in differences models with panel data.  _Quantitative Economics_ , 10(4):1579–1618, 2019\.  https://doi.org/10.3982/QE935.  URL https://onlinelibrary.wiley.com/doi/abs/10.3982/QE935.
  * Card and Krueger (1993) David Card and Alan B Krueger.  Minimum wages and employment: A case study of the fast food industry in new jersey and pennsylvania.  Working Paper 4509, National Bureau of Economic Research, October 1993\.  URL http://www.nber.org/papers/w4509.
  * Card and Krueger (2000) David Card and Alan B. Krueger.  Minimum wages and employment: A case study of the fast-food industry in new jersey and pennsylvania: Reply.  _The American Economic Review_ , 90(5):1397–1420, 2000.  ISSN 00028282.  URL http://www.jstor.org/stable/2677856.
  * Cuturi (2013) M. Cuturi.  Sinkhorn distances: Lightspeed computation of optimal transport.  In _Advances in Neural Information Processing Systems_ , pages 2292–2300, 2013.
  * Deshpande et al. (2019) Ishan Deshpande, Yuan-Ting Hu, Ruoyu Sun, Ayis Pyrros, Nasir Siddiqui, Sanmi Koyejo, Zhizhen Zhao, David Forsyth, and Alexander G. Schwing.  Max-sliced wasserstein distance and its use for gans.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2019.
  * Dhouib et al. (2020) Sofien Dhouib, Ievgen Redko, Tanguy Kerdoncuff, Rémi Emonet, and Marc Sebban.  A swiss army knife for minimax optimal transport.  In _International Conference on Machine Learning_ , pages 2504–2513. PMLR, 2020.
  * Donald and Lang (2007) Stephen G Donald and Kevin Lang.  Inference with Difference-in-Differences and Other Panel Data.  _The Review of Economics and Statistics_ , 89(2):221–233, 05 2007.  ISSN 0034-6535.  10.1162/rest.89.2.221.  URL https://doi.org/10.1162/rest.89.2.221.
  * Gangbo (1999) Wilfrid Gangbo.  The monge mass transfer problem and its applications.  _Contemporary Mathematics_ , 226:79–104, 1999.
  * Gangl (2010) Markus Gangl.  Causal inference in sociological research.  _Annual Review of Sociology_ , 36(1):21–47, 2010\.  10.1146/annurev.soc.012809.102702.  URL https://doi.org/10.1146/annurev.soc.012809.102702.
  * Glass et al. (2013) Thomas A. Glass, Steven N. Goodman, Miguel A. Hernán, and Jonathan M. Samet.  Causal inference in public health.  _Annual Review of Public Health_ , 34(1):61–75, 2013.  10.1146/annurev-publhealth-031811-124606.  URL https://doi.org/10.1146/annurev-publhealth-031811-124606.  PMID: 23297653.
  * Gruber (1994) Jonathan Gruber.  The incidence of mandated maternity benefits.  _The American Economic Review_ , 84(3):622–641, 1994.  ISSN 00028282.  URL http://www.jstor.org/stable/2118071.
  * Hwang et al. (2023) Inwoo Hwang, Yunhyeok Kwak, Yeon-Ji Song, Byoung-Tak Zhang, and Sanghack Lee.  On discovery of local independence over continuous variables via neural contextual decomposition.  In _2nd Conference on Causal Learning and Reasoning_ , 2023.  URL https://openreview.net/forum?id=-aFd28Uy9td.
  * Immer et al. (2023) Alexander Immer, Christoph Schultheiss, Julia E Vogt, Bernhard Schölkopf, Peter Bühlmann, and Alexander Marx.  On the identifiability and estimation of causal location-scale noise models.  In _Proceedings of the 40th International Conference on Machine Learning_ , ICML’23. JMLR.org, 2023.
  * Keel et al. (1988) Lee H Keel, SP Bhattacharyya, and Jo W Howze.  Robust control with structure perturbations.  _IEEE Transactions on Automatic Control_ , 33(1):68–78, 1988.
  * Kingma and Ba (2015) Diederik Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In _International Conference on Learning Representations (ICLR)_ , San Diega, CA, USA, 2015.
  * Kladny et al. (2023) Klaus-Rudolf Kladny, Julius von Kügelgen, Bernhard Schölkopf, and Michael Muehlebach.  Deep backtracking counterfactuals for causally compliant explanations, 2023.
  * Kleinberg and Hripcsak (2011) Samantha Kleinberg and George Hripcsak.  A review of causal inference for biomedical informatics.  _Journal of Biomedical Informatics_ , 44(6):1102–1112, 2011.  ISSN 1532-0464.  https://doi.org/10.1016/j.jbi.2011.07.001.  URL https://www.sciencedirect.com/science/article/pii/S1532046411001195.
  * Lechner (2011) Michael Lechner.  The estimation of causal effects by difference-in-difference methods.  _Foundations and Trends® in Econometrics_ , 4(3):165–224, 2011.  ISSN 1551-3076.  10.1561/0800000014.  URL http://dx.doi.org/10.1561/0800000014.
  * Lu and Rosenbaum (2004) Bo Lu and Paul R Rosenbaum.  Optimal pair matching with two control groups.  _Journal of Computational and Graphical Statistics_ , 13(2):422–434, 2004.  10.1198/1061860043470.  URL https://doi.org/10.1198/1061860043470.
  * McCann (1995) Robert J McCann.  Existence and uniqueness of monotone measure-preserving maps.  1995\.
  * Morimoto and Doya (2001) Jun Morimoto and Kenji Doya.  Robust reinforcement learning.  _Advances in neural information processing systems_ , pages 1061–1067, 2001.
  * Olden and Møen (2022) Andreas Olden and Jarle Møen.  The triple difference estimator.  _The Econometrics Journal_ , 25(3):531–553, 03 2022.  ISSN 1368-4221.  10.1093/ectj/utac010.  URL https://doi.org/10.1093/ectj/utac010.
  * Panaganti et al. (2022) Kishan Panaganti, Zaiyan Xu, Dileep Kalathil, and Mohammad Ghavamzadeh.  Robust reinforcement learning using offline data.  In _Advances in Neural Information Processing Systems_ , volume 35, pages 32211–32224, 2022.
  * Paty and Cuturi (2019) François-Pierre Paty and Marco Cuturi.  Subspace robust Wasserstein distances.  In _Proceedings of the 36th International Conference on Machine Learning_ , pages 5072–5081, 2019.
  * Pearl (2009) Judea Pearl.  Causal inference in statistics: An overview.  _Statistics Surveys_ , 3(none):96 – 146, 2009\.  10.1214/09-SS057.  URL https://doi.org/10.1214/09-SS057.
  * Peyré and Cuturi (2019) Gabriel Peyré and Marco Cuturi.  Computational optimal transport.  _Foundations and Trends® in Machine Learning_ , 11(5-6):355–607, 2019.
  * Rabin et al. (2011) Julien Rabin, Gabriel Peyré, Julie Delon, and Marc Bernot.  Wasserstein barycenter and its application to texture mixing.  In _International Conference on Scale Space and Variational Methods in Computer Vision_ , pages 435–446, 2011.
  * Roth and Sant’Anna (2023) Jonathan Roth and Pedro H. C. Sant’Anna.  When is parallel trends sensitive to functional form?  _Econometrica_ , 91(2):737–747, 2023.  https://doi.org/10.3982/ECTA19402.  URL https://onlinelibrary.wiley.com/doi/abs/10.3982/ECTA19402.
  * Roth et al. (2022) Jonathan Roth, Pedro H. C. Sant’Anna, Alyssa Bilinski, and John Poe.  What’s Trending in Difference-in-Differences? A Synthesis of the Recent Econometrics Literature.  Papers 2201.01194, arXiv.org, January 2022.  URL https://ideas.repec.org/p/arx/papers/2201.01194.html.
  * Santambrogio (2015) Filippo Santambrogio.  _Optimal transport for applied mathematicians_.  Birkäuser, 2015.
  * Sauter et al. (2023) Andreas W.M. Sauter, Erman Acar, and Vincent Francois-Lavet.  A meta-reinforcement learning algorithm for causal discovery.  In _2nd Conference on Causal Learning and Reasoning_ , 2023.  URL https://openreview.net/forum?id=p6NnDqJM_jL.
  * Schölkopf (2022) Bernhard Schölkopf.  _Causality for Machine Learning_ , page 765–804.  Association for Computing Machinery, New York, NY, USA, 1 edition, 2022\.  ISBN 9781450395861.  URL https://doi.org/10.1145/3501714.3501755.
  * Snow (1854) John Snow.  The cholera near golden square, and at deptford.  _Medical Times and Gazette_ , 9:321–322, 1854.
  * Snow (1855) John Snow.  On the mode of communication of cholera.  page 162.1, 1855.
  * Sofer et al. (2016) Tamar Sofer, David B. Richardson, Elena Colicino, Joel Schwartz, and Eric J. Tchetgen Tchetgen.  On Negative Outcome Control of Unobserved Confounding as a Generalization of Difference-in-Differences.  _Statistical Science_ , 31(3):348 – 361, 2016\.  10.1214/16-STS558.  URL https://doi.org/10.1214/16-STS558.
  * Torous et al. (2021) William Torous, Florian Gunsilius, and Philippe Rigollet.  An optimal transport approach to causal inference.  2021\.  10.48550/ARXIV.2108.05858.  URL https://arxiv.org/abs/2108.05858.
  * Tu et al. (2022) Ruibo Tu, Hedvig Kjellstrom, Kun Zhang, and Cheng Zhang.  Optimal transport for causal discovery.  In _ICLR 2022_ , April 2022.  URL https://www.microsoft.com/en-us/research/publication/optimal-transport-for-causal-discovery-2/.
  * Villani (2003) Cédric Villani.  _Topics in Optimal Transportation_.  American Mathematical Society, 2003.
  * Villani (2008) Cédric Villani.  _Optimal transport: Old and New_.  Springer, 2008.
  * Wildberger et al. (2023) Jonas Bernhard Wildberger, Siyuan Guo, Arnab Bhattacharyya, and Bernhard Schölkopf.  On the interventional kullback-leibler divergence.  In Mihaela van der Schaar, Cheng Zhang, and Dominik Janzing, editors, _Proceedings of the Second Conference on Causal Learning and Reasoning_ , volume 213 of _Proceedings of Machine Learning Research_ , pages 328–349. PMLR, 11–14 Apr 2023.  URL https://proceedings.mlr.press/v213/wildberger23a.html.
  * Xu et al. (2009) Huan Xu, Constantine Caramanis, and Shie Mannor.  Robustness and regularization of support vector machines.  _Journal of machine learning research_ , 10(7), 2009.

##  Appendix A Details on experiment settings

###  A.1 The illustrative examples

The latent distributions in two examples are as follows. In the bivariate
Gamma example, the first dimension of $\nu^{\texttt{C}}$ is a Gamma
distribution with shape $2$ and scale $3$, while the second dimension is Gamma
with shape $3$ and scale $2$. For $\nu^{\texttt{T}}$, the first and second
dimensions are reversed. In the Gaussian mixture example, the first and second
dimensions of $\nu^{\texttt{C}}$ are $0.5\mathcal{N}(1,1)+0.5\mathcal{N}(5,1)$
and $0.5\mathcal{N}(2,1)+0.5\mathcal{N}(4,1)$, respectively. For
$\nu^{\texttt{T}}$, the first and second dimensions are reversed.

The production functions $h_{0}$ and $h_{1}$ are
$h_{0}(u)=\big{(}\begin{smallmatrix}1&0.5\\\ 0.5&1\end{smallmatrix}\big{)}u$
and $h_{1}(u)=\big{(}\begin{smallmatrix}1&-0.5\\\
-0.5&1\end{smallmatrix}\big{)}u$ for a length-$2$ vector $u$. These functions
are co-monotone, and thus the natural drift map f is identifiable (Torous et
al., 2021).

###  A.2 Settings for the experiments with varying $d$

We discuss the settings for the production functions $h_{0}$ and $h_{1}$. When
$d\geq 2$, for the natural drift map f to be identifiable, the functions
$h_{0}$ and $h_{1}$ need to be co-monotone (Torous et al., 2021), i.e.,

$\langle h_{0}(x)-h_{0}(y),h_{1}(x)-h_{1}(y)\rangle\geq 0,\quad\forall
x,y\in\mathbb{R}^{d}.$

When $h_{0}(x)=\mathbf{H}_{0}x$ and $h_{1}(x)=\mathbf{H}_{1}x$ for
$\mathbb{R}^{d\times d}$ matrices $\mathbf{H}_{0}$ and $\mathbf{H}_{1}$, this
condition is satisfied if $\mathbf{H}_{0}^{T}\mathbf{H}_{1}$ is positive semi-
definite. We generate one matrix $\mathbf{H}_{0}$ as a $d\times d$ matrix
where each off-diagonal entry is uniformly distributed in $(0,1)$ and the
diagonal entries are $1$. Note that $\mathbf{H}_{0}$ generated this way is
almost surely invertible. We then generate a diagonal matrix $\mathbf{B}$
where each diagonal entry is uniformly distributed in $(0,1)$. We then let
$\mathbf{H}_{1}=(\mathbf{H}_{0}^{-1})^{T}\mathbf{B}$. This will ensure that
$\mathbf{H}_{0}^{T}\mathbf{H}_{1}$ is equal to $\mathbf{B}$ and thus positive
semi-definite. We then fix the pair ($\mathbf{H}_{0}$, $\mathbf{H}_{1}$) and
then generate datasets using this pair.

The latent distributions are as follows. For the control group, each dimension
independently follows a Gamma distribution with shape $2$ and scale $3$. For
the treatment group, each dimension independently follows a Gamma distribution
with shape $3$ and scale $2$.

###  A.3 CK data

The dataset is downloaded from
https://davidcard.berkeley.edu/data_sets/njmin.zip. We also include the
following covariates into the analysis: HRSOPEN, OPEN, NMGRS, NREGS, INCTIME,
PSODA, and PENTREE, and estimate the $9$-dimensional counterfactual
distribution. We process the data by removing samples that contain missing
values at any covariates. The final numbers of samples after this pre-
processing are $57$ for the control group and $220$ for the treatment group.

##  Appendix B Additional results with different hyperparameters of Sinkhorn

In Sinkhorn, an entropic regularization term $\lambda\Omega(T)$, where
$\Omega(T)$ is the entropy of the transportation plan $T$, is added to the
objective function of OT. The hyperparameter $\lambda\geq 0$ controls the
strength of the regularization. The results reported in Figs. 4 and 5 are
obtained using $\lambda=30$. We report in the following figures two more cases
when $\lambda=10$ and $\lambda=90$.

Figure B.1: Additional results of Sinkhorn for the case of varying the number
of samples $n$ while fixing the dimension $d=2$. Figure B.2: Additional
results of Sinkhorn for the case of varying the dimension $d$ while fixing
number of samples $n=5000$.

##  Appendix C Brief review

In this section, we further give brief review for some technical details on
optimal transport (OT) which are used in our work.

#### Kantorovich formulation of OT.

Given two probability distributions $\mu,\nu$ with a cost function $c$, the
Kantorovich formulation of OT is as follow:

$\text{OT}(\mu,\nu)=\inf_{\pi\in\Pi(\mu,\nu)}\int_{\mathbb{R}^{d}\times\mathbb{R}^{d}}\pi(x,y)c(x,y)d\mu(x)d\nu(y),$

where $\pi$ is known as the transport plan, and $\Pi(\mu,\nu)$ is the set of
all probability distributions on the product space
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ such that its first and second marginals
equal to $\mu,\nu$ respectively.

#### Optimal condition for OT with $d=1$.

If $\mu,\nu$ are two probability measures supported on $\mathbb{R}$, and $\mu$
is atomless (i.e., $\mu$ is absolutely continuous with respect to the Lebesgue
measure), then there exists at least a transport map $T$ such that
$T_{\sharp}\mu=\nu$ (Santambrogio, 2015, Lemma 1.27). With quadratic cost, the
transport map will be the derivative of a convex function, i.e., a
nondecreasing map (see (Santambrogio, 2015, Remark 1.23) and (Gangbo, 1999,
§2) for further details).

#### Optimal condition for OT with $d>1$.

It is also known as Brenier theorem (Brenier, 1991).

Given two probability measures $\mu,\nu$ supported on the $\mathbb{R}^{d}$
space such that $\mu$ is absolutely continuous with respect to the Lebesgue
measure. Then, for all possible map
$\widetilde{T}:\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$ such that
$\widetilde{T}_{\sharp}\mu=\nu$, there is an unique Brenier map $T$ which is
the gradient of a convex function, and $T$ is the optimal transport in the
following sense: the Kantorich formulation of OT between $\mu$ and $\nu$
admits a unique optimal transportation plan $\pi^{*}$ such that
$(x,y)\sim\pi^{*}$ if and only if $x\sim\mu$ and $y=T(x)$, $\mu$-a.s. For
further details, please see (Santambrogio, 2015, Theorem 1.17) and (Brenier,
1991; McCann, 1995).

##  Appendix D Further discussions

#### About the max-min robust OT.

Max-sliced Wasserstein (Deshpande et al., 2019) is similar to our proposed
method, i.e., they consider the sphere for the uncertainty set $\Omega$ of
projections. However, when $\Omega$ is a sphere, the problem becomes non-
convex and non-smooth. The entropic regularized OT (Cuturi, 2013) is a popular
approach to reduce the computation of OT into quadratic computation
complexity, but it comes with a trade-off for a dense optimal transport plan.

#### About the causal model.

The causal model described in this paper is the same as that of Athey and
Imbens (2006) when $d=1$ and that of Torous et al. (2021) for a general $d$.
When $d=1$, it is a non-linear extension of the classical DiD model (Angrist
and Krueger, 1999; Blundell and Macurdy, 1999). Another notable non-linear
extension of this DiD model is provided in Abadie (2005). Some stronger
assumptions to model the natural drift have been proposed in, for example,
Callaway and Li (2019), Roth and Sant’Anna (2023), and Bonhomme and Sauder
(2011).

#### About the use of machine learning methods in causal inference and causal
discovery

Recently, Tu et al. (2022) provided the first method that uses optimal
transport for solving causal discovery, i.e., the task of uncovering the graph
that represents the dependence relationships between variables. Akbari et al.
(2023) generalized the method to higher dimensions and different noise
settings. Other machine learning tools and frameworks have also been used in
causal inference/causal discovery, such as neural networks (Hwang et al.,
2023; Ashman et al., 2023; Immer et al., 2023; Kladny et al., 2023), the
Kullback-Leibler divergence (Wildberger et al., 2023), and reinforcement
learning and meta learning (Sauter et al., 2023).


# Hyper-Auxeticity and the Volume Phase Transition of Polymer Gels

Andrea Ninarello<EMAIL_ADDRESS>CNR Institute of Complex Systems,
Uos Sapienza, Piazzale Aldo Moro 2, 00185, Roma, Italy  Department of Physics,
Sapienza University of Rome, Piazzale Aldo Moro 2, 00185 Roma, Italy  Emanuela
Zaccarelli<EMAIL_ADDRESS>CNR Institute of Complex Systems, Uos
Sapienza, Piazzale Aldo Moro 2, 00185, Roma, Italy  Department of Physics,
Sapienza University of Rome, Piazzale Aldo Moro 2, 00185 Roma, Italy

###### Abstract

Thermoresponsive hydrogels exhibit reversible deswelling at the Volume Phase
Transition (VPT), associated to a minimum of the Poisson’s ratio $\nu$. Recent
numerical investigations uncovered the occurrence of a Hyper-Auxetic
Transition (HAT) ($\nu=-1$) for low-crosslinked hydrogels at low temperature,
accompanied by a critical-like behavior. Here, we perform extensive numerical
simulations to unveil the relation between these two transitions. We find that
the HAT occurs at different temperatures $T$ up to a maximum value of the
crosslinker concentration $c$, thus being clearly distinct from the VPT,
taking place at fixed $T$ and for all $c$. Our results provide novel
fundamental insights on the interplay between network collapse and mechanical
properties of these fascinating materials.

Chemical hydrogels are covalently crosslinked polymer chains dispersed in a
water-based solvent, offering versatility in synthesis and functionality,
driving both scientific inquiry and practical applications [1, 2].
Thermoresponsive hydrogels, like those based on Poly(N-isopropylacrylamide)
(pNIPAM), have gathered significant attention for their reversible swelling
behavior induced by a change of temperature $T$. This is due to a variation of
the affinity between polymer and solvent, usually water, which decreases,
becoming more hydrophobic, upon increasing $T$ [3]. Experimental studies on
neutral and charged pNIPAM hydrogels have elucidated this phenomenon,
revealing a Volume Phase Transition (VPT) and its theoretical underpinnings
within the Flory-Rehner theory of swelling [4, 5].
Mechanical instabilities are found to accompany the transition [6, 7, 8],
showcasing intriguing elastic properties and the onset of auxetic response,
where materials expand perpendicularly with respect to the applied strain,
yielding a negative Poisson’s ratio that is a signature of auxetic behavior
[9, 10, 11]. Auxetic materials often captivated the interest of researchers
for their theoretical implications and practical applications, arising either
from material topology or from the proximity to a thermodynamic transition
[12, 13, 14, 15, 16, 17, 18, 19, 20].

In general, when approaching the lowest limit of mechanically stable
unconstrained solids ($\nu=-1$), a condition that we name as hyper-auxeticity,
the bulk modulus tends to zero, suggesting the presence of very large volume
fluctuations from a thermodynamic viewpoint [15]. However, according to the
mean field theory of polymer networks, the relationship $K+4/3G=0$ is expected
to hold at the VPT [21, 22], implying a positive $\nu$ within linear
elasticity theory, as $\nu=(3K-2G)/[2(3K+G)]$. This prediction is however
inconsistent with experimental observations. While simulations could aid in
reconciling this disparity, so far most computational studies concerning the
VPT of hydrogels have overlooked the interactions between polymers and
solvents and the link between thermodynamic and mechanical response [23, 24,
25, 26].

Recent numerical works [27, 28] have uncovered hyper-auxetic behavior of
polymer networks in good solvent, hence purely repulsive systems. This
condition is met under a slight tension and by reduction of the crosslinker
concentration $c$, leading to the occurrence of a thermodynamic criticality,
that is referred as Hyper-Auxetic Transition (HAT) [28]. However, in that
case, the role played by the attraction and the influence of the deswelling
phenomenon on the elastic properties remained unexplored.

In the present work, we fill this gap by performing extensive computer
simulations to investigate the occurrence of the HAT in polymer networks
experiencing temperature-induced shrinking in a wide range of phase diagram
and for several values of the crosslinker concentration. To this aim, we rely
on a widespread bead-spring model of polymer networks [29, 30] and we model
implicitly the presence of the solvent via a solvophobic interaction [3], that
was previously shown to quantitatively reproduce the features of the networks
across the VPT [32].

Our findings reveal the presence of an attraction-induced HAT consistent with
the Ising universality class, which vanishes with increasing $c$. In addition,
we are able to clearly distinguish the solvent-mediated network collapse
driven by the VPT and the hyper-auxetic criticality, which at low enough $c$
takes place in different regions of the phase diagram. This is in contrast to
the VPT, which always occurs at a fixed temperature, that is dictated by the
underlying monomer-solvent affinity. Our study thus establishes a clear
relationship between thermodynamic and mechanical properties of polymer
networks, paving the way to their fine control in experiments.

We simulate hydrogels using the well-established Kremer-Grest model for
monomer-monomer interactions, where beads of diameter $\sigma$ experience a
Weeks-Chandler-Andersen (WCA) soft repulsion and are linked by finitely
extensible nonlinear elastic (FENE) springs [29, 30]. The network comprises
polymer strands connected by four-valence crosslinkers. The majority of
initial configurations are prepared by adding monomers at the equilibrium bond
length in a diamond-like crystalline structure (ordered systems), while for
one case we also perform self-assembly of patchy particles (disordered
systems). This latter method has been recently envisaged to create low-density
disordered polymer networks [4, 34, 35]. We implicitly account for the solvent
by introducing an attractive term between all polymer beads in the network,
known as the solvophobic potential [3], where the strength of the monomer-
monomer attractive interactions is controlled by a parameter $\alpha$. We
anticipate our findings to align with those of explicit solvent simulations,
as demonstrated previously with microgel particles [36].

We perform NPT molecular dynamics simulations of the hydrogels at different
$c$, $P$ and $\alpha$ values using the LAMMPS package [37] either in
unperturbed conditions or under uniaxial strain. We then calculate the bulk
modulus directly from volume fluctuations as
$K=\frac{k_{b}T<V>}{<V^{2}>-<V>^{2}}$, while the Young modulus and the
Poisson’s ratio are computed though their definitions in strain-stress
simulations ($Y=\frac{\sigma_{\parallel}}{\lambda_{\parallel}}$,
$\nu=-\frac{\lambda_{\perp}}{\lambda_{\parallel}}$), where
$\sigma_{\parallel,\perp}$ and $\lambda_{\parallel,\perp}$ are respectively
the stress and the strain either parallel or perpendicular with respect to the
deformation axis. Strain-stress simulations consist in uniaxial deformations
of different amounts, followed by equilibrium simulations spanning $10^{6}$
steps. We perform deformations in each of the three spatial dimensions in
separate runs and then average the results. The time unit is deﬁned as
$\tau=\sqrt{m\sigma^{2}/\epsilon}$, where $m$ is the monomer mass and
$\epsilon$ sets the energy scale. In addition to the hydrogels, we also
perform NVT simulations of chains of various lengths interacting with the same
potential, in order to unveil the role of the VPT. Further details about the
simulations can be found in the SI.

Figure 1: Hyper-auxeticity in solvent-mediated networks: (a) Equation of state
representing the negative pressure $P$ as a function of $\rho$ for different
attraction strengths $\alpha$ as indicated in the legend for the $c=1\%$
system. (b) Poisson’s ratio $\nu$ as a function of the pressure for the
$c=1\%$ (blue circles) and $1\%_{\text{diso}}$ (orange squares) systems
respectively at $\alpha=0.5$ and $0.425$. (c) Normalized probability
distribution $\Pi(M)$ of the order parameter $M=\rho+se_{X}$, using either
total energy $e_{T}$ (blue line) or non-bonded particles energy $e_{nb}$
(orange line). Curves are rescaled to zero-mean and unitary variance. Results
are for the $c=1\%$ diamond network at $\alpha_{HAT}=0.5064$, reweighted to
$P=9.319\times 10^{-4}$, and using the mixing parameter $s=-3.2,-0.8$ for the
two energies, respectively.

We start by investigating the onset of auxetic behavior in the presence of
attraction. For this purpose, we analyze the equation of state of the ordered
system with $1\%$ crosslinking, varying tension (or negative pressure $-P$)
for different temperatures (quantified by the attraction strength $\alpha$).
The corresponding findings are reported in Fig. 1(a). Notably, as the
effective attraction between particles increases, the change in the density of
the network becomes more abrupt, leading to a discontinuous transition between
two states, characterized by a low and a high density, respectively, when
$\alpha>0.5$. Given the resemblance to a classical liquid-gas transition, we
examine in more detail the system behavior at $\alpha=0.5$ for different
pressures.
Our analysis of the elastic moduli reveals a non-monotonic behavior of the
bulk $K$ and Young $Y$ moduli as pressure changes (see SI) that reflects the
behavior of the Poisson’s ratio $\nu$, reported in Fig.1(b). First, $\nu$
decreases with decreasing $P$, reaching negative values, and then it increases
again at larger tensions. Notably, these two behaviors are distinguished by a
minimum where $\nu=-1$, marking the mechanical stability limit [38, 39]. We
identify this condition as a Hyper-Auxetic point, that in addition to the
mechanical instability, is characterized by huge density fluctuations,
reminiscent of critical-like behavior. To quantify this behavior in the
hypothesis of the occurrence of a second-order critical point, we follow the
typical approach used to study gas-liquid phase separation in atomic and
molecular systems [40, 41]. We thus monitor the behavior of the characteristic
order parameter $M$, which combines density and energy of the system, akin to
approaches used in liquid-gas transitions [40]. Hence, $M=\rho+se_{X}$, where
$s$ is a mixing parameter and $X$ identifies the type of energy employed in
the analysis. In particular, we either employ the total energy $e_{t}$ or the
energy of the non-bonded particles $e_{nb}$ only, as in our previous studies
[27, 28]. By examining the probability distribution of this order parameter
rescaled to zero-mean and unitary variance $\Pi(M)$, we compare it to the
theoretical Ising behavior and find a good agreement for the critical
temperature $\alpha_{HAT}^{1\%}=0.5064$, between the predicted and calculated
histograms with both observables, as shown in Fig.1(c). We ascribe the
differences between the results obtained with either $e_{t}$ or $e_{nb}$ to
the intrinsic difficulties in finding these distributions in computational
studies and to the weak correlation between energy and density for the
considered low-density polymer network. Our results thus suggest that the HAT
transition observed in the presence of attraction is of the Ising universality
class, deepening the analogy to a gas-liquid phase separation.

To provide robustness to our findings, we perform a similar analysis also for
the disordered $c=1\%$ network ($1\%_{\text{diso}}$). We thus evaluate the
equation of state (not shown), establishing the critical-like temperature
$\alpha_{HAT}^{1\%_{\text{diso}}}=0.425$, that is notably distinct from the
corresponding diamond network one. We then compute the elastic properties,
confirming the occurrence of the HAT where $\nu=-1$ at an intermediate
pressure, as also shown in Fig.1(b). Similarly to our previous studies at
$\alpha=0$ [27, 28], we thus confirm that the topology of the network does not
play a role for the existence of the HAT, ruling out a possible geometric
origin.

We now aim to compare these results with previous measurements on hydrogels,
reporting the occurrence of auxetic behavior [9]. Indeed, at first sight, it
may appear that the Hyper-Auxetic Transition (HAT) is merely an amplification
of the negative Poisson’s ratio observed in experiments [9]. However, a more
quantitative analysis is needed to verify this assumption, specifically
examining the relationship between the two transitions: the HAT and the VPT.
To address this issue, we report the hydrogel swelling curve for the $c=1\%$
diamond network in Fig. 2, where the cubic root of volume $V^{1/3}$, analogous
to an effective size of the hydrogel, is plotted as a function of $\alpha$,
which acts as an effective temperature, for different studied pressures.
First, we focus on $P=0$, that is the case representative of microgels in
experiments, and observe that, as usual, the effective size smoothly decreases
with increasing $\alpha$, showing an inflection point at the VPT, namely at
$\alpha_{VPT}\simeq 0.65$ [4]. However, upon decreasing $P$, we find that
deswelling happens in two steps. At first, there is a decrease which
progressively hits the HAT point at $\alpha_{HAT}\sim 0.5$, above which
$V^{1/3}$ becomes discontinuous. Representative snapshots of the hydrogel
around the HAT are also reported in Fig. 2, showing the spinodal-like
character of the transition. Indeed, the high density state is far from
homogeneous and a further compaction only happens at larger $\alpha$. This
second deswelling step finally recovers the VPT behavior at high temperatures,
where for all pressures the deswelling curves are virtually indistinguishable
from one another.

These results clearly indicate that the HAT driving the auxetic response of
the system is a distinct phenomenon from the VPT and we now aim to unveil
their underlying nature and mutual interplay in more detail.

Figure 2: Top: Swelling curve for the $c=1\%$ ordered network. Cubic root of
the volume as a function of $\alpha$ for different $P$ as indicated in the
legend. The vertical region indicates the location of the VPT. Bottom:
Snapshots of the $c=1\%$ system at $P=-9.05\times 10^{-4}$ around the HAT, for
$\alpha=0.48$ (left) and $0.56$ (right), respectively. Monomers and
crosslinkers are shown in violet and dark blue, respectively.

Figure 3: HAT and VPT. (a) Poisson’s ratio as a function of $\alpha$ for
ordered systems with $c=1,2,3,5,7.5\%$ and corresponding $P=9.25\times
10^{-4},3.5\times 10^{-3},7.9\times 10^{-3},1.8\times 10^{-2},4.25\times
10^{-2}$. The horizontal dashed black line indicates the hyper-auxetic
threshold, when $\nu_{min}=-1$, while the vertical box marks $\alpha_{c}$,
i.e., the occurrence of the VPT. (b) Top. The phase separation effective
temperature of chains $\alpha_{c}$ is reported as a function of chain length
$l$. Bottom. Effective temperature, corresponding to the minimum in $\nu$ as a
function of $c$. Filled circles indicate HAT points where $\nu=-1$, while
stars represent points where $\nu\gtrsim-1$. Both top and bottom curves tend
to the VPT one, $\alpha_{VPT}$, as $c$ and $l$ increase, as signaled by the
horizontal line.

To this aim, we first examine the occurrence of the HAT for various
crosslinker concentrations. We explore diamond networks with $c=0.5,1,2,3,5,$
and $7.5\%$ and investigate their elastic behavior in order to identify the
conditions corresponding to the minimum Poisson’s ratio $\nu_{min}$. We first
retrieve the corresponding pressure $P_{min}$, that, as expected, increases in
absolute value for stiffer networks. Then, we calculate $\nu$ as a function of
$\alpha$ for each network in order to determine whether the system is able to
reach $\nu_{min}=-1$ independently of crosslinker concentration. The results
of this extensive analysis are summarized in Fig. 3(a), where only some of the
investigated $c$ values are reported to improve visualization. Our findings
reveal the occurrence of hyper-auxeticity at a given $\alpha_{HAT}$,
associated to $\nu_{min}\approx-1$, up to crosslinker concentrations $c\simeq
3\%$, while networks with higher $c$ are still found to exhibit a minimum in
$\nu$, at a characteristic temperature that we call $\alpha_{min}$, which
however does not reach the HAT condition any longer, since $\nu_{min}>-1$.
Indeed, in analogy with previous results obtained for $\alpha=0$ [27] and for
microgels [42], we find that $\nu_{min}$ increases with the stiffness of the
network, approaching the HAT condition only at sufficiently low $c$. It is now
interesting to look at the $\alpha_{HAT}$ and $\alpha_{min}$ dependence on
crosslinker concentration, which is reported in Fig. 3(b). It is evident that
$\alpha_{HAT}$ continuously turn to $\alpha_{min}$, progressively increasing
with the system connectivity and finally tending to $\alpha_{VPT}\sim 0.65$ at
high $c$. Hence, while the HAT terminates at a finite $c$ value, we provide
evidence that its echo, characterized by the presence of a minimum in $\nu$ at
even larger $c$, converges to the temperature of the VPT as crosslinker
concentration increases. This is why a minimum in $\nu$ is detected even in
microgels [43, 11, 42], a feature here confirmed for hydrogels at $P=0$, as
also shown for the $c=3\%$ network in Fig. 3(a).
Interestingly, we previously detected the existence of the HAT at $\alpha=0$
in the $c=0.35\%$ network [27] in the absence of attraction. Although in that
case we found comparable phenomenology with the present investigation, here we
notice that the two phenomena appear to be continuous from the behavior of
$\alpha_{HAT}$ in Fig. 3(b). However, the underlying order parameter in that
case was found to be strictly $\rho+e_{nb}$, so that differences between the
attractive and purely repulsive case need to be further investigated.
Having discussed the HAT, we now address the underlying nature of the VPT. We
thus consider simple polymer chains featuring the same interactions as the
hydrogels and evaluate the occurrence of gas (chain poor)-liquid (chain rich)
phase separation at a given critical temperature $\alpha_{c}$, as detailed in
the SI, for systems of different chain lengths $l$ from monomers ($l=1$) to
very long chains ($l\sim 500$). We find that, as expected, the various systems
phase separate at a temperature which depends on $l$. In particular, a larger
temperature is needed for monomers and small chains, while connectivity
enhances the tendency to phase separate, thus lowering the critical
temperature as $l$ increases. Remarkably, for long enough chains, i.e.
$N\gtrsim 100$, we find that phase separation occurs at a constant value
$\alpha_{c}\sim\alpha_{VPT}$. Given the presence of very long chains in any
microgel or hydrogel realization, either in experiments or in simulations,
this finding explains the almost universal value of $\alpha_{VPT}$ and of the
corresponding temperature in experiments, that occurs independently of $c$ for
PNIPAM-based systems in the absence of co-monomers. This is simply due to the
solvophobic interactions between polymer and water, which drives the phase
separation and is mitigated by the presence of crosslinks within the network,
which transforms the critical-like behavior into a smooth and reversible
deswelling process.

In summary, the present study elucidates the multifaceted relationship between
critical swelling, hyper-auxeticity, and Volume Phase Transition in
thermoresponsive polymer networks. Through extensive numerical simulations in
equilibrium and under uniaxial deformation, we demonstrate that the Hyper-
Auxetic Transition represents a clearly distinct phenomenon from the $VPT$,
that is essentially driven by the phase separation of polymer chains, echoing
the coil-to-globule transition of single chains. By analyzing the swelling
curves of the hydrogels, we are indeed able to discern a distinct two-step
process at slightly negative pressures, whereby an initial compression leads
to the onset of a hyper-auxetic point before a further volume decrease takes
place, toward the generic $VPT$ regime. Furthermore, our investigation
spanning a wide range of crosslinker concentrations uncovers hyper-auxetic
behavior up to a finite $c\simeq 3\%$, while higher crosslinker concentrations
exhibit increasing stiffness without reaching the mechanical instability. The
convergence of the Hyper-Auxetic Transition with the $VPT$ as system
connectivity increases underscores the complex nature of polymer network
behavior. Future theoretical work will have to be devoted to the understanding
of the link between the hyper-auxetic transition found in this work in
presence of attractive (solvophobic) interactions and the one found at
$\alpha=0$ [27], where only entropy plays a role, in order to shed light on
the fundamental nature of this transition from the statistical mechanics point
of view. Notwithstanding this open theoretical problem, the present study
provides a useful framework for the understanding of previous and future
experimental observations thanks to the establishment of a clear link between
thermodynamic and mechanical properties of polymer networks. By clarifying the
conditions under which the HAT occurs, the present results call for further
experimental and computational exploration of elastic and thermodynamic
properties of polymer networks, both at negative and at positive pressures,
paving the way for the development of novel materials with controlled auxetic
properties.

We acknowledge funding from ICSC – Centro Nazionale di Ricerca in High
Performance Computing, Big Data and Quantum Computing, funded by European
Union – NextGenerationEU - PNRR, Missione 4 Componente 2 Investimento 1.4 and
CINECA-ISCRA for computational resources.

## References

  * DeRossi  _et al._ [1991] D. DeRossi, K. Kajiwara, Y. Osada, and A. Yamauchi, eds.,  _Polymer gels_ , 1991st ed. (Springer, New York, NY, 1991).
  * Loh and Scherman [2012] X. J. Loh and O. A. Scherman, eds.,  _Polymeric and self assembled hydrogels_ , Monographs in supramolecular chemistry (Royal Society of Chemistry, Cambridge, England, 2012).
  * Koetting  _et al._ [2015] M. C. Koetting, J. T. Peters, S. D. Steichen, and N. A. Peppas, Materials Science and Engineering: R: Reports 93, 1–49 (2015).
  * Onuki [1993] A. Onuki, Theory of phase transition in polymer gels, in  _Advances in Polymer Science_ (Springer Berlin Heidelberg, 1993) p. 63–121.
  * Dimitriyev  _et al._ [2019] M. S. Dimitriyev, Y.-W. Chang, P. M. Goldbart, and A. Fernández-Nieves, Nano Futures 3, 042001 (2019).
  * Tanaka  _et al._ [1987] T. Tanaka, S.-T. Sun, Y. Hirokawa, S. Katayama, J. Kucera, Y. Hirose, and T. Amiya, Nature 325, 796–798 (1987).
  * Matsuo and Tanaka [1992] E. S. Matsuo and T. Tanaka, Nature 358, 482–485 (1992).
  * Chang  _et al._ [2018] Y.-W. Chang, M. S. Dimitriyev, A. Souslov, S. V. Nikolov, S. M. Marquez, A. Alexeev, P. M. Goldbart, and A. Fernández-Nieves, Physical Review E 98, 10.1103/physreve.98.020501 (2018).
  * Hirotsu [1991] S. Hirotsu, The Journal of Chemical Physics 94, 3949–3957 (1991).
  * Li  _et al._ [1993] C. Li, Z. Hu, and Y. Li, Physical Review E 48, 603–606 (1993).
  * Boon and Schurtenberger [2017] N. Boon and P. Schurtenberger, Physical Chemistry Chemical Physics 19, 23740–23746 (2017).
  * Lakes [1987] R. Lakes, Science 235, 1038–1040 (1987).
  * Evans [1991] K. E. Evans, Endeavour 15, 170–174 (1991).
  * Xinchun and Lakes [2007] S. Xinchun and R. S. Lakes, physica status solidi (b) 244, 807–807 (2007).
  * Greaves  _et al._ [2011] G. N. Greaves, A. L. Greer, R. S. Lakes, and T. Rouxel, Nature Materials 10, 823–837 (2011).
  * Nicolaou and Motter [2012] Z. G. Nicolaou and A. E. Motter, Nature Materials 11, 608–613 (2012).
  * Ciarletta  _et al._ [2013] M. Ciarletta, M. Fabrizio, and V. Tibullo, Mathematical Methods in the Applied Sciences 37, 2864–2871 (2013).
  * Pigłowski  _et al._ [2017] P. M. Pigłowski, J. W. Narojczyk, K. W. Wojciechowski, and K. V. Tretiakov, Soft Matter 13, 7916–7921 (2017).
  * Reid  _et al._ [2018] D. R. Reid, N. Pashine, J. M. Wozniak, H. M. Jaeger, A. J. Liu, S. R. Nagel, and J. J. de Pablo, Proceedings of the National Academy of Sciences 115, 10.1073/pnas.1717442115 (2018).
  * Hanifpour  _et al._ [2018] M. Hanifpour, C. F. Petersen, M. J. Alava, and S. Zapperi, The European Physical Journal B 91, 10.1140/epjb/e2018-90073-1 (2018).
  * Tanaka  _et al._ [1977] T. Tanaka, S. Ishiwata, and C. Ishimoto, Physical Review Letters 38, 771–774 (1977).
  * Tanaka [1978] T. Tanaka, Physical Review Letters 40, 820–823 (1978).
  * Escobedo and de Pablo [1996] F. A. Escobedo and J. J. de Pablo, The Journal of Chemical Physics 104, 4788–4801 (1996).
  * Escobedo and de Pablo [1997] F. A. Escobedo and J. J. de Pablo, The Journal of Chemical Physics 106, 793–810 (1997).
  * Escobedo and De Pablo [1997] F. Escobedo and J. De Pablo, Molecular Physics 90, 437–443 (1997).
  * Jha  _et al._ [2011] P. K. Jha, J. W. Zwanikken, F. A. Detcheverry, J. J. de Pablo, and M. Olvera de la Cruz, Soft Matter 7, 5965 (2011).
  * Ninarello  _et al._ [2022] A. Ninarello, J. Ruiz-Franco, and E. Zaccarelli, Nature Communications 13, 10.1038/s41467-022-28026-z (2022).
  * Ninarello  _et al._ [2023] A. Ninarello, J. Ruiz-Franco, and E. Zaccarelli, The Journal of Chemical Physics 159, 10.1063/5.0178409 (2023).
  * Grest and Kremer [1986] G. S. Grest and K. Kremer, Physical Review A 33, 3628–3631 (1986).
  * Grest and Kremer [1990] G. S. Grest and K. Kremer, Macromolecules 23, 4994–5000 (1990).
  * Soddemann  _et al._ [2001] T. Soddemann, B. Dünweg, and K. Kremer, The European Physical Journal E 6, 409–419 (2001).
  * Ninarello  _et al._ [2019] A. Ninarello, J. J. Crassous, D. Paloli, F. Camerin, N. Gnan, L. Rovigatti, P. Schurtenberger, and E. Zaccarelli, Macromolecules 52, 7584–7592 (2019).
  * Gnan  _et al._ [2017] N. Gnan, L. Rovigatti, M. Bergman, and E. Zaccarelli, Macromolecules 50, 8777–8786 (2017).
  * Sorichetti  _et al._ [2021] V. Sorichetti, V. Hugouvieux, and W. Kob, Macromolecules 54, 8575–8589 (2021).
  * Sorichetti  _et al._ [2023] V. Sorichetti, A. Ninarello, J. Ruiz-Franco, V. Hugouvieux, E. Zaccarelli, C. Micheletti, W. Kob, and L. Rovigatti, The Journal of Chemical Physics 158, 10.1063/5.0134271 (2023).
  * Camerin  _et al._ [2018] F. Camerin, N. Gnan, L. Rovigatti, and E. Zaccarelli, Scientific Reports 8, 10.1038/s41598-018-32642-5 (2018).
  * Plimpton [1995] S. Plimpton, Journal of Computational Physics 117, 1–19 (1995).
  * Milton [1992] G. W. Milton, Journal of the Mechanics and Physics of Solids 40, 1105–1137 (1992).
  * Lakes [1993] R. Lakes, Advanced Materials 5, 293–296 (1993).
  * Wilding and Bruce [1992] N. B. Wilding and A. D. Bruce, Journal of Physics: Condensed Matter 4, 3087–3108 (1992).
  * Debenedetti  _et al._ [2020] P. G. Debenedetti, F. Sciortino, and G. H. Zerze, Science 369, 289–292 (2020).
  * Rovigatti  _et al._ [2019] L. Rovigatti, N. Gnan, A. Ninarello, and E. Zaccarelli, Macromolecules 52, 4895–4906 (2019).
  * Voudouris  _et al._ [2013] P. Voudouris, D. Florea, P. van der Schoot, and H. M. Wyss, Soft Matter 9, 7158 (2013).

## Supplementary Information for: Hyper-Auxeticity and the Volume Phase
Transition of Polymer Gels

### Additional Methods

We perform Molecular Dynamics simulations of polymer networks composed of
monomers that interact through the Kremer-Grest potential. Excluded volume for
all particles are given by the Weeks-Chandler-Andersen potential: [1]

$V_{WCA}\left(r\right)=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]+\epsilon\
\ \ if\ \ \ r\leq 2^{1/6}\sigma,\\\ $ (1)

where $\sigma$ is the monomer diameter, which sets the unit of length, and
$\epsilon$ controls the energy scale. Defining $m$ as the mass of the
particles, the unit time of our simulations is defined as
$\tau=\sqrt{m\sigma^{2}/\epsilon}$. Chemical bonds between connected monomers
are modeled by a FENE potential $V_{FENE}\left(r\right)$ [2]:

$V_{FENE}\left(r\right)=-\epsilon
k_{F}R^{2}_{0}\ln\left[1-\left(\frac{r}{R_{0}\sigma}\right)\right]\ \ \ if\ \
\ r<R_{0}\sigma,\\\ $ (2)

where $k_{F}=15$ is the spring constant and $R_{0}=1.5$ is the maximum
extension of the bond. To simulate the presence of the solvent and the polymer
progressively going from a hydrophilic to a hydrophobic condition with
increasing temperature, we employ a so-called solvophobic potential [3] of the
form:

$V_{\alpha}(r)=\begin{cases}-\epsilon\alpha\ \ &\ r\leq
2^{\frac{1}{6}}\sigma,\\\
\frac{1}{2}\alpha\epsilon\big{[}\cos(\gamma(\frac{r}{\sigma})^{2}+\beta)-1\big{]}\
\ &\ 2^{\frac{1}{6}}\sigma<r<R_{0}\sigma,\end{cases}$ (3)

where $\gamma=(\pi(2.25-2^{1/3}))^{-1}$, $\beta=2\pi-2.25\gamma$, $\epsilon$
is the unit of energy and the parameter $\alpha$ controls the strength of the
monomer-monomer attractive interactions.

As anticipated in the main text, we employ networks with two different
topologies. In one case, we build networks with an ordered structure based on
a diamond-like lattice by crosslinking equal length chains where monomers are
placed at the equilibrium distance of the FENE potential. The crosslinkers
concentration is therefore directly determined by the chain length $l$ through
the relationship: $c=1/(2l+1)$. The total number of monomers $N$ of the
networks is $N=\frac{N_{c}}{c}$ where $N_{c}$ is the number of crosslinkers.
In particular, we use $N_{c}=64$ and $l=99,49,25$, respectively for the
$c=0.5\%,1\%,2\%$ systems, while $N_{c}=216$ and $l=16,10$ for the $c=3\%,5\%$
systems, and finally $N_{c}=512$ and $l=6$ for the $c=7.5\%$ system.

In the disordered case, we considered a network with $N=4384$ total number of
monomers, of which $N_{c}=49$ crosslinkers. This is assembled through a
recently developed technique [4], which makes us able to produce low density
computational polymer networks. In this method, two and four valence patchy
particles that mimic respectively monomers and crosslinkers are mixed in
appropriate proportions in a cubic box with periodic boundary conditions to
obtain given connectivity amount. The assembly is performed via NVT molecular
dynamics simulations at low temperature $T=0.3$ and density $\rho=2.85\times
10^{-2}$ with the OxDNA package [5, 6]. We wait until the $99.9\%$ of possible
bonds are satisfied, then remove all particles non belonging to the biggest
cluster and all dangling chains in order to have a fully connected network
without defects. Once the assembly is ready, we then substitute the patchy
reversible bonds with the bead-spring model.

### Bulk and Young moduli

Figure 4: Bulk (a) and Young (b) moduli as a function of $\alpha$ for
$c=1,2,3,5,7.5\%$ and corresponding $P=9.25\times 10^{-4},3.5\times
10^{-3},7.9\times 10^{-3},1.8\times 10^{-2},4.25\times 10^{-2}$.

As discussed in the main text, the behavior of the bulk $K$ modulus for our
hydrogels suggests the presence of an underlying criticality. In particular, a
vanishing bulk modulus is related with the presence of critical-like
fluctuations of the system volume, hinting at the occurrence of a second order
transition as observed for instance in martensitic transformations [7]. In
Fig. 4(a) we report the bulk modulus as a function of $\alpha$ for different
$c$ values, corresponding to state points shown in Fig.3(a) of the main text.
From low $\alpha$ results included in Fig. 4(a) , it can be seen that the
elasticity of the network highly depends on crosslinker concentration, with
values of $K$ spanning two orders of magnitude for increasing $c$. However,
increasing $\alpha$, $K$ shows a dramatic decrease with a minimum reaching
very low values as $c$ decreases, up to about five orders of magnitude with
respect to the high compact state observed at large $\alpha$ when all networks
behave in the same way. In addition, we report in Fig. 4(b) the corresponding
behavior of the Young modulus $Y$. Despite this is not expected to show
critical behavior, we observe the same qualitative trend also present in $K$,
with a clear minimum also developing, becoming more and more pronounced as $c$
decreases. Increasing $\alpha$ to values above the VPT, we find that $Y$
increases and remains distinct for different $c$ values, differently from $K$.

### Phase separation of polymer chains

To allow the comparison with the network results, we also explore chain-based
systems of different lengths $l$. We thus perform MD simulations of chains of
equal length $l=1,2,5,10,49,156,500$ and varying total number of chains
varying in the range between 10000 and 20000.

To allow a quick scan of the phase diagram in order to look for phase
separation, for each chain length we perform several NVT simulations at
different attractive strengths $\alpha$ and system volume. After allowing the
system to equilibrate for $2\times 10^{6}$ time steps, we generate multiple
equilibrium configurations by running simulations for $10^{7}$ time steps.
These configurations are used to calculate the structure factor $S(k)$ of the
network, that is computed removing the (almost negligible) contribution of the
single chain form factor $S_{1}(k)$. Subsequently, we extract the value at
$S(k)$ at the smallest simulated value of $k$, namely
$k_{min}=\frac{2\pi}{L}$, with $L$ representing the box side. In this way, we
have access to a proxy for the isothermal compressibility $\kappa_{T}$ of the
system (that is the inverse of the bulk modulus for network systems), since
$S(k\to 0)$ is directly related to the system through the equation
$\kappa_{T}=\frac{S(0)}{\rho K_{B}T}$ [8]. It is important to note that
$\kappa_{T}$ theoretically approaches infinity at phase separation in
infinite-size systems, but for finite systems, density fluctuations remain
finite, although showing a significant growth.

Figure 5: Value of the structure factor at the minimum accessible wave vector
$S(k_{min})$ as a function of $\alpha$, rescaled by the single chain form
factor at the same point $S_{1}(k_{min})$. (a) Results for the $l=49$ system
for different volumes increasing from bottom to top in the interval $V=1\times
10^{3}-1\times 10^{4}$. (b) Results for systems with various chain lengths as
indicated in the legend at the volume corresponding to the steepest variation.
Vertical lines highlight the inflection points that identify the transition
temperature $\alpha_{c}$ for each $l$.

We report $S(k_{min})/S_{1}(k_{min})$ as a function of $\alpha$ for several
studied isochores for the $l=49$ system in Fig. 5(a) as a representative
example of our analysis. By increasing $\alpha$, a clear increase of
$S(k_{min})$ is observed, which goes through an inflection point at decreasing
value of $\alpha$ with increasing volume. At large volumes, the behavior
saturates, collapsing onto the same curve, with a value of
$S(k_{min})/S_{1}(k_{min})$ jumping by more than two orders of magnitude for a
slight increase of $\alpha$. We identify this inflection point with the
critical value $\alpha_{c}$ that is required to induce phase separation. As we
aim to compare these results with the network system ones, it is particularly
relevant to perform a similar analysis at varying the chain length. This is
reported in Fig. 5(b) where we report $S(k_{min})/S_{1}(k_{min})$ for each
studied $l$ at the volume where the phase-separation temperature $\alpha_{c}$
is smallest. We thus observe a clear decrease of $\alpha_{c}$ with increasing
$l$, approaching the value of the VPT ($\alpha_{VPT}\sim 0.65$) for $l>100$,
as shown in the main text in Fig. 3(b). Since networks generally contains
mostly longer chains, whose effect is dominating on the phase behavior, we
thus conclude that the VPT is an effect purely originating from the
solvophobic attraction between the monomers as $\alpha$ increases and its
universal value $\alpha_{VPT}$ is due to the dominance of long chains in the
systems, both for hydrogels and for microgels.

## References

  * Weeks  _et al._ [1971] J. D. Weeks, D. Chandler, and H. C. Andersen, Role of repulsive forces in determining the equilibrium structure of simple liquids, The Journal of Chemical Physics 54, 5237–5247 (1971).
  * Kremer and Grest [1990] K. Kremer and G. S. Grest, Dynamics of entangled linear polymer melts: a molecular-dynamics simulation, The Journal of Chemical Physics 92, 5057–5086 (1990).
  * Soddemann  _et al._ [2001] T. Soddemann, B. Dünweg, and K. Kremer, The European Physical Journal E 6, 409–419 (2001).
  * Gnan  _et al._ [2017] N. Gnan, L. Rovigatti, M. Bergman, and E. Zaccarelli, In silico synthesis of microgel particles, Macromolecules 50, 8777–8786 (2017).
  * Rovigatti  _et al._ [2014] L. Rovigatti, P. Šulc, I. Z. Reguly, and F. Romano, A comparison between parallelization approaches in molecular dynamics simulations on gpus, Journal of Computational Chemistry 36, 1–8 (2014).
  * Poppleton  _et al._ [2020] E. Poppleton, J. Bohlin, M. Matthies, S. Sharma, F. Zhang, and P. Šulc, Design, optimization and analysis of large dna and rna nanostructures through interactive visualization, editing and molecular simulation, Nucleic Acids Research 48, e72–e72 (2020).
  * Dong  _et al._ [2010] L. Dong, D. S. Stone, and R. S. Lakes, Softening of bulk modulus and negative poisson ratio in barium titanate ceramic near the curie point, Philosophical Magazine Letters 90, 23–33 (2010).
  * Hansen and McDonald [2013] J.-P. Hansen and I. R. McDonald, Introduction, in  _Theory of Simple Liquids_ (Elsevier, 2013).