A simple approach to construct confidence bands for a regression function with incomplete data

TL;DR

This paper introduces a new kernel-based method for constructing asymptotically valid confidence bands for regression functions with incomplete data, addressing challenges of unknown data distribution and missing responses.

Contribution

It proposes a novel kernel-type estimator and derives its limiting distribution to enable accurate confidence band construction under missing data conditions.

Findings

The method achieves asymptotically correct coverage probabilities.

Numerical studies show good finite-sample performance.

The approach handles unknown density and variance functions.

Abstract

A long-standing problem in the construction of asymptotically correct confidence bands for a regression function $m(x)=E[Y|X=x]$, where $Y$ is the response variable influenced by the covariate $X$, involves the situation where $Y$ values may be missing at random, and where the selection probability, the density function $f(x)$ of $X$, and the conditional variance of $Y$ given $X$ are all completely unknown. This can be particularly more complicated in nonparametric situations. In this paper we propose a new kernel-type regression estimator and study the limiting distribution of the properly normalized versions of the maximal deviation of the proposed estimator from the true regression curve. The resulting limiting distribution will be used to construct uniform confidence bands for the underlying regression curve with asymptotically correct coverages. The focus of the current paper is on the case where $X\in \mathbb{R}$. We also perform numerical studies to assess the finite-sample performance of the proposed method. In this paper, both mechanics and the theoretical validity of our methods are discussed.

Figures

Tables

Table 1. Table 1: Coverage and the average area of confidence bands (averaged over 3000 bands) for the regression function m ​ ( x ) 𝑚 𝑥 m(x) ; here κ = 10 − 3 𝜅 superscript 10 3 \kappa=10^{-3} . Nominal coverage probability = 90%.

		$n=200$		$n=500$		$n=1000$
	Missing =	50%	25%	50%	25%	50%	25%
Method
Via (18), with	Coverage =	0.872	0.878	0.891	0.901	0.904	0.905
${ϵ}_i\sim$Unif$(-\kappa ,\kappa )$	(Area) =	1.884453	1.440004	1.416213	1.150745	1.270975	0.820889
Via (18), with	Coverage =	0.872	0.878	0.891	0.901	0.904	0.905
${ϵ}_1,\dots ,{ϵ}_n=0$	(Area) =	1.884451	1.440004	1.416212	1.150744	1.270973	0.820827
Using complete	Coverage =	0.726	0.794	0.731	0.836	0.778	0.859
cases only	(Area) =	1.216772	1.257136	1.040574	0.954291	0.841123	0.713148

Table 2. Table 2: Coverage and the average area of confidence bands (averaged over 3000 bands) for the regression function m ​ ( x ) 𝑚 𝑥 m(x) ; here κ = 10 − 3 𝜅 superscript 10 3 \kappa=10^{-3} . Nominal coverage probability = 95%.

		$n=200$		$n=500$		$n=1000$
	Missing =	50%	25%	50%	25%	50%	25%
Method
Via (18), with	Coverage =	0.932	0.938	0.947	0.957	0.949	0.954
${ϵ}_i\sim$Unif$(-\kappa ,\kappa )$	(Area) =	2.101318	1.658559	1.613458	1.300233	1.423841	0.926563
Via (18), with	Coverage =	0.932	0.938	0.947	0.957	0.949	0.954
${ϵ}_1,\dots ,{ϵ}_n=0$	(Area) =	2.101317	1.658558	1.613456	1.300233	1.423840	0.926561
Using complete	Coverage =	0.841	0.887	0.848	0.918	0.872	0.928
cases only	(Area) =	1.451007	1.443447	1.194201	1.085454	0.955346	0.808382

Equations

$$m_n(x)=\sum _{i=1}^n{Y}_{i}K((x-X_i)/h_n)÷\sum _{i=1}^{n}K((x-X_i)/h_n),$$

$$P\{\Delta =1|X=x,Y=y\}=P\{\Delta =1|X=x\}=E[\Delta |X=x]=:p(x),$$

$${\hat{f}}_n(x)=\frac{1}{n{h}_n}\sum _{i=1}^{n}K((x-X_i)/h_n)\text{mbox}and{\sigma }_n^2(x)=\frac{{\sum }_{i=1}^n{Y}_i^{2}K((x-X_i)/h_n)}{{\sum }_{i=1}^{n}K((x-X_i)/h_n)}-m_n^2(x)$$

$$P\left\{\sqrt{2\delta\log n}\left(\sqrt{\frac{nh_{n}}{c_{K}}}\,\,\sup_{x\in[0,1]}\,\sqrt{\frac{\widehat{f}_{n}(x)}{\sigma^{2}_{n}(x)}}~{}\Big{|}m_{n}(x)-m(x)\Big{|}-d_{n}\right)\,\leq\,y\right\}\rightarrow\exp\left(-2\,e^{-y}\right),$$

$$d_n$$

$$C_1=\frac{1}{2c_K}\left(K^2(A)+K^2(-A)\right)\text{mbox}and{C}_2=\frac{1}{2c_K}\int [K^{\prime }(t){]}^{2}dt.$$

$$m_n(x)\pm {\left(\frac{c_K\cdot {\sigma }_n^2(x)}{n{h}_n\cdot {\hat{f}}_n(x)}\right)}^{1/2}\left(\frac{x^{(\alpha )}}{\sqrt{2\delta \log n}}+d_n\right),x\in [0,1],$$

$$P\left\{m(x)\in m_n(x)\pm {\left(\frac{c_K\cdot {\sigma }_n^2(x)}{n{h}_n\cdot {\hat{f}}_n(x)}\right)}^{1/2}\left(\frac{x^{(\alpha )}}{\sqrt{2\delta \log n}}+d_n\right),x\in [0,1]\right\}⟶1-\alpha ,$$

$${\overline{m}}_n(x):=\sum _{i=1}^n{\Delta }_i{Y}_{i}K((x-X_i)/h_n)÷\sum _{i=1}^n{\Delta }_{i}K((x-X_i)/h_n).$$

