\section{Introduction} \label{sec:1}
A statistical depth function provides an ordering of the points such that the points close to the centre of the distribution have higher depth values than those that are away from the centre. Several depth functions are defined in the literature, e.g., the halfspace depth \citep{tukey1975mathematics,donoho1992breakdown}, the simplicial depth \citep{liu1990notion}, the spatial depth \citep{vardi2000multivariate,serfling2002depth}, the projection depth \citep{zuo2000general,zuo2000structural}, and many others.
Depth functions have been used earlier for various purposes like detecting outliers \citep{chen2009outlier}, clustering \citep{jornsten2004clustering} and classification \citep{ghosh2005maximum,dutta2012robust,li2012dd}.
We develop nonparametric regression methods using statistical depth functions. The response in our setup is multivariate, while the covariate is a random element in a separable metric space, which includes the cases of finite dimensional as well as infinite dimensional or function-valued covariates.
The traditional nonparametric regression methods estimate the conditional mean or median (see, e.g., \cite{hardle1990applied}, \cite{fan1996local}, \cite{wand1994kernel}, etc.), and provide information about only the center of the conditional distribution of the response given the covariate. These methods are not useful if one is interested in other features of the conditional distribution, like the conditional spread and heteroscedasticity in the sample. The methods we develop will provide simultaneous information about the center as well as other parts of the conditional distribution, and help us construct measures of conditional spread of the response.

Our methods based on statistical depth are inspired from Tukey's boxplot \citep{tukey1975mathematics,tukey1977exploratory}. The boxplot for a univariate random variable provides information about the center of the distribution through the median as well as the spread and the skewness through the first and third quartiles. Let us consider the following two examples, which demonstrate the inadequacy of regression procedures that consider only the center of the conditional distribution.

\emph{Example 1}: In our first example, we consider the BudgetItaly Data, which is available in the Ecdat package in R. This dataset contains the following information on budget shares of 1729 Italian households for the period 1973 to 1992: the shares of food expenditure, housing and fuels expenditure and miscellaneous expenditure in the total budget along with the corresponding food price, housing and fuels price and miscellaneous price. We consider the variable pair of food share and housing and fuels share as the response, and the variable pair of food price and housing and fuels price as the covariate.

To carry out nonparametric regression, we construct a neighborhood consisting of 50 nearest covariate values for each covariate value, based on the euclidiean metric. The collection of the corresponding response values for a covariate value represents the conditional distribution of the response given that covariate value. We call these response values the local response values for the covariate. We choose 4 covariate values, and plot their local response values and the marginal boxplots for the local response values in \autoref{fig:paper3ex1}.
The selected covariate values are mentioned over the corresponding columns in the plot.
\begin{figure}
	\centering
	\includegraphics[width=1\linewidth]{"paper3ex1"}
	\caption{The local response values (1st row), the local boxplots for the 1st (2nd row) and the 2nd coordinates (3rd row) of the response variable in the BudgetItaly Data.}
	\label{fig:paper3ex1}
\end{figure}
The local boxplots, which are the boxplots for the first and the second coordinates of the local response values, estimate the conditional boxplots of the two coordinates of the response given the covariate. They clearly reflect the varying conditional spread and skewness of the response over the covariate values, which could not be detected from mean or median regression. But as indicated by the scatter plots of the local response values, the two coordinates are correlated. Consequently, it is desirable to investigate the features of the joint distribution of the two coordinates given the covariate, which cannot be done using the local marginal boxplots.

\emph{Example 2}: Our second example deals with the Tecator Data, which is available in the caret package in R. This dataset contains the percentage values of moisture, fat and protein contents of 215 meat samples along with their absorbance spectra in the wavelength range 850--1050 nm. The moisture, the fat and the protein contents were measured by analytical chemistry, while a Tecator Infratec Food and Feed Analyzer was used to record the absorbance spectrum. Being able to predict the nutritional contents of a meat sample from its absorbance spectra is economically beneficial since obtaining the spectra is relatively cheaper.
We consider the pair of fat and protein contents as the response, and the curve of absorbance spectra as the covariate. So, the response is bivariate and the covariate is functional. The covariate is considered to be a random element in the $ L_2 $ space. We do the same analysis as in \emph{Example 1}, with the euclidean metric replaced by the $ L_2 $ metric and changing the nearest neighbor size to 30. We plot the results, along with the chosen covariate curves, in \autoref{fig:paper3ex2}.
\begin{figure}
	\centering
	\includegraphics[width=1\linewidth]{"paper3ex2"}
	\caption{The selected covariate curves (1st row), the local response values (2nd row), the local boxplots for the 1st (3rd row) and the 2nd coordinates (4th row) of the response variable in the Tecator Data.}
	\label{fig:paper3ex2}
\end{figure}
Here again, the local boxplots reflect the varying conditional spread and skewness of the response over the covariate curves, and the scatter plots reflect the dependence between the two coordinates of the local response values. This again implies the importance of investigating the features of the bivariate conditional distribution of the response given the covariate curves.

The preceding two examples demonstrate the usefulness of the boxplots in yielding information about univariate data, as well as their shortcomings when the data is not univariate.
We develop depth based regression methods, which is applicable for multivariate responses.
Earlier attempts to generalize boxplot using statistical depth functions in unconditional setup can be found in \cite{rousseeuw1999bagplot} and \cite{liu1999multivariate} for bivariate data and in \cite{sun2011functional} for functional data.
For a univariate distribution, the centre of the boxplot, which is the median, is the point having the largest depth value, and the box contains 50\% of the observations which have higher depth values than the points outside the box. These basic features of the boxplot are utilized to generalize boxplots for multivariate conditional distributions in the following section. Such a boxplot based regression has fundamental connections with the well-known quantile regression \citep{koenker2005quantile}.

\section{Conditional depth and central regions} \label{subsec:2.1}
Let $ \mathbf{Y} $ be a random vector in $ \mathbb{R}^p $ and $ \mathbf{X} $ be a random element in a complete separable metric space $ ( \mathcal{C}, d ) $. Let $ \mu( \cdot \midil \mathbf{z} ) $ be the conditional probability measure of $ \mathbf{Y} $ given $ \mathbf{X} = \mathbf{z} $, and $ \rho( \cdot \midil \mathbf{z} ) $ be a conditional depth function on the response space $ \mathbb{R}^p $ related to $ \mu( \cdot \midil \mathbf{z} ) $. Values of the depth functions usually lie between 0 and 1 and consequently it is uniformly bounded. Let $ \mathbf{x} $ be a fixed element in $ \mathcal{C} $. Define $ D( \alpha \midil \mathbf{x} ) = \{ \mathbf{y} \in \mathbb{R}^p \midil \rho( \mathbf{y} \midil \mathbf{x} ) \ge \alpha \} $ for $ \alpha \in \mathbb{R} $. $ D( \alpha \midil \mathbf{x} ) $ is called the conditional $ \alpha $-trimmed region of $ \mathbf{Y} $ given $ \mathbf{X} = \mathbf{x} $ corresponding to the conditional depth $ \rho( \cdot \midil \mathbf{x} ) $ (cf. \cite{zuo2000general,zuo2000structural}).
For $ 0 \le r < 1 $, let $ \alpha( r ) = \sup\{ \alpha \midil \mu( D( \alpha \midil \mathbf{x} ) \midil \mathbf{x} ) \ge r \} $, which is finite when $ \rho( \mathbf{y} \midil \mathbf{x} ) $ is uniformly bounded over $ \mathbf{y} $. We define the set $ D( \alpha( r ) \midil \mathbf{x} ) $ as the conditional $ 100 r $\% central region of $ \mathbf{Y} $ given $ \mathbf{X} = \mathbf{x} $ with respect to the conditional depth $ \rho( \cdot \midil \mathbf{x} ) $.
Clearly, $ \mu( D( \alpha( r ) \midil \mathbf{x} ) \midil \mathbf{x} ) \ge r $, and $ \mu( D( \alpha( r ) \midil \mathbf{x} ) \midil \mathbf{x} ) = r $ whenever $ \mu( \{ \mathbf{y} \in \mathbb{R}^p \midil \rho( \mathbf{y} \midil \mathbf{x} ) = \alpha( r ) \} \midil \mathbf{x} ) = 0 $. Further, for any $ \mathbf{y}_1 \in D( \alpha( r ) \midil \mathbf{x} ) $ and $ \mathbf{y}_2 \not\in D( \alpha( r ) \midil \mathbf{x} ) $, $ \rho( \mathbf{y}_1 \midil \mathbf{x} ) \ge \rho( \mathbf{y}_2 \midil \mathbf{x} ) $.
When the response is univariate, the conditional $ 50 $\% central region corresponds to the box in the conditional boxplot of the response.

