Multivariate incremental effects for continuous treatments: Studying the health effects of environmental mixtures

Let our observed data consist of $n$ independent and identically distributed samples $\{\boldsymbol{Z}_{i}\}_{i=1}^{n}$, drawn from some underlying distribution $P_{0}$. Each observation $\boldsymbol{Z}_{i}=(\boldsymbol{X}_{i},\boldsymbol{W}_{i},Y_{i})$ is composed of a $p$-dimensional vector of covariates $\boldsymbol{X}_{i}\in\mathcal{X}\subseteq\mathbb{R}^{p}$, a $q$-dimensional vector representing the environmental exposures or treatments¹¹1Note that we use the terms treatment and exposure interchangeably throughout the manuscript. Also note that environmental mixture simply refers to a vector of environmental exposures. $\boldsymbol{W}_{i}=(W_{i1},\dots,W_{iq})\in\mathcal{W}\subseteq\mathbb{R}^{q}$, and a scalar outcome of interest $Y_{i}\in\mathbb{R}$. We denote the conditional density of the exposure mixture given the covariates as $f(\boldsymbol{w}\mid\boldsymbol{x})$.

To formally define causal effects, we operate within the potential outcomes framework. For each individual $i$ and any potential exposure vector $\boldsymbol{w}\in\mathcal{W}$, we let $Y_{i}(\boldsymbol{w})$ denote the potential outcome that would have been observed had individual $i$ received exposure level $\boldsymbol{w}$. This notation relies on the Stable Unit Treatment Value Assumption (SUTVA, (Rubin, 1980)), which comprises two key principles:

1. 1.

   No Interference: The potential outcome for one individual is unaffected by the exposure assignments of other individuals. That is, $Y_{i}(\boldsymbol{w})$ depends only on the exposure $\boldsymbol{w}$ assigned to individual $i$.

2. 2.

   No multiple versions of treatment: Treatment is well-defined in the sense that there are not two distinct treatments that lead to the same value of $\boldsymbol{W}_{i}$.

Importantly, this assumption ensures that an individual’s observed outcome corresponds to their potential outcome under their observed exposure. Formally, if an individual $i$ is observed to have exposure $\boldsymbol{W}_{i}$, their observed outcome is $Y_{i}=Y_{i}(\boldsymbol{W}_{i})$.

In this work, we extend the framework of exponential tilting incremental causal effects to the multivariate treatment setting. This type of stochastic treatment was first proposed for binary treatments (Kennedy, 2019) and later generalized to single continuous treatments (Díaz and Hejazi, 2020; Schindl et al., 2024). We adapt this formulation to define interventions on the entire $q$-dimensional exposure vector, $\boldsymbol{W}$. Given the conditional density of the exposure mixture, $f(\boldsymbol{w}\mid\boldsymbol{x})$, we define an exponentially tilted interventional density, $g_{\boldsymbol{\delta}}(\boldsymbol{w}\mid\boldsymbol{x})$, indexed by a user-specified vector $\boldsymbol{\delta}\in\mathbb{R}^{q}$:

Here, the denominator is a normalizing constant that ensures $g_{\boldsymbol{\delta}}$ integrates to one, and $\boldsymbol{\delta}$ is a vector determining both the direction and magnitude of how the natural density is shifted.

Our causal estimand of interest, the incremental effect, is the expected potential outcome under the stochastic intervention defined by the tilted density $g_{\boldsymbol{\delta}}$. We denote this estimand as $\psi(\boldsymbol{\delta})$:

This represents the population average outcome if, for covariates $\boldsymbol{X}=\boldsymbol{x}$, each individual’s exposure were a random draw from the shifted distribution $g_{\boldsymbol{\delta}}(\boldsymbol{w}\mid\boldsymbol{x})$. In order to identify this quantity from the observed data, we must make a standard no unmeasured confounding assumption: for all $\boldsymbol{w}\in\mathcal{W}$,

We develop sensitivity analysis approaches to assess violations of this assumption in Section 5. Under SUTVA and the no unmeasured confounding assumption, this causal quantity is identified from the observed data distribution via:

where $P(\boldsymbol{x})$ denotes the marginal probability measure of the covariates $\boldsymbol{X}$. Note that because the support of $g_{\boldsymbol{\delta}}$ is identical to the support of $f$, we do not have to invoke any positivity assumptions and the intervention does not require us to estimate outcomes for exposure combinations that are never observed in the data.

The parameter $\boldsymbol{\delta}$ has a similar interpretation as the constant gradient of the log-likelihood ratio between the interventional and observational densities as in (Schindl et al., 2024):

This means that each component, $\delta_{j}$, quantifies the change in this log density ratio for an infinitesimal increase in the $j$-th exposure component, $w_{j}$. Intuitively, setting $\boldsymbol{\delta}=(1,0,\dots,0)^{\top}$ defines an intervention that tilts the distribution to favor higher values of the first exposure, $W_{1}$.

Let $\mu(\boldsymbol{x},\boldsymbol{w})=\mathbb{E}[Y\mid\boldsymbol{X}=\boldsymbol{x},\boldsymbol{W}=\boldsymbol{w}]$. Following the same derivation as (Schindl et al., 2024) but generalized to the multivariate setting, we get the following formula for the efficient influence function.

Note the subscript $g_{\boldsymbol{\delta}}$ denotes expectations with respect to the tilted exposure density. The asymptotic variance of any regular and asymptotically linear (RAL) estimator for $\psi(\boldsymbol{\delta})$ is equal to the variance of its EIF, $\mathrm{Var}(\varphi(\boldsymbol{Z};\boldsymbol{\delta}))$. A critical observation from the structure of $\varphi$ is the repeated appearance of the density ratio, $g_{\boldsymbol{\delta}}(\boldsymbol{W}\mid\boldsymbol{X})/f(\boldsymbol{W}\mid\boldsymbol{X})$. When this ratio exhibits high variability, it inflates the variance of the first two components of the EIF, leading to less precise estimates of the causal effect. Therefore, to improve statistical efficiency for a given intervention strength, we could select an intervention direction $\boldsymbol{\delta}$ that minimizes the variance of this density ratio, $\mathrm{Var}(g_{\boldsymbol{\delta}}/f)$. The third term, on the other hand, represents the deviation between the average post-intervention effect for a specific covariate group and the overall average post-intervention effect. It quantifies how much higher or lower the expected intervention effect for an individual with covariates $\boldsymbol{X}$ is compared to the overall population average.

To gain analytical insight into this variance, we can let the exposures follow a multivariate normal distribution, $f(\boldsymbol{w}\mid\boldsymbol{x})\sim\mathcal{N}(\boldsymbol{\mu_{x}},\boldsymbol{\Sigma})$. Under this assumption, we have two key results. First, the tilted distribution $g_{\boldsymbol{\delta}}(\boldsymbol{w}\mid\boldsymbol{x})$ is also normal, with a shifted mean:

This provides a clear interpretation of the intervention: it is a shift of the mean of the exposure distribution in the direction $\boldsymbol{\Sigma}\boldsymbol{\delta}$. Second, the variance of the density ratio has a simple, closed-form expression:

Minimizing this variance is equivalent to minimizing the quadratic form $\boldsymbol{\delta}^{\top}\boldsymbol{\Sigma}\boldsymbol{\delta}$.

The question of interest thus becomes an optimization problem: which direction $\boldsymbol{\delta}$ minimizes $\boldsymbol{\delta}^{\top}\boldsymbol{\Sigma}\boldsymbol{\delta}$ subject to a constraint on the “strength” of the intervention? While various constraints are possible (e.g., fixing $\boldsymbol{\delta}^{\top}\boldsymbol{\delta}$), a particularly meaningful constraint is to fix the distance between the original and tilted distributions, for which the 2-Wasserstein distance provides a natural metric. We discuss the choice of this metric in the following section. For the multivariate normal case with the same covariance matrix, the squared Wasserstein distance has a simple closed-form expression as $d_{W}^{2}(f,g_{\boldsymbol{\delta}})=||(\boldsymbol{\mu_{x}}+\boldsymbol{\Sigma}\boldsymbol{\delta})-\boldsymbol{\mu_{x}}||^{2}=||\boldsymbol{\Sigma}\boldsymbol{\delta}||^{2}=\boldsymbol{\delta}^{\top}\boldsymbol{\Sigma}^{2}\boldsymbol{\delta}$.