$$\frac{E[(\Delta Y|X)}{E(\Delta |X)}\stackrel{\text{mbox}a.s.}{=}\frac{E[E(\Delta Y|X,Y)|X]}{E(\Delta |X)}\stackrel{\text{mbox}by(\text{ref}MAR)}{=}\frac{E[YE(\Delta |X)|X]}{E(\Delta |X)}=E(Y|X),$$

$${\breve{m}}_n(x)=\sum _{i=1}^n\frac{{\Delta }_i{Y}_i}{p(X_i)}K\left(\frac{x-X_i}{h_n}\right)÷\sum _{i=1}^{n}K\left(\frac{x-X_i}{h_n}\right),$$

$$E[\frac{\Delta Y}{p(X)}\big{|}X]\stackrel{{\scriptstyle\mbox{\tiny a.s.}}}{{=}}\frac{1}{p(X)}E[E(\Delta Y|X,Y)|X]\stackrel{{\scriptstyle\mbox{\tiny by (\ref{MAR})}}}{{=}}\frac{1}{p(X)}E[Yp(X)|X]\stackrel{{\scriptstyle\mbox{\tiny a.s.}}}{{=}}E(Y|X).$$

$${\tilde{m}}_n(x)=\sum _{i=1}^n\left\{\left[\frac{{\Delta }_i{Y}_i}{p(X_i)}+{ϵ}_i\right]K\left(\frac{x-X_i}{h_n}\right)\right\}÷\sum _{i=1}^{n}K\left(\frac{x-X_i}{h_n}\right),$$

$${\hat{m}}_n(x)=\sum _{i=1}^n\left\{\left[\frac{{\Delta }_i{Y}_i}{{\hat{p}}_n(X_i)}+{ϵ}_i\right]K\left(\frac{x-X_i}{h_n}\right)\right\}÷\sum _{i=1}^{n}K\left(\frac{x-X_i}{h_n}\right),$$

$${\hat{p}}_n(x)=\sum _{i=1}^n\left\{({\Delta }_i+{ϵ}_i)K\left(\frac{x-X_i}{{\lambda }_n}\right)\right\}÷\sum _{i=1}^{n}K\left(\frac{x-X_i}{{\lambda }_n}\right),$$

$$Y^{*}=\frac{\Delta Y}{p(X)}+\epsilon~{}~{}~{}~{}\mbox{and put}~{}~{}~{}~{}~{}~{}Z=Y^{*}-E[Y^{*}|X]~{}\Big{(}\stackrel{{\scriptstyle\mbox{\tiny a.s.}}}{{=}}Y^{*}-m(X)\Big{)}.$$

$${\hat{\sigma }}_{{\hat{p}}_n}^2(x)=\left\{\sum _{i=1}^n{\left[\frac{{\Delta }_i{Y}_i}{{\hat{p}}_n(X_i)}+{ϵ}_i\right]}^{2}K\left(\frac{x-X_i}{h_n}\right)÷\sum _{i=1}^{n}K\left(\frac{x-X_i}{h_n}\right)\right\}-{\left({\hat{m}}_n(x)\right)}^2,$$

$$P\left\{\sqrt{2\delta\log n}\left(\sqrt{\frac{nh_{n}}{c_{K}}}\,\,\sup_{x\in[0,1]}\,\sqrt{\frac{\widehat{f}_{n}(x)}{\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)}}~{}\,\bigg{|}\widehat{m}_{n}(x)-m(x)\bigg{|}-d_{n}\right)\leq y\right\}\rightarrow\exp\left(-2\,e^{-y}\right),~{}$$

$${\hat{m}}_n(x)\pm {\left(\frac{c_K\cdot {\hat{\sigma }}_{{\hat{p}}_n}^2(x)}{n{h}_n\cdot {\hat{f}}_n(x)}\right)}^{1/2}\left(\frac{x^{(\alpha )}}{\sqrt{2\delta \log n}}+d_n\right),x\in [0,1],$$

$$Y=\sin \left(\pi [X^4+e^{\cos (X)}]\right)+\sigma (X)\cdot Z,\text{mbox}withZ\sim N(0,1),$$

$$U_n$$

$$U=\exp \{-2\exp (-U_n)\}$$

$$V_n$$

$$\sqrt{n{h}_n\log n}\underset{x\in [0,1]}{\sup }\left|{\hat{m}}_n(x)-{\tilde{m}}_n(x)\right|{⟶}^{p}0.$$

$$\left|{\hat{m}}_n(x)-{\tilde{m}}_n(x)\right|$$

$$\underset{x\in [0,1]}{\sup }\left|{\hat{m}}_n(x)-{\tilde{m}}_n(x)\right|$$

$${\tilde{\sigma }}_n^2$$

$${\sigma }^2(x)$$

$$\underset{x\in [0,1]}{\sup }\left|{\hat{\sigma }}_{{\hat{p}}_n}^2(x)-{\tilde{\sigma }}_n^2(x)\right|$$

$$\underset{x\in [0,1]}{\sup }\left|{\tilde{\sigma }}_n^2(x)-{\sigma }^2(x)\right|$$

$$\left|{\hat{\sigma }}_{{\hat{p}}_n}^2(x)-{\tilde{\sigma }}_n^2(x)\right|$$

Full text

Keywords: Kernel regression; incomplete data; confidence bands.

1 Introduction

Nonparametric regression estimation has important applications in both statistical estimation theory and classical statistical inferential procedures such as tests of hypotheses or the construction of uniform confidence bands for a true regression function. Confidence bands, in particular, provide insight into the variability of the estimators of the entire regression curve and can also be used to study and investigate certain global features, such as the shape, of the true curve. A very standard procedure to construct asymptotically correct confidence bands for a regression function, over a connected compact set, is based on the limiting distribution of the properly normalized versions of the maximal deviation of the regression estimator from the true regression curve. Results along these lines include the work of Johnston (1982) who constructed confidence bands based on kernel regression estimators; Härdle (1989) established uniform confidence bands for M-smoothers; Eubank and Speckman (1993) proposed bias-corrected confidence bands based on nonparametric kernal regression with fixed design points, whereas Xia (1998) considered random design points under dependence. Additionally, bootstrap confidence bands based on nonparametric regression have been proposed by Neumann and Polzehl (1998), Claeskens and Keilegom (2003), and Song et al. (2012). Härdle and Song (2010) constructed uniform confidence bands for a quantile regression curve with a one-dimensional predictor, whereas Cai et al. (2014) constructed adaptive confidence bands based on nonparametric regression functions. Massé and Meiniel (2014) developed adaptive confidence bands for the case of nonparametric fixed design regression models, and Proksch (2016) developed uniform confidence bands in a nonparametric regression setting with deterministic and multivariate predictor. Another related result is that of Gu and Yang (2015).

The papers cited above as well as most of the results in the literature deal with the cases where the data are fully observable. The focus of this paper is on the realistic case where the response variable could be unobservable or missing. More specifically, let $(X,Y)$ be a random pair with the cumulative distribution function (cdf) $G(x,y)$. Here, $X$ is a $d$-dimensional random vector of covariates and $Y$ is the response variable influenced by $X$. Given the independent and identically distributed (iid) data $(X_{1},Y_{1}),\dots,(X_{n},Y_{n})$ from $G(x,y)$, let $m_{n}(x)$ be the Nadaraya-Watson (Nadaraya (1970), Watson (1964)) kernel regression estimator of the regression function $m(x):=E[Y|X=x]$, i.e.,

[TABLE]

