Adaptive Inference in Multivariate Nonparametric Regression Models Under Monotonicity

Adaptive Inference in Multivariate Nonparametric Regression Models Under Monotonicity

Introduction

We consider the problem of inference on a regression function at a point under the nonparametric regression model

where $f$ is assumed to lie in a H\"older class with exponent $\gamma \in (0,1]$. Procedures based on $\gamma$ is conservative (or suboptimal) when the true regression function in fact lies in a smoother H\"older class with $\gamma' > \gamma$. Adaptive procedures try to overcome this issue by automatically adjusting to the (unknown) underlying smoothness class. However, unlike in the case of estimation, where adaptation to the unknown smoothness class is in general possible with an additional logarithmic term (lepskii1991ProblemAdaptiveEstimationa), adaptation is impossible in the case of inference without further restrictions on the function class (low1997nonparametric) .

Two shape restrictions that can be used to overcome this impossibility have been discussed in the literature, convexity and monotonicity. In this paper, we impose monotonicity on the regression function to construct a CI that adapts to the underlying smoothness of the regression function. The main difference with other papers that consider adaptation under a monotonicity condition (cai2013AdaptiveConfidenceIntervals; armstrong2015AdaptiveTestingRegression) is our general treatment of the dimension of $x_{i}$. To our knowledge, this is the first paper to construct adaptive CIs, under a multivariate nonparametric regression setting.

We consider coordinate-wise monotonicity with respect to all or some of the coordinates. A function $f$ is coordinate-wise monotone with respect to $\mathcal{V} \subseteq \{1, \dots, k\}$ if $x_{j}\geq z_{j}$ for all $j \in\mathcal{V}$ and $x_{j}= z_{j}$ for all $j\notin\mathcal{V}$ imply $f(z) \geq f(z)$. The minimax expected length of a CI over the H\"older class with exponent $\gamma$ converges to $0$ at the well-known rate of $n^{-1/(2 + k/\gamma)}$. When the regression function is monotone in all variables, i.e., $\mathcal{V} = \{1, \dots, k\}$, we can construct a CI that achieves this minimax rate over all $\gamma \in (0,1]$ just as in the univariate case. Also, again as in the univariate case, if the regression is not monotone to any of the variables so that $\mathcal{V} = \emptyset$, there is no scope for adaptation.

An interesting case is when the function is monotone with respect to only some of the variables so that $k_{+} := \lvert \mathcal{V} \rvert < k$, which can arise due to the multivariate nature of the problem. In this case, we show that for a CI that maintains coverage over the H\"older class with exponent $\gamma$, the minimax expected length over a smoother class $\gamma' > \gamma$ converges to $0$ at the rate $n^{-1/(2 + k_{+}/\gamma' + (k-k_{+})/\gamma)}$. The denominator of the exponent can be written as $2 + k/\gamma - k_{+}( 1/\gamma- 1/\gamma')$. This is the sum of a term that comes from the minimax rate over $\gamma$, $2 + k/\gamma$, and $- k_{+}( 1/\gamma- 1/\gamma')$. In this sense, $k_{+}( 1/\gamma- 1/\gamma')$ exactly quantifies the possible gain from monotonicity, indicating larger gains if the regression function is monotone in more variables and/or smoother.

We propose a CI that obtains this minimax rate (of adaptation) for a sequence of H\"older exponents $\{\gamma_{j}\}_{j=1}^{J} \subset (0, 1]$. While the method provided by cai2004adaptation can be used to construct such a CI, we provide an alternative method that builds upon the one-sided CI proposed by armstrong2018optimal. Their one-sided CI “directs power” to a smoother class while maintaining coverage over a larger class of functions. Our CI is constructed by combining the lower and upper versions of their one-sided CI to create a two-sided CI, and then taking the intersection of a sequence of such two-sided CIs that direct power to each $\gamma_{j}$. An appropriate Bonferroni correction is used to obtain correct coverage. This CI can be used in more general nonparametric regression settings, as long as the parameter of interest is a linear functional of the regression function and the regression functions lies in a convex function class.

While the proposed CI obtains the minimax length over $\gamma_{j}$ for each $j$ up to a constant factor that does not depend on the sample size, this constant does depend on the number of parameter spaces $J$ the CI adapts to. This is in contrast with the CI of cai2004adaptation, which gives a multiplicative constant that does not depend on $J$. However, the multiplicative constant of our CI grows slowly with $J$ at a $(\log J)^{1/2}$ rate, and is smaller than the constant given by cai2004adaptation for any reasonable specification of $J$. Even if one wishes to adapt to $J= 10^{3}$ parameter spaces, our CI obtains the minimax expected length of each parameter space within a multiplicative constant of $4.14$, whereas this constant is $16$ for the CI by cai2004adaptation. A simulation study confirms that our CI can be significantly shorter in practice as well. Nonetheless, the uniform constant that cai2004adaptation obtain is theoretically attractive and allows one to adapt to the continuum of H\"older exponents $(0,1]$ in this context.\\

Related literature. An adaptation theory for CIs in a nonparametric regression setting was developed by cai2004adaptation. cai2013AdaptiveConfidenceIntervals provide a procedure for constructing adaptive CIs that adapt to each individual function under monotonicity and convexity. armstrong2015AdaptiveTestingRegression provides an inference method for the regression function at a point, possibly on the boundary of the support, that adapts to the underlying H\"older classes under a monotonicity assumption. As noted earlier, the main difference of our paper is that we consider a multivariate regression setting where there is no restriction on the dimension of the independent variable as long as it is fixed and finite. The adaptation theory for CIs builds upon the more classical minimax theory for CIs, which has been developed in Donoho1994 and low1997nonparametric. cai2012MinimaxAdaptiveInference provides an excellent review on the theory of minimax and adaptive CIs, along with the minimax and adaptive estimation problems.

