Variable Selection for Additive Quantile Regression with Nonlinear Interaction Structures

Abstract

In high-dimensional data analysis, main effects and interaction effects often coexist, especially when complex nonlinear relationships are present. Effective variable selection is crucial for avoiding the curse of dimensionality and enhancing the predictive performance of a model. In this paper, we introduce a nonlinear interaction structure into the additive quantile regression model and propose an innovative penalization method. This method considers the complexity and smoothness of the additive model and incorporates heredity constraints on main effects and interaction effects through an improved regularization algorithm under marginality principle. We also establish the asymptotic properties of the penalized estimator and provide the corresponding excess risk. Our Monte Carlo simulations illustrate the proposed model and method, which are then applied to the analysis of Parkinson’s disease rating scores and further verify the effectiveness of a novel Parkinson’s disease (PD) treatment.

1. Introduction

Quantile regression [1] has become a widely used tool in both empirical studies and the theoretical analysis of conditional quantile functions. Additive models offer greater flexibility by expressing linear predictors as the sum of nonparametric functions of each covariate, which generally results in lower variance compared to fully nonparametric models. However, many practical problems require the consideration of interactions between covariates, an aspect that has been extensively explored in linear and generalized linear models [2,3], where incorporating interactions has been shown to improve prediction accuracy. Despite this, there is limited research on additive quantile regression models that explicitly account for interaction structures, particularly in the context of variable selection.

While nonparametric additive quantile regression models have been extensively studied in the literature [4,5,6,7,8], with recent advancements in the analysis of longitudinal data [9] and dynamic component modeling [10]. Specifically, Ref. [9] explored a partially linear additive model for longitudinal data within the quantile regression framework, utilizing quadratic inference functions to account for within-subject correlations. In parallel, Ref. [10] proposed a quantile additive model with dynamic component functions, introducing a penalization-based approach to identify non-dynamic components. Despite these advances, the integration of interaction effects remains relatively underexplored in the context of nonparametric additive quantile regression models. Building upon the existing literature, this paper considers an additive quantile regression model with a nonlinear interaction structure.

When p, the number of main effect terms and interaction effect terms, is large, many additive methods are not feasible since their implementation requires storing and manipulating the entire $O\left(p^2\right)\times n$ design matrix; thus, variable selection becomes necessary. Refs. [11,12] proposed component selection and smoothing operator (COSSO) and an adaptive version of the COSSO algorithm (ACOSSO) to fit the additive model with interaction structures, respectively. But these methods violated the heredity condition, whereas [13] took the heredity constraint into account and proposed a penalty $\lambda ({\sum }_{j=1}^p\parallel f_j{\parallel }_2+{\sum }_{j=1}^p{\sum }_{k=j+1}^p\parallel f_{jk}{\parallel }_2)$ on the empirical $L_2$ norm of the main effects and interaction effects. This method shrinks the interactions, depending on the main effects that are already present in the model [9,10].

Several regularization penalties, especially Lasso, have been widely used to shrink the coefficients of main effects and interaction effects to achieve a sparse model. However, existing methods treated interaction terms and main terms similarly, which may select an interaction term but not the corresponding main terms, thus making them difficult to interpret in practice. There is a natural hierarchy among the variables in the model with interaction structures. Refs. [14,15] proposed a two-stage Lasso method to select important main and interaction terms. However, this approach is inefficient, as the solution path in the second stage heavily depends on the selection results from the first stage. To address this issue, several regularization methods [16,17,18] have been proposed that employ special reparametrizations of regression coefficients and penalty functions under the hierarchical principle. Ref. [19] introduced strong and weak heredity constraints into models with interactions. They proposed a Lasso penalty with convex constraints that produces sparse coefficients while ensuring that the strong or weak heredity is satisfied. However, existing algorithms often suffer from slow convergence, even with a moderate number of predictors. Ref. [20] proposed a group-regularized estimation method under both strong and weak heredity constraints and developed a computational algorithm that guarantees the convergence of the iterates. Ref. [21] extended the linear interaction model to a quadratic model, accounting for all two-term interactions, and proposed the Regularization Path Algorithm under the Marginality Principle (RAMP), which efficiently computes the regularization solution path while maintaining strong or weak heredity constraints.

Recent advances in Bayesian hierarchical modeling for interaction selection have been driven by innovative prior constructions. Ref. [22] developed mixture priors that link interaction inclusion probabilities to main effect strengths, automatically enforcing heredity constraints. Ref. [23] introduced structured shrinkage priors using hierarchical Laplace distributions to facilitate group-level interaction selection. Ref. [24] incorporated heredity principles into SSVS frameworks through carefully designed spike-and-slab priors. These methodological developments have enhanced the ability to capture complex dependency structures in interaction selection while preserving Bayesian advantages, with demonstrated applications across various domains, including spatiotemporal analysis.

However, all these works assumed that the effects of all covariates can be captured in a simple linear form, which does not always hold in practice. To handle this issue, several authors have proposed nonparametric and semiparametric interaction model. For example, Refs. [25,26] considered the variable selection and estimation procedure of single-index models with interactions for functional data and longitudinal data. Some recent work on semiparametric single-index models with interactions included [27,28], among others.

Driven by both practical needs and theoretical considerations, this paper makes the following contributions: (1) We propose a method for variable selection in quantile regression models that incorporates nonlinear interaction structures; (2) We modify the RAMP algorithm to extend its applicability to additive models with nonlinear interactions, which enhances its capability to adeptly handle complex nonlinear interactions while also preserving the hierarchical relationship between main effects and interaction effects throughout the variable selection process. This work addresses the challenges posed by complex nonlinear interactions while ensuring that the inherent hierarchical structure is maintained, thereby providing a robust framework for more accurate and interpretable models.

The rest of this paper is organized as follows. In Section 2, we present the additive quantile regression with nonlinear interaction structures and discuss the properties of the oracle estimator. In Section 3, we present a sparsity smooth penalized method fitted by regularization algorithm under marginality principle and derive its oracle property. In Section 4, simulation studies are provided to demonstrate the outperformance of the proposed method. We illustrate an application of the proposed method on a real dataset in Section 5, and a conclusion is presented in Section 6.

2. Additive Quantile Regression with Nonlinear Interaction Structures

Suppose that $\{(Y_i,{\mathbf{x}}_i):i=1,\dots ,n\}$ is an independent and identically distributed (iid) sample, where ${\mathbf{x}}_i=(x_{i1},\dots ,x_{ip})$ is a p-dimensional vector of covariates. The $\tau$th ($0<\tau <1$) conditional quantile of $Y_i$ given ${\mathbf{x}}_i$ is defined as $Q_{Y_i|{\mathbf{x}}_i}\left(\tau \right)=inf\{t:F\left(t\right|{\mathbf{x}}_i)\ge \tau \}$, where $F(·|{\mathbf{x}}_i)$ is the conditional distribution function of $Y_i$ given ${\mathbf{x}}_i$. We consider the following additive nonlinear interaction model for the conditional quantile function:

(1)$Q_{Y_i|{\mathbf{x}}_i}\left(\tau \right)=\sum _{j=1}^p{f}_j\left(x_{ij}\right)+\sum _{1\le j<k\le p}{f}_{jk}(x_{ij},x_{ik}),$

where the unknown real-valued two-variate functions $f_{jk}(x_{ij},x_{ik})$ represent the pairwise interaction terms between $x_{ij}$ and $x_{ik}$. We assume that $f_{jk}(a,b)=f_{kj}(b,a)$ for all a and b and that all $j\ne k$ and $x_{ij}\in [0,1]$ for all i and j. Let ${ϵ}_i=Y_i-Q_{Y_i|{\mathbf{x}}_i}\left(\tau \right)$; then, ${ϵ}_i$ satisfies $P({ϵ}_i\le 0|{\mathbf{x}}_i)=\tau$ and we may also write $Y_i={\displaystyle \sum _{j=1}^p}{f}_j\left(x_{ij}\right)+{\displaystyle \sum _{1\le j<k\le p}}{f}_{jk}(x_{ij},x_{ik})+{ϵ}_i$, where ${ϵ}_i$ is the random error. For identification, we assume the $\tau$ quantile of $f_j\left(x_{ij}\right)$ and $f_{jk}(x_{ij},x_{ik})$ for $1\le j<k\le p$ to be zero (see [29]).

Let $\{{\phi }_1,{\phi }_2,\cdots \}$ denote the preselected orthonormal basis with respect to the Lebesgue measure on the unit interval. We center all the candidate functions, so the basis functions are centered as well. We approximate each function $f_j\left(x_{ij}\right)$ by a linear combination of the basis functions, i.e.,

$f_j\left(x_{ij}\right)\approx \sum _{l=1}^{L_n}{\phi }_{jl}\left(x_{ij}\right){\mathbf{\beta }}_{jl},$

where $\{{\phi }_{jl},l=1,\dots ,L_n\}$s are centered orthonormal bases. Let ${\mathsf{\Psi }}_{i,j}={({\phi }_{j1}\left(x_{ij}\right),\dots ,{\phi }_{j{L}_n}\left(x_{ij}\right))}^{\top }$; the main terms can be expressed as $f_j\left(x_{ij}\right)={\mathbf{\Psi }}_{i,j}^{\top }{\beta }_j$, where ${\mathit{\beta }}_j={({\beta }_{j1},\cdots ,{\beta }_{j{L}_n})}^{\top }$ denotes the $L_n$-dimensional vector of basis coefficients for the j-th main terms. We consider the tensor product of the interaction terms $f_{jk}$ as a surface that can be approximated by

$f_{jk}(x_{ij},x_{ik})\approx \sum _{s=1}^{L_n}\sum _{d=1}^{L_n}{\phi }_{js}\left(x_{ij}\right){\phi }_{kd}\left(x_{ik}\right){\beta }_{jk,sd}.$

To simplify the expression, it is computationally efficient to re-express the surface in vector notation as $f_{jk}(x_{ij},x_{ik})={\mathsf{\Phi }}_{jk}^{\top }{\mathit{\beta }}_{jk}$, where ${\mathsf{\Phi }}_{i,jk}={\mathbf{\Psi }}_{i,j}\otimes {\mathbf{\Psi }}_{i,k}$, and ${\mathit{\beta }}_{jk}=vec\left({\mathsf{\Gamma }}_{jk}^{\top }\right)$, where ${\mathsf{\Gamma }}_{jk}=\left[{\beta }_{jk,sd}\right],s=1,\dots ,L_n,d=1,\dots ,L_n$ is the matrix with elements ${\beta }_{jk,sd}$. Here, ⊗ denotes the Kronecker product. We can now express (1) as

(2)$Q_{Y_i|{\mathbf{x}}_i}\left(\tau \right)=\sum _{j=1}^p{\mathbf{\Psi }}_{i,j}^{\top }{\mathit{\beta }}_j+\sum _{1\le j<k\le p}{\mathsf{\Phi }}_{i,jk}^{\top }{\mathit{\beta }}_{jk}$

2.1. Oracle Estimator

For high-dimensional inference, it is often assumed that p is large, but only a small fraction of the main effects and interaction terms are present in the true model. Let $\mathcal{M}=\{j:f_j\ne 0\mathrm{for}1\le j\le p\}$ and $\mathcal{I}=\{(j,k):1\le j<k\le p\mathrm{and}{f}_{jk}\ne 0\}$ be the index set of nonzero main effects and interaction effects, and $q=\left|\mathcal{M}\right|$ and $s=\left|\mathcal{I}\right|$ be the cardinality of $\mathcal{M}$ and $\mathcal{I}$, respectively. That is, we assume that the first q of ${\mathit{\beta }}_j,j=1,\dots ,p$ are nonzero and the remaining components are zero. Hence, we can write ${\mathit{\beta }}_0={({\mathit{\beta }}_{01}^{\top },{\mathbf{0}}_{(p-q)L_n}^{\top })}^{\top }$, where ${\mathit{\beta }}_{01}={({\mathit{\beta }}_1^{\top },\dots ,{\mathit{\beta }}_q^{\top })}^{\top }$. Further, we assume that the first s of pairwise interaction ${\mathit{\beta }}_{jk}$ corresponding to ${\mathit{\beta }}_{01}$ are nonzero and the remaining components are zero. Thus, we write ${\mathit{\beta }}_{00}={({\mathit{\beta }}_{001}^{\top },{\mathbf{0}}_{{\displaystyle \frac{p(p-1)-s(s-1)}{2}}{L}_n^2}^{\top })}^{\top }$, where ${\mathit{\beta }}_{001}={({\mathit{\beta }}_{\left(12\right)}^{\top },\dots ,{\mathit{\beta }}_{\left(1s\right)}^{\top },\dots ,{\mathit{\beta }}_{\left(\right(s-1\left)1\right)}^{\top },\dots ,{\mathit{\beta }}_{\left(\right(s-1\left)s\right)}^{\top })}^{\top }$. Here, ${\mathit{\beta }}_{\left(jk\right)}$ denotes the component of ${\mathit{\beta }}_{jk}$ with $\left\{\right(j,k):(j,k)\in \mathcal{I},j\in \mathcal{M},k\in \mathcal{M}\}$ or $\left\{\right(j,k):(j,k)\in \mathcal{I};j\in \mathcal{M},\mathrm{or}k\in \mathcal{M}\}$. Let $\mathsf{\Pi }={({\mathsf{\Pi }}_1,\dots ,{\mathsf{\Pi }}_n)}^{\top }$ be an n by $p{L}_n+(p(p-1)/2)L_n^2$ design matrix with ${\mathsf{\Pi }}_i={({\mathbf{\Psi }}_i^{\top },{\mathsf{\Phi }}_i^{\top })}^{\top }$, where ${\mathbf{\Psi }}_i={({\mathbf{\Psi }}_{i,1}^{\top },\dots ,{\mathbf{\Psi }}_{i,p}^{\top })}^{\top }\in R^{p{L}_n}$ and ${\mathsf{\Phi }}_i={({\mathsf{\Phi }}_{i,12}^{\top },\dots ,{\mathsf{\Phi }}_{i,1p}^{\top },\dots ,{\mathsf{\Phi }}_{i,(p-1)1}^{\top },\dots ,{\mathsf{\Phi }}_{i,(p-1)p}^{\top })}^{\top }\in R^{p(p-1)/2L_n^2}$. Likewise, we write ${\mathsf{\Pi }}_{{\mathcal{A}}_i}={({\mathbf{\Psi }}_{{\mathcal{M}}_i}^{\top },{\mathsf{\Phi }}_{{\mathcal{I}}_i}^{\top })}^{\top }\in R^{q{L}_n+s{L}_n^2}$, where ${\mathbf{\Psi }}_{{\mathcal{M}}_i}$ is the subvector consisting of the first $q{L}_n$ elements of ${\mathbf{\Psi }}_i$ with the active covariates and ${\mathsf{\Phi }}_{{\mathcal{I}}_i}$ is the subvector consisting of the first $s{L}_n^2$ elements of ${\mathsf{\Phi }}_i$ with the active interaction covariates. We first investigate the estimator, which is obtained when the index sets $\mathcal{M}$ and $\mathcal{I}$ are known in advance, which we refer to as the oracle estimator. Now, we consider the quantile regression with the oracle information. Let

(3)$\left({\widehat{\mathit{\beta }}}_{01},{\widehat{\mathit{\beta }}}_{001}\right)=arg\underset{{\mathit{\beta }}_{01},{\mathit{\beta }}_{001}}{min}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\rho }_{\tau }\left(Y_i-{\mathbf{\Psi }}_{{\mathcal{M}}_i}^{\top }{\mathit{\beta }}_{01}-{\mathsf{\Phi }}_{{\mathcal{I}}_i}^{\top }{\mathit{\beta }}_{001}\right).$

The oracle estimators for ${\mathit{\beta }}_0$ and ${\mathit{\beta }}_{00}$ are ${\left({\widehat{\mathit{\beta }}}_{01}^{\top },{\mathbf{0}}_{(p-q)L_n}^{\top }\right)}^{\top }$ and ${\left({\widehat{\mathit{\beta }}}_{001},{\mathbf{0}}_{{\displaystyle \frac{p(p-1)-s(s-1)}{2}}{L}_n^2}^{\top }\right)}^{\top }$, respectively. Accordingly, the oracle estimators for nonparametric function $f_j\left(x_{ij}\right)$ and $f_{ij}(x_{ij},x_{ik})$ are ${\widehat{f}}_j\left(x_{ij}\right)={\mathbf{\Psi }}_{{\mathcal{M}}_i}^{\top }{\widehat{\mathit{\beta }}}_{01}$ and ${\widehat{f}}_j(x_{ij},x_{ik})={\mathsf{\Phi }}_{{\mathcal{I}}_i}^{\top }{\widehat{\mathit{\beta }}}_{001}$, respectively.

2.2. Asymptotic Properties

We next present the asymptotic properties of the oracle estimators. Similar to [13], we assume that all true main effects belong to the Sobolev space of order two: ${\sum }_{l=1}^{\infty }{\left({\beta }_{jl}\right)}^2{l}^4<C^2$, and impose the same requirement on each true interaction function ${\sum }_{s=1}^{\infty }{\left({\beta }_{jk,sd}\right)}^2{s}^4<C^2$ and ${\sum }_{d=1}^{\infty }{\left({\beta }_{jk,sd}\right)}^2{d}^4<C^2$ for each $j,k,1\le j<k\le p$, with C being some constant.

