Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
103,482 characters · 17 sections · 78 citation commands
Inference on Optimal Policy Values and Other Irregular Functionals via Softmax Smoothing
Individualized treatment policies, or mappings from a unit's observed characteristics to an assigned treatment, are deployed in a wide range of decision-making problems. In medicine, individualized policies may be used to maximize the expected outcome or well-being of a patient qian2011performance, xu2022estimating. In advertising, such policies may be employed for maximizing patron engagement or revenue bottou2013counterfactual. Motivated by such applications, a large literature has developed that focuses on using observational data to learn a treatment policy with a high value, i.e.\ a large outcome averaged over a covariate distribution athey2021policy,kitagawa2018should, kitagawa2022treatment, hirano2009asymptotics.
Before investing resources in the construction and deployment of a highly individualized treatment rule, one may first wish to quantify how much can be gained from policy learning. The natural estimand for assessing this is the value of the optimal treatment policy. Let $Z = (X, A, Y)$ denote a observation, where $X$ represents a unit's covariates, $A \in [N]$ denotes a categorical treatment, $Y(a)$ a unit's potential outcome associated with treatment $a \in [N]$, and $Y = Y(A)$ the observed outcome. Letting $\pi(X) \in [N]$ denote a arbitrary policy that assigns a unit to a treatment, the value of the optimal treatment policy is given by \[ V^\ast := \max_\pi V(\pi), \;\; \text{where} \;\; V(\pi) := \mathbb{E}[Y(\pi(X))]. \] Inference on $V^\ast$ can be used to determine the value of developing individualized policies in several ways. For instance, if an advertiser wanted to determine if some baseline treatment policy $\pi_0$ (perhaps a rule-based policy) could be improved through the use of machine learning (ML) algorithms, he could estimate the policy value gap $V^\ast - V(\pi_0)$. Likewise, if a teacher wanted to determine if a single lesson plan was a good fit for all students, she could estimate the value of personalization, i.e.\ the quantity $V^\ast - \max_{a \in [N]}V(\pi_a)$, where $\pi_a(X) = a$ denotes a constant treatment policy. These examples motivate the need for simple, computationally efficient inferential procedures for $V^\ast$.
At first glance, because $V^\ast$ is just a scalar value, it may appear that estimating $V^\ast$ should be easier than learning a good individualized treatment policy $\pi(X)$. However, because the target estimand $V^\ast$ is irregular, standard approaches for constructing confidence intervals fail to give valid coverage. In more detail, under standard causal assumptions, $V^\ast$ can be re-written in terms of the outcome regression (i.e.\ Q-function) $Q^\ast(a, x) = \mathbb{E}[Y \mid X = x, A = a]$ as
When the argument maximizing set $\arg\max_{\ell \in [N]}Q^\ast(\ell, X)$ is non-unique with positive probability, the functional $Q \mapsto \mathbb{E}\left[\max_{\ell} Q(\ell, X)\right]$ is not Gateaux/path-wise differentiable at $Q^\ast$. Classical doubly-robust and targeted maximum likelihood estimators require the Gateaux differentiability of the underlying functional, and are thus inapplicable when multiple treatments offer the same expected reward. Restricting ourselves to the case of binary treatments $A \in \{0, 1\}$ and letting $\tau^\ast(X) := Q^\ast(1, X) - Q^\ast(0, X)$ denote the conditional average treatment effect (i.e.\ CATE), these classical estimation approaches are inapplicable precisely when $\mathbb{P}(\tau^\ast(X) = 0) > 0$, i.e.\ when there is a positive probability treatment non-response in the population.\footnote{Throughout this paper, when we say “non-response to treatment,” we mean that multiple treatments yield the same expected outcome. A stronger notion of non-response would be that the potential outcomes are identical, i.e.\ that $Y(1) = Y(0)$ with positive probability in the binary treatment setting.}
This irregularity has forced researchers to develop estimators that are tailor-made for estimating $V^\ast$. Some works assume that the optimal action is almost surely unique semenova2023aggregated and others make parametric assumptions on the outcome regression (e.g.\ linearity) to enable inference laber2014dynamic, goldberg2014comment, chakraborty2010inference. While these assumptions may be reasonable in some settings, they can be violated in domains such as medicine or advertising where treatments can be entirely ineffective and response curves may be highly non-parametric. The most flexible estimation approaches, due to luedtke2016statistical and shi2020breaking, allow for a positive probability of non-response to treatment and also permit non-parametric outcome regressions. This flexibility comes at the cost of refitting a number of nuisance models that grows with the sample size. Refitting nuisances many times may be possible when nuisance models are parametric or datasets are small. However, when datasets are large or black-box learners such as neural networks are used, training a growing number of nuisance models may not be feasible. We discuss these approaches amongst others in greater detail in the related work section (Section (ref)). There is thus a need for computationally efficient estimators that allow non-parametric outcome regressions and non-response to treatment.
In this paper, we show that careful smoothing can be used to construct asymptotically normal estimates for both the value of the optimal treatment policy and also a broader class of maximum-type estimands. Defining the softmax function $\mathbf{sm}^\beta : \mathbb{R}^N \rightarrow \mathbb{R}$ by
our approach is to replace the target $V^\ast$ outlined in Equation (ref) with a smoothed surrogate: \[ \underbrace{V^\ast := \mathbb{E}\left[\max_\ell\{Q^\ast(\ell, X)\}\right]}_{\text{Non-differentiable functional}} \quad\Longrightarrow\quad \underbrace{V^\beta := \mathbb{E}\left[\mathbf{sm}^\beta_\ell\{Q^\ast(\ell, X)\}\right]}_{\text{Smoothed approximation}}. \] The functional $Q \mapsto \mathbb{E}\left[\mathbf{sm}^\beta_\ell Q(\ell, X)\right]$ is twice Gateaux differentiable, thus enabling the use of semi-parametric de-biasing techniques for performing inference on $V^\beta$. We enumerate our contributions in greater detail below.
Two theoretical insights drive our analysis. Our first insight is that the particular choice of smoothing function is critically important in the estimation of irregular parameters. Several works have historically leveraged the softplus function $\mathbf{sp}^\beta\{u_1, \dots, u_N\} := \frac{1}{\beta}\log\left(\sum_{i = 1}^N e^{\beta u_i}\right)$ to construct smooth surrogates for non-differentiable functions levis2023covariate, goldberg2014comment. This smoother possesses a notable defect: when $u_1, \dots, u_N$ are equal, the soft-plus function differs from the true maximum by exactly $\log(N)/\beta$. In general, the magnitude of this bias is too large to permit asymptotically normal inference around the un-smoothed policy value $V^\ast$. In contrast, the softmax exhibits zero bias in presence of ties, which will enable inference even when argument maximizing set contains multiple elements.
Our second insight is that the softmax function exhibits fast bias decay under mild distributional assumptions on the sub-optimality gaps $ \Delta_k$. If these gaps admit a density bounded in the sense of Informal Theorem (ref), the bias of softmax smoothing will decay as $|V^\ast - V^\beta| = O\left(\beta^{-(1 + \delta)}\right)$ (Lemma (ref)). The provides a window in which one can select $\beta_n$ to obtain valid inference. This density assumption is analogous to the soft-margin condition considered in luedtke2016statistical and shi2020breaking, and thus allows for an arbitrary large probability of treatment non-response.\footnote{The classical soft-margin condition states that $\mathbb{P}(0 < \Delta_k \leq t) \lesssim t^{\delta}$ for $t \in (0, t_0]$ for some $\delta, t_0 > 0$.} When $\delta = 1$ and there are binary treatments, this is equivalent to the CATE have bounded density near zero.
Lastly, we note that the softmax function described above is commonly used in the machine learning literature for tasks such as soft Q-learning schulman2017equivalence, nachum2017bridging, haarnoja2017reinforcement,haarnoja2018soft. However, the goal in this literature is not inference, but rather regret minimization, for which a naive bias bound of $O(1/\beta)$ actually suffices. For inferential tasks, naive bias bounds generally lead to pessimistic confidence intervals for smoothed estimates levis2023covariate, zhang2024winners. Our work indicates that, with appropriate care, tools that are commonly leveraged in machine learning tasks may be applicable to broad classes of statistical problems.
\paragraph{Inference on Optimal Policy Values:} Many works have focused on performing inference on optimal policy values under parametric restrictions on outcome regressions/Q-functions. laber2014dynamic assume linear Q-functions and propose estimating the optimal policy value in a dynamic regime using Q-learning. Instead of smoothing, they construct pessimistic confidence intervals in regions of non-regularity (i.e.\ regions of non-responders) via a union bound and bootstrapping methods. goldberg2014comment consider the same parametric restrictions and show how to use a softplus function to approximate $\max\{a, b\}$. These authors also require conservative estimation when the probability of non-response is positive. We emphasize that our estimator is asymptotically normal even when non-response to treatment occurs, and that our results hold in both non-parametric and semi-parametric settings. chakraborty2013inference propose a method based on $m$-out-of-$n$ bootstrapping to build confidence intervals on structural parameters, but assume Q-functions are linear and achieve width $m$ confidence intervals instead of width $n$ (here, $m = o(n)$).
Other authors have considered settings where the Q-functions/conditional average treatment effects (CATEs) are allowed to be fully non-parametric. For instance, semenova2023aggregated analyze a first-order de-biased estimator estimator for the optimal policy value. However, this estimator requires that nuisances be estimated at $o(n^{-1/4})$ rates in the $L^\infty$ norm, which may be restrictive when off-the-shelf ML estimators like gradient boosted trees or neural networks are used. Further, they require the reward maximizing action to be unique almost surely. We also mention the work of park2024debiased\footnote{This cited paper was released after the first draft of this paper chen2023inference}, who use smoothing but obtain slower, $\sqrt{n/\beta_n^2}$ rates of convergence, where $\beta_n$ grows to infinity (Proposition 1 of their work). In contrast, our Theorem (ref) establishes $\sqrt{n}$ asymptotic linearity of the softmax smoothed estimator, which results in tighter confidence intervals. Our results also apply to other irregular functionals of interest.
The two works most closely related to the present paper in are those of shi2020breaking and luedtke2016statistical, which are discussed above. These works are very general, describing estimators that allow for non-parametric outcome regressions, dynamic treatment regimes, and positive probabilities of treatment non-response. However, both approaches require training a number of nuisance models that grows with the number of samples. Our estimator requires only fitting a fixed number of nuisances, but at the cost of selecting a smoothing parameter. There is thus a fundamental tradeoff between our estimator and the two mentioned above. Because these methods fit a growing number of nuisance models, they require stronger nuisance error control than our estimator. In our setting, it suffices that the $L^2$ nuisance errors vanish in probability, which is the standard in semi-parametric literature chernozhukov2018double. On the other hand, these works require either $L^2$ convergence of the nuisance estimates or high probability bounds on the $L^2$ error (see bibaut2020sufficient for a discussion of this in the context of the estimator of luedtke2016statistical).
\paragraph{Policy Learning:}
In the causal inference and reinforcement learning (RL) literature, a large body of work has developed related to learning an optimal or near-optimal individualized treatment policy. These include approaches based on Q-learning zhang2012robust, watkins1992q, watkins1989learning, qian2011performance,zhao2009reinforcement, zhao2011reinforcement, moodie2012q, A-learning (or advantage learning) shi2024statistically, murphy2003optimal, moodie2007demystifying, robins2004optimal, empirical risk minimization athey2021policy, luedtke2020performance, kallus2018confounding, and more. We particularly emphasize works relating to empirical risk minimization athey2021policy, kitagawa2018should, kitagawa2022treatment, which describe how to learn policies with low regret, which is defined as the difference in value between the best policy in the class and the learned policy manski2004statistical, hirano2009asymptotics, luedtke2020performance. We also note works on optimal policy learning in partially-identified settings olea2023decision, kitagawa2023treatment, pu2021estimating, which may be of particular importance when one does not believe commonly-made conditional exogeneity assumptions. We emphasize that the algorithms in the works listed above do not discuss inference on the value of the optimal treatment policy. Also relevant to our work are algorithms from the RL literature based on soft-Q learning garg2021iq, nachum2017bridging, schulman2017equivalence, which use entropy regularization to learn a near-optimal Boltzmann policy. Again, these algorithms do not provide a means for a learner to perform inference on optimal policy values.
\paragraph{Causal Inference:}
Lastly, we touch on tools from the causal inference and semi-parametric statistics that are relevant to our work, including work on semi-parametric efficiency and estimation kosorok2008introduction, newey1994asymptotic, van2000asymptotic,robins1995analysis,robins1995semiparametric,bickel1993efficient, levit1976efficiency, missing data and doubly-robust estimation robins2004optimal,bang2005doubly,tsiatis2006semiparametric, and targeted maximum likelihood estimation van2006targeted, van2011targeted. Most closely related to our work is the literature of double/de-biased machine learning bang2005doubly, chernozhukov2018double, which advocates performing generic first-order corrections and cross-fitting to reduce bias in estimation problems involving nuisance components. These approaches are commonly leveraged in both causal inference wager2018estimation, athey2019generalized,runge2023causal, causal machine learning foster2023orthogonal, whitehouse2024orthogonal, lan2025meta, van2007super, and automatic de-biased learning methods hirshberg2021augmented, chernozhukov2022automatic, chernozhukov2022riesznet, chernozhukov2022nested.
We let $P_V$ be distribution of a arbitrary random variable $V$. We let $\mathbb{E}[\cdots]$ and $\mathbb{P}(\cdots)$ denote expectation and probability respectively over all sources of randomness. For arbitrary independent random variables $U$ and $V$ with distributions $P_U$ and $P_V$ and an arbitrary function $m(u, v)$, we let $\mathbb{E}_V[m(U, V)] := \int m(U, v)P_V(d v) = \mathbb{E}_{U,V}[m(U, V) \mid U]$. In this case $\mathbb{E}[m(U, V)] = \mathbb{E}_{U, V}[m(U, V)] = \int\int m(u, v) P_U(du) P_V(dv)$. Given samples $Z_1, \dots, Z_n$, we let $\mathbb{P}_n := \frac{1}{n}\sum_{m = 1}^n \delta_{Z_m}$ denote the corresponding empirical distribution, $\mathbb{P}_n f(Z) := \frac{1}{n}\sum_m f(Z_m)$, and $\mathbb{G}_n := \sqrt{n}(\mathbb{P}_n - \mathbb{E}_Z)$.
For a given integer $N > 0$, we let $[N] := \{1, 2, \dots, N\}$. For $p \in [1, \infty)$, a random variable $V \in \mathcal{V}$ with distribution $P_V$, and a function $f : \mathcal{V} \rightarrow \mathbb{R}^d$, we let $\|f\|_{L^p(P_V)} := \left(\mathbb{E}_V\|f(V)\|_p^p\right)^{1/p}$, where $\|a\|_p := \left(\sum_i |a_i|^p\right)^{1/p}$ for $a \in \mathbb{R}^d$. For $p = \infty$, we let $\|f\|_{L^\infty(P_V)} := \operatorname*{ess\,sup} |f(V)|$. We define the $L^p$ space of functions by $L^p(P_V) := \left\{f : \mathcal{V} \rightarrow \mathbb{R}^d \text{ s.t. } \|f\|_{L^p(P_V)} < \infty\right\}$. The dimension $d$ of the range of functions will either be clear from context or made explicit via the notation $L^p(P_V; \mathbb{R}^d)$.
Given a Banach space $B$, $T : B \rightarrow \mathbb{R}^d$, and $f^\ast, g \in B$, we define the Gateaux derivative of $T$ at $f^\ast$ in the direction $g$ as $D_f T(f^\ast)(g) := \frac{\partial}{\partial t}T(f^\ast + tg) \vert_{t = 0}$. We define the second Gateaux derivative as $D_{f}^2T(f^\ast)(g) := \frac{\partial^2}{\partial t^2}T(f^\ast + tg)\vert_{t = 0}$. If we have a bi-variate map $T : B\times B \rightarrow \mathbb{R}^d$ and $f_1^\ast, f_2^\ast, g_1, g_2 \in B$, we let the cross Gateaux derivative be defined as $D_{f_1, f_2}T(f_1^\ast, f_2^\ast)(g_1, g_2) := \frac{\partial^2}{\partial s \partial t}T(f_1^\ast + tg_1, f_2^\ast + s g_2)\vert{t = s = 0}$. For $f : \mathbb{R}^d \rightarrow \mathbb{R}$, we let $\partial_{x_i} f(x)$ denote the $i$th partial derivative, $\nabla_x f(x)$ the gradient, and $\nabla_x^2 f(x)$ the Hessian, when said objects exist. If $f : \mathbb{R}^d \rightarrow \mathbb{R}^m$, we will let $\partial_x f(x)$ denote the Jacobian. When denote the derivative of composed functions with respect to the argument of the outer function, we use the notation $\partial_i f(g(x)) := \partial_{u_i} f(u) \vert_{u = g(x)}$, $\nabla_u f(g(x)) := \nabla_u f(u)\vert_{u = g(x)}$, and $\nabla_u^2 f(g(x)) = \nabla_u^2 f(u)\vert_{u = g(x)}$ for convenience unless otherwise noted. Lastly, for functions $f, g : \mathcal{V} \rightarrow \mathbb{R}^{d}$, we define the bracket between $f$ and $g$ as $[f, g] := \{h : \mathcal{Z} \rightarrow \mathbb{R} : \exists \lambda : \mathcal{X} \rightarrow [0, 1] \text{ s.t. } h(z) = \lambda(z)f(z) + (1 - \lambda(z))g(z)\}$.
In this section, we formalize the problem of performing inference on the value of the optimal treatment policy. We present our assumptions on the underlying data-generating process and identify the optimal value both in terms of outcome regressions and blip effects. The latter identification will be used for performing inference under semi-parametric restrictions (see Section (ref)). In the main body of the paper we just focus on static treatment regimes, i.e.\ settings where there is a single round of treatment. We generalize this discussion to dynamic (namely, two-stage) treatment regimes in Appendix (ref). We conclude by discussing softmax smoothing, the main density assumption made throughout this work, and the bias introduced by smoothing. To streamline our presentation, we defer discussion of the case of general irregular functionals to Section (ref).
We assume observations are of the form $Z = (X, A, Y)$, where $X \in \mathcal{X}$ represents an individual's covariates, $A \in [N] := \{1, 2, \dots, N\}$ represents a treatment coming from a fixed set of options, and $Y \in \mathbb{R}$ represents an observed outcome. We make the following assumptions on the data generating process.
In what follows, we do not assume that the propensity $p^\ast(a \mid x)$ is known. By a treatment policy, we mean a mapping $\pi : \mathcal{X} \rightarrow [N]$ taking an individual's covariates to a categorical treatment. Under the above assumption, straightforward manipulation yields that $V^\ast := \max_{\pi}\mathbb{E}[Y(\pi(X))]$ is identified as $V^\ast = \mathbb{E}\left[\max_\ell\{Q^\ast(\ell, X)\}\right]$, which is precisely the characterization provided in Equation (ref).
We can also characterize the value of the optimal treatment policy via blip effects. The blip effect $\gamma^\ast(a, x) := Q^\ast(a, x) - Q^\ast(a^\ast, x)$ gives the expected value of an action $a \in [N]$ for an individual with covariates $x \in \mathcal{X}$ relative to some baseline or control action $a^\ast \in [N]$. Given that one can write the value of a policy $\pi$ relative to the observational policy as
it follows that the value of an optimal treatment policy can also be characterized as
We use this identification in Section (ref), in which we discuss the estimation of the optimal policy value under semi-parametric restrictions on the blip effects.
As noted in the introduction, the mapping $Q \mapsto \mathbb{E}\left[\max_{\ell}Q(\ell, X)\right]$ is not generally Gateaux differentiable at $Q^\ast$ when there are non-responders to treatment, i.e.\ when \[ \mathbb{P}\left(\left|\arg\max_{\ell \in [N]}Q^\ast(\ell, X)\right| > 1\right) > 0. \] In this setting, such laws are generally termed to be irregular or nonregular shi2020breaking, luedtke2016statistical. When the law is irregular, hirano2012impossibility show that regular asymptotically linear (RAL) estimators fail to exist in general. Consequently, standard notions of semi-parametric efficiency may fail to be applicable. Likewise, de-biased or one-step estimators chernozhukov2018double, van2006targeted rely on the Gateaux differentiability of the underlying functional, and thus cannot be used.
To work around this irregularity, we defined the smoothed surrogate $V^\beta := \mathbb{E}\left[\mathbf{sm}^\beta_\ell\{Q^\ast(\ell, X)\}\right]$. Since the softmax function $\mathbf{sm}^\beta\{u\}$ is twice continuously differentiable with respect to $u$, the above surrogate is twice Gateaux differentiable with respect to $Q$ for any fixed $\beta > 0$. This enables the estimation of $V^\beta$ via standard de-biased estimators chernozhukov2018double, van2006targeted. We ultimately want to obtain root-$n$ inference for $V^\ast$ through a de-biased estimator for $V^\beta$. To do this, we need to be able to select the smoothing parameter $\beta$ as a function of the sample size in a way to ensure the bias decays like $|V^\beta - V^\ast| = o(n^{-1/2})$. The bias of softmax smoothing is governed by the distribution of the sub-optimality gaps \[ \Delta_k := \max_\ell\{Q^\ast(\ell, X)\} - Q^\ast(k, X), \] which are small when action $k$ is close to optimal, and large otherwise. Without further assumptions on these gaps, the worst-case bias of smoothing decays as $|V^\beta - V^\ast| = \Theta\left(\tfrac{1}{\beta}\right)$. In this regime, to ensure $|V^\beta - V^\ast| = o(n^{-1/2})$, one would need to take $\beta = \omega(n^{1/2})$. In the sequel, we will see that taking $\beta$ to be this large would require the learner to estimate the unknown regression $Q^\ast$ at super-parametric $o_\mathbb{P}(n^{-1/2})$ rates, which is generally impossible.
Our insight is that, under a mild assumption of the density of the gaps $\Delta_k$, the bias will shrink rapidly as the smoothing parameter grows. Our main assumption formally codifies this assumption on the sub-optimality gaps, and is analogous to the assumptions considered in luedtke2016statistical and shi2020breaking (which are instead stated in terms of the CDF of the sub-optimality gaps).
When Assumption (ref) is applied to the $\Delta_k$, it allows the maximizing index of $Q^\ast(X) = (Q^\ast(1, X), \dots, Q^\ast(N, X))$ to be non-unique with positive probability. The parameter $\delta > 0$ governs the density of near ties: $\delta > 1$ indicates the density of $\Delta_k$ is small near zero, $\delta = 1$ corresponds to a bounded density, and $\delta \in (0, 1)$ indicates the $\Delta_k$ may be heavy-tailed. These cases are illustrated in Figure (ref). Under this assumption, we can prove a bound on the bias of softmax smoothing, which we prove in Appendix (ref).
To gain visual intuition for the above lemma, we can look at Panel (b) of Figure (ref). This panel illustrates a softmax smoothing approximation for the function $\max\{0, x\}$. From the figure, it is clear that the bias from smoothing is most pronounced near (but not at) zero. Further, as $\beta$ grows, the region in which the bias is pronounced shrinks. Thus, if the gaps are unlikely to be very close to the origin (in the sense of Assumption (ref)), we should be able to get faster than the naive $O\left(\tfrac{1}{\beta}\right)$ bias control.
In this section, we propose a de-biased estimator for the value of the optimal treatment policy based on softmax smoothing. Since $Q \mapsto \mathbb{E}\left[\max_\ell\{Q(\ell, X)\}\right]$ was not generally Gateaux differentiable, we proposed replacing it by a smoothed alternative $V^\beta := \mathbb{E}\big[\mathbf{sm}^\beta_\ell\{Q^\ast(\ell, X)\}\big]$. Our first goal is to provide a de-biased or Neyman orthogonal score for $V^\beta$ chernozhukov2018double.
Typically, $\theta_0$ is the unique solution to the moment equation $0 = \mathbb{E}\left[m(Z; \theta_0, g^\ast)\right]$, and one regularly has $g^\ast, \omega \in \mathcal{G} := L^2(P_W; \mathbb{R}^p)$, where $p \geq 1$ and $W \subset Z$ denotes a subset of features. We will always assume this to be the case unless otherwise stated. If we are interested in a parameter $\theta_0 = \mathbb{E}[m(Z; g^\ast)]$ instead, orthogonality reduces to the condition that $D_g\mathbb{E}[m(Z; g^\ast)](\omega) = 0$. Heuristically, a score is Neyman orthogonal if any errors in nuisance estimation are only second order in nature. A standard example of a Neyman orthogonal score is the doubly-robust pseudo-outcome leveraged in the estimation of average treatment effects kennedy2016semiparametric, kennedy2024semiparametric. The following proposition describes a Neyman orthogonal score for $V^\beta$ for any fixed $\beta > 0$. We provide a proof in Appendix (ref).
In general, the solution $V^\beta$ to the smoothed, de-biased score is not the same as the optimal policy value $V^\ast$. However, by appropriately selecting an increasing sequence of smoothing parameters $(\beta_n)_{n \geq 1}$ and exploiting the bias properties outlined in Lemma (ref), one can use the above score to obtain asymptotically linear (and hence normal) estimates of $V^\ast$. The following theorem describes an estimator for $V^\ast$ in which nuisance estimates are independent of the data. This sample-splitting approach is just for ease of exposition, and one can show that a natural cross-fitting variant of Theorem (ref) is true as well (see Remark (ref) below).
The above theorem provides an asymptotically valid range in which the learner can select a smoothing parameter. On one hand, we must smooth sufficiently aggressively to make the bias negligible. This is captured through the requirement that $\beta_n = \omega\left(n^{\frac{1}{2(1 + \delta)}}\right)$. On the other hand, $\beta_n$ must grow sufficiently slowly to ensure the estimating the unknown regression is feasible. This emerges through the condition that $\|\widehat{Q} - Q^\ast\|_{L^2(P_W)} = o_\mathbb{P}(\beta_n^{-1/2}n^{-1/4})$. This places an implicit upper bound of $\beta_n = o(n^{1/2})$ on the smoothing parameter: if $\beta_n = \Omega(n^{1/2})$, then we would need $\|\widehat{Q} - Q^\ast\|_{L^2(P_W)} = o_\mathbb{P}(n^{-1/2})$, which is generally not possible. We note that such “faster than $n^{-1/4}$” rates on the regression $Q^\ast$ are also required in luedtke2016statistical and shi2020breaking --- see Remark (ref) below.
We can interpret the asymptotic behavior of our estimator in the context of the results of luedtke2016statistical and van2014targeted. When the maximizing coordinate of $Q^\ast(X) := (Q^\ast(1, X), \dots, Q^\ast(N, X))$ is unique almost surely, van2014targeted remarkably show that performing inference on $V^\ast$ is asymptotically equivalent to being directly given the (almost surely unique) optimal policy $\pi^\ast(x) := \arg\max_\ell Q^\ast(\ell, X)$ and performing policy evaluation. More formally, they show that the efficient influence function for $V^\ast$ is given by \[ \rho^{\mathrm{reg}}_V(Z) := \underbrace{\max_\ell\{Q^\ast(\ell, X)\}}_{= Q^\ast(\pi^\ast(X), X)} + \frac{\mathbbm{1}\{A = \pi^\ast(X)\}}{p^\ast(\pi^\ast(X) \mid X)} \left\{Y - Q^\ast(A, X)\right\} - V^\ast. \] This exactly aligns with our influence function $\rho_V(Z)$ defined in Equation (ref), since in the regular setting $\mathfrak{s}_k(Q^\ast(X)) = 1$ if and only if $k = \arg\max_\ell Q^\ast(\ell, X)$. In particular, this shows our estimator is semi-parametric efficient in this setting.
When the maximizing coordinate of $Q^\ast = (Q^\ast(1, X), \dots, Q^\ast(N, X))$ is non-unique with positive probability, efficiency results in the standard semi-parametric sense are not possible hirano2012impossibility. Further, there is no longer an almost surely unique optimal treatment policy. Any (potentially randomized) treatment policy $\pi$ satisfying $\pi(X) \in \arg\max_\ell Q^\ast(\ell, X)$ almost surely will attain the optimal policy value. Note that for any fixed $u \in \mathbb{R}^N$, the vector $(\mathfrak{s}_1(u), \dots, \mathfrak{s}_N(u))$ defines a probability distribution on $[N]$ that samples an action from the maximizing set $\arg\max_\ell\{u_1, \dots, u_N\}$ uniformly at random. Thus, the influence function $\rho_V(Z)$ corresponds to doing policy evaluation on the randomized policy that selects an optimal action uniformly at random when there are ties.
We do not prove Theorem (ref) in the main body of the paper. We provide our proof in Appendix (ref), which follows as a special case of our result on more general irregular functionals (Section (ref)). At a high level, the proof follows from a careful balancing of the smoothing bias and the second order nuisance estimation errors. In our arguments, we leverage both Lemma (ref) and various analytical properties of the softmax function, which are noted in Lemma (ref). We now provide a very brief proof sketch illustrating how one should select $\beta_n$ to ensure asymptotic normality.
The following corollary allows us to build confidence intervals using the quantiles of the standard normal distribution. Again, we prove said corollary in Appendix (ref). The proof largely boils down to showing the consistency of the plug-in variance estimate and then applying the continuous mapping theorem.
In the previous subsection, we considered performing inference on the value of the optimal treatment policy in a setting where we made no structural assumptions on the outcome regressions/Q-functions. As noted in Remark (ref), there is a price for smoothing: the learner needed to estimate the regressions at slightly faster than the typical $o_\mathbb{P}(n^{-1/4})$ required in causal inference. While existing works such as shi2020breaking and luedtke2016statistical also suffer from this problem, this requirement may be unpalatable when one believes the regressions are particularly hard to estimate.
To allow for more flexible nuisance estimation rates, we now consider a setting where we place semi-parametric restrictions on the underlying outcome regressions, namely through the assumption of linearity of the blip effects in a known feature embedding. Importantly, this is strictly weaker than assuming the Q-function $Q^\ast(a, x)$ is itself linear (which was the assumption made in works such as chakraborty2010inference, goldberg2014comment, laber2014dynamic). This modeling allows the expected baseline response to some fixed treatment to be arbitrarily complicated.
Recalling the identification formula for $V^\ast$ in terms of blip effects from Section (ref), we see that under Assumption (ref) the value of the optimal treatment policy can be written as
If we are able to obtain a $\sqrt{n}$-consistent estimate $\widehat{\theta}$ of the structural parameter $\theta_0$, we can follow the softmax smoothing strategy from the previous section to perform inference on $V^\ast$ via a plug-in estimate. Thus, we start by by identifying $\theta_0$ in terms of a Neyman orthogonal score.
Under standard convergence and regularity assumptions, we can use the score $\Psi$ to perform inference on the structural parameter $\theta_0$. Notably, we will be able to obtain asymptotically linear estimates for $\theta_0$ only assuming the regression $\widehat{\mu}$ converges to to $\mu^\ast$ in $L^2$ at $o_\mathbb{P}(n^{-1/4})$ rates. This is uniformly slower than the rates required on the outcome regression in the previous section. The proof of the following theorem can be found in Appendix (ref)
Using Theorem (ref), we can directly obtain asymptotically linear estimates for the value $V^\ast$ as well. In particular, because the estimates for $\theta_0$ provided by the aforementioned theorem are asymptotically linear, we do not need to construct a Neyman orthogonal score for the value of the optimal treatment policy. We can simply plug-in our estimates and appropriately account for the additional introduced variance. This is detailed in the following corollary.
We prove Corollary (ref) in Appendix (ref) by checking the conditions of Theorem (ref), our smoothed central limit theorem proved in Appendix (ref). Further, we extend our estimator of the value of the optimal treatment policy under semi-parametric assumptions on the blip effect to dynamic settings in Appendix (ref).
Up to this point, we have focused on performing inference on the value of an optimal treatment policy. We now show that softmax smoothing can be applied more broadly to classes of irregular functionals that are specified as the expected point-wise maximum of a collection of sufficiently regular scores. This framework will not only subsume the results in the previous section, but will also be applicable in many other settings as well. Throughout this section, we assume that observations take the generic form $Z$ and that we have $X \subset W \subset Z$ for a generalized set of covariates $X$ and an extended set of features $W$. Often, we will have $W = (X, A)$. Here, by $X \subset W \subset Z$, we mean that $X$ consists of a subset of the coordinates of a vector $W$, and likewise $W$ consists of a sbuset of the coordinates of $Z$. We let $\mathcal{X}$, $\mathcal{W}$, and $\mathcal{Z}$ denote the measurable spaces in which $X$, $W$, and $Z$ take their values. The goal is to perform inference on the maximum of the scores:
Here, for each $k \in [N]$, $\psi_k$ is a score depending on the data vector $X$ and a nuisance component $g_k : \mathcal{W} \rightarrow \mathbb{R}^{d_k}$, and $g_k^\ast(w) \in L^2(P_{W}; \mathbb{R}^{d_k})$ denotes some true, unknown nuisance function. We let $\psi(X; g) := (\psi_1(X; g_1), \dots, \psi_N(X; g_N))$ for convenience throughout, where $g = (g_1, \dots, g_N)$. In the example of the value of an optimal treatment policy, for each $k \in [N]$, we have $W = (A, X)$, $g_k^\ast(a, x) = Q^\ast(a, x)= \mathbb{E}[Y \mid X =x, A = a]$, and $\psi_k(x; g_k) = g_k(k, x)$. We describe how to instantiate this framework with respect to two other examples, conditional Balke and Pearl bounds levis2023covariate and $L^1$ calibration error gupta2022post, in the sequel. Throughout this section, we make the following assumption on the constituent scores $\psi_1(x; g_1), \dots, \psi_N(x; g_N)$.
In words, the first assumption just states that each $g_k^\ast$ is a regression of some random variables $U_k$ onto the extended set of covariates $W$. As noted above, this assumption is trivially satisfied for the value of the optimal treatment policy. The first part of the second assumption says that each score is affine in the nuisance component $g_k$. Given $\psi_k(x; Q) = Q(k, x)$, this is trivially satisfied for the value of the optimal treatment policy. The second part ensures that each constituent score is sufficiently regular to exchange derivatives with expectations. In particular, it will hold whenever
for an integer $P_k$, constants $c_{k, 0}, \dots, c_{k, P_k} \in \mathbb{R}$, indices $i_1, \dots, i_p \in [d_k]$, and values $a_1, \dots, a_{P_k} \in \mathcal{A}$. This is the case for the optimal policy value and all examples in the sequel. Finally, the third assumption states that the scores are sufficiently regular to admit a Riesz representer conditional on $X$. We have $\zeta_k^\ast(a, x) = \frac{\mathbbm{1}\{a = k\}}{p^\ast(k \mid x)}$ for the optimal policy value. We derive the representers for conditional Balke and Pearl bounds and the $L^1$ calibration error in the sequel.
While estimands such as the one outlined in Equation (ref) have been considered by semenova2023aggregated, the existing results have several key limitations. First, all convergence rates of nuisance estimates must be $o(n^{-1/4})$ in the $L^\infty(P_Z)$ norm, which can be a strong requirement when ML learners are used. Second, the argument maximizer is assumed to be unique almost surely, e.g.\ in the case covered in the previous sections the probability of non-response to treatment is assumed to be zero. In this section, we use softmax smoothing to construct confidence intervals (a) only assuming in probability rates of convergence for nuisance estimates and (b) allowing for an arbitrarily large probability of the argument maximizer being non-unique. Once again, the first step in our recipe is to construct a softmax approximation to the main objective. To this end, we define \[ V^\beta := \mathbb{E}\left[\mathbf{sm}^\beta_\ell \psi_\ell(X; g_\ell^\ast)\right]. \] As before, the quantity $\mathbb{E}[\mathbf{sm}^\beta_\ell \psi_\ell(X; g_\ell)]$ is, without modification, highly sensitive to misestimation in the nuisance components $g_1, \dots, g_N$. Thus, we should subtract a first order correction/de-biasing term to arrive at a Neyman orthogonal score. The following proposition accomplishes precisely this --- we provide a proof in Appendix (ref).
Using the above Neyman orthogonal score, we can arrive at a general asymptotic linearity/normality result that allows one to construct confidence intervals for broad classes of irregular statistical parameters. Of particular importance will be the following, limiting score \[ \Psi^\ast(Z; g, \alpha) = \max_\ell\{\psi_\ell(X; g_\ell)\} + \sum_{\ell = 1}^N \alpha_\ell(W)^\top(U_\ell - g_\ell(W)), \] where the limiting $k$th nuisance $\alpha^\ast_k$ is given as the point-wise limit $\alpha^\ast_k(W) := \lim_{\beta \rightarrow \infty}\alpha^\beta_k(W)$. This limit is precisely the quantity \[ \alpha^\ast_k(W) = \mathfrak{s}_k(\psi_\ell(X; g_\ell^\ast))\zeta^\ast_k(W) \] where, as before, $\mathfrak{s}_k(u) := \frac{\mathbbm{1}\{k \in \arg\max_\ell\{u_1, \dots, u_N\}\}}{|\arg\max_\ell\{u_1, \dots, u_N\}|}$. This score is reminiscent of the influence function discussed following Theorem (ref) (up to the centering at $V^\ast$), just generalized to a broader class of irregular parameters. In particular, the same asymptotic interpretation holds for $\Psi^\ast$ --- in the limit of large $\beta$, the estimator breaks ties via uniform randomization. We now state our main result, which we prove in Appendix (ref)
Most of the assumptions above are commonly made in the causal inference literature. The most peculiar of the assumptions is that of score continuity (Assumption 3), as the quantity on the left-hand side only involves the random covariates $X$, whereas the quantity of the right involves the extended set of features $W$. In the context of Equation (ref), this assumption ensures that each value $a_{k, p} \in \mathcal{A}$ in the score $\psi_k$ occurs with conditional probability bounded away from zero given $X$. In the example of the value of the optimal treatment policy, we have $\psi_k(X; Q) = Q(k, X)$, and so mean-squared continuity is implied by the assumption of strong positivity. A similar positivity assumption will ensure mean-squared continuity in the example of conditional Balke and Pearl bounds discussed below. For the other example, that of $L^1$ calibration error, we will have $W = X$, and so the assumption will follow trivially.
In the above, we have chosen to assume the $\psi_k(X; g_k)$ are affine for ease of exposition. In principle, one could generalize our result to non-affine scores that are sufficiently regular. More specifically, one would need to assume the $\ell_{\nu, x}^k(\epsilon)$ as defined above are twice continuously differentiable and that there is another integrable function $G_\nu^{k}(X)$ such that $\sup_{|\epsilon|\leq 1}\|\partial_\epsilon^2 \ell^k_{\nu, X}(\epsilon)\|_2 \leq G_{\nu}^k(X)$. The first assumption allows for twice Gateaux differentiability, and the second envelope enables the exchange of second derivatives with expectations. One would also need to assume conditions analogous to mean-squared continuity for the first and second Gateaux derivatives of $\psi_k(X; g)$. Given that the examples discussed in the paper all fall into the setting of affine scores, we avoid introducing additional, technical assumptions.
As before, the above theorem yields a corollary that allows for the practical construction of asymptotically valid confidence intervals using a plug-in estimate of the variance.
Throughout this work, we have focused on $V^\ast := \mathbb{E}[\max_\ell\{Q^\ast(\ell,X)\}]$ as an example of an irregular parameter. We now discuss another example focused on the partial identification of the average treatment effect (ATE) in instrumental variable settings. We largely follow the exposition given in levis2023covariate, and recommend balke1994counterfactual, balke1995probabilistic, balke1997bounds for further exposition and history. In general, one can consider both upper and lower bounds on the average treatment effect --- we solely discuss lower bounds as results for upper bounds are exactly analogous. The main contributions of this application in comparison to the results of levis2023covariate are that (a) our main theorem allows one to select $(\beta_n)_{n \geq 1}$ to guarantee asymptotic normality around the un-smoothed, population Balke and Pearl bounds (defined below) as opposed to smoothed analogues, and (b) the theorem holds even when the argument maximizer used in defining the Balke and Pearl bounds is non-unique with arbitrarily large probability.
For this example, observations are of the form $Z = (X, D, V, Y)$, $X \in \mathcal{X}$ represent covariates, $D \in \{0, 1\}$ represents a binary treatment, $V \in \{0, 1\}$ a binary instrument, and $Y \in \{0, 1\}$ a binary outcome. For simplicity, we let $W = (X, V)$ denote the tuple consisting of covariates and instrument. We make the following assumptions on the vector of observations $Z$.
Because of the second assumption, it makes sense to define the natural potential outcomes $Y(d) := Y(0, d) = Y(1, d)$. The goal is to provide a lower bound on the average treatment effect $\mathbb{E}[Y(1) - Y(0)]$, which in general is not point identified. We use a the construction of levis2023covariate to fit this example into the framework of general irregular functionals.
For any $d, y \in \{0, 1\}$, define the nuisance function $q_{yd}^\ast(x, v) := \mathbb{P}(Y = y, D = d \mid X = x, V = v)$. Next, define the scores $\psi_1(x; q), \dots, \psi_8(x; q)$ as follows:
In the above, we write $q := (q_{yd} : (y, d) \in \{0, 1\}^2)$ and let $q^\ast$ denote the full-vector of true, unknown regressions. While each score only depends on a subset of the regressions, we let the scores depend on the full vector for simplicity. levis2023covariate show that one has the following lower bound on the ATE: \[ V^\ast := \mathbb{E}[\max_{\ell = 1}^8 \psi_\ell(X; q^\ast)] \leq \mathbb{E}[Y(1) - Y(0)]. \] Our goal is to show that our smoothing approach can be used to perform inference on the lower Balke and Pearl bound $V^\ast$. In the following corollary, for each $\ell \in [8]$, we leverage the representer $\zeta^\ast_\ell(V,X) \in \mathbb{R}^4$ whose coordinates are defined as
where $\mathrm{sgn}_\ell(y, d, v) = 1$ if $q_{yd}(x, v)$ appears in the definition of $\psi_\ell$ with positive sign, $\mathrm{sgn}_\ell(y, d, v) = -1$ if it appears with negative sign, and $\mathrm{sgn}_\ell(y, d, v) =0$ otherwise.
We prove Corollary (ref) by first showing that, under Assumption (ref), the constituent scores satisfy Assumption (ref). We then check the relevant conditions of Theorem (ref). We provide a proof in Appendix (ref).
We now consider an application of Theorem (ref) related to the calibration literature. Calibration is a notion of consistency for an estimator $\theta$ that establishes the trustworthiness of its predictions kumar2019verified, gupta2022post. Formally, perfect calibration is defined as follows.
We use $O$ in this section to represent observed covariates to prevent notational collision with $X$ as used in Theorem (ref) and Assumption (ref). In words, $\theta$ is perfectly calibrated if, on average, whenever it makes a prediction $\theta(O) = x$, the true outcome is also $x$. Given most estimators will not be perfectly calibrated in practice, various notions of approximate calibration have been developed in the calibration literature. Perhaps the most commonly used approximation is the $L^1$ calibration error, which aims to capture the “distance” from perfect calibration by measuring the absolute deviations from perfect calibration on each level set of $\theta$ naeini2015obtaining, gupta2022post. Other analogous notions, such as $L^2$ calibration error kumar2019verified, whitehouse2024orthogonal, lee2023t, van2023causal, have been proposed in the literature as well, but we do not define those in this paper. The $L^1$ calibration error is defined as follows.
The goal of this subsection is to use Theorem (ref) to construct asymptotically-valid confidence intervals for the $L^1$ calibration error of a fixed (i.e.\ non-random) model $\theta : \mathcal{O} \rightarrow \mathbb{R}$. Estimating the $L^1$ calibration error is understood to be a hard problem in the calibration literature, and gupta2022post note that it is impossible to construct unbiased estimates in general. It is also non-trivial to construct asymptotically-normal estimate for $\mathrm{Cal}_1(\theta)$. First, the calibration function $\chi^\ast$, which is the projection of $Y$ onto the space of $\theta(O)$-measurable functions, is an unknown nuisance component and thus must be estimated from data. Second, the $L^1$ calibration error is defined as the expected absolute deviation between $\theta$ and $\chi^\ast$, and is thus not differentiable when both are equal to zero. This could happen when the estimator is perfectly calibrated, or when it is calibrated on certain subsegments of the population. In sum, performing inference on $\mathrm{Cal}_1(\theta)$ is a semi-parametric problem that falls into the scope of Theorem (ref). To prove our result, we just make the following assumptions.
With the above assumption, we can now present the main result of this subsection. For this example, we take $X = W := \theta(O)$, and take the outcome to be simply $Y$. The constituent scores are $\psi_1(x; \chi), \psi_2(x; \chi)$, defined respectively as \[ \psi_1(x; \chi) := \chi(x) - x \quad \text{and} \quad \psi_2(x; \chi) := x - \chi(x). \] One can check the identity \[ \mathrm{Cal}_1(\theta) := \mathbb{E}|\theta(O) - \chi^\ast(\theta(O))| = \mathbb{E}\big[\max\{\psi_1(\theta(O); \chi^\ast), \psi_2(\theta(O); \chi^\ast)\}\big]. \] The following corollary follows in a straightforward manner from Theorem (ref). We prove the following result in Appendix (ref).
We now evaluate the coverage of our softmax smoothing methods from Section (ref) in both non-parametric and semi-parametric settings. We also also measure the sensitivity of coverage to the choice of smoothing parameter $\beta$. Analyzing this sensitivity is important for practical applications of our estimator, as the exponent $\delta$ that governs the density of the sub-optimality gaps may not be known to the user.
Throughout this section, we study a setting with a binary treatment for simplicity. To align with the majority of the causal inference literature, we let our action set in this case be denoted by $\{0, 1\}$ and let $\tau^\ast(x) := \mathbb{E}[Y(1) - Y(0) \mid X = x] = Q^\ast(1, x) - Q^\ast(0, x)$ denote the CATE. We start by describing data generating processes in both the non-parametric and semi-parametric settings which ensure Assumption (ref) holds.
\paragraph{Non-Parametric Data Generating Process}
We describe how a single sample $Z := (X, A, Y)$ is generated in the non-parametric setting. We first define the abstract data-generating process, and then specify the parameter values used in our experiments. Let
where all random variables are mutually independent. Fix parameters $\pi_0 \in (0,1)$, $\rho \in (0,1)$, $\delta > 0$, and $t_0 > 0$. In our applications, we take $d = 3, \pi_0 = 0.2, \rho = 0.3$, and $t_0 = 0.3$. Further, define \[ U := \frac{C_1 - \pi_0}{1-\pi_0}. \] We let the baseline regression under control be given by \[ Q^\ast(0,X) := 1.0 + 0.5W_1, \] and define the CATE $\tau^\ast(X)$ by
where $h:\mathbb{R} \to \mathbb{R}_{\geq 0}$ is a non-negative function, which we take to be $h(W_1) := 0.5|W_1|$. We then define \[ Q^\ast(1,X) := Q^\ast(0,X) + \tau^\ast(X). \] Treatments are sampled in the following manner: \[ A \mid X \sim \mathrm{Bern}(p^\ast(X)),\;\; \text{where } p^\ast(X) = \frac{1}{1 + \exp\{-0.75(-0.2 + W_1)\}}, \] and outcomes are generated as \[ Y := Q^\ast(0,X) + A\cdot \tau^\ast(X) + \sigma\epsilon, \qquad \epsilon \sim \mathcal{N}(0,1), \] where $\epsilon$ is independent of $(X,A)$ and we take $\sigma = 0.3$ in our settings.
First, because $\tau^\ast(X) = 0$ whenever $C_1 \leq \pi_0$, we have $\mathbb{P}(\tau^\ast(X)=0)=\pi_0>0$, and thus there is a positive probability of non-response to treatment. Second, one can check that the random variable $|\tau^\ast(X)|$ admits a Lebesgue density $f_\tau(u)$ on the interval $(0, t_0)$ that satisfies $f_\tau(u) \lesssim u^{\delta - 1}$. That is, the data-generating process satisfies the density condition from Assumption (ref) with parameter $\delta$.
\paragraph{Semi-Parametric Setting} For the semi-parametric setting, a single sample $Z := (X, A, Y)$ is generated as follows.
where we set the structural parameter as $\theta_0 = 1.0$. In this setting, there is zero-probability of non-response, and because the blip effect using $A = 0$ as a reference is precisely the CATE, we see that $\tau^\ast(X) = \theta_0(X_2 + 1.0)$ and thus Assumption (ref) is satisfied with $\delta = 1.0$ (corresponding to a bounded density).
\paragraph{Evaluation of Non-Parametric Estimator} We start by describing how we measure coverage in a single run of our method. We then discuss how we aggregate our experimental results over independent runs to obtain estimates of coverage with corresponding confidence intervals. We generate $n = 5,000$ independent samples from the data generating process, where we let $\delta^{\mathrm{true}} \equiv \delta \in \{0.75, 1.0, 1.25\}$. We then compute the estimator outlined in Theorem (ref) as follows. We use $K$-fold cross-fitting (see Remark (ref)) with $K = 5$ and stratified splitting based on the treatment $A$. Using the data outside of each fold $k \in \{1, 2, \dots, 5\}$, we separately train nuisance models $\widehat{Q}^{(-k)}(0, X), \widehat{Q}^{(-k)}(1, X)$, and $\widehat{p}^{(-k)}(X)$ predicting $Q^\ast(0, X), Q^\ast(1, X),$ and $p^\ast(X)$ respectively. In each run, we estimate the Q-functions $Q^\ast(0, X)$ and $Q^\ast(1, X)$ using the LGBMRegressor class from the module LightGBM ke2017lightgbm and we estimate the propensity $p^\ast(X)$ using the LogisticRegression class from the Scikit-learn module pedregosa2011scikit with $\ell_2$ regularization.\footnote{We use the following parameter settings for LGBMRegressor \{`metric' : `rmse', `learning_rate': 0.05, `num_leaves' : 2**3, `n_estimators' : 250, `max_depth' : 3, `verbose' : -1, `random_state' : 123\} . For LogisticRegression, we use the parameter settings `random_state' : 123 and otherwise use default settings.} For each $\delta^{\mathrm{true}}$, we consider a grid of smoothing parameters to measure sensitivity of the coverage of the smoothing estimator. For each assumed parameter \[ \delta^{\mathrm{ass}} \in \{\delta^{\mathrm{true}} - 0.5, \delta^{\mathrm{true}} - 0.25, \delta^{\mathrm{true}} - 0.1, \delta^{\mathrm{true}}, \delta^{\mathrm{true}} + 0.1, \delta^{\mathrm{true}} + 0.25, \delta^{\mathrm{true}} + 0.5\} \] we compute the softmax-smoothing estimator with
We produce corresponding confidence intervals using Corollary (ref). We multiply by $\log\log(n)$ to ensure $\beta_n = \omega\left(n^{\frac{1}{2(1 + \delta^{\mathrm{ass}})}}\right)$, which is needed when $\delta^{\mathrm{ass}} = \delta^{\mathrm{true}}$ per Theorem (ref) to maintain coverage. For each setting of parameters $(n, \delta^{\mathrm{true}}, \delta^{\mathrm{ass}})$ (which implicitly defines $\beta_n$) outlined above, we repeat the measurement process $M = 250$ times and report 95% Wilson Binomial confidence intervals (CIs).
\paragraph{Evaluation of Semi-Parametric Estimator} In the semi-parametric setting, we estimate the structural parameter $\theta_0 = 1.0$ using the moment condition estimator outlined in Theorem (ref), again using $K$-fold cross-fitting with stratified splitting and $K = 5$. All nuisances are estimated using the LGBMRegressor class with the same hyperparameter settings as above. We then use the estimate of the structural parameter to compute an estimate of the optimal policy value using the estimator in Corollary (ref). In this setting, we consider the grid $\delta^{\mathrm{ass}} \in \{0.5, 0.75, 0.9, 1.0, 1.1, 1.25, 1.5\}$ and compute an estimator for each corresponding value of $\beta_n$ as noted in Equation (ref). Again, for each setting of $\delta^{\mathrm{ass}}$, we repeat the measurement process $M = 250$ times and produce Wilson CIs as above.
\paragraph{Non-Parametric Estimator} We display a summary of the coverage of our estimator in the first three columns of Table (ref) and more extensive measurements of sensitivity in Figure (ref). Table (ref), which assume $\delta^{\mathrm{ass}} \in \{\delta^{\mathrm{true}} - 0.1, \delta^{\mathrm{true}}, \delta^{\mathrm{true}} + 0.1\}$, shows our estimator obtains close to 95% coverage across settings of $\delta^{\mathrm{true}}$ even when $\delta^{\mathrm{ass}}$ is incorrect. Figure (ref), which displays a wider range of $\delta^{\mathrm{ass}}$, seems to indicate that having $\delta^{\mathrm{ass}} \gg \delta^{\mathrm{true}}$ can lead to incorrect coverage. This makes sense in the context of our theoretical results, as when $\delta^{\mathrm{ass}} > \delta^{\mathrm{true}}$ (over-smoothing), the selected $\beta_n$ will satisfy $\beta_n = o\left(n^{\frac{1}{2(1 + \delta^{\mathrm{true}})}}\right)$, and thus our asymptotic normality results will fail to hold. We note that in settings where $\delta^{\mathrm{ass}} \leq \delta^{\mathrm{true}}$ (under-smoothing) we generally observe valid coverage as well. Overall, our experimental results seem to indicate that it is advantageous to under-smooth (i.e.\ to select $\beta_n \gg n^{\frac{1}{2(1 + \delta)}}$) if there is uncertainty about the behavior of the density of the CATE near zero. This indicates that the bias from smoothing is likely the limiting factor in maintaining valid coverage. While under-smoothing requires the $L^2$ error of the regression estimate to decay as $o_\mathbb{P}(\beta_n^{-1/2}n^{-1/4})$, flexible learners may readily be able to obtain such rates for reasonable DGPs.
\paragraph{Semi-Parametric Estimator} In the semi-parametric setting, due to the assumed linearity of the blip effects, there is no real trade-off between first-order biases and second order nuisance errors. That is, the learner just needs to select $\beta_n = o(n^{1/2})$ and $\beta_n = \omega\left(n^{\frac{1}{2(1 + \delta^{\mathrm{true}})}}\right)$ in order to obtain valid asymptotic coverage. Our coverage findings, illustrated in the final column of Table (ref) and in the final image of Figure (ref), indicate a general lack of sensitivity to the choice of the assumed density parameter $\delta^{\mathrm{ass}}$. Interestingly, even under-smoothing does not seem to result in invalid coverage.
In this paper, we considered the problem of constructing confidence intervals for the value of the optimal treatment policy. This is known to be a difficult problem in semi-parametric inference, as the non-Gateaux differentiability of the target functional prevents the application of standard, de-biased estimators luedtke2016statistical, shi2020breaking, laber2014dynamic, chakraborty2013inference, semenova2023aggregated. Our approach is intuitive and computationally lightweight, replacing the non-differentiable maximum function with a smoothed analogue defined in terms of the softmax function. Through a careful analysis of this smoother, we show one can construct estimators for the optimal policy value that allow for an arbitrarily large probability of non-response to a treatment. This estimator naturally generalizes to other important examples of irregular parameters, such as conditional Balke and Pearl bounds on the average treatment effect levis2023covariate and the $L^1$ calibration error of an ML model gupta2022post. In sum, our work shows that the careful choice and analysis of a smoothing function can turn irregular problems into ones that can again be handled by first-order de-biasing methods.
While the results presented in this paper are quite general, many important open questions remain. First, throughout our work we assume that the sequence $(\beta_n)_{n \geq 1}$ is fixed and chosen according to the behavior of the sub-optimality gaps in a small neighborhood of zero. Although our experiments demonstrated some lack of sensitivity to the choice of smoothing parameter, it would be of great theoretical and practical interest to develop data-adaptive methods for selecting $\beta_n$, if this is possible. Second, while we established efficiency in the setting where the optimal treatment is almost surely unique, it would be interesting to reason about efficiency when non-response occurs with positive probability. In this setting, hirano2012impossibility show that regular asymptotically linear estimators do not generally exist for such target parameters. Since semi-parametric efficiency is most often described in terms of such estimators, one would need to develop an analogue of efficiency theory for this irregular setting. We view this as an important direction as it would facilitate the direct comparison of estimators in irregular regimes.