To estimate $ \rho( \cdot \midil \mathbf{x} ) $, $ D( \alpha \midil \mathbf{x} ) $ and $ D( \alpha( r ) \midil \mathbf{x} ) $ based on a random sample $ ( \mathbf{X}_1, \mathbf{Y}_1 ), \allowbreak ..., ( \mathbf{X}_n, \mathbf{Y}_n ) $, we adopt a nonparametric regression procedure. Let $ W_{i,n}( \mathbf{x} ) $ be the weight on the observation pair $ ( \mathbf{X}_i, \mathbf{Y}_i ) $, $ i = 1, ..., n $, $ W_{i,n}( \mathbf{x} ) \ge 0 $ for each $ i $ and $ \sum_{i=1}^{n} W_{i,n}( \mathbf{x} ) = 1 $. The sample conditional probability measure of $ \mathbf{Y} $ given $ \mathbf{X} = \mathbf{x} $ is defined as
\begin{align*}
\mu_n( B \midil \mathbf{x} ) = \sum_{i=1}^{n} I( \mathbf{Y}_i \in B ) W_{i,n}( \mathbf{x} ) ,
\end{align*}
where $ B $ is any Borel set. The conditional sample depth function $ \rho_n( \cdot \midil \mathbf{x} ) $ is defined based on $ \mu_n( \cdot \midil \mathbf{x} ) $ in the same way as the conditional population depth function $ \rho( \cdot \midil \mathbf{x} ) $ is defined based on $ \mu( \cdot \midil \mathbf{x} ) $. The conditional sample $ \alpha $-trimmed region of $ \mathbf{Y} $ given $ \mathbf{X} = \mathbf{x} $ is defined as $ D_n( \alpha \midil \mathbf{x} ) = \{ \mathbf{y} \in \mathbb{R}^p \midil \rho_n( \mathbf{y} \midil \mathbf{x} ) \ge \alpha \} $. The conditional sample $ 100 r $\% central region of $ \mathbf{Y} $ given $ \mathbf{X} = \mathbf{x} $ is $ D_n( \alpha_n( r ) \midil \mathbf{x} ) $, where $ \alpha_n( r ) = \sup\{ \alpha \midil \mu_n( D_n( \alpha \midil \mathbf{x} ) \midil \mathbf{x} ) \ge r \} $.
In an unconditional setup, a version of the central region was defined in \citet[p.~788]{liu1999multivariate}, and \cite{chowdhury2016nonparametric} defined conditional maximal depth sets based on the spatial depth. However, neither \cite{liu1999multivariate} nor \cite{chowdhury2016nonparametric} investigated theoretical properties like asymptotic consistency of the sample central regions and maximal depth sets.

We now describe the conditional versions of several well-known depth functions. The conditional halfspace depth \citep{tukey1975mathematics,donoho1992breakdown} of $ \mathbf{Y} $ given $ \mathbf{X} = \mathbf{x} $ is defined as $ \rho( \mathbf{y} \midil \mathbf{x} ) = \inf\{ \mu( \{ \mathbf{v} \in \mathbb{R}^p \midil \mathbf{u}^t \mathbf{v} \ge \mathbf{u}^t \mathbf{y} \} \midil \mathbf{x} ) \midil \mathbf{u} \in \mathbb{R}^p \} $. Its estimate is $ \rho_n( \mathbf{y} \midil \mathbf{x} ) = \inf\{ \mu_n( \{ \mathbf{v} \in \mathbb{R}^p \midil \mathbf{u}^t \mathbf{v} \ge \mathbf{u}^t \mathbf{y} \} \midil \mathbf{x} ) \midil \mathbf{u} \in \mathbb{R}^p \} $.
The conditional Spatial Depth \citep{vardi2000multivariate,serfling2002depth} is defined as $ \rho( \mathbf{y} \midil \mathbf{x} ) = 1 - \| E[ \| \mathbf{y} - \mathbf{Y} \|^{-1} ( \mathbf{y} - \mathbf{Y} ) \midil \mathbf{X} = \mathbf{x} ] \| $. The estimate $ \rho_n( \mathbf{y} \midil \mathbf{x} ) $ is given by
\begin{align*}
\rho_n( \mathbf{y} \midil \mathbf{x} ) = 1 - \left\| \sum_{i=1}^{n} \| \mathbf{y} - \mathbf{Y}_i \|^{-1} ( \mathbf{y} - \mathbf{Y}_i ) W_{i,n}( \mathbf{x}  ) \right\| .
\end{align*}
The conditional Projection Depth \citep{zuo2000general,zuo2000structural} is defined as
\begin{align*}
\rho( \mathbf{y} \midil \mathbf{x} ) = \left[ 1 + \sup_{\| \mathbf{u} \| = 1} \frac{| \mathbf{u}^t \mathbf{y} - M( \mathbf{u}^t \mathbf{Y} \midil \mathbf{x} ) |}{M( | \mathbf{u}^t \mathbf{Y} - M( \mathbf{u}^t \mathbf{Y} \midil \mathbf{x} ) | \midil \mathbf{x} )} \right]^{-1} ,
\end{align*}
where $ M( \mathbf{u}^t \mathbf{Y} \midil \mathbf{x} ) $ and $ M( | \mathbf{u}^t \mathbf{Y} - M( \mathbf{u}^t \mathbf{Y} \midil \mathbf{x} ) | \midil \mathbf{x} ) $ are the conditional medians of $ \mathbf{u}^t \mathbf{Y} $ and $ | \mathbf{u}^t \mathbf{Y} - M( \mathbf{u}^t \mathbf{Y} \midil \mathbf{x} ) | $ given $ \mathbf{X} = \mathbf{x} $, respectively.
The conditional sample projection depth is
\begin{align*}
\rho_n( \mathbf{y} \midil \mathbf{x} ) = \left[ 1 + \sup_{\| \mathbf{u} \| = 1} \frac{| \mathbf{u}^t \mathbf{y} - M_n( \mathbf{u}^t \mathbf{Y} \midil \mathbf{x} ) |}{M_n( | \mathbf{u}^t \mathbf{Y} - M_n( \mathbf{u}^t \mathbf{Y} \midil \mathbf{x} ) | \midil \mathbf{x} )} \right]^{-1} ,
\end{align*}
where $ M_n( \mathbf{u}^t \mathbf{Y} \midil \mathbf{x} ) $ and $ M_n( | \mathbf{u}^t \mathbf{Y} - M_n( \mathbf{u}^t \mathbf{Y} \midil \mathbf{x} ) | \midil \mathbf{x} ) $ are the sample analogues of $ M( \mathbf{u}^t \mathbf{Y} \midil \mathbf{x} ) $ and $ M( | \mathbf{u}^t \mathbf{Y} - M( \mathbf{u}^t \mathbf{Y} \midil \mathbf{x} ) | \midil \mathbf{x} ) $, respectively.
Let $ \mathbf{Y}^{(1)}, ..., \mathbf{Y}^{(p+1)} $ be independent copies of $ \mathbf{Y} \in \mathbb{R}^p $, and let $ \Delta\left( \mathbf{Y}^{(1)}, ..., \mathbf{Y}^{(p+1)} \right) $ denote the $ p $-dimensional simplex formed from $ \mathbf{Y}^{(1)}, ..., \mathbf{Y}^{(p+1)} $. The conditional simplicial depth \citep{liu1990notion} is defined as
\begin{align*}
	\rho( \mathbf{y} \midil \mathbf{x} ) = P\left[ \mathbf{y} \in \Delta\left( \mathbf{Y}^{(1)}, ..., \mathbf{Y}^{(p+1)} \right) \mid \mathbf{X} = \mathbf{x} \right] .
