Estimation and Uniform Inference in Sparse High-Dimensional Additive Models

Estimation and Uniform Inference in Sparse High-Dimensional Additive Models

Introduction

Nonparametric regression offers a way to estimate the relationship $f$ between a target variable $Y$ and input variables $X=(X_1,\ldots,X_p)^T$ without imposing restrictive functional assumptions. This relationship is expressed as \[ Y=f(X_1,\ldots,X_p) + \varepsilon, \] where $\varepsilon$ denotes the random error term satisfying ${\mathrm{E}}[\varepsilon|X]=0$. However, in settings with a large number of regressors $p$, especially in cases where $p$ exceeds the number of observations $n$, the well-known curse of dimensionality makes it practically impossible to estimate the regression function $f(X_1, \ldots, X_p)$. A popular approach to overcome this limitation of nonparametric estimation in practice is to impose an additive structure of the regression function, leading to the following additive models

where $f_j(\cdot)$, $j=1,\ldots,p$, are smooth univariate functions. However, even within this framework, challenges can arise when incorporating a large number of regressors or components, or numerous sieve terms to model the component functions $f_j$. While in many studies, the dimension of the problem is fixed, we follow the recent literature and allow the number of parameters to grow asymptotically with and exceed the sample size. This approach provides a reasonable approximation for various empirical settings involving complex data. The concept of additive models has its origins in the works of friedman1981, stone1985 and hastie1990. Readers seeking an introduction to additive models will find useful textbook treatments in hastie1990 and wood2017. A seminal contribution to the analysis of semi- and nonparametric models, including additive models, was made by Robinson1988. Since then, numerous studies focusing on nonparametric estimation of these models have been published, with contributions from LintonNielsen, Andrews1990, FanJiang, Horowitz2004, HorowitzMammen2006, Mammen2006, Horowitz2012 and many more.

Traditionally, the literature has concentrated on estimation and inference in low-dimensional settings where $p$ is fixed. Recent years, however, have witnessed considerable progress in the understanding and analysis of high-dimensional additive models that allow the number of components to grow with the sample size. For example, the theoretical literature has provided insights into estimation rates in high-dimensional additive models, as in the work of sardy2004, lin2006 and many others ravikumar2009, meier2009, huang2010, koltchinskii2010, kato2012two, petersen2016, lou2016.

The theoretical results in high-dimensional settings rely on a sparsity assumption, which requires that only a small number $s$ of coefficients or components are relevant, meaning they are non-zero. While the number of parameters or covariates exceed the sample size, or the number of parameters can increase along with the sample size, sparsity requires that the number of relevant parameters remains smaller than the sample size. However, it may still increase as the sample size grows, thus introducing additional structure. While high-dimensional additive models have been studied extensively, there has been limited focus on valid inference within them, particularly in terms of constructing valid hypothesis tests or confidence regions. Approaches to construct confidence bands have been proposed by haerdle1989, sun1994, fan2000, claeskens2003 and zhang2010, but only in the widely studied fixed-dimensional setting. Only recently have new results been derived regarding valid inference in high-dimensional additive models. In the remainder of the introduction section, we will review these recent advancements and highlight our contribution to this emerging literature.\\

The primary aim of our paper is to provide a method for constructing uniformly valid inference and confidence bands in sparse high-dimensional models in the sieve framework. In doing so, we contribute to the growing literature on high-dimensional inference in additive models, especially that on debiased/double machine learning. The double machine learning approach belloni2014, Chernozhukov2018 offers a general framework for uniformly valid inference in high-dimensional settings. Similar methods, such as those proposed by vandegeer2014 and zhang2014, have also produced valid confidence intervals for low-dimensional parameters in high-dimensional linear models. These studies are based on the so-called debiasing approach, which provides an alternative framework for valid inference. The framework entails a one-step correction of the lasso estimator, resulting in an asymptotically normally distributed estimator of the low-dimensional target parameter. For a survey on post-selection inference in high-dimensional settings and its generalizations, we refer to spindler2015.\\ In research closely related to ours, Kozbur2021 proposes a post-nonparametric double selection approach for a scalar functional of a single component. We consider the same setting as Kozbur2021, which encompasses a more general additively separable model \[ Y= f_{1}(X_1)+ f_{-1}(X_2,\dots,X_p) + \varepsilon, \] which includes the additive model \[ Y= f_1(X_1) + \ldots + f_p(X_p) + \varepsilon. \] Kozbur2021 focuses on inference on functionals of the form $\theta=a(f_1)$ resulting in pointwise confidence intervals based on a penalized series estimator. Although many of these functionals are of interest, Kozbur2021 approach cannot be easily extended to uniform inference because it considers only differentiable functionals. In contrast, our framework allows us to extend these results and clarify the underlying assumptions. First, building on recent results for inference on high-dimensional target parameters by belloni2018uniformly and Belloni2014b, we are able to construct uniformly valid confidence bands for the whole function $f_1$. Second, the post-nonparametric double selection approach in Kozbur2021 relies on very technical assumptions concerning model selection properties of lasso in Assumption 7 that are difficult to verify (see also the discussion in lu2019 and Example 2 in Kozbur2021). We have much simpler and less technical assumptions since we do not need to rely on model selection properties of lasso. We only make use of the rate of convergence of our lasso estimate to orthogonalize our estimation procedure. In this context, we provide tailored results for uniform lasso estimation that are necessary for establishing valid inference.\\

