Robust subgroup-classifier learning and testing in change-plane regressions

Li et al., (2021) considered the change-plane model (1) with Gaussian errors. Also, Zhang et al., (2022) investigated a quantile regression with a change plane and derived the asymptotic normalities for the grouping difference parameter and the grouping parameter. However, although quantile or median regression models require no Gaussian or sub-Gaussian assumption, they essentially estimate the conditional quantile or median regression instead of the conditional mean regression. If the mean regression is of interest in practice, then these procedures are not feasible unless the error distribution is symmetric around zero, which may be too strong to cause the misspecification problem. See Fan et al., 2017b for some examples that demonstrate the distinction between conditional mean regression and conditional quantile or median regression.

Linear regression models with heavy-tailed errors are prevalent in the literature. Fan et al., 2017b proposed a robust estimator of high-dimensional mean regression in the absence of asymmetry and with light tail assumptions. Zhou et al., (2018) provided a robust M-estimation procedure with applications to dependence-adjusted multiple testing. Sun et al., (2020) and Wang et al., (2021) studied adaptive Huber regression for linear regression models with heavy-tailed errors. Chen and Zhou, (2020) investigated robust inference via multiplier bootstrap in multiple response regression models, constructing robust bootstrap confidence sets and addressing large-scale simultaneous hypothesis testing problems.

Mukherjee et al., (2022) studied the change-plane problem under heavy-tailed errors when $\boldsymbol{\alpha}=0$ and ${\boldsymbol{X}}=1$, which is a special case of the change-plane model (1), and they left the general change-plane model with heavy-tailed errors for future work. Herein, a new robust procedure is introduced to estimate parameters and consequently to learn the subgroup classifier. Figure 1 shows boxplots of the estimation errors of the parameter $\boldsymbol{\theta}=(\boldsymbol{\alpha}^{\rm T},\boldsymbol{\beta}^{\rm T})^{\rm T}$ with $L_{2}$ norm and the accuracies of the estimated subgroup classifier, where the $L_{2}$ norm of the estimation errors is defined as $\|\hat{\boldsymbol{\theta}}_{\tau,h}-\boldsymbol{\theta}^{*}\|$ with the robust estimator $\hat{\boldsymbol{\theta}}_{\tau,h}$ of the true grouping parameter $\boldsymbol{\theta}^{*}$, and the accuracy is defined as $\mbox{ACC}=1-n^{-1}\sum_{i=1}^{n}\left|{\boldsymbol{1}}({\boldsymbol{U}}_{i}^{\rm T}\hat{\boldsymbol{\gamma}}_{\tau,h}\geq 0)-{\boldsymbol{1}}({\boldsymbol{U}}_{i}^{\rm T}\boldsymbol{\gamma}^{*}\geq 0)\right|$ with the robust estimator $\hat{\boldsymbol{\gamma}}_{\tau,h}$ of the true parameter $\boldsymbol{\gamma}^{*}$. Here, the settings are $(p,q,r)=(3,3,3)$ and $n=(200,400,600)$ with 1000 repetitions, the heavy-tailed errors are generated from the Pareto distribution $Par(2,1)$ with shape parameter 2 and scale parameter 1, and $X_{1}=Z_{1}=1$ and $(X_{2},\cdots,X_{p})^{\rm T}=(Z_{2},\cdots,Z_{p})^{\rm T}$ and $(U_{2},\cdots,U_{r})^{\rm T}$ are generated independently from multivariate normal distributions $N({\boldsymbol{0}}_{p-1},\sqrt{3}{\boldsymbol{I}}_{p-1})$ and $N({\boldsymbol{0}}_{r-1},\sqrt{3}{\boldsymbol{I}}_{r-1})$, respectively; see Section 4 for details. As used by Zhang et al., (2022), the smooth function $K(u)=\{1+\exp(-u)\}^{-1}$ with smoothness parameter $h=\sqrt{\log(n)/n}$ is chosen, and the proposed robust estimation procedure (AHu) is compared with the method based on OLS (Li et al., 2021). Figure 1 sends the important message that in the presence of heavy tails, compared with the existing method (Li et al., 2021), the proposed robust estimators not only reduce the estimation error dramatically but also improve significantly the accuracy of the subgroup classifier.

Another goal of this paper is to test for the existence of subgroups, i.e.,

$\displaystyle H_{0}:\boldsymbol{\beta}={\boldsymbol{0}}\quad versus\quad H_{1}:\boldsymbol{\beta}\neq{\boldsymbol{0}}.$

(2)

Note that the grouping parameter $\boldsymbol{\gamma}$ is not identifiable under the null hypothesis.

The classic Wald-type test or score-based test is powerful in standard test problems when there is no identifiability problem in both the null and alternative hypotheses, but these common procedures are not feasible when nuisance parameters are present. Andrews and Ploberger, (1994) and Andrews and Ploberger, (1995) studied the weighted average exponential form, which was originally introduced by Wald, (1943). Davies, (1977) investigated well the supremum of the squared score test (SST) statistic for mixture models, which was applied by Song et al., (2009) and Kang et al., (2017) to semiparametric models in censoring data.

