Statistical Inference of Optimal Allocations I: Regularities and their Implications

Statistical Inference of Optimal Allocations 1: Regularities and their Implications

Introduction

Characterizing individual heterogeneity is crucial in the optimal allocation of scarce resources especially as a consequence of the rapid digitization of social-economic data. Examples include medical treatment assignments when doctors adopt different therapeutic approaches to different patients, pretrial bails when judges decide where defendants will await trials, college admissions when committees assign students to various educational programs. In order to address the heterogeneity inherent in such problems, where the potential effects of a treatment on an outcome vary among individuals, personalized allocation is more of a norm than an exception. Personalized allocation is compelling in the presence of the inherent individual heterogeneity in the potential effects of a treatment on an outcome.

Based on past experimental or observational data, the core objective of optimal allocation problem is to find a strategy that can perform well out of sample. An effective approach advocated by the literature is to make use of information provided by prediction algorithms. See for example QJEbail, agrawal2019artificial and babina2024artificial. A common decision framework is to form a nonrandom treatment rule using an indicator function that benchmarks a first step regression estimate of individual level treatment effect against certain thresholds, representing a plug-in strategy in a weighted classification problem.

Such a plug-in strategy is based on solving a population optimal allocation problem. For a joint distribution of $\(X, Y\right)$ where $Y = \(Y_{j}, j=0,\ldots,J\right)$, consider the social welfare potential function\footnote{Appendix (ref) provides a mathematical justification for the term “potential function”.}

In the above, $\Phi$ is the set of multi-critical functions, i.e. $\phi = \left(\phi_{j}, j=0,\ldots,J\right)$, $\phi_j: \Omega\rightarrow \mathbb{R}$, with $\left(\Omega, F, \mu\right)$ being a Borel probability space, such that $\phi_j \in L^\infty\left(\mu\right), 0 \leq \phi_j\leq 1, \forall j\in \{0,\ldots,J\}$, and such that $\sum_{j=0}^J \phi_j\(x\right) \equiv 1$ for all $x \in \Omega$; $g\left(\cdot\right) = \(g_j\left(\cdot\right), j=0,\ldots,J\right)$ and $g_j\(x\right) = \mathbb{E}\(Y_j \vert X = x\right)$. Under very general conditions, the problem posed in (ref) admits an intuitive solution where the $j$th-category receives all allocation when $\mathbb{E}\left(\lambda_j Y_j \vert X\right)$ is the unique maximum. If there exists some $l$ such that $\mathbb{E}\left(\lambda_j Y_j \vert X\right) = \mathbb{E}\left(\lambda_l Y_l \vert X\right)$, then all allocation can be divided among the maximal indexes.

In this paper, we are interested in differentiating $\gamma\left(\lambda, g\right)$ in (ref) functionally with respect to both $\lambda$ and $g\left(\cdot\right)$. The solution can provide inference methods for the value of optimal treatment allocation in multiple applications. In section (ref), we derive the asymptotic distributions and justify bootstrapping for binary optimal allocation subject to a resource constraint and for the receiver operating characteristic (ROC) curve in binary classification. Section (ref) provides asymptotic regret bounds and the asymptotic distribution of the value function for unconstrained multi-class optimal allocation.

Our paper attempts to synthesize several strands of literature. It is widely-known in the plug-in classification literature that the error of the second step classification is smaller than that of the first step regression devroye1996probabilistic, audibert2007fast. Our paper builds on seminal works by chernozhukov2018sorted and chen2003estimation to provide conditions supporting $\sqrt{n}$ value function convergence rates and tight regret bounds in multiclass and constrained optimal allocations, covering settings both for classical empirical processes and for recent advances in Neyman orthogonality and sample splitting chernozhukov2018double, chernozhukov2022locally. Subsequent works following manski2004statistical have defined regret as the difference between the social welfare potential and the expected utility achieved by an estimated decision rule.

In the binary case of $J=1$, $\lambda = \left(\lambda_0, \lambda_1\right)$. Consider the finite dimensional derivatives of $\gamma\left(\lambda, g\right) =\gamma\left(\lambda_0, \lambda_1, g\right)$ with respect to $\lambda_0$. Assuming $\mathbb{P}\{\lambda_1 g_1\(X\right)=\lambda_0 g_0\(X\right)\}=0$, we can write

Given the definition of $\gamma\left(\lambda_0,\lambda_1, g\right)$, we expect based on intuitive reasoning that

In a rigorous proof of (ref) using direct calculation, the terms of major difficulties are

We will see later in section (ref) that under mild conditions, the two terms in (ref) cancel out because they converge to the same integral on the level set $\{x: \lambda_0 g_0\(x\right)=\lambda_1 g_1\(x\right)\}$ with opposing signs. As a consequence, the desired partial derivatives exist and take the displayed form (ref).

Our approach in section (ref) is greatly inspired by the integration on manifold approach in chernozhukov2018sorted. In particular, we reapproach and generalize their results based on the concept of Hausdorff measure and powerful area and coarea integration formulas from geometric measure theory. We show that the Hadamard derivative of $\gamma\left(\lambda, g\right)$ contains terms that are integrals under certain Hausdorff measure. The Hadamard differentiability results are derived under mild primitive conditions instead of being directly assumed. See for example sherman1993limiting and luedtke2016optimal.

In the general multi-class setting, the social welfare potential function in (ref) is not only subadditive but also positive homogeneous of degree one with respect to both $\lambda$ and $g\left(\cdot\right)$, respectively. Consequently, while $\gamma\left(\lambda, g\right)$ is not necessarily jointly convex in $\left(\lambda, g\right)$, $\gamma\left(\lambda, g\left(\cdot\right)\right)$ is convex in $g\left(\cdot\right)$ given $\lambda$ and is convex in $\lambda$ given $g\left(\cdot\right)$. Interestingly, this simple observation seems to have gone unnoticed by the existing literature. As discussed in section (ref), profound results from geometric functional analysis imply that under nearly no assumptions, the social welfare potential function $\gamma\left(\lambda, g\right)$ is Fr\'{e}chet differentiable on a large set and admits an application of the envelope theorem.

The implied envelope theorem allows us to verify the degeneracy of the first order Hadamard derivative of a multi-discrete social welfare potential function. In section (ref), we show that the primitive conditions in section (ref) can also be used to generalize fast convergence rates for the regret of binary plug-in classifiers audibert2007fast to multi-discrete allocation problems. The same methodology can be combined with recent developments of double/debiased machine learning chernozhukov2018double,chernozhukov2022locally to speed up the estimation of $\gamma\left(\lambda, g\right)$.

Literature

A multidisciplinary literature has studied the construction of optimal decision rules from experimental and observational data in the forms of treatment allocation, policy learning and individual treatment rule (ITR). A large body of the literature focuses on empirical risk minimization (ERM). See for example manski2004statistical, qian2011performance, zhao2012estimating, swaminathan2015batch, zhou2017residual, kitagawa2018should, rai2018statistical, luedtke2020performance, athey2021policy, mbakop2021model, zhou2023offline, ben2024policy and ai2024data, among others. An important objective is to provide probabilistic bounds for

where $\phi^{*}$ solves the population optimization program (ref). The function $F\(X, D, Y\right)$ usually takes the inverse propensity weighting (IPW) or doubly robust (DR) form and $\Phi_{0} \subset \Phi$ is not too complex in an entropy sense. See for example kitagawa2018should and athey2021policy.