In a recent study based on the previously mentioned debiasing approach by zhang2014, gregory2016 propose an estimator for the first component $f_1$ in a high-dimensional additive model in which the number of additive components $p$ may increase with the sample size. The estimator is constructed in two steps. The first step entails constructing an undersmoothed estimator based on near-orthogonal projections with a group lasso bias correction. A debiased version of the first step estimator is used to generate pseudo responses $\hat{Y}$. These pseudo responses are then used in the second step, which involves applying a smoothing method to a nonparametric regression problem with $\hat{Y}$ and covariates $X_1$. Under sparsity assumptions concerning the number of nonzero additive components, the so-called oracle property is shown. Accordingly, the proposed estimator in gregory2016 is asymptotically equivalent to the oracle estimator, which is based on the true functions $f_2,\dots,f_p$. The asymptotic properties of the oracle estimator are well understood and carry over to the proposed debiasing estimate, including the methodology for constructing uniformly valid confidence intervals for $f_1$. Importantly, gregory2016 do not explicitly focus on inference, and their analysis requires much stronger assumptions to obtain the oracle property. For example, these assumptions include normally distributed errors independent of $X$, as well as a bounded support of $X$. Similar to our framework, a large set of basis functions is chosen in their work, such as polynomials or splines, to approximate the components $f_1$ and $f_{-1}$. One distinctive feature of our work compared to gregory2016 is that we allow the degree of approximating functions to grow to infinity with increasing sample size.\\

A procedure explicitly addressing the construction of uniformly valid confidence bands for the components in high-dimensional additive models has been developed by lu2019. The authors emphasize that achieving uniformly valid inference in these models is challenging due to the difficulty of directly generalizing the ideas from the fixed-dimensional case. Whereas confidence bands in the low-dimensional case are mostly built using kernel methods, the estimators for high-dimensional sparse additive models typically rely on sieve estimators based on dictionaries. To derive their results, lu2019 combine both kernel and sieve methods to draw upon the advantages of each, resulting in a kernel-sieve hybrid estimator. This is a two-step estimator with tuning parameters for kernel estimation and sieves estimation, such as the bandwidth and penalization levels, which must be chosen by cross-validation. Because of the local structure of the hybrid estimator, the framework of lu2019 differs from ours in that they consider an additive local approximation model with sparsity (ATLAS), in which they need to impose a local sparsity structure.

The advantage of our proposed estimator is that we do not have to leave the sieves framework while establishing the uniform validity of the resulting confidence bands. Interestingly, lu2019 conclude that "it is challenging to study the uniform confidence band through pure sieve-type approaches" (p. 4). However, this is exactly our main contribution: achieving uniform inference within a sieve framework. We accomplish this by framing the problem as one of high-dimensional Z-estimation, utilizing recent results from belloni2018uniformly. Additionally, we offer a theory-driven approach for choosing the penalization level for the lasso estimation steps, eliminating the need for computationally intense cross-validation procedures. Similar to gregory2016, lu2019 assume normally distributed errors that are independent of $X$. In contrast, our model framework allows us to dispense with the normality assumption and only requires the existence of the fourth moment of the error term. Moreover, our main results are also compatible with a heteroskedastic error. Additionally, we can overcome the requirement outlined in lu2019 that the number of relevant regressors is bounded. Instead, our model allows the number of relevant regressors to grow to infinity with increasing sample size.

Organization of the Paper

The paper is organized as follows. Section (ref) introduces and motivates the main regression problem in a high-dimensional additive model. Section (ref) presents the estimation method. In Section (ref), the main result is provided. Section (ref) presents a simulation study, highlighting the small sample properties and implementation of our proposed method. Section (ref) provides the proof of the main theorem. The Appendix includes additional technical material. Appendix (ref) presents a general result for uniform inference on a high-dimensional linear functional. Appendix (ref) provides results in terms of uniform lasso estimation rates in high-dimensions which might be of independent interest. Computational details and additional simulation results are presented in Appendix (ref) and Appendix (ref).

Notation

Throughout the paper, we consider a random element $W$ taken from a common probability space $(\Omega,\mathcal{A},P)$. We denote by $P\in \mathcal{P}_n$ a probability measure out of a large class of probability measures, which may vary with the sample size (because the model is allowed to change with $n$). By $\mathbb{P}_n$ we denote the empirical probability measure. Additionally, let $\mathbb{E}$ and $\mathbb{E}_n$ be the expectations with regard to $P$, and $\mathbb{P}_n$, respectively. $\mathbb{G}_n(\cdot)$ denotes the empirical process