\end{align*}
The estimate $ \rho_n( \mathbf{y} \midil \mathbf{x} ) $ is defined as
\begin{align*}
	\rho_n( \mathbf{y} \midil \mathbf{x} ) = \frac{\sum_{1 \le i_1 < ... < i_{p+1} \le n} I\left( \mathbf{y} \in \Delta\left( \mathbf{Y}_{i_1}, ..., \mathbf{Y}_{i_{p+1}} \right) \right) W_{i_1,n}( \mathbf{x} ) ... W_{i_{p+1},n}( \mathbf{x} )}{\sum_{1 \le i_1 < ... < i_{p+1} \le n} W_{i_1,n}( \mathbf{x} ) ... W_{i_{p+1},n}( \mathbf{x} )} .
\end{align*}

The weights $ \{  W_{i,n}( \mathbf{x} ) \} $ are constructed based on the covariate values $ \mathbf{X}_1, ..., \mathbf{X}_n $. There are several methods of selecting such weights. In the kernel regression method, we choose a kernel function $ K( \cdot ) $ and a bandwidth $ h $, and the weight $ W_{i,n}( \mathbf{x} ) $ is
\begin{align*}
W_{i,n}( \mathbf{x} ) = \frac{K( h^{-1} d( \mathbf{x}, \mathbf{X}_i ) )}{\sum_{i=1}^{n} K( h^{-1} d( \mathbf{x}, \mathbf{X}_i ) )} ,
\end{align*}
where $ d( \cdot, \cdot ) $ is the metric in the covariate space. This leads to a Nadaraya-Watson type kernel estimate \citep{nadaraya1964estimating,watson1964smooth}. In the nearest neighbor method, we choose a positive integer $ k $ for the number of nearest neighbors to be considered, and define
\begin{align*}
h( \mathbf{x}, k, n ) = \min\left\{ h \mid \sum_{i=1}^{n} I( d( \mathbf{x}, \mathbf{X}_i ) \le h ) \ge k \right\} .
\end{align*}
The weight $ W_{i,n}( \mathbf{x} ) $ in this case is
\begin{align*}
W_{i,n}( \mathbf{x} ) = \frac{I( d( \mathbf{x}, \mathbf{X}_i ) \le h( \mathbf{x}, k, n ) )}{\sum_{i=1}^{n} I( d( \mathbf{x}, \mathbf{X}_i ) \le h( \mathbf{x}, k, n ) )} .
\end{align*}

\subsection{Median and trimmed mean regression} \label{subsec:2.2}
The depth-based conditional median $ \mathbf{m}( \mathbf{x} ) $ of $ \mathbf{Y} $ given $ \mathbf{X} = \mathbf{x} $ with respect to the conditional depth $ \rho( \cdot \midil \mathbf{x} ) $ is a point such that $ \rho( \mathbf{m}( \mathbf{x} ) \midil \mathbf{x} ) \ge \rho( \mathbf{y} \midil \mathbf{x} ) $ for every $ \mathbf{y} $. Note that $ \mathbf{m}( \mathbf{x} ) $ may not be unique. For a univariate response and the usual depth functions like the conditional halfspace depth, the point $ \mathbf{m}( \mathbf{x} ) $ becomes the usual conditional median, which can be viewed as the statistical center of the conditional boxplot, and the set $ D( \alpha( r ) \midil \mathbf{x} ) $ becomes the conditional interquartile interval for $ r = 0.5 $, which corresponds to the box in the conditional boxplot as already noted. $ \mathbf{m}( \mathbf{x} ) $ along with the conditional central region $ D( \alpha( r ) \midil \mathbf{x} ) $ yields information about the centre and the spread of the conditional distribution of the response.
The sample conditional median $ \mathbf{m}_n( \mathbf{x} ) $ is a point such that $ \rho_n( \mathbf{m}_n( \mathbf{x} ) \midil \mathbf{x} ) \ge \rho_n( \mathbf{y} \midil \mathbf{x} ) $ for every $ \mathbf{y} $.

Trimmed means based on depth functions in an unconditional setup were investigated earlier in \cite{donoho1992breakdown}, \cite{liu1999multivariate}, \cite{zuo2006multidimensional}, \cite{masse2009multivariate}, etc.
The conditional $ 100 r $\% trimmed mean $ \mathbf{m}( r \midil \mathbf{x} ) $ of $ \mathbf{Y} $ given $ \mathbf{X} = \mathbf{x} $ can be defined as $ \mathbf{m}( r \midil \mathbf{x} ) = [ \int \mathbf{y} I( \mathbf{y} \in D( \alpha( 1 - r ) \midil \mathbf{x} ) ) \mu( d \mathbf{y} \midil \mathbf{x} ) ] / \mu( D( \alpha( 1 - r ) \midil \mathbf{x} ) \midil \mathbf{x} ) $.
Unlike the conditional median, the conditional trimmed mean is always unique.
For a real valued response, the depth-based conditional trimmed mean coincides with the usual conditional trimmed mean.
The sample conditional $ 100 r $\% trimmed mean $ \mathbf{m}_n( r \midil \mathbf{x} ) $ is $ \mathbf{m}_n( r \midil \mathbf{x} ) = [ \int \mathbf{y} I( \mathbf{y} \in D_n( \alpha_n( 1 - r ) \midil \mathbf{x} ) ) \mu_n( d \mathbf{y} \midil \mathbf{x} ) ] \allowbreak / \mu_n( D_n( \alpha_n( 1 - r ) \midil \mathbf{x} ) \midil \mathbf{x} ) $.
Conditional trimmed means can be used to detect conditional skewness in the data. When the distribution is symmetric, the trimmed means and the median coincide, whereas for a skewed distribution, the trimmed means lie away from the median.

\subsection{Data demonstration} \label{subsec:2.3}
We demonstrate the conditional central regions, the conditional medians and the conditional trimmed means in the two real data examples considered in \autoref{sec:1}.
We choose the weights $ \{ W_{i,n}( \mathbf{x} ) \} $ in a way which ensures the asymptotic consistency of the estimates. For this, we have employed the nearest neighbor approach, where the integer $ k = k_n $ is
\begin{align}
k_n = \lfloor ( \log n )^2 \rfloor + 1 ,
\label{eq:2.1}
\end{align}
$ \lfloor r \rfloor $ being the largest integer less than or equal to the real number $ r $. So,
\begin{align*}
h( \mathbf{x}, k_n, n ) = \min\left\{ h \mid \sum_{i=1}^{n} I( d( \mathbf{x}, \mathbf{X}_i ) \le h ) > ( \log n )^2 \right\} .
\end{align*}
This approach is equivalent to the kernel method of choosing the weights when the kernel function is $ K( u ) = I( 0 \le u \le 1 ) $, and bandwidth $ h = h( \mathbf{x}, k_n, n ) $.
We shall see in \autoref{sec:5} that such a choice ensures the asymptotic consistency of the estimates.
In the demonstrations, we use the conditional halfspace depth.
With our choice of the weights, the conditional sample depth $ \rho_n( \cdot \midil \mathbf{x} ) $ becomes the corresponding sample depth based on the local response values of $ \mathbf{x} $, which we compute using algorithms by \cite{rousseeuw1996algorithm} and \cite{rousseeuw1998computing}.