All the aforementioned testing methods are optimal tests based on the weighted average power criterion. However, because these optimal tests take the weighted exponential average of the classical tests over the grouping parametric space $\Theta_{\gamma}$, they may not only not perform well in practice when the dimension of $\Theta_{\gamma}$ is large but also give rise to a heavy-burden calculation of the p-value or the critical value. Instead, Liu et al., (2024) introduced a new test statistic by taking the weighted average of the SST (WAST) over $\Theta_{\gamma}$ and removing both the inverse of the covariance and the cross-interaction terms to overcome the drawbacks of SST. Thanks to its closed form, WAST achieves more-accurate type-I errors and significantly improved power and hence dramatically reduced computational time as a byproduct.

All the aforementioned test procedures require the important assumption of Gaussian or sub-Gaussian errors in the change-plane models, none of which apply to heavy-tailed data sets. Therefore, proposed herein is a robust test procedure based on WAST (Liu et al., 2024), called robust WAST (RWAST). Figure 2 shows the power curves of the proposed RWAST, WAST as introduced by Liu et al., (2024), and SST as considered by Davies, (1977); Kang et al., (2017). A total of 1000 bootstrap samples is set, and the other settings are the same as those in Section 1.1. With the nominal significance level $\alpha=0.05$, the type-I errors for these three methods are $(\mbox{rwast},\mbox{olsw},\mbox{olss})=(0.042,0.065,0.031)$ for $n=200$, $(\mbox{rwast},\mbox{olsw},\mbox{olss})=(0.043,0.047,0.118)$ for $n=400$, and $(\mbox{rwast},\mbox{olsw},\mbox{olss})=(0.033,0.047,0.135)$ for $n=600$. It follows that RWAST controls the type-I errors well, while those of SST based on quadratic loss are larger and deviate far from 0.05. Figure 2 shows that in the presence of heavy tails, the proposed RWAST achieves larger power in comparison with WAST based on the ordinary quadratic loss.

In summary, from the demonstration examples in Section 1.1 and Section 1.2, compared with the existing nonrobust methods for heavy-tailed data, the proposed robust estimation procedure is characterized by lower estimation errors and higher accuracy, and the robust test procedure has more-accurate type-I errors and larger power.

The main contributions of this paper are summarized as follows. First, the robust estimation procedure for the change-plane model (1) with heavy-tailed errors is investigated carefully. The proposed robust estimator adapts to the sample size, the robustification parameter in the Huber loss, the smoothness parameter when approximating the indicator function in the subgroup classifier, and the moments of errors. The sacrifices made in pursuit of robustness and smoothness are analyzed theoretically, with the bias involving the robustification parameter arising from the pursuit of robustness, and the one involving the smoothness parameter arising from the approximation to the indicator function. The nonasymptotic properties for parameters $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ are established, as well as those for the grouping parameter. The theoretical results reveal that the proposed estimators of the grouping difference parameter and baseline parameter have Gaussian-type deviations (Devroye et al., 2016). Also provided is the nonasymptotic Bahadur representation of the proposed robust estimators, which is convenient for deriving the classical asymptotic results needed for statistical inference such as hypothesis tests and constructing confidence regions. Extensive simulation studies show that the proposed robust estimation procedure is superior to its several competitors.

Second, for the change-plane model with heavy-tailed errors, RWAST is proposed, which makes full use of the data’s heavy-tailed information and overcomes the drawbacks of loss of power in practice and the heavy computational burden of SST when the dimension of the grouping parameter is large. The asymptotic distributions of the proposed RWAST under the null and alternative hypotheses are established based on the theory of degenerate U-statistics. As with exponential average tests, the proposed asymptotic distributions are not standard (e.g., the normal or $\chi^{2}$ distribution), so a novel bootstrap method that is easily implemented and theoretically guaranteed is introduced to mimic the critical value or p-value. Comprehensive simulation studies conducted with finite sample sizes and for various heavy-tailed error distributions show the excellent performance of the proposed RWAST, which improves the power significantly and reduces the computational burden dramatically.

In summary, a novel robust estimator is proposed that adapts to the sample size, dimension, robustification parameter, moments, and smoothness parameter in pursuit of the optimal tradeoff among bias, robustness, and smoothness. To the best of the author’s knowledge about change-plane analysis, the literature contains no nonasymptotic results with sub-Gaussian tails for parameters $\boldsymbol{\theta}$ and $\boldsymbol{\eta}$, and no nonasymptotic results with sub-exponential tails for the Bahadur representation of these parameters. Furthermore, a robust test procedure is proposed that improves on WAST in the change-plane model with heavy-tailed errors.

Here, some useful notation is introduced for convenience of expression. For a vector ${\boldsymbol{v}}\in\mathbb{R}^{d}$ and a square matrix $A=(a_{ij})\in\mathbb{R}^{d\times d}$, denote by $\|{\boldsymbol{v}}\|$ the Euclidean norm of ${\boldsymbol{v}}$, by $\mbox{trace}(A)=\sum_{i=1}^{d}a_{ii}$ the trace of $A$; $\|{\boldsymbol{v}}\|_{A}^{2}=\sum_{i,j}a_{ij}v_{i}v_{j}$ and ${\boldsymbol{v}}^{\otimes 2}={\boldsymbol{v}}{\boldsymbol{v}}^{\rm T}$. Denote by $\|A\|_{p}=\sup\{\|A{\boldsymbol{x}}\|_{p}:{\boldsymbol{x}}\in\mathbb{R}^{d},\|{\boldsymbol{x}}\|_{p}=1\}$ the induced operator norm for a matrix $A=(a_{ij})\in\mathbb{R}^{m\times d}$.