While the focus of this paper is on adaptive CIs, there are other forms of confidence sets that are of interest in the context of nonparametric regression setting. Adaptive confidence balls have been considered in genovese2005ConfidenceSetsNonparametric, cai2006AdaptiveConfidenceBalls and robins2006AdaptiveNonparametricConfidencea. An adaptation theory for confidence bands has been considered in, for example, dumbgen1998NewGoodnessoffitTests, genovese2008AdaptiveConfidenceBands, and cai2014AdaptiveConfidenceBands. In the context of density estimation, adaptive confidence bands have also been considered in hengartner1995FiniteSampleConfidenceEnvelopes, gine2010ConfidenceBandsDensitya, and hoffmann2011AdaptiveInferenceConfidencea.

Recently, there has been interest in isotonic regression in general dimensions. The monotonicity condition imposed in such models is the same as the one we impose here with $\mathcal{V} = \{1, \dots, k\}$. han2019IsotonicRegressionGeneral derive minimax rates for the least squares estimation problem. deng2020ConfidenceIntervalsMultiple provide a method for constructing CIs at a point based on block max-min and min-max estimators.\\

Outline. Section (ref) describes the nonparametric regression model and the function class we consider. Section (ref) introduces the notion of adaptivity in more detail and describes our procedure for constructing adaptive CIs. Section (ref) presents the main result of the paper, the minimax rate of adaptation, and an adaptive CI that obtains this rate by solving the corresponding modulus problem. Section (ref) provides a simulation study, and Section (ref) illustrates our method in the context of production function estimation.

Any proof omitted in the main text can be found in the appendix. Appendix (ref) collects the proofs for lemmas and corollaries. Appendix (ref) contains the proof for our main theoretical result, Theorem (ref).

Nonparametric Regression Under Monotonicity

We observe $\left\{ \left(y_{i},x_{i}\right)\right\} _{i=1}^{n}$ and consider a nonparametric regression model,

where $x_{i}\mathcal{\in X}\subset\mathbb{R}^{k}$ is a (fixed) regressor, $f:\mathbb{R}^{k}\to\mathbb{R}$ is the unknown regression function that lies in some function class $\mathcal{F}$, and $u_{i}$'s are independent with $u_{i}\sim N(0,\sigma^{2}(x_{i}))$ and $\sigma^{2}(\cdot)$ known. The parameter of interest is $f(x_{0})$. For the rate results provided in Section (ref), we require that $x_{0} \in \mathrm{Int}\, \mathcal{X}$. However, we note that the solution to the modulus problem given in Section (ref) does not depend on whether $x_0$ is on the boundary or not. Without loss of generality, we normalize $x_{0}$ to be $0$.

We take the $\mathcal{F}$ to be the class of functions that are H\"older continuous and nondecreasing in all or some of the variables. Let $\Lambda(\gamma,C)$ denote the set of functions from $\mathbb{R}^{k}$ to $\mathbb{R}$ that are H\"older continuous with H\"older constants $(\gamma,C)$,

where $\mathcal{F}\left(\mathbb{R}^{k}, \mathbb{R}\right)$ is the set of functions from $\mathbb{R}^{k}$ to $\mathbb{R}$, $\gamma\in[0,1],$ $C \geq 0$ and $\left\Vert \cdot\right\Vert $ is a norm on $\mathbb{R}^{k}$. For notational simplicity, we omit the dependence of the function class on the choice of the norm $\lVert \cdot \rVert$. We impose the following restriction that $\lVert \cdot \rVert$ is monotone in the magnitude of each element, which is satisfied by most norms used in practice. such as the $\ell_{p}$ norm or a weighted version of it. We discuss the relationship between this assumption and the monotonicity of the regression function in Remark (ref).

We now define the (coordinate-wise) monotone H\"older class. For a subset of the covariate indices $\mathcal{V}\subset\left\{ 1,\dots,k\right\} $, write

This is the set of H\"older continuous functions that are nondecreasing, coordinate-wise, with respect to the $j$th element for $j \in \mathcal{V}$. Define $k_{+}:=\left\lvert \mathcal{V}\right\rvert .$ By a relabeling argument, it is without loss of generality to write $\mathcal{V}:=\left\{ 1,\dots,k_{+}\right\} $. If $k_{+}=k$, then $\Lambda_{+,\mathcal{V}}(\gamma,C)$ is the set of nondecreasing and H\"older continuous functions where the monotonicity is with respect to the coordinate-wise partial ordering on $\mathbb{R}^{k}$.

Adaptive Confidence Intervals

Notion of Adaptivity

In this section, we discuss the problem of inference for a general linear functional of the regression function, $Lf$. Consider a sequence of convex parameter spaces $\mathcal{F}_{1}$, \dots, $\mathcal{F}_{J}$, with the requirement that $\mathcal{F}_{j} \subset \mathcal{F}_{J}$ for all $j \leq J$. Note that the parameter spaces are not necessarily nested, but there is a largest convex parameter space that nests all the other parameter spaces. Here, $\mathcal{F}_{J}$ reflects a conservative choice of the parameter space where the researcher believes the true regression function to lie in. Hence, the CI we construct will be required to maintain correct coverage over this space. An adaptive CI maintains this correct coverage over the largest parameter space $\mathcal{F}_{J}$ while having good performance (e.g. shorter expected length) when the true function happens to lie in the smaller parameter space $\mathcal{F}_{j}$, simultaneously for all $j \leq J$.

Then, a natural question is how well a CI that maintains coverage over $\mathcal{F}_{J}$ can perform over $\mathcal{F}_{j}$, which is one of the main questions that cai2004adaptation raise and address in detail in the context of two-sided CIs. The case of one-sided CIs has been considered by armstrong2018optimal, along with other questions.

Two-sided Adaptive CIs

Let $\mathcal{I}_{\alpha,2}^{J}$ denote the set of all two-sided CIs that have coverage at least $1-\alpha$ over $\mathcal{F}_{J}.$ Following cai2004adaptation, the performance criterion we consider for two-sided CIs is the worst-case expected length. That is, the performance of a CI, $CI$, over the parameter space $\mathcal{F}_{j}$ is measured by $\sup_{f\in\mathcal{F}_{j}}\operatorname{\mathbf{E}}_{f}\mu(CI)$ with smaller values of this quantity meaning better performance. Here, $\operatorname{\mathbf{E}}_{f}$ denotes the expectation when the true regression function is $f$ and $\mu$ is the Lebesgue measure on the real line. Then, the shortest possible worst-case expected length a CI can achieve over $\mathcal{F}_{j}$ (while maintaining correct coverage over $\mathcal{F}_{J}$) is characterized by the quantity

