Nonparametric modal regression

1 Introduction

Modal regression (Sager and Thisted, 1982; Lee, 1989; Yao et al., 2012; Yao and Li, 2014)

is an alternate approach to the usual regression methods for exploring the relationship between a response variable

Y

and a predictor variable

X

. Unlike conventional regression, which is based on the conditional mean of

Y

given

X = x

, modal regression estimates conditional modes of

Y

given

X = x

Why would we ever use modal regression in favor a conventional regression method? The answer, at a high-level, is that conditional modes can reveal structure that is missed by the conditional mean. Figure 1 gives a definitive illustration of this point: we can see that, for the data examples in question, the conditional mean both fails to capture the major trends present in the response, and produces unecessarily wide prediction bands. Modal regression is an improvement in both of these regards (better trend estimation, and narrower prediction bands). In this paper, we rigorously study and develop its properties.

Figure 1:
Examples of modal regression versus a common nonparametric regression estimator, local linear regression. In the top row, we show local regression estimate and its associated 95% prediction bands alongside the modal regression and its 95% prediction bands. The bottom row does the same for a different data example. The local regression method fails to capture the structure, and produces prediction bands that are too wide.

Modal regression has been used in transportation (Einbeck and Tutz, 2006), astronomy (Rojas, 2005), meteorology (Hyndman et al., 1996) and economics (Huang and Yao, 2012; Huang et al., 2013). Formally, the conditional (or local) mode set at $x$ is defined as

M (x) = {y : \frac{\partial}{\partial y} p (y | x) = 0, \frac{\partial^{2}}{\partial y^{2}} p (y | x) < 0},

(1)

where $p (y | x) = p (x, y) / p (x)$ is the conditional density of $Y$ given $X = x$ . As a simplification, the set $M (x)$ can be expressed in terms of the joint density alone:

M (x) = {y : \frac{\partial}{\partial y} p (x, y) = 0, \frac{\partial^{2}}{\partial y^{2}} p (x, y) < 0} .

(2)

At each $x$ , the local mode set $M (x)$ may consist of several points, and so $M (x)$ is in general a multi-valued function. Under appropriate conditions, as we will show, these modes change smoothly as $x$ changes. Thus, local modes behave like a collection of surfaces which we call modal manifolds.

We focus on a nonparametric estimate of the conditional mode set, derived from a kernel density estimator (KDE):

{ˆ M}_{n} (x) = {y : \frac{\partial}{\partial y} {ˆ p}_{n} (x, y) = 0, \frac{\partial^{2}}{\partial y^{2}} {ˆ p}_{n} (x, y) < 0},

(3)

where ${ˆ p}_{n} (x, y)$ is the joint KDE of $X, Y$ . Scott (1992) proposed this plug-in estimator for modal regression, and Einbeck and Tutz (2006) proposed a fast algorithm. We extend the work of these authors by giving a thorough treatment and analysis of nonparametric modal regression. In particular, our contributions are as follows.

We study the geometric properties of modal regression.
We prove consistency of the nonparametric modal regression estimator, and furthermore derive explicit convergence rates, with respect to various error metrics.
We propose a method for constructing confidence sets, using the bootstrap, and prove that it has proper asymptotic coverage.
We propose a method for constructing prediction sets, based on plug-in methods, and prove that the population prediction sets from this method can be smaller than those based on the regression function.
We propose a rule for selecting the smoothing bandwidth of the KDE based on minimizing the size of prediction sets.
We draw enlightening comparisons to mixture regression (which suggests a clustering method using modal regression) and to density ridge estimation.

We begin by reviewing basic properties of modal regression and recalling previous work, in Section 2. Sections 3 through 8 then follow roughly the topics described in items 1–6 above. In Section 9 we end with some discussion. Simple R code for the modal regression and some simulation data sets used in this paper can be found at http://www.stat.cmu.edu/~yenchic/ModalRegression.zip.

2 Review of modal regression

Consider a response variable $Y \in K \subseteq R$ and covariate or predictor variable $X \in D \subseteq R^{d}$ , where $D$ is a compact set. A classic take on modal regression (Sager and Thisted, 1982; Lee, 1989; Yao and Li, 2014) is to assume a linear model

M o d e (Y | X = x) = β_{0} + β^{T} x,

where $β_{0} \in R$ , $β \in R^{d}$ are unknown coefficients, and $M o d e (Y | X = x)$ denotes the (global) mode of $Y$ given $X = x$ . Nonparametric modal regression is more flexible, because it allows $M (x)$ to be multi-valued, and also it models the components of $M (x)$ as smooth functions of $x$ (not necessarily linear). As another nonlinear generalization of the above model, Yao et al. (2012) propose an interesting local polynomial smoothing method for mode estimation; however, they focus on the global mode $M o d e (Y | X = x)$ rather than $M (x)$ , the collection of all conditional modes.

The estimated local mode set ${ˆ M}_{n} (x)$ in (3) from Scott (1992) relies on an estimated joint density function ${ˆ p}_{n} (x, y)$ , most commonly computed using a KDE. Let $(X_{1}, Y_{1}), \dots, (X_{n}, Y_{n})$ be the observed data samples. Then the KDE of the joint density $p (x, y)$ is