Denote by $\boldsymbol{\mathrm{P}}$ the ordinary probability measure such that $\boldsymbol{\mathrm{P}}f=\int fd\boldsymbol{\mathrm{P}}$ for any measurable function $f$, by $\mathbb{P}_{n}$ the empirical measure of a sample of random elements from $\boldsymbol{\mathrm{P}}$ such that $\mathbb{P}_{n}f=n^{-1}\sum_{i=1}^{n}f({\boldsymbol{V}}_{i})$, and by ${\mathbb{G}}_{n}$ the empirical process indexed by a class ${\cal F}$ of measurable functions such that ${\mathbb{G}}_{n}f=\sqrt{n}(\mathbb{P}_{n}-\boldsymbol{\mathrm{P}})f$ for any $f\in{\cal F}$. Let $L^{p}(Q)$ be the space of all measurable functions $f$ such that $\|f\|_{Q,p}:=(Q|f|^{p})^{1/p}<\infty$, where $p\in[1,\infty)$ and $(Q|f|^{p})^{1/p}$ denotes the essential supremum when $p=\infty$. Let $N(\epsilon,{\cal F},\|\cdot\|_{Q,2})$ be an $\epsilon$-covering number of ${\cal F}$ with respect to the $L^{2}(Q)$ seminorm $\|\cdot\|_{Q,2}$, where ${\cal F}$ is a class of measure functions and $Q$ is finite discrete.

The remainder of this paper is organized as follows. Section 2 provides the robust estimators for the grouping difference parameter and grouping parameter as well as the subgroup classifier, and theorems reveal that these estimators achieve Gaussian-type deviations. Also derived is the Bahadur representation of the robust estimators, and it is shown that the remainder of the Bahadur representation achieves sub-Gaussian tails. Section 3 presents the RWAST statistic and establishes its limiting distributions under the null and alternative hypotheses. Section 4 reports the results of simulation studies conducted to evaluate the finite-sample performance of the proposed methods with competitors in the change-plane models with heavy-tailed errors. The performance of the proposed methods is illustrated further by applying them to a medical dataset in Section 5. Finally, Section 6 concludes with remarks and further extensions. The proofs are provided in the Supplementary Material, and an R package named “wasthub” is available at https://github.com/xliusufe/wasthub.

Rewrite model (1) as

$\displaystyle y_{i}={\boldsymbol{X}}_{i}^{\rm T}\boldsymbol{\alpha}+{\boldsymbol{Z}}_{i}^{\rm T}\boldsymbol{\beta}{\boldsymbol{1}}(U_{1i}+{\boldsymbol{U}}_{2i}^{\rm T}\boldsymbol{\eta}\geq 0)+\epsilon_{i},$

(6)

where ${\boldsymbol{U}}_{i}=(U_{1i},{\boldsymbol{U}}_{2i}^{\rm T})^{\rm T}$, $\boldsymbol{\eta}=\gamma_{1}^{-1}\boldsymbol{\gamma}_{-1}$ with $\boldsymbol{\gamma}_{-1}=(\gamma_{2},\cdots,\gamma_{r})^{\rm T}$. To avoid the identifiability problem for $\boldsymbol{\eta}$, $\boldsymbol{\beta}\neq{\boldsymbol{0}}$ is assumed in this section. Denote $\boldsymbol{\zeta}=(\boldsymbol{\alpha}^{\rm T},\boldsymbol{\beta}^{\rm T},\boldsymbol{\eta}^{\rm T})^{\rm T}\in\Theta_{\zeta}$, where $\Theta_{\zeta}=\Theta_{\alpha}\times\Theta_{\beta}\times\Theta_{\eta}$ is the product space.

Because the indicator function ${\boldsymbol{1}}(U_{1}+{\boldsymbol{U}}_{2}^{\rm T}\boldsymbol{\eta}\geq 0)$ is not differentiable, it is natural to approximate it by a smooth function $K(u)$ satisfying

$\displaystyle\lim_{u\rightarrow+\infty}K(u)=1~{}\mbox{and}~{}\lim_{u\rightarrow-\infty}K(u)=0.$

Note that this smooth function characterizes the cumulative distribution function instead of the density function; see Seo and Linton, (2007); Li et al., (2021); Zhang et al., (2022); Mukherjee et al., (2020) for more details. The literature contains many commonly used smooth functions, such as the cumulative distribution function of standard normal distribution $K(u)=\Phi(u)$, the sigmoid function $K(u)=\{1+\exp(-u)\}^{-1}$, and the mixture of the cumulative distribution function and density of standard normal distribution $K(u)=\Phi(u)+u\phi(u)$. Thus, model (6) can be approximated by

