Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
121,227 characters · 30 sections · 80 citation commands
Uniform Estimation and Inference for Nonparametric Partitioning-Based M-Estimators
\pagenumbering{roman}
Keywords: nonparametric estimation and inference, series methods, partitioning estimators, quantile regression, nonlinear regression, robust regression, generalized linear models, uniform distribution theory.
\thispagestyle{empty}
\setcounter{tocdepth}{2}
\ifthenelse{\boolean{standalone}}{
\listoftheorems[ title={List of important statements}, ignore={lemma,remark,remark*}, numwidth=3em, ]
\listoftheorems[ title={List of lemmas}, ignoreall, show={lem,lemma}, numwidth=3em, ] }
\pagenumbering{arabic}
\pagestyle{plain}
Let $(y_1,\bx_1),(y_2,\bx_2),\cdots,(y_n,\bx_n)$ be independent and identically distributed copies of the random vector $(Y,\bX)\in\mathcal{Y}\times\mathcal{X} \subseteq \mathbb{R}\times\mathbb{R}^d$. Given a loss function $\rho\colon \mathcal{Y}\times \mathcal{E} \times\mathcal{Q} \to \mathbb{R}$ with $\mathcal{E}\subseteq\mathbb{R}$ an open connected set and $\mathcal{Q}\subseteq\mathbb{R}$ a connected compact set, and $\eta\colon \mathbb{R} \to \mathcal{E}$ a strictly monotonic transformation function, consider the functional parameter $\mu_{0}\colon \mathcal{X}\times\mathcal{Q} \to \mathbb{R}$ satisfying
where the minimization is over the space of measurable functions from $\mathcal{X}$ to $\mathbb{R}$. In particular, we assume that the (local) minimum is achieved, which is true in most cases. This setup covers settings of interest in nonparametric statistics, econometrics, and data science, including generalized linear models, robust nonlinear regression, and generalized conditional quantile regression. In practice, the parameter of interest may be $\mu_{0}$ itself, or otherwise specific transformations thereof such as $\eta(\mu_0(\cdot, \cdot))$ or its partial derivatives. This paper presents uniform over $\mathcal{X}\times\mathcal{Q}$ estimation and inference results for $\mu_{0}$, and transformations thereof, based on nonparametric partitioning-based $M$-estimation.
The series (or sieve) nonparametric partitioning-based $M$-estimator is
where $\mathcal{B} \subseteq \mathbb{R}^K$ is the feasible set of the optimization problem, and $\bx \mapsto \bp(\bx) = \bp(\bx; \Delta, m) = \big( p_1(\bx; \Delta, m), \ldots, p_K(\bx; \Delta, m) \big)\trans$ is a dictionary of $K$ locally supported basis functions of order $m$ based on a quasi-uniform partition $\Delta = \left\{ \delta_l: 1 \leq l \leq \kappa \right\}$ containing a collection of open disjoint polyhedra in $\mathcal{X}$ such that the closure of their union covers $\mathcal{X}$.
Gyorfi-etal_2002_book give a textbook introduction to the partitioning-based estimation literature. In this nonparametric framework, the Haar basis
corresponds to the canonical basis with $m=1$ and $K=\kappa$, and is an essential building block for the construction of other basis functions. Since the Haar basis is “unconnected” across cells (i.e., each basis function is supported on a single cell), the resulting estimator in (ref) reduces to $K$ separate M-estimators, each only using observations with $\bx_i\in\delta_k$, for $k=1,\ldots,K$. For estimation and inference, “small” cells decrease bias but increase variance, while “large” cells have the opposite effect.
A natural generalization is the piecewise polynomial basis
where the vector $\br_m(\bx)$ contains the unique terms of an $(m-1)$th-degree polynomial expansion based on $\bx$ and $\otimes$ is the Kronecker product operator, and thus $K = \frac{(m+d-1)!}{(m-1)!d!}\kappa$. The resulting piecewise polynomial fit within each cell gives more flexible approximation, thus decreasing bias, but the estimation approach remains unconnected.
Since the piecewise polynomial estimator may be discontinuous over $\mathcal{X}$, it is sometimes preferred to impose smoothness restrictions across cells: for example, if the partition $\Delta$ admits a tensor product representation with equal number of partitions along the $d$ axes, then the Splines basis is
where $\be_k$ denotes the $k$-th unit vector $(1\leq k\leq d)$, and $\bT_{s}$ denotes a transformation matrix that ensures the estimator $\bx\mapsto\widehat{\mu}(\bx, q)$ is $(s-1)$-times continuously differentiable ($s<m$) over $\mathcal{X}$, and thus $K=((m-s)\kappa^{1/d}+s)^d$. Due to the global smoothness restrictions, the spline basis is no longer unconnected, and the resulting estimator in (ref) cannot be reduced to separate local estimators. Other spline constructions on more general partitioning schemes are available, and compactly supported wavelets are yet another example of a local basis constructed recursively out of the Haar basis. See Belloni-Chernozhukov-Chetverikov-Kato_2015_JoE, Cattaneo-Farrell_2013_JoE, Cattaneo-Farrell-Feng_2020_AOS, and Chen-Christensen_2015_JOE for more discussion on these and other basis of approximation. Furthermore, partitioning-based estimation naturally arises in the recursive partitioning literature Devroye-etal2013_book,Zhang-Singer_2010_Book.
To enable good statistical performance, we need to restrict the partition of $\mathcal{X}$, and the local basis constructed on it. The first assumption concerns the regularity of the cells in the partition. Let $a_n\lesssim b_n$ denote $\limsup_{n\to\infty} |a_n/b_n|<\infty$.
Assumption (ref) requires the partition $\Delta$ be quasi-uniform: the elements in the partition $\Delta$ do not differ too much in size asymptotically. As a consequence, we can use the maximum diameter $h$ as a universal measure of mesh sizes. The next assumption requires the basis be “locally supported”, non-collinear, and bounded in a proper sense. A function $p(\cdot)$ on $\mathcal{X}$ is said to be active on $\delta\in\Delta$ if it is not identically zero on $\delta$; we also employ standard multi-index notation (see Section (ref) for details).
In Assumption (ref), condition (i) implies that each basis function in $\bp(\bx)$ is supported by a region consisting of a finite number of cells in $\Delta$ (independent of $n$). Then, as $\kappa\to\infty$, all basis functions are locally supported relative to the whole support of the data. Condition (ii) can be read as “non-collinearity” of the basis functions in $\bp(\bx)$. Since local support condition has been imposed, it suffices to require the basis functions are not too collinear “locally”. Condition (iii) controls the magnitude of the local basis in a uniform sense.
Assumptions (ref) and (ref) implicitly relate the number of approximating series terms, the number of cells in $\Delta$, and the maximum mesh size: $K\asymp \kappa\asymp h^{-d}$, where $a_n\asymp b_n$ means $a_n\lesssim b_n$ and $b_n\lesssim a_n$. Under appropriate assumptions on the statistical model (Assumption (ref) in Section (ref)), the parameter $m$ will control how well $\mu_0$ can be approximated by linear combinations of the local basis (via Assumption (ref) in Section (ref)). We consider large sample approximations where $d$ and $m$ are fixed constants, and $\kappa\to\infty$ (and thus $K\asymp h^{-d}\to\infty$) as $n\to\infty$. As a consequence, appropriate choices of $\Delta$ and $\bp(\cdot)$ will enable valid nonparametric approximations of $\mu_{0}$, and transformations thereof, in large samples.
In practice, the parameter (ref) and its associated plug-in estimator (ref) may not be unique (e.g., when the objective function is not convex and hence several local minima may exist). In such cases the interpretation of the estimator and its probability limit may depend on the specific (algorithmic) implementation used. This paper does not study these additional complications, but rather assumes that the estimator (ref) has been computed, and then relies on the assumptions in Section (ref) concerning the data generating process to study the large sample statistical properties of the partitioning-based $M$-estimator and transformations thereof.
The objective function in (ref) may not be convex. To address this challenge, we first provide primitive conditions for uniform over $\mathcal{X}\times\mathcal{Q}$ (and mean square) consistency of the partitioning-based estimator $\widehat{\mu}(\bx, q)$, taking explicitly into account whether the loss function $\theta\mapsto\rho(y,\eta(\theta);q)$ is convex: setting $\mathcal{B}=\mathbb{R}^K$ if it is convex, or otherwise $\mathcal{B}=\{\bb\in\mathbb{R}^K:\|\bb\|_\infty\leq R\}$ for some large enough fixed constant $R>0$, we establish $\sup_{q \in \mathcal{Q}}\big\| \widehat{\boldsymbol{\beta}}(q) - \bbeta_0(q) \big\|_{\infty} = o_{\P}(1)$, where $\|\cdot\|_{\infty}$ denotes the $\ell^\infty$-norm, and $\bbeta_0\colon \mathcal{Q} \to \mathbb{R}^K$ denotes coefficients such that $\bbeta_0(q)\trans \bp$ approximates $\mu_{0}$ well enough uniformly over $\mathcal{X}\times\mathcal{Q}$ (Assumption (ref) in Section (ref)). For the non-convex case, the resulting “fixed box” constrained optimization is arguably a mild assumption in practice, and may be justified in theory under different regularity conditions. These results are presented in Section (ref).
Taking the uniform consistency of the partitioning-based estimator as given, and hence being agnostic about the shape of the objective function and other optimization-related aspects, we establish three theoretical results for the partitioning-based series $M$-estimator in (ref):
These results allow for a large class of possibly non-smooth loss functions. In addition, we precisely characterize how the degrees of smoothness of $\rho$ and $\eta$ affect the order of the remainder in the uniform Bahadur representation for $\widehat{\mu}$, its convergence rates, and the validity of the associated uniform inference procedures. Results (i) and (ii) are presented in Section (ref), while results (iii) are presented in Sections (ref) and (ref).
Section (ref) introduces four examples: Generalized Conditional Quantile Regression, Generalized Conditional Distribution Regression, Generalized $L_p$ Regression Estimation, and Maximum Likelihood Logistic Regression. These examples are used to both motivate our high-level assumptions and demonstrate the broad applicability of our uniform estimation and inference results. Our most general results cover other applications such as nonparametric partitioning-based (quasi-maximum likelihood) Poisson regression, censored and truncated regression, as well as Tukey and Huber regression. Section (ref) presents the slightly simplified high-level technical assumptions used throughout the paper, but their most general form is given in the supplemental appendix to streamline the presentation. Section (ref) demonstrates how our general sufficient conditions are verified for each of our motivating examples.
Section (ref) discusses how our results can be extended to cover other parameters of interest, while Section (ref) concludes. The supplemental appendix reports simulation evidence, collects all the technical proofs, presents other theoretical results that may be of independent interest. In particular, our more general theoretical results (i) allow for $\mathcal{Q}$ to be a set of vectors rather than scalars, which can be useful in other examples beyond those studied in this paper; and (ii) consider more complex (VC-type) classes of loss and transformation functions, thereby covering a broader class of settings than those studied herein, but at the cost of additional, cumbersome notation and technicalities. In addition, the supplemental appendix presents new strong approximation results for a class of $K$-dimensional linear stochastic processes indexed by $\mathcal{X}\times\mathcal{Q}$ under standard complexity and smoothness conditions, leveraging a conditional Strassen’s Theorem chen2020jackknife,monrad1991nearby and generalizing prior Yurinskii's coupling results in the literature Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE,yurinskii1978error.
Our paper contributes to the literature on nonparametric curve estimation and inference, focusing in particular on series (or sieve) partitioning-based methods. See, for example, Eggermont-LaRiccia_2009_Book and Gyorfi-etal_2002_book for textbook introductions. This literature is mature and well-developed for the special case of a square loss function $\rho(y,\eta(\theta);q) = (y - \theta)^2$ with identity transformation $\eta(u)=u$, and hence not a function of $q\in\mathcal{Q}$. See, for example, Belloni-Chernozhukov-Chetverikov-Kato_2015_JoE, Cattaneo-Farrell_2013_JoE, Cattaneo-Farrell-Feng_2020_AOS, Cattaneo-Crump-Farrell-Feng_2024_AER, Chen-Christensen_2015_JOE, Huang_2003_AoS, Zhou-Shen-Wolfe_1998_AoS for pointwise and uniform over $\mathcal{X}$ estimation and inference results at different levels of generality, and with increasingly weaker technical conditions. This strand of the literature explicitly exploits the special structure, which leads to a closed-form solution of the estimator in (ref), and hence results are often obtained under minimal assumptions and technical regularity conditions. To be more precise, up to $\polylog(n)$ terms and mild regularity conditions, Cattaneo-Farrell-Feng_2020_AOS show that the minimal requirement $K/n\to0$ is (necessary and) sufficient for rate-optimal convergence rates for any $d\geq1$, and for strong approximations uniformly over $\mathcal{X}$ when $d=1$. They also establish valid strong approximations uniformly over $\mathcal{X}$ for $d>1$ under the requirement $K^3/n\to0$, up to $\polylog(n)$ terms and mild regularity conditions.
Despite aiming for generality, that is, allowing for a large class of loss functions with different levels of smoothness and a non-identity transformation function, this paper establishes rate-optimal uniform over both $\mathcal{X}$ and $\mathcal{Q}$ estimation results under the same minimal assumption $K/n\to0$ for unconnected bases, and under the slightly stronger assumption $K^2/n\to0$ for general partitioning-based estimators. Furthermore, we establish valid uniform over both $\mathcal{X}$ and $\mathcal{Q}$ inference under the same condition $K^3/n\to0$, leveraging a new strong approximation result given in the supplemental appendix. Compared to the prior literature focusing on the special case of square loss and identity transformation, we are able to achieve the same best known (in some cases rate-optimal) estimation and inference results, under the same (in some cases minimal) side rate restrictions and conditions on the partitioning-based method (Assumptions (ref) and (ref)). Furthermore, our results on uniform consistency disentangling convex and non-convex loss functions (Section (ref)), rate-optimal uniform Bahadur representation and convergence capturing explicitly the smoothness degree of the loss function (Section (ref)), and uniform feasible inference validity (Sections (ref) and (ref)), are necessarily new relative to prior work studying the special case of least square partitioning-based methods.
Going beyond square loss and identity transformation, there are only a handful of results available in the literature. The closest antecedent is Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE, who consider nonparametric conditional quantile series regression estimation and inference uniformly over $\mathcal{X}\times\mathcal{Q}$ with $\eta(u)=u$, and under the side rate restriction $K^4/n\to0$, up to $\polylog(n)$ terms, and other regularity conditions. As a comparison, for the special case of nonparametric quantile regression (Example (ref) below), this paper allows for a non-identity (inverse) link function $\eta(\cdot)$, and establishes convergence rates under the minimal condition $K/n\to0$ for piecewise polynomials, and the improved condition $K^2/n\to0$ for connected bases, while for uniform inference we require the weaker condition $K^3/n\to0$, in all cases up to $\polylog(n)$ terms. We also weaken other assumptions imposed in Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE: see Example (ref) in Sections (ref) and (ref), and Sections (ref) and (ref). On the other hand, it is worth noting that Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE also consider generic dimension-increasing covariates, while our paper focuses exclusively on partitioning-based local basis.
Our contributions can also be compared to recent work on nonparametric M-estimation employing other smoothing techniques. For example, Kong-Linton-Xia_2013_ET considers local polynomial methods, and Shang-Cheng_2013_AOS considers smoothing spline methods. Sections (ref) and (ref) give a more detailed comparison with prior literature, and explain precisely how our general results are either on par with or improve upon prior work. In a nutshell, we present estimation and inference results for partitioning-based $M$-estimators that (i) allow for a large class of possibly non-smooth loss functions, (ii) are uniformly valid over both $\mathcal{X}$ and $\mathcal{Q}$, (iii) achieve the best known (in some cases rate-optimal) convergence rates, and (iv) require substantially weaker (in some cases minimal) side rate restrictions and regularity conditions. Our results, for example, permit the use of Haar $M$-estimators for estimation and inference, which were ruled out by prior work.
We employ standard notation in probability, statistics and empirical process theory bhatia2013matrix,dudley2014uniform,Kallenberg2021,vaart96_weak_conver_empir_proces. For any vector $\ba=(a_1, \cdots, a_M)\in\mathbb{R}^M$, we write $\|\ba\|=(\sum_{j=1}^Ma_j^2)^{1/2}$ and $\|\ba\|_\infty=\max_{1\leq j\leq M} |a_j|$. For any real function $f$ depending on $d$ variables $(t_1, \ldots, t_d)$ and any vector $\bv=(v_1, \cdots, v_d)$ of nonnegative integers, denote $f^{(\bv)} = \frac{\partial^{\multindnorm{\bv}}}{\partial t_1^{v_1} \ldots \partial t_d^{v_d}} f$ where $\multindnorm{\bv} = \sum_{k = 1}^d v_j$. For functions that depend on $(\bx,q)$, the multi-index derivative notation is taken with respect to the first argument $\bx$, unless otherwise noted. We say a function $f$ is $\alpha$-H\"{o}lder on a set $\mathcal{I}$ if for some constant $C>0$ and $\alpha>0$, $|f(\bx_1)-f(\bx_2)|\leq C\|\bx_1-\bx_2\|^{\alpha}$ for any $\bx_1,\bx_2\in\mathcal{I}$. For any two numbers $a$ and $b$, $a\vee b=\max\{a,b\}$, and $a\wedge b=\min\{a,b\}$. Let $\E_n[g(x_i)]=\frac{1}{n}\sum_{i=1}^ng(x_i)$ and $\mathbb{G}_n[g(x_i)]=\frac{1}{\sqrt{n}}\sum_{i=1}^n(g(x_i)-\E[g(x_i)])$. For sequences, $a_n=O(b_n)$ or $a_n\lesssim b_n$ denotes $\limsup_n |a_n/b_n|$ is finite, $a_n=O_\P(b_n)$ denotes $\limsup_{\epsilon\to\infty} \limsup_{n\to\infty} \P[|a_n/b_n| \geq \epsilon ] = 0$, $a_n = o(b_n)$ denotes $a_n/b_n\to 0$, and $a_n = o_\P(b_n)$ denotes $a_n/b_n\to_\P 0$, where $\to_\P$ is convergence in probability. Limits are taken as $n \to \infty$, and the dependence on $n$ is often suppressed, e.\,g. $K = K_n$. Also, we say a random variable $\xi$ is sub-Gaussian conditional on $\bX$ if for some constant $\sigma^2>0$, $\P(|\xi|\geq t|\bX=\bx)\leq 2\exp(-t^2/\sigma^2)$ for all $t\geq 0$ and $\bx\in\mathcal{X}$.
We discuss four motivating examples of interest covered by our theoretical results. Section (ref) demonstrates how our high-level assumptions, introduced in Section (ref), are verified for these examples in order to obtain uniform estimation and inference results; the supplemental appendix collects omitted details.
Our first example generalizes the work of Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE, who studies the large sample properties of nonparametric conditional quantile series regression with $\eta(u)=u$. We allow for non-identity transformation under substantially weaker technical conditions.
chernozhukov-et-al_2013_ECMA obtains large sample estimation and inference results for parametric ($K$ fixed) generalized conditional distribution regression, and applies them to counterfactual analysis and causal inference. The following example discusses a nonparametric partitioning-based generalized conditional distribution regression estimator.
Lai-Lee_JASA_2005 studies $L_p$ regression estimation with identity transformation $\eta(\cdot)$ in a parametric setting ($K$ fixed). The next example considers a class of nonparametric partitioning-based generalized $L_p$ regression with $p\in[1,2]$, covering the full interpolation between nonparametric generalized median regression ($p=1$) and nonparametric nonlinear least squares regression ($p=2$), with possibly non-identity $\eta(\cdot)$.
The final example considers (nonparametric) Generalized Linear Models Mccullagh_2019_book. For specificity, we focus on (quasi-)maximum likelihood logistic regression, but our results cover many other examples within this class such as regression models with limited dependent variables (e.g., Poisson, fractional, censored and truncation regression).
The four examples illustrate distinct settings from a technical perspective. In Example (ref) uniformity over $\mathcal{X}\times\mathcal{Q}$ is of interest, and the loss function is non-smooth as a function of $\bx\in\mathcal{X}$ but smooth as a function of $q\in\mathcal{Q}$. Example (ref) is the antithesis of Example (ref) because now the loss function is smooth as a function of $\bx\in\mathcal{X}$ and non-smooth as a function of $q\in\mathcal{Q}$, while uniformity over $\mathcal{X}\times\mathcal{Q}$ is still of interest. In Example (ref) only uniformity over $\mathcal{X}$ is of interest because $q\in\mathcal{Q}$ is not present in the loss function, but its smoothness depends on $p\in[0,1]$; the a.\,e. derivative of $\eta\mapsto\rho(y,\eta)$ ranges from discontinuous ($p=1$), to H\"{o}lder continuous ($p\in(1,2)$), to linear ($p=2$). Likewise, Example (ref) only involves uniformity over $\mathcal{X}$ because $q\in\mathcal{Q}$ is not present in the loss function, but now the loss function is smooth and well-behaved; this last example serves as a benchmark for our theoretical development. All of the examples above have a convex loss function when $\eta(u)=u$, but can be non-convex when $\eta(\cdot)$ is not the identity function.
Our theoretical results cover other examples. For instance, Tukey and Huber regression are popular methods in robust statistics, and our theory allows for their generalizations to nonparametric partitioning-based uniform estimation and inference. Specifically, Tukey regression employs the loss function $\rho(y,\eta;q) = q^2 (1 - [1 - (y-\eta)^2/q^2]^3) \I(|y-\eta|\leq q) + q^2 \I(|y-\eta|>q)$, while Huber regression uses the loss function $\rho(y,\eta;q) = (y-\eta)^2\I(|y-\eta|\leq q)+ q(2|y-\eta|-q)\I(|y-\eta|>q)$, where $q$ is treated as a tuning parameter that balances the robustness and the bias of the estimation. We do not discuss these and other examples for brevity.
Our theoretical work proceeds under Assumptions (ref) and (ref) on the partitioning-based estimation framework, three assumptions concerning the data generating process and the loss function, and a final assumption linking the statistical model and partition-based approximation.
Assumption (ref) imposes standard conditions from the nonparametric regression literature, including basic support and smoothness restrictions. Minimal additional regularity is imposed to accommodate uniformity over $q\in\mathcal{Q}$, and different types of conditional distributions of $Y|\bX$ (e.\,g., absolutely continuous, discrete or mixed) are allowed. We consider continuously distributed covariates for simplicity, but with additional notation, and by appropriate modification of our assumptions and proof, it is possible to accommodate $\bx_i$ with continuous and discrete components.
The next assumption requires regularity conditions on the loss and transformation functions. Define $B_{q}(\bx)=\{\zeta: |\zeta-\mu_0(\bx,q)|\leq r\}$ for some fixed (small enough) constant $r>0$, which is a “ball” around the true value $\mu_0(\bx,q)$ with radius $r$.
Assumption (ref) is carefully crafted to accommodate all the examples discussed in Section (ref), and many others. Part (i) allows for different degrees of smoothness in the loss function, assuming only absolute continuity (with respect to $\eta$). It also makes clear that $q$ is scalar, which is assumed only to simplify the notation; see the supplemental appendix for the general case $\mathcal{Q}\subseteq\mathbb{R}^{d_\mathcal{Q}}$ with $d_\mathcal{Q}\geq1$. Part (ii) formalizes the idea that $\mu_0(\bx,q)$ may not be a unique (global) minimizer in (ref), and consequently it is only required to be a root of the (conditional) first-order condition; the rest of the assumptions in that part are mild regularity conditions. In some applications, $\mu_0(\bx,q)$ can be the unique minimizer; see, for example, Lin-Kulasekera_2007_Biometrika, matzkin2007nonparametric, and references therein.
Part (iii) of Assumption (ref) imposes additional structure on the a.\,e. first derivative of the loss function, allowing for all types of outcome data (discrete, mixed, and continuous) and rescalings emerging in some of the motivating examples. Importantly, this part characterizes precisely the role of (H\"older) smoothness, which is controlled by the parameter $\alpha\in(0,1]$. We illustrate the full power of this general assumption in Section (ref), where $\alpha=1$ in Examples (ref) and (ref), $\alpha=p-1$ in Example (ref) when $p>1$, and $\alpha=1$ in Example (ref). The strict monotonicity condition on $\eta(\cdot)$ is satisfied by usual (inverse) link functions used in generalized linear models. Finally, part (iv) of Assumption (ref) collects mild regularity conditions on the smoothed-out a.\,e. derivative of the loss function.
Assumptions (ref) and (ref) have restricted basic aspects of the statistical model, imposing standard support, moment, and smoothness conditions, in addition to other minimal structure required on the loss and transformation functions. These conditions are sufficient for pointwise estimation and inference, but more is needed for uniform over $\mathcal{X}\times\mathcal{Q}$ results. In the supplemental appendix, our theoretical results are established under one more condition that governs the complexity of the loss function and related function classes. To avoid a long list of complexity bounds, we present a more restrictive but simpler assumption motivated by the examples discussed in Section (ref): we consider a loss function $\rho(y, \eta; q)$ that can be expressed as a linear combination of certain “simple” functions. See Section \saref{sec:simpler-generalization-of-examples} for omitted details.
This third assumption imposes an additional weak monotonicity condition on $\mu_0(\bx, q)$ as a function of $q$, which is compatible with all the examples in Section (ref). Finally, the key restriction emerging from Assumption (ref) is on the structure of the loss function, which allows for linear combinations of smooth loss functions of $y$ and $\eta$, and non-smooth loss functions involving indicator functions of either $y$ and $\eta$, or $y$ and $q$. These restrictions are still general enough to cover the four motivating examples: the loss function in Example (ref) is a combination of Type I and Type III functions with $f_1$ and $f_3$ being linear functions of $y$; the loss function in Example (ref) is a combination of Type II and Type IV functions; the loss function in Example (ref) is of Type IV for $p>1$, and of the same type as Example (ref) when $p=1$ (median regression); and the loss function in Example (ref) is usually a Type IV function. See Section (ref) for details.
Our final assumption concerns the approximation power of the basis $\bp(\cdot)$ in connection with the underlying functional parameter.
The vector $\bbeta_0(q)$ can be viewed as a pseudo-true value, and does not have to be unique. The existence of such $\bbeta_0(q)$ can be established using approximation theory or related methods, and necessarily depends on the specific underlying structure of the statistical model (determining $\mu_0(\bx,q)$) and the partitioning-based method (determining $\bp^{(\boldsymbol{\varsigma})}(\bx)$). For standard local bases, the assumption can be verified by imposing smoothness conditions on $(\bx,q)\mapsto\mu_0^{(\boldsymbol{\varsigma})}(\bx,q)$ (see Assumption (ref)(iv)). For more discussion, see Belloni-Chernozhukov-Chetverikov-Kato_2015_JoE, Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE, Cattaneo-Farrell_2013_JoE, Cattaneo-Farrell-Feng_2020_AOS, Chen-Christensen_2015_JOE, Huang_2003_AoS, and references therein.
We show that the partitioning-based $M$-estimator is consistent, which is the starting point for establishing its main point estimation and inference asymptotic properties. We endeavor to impose the weakest possible conditions, which requires careful consideration of the specific shape of the loss function in (ref): we consider two cases, either the loss function $\rho(y,\eta(\theta);q)$ is convex with respect to $\theta$ or not; in the latter case, we will have to restrict the feasibility region $\mathcal{B}$.
For the case of convex $\theta\mapsto\rho(y,\eta(\theta);q)$, consistency can be established for general unconstrained estimators ($\mathcal{B}=\mathbb{R}^K$ in (ref)) under mild conditions. The proof is deferred to the supplemental appendix (Lemma \saref{lem:holder-consistency}).
Lemma (ref) shows that the function estimator $\widehat\mu$ is uniform-in-$q$ consistent for the true value $\mu_0$ in both $L_2$-norm and $\sup$-norm over $\mathcal{X}$, whereas for the derivative estimator $\widehat\mu^{(\bv)}$ (with $|\bv|>0$) the lemma only provides a bound on its deviation from the estimand $\mu_0^{(\bv)}$. Technically, all we need from this lemma to establish the Bahadur representation later is the uniform-in-$q$ consistency of the coefficients estimator $\widehat\bbeta(q)$ for the pseudo-true coefficients $\bbeta_0(q)$ in $\sup$-norm, i.e., $\|\widehat\bbeta(q)-\bbeta_0(q)\|_\infty=o_\P(1)$, which is immediate from (ref), the uniform-in-$q$ consistency in the Euclidean norm.
Two kinds of rate restrictions are imposed in Lemma (ref), depending on the moment condition assumed for the generalized residual $\psi(y_i,\eta(\mu_0(\bx_i,q));q)$. In the best case when the residual has a sub-Gaussian envelope, we need $1/(nh^{2d}) \asymp K^2/n = o(1)$, up to $\polylog(n)$ terms, while in the worst case when the envelope of $\psi(y_i,\eta(\mu_0(\bx_i,q));q)$ has a bounded $\nu$-th moment with $\nu$ close to $2$, we roughly need $1/(nh^{4d}) \asymp K^4/n = o(1)$, up to $\polylog(n)$ terms.
A feature of Lemma (ref) is that no constraints are imposed on the coefficients in the optimization procedure, which allows the estimation space to be, for example, piecewise polynomials. In contrast, many studies of series (or sieve) methods restrict the functions in the estimation space to satisfy certain smoothness conditions, e.g., Lipschitz continuity, to derive the uniform consistency chernozhukov-imbens-newey_2007_JoE.
Consider the case when the loss $\rho(y,\eta(\theta); q)$ is possibly non-convex with respect to $\theta$. This setting naturally arises, for example, in nonlinear regression when $\rho(y,\eta(\theta);q) = (y - \eta(\theta))^2$ with $\eta(\cdot)$ non-identity: while $\eta\mapsto\rho(y,\eta;q)$ is a square loss function, hence convex, introducing a transformation function $\eta$ such as the (inverse) logistic link will often make $\theta\mapsto\rho(y,\eta(\theta);q)$ non-convex.
A proof of consistency for the unconstrained estimator in (ref) with a non-convex loss function is not available, but we are able to establish consistency of a regularized $M$-estimator. Specifically, we add a fixed “box” constraint: for some fixed constant $R>0$,
In the supplemental appendix we show that the pseudo-true coefficients $\bbeta_0(q)$ from Assumption (ref) are bounded in $\sup$-norm by a universal constant: $\sup_{q\in\mathcal{Q}}\|\bbeta_0(q)\|_\infty \lesssim 1$ (because $\bp(\bx)\trans \bbeta_0(q)$ has to be close to $\mu_0(\bx, q)$ which is uniformly bounded). Therefore, we can always choose a sufficiently large constant $R$ in the optimization procedure, making the box constraint set contain $\bbeta_0(q)$ as an interior point. The following lemma, proven in the supplemental appendix (Lemma \saref{lem:consistency-nonconvex}), establishes consistency of the constrained estimator.
Compared to Lemma (ref), two additional restrictions are imposed in this lemma. The first one, $R\geq 2\sup_{q\in\mathcal{Q}}\|\bbeta_0(q)\|_\infty$, can be theoretically justified by Lemma \saref{lem:assumptions-of-nonconvex-consistency-are-reasonable} in the supplemental appendix, and in practice a large enough $R$ is recommended. The other restriction concerns a lower bound for $\Psi_1$, which implies that the (population) loss function is strongly convex in a neighborhood of the true value $\eta(\mu_0(\bx,q))$, making the (constrained) minimizer well defined. (This condition does not\/ break because of the shape of $\eta(\cdot)$, in contrast with the convexity of $\rho(y, \eta(\theta); q)$ in $\theta$.) The other conditions in this lemma are the same as those in the convex case, and thus all improvements discussed before also apply to this second consistency result.
In the supplemental appendix we provide additional consistency results for two special cases:
The first case covers the usual square loss function with identity transformation. The second case covers partitioning-based $M$-estimation using the (Haar and) piecewise polynomial basis. Notably, in these two special cases, the consistency result $\|\widehat\bbeta(q)-\bbeta_0(q)\|_\infty=o_\P(1)$ is established for any $m$ and $d$, so the requirement $m>d/2$ imposed in Lemmas (ref) and (ref) is not needed. Furthermore, in these two cases, it only requires the minimal side rate restrictions $1/(nh^d) \asymp K/n = o(1)$ in the sub-Gaussian case, and $1/(nh^{\frac{\nu}{\nu-1}d}) \asymp K^{\frac{\nu}{\nu-1}}/n = o(1)$ in the bounded $\nu$-th moment case, up to $\polylog(n)$ terms. See Section \saref{sec:consistency} in the supplemental appendix for more details.
These restrictions imposed in Lemmas (ref) and (ref), and the associated results in the supplemental appendix for special cases, are either comparable to or improve upon the existing literature. In the special case of square loss function and identity transformation, uniform (over $\bx\in\mathcal{X}$) consistency of the partitioning-based estimator is essentially automatic due to the intrinsic closed-form and linearity of the estimator. Nevertheless, compared to the best known result in that strand of the literature Cattaneo-Farrell_2013_JoE,Belloni-Chernozhukov-Chetverikov-Kato_2015_JoE,Cattaneo-Farrell-Feng_2020_AOS, our general results are essentially on par in terms of side rate restrictions and regularity conditions. For example, in the sub-Gaussian case, and up to $\polylog(n)$ terms, the best side rate restriction in that literature requires $1/(nh^{d}) \asymp K/n = o(1)$, while Lemmas (ref) and (ref) require $1/(nh^{2d}) \asymp K^2/n = o(1)$, and our improved results in the supplemental appendix for unconnected bases require $1/(nh^{d}) \asymp K/n = o(1)$. Therefore, our results are on par with the best available results for series-based least squares regression Cattaneo-Farrell_2013_JoE,Belloni-Chernozhukov-Chetverikov-Kato_2015_JoE,Cattaneo-Farrell-Feng_2020_AOS, despite being able to cover a large class of $M$-estimation settings such as piecewise-polynomial-based quantile, nonlinear, or robust regression.
In the case of quantile regression with tensor-product $B$-splines, Corollary 1 of Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE implies that the $L_2$-consistency (ref) can be obtained under $1/(nh^{2d}) \asymp K^2/n = o(1)$, and their Corollary 2 implies that the uniform consistency (ref) can be obtained under $1/(nh^{4d}) \asymp K^4/n = o(1)$. In contrast, since the generalized residual from quantile regression has a sub-Gaussian envelope, we only require $1/(nh^{2d}) \asymp K^2/n = o(1)$ to establish both kinds of consistency. Moreover, when an unconnected basis is used, or when the loss function is strongly convex and smooth (e.g., the square loss for mean regression), we establish consistency under the weakest possible restriction: $1/(nh^d) \asymp K/n = o(1)$, up to $\polylog(n)$ terms. It remains an open question whether it is possible to establish consistency under the weakest condition $1/(nh^d) \asymp K/n = o(1)$ for general partitioning-based $M$-estimators.
The Bahadur representation is
where
The following theorem takes the sup-norm consistency of the coefficient estimators $\widehat\bbeta(q)$ as a high-level assumption, and thus avoids imposing any of the specific sufficient conditions discussed in Section (ref). The proof is provided in the supplemental appendix (Theorem \saref{th:bahadur-repres}).
The Bahadur representation (ref) applies to the case where the “derivative” $\psi(\cdot,\cdot;q)$ of the loss function may be discontinuous. One typical example is quantile regression (Example (ref)), where the “derivative” $\psi(y,\eta;q)=\I(y-\eta<0)-q$, as a function of $(y-\eta)$, is piecewise constant with a jump at zero. In this case we can let $\alpha=1$, and (ref) implies that the order of the remainder in the Bahadur representation for partitioning-based quantile regression is $O(h^{-|\bv|}(nh^{d})^{-3/4}+h^{m-|\bv|})$, up to $\polylog(n)$ terms. Another example is $L_p$ regression with $p\in(1, 2)$, Example (ref), where the derivative of the loss function is $\psi(y,\eta)\equiv \psi(y-\eta)=p|y-\eta|^{p-1}\sgn(\eta-y)$ with $\sgn(\cdot)$ denoting the sign function. As a function of $(y-\eta)$, $\psi(\cdot)$ is $\alpha$-H\"{o}lder on $[0,\infty)$ or $(-\infty,0]$ for all $\alpha\in(0, p-1]$ but not for $\alpha>p-1$. Thus, (ref) applies with the order of the remainder depending on $p$, which is the same as that for quantile regression when $p\geq 3/2$.
On the other hand, the Bahadur representation (ref) applies to the case where the “derivative” of the loss is a continuous function of $(y-\eta)$. Nonlinear least squares regression (Example (ref) and quasi-maximum likelihood estimation of generalized linear models (Example (ref)) fall into this category with the H\"{o}lder parameter $\alpha=1$. In such cases, (ref) implies that the order of the remainder in the Bahadur representation is $O(h^{-|\bv|}(nh^{d})^{-1}+h^{m-|\bv|})$, up to $\polylog(n)$ terms, which is a tighter upper bound than that implied by (ref). See Section (ref) and the supplemental appendix for more details.
In both cases, the remainder of the Bahadur representation consists of two terms. The last term $h^{m-|\bv|}$ corresponds to the error from approximating the function $\mu_0$ using the partitioning basis (cf. Assumption (ref)), whereas the first term arises from the (potential) nonlinearity underlying the $M$-estimation, and reflects explicitly the role of non-smoothness of the loss function. When the “derivative” of the loss function has discontinuity points, the order of the remainder in (ref) is greater than that in the continuous case (ref); with a smaller H\"{o}lder parameter $\alpha$, the order of the remainder in both cases could increase.
The uniform Bahadur representations (Theorem (ref)) can be used to establish convergence rates for the general partitioning-based $M$-estimators. We focus first on uniform convergence over $\bx\in\mathcal{X}$ and $q\in\mathcal{Q}$.
By setting $h\asymp \big(\frac{\log n}{n}\big)^{\frac{1}{2m+d}}$, Corollary (ref) implies that the partitioning-based $M$-estimator can achieve the uniform convergence rate $\big(\frac{\log n}{n}\big)^{\frac{m}{2m+d}}$. This matches the optimal rate of convergence in $\sup$-norm for nonparametric estimators of the conditional mean Stone_1982_AoS and conditional quantiles chaudhuri1991global. In this sense, the rate of convergence in Corollary (ref) is optimal and cannot be further improved at our level of generality.
Theorem (ref) can also be used to obtain the mean square convergence rate of the partitioning-based $M$-estimator uniformly-in-$q$.
By setting $h\asymp n^{-\frac{1}{2m+d}}$, Corollary (ref) implies that the partitioning-based $M$-estimator can also achieve the $L_2$ convergence rate $n^{-\frac{m}{2m+d}}$, uniformly over $\mathcal{Q}$, thereby matching the optimal rate of convergence in $L_2$-norm for nonparametric estimators of conditional means Stone_1980_AoS and conditional quantiles chaudhuri1991global.
The convergence rates in (ref) and (ref) capture two contributions: the first term reflects the variance of the estimator, while the second term arises from the error of approximating the unknown $\mu_0$ by the partitioning basis. In the case of Corollary (ref), it is possible to further leverage Theorem (ref) to obtain a precise first-order asymptotic approximation for the integrated mean square error of the partitioning-based $M$-estimator, uniformly over $\mathcal{Q}$, which in turn could be used to develop plug-in asymptotically optimal rules for selecting $K\asymp h^{-d}$. See, for example, Theorem 4.2 in Cattaneo-Farrell-Feng_2020_AOS for a similar result in the special case of square loss and identity transformation functions. We do not pursue this result here for brevity.
To our knowledge, this paper is the first to establish uniformly valid Bahadur representations for partitioning-based $M$-estimators at the level of generality allowed in Theorem (ref), and the implied convergence rates in Corollaries (ref) and (ref). The restriction on the tuning parameter $h$ required by the theorem is seemingly minimal: when the envelope of the generalized residual $\psi(y_i,\eta(\mu_0(\bx_i,q));q)$ is sub-Gaussian (or its $\nu$-th moment is bounded with a large $\nu$), we roughly only need $1/(nh^d) \asymp K/n = o(1)$, up to $\polylog(n)$ terms. Having noted this, verification of the high-level consistency assumption $\|\widehat\bbeta(q)-\bbeta_0(q)\|_\infty=o_\P(1)$ in the sub-Gaussian case may require a more stringent condition on $h$, as discussed in Section (ref). In the best scenario (e.g., an unconnected basis is used), the minimal restriction $1/(nh^d) \asymp K/n = o(1)$ suffices, while in the worst scenario we need at most $1/(nh^{2d}) \asymp K^2/n = o(1)$, up to $\polylog(n)$ terms.
The rest of this section discusses precisely how our results improve on prior literature.
The usual mean regression is a special case of our general setup where $\rho(\cdot,\cdot)$ is the square loss, $\eta(\cdot)$ is the identity link, and $\mathcal{Q}$ is a singleton. Bahadur representations for this special case were established by Belloni-Chernozhukov-Chetverikov-Kato_2015_JoE and Cattaneo-Farrell-Feng_2020_AOS. Since the derivative of the square loss for mean regression is linear, the first term in (ref) or (ref) does not show up in the uniform linearization of least squares series estimators. See, for example, Lemma SA-4.2 of Cattaneo-Farrell-Feng_2020_AOS; $R_{1n,q}$ defined therein has been implicitly included in the leading variance term in (ref) above. Theorem (ref) substantially extends these prior results to other nonlinear settings, under minimal additional conditions.
Finally, Corollaries (ref) and (ref) demonstrate the convergence rate optimality of general partitioning-based series $M$-estimation, recovering in particular known results for mean regression Belloni-Chernozhukov-Chetverikov-Kato_2015_JoE,Cattaneo-Farrell-Feng_2020_AOS under essentially the same minimal conditions.
Theorem (ref) improves upon prior theoretical results for nonparametric series quantile regression estimators. The most recent advance in this literature is due to Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE, which establishes a uniform linear approximation for general series-based quantile regression estimators. In comparison, we exploit the “local support” feature of the partitioning basis, and make improvements in (at least) four aspects. To summarize these improvements without additional cumbersome notation, we set $\bv=\bm{0}$ and ignore the smoothing bias $h^m$ in the Bahadur approximation remainders.
First, Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE shows that the order of the remainder in the Bahadur representation is $O((nh^d)^{-3/4} h^{-d/2})$, up to $\polylog(n)$ terms (see proofs of Theorem 2 and Corollary 2 therein for details). In contrast, Theorem (ref) implies that the remainder in the Bahadur representation for partitioning-based quantile regression estimators is $O((nh^d)^{-3/4})$, up to $\polylog(n)$ terms, which is not only a much tighter bound but also matches the optimal parametric bound when taking $nh^d$ as the effective sample size.\sloppy
Second, the rate restriction $1/(nh^{4d}) \asymp K^4/n = o(1)$ is required for $B$-spline-based estimators in Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE. In contrast, the restriction on $h$ in Theorem (ref) depends on the tail behavior of the generalized residuals and becomes weaker as $\nu$ gets larger. In the best case (the residuals have a sub-Gaussian envelope) we only need the seemingly weakest restriction $1/(nh^d) \asymp K/n = o(1)$, up to $\polylog(n)$ terms, along with the consistency condition for $\widehat\bbeta(q)$. Recall that in the sub-Gaussian scenario we need at worst $1/(nh^{2d}) \asymp K^2/n =o(1)$, up to $\polylog(n)$ terms, to satisfy the consistency requirement.
Third, the restriction $h^{m-d}=o(n^{-\varepsilon})$ for some $\varepsilon>0$ in Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE implicitly requires the smoothness $m$ of the conditional quantile function be greater than the dimensionality $d$ of the covariates. In contrast, the proof of Theorem (ref) does not need such a restriction, though a weaker condition $m>d/2$ might be needed to verify the consistency condition on $\widehat{\bbeta}(q)$; see Lemmas (ref) and (ref). Furthermore, when an unconnected basis (e.g., piecewise polynomials) is used for approximation, the condition $m>d/2$ is unnecessary for consistency, and thus we have no constraint on the relation between smoothness $m$ and dimensionality $d$; see Section (ref).
Fourth, compared to Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE, we allow for a possibly non-identity link. Introducing a link function may lead to non-convexity of the loss $\rho(y,\eta(\theta);q)$ with respect to $\theta$, making the usual proof strategies for consistency and Bahadur representation under convexity inapplicable. For example, non-convex quantile regression is covered in Theorem (ref) by virtue of our general consistency results in Lemma (ref).
All of the aforementioned improvements are practically relevant. For example, they accommodate univariate quantile regression using the piecewise constant basis with the IMSE-optimal choice of the mesh size $h$ (in this case $h\asymp n^{-1/3}$ and $m=d$), which was theoretical excluded in prior literature.
Finally, Corollaries (ref) and (ref) establish the optimal rate of convergence for general partitioning-based series $M$-estimators, which substantially improve on prior work on quantile series regression in particular. More specifically, the conditions on the mesh size $h$, the smoothness $m$, and the dimensionality $d$ in both corollaries are weaker than in prior work. In the best case (e.g., an unconnected basis is used and a sub-Gaussian envelope for residuals exists), we only require the seemingly minimal restriction $1/(nh^d) \asymp K/n =o(1)$, up to $\polylog(n)$ terms, and an arbitrary relation between $m$ and $d$ is permitted. For comparison, in the special case of quantile regression, Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE shows that series estimators can achieve the fastest possible uniform-in-$q$ rate in both $L_2$-norm and $\sup$-norm (see Comments 3 and 4 therein), but under more stringent conditions: $\eta$ is the identity function, $m>d$, and $1/(nh^{4d}) \asymp K^4/n = o(n^{-\varepsilon})$ for some $\varepsilon>0$ (see their Corollary 2). Such conditions exclude, e.g., the IMSE-optimal choice of $h$ for the Haar basis when $d=1$ (since $h\asymp n^{-1/3}$ and $m=d=1$), or piecewise linear fit when $d=2$ (since $m=d=2$).
Kong-Linton-Xia_2013_ET establishes a similar Bahadur representation for kernel-based $M$-estimators using weakly stationary time series data. They consider a special case of our setup in Assumption (ref): their loss function class $\mathcal{Q}$ is a singleton, $\eta$ is an identity function, and the “derivative” of the loss can be written as $\psi(y,\eta)\equiv\psi(y-\eta)$ and is assumed to be piecewise Lipschitz continuous. In a comparable cross-sectional context with $\alpha=1$ and $\bv=\bm{0}$, the order of the remainder in our Bahadur representation (ref) is $O((nh^d)^{-3/4})$, up to $\polylog(n)$ terms, and thus Theorem (ref) matches their approximation error up to a minor difference in $\log n$ terms. Taking $nh^d$ to be the effective sample size, the approximation rate can not be further improved at this level of generality, and hence Theorem (ref) establishes that the partitioning-based series $M$-estimator in (ref) can achieve the same best Bahadur approximation as local polynomial kernel methods, up to $\polylog(n)$ terms.
Furthermore, compared to Kong-Linton-Xia_2013_ET or other similar contributions in the literature, Theorem (ref) exhibits (at least) two novel features. First, the Bahadur representations (ref) and (ref) hold uniformly not only over the evaluation point $x\in\mathcal{X}$, but also over the loss function index $q\in\mathcal{Q}$, which may be important, for example, to study simultaneous quantile regression where the entire conditional quantile process may be of interest. Second, Theorem (ref) also covers the more general setup where the “derivative” function may exhibit different degrees of smoothness, reflected by discontinuity points and/or the H\"{o}lder parameter $\alpha$, or admits a more complex structure so that $\psi(y,\eta;q)$ cannot be written as $\psi(y-\eta;q)$. Thus, we cover more examples such as distribution regression (Example (ref)) and $L_p$ regression with $p\in(1,2)$ (Example (ref)). Finally, Kong-Linton-Xia_2013_ET does not discuss convergence rates as we do in Corollaries (ref) and (ref).
In the context of nonparametric penalized smoothing spline methods, Shang-Cheng_2013_AOS also establishes a uniform Bahadur representation (and other results) that can be compared to Theorem (ref). However, their paper imposes more stringent assumptions and hence cover a smaller class of settings: using our notation, they assume that (i) $\mathcal{Q}$ is a singleton so their uniformity is only over $\mathcal{X}$; (ii) $d=1$ so they consider only scalar covariate $\bx_i$; and (iii) $\rho(\cdot,\eta(\cdot))$ is smooth so they rule out many important examples such as quantile regression, and Tukey and Huber regression. Furthermore, their results do not take explicit advantage of specific moment and boundedness conditions, or the structure of the nonparametric estimator, and instead impose the generic side condition $nh^2\to\infty$, which is comparable to our condition $K^2/n\to\infty$, up to $\polylog(n)$ terms. Most importantly, in the closest comparable case ($d=1$, $\alpha=1$, and $\bv=\bm{0}$), and only focusing on the variance component for simplicity, the order of the remainder in their uniform Bahadur representation (a combination of Theorem 3.4 and Lemma 3.1 in Shang-Cheng_2013_AOS) is $O((nh)^{-1} h^{-(6m-1)/(4m)})$, while (ref) in Theorem (ref) gives the optimal result $O((nh)^{-1})$, thereby demonstrating a substantial improvement over their result. Finally, as for convergence rates, Proposition 3.3 in Shang-Cheng_2013_AOS and our Corollary (ref) are essentially equivalent, both delivering optimal mean square convergence. They do not explicitly discuss uniform convergence rates as we do in Corollary (ref).
The uniform Bahadur representations in Theorem (ref) can also be leveraged to establish uniform distribution theory for $\widehat\mu^{(\bv)}$. The infeasible conditional variance of the estimator can be written as
where
Accordingly, we define a feasible variance estimator as
where $\widehat\bQ_q$ and $\widehat\bSigma_q$ are some estimators of $\bar\bQ_q$ and $\bar\bSigma_q$, respectively, which are consistent in a sense described below. Therefore, $\widehat\Omega_{\bv}(\bx,q)$ is an estimator of the infeasible conditional variance $\bar\Omega_{\bv}(\bx,q)$.
Statistical inference on $\mu_0^{(\bv)}$ usually relies on the following $t$-statistic process:
where we drop the dependence of $T(\cdot)$ on $\bv$ for simplicity.
Employing Theorem (ref), or more precise arguments under slightly weaker conditions, it is easy to show that $T(\bx,q)$ converges in distribution to $\mathsf{N}(0,1)$ for each $(\bx,q)\in\mathcal{X}\times\mathcal{Q}$. However, the stochastic process $(T(\bx,q):(\bx,q)\in\mathcal{X}\times\mathcal{Q})$ is generally not asymptotically tight and, therefore, does not converge weakly in $\ell^\infty(\mathcal{X}\times\mathcal{Q})$, where $\ell^\infty(\mathcal{X}\times\mathcal{Q})$ denotes the set of all (uniformly) bounded real functions on $\mathcal{X} \times \mathcal{Q}$ equipped with uniform norm vaart96_weak_conver_empir_proces. Nevertheless, we can construct a Gaussian process, in a possibly enlarged probability space, that approximates the entire process $T(\cdot)$ sufficiently fast, which can then be used to approximate the finite sample distribution of the function $M$-estimator $\widehat{\mu}^{(\bv)}(\cdot)$.
More precisely, under some mild consistency conditions on $\widehat\Omega_{\bv}(\bx,q)$, our Theorem (ref) guarantees that $\sup_{q\in\mathcal{Q}} \sup_{\bx\in\mathcal{X}} \big| T(\bx,q) - t(\bx,q) \big| \to_\P 0$ sufficiently fast, where
It follows that, conditional on $\{\bx_i\}_{i=1}^n$, the randomness of $t(\bx,q)$ comes exclusively from the $K$-dimensional vector $\mathbb{G}_n[\bp(\bx_i)\eta^{(1)}(\mu_0(\bx_i,q)) \psi(y_i,\eta(\mu_0(\bx_i,q));q)]$. Thus, our proof strategy is to further “discretize” this vector with respect to $q\in\mathcal{Q}$, and then apply Yurinskii's coupling yurinskii1978error to construct a (conditional) Gaussian process that is close to the original $t$-statistic process $T(\bx,q)$ uniformly over both $\bx\in\mathcal{X}$ and $q\in\mathcal{Q}$. Our construction leverages a conditional Strassen’s theorem chen2020jackknife to generalize prior coupling results Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE. See Section \saref{sec:strong-approximation} in the supplemental appendix for details.
Our strong approximation approach is formalized in the next theorem. We employ high-level conditions to ease the exposition, but those conditions can be verified using Corollaries (ref) and (ref), and Theorem (ref), as well as using the more general results in the supplemental appendix. Let $r_{\tt UC}$, $r_{\tt BR}$, $r_{\tt VC}$, and $r_{\tt SA}$ be positive non-random sequences as $n\to\infty$. The proof is available in the supplemental appendix (Theorem \saref{th:strong-approximation-binscatter-yurinski}).
The speed of strong approximation in Theorem (ref) is determined by four factors: the uniform convergence rate $r_{\tt UC}$, the order of the remainder in the Bahadur representation $r_{\tt BR}$, the convergence rate $r_{\tt VC}$ of the variance estimator $\widehat\Omega_{\bv}$, and the strong approximation rate $r_{\tt SA}$. Therefore, our strong approximation results are established at a high level of generality, building on our prior theoretical results: Corollary (ref) for $r_{\tt UC}$, and Theorem (ref) for $r_{\tt BR}$, while $r_{\tt VC}$ is a high-level condition that needs to be verified on a case-by-case basis. See Sections (ref) and (ref) for more discussion.
With respect to the strong approximation rate, Theorem (ref) lays down two versions of lower bounds on $r_{\tt SA}$, depending on the tail behavior of the generalized residuals. Such restrictions may not be optimal, but are still weak enough to cover almost all partition size choices commonly used in practice. In particular, the restriction on $r_{\tt SA}$ in Theorem (ref) allows for the MSE-optimal choice $h\asymp n^{-\frac{1}{2m+d}}$ in all cases except the unidimensional Haar basis approximation ($m=d=1$); there is also room for undersmoothing in order to make the smoothing bias negligible in all cases but $m=d=1$. The strong approximation for one dimensional partitioning-based series estimators in the special case of square loss and identity transformation functions was studied in Cattaneo-Farrell-Feng_2020_AOS,Cattaneo-Crump-Farrell-Feng_2024_AER via a different coupling strategy, which delivered tighter approximation results allowing for an MSE-optimal choice of $h$. We conjecture those techniques could be adapted to cover the case $m=d=1$ for general partitioning-based $M$-estimator in (ref), but we do not pursue this line of research here because it would require a different theoretical treatment.
Theorem (ref) is the first to establish strong approximation results for general partitioning-based $M$-estimators at the level of generality considered in this paper. In the prior literature, similar results are usually available only in specific scenarios such as least squares regression or quantile regression. To be more precise, in the least squares context ($\mathcal{Q}$ is a singleton), Cattaneo-Farrell-Feng_2020_AOS establishes uniform inference theory for univariate regression ($d=1$) and multivariate regression ($d>1$) separately via different strong approximation methods. In particular, when $d>1$, the same Yurinskii coupling technique is employed to obtain strong approximation for $t$-statistic processes, leading to similar rate restrictions on $h$. Theorem (ref) is a substantial generalization of results therein, not just covering other loss functions, but also providing distributional approximation uniformly over the loss function index $q\in\mathcal{Q}$.
In the quantile regression context, Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE provides two strong approximations for general series-based estimators. When $B$-splines are used, their first strategy relies on a pivotal coupling, imposing $1/(nh^{10 d}) = o(n^{-\varepsilon})$ and $h^{m-d}=o(n^{-\varepsilon})$ for some constant $\varepsilon>0$ (see Theorem 11 therein), while the second strategy uses a Gaussian coupling (as in this paper), imposing $1/(nh^{4d\vee(2+3d)}) = o(n^{-\varepsilon})$ and $h^{m-d}=o(n^{-\varepsilon})$ (see Theorem 12 and Comment 13 therein). In comparison, our Theorem (ref) requires weaker conditions on the tuning parameter $h$ and the relation between the smoothness $m$ and the dimensionality $d$. Specifically, we assume $1/(nh^{3d})=o(1)$ up to $\polylog(n)$ terms for a valid approximation. This improvement is practically relevant: for example, it allows for Gaussian approximation of linear-spline-based univariate quantile regression estimators with the MSE-optimal mesh size $h\asymp n^{-1/5}$. In addition, our general strategy to verify the consistency condition on $\widehat\bbeta(q)$ in Theorem (ref) only requires $m>d/2$ (not required for the special case of unconnected basis), which is weaker than $m>d$ as implicitly assumed in Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE. In practice, this improvement can accommodate, for example, the use of cubic splines for trivariate quantile regression.
Finally, Shang-Cheng_2013_AOS establishes uniform inference results for nonparametric penalized smoothing spline M-estimators. As mentioned before, their work is more specialized because they assume that $d=1$, $\mathcal{Q}$ is a singleton, and $\rho(\cdot,\eta(\cdot))$ is smooth. Furthermore, their approach to constructing valid confidence bands and related uniform inference methods relies on approximating the suprema of the stochastic process directly via extreme value theory Shang-Cheng_2013_AOS, which leads to substantially slower approximation rates and requires stronger assumption and side rate restrictions; see Belloni-Chernozhukov-Chetverikov-Kato_2015_JoE and Cattaneo-Farrell-Feng_2020_AOS for more discussion in the context of nonparametric least squares series estimation. In contrast, Theorem (ref), and our related uniform inference methods, provide a pre-asymptotic approximation with better finite sample properties, faster approximation rates, and weaker regularity conditions.
The (conditional) Gaussian process $(Z(\bx,q):(\bx,q)\in\mathcal{X}\times\mathcal{Q})$ in Theorem (ref) is still infeasible since its covariance structure contains unknowns. This section establishes the validity of a generic plug-in method to construct a feasible version of $Z(\bx,q)$. We later employ this result in Section (ref) to develop feasible uniform inference in the context of our motivating examples.
The core idea behind the plug-in method is to estimate the covariance structure of $Z(\bx,q)$ and then simulate its feasible version $\widehat{Z}(\bx,q)$, a Gaussian process conditional on the data. If the covariance estimate converges to the true covariance sufficiently fast, $\widehat{Z}(\bx,q)$ will be “close” to a copy of $Z(\bx,q)$. The covariance structure of the process $Z(\bx,q)$ in Theorem (ref) is
for all $(\bx,q),(\tilde{\bx},\tilde{q})\in\mathcal{X}\times\mathcal{Q}$, where
with
Given context-specific estimates $\widehat{\bQ}_q$ and $\bx \mapsto \widehat{S}_{q, \tilde{q}}(\bx)$, we can put
and $\widehat{\Omega}_{\bv}(\bx,q) = \bp^{(\bv)}(\bx)\trans\widehat\bQ_q^{-1}\widehat\bSigma_{q, q} \widehat\bQ_q^{-1}\bp^{(\bv)}(\bx)$ as above. Section (ref) illustrates how the estimates $\widehat{\bQ}_q$ and $\bx \mapsto \widehat{S}_{q, \tilde{q}}(\bx)$ can be constructed in specific examples. Then, a feasible Gaussian approximation $\widehat{Z}(\bx,q)$ can be constructed as a mean-zero Gaussian process conditional on the data $\bD_n = ((y_1,\bx_1),\cdots,(y_n,\bx_n))$ with conditional covariance structure
for all $(\bx,q),(\tilde{\bx},\tilde{q})\in\mathcal{X}\times\mathcal{Q}$.
The following theorem establishes the validity of the plug-in approach.
Once we have a feasible process $\widehat{Z}(\bx,q)$ that is “close” to a copy of $Z(\bx,q)$ uniformly over $\mathcal{X}\times\mathcal{Q}$, conditional on the data, then $\widehat{Z}(\bx,q)$ can be used to conduct inference on the entire function $\mu_0(\bx,q)$, and functionals thereof. For example, our strong approximation results can be converted to convergence of the Kolmogorov distance between the distributions of $\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}|T(\bx,q)|$ and its feasible Gaussian approximation $\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}|\widehat{Z}(\bx,q)|$. See Theorem \saref{th:k-s-distance} in the supplemental appendix for the formal result.
Furthermore, Theorem \saref{th:confidence-bands} in the supplemental appendix establishes the asymptotic validity of the uniform confidence band for $\mu_0^{(\bv)}$ given by
with $\mathfrak{c}_{1-\alpha}$ satisfying \[\P\Big(\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}|\widehat{Z}(\bx,q)|\leq \mathfrak{c}_{1-\alpha} \Big| \bD_n \Big) =1-\alpha+o_\P(1), \] provided the smoothing (or misspecification) bias relative to the standard error of the estimator is small, which could be achieved by undersmoothing, bias correction hall1992effect, simply ignoring the bias hall2001bootstrapping, robust bias correction Calonico-Cattaneo-Farrell_2018_JASA, Calonico-Cattaneo-Farrell_2022_Bernoulli, or the Lepskii's method lepskii1992asymptotically,birge2001alternative, among other possibilities. Thus, under regularity conditions, the confidence band (ref) covers $\mu_0^{(\bv)}$ with probability approximately $1-\alpha$ in large samples, that is,
We verify assumptions in our four motivating examples (Section (ref)). Assumptions (ref) and (ref) concern the partitioning-based methodology itself, Assumption (ref) are already primitive conditions on the data generating process, and Assumption (ref) can be verified for usual local bases (e.g., piecewise polynomials and splines) when we assume the functional parameter $\mu_0$, such as the conditional quantile function in Example 1 or conditional distribution function in Example 2, is smooth enough (see Assumption (ref)(iv) for details). Thus, we focus attention on two major issues that remain: (i) how the high-level conditions imposed in Assumptions (ref) and (ref), and Condition (ref) in Theorem (ref), can be verified under intuitive primitive assumptions; and (ii) how to implement uniform inference based on our theory in Section (ref).
This example considers generalized conditional quantile regression with a possibly non-identity link: $\rho(y,\eta;q)=(q-\I(y<\eta))(y-\eta)$, where $q\in\mathcal{Q}$ denotes the quantile position. Thus, let $\eta(\mu_0(\bx, q))$ be the conditional $q$-quantile of $Y$ given $\bX = \bx$; we verify in the supplemental appendix that such $\mu_0$ solves (ref). For this example, the following simple proposition, proven in the supplemental appendix (Proposition \saref{prop:gl-holder-quantile-regression}), gives sufficient conditions to verify the general Assumptions (ref) and (ref), and Condition (ref) in Theorem (ref).
The additional conditions in this proposition are primitive and easy-to-interpret, only restricting the conditional density of $Y$ given $\bX$ to be bounded and smooth in a mild sense. Our assumptions are on par with or are weaker than those imposed in Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE, despite the high level of generality of our theoretical results.
We can implement uniform inference following the plug-in method described in Section (ref). In this context $S_{q, \tilde{q}}(\bx) = q \wedge \tilde{q} - q \tilde{q}$ is known and constant in $\bx$, so a natural plug-in estimator of $\bar\bSigma_{q, \tilde{q}}$ is
On the other hand, the matrix $\bar\bQ_q$
depends on the unknown conditional density $f_{Y|X}$, and a plug-in estimator is not immediately available. However, many estimation strategies have been proposed in the literature Koenker_2005_book. We do not recommend a particular choice, but rather any estimator satisfying the mild convergence rate requirement in Condition (iii) of Theorem (ref) may be used.
The loss function is $\rho(y, \eta;q) = ( \I(y \leq q) - \eta )^2$ with a possibly non-identity inverse link function $\eta(\cdot)$. The derivative function is $\psi(y,\eta;q)=-2(\I(y\leq q)-\eta)$. The following proposition, proven in the supplemental appendix (Proposition \saref{prop:gl-holder-distribution-regression}), verifies our high-level assumptions under mild regularity conditions on the conditional distribution function of $Y$ given $\bX$.
The implementation of uniform inference follows the plug-in method described in Section (ref). To construct the prerequisite estimators, in this case $S_{q, \tilde{q}}(\bx_i) = 4 F_{Y|X}(q \wedge \tilde{q} | \bx_i) \big(1 - F_{Y|X}(q \vee \tilde{q} | \bx_i)\big)$. Therefore, a simple plug-in estimator of $\bar\bSigma_{q,\tilde{q}}$ is
In addition, a plug-in estimator of the matrix $\bar\bQ_q$ is $\widehat\bQ_q=2\E_n[(\eta^{(1)}(\widehat\mu(\bx_i,q)))^2 \bp(\bx_i)\bp(\bx_i)\trans]$.
The loss function is $\rho(y, \eta) = | y - \eta |^p$, $p \in(1,2]$ with a possibly non-identity link. The case $p=1$ is equivalent to quantile (median) regression discussed previously. The derivative function is $\psi(y,\eta)\equiv\psi(y-\eta)=p|y-\eta|^{p-1}\sgn(\eta-y)$. In this example the family $\mathcal{Q}$ of the loss functions is a singleton, and hence the dependence on the index $q$ can be dropped to simplify notation.
The following proposition, proven in the supplemental appendix (Proposition \saref{prop:lp-regression}), provides a set of simple regularity conditions that ensure our general theory can be applied to study generalized $L_p$ regression estimation and inference.
For implementation, we follow the plug-in method in Section (ref). Since $\mathcal{Q}$ is a singleton, dependence on $q$ can be dropped. Direct plug-in choices for estimating the prerequisite matrices take the form \[ \widehat\bQ=\E_n[\bp(\bx_i)\bp(\bx_i)\trans\widehat\Psi_{1, i} [\eta^{(1)}(\widehat\mu(\bx_i))]^2]\quad \text{and}\quad \widehat\bSigma=\E_n[\bp(\bx_i)\bp(\bx_i)\trans\psi(\widehat\epsilon_i)^2[\eta^{(1)}(\widehat\mu(\bx_i))]^2], \] where $\widehat\epsilon_i=y_i-\eta(\widehat\mu(\bx_i))$ and $\widehat\Psi_{1, i}$ is some estimator of the function $\Psi_1(\bx_i,\eta(\mu_0(\bx_i)))$. In $L_p$ regression with $p\in(1,2]$, $\Psi_1(\bx,\eta)=p(p-1)\E[|Y-\eta|^{p-2}\sgn(\eta-Y)|\bX=\bx]$, and therefore a simple plug-in choice is $\widehat\Psi_{1, i} =p(p-1)|y_i-\eta(\widehat\mu(\bx_i))|^{p-2}\sgn(\eta(\widehat\mu(\bx_i))-y_i)$. As an alternative, bootstrap-based inference could be used.
For this final example, the loss function is $\rho(y,\eta)=-y\log \eta-(1-y)\log (1-\eta)$, the inverse link function is $\eta(\theta)=1/(1+e^{-\theta})$, and the derivative function is $\psi(y,\eta)=- y / \eta + (1-y) / (1-\eta)$, and the loss function does not depend on $q\in\mathcal{Q}$. The following proposition, proven in the supplemental appendix (Proposition \saref{prop:gl-holder-logistic-regression}), gives simple primitive conditions verifying the high-level assumptions for our general theoretical results.
It is easy to construct a feasible Gaussian process $\widehat{Z}(\bx)$ conditional on the data $\bD_n$ with covariance structure (ref). Standard choices are
where $\widehat{\eta}_i=\eta(\widehat\mu(\bx_i))$ and $\widehat\epsilon_i=y_i-\widehat\eta_i$. See Section (ref) for more discussion.
We focused on uniform estimation and inference for the unknown function $\mu_0$ and derivatives thereof. However, the parameter of interest may be other linear or nonlinear transformations of $\mu_0$. For example, in generalized linear models usually the goal is to estimate the function $\eta(\mu_0(\bx, q))$, or the marginal effect of a covariate on that function $\frac{\partial}{\partial x_k}\eta(\mu_0(\bx,q))=\eta^{(1)}(\mu_0(\bx,q))\mu_0^{(\be_k)}(\bx,q)$. Furthermore, in treatment effect and causal inference settings Abadie-Cattaneo_2018_ARE, interest often lies in differences of such estimands across two or more subgroups: for two treatment levels $j=1,2$, $\eta(\mu_2(\bx,q))-\eta(\mu_1(\bx,q))$ can be interpreted as a mean, quantile, or other conditional (on $(\bx,q)$) treatment effect, where $\mu_j(\bx,q)$ is estimated using separately the subsample of, say, control ($j=1$) and treated ($j=2$) units. Our results can be applied to all these cases of practical interest with minimal additional effort.
We showcase the generality of our theory by briefly discussing uniform inference on the transformed function $\eta(\mu_0(\bx,q))$, its first derivative, and differences thereof across subgroups. Given the partitioning-based $M$-estimators $\widehat\mu(\bx,q)$ and $\widehat\mu_j(\bx,q)$, $j=1,2$, where $\widehat\mu_j$ is constructed using only data from the subsample $j$ of the full sample, we can immediately plug in to form the desired estimators.
Uniform consistency of the three estimators follows from uniform consistency of $\widehat\mu(\bx,q)$ (Corollary (ref)) because the transformation function $\eta$ is twice continuously differentiable. A Bahadur representation for each of the transformation estimators can be established via Theorem (ref) and a Taylor expansion. For example, for the level estimator,
with
and for the marginal effect of the $k$th covariate,
with
where the approximation remainders from the Taylor expansion, and their uniform rates $r_{\tt LE}$ and $r_{\tt ME}$, are precisely characterized in the supplemental appendix (Theorem \saref{th:other-parameters}). The conditional treatment effect estimator is simply a difference of two level estimators, each employing a disjoint sub-sample, and therefore it follows directly that
with $\mathsf{L}_{\tt CTE}(\bx,q) = \mathsf{L}_{{\tt LE},2}(\bx,q) - \mathsf{L}_{{\tt LE},1}(\bx,q)$ with $\mathsf{L}_{{\tt LE},j}(\bx,q)$ denoting the Bahadur approximation $\mathsf{L}_{\tt LE}(\bx,q)$ but when only using the sub-sample $j$.
Given the uniform Bahadur representations for each of the transformation estimators, strong approximations of their corresponding $t$-statistic processes can be constructed as in Section (ref). For example, conditional on $\bX_n$, the stochastic process $(\mathsf{L}_{\tt LE}(\bx,q) : (\bx,q)\in\mathcal{X}\times\mathcal{Q})$ has mean zero and variance $|\eta^{(1)}(\mu_0(\bx, q))|^2\bar\Omega_{\bm{0}}(\bx,q)/n$. Then, applying our strong approximation strategy, we can construct a conditional Gaussian process $Z_{\tt LE}(\bx, q)$ that approximates the $t$-statistic process of $\eta(\widehat{\mu}(\bx,q))$:
with strong approximation rate $r_{\tt SALE}$ as in Theorem (ref). Similarly, we can also construct a conditional Gaussian process $Z_{\tt ME}(\bx,q)$ that approximates the $t$-statistic process of the marginal effect estimator $\frac{\partial}{\partial x_k}\eta(\widehat\mu(\bx,q))$:
with strong approximation rate $r_{\tt SAME}$ as in Theorem (ref). These results are formalized in the supplemental appendix (Theorem \saref{th:other-parameters}). An analogous result holds for the conditional treatment effect estimator.
Finally, for implementation we can construct feasible processes to approximate $Z_{\tt LE}(\bx, q)$ and $Z_{\tt ME}(\bx, q)$ via the plug-in method discussed in Section (ref), and illustrated in Section (ref), which then can be employed to approximate the distributions of the entire level process $(\eta(\widehat\mu(\bx,q)):(\bx,q)\in\mathcal{X}\times\mathcal{Q})$, marginal effect process $(\frac{\partial}{\partial x_k}\eta(\widehat\mu(\bx,q)):(\bx,q)\in\mathcal{X}\times\mathcal{Q})$, and conditional treatment effect process $(\eta(\widehat\mu_2(\bx,q)) - \eta(\widehat\mu_1(\bx,q)) :(\bx,q)\in\mathcal{X}\times\mathcal{Q})$.
This paper investigated the asymptotic properties of a large class of nonparametric partitioning-based M-estimators, allowing for different degrees of non-smoothness in the loss function and a possibly non-identity monotonic transformation function. Our main theoretical results include uniform consistency for convex and non-convex objective functions, uniform Bahadur representations with optimal remainder under appropriate conditions, uniform and mean square convergence rates achieving optimal approximation under appropriate conditions, uniform strong approximation methods under general conditions, and uniform inference methods via plug-in approximations. We illustrated our general theory with four examples, and demonstrated how our results improve on prior literature, in many cases requiring minimal side rate restrictions on tuning parameters and achieving rate-optimal approximation rates. The supplemental appendix collects further theoretical results and generalizations that may be of independent interest. In future work, we plan to investigate optimal tuning parameter selection, including random partitioning schemes, and the validity of bootstrap-based approximations.
We thank Richard Crump, Max Farrell, Will Underwood, and Rae Yu for helpful comments and discussions. We also thank the Co-Editor, Associate Editor, and two reviewers for their comments. Cattaneo gratefully acknowledges financial support from the National Science Foundation through grants DMS-2210561 and SES-2241575. Feng gratefully acknowledges financial support from the National Natural Science Foundation of China (NSFC) through grants 72203122.
\begingroup \endgroup
\ifthenelse{\boolean{standalone}}{