To establish the asymptotic properties of the estimators, the following regularity conditions are needed:

(A1). The distribution of $\mathbf{x}$ is absolutely continuous with density $f\left(x\right)$, which is bounded away from 0;

(A2). The conditional distribution $F_{\epsilon |\mathbf{x}}(·|x)$ of the error $\epsilon$, given $\mathbf{x}=x$, has a density $f_{\epsilon |\mathbf{x}}(·|x)$, which satisfies the following two conditions:

(1). $\underset{\epsilon ,x}{sup}{f}_{\epsilon |\mathbf{x}}\left(\epsilon \right|x)<\infty$;

(2). There exists positive constants $b_1$ and $b_2$ such that $\underset{x}{inf}\underset{\left|\epsilon \right|<b_1}{inf}{f}_{\epsilon |\mathbf{x}}\left(\epsilon \right|x)\ge b_2$;

(A3). Basis function ${\phi }_1,\cdots ,{\phi }_{L_n}$ satisfies the following:

(1). $L_n\approx n^{1/(2r+1)},r>1/2$;

(2). $\underset{\mathbf{x}\in {[0,1]}^p}{sup}{\parallel \mathsf{\Pi }\parallel }_2=O\left(L_n\right)$ and $H_n={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{\mathsf{\Pi }}_i{\mathsf{\Pi }}_i^{\top }$ are uniformly bounded from 0 to ∞;

(3). There is a vector $\mathit{\gamma }={({\mathit{\beta }}_0^{\top },{\mathit{\beta }}_{00}^{\top })}^{\top }$ such that ${sup}_{\mathbf{x}\in {[0,1]}^p}|m\left(\mathbf{x}\right)-{\mathsf{\Pi }}^{\top }\mathit{\gamma }|=O\left(L_n^{-r}\right)$, where $m\left(\mathbf{x}\right)={\displaystyle \sum _{j=1}^p}{f}_j\left({\mathrm{x}}_j\right)+{\displaystyle \sum _{1<j<k\le p}}{f}_{jk}({\mathrm{x}}_j,{\mathrm{x}}_k)$ and ${\mathrm{x}}_j={({\mathrm{x}}_{1j},\dots ,{\mathrm{x}}_{nj})}^{\top }$;

(A4). $q=O\left(n^{C_1}\right)$ for some $C_1<{\displaystyle \frac{1}{3}}$.

It should be noticed that (A1) and (A2) are typical assumptions in nonparametric quantile regression (see [30]). Condition (3) in (A3) requires that the true component function $f_j\left({\mathrm{x}}_j\right)$ and $f_{jk}({\mathrm{x}}_j,{\mathrm{x}}_k)$ can be uniformly well-approximated by the basis ${\phi }_{jl}\left(x\right),l=1,\dots ,L_n$ and ${\phi }_{js}\left(x\right){\phi }_{kd}\left(x\right),s=1,\dots ,L_n,d=1,\dots ,L_n$. Finally, Condition (A4) is set to control the model size.

Theorem 1. 

Let ${\mathit{\gamma }}_{01}={({\mathit{\beta }}_{01}^{\top },{\mathit{\beta }}_{001}^{\top })}^{\top }$. Assume that Conditions (A1)–(A4) hold. Then, the following are true: (a) 

$\parallel {\widehat{\mathit{\gamma }}}_{01}-{\mathit{\gamma }}_{01}{\parallel }_2=o_p\left(n^{-1/2}{L}_n\right)$;

(b) 

$\parallel \widehat{m}{\left(\mathbf{x}\right)-m\left(\mathbf{x}\right)\parallel }_{L_2}=O\left(n^{-(2r-1)/2(2r+1)}\right),$ where $m\left({\mathbf{x}}_i\right)={\displaystyle \sum _{j=1}^p}{f}_j\left(x_{ij}\right)+{\displaystyle \sum _{1\le j<k\le p}}{f}_{jk}(x_{ij},x_{ik})$.

Theorem 1 summarizes the rate of convergence for the oracle estimator, which is the same as the optimal convergence rate in the additive model [30] due to the merit of the additive structure. The proof of Theorem 1 can be found in Appendix A.1.

3. Penalized Estimation for Additive Quantile Regression with Nonlinear Interaction Structures

3.1. Penalized Estimator

Model (2) contains a large number of main effects $f_j\left(x_{ij}\right)$ and interaction effects $f_{jk}(x_{ij},x_{ik})$, which will significantly increase the computational burden, especially in cases where p (the number of features) is large. It is necessary to use penalty methods to ensure computational efficiency and prevent overfitting. Applying sparsity penalties (such as Lasso or grouped Lasso) to each $f_j$ can help to automatically select important variables and reduce the complexity of the model.

However, when using a large number of basis functions to fit complex relationships, a simple sparse penalty may oversimplify the model, ignore potential patterns in the data, and introduce unstable fluctuations. Therefore, in addition to sparse penalties, it is also necessary to introduce smoothness penalties (such as second-order derivative penalties) to avoid drastic fluctuations in the function and ensure smooth and stable fitting results. This method helps to control model complexity, reduce overfitting, and improve the model’s generalization ability.

Therefore, we minimize the following penalized objective function for $({\mathit{\beta }}_0,{\mathit{\beta }}_{00})$:

(4)${\displaystyle \frac{1}{n}}\sum _{i=1}^n{\rho }_{\tau }\left(Y_i-\sum _{j=1}^p{\mathbf{\Psi }}_{i,j}^{\top }{\mathit{\beta }}_j+\sum _{1\le j<k\le p}{\mathsf{\Phi }}_{i,jk}^{\top }{\mathit{\beta }}_{jk}\right)+P\left(f\right),$

where

$P\left(f\right)={\lambda }_1\left[\sum _{j=1}^p{\left\{\parallel f_j\left({\mathrm{x}}_j\right){\parallel }_n^2+{\lambda }_2{J}_2^1\left(f_j\left({\mathrm{x}}_j\right)\right)\right\}}^{1/2}+\sum _{1\le j<k\le p}{\left\{\parallel f_{jk}({\mathrm{x}}_j,{\mathrm{x}}_k){\parallel }_n^2+{\lambda }_2{J}_2^1\left(f_{jk}({\mathrm{x}}_j,{\mathrm{x}}_k)\right)\right\}}^{1/2}\right],$

is a sparsity smooth penalty function, and the first and second terms in the penalty function penalize the main and interaction effects, respectively. We employ ${\parallel f\parallel }_n^2={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{f}_i^2$ as a sparsity penalty, which encourages sparsity at the function level, and $J_2^1\left(f\right)=\int f^{″2}\left(x\right)dx$ as a roughness penalty, which controls the model’s complexity by ensuring that the estimated function remains smooth while fitting the data, thus preventing high-frequency fluctuations caused by an excessive number of basis functions. The two tuning parameters ${\lambda }_1,{\lambda }_2\ge 0$ control the amount of penalization.

It can be written that $J_2^1\left(f_j\left({\mathrm{x}}_j\right)\right)={\mathit{\beta }}_j^{\top }{\mathsf{\Omega }}_j{\mathit{\beta }}_j$, where ${\mathsf{\Omega }}_j$ is a band diagonal matrix of known coefficients and the $(l_1,l_2)$th element can be expressed as ${\mathsf{\Omega }}_{j,l_1,l_2}=\int {\phi }_{j{l}_1}^{″}\left({\mathrm{x}}_j\right){\phi }_{j{l}_2}^{″}\left({\mathrm{x}}_j\right)dx,l_1,l_2\in \{1,\dots ,L_n\}$—see [31] for details. According to [32], the penalty $J_2^1\left(f_{jk}({\mathrm{x}}_j,{\mathrm{x}}_k)\right)$ can be represented as

${\mathit{\beta }}_{jk}^{\top }({\mathsf{\Lambda }}_{jj}\otimes I_{L_n}+I_{L_n}\otimes {\mathsf{\Lambda }}_{kk}){\mathit{\beta }}_{jk}\triangleq {\mathit{\beta }}_{jk}^{\top }{\mathsf{\Lambda }}_{jk}{\mathit{\beta }}_{jk},$

where ${\mathsf{\Lambda }}_{jj}=\int f_j^{″2}d{x}_j$ and ${\mathsf{\Lambda }}_{kk}=\int f_k^{″2}d{x}_k$, with $f_j^{″}$ and $f_k^{″}$ being the second derivatives of $f({\mathrm{x}}_j,{\mathrm{x}}_k)$ with respect to ${\mathrm{x}}_j$ and ${\mathrm{x}}_k$, respectively. Hence, the penalty function can be rewritten as

$P\left(f\right)={\lambda }_1\left(\sum _{j=1}^p\sqrt{{\mathit{\beta }}_j^{\top }{M}_j{\mathit{\beta }}_j}+\sum _{1\le j<k\le p}\sqrt{{\mathit{\beta }}_{jk}^{\top }{K}_{jk}{\mathit{\beta }}_{jk}}\right),$

where $M_j={\displaystyle \frac{1}{n}}{\mathit{\psi }}_j^{\top }{\mathit{\psi }}_j+{\lambda }_2{\mathsf{\Omega }}_j$ and $K_{jk}={\displaystyle \frac{1}{n}}{\mathit{\varphi }}_{jk}^{\top }{\mathit{\varphi }}_{jk}+{\lambda }_2{\mathsf{\Lambda }}_{jk}$, with ${\mathit{\psi }}_j$ denoting the $n\times L_n$ matrix with the $(i,k)$-th entry given by ${\phi }_{jk}\left(x_{ij}\right)$, and ${\mathit{\varphi }}_{jk}$ representing the $n\times L_n$ matrix with the $(i,l)$-th entry given by ${\phi }_{jl}\left(x_{ij}\right){\phi }_{kl}\left(x_{ik}\right)$. Similar to [33], we decompose $M_j=R_j^{\top }{R}_j$ and $K_{jk}=Q_{jk}^{\top }{Q}_{jk}$ for some quadratic $L_n\times L_n$ matrix $R_j$ and $L_n^2\times L_n^2$ matrix $Q_{jk}$. Then, Model (4) can be represented as

(5)${\displaystyle \frac{1}{n}}\sum _{i=1}^n\rho \left(Y_i-\sum _{j=1}^p{\mathbf{\Psi }}_{i,j}^{\top }{\mathit{\beta }}_j-\sum _{1\le j<k\le p}{\mathsf{\Phi }}_{i,jk}^{\top }{\mathit{\beta }}_{jk}\right)+{\lambda }_1\left(\sum _{j=1}^p\parallel R_j{\mathit{\beta }}_j{\parallel }_2+\sum _{1\le j<k\le p}{\parallel Q_{jk}{\mathit{\beta }}_{jk}\parallel }_2\right).$

However, the above penalty treats the main and interaction effects similarly. That is, an entry of an interaction into the model generally adds more predictors than an entry of a main effect, which usually demands high computational cost and is hard to interpret for more interactions. To deal with this difficulty, we propose to add a set of heredity restrictions to produce sparse interaction models that an interaction only be included in a model if one or both variables are marginally important.

Naturally, $T\subseteq S^{\circ 2}.$ Our target is to estimate S and T from the data consistently or sign consistently, like Th1 in [21]

There are two types of heredity restrictions, which iare called strong and weak hierarchy:

$\mathrm{Strong}\mathrm{hierarchy}:{\widehat{\mathit{\beta }}}_{jk}\ne 0\mathrm{and}{\widehat{\mathit{\beta }}}_j{\widehat{\mathit{\beta }}}_k\ne 0$

$\mathrm{Weak}\mathrm{hierarchy}:{\widehat{\mathit{\beta }}}_{jk}\ne 0\mathrm{and}max\{\parallel {\widehat{\mathit{\beta }}}_j{\parallel }_0,\parallel {\widehat{\mathit{\beta }}}_k{\parallel }_0\}\ne 0.$

To ensure the strong and weak hierarchy conditions for ${\widehat{\mathit{\beta }}}_{jk}$ and ${\widehat{\mathit{\beta }}}_j,{\widehat{\mathit{\beta }}}_k$, we structure the following penalties to enforce the strong hierarchy:

(6)$\begin{array}{ccc}& & Q^P({\mathit{\beta }}_j,{\mathit{\beta }}_{jk})={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\rho }_{\tau }(Y_i-\sum _{j=1}^p{\mathbf{\Psi }}_{i,j}^{\top }{\mathit{\beta }}_j-\sum _{1\le j<k\le p}{\mathsf{\Phi }}_{i,jk}^{\top }{\mathit{\beta }}_{jk})\hfill \\ & & +{\lambda }_1\left({\left(\sum _{j=1}^p\parallel R_j{\mathit{\beta }}_j{\parallel }_2+\sum _{1\le k<j\le p}{I}_{{\beta }_{jk\ne 0}}(\parallel R_j{\mathit{\beta }}_j{\parallel }_2+\parallel R_k{\mathit{\beta }}_k{\parallel }_2)\right)}^{1/2}+\sum _{1\le k<j\le p}{\parallel Q_{jk}{\mathit{\beta }}_{jk}\parallel }_2\right).\hfill \end{array}$

and weak hierarcy:

(7)$\begin{array}{ccc}& & Q^P({\mathit{\beta }}_j,{\mathit{\beta }}_{jk})={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\rho }_{\tau }(Y_i-\sum _{j=1}^p{\mathbf{\Psi }}_{i,j}^{\top }{\mathit{\beta }}_j-\sum _{1\le j<k\le p}{\mathsf{\Phi }}_{i,jk}^{\top }{\mathit{\beta }}_{jk})\hfill \\ & & +{\lambda }_1\left({\left(\sum _{j=1}^p\parallel R_j{\mathit{\beta }}_j{\parallel }_2+\sum _{1\le k<j\le p}{I}_{{\beta }_{jk\ne 0}}min(\parallel R_j{\mathit{\beta }}_j{\parallel }_2,\parallel R_k{\mathit{\beta }}_k{\parallel }_2)\right)}^{1/2}+\sum _{1\le k<j\le p}{\parallel Q_{jk}{\mathit{\beta }}_{jk}\parallel }_2\right).\hfill \end{array}$

These structures ensure hierarchy through the introduction of an indicative penalty. Specifically, the indicator function $I_{{\mathit{\beta }}_{jk\ne 0}}$ activates the penalty only when ${\mathit{\beta }}_{jk}\ne 0$, thus enforcing the strong condition. This ensures that both corresponding main effects must be present if the interaction term is non-zero. Additionally, the term $min(\parallel R_j{\mathit{\beta }}_j{\parallel }_2,\parallel R_k{\mathit{\beta }}_k{\parallel }_2)$ ensures that the penalty is applied to the smaller of the two norms, thereby enforcing the weaker condition. This means that at least one of the corresponding main effects must be present for an interaction term to be included.

3.2. Algorithm

The regularization path algorithm under the marginality principle (RAMP, [21]) is a method for variable selection in high-dimensional quadratic regression models, designed to maintain the hierarchical structure between main effects and interaction effects during the selection process. However, the RAMP algorithm primarily focuses on linear interaction terms and considers interactions of variables with themselves, making it unsuitable for directly handling nonlinear interaction terms.

To address this limitation, we modified the RAMP algorithm, extending its application to additive models that include nonlinear interaction models. Through this enhancement, we can not only handle complex nonlinear interactions but also continue to preserve the hierarchical structure between main effects and interaction effects during variable selection. Specifically, the detailed steps of the modified algorithm are as follows.

Let $\mathcal{P}=\{1,2,\dots ,p\}$ and $\mathcal{Q}=\left\{\right(j,k):1\le j\le k\le p\}$ be the index set for main effects and interaction effects, respectively. For an index set $\mathcal{A}\subset \mathcal{P}$, define ${\mathcal{A}}^{\circ 2}=\mathcal{A}\circ \mathcal{A}=\{(j,k):j\le k;j,k\in \mathcal{A}\}\subset \mathcal{Q}$, and $\mathcal{A}\circ \mathcal{P}=\left\{\right(j,k):j\le k;\mathrm{or}j,k\in \mathcal{A}\}\subset \mathcal{Q}$.

We fix a sequence of values ${\left\{{\lambda }_{1,s}\right\}}_{s=1}^S$ between ${\lambda }_{1,max}$ and ${\lambda }_{1,min}$. Following [21], we set ${\lambda }_{1,max}=n^{-1}max|X^{\top }y|$ and ${\lambda }_{1,min}=\zeta {\lambda }_{1,max}$ with some small $\zeta >0$. At step $s-1$, we denote the current active main effects set as ${\mathcal{P}}_{s-1}$, and the interaction effects are set as ${\mathcal{Q}}_{s-1}$. Define ${\mathcal{H}}_{s-1}$ as the parent set of ${\mathcal{Q}}_{s-1}$, which contains the main effects that have at least one interaction effect in ${\mathcal{Q}}_{s-1}$. Set ${\mathcal{H}}_{s-1}^c=\mathcal{P}-{\mathcal{H}}_{s-1}$. Then, the algorithm is as detailed below:

Step 1. Generate a decreasing sequence ${\lambda }_{1,max}={\lambda }_{1,1}>{\lambda }_{1,2}>\dots ,{\lambda }_{1,S}={\lambda }_{1,min}$, and ${\lambda }_{2,s}={\lambda }_{1,s}^2$ for $s=1,\dots ,S$.

We use the warm start strategy, where the solution for $s-1$ is used as a starting value for s;

Step 2. Given ${\mathcal{P}}_{s-1},{\mathcal{Q}}_{s-1},{\mathcal{H}}_{s-1}$, add the possible interactions among main effects in ${\mathcal{P}}_{s-1}$ to the current model. Then, we minimize the penalty loss function

$\begin{array}{ccc}& & {\displaystyle \frac{1}{n}}\sum _{i=1}^n{\rho }_{\tau }\left(Y_i-\sum _{j\in \mathcal{P}}{\mathbf{\Psi }}_{i,j}^{\top }{\mathit{\beta }}_j-\sum _{(j,k)\in {\mathcal{P}}_{s-1}^{\circ 2}}{\mathsf{\Phi }}_{i,jk}^{\top }{\mathit{\beta }}_{jk}\right)\hfill \\ & & +{\lambda }_{1,s}\left(\sum _j{I}_{j\in {\mathcal{H}}_{s-1}^c}\parallel R_j{\mathit{\beta }}_j{\parallel }_2+\sum _{j,k}{I}_{(j,k)\in {\mathcal{P}}_{s-1}^{\circ 2}}\parallel Q_{jk}{\mathit{\beta }}_{jk}{\parallel }_2\right)\hfill \end{array}$

where the penalty is imposed on the candidate interaction effects and ${\mathcal{H}}_{s-1}^c$, which contains the main effects not enforced by the strong heredity constraint. The penalty encourages smoothness and sparsity in the estimated functional components;

Step 3. Record ${\mathcal{P}}_s,{\mathcal{Q}}_s$, and ${\mathcal{H}}_s$ according to the above solution. Add the corresponding effects from ${\mathcal{Q}}_s$ into ${\mathcal{P}}_s$;

Step 4. Calculate the quantile estimation based on the current model

$({\widehat{\mathit{\beta }}}_j^{\left(s\right)},{\widehat{\mathit{\beta }}}_{jk}^{\left(s\right)})=argmin{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\rho }_{\tau }\left(Y_i-\sum _{j\in {\mathcal{P}}_s}{\mathbf{\Psi }}_{i,j}^{\top }{\mathit{\beta }}_j-\sum _{(j,k)\in {\mathcal{Q}}_s}{\mathsf{\Phi }}_{i,jk}^{\top }{\mathit{\beta }}_{jk}\right)$

Step 5. Repeat steps 2–4 S times and determine the active sets ${\mathcal{M}}_s$ and ${\mathcal{I}}_s$ according to ${\lambda }_s$, which minimizes GIC.

Ref. [34] proposed the following GIC criterion:

$\begin{array}{c}\hfill GIC\left({\lambda }_{1s}\right)=log\left(\sum _{i=1}^n{\rho }_{\tau }(Y_i-\sum _{j\in {\mathcal{P}}_s}{\mathbf{\Psi }}_{i,j}^{\top }{\widehat{\mathit{\beta }}}_j^{\left(s\right)}-\sum _{(j,k)\in {\mathcal{Q}}_s}{\mathsf{\Phi }}_{i,jk}^{\top }{\widehat{\mathit{\beta }}}_{jk}^{\left(s\right)})\right)+C_n{\kappa }_{s}d{f}_{\lambda },\end{array}$

where ${\widehat{\mathit{\beta }}}_j^{\left(s\right)}$ and ${\widehat{\mathit{\beta }}}_{jk}^{\left(s\right)}$ are the penalized estimators obtained by step 4, ${\kappa }_s=|{\mathcal{P}}_s|+|{\mathcal{Q}}_s|$ is the cardinality of the index set of nonzero coefficients in main effects and interaction effects, $d{f}_{\lambda }$ represents the degrees of freedom of the model, and $C_n$ is a sequence of positive constants diverging to infinity as n increases. We take $C_n=log(log\left(n\right))$ in our simulation studies and real data analysis.

Following the same strategy, under the weak heredity condition, we use the set ${\mathcal{P}}_{s-1}\circ \mathcal{P}$ instead of ${\mathcal{P}}_{s-1}^{\circ 2}$ and solve the following optimization problem:

$\begin{array}{ccc}& & {\displaystyle \frac{1}{n}}\sum _{i=1}^n{\rho }_{\tau }\left(Y_i-\sum _{j\in \mathcal{P}}{\mathbf{\Psi }}_{i,j}^{\top }{\mathit{\beta }}_j-\sum _{(j,k)\in {\mathcal{P}}_{s-1}\circ \mathcal{P}}{\mathsf{\Phi }}_{i,jk}^{\top }{\mathit{\beta }}_{jk}\right)\hfill \\ & & +{\lambda }_{1,s}\left(\sum _j{I}_{j\in {\mathcal{H}}_{s-1}^c}\parallel R_j{\mathit{\beta }}_j{\parallel }_2+\sum _{(j,k)}{I}_{(j,k)\in {\mathcal{P}}_{s-1}\circ \mathcal{P}}\parallel Q_{jk}{\mathit{\beta }}_{jk}{\parallel }_2\right),\hfill \end{array}$

with respect to ${\mathit{\beta }}_j$ and ${\mathit{\beta }}_{jk}$. That is, an interaction effect can enter the model for selection if at least one of its parents are selected in a previous step.

3.3. Asymptotic Theory

Let ${\Xi }_n\left(\lambda \right)$ be the set of local minima of $Q^P({\mathit{\beta }}_j,{\mathit{\beta }}_{jk})$. The following result shows that, with probability approaching one, the oracle estimator belongs to the set ${\Xi }_n\left(\lambda \right)$.

Theorem 2. 

Assume that conditions (A1)–(A4) are satisfied and consider the penalized function under strong heredity. Let $\widehat{\mathit{\gamma }}={({\widehat{\mathit{\beta }}}_0^{\top },{\widehat{\mathit{\beta }}}_{00}^{\top })}^{\top }$ be the oracle estimator. If ${\lambda }_1=o\left(1\right)$ and $n^{-1}{L}_n^3{\lambda }^{-1}\to 0$ as $n\to \infty$, then

$P(\widehat{\mathit{\gamma }}\in {\Xi }_n\left(\lambda \right))\to 1,\mathrm{as}n\to \infty .$

Remark 1. 

The penalty term ${\lambda }_1$ decays to zero as $n\to \infty$, allowing the oracle estimator to dominate. The rate condition $n^{-1}{L}_n^3{\lambda }^{-1}\to 0$ guarantees that the approximation error from basis expansion $L_n$ is controlled by the penalty strength.

The proof of Theorem 2 can be found in Appendix A.2.

Define ${\mathcal{R}}_n\left(m\right)={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{\rho }_{\tau }(Y_i-m\left({\mathbf{x}}_i\right))$ and $\mathcal{R}\left(m\right)=\mathbb{E}\left({\rho }_{\tau }(Y-m\left(\mathbf{x}\right))\right)$ as the empirical risk and predictive risk of quantile loss, respectively. And

$\begin{array}{ccc}\hfill {\mathcal{M}}_n& =& \{m:m\left(x\right)=\sum _{j=1}^p{\mathit{\psi }}_j^{\top }\left(x_j\right){\mathit{\beta }}_j+\sum _{1<j<k\le p}{\mathit{\varphi }}_{jk}^{\top }(x_j,x_k){\mathit{\beta }}_{jk}:\hfill \\ & & E\left({\mathit{\psi }}_j\right)=0,E\left({\mathit{\varphi }}_{jk}\right)=0,E\left({\mathit{\psi }}_j^2\right)=1,E\left({\mathit{\varphi }}_{jk}^2\right)=1\}\hfill \end{array}$

is a functional class. Let $\tilde{m}=\underset{m\in {\mathcal{M}}_n}{argmin}\mathcal{R}\left(m\right)$ denote the predictive oracle, i.e., the minimizer of the predictive risk over ${\mathcal{M}}_n$, and $\widehat{m}$ represent the minimization of Equation (6) over ${\mathcal{M}}_n$. Following [35], we say that an estimator $\widehat{m}$ is persistent (risk consistent) relative to a class of functions ${\mathcal{M}}_n$ if $\mathcal{R}\left(\widehat{m}\right)-\mathcal{R}\left(\tilde{m}\right)\stackrel{P}{⟶}0$. Then, we have the following result.

Theorem 3. 

Under conditions (A1)–(A4), if $p=o\left(n\right)$, then for some constants $t>0$,

$P\left(\mathcal{R}\left(\widehat{m}\right)-\mathcal{R}\left(\tilde{m}\right)\le (2\sqrt{2t^2}+10t^2/3)n^{-(2r-1)/2(2r+1)}\right)\ge 1-exp(-n{t}^2).$

The above theorem establishes the convergence rate in terms of the excess risk of the estimator, which shows the prediction accuracy and consistency of the proposed estimator. The proof of Theorem 3 can be found in Appendix A.3.

4. Simulation Studies

We conduct extensive Monte Carlo simulations to assess the performance of our proposed method in comparison with two established approaches, hirNet [19] and RAMP [21], across five specific quantile levels: $\tau \in 0.1,0.25,0.5,0.75,0.90$. In each simulation scenario, we generate 500 independent training datasets with two sample sizes, $n=300$ and $n=100$, and $p=15$ main effects, resulting in a total of $15+15\times 14/2=121$ potential terms (15 main effects and 105 pairwise interactions). All methods are implemented using pseudo-spline basis functions to ensure consistent comparison. For regularization parameters, we adopt the relationship ${\lambda }_2={\lambda }_1^2$ based on [13], which provides superior empirical performance.

The covariates $x_{ij}$ are generated independently from a uniform distribution $U(0,1)$. We specify five nonlinear main effect functions that are standardized before model fitting:

$f_1\left(x\right)=10x,f_2\left(x\right)={\displaystyle \frac{1}{1+x}}+20x,$

$f_3\left(x\right)=20\sin \left(x\right),f_4\left(x\right)=20\exp \left(x\right),f_5\left(x\right)=10x^2.$

Hence, each $f_j$ is standardized and the interaction functions are generated by multiplying together the standardized main effects,

$f_{jk}(x_j,x_k)=f_j\left(x_j\right)\times f_k\left(x_k\right),1\le j<k\le p.$

The response variables are generated through quantile-specific models that incorporate either strong or weak heredity constraints. Under strong heredity, interactions are only included when all parent main effects are active, while weak heredity requires just one parent main effect.

Example 1. 

Strong heredity. Building on the quantile-specific sparsity framework established by [36], we incorporate strong heredity constraints into our modeling approach. This integration ensures two key aspects: (1) the set of active predictors varies depending on the quantile, and (2) all selected variables preserve strict hierarchical relationships. The resulting nonlinear data-generating process is as follows:

$\tau =0.1$ (Sparse):

$Y_i=f_1\left(x_{i1}\right)+f_2\left(x_{i2}\right)+f_3\left(x_{i3}\right)+{\epsilon }_i.$

$\tau =0.25,0.50,0.75$ (Dense):

$Y_i=f_1\left(x_{i1}\right)+f_2\left(x_{i2}\right)+f_3\left(x_{i3}\right)+f_4\left(x_{i4}\right)+f_5\left(x_{i5}\right)+f_{12}(x_{i1},x_{i2})+f_{13}(x_{i1},x_{i3})+{\epsilon }_i.$

$\tau =0.9$ (Sparse):

$Y_i=f_1\left(x_{i1}\right)+f_3\left(x_{i3}\right)+f_5\left(x_{i5}\right)+{\epsilon }_i.$

Example 2. 

Weak heredity. In contrast to strong heredity, we also incorporate weak heredity constraints into our modeling approach by extending the quantile-specific sparsity framework. The resulting nonlinear data-generating process is as follows:

$\tau \in 0.1,0.9$ (Sparse):

$Y_i=f_1\left(x_{i1}\right)+f_2\left(x_{i2}\right)+f_3\left(x_{i3}\right)+{\epsilon }_i.$

$\tau =0.25,0.50,0.75$ (Dense):

$Y_i=f_1\left(x_{i1}\right)+f_2\left(x_{i2}\right)+f_3\left(x_{i3}\right)+f_{14}(x_{i1},x_{i4})+f_{15}(x_{i1},x_{i5})+{\epsilon }_i.$

To assess robustness across different error conditions, we consider four error distributions:

1.. Gaussian: The error terms ${\epsilon }_i$ follow a normal distribution with mean 0 and variance 1, i.e., ${\epsilon }_i\sim N(0,1)$;

2.. Heavy-tailed: The error terms ${\epsilon }_i$ follow a Student’s t-distribution with 3 degrees of freedom, i.e., ${\epsilon }_i\sim t\left(3\right)$;

3.. Skewed: The error terms ${\epsilon }_i$ follow a chi-squared distribution with 2 degrees of freedom, i.e., ${\epsilon }_i\sim \chi \left(2\right)$;

4.. Heteroscedasticity [37]: The error terms ${\epsilon }_i$ are heteroscedastic and modeled as follows:

$\tau =0.1:{\epsilon }_i=(x_{i1}+x_{i2})e_i$;

$\tau \in 0.25,0.50,0.75:{\epsilon }_i=(x_{i1}+x_{i5})e_i$;

$\tau =0.9\left(strongheredity\right):{\epsilon }_i=(x_{i1}+x_{i3})e_i$;

$\tau =0.9\left(weakheredity\right):{\epsilon }_i=(x_{i1}+x_{i4})e_i$.

Recall that $\mathcal{M}=\{j:f_j\ne 0\}$ and $\mathcal{I}=\{(j,k):f_{jk}\ne 0\}$ are the active sets of main and interaction effects. For each example, we run $R=500$ Monte Carlo simulations and denote the estimated subsets as ${\widehat{\mathcal{M}}}^{\left(r\right)}$ and ${\widehat{\mathcal{I}}}^{\left(r\right)}$, $r=1,\dots ,R$. We evaluate the performance on variable selection based on the following criteria:

(1). True positive rate and false positive rate of main effects (mTPR, mFDR);

(2). True positive rate and false positive rate of interaction effects (iTPR, iFDR);

(3). Main effects coverage percentage (${\mathcal{P}}_m$): $R^{-1}{\sum }_{r=1}^{R}I(\mathcal{M}\in {\widehat{\mathcal{M}}}^{\left(r\right)})$;

(4). Interaction effects coverage percentage (${\mathcal{P}}_i$): $R^{-1}{\sum }_{r=1}^{R}I(\mathcal{I}\in {\widehat{\mathcal{I}}}^{\left(r\right)})$;

(5). Coverage rate of a single variable in the main effects and interaction effects;

(6). Model size (size): $R^{-1}{\sum }_{r=1}^R\left(\right|{\widehat{\mathcal{M}}}^{\left(r\right)}|+|{\widehat{\mathcal{I}}}^{\left(r\right)}\left|\right)$;

(7). Root mean squared error (RMSE): $R^{-1}{\sum }_{r=1}^R\sqrt{{\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{\left({\rho }_{\tau }(Y_i-{\widehat{Q}}_{\tau }^{\left(r\right)}\left(Y_i\right|{\mathbf{x}}_i))\right)}^2}$, where ${\widehat{Q}}_{\tau }^{\left(r\right)}\left(Y_i\right|{\mathbf{x}}_i)$ is obtained by the estimated coefficients ${\widehat{\mathit{\beta }}}_j^{\left(r\right)}\in {\widehat{\mathcal{M}}}^{\left(r\right)}$ and ${\widehat{\mathit{\beta }}}_{jk}^{\left(r\right)}\in {\widehat{\mathcal{I}}}^{\left(r\right)}$; R is the number of iterations.

The results under the strong heredity for $n=100$ and $n=300$ are shown in Table 1 and Table 2, respectively. In the context of small sample data, the performance of the proposed method is a little worse than the other two methods in terms of true positive rates and coverage percentages of a single variable in the main and interaction effects. This is because the high complexity and degrees of freedom in nonparametric additive interaction models make them more prone to overfitting, leading to an inaccurate reflection of the underlying data structure and, consequently, a lower true positive rate (TPR) in variable selection. In contrast, methods like hirNet and RAMP, by adopting regularization and relatively simple linear structures, can more effectively prevent overfitting. Although this may come at the cost of some loss in model interpretability.

In large sample scenarios, our proposed method not only demonstrates higher true positive rates and broader coverage of main effects and interaction effects but also significantly reduces false positive rates and generates a more parsimonious model. These results indicate that the method can effectively handle nonlinear additive interaction models, accurately capturing complex relationships in the data while maintaining a low false discovery rate and higher model simplicity.

The results from weak heredity are shown in Table 3 and Table 4. In small sample cases, the hirNet method is effective at capturing both main and interaction effects but suffers from a high rate of false selections, leading to increased model complexity. On the other hand, the RAMP method selects fewer variables, reducing false selections but also under-selecting important variables, which negatively affects predictive accuracy. In contrast, the proposed method in this paper offers a better balance, with higher precision and stability in selecting both main and interaction effects, while significantly reducing false selections. This method results in a moderate model size and the lowest prediction error, demonstrating an optimal trade-off between model complexity and predictive performance.

In large sample cases, the proposed method excels in identifying both main and interaction effects, with the main effect true positive rate (mTPR) and interaction effect true positive rate (iTPR) approaching 1.000, indicating a nearly perfect identification of true positives. This method also maintains low mFDR and iFDR, minimizing false positives. In contrast, although the hirNet methods perform better in terms of mTPR and iTPR, they tend to include more noisy variables in the model, resulting in higher RMSE values and lower predictive accuracy.

From the perspective of quantile-specific sparsity, the proposed method demonstrates varying levels of sparsity across different quantiles (e.g., $\tau \in 0.1,0.25,0.5,0.75,0.90$). At $\tau =0.1$ and $\tau =0.9$ quantiles, fewer variables are selected, indicating higher sparsity, while, at $\tau =0.25,0.50,0.9$, more variables and interaction effects are included, showing that the model adapts its selection based on the distribution of the response variable. Simulation results further show that, when the $\tau =0.9$, DGP consists only of main effects and our method achieves near-perfect variable selection: the TPR for main effects approaches 1 (all true effects retained) and the FDR is close to 0 (almost no false positives), leading to model sizes nearly identical to the true DGP and minimal prediction errors, as reflected in the low RMSE. These results are consistent with [38]’s work, which provides theoretical support for this outcome. Moreover, the RMSE patterns suggest that the method maintains strong predictive performance even when more variables are selected, thereby confirming the effectiveness of the quantile-specific sparsity approach.

5. Applications

Parkinson’s disease (PD) is a common degenerative disease of the nervous system that leads to shaking, stiffness, and difficulty with walking, balance, and coordination. PD symptom monitoring with the Unified Parkinson’s Disease Rating Scale (UPDRS) is very important but costly and logistically inconvenient due to the requirement of patients’ presence in clinics and time-consuming physical examinations by trained medical staff. A new treatment technology based on the rapid and remote replication of UPDRS assessments with more than a dozen biomedical voice measurement indicators has been proposed recently. The dataset consists of 5875 UPDRS from patients with early-stage PD and 16 corresponding biomedical voice measurement indicators. Reference [39] analyzed a PD dataset and mapped the clinically relevant properties from speech signals of PD patients to UPDRS using linear regression and nonlinear regression models, with the aim to verify the feasibility of frequent, remote, and accurate UPDRS tracking and effectiveness in telemonitoring frameworks that enable large-scale clinical trials into novel PD treatments. However, existing analyses of the UPDRS scores from the 16 voice measures with a model including only the main effects may fail to reflect the relationship between the UPDRS and 16 voice measures. For example, fundamental frequency and amplitude can work together on voice and thus may have an interactive effect on UPDRS. In this situation, it is necessary to consider the model with pairwise interactions. Furthermore, the effect on UPDRS from these indicators may be complicated and needs to be investigated further, including their nonlinear dependency, even nonparametric relation, so that both main effect and interactive effect may exist. By noting the typical skewness of the distribution of the UPDRS and the asymmetry of its $\tau$-quantile-level and ($1-\tau$)-quantile-level score $(0<\tau <1$) in Figure 1 (left), the distribution of the UPDRS tends to asymmetric and heavy-tailed, so that one should examine how any quantile, including the median (rather than the mean) of the conditional distribution, is affected by changes in the 16 voice measures and then make a reliable decision about the new treatment technology. Therefore, the proposed additive quantile regression with nonlinear interaction structure is applied for the analysis of the data.

The 16 biomedical voice measures in the PD data include several measurements of variation in fundamental frequency (Jitter, Jitter (Abs), Jitter:RAP, Jitter:PPQ5, and Jitter:DDP), several measurements of variation in amplitude (Shimmer, Shimmer (dB), Shimmer:APQ3, Shimmer:APQ5, Shimmer:APQ11, and Shimmer:DDA), two measurements of the ratio of noise to tonal components in the voice (NHR and HNR), a nonlinear dynamical complexity measurement (RPDE), signal fractal scaling exponent (DFA), and a nonlinear measure of fundamental frequency variation (PPE). Our interest is to identify the biomedical voice measurements and their interactions that may effectively affect the UPDRS (Unified Parkinson Disease Rating Scale) in Parkinson’s disease patients. Figure 1 (right) shows the correlations among 16 biomedical voice measurements. As we can see from Figure 1 (right), except for variable DFA, there are strong correlations among the covariates, which makes variable selection more challenging. Therefore, the proposed additive quantile regression model with nonlinear interaction structures may provide a complete picture on how those factors and their interactions affect the distribution of UPDRS.

We first randomly generated 100 partitions of the data into training and test sets. For each training set, we select 5475 observations as training set and the remaining 400 observations as the test set. Normalization is conducted for the dataset. Also, quantile regression is used to fit the model selected by the proposed method. For comparison, we also use least square to fit the model selected by hirNet and RAMP on the training set. Finally, we evaluate the performance of these models using the fixed active sets on the corresponding test sets.

Table 5 summarizes the covariates selected by different methods. The proposed quantile-specific method reveals distinct covariate patterns across different UPDRS severity levels. At lower quantiles ($\tau =0.1$ and $\tau =0.25$), the model highlights DFA and its interactions (e.g., HNRDFA, DFAPPE), suggesting that these acoustic features may be particularly relevant in the early stages of the disease. The middle quantile ($\tau =0.5$) incorporates additional speech markers (APQ11, RPDE), which align with conventional methods, while higher quantiles ($\tau =0.75$ and $\tau =0.9$) progressively simplify to core features (HNR, dB), indicating reduced covariate dependence in advanced stages. This dynamic selection demonstrates how quantile-specific modeling can capture stage-dependent biological mechanisms while maintaining model parsimony, selecting just 3–5 main effects and 1–2 interactions per quantile, compared to hitNet’s fixed 12-variable approach. The shifting importance of HNR and DFA interactions across quantiles particularly underscores their complex, nonlinear relationships with symptom progression.

Figure 2 and Figure 3 show the estimation results of main effects and interaction effects at $\tau =0.50$. The solid lines represent the estimated effects, while the dashed lines indicate the $95\%$ confidence intervals obtained through a simulation-based approach. Specifically, for the chosen model, we conducted 100 repeated simulations on the dataset to derive point-wise confidence intervals for each covariate’s effect.

From Figure 2, it can be observed that, as APQ11 increases, the UPDRS score initially rises and then declines, peaking at a specific value. This inverted U-shaped relationship suggests that moderate amplitude variations are linked to more severe symptoms. For HNR, the UPDRS score decreases as the noise ratio increases, indicating that higher noise levels correlate with disease worsening. The confidence intervals around these trends confirm their statistical significance across the simulated datasets. RPDE shows a U-shaped curve, meaning that intermediate complexity levels are associated with more severe symptoms, while extreme values (high or low) indicate milder symptoms. The confidence intervals further support this nonlinear relationship, demonstrating its robustness to the data. DFA has a negative correlation with the UPDRS score: as fundamental frequency variability increases, the UPDRS score decreases, suggesting that greater variability is linked to symptom relief. The narrow confidence intervals around this trend highlight its consistency across the simulations.

Figure 3 illustrates the interaction effect between HNR and RPDE on the UPDRS score. This shows that, as HNR increases, the impact on the UPDRS score becomes more negative, especially when RPDE is at lower values, indicating a worsening in symptoms with higher noise to harmonic ratio. Conversely, for higher values of RPDE, the interaction effect stabilizes, suggesting less influence on the UPDRS score. The plot reveals a valley-like structure around the center where both HNR and RPDE are near zero, signifying minimal interaction effect. Notably, the most significant interactions occur at extreme values of either HNR or RPDE, highlighting the complex, nonlinear relationship between these variables and Parkinson’s disease symptom severity. The estimation results of the main effects and interaction effects at $\tau =0.1,0.25,0.75,0.9$ are presented in Appendix B.

From the results in Appendix B, it can be seen that, although the main effects selected at $\tau =0.25$ and $\tau =0.75$ are the same, their impacts on UPDRS scores differ due to nonlinear relationships or complex interactions between variables. This highlights the advantage of quantile regression, which not only captures the average trend of the dependent variable but also reveals how these effects vary across different quantiles. This method provides a more comprehensive understanding of the complex relationships between variables and helps to analyze the specific effects of variables on UPDRS under different conditions.

From the perspective of quantile-specific sparsity, the proposed method demonstrates varying levels of covariate selection across different quantiles. At lower quantiles ($\tau =0.1$ and $\tau =0.25$), the model tends to select fewer covariates, focusing primarily on HNR, DFA, and PPE, which suggests a higher degree of sparsity. As the quantile increases to $\tau =0.50$, the model selects more covariates, including APQ11 and RPDE, indicating a reduction in sparsity. At higher quantiles ($\tau =0.75$ and $\tau =0.9$), the model again shows increased sparsity by selecting only a few key covariates such as HNR, DFA, and PPE. This pattern reflects the adaptability of the proposed method in capturing the specific characteristics of the data distribution at different quantiles, thereby achieving optimal sparsity tailored to each quantile level.

Table 6 compares the average RMSE and model sizes across 100 datasets for three methods: the proposed quantile-based approach (evaluated at $\tau =0.1,0.25,0.5,0.75,0.9$), RAMP, and hirNet. The results show that our method outperforms the others in predictive accuracy, particularly at $\tau =0.75$, where it achieves the lowest RMSE of 0.858, outperforming RAMP (0.878) and hirNet (0.887). Additionally, the proposed method uses much smaller models, with only 2.3–5.0 variables, compared to hirNet’s larger model with 18.39 variables. This reduction in model size highlights the method’s efficiency without sacrificing predictive power. The method also performs consistently well across different quantiles, with a particularly strong result at $\tau =0.25$ (RMSE = 0.873), though the RMSE at $\tau =0.5$ is slightly higher (1.00), showing some variation in performance across quantiles. Furthermore, models that include interaction terms perform better than those with only main effects, emphasizing the importance of interaction effects in the model. In conclusion, the proposed method strikes an optimal balance between prediction accuracy, model simplicity, and computational efficiency, making it ideal for applications that require both interpretability and strong predictive performance.

6. Conclusions

We explore the variable selection for additive quantile regression with interactions. We fit the model using the sparsity smooth penalty function and add the regularization algorithm under the marginality principle to the additive quantile regression with interactions, which can select main and interaction effects simultaneously while keeping either the strong or weak heredity constraints. We demonstrate theoretical properties of the proposed method for additive quantile regression with interactions. Simulation studies demonstrate good performance of the proposed model and method.

Also, we applied the proposed method to the Parkinson’s disease (PD) data; our method successfully validates the important finding in the literature that frequent, remote, and accurate UPDRS tracking as a novel PD treatment could be effective in telemonitoring frameworks for PD symptom monitoring and large-scale clinical trials.

Author Contributions

Y.B.: Conceptualization, Formal analysis, Methodology, Software, Writing—original draft, Data curation; J.J.: Formal analysis, Methodology, Writing—review and editing; M.T.: Supervision, Funding acquisition. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The dataset used in this study is publicly available at https://www.worldbank.org/en/home (accessed on 1 June 2023).

Conflicts of Interest

The authors declare no conflictw of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 (Left) Histogram and density curve of UPDRS; (right) Correlation between voice measurements.

Figure 2 Estimated main effects for the APQ11, HNR, RPDE, and DFA variables at $\tau =0.5$.

Figure 3 Estimated interaction term for the HNR and RPDE variables at $\tau =0.5$.

Table 1

Selection and estimation results for strong heredity with $n=100$.

Method	Main								Inter
	${\mathcal{P}}_{\mathbf{m}}$	${\mathrm{x}}_1$	${\mathrm{x}}_2$	${\mathrm{x}}_3$	${\mathrm{x}}_4$	${\mathrm{x}}_5$	mTPR	mFDR	${\mathcal{P}}_{\mathbf{i}}$	${\mathrm{x}}_1{\mathrm{x}}_2$	${\mathrm{x}}_1{\mathrm{x}}_3$	iTPR	iFDR	Size	RMSE
$N(0,1)$
$hirNet$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.045	1.000	1.000	1.000	1.000	0.710	12.140	6.056
$RAMP$	0.920	1.000	1.000	1.000	1.000	0.920	0.984	0.046	1.000	1.000	1.000	1.000	0.382	8.400	7.200
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	-	-	1.000	0.038	-	-	-	-	-	3.420	1.030
$proposed(\tau =0.25)$	0.400	1.000	1.000	1.000	0.920	0.400	0.864	0.068	0.340	0.920	0.920	0.350	0.426	7.580	4.251
$proposed(\tau =0.50)$	0.700	1.000	1.000	1.000	0.980	0.700	0.936	0.052	0.940	0.980	0.940	0.960	0.200	7.360	2.919
$proposed(\tau =0.75)$	0.640	1.000	1.000	1.000	1.000	0.640	0.928	0.021	0.980	0.980	1.000	0.990	0.232	7.320	3.216
$proposed(\tau =0.9)$	1.000	1.000	-	1.000	-	1.000	1.000	0.044	-	-	-	-	-	3.920	1.042
$t\left(3\right)$
$hirNet$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.110	1.000	1.000	1.000	1.000	0.727	12.960	7.392
$RAMP$	0.900	1.000	1.000	1.000	1.000	0.900	0.980	0.050	1.000	1.000	1.000	1.000	0.363	8.300	8.978
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	-	-	1.000	0.056	-	-	-	-	-	3.480	2.994
$proposed(\tau =0.25)$	0.380	1.000	1.000	1.000	0.940	0.380	0.864	0.056	0.520	0.880	0.880	0.550	0.388	7.440	7.718
$proposed(\tau =0.50)$	0.620	1.000	1.000	1.000	1.000	0.620	0.924	0.057	0.940	0.960	0.960	0.960	0.213	7.360	6.033
$proposed(\tau =0.75)$	0.660	1.000	1.000	1.000	1.000	0.660	0.932	0.025	1.000	1.000	1.000	1.000	0.253	7.460	6.146
$proposed(\tau =0.9)$	0.960	0.960	-	1.000	-	1.000	0.986	0.019	-	-	-	-	-	3.480	2.827
${\chi }^2\left(2\right)$
$hirNet$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.128	1.000	1.000	1.000	1.000	0.735	13.300	8.663
$RAMP$	0.860	1.000	1.000	1.000	1.000	0.860	0.972	0.061	1.000	1.000	1.000	1.000	0.425	8.680	10.217
$proposed(\tau =0.1)$	0.980	0.980	1.000	1.000	-	-	0.993	0.032	-	-	-	-	-	3.160	3.931
$proposed(\tau =0.25)$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	0.1000	0.820	0.820	0.100	0.411	7.180	7.638
$proposed(\tau =0.50)$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	0.940	1.000	0.940	0.000	0.210	7.260	8.583
$proposed(\tau =0.75)$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	0.180	0.307	7.640	8.942
$proposed(\tau =0.1)$	0.860	0.860	-	1.000	-	1.000	0.953	0.046	-	-	-	-	-	3.900	6.049
$\sigma \left({\mathbf{x}}_{\mathbf{i}}\right)e_i$
$hirNet$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.031	1.000	1.000	1.000	1.000	0.717	12.240	6.199
$RAMP$	0.900	1.000	1.000	1.000	1.000	0.900	0.980	0.057	1.000	1.000	1.000	1.000	0.411	8.600	7.441
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	-	-	1.000	0.026	-	-	-	-	-	3.380	1.306
$proposed(\tau =0.25)$	0.404	1.000	1.000	1.000	0.957	0.404	0.872	0.072	0.914	0.936	0.978	0.957	0.403	7.914	2.893
$proposed(\tau =0.50)$	0.640	1.000	1.000	1.000	0.980	0.6400	0.924	0.057	0.940	0.980	0.940	0.960	0.219	7.380	2.355
$proposed(\tau =0.75)$	0.660	1.000	1.000	1.000	1.000	0.660	0.932	0.029	0.980	0.980	1.000	0.990	0.232	7.380	2.706
$proposed(\tau =0.9)$	1.000	1.000	-	1.000	-	1.000	1.000	0.063	-	-	-	-	-	4.080	1.344

Table 2

Selection and estimation results for strong heredity with $n=300$.

Method	Main								Inter
	${\mathcal{P}}_{\mathbf{m}}$	${\mathrm{x}}_1$	${\mathrm{x}}_2$	${\mathrm{x}}_3$	${\mathrm{x}}_4$	${\mathrm{x}}_5$	mTPR	mFDR	${\mathcal{P}}_{\mathbf{i}}$	${\mathrm{x}}_1{\mathrm{x}}_2$	${\mathrm{x}}_1{\mathrm{x}}_3$	iTPR	iFDR	Size	RMSE
$N(0,1)$
$hirNet$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.695	11.560	6.158
$RAMP$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.629	10.840	2.058
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	-	-	1.000	0.000	-	-	-	-	-	3.000	0.978
$proposed(\tau =0.25)$	1.000	1.000	1.000	1.000	1.000	1.00	1.000	0.000	1.000	1.000	1.000	1.000	0.091	7.200	2.520
$proposed(\tau =0.50)$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.000	7.000	2.472
$proposed(\tau =0.75)$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.038	7.080	1.120
$proposed(\tau =0.9)$	1.000	1.000	-	1.000	-	1.000	1.000	0.000	-	-	-	-	-	3.040	0.958
$t\left(3\right)$
$hirNet$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.696	11.580	8.741
$RAMP$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.635	11.320	4.641
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	-	-	1.000	0.000	-	-	-	-	-	3.020	3.765
$proposed(\tau =0.25)$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.074	7.160	4.684
$proposed(\tau =0.50)$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.000	7.000	4.664
$proposed(\tau =0.75)$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.029	7.060	4.786
$proposed(\tau =0.1)$	1.000	1.000	-	1.000	-	1.000	1.000	0.000	-	-	-	-	-	3.060	3.658
${\chi }^2\left(2\right)$
$hirNet$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.698	11.660	9.979
$RAMP$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.375	1.000	1.000	1.000	1.000	0.646	11.660	5.843
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	-	-	1.000	0.000	-	-	-	-	-	3.000	4.913
$proposed(\tau =0.25)$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.056	7.120	5.702
$proposed(\tau =0.50)$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.029	7.060	5.777
$proposed(\tau =0.75)$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.056	7.120	5.787
$proposed(\tau =0.9)$	0.940	0.940	-	1.000	-	1.000	0.980	0.000	-	-	-	-	-	3.100	6.053
$\sigma \left({\mathbf{x}}_{\mathbf{i}}\right)e_i$
$hirNet$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.696	11.660	6.282
$RAMP$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.100	1.000	1.000	1.000	1.000	0.625	10.900	2.181
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	-	-	1.000	0.000	-	-	-	-	-	3.040	1.346
$proposed(\tau =0.25)$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.090	7.200	2.777
$proposed(\tau =0.50)$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.010	7.020	2.693
$proposed(\tau =0.75)$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.000	1.000	1.000	1.000	1.000	0.010	7.020	2.782
$proposed(\tau =0.9)$	1.000	1.000	-	1.000	-	1.000	1.000	0.000	-	-	-	-	-	3.080	1.357

Table 3

Selection and estimation results for weak heredity with $n=100$.

Method	Main								Inter
	${\mathcal{P}}_{\mathbf{m}}$	${\mathrm{x}}_1$	${\mathrm{x}}_2$	${\mathrm{x}}_3$	mTPR	mFDR	${\mathcal{P}}_{\mathbf{i}}$	${\mathrm{x}}_1{\mathrm{x}}_4$	${\mathrm{x}}_1{\mathrm{x}}_5$	iTPR	iFDR	Size	RMSE
$N(0,1)$
$hirNet$	0.980	1.000	1.000	1.000	1.000	0.446	1.000	1.000	1.000	1.000	0.781	14.58	31.954
$RAMP$	0.060	1.000	0.260	0.240	0.500	0.050	1.000	1.000	1.000	1.000	0.619	7.140	67.840
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	1.000	0.047	-	-	-	-	-	3.500	1.067
$proposed(\tau =0.25)$	0.340	1.000	0.580	0.480	0.686	0.507	0.720	1.000	0.720	0.860	0.505	7.860	48.135
$proposed(\tau =0.50)$	0.320	1.000	0.520	0.500	0.673	0.504	0.740	1.000	0.740	0.870	0.462	7.480	42.688
$proposed(\tau =0.75)$	0.340	1.000	0.520	0.520	0.680	0.507	0.780	1.000	0.780	0.890	0.476	7.780	41.641
$proposed(\tau =0.9)$	1.000	1.000	1.000	1.000	1.000	0.062	-	-	-	-	-	4.000	1.058
$t\left(3\right)$
$hirNet$	1.000	1.000	1.000	1.000	1.000	0.450	1.000	1.000	1.000	1.000	0.788	14.920	32.742
$RAMP$	0.060	1.000	0.280	0.240	0.506	0.037	1.000	1.000	1.000	1.000	0.619	7.140	69.608
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	1.000	0.006	-	-	-	-	-	7.060	1.710
$proposed(\tau =0.25)$	0.340	1.000	0.500	0.520	0.673	0.502	0.760	1.000	0.760	0.880	0.443	7.380	46.631
$proposed(\tau =0.50)$	0.360	1.000	0.540	0.540	0.693	0.497	0.720	1.000	0.720	0.860	0.459	7.440	32.514
$proposed(\tau =0.75)$	0.360	1.000	0.560	0.540	0.700	0.492	0.720	1.000	0.720	0.860	0.679	7.680	52.425
$proposed(\tau =0.9)$	1.000	1.000	1.000	1.000	1.000	0.010	-	-	-	-	-	7.000	1.765
${\chi }^2\left(2\right)$
$hirNet$	0.980	1.000	1.000	1.000	1.000	0.468	1.000	1.000	1.000	1.000	0.794	15.360	33.515
$RAMP$	0.040	1.000	0.240	0.220	0.486	0.087	1.000	1.000	1.000	1.000	0.628	7.300	69.560
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	1.000	0.051	-	-	-	-	-	6.940	2.821
$proposed(\tau =0.25)$	0.280	1.000	0.500	0.400	0.633	0.515	0.620	0.980	0.620	0.810	0.683	7.280	41.928
$proposed(\tau =0.50)$	0.340	1.000	0.540	0.460	0.666	0.504	0.640	0.980	0.640	0.810	0.480	7.380	42.535
$proposed(\tau =0.75)$	0.360	1.000	0.580	0.500	0.693	0.500	0.680	0.980	0.680	0.830	0.522	7.920	58.314
$proposed(\tau =0.9)$	1.000	1.000	1.000	1.000	1.000	0.011	-	-	-	-	-	6.166	2.929
$\sigma \left({\mathbf{x}}_{\mathbf{i}}\right)e_i$
$hirNet$	1.000	1.000	1.000	1.000	1.000	0.444	1.000	1.000	1.000	1.000	0.784	14.700	32.115
$RAMP$	0.000	0.300	0.380	0.120	0.266	0.259	1.000	1.000	0.300	0.650	0.821	8.360	13.755
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	1.000	0.124	-	-	-	-		7.775	0.681
$proposed(\tau =0.25)$	0.333	1.000	0.562	0.500	0.687	0.505	0.750	1.000	0.750	0.875	0.478	7.687	24.665
$proposed(\tau =0.50)$	0.354	1.000	0.604	0.479	0.694	0.502	0.750	1.000	0.750	0.875	0.481	7.750	26.570
$proposed(\tau =0.75)$	0.354	1.000	0.625	0.520	0.715	0.500	0.791	1.000	0.791	0.895	0.500	8.083	29.956
$proposed(\tau =0.9)$	1.000	1.000	1.000	1.000	1.000	0.013	-	-	-	-	-	6.560	0.760

Table 4

Selection and estimation results for weak heredity with $n=300$.

Method	Main								Inter
	${\mathcal{P}}_{\mathbf{m}}$	${\mathrm{x}}_1$	${\mathrm{x}}_2$	${\mathrm{x}}_3$	mTPR	mFDR	${\mathcal{P}}_{\mathbf{i}}$	${\mathrm{x}}_1{\mathrm{x}}_4$	${\mathrm{x}}_1{\mathrm{x}}_5$	iTPR	iFDR	Size	RMSE
$N(0,1)$
$hirNet$	1.000	1.000	1.000	1.000	1.000	0.401	0.160	0.900	0.160	0.530	0.047	19.600	32.574
$RAMP$	1.000	1.000	0.700	0.720	1.000	0.000	1.000	1.000	1.000	1.000	0.695	9.360	68.139
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	1.000	0.000	-	-	-	-	-	3.620	0.908
$proposed(\tau =0.25)$	1.000	1.000	1.000	1.000	1.000	0.400	1.000	1.000	1.000	1.000	0.107	7.280	3.944
$proposed(\tau =0.50)$	1.000	1.000	1.000	1.000	1.000	0.402	1.000	1.000	1.000	1.000	0.065	7.160	3.288
$proposed(\tau =0.75)$	1.000	1.000	1.000	1.000	1.000	0.400	1.000	1.000	1.000	1.000	0.137	7.360	3.912
$proposed(\tau =0.9)$	1.000	1.000	1.000	1.000	1.000	0.000	-	-	-	-	-	3.033	0.921
$t\left(3\right)$
$hirNet$	1.000	1.000	1.000	1.000	1.000	0.400	1.000	1.000	1.000	1.000	0.676	11.200	32.903
$RAMP$	1.000	1.000	0.700	0.680	1.000	0.000	1.000	1.000	1.000	1.000	0.696	9.360	70.609
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	1.000	0.000	-	-	-	-	-	3.320	3.572
$proposed(\tau =0.25)$	1.000	1.000	1.000	1.000	1.000	0.400	1.000	1.000	1.000	1.000	0.074	7.180	5.371
$proposed(\tau =0.50)$	1.000	1.000	1.000	1.000	1.000	0.400	1.000	1.000	1.000	1.000	0.056	7.160	5.568
$proposed(\tau =0.75)$	1.000	1.000	1.000	1.000	1.000	0.400	0.980	1.000	0.980	0.990	0.091	7.200	5.145
$proposed(\tau =0.9)$	1.000	1.000	1.000	1.000	1.000	0.000	-	-	-	-	-	3.200	3.560
${\chi }^2\left(2\right)$
$hirNet$	1.000	1.000	1.000	1.000	1.000	0.4000	1.000	1.000	1.000	1.000	0.679	11.260	34.494
$RAMP$	0.000	1.000	0.760	0.600	0.666	0.000	1.000	1.000	1.000	1.000	0.682	9.060	72.785
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	1.000	0.011	-	-	-	-	-	4.766	4.723
$proposed(\tau =0.25)$	1.000	1.000	1.000	1.000	1.000	0.400	0.900	1.000	0.900	0.950	0.251	7.620	5.462
$proposed(\tau =0.50)$	1.000	1.000	1.000	1.000	1.000	0.402	0.940	1.000	0.940	0.970	0.163	7.400	5.814
$proposed(\tau =0.75)$	1.000	1.000	1.000	1.000	1.000	0.400	0.920	1.000	0.920	0.960	0.150	7.340	6.005
$proposed(\tau =0.9)$	1.000	1.000	1.000	1.000	1.000	0.000	-	-	-	-	-	3.566	4.820
$\sigma \left({\mathbf{x}}_{\mathbf{i}}\right)e_i$
$hirNet$	1.000	1.000	1.000	1.000	1.000	0.402	1.000	1.000	1.000	1.000	0.673	11.140	30.643
$RAMP$	0.840	1.000	0.920	0.900	0.940	0.241	1.000	1.000	0.860	0.930	0.845	15.860	7.298
$proposed(\tau =0.1)$	1.000	1.000	1.000	1.000	1.000	0.000	-	-	-	-	-	3.440	1.055
$proposed(\tau =0.25)$	1.000	1.000	1.000	1.000	1.000	0.400	1.000	1.000	1.000	1.000	0.010	7.020	3.950
$proposed(\tau =0.50)$	1.00	1.000	1.000	1.000	1.000	0.400	1.000	1.000	1.000	1.000	0.029	7.060	3.948
$proposed(\tau =0.75)$	1.000	1.000	1.000	1.000	1.000	0.400	1.000	1.000	1.000	1.000	0.029	7.060	3.948
$proposed(\tau =0.9)$	1.000	1.000	1.000	1.000	1.000	0.000	-	-	-	-	-	3.166	1.052

Table 5

The covariates selected by the proposed methods.

Method	Covariates
RAMP	APQ11	HNR	RPDE	DFA	PPE
	HNR*RPDE	HNR*DFA	DFA*PPE
hirNet	PPQ5	APQ11	NHR	HNR	RPDE
	DFA	PPE	APQ11*PPE	APQ11*DFA	HNR*RPDE
	RPDE*PPE	DFA*PPE
Proposed method
$\tau =0.1$	HNR	RPDE	DFA	HNR*DFA
$\tau =0.25$	HNR	DFA	PPE	DFA*PPE
$\tau =0.5$	APQ11	HNR	RPDE	DFA	HNR*RPDE
$\tau =0.75$	HNR	DFA	PPE	HNR*DFA	HNR*PPE
$\tau =0.9$	dB	HNR

The symbol * represents the interaction between the two main effects.

Table 6

The average RMSE over the 100 data sets.

Method	RMSE	Size	Method	Size	RMSE	Method	$\mathit{\tau }$	Size	RMSE
							$\tau =0.1$	3	0.891
							$\tau =0.25$	4.7	0.873
RAMP	0.878	9.5	hirNet	0.887	18.39	proposed method	$\tau =0.5$	4.9	1.00
							$\tau =0.75$	5.0	0.858
							$\tau =0.9$	2.3	0.971

Appendix A. Proofof the Theorem

Appendix A.1. Proof of Theorem 1

Throughout the proof, if $v={(v_1,\dots ,v_k)}^{\top }$ is a vector, we use the norms ${\parallel v\parallel }_2=\sqrt{{\sum }_{j=1}^k{v}_j^2}$ and ${\parallel v\parallel }_{\infty }={max}_j\left|v_j\right|$. For a function f on $[0,1]$, we denote its $L_2\left(P\right)$ norm by ${\parallel f\parallel }_{L_2}=\sqrt{{\int }_0^1}{f}^2\left(x\right)dP\left(x\right)=\sqrt{E{\left(f\right)}^2}$. We write ${\mathit{\gamma }}_{01}={({\mathit{\beta }}_{01}^{\top },{\mathit{\beta }}_{001}^{\top })}^{\top }$ and write ${\tilde{\mathit{\gamma }}}_{01}={({\tilde{\mathit{\beta }}}_{01}^{\top },{\tilde{\mathit{\beta }}}_{001}^{\top })}^{\top }$ in the same fashion. Let $\mathit{\theta }=\sqrt{n}({\tilde{\mathit{\gamma }}}_{01}-{\mathit{\gamma }}_{01})$. Note that we can express ${\widehat{Q}}_{Y_i|{\mathbf{x}}_i}\left(\tau \right)={\mathbf{\Psi }}_{{\mathcal{M}}_i}^{\top }{\widehat{\mathit{\beta }}}_{01}+{\mathsf{\Phi }}_{{\mathcal{I}}_i}^{\top }{\widehat{\mathit{\beta }}}_{001}$; alternatively, we can express it as ${\widehat{Q}}_{Y_i|{\mathbf{x}}_i}\left(\tau \right)={\mathbf{\Psi }}_{{\mathcal{M}}_i}^{\top }{\tilde{\mathit{\beta }}}_{01}+{\mathsf{\Phi }}_{{\mathcal{I}}_i}^{\top }{\tilde{\mathit{\beta }}}_{001}$. By the identifiability of the model, we must have ${\widehat{\mathit{\gamma }}}_{01}={\tilde{\mathit{\gamma }}}_{01}$.

Notice that$${\displaystyle \frac{1}{n}}\sum _{i=1}^n{\rho }_{\tau }\left(Y_i-{\mathbf{\Psi }}_{{\mathcal{M}}_i}^{\top }{\tilde{\mathit{\beta }}}_{01}+{\mathsf{\Phi }}_{{\mathcal{I}}_i}^{\top }{\tilde{\mathit{\beta }}}_{001}\right)={\displaystyle \frac{1}{n}}\sum _{i=1}^n{\rho }_{\tau }\left({ϵ}_i-{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }-U_{ni}\right)$$ where ${\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}=n^{-1/2}{\mathsf{\Pi }}_{{\mathcal{A}}_i}$ and $U_{ni}={\mathsf{\Pi }}_{{\mathcal{A}}_i}^{\top }{\tilde{\mathit{\gamma }}}_{01}-m\left({\mathbf{x}}_i\right)$. Define the minimizers under the transformation as$$\widehat{\mathit{\theta }}=arg\underset{\mathbf{\theta }}{min}{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\rho }_{\tau }\left({ϵ}_i-{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }-U_{ni}\right).$$

Let $a_n$ be a sequence of positive numbers and define$$Q_i\left(a_n\right)={\rho }_{\tau }({ϵ}_i-a_n{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }-U_{ni}).$$

Define$$D_i(\mathit{\theta },a_n)=Q_i\left(a_n\right)-Q_i\left(0\right)-E\left[Q_i\left(a_n\right)-Q_i\left(0\right)|{\mathbf{x}}_i\right]+a_n{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }{\phi }_{\tau }\left({ϵ}_i\right),$$ and$${\tilde{Q}}_i(\mathit{\theta },a_n)=Q_i\left(a_n\right)-Q_i\left(0\right)+a_n{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }{\phi }_{\tau }\left({ϵ}_i\right),$$ where ${\phi }_{\tau }\left({ϵ}_i\right)=\tau -I({ϵ}_i<0)$.

Lemma A1. 

Let $q_n=q{L}_n+s{L}_n^2$. If Conditions (A1)–(A4) are satisfied, then, for any positive constant L, $$q_n^{-1}\underset{\parallel \mathit{\theta }\parallel \le L}{sup}\left|D_i(\mathit{\theta },\sqrt{q_n})\right|=o_p\left(1\right).$$

Proof. 

Note that$$\begin{array}{c}\hfill D_i(\mathit{\theta },\sqrt{q_n})={\tilde{Q}}_i(\mathit{\theta },\sqrt{q_n})-E\left[{\tilde{Q}}_i(\mathit{\theta },\sqrt{q_n})\right].\end{array}$$

Using Knight’s identity$${\rho }_{\tau }(r-s)-{\rho }_{\tau }\left(r\right)=-s\{\tau -I(r\le 0)\}+{\int }_0^s\{I(r\le t)-I(r\le 0)\}dt,$$ we have$$\begin{array}{c}\hfill \begin{array}{cc}{\tilde{Q}}_i(\mathit{\theta },\sqrt{q_n})\hfill & ={\rho }_{\tau }\left({ϵ}_i-\sqrt{q_n}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }-U_{ni}\right)-{\rho }_{\tau }\left({ϵ}_i-U_{ni}\right)+\sqrt{q_n}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }{\phi }_{\tau }\left({ϵ}_i\right)\hfill \\ & {\displaystyle ={\int }_0^{\sqrt{q_n}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }}\left[I({ϵ}_i-U_{ni}<t)-I({ϵ}_i-U_{ni}<0)\right]dt.}\hfill \end{array}\end{array}$$

Therefore, $Var\left(D_i(\mathit{\theta },\sqrt{q_n})\right)=Var\left({\tilde{Q}}_i(\mathit{\theta },\sqrt{q_n})\right)\le E\left({\tilde{Q}}_i^2(\mathit{\theta },\sqrt{q_n})\right)$. We have$$\begin{array}{ccc}& & \sum _{i=1}^{n}E\left({\tilde{Q}}_i^2(\mathit{\theta },\sqrt{q_n})|{\mathbf{x}}_i\right)\hfill \\ & \le & {\displaystyle C{q}_n{n}^{-1/2}\sum _{i=1}^n{\int }_0^{\sqrt{q_n}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }}{F}_i(t+U_{ni})-F_i\left(U_{ni}\right)dt}\hfill \\ & \le & {\displaystyle C{q}_n{n}^{-1/2}\sum _{i=1}^n{\int }_0^{\sqrt{q_n}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }}(f\left(0\right)+o\left(1\right))(t+o\left(t^2\right))d2}\hfill \\ & \le & C{q}_n^2{n}^{-1/2}\left[{\mathit{\theta }}^{\top }\left(\sum _{i=1}^n{f}_i\left(0\right){\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\right)\mathit{\theta }\right](1+o\left(1\right))\hfill \\ & \le & C{q}_n^2{n}^{-1/2}\left[{\parallel \mathit{\theta }\parallel }_2^2{\lambda }_{max}\left(n^{-1}{\mathsf{\Pi }}_{\mathcal{A}}^{\top }{B}_n{\mathsf{\Pi }}_{\mathcal{A}}\right)\right]\hfill \\ & \le & C{q}_n^2{n}^{-1/2}(1+o\left(1\right)).\hfill \end{array}$$ for some positive constant C, where $B_n=diag(f_1\left(0\right),\dots ,f_n\left(0\right))$ is an $n\times n$ diagonal matrix with $f_i\left(0\right)$ denoting the conditional density function of ${ϵ}_i$ given ${\mathbf{x}}_i$ evaluated at zero. Therefore, ${\displaystyle \sum _{i=1}^n}Var\left\{D_i(\mathbf{\theta },\sqrt{q_n})\right\}\le C{q}_n^2{n}^{-1/2}$ for some positive constant C and all n sufficiently large. By Bernstein’s inequality, for all n sufficiently large,$$\begin{array}{c}\hfill \begin{array}{cc}& P\left(q_n^{-1}|{\displaystyle \sum _{i=1}^n}{D}_i(\mathbf{\theta },\sqrt{q_n})|>\nu |{\mathbf{x}}_i\right)\hfill \\ & \le q_n^{-1}exp\left(-{\displaystyle \frac{{\nu }^2}{C{q}_n^2{n}^{-1/2}+C\nu n^{-1/2}}}\right)\hfill \\ & \le q_n^{-1}exp(-C{n}^{1/2}{q}_n^{-2})\hfill \end{array}\end{array}$$ which converges to 0 as $n\to \infty$ by Conditions (A3) and (A4). Note that the upper bound does not depend on ${\mathbf{x}}_i$, so the above bound also holds unconditionally. □

Lemma A2. 

Suppose that Conditions (A1)–(A4) hold. Then, for any sequence $\left\{b_n\right\}$ with $1\le b_n\le L_n^{\xi /10},0<\xi <(r-1/2)/(2r+1)$, we have $$\begin{array}{cc}\hfill \underset{{\mathit{\theta }}^{\top }{\tilde{\mathsf{\Pi }}}_{\mathcal{A}}^{\top }{\tilde{\mathsf{\Pi }}}_{\mathcal{A}}\mathit{\theta }\le b_n^2{L}_n}{sup}|\sum _{i=1}^n& \rho ({ϵ}_i-{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }-U_{ni})-\rho ({ϵ}_i-U_{ni})+{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }(\tau -I({ϵ}_i<0))\hfill \\ \hfill & -E_{{ϵ}_i|{\mathbf{x}}_i}(\rho ({ϵ}_i-{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }-U_{ni})-\rho ({ϵ}_i-U_{ni}))|=o_p\left(L_n\right).\hfill \end{array}$$

Using the similar arguments as described to prove Lemma 3.2 in [40], Lemma A2 can be proven.

Proof. 

For Theorem 1(1), we first prove that, $\forall \eta >0$, there exists a $C>0$ such that$$P\left(\underset{{\parallel \mathit{\theta }\parallel }_2\le L}{inf}{q}_n^{-1}\sum _{i=1}^n(Q_i\left(\sqrt{q_n}\right)-Q_i\left(0\right))>0\right)\ge 1-\eta .$$

Note that$$\begin{array}{c}\hfill \begin{array}{cc}{q}_n^{-1}{\displaystyle \sum _{i=1}^n}(Q_i\left(\sqrt{q_n}\right)-Q_i\left(0\right))\hfill & =q_n^{-1}{\displaystyle \sum _{i=1}^n}{D}_i(\mathit{\theta },\sqrt{q_n})+q_n^{-1}{\displaystyle \sum _{i=1}^n}E\left[Q_i\left(\sqrt{q_n}\right)-Q_i\left(0\right)\right]\hfill \\ & -q_n^{-1/2}{\displaystyle \sum _{i=1}^n}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }{\phi }_{\tau }\left({ϵ}_i\right)\hfill \\ & \triangleq G_{n1}+G_{n2}+G_{n3},\hfill \end{array}\end{array}$$ where the definition of $G_{ni},i=1,2,3$ is clear from the context. First, we can see that $\underset{{\parallel \mathit{\theta }\parallel }_2\le L}{inf}\left|G_{n1}\right|=O_p\left(1\right)$ by Lemma A1.