for a class of suitably measurable functions $\mathcal{F}:\mathcal{W}\to \mathbb{R}$. $\|\cdot\|_{P,q}$ denotes the $L^q(P)$-norm. Furthermore, $\|v\|_1=\sum_{l=1}^p |v_l|$ denotes the $\ell_1$-norm, $\|v\|_2=\sqrt{v^Tv}$ the $\ell_2$-norm and $\|v\|_0$ equals the number of non-zero components of a vector $v\in\mathbb{R}^p$. We define $v_{-l}:=(v_1,\dots,v_{l-1},v_{l+1}\dots,v_p)^T\in\mathbb{R}^{p-1}$ for any $1\le l\le p$. $\|v\|_\infty=\sup_{l=1,\dots,p }|v_l|$ denotes the $\sup$-norm. Let $c$ and $C$ denote positive constants that are independent of $n$ with values that may change at each appearance. The notation $a_n\lesssim b_n$ means $a_n\le Cb_n$ for all $n$ and some $C$. Furthermore, $a_n=o(1)$ denotes that there exists a sequence $(b_n)_{\ge 1}$ of positive numbers such that $|a_n|\le b_n$ for all $n$, where $b_n$ is independent of $P\in\mathcal{P}_n$ for all $n$, and $b_n$ converges to zero. Finally, $a_n=O_P(b_n)$ means that for any $\epsilon>0$ there exists a $C$ such that $P(a_n>Cb_n)\le\epsilon$ for all $n$.

Setting

Motivation and Illustration

Before introducing the formal setting in Section (ref), we would like to motivate the basic ideas with a simplified example. We consider an additive model with two covariables

Our goal is to perform valid inference on the object of interest $f_1$. In other words, we would like to provide a uniform confidence band for $f_1$ as illustrated in Figure (ref). Hence, we treat $f_2$ in (ref) as a nuisance function. Next, we assume that it is possible to represent the two components through an approximately linear representation. For the first component, the representation is given by $$f_1(X_1) = \theta_0^T g(X_1) + b_1(X_1).$$ Here, $g(X_1)$ is a basis (e.g., a spline basis, sieve terms or a polynomial series) consisting of $d_1$ terms, $\theta_0$ is the corresponding coefficient vector and $b_1(X_1)$ is an approximation error. For the second component, we assume an analogous representation

where the basis $h(X_2)$ consists of $d_2$ terms. We allow the dimensions $d_1$ and $d_2$ to grow with the sample size in order to establish the asymptotic results in high-dimensional settings. Accordingly, the number of components $p$ can also grow with the sample size.\footnote{However, for the ease of notation we will later subsume them in the component $f_{-1}(X_2,\dots,X_p)$.} To derive a uniformly valid confidence band, estimation and inference of

is required. In a naive approach to estimate the first component $\theta_{0,1}$ of the vector $\theta_0$, one might employ lasso or other machine learning methods to select the relevant components of the basis expansion terms in $\theta_0$ and $\beta_0$ in the regression model ((ref)). A possible second step could involve estimating a regression of the dependent variable on all components that have been selected in the lasso estimation step. The final estimator for the first component $\theta_{0,1}$ could then be obtained from this regression, and this procedure could be repeated iteratively for all other components in $f_1(x_1)$.

However, this approach is problematic because it generally fails to yield valid results. In other words, using the lasso estimator necessitates accounting for variable selection, leading to the non-standard problem of post-selection inference. Intuitively, the naive procedure based on a single lasso estimation step could result in selection mistakes that invalidate the resulting confidence interval for $\theta_{0,1}$. These critical selection mistakes will not occur for variables with high predictive power for the dependent variable $Y$, but rather refer to variables that are potentially highly correlated with the basis expansion term $g_1()$. This could lead to an omitted variable bias, preventing the estimator from asymptotically converging to a normal distribution. To overcome this limitation of post-selection inference, the statistical literature has developed the double machine learning and the debiasing approach, as mentioned in the introduction. To address the potential bias introduced by lasso estimation, belloni2014 propose including an auxiliary regression for the corresponding covariate of the target parameter. Here, we consider

where $\nu$ is an error term and $Z_{-1}$ is defined as $$Z_{-1}=(g_2(X_1),\dots,g_{d_1}(X_1),h_1(X_{-1}),\dots,h_{d_2}(X_{-1}))^T.$$ Later, we will also allow for an approximation error in this equation. belloni2014 propose including in the final regression not only the covariates selected in the first step of the naive approach but also augmenting this set of variables with Lasso-selected regressors from the auxiliary regression. This procedure is equivalent to constructing a so-called Neyman orthogonal moment function with respect to the nuisance part. This is essential for ensuring valid post-selection inference for the first component of the vector $\theta_0$. In Section (ref), we will provide more details about this property. Heuristically, the additional regression step in Equation (ref) will lead to robustness against moderate selection mistakes. It can be shown formally, that this procedure implements an orthogonal moment equation $$\mathbb{E}\left[\psi_1(W,\theta_{0,1},\eta_{0,1})\right]=0,$$ where the first component of $\theta_0$ is our target parameter and all other involved parameters are considered nuisance parameters, denoted $\eta_{0,1}$. belloni2014 developed an approach for valid inference for one parameter. In high-dimensional additive models, the major technical challenge arises from the need to conduct inference for the potentially high-dimensional vector $\theta_0$. In other words, the number of elements in $\theta_0$ for which we would like to construct a valid confidence region is allowed to grow with the sample size. Each component of $\theta_0$, $\theta_{0,l}$ with $\l=1, \ldots, d_1$, is determined by an orthogonal moment condition and we will demonstrate how uniformly valid confidence bands can be constructed by embedding into a high-dimensional Z-estimation framework. Finally, we illustrate how the estimation of $\theta_0$ can be translated into uniformly valid confidence bands for the target function $f_1\approx \theta_0^T g(\cdot)$ using a multiplier bootstrap procedure.\\ Figure (ref) illustrates the use of our estimation procedure, which is described in Section (ref), by providing a preview of our simulation studies in Section (ref). Our estimation methods yield an unbiased estimate for the target component $f_1$, as well as corresponding $(1-\alpha)$ confidence interval that is valid over a compact interval of $X_1$. In the next section, we consider a more general additively separable model and introduce a general formulation of the underlying problem.

