The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
121,227 characters
Uniform Estimation and Inference for Nonparametric Partitioning-Based M-Estimators
\pagenumbering{roman}
\title{Uniform Estimation and Inference\\ for Nonparametric Partitioning-Based M-Estimators
\bigskip }
\author{Matias D. Cattaneo\thanks{Department of Operations Research and Financial Engineering, Princeton University.} \and
Yingjie Feng\thanks{School of Economics and Management, Tsinghua University.} \and
Boris Shigida\thanks{Department of Operations Research and Financial Engineering, Princeton University.}
}
\maketitle
\begin{abstract}
This paper presents uniform estimation and inference theory for a large class of nonparametric partitioning-based M-estimators. The main theoretical results include:
(i) uniform consistency for convex and non-convex objective functions;
(ii) rate-optimal uniform Bahadur representations;
(iii) rate-optimal uniform (and mean square) convergence rates;
(iv) valid strong approximations and feasible uniform inference methods; and
(v) extensions to functional transformations of underlying estimators.
Uniformity is established over both the evaluation point of the nonparametric functional parameter and a Euclidean parameter indexing the class of loss functions. The results also account explicitly for the smoothness degree of the loss function (if any), and allow for a possibly non-identity (inverse) link function. We illustrate the theoretical and methodological results in four examples: quantile regression, distribution regression, $L_p$ regression, and Logistic regression. Many other possibly non-smooth, nonlinear, generalized, robust M-estimation settings are covered by our results. We provide detailed comparisons with the existing literature and demonstrate substantive improvements: we achieve the best (in some cases optimal) known results under improved (in some cases minimal) requirements in terms of regularity conditions and side rate restrictions. The supplemental appendix reports complementary technical results that may be of independent interest, including a novel uniform strong approximation result based on Yurinskii's coupling.
\end{abstract}
\textit{Keywords}: nonparametric estimation and inference, series methods, partitioning estimators, quantile regression, nonlinear regression, robust regression, generalized linear models, uniform distribution theory.
\thispagestyle{empty}
\clearpage
\setcounter{tocdepth}{2}
\tableofcontents
\ifthenelse{\boolean{standalone}}{}{
\clearpage
\listoftheorems[
title={List of important statements},
ignore={lemma,remark,remark*},
numwidth=3em,
]
\clearpage
\listoftheorems[
title={List of lemmas},
ignoreall,
show={lem,lemma},
numwidth=3em,
]
}
\clearpage
\pagenumbering{arabic}
\pagestyle{plain}
\section{Introduction}\label{sec: Introduction}
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
\begin{equation}\label{eq: Rho problem}
\mu_{0}(\cdot,q) \in \argmin_{\mu\in\mathscr{M}} \E\big[\rho\big(y_{i} , \eta\big(\mu(\bx_{i})\big); q\big)\big],
\end{equation}
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.
\subsection{Partitioning-Based Methodology}
The series (or sieve) nonparametric partitioning-based $M$-estimator is
\begin{equation}\label{eq: M estimator}
\widehat{\mu}(\bx, q)=\bp(\bx)\trans
\widehat{\bbeta}(q),\qquad
\widehat\bbeta(q) \in \argmin_{\bb\in\mathcal{B}} \sum_{i=1}^n\rho\big(y_i,\eta(\bp(\bx_i)\trans \bb);q\big),
\end{equation}
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}$.
\cite{Gyorfi-etal_2002_book} give a textbook introduction to the partitioning-based estimation literature. In this nonparametric framework, the \textit{Haar} basis
\begin{align*}
\bp_{\mathtt{H}}(\bx) = \big( \I(\bx\in \delta_1), \ldots, \I(\bx\in \delta_{\kappa}) \big)\trans
\end{align*}
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 \eqref{eq: M estimator} 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 \textit{piecewise polynomial} basis
\begin{align*}
\bp_{\mathtt{P}}(\bx) = \bp_{\mathtt{H}}(\bx) \otimes \br_m(\bx),
\end{align*}
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 \textit{Splines} basis is
\begin{align*}
\bp_{\mathtt{S}}(\bx)=\otimes_{k=1}^d \bT_s\bp_{\mathtt{P}}(\be\trans_k\bx),
\end{align*}
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 \eqref{eq: M estimator} 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 \cite{Belloni-Chernozhukov-Chetverikov-Kato_2015_JoE}, \cite{Cattaneo-Farrell_2013_JoE}, \cite{Cattaneo-Farrell-Feng_2020_AOS}, and \cite{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 \citep{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$.
\begin{assumption}[Quasi-uniform partition]\label{Assumption: gl-quasi-uniform-partition}
The ratio of the sizes of inscribed and circumscribed balls of each $\delta \in \Delta$ is bounded away from zero uniformly in $\delta \in \Delta$, and
\begin{equation*}
\frac{\max \{\operatorname{diam}(\delta): \delta \in \Delta\}}{\min \{\operatorname{diam}(\delta): \delta \in \Delta\}} \lesssim 1
\end{equation*}
where $\operatorname{diam}(\delta)$ denotes the diameter of $\delta$. Further, for $h = \max \{\operatorname{diam}(\delta): \delta \in \Delta\}$, assume $h = o(1)$ and $\log (1/h)\lesssim \log n$, as $n\to\infty$.
\end{assumption}
Assumption \ref{Assumption: gl-quasi-uniform-partition} 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 \textit{active} on $\delta\in\Delta$ if it is not identically zero on $\delta$; we also employ standard multi-index notation (see Section \ref{sec: Notation} for details).
\begin{assumption}[Local basis]\label{Assumption: gl-local-basis}\leavevmode$ $
\begin{enumerate}[label=\emph{(\roman*)}]
\item For each basis function $p_k$, $k=1, \ldots, K$, the union of elements of $\Delta$ on which $p_k$ is active is a connected set, denoted by $\mathcal{H}_k$. For all $k=1, \ldots, K$, both the number of elements of $\mathcal{H}_k$ and the number of basis functions which are active on $\mathcal{H}_k$ are bounded by a constant.
\item For any $\ba=\left(a_1, \ldots, a_K\right)\trans \in \mathbb{R}^K$, $a_k^2 h^d \lesssim \ba\trans \int_{\mathcal{H}_k} \bp(\bx) \bp(\bx)\trans d \bx \ \ba$ for $k=1, \ldots, K$.
\item Let $\multindnorm{\bv} < m$. There exists an integer $\varsigma \in[\multindnorm{\bv}, m)$ such that, for all $\boldsymbol{\varsigma},\multindnorm{\boldsymbol{\varsigma}} \leq \varsigma$,
\begin{equation*}
h^{-\multindnorm{\boldsymbol{\varsigma}}}
\lesssim \inf _{\delta \in \Delta} \inf _{\bx \in \closure(\delta)}\big\|\bp^{(\boldsymbol{\varsigma})}(\bx)\big\|
\leq \sup _{\delta \in \Delta} \sup _{\bx \in \closure(\delta)}\big\|\bp^{(\boldsymbol{\varsigma})}(\bx)\big\| \lesssim h^{-\multindnorm{\boldsymbol{\varsigma}}},
\end{equation*}
where $\closure(\delta)$ is the closure of $\delta$.
\end{enumerate}
\end{assumption}
In Assumption \ref{Assumption: gl-local-basis}, 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{Assumption: gl-quasi-uniform-partition} and \ref{Assumption: gl-local-basis} 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{Assumption: dgp-holder} in Section \ref{sec: Assumptions}), the parameter $m$ will control how well $\mu_0$ can be approximated by linear combinations of the local basis (via Assumption \ref{Assumption: gl-approximation-error} in Section \ref{sec: Assumptions}). 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 \eqref{eq: Rho problem} and its associated plug-in estimator \eqref{eq: M estimator} 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 \eqref{eq: M estimator} has been computed, and then relies on the assumptions in Section \ref{sec: Assumptions} concerning the data generating process to study the large sample statistical properties of the partitioning-based $M$-estimator and transformations thereof.
\subsection{Summary of Contributions}
The objective function in \eqref{eq: M estimator} 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{Assumption: gl-approximation-error} in Section \ref{sec: Assumptions}). 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{sec: consistency}.
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 \eqref{eq: M estimator}:
\begin{enumerate}[label=(\roman*)]
\item rate-optimal Bahadur representation uniformly over $\mathcal{X}\times\mathcal{Q}$,
\item rate-optimal convergence rates in mean square and uniformly over $\mathcal{X}\times\mathcal{Q}$, and
\item valid strong approximation and feasible distribution theory uniformly over $\mathcal{X}\times\mathcal{Q}$.
\end{enumerate}
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{sec: Bahadur representation}, while results (iii) are presented in Sections \ref{sec: uniform inference} and \ref{sec: Feasible Uniform Inference}.
Section \ref{sec: Motivating Examples} introduces four examples: \textit{Generalized Conditional Quantile Regression}, \textit{Generalized Conditional Distribution Regression}, \textit{Generalized $L_p$ Regression Estimation}, and \textit{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{sec: Assumptions} 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{sec: Examples} demonstrates how our general sufficient conditions are verified for each of our motivating examples.
Section \ref{sec: Extensions} discusses how our results can be extended to cover other parameters of interest, while Section \ref{sec: conclusion} 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 \citep{chen2020jackknife,monrad1991nearby} and generalizing prior Yurinskii's coupling results in the literature \citep{Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE,yurinskii1978error}.
\subsection{Related Literature}
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, \cite{Eggermont-LaRiccia_2009_Book} and \cite{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, \cite{Belloni-Chernozhukov-Chetverikov-Kato_2015_JoE}, \cite{Cattaneo-Farrell_2013_JoE}, \cite{Cattaneo-Farrell-Feng_2020_AOS}, \cite{Cattaneo-Crump-Farrell-Feng_2024_AER}, \cite{Chen-Christensen_2015_JOE}, \cite{Huang_2003_AoS}, \cite{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 \eqref{eq: M estimator}, 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, \cite{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 \textit{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 \textit{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{Assumption: gl-quasi-uniform-partition} and \ref{Assumption: gl-local-basis}). Furthermore, our results on uniform consistency disentangling convex and non-convex loss functions (Section \ref{sec: consistency}), rate-optimal uniform Bahadur representation and convergence capturing explicitly the smoothness degree of the loss function (Section \ref{sec: Bahadur representation}), and uniform feasible inference validity (Sections \ref{sec: uniform inference} and \ref{sec: Feasible Uniform Inference}), 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 \cite{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{example: Generalized Conditional Quantile Regression} 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 \cite{Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE}: see Example \ref{example: Generalized Conditional Quantile Regression} in Sections \ref{sec: Motivating Examples} and \ref{sec: quantile reg}, and Sections \ref{sec: Comparison with Existing Results -- Bahadur} and \ref{sec: Comparison with Existing Results -- SA}. On the other hand, it is worth noting that \cite{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, \cite{Kong-Linton-Xia_2013_ET} considers local polynomial methods, and \cite{Shang-Cheng_2013_AOS} considers smoothing spline methods. Sections \ref{sec: Comparison with Existing Results -- Bahadur} and \ref{sec: Comparison with Existing Results -- SA} 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.
\subsection{Notation}\label{sec: Notation}
We employ standard notation in probability, statistics and empirical process theory \citep{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}$.
\section{Examples}\label{sec: Motivating Examples}
We discuss four motivating examples of interest covered by our theoretical results. Section \ref{sec: Examples} demonstrates how our high-level assumptions, introduced in Section \ref{sec: Assumptions}, 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 \cite{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.
\begin{example}[Generalized Conditional Quantile Regression]\label{example: Generalized Conditional Quantile Regression}
The quantile regression loss function is
\[ \rho(y,\eta;q)=(q-\mathds{1}(y<\eta))(y-\eta),\]
where $q\in\mathcal{Q}=[\varepsilon_0,1-\varepsilon_0]$ denotes the quantile position, with $\varepsilon_0>0$. Then, $\eta(\mu_0(\bx,q))$ is the $q$-th conditional quantile function of $Y$ given $\bX=\bx$, and the partitioning-based quantile regression estimator is $\eta(\widehat{\mu}(\bx,q))$ as defined in \eqref{eq: M estimator}. In the classical case, $\eta(\cdot)$ is the identity function, but our theory accommodates other transformations. Interest lies on the quantile process estimator $(\eta(\widehat{\mu}(\bx,q)):(\bx,q)\in\mathcal{X}\times\mathcal{Q})$, which can be used to characterize heterogeneous effects of covariates on the outcome distribution and to conduct specification testing. See Section \ref{sec: quantile reg} for our main results, and Section \ref{sec: Extensions} for results on transformations.
\end{example}
\cite{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.
\begin{example}[Generalized Conditional Distribution Regression]\label{example: Conditional Distribution Regression}
Non-linear least squares conditional distribution regression employs
\[\rho(y,\eta;q)= ( \mathds{1}(y \leq q) - \eta)^2,\]
where we can, for example, use the complementary log-log link $\eta(\theta)=1-\exp(-\exp(\theta))$. Estimands of interest are $\eta(\mu_0(\bx,q))$, which corresponds to the conditional distribution function of $Y$ given $\bX=\bx$ (i.e., $F_{Y|\bX}(q|\bx)=\E[\mathds{1}(Y \leq q)|\bX=\bx]$), and derivatives thereof. Uniform estimation and inference results based on $(\eta(\widehat{\mu}(\bx,q)):(\bx,q)\in\mathcal{X}\times\mathcal{Q})$ are useful for a variety of purposes, including heterogeneous treatment effect estimation and specification testing. See Section \ref{sec: distribution reg} for our main results, and Section \ref{sec: Extensions} for transformations.
\end{example}
\cite{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)$.
\begin{example}[Generalized \texorpdfstring{$L_p$}{Lp} Regression]\label{example: $L_p$ regression}
The (possibly nonlinear) $L_p$ regression estimator is defined by taking
\[\rho(y, \eta)=|y-\eta|^{p},\]
for a fixed $p>0$. In particular, $p=2$ leads to nonlinear least squares, and $p=1$ leads to generalized least absolute deviations, when a non-identity transformation function $\eta(\cdot)$ is used. The estimand of interest is usually the transformed regression function $\eta(\mu_0(\bx))$, which needs to be interpreted in context. Our general theory applies to any choice $p\in[1,2]$, delivering uniform estimation and inference methods based on $(\widehat{\mu}(\bx):\bx\in\mathcal{X})$, and transformations thereof. See Section \ref{sec: lp reg} for our main results, and Section \ref{sec: Extensions} for transformations.
\end{example}
The final example considers (nonparametric) Generalized Linear Models \citep{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).
\begin{example}[Generalized Linear Models]\label{example: Generalized linear models}
The classical logistic regression model, or binary classification with sigmoid (inverse) link, employs
\[\rho(y,\eta)=-y\log \eta-(1-y)\log (1-\eta),\qquad \eta(\theta)=\exp(\theta)/(1+\exp(\theta)),\]
with $\mathcal{Y}=\{0,1\}$. The estimand $\eta(\mu_0(\bx))$ characterizes the conditional probability of $Y=1$ given $\bX=\bx$. See Section \ref{sec: logit est} for uniform estimation and inference methods based on $(\widehat{\mu}(\bx):\bx\in\mathcal{X})$, and Section \ref{sec: Extensions} for transformations thereof. Furthermore, our results also cover other related quasi-maximum likelihood (and non-linear least squares) problems such as fractional regression where $\mathcal{Y}=[0,1]$.
\end{example}
The four examples illustrate distinct settings from a technical perspective. In Example \ref{example: Generalized Conditional Quantile Regression} 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{example: Conditional Distribution Regression} is the antithesis of Example \ref{example: Generalized Conditional Quantile Regression} 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{example: $L_p$ regression} 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{example: Generalized linear models} 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.
\section{Assumptions}\label{sec: Assumptions}
Our theoretical work proceeds under Assumptions \ref{Assumption: gl-quasi-uniform-partition} and \ref{Assumption: gl-local-basis} 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.
\begin{assumption}[Data Generating Process]
\label{Assumption: dgp-holder}\leavevmode
\begin{enumerate}[label=\emph{(\roman*)}]
\item $((y_i, \bx_i) \colon 1 \leq i \leq n)$ is a random sample satisfying~\eqref{eq: Rho problem}.
\item The distribution of $\bx_i$ admits a Lebesgue density $f_{X}(\cdot)$ which is continuous and bounded away from zero on support $\mathcal{X} \subset \mathbb{R}^d$, where $\mathcal{X}$ is the closure of an open, connected and bounded set.
\item The conditional distribution of $y_{i}$ given $\bx_i$ admits a conditional density $f_{Y|X}(y | \bx)$ with support $\mathcal{Y}_{\bx}$ with respect to some (sigma-finite) measure $\mathfrak{M}$, and $\underset{\bx \in \mathcal{X}}{\sup}\, \underset{y \in \mathcal{Y}_{\bx}}{\sup}\, f_{Y|X}(y \mid \bx) < \infty$.
\item $\bx\mapsto\mu_{0}(\bx,q)$ is $m \geq 1$ times continuously differentiable for every $q\in\mathcal{Q}$, $\bx\mapsto\mu_{0}(\bx,q)$ and its derivatives of order no greater than $m$ are bounded uniformly over $(\bx,q) \in \mathcal{X}\times\mathcal{Q}$, and
\begin{align*}
\sup_{\bx\in\mathcal{X}}\sup_{q_1\neq q_2} \frac{|\mu_{0}(\bx,q_1) - \mu_{0}(\bx,q_2)|}{|q_1 - q_2|}
\lesssim 1.
\end{align*}
\end{enumerate}
\end{assumption}
Assumption \ref{Assumption: dgp-holder} 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$.
\begin{assumption}[Loss Function]\label{Assumption: loss function}\leavevmode
\begin{enumerate}[label=\emph{(\roman*)}]
\item Let $\mathcal{Q}\subseteq\mathbb{R}$ be a connected compact set. For each $q \in \mathcal{Q}$, $y \in \mathcal{Y}$, and some open connected subset $\mathcal{E}$ of $\mathbb{R}$ not depending on $y$, $\eta\mapsto\rho(y,\eta;q)$ is absolutely continuous on closed bounded intervals within $\mathcal{E}$, and admits an a.\,e. derivative $\psi(y,\eta; q)$.
\item The first-order optimality condition $\E[\psi(y_i, \eta(\mu_{0}(\bx_i,q));q)| \bx_i] =0$ holds; the function $\E[\psi(y_i, \eta(\mu_{0}(\bx_i,q));q)^{2} | \bx_i = \bx]$ is continuous in both arguments $(\bx,q) \in \mathcal{X}\times\mathcal{Q}$, bounded away from zero, and Lipschitz in $q$ uniformly in $\bx$; there is a positive measurable envelope function $\overline{\psi}(\bx_i, y_i)$ such that $\sup_{q\in\mathcal{Q}} |\psi(y_i, \eta(\mu_{0}(\bx_i,q));q)| \leq \overline{\psi}(\bx_i, y_i)$ with $\sup_{\bx \in \mathcal{X}}\E[\overline{\psi}(\bx_i, y_i)^{\nu} | \bx_i = \bx] < \infty$ for some $\nu > 2$.\sloppy
\item\label{shape of derivative} $\eta(\cdot)$ is strictly monotonic and twice continuously differentiable. Furthermore, for some fixed constant $\alpha \in (0, 1]$, for any $(\bx,q) \in \mathcal{X}\times\mathcal{Q}$, and a pair of points $\zeta_1, \zeta_2 \in B_{q}(\bx)$, $\psi(\cdot)$ satisfies the following (constants hidden in $\lesssim$ do not depend on $\bx$, $q$, $\zeta_1$, $\zeta_2$):
\begin{itemize}
\item if $\mathfrak{M}$ is the Lebesgue measure, then
\begin{align*}
&\sup_{\lambda \in [0, 1]} \sup_{y \not\in [\eta(\zeta_1)\wedge \eta(\zeta_2), \eta(\zeta_1)\vee \eta(\zeta_2)]}\; | \psi ( y, \eta(\zeta_1 + \lambda(\zeta_2 - \zeta_1));q) -
\psi ( y, \eta(\zeta_2);q) | \lesssim | \zeta_1 - \zeta_2 |^{\alpha},\\
&\sup_{\lambda \in [0, 1]} \sup_{y \in [\eta(\zeta_1)\wedge \eta(\zeta_2), \eta(\zeta_1)\vee \eta(\zeta_2)]}\, | \psi ( y, \eta(\zeta_1 + \lambda(\zeta_2 - \zeta_1));q ) - \psi ( y, \eta(\zeta_2); q) | \lesssim 1;
\end{align*}
\item if $\mathfrak{M}$ is \emph{not} the Lebesgue measure, then
\begin{align*}
\sup_{\lambda \in [0, 1]} \sup_{y\in\mathcal{Y}}| \psi(y,\eta(\zeta_1 + \lambda(\zeta_2 - \zeta_1));q) - \psi(y,\eta(\zeta_2);q)|\lesssim |\zeta_1-\zeta_2|^\alpha.
\end{align*}
\end{itemize}
\item $\Psi(\bx,\eta;q)= \E[\psi(y_{i}, \eta;q) | \bx_i=\bx]$ is twice continuously differentiable with respect to $\eta$,
\begin{align*}
\sup_{\bx \in \mathcal{X}, q \in \mathcal{Q}} \sup_{\zeta \in B_{q}(\bx)} |\Psi_k(\bx, \eta(\zeta);q)| < \infty, \quad \Psi_k(\bx, \eta; q)= \frac{\partial^k}{\partial \eta^k} \Psi(\bx, \eta; q),\quad k=1,2,
\end{align*}
and
\begin{align*}
&\inf_{\bx \in \mathcal{X}, q \in \mathcal{Q}} \inf_{\zeta \in B_{q}(\bx)} \Psi_1(\bx, \eta(\zeta);q)(\eta^{(1)}(\zeta))^{2} > 0.
\end{align*}
\end{enumerate}
\end{assumption}
Assumption \ref{Assumption: loss function} is carefully crafted to accommodate all the examples discussed in Section \ref{sec: Motivating Examples}, 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 \eqref{eq: Rho problem}, 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, \cite{Lin-Kulasekera_2007_Biometrika}, \cite{matzkin2007nonparametric}, and references therein.
Part (iii) of Assumption \ref{Assumption: loss function} 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{sec: Examples}, where $\alpha=1$ in Examples \ref{example: Generalized Conditional Quantile Regression} and \ref{example: Conditional Distribution Regression}, $\alpha=p-1$ in Example \ref{example: $L_p$ regression} when $p>1$, and $\alpha=1$ in Example \ref{example: Generalized linear models}. 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{Assumption: loss function} collects mild regularity conditions on the smoothed-out a.\,e. derivative of the loss function.
Assumptions \ref{Assumption: dgp-holder} and \ref{Assumption: loss function} 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{sec: Motivating Examples}: 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.
\begin{assumption}[Simplified Setup]\label{Assumption: simple setup} $ $
\begin{enumerate}[label=\emph{(\roman*)}]
\item $q\mapsto\mu_0(\bx, q)$ is non-decreasing.
\item $\rho(y,\eta;q)=\sum_{j=1}^4\rho_j(y,\eta;q)$, where the functions $\rho_j(\cdot)$ are of the following types:
\begin{itemize}
\item Type I: $\rho_1(y,\eta;q)=(f_1(y)+D_1\eta)\I(y\leq \eta)$,
\item Type II: $\rho_2(y,\eta;q)=(f_2(y)+D_2\eta)\I(y\leq q)$,
\item Type III: $\rho_3(y,\eta;q)=(f_3(y)+D_3\eta)q$,
\item Type IV: $\rho_4(y,\eta;q)=\mathcal{T}(y,\eta)$,
\end{itemize}
where $f_{j}(\cdot)$ are fixed continuous functions, $D_j$ are universal constants, and $\mathcal{T}(\cdot)$ is a fixed continuous function.
\item $\eta\mapsto\E[\tau(y_i,\eta)|\bx_i=\bx]$ is differentiable, where $\tau(y,\eta) = \frac{\partial}{\partial \eta}\mathcal{T}(y,\eta)$, $\tau(y,\eta)$ and $\frac{\partial}{\partial \eta}\E[\tau(y_i,\eta)|\bx_i=\bx]$ are continuous in their arguments and $\alpha$-H\"{o}lder continuous ($\alpha\in(0,1]$) in $\eta$ for $\eta$ in any fixed compact subset of $\mathcal{E}$ with the H\"{o}lder constants independent of $(y,\bx)$.
\item $\sup_{q\in\mathcal{Q}} |\tau(y,\eta(\mu_{0}(\bx,q)))|\leq \bar\tau(\bx,y)$ with $\sup_{\bx \in \mathcal{X}} \E[\bar{\tau}(\bx_i, y_i)^{\nu} | \bx_i = \bx] < \infty$ for some $\nu > 2$.
\item If $D_1 \neq 0$, then $F_{Y|X}$ is differentiable with a Lebesgue density $f_{Y|X}$, and $f_{Y|X}$ is continuous in both arguments and $y\mapsto f_{Y|X}(y|\bx)$ is continuously differentiable.
\end{enumerate}
\end{assumption}
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{sec: Motivating Examples}. Finally, the key restriction emerging from Assumption \ref{Assumption: simple setup} 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{example: Generalized Conditional Quantile Regression} 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{example: Conditional Distribution Regression} is a combination of Type II and Type IV functions; the loss function in Example \ref{example: $L_p$ regression} is of Type IV for $p>1$, and of the same type as Example \ref{example: Generalized Conditional Quantile Regression} when $p=1$ (median regression); and the loss function in Example \ref{example: Generalized linear models} is usually a Type IV function. See Section \ref{sec: Examples} for details.
Our final assumption concerns the approximation power of the basis $\bp(\cdot)$ in connection with the underlying functional parameter.
\begin{assumption}[Approximation Error]\label{Assumption: gl-approximation-error}
There exists a vector of coefficients $\bbeta_0(q) \in \mathbb{R}^K$ such that for all $\boldsymbol{\varsigma}$ satisfying $\multindnorm{\boldsymbol{\varsigma}} \leq \varsigma$ in Assumption~\ref{Assumption: gl-local-basis},
\begin{equation*}
\sup_{q\in\mathcal{Q}}\sup_{\bx \in \mathcal{X}}
\big| \mu_0^{(\boldsymbol{\varsigma})}(\bx,q) -
\bbeta_0(q) \trans \bp^{(\boldsymbol{\varsigma})}(\bx) \big|
\lesssim h^{m - \multindnorm{\boldsymbol{\varsigma}}}.
\end{equation*}
\end{assumption}
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{Assumption: dgp-holder}(iv)). For more discussion, see \cite{Belloni-Chernozhukov-Chetverikov-Kato_2015_JoE}, \cite{Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE}, \cite{Cattaneo-Farrell_2013_JoE}, \cite{Cattaneo-Farrell-Feng_2020_AOS}, \cite{Chen-Christensen_2015_JOE}, \cite{Huang_2003_AoS}, and references therein.
\section{Consistency}\label{sec: consistency}
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 \eqref{eq: M estimator}: 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}$.
\subsection{Convex Loss Function}
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 \eqref{eq: M estimator}) under mild conditions. The proof is deferred to the supplemental appendix (Lemma~\saref{lem:holder-consistency}).
\begin{lem}[Consistency, convex case]\label{lem:consistency, convex}
Suppose that Assumptions \ref{Assumption: gl-quasi-uniform-partition}--\ref{Assumption: gl-approximation-error} hold,
$\rho(y, \eta(\theta); q)$ is convex with respect to $\theta$ with left or right derivative $\psi(y, \eta(\theta); q) \eta^{(1)}(\theta)$, $\mathcal{B}=\mathbb{R}^K$ in \eqref{eq: M estimator}, and $m > d / 2$. Furthermore, assume that one of the following two conditions holds:
\begin{enumerate}[label=\emph{(\roman*)}]
\item $\frac{( \log n )^{\frac{\nu}{\nu - 1}}}{n h^{\frac{2 \nu}{\nu - 1} d}} = o(1)$, or
\item $\frac{(\log n)^{3/2}}{nh^{2d}}=o(1)$ and $\overline\psi(\bx_i,y_i)$ is sub-Gaussian conditional on $\bx_i$.
\end{enumerate}
Then
\begin{align}
\sup_{q \in \mathcal{Q}}\; \big\| \widehat{\bbeta}(q) - \bbeta_0(q) \big\| &= o_{\P}(1),\label{eq: consistency, beta}\\
\sup_{\bx \in \mathcal{X}} \sup_{q \in \mathcal{Q}}\; \big|\widehat{\mu}^{(\bv)}(\bx,q) - \mu_{0}^{(\bv)}(\bx,q) \big| &= o_{\P}(h^{- \multindnorm{\bv}}),\label{eq: consistency, sup}\\
\sup_{q \in \mathcal{Q}}\; \int \big(\widehat\mu^{(\bv)}(\bx,q) - \mu_{0}^{(\bv)}(\bx,q) \big)^2 f_X(\bx)\, \mathrm{d} \bx &= o_{\P} ( h^{d - 2\multindnorm{\bv}} )\label{eq: consistency, l2}.
\end{align}
\end{lem}
Lemma~\ref{lem:consistency, convex} 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 \eqref{eq: consistency, beta}, the uniform-in-$q$ consistency in the Euclidean norm.
Two kinds of rate restrictions are imposed in Lemma \ref{lem:consistency, convex}, 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{lem:consistency, convex} is that \textit{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 \citep[e.g.,][]{chernozhukov-imbens-newey_2007_JoE}.
\subsection{Non-Convex Loss Function}
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 \eqref{eq: M estimator} 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 \textit{fixed} ``box'' constraint: for some fixed constant $R>0$,
\begin{equation*}
\widehat{\bbeta}(q) \in \argmin_{\| \bb \|_{\infty} \leq R}\; \sum_{i=1}^n
\rho ( y_i, \eta(\bp(\bx_i) \trans \bb);q ).
\end{equation*}
In the supplemental appendix we show that the pseudo-true coefficients $\bbeta_0(q)$ from Assumption \ref{Assumption: gl-approximation-error} 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.
\begin{lem}[Consistency, non-convex case]\label{lem:consistency, nonconvex}
Suppose that Assumptions \ref{Assumption: gl-quasi-uniform-partition}--\ref{Assumption: gl-approximation-error} hold, $\mathcal{B}=\{\bb\in\mathbb{R}^K: \|\bb\|_\infty\leq R\}$ with $R\geq 2\sup_{q\in\mathcal{Q}}\|\bbeta_0(q)\|_\infty$ in \eqref{eq: M estimator},
$m > d / 2$, and that there exists some constant $c>0$ such that $\inf \Psi_{1}(\bx,\zeta;q)>c$, where the infimum is over $\bx \in \mathcal{X}$, $q \in \mathcal{Q}$, $\zeta$ between $\eta(\bp(\bx)\trans\bbeta)$ and $\eta( \mu_0(\bx, q))$,
and $\bbeta \in \mathcal{B}$.
Furthermore, assume one of the following two conditions holds:
\begin{enumerate}[label=\emph{(\roman*)}]
\item $\frac{\left( \log n \right)^{\frac{\nu}{\nu - 1}}}{n h^{\frac{2 \nu}{\nu - 1} d}} = o(1)$, or
\item $\frac{(\log n)^{3/2}}{nh^{2d}}=o(1)$ and $\overline\psi(\bx_i,y_i)$ is sub-Gaussian conditional on $\bx_i$.
\end{enumerate}
Then \eqref{eq: consistency, beta}, \eqref{eq: consistency, sup}, and \eqref{eq: consistency, l2} hold.
\end{lem}
Compared to Lemma \ref{lem:consistency, convex}, 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 \textit{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.
\subsection{Weaker Conditions for Special Cases}\label{sec: Weaker Conditions for Special Cases}
In the supplemental appendix we provide additional consistency results for two special cases:
\begin{itemize}
\item the loss function is strongly convex and smooth (i.e., the second ``derivative'' of $\rho(y, \eta(\cdot); q)$ is bounded and bounded away from zero), or
\item an unconnected basis (i.e., each basis function is supported on a single cell of $\Delta$) is employed.
\end{itemize}
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{lem:consistency, convex} and \ref{lem:consistency, nonconvex} 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.
\subsection{Comparison with Existing Results}\label{sec: Comparison with Existing Results -- Consistency}
These restrictions imposed in Lemmas \ref{lem:consistency, convex} and \ref{lem:consistency, nonconvex}, 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 \citep{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{lem:consistency, convex} and \ref{lem:consistency, nonconvex} 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 \citep{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 \cite{Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE} implies that the $L_2$-consistency \eqref{eq: consistency, l2} can be obtained under $1/(nh^{2d}) \asymp K^2/n = o(1)$, and their Corollary 2 implies that the uniform consistency \eqref{eq: consistency, sup} 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.
\section{Bahadur Representation and Convergence Rates}\label{sec: Bahadur representation}
The Bahadur representation is
\begin{align*}
\mathsf{L}^{(\bv)}(\bx,q)
= -\bp^{(\bv)}(\bx)\trans \bQ_{0,q}^{-1}
\E_{n}\big[\bp(\bx_i) \eta^{(1)}(\mu_0(\bx_i,q)) \psi(y_i,\eta(\mu_{0}(\bx_i,q)); q)\big]
\end{align*}
where
\begin{align*}
\bQ_{0,q} = \E\big[\bp(\bx_i)\bp(\bx_i)\trans \Psi_1(\bx_i,\eta(\mu_0(\bx_i,q));q)[\eta^{(1)}(\mu_0(\bx_i,q))]^2\big].
\end{align*}
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{sec: consistency}. The proof is provided in the supplemental appendix (Theorem~\saref{th:bahadur-repres}).
\begin{thm}[Bahadur representation]\label{theorem: bahadur-repres}
Suppose that Assumptions \ref{Assumption: gl-quasi-uniform-partition}--\ref{Assumption: gl-approximation-error} hold.
Furthermore, assume the following four conditions:
\begin{enumerate}[label=\emph{(\roman*)}]
\item $\sup_{q \in \mathcal{Q}}\big\| \widehat{\boldsymbol{\beta}}(q) - \bbeta_0(q) \big\|_{\infty} = o_{\P}(1)$;
\item there exists a fixed constant $c>0$ such that
$\{\bb\in\mathbb{R}^K:\|\bb-\bbeta_0(q)\|_\infty\leq c, q\in\mathcal{Q}\}\subseteq\mathcal{B}$;
\item $\frac{(\log n)^{d + 2}}{nh^d}=o(1)$;
\item either $\frac{(\log n)^d}{n^{1-2/\nu} h^d} = o(1)$ or
$\overline{\psi}(\bx_i,y_i)$ is sub-Gaussian conditional on $\bx_i$.
\end{enumerate}
Then
\begin{small}
\begin{align}\label{eq: Bahadur rep, disc}
\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}
\big|\widehat{\mu}^{(\bv)}(\bx,q) - \mu_{0}^{(\bv)}(\bx,q) - \mathsf{L}^{(\bv)}(\bx,q) \big|
\lesssim_{\mathbb{P}}
h^{-\multindnorm{\bv}}
\Big(\frac{(\log n)^d}{nh^d}\Big)^{\frac{1}{2}+(\frac{\alpha}{2}\wedge \frac{1}{4})} \log n +
h^{m - \multindnorm{\bv}}.
\end{align}
\end{small}
If, in addition, $\sup_{y\in\mathcal{Y},q\in\mathcal{Q}} | \psi(y,\eta(\zeta_1);q) - \psi(y,\eta(\zeta_2);q)| \lesssim |\zeta_1-\zeta_2|^\alpha$ without any restrictions on $y$ in Assumption~\ref{Assumption: loss function}\ref{shape of derivative}, then
\begin{small}
\begin{align}\label{eq: Bahadur rep, cont}
\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}
\big|\widehat\mu^{(\bv)}(\bx,q) - \mu_{0}^{(\bv)}(\bx,q) - \mathsf{L}^{(\bv)}(\bx,q) \big|
\lesssim_{\mathbb{P}} h^{- \multindnorm{\bv}} \Big(\frac{(\log n)^d}{n h^d}\Big)^{\frac{1 + \alpha}{2}} \log n +
h^{m - \multindnorm{\bv}}.
\end{align}
\end{small}
\end{thm}
The Bahadur representation \eqref{eq: Bahadur rep, disc} 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{example: Generalized Conditional Quantile Regression}), 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 \eqref{eq: Bahadur rep, disc} 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{example: $L_p$ regression}, 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, \eqref{eq: Bahadur rep, disc} 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 \eqref{eq: Bahadur rep, cont} applies to the case where the ``derivative'' of the loss is a continuous function of $(y-\eta)$. Nonlinear least squares regression (Example \ref{example: Conditional Distribution Regression} and quasi-maximum likelihood estimation of generalized linear models (Example \ref{example: Generalized linear models}) fall into this category with the H\"{o}lder parameter $\alpha=1$. In such cases, \eqref{eq: Bahadur rep, cont} 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 \eqref{eq: Bahadur rep, disc}. See Section \ref{sec: Examples} 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{Assumption: gl-approximation-error}), 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 \eqref{eq: Bahadur rep, disc} is greater than that in the continuous case \eqref{eq: Bahadur rep, cont}; with a smaller H\"{o}lder parameter $\alpha$, the order of the remainder in both cases could increase.
\subsection{Rates of Convergence}\label{sec: rates of convergence}
The uniform Bahadur representations (Theorem \ref{theorem: bahadur-repres}) 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}$.
\begin{coro}[Uniform Rate of Convergence]\label{coro: convergence}
Suppose that Assumptions \ref{Assumption: gl-quasi-uniform-partition}--\ref{Assumption: gl-approximation-error} and the four conditions (i)--(iv) in Theorem \ref{theorem: bahadur-repres} hold. Furthermore, assume one of the following two conditions holds:
\begin{enumerate}[label=\emph{(\roman*)}]
\item $\frac{(\log n)^{d+\frac{d+1}{\alpha\wedge 0.5}}}{nh^d}=O(1)$, and $h^{(\alpha\wedge 0.5)m}(\log n)^{0.5 d}=O(1)$, or
\item the additional condition for \eqref{eq: Bahadur rep, cont} holds, $\frac{(\log n)^{d+\frac{d+1}{\alpha}}}{nh^d}=O(1)$, and $h^{\alpha m}(\log n)^{0.5d}=O(1)$.
\end{enumerate}
Then
\begin{equation}\label{eq: uniform convergence}
\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}
\big|\widehat{\mu}^{(\bv)}(\bx,q)-\mu_0^{(\bv)}(\bx,q)\big|
\lesssim_\P
h^{-|\bv|}\sqrt{\frac{\log n}{nh^d}}+h^{m-|\bv|}.
\end{equation}
\end{coro}
By setting $h\asymp \big(\frac{\log n}{n}\big)^{\frac{1}{2m+d}}$, Corollary \ref{coro: convergence} 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 \citep{Stone_1982_AoS} and conditional quantiles \citep{chaudhuri1991global}. In this sense, the rate of convergence in Corollary \ref{coro: convergence} is optimal and cannot be further improved at our level of generality.
Theorem \ref{theorem: bahadur-repres} can also be used to obtain the mean square convergence rate of the partitioning-based $M$-estimator uniformly-in-$q$.
\begin{coro}[Mean Square Rate of Convergence]\label{coro: Mean Square Rate of Convergence}
Suppose that Assumptions \ref{Assumption: gl-quasi-uniform-partition}--\ref{Assumption: gl-approximation-error} and the four conditions (i)--(iv) in Theorem \ref{theorem: bahadur-repres} hold. Furthermore, assume one of the following two conditions holds:
\begin{enumerate}[label=\emph{(\roman*)}]
\item $\frac{(\log n)^{d+\frac{d+2}{\alpha\wedge 0.5}}}{nh^d}=o(1)$ and $h^{(\alpha\wedge 0.5)m}(\log n)^{\frac{d+1}{2}}=o(1)$, or
\item the additional condition for \eqref{eq: Bahadur rep, cont} holds, $\frac{(\log n)^{d+\frac{d+2}{\alpha}}}{nh^d}=o(1)$, and $h^{\alpha m}(\log n)^{\frac{d+1}{2}}=o(1)$.
\end{enumerate}
Then
\begin{equation}\label{eq: l2 convergence}
\sup_{q\in\mathcal{Q}}\int_{\mathcal{X}} \big(\widehat{\mu}^{(\bv)}(\bx,q)-\mu_0^{(\bv)}(\bx,q)\big)^2 f_X(\bx)d\bx
\lesssim_\P \frac{1}{nh^{d+2|\bv|}}+h^{2(m-|\bv|)}.
\end{equation}
\end{coro}
By setting $h\asymp n^{-\frac{1}{2m+d}}$, Corollary \ref{coro: Mean Square Rate of Convergence} 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 \citep{Stone_1980_AoS} and conditional quantiles \citep{chaudhuri1991global}.
The convergence rates in \eqref{eq: uniform convergence} and \eqref{eq: l2 convergence} 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{coro: Mean Square Rate of Convergence}, it is possible to further leverage Theorem \ref{theorem: bahadur-repres} 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 \cite{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.
\subsection{Comparison with Existing Results}\label{sec: Comparison with Existing Results -- Bahadur}
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{theorem: bahadur-repres}, and the implied convergence rates in Corollaries \ref{coro: convergence} and \ref{coro: Mean Square Rate of Convergence}. 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{sec: consistency}. 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.
\subsubsection*{Mean Regression}
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 \cite{Belloni-Chernozhukov-Chetverikov-Kato_2015_JoE} and \cite{Cattaneo-Farrell-Feng_2020_AOS}. Since the derivative of the square loss for mean regression is linear, the first term in \eqref{eq: Bahadur rep, disc} or \eqref{eq: Bahadur rep, cont} does not show up in the uniform linearization of least squares series estimators. See, for example, Lemma SA-4.2 of \cite{Cattaneo-Farrell-Feng_2020_AOS}; $R_{1n,q}$ defined therein has been implicitly included in the leading variance term in \eqref{eq: Bahadur rep, cont} above. Theorem \ref{theorem: bahadur-repres} substantially extends these prior results to other nonlinear settings, under minimal additional conditions.
Finally, Corollaries \ref{coro: convergence} and \ref{coro: Mean Square Rate of Convergence} demonstrate the convergence rate optimality of general partitioning-based series $M$-estimation, recovering in particular known results for mean regression \citep{Belloni-Chernozhukov-Chetverikov-Kato_2015_JoE,Cattaneo-Farrell-Feng_2020_AOS} under essentially the same minimal conditions.
\subsubsection*{Quantile Regression}
Theorem \ref{theorem: bahadur-repres} improves upon prior theoretical results for nonparametric series quantile regression estimators. The most recent advance in this literature is due to \cite{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, \cite{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{theorem: bahadur-repres} 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 \cite{Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE}. In contrast, the restriction on $h$ in Theorem \ref{theorem: bahadur-repres} 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 \cite{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{theorem: bahadur-repres} 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{lem:consistency, convex} and \ref{lem:consistency, nonconvex}. 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 \textit{no} constraint on the relation between smoothness $m$ and dimensionality $d$; see Section \ref{sec: Weaker Conditions for Special Cases}.
Fourth, compared to \cite{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{theorem: bahadur-repres} by virtue of our general consistency results in Lemma \ref{lem:consistency, nonconvex}.
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{coro: convergence} and \ref{coro: Mean Square Rate of Convergence} 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, \cite{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$).
\subsubsection*{Other Nonparametric Smoothing Methods}
\cite{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{Assumption: loss function}: 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 \eqref{eq: Bahadur rep, disc} is $O((nh^d)^{-3/4})$, up to $\polylog(n)$ terms, and thus Theorem \ref{theorem: bahadur-repres} 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{theorem: bahadur-repres} establishes that the partitioning-based series $M$-estimator in \eqref{eq: M estimator} can achieve the same best Bahadur approximation as local polynomial kernel methods, up to $\polylog(n)$ terms.
Furthermore, compared to \cite{Kong-Linton-Xia_2013_ET} or other similar contributions in the literature, Theorem \ref{theorem: bahadur-repres} exhibits (at least) two novel features. First, the Bahadur representations \eqref{eq: Bahadur rep, disc} and \eqref{eq: Bahadur rep, cont} 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 \textit{entire} conditional quantile process may be of interest. Second, Theorem \ref{theorem: bahadur-repres} 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{example: Conditional Distribution Regression}) and $L_p$ regression with $p\in(1,2)$ (Example \ref{example: $L_p$ regression}). Finally, \cite{Kong-Linton-Xia_2013_ET} does not discuss convergence rates as we do in Corollaries \ref{coro: convergence} and \ref{coro: Mean Square Rate of Convergence}.
In the context of nonparametric penalized smoothing spline methods, \cite{Shang-Cheng_2013_AOS} also establishes a uniform Bahadur representation (and other results) that can be compared to Theorem \ref{theorem: bahadur-repres}. 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 \cite{Shang-Cheng_2013_AOS}) is $O((nh)^{-1} h^{-(6m-1)/(4m)})$, while \eqref{eq: Bahadur rep, cont} in Theorem \ref{theorem: bahadur-repres} gives the optimal result $O((nh)^{-1})$, thereby demonstrating a substantial improvement over their result. Finally, as for convergence rates, Proposition 3.3 in \cite{Shang-Cheng_2013_AOS} and our Corollary \ref{coro: Mean Square Rate of Convergence} are essentially equivalent, both delivering optimal mean square convergence. They do not explicitly discuss uniform convergence rates as we do in Corollary \ref{coro: convergence}.
\section{Strong Approximation}\label{sec: uniform inference}
The uniform Bahadur representations in Theorem \ref{theorem: bahadur-repres} can also be leveraged to establish uniform distribution theory for $\widehat\mu^{(\bv)}$. The infeasible conditional variance of the estimator can be written as
\begin{align*}
\bar\Omega_{\bv}(\bx,q) = \bp^{(\bv)}(\bx)\trans\bar\bQ_q^{-1}\bar\bSigma_q\bar\bQ_q^{-1}\bp^{(\bv)}(\bx),
\end{align*}
where
\begin{align*}
\bar\bQ_{q} &= \E_n\big[\bp(\bx_i)\bp(\bx_i)\trans \Psi_1(\bx_i,\eta(\mu_0(\bx_i,q));q)[\eta^{(1)}(\mu_0(\bx_i,q))]^2\big], \qquad \text{and}\\
\bar\bSigma_q &= \E_n\big[\bp(\bx_i)\bp(\bx_i)\trans\E[\psi(y_i,\eta(\mu_0(\bx_i,q));q)^2|\bx_i]
[\eta^{(1)}(\mu_0(\bx_i,q))]^2\big]
\end{align*}
Accordingly, we define a feasible variance estimator as
\begin{align*}
\widehat{\Omega}_{\bv}(\bx,q) = \bp^{(\bv)}(\bx)\trans\widehat\bQ_q^{-1}\widehat\bSigma_q\widehat\bQ_q^{-1}\bp^{(\bv)}(\bx),
\end{align*}
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:
\begin{align*}
T(\bx,q) = \frac{\widehat\mu^{(\bv)}(\bx,q)-\mu_0^{(\bv)}(\bx,q)}
{\sqrt{\widehat{\Omega}_{\bv}(\bx,q)/n}}, \qquad (\bx,q)\in\mathcal{X}\times\mathcal{Q},
\end{align*}
where we drop the dependence of $T(\cdot)$ on $\bv$ for simplicity.
Employing Theorem \ref{theorem: bahadur-repres}, 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 \citep{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{theorem: bahadur-repres} 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
\begin{align*}
t(\bx,q) = -\frac{\bp^{(\bv)}(\bx)\trans\bar\bQ_q^{-1}}{\sqrt{\bar\Omega_{\bv}(\bx,q)}}
\mathbb{G}_n\big[\bp(\bx_i)\eta^{(1)}(\mu_0(\bx_i,q)) \psi(y_i,\eta(\mu_0(\bx_i,q));q)\big].
\end{align*}
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 \citep{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 \citep[Theorem B.2]{chen2020jackknife} to generalize prior coupling results \citep[Lemma 36]{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{coro: convergence} and \ref{coro: Mean Square Rate of Convergence}, and Theorem \ref{theorem: bahadur-repres}, 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}).
\begin{thm}[Strong approximation]\label{thm: strong approximation}
Suppose that Assumptions \ref{Assumption: gl-quasi-uniform-partition}--\ref{Assumption: gl-approximation-error} with $\nu\geq3$ hold. Furthermore, assume the following four conditions hold:
\begin{enumerate}[label=\emph{(\roman*)}]
\item $\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}} |\widehat\mu^{(\bv)}(\bx,q)-\mu_0^{(\bv)}(\bx,q)| \lesssim_\P h^{-|\bv|}r_{\tt UC}$.
\item $\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}
\big|\widehat{\mu}^{(\bv)}(\bx,q) - \mu_0^{(\bv)}(\bx,q) - \mathsf{L}^{(\bv)}(\bx,q) \big|
\lesssim_\P h^{-|\bv|}r_{\tt BR}$, with $\frac{\log n}{nh^d} \lesssim r_{\tt BR}$.
\item $\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}
\big|\widehat\Omega_{\bv}(\bx,q)-\bar\Omega_{\bv}(\bx,q)\big|
\lesssim_\P h^{-2|\bv|-d}r_{\tt VC}$, with $r_{\tt VC} = o(1)$.
\item \label{eq: strong approx, lipschitz condition}
$\E\big[\big| \psi(y_i,\eta(\mu_0(\bx_i,q));q) \eta^{(1)}(\mu_{0}(\bx_i,q)) - \psi(y_i, \eta(\mu_{0}(\bx_i,\tilde{q})); \tilde{q}) \eta^{(1)}(\mu_{0}(\bx_i,\tilde{q})) \big|^2 \big | \bx_i \big] \lesssim | q - \tilde{q} |$, for all $q, \tilde{q} \in \mathcal{Q}$.
\end{enumerate}
Then (provided the probability space is rich enough) there exists a stochastic process $(Z(\bx,q):(\bx,q)\in\mathcal{X}\times\mathcal{Q})$ such that, conditional on $\bX_n=(\bx_1,\cdots,\bx_n)$, $Z$ is a mean-zero Gaussian process with $\E[Z(\bx,q)Z(\tilde\bx,\tilde{q})|\bX_n]=\E[t(\bx,q)t(\tilde\bx,\tilde{q})|\bX_n]$ for all $(\bx,q),(\tilde{\bx},\tilde{q})\in\mathcal{X}\times\mathcal{Q}$, and
\begin{align*}
\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}
\big|T(\bx,q)-Z(\bx,q)\big| \lesssim_\P r_{\tt SA} + \sqrt{nh^d}(r_{\tt UC} \; r_{\tt VC} + r_{\tt BR}).
\end{align*}
where $r_{\tt SA}=o(1)$ is any positive sequence satisfying
\[\Big(\frac{1}{nh^{3d}}\Big)^{\frac{1}{10}}\sqrt{\log n} + \frac{\log n}{\sqrt{n^{1-2/\nu}h^{d}}} = o(r_{\tt SA}).\]
Furthermore, if $\overline\psi(\bx_i,y_i)$ is sub-Gaussian conditional on $\bx_i$, then the same result holds with any positive sequence $r_{\tt SA}=o(1)$ satisfying
\[\Big(\frac{1}{nh^{3d}}\Big)^{\frac{1}{10}}\sqrt{\log n} + \frac{(\log n)^{3/2}}{\sqrt{nh^d}} = o(r_{\tt SA}).\]
\end{thm}
The speed of strong approximation in Theorem \ref{thm: strong approximation} 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{coro: convergence} for $r_{\tt UC}$, and Theorem \ref{theorem: bahadur-repres} 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{sec: Feasible Uniform Inference} and \ref{sec: Examples} for more discussion.
With respect to the strong approximation rate, Theorem \ref{thm: strong approximation} 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{thm: strong approximation} 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 \cite{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 \eqref{eq: M estimator}, but we do not pursue this line of research here because it would require a different theoretical treatment.
\subsection{Comparison with Existing Results}\label{sec: Comparison with Existing Results -- SA}
Theorem \ref{thm: strong approximation} 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), \cite{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{thm: strong approximation} 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, \cite{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{thm: strong approximation} 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{theorem: bahadur-repres} only requires $m>d/2$ (not required for the special case of unconnected basis), which is weaker than $m>d$ as implicitly assumed in \cite{Belloni-Chernozhukov-Chetverikov-FernandezVal_2019_JoE}. In practice, this improvement can accommodate, for example, the use of cubic splines for trivariate quantile regression.
Finally, \cite{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 \citep[Theorem 5.1]{Shang-Cheng_2013_AOS}, which leads to substantially slower approximation rates and requires stronger assumption and side rate restrictions; see \cite{Belloni-Chernozhukov-Chetverikov-Kato_2015_JoE} and \cite{Cattaneo-Farrell-Feng_2020_AOS} for more discussion in the context of nonparametric least squares series estimation. In contrast, Theorem \ref{thm: strong approximation}, and our related uniform inference methods, provide a pre-asymptotic approximation with better finite sample properties, faster approximation rates, and weaker regularity conditions.
\section{Feasible Plug-in Uniform Inference}\label{sec: Feasible Uniform Inference}
The (conditional) Gaussian process $(Z(\bx,q):(\bx,q)\in\mathcal{X}\times\mathcal{Q})$ in Theorem \ref{thm: strong approximation} is still infeasible since its covariance structure contains unknowns. This section establishes the validity of a generic \textit{plug-in method} to construct a feasible version of $Z(\bx,q)$. We later employ this result in Section \ref{sec: Examples} 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{thm: strong approximation} is
\begin{align*}
\E[Z(\bx,q)Z(\tilde\bx,\tilde{q})|\bX_n]
= \frac{\bp^{(\bv)}(\bx)\trans\bar{\bQ}_q^{-1} \bar\bSigma_{q,\tilde{q}}\bar{\bQ}_{\tilde{q}}^{-1}\bp^{(\bv)}(\tilde\bx)}
{\sqrt{\bar{\Omega}_{\bv}(\bx,q)\bar{\Omega}_{\bv}(\tilde\bx,\tilde{q})}},
\end{align*}
for all $(\bx,q),(\tilde{\bx},\tilde{q})\in\mathcal{X}\times\mathcal{Q}$, where
\begin{align*}
\bar\bSigma_{q,\tilde{q}}
= \E_n\big[S_{q, \tilde{q}}(\bx_i) \eta^{(1)}(\mu_0(\bx_i, q)) \eta^{(1)}(\mu_0(\bx_i, \tilde{q})) \bp(\bx_i) \bp(\bx_i)\trans]
\end{align*}
with
\begin{align*}
S_{q, \tilde{q}}(\bx)
= \E\big[\psi(y_i, \eta(\mu_0(\bx_i, q)); q) \psi(y_i, \eta(\mu_0(\bx_i, \tilde{q})); \tilde{q}) \big| \bx_i = \bx \big].
\end{align*}
Given context-specific estimates $\widehat{\bQ}_q$ and $\bx \mapsto \widehat{S}_{q, \tilde{q}}(\bx)$, we can put
\begin{align*}
\widehat\bSigma_{q,\tilde{q}}
= \E_n\big[\widehat{S}_{q, \tilde{q}}(\bx_i) \eta^{(1)}(\widehat{\mu}(\bx_i, q)) \eta^{(1)}(\widehat{\mu}(\bx_i, \tilde{q})) \bp(\bx_i) \bp(\bx_i)\trans\big],
\end{align*}
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{sec: Examples} 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
\begin{align}\label{eq: feasible cov}
\E[\widehat{Z}(\bx,q)\widehat{Z}(\tilde{\bx},\tilde{q})|\bD_n]
= \frac{\bp^{(\bv)}(\bx)\trans \widehat\bQ_q^{-1} \widehat\bSigma_{q,\tilde{q}}\widehat\bQ_{\tilde{q}}^{-1}\bp^{(\bv)}(\tilde\bx)}
{\sqrt{\widehat\Omega_{\bv}(\bx,q)\widehat\Omega_{\bv}(\tilde\bx,\tilde{q})}},
\end{align}
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.
\begin{thm}[Feasible Plug-in Strong Approximation]\label{thm: Feasible Plug-in Strong Approximation}
Suppose the assumptions in Theorem \ref{thm: strong approximation} hold. Furthermore, assume the following three conditions hold:
\begin{enumerate}[label=\emph{(\roman*)}]
\item $\|\widehat{\bQ}_{q} - \bar{\bQ}_{q}\|_{\infty} \lesssim_{\P} h^d r_{\mathtt{Q}}$ and $\|\widehat{\bQ}^{-1}_{q}\|_{\infty} \lesssim_{\P} h^{-d}$, with $r_{\mathtt{Q}} = o(1)$.
\item $S_{q,\tilde{q}}(\bx)$ is continuous for all $q, \tilde{q}$, and $\sup_{\bx, q, q_1 \neq q_2} \frac{|S_{q,q_1}(\bx) - S_{q,q_2}(\bx)|}{|q_1 - q_2|} \lesssim 1$.
\item $\sup_{q, \tilde{q}, \bx} |\widehat{S}_{q,\tilde{q}}(\bx) - S_{q,\tilde{q}}(\bx)| \lesssim_{\P} r_{\mathtt{S}}$, with $r_{\mathtt{S}} = o(1)$.
\end{enumerate}
Then (provided the probability space is rich enough) there exists a mean-zero Gaussian process, conditional on $\bD_n$, $(\widehat{Z}(\bx,q):(\bx,q)\in\mathcal{X}\times\mathcal{Q})$ satisfying \eqref{eq: feasible cov} and
\begin{align*}
\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}
\big| \widehat{Z}(\bx,q) - Z^{\star}(\bx,q) \big| \lesssim_\P
[(r_{\mathtt{UC}} + r_{\mathtt{S}})^{1 / 4} + r_{\mathtt{Q}} + r_{\mathtt{VC}}] \sqrt{\log n},
\end{align*}
where $(Z^{\star}(\bx,q):(\bx,q)\in\mathcal{X}\times\mathcal{Q})$ is a process such that, conditional on $\bX_n$, $Z(\cdot)$ and $Z^{\star}(\cdot)$ have the same (conditional) distribution, and $Z^{\star}(\cdot)$ is (conditionally) independent of $(y_1,\ldots,y_n)$.
\end{thm}
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
\begin{align}\label{eq: confidence band}
\mathsf{CB}_{1-\alpha}(\bx,q) = \Big[ \widehat{\mu}^{(\bv)}(\bx,q) \pm \mathfrak{c}_{1-\alpha}\sqrt{\widehat\Omega_{\bv}(\bx,q)/n} \Big]
\end{align}
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 \citep{hall1992effect}, simply ignoring the bias \citep{hall2001bootstrapping}, robust bias correction
\citep{Calonico-Cattaneo-Farrell_2018_JASA,
Calonico-Cattaneo-Farrell_2022_Bernoulli},
or the Lepskii's method \citep{lepskii1992asymptotically,birge2001alternative},
among other possibilities. Thus, under regularity conditions, the confidence band \eqref{eq: confidence band} covers $\mu_0^{(\bv)}$ with probability approximately $1-\alpha$ in large samples, that is,
\begin{align*}
\lim_{n\to\infty} \P \big( \mu_0^{(\bv)}(\bx,q) \in \mathsf{CB}_{1-\alpha}(\bx,q), \;\;\text{for all } (\bx,q)\in\mathcal{X}\times\mathcal{Q} \big) = 1-\alpha.
\end{align*}
\section{Verification of Assumptions in Examples}\label{sec: Examples}
We verify assumptions in our four motivating examples (Section \ref{sec: Motivating Examples}). Assumptions \ref{Assumption: gl-quasi-uniform-partition} and \ref{Assumption: gl-local-basis} concern the partitioning-based methodology itself, Assumption \ref{Assumption: dgp-holder} are already primitive conditions on the data generating process, and Assumption \ref{Assumption: gl-approximation-error} 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{Assumption: dgp-holder}(iv) for details).
Thus, we focus attention on two major issues that remain: (i) how the high-level conditions imposed in Assumptions \ref{Assumption: loss function} and \ref{Assumption: simple setup}, and Condition \emph{\ref{eq: strong approx, lipschitz condition}} in Theorem \ref{thm: strong approximation}, can be verified under intuitive primitive assumptions; and (ii) how to implement uniform inference based on our theory in Section \ref{sec: Feasible Uniform Inference}.
\subsection{Example \ref{example: Generalized Conditional Quantile Regression}: Generalized Conditional Quantile Regression}\label{sec: quantile reg}
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 \eqref{eq: Rho problem}. 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{Assumption: loss function} and \ref{Assumption: simple setup}, and Condition \emph{\ref{eq: strong approx, lipschitz condition}} in Theorem \ref{thm: strong approximation}.
\begin{prop}[Quantile Regression]\label{lem: quantile reg}
Suppose Assumption~\ref{Assumption: dgp-holder} holds with $\mathcal{Q}=[\varepsilon_0, 1-\varepsilon_0]$ for some $\varepsilon_0\in(0,0.5)$, the loss is $\rho(y,\eta;q) = (q-\I(y<\eta))(y-\eta)$, the first moment of $Y$ is finite $\E[|Y|] < \infty$. Assume further that $\eta(\cdot)\colon \mathbb{R} \to \mathcal{E}$ is strictly monotonic and twice continuously differentiable with $\mathcal{E}$ an open connected subset of $\mathbb{R}$ containing the conditional $q$-quantile of $Y|\bX=\bx$, given by $\eta(\mu_{0}(\bx, q))$ for all $(\bx,q)$; $f_{Y|X}(\eta(\mu_{0}(\bx,q)) | \bx)$ is bounded away from zero uniformly over $q \in \mathcal{Q}$ and $\bx \in \mathcal{X}$, and the derivative of $y \mapsto f_{Y|X}(y | \bx)$ is continuous and bounded in absolute value from above uniformly over $y \in \mathcal{Y}_{\bx}$ and $\bx \in \mathcal{X}$.
Then Assumptions \ref{Assumption: loss function}--\ref{Assumption: simple setup} and Condition \emph{\ref{eq: strong approx, lipschitz condition}} in Theorem \ref{thm: strong approximation} hold.
\end{prop}
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 \cite{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{sec: Feasible Uniform Inference}. 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
\begin{align*}
\widehat\bSigma_{q,\tilde{q}} = (q \wedge \tilde{q} - q \tilde{q}) \E_n\big[\eta^{(1)}(\widehat{\mu}(\bx_i, q)) \eta^{(1)}(\widehat{\mu}(\bx_i, \tilde{q})) \bp(\bx_i) \bp(\bx_i)\trans\big].
\end{align*}
On the other hand, the matrix $\bar\bQ_q$
\begin{align*}
\bar\bQ_q = \E_n\big[\bp(\bx_i)\bp(\bx_i)\trans f_{Y|X}(\eta(\mu_0(\bx_i,q))|\bx_i) [\eta^{(1)}(\mu_0(\bx_i,q))]^2\big]
\end{align*}
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 \citep{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{thm: strong approximation} may be used.
\subsection{Example \ref{example: Conditional Distribution Regression}: Generalized Conditional Distribution Regression}\label{sec: distribution reg}
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$.
\begin{prop}[Distribution Regression]\label{lem: distribution reg}
Let $\mathcal{Q} = [-A, A]$ for some $A > 0 $. Suppose that Assumption \ref{Assumption: dgp-holder} holds, the loss is $\rho(y, \eta;q) = \left( \I(y \leq q) - \eta \right)^2$, $\eta(\cdot)\colon \mathbb{R} \to (0,1)$ is strictly monotonic and twice continuously differentiable, $\bx \mapsto F_{Y|X}(q|\bx)$ is a continuous function, and $F_{Y|X}(q|\bx) = \eta(\mu_0(\bx, q))$ lies in a compact subset of $(0,1)$ for all $q\in\mathcal{Q}$ and $\bx\in\mathcal{X}$ (this subset does not depend on $q$ and $\bx$).
Then Assumptions \ref{Assumption: loss function}--\ref{Assumption: simple setup} and Condition \emph{\ref{eq: strong approx, lipschitz condition}} in Theorem \ref{thm: strong approximation} hold.
\end{prop}
The implementation of uniform inference follows the plug-in method described in Section~\ref{sec: Feasible Uniform Inference}. 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
\begin{align*}
\widehat\bSigma_{q,\tilde{q}}
= 4\E_n\big[\bp(\bx_i)\bp(\bx_i)\trans\eta(\widehat{\mu}(\bx_i,q\wedge\tilde{q}))(1-\eta(\widehat\mu(\bx_i,q\vee\tilde{q})))\eta^{(1)}(\widehat\mu(\bx_i,q))\eta^{(1)}(\widehat\mu(\bx_i,\tilde{q}))\big].
\end{align*}
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]$.
\subsection{Example \ref{example: $L_p$ regression}: Generalized \texorpdfstring{$L_p$}{Lp} Regression}\label{sec: lp reg}
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.
\begin{prop}[\texorpdfstring{$L_p$}{Lp} Regression]\label{lem: lp reg}
Suppose that Assumption \ref{Assumption: dgp-holder} holds with the loss function $\rho(y, \eta) = \left| y - \eta \right|^p$, $p \in(1,2]$, and $\eta(\cdot)\colon \mathbb{R} \to \mathcal{E}$ is strictly monotonic and twice continuously differentiable with $\mathcal{E}$ an open connected subset of $\mathbb{R}$. Denoting by $a_l$ and $a_r$ the left and right ends of $\mathcal{E}$ respectively (possibly $\pm \infty$), assume that $\int_{\mathbb{R}} \psi(y;a_l) f_{Y |X}(y|x)dy < 0$ if $a_l$ is finite, and $\int_{\mathbb{R}} \psi(y;a_r)f_{Y |X}(y|x)dy > 0$ if $a_r$ is finite. Also assume that $\E[|Y|^{\nu (p - 1)}]<\infty$ for some $\nu > 2$, and that $\bx \mapsto f_{Y|X}(y | \bx)$ is continuous for any $y \in \mathcal{Y}$. In addition, assume that $\eta \mapsto \int_{\mathbb{R}} |\eta - y |^{p - 1} \sgn(\eta - y) f_{Y|X}(y|\bx)\,\mathrm{d}y$ is twice continuously differentiable with derivatives $\frac{\mathrm{d}^j}{\mathrm{d} \eta^j} \int_{\mathbb{R}} |\eta - y |^{p - 1} \sgn(\eta - y) f_{Y|X}(y|\bx)\,\mathrm{d}y = \int_{\mathbb{R}} |\eta - y |^{p - 1} \sgn(\eta - y) \frac{\partial^j}{\partial y^j} f_{Y|X}(y|\bx)\,\mathrm{d}y$ for $j \in \{1, 2\}$.
Moreover, the function $\int_{\mathbb{R}} |\eta(\zeta) - y |^{p - 1} \sgn(\eta(\zeta) - y) \frac{\partial}{\partial y}f_{Y|X}(y|\bx)\,\mathrm{d}y$ is bounded and bounded away from zero uniformly over $\bx \in \mathcal{X}$ and $\zeta \in B(\bx)$ with $B(\bx)=\{\zeta:|\zeta-\mu_0(\bx)|\leq r\}$ for some $r>0$, and the function $\int_{\mathbb{R}} | \eta(\zeta) - y |^{p - 1} \sgn(\eta(\zeta) - y) \frac{\partial^2}{\partial y^2}f_{Y|X}(y|\bx)\,\mathrm{d}y$ is bounded in absolute value uniformly over $\bx \in \mathcal{X}$ and $\zeta \in B(\bx)$.
Then Assumptions \ref{Assumption: loss function}--\ref{Assumption: simple setup} and Condition \emph{\ref{eq: strong approx, lipschitz condition}} in Theorem \ref{thm: strong approximation} hold.\sloppy
\end{prop}
For implementation, we follow the plug-in method in Section~\ref{sec: Feasible Uniform Inference}. 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.
\subsection{Example \ref{example: Generalized linear models}: Logistic Regression}\label{sec: logit est}
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.
\begin{prop}[Logit Estimation]\label{lem: logit est}
Suppose that Assumption \ref{Assumption: dgp-holder} holds with the loss function $\rho(y,\eta)=-y\log \eta-(1-y)\log (1-\eta)$ and the inverse link $\eta(\theta)=1/(1+e^{-\theta})$; $\mathcal{Y}=\{0,1\}$; $\P(Y = 1 | \bX = \bx)$ is continuous and lies in the interval $(0, 1)$
for all $\bx \in \mathcal{X}$.
Then Assumptions \ref{Assumption: loss function}--\ref{Assumption: simple setup} and Condition \emph{\ref{eq: strong approx, lipschitz condition}} in Theorem \ref{thm: strong approximation} hold.
\end{prop}
It is easy to construct a feasible Gaussian process $\widehat{Z}(\bx)$ conditional on the data $\bD_n$ with covariance structure \eqref{eq: feasible cov}. Standard choices are
\begin{align*}
\widehat\bQ = \E_n[\bp(\bx_i)\bp(\bx_i)\trans\widehat{\eta}_i(1-\widehat{\eta}_i)]
\qquad\text{and}\qquad
\widehat\bSigma = \E_n[\bp(\bx_i)\bp(\bx_i)\trans\widehat\epsilon_i^2],
\end{align*}
where $\widehat{\eta}_i=\eta(\widehat\mu(\bx_i))$ and $\widehat\epsilon_i=y_i-\widehat\eta_i$. See Section \ref{sec: Feasible Uniform Inference} for more discussion.
\section{Other Parameters of Interest}\label{sec: Extensions}
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 \citep[and references therein]{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.
\begin{itemize}
\item \textit{Level Estimator}: $\eta(\widehat\mu(\bx,q))$.
\item \textit{Marginal Effect Estimator}: $\eta^{(1)}(\widehat\mu(\bx,q))\widehat\mu^{(\be_k)}(\bx,q)$.
\item \textit{Conditional Treatment Effect Estimator}: $\eta(\widehat\mu_1(\bx,q))-\eta(\widehat\mu_2(\bx,q))$.
\end{itemize}
Uniform consistency of the three estimators follows from uniform consistency of $\widehat\mu(\bx,q)$ (Corollary \ref{coro: convergence}) because the transformation function $\eta$ is twice continuously differentiable. A Bahadur representation for each of the transformation estimators can be established via Theorem \ref{theorem: bahadur-repres} and a Taylor expansion. For example, for the level estimator,
\begin{align*}
\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}
\big| \eta(\widehat{\mu}(\bx,q))-\eta(\mu_0(\bx,q)) - \mathsf{L}_{\tt LE}(\bx,q) \big|
\lesssim_\P r_{\tt LE}
\end{align*}
with
\begin{align*}
\mathsf{L}_{\tt LE}(\bx,q)
= -\eta^{(1)}(\mu_0(\bx, q)) \bp(\bx)\trans \bQ_{0,q}^{-1}
\E_{n}\big[\bp(\bx_i) \eta^{(1)}(\mu_{0}(\bx_i,q)) \psi(y_i, \eta(\mu_{0}(\bx_i,q));q)\big],
\end{align*}
and for the marginal effect of the $k$th covariate,
\begin{align*}
\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}
\big| \eta^{(1)}(\widehat{\mu}(\bx,q))\widehat{\mu}^{(\be_k)}(\bx,q)
- \eta^{(1)}(\mu_0(\bx,q))\mu_0^{(\be_k)}(\bx,q)
- \mathsf{L}_{\tt ME}(\bx,q) \big|
\lesssim_\P r_{\tt ME}
\end{align*}
with
\begin{align*}
\mathsf{L}_{\tt ME}(\bx,q)
= -\eta^{(1)}(\mu_0(\bx, q))
\bp^{(\be_k)}(\bx)\trans \bQ_{0,q}^{-1}
\E_{n}\big[\bp(\bx_i) \eta^{(1)}(\mu_{0}(\bx_i,q)) \psi(y_i, \eta(\mu_{0}(\bx_i,q));q)\big],
\end{align*}
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
\begin{align*}
\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}
\big| (\eta(\widehat\mu_2(\bx,q))-\eta(\widehat\mu_1(\bx,q)))
- (\eta(\mu_2(\bx,q))-\eta(\mu_1(\bx,q)))
- \mathsf{L}_{\tt CTE}(\bx,q) \big|
\lesssim_\P r_{\tt LE}
\end{align*}
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{sec: uniform inference}. 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))$:
\begin{align*}
\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}
\Big| \frac{\eta(\widehat{\mu}(\bx,q))-\eta(\mu_0(\bx,q))}
{|\eta^{(1)}(\mu_0(\bx,q))|\sqrt{\bar\Omega_{\bm{0}}(\bx,q)/n}}
- Z_{\tt LE}(\bx, q) \Big|
\lesssim_\P r_{\tt SALE}
\end{align*}
with strong approximation rate $r_{\tt SALE}$ as in Theorem \ref{thm: strong approximation}. 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))$:
\begin{align*}
\sup_{q\in\mathcal{Q}}\sup_{\bx\in\mathcal{X}}
\Big| \frac{\eta^{(1)}(\widehat{\mu}(\bx,q)) \widehat{\mu}^{(\be_k)}(\bx,q)
- \eta^{(1)}(\mu_0(\bx,q))\mu_0^{(\be_k)}(\bx, q)}
{|\eta^{(1)}(\mu_0(\bx,q))|\sqrt{\bar\Omega_{\be_k}(\bx,q)/n}}
- Z_{\tt ME}(\bx, q) \Big|
\lesssim_\P r_{\tt SAME}
\end{align*}
with strong approximation rate $r_{\tt SAME}$ as in Theorem \ref{thm: strong approximation}. 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{sec: Feasible Uniform Inference}, and illustrated in Section \ref{sec: Examples}, 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})$.
\section{Conclusion}\label{sec: conclusion}
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.
\section*{Acknowledgments}
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.
\section*{References}
\begingroup
\bibliography{bib}
\bibliographystyle{jasa}
\endgroup
\ifthenelse{\boolean{standalone}}{}{
\clearpage