If our goal then is to find a shift that we can estimate with a high degree of efficiency, we can solve the following:

for some constant $c$ defining the size of intervention. It is straightforward to show that this expression is minimized when $\boldsymbol{\delta}$ is chosen to be proportional to the eigenvector of $\boldsymbol{\Sigma}$ corresponding to its largest eigenvalue, $\lambda_{\max}$. Hence our finding demonstrates that for a fixed intervention strength (as measured by the Wasserstein distance), the most statistically efficient causal effect to estimate, at least in terms of this density ratio, is the one that corresponds to shifting the exposure mixture along its primary axis of variation. Note that this result only holds exactly under a multivariate normal distribution, but we have seen empirically that this choice of $\boldsymbol{\delta}$ typically leads to efficient estimates of $\psi(\boldsymbol{\delta})$ under a variety of exposure distributions.

A primary motivation for developing our framework is to answer policy-relevant questions, such as identifying the optimal way to modify an environmental exposure mixture to achieve the greatest public health benefit. This naturally leads to an optimization problem: finding the intervention vector $\boldsymbol{\delta}$ that maximizes (or minimizes) the causal estimand $\psi(\boldsymbol{\delta})$. However, a naive comparison across all possible $\boldsymbol{\delta}$ vectors is misleading. For any intervention direction that yields a beneficial effect, one could simply increase the intervention’s strength, for example, by scaling the magnitude of $\boldsymbol{\delta}$ to produce an arbitrarily larger (or smaller) value of $\psi(\boldsymbol{\delta})$. This would invariably lead to trivial solutions that correspond to extreme, unrealistic shifts in the exposure distribution. A more relevant policy question is “for a fixed amount of interventional effort, what is the best direction to apply that effort?”. To formalize this, we must first establish a fair basis for comparison, constraining our search to a set of interventions that are of the same size, which captures the actual movement of the exposure values, not just the magnitude of the parameter $\boldsymbol{\delta}$. Note that the notion of fairness here is with respect to the size of the intervention being applied, which differs from typical notions of fairness in causal inference or related fields, such as in recent work that defines a fairness criterion based on whether an estimand preserves the ordinal ranking of effects across all covariate subgroups.(McClean et al., 2024)

In order to establish our notion of fairness between interventions, we can use the 2-Wasserstein distance, $d_{W}(f,g_{\boldsymbol{\delta}})$, which is widely studied and serves this purpose (Panaretos and Zemel, 2019). Intuitively, the Wasserstein distance measures the minimum cost of transporting the probability mass of one distribution to match another, akin to the cost of moving a pile of dirt. By fixing the Wasserstein distance, we ensure that the total amount of shift between the old and new distributions is fixed. The optimization problem thus becomes a meaningful search for the $\boldsymbol{\delta}$ that minimizes $\psi(\boldsymbol{\delta})$ among all possible intervention directions of a comparable magnitude. This distributional fairness notion is widely used in both the operations and statistics research literatures. (Mohajerin Esfahani and Kuhn, 2018; Blanchet and Murthy, 2019; Duchi and Namkoong, 2021; Gao and Kleywegt, 2023)

For analytical tractability, we rely on a well-established result providing a formula for the squared 2-Wasserstein distance $d_{W}^{2}$ based on the means ($\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2}$) and covariance matrices ($\boldsymbol{\Sigma}_{1},\boldsymbol{\Sigma}_{2}$) of two distributions $(P_{1},P_{2})$, which is referred to as the Gelbrich formula (Gelbrich, 1990). Crucially, this formula serves as a general lower bound for the true squared 2-Wasserstein distance between any two probability measures on $\mathbb{R}^{q}$ with finite second moments. Specifically, we have that