where 0/0 := 0 by convention. Here $K:\mathbb{R}^{d}\to\mathbb{R}_{+}$ is the kernel used with the bandwidth $h_{n}$. However, our focus in this paper is on the case of $d=1$. Now, for various reasons, some of the $Y_{i}$’s may be unavailable or missing from the data. Missing data are common in opinion polls, survey data, mail questionnaires, data collected in medical research and other scientific studies. In this paper we consider the case where $Y$ may be Missing At Random (MAR). More specifically, let $\Delta$ be the Bernoulli random variable defined as $\Delta=0$ if $Y$ is missing and $\Delta=1$, otherwise. Then, the MAR assumption states:

[TABLE]

i.e., the probability that $Y$ is missing does not depend on $Y$ itself. For more on this and other types of missing probability mechanism see, for example, Little and Rubin (2002). Here, the missing probability mechanism $p(x)$, also called the selection probability, is assumed to be completely unknown. In the rest of this paper, the iid data will be represented by $(X_{1},Y_{1},\Delta_{1}),\dots,(X_{n},Y_{n},\Delta_{n})$. Some work has been done for the simpler problem of constructing confidence intervals for $m(x_{0})=E[Y|X=x_{0}]$, where $x_{0}$ is a given point in $\mathbb{R}^{d}$; see, for example, Qin et al. (2014) as well as Lei and Qin (2011). However, to the best of our knowledge, a long-standing problem in constructing uniform confidence bands for the regression function $m(x)$, over compact sets, involves the situation where the response variable $Y$ may be missing at random and the function $p(x)$ in (2), the density function $f(x)$ of $X$, and the conditional variance $\sigma_{0}^{2}(x)=E[Y^{2}|X=x]-m^{2}(x)$ are all completely unknown. Of course, it should be possible to form asymptotically correct uniform confidence bands for $m(x)$ under the restrictive assumptions that $p(x)$, $f(x)$, and $\sigma_{0}^{2}(x)$, (or certain functions of these quantities) are known. However, since such assumptions are unrealistic and not warranted in practice, they will not be pursued in this paper.

In passing we also note that in the case of censored data, Hollander et al. (1997) proposed confidence bands for survival functions based on the empirical likelihood method. Li and van Keilegom (2002) constructed confidence bands for the conditional survival function under random censorship. Wang and Shen (2008) constructed confidence bands of a conditional survival function when the censoring indicators are missing at random. Mondal and Subramanian (2016) developed simultaneous confidence bands for Cox regression in a semiparametric random censorship setup. In another closely related result, Wang and Qin (2010) constructed empirical likelihood confidence bands for a distribution function with missing responses.

In the next section we propose a new kernel-type regression estimator with missing response variables. Our main result in Theorem 2 deals with the limiting distribution of the properly normalized versions of the maximal deviation of the proposed estimator from the true regression curve. Theorem 2 may be viewed as a counterpart of the classical result of Liero (1982) for the setup with no missing data; see Theorem 1 in Section 2.1. Our results will be used to develop a new effective procedure for constructing uniform confidence bands for a regression function in the presence of missing response variables with asymptotically correct coverages. Our numerical results also confirm the finite-sample effectiveness of our procedures.

2 Main results

2.1 Preliminaries and the background tools

To provide the necessary background tools, let $m_{n}(x)$ be the kernel regression estimator defined in (1). Also, let

[TABLE]

be, respectively, the kernel estimators of the density $f$ of $X$ and the conditional variance $\sigma_{0}^{2}(x)=E[Y^{2}|X=x]-m^{2}(x)$. The limiting distribution of the properly normalized versions of the statistic   \sup_{x\in[0,1]}\,\sqrt{\widehat{f}_{n}(x)/\sigma^{2}_{n}(x)}~{}\big{|}m_{n}(x)-m(x)\big{|} have been studied by many authors; see, for example Wandl (1980), Johnston (1982), Liero (1982), and Härdle (1990). In passing, we also note that the interval $[0,1]$ may be any connected compact subset of the interior of the support of $f$. For the case where $X$ is a $d$-dimensional vector, one may refer to the results of Konakov and Piterbarg (1984) and those of Muminov (2011, 2012). To state our proposed estimators and results, we first state a number of classical assumptions, some of which will also be used throughout this paper. These assumptions are virtually all the same as those in Liero (1982).

Assumption (A). The random pair $(X,Y)$ has a probability density function (pdf), $g(x,y)$, with respect to the Lebesgue measure. The random variable $Y$ is almost surely bounded, i.e., $P\{B_{1}\leq Y\ \leq B_{2}\}=1$ for constants $-\infty<B_{1}<B_{2}<\infty$.

Assumption (B). The pdf of $X$, $f$, is strictly positive on $[0,1]$ and vanishes outside of a finite interval $[a,b]$, where $[0,1]\subset(a,b)$.

Assumption (C). The functions $m(x)$, $f(x)$, and $\sigma_{0}^{2}(x)=E[(Y-m(X))^{2}|X=x]$ are twice differentiable with bounded derivatives. Furthermore, $\sigma_{0}^{2}(x)$ is strictly positive on $[0,1]$.

To state the next assumption, put $Z=Y-m(X)$ and let $\widetilde{G}(x,z)$ and $\widetilde{g}(x,z)$ be the cdf and the pdf of the vector $(X,Z)$. Also, let $F$ be the cdf of $X$ and define $H(z|x)$ and $h(z|x)$ to be the conditional cdf and the conditional pdf of $Z$ given $X$, respectively.

Assumption (D).  $\widetilde{g}^{1/2}(x,z)$ is differentiable with respect to both $x$ and $z$, and the partial derivatives are bounded. Furthermore, the inverse functions $H^{-1}$ and $F^{-1}$ of $H$ and $F$ exist and \frac{\partial}{\partial x}H^{-1}\big{(}z|F^{-1}(x)\big{)} and \frac{\partial}{\partial z}H^{-1}\big{(}z|F^{-1}(x)\big{)} are bounded.

Assumption (E). The kernel $K$ is a density function and has a bounded support $[-A,A]$ for some $A>0$. Furthermore, $K$ is continuously differentiable and satisfies $\int x\,K(x)\,dx=0$.

We have the following classical result (see, for example, Liero (1982)).

Theorem 1

Let $h_{n}=n^{-\delta}$, $\frac{1}{5}<\delta<\frac{1}{3}$, and suppose that assumptions (A)-(E) hold. Then

[TABLE]

as $n\to\infty$, where $c_{K}=\int K^{2}(t)\,dt$ and

[TABLE]

with

[TABLE]

The result in (4) can be used to construct confidence bands for $m(x)$. In fact, in light of (4),

[TABLE]

represents an asymptotic $(1-\alpha)100$% confidence band for the regression function $m(x)$ in the sense that, as $n\to\infty$,

[TABLE]