A statistical difficulty in the optimal allocation problem arises from the indicator function in the population solution qian2011performance. To address such a problem, we build on a divergent literature motivated by very different purposes. To restore monotonicity in conditional quantile estimation, the functional derivative of the sorting operator in the univariate case is studied by chernozhukov2010quantile and subsequently extended to the multivariate cases by chernozhukov2018sorted using calculus on manifold techniques. kim1990cube hinted at a rudimentary form of calculus on manifold in deriving the large sample properties of cube root consistent estimators. sasaki2015quantile incorporated Hausdorff measure tools in fluid mechanics to characterize the information content of quantile partial derivatives for general structural functions.

Optimal allocation problems appear in different but related settings. Examples include asymptotically minimax optimal decision under the limits of experiments framework hirano2009asymptotics, optimal decision under minimax regret stoye2009minimax, tetenov2012statistical, ben2024policy, uniform confidence interval armstrong2023inference, allocation with spillover effects kitagawa2023individualized, kitagawa2023should, viviano2025policy, as well as binary constrained optimal allocation bhattacharya2012inferring, luedtke2016optimal for which our general functional Hadamard derivative results in section (ref) are applicable.

Closely related is the receiver operating characteristic (ROC) curve widely used in binary classification. In the economics and finance literature, the plug-in ROC curves are used for example in berge2011evaluating, vallee2019marketplace, sadhwani2021deep, feng2025statisticalMS to evaluate prediction and decision making performance. A previous literature studying the inference of the ROC curves, including hsieh1996nonparametric, lloyd1998using, li1999semiparametric, hall2004nonparametric and bertail2008bootstrapping, did not account for the sampling error in the estimated propensity score. An exception is luckett2021receiver which derives uniform asymptotics for the estimated ROC curve using a support vector machine (SVM) classifier. Unlike previous studies, we derive the asymptotic properties of the ROC curves obtained by plugging-in the propensity score estimators, which is consistent with the prevalent empirical usage of predictive classifiers.

The methodology in this paper also has an intriguing connection with the faster convergence phenomenon in the classification literature devroye1996probabilistic, audibert2007fast. Consider, especially, the margin assumption (MA) evoked by works including mammen1999smooth, tsybakov2004optimal, boucheron2005theory, audibert2007fast, kitagawa2018should, luedtke2020performance and semenova2023debiased:

for fixed $c$ and some constant $C > 0$ and $\alpha \geq 0$. Our conditions in section (ref) directly imply the margin assumption with $\alpha \geq 1$. This connection can be further combined with the recent double/debiased machine learning approach chernozhukov2018double, chernozhukov2022locally to form a transparent cross-validation and plug-in based method to complement the optimal allocation literature.

Hadamard differentiability

Following the seminal work by chernozhukov2018sorted, we consider (ref) as a functional mapping and are interested in the infinitesimal change in $\gamma\left(\lambda, g\left(\cdot\right)\right)$ induced by a marginal change in $\left(\lambda, g\left(\cdot\right)\right)$. More precisely, we derive the Hadamard derivative of $\gamma\left(\lambda, g\left(\cdot\right)\right)$ with respect to $\left(\lambda, g\left(\cdot\right)\right)$, which justifies and generalizes the calculations in (ref). We assume that the underlying distribution of $X$, denoted as $\mu$, is absolutely continuous with density function $f\(x\right)$ with respect to $\mathcal{L}_{n}$, the $n$-dimensional Lebesgue measure. Then we write

Typically, the optimizing policy function takes the form of $\phi_j\(x\right) = \prod_{l \neq j} 1\left(\lambda_j g_j\(x\right) > \lambda_l g_l\(x\right) \right)$. For simplicity we abbreviate $\lambda_j g_j\(x\right)$ as $g_j\(x\right)$. The main object of study is

which can be generalized to a form of the following sorting operator (as defined by chernozhukov2018sorted) for vector-valued function $h\(x\right)$ and vector-valued cutoff threshold $c$,

The number of terms in the product should be less than $n$, the dimension of $X$.

Hausdorff measure and integration formulas

The following discussions are based on standard references by federer2014geometric and evans2018measure. In order to extend the analysis in chernozhukov2018sorted, this section first defines extrinsically the notion of length, area and volume. For this purpose, the Carath\'{e}odory criterion is used to construct a measure space from a $\sigma$-subadditive set function (outer measure) on power set of a space. More specifically, we introduce the series of Hausdorff's outer measures (See definition (ref) for outer measure) which does not require any local parameterization or regularity assumption on these sets.

In $\mathbb{R}^{n}$, $\mathcal{H}^{*}_{n} = \mathcal{L}^{*}_{n}$, $\mathcal{H}_{n} = \mathcal{L}_{n}$, where $\mathcal{L}^{*}_n$ denotes the Lebesgue outer measure and $\mathcal{L}_{n}$ denotes the Lebesgue measure on $\mathbb{R}^{n}$. $\mathcal{H}^{*}_{0}\(E\right)$ measures the cardinality of the set $E$:

We focus on the Hausdorff measures $\mathcal{H}_{k}$ defined on Euclidean spaces $\mathbb{R}^{n}$ for $k \leq n$.

The following Theorem (ref) generalizes the conventional integration change of variable formula for invertible mapping to allow for non-invertible mapping and an increase in dimension.

Furthermore, the next curving Fubini-Tonelli theorem allows for slicing among curving surfaces and mapping from a higher dimension to a lower dimension.

Hadamard derivatives of the sorting operator

In this section, we present our mathematics results on the Hadamard differentiability of the sorting operator in (ref), where we intend to compute

Because the main term that needs to be accounted for is

we focus on the following simplified form of the sorting operator with vector-valued $c$ and $h\(x\right)$:

Without loss of generality, the $g_j\(x\right) f\(x\right)$ term in (ref) is replaced by $f\(x\right)$ in (ref) above.

First, in the special case of the sorting operator where $h\(x\right)$ is a scalar function, consider the derivative with respect to $c$.

Let $h$ be $C^{1}$ or Lipschitz, with $\Vert \nabla h\(x\right)\Vert > 0, \mathcal{L}_{n} \; a.e. \; x.$ Suppose also that $f$ is an integrable function. By the coarea formula Theorem (ref),

Further, by the Lebesgue differentiation theorem, the conventional derivative is given by

for $\mathcal{L}_{1} \; a.e. \; c$. The right hand side above is also a derivative in the distributional (weak) sense.

To verify Hadamard differentiability, we calculate the derivative of the sorting operator,

where $\mathcal \mathcal{H}_{n} \to H$ and $t_{n} \downarrow 0$ as $n \to \infty$. This limit is closely related to the Hausdorff measure $\mathcal{H}_{k}$ along the {\it surface} $\{x: h\(x\right) = c\}$, which defines the {\it volume} of $k$-dimensional surfaces in $n$-dimensional spaces, $k \in \{1, \ldots, n\}$. The important insight of chernozhukov2018sorted is a representation by an integral under $\mathcal{H}_{n - 1}$.

Let $E \subset \mathbb{R}^{n}$ be an open set, $k \in \{1, 2, \ldots, n\}$, $h: E \mapsto \mathbb{R}^{k}$ be a $C^{1}$ function, $f: E \mapsto \mathbb{R}$ be a continuous function. Let $C\(E, \mathbb{R}^{k}\right)$ denote the space of $\mathbb{R}^{k}$-valued continuous functions on $E$ equipped with the sup-norm. We need the following technical assumptions about $h$ and $f$.

