Estimation of Optimal Individualized Treatment Rules Using a Covariate-Specific Treatment Effect Curve with High-dimensional Covariates

TL;DR

This paper introduces a semi-parametric method to estimate covariate-specific treatment effects and optimal individualized treatment rules in high-dimensional settings, enabling personalized medicine with confidence quantification.

Contribution

It proposes a flexible, high-dimensional covariate modeling approach with confidence bands for treatment effect estimation and personalized treatment decision-making.

Findings

Accurately estimates covariate-specific treatment effects in high-dimensional data.

Provides simultaneous confidence bands for treatment effect curves.

Enables identification of patient subgroups benefiting from specific treatments.

Abstract

With a large number of baseline covariates, we propose a new semi-parametric modeling strategy for heterogeneous treatment effect estimation and individualized treatment selection, which are two major goals in personalized medicine. We achieve the first goal through estimating a covariate-specific treatment effect (CSTE) curve modeled as an unknown function of a weighted linear combination of all baseline covariates. The weight or the coefficient for each covariate is estimated by fitting a sparse semi-parametric logistic single-index coefficient model. The CSTE curve is estimated by a spline-backfitted kernel procedure, which enables us to further construct a simultaneous confidence band (SCB) for the CSTE curve under a desired confidence level. Based on the SCB, we find the subgroups of patients that benefit from each treatment, so that we can make individualized treatment selection. The innovations of the proposed method are three-fold. First, the proposed method can quantify variability associated with the estimated optimal individualized treatment rule with high-dimensional covariates. Second, the proposed method is very flexible to depict both local and global associations between the treatment and baseline covariates in the presence of high-dimensional covariates, and thus it enjoys flexibility while achieving dimensionality reduction. Third, the SCB achieves the nominal confidence level asymptotically, and it provides a uniform inferential tool in making individualized treatment decisions.

Figures

Tables

Table 1. Table 1: Simulation results for example 1-3. Oracle estimator is obtained when true index sets are given. True functions are x ​ ( 1 − x ) , x 2 , − exp ⁡ ( x ) / 1.5 𝑥 1 𝑥 superscript 𝑥 2 𝑥 1.5 x(1-x),x^{2},-\exp(x)/1.5 . MSE and MAE represent the mean squared error and the mean absolute error of the estimator for g 1 ​ ( ⋅ ) subscript 𝑔 1 ⋅ g_{1}(\cdot) , and CP represents the coverage probability of its SCBs.

		Model 1				Model 2				Model 3
$(n,p)$		MSE	MAE	CP		MSE	MAE	CP		MSE	MAE	CP
(500, Oracle)		0.230	0.364	0.880		0.217	0.336	0.861		0.389	0.400	0.840
(750, Oracle)		0.125	0.276	0.925		0.113	0.251	0.905		0.143	0.268	0.905
(1000, Oracle)		0.095	0.23	0.947		0.088	0.230	0.935		0.083	0.220	0.940
(500, 10)		0.280	0.494	0.874		0.325	0.428	0.889		0.405	0.461	0.858
(750, 10)		0.169	0.318	0.897		0.151	0.293	0.924		0.329	0.472	0.885
(1000, 10)		0.112	0.258	0.939		0.095	0.243	0.929		0.099	0.241	0.931
(500, 50)		0.530	0.461	0.871		0.506	0.576	0.870		0.679	0.621	0.863
(750, 50)		0.258	0.371	0.930		0.331	0.426	0.903		0.280	0.393	0.920
(1000, 50)		0.100	0.248	0.957		0.129	0.280	0.955		0.099	0.242	0.959
(500, 100)		0.698	0.594	0.836		0.525	0.566	0.826		0.640	0.660	0.862
(750, 100)		0.153	0.305	0.906		0.452	0.477	0.896		0.141	0.291	0.907
(1000, 100)		0.108	0.230	0.947		0.188	0.338	0.920		0.093	0.237	0.940
(500, 500)		0.718	0.789	0.826		0.642	0.777	0.835		0.683	0.519	0.831
(750, 500)		0.413	0.503	0.857		0.434	0.453	0.860		0.524	0.444	0.854
(1000, 500)		0.123	0.249	0.930		0.174	0.319	0.925		0.143	0.287	0.912

Table 2. Table 2: Simulation results for example 4-5 (mis-specified cases) with p = 50 𝑝 50 p=50 and n = 1000 𝑛 1000 n=1000 . MSE and MAE represent the mean squared error and the mean absolute error of the estimator for g 1 ​ ( ⋅ ) subscript 𝑔 1 ⋅ g_{1}(\cdot) , and CP represents the coverage probability of its SCBs.

	Model 4				Model 5
	MSE	MAE	CP		MSE	MAE	CP
c = 0	0.100	0.248	0.957		0.100	0.248	0.957
c=0.1	0.138	0.212	0.933		0.161	0.293	0.938
c=0.2	0.136	0.235	0.931		0.266	0.338	0.914
c=0.5	0.208	0.341	0.904		0.290	0.415	0.895

Table 3. Table 3: Comparisons between the pointwise confidence interval and the SCBs for the empirical coverage rate of a group of randomly generated new patients with group size 1, 10, 100, 200, and 500, respectively.

	group size	1	10	100	200	500
Model 1	CI	0.962	0.806	0.621	0.561	0.554
	SCBs	1.000	0.991	0.952	0.948	0.935
Model 2	CI	0.959	0.911	0.835	0.771	0.773
	SCBs	1.000	1.000	0.963	0.951	0.954
Model 3	CI	0.961	0.853	0.694	0.661	0.668
	SCBs	1.000	0.981	0.952	0.934	0.937

Table 4. Table 4: Model selection results for example 1-3. FP(false positive): zero is estimated as non-zero; FN(false negative): non-zero is estimated as zero; C(correctly selected): all relevant covariates are selected; IC(incorrectly selected): some relevant covariates are not selected.

		Model 1					Model 2					Model 3
$(n,p)$		FPR	FNR	C	IC		FPR	FNR	C	IC		FPR	FNR	C	IC
(500, 10)		0.282	0.025	0.885	0.115		0.274	0.019	0.904	0.095		0.381	0.015	0.909	0.090
(750, 10)		0.206	0.008	0.956	0.043		0.170	0.007	0.960	0.039		0.280	0.018	0.941	0.058
(1000, 10)		0.162	0.000	1.000	0.000		0.110	0.005	0.994	0.005		0.217	0.002	0.989	0.010
(500, 50)		0.271	0.042	0.785	0.214		0.172	0.181	0.564	0.435		0.165	0.236	0.535	0.465
(750, 50)		0.108	0.041	0.823	0.176		0.178	0.025	0.892	0.107		0.094	0.125	0.563	0.436
(1000, 50)		0.070	0.009	0.952	0.047		0.111	0.005	0.975	0.025		0.082	0.065	0.788	0.211
(500, 100)		0.063	0.093	0.625	0.375		0.107	0.084	0.605	0.395		0.093	0.087	0.670	0.330
(750, 100)		0.055	0.065	0.703	0.296		0.098	0.044	0.688	0.311		0.040	0.106	0.681	0.318
(1000, 100)		0.031	0.032	0.845	0.154		0.053	0.037	0.830	0.170		0.021	0.061	0.803	0.196
(500, 500)		0.114	0.138	0.471	0.519		0.119	0.164	0.223	0.777		0.104	0.166	0.243	0.757
(750, 500)		0.974	0.096	0.540	0.460		0.069	0.138	0.294	0.706		0.059	0.132	0.302	0.698
(1000, 500)		0.030	0.062	0.609	0.391		0.039	0.102	0.413	0.587		0.049	0.094	0.449	0.551

Table 5. Table 5: Estimates of β 1 subscript 𝛽 1 \beta_{1} and the corresponding standard errors (in parentheses). First row: using the unpenalized model; second row: using the penalized model.

Variable	Gender	Height	Weight	Huozheng	Fengzheng	Xueyu	Tanshi	Qixu	Yuzheng	Xuexu	Yinxu	Yangxu
${\beta }_1$	0.88	0.10	-0.04	-0.06	0.08	-0.09	-0.08	0.04	0.18	-0.24	0.30	-0.03
	(0.23)	(0.01)	(0.01)	(0.03)	(0.02)	(0.01)	(0.03)	(0.11)	(0.04)	(0.05)	(0.07)	(0.02)
Penalizd ${\beta }_1$	0.70	0.11	-0.05	-0.04	0.14	-0.15	0.00	0.00	0.21	-0.38	0.49	-0.14
	(0.31)	(0.02)	(0.01)	(0.02)	(0.03)	(0.01)	(0.00)	(0.00)	(0.05)	(0.07)	(0.09)	(0.06)

Equations

$$\mathrm{CSTE}(X)=\mathrm{logit}(E(Y(1)|X))-\mathrm{logit}(E(Y(0)|X)),$$