where $x^{(\alpha)}=-(\log\log\left(\frac{1}{1-\alpha}\right)-\log 2)$ is the unique solution of the equation $\exp\{-2\exp(-x)\}=1-\alpha$. Furthermore, one can perform the test of hypothesis $H_{0}:m=m_{0}$, based on the statistic t_{n}:=\sup_{x\in[0,1]}\,\big{(}\widehat{f}_{n}(x)/\sigma^{2}_{n}(x)\big{)}^{1/2}~{}\big{|}m_{n}(x)-m_{0}(x)\big{|}, and reject $H_{0}$, at the significance level $\alpha$, if $t_{n}\,>~{}\sqrt{c_{K}/(nh_{n})}\,\{(2\delta\log n)^{-1/2}\left(\log 2-\log\log\left(\frac{1}{1-\alpha}\right)\right)+d_{n}\}.$

2.2 The proposed estimator

When the response variable $Y_{i}$ may be missing at random according to (2), a very simple counterpart of (1) is the estimator that uses the complete cases only, i.e., the estimator

[TABLE]

The estimator in (10) is in a sense the right estimator. To appreciate this, observe that upon dividing the numerator and the denominator of the right hand side of (10) by the quantity $\sum_{i=1}^{n}K((x-X_{i})/h_{n})$, the estimator $\overline{m}_{n}(x)$ becomes the ratio of the kernel regression estimator of $E(\Delta Y|X=x)$ and the kernel regression estimator of $E(\Delta|X=x)$. Since

[TABLE]

the estimator (10) is indeed a correct kernel-type regression estimator of $m(x)=E(Y|X=x)$. Despite its simplicity, there are no results available in the literature for the maximal deviations of $\overline{m}_{n}(x)$ similar to that in (4). Of course, it may be possible to establish results such as (4) for the maximal deviation of $\overline{m}_{n}(x)$ under the restrictive assumptions that $\sigma_{0}^{2}(x)$ or $p(x)=E[Y|X=x]$ are known, but in this paper we do not impose such assumptions.

In what follows, we propose a kernel regression estimator of $m(x)$ where the presence of missing values is handled via a Horvitz-Thompson-type inverse weighting approach (Horvitz and Thompson (1952)). To motivate our proposed estimator, we first consider the simple but unrealistic case where the selection probability $p(x)=E[\Delta|X=x]$ is known. Now define

[TABLE]

which is the kernel regression estimator of E[\frac{\Delta Y}{p(X)}\big{|}X=x]=E[Y|X=x], where we have used the fact that

[TABLE]

Therefore, when $p(x)$ is known, the estimator in (11) is the kernel regression estimator of $E[Y|X=x]$. When $p(x)$ is unknown it can be replaced by an estimator $\widehat{p}_{n}(x)$; here we have in mind kernel regression estimators of $p(x)$. However, regardless of whether $p(x)$ or $\widehat{p}_{n}(x)$ is used in (11), the maximal deviation of $\breve{m}_{n}(x)$ from $m(x)$ cannot be expected to yield the conclusion of Theorem 1, as given by (4). This is because $\breve{m}_{n}(x)$ in (11) is the kernel regression estimator of $E[Y^{*}|X=x]$, where $Y^{*}:=\frac{\Delta Y}{p(X)}$ does not always have a density with respect to the Lebesgue measure ($P\{Y^{*}=0\}\neq 0$ because $P\{\Delta=0\}>0$), which violates the assumption that the response variable ($Y^{*}$ in this case) has a pdf. To rectify this difficulty we start by artificially adding to $\frac{\Delta Y}{p(X)}$ a zero-mean continuous random variable $\epsilon$, whose pdf has a finite support, and where $\epsilon$ is independent of $(X,Y,\Delta)$. Clearly, the independence of $\epsilon$ and $X$ combined with the MAR assumption in (2) yield E\big{[}\frac{\Delta Y}{p(X)}+\epsilon\big{|}X=x\big{]}=E[Y|X=x]=m(x). The choice of the distribution of $\epsilon$ will be discussed later in Remark 1. Now, let $\epsilon_{1},\dots,\epsilon_{n}$ be iid copies of $\epsilon$, independent of the data $(X_{i},Y_{i},\Delta_{i}),~{}i=1,\dots,n$, and consider the following revised version of (11)

[TABLE]

which is the kernel regression estimator of $m(x)$. Motivated by the naive estimator in (13) and the fact that $p(x)=E[\Delta|X=x]=E[\Delta+\epsilon|X=x]$, our proposed kernel-type regression estimator of $m(x)$ is given by

[TABLE]

where, for technical reasons that will be discussed under Remark 1, we propose to use the following kernel-type estimator of $p(x)$

[TABLE]

instead of the usual kernel estimator $\widehat{p}_{n}(x)=\sum_{i=1}^{n}\Delta_{i}K((x-X_{i})/\lambda_{n})\div\sum_{i=1}^{n}K((x-X_{i})/\lambda_{n})$. Here, $\lambda_{n}\,(\neq h_{n})$ is the smoothing parameter of the kernel. Remark 1 below discusses the use of the random variable $\epsilon$, its justification in the literature, and the choice of its distribution from both theoretical and applied points of view. Next, to establish the limiting distribution of the maximal deviation of $\widehat{m}_{n}(x)$ from $m(x)$, let

[TABLE]

Let $\breve{G}(x,z)$ and $\breve{g}(x,z)$ be the joint cdf and the joint pdf of $(X,Z)$, respectively. Also, let $Q(z|x)$ and $q(z|x)$ be the conditional cdf and the conditional pdf of $Z$ given $X$ and consider the following counterpart of assumption (D):

Assumption (D′).**  $\breve{g}^{1/2}(x,z)$ is differentiable with respect to both $x$ and $z$, and the partial derivatives are bounded. Furthermore, the inverse functions $Q^{-1}$ and $F^{-1}$ of $Q$ and $F$ exist and \frac{\partial}{\partial x}Q^{-1}\big{(}z|F^{-1}(x)\big{)} and \frac{\partial}{\partial z}Q^{-1}\big{(}z|F^{-1}(x)\big{)} are bounded.

Assumption (E′).** The kernel $K$ satisfies Assumption (E). Additionally, $K(x)=K(-x)$ and $K(x)\to 0$, as $|x|\to\infty$.

Regarding the selection probability $p(x)=E[\Delta|X=x]$ we assume

Assumption (F). The selection probability $p(x)=P[\Delta=1|X=x]$ given by (2) is twice differentiable with bounded derivatives. Also, $\inf_{x\in[a,b]}\,p(x)=:\,p_{0}>0$, for some $p_{0}\in(0,1]$, where $[a,b]$ is as in Assumption (B).

Assumption (G). The random variables $\epsilon,\epsilon_{1},\dots,\epsilon_{n}$ are iid zero-mean bounded random variables with a density function that vanishes off the interval $(a_{0},b_{0})$, for some $-\infty<a_{0}<b_{0}<\infty$. Also, $\epsilon_{i}$’s are independent of the data $(X_{i},Y_{i},\Delta_{i})$.