\begin{figure}
	\centering
	\includegraphics[width=0.9\linewidth]{"paper3sliceplotrealbudgetitaly1a"}
	\caption{The local response values (1st row), the local boxplots for the 1st (2nd row) and the 2nd coordinates (3rd row) and the conditional 50\% central regions along with conditional medians (circle) and conditional 10\% trimmed means (cross) (4th row) in the BudgetItaly Data.}
	\label{fig:paper3sliceplotrealbudgetitaly}
\end{figure}
\begin{figure}
	\centering
	\includegraphics[width=0.9\linewidth]{"paper3sliceplotrealtecator1a"}
	\caption{The selected covariate curves (1st row), the local response values (2nd row), the local boxplots for the 1st (3rd row) and the 2nd coordinates (4th row) and the conditional 50\% central regions along with conditional medians (circle) and conditional 10\% trimmed means (cross) (5th row) in the Tecator Data.}
	\label{fig:paper3sliceplotrealtecator}
\end{figure}

We present the 50\% conditional central regions for the conditional halfspace depth corresponding to four selected covariate values in the BudgetItaly Data in \autoref{fig:paper3sliceplotrealbudgetitaly}, along with the scatter plots of the local response values and the corresponding local boxplots. The circles and the crosses inside the conditional central regions denote the conditional medians and the 10\% conditional trimmed means, respectively. Similar plots for the Tecator Data are presented in \autoref{fig:paper3sliceplotrealtecator}.

The scatter plots of the local response values in \autoref{fig:paper3sliceplotrealbudgetitaly} and \autoref{fig:paper3sliceplotrealtecator} indicate the correlation between the two coordinates of the respective response variables. The correlation is higher in the case of the Tecator Data. The local boxplots and the conditional central regions both demonstrate the heteroscedasticity present in the two datasets. The local boxplots and the conditional medians and the conditional trimmed means reflect the variation of conditional skewness of the response over the covariate values in both the datasets.

\section{Conditional Spread and Heteroscedasticity} \label{sec:3}
We define a measure of conditional spread based on the conditional central region, and develop a test for heteroscedasticity in this section.
In \cite{dette1998testing}, \cite{dette2002consistent} and \cite{dette2010robust}, some tests of heteroscedasticity were considered in a nonparametric regression setup with real valued response and covariate. Some tests of heteroscedasticity were also proposed in the case of parametric regression with multivariate response \citep{gupta1984distribution,holgersson2004testing}. Our test of heteroscedasticity is developed in a nonparametric setup, where the response is multivariate and the covariate lies in a complete separable metric space. This covers the cases of real valued response as well as finite and infinite dimensional covariates. To the best of our knowledge, there is no other nonparametric test of heteroscedasticity for a multivariate response as in our case.

A large conditional central region indicates a high spread of the conditional distribution of the response.
In \cite{liu1999multivariate}, a measure of spread based on central regions in an unconditional setup was considered. The analogous measure of conditional spread in our regression setup is
$ S( r \midil \mathbf{x} ) = \text{volume}( D( \alpha( r ) \midil \mathbf{x} ) ) = \lambda( D( \alpha( r ) \midil \mathbf{x} ) ) $,
where $ \lambda( \cdot ) $ is the Lebesgue measure in $ \mathbb{R}^p $.
We can estimate $ S( r \midil \mathbf{x} ) $ by its sample analogue $ S_n( r \midil \mathbf{x} ) = \text{volume}( D_n( \alpha_n( r ) \midil \mathbf{x} ) ) = \lambda( D_n( \alpha_n( r ) \midil \mathbf{x} ) ) $. In \autoref{sec:5}, we show that $ S_n( r \midil \mathbf{x} ) $ is a consistent estimate of $ S( r \midil \mathbf{x} ) $. In \cite{liu1999multivariate}, a different estimate of $ S( r \midil \mathbf{x} ) $ was considered, which is the volume of the convex hull of the sample observations lying in $ D_n( \alpha_n( r ) \midil \mathbf{x} ) $. But this estimate may not be consistent when the conditional central region $ D( \alpha( r ) \midil \mathbf{x} ) $ is not convex. The measure $ S_n( r \midil \mathbf{x} ) $ is computationally intensive, which causes problems when using it for our test of heteroscedasticity based on a permutation procedure. For this reason, we consider an alternative spread measure $ \Delta( r \midil \mathbf{x} ) $ in our data analysis.
The measure $ \Delta( r \midil \mathbf{x} ) $ of conditional spread of $ \mathbf{Y} $ given $ \mathbf{X} = \mathbf{x} $ is defined as the diameter of the set $ D( \alpha( r ) \midil \mathbf{x} ) $, i.e., $ \Delta( r \midil \mathbf{x} ) = \sup\{ \| \mathbf{y}_1 - \mathbf{y}_2 \| \midil \mathbf{y}_1, \mathbf{y}_2 \in D( \alpha( r ) \midil \mathbf{x} ) \} $.
For a real valued response, $ \Delta( 0.5 \midil \mathbf{x} ) $ coincides with the conditional interquartile range.
We estimate $ \Delta( r \midil \mathbf{x} ) $ by
$ \Delta_n( r \midil \mathbf{x} ) = \max\{ \| \mathbf{Y}_i - \mathbf{Y}_j \| \midil \rho_n( \mathbf{Y}_i \midil \mathbf{x} ) \ge \alpha_n( r ),\, \rho_n( \mathbf{Y}_j \midil \mathbf{x} ) \ge \alpha_n( r );\, i,j = 1, ..., n \} $.
In \autoref{sec:5}, we show that the estimate $ \Delta_n( r \midil \mathbf{x} ) $ is asymptotically consistent.
In \cite{chowdhury2016nonparametric}, a similar measure of conditional spread using the kernel estimate of the conditional spatial depth was considered. But they did not investigate its asymptotic consistency, which now follows as a special case of \autoref{thm:3} in \autoref{sec:5} for finite dimensional responses.

\subsection{Data Demonstration}
We now demonstrate the conditional spread measure in the BudgetItaly Data and the Tecator Data considered in \emph{Example 1} and \emph{Example 2}, respectively, in \autoref{sec:1}. For the BudgetItaly Data, we take the pair of food share and housing and fuels share as the response, and the variable pair of food price and housing and fuels price as the covariate.
For the Tecator Data, we consider two regression problems. In the first case, we take the pair of fat and protein contents as a bivariate response, and in the second case, we take the triplet of moisture, fat and protein contents as a trivariate response. The covariate is same in both cases, namely, the curve of absorbance spectra.
For both the datasets, we compute $ \Delta_n( r \midil \mathbf{x} ) $ with $ r = 0.5 $ and $ \mathbf{x} $ varying over the sample covariate values $ \mathbf{X}_i, i = 1, ..., n $. The choice of weights and the depth function are the same as in \autoref{subsec:2.3}. We compute the scores of the observations with respect to the first and the second principal components of the sample dispersion of the covariate, and denote the two scores as $ P_1 $ and $ P_2 $, respectively. For the BudgetItaly Data, we plot $ \Delta_n( r \midil \mathbf{x} ) $ against $ P_1 $ and $ P_2 $ in \autoref{fig:paper3spreadbudgetitaly}.
\begin{figure}
	\centering
	\includegraphics[width=1\linewidth]{"paper3spreadbudgetitaly2d1"}
	\caption{Plot of $ \Delta_n( r \midil \mathbf{x} ) $ with $ r = 0.5 $ against $ P_1 $ and $ P_2 $ in the BudgetItaly Data.}
	\label{fig:paper3spreadbudgetitaly}