Furthermore, this bound becomes an exact equality for any two distributions belonging to the same family of elliptically contoured distributions, a class which notably includes all multivariate normal distributions as well as uniform distributions on ellipsoids. The tractability of the Gelbrich formula has invited various researchers to use it as a surrogate for the 2-Wasserstein distance, a method frequently employed to overcome the computational complexity associated with the Wasserstein metric (Kuhn et al., 2019; Hakobyan and Yang, 2024), since the empirical 2-Wasserstein distance lacks a closed-form formula from data and can only be approached by numerical methods (Panaretos and Zemel, 2020). Moreover, the lower bound it provides has been shown to be tight in a fairly general situation against an upper bound derived for the 2-Wasserstein distance (Biswas and Mackey, 2024; Papp and Sherlock, 2024), and is therefore generally close to the 2-Wasserstein distance (Nguyen et al., 2021; Ye et al., 2024). Therefore, we believe our use of the Gelbrich formula as an approximation formula for the true 2-Wasserstein distance is reasonable.

This provides a computationally tractable and well-justified measure of intervention size. Accordingly, we measure intervention size through the Gelbrich formula applied to the marginal baseline and tilted laws of $\boldsymbol{W}$, and write the resulting quantity as $G(\boldsymbol{\delta})$. Comparing interventions on the level set $G(\boldsymbol{\delta})=c^{2}$ puts them on a common scale. Under a multivariate normal model this coincides with the squared 2-Wasserstein distance, while in more general settings it remains a rigorous lower bound. The optimization problem thus becomes a search for the $\boldsymbol{\delta}$ that minimizes $\psi(\boldsymbol{\delta})$ among all possible intervention directions of a comparable size, allowing us to disentangle the direction of an intervention from its size and enabling a principled exploration of which changes to an environmental exposure mixture are most beneficial or harmful.

Having established a principled method for comparing interventions of equivalent magnitude, we are now positioned to define our estimands of interest. The framework of incremental effects, when extended to a multivariate setting, moves beyond simply estimating the effect of a single, pre-specified shift or simply examining how the causal effect depends on shift size. Instead, it allows us to explore the entire space of potential interventions to identify those that are most impactful.

For a fixed intervention vector $\boldsymbol{\delta}$, the target estimand is the expected outcome under the tilted exposure distribution:

where the outer expectation is over the marginal distribution of covariates $\boldsymbol{X}$. Therefore, a direct approach to estimating $\psi(\boldsymbol{\delta})$ is through a plug-in procedure. This method involves replacing each component in the expression above with a corresponding empirical estimate, which leads to the plug-in estimator:

However, this estimator faces certain practical and theoretical challenges. One issue is that the estimator’s performance is highly dependent on an accurate estimate of the multivariate conditional density $\hat{f}(\boldsymbol{w}|\boldsymbol{x})$, which can be difficult to estimate for multivariate exposures, and will likely have slower convergence rates. Furthermore, this approach does not leverage the structure of the efficient influence function and, as a result, is generally not statistically efficient. These limitations motivate alternative estimators with superior statistical properties.

Our primary estimand, the optimal policy shift $\boldsymbol{\delta}^{*}_{c}$, is defined as the solution to

where $\mathcal{M}:=\{\boldsymbol{\delta}\in\mathbb{R}^{d}:G(\boldsymbol{\delta})=c^{2}\}$ and $G$ is the Gelbrich quantity defined in Section 2.4. We assume $c^{2}$ is a regular value of $G$ (equivalently, $\nabla G(\boldsymbol{\delta})\neq\boldsymbol{0}$ for all $\boldsymbol{\delta}\in\mathcal{M}$), so that $\mathcal{M}$ is a smooth embedded hypersurface. This regularity holds generically and is verified for the Gelbrich constraint in Appendix B. Standard Euclidean optimization algorithms such as gradient descent are not directly applicable as a gradient step taken from a point on the manifold will likely lead to a point outside of the feasible set.

To properly solve this problem, we must recognize that the constraint set forms a smooth, curved space known as a Riemannian manifold. Optimization on manifolds requires specialized techniques that generalize concepts from Euclidean optimization. The core idea is to perform optimization steps within the tangent space at each point on the manifold, which is a local linear approximation of the manifold, and then map the result back onto the manifold itself. For our specific constrained problem, we employ a Riemannian Broyden Fletcher Goldfarb Shanno (BFGS) algorithm. This is a powerful quasi-Newton method adapted for optimization on manifolds, which generally offers faster convergence than simpler first-order methods. In practice, the method is implemented with two computationally light ingredients: a projection-based retraction back to the level set, and an inexpensive projection-based vector transport between tangent spaces. The key components of this algorithm for our problem are:

1. 1.

   Preliminary definitions: A definition of the tangent space is to take any small smooth curve on $\mathcal{M}$

   $$\gamma:(-\varepsilon,\varepsilon)\rightarrow\mathcal{M},\quad\gamma(0)=\boldsymbol{\delta}.$$

   Then

   $$\gamma^{\prime}(0)\in\mathbb{R}^{d},$$

   is a valid tangent vector. The tangent space at a point $\boldsymbol{\delta}$ is the collection of all such vectors:

   $$T_{\delta}\mathcal{M}=\left\{\gamma^{\prime}(0):\gamma(0)=\delta,\gamma(t)\in\mathcal{M}\text{ for all }t\right\}$$

   For a level set $\mathcal{M}=\{\boldsymbol{\delta}:G(\boldsymbol{\delta})=c^{2}\}$, we also have the equivalent characterization $T_{\boldsymbol{\delta}}\mathcal{M}=\text{Null}(\nabla G(\boldsymbol{\delta})^{\top})$. The BFGS (Broyden–Fletcher–Goldfarb–Shanno) algorithm is a highly effective quasi-Newton method, which finds an optimum by maintaining and iteratively updating an approximation $\mathbf{B}$ of the inverse Hessian matrix. In Euclidean space $\mathbb{R}^{n}$, the standard update formula is:

   $$\mathbf{B}_{k+1}=\mathbf{B}_{k}+\frac{\mathbf{b}_{k}\mathbf{b}_{k}^{\top}}{\mathbf{b}_{k}^{\top}\mathbf{a}_{k}}-\frac{(\mathbf{B}_{k}\mathbf{a}_{k})(\mathbf{B}_{k}\mathbf{a}_{k})^{\top}}{\mathbf{a}_{k}^{\top}\mathbf{B}_{k}\mathbf{a}_{k}}$$

   Where:

   • •

     $\mathbf{B}_{k+1}$ is the updated approximation of the inverse Hessian, which is the target of the update.

   • •

     $\mathbf{B}_{k}$ is the current approximation of the inverse Hessian.

   • •

     $\mathbf{a}_{k}$ is the change in position (step), defined as $\mathbf{a}_{k}=\boldsymbol{\delta}_{k+1}-\boldsymbol{\delta}_{k}$.

   • •

     $\mathbf{b}_{k}$ is the change in gradient, defined as $\mathbf{b}_{k}=\nabla\psi(\boldsymbol{\delta}_{k+1})-\nabla\psi(\boldsymbol{\delta}_{k})$

2. 2.

   Riemannian Gradient: The standard (Euclidean) gradient, $\nabla\psi(\boldsymbol{\delta})$, must be projected onto the tangent space at the current point $\boldsymbol{\delta}$ to find the direction of steepest ascent along the manifold. Let $\nabla G(\boldsymbol{\delta})$ denote the Euclidean normal to the level set at $\boldsymbol{\delta}$, and let $\boldsymbol{n}(\boldsymbol{\delta})=\nabla G(\boldsymbol{\delta})/\|\nabla G(\boldsymbol{\delta})\|$ be the corresponding unit normal vector. The orthogonal projection matrix onto the tangent space is defined as $P_{\boldsymbol{\delta}}:=I-\boldsymbol{n}(\boldsymbol{\delta})\boldsymbol{n}(\boldsymbol{\delta})^{\top}$. Consequently, the Riemannian gradient under the induced Euclidean metric is given by

   $$\operatorname{grad}\psi(\boldsymbol{\delta})=P_{\boldsymbol{\delta}}\nabla\psi(\boldsymbol{\delta}).$$