{ˆ p}_{n} (x, y) = \frac{1}{n h^{d + 1}} n \sum i = 1 K (\frac{∥ x - X_{i} ∥}{} h) K (\frac{y - Y_{i}}{h}) .

(4)

Here $K$ is a symmetric, smooth kernel function, such as the Gaussian kernel (i.e., $K (u) = e^{- u^{2} / 2} / \sqrt{2 π}$ ), and $h > 0$ is the smoothing bandwidth. For simplicity, we have used the same kernel function $K$ and bandwidth $h$ for the covariates and the response, but this is not necessary. For brevity, we will write the estimated modal set as

{ˆ M}_{n} (x) = {y : {ˆ p}_{y, n} (x, y) = 0, {ˆ p}_{y y, n} (x, y) < 0},

(5)

where the subscript notation denotes partial derivatives, as in $f_{y} = \partial f (x, y) / \partial y$ and $f_{y y} = \partial^{2} f (x, y) / \partial y^{2}$ .

In general, calculating ${ˆ M}_{n} (x)$ can be challenging, but for special kernels, Einbeck and Tutz (2006) proposed a simple and efficient algorithm for computing local mode estimates, based on the mean-shift algorithm (Cheng, 1995; Comaniciu and Meer, 2002). A related approach can be found in Yao (2013), where the authors consider a mode hunting algorithm based on the EM algorithm. For a discussion of how the mean-shift and EM algorithms relate, see Carreira-Perpiñán (2007). In general, mean-shift algorithms can be derived for any KDEs with radially symmetric kernels (Comaniciu and Meer, 2002), but for simplicity we assume here that $K$ is Gaussian. The partial mean-shift algorithm of Einbeck and Tutz (2006), to estimate conditional modes, is described in Algorithm 1.

Input: Data samples

D = {(X_{1}, Y_{1}), \dots, (X_{n}, Y_{n})}

, bandwidth

h

. (The kernel

K

is chosen to be Gaussian.)

Initialize mesh points $M \subseteq R^{d + 1}$ (a common choice is $M = D$ ).
For each $(x, y) \in M$ , fix $x$ , and update $y$ using the following iterations until convergence:

$y ⟵ \frac{n \sum i = 1 Y_{i} K (\frac{∥ x - X_{i} ∥}{} h) K (\frac{y - Y_{i}}{h})}{n \sum i = 1 K (\frac{∥ x - X_{i} ∥}{h}) K (\frac{y - Y_{i}}{h})} .$ (6)

Output: The set

M^{\infty}

, containing the points

(x, y^{\infty})

, where

x

is a predictor value as fixed in

M

, and

y^{\infty}

is the corresponding limit of the mean-shift iterations (6).

Algorithm 1 Partial mean-shift

A straightforward calculation shows that the mean-shift update (6) is indeed a gradient ascent update on the function $f (y) = {ˆ p}_{n} (x, y)$ (for fixed $x$ ), with an implicit choice of step size. Because this function $f$ is generically nonconcave, we are not guaranteed that gradient ascent will actually attain a (global) maximum, but it will converge to critical points under small enough step sizes (Arias-Castro et al., 2013).

3 Geometric properties

In this section, we study the geometric properties of modal regression. Recall that $M (x)$ is a set of points at each input $x$ . We define the modal manifold collection as the union of these sets over all inputs $x$ ,

S = {(x, y) : x \in D, y \in M (x)} .

(7)

By the implicit function theorem, the set $S$ has dimension $d$ ; see Figure 2 for an illustration with $d = 1$ (univariate $x$ ).

We will assume that the modal manifold collection $S$ can be factorized as

S = S_{1} \cup \dots \cup S_{K},

(8)

where each $S_{j}$ is a connected manifold that admits a parametrization

S_{j} = {(x, m_{j} (x)) : x \in A_{j}},

(9)

for some function $m_{j} (x)$ and open set $A_{j}$ . For instance, in Figure 2, each $S_{j}$ is a connected curve. Note that $A_{1}, \dots, A_{K}$ form an open cover for the support $D$ of $X$ . We call $S_{j}$ the $j$ th modal manifold, and $m_{j} (x)$ the $j$ th modal function. By convention we let $m_{j} (x) = \emptyset$ if $x \notin A_{j}$ , and therefore we may write

M (x) = {m_{1} (x), \dots, m_{K} (x)} .

(10)

That is, at any $x$ , the values among $m_{1} (x), \dots, m_{K} (x)$ that are nonempty give local modes.

Under weak assumptions, each $m_{j} (x)$ is differentiable, and so is the modal set $M (x)$ , in a sense. We discuss this next.

Lemma 1 (Derivative of modal functions).

Assume that $p$ is twice differentiable, and let $S = {(x, y) : x \in D, y \in M (x)}$ be the modal manifold collection. Assume that $S$ factorizes according to (7), (8). Then, when $x \in A_{j}$ ,

\nabla m_{j} (x) = - \frac{p_{y x} (x, m_{j} (x))}{p_{y y} (x, m_{j} (x))},

(11)

where $p_{y x} (x, y) = \nabla_{x} \frac{\partial}{\partial y} p (x, y)$ is the gradient over $x$ of $p_{y} (x, y)$ .

Proof. Since we assume that $x \in A_{j}$ , we have $p_{y} (x, m_{j} (x)) = 0$ by definition. Taking a gradient over $x$ yields