Here, the conditions $K(x)=K(-x)$ and $K(x)\to 0$, as $|x|\to\infty$ that appear under Assumption (E*′) are as in Mack and Silverman (1982), whereas the second part of assumption (F) is standard in missing data literature and essentially amounts to requiring $Y$ to be observable with a non-zero probability for each $X=x$. Assumption (D′*) is the counterpart of assumption (D). Next, let $\widehat{f}_{n}(x)$ be the kernel density estimator defined in (3), and let

[TABLE]

be the kernel regression estimator of the conditional variance \sigma^{2}(x)=E\big{[}\big{(}\frac{\Delta Y}{p(X)}+\epsilon\big{)}^{2}\,\big{|}\,X=x\big{]}-\big{(}E\big{[}\frac{\Delta Y}{p(X)}+\epsilon\,\big{|}\,X=x\big{]}\big{)}^{2}~{}\stackrel{{\scriptstyle\mbox{\tiny(*)}}}{{=}}E\big{[}\big{(}\frac{\Delta Y}{p(X)}+\epsilon\big{)}^{2}\,\big{|}\,X=x\big{]}-m^{2}(x), where $(*)$ follows from (2), (12), and assumption (G). Then we have the following result.

Theorem 2

Let $\widehat{m}_{n}(x)$, $\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)$, and $\widehat{f}_{n}(x)$ be as in (14), (17), (3), respectively. Let $h_{n}=n^{-\delta}$ and $\lambda_{n}=n^{-\beta}$ for any $\delta$ and $\beta$ satisfying $1/5<\beta<\delta<1/3$. Then, under assumptions (A)-(C), (D′), (E′), (F), and (G)

[TABLE]

as $n\to\infty$, where $c_{K}=\int K^{2}(t)\,dt$ and $d_{n}$ is as in (7).

This result immediately yields the following asymptotic (1$$-$$\alpha)100% confidence bands for $m(x)$ when the response variable may be missing at random:

[TABLE]

where $x^{(\alpha)}=-(\log\log\left(\frac{1}{1-\alpha}\right)-\log 2)$.

Remark 1

The main reason for employing the artificial variables $\epsilon_{1},\dots,\epsilon_{n}$ in our methodology above is purely technical (and not quite necessary in numerical studies). The use of artificial or contrived variables in statistical estimation and inference is not new and, in fact, has a long history in the literature. Some classical examples include the problem of nearest neighbor classification when the $d$-dimensional covariate vectors do not have a pdf in which case the dimension is artificially increased to $d+1$ by including an additional random variable $\epsilon$ that has a pdf. This helps to establish the strong consistency of the nearest neighbor classifier and also works as a tie-breaking procedure (see, for example, Devroye et al. (1996, pp 175-176) or Györfi et al. (2002, p. 245)). Another, perhaps more important, example of the use of artificial random variables is related to the weighted bootstrap approximation; see, for example, Mason and Newton (1992), Praestgaard and Wellner (1993), Janssen and Pauls (2003), Janssen (2005), Horváth et al. (2000), Horváth (2000), Burke (1998, 2000), Kojadinovic and Yan (2012), Kojadinovic, Yan, and Holmes (2011), and Mojirsheibani and Pouliot (2017) among others. Since the conditional variance $\sigma^{2}(x)=E[(\frac{\Delta Y}{p(X)}+\epsilon)^{2}|X=x]-m^{2}(x)=[p(X)]^{-1}E[Y^{2}|X=x]+E(\epsilon^{2})-m^{2}(x)$, if $\epsilon$ is chosen to have a large variance, the estimator $\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)$ in (17) can be expected to be inflated. It therefore makes sense to choose $\epsilon$ to have a small variance. In fact, in our numerical work (Section 3), we chose $\epsilon\sim$ Unif$\,[-a,a]$, where $a=10^{-3}$. However, our numerical results also show that one can actually replace $\epsilon_{i}$ by zero in (14) and (15), which is intuitively more appealing, and still expect to see the same numerical results. In other words, the presence of $\epsilon_{i}$ in (14) and (17) is only for theoretical purposes.

3 Numerical results

In this section we carry out some simulation studies to evaluate (numerically) the finite-sample performance of the methods discussed in this paper. The results show that, in general, the proposed estimator performs well. We also take a close look at the performance of the complete-case estimator that is constructed based on the complete cases only. More specifically, in what follows we consider random samples $(X_{1},Y_{1}),\dots,(X_{n},Y_{n})$ of sizes $n=200,500,$ and $1000$ from the model

[TABLE]

where $X\sim N(0.5,1)$ is independent of $Z$, and $\sigma^{2}(x)=1+e^{-(x+2)}$ represents the variance function $E(Y^{2}|X=x)-[E(Y|X=x)]^{2}$. Here, $Y_{i}$ can be missing at random based on one of the following two logistic missing probability models for the function $\pi$ defined in (2):

Model A.   $p(x):=P\{\Delta=1|X=x\}=\exp(1-2x)/\{1+\exp(1-2x)\}$.

Model B.   $p(x):=P\{\Delta=1|X=x\}=\exp(1+0.2x)/\{1+\exp(1+0.2x)\}$,

where $\Delta$=0 if $Y$ is missing (and $\Delta$=1, otherwise). These missing probability mechanisms yield roughly 50% missing data under Model A and about 25% for Model B. Next, for our kernel estimators $\widehat{p}_{n}(x)$ and $\widehat{m}_{n}(x)$ in (15) and (14) and their smoothing parameters $h_{n}=n^{-\delta}$ and $\lambda_{n}=n^{-\beta}$ (subject to $1/5<\beta<\delta<1/3$), we employed the cross-validation approach of Racine and Li (2004) with the Epanechnikov kernel $K(u)=(3/4)(1-u^{2})\cdot I\{|u|\leq 1\}$; this is implemented in the R package “np” (Hayfield and Racine 2008). As for the choice of $\epsilon_{i}$’s that appear in (15) and (14), we considered $\epsilon_{i}\sim\mbox{Unif}(-\kappa\,,\,\kappa)$, $\kappa=10^{-3}$, but we have also considered the more appealing and practical choice of $\epsilon_{i}=0$, $i=1,\dots,n$. Next, to evaluate the performance of various estimators numerically, we computed the statistic

[TABLE]

for sample sizes $n=200,\,500,1000$, and each of the missing probability models A and B, where $\widehat{m}_{n}(x)$, $\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)$, $\widehat{f}_{n}(x)$, $c_{K}$, and $d_{n}$ are as in Theorem 2. In practice, to compute the supremum functional in (19), we used the maximum of \sqrt{\widehat{f}_{n}(x)/\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)}~{}\big{|}\widehat{m}_{n}(x)-m(x)\big{|} over a grid of 200 equally spaced values of $x$ in the interval [0, 1]. Our initial pilot study shows that increasing the grid size to as large as 500 does not make any noticeable changes. Next we note that if $n$ is not very small then by Theorem 2 the quantity

[TABLE]

should be approximately a Unif [0,1] random variable. Repeating the entire above process a total of 3000 times yields $U_{1},\dots,U_{3000}$. We also constructed the above statistic based on the complete cases only, i.e.,