For $G_{n2}$, we have$$\begin{array}{ccc}& & q_n^{-1}\sum _{i=1}^{n}E\{Q_i\left(\sqrt{q_n}\right)-Q_i\left(0\right)\}\hfill \\ & =& q_n^{-1}\sum _{i=1}^{n}E\left({\rho }_{\tau }({ϵ}_i-\sqrt{q_n}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }-U_{ni})-{\rho }_{\tau }({ϵ}_i-U_{ni})\right)\hfill \\ & =& q_n^{-1}\sum _{i=1}^{n}E(-\sqrt{q_n}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }(\tau -I({ϵ}_i-U_{ni}\le 0))\hfill \\ & & {\displaystyle +{\int }_0^{\sqrt{q_n}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }}\left\{I({ϵ}_i-U_{ni}<s)-I({ϵ}_i-U_{ni}<0)\right\}ds|{\mathbf{x}}_i)}\hfill \\ & =& -q_n^{-1/2}\sum _{i=1}^{n}E\left({\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }(\tau -I({ϵ}_i-U_{ni}\le 0))\right)\hfill \\ & & {\displaystyle +q_n^{-1}\sum _{i=1}^{n}E\left({\int }_{U_{ni}}^{\sqrt{q_n}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }+U_{ni}}\left\{I({ϵ}_i<s)-I({ϵ}_i<0)\right\}ds|{\mathbf{x}}_i\right)}\hfill \\ & \triangleq & W_{n1}\left(\mathit{\theta }\right)+W_{n2}\left(\mathit{\theta }\right),\hfill \end{array}$$ where the definition of $W_{ni}\left(\mathit{\theta }\right),i=1,2$ is clear from the context. Note that $|F_{ϵ|\mathbf{x}}\left(0\right|\mathbf{x})-F_{ϵ|\mathbf{x}}\left(U_{ni}\right|\mathbf{x})|\le B|U_{ni}|$ for all $\mathbf{x}$, where B is the constant in the assumption (A2). Let $U_n={(u_{n1},\dots ,u_{nn})}^{\top }$. By Condition (A3), we have $\parallel U_n{\parallel }_2=O\left(L_n^{-r}\right)$. Consequently, we can take a constant $M>0$ such that ${sup}_{{\parallel \mathit{\theta }\parallel }_2\le L}|W_{n1}\left(\mathit{\theta }\right)|\le M{q}_n^{-1/2}\parallel n^{-1/2}{\mathit{\theta }}^{\top }{\mathsf{\Pi }}_{\mathcal{A}}^{\top }{\parallel }_2\parallel U_n{\parallel }_2=O_p\left(q_n^{-1/2}{n}^{-1/2}{L}_n^{-r}\right){\parallel \mathit{\theta }\parallel }_2=O_p{(\parallel \mathit{\theta }\parallel }_2)$ by Condition (A3) and Lemma A2. For $W_{n2}\left(\mathit{\theta }\right)$, we have$$\begin{array}{ccc}\hfill W_{n2}\left(\mathit{\theta }\right)& =& {\displaystyle q_n^{-1}\sum _{i=1}^n{\int }_{U_{ni}}^{\sqrt{q_n}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }+U_{ni}}{f}_i\left(0\right)sds(1+o\left(1\right))}\hfill \\ & =& q_n^{-1}\sum _{i=1}^n{f}_i\left(0\right)\left[{\displaystyle \frac{1}{2}}{q}_n{\left({\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }\right)}^2+U_{ni}\sqrt{q_n}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }\right](1+o\left(1\right))\hfill \\ & =& C{\mathit{\theta }}^{\top }\left(n^{-1}\sum _{i=1}^n{f}_i\left(0\right){\mathsf{\Pi }}_{{\mathcal{A}}_i}{\mathsf{\Pi }}_{{\mathcal{A}}_i}^{\top }\right)\mathit{\theta }\times (1+o\left(1\right))+q_n^{-1/2}\sum _{i=1}^n{f}_i\left(0\right)U_{ni}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }\hfill \\ & =& C{\mathit{\theta }}^{\top }{K}_n\mathit{\theta }\times (1+o\left(1\right))+q_n^{-1/2}\sum _{i=1}^n{f}_i\left(0\right)U_{ni}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }\hfill \end{array}$$ where $K_n={\displaystyle \frac{1}{n}}{\mathsf{\Pi }}_{\mathcal{A}}^{\top }{B}_n{\mathsf{\Pi }}_{\mathcal{A}}$. Based on Condition (A3), there exists a finite constant $c>0$, such that $C{\mathit{\theta }}^{\top }{K}_n\mathit{\theta }\times (1+o\left(1\right))\ge c{\parallel \mathit{\theta }\parallel }_2^2$ with the probability approaching one. Combining Condition (A2) and the Cauchy–Schwarz inequality, we obtain$$\begin{array}{ccc}\hfill q_n^{-1/2}\sum _{i=1}^n{f}_i\left(0\right)U_{ni}{\tilde{\mathsf{\Pi }}}_{{\mathcal{A}}_i}^{\top }\mathit{\theta }& =& q_n^{-1/2}{n}^{-1/2}{\mathit{\theta }}^{\top }{\mathsf{\Pi }}_{\mathcal{A}}^{\top }{B}_n{U}_n\hfill \\ & \le & q_n^{-1/2}\parallel n^{-1/2}{\mathit{\theta }}^{\top }{\mathsf{\Pi }}_{\mathcal{A}}^{\top }{\parallel }_2·{\parallel B_n{U}_n\parallel }_2\hfill \\ & =& O_p\left(q_n^{-1/2}{n}^{-1/2}{L}_n^{-r}\right){\parallel \mathit{\theta }\parallel }_2=O_p{(\parallel \mathit{\theta }\parallel }_2).\hfill \end{array}$$