3. 3.

   Retraction: Even though we are moving based on a Riemannian gradient step that is in the current point’s tangent space (which is essentially a linear approximation of the manifold at the current point), we are still technically moving out of the manifold. After finding a search direction $\boldsymbol{v}$, we need a retraction to map the point from the tangent space back onto the manifold. For the level-set manifold $\mathcal{M}$, we take the projection retraction $R_{\boldsymbol{\delta}}(\boldsymbol{\xi})=\Pi(\boldsymbol{\delta}+\boldsymbol{\xi})$, where $\Pi$ is the orthogonal projection onto $\mathcal{M}$, which is a canonical construction for embedded manifolds. In our code, $\Pi$ is computed approximately by a single normal-direction correction, which can be interpreted as a projection-like update that preserves the local first-order accuracy required of a retraction while substantially reducing computational cost.

4. 4.

   Vector Transport: Vector transport is required to compare tangent vectors that live in different tangent spaces across iterates. For our embedded level-set manifold, we use the simple projection transport:

   $$\widetilde{\mathcal{T}}_{\boldsymbol{\delta}\to\boldsymbol{\delta}_{+}}(\boldsymbol{\zeta}):=P_{\boldsymbol{\delta}_{+}}\boldsymbol{\zeta},$$

   namely, orthogonally projecting an ambient vector onto the new tangent space. Such a transport is generally not an isometry, even though it is often computationally attractive. This transport is also used to move tangent-space quantities between iterates when forming the RBFGS update.

To summarize, the iterative process of the Riemannian BFGS algorithm is outlined as follows.

When the objective is to minimize $\psi(\boldsymbol{\delta})$ over the Gelbrich level set, convergence guarantees for RBFGS depend on the line-search conditions, the regularity of the objective along retraction curves, and the choice of vector transport. Under the assumptions stated in Appendix B, the global convergence result in Theorem 4.2 of Huang et al. (2018) yields $\liminf_{k\to\infty}\|\operatorname{grad}\psi(\boldsymbol{\delta}_{k})\|=0$ for the RBFGS scheme equipped with a suitable isometric vector transport and weak Wolfe line search. Appendix B verifies the geometric regularity of the Gelbrich level set, establishes objective regularity along the retraction curve, and invokes the existence of a weak Wolfe step. The proof therefore uses a fully rigorous projection retraction and the scaled isometric vector transport $\mathcal{T}^{S}$ (Huang et al., 2015) together with an explicit boundary condition on $G$ that yields compactness of the relevant level set. Generally, however, these stricter choices are made for the convergence proof. In the actual implementation we use a computationally cheaper projection-type retraction and the orthogonal projection transport $\widetilde{\mathcal{T}}$, together with safeguard checks on the RBFGS update. This separation between the proof device and the working algorithm is standard in large-scale Riemannian optimization (Ring and Wirth, 2012; Huang et al., 2018). In our numerical work, the simplified implementation converged stably to a local solution while substantially reducing computation time.

In embedded-manifold settings, the projection transport $\boldsymbol{\zeta}\mapsto P_{\boldsymbol{\delta}_{+}}\boldsymbol{\zeta}$ can significantly reduce computational time without degrading practical convergence behavior. This phenomenon is well documented in the RBFGS experiments of Qi et al. (2010), where projection-based transport on the Stiefel manifold reduces wall-clock time from $304.0$ seconds to $24.0$ seconds and also decreases the iteration count from $175$ to $83$ in the Procrustes problem. Moreover, Huang et al. (2018) compare combinations that satisfy sufficient conditions used in global convergence theory with cheaper combinations that do not, and find that the numbers of function and gradient evaluations (as well as vector transports) are not significantly affected by the choice of retraction and transport; as a result, lower-complexity retractions and transports can be markedly faster in terms of computational time. In our applications, the simplified implementation described above converges reliably and is computationally fast.

In this section, we establish a minimax lower bound for the incremental effect

under a flexible nonparametric model. Minimax lower bounds benchmark what is statistically achievable without imposing additional structure: they show that no estimator can attain uniformly smaller risk over a given model class. Before presenting the main results, we must first define a number of important terms. First, we refer to the centered version of the exposures as

To analyze the asymptotic variance of the effect difference $\theta(\boldsymbol{\delta})$, we decompose the problem into an incremental component $\boldsymbol{h}(\boldsymbol{Z})$ and a baseline influence function $\varphi_{0}(\boldsymbol{Z})$. We define