Following cai2004adaptation, we say a CI is adaptive if it achieves $L_{j,J}^{\ast}$ for all $j \leq J$ up to a multiplicative constant that does not depend on the sample size. Let $z_{q}$ denote the $q$--quantile of the standard normal distribution. cai2004adaptation show that $L_{j,J}^{\ast} \asymp\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J}),$ with $\asymp$ denoting asymptotic equivalence\footnote{We write $a_{n}\asymp b_{n}$ if $ 0<\underset{n\to\infty}{\lim\inf}({a_{n}/}{b_{n}})\leq\underset{n\to\infty}{\lim\sup}({a_{n}}/{b_{n}})<\infty.$} and $\omega_{+}(\delta,\mathcal{F}_{j},\mathcal{F}_{J})$ is the between class modulus of continuity defined as

for $\delta \geq 0.$\footnote{Note that the definition is slightly different with cai2004adaptation due to the $\sigma(x_{i})$ term that appears in the denominator of the summand. This is because we divide both sides of (ref) by the (known) $\sigma(x_{i})$ to convert the model into the same form as that of cai2004adaptation.} In general, $\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J})$ is more tractable than $L_{j,J}^{\ast}$, and thus the strategy is to construct a CI that has worst case length over $\mathcal{F}_{j}$ bounded by $\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J})$, up to a multiplicative constant. We refer to the rate at which $\omega_{+}(z_{1-\alpha}, \mathcal{F}_{j}, \mathcal{F}_{J}) $ converges to $0$ as the minimax rate of adaptation (of $\mathcal{F}_{j}$ over $\mathcal{F}_{J}$). If $\mathcal{F}_{j}=\mathcal{F}_{J}$, this is the minimax rate over $\mathcal{F}_{J}$, which is the fastest rate at which the worst-case expected length over $\mathcal{F}_{J}$ of a CI that maintains correct coverage over the same space $\mathcal{F}_{J}$ can achieve.

One-sided CIs

While our main focus is on adaptive two-sided CIs, the construction of our adaptive CI relies heavily on the one-sided CI proposed by armstrong2018optimal. Hence, we briefly describe the notion of adaptivity in the context of one-sided CIs. For one-sided CIs, we follow armstrong2018optimal and consider the $\beta$th quantile of excess length as the performance criterion. More specifically, for a one sided lower CI, $[\hat{c},\infty)$, we denote the $\beta$th quantile of the excess length at $f$ as $q_{\beta,f}(Lf-\hat{c})$, where $q_{\beta,f}(\cdot)$ denotes the $\beta$th quantile function when the true regression function is $f$. Under this criterion, the best possible performance over $\mathcal{F}_{j}$ is quantified by

where $\mathcal{I}_{\alpha, \ell}^{J}$ denotes the set of all one-sided lower CIs that have coverage at least $1-\alpha$ over $\mathcal{F}_{J}.$ armstrong2018optimal showed that $\ell_{j,J}^{\ast}=\omega(z_{1-\alpha}+z_{\beta},\mathcal{\mathcal{F}}_{J},\mathcal{\mathcal{F}}_{j})$, where $\omega(z_{1-\alpha}+z_{\beta},\mathcal{\mathcal{F}}_{J},\mathcal{\mathcal{F}}_{j})$ is the ordered class modulus of continuity defined as

for any $\delta \geq 0$ and $j,k \leq J$. We refer to the optimization problem in the definition as the ordered modulus problem. Naturally, an analogous result holds for upper one-sided CIs so that $u_{j,J}^{\ast}=\omega(z_{1-\alpha}+z_{\beta},\mathcal{\mathcal{F}}_{j},\mathcal{\mathcal{F}}_{J})$, where

with $\mathcal{I}_{\alpha, u}^{J}$ denoting the set of all one-sided upper CIs that have coverage at least $1-\alpha$ over $\mathcal{F}_{J}.$

We say a one-sided lower CI, $[\hat{c}^{\ast},\infty)$, is adaptive if there exists some $c>0$ that does not depend on $n$ such that

for all $j \leq J$, and similarly for one-sided upper CIs.

Modes of Adaptation

Note that it must be the case that $\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{j}) \leq \omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J})$ (and similarly for the ordered moduli) because $ \omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{j})$ takes the supremum over a smaller set. However, if it happens to be the case that $\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{j})\asymp\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J})$, an adaptive CI, $CI^{\ast}$, satisfies $\sup_{f\in\mathcal{F}_{j}}\operatorname{\mathbf{E}}_{f}\mu(CI^{\ast})\leq\overline{c}\,L_{j,j}^{\ast}$ for all $j\leq J$. cai2004adaptation define such CI to be strongly adaptive. This is an ideal case because we obtain $L_{j,j}^{\ast}$, up to a multiplicative constant, which is the minimax length we could have achieved if we “knew” that our true regression function lied in the smaller class $\mathcal{F}_{j}$ (i.e., if we made a stronger assumption that the true regression function lies in this smaller class). While adaptive CIs exist in general, strong adaptation is possible only when $\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J})\asymp\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{j})$ for all $j \leq J$. This is not a property of a given procedure, but of the given statistical model.

The least desirable case is when $\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J})\asymp\omega_{+}(z_{1-\alpha},\mathcal{F}_{J},\mathcal{F}_{J})$, because this leaves no scope of adaptation. An intermediate case is when