[TABLE]

where $\overline{m}_{n}(x)$, $\bar{f}_{n}(x)$, $\bar{\sigma}^{2}_{n}(x)$ are the estimates of $m(x)$, $f(x)$, and $\sigma^{2}(x)$ based on the complete cases only; see (10) for the definition of the estimator $\overline{m}_{n}(x)$. If we put $V=\exp\{-2\exp(-V_{n})\}$ then the 3000 Monte Carlo runs yield $V_{1},\dots,V_{3000}$. Figure 1 gives plots of the empirical distribution functions of $U_{1},\dots,U_{3000}$ and $V_{1},\dots,V_{3000}$ for different sample sizes, different missing proportions, and the two choices of $\epsilon_{i}$’s. We have also included the $45^{\circ}$ line, which is the CDF of the Unif$\,[0,1]$ random variable.

Comparing plots (a) and (b) in Figure 1, we see that the proposed estimator performs much better than the one based on complete cases when we have 50% missing data and $n=1000$; this is shown by the fact that the empirical CDF of $U_{1},\dots,U_{3000}$ (which corresponds to the proposed estimator) is much closer to the $45^{\circ}$ line. In fact, as Figure 1 shows, the proposed estimator performs better at both 25% and 50% missing rates and for all sample sizes. Of course, the performance of the proposed estimator improves as $n$ increases, confirming the conclusion of Theorem 2. From a practical point of view, the presence of the artificial random variables $\epsilon_{1},\dots,\epsilon_{n}$ in the definition of the estimators (14) and (15) can be viewed as a nuisance. This is only needed as a technical tool and, in practice, one can take $\epsilon_{i}=0$. To confirm this, we also carried out the same simulation study with the more realistic choice of $\epsilon_{i}=0$, $i=1,\dots,n$ and the results were indistinguishable. Plots (m) to (x) in Figure 2 correspond to this setup. As Figure 2 shows, the results are virtually identical.

Next, we used our Monte Carlo simulation results to construct 90% and 95% uniform confidence bands for the regression function $m(x)$, over the set [0, 1]. This resulted in 3000 confidence bands for the regression function $m(x)$ for each sample size $n$ (= 200, 500, 1000) and each missing proportion (25% and 50%). The results summarized in Table 1 are for 90% confidence bands.

Here, coverage is computed as the proportion of confidence bands (out of 3000) that actually captured the true regression function $m(x)$ in the interval $[0,1]$. Table 1 also gives the average area of the 3000 confidence bands constructed under each setup. We make several observations here: (i) The table shows that, numerically, it makes no difference as to whether we consider $\epsilon_{i}\sim\,$Unif$(-10^{-3},10^{-3})$ or the more intuitive and realistic choice of $\epsilon_{1},\dots,\epsilon_{n}=0$ in the estimators (14) and (15). This is clearly evident by the equal coverages in rows one and two of the Table 1 as well as the reported average areas (which are virtually the same up to 6 decimal places). As mentioned earlier, the presence of $\epsilon_{i}$’s in (14) and (15) are only for technical reasons. (ii) Table 1 also shows that the coverage of the bands based on complete cases can reduce far more noticeably as the missing rate changes from 50% to 25% than that of the proposed method. In other words, the proposed method is not as heavily influenced by the amount of missing data; this can be particularly important when a much larger proportion of the data is missing. (iii) We also note that the average areas of the bands based on our proposed method is somewhat higher than those based on complete cases. However, this does not mean that our bands are unnecessarily “wider” than what the theory suggests. In fact, as Theorem 2 as well as Figures 1 and 2 show, the proposed bands are precisely those that are supported by the theory. Table 2 gives the corresponding results for 95% confidence bands. The conclusions are virtually the same those of Table 1.

4 Concluding remarks

In this article, we have proposed a kernel-type method to construct asymptotically correct uniform confidence bands for an unknown regression function $m(x)=E[Y|X=x]$, over compact sets, where the response variable $Y$ may be missing at random. The proposed method is fully nonparametric in that the selection probability, the density function $f(x)$ of $X$, and the conditional variance of $Y$ given $X$ are all completely unknown. The proposed method is quite straightforward to implement and has good asymptotic properties. Furthermore, our numerical work shows that the proposed method has good finite-sample performance. As explained in Remark 1, the presence of the artificial variables, i.e., $\epsilon_{i}$’s, in our methodology and theoretical results are purely for technical reasons and, in practice, such variables can be taken to be zero.

5 Appendix: Proof of Theorem 2

We prove Theorem 2 in a number of steps.

STEP 1. Let $\widehat{m}_{n}(x)$ and $\widetilde{m}_{n}(x)$ be as in (14) and (13), respectively. Then, as $n\to\infty$, we have

[TABLE]

To show this, put $r_{n}(x)=\sum_{i=1}^{n}|\Delta_{i}Y_{i}|\,K((x-X_{i})/h_{n})/\sum_{i=1}^{n}K((x-X_{i})/h_{n})$ and observe that in view of Assumption (E)

[TABLE]

where $\mathbb{I}_{B}$ denotes the indicator of a set $B$. But, $r_{n}(x)\leq\max(|B_{1}|,|B_{2}|)$, by the second part of Assumption (A). Furthermore, in view of assumptions (A), (B), (C), (E*′*), and Theorem B of Mack and Silverman (1982) one has \sup_{x\in[0,1]}|\widehat{p}_{n}(x)-p(x)|=O_{p}\big{(}\sqrt{\log n/(n\lambda_{n})}\,\big{)}. Therefore,

[TABLE]

where we have used the fact that $p_{0}\leq\lim_{n\to\infty}\inf_{x}\,\widehat{p}_{n}(x)\leq 1$ (which follows by noticing that $-\sup_{x}|\widehat{p}_{n}(x)-p(x)|+p_{0}\leq\inf_{x}\widehat{p}_{n}(x)\leq\sup_{x}|\widehat{p}_{n}(x)-p(x)|+1$ and then taking the limit, as $n\to\infty$); here, the infimums are taken over the set [-Ah_{n}\,,\,\mbox{1+$$A}h_{n}]. Now (20) follows from (5) together with the fact that $(\log n)^{2}h_{n}/\lambda_{n}=(\log n)^{2}\,n^{\beta-\delta}\to 0$, as $n\to\infty$ (because $\beta<\delta$).

STEP 2. Define the quantity

[TABLE]

where $\widetilde{m}_{n}(x)$ is as in (13), and put

[TABLE]

Also, let $\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)$ be as in (17). Then we have

[TABLE]

where $\sigma^{2}(x)$ is as in (23). To establish (24) and (25), first observe that

[TABLE]

