Doubly robust estimation of causal effects for random object outcomes with continuous treatments

Keywords: Causal inference, Continuous treatment, Fr{\'e}chet regression, Semiparametric efficiency, Air pollution and mortality.

Introduction

Causal inference is pivotal in statistics and scientific research, enabling the identification and estimation of cause-and-effect relationships beyond associations (e.g., rubi:74, rubi:05, holland1986statistics, pearl:16). Going beyond randomized controlled trials coln:24, causal inference provides a powerful toolkit to analyze observational data and infer causal relationships, even in the presence of confounding factors, biases, and noise. By uncovering the underlying mechanisms driving observed phenomena, causal inference approaches facilitate better decision-making and more effective interventions across a broad range of disciplines.

However, most existing work in causal inference has focused on investigating causal effects within linear spaces, particularly within the Euclidean space $\mathbb{R}^p$. In contrast, the rise of complex non-Euclidean data, taking values in general metric spaces that often lack inherent linear structures, has become increasingly prominent in real-world applications. Instances of such data, referred to as “random objects”, include diverse forms such as images, shapes, networks, or life tables marr:14. Other notable examples include symmetric positive definite matrices, networks, spherical surface data, and Riemannian manifolds, among others. Given the metric space nature of the data, conventional statistical concepts such as sample or population means, defined as averages or expected values, do not readily apply and necessitate substitution with notions like barycenters or Fr\'echet means frec:48. In many modern applications, observed data either inherently manifests as or can be abstracted into such complex, non-Euclidean random objects. Often, the primary interest lies in understanding the causal effect on the random objects themselves. Consequently, there is a growing recognition that such applications are suitably characterized using non-Euclidean random objects. Modeling these as metric space-valued stochastic processes preserves their shape and geometry, providing richer information than scalar or vector summaries and necessitating new approaches to causal inference.

A specific application that motivated our work is an environmental study, in which we examine the causal relation between the age-specific mortality distribution and the annual exposure to fine particles (with an aerodynamic diameter of 2.5$\mu$m\ or smaller), denoted as \pmfine\ \ across the U.S.. The National Ambient Air Quality Standards (NAAQS) for \pmfine\ \ set the current primary standard for annual average at $9$ $\mu$g/m$^3$. While a body of literature concluded that exposure to \pmfine\ \ increases the risk of premature death among older adults wu:20, jose:23, the aggregate mortality rate is a scalar random variable that often fails to capture the age-specific mortality of the given region. Our interest is to summarize how the distribution of age-at-death, which is a metric-space-valued random element, can be causally explained by continuously distributed \pmfine\ \ exposure in the presence of confounders while developing theoretical guarantees for the proposed estimator. This presents a significant challenge since the space of distribution-valued random variables lacks an inherent linear structure; as such, basic algebraic operations such as addition or scalar multiplication are not well defined in the space of distributions. However, one can consider the space of distributions, represented as quantile functions, CDF, or density functions, to be a metric space equipped with an appropriate metric such as the Wasserstein or Fisher-Rao metric deli:17, le:17, pana:19.

In this paper, we develop an inverse probability weighted and a double-debiased estimation approaches for causal effects with non-Euclidean outcomes in the presence of moderately high-dimensional vector-valued confounding variables. We refer to this problem as CTROCIN (read as C-trocin), which stands for Continuous Treatment, Random Object Causal Inference. The regime of binary treatment and random object outcomes, referred to as BTROCIN, has been studied very recently by lin:23, kuri:24. In the case of BTROCIN, lin:23 developed a doubly robust estimation approach, which is limited to distributional outcomes, while kuri:24 considered a specific, but often restrictive, modeling framework for outcomes in geodesic spaces. However, the general methodologies and theoretical foundations for CTROCIN remain underdeveloped. Extending from binary to continuous treatments introduces additional complexities, particularly in interpreting treatment effects across varying levels, due to potential selection bias. To address this challenge, we use dose density weights to estimate average causal responses. We propose a nonparametric doubly-debiased inference approach for non-Euclidean outcomes in general metric spaces that allow for embedding into some underlying Hilbert Spaces under the assumption of unconfoundedness given observed covariates.

The main contributions of this paper are as follows:

Section 7 includes a few concluding remarks. The additional numerical results, technical details, and complete proofs are presented in the supplement.

Basic elements of CTROCIN

We consider the continuous evolution of the random object outcome $Y$ in response to a continuously varying treatment $T$, in the presence of possibly high-dimensional confounders $X$, in observational studies. The key ingredients of CTROCIN are summarized as follows: $(\Omega, \ca F, P)$, the probability space; $(\ca T, \ca F \lo T)$, the treatment space, with $T: \Omega \to \ca T\subset \mathbb R$ being the treatment; $(\ca X, \ca F \lo X)$, the covariate space, with $X: \Omega \to \ca X\subset \mathbb R^p$ being the covariate; $(\ca Y, d\lo Y)$, the metric space for the cross-sectional outcome $Y \lo t$ under treatment $t \in \ca T$; $\ca M \lo Y$, the space of functions $\ca T \to \ca Y$, with $Y: \Omega \to \ca M \lo Y$ being the outcome function, denoted by $Y =\{ Y\lo t: t\in \ca T \}$. We use $Y_t$ to represent the random-object-valued response at treatment level $t$. The random elements $X,T,Y$ are measurable with respect to (w.r.t.) $\ca F/ \ca F \lo X$, $\ca F/ \ca F \lo T$, and $\ca F / \ca F \lo Y$, respectively, where $\ca F \lo X, \ca F\lo T,$ and ${\cal{\cal F}}\lo Y$ are the Borel $\sigma$-fields on $\ca X,\ca T,$ and $\ca Y$, induced by their respective metrics. We denote $Z =(Y, T, X) \in \mathcal{Z}:= \mathcal{Y}\times \mathcal{T} \times \mathcal{X}$ from a population $\mathcal{P}$ with the CDF $F_Z$. By construction, the Stable Unit Treatment Values Assumption (SUTVA) holds: each subject’s potential outcome is unaffected by the treatment assignments of other subjects, and each treatment level is well-defined without hidden variations that could lead to different potential outcomes for the same unit.

\def\tfrac#1#2{\textstyle{\frac{#1}{#2}}}

The potential outcome function $Y$ in a causal setting is a function of $T$, $X$, and random unobserved heterogeneity. For Euclidean outcomes, a quantity of interest is often a summary measure of the potential outcome distribution that reflects the change due to the continuous treatment or exposure under the assumption of unconfoundedness (ref). A widely used measure is the average dose response or exposure-response-function (ERF) defined as $E(Y\lo t)$, where $Y\lo t$ denotes the potential outcome under the hypothetical treatment value $T=t$, and the expectation is taken over the distributions of $(X,\epsilon)$, with $\epsilon$ being the unobserved noise. When the outcome is a metric space-valued random object, however, linear functionals such as expectation or the usual additive error structure are not well-defined. Thus, understanding the effect of continuous treatment on random object responses requires leveraging the underlying geometry of the metric space. As such, the definition of an average or expected value is replaced by barycenters or Fr\'echet means frec:48. For any random variable $U$ in a metric space $(\ca M, d)$, its Fr\'echet mean is defined as $E\lo \oplus(U) := \operatorname*{argmin} \{ E[ d\hi 2(U, u)] : u \in (\ca M, d)\}.$ Accordingly, for any potential outcome $Y\lo t$ at treatment level $t\in \ca T$ that takes value in a metric space $(\ca Y,d\lo Y)$, the central tendency, interpreted as the `expected' potential distribution, is defined as the Fr\'echet mean

Henceforward, we will call this the Fr\'echet exposure-dose function (FERF)

A central challenge in causal effect estimation is that we do not observe $(X, T, \{ Y \lo t: t \in \ca T\})$ for each subject. Instead, we observe only $(X, T, Y \lo T)$ for each subject, where $Y \lo T$ is the cross-sectional response at the observed treatment $T$. It is helpful to view $Y\lo T$ as $\int Y\lo t d\delta_T(t)$, where $\delta_a$ is a Dirac measure at $a\in\ca T$, satisfying $\delta_a(B) = 1$ if $a\in B$ and $\delta_a(B) = 0$ otherwise. From the property of the Dirac measure that $\int f(t) d\delta_a(t) = f(a)$, the expression for $Y\lo T = \int Y \delta_T$ follows. This formulation separates the random outcome function $Y$ from the random treatment variable $T$. In particular, this facilitates the derivation of the influence function and its semiparametric efficiency in Section 3 before we describe our proposed doubly robust estimator. We refer to the map $(X, T, Y) \mapsto (X, T, Y\lo T) = (X, T, \int_{\scriptscriptstyle{\ca T}}Y \delta_T )$ as the observation mapping, since only the image of this map is observed.

In general, for any given $t \in \ca T$, $E(Y \lo T|T = t)$ does not give us an unbiased estimate of $E(Y \lo t)$, because $T$ may be affected by other factors, known as confounders, that also affect $T$ and $Y$. To account for confounding, we introduce the covariate $X$, which we assume contains all the confounders: after conditioning on $X$, the response $Y$ no longer depends on $T$. This assumption, known as ignorability, conditional independence, or selection of observables, is standard in the causal inference literature rubi:74, rubi:05:

This assumption asserts that, conditional on observables, the treatment assignment is conditionally exogenous or behaves as if it were randomized. It implies that the observational study is similar to a randomized controlled trial that facilitates valid estimation of causal effects. To achieve valid inference, we adopt a doubly debiased machine learning approach that leverages a doubly robust moment function combined with cross-fitting, as long as the response space admits a suitable embedding into a Hilbert space. Our approach is fully nonparametric, imposing no distributional or functional form assumptions on $T$, $X$, or $\epsilon$.

Embedding of a metric space into a Hilbert space

Embedding metric spaces into simpler and more structured spaces that have low distortion plays an important role in the analysis of random objects, and such embeddings have widespread applications across fields. Proposition 3 of sejd:12 implies that whenever $d\lo Y$ is a semi-metric of negative type, there exists a Hilbert space $\ca H$ and an injective map, say $f : \ca Y \to \ca H$ with $d\lo Y\hi 2(Y_1,Y_2) = \|f(Y_1) - f(Y_2)\|_{\scriptscriptstyle{\ca H}}\hi 2,$ for any $Y_1, Y_2 \in \ca Y$. Thus, if the metric space $(\ca Y,d\lo Y)$ where the outcome function takes values in is of strong negative type, the existence of an isometric continuous embedding from $\ca Y $ to an underlying Hilbert space $\ca H$ is guaranteed. Here, a space $(M,\rho)$ with a semi-metric $\rho$ is of negative type if for all $n \geq 2$, $z_1, z_2, \dots,z_n \in M$ and $\alpha_1,\alpha_2,\dots,\alpha_n \in \mathbb R$, with $\sum_{i=1}^n \alpha\lo i = 0,$ one has $\sum_{i=1}^n\sum_{j=1}^n \alpha\lo i \alpha_j \rho(z\lo i,z_j) \leq 0$. Every separable Hilbert space is of strong negative type. An explicit form of such a continuous, injective, isometric map is given in the following theorem.

Next, we discuss special cases for a metric space to be embeddable in a Hilbert space. The constructions of the Hilbert space embeddings for commonly observed random object data, namely distributional objects, SPD matrices, compositional data, and phylogenetic trees, are presented in the supplement, with more detailed literature review and discussion along with examples.

Embedding for the space of probability distribution: Let $(\ca P,d\lo Y)$ denote the space of probability distributions on a measurable space $(\ca X, \ca B_{\scriptscriptstyle{X}})$, and let $\kappa:\ca X\times \ca X\to \mathbb R$ be a measurable, positive definite kernel with associated RKHS $\ca H_{\scriptscriptstyle{\kappa}}$, such that $\sup_{\scriptscriptstyle{x\in \ca X}} \kappa(x,x)<\infty$. The kernel mean embedding of the probability measure $P$ is defined as the Bochner integral of $\kappa(\cdot,x)$ w.r.t. $P$, that is, $\rho: P\mapsto \int \kappa(\cdot,x) dP(x), \text{ for } P\in \ca P,$ and $\rho$ is a continuous injective map if the kernel is characteristic fuku:04, meaning that $\rho(P)$ uniquely represents $P\in \ca P$, preserving all information. Examples of characteristic kernels on $X = \mathbb R^d$ include Gaussian, Mat\'ern, and Laplace kernels srip:10. Note that the space of the probability distribution is convex; thus, for any $P,Q\in \ca P$ $\lambda P+ (1-\lambda)Q \in \ca P$ for $0<\lambda<1$. Combined with the linearity of the integral operation, this yields that $\rho(\ca P)$ is convex. Furthermore, for continuous and bounded kernel functions, $\rho(\ca P)$ is closed by the portmanteau lemma.

Embedding for the space of SPD matrices and networks: The cone of $K\times K$ symmetric positive semi-definite matrices, $\mathcal{S}_K$, equipped with a suitable choice of metric, such as the Frobenius metric, log-Euclidean metric arsi:07, the power metric family dryd:10,pigo:14, tava:19, the Log-Cholesky metric lin:19, the Bures-Wasserstein metric taka:11, and so on induce a Riemannian manifold structure on $\mathcal{S}_K$ bhat:09.

Embedding for the finite-dimensional Riemannian manifold: In the context of Riemannian manifolds, mapping data into a Hilbert space is well-studied. Embedding into a Reproducing Kernel Hilbert Space (RKHS) can be achieved by using heat kernels bera:94, chu:22 or Gaussian RBF kernels jaya:15, jaya:16, with appropriate adjustments for the curvature of the manifold. Figure (ref) shows an illustration of a Hilbert space embedding for compositional data situated on the surface of a sphere, using a Legendre polynomial embedding map.

Embedding for Phylogenetic trees: Tree space bill:01 is an example that may not admit a Riemannian structure. Phylogenetic trees are widely used in evolutionary biology to represent the ancestral relationships among a set of organisms, and a vector space embedding of tree space is possible using “tropical geometry” (e.g., song:11).

Fr\'echet exposure-dose response function (FERF)

The existence and uniqueness of the Fr\'echet means depend on the nature of the space, as well as the metric considered. For example, in the case of Euclidean responses, the Fr\'echet means coincide with the usual means for random vectors with finite second moments. In the case of Riemannian manifolds, the existence, uniqueness, and convexity of the center of mass are guaranteed afsa:11, penn:18. In a space with a negative or zero curvature, or a Hadamard space, unique Fr\'echet means are also shown to exist (see e.g. stur:03).

An efficient way to compute the Fr\'echet means is to embed the metric space into a Hilbert space, as discussed in the previous subsection, and use the Riesz representation theorem to compute the expectation within that space. This lends more structure to the abstract metric space, enabling the use of Hilbert space geometry, such as the inherent direction and interpretability of linear space, to support nonparametric or functional causal inference.

The following result, which we rely on heavily in the subsequent development, seems not to have been recorded in the literature to the best of our knowledge. So we formally state it here and rigorously prove it in the Supplementary Material. For a random element $V$ taking values in a Hilbert space $\ca H$, we define its expectation as the Riesz representation of the linear functional that maps $f \in {\ca H}$ to $E \langle f, V \rangle _{\scriptscriptstyle{\ca H}} \in \mathbb R$. This linear functional is bounded if $E \| V \| _{\scriptscriptstyle{\ca H}} < \infty$.

Hereafter, we assume the form of the continuous injective map $\rho$ is known. The object responses $Y\lo t$ are thus embedded in the Hilbert space $\ca H$, and the effective outcomes are denoted as $V\lo t=\rho(Y\lo t)$ for each $t\in \ca T$. Let $\ca M \lo V$ be a space of $\ca H$-valued functions defined on $\ca T$. Assume, for each $\omega \in \Omega$, the function $V (\omega)$ defined by $t \mapsto V \lo t(\omega)$ is a member of $\ca M \lo V$. Thus, the mapping $V: \omega \mapsto V(\omega)$ from $\Omega$ to $\ca M \lo V$ defines a random element in $\ca M \lo V$. In this notation system, $\{\ca H, V, \ca M \lo V\}$ are counterparts of $\{\ca Y, Y, \ca M \lo Y\}$ after the Hilbert-space embedding. This embedding $\rho$ allows us to compute $E \lo \oplus(Y \lo t)$ through operations in the Hilbert space. Under Assumptions (ref) and (ref), the FERF, defined as the Fr\'echet mean of the responses $Y\lo t \in (\ca Y,d\lo Y)$ in (ref) can be equivalently written as

The significance of this relation is that $\rho (Y \lo t)$ is a Hilbert-space-valued random element, and its expectation can be computed by Riesz representation and thus can be estimated. In fact, we can rewrite the FERF estimand in a computationally cheap and interpretable way under the Hilbert space embedding as follows. Our strategy is to estimate $E(V \lo t)$ in a Hilbert space and then transform the result by $\rho ^{\scriptscriptstyle{\scriptscriptstyle -1}}$ to estimate $\beta\lo t= E \lo \oplus (Y \lo t)$.

In the following, we use $m(t, X)$ to denote $E \lo \oplus ( Y \lo t | X = x)$, and call it the Fr\'echet conditional ERF. Let us denote $\gamma(t,x):= \rho(m(t,x))$. Further, the causal effects map between two different levels of treatments $t$ and $t'$ can be quantified as $\Delta_{\scriptscriptstyle{t,t'}}:= \mathbb{E}(\|V\lo t - V_{\scriptscriptstyle{t'}}\|_{\scriptscriptstyle{\ca H}})$, again via the embedding map $\rho$.

