© Wuhan University 2025

0 Introduction

Subgroups are usually defined based on biomarkers such as blood sugar, blood pressure, gene or protein expression levels, etc^[1]. When evaluating the efficacy of new therapy on patients in clinical studies, these characteristic differences may lead to different responses of the same therapy to different patients^[2]. Therefore, when physicians assess therapeutic effect, they should not only consider the average effect across the overall population, but also the heterogeneous effects within subgroups. Many methods for subgroup identification have proposed in research literature, such as method based on the tree structure^[3], global model^[4], clustering analysis^[5], etc. In clinical practice, a straightforward and easily interpretable approach is often preferred, which involves dividing subgroups according to whether a single continuous biomarker exceeds a certain threshold. This is the threshold model we focus on.

However, in the threshold model described above, a single variable is used to classify subgroups, which can result in a limited representation of subgroup characteristics and may be compromised by the inappropriate choice of the threshold variable. Numerous studies suggest that subgroups may arise from the combined influence of multiple variables, for instance, Jiang et al^[6] demonstrated that multiple single nucleotide polymorphisms (SNPs) can frequently alter the risk of developing specific diseases, whereas any single SNP alone lacks this capability; Fan et al^[7] showed that a risk score defined as a function of multiple predictors is useful for identifying subgroups of AIDS patients who benefit more from treatment. Similarly, Vander Weele et al^[8] discussed the importance of using multiple variables to identify optimal treatment subgroups. He et al^[9] and Wei et al^[10] extended this model to survival and longitudinal data, respectively. While they utilized penalization methods to identify significant biomarkers, they did not perform variable selection for the included covariates.

Moreover, the majority of current studies on threshold model primarily focus on traditional mean regression. However, sometimes we place greater emphasis on sample data at other quantiles, for example, we prefer to conduct heterogeneity analysis on patients with worse condition and provide targeted medical solutions. Consequently, we consider introducing quantile regression to analyze the threshold model. The threshold quantile regression was initially introduced by Cai et al^[11] and applied to quantile self-exciting threshold autoregressive time series models. Galvao et al^[12] studied the threshold quantile autoregressive processes. Furthermore, Lee et al^[13] developed the general tests for threshold effects in regression models, Su and Xu^[14] presented a systematic estimation and inference procedure for threshold quantile regression. On this basis, Zhang et al^[15] proposed a more general single-index threshold quantile regression model. However, there is no existing literature that addresses the variable selection for covariates and biomarkers in single-index threshold quantile regression.

In practical applications, incorporating excessive number of biomarkers for subgroup division may diminish the interpretability of the subgroups, and including excessive number of covariates may also introduce multicollinearity. Consequently, we prefer to select the variables most pertinent to the response from the model's explanatory variables. Common variable selection methods include those based on penalties. This method, initially proposed by Tibshirani in 1996^[16], is known as Lasso estimation, which uses an L1 norm penalty. Subsequent researches include: Fused Lasso penalty^[17], Adaptive Lasso^[18], Graphical Lasso^[19] and so on. Furthermore, as Lasso compresses all coefficients equally, the resulting estimator is often biased. To this end, Fan et al^[20] introduced the smoothly clipped absolute deviation (SCAD) penalty and proved that the estimation result satisfied the oracle property.

Our paper extends the single-index threshold quantile regression^[15] to longitudinal data and introduces the SCAD penalty for variable selection of covariates and biomarkers. To this end, we propose an efficiency algorithm similar to that proposed by Zhang et al^[21], which transforms the problem of estimating the covariate regression coefficients and threshold parameters into a penalized linear quantile regression. Furthermore, we obtain a degraded linear quantile regression through pseudo observations. In addition, the penalty for covariates and biomarkers should be performed separately. We combine the two-step grid search algorithm^[22] to select the compression parameters that minimize the Bayesian information criterion (BIC), and then substitute the pre-selected compression parameters into the iterative algorithm for parameter estimation. Finally, random simulations and example analysis illustrate the practicality of our method.

1 Models and Estimation

We assume that $y_{ij}$ is the response variable of the i-th individual at the j-th observation, $i=\mathrm{1},\cdots ,n$, $j=\mathrm{1},\cdots ,m$, total observations $N=nm$, ${\mathit{z}}_{ij}={\left(z_{\mathrm{1},ij},\cdots ,z_{p,ij}\right)}^{\mathrm{T}}$ is the covariates with dimension of $p$, where $z_{\mathrm{1},ij}=\mathrm{1}$, ${\tilde{\mathit{z}}}_{ij}$ is a subset of ${\mathit{z}}_{ij}$, and x_i=${\left(x_{i\mathrm{0}},\cdots ,x_{id}\right)}^{\mathrm{T}}$ is the baseline biomarker including the intercept term. We consider the following single-index threshold quantile regression with$\tau \in \left(\mathrm{0,1}\right):y_{ij}={\mathit{z}}_{ij}^{\mathrm{T}}\mathit{\beta }\left(\tau \right)+{\tilde{\mathit{z}}}_{ij}^{\mathrm{T}}\mathit{\delta }\left(\tau \right)I\left({\mathit{x}}_i^{\mathrm{T}}\mathit{\gamma }\left(\tau \right)>\mathrm{0}\right)+{\epsilon }_{ij},$

$i=\mathrm{1},\cdots ,n,j=\mathrm{1},\cdots ,m.$

The regression coefficients $\mathit{\beta }\left(\tau \right)$ characterize the baseline effect of ${\mathit{z}}_{ij}$ at the $\tau$-quantile of $y_{ij}$, and the regression coefficients $\mathit{\delta }\left(\tau \right)$ describe the heterogeneity of ${\tilde{\mathit{z}}}_{ij}$ within the subgroup, the threshold coefficients $\mathit{\gamma }\left(\tau \right)$ determine the subgroup based on whether the linear combination ${\mathit{x}}_i^{\mathrm{T}}\mathit{\gamma }\left(\tau \right)$ exceeds 0. We consider time-invariant biomarkers to prevent time-varying biomarkers from dividing repeated observations of the same sample into different subgroups due to the influence of treatment. In addition, the error terms $\left({\epsilon }_{i\mathrm{1}},\cdots ,{\epsilon }_{im}\right)$ and $\left({\epsilon }_{k\mathrm{1}},\cdots ,{\epsilon }_{km}\right)$ are typically assumed to be independent for any $i\ne k$ stemming from the index set $\left\{\mathrm{1},\cdots ,n\right\}$, but correlated within $\left({\epsilon }_{i\mathrm{1}},\cdots ,{\epsilon }_{im}\right)$ for any fixed $i\in \left\{\mathrm{1},\cdots ,n\right\}$. Assume $P\left({\epsilon }_{ij}<\mathrm{0}|{\mathit{z}}_{ij},{\mathit{x}}_i\right)=\tau .$