Formal Setting

Consider the following nonparametric additively separable model $$Y=f(X)+\varepsilon=f_1(X_1)+f_{-1}(X_{-1}) + \varepsilon$$ with $\mathbb{E}[\varepsilon|X]=0$ and $c\le \operatorname{Var}(\varepsilon|X)\le C$. Let the scalar response $Y$ and features $X=(X_1,\dots,X_p)$ take values in $\mathcal{Y}$ and in $\mathcal{X}=\mathcal{X}_1\times\dots\times\mathcal{X}_p$, respectively. We assume the observation of $n$ i.i.d. copies $(W^{(i)})_{i=1}^n= (Y^{(i)},X^{(i)})_{i=1}^n$ of $W=(Y,X)$, with the number of covariates $p$ allowed to grow with the sample size $n$. To ensure identifiability, we assume $\mathbb{E}[f_{1}(X_{1})]=0$. Our aim is to construct uniformly valid confidence regions for the first nonparametric component of the regression function. In other words, we want to find functions $\hat{l}(x)$ and $\hat{u}(x)$ fulfill

Here, $I\subseteq\mathcal{X}_1$ is a bounded interval of interest in which we want to conduct inference. We approximate $f_1$ and $f_{-1}$ using linear combinations of approximating functions $g_1,\dots,g_{d_1}$ and $h_1,\dots,h_{d_2}$, respectively. Define $$g(x_1):=(g_1(x_1),\dots,g_{d_1}(x_1))^T$$ for $x_1\in\mathbb{R}^{}$ and $$h(x_{-1}):=(h_1(x_{-1}),\dots,h_{d_2}(x_{-1}))^T$$ for $x_{-1}\in\mathbb{R}^{p-1}$. It is important to note that we allow the number of approximating functions $d_1$ and $d_2$ to increase with sample size. Assume that the approximations are given by

where $\theta_{0,l}\in\Theta_l$ and analogously

where $b_1$ and $b_2$ denote the error terms. Additionally, it is convenient to define the union of approximation functions $$z(x):=(g_1(x_1),\dots,g_{d_1}(x_1),h_1(x_{-1}),\dots,h_{d_2}(x_{-1}))^T$$ for $x \in\mathbb{R}^p$, where we abbreviate $Z:=z(X)$. As a common assumption in high-dimensional statistics, we will impose sparsity in the linear approximation of $f$, meaning that only a relatively small number of the coefficients in $(\theta_0,\beta_0)$, denoted $s_1$, are non-zero which implies that only a subset of approximation functions in $Z$ is needed to approximate $f_1$ and $f_{-1}$ sufficiently well. Of course, the identity of these relevant approximation functions is not known. Sparsity can occur because of two different reasons: Either we have many irrelevant covariables $X_2,\ldots,X_p$ in our model or we have sparse representations of the regression functions. For each element $g_l$ of $g$, we also consider

with $\mathbb{E}[\nu^{(l)}Z_{-l}] = 0$, such that $(\gamma^{(l)}_0)^TZ_{-l}$ is the linear projection of $g_l(X_1)$ onto $Z_{-l}$. Here, $\gamma^{(l,1)}_0$ denotes the sparse part of the coefficient vector $\gamma^{(l)}_0$. The auxiliary regression ((ref)) is used to construct an orthogonal score function for valid inference in a high-dimensional setting, as described in Section (ref). Estimating $$f_1(\cdot)\approx \theta_0^Tg(\cdot)$$ can be recast into a general Z-estimation problem of the form $$\mathbb{E}\left[\psi_l(W,\theta_{0,l},\eta_{0,l})\right]=0,\quad l\in1,\dots, d_1$$ with target parameter $\theta_0$ where the score functions are defined by

Here, $\eta=(\eta^{(1)},\eta^{(2)},\eta^{(3)})^T$ with $\eta^{(1)}\in \mathbb{R}^{d_1+d_2-1},\eta^{(2)}\in \mathbb{R}^{d_1+d_2-1}$ and $\eta^{(3)}\in \ell^{\infty}(\mathbb{R}^p)$ are nuisance functions. The true nuisance parameter $\eta_{0,l}$ is given by