so that the minimax rate of adaptation is better than the worst-case minimax rate over $\mathcal{F}_{J}$ but not as good as the minimax rate over $\mathcal{F}_{j}$.\footnote{For positive sequences $\left\{ a_{n}\right\} $ and $\left\{ b_{n}\right\} $, we write $a_{n} \prec b_{n}$ if $\underset{n\to\infty}{\lim\inf}({b_{n}/}{a_{n}}) = \infty$.} That is, one can do better than simply taking the most conservative parameter space as the true space but not quite as good as knowing that the true function actually lies in the smaller parameter space. Hence, the minimax adaptation rate plays an important role in determining whether sharp adaptation is possible. In Section (ref), we derive the minimax rates of adaptation under the model given in Section (ref).

Construction of Adaptive CIs

cai2004adaptation provide a general method of constructing adaptive CIs of $Lf$ under the general model (ref). Here, we provide an alternative method that is intuitive and gives smaller constants in the case of non-nested parameter spaces.\footnote{For a given adaptive CI, $CI^{\ast}$, we refer to the positive number $c$ (that does not depend on $n$) such that $ \sup_{f\in\mathcal{F}_{j}}E\mu(CI^{\ast})\leq c\,\omega_{+}\left(z_{\alpha},\mathcal{F}_{j},\mathcal{F}_{J}\right), $ as the “constant” of $CI^{\ast}$.} For the nested case, the CI of cai2004adaptation has a bounded constant even as $J \to \infty$, which is an attractive theoretical property. For the CI we propose, the constant will grow with $J$ in general. In practice, however, one can only adapt to finitely many parameter spaces due to computational constraints. The proposed procedure gives a smaller constant than that of cai2004adaptation even for unrealistically large values of $J$ (e.g., $J = 10^{10}$).

The main building block for our adaptive CI is the minimax one-sided CI proposed by armstrong2018optimal, which relies on the ordered modulus. We say that $(f_{j}, f_{k}) \in \mathcal{F}_{j} \times \mathcal{F}_{k} $ is a solution to $\omega(\delta, \mathcal{F}_{j}, \mathcal{F}_{k})$ if $(f_{j}, f_{k})$ solves the optimization problem corresponding to $\omega(\delta, \mathcal{F}_{j}, \mathcal{F}_{k})$. Let $(f_{J,{\delta}}^{*,Jj}, g_{j,{\delta}}^{*,Jj})\in\mathcal{F}_{J}\times\mathcal{F}_{j}$ be a solution to the ordered modulus $\omega\left({\delta},\mathcal{F}_{J},\mathcal{F}_{j}\right),$ and define the estimator

where $\omega'(\cdot, \mathcal{F}_{J}, \mathcal{F}_{j})$ is the derivative of $\omega(\cdot, \mathcal{F}_{J}, \mathcal{F}_{j})$. Based on this estimator, define a lower one-sided CI by subtracting the maximum bias and an appropriately scaled normal quantile:

The following theorem from armstrong2018optimal shows that for a specific choice of $\delta$, this CI is optimal in the sense that it achieves $\ell^{\ast}_{j,J}.$

The excess length $Lf-\hat{c}_{\alpha, \underline{\delta}}^{\ell,j}$ follows a Gaussian distribution because it is a affine transformation of the data, which follows a Gaussian distribution by assumption. Hence, the median and mean of the excess length are the same. Taking $\beta=1/2$, we can replace $q_{f,\beta}$ with the expectation under $f$, which gives

where we define $\hat{c}_{\alpha}^{\ell,j}:= \hat{c}_{\alpha, z_{1-\alpha}}^{\ell,j}.$ Likewise, we can define an optimal upper one-sided CI $(-\infty, \hat{c}_{\alpha, \underline{\delta}}^{\ell,j}]$ such that

where the precise definition of $\hat{c}_{\alpha, \underline{\delta}}^{\ell,j}$ is given in Appendix (ref). Similarly, let $\hat{c}_{\alpha}^{u,j}$ denote the upper counterpart of $\hat{c}_{\alpha}^{\ell,j}.$

Using the optimal one-sided CIs, we first show how a naive Bonferroni procedure leads to a two-sided adaptive CI. We then provide a method that improves upon this naive Bonferroni CI by taking into account the correlation among the CIs. The naive Bonferroni CI is defined as

This has coverage at least $1-\alpha$ over $\mathcal{F}_{J}$ because each $[\hat{c}_{\alpha/2J}^{\ell,j},\hat{c}_{\alpha/2J}^{u,j}]$ has coverage $1 - \alpha/J$ over $\mathcal{F}_{J}$ and $CI_{\alpha}^{Bon,J}$ is simply the intersection of such CIs. The following theorem shows that this CI is indeed adaptive.

The constant $2z_{1-\frac{\alpha}{2J}}/z_{1-\frac{\alpha}{2}}$ increases with the number of parameter spaces $J$.\footnote{The constant, $z_{1-\frac{\alpha}{2J}}/z_{1-\frac{\alpha}{2}}$, grows with $J$ at the rate $(\log J)^{1/2}$. This is the same rate that cai2004adaptation find in their analysis of the case with non-nested parameter spaces. Their constant is at least eight times greater than what we provide here, but does not require that the largest space in consideration is convex.} On the other hand, the constant given in cai2004adaptation is 16 and thus does not depend on the number of parameter spaces. However, we note that $2z_{1-\frac{\alpha}{2J}}/z_{1-\frac{\alpha}{2}}$ is not too large, in fact smaller than $16$, for reasonable specifications of $J.$ For example, when $\alpha=0.05$ and $J=50$, we get $2z_{1-\frac{\alpha}{2J}}/z_{1-\frac{\alpha}{2}}\approx3.36$, which is considerably smaller than the constant given in cai2004adaptation. Even for unrealistically large $J$ such as $J=10^{10}$, we have $2z_{1-\frac{\alpha}{2J}}/z_{1-\frac{\alpha}{2}}<8$, which is still less than half of the constant given by cai2004adaptation. Simulation results given in Section (ref) confirm that not only the upper bound, but also the actual length itself is often much shorter for our CI.