Estimation

In this section ,we introduce two estimators of the CITROCIN: one based on inverse probability weighting (IPW), and the other on doubly robust estimation (DR). We assume that $\rho$ is a known bijective transformation satisfying Assumptions (ref)--(ref). After the transformation, we have an i.i.d. sample $\{X_i, T_i, V_i\}_{i=1}^n$. As mentioned earlier, our strategy is to estimate $E(V \lo t)$ and then employ the relation $E \lo \oplus (Y \lo t) = \rho ^{\scriptscriptstyle{\scriptscriptstyle -1}} ( E (V \lo t))$ to estimate $ E \lo \oplus (Y \lo t)$. Since we observe $V \lo T$ but not $V \lo t$, our estimator must be based on $V \lo T$ instead of $V \lo t$, and this is the main challenge in causal estimation.

Inverse probability weighting approach for CTROCIN

This subsection constructs a causally unbiased estimate of $\rho(\beta\lo t) = E(V \lo t)$ by inverse probability weighting (IPW); see, for example, imbe:15 and pearl:16. Since $T$ is continuous, the probability of observing $T = t$ is zero. Thus, combining information from nearby points is necessary. We approach this by first approximating $\beta \lo t$ as follows. Let

where the integral is taken as Brochner integral hsing2015theoretical, $K \lo h (u) = k (u / h) / h$, $h$ is a positive constant, and $k$ is a probability density function defined on $\ca T$. By construction, $\lim _{\scriptscriptstyle{h \to 0}} \vartheta \lo t (h) = \beta \lo t. $ So, to construct a consistent estimate of $\beta \lo t$, we let $h\to 0$ as $n \to \infty$.