\end{figure}
For the Tecator Data, we plot $ \Delta_n( r \midil \mathbf{x} ) $ for both the bivariate response and the trivariate response against $ P_1 $ and $ P_2 $ in \autoref{fig:paper3spreadtecator}.
\begin{figure}
	\centering
	\includegraphics[width=1\linewidth]{"paper3spreadtecator2d1"}
	\includegraphics[width=1\linewidth]{"paper3spreadtecator3d1"}
	\caption{Plots of $ \Delta_n( r \midil \mathbf{x} ) $ with $ r = 0.5 $ against $ P_1 $ and $ P_2 $ for bivariate and trivariate responses in the Tecator Data.}
	\label{fig:paper3spreadtecator}
\end{figure}
We notice that the indication of heteroscedasticity in the BudgetItaly Data is not so strong in \autoref{fig:paper3spreadbudgetitaly}, compared to that for the Tecator Data in \autoref{fig:paper3spreadtecator}. In the Tecator Data, we can see clear pattern in the plots of $ \Delta_n( r \midil \mathbf{x} ) $ against $ P_1 $ for both the bivariate and the trivariate responses in \autoref{fig:paper3spreadtecator}, which indicates the presence of heteroscedasticity. The plots of $ \Delta_n( r \midil \mathbf{x} ) $ against $ P_2 $ does not show much heteroscedasticity.

\subsection{A Test for Heteroscedasticity}
We now develop a nonparametric test for heteroscedasticity based on the conditional central region. In the presence of heteroscedasticity, $ \Delta( r \midil \mathbf{x} ) $ will vary with $ \mathbf{x} $. To capture the variation of $ \Delta( r \midil \mathbf{x} ) $ over $ \mathbf{x} $, we compute $ \Delta_n( r \midil \mathbf{X}_i ) $ for $ i = 1, ..., n $. Our hypotheses are
\begin{center}
\begin{minipage}[c]{.5\textwidth}
	\begin{itemize}
	\item[$ H_0 : $] $ \Delta( r \midil \mathbf{x} ) $ is constant over $ \mathbf{x} $,
	\item[$ H_A : $] $ \Delta( r \midil \mathbf{x} ) $ varies with $ \mathbf{x} $.
	\end{itemize}
\end{minipage}
\end{center}

We define our test-statistic as
\begin{align*}
	\mathbf{T}_n = \frac{1}{n} \sum_{i=1}^{n} \left[ \Delta_n( r \midil \mathbf{X}_i ) - \left( \frac{1}{n} \sum_{j=1}^{n} \Delta_n( r \midil \mathbf{X}_j ) \right) \right]^2 .
\end{align*}
Large values of $ \mathbf{T}_n $ bear evidence against $ H_0 $.
The p-value of the test is computed using a permutation procedure. Let $ \mathbf{S}_n $ be the set of all permutations of the integers $ 1, ..., n $, defined as $ \mathbf{S}_n = \{ \sigma \midil \sigma : \{ 1, ..., n \} \to \{ 1, ..., n \}, \, \sigma \text{ is 1-1 and onto} \} $. Consider the set $ \mathbb{S}_n = \{ ( \mathbf{X}_1, \mathbf{Y}_{\sigma( 1 )} ), ..., ( \mathbf{X}_n, \mathbf{Y}_{\sigma( n )} ) \allowbreak \midil \sigma \in \mathbf{S}_n \} $, the set of all permuted samples, where the response is freely permuted.
We compute the value of $ \mathbf{T}_n $ for all the permuted samples in $ \mathbb{S}_n $, and the empirical distribution of those values can be taken as an approximation of the null distribution of $ \mathbf{T}_n $. The p-value is computed as the proportion of those values of $ \mathbf{T}_n $ which are larger than the actually observed value of $ \mathbf{T}_n $. In practice, the number of all permutations is too large for even moderate sample sizes, and we consider a fixed number of random permutations of $ 1, ..., n $ to calculate the p-value based on that.
\begin{table}[h]
\caption{p-values for the test of heteroscedasticity based on several conditional depths in the BudgetItaly Data and the Tecator Data}
\begin{center}
\begin{tabular}{cccccc}
	\hline
	Data		& Response type		& Halfspace		& Spatial		& Projection	& Simplicial	\\\hline
	BudgetItaly	& Bivariate			& 0.018			& 0.006			& 0				& 0.046			\\\hline
	Tecator		& Bivariate			& 0				& 0				& 0				& 0				\\
	Tecator		& Trivariate		& 0.042			& 0				& 0.002			& 0				\\\hline
\end{tabular}
\end{center}
\label{table1}
\end{table}
We present the computed p-values for the test based on 500 random permutations in the BudgetItaly Data and the Tecator Data in \autoref{table1}. Overall, the p-values indicate strong presence of heteroscedasticity in both the data. We note that, though the plots of the conditional spread measure in the BudgetItaly Data in \autoref{fig:paper3spreadbudgetitaly} did not clearly indicate the presence of heteroscedasticity, our test detects the presence of heteroscedasticity in the data. The p-values for the Tecator Data are consistent with the plots in \autoref{fig:paper3spreadtecator} as both indicate presence of heteroscedasticity.

\subsection{Level and Power study for the test of heteroscedasticity} \label{supsec:1}
We present a level and power study for our test of heteroscedasticity. We consider four regression models. Let $ \Sigma_p = (( \sigma_{ij} ))_{p \times p} $ with $ \sigma_{ij} = 0.5 + 0.5 I( i = j ) $, and $ a $ be a positive number. In the first model, we consider a trivariate covariate and a bivariate response. The trivariate covariate $ \mathbf{X} = ( \mathbf{X}^{(1)}, \mathbf{X}^{(2)}, \mathbf{X}^{(3)} ) $ is such that $ \mathbf{X}^{(1)} $, $ \mathbf{X}^{(2)} $ and $ \mathbf{X}^{(3)} $ have independent $ Uniform[ 0, 1.5 ] $ distribution. The conditional distribution of the response $ \mathbf{Y} $ given covariate $ \mathbf{X} $ is bivariate normal with mean vector $ ( 0, 0 ) $ and dispersion matrix $ (1 + a ( \mathbf{X}^{(1)} \mathbf{X}^{(2)} \mathbf{X}^{(3)} ) ) \Sigma_2 $. In the second model, the covariate is the same as in the first model, but the response $ \mathbf{Y} $ is trivariate, with its conditional distribution given $ \mathbf{X} $ being trivariate normal with mean vector $ ( 0, 0, 0 ) $ and dispersion matrix $ (1 + a ( \mathbf{X}^{(1)} \mathbf{X}^{(2)} \mathbf{X}^{(3)} ) ) \Sigma_3 $.
\begin{table}[h]
\caption{Estimated powers for bivariate and trivariate responses with trivariate covariate}
\begin{center}
\begin{tabular}{cc|ccccc|ccccc}
	\hline
	\multicolumn{2}{c|}{}	& \multicolumn{5}{c|}{Bivariate $ \mathbf{Y} $, trivariate $ \mathbf{X} $}
							& \multicolumn{5}{c}{Trivariate $ \mathbf{Y} $, trivariate $ \mathbf{X} $}\\\hline
	$ n $					& level	& $ a = 0 $	& $ a = 2 $	& $ a = 4 $	& $ a = 6 $	& $ a = 8 $
									& $ a = 0 $	& $ a = 2 $	& $ a = 4 $	& $ a = 6 $	& $ a = 8 $\\\hline
	100						& 5\%	& 0.052		& 0.322		& 0.576		& 0.724		& 0.800
									& 0.054		& 0.098		& 0.184		& 0.254		& 0.348	\\
	100						& 1\%	& 0.008		& 0.142		& 0.318		& 0.486		& 0.574
									& 0.006		& 0.022		& 0.070		& 0.114		& 0.170	\\\hline
	200						& 5\%	& 0.046		& 0.718		& 0.948		& 0.986		& 0.996
									& 0.052		& 0.330		& 0.588		& 0.760		& 0.858	\\
	200						& 1\%	& 0.008		& 0.500		& 0.844		& 0.936		& 0.976
									& 0.010		& 0.146		& 0.390		& 0.536		& 0.660	\\\hline
	300						& 5\%	& 0.044		& 0.896		& 0.998		& 1.000		& 1.000
									& 0.070		& 0.626		& 0.936		& 0.992		& 0.998	\\
	300						& 1\%	& 0.014		& 0.734		& 0.976		& 1.000		& 1.000
									& 0.012		& 0.374		& 0.778		& 0.924		& 0.980	\\\hline