$\displaystyle y_{i}={\boldsymbol{X}}_{i}^{\rm T}\boldsymbol{\alpha}+{\boldsymbol{Z}}_{i}^{\rm T}\boldsymbol{\beta}K_{h}(U_{1i}+{\boldsymbol{U}}_{2i}^{\rm T}\boldsymbol{\eta})+\epsilon_{i},$

(7)

where $K_{h}(u)=K(u/h)$, and $h$ is a predetermined tuning parameter associated with $n$ satisfying $\lim_{n\rightarrow\infty}h=0$ , called the smoothness parameter.

For any $\tau>0$, let $\boldsymbol{\zeta}_{\tau,h}^{*}$ be the minimizer defined as

$\displaystyle\boldsymbol{\zeta}_{\tau,h}^{*}=\mathop{\rm argmin}_{\boldsymbol{\zeta}\in\Theta_{\zeta}}\boldsymbol{\mathrm{P}}L_{\tau}(y-{\boldsymbol{X}}^{\rm T}\boldsymbol{\alpha}-{\boldsymbol{Z}}^{\rm T}\boldsymbol{\beta}K_{h}(U_{1}+{\boldsymbol{U}}_{2}^{\rm T}\boldsymbol{\eta})),$

(8)

which approximates the minimizer

$\displaystyle\boldsymbol{\zeta}_{\tau}^{*}=\mathop{\rm argmin}_{\boldsymbol{\zeta}\in\Theta_{\zeta}}\boldsymbol{\mathrm{P}}L_{\tau}(y-{\boldsymbol{X}}^{\rm T}\boldsymbol{\alpha}-{\boldsymbol{Z}}^{\rm T}\boldsymbol{\beta}{\boldsymbol{1}}(U_{1}+{\boldsymbol{U}}_{2}^{\rm T}\boldsymbol{\eta}\geq 0)).$

(9)

As did Sun et al., (2020), $\boldsymbol{\zeta}_{\tau}^{*}$ is called the Huber coefficient, which usually distinguishes from the true parameter $\boldsymbol{\zeta}^{*}$. Measured by $\|\boldsymbol{\zeta}_{\tau}^{*}-\boldsymbol{\zeta}^{*}\|$, the Huber error is caused by the robustification for the heavy-tailed errors, while the distance $\|\boldsymbol{\zeta}_{\tau,h}^{*}-\boldsymbol{\zeta}^{*}\|$ is a consequence of both robustification and smoothness. Theorem 2 reveals that $\|\boldsymbol{\zeta}_{\tau,h}^{*}-\boldsymbol{\zeta}^{*}\|$ is controlled by both $\tau$ and $h$, with $h$ playing the role of the bandwidth in the nonparametric area.

Minimizing the empirical loss in (8) produces the robust estimator of interest, i.e.,

$\displaystyle\hat{\boldsymbol{\zeta}}_{\tau,h}=\mathop{\rm argmin}_{\boldsymbol{\zeta}\in\Theta_{\zeta}}\mathbb{P}_{n}L_{\tau}(y-{\boldsymbol{X}}^{\rm T}\boldsymbol{\alpha}-{\boldsymbol{Z}}^{\rm T}\boldsymbol{\beta}K_{h}(U_{1}+{\boldsymbol{U}}_{2}^{\rm T}\boldsymbol{\eta})).$

(10)

From (8), (9), and (10), the total estimation error $\|\hat{\boldsymbol{\zeta}}_{\tau,h}-\boldsymbol{\zeta}^{*}\|$ can be decomposed into three parts, i.e.,

$\displaystyle\underbrace{\|\hat{\boldsymbol{\zeta}}_{\tau,h}-\boldsymbol{\zeta}^{*}\|}_{\mbox{Total error}}\leq\underbrace{\|\hat{\boldsymbol{\zeta}}_{\tau,h}-\boldsymbol{\zeta}_{\tau,h}^{*}\|}_{\mbox{estimation error}}+\underbrace{\|\boldsymbol{\zeta}_{\tau,h}^{*}-\boldsymbol{\zeta}_{\tau}^{*}\|}_{\mbox{smooth error}}+\underbrace{\|\boldsymbol{\zeta}_{\tau}^{*}-\boldsymbol{\zeta}^{*}\|}_{\mbox{Huber error}}.$

(11)

It is natural to use the alternating strategy to obtain the estimate, denoted by $\hat{\boldsymbol{\zeta}}=(\hat{\boldsymbol{\alpha}}^{\rm T},\hat{\boldsymbol{\beta}}^{\rm T},\hat{\boldsymbol{\eta}}^{\rm T})^{\rm T}$. Specifically, the parameters $(\boldsymbol{\alpha}^{\rm T},\boldsymbol{\beta}^{\rm T})^{\rm T}$ and $\boldsymbol{\eta}$ can be estimated iteratively as follows. For given $\boldsymbol{\eta}^{(k)}$, $\boldsymbol{\alpha}^{(k+1)}$ and $\boldsymbol{\beta}^{(k+1)}$ are obtained by minimizing