Subsequently, we consider a kernel function $K_h$ satisfying the following integrability assumptions

For an $a \in \ca T$, let $\delta \lo a$ be the Dirac measure at $a$. Let $w: \ca T \to \mathbb R$ be an arbitrary nonnegative measurable function on $\ca T$. In particular, it could be the function $s \mapsto K \lo h ( s - t)$ for the fixed $t$ in $V \lo t$. Our IPW estimator is based on the following population-level result. The key rule we have to obey when constructing a causal estimator of $E(V \lo t)$ is that we are only allowed to use $V \lo T$, as $V \lo t$ is not observed for any fixed $t \in \ca T$. Note that, in terms of the Dirac measure $\delta \lo a$, $V \lo T$ can be rewritten as the integral form $\int _{\scriptscriptstyle{\ca T}} V \lo s d \delta \lo T (s)$. The benefit of using this alternative expression is that it separates $V$ from $T$, allowing us to apply the conditional independence $Y \;\, \rule[0em]{.03em}{.65em} \hspace{-.41em} \rule[-.02em]{.65em}{.03em} \hspace{-.41em} \rule[0em]{.03em}{.65em}\;\, T |X$ in a more intuitive fashion.

The point of this equality is that the left-hand side depends on $V \lo T$ but not $V \lo t$; whereas the right-hand side does depend on $V \lo t$. The IPW estimator is based on the mimicry of the left-hand side of the above equality at the sample level. In fact, if we choose $w(t)$ to be the kernel function $K \lo h ( s- t)$, then, by letting $h \to 0$, we can prove the following limit form of the above equality in Proposition (ref). The conditional density $f_{T|X}$, a.k.a. the General Propensity Score (GPS), plays a central role. The following assumption on the smoothness of the GPS is used in the subsequent sections whenever needed.