\end{tabular}
\end{center}
\label{table:1}
\end{table}
In the third and the fourth regression models, we consider a functional covariate $ \mathbf{X} $, such that $ \mathbf{X}( t ) = \mathbf{B} e^t $, with $ \mathbf{B} \sim Uniform[ 0, 1 ] $ and $ t \in [ 0, 1 ] $. $ \mathbf{X} $ is considered as a random element in $ L_2[ 0, 1 ] $. In the third model, we take the conditional distribution of the response $ \mathbf{Y} $ given $ \mathbf{X} $ to be bivariate normal with mean vector $ ( 0, 0 ) $ and dispersion matrix $ ( 1 + a \| \mathbf{X} \|_2 ) \Sigma_2 $, where $ \| \cdot \|_2 $ is the $ L_2 $-norm. In the fourth model, the conditional distribution of the response $ \mathbf{Y} $ given $ \mathbf{X} $ is trivariate normal with mean vector $ ( 0, 0, 0 ) $ and dispersion matrix $ ( 1 + a \| \mathbf{X} \|_2 ) \Sigma_3 $.
We estimate the power of the test of heteroscedasticity based on conditional halfspace depth for different sample sizes and different values of the constant $ a $ based on 500 independent replications, with the number of random permutations in each replication being also 500. The estimated powers for the first and the second models are presented in \autoref{table:1}, while those for the third and the fourth models are presented in \autoref{table:2}.
\begin{table}
\caption{Estimated powers for bivariate and trivariate responses with functional covariate}
\begin{center}
\begin{tabular}{cc|ccccc|ccccc}
	\hline
	\multicolumn{2}{c|}{}	& \multicolumn{5}{c|}{Bivariate $ \mathbf{Y} $, functional $ \mathbf{X} $}
							& \multicolumn{5}{c}{Trivariate $ \mathbf{Y} $, functional $ \mathbf{X} $}\\\hline
	$ n $					& level	& $ a = 0 $	& $ a = 2 $	& $ a = 4 $	& $ a = 6 $	& $ a = 8 $
									& $ a = 0 $	& $ a = 2 $	& $ a = 4 $	& $ a = 6 $	& $ a = 8 $\\\hline
	100						& 5\%	& 0.056		& 0.706		& 0.740		& 0.752		& 0.760
									& 0.050		& 0.316		& 0.376		& 0.400		& 0.400\\
	100						& 1\%	& 0.006		& 0.436		& 0.458		& 0.476		& 0.488
									& 0.004		& 0.110		& 0.150		& 0.172		& 0.172\\\hline
	200						& 5\%	& 0.038		& 0.982		& 0.996		& 0.994		& 0.994
									& 0.052		& 0.750		& 0.786		& 0.796		& 0.800\\
	200						& 1\%	& 0.008		& 0.896		& 0.916		& 0.928		& 0.928
									& 0.004		& 0.512		& 0.560		& 0.588		& 0.586\\\hline
	300						& 5\%	& 0.052		& 0.998		& 0.998		& 1.000		& 1.000
									& 0.048		& 0.934		& 0.958		& 0.966		& 0.972\\
	300						& 1\%	& 0.008		& 0.984		& 0.992		& 0.996		& 0.994
									& 0.008		& 0.722		& 0.772		& 0.786		& 0.786\\\hline
\end{tabular}
\end{center}
\label{table:2}
\end{table}
Note that the estimated power corresponding to $ a = 0 $ is nothing but the estimated size of the test. The observed power increases with the increment in the value of $ a $, which is a parameter that controls the extent of heteroscedasticity in the model. Evidently, the power of the test tends to 1 as the sample size increases.

\section{Asymptotic Results} \label{sec:5}
Recall the definitions of $ \mu( \cdot \midil \mathbf{z} ) $, $ \mu_n( \cdot \midil \mathbf{x} ) $, $ \rho( \cdot \midil \mathbf{x} ) $ and $ \rho_n( \cdot \midil \mathbf{x} ) $ from \autoref{subsec:2.1}.
We shall first establish the uniform asymptotic consistency for some standard conditional depth functions.
We need the following conditions on the weights $ W_{i,n}( \mathbf{x} ) $ and the conditional probability $ \mu( \cdot \midil \mathbf{z} ) $ of $ \mathbf{Y} $ given $ \mathbf{X} = \mathbf{z} $.
\begin{enumerate}[label=C(i.\alph*), ref=C(i.\alph*)]
\item \label{a1}
$ ( \log n ) \sum_{i=1}^{n} W_{i,n}^2( \mathbf{x} ) \stackrel{a.s.}{\longrightarrow} 0 $ and $ ( \log n ) \max_{1 \le i \le n} W_{i,n}( \mathbf{x} ) \stackrel{a.s.}{\longrightarrow} 0 $ as $ n \to \infty $,
and for any $ \delta > 0 $, $ \sum_{i=1}^{n} W_{i,n}( \mathbf{x} ) I( d( \mathbf{x}, \mathbf{X}_i ) \ge \delta ) \allowbreak \stackrel{a.s.}{\longrightarrow} 0 $ as $ n \to \infty $.

\item \label{a1b}
There is a collection of indices $ \mathbf{S}_n \subset \{ 1, 2, ..., n \} $ with cardinality $ k_n $ such that for any $ \delta > 0 $, $ \mathbf{S}_n \subset \{ i \midil d( \mathbf{x}, \mathbf{X}_i ) < \delta , \, i = 1, ..., n \} $ for all sufficiently large $ n $. Also,
\begin{align*}
& \frac{\log n}{k_n} \to 0
\text{ as } n \to \infty , \\
& \frac{\min_{i \in \mathbf{S}_n} W_{i,n}( \mathbf{x} )}{\max_{i \in \mathbf{S}_n} W_{i,n}( \mathbf{x} )}
\text{ is \emph{almost surely} bounded away from 0 as } n \to \infty , \\
& \text{and } \frac{\sum_{\{ i_1, ..., i_{p+1} \} \not\subset \mathbf{S}_n} W_{i_1,n}( \mathbf{x} ) ... W_{i_{p+1},n}( \mathbf{x} )}{\sum_{\{ j_1, ..., j_{p+1} \} \subset \mathbf{S}_n} W_{j_1,n}( \mathbf{x} ) ... W_{j_{p+1},n}( \mathbf{x} )}
\stackrel{a.s.}{\longrightarrow} 0
\text{ as }
n \to \infty .
\end{align*}
\end{enumerate}
\begin{enumerate}[label=C(\roman*), ref=C(\roman*)]
\setcounter{enumi}{1}
\item \label{assume:2}
$ \mu( \cdot \midil \mathbf{z} ) \stackrel{w}{\longrightarrow} \mu( \cdot \midil \mathbf{x} ) $ as $ \mathbf{z} \to \mathbf{x} $.
\end{enumerate}
Suppose that $ P[ d( \mathbf{x}, \mathbf{X} ) \le h ] > 0 $ for all $ h > 0 $.
Then, in the nearest neighbor method of constructing weights $ W_{i,n}( \mathbf{x} ) $, conditions \ref{a1} and \ref{a1b} hold when $ k_n = \lfloor ( \log n )^2 \rfloor + 1 $.
In the kernel method, suppose that the kernel function $ K( \cdot ) $ satisfies $ l I( 0 \le s \le 1 ) \le K( s ) \le u I( 0 \le s \le 1 ) $ for some constants $ 0 < l, u < \infty $.
Depending on the distribution of the covariate $ \mathbf{X} $, we can choose a sequence of bandwidths $ \{ h_n \} $ such that $ h_n \to 0 $ and $ [ n P[ d( \mathbf{x}, \mathbf{X} ) \le h_n ] ]^{-1} ( \log n ) \to 0 $ as $ n \to \infty $.
For this choice of $ \{ h_n \} $ and the kernel function, it can be shown using the Bernstein inequality (see \citet[p.~95, Lemma A]{serfling2009approximation}) and the Borel-Cantelli Lemma that conditions \ref{a1} and \ref{a1b} are satisfied.

