Selecting Biomarkers for building optimal treatment selection rules using Kernel Machines

TL;DR

This paper develops methods for selecting optimal biomarker combinations for treatment decisions by incorporating measurement costs into a feature selection framework using kernel machines, improving model efficiency and cost-effectiveness.

Contribution

It introduces a novel feature selection approach with cost considerations for treatment-rule optimization using kernel machine techniques, addressing cost and performance trade-offs.

Findings

Feature selection improves treatment rule accuracy

Incorporating biomarker costs reduces measurement expenses

Methods outperform existing approaches in simulations and real data

Abstract

Optimal biomarker combinations for treatment-selection can be derived by minimizing total burden to the population caused by the targeted disease and its treatment. However, when multiple biomarkers are present, including all in the model can be expensive and hurt model performance. To remedy this, we consider feature selection in optimization by minimizing an extended total burden that additionally incorporates biomarker measurement costs. Formulating it as a 0-norm penalized weighted classification, we develop various procedures for estimating linear and nonlinear combinations. Through simulations and a real data example, we demonstrate the importance of incorporating feature-selection and marker cost when deriving treatment-selection rules.

Figures

Tables

Table 1. Table 1: Setting 1(i) (Linear): Setting with disease prevalence approximately 0.30 and 0.20 among the untreated and treated, respectively (Number of markers 27, number of significant markers 2) - Performance scores: (a) Monte Carlo mean of proportion of correct markers chosen, (b) Monte Carlo mean of proportion of incorrect markers chosen, (c) Monte Carlo mean of θ 𝜃 \theta with (i) δ 1 = 0 subscript 𝛿 1 0 \delta_{1}=0 and δ 2 = 0 subscript 𝛿 2 0 \delta_{2}=0 , (ii) with δ 1 = 0 subscript 𝛿 1 0 \delta_{1}=0 and δ 2 = 0.0001 subscript 𝛿 2 0.0001 \delta_{2}=0.0001 (iii) with δ 1 = 0 subscript 𝛿 1 0 \delta_{1}=0 and δ 2 = 0.001 subscript 𝛿 2 0.001 \delta_{2}=0.001 for different feature selection methods .

Setting 1 (Underlying model: Linear)
	${\delta }_1=0$, ${\delta }_2=0$			${\delta }_1=0$, ${\delta }_2=0.0001$			${\delta }_1=0$, ${\delta }_2=0.001$
Linear Methods	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$
	correct)	incorrect)		correct)	incorrect)		correct)	incorrect)
SCAD	1	0.17	0.1048	1	0.16	0.1069	1	0.01	0.1045
elastic SCAD	1	0.20	0.1076	1	0.18	0.1128	1	0.11	0.1145
$L_1$ WSVM	1	0.31	0.1107	1	0.29	0.1123	1	0.16	0.1109
$L_0$ WSVM	1	0.26	0.1114	1	0.16	0.1105	1	0.08	0.1127
FSV	1	0.32	0.1074	1	0.22	0.1096	1	0.11	0.1075
AROM	0.97	0.10	0.1059	0.97	0.09	0.1079	0.97	0.05	0.1116
riskRFE (L)	1	0.03	0.1016	1	0.03	0.1019	1	0.03	0.1044
LASSO	1	0.29	0.1065	1	0.29	0.1074	1	0.29	0.1158
W-SVM Linear no sel.†	1	1	0.1178	1	1	0.1205	1	1	0.1448
	${\delta }_1=0$, ${\delta }_2=0$			${\delta }_1=0$, ${\delta }_2=0.0001$			${\delta }_1=0$, ${\delta }_2=0.001$
Nonlinear Methods	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$
	correct)	incorrect)		correct)	incorrect)		correct)	incorrect)
KNIFE	1	0.27	0.1075	1	0.26	0.1101	1	0.12	0.1120
$L_1$ KP-WSVM	1	0.24	0.1080	1	0.24	0.1094	1	0.09	0.1082
$L_0$ KP-WSVM	1	0.32	0.1113	1	0.27	0.1102	1	0.17	0.1127
riskRFE (G)	1	0.03	0.1021	1	0.03	0.1024	1	0.03	0.1049
W-SVM Gauss no sel.†	1	1	0.1202	1	1	0.1229	1	1	0.1472

Table 2. Table 2: Setting 1(ii) (Nonlinear): Setting with disease prevalence approximately 0.18 and 0.17 among the untreated and treated, respectively (Number of markers 27, number of significant markers 2) - Performance scores: (a) Monte Carlo mean of proportion of correct markers chosen, (b) Monte Carlo mean of proportion of incorrect markers chosen, (c) Monte Carlo mean of θ 𝜃 \theta with (i) δ 1 = 0 subscript 𝛿 1 0 \delta_{1}=0 and δ 2 = 0 subscript 𝛿 2 0 \delta_{2}=0 , (ii) with δ 1 = 0 subscript 𝛿 1 0 \delta_{1}=0 and δ 2 = 0.0001 subscript 𝛿 2 0.0001 \delta_{2}=0.0001 (iii) with δ 1 = 0 subscript 𝛿 1 0 \delta_{1}=0 and δ 2 = 0.001 subscript 𝛿 2 0.001 \delta_{2}=0.001 for different feature selection methods.

Setting 1 (Underlying model: Nonlinear)
	${\delta }_1=0$, ${\delta }_2=0$			${\delta }_1=0$, ${\delta }_2=0.0001$			${\delta }_1=0$, ${\delta }_2=0.001$
Linear Methods	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$
	correct)	incorrect)		correct)	incorrect)		correct)	incorrect)
SCAD	0.91	0.42	0.1296	0.89	0.32	0.1292	0.85	0.15	0.1322
elastic SCAD	0.96	0.36	0.1289	0.96	0.33	0.1291	0.95	0.23	0.1340
$L_1$ WSVM	0.94	0.43	0.1281	0.94	0.41	0.1296	0.90	0.29	0.1358
$L_0$ WSVM	0.88	0.30	0.1247	0.86	0.29	0.1269	0.82	0.16	0.1314
FSV	0.89	0.32	0.1274	0.90	0.26	0.1266	0.88	0.19	0.1328
AROM	0.91	0.25	0.1269	0.92	0.23	0.1263	0.89	0.16	0.1317
riskRFE (L)	0.85	0.10	0.1223	0.85	0.10	0.1228	0.85	0.10	0.1266
LASSO	0.79	0.19	0.1320	0.79	0.19	0.1326	0.79	0.19	0.1383
W-SVM Linear no sel.†	1	1	0.1349	1	1	0.1376	1	1	0.1619
	${\delta }_1=0$, ${\delta }_2=0$			${\delta }_1=0$, ${\delta }_2=0.0001$			${\delta }_1=0$, ${\delta }_2=0.001$
Nonlinear Methods	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$
	correct)	incorrect)		correct)	incorrect)		correct)	incorrect)
KNIFE	0.99	0.35	0.1156	0.97	0.30	0.1147	0.95	0.21	0.1196
$L_1$ KP-WSVM	0.96	0.33	0.1130	0.96	0.31	0.1126	0.93	0.20	0.1183
$L_0$ KP-WSVM	0.97	0.45	0.1191	0.97	0.42	0.1196	0.95	0.35	0.1245
riskRFE (G)	0.88	0.11	0.1178	0.88	0.11	0.1183	0.88	0.11	0.1222
W-SVM Gauss no sel.†	1	1	0.1356	1	1	0.1383	1	1	0.1626

Table 3. Table 3: Setting 2(i) (Linear): Setting with disease prevalence approximately 0.22 among both the untreated and treated (Number of markers 53, number of significant markers 3) - Performance scores: (a) Monte Carlo mean of proportion of correct markers chosen, (b) Monte Carlo mean of proportion of incorrect markers chosen, (c) Monte Carlo mean of θ 𝜃 \theta with (i) δ 1 = 0 subscript 𝛿 1 0 \delta_{1}=0 and δ 2 = 0 subscript 𝛿 2 0 \delta_{2}=0 for different feature selection methods.

Setting 2 (Underlying model: Linear)
	${\delta }_1=0$, ${\delta }_2=0$			${\delta }_1=0$, ${\delta }_2=0.0001$			${\delta }_1=0$, ${\delta }_2=0.001$
Linear Methods	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$
	correct)	incorrect)		correct)	incorrect)		correct)	incorrect)