$$\mathrm{CSTE}(X)=\mathrm{logit}(E(Y|X,Z=1))-\mathrm{logit}(E(Y|X,Z=0)).$$

$$\mathrm{logit}(\mu (X,Z))=g_1(X^{\top }{\beta }_1)\cdot Z+g_2(X^{\top }{\beta }_2),$$

$$\mathrm{CSTE}(X)=\mathrm{logit}(\mu (X,1))-\mathrm{logit}(\mu (X,0))=g_1(X^{\top }{\beta }_1).$$

$$\frac{1}{n}\sum _{i=1}^n\log \{1+\exp (g_1(X_i^{\top }{\beta }_1)Z_i+g_2(X_i^{\top }{\beta }_2))\}-\frac{1}{n}\sum _{i=1}^n{Y}_i\{g_1(X_i^{\top }{\beta }_1)Z_i+g_2(X_i^{\top }{\beta }_2)\}.$$

$$\frac{1}{n}\sum _{i=1}^n\log \{1+\exp (g_1(X_i^{\top }{\beta }_1)Z_i+g_2(X_i^{\top }{\beta }_2))\}-\frac{1}{n}\sum _{i=1}^n{Y}_i\{g_1(X_i^{\top }{\beta }_1)Z_i+g_2(X_i^{\top }{\beta }_2)\}+\sum _{k=1}^2\sum _{j=2}^{p}p({\beta }_{kj},\lambda ),$$

$$\frac{1}{n}\sum _{i=1}^n\log \{1+\exp (B_1(X_i^{\top }{\beta }_1{)}^{\top }{\delta }_1{Z}_i+B_2(X_i^{\top }{\beta }_2{)}^{\top }{\delta }_2)\}$$

$$-\frac{1}{n}\sum _{i=1}^n{Y}_i\{B_1(X_i^{\top }{\beta }_1{)}^{\top }{\delta }_1{z}_i+B_2(X_i^{\top }{\beta }_2{)}^{\top }{\delta }_2\}+\sum _{k=1}^2\sum _{j=2}^{p}p({\beta }_{kj},\lambda ),$$

$$\mathrm{logit}(\mu(X,Z))\approx B(X^{\top}\beta_{k})^{\top}\delta_{1}Z+B(X^{\top}\beta_{k})^{\top}\delta_{2}=\left(ZB(X^{\top}\beta_{k})^{\top},B(X^{\top}\beta_{k})^{\top}\right)\left(\begin{array}[]{c}\delta_{1}\\ \delta_{2}\end{array}\right).$$

$${\tilde{g}}_k(X^{\top }{\beta }_k)\approx {\tilde{g}}_k(X^{\top }{\beta }_k^{old})+{\tilde{g}}_k^{\prime }(X^{\top }{\beta }_k^{old})X^{\top }J({\beta }_k^{old})({\beta }_{k,-1}-{\beta }_{k,-1}^{old}),$$

$$Z{\tilde{g}}_1(X^{\top }{\beta }_1^{old})+{\tilde{g}}_2(X^{\top }{\beta }_2^{old})-Z{\tilde{g}}_1^{\prime }(X^{\top }{\beta }_1^{old})X^{\top }J({\beta }_1^{old}){\beta }_{1,-1}^{old}-{\tilde{g}}_2^{\prime }(X^{\top }{\beta }_2^{old})X^{\top }J({\beta }_2^{old}){\beta }_{2,-1}^{old}.$$

$$l(a,b,X)=-\frac{1}{n}\sum _i[Y_i(a{Z}_i+{\hat{g}}_2(X_i^{\top }{\hat{\beta }}_2))-\log (1+\exp (a{Z}_i+{\hat{g}}_2(X_i^{\top }{\hat{\beta }}_2)))]K_h(X_i^{\top }{\hat{\beta }}_1-X^{\top }{\hat{\beta }}_1),$$

$$\mathrm{CSTE}(X_0)\approx {\hat{g}}_{1,sbk}(X_0^{\top }{\hat{\beta }}_1)=\hat{a}(X_0).$$

$$E(Y|X,Z)=\mu (X,Z)=S(g_1(X^{\top }{\beta }_1)Z+g_2(X^{\top }{\beta }_2)).$$

$$\mathcal{M}=\left\{g(x,z)=g_1(x^{\top }{\beta }_1^0)z+g_2(x^{\top }{\beta }_2^0),E\{g_k(X^{\top }{\beta }_k^0){\}}^2<\infty \right\}.$$

$${\mathbb{P}}_{\mathcal{M}}(U)=\arg \underset{g\in \mathcal{M}}{\min }E\left[w^0\{U-g(X,Z){\}}^2\right],$$

$${\mathbb{P}}_{\mathcal{M}}(U)=\{{\mathbb{P}}_{\mathcal{M}}(U_1),...,{\mathbb{P}}_{\mathcal{M}}(U_d){\}}^{\top }.$$

$$\Omega_{n}=\frac{1}{n}\sum_{i=1}^{n}J(\beta_{(1)}^{0})^{\top}\left(\begin{array}[]{c}g_{1}^{\prime}(X_{i,11}^{\top}\beta_{11}^{0})\widetilde{X}_{i,11}Z_{i}\\ g_{2}^{\prime}(X_{i,21}^{\top}\beta_{21}^{0})\widetilde{X}_{i,12}\end{array}\right)\left(\begin{array}[]{c}g_{1}^{\prime}(X_{i,11}^{\top}\beta_{11}^{0})\widetilde{X}_{i,11}Z_{i}\\ g_{2}^{\prime}(X_{i,21}^{\top}\beta_{21}^{0})\widetilde{X}_{i,12}\end{array}\right)^{\top}J(\beta_{(1)}^{0}),$$

$$\Phi_{n}=\frac{1}{n}\sum_{i=1}^{n}\sigma^{2}(X_{i},Z_{i})J(\beta_{(1)}^{0})^{\top}\left(\begin{array}[]{c}g_{1}^{\prime}(X_{i,11}^{\top}\beta_{11}^{0})\widetilde{X}_{i,11}Z_{i}\\ g_{2}^{\prime}(X_{i,21}^{\top}\beta_{21}^{0})\widetilde{X}_{i,21}\end{array}\right)\left(\begin{array}[]{c}g_{1}^{\prime}(X_{i,11}^{\top}\beta_{11}^{0})\widetilde{X}_{i,11}Z_{i}\\ g_{2}^{\prime}(X_{i,21}^{\top}\beta_{21}^{0})\widetilde{X}_{i,21}\end{array}\right)^{\top}J(\beta_{(1)}^{0}),$$

$${\sigma }^2(u)=E[\mu (u,Z)(1-\mu (u,Z))|Z=1],D(u)=f(u){\sigma }^2(u),v^2(u)=||K|{|}_2^{2}f(u){\sigma }^2(u),$$

$$Q_h(\alpha )=a_h+a_h^{-1}[\log (\sqrt{C_K}/(2\pi ))-\log (-\log \sqrt{1-\alpha })]$$

$$\underset{n\to \infty }{\lim }P\left(\underset{x\in \mathcal{C}}{\sup }\left|\frac{{\hat{g}}_{1,sbk}(x^{\top }{\hat{\beta }}_1)-g_1(x^{\top }{\beta }_1^0)}{{\sigma }_n(x^{\top }{\hat{\beta }}_1)}\right|\le Q_h(\alpha )\right)=1-\alpha ,$$

Full text

Keywords: Personalized medicine; Semiparametric model; High-dimensional covariates; Optimal treatment selection.

1 Introduction

Personalized medicine aims to tailor medical treatments according to patient characteristics, and it has gained much attention in modern biomedical research. The success of personalized medicine crucially depends on the development of reliable statistical tools for estimating an optimal treatment regime given the data collected from clinical trials or observational studies (KL19). In the literature, there are two general statistical approaches and their hybrids for deriving an optimal Individualized Treatment Rule (ITR) based on clinical trial data. The first general approach targets at direct optimization of the population average outcome under an ITR. Several independent research groups have proposed a new framework of deriving the ITR that recasts the problem of maximizing treatment benefit as a weighted classification problem (Zhang2012a; Rubin2012; Zhao2012; ZhouKosorork2015; ChenKosork2016). Several semi-parametric and non-parametric methods have been further proposed under this framework (Zhao2012; Huang2014; Huang2015; Laber2015; Zhu2018), and they are more robust against model misspecification than the parametric approaches. Under this framework, the convergence rates of the resulting estimators were thoroughly studied in Zhao2012 and zhang2018, among others, but it is generally difficult to build statistical inference upon the estimated ITRs (laberqian2019). To solve this problem, jslz2019 proposed an entropy learning method, and wager2018 established asymptotic normality of a causal forest estimator.