Next, we calculate the covariance structure of $\boldsymbol{h}(\boldsymbol{Z})$. A key part of this variance arises from the outcome regression. To ensure this term represents a proper covariance (centered moment), we compute the raw second moment and subtract the outer product of the mean:

Then, letting $\boldsymbol{\Sigma}_{\tilde{\boldsymbol{W}}\tilde{\boldsymbol{W}}^{\top},\,\varepsilon}:=\operatorname{E}\big[\operatorname{Var}(Y\mid\boldsymbol{X},\boldsymbol{W})\,\tilde{\boldsymbol{W}}\tilde{\boldsymbol{W}}^{\top}\big]$ be the residual variance component, we show in the appendix that we can write the full covariance matrix of $\boldsymbol{h}(\boldsymbol{Z})$ as

Finally, estimating the incremental effect $\theta(\boldsymbol{\delta})$ requires distinguishing the variation in the shift direction from the variation in the baseline. Since the gradient component $\boldsymbol{h}(\boldsymbol{Z})$ and the baseline influence $\varphi_{0}(\boldsymbol{Z})$ are typically correlated, the fundamental difficulty of estimating their difference is governed by the variability of $\boldsymbol{h}$ that cannot be explained by $\varphi_{0}$. Mathematically, this corresponds to the residual variance of $\boldsymbol{h}$ after linearly projecting it onto the space spanned by $\varphi_{0}$, yielding the Schur–complement covariance

The next result records how the efficiency bound depends on $\boldsymbol{\delta}$.

Lemma 1 makes explicit that the nonparametric efficiency bound (the variance of the efficient influence function) grows quadratically in the tilt magnitude, with geometry determined by the matrix $\boldsymbol{H}$; this geometry blends the noise component $\boldsymbol{\Sigma}_{\tilde{\boldsymbol{W}}\tilde{\boldsymbol{W}}^{\top},\,\varepsilon}$ and the outcome–regression component $\boldsymbol{\Sigma}_{\mu,\mathrm{full}}$. Since $\boldsymbol{H}=\boldsymbol{\Sigma}_{\tilde{\boldsymbol{W}}\tilde{\boldsymbol{W}}^{\top},\,\varepsilon}+\boldsymbol{\Sigma}_{\mu,\mathrm{full}}\succeq\boldsymbol{\Sigma}_{\tilde{\boldsymbol{W}}\tilde{\boldsymbol{W}}^{\top},\,\varepsilon}$, the lower bound immediately implies

As the direction of $\boldsymbol{\delta}$ changes, the bound can increase substantially showing how the most efficient shifts are those proportional to the first eigenvector of $\boldsymbol{\Sigma}_{\tilde{\boldsymbol{W}}\tilde{\boldsymbol{W}}^{\top},\,\varepsilon}$. Note that this matrix is not the same as the covariance matrix of the exposures, but is closely related to it, and therefore this result confirms our findings of Section 2 about which directions are most efficient to estimate. Now, we derive minimax bounds for the estimation error of $\theta(\boldsymbol{\delta})$.

Theorem 2 shows that the best possible root mean–squared error obeys

revealing an effective sample size of order $n_{\mathrm{eff}}\asymp n/(\boldsymbol{\delta}^{\top}\boldsymbol{\Gamma}\boldsymbol{\delta})$. Hence larger tilts and directions aligned with high–variance components of $\boldsymbol{\Gamma}$ are intrinsically harder, regardless of the estimation strategy. Conversely, if $\boldsymbol{\delta}$ lies near a low–variance direction of $\boldsymbol{\Gamma}$, faster convergence is attainable. Together with Lemma 1, the theorem clarifies why one must account for $\boldsymbol{\delta}$ when assessing estimation difficulty and designing procedures: the geometry induced by $\boldsymbol{\Gamma}$ determines how precision scales with both the size and the direction of the tilt.

Under mild regularity conditions (boundedness, i.i.d. sampling, and a fixed finite tilt), we establish a finite-$\boldsymbol{\delta}$ central limit theorem for our cross-fitted one-step estimator. At first sight, the efficient influence function suggests that asymptotic linearity requires controlling the $\boldsymbol{\delta}$-specific nuisance functions