SCAD	0.98	0.31	0.0557	0.97	0.30	0.0582	0.95	0.01	0.0545
elastic SCAD	1	0.30	0.0599	1	0.23	0.0588	1	0.13	0.0642
$L_1$ WSVM	1	0.45	0.0605	1	0.33	0.0574	1	0.14	0.0634
$L_0$ WSVM	0.97	0.23	0.0563	0.96	0.17	0.0565	0.95	0.01	0.0542
FSV	0.95	0.36	0.0594	0.94	0.24	0.0627	0.95	0.12	0.0630
AROM	0.94	0.06	0.0562	0.94	0.04	0.0549	0.94	0.02	0.0598
riskRFE (L)	0.97	0.02	0.0519	0.97	0.02	0.0523	0.97	0.02	0.0558
LASSO	1	0.17	0.0555	1	0.17	0.0567	1	0.17	0.0670
W-SVM Linear no sel. †	1	1	0.0686	1	1	0.0739	1	1	0.1216
	${\delta }_1=0$, ${\delta }_2=0$			${\delta }_1=0$, ${\delta }_2=0.0001$			${\delta }_1=0$, ${\delta }_2=0.001$
Nonlinear Methods	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$
	correct)	incorrect)		correct)	incorrect)		correct)	incorrect)
KNIFE	0.98	0.24	0.0557	0.96	0.16	0.0567	0.96	0.03	0.0561
$L_1$ KP-WSVM	0.98	0.19	0.0550	0.98	0.08	0.0539	0.97	0.03	0.0560
$L_0$ KP-WSVM	0.99	0.20	0.0606	0.99	0.17	0.0572	0.99	0.09	0.0630
riskRFE (G)	0.98	0.02	0.0505	0.98	0.02	0.0509	0.98	0.02	0.0545
W-SVM Gauss no sel. †	1	1	0.0681	1	1	0.0734	1	1	0.1211

Table 4. Table 4: Setting 2(ii) (Nonlinear): Setting with disease prevalence approximately 0.30 and 0.18 among the untreated and treated, respectively (Number of markers 53, number of significant markers 3) - Performance scores: (a) Monte Carlo mean of proportion of correct markers chosen, (b) Monte Carlo mean of proportion of incorrect markers chosen, (c) Monte Carlo mean of θ 𝜃 \theta with (i) δ 1 = 0 subscript 𝛿 1 0 \delta_{1}=0 and δ 2 = 0 subscript 𝛿 2 0 \delta_{2}=0 , (ii) with δ 1 = 0 subscript 𝛿 1 0 \delta_{1}=0 and δ 2 = 0.0001 subscript 𝛿 2 0.0001 \delta_{2}=0.0001 (iii) with δ 1 = 0 subscript 𝛿 1 0 \delta_{1}=0 and δ 2 = 0.001 subscript 𝛿 2 0.001 \delta_{2}=0.001 for different feature selection methods.

Setting 2 (Underlying model: Nonlinear)
	${\delta }_1=0$, ${\delta }_2=0$			${\delta }_1=0$, ${\delta }_2=0.0001$			${\delta }_1=0$, ${\delta }_2=0.001$
Linear Methods	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$
	correct)	incorrect)		correct)	incorrect)		correct)	incorrect)
SCAD	0.84	0.65	0.1578	0.74	0.39	0.1588	0.68	0.24	0.1693
elastic SCAD	0.88	0.63	0.1607	0.81	0.47	0.1641	0.77	0.39	0.1821
$L_1$ WSVM	0.76	0.47	0.1597	0.71	0.37	0.1610	0.58	0.18	0.1696
$L_0$ WSVM	0.74	0.33	0.1576	0.69	0.27	0.1621	0.63	0.15	0.1685
FSV	0.83	0.64	0.1597	0.70	0.40	0.1630	0.69	0.33	0.1762
AROM	0.74	0.33	0.1564	0.68	0.24	0.1584	0.65	0.17	0.1653
riskRFE (L)	0.64	0.22	0.1564	0.64	0.22	0.1577	0.64	0.22	0.1693
LASSO	0.59	0.15	0.1636	0.59	0.15	0.1645	0.59	0.15	0.1729
W-SVM Linear no sel.†	1	1	0.1674	1	1	0.1727	1	1	0.2204
	${\delta }_1=0$, ${\delta }_2=0$			${\delta }_1=0$, ${\delta }_2=0.0001$			${\delta }_1=0$, ${\delta }_2=0.001$
Nonlinear Methods	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$
	correct)	incorrect)		correct)	incorrect)		correct)	incorrect)
KNIFE	0.79	0.31	0.1435	0.74	0.21	0.1388	0.71	0.07	0.1304
$L_1$ KP-WSVM	0.81	0.53	0.1528	0.62	0.27	0.1513	0.50	0.10	0.1455
$L_0$ KP-WSVM	0.82	0.53	0.1535	0.81	0.47	0.1567	0.77	0.38	0.1642
riskRFE (G)	0.65	0.09	0.1369	0.65	0.09	0.1375	0.65	0.09	0.1432
W-SVM Gauss no sel.†	1	1	0.1645	1	1	0.1698	1	1	0.2175

Table 5. Table 5: Setting 3 (Nonlinear): Setting with disease prevalence approximately 0.53 and 0.17 among the untreated and treated, respectively (Number of markers 27, number of significant markers 2) - Performance scores: (a) Monte Carlo mean of proportion of correct markers chosen, (b) Monte Carlo mean of proportion of incorrect markers chosen, (c) Monte Carlo mean of θ 𝜃 \theta with (i) δ 1 = 0 subscript 𝛿 1 0 \delta_{1}=0 and δ 2 = 0 subscript 𝛿 2 0 \delta_{2}=0 for different feature selection methods.

Setting 3 (Underlying model: Nonlinear)
	${\delta }_1=0$, ${\delta }_2=0$			${\delta }_1=0$, ${\delta }_2=0.0001$			${\delta }_1=0$, ${\delta }_2=0.001$
Linear Methods	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$
	correct)	incorrect)		correct)	incorrect)		correct)	incorrect)
SCAD	NA	NA	NA	NA	NA	NA	NA	NA	NA
elastic SCAD	NA	NA	NA	NA	NA	NA	NA	NA	NA
$L_1$ WSVM	0.19	0.18	0.1744	0.17	0.20	0.1749	0.16	0.19	0.1795
$L_0$ WSVM	1	0.99	0.1746	1	1	0.1773	1	1	0.2015
FSV	0.91	0.95	0.1745	0.50	0.52	0.1746	0.50	0.52	0.1873
AROM	NA	NA	NA	NA	NA	NA	NA	NA	NA
riskRFE (L)	0.22	0.29	0.1749	0.22	0.29	0.1757	0.22	0.29	0.1826
LASSO	0.05	0.14	0.1776	0.05	0.14	0.1780	0.05	0.14	0.1812
W-SVM Linear no sel.†	1	1	0.1746	1	1	0.1773	1	1	0.2015
	${\delta }_1=0$, ${\delta }_2=0$			${\delta }_1=0$, ${\delta }_2=0.0001$			${\delta }_1=0$, ${\delta }_2=0.001$
Nonlinear Methods	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$	$𝔼$(Prop.	$𝔼$(Prop.	$𝔼\left(\theta \right)$
	correct)	incorrect)		correct)	incorrect)		correct)	incorrect)
KNIFE	0.42	0.69	0.1765	0.39	0.68	0.1799	0.40	0.71	0.1958
$L_1$ KP-WSVM	0.88	0.32	0.1100	0.88	0.28	0.1060	0.88	0.26	0.1127
$L_0$ KP-WSVM	1	0.88	0.1621	1	0.89	0.1667	1	0.88	0.1860
riskRFE (G)	0.37	0.08	0.1349	0.37	0.08	0.1352	0.37	0.08	0.1376
W-SVM Gauss no sel.†	1	1	0.1761	1	1	0.1788	1	1	0.2031

Table 6. Table 6: Cross-validated treatment-selection performance of various approaches for making treatment recommendation in the RV144 trial: Disease cost + + vaccine cost + + SNP cost - expected disease rate per 1000 individuals with added cost of vaccine ( δ 1 = 0.001 subscript 𝛿 1 0.001 \delta_{1}=0.001 ) and 3 different costs for marker evaluation ( δ 2 ∈ { 10 − 6 , 5 × 10 − 6 , 10 − 5 } subscript 𝛿 2 superscript 10 6 5 superscript 10 6 superscript 10 5 \delta_{2}\in\{10^{-6},5\times 10^{-6},10^{-5}\} for each of the 28 SNPs evaluated) for treatment selection with various methods of marker selection.

	Disease cost	Disease cost	Disease cost
Methods	$+$ vaccine cost	$+$ vaccine cost	$+$ vaccine cost
	$+$SNP cost	$+$SNP cost	$+$SNP cost
	${\delta }_2={10}^{-6}$	${\delta }_2=5\times {10}^{-6}$	${\delta }_2={10}^{-5}$