where $\beta_0^{(l)}$ is defined as $$\beta_0^{(l)}:=(\theta_{0,1},\dots,\theta_{0,l-1},\theta_{0,l+1},\dots\theta_{0,d_1},\beta_{0,1},\dots,\beta_{0,d_2})^T.$$ Essentially, the index $l$ determines which coefficient is not contained in $\beta_0^{(l)}$. Formally, we will assume sparsity of $\beta_0^{(l)}$, i.e. $\|\beta_0^{(l)}\|_0\le s_1$ and we also assume sparsity in the auxiliary regression, i.e. $\|\gamma_0^{(l,1)}\|_0\le s_2$ for all $l=1,\ldots,d_1$. These assumptions are discussed in detail in Comment (ref) and Comment (ref). The third part of the nuisance functions captures the error made by the approximation of $f_1$ and $f_{-1}$, which is independent from $l$. Therefore, we sometimes omit $l$.

Estimation

In this section, we describe our estimation method and how the uniform valid confidence bands are constructed. First, nuisance functions are estimated using lasso regressions. Subsequently, they are plugged into the moment conditions and solved for the target parameters, which yield an estimate $\hat{f}_1$ for the first component in the additive regression model. The lower and upper curve of the confidence bands are then determined using the estimated covariance matrix and a critical value, which is obtained through a multiplier bootstrap procedure. The technical details for the estimation are given in this section.\\

Estimation and Algorithm

Let $ g(x)=\left(g_1(x),\dots,g_{d_1}(x)\right)^T\in\mathbb{R}^{d_1\times 1}, $ and

for some vectors $\theta=(\theta_1,\dots,\theta_{d_1})^T$ and $\eta=(\eta_1,\dots,\eta_{d_1})^T.$ For each $l=1,\dots,d_1$, let $\hat{\eta}_l=\left(\hat{\eta}_l^{(1)},\hat{\eta}_l^{(2)},\hat{\eta}_l^{(3)}\right)$ be an (lasso) estimator of the nuisance function. The estimator $\hat{\theta}_0$ of the target parameter $\theta_0=(\theta_{0,1},\dots,\theta_{0,d_1})^T$ is defined as the solution of

where $\epsilon_{n}=O\left(\delta_n\varsigma_n^{-1/2}n^{-1/2}\right)$ is the numerical tolerance. The sequences of positive constants, $\delta_n$ and $\varsigma_n$, are defined in Assumption A.(ref) and Assumption A.(ref). Since the score is linear, the explicit solution is given by $$\hat{\theta}_{0,l}=-\frac{\mathbb{E}_n[\psi^b_l(W,\hat{\eta}_l)]}{\mathbb{E}_n[\psi^{a}_l(W,\hat{\eta}^{(2)}_l)]}.$$ Therefore, (ref) obviously holds if we use the estimator $\hat{\theta}_{0,l}$. We consider the representation in (ref) because we rely on more general results derived in Appendix (ref) to analyze the theoretical properties of our new estimation procedure. Lastly, the target function $f_1(\cdot)$ can be estimated by

Define the Jacobian matrix

with

for all $l=1,\dots,d_1$. Observe that $\Sigma_{\varepsilon\nu}:=\mathbb{E}\big[\psi(W,\theta_0,\eta_0)\psi(W,\theta_0,\eta_0)^T\big]$ is the covariance matrix of $\varepsilon\nu:=(\varepsilon\nu^{(1)},\dots,\varepsilon\nu^{(d_1)})$. Define the approximate covariance matrix

with

The approximate covariance matrix can be estimated by replacing every expectation with its empirical analog and plugging in the estimated parameters

This estimated covariance matrix can be used to construct the confidence bands

where $c_{\alpha}$ is a critical value determined by the following standard multiplier bootstrap method introduced in chernozhukov2013gaussian. Define $$\hat{\psi}_x(\cdot):=(g(x)^T\hat{\Sigma}_n g(x))^{-1/2}g(x)^T\hat{J}^{-1}\psi(\cdot,\hat{\theta}_{0},\hat{\eta}_{0})$$ and let

where $(\xi^{(i)})_{i=1}^n$ are independent standard normal random variables (especially independent from the data $(W^{(i)})_{i=1}^n$). The multiplier bootstrap critical value $c_\alpha$ is given by the $(1-\alpha)$-quantile of the conditional distribution of $\sup_{x\in I}|\hat{\mathcal{G}}_x|$ given $(W^{(i)})_{i=1}^n$. This estimation procedure can be summarized in Algorithm (ref).

Practical Considerations: Estimation of Sparse High-Dimensional Additive Models

In addition to the significance level $\alpha$, the number of bootstrap repetitions and an interval $I$ for inference, it is necessary to specify a dictionary of approximating functions $g_1, \ldots, g_{d_1}$ for $f_1$ and $h_1, \ldots, h_{d_2}$ for $f_{-1}$ as inputs to Algorithm (ref). In our theoretical results and our simulation experiments, we focus on B-splines. From a practical point of view, B-splines are specified by parameters for a polynomial degree and the number of knots. In our simulation study, we use cubic B-splines. The number of knots plays a role for the quality of the approximation for $f_1$ as well as the width of the confidence bands. For example, a higher number of knots leads to a higher number of basis functions and corresponding number of fitted coefficients. In our simulation study, we observed that our results are relatively robust with regard to the number of knots and only change moderately with changes of the specification. We base the choice of the number of knots on preliminary evidence in terms of the approximation quality for $f_1$. We also investigated the use of a cross-validated choice of the number of knots, which would not be backed by our theoretical framework, but could serve as a practical guidance. Our results, which are available on request, indicate that a cross-validated choice of the B-spline specification tends to lead to conservative results, i.e., wide confidence intervals with higher coverage than the nominal level.