Assumption \ref{assume:2} holds in many common situations that we encounter. As an example, consider the location-scale model defined by $ \mathbf{Y} = l( \mathbf{X} ) + s( \mathbf{X} ) \mathbf{G} $. Here, $ \mathbf{G} $ is a random vector independent of $ \mathbf{X} $, and the functions $ l( \cdot ) : \mathcal{C} \to \mathbb{R}^p $ and $ s( \cdot ) : \mathcal{C} \to \mathbb{R} $ are both continuous at $ \mathbf{x} $. In such a setup, $ \mu( \cdot \midil \mathbf{z} ) \stackrel{w}{\longrightarrow} \mu( \cdot \midil \mathbf{x} ) $ as $ \mathbf{z} \to \mathbf{x} $.
Under \ref{a1} and \ref{assume:2}, it follows that
\begin{align}
\mu_n( \cdot \midil \mathbf{x} ) \stackrel{w}{\longrightarrow} \mu( \cdot \midil \mathbf{x} )
\quad
\text{\emph{almost surely} as }
n \to \infty .
\label{assume:1}
\end{align}
See Lemma \ref{lemma1a} for the proof of \eqref{assume:1}.
The following theorem states the uniform consistency for several depth functions.
\begin{theorem} \label{thm:1a}
Let $ \rho( \mathbf{y} \midil \mathbf{x} ) $ be any of the four conditional depth functions described in \autoref{subsec:2.1}. Suppose that there is a neighborhood of $ \mathbf{x} $ such that for all $ \mathbf{z} $ in that neighborhood, the conditional distribution of $ \mathbf{Y} $ given $ \mathbf{X} = \mathbf{z} $ has a continuous positive density $ f( \cdot \midil \mathbf{z} ) $ that is continuous in $ \mathbf{z} $.
Then, under conditions \ref{a1}, \ref{a1b} and \ref{assume:2},
\begin{align}
\sup_{ \mathbf{y} \in \mathbb{R}^p } | \rho_n( \mathbf{y} \midil \mathbf{x} ) - \rho( \mathbf{y} \midil \mathbf{x} ) | \stackrel{P}{\longrightarrow} 0
\quad
\text{as }
n \to \infty .
\label{a3}
\end{align}
\end{theorem}
Examples of continuous positive conditional densities $ f( \cdot \midil \mathbf{z} ) $ that are continuous in $ \mathbf{z} $ include the location-scale model: $ \mathbf{Y} = l( \mathbf{X} ) + s( \mathbf{X} ) \mathbf{G} $, where $ \mathbf{G} $ is a random vector independent of $ \mathbf{X} $ and has a continuous positive density on $ \mathbb{R}^p $, and the functions $ l( \cdot ) : \mathcal{C} \to \mathbb{R}^p $ and $ s( \cdot ) : \mathcal{C} \to \mathbb{R} $ are both continuous in a neighborhood of $ \mathbf{x} $.

\subsection{Asymptotic consistency of conditional central regions and related measures}
\label{subsec:5.1}
We now proceed to state the asymptotic consistency of the conditional central regions, medians and trimmed means, and spread measures. Recall the definitions of $ D( \cdot \midil \mathbf{x} ) $, $ D_n( \cdot \midil \mathbf{x} ) $, $ \alpha( r ) $, $ \alpha_n( r ) $, $ \mathbf{m}( \mathbf{x} ) $, $ \mathbf{m}_n( \mathbf{x} ) $, $ \mathbf{m}( r \midil \mathbf{x} ) $ and $ \mathbf{m}_n( r \midil \mathbf{x} ) $ from \autoref{subsec:2.1}. We need the following conditions.
\begin{enumerate}[label=C(\roman*), ref=C(\roman*)]
\setcounter{enumi}{2}
\item \label{a4}
	$ \mu( \{ \mathbf{y} \midil \rho( \mathbf{y} \midil \mathbf{x} ) = \alpha \} \midil \mathbf{x} ) = 0 $ for all $ \alpha $.
	
\item \label{a5}
	For any $ 0 \le \gamma_1 < \gamma_2 \le \rho( \mathbf{m}( \mathbf{x} ) \midil \mathbf{x} ) $, $ \mu( \{ \mathbf{y} \midil \gamma_1 \le \rho( \mathbf{y} \midil \mathbf{x} ) \le \gamma_2 \} \midil \mathbf{x} ) > 0 $.
	
\item \label{a2}
	$ \rho( \mathbf{y} \midil \mathbf{x} ) $ is a continuous function of $ \mathbf{y} $, and $ \rho( \mathbf{y} \midil \mathbf{x} ) \to 0 $ as $ \| \mathbf{y} \| \to \infty $.
	
\item \label{a6}
	Define $ D_0( \alpha \midil \mathbf{x} ) = \{ \mathbf{y} \midil \rho( \mathbf{y} \midil \mathbf{x} ) > \alpha \} $. Then, the closure of $ D_0( \alpha \midil \mathbf{x} ) $ is $ D( \alpha \midil \mathbf{x} ) $ for any $ 0 < \alpha < \rho( \mathbf{m}( \mathbf{x} ) \midil \mathbf{x} ) $.
\end{enumerate}
Assumptions \ref{a4} and \ref{a5} imply that the distribution function of the random variable $ \rho( \mathbf{Y} \midil \mathbf{x} ) $ is continuous and strictly increasing. Assumptions \ref{a2} and \ref{a6} are related to the smoothness of the conditional depth and the corresponding central regions.

Suppose $ \mu( \cdot \midil \mathbf{x} ) $ has a probability density with a convex support. Then, from the proof of Lemma 6.1 in \cite{donoho1992breakdown}, we get that \ref{a2} is satisfied for the conditional halfspace depth. From the arguments in the proof of Lemma 6.3 in \cite{donoho1992breakdown}, we get that $ D( \alpha \midil \mathbf{x} ) $ is convex for all $ 0 < \alpha < \rho( \mathbf{m}( \mathbf{x} ) \midil \mathbf{x} ) $, and consequently, \ref{a4}, \ref{a5} and \ref{a6} are also satisfied for the conditional halfspace depth.