The naive CI given in (ref) does not take into account the possible correlation among the CIs that we take the intersection of. However, if parameter spaces are “close” to each other, the corresponding CIs will be correlated, implying that there is room for improvement over the Bonferroni procedure. Consider the CIs of the form $CI^{\tau,\mathcal{J}}=\cap_{j=1}^{J}[\hat{c}_{\tau}^{\ell,j},\hat{c}_{\tau}^{u,j}]$. If we take $\tau = \alpha/(2J)$, this is precisely the CI given in (ref). The CI that gives the smallest constant among CIs of such forms is $CI^{\tau^{\ast},\mathcal{J}},$ where $\tau^{\ast}$ is the largest possible $\tau$ such that $CI^{\tau,\mathcal{J}}$ has correct coverage over $\mathcal{F}_{J}$:

We know that $\tau = \alpha/(2J)$ satisfies the constraint, and also that any $ \tau > \alpha$ does not because then $[\hat{c}_{\tau}^{\ell,j}, \infty)$ will have coverage probability $1-\tau < 1 - \alpha$. Hence, we can restrict $\tau$ to lie in $[\alpha/(2J), \alpha]$.

However, the coverage probability $ \inf_{f \in \mathcal{F}_{J}}\operatorname{\mathbf{P}} _{f}(Lf \in CI^{\tau,\mathcal{J}})$ is unknown in general, rendering $CI^{\tau^{\ast},\mathcal{J}}$ infeasible. Instead, we replace this coverage probability with a lower bound that we can calculate either analytically or via simulation. Then, we take $\tau^{\ast}$ as the largest value that makes this lower bound at least $1-\alpha$. As we show later, using $\tau^{\ast}$ rather than $\alpha/(2J)$ can only make the resulting CI shorter.

Let $(V(\tau)', W(\tau)')'$ be a centered Gaussian random vector with unit variance. The covariance terms for $V(\tau)=\left(V_{1}(\tau),...,V_{J}(\tau)\right)' $ is given by \[ \text{Cov}\left(V_{j}(\tau),V_{\ell}(\tau)\right)=\frac{1}{z_{1-\tau}^{2}}\sum_{i=1}^{n}\big(g_{j,z_{1-\tau}}^{*, Jj}(x_{i})-f_{J,z_{1-\tau}}^{*, J j}(x_{i})\big)\big(g_{\ell,z_{1-\tau}}^{*, J\ell}(x_{i})-f_{J,z_{1-\tau}}^{*, J \ell}(x_{i})\big). \] Likewise, the covariance terms for $ W(\tau)=\big(W_{1}(\tau),...,W_{J}(\tau)\big)'$ is given by \[ \text{Cov}\big(W_{j}(\tau),W_{\ell}(\tau)\big)=\frac{1}{z_{1-\tau}^{2}}\sum_{i=1}^{n}\big(g_{j,z_{1-\tau}}^{*, jJ}(x_{i})-f_{J,z_{1-\tau}}^{*, jJ}(x_{i})\big)\big(g_{\ell,z_{1-\tau}}^{*, \ell J}(x_{i})-f_{J,z_{1-\tau}}^{*, \ell J}(x_{i})\big). \] Finally, the covariance terms across $V(\tau)$ are $W(\tau)$ given as \[ \text{Cov}\big(V_{j}(\tau),W_{\ell}(\tau)\big)=\frac{1}{z_{1-\tau}^{2}}\sum_{i=1}^{n}\big(g_{j,z_{1-\tau}}^{*, Jj}(x_{i})-f_{J,z_{1-\tau}}^{*, Jj}(x_{i})\big)\big(g_{\ell,z_{1-\tau}}^{*, \ell J}(x_{i})-f_{J,z_{1-\tau}}^{*, \ell J}(x_{i})\big). \] This Gaussian random vector can be used to tune the critical value, as the following lemma implies.

Such a $\tau^{\ast}$ always exists because the inequality (ref) holds with $\tau = \alpha/(2J)$ due to the union bound. A solution $\tau^{*}$ can be found via numerical simulation. By construction, its length will be also bounded by ((ref)). In Section (ref), we show that as $n \to \infty$ the distribution of $(V(\tau)', W(\tau)')'$ does not depend on $\tau$, under our setting of $Lf = f(0)$ with $f$ belonging to a H\"older class. Hence, finding $\tau^{*}$ boils down to simply finding the $1-\alpha$ quantile of the maximum of a Gaussian vector in this case.

Adaptive Inference for $f(0)$

In this section, we provide an adaptive inference procedure for $f(0)$. To construct the adaptive CI introduced in Section (ref), we first solve the corresponding modulus problem. By using this solution to the modulus problem, we derive the minimax rate of adaptation. Finally, we provide a CI that obtain this rate, using the method described in Section (ref).

Solution to the Modulus Problem

Let $\Lambda_{+,\mathcal{V}}(\gamma_{j},C_{j})\subset\Lambda_{+,\mathcal{V}}(\gamma_{J},C_{J})$ with $\gamma_{j}\geq\gamma_{J}$ and $C_{j}\leq C_{J}$. To construct the adaptive CI, we first calculate the ordered moduli, $\omega\left(\delta,\Lambda_{+,\mathcal{V}}\left(\gamma_{j},C_{j}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{J},C_{J}\right)\right)$ and $\omega\left(\delta,\Lambda_{+,\mathcal{V}}\left(\gamma_{J},C_{J}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{j},C_{j}\right)\right),$ for each $j=1,\dots,J.$ For notational simplicity, we consider the case with $J=2$ and solve $\omega_{+}\left(\delta,\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right)\right)$, from which the general solution follows immediately.

Recall the definition of the ordered modulus of continuity

with the maximized value denoted by $\omega\left(\delta,\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right)\right).$ It is convenient to solve the inverse modulus problem instead, which is defined as

for $b>0$, with the square root of the maximized value denoted by the inverse (ordered) modulus $ \omega^{-1}\left(b,\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right)\right).$ We provide a closed form solution for the this inverse problem, from which we can recover the solution to the original problem by finding $b$ such that $\omega^{-1}\left(b,\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right)\right)=\delta$. Note that this is simply a search problem on the positive real line.