The second general approach is to model the difference in average outcomes between two treatment groups conditional on predictive biomarkers that provide information about the effect of a therapeutic intervention. There are hybrids from the two approaches using mean-modeling (Zhang2012a; taylor2015; zhang2018c; luckett), but the hybrid methods also inherit the difficulty of developing inferential tools from the first general approach. In this paper, we focus on the idea of the second approach. With a single predictive biomarker, several authors have proposed to plot the estimated conditional treatment difference against the biomarker’s values or its percentiles (e.g., a covariate-specific treatment effect (CSTE) curve) obtained from nonparametric smoothing, along with its confidence bands, to derive an optimal ITR (zhou2012bate; ma2014treatment; Janes2014; CSTE). It provides a direct and effective visual way of using the predictive biomarker for deriving an ITR, but the nonparametric smoothing method suffers from “curse of dimensionality”. With multiple biomarkers, Qian2011 considered a parametric conditional mean model that involves higher order terms, and employed a penalized method for estimation. doi:10.1093/biostatistics/kxq060 proposed a two-step method, and Kang2014 devised a boosting iterative algorithm to reduce the bias caused by the model misspecification of a working model. Moreover, shi2018 proposed a penalized multi-stage A-learning, and zhu2019 developed a high-dimensional Q-learning method to simultaneously estimate the optimal dynamic treatment regimes and select important variables. The existing methods try to either provide a flexible modeling strategy by relaxing the restrictive structural assumption embodied in parametric models, or handle high-dimensional covariates, but not both.

In our paper, we aim at developing a new method that estimates the outcome model in a flexible way as well as selecting important variables simultaneously, when a large number of baseline characteristics such as genetic variables are present. Moreover, we would like to provide a uniform inferential procedure that takes into account uncertainty associated with the estimated optimal ITR, and graphically explores heterogeneity of treatment effects based on varying levels of biomarkers. To achieve these goals, we propose a generalization of the CSTE curve suitable for high dimensional baseline covariates. The CSTE curve represents the predictive ability of covariates in evaluating whether a patient responds better to one treatment over another. It depicts the heterogeneous treatment effects on the outcome through an unknown function of covariates, and thus it can visually represent the magnitude of the predictive ability of the covariates.

We propose a penalized semi-parametric modeling approach for estimating the CSTE curve and selecting variables simultaneously. The proposed semi-parametric model enjoys flexibility while achieving dimensionality reduction. It is motivated by the logistic varying-coefficient model considered in CSTE for treatment selection with one covariate, and it meets the immediate needs from modern biomedical studies which can have a large number of baseline covariates. To make use of all covariates, one simple but effective way is to derive a weighted linear combination of all covariates as a summary predictor, and model the CSTE curve as an unknown function of this summary predictor. The weight or the coefficient for each covariate represents how important the covariate is for the prediction of the outcome. As a result, our proposed model involves two sets of high-dimensional coefficients fed into additive unknown functions for the two treatment groups. The development of the estimation procedure and the associated statistical properties is challenging and it needs different tools from the high-dimensional single-index model (RADCHENKO2015266). Moreover, we establish different convergence rates for the estimators of the high-dimensional coefficients and the estimators for the unknown additive functions, respectively. This new theoretical result makes it possible to further construct a simultaneous confidence band (SCB) for the CSTE curve based on the asymptotic extreme value distribution of a spline-backfitted kernel estimator (MR2382655; doi:10.1080/01621459.2013.763726; Zheng2016) for the unknown function of interest, while RADCHENKO2015266 and other related works only provide a nonseparable convergence rate for the estimators of the coefficients and the unknown function in single-index models. Based on the SCB, we identify the subgroups of patients that benefit from each treatment, and the proposed method is flexible enough to depict both local and global associations between the treatment and baseline covariates.

The rest of the paper is organized as follows. In Section 2, we introduce the CSTE curve and the proposed semi-parametric logistic single-index coefficient model. Section 3 presents the estimation procedure and the asymptotic properties of the proposed estimators. Section 4 illustrates the application of the CSTE curves and the SCBs for treatment selection. In Section 5 we evaluate the finite sample properties of the proposed method via simulation studies, while Section 6 illustrates the usefulness of the proposed method through the analysis of a real data example. A discussion is given in Section 7. All technical proofs are relegated to the online Supplemental Materials.

2 Methodology

We consider a sample of $n$ subjects, a binary treatment, denoted by $Z_{i}=1$ if the subject $i$ is assigned to treatment and $Z_{i}=0$ otherwise, a $p$-dimensional vector of covariates, denoted by $X_{i}$, and binary-valued outcomes, denoted by $Y_{i}$. Let $\left(Y_{i},Z_{i},X_{i}^{\top}\right),i=1,...,n$, be independent and identically distributed (i.i.d.) copies of $\left(Y,Z,X\right)$. The goal is to estimate the optimal treatment regime using the observed data. We consider the CSTE curve given in CSTE, which has the following form:

[TABLE]

where $\mathrm{logit}(u)=\log(u)-\log(1-u)$, and $Y(1)$ and $Y(0)$ denote the potential outcomes if the active treatment and the control treatment are received, respectively (rubin2005). Moreover, $Y\equiv ZY(1)+(1-Z)Y(0)$. Under the unconfoundeness assumption such that $(Y(0),Y(1))\bot Z|X$, the CSTE curve can be re-expressed as

[TABLE]

Denoting $\mu(X,Z)=E(Y|X,Z)$, we model the logarithm of odds ratio as

[TABLE]

where $g_{1}(\cdot)$ and $g_{2}(\cdot)$ are unknown single-valued functions of $p$ variables, $\beta_{1}=(\beta_{11},\dots,\beta_{1p})$ and $\beta_{2}=(\beta_{21},\dots,\beta_{2p})$ are two $p$-vectors of unknown parameters. We see

[TABLE]

The two sets of high-dimensional coefficients and the two unknown functions in (1) are chosen to simultaneously maximize the log-likelihood function of the binomial distribution. The proposed model is a flexible semi-parametric model and is robust against model misspecification. We call it sparse logistic single index coefficient model (SLSICM). Our SLSICM contains the varying coefficient (VC) model considered in CSTE as a special case. We have the same form when $p=1$. As an extension of the VC model, the varying-index coefficient model considered in VICM can be applied to the cases with several covariates and continuous responses. Moreover, an important and related work, song2017, proposed a semi-parametric model in which $g_{1}$ has a single-index structure and $g_{2}$ is a pure nonparametric function of $X$. This model suffers from the curse of dimensionality when the number of covariates is large. In our paper, we allow the number of covariates to be much larger than sample size. As a result, our proposed SLSICM involves two sets of high-dimensional unknown coefficients built into additive unknown functions with a nonlinear link function, and it is considered as a hybrid of the high-dimensional single index model (HSIM) and the nonparametric additive model. Developing the estimation procedure and the statistical theories for SLSICM is challenging, and it needs different tools from HSIM. It is worth noting that estimation of HSIM has been studied in the past several years, and most existing works use a sliced inverse regression approach to estimate the index coefficients in HSIM under a linearity condition on the covariates; see JL2014 and NLC2016, among others. This method is not applicable to our setting, and it does not directly estimate the unknown function. RADCHENKO2015266 proposed to estimate the coefficients and the unknown function in HSIM jointly, but they only provide a nonseparable convergence rate for the resulting estimators of the coefficients and the unknown function. In our proposed SLSICM, we establish different convergence rates for the estimators of the coefficients $\beta_{1}$ and $\beta_{2}$ and the estimators of the unknown functions $g_{1}$ and $g_{2}$, respectively. This new theoretical result makes it possible to further construct the asymptotic SCB based on the asymptotic extreme value distribution of a spline-backfitted kernel estimator for the CSTE curve.

Unlike parametric models, the parameter vectors $\beta_{k}$ are not identified without further assumptions. For the purpose of model identification, we assume that $\beta_{k}$ for $k=1,2$ belong to the parameter space $\Theta=\{\beta=(\beta_{1}^{\top},\beta_{2}^{\top})^{\top}:\|\beta_{k}\|=1,\beta_{k1}>0,\beta_{k}\in R^{p},k=1,2\},$ where $\|\cdot\|$ denotes the $l_{2}$ norm of a vector. This restriction together with Assumption 2 given in Section 3.2 will make the model (1) identifiable. We will provide a formal proof of the model identifiability in Section 1.1 of the Supplemental Materials. Based on the constraint that $\|\beta_{k}\|=1$, we eliminate the first component in $\beta_{k}$ and obtain the resulting parameter space: $\Theta_{-1}=\{\beta_{k,-1}=(\beta_{k2},\dots,\beta_{kp})^{\top}:\sum\nolimits_{j>1}\beta_{kj}^{2}<1,k=1,2\}$. Let $\beta_{k1}=\sqrt{1-\sum\nolimits_{j>1}\beta_{kj}^{2}}$. The derivative with respect to the coefficients $(\beta_{k2},\dots,\beta_{kp})^{\top}$ can be easily obtained using the chain rule under the above parameter space. For high-dimensional problems (with a large number of covariates), $p$ can be much larger than $n$ but only a small number of covariates are important or relevant for treatment selection. To this end, we assume that the number of nonzero elements increases as $n$ increases, but it is much smaller than $n$. Without loss of generality, we assume that only the first $s_{k}=s_{kn}$ components of $\beta_{k}$ are non-zeros, i.e., we can write the true values as $\beta_{k}$ $=(\beta_{k1},\dots,\beta_{ks_{k}},0,\dots,0)^{\top}$. For model identifiability, we require one component in $\beta_{k}$ to be nonzero. Without loss of generality, we let the first component $\beta_{k1}$ to be nonzero.