Based on the provided input, the algorithm proceeds with various nuisance estimation steps. Due to the large number of coefficients, we base estimation on the lasso or post-lasso, i.e., a lasso estimator that is followed by an ordinary least squares estimation step belloni:2013. The performance of regularized estimators depends on the choice of the corresponding penalty parameters. We base the choice of the lasso penalty on theoretical arguments as provided in Appendix (ref). In our software implementation, we use the rigorous lasso learner as provided by the R package hdm hdm. We also experimented with a cross-validated choice of the lasso penalty as provided by the glmnet package glmnet:2010. However, we find that the cross-validated penalty choice was more costly under computational consideration and also less stable. In our simulation experiments, we observe that the estimation performance of the lasso learner has important implications on the accuracy of our point estimator and, accordingly, on the coverage results of the confidence bands. High-dimensional additive models are characterized by a large number of parameters, which have to be estimated with a possibly small number of observations. In our simulation results in Section (ref), we find that the theory-based lasso learner exhibits a good performance in selecting the sparse parameters even in highly correlated settings.

Once, the lasso estimation has been performed, the corresponding residuals are plugged into the variance-covariance matrix. This, in turn, is used to construct the simultaneous confidence bands via a multiplier bootstrap procedure. The latter is based on a random perturbation of the score function, for example, by an i.i.d. standard normal distributed random variable. This procedure is very appealing from a computational point of view as it does not require resampling and reestimation of the parameters, as for example in classical bootstrap procedures. It is generally recommended to use a large number of bootstrap repetitions, $B\ge 500$.

Main Results

Now, we proceed to specifying the conditions required to create uniformly valid confidence bands using Algorithm (ref). To represent $f_1$ and $f_{-1}$ by their approximations in ((ref)) and ((ref)), we need to choose an appropriate set of approximating functions $g=(g_1,\dots,g_{d_1})$ and $h=(h_1,\dots,h_{d_2})$, respectively. In this context, let $\bar{d}_n:=\max(d_1,d_2,n,e)$ and $C$ be a strictly positive constant independent of $n$ and $l$, where $e$ in $\bar{d}_n$ denotes the Euler's number. $(\tilde{A}_n)_{n \geq 1}$ denotes a sequence of positive constants, possibly growing to infinity with $\tilde{A}_n \geq n$ for all $n$. For more details, we refer to Appendix (ref). The number of non-zero coefficients in Equation ((ref)) and ((ref)), is given by the sparsity index $s=\max(s_1,s_2)$. Additionally, we set $t_1:=\sup_{x\in I} \|g(x)\|_0\le d_1$. The definition of $t_1$ is helpful if the functions $g_l$, $l=1,\dots,d_1$, are local, meaning that for any point $x$ in $I$ there are at most $t_1<<d_1$ non-zero functions. Furthermore, $g(I)$ denotes the image of the approximation functions concerning the interval of interest $I$.\\ \\ The following assumptions are valid uniformly in $n\ge n_0$ and $P\in\mathcal{P}_n$:

Assumption A.(ref)$(i)$ contains regularity conditions on $g$. We assume that the supremem of the $\ell_2$-norm of $g(x)$ is bounded, but the infimum is allowed to decrease with sample size, affecting the growth conditions in A.(ref)$(v)$. Assumption A.(ref)$(ii)$ is a condition on the approximation error and is discussed further in Comment (ref). This assumption is mild because the number of approximating functions may increase with sample size. Lastly, Assumption A.(ref)$(iii)$ ensures that any violation of the exact Neyman Orthogonality due to approximation errors is negligible. It is worth noting that if $b_1(X_{1})$ and $b_2(X_{-1})$ are measurable with respect to $Z_{-l}$ (for example, in a linear approximate sparse setting for the conditional expectation) the exact Neyman Orthogonality holds.

Assumptions A.(ref)$(i)$ and $(ii)$ are regularity and sparsity conditions, permitting the number of nonzero regression coefficients $s_1$ and $s_2$ to grow to infinity with increasing sample size. A detailed comment on the sparsity condition is given in Comment (ref). Assumption A.(ref)$(iii)$ contains tail conditions on the approximating functions (and therefore on the original variables) as well as for the error term. Assumption A.(ref)$(iv)$ is an eigenvalue condition that restricts the correlation between the basis elements (and therefore between the original variables), as well as the correlation between the error term $\nu_l$ from the auxiliary regressions. The key innovation of A.(ref)$(iv)$ is that we allow the eigenvalues to vanish as $n$ increases. This is essential if we use an increasing number of local approximation functions $g_l$, $l=1,\dots,d_1$, which implies that the support and hence the variance of each $g_l$ is vanishing with increasing sample size. Typically, we have $\varsigma_n\approx d_1$ as we will illustrate in Comment (ref) for B-Splines. This is in line with the Assumption 5 in kozbur2015 which tackles the same problem of local estimators where $\zeta_0(d_1)=d_1^{1/2}$ is the rate to orthogonalize B-Splines. A related assumption on the restricted eigenvalue can also be found in Assumption A3 in lu2019 where the minimum eigenvalue is of order $O(d_1^{-1})$ as well. Lastly, Assumption A.(ref)$(v)$ provides the growth conditions. These are given in general terms and depend on the choice of the approximation functions. Choosing B-Splines simplifies the growth conditions significantly as we will also outline in Comment (ref).\\