We next evaluate $G_{n3}$, noting that $E\left(G_{n3}\right)=0$ and$$E\left(G_{n3}^2\right)\le C{q}_n^{-1}E\left({\mathit{\theta }}^{\top }\left(n^{-1}\sum _{i=1}^n{\mathsf{\Pi }}_{{\mathcal{A}}_i}{\mathsf{\Pi }}_{{\mathcal{A}}_i}^{\top }\right)\mathit{\theta }\right)=O(q_n^{-1}{\parallel \mathit{\theta }\parallel }_2^2).$$

Therefore, $G_{n3}=O(q_n^{-1}{\parallel \mathit{\theta }\parallel }_2^2)$. Hence, for L sufficiently large, the quadratic term will dominate and $q_n^{-1}{\sum }_{i=1}^n(Q_i\left(\sqrt{q_n}\right)-Q_i\left(0\right))$ has asymptotically a lower bound $c{L}^2$. By convexity, this implies $\parallel \widehat{\mathit{\theta }}{\parallel }_2=O_p\left(\sqrt{q_n}\right)$. From the definition of $\widehat{\mathit{\theta }}$, it follows that $\parallel {\widehat{\mathit{\gamma }}}_{01}-{\mathit{\gamma }}_{01}{\parallel }_2=O_p\left(\sqrt{n^{-1}{q}_n}\right)$. This completes the proof of Theorem 1 (a).