To characterize the solution to (ref), we show two simple lemmas about the properties of the class $\Lambda_{+,\mathcal{V}}\left(\gamma,C\right)$. For $z=(z_{1},\dots,z_{k})\in\mathbb{R}^{k}$, define

\[ \left(z\right)_{\mathcal{V}+}=

\] and $\left(z\right)_{\mathcal{V}-}=\left(-z\right)_{\mathcal{V}+}.$

The following lemma asserts that the class of functions we consider is closed under the maximum operator.

The next lemma can be used to establish the solutions to the problem ((ref)). This is a generalization of Proposition 4.1 of beliakov2005monotonicity, which gives the same result for the special case of $\gamma=1$.

We are now ready to characterize the solution to the inverse modulus problem ((ref)). For $r\in\mathbb{R}$, define $\left(r\right)_{+}:= \max\left\{ r,0\right\} $.

The following corollary states an analogous result regarding the inverse modulus $\omega^{-1}\left(b,\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right)\right)$.

Minimax Rate of Adaptation

Using this solution to the inverse modulus, we derive the rate of convergence of the between class of modulus, which characterizes how fast the worst-case expected length of the adaptive CIs can go to 0 as $n\rightarrow\infty$. We derive the rates under the assumption that the sequence of design points $\{x_{i}\}_{i=1}^{\infty}$ is a realization of a sequence of independent and identically distributed random vectors $\{ X_{i} \}_{i=1}^{\infty}$ drawn from a distribution that satisfies some mild regularity conditions. This gives an intuitive restriction on the design points, and also shows that the result applies under random design points as well.\footnote{Consider the model $y_{i}=f(X_{i})+\varepsilon_{i}$, for $i=1,\dots,n$, with the $X_{i}\overset{\mathrm{i.i.d.}}{\sim}p_{X}$ with $\varepsilon_{i}|X_{i}\sim N(0,\sigma^{2}(X_{i}))$. Then, conditional on $\{ X_{i} \}_{i=1}^{n}=\{x_{i}\}_{i=1}^{n}$, this model is equivalent with our model.} Define $r(\gamma_{1}, \gamma_{2}) = ({2+k_{+}/\gamma_{1}+(k-k_{+})/\gamma_{2}})^{-1}$. The following theorem fully characterizes the minimax rate of adaptation.

Theorem (ref) shows how the monotonicity restriction plays a role in determining the minimax rates of adaptation to H\"older coefficients under the multivariate nonparametric regression setting. When $k_{+}=k$, the minimax rate of adaptation is $n^{-\frac{1}{2+k/\gamma_{1}}}$, which equals the minimax convergence rate over $\omega\left(\delta,\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right)\right)$. This shows that strong adaptation is possible if the regression function is monotone with respect to all the variables, just like in the univariate case. On the other hand, when $k_{+}=0$, the rate becomes $n^{-\frac{1}{2+k/\gamma_{2}}}$, consistent with the previous findings that there is no scope of adaptation for general H\"older classes without any shape constraint. Importantly, Theorem (ref) characterizes the convergence rate for the case where $0<k_{+}<k$, where it gives an intuitive intermediate rate between the two extreme.

Construction of the Adaptive CI

Here, we give the explicit formula of the CIs for our parameters of interest, now that we have derived the form of the moduli of continuity and the solutions to the modulus problems in the previous section. We first consider $L_{0}f$. Before stating the result, it is convenient to define the following functions

The first terms in the formula of $\hat{L}_{\delta}^{\ell,j}$ and $\hat{L}_{\delta}^{u,j}$ are the random terms linear in $y_{i}$ while the remaining terms are non-random fixed terms. If $\mathcal{V}=\{1,...,k\}$ (so the function is monotone in every coordinate), the random terms can be viewed as a kernel estimator with a data-dependent bandwidth. Too see this, if we define \[ k(x)=\left[1-C_{j}\left\lVert \left(x\right)_{\mathcal{V}-}\right\rVert ^{\gamma_{j}}-C_{J}\left\lVert \left(x\right)_{\mathcal{V}+}\right\rVert ^{\gamma_{J}}\right]_{+}, \] and \[ h_{mn}\left(x\right)=

\] we have \[ \frac{\sum_{i=1}^{n}D_{Jj,\delta}\left(x_{i}\right)y_{i}}{\sum_{i=1}^{n}D_{Jj,\delta}\left(x_{i}\right)}=\frac{\sum_{i=1}^{n}k\left(x_{1i}/h_{1n}\left(x_{i}\right),...,x_{ki}/h_{kn}\left(x_{i}\right)\right)y_{i}}{\sum_{i=1}^{n}k\left(x_{1i}/h_{1n}\left(x_{i}\right),...,x_{ki}/h_{kn}\left(x_{i}\right)\right)}. \] Hence, the CI can be considered to be based on a Nadaraya-Watson type estimator, correcting for the bias.

As described in Section (ref), the proposed CI is given by $ \cap_{j=1}^{J} [\hat{c}_{z_{1-\tau^{\ast}}}^{\ell,j}, \hat{c}_{z_{1-\tau^{\ast}}}^{\ell,j}]$, where $\tau^{\ast}$ is defined in Lemma (ref). Here, we show that the distribution of $(V(\tau)',W(\tau)')'$ does not depend on $\tau$ as $n \to \infty$. The implication of this invariance with respect to $\tau$, is that calculating $\tau^{\ast}$ boils down to calculating the quantile of the maximum of Gaussian vectors. The variance matrix of this limiting Gaussian random vector is known, and thus the said quantile can be easily simulated. Moreover, when $\gamma_{1} = \cdots = \gamma_{J}$ so that the parameters spaces differs only in $C_{j}$, $\tau^{\ast}$ can be shown to be bounded away from zero by a constant that does not depend on $J$, for large $n$. Hence, the constant of the CI does not grow to infinity as $J \to \infty$ in this case.\footnote{This is especially useful when one wishes to adapt to $C$ while keeping $\gamma$ fixed. For example, kwon2020InferenceRegressionDiscontinuity take $\gamma_{j} = 1$ and consider the problem of adapting to the Lipschitz constant in a regression discontinuity setting.}