For the recognizability of $\mathit{\gamma }\left(\tau \right)$, reference Zhang et al^[15] and Horowitz^[23], assuming that at least one threshold variable has a nonzero coefficient, and the corresponding probability distribution is absolutely continuous with respect to the Lebesgue measure conditional on the other threshold variables. The variables in ${\mathit{x}}_i$ can be rearranged so that $x_{\mathrm{1}i}$ meets this condition, then the model can be rewritten as:

$Q_{y_{ij}}\left(\tau |{\mathit{z}}_{ij},{\mathit{x}}_i\right)={\mathit{z}}_{ij}^{\mathrm{T}}\mathit{\beta }\left(\tau \right)+{\tilde{\mathit{z}}}_{ij}^{\mathrm{T}}\mathit{\delta }\left(\tau \right)I\left(x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}\mathit{\psi }\left(\tau \right)>\mathrm{0}\right)$

where ${\mathit{x}}_i={\left(x_{\mathrm{1}i},{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}\right)}^{\mathrm{T}}$,$\mathit{\psi }\left(\tau \right)=\frac{{\mathit{\gamma }}_{-\mathrm{1}}\left(\tau \right)}{{\gamma }_{\mathrm{1}}\left(\tau \right)}$ with ${\mathit{\gamma }}_{\mathrm{1}}\left(\tau \right)$ and ${\mathit{\gamma }}_{-\mathrm{1}}\left(\tau \right)$ corresponding to the coefficients of $x_{\mathrm{1}i}$ and ${\mathit{x}}_{\mathrm{2}i}$, respectively. If ${\gamma }_{\mathrm{1}}\left(\tau \right)$ is positive, the normalization without altering the coefficients of $\mathit{\beta }\left(\tau \right)$ and $\mathit{\delta }\left(\tau \right)$. If ${\gamma }_{\mathrm{1}}\left(\tau \right)$ is negative, we should redefined $\mathit{\beta }\left(\tau \right)$ and $\mathit{\delta }\left(\tau \right)$ to ensure the consistency of the model^[15]. Moreover, for the convenience of expression, we remove $\left(\tau \right)$ in the parameters below.

In order to estimate all parameters $\mathit{\theta }={\left({\mathit{\beta }}^{\mathrm{T}},{\mathit{\delta }}^{\mathrm{T}},{\mathit{\psi }}^{\mathrm{T}}\right)}^{\mathrm{T}}$, a computationally efficient method is to ignore the possible correlation between repeated observations within the longitudinal data, that is, to minimize the following objective function under the independence framework:

${\mathcal{l}}_n\left(\theta \right)=\frac{\mathrm{1}}{nm}{\displaystyle \sum _{i=\mathrm{1}}^n}{\displaystyle \sum _{j=\mathrm{1}}^m}{\rho }_{\tau }\{y_{ij}-{\mathit{z}}_{ij}^{\mathrm{T}}\mathit{\beta }-{\tilde{\mathit{z}}}_{ij}^{\mathrm{T}}\mathit{\delta }I\left(x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}\mathit{\psi }>\mathrm{0}\right)\}$

where ${\rho }_{\tau }\left(u\right)=\left\{\tau -I\left(u<\mathrm{0}\right)\right\}u$ is the quantile loss function. In order to address the challenge of discontinuity in the objective function, similar to Seo and Linton^[24], we employ an integrated kernel function to smooth this characteristic function, for example, the cumulative distribution function of the standard normal distribution $\Phi \left(\cdot \right)$, then we can derive the smoothing objective function as following:

$S_n\left(\mathit{\theta };h_n\right)=\frac{\mathrm{1}}{nm}{\displaystyle \sum _{i=\mathrm{1}}^n}{\displaystyle \sum _{j=\mathrm{1}}^m}{\rho }_{\tau }\left\{y_{ij}-{\mathit{z}}_{ij}^{\mathrm{T}}\mathit{\beta }-{\tilde{\mathit{z}}}_{ij}^{\mathrm{T}}\mathit{\delta }\Phi \left(\frac{x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}\mathit{\psi }}{h_n}\right)\right\}$

In refer to Zhang et al^[15], this article selects the bandwidth $h_n=\mathrm{0.5}\mathrm{l}\mathrm{o}\mathrm{g}(N)/\sqrt[]{N}$.

Furthermore, it is hoped to select important covariates ${\mathit{z}}_{ij}$ and threshold variables ${\mathit{x}}_{\mathrm{2}i}$. In practice, a few exposure variables of interest are included in ${\tilde{\mathit{z}}}_{ij}$, such as treatment group indices or drug doses, thus variable selection for interaction terms is not considered. Consequently, the objective function, which includes corresponding penalty terms, is formulated as:

$Q_n\left(\mathit{\theta };h_n\right)=S_n\left(\mathit{\theta };h_n\right)+{\displaystyle \sum _{l=\mathrm{1}}^p}{p}_{{\lambda }_{\mathrm{1}}}\left(\left|{\beta }_l\right|\right)+{\displaystyle \sum _{l=\mathrm{1}}^d}{p}_{{\lambda }_{\mathrm{2}}}\left(\left|{\psi }_l\right|\right)$

where $p_{\lambda }\left(\cdot \right)$ is the penalty function, and ${\lambda }_{\mathrm{1}}$ and $\lambda {}_{\mathrm{2}}{}^{}$ are the shrinkage adjustment parameters that control the partial linear regression parameters and threshold parameters, respectively. We consider the following non-convex SCAD penalty function:

$\begin{array}{l}{p}_{\lambda }\left(\left|\alpha \right|\right)=\lambda \left|\alpha \right|I\left(\mathrm{0}<\left|\alpha \right|<\lambda \right)\\ +\frac{(a^{\mathrm{2}}-\mathrm{1}){\lambda }^{\mathrm{2}}-{(\left|\alpha \right|-a\lambda )}^{\mathrm{2}}}{\mathrm{2}(a-\mathrm{1})}I\left(\lambda <\left|\alpha \right|<a\lambda \right)\\ +\frac{(a^{\mathrm{2}}+\mathrm{1}){\lambda }^{\mathrm{2}}}{\mathrm{2}}I\left(\left|\alpha \right|\ge a\lambda \right),\end{array}$

when $\lambda >\mathrm{0}$, and $a=\mathrm{3.7}$. SCAD penalty approaches the optimal Bayesian risk^[20], therefore, the nuisance parameters involved in the model only are ${\lambda }_{\mathrm{1}}>\mathrm{0}$ and ${\lambda }_{\mathrm{2}}>\mathrm{0}$.