As mentioned earlier, we choose the weight function $w(t)$ in Proposition (ref) to be a kernel function $k ((s-t)/h)/h$ for some probability density function on $\ca T$.

Doubly robust approach for CTROCIN

Inverse probability weighting described in the last subsection is essentially a moment estimator that allows re-expressing the moment $E(\int _{\scriptscriptstyle{\ca T}} w(t) V \lo t d t )$ in terms of the observable variable $V \lo T$. An alternative approach is the doubly robust estimator, which is semiparametrically efficient under regularity conditions. In this section, we develop such an estimate. Similar to the IPW case, we first target $\vartheta \lo t ( h) = E(\int _{\scriptscriptstyle{\ca T}} K \lo h (s - t) V \lo s d s )$, and then let $h$ go to 0 to estimate $E (V \lo t)$. The next theorem gives the semiparametrically efficient influence function for estimating the real-valued parameter $E(\int _{\scriptscriptstyle{\ca T}} w(s) V \lo s d s )$, where $w(s)$ is a general weighting function, and when the infinite-dimensional nuisance parameters -- $f \lo X$, $f _{\scriptscriptstyle{T|X}}$ and $\{f _{\scriptscriptstyle{V \lo t |X}}: t \in \ca T \}$ --- are completely unknown.

In the binary treatment setting where $\ca T = \{0,1\}$, the potential outcome framework closely aligns with the missing at random (MAR) problem, allowing the average treatment effect (ATE) at two levels to be formulated as a semiparametric estimand. Here, ATE is the parametric component, while functions $f\lo X, f_{\scriptscriptstyle{T|X}}, f _{\scriptscriptstyle{V \lo 0|X}}, f _{\scriptscriptstyle{V \lo 1 | X}}$ are the nonparametric components or the infinite-dimensional nuisance parameters. This setting admits only one influence function in this scenario, so it is semiparametrically efficient. This theory easily extends to the case of finite treatment levels $k$, where standard semiparametric theory can be used to derive the efficient influence function, efficiency bound, and corresponding estimator. However, extending this framework to the continuous treatment setting is nontrivial, as standard semiparametric theory does not apply directly. To address this gap, the key insight of our proof of Theorem (ref) is to introduce the random Dirac measure $\delta \lo T$ to separate the random function $V$ and the random vector $(X, T)$.