Simulation Results

In this section, we compare the performances of the adaptive CI of cai2004adaptation and the adaptive CI constructed using the naive Bonferroni procedure described in Section (ref). As a benchmark, we also provide the lengths of the shortest fixed length confidence intervals of Donoho1994, referred to as minimax CIs. We consider inference for $f(0)$, given some regression function $f$. We consider the case where the researcher is uncertain about the value of the H\"older exponent $\gamma$, and thus tries to adapt to its value.

First, we construct adaptive CIs with respect to two smoothness parameters $\left(\gamma_{1},\gamma_{2}\right)=\left(1,10^{-3}\right)$ while fixing $C = 1$, which gives $J = 2$. We vary $n$ over $\{10^2,\ 5 \times 10^2,\ 10^3,\ 5 \times 10^3,\ 10^4 \}$ to investigate the rate of adaptation as the sample size grows. The true regression function is over $\mathbb{R}^2$ and given by either $f_{1}$ or $f_{2}$, defined as

By construction, we have $f_{j}\in\Lambda_{+,\mathcal{V}}(\gamma_{j},1)$. The covariates are drawn from a uniform distribution over $[-1/(2\sqrt{2}), 1/(2\sqrt{2})]^{2}$, and the noise terms, $\{u_{i}\}_{i=1}^{n}$, are drawn from a standard normal distribution. The outcome variable is given as $y_{i} = f(x_{i}) + u_{i}$, for $f \in \{f_1, f_2\}$. We fix the draw of $\left\{ x_{i} \right\}_{i=1}^{n}$ within each simulation iteration. We run 500 iterations to calculate the average lengths and coverage probabilities of CIs. The nominal coverage probability is $.95$ for all CIs.

Table (ref) shows the results for the case where $f=f_{1}$. Each column corresponds to 1) our proposed (naive) Bonferroni adaptive procedure (AdaptBonf), 2) the adaptive CI of cai2004adaptation (CL, henceforth), 3) the minimax CI with respect to $\Lambda_{+,\mathcal{V}}(\gamma_2, 1)$, and 4) the minimax CI with respect to $\Lambda_{+,\mathcal{V}}(\gamma_1, 1)$. Regarding the last two minimax procedures, we refer to them as the “conservative minimax CI” and the “oracle minimax CI”, respectively. Note that the oracle minimax CI is an optimal benchmark, which is only feasible when we actually know the true regression function is in the smaller parameter space $\Lambda_{+,\mathcal{V}}(\gamma_1, 1)$. In Table (ref), the average lengths of both adaptive confidence intervals decrease considerably as $n$ increases from 100 to 10,000. In comparison, the length of the conservative minimax CI (column 3) decreases only about 28% for the same change in the sample size. This shows the lengths of the adaptive confidence intervals decrease more sharply when the true function is smooth, as predicted by the theory.

To compare the performances of different adaptive inference procedures, note that the average lengths of the CI of CL adapting to the H\"older exponents (column 2) are often wider than the conservative minimax CI (column 3). When $n = 100$, the former is more than three times wider than the latter, and the adaptive procedure starts to dominate the minimax procedure only when $n$ is greater than 5,000. In comparison, our proposed Bonferroni adaptive procedure (column 1) yields shorter CIs than those by CL, as predicted in Section (ref). To compare the Bonferroni adaptive CI with the conservative minimax CI, the lengths of the former are always exceeded by those of the minimax CI, even for the relatively small sample size of $n=100$. Moreover, the length of the adaptive CI becomes only 20% of the length of the conservative minimax CI for the sample size of $n=10^4$. The Bonferroni procedure also performs well even when compared to the infeasible oracle minimax CI (column 4), with the length of the former only 13% wider than the latter when $n = 10^4$. This demonstrates the strong adaptivity property of the adaptive procedure when the regression function is monotone with respect to all variables, as shown in Section (ref).

Table (ref) demonstrates the analogous simulation results when $f=f_{2}$. In this case, the minimax CI with respect to $\Lambda_{+,\mathcal{V}}(\gamma_2, 1)$ (column 3) is referred to as the oracle minimax CI. While the lengths of the oracle minimax procedure are considerably shorter than the CIs of CL for various values of $n$, the performance of the Bonferroni CIs almost matches that of the oracle minimax procedure. Especially, the performance of the Bonferroni adaptive procedure becomes extremely close to the oracle minimax procedure when $n$ is greater than 500.

Table (ref) shows the coverage probabilities of adaptive CIs for both of the cases when $f = f_1$ and $f = f_2$. While all the CIs achieve the correct coverage, none of those CIs exactly achieves the nominal coverage of $.95$, reflecting the conservative nature of the adaptive CIs. We can see that the adaptive procedure of CL is particularly conservative, almost always yielding 100% coverage probabilities.

So far we considered adapting to the smoothness parameters at two extremes, $\gamma\in(0.001,1)$. Since the multiplicative constant for the Bonferroni procedure increases with $J$, a concern is that the performance of the Bonferroni procedure relative to the CL procedure might get worse when $J$ is larger. To investigate the possibility, we consider adapting to a wider set of parameters, $\{\gamma_j\}_{j = 1}^6$, where $\gamma_j = 1 - (j - 1)/5$ for $j = 1,...,5$ and $\gamma_6 = 10^{-3}$. Moreover, rather than taking the extreme value of $\gamma$ as the true parameter, we consider the case where $\gamma$ takes an intermediate value, $\gamma=1/2$. The true regression function is given by \[ f_3(x_{1},x_{2})=\left\lVert \left(x_{1},x_{2}\right)_{\mathcal{V}+}\right\rVert _{2}^{1/2}, \quad \mathcal{V} = \{1,2\}, \] so that $f_3 \in\Lambda_{+,\mathcal{V}}(1/2,1)$.