Since the objective function $Q_n\left(\mathit{\theta };h_n\right)$ is non-convex and not differentiable everywhere, we use the local linear approximation algorithm of Zou and Li^[25] to transform the non-convex SCAD penalty into a convex approximation form. Specifically,

$p_{{\lambda }_{\mathrm{1}}}\left(\left|{\beta }_l\right|\right)\approx p_{{\lambda }_{\mathrm{1}}}\left(\left|{\widehat{\beta }}_l^{(\mathrm{0})}\right|\right)+p_{{\lambda }_{\mathrm{1}}}^{\text{'}}\left(\left|{\widehat{\beta }}_l^{(\mathrm{0})}\right|\right)\left(\left|{\beta }_l\right|-\left|{\widehat{\beta }}_l^{(\mathrm{0})}\right|\right),\mathrm{f}\mathrm{o}\mathrm{r}{\beta }_l\approx {\widehat{\beta }}_l^{(\mathrm{0})}$

$p_{{\lambda }_{\mathrm{2}}}\left(\left|{\psi }_l\right|\right)\approx p_{{\lambda }_{\mathrm{2}}}\left(\left|{\widehat{\psi }}_l^{(\mathrm{0})}\right|\right)+p_{{\lambda }_{\mathrm{2}}}^{\text{'}}\left(\left|{\widehat{\psi }}_l^{(\mathrm{0})}\right|\right)\left(\left|{\psi }_l\right|-\left|{\widehat{\psi }}_l^{(\mathrm{0})}\right|\right),\mathrm{f}\mathrm{o}\mathrm{r}{\psi }_l\approx {\widehat{\psi }}_l^{(\mathrm{0})}$

Moreover, even if we employ a profile estimation similar to Zhang et al^[15] to estimate the slope ${\left({\mathit{\beta }}^{\mathrm{T}},{\mathit{\delta }}^{\mathrm{T}}\right)}^{\mathrm{T}}$ and threshold parameters $\mathit{\psi }$ separately, the objective function is still non-convex with respect to $\mathit{\psi }$ with fixing ${\left({\mathit{\beta }}^{\mathrm{T}},{\mathit{\delta }}^{\mathrm{T}}\right)}^{\mathrm{T}}$. Owing to the slow computation speed of genetic algorithms and the sensitivity of the Nelder-Mead method to initial values, therefore, we employ the first-order Taylor linear approximation to $\Phi$,

$\Phi \left(\frac{x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}\mathit{\psi }}{h_n}\right)\approx \Phi \left(\frac{x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}{\widehat{\mathit{\psi }}}^{(\mathrm{0})}}{h_n}\right)+{\Phi }^{\text{'}}\left(\frac{x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}{\widehat{\mathit{\psi }}}^{(\mathrm{0})}}{h_n}\right)\frac{{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}}{h_n}(\mathit{\psi }-{\widehat{\mathit{\psi }}}^{(\mathrm{0})}),$

where ${\Phi }^{\text{'}}\left(\cdot \right)$ indicates the probability density function of the standard normal distribution. Accordingly, with initial estimation for $\widehat{\mathit{\delta }}$ and $\left(\widehat{\mathit{\beta }},\widehat{\mathit{\psi }}\right)$, we can derive a new objective function $\tilde{Q}\left(\mathit{\beta },\mathit{\psi }\right)$ with respect to $\mathit{\beta }$ and $\mathit{\psi }$:

$\tilde{Q}\left(\mathit{\beta },\mathit{\psi }\right)={\displaystyle \sum _{i=\mathrm{1}}^n}{\displaystyle \sum _{j=\mathrm{1}}^m}{\rho }_{\tau }\left({\tilde{y}}_{ij}-{\mathit{z}}_{ij}^{\mathrm{T}}\mathit{\beta }-{\tilde{\mathit{x}}}_i^{\mathrm{T}}\mathit{\psi }\right)$

$+N{\displaystyle \sum _{l=\mathrm{1}}^p}{p}_{{\lambda }_{\mathrm{1}}}^{\text{'}}(\left|{\widehat{\beta }}_l\right|)\left|{\beta }_l\right|+N{\displaystyle \sum _{l=\mathrm{1}}^d}{p}_{{\lambda }_{\mathrm{2}}}^{\text{'}}(\left|{\widehat{\psi }}_l\right|)\left|{\psi }_l\right|$

and

${\tilde{y}}_{ij}=y_{ij}-{\tilde{\mathit{z}}}_{ij}^{\mathrm{T}}\widehat{\mathit{\delta }}\Phi \left(\frac{x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}\widehat{\mathit{\psi }}}{h_n}\right)+{\tilde{\mathit{z}}}_{ij}^{\mathrm{T}}\widehat{\mathit{\delta }}{\Phi }^{\text{'}}\left(\frac{x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}\widehat{\mathit{\psi }}}{h_n}\right)\frac{{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}}{h_n}\widehat{\mathit{\psi }},$

${\tilde{\mathit{x}}}_i={\tilde{\mathit{z}}}_{ij}^{\mathrm{T}}\widehat{\mathit{\delta }}{\Phi }^{\text{'}}\left(\frac{x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}\widehat{\mathit{\psi }}}{h_n}\right)\frac{{\mathit{x}}_{\mathrm{2}i}}{h_n}.$

Now, the estimation of $\mathit{\beta }$ and $\mathit{\psi }$ is essentially a penalized linear quantile regression problem, specific iterative algorithm is as follows:

Step 0   Initialize $\left({\widehat{\mathit{\beta }}}^{(\mathrm{0})},{\widehat{\mathit{\psi }}}^{(\mathrm{0})}\right)$.

Step 1   Given $\left({\widehat{\mathit{\beta }}}^{(k-\mathrm{1})},{\widehat{\mathit{\psi }}}^{(k-\mathrm{1})}\right)$, we have:

${\widehat{\mathit{\delta }}}^{(k)}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}{\displaystyle \sum _{i=\mathrm{1}}^n}{\displaystyle \sum _{j=\mathrm{1}}^m}{\rho }_{\tau }\{y_{ij}-{\mathit{z}}_{ij}^{\mathrm{T}}{\widehat{\mathit{\beta }}}^{(k-\mathrm{1})}-{\tilde{\mathit{z}}}_{ij}^{\mathrm{T}}\mathit{\delta }\Phi \left(\frac{x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}{\widehat{\mathit{\psi }}}^{(k-\mathrm{1})}}{h_n}\right)\}$

Step 2   Given ${\widehat{\mathit{\delta }}}^{(k)}$, the estimation $\left({\widehat{\mathit{\beta }}}^{(k)},{\widehat{\mathit{\psi }}}^{(k)}\right)$ is derived by minimizing the following objective function:

$\begin{array}{l}\tilde{Q}\left(\mathit{\beta },\mathit{\psi }\right)={\displaystyle \sum _{i=\mathrm{1}}^n}{\displaystyle \sum _{j=\mathrm{1}}^m}{\rho }_{\tau }\left({\tilde{y}}_{ij}-{\mathit{z}}_{ij}^{\mathrm{T}}\mathit{\beta }-{\tilde{\mathit{x}}}_i^{\mathrm{T}}\mathit{\psi }\right)\\ +N{\displaystyle \sum _{l=\mathrm{1}}^p}{p}_{{\lambda }_{\mathrm{1}}}^{\text{'}}(\left|{\widehat{\beta }}_l^{(k-\mathrm{1})}\right|)\left|{\beta }_l\right|\\ +N{\displaystyle \sum _{l=\mathrm{1}}^d}{p}_{{\lambda }_{\mathrm{2}}}^{\text{'}}\left(\left|{\widehat{\psi }}_l^{(k-\mathrm{1})}\right|\right)\left|{\psi }_l\right|,\end{array}$

and

$\begin{array}{l}{\tilde{y}}_{ij}=y_{ij}-{\tilde{\mathit{z}}}_{ij}^{\mathrm{T}}{\widehat{\mathit{\delta }}}^{(k)}\Phi \left(\frac{x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}{\widehat{\mathit{\psi }}}^{(k-\mathrm{1})}}{h_n}\right)\\ +{\tilde{\mathit{z}}}_{ij}^{\mathrm{T}}{\widehat{\mathit{\delta }}}^{(k)}{\Phi }^{\text{'}}\left(\frac{x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}{\widehat{\mathit{\psi }}}^{(k-\mathrm{1})}}{h_n}\right)\frac{{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}}{h_n}{\widehat{\mathit{\psi }}}^{(k-\mathrm{1})},\end{array}$

set ${\tilde{\mathit{x}}}_i={\tilde{\mathit{z}}}_{ij}^{\mathrm{T}}{\widehat{\mathit{\delta }}}^{(k)}{\Phi }^{\text{'}}\left(\frac{x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}{\widehat{\mathit{\psi }}}^{(k-\mathrm{1})}}{h_n}\right)\frac{{\mathit{x}}_{\mathrm{2}i}}{h_n}$, repeat step 1 and step 2 until the L2 norm of adjacent estimation for all parameters are less than the tolerance error.

Finally, according to the facts ${\rho }_{\tau }(cv)=c{\rho }_{\tau }(v)(c>\mathrm{0})$ and $\left|{\beta }_l\right|={\rho }_{\tau }\left({\beta }_l\right)+{\rho }_{\tau }\left(-{\beta }_l\right)$, set $N{p}_{{\lambda }_{\mathrm{1}}}^{\text{'}}(\left|{\widehat{\beta }}_l^{(k-\mathrm{1})}\right|)=c_l$, and $N{p}_{{\lambda }_{\mathrm{2}}}^{\text{'}}(\left|{\psi }_l^{(k-\mathrm{1})}\right|)=d_l$, then, we rewrite the first and second penalty terms as follows:

${\displaystyle \sum _{l=\mathrm{1}}^p}({\rho }_{\tau }(c_l{\beta }_l)+{\rho }_{\tau }(-c_l{\beta }_l)),{\displaystyle \sum _{l=\mathrm{1}}^d}({\rho }_{\tau }(d_l{\psi }_l)+{\rho }_{\tau }(-d_l{\psi }_l)).$

Furthermore, we construct the "unpenalized" linear quantile regression through incorporating additional pseudo observations. The augmented form of the defining variables is as follows:

${\tilde{y}}_l^{+}=\{\begin{array}{l}{\tilde{y}}_l,l=\mathrm{1},\cdots ,N\\ \mathrm{0},l=N+\mathrm{1},\cdots ,N+\mathrm{2}p+\mathrm{2}d\end{array}$

${\tilde{\mathit{z}}}_l^{+}=\{\begin{array}{l}{z}_l,l=\mathrm{1},\cdots ,N,\\ \left(\mathrm{0},\cdots ,\mathrm{0},c_l,\mathrm{0},\cdots ,\mathrm{0}\right),l=N+\mathrm{1},\cdots ,N+p,\\ \left(\mathrm{0},\cdots ,\mathrm{0},-c_l,\mathrm{0},\cdots ,\mathrm{0}\right),l=N+p+\mathrm{1},\cdots ,N+\mathrm{2}p,\\ \left(\mathrm{0},\cdots ,\mathrm{0}\right),l=N+\mathrm{2}p+\mathrm{1},\cdots ,N+\mathrm{2}p+\mathrm{2}d,\end{array}$

${\tilde{\mathit{x}}}_l^{+}=\{\begin{array}{l}{x}_l,l=\mathrm{1},\cdots ,N,\\ \left(\mathrm{0},\cdots ,\mathrm{0}\right),l=N+\mathrm{1},\cdots ,N+\mathrm{2}p,\\ \left(\mathrm{0},\cdots ,\mathrm{0},d_l,\mathrm{0},\cdots ,\mathrm{0}\right),l=N+\mathrm{2}p+\mathrm{1},\cdots ,N+\mathrm{2}p+d,\\ \left(\mathrm{0},\cdots ,\mathrm{0},-d_l,\mathrm{0},\cdots ,\mathrm{0}\right),l=N+\mathrm{2}p+d+\mathrm{1},\cdots ,N+\mathrm{2}p+\mathrm{2}d,\end{array}$

then the objective function in Step 2 degenerates into linear quantile regression:${\displaystyle \sum _{l=\mathrm{1}}^{N+\mathrm{2}p+\mathrm{2}d}}{\rho }_{\tau }\left({\tilde{y}}_l^{+}-{{\tilde{\mathit{z}}}_l^{+}}^{\mathrm{T}}\mathit{\beta }-{{\tilde{\mathit{x}}}_l^{+}}^{\mathrm{T}}\mathit{\psi }\right)$,and proceed with solving it.

In practice, BIC is frequently employed in variable selection problems for tuning parameter selection,

$\begin{array}{l}\mathrm{B}\mathrm{I}\mathrm{C}\left({\lambda }_{\mathrm{1}},{\lambda }_{\mathrm{2}}\right)=\mathrm{l}\mathrm{o}\mathrm{g}\left({\displaystyle \sum _{i=\mathrm{1}}^n}{\displaystyle \sum _{j=\mathrm{1}}^m}{\rho }_{\tau }\left\{y_{ij}-{\mathit{z}}_{ij}^{\mathrm{T}}{\widehat{\mathit{\beta }}}_{{\lambda }_{\mathrm{1}}}-{\tilde{\mathit{z}}}_{ij}^{\mathrm{T}}\mathit{\delta }\Phi \left(\frac{x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}{\widehat{\mathit{\psi }}}_{{\lambda }_{\mathrm{2}}}}{h_n}\right)\right\}\right)\\ +\frac{\mathrm{l}\mathrm{o}\mathrm{g}(N)}{\mathrm{2}N}{d}_{\lambda },\end{array}$ $d_{\lambda }$ represents the total count of non-zero coefficients in both the slope and threshold components, while $N$ denotes the total number of observations. Following Ruppert and Carroll^[22], we employ a two-step grid search algorithm to derive the optimal tuning parameters $\left({\lambda }_{\mathrm{1}},{\lambda }_{\mathrm{2}}\right)$.

2 Simulation

In order to assess the performance of our proposed estimation in finite samples, we consider the following data generation model:

$y_{ij}={\mathit{z}}_{ij}^{\mathrm{T}}{\mathit{\beta }}_{\mathrm{0}}+{\tilde{\mathit{z}}}_{ij}^{\mathrm{T}}{\mathit{\delta }}_{\mathrm{0}}I\left\{x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}{\mathit{\psi }}_{\mathrm{0}}+\mathrm{0.1}{\Phi }^{-\mathrm{1}}\left(U_i\right)>\mathrm{0}\right\}+(\mathrm{1}+\mathrm{0.1}{z}_{\mathrm{2},ij}){\epsilon }_{ij},$

covariates ${\mathit{z}}_{ij}={\left(\mathrm{1},z_{\mathrm{1},ij},\cdots ,z_{p-\mathrm{1},ij}\right)}^{\mathrm{T}}$, and $z_{\mathrm{1},ij}\sim \mathrm{e}\mathrm{x}\mathrm{p}(\mathrm{1})$, $z_{\mathrm{2},ij}\sim N\left(-\mathrm{1,1}\right)$, the remaining covariates $z_{\mathrm{3},ij},\cdots ,z_{p-\mathrm{1},ij}$ are independent and conform to uniform distribution $U(\mathrm{0,1})$; interaction variable ${\tilde{\mathit{z}}}_{ij}={\left(\mathrm{1},z_{\mathrm{1},ij}\right)}^{\mathrm{T}}$, threshold variables $x_{\mathrm{1}i}\sim N\left(\mathrm{0,1}\right)$, component $x_{\mathrm{2}i}^{\mathrm{*}}\sim N\left(\mathrm{1,1}\right)$ within ${\mathit{x}}_{\mathrm{2}i}={\left(\mathrm{1},x_{\mathrm{2}i}^{\mathrm{*}},\cdots ,x_{di}^{\mathrm{*}}\right)}^{\mathrm{T}}$, the remaining covariates $x_{\mathrm{3}i}^{\mathrm{*}},\cdots ,x_{di}^{\mathrm{*}}$ are independent and conform to uniform distribution $U(\mathrm{0,1})$ too; setting ${\mathit{\beta }}_{\mathrm{0}}={\left(\mathrm{1,1},\mathrm{1,0},\cdots ,{\mathrm{0}}_{p-\mathrm{3}}\right)}^{\mathrm{T}}$, ${\mathit{\delta }}_{\mathrm{0}}={\left(\mathrm{1,1}\right)}^{\mathrm{T}}$, and ${\mathit{\psi }}_{\mathrm{0}}={\left(-\mathrm{1,1},\mathrm{0},\cdots ,{\mathrm{0}}_{d-\mathrm{2}}\right)}^{\mathrm{T}}$. However, in our data generation model, the true value of $\left({\beta }_{\mathrm{01}},{\beta }_{\mathrm{03}},{\psi }_{\mathrm{01}}\right)$ varies with the quantile. The random error term ${\epsilon }_{ij}=\xi +e_{ij}$, where $\xi$ ensures $P\left({\epsilon }_{ij}<\mathrm{0}|z_{ij},x_i\right)=\tau$, we take into account two error distributions:

Case 1 The random error vector conforms to multivariate normal distribution, that is, ${\mathit{\epsilon }}_i={\left({\epsilon }_{i\mathrm{1}},\cdots ,{\epsilon }_{im}\right)}^{\mathrm{T}}$$\sim \mathrm{M}\mathrm{V}\mathrm{N}\left(\mathrm{0},\mathit{\Sigma }\right)$, and the covariance matrix $\mathit{\Sigma }$ follows $\mathrm{A}\mathrm{R}$-1 correlation structure with correlation coefficient $\rho =\mathrm{0.3}$, where the partial coefficients of the model are set to

$\left({\beta }_{\mathrm{01}},{\beta }_{\mathrm{03}},{\psi }_{\mathrm{01}}\right)=\left(\mathrm{1}+{\Phi }^{-\mathrm{1}}\left(\tau \right),\mathrm{1}+\mathrm{0.1}{\Phi }^{-\mathrm{1}}\left(\tau \right),\mathrm{0.1}{\Phi }^{-\mathrm{1}}\left(\tau \right)-\mathrm{1}\right).$

Case 2 The random error vector conforms to multivariate $t$ distribution, that is, ${\mathit{\epsilon }}_i={\left({\epsilon }_{i\mathrm{1}},\cdots ,{\epsilon }_{im}\right)}^{\mathrm{T}}$

$\sim \mathrm{M}\mathrm{V}{\mathrm{T}}_{\mathrm{3}}\left(\mathrm{0},\mathit{\Sigma }\right)$, and the covariance matrix shares the same structure, where $\left({\beta }_{\mathrm{01}},{\beta }_{\mathrm{03}},{\psi }_{\mathrm{01}}\right)=(\mathrm{1}+qt\left(\tau ,\mathrm{d}\mathrm{f}=\mathrm{3}\right),\mathrm{1}+\mathrm{0.1}\times qt\left(\tau ,\mathrm{d}\mathrm{f}=\mathrm{3}\right),$$\mathrm{0.1}{\Phi }^{-\mathrm{1}}\left(\tau \right)-\mathrm{1})$, $qt(\tau ,$$\mathrm{d}\mathrm{f}=\mathrm{3})$ is the $\tau$ quantile of the $t$ distribution with degrees of freedom (df) 3.

We will assess the performance of our estimation across varying numbers of parameters $(p=d=\mathrm{10}$ or $p=d=\mathrm{20})$ or different quantiles $\left(\tau =\mathrm{0.25,0.5,0.75}\right)$, each random experiment yields $\left(n=\mathrm{400}\mathrm{o}\mathrm{r}\mathrm{800}\right)$ samples, with each sample being observed $\left(m=\mathrm{4}\right)$ times, and the simulation is repeated 500 times. We assess the model's performance using four criterias: 1) True Negative (TN): This represents the average number of zero parameters that are correctly classified as zero by the model. A higher TN indicates better model performance; 2) False Negative (FN): The average number of non-zero parameters erroneously classify as zero. A lower FN indicates better model performance; 3) Correct (%): This represents the percentage of accurate true models with non-zero parameters that are correctly identified. The closer this percentage is to 100, the more accurate the model estimation effect is; 4) MSE: It represents the mean square error of the parameter estimation, smaller MSE indicates that the estimated value $\widehat{\theta }$ is closer to the true value ${\theta }_{\mathrm{0}}$, and the better estimation effect of the model. By comparing the mean square error O.MSE of the oracle estimator with the mean square error P.MSE of our method, the P.MSE is closer to O.MSE means the better estimation effect of our approach.

Tables 1 and 2 present the variable selection outcomes of the single-index threshold quantile regression when the error terms follow a multivariate normal distribution or a multivariate T(3) distribution, respectively. In an ideal scenario, the component of TN for the linear segment is $p-\mathrm{3}$, and for the threshold segment it is $d-\mathrm{2}$. The oracle TN in scenario $p=d=\mathrm{10}$ and scenario $p=d=\mathrm{20}$ is 15 or 35, respectively. In all settings, our method can compress the true zero coefficients to zero and correctly identify all relevant covariates (FN=0), with P.MSE and the O.MSE is very close, which indicate that the model estimation is very effective.

Tables 3 and 4 show the estimators of the identified nonzero parameters by our method when the error terms follow a multivariate normal distribution or a multivariate T(3) distribution, respectively. The Biases for all nonzero parameters remain minimal across various parameter counts $(p=d=\mathrm{10}$ or $p=d=\mathrm{20})$ or quantiles $\left(\tau =\mathrm{0.25,0.5,0.75}\right)$, and all metrics: Biases, SDs, and RMSEs decrease as the sample size increases.

We further conduct sensitivity comparisons for two bandwidths $r=\mathrm{l}\mathrm{o}\mathrm{g}\left(N\right)/\sqrt[]{N}$ and $r^{\mathrm{a}\mathrm{d}\mathrm{j}}=\mathrm{1}/\sqrt[\mathrm{5}]{N}$ under $p=d=\mathrm{10}$, $\tau =\mathrm{0.5}$, while keeping the remaining settings as the same as before. The RMSEs of the proposed estimators in Fig. 1 show that these two bandwidths do not exert significant effects on the performance of the proposed method

Fig. 1 The sensitivity diagram of parameter estimation under the bandwidth rates $\mathit{r}$ and ${\mathit{r}}^{\mathbf{a}\mathbf{d}\mathbf{j}}$ with $\mathit{p}=\mathit{d}=\mathrm{10}$ and $\mathit{\tau }=\mathrm{0.5}$

3 Real Data Analysis

We apply the proposed method to analyze the PA.3 trial dataset conducted by the Canadian Cancer Group. In this trial, patients with locally advanced or metastatic pancreatic cancer were randomly assigned to the erlotinib plus gemcitabine group or the gemcitabine monotherapy group. In the primary analysis, the survival rate of the erlotinib plus gemcitabine group was significantly improved compared with the gemcitabine monotherapy group^[26]. Shultz et al^[27] found that when each biomarker was evaluated individually, none of them was significantly associated with erlotinib. However, He et al^[9] selected CA19-9 and AXL from 47 biomarkers to define the treatment-sensitive subgroup, and found that the erlotinib plus gemcitabine could prolong patient's survival time compared with gemcitabine monotherapy. Similarly, Wei et al^[10] selected CA19-9 and AXL from 6 biomarkers to define a treatment-sensitive subgroup and found that erlotinib plus gemcitabine improved patient's quality of life (QoL) scores compared with gemcitabine monotherapy. However, these studies only focused on selecting biomarkers for defining subgroups, not covariates. Moreover, none of these studies investigated whether distinct subgroup effects exist in patients across various QoL levels. Therefore, we use our method to select relevant covariates and biomarkers, and assess the differential therapeutic effect of erlotinib in subgroups defined by linear combinations of these biomarkers, as well as the differences in effects at different quantile levels of the outcomes.

A total of 377 patients were included in the analysis, of whom 196 received erlotinib treatment, and their quality of life was assessed using the global QoL score. For each patient in the trial, QoL scores were recorded at the baseline, as well as 4, 8, 12, and 16 weeks following the treatment. We focus on the changes in patient's QoL scores relative to the baseline (cQoL), and missing values have imputed using multiple imputation. The Q-Q plot of the imputed variable $y=\mathrm{c}\mathrm{Q}\mathrm{o}\mathrm{L}$ and the Shapiro-Wilk normality test indicate that the data do not conform to a normal distribution. The calculated skewness is 0.139 3, indicating a right skewed distribution, so that quantile regression is considered to be more suitable than the mean regression. The covariates included in the analysis are: baseline QoL score (bQoL), time point (Time), baseline age (Age), treatment index B (if the patient received erlotinib plus gemcitabine, B = 1; if the patient received gemcitabine alone, B=0), and pain level (Painmeas). The twelve biomarkers included in the analysis include proteins with lower overall survival and baseline expression levels associated with metastatic diseases, specifically: AXL: AXL receptor tyrosine kinase; CA 19-9: Carbohydrate antigen 19-9; IL8: Interleukin-8; CEA: Carcinoembryonic antigen; MUC-1: Mucin 1, cell surface associate; PDGFRA: Platelet-derived growth factor receptor; EGFR: Epidermal growth factor receptor; PDK1: Pyruvate dehydrogenase kinase 1; BMP-2: Bone morphogenetic protein 2; HER2: Erb-b2 receptor tyrosine kinase 2; PF4: Platelet factor 4; GAS6: Growth arrest-specific protein 6. We standardize the 12 biomarkers and impute missing values using the median.

It is worth noting that our model requires the development of significant biomarker$x_{\mathrm{1}}$and assumes its coefficient to be nonzero. With prior knowledge, we should primarily consider results from expert opinions or existing literature. In terms of analyzing the PA.3 dataset, existing literature has selected CA19-9 and AXL as significant biomarkers in mean regression^[10] and Cox regression^[9] analysis based on threshold models for $y$. Hence, we may consider utilizing CA19-9 or AXL as significant threshold variable for defining subgroups, with one designate as $x_{\mathrm{1}}$ and the rest store in $x_{\mathrm{2}}$, we can also compare the differences in subgroup analysis when different biomarker is designated as $x_{\mathrm{1}}$. We construct the following single-index threshold quantile regression to fit the data:

$\begin{array}{l}\mathrm{c}\mathrm{Q}\mathrm{o}{\mathrm{L}}_{ij}={\beta }_{\mathrm{1}}+{\beta }_{\mathrm{2}}\mathrm{b}\mathrm{Q}\mathrm{o}{\mathrm{L}}_i+{\beta }_{\mathrm{3}}\mathrm{T}\mathrm{i}\mathrm{m}{\mathrm{e}}_{ij}+{\beta }_{\mathrm{4}}\mathrm{A}\mathrm{g}{\mathrm{e}}_i+{\beta }_{\mathrm{5}}{B}_i\\ +{\beta }_{\mathrm{6}}\mathrm{P}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{m}\mathrm{e}\mathrm{a}{\mathrm{s}}_i+{\delta }_{\mathrm{1}}I\left(x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}\mathit{\psi }>\mathrm{0}\right)\\ +{\delta }_{\mathrm{2}}{B}_{i}I\left(x_{\mathrm{1}i}+{\mathit{x}}_{\mathrm{2}i}^{\mathrm{T}}\mathit{\psi }>\mathrm{0}\right)+{\epsilon }_{ij},i=\mathrm{1},\cdots ,\mathrm{377},j=\mathrm{1},\cdots ,\mathrm{4}.\end{array}$

Set $y=\mathrm{c}\mathrm{Q}\mathrm{o}\mathrm{L}$,$\mathit{z}=(\mathrm{1},\mathrm{b}\mathrm{Q}\mathrm{o}\mathrm{L},\mathrm{T}\mathrm{i}\mathrm{m}\mathrm{e},\mathrm{A}\mathrm{g}\mathrm{e},B,\mathrm{P}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}\mathrm{s})$, $\tilde{\mathit{z}}=(\mathrm{1},B)$, $x_{\mathrm{1}}=\mathrm{C}\mathrm{A}\mathrm{19}$-$\mathrm{9}\mathrm{o}\mathrm{r}\mathrm{A}\mathrm{X}\mathrm{L}$, $\mathit{x}=(\mathrm{1},\mathrm{C}\mathrm{A}\mathrm{19}$-$\mathrm{9},\mathrm{A}\mathrm{X}\mathrm{L},\mathrm{I}\mathrm{L}\mathrm{8},\mathrm{C}\mathrm{E}\mathrm{A}$,

$\mathrm{M}\mathrm{U}\mathrm{C}$-$\mathrm{1},\mathrm{P}\mathrm{D}\mathrm{G}\mathrm{F}\mathrm{R}\mathrm{A},\mathrm{E}\mathrm{G}\mathrm{F}\mathrm{R},\mathrm{P}\mathrm{D}\mathrm{K}\mathrm{1},\mathrm{I}\mathrm{P}$-$\mathrm{10},\mathrm{H}\mathrm{E}\mathrm{R}\mathrm{2},\mathrm{P}\mathrm{F}\mathrm{4},$$\mathrm{G}\mathrm{A}\mathrm{S}\mathrm{6}),$

${\mathit{x}}_{\mathrm{2}}$ represents the part of $\mathit{x}$ minus $x_{\mathrm{1}}$, assuming that the error term satisfies $P\left({\epsilon }_{ij}<\mathrm{0}|{\mathit{z}}_{ij},{\mathit{x}}_i\right)=\tau$.

In order to conduct subgroup analysis of outcome variables at different levels as comprehensively as possible, we consider 7 quantile levels $\tau =\mathrm{0.2,0.3},\cdots ,\mathrm{0.8}$, and the same two-step grid search method as in the simulation study is used to select the tuning parameters ${\lambda }_{\mathrm{1}}$ and ${\lambda }_{\mathrm{2}}$. Table 5 presents the point estimator, bootstrap standard error (BS-SE), 95% confidence interval (95%-CI), and p-value of each parameter obtained at the above quantile levels with fixing $x_{\mathrm{1}}=\mathrm{A}\mathrm{X}\mathrm{L}$.

From Table 5, it can be seen that in terms of variable selection, our method sets the coefficients for observation time (Time) and treatment index variable (B) to 0 across all quantiles ranging from 0.2 to 0.8, additionally, it sets the coefficients for age (Age) to 0 at quantiles 0.4 to 0.8, while the effect at quantiles 0.2 to 0.3 have not statistically significant. Thus, we conclude that the effects of observation time (Time) and age (Age) on patient's QoL scores are not significant. However, the selection of treatment effect (B) as an unrelated variable may be due to the masking of the interaction effect between treatment and subgroups, resulting in an insignificant effect, which requires further confirmation in stratified analysis. In addition, CA19-9 is selected as an important biomarker at all quantiles, and the coefficients of the remaining 10 biomarkers are compressed to 0. The biomarkers we selected are consistent with the results of He et al^[9] and Wei et al^[10].

In terms of subgroup analysis, within the mid or low percentiles (0.2-0.6), the p-value of the interaction coefficient ${\delta }_{\mathrm{2}}$ between treatment and subgroup is less than 0.05, this demonstrates a significant difference in therapeutic effect between the two subgroups at the 0.05 significance level. However, within the high percentile range (0.7-0.8), the corresponding p-value is met with $p>\mathrm{0.05}$, this demonstrates that the difference in therapeutic effect between the two subgroups is not significant at the 0.05 significance level, the possible reason is that there is no treatment sensitive subgroup, and further exploration is required through stratified analysis.

Next, we define the biomarker positive subgroup and the biomarker negative subgroup according to linear combination of fixing biomarker AXL, intercept and selected biomarker CA19-9, corresponding coefficients $\left({\psi }_{\mathrm{0}},{\psi }_{\mathrm{1}}\right)$ as shown in Table 5. Linear quantile regression has conducted on all covariates separately for the samples in two subgroups, the estimated therapeutic effect and their 95% confidence intervals at different quantiles present in Fig. 2, and the confidence intervals are derived from rank tests. Figure 2 reveals that there is significant heterogeneity in therapeutic effect between the two subgroups. Specifically, at all quantiles, within the biomarker positive subgroup, the therapeutic effect are all significant positive values, and the 95% confidence intervals mostly fall above the zero line, which indicate that erlotinib has a positive impact on improving QoL scores. Therefore, we believe that at fixed $x_{\mathrm{1}}=\mathrm{A}\mathrm{X}\mathrm{L}$, the identified biomarker positive subgroup is the treatment sensitive subgroup, in which patients receiving treatment of erlotinib plus gemcitabine have significantly improved quality of life scores compared to treatment with gemcitabine alone. However, within the biomarker negative subgroup, the 95% confidence interval for therapeutic effect remains zero, which indicates that the results are not significant. Therefore, we believe that in this subgroup, patients do not benefit more from erlotinib plus gemcitabine treatment compared to treatment with gemcitabine alone.

Fig. 2 The 95% confidence intervals of therapeutic effect when fixing ${\mathit{x}}_{\mathrm{1}}=\mathbf{A}\mathbf{X}\mathbf{L}$

Moreover, Table 6 presents the same four estimators at the corresponding quantile level with fixing $x_{\mathrm{1}}=\mathrm{C}\mathrm{A}\mathrm{19}$-9. As can be seen from Table 6, the corresponding main results of variable selection and subgroup analysis are similar to those in Table 5 with fixing $x_{\mathrm{1}}=\mathrm{A}\mathrm{X}\mathrm{L}$.

Similarly, we also define the biomarker positive subgroup and the biomarker negative subgroup according to linear combination of fixing biomarker CA19-9. The corresponding estimated therapeutic effect and their 95% confidence intervals at different quantiles present in Fig. 3. Figure 3 shows that there is still significant heterogeneity in the therapeutic effect between the two subgroups. But unlike in the case of fixing $x_{\mathrm{1}}=\mathrm{A}\mathrm{X}\mathrm{L}$, when $x_{\mathrm{1}}=\mathrm{C}\mathrm{A}\mathrm{19}$-9 is fixed, the treatment sensitive subgroup with significant therapeutic effect is the identified biomarker negative subgroup, patients in this subgroup who receive treatment with erlotinib plus gemcitabine will significantly improve their overall quality of life score, however, within the biomarker positive subgroup, patients did not benefit more from the combination therapy of erlotinib and gemcitabine. This indicates that the 0-1 values of the subgroup index variables have reversed in both scenarios, that is, the interaction coefficient ${\delta }_{\mathrm{2}}$ between treatment and subgroups is negative when $x_{\mathrm{1}}=\mathrm{C}\mathrm{A}\mathrm{19}$-9, and positive when $x_{\mathrm{1}}=\mathrm{A}\mathrm{X}\mathrm{L}$. It also explains why the coefficient ${\beta }_{\mathrm{5}}$ of the treatment term is selected as a significant correlation variable when $x_{\mathrm{1}}=\mathrm{C}\mathrm{A}\mathrm{19}$-9, while it is compressed to 0 when $x_{\mathrm{1}}=\mathrm{A}\mathrm{X}\mathrm{L}$.

Fig. 3 The 95% confidence intervals of therapeutic effect when fixing${\mathit{x}}_{\mathrm{1}}=\mathbf{C}\mathbf{A}\mathrm{19}$-9

In summary, combining the results of subgroup analysis for setting $x_{\mathrm{1}}=\mathrm{A}\mathrm{X}\mathrm{L}$ and $x_{\mathrm{1}}=\mathrm{C}\mathrm{A}\mathrm{19}$-9, we can see that our method can identify treatment sensitive subgroups. In this subgroup, patients receiving treatment of erlotinib plus gemcitabine are more beneficial in improving their QoL score compared with using gemcitabine alone, and erlotinib is effective for patients with different QoL levels.

Finally, when fixing $x_{\mathrm{1}}=\mathrm{A}\mathrm{X}\mathrm{L}$ or $x_{\mathrm{1}}=\mathrm{C}\mathrm{A}\mathrm{19}$-9 separately, we also want to compare whether the defined subgroups are consistent. We draw Fig. 4 to visualize the overlap of the subgroups divided in both cases, when fixing $x_{\mathrm{1}}=\mathrm{A}\mathrm{X}\mathrm{L}$, the treatment sensitive subgroup corresponding to the biomarker positive subgroup is defined as group 1; the treatment insensitive subgroup corresponding to the biomarker negative subgroup is defined as group 2. Conversely, when fixing $x_{\mathrm{1}}=\mathrm{C}\mathrm{A}\mathrm{19}$-9, the treatment sensitive subgroup corresponding to the biomarker negative subgroup is defined as group 1; the treatment insensitive subgroup corresponding to the biomarker positive subgroup is defined as group 2. Figure 4 shows the intersection and union of group 1 and group 2 in two different scenarios, the orange bars corresponding to 11 indicate that they are all assigned to group 1 in both scenarios, the blue bars corresponding to 22 signify that they are all assigned to group 2, 12 and 21 represent the number of individuals assigned to distinct subgroups in each scenario, with no intersection, and are uniformly displayed with gray bars.

Fig. 4 Schematic diagram of subgroup overlap for two different settings

In order to clearly display the proportion of patients divided into different subgroups in both scenarios, we set the highest count point of each subgraph to 100 and compare the height of the orange and gray bars. It can be found that the overlap of the treatment sensitive subgroups identified in both scenarios is very high, which indicates that the application of our method for subgroup identification has good robustness in the selection of $x_{\mathrm{1}}$.

4 Conclusion

This article considers a single-index threshold quantile regression for subgroup analysis, and selects covariates and biomarkers for defining subgroups, an efficient algorithm is proposed to smooth and locally linearly approximate the subgroup index function and non-convex penalty, respectively. Based on pseudo observations, the corresponding estimation problem degenerates into linear quantile regression, and all unknown parameters are iteratively solved through a two-step estimation process. Numerical simulations demonstrate that the proposed algorithm performs well with a moderate number of variables. Analysis of real data from the PA.3 trial shows that it can distinguish patient between treatment sensitive and treatment insensitive subgroups.

However, our model only considers the case of one threshold, so patient can only be divided into two subgroups. Furthermore, multiple thresholds can be introduced to divide subjects into multiple subgroups with different covariate effects, for related studies, please refer to Li et al^[28-29]. In addition, the recognizability constraint in our model requires specifying a threshold variable $x_{\mathrm{1}}$ with non-zero coefficient. In practical applications, this variable can be determined through professional knowledge. However, when prior knowledge is not available, it may be achieved by extending the score type specification test method designed by Zhang et al^[30] for selecting the non-zero threshold variables. In addition, considering the internal correlation of longitudinal data and conducting variable selection and parameter estimation under certain condition $p>n$ is our future research direction.

References

All Tables

All Figures

	Fig. 1 The sensitivity diagram of parameter estimation under the bandwidth rates $\mathit{r}$ and ${\mathit{r}}^{\mathbf{a}\mathbf{d}\mathbf{j}}$ with $\mathit{p}=\mathit{d}=\mathrm{10}$ and $\mathit{\tau }=\mathrm{0.5}$
In the text

	Fig. 2 The 95% confidence intervals of therapeutic effect when fixing ${\mathit{x}}_{\mathrm{1}}=\mathbf{A}\mathbf{X}\mathbf{L}$
In the text

	Fig. 3 The 95% confidence intervals of therapeutic effect when fixing${\mathit{x}}_{\mathrm{1}}=\mathbf{C}\mathbf{A}\mathrm{19}$-9
In the text

	Fig. 4 Schematic diagram of subgroup overlap for two different settings
In the text