Theorem (ref) provides uniform confidence bands for the whole target component $f_1$. Particularly, with probability $1-\alpha$, we have $$\sup_x|\hat{f}_1(x)-f_1(x)|^2=O\left(\frac{sup_x (g(x)^T\hat{\Sigma}_ng(x))c_\alpha^2}{n}\right)=O\left(\frac{\varsigma_n c_\alpha^2}{n}\right)$$ by Assumption A.(ref)$(iv)$. Hence, the rate of convergence/the width of the confidence band strongly depends on the parameter $\varsigma_n$ and therefore on the number of approximating functions $d_1$ where typically $\varsigma_n\approx d_1$ for local approximating functions.

Simulation Results

We assess the empirical performance of our estimator through a series of simulation experiments based on the settings outlined in gregory2016, which in turn build on the work of meier2009. In our first set of simulations, we evaluate the empirical coverage of our simultaneous confidence bands replicating the data generating process from gregory2016 and meier2009, which incorporates a homoskedastic error term. Additionally, because our estimation framework is also compatible with heteroskedastic error terms, we extend the simulation study to include a heteroskedastic version of the aforementioned design.

Furthermore, we investigate the performance of our estimator and confidence bands in scenarios involving interactions among the nuisance functions $f_{-j}(x_{-j})$. These scenarios are challenging in two regards: First, modeling interactions requires a large number of interaction terms, quickly leading to very high-dimensional settings. Second, specifying interactions among various functions in empirical applications is challenging because of the risk that the chosen model will not precisely match the true underlying data generating process. By exploiting a sparse structure, our theoretical framework can tolerate a moderate degree of misspecification in the nuisance functions. We introduce two interacted designs -- one correctly specified and one misspecified--providing a comprehensive assessment of our methodology. We consider the finite sample performance of our estimator in a high-dimensional additive model of the form

with $i = 1, \ldots, n$ and $j = 1, \ldots, p$. The functions $f_j(x_j)$ are nonlinear for $j=1, 2, 4$ and linear for $j=3$, as can be recognized from their definition in Table (ref). The functions $f_5(x_5), \ldots, f_p(x_p)$ are null components, thus implementing a sparse setting. In all of our considered settings, the regressors $X$ are marginally uniformly distributed on an interval, $I = [-2.5, 2.5]$ with correlation matrix $\Sigma$ with $\Sigma_{k,l}=\rho^{|k-l|}$, $1 \le k,l \le p$. We consider the values $\rho = 0$ and $\rho = 0.5$ because they implement the weakest and strongest correlation structures from the simulations in gregory2016. We generate data sets for scenarios with $n \in \{100, 1000\}$ observations and $p \in \{50, 150\}$ explanatory variables $X_1, \ldots, X_p$.

First, we consider a homoskedastic setting with a standard normally distributed error term $\epsilon \sim N\left(0, 1\right)$, which is identical to the design in gregory2016. Second, we define a comparable heteroskedastic setting by specifying an error term $\epsilon_j\sim N\left(0, \sigma_j(x_j)\right)$ with $\sigma_j(x_j) = \underline{\sigma} \cdot (1 + |x_j|)$ and $\underline{\sigma} = \sqrt{\frac{12}{67}}$. This value of $\underline{\sigma}$ ensures a signal-to-noise ratio that is comparable to the homoskedastic setting. Third, we implement two different interacted scenarios, in which we focus on the estimation of the function $f_1(x_1)$. In the correctly specified interacted scenario, which we will later refer to as Scenario $I$, the specification of the estimated sparse additive model is identical to that of the underlying data generating process. In this model, the function $f_2(x_2)$ is interacted with all other components except for $f_1(x_1)$, making the scenario suitable for our estimation approach. Including the interactions of all $p-1$ nuisance functions leads to a high-dimensional setting. For example, if one specifies cubic B-splines with 3 knots in a setting with $p=150$ covariates, the overall dimension of the estimated model is $d_1+d_2 = 6,228$. In the data generating process of the misspecified interacted Scenario $II$, all non-sparse functions are interacted with each other, except for $f_1(x_1)$. In contrast, the model that is estimated empirically only specifies a basic additive model without any interaction terms. Hence, the interacted Scenario $II$ implements a situation with considerable model misspecification.