Treat all	7.44	7.44	7.44
Treat none	8.54	8.54	8.54
SCAD	7.21	7.28	7.37
elastic SCAD	7.34	7.44	7.56
$L_1$ WSVM	7.24	7.32	7.42
$L_0$ WSVM	7.31	7.37	7.44
FSV	7.29	7.35	7.43
AROM	6.97	7.03	7.11
riskRFE (L)	7.15	7.17	7.19
LASSO	7.53	7.56	7.60
W-SVM Linear no sel.†	7.42	7.55	7.71
KNIFE	7.33	7.38	7.44
$L_1$ KP-WSVM	7.39	7.46	7.55
$L_0$ KP-WSVM	7.32	7.42	7.55
riskRFE (G)	7.35	7.40	7.46
W-SVM Gauss no sel.†	7.44	7.57	7.73

Equations

$$\theta$$

$$=\sum _{a=0}^{1}E\left[Y(a)\times I\{A(\text{X})=a\}\right]+{\delta }_1\times E\{A(\text{X})\}+{\delta }_2\times \dim (\text{X}).$$

$$f_0=\arg \underset{\text{X}\subseteq {\text{X}}^p}{\min }\underset{f\in {\mathcal{F}}_{\text{X}}}{\min }\theta$$

$$\theta$$

$$=E(Y|T=0)-E\left[A(\text{X})\times \{Ris{k}_0(\text{X})-Ris{k}_1(\text{X})-{\delta }_1\}\right]+{\delta }_2\times \dim (\text{X})$$

$$\hat{f}=\arg \underset{\text{X}\subseteq {\text{X}}^p}{\min }\underset{f\in {\mathcal{F}}_{\text{X}}}{\min }\left\{\sum _{i=1}^n\{Y_i\times \frac{T_i}{n_1}-Y_i\times \frac{(1-T_i)}{n_0}+\frac{{\delta }_1}{n}\}\times I\{f({\text{X}}_i)\le 0\}+{\delta }_2\times \dim (\text{X})\right\}.$$

$$\sum _{i=1}^n{W}_{i}I\left\{f({\text{X}}_i)\le 0\right\}/n+{\delta }_2\times \dim (\text{X}),$$

$$\sum _{i=1}^n|W_i|I\left\{\text{sgn}(f({\text{X}}_i))\ne \text{sgn}(W_i)\right\}/n+{\delta }_2\times \dim (\text{X})$$

$$\underset{\text{X}\subseteq {\text{X}}^p}{\min }\underset{f\in {\mathcal{F}}_{\text{X}}}{\min }\sum _{i=1}^n|W_i|\max (1-\text{sgn}(W_i)\times f({\text{X}}_i),0)/n+{\delta }_2\times \dim (\text{X})$$

$$\underset{f\in {\mathcal{H}}_{\text{X}}}{\min }\sum _{i=1}^n|W_i|\max (1-\text{sgn}(W_i)\times f({\text{X}}_i),0)/n+\lambda \Vert f{\Vert }_{{\mathcal{H}}_{\text{X}}}^2$$

$$\underset{\mathbf{\alpha }}{\max }\sum _{i=1}^n{\alpha }_i-\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{\alpha }_i{\alpha }_j{Y}_i^{\ast }{Y}_j^{\ast }{k}_{\text{X}}({\text{x}}_i,{\text{x}}_j)\text{subject to }\sum _{i=1}^n{Y}_i^{\ast }{\alpha }_i=0,0\le {\alpha }_i\le C,$$

$$\underset{\mathbf{\beta },{\beta }_0}{\min }\sum _{i=1}^n{q}_i\max \{1-Y_i^{\ast }(⟨X_i,\mathbf{\beta }⟩+{\beta }_0),0\}/n+\lambda \sum _{j=1}^p|{\beta }_j|.$$

$$\underset{\mathbf{\beta },{\beta }_0}{\min }\sum _{i=1}^n{q}_i\max \{1-Y_i^{\ast }(⟨X_i,\mathbf{\beta }⟩+{\beta }_0),0\}/n+\sum _{j=1}^p{p}_{\lambda }(|{\beta }_j|),$$

$$\underset{\mathbf{\beta },{\beta }_0}{\min }\sum _{i=1}^n{q}_i\max \{1-Y_i^{\ast }(⟨X_i,\mathbf{\beta }⟩+{\beta }_0),0\}/n+\sum _{j=1}^p{p}_{{\lambda }_1}(|{\beta }_j|)+{\lambda }_2\Vert \mathbf{\beta }{\Vert }^2.$$

$$\underset{\mathbf{\beta },{\beta }_0,\mathbf{v}}{\min }$$

$$+\sum _j\alpha e^{-\alpha v_j^{\ast }(v_j-v_j^{\ast })}\text{subject to }-\mathbf{v}\le \mathbf{\beta }\le \mathbf{v},$$

$$\underset{\mathbf{\beta },{\beta }_0}{\min }\sum _{i=1}^n{q}_i\max \{1-Y_i^{\ast }(⟨z_i\ast X_i,\mathbf{\beta }⟩+{\beta }_0),0\}/n+\lambda \sum _{j=1}^p|{\beta }_j{|}^2$$

$$({\mathbf{\beta }}^{(t)},{\beta }_0^{(t)})=\arg \underset{\mathbf{\beta },{\beta }_0}{\min }\sum _{i=1}^n{q}_i{\xi }_i/n+\lambda {\mathbf{\beta }}^T{\Lambda }^{(t-1)}\mathbf{\beta },$$

$$\underset{{\beta }_0,\mathbf{\alpha },\mathbf{w}}{\min }$$

$$0\le w_j<1\text{ for all }j=1,\dots ,p.$$

$$\underset{{\beta }_0,\mathbf{\alpha },\mathbf{w}}{\min }$$

$$w_j\ge 0\forall j=\{1,\dots ,n\}.$$

Full text

1 Introduction.

A considerable amount of recent biometric research is being conducted in the ‘personalized medicine’ framework, because it has been well accepted now that heterogeneity can exist among individual subjects’ response to treatment in many disease settings. The characteristics contributing to this heterogeneity may include patient demographics, genetic/genomic information or other biological markers, henceforth referred to as treatment-selection biomarkers. These biomarkers can be effectively used to select optimal therapies for individuals in order to optimize this clinical outcome. Also it is important to remember that a single biomarker may not sufficiently explain this heterogeneity, and multiple biomarkers may need to be combined to build the correct statistical framework to optimize the process of treatment selection.

The direct approach to identifying these optimal marker combinations in treatment selection involves parametric modeling of the disease risk conditional on biomarkers, treatment assignment and other baseline patient characteristics, and recommending treatment assignments based on whether the predicted risk under treatment is lower than the predicted risk under no treatment. This framework was first introduced by Song and Pepe (2004), and studied extensively in Foster et al. (2011); Qian and Murphy (2011); Lu et al. (2013). Treatment selection rules based on parametric risk models rely heavily on the correct specification of this model, which is often challenging given the complexity of biological mechanisms. An alternate, much more robust approach is to build optimization algorithms to minimize (or maximize) a desired criterion. This criterion, often called the objective function, is formulated based on relevant goals for treatment selection, and the optimization algorithm tries to find the best marker combination that optimizes it. These methods (also called indirect approaches) can be made completely nonparametric and assumption free and much more robust to model misspecifications. Zhang et al. (2012a, b) proposed finding the optimal marker combination within a pre-specified class by optimizing an estimator of the overall population mean outcome. Zhao et al. (2012) approached the same as an outcome-weighted learning problem and derived optimal treatment-selection rules using a weighted support vector machine method.

In Huang and Fong (2014), the authors proposed a new method to identify linear and nonlinear marker combinations by directly optimizing a targeted criterion function in the manner of Zhang et al. (2012a) and Zhao et al. (2012). However, the targeted criterion differed from others as adverse side-effects and/or cost incurred by the treatment were considered along with the event rates of the targeted disease in establishing the objective function. This was a crucial extension, as reducing safety events or cost of administering the treatment is often a valid policy goal (see Vickers et al., 2007). The authors formulated the optimization problem as minimization of a weighted sum of 0-1 loss, and used the ramp loss as an approximation of the 0-1 loss. In this article, we extend the idea of Huang and Fong (2014) in creating an augmented target criterion that takes on the additional challenge of controlling for the number of markers in the risk model as well. This is crucial for two vital reasons: (a) Measuring biomarkers can be expensive with respect to both time and money, and possibly invasive too. It is therefore of significant interest to limit the number of biomarkers that need to be collected for an individual to selecting his/her optimal treatment. (b) As discussed earlier, our original problem requires us to find optimal biomarker combinations to explain the disease response heterogeneity in individuals, but finding this optimal combination can be difficult in the presence of redundant markers that do not contribute to treatment selection. This may lead to overfitting and result in overly complex selection rules that have poor prediction performance. One way to deal with this ‘curse of dimensionality’ is through marker selection. Recently a lot of attention has been directed to effect modifier selection in precision medicine. For example, in Zhao (2017), the authors studied selective inference in effect modification models via LASSO; in Liang (2017), the authors constructed sparse decision rules in the context of concordance-assisted learning; and in Shi (2018), the authors proposed a penalized multi-stage A-learning algorithm for deriving the optimal dynamic treatment regime. Although these feature selection methods are all relevant to the broad genre of precision medicine, each of them target a very specific problem, which are quite different from our objective of penalized optimization of the single stage targeted objective function, where still only limited research exists: some insight into marker selection in treatment recommendation problems were briefly studied by Huang (2016) and Zhou et al. (2015), where both targeted minimization of disease rate in treatment selection.