$\displaystyle\mathbb{P}_{n}L_{\tau}(y-{\boldsymbol{X}}^{\rm T}\boldsymbol{\alpha}-{\boldsymbol{Z}}^{\rm T}\boldsymbol{\beta}K_{h}(U_{1}+{\boldsymbol{U}}_{2}^{\rm T}\boldsymbol{\eta}^{(k)})),$

and for given $\boldsymbol{\alpha}^{(k+1)}$ and $\boldsymbol{\beta}^{(k+1)}$, $\boldsymbol{\gamma}^{(k+1)}$ is estimated by minimizing

$\displaystyle\mathbb{P}_{n}L_{\tau}(y-{\boldsymbol{X}}^{\rm T}\boldsymbol{\alpha}^{(k+1)}-{\boldsymbol{Z}}^{\rm T}\boldsymbol{\beta}^{(k+1)}K_{h}(U_{1}+{\boldsymbol{U}}_{2}^{\rm T}\boldsymbol{\eta})).$

Iterating these two maximizers leads to the desired robust estimator. The above alternating strategy is summarized in Algorithm A in Appendix A of the Supplementary Material.

This section begins with assumptions needed to establish the nonasymptotic properties. Let ${\bar{\boldsymbol{V}}}$ be ${\boldsymbol{V}}$ by removing $U_{1}$, i.e., ${\bar{\boldsymbol{V}}}=(y,{\boldsymbol{X}},{\boldsymbol{Z}},{\boldsymbol{U}}_{2}))$, and $\bar{\omega}(\boldsymbol{\eta})=U_{1}+{\boldsymbol{U}}_{2}^{\rm T}\boldsymbol{\eta}$ and $\bar{\omega}=\bar{\omega}(\boldsymbol{\eta}^{*})$.

1. (A1)

   The conditional random vectors ${\mathrm{E}}({\boldsymbol{X}}|{\boldsymbol{U}}_{2})\in\mathbb{R}^{p}$ and ${\mathrm{E}}({\boldsymbol{Z}}|{\boldsymbol{U}}_{2})\in\mathbb{R}^{q}$ given ${\boldsymbol{U}}_{2}$ are sub-Gaussian, and ${\boldsymbol{U}}_{2}$ is sub-Gaussian. There is a universal constant $K_{1}>0$ satisfying $\|{\mathrm{E}}({\boldsymbol{X}}|{\boldsymbol{U}}_{2})\|_{\psi_{2}}\leq K_{1}$, $\|{\mathrm{E}}({\boldsymbol{Z}}|{\boldsymbol{U}}_{2})\|_{\psi_{2}}\leq K_{1}$, and $\|{\boldsymbol{U}}_{2}\|_{\psi_{2}}\leq K_{1}$. For any ${\boldsymbol{u}}_{2}\in\mathbb{R}^{r-1}$, the matrices ${\mathrm{E}}({\boldsymbol{X}}{\boldsymbol{X}}^{\rm T}|{\boldsymbol{U}}_{2}={\boldsymbol{u}}_{2})$, ${\mathrm{E}}({\boldsymbol{X}}{\boldsymbol{X}}^{\rm T}|{\boldsymbol{U}}_{2}={\boldsymbol{u}}_{2})$, and ${\mathrm{E}}({\boldsymbol{U}}_{2}{\boldsymbol{U}}_{2}^{\rm T})$ are uniformly positive definite, and there is a universal constant $K_{0}>0$ satisfying $\lambda_{\min}({\mathrm{E}}({\boldsymbol{X}}{\boldsymbol{X}}^{\rm T}|{\boldsymbol{U}}_{2}={\boldsymbol{u}}_{2}))\geq K_{0}$, $\lambda_{\min}({\mathrm{E}}({\boldsymbol{Z}}{\boldsymbol{Z}}^{\rm T}|{\boldsymbol{U}}_{2}={\boldsymbol{u}}_{2}))\geq K_{0}$, and $\lambda_{\min}({\mathrm{E}}({\boldsymbol{U}}_{2}{\boldsymbol{U}}_{2}^{\rm T}))\geq K_{0}$.

2. (A2)

   The error variable $\epsilon$ is independent of $({\boldsymbol{X}}^{\rm T},{\boldsymbol{Z}}^{\rm T},{\boldsymbol{U}}^{\rm T})^{\rm T}$ and satisfies ${\mathrm{E}}(\epsilon)=0$ and ${\mathrm{E}}(|\epsilon|^{2+\delta})=M_{\delta}<\infty$ with $\delta\geq 0$. Denote $M_{0}={\mathrm{E}}(|\epsilon|^{2})$.

3. (A3)

   $0<{\mathrm{E}}[{\boldsymbol{1}}({\boldsymbol{U}}^{\rm T}\boldsymbol{\gamma}\geq 0)]<1$ for any $\boldsymbol{\gamma}\in\Theta_{\gamma}$, and there is a constant $\delta_{u_{2}}>0$ satisfying $\sup_{\boldsymbol{\eta}\in\mathbb{S}^{r-1}}{\mathrm{E}}[{\boldsymbol{U}}_{2}^{\rm T}\boldsymbol{\eta}{\boldsymbol{1}}({\boldsymbol{U}}_{2}^{\rm T}\boldsymbol{\eta}\geq 0)]\geq\delta_{u_{2}}$.