The proof of Theorem 1 (b) is immediate since $\parallel \widehat{m}{\left(\mathbf{x}\right)-m\left(\mathbf{x}\right)\parallel }_{L_2}={\parallel {\widehat{\mathit{\gamma }}}_{01}-{\mathit{\gamma }}_{01}\parallel }_2$ and ${sup}_{\mathbf{x}\in {[0,1]}^p}{\parallel \mathsf{\Pi }\parallel }_2=O\left(L_n\right)$. □

Appendix A.2. Proof of Theorem 2

Proof. 

Note that $\widehat{\mathit{\gamma }}={({\widehat{\mathit{\beta }}}_0^{\top },{\widehat{\mathit{\beta }}}_{00}^{\top })}^{\top }$ is the oracle estimator. Our goal is to prove that $\widehat{\mathit{\gamma }}$ is a local minimizer of $Q^P({\mathit{\beta }}_j,{\mathit{\beta }}_{jk})$. The gradient function of ${\sum }_{i=1}^n{\rho }_{\tau }(Y_i-{\sum }_{j=1}^p{\mathbf{\Psi }}_{i,j}^{\top }{\mathit{\beta }}_j-{\sum }_{1\le j<k\le p}{\mathsf{\Phi }}_{i,jk}^{\top }{\mathit{\beta }}_{jk})$ is not applicable in the proof because the check loss function ${\rho }_{\tau }$ is not differentiable at zero. We derive it directly from a certain lower bound of the difference of two check loss functions.