Assumption (ref) guarantees that $K_f \cap \{x: h\(x\right) = c'\}$ is close to $K_f \cap \{x: h\(x\right) = c\}$ when $c'$ is close to $c$. The purpose of the Assumption (ref) is to guarantee that $h$ is regular in a neighborhood of $K_f$ in order to apply the implicit function theorem.

Embedded in the proof of Theorem (ref) is the continuity of the derivative of the sorting operator, expressed in the following Proposition (ref), which handles the general case when $k \geq 1$.

We next consider a variant of Theorem (ref) under Lipschitz conditions.

The next theorem (ref) is a nonsmooth version of Theorem (ref).

The assumptions employed in Theorem (ref) play several roles. Assumption (ref) provides compactness and Assumption (ref) guarantees $h$ is regular. The compactness Assumption (ref) ensures that the set $\{x: h\(x\right)=c\} \cap K_f$ in the area of integration is not too large and changes continuously with $c$, meaning that each open covering of $\{x: h\(x\right)=c\} \cap K_f$ covers $\{x: h\(x\right)=c'\} \cap K_f$ for all $c'$ close enough to $c$. Assumption (ref) requires each $x \in E$ to be a regular point of $h\left(\cdot\right)$. Under this assumption, $h\left(\cdot\right)$ can be expressed by the implicit function theorem to facilitate the change of variables.

In addition, Assumption (ref) requires each level set $\{x: h\(x\right)=c\}$ to be sufficiently regular in the sense of having a small number of nonsmooth points and the integrand to be almost everywhere continuous. Under this assumption, derivatives can be taken almost everywhere on the level set. Assumption (ref) can be replaced by continuous differentiability of $h$ and continuity of $f$.

Alternative assumptions are possible. Consider Theorem (ref). Therein, Assumption (ref), which requires compactness of the support of $f\left(\cdot\right)$, can be replaced by the next Assumption (ref) which requires that local union of level sets of $h\left(\cdot\right)$ is compact.

In addition, when $k = 1$ and when we only need Hadamard differentiability on a singleton $\{c\}$, where $c \in h\(E\right)$, Assumption (ref) can be replaced by the following

Binary treatment allocation and the ROC curve

In this section, we derive uniform asymptotic distributions for a binary treatment allocation value function under a resource constraint. Using Hadamard differentiability results from section (ref), we directly and concisely analyze a plug-in two-step estimator of the value function. The functional delta methods alone (see for example van2023weak and fang2019inference) are sufficient to obtain asymptotic results. An important result of this section is showing consistency of a computationally feasible bootstrap procedure for the plug-in two-step ROC estimator commonly used in computer science and other fields.

Binary constrained allocation

In population, binary optimal constrained treatment allocation takes the form of

where $\Phi$ is the set of critical functions with $J = 1$. By the Neyman-Pearson lemma (see for example Theorem 3.2.1 in lehmann2022testing), solutions of (ref) satisfy

where $k$ is a constant to be determined, $g_{j}\(x\right) = \mathbb{E}\(Y_{j}\vert X = x\right)$, $c_{j}\(x\right) = \mathbb{E}\(Z_{j}\vert X = x\right)$, $j=0,1$.

The constrained population programming in (ref) uses observed covariates to allocate a binary treatment among a population when a budget constraint limits the size of the treated group. In this context, $Y_{1}$ and $Y_{0}$ are the welfare measures with and without treatment; $Z_{1}$ and $Z_{0}$ are the respective costs of treatment and control. For a typical economic problem, the cost of treatment is usually no less than the cost of control. For example, bhattacharya2012inferring study the efficient provision of anti-malaria bed net subsidies using a set of features $X$ that include the presence of a child under 10, the wealth per capita and the ownership of a bank account. Here, $Y_{1}$ and $Y_{0}$ are indexed by whether a household is covered by an insecticide-treated net; $Z_{1}$ is the cost of a redeemed net coupon and $Z_{0}$ is taken to be zero. In elliott2013predicting, $\(Y_1, Y_0, Z_1, Z_0\right)$ represent the utility of a true positive, true negative, false positive and false negative, respectively, in a classification problem. In luedtke2016optimal, the resource constraint takes the form of an upper bound on the fraction of the treated population. Their setting is equivalent to $Z_1=1$ and $Z_0=0$. In this paper, we therefore assume $Z_{1} \geq Z_{0}$. This section first considers full observability of $\(Y_1, Y_0, Z_1, Z_0, X\right)$.

Following empirical process notation conventions, let $Qf = \int fdQ$. We then express the value and constrain functions of (ref) as

where level set $\left\{x: \(g_{1}\(x\right) - g_{0}\(x\right)\right) = k\(c_{1}\(x\right) - c_{0}\(x\right)\right)\right\}$ is of zero probability mass. The Hadamard differentiability results in section (ref) allow us to derive the asymptotic distribution of an estimator for (ref) by the functional delta method under the following assumptions. We largely adopt the styles in chernozhukov2018sorted.

Theorem (ref) simplifies under correct specification when $\Delta^{\star}\left(\cdot\right) = \Delta\left(\cdot\right)$. In this case, $k = \frac{f_{\beta}\(Q, \Delta^{\star}, k\right)}{f_{\alpha}\(Q, \Delta^{\star},k\right)}$, since on $\Delta^{\star}\(x; k\right) = 0$, $\Delta\(x; k\right) = \(g_{1}\(x\right) - g_{0}\(x\right) - k \(c_{1}\(x\right) - c_{0}\(x\right)\right) \right) =0$,

Then the asymptotic distribution in (ref) simplifies to a Gaussian process in $\alpha$:

In parametric models where $\Delta\(x\right) = \Delta\(x; \theta\right)$, $r_n = \sqrt{n}$ and $t_Q = 1$. The covariance kernel for (ref) over $\left(\alpha, \alpha'\right)$ can be consistently estimated by its sample analog. Appendix (ref) analyzes the implications of parametric model misspecification.

An alternative model-based representation of (ref) is

Corresponding to (ref), the measure $Q$ for $\(Y_{1}, Y_{0}, Z_{1}, Z_{0}, X\right)$ in the functional expression of (ref) can be changed to the measure $\mu$ for $X$:

When the conditional data generating process of $y_1, y_0, z_1, z_0$ given $x$ is a parametric model $f\(y_1, y_0, z_1, z_0 \vert x; \theta\right)$, we denote $\hat \Delta\(x\right) = \Delta\(x; \hat\theta\right)$. A plug-in estimator of (ref) is then

We will call $\beta\left(\mathbb{Q}_{n}, \hat{\Delta}, \alpha\right)$ in (ref) an orthogonal estimator since its first order derivative cancels out under correct specification, and $\beta\left(\hat\mu, \hat\theta, \alpha\right)$ in (ref) a parametric model-based estimator.

Typically, $\hat\theta$ is the maximum likelihood estimator: $\hat\theta = \mathop{\arg\max}\limits_{\theta\in\Theta} \sum_{i=1}^n \log f\(y_{i1}, y_{i0}, z_{i1}, z_{i0} \vert x; \theta\right).$ Let $\omega = \{y_1, y_0, z_1, z_0\}$ and denote the score function $s_{\theta}\left(\omega \vert x; \theta\right) = \frac{\partial}{\partial \theta} \log f\left(\omega \vert x, \theta\right)$. It is a standard result from maximum likelihood estimation that

See for example newey_mcfadden section 3 or van2000asymptotic chapter 5 for sufficient regularity conditions for the asymptotic linear representation.

Under correct specification, the parametric model-based estimator has an efficiency advantage over the orthogonal estimator. We have the following efficiency comparison proposition.

The population value function $\beta\left(\mu, \Delta, \alpha\right)$ of the constrained optimization problem in (ref) satisfies concavity and monotonicity shape restrictions in the arguments as in the general framework developed by chen2021shape. The orthogonal estimator $\beta\left(\mathbb{Q}_{n}, \hat{\Delta}, \alpha\right)$ is only monotonically increasing but not necessarily concave. When shape restrictions like concavity and monotonicity are applicable, shape-enforcing operators proposed by chen2021shape can be used to improve estimation and inference ex post. In particular, chen2021shape demonstrate that an estimator transformed by the shape-enforcing operator is more accurate, and confidence intervals transformed by the shape-enforcing operator have greater coverage and shorter length.

The double Legendre-Fenchel transform and the “pooled adjacent violators algorithm” as described in chen2021shape can be used to enforce concavity of the orthogonal estimator $\beta\left(\mathbb{Q}_{n}, \hat{\Delta}, \alpha\right)$. An alternative approach of concavification is via the Lagrangian dual. In particular, the model based estimator $\beta\left(\hat\mu, \hat{\theta}, \alpha\right)$ is equivalent to the Lagrangian dual

where $\hat{g}_{j}\left(\cdot\right) = g_{j}\left(\cdot, \hat{\theta}\right)$, $\hat{c}_{j}\left(\cdot\right) = c_{j}\left(\cdot, \hat{\theta}\right)$ and $j = 0, 1$. Since the Lagrangian dual minimization problem is affine with respect to $\alpha$, $\beta\left(\hat\mu, \hat{\theta}, \alpha\right)$ can be considered an implicit concavification of $\beta\left(\mathbb{Q}_{n}, \hat{\Delta}, \alpha\right)$. Both the double Legendre-Fenchel transform and the Lagrangian dual $\mathcal{L}\left(\hat{\mu}, \hat{\theta}, \alpha\right)$ can be generalized to the problem of multi-constraints.

Under misspecification, the asymptotic distribution of the parametric orthogonal estimator $\beta\left(\mathbb{Q}_{n}, \hat{\theta}, \alpha\right)$ is described in Appendix (ref), where the integrals with respect to Hausdorff measure in (ref) do not vanish in general.

In causal inference, $Y_j, Z_j, j \in \{0,1\}$ are not fully observed. Only $Y=D Y_1 + \left(1-D\)Y_0$ and $Z=D Z_1 + \left(1-D\)Z_0$ are observed, where $D \in \{0,1\}$ and unconfoundedness holds:

The unconfoundedness assumption (ref) allows the model based constrained value function (ref) to be identified using $g_1\(x\right) = \mathbb{E}\(Y \vert D=1, x\right)$, $g_0\(x\right) = \mathbb{E}\(Y \vert D=0, x\right)$, $c_1\(x\right) = \mathbb{E}\(Z \vert D=1, x\right)$, $c_0\(x\right) = \mathbb{E}\(Z \vert D=0, x\right)$. When the conditional data generating processes are parametric models of $f\(y, z \vert x, d; \theta\right)$, denoted as

a parametric model based estimator is given in (ref) such that $f\(y_1 \vert x; \theta\right)$, $f\(z_1 \vert x; \theta\right)$, $f\(y_0 \vert x; \theta\right)$ and $f\(z_0 \vert x; \theta\right)$ are the marginal densities implied by (ref), and such that $\hat\theta$ is the conditional MLE:

Under standard regularities, the asymptotic linear representation as in (ref) continues to hold with

In addition to the model based estimator (ref), an orthogonal doubly robust value function estimator under unconfoundedness can be formed by re-expressing (ref) as

In the above, for $j \in \{0, 1\}$, $p_{1}\left(\cdot\right) \coloneqq p\left(\cdot\right)$ and $p_{0}\left(\cdot\right) \coloneqq 1 - p\left(\cdot\right)$, we denote $\Gamma_{j}\(x, y, d\right) \coloneqq g_{j}\(x\right) + \frac{d}{p_{j}\(x\right)}\(y - g_{j}\(x\right)\right)$ and $\Xi_{j}\(x, z, d\right) = c_{j}\(x\right) + \frac{d}{p_{j}\(x\right)}\(z - c_{j}\(x\right)\right)$. Under correct specification, and assume that $\epsilon < p\left(\cdot\right) < 1-\epsilon$, for some constant $\epsilon > 0$, an adaption to Theorem (ref) shows that

where for $\phi^{\star}\(x; \alpha\right) \coloneqq 1\left(\Delta^{\star}\(x; k\(Q, \Delta^{\star}, \alpha\right)\right) > 0\right)$, we write

The ROC curve

A special case of (ref) is when $Y_{1}=Y \in \{0, 1\}$, $Y_{0} = Z_{0} = 0$, $Z_{1} = 1 - Y_{1}$, $g_{0}\(x\right) = c_{0}\(x\right) \equiv 0$, $g_{1}\(x\right) = p\(x\right)$, $c_{1}\(x\right) = 1 - p\(x\right)$ where $p\(x\right) = \mathbb{E}\(Y \vert\ X = x\right)$. This special case, i.e.

is closely related to the population receiver operating characteristic (ROC) curve defined by

The objective in (ref) is to maximize the power $\beta$ while controlling the size of the test $\alpha$. This optimization problem is exactly the same as that of the Uniformly Most Powerful Tests under simple null and alternative hypotheses. See for example section 3 of lehmann2022testing for rigorous definitions and implications of these statistical testing concepts. In this paper, we adopt the more commonly employed terminology in the computer science, biostatistics and economics literatures to refer to $\beta$ and $\alpha$ as true positive rate (TPR) and false positive rate (FPR). In berge2011evaluating, $Y_1 = Y$ is interpreted as the state of the economy ($Y = 0$ for expansion and $Y = 1$ for recession). The ROC curve is also used in the finance loan literature vallee2019marketplace,sadhwani2021deep to evaluate the tradeoff between earned interest and lost principal, where $Y=0$ indicates loan fulfillment and $Y=1$ indicates loan default. The abnormal birth detection application of feng2025statisticalMS uses $Y$ as the indicator of abnormal birth ($Y = 0$ for normal birth and $Y = 1$ for abnormal birth).

The analog of (ref) for the constraint and value functions of (ref) is

In order to adapt Theorem (ref) to the ROC curve in (ref) and (ref), define similarly to (ref),

Different from (ref), (ref) has denominators and $\Delta\(x; k\right)$ is replaced by $p\(x\right) - k$. We also define

The following corollary adapts Theorem (ref) to the commonly used orthogonal estimator of (ref).

Under correct specification, $p^\star\left(\cdot\right) = p\left(\cdot\right)$ and $\frac{ k }{ 1- k } = \frac{Q y}{1- Qy} \frac{f_{\beta}\(Q, p^{\star}, k\right)}{f_{\alpha}\(Q, p^{\star}, k\right)}$. This is because, as $\left(1-k\right) p\(x\right) - k \left(1-p\(x\right)\right) = 0$ on $p^\star\(x\right) = p\(x\right)= k$,

Furthermore, the two integrals on the level set in Corollary (ref) cancel out since

Under correct specification (ref) simplifies to

In parametric models of $p\(x\right)$, where $r_n = \sqrt{n}$ and $t_Q = 1$, (ref) has a Gaussian covariance kernel

Similar to (ref) and (ref), a model based representation of the ROC curve (ref) is

The corresponding estimator is

Typically, $\hat\theta = \arg\max_{\theta\in\Theta} \sum_{i=1}^n \[y_i \log p\(x_i; \theta\right) + \left(1-y_i\right) \log \left(1-p\(x_i; \theta\right)\right)\right]$ is the maximum likelihood estimator.

When the parametric model is correctly specified, the model based estimator is more efficient than the orthogonal estimator, as formalized in the next proposition and illustrated in a synthetic data example in section (ref).

Misspecified parametric models estimated by the orthogonal estimator are discussed in Appendix (ref). The next Corollary (ref) validates bootstrap consistency under general specification. For correctly specified parametric models, the first step model only needs to be trained once and does not need to be re-estimated for every bootstrapped sample.

Let $M = \(M_{1}, \ldots, M_{n}\right)$ be the vector of bootstrap weights. For a random element $\tilde{\mathbb{Z}}_{n} \coloneqq \tilde{\mathbb{Z}}_{n}\left(\{\(X_{i}, Y_{i}\right)\}^{n}_{i = 1}, M\right)$ on normed space $\mathcal{X}$, we say its bootstrap law is consistent for a tight process $\mathbb{G}$ and denote it as $\tilde{\mathbb{Z}}_{n} \overset{\mathbb{P}}{\underset{\mathbb{W}}{\rightsquigarrow}} \mathbb{G}$ if in outer probability\footnote{ Under measurability, the consistency definition can be simplified to $\sup_{l \in BL_{1}\left(\mathcal{X}\right)} \left\vert \mathbb{E}_{M}l\left(\tilde{\mathbb{Z}}_{n}\right) - \mathbb{E}l\left(\mathbb{G}\right)\right\vert = o_{\mathbb{P}}\left(1\right).$ }

In the above $\mathbb{E}_{M}$ denotes the conditional expectation with respect to $M$ given the data $\{\(X_{i}, Y_{i}\right)\}^{n}_{i = 1}$, and $BL_{1}\left(\mathcal{X}\right)$ denotes the set of real-valued functions on $\mathcal{X}$ with Lipschitz constant and supremum norm both bounded by $1$.

Consider i.i.d. data $\{\(Y_i, W_i, X_i\right), i=1,\ldots,n\}$, where $Y_i$ is the true outcome, $W_i$ is the human decision, and $X_i$ are observable features. An estimator of the human TPR/FPR pair is

Human decision making quality can be evaluated by testing whether the population TPR/FPR pair of a human decision maker, denoted $\left(\alpha_{H}, \beta_{H}\right)$, lies on, above or below a population ROC curve $\beta\left(\alpha\right)$ of a machine learning model. In the above,

The null and alternative hypotheses can be formulated as

where $H_A$ can be further categorized into two one-sided directional alternatives $H_{A,h}: \beta\left(\alpha_{H}\right) < \beta_{H}$ and $H_{A,m}: \beta\left(\alpha_{H}\right) > \beta_{H}$, as in vuong1989likelihood. It follows from the Lindeberg-L\'{e}vy central limit theorem that for $\theta_{H} = \left(\alpha_{H}, \beta_{H}\right)'$, $\hat{\theta}_{H} = \left(\hat{\alpha}_{H}, \hat{\beta}_{H}\right)'$ and $p = Q Y$,

When the ROC $\beta\left(\alpha\right)$ is fully known,\footnote{ For example, feng2025statisticalMS assume a known ROC because the sample size of the data set used to train the machine learning model is by orders of magnitude larger than the sample size for each individual decision maker. } the hypothesis (ref) can be tested using the statistics $\beta\left(\hat{\alpha}_{H}\right) - \hat{\beta}_{H}$. By the delta method, under the null hypothesis $\beta\left(\alpha_{H}\right) = \beta_{H}$ in (ref),

A more complex situation is when the machine ROC curve $\hat\beta\left(\alpha\right)$ and the human TPR/FPR pair $\left(\hat{\alpha}_{H}, \hat{\beta}_{H}\right)$ are estimated from the same dataset. In such a case, a revised test statistic for the hypothesis in (ref) is then $\hat\beta\left(\hat{\alpha}_{H}\right) - \hat{\beta}_{H} = \hat{\beta}\left(\hat{k}\left(\hat{\alpha}_{H}\right)\right) - \hat{\beta}_{H}$, where $\hat\beta\left(\cdot\right)$ can be any one of orthogonal ROC $\beta\left(\mathbb Q_n, \hat p\left(\cdot\right), \cdot\right)$, parametric model based ROC $\beta\left(\hat\mu, \hat\theta, \cdot\right)$ and parametric orthogonal ROC $\beta\left(\mathbb Q_n, \hat \theta, \cdot\right)$. Similarly, $\hat{k}\left(\cdot\right)$ can be the corresponding orthogonal $k\left(\mathbb Q_n, \hat p\left(\cdot\right), \cdot\right)$, parametric model based $k\left(\hat\mu, \hat\theta, \cdot\right)$ and parametric orthogonal $k\left(\mathbb Q_n, \hat \theta, \cdot\right)$ in (ref) and (ref). The following Hadamard differentiation result applies to these estimators.

When $p\(x\right)$ is modeled parametrically as $p\(x; \theta\right)$ and $\sqrt{n}\left(\hat{\theta} - \theta^{\star}\right)$ has an asymptotic linear representation by $\kappa\(x, y\right)$, by the delta method, $H\(x\right)$ will be determined by $\frac{\partial}{\partial\theta}p\(x; \theta^{\star}\right)$ and $\kappa\(x_i, y_i\right)$ giving rise to an asymptotic linear representation. In each case, the asymptotic variance can be consistently estimated by sample analogs. Bootstrap also provides valid inference. The asymptotic distribution and verification of bootstrapped confidence interval of $\sqrt{n}\left[\left(\hat{\beta}\left(\hat{\alpha}_{H}\right) - \hat{\beta}_{H}\right) - \left(\beta\left(\alpha_{H}\right) - \beta_{H}\right)\right]$ is obtained by combining the standard term $\hat{\beta}_{H} - \beta_{H}$ with $\hat\beta\left(\hat{\alpha}_{H}\right) - \beta\left(\alpha_{H}\right)$ handled by (ref). Section (ref) showcases a synthetic data example.

Fast convergence rates

In this section we build on the previous analysis to generalize the fast learning rates for plug-in classifiers in devroye1996probabilistic and audibert2007fast to multiclassification and multi-discrete allocation problems. The fast learning rate phenomena can also be found in recent works by semenova2023debiased and liu2024inference. Following chernozhukov2018double, chernozhukov2022locally and foster2023orthogonal, Neyman-orthogonalization and sampling splitting can be combined with the plug-in scheme to facilitate statistical inference.

Second order Hadamard differentiability

Our generalizations are motivated by the partial convexity of the social welfare potential function $\gamma\left(\cdot, \cdot\right)$. Consider, for a given $g\left(\cdot\right)$,

By construction, $e\left(\lambda'; \phi^{*}\left(\cdot; \lambda, g\right)\right)$ achieves its maximum of zero at $\lambda'=\lambda$. When there is no ambiguity, we also denote $\phi^{*}\left(\cdot\right) \coloneqq \phi^{*}\left(\cdot; \lambda, g\right)$. Continuous convex function is locally Lipschitz and (locally) Lipschitz function is differentiable almost everywhere by the Rademacher theorem (ref). Due to the convexity of $\gamma\left(\lambda', g\right)$ in $\lambda'$ given $g\left(\cdot\right)$, $\gamma\left(\lambda', g\right)$ is differentiable in $\lambda'$ given $g\left(\cdot\right)$ in a set $\Lambda$ such that its complement $\Lambda^c$ is a Lebesgue null set. Then for each $\lambda \in \Lambda$ and for each selection $\phi^{*}$ given this choice of $\lambda$, $e\left(\lambda'; \phi^{*}\right)$ is differentiable in $\lambda'$ at $\lambda'=\lambda$. Since zero must be in the set of supgradients at a point of maximum of any function, there is $0 \in \nabla^+ e\left(\lambda'; \phi^{*}\right)\big\vert_{\lambda'=\lambda}$ and thus $\nabla e\left(\lambda'; \phi^{*}\right)\big\vert_{\lambda'=\lambda} = 0$, implying that

We expect such envelope theorem like statement remains true with functional derivatives of

at $\left(\lambda', g'\right) = \left(\lambda, g\right)$ (here, $\phi^{*}\left(\cdot\right) \coloneqq \phi^{*}\left(\cdot; \lambda, g\right) $). Theorem (ref) shows Fr\'{e}chet differentiability and verifies (ref). Generically, under weak assumptions,

By construction, the remainder term of $o\left(\Vert \theta\Vert\right)$ in the above expression is

The magnitude of this term can be investigated by relating the Hadamard differentiability framework in section (ref) to the margin assumption (MA) in the plug-in classifier literature.

Consider the estimation of social welfare potential function under full observability. Let $\mu$ denote the distribution of $X$ and $Q$ the joint distribution of $(y,x).$ Define the following notation:

To contrast the application of measure $\mu$ with that of measure $Q$, we first consider the perturbation of $\gamma\left(\lambda, g(\cdot),\mu\right)$ to $\gamma\left(\lambda, g(\cdot)+t_n\,H_n(\cdot),\mu+t_n{\mathcal P}_n\right)$ where $t_n$ approaches zero and $H_n(\cdot)$ and ${\mathcal P}_n$ provide directions of the perturbation of the conditional expectation function $g(\cdot)$ and of the measure $\mu\left(\cdot\right)$ over ${\mathcal X},$ respectively. We consider the following decomposition:

The first term $\left(1\right)$ is the second order difference “empirical process” term that should vanish under suitable continuity conditions. The third term $\left(3\right)$ is a standard “summation of random variables” term. The second term $\left(2\right)$ can be further decomposed as

where for brevity we omit the argument of functions $g(\cdot)$ and $H(\cdot)$. If we use $\gamma\left(\lambda, g, Q\right)$ instead of $\gamma\left(\lambda, g, \mu\right)$, then we will not have \textcircled{{2}} and \textcircled{{3}} terms. In this sense, using $Y$ instead of $g$ provides a natural Neyman-orthogonalization of the estimator. We show that under Hadamard differentiability, the term \textcircled{{1}} is $o\(t_n\right)$, leading to a simplified form of the Hadamard derivative of $\gamma\left(\lambda, g, Q\right)$.

Compare Lemma (ref) with Theorem (ref). By (ref) in Theorem (ref), the first order influence of pertubing the second variable $g$ will not come from the indicator term. Equation (ref) in Lemma (ref) derives an analogous Hadamard differentiability result. Neyman-orthogonalization changes the remainder term (ref) to $\mathbb{E}\left[\sum_{j}\lambda_{j}\left(\phi^{*}_{j}\(X; \lambda, g + \theta\right) - \phi^{*}_{j}\(X; \lambda, g\right)\)g_{j}\(X\right)\right].$ Next, we show under a strict MA in (ref) with $\alpha > 1$, this term is second order degenerate.

In the following we adopt a convention for the optimal policy in (ref) such that all the weight among the maximal indexes is allocated to the smallest index member:

We introduce this notation in particular to handle the case $\lambda_{j}\hat{g}_{j} = \lambda_{l}\hat{g}_{l}$ when $\hat{g}$ is estimated. For $k\neq l, k,l=0,\ldots,J$, we also write

The main drawback of Hadamard differentiability (functional delta method) is that it does not fully utilize the convergence rate of the first step estimator. Since Proposition (ref) guarantees $\alpha = 1$, the requirement of $\alpha > 1$ in Proposition (ref) is less satisfactory. In section (ref) we overcome this drawback using sample splitting and Neyman orthogonality. In the following we provide two additional results to show more connection between the theory in section (ref) and the MA.

Recall that the push-forward measure $T_{\#}\mu$ generated by a Borel measurable map $T$ from the measure $\mu$ is defined by $T_{\#}\mu\[A\right] = \mu\[T^{-1}\(A\right)\right]$. The original margin assumption in (ref) can be extended to the following:

A definition of a margin controlled random variable is to impose local H\"{o}lder-continuity on its distribution function with a universal exponent $\alpha$.

Multi discrete allocation problem

In this section we apply the fast convergence rate to social welfare functions with multiple treatments, generalizing the arguments in devroye1996probabilistic and audibert2007fast.

Recall that in a general discrete allocation problem, we define welfare under policy $\phi\left(\cdot\right)$ as

where $g_j\(x\right) = \mathbb{E}\(Y_j \vert X = x\right)$ and $\phi_j\(x\right) \geq 0, \forall j$, $\sum_{j=0}^J \phi_j\(x\right) = 1$. A special case of (ref) is the weighted multiclassification accuracy when $Y_j \in \{0,1\}$, $\sum_{j=0}^J Y_j =1$, and when $\phi_j\(x\right), j=0,\ldots,J$ is a multiclassification allocation.\footnote{The binary treatment case in audibert2007fast is a special of (ref) where $g_1\(x\right)=p(x)$, the propensity score, $g_0\(x\right)=1-p\(x\right)$, $\phi_1\(x\right) = 1\(p(x)>\tau\right)$, $\phi_0\(x\right) = 1\(p(x) \leq \tau\right)$, $\lambda_1=\tau$ and $\lambda_0=1-\tau$.}

In classification problems, $Y_j, j=0,\ldots,J$ are directly observed. In contrast, in optimal treatment allocation problem, $Y_j$ is typically observed only when $d_j=1$, where $d_j$ is the indicator of whether the $j$th treatment is taken, such that $d_j \in \{0,1\}$ and $\sum_{j=0}^J d_j = 1$. We assume the standard unfoundedness condition, i.e., $Y_{j}\perp\!\!\!\perp D_{j} \vert X, \forall j=0,\ldots,J$. Then $g_j\(x\right) = \mathbb{E}\(Y_j \vert D_j =1, x\right)$ can be estimated based on the subsample where $d_j = 1$.

The envelope-like Theorem (ref) essentially states that generically there is no first order impact of estimating $\phi^*\left(\cdot; \lambda, g\right)$ in $\gamma\left(\lambda, g\right) = \mathbb{E} \left[\sum_j \phi_j^*\(X; \lambda, g\right) g_j\(X\right)\right]$ by $\phi^{*}\left(\cdot; \lambda, \hat{g}\right)$. The fast convergence rate of regret to zero is a consequence of the second order effect of estimating $\phi^*\left(\cdot; \lambda, g\right)$ by $\phi^{*}\left(\cdot; \lambda, \hat{g}\right)$. We use the same convention as (ref) in the following expressions for the optimal and feasible policies:

The regret of using the feasible policy for welfare is:

By definition, $\phi_j^{*}\(x\right)$ is either $0$ or $1$. Furthermore, for each $x$, one and only one of $\phi_j^{*}\(x\right)$ out of $j=0,\ldots,J$ takes the value $1$. The same is true for $\hat\phi_j\(x\right)$ for each $x$. Therefore we can write

where the second inequality uses the simplified notation (ref). The first inequality of (ref) is due to the fact that when $\phi_i^{*}\(x\right) = \hat\phi_j\(x\right) =1$,

Following Example 1 in chen2003estimation, for $\delta_n = \sup_{x\in\mathcal X} \max_{j=0,\ldots,J} \vert \lambda_j g_j\(x\right) - \lambda_j \hat{g}_j\(x\right) \vert$,

Then we can further bound the regret as

The conditions in section (ref) imply that for $k\neq l$, $\lambda_k g_k\(x\right) - \lambda_l g_l\(x\right)$ has a bounded density at zero, leading to $\mathbb{P}\left(-2 \delta_n \leq \lambda_k g_k\(X\right) - \lambda_l g_l\(X\right)\leq 2 \delta_n\right) \leq C \delta_n$ and $\gamma\left(\lambda, g, \phi^{*}\right) -\gamma\left(\lambda, g, \hat\phi\right) \leq C \delta_n^2.$

Nonparametric estimators $\hat{g}_j\(x\right)$ typically achieve the convergence rate $\delta_n = O_{\mathbb{P}}\(n^{-\frac{\beta}{2\beta+d}} \log n\right)$ where $d$ is the dimension of $x$, and $\beta$ a smoothness parameter, implying that

Whenever $\beta > d/2$, the convergence rate of $\gamma\left(\lambda, g, \phi^{*}\right) -\gamma\left(\lambda, g, \hat\phi\right)$ becomes $O_{\mathbb{P}}\left( n^{-\frac{1}{2}}\right).$

Inference on social welfare potential function

The unknown social welfare potential function $\beta_0 = \gamma\left(\lambda, g\right) = \mathbb{E}\left[\sum_{j}\lambda_j \phi^{*}_{j}\(X\)g_{j}\(X\right)\right]$ can be estimated by its sample analog $\frac{1}{n'} \sum_{i=1}^{n'} \sum_{j} \lambda_j \hat{\phi}_{j}\(X_{i}\right)\hat g_j\(X_{i}\right)$ for $n' = O\(n\right)$. On the one hand, under full observability, the Neyman-orthogonalized estimator $\frac{1}{n'} \sum_{i=1}^{n'} \sum_{j} \lambda_j \hat{\phi}_{j}\(X_{i}\)Y_{ij}$ can be combined with a sample splitting scheme for estimating $\hat{\phi}_{j}\(X_{i}\right)$ to achieve desirable statistical properties. On the other hand, for a typical treatment allocation problem, $Y_{ij}$ are not fully observable. We only observe $Y_{i} = \sum_{j}D_{ij}Y_{ij}$ and need to use the treatment status variables $D_{ij}$ to construct a Neyman-orthogonalized estimator.

We note that luedtke2016statistical and luedtke2020performance utilize a formula similar to (ref) for the binary case under empirical process settings. To apply broadly to machine learning methods, Theorem (ref) and Theorem (ref) combine the fast convergence rate of a plug-in classifier with the insighful techniques of chernozhukov2018double, chernozhukov2022locally and foster2023orthogonal based on Neyman orthogonalization and sample splitting. Inspired by semenova2023debiased, the asymptotic distributions of (ref) and (ref) can be estimated through bootstrap.

Simulation and empirical application

Confidence intervals and hypothesis testing of ROC curves

To evaluate the finite sample performance of our asymptotic results on the ROC curves, we first study the simplified asymptotic distribution of the ROC curves under correct specification. Next, we consider a synthetic data example and illustrate an important application of comparing a machine generated ROC curve with a human decision maker.

In our first simulation, we generate data according to a logit model:

The Epanechnikov distribution with the parameter $c = 1$ has the density function

We examine the empirical coverage frequencies and average lengths of three types of confidence intervals (hereafter abbreviated as CI) for the ROC curve of the correctly specified logit model. The first $CI_{bootstrap}$ is formed via the orthogonal estimator using the maximum likelihood estimator (MLE) of the logit model and the bootstrapping method in Corollary (ref). The second $CI_{theoretical}$ is formed via the MLE orthogonal estimator and the analytic asymptotic distribution in (ref). The third $CI_{model \; based}$ is formed by bootstrapping the model based ROC curve described in Proposition (ref). To calculate the $CI_{bootstrap}$ under correctly specification we do not need to reestimate the logit model for each bootstrap sample. The constructions of the ROC estimators and three types of CIs are summarized in Algorithm (ref).

We use $R = 1000$ bootstrap repetitions each in 1000 Monte Carlo simulations for sample sizes $n \in \{1000, 2000, 5000, 10000\}$. We evaluate the pointwise CI at FPR level between $0.3$ and $0.7$ with stepsize $0.05$. Since the analytic expression for the population ROC curve is difficult to derive, we calculate an approximation by simulating a large sample size of $5 \cdot 10^{6}$ from the data generating model. As shown in Table (ref), the coverage is very close to the nominal level $95\%$. The average lengths of $CI_{bootstrap}$ and $CI_{theoretical}$ are nearly the same, while $CI_{model \; based}$ is saliently the shortest.

Next, we consider a problem of comparing algorithms with decision making of “doctors” using synthetic data. Consider i.i.d. observations $\{\(Y_i, W_i, X_i\right), i = 1, \ldots, n\}$, such that $Y_i$ is the true outcome, $W_i$ is the doctor's diagnosis, and $X_i$ are observable features. To compare an algorithm with a doctor, we test the null hypothesis that the human TPR/FPR pair lies on the machine ROC curve against the alternative hypotheses that the human pair lies above or below the ROC curve. We generate synthetic data according to the following model

where $I_{4}$ is the $4$-dimensional identity matrix. The $W^{E}_{i}$ decisions are made by an ideal “experienced doctor” who knows the data generating process, while the $W^{I}_{i}$ decisions are made by an “inexperienced doctor” who only utilizes $X_{1}$ and $X_{3}$. Theorem (ref) is applicable to this model.

A human decision maker (e.g. doctor, judge and admission officer) is unlikely to make a large number of decisions. We set the simulation sample size $n$ to $\{100, 250, 500, 1000, 2000, 5000, 10000\}$. We first use Algorithm (ref) with $R = 1000$ to form a CI for the difference $\beta\left(\alpha_{doc}\right) - \beta_{doc}$. We then determine whether $0$ is within, above or below the interval. In the second or third cases, we reject the hypothesis of equal performance between the doctor and the machine model in favor of either the doctor or the model.

Table (ref) shows the empirical rejection frequencies for nominal 5% and 10% tests. In the first two entries of both panels, the “experienced doctor” and the “inexperienced doctor” are compared with a misspecified logit model using only $\(X_{1}, X_{2}, X_{3}\right)$. Additionally, in the last entry of both panels labelled “model-like case”, we compare the same “experienced doctor” with a correctly specified and estimated logit model using all of $\(X_{1}, X_{2}, X_{3}, X_{4}\right)$. When the sample size is only $100$, the test has low power but rarely rejects in the wrong direction. A high power level is achieved when the sample size reaches $500$ or $1000$. The “model-like case” represents the null hypothesis of equal performance between the “experienced doctor” and the model. As shown in Table (ref), the empirical size of the test is close to the nominal level. The simulation results validate the finite sample quality of our asymptotic approximation theory.

Real data application and simulation

In this section, we apply the plug-in multi-class classification method in section (ref) to a voting dataset of the August 2006 primary election in Michigan. This dataset originates from the field experiment reported in gerber2008social.\footnote{We download the dataset from \url{https://github.com/gsbDBI/ExperimentData/tree/master/Social}.} The optimal voting incentive policy has been studied in the preprint version of zhou2023offline using their decision-tree search method.\footnote{See the preprint version of zhou2023offline at \url{https://arxiv.org/pdf/1810.04778}.}

The voting dataset contains 180002 electorates from households spanning Michigan. Following zhou2023offline, we use $10$ features in the dataset: year of birth, sex, household size, city, turnout for 2000, 2002 and 2004 primary and general elections. We experiment with applying both the multilayer perceptron (MLP) algorithm and the causal forest (CF) algorithm (developed in wager2018estimation and athey2019generalized) to the voting data set.

There are four treatment groups and a control group in the voting field experiment. Hereafter, the control group will be called Control, and the four treatment groups will be called Civic Duty, Hawthorne, Self and Neighbors. See gerber2008social for the specific meaning of those treatments. Electorates are assigned to the five groups at random with probabilities $\left(\frac{5}{9}, \frac{1}{9}, \frac{1}{9}, \frac{1}{9}, \frac{1}{9}\right)$.

First, we generate a synthetic dataset similar to the voting experiment in order to validate the multiplier bootstrapping method described in Corollary (ref). The data generating process is a MLP model with $1$ hidden layer of size $8$, ReLU activation and a softmax function (i.e. standard logistic function). The outcome is then generated according to

where the treatment assignment probabilities are independent of the characteristics $X_i$. We let $X_{i} \sim N\left(\mu, \Sigma\right)$, where $\mu$ and $\Sigma$ are specified in Appendix (ref). The weights $W_{1}$, $W_{2}$, $b_{1}$ and $b_{2}$ are randomly drawn. The values of the weights matrices are also reported in Appendix (ref).

Using this data generating process, we first draw a large sample of size $10^{7}$ to compute a close approximation of the social welfare potential $\beta_{0}$ around $0.311$ . An overparameterized MLP with $2$ hidden layers of size (48, 24) is used in estimation. We then use the multiplier bootstrap method ($1000$ repetitions) in section (ref) to construct CI based on a sample splitting scheme of $20000$ training and $5000$ testing observations. To the best of our knowledge, there is no current guarantee of the almost sure uniform convergence rate for neural networks. However, inspired by the theory in bos2022convergence, it is reasonable to expect a fast uniform convergence rate of the MLP classifier on a large probability set of $X_{i}$. In $1000$ simulations, our multiplier bootstrap CI has an empirical coverage rate of $0.939$, close to the nominal $0.95$ level. The average length of the CI is around $0.0423$.

Next, two machine learning methods are applied to this data. The first method is a MLP with three hidden layers of size $\left(64, 64, 64\right)$ and ReLU activation (i.e. $\max\{0, x\}$). The second method is a CF with $500$ trees where the minimum number of samples for each group in each leaf is equal to $600$. We use the Python package scikit-learn scikit-learn for the MLP and the Python implement from the EconML package econml for the CF. The treatment group of each electorate is concatenated to the MLP input features as a $1$-dimensional variable taking value $0$ to $4$ (in the order listed above).

We follow zhou2023offline and use $5$-fold cross validation. We report the estimated difference and CI of the difference between the plug-in policy and the control group based on (ref). Note that (ref) is Neyman orthogonal and amenable to a sample splitting scheme described in Algorithm (ref) with $R = 1000$ and $K = 5$.

Both the MLP and the CF assign nearly all electorates to the Neighbors group. Meanwhile, the two models give value differences of $0.0865$ and $0.0837$ compared to control. The CIs for the true $\beta_{0}$ are reported in Table (ref). These findings are similar to the value difference around $0.0864$ using the ERM policy of zhou2023offline. Their preprint version reported a value difference of $0.082$ for the same ERM policy. The difference is due to the randomness from sample splitting and model training. The policy gain results need to be interpreted with caution. As shown in Table (ref), the CIs for the policy gain between the MLP / CF and the ERM policy include $0$. Given that the performance of our proposed method is statistically on par with ERM, the method that we propose provides a viable computationally practical alternative approach to ERM policy learning.

Conclusion

In this paper, we explore a functional differentiability approach for a class of statistical optimal allocation problems. Inspired by chernozhukov2018sorted, Hadamard differentiability is facilitated by a study of the general properties of the sorting operator. We provide a derivation that makes use of the concept of Hausdorff measure and the area and coarea formulas from geometric measure theory. Based on our general Hadamard differentiability results, in section (ref) we derive the asymptotic properties of the plug-in estimators for both the value function process of a binary constrained optimal allocation problem and the ROC curve from the functional delta method. When the first step propensity score model is correctly specified, a computationally feasible bootstrap procedure is validated for the plug-in ROC estimator.

Importantly, we build on the partial convexity of the social welfare potential function, i.e. the value function of the optimal allocation problem, to demonstrate the degeneracy of the first order derivative of the social welfare potential function with respect to the policy. These intriguing results provide us with insights when we combine techniques from the literature of nonsmooth method of moment, plug-in classification and the recent development of double / debiased machine learning to develop a debiased estimator of the social welfare potential function in section (ref). Here, the conditions required for Hadamard differentiability validate the margin assumption, which leads to a faster convergence rate.

Acknowledgments

We thank the editor and two anonymous referees for insightful comments. We also acknowledge helpful and encouraging conversations with collaborators and colleagues including Chunrong Ai, Timothy Armstrong, Haoge Chang, Qihui Chen, Xiaohong Chen, Victor Chernozhukov, Michael Fan, Yanqin Fan, Yue Fang, Yingjie Feng, Robin Han, Toru Kitagawa, Michael P. Leung, Hongjun Li, Jessie Li, Ruixuan Liu, Ye Luo, Thomas MaCurdy, Chen Qiu, Shuyang Sheng, Zhentao Shi, Kurt Sweat, Liangjun Su, Ke Tang, Guanyi Wang, Haitian Xie and Ping Yu, as well as participants at various seminars and conferences. We are most grateful to Yuhao Xue for enlightening us with his knowledge of geometry measure theory and to Yichuan Zhang for research assistance. We acknowledge funding support from the National Science Foundation (SES 1658950 to Han Hong).

\FloatBarrier \phantomsection \addcontentsline{toc}{section}{References}

\setstretch{1.25} {{