4. (A4)

   For almost every ${\boldsymbol{u}}_{2}$, the density of $U_{1}$ conditional on ${\boldsymbol{U}}_{2}={\boldsymbol{u}}_{2}$ is everywhere positive. The conditional density $f_{\varpi|\bar{{\boldsymbol{v}}}}(\varpi)$ of $\bar{\omega}$ given ${\bar{\boldsymbol{V}}}$ has continuous derivative, and there is a constant $\kappa_{f}>0$ such that $f_{\varpi|\bar{{\boldsymbol{v}}}}(\varpi)$ and $|f_{\varpi|\bar{{\boldsymbol{v}}}}^{\prime}(\varpi)|$ are uniformly bounded from above by $\kappa_{f}$ over $(\varpi,\bar{{\boldsymbol{v}}})$. $f_{\varpi|\bar{{\boldsymbol{v}}}}(0)$ is uniformly bounded from below by $\delta_{f_{0}}$ over $\bar{{\boldsymbol{v}}}$, and $F_{\varpi|\bar{{\boldsymbol{v}}}}(0)$ is uniformly bounded from above by $\kappa_{F}$ over $\bar{{\boldsymbol{v}}}$, where $\delta_{f_{0}}>0$ and $0<\kappa_{F}<1$ are constants.

5. (A5)

   The smooth function $K(\cdot)$ is twice differentiable and $K(-t)=1-K(t)$. $K^{\prime}(\cdot)$ is symmetric around zero. Moveover, there is a universal constant $\kappa_{k}>0$ satisfying $\max\left\{\sup_{t\in\mathbb{R}}|K^{\prime}(t)|,\sup_{t\in\mathbb{R}}|K^{\prime\prime}(t)|,\int|K^{\prime}(t)|^{j}dt,\int|K^{\prime\prime}(t)|^{j}dt,\int|t||K^{\prime}(t)|^{j}dt\right\}\leq\kappa_{k}$ with $j=1,2$.

Theorem 1 states that the Huber error is of order $\tau^{-(1+\delta)}$ for $\|\boldsymbol{\theta}_{\tau}^{*}-\boldsymbol{\theta}^{*}\|$ and of order $\tau^{-(2+2\delta)}$ for $\|\boldsymbol{\eta}^{*}-\boldsymbol{\eta}_{\tau}^{*}\|$. As $\tau$ tends to infinity, because the Huber loss becomes the ordinary quadratic loss, the Huber error vanishes as expected. The next theorem studies the smooth error and Huber loss together.

Theorem 2 reveals that $\|\boldsymbol{\theta}_{\tau,h}^{*}-\boldsymbol{\theta}^{*}\|$ and $\|\boldsymbol{\eta}_{\tau,h}^{*}-\boldsymbol{\eta}^{*}\|$ are associated with both the Huber error and the smooth error. The smoothness parameter $h$ can be of order $\tau^{-(1+\delta)}$, and because $K_{h}(t)$ approximates ${\boldsymbol{1}}(t\geq 0)$ as $h$ tends to zero, the upper bounds in Theorem 2 are of the same order as those in Theorem 1 when $h\rightarrow 0$. Thus, the deviations $\|\boldsymbol{\theta}_{\tau,h}^{*}-\boldsymbol{\theta}^{*}\|$ and $\|\boldsymbol{\eta}_{\tau,h}^{*}-\boldsymbol{\eta}^{*}\|$ are sacrifices in pursuit of robustification and smoothness. The next theorem provides the exponential-type deviations for the baseline and grouping difference parameters and grouping parameter.

With appropriate choice of $\tau$ with $\delta=0$, such as $\tau=O((n/t)^{1/2})$ with $t=\log(n)$, the nonasymptotic property of the sub-Gaussian estimator $\hat{\boldsymbol{\theta}}_{\tau,h}$ in Theorem 3 demonstrates that the deviation $\|\hat{\boldsymbol{\theta}}_{\tau,h}-\boldsymbol{\theta}^{*}\|$ adapts to the sample size, dimension, robustification parameter $\tau$, and moments in pursuit of the optimal tradeoff between bias and robustness, and the deviation $\|\hat{\boldsymbol{\eta}}_{\tau,h}-\boldsymbol{\eta}^{*}\|$ adapts to the extra smooth parameter $h$ to achieve smoothness. Adaptation to the robustification parameter $\tau$ is caused by pursuing the robustness for linear regression with heavy-tailed errors. The deviations $\|\hat{\boldsymbol{\theta}}_{\tau,h}-\boldsymbol{\theta}^{*}\|$ and $\|\hat{\boldsymbol{\eta}}_{\tau,h}-\boldsymbol{\eta}^{*}\|$ coincide with the smoothed OLS estimator when $\tau\rightarrow\infty$, because $a(n,\tau)=\sqrt{3\nu_{0}(p+q+r+2t)/n}$ as $\tau\rightarrow\infty$. The next theorem derives the nonasymptotic Bahadur representation of the robust estimators given in (10).