Suppose that there exists an index $j_0\in {\mathcal{M}}^c,$ $k_0\in {\mathcal{M}}^c$ and $(j_0,k_0)\in {\mathcal{I}}^c$, such that $f_{j_0}\ne 0$ and $f_{j_0{k}_0}\ne 0$. That is, ${\widehat{\beta }}_{j_0}\ne 0$ and ${\widehat{\beta }}_{j_0{k}_0}\ne 0$. Let ${\widehat{\mathit{\gamma }}}^{*}={({\widehat{\mathit{\beta }}}_0^{*\top },{\widehat{\mathit{\beta }}}_{00}^{*\top })}^{\top }$ be the vector obtained with ${\widehat{\beta }}_{j_0}\ne 0$ and ${\widehat{\beta }}_{j_0{k}_0}\ne 0$ being replaced by 0. Since ${\rho }_{\tau }\left(u\right)-{\rho }_{\tau }\left(v\right)\ge (\tau -I(v\le 0))(u-v)$ for any $u,v\in R$, then$$\begin{array}{ccc}& & Q^P({\widehat{\mathit{\beta }}}_j,{\widehat{\mathit{\beta }}}_{jk})-Q^P({\widehat{\mathit{\beta }}}_j^{*},{\widehat{\mathit{\beta }}}_{jk}^{*})\hfill \\ & \ge & -{\displaystyle \frac{1}{n}}\sum _{i=1}^n(\tau -I(Y_i\le {\mathsf{\Pi }}_i^{\top }{\widehat{\mathit{\gamma }}}^{*})){\mathsf{\Pi }}_i^{\top }(\widehat{\mathit{\gamma }}-{\widehat{\mathit{\gamma }}}^{*})\hfill \\ & & +{\lambda }_1({\left(\sum _{j=1}^p\parallel R_j{\parallel }_2\parallel {\mathit{\beta }}_{j_0}{\parallel }_2+\sum _{1\le k<j\le p}{I}_{{\beta }_{jk\ne 0}}(\parallel R_j{\parallel }_2\parallel {\mathit{\beta }}_{j_0}{\parallel }_2+\parallel R_k{\parallel }_2\parallel {\mathit{\beta }}_{k_0}{\parallel }_2)\right)}^{1/2}\hfill \\ & & +\sum _{1\le k<j\le p}\parallel Q_{jk}{\parallel }_2{\parallel {\mathit{\beta }}_{j_0{k}_0}\parallel }_2)\hfill \\ & =& -{\displaystyle \frac{1}{n}}\sum _{i=1}^n(\tau -I({ϵ}_i\le 0)){\mathsf{\Pi }}_i^{\top }(\widehat{\mathit{\gamma }}-{\widehat{\mathit{\gamma }}}^{*})-{\displaystyle \frac{1}{n}}\sum _{i=1}^n(I({ϵ}_i\le 0)-I({ϵ}_i\le r_{ni})){\mathsf{\Pi }}_i^{\top }(\widehat{\mathit{\gamma }}-{\widehat{\mathit{\gamma }}}^{*})\hfill \\ & & +{\lambda }_1({\left(\sum _{j=1}^p\parallel R_j{\parallel }_2\parallel {\mathit{\beta }}_{j_0}{\parallel }_2+\sum _{1\le k<j\le p}{I}_{{\beta }_{jk\ne 0}}(\parallel R_j{\parallel }_2\parallel {\mathit{\beta }}_{j_0}{\parallel }_2+\parallel R_k{\parallel }_2\parallel {\mathit{\beta }}_{k_0}{\parallel }_2)\right)}^{1/2}\hfill \\ & & +\sum _{1\le k<j\le p}\parallel Q_{jk}{\parallel }_2{\parallel {\mathit{\beta }}_{j_0{k}_0}\parallel }_2)\hfill \\ & \ge & -{\displaystyle \frac{1}{n}}\parallel \sum _{i=1}^n(\tau -I({ϵ}_i\le 0)){\mathsf{\Pi }}_i\parallel \parallel \widehat{\mathit{\gamma }}-{\widehat{\mathit{\gamma }}}^{*}\parallel -{\displaystyle \frac{1}{n}}\parallel \sum _{i=1}^n(I({ϵ}_i\le 0)-I({ϵ}_i\le r_{ni})){\mathsf{\Pi }}_i\parallel \parallel \widehat{\mathit{\gamma }}-{\widehat{\mathit{\gamma }}}^{*}\parallel \hfill \\ & & +{\lambda }_1({\left(\sum _{j=1}^p\parallel R_j{\parallel }_2\parallel {\mathit{\beta }}_{j_0}{\parallel }_2+\sum _{1\le k<j\le p}{I}_{{\beta }_{jk\ne 0}}(\parallel R_j{\parallel }_2\parallel {\mathit{\beta }}_{j_0}{\parallel }_2+\parallel R_k{\parallel }_2\parallel {\mathit{\beta }}_{k_0}{\parallel }_2)\right)}^{1/2}\hfill \\ & & +\sum _{1\le k<j\le p}\parallel Q_{jk}{\parallel }_2{\parallel {\mathit{\beta }}_{j_0{k}_0}\parallel }_2)\hfill \\ & \ge & -T_{n1}-T_{n2}+T_{n3}\hfill \end{array}$$ where $r_{ni}=U_{ni}-{\mathsf{\Pi }}_i^{\top }(\widehat{\mathit{\gamma }}-{\widehat{\mathit{\gamma }}}^{*})$. By simple calculation, one has that $\parallel {\sum }_{i=1}^n(\tau -I({ϵ}_i\le 0)){\mathsf{\Pi }}_i{\parallel }_2=O_p\left(n^{1/2}{L}_n\right)$. Therefore, $T_{n1}=O_p\left(n^{-1}{L}_n^2\right)$. For $T_{n2}$, from Conditions (A2) and (A3), we have$$\begin{array}{ccc}& & E{\left[\sum _{i=1}^n(I({ϵ}_i\le 0)-I({ϵ}_i\le r_{ni})){\mathsf{\Pi }}_i\right]}^2\hfill \\ & \le & n\sum _{i=1}^{n}E{\left[(I({ϵ}_i\le 0)-I({ϵ}_i\le r_{ni})){\mathsf{\Pi }}_i\right]}^2\hfill \\ & \le & n\sum _{i=1}^{n}E\left[{\mathsf{\Pi }}_i^{\top }{\mathsf{\Pi }}_i|I({ϵ}_i-r_{ni}\le 0)-I({ϵ}_i\le 0)|\right]\hfill \\ & \le & n\sum _{i=1}^{n}E\left[s_{\left(n\right)}^{2}I(0\le |{ϵ}_i|\le |r_{ni}\left|\right)\right]\hfill \\ & \le & L_n\sum _{i=1}^n{\int }_{-|r_{ni}|}^{|r_{ni}|}{f}_i\left(s\right)ds\hfill \\ & =& O_p\left(n^{1/2}{L}_n^3\right)\hfill \end{array}$$ where $s_{\left(n\right)}={max}_i{\parallel \phi \left({\mathbf{x}}_i\right)\parallel }_2\le c_1{n}^{-1/2}{L}_n$ for some positive constant $c_1$. The last equality follows by observing that $\underset{1\le i\le n}{max}|r_{ni}|\le O\left(L_n^{-r}\right)+L_n{\parallel \widehat{\mathit{\gamma }}-{\widehat{\mathit{\gamma }}}^{*}\parallel }_2=O_p\left(n^{-1/2}{L}_n^2\right)$. This implies that $\parallel {\sum }_{i=1}^n(\tau -I({ϵ}_i\le 0)){\mathsf{\Pi }}_i{\parallel }_2=O_p\left(n^{1/2}{L}_n^3\right)$. Therefore, $T_{n2}=O_p\left(n^{-1/2}{L}_n^3\right)$. By $T_{n1}$, $T_{n2}$, and $n^{-1}{L}_n^3{\lambda }_1^{-1}\to 0$, we have that $T_{n3}$ will dominate. Therefore, $Q^P({\widehat{\mathit{\beta }}}_j,{\widehat{\mathit{\beta }}}_{jk})-Q^P({\widehat{\mathit{\beta }}}_j^{*},{\widehat{\mathit{\beta }}}_{jk}^{*})>0$ with probability tending to one, which contradicts the fact that $\widehat{\mathit{\gamma }}$ is the minimizer of $Q^P({\mathit{\beta }}_j,{\mathit{\beta }}_{jk})$. Similarly, the same results are proven under weak heredity. This completes the proof of Theorem 2. □