But, by the second part of assumption (A), the second part of assumption (F), and the boundedness of the support of the distribution of $\epsilon$, we immediately find $\sup_{x\in[0,1]}|\widetilde{m}_{n}(x)|=O_{p}(1)$. Thus, by (20), $\sup_{x\in[0,1]}|\widehat{m}_{n}(x)|\leq\sup_{x\in[0,1]}|\widehat{m}_{n}(x)-\widetilde{m}_{n}(x)|+\sup_{x\in[0,1]}|\widetilde{m}_{n}(x)|=o_{p}\left((nh_{n}\log n)^{-1/2}\right)+O_{p}(1)=O_{p}(1)$. Therefore, by (20),

[TABLE]

Next, observe that

[TABLE]

which follows from the last part of assumption (F), the fact that $p_{0}^{2}\leq\lim_{n\to\infty}\inf_{x}\,\widehat{p}^{2}(x)\leq\sup_{x}\,p^{2}(x)=1$, together with Theorem B of Mack and Silverman (1982), where, the infimum is taken over the set [-Ah_{n}\,,\,\mbox{1+$$A}h_{n}]. Therefore,

[TABLE]

which follows because, in view of assumption (A), the first supremum term on the right side of the above inequality is bounded by $\max(B_{1}^{2},B_{2}^{2})$. Similarly, we have $\sup_{x\in[0,1]}|V_{n,2}(x)|=O_{p}(\sqrt{\log n/(n\lambda_{n})}\,)$, from which (24) follows. The proof of (25) is rather straightforward and goes as follows

[TABLE]

STEP 3. Let $\widehat{f}_{n}(x)$, $\widehat{m}_{n}(x)$, $\widetilde{m}_{n}(x)$, $\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)$, and $\widetilde{\sigma}^{2}_{n}(x)$ be as in (3), (14), (13), (17), and (22), respectively, and write

[TABLE]

where

[TABLE]

To deal with the supremum on the r.h.s of (26), we note that $\widetilde{m}_{n}(x)$ and $\widetilde{\sigma}_{n}^{2}(x)$ that appear in this supremum term are, respectively, the kernel regression estimator of $E(Y^{*}|X=x)$ and the kernel estimator of the conditional variance of $Y^{*}$, as given by (23), based on the iid “data” $(X_{i},Y^{*}_{i}),$ $i=1,\dots,n$, where $Y^{*}$ is given by (16). It is straightforward to see that when assumptions (A) and (F) hold then $P\{B^{*}_{1}\leq Y^{*}\leq B^{*}_{2}\}=1$, where $B_{1}^{*}=\min(0,B_{1})+a_{0}$  and $B_{2}^{*}=\frac{B_{2}}{p_{0}}+b_{0}$, with the constants $B_{1}$ and $B_{2}$ as in assumption (A), and where $a_{0}$ and $b_{0}$ are as in assumption (G). Also, in view of assumptions (A) and (G), the random vector $(X,Y^{*})$ has a pdf. Therefore, when assumption (A) holds for the distribution of $(X,Y)$ then, because of asumption (F), it also holds for the distribution of $(X,Y^{*})$. Similarly, if $\sigma_{0}^{2}(x)$ satisfies assumption (C) then so does $\sigma^{2}(x)$ (in view of assumption (F)); to show this, simply observe that in view of (2) we have $\sigma^{2}(x)\stackrel{{\scriptstyle\mbox{\tiny via (\ref{SIG})}}}{{=}}[(p(x))^{-1}-1]E(Y^{2}|X=x)+E(\epsilon^{2})+\sigma_{0}^{2}(x)$. Therefore, as a consequence of Theorem 1, under assumptions (A), (B), (C), (D*′*), (E), (F), and (G),

[TABLE]

where $P(Y\leq y)=\exp\left\{-2\exp(-y)\right\}$, $y\in\mathbb{R}$, $c_{K}=\int K^{2}(t)\,dt$, and $d_{n}$ is as in (7). Therefore to prove Theorem 2, it is sufficient to show that \sqrt{nh_{n}\log n}\,\big{(}R_{n}(1)+R_{n}(2)\big{)}\to^{p}0, as $n\to\infty$. First we show that $\sqrt{nh_{n}\log n}\,|R_{n}(2)|\to^{p}0$. To show this, observe that by (28)

[TABLE]

Furthermore \big{|}\sup_{x}\sqrt{\widetilde{\sigma}^{2}_{n}(x)/\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)}-1\big{|}\leq\sup_{x}\big{|}\sqrt{\widetilde{\sigma}^{2}_{n}(x)/\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)}-1\big{|}\leq\frac{\sup_{x}\big{|}\widetilde{\sigma}^{2}_{n}(x)-\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)\big{|}}{\inf_{x}\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)}. But by (24), $\sup_{x\in[0,1]}\left|\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)-\widetilde{\sigma}^{2}_{n}(x)\right|=O_{p}\left(\sqrt{\log n/(n\lambda_{n})}\right)+o_{p}\left((nh_{n}\log n)^{-1/2}\right)$. We also note that \,\inf_{x}\,\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)\geq\inf_{x}\,\left\{\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)-\sigma^{2}(x)\right\}+\inf_{x}\,\sigma^{2}(x)\geq-\sup_{x}\big{|}\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)-\sigma^{2}(x)\big{|}+\inf_{x}\,\sigma^{2}(x), where $\sigma^{2}(x)$ is as in (23). Similarly, observe that \inf_{x}\,\widehat{\sigma}_{n}^{2}(x)\leq\sup_{x}\big{|}\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)-\sigma^{2}(x)\big{|}+\sup_{x}\,\sigma^{2}(x). Thus we have

[TABLE]

Now, in view of (24) and (25), and upon taking the limit in the above chain of inequalities, as $n\to\infty$, we find $0<\lim_{n\to\infty}\inf_{x}\,\widehat{\sigma}_{n}^{2}(x)<\infty$, which yields

[TABLE]

This in conjunction with (29) imply that $\sqrt{nh_{n}\log n}\,|R_{n}(2)|\to^{p}0$, where $R_{n}(2)$ is as in (27). Next, observe that \sup_{x}\big{|}\widehat{f}_{n}(x)/\widehat{\sigma}_{n}^{2}(x)\big{|}\leq\big{[}\sup_{x}\big{|}\widehat{f}_{n}(x)-f(x)\big{|}+\sup_{x}f(x)\big{]}/\inf_{x}\widehat{\sigma}^{2}_{\widehat{p}_{n}}(x)=O_{p}(1), which follows because $\sup_{x}|\widehat{f}_{n}(x)-f(x)|=o_{p}(1)$ and by the fact that $0<\lim_{n\to\infty}\inf_{x}\,\widehat{\sigma}_{n}^{2}(x)<\infty$ (as shown above). Combining these results, we have

[TABLE]

Putting the above results together, we have $\sqrt{nh_{n}\log n}\,|R_{n}|=o_{p}(1)$. Theorem 2 now follows from this together with (26), (27), and (28).

$\Box$

References