Theorem 4 shows that the remainder of the Bahadur representation of $\|\hat{\boldsymbol{\theta}}_{\tau,h}-\boldsymbol{\theta}^{*}\|$ achieves the rate $h^{1/2}a(n,\tau)$, which is the same as that for $\|\hat{\boldsymbol{\eta}}_{\tau,h}-\boldsymbol{\eta}^{*}\|$. Because of the rate restriction $h^{2}n/(p+q+r)\rightarrow\infty$ and $h=o(1)$, the remainder of the Bahadur representation in inequality (15) does not exhibit subexponential behavior as considered by Sun et al., (2020); Chen and Zhou, (2020). This reason is that there is a change plane involved a smooth function in model (1). To the best of the authors’ knowledge, this is the first time that this type of nonasymptotic Bahadur representation has been reported in the literature, especially for the robust estimator $\hat{\boldsymbol{\eta}}_{\tau,h}$ of the grouping parameter, with previous studies reporting only polynomial-type deviation bounds; see Liu et al., (2024); Zhang et al., (2022). It is convenient to derive the classical asymptotic results from the Bahadur representation, and Theorems 3 and 4 show that the robustification parameter $\tau$ and the smoothness parameter $h$ play the same role as bandwidth in constructing classical nonparametric estimators.

The theoretical properties in Section 2.2 guarantee that the robust estimation performs well with appropriate choices of the robustification parameter $\tau$ and the smoothness parameter $h$. Because the robustification parameter $\tau$ is treated as a tuning parameter to balance bias and robustness, it is natural to consider using the cross-validation (CV) method to select an appropriate $\tau$ in practice. However, as noted by Chen and Zhou, (2020) and Wang et al., (2021), because $M_{\delta}=\textrm{E}(|\epsilon|^{2+\delta})$ is typically unknown in practice, its empirical OLS estimator $\hat{M}_{\delta}=(n-p-q)^{-1}\sum_{i=1}^{n}(y_{i}-{\boldsymbol{X}}_{i}^{\rm T}\hat{\boldsymbol{\alpha}}_{\tau,h}-{\boldsymbol{Z}}_{i}^{\rm T}\hat{\boldsymbol{\beta}}_{\tau,h}{\boldsymbol{1}}({\boldsymbol{U}}_{i}^{\rm T}\hat{\gamma}_{\tau,h}\geq 0))^{2+\delta}$ is poor when the errors are heavy-tailed. Instead, there are two good alternatives (Chen and Zhou, 2020): (i) an adaptive technique based on Lepski’s method Lepskii, (1992) and (ii) a Huber-type method by solving a so-called censored equation Hahn et al., (1990); see Chen and Zhou, (2020); Wang et al., (2021) for details. For the smoothness parameter $h$, the CV method is always a natural choice for selecting an appropriate one. Alternatively, from the theoretical conditions on the smoothness parameter $h$, a rule of thumb $h_{n}=c_{h}\hat{\sigma}_{u}\log(n)/\sqrt{n}$ suggested by Seo and Linton, (2007), Li et al., (2021), and Zhang et al., (2022) is recommended by considering the computation reduction, where $c_{h}$ is a constant and $\hat{\sigma}_{u}=\sqrt{(n-r)^{-1}\sum_{i=1}^{n}(U_{1i}+{\boldsymbol{U}}_{2i}^{\rm T}\hat{\boldsymbol{\eta}}^{ols})^{2}}$ is the estimated standard deviation of ${\boldsymbol{U}}^{\rm T}\boldsymbol{\gamma}$, with $(\hat{\boldsymbol{\alpha}}^{ols},\hat{\boldsymbol{\beta}}^{ols},\hat{\boldsymbol{\eta}}^{ols})$ being the estimator of $(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\eta})$ by using ordinary quadratic loss.

Attention now turns to the implementation for subgroup-classifier learning, with the parameters $\boldsymbol{\theta}=(\boldsymbol{\alpha}^{\rm T},\boldsymbol{\beta}^{\rm T})^{\rm T}$ and $\boldsymbol{\eta}$ estimated iteratively as follows. Let $\ell_{\tau}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\eta})$ be the loss function for model (1), i.e.,

$\displaystyle\ell_{\tau}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\eta})=\sum_{i=1}^{n}L_{\tau}\left(y_{i}-{\boldsymbol{X}}_{i}^{\rm T}\boldsymbol{\alpha}-{\boldsymbol{Z}}_{i}^{\rm T}\boldsymbol{\beta}{\boldsymbol{1}}(U_{1i}+{\boldsymbol{U}}_{2i}^{\rm T}\boldsymbol{\eta}\geq 0)\right),$

and let the smoothed loss function be

$\displaystyle\tilde{\ell}_{\tau,h}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\eta})=\sum_{i=1}^{n}L_{\tau}\left(y_{i}-{\boldsymbol{X}}_{i}^{\rm T}\boldsymbol{\alpha}-{\boldsymbol{Z}}_{i}^{\rm T}\boldsymbol{\beta}K_{h}(U_{1i}+{\boldsymbol{U}}_{2i}^{\rm T}\boldsymbol{\eta}\geq 0)\right).$

(16)