Table (ref) displays the simulation results corresponding to this specification. Each column corresponds to 1) our proposed Bonferroni adaptive procedure, 2) the adaptive CI of CL, 3) the minimax CI with respect to $\Lambda_{+,\mathcal{V}}(\gamma_6, 1)$, and 4) the minimax CI with respect to $\Lambda_{+,\mathcal{V}}(1/2, 1)$. As before, we refer to the last two CIs as the conservative minimax CI and the oracle minimax CI, respectively. We observe the same pattern as in the case of adapting to two parameters---adaptive CIs shrink faster than the conservative minimax CI as the sample size increases, and the Bonferroni adaptive CIs are shorter than the ones of CL. While the ratio of the length of the Bonferroni CI to that of the CI of CL is larger in this case compared to the case where $J = 2$, especially when $n$ is large, the Bonferroni CI is still more than 50 % narrower than the CI of CL, and not much wider than the oracle minimax CI.

Empirical Illustration

In this section, we apply our procedure to the production function estimation problem for the Chinese chemical industry. Specifically, we use the firm-level data of jacho2010identification for the year 2001, which was also used by horowitz2017nonparametric to illustrate their method of constructing the uniform confidence band for the production function under shape restrictions.

In the dataset, the dependent variable is the logarithm of value-added real output ($y$), and the explanatory variables are the logarithms of the net value of the real fixed asset ($k$) and the number of employees ($\ell$). After removing the outliers for $y,k$ and $\ell$, the remaining sample size was $n=1,636$.\footnote{We used the conventional way of outlier detection, removing the observations that are greater than the third quantile plus IQR times 1.5, or less than the first quantile minus IQR times 1.5. Our resulting sample size is close to horowitz2017nonparametric, who have $n=1,638$. } Table (ref) shows the brief summary of the variables used in our analysis. We are interested in construction of the confidence interval for $f(k_{0},\ell_{0}) := \operatorname{\mathbf{E}} \left[y|k=k_{0},\ell=\ell_{0}\right]$. We take $\left(k_{0},\ell_{0}\right)$ to be medians of each variable.

The first step is to estimate the variance of the error term. We assume homoskedastic errors for simplicity. The variance estimator is defined as \[ \hat{\sigma}^{2}=\frac{\sum_{i=1}^{n}\left(y_{i}-\hat{r}(k_{i},\ell_{i})\right)}{n-2\nu_{1}+\nu_{2}}, \] where $\hat{r}(k_{i},\ell_{i})$ is the estimator for the conditional mean using kernel regression, $\nu_{1}=\text{tr}(L)$, $\nu_{2}=\text{tr}(L'L)$, where $L$ is the weight matrix for the kernel estimator. Refer to wassermann2006all for a justification for this variance estimator. We used the Gaussian kernel with the bandwidth chosen by expected Kullback-Leibler cross validation as in hurvich1998smoothing.

For the function space, we consider adapting to a sequence of parameter spaces $\{\Lambda_{+,\mathcal{V}}(\gamma_{j},C) \}_{j = 1}^6$ with $\gamma_j = 1 - (j - 1)/5$ for $j = 1,...,5$ and $\gamma_6 = 10^{-3}$. We take $\mathcal{V}=\{1,2\}$, assuming that the production function is nondecreasing in both fixed assets and labor, which is consistent with economic theory. To make $\Lambda_{+,\mathcal{V}}(\gamma_{j},C) \subset \Lambda_{+,\mathcal{V}}(\gamma_{6},C)$ hold for all $j = 1,...,5$, we only use observations in a restricted support, and the effective sample size is given by $n_{\text{eff}} = 272$.

For the norm, we use the Euclidean norm weighted by the inverse of the standard deviation of each input, $\left\lVert (k,\ell)\right\rVert =(({k}/{s_{k}})^{2}+({\ell}/{s_{\ell}})^{2})^{1/2}$ where $s_{k}$ and $s_{\ell}$ are standard deviations of $k$ and $\ell$, respectively. We take conservative values of $C$ by setting \[ C=2\times\underset{(i,j)\in\left\{ 1,...,n_{\text{eff}}\right\} ^{2}}{\max }\frac{\left\lvert y_{j}-y_{i}\right\rvert }{\left\lVert (k_{j},\ell_{j})-(k_{i},\ell_{i})\right\rVert ^{\gamma_{6}}}. \]

We compare different procedures to construct CIs. The methods in comparison are the minimax CI with respect to the largest space $\Lambda_{+,\mathcal{V}}(\gamma_{6},C)$ (row 1), the restricted minimax CI with respect to the smallest space $\Lambda_{+,\mathcal{V}}(\gamma_{1},C)$ (row 2), the adaptive Bonferroni CI adapting to $\{\gamma_j\}_{j = 1}^6$ (row 3), the same adaptive CI, but taking into account the correlations between different CIs (fourth row), and the adaptive CI of cai2004adaptation (henceforth CL) adapting to $\{\gamma_j\}_{j = 1}^6$ (fifth row). Note that all the CIs maintain correct coverage over the largest space $\Lambda_{+,\mathcal{V}}(\gamma_{6},C)$, except for the second one, which is valid only over the smallest space $\Lambda_{+,\mathcal{V}}(\gamma_{1},C)$. We refer to the first minimax CI as the conservative minimax CI.

Table (ref) demonstrates the 95% confidence intervals for $f(k_{0},\ell_{0})$ produced by different inference methods. First of all, the lengths of the adaptive Bonferroni CIs are much shorter than the conservative minimax CI, while the procedure of CL yields a wider CI, almost as long as the conservative minimax CI. We can also observe that the adaptive Bonferroni CI using the calibrated value of $\tau^*$ (fourth row) is relatively narrower than its naive version taking $\tau = 0.05/2J$ (third row). Lastly, while the length of the second minimax procedure (second row) is the shortest, it is only valid when we are confident that the true regression function is in the smallest function space we consider, $ \Lambda_{+,\mathcal{V}}(\gamma_{1},C)$. Together with the simulation results in the previous section, our empirical analysis demonstrates the advantage of using an adaptive procedure when the monotonicity restriction is plausible as well as good finite sample performance of our proposed Bonferroni adaptive procedure.