Burke, M.: A Gaussian bootstrap approach to estimation and tests in Asymptotic Methods in Probability and Statistics. E. (eds.) B. Szyszkowicz, pp. 697-706. North-Holland, Amsterdam (1998)

Burke, M.: Multivariate tests-of-fit and uniform confidence bands using a weighted bootstrap. Statist. Probab. Lett. 46, 13-20 (2000)

Cai, T., Low, M., Zongming, M.: Adaptive confidence bands for nonparametric regression functions. J. Amer. Statist. Assoc. 109, 1054-1070 (2014)

Claeskens, G., Van Keilegom, I.: Bootstrap confidence bands for regression curves and their derivatives. Ann. Statist. 31, 1852-1884 (2003)

Devroye, L., Györfi, L., Lugosi, G.: A probabilistic theory of pattern recognition. Springer-Verlag, New York (1996)

Eubank, R.L., Speckman, P.L.: Confidence Bands in Nonparametric Regression. J. Amer. Statist. Assoc. 88 (424), 1287-1301 (2012)

Gu, L., Yang, L.: Oracally efficient estimation for single-index link function with simultaneous confidence band. Electron. J. Stat. 9, 1540-1561 (2015)

Györfi, L., Kohler, M., Krzyżak, A., Walk, H.: A distribution-free theory of nonparametric regression. Springer-Verlag, New York (2002)

Härdle, W.: Asymptotic maximal deviation of M-smoothers. J. Multivariate Anal. 29, 163-179 (1989)

Härdle, W.: Applied Nonparametric Regression. Cambridge University Press (1990)

Härdle, W., Song, S.: Confidence bands in quantile regression. Econometric Theory 26 (4), 1-22 (2010)

Hollander, M., McKeague., I.W., Yang, J.: Likelihood ratio-based confidence bands for survival functions. J. Amer. Statist. Assoc. 92, 215-227 (1997)

Horváth, L. Approximations for hybrids of empirical and partial sums processes. J. Statist. Plann. Inference. 88:1-18 (2000)

Horváth, L., Kokoszka, P., Steinebach, J.: Approximations for weighted bootstrap processes with an application. Statist. Probab. Lett. 48, 59-70 (2000)

Horvitz, D.G., Thompson D.J. A generalization of sampling without replacement from a finite universe. J. Amer. Statist. Assoc. 47, 663-685 (1952)

Janssen, A.: Resampling Student’s t-type statistics. Ann. Inst. Statist. Math. 57, 507-529 (2005)

Janssen, A., Pauls, T.: How do bootstrap and permutation tests work? Ann. Statist. 31, 768-806 (2003)

Johnston, G.J.: Probabilities of maximal deviations for nonparametric regression function estimates. J. Multivariate Anal. 12, 402-414 (1982)

Kojadinovic, I., Yan., J.: Goodness-of-fit testing based on a weighted bootstrap: A fast large-sample alternative to the parametric bootstrap. Canad. J. Statist. 40, 480-500 (2012)

Kojadinovic, I., Yan, J., Holmes, M.: Fast large-sample goodness-of-fit for copulas. Statist. Sinica 21, 841-871 (2011)

Konakov, V.D., Piterbarg, V.I. : On the convergence rate of maximal deviation distribution. J. Multivariate Anal. 15, 279-294 (1984)

Liero, H.: On the maximal deviation of the kernel regression function estimate. Series Statistics 13, 171-182 (1982)

Lei, Q., Qin, Y.: Confidence intervals for nonparametric regression functions with missing data: multiple design case. J. Syst. Sci. Complex. 24, 1204-1217 (2011)

Little, R.J.A., Rubin, D.B: Statistical analysis with missing data. Wiley, New York (2002)

Li, G., van Keilegom, I.: Likelihood ratio confidence bands in nonparametric regression with censored data, Scand. J. Statist. 2, 547-562 (2002)

Mack, Y.P., Silverman, Z.: Weak and strong uniform consistency of kernel regression estimates. Z. Wahrsch. Verw. Gebiete 61, 405-415 (1982)

Mason, D.M., Newton, M.A.: A rank statistics approach to the consistency of a general bootstrap. Ann. Statist. 20, 1611-1624 (1992)

Massé, P., Meiniel, W.: Adaptive confidence bands in the nonparametric fixed design regression model. J. Nonparametr. Stat. 26, 451-469 (2014)

Mojirsheibani, M., Pouliot, W.: Weighted bootstrapped kernel density estimators in two sample problems.J. Nonparametr. Stat. 29, 61-84 (2017)

Mondal, S., Subramanian, S.: Simultaneous confidence bands for Cox regression from semiparametric random censorship. Lifetime Data Anal. 22, 122-144 (2016)

Muminov, M.S.: On the limit distribution of the maximum deviation of the empirical distribution density and the regression function. I. Theory Probab. Appl. 55, 509-517 (2011)

Muminov, M.S.: On the limit distribution of the maximum deviation of the empirical distribution density and the regression function II. Theory Probab. Appl. 56, 155-166 (2012)

Nadaraya, E.A.: Remarks on nonparametric estimates for density functions and regression curves. Theory Probab. Appl. 15, 134-137 (1970)

Neumann, M.H., Polzehl, J.: Simultaneous bootstrap confidence bands in nonparametric regression. J. Nonparametr. Stat. 9, 307-333 (1998)

Praestgaard, J., Wellner, J.A.: Exchangeably weighted bootstraps of the general empirical process. Ann. of Probab. 21, 2053-2086 (1993)

Proksch, K.: On confidence bands for multivariate nonparametric regression. Ann. Inst. Statist. Math. 68, 209-236 (2016)

Qin, Y., Qiu, T., Lei, Q.: Confidence intervals for nonparametric regression functions with missing data. Comm. Statist. Theory Methods. 43, 4123-4142 (2014)

Racine, J., Li, Q.: Cross-validated local linear nonparametric regression. Statist. Sinica 14, 485-512 (2004)

Hayfield, T., Racine, J.: Nonparametric econometrics: The np package. Journal of statistical software (2008)

Song, S., Ritov, Y., Härdle, W.: Bootstrap confidence bands and partial linear quantile regression. J. Multivariate Anal. 107, 244-262 (2012)

Wandl, H.: On kernel estimation of regression functions. Wissenschaftliche Sitzungen zur Stochastik, WSS-03, Berlin. (1980)

Wang, Q., Shen, J.: Estimation and confidence bands of a conditional survival function with censoring indicators missing at random. J. Multivariate Anal. 99, 928-948 (2008)

Wang, Q., Qin, Y.: Empirical likelihood confidence bands for distribution functions with missing responses. J. Statist. Plann. Inference 140, 2778-2789 (2010)

Watson, G.S.: Smooth regression analysis. Sankhya Ser. A 26, 359-372 (1964)

Xia, Y.: Bias-corrected confidence bands in nonparametric regression. J. R. Stat. Soc. Ser. B. Stat. Methodol. 60, 797-811 (1998)