Since our goal is to estimate $\rho ( \beta \lo t) = \vartheta \lo t = E (V \lo t)$ instead of $\vartheta \lo t (h)$, we further derive below the limiting (as $h \to 0$) moment condition derived from the efficient score ((ref)).

Mimicking the limit form of the efficient influence function given in Corollary (ref), we now propose the following kernel-based doubly-debiased machine learning estimator as

Again, note that the observable $V \lo T$, rather than the unobservable $V \lo t$, appears in the above equation. The two sample-level conditional expectation $\widehat E (V \lo T |X, T )$ are computed by performing nonparametric regression: ${ \operatorname*{argmin}_{v \in \rho ( \ca Y) } \ \widehat E (\| V \lo i - v \| _{\scriptscriptstyle{\ca H}} \hi 2 | X \lo i , T \lo i )} $ by smoothing spline, kernel regression, or RKHS. The estimation strategies for these infinite-dimensional nuisance parameters are discussed in more detail in Section 4.

The above proposition shows that

from which it is easily deduced that our proposed estimator $\hat \vartheta _{\scriptscriptstyle{t, DR}}$ is consistent even if $\gamma (t,x)$ is misspecified as $\tilde \gamma (x,t)$. The proof requires the structures of the latent Hilbert space, the law of iterated expectation, and the conditional independence assumption, and can be found in Supplementary Material.

Cross-fitting estimator for CTROCIN

While the doubly robust estimator is useful, it requires strong and unverifiable assumptions on the metric space, such as the Donsker property for the function class of the outcome regression. To circumvent these constraints while preserving the desired asymptotic properties, we adopt a sample splitting strategy proposed in cher:18 and cola:22. Building on this idea, let $\hat{\gamma}_{l}(t,x)$ and $\hat{f}_{l}(t,x)$ be suitable estimators based on $(\hat{V}_i, T_i,X_i)_{i=1}^n$ for $\gamma(t,x)$ and $f_{T|X}(t|x)$, respectively, and we propose a kernel-based estimator that utilizes the following double-debiased moment function and a cross-fitting strategy.

Theory

In this section, we first describe the estimation of the auxiliary quantities involved in (ref) and then derive results regarding the asymptotic distribution for the proposed estimator in (ref). Before proceeding, we define the following norms:

First, the outcome curves $\hat{V}_i$ (or $\hat Y \lo i$) need to be constructed from the discrete observations on $V_i$ (or $Y \lo i$), $i=1,\dots, n.$ About these estimates, we make the following convergence assumptions.

Proposition (ref) ensures consistent estimation of the outcome trajectories from the data that are not fully observed. Next, we establish the asymptotic distribution of the proposed estimate $\hat \vartheta \lo t$. To do so, we make some key assumptions about the convergence rates of the estimated GPS and conditional expectation, and so on, and describe some existing estimation procedures of $f _{\scriptscriptstyle{T|X}}$ and $\gamma (t,x)$.

We further require the following assumption:

We discuss the estimation of two nuisance parameters. For the consistent estimation of the quantities involved in (ref), we can use any suitable estimates of the conditional Fr\'echet mean dose-response $m(t,X)$ and the GPS $f_{T|X}(t|x)$. For example, one can employ any local or global object regression method to estimate the above using data in the observations not present in the $l$-th bin (as described in Step 2 of the algorithm), $l\in\{1,\dots, L\}$. We discuss some available options in the literature as follows:

To simplify our understanding in the transformed random variables now taking values in a Hilbert space, we can also view the transformed conditional mean $\gamma(t,x)= \rho(m(t,X))$ as a unique Riesz representation. As the least squares regression fits nicely into the Hilbert space setting, the conditional mean can be perceived as the best linear predictor using the orthonormal basis of the space once the transformation is done rosi:01.

The estimation of the GPS $f_{T|X}(t|x)$, on the other hand, is a well-studied problem: available estimating procedures include a re-weighted Nadaraya-Watson or locally linear estimator fan:96, orthogonal series estimators whit:58,wats:69, penalized quantile regression methods bell:19, cola:22, neural networks mcca:13,roth:19, and so on. Essentially, any estimators that satisfy Assumption (ref) would be a candidate for the nuisance parameters estimation from the underlying outcome curves $V_i$, $i=1,\dots,n$. To tie everything together, the conditional Fr\'echet mean estimator of $\gamma(t,x)$ based on the unobserved underlying quantities $V_i$'s and that based on the estimated trajectories $\hat{V}_i$'s are required to be asymptotically close:

Furthermore, the following technical condition is required for deriving the asymptotically linear representation and asymptotic normality for the proposed DML cross-fitting estimator.

The above theorem yields an asymptotic inferential framework for the FERF and can also be used to quantify the average causal effect between two treatment levels. We define the causal effect map between treatment levels $t$ and $t'$ as $\Delta_{\scriptscriptstyle{tt'}} = \rho(\beta\lo t) - \rho(\beta_{\scriptscriptstyle{t'}})$. An exact finite sample guarantee for the uncertainty quantification of the estimators $\hat{\Delta}_{\scriptscriptstyle{tt'}}$ can be obtained using the adaptive HulC method by kuch:23 to construct confidence regions for the contrast $\Delta_{\scriptscriptstyle{tt'}}$, with the following implementation.

Let $\{S_b\}_{b=1}^{B}$ be a (random) partition of $\{1, \dots, n\}$ into $B$ subsets, and the estimators $\hat{\Delta}_{\scriptscriptstyle{tt'}}^{\scriptscriptstyle{b}} := \{\hat \vartheta_{\scriptscriptstyle{b;t}}) - \hat \vartheta_{\scriptscriptstyle{b;t'}})\}_{b=1}^{B}$ for $ \Delta_{\scriptscriptstyle{tt'}} =\rho(\beta\lo t) - \rho(\beta_{\scriptscriptstyle{t'}})$ computed for each subsample $\{(V_i, T_i, X_i) : i \in S_b\}_{b=1}^{B}$. Define the maximum median bias of the estimators $\{\hat{\Delta}_{\scriptscriptstyle{tt'}}^{\scriptscriptstyle{b}}\}_{b=1}^{B}$ for $\Delta_{\scriptscriptstyle{tt'}}$ as

\[ \Delta := \max_{1 \leq b \leq B} \left\{ 0, \frac{1}{2} - \min \left\{ P(U_b \geq 0), P(U_b \leq 0) \right\} \right\}, \] where $U_b = \hat{\Delta}_{\scriptscriptstyle{tt'}}^{\scriptscriptstyle{b}} - \Delta_{\scriptscriptstyle{tt'}}$. Adopting the HulC algorithm, we construct a confidence interval with coverage probability $1 - \alpha$ for $\Delta_{\scriptscriptstyle{tt'}}$ as follows:

From Theorem 1 in kuch:23, the coverage probability of this confidence interval is guaranteed for finite samples, i.e., $ P\left( \Delta_{\scriptscriptstyle{tt'}} \in \hat{C}_{\alpha, \Delta} \right) \geq 1 - \alpha.$

Simulation studies

The four estimators -- Outcome Regression (OR), Inverse Probability Weighting (IPW), Doubly-Robust (DR), and Doubly-Debiased Cross Fitting (CF) -- are introduced in logical progression in Section 3. To assess their performance, we conduct simulations across various settings involving different types of metric space-valued responses: univariate distributions equipped with the Wasserstein metric, covariance matrices equipped with the Frobenius metric, and six-dimensional compositional data on the positive segment of the sphere $\ca S\hi 6 \lo +$ equipped with the geodesic metric on $\ca S\hi 6$. It is important to note that we do not have any existing methods to compare with since the CTROCIN framework is mostly unexplored as of now. We consider sample sizes $n= 50,200,$ and $1000$ and report the average and standard deviation of leave-one-out mean squared error (MSE) over $B= 100$ Monte Carlo replications. For the $b^{\text{th}}$ simulation, the MSE is given by \[ \text{MSE}^{\scriptscriptstyle{(b)}} = \frac{1}{n} \sum_{i=1}^n d^2(Y_i, \hat{Y}_{(i)}), \] where $Y_i$ is the observed response for the $i^{\text{th}}$ sample, $\hat{Y}_{(i)}$ is the predicted value for $Y_i$ obtained by fitting the model on the remaining $n-1$ observations and predicting $Y_i$. To compute the MSE efficiently, each response $Y_i$ is first embedded into a Hilbert space via a known isometric map $\rho$, resulting in $V_i = \rho(Y_i)$. This transformation allows MSE to be computed as squared distances using the inner product in the embedding space.

In all simulations, the confounder or pre-exposure covariates $X$ and the treatment variable conditioned on $X$ are generated as follows: We generate six pre-exposure covariates $(X_1, X_2,\dots, X_6)$ as a combination of continuous and categorical variables: $$X_1,\dots,X_4 \sim N(0,I_4),\ X_5 \sim V\{-2,2\},\ X_6 \sim U(-3,3),$$ where $N(0,I_4)$ denotes a multivariate normal distribution, $V\{-2,2\}$ denotes a discrete uniform distribution, and $U(-3,3)$ denotes a uniform distribution. We generate $T$ using three different specifications of the GPS model, all relying on the cardinal function $r(X) = - 0.8 + (0.1,0.1,-0.1,0.2,0.1,0.1)^\top X$. The coefficients of the cardinal function $r(X)$ are modified from wu2024matching. Specifically, we consider the following scenarios:

Scenario 1 serves as the baseline, where the exposure $T$ is generated as a linear function of the confounders, and the residuals are normally distributed without extreme values. Scenario 2 introduces heavy-tailed behavior by generating $T$ from a t-distribution, leading to extreme values and, consequently, extreme GPS values. Scenario 3 can be seen as a variant of scenario 1 by incorporating a more complex data-generating process and deliberately misspecifying the GPS model, thereby testing robustness under model misspecification.

For space considerations, we only present simulation results for distributional responses. The additional simulation results for SPD matrix objects and compositional data taking values on the surface of the sphere $S^2 \subset \mathbb R^3$ can be found in the supplement.

We generate $Y$ from an outcome model that is assumed to depend on the treatment and confounders. To this end, we consider two settings as follows.

Setting (A) assumes Gaussian outcomes, while Setting (B) introduces greater complexity by generating non-Gaussian distributions. In Setting (B), distribution parameters are first sampled as in Setting (A), but the resulting distributions are then “transported” in Wasserstein space following pete:19 and chen:22.

The simulated outcomes, represented as quantile functions, are embedded in the Hilbert space $L\hi2 [0,1]$ and taken as the effective outcomes. The estimated GPS $\hat{f}_{t|X}$ is computed by a cross-validated Super Learner ensemble algorithm, implemented by the R package SuperLearner with the extreme gradient boosting machines algorithm SL.xgboost. The other outcome regression function $\gamma(t,x)$ is estimated using smoothing splines via the R function smooth.spline, employing its default generalized cross-validation for tuning.

Table (ref) compares the performance of the OR, IPW, DR, and CF estimators. Overall, all estimators successfully recover the true outcome distribution. The MSE of the doubly robust estimator decreases as the sample size increases. For more complicated data-generating mechanisms (e.g., combinations of scenario 2 for the propensity score model and outcome generation model B), the MSEs are generally higher. Furthermore, neither the OR nor IPW estimator demonstrates double robustness under model misspecification, displaying large MSEs with high variances even with a sample size of $n =1000$. The DR and CF estimators offer improved estimations. These results are in line with our theoretical analysis.

Table (ref) shows the average width of the $95\%$ confidence band obtained from the asymptotic properties of the CF estimate under various simulation scenarios and varying sample sizes, again over $100$ Monte Carlo simulation runs. The width of the confidence band generally narrows with increasing sample size.

Air pollution and health application

The key scientific question in air pollution epidemiology studies is to assess whether and to what extent exposure to air pollution is causally linked to adverse health outcomes. Specifically, we aim to apply our proposed doubly debiased approach to estimate the causal FERF of long-term \pmfine\ \ exposure on all-cause mortality. The response is the age-at-death densities, which reside in the metric space of univariate distributions. Equipped with the Wasserstein-2 metric, this space of distributions can be isometrically embedded into the Hilbert space $L^2[0,1]$. The mortality data for various counties in the United States are obtained from the Centers for Disease Control and Prevention (CDC) website. We use the data for the year $2012$, which is available in the form of cohort life tables over the age interval $[0,110]$ for each of $n = 2392$ counties in the U.S.. For each county, the life table data that correspond to histograms are smoothed with a bin width of five years by adding a smoothing step using the R package fr\'echet fr:package with a bandwidth of 2 years. Thus, a sample of $n =2392$ age-at-death densities is obtained, which are then embedded in the underlying Hilbert space $L^2[0,1]$ to produce the outcomes of interest.

The treatment data is collected as the county-level long-term exposure to \pmfine\ \ (averaged from $2000$ to $2016$) from an established exposure prediction model. The resources used are discussed below. The estimated daily concentrations of \pmfine\ \ are available for 1-km square grids for the contiguous USA between $2000$ and $2016$ di2019ensemble. These measurements are generated using an ensemble-based model that fuses predictions from three machine learning methods: a random forest regression, a gradient boosting machine, and a neural network. Each model uses more than 100 predictor variables derived from satellite data, land-use data, weather measurements, and output from chemical transport model (CTM) simulations. The ensemble was trained on daily \pmfine\ \ concentrations measured at $2,156$ U.S. EPA monitoring sites, with an average cross-validated $R^2$ of $0.86$ for daily \pmfine\ \ predictions and $0.89$ for annual \pmfine\ \ predictions, indicating excellent performance. These predictions have been used in previous high-impact studies evaluating \pmfine\ \ and mortality jose:23, wei2021emulating. In short, CTM and satellite data are combined to estimate a high-resolution \pmfine\ \ surface across the whole United States. This surface is bias-corrected for ground-monitor \pmfine\ \ observations using a geographically weighted regression. We aggregated these levels spatially by averaging the values for all grid points within a county to obtain the temporally averaged \pmfine\ \ values (2000–2016) at the county level by averaging estimated \pmfine\ \ values within a given county.

We first estimate the GPS by using an Extreme Gradient Boosting Machine (GBM) (i.e., a single learner in the Super Learner algorithm) chen2016xgboost, zhu2015boosting, with the county-level \pmfine\ \ exposure as the dependent variable and 19 zip code-level potential confounders as independent variables. The latter include population demographics (%Female, %Black, %Hispanic populations, Population density), health-related (Mean BMI, % Ever Smoked), and educational (% Below High School Education) information regarding socioeconomic status (Median Home Value, Median Household Income, % Owner-occupied Housing), and meteorological information (Summer and Winter minimum and maximum temperature, humidity). The extreme GBM learner is desirable for the estimation of the infinite-dimensional nuisance parameter $f_{T|X}$ for 19-dimensional confounder $X$, as described, since the method is flexible, computationally feasible on a large dataset, and achieves better covariate balance compared to a linear regression model on the complex application data.

Next, for estimating the outcome regression model $E(\rho(Y_t)|T=t, X)$, we implemented a global linear Fr\'echet regression pete:19, bhat:23 model using the Hilbertian quantile functions $\rho(Y_t)$ across $n= 2392$ counties as responses and the confounders $X$ as predictors, at every given level $t$ of the treatment, \pmfine\ \ exposure. Finally, we implemented the cross-fitting estimator in (ref) using $L = 100$. The left panel of Figure (ref) shows the changes in the potential age-at-densities for varying levels of the treatment. The counterfactual densities are color-coded such that blue to red indicates smaller to larger values of the \pmfine\ \ exposures. We find that a lower treatment level can be causally linked with left-shifted age-at-death distributions in general, while higher treatment levels result in a shift of the mode of the age-at-death toward the right. We also plot the $95\%$ pointwise confidence band at three different levels of the \pmfine\ \ exposure, namely, at the $10^{\text{th}}, 50^{\text{th}}$, and $90^{\text{th}}$ percentiles of the treatment values (see the right panel of Figure (ref)). However, the bands are very narrow, which makes the interpretation of the potential outcome distribution rather obscure, possibly indicating a heterogeneous group variation in the treatment effects.

To examine heterogeneity in the causal effect of \pmfine\ \ exposure, we divide the sample into four disjoint groups based on socio-economic characteristics, following jose:23. Specifically, we stratify counties by combinations of high or low income and high or low percentage of Black residents. The CF estimator is implemented for four groups separately. The resulting counterfactual age-at-death distributions exhibit different shapes and structures across these four groups and resonate with the finding of jose:23. As shown in Figure (ref), lower \pmfine\ \ exposure is causally associated with lower mortality in the full population, but marginalized subpopulations appear to benefit more as the \pmfine\ \ levels decrease (see Figure (ref)). For example, the mode of the distribution shifts toward the top-right when the \pmfine\ \ levels are lower for the low-income, high-black% group, indicating a higher longevity corresponding to lower \pmfine\ \ levels. The child mortality is also lower with a lower \pmfine\ \ level. The effect is not so evident for the high-income-low black group.

We also compute the causal effect map for pairwise comparisons of three treatment levels (at 5%, 50%, and 95% values of the treatment levels, respectively). We evaluate a $95\%$ pointwise confidence band obtained from the asymptotic distribution of our proposed estimate using Theorem 2. We also implement the $95\%$ pointwise confidence band using the HulC method kuch:23. Figure (ref) shows the causal effect maps for the contrast of $t= 90\%- t' = 50\%$, $t= 50\%- t' = 10\%$, and $t= 90\%- t' = 10\%$, where for each contrast level, the confidence band obtained from the two methods are overlaid. The confidence band derived from the asymptotic distribution (in orange) is tighter than the general method of minimizing maximum bias, as designed in kuch:23 (in blue).

Figure (ref) shows the $95\%$ confidence bands for each group. As a conclusion, we might infer that the groups with higher income-higher Black%, low income-low Black%, and low income-higher Black% populations may benefit more from lower \pmfine\ \ levels than higher income low Black% groups. These findings underscore the importance of considering racial identity and income when assessing health inequities.

Conclusion

We proposed a method for estimating the causal effect of a continuous treatment on random object response, assuming that the metric space for the response can be embedded in a Hilbert space. This covers many commonly observed random objects such as distributional data, SPD matrices, data on a Riemannian manifold, etc. However, this embedding assumption imposes certain limitations on the method. For example, in certain cases, it may be impossible to embed the metric space or the form of the embedding map may be unknown. This occurs, for example, in spaces such as phylogenetic trees with the BHV or hyperbolic metric bill:01, mata:24. Moreover, the embedding map is not necessarily bijective, making it non-trivial to project back to the original metric space. This highlights the need for intrinsic methods, in which all model components and fits are defined directly within the metric space schotz2021frechet, bhat:23. Alternatively, we can consider a general method for CTROCIN without Hilbert space embedding by focusing on the metric $d (Y,y)$ rather than the random object $Y$ itself. This insight provides a useful paradigm for performing Fr\'echet regression in settings where embedding into a Hilbert space is not feasible. In particular, we can define a causally unbiased estimate of the metric $d \hi 2 (Y \lo t, y)$ for each $y \in \ca Y$ using the IPW or doubly robust estimate and allows for a semiparametric efficient (a.k.a. doubly robust) estimate. This calls for a future research agenda in the intersection of random object data analysis and causal inference.

{

}