Abstract

When estimating causal effects from observational studies, researchers often need to adjust for many covariates to deconfound the non-causal relationship between exposure and outcome, and often many covariates are discrete. The behavior of commonly used estimators in the presence of many discrete covariates is not well understood, since standard approaches often employ structural assumptions such as smoothness, which do not apply in discrete settings. In this work, we study estimation of causal effects in a model where the covariates required for confounding adjustment are discrete but high-dimensional, meaning the number of categories \(d\) can be comparable to or even larger than sample size \(n\). Specifically, we show the mean squared error of commonly used regression, weighting and doubly robust estimators is bounded by \(\frac{d^2}{n^2}+\frac{1}{n}\). We then prove that the minimax lower bound for estimating the average treatment effect is of order \(\frac{d^2}{n^2 \log^2 n}+\frac{1}{n}\), which characterizes the fundamental difficulty of causal effect estimation in the high-dimensional discrete setting, and shows the estimators mentioned above are rate-optimal up to log factors. Finally we consider two other kinds of structure that can be exploited: effect homogeneity, and prior knowledge of the covariate distribution. We propose new estimators that enjoy faster convergence rates here, of order \(\frac{d}{n^2} + \frac{1}{n}\), thus achieving consistency in a broader regime. The results are illustrated empirically via simulation studies and a real data example.

Keywords: average treatment effects, categorical data, effect homogeneity, high-dimensionality, minimax lower bounds.

1 Introduction↩︎

To draw causal conclusions from observational studies, researchers typically have to measure and adjust for many covariates (e.g., any confounders that could affect treatment assignment or the outcome of interest). Numerous approaches exist for estimating causal effects from such data. Common methods include outcome modeling [1], [2], inverse propensity score weighting [3]–[5], and semiparametric methods (i.e., doubly robust estimation, or augmented inverse probability weighting, or targeted or double/debiased machine learning) [6]–[8]. Notably, the latter semiparametric methods are consistent for the average treatment effect (ATE) even when either the outcome regression or propensity score is misspecified, yielding a robust tool for estimating ATEs under possible model misspecification. More generally, the bias of these methods is a product of errors in outcome regression and propensity score estimation, thus allowing parametric rates of convergence (and semiparametric efficiency) even when nuisance functions are estimated at slower nonparametric rates. Sample splitting or cross-fitting [7], [9], [10] are often used to prevent overfitting the nuisance functions and avoid empirical process conditions [8]. In the classic setting where the covariate dimension is held fixed, while the sample size grows to infinity, the aforementioned methods can enjoy favorable properties (e.g., \(\sqrt{n}\)-consistency and asymptotic normality) and performance in applications. However, when the number of covariates needed for confounding adjustment is large, some new and important problems arise.

Indeed, recently there has been much work in the causal inference literature focused on understanding the behavior and properties of various estimators in the high-dimensional setting. [11] derived the excess variance of commonly used estimators due to first-stage nuisance estimation with high-dimensional covariates by assuming a linear outcome model and known propensity score. They also illustrated the inflation of variance in simulations, showing even the doubly robust estimator may not achieve the efficiency bound when the covariates included are high-dimensional. [12] established a novel central limit theorem for the doubly robust estimator, when the outcome regression model and propensity score follow high-dimensional generalized linear models, without assumptions on sparsity, but under a stylized setting where covariates are normally distributed and the covariate dimension \(d\) is of the same order as sample size \(n\). They also showed estimates obtained by permuting the folds in cross-fitting are asymptotically correlated in the high-dimensional regime. [13] proposed a novel debiased method for missing data models in the \(n < d\) setting, where ordinary least square estimation is not feasible. However, to help overcome the curse of dimensionality and achieve non-trivial rates of convergence, most of the work in this literature combines the debiased machine learning/doubly robust techniques with additional structural assumptions on the nuisance functions, such as the aforementioned linearity, or smoothness [14]–[16], or sparsity [7], [17]–[19]. We refer to [20]–[28] and others for more relevant discussion on causal inference with high-dimensional data.

Our work enriches the high-dimensional causal literature under a different structural assumption, i.e., that the covariates are discrete. This setting is of interest for several reasons. First, it can help uncover new phenomena in more general high-dimensional regimes where the dimension \(d\) can be comparable to or larger than sample size \(n\), and motivate further exploration. The discrete covariate setting can also be viewed as an interesting base case, with crucial implications for continuous data, as has been seen in semiparametric efficiency bounds [29], bandit problems [30], and hypothesis testing [31], for example. Further, discrete high-dimensional covariates often arise in practice, for example in applied medical and health policy research, e.g., when adjusting for International Classification of Diseases (ICD) codes [32]. ICD codes include over 100,000 indicators, so their cartesian product, together with other common demographic covariates (e.g., sex, race, education level), induces a huge number of categories, certainly corresponding to a high-dimensional regime. The structured discrete data, including graphs, texts and images, is also ubiquitous in the machine learning literature, which commonly has high-dimensional representations. It is often important to take these non-numerical variables into account when estimating treatment effects [33], [34]. Our work establishes theoretical foundations for adjusting for high-dimensional discrete covariates and helps practitioners evaluate the reliability of their estimates given the sample size and covariate dimension.

Our work is also related to [35]–[39], and others, who considered estimating simpler functionals (e.g., entropy, support size) in high-dimensional discrete models. In this work, important breakthroughs have been made in terms of both improved estimation methodologies and minimax fundamental limits. Using tools in polynomial approximation theory [40], it has been established that a best-polynomial-approximation estimator usually enjoys the so-called “effective sample size enlargement" property, meaning its behavior with \(n\) samples resembles that of an MLE (the naive plug-in estimator) with \(n \log n\) samples. Lower bounds are established by a moment matching technique. Upper bounds and lower bounds using moment matching usually coincide here since, in the optimization sense, moment matching is the dual problem of best polynomial approximation where strong duality holds. The readers are referred to [38], [41] for more detailed discussion on this duality phenomenon. Estimation of causal effects under the discrete covariate framework has remained largely unexplored, but our work helps bridge this gap.

Specifically, in this paper we study the treatment effect estimation problem when adjusting for discrete covariates with possibly more categories than samples. Our four main contributions are summarized as follows:

First, we provide finite-sample bounds on the mean squared error of commonly used regression, weighting and doubly robust estimators of ATE estimation. Our results imply \(d/n \rightarrow 0\) is a sufficient and necessary condition for these estimators to be consistent under positivity assumption, where \(d\) is the number of categories and \(n\) is the sample size.
We then characterize the fundamental limits of ATE estimation under a positivity assumption. The minimax lower bound is (in terms of mean squared error) of order \(\frac{d^2}{n^2 \log^2 n} + \frac{1}{n}\), which shows that commonly used estimators are minimax optimal up to log factors. Moreover, \(d=o(n \log n)\) is a necessary condition for consistent ATE estimation.
Next, we explore the role of effect homogeneity in ATE estimation, showing that faster rates are achievable when the treatment effects are more homogeneous across categories. In fact, here consistent estimation is possible if \(d/n^2 \rightarrow 0\).
Finally we study estimation of treatment effects given prior knowledge of the covariate distribution. We first provide a negative result, showing the covariate distribution may not help improve the rate for ATEs. We then consider a variance-weighted average treatment effect, and derive its faster estimation rate, achieved by a second-order estimator. As in the homogeneous effects setting, our second-order estimator is consistent here if \(d/n^2 \rightarrow 0\).

The structure of our paper is as follows: After introducing the data-generating process, causal assumptions and notation in Section 2, in Section 3 we study the properties of commonly used regression, weighting and doubly robust estimators in estimating the average treatment effects (ATE). We first show the numerical equivalence between these three estimators and the plug-in estimator. We then provide finite-sample bounds on their mean squared error and characterize the sufficient and necessary conditions for them to be consistent in Section 3.1. In Section 3.2 we study the minimax lower bound and sample complexity of ATE estimation in the high-dimensional setting, based on the moment matching method borrowed from the theoretical computer science literature. Next, we consider two additional structures that one can exploit in the high-dimensional problem: effect homogeneity (Section 4) and prior knowledge of covariate distribution (Section 5). We propose novel estimators that can take advantage of these structures and enjoy faster convergence rates, making consistent estimation possible in the regime where conventional estimators examined in Section 3 fail to be consistent. Finally in Section 6 we perform numerical experiments to verify our results in previous sections. All the proofs and additional complementary results are presented in the appendix. To the best of our knowledge, our work is the first in the literature that directly analyzes the theoretical properties of different estimators of causal effects, provides a minimax lower bound for ATE estimation, and explores additional structures that we can exploit to achieve faster convergence rates, in the high-dimensional discrete setting.

2 Setup and Notation↩︎

In this section, we first introduce the data-generating process in the covariate setting and characterize the distributions of several counting statistics that play an important role in estimating causal effects. Identification assumptions of causal estimands and additional notation are further introduced with discussion.

2.1 Data Generating Process↩︎