3 Estimation and Theory

3.1 Algorithm

In this subsection, we discuss the estimation of model (1). We minimize the negative log-likelihood function simultaneously with respect to the parameters $\beta_{k}$ and the functions $g_{k}$ for $(k=1,2)$. The SLSICM has the form $\mathrm{logit}(\mu(X,Z))=g_{1}(X^{\top}\beta_{1})\cdot Z+g_{2}(X^{\top}\beta_{2}).$ Therefore, we seek the minimizer of the following negative log-likelihood function given $(X_{i},Y_{i},Z_{i}),i=1,\dots,n$:

[TABLE]

To overcome the problem of high-dimensional covariates, we exploit the sparsity through parameter regularization. With the sparsity constraint of $\beta_{k}$’s, we minimize the following penalized negative log-likelihood

[TABLE]

where $p(\cdot)$ is a penalty function with a tuning parameter $\lambda$ that controls the level of sparsity in $\beta_{k}$, $k=1,2$.

The functions $g_{k}(\cdot)$, $k=1,2$, are unspecified and are estimated using B-splines regression. Next, we introduce the B-splines that will be used to approximate the unknown functions. For $k=1,2$, we assume the support of $g_{k}(\cdot)$ is $[\inf_{X}(X^{\top}\beta_{k}),\sup_{X}(X^{\top}\beta_{k})]=[a_{k},b_{k}].$ Let $a_{k}=t_{0,0}<t_{1,k}<\cdots<t_{N_{k},k}<b_{k}=t_{N_{k}+1}$ be an equally-spaced partition of $[a_{k},b_{k}]$, called interior knots. Then, $[a_{k},b_{k}]$ is divided into subintervals $I_{\ell,k}=[t_{\ell,k},t_{\ell+1,k})$, $0\leq\ell\leq N_{k}-1$ and $I_{N_{k}}=[t_{N_{k}},t_{N_{k}+1}]$, satisfying $\max_{0\leq\ell\leq N_{k}}|t_{\ell+1,k}-t_{\ell,l}|/\min_{0\leq\ell\leq N_{k}}|t_{\ell+1,k}-t_{\ell,k}|\leq M$ uniformly in $n$ for some constant $0<M<\infty$, where $N_{k}\equiv N_{k,n}$ increases with the sample size $n$. We write the normalized B spline basis of this space (MR1900298) as $B_{k}(u_{k})=\{B_{\ell,k}(u_{k}):1\leq\ell\leq L_{n,k}\}^{\top}$, where the number of spline basis functions is $L_{n,k}=L_{k}=N_{k}+q_{k}$, and $q_{k}$ is the spline order. For computational convenience, we let $N_{k}=N$ and $q_{k}=q$ so that $L_{k}=L$. In practice, cubic splines with order $q=4$ are often used. By the result in MR1900298, the nonparametric function can be approximated well by a spline function such that $g_{k}(X^{\top}\beta_{k})\approx B_{k}(X^{\top}\beta_{k})^{\top}\delta_{k}$ for some $\delta_{k}\in R^{L}$, $k=1,2.$ For notational simplicity, we write $B_{k}(X^{\top}\beta_{k})$ as $B(X^{\top}\beta_{k})$. Therefore, the estimates of the unknown index parameters $\beta_{k}$ and the spline coefficients $\delta_{k}$ are the minimizers of

[TABLE]

Denote $U_{k}=X^{\top}\beta_{k}$ for $k=1,2$. For notational simplicity, we use $U$ to denote $U_{k}$ by suppressing $k$. To allow for possibly unbounded support $\mathcal{U}$ of $U$, according to the method given in CC15, we can weight the B-spline basis functions by a sequence of nonnegative weighting functions $w_{n}:\mathcal{U}\rightarrow\{0,1\}$ given by $w_{n}(u)=1$ if $u\in\mathcal{D}_{n}$ and $w_{n}(u)=0$ otherwise, where $\mathcal{D}_{n}\subseteq\mathcal{U}$ is compact, convex and has nonempty interior and $\mathcal{D}_{n}\subseteq\mathcal{D}_{n+1}$ for all $n$. For the choices of $\mathcal{D}_{n}$, we refer to CC15 for the discussions. When $\mathcal{U}$ is compact, we simply set $w_{n}(u)=1$ for all $u\in\mathcal{U}$. Then we obtain a set of new B-spline basis functions: $B_{k}^{w}(u_{k})=\{B_{\ell,k}(u_{k})w_{n}(u_{k}):1\leq\ell\leq L_{n,k}\}^{\top}$, and replace $B_{k}(u_{k})$ by $B_{k}^{w}(u_{k})$ in the above minimization problem for obtaining the estimators $\widehat{\beta}_{k}$ and $\widehat{\delta}_{k}$. As given in CC15, the asymptotic properties including consistency and asymptotic distribution for the estimator of the unknown function $g_{1}(u)$ still hold for $u\in\mathcal{D}_{n}$.

In principle, we would like to obtain the estimates of $\beta_{k}$ and $\delta_{k}$ by minimizing the penalized negative log-likelihood given in (4). In practice, this minimization is different to achieve, and thus an iterative algorithm (WE83) is often applied. To this end, we obtain the estimates of $\widehat{\beta}_{k}$ and $\widehat{\delta}_{k}$ through an iterative algorithm described as follows, although our theoretical properties given in Section 3.2 are established for the estimators which minimize (4).

Step 1: Given $\beta_{k}$, the solution of $\delta_{k}$ is easily obtained. Reexpressing the model gives

[TABLE]

This can be viewed as a logistic regression model using $\left(ZB(X^{\top}\beta_{k})^{\top},B(X^{\top}\beta_{k})^{\top}\right)^{\top}$ as the regressors without intercept term.

Step 2: Given $\delta_{k}$, it remains to find the solution that minimizes (3) with respect to $\beta_{k}$. Let $\beta_{k}^{old}$ and $\beta_{k,-1}^{old}$ be the current estimates for $\beta_{k}$ and $\beta_{k,-1}$, respectively. Let $\widetilde{g}_{k}(X^{\top}\beta_{k})=B_{k}(X^{\top}\beta_{k})^{\top}\delta_{k}$. We approximate

[TABLE]

where $J(\beta_{k})=\partial\beta_{k}/\partial\beta_{k,-1}=(-\beta_{k,-1}/\sqrt{1-\rVert\beta_{k,-1}\lVert_{2}^{2}},I_{p-1})^{\top}$ is the Jacobian matrix of size $p$ by $p-1$. To obtain the sparse estimates of $\beta_{k}$, we carry out a regularized logistic regression with $\left[\{Z\widetilde{g}_{1}^{\prime}(X^{\top}\beta_{1}^{old})J(\beta_{1}^{old})^{\top}X\}^{\top},\{\widetilde{g}_{2}^{\prime}(X^{\top}\beta_{2}^{old})J(\beta_{2}^{old})^{\top}X\}^{\top}\right]^{\top}$ as the regressors with a known intercept term given as

[TABLE]

This produces an updated vector $\beta_{k,-1}^{new}$. Then we set $\beta_{k}^{new}=(\sqrt{1-||\beta_{k,-1}^{new}||^{2}},(\beta_{k,-1}^{new})^{\top})^{\top}$, for $k=1,2$. Steps 1 and 2 are repeated until convergence.

We obtain the initial value of $\beta_{k}$ through fitting a regularized logistic regression by assuming that $g_{k}(X^{\top}\beta_{k})=X^{\top}\beta_{k}$. In step 2, we use the coordinate descent algorithm (MR2810396) to fit the regularized regression. Moreover, we choose to use the non-convex penalties such as MCP and SCAD which induce nearly unbiased estimators. The MCP (MR2604701) has the form $p_{\gamma}(t,\lambda)=\lambda\int_{0}^{t}(1-x/(\gamma\lambda))_{+}dx,\gamma>1$ and the SCAD (MR1946581) penalty is $p_{\gamma}(t,\lambda)=\lambda\int_{0}^{t}\min\{1,(\gamma-x/\lambda)_{+}/(\gamma-1)\}dx,\gamma>2,$ where $\gamma$ is a parameter that controls the concavity of the penalty functions. In particular, both penalties converge to the $L_{1}$ penalty as $\gamma\rightarrow\infty$. We put $\gamma$ in the subscript to indicate the dependence of these penalty functions on it. In practice, we treat $\gamma$ as a fixed constant. The B-spline basis functions and their derivatives are calculated using the bsplineS function in R package fda.

From the above algorithm, we obtain the spline estimators of the functions $g_{1}(\cdot)$ and $g_{2}(\cdot)$. However, the spline estimator only has convergence rates but its asymptotic distribution is not available in the additive model settings with multiple unknown functions (S85), so no measures of confidence can be assigned to the estimators for conducting statistical inference(MR2382655). The spline-backfitted kernel (SBK) estimator is designed to overcome this issue for generalized additive models (doi:10.1080/01621459.2013.763726), which combines the strengths of kernel and spline smoothing, is easy to implement and has asymptotic distributions. Denote the SBK estimator of $g_{1}(\cdot)$ as $\widehat{g}_{1,sbk}(\cdot)$. We obtain $\widehat{g}_{1,sbk}(X_{0}^{\top}\widehat{\beta}_{1})$ for a new input vector $X_{0}$ as follow:

Step 3: Given the spline estimate $\widehat{g}_{2}(X_{i}^{\top}\widehat{\beta}_{2})$, the loss criterion for a local linear logistic regression can be expressed as the following negative quasi-likelihood function:

[TABLE]

where $K_{h}(\cdot)$ is a kernel function with bandwidth $h$. We obtain the estimate $\widehat{a}(X_{0})$ by minimizing the above loss function. Then the predicted CSTE value at $X_{0}$ is

[TABLE]

Based on the above SBK estimator of CSTE, we can construct a SCB which is used for optimal treatment selection. The details will be discussed in Section 4.

3.2 Asymptotic Analysis

For any positive sequences $\{a_{n}\}$ and $\{b_{n}\}$, let $a_{n}\asymp b_{n}$ denote $\lim_{n\rightarrow\infty}a_{n}b_{n}^{-1}=C$ for a constant $0<C<\infty$ and $a_{n}\ll b_{n}$ denote $\lim_{n\rightarrow\infty}a_{n}b_{n}^{-1}=0$. For a vector $a=(a_{1},\dots,a_{p})^{\top}\in R^{p},$ denote $\|a\|=(\sum_{l=1}^{p}a_{l}^{2})^{1/2}$ and $\|a\|_{\infty}=\max_{l}|a_{l}|$. For a matrix $A=(A_{ij})$, denote $\|A\|=\max_{\|\zeta\|=1}\|A\zeta\|$, $\|A\|_{\infty}=\max_{i}\sum_{j}|A_{ij}|$, and $\|A\|_{2,\infty}=\max_{i}||A_{i}||$, where $A_{i}$ is the $i$th row. For a symmetric matrix $A$, let $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ be the smallest and largest eigenvalues of $A$, respectively. We assume that the nonsparsity size $s=\max(s_{1},s_{2})\ll n$ and the dimensionality satisfies $\log p=O(n^{\alpha})$ for some $\alpha\in(0,1)$. Denote $\beta_{k1}=(\beta_{k1},\dots,\beta_{ks_{k}})^{\top},\beta_{k1,-1}=(\beta_{k2},\dots,\beta_{ks_{k}})^{\top},\beta_{k2}=(\beta_{k(s_{k}+1)},\dots,\beta_{kp})^{\top}.$ Then we write $\beta_{(1)}=(\beta_{11}^{\top},\beta_{21}^{\top})^{\top},\beta_{(2)}=(\beta_{12}^{\top},\beta_{22}^{\top})^{\top},\beta_{(1),-1}=(\beta_{11,-1}^{\top},\beta_{21,-1}^{\top})^{\top}.$ Denote the Jacobian matrix as $J(\beta_{k1})=\partial\beta_{k1}/\beta_{k1,-1},k=1,2$, and $J(\beta_{(1)})=\partial\beta_{(1)}/\partial\beta_{(1),-1}=diag(J(\beta_{11}),J(\beta_{21})),$ which is a block diagonal matrix. We use the superscript ‘0’ to represent the true values.

Denote the first $s_{k}(1\leq k\leq 2)$ components of $X_{i}$ as $X_{i,k1}=(x_{ij},1\leq j\leq s_{k})^{\top}$ and the last $p-s_{k}$ components as $X_{i,k2}=(x_{ij},s_{k}<j\leq p)^{\top}.$ Denote $S(x)=(1+e^{-x})^{-1}$ as the sigmoid function. Then the true expected value of the response given $(X,Z)$ is

[TABLE]

Define the space $\mathcal{M}$ as a collection of functions with finite L2 norm on $\mathscr{C}\times\{0,1\}$ by

[TABLE]

For a given random variable $U$, define its projection onto the space $\mathcal{M}$ as

[TABLE]

where $w^{0}=\pi^{0}(1-\pi^{0})$ and $\pi^{0}=\mu(X,Z)$. For a vector $U=\{U_{1},...,U_{d}\}$, let

[TABLE]

Moreover, we define

[TABLE]

where $\widetilde{X}_{i,kv}=X_{i,kv}-\mathbb{P}_{\mathcal{M}}(X_{i,kv})$ for $k,v=1,2$, and

[TABLE]

where $\sigma^{2}(X,Z)=E[\{Y-S(g(X,Z))\}^{2}|X,Z]$ .

To establish asymptotic properties, we need the following regularity conditions.

Assumption 1**.**

The penalty function $p_{\gamma}(t,\lambda)$ is a non-decreasing symmetric function and concave on $[0,\infty)$. For some constant $a>0$, $\rho(t)=\lambda^{-1}p_{\gamma}(t,\lambda)$ is a constant for all $t\geq a\lambda$ and $\rho(0)=0$. $\rho^{\prime}(t)$ exists and is continuous except for a finite number of $t$ and $\rho^{\prime}(0+)=1$.

Assumption 2**.**

For $k=1,2$, for any given $\beta_{k}\in$ $\Theta$, $g_{k}$ is a nonconstant function on $\{x\in{\mathcal{X}\}}$, where $\mathcal{X}$ is compact, convex and has nonempty interior, and $g_{k}\in\mathcal{H}_{r}$ for some $r>1$, where $\mathcal{H}_{r}$ is the collection of all functions such that the $q^{\text{th}}$ order derivative satisfies the Hölder condition of order $\gamma$ with $r\equiv q+\gamma$ and $q\geq 1$, i.e. for any $\phi\in\mathcal{H}_{r}$, there is a $C_{0}\in(0,\infty)$ such that $|\phi^{(q)}(u_{1})-\phi^{(q)}(u_{2})|\leq C_{0}|u_{1}-u_{2}|^{\gamma}$ for all $u_{1}$ and $u_{2}$ in the support of $g_{k}$.

Assumption 3**.**

For $1\leq j\leq s_{k}$, $E(X_{j}^{2+2(\varkappa+1)})\leq C_{\varkappa,k}$ for some constant $\varkappa\in(0,1)$. For $s_{k}\leq j\leq p-s_{k}$, $E(|X_{j}|^{2+\varrho})\leq C_{\varrho,k}$, for some $\varrho>(8/3)(1-\alpha)^{-1}-2$ and $\varrho\geq 2$.

Assumption 4**.**

There exist $c,c^{\prime},c_{k}\in(0,\infty)$ such that $\lambda_{\min}(\Omega_{n})\geq c$, $\lambda_{\min}(\Phi_{n})\geq c^{\prime}$ almost surely, and $\|E(\widetilde{X}_{k2}\widetilde{X}_{k1}^{\top})\|_{2,\infty}\leq c_{k}$, where $\widetilde{X}_{kv}=X_{kv}-\mathbb{P}_{\mathcal{M}}(X_{kv})$, for $k,v=1,2$.

Assumption 5**.**

Let $w_{n}=2^{-1}\min\{|\beta_{kj}^{0}|:2\leq j\leq s,k=1,2\}.$ Assume that $\max((s/n)^{1/2},\newline L^{1-r}+L^{3/2}n^{\alpha/2-1/2})\ll\lambda\ll w_{n}$.

Assumption 1 is a typical condition on the penalty function, see FLv11. The concave penalties such as SCAD and MCP satisfy Assumption 1. Assumption 2 is a typical smoothness condition on the unknown nonparametric function, see for instance Condition (C3) in 6cfe1e56953040a0bb6a2f0f4c1310c4. Moreover, we require that the unknown functions be nonconstant functions in order to identify the index parameters $\beta_{k}$. Assumption 3 is required for the covariates, see Condition (A5) in MR2780816. Moreover, the design matrix needs to satisfy Assumption 4. A similar condition can be found in FLv11. Assumption 5 assumes that half of the minimum nonzero signal in $\beta_{k}^{0}$ is bounded by some thresholding value, which is allowed to go to zero as $n\rightarrow\infty$. This assumption is needed for variable selection consistency established in Theorem 1.

Denote $\beta_{-1}=(\beta_{11,-1}^{\top},\beta_{12}^{\top},\beta_{21,-1}^{\top},\beta_{22}^{\top})^{\top}$. Let $s=\max(s_{1},s_{2})$. Theorem 1 establishes the consistency for the parameters in model (1).

Theorem 1**.**

Under Assumptions (A1)-(A5), and $\alpha\in(0,(2r-5/2)/(2r-1))$, $L^{1-r}s^{1/2}=o(1)$, $n^{\alpha-1}(L^{3}+s)=o(1)$, $n^{1/2}L^{1-2r}=o(1)$, $\log(n)(s^{1/2}+L^{1/2})L^{3/2}n^{-1/2}=o(1)$, there exists a strict local minimizer $\widehat{\beta}_{-1}=(\widehat{\beta}_{11,-1}^{\top},\widehat{\beta}_{12}^{\top},\widehat{\beta}_{21,-1}^{\top},\widehat{\beta}_{22}^{\top})^{\top}$ of the loss function given in (3) such that $\widehat{\beta}_{k2}=0$ for $1\leq k\leq 2$ with probability approaching $1$ as $n\rightarrow\infty$, and $\|\widehat{\beta}_{-1}-\beta_{-1}^{0}\|=O_{p}(\sqrt{s/n})$.

**Remark. **Based on the assumption given in Theorem 1, the number of spline basis functions needs to satisfy $n^{1/\{2(2r-1)\}}\ll L\ll\min\{n^{1/4}\{\log(n)\}^{-1},n^{(1-\alpha)/3}\}$.

The following Theorem presents the convergence rate for the spline estimator of the unknown functions.

Theorem 2**.**

Under conditions given in Theorem 1, we have $n^{-1}\sum_{i=1}^{n}\{\widehat{g}_{k}(X_{i}^{\top}\widehat{\beta}_{k})-{g}_{k}(X_{i}^{\top}\beta_{k}^{0})\}^{2}=O_{p}(L^{-2r}+L/n+s/n)$, for $k=1,2$.

To estimate the SCB, we need the following assumptions, see Assumptions (A4) and (A6) in Zheng2016.

Assumption 6**.**

Let $r=2$. The kernel function $K$ is symmetric probability density function supported on $[-1,1]$ and has bounded derivative. The bandwidth $h$ satisfies $h=h_{n}=o(n^{-1/5}(\log n)^{-1/5})$ and $h^{-1}=O(n^{1/5}(\log n)^{\delta})$ for some constant $\delta>1/5$.

Assumption 7**.**

The joint density of $X^{\top}\beta_{1}^{0}$ and $X^{\top}\beta_{2}^{0}$ is a bounded and continuous function. The marginal probability density functions have continuous derivatives and the same bound as the joint density. When $\beta_{1}^{0}=\beta_{2}^{0}=\beta^{0}$, the probability density function of $X^{\top}\beta^{0}$ is bounded, continuous, and has continuous derivative.

We borrow some notations in Zheng2016:

[TABLE]

where $\mu(u,Z)=S(g_{1}(u)Z+g_{2}(X^{\top}\beta_{2}^{0}))$ and $f(u)$ is the density function of $X^{\top}\beta_{1}^{0}$. Define the quantile function

[TABLE]

for any $\alpha\in(0,1)$, where $a_{h}=\sqrt{-2\log h}$ and $C_{K}=||K^{\prime}||_{2}^{2}/||K||_{2}^{2}$. Without loss of generality, assume $x^{\top}\beta_{1}^{0}$ is in the range $[0,1]$ and let $\mathscr{C}$ be a set of $x$ such that $\mathscr{C}=\{x:~{}x^{\top}\beta_{1}^{0}\in[h,1-h],x\in R^{p}\}$. Theorem 3 is an adaptation from Theorem 1 in Zheng2016. It provides a method to construct the SCB for $g_{1}$.

Theorem 3**.**

Under Assumptions 1-7, and $\alpha\in(0,2/5)$, $sn^{-1/10}(\log n)^{-3/5}=O(1)$ and $n^{1/5}\ll L\ll\min\{n^{1/4}\{\log(n)\}^{-1},n^{(1-\alpha)/3}\}$, we have

[TABLE]

where $\sigma_{n}(x^{\top}\beta_{1})=n^{-0.5}h^{-0.5}v(x^{\top}\beta_{1})/D(x^{\top}\beta_{1})$.

**Remark 1. **Based on the proof for Theorem 3, we can readily obtain that

$\sigma_{n}(x^{\top}\widehat{\beta}_{1})^{-1}\{\widehat{g}_{1,sbk}(x^{\top}\widehat{\beta}_{1})-g_{1}(x^{\top}\beta_{1}^{0})\}\rightarrow\mathcal{N}(0,1)$ in distribution, for given $x\in\mathscr{C}$. Thus, the $100(1-\alpha)\%$ pointwise confidence interval for $g_{1}(x^{\top}\beta_{1}^{0})$ is $\widehat{g}_{1,sbk}(x^{\top}\widehat{\beta}_{1})\pm Z_{1-\alpha/2}\sigma_{n}(x^{\top}\widehat{\beta}_{1})$, where $Z_{1-\alpha/2}$ is the $(1-\alpha/2)100\%$ quantile of the standard normal distribution.

**Remark 2. **For details of implementations of kernel and spline estimation for functions in (7), we refer to Zheng2016. As suggested in Zheng2016, we use a data-driven under smoothing bandwidth $h=h_{opt}(\log n)^{-1/4}$ and $h_{opt}$ is given in Zheng2016. We let the number of spline interior knots be $\left\lfloor n^{1/5}(\log n)\right\rfloor+1$. The $100(1-\alpha)\%$ SCB for $g_{1}(x^{\top}\beta_{1}^{0})$ is $\widehat{g}_{1,sbk}(x^{\top}\widehat{\beta}_{1})\pm\sigma_{n}(x^{\top}\widehat{\beta}_{1})Q_{h}(\alpha).$

**Remark 3. ** From Theorem 3, we see that the width of the SCB has the order of $\sqrt{\log(h)n/h}$. The width of the pointwise confidence interval given in Remark 1 has the order of $\sqrt{n/h}$. Based on the assumption for the bandwidth $h$, the SCB is wider than the pointwsie confidence interval asymptotically. In deriving the asymptotic distribution of the SBK estimator for the unknown function, it also involves the estimation error from the parameter estimation, which is given in Theorem 1. We assume that $s$ increases with $n$ in a slow rate, so that this estimation error is negligible in the construction of the asymptotic SCB. Our interest focuses on constructing the SCB for the CSTE curve, based on which we make treatment recommendations to a group of new patients. For conducting inference for the parameters, it can be a future interesting research topic to explore.

The theoretical result in Theorem 3 is obtained from the asymptotic extreme value distribution of the SBK estimator for the CSTE curve. This is the key result for obtaining our asymptotic SCB, which achieves the nominal confidence level asymptotically. It is worth noting that we provide different convergence rates for the estimators of the high-dimensional coefficients and the estimators for the unknown functions in Theorems 1 and 2, respectively. This theoretical result enables us to further derive the asymptotic distribution of $\widehat{g}_{1,sbk}$, when the number of important covariates satisfies certain order given in Theorem 3. The SCB provides an important uniform inferential tool for identifying subgroups of patients that benefit from each treatment and make treatment recommendations to the subgroups, while the pointwise confidence interval can only provide treatment selection for a given patient.

4 Treatment Selection via Confidence Bands

In this section, attention is focused on making individualized treatment decision rule for patients. We provide an example to illustrate how to select the optimal treatment based on the CSTE curve and its SCBs. The goal of treatment selection is to find which group of patients will benefit from new treatment based on their covariates. By the definition of CSTE curve, if we assume that the outcome of interest is death, the CSTE curve is the odds ratio of the treatment effect in reducing the probability of death. That is, a positive ${\mathrm{CSTE}}(X)$ value means that the patients will not benefit from new treatment since they may have a higher death rate than patients who receive old treatment. We define the cutoff points as the places where the upper and lower confidence bands equal to 0. Based on these cutoff points, we are able to identify the regions with the positive and negative values of ${\mathrm{CSTE}}(X)$, respectively, so that it will guide us to select the best treatment for a future patient. To summarize, this treatment selection method consists of the following steps:

Step 1**.**

Use the existing dataset to obtain the SBK estimator $\widehat{g}_{1,sbk}(x^{\top}\widehat{\beta}_{1})$ of the CSTE curve and the corresponding SCBs $\widehat{g}_{1,sbk}(x^{\top}\widehat{\beta}_{1})\pm\sigma_{n}(x^{\top}\widehat{\beta}_{1})Q_{h}(\alpha)$ at a given confidence level $100(1-\alpha)\%$.

Step 2**.**

Identify the cutoff points for $x^{\top}\widehat{\beta}_{1}$ and the regions of positive and negative values of the constructed SCBs, as illustrated in Figure 1.

Step 3**.**

For a new patient, we first calculate $\tilde{x}^{\top}\widehat{\beta}_{1}$, where $\tilde{x}$ is the observed value of the covariates for this new patient. Then we make treatment recommendation for this patient based on what region the value of $x^{\top}\widehat{\beta}$ falls into.

We use the following example to illustrate the method of using the SCBs to select optimal treatment for patients. We assume that $X^{\top}\hat{\beta}_{1}$ has a range of $(-4,2)$. In Figure 1, the solid line represents a CSTE curve and dashed lines above and below the curve are the corresponding 95% SCBs. We assume that the outcome variable $Y$ is the indicator of death. As shown in Figure 1, the CSTE is decreasing when the value of $X^{\top}\hat{\beta}_{1}$ is from -4 to -1.6 and it is increasing when the value of $X^{\top}\hat{\beta}_{1}$ is from -1.6 to 2. In general, when the $X^{\top}\hat{\beta}_{1}$ value of a patient is within the range of (-4,-1.6), a larger value implies that the patient more likely benefits from the new treatment than from the old treatment. On the other hand, if the patient’s $X^{\top}\hat{\beta}_{1}$ value falls into (-1.6,2), the new treatment is more beneficial when a smaller value of $X^{\top}\hat{\beta}_{1}$ is observed. Moreover, the results in Figure 1 also imply that we are 95% confident that a new patient with the $X^{\top}\hat{\beta}_{1}$ value falling into the region $[-4,a]$ or $[d,2]$ will benefit from the old treatment, since the lower bands are all above zero in this region. If $X_{0}^{\top}\hat{\beta}_{1}$ falls into region $[b,c]$, then we have 95% confidence that the patient should receive the new treatment, as the upper bands are all below zero in this region. For the intervals $[a,b]$ and $[c,d]$, zero is covered in the confidence band indicating that there is no significant difference between the effects of the new and old treatments.

It is worth nothing that the regions of indifference can also be constructed based on a parametric model using the limiting distribution of the parametric estimators; see more discussions in wu2016set. However, this approach is only directly applicable to the settings with low-dimensional covariates and to each given patient instead of a group of patients. CSTE define a modified version of the CSTE curve using the quantile of the covariate to compare covariate’s capacities for predicting responses to a treatment. Then they use the “best” covariate as a guidance to select treatment. This may be time-consuming when there is a large number of covariates. Our method selects relevant covariates automatically in the estimation procedure and combines information from all covariates through a weighted combination of those covariates with the weights estimated from the data. The weight reflects how important the corresponding covariate is for predicting the response. As a result, our method is more convenient and flexible in selecting optimal treatment.

5 Simulation Study

In this section, we investigate the finite-sample performance of our proposed method via simulated datasets. We run all simulations in R in a linux cluster. We consider the following examples:

Example 1**.**

${\rm logit}(\mu(X,Z))=X^{\top}\beta_{1}(1-X^{\top}\beta_{1})Z+\exp(X^{\top}\beta_{2}).$**

Example 2**.**

${\rm logit}(\mu(X,Z))=(X^{\top}\beta_{1})^{2}Z+\sin(\pi X^{\top}\beta_{2}/2).$**

Example 3**.**

${\rm logit}(\mu(X,Z))=-\exp(X^{\top}\beta_{1})Z/1.5+(X^{\top}\beta_{2})^{2}.$**

Example 4**.**

${\rm logit}(\mu(X,Z))=X^{\top}\beta_{1}(1-X^{\top}\beta_{1})Z+\exp(X^{\top}\beta_{2})+c\cdot\sin(\pi X^{\top}\beta_{3}/2).$**

Example 5**.**

${\rm logit}(\mu(X,Z))=[X^{\top}\beta_{1}(1-X^{\top}\beta_{1})+c\cdot\sin(\pi X^{\top}\beta_{3}/2)]Z+\exp(X^{\top}\beta_{2})$**

The simulated data are generated as follows: the outcome $Y$ is sampled from a binomial distribution with probability of success equals to $\mu(x,z)$; the covariates $X$ are generated from a truncated multivariate normal distribution with mean vector 0, covariance matrix with $\Sigma_{ij}=0.5^{|i-j|}$, and each covariate is truncated by $(-2,2)$; the binary scale covariate $Z$ is sampled from $Binomial(1,0.5)$, which means that each subject is randomly assigned to either control or treatment group. We set $\beta_{1}=(1,1,1,0,\dots,0)^{\prime}/\sqrt{3}$, $\beta_{2}=(1,-2,0,\dots,0)^{\prime}/\sqrt{5}$, and $\beta_{3}=(1,1,1,0,\dots,0)^{\prime}/\sqrt{3}$. For examples 1-3, we set $p=10,50,100,500$, and sample size $n=500,750,1000$. For examples 4-5, we consider consider that $g_{2}$ and $g_{1}$ are mis-specified respectively. We choose $p=50$, $n=1000$, and $c=0.1,0.2,0.5$. For each pair of $n$ and $p$, we repeat the simulations $J=300$ times.

To obtain sparse estimates of $\beta_{k}$ for high-dimensional cases, we choose SCAD as the penalty function and let $\gamma=3.7$ (MR1946581). The optimal tuning parameter $\lambda$ is chosen from a geometrically increasing sequence of 30 parameters by minimizing the modified Bayesian information criterion (wang2009shrinkage): $BIC_{\lambda}=-2\cdot\log(Loss)+df_{\lambda}\cdot\log(n)\cdot C(p),$ where $Loss$ is the loss function in (3), $df_{\lambda}$ is the number of non-zero elements in $\beta_{k}$, $C(p)=C(p)=\log\log(p)$.

We evaluate the following metrics for the non-parametric function $g_{1}(\cdot)$ based on 200 equally spaced grid points $t_{1},\dots,t_{200}$ on the range of $\hat{\eta}_{1}=X^{\top}\hat{\beta}_{1}$ : the mean squared error (MSE) of its estimator; the mean absolute error (MAE) of its estimator; the average coverage probability (CP) of its SCBs. Since the models in Example 1-3 are correcly specified, we compute the average number of parameters that are incorrectly estimated to be non-zero, the average number of parameters that are incorrectly estimated as zero; the proportions that all relevant covariates are correctly selected and the proportions that some relevant covariates are not selected.

We define that the oracle estimator of $g_{1}(\cdot)$ is obtained when the true indexes of non-zero components in $\beta_{1}$ and $\beta_{2}$ are given in example 1-3. The results are summarized in Table 1 and 2. In Table 1, we see that the coverage probabilities are slightly less than 95% but close to 90% when the sample size is 750. When the sample is 1000, the empirical coverage is close to the nominal 95% confidence level. The mean square error and mean absolute error also decrease as sample size $n$ increases. This indicates that the estimates of confidence bands become more accurate as the sample size increases. Table 1 shows that the proposed asymptotic SCB performs well when the sample size is relatively large. However, with small sample size, a bootstrap procedure is recommended to construct the SCB. This can be a future research topic to explore. In Table 2, when $c$ is small (i.e., the degree of misspecification is small), the numerical results are comparable to those for $c=0$ with no misspecification. However, the performance deteriorates when $c$ becomes large. Moreover, the method is more robust to the misspecification of $g_{2}$ than that of $g_{1}$. The model selection results are summarized in Table 4. The false positive and false negative are decreasing when $n$ increases. It is worth noting that in practice we often need to make a reliable treatment recommendation strategy to a group of patients instead of a specific given patient, and the SCBs can serve the former purpose well. In Table 3, we also compare the performance of the pointwise confidence interval (CI) with that of the SCB for a group of randomly generated new patients with group size 1, 10, 100, 200, 500, respectively. The CI and SCB are constructed based on the data with sample $n=1000,p=10$. We see that the pointwise CI has low empirical coverage rates when the group size is greater than 1. The pointwise CI can work well when we only have one patient. However, when we have a group of patients with different covariate values, we can make an incorrect treatment decision based on the pointwise CI for these patients.

6 A Real Data Example

In this section, we present and discuss the results of applying the procedure described in previous sections to a real data set. We illustrate the applications of the CSTE curve in a real-world example and provide some insight into the interpretation of the results. The goal is to analyze the effect of Zhengtianwan in the treatment of migraines. This is a multicenter, randomized, double-blind, placebo-controlled trial on the effectiveness of Zhengtian pill on treating patients with migraines. A migraine is a common neurological disorder characterized by recurrent headache attacks. Migraine treatment involves acute and prophylactic therapy. The objective of this study is to evaluate the efficacy and safety of Zhengtian Pill for migraine prophylaxis. Zhengtian Pill, a Chinese Patented Medicine approved by the State Food and Drug Administration of China in 1987, has been used in clinical practice for more than 20 years in China to stimulate blood circulation, dredge collaterals, alleviate pains, and nourish the liver. To evaluate the effectiveness of Zhengtian Pill in preventing the onset of migraine attacks, a large-scale, randomized and prospective clinical study was conducted. Eligible patients were monitored during a baseline period of four weeks, during which the headache characteristics were recorded as baseline data. In this period, any use of migraine preventive medications was prohibited. After the baseline period, a 12-week treatment period and four-week follow-up period were carried out. Patients were requested to keep a headache diary throughout the whole study period, from which investigators were able to extract detailed information of migraine attacks including migraine days, frequency, duration, and intensity as well as the use of acute medication during the study period. The outcome measures were evaluated at 4, 8, and 12 weeks, and during the follow-up period. The patients who met the inclusion criteria were randomly assigned into the experimental group and control group in a 1:1 ratio using a computer-generated stochastic system.

In our analysis, the response variable $Y$ is a binary outcome indicating if the number of days that headaches occur has decreased 8 weeks after patients were treated. $Z$ is another indicator: $Z=1$ means the subject is assigned in the experimental group; $Z=0$ means in the control group. The covariates $x_{1}$ to $x_{3}$ are gender (0 for male, 1 for female), height and body weight, $x_{4}$ to $x_{12}$ are overall scores of Traditional Chinese Medicine Syndromes (Huozheng, Fengzheng, Xueyu, Tanshi, Qixu, Yuzheng, Xuexu, Yinxu, and Yangxu) for migraine. Those scores, with a range of 0 to 20, are commonly used to describe the severity of migraine symptoms in Traditional Chinese Medicine. The higher of the score means that the syndrome is more severe. All covariates are centered and standardized as input. We have 204 observations where 99 are in the experimental group and 105 are in the control group. The purpose of this exercise is to model the odds ratio as a function of those covariates. The number of covariates $p=12$ might be regarded as small, so we estimate the CSTE curve using the algorithm in Section 3.1 with and without model selection. We use SCAD as our penalty function and the optimal tuning parameter is selected via the modified BIC criterion. The corresponding SCBs for the CSTE curve are calculated. The results are not intended as definitive analyses of these data.

Table 5 summarizes the point estimates of $\beta_{1}$ and the corresponding standard errors. Figure 2 shows two estimated CSTE curve and their SCBs: (a) using all 12 variables; (b) using all variables except $x_{7}$ and $x_{8}$ (not selected). To aid in interpretation for each covariate, we depict each covariate versus $g_{1}$ where other covariates are projected onto their mean values in Figure 3. As we can see, Huozheng has a monotonic dependence on the CSTE, but the corresponding relationships for other overall scores are highly nonlinear; most of them have a quadratic appearance. On the basis of these Figures, we can conclude that the odds ratio does depend on the linear combination of the covariates in a nonlinear manner.

When all variables are used to estimate the curve, the two cutoff points are $c_{1}=-0.502,c_{2}=2.182$. The estimates of the overall rating of biomarker values are $x^{\top}\hat{\beta}_{1}$ where the majority of the points ($>$ 95%) fall into the interval $(-1.2,5)$. Two cutoff points divide this interval into three parts. Since the response variable $y=1$ represents headache improvements after 8 weeks treatment, higher CSTE value means that the patient more likely benefits from the treatment. Suppose a new patient with biomarker value $x_{0}$. When $x^{\top}\hat{\beta}_{1}$ falls into $[-1.2,c_{1}]$ and $[c_{2},5]$, the treatment does not have a significant effect. When it falls into $(c_{1},c_{2})$, then this treatment can improve the migraine of this patient. In the Zhengtianwan data, there are 37.7% of patients fall into the interval. On the other hand, when $x_{7},x_{8}$ are removed from the model, we obtain a new estimate $\hat{\beta}_{1,ms}$ of $\beta_{1}$. As we can see from Figure 2 (b), the positive region of the CSTE curve becomes wider than that in (a), and the two cutoff points become -0.811 and 2.366. The majority of $x^{\top}\hat{\beta}_{1,ms}$ falls into $-(2.5,5)$. When $x^{\top}\hat{\beta}_{1,ms}$ is in the range of $[-0.811,2.366]$, the patient should receive Zhengtianwan treatment and 50.98% of all subjects in the training dataset fall into this interval.

The number of covariates is not very large in this dataset, so it is possible to fit a sparse linear model with high order terms and make treatment recommendation based on the parametric estimate of the CSTE curve. To this end, we model $g_{k}$ for $k=1,2$ as a multivariate polynomial function of the covariates $X$ with degree of two, i.e. $g_{k}(X)=\alpha_{0,k}+\sum_{j}\alpha_{j,k}X_{j}+\sum_{j\leq j^{\prime}}\alpha_{jj^{\prime},k}X_{j}X_{j^{\prime}}.$ We obtain the estimates ($\hat{\alpha}_{0,k}$,$\hat{\alpha}_{j,k}$,$\hat{\alpha}_{jj^{\prime},k}$) of the unknown coefficients ($\alpha_{0,k}$,$\alpha_{j,k}$,$\alpha_{jj^{\prime},k}$) by fitting regularized $L_{1}$ logistic regression with the tuning parameter selected by 3-fold cross-validation. For a new patient with the observed covariates $\tilde{x}$, we obtain the predicted value for the CSTE curve as $\hat{g_{1}}(\tilde{x})=\hat{\alpha}_{0,1}+\sum_{j}\hat{\alpha}_{j,1}\tilde{x}_{j}+\sum_{j\leq j^{\prime}}\hat{\alpha}_{jj^{\prime},1}\tilde{x}_{j}\tilde{x}_{j^{\prime}}$. For this parametric modeling method, we can only use the predicted value for $g_{1}$ to make treatment recommendation. Thus, based on the rule that a subject would benefit from the treatment if $\hat{g}_{1}$ is positive, $67.64\%$ of the subjects should receive the treatment. We see that this percentage is higher than that obtained by our proposed method, as the parametric method only uses the point estimate for treatment selection. The construction of confidence intervals for $g_{1}$ by the parametric method is very challenging, as it involves finding the joint asymptotic distribution of ($\hat{\alpha}_{0,k}$,$\hat{\alpha}_{j,k}$,$\hat{\alpha}_{jj^{\prime},k}$) after model selection, which is a difficult problem. Moreover, SCBs are not available based on the parametric modeling method. As discussed and shown in the simulation studies given in Section 5, SCBs are more reliable than the pointwise CIs for treatment recommendation for a group of patients.

7 Discussion

Both the simulation and real-world studies in Sections 5 and 6 suggest that the modeling procedure for CSTE curve can successfully detect and model complicated non-linear relationships between binary response and high-dimensional covariates. In practice, the non-linear dependencies we suggest are not characteristic of all situations. We adapt the spline-backfitted kernel smoothing to construct the SCBs for the non-linear functions to choose the optimal treatment. Moreover, the SCBs can be used to verify the presence of non-linear relationships as well.

Our model is motivated by the desire to provide an individualized decision rule for patients along with the ability to deal with high-dimensional covariates when the outcome is binary. The semi-parametric modeling approach can be viewed as a generalization of the CSTE curve with one covariate proposed in CSTE in the sense that the odds ratio depends on a weighted linear combination of all covariates. Although we consider a single decision with two treatment options, our model can be readily generalized to multiple treatment arms.

To identify model (1), we require that the unknown functions $g_{1}\left(\cdot\right)$ and $g_{2}\left(\cdot\right)$ be nonconstant functions on their supports. Otherwise, the index parameters $\beta_{1}$ and $\beta_{2}$ are not identified. That is, if $g_{k}\left(\beta_{k}^{\top}x\right)=a_{k}$ for all $x\in{\cal X}$, then $\beta_{k}$ can take an arbitrary value. Our treatment recommendation method is proposed based on the prior information that at least one baseline covariate is useful, and the treatment effect is not zero for some values of the covariates. Based on this prior information, our goal is to estimate the optimal treatment regimes and provide recommendations for individual patients based on their observed baseline covariates. For example, Figures 1-2 show that the estimated values of $g_{1}\left(\cdot\right)$ are nonzero in some regions, so that one can apply our method to identify those regions and then make treatment recommendations to patients. In practice, if we conjecture that all covariates are irrelevant, then for an arbitrarily given value of $\beta_{k}$, $g_{k}\left(\beta_{k}^{\top}x\right)=a_{k}$ for all $x\in{\cal X}$. We can conduct a test for this conjecture by choosing a value for $\beta_{k}$, estimating the unknown function $g_{k}\left(\beta_{k}^{\top}x\right)$ and then testing whether the function is a constant or not. Likewise, we can use this approach to testing whether $g_{1}\left(\beta_{1}^{\top}x\right)=0$ for all $x\in{\cal X}$, i.e. whether the treatment has effect at all. If they are rejected, next we can fit our model to more precisely find the optimal value of $\beta_{k}$. Since this paper focuses on making individualized treatment recommendations to patients based on their available covariates, we leave this interesting topic as a future work to explore.

SUPPLEMENTARY MATERIAL

Supplementary material includes the technical proofs for all the theoretical results.