0 = \nabla_{x} p_{y} (x, m_{j} (x)) = p_{y x} (x, m_{j} (x)) + p_{y y} (x, m_{j} (x)) \nabla m_{j} (x) .

After rearrangement,

\nabla m_{j} (x) = - \frac{p_{y x} (x, m_{j} (x))}{p_{y y} (x, m_{j} (x))} .

$□$

Lemma 1 links the geometry of the joint density function to the smoothness of the modal functions (and modal manifolds). The formula (11) is well-defined as long as $p_{y y} (x, m_{j} (x))$ is nonzero, which is guaranteed by the definition of local modes. Thus, when $p$ is smooth, each modal manifold is also smooth.

To characterize smoothness of $M (x)$ itself, we require a notion for smoothness over sets. For this, we recall the Hausdorff distance between two sets $A, B$ , defined as

H a u s (A, B) = inf {r : A \subseteq B \oplus r, B \subseteq A \oplus r},

where $A \oplus r = {x : d (x, A) \leq r}$ with $d (x, A) = {inf}_{y \in A} ∥ x - y ∥$ .

Theorem 2 (Smoothness of the modal manifold collection).

Assume the conditions of Lemma 1. Assume furthermore all partial derivatives of $p$ are bounded by $C$ , and there exists $λ_{2} > 0$ such that $p_{y y} (x, y) < - λ_{2}$ for all $y \in M (x)$ and $x \in D$ . Then