Since both objects are indexed by the tilt parameter $\boldsymbol{\delta}$, and $m_{\boldsymbol{\delta}}$ is a nonlinear functional of the observed-data law, it is not obvious how to guarantee the $L_{2}$ rates needed for one-step asymptotics.

Our analysis shows that, for any fixed finite tilt $\boldsymbol{\delta}$, it is in fact sufficient to estimate only two familiar observed-data nuisance components:

namely the outcome regression and the (generalized) propensity score. In the Appendix we show that, under bounded support for $\boldsymbol{W}$ and the finite-tilt condition $\|\boldsymbol{\delta}\|\leq\Delta$, the maps $(\mu,f)\mapsto r_{\boldsymbol{\delta}}$ and $(\mu,f)\mapsto m_{\boldsymbol{\delta}}$ are Lipschitz in $L_{2}(P)$: there exist fixed finite constants $C_{1}(\Delta),C_{2}(\Delta)$ such that for any estimators $(\widehat{\mu},\widehat{f})$,

where $(\widehat{m}_{\boldsymbol{\delta}},\widehat{r}_{\boldsymbol{\delta}})$ are obtained from $(\widehat{\mu},\widehat{f})$ by the same formulas as above; see Lemma 9. Consequently, the product condition

which ensures that the second-order remainder is negligible, is implied by the more interpretable requirement that both $\widehat{\mu}$ and $\widehat{f}$ converge at rate $n^{-1/4}$ in $L_{2}(P)$; see Corollary 3. These $L_{2}$ rates for $(\widehat{\mu},\widehat{f})$, and hence for $(\widehat{m}_{\boldsymbol{\delta}},\widehat{r}_{\boldsymbol{\delta}})$, can be guaranteed under standard smoothness and complexity conditions by the highly adaptive lasso (van der Laan, 2017), which attains near-optimal convergence rates for a broad nonparametric class.

Theorem 3 provides the basis for asymptotic inference, which proceeds in a standard influence-function manner. Let $\boldsymbol{Z}_{i}=(\boldsymbol{X}_{i},\boldsymbol{W}_{i},Y_{i})$, and define the empirical influence values for $\psi(\boldsymbol{\delta})$ by

where $\widehat{r}_{\boldsymbol{\delta}}$ and $\widehat{m}_{\boldsymbol{\delta}}$ are constructed from the cross-fitted estimators $(\widehat{\mu},\widehat{f})$ as above. A consistent estimator of the asymptotic variance is the sample second moment

and analogously for the incremental effect we use

These yield Wald-type confidence intervals of the form

with a similar process for $\psi(\boldsymbol{\delta})$. When inference is required for several tilts $\boldsymbol{\delta}_{1},\ldots,\boldsymbol{\delta}_{J}$, a joint covariance estimator is obtained by the empirical covariance matrix of the vector

which enables simultaneous confidence intervals or Wald tests via a multivariate normal approximation.

Let $f(\boldsymbol{w}\mid\boldsymbol{V})$ denote the true conditional exposure density given the full adjustment set. The parameter $\psi(\boldsymbol{\delta})$ can be written as a linear functional of the long regression evaluated under the observed distribution of $(\boldsymbol{V},\boldsymbol{W})$:

Analogously, the short estimand admits the representation

Under mild regularity conditions ensuring square-integrability of the relevant Riesz representers, Chernozhukov et al. (2022) shows that the bias admits the exact representation

where

This identity isolates two necessary sources of bias: the additional outcome variation explained by $\boldsymbol{U}$ beyond $(\boldsymbol{X},\boldsymbol{W})$ through $\Delta_{\mu}$, and the additional information about the exposure distribution provided by $\boldsymbol{U}$ beyond $\boldsymbol{X}$ through $\Delta_{\alpha}$. By Cauchy–Schwarz,

To express the bound in terms of identifiable scale components and partial $R^{2}$-type sensitivity parameters, define the following identifiable parameters

We parameterize the outcome component by the nonparametric partial $R^{2}$ of $\boldsymbol{U}$ with $Y$ given $(\boldsymbol{X},\boldsymbol{W})$,

For the treatment component, we use the relative increase in $L_{2}$ variation of the Riesz representer induced by conditioning on $\boldsymbol{U}$:

Combining these definitions yields the sensitivity bound