In the simulation, we use the estimator and the multiplier bootstrap procedure we proposed in Section (ref) to generate predictions $\hat{f}_j(x_j)$ for the function $f_j(x_j)$ and construct simultaneous confidence bands that are defined in terms of $\hat{l}_j(x_j)$ and $\hat{u}_j(x_j)$. We approximate the functions $f_j(x_j)$ in the additive model by cubic B-splines. Estimation is performed using post-lasso with a theory-based choice of the penalty level as implemented in the R package hdm hdm. Further details related to the implementation in the simulation study can be found in Appendix (ref).\\ Table (ref) presents the empirical coverage achieved by the estimated simultaneous $95\%$-confidence bands in the homoskedastic and heteroskedastic scenario. The results are obtained in $R = 500$ repetitions and the confidence bands are constructed over an interval of values of $x_j$, $I = [-2,2]$. A confidence band is considered to cover the function $f_j(x_j)$ if it entirely contains the true function, or, stated more formally, if for all values of $x_j \in I$ it holds that $\hat{l}_j(x_j)\le f_j(x_j) \le \hat{u}_j(x_j)$.\\

The results can be interpreted as evidence supporting the validity of our inference method in high-dimensional additive models in both homoskedastic and heteroskedastic settings. The coverage of the confidence bands is close to the nominal level across all settings, in particular when the sample size is large ($n=1000$) or if the correlation structure among the covariates is weak ($\rho=0$). We find that a precise approximation of the target function $f_j(x_j)$ through an appropriate specification of the employed B-Splines is essential for the coverage results of the confidence bands.

Figures (ref) to (ref) illustrate the average results for the estimators of $f_1(x_1), \ldots, f_5(x_5)$, along with the corresponding joint confidence bands. The green dashed line and the green shaded areas represent the results for the settings with $\rho = 0$, and the blue curve and shaded areas refer to the setting with $\rho = 0.5$. The red solid line shows the true value of the function $f_j(x_j)$. The corresponding results for the heteroskedastic settings are very similar and, hence, made available in Appendix (ref). The plots depicting the average results provide additional insights into the quality of estimation across different settings. First, they illustrate that the estimation benefits substantially from larger samples. In samples with $n=100$ the confidence bands exhibit noticeable variability and width, whereas in larger samples they become considerably narrower. It is important to note that this observation is expected. For example, in the high-dimensional settings with $n=100$ and $p=150$, cubic B-splines with 3 knots would lead to an overall dimension of $d_1 + d_2 = 900$, which is challenging in terms of estimation. Second, comparing the results across different values of $\rho$ helps elucidate the challenge in terms of approximation quality. For example, Figure (ref) shows the results for $f_4(x_4)$ in the homoskedastic design. Whereas the approximation appears to be very accurate for $\rho = 0$ in the setting with $n=100$, it is considerably worse for $\rho=0.5$. However, in larger samples, the quality of the approximation improves, which is in line with the results of gregory2016. Similar observations can be made for $f_3(x_3)$ and $f_5(x_5)$. As in the previously considered settings, the strength of the correlation of the regressors $X$ plays an important role among the performance of the estimator $\hat{f}_1(x_1)$.

In addition to the homoskedastic and heteroskedastic designs, we investigate the performance of our estimator and confidence bands in interacted simulation settings. The results in Tables (ref) and (ref) and Figures (ref) and (ref) provide simulation evidence in very high-dimensional settings -- both with a correctly (Scenario $I$) and an incorrectly (Scenario $II$) specified additive model. As a consequence of the higher dimensionality and the increased noise in the estimation problem, the confidence intervals are wider than in the basic settings considered previously. In all cases and irrespective of the model misspecification in Scenario $II$, the empirical coverage is close to the nominal level. As in the settings considered before, the strength of the correlation among the regressors $X$ makes estimation of $f_1(x_1)$ more challenging (see for example the top panels in Figure (ref)).

In general, the performance of the estimator and the confidence bands depends on the specification of the cubic B-splines used to approximate the target functions (in terms of the number of knots). In our simulations, we observe that the quality of estimation and width of the confidence bands change only moderately when varying the number of knots. Overall, we conclude that our results are relatively robust to the specifications of the B-splines. We find that the performance of the point estimator for the target component $f_j(x_j)$ relies on variable selection through the lasso estimator. In line with our expectations, the lasso estimator is better at selecting the sparse components of the B-splines in larger samples and settings with weaker correlation structure, i.e., $\rho = 0$. In our simulations, we chose the specification of the B-splines estimator based on preliminary evidence. We also investigated the performance of a cross-validated choice (results available upon request). While the latter lacks theoretical justification, we found it to be relatively conservative, resulting in wide confidence bands. During our simulation experiments, we also investigated the performance of an alternative lasso learner that is based on a cross-validated choice of the penalty term. We found that the cross-validated penalty choice was computationally more expensive and inferior in terms of the estimation performance, as indicated by selecting very few of the sparse components. Overall, we interpret the simulation results as supportive evidence for our estimation approach. The approximation and coverage results for our estimator and confidence bands across different simulation settings are comparable to those provided in gregory2016. As in gregory2016, the approximation quality of our estimator depends on the underlying correlation structure of the covariates $X$, the parametrization of the B-splines, and the variable selection performance of lasso. In their simulation study, gregory2016 report pointwise coverage results, whereas our results refer to the coverage of the true function $f_j(x_j)$ over all points of $x_j$ in an interval $I$ by a simultaneous confidence band.

Proofs