lim | ϵ | \to 0 \frac{H a u s (M (x), M (x + ϵ))}{| ϵ |} \leq max j = 1, \dots, K ∥ m_{j}^{'} (x) ∥ \leq \frac{C}{λ_{2}} < \infty .

(12)

The proof of this result follows directly from Lemma 1 and the definition of Hausdorff distance, so we omit it. Since $M (x)$ is a multi-valued function, classic notions of smoothness cannot be applied, and Theorem 2 describes its smoothness in terms of Hausdorff distance. This distance can be thought of as a generalized $ℓ_{\infty}$ distance for sets, and Theorem 2 can be interpreted as a statement about Lipschitz continuity with respect to Hausdorff distance.

Modal manifolds can merge or bifurcate as $x$ varies. Interestingly, though, the merges or bifurcations do not necessarily occur at points of contact between two modal manifolds. See Figure 3 for an example with $d = 1$ . Shown is a modal curve (manifold with $d = 1$ ), starting at $x = 0$ and stopping at about $x = 0.35$ , which leaves a gap between itself and the neighboring modal curve. We take a closer look at the joint density contours, in panel (c), and inspect the conditional density along four slices $X = x_{1}, \dots, x_{4}$ , in panel (d). We see that after $X = x_{2}$ , the conditional density becomes unimodal and the first (left) mode disappears, as it slides into a saddle point.

A remark about the uniqueness of the modal manifold collection $S$ in (8): this factorization is unique if the second derivative $p_{y y} (x, y)$ is uniformly bounded away from zero. This will later become one of our assumptions (assumption (A3)) in the theoretical analysis of Section 4. Note that in the left panel of Figure 2, the collection $S$ is uniquely defined, while this is not the case in the right panel (at the points of intersection between curves, the density $p$ has vanishing second derivatives with respect to $y$ ).

Lastly, the population quantities defined above all have sample analogs. For the estimate ${ˆ M}_{n} (x)$ , we define

{ˆ S}_{n} = {(x, y) : y \in {ˆ M}_{n} (x), x \in R} = {ˆ S}_{1} \cup \dots \cup {ˆ S}_{ˆ K},

(13)

where each ${ˆ S}_{j}$ is a connected manifold, and $ˆ K$ is the total number. We also define ${ˆ m}_{j} (x)$ in a similar way for $j = 1, \dots, ˆ K$ . Thus, we can write

{ˆ M}_{n} (x) = {{ˆ m}_{1} (x), \dots, {ˆ m}_{ˆ K} (x)} .

(14)

In practice, determining the manifold memberships and the total number of manifolds $ˆ K$ is not trivial. In principle, the sample manifolds ${ˆ S}_{1}, \dots, {ˆ S}_{ˆ K}$ are well-defined in terms of the sample estimate ${ˆ M}_{n} (x)$ ; but even with a perfectly convergent mean-shift algorithm, we would need to run mean-shift iterations at every input $x$ in the domain $D$ to determine these manifold components. Clearly this is not an implementable strategy. Thus from the output of the mean-shift algorithm over a finite mesh, we usually employ some type of simple post-processing technique to determine connectivity of the outputs and hence the sample manifolds. This is discussed further in Section 7.

4 Asymptotic error analysis

In this section we present asymptotic results about the convergence of the estimated modal regression set ${ˆ M}_{n} (x)$ to the underlying modal set $M (x)$ . Let ${B C}^{k} (C)$ denote the collection of $k$ times continuously differentiable functions with all partial derivatives bounded in absolute value by $C$ . (The domain of these functions should be clear from the context.) Given a kernel function $K : R \to R$ , denote the collection of functions

K = {v \mapsto K^{(α)} (\frac{z - v}{h}) : z \in R, h > 0, α = 0, 1, 2},

where $K^{(α)}$ denotes the $α$ -th order derivative of $K$ .

Our assumptions are as follows.

The joint density $p \in {B C}^{4} (C_{p})$ for some $C_{p} > 0$ .
The collection of modal manifolds $S$ can be factorized into $S = S_{1} \cup \dots \cup S_{K},$ where each $S_{j}$ is a connected curve that admits a parametrization $S_{j} = {(x, m_{j} (x)) : x \in A_{j}}$ for some $m_{j} (x)$ , and $A_{1}, \dots, A_{K}$ form an open cover for the support $D$ of $X$ .
There exists $λ_{2} > 0$ such that for any $(x, y) \in D \times K$ with $p_{y} (x, y) = 0$ , $| p_{y y} (x, y) | > λ_{2}$ .
The kernel function $K \in {B C}^{2} (C_{K})$ and satifies

$\int_{R} (K^{(α)})^{2} (z) d z < \infty, \int_{R} z^{2} K^{(α)} (z) d z < \infty,$

for $α = 0, 1, 2$ .
The collection $K$ is a VC-type class, i.e. there exists $A, v > 0$ such that for $0 < ϵ < 1$ ,

$sup Q N (K, L_{2} (Q), C_{K} ϵ) \leq {(\frac{A}{ϵ})}^{v},$

where $N (T, d, ϵ)$ is the $ϵ$ -covering number for a semi-metric space $(T, d)$ and $Q$
is any probability measure.

Assumption (A1) is an ordinary smoothness condition; we need fourth derivatives since we need to bound the bias of second derivatives. The assumption (A2) is to make sure the collection of all local modes can be represented as finite collection of manifolds. (A3) is a sharpness requirement for all critical points (local modes and minimums); and it excludes the case that the modal manifolds bifurcate or merge, i.e., it excludes cases such as the right panel of Figure 2. Similar conditions appear in Romano (1988); Chen et al. (2014b)

for estimating density modes. Assumption (K1) is assumed for the kernel density estimator to have the usual rates for its bias and variance. (K2) is for the uniform bounds on the kernel density estimator; this condition can be found in

Giné and Guillou (2002); Einmahl et al. (2005); Chen et al. (2014c). We study three types of error metrics for regression modes: pointwise, uniform, and mean integrated squared errors. We defer all proofs to LABEL:suppA.

First we study the pointwise case. Recall that ${ˆ p}_{n}$ is the KDE in (4) of the joint density based on $n$ samples, under the kernel $K$ , and ${ˆ M}_{n} (x)$ is the estimated modal regression set in (5) at a point $x$ . Our pointwise analysis considers

Δ_{n} (x) = H a u s ({ˆ M}_{n} (x), M (x)),

the Hausdorff distance between ${ˆ M}_{n} (x)$ and $M (x)$ , at a point $x$ . For our first result, we define the quantities:

	$∥ {ˆ p}_{n} - p ∥_{\infty}^{(0)}$	$= sup x, y ∥ ˆ p (x, y) - p (x, y) ∥$
	$∥ {ˆ p}_{n} - p ∥_{\infty}^{(1)}$	$= sup x, y ∥ {ˆ p}_{y} (x, y) - p_{y} (x, y) ∥$
	$∥ {ˆ p}_{n} - p ∥_{\infty}^{(2)}$	$= sup x, y ∥ {ˆ p}_{y y} (x, y) - p_{y y} (x, y) ∥$
	$∥ {ˆ p}_{n} - p ∥_{\infty, 2}^{*}$	$= max {∥ {ˆ p}_{n} - p ∥_{\infty}^{(0)}, ∥ {ˆ p}_{n} - p ∥_{\infty}^{(1)}, ∥ {ˆ p}_{n} - p ∥_{\infty}^{(2)}} .$

Theorem 3 (Pointwise error rate).

Assume (A1-3) and (K1-2). Define a stochastic process $A_{n} (x)$ by

A_{n} (x) = {\begin{matrix} \frac{1}{Δ_{n} (x)} ∣ ∣ Δ_{n} (x) - {max}_{z \in M (x)} {| p_{y y}^{- 1} (x, z) | | {ˆ p}_{y, n} (x, z) |} ∣ ∣ & if Δ_{n} (x) > 0 0 & if Δ_{n} (x) = 0 . \end{matrix}

Then, when

∥ {ˆ p}_{n} - p ∥_{\infty, 2}^{*} = max {∥ {ˆ p}_{n} - p ∥_{\infty}^{(0)}, ∥ {ˆ p}_{n} - p ∥_{\infty}^{(1)}, ∥ {ˆ p}_{n} - p ∥_{\infty}^{(2)}}

is sufficiently small, we have

sup x \in D A_{n} (x) = O_{P} (∥ {ˆ p}_{n} - p ∥_{\infty, 2}^{*}) .

Moreover, at any fixed $x \in D$ , when $\frac{n h^{d + 5}}{log n} \to \infty$ and $h \to 0$ ,

Δ_{n} (x) = O (h^{2}) + O_{P} (\sqrt{\frac{1}{n h^{d + 3}}}) .

The proof is in the supplementary material (Chen et al., 2015). This shows that if the curvature of the joint density function along $y$ is bounded away from $0$ , then the error can be approximated by the error of ${ˆ p}_{y, n} (x, z)$ after scaling. The rate of convergence follows from the fact that ${ˆ p}_{y, n} (x, z)$ is converging to $0$ at the same rate. Note that as $z$ is a conditional mode, the partial derivative of the true density is $0$ . We defined $A_{n} (x)$ as above since $Δ_{n} (x) = 0$ implies ${max}_{z \in M (x)} | {ˆ p}_{y, n} (x, z) | = 0$ , so that the ratio would be ill-defined if $Δ_{n} (x) = 0$ . Also, the constraints on $h$ in the second assertion ( $\frac{n h^{d + 5}}{log n} \to \infty$ and $h \to 0$ ) are to ensure $∥ {ˆ p}_{n} - p ∥_{\infty, 2}^{*} = o_{P} (1)$ .

For our next result, we define the uniform error

Δ_{n} = sup x \in D Δ_{n} (x) = sup x \in D H a u s ({ˆ M}_{n} (x), M (x)) .

This is an $ℓ_{\infty}$ type error for estimating regression modes (and is also closely linked to confidence sets; see Section 5).

Theorem 4 (Uniform error rate).

Assume (A1-3) and (K1-2). Then as $\frac{n h^{d + 5}}{log n} \to \infty$ and $h \to 0$ ,

Δ_{n} = O (h^{2}) + O_{P} (\sqrt{\frac{log n}{n h^{d + 3}}}) .

The proof is in the supplementary material (Chen et al., 2015). Compared to the pointwise error rate in Theorem 3, we have an additional $\sqrt{log n}$ factor in the second term. One can view this as the price we need to pay for getting an uniform bound over all points. See Giné and Guillou (2002); Einmahl et al. (2005) for similar findings in density estimation.

The last error metric we consider is the mean integrated squared error (MISE), defined as

M I S E ({ˆ M}_{n}) = E (\int_{x \in D} Δ_{n}^{2} (x) d x) .

Note that the MISE is a nonrandom quantity, unlike first two error metrics considered.

Theorem 5 (MISE rate).

Assume (A1-3) and (K1-2). Then as $\frac{n h^{d + 5}}{log n} \to \infty$ and $h \to 0$ ,

M I S E ({ˆ M}_{n}) = O (h^{2}) + O (\sqrt{\frac{1}{n h^{d + 3}}}) .

The proof is in the supplementary material (Chen et al., 2015). If we instead focus on estimating the regression modes of the smoothed joint density $˜ p (x, y) = E ({ˆ p}_{n} (x, y))$ , then we obtain much faster convergence rates. Let $˜ M (x) = E ({ˆ M}_{n} (x))$ be the smoothed regression modes at $x \in D$ . Analogously define

	${˜ Δ}_{n} (x)$	$= H a u s ({ˆ M}_{n} (x), ˜ M (x))$
	${˜ Δ}_{n}$	$= sup x \in D {˜ Δ}_{n} (x)$
	$˜ M I S E ({ˆ M}_{n})$	$= E (\int_{x \in D} {˜ Δ}_{n}^{2} (x) d x) .$

Corollary 6 (Error rates for smoothed conditional modes).

Assume (A1-3) and (K1-2). Then as $\frac{n h^{d + 5}}{log n} \to \infty$ and $h \to 0$ ,

	$\sqrt{n h^{d + 3}} sup x \in D ∣ ∣ {˜ Δ}_{n} (x)$	$- max z \in ˜ M (x) {{˜ p}_{y y}^{- 1} (x, z) {ˆ p}_{y, n} (x, z)} ∣ ∣ = O_{P} (ϵ_{n, 2})$
	${˜ Δ}_{n} (x)$	$= O_{P} (\sqrt{\frac{1}{n h^{d + 3}}})$
	${˜ Δ}_{n}$	$= O_{P} (\sqrt{\frac{log n}{n h^{d + 3}}})$
	$˜ M I S E ({ˆ M}_{n})$	$= O (\sqrt{\frac{1}{n h^{d + 3}}}),$

where $ϵ_{n, 2} = {sup}_{x, y} | {ˆ p}_{y y, n} (x, y) - {˜ p}_{y y} (x, y) | = {sup}_{x, y} | {ˆ p}_{y y, n} (x, y) - E ({ˆ p}_{y y, n} (x, y)) |$ .

5 Confidence sets

In an idealized setting, we could define a confidence set at $x$ by

{ˆ C}_{n}^{0} (x) = {ˆ M}_{n} (x) \oplus δ_{n, 1 - α} (x),

where

P (Δ_{n} (x) > δ_{n, 1 - α} (x)) = α .

By construction, we have $P (M (x) \in {ˆ C}_{n}^{0} (x)) = 1 - α$ . Of course, the distribution of $Δ_{n} (x)$ is unknown, but we can use the bootstrap (Efron, 1979) to estimate $δ_{n, 1 - α} (x)$ .

Given the observed data samples $(X_{1}, Y_{1}), \dots, (X_{n}, Y_{n})$ , we denote a bootstrap sample as $(X_{1}^{*}, Y_{1}^{*}), \dots, (X_{n}^{*}, Y_{n}^{*})$ . Let ${ˆ M}_{n}^{*} (x)$ be the estimated regression modes based on this bootstrap sample, and

{ˆ Δ}_{n}^{*} (x) = H a u s ({ˆ M}_{n}^{*} (x), {ˆ M}_{n} (x)) .

We repeat the bootstrap sampling $B$ times to get ${ˆ Δ}_{1, n}^{*} (x), \dots, {ˆ Δ}_{B, n}^{*} (x)$ . Define ${ˆ δ}_{n, 1 - α} (x)$ by

\frac{1}{B} B \sum j = 1 I ({ˆ Δ}_{j, n}^{*} (x) > {ˆ δ}_{n, 1 - α} (x)) = α .

Our confidence set for $M (x)$ is then given by

{ˆ C}_{n} (x) = {ˆ M}_{n} (x) \oplus {ˆ δ}_{n, 1 - α} (x) .

Note that this is a pointwise confidence set, at $x \in D$ .

Alternatively, we can use $Δ_{n} = {sup}_{x \in D} Δ_{n} (x)$ to build a uniform confidence set. Define $δ_{n, 1 - α}$ by

P (M (x) \subseteq {ˆ M}_{n} (x) \oplus δ_{n, 1 - α}, \forall x \in D) = 1 - α .

As above, we can use bootstrap sampling to form an estimate ${ˆ δ}_{n, 1 - α}$

, based on the quantiles of the bootstrapped uniform error metric

{ˆ Δ}_{n}^{*} = sup x \in D H a u s ({ˆ M}_{n}^{*} (x), {ˆ M}_{n} (x)) .

Our uniform confidence set is then

{ˆ C}_{n} = {(x, y) : x \in D, y \in {ˆ M}_{n} (x) \oplus {ˆ δ}_{n, 1 - α}} .

In practice, there are many possible flavors of the bootstrap that are applicable here. This includes the ordinary (nonparametric) bootstrap, the smoothed bootstrap and the residual bootstrap. See Figure 4 for an example with the ordinary bootstrap.

We focus on the asymptotic coverage of uniform confidence sets built with the ordinary bootstrap. We consider coverage of the smoothed regression mode set $˜ M (x)$ (to avoid issues of bias), and we employ tools developed in Chernozhukov et al. (2014a, b); Chen et al. (2014c).

Figure 4: An example with pointwise (left) and uniform (right) confidence sets. The significance level is 90%.

Consider a function space $F$ defined as

F = {(u, v) \mapsto f_{x, y} (u, v) : f_{x, y} (u, v) = {˜ p}_{y y}^{- 1} (x, y) \cdot K (\frac{∥ x - u ∥}{h}) K^{(1)} (\frac{y - v}{h}), x \in D, y \in ˜ M (x)},

(15)

and let $B$ be a Gaussian process defined on $F$ such that

Cov (B (f_{1}), B (f_{2})) = E (f_{1} (X_{i}, Y_{i}) f_{2} (X_{i}, Y_{i})) - E (f_{1} (X_{i}, Y_{i})) E (f_{2} (X_{i}, Y_{i})),

(16)

for all $f_{1}, f_{2} \in F$ .

Theorem 7.

Assume (A1-3) and (K1-2). Define the random variable

B = \frac{1}{\sqrt{h^{d + 3}}} {sup}_{f \in F} | B (f) |

. Then as

\frac{n h^{d + 5}}{log n} \to \infty

h \to 0

sup t \geq 0 ∣ ∣ P (\sqrt{n h^{d + 3}} {˜ Δ}_{n} < t) - P (B < t) ∣ ∣ = O ⎛ ⎝ {(\frac{{log}^{7} n}{n h^{d + 3}})}^{1 / 8} ⎞ ⎠ .

The proof is in the supplementary material (Chen et al., 2015). This theorem shows that the smoothed uniform discrepancy ${˜ Δ}_{n}$ is distributed asymptotically as the supremum of a Gaussian process. In fact, it can be shown that the two random variables ${˜ Δ}_{n}$ and $B$ are coupled by

∣ ∣ \sqrt{n h^{d + 3}} {˜ Δ}_{n} - B ∣ ∣ = O_{P} ⎛ ⎝ {(\frac{{log}^{7} n}{n h^{d + 3}})}^{1 / 8} ⎞ ⎠ .

Now we turn to the limiting behavior for the bootstrap estimate. Let $D_{n} = {(X_{1}, Y_{1}), \dots, (X_{n}, Y_{n})}$ be the observed data and denote the bootstrap estimate by

{ˆ Δ}_{n}^{*} = sup x \in D H a u s ({ˆ M}_{n}^{*} (x), {ˆ M}_{n} (x))

where ${ˆ M}_{n}^{*} (x)$ is the the bootstrap regression mode set at $x$ .

Theorem 8 (Bootstrap consistency).

Assume conditions (A1-3) and (K1-2). Also assume that $\frac{n h^{d + 5}}{log n} \to \infty, h \to 0$ . Define $B = \frac{1}{\sqrt{h^{d + 3}}} {sup}_{f \in F} | B (f) |$ . There exists $X_{n}$ such that $P (X_{n}) \geq 1 - O (\frac{1}{n})$ and, for all $D_{n} \in X_{n}$ ,

sup t \geq 0 ∣ ∣ P (\sqrt{n h^{d + 3}} {ˆ Δ}_{n}^{*} < t ∣ ∣ D_{n}) - P (B < t) ∣ ∣ = O_{P} ⎛ ⎝ {(\frac{{log}^{7} n}{n h^{d + 3}})}^{1 / 8} ⎞ ⎠ .

The proof is in the supplementary material (Chen et al., 2015). Theorem 8 shows that the limiting distribution for the bootstrap estimate ${ˆ Δ}_{n}^{*}$ is the same as the limiting distribution of ${˜ Δ}_{n}$ (recall Theorem 7) with high probability. Using Theorems 7 and 8, we conclude the following.

Corollary 9 (Uniform confidence sets).

Assume conditions (A1-3) and (K1-2). Then as $\frac{n h^{d + 5}}{log n} \to \infty$ and $h \to 0$ ,

P (˜ M (x) \subseteq {ˆ M}_{n} (x) \oplus {ˆ δ}_{n, 1 - α}, \forall x \in D) = 1 - α + O ⎛ ⎝ {(\frac{{log}^{7} n}{n h^{d + 3}})}^{1 / 8} ⎞ ⎠ .

6 Prediction sets

Modal regression can be also used to construct prediction sets. Define

	$ϵ_{1 - α} (x)$	$= inf {ϵ \geq 0 : P (d (Y, M (X)) > ϵ \| X = x) \leq α}$
	$ϵ_{1 - α}$	$= inf {ϵ \geq 0 : P (d (Y, M (X)) > ϵ) \leq α} .$

Recall that $d (x, A) = {inf}_{y \in A} | x - y |$ for a point $x$ and a set $A$ . Then

	$P_{1 - α} (x)$	$= M (x) \oplus ϵ_{1 - α} (x) \subseteq R$
	$P_{1 - α}$	$= {(x, y) : x \in D, y \in M (x) \oplus ϵ_{1 - α}} \subseteq D \times R$

are pointwise and uniform prediction sets, respectively, at the population level, because

	$P (Y \in P_{1 - α} (x) \| X = x)$	$\geq 1 - α$
	$P (Y \in P_{1 - α})$	$\geq 1 - α .$

At the sample level, we use a KDE of the conditional density ${ˆ p}_{n} (y | x) = {ˆ p}_{n} (x, y) / {ˆ p}_{n} (x)$ , and estimate $ϵ_{1 - α} (x)$ via

{ˆ ϵ}_{1 - α} (x) = inf {ϵ \geq 0 : \int_{{ˆ M}_{n} (x) \oplus ϵ} {ˆ p}_{n} (y | x) d y \geq 1 - α} .

An estimated pointwise prediction set is then

{ˆ P}_{1 - α} (x) = {ˆ M}_{n} (x) \oplus {ˆ ϵ}_{1 - α} (x) .

This has the proper pointwise coverage with respect to samples drawn according to ${ˆ p}_{n} (y | x)$ , so in an asymptotic regime in which ${ˆ p}_{n} (y | x) \to p_{n} (y | x)$ , it will have the correct coverage with respect to the population distribution, as well.

Similarly, we can define

{ˆ ϵ}_{1 - α} = Q u a n t i l e ({d (Y_{i}, ˆ M_{n} (X_{i})) : i = 1, \dots, n}, 1 - α),

(17)

the $(1 - α)$ quantile of $d (Y_{i}, {ˆ M}_{n} (X_{i}))$ , $i = 1, \dots, n$ , and then the estimated uniform prediction set is

{ˆ P}_{1 - α} = {(x, y) : x \in D, y \in {ˆ M}_{n} (x) \oplus {ˆ ϵ}_{1 - α}} .

(18)

The estimated uniform prediction set has proper coverage with respect to the empirical distribution, and so certain conditions, it will have valid limiting population coverage.

6.1 Bandwidth selection

Prediction sets can be used to select the smoothing bandwidth of the underlying KDE, as we describe here. We focus on uniform prediction sets, and we will use a subscript $h$ throughout to denote the dependence on the smoothing bandwidth. From its definition in (18), we can see that the volume (Lebesgue measure) of the estimated uniform prediction set is

V o l ({ˆ P}_{1 - α, h}) = {ˆ ϵ}_{1 - α, h} \int_{x \in D} {ˆ K}_{h} (x) d x,

where ${ˆ K}_{h} (x)$ is the number of estimated local modes at $X = x$ , and ${ˆ ϵ}_{1 - α, h}$ is as defined in (17). Roughly speaking, when $h$ is small, ${ˆ ϵ}_{1 - α, h}$ is also small, but the number of estimated manifolds is large; on the other hand, when $h$ is large, ${ˆ ϵ}_{1 - α, h}$ is large, but the number of estimated manifolds is small. This is like the bias-variance trade-off: small $h$ corresponds to less bias ( ${ˆ ϵ}_{1 - α, h}$ ) but higher variance (number of estimated manifolds).

Our proposal is to select $h$ by

h^{*} = a r g m i n h \geq 0 V o l ({ˆ P}_{1 - α, h}) .

Figure 5 gives an example this rule when $α = 0.05$ , i.e., when minimizing the size of the estimated 95% uniform prediction set. Here we actually use cross-validation to obtain the size of the prediction set; namely, we use the training set to estimate the modal manifolds and then use the validation set to estimate the width of prediction set. This helps us to avoid overfitting. As can be seen, there is a clear trade-off in the size of the prediction set versus $h$ in the left plot. The optimal value $h^{*} = 0.07$ is marked by a vertical line, and the right plot displays the corresponding modal regression estimate and uniform prediction set on the data samples.

In the same plot, we also display a local regression estimate and its corresponding 95% uniform prediction set. We can see that the prediction set from the local regression method is much larger than that from modal regression. (To even the comparison, the bandwidth for the local linear smoother was also chosen to minimize the size of the prediction set.) This illustrates a major strength of the modal regression method: because it is not constrained to modeling conditional mean structure, it can produce smaller prediction sets than the usual regression methods when the conditional mean fails to capture the main structure in the data. We investigate this claim theoretically, next.

Figure 5: An example of bandwidth selection based on the size of the prediction sets.

6.2 Theory on the size of prediction sets

We will show that, at the population level, prediction sets from modal regression can be smaller than those based on the underlying regression function $μ (x) = E (Y | X = x)$ . Defining

	$η_{1 - α} (x)$	$= inf {η \geq 0 : P (d (Y, μ (X)) > η \| X = x) \leq α}$
	$η_{1 - α}$

pointwise and uniform prediction sets based on the regression function are

	$R_{1 - α} (x)$	$= μ (x) \oplus η_{1 - α} (x) \subseteq R$
	$R_{1 - α}$	$= {(x, μ (x) \oplus η_{1 - α}) : x \in D} \subseteq D \times R,$

respectively.

For a pointwise prediction set $A (x)$ , we write $l e n g t h (A (x))$ for its Lebesgue measure on $R$ ; note that in the case of modal regression, this is somewhat of an abuse of notation because the Lebesgue measure of $A (x)$ can be a sum of interval lengths. For a uniform prediction set $A$ , we write $V o l (A)$ for its Lebesgue measure on $D \times R$ .

We consider the following assumption.

The conditional density satisfies

$p (y | x) = K (x) \sum j = 1 π_{j} (x) ϕ (y; μ_{j} (x), σ_{j}^{2} (x))$

with $μ_{1} (x) < μ_{2} (x) < \dots < μ_{K (x)} (x)$ by convention, and $ϕ (\cdot; μ, σ^{2})$ denoting the Gaussian density function with mean $μ$ and variance $σ^{2}$ .

The assumption that the conditional density can be written as a mixture of Gaussians is only used for the next result. It is important to note that this is an assumption made about the population density, and does not reflect modeling choices made in the sample. Indeed, recall, we are comparing prediction sets based on the modal set $M (x)$ and the regression function $μ (x)$ , both of which use true population information.

Before stating the result, we must define several quantities. Define the minimal separation between mixture centers

Δ_{min} (x) = min {| μ_{i} (x) - μ_{j} (x) | : i \neq j}

and

	$σ_{max}^{2} (x) = max j = 1, \dots, K (x) σ_{j}^{2} (x),$
	$π_{max} (x) = max j = 1, \dots, K (x) π_{j} (x), π_{min} (x) = min j = 1, \dots, K (x) π_{j} (x) .$

Also define

Δ_{min} = inf x \in D Δ_{min} (x), σ_{max}^{2} = sup x \in D σ_{max}^{2} (x),

and

π_{max} = sup x \in D π_{max} (x), π_{min} = inf x \in D π_{min} (x),

and

¯ ¯¯¯ ¯ K = \frac{\int_{x \in D} K (x) d x}{\int_{x \in D} d x}, K_{min} =

Nonparametric modal regression

Modal clustering asymptotics with applications to bandwidth selection

Modal-set estimation with an application to clustering

An implicit function learning approach for parametric modal regression

Modal clustering of matrix-variate data

Asymptotic Confidence Regions for Density Ridges

Optimal Sampling Density for Nonparametric Regression

Kernel Selection for Modal Linear Regression: Optimal Kernel and IRLS Algorithm

Code Repositories

ModalRegression

1 Introduction

2 Review of modal regression

3 Geometric properties

Lemma 1 (Derivative of modal functions).

Theorem 2 (Smoothness of the modal manifold collection).

4 Asymptotic error analysis

Theorem 3 (Pointwise error rate).

Theorem 4 (Uniform error rate).

Theorem 5 (MISE rate).

Corollary 6 (Error rates for smoothed conditional modes).

5 Confidence sets

Theorem 7.

Theorem 8 (Bootstrap consistency).

Corollary 9 (Uniform confidence sets).

6 Prediction sets

6.1 Bandwidth selection

6.2 Theory on the size of prediction sets

Nonparametric modal regression

Related Research

Modal clustering asymptotics with applications to bandwidth selection

Modal-set estimation with an application to clustering

An implicit function learning approach for parametric modal regression

Modal clustering of matrix-variate data

Asymptotic Confidence Regions for Density Ridges

Optimal Sampling Density for Nonparametric Regression

Kernel Selection for Modal Linear Regression: Optimal Kernel and IRLS Algorithm

Code Repositories

ModalRegression

1 Introduction

2 Review of modal regression

3 Geometric properties

Lemma 1 (Derivative of modal functions).

Theorem 2 (Smoothness of the modal manifold collection).

4 Asymptotic error analysis

Theorem 3 (Pointwise error rate).

Theorem 4 (Uniform error rate).

Theorem 5 (MISE rate).

Corollary 6 (Error rates for smoothed conditional modes).

5 Confidence sets

Theorem 7.

Theorem 8 (Bootstrap consistency).

Corollary 9 (Uniform confidence sets).

6 Prediction sets

6.1 Bandwidth selection

6.2 Theory on the size of prediction sets