Appendix A.3. Proof of Theorem 3

Lemma A3. 

Let $z_1,\dots ,z_n$ be independent random variable with values in some space $\mathcal{Z}$ and let Γ be a class of real-valued functions on $\mathcal{Z}$, satisfying, for some positive constants ${\eta }_n$ and ${\tau }_n$, $${\parallel \gamma \left(z\right)\parallel }_2\le {\eta }_n\mathit{and}{\displaystyle \frac{1}{n}}\sum _{i=1}^{n}var\left(\gamma \left(z_i\right)\right)\le {\tau }_n^2,\forall \gamma \in \mathsf{\Gamma }.$$Define $\mathbf{Z}:={sup}_{\gamma \in \mathsf{\Gamma }}\left|{\displaystyle \frac{1}{n}}{\sum }_{i=1}^n(\gamma \left(z_i\right)-E\gamma \left(z_i\right))\right|$. Then, for $t>0$, $$P\left(\mathbf{Z}\le E\left(\mathbf{Z}\right)+t\sqrt{2({\tau }_n^2+2{\eta }_{n}E\left(\mathbf{Z}\right))}+{\displaystyle \frac{2{\eta }_n{t}^2}{3}}\right)\le exp(-n{t}^2).$$

For the details of this proof, see [41].

Proof. 

Define $\mathsf{\Gamma }=\{\gamma \left(z\right):\gamma \left(z\right)=\rho (Y-\widehat{m}\left(x\right))-\rho (Y-\tilde{m}\left(x\right))\}$. We can write $[\mathcal{R}\left(\widehat{m}\right)-\mathcal{R}\left(\tilde{m}\right)]-[{\mathcal{R}}_n\left(\widehat{m}\right)-{\mathcal{R}}_n\left(\tilde{m}\right)]=E\gamma \left(z\right)-{\displaystyle \frac{1}{n}}{\sum }_{i=1}^n\gamma \left(z_i\right),\gamma \in \mathsf{\Gamma }$. By Lemma A3, we have$$\mathbf{Z}\le E\left(\mathbf{Z}\right)+\sqrt{2t^2({\tau }_n^2+2{\eta }_{n}E\left(\mathbf{Z}\right))}+{\displaystyle \frac{2{\eta }_n{t}^2}{3}}.$$ with probability at least $1-exp(-n{t}^2)$ for $t>0$. Based on the subadditivity and inequality $\sqrt{xy}\le (x+y)/2,x,y\ge 0$, we have$$\sqrt{2t^2({\tau }_n^2+2{\eta }_{n}E\left(\mathbf{Z}\right))}\le \sqrt{2t^2{\tau }_n^2}+2\sqrt{t^2{\eta }_{n}E\left(\mathbf{Z}\right)}\le \sqrt{2t^2{\tau }_n^2}+2E\left(\mathbf{Z}\right)+t^2{\eta }_n.$$

Then, with probability at least $1-exp(-n{t}^2)$, we have(A1)$$\mathbf{Z}\le 3E\left(\mathbf{Z}\right)+\sqrt{2t^2{\tau }_n^2}+{\displaystyle \frac{5t^2{\eta }_n}{3}}.$$

Let $\gamma \left(x\right)$ be the collection of all differences $\rho (Y-\widehat{m}\left(x\right))-\rho (Y-\tilde{m}\left(x\right))$. The bracketing number $N_{\left[\right]}(\delta ,{\mathcal{M}}_n)$ is the minimum number of $\epsilon$-brackets $[l_i,u_i]$ over ${\mathcal{M}}_n$, where $\parallel u_j-l_j{\parallel }_{\infty }\le ϵ,1\le j\le k$. We can construct 2$\epsilon$-brackets over $\gamma \left(x\right)$ by taking difference $[l_i-u_j,u_i-l_j]$ for the upper and lower bounds. Therefore, the bracketing numbers $N_{\left[\right]}(ϵ,\gamma \left(x\right))$ are bounded by the squares of the bracketing numbers $N_{\left[\right]}(ϵ/2,{\mathcal{M}}_n)$. For $\delta >0$, by theorem 19.5 in [42], there exists a finite number $a\left(\delta \right)$ such that$$E\left(\underset{i}{sup}\underset{\gamma \in \mathsf{\Gamma }}{sup}|{\displaystyle \frac{1}{n}}\sum _{i=1}^n\gamma \left(x_i\right)-E\gamma \left(z\right)|\right)\lesssim {\displaystyle \frac{J_{\left[\right]}(\delta ,{\mathcal{M}}_n)}{\sqrt{n}}}+M1\{M>a\left(\delta \right)\sqrt{n}\}.$$ where $J_{\left[\right]}(\delta ,{\mathcal{M}}_n)={\int }_0^{\delta }\sqrt{log{N}_{\left[\right]}(ϵ,{\mathcal{M}}_n)}dϵ$ is the bracketing integral. The envelope function M can be taken as equal to the supremum of the absolute values of the upper and lower bounds of finitely many brackets over ${\mathcal{M}}_n$. Based on Theorem 19.5 in [42], the second term on the right is bounded by $a{\left(\delta \right)}^{-1}P{M}^{2}1\{M>a\left(\delta \right)\sqrt{n}\}$ and hence converges to zero as $n\to \infty$. Given $K>0$, there exists a constant K such that, for Sobolev space $S_C^2\left([0,1]\right)$,$$log{N}_{\left[\right]}(\delta ,S_C^2\left([0,1]\right))\le K{\left({\displaystyle \frac{1}{\delta }}\right)}^{1/2},$$Note that, ${\sum }_{j\in \mathcal{M}}{\mathit{\phi }}_j^{\top }{\mathit{\beta }}_j\in {\mathcal{M}}_n\left(L_n\right)$ and ${\sum }_{\left(jk\right)\mathcal{I}}{\mathit{\varphi }}_{jk}^{\top }{\mathit{\beta }}_{jk}\in {\mathcal{M}}_n\left(L_n\right)$, where ${\mathit{\phi }}_j\in S_C^2\left([0,1]\right)$ and ${\mathit{\varphi }}_{jk}\in S_C^2\left({[0,1]}^2\right)$. Reference [43] implies that the bracketing integral of Sobolev space $S_C^2\left([0,1]\right)$ and $S_C^2\left({[0,1]}^2\right)$ is bounded. Then,$$J_{\left[\right]}(\delta ,{\mathcal{M}}_n)=O({(logp)}^{1/2}+{(2logp)}^{1/4})=O\left(\sqrt{logp}\right).$$

By the convexity of ${\rho }_{\tau }$, there exists $C>0$ such that $E{\left(\gamma \left(z\right)\right)}^2\le 2{\parallel \widehat{m}\left(x\right)-\tilde{m}\left(x\right)\parallel }_{L_2}^2$ and ${\parallel \gamma \parallel }_{\infty }\le C\parallel \widehat{m}\left(x\right)-\tilde{m}{\left(x\right)\parallel }_{\infty }\le {\parallel \widehat{m}\left(x\right)-\tilde{m}\left(x\right)\parallel }_{L_2}$. Let ${\tau }_n^2=\parallel \widehat{m}\left(x\right)-\tilde{m}{\left(x\right)\parallel }_{L_2}^2\le \parallel \widehat{m}{\left(x\right)-m\left(x\right)\parallel }_{L_2}^2+{\parallel \tilde{m}\left(x\right)-m\left(x\right)\parallel }_{L_2}^2=O\left(n^{-(2r-1)/(2r+1)}\right)$ and ${\eta }_n=\parallel \widehat{m}\left(x\right)-\tilde{m}{\left(x\right)\parallel }_{L_2}\le \parallel \widehat{m}{\left(x\right)-m\left(x\right)\parallel }_{L_2}+{\parallel \tilde{m}\left(x\right)-m\left(x\right)\parallel }_{L_2}=O\left(n^{-(2r-1)/2(2r+1)}\right)$. Then, Equation (A1) implies that$$\begin{array}{ccc}\hfill \mathbf{Z}& \le & n^{-1/2}\sqrt{logp}+\sqrt{2t^2}\parallel \widehat{m}\left(x\right)-\tilde{m}{\left(x\right)\parallel }_{L_2}+{\displaystyle \frac{5t^2}{3}}{\parallel \widehat{m}\left(x\right)-\tilde{m}\left(x\right)\parallel }_{L_2}\hfill \\ & \le & \sqrt{logp}{n}^{-1/2}+(\sqrt{2t^2}+5t^2/3)n^{-(2r-1)/2(2r+1)}\hfill \\ & \le & (2\sqrt{2t^2}+10t^2/3)n^{-(2r-1)/2(2r+1)}\hfill \end{array}$$ with probability at least $1-exp(-n{t}^2)$ as $n\to \infty$. From the above, we have $R\left(\widehat{m}\right)-R\left(\tilde{m}\right)\le {\mathcal{R}}_n\left(\widehat{m}\right)-{\mathcal{R}}_n\left(\tilde{m}\right)+(2\sqrt{2t^2}+10t^2/3)n^{-(2r-1)/2(2r+1)}$ and that there exists $D_n$ such that $P(\mathcal{R}\left(\widehat{m}\right)-\mathcal{R}\left(\tilde{m}\right)\le D_n)=P({\mathcal{R}}_n\left(\widehat{m}\right)-{\mathcal{R}}_n\left(\tilde{m}\right)+(2\sqrt{2t^2}+10t^2/3)n^{-(2r-1)/2(2r+1)})$. Under the assumption that the regression function is bounded, it follows, for ${\mathit{\phi }}_j\in S_C^2\left([0,1]\right)$ and ${\mathit{\varphi }}_{jk}\in S_C^2\left({[0,1]}^2\right)$, that$${\left(\sum _{j=1}^p\parallel {\mathit{\beta }}_j{\parallel }_2+\sum _{1\le k<j\le p}{I}_{{\beta }_{jk\ne 0}}(\parallel {\mathit{\beta }}_j{\parallel }_2+\parallel {\mathit{\beta }}_k{\parallel }_2)\right)}^{1/2}+\sum _{1\le k<j\le p}{\parallel {\mathit{\beta }}_{jk}\parallel }_2\le L_n$$ where $L_n=o\left({[n/log\left(n\right)]}^{1/4}\right)$. Reference [13] also shows that the Lasso is persistent when p grows polynomially in n. Furthermore, according to the definition of $\widehat{m}\left(x\right)$, we have ${\mathcal{R}}_n\left(\widehat{m}\right)-{\mathcal{R}}_n\left(\tilde{m}\right)\le {\lambda }_1{L}_n$ and $\mathcal{R}\left(\widehat{m}\right)-\mathcal{R}\left(\tilde{m}\right)={\lambda }_1{L}_n+(2\sqrt{2t^2}+10t^2/3)n^{-(2r-1)/2(2r+1)}$. Let $D_n={\lambda }_1{L}_n+(2\sqrt{2t^2}+10t^2/3)n^{-(2r-1)/2(2r+1)}$. Since ${\mathcal{R}}_n\left(\widehat{m}\right)-{\mathcal{R}}_n\left(\tilde{m}\right)\le {\lambda }_1{L}_n$ always holds, the probability of $P(\mathcal{R}\left(\widehat{m}\right)-\mathcal{R}\left(\tilde{m}\right)\le D_n)=P({\mathcal{R}}_n\left(\widehat{m}\right)-{\mathcal{R}}_n\left(\tilde{m}\right)+(2\sqrt{2t^2}+10t^2/3)n^{-(2r-1)/2(2r+1)}\le {\lambda }_1{L}_n+(2\sqrt{2t^2}+10t^2/3)n^{-(2r-1)/2(2r+1)})$ is 1. Note that ${\lambda }_1{L}_n+(2\sqrt{2t^2}+10t^2/3)n^{-(2r-1)/2(2r+1)}=O\left((2\sqrt{2t^2}+10t^2/3)n^{-(2r-1)/2(2r+1)}\right)$. Therefore, $\mathcal{R}\left(\widehat{m}\right)-\mathcal{R}\left(\tilde{m}\right)\le (2\sqrt{2t^2}+10t^2/3)n^{-(2r-1)/2(2r+1)}$ with probability at least $1-exp(-n{t}^2)$. Similarly, the same results are proved under weak heredity. □