Now suppose that $ \mu( \cdot \midil \mathbf{x} ) $ has a positive probability density on $ \mathbb{R}^p $. Then, the assumptions \ref{a4}, \ref{a5} and \ref{a2} hold for the conditional spatial depth. Note that when assumption \ref{a2} holds, assumption \ref{a6} is also satisfied if the conditional depth has no local maximum. From the arguments in section 2 in \cite{chowdhury2016nonparametric}, we get that the conditional spatial median is unique. Using arguments similar to those in the proofs of Lemma 2.5 and Lemma 2.6 in the supplement of \cite{chowdhury2016nonparametric}, one can show that the conditional spatial depth $ \rho( \mathbf{y} \midil \mathbf{x} ) $ has non-zero Fr\'{e}chet derivative with respect to $ \mathbf{y} $ for $ \mathbf{y} $ not being the conditional spatial median. Consequently, the conditional spatial depth cannot have any local maximum and hence it satisfies assumption \ref{a6}.

From Theorem 3.4 in \cite{zuo2000structural}, it follows that it is sufficient for $ \mu( \cdot \midil \mathbf{x} ) $ to have an elliptically symmetric probability density on $ \mathbb{R}^p $ for the associated conditional projection depth and conditional simplicial depth to satisfy assumptions \ref{a4}, \ref{a5}, \ref{a2} and \ref{a6}.

In the theorem below, we establish an asymptotic consistency result for the conditional sample central regions.
Note that the Hausdorff distance $ d_H( A, B ) $ between two nonempty closed and bounded subsets of a metric space is defined as $ d_H( A, B ) = \inf\{ \epsilon \midil A \subseteq B^\epsilon , B \subseteq A^\epsilon \} $, where $ A^\epsilon $ and $ B^\epsilon $ denote the $ \epsilon $-neighborhoods of $ A $ and $ B $, respectively (see, e.g., \citet[p.~281]{munkres2000topology}).
\begin{theorem} \label{thm:1}
Suppose \ref{a4}, \ref{a5}, \eqref{assume:1} and \eqref{a3} are satisfied. Then, for any $ \epsilon > 0 $ and any $ 0 < r < 1 $,
\begin{align*}
P\left[ D( \alpha( r ) + \epsilon \midil \mathbf{x} ) \subseteq D_n( \alpha_n( r ) \midil \mathbf{x} ) \subseteq D( \alpha( r ) - \epsilon \midil \mathbf{x} ) \right] \to 1
\quad \text{as } n \to \infty .
\end{align*}
If in addition, conditions \ref{a2} and \ref{a6} hold, then
$ d_H( D_n( \alpha_n( r ) \midil \mathbf{x} ) , D( \alpha( r ) \midil \mathbf{x} ) ) \stackrel{P}{\longrightarrow} 0 $
as $ n \to \infty $ for any $ 0 < r < 1 $.
\end{theorem}

The conditional $ \alpha $-depth contour $ \delta( \alpha \midil \mathbf{x} ) $ of $ \mathbf{Y} $ given $ \mathbf{X} = \mathbf{x} $ is defined as the boundary of $ D( \alpha \midil \mathbf{x} ) $, and for $ 0 < r < 1 $, the conditional $ 100 r $\% central region contour of $ \mathbf{Y} $ given $ \mathbf{X} = \mathbf{x} $ is defined as $ \delta( \alpha( r ) \midil \mathbf{x} ) $. The sample counterparts of $ \delta( \alpha \midil \mathbf{x} ) $ and $ \delta( \alpha( r ) \midil \mathbf{x} ) $ are denoted as $ \delta_n( \alpha \midil \mathbf{x} ) $ and $ \delta_n( \alpha_n( r ) \midil \mathbf{x} ) $, respectively. A consequence of \autoref{thm:1}, stated below in Corollary \ref{coro1}, is that for any band around $ \delta( \alpha( r ) \midil \mathbf{x} ) $, which may be of arbitrarily small width, $ \delta_n( \alpha_n( r ) \midil \mathbf{x} ) $ eventually lies inside that band with high probability for large $ n $. Since the contours determine the shapes of the central regions, this implies that the shapes of the sample central regions are good approximations of their population counter-parts in large samples.
\begin{corollary} \label{coro1}
Suppose \ref{a4}, \ref{a5}, \eqref{assume:1} and \eqref{a3} are satisfied. Then, for any $ \epsilon > 0 $,
\begin{align*}
P\left[ \delta_n( \alpha_n( r ) \midil \mathbf{x} ) \subseteq \{ \mathbf{y} \in \mathbb{R}^p \midil \alpha( r ) - \epsilon \le \rho( \mathbf{y} \midil \mathbf{x} ) < \alpha( r ) + \epsilon \} \right] \to 1
\quad \text{as } n \to \infty .
\end{align*}
\end{corollary}

Convergence of sample central regions were studied earlier in the unconditional setup by \cite{he1997convergence} and \cite{zuo2000structural}. In \cite{he1997convergence}, the authors assumed that the population depth contours are elliptic in nature in order to establish the convergence of the contours of the central regions. In \cite{zuo2000structural}, the authors restricted themselves to elliptic distributions for proving the convergence of the central regions in the Hausdorff distance.

The next theorem states the asymptotic consistency of the sample conditional median and the sample conditional trimmed means.
\begin{theorem} \label{thm:2}
Let $ \mathbf{M}( \mathbf{x} ) $ be the set of all conditional medians corresponding to the population conditional depth function $ \rho( \mathbf{y} \midil \mathbf{x} ) $.
Suppose that conditions \ref{a4}--\ref{a6} are satisfied and \eqref{assume:1} and \eqref{a3} hold. Then,
for any sequence of sample conditional medians $ \{ \mathbf{m}_n( \mathbf{x} ) \} $,
$ \inf_{ \mathbf{m} \in \mathbf{M}( \mathbf{x} ) } \| \mathbf{m}_n( \mathbf{x} ) - \mathbf{m} \|
\stackrel{P}{\longrightarrow} 0 $
as $ n \to \infty $.
Also, $ \mathbf{m}_n( r \midil \mathbf{x} ) \stackrel{P}{\longrightarrow} \mathbf{m}( r \midil \mathbf{x} ) $ as $ n \to \infty $ for any $ 0 < r < 1 $.
\end{theorem}
We state the asymptotic consistency of the conditional spread estimates $ S_n( r \midil \mathbf{x} ) $ and $ \Delta_n( r \midil \mathbf{x} ) $ in the theorem below.
\begin{theorem} \label{thm:3}
(a) Suppose \ref{a4}, \ref{a5}, \eqref{assume:1} and \eqref{a3} hold, and
$ \mu( \cdot \midil \mathbf{x} ) $ has a continuous positive density on $ \mathbb{R}^p $.
Then,
$ S_n( r \midil \mathbf{x} ) \stackrel{P}{\longrightarrow} S( r \midil \mathbf{x} ) $
as $ n \to \infty $.

(b) Let $ 0 < r < 1 $, and suppose that for any open set $ G $ with $ G \cap D( \alpha( r ) \midil \mathbf{x} ) \neq \emptyset $, we have $ P\left[ \mathbf{Y} \in G \cap D( \alpha( r ) \midil \mathbf{x} ) \right] > 0 $.
Then, under \ref{a4}--\ref{a6}, \eqref{assume:1} and \eqref{a3}, $ \Delta_n( r \midil \mathbf{x} ) \stackrel{P}{\longrightarrow} \Delta( r \midil \mathbf{x} ) $ as $ n \to \infty $.
\end{theorem}

The first condition in part (b) in the above theorem is trivially satisfied for any depth function when the support of the distribution of $ \mathbf{Y} $ is the whole response space.
On the other hand, for particular depths like the conditional halfspace depth, the conditional central regions are bounded convex sets. So, if the support of the distribution of $ \mathbf{Y} $ is a bounded convex set, then also that condition is satisfied for the conditional halfspace depth.