For given $\boldsymbol{\eta}^{(k)}$, one obtains $\boldsymbol{\alpha}^{(k+1)}$ and $\boldsymbol{\beta}^{(k+1)}$ by minimizing the smoothed loss function

$\displaystyle(\boldsymbol{\alpha}^{(k+1)},\boldsymbol{\beta}^{(k+1)})=\mathop{\rm argmin}_{\boldsymbol{\alpha}\in\Theta_{\alpha},\boldsymbol{\beta}\in\Theta_{\beta}}\tilde{\ell}_{\tau,h}(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\eta}^{(k)}),$

and for given $\boldsymbol{\alpha}^{(k+1)}$ and $\boldsymbol{\beta}^{(k+1)}$, one estimates $\boldsymbol{\eta}^{(k+1)}$ by

$\displaystyle\boldsymbol{\eta}^{(k+1)}=\mathop{\rm argmin}_{\boldsymbol{\eta}\in\Theta_{\eta}}\tilde{\ell}_{\tau,h}(\boldsymbol{\alpha}^{(k+1)},\boldsymbol{\beta}^{(k+1)},\boldsymbol{\eta}).$

Iterating these two minimizers leads to the desired robust estimators. These are summarized with the multiplier bootstrap calibration in Algorithm 1 in Appendix A of the Supplementary Material, which provides the strategy for estimating the confidence intervals for the estimators $\hat{\boldsymbol{\alpha}}$, $\hat{\boldsymbol{\beta}}$, and $\hat{\boldsymbol{\eta}}$.

Note that herein, $\boldsymbol{\alpha}^{(k+1)}$ and $\boldsymbol{\beta}^{(k+1)}$ for given $\boldsymbol{\eta}^{(k)}$ are obtained by a robust data-adaptive method proposed by Wang et al., (2021). Specifically, $\boldsymbol{\theta}=(\boldsymbol{\alpha}^{\rm T},\boldsymbol{\beta}^{\rm T})^{\rm T}$ is estimated and $\tau$ is calibrated simultaneously by solving the following system of equations:

$\displaystyle\left\{\begin{array}[]{l}\sum_{i=1}^{n}\psi_{\tau}(y_{i}-{\boldsymbol{W}}_{i}^{\rm T}\boldsymbol{\theta}){\boldsymbol{W}}_{i}=0,\\ (\tau^{2}n)^{-1}\sum_{i=1}^{n}\min\{(y_{i}-{\boldsymbol{W}}_{i}^{\rm T}\boldsymbol{\theta})^{2},\tau^{2}\}-n^{-1}(d+z)=0,\end{array}\right.$

where $d=p+q-1$, ${\boldsymbol{W}}_{i}=({\boldsymbol{X}}_{i}^{\rm T},{\boldsymbol{Z}}_{i}^{\rm T}{\boldsymbol{1}}(U_{1i}+{\boldsymbol{U}}_{2i}^{\rm T}\boldsymbol{\eta}^{(k)}\geq 0))^{\rm T}$, and $z=\log(n)$ as suggested by Wang et al., (2021). The initial values $\boldsymbol{\theta}^{(0)}=\boldsymbol{\theta}^{(ols)}$ and $\tau^{(0)}=\hat{\sigma}_{\epsilon}\sqrt{n/(d+z)}$ are set using the ordinary quadratic loss, where $\hat{\sigma}_{\epsilon}^{2}=(n-p-q-r)^{-1}\sum_{i=1}^{n}(y-{\boldsymbol{X}}_{i}^{\rm T}\hat{\boldsymbol{\alpha}}^{ols}-{\boldsymbol{Z}}_{i}^{\rm T}\hat{\boldsymbol{\beta}}^{ols}{\boldsymbol{1}}(U_{1i}+{\boldsymbol{U}}_{2i}^{\rm T}\hat{\boldsymbol{\eta}}^{ols}\geq 0))^{2}$, with $(\hat{\boldsymbol{\alpha}}^{ols},\hat{\boldsymbol{\beta}}^{ols},\boldsymbol{\eta}^{ols})$ being the estimator of $(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\eta})$.

Under the null hypothesis, model (1) reduces to the ordinary linear regression model with heavy-tailed errors, i.e.,

$\displaystyle y_{i}={\boldsymbol{X}}_{i}^{\rm T}\boldsymbol{\alpha}+\epsilon_{i}.$

(18)

Parametric estimation in model (18) is well-studied in the literature; see Huber, (1964, 1973); Fan et al., 2017b ; Sun et al., (2020); Wang et al., (2021); Han et al., (2022); Chen and Zhou, (2020); Zhou et al., (2018), among others. Let $\hat{\boldsymbol{\alpha}}_{\tau}$ be the estimate of $\boldsymbol{\alpha}_{\tau}$ under the null hypothesis, i.e.,

$\displaystyle\hat{\boldsymbol{\alpha}}_{\tau}=\mathop{\rm argmin}_{\boldsymbol{\alpha}\in\Theta_{\alpha}}\mathbb{P}_{n}L_{\tau}(y-{\boldsymbol{X}}^{\rm T}\boldsymbol{\alpha}).$

(19)