Appendix B. Estimated Results at τ = 0.1, 0.25, 0.75, 0.9 in Applications

Figure A1 illustrates the effects of HNR, RPDE, and DFA on the UPDRS score at $\tau =0.1$. As HNR increases, the UPDRS score significantly decreases, showing a nonlinear relationship, where higher noise-to-harmonic ratios are associated with milder symptoms. RPDE exhibits a clear U-shaped curve, indicating that extreme values correspond to more severe symptoms, while intermediate values are associated with milder symptoms, reflecting a nonlinear effect pattern. DFA shows a negative correlation, with UPDRS scores decreasing as DFA increases, again revealing a nonlinear trend, suggesting that increased fundamental frequency variability may help to alleviate symptoms.

Figure A2 shows the interaction effect between HNR and DFA at $\tau =0.1$, revealing a nonlinear relationship. When HNR is low and DFA is high, the UPDRS score is higher; as HNR increases, the UPDRS score decreases. However, when HNR is high and DFA is low, the score increases again. This suggests that the effects of HNR and DFA on the UPDRS score are complex and interdependent, highlighting the importance of understanding and analyzing the interactions between these variables.

Figure A1 Estimated main effects for the HNR, RPDE, and DFA variables at $\tau =0.10$.

Figure A2 Estimated interaction between HNR and DFA at $\tau =0.1$.

Figure A3 illustrates the estimated main effect trends of three acoustic features—HNR, DFA, and PPE—at $\tau =0.25$. Under this condition, as HNR increases, its impact on the analysis first rises slightly and then drops rapidly, indicating that the purity of the speech signal influences recognition but with diminishing returns beyond a certain point. DFA reveals an optimal fluctuation pattern, after which its ability to reflect health status declines. PPE shows that the randomness and complexity of the speech signal also have an optimal range, beyond which the effect weakens. Overall, the variations in these features highlight their potential applications in fields such as disease diagnosis, where analyzing these acoustic characteristics can assist in medical diagnostics and improve accuracy.

Figure A4 illustrates the estimated interaction between DFA and PPE at $\tau =0.25$, presented through a 3D surface plot. As DFA and PPE values vary, their interaction significantly impacts the overall effect. Specifically, when DFA is at a low level, the overall effect initially increases slowly and then drops rapidly as PPE increases. In contrast, when DFA is high, this trend becomes more gradual, indicating a more complex interaction pattern. This suggests that the influence of their interaction varies significantly across different DFA and PPE combinations. For instance, in certain combinations, their interaction can more accurately reflect changes in health status within speech signals, which is important for disease diagnosis and improving speech processing technologies. Overall, the figure highlights the nonlinear interaction between DFA and PPE and its potential applications.

Figure A3 Estimated main effects for the HNR, DFA, and PPE variables at $\tau =0.25$.

Figure A4 Estimated interaction between DFA and PPE at $\tau =0.25$.

Figure A5 illustrates the estimated main effects of three variables HNR, DFA, and PPE at $\tau =0.75$. From the figure, it can be observed that, as the HNR value increases, its impact on the overall effect first rises and then rapidly decreases. The change in DFA values shows a similar trend, but the decline is more gradual. Meanwhile, the PPE value exhibits a distinct peak before gradually decreasing. These variations in the features highlight their significance in practical applications, particularly in disease diagnosis and speech processing fields.

Figure A6 illustrates the estimated interaction effects between HNR and DFA, as well as between HNR and PPE at $\tau =0.75$. From the figure, it can be observed that the interaction between HNR and DFA exhibits a complex nonlinear relationship, with the overall effect first increasing and then decreasing as the HNR value rises, showing a distinct fluctuation pattern. Similarly, the interaction between HNR and PPE also displays a complex pattern, but with more pronounced changes, particularly when the PPE value is high, leading to more significant variations in the overall effect.

Figure A5 Estimated main effects for the HNR and DFA variables at $\tau =0.75$.

Figure A6 Estimated interaction effects between HNR and DFA, as well as the interaction effect between HNR and PPE variables at $\tau =0.75$.

Figure A7 provides the effects of dB and HNR on the UPDRS score at $\tau =0.9$. The main effect of dB presents a U-shaped curve, indicating that both extremely high and low dB values lead to higher UPDRS scores, with lower scores at intermediate values. This emphasizes the importance of moderate changes in sound levels and their nonlinear effects. HNR also shows an inverted U-shaped curve, indicating that moderate HNR values are associated with the most severe symptoms, with this nonlinear relationship being particularly prominent at higher quantiles. These findings highlight the complex nonlinear effects of different variables on the UPDRS score, providing valuable insights for understanding Parkinson’s disease symptoms.

Figure A7 Estimated main effects for the dB and HNR variables at $\tau =0.90$.

References

1. Koenker, R.; Bassett, G. Regression quantiles. Econometrica; 1978; 46, pp. 33-50. [DOI: https://dx.doi.org/10.2307/1913643]

2. Nelder, J.A. A reformulation of linear models. J. R. Stat. Soc. Ser. A (Gen.); 1977; 140, pp. 48-63. [DOI: https://dx.doi.org/10.2307/2344517]

3. McCullagh, P.; Nelder, J. Monographs on Statistics and Applied Probability; Chapman & Hall: London, UK, 1989.

4. De Gooijer, J.G.; Zerom, D. On additive conditional quantiles with high-dimensional covariates. J. Am. Stat. Assoc.; 2003; 98, pp. 135-146. [DOI: https://dx.doi.org/10.1198/016214503388619166]

5. Cheng, Y.; De Gooijer, J.G.; Zerom, D. Efficient estimation of an additive quantile regression model. Scand. J. Stat.; 2011; 38, pp. 46-62. [DOI: https://dx.doi.org/10.1111/j.1467-9469.2010.00706.x]

6. Lee, Y.K.; Mammen, E.; Park, B.U. Backfitting and smooth backfitting for additive quantile models. Ann. Stat.; 2010; 38, pp. 2857-2883. [DOI: https://dx.doi.org/10.1214/10-AOS808]

7. Horowitz, J.L.; Lee, S. Nonparametric estimation of an additive quantile regression model. J. Am. Stat. Assoc.; 2005; 100, pp. 1238-1249. [DOI: https://dx.doi.org/10.1198/016214505000000583]

8. Sherwood, B.; Maidman, A. Additive nonlinear quantile regression in ultra-high dimension. J. Mach. Learn. Res.; 2022; 23, pp. 1-47.

9. Zhao, W.; Li, R.; Lian, H. Estimation and variable selection of quantile partially linear additive models for correlated data. J. Stat. Comput. Simul.; 2024; 94, pp. 315-345. [DOI: https://dx.doi.org/10.1080/00949655.2023.2245098]

10. Cui, X.; Zhao, W. Pursuit of dynamic structure in quantile additive models with longitudinal data. Comput. Stat. Data Anal.; 2019; 130, pp. 42-60. [DOI: https://dx.doi.org/10.1016/j.csda.2018.08.017]

11. Lin, Y.; Zhang, H.H. Component selection and smoothing in multivariate nonparametric regression. Ann. Stat.; 2006; 34, pp. 2272-2297. [DOI: https://dx.doi.org/10.1214/009053606000000722]

12. Storlie, C.B.; Bondell, H.D.; Reich, B.J.; Zhang, H.H. Surface estimation, variable selection, and the nonparametric oracle property. Stat. Sin.; 2011; 21, pp. 679-705. [DOI: https://dx.doi.org/10.5705/ss.2011.030a] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21603586]

13. Radchenko, P.; James, G.M. Variable selection using adaptive nonlinear interaction structures in high dimensions. J. Am. Stat. Assoc.; 2010; 105, pp. 1541-1553. [DOI: https://dx.doi.org/10.1198/jasa.2010.tm10130]

14. Efron, B.; Hastie, T.; Johnstone, I.; Tibshirani, R. Least angle regression. Ann. Stat.; 2004; 32, pp. 407-499. [DOI: https://dx.doi.org/10.1214/009053604000000067]

15. Wu, T.T.; Chen, Y.F.; Hastie, T.; Sobel, E.; Lange, K. Genome-wide association analysis by lasso penalized logistic regression. Bioinformatics; 2009; 25, pp. 714-721. [DOI: https://dx.doi.org/10.1093/bioinformatics/btp041]

16. Zhao, P.; Rocha, G.; Yu, B. The composite absolute penalties family for grouped and hierarchical variable selection. Ann. Stat.; 2009; 37, pp. 3468-3497. [DOI: https://dx.doi.org/10.1214/07-AOS584]

17. Yuan, M.; Joseph, V.R.; Zou, H. Structured variable selection and estimation. Ann. Appl. Stat.; 2009; 3, pp. 1738-1757. [DOI: https://dx.doi.org/10.1214/09-AOAS254]

18. Choi, N.H.; Li, W.; Zhu, J. Variable selection with the strong heredity constraint and its oracle property. J. Am. Stat. Assoc.; 2010; 105, pp. 354-364. [DOI: https://dx.doi.org/10.1198/jasa.2010.tm08281]

19. Bien, J.; Taylor, J.; Tibshirani, R. A lasso for hierarchical interactions. Ann. Stat.; 2013; 41, pp. 1111-1141. [DOI: https://dx.doi.org/10.1214/13-AOS1096]

20. She, Y.; Wang, Z.F.; Jiang, H. Group regularized estimation under structural hierarchy. J. Am. Stat. Assoc.; 2018; 113, pp. 445-454. [DOI: https://dx.doi.org/10.1080/01621459.2016.1260470]

21. Hao, N.; Feng, Y.; Zhang, H.H. Model Selection for High Dimensional Quadratic Regression via Regularization. J. Am. Stat. Assoc.; 2018; 113, pp. 615-625. [DOI: https://dx.doi.org/10.1080/01621459.2016.1264956]

22. Liu, C.; Ma, J.; Amos, C.I. Bayesian variable selection for hierarchical gene–environment and gene–gene interactions. Hum. Genet.; 2015; 134, pp. 23-36. [DOI: https://dx.doi.org/10.1007/s00439-014-1478-5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25154630]

23. Yang, Y.; Basu, S.; Zhang, L. A Bayesian hierarchical variable selection prior for pathway-based GWAS using summary statistics. Stat. Med.; 2020; 39, pp. 724-739. [DOI: https://dx.doi.org/10.1002/sim.8442] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31777110]

24. Kim, J.E.A. Bayesian variable selection with strong heredity constraints. J. Korean Stat. Soc.; 2018; 47, pp. 314-329. [DOI: https://dx.doi.org/10.1016/j.jkss.2018.03.003]

25. Li, Y.; Wang, N.; Carroll, R.J. Generalized functional linear models with semiparametric single-index interactions. J. Am. Stat. Assoc.; 2010; 105, pp. 621-633. [DOI: https://dx.doi.org/10.1198/jasa.2010.tm09313]

26. Li, Y.; Liu, J.S. Robust variable and interaction selection for logistic regression and multiple index models. J. Am. Stat. Assoc.; 2019; 114, pp. 271-286. [DOI: https://dx.doi.org/10.1080/01621459.2017.1401541]

27. Liu, Y.; Li, Y.; Carroll, R.J. Predictive functional linear models with diverging number of semiparametric single-index interactions. J. Econom.; 2021; 230, pp. 221-239. [DOI: https://dx.doi.org/10.1016/j.jeconom.2021.03.010]

28. Liu, H.; You, J.; Cao, J. A Dynamic Interaction Semiparametric Function-on-Scalar Model. J. Am. Stat. Assoc.; 2021; 118, pp. 360-373. [DOI: https://dx.doi.org/10.1080/01621459.2021.1933496]

29. Yu, K.; Lu, Z. Local linear additive quantile regression. Scand. J. Stat.; 2004; 31, pp. 333-346. [DOI: https://dx.doi.org/10.1111/j.1467-9469.2004.03_035.x]

30. Noh, H.; Lee, E.R. Component selection in additive quantile regression models. J. Korean Stat. Soc.; 2014; 43, pp. 439-452. [DOI: https://dx.doi.org/10.1016/j.jkss.2014.01.002]

31. Wood, S.N. P-splines with derivative based penalties and tensor product smoothing of unevenly distributed data. Stat. Comput.; 2017; 27, pp. 985-989. [DOI: https://dx.doi.org/10.1007/s11222-016-9666-x]

32. Eilers, P.H.C.; Marx, B.D. Flexible smoothing with B-splines and penalties. Stat. Sci.; 1996; 11, pp. 89-102. [DOI: https://dx.doi.org/10.1214/ss/1038425655]

33. Meier, L.; van de Geer, S.; Buehlmann, P. High-dimensional additive modeling. Ann. Stat.; 2009; 37, pp. 3779-3821. [DOI: https://dx.doi.org/10.1214/09-AOS692]

34. Fan, Y.; Tang, C.Y. Tuning parameter selection in high dimensional penalized likelihood. J. R. Stat. Soc. Ser. B; 2013; 75, pp. 531-552. [DOI: https://dx.doi.org/10.1111/rssb.12001]

35. Greenshtein, E.; Ritov, Y. Persistence in high-dimensional linear predictor selection and the virtue of overparametrization. Bernoulli; 2004; 10, pp. 971-988. [DOI: https://dx.doi.org/10.3150/bj/1106314846]

36. Jiang, L.; Bondell, H.; Wang, H. Interquantile shrinkage and variable selection in quantile regression. Comput. Stat. Data Anal.; 2014; 69, pp. 208-219. [DOI: https://dx.doi.org/10.1016/j.csda.2013.08.006]

37. Zou, H.; Ming, Y. Composite quantile regression and the oracle model selection theory. Ann. Stat.; 2008; 36, pp. 1108-1126. [DOI: https://dx.doi.org/10.1214/07-AOS507]

38. Kohns, D.; Szendrei, T. Horseshoe prior Bayesian quantile regression. J. R. Stat. Soc. Ser. C Appl. Stat.; 2024; 73, pp. 193-220. [DOI: https://dx.doi.org/10.1093/jrsssc/qlad091]

39. Tsanas, A.; Little, M.A.; McSharry, P.E.; Ramig, L.O. Accurate telemonitoring of parkinson’s disease progression by noninvasive speech tests. IEEE Trans. Bio-Med. Eng.; 2010; 57, pp. 884-893. [DOI: https://dx.doi.org/10.1109/TBME.2009.2036000]

40. He, X.M.; Shi, P. Convergence rate of B-spline estimators of nonparametric conditional quantile functions. J. Nonparametr. Stat.; 1994; 3, pp. 299-308. [DOI: https://dx.doi.org/10.1080/10485259408832589]

41. Bousquet, O. A Bennett concentration inequality and its application to suprema of empirical processes. Comptes Rendus Math.; 2002; 334, pp. 495-500. [DOI: https://dx.doi.org/10.1016/S1631-073X(02)02292-6]

42. Vaart, A. Asymptotic Statistics; Cambridge University Press: Cambridge, UK, 1990.

43. Birman, M.; Solomjak, M.Z. Piecewise-polynomial approximations of functions of the classes. Math. USSR-Sb.; 1967; 2, pp. 295-317. [DOI: https://dx.doi.org/10.1070/SM1967v002n03ABEH002343]