Suppose we observe \(n\) i.i.d. copies of \(\mathbf{Z}=(X,A,Y)\) where \(X\) is the discrete covariate with \(d\) categories, \(A \in\{0,1\}\) is the binary treatment and \(Y \in\{0,1\}\) is the binary outcome. On the population level, the covariate \(X\) has a categorical distribution on \([d]=\{1,\dots,d\}\) with \[\mathbb{P}(X=k) = p_k, \,1 \leq k \leq d.\] Let \({\boldsymbol{p}}=(p_1,\dots,p_d)\) be the probability vector. In real applications one may observe multiple discrete covariates. For example, we may observe \(K\) binary covariates and it is easy to see that one could encode these binary variables as one categorical variable with \(d=2^K\) categories. Hence we will assume a single discrete covariate \(X\) in our problem. Given the covariate \(X=k\), the treatment \(A\) follows a Bernoulli distribution with parameter \(\pi_k\) \[A \mid X=k\, \sim \text{Bernoulli} (\pi_k) \text{ with } \pi_k =\mathbb{P}(A=1 \mid X=k).\] Conditioned on the covariate \(X=k\) and the treatment \(A=a\), \(Y\) has a Bernoulli distribution with parameter \(\mu_{ak}\) \[Y\mid A=a, X=k \, \sim \text{Bernoulli} (\mu_{a k}) \text{ with } \mu_{a k} =\mathbb{P}(Y=1 \mid X=k, A=a).\] \(\pi_k\) and \(\mu_{ak}\) are the propensity score and regression functions, respectively. We further denote \(q_{a k}=\mathbb{P}(X=k, A=a, Y=1)=p_k\pi_k^a\left(1-\pi_k\right)^{1-a}\mu_{a k}\) and \(w_k=\mathbb{P}(X=k, A=1)= p_k \pi_k\). Based on a sample of size \(n\), define the following empirical average estimators of the parameters: \[n \widehat{q}_{a k}= \, \sum_{i=1}^n I\left(X_i=k, A_i=a, Y_i=1\right) \sim \operatorname{Bin}\left(n, q_{a k}\right),\] \[n \widehat{w}_k =\, \sum_{i=1}^n I\left(X_i=k, A_i=1\right) \sim \operatorname{Bin}\left(n, w_k\right),\] \[n \widehat{p}_k =\, \sum_{i=1}^n I\left(X_i=k\right) \sim \operatorname{Bin}\left(n, p_k\right)\] for each \(1 \leq k \leq d\). The corresponding empirical average estimates of the propensity score \(\pi_k\) and regression function \(\mu_{ak}\) are \[\label{eq:mle-nuisance} \begin{align} \widehat{\pi}_{k} = &\, \frac{\widehat{w}_k}{\widehat{p}_k} = \frac{\#\{i: X_i=k, A_i=1\}}{\#\{i: X_i=k\}}, \\ \widehat{\mu}_{1k} = &\, \frac{\widehat{q}_{1k}}{\widehat{w}_k} = \frac{\#\{i: X_i=k, A_i=1, Y_i=1\}}{\#\{i: X_i=k,A_i=1\}}, \\ \widehat{\mu}_{0k} = &\, \frac{\widehat{q}_{0k}}{\widehat{p}_k-\widehat{w}_k} = \frac{\#\{i: X_i=k, A_i=0, Y_i=1\}}{\#\{i: X_i=k,A_i=0\}}, \end{align}\tag{1}\] where we define \(0/0=0\) whenever both the numerator and denominator are zero. This may happen when no individual in \(k\)-th category is assigned to treatment (\(\widehat{w}_k = 0\)) and hence the response \(Y\) under treatment is unavailable (\(\widehat{q}_{1k} = 0\)). Let \(\mathbf{X}^n = (X_1 ,\dots, X_n)\) and \(\mathbf{A}^n = (A_1, \dots, A_n)\) be the collection of covariates and treatments in the sample, respectively. According to the sampling schemes we have \[(n\widehat{p}_1, \dots,n\widehat{p}_d ) \, \sim \text{Multinomial} (n, p_1, \dots, p_d),\] \[n \widehat{w}_k=n \widehat{p}_k \widehat{\pi}_k \mid \mathbf{X}^n \sim \operatorname{Bin}\left(n \widehat{p}_k, \pi_k\right),\] \[\begin{align} &\,n \widehat{q}_{a k}= n \widehat{p}_k \left[a \widehat{\pi}_k+(1-a)\left(1-\widehat{\pi}_k\right)\right] \widehat{\mu}_{a k} \mid \mathbf{X}^n, \mathbf{A}^n \\ &\,\sim \operatorname{Bin}\left(n \widehat{p}_k \left[a \widehat{\pi}_k+(1-a)\left(1-\widehat{\pi}_k\right)\right] , \mu_{a k}\right). \end{align}\] Moreover, we have \(\widehat{w}_k \mathpalette{\independenT}{\perp}\widehat{w}_{\ell} \mid \mathbf{X}^n\) for \(k \neq \ell\) since conditioned on the number of samples in each category, the treatment assignment within different categories \(X=k,\ell (k \neq \ell)\) proceeds independently. Similarly, conditioned on the number of samples in each category and treatment assignment, the outcomes within different categories are conditionally independent, i.e., \(\widehat{q}_{a k} \mathpalette{\independenT}{\perp}\widehat{q}_{a \ell} \mid \mathbf{X}^n, \mathbf{A}^n\) for \(k \neq \ell\).

2.2 Causal Assumptions and Other Notation↩︎

To properly define and identify causal estimands of interests, we rely on the potential outcome framework [42], [43] and additional identification assumptions to connect counterfactual outcomes with the observed data. We use the random variable \(Y^a\) to denote the potential/counterfactual outcome we would have observed had a subject received treatment \(A = a\), which may be contrary to the observation \(Y\). The average treatment effect (ATE) \(\psi\) is defined as \[\psi = \mathbb{E}\left[Y^1-Y^0\right].\] The following identification assumptions are often imposed to identify \(\psi\) as a functional of the observed data distribution \(\mathbb{P}\).

Assumption 1. Consistency: \(Y = Y^a \text{\, if \,} A=a.\)

Assumption 2. Positivity: For any \(k \in [d]\), \(\pi_k \in [\epsilon, 1-\epsilon]\) for some constant \(\epsilon \in (0,1/2)\).

Assumption 3. No unmeasured confounding: \(Y^a \mathpalette{\independenT}{\perp}A \mid X\) for \(a=0,1\).

Under Assumptions 1–3, the ATE can be identified as \[\label{eq:ATE} \begin{align} \psi = &\, \mathbb{E}\left[\mathbb{E}(Y\mid X,A=1)-\mathbb{E}(Y\mid X,A=0)\right] \\ = &\, \sum_{k=1}^d p_k(\mu_{1k}- \mu_{0k}) =\sum_{k=1}^d p_k\left(\frac{q_{1k}}{w_k}- \frac{q_{0k}}{p_k-w_k}\right). \end{align}\tag{2}\] We refer the readers to [2] for detailed discussion on these identification assumptions. From this point forward, \(\psi\) will denote the functional of observed data distribution in 2 , which is equal to ATE when Assumptions 1–3 hold. If Assumption 1–3 are violated, the functional \(\psi\) should only be viewed as the expected difference in the regression functions between the treatment and control group, which may not represent a causal effect. Nonetheless, all our results still hold for \(\psi\) in 2 under only the positivity assumption.

In this paper we will consider the following model in which the number of categories for the covariate is at most \(d\): \[\label{eq:model-D} \mathcal{D}(\epsilon) = \left\{ \sum_{k=1}^d p_k = 1, 0 \leq p_k\leq 1 , \epsilon \leq \pi_k \leq 1-\epsilon, 0\leq \mu_{0k}, \mu_{1k} \leq 1 ,\, \forall k \in [d] \right \}.\tag{3}\] We will always assume positivity in addition to the basic bounds on model parameters. The binary nature of the outcome variable \(Y\) is assumed primarily for simplicity and is not essential to our analysis. Our results on rates remain valid as long as \(\mathbb{E}[Y^2 \mid X = k, A = a] \leq C\) for some constant \(C > 0\). This condition holds if, within each covariate level \(X = k\) and treatment arm \(A = a\), the outcome \(Y\) has bounded variance. For example, this is satisfied when \(Y\) is bounded or sub-Gaussian with a uniformly bounded sub-Gaussian norm given \(X=k\) and \(A=a\).

Remark 1. It is worth noting that in the high-dimensional regime where \(d\) can be large (the regime we are mainly interested in), the positivity assumption 2 can impose additional restriction on the observed distribution [44]. However, when positivity is violated, the ATE is not identified and may not be an appropriate estimand to focus on. This together with other possible issues due to violation of positivity are tangential to our main points and thus we may proceed assuming positivity holds. Some of our results characterize the rate in terms of \(\epsilon\) explicitly (e.g., Theorem 1) and still provide meaningful implications under weak positivity assumption where \(\epsilon = \epsilon_n\) shrinks to zero.

For a (possibly random) function \(f\) of the observation \(\mathbf{Z}= (X,A,Y)\), we use \(\mathbb{P}_n[f(\mathbf{Z})]\) to denote the sample average \(\frac{1}{n} \sum_{i=1}^n f(\mathbf{Z}_i)\), \(\mathbb{P}[f(\mathbf{Z})] = \int f(\mathbf{z}) d\mathbb{P}(\mathbf{z})\), and \(\|f\|_2\) to denote the \(L_2\)-norm \(\left[\int f^2(\mathbf{z}) d \mathbb{P}(\mathbf{z})\right]^{1/2}\), where all expectations are only taken with respect to the randomness of \(\mathbf{Z}\). For a bivariate function \(g\) on \((\mathbf{Z}_1,\mathbf{Z}_2)\), let \(\mathbb{U}_n[g(\mathbf{Z}_1,\mathbf{Z}_2)]\) denote the second-order U-statistic measure \(\frac{1}{n(n-1)} \sum_{i \neq j} g(\mathbf{Z}_i,\mathbf{Z}_j)\). For two sequences \(\{a_n\}\) and \(\{b_n\}\), we use \(a_n \lesssim b_n\) to denote that there exists a constant \(C>0\) such that \(a_n \leq C b_n\) when \(n\) is sufficiently large, and \(a_n \gtrsim b_n\) means \(b_n \lesssim a_n\). We use \(a\vee b\) and \(a \wedge b\) to denote the maximum and minimum of \(a\) and \(b\), respectively.

3 Average Treatment Effects↩︎

In this section, we focus on the theoretical properties of commonly used estimators of ATE \[\psi = \mathbb{E}\left[\mathbb{E}(Y\mid X,A=1)-\mathbb{E}(Y\mid X,A=0)\right] = \sum_{k=1}^d p_k(\mu_{1k}- \mu_{0k}) .\] A simple plug-in-style estimator based on empirical average estimates of model parameters \((p_k,\pi_k,\mu_{1k}, \mu_{0k})\) in 1 is \[\label{eq:plugin} \widehat{\psi} = \sum_{k=1}^d \widehat{p}_k \left(\widehat{\mu}_{1k}-\widehat{\mu}_{0k} \right)= \sum_{k=1}^d \widehat{p}_k\left(\frac{\widehat{q}_{1k}}{\widehat{w}_k}- \frac{\widehat{q}_{0k}}{\widehat{p}_k-\widehat{w}_k}\right),\tag{4}\] We also consider three popular estimators in the literature to estimate ATE, namely the outcome regression estimator \(\widehat{\psi}_{\text{reg}}\) [1], inverse probability weighting estimator \(\widehat{\psi}_{\text{ipw}}\) [3], [4] and doubly robust estimator \(\widehat{\psi}_{\text{dr}}\) [45], [46] defined as follows: \[\label{eq:ate-ests} \begin{align} \widehat{\psi}_{\text{reg}} = &\, \mathbb{P}_n\left[\widehat{\mu}_{1X} - \widehat{\mu}_{0X} \right], \\ \widehat{\psi}_{\text{ipw}} = &\, \mathbb{P}_n\left[\frac{AY}{\widehat{\pi}_X} - \frac{(1-A)Y}{1-\widehat{\pi}_X} \right], \\ \widehat{\psi}_{\text{dr}} = &\, \mathbb{P}_n\left[\frac{A(Y-\widehat{\mu}_{1X})}{\widehat{\pi}_X} + \widehat{\mu}_{1X} - \frac{(1-A)(Y-\widehat{\mu}_{0X})}{1-\widehat{\pi}_X}- \widehat{\mu}_{0X}\right], \end{align}\tag{5}\] where again in the IPW and DR estimator we define \(0/0=0\) whenever it occurs. In this paper, we will focus on the observational study setting where propensity score is unknown and needs to be estimated. In randomized experiments where the treatment process is known, the IPW estimator with known propensity scores is unbiased and \(\sqrt{n}\)-consistent.

Surprisingly, all these three commonly used estimators of ATE are numerically equivalent to the simple plug-in estimator 4 in the discrete covariate setting, if we use the same sample to estimate the nuisance parameters in 1 and take average over in 5 . This property unique to the discrete covariate setting is presented in regression coefficients estimation with missing covariates in [47]. Recently similar equivalence in ATE estimation is shown in [48] when the correct parametric model for propensity score is specified. The discrete covariate setting can be viewed as a special case where the design matrix only contains indicators specifying the category membership of each sample. We restate their result in the discrete covariate setting and emphasize that such equivalence still holds in the high-dimensional setting where some categories may not have treated/untreated samples observed, as long as we define \(0/0 = 0\) whenever it appears.

Proposition 1. Supposed \(\mathcal{D}=\{(X_i,A_i,Y_i), 1 \leq i \leq n\}\) is a sample of size \(n\) and \(X \in [d]\) is discrete. If the nuisance estimators are the empirical averages defined in equation 1 using \(\mathcal{D}\) and we take average over \(\mathcal{D}\) in 5 . Then we have \[\widehat{\psi}_{\text{reg}} = \widehat{\psi}_{\text{ipw}} = \widehat{\psi}_{\text{dr}} = \widehat{\psi}.\] So the regression, weighting and doubly robust estimators are numerically equivalent in the discrete covariate setting.

The numerical equivalence in Proposition 1 does not necessarily hold when the covariates have continuous components and smoothing is applied to construct estimates of propensity scores and regression functions. One explanation is that in the discrete case, the arguments used to show \(\mathbb{E}[{\mu}_{1X} ] = \mathbb{E}[AY/\pi_X]=\psi_1\) (See e.g., Chapter 2 in [2]) can be applied to the empirical distribution \(\mathbb{P}_n\) (i.e. replace every expectation and conditional expectation in the argument with sample averages) to show \(\mathbb{P}_n [\widehat{\mu}_{1X}] = \mathbb{P}_n[AY/\widehat{\pi}_X]=\widehat{\psi}\).

Proposition 1 shows that we only need to consider the properties of plug-in-style estimator \(\widehat{\psi}\) and all the results hold for estimators in 5 as well. In the following discussion, we will first derive the upper bound on the mean squared error of \(\widehat{\psi}\) in Section 3.1. Then we characterize the minimax lower bound in Section 3.2.

3.1 Upper Bound↩︎

In this section, we study the behavior of \(\widehat{\psi}\) in a potentially high-dimensional regime and characterize its mean squared error in terms of \((n,d)\). In the low-dimensional regime where \(d\) is fixed, the empirical average estimators of nuisance functions 1 have parametric convergence rates. One can expect this ideal property to propagate to plug-in-style estimator \(\widehat{\psi}\), which crucially depends on the empirical average estimates. We provide a central limit theorem of \(\widehat{\psi}\) in the appendix when \(d\) is fixed to avoid distracting the readers from the high-dimensional regime of main interest. In the following discussions, our main results hold non-asymptotically for all \(n,d\) sufficiently large. We will focus on the point estimation theory in this work and characterize the condition under which \(\widehat{\psi}\) is consistent. Let \(\psi_a = \mathbb{E}[Y^a] = \mathbb{E}[\mathbb{E}(Y\mid X,A=a)]= \sum_{k=1}^d p_k \mu_{ak}\) be the mean of potential outcome \(Y^a\) and \(\widehat{\psi}_a = \sum_{k=1}^d \widehat{p}_k \widehat{\mu}_{ak}\) be the corresponding plug-in-style estimator. We first summarize the exact bias of \(\widehat{\psi}_a\) in the following proposition.

Proposition 2. The exact bias of \(\widehat{\psi}_a\) is \[\begin{align} \mathbb{E}\left[\widehat{\psi}_1 - \psi_1\right] =&\, -\sum_{k=1}^d \mu_{1k}p_k(1-\pi_k)(1-p_k \pi_k)^{n-1}, \\ \mathbb{E}\left[\widehat{\psi}_0 - \psi_0\right] =&\, -\sum_{k=1}^d \mu_{0k}p_k\pi_k(1-p_k +p_k\pi_k)^{n-1}. \\ \end{align}\] Hence the exact bias of \(\widehat{\psi}=\widehat{\psi}_1 - \widehat{\psi}_0\) is \[\begin{align} &\,\mathbb{E}\left[\widehat{\psi} - \psi\right] \\ = &\,- \sum_{k=1}^d \left[ \mu_{1k}p_k(1-\pi_k)(1-p_k \pi_k)^{n-1} - \mu_{0k}p_k\pi_k(1-p_k +p_k\pi_k)^{n-1}\right] \end{align}\]

From the proof of Proposition 2, the bias of \(\widehat{\psi}_a\) comes from categories with no subjects receiving treatment \(A=a\), and hence it is not possible to obtain an unbiased estimator of regression functions \(\mu_{ak}\) within these categories. It is unclear from the exact bias how fast \(d\) can grow with \(n\) while still having vanishing bias. To make the dependency of bias on \((n,d)\) explicit, we maximize the bias over model \(\mathcal{D}(\epsilon)\) and characterize the worst-case bias in the next proposition. In the following discussion of this section, we will mainly consider the functional \(\psi_1\), with the understanding that analogous arguments can be applied to \(\psi_0\) to obtain the same rate/property.

Proposition 3. The worst-case bias of \(\widehat{\psi}_1\) is \[\sup_{\mathbb{P}\in \mathcal{D(\epsilon)}}\big |\mathbb{E}_\mathbb{P}[\widehat{\psi}_1 - \psi_1]\big| = \sup_{{\boldsymbol{p}}} \sum_{k=1}^d p_k (1-\epsilon)(1-\epsilon p_k )^{n-1}.\] As a consequence, we have the following bounds on the worst-case bias: \[\frac{1-\epsilon}{2e}\left(\frac{d-1}{\epsilon n}\wedge 1\right) \leq \sup_{\mathbb{P}\in \mathcal{D(\epsilon)}} \big|\mathbb{E}_\mathbb{P}[\widehat{\psi}_1 - \psi_1]\big| \leq \left( \frac{1-\epsilon}{\epsilon} \right) \frac{d}{n}.\]

Proposition 3 shows that a sufficient and necessary condition for the bias to vanish as \(n \rightarrow \infty\) in the worst case over model class \(\mathcal{D}(\epsilon)\) is \(d/n \rightarrow 0\). The following theorem further characterizes a bound on the variance of \(\widehat{\psi}_1\) and arrives at the MSE of the plug-in-style estimator.

Theorem 1. The variance of \(\widehat{\psi}_1\) is upper bounded as \[\sup_{\mathbb{P}\in \mathcal{D}(\epsilon)} \operatorname{Var}_{\mathbb{P}}(\widehat{\psi}_1) \leq \frac{C}{\epsilon n},\] where \(C>0\) is an absolute constant. Hence the MSE is upper bounded as \[\sup_{\mathbb{P}\in \mathcal{D}(\epsilon)} \mathbb{E}_\mathbb{P}\left[\left(\widehat{\psi}_1 - \psi_1\right)^2\right] \leq \frac{(1-\epsilon)^2}{\epsilon^2}\frac{d^2}{n^2} + \frac{C}{\epsilon n}.\]

Under positivity condition 2, the marginal probability for a subject to receive treatment is \(\mathbb{P}(A=1) \geq \epsilon\). Thus on average the total number of treated samples is lower bounded by \(\epsilon n\). Since the denominator of the bound on variance is exactly \(\epsilon n\), intuitively \(\epsilon n\) acts as the effective sample size for estimating \(\psi_1\). The results in Proposition 3 and Theorem 1 imply:

The plug-in estimator \(\widehat{\psi}_1\) is consistent in non-classical high-dimensional regimes where \(d\rightarrow \infty\) as long as \(d/n \rightarrow 0\) holds.
\(d/n \rightarrow 0\) is also a necessary condition for the worst-case bias to vanish as \(n \rightarrow \infty\). Thus consistency of \(\widehat{\psi}_1\) in the high-dimensional regime \(n \lesssim d\) is not achievable without further assumptions.

The intuition is that with \(n\) samples and \(d\) categories, there are \(n/d\) samples within each category on average. Without further assumptions, one needs to consistently estimate the regression functions \(\mu_{1k}\) in each category to achieve consistent estimation of \(\psi_1 = \sum_k p_k \mu_{1k}\). And consistent estimation of \(\mu_{1k}\) requires infinite samples assigned to each category, i.e. \(n/d \rightarrow \infty\).

Remark 2. When the covariate distribution is sparse in the sense that the support size \(s := |{k : p_k > 0}| \ll d\), we can replace \(d\) with \(s\) in the bound of Theorem 1, yielding a rate of \(s^2/n^2 + 1/n\). Thus, consistency is achieved as long as \(s/n \to 0\).

Remark 3. Note that in the estimators above (including \(\widehat{\psi}_{\text{reg}}, \widehat{\psi}_{\text{ipw}}, \widehat{\psi}_{\text{dr}}\), and \(\widehat{\psi}\)), the same sample is used to estimate the nuisance functions and to compute the final estimator. Alternatively, one can adopt a sample splitting strategy, as in the double machine learning literature [7], where a separate sample is used for nuisance estimation. In this case, the numerical equivalence in Proposition 1 no longer holds. However, sample splitting does not necessarily improve the convergence rate or guarantee consistency in the high-dimensional regime where \(n \lesssim d\). Specifically, assume the nuisance functions are estimated from an independent sample of size \(n\) using empirical averages, which are natural estimators in the discrete covariate setting without additional structural assumptions. Following the proof of Theorem 7, we have \[\|\widehat{\mu}_{aX}-\mu_{aX}\|_2=O_{\mathbb{P}}\left(\sqrt{\frac{d}{n}} \right), \, \|\widehat{\pi}_{X}-\pi_{X}\|_2=O_{\mathbb{P}}\left(\sqrt{\frac{d}{n}} \right),\] where for a function \(g_X\), we define \(\|g_X\|_2^2 = \sum_k p_k g_k^2\). The squared errors of the regression estimator \(\widehat{\psi}_{\text{reg}}\), the IPW estimator \(\widehat{\psi}_{\text{ipw}}\), and the doubly robust estimator \(\widehat{\psi}_{\text{dr}}\) are therefore dominated by \(\|\widehat{\mu}_{aX}-\mu_{aX}\|_2^2=O_{\mathbb{P}}(d/n)\), \(\|\widehat{\pi}_{X}-\pi_{X}\|_2^2=O_{\mathbb{P}}(d/n)\), and \(\|\widehat{\mu}_{aX}-\mu_{aX}\|_2^2\|\widehat{\pi}_{X}-\pi_{X}\|_2^2=O_{\mathbb{P}}(d^2/n^2)\), respectively. As a result, consistency of these estimators still requires \(d/n \to 0\), even when sample splitting is employed.

However, our results in this section do not preclude the existence of other consistent estimates. In order to conclusively determine whether consistent estimates exist in the high-dimensional regime, one needs to characterize the minimax lower bound for \(\psi_1\).

3.2 Minimax Lower Bound↩︎

In Section 3.1 we showed that \(d/n \rightarrow 0\) is a sufficient and necessary condition for the plug-in estimator \(\widehat{\psi}\) to be consistent. In this section, we study the existence of consistent estimators in the high-dimensional regime by considering the minimax lower bound for the mean of the regression function \(\psi_1 = \mathbb{E}[\mathbb{E}(Y\mid X, A=1)]\), with the understanding that similar arguments show \(\psi_0=\mathbb{E}[\mathbb{E}(Y\mid X, A=0)]\) and ATE share the same minimax rate as \(\psi_1\). Recall the model we consider is: \[\mathcal{D}(\epsilon) = \left\{ 0 \leq p_k\leq 1 , \epsilon \leq \pi_k \leq 1-\epsilon, 0\leq \mu_{0k}, \mu_{1k} \leq 1 ,\, \forall k \in [d] \right \},\] which corresponds to the setting where statisticians have knowledge on a set of possible values of \(X\) but some categories may have zero proportion.

By Theorem 1, the MSE of plug-in-style estimator satisfies \[\label{eq:plugin-p} \mathbb{E}_\mathbb{P}\left[\left(\widehat{\psi}_1 - \psi_1\right)^2\right] \lesssim \frac{d^2}{n^2} + \frac{1}{ n}, \forall \, \mathbb{P}\in \mathcal{D} (\epsilon),\tag{6}\] which shows a sufficient condition for the plug-in estimator to be consistent is \(d/n \rightarrow 0\), i.e. we require sample size \(n \gg d\) to achieve consistency. The following theorem characterizes the minimax lower bound of the estimation error of \(\psi_1\) over model class \(\mathcal{D} (\epsilon)\).

Theorem 2. In the regime \(d \lesssim n \log n\), we have \[R^*(d,n) := \inf_{\widehat{\psi}_1} \sup_{\mathbb{P}\in \mathcal{D} (\epsilon)} \mathbb{E}_{\mathbb{P}}\left[\left(\widehat{\psi}_1 - \psi_1\right)^2\right] \gtrsim \left( \frac{d}{n \log n}\right)^2 + \frac{1}{n},\] where the constant behind “\(\gtrsim\)" depends on \(\epsilon\).

The proof adopts the idea of moment matching commonly used in deriving minimax lower bounds for functionals of discrete distributions [37]–[39]. Specifically, we need to show that there exist two probability measures \(\mu_0, \mu_1\) over the tuple \((p, \pi, \mu )\) satisfying \[\begin{align} a_0 := \mathbb{E}_{\mu_0}[p] &= \mathbb{E}_{\mu_1}[p] \leq \frac{1}{d-1}, \tag{7} \\ \mathbb{E}_{\mu_0}\left[ p^i (p \pi)^j (p \pi \mu)^k \right] &= \mathbb{E}_{\mu_1} \left[ p^i (p \pi)^j (p \pi \mu)^k \right], \quad \forall i, j ,k \geq 0, \,i+j+k \leq K, \tag{8} \\ |\mathbb{E}_{\mu_0}[p \mu] - \mathbb{E}_{\mu_1}[p \mu]| &\gtrsim \frac{1}{n \log n},\tag{9} \end{align}\] with \(K \asymp \log n\). Once the measures \((\mu_0, \mu_1)\) are constructed, we construct two “difficult" hypotheses \((H_0, H_1)\), where under the hypothesis \(H_u\) with \(u\in \{0,1\}\), let \((p_k, \pi_k, \mu_{1k})\overset{\text{i.i.d.}}{\sim} \mu_u\) for \(1\leq k \leq d-1\), and \((p_{d}, \pi_{d}, \mu_{1d})=(1-(d-1)a_0, \epsilon, 0)\). On a high level, the mean constraint 7 ensures that \({\boldsymbol{p}}\) is an approximate probability measure with high probability, the moment matching constraint 8 ensures that \(H_0\) and \(H_1\) are statistically indistinguishable based on the observations \((\mathbf{X}^n, \mathbf{A}^n, \mathbf{Y}^n)\), and the separation condition 9 ensures that the value of \(\psi_1 = \sum_{k=1}^d p_k \mu_{1k}\) is separated apart by an amount of \(\Omega(d/(n\log n))\) under hypotheses \(H_0\) and \(H_1\). Based on these intuitions, the minimax lower bound in Theorem 2 follows from the method of fuzzy hypothesis [49].

To show the existence of measures \((\mu_0, \mu_1)\) satisfying 7 , 8 , and 9 , by the duality of moment matching and best polynomial approximation (cf. e.g. [50]), the maximum separation in 9 subject to the moment constraint 8 is characterized by \[\begin{align} \label{eq:3D-approx-error} \inf_{Q\in \mathbb{R}[x,y,z], \deg Q\le K} \sup_{(p,p\pi,p\pi\mu)\in D} |Q(p,p\pi,p\pi\mu) - p\mu|, \end{align}\tag{10}\] where the approximation domain \(D\) is a 3-dimensional polytope given by \[\begin{align} D = \{(x,y,z): 0\le x\le 1, \epsilon x\le y\le (1-\epsilon)x, 0\le z\le y \}. \end{align}\] Characterizing the best 3D polynomial approximation error 10 is generally very challenging [51], and we lower bound 10 via a proper one-dimensional subproblem, carefully chosen as \[\label{eq:1D-approx} \inf_{P \in \text{span}\{1/x, 1, x, \dots, x^{K}\} } \max_{x \in [c/K^2,1]} \left | \frac{x}{x+c/K^2}-P(x) \right|.\tag{11}\] The intuition behind the subproblem 11 is a one-dimensional trajectory \(x\mapsto (x,\epsilon x+c',\epsilon x)\in D\) (roughly) parallel to an edge of \(D\), which is further motivated by the approximation-theoretic results in [52] and [53]. Lower bounding 11 using machineary in [52] then resolves 8 and 9 . As for the mean constraint 7 , we overcome it by a change-of-measure trick, which is the reason behind the additional basis \(1/x\) and the constraint \(x\ge c/K^2\) in 11 . We defer the details to the appendix.

Theorem 2 shows in the regime \(d \gtrsim n \log n\), a consistent estimator of \(\psi_1\) does not exist over model class \(\mathcal{D} (\epsilon)\). Compared with the upper bound in 6 , we conclude that up to log-factors the plug-in estimator \(\widehat{\psi}_1\) is minimax optimal and one needs sample size at least of order \(d\) to consistently estimate \(\psi_1\) over model class \(\mathcal{D}(\epsilon)\). In applications, if the observed number of categories of the discrete covariate \(X\) is comparable to sample size \(n\), then we should interpret the estimated ATE carefully since it is possible that our estimator is not consistent in that high-dimensional regime. It is worth noting that the lower bound in Theorem 2 does not exactly match the upper bound provided by the plug-in-style estimator with the difference being a log-factor. It is possible that some estimator of \(\psi_1\) based on polynomial approximation, which further reduces the bias, enjoys the “effective sample size enlargement" property [37]–[39] and achieves the minimax lower bound. We leave the exploration of such estimators to future investigation.

In the following two sections, we consider how effect homogeneity and prior knowledge of the covariate distribution can allow faster rates and consistent estimation of causal effects in the regime \(n \lesssim d\) under certain scaling conditions on \(n\) and \(d\).

4 Role of Effect Homogeneity↩︎

In this section, we study the role of effect homogeneity in consistent estimation of treatment effects in the regime \(n \lesssim d\). Let \(\tau_k = \mathbb{E}[Y|X=k, A=1] - \mathbb{E}[Y|X=k, A=0] = \mu_{1k}-\mu_{0k}\) be the conditional average treatment effect (CATE) in the \(k\)-th category. In the high-dimensional regime where \(d=d_n\) can grow with \(n\), \(\tau_k\)’s should be viewed as a triangular array \(\{\tau_{nk}, 1\leq k\leq d_n, n\geq 1\}\) and we will slightly abuse the notation to denote \(\tau_k = \tau_{nk}\). Homogeneous effects (i.e. when \(\tau_k = \psi\) for all \(k \in [d]\)) can be helpful in terms of estimation. Intuitively, since the treatment effects are the same within each level of covariate, there is no need to consistently estimate CATE \(\tau_k\) for all possible \(d\) categories. One only needs to estimate the CATE within a few categories accurately and due to effect homogeneity, these estimated CATEs generalize to other levels of covariate, which potentially reduces the sample size required for consistent estimation. In this work, we adopt a novel form of approximate effect homogeneity via parameter capturing the extent of heterogeneity. Mathematically, denote \[\label{eq:homogeneity} \sigma_n := \max_{1 \leq k \leq d_n} |\tau_k -\psi|\tag{12}\] as the maximal effect heterogeneity. Note that \(\sigma_n=0\) corresponds to the constant conditional average treatment effects and \(\sigma_n=2\) imposes no restriction on the model class. Our parameterization can interpolate between these two extremes as \(\sigma_n\) varies.

4.1 Upper Bound↩︎

The estimator we propose to take advantage of the effect homogeneity is \[\label{eq:homo-est} \widehat{\tau} = \frac{\sum_{k=1}^d \widehat{t}_k \widehat{\tau}_k}{\sum_{k=1}^d \widehat{t}_k}\tag{13}\] where \(\widehat{t}_k = \widehat{p}_kI( 0 < \widehat{\pi}_k <1)\) is the indicator of whether both treated and untreated samples are observed in the \(k\)-th level and \(\widehat{\tau}_k = \widehat{\mu}_{1k}-\widehat{\mu}_{0k}\). The idea is to only estimate the CATE \(\tau_k\) within those categories with both treated and untreated units and take a weighted average over such categories. We restrict our attention to these categories so that unbiased estimation of \((\mu_{1k}, \mu_{0k})\) (and hence \(\tau_k\)) is possible. We again define \(0/0=0\) if \(\sum_{k=1}^d \widehat{t}_k = 0\), i.e. there is no such category that contains both treated and untreated units. Clearly on the event \(\left\{\sum_{k=1}^d \widehat{t}_k = 0\right\}\), we cannot obtain any information on treatment effects from \(\widehat{\tau}\). Under appropriate scaling conditions, the probability of this “adverse" event converges to \(0\) even in the case \(n \lesssim d\), as summarized in the following lemma.

Lemma 1. The chance of having every category consist of either all treated or all untreated units (i.e., no collisions of any treated and untreated units at any category) is upper-bounded as \[\mathbb{P}\left(\sum_{k=1}^d \widehat{t}_k = 0\right) \leq 2 \exp\left(-C (\epsilon) \frac{n^2}{n \vee d} \right)\] where \(C(\epsilon)> 0\) is a constant depending on \(\epsilon\).

This lemma is not only critical to our analysis of \(\widehat{\tau}\) but also of independent interests in the literature of applied probability. It is closely related to the birthday problem [54] and can be viewed as the occupancy problem [55], [56] under a different sampling scheme. In the classic occupancy problem, the number of treated and untreated units are fixed first and then they are assigned to different categories of \(X\) with probability vectors \({\boldsymbol{p}}^{\text{T}}\) and \({\boldsymbol{p}}^{\text{C}}\) separately. While in our setting people’s covariates are first sampled, following which their treatment assignments are determined. Importantly, Lemma 1 shows the probability that the denominator of \(\widehat{\tau}\) is \(0\) vanishes even in the high-dimensional regime \(n \lesssim d\), as long as \(d/n^2 \rightarrow 0\). With Lemma 1 we can derive the following theorem bounding the MSE of \(\widehat{\tau}\).

Theorem 3. The bias and variance of \(\widehat{\tau}\) are bounded as \[|\mathbb{E}[\widehat{\tau}-\psi]| \leq \sigma_n + 2 \exp\left( -C(\epsilon)\frac{n^2}{n \vee d} \right),\] \[\operatorname{Var}(\widehat{\tau}) \lesssim \sigma_n^2 + \frac{ d}{n^2} +\frac{1}{n},\] where \(C(\epsilon)\) and the constant behind “\(\lesssim\)" depend on \(\epsilon\). As a consequence, the estimator \(\widehat{\tau}\) is consistent if \(\sigma_n \rightarrow 0\) and \(d/n^2 \rightarrow 0\) in the regime \(n \lesssim d\).

Theorem 3 shows that if asymptotic effect homogeneity holds (i.e. \(\sigma_n \rightarrow 0\)), then the estimator \(\widehat{\tau}\) is consistent in the regime \(d/n^2 \rightarrow 0\). Compared with the rate condition \(d/n \rightarrow0\) required for consistency by the plug-in estimator \(\widehat{\psi}\), \(\widehat{\tau}\) has a faster convergence rate and enables us to achieve consistency in a wider regime under asymptotic effect homogeneity.

4.2 Minimax Lower Bound↩︎

In this section, we establish a matching minimax lower bound under exact effect homogeneity. Specifically, we consider the model class \[\mathcal{H}(\epsilon):=\left\{p_k=\frac{1}{d},\;\epsilon \le \pi_k \le 1-\epsilon,\;0 \le \mu_{0k},\mu_{1k}\le 1,\;\tau_k=\tau_{k'},\, \forall k,k'\in[d]\right\}.\] By Theorem 3, the estimator \(\widehat{\tau}\) attains the upper bound \[\frac{d}{n^2}+\frac{1}{n}.\] The following theorem shows that this rate is minimax optimal over \(\mathcal{H}(\epsilon)\).

Theorem 4. In the regime \(d \lesssim n^2\), we have \[\inf_{\widehat{\psi}} \sup_{\mathbb{P}\in \mathcal{H}(\epsilon)}\mathbb{E}_{\mathbb{P}}\left[\left(\widehat{\psi} - \psi\right)^2\right] \gtrsim \frac{d}{n^2} + \frac{1}{n},\] where the constant behind “\(\gtrsim\)" depends on \(\epsilon\).

Theorem 4 implies that consistency over the model class \(\mathcal{H}(\epsilon)\) cannot be achieved when \(d \gtrsim n^2\). Intuitively, in this regime, with high probability, very few categories contain repeated observations, and in particular very few contain both treated and untreated units. Moreover, one can construct pairs of distributions whose one-observation marginals within each category are identical, so that categories observed only once carry no information for distinguishing them. As a result, only within-category collisions are informative, and when \(d \gtrsim n^2\) there are too few such collisions to permit consistent estimation.

Theorem 3 establishes that \(\sigma_n \to 0\) is a sufficient condition for consistency when \(d/n^2 \to 0\). We regard \(\sigma_n \rightarrow 0\) as an interpretable structural assumption capturing asymptotic effect homogeneity. In addition, the lower bound in Theorem 2 suggests that some such restriction is necessary to extend consistency beyond the general \(d=o(n\log n)\) regime. We leave a sharp characterization of the minimax dependence on \(\sigma_n\) to future work.

5 Prior Knowledge of Covariate Distribution↩︎

In this section, we consider a different structure that could be exploited in causal effects estimation: the covariate distribution is known or can be estimated at a fast rate. The role of the covariate distribution has been studied in both causal functional estimation (e.g., expected conditional covariance (ECC); see [9]) and function estimation (e.g., CATE; see [15]), where incorporating knowledge of the covariate density can lead to improved convergence rates. When the covariate is discrete, faster rate of ATE estimation is not achievable with information on covariate distribution, as shown in Section 5.1. The fundamental difficulty is that the measure underlying ATE estimation is the product of covariate distribution and propensity score, whereas propensity score is unknown in our observational study setting. We thus switch attention to variance-weighted average treatment effects (WATE) in Section 5.2, which itself also receives an increasing amount of attention in the literature. The underlying measure of WATE is the covariate distribution and a faster rate of estimation is achieved with covariate information.

5.1 Average Treatment Effects↩︎

In this section, we introduce a minimax lower bound for \(\psi_1= \mathbb{E}[\mathbb{E}(Y|X,A=1)]\) under a uniform covariate distribution \({\boldsymbol{p}}= (1/d,\cdots,1/d)\). We show that a significant improvement on the rate of ATE estimation, such as the \(d/n^2\) rate discussed in Section 4, cannot be attained. This highlights the limitations of leveraging covariate distribution knowledge in improving the efficiency of ATE estimation. Consider the following model class where the covariate distribution is uniform over \([d]\): \[\mathcal{D}^U(\epsilon) = \left\{ p_k = 1/d , \epsilon \leq \pi_k \leq 1-\epsilon, 0\leq \mu_{0k}, \mu_{1k} \leq 1 ,\, \forall 1\leq k \leq d \right \}.\] The minimax lower bound of \(\psi_1\) over model class \(\mathcal{D}^U(\epsilon)\) is summarized in the following theorem.

Theorem 5. For any fixed constant \(\beta \in (0,1)\), in the regime \(n \lesssim d^{1-\beta}\) we have \[\inf_{\widehat{\psi}_1} \sup_{\mathbb{P}\in \mathcal{D}^U (\epsilon)} \mathbb{E}_{\mathbb{P}}\left[\left(\widehat{\psi}_1 - \psi_1\right)^2\right] \gtrsim 1,\] where the constant behind “\(\gtrsim\)" depends on \(\beta, \epsilon\).

Similar to the proof of Theorem 2, the proof of Theorem 5 relies on the method of fuzzy hypotheses [49] and the characterization of the following best 2D polynomial approximation error: \[\begin{align} \label{eq:2D-approx-error} \inf_{Q\in \mathbb{R}[x,y], \deg Q\le L} \max_{\pi\in [\epsilon, 1-\epsilon], \mu\in [0,1]} |Q(\pi, \pi\mu) - \mu|. \end{align}\tag{14}\] Again, we lower bound the above quantity via a suitable \(1\)-dimensional subproblem, which simply sets \(\mu = \epsilon/\pi\). The main difference from the proof of Theorem 2 is that, instead of choosing \(L \asymp \log n\), here we choose \(L = O(1)\) to be a large constant. It turns out that the approximation error 14 is lower bounded by a constant \(c=c(\epsilon,L)>0\), and the two fuzzy hypotheses are statistically indistinguishable as long as \(n\lesssim d^{1-\beta}\).

Theorem 5 establishes that the minimax rate for estimating \(\psi_1\) is lower bounded away from zero within the regime \(n\lesssim d^{1-\beta}\), provided the model class includes \(\mathcal{D}^U(\epsilon)\) as a subclass. This result indicates that enhancing the estimation rate of \(\psi_1\) by incorporating knowledge of the covariate distribution may not yield significant improvements when \(n \lesssim d\). Specifically, achieving consistency in the scenario where \(d/n^2 \rightarrow 0\)—analogous to the situation of effect homogeneity discussed in Section 4—is unattainable solely with insights into the covariate distribution.

The negative result in Theorem 5 can be explained by the fact that the underlying measure of ATE estimation is the product of covariate distribution and propensity score, as shown in analyzing its second-order estimator [14], [57]. In observational studies where the propensity score is unknown, it’s hard to obtain nearly unbiased estimates of \(p_k \mu_{1k}\) (the summand in \(\psi_1\)) with the observed data. One may write \(p_k \mu_{1k} = p_k \pi_k \mu_{1k}/\pi_k\) with \(p_k \pi_k \mu_{1k}\) unbiasedly estimable. The extra term \(1/\pi_k\) could be approximated by a polynomial with possibly diverging degrees \[\frac{1}{\pi_k} \approx \sum_{j=0}^L (1-\pi_k)^j,\] where \(L= L_n\) determines the approximation error. Since the product \((p_k \pi_k)^j\) can be estimated unbiasedly and \(p_k\) is known, this strategy could reduce the bias. However, this approach suffers from a large variance, particularly in situations where \(n \lesssim d\). In fact, our construction in the proof of Theorem 5 indeed replies on the difficulty of approximating \(1/\pi_k\) with a polynomial of \(\pi_k\), thereby highlighting the fundamental barrier in ATE estimation even with information on the covariate distribution.

5.2 Variance-weighted Average Treatment Effects↩︎

In this section we switch our attention to a popular variance-weighted average treatment effect estimand (WATE), defined as \[\theta= \frac{\mathbb{E}[\operatorname{Cov}(Y,A\mid X)]}{\mathbb{E}[\operatorname{Var}(A\mid X)]} = \frac{\mathbb{E}[\operatorname{Var}(A\mid X )\tau_X]}{\mathbb{E}[\operatorname{Var}(A\mid X)]},\] where the weight \(\operatorname{Var}(A\mid X)\) minimizes the asymptotic variance among all weighted average treatment effects with known weight under homoskedasticity [58]. Another interpretation of \(\theta\) in terms of robustness is that if we assume a partial linear (homogeneous effect) model for the outcome and use E-estimator [59] to estimate the parametric component, under model misspecification (i.e. the partial linear model is wrong and effect of treatment is heterogeneous), one can show the E-estimator converges to WATE [60]. The WATE has also been derived under different frameworks recently. [61] showed the marginal interventional effect under incremental propensity score intervention [62] coincides with WATE. The expected conditional covariance also appears in conditional independence test [63] and can be interpreted as a causal effect under a stochastic intervention [64]. In contrast to ATE, the underlying measure of \(\theta\) is the covariate distribution alone, making improvements from knowledge of covariate distribution possible.

The strategy of estimating \(\theta\) is to deal with the numerator and the denominator separately. Let \(\eta = \mathbb{E}[\operatorname{Cov}(Y,A\mid X)]\) and \(\rho = \mathbb{E}[\operatorname{Var}(A\mid X)]\). Further denote the regression function in \(k\)-th category as \(\mu_k : = \mathbb{E}[Y|X=k]\). The first-order and second-order influence functions of \(\eta\) under a nonparametric model, if they exist, have form \[\varphi_1(\mathbf{Z}) = (A-\pi_X)(Y-\mu_X),\] \[\varphi_2(\mathbf{Z}_1,\mathbf{Z}_2) = -\frac{(A_1-\pi_{X_1})I(X_1=X_2)(Y_2-\mu_{X_2})}{p_{X_1}}.\] Similar forms of influence functions of \(\rho\) can be obtained by replacing \((Y,\mu)\) with \((A,\pi)\). However, in general settings where covariates have continuous components, the second-order influence function does not exist for \(\eta\) [14]. Intuitively, for two independent samples \(X_1, X_2\) from a continuous distribution, we have \(X_1 \neq X_2\) with probability \(1\) and hence \(\varphi_2(\mathbf{Z}_1,\mathbf{Z}_2)=0\). So the “second-order influence function" is always \(0\) when covariates have continuous components and cannot help us improve the estimator. On the other hand, when \(X\) is discrete it is possible that \(X_1, X_2\) fall into the same category and \(\varphi_2(\mathbf{Z}_1,\mathbf{Z}_2) \neq 0\). Then the idea is to use second-order estimator based on \(\varphi_1, \varphi_2\) to simultaneously correct for the first-order and second-order bias of plug-in-style estimator \(\widehat{\eta} = \mathbb{P}_n[YA - \widehat{\mu}_X \widehat{\pi}_X]\). Consider the following second-order estimator of \(\eta\) and \(\rho\), \[\label{eq:high-order-eta} \widehat{\eta} = \mathbb{P}_n[(A-\widehat{\pi}_X)(Y-\widehat{\mu}_X)] - \mathbb{U}_n \left[ \frac{(A_1-\widehat{\pi}_{X_1})I(X_1=X_2)(Y_2-\widehat{\mu}_{X_2})}{\widehat{p}_{X_1}}\right],\tag{15}\] \[\label{eq:high-order-rho} \widehat{\rho} = \mathbb{P}_n\left[(A-\widehat{\pi}_X)^2\right] - \mathbb{U}_n \left[ \frac{(A_1-\widehat{\pi}_{X_1})I(X_1=X_2)(A_2-\widehat{\pi}_{X_2})}{\widehat{p}_{X_1}}\right],\tag{16}\] where \(\mathbb{P}_n,\mathbb{U}_n\) are empirical and U-statistic measures. The following theorem establishes the estimation guarantees of second-order estimators when the nuisance functions are estimated from a separate independent sample agnostically, i.e., we do not specify the way to estimate \(\mu\) and \(\pi\).

Theorem 6. Suppose the nuisance functions \(\pi_k,\mu_k\) and covariate distribution \(p_k\) are estimated from a separate independent sample \(D\) as \(\widehat{\pi}_k,\widehat{\mu}_k, \widehat{p}_k\) with \(\widehat{\mu}_k, \widehat{\pi}_k \in [0,1]\). Then for second-order estimators in 15 –16 we have \[\begin{align} |\mathbb{E}\left[\widehat{\eta}\mid D\right]-\eta| \leq&\, \|\widehat{\mu}-\mu\|_2 \|\widehat{\pi}-\pi\|_2 \max_{k} \left| 1-\frac{p_k}{\widehat{p}_k} \right|, \\ \operatorname{Var}\left(\widehat{\eta}\mid D\right) \lesssim &\, \frac{1}{n} \left(1+ \max_k \left| 1-\frac{p_k}{\widehat{p}_k} \right|\right)^2 + \frac{d}{n^2}\left(1+ \max_k \left| 1-\frac{p_k}{\widehat{p}_k} \right|\right)^2, \\ |\mathbb{E}\left[\widehat{\rho}\mid D \right]-\rho| \leq&\, \|\widehat{\pi}-\pi\|_2^2 \max_{k} \left| 1-\frac{p_k}{\widehat{p}_k} \right|, \\ \operatorname{Var}\left(\widehat{\rho}\mid D\right) \lesssim &\, \frac{1}{n} \left(1+ \max_k \left| 1-\frac{p_k}{\widehat{p}_k} \right|\right)^2 + \frac{d}{n^2}\left(1+ \max_k \left| 1-\frac{p_k}{\widehat{p}_k} \right|\right)^2, \end{align}\] where recall for a function \(f = f(\mathbf{Z})\) the \(L_2\)-norm is defined as \(\|f\|_2^2 = \int f^2(\mathbf{z}) d \mathbb{P}(\mathbf{z})\) and the constant behind “\(\lesssim\)" is an absolute constant.

Theorem 6 shows the advantage of second-order estimators: after correcting for the first-order and second-order bias, the conditional bias of \(\widehat{\eta}\) depends on the product of estimation error of \(\mu, \pi, p\) and is “third-order small". In the following theorem, we parameterize the estimation rate of covariate distribution to show the consistency of second-order estimators in the high-dimensional regime \(n \lesssim d\) with different choices of nuisance estimators \(\widehat{\pi}, \widehat{\mu}\).

Theorem 7. Suppose the nuisance functions \(\pi_k,\mu_k\) and covariate distribution \(p_k\) are estimated from a separate independent sample \(D\) as \(\widehat{\pi}_k,\widehat{\mu}_k, \widehat{p}_k\). Further assume the estimated covariate distribution \(\widehat{p}_k\) satisfies \[\max_{1 \leq k \leq d} \left | 1-\frac{p_k}{\widehat{p}_k} \right|\leq \xi_n < 1.\] If nuisance functions \(\pi_k, \mu_k\) are estimated from a sample of size \(n\) using empirical averages, then the second-order estimators of \(\eta\) and \(\rho\) in 15 –16 satisfy \[\mathbb{E}\left[ (\widehat{\eta}-\eta)^2 \right] \lesssim \xi_n^2 \frac{d\wedge n}{n} + \frac{d }{n^2}+ \frac{1}{n},\] \[\mathbb{E}\left[ (\widehat{\rho}-\rho)^2 \right] \lesssim \xi_n^2 \frac{d \wedge n}{n} + \frac{d}{n^2}+ \frac{1}{n}.\] If we set the nuisance functions \(\widehat{\pi}_k, \widehat{\mu}_k\) as \(0\), then the second-order estimators of \(\eta\) and \(\rho\) satisfy \[\mathbb{E}\left[ (\widehat{\eta}-\eta)^2 \right] \lesssim \xi_n^2 + \frac{d }{n^2}+ \frac{1}{n},\] \[\mathbb{E}\left[ (\widehat{\rho}-\rho)^2 \right] \lesssim \xi_n^2 + \frac{d}{n^2}+ \frac{1}{n}.\] where the constant behind “\(\lesssim\)" is an absolute constant.

As a consequence, under positivity \(\epsilon \leq \pi_k \leq 1-\epsilon\) and assume \(\widehat{\rho} \geq \epsilon(1-\epsilon)\), the estimator \(\widehat{\theta} = \widehat{\eta}/\widehat{\rho}\) satisfies \[\begin{align} \mathbb{E}\left[ (\widehat{\theta}-\theta)^2 \right] \lesssim &\, \mathbb{E}\left[ (\widehat{\eta}-\eta)^2 \right] + \mathbb{E}\left[ (\widehat{\rho}-\rho)^2 \right]\\ \lesssim &\, \xi_n^2 + \frac{d}{n^2}+ \frac{1}{n}, \end{align}\] where the constant behind “\(\lesssim\)" depends on \(\epsilon\).

We note that Theorem 7 implies a faster convergence rate than \(d^2/n^2 + 1/n\) when \(\xi_n\) is small, which is hard to achieve for estimators \(\widehat{p}_k\)’s constructed from a sample of size \(n\) in the high-dimensional regime \(n \lesssim d\). For example, in the uniform case \(p_k = 1/d, 1\leq k \leq d\), one can use Chernoff bound to show that with probability \(1-\delta\), \[\max_{1\leq k \leq d}|d\widehat{p}_k-1| \leq \sqrt{\frac{3d \log(2d/\delta)}{n}}.\] Thus we need a sample size of order \(d \log d\) to guarantee small uniform estimation errors of \(p_k\), which is not achievable in the case \(n \lesssim d\). The improvements on the convergence rate of \(\widehat{\eta}\) should come from prior knowledge/assumptions on covariate distribution. The results imply that consistency is still possible in the regime \(n \lesssim d\) when such information on covariate distribution is available. One example is that the covariate distribution is known to statisticians, then we have \(\widehat{p}_k = p_k,\xi_n = 0\) and the estimators are consistent as long as \(d/n^2 \rightarrow 0\). Another useful setting is semi-supervised causal inference [65]–[68]: there is a large number (\(\gg n\)) of individuals in the database who were not selected into the randomized trial/observational study, hence treatment and outcome are not available for them but one can use the database to estimate the covariate distribution very accurately thanks to the large sample size. Then \(\xi_n\) is small and the estimation error of \(\widehat{\eta}\) will be small according to Theorem 7.

When the covariate distribution can be estimated very well and \(\xi_n\) is small, we can use arbitrary (bounded) nuisance estimators and still maintain consistency. As shown in Theorem 7, simply setting \(\widehat{\pi}, \widehat{\mu}\) as \(0\) leads to consistent estimation as long as \(\xi_n\) vanishes and \(d/n^2 \rightarrow 0\). Using empirical average estimators has advantage over setting \(\widehat{\pi}, \widehat{\mu}\) as \(0\) only when \(d \ll n\). The interpretation is that when \(\xi_n\) is small, the conditional bias in Theorem 6 is negligible regardless of the convergence rate of \(\widehat{\pi}, \widehat{\mu}\). The conditional variance, whose order is independent of the nuisance estimation rate, is the dominant term in the MSE.

We conclude this section with additional discussion on second-order estimation of \(\psi_1 = \mathbb{E}[\mathbb{E}(Y\mid X,A=1)]\) in the discrete covariate setting. The forms of first-order and second-order influence functions are \[\phi_1(\mathbf{Z}) = \frac{A(Y-\mu_{1X})}{\pi_X}+ \mu_{1X},\] \[\phi_2(\mathbf{Z}_1,\mathbf{Z}_2) = -\left(\frac{A_1}{\pi_{X_1}}-1 \right)\frac{I(X_1=X_2)}{p_{X_1}\pi_{X_2}}A_2(Y_2-\mu_{1X_2}).\] The second-order estimator is then \[\label{eq:high-order-ate} \mathbb{P}_n\left[\widehat{\phi}_1(\mathbf{Z})\right] + \mathbb{U}_n\left[\widehat{\phi}_2(\mathbf{Z}_1,\mathbf{Z}_2)\right].\tag{17}\] Assume we use a separate independent sample \(D\) to estimate nuisance functions \(\pi, \mu_1\) and covariate distribution \(p\), the conditional bias of second-order estimator 17 is \[\mathbb{E}\left[ (\widehat{\mu}_{1X}-\mu_{1X})\left(1-\frac{\pi_X}{\widehat{\pi}_X} \right)\left(1-\frac{\pi_X p_X}{\widehat{\pi}_X \widehat{p}_X} \right) \mid D\right].\] In order to obtain a faster rate by making use of information on covariate distribution, one also needs an accurate estimator of propensity score \(\pi\) to ensure \[\label{eq:product-pi-p} \max_{1 \leq k \leq d}\bigg |1-\frac{\pi_k p_k}{\widehat{\pi}_k \widehat{p}_k} \bigg |\tag{18}\] is small. In observational studies where the treatment process is unknown, it is hard to estimate the product of propensity score and covariate distribution well and make 18 small in the regime \(n \lesssim d\), due to the same reason discussed after Theorem 7. This further demonstrates the difficulty in ATE estimation with solely information on covariate distribution.

The fact that covariate distribution improves the estimation rates for ECC and WATE, but not for ATE, can be understood through the framework of [9]. In their analysis, the improved convergence rates achieved by higher-order estimators rely on knowledge of the so-called “\(g\)-function” (defined in their equation (3.3)), which is assumed to be known or sufficiently smooth. For the expected conditional covariance (ECC) \(\eta\), this \(g\)-function simplifies to the covariate density. Hence, knowledge of the covariate distribution directly contributes to improved estimation rates. In contrast, for the covariate-adjusted mean in the treated group \(\psi_1\), the \(g\)-function involves the product of the covariate density and the propensity score. Therefore, faster rates can only be achieved when information is available about this product, not the covariate distribution alone. In practice, WATE is also a useful and interpretable summary measure of treatment effects that is often important to estimate and report.

6 Simulation Study↩︎

In this section, we provide numerical results to verify the theoretical properties of the estimators discussed in Section 3–5. Consider the following data generating process: \(X\) has a categorical distribution with uniform probability \({\boldsymbol{p}}=(1/d,\dots, 1/d)^{\top} \in \mathbb{R}^d\). For each category \(X=k\), the propensity score \(\pi_k = 1/2\) and the conditional means of outcome within treated and untreated groups are \(\mu_{1k}=1/2, \mu_{0k}=1/4\), respectively. Hence the treatment effects are homogeneous across different levels (i.e. \(\tau_k = 1/4\) for all \(k \in [d]\)) in our simulation setting and expected conditional covariance between treatment and outcome \(\eta=1/16\) (considered in Section 5). For each estimator we generate data containing \(n\) observations with \(n \in \{100,500,1000,5000,10000 \}\), apply the proposed method, and compute the estimated RMSE. Such procedure is repeated \(M=500\) times and the average RMSE is reported in each plot.

We first consider the plug-in estimator discussed in Section 3. For each sample size \(n\) let the number of categories \(d=\lfloor n^{\gamma} \rfloor\) with \(\gamma \in \{0.5, 0.55, 0.6,\\ \dots, 1.5\}\). The idea is to evaluate how large the order of \(d\) can be to maintain consistency [69]. The average treatment effect is \(\psi=1/4\) and we estimate RMSE as \[\widehat{\text{RMSE}} =\left[ \frac{1}{M}\sum_{m=1}^M (\widehat{\psi}^m-\psi)^2 \right]^{1/2},\] where \(\widehat{\psi}^m\) is the estimated ATE from \(m\)-th repetition. The relationship between RMSE and \(\gamma\) is summarized in Figure 1. For a fixed sample size \(n\), as \(\gamma\) increases (i.e. \(d\) increases) the estimated RMSE also increases as expected. We see a clear phase transition in the plot: in the region \(\gamma < 0.8\), the RMSE is quite stable of order \(1/\sqrt{n}\); around \(\gamma=1\) the RMSE increases drastically, indicating the plug-in estimator starts to have larger error. This corresponds to our theoretical results in Section 3: the plug-in estimator is consistent if and only if \(d/n \rightarrow 0\).

Figure 1: RMSE V.S. \gamma for the plug-in estimator, where \gamma controls the order of d. — Figure 1: RMSE V.S. \(\gamma\) for the plug-in estimator, where \(\gamma\) controls the order of \(d\).

Then we consider the estimator 13 proposed in Section 4 under effect homogeneity. Again for each fixed sample size \(n\) we let \(d=\lfloor n^{\gamma} \rfloor\) with \(\gamma \in \{0.5, 0.55, 0.6, \dots, 2\}\). Note that the treatment effects are homogeneous in our setting (i.e. \(\sigma_n = 0\)) and we expect \(\widehat{\tau}\) in equation 13 to be consistent in a wider regime.

The relationship between RMSE and \(\gamma\) is summarized in Figure 2. The phase transition occurs in the region \(\gamma > 1\) (instead of at \(\gamma=1\)) and the estimator in 13 has a smaller error in the region \(1 < \gamma < 1.25\) compared with the plug-in estimator when \(n \geq 1000\). For instance, in the case \(n=1000, \gamma = 1\) the plug-in estimator has RMSE around \(0.08\) while the estimator under effect homogeneity has RMSE around \(0.05\). This coincides with our theoretical results in Section 4 that the estimator in equation 13 has a faster rate and is consistent in a wider regime.

Figure 2: RMSE V.S. \gamma for the estimator under effect homogeneity, where \gamma controls the order of d. — Figure 2: RMSE V.S. \(\gamma\) for the estimator under effect homogeneity, where \(\gamma\) controls the order of \(d\).

To further understand the order of the RMSE, for \(n \in \{1000, 10000\}\) we include the theoretical upper bound on RMSE, which is of order \(C_1\sqrt{d/n^2} = C_1 n^{\gamma/2-1}\), in the plot as a benchmark. Here the constant \(C\) is chosen as \(1.5\) to fit the empirical RMSE curve. The results are summarized in Figure 3. When \(\gamma < 1.5\), the empirical RMSE fits the theoretical upper bound quite well. When \(\gamma>1.5\) the empirical RMSE starts to deviate from the theoretical bound. In our experiments we found when \(\gamma > 1.5\) and \(d\) is large, the denominator in the estimator 13 is usually small (i.e., only a few categories have both treated and untreated samples) and the variation is large, which may explain the deviation of the empirical RMSE from the theoretical bound.

Figure 3: Comparison of theoretical and empirical order of RMSE

Finally, we evaluate the performance of the second-order estimator in 15 of expected conditional covariance \(\eta\) (the estimator 16 of expected conditional variance of treatment is expected to have similar performance). In our setting \(\eta=1/16\) and for each fixed sample size \(n\), let \(d=\lfloor n^{\gamma} \rfloor\) with \(\gamma \in \{0.5, 0.55, 0.6, \dots, 2\}\). We set the estimates of probabilities \(\widehat{p}_k\)’s as the true \(p_k = 1/d\) and hence \(\xi_n = 0\) in Theorem 7. The estimates \(\widehat{\pi}_k, \widehat{\mu}_k\)’s are all set to \(0\). The results are summarized in Figure 4. Similar to Figure 2, the estimated RMSE is quite stable in the region \(\gamma < 1.25\) and the phase transition seems to happen in the region \(\gamma > 1\).

Figure 4: RMSE V.S. \gamma for the second-order estimator 15 using true covariate distribution {\boldsymbol{p}}, where \gamma controls the order of d. — Figure 4: RMSE V.S. \(\gamma\) for the second-order estimator 15 using true covariate distribution \({\boldsymbol{p}}\), where \(\gamma\) controls the order of \(d\).

We also plot the relationship between empirical RMSE and theoretical bound \(C_2\sqrt{d/n^2} =C_2 n^{\gamma/2-1}\) with \(C_2 = 0.5\) in Figure 5. We see the empirical RMSE matches the theoretical bound very well.

Figure 5: Comparison of theoretical and empirical order of RMSE

7 Real Data Analysis↩︎

In this section, we employ the methods presented in Sections 3 and 4 to investigate the impact of 401(k) participation on net financial assets, utilizing data from the study by [70]. Since 401(k) eligibility is not randomly assigned, the effects must be assessed within an observational study framework. Following the strategy outlined in [7], [71], [72], 401(k) eligibility can be treated as exogenous after adjusting for confounding variables related to job choice that might be correlated with eligibility, such as income and education level.

The outcome variable \(Y\) represents net financial assets, which includes IRA balances, 401(k) balances, checking accounts, U.S. savings bonds, other interest-earning accounts in banks and financial institutions, other interest-earning assets (such as personally held bonds), stocks, and mutual funds, less nonmortgage debt. The treatment variable \(A\) is an indicator of 401(k) eligibility, where \(A=1\) represents eligibility. The confounding variables \(X\) include age, income, family size, years of education, marital status, two-earner status, defined benefit pension status, IRA participation, and home ownership. In our analysis, the variables age, income, family size, and years of education are discretized into categorical variables with four levels each. The sample size is \(n=9915\), and the number of categories, considering a saturated model with all possible interactions, is \(d=8196\).

The treatment effect estimated using the plug-in-style estimator discussed in Section 3 is \(7802\). The 95% confidence intervals, based on bootstrap and asymptotic normality in the appendix, are \([5352, 9559]\) and \([6212, 10238]\), respectively. It is important to note that, given the sample size \(n=9915\) and number of categories \(d = 8196\), these confidence intervals may not achieve the nominal coverage probability [73]; thus, they are presented for reference. Simulation results illustrating bootstrap-based variance estimation in high-dimensional settings with discrete covariates are provided in the appendix. These results indicate that 401(k) eligibility significantly increases the net financial assets. This finding aligns with the results of [7], [71], [72]. If researchers are confident that the covariates account for all relevant confounders, our estimate can be interpreted causally; otherwise, it should be understood as the expected difference between the covariate-adjusted regression functions for the treatment and control groups.

Additionally, we identified 993 categories in our data where positivity may be violated, meaning there is no overlap between treatment and control within those categories. The estimator in 13 , which only considers categories with both treated and control units, yields an estimate of \(9149\), with a 95% bootstrap confidence interval of \([6164,11715]\). This confidence interval is comparable in length to that of the plug-in-style estimator. It is possible that the treatment effects are heterogeneous, suggesting that the faster convergence rate of the estimator in 13 under effect homogeneity may not be fully achievable with the given data.

8 Discussion↩︎

In this paper, we studied the treatment effects estimation problem in the context of high-dimensional discrete covariates. Theoretical properties of commonly used regression, weighting and doubly robust estimators are examined in this non-classic regime. We also evaluated the role of effect homogeneity and covariate distribution in treatment effects estimation and proposed estimators that can properly take advantage of these structures and achieve faster convergence rates. Finally, we explored the fundamental limits of treatment effects estimation and showed consistent estimation of ATE is a difficult task on high-dimensional data. The discrete covariate setting is not only an interesting base case but also informative for the general dataset with continuous components. We hope our work can help researchers appropriately understand and interpret the treatment effects estimated from datasets with many covariates.

There are several possible extensions for future work. In this paper, we borrowed the moment matching tools from the theoretical computer science literature [37]–[39] to study the minimax lower bound on ATE estimation. It would be interesting to explore how to construct estimators of ATE based on polynomial approximation theory and realize the “effective sample size enlargement" phenomenon, which is feasible in entropy and support size estimation. Moreover, the minimax lower bound for ATE we proved is an initial result that does not exactly match the upper bound. More effort is needed to come up with new construction and tighten the lower bound if possible. Examining how other structural assumptions can help us achieve faster estimation rates is also an interesting topic. For example, although the covariate may take many categories, the conditional means \(\mu_{ak}\) may be constant over a small number of (known or unknown) groups: \[\mu_{ak} = \mu_{a,g(k)}, \qquad g(k) \in [G], \qquad G \ll d.\] When the grouping map \(g(\cdot)\) is known, the problem reduces from estimating \(2d\) stratum-specific means to estimating only \(2G\) group-level means, and Theorem 1 then yields the corresponding rate \(G^2/n^2 + 1/n\) (with \(d\) replaced by \(G\)). When the grouping is unknown, one can encourage a small number of distinct values among the \(\mu_{ak}\) by solving \[\min_{\mu_{a1},\dots,\mu_{ad}} \sum_{k=1}^d \bigl(\widehat{\mu}_{ak} - \mu_{ak} \bigr)^2 \;+\; \lambda \sum_{k < \ell} \bigl|\mu_{ak} - \mu_{a\ell}\bigr|,\] motivated by the literature on convex clustering with fusion penalties [74]–[78]. The fusion term shrinks stratum-specific means toward one another and can therefore recover, or approximately recover, latent grouping structure when many of the \(\mu_{ak}\) are equal or nearly equal. The resulting estimates of \(\mu_{ak}\) can then be substituted into the plug-in ATE estimator considered in Section 3. Similarly, sparsity may also arise in the effect modifiers, in the sense that the category-specific treatment effects \[\tau_k = \mu_{1k} - \mu_{0k}\] are constant across a small number of groups. In this case, one can regularize the \(\tau_k\) directly by solving \[\min_{\tau_1,\dots,\tau_d} \sum_{k=1}^d \bigl(\widehat{\tau}_{k} - \tau_{k} \bigr)^2 \;+\; \lambda \sum_{k < \ell} \bigl|\tau_{k} - \tau_{\ell}\bigr|,\] which encourages a small number of distinct values among the \(\tau_k\).

Other extensions could involve the estimation of different causal estimands in a similar discrete setting, including time-varying treatment effects, generalizability and transportability, optimal treatment regimes, instrumental variable and more. Developing a comprehensive understanding on the properties of popular estimators and fundamental limits of different causal functionals in the discrete setting could be an interesting avenue left for future investigation.

Appendix

9 Asymptotic Normality of the Plug-in Estimator↩︎

In this section, we provide a central limit theorem for the plug-in estimator when the number of categories \(d\) is fixed as a constant, summarized in the following theorem.

Theorem 8. Supposed \(X \in [d]\) is discrete with \(d\) fixed and the nuisance estimators are the empirical averages defined in 1 . Then we have \[\sqrt{n}\left(\widehat{\psi}-\psi\right) \stackrel{d}{\rightarrow} N\left(0, \operatorname{Var}(\varphi(\mathbf{Z}))\right),\] where \[\varphi(\mathbf{Z}) = \mu_{1X}-\mu_{0X}+\left (\frac{A}{\pi_X}-\frac{1-A}{1-\pi_X}\right )\left (Y-\mu_{AX}\right )\] is the first-order influence function of ATE \(\psi\) under a nonparametric model.

Theorem 8 implies that when the covariate is discrete with fixed dimension \(d\), the plug-in-style estimator \(\widehat{\psi}\) is \(\sqrt{n}\)-consistent and asymptotically normal. Hence in low-dimensional problems, the plug-in estimator \(\widehat{\psi}\) enjoys appealing properties and we can construct confidence intervals and perform statistical tests on ATE based on Theorem 8. It is worth noting that asymptotic normality also holds in the regime \(d=o(\sqrt{n})\). However, truncating the propensity score estimates at \(\epsilon\) and \(1-\epsilon\) is required to avoid instability induced by imprecise estimation of \(\pi_k\) in the high-dimensional regime. Moreover, the empirical influence function \(\widehat{\varphi}\) belongs to a Donsker class when \(d\) is fixed since it could be expressed as finite-dimensional parametric models (the dimension depends on \(d\)). As \(d=d_n\) grows with \(n\), they may not belong to a Donsker class and sample splitting is required to control the empirical process term.

10 Additional Simulation Results↩︎

In this section, we present additional simulation results on the asymptotic normality of the estimators studied in the main text, as well as on the accuracy of bootstrap-based variance estimation.

10.1 Asymptotic Normality↩︎

We first illustrate the asymptotic normality of the estimators considered in the paper, including the plug-in estimator \(\widehat{\psi}_1\) from Section 3, the estimator \(\widehat{\tau}\) exploiting effect homogeneity from Section 4, and the second-order estimator \(\widehat{\eta}\) from Section 5. The data-generating process is the same as in Section 6. For each repetition \(m \in [M]\), we compute a standardized version of the estimator. For example, the standardized plug-in estimator is \((\widehat{\psi}_1^{(m)}-\psi_1)/\widehat{\mathrm{sd}}(\widehat{\psi}_1)\), where \(\widehat{\mathrm{sd}}(\widehat{\psi}_1)\) denotes the empirical standard deviation of \(\widehat{\psi}_1\) across the \(M\) repetitions. Throughout, we center the estimators at the true parameter values, rather than at their Monte Carlo means. We then plot the estimated densities of the standardized estimators for different choices of \(\gamma\), recalling that \(d=\lfloor n^\gamma\rfloor\). We focus on \(n=1000\) and \(n=10{,}000\). The results are summarized in Figures 6–8.

Figure 6: Density of standardized plug-in estimator \widehat{\psi}_1. The black dashed line is the standard normal density. — Figure 6: Density of standardized plug-in estimator \(\widehat{\psi}_1\). The black dashed line is the standard normal density.

Figure 7: Density of the estimator \widehat{\tau} exploiting effect homogeneity. The black dashed line is the standard normal density. — Figure 7: Density of the estimator \(\widehat{\tau}\) exploiting effect homogeneity. The black dashed line is the standard normal density.

Figure 8: Density of the second-order estimator \widehat{\eta}. The black dashed line is the standard normal density. — Figure 8: Density of the second-order estimator \(\widehat{\eta}\). The black dashed line is the standard normal density.

In Figure 6, when \(\gamma \in \{0.3, 0.5, 0.7\}\), the estimated densities of the standardized plug-in estimator closely match the standard normal distribution, consistent with the asymptotic normality established in Theorem 8. When \(\gamma=0.9\), the density remains qualitatively bell-shaped, but it is no longer centered at zero because we center the estimator using the true parameter \(\psi_1\) rather than its Monte Carlo mean. In this high-dimensional regime the plug-in estimator is biased, so centering at \(\psi_1\) induces a visible shift in the standardized distribution.

In Figures 7–8, when \(\gamma \in \{0.4, 0.8, 1.2\}\), the estimated densities of the standardized estimators \(\widehat{\tau}\) and \(\widehat{\eta}\) closely match the standard normal distribution. Since these estimators can attain faster convergence rates under the additional structure they exploit, the standardized distributions remain approximately normal even at \(\gamma=1.2\). This contrasts with Figure 6, where the plug-in estimator \(\widehat{\psi}_1\) is biased in the high-dimensional regime and its standardized distribution is noticeably shifted when \(\gamma=0.9\). When \(\gamma=1.6\), the densities deviate from normality for \(n=1000\) but move closer to a standard normal shape for \(n=10{,}000\), suggesting that a CLT may still hold in this more challenging regime. We leave a formal analysis of asymptotic normality for \(\widehat{\tau}\) and \(\widehat{\eta}\) as an interesting direction for future work.

10.2 Bootstrap Variance Estimation↩︎

We then examine the accuracy of bootstrap-based variance estimation in our high-dimensional setting with discrete covariates. We focus on estimating the variance of the plug-in estimator \(\widehat{\psi}_1\). The data-generating process is the same as in Section 6. For each replication \(m \in [M]\), we estimate the variance of \(\widehat{\psi}_1^{(m)}\) using the nonparametric bootstrap with \(B=500\) resamples. We then compare the mean and median of the resulting bootstrap variance estimates to the empirical variance of \(\widehat{\psi}_1^{(m)}\) across \(M=500\) replications, for each \(\gamma \in \{0.5,0.6,\ldots,1.5\}\). The results are summarized in Figure 9.

Figure 9: Comparison between estimated variance from bootstrap with the empirical variance

In Figure 9, the mean and median of the bootstrap variance estimates are generally smaller than the empirical variance of \(\widehat{\psi}_1\) across replications, which serves as a consistent estimate of the true sampling variance. This pattern holds for both \(n=1000\) and \(n=10{,}000\). Overall, these results suggest that, in our simulated high-dimensional discrete-covariate setting, the nonparametric bootstrap can systematically underestimate the variance of \(\widehat{\psi}_1\). When \(d\) is large relative to \(n\), many categories have \(0\) or very small counts. The sampling distribution of plug-in estimators is then strongly influenced by whether a stratum appears at all and whether it has both treated/control observations. A bootstrap resample is drawn from the observed data, so it cannot recreate “new” strata that were absent in the original sample and tends to under-represent the variability coming from these appearance/disappearance events. Developing more robust variance estimators in this regime is an interesting direction for future work.

11 Proof of Main Results↩︎

11.1 Proof of Proposition 1 ↩︎

Proof. We will prove three estimators for \(\psi_1 = \mathbb{E}[Y^1] = \sum_{k=1}^d {p}_k {\mu}_{1k}\) are equal and similar results hold for \(\mathbb{E}[Y^0] = \sum_{k=1}^d {p}_k {\mu}_{0k}\). First consider the regression estimator \[\begin{align} \widehat{\psi}_{1,\text{reg}} = &\, \frac{1}{n} \sum_{i=1}^n \widehat{\mu}_{1X_i} \\ =&\, \frac{1}{n} \sum_{i=1}^n \sum_{k=1}^d I(X_i = k) \widehat{\mu}_{1X_i} \\ =&\, \frac{1}{n} \sum_{i=1}^n \sum_{k=1}^d I(X_i = k) \widehat{\mu}_{1 k} \\ = &\, \sum_{k=1}^d \left(\frac{1}{n} \sum_{i=1}^n I(X_i = k) \right) \widehat{\mu}_{1 k} \\ = &\, \sum_{k=1}^d \widehat{p}_k\widehat{\mu}_{1 k} = \widehat{\psi}, \end{align}\] where the second equation follows from the fact \(\sum_{k=1}^d I(X_i=k)=1\). In fact, given \(n\) samples, we have \(\sum_{k=1}^d I(X_i=k)=\sum_{k: \widehat{p}_k >0} I(X_i=k)=1\). With this in mind, for the inverse probability weighting estimator we have \[\begin{align} \widehat{\psi}_{1,\text{ipw}} = &\, \frac{1}{n}\sum_{i=1}^n \frac{A_i Y_i}{\widehat{\pi}_{X_i}} \\ = &\, \frac{1}{n}\sum_{i=1}^n \sum_{k:\widehat{p}_k>0} I(X_i=k)\frac{A_i Y_i}{\widehat{\pi}_{X_i}} \\ = &\, \frac{1}{n}\sum_{i=1}^n \sum_{k:\widehat{p}_k>0} I(X_i=k) \frac{A_i Y_i}{\widehat{\pi}_{k}} \\ = &\, \sum_{k:\widehat{p}_k>0} \frac{1}{n}\sum_{i=1}^n \frac{I(X_i=k)A_i Y_i}{\widehat{\pi}_{k}} \\ = &\, \sum_{k:\widehat{p}_k>0} \widehat{p}_k \frac{1}{n}\sum_{i=1}^n \frac{I(X_i=k)A_i Y_i}{\widehat{p}_k\widehat{\pi}_{k}} \\ = &\, \sum_{k:\widehat{p}_k>0} \widehat{p}_k \frac{\widehat{q}_{1k}}{\widehat{w}_k} \\ = &\, \sum_{k:\widehat{p}_k>0} \widehat{p}_k \widehat{\mu}_{1k} \\ = &\, \sum_{k=1}^d \widehat{p}_k \widehat{\mu}_{1k} = \widehat{\psi}, \end{align}\] where \(\widehat{q}_{1k} = \frac{1}{n}\sum_{i=1}^n I(X_i=k)A_i Y_i\) and \(\widehat{q}_{1k}/\widehat{w}_k = \widehat{\mu}_{1k}\) by definition in 1 . Note that all the equations still hold when some categories have no treated samples (i.e. \(\widehat{w}_k=0\)) since we define \(0/0=0\) whenever it appears. The last equation holds since the categories with \(\widehat{p}_k=0\) do not contribute to the estimation. Finally, we consider the doubly robust estimator. Note that \[\begin{align} &\, \mathbb{P}_n \left [ \frac{A\widehat{\mu}_{1X}}{\widehat{\pi}_X} \right] \\ =&\, \frac{1}{n} \sum_{i=1}^n \sum_{k:\widehat{p}_k>0} I(X_i=k) \frac{A_i \widehat{\mu}_{1X_i}}{\widehat{\pi}_{X_i}} \\ =&\, \sum_{k:\widehat{p}_k>0} \left( \frac{\sum_{i=1}^n A_i I(X_i=k)}{n \widehat{p}_k \widehat{\pi}_k} \right) \widehat{p}_k \widehat{\mu}_{1k}. \end{align}\] Note that \(n \widehat{p}_k \widehat{\pi}_k = \sum_{i=1}^n A_i I(X_i=k)\). If \(n \widehat{p}_k \widehat{\pi}_k>0\) then \[\frac{\sum_{i=1}^n A_i I(X_i=k)}{n \widehat{p}_k \widehat{\pi}_k} = 1.\] If \(n \widehat{p}_k \widehat{\pi}_k>0\) by our definition on \(0/0=0\) \[\frac{\sum_{i=1}^n A_i I(X_i=k)}{n \widehat{p}_k \widehat{\pi}_k} = 0, \,\widehat{\mu}_{1k} = \frac{\widehat{q}_{1k}}{\widehat{w}_k} = 0\] and hence we always have \[\left( \frac{\sum_{i=1}^n A_i I(X_i=k)}{n \widehat{p}_k \widehat{\pi}_k} \right) \widehat{p}_k \widehat{\mu}_{1k} = \widehat{p}_k \widehat{\mu}_{1k}.\] We conclude \[\mathbb{P}_n \left [ \frac{A\widehat{\mu}_{1X}}{\widehat{\pi}_X} \right] = \sum_{k: \widehat{p}_k >0} \widehat{p}_k \widehat{\mu}_{1k} = \sum_{k=1}^d \widehat{p}_k \widehat{\mu}_{1k}.\] This together with \[\mathbb{P}_n \left[\frac{AY}{\widehat{\pi}_X} \right] = \widehat{\psi}\] shows \[\widehat{\psi}_{1,\text{dr}} = \mathbb{P}_n \left[ \frac{A(Y-\widehat{\mu}_{1X})}{\widehat{\pi}_X} + \widehat{\mu}_{1X}\right] = \widehat{\psi}.\] ◻

11.2 Proof of Proposition 2 ↩︎

Proof. We focus on the bias of \(\widehat{\psi}_1\) and \(\widehat{\psi}_0\) can be similarly analyzed. We may rewrite the plug-in estimator as \[\widehat{\psi}_1 = \sum_{k=1}^d \frac{\widehat{q}_{1k}\widehat{p}_k}{\widehat{w}_k} I(\widehat{w}_k > 0)\] to emphasize the definition \(0/0=0\) in each term. On the event \(\{ n\widehat{w}_k = 0\}\), we automatically have \(\widehat{q}_{1k} = 0\). On the event \(\{ n\widehat{w}_k > 0\}\), \[\mathbb{E}[n\widehat{q}_{1k} \mid\mathbf{X}^n, \mathbf{A}^n] = n \widehat{w}_k \mu_{1k} = n\frac{q_{1k}}{w_{k}}\widehat{w}_k\] since \(n\widehat{q}_{1k} \mid \mathbf{X}^n, \mathbf{A}^n \sim \text{B}(n\widehat{w}_k, \mu_{1k})\). Hence we have \[\mathbb{E}\left[ \frac{\widehat{q}_{1k}\widehat{p}_k}{\widehat{w}_k} I(\widehat{w}_k > 0) \right] = \mathbb{E}\left[ \frac{\widehat{p}_k}{\widehat{w}_k} I(\widehat{w}_k > 0) \mathbb{E}(\widehat{q}_{1k}\mid\mathbf{X}^n, \mathbf{A}^n) \right] = \mathbb{E}[\widehat{p}_k I(\widehat{w}_k>0)]\frac{q_{1k}}{w_k}.\] The bias of each individual term is \[\begin{align} &\, \mathbb{E}\left[ \frac{\widehat{q}_{1k}\widehat{p}_k}{\widehat{w}_k} I(\widehat{w}_k > 0) \right]-\frac{q_{1k}p_k}{w_k} \\ = &\, \mathbb{E}[\widehat{p}_k I(\widehat{w}_k>0)] \frac{q_{1k}}{w_k} - \frac{q_{1k}p_k}{w_k} \\ = &\, -\mathbb{E}[\widehat{p}_k I(\widehat{w}_k=0)]\frac{q_{1k}}{w_{k}} \end{align}\] Further note that \(n\widehat{w}_k \mid \mathbf{X}^n \sim \text{B}(n\widehat{p}_k, \pi_k)\), \[\begin{align} &\, \mathbb{E}[\widehat{p}_k I(\widehat{w}_k=0)] \\ = &\, \mathbb{E}[\widehat{p}_k I(n\widehat{w}_k=0)] \\ =&\, \mathbb{E}[\widehat{p}_k \mathbb{P}(n\widehat{w}_k =0\mid\mathbf{X}^n)] \\ = &\, \mathbb{E}[\widehat{p}_k (1-\pi_k)^{n\widehat{p}_k}]. \end{align}\] We then use the fact \(\mathbb{E}[Vc^V] = npc(1-p+cp)^{n-1}\) for \(V \sim \text{Bin}(n,p)\) and obtain \[\begin{align} &\, \mathbb{E}[\widehat{p}_k (1-\pi_k)^{n\widehat{p}_k}] \\ = &\, \frac{1}{n}\mathbb{E}[n\widehat{p}_k (1-\pi_k)^{n\widehat{p}_k}] \\ = &\, \frac{1}{n} \times n p_k (1-\pi_k) [1-p_k+(1-\pi_k)p_k]^{n-1} \\ = &\, p_k (1-\pi_k)(1-p_k \pi_k)^{n-1} \end{align}\] Hence the bias for individual terms is \[\mathbb{E}\left[ \frac{\widehat{q}_{1k}\widehat{p}_k}{\widehat{w}_k} I(\widehat{w}_k > 0) \right]-\frac{q_{1k}p_k}{w_k}= -\mu_{1k}p_k (1-\pi_k)(1-p_k \pi_k)^{n-1}\] and we conclude \[\mathbb{E}[\widehat{\psi}_1 - \psi_1] = -\sum_{k=1}^d \mu_{1k}p_k (1-\pi_k)(1-p_k \pi_k)^{n-1}.\] Similarly (or by symmetry) for \(\widehat{\psi}_0\) we have \[\mathbb{E}[\widehat{\psi}_0 - \psi_0] = -\sum_{k=1}^d \mu_{0k}p_k \pi_k(1-p_k+p_k \pi_k)^{n-1}\] ◻

11.3 Proof of Proposition 3 ↩︎

Proof. Recall that we use bold letters \(({\boldsymbol{p}}, \boldsymbol{\pi}, \boldsymbol{\mu}_1)\) to denote the vectors \((p_k, \pi_k, \mu_{1k})_{k=1}^d\). We have \[\begin{align} \sup_{\mathbb{P}\in \mathcal{D}(\epsilon)}|\mathbb{E}_{\mathbb{P}}[\widehat{\psi}_1 - \psi_1]| =&\, \sup_{{\boldsymbol{p}}, \boldsymbol{\pi}, \boldsymbol{\mu}_1} \sum_{k=1}^d \mu_{1k}p_k (1-\pi_k)(1-p_k \pi_k)^{n-1} \\ =&\, \sup_{ {\boldsymbol{p}}, \boldsymbol{\pi}} \sum_{k=1}^d p_k (1-\pi_k)(1-p_k \pi_k)^{n-1} \\ = &\, \sup_{ {\boldsymbol{p}}} \sum_{k=1}^d p_k (1-\epsilon)(1-\epsilon p_k )^{n-1}, \end{align}\] where the first equation follows from \(\mu_{1k} \leq 1\), and the second follows from \(\pi_k \geq \epsilon\). The upper bound then follows from \[\max_{p\ge 0} p(1-\epsilon p)^{n-1} = \frac{1}{\epsilon n}\left(1-\epsilon\cdot \frac{1}{\epsilon n}\right)^{n-1} < \frac{1}{\epsilon n},\] where the maximizer \(p^\star = 1/(\epsilon n)\) follows from simple differentiation. The lower bound follows from the enumeration of three different scenarios:

If \(n\epsilon < 1\), choosing \({\boldsymbol{p}}= (1,0,\cdots,0)\) gives the lower bound \[\begin{align} (1-\epsilon)\cdot (1-\epsilon)^{n-1} \ge (1-\epsilon) \left(1-\frac{1}{n}\right)^{n-1} \ge \frac{1-\epsilon}{e}. \end{align}\]
If \(1\le n\epsilon<d-1\), choosing \({\boldsymbol{p}}= (1/(n\epsilon), \cdots, 1/(n\epsilon), 1-\lfloor n\epsilon \rfloor / (n\epsilon), 0, \cdots, 0)\) gives the lower bound \[\begin{align} \lfloor n\epsilon \rfloor \cdot \frac{1-\epsilon}{n\epsilon}\left(1-\frac{1}{n}\right)^{n-1} \ge \frac{n\epsilon}{2}\cdot \frac{1-\epsilon}{en\epsilon} \ge \frac{1-\epsilon}{2e}. \end{align}\]
If \(n\epsilon \ge d-1\), choosing \({\boldsymbol{p}}= (1/(n\epsilon), \cdots, 1/(n\epsilon), 1-(d-1)/(n\epsilon))\) gives the lower bound \[\begin{align} (d-1)\cdot \frac{1-\epsilon}{n\epsilon}\left(1-\frac{1}{n}\right)^{n-1} \ge \frac{1-\epsilon}{e}\cdot \frac{d-1}{n\epsilon}. \end{align}\]

◻

11.4 Proof of Theorem 1 ↩︎

Proof. \[\label{eq:var-psi1} \operatorname{Var}(\widehat{\psi}_1) = \mathbb{E}[\operatorname{Var}(\widehat{\psi}_1 \mid \mathbf{X}^n, \mathbf{A}^n)] + \operatorname{Var}(\mathbb{E}[\widehat{\psi}_1\mid\mathbf{X}^n, \mathbf{A}^n]).\tag{19}\] We start with the first term. By the conditional independence of \(\widehat{q}_{11}, \dots, \widehat{q}_{1d}\) given \(\mathbf{X}^n, \mathbf{A}^n\), we have \[\begin{align} &\, \operatorname{Var}(\widehat{\psi}_1 \mid \mathbf{X}^n, \mathbf{A}^n)\\ = &\, \operatorname{Var}\left( \sum_{k=1}^d \frac{\widehat{q}_{1k}\widehat{p}_k I(\widehat{w}_k>0)}{\widehat{w}_k} \mid \mathbf{X}^n, \mathbf{A}^n \right) \\ =&\, \sum_{k=1}^d \frac{\widehat{p}_k^2 I(\widehat{w}_k >0)}{\widehat{w}_k^2} \operatorname{Var}(\widehat{q}_{1k}\mid\mathbf{X}^n, \mathbf{A}^n) \\ = &\, \sum_{k=1}^d \frac{\widehat{p}_k^2 I(\widehat{w}_k >0)}{\widehat{w}_k^2} \frac{\widehat{w}_k \mu_{1k}(1-\mu_{1k})}{n} \\ = &\, \frac{1}{n} \sum_{k=1}^d \frac{\widehat{p}_k^2 I(\widehat{w}_k >0)}{\widehat{w}_k}\mu_{1k}(1-\mu_{1k}) \end{align}\] where we use the fact \(n \widehat{q}_{1k} \mid \mathbf{X}^n, \mathbf{A}^n \sim \text{Bin}(n\widehat{w}_k, \mu_{1k})\). We need the following lemma to proceed the analysis.

Lemma 2. (Lemma A.2 in [79])If \(X \sim B(n,p)\), then \[\mathbb{E}\left\{\frac{I(X>0)}{X}\right\} \leq \frac{2}{(n+1) p}.\]

Note that \(n\widehat{w}_k \sim \text{Bin}(n\widehat{p}_k, \pi_k) \mid \mathbf{X}^n\) and apply Lemma 2, we have \[\begin{align} &\, \mathbb{E}[\operatorname{Var}(\widehat{\psi}_1 \mid \mathbf{X}^n, \mathbf{A}^n)] \\ = &\, \frac{1}{n} \sum_{k=1}^d \mathbb{E}\left[ \frac{\widehat{p}_k^2 I(\widehat{w}_k >0)}{\widehat{w}_k} \right]\mu_{1k}(1-\mu_{1k}) \\ = &\, \sum_{k=1}^d \mathbb{E}\left\{ \widehat{p}_k^2 \mathbb{E}\left[ \frac{ I(\widehat{w}_k >0)}{n\widehat{w}_k} \mid \mathbf{X}^n\right] \right\}\mu_{1k}(1-\mu_{1k}) \\ \leq &\, \sum_{k=1}^d \mathbb{E}\left\{ \widehat{p}_k^2 \frac{2}{(n \widehat{p}_k +1)\pi_k} \right\}\mu_{1k}(1-\mu_{1k})\\ \leq &\, \frac{2}{n} \sum_{k=1}^d \mathbb{E}\left\{ \frac{\widehat{p}_k}{\pi_k} \right\}\mu_{1k}(1-\mu_{1k}) \\ \leq &\, \frac{1}{2n\epsilon} \sum_{k=1}^d \mathbb{E}[\widehat{p}_k] \\ = &\, \frac{1}{2n\epsilon}, \end{align}\] where we use the fact \(\mu_{1k}(1-\mu_{1k}) \leq 1/4\). Then we evaluate the second term in 19 \[\mathbb{E}[\widehat{\psi}_1 \mid \mathbf{X}^n, \mathbf{A}^n] = \sum_{k=1}^d \mathbb{E}\left[ \frac{\widehat{q}_{1k}\widehat{p}_k}{\widehat{w}_k}I(\widehat{w}_k >0) \mid \mathbf{X}^n, \mathbf{A}^n \right] = \sum_{k=1}^d \widehat{p}_kI(\widehat{w}_k >0) \mu_{1k}.\] By further using the property of conditional variance we have \[\label{eq:var-psi195cond} \operatorname{Var}(\mathbb{E}[\widehat{\psi}_1 \mid \mathbf{X}^n, \mathbf{A}^n]) = \mathbb{E}\left[\operatorname{Var}\left(\sum_{k=1}^d \widehat{p}_kI(\widehat{w}_k >0) \mu_{1k} \mid \mathbf{X}^n \right)\right] + \operatorname{Var}\left(\mathbb{E}\left[\sum_{k=1}^d \widehat{p}_kI(\widehat{w}_k >0) \mu_{1k} \mid \mathbf{X}^n \right]\right)\tag{20}\] We analyze the first term in 20 . By the conditional independence of \(n\widehat{w}_1 , \dots, n \widehat{w}_d\) given \(\mathbf{X}^n\) we have \[\begin{align} &\, \operatorname{Var}\left(\sum_{k=1}^d \widehat{p}_kI(\widehat{w}_k >0) \mu_{1k} \mid \mathbf{X}^n \right)\\ =&\, \sum_{k=1}^d \widehat{p}_k^2 \mu_{1k}^2 \operatorname{Var}(I(\widehat{w}_k > 0)\mid\mathbf{X}^n)\\ \leq &\, \sum_{k=1}^d \widehat{p}_k^2 \mu_{1k}^2 (1-\pi_k)^{n\widehat{p}_k} \\ \leq &\, \sum_{k=1}^d \widehat{p}_k^2 \mu_{1k}^2 \frac{1}{\pi_k(1+n\widehat{p}_k)} \\ \leq &\, \frac{1}{\epsilon n} \sum_{k=1}^d \widehat{p}_k \\ = &\, \frac{1}{\epsilon n}, \end{align}\] where we use the fact \[(1-x)^n \leq \frac{1}{x(1+n)}, \, \forall \, 0\leq x \leq 1.\] Hence we have \[\mathbb{E}\left[\operatorname{Var}\left(\sum_{k=1}^d \widehat{p}_kI(\widehat{w}_k >0) \mu_{1k} \mid \mathbf{X}^n \right)\right] \leq \frac{1}{\epsilon n}.\] For the second term in 20 , we have \[\begin{align} &\, \mathbb{E}\left[\sum_{k=1}^d \widehat{p}_kI(\widehat{w}_k >0) \mu_{1k} \mid \mathbf{X}^n \right] \\ = &\, \sum_{k=1}^d \widehat{p}_k \mu_{1k}(1-\mathbb{P}(n\widehat{w}_k=0\mid\mathbf{X}^n)) \\ = &\, \sum_{k=1}^d \mu_{1k} \widehat{p}_k [1-(1-\pi_k)^{n\widehat{p}_k}] \end{align}\] So we need to evaluate \[\begin{align} &\, \operatorname{Var}\left(\mathbb{E}\left[\sum_{k=1}^d \widehat{p}_kI(\widehat{w}_k >0) \mu_{1k} \mid \mathbf{X}^n \right]\right) \\ = &\, \operatorname{Var}\left( \sum_{k=1}^d \mu_{1k} \widehat{p}_k [1-(1-\pi_k)^{n\widehat{p}_k}]\right) \\ = &\, \sum_{k=1}^d \mu_{1k}^2 \operatorname{Var}\left(\widehat{p}_k [1-(1-\pi_k)^{n\widehat{p}_k}]\right) + \sum_{i\neq j} \mu_{1i}\mu_{1j}\operatorname{Cov}\left(\widehat{p}_i [1-(1-\pi_i)^{n\widehat{p}_i}], \widehat{p}_j [1-(1-\pi_j)^{n\widehat{p}_j}]\right). \end{align}\] We note that \[\mathbb{E}[\widehat{p}_k [1-(1-\pi_k)^{n\widehat{p}_k}]] = p_k - p_k(1-\pi_k)(1-\pi_k p_k)^{n-1}\] and \[\mathbb{E}\left\{\widehat{p}_k^2 [1-(1-\pi_k)^{n\widehat{p}_k}]^2\right\} \leq \mathbb{E}[\widehat{p}_k^2] = \frac{p_k(1-p_k)}{n} + p_k^2.\] This yields the following bound \[\begin{align} &\, \operatorname{Var}(\widehat{p}_k [1-(1-\pi_k)^{n\widehat{p}_k}]) \\ \leq &\, \frac{p_k(1-p_k)}{n} + 2p_k^2 (1-\pi_k)(1-\pi_k p_k)^{n-1} \\ \leq &\, \frac{p_k(1-p_k)}{n} + 2p_k^2 (1-\pi_k) \frac{1}{n\pi_k p_k} \\ \leq &\, \frac{p_k(1-p_k)}{n} + \frac{2p_k}{n \epsilon} \\ \leq &\, \left( 1 + \frac{2}{\epsilon} \right) \frac{p_k}{n}. \end{align}\] So we have \[\sum_{k=1}^d \mu_{1k}^2 \operatorname{Var}\left(\widehat{p}_k [1-(1-\pi_k)^{n\widehat{p}_k}]\right) \leq \left( 1 + \frac{2}{\epsilon} \right) \frac{1}{n}.\] For the covariance part, the computations are more involved. We need to compute \[\begin{align} &\, \mathbb{E}\left\{ \widehat{p}_i[1-(1-\pi_i)^{n\widehat{p}_i}] \widehat{p}_j [1-(1-\pi_j)^{n\widehat{p}_j}] \right\} \\ = &\, \mathbb{E}[\widehat{p}_i \widehat{p}_j] -\mathbb{E}[\widehat{p}_i\widehat{p}_j (1-\pi_i)^{n\widehat{p}_i}] -\mathbb{E}[\widehat{p}_i\widehat{p}_j (1-\pi_j)^{n\widehat{p}_j}] + \mathbb{E}[\widehat{p}_i\widehat{p}_j (1-\pi_i)^{n\widehat{p}_i} (1-\pi_j)^{n\widehat{p}_j}] \end{align}\] For three-dimensional multinomial distribution \((X_1, X_2, X_3)\sim\) Multinomial\((n, p_1, p_2, p_3)\) the probability generating function is \[\mathbb{E}[z_1^{X_1}z_2^{X_2}z_3^{X_3}] = (p_1z_1 + p_2z_2 + p_3z_3)^n.\] From this formula and differentiation we have \[\mathbb{E}[X_1X_2 z_1^{X_1}z_2^{X_2}] = n(n-1)p_1 z_1 p_2 z_2 (p_1z_1 + p_2z_2 + p_3)^{n-2}.\] We can obtain the expectations appearing in the covariance as: \[\mathbb{E}[\widehat{p}_i \widehat{p}_j] = \frac{n-1}{n}p_i p_j\] \[\mathbb{E}[\widehat{p}_i\widehat{p}_j (1-\pi_i)^{n\widehat{p}_i}] = \frac{n-1}{n}p_i p_j (1-\pi_i)(1-\pi_i p_i)^{n-2}\] \[\mathbb{E}[\widehat{p}_i\widehat{p}_j (1-\pi_i)^{n\widehat{p}_i} (1-\pi_j)^{n\widehat{p}_j}] = \frac{n-1}{n} p_i p_j (1-\pi_i)(1-\pi_j)(1-\pi_i p_i - \pi_j p_j)^{n-2}\] Now the covariance of pair \((i,j)\) is \[\begin{align} &\, \operatorname{Cov}\left(\widehat{p}_i [1-(1-\pi_i)^{n\widehat{p}_i}], \widehat{p}_j [1-(1-\pi_j)^{n\widehat{p}_j}]\right)\\ = &\, \mathbb{E}\left[ \widehat{p}_i[1-(1-\pi_i)^{n\widehat{p}_i}] \widehat{p}_j [1-(1-\pi_j)^{n\widehat{p}_j}] \right] - \mathbb{E}\left[ \widehat{p}_i[1-(1-\pi_i)^{n\widehat{p}_i}] \right] \mathbb{E}\left[ \widehat{p}_j [1-(1-\pi_j)^{n\widehat{p}_j}]\right] \\ = &\, \frac{n-1}{n}p_i p_j - \frac{n-1}{n}p_i p_j (1-\pi_i)(1-\pi_i p_i)^{n-2} - \frac{n-1}{n}p_i p_j (1-\pi_j)(1-\pi_j p_j)^{n-2} \\ &\, +\frac{n-1}{n} p_i p_j (1-\pi_i)(1-\pi_j)(1-\pi_i p_i - \pi_j p_j)^{n-2} \\ &\, - [p_i - p_i(1-\pi_i)(1-\pi_i p_i)^{n-1}][p_j - p_j(1-\pi_j)(1-\pi_j p_j)^{n-1}] \\ =&\, -\frac{1}{n}p_ip_j + p_i p_j (1-\pi_i)(1-\pi_ip_i)^{n-2}\left( \frac{1}{n} - \pi_i p_i \right) + p_i p_j (1-\pi_j)(1-\pi_jp_j)^{n-2}\left( \frac{1}{n} - \pi_j p_j \right) \\ &\, + \frac{1}{n}p_i p_j (1-\pi_i)(1-\pi_j)\left[(n-1)(1-\pi_ip_i-\pi_jp_j)^{n-2} - n(1-\pi_ip_i)^{n-1}(1-\pi_jp_j)^{n-1} \right] \end{align}\] where in the last equation we combine the four terms in \(\mathbb{E}\left[ \widehat{p}_i[1-(1-\pi_i)^{n\widehat{p}_i}] \widehat{p}_j [1-(1-\pi_j)^{n\widehat{p}_j}] \right]\) and \(\mathbb{E}\left[ \widehat{p}_i[1-(1-\pi_i)^{n\widehat{p}_i}] \right] \mathbb{E}\left[ \widehat{p}_j [1-(1-\pi_j)^{n\widehat{p}_j}]\right]\) correspondingly. We proceed as taking summations \[\begin{align} &\,\bigg| \sum_{i \neq j} \mu_{1i} \mu_{1j} \operatorname{Cov}\left(\widehat{p}_i [1-(1-\pi_i)^{n\widehat{p}_i}], \widehat{p}_j [1-(1-\pi_j)^{n\widehat{p}_j}]\right) \bigg| \\ \leq &\,\sum_{i \neq j} \frac{1}{n}p_i p_j + 2 \sum_{i \neq j}p_i p_j (1-\pi_i)(1-\pi_i p_i)^{n-2}\bigg| \frac{1}{n}-\pi_i p_i \bigg| \\ &\, + \sum_{i \neq j} \frac{1}{n}p_i p_j (1-\pi_i)(1-\pi_j)\big |(n-1)(1-\pi_ip_i-\pi_jp_j)^{n-2} - n(1-\pi_ip_i)^{n-1}(1-\pi_jp_j)^{n-1} \big | \end{align}\] For each individual term in the expression above, we have \[\sum_{i \neq j} \frac{1}{n}p_i p_j = \sum_{i} \frac{1}{n}p_i (1-p_i) \leq \frac{1}{n}.\] \[\begin{align} &\, \sum_{i \neq j}p_i p_j (1-\pi_i)(1-\pi_i p_i)^{n-2}\bigg| \frac{1}{n}-\pi_i p_i\bigg|\\ \leq &\, \sum_{i \neq j} \frac{1}{n} p_i p_j + \sum_{i \neq j} p_i p_j (1-\pi_i)(1-\pi_i p_i)^{n-2} \pi_i p_i\\ \leq &\, \frac{1}{n} + \sum_{i \neq j} p_i p_j (1-\pi_i)\frac{1}{\pi_i p_i (n-1)} \pi_i p_i \\ \leq &\, \frac{1}{n}+ \frac{1}{n-1}. \end{align}\] For the last term we need an auxiliary lemma as follows:

Lemma 3. For \(n \geq 3\), consider the function \[f_2(x,y) = (n-1)(1-x-y)^{n-2} - n(1-x)^{n-1}(1-y)^{n-1}\] defined on the triangle \(D = \{ x \geq 0, y \geq 0, x+y\leq 1\}\). Then we have \[|f_2(x,y)| \leq 1.\]

For the last term we invoke Lemma 3 and have \[\begin{align} &\, \sum_{i \neq j} \frac{1}{n}p_i p_j (1-\pi_i)(1-\pi_j)\big |(n-1)(1-\pi_ip_i-\pi_jp_j)^{n-2} - n(1-\pi_ip_i)^{n-1}(1-\pi_jp_j)^{n-1} \big |\\ \leq &\, \sum_{i \neq j} \frac{1}{n}p_i p_j \leq \frac{1}{n}. \end{align}\] We thus showed \[\big| \sum_{i \neq j} \mu_{1i} \mu_{1j} \operatorname{Cov}\left(\widehat{p}_i [1-(1-\pi_i)^{n\widehat{p}_i}], \widehat{p}_j [1-(1-\pi_j)^{n\widehat{p}_j}]\right) \big| \leq \frac{4}{n} + \frac{2}{n-1}.\] \[\operatorname{Var}\left(\mathbb{E}\left[\sum_{k=1}^d \widehat{p}_kI(\widehat{w}_k >0) \mu_{1k} \bigg | \mathbf{X}^n \right]\right) \leq \left( 1 + \frac{2}{\epsilon} \right) \frac{1}{n} + \frac{4}{n} + \frac{2}{n-1}.\] \[\operatorname{Var}(\mathbb{E}[\widehat{\psi}_1 \mid \mathbf{X}^n, \mathbf{A}^n]) \leq \left( 1 + \frac{3}{\epsilon} \right) \frac{1}{n} + \frac{4}{n} + \frac{2}{n-1}.\] \[\operatorname{Var}(\widehat{\psi}_1) \leq \left( 1 + \frac{7}{2\epsilon} \right) \frac{1}{n} + \frac{4}{n} + \frac{2}{n-1}.\] ◻

11.5 Proof of Theorem 2 ↩︎

Proof. The proof strategy follows from the moment matching method commonly used in theoretical computer science literature [37]–[39]. We first define a “relaxed" model class with proportion vector \({\boldsymbol{p}}= (p_1, p_2,\dots)\) being approximately a probability vector. Mathematically, for \(\delta>0\) define \[\label{eq:model-class-relaxed} \mathcal{D}(\epsilon, \delta) = \left\{ \left|\sum_{k=1}^d p_k -1\right| \leq \delta, p_k \in \left[0,1\right], \pi_k \in [\epsilon, 1-\epsilon], \forall k \in [d] \right\},\tag{21}\] where we relax the assumption \(\sum_k p_k=1\) to \(\left|\sum_k p_k -1\right| \leq \delta\) so that we can set \(p_k\)’s to be random variables in our construction and use the method of fuzzy hypotheses [49], [80]. Under \(\mathbb{P}\in \mathcal{D}(\epsilon, \delta)\), we assume \(n\widehat{p}_k \sim \text{Poi}(np_k)\) and \((n\widehat{p}_1, n\widehat{p}_2, \dots)\) are independent, i.e. we again rely on a Poisson sampling model to prove our results. The treatment assignment \(A\) and outcome \(Y\) have the same distribution as in Section 2.1 conditioned on the category each sample falls into. Define the minimax lower bound over \(\mathcal{D}(\epsilon, \delta)\) as \[\widetilde{R}^*(d,n,\delta) := \inf_{\widehat{\psi}_1} \sup_{\mathbb{P}\in \mathcal{D}(\epsilon, \delta)} \mathbb{E}_{\mathbb{P}}\left[\left(\widehat{\psi}_1 - \psi_1\right)^2\right]\] The following lemma allows us to relate \({R}^*(d,n)\) with \(\widetilde{R}^*(d,n,\delta)\).

Lemma 4. For any \(\delta \in [0,1/3)\), \[{R}^*(d,n/2) \geq \frac{1}{2} \widetilde{R}^*(d,n,\delta) - \exp(-n/50)-\delta^2.\]

We then present an auxiliary result characterizing the prior distribution on the parameters in the method of fuzzy hypotheses. In the following proof of this section, we will set \(\mu_{0k}=0\) and abbreviate \(\mu_{1k}\) as \(\mu_k\), which is different from the mean of outcome in \(k\)-th category \(\mathbb{E}[Y|X=k]\) as in Section 5.

Lemma 5. There exists constants \(c, c' >0\) such that for any constants \(c_1, c_2, c_3>0\) satisfying \(\frac{c c_1 c_3}{c_2^2} \leq 1\) and \(d \leq c_3 n \log n\), there exist two distributions \(\mu_0\) and \(\mu_1\) on \((p, \pi, \mu )\) satisfying the following properties:

\(\mu_0\) a.s. and \(\mu_1\) a.s. \[0 \leq p \leq \frac{c_1 \log n}{n}, \, \epsilon \leq \pi \leq 1-\epsilon, \, 0 \leq \mu \leq 1;\]
For \(i, j ,k \geq 0\) and \(i+j+k \leq 3K\) with \(K = c_2 \log n\), \[\mathbb{E}_{\mu_0}\left[ p^i (p \pi)^j (p \pi \mu)^k \right] = \mathbb{E}_{\mu_1} \left[ p^i (p \pi)^j (p \pi \mu)^k \right].\]
\[a_0 : =\mathbb{E}_{\mu_0}[p] = \mathbb{E}_{\mu_1}[p] \leq \frac{1}{d}.\]
\[|\mathbb{E}_{\mu_0}[p \mu] - \mathbb{E}_{\mu_1}[p \mu]| \geq \frac{c_4}{n \log n},\] where \(c_4 = \frac{cc'c_1}{c_2^2}\).

The proof of Lemma 4 and 5 are provided in Section 12. Note that in the construction of \(\mu_i\), \(\pi\) and \(\mu\) are actually functions of \(p\). Under null hypothesis \(P\), let \((p_k, \pi_k, \mu_k)\) i.i.d. \(\sim \mu_0, 1\leq k \leq d\). We add one more category with \(p_{d+1}=1-d a_0, \pi_{d+1}=\epsilon, \mu_{d+1}=0\). Under the alternative hypothesis \(P'\), let \((p_k', \pi_k', \mu_k')\) i.i.d. \(\sim \mu_1, 1\leq k \leq d\) and \(p_{d+1}'=1-d a_0, \pi_{d+1}'=\epsilon, \mu_{d+1}'=0\). Obviously, adding one more category will not affect the final rate. The sufficient statistics for \(\psi_1\) are \[\#\{i:X_i=k,A_i=1,Y_i=1 \}, \#\{i:X_i=k,A_i=1,Y_i=0 \}, \#\{i:X_i=k,A_i=0 \}, 1\leq k \leq d+1.\] By the property of Poisson distribution these counting statistics are independent under the Poisson sampling model and under the null hypothesis \(P\) (given \(\{p_1, \dots, p_d\}\), which is equivalent to conditioning on \(\{(p_k, \pi_k, \mu_k ), 1\leq k \leq d\}\) since \(\pi, \mu\) are functions of \(p\)) \[N_{k11} : = \#\{i:X_i=k,A_i=1,Y_i=1 \} \sim \text{Poi} (np_k \pi_k\mu_{k}) ,\] \[N_{k10} : = \#\{i:X_i=k,A_i=1,Y_i=0 \} \sim \text{Poi} (np_k \pi_k(1-\mu_{k})) ,\] \[N_{k0} := \#\{i:X_i=k,A_i=0 \} \sim \text{Poi} (np_k (1-\pi_k)) .\] Similarly under the alternative hypothesis \(P'\), \[N_{k11}' : = \#\{i:X_i=k,A_i=1,Y_i=1 \} \sim \text{Poi} (np_k' \pi_k'\mu_{k}') ,\] \[N_{k10}' : = \#\{i:X_i=k,A_i=1,Y_i=0 \} \sim \text{Poi} (np_k' \pi_k'(1-\mu_{k}')) ,\] \[N_{k0}' := \#\{i:X_i=k,A_i=0 \} \sim \text{Poi} (np_k' (1-\pi_k')) .\] Denote \(\mathbf{N}_k = (N_{k11}, N_{k10}, N_{k0})\) as the sufficient statistics in \(k\)-th category and \(\mathbf{N}=(\mathbf{N}_1,\dots, \mathbf{N}_{d}, \mathbf{N}_{d+1})\) as the collection of these sufficient statistics. Define \(\mathbf{N}'\) similarly for the alternative hypothesis. The total variation distance between the marginal distribution of \(\mathbf{N}\) and \(\mathbf{N}'\) (marginalize over the distribution of \(\{(p_k, \pi_k, \mu_k ), 1\leq k \leq d\}\) and \(\{(p_k', \pi_k', \mu_k' ), 1\leq k \leq d\}\)) can be bounded by triangle inequality as (since \(\mathbf{N}_k\) depends only on \((p_k, \pi_k, \mu_k )\) for each \(k\) and \(\{(p_k, \pi_k, \mu_k ), 1\leq k \leq d\}\) are independent, \((\mathbf{N}_1,\dots, \mathbf{N}_d, \mathbf{N}_{d+1})\) are also independent) \[\text{TV}(\mathbf{N},\mathbf{N}') \leq\, \sum_{k=1}^d \text{TV}(\mathbf{N}_k,\mathbf{N}_k').\] Note that the marginal distributions of components of \(\mathbf{N}_k\) are not independent since \(N_{k11}, N_{k10}, N_{k0}\) all depend on \((p_k, \pi_k, \mu_k )\). Conditioned on \((p_k, \pi_k, \mu_k )\), \((N_{k11}, N_{k10}, N_{k0})\) are conditionally independent with each being a Poisson distribution. By definition of total variation distance, \[\begin{align} &\, \text{TV}(\mathbf{N}_k,\mathbf{N}_k')\\ = &\, \frac{1}{2} \sum_{i,j,\ell=0}^{\infty} \big | \mathbb{P}(N_{k11}=i, N_{k10}=j, N_{k0}=\ell)- \mathbb{P}(N_{k11}'=i, N_{k10}'=j, N_{k0}'=\ell)\big |. \end{align}\] Conditioning on \((p_k, \pi_k, \mu_k )\) we have \[\begin{align} &\,\mathbb{P}(N_{k11}=i, N_{k10}=j, N_{k0}=\ell)\\ = & \, \mathbb{E}\left [ \exp(-np_k \pi_k \mu_k) \frac{(np_k \pi_k \mu_k)^i}{i!} \exp(-np_k \pi_k(1- \mu_k)) \frac{[np_k \pi_k(1- \mu_k)]^j}{j!} \right.\\ &\, \left.\exp(-np_k (1-\pi_k)) \frac{[np_k(1-\pi_k)]^{\ell}}{\ell !} \right ] \\ = &\, \mathbb{E}\left[ \exp(-n p_k) \frac{(np_k \pi_k \mu_k)^i[np_k \pi_k(1- \mu_k)]^j[np_k(1-\pi_k)]^{\ell}}{i!j!\ell!} \right] \\ = & \, \mathbb{E}\left[ \sum_{t=0}^{\infty} \frac{(-np_k)^t(np_k \pi_k \mu_k)^i[np_k \pi_k(1- \mu_k)]^j[np_k(1-\pi_k)]^{\ell}}{i!j!\ell!t!} \right]. \end{align}\] Hence the total variation distance can be written as \[\begin{align} &\, \text{TV}(\mathbf{N}_k,\mathbf{N}_k') \\ = &\, \frac{1}{2} \sum_{i,j,\ell} \left | \mathbb{E}\left[ \sum_{t=0}^{\infty} \frac{(-np_k)^t(np_k \pi_k \mu_k)^i[np_k \pi_k(1- \mu_k)]^j[np_k(1-\pi_k)]^{\ell}}{i!j!\ell!t!} \right] \right. \\ &\, \left.- \mathbb{E}\left[ \sum_{t=0}^{\infty} \frac{(-np_k')^t(np_k' \pi_k' \mu_k')^i[np_k' \pi_k'(1- \mu_k')]^j[np_k'(1-\pi_k')]^{\ell}}{i!j!\ell!t!} \right] \right |. \end{align}\] Note that \[(np_k)^t(np_k \pi_k \mu_k)^i[np_k \pi_k(1- \mu_k)]^j[np_k(1-\pi_k)]^{\ell}\] is a polynomial of \((p_k, p_k \pi_k, p_k\pi_k \mu_k)\) with degree \(t+i+j+\ell\). By moment matching property \(1\) in Theorem 5, we have \[\begin{align} &\, \text{TV}(\mathbf{N}_k,\mathbf{N}_k') \\ \leq &\, \frac{1}{2} \left\{ \sum_{i+j+\ell+t > 3K} \mathbb{E}\left[ \frac{(np_k)^t(np_k \pi_k \mu_k)^i[np_k \pi_k(1- \mu_k)]^j[np_k(1-\pi_k)]^{\ell}}{i!j!\ell!t!} \right] \right. \\ &\, \left. + \mathbb{E}\left[ \frac{(np_k')^t(np_k' \pi_k' \mu_k')^i[np_k' \pi_k'(1- \mu_k')]^j[np_k'(1-\pi_k')]^{\ell}}{i!j!\ell!t!} \right] \right\} . \end{align}\] Note that \[\begin{align} &\,\frac{(np_k)^t(np_k \pi_k \mu_k)^i[np_k \pi_k(1- \mu_k)]^j[np_k(1-\pi_k)]^{\ell}}{i!j!\ell!t!} \\ = &\, \exp(2np_k) \exp(-np_k) \frac{(np_k)^t}{t!} \exp(-np_k \pi_k \mu_k) \frac{(np_k \pi_k \mu_k)^i}{i!} \\ &\, \exp(-np_k \pi_k(1- \mu_k)) \frac{[np_k \pi_k(1- \mu_k)]^j}{j!} \exp(-np_k (1-\pi_k)) \frac{[np_k(1-\pi_k)]^{\ell}}{\ell !} \\ = & \, \exp(2np_k)\mathbb{P}(V_1 = t, V_2 = i, V_3 = j, V_4 = \ell \mid p_k, \pi_k, \mu_k), \end{align}\] where conditioning on \((p_k, \pi_k, \mu_k)\), \[\begin{align} V_1 \sim &\, \text{Poi}(n p_k) ,\\ V_2 \sim &\, \text{Poi}(n p_k \pi_k \mu_k) ,\\ V_3 \sim &\, \text{Poi}(n p_k \pi_k (1-\mu_k)) ,\\ V_4 \sim &\, \text{Poi}(n p_k (1-\pi_k)) , \end{align}\] and \((V_1, V_2, V_3, V_4)\) are independent. Further let \(V= V_1 + V_2 + V_3 + V_4 \sim \text{Poi}(2n p_k)\) given \((p_k, \pi_k, \mu_k)\). Thus we have \[\begin{align} &\, \sum_{i+j+\ell+t > 3K} \mathbb{E}\left[ \frac{(np_k)^t(np_k \pi_k \mu_k)^i[np_k \pi_k(1- \mu_k)]^j[np_k(1-\pi_k)]^{\ell}}{i!j!\ell!t!} \right]\\ = &\, \mathbb{E}[\exp(2np_k) \mathbb{P}(V_1+V_2+V_3+V_4 > 3K \mid p_k, \pi_k,\mu_k)] \\ = &\, \mathbb{E}[\exp(2np_k) \mathbb{P}(V > 3K \mid p_k, \pi_k,\mu_k)] \\ = &\, \mathbb{E}\left[\exp(2np_k) \sum_{\ell > 3K} \exp(-2np_k) \frac{(2np_k)^{\ell}}{\ell!}\right] \\ = &\, \mathbb{E}\left[\sum_{\ell > 3K} \frac{(2np_k)^{\ell}}{\ell!} \right] \\ \leq &\, \sum_{\ell > 3K} \frac{(2c_1 \log n)^{\ell}}{\ell!}\\ = & \, \exp(2c_1 \log n) \sum_{\ell > 3K} \exp(-2c_1 \log n) \frac{(2c_1 \log n)^{\ell}}{\ell!} \\ = &\, n^{2c_1} \mathbb{P}(W>3K), \end{align}\] where \(W \sim \text{Poi}(2c_1 \log n)\). Apply the following Chernoff bound: For \(L>eM\) \[\label{eq:chernoff-poi} \mathbb{P}(\text{Poi}(M)>L) \leq \exp(-M) \left(\frac{eM}{L} \right)^L.\tag{22}\] When \(3K > 2e c_1 \log n\), we have \[\mathbb{P}(W>3K) \leq \exp(-2c_1 \log n) \left(\frac{2ec_1}{3c_2} \right)^{3c_2 \log n},\] \[n^{2c_1} \mathbb{P}(W>3K) \leq \left(\frac{2ec_1}{3c_2} \right)^{3c_2 \log n} = n^{3c_2 \log \left(\frac{2ec_1}{3c_2} \right)} \leq n^{-2} \leq \frac{c_3 \log n}{nd}\] as long as we choose constant \(c_1, c_2\) satisfying \(3c_2 \log \left(\frac{3c_2}{2ec_1} \right) \geq 2\) and \(c_3\) is the constant in Lemma 5. Thus the total variation distance is bounded as (similar inequality holds under the alternative hypothesis) \[\text{TV}(\mathbf{N}_k,\mathbf{N}_k') \leq \frac{c_3 \log n}{nd}\] \[\text{TV}(\mathbf{N},\mathbf{N}') \leq \frac{c_3 \log n}{n}.\] The functional separation is \[\psi_1(P) - \psi_1(P') = \sum_{k=1}^d (p_k \mu_k - p_k' \mu_k'),\] \[\left |\mathbb{E}\left[\psi_1(P) - \psi_1(P') \right] \right | = d | \mathbb{E}[p\mu - p' \mu'] | \geq \frac{c_4 d}{n \log n}:= q.\] Consider the following events: \[E=\left\{\bigg| \sum_{k=1}^d p_k -da_0\bigg|\leq \delta, \bigg | \sum_{k=1}^d p_k \mu_k - d\mathbb{E}_{\mu_0}[p \mu]\bigg|\leq q/4 \right\},\] \[E'=\left\{\bigg| \sum_{k=1}^d p_k' -da_0\bigg|\leq \delta, \bigg | \sum_{k=1}^d p_k' \mu_k' -d \mathbb{E}_{\mu_1}[p' \mu']\bigg|\leq q/4 \right\},\] By Chebyshev’s inequality, we have \[\begin{align} \mathbb{P}(E^c) \leq &\, \mathbb{P}\left(\bigg| \sum_{k=1}^d p_k -da_0\bigg|> \delta\right) \\ &\, + \mathbb{P}\left(\bigg | \sum_{k=1}^d p_k \mu_k - d\mathbb{E}_{\mu_0}[p \mu]\bigg| > q/4 \right) \\ \leq &\, \frac{d \operatorname{Var}_{\mu_0}(p)}{ \delta^2} + \frac{16 d \operatorname{Var}_{\mu_0}(p \mu)}{q^2} \\ \leq &\, \frac{c_1^2 d (\log n)^2}{ \delta^2 n^2} + \frac{16 c_1^2 d (\log n)^2}{n^2 q^2} \\ = &\, \frac{1}{16} + \frac{16c_1^2 (\log n)^4}{ c_4^2 d}, \end{align}\] where the third inequality follows from the bound \(|p\mu| \leq |p|\leq c_1 \log n/n\) and the last equation follows by setting \(\delta = \frac{4 c_1 \sqrt{d} \log n }{n}\). Similarly, we have \[\mathbb{P}(E'^c) \leq \frac{1}{16} + \frac{16c_1^2 (\log n)^4}{ c_4^2 d}.\] We put the following prior distributions induced by \(\{(p_k, \pi_k, \mu_k), 1 \leq k \leq d\}\) and \(\{(p_k', \pi_k', \mu_k'), 1 \leq k \leq d\}\) on \(P\) and \(P'\), respectively: \[\pi \stackrel{d}{=} P \mid E, \pi' \stackrel{d}{=} P' \mid E'.\] Note that under \(\pi, \pi'\), \[|\psi_1(P)-\psi_1(P')|\geq q/2.\] By triangle inequality, the total variation distance of the sufficient counting statistics \(\mathbf{N}\) and \(\mathbf{N}'\) under two priors is bounded by \[\begin{align} \text{TV}\left( {\mathbf{N}\mid E}, {\mathbf{N}' \mid E'} \right) \leq &\, \text{TV}\left( {\mathbf{N}\mid E}, {\mathbf{N}} \right) + \text{TV}({\mathbf{N}},{\mathbf{N}'} ) + \text{TV}\left( {\mathbf{N}' \mid E'}, {\mathbf{N}'} \right)\\ \leq &\, \mathbb{P}(E^c) + \mathbb{P}(E'^c) + \text{TV}({\mathbf{N}},{\mathbf{N}'} )\\ \leq &\, \frac{1}{8} + \frac{32c_1^2 (\log n)^4}{ c_4^2 d} + \frac{c_3 \log n}{n} . \end{align}\] By method of fuzzy hypotheses (Section 2.7.4 of [49]) we conclude \[\begin{align} &\, \widetilde{R}^*(d,n,\delta) \\ \geq &\, \frac{q^2}{32} \left( 1 - \text{TV}\left( {\mathbf{N}\mid E}, {\mathbf{N}' \mid E'} \right) \right).\\ \geq &\, \frac{q^2}{32} \left( \frac{7}{8} - \frac{32c_1^2 (\log n)^4}{ c_4^2 d} - \frac{c_3 \log n}{n} \right). \end{align}\] Hence in the regime \(d \gtrsim (\log n)^4\), we have \[\widetilde{R}^*(d,n,\delta) \gtrsim \frac{d^2}{(n \log n)^2}.\] By Lemma 4 we conclude \[\begin{align} {R}^*(d,n) \gtrsim &\, \frac{d^2}{(n \log n)^2} - \exp(-n/50)- \frac{d (\log n)^2}{n^2}\\ \gtrsim &\, \frac{d^2}{(n \log n)^2}. \end{align}\] The overall requirements on the constants are \[\frac{c c_1 c_3}{c_2^2} \leq 1, \, 3c_2 \log \left(\frac{3c_2}{2ec_1} \right) \geq 2.\] Clearly for \(c, c_3 >0\), one can choose \(c_1\) sufficiently small and \(c_2\) sufficiently large to satisfy these conditions.

The lower bound \(1/n\) can be proved using a two-point method. Without loss of generality, assume that \(n \geq 8\). Under the null \(P\) and alternative hypothesis \(P'\), set \[p_k = \frac{1}{d}, \pi_k = \frac{1}{2}, \mu_{0k}=0, \mu_{1k}(P) = \frac{1}{2}, \mu_{1k}(P') = \frac{1}{2} + \delta,\] with \(\delta = 1/\sqrt{n}\), i.e., the covariate distribution and propensity score are the same. By Le Cam’s two-point method we have \[\label{eq:two-point} \inf_{\widehat{\psi}_1} \sup_{\mathbb{P}\in \mathcal{D} (\epsilon)} \mathbb{E}_{\mathbb{P}}\left[\left(\widehat{\psi}_1 - \psi_1\right)^2\right] \geq \frac{1}{4}(\psi_1(P)-\psi_1(P'))^2 \exp (-n D(P \| P')).\tag{23}\] Note that the functional separation is \[|\psi_1(P)-\psi_1(P')| = \delta = \frac{1}{\sqrt{n}}.\] Under the null hypothesis \(P\), we have \[P(X=k, A=a, Y=y) = \left\{ \begin{array}{rl} \frac{1}{4d} & \text{if } a=1, y=1,\\ \frac{1}{4d} & \text{if } a=1, y=0,\\ \frac{1}{2d} & \text{if } a=0, y=0. \end{array} \right.\] Under the alternative hypothesis \(P'\), we have \[P'(X=k, A=a, Y=y) = \left\{ \begin{array}{rl} \frac{1}{2d}\left( \frac{1}{2}+\delta \right) & \text{if } a=1, y=1,\\ \frac{1}{2d} \left( \frac{1}{2}-\delta \right) & \text{if } a=1, y=0,\\ \frac{1}{2d} & \text{if } a=0, y=0. \end{array} \right.\] The K-L divergence between \(P\) and \(P'\) is \[D(P \| P') = \sum_{k=1}^d \left[- \frac{1}{4d} \log(1+2\delta) - \frac{1}{4d} \log(1-2\delta) \right] = -\frac{1}{4} \log (1-4\delta^2).\] Using the inequality \[\log(1-x) \geq -2x, x \in [0,1/2],\] we have \[D(P \| P') =-\frac{1}{4} \log \left(1-\frac{4}{n} \right) \leq \frac{2}{n}.\] Plug in the functional separation and bound on K-L divergence into 23 , we conclude \[\inf_{\widehat{\psi}_1} \sup_{\mathbb{P}\in \mathcal{D} (\epsilon)} \mathbb{E}_{\mathbb{P}}\left[\left(\widehat{\psi}_1 - \psi_1\right)^2\right] \gtrsim \frac{1}{n}.\] ◻

11.6 Proof of Lemma 1 ↩︎

Proof. By the definition of \(\widehat{t}_k = I(\widehat{p}_k >0, 0 < \widehat{\pi}_k < 1)\) we have \[\begin{align} \mathbb{P}\left(\sum_{j=1}^d \widehat{t}_j=0\right) & =\mathbb{P}\left(\widehat{w}_k \in\left\{0, \widehat{p}_k\right\} \forall k \text{ with } \widehat{p}_k>0\right) \\ & =\mathbb{E}\left\{\mathbb{P}\left(\widehat{w}_k \in\left\{0, \widehat{p}_k\right\} \forall k \text{ with } \widehat{p}_k>0 \mid \mathbf{X}^n\right)\right\} \\ & =\mathbb{E}\left\{\prod_{k: \widehat{p}_k > 0} \mathbb{P}\left(\widehat{w}_k \in\left\{0, \widehat{p}_k\right\} \mid \mathbf{X}^n\right)\right\} \\ & =\mathbb{E}\left(\prod_{k=1}^d\left[\left\{\left(1-\pi_k\right)^{n \widehat{p}_k}+\pi_k^{n \widehat{p}_k}\right\} I\left(\widehat{p}_k>0\right)+I\left(\widehat{p}_k=0\right)\right]\right) \\ & =\mathbb{E}\left[\prod_{k=1}^d\left\{\left(1-\pi_k\right)^{n \widehat{p}_k}+\pi_k^{n \widehat{p}_k}-I\left(\widehat{p}_k=0\right)\right\}\right] \end{align}\] The proof relies on the poissonization technique to bound the above expectation of the product. Poissonization allows us to replace \(\left(n \widehat{p}_1, \ldots, n \widehat{p}_d\right) \sim \operatorname{Multinomial}\left(n, p_1, \ldots, p_d\right)\) with \(n\widehat{p}_k \sim \text{Poisson}(n p_k)\) and \(n\widehat{p}_1, \dots, n\widehat{p}_d\) are independent. The following lemma connects the expectation in the multinomial case with that in the independent Poisson case.

Lemma 6 (Theorem 5.10 in [81]).

Let \(\mathbf{X}^n \in \mathbb{R}^d \sim \text{Multinomial}\left(n, p_1, \ldots, p_d\right)\), \(\mathbf{Y}^n \in \mathbb{R}^d\) and \(Y_{i}^n \sim \text{Poisson}(np_i)\) and \(Y_{1}^n, \dots, Y_{d}^n\) are independent. Consider a non-negative function \(f(x_1, \dots, x_d)\), if \(\mathbb{E}[f(X_1^n,\dots, X_d^n)]\) is monotonely non-increasing with \(n\), then \[\mathbb{E}[f(X_1^n,\dots, X_d^n)] \leq 2 \mathbb{E}[f(Y_1^n,\dots, Y_d^n)].\]

The proof is left as an exercise in [81] and we include it in Appendix 12.4. Let \[\begin{align} a_n =&\, \mathbb{E}\left[\prod_{k=1}^d\left\{\left(1-\pi_k\right)^{n \widehat{p}_k}+\pi_k^{n \widehat{p}_k}-I\left(\widehat{p}_k=0\right)\right\}\right] \\ =&\, \mathbb{E}\left[\prod_{k=1}^d\left\{\left(1-\pi_k\right)^{\sum_{i=1}^n I(X_i=k)}+\pi_k^{\sum_{i=1}^n I(X_i=k)}-I\left(\sum_{i=1}^n I(X_i=k) = 0\right)\right\}\right]. \end{align}\]

First we verify the monotonicity of \(a_n\). We claim \[\begin{align} &\, \left(1-\pi_k\right)^{\sum_{i=1}^{n+1} I(X_i=k)}+\pi_k^{\sum_{i=1}^{n+1} I(X_i=k)}-I\left(\sum_{i=1}^{n+1} I(X_i=k) = 0\right) \\ \leq &\, \left(1-\pi_k\right)^{\sum_{i=1}^{n} I(X_i=k)}+\pi_k^{\sum_{i=1}^{n} I(X_i=k)}-I\left(\sum_{i=1}^{n} I(X_i=k) = 0\right). \end{align}\] On the event \(\left\{\sum_{i=1}^{n+1 } I(X_i=k) = 0\right\}\) then \(\sum_{i=1}^{n } I(X_i=k) = 0\) and \[\left(1-\pi_k\right)^{\sum_{i=1}^{n+1} I(X_i=k)}+\pi_k^{\sum_{i=1}^{n+1} I(X_i=k)}-I\left(\sum_{i=1}^{n+1} I(X_i=k) = 0\right) = 1\] \[\left(1-\pi_k\right)^{\sum_{i=1}^{n} I(X_i=k)}+\pi_k^{\sum_{i=1}^{n} I(X_i=k)}-I\left(\sum_{i=1}^{n} I(X_i=k) = 0\right) = 1.\] On the event \(\left\{\sum_{i=1}^{n} I(X_i=k) = 0, X_{n+1} = k\right\}\), \[\left(1-\pi_k\right)^{\sum_{i=1}^{n+1} I(X_i=k)}+\pi_k^{\sum_{i=1}^{n+1} I(X_i=k)}-I\left(\sum_{i=1}^{n+1} I(X_i=k) = 0\right) = 1-\pi_k + \pi_k = 1\] \[\left(1-\pi_k\right)^{\sum_{i=1}^{n} I(X_i=k)}+\pi_k^{\sum_{i=1}^{n} I(X_i=k)}-I\left(\sum_{i=1}^{n} I(X_i=k) = 0\right) = 1.\] On the event \(\left\{\sum_{i=1}^n I(X_i=k) > 0 \right\}\), \[\begin{align} &\, \left(1-\pi_k\right)^{\sum_{i=1}^{n+1} I(X_i=k)}+\pi_k^{\sum_{i=1}^{n+1} I(X_i=k)}-I\left(\sum_{i=1}^{n+1} I(X_i=k) = 0\right) \\ = &\, \left(1-\pi_k\right)^{\sum_{i=1}^{n+1} I(X_i=k)}+\pi_k^{\sum_{i=1}^{n+1} I(X_i=k)} \\ \leq &\, \left(1-\pi_k\right)^{\sum_{i=1}^{n} I(X_i=k)}+\pi_k^{\sum_{i=1}^{n} I(X_i=k)} \\ = &\, \left(1-\pi_k\right)^{\sum_{i=1}^{n} I(X_i=k)}+\pi_k^{\sum_{i=1}^{n} I(X_i=k)}-I\left(\sum_{i=1}^{n} I(X_i=k) = 0\right). \end{align}\] Hence the claim is verified and \(a_n\) is non-increasing. Apply Lemma 6 we can now assume \(n\widehat{p}_k \sim \text{Poisson}(np_k)\) and \(n\widehat{p}_1,\dots, n\widehat{p}_d\) are independent (with an additional factor 2) \[\begin{align} &\, \mathbb{P}\left(\sum_{j=1}^d \widehat{t}_j=0\right) \\ \leq &\, 2 \mathbb{E}\left[\prod_{k=1}^d\left\{\left(1-\pi_k\right)^{n \widehat{p}_k}+\pi_k^{n \widehat{p}_k}-I\left(\widehat{p}_k=0\right)\right\}\right] \\ =&\, 2 \prod_{k=1}^d \mathbb{E}\left[\left(1-\pi_k\right)^{n \widehat{p}_k}+\pi_k^{n \widehat{p}_k}-I\left(\widehat{p}_k=0\right)\right] \\ =&\, 2 \prod_{k=1}^d \left[ \exp(-n\pi_k p_k) + \exp(-n(1-\pi_k)p_k) -\exp(-np_k) \right] \\ \leq &\, 2 \prod_{k=1}^d \left[ \exp(-\epsilon n p_k) + \exp(-(1-\epsilon)n p_k) -\exp(-np_k) \right]. \end{align}\] The first equation follows from independence and second one follows from probability generating function of Poisson distribution and the last inequality follows since \(f_3(x) = \exp(-cx) + \exp(-c(1-x))\) (for constant \(c>0\)) is decreasing on \([0,1/2]\) and increasing on \([1/2, 1]\). For simplicity define \[Z_k = \exp(-\epsilon n p_k) + \exp(-(1-\epsilon)n p_k) -\exp(-np_k) .\] 1. In the first case \(\epsilon n p_k \geq 2\log 2\) we have (recall \(0 < \epsilon < 1/2\)) \[Z_k \leq 2\exp(-\epsilon n p_k) \leq \exp\left( -\frac{\epsilon n p_k}{2} \right).\] 2. In the second case \(n p_k \leq \epsilon\), we use the following inequalities for \(t \geq 0\): \[\begin{align} & \exp (-t) \leq 1-t+t^2 / 2, \\ & \exp (-t) \geq 1-t+t^2 / 2-t^3 / 6 \end{align}\] and obtain \[\begin{align} Z_k \leq &\, 1-\epsilon n p_k + \frac{\epsilon^2 n^2 p_k^2}{2} + 1-(1-\epsilon)n p_k + \frac{(1-\epsilon)^2 n^2 p_k^2}{2} \\ &\, - \left( 1-n p_k + \frac{n^2 p_k^2}{2} - \frac{n^3 p_k^3}{6} \right) \\ = &\, 1-\epsilon(1-\epsilon)n^2 p_k^2 + \frac{n^3 p_k^3}{6} \\ \leq & \, 1- \left( \frac{5}{6}\epsilon -\epsilon^2 \right)n^2 p_k^2\\ \leq & \, 1-\frac{\epsilon n^2 p_k^2}{3} \\ \leq &\, \exp \left( -\frac{\epsilon n^2 p_k^2}{3} \right). \end{align}\] 3. In the last case \(\epsilon < n p_k < \frac{2 \log 2}{\epsilon}\), since \(f_4(x) = \exp(-\epsilon x) + \exp(-(1-\epsilon)x) -\exp(-x)\) is monotonely non-increasing on \([0,+\infty]\) (one can check this by taking the first-order derivative easily), we have \[Z_k \leq \exp(-\epsilon^2) + \exp(-\epsilon(1-\epsilon)) - \exp(-\epsilon) = \exp (-C_1(\epsilon))\] where \(C_1(\epsilon) = -\log(\exp(-\epsilon^2) + \exp(-\epsilon(1-\epsilon)) - \exp(-\epsilon)) \in(0, \epsilon^2)\). Note that \[1 > \frac{\epsilon n p_k}{2 \log 2}.\] We have \[\exp (-C_1(\epsilon)) \leq \exp (-C_2(\epsilon) n p_k)\] where \[C_2(\epsilon) = \frac{C_1(\epsilon)\epsilon}{2 \log 2} < \frac{\epsilon ^3}{2 \log 2} < \frac{\epsilon }{2}\] Hence we can combine the first case with the third case as \[Z_k \leq \exp (-C_2(\epsilon) n p_k).\] Let \(I_1\) include the index \(k \in [d]\) in the first and third case, \(I_2\) include the index \(k\) in the second case. Denote \(S_1 = \sum_{i\in I_1} p_i, T_2 = \sum_{i \in I_2}p_i^2\). We thus have \[\prod_{k=1}^d Z_k \leq \exp(-C_2(\epsilon)nS_1) \exp \left(-\frac{\epsilon n^2 T_2}{3} \right).\] In the case \(n \geq d\) we have \[1-S_1 = \sum_{i \in I_2} p_i \leq \frac{\epsilon |I_2|}{n} \leq \frac{\epsilon d}{n} \leq \epsilon,\] i.e. \(S_1 > 1-\epsilon\) and we have \[\prod_{k=1}^d Z_k \leq \exp(-C_2(\epsilon)(1-\epsilon)n) \leq \exp\left(- \frac{C_2(\epsilon)n}{2 } \right).\] In the case \(n < d\), if \(S_1 \geq 1/2\) then we also have \[\prod_{k=1}^d Z_k \leq \exp\left(- \frac{C_2(\epsilon)n}{2 } \right).\] In the case \(S_1 < 1/2\), then by Cauchy-Schwarz inequality we have \[\frac{1}{4} \leq (1-S_1)^2 = \left(\sum_{i \in I_2} p_i \right)^2 \leq |I_2|T_2 \leq d T_2,\] i.e. \(T_2 \geq \frac{1}{4d}\), this further implies \[\prod_{k=1}^d Z_k \leq \exp \left(- \frac{\epsilon n^2}{12d} \right).\] So we conclude that \[\prod_{k=1}^d Z_k \leq \max \left\{ \exp\left(- \frac{C_2(\epsilon)n}{2 } \right), \exp \left(- \frac{\epsilon n^2}{12d} \right) \right \} \leq \exp\left(-C (\epsilon) \frac{n^2}{n \vee d} \right)\] where \[C(\epsilon) = \min \left( \frac{C_2(\epsilon)}{2}, \frac{\epsilon}{12} \right).\] ◻

11.7 Proof of Theorem 3 ↩︎

Proof. We first bound the bias. By conditioning on \(\mathbf{X}^n,\mathbf{A}^n\) and noting \(\mathbb{E}[\widehat{\mu}_{1k} \mid \mathbf{X}^n, \mathbf{A}^n] = \mu_{1k} I(\widehat{p}_k \widehat{\pi}_k >0)\), \(\mathbb{E}[\widehat{\mu}_{0k} \mid \mathbf{X}^n, \mathbf{A}^n] = \mu_{0k} I(\widehat{p}_k (1-\widehat{\pi}_k) >0)\) we have \[\mathbb{E}\left[\widehat{\tau} \mid\mathbf{X}^n, \mathbf{A}^n\right] = \frac{\sum_{k=1}^d \widehat{t}_k \tau_k}{\sum_{k=1}^d \widehat{t}_k} I\left(\sum_{j=1}^d \widehat{t}_j>0\right).\] Note that \[\psi = \psi I\left(\sum_{j=1}^d \widehat{t}_j>0\right) + \psi I\left(\sum_{j=1}^d \widehat{t}_j=0\right) = \frac{\sum_{k=1}^d \widehat{t}_k \psi}{\sum_{k=1}^d \widehat{t}_k} I\left(\sum_{j=1}^d \widehat{t}_j>0\right) + \psi I\left(\sum_{j=1}^d \widehat{t}_j=0\right).\] Hence the bias can be expressed as \[\begin{align} &\,\mathbb{E}\left[\widehat{\tau}-\psi\right] \\ =&\, \mathbb{E}\left[ \frac{\sum_{k=1}^d \widehat{t}_k (\tau_k-\psi)}{\sum_{k=1}^d \widehat{t}_k} I\left(\sum_{j=1}^d \widehat{t}_j>0\right) \right] - \psi \mathbb{P}\left(\sum_{j=1}^d \widehat{t}_j=0\right) \end{align}\] Using \(\sigma_n\) in 12 to parameterize the effect homogeneity and Lemma 1, we have \[\begin{align} &\,\left |\mathbb{E}\left[\widehat{\tau}-\psi\right] \right| \\ \leq &\, \mathbb{E}\left[ \frac{\sum_{k=1}^d \widehat{t}_k |\tau_k-\psi|}{\sum_{k=1}^d \widehat{t}_k} I\left(\sum_{j=1}^d \widehat{t}_j>0\right) \right]+ |\psi| \mathbb{P}\left(\sum_{j=1}^d \widehat{t}_j=0\right) \\ \leq &\, \sigma_n + 2 \exp\left(-C (\epsilon) \frac{n^2}{n \vee d} \right). \end{align}\] We then consider the variance. By the property of conditional variance, we have \[\label{eq:var-homogeneity} \operatorname{Var}(\widehat{\tau}) = \operatorname{Var}\left(\mathbb{E}\left[\widehat{\tau} \mid\mathbf{X}^n, \mathbf{A}^n\right] \right) + \mathbb{E}\left[ \operatorname{Var}\left(\widehat{\tau} \mid \mathbf{X}^n, \mathbf{A}^n\right) \right].\tag{24}\] Rewrite \[\mathbb{E}\left[\widehat{\tau} \mid\mathbf{X}^n, \mathbf{A}^n\right] = \frac{\sum_{k=1}^d \widehat{t}_k \tau_k}{\sum_{k=1}^d \widehat{t}_k} I\left(\sum_{j=1}^d \widehat{t}_j>0\right) = \frac{\sum_{k=1}^d \widehat{t}_k (\tau_k-\psi)}{\sum_{k=1}^d \widehat{t}_k} I\left(\sum_{j=1}^d \widehat{t}_j>0\right) + \psi I\left(\sum_{j=1}^d \widehat{t}_j>0\right).\] For a bounded random variable \(X\) satisfying \(m \leq X \leq M\) we have \(\operatorname{Var}(X) \leq (M-m)^2/4\). Note that \[\left| \frac{\sum_{k=1}^d \widehat{t}_k (\tau_k-\psi)}{\sum_{k=1}^d \widehat{t}_k} I\left(\sum_{j=1}^d \widehat{t}_j>0\right) \right| \leq \sigma_n\] and we have \[\operatorname{Var}\left(\frac{\sum_{k=1}^d \widehat{t}_k (\tau_k-\psi)}{\sum_{k=1}^d \widehat{t}_k} I\left(\sum_{j=1}^d \widehat{t}_j>0\right) \right) \leq \sigma_n^2.\] Further notice \[\operatorname{Var}\left( \psi I\left(\sum_{j=1}^d \widehat{t}_j>0\right) \right) \leq \psi^2 \mathbb{P}\left(\sum_{j=1}^d \widehat{t}_j=0 \right) \leq 2 \exp\left(-C (\epsilon) \frac{n^2}{n \vee d} \right).\] The covariance can be expressed as \[\begin{align} &\, \operatorname{Cov} \left( \frac{\sum_{k=1}^d \widehat{t}_k (\tau_k-\psi)}{\sum_{k=1}^d \widehat{t}_k} I\left(\sum_{j=1}^d \widehat{t}_j>0\right), \psi I\left(\sum_{j=1}^d \widehat{t}_j>0\right)\right) \\ = &\, \psi \mathbb{E}\left[\frac{\sum_{k=1}^d \widehat{t}_k (\tau_k-\psi)}{\sum_{k=1}^d \widehat{t}_k} I\left(\sum_{j=1}^d \widehat{t}_j>0\right) \right] -\psi \mathbb{E}\left[\frac{\sum_{k=1}^d \widehat{t}_k (\tau_k-\psi)}{\sum_{k=1}^d \widehat{t}_k} I\left(\sum_{j=1}^d \widehat{t}_j>0\right) \right] \mathbb{P}\left(\sum_{j=1}^d \widehat{t}_j>0 \right) \\ =&\, \psi \mathbb{E}\left[\frac{\sum_{k=1}^d \widehat{t}_k (\tau_k-\psi)}{\sum_{k=1}^d \widehat{t}_k} I\left(\sum_{j=1}^d \widehat{t}_j>0\right) \right] \mathbb{P}\left(\sum_{j=1}^d \widehat{t}_j=0 \right) \end{align}\] and hence we can bound the covariance as follows: \[\begin{align} &\, \left|\operatorname{Cov} \left( \frac{\sum_{k=1}^d \widehat{t}_k (\tau_k-\psi)}{\sum_{k=1}^d \widehat{t}_k} I\left(\sum_{j=1}^d \widehat{t}_j>0\right), \psi I\left(\sum_{j=1}^d \widehat{t}_j>0\right)\right) \right|\\ \leq &\, 2\sigma_n \exp\left(-C (\epsilon) \frac{n^2}{n \vee d} \right). \end{align}\] Combining these three terms above we have \[\label{eq:var-cond-exp} \operatorname{Var}\left(\mathbb{E}\left[\widehat{\tau}\mid\mathbf{X}^n, \mathbf{A}^n\right] \right) \leq \sigma_n^2 + 2 \exp\left(-C (\epsilon) \frac{n^2}{n \vee d} \right) + 4\sigma_n \exp\left(-C (\epsilon) \frac{n^2}{n \vee d} \right).\tag{25}\] To complete the proof we need to bound \[\mathbb{E}\left[ \operatorname{Var}\left(\widehat{\tau} \mid \mathbf{X}^n, \mathbf{A}^n\right) \right].\] Note that conditioning on \(\mathbf{X}^n,\mathbf{A}^n\) the estimated \(\{\widehat{\tau}_k, 1\leq k \leq d \}\) are independent with \[\operatorname{Var}(\widehat{\tau}_k\mid\mathbf{X}^n, \mathbf{A}^n) = \frac{\mu_{1k}(1-\mu_{1k})}{n\widehat{p}_k \widehat{\pi}_k} I(\widehat{p}_k \widehat{\pi}_k>0) + \frac{\mu_{0k}(1-\mu_{0k})}{n\widehat{p}_k (1-\widehat{\pi}_k)}I\left(\widehat{p}_k (1-\widehat{\pi}_k)>0\right).\] We have \[\begin{align} \operatorname{Var}(\widehat{\tau}\mid\mathbf{X}^n, \mathbf{A}^n) = &\, \frac{I\left(\sum_{j=1}^d \widehat{t}_j > 0\right)}{\left(\sum_{j=1}^d \widehat{t}_j\right)^2} \sum_{k=1}^d \widehat{p}_k^2 I(0 < \widehat{\pi}_k < 1) \left( \frac{\mu_{1k}(1-\mu_{1k})}{n\widehat{p}_k \widehat{\pi}_k} +\frac{\mu_{0k}(1-\mu_{0k})}{n\widehat{p}_k (1-\widehat{\pi}_k)} \right)\\ \leq &\, \frac{1}{4}\frac{I\left(\sum_{j=1}^d \widehat{t}_j > 0\right)}{\left(\sum_{j=1}^d \widehat{t}_j\right)^2} \sum_{k=1}^d \widehat{p}_k^2 I(0 < \widehat{\pi}_k < 1) \left( \frac{1}{n\widehat{p}_k \widehat{\pi}_k} +\frac{1}{n\widehat{p}_k (1-\widehat{\pi}_k)} \right) \end{align}\] Denote \(\mathbf{I}= (I_1,\dots, I_d)\) with \(I_k = I(0 < \widehat{\pi}_k < 1)\) and further note that \[\begin{align} &\,\frac{1}{4} \mathbb{E}\left[ \frac{I\left(\sum_{j=1}^d \widehat{t}_j > 0\right)}{\left(\sum_{j=1}^d \widehat{t}_j\right)^2} \sum_{k=1}^d \frac{\widehat{p}_k^2 I(0 < \widehat{\pi}_k < 1)}{n\widehat{p}_k \widehat{\pi}_k} \right]\\ = &\, \frac{1}{4} \mathbb{E}\left[ \frac{I\left(\sum_{j=1}^d \widehat{t}_j > 0\right)}{\left(\sum_{j=1}^d \widehat{t}_j\right)^2} \sum_{k=1}^d \widehat{p}_k^2 I(0 < \widehat{\pi}_k < 1) \mathbb{E}\left(\frac{1}{n\widehat{p}_k \widehat{\pi}_k} \mid \mathbf{X}^n, \mathbf{I}\right) \right] \\ = &\, \frac{1}{4} \mathbb{E}\left[ \frac{I\left(\sum_{j=1}^d \widehat{t}_j > 0\right)}{\left(\sum_{j=1}^d \widehat{t}_j\right)^2} \sum_{k=1}^d \widehat{p}_k^2 I(0 < \widehat{\pi}_k < 1) \mathbb{E}\left(\frac{1}{n\widehat{p}_k \widehat{\pi}_k} \mid \mathbf{X}^n, I_k \right) \right], \end{align}\] where the last equation follows from the independence of treatment assignment within each category after covariates \(\mathbf{X}^n\) are sampled. It is easy to see from Lemma 2 that for a Binomial random variable \(V \sim \text{B}(n,p)\) we have \[\mathbb{E}\left[\frac{1}{V} \mid 0 < V < n\right] \leq \frac{2}{(n+1)p\left[1-p^n-(1-p)^n\right]}.\] Thus we have (assume \(n \geq 2\)) \[\begin{align} &\,\widehat{p}_k^2 I(0 < \widehat{\pi}_k < 1) \mathbb{E}\left(\frac{1}{n\widehat{p}_k \widehat{\pi}_k} \mid \mathbf{X}^n, I_k \right)\\ \leq &\, \frac{2\widehat{p}_k^2 I(0 < \widehat{\pi}_k < 1)}{(n\widehat{p}_k+1)\pi_k\left[1-\pi_k^n-(1-\pi_k)^n\right]}\\ \leq &\, \frac{2\widehat{p}_k^2 I(0 < \widehat{\pi}_k < 1)}{(n\widehat{p}_k+1)\epsilon\left[1-\epsilon^n-(1-\epsilon)^n\right]}\\ \leq &\, \frac{2\widehat{p}_k^2 I(0 < \widehat{\pi}_k < 1)}{(n\widehat{p}_k+1)\epsilon\left[1-\epsilon^2-(1-\epsilon)^2\right]}\\ = &\, \frac{\widehat{p}_k^2 I(0 < \widehat{\pi}_k < 1)}{(n\widehat{p}_k+1)\epsilon^2(1-\epsilon)} \\ \leq &\, \frac{\widehat{p}_k I(0 < \widehat{\pi}_k < 1)}{\epsilon^2(1-\epsilon)n}. \end{align}\] Sum these terms up, we have \[\frac{1}{4} \mathbb{E}\left[ \frac{I\left(\sum_{j=1}^d \widehat{t}_j > 0\right)}{\left(\sum_{j=1}^d \widehat{t}_j\right)^2} \sum_{k=1}^d \widehat{p}_k^2 I(0 < \widehat{\pi}_k < 1) \mathbb{E}\left(\frac{1}{n\widehat{p}_k \widehat{\pi}_k} \mid \mathbf{X}^n, I_k \right) \right] \leq \frac{1}{4\epsilon^2(1-\epsilon)} \mathbb{E}\left[ \frac{I\left(\sum_{j=1}^d \widehat{t}_j > 0\right)}{n\sum_{j=1}^d \widehat{t}_j} \right]\] Similarly, we can show \[\frac{1}{4} \mathbb{E}\left[ \frac{I\left(\sum_{j=1}^d \widehat{t}_j > 0\right)}{\left(\sum_{j=1}^d \widehat{t}_j\right)^2} \sum_{k=1}^d \frac{\widehat{p}_k^2 I(0 < \widehat{\pi}_k < 1)}{n\widehat{p}_k (1-\widehat{\pi}_k)} \right]\leq \frac{1}{4\epsilon(1-\epsilon)^2} \mathbb{E}\left[ \frac{I\left(\sum_{j=1}^d \widehat{t}_j > 0\right)}{n\sum_{j=1}^d \widehat{t}_j} \right]\] So we have the following bound on the expected conditional variance: \[\mathbb{E}\left[\operatorname{Var}(\widehat{\tau}\mid\mathbf{X}^n, \mathbf{A}^n)\right]\leq \frac{1}{4\epsilon^2(1-\epsilon)^2} \mathbb{E}\left[ \frac{I\left(\sum_{j=1}^d \widehat{t}_j > 0\right)}{n\sum_{j=1}^d \widehat{t}_j} \right]\] Define \(h(\widehat{{\boldsymbol{p}}}, \widehat{\pi})\) \[h(\widehat{{\boldsymbol{p}}}, \widehat{\boldsymbol{\pi}}) = \begin{cases} 1 &\quad\text{if} \sum_{j=1}^d \widehat{t}_j =0\\ \frac{1}{\sum_j n\widehat{p}_j I(0 < \widehat{\pi}_j < 1)} &\quad\text{if} \sum_{j=1}^d \widehat{t}_j >0 . \\ \end{cases}\] Note that \(\sum_j n\widehat{p}_j I(0 < \widehat{\pi}_j < 1)\) is the number of subjects in categories with both treated and untreated units, thus it is an integer and will not decrease as we collect more samples. As a result, \(\mathbb{E}[h(\widehat{{\boldsymbol{p}}}, \widehat{\boldsymbol{\pi}})]\) is non-increasing in \(n\) and by Lemma 6 we have \[\label{eq:bound-the-inverse} \begin{align} &\, \mathbb{E}\left[ \frac{I\left(\sum_{j=1}^d \widehat{t}_j > 0\right)}{n\sum_{j=1}^d \widehat{t}_j} \right]\\ \leq &\, \mathbb{E}[h(\widehat{{\boldsymbol{p}}}, \widehat{\boldsymbol{\pi}})]\\ \leq &\,2\mathbb{E}_{\text{Poi}}[h(\widehat{{\boldsymbol{p}}}, \widehat{\boldsymbol{\pi}})] \\ =&\, 2\mathbb{P}\left(\sum_j n\widehat{p}_j I(0 < \widehat{\pi}_j < 1) = 0\right) + 2 \mathbb{E}_{\text{Poi}} \left[ \frac{I\left(\sum_j n\widehat{p}_j I(0 < \widehat{\pi}_j < 1)\geq 1\right)}{\sum_j n\widehat{p}_j I(0 < \widehat{\pi}_j < 1)}\right]\\ \leq &\, 2 \exp\left(-C(\epsilon)\frac{n^2}{n \vee d} \right) + 4\mathbb{E}_{\text{Poi}} \left[ \frac{1}{1+\sum_j n\widehat{p}_j I(0 < \widehat{\pi}_j < 1)} \right], \end{align}\tag{26}\] where we use the fact that \(\sum_j n\widehat{p}_j I(0 < \widehat{\pi}_j < 1)\) is an integer and the last inequality follows from Lemma 1 and the inequality \[\frac{I(x \geq 1)}{x} \leq \frac{2}{1+x}, \, x\geq 0.\] \(\mathbb{E}_{\text{Poi}}\) means the components of \((n\widehat{p}_1, \dots, n\widehat{p}_d)\) are independent and \(n\widehat{p}_k \sim \text{Poi}(np_k)\). Let \(W = \sum_j n\widehat{p}_j I(0 < \widehat{\pi}_j < 1)\) and we will derive a tail bound for \(1/(1+W)\) by bounding its MGF \(\mathbb{E}_{\text{Poi}}[\exp(-c_1W)]\) for an absolute constant \(c_1>0\) (e.g. can be taken as \(1/2\)). Note that given \(\widehat{p}_k\), \(n\widehat{p}_k \widehat{\pi}_k\) is a binomial variable and \[\mathbb{P}(0 < \widehat{\pi}_k < 1 \mid \widehat{p}_k) = 1-\pi_k^{n\widehat{p}_k} - (1-\pi_k)^{n\widehat{p}_k}\] when \(\widehat{p}_k>0\). Denote \(q_k(j):= 1-\pi_k^{j} - (1-\pi_k)^{j}\). We have \[\begin{align} &\,\mathbb{E}_{\text{Poi}}[\exp(-c_1n \widehat{p}_k I(0<\widehat{\pi}_k <1))] \\ \leq &\, \exp(-np_k)(1+np_k) +\sum_{j=2}^\infty \mathbb{E}_{\text{Poi}}[\exp(-c_1n \widehat{p}_k I(0<\widehat{\pi}_k <1))\mid n\widehat{p}_k = j] \exp(-np_k) \frac{(np_k)^j}{j!}\\ =&\, \exp(-np_k)(1+np_k) +\sum_{j=2}^\infty \left[ \exp(-c_1j)q_k(j)+1-q_k(j) \right] \exp(-np_k) \frac{(np_k)^j}{j!}. \end{align}\] For \(c_2>0\) another constant to be fixed, we divide into two cases:

Case 1: \(np_k \leq c_2\), we have \[\begin{align} &\,\mathbb{E}_{\text{Poi}}[\exp(-c_1n \widehat{p}_k I(0<\widehat{\pi}_k <1))]\\ \leq &\, 1-\sum_{j=2}^\infty q_k(j)(1-\exp(-c_1j))\exp(-np_k)\frac{(np_k)^j}{j!}\\ \leq &\, 1-q_k(2)(1-\exp(-2c_1))\exp(-np_k)\frac{(np_k)^2}{2} \\ \leq &\, 1-\epsilon (1-\epsilon)(1-\exp(-2c_1))\exp(-c_2){(np_k)^2} \\ \leq &\, \exp(-c_3 n^2p_k^2), \end{align}\] where \(c_3 = \epsilon (1-\epsilon)(1-\exp(-2c_1))\exp(-c_2)\).

Case 2: \(np_k > c_2\), we have \[\begin{align} &\,\mathbb{E}_{\text{Poi}}[\exp(-c_1n \widehat{p}_k I(0<\widehat{\pi}_k <1))]\\ \leq &\, \exp(-np_k)(1+np_k) +\sum_{j=2}^\infty \exp(-c_1j)\exp(-np_k) \frac{(np_k)^j}{j!} \\ &\,+\sum_{j=2}^\infty (1-q_k(j))(1-\exp(-c_1j)) \exp(-np_k) \frac{(np_k)^j}{j!}\\ \end{align}\] The summation of first two terms is equal to \[\begin{align} &\,\exp(-np_k)(1+np_k) +\exp(-np_k) \left [\exp\{\exp(-c_1)np_k\}-1-\exp(-c_1)np_k \right]\\ =&\, \exp\{-(1-\exp(-c_1))np_k\} +np_k \exp(-np_k)(1-\exp(-c_1)). \end{align}\] The third term can be bounded as \[\begin{align} &\,\sum_{j=2}^\infty (1-q_k(j))(1-\exp(-c_1j)) \exp(-np_k) \frac{(np_k)^j}{j!}\\ \leq &\, \sum_{j=2}^\infty (1-q_k(j)) \exp(-np_k) \frac{(np_k)^j}{j!}\\ = &\, \sum_{j=2}^\infty (\pi_k^j + (1-\pi_k)^j) \exp(-np_k) \frac{(np_k)^j}{j!} \\ \leq &\,\sum_{j=2}^\infty (\epsilon^j + (1-\epsilon)^j) \exp(-np_k) \frac{(np_k)^j}{j!} \\ \leq &\,\sum_{j=0}^\infty (\epsilon^j + (1-\epsilon)^j) \exp(-np_k) \frac{(np_k)^j}{j!} \\ =&\, \exp(-(1-\epsilon)np_k) + \exp(-\epsilon np_k). \end{align}\] Hence we have \[\begin{align} &\,\mathbb{E}_{\text{Poi}}[\exp(-c_1n \widehat{p}_k I(0<\widehat{\pi}_k <1))]\\ \leq &\, \exp\{-(1-\exp(-c_1))np_k\} +np_k \exp(-np_k)(1-\exp(-c_1)) + \exp(-(1-\epsilon)np_k) + \exp(-\epsilon np_k)\\ \leq &\, \exp(-c_4np_k) \end{align}\] where we take \(c_2>0\) sufficiently large and \(c_4>0\) sufficiently small (both depend on \(\epsilon\)) such that \[\exp\{-(1-\exp(-c_1))x\} +x \exp(-x)(1-\exp(-c_1)) + \exp(-(1-\epsilon)x) + \exp(-\epsilon x)\leq \exp(-c_4x) , \, x\geq c_2.\] Denote the indices corresponding to case 1 as \(I_1\) and case 2 as \(I_2\). Now we have \[\begin{align} &\,\mathbb{E}_{\text{Poi}}[\exp(-c_1W)]\\ = &\, \prod_{k=1}^d \mathbb{E}_{\text{Poi}}[\exp(-c_1n\widehat{p}_k I(0 < \widehat{\pi}_k < 1))]\\ \leq &\, \exp \left\{-\left(c_3 \sum_{k \in I_1}n^2p_k^2 + c_4 \sum_{k \in I_2}np_k \right)\right\}\\ = &\, \exp \left\{-\left(c_3 n^2T_1 + c_4 nS_2 \right)\right\}, \end{align}\] where \(T_1 = \sum_{k \in I_1}p_k^2, S_2 = \sum_{k \in I_2}p_k\). In the case \(d \leq n/(2c_2)\), we have \[1-S_2 = \sum_{k \in I_1}p_k \leq \frac{c_2|I_1|}{n} \leq \frac{c_2d}{n} \leq \frac{1}{2}.\] So \(S_2 \geq 1/2\) and we have \[\mathbb{E}_{\text{Poi}}[\exp(-c_1W)] \leq \exp(-c_4n/2)\] In the case \(d > n/(2c_2)\), if \(S_2 \geq 1/2\), the above inequality still holds. If \(S_2 < 1/2\), by Cauchy-Schwarz inequality we have \[\frac{1}{4} \leq (1-S_2)^2=\left(\sum_{k \in I_1}p_k \right)^2 \leq |I_1|T_1 \leq dT_1.\] Hence \(T_1 \geq \frac{1}{4d}\) and \[\mathbb{E}_{\text{Poi}}[\exp(-c_1W)] \leq \exp\left(-c_3\frac{n^2}{4d} \right)\] Thus we always have \[\mathbb{E}_{\text{Poi}}[\exp(-c_1W)] \leq \max\left\{\exp(-c_4n/2), \exp\left(-c_3\frac{n^2}{4d} \right)\right\} \leq \exp \left(-c_5 \frac{n^2}{n \vee d}\right).\] Finally for a small constant \(c_6>0\) such that \(c_1c_6 < c_5\) we have \[\begin{align} &\,\mathbb{P}_{\text{Poi}}\left(W\leq \frac{c_6n^2}{n \vee d} \right)\\ \leq &\, \exp\left(\frac{c_1c_6n^2}{n\vee d}\right) \mathbb{E}_{\text{Poi}}[\exp(-c_1W)] \\ \leq &\, \exp\left(\frac{c_1c_6n^2}{n\vee d}-\frac{c_5n^2}{n\vee d}\right) \\ \leq &\, \exp\left(\frac{-c_7n^2}{n\vee d}\right) \end{align}\] Hence we have the following bound \[\mathbb{E}_{\text{Poi}}\left[\frac{1}{1+W} \right] \leq \frac{1}{c_6} \frac{n \vee d}{n^2} + \mathbb{P}_{\text{Poi}}\left(W\leq \frac{c_6n^2}{n \vee d}\right) \\ \leq \frac{1}{c_6} \frac{n \vee d}{n^2} + \exp\left(\frac{-c_7n^2}{n\vee d}\right).\] Plug into 26 , we conclude \[\mathbb{E}\left[ \frac{I\left(\sum_{j=1}^d \widehat{t}_j > 0\right)}{n\sum_{j=1}^d \widehat{t}_j} \right] \lesssim \exp\left(-C'(\epsilon)\frac{n^2}{n\vee d}\right) + \frac{n \vee d}{n^2},\] which is also a bound on expected conditional variance. Combining this bound with 25 we have \[\operatorname{Var}(\widehat{\tau}) \lesssim \sigma_n^2 + \exp\left(-C'(\epsilon)\frac{n^2}{n\vee d}\right) + \frac{n \vee d}{n^2}.\] Thus the mean squared error of \(\widehat{\tau}\) can be bounded as \[\mathbb{E}[(\widehat{\tau}-\psi)^2] \lesssim \sigma_n^2 + \exp\left(-C'(\epsilon)\frac{n^2}{n\vee d}\right) + \frac{n \vee d}{n^2}.\] ◻

11.8 Proof of Theorem 4 ↩︎

Proof. The parametric rate \(1/n\) can be similarly achieved by using the construction in proof of Theorem 2. We focus on proving the \(d/n^2\) lower bound when \(d \lesssim n^2\).

By Lemma 4, we can again focus on a Poissonized experiment where the total sample size \(N \sim \mathrm{Poi}(n)\). Since \(p_k = 1/d\), this is equivalent to having independent category counts \[N_k \sim \mathrm{Poi}(\lambda), \qquad \lambda := \frac{n}{d}, \qquad 1 \le k \le d.\] We then define the null and alternative distributions. Fix \(q \in (0,q_0]\), where \[q_0 := \min\left\{\frac{1}{8},\, 8\Bigl(\frac{1}{2}-\epsilon\Bigr)^2\right\}, \qquad r := \sqrt{\frac{q}{8}}.\] Define the null model \(P_0 \in \mathcal{H}_d(\epsilon)\) by \[p_k = \frac{1}{d}, \qquad \pi_k = \frac{1}{2}, \qquad \mu_{0k} = \mu_{1k} = \frac{1}{2}, \qquad 1 \le k \le d.\] Then \(\psi(P_0)=0\). For the alternative, let \(\xi_1,\dots,\xi_d\) be i.i.d.Rademacher random variables. Given \(\boldsymbol{\xi}=(\xi_1,\dots,\xi_d)\), define the model \(P_\xi\) by \[p_k = \frac{1}{d}, \qquad \pi_k = \frac{1}{2} + r \xi_k,\] \[\mu_{0k} = \frac{1}{2} - r \xi_k - \frac{q}{4}, \qquad \mu_{1k} = \frac{1}{2} - r \xi_k + \frac{q}{4}, \qquad 1 \le k \le d.\] Since \(r \le \frac{1}{2}-\epsilon\) and \(r + q/4 \le 1/2\), every realization \(P_\xi\) belongs to \(\mathcal{H}_d(\epsilon)\). Moreover, \[\mu_{1k} - \mu_{0k} = \frac{q}{2} \qquad \text{for all } k \in [d],\] so that \[\psi(P_\xi) = \frac{q}{2} \qquad \text{almost surely}.\]

With this construction, we will see that one observation from a category is completely uninformative on average. Order the four cells as \[(A,Y)\in\{(1,1),(1,0),(0,1),(0,0)\},\] and write \[u_0 := \left(\frac{1}{4},\frac{1}{4},\frac{1}{4},\frac{1}{4}\right).\] Under the null model \(P_0\), the cell probabilities are exactly \(u_0\).

Under the alternative \(P_\xi\), for category \(k\) a direct calculation gives the following cell probabilities \[u_{\xi_k} = \left( \frac{1}{4}+2r^3\xi_k,\; \frac{1}{4}+(r-2r^3)\xi_k,\; \frac{1}{4}+(-r+2r^3)\xi_k,\; \frac{1}{4}-2r^3\xi_k \right).\] Equivalently, \[u_{\xi_k} = u_0 + \xi_k \Delta, \qquad \Delta := \bigl(2r^3,\;r-2r^3,\;-r+2r^3,\;-2r^3\bigr).\] Hence \[\frac{1}{2}u_{+1} + \frac{1}{2}u_{-1} = u_0.\] Therefore, after averaging over the latent sign \(\xi_k\), the law of one observation from category \(k\) is exactly the same under the alternative prior and the null. Thus categories observed zero or one time contain no information for distinguishing the null from the alternative; only repeated observations within a category can separate them.

We then compute the \(\chi^2\) divergence between the null and marginalized alternative \(P_1 = \mathbb{E}[P_{\xi}]\). Under the null \(P_0\), the four cell counts in category \(k\), \[\mathbf{N}_k = (N_{k,11},N_{k,10},N_{k,01},N_{k,00}),\] are independent with law \[Q_0 := \bigotimes_{j=1}^4 \mathrm{Poi}(\lambda/4).\] Under the alternative prior, conditional on \(\xi_k=\pm 1\), the category law is \[Q_\pm := \bigotimes_{j=1}^4 \mathrm{Poi}\!\bigl(\lambda(u_{0j}\pm \Delta_j)\bigr),\] and after averaging over \(\xi_k\), \[Q_1 := \frac{1}{2}Q_+ + \frac{1}{2}Q_-.\] Let \[s_q := \|\Delta\|_2^2.\] A direct computation gives \[s_q = (2r^3)^2 + (r-2r^3)^2 + (-r+2r^3)^2 + (-2r^3)^2 = \frac{q}{4} - \frac{q^2}{8} + \frac{q^3}{32} \le \frac{q}{4}.\] We now compute \(\chi^2(Q_1 \|Q_0)\). Since \[u_{0j}\pm \Delta_j = \frac{1}{4}(1\pm 4\Delta_j),\] the likelihood ratios satisfy \[\frac{dQ_\pm}{dQ_0}(\mathbf{N}_k) = \prod_{j=1}^4 (1\pm 4\Delta_j)^{N_{k,j}},\] because \(\sum_{j=1}^4 \Delta_j = 0\). Therefore \[\mathbb{E}_{Q_0}\!\left[\left(\frac{dQ_{\pm}}{dQ_0}\right)^2\right] = \prod_{j=1}^4 \exp\!\left( \frac{\lambda}{4}\bigl[(1\pm 4\Delta_j)^2 - 1\bigr] \right) = \exp(4\lambda s_q).\] Similarly, \[\mathbb{E}_{Q_0}\!\left[ \frac{dQ_+}{dQ_0}\frac{dQ_-}{dQ_0} \right] = \prod_{j=1}^4 \exp\!\left( \frac{\lambda}{4}\bigl[(1+4\Delta_j)(1-4\Delta_j)-1\bigr] \right) = \exp(-4\lambda s_q).\] Hence \[1+\chi^2(Q_1\|Q_0) = \mathbb{E}_{Q_0}\!\left[\left(\frac{dQ_1}{dQ_0}\right)^2\right] = \frac{1}{4} \left( e^{4\lambda s_q}+2e^{-4\lambda s_q}+e^{4\lambda s_q} \right) = \cosh(4\lambda s_q).\] Since \(4\lambda s_q \le \lambda q\), and later we will choose \(q\) so that \(\lambda q \le 1\), it follows that \[\chi^2(Q_1\|Q_0) = \cosh(4\lambda s_q)-1 \leq C' \lambda^2 s_q^2 \leq C(\lambda q)^2\] for some universal constant \(C>0\).

Because the categories are independent under Poissonization and the latent signs \(\xi_k\) are independent, the cell counts’ laws satisfy \[{P}_0 = Q_0^{\otimes d}. \qquad {P}_1 = Q_1^{\otimes d},\] Therefore \[1+\chi^2({P}_1\|{P}_0) = \bigl(1+\chi^2(Q_1\|Q_0)\bigr)^d \le \exp\!\bigl(Cd\lambda^2 q^2\bigr) = \exp\!\left(C\frac{n^2q^2}{d}\right).\] Now choose \[q^2 = c_0\frac{d}{n^2},\] where \(c_0>0\) is a sufficiently small constant depending only on \(\epsilon\), chosen so that \(q\le q_0\) and \(\lambda q \le 1\). Then the \(\chi^2\) divergence is bounded away from infinity. By method of fuzzy hypothesis, we conclude that \[\inf_{\widehat{\psi}} \sup_{\mathbb{P}\in \mathcal{H}(\epsilon)}\mathbb{E}_{\mathbb{P}}\left[\left(\widehat{\psi} - \psi\right)^2\right] \gtrsim q^2 \asymp \frac{d}{n^2}.\] ◻

11.9 Proof of Theorem 5 ↩︎

Proof. By the same argument in Lemma 4, we only need to prove the result under the Poisson sampling model, where \((n \widehat{p}_1,\dots,n \widehat{p}_d) \text{ i.i.d} \sim \text{Poi}(\lambda)\) with \(\lambda = \frac{n}{d}\). Let \(\nu_0, \nu_1\) be two prior distributions on propensity score \(\pi\) defined on \([\epsilon, 1-\epsilon]\) satisfying

\(\mathbb{E}_{\nu_0}[\pi^j] = \mathbb{E}_{\nu_1}[\pi^j], 1\leq j \leq 2L, L = \lfloor \frac{1}{\beta} \rfloor +1\).
\(\mathbb{E}_{\nu_0}\left[\frac{1}{\pi} \right] - \mathbb{E}_{\nu_1}\left[\frac{1}{\pi} \right] = c_1(\beta, \epsilon)>0.\)

The existence of \(\nu_0, \nu_1\) follows from the duality between moment matching and best polynomial approximation. In fact, one can show \[\label{eq:const-147x} c_1(\beta, \epsilon) = 2 \mathbb{E}_{2L} \left( \frac{1}{x}; [\epsilon, 1-\epsilon] \right),\tag{27}\] where \(E_n(f;S)\) is the best polynomial (with order no greater than \(n\)) approximation error of \(f\) on the interval \(S\). Since \(L\) depends on \(\beta\), the RHS of 27 is a constant only depending on \(\beta,\epsilon\).

Under null hypothesis \(P\), let \((\pi_1, \dots, \pi_d)\) i.i.d. \(\sim \nu_0\) and \(\mu_k = \epsilon /\pi_k\). The sufficient statistics for \(\psi_1\) (conditioning on \((\pi_1, \dots, \pi_d)\)) are \[N_{k11} : = \#\{i:X_i=k,A_i=1,Y_i=1 \} \sim \text{Poi} (\epsilon \lambda) ,\] \[N_{k10} : = \#\{i:X_i=k,A_i=1,Y_i=0 \} \sim \text{Poi} \left(\lambda(\pi_k-\epsilon)\right) ,\] \[N_{k0} := \#\{i:X_i=k,A_i=0 \} \sim \text{Poi} \left(\lambda (1-\pi_k)\right)\] with \(N_{k11}, N_{k10}, N_{k0}\) conditionally independent. Similarly under the alternative hypothesis \(P'\), let \((\pi_1', \dots, \pi_d')\) i.i.d. \(\sim \nu_1\) and \(\mu_k' = \epsilon /\pi_k'\). The sufficient statistics for \(\psi_1\) (conditioning on \((\pi_1', \dots, \pi_d')\)) are \[N_{k11}' : = \#\{i:X_i=k,A_i=1,Y_i=1 \} \sim \text{Poi} (\epsilon \lambda) ,\] \[N_{k10}' : = \#\{i:X_i=k,A_i=1,Y_i=0 \} \sim \text{Poi} \left(\lambda(\pi_k'-\epsilon)\right) ,\] \[N_{k0}' := \#\{i:X_i=k,A_i=0 \} \sim \text{Poi} \left(\lambda (1-\pi_k')\right) .\] Using the same notation as in proof of Theorem 2, the total variation distance between the sufficient statistics is \[\begin{align} &\, \text{TV}(\mathbf{N}_k,\mathbf{N}_k')\\ = &\, \frac{1}{2} \sum_{i,j,\ell=0}^{\infty} \big | \mathbb{P}(N_{k11}=i, N_{k10}=j, N_{k0}=\ell)- \mathbb{P}(N_{k11}'=i, N_{k10}'=j, N_{k0}'=\ell)\big |. \end{align}\] Conditioning on \(\pi_k\) we have \[\begin{align} &\,\mathbb{P}(N_{k11}=i, N_{k10}=j, N_{k0}=\ell)\\ = & \, \mathbb{E}\left [ \exp(-\epsilon \lambda ) \frac{(\epsilon \lambda)^i}{i!} \exp(-\lambda(\pi_k-\epsilon)) \frac{[\lambda(\pi_k-\epsilon)]^j}{j!} \right.\\ &\, \left.\exp(-\lambda (1-\pi_k)) \frac{[\lambda(1-\pi_k)]^{\ell}}{\ell !} \right ] \\ = &\, \mathbb{E}\left[ \exp(-\lambda) \frac{(\epsilon \lambda)^i[\lambda(\pi_k-\epsilon)]^j[\lambda(1-\pi_k)]^{\ell}}{i!j!\ell!} \right]. \\ \end{align}\] Hence the total variation distance can be written as \[\begin{align} &\, \text{TV}(\mathbf{N}_k,\mathbf{N}_k') \\ = &\, \frac{1}{2} \sum_{i,j,\ell} \left | \mathbb{E}\left[ \exp(-\lambda) \frac{(\epsilon \lambda)^i[\lambda(\pi_k-\epsilon)]^j[\lambda(1-\pi_k)]^{\ell}}{i!j!\ell!} \right] \right. \\ &\, \left.- \mathbb{E}\left[ \exp(-\lambda) \frac{(\epsilon \lambda)^i[\lambda(\pi_k'-\epsilon)]^j[\lambda(1-\pi_k')]^{\ell}}{i!j!\ell!} \right] \right | \\ = &\, \frac{1}{2} \sum_{j,\ell} \left | \mathbb{E}\left [ \exp(-(1-\epsilon)\lambda) \frac{[\lambda(\pi_k-\epsilon)]^j[\lambda(1-\pi_k)]^{\ell}}{j!\ell!}\right]\right. \\ &\, \left. -\mathbb{E}\left [ \exp(-(1-\epsilon)\lambda) \frac{[\lambda(\pi_k'-\epsilon)]^j[\lambda(1-\pi_k')]^{\ell}}{j!\ell!}\right] \right|, \end{align}\] where the last equation follows from \(\sum_i \exp(-\epsilon \lambda) \frac{(\epsilon \lambda)^i}{i!} = 1\). Since \(\mathbb{E}[\pi_k^j] = \mathbb{E}[\pi_k'^j], 1 \leq j \leq 2L\), we have \[\begin{align} &\, \text{TV}(\mathbf{N}_k,\mathbf{N}_k') \\ = &\, \frac{1}{2} \sum_{j+\ell>2L} \left | \mathbb{E}\left [ \exp(-(1-\epsilon)\lambda) \frac{[\lambda(\pi_k-\epsilon)]^j[\lambda(1-\pi_k)]^{\ell}}{j!\ell!}\right]\right. \\ &\, \left. -\mathbb{E}\left [ \exp(-(1-\epsilon)\lambda) \frac{[\lambda(\pi_k'-\epsilon)]^j[\lambda(1-\pi_k')]^{\ell}}{j!\ell!}\right] \right| \\ \leq &\, \frac{1}{2} \left \{ \sum_{j+\ell>2L} \mathbb{E}\left [ \exp(-(1-\epsilon)\lambda) \frac{[\lambda(\pi_k-\epsilon)]^j[\lambda(1-\pi_k)]^{\ell}}{j!\ell!}\right] \right. \\ &\, \left. +\mathbb{E}\left [ \exp(-(1-\epsilon)\lambda) \frac{[\lambda(\pi_k'-\epsilon)]^j[\lambda(1-\pi_k')]^{\ell}}{j!\ell!}\right] \right \}. \\ \end{align}\] Since given \(\pi_k\), \(N_{k11} \sim \text{Poi}(\lambda(\pi_k-\epsilon)), N_{k10} \sim \text{Poi}(\lambda(1-\pi_k)), N_{k11}+N_{k10} \sim \text{Poi}(\lambda(1-\epsilon))\), we have \[\begin{align} &\,\sum_{j+\ell>2L} \mathbb{E}\left [ \exp(-(1-\epsilon)\lambda) \frac{[\lambda(\pi_k-\epsilon)]^j[\lambda(1-\pi_k)]^{\ell}}{j!\ell!}\right]\\ = &\, \mathbb{E}[\mathbb{P}(N_{k11}+N_{k10}>2L|\pi_k)] \\ \leq&\, \exp(-\lambda(1-\epsilon)) \left( \frac{e\lambda(1-\epsilon)}{2L} \right)^{2L}, \end{align}\] where the last inequality follows from Chernoff bound 22 . Hence in the regime \(n \lesssim d^{1-\beta}\) with the choice \(L=\lfloor 1/\beta \rfloor +1\), the total variation distance is bounded as \[\text{TV}(\mathbf{N}_k,\mathbf{N}_k') \leq \left( \frac{en(1-\epsilon)}{2dL} \right)^{2L} \lesssim d^{-2\beta L} = o(1/d^2).\] Since \((\mathbf{N}_1, \dots, \mathbf{N}_d)\) are independent, we have \[\text{TV}(\mathbf{N},\mathbf{N}') \leq \sum_{k=1}^d \text{TV}(\mathbf{N}_k,\mathbf{N}_k') = o(1/d) \rightarrow 0.\] The functional separation between the null and the alternative is \[\psi_1(P) - \psi_1(P') = \sum_{k=1}^d \frac{\epsilon}{d} \left( \frac{1}{\pi_k} - \frac{1}{\pi_k'} \right)\] with expectation \[\mathbb{E}[\psi_1(P) - \psi_1(P') ] = \epsilon \mathbb{E}\left[ \frac{1}{\pi_k} - \frac{1}{\pi_k'} \right]= \epsilon c_1( \beta, \epsilon) := c_2(\beta, \epsilon) = c_2.\] Define two events \[E =\left \{ \left| \frac{1}{d}\sum_{k=1}^d \frac{1}{\pi_k} - \mathbb{E}\left[\frac{1}{\pi_k}\right] \right| \leq \frac{c_2}{4} \right\},\] \[E' =\left \{ \left| \frac{1}{d}\sum_{k=1}^d \frac{1}{\pi_k'} - \mathbb{E}\left[\frac{1}{\pi_k'}\right] \right| \leq \frac{c_2}{4} \right\}.\] By Chebyshev’s inequality, we have \[\mathbb{P}(E^c) \leq \frac{16}{c_2^2d} \operatorname{Var}\left(\frac{1}{\pi_k} \right) \leq \frac{16}{\epsilon^2c_2^2d}.\] \[\mathbb{P}(E'^c) \leq \frac{16}{c_2^2d} \operatorname{Var}\left(\frac{1}{\pi_k'} \right) \leq \frac{16}{\epsilon^2c_2^2d}.\] We put the following prior distributions induced by \(\{ \pi_k, 1 \leq k \leq d\}\) and \(\{ \pi_k', 1 \leq k \leq d\}\) on \(P\) and \(P'\), respectively: \[\pi \stackrel{d}{=} P \mid E, \pi' \stackrel{d}{=} P' \mid E'.\] Note that under \(\pi, \pi'\), \[|\psi_1(P)-\psi_1(P')|\geq c_2/2.\] By triangle inequality, the total variation distance of the sufficient counting statistics \(\mathbf{N}\) and \(\mathbf{N}'\) under two priors is bounded by \[\begin{align} \text{TV}\left( {\mathbf{N}\mid E}, {\mathbf{N}' \mid E'} \right) \leq &\, \text{TV}\left( {\mathbf{N}\mid E}, {\mathbf{N}} \right) + \text{TV}({\mathbf{N}},{\mathbf{N}'} ) + \text{TV}\left( {\mathbf{N}' \mid E'}, {\mathbf{N}'} \right)\\ \leq &\, \mathbb{P}(E^c) + \mathbb{P}(E'^c) + \text{TV}({\mathbf{N}},{\mathbf{N}'} )\\ \leq &\, \frac{32}{ \epsilon^2 c_2^2 d} + \text{TV}({\mathbf{N}},{\mathbf{N}'} ) \rightarrow 0 . \end{align}\] By method of fuzzy hypotheses, we conclude \[\begin{align} &\, \inf_{\widehat{\psi}_1} \sup_{\mathbb{P}\in \mathcal{D}^U (\epsilon)} \mathbb{E}_{\mathbb{P}}\left[\left(\widehat{\psi}_1 - \psi_1\right)^2\right] \\ \geq &\, \frac{c_2^2}{32} \left( 1 - \text{TV}\left( {\mathbf{N}\mid E}, {\mathbf{N}' \mid E'} \right) \right).\\ \geq &\, c(\beta, \epsilon), \end{align}\] for some constant \(c(\beta, \epsilon)\) when \(d\) is large enough in the regime \(n \lesssim d^{1-\beta}\). ◻

11.10 Proof of Theorem 6 ↩︎

Proof. Note that \(\eta = \mathbb{E}[(A-\pi_X)(Y-\mu_X)]\) The conditional bias is (Let \(\mathbf{Z}=(X,A,Y), \mathbf{Z}_1=(X_1,A_1,Y_1), \mathbf{Z}_2=(X_2,A_2,Y_2)\) be samples independent of \(D\)) \[\begin{align} &\, \mathbb{E}[\widehat{\eta}\mid D]-\eta \\ =&\, \mathbb{E}[(A-\widehat{\pi}_X)(Y-\widehat{\mu}_X)\mid D] - \mathbb{E}[(A-\pi_X)(Y-\mu_X)\mid D] \\ &\,- \mathbb{E}\left[\frac{(A_1-\widehat{\pi}_{X_1})I(X_1=X_2)(Y_2 - \widehat{\mu}(X_2))}{\widehat{p}_{X_1}} \mid D\right] \end{align}\] Note that \[\begin{align} &\, \mathbb{E}[(A-\widehat{\pi}_X)(Y-\widehat{\mu}_X)\mid D] - \mathbb{E}[(A-\pi_X)(Y-\mu_X)\mid D] \\ = &\, \mathbb{E}[(A-\widehat{\pi}_X)(Y-\widehat{\mu}_X)\mid D] - \mathbb{E}[(A-\widehat{\pi}_X)(Y-{\mu}_X)\mid D] \\ &\, +\mathbb{E}[(A-\widehat{\pi}_X)(Y-{\mu}_X)\mid D] - \mathbb{E}[(A-{\pi}_X)(Y-{\mu}_X)\mid D] \\ = &\, \mathbb{E}[(A-\widehat{\pi}_X)(\mu_X-\widehat{\mu}_X)\mid D] + \mathbb{E}[(\pi_X-\widehat{\pi}_X)(Y-{\mu}_X)\mid D] \\ = &\, \mathbb{E}[(\pi_X-\widehat{\pi}_X)(\mu_X-\widehat{\mu}_X)\mid D], \end{align}\] the last equation follows from conditioning on \(X\). By conditioning on \(X_1\) we have \[\mathbb{E}\left[ I(X_2=X_1)(Y_2 - \widehat{\mu}(X_2))\mid D,X_1\right] = p_{X_1}(\mu_{X_1}-\widehat{\mu}_{X_1}).\] Hence condition on \(\mathbf{Z}_1\) we have \[\begin{align} &\, \mathbb{E}\left[\frac{(A_1-\widehat{\pi}_{X_1})I(X_1=X_2)(Y_2 - \widehat{\mu}(X_2))}{\widehat{p}_{X_1}} \mid D\right]\\ = &\, \mathbb{E}\left[ \frac{p_{X_1}}{\widehat{p}_{X_1}}(A_1-\widehat{\pi}_{X_1})(\mu_{X_1}-\widehat{\mu}_{X_1}) \mid D\right] \\ = &\, \mathbb{E}\left[ \frac{p_{X_1}}{\widehat{p}_{X_1}}(\pi_{X_1}-\widehat{\pi}_{X_1})(\mu_{X_1}-\widehat{\mu}_{X_1}) \mid D\right]. \end{align}\] Hence the conditional bias is \[\label{eq:high-order-bias} \mathbb{E}[\widehat{\eta}\mid D]-\eta = \mathbb{E}\left[ (\widehat{\mu}_{X}-\mu_{X})(\widehat{\pi}_{X} - \pi_X)\left (1-\frac{p_{X}}{\widehat{p}_X} \right) \mid D \right].\tag{28}\] By Cauchy-Schwarz inequality one can bound the conditional bias as \[|\mathbb{E}[\widehat{\eta}\mid D]-\eta| \leq \|\widehat{\mu}-\mu\|_2 \|\widehat{\pi}-\pi\|_2 \max_{k} \left| 1-\frac{p_k}{\widehat{p}_k} \right|\] We next bound the variance. Let \[h_2(\mathbf{Z}_1,\mathbf{Z}_2) = (A_1 - \widehat{\pi}_{X_1})(Y_1-\widehat{\mu}_{X_1}) - \frac{(A_1-\widehat{\pi}_{X_1})I(X_1=X_2)(Y_2-\widehat{\mu}_{X_2})}{\widehat{p}_{X_1}}\] and note that the estimator can be expressed as \[\widehat{\eta} = \mathbb{U}_n[h_2(\mathbf{Z}_1,\mathbf{Z}_2)].\] We will use the variance of U-statistics (Lemma 6 in [14]) to bound the conditional variance of \(\widehat{\eta}\). Let \[h_1(\mathbf{Z}_1) = \mathbb{E}[h_2(\mathbf{Z}_1,\mathbf{Z}_2)\mid\mathbf{Z}_1] = (A_1 - \widehat{\pi}_{X_1})(Y_1-\widehat{\mu}_{X_1}) - \frac{(A_1-\widehat{\pi}_{X_1})p_{X_1}(\mu_{X_1}-\widehat{\mu}_{X_1})}{\widehat{p}_{X_1}}.\] \[\begin{align} \sigma_1^2 =&\, \mathbb{E}[h_1^2(\mathbf{Z}_1)\mid D] \\ \leq &\, 2 \left( \mathbb{E}[(A_1 - \widehat{\pi}_{X_1})^2(Y_1-\widehat{\mu}_{X_1})^2\mid D] + \mathbb{E}\left[ \frac{(A_1-\widehat{\pi}_{X_1})^2p_{X_1}^2(\mu_{X_1}-\widehat{\mu}_{X_1})^2}{\widehat{p}_{X_1}^2} \mid D \right] \right) \\ \lesssim &\, 1 + \mathbb{E}\left[ \frac{p_{X}^2}{\widehat{p}_X^2} \mid D \right] \\ \lesssim &\, 1 + \left(1+\max_{k} \left| 1-\frac{p_k}{\widehat{p}_k} \right|\right)^2 \\ \lesssim &\, \left(1+\max_{k} \left| 1-\frac{p_k}{\widehat{p}_k} \right|\right)^2. \end{align}\] \[\begin{align} \sigma_2^2 =&\, \mathbb{E}[h_2^2(\mathbf{Z}_1,\mathbf{Z}_2)\mid D] \\ \leq &\, 2 \left( \mathbb{E}[(A_1 - \widehat{\pi}_{X_1})^2(Y_1-\widehat{\mu}_{X_1})^2\mid D] + \mathbb{E}\left[ \frac{(A_1-\widehat{\pi}_{X_1})^2I(X_1=X_2)(Y_2-\widehat{\mu}_{X_2})^2}{\widehat{p}_{X_1}^2} \mid D \right] \right) \\ \lesssim &\, 1 + \mathbb{E}\left [ \frac{I(X_1=X_2)}{\widehat{p}_{X_1}^2} \mid D \right] \\ \lesssim &\, 1 + \sum_{k=1}^d \frac{p_k^2}{\widehat{p}_k^2} \\ \leq &\, 1 + d\left(1+\max_{k} \left| 1-\frac{p_k}{\widehat{p}_k} \right|\right)^2 \\ \lesssim&\, d\left(1+\max_{k} \left| 1-\frac{p_k}{\widehat{p}_k} \right|\right)^2. \end{align}\] By Lemma 6 in [14] the conditional variance is upper bounded as \[\operatorname{Var}(\widehat{\eta}\mid D) \lesssim \frac{1}{n}\left(1+\max_{k} \left| 1-\frac{p_k}{\widehat{p}_k} \right|\right)^2 + \frac{d}{n^2}\left(1+\max_{k} \left| 1-\frac{p_k}{\widehat{p}_k} \right|\right)^2.\] ◻

11.11 Proof of Theorem 7 ↩︎

Proof. We will use the following equation for MSE \[\label{eq:mse-decomposition} \begin{align} &\, \mathbb{E}[(\widehat{\eta}-\eta)^2] \\ =&\, \mathbb{E}\left\{ \mathbb{E}[(\widehat{\eta}-\eta)^2\mid D] \right\} \\ = &\, \mathbb{E}[\operatorname{Var}(\widehat{\eta}\mid D)] + \mathbb{E}\{(\mathbb{E}[\widehat{\eta}\mid D] - \eta)^2\}. \end{align}\tag{29}\] For the conditional bias derived in 28 we have (using the property \((\mathbb{E}[X])^2 \leq \mathbb{E}[X^2]\)) \[\begin{align} &\,\mathbb{E}\{(\mathbb{E}[\widehat{\eta}\mid D] - \eta)^2\}\\ \leq &\, \mathbb{E}\left \{ \mathbb{E}\left [ (\widehat{\mu}_{X}-\mu_{X})^2(\widehat{\pi}_{X} - \pi_X)^2\left (1-\frac{p_{X}}{\widehat{p}_X} \right)^2 \mid D\right] \right\} \\ = &\, \mathbb{E}\left[ \sum_{k=1}^d p_k (\widehat{\mu}_{k}-\mu_{k})^2(\widehat{\pi}_{k} - \pi_k)^2\left (1-\frac{p_{k}}{\widehat{p}_k} \right)^2\right] \\ = &\, \sum_{k=1}^d \mathbb{E}\left[ p_k (\widehat{\mu}_{k}-\mu_{k})^2(\widehat{\pi}_{k} - \pi_k)^2\left (1-\frac{p_{k}}{\widehat{p}_k} \right)^2\right] \\ \leq&\, \xi_n^2 \sum_{k=1}^d p_k\mathbb{E}[(\widehat{\pi}_k - \pi_k)^2], \end{align}\] where in the last inequality we use the bound \[\max_k \left |1-\frac{p_k}{\widehat{p}_k} \right| \leq \xi_n.\] and \((\widehat{\mu}_k - \mu_k)^2 \leq 1\). A naive bound \[\sum_{k=1}^d p_k\mathbb{E}[(\widehat{\pi}_k - \pi_k)^2] \leq 1\] holds for both empirical average estimator \(\widehat{\pi}_k, \widehat{\mu}_k\) or simply letting \(\widehat{\pi}_k = \widehat{\mu}_k =0\). For empirical average estimator \(\widehat{\pi}_k\) from a sample of size \(n\) we can derive an alternative bound. By the property of conditional variance, we have \[\operatorname{Var}(\widehat{\pi}_k) = \mathbb{E}[\operatorname{Var}(\widehat{\pi}_k \mid \mathbf{X}^n)] + \operatorname{Var}(\mathbb{E}[\widehat{\pi}_k \mid \mathbf{X}^n]).\] Recall \(\mathbb{E}[\widehat{\pi}_k \mid \mathbf{X}^n] = \pi_k I(\widehat{p}_k>0), \operatorname{Var}(\widehat{\pi}_k \mid \mathbf{X}^n) = \frac{\pi_k(1-\pi_k)}{n\widehat{p}_k} I(\widehat{p}_k>0)\) we obtain \[\operatorname{Var}(\widehat{\pi}_k) \leq \frac{1}{4} \mathbb{E}\left[\frac{I(\widehat{p}_k>0)}{n\widehat{p}_k}\right] + \pi_k^2 (1-p_k)^n \leq \frac{1}{2(n+1)p_k} + (1-p_k)^n,\] \[(\mathbb{E}[\widehat{\pi}_k ] - \pi_k)^2 = \pi_k^2 (1-p_k)^{2n} \leq (1-p_k)^{2n}.\] Thus we have \[\mathbb{E}[(\widehat{\pi}_k - \pi_k)^2] \leq \frac{1}{2(n+1)p_k} + (1-p_k)^n + (1-p_k)^{2n}\] \[\sum_{k=1}^d p_k\mathbb{E}[(\widehat{\pi}_k - \pi_k)^2] \leq \frac{d}{2(n+1)} + \sum_{k=1}^d p_k(1-p_k)^n +\sum_{k=1}^d p_k(1-p_k)^{2n} .\] Let \(f_6(x) = x(1-x)^n, x\in [0,1]\) and \(f_6'(x) = (1-nx)(1-x)^{n-1}\). Hence \(f_6(x) \leq f_6(1/n) < 1/n\), which implies \[\sum_{k=1}^d p_k(1-p_k)^n \leq \frac{d}{n}, \sum_{k=1}^d p_k(1-p_k)^{2n} \leq \frac{d}{2n}.\] We conclude that \[\sum_{k=1}^d p_k\mathbb{E}[(\widehat{\pi}_k - \pi_k)^2] \leq \frac{2d}{n},\] Combining this bound the naive bound \[\sum_{k=1}^d p_k\mathbb{E}[(\widehat{\pi}_k - \pi_k)^2] \leq 1\] we have \[\mathbb{E}\{(\mathbb{E}[\widehat{\eta}\mid D] - \eta)^2\} \lesssim \xi_n^2 \frac{n \wedge d}{n}.\] The bound on conditional variance in proof of Theorem 6 can be reduced to \[\operatorname{Var}(\widehat{\eta}\mid D) \lesssim \frac{1}{n} + \frac{d}{n^2}.\] The proof is completed by combining the bounds on conditional bias with variance as in 29 . ◻

11.12 Proof of Theorem 8 ↩︎

Proof. We will show the plug-in estimator of \(\psi_1 = \sum_{k}p_k \mu_{1k}\) satisfies \[\widehat{\psi}_1 - \psi_1 = (\mathbb{P}_n - \mathbb{P}) [\varphi_1(\mathbf{Z})] + o_{\mathbb{P}}(1/\sqrt{n}),\] where \[\varphi_1(\mathbf{Z}) = \frac{A(Y-\mu_{1X})}{\pi_X}+\mu_{1X}.\] Similar argument can be applied to \(\psi_0 = \sum_{k}p_k \mu_{0k}\). By Proposition 1 we will use the doubly robust form of \(\widehat{\psi} = \mathbb{P}_n[\widehat{\varphi}_1(\mathbf{Z})]\). Using the following decomposition \[\widehat{\psi}_1 - \psi_1 = (\mathbb{P}_n - \mathbb{P}) [\varphi_1(\mathbf{Z})] + (\mathbb{P}_n - \mathbb{P}) [\widehat{\varphi}_1(\mathbf{Z})-\varphi_1(\mathbf{Z})] + \mathbb{P}[\widehat{\varphi}_1(\mathbf{Z})-\varphi_1(\mathbf{Z})].\] Step1: Bound \(\mathbb{P}[\widehat{\varphi}_1(\mathbf{Z})-\varphi_1(\mathbf{Z})]\). We first show \(\mathbb{P}[\widehat{\varphi}_1(\mathbf{Z})-\varphi_1(\mathbf{Z})] = \int \widehat{\varphi}_1(\mathbf{z})-\varphi_1(\mathbf{z}) d \mathbb{P}(\mathbf{z}) = o_{\mathbb{P}}(1/\sqrt{n})\). By direct calculations, one can show \[|\mathbb{P}[\widehat{\varphi}_1(\mathbf{Z})-\varphi_1(\mathbf{Z})] |= \bigg |\mathbb{P}\left[ (\widehat{\mu}_{1X}-\mu_{1X})\left(1-\frac{\pi_X}{\widehat{\pi}_X} \right) \right] \bigg |\leq \|\widehat{\mu}_{1X}-\mu_{1X}\|_2 \left \| 1-\frac{\pi_X}{\widehat{\pi}_X}\right \|_2.\] Note that since \(n\widehat{p}_k \widehat{\pi}_k \widehat{\mu}_{1k}| \mathbf{X}^n, \mathbf{A}^n \sim \text{Bin}(n\widehat{p}_k \widehat{\pi}_k, \mu_{1k})\), we have \[\mathbb{E}[\widehat{\mu}_{1k}\mid\mathbf{X}^n,\mathbf{A}^n] = \mu_{1k} I(\widehat{p}_k \widehat{\pi}_k>0), \, \operatorname{Var}(\widehat{\mu}_{1k}\mid\mathbf{X}^n, \mathbf{A}^n) = \frac{\mu_{1k}(1-\mu_{1k})I(\widehat{p}_k \widehat{\pi}_k>0)}{n\widehat{p}_k \widehat{\pi}_k}.\] The bias of \(\widehat{\mu}_{1k}\) is \[\mathbb{E}[\widehat{\mu}_{1k}-\mu_{1k}] = -\mu_{1k}\mathbb{P}(\widehat{p}_k \widehat{\pi}_k=0).\] By conditioning on \(\mathbf{X}^n\), we have \[\mathbb{P}(\widehat{p}_k \widehat{\pi}_k=0) = \mathbb{E}[\mathbb{P}(\widehat{p}_k \widehat{\pi}_k=0\mid \mathbf{X}^n)] = \mathbb{E}[(1-\pi_k)^{n\widehat{p}_k}] = (1-p_k \pi_k)^n \leq (1-\epsilon p_k)^n,\] where we use the fact \(\mathbb{E}[c^V] = (1-p+pc)^n\) for \(V \sim B(n,p)\). Similarly, we have \[\operatorname{Var}(\mathbb{E}[\widehat{\mu}_{1k}\mid\mathbf{X}^n,\mathbf{A}^n]) = \mu_{1k}^2 \mathbb{P}(\widehat{p}_k \widehat{\pi}_k=0)\mathbb{P}(\widehat{p}_k \widehat{\pi}_k>0) \leq (1-\epsilon p_k)^n.\] For the expected conditional variance, \[\mathbb{E}\left[ \operatorname{Var}(\widehat{\mu}_{1k}\mid\mathbf{X}^n, \mathbf{A}^n) \right] \leq \frac{1}{4}\mathbb{E}\left[\frac{I(\widehat{p}_k \widehat{\pi}_k>0)}{n\widehat{p}_k \widehat{\pi}_k} \right]\] Note that \(n\widehat{p}_k \widehat{\pi}_k \sim B(n,p_k \pi_k)\), by lemma 2 we have \[\mathbb{E}\left[ \operatorname{Var}(\widehat{\mu}_{1k}\mid\mathbf{X}^n, \mathbf{A}^n) \right] \leq \frac{1}{2(n+1)p_k\pi_k}\leq \frac{1}{2\epsilon(n+1)p_k}.\] Thus the variance of \(\widehat{\mu}_{1k}\) is bounded by \[\operatorname{Var}(\widehat{\mu}_{1k}) \leq (1-\epsilon p_k)^n + \frac{1}{2\epsilon(n+1)p_k}.\] We conclude \[\mathbb{E}[(\widehat{\mu}_{1k}-\mu_{1k})^2] \leq (1-\epsilon p_k)^n + \frac{1}{2\epsilon(n+1)p_k}+ (1-\epsilon p_k)^{2n} = O(1/n)\] since when \(d\) is fixed, \(p_k\)’s are considered as fixed. Hence \[\|\widehat{\mu}_{1X}-\mu_{1X}\|_2^2 = \sum_{k}p_k (\widehat{\mu}_{1k}-\mu_{1k})^2 = O_{\mathbb{P}}(1/n),\] \[\|\widehat{\mu}_{1X}-\mu_{1X}\|_2 = O_{\mathbb{P}}(1/\sqrt{n}).\] Similarly \(\mathbb{E}[\widehat{\pi}_k \mid \mathbf{X}^n] = \pi_k I(\widehat{p}_k>0), \operatorname{Var}(\widehat{\pi}_k \mid \mathbf{X}^n) = \frac{\pi_k(1-\pi_k)}{n\widehat{p}_k} I(\widehat{p}_k>0)\), we obtain \[\operatorname{Var}(\widehat{\pi}_k) \leq \frac{1}{4} \mathbb{E}\left[\frac{I(\widehat{p}_k>0)}{n\widehat{p}_k}\right] + \pi_k^2 (1-p_k)^n \leq \frac{1}{2(n+1)p_k} + (1-p_k)^n,\] \[(\mathbb{E}[\widehat{\pi}_k ] - \pi_k)^2 = \pi_k^2 (1-p_k)^{2n} \leq (1-p_k)^{2n}.\] Thus we have \[\mathbb{E}[(\widehat{\pi}_k - \pi_k)^2] \leq \frac{1}{2(n+1)p_k} + (1-p_k)^n + (1-p_k)^{2n} = O(1/n)\] \[(\widehat{\pi}_k-\pi_k)^2 = O_{\mathbb{P}}(1/n).\] Since \(\pi_k \geq \epsilon\), this shows \(1/\widehat{\pi}_k = O_{\mathbb{P}}(1)\) and \[\frac{(\widehat{\pi}_k-\pi_k)^2}{\widehat{\pi}_k^2} =O_{\mathbb{P}}(1/n).\] We conclude \[\left \| 1-\frac{\pi_X}{\widehat{\pi}_X}\right \|_2^2 = \sum_{k} p_k \frac{(\widehat{\pi}_k-\pi_k)^2}{\widehat{\pi}_k^2}=O_{\mathbb{P}}(1/n),\] \[\left \| 1-\frac{\pi_X}{\widehat{\pi}_X}\right \|_2=O_{\mathbb{P}}(1/\sqrt{n}).\] This shows \[\mathbb{P}[\widehat{\varphi}_1(\mathbf{Z})-\varphi_1(\mathbf{Z})] = O_{\mathbb{P}}(1/n)=o_{\mathbb{P}}(1/\sqrt{n}).\] Step2: Bound \((\mathbb{P}_n - \mathbb{P}) [\widehat{\varphi}_1(\mathbf{Z})-\varphi_1(\mathbf{Z})]\).We then show \[(\mathbb{P}_n - \mathbb{P}) [\widehat{\varphi}_1(\mathbf{Z})-\varphi_1(\mathbf{Z})] = o_{\mathbb{P}}(1/\sqrt{n}).\] Since \(X\) is discrete we can write the nuisance functions \((\widehat{\pi}, \widehat{\mu}_1)\) as saturated linear models, i.e. \[{\pi}_x=\pi(x ; \boldsymbol{\alpha})=\boldsymbol{\alpha}^{\top} \boldsymbol{w},\] where \(\boldsymbol{w}^{\top}=\{I(x=1), \ldots, I(x=d)\} \in\{0,1\}^{d}\) and \[{\mu}_{1x}=\mu_1(x ; \boldsymbol{\beta})=\boldsymbol{\beta}^{\top} \boldsymbol{w}.\] Here the components of \(\boldsymbol{\alpha}\) and \(\boldsymbol{\beta}\) are simply propensity scores and regression functions within different categories and \(\|\boldsymbol{w}\|_2=1\). Define the function class \(\mathcal{F} = \{\varphi_{1}(\mathbf{z};\boldsymbol{\gamma}), \boldsymbol{\gamma}= (\boldsymbol{\alpha}, \boldsymbol{\beta}), \alpha_k \in [\epsilon/2, 1], \beta_k \in [0, 1]\}\). Since propensity scores are lower bounded by \(\epsilon /2\) for \(\boldsymbol{\alpha}\in \mathcal{F}\), one can show there is a constant \(K\) that depends on \(\epsilon\) such that \[|\varphi_{1}(\mathbf{z};\boldsymbol{\gamma}_1) - \varphi_{1}(\mathbf{z};\boldsymbol{\gamma}_2)| \leq K\|\boldsymbol{\gamma}_1 - \boldsymbol{\gamma}_2\|_2.\] Since \(d\) is fixed as constant, by example 19.7 in [82], \(\mathcal{F}\) is Donsker. Let \[\widetilde{\pi}_k = \widehat{\pi}_k I(\widehat{\pi}_k \geq \epsilon/2) + \frac{\epsilon}{2} I(\widehat{\pi}_k <\epsilon/2), \, \widetilde{\mu}_{1k} = \widehat{\mu}_{1k},\] so that we truncate the estimated propensity score \(\widehat{\pi}_k\) to obtain \(\widetilde{\pi}_k\). Let \[\widetilde{\varphi}_1(\mathbf{Z}) = \varphi_1(\mathbf{Z};\widetilde{\boldsymbol{\gamma}}) = \frac{A(Y-\widetilde{\mu}_{1X})}{\widetilde{\pi}_X}+\widetilde{\mu}_{1X}.\] Clearly \(\widetilde{\varphi}_{1}(\cdot) = \varphi_1(\cdot, \widetilde{\boldsymbol{\gamma}}) \in \mathcal{F}\). We have \[\|\widetilde{\boldsymbol{\gamma}}- \boldsymbol{\gamma}\|_2^2 =\|\widetilde{\boldsymbol{\alpha}} -\boldsymbol{\alpha}\|_2^2+\|\widetilde{\boldsymbol{\beta}}-\boldsymbol{\beta}\|_2^2 = \sum_{k=1}^d \left[(\widetilde{\pi}_k-\pi_k)^2 + (\widetilde{\mu}_{1k}-\mu_{1k})^2 \right] .\] Since we assume \(\pi_k \geq \epsilon\), truncating \(\widehat{\pi}_k\) yields smaller error and hence \(|\widetilde{\pi}_k - \pi_k|\leq |\widehat{\pi}_k - \pi_k|\). By the consistency of \(\widehat{\pi}_k\) and \(\widehat{\mu}_{1k}\) shown above, we have \[\|\widetilde{\boldsymbol{\gamma}}- \boldsymbol{\gamma}\|_2^2 \leq \sum_{k=1}^d \left[(\widehat{\pi}_k-\pi_k)^2 + (\widehat{\mu}_{1k}-\mu_{1k})^2 \right] = o_{\mathbb{P}}(1)\] again because \(d\) is fixed. Thus \[\|\widetilde{\varphi}_1-\varphi\|_2 \leq K \|\widetilde{\boldsymbol{\gamma}}- \boldsymbol{\gamma}\| = o_{\mathbb{P}}(1)\] and Lemma 19.24 in [82] shows \[\label{eq:varphitilde-neg} (\mathbb{P}_n - \mathbb{P}) [\widetilde{\varphi}_1(\mathbf{Z})-\varphi_1(\mathbf{Z})] = o_{\mathbb{P}}(1/\sqrt{n}).\tag{30}\] Now consider \[(\mathbb{P}_n - \mathbb{P}) [\widehat{\varphi}_1(\mathbf{Z})-\widetilde{\varphi}_1(\mathbf{Z})].\] By strong law of large numbers, we have \[\widehat{\pi}_k = \frac{\sum_{i=1}^n I(X_i=k,A_i=1)}{\sum_{i=1}^nI(X_i=k)} \rightarrow \pi_k\] for any \(k \in [d]\) almost surely. Thus for almost every \(\omega \in \Omega\) (sample space), one can find a \(N(\omega,\epsilon,d) \in \mathbf{N}^+\) such that for all \(n \geq N(\omega, \epsilon,d)\), we have \(|\widehat{\pi}_k-\pi_k| \leq \epsilon/2\) for any \(k \in [d]\). This together with \(\pi_k \geq \epsilon\) shows \(\widehat{\pi}_k \geq \epsilon/2\) and hence \(\widehat{\pi}_k = \widetilde{\pi}_k\) holds for all \(k\) when \(n \geq N(\omega,\epsilon,d)\). This implies \[\widetilde{\varphi}_1(\mathbf{z}) \equiv \widehat{\varphi}_1(\mathbf{z}),\] \[(\mathbb{P}_n - \mathbb{P}) [\widehat{\varphi}_1(\mathbf{Z})-\widetilde{\varphi}_1(\mathbf{Z})]=0\] when \(n \geq N(\omega,\epsilon,d)\). This clearly implies \[\sqrt{n}(\mathbb{P}_n - \mathbb{P}) [\widehat{\varphi}_1(\mathbf{Z})-\widetilde{\varphi}_1(\mathbf{Z})] \rightarrow 0\] almost surely since the left-hand side is exactly \(0\) when \(n\) is large enough. We conclude \[\label{eq:diff-varphi} (\mathbb{P}_n - \mathbb{P}) [\widehat{\varphi}_1(\mathbf{Z})-\widetilde{\varphi}_1(\mathbf{Z})] = o_{\mathbb{P}}(1/\sqrt{n})\tag{31}\] Combining 30 and 31 , we have \[(\mathbb{P}_n - \mathbb{P}) [\widehat{\varphi}_1(\mathbf{Z})-\varphi_1(\mathbf{Z})] = o_{\mathbb{P}}(1/\sqrt{n}).\] ◻

12 Proof of Auxiliary Results↩︎

12.1 Proof of Lemma 3 ↩︎

Proof. First consider the boundary. If \(y=0\) then the function \(f_2\) is reduced to \[g_1(x) = (n-1)(1-x)^{n-2} - n(1-x)^{n-1}, 0 \leq x \leq 1,\] \[g_1'(x) = (n-1)(1-x)^{n-3}(2-nx).\] So \(|g_1(x)|\leq \max (|g_1(0)|, |g_1(1)|, |g_1(2/n)|)\) and \(g_1(0) = -1, g_1(1) = 0, g_1(2/n) = \left(\frac{n-2}{n}\right)^{n-2} \leq 1\). This shows \(|g_1(x)| \leq 1\). On the boundary \(x=0\) similar arguments hold. Now consider the boundary \(x+y=1\), the function \(f_2\) is reduced to \[|f_2(x,y)| = nx^{n-1}(1-x)^{n-1} \leq \frac{n}{4^{n-1}} \leq 1\] Now we consider the interior of the triangle. By taking the partial derivative (or by noting that \(x,y\) have symmetric roles in the function \(f_2\)) we see the maximizer of \(|f_2(x,y)|\) must lie on the line \(x=y\). So we define \[g_2(x) = (n-1)(1-2x)^{n-2} - n(1-x)^{2n-2}, 0 \leq x \leq 1/2\] and only need to show \(|g_2(x)| \leq 1\). Let \(X_0\) be the set of stationary points of \(g_2\) on \([0,1/2]\). Any stationary point \(x_0 \in X_0\) must satisfy \(g_2'(x_0) = 0\), i.e. \[n(1-x_0)^{2n-3} = (n-2)(1-2x_0)^{n-3}.\] So for any \(x_0 \in X_0\), \[g_2(x_0) = (n-1)(1-2x_0)^{n-2} - (n-2)(1-2x_0)^{n-3}(1-x_0) = (1-2x_0)^{n-3}(1-nx_0).\] We then define a new function to show \(|g_2(x_0)| \leq 1\) for any \(x_0 \in X_0\). Let \[g_3(x) = (1-2x)^{n-3}(1-nx), 0 \leq x \leq 1/2\] \[g_3'(x) = (n-2)(2nx-3)(1-2x)^{n-4}.\] So we have \(|g_3(x)| \leq \max(|g_3(0)|,|g_3(1/2)|, |g_3(3/2n)|)\) and \(g_3(0) = 1, g_3(1/2) = 0, |g_3(3/2n)| = |\frac{1}{2}(\frac{n-3}{n})^{n-3}| \leq 1/2\). This shows \(|g_3(x)| \leq 1\) on \([0,1/2]\), which implies \(|g_2(x_0)| \leq 1\) for any \(x \in X_0\). Note that \(g_2(0)=-1, g_2(1/2) = -n/4^{n-1}\) We conclude \[|g_2(x)| \leq \max\left( \max_{x_0 \in X_0 }|g_2(x)|, 1, n/4^{n-1} \right) \leq 1\] ◻

12.2 Proof of Lemma 4 ↩︎

Proof. For any \(\gamma > 0\), let \(\widehat{\psi}_1(n)\) be a near-minimax optimal estimator of \(\psi_1({\boldsymbol{p}})\) for fixed sample \(n\), i.e., \[\sup_{\mathbb{P}\in \mathcal{D}(\epsilon)} \mathbb{E}_{\mathbb{P}}\left[\left(\widehat{\psi}_1(n) - \psi_1({\boldsymbol{p}})\right)^2\right] \leq \gamma + R^*(d,n).\] Note we emphasize the dependency of \(\psi_1\) on \({\boldsymbol{p}}=(p_1,\dots)\). Under a Poisson-sampling model \(\mathbb{P}\in \mathcal{D}(\epsilon, \delta)\), let \(n' = \sum_{k}n\widehat{p}_k \sim \text{Poi}(n\sum_{k} p_k)\) and construct an estimator \(\widetilde{\psi}_1 = \widehat{\psi}_1(n')\). Note that conditioned on \(n'=m, (n\widehat{p}_1, \dots,) \sim \text{Multinomial}\left(m, \frac{{\boldsymbol{p}}}{\sum_{k}p_k}\right)\). By definition, the model with probability vector \(\frac{{\boldsymbol{p}}}{\sum_{k}p_k}\) and the same propensity score, regression functions as \(\mathbb{P}\) is in \(\mathcal{D}(\epsilon)\). Under the Poisson-sampling model we have \[\begin{align} &\, \mathbb{E}\left[ \left( \widetilde{\psi}_1 -\psi_1\left( \frac{{\boldsymbol{p}}}{\sum_{k}p_k}\right) \right)^2 \right] \\ = &\, \sum_{m=0}^{\infty} \mathbb{E}\left[ \left( \widetilde{\psi}_1 -\psi_1\left( \frac{{\boldsymbol{p}}}{\sum_{k}p_k}\right) \right)^2 \mid n'=m\right] \mathbb{P}(n'=m)\\ \leq & \, \sum_{m=0}^{\infty} R^*(d,m)\mathbb{P}(n'=m) + \gamma. \end{align}\] Since \(n \mapsto R^*(d, n)\) is non-increasing for fixed \(d\) and \(R^*(d, n) \leq 1\) , we have \[\begin{align} &\,\mathbb{E}\left[ \left( \widetilde{\psi}_1 -\psi_1\left( \frac{{\boldsymbol{p}}}{\sum_{k}p_k}\right) \right)^2 \right] \\ \leq &\, \sum_{m \geq n / 2} R^*(d, m) \mathbb{P}\left[n^{\prime}=m\right]+ \mathbb{P}\left[n^{\prime} \leq \frac{n}{2}\right] + \gamma \\ \leq &\, R^*(d, n/2) + \exp(-n/50)+\gamma. \end{align}\] The last inequality follows from Chernoff bound and \(|\sum_{k}p_k-1|\leq \delta < 1/3\). The difference between \(\psi_1\left( \frac{{\boldsymbol{p}}}{\sum_{k}p_k}\right)\) and \(\psi_1\left( {{\boldsymbol{p}}}\right)\) is \[\begin{align} &\, \bigg |\psi_1\left( \frac{{\boldsymbol{p}}}{\sum_{k}p_k}\right) - \psi_1\left( {{\boldsymbol{p}}}\right) \bigg |\\ = &\, \big|\sum_{k}p_k-1\big| \frac{\sum_k p_k \mu_{1k}}{\sum_k p_k} \\ \leq &\, \delta. \end{align}\] The last inequality follows since \(\frac{\sum_k p_k \mu_{1k}}{\sum_k p_k}\) is a weighted average of quantities bounded by 1. Thus we have \[\begin{align} &\, \frac{1}{2} \mathbb{E}\left[ \left( \widetilde{\psi}_1 -\psi_1\left( {{\boldsymbol{p}}}\right) \right)^2 \right]\\ \leq &\, \mathbb{E}\left[ \left( \widetilde{\psi}_1 -\psi_1\left( \frac{{\boldsymbol{p}}}{\sum_{k}p_k}\right) \right)^2 \right] +\left(\psi_1\left( \frac{{\boldsymbol{p}}}{\sum_{k}p_k}\right) - \psi_1\left( {{\boldsymbol{p}}}\right) \right)^2 \\ \leq &\, R^*(d, n/2) + \exp(-n/50)+\gamma + \delta^2. \end{align}\] Take supremum over \(\mathbb{P}\in D(\epsilon, \delta)\) and since \(\gamma\) is arbitrary, we have \[\frac{1}{2}\widetilde{R}^*(d, n, \delta) \leq R^*(d, n/2) + \exp(-n/50) + \delta^2.\] ◻

12.3 Proof of Lemma 5 ↩︎

Proof. We claim it suffices to find probability measures \(\omega_0, \omega_1\) on \([c/K^2,1]\) such that

\[\label{eq:matched-moments} \mathbb{E}_{\omega_0}[X^l] = \mathbb{E}_{\omega_1}[X^l] \text{ for all } l=-1,0,1,\dots, 3K.\tag{32}\]
\[\left | \mathbb{E}_{\omega_0}\left[\frac{X}{X+c/K^2} \right]- \mathbb{E}_{\omega_1}\left[\frac{X}{X+c/K^2} \right] \right| \geq c'.\]

Note that here we use \(X\) to denote a different random variable from the covariate in the maintexts. We first show the claim leads to the results in Lemma 5. With \(\omega_i\) (\(i=0,1\)) we define \(\widetilde{\omega}_i\) supported on \(\{0\} \cup [c/K^2,1]\) such that \[\widetilde{\omega}_i (dx) = \frac{c}{K^2x}\omega_i(dx) + \left(1-\mathbb{E}\left[\frac{c}{K^2X} \right] \right) \delta_0(dx).\] And for \(X \sim \widetilde{\omega}_i\), let \(\nu_i\) be the distribution of \(\frac{c_1 \log n}{n}X\). We let \(p \sim \nu_i\), \[\pi = \left\{ \begin{array}{rl} \epsilon \left(1+\frac{c_1 \log n}{np} \frac{c}{K^2} \right) & \text{if } p>0,\\ \epsilon & \text{if } p=0 \end{array} \right.\] and \(\mu = \epsilon /\pi\). Note that \(\pi, \mu\) are both functions of \(p\). We then verify the properties in Lemma 5 with joint distribution \(\mu_i\) defined above.

\(\text{supp}(\nu_i) \subseteq \{0\} \cup \frac{c_1 \log n}{n}[c/K^2,1]\) implies the range of \(p\).
For \(i,j,k \geq 0\) and \(1\leq i+ j +k \leq 3K\) we have \[\begin{align} &\, \mathbb{E}_{\mu_0}[p^i (p \pi)^j (p \pi \mu)^k] \\ = &\, \mathbb{E}_{\mu_0}[p^i (p \pi)^j (p \pi \mu)^k I(p>0)] \\ = &\, \mathbb{E}_{\nu_0}\left[ p^{i+k} \epsilon^{j+k} \left( p+\frac{cc_1 \log n}{nK^2} \right)^j I(p>0) \right] \\ = &\, \mathbb{E}_{\widetilde{\omega}_0} \left[ \left(\frac{c_1 \log n}{n} \right)^{i+j+k} \epsilon^{j+k} X^{i+k} \left(X + \frac{c}{K^2} \right)^j I(X>0) \right] \\ = &\, \int \left(\frac{c_1 \log n}{n} \right)^{i+j+k} \epsilon^{j+k} \frac{c}{K^2} x^{i+k-1} \left(x+ \frac{c}{K^2} \right)^j \omega_0 (dx) \\ = &\, \int \left(\frac{c_1 \log n}{n} \right)^{i+j+k} \epsilon^{j+k} \frac{c}{K^2} x^{i+k-1} \left(x+ \frac{c}{K^2} \right)^j \omega_1 (dx) \\ = &\, \mathbb{E}_{\mu_1}[p^i (p \pi)^j (p \pi \mu)^k]. \end{align}\] The first four equations follow from the definition of distributions above and the fifth follows from 32 .
\[\begin{align} & \, \mathbb{E}_{\mu_i}[p]\\ = &\, \int_{p>0} p \nu_i(dp) \\ = &\, \int_{x>0}\frac{c_1 \log n}{n} x \widetilde{\omega}_i (dx)\\ = &\, \int_{x>0}\frac{c_1 \log n}{n} \frac{c}{K^2} \widetilde{\omega}_i (dx)\\ = &\, \frac{c c_1}{c_2^2} \frac{1}{n \log n} \leq \frac{c c_1 c_3}{c_2^2 d} \leq \frac{1}{d} \end{align}\] as long as \(c c_1 c_3/c_2^2 \leq 1\).
\[\begin{align} &\, |\mathbb{E}_{\mu_0}[p \mu] - \mathbb{E}_{\mu_1}[p \mu]|\\ = &\, \left | \mathbb{E}_{\nu_0}\left[ \frac{p^2}{p+ \frac{c_1 \log n}{n}\frac{c}{K^2} } \right] - \mathbb{E}_{\nu_1}\left[ \frac{p^2}{p+ \frac{c_1 \log n}{n}\frac{c}{K^2} } \right]\right| \\ =&\, \frac{c_1 \log n}{n} \left | \mathbb{E}_{\widetilde{\omega}_0}\left[ \frac{X^2}{X+ \frac{c}{K^2} } \right] - \mathbb{E}_{\widetilde{\omega}_1}\left[ \frac{X^2}{X+ \frac{c}{K^2} } \right]\right|\\ = &\, \frac{c c_1 \log n}{n K^2} \left | \mathbb{E}_{{\omega}_0}\left[ \frac{X}{X+ \frac{c}{K^2} } \right] - \mathbb{E}_{{\omega}_1}\left[ \frac{X}{X+ \frac{c}{K^2} } \right]\right| \\ \geq &\, \frac{c c' c_1}{c_2^2 n \log n}. \end{align}\] And we define \(c_4 = cc' c_1/c_2^2\).

Thus we only need to prove the claim. By duality of polynomial approximation and moment matching (see e.g., Lemma 19 in [83]) it suffices to prove \[\inf_{P \in \text{span}\{1/x, 1, x, \dots, x^{3K}\} } \max_{x \in [c/K^2,1]} \left | \frac{x}{x+c/K^2}-P(x) \right| \geq c'/2.\] Let \[\label{eq:poly-approximation} E_n(f; S) = \min_{\text{deg}(P) \leq n} \max_{x \in S} |f(x)-P(x)|\tag{33}\] be the best polynomial approximation error (with polynomial’s degree smaller than \(n\)) of \(f\) on a set \(S\). We have the following lemma.

Lemma 7. There exists \(c_0, c' >0\) such that \[\liminf_{n \rightarrow \infty} \inf_{\alpha \in \mathbb{R}} E_n \left(\frac{x}{x+c_0n^{-2}} + \frac{\alpha}{x}; [c_0 n^{-2},1] \right) \geq c' >0.\]

Let \(c = c_0/9\), we then have for any \(P \in \text{span}\{1/x, 1, x,\dots, x^{3K} \}\), write \(P(x) = \frac{\alpha_{-1}}{x}+ \sum_{k=0}^{3K}\alpha_k x^k\) and let \(P_{\geq}(x) = \sum_{k=0}^{3K}\alpha_k x^k\) \[\begin{align} & \,\max_{x \in [c/K^2,1]} \left| \frac{x}{x+c/K^2} -P(x) \right| \\ \geq &\, \inf_{\alpha \in \mathbb{R}} \max_{x \in [c/K^2,1]} \left| \frac{x}{x+c/K^2} + \frac{\alpha}{x} -P_{\geq}(x) \right| \\ = &\, \inf_{\alpha \in \mathbb{R}} \max_{x \in [c_0/(3K)^2,1]} \left| \frac{x}{x+c_0/(3K)^2} + \frac{\alpha}{x} -P_{\geq}(x) \right| \\ \geq &\, \inf_{\alpha \in \mathbb{R}} E_{3K} \left( \frac{x}{x+c_0/(3K)^2} + \frac{\alpha}{x}; [c_0 / (3K)^2, 1] \right) \\ \geq &\, \frac{c'}{2}, \end{align}\] as \(K = c_2 \log n\) is large enough.

We then prove Lemma 7. By Theorem 7.2.4 in [52] (and use the notation there), we have \[\label{eq:moduli-smooth} \omega_{\varphi}^1 \left(f,\frac{1}{n} \right)_{\infty} \leq \frac{M}{n}\sum_{\ell=0}^n E_{\ell}(f;[0,1]),\tag{34}\] where \[\omega_{\varphi}^1 \left(f,\frac{1}{n} \right)_{\infty} = \sup_{0 < h \leq 1/n} \sup_{x,x+h\varphi(x) \in [0,1]} |f(x+h\varphi(x))-f(x)|, \, \varphi(x) = \sqrt{x(1-x)}\] and \(M\) is an absolute constant. Let \(f_{\alpha}(x) = \frac{x}{x+c_0 n^{-2}} + \frac{\alpha}{x}\) and \(\widetilde{f}_{\alpha}(x) = f_{\alpha}\left(c_0n^{-2} + (1-c_0n^{-2}) x\right), x\in[0,1]\). Denote \(m = \lceil n/\sqrt{c_0} \rceil\). Here \(c_0\) is a constant to be chosen. Consider two cases of \(\alpha\) separately:

Case 1: \(|\alpha| \leq c_0 / n^2\). By 34 we have \[\begin{align} &\,E_n(f_{\alpha}; [c_0n^{-2},1]) \\ = &\,E_n(\widetilde{f}_{\alpha}; [0,1]) \\ \geq &\, \frac{1}{m-n+1} \sum_{\ell =n}^m \mathbb{E}_{\ell} (\widetilde{f}_{\alpha}; [0,1]) \\ \geq &\, \frac{1}{m} \sum_{\ell =n}^m \mathbb{E}_{\ell} \left(\widetilde{f}_{\alpha}; [0,1] \right) \\ \geq &\, \frac{1}{M} \omega_{\varphi}^1 \left(\widetilde{f}_{\alpha},\frac{1}{m} \right)_{\infty} - \frac{1}{m} \sum_{\ell =0}^{n-1} \mathbb{E}_{\ell} \left(\widetilde{f}_{\alpha}; [0,1] \right) \\ \geq &\, \frac{1}{M} \omega_{\varphi}^1 \left(\widetilde{f}_{\alpha},\frac{1}{m} \right)_{\infty} - \frac{n}{m} \mathbb{E}_{0} \left(\widetilde{f}_{\alpha}; [0,1] \right) \\ \geq &\, \frac{1}{M} \omega_{\varphi}^1 \left(\widetilde{f}_{\alpha},\frac{1}{m} \right)_{\infty} - 2 \sqrt{c_0}, \end{align}\] where the first and the fourth inequality follow from the monotonity of \(E_n\), the third inequality applies 34 and the last one follows from \(|\widetilde{f}_{\alpha}| \leq 2\). For the first term, we have \[\begin{align} &\, \omega_{\varphi}^1 \left(\widetilde{f}_{\alpha},\frac{1}{m} \right)_{\infty} \\ \geq &\, \sup_{t_1, t_2 \in [3,4]} \left | \widetilde{f}_{\alpha}\left(\frac{t_1}{m^2} \right) -\widetilde{f}_{\alpha}\left(\frac{t_2}{m^2} \right) \right| \\ = &\, \sup_{t_1, t_2 \in [3,4]} \left | {f}_{\alpha}\left(\frac{c_0}{n^2} + \left(1-\frac{c_0}{n^2} \right) \frac{t_1}{m^2} \right) -{f}_{\alpha}\left(\frac{c_0}{n^2} + \left(1-\frac{c_0}{n^2} \right) \frac{t_2}{m^2} \right) \right|, \\ \end{align}\] where the inequality follows from the fact that when \(m\) is sufficiently large, \[\frac{t_2-t_1}{m^2} \leq \frac{1}{m} \sqrt{\frac{t_2}{m} \left(1-\frac{t_2}{m} \right)}\] holds for all \(t_1 < t_2\) and \(t_1,t_2 \in [3,4]\). Note that as \(m \rightarrow \infty\) \[1+ \left(\frac{n^2}{c_0}-1 \right)\frac{3}{m^2} \rightarrow 4,\] \[1+ \left(\frac{n^2}{c_0}-1 \right)\frac{4}{m^2} \rightarrow 5.\] Thus we have \[\begin{align} & \, \sup_{t_1, t_2 \in [3,4]} \left | {f}_{\alpha}\left(\frac{c_0}{n^2} + \left(1-\frac{c_0}{n^2} \right) \frac{t_1}{m^2} \right) -{f}_{\alpha}\left(\frac{c_0}{n^2} + \left(1-\frac{c_0}{n^2} \right) \frac{t_2}{m^2} \right) \right|\\ \geq &\, \sup_{t_1, t_2 \in [3.5, 5]} \left | {f}_{\alpha}\left(\frac{c_0}{n^2} t_1 \right) - {f}_{\alpha}\left(\frac{c_0}{n^2} t_2 \right)\right| \\ \geq &\, \sup_{t \in[3.5,4]} \left | {f}_{\alpha}\left(\frac{c_0}{n^2} \left(t+ \frac{1}{2} \right) \right) - {f}_{\alpha}\left(\frac{c_0}{n^2} t \right)\right| \\ = &\, \sup_{t \in[3.5,4]} \left | \frac{t+1/2}{t+3/2} - \frac{t}{t+1} + \frac{\alpha n^2}{c_0} \left(\frac{1}{t+1/2} -\frac{1}{t} \right) \right |\\ \geq &\, \inf_{\beta} \sup_{t \in[3.5,4]} \left | \frac{t+1/2}{t+3/2} - \frac{t}{t+1} + \beta \left(\frac{1}{t+1/2} -\frac{1}{t} \right) \right |\\ \geq &\, a_1, \end{align}\] where \(a_1\) is a positive constant independent of \(\alpha, n, c_0\) since \(\frac{t+1/2}{t+3/2} - \frac{t}{t+1}\) and \(\frac{1}{t+1/2} -\frac{1}{t}\) are linearly independent. Hence we have \[E_n(f_{\alpha}; [c_0n^{-2},1]) \geq \frac{a_1}{M}-2\sqrt{c_0}, \, |\alpha|\leq \frac{c_0}{n^2}.\] Case 2: \(|\alpha| >{c_0}/{n^2}\), similar to the proof in Case 1 we have \[\begin{align} &\, E_n(f_{\alpha}; [c_0n^{-2},1]) \\ \geq &\, \frac{1}{M} \omega_{\varphi}^1 \left(\widetilde{f}_{\alpha},\frac{1}{m} \right)_{\infty} - \sqrt{c_o} \mathbb{E}_{0} \left(\widetilde{f}_{\alpha}; [0,1] \right) \\ \geq &\, \frac{1}{M} \omega_{\varphi}^1 \left(\widetilde{f}_{\alpha},\frac{1}{m} \right)_{\infty} - \sqrt{c_o} \left(1+ \frac{|\alpha|n^2}{c_0} \right). \end{align}\] The second inequality follows from \(|f_{\alpha}| \leq \left(1 + \frac{|\alpha|n^2}{c_0} \right)\). And for the first term by the same argument in Case 1, \[\begin{align} &\, \omega_{\varphi}^1 \left(\widetilde{f}_{\alpha},\frac{1}{m} \right)_{\infty}\\ \geq &\, \sup_{t \in[3.5,4]} \left | {f}_{\alpha}\left(\frac{c_0}{n^2} \left(t+ \frac{1}{2} \right) \right) - {f}_{\alpha}\left(\frac{c_0}{n^2} t \right)\right| \\ = &\, \sup_{t \in[3.5,4]} \left | \frac{t+1/2}{t+3/2} - \frac{t}{t+1} + \frac{\alpha n^2}{c_0} \left(\frac{1}{t+1/2} -\frac{1}{t} \right) \right |\\ \geq &\, \frac{\alpha n^2}{c_0} \inf_{\beta \in \mathbb{R}} \sup_{t \in[3.5,4]} \left | \beta \left(\frac{t+1/2}{t+3/2} - \frac{t}{t+1}\right) + \frac{1}{t+1/2} -\frac{1}{t} \right |\\ \geq & \, \frac{\alpha n^2}{c_0} a_2, \end{align}\] where \(a_2\) is a positive constant independent of \(\alpha, n, c_0\) again since \(\frac{t+1/2}{t+3/2} - \frac{t}{t+1}\) and \(\frac{1}{t+1/2} -\frac{1}{t}\) are linearly independent. Hence we have for \(|\alpha|> c_0/n^2\) (choose \(c_0\) such that \(\frac{a_2}{M}-\sqrt{c_0}>0\)) \[\begin{align} &\, E_n(f_{\alpha}; [c_0n^{-2},1]) \\ \geq &\, \frac{\alpha n^2}{c_0 M}a_2 -\sqrt{c_o} \left(1+ \frac{|\alpha|n^2}{c_0} \right) \\ \geq &\, \inf_{|\alpha|> c_0/n^2} \frac{|\alpha|n^2}{c_0} \left(\frac{a_2}{M}-\sqrt{c_0} \right) -\sqrt{c_0} \\ \geq & \, \frac{a_2}{M}-2\sqrt{c_0}. \end{align}\] Combining Case 1 and Case 2, we conclude for any \(\alpha \in \mathbb{R}\) \[E_n(f_{\alpha}; [c_0n^{-2},1]) \geq \min \left\{ \frac{a_1}{M} - 2\sqrt{c_0}, \frac{a_2}{M} - 2\sqrt{c_0}\right\}.\] Choosing \(c_0>0\) sufficiently small completes the proof. ◻

12.4 Proof of Lemma 6 ↩︎

Proof. It is known that the conditional distribution of \(Y^n\) given the summation of its components is multinomial: \[\label{eq:poisson-multinomial} \mathbf{Y}^n \mid \sum_{k=1}^d Y_k^n = m \stackrel{d}{=} \mathbf{X}^m\tag{35}\] Hence we have \[\begin{align} & \, \mathbb{E}[f(Y_1^n,\dots, Y_d^n)]\\ \geq&\, \sum_{m=0}^n \mathbb{E}\left[f(Y_1^n,\dots, Y_d^n) | \sum_{k=1}^d Y_k^n = m\right] \mathbb{P}\left(\sum_{k=1}^d Y_k^n = m\right) \\ = &\, \sum_{m=0}^n \mathbb{E}\left[f(X_1^m,\dots, X_d^m)\right] \mathbb{P}\left(\sum_{k=1}^d Y_k^n = m\right) \\ \geq &\, \mathbb{E}\left[f(X_1^n,\dots, X_d^n)\right] \mathbb{P}\left(\sum_{k=1}^d Y_k^n \leq n\right) \\ \geq &\, \frac{1}{2} \mathbb{E}\left[f(X_1^n,\dots, X_d^n)\right] \end{align}\] where the first inequality follows from truncating the summation, the first equation follows from 35 , the second inequality follows from monotonicity of \(\mathbb{E}[f(X_1^n,\dots, X_d^n)]\) and the last inequality follows from \(\sum_{k=1}^d Y_k^n \sim \text{Poisson}(n)\) and \(\mathbb{P}(\sum_{k=1}^d Y_k^n \leq n) \geq 1/2\). ◻

References↩︎

[1]

D. B. Rubin, “Using multivariate matched sampling and regression adjustment to control bias in observational studies,” Journal of the American Statistical Association, vol. 74, no. 366a, pp. 318–328, 1979.

[2]

M. A. Hernán and J. M. Robins, “Causal inference.” CRC Boca Raton, FL, 2010.

[3]

P. R. Rosenbaum and D. B. Rubin, “The central role of the propensity score in observational studies for causal effects,” Biometrika, vol. 70, no. 1, pp. 41–55, 1983.

[4]

J. Hahn, “On the role of the propensity score in efficient semiparametric estimation of average treatment effects,” Econometrica, pp. 315–331, 1998.

[5]

K. Hirano, G. W. Imbens, and G. Ridder, “Efficient estimation of average treatment effects using the estimated propensity score,” Econometrica, vol. 71, no. 4, pp. 1161–1189, 2003.

[6]

M. J. van der Laan, J. M. Robins, M. J. van der Laan, and J. M. Robins, “Unified approach for causal inference and censored data,” Unified Methods for Censored Longitudinal Data and Causality, pp. 311–370, 2003.

[7]

V. Chernozhukov et al., “Double/debiased machine learning for treatment and structural parameters.” Oxford University Press Oxford, UK, 2018.

[8]

E. H. Kennedy, “Semiparametric doubly robust targeted double machine learning: A review,” arXiv preprint arXiv:2203.06469, 2022.

[9]

J. Robins, L. Li, E. Tchetgen, A. van der Vaart, et al., “Higher order influence functions and minimax estimation of nonlinear functionals,” in Probability and statistics: Essays in honor of david a. freedman, vol. 2, Institute of Mathematical Statistics, 2008, pp. 335–422.

[10]

W. Zheng and M. J. van der Laan, “Asymptotic theory for cross-validated targeted maximum likelihood estimation,” 2010.

[11]

S. Yadlowsky, “Explaining practical differences between treatment effect estimators with high dimensional asymptotics,” arXiv preprint arXiv:2203.12538, 2022.

[12]

K. Jiang, R. Mukherjee, S. Sen, and P. Sur, “A new central limit theorem for the augmented IPW estimator: Variance inflation, cross-fit covariance and beyond,” arXiv preprint arXiv:2205.10198, 2022.

[13]

M. Celentano and M. J. Wainwright, “Challenges of the inconsistency regime: Novel debiasing methods for missing data models,” arXiv preprint arXiv:2309.01362, 2023.

[14]

J. Robins, L. Li, E. Tchetgen, and A. W. van der Vaart, “Quadratic semiparametric von mises calculus,” Metrika, vol. 69, pp. 227–247, 2009.

[15]

E. H. Kennedy, “Towards optimal doubly robust estimation of heterogeneous causal effects,” Electronic Journal of Statistics, vol. 17, no. 2, pp. 3008–3049, 2023, doi: 10.1214/23-EJS2157.

[16]

V. Semenova and V. Chernozhukov, “Debiased machine learning of conditional average treatment effects and other causal functions,” The Econometrics Journal, vol. 24, no. 2, pp. 264–289, 2021.

[17]

A. Belloni, V. Chernozhukov, I. Fernandez-Val, and C. Hansen, “Program evaluation and causal inference with high-dimensional data,” Econometrica, vol. 85, no. 1, pp. 233–298, 2017.

[18]

S. Athey, “The impact of machine learning on economics,” in The economics of artificial intelligence: An agenda, University of Chicago Press, 2018, pp. 507–547.

[19]

J. Bradic, V. Chernozhukov, W. K. Newey, and Y. Zhu, “Minimax semiparametric learning with approximate sparsity,” arXiv preprint arXiv:1912.12213, 2019.

[20]

M. H. Maathuis, M. Kalisch, and P. Bühlmann, “Estimating high-dimensional intervention effects from observational data,” The Annals of Statistics, vol. 37, no. 6A, pp. 3133–3164, 2009.

[21]

W. Lin, “Agnostic notes on regression adjustments to experimental data: Reexamining Freedman’s critique,” The Annals of Applied Statistics, vol. 7, no. 1, pp. 295–318, 2013, doi: 10.1214/12-AOAS583.

[22]

Q. Zhao, “Topics in causal and high dimensional inference,” PhD thesis, Stanford University, 2016.

[23]

A. Chakrabortty, P. Nandy, and H. Li, “Inference for individual mediation effects and interventional effects in sparse high-dimensional causal graphical models,” arXiv preprint arXiv:1809.10652, 2018.

[24]

S. Ma, L. Zhu, Z. Zhang, C.-L. Tsai, and R. J. Carroll, “A robust and efficient approach to causal inference based on sparse sufficient dimension reduction,” Annals of statistics, vol. 47, no. 3, p. 1505, 2019.

[25]

L. Lei and P. Ding, “Regression adjustment in completely randomized experiments with a diverging number of covariates,” Biometrika, vol. 108, no. 4, pp. 815–828, 2021.

[26]

J. Antonelli, G. Papadogeorgou, and F. Dominici, “Causal inference in high dimensions: A marriage between bayesian modeling and good frequentist properties,” Biometrics, vol. 78, no. 1, pp. 100–114, 2022.

[27]

D. Tang, D. Kong, W. Pan, and L. Wang, “Ultra-high dimensional variable selection for doubly robust causal inference,” Biometrics, vol. 79, no. 2, pp. 903–914, 2023.

[28]

J.-H. Du, Z. Zeng, E. H. Kennedy, L. Wasserman, and K. Roeder, “Causal inference for genomic data with multiple heterogeneous outcomes,” arXiv preprint arXiv:2404.09119, 2024.

[29]

G. Chamberlain, “Efficiency bounds for semiparametric regression,” Econometrica: Journal of the Econometric Society, pp. 567–596, 1992.

[30]

R. Kleinberg, “Nearly tight bounds for the continuum-armed bandit problem,” Advances in Neural Information Processing Systems, vol. 17, 2004.

[31]

S. Balakrishnan and L. Wasserman, “HYPOTHESIS TESTING FOR DENSITIES AND HIGH-DIMENSIONAL MULTINOMIALS,” The Annals of Statistics, vol. 47, no. 4, pp. 1893–1927, 2019.

[32]

W. H. Organization, International statistical classification of diseases and related health problems: Alphabetical index, vol. 3. World Health Organization, 2004.

[33]

L. Yao, S. Li, Y. Li, H. Xue, J. Gao, and A. Zhang, “On the estimation of treatment effect with text covariates,” in International joint conference on artificial intelligence, 2019.

[34]

K. A. Keith, D. Jensen, and B. O’Connor, “Text and causal inference: A review of using text to remove confounding from causal estimates,” arXiv preprint arXiv:2005.00649, 2020.

[35]

L. Paninski, “Estimation of entropy and mutual information,” Neural computation, vol. 15, no. 6, pp. 1191–1253, 2003.

[36]

G. Valiant and P. Valiant, “A CLT and tight lower bounds for estimating entropy.” in Electron. Colloquium comput. complex., 2010, vol. 17, p. 179.

[37]

J. Jiao, K. Venkat, Y. Han, and T. Weissman, “Minimax estimation of functionals of discrete distributions,” IEEE Transactions on Information Theory, vol. 61, no. 5, pp. 2835–2885, 2015.

[38]

Y. Wu and P. Yang, “Minimax rates of entropy estimation on large alphabets via best polynomial approximation,” IEEE Transactions on Information Theory, vol. 62, no. 6, pp. 3702–3720, 2016.

[39]

Y. Wu and P. Yang, “Chebyshev polynomials, moment matching, and optimal estimation of the unseen,” The Annals of Statistics, vol. 47, no. 2, pp. 857–883, 2019.

[40]

A. F. Timan, Theory of approximation of functions of a real variable. Elsevier, 2014.

[41]

D. G. Luenberger, Optimization by vector space methods. John Wiley & Sons, 1997.

[42]

D. B. Rubin, “Estimating causal effects of treatments in randomized and nonrandomized studies.” Journal of educational Psychology, vol. 66, no. 5, p. 688, 1974.

[43]

J. Splawa-Neyman, D. M. Dabrowska, and T. P. Speed, “On the application of probability theory to agricultural experiments. Essay on principles. Section 9.” Statistical Science, pp. 465–472, 1990.

[44]

A. D’Amour, P. Ding, A. Feller, L. Lei, and J. Sekhon, “Overlap in observational studies with high-dimensional covariates,” Journal of Econometrics, vol. 221, no. 2, pp. 644–654, 2021.

[45]

J. M. Robins, A. Rotnitzky, and L. P. Zhao, “Estimation of regression coefficients when some regressors are not always observed,” Journal of the American statistical Association, vol. 89, no. 427, pp. 846–866, 1994.

[46]

D. O. Scharfstein, A. Rotnitzky, and J. M. Robins, “Adjusting for nonignorable drop-out using semiparametric nonresponse models,” Journal of the American Statistical Association, vol. 94, no. 448, pp. 1096–1120, 1999.

[47]

C. Wang, S.-M. Lee, and E. C. Chao, “Numerical equivalence of imputing scores and weighted estimators in regression analysis with missing covariates,” Biostatistics, vol. 8, no. 2, pp. 468–473, 2007.

[48]

T. Słoczyński, S. D. Uysal, and J. M. Wooldridge, “Covariate balancing and the equivalence of weighting and doubly robust estimators of average treatment effects,” arXiv preprint arXiv:2310.18563, 2023.

[49]

A. B. Tsybakov, Introduction to nonparametric estimation, 1st ed. 2009. New York, NY: Springer New York, 2009.

[50]

O. Lepski, A. Nemirovski, and V. Spokoiny, “On estimation of the l r norm of a regression function,” Probability theory and related fields, vol. 113, pp. 221–253, 1999.

[51]

J. R. Rice, “Tchebycheff approximation in several variables,” Transactions of the American Mathematical Society, vol. 109, no. 3, pp. 444–466, 1963.

[52]

Z. Ditzian and V. Totik, Moduli of smoothness, vol. 9. Springer Science & Business Media, 2012.

[53]

V. Totik, Polynomial approximation on polytopes, vol. 232. American Mathematical Society, 2014.

[54]

M. L. Clevenson and W. Watkins, “Majorization and the birthday inequality,” Mathematics Magazine, vol. 64, no. 3, pp. 183–188, 1991.

[55]

M. C. Wendl, “Collision probability between sets of random variables,” Statistics & probability letters, vol. 64, no. 3, pp. 249–254, 2003.

[56]

T. Nakata, “The number of collisions for the occupancy problem with unequal probabilities,” Advances in Applied Probability, vol. 46, no. 1, pp. 168–185, 2014.

[57]

Z. Zeng, E. H. Kennedy, L. M. Bodnar, and A. I. Naimi, “Efficient generalization and transportation,” arXiv preprint arXiv:2302.00092, 2023.

[58]

R. K. Crump, V. J. Hotz, G. Imbens, and O. Mitnik, “Moving the goalposts: Addressing limited overlap in the estimation of average treatment effects by changing the estimand.” National Bureau of Economic Research Cambridge, Mass., USA, 2006.

[59]

J. M. Robins, S. D. Mark, and W. K. Newey, “Estimating exposure effects by modelling the expectation of exposure conditional on confounders,” Biometrics, pp. 479–495, 1992.

[60]

S. Vansteelandt and R. M. Daniel, “On regression adjustment for the propensity score,” Statistics in medicine, vol. 33, no. 23, pp. 4053–4072, 2014.

[61]

X. Zhou and A. Opacic, “Marginal interventional effects,” arXiv preprint arXiv:2206.10717, 2022.

[62]

E. H. Kennedy, “Nonparametric causal effects based on incremental propensity score interventions,” Journal of the American Statistical Association, vol. 114, no. 526, pp. 645–656, 2019.

[63]

R. D. Shah and J. Peters, “The hardness of conditional independence testing and the generalised covariance measure,” The Annals of Statistics, vol. 48, no. 3, pp. 1514–1538, 2020, doi: 10.1214/19-AOS1857.

[64]

I. Dı́az, “Non-agency interventions for causal mediation in the presence of intermediate confounding,” Journal of the Royal Statistical Society Series B: Statistical Methodology, p. qkad130, 2023.

[65]

A. Chakrabortty, G. Dai, and E. T. Tchetgen, “A general framework for treatment effect estimation in semi-supervised and high dimensional settings,” arXiv preprint arXiv:2201.00468, 2022.

[66]

Y. Zhang, A. Chakrabortty, and J. Bradic, “Semi-supervised causal inference: Generalizable and double robust inference for average treatment effects under selection bias with decaying overlap,” arXiv preprint arXiv:2305.12789, 2023.

[67]

N. Kallus and X. Mao, “On the role of surrogates in the efficient estimation of treatment effects with limited outcome data,” arXiv preprint arXiv:2003.12408, 2020.

[68]

Z. Zeng et al., “Continuous treatment effects with surrogate outcomes,” arXiv preprint arXiv:2402.00168, 2024.

[69]

Y. He, B. Meng, Z. Zeng, and G. Xu, “On the phase transition of wilks’ phenomenon,” Biometrika, vol. 108, no. 3, pp. 741–748, 2021.

[70]

V. Chernozhukov and C. Hansen, “The effects of 401 (k) participation on the wealth distribution: An instrumental quantile regression analysis,” Review of Economics and statistics, vol. 86, no. 3, pp. 735–751, 2004.

[71]

J. M. Poterba and S. F. Venti, “401 (k) plans and tax-deferred saving,” in Studies in the economics of aging, University of Chicago Press, 1994, pp. 105–142.

[72]

J. M. Poterba, S. F. Venti, and D. A. Wise, “Do 401 (k) contributions crowd out other personal saving?” Journal of Public Economics, vol. 58, no. 1, pp. 1–32, 1995.

[73]

N. E. Karoui and E. Purdom, “Can we trust the bootstrap in high-dimensions? The case of linear models,” Journal of Machine Learning Research, vol. 19, no. 5, pp. 1–66, 2018, [Online]. Available: http://jmlr.org/papers/v19/17-006.html.

[74]

T. D. Hocking, A. Joulin, F. Bach, and J.-P. Vert, “Clusterpath an algorithm for clustering using convex fusion penalties,” in 28th international conference on machine learning, 2011, p. 1.

[75]

E. C. Chi and K. Lange, “Splitting methods for convex clustering,” Journal of Computational and Graphical Statistics, vol. 24, no. 4, pp. 994–1013, 2015.

[76]

K. M. Tan and D. Witten, “Statistical properties of convex clustering,” Electronic journal of statistics, vol. 9, no. 2, p. 2324, 2015.

[77]

P. Radchenko and G. Mukherjee, “Convex clustering via l 1 fusion penalization,” Journal of the Royal Statistical Society Series B: Statistical Methodology, vol. 79, no. 5, pp. 1527–1546, 2017.

[78]

D. Sun, K.-C. Toh, and Y. Yuan, “Convex clustering: Model, theoretical guarantee and efficient algorithm,” Journal of Machine Learning Research, vol. 22, no. 9, pp. 1–32, 2021.

[79]

L. Devroye, L. Györfi, and G. Lugosi, A probabilistic theory of pattern recognition, vol. 31. Springer Science & Business Media, 2013.

[80]

L. Le Cam, Asymptotic methods in statistical decision theory. Springer Science & Business Media, 2012.

[81]

M. Mitzenmacher and E. Upfal, Probability and computing: Randomization and probabilistic techniques in algorithms and data analysis. Cambridge university press, 2017.

[82]

A. W. Van der Vaart, Asymptotic statistics, vol. 3. Cambridge university press, 2000.

[83]

Y. Han, J. Jiao, T. Weissman, and Y. Wu, “Optimal rates of entropy estimation over lipschitz balls,” Annals of Statistics, vol. 48, no. 6, 2020.

Causal inference with high-dimensional discrete covariates

Abstract

1 Introduction↩︎

2 Setup and Notation↩︎

2.1 Data Generating Process↩︎

2.2 Causal Assumptions and Other Notation↩︎

3 Average Treatment Effects↩︎

3.1 Upper Bound↩︎

3.2 Minimax Lower Bound↩︎

4 Role of Effect Homogeneity↩︎

4.1 Upper Bound↩︎

4.2 Minimax Lower Bound↩︎

5 Prior Knowledge of Covariate Distribution↩︎

5.1 Average Treatment Effects↩︎

5.2 Variance-weighted Average Treatment Effects↩︎

6 Simulation Study↩︎

7 Real Data Analysis↩︎

8 Discussion↩︎

9 Asymptotic Normality of the Plug-in Estimator↩︎

10 Additional Simulation Results↩︎

10.1 Asymptotic Normality↩︎

10.2 Bootstrap Variance Estimation↩︎

11 Proof of Main Results↩︎

11.1 Proof of Proposition 1↩︎

11.2 Proof of Proposition 2↩︎

11.3 Proof of Proposition 3↩︎

11.4 Proof of Theorem 1↩︎

11.5 Proof of Theorem 2↩︎

11.6 Proof of Lemma 1↩︎

11.7 Proof of Theorem 3↩︎

11.8 Proof of Theorem 4↩︎

11.9 Proof of Theorem 5↩︎

11.10 Proof of Theorem 6↩︎

11.11 Proof of Theorem 7↩︎

11.12 Proof of Theorem 8↩︎