To address these issues along with the original goals of treatment selection, we reformulate the criterion for treatment selection by associating a cost for measuring each biomarker in creating the optimal rule. The optimization problem can be written as minimization of a weighted sum of 0-1 loss, but with a $L_{0}$ penalty added for the number of markers in the model. In this article, we adopt the usual hinge loss convex relaxation of the 0-1 loss and perform a comprehensive investigation and comparison of various algorithms for solving this optimization problem. In linear support vector machines, $L_{0}$ feature extraction is a well-researched problem (see Bradley and Mangasarian, 1998; Weston et al., 2003; Mangasarian, 2006; Huang et al., 2010), although it is a much more challenging problem in nonlinear support vector machines (see Mangasarian and Wild, 2007). However, it is worth noting that our setting differs from a simple classification format in two vital aspects: (a) although the treatment selection objective can be rewritten into a (weighted) classification problem (as shown in Section 2), it is still in essence a fundamentally different problem from classification, and feature selection techniques in SVMs have not been studied under this context, and (b) weighted SVM is a more complicated optimization problem than the standard SVM, where the constraint on each support vector varies according to the weight associated with it, and research into feature extraction under this setting has also been fairly limited till now. Some work has been done to extend the SCAD penalty with the weighted linear support vector machines with special forms of such weights (see Jung, 2013), but beyond that, there hasn’t been any targeted investigation of such as per our knowledge. These two reasons make these explorations vitally important.

Thus, the biggest contribution of this article is to combine the methodological advances in two different areas of research, advances in indirect approaches to treatment selection through kernelized methods and those in feature selection techniques in nonparametric statistical learning methods like SVMs. We believe that combining these two areas of research is a novel methodological formulation in itself. Additionally each of these feature selection methods has been translated from the standard SVM framework to the weighted SVM formulation to match the proposed objective. Moreover, adding a penalty for the number of markers in the model redefines the objective function as well, which has been used in conjunction with these feature selection techniques - for example, in choosing tuning parameters for a given feature selection method, we propose the use of a generalized cross validation (GCV) technique that utilizes the whole objective criterion including the penalty for the number of markers.

The article continues in the following manner: In Section 2, we establish the problem of minimizing the total disease, treatment, and marker measurement cost in treatment-selection, and discuss marker selection in conjunction with treatment selection. We briefly discuss various methods built for support vector machines, adapted specifically for the weighted classification setup, for deriving the best linear and nonlinear marker combinations in order to optimize the desired objective function. In Section 3, we set up different simulation examples to test the strength of these linear and nonlinear feature selection methods, modified for our setting, and discuss the results. Then in Section 4, we illustrate the application of our methods using a real data example from an HIV vaccine trial. And finally, in Section 5, we discuss our findings regarding the use of these feature selection methods and the impact of incorporating marker measurement cost into treatment selection, and link to additional materials are presented in the Supplementary Section 6.

2 Methods.

In this article, we consider the problem of finding optimal treatment selection rules for a binary clinical outcome $Y$ (0 for nondiseased and 1 for diseased) based on a set of $p\geq 1$ candidate markers collected from an individual’s biological characteristics (the solution we propose is applicable to clinical outcome measured in continous scale as well). Denote this set of candidate markers by $\textbf{X}_{p}=\{X_{1},\dots,X_{p}\}$. For a given subset of markers $\textbf{X}\subseteq\textbf{X}_{p}$, we consider marker based treatment selection rules of the form $A(\textbf{X})=I(f(\textbf{X})>0)$, where $f\in\mathcal{F}_{\textbf{X}}$ and $\mathcal{F}_{\textbf{X}}$ denotes a class of functions spanning over X, and $I(\cdot)$ is the indicator function. The above then translates to the rule: $A=0$ for not treating, and $A=1$ for treating. The treatment-selection benefit of a decision rule $A(\textbf{X})$ can be quantified by $E_{A(\textbf{X})}(Y)$, the expected disease rate in the population as a result of treatment selection based on $A(\textbf{X})$ (Song and Pepe, 2004; Qian and Murphy, 2011). This measure has been widely accepted as a crucial metric in recent literature (Zhao et al., 2012; Zhang et al., 2012a). Although $E_{A(\textbf{X})}(Y)$ characterizes the burden of disease upon the population, one may be concerned with additional burden associated with a treatment regimen, such as its side effects and/or the monetary cost of its implementation. Hence it might be preferable to search for treatment-selection strategies that take these aspects into consideration as well. Huang and Fong (2014) proposed to incorporate additional burden associated with a treatment regimen, such as due to its side effects of monetary cost, by pre-specifying a treatment/disease harm ratio such that each burden type can be put on the same scale. Following the decision-theoretic framework of Vickers et al. (2007), let $\delta_{1}$ be a pre-specified ratio of the burden per treatment relative to the burden per disease event, and let $Y(1)$ and $Y(0)$ indicate the potential disease outcome if a subject were to receive or not receive the treatment. Then the total burden due to disease and treatment for $A(\textbf{X})$, represented in the unit of burden per disease event (see Huang and Fong, 2014) is given as $E_{A(\textbf{X})}(Y)+E\{\delta_{1}\times A(\textbf{X})\}$ $=\sum_{a=0}^{1}E\left[Y(a)\times I\{A(\textbf{X}=a)\}\right]+\delta_{1}\times E\{A(\textbf{X})\}$. And the best treatment-selection rule is derived as the one that minimizes this total burden.

In this article we further take into consideration the cost of measuring biomarkers in deriving the optimal treatment-selection rule, under the expectation that inclusion of the burden of measuring biomarkers in the criterion function can lead to the derivation of more cost-effective treatment-selection rules in the sense of achieving desired public health impact under the guidance of a parsimonious biomarker panel. We make the following assumptions: (i) measurement of each biomarker induces roughly equal cost, and (ii) the cost of measuring one biomarker is $\delta_{2}$ times the burden per disease event. Then the total burden due to disease, treatment, and biomarker measurement for a treatment-selection rule $A(\textbf{X})$ can be represented in the unit of the burden per disease event as

[TABLE]

Note: Although we assume a roughly equal cost for each biomarker for simplicity in this paper, it is straightforward to extend it to a setting where each biomarker has a different cost. Biomarkers which are absolutely essential to collect can be considered to have no cost at all, or there may be several groups of biomarkers, such that each group has a different cost depending on its burden and importance for the population of interest. In such situations, the penalty $\delta_{2}\times\dim(\textbf{X})$ can be replaced by $\sum_{i=1}^{\dim(\textbf{X})}\delta_{2,i}$ where $\delta_{2,i}\geq 0$ is the cost associated with the biomarker $X_{i}$, and our methods can still be applied as described with only minor modifications.

We propose to derive an optimal treatment-selection rule by minimizing this quantity $\theta$. Suppose there exists an ‘optimal’ or ‘correct’ set of $p_{0}$ biomarkers $\textbf{X}_{0}$ within the list of biomarkers measured, such that $\textbf{X}_{0}=\{X^{\ast}_{1},\dots,X^{\ast}_{p_{0}}\}\in\textbf{X}_{p}$, which when combined through an oracle rule $f_{0}$ minimizes $\theta$. Thus, we can write the above problem as,

[TABLE]

Note that under this setup, $f_{0}\in\mathcal{F}_{\textbf{X}_{0}}$. In practice, it is important to have a sensible way to specify the values or the range of values for $\delta_{1}$ and $\delta_{2}$. Choice of these cost ratios can be facilitated based on information of the monetary cost for controlling the targeted disease, for applying the treatment, and for biomarker measurement, as in our data example of making recommendation for HIV vaccine, presented later in Section 4.

Now imagine data from a two-arm randomized trial with the treatment indicator $T$ taking values 0 and 1, to refer being untreated and treated, respectively. Let $n_{0}$ and $n_{1}$ indicate the number of subjects in the untreated and treated arms, respectively. Thus we have i.i.d. samples of the form $\{Y_{i},\textbf{X}_{i},T_{i}\}$ for $i=1,\dots,n$ with $n=n_{0}+n_{1}$. As in Huang and Fong (2014), we assume: (i) stable unit treatment value (SUTVA) (Rubin, 1980) and consistency: $Y(0)$, $Y(1)$ of one subject is independent of the treatment assignments of other subjects, and given the treatment a subject actually received, a subject’s potential outcomes equal the observed outcomes; (ii) ignorable treatment assignments assumption: T\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}Y(0),Y(1)|\textbf{X}. Assumption (i) is plausible in trials where participants do not interact with one another and assumption (ii) is ensured by randomization. Under the above assumptions, it can be shown that,

[TABLE]

where $Risk_{0}(\textbf{X})=P(Y=1|\textbf{X},T=0)$ and $Risk_{1}(\textbf{X})=P(Y=1|\textbf{X},T=1)$ are the risk of $Y$ conditional on X among the untreated and treated, respectively.

Huang and Fong (2014) showed that when $\delta_{2}=0$, an optimal rule $A(\textbf{X})$ can be obtained as $A(\textbf{X})=1$ if $Risk_{0}(\textbf{X})-Risk_{1}(\textbf{X})>\delta_{1}$, and $A(\textbf{X})=0$ otherwise. But when $\delta_{2}$ is positive, such a strategy alone wouldn’t exactly work as we need to also control the number of markers in the model. Since the quantity $\delta_{2}$ imposes an $L_{0}$ penalty on the set of markers X, one ad hoc method may be to use standard regression-based methods with $L_{0}$/$L_{1}$ penalization to estimate the quantities $Risk_{0}(\textbf{X})$ and $Risk_{1}(\textbf{X})$ and then derive the optimal treatment-selection rule as $A(\textbf{X})=I\{Risk_{1}(\textbf{X})-Risk_{0}(\textbf{X})>\delta_{1}\}$. But note that the above procedure may lead to suboptimal rules with respect to our goal of minimizing $\theta$.

Alternatively, following the strategy in Zhang et al. (2012a), Zhao et al. (2012) and Huang and Fong (2014), we can consider a class of rules for treatment recommendation based on functionals of the form $f(\textbf{X})$ (belonging to some functional class $\mathcal{F}_{\textbf{X}}$ on X), and a given threshold (usually 0). In particular, we let $A(\textbf{X})=I\{f(\textbf{X})>0\}$ with $I(\cdot)$ the indicator function, $f(\textbf{X})=b+g(\textbf{X})$ with $g(\textbf{X})$ a function of markers X. Then, assuming randomization does not depend on X for simplicity, $\theta$ can be expressed as $\frac{E\left[Y\times T\times I\{f(\textbf{X})>0\}\right]}{P(T=1)}+\frac{E\left[Y\times(1-T)\times I\{f(\textbf{X})\leq 0\}\right]}{P(T=0)}+\delta_{1}\times P\{f(\textbf{X})>0\}+\delta_{2}\times\dim(\textbf{X})$. The optimal $f(\textbf{X})$ can be found by minimizing the empirical estimate of $\theta$, that is,

[TABLE]

Therefore, we can formulate this problem as the minimization of the sum of a weighted sum of 0-1 loss and a term proportional to the number of biomarkers. That is, $\hat{f}$ can be found as the minimizer of

[TABLE]

with the case-specific weight $W_{i}=W_{1i}=-\left\{\frac{Y_{i}\times T_{i}}{n_{1}/n}-\frac{Y_{i}\times(1-T_{i})}{n_{0}/n}+\delta_{1}\right\}$. Other types of weights such as control-specific or case-control mixture weights or their robust substitutes can also be adopted and are discussed in Huang and Fong (2014). For example, the robust substitute for the case-only weights $W^{rob}_{1i}=-\left\{\frac{Y_{i}\times T_{i}}{n_{1}/n}-\frac{Y_{i}\times(1-T_{i})}{n_{0}/n}+\frac{\pi-T_{i}}{\pi}\widehat{Risk}_{0}+\frac{\pi-T_{i}}{1-\pi}\widehat{Risk}_{1}+\delta_{1}\right\}$, where $\widehat{Risk_{0}}(\textbf{X})$ and $\widehat{Risk_{1}}(\textbf{X})$ are risk estimates obtained from a working model and $\pi=P(T=1)$. It is worthwhile to note that since the minimization of (2) is equivalent to the minimization of $E[I\{f(\textbf{X})\leq 0\}\times\{Risk_{0}(\textbf{X})-Risk_{1}(\textbf{X})-\delta_{1}\}]$, any consistent estimate of $Risk_{0}(\textbf{X})-Risk_{1}(\textbf{X})-\delta_{1}$ can also serve as weights in (4).

Note that if we are only interested in combining a set of candidate biomarkers without performing any variable selection, then it is not necessary to include the biomarker measurement cost in the criterion function as in Huang and Fong (2014). However, in the presence of a large number of biomarkers, incorporation of marker cost will impact the features selected in the estimated rule.

2.1 Treatment selection as a weighted Support Vector Machines problem.

In this section, we consider the minimization of the regularized loss function (4), conditional on a pre-specified set of weights $W$. Since $\sum_{i=1}^{n}W_{i}I\left\{f(\textbf{X}_{i})\leq 0\right\}/n\propto\sum_{i=1}^{n}|W_{i}|I\left\{\text{sgn}(f(\textbf{X}_{i}))\neq\text{sgn}(W_{i})\right\}/n$, (4) can be reformulated as

[TABLE]

Note that (5) is a weighted classification problem where $Y^{\ast}_{i}=\text{sgn}(W_{i})=\{-1,1\}$ is the true binary classes, $\text{sgn}\{f(\textbf{X}_{i})\}$ is the predicted binary class based on X, and $|W_{i}|$ is the subject specific weight. This type of problem can be resolved using the weighted support vector machine (Lin and Wang, 2002), based on $\textbf{X}_{i}$, $Y_{i}$, and $|W_{i}|$. Since minimization of a weighted sum of 0-1 loss is non-convex and intractable, we look to replace the 0-1 loss with a convex surrogates; the hinge loss is often used in this context. The hinge loss, given as $h(u)=\max(0,1-u)$ has been proven to be a useful surrogate to the classification loss, such that the term $|W_{i}|I\left\{\text{sgn}(f(\textbf{X}_{i}))\neq\text{sgn}(W_{i})\right\}$ in (5) is replaced by $|W_{i}|h\left(f(\textbf{X}_{i})\times\text{sgn}(W_{i})\right)$. It is worth noting that the hinge loss penalizes departure of a decision rule from its observed class label, based on the extent of this departure, which the classification loss fails to do. The hinge-loss SVM formulation of the above problem is given as,

[TABLE]

This is the $L_{0}$ penalized weighted support vector machines framework. In SVMs, the functional space $\mathcal{F}_{\textbf{X}}$ is generally restricted to be a reproducing kernel Hilbert space $\mathcal{H}_{\textbf{X}}(k_{\textbf{x}})$, represented uniquely by its kernel $k_{\textbf{x}}$. Before considering the solution to (6), we first review the weighted support vector machines formulation from Lin and Wang (2002), where the square of the Hilbert space norm $\|\cdot\|_{\mathcal{H}_{\textbf{X}}}$ is used instead of the $L_{0}$ norm. For a given subspace $\textbf{X}\in\textbf{X}^{p}$, we define it as:

[TABLE]

In linear support vector machines, $\mathcal{F}^{\text{lin}}_{\textbf{X}}=\{f_{\beta,b_{0}}:f_{\beta,b_{0}}(\textbf{x})=\left\langle\beta,\textbf{x}\right\rangle+b_{0},\>\textbf{x}\in\textbf{X},\>\beta\in\mathbb{R}^{\dim(\textbf{X})},\>b\in\mathbb{R}\}$ is used for optimization, an RKHS with the Euclidean inner product as its kernel function, $k_{\textbf{x}}(\textbf{x}_{1},\textbf{x}_{2})=\left\langle\textbf{x}_{1},\textbf{x}_{2}\right\rangle$. The Hilbert norm $\|f\|^{2}_{\mathcal{H}_{\textbf{X}}}$ in this case becomes the usual Euclidean $L_{2}$ norm on the linear combination weights $\|\beta\|^{2}$, and the estimated decision functions $\hat{f}$ can be expressed as a linear combination of the input marker set X. However, as many optimal marker combination for treatment selection may not be among linear combinations of the input markers, it is important to extend $\mathcal{F}_{\textbf{X}}$ to include more complex nonlinear functions. This can be achieved by considering transformations of the feature space through feature maps of the form $\phi(\textbf{x})$, the appropriate RKHS for which is the one with the kernel $k_{\textbf{x}}$ satisfying $k_{\textbf{x}}(\textbf{x}_{1},\textbf{x}_{2})=\left\langle\phi(\textbf{x}_{1}),\phi(\textbf{x}_{2})\right\rangle$. Examples of popular nonlinear kernels include the polynomial kernel of $d^{th}$ degree $k_{\textbf{x}}(\textbf{x}_{i},\textbf{x}_{j})=(1+\left\langle\textbf{x}_{i},\textbf{x}_{j}\right\rangle)^{d}$ and the radial basis function (RBF) kernel $k_{\textbf{x}}(\textbf{x}_{i},\textbf{x}_{j})=\exp\left(-\gamma\|\textbf{x}_{i}-\textbf{x}_{j}\|^{2}\right)$ with $\gamma$ as a tuning parameter. The resulting SVM solution at a given covariate vector $\textbf{x}_{0}$ is given as $f(\textbf{x}_{0})=\sum_{j=1}^{n}\alpha_{j}k(\textbf{x}_{0},\textbf{x}_{j})+\beta_{0}$, where $\alpha_{i}\in\mathbb{R}$’s are the trained weights on the support vectors $k(\cdot,\textbf{x}_{j})$ and $b_{0}\in\mathbb{R}$ is the estimated global constant.

2.2 Weighted versus unweighted support vector machines.

The weighted support vector machines is a well-known technique for classification. First proposed by Lin and Wang (2002), who called it the “fuzzy support vector machines” or FSVM, it has been further applied and studied in subsequent works (Fan and Ramamohanarao, 2005; Yang et al., 2007; Zhao et al., 2012, see for example). In essence, the weighted support vector machines becomes a very important tool in classification, when some training points are more important than others for the given problem. The main difference in solving the weighted SVM and the standard unweighted SVM lie in the constraints that we put on the support vectors in the dual formulation of the problem. To see this, let us write down the dual of the unweighted SVM below and note the changes that can transform it into a weighted SVM problem. Let us denote $Y^{\ast}_{i}=\text{sgn}(W_{i})$ and $q_{i}=|W_{i}|$, and then see that the unweighted squared Hilbert loss penalized SVM is given as $\min_{f\in\mathcal{H}_{\textbf{X}}}\sum_{i=1}^{n}\max(1-Y^{\ast}_{i}f(\textbf{X}_{i}),0)/n+\lambda\|f\|^{2}_{\mathcal{H}_{\textbf{X}}}$, the dual of which is given as,

[TABLE]

for $C=1/(2\lambda n)$. The dual of (7) differs from (8) only through the constraints that we put on the $\alpha_{i}$s. In the unweighted SVM, the upper-bound for each of $\alpha_{i}$ is the fixed constant $C$, while in the weighted SVM, the upper bound for each $\alpha_{i}$ is further multiplied by the quantity $q_{i}$ and thus becomes $C\times q_{i}$ (that is, $0\leq\alpha_{i}\leq Cq_{i}$). Intuitively it just means that $\alpha_{i}$, the coefficient associated with a given sample, is each constrained differently according to their weight.

2.3 Identifying optimal biomarker combinations in the linear space.

As noted before, variable selection is a key motivation in our pursuit to solve (6), which isn’t an inherent feature of the support vector machines framework of (7). The Hilbert norm provides some control on the overcomplexification of the estimated functions, but without performing any variable selection. Many authors have proposed to use the $L_{1}$ and $L_{0}$ penalty either as a replacement or in conjunction with the $L_{2}$ norm in the linear SVMs. First of all, note that solving (6) is combinatorially a very hard problem (see Amaldi and Kann, 1998), but the zero norm is directly related to finding minimal subsets and can provide optimal feature extraction if properly implemented. There are however many feature extraction algorithms for the unweighted linear support vector machines that do not rely on the $L_{0}$ norm for selection (see for example Mangasarian, 2006; Zhang et al., 2006). To this effect, in this article, we also look at alternate ways to perform feature selection in support vector machines, instead of focusing solely on solving (6). Here, we build our optimization approaches on some of the most commonly used algorithms for unweighted SVMs. That is, we modify each of them to cater to the weighted version of the algorithm, and we will see how it paves the way for directly or indirectly solving (6) in the linear space. We now briefly introduce each method below, along with their modified objective function under the weighted SVM setup, while a more detailed description of each is reserved for the Supplementary Section (see Web Appendix LABEL:sec:appA).

$L_{1}$ weighted support vector machines ($L_{1}$ WSVM), that utilizes the $L_{1}$ norm for penalization instead of the squared Hilbert norm, built upon Mangasarian (2006), with objective function,

[TABLE] 2. 2.

SCAD applies the smoothly clipped absolute deviation penalty (Zhang et al., 2006) to the WSVM objective function ,

[TABLE]

$p_{\lambda}(|\boldsymbol{\beta}|)=\lambda|\boldsymbol{\beta}|1(|\boldsymbol{\beta}|\leq\lambda)-\frac{|\boldsymbol{\beta}|^{2}-2a\lambda|\boldsymbol{\beta}|+\lambda^{2}}{2(a-1)}1(\lambda<|\boldsymbol{\beta}|\leq a\lambda)+\frac{(a+1)\lambda^{2}}{2}1(|\boldsymbol{\beta}|>a\lambda)$. 3. 3.

elastic SCAD, that uses a mixture of the SCAD penalty and the $L_{2}$ norm,

[TABLE] 4. 4.

feature selection concave (FSV), built on a concave minimization algorithm originally proposed by Bradley and Mangasarian (1998), achieves approximate $L_{0}$ penalization in weighted SVMs through solving,

[TABLE] 5. 5.

approximation of the zero norm minimization (AROM), that achieves approximate $L_{0}$ penalization (see Weston et al., 2003) by solving the following objective function iteratively. It initializes a vector $\boldsymbol{z}=(1,\dots,1)$, and at each successive step, it resets $\boldsymbol{z}$ as $\boldsymbol{z}\times\boldsymbol{\beta}$,

[TABLE]

where $\boldsymbol{x}\ast\boldsymbol{w}=\left(x_{1}w_{1},\dots,x_{p}w_{p}\right)$. 6. 6.

$L_{0}$ weighted support vector machines ($L_{0}$ WSVM), built on the works of Huang et al. (2010), that achieves $L_{0}$ norm in weighted support vector machines through an iterative scheme, by solving the following objective function at the $t^{th}$ step,

[TABLE]

subject to $Y^{\ast}_{i}(\left\langle X_{i}\boldsymbol{\beta}\right\rangle+\beta_{0})\geq 1-\xi_{i},\>\xi_{i}\geq 0,\>i=1,\dots,n$, where $\Lambda^{(t-1)}=diag(1/|\beta^{(t-1)}_{1}|^{2},\dots,1/|\beta^{(t-1)}_{p}|^{2})$.

2.4 Identifying optimal biomarker combinations in the nonlinear space.

In this section, we consider the derivation of best biomarker combinations in the nonlinear space, as most optimal marker combinations for treatment-selection rules may not be among the linear combinations of the input markers. Feature selection is a fairly straightforward procedure for linear SVM classifiers, but it is a much more challenging problem in nonlinear support vector machines (see Mangasarian and Wild, 2007). For example, in the linear support vector machines, we can use the $L_{1}$ or $L_{0}$ norm instead of the Hilbert norm. However, when similar techniques are applied to nonlinear SVM classifiers, we have reduction in the number of support vectors but not in the number of input space features (Fung and Mangasarian, 2004). So although the dimensionality of the transformed space is reduced, it does not provide any direct reduction of input space features. A few procedures have recently been developed however to deal with feature selection in the original input space under nonlinear feature maps (see Dasgupta et al., 2013; Allen, 2013), but not a lot of them have been developed with the goal of $L_{0}$ penalization. We will investigate some of these newly developed methods here, but under the more generalized setting of weighted support vector machines. We briefly introduce each method below, along with their modified objective function under the weighted SVM setup, while a more detailed description of each is reserved for the Supplementary Section (see Web Appendix LABEL:sec:appB).

risk recursive feature elimination (riskRFE), a newly developed, powerful wrapper technique based on recursive computation of the learning function (see Dasgupta et al., 2013), that works both in the linear and the nonlinear space (see Figure 1 for a flow chart of the algorithm). 2. 2.

kernel iterative feature extraction (KNIFE), that achieves feature selection in nonlinear SVMs, by optimizing a feature-regularized loss function (by weighting features within the kernel), iteratively (see Allen, 2013). For $k_{\boldsymbol{w}}\left(\textbf{x}_{1},\textbf{x}_{2}\right)=k\left(\boldsymbol{w}\ast\textbf{x}_{1},\boldsymbol{w}\ast\textbf{x}_{2}\right)$, the weighted version of the algorithm solves,

[TABLE] 3. 3.

kernel penalized weighted support vector machines (KP-WSVM), built on the works of Maldonado et al. (2011) for feature selection in nonlinear SVMs, that relies on penalizing each feature’s use in the dual formulation.

[TABLE]

The approximated $L_{0}$ penalty in the above objective function can also be replaced by a $L_{1}$ norm.

3 A simulation study

In this section, we investigate the performance of the aforementioned feature selection methods in conjunction with the weighted support vector machines algorithms in minimizing the $L_{0}$ penalized sum of the weighted classification loss as in (5) in various simulation settings that we describe below.

3.1 The simulation setup

We consider data from 1:1 randomized trial of size $n=500$ for studying the performance of our proposed strategy. The linear feature selection methods (riskRFE (L), $L_{1}$ WSVM, SCAD/eSCAD, FSV, AROM, $L_{0}$ WSVM) are used in conjunction with the linear kernel based weighted SVM, while the nonlinear marker selection methods (riskRFE (G), KNIFE, $L_{0}$ KP-WSVM, $L_{1}$ KP-WSVM) are used in conjunction with the Gaussian RBF kernel based weighted SVM. For each simulation setting, we restrict ourselves to when the treatment/disease burden ratio $\delta_{1}$ is equal to 0, but consider three different values for the marker cost $\delta_{2}$, which is either $=0/0.0001/0.001$ units. The feature selection methods considered here can be broadly categorized into two groups - (a) Embedded (penalized) methods, which have a separate tuning parameter for penalizing the $L_{1}$ or $L_{0}$ cost, including SCAD, eSCAD, $L_{1}$ WSVM, $L_{0}$ WSVM, FSV, KNIFE and KP-WSVMs. For these methods, $\delta_{2}$, the ratio of the cost of measuring one biomarker to the burden per disease event, is utilized directly for choosing the optimum value of the respective tuning parameter(s). Although AROM also belongs to this class, it achieves $L_{0}$ penalization through an iterative fitting of the $L_{2}$ norm SVM. Thus in this case, $\delta_{2}$ is used for choosing the optimum amount of $L_{2}$ tuning; and (b) Wrappers that have no separate penalization for the $L_{1}$ or $L_{0}$ cost (both riskRFE(L) and riskRFE(G)); thus are not tuned on the value of $\delta_{2}$. The other global SVM parameters, such as the width of the Gaussian kernel in nonlinear SVMs, are tuned based on the prediction performance of the weighted support vector machines in the full model. In this simulation study, we only consider the double-robustness weights of Huang and Fong (2014) for obtaining $W_{i}$ for each individual. As a comparative method, we use linear logistic regression model with LASSO and weighted SVM without feature selection to find the optimal treatment selection rule. In logistic regression, the $P(Y=1|\textbf{X},T)$ is modeled as a function of the main effects of the treatment and the biomarkers, and the interaction effects between the treatment and the biomarkers, with the amount of $L_{1}$ penalization tuned globally through usual cross validation. This is also the risk model used for constructing the double-robustness weights.

In each Monte Carlo simulation, we follow a 5-fold cross-validation procedure to identify the optimal tuning parameters in the model for a given method. As mentioned before, we perform this tuning in a two-step procedure - the global SVM parameters are tuned first using the full model in an weighted SVM analysis. The Hilbert Space norm of the estimated function (which in linear kernel becomes the $L_{2}$ norm) is tuned on a grid of values lying between $(0.0001,100)$, while the width of the Gaussian kernel, $\gamma$, is tuned over a grid of values lying between $(0.001,10)$. Then for each of the embedded methods, we use cross validation to tune for the additional penalization parameter(s) based on the performance of the method in question, using the selected global SVM parameters from step one. Each of these parameters was tuned on a grid of a viable range of values, and the value obtaining the optimal GCV performance is chosen for that method. This two-step procedure is followed to gain computation time, and also because the $L_{2}$ norm allows shrinkage of the estimated effects, rather than controlling sparsity directly.

The optimal treatment-selection rule $\hat{A}(\hat{\textbf{X}}_{0})$, where $\hat{\textbf{X}}_{0}\subseteq\textbf{X}_{p}$ is the optimal set of markers chosen by a given method, is then estimated from the training sample, after retuning the global SVM parameters in the model with the selected markers. To evaluate the performance of the estimated rule, a test set of $n=5000$ is generated in each simulation run, based on which we estimate $\theta$ as $\hat{\theta}=\sum_{i=1}^{n}\left\{(1-T_{i})\times Y_{i}\times(1-\hat{A}(\hat{\textbf{X}}_{0i}))/n_{0}+T_{i}\times Y_{i}\times\hat{A}(\hat{\textbf{X}}_{0i})/n_{1}+\delta_{1}\hat{A}(\hat{\textbf{X}}_{0i})\right\}+\delta_{2}|\hat{\textbf{X}}_{0}|$. We evaluate the performance of each method over 100 Monte Carlo runs. In our results, we present the estimated version of the quantity $\theta$ that we evaluate from these Monte Carlo runs. We compare performance of different methods in the following three settings.

In the first setting, we have in total a list of 27 markers, of which only 2 are significant in explaining treatment-marker interactions ($|\textbf{X}_{p}|=27,\>|\textbf{X}_{0}|=2$). The significant markers, $(X_{1},X_{2})$ are generated from the multivariate normal distribution \sim N\left(\left(\begin{array}[]{c}0\\ 0\end{array}\right),\left(\begin{array}[]{cc}1&0.2\\ 0.2&1\end{array}\right)\right). Each of the rest is generated independently from a $N(0,1)$ distribution. We consider two subcases under this setting: (i) a ‘linear’ underlying model, given by $logitP(Y=1|X_{1},X_{2},T)=-1.5-1.5X_{1}-1.25X_{2}+2X_{1}T+1.5X_{2}T$, with disease prevalence approximately 0.3 and 0.2 among the untreated and treated, respectively; and (ii) a polynomial ‘nonlinear’ underlying model $logitP(Y=1|X_{1},X_{2},T)=-1.5+0.2X_{1}-0.2X_{2}-3T-X_{1}T-X_{2}T+X_{1}^{2}T+X_{2}^{2}T$, with disease prevalence approximately 0.18 and 0.17 among the untreated and treated, respectively.

In the second setting, we have in total a list of 53 markers, of which only 3 are significant in explaining treatment-marker interactions ($|\textbf{X}_{p}|=53,\>|\textbf{X}_{0}|=3$). The significant markers, $(X_{1},X_{2},X_{3})$ are generated from the multivariate normal distribution \sim N\left(\left(\begin{array}[]{c}0\\ 0\\ 0\end{array}\right),\left(\begin{array}[]{ccc}1&0.2&0.2\\ 0.2&1&0.2\\ 0.2&0.2&1\end{array}\right)\right). Each of the rest is generated independently from a $N(0,1)$distribution with mean 0 and standard deviation 1. Again we consider two subcases: (i) a ‘linear’ underlying model, given by $logitP(Y=1|X_{1},X_{2},X_{3},T)=-2+X_{1}+0.75X_{2}+X_{3}-2X_{1}T-1.5X_{2}T-2X_{3}T$, with disease prevalence approximately 0.22 among both the untreated and treated groups; and (ii) a ‘nonlinear’ underlying model $logitP(Y=1|X_{1},X_{2},X_{3},T)=-0.8-3T-X_{1}T-X_{2}T+X_{1}^{2}T+X_{2}^{2}T-\sqrt{|X_{3}|}T$, with disease prevalence approximately 0.30 and 0.18 among the untreated and treated, respectively.

In the third setting, we consider a complex nonlinear setup, with a total of 27 markers, of which only 2 are significant in explaining treatment-marker interactions. All of the markers are generated from the uniform distribution, that is, $X_{i}\sim U(-2,2),\>i=1,\dots,27$. The model for treatment marker interaction considers a situation where the treatment has benevolent effect for an individual only if his/her marker values $X_{1}$ and $X_{2}$ lie in a specific region of the covariate space, without which the treatment might yield a harmful result. This is achieved by dividing the covariate space spanned by $X_{1}$ and $X_{2}$, the square $\{-2\leq X_{1}\leq 2,-2\leq X_{2}\leq 2\}$, into two regions: (i) a ‘harmful’ zone given by the smaller square $\{-1\leq X_{1}\leq 1,-1\leq X_{2}\leq 1\}$, denoted as $H$; and (ii) a ‘beneficial’ zone given by the region between the two concentric squares, denoted as $B$ (Figure 2 plots $X_{1}$ and $X_{2}$ marker values for 10000 such hypothetical individuals). Thus, the probability that the treatment is benevolent for the population is 0.75. The probability of disease for an individual in region $B$ is 0 if the individual receives treatment, but 0.6 otherwise, and they are reversed for an individual from $H$. The global probability for disease in the population is thus fixed at 0.3. The disease prevalences are approximately around 0.53 and 0.17 among the untreated and treated, respectively.

3.2 Results

The results for the simulation exercise are summarized in Tables 1-5. We give an overview of these results below, while a more detailed discussion of each of these settings is provided in the Supplementary Materials (see Web Appendix LABEL:sec:appC).

While most biomarker-based treatment-selection rules result in reduction of the total cost compared to the optimal strategy between treating all or treating none, methods that allow for feature selection lead to substantial improvement compared to the weighted SVM method without feature selection. 2. 2.

Embedded feature selection methods that involve tuning by marker measure cost $\delta_{2}$ in general show an decreasing trend in sensitivity and an increasing trend in specificity with increasing values of $\delta_{2}$, which is not the case for wrapper methods (such as riskRFE) or LASSO, that do not rely on $\delta_{2}$ for tuning. As a result, for the latter type of methods, the total cost $E(\theta)$ necessarily increases with increasing $\delta_{2}$, the cost of measuring a marker; in contrast, for methods that involve $\delta_{2}$ in tuning, the total cost $E(\theta)$ can often decrease with increasing $\delta_{2}$, especially when the improvement in the specificity score for a given method is more substantial than the decline in its sensitivity score. 3. 3.

In presence of a linear trend, best performing methods are typically among feature selection methods with linear kernel, yet feature selection methods with nonlinear kernel can have comparable performance; in presence of nonlinear trends, feature selection methods with nonlinear kernel can have substantial improvement over linear methods. 4. 4.

Relative performance of various methods vary with settings. In general, riskRFE performs really well when the cost of marker measurement is low. Apart from it, $L_{1}$ WSVM, SCAD and AROM are the best performing linear methods, especially when the cost of marker measurement is high, while $L_{1}$ KP-WSVM stands out among the nonlinear feature selection methods, with robust performance across settings.

We also evaluate these algorithms in a few additional settings. We compare their performance against that of the Decision List method of Zhang et al. (2015). We also consider a setting where the true optimal decision rule is not sparse but the effect of some biomarkers on the optimal treatment decision are so small that given the cost consideration, a sparser decision rule is more optimal in terms of the criterion $\theta$. For space constraint, these additional settings are discussed in detail in the Supplementary Materials (Web Appendix LABEL:sec:appD), and the results are provided in Supplementary Tables LABEL:tab:set4_decision, LABEL:tab:set4_linear and LABEL:tab:set5_linear.

4 Real data analysis

We now use an example from the RV144 Thailand HIV vaccine trial to examine the performance of the methods discussed above for selecting and combining markers. RV144 is one of the first vaccine trials that showed a significant positive effect of vaccine in preventing HIV infection. It included 16,402 participants, aged between 18 and 30 years, randomized 1:1 between vaccine and placebo (Rerks-Ngarm et al., 2009). A followup hostgenetic study was conducted to measure the effect of genotypes of Fc receptor genes on vaccine efficacy, and to that effect 190 single nucleotide polymorphisms (SNPs) (including five Fc-$\gamma$ and one Fc-$\alpha$ receptors) were genotyped on 125 cases (74 placebo recipients and 51 vaccine recipients), and 225 controls (20 placebo recipients and 205 vaccine recipients), of which 28 SNPs were selected based on Hardy-Weinberg equilibrium by Li et al. (2014), each categorized into a binary variable, to study association of each with vaccine efficacy. Here we consider all 28 SNPs as the expensive candidate biomarkers to select from. In addition, age, gender, and baseline behaviorial risk for vaccine are combined with SNPs for treatment recommendation as in Huang (2016); no penalty is put on those baseline variables as they are readily available from the trial. Lifetime HIV treatment cost is estimated to be around $370K according to CDC (https://www.cdc.gov/hiv/programresources/guidance/costeffectiveness/index.html). Vaccines typically cost a few hundred dollars and have minimal side effects, so we set$\delta_{1}=0.001$. The average cost of a SNP evaluation varies between$ 2-4, and hence we consider a range of values for $\delta_{2}\in\{10^{-6},5\times 10^{-6},10^{-5}\}$ based on the burden of measuring a SNP with respect to the treatment cost (as $2/370000\approx 5\times 10^{-6}$). We compare the performance of the linear and nonlinear feature selection methods to select the optimal subset of markers for the purposes of vaccine recommendation for an individual. We also use LASSO as a method of selecting markers and estimating the treatment recommendation rule. A fivefold CV is performed to select $\lambda$, the parameter which controls for overfitting in the estimated WSVM decision function, and we do the same to select the width of the Gaussian kernel, whenever it was used. We also employed a fivefold CV to choose the other tuning parameters, specific to each of the feature selection methods employed. To compute the expected disease rate for marker selection for each of these methods, we perform a cross validation procedure by splitting the data into five random folds, and using four folds for training the WSVM (or LASSO) procedures with the selected markers, and using the remaining one for testing. The procedure is repeated 100 times and the average disease rate is computed.

Table 6 shows the estimated performance of different selection strategies along with the strategy of treating none and treating all. Those two lead to an estimated HIV infection rate of 8.54 and 6.44 per 1000 persons, respectively, consistent with the positive vaccine efficacy we observed in the RV144 trial. From Table 6, it is clear that linear weighted support vector machines (with or without marker selection) does a better job of treatment recommendation compared to the nonlinear weighted support vector machines (with or without marker selection) with the Gaussian kernel in this particular example. It can be seen that both the linear and nonlinear WSVM without any selection perform at par with the strategy of treating all, even when we consider a biomarker cost, when it is relatively low, but gets worse with higher values of $\delta_{2}$. Most of the linear selection methods yield a total cost at least as good as the strategy of treating all, even when $\delta_{2}$ is high. Among the linear methods, only elastic SCAD achieves a total cost worse than the strategy of treating all when $\delta_{2}$ is high. On the other hand, AROM is outperforming every other strategy, followed closely by riskRFE (L). Overall we can conclude that marker selection is an important formulation for treatment recommendation using the weighted support vector machines framework.

5 Concluding remarks

In this article, we developed a new framework to incorporate marker measurement cost into treatment regime identification, with a pre-determined cost for inclusion of each biomarker in the model. We extended several different marker selection methods to apply to the weighted support vector machines setting, encompassing both the linear and the nonlinear space, in order to derive the optimal treatment-selection rules that minimize the total cost due to disease, treatment, and marker measurements. We investigated their performance in a number of different setups through a detailed simulation study, and also in a real data scenario, the RV144 HIV vaccine trial. We showed that in presence of a large number of candidate biomarkers, marker selection is essential for deriving cost-effective treatment-selection rules that effectively reduces disease and treatment burdens to the population while avoiding the burdens of collecting information on irrelevant biomarkers. We showed that marker selection also reduced the chance of obtaining overfitted treatment recommendation rules, which can result in lower test disease rate predictions than when we use the full model for treatment recommendation. It is worthwhile to note here that the indirect approaches can sometimes be inefficient compared with direct approaches, especially when the underlying disease risk model is correctly specified in the direct approaches. Moreover, inference for the estimated optimal regime can become challenging under indirect approaches. However, these methods are still appealing and provide a complementary alternative to direct approaches because of their robustness to model misspecifications. In fact, comparison of our algorithms with logistic regression with LASSO showed that in some settings, even when the risk model assumptions hold (as was the case in most of our simulation settings), feature selection using the indirect approaches can lead to treatment-selection rules with performance comparable to or even better than that based on direct approaches. Among the various marker selection technique employed in the simulations, riskRFE stood out as the best performing method when cost of measuring marker measurement cost was low, while AROM, SCAD and $L_{0}$ WSVM were the best performing linear feature selection methods for higher costs of marker selection, depending on the setting. On the other hand, $L_{1}$ KP-WSVM was clearly the best performing nonlinear feature selection method in most of the settings. In general, we showed that using nonlinear feature selection methods can perform substantially better than linear feature selection methods in presence of nonlinear patterns, at a minimal cost of efficiency in presence of linear patterns. One last thing to note here is that the methodology put forward in this article typically invoke the randomized setting. However, if the assumptions of (i) SUTVA and (ii) ignorable treatment assignment (see Page 5 in Section 2 for more details) hold in an observational setting too, our proposed methods will continue to be applicable there, although assumption (ii) is typically not verifiable without randomization.

6 Supplementary Materials

Additional materials referenced in Sections 2 and 3 are available with this paper at the Journal of Royal Statistical Society Series C web page hosted at the Wiley Online Library.

7 Acknowledgments

The work was supported by NIH grant R01 GM106177-01. The authors thank the participants, investigators, and sponsors of the RV144 trials. The authors thank Dr. Dan Geraghty from the Fred Hutchinson Cancer Center for generating the genetics data, and Drs. Sue Li and Peter Gilbert for pre-processing and preliminary analysis of the SNP data. The views expressed are those of the authors and should not be construed to represent the positions of the U.S. Army or the Department of Defense.