Feature Screening in Large Scale Cluster Analysis

TL;DR

This paper introduces a fast feature screening method for large-scale clustering, effectively discarding non-informative features in massive datasets using a fusion penalization approach, supported by theoretical bounds and empirical validation.

Contribution

The paper presents a novel, efficient feature screening technique for clustering large datasets, with theoretical guarantees and demonstrated effectiveness in high-dimensional applications.

Findings

High probability of perfect screening of noise features

Competitive performance in simulation studies

Effective application to single-cell gene expression data

Abstract

We propose a novel methodology for feature screening in clustering massive datasets, in which both the number of features and the number of observations can potentially be very large. Taking advantage of a fusion penalization based convex clustering criterion, we propose a very fast screening procedure that efficiently discards non-informative features by first computing a clustering score corresponding to the clustering tree constructed for each feature, and then thresholding the resulting values. We provide theoretical support for our approach by establishing uniform non-asymptotic bounds on the clustering scores of the "noise" features. These bounds imply perfect screening of non-informative features with high probability and are derived via careful analysis of the empirical processes corresponding to the clustering trees that are constructed for each of the features by the associated clustering procedure. Through extensive simulation experiments we compare the performance of our proposed method with other screening approaches, popularly used in cluster analysis, and obtain encouraging results. We demonstrate empirically that our method is applicable to cluster analysis of big datasets arising in single-cell gene expression studies.

Tables

Table 1. Table 1: False Negatives and False Positive rates for Simulation Experiment I. Here, ℐ S = { 1 , … , 5 } subscript ℐ 𝑆 1 … 5 \mathcal{I}_{S}=\{1,\ldots,5\} , p = 50 𝑝 50 p=50 . The numbers in parenthesis are standard errors over 50 50 50 repetitions.

		$n=200$		$n=1000$		$n=2500$
		Avg FN	Avg FP	Avg FN	Avg FP	Avg FN	Avg FP
	0.05	0.06 (0.03)	42.38 (0.22)	0.10 (0.04)	21.26 (0.53)	0.14 (0.06)	6.08 (0.33)
COSCI	0.08	0.22 (0.06)	32.88 (0.35)	0.24 (0.06)	10.24 (0.41)	0.34 (0.08)	2.22 (0.23)
with	0.1	0.34 (0.07)	26.96 (0.43)	0.34 (0.07)	7.14 (0.33)	0.40 (0.08)	1.20 (0.17)
${\alpha }_0$	0.12	0.52 (0.08)	22.12 (0.41)	0.50 (0.1)	5.16 (0.27)	0.48 (0.09)	0.82 (0.11)
fixed	0.15	0.74 (0.1)	16.16 (0.48)	0.76 (0.1)	3.28 (0.19)	0.64 (0.09)	0.50 (0.10)
	0.2	1.12 (0.1)	9.36 (0.40)	1.04 (0.1)	1.68 (0.15)	0.92 (0.09)	0.30 (0.08)
	Data driven	1.8 (0.06)	0.54 (0.13)	0.88 (0.1)	2.14 (0.17)	0.36 (0.08)	1.28 (0.10)
	SpKM	2.96 (0.24)	16.04 (2.28)	3.26 (0.26)	13.28 (2.27)	3.9 (0.41)	7.00 (4.19)
Other	SpHC	3.34 (0.21)	16.86 (1.15)	2.08 (0.19)	22.54 (0.77)	1.6 (0.40)	24.5 (2.02)
methods	SAS	1.46 (0.14)	19.66 (1.15)	1.70 (0.18)	19.52 (1.29)	1.86 (0.20)	17.04 (1.25)
	Ex. Mass	2.3 (0.09)	0.00 (0.00)	2.00 (0.00)	0.00 (0.00)	2.00 (0.00)	0.00 (0.00)
	IF-PCA	0.32 (0.07)	17.64 (0.93)	0.10 (0.04)	18.08 (0.86)	0.10 (0.04)	17.24 (0.99)

Table 2. Table 2: False Negatives and False Positive rates for Simulation Experiment II. Here, ℐ S = { 1 , … , 6 } subscript ℐ 𝑆 1 … 6 \mathcal{I}_{S}=\{1,\ldots,6\} , p = 100 𝑝 100 p=100 . The numbers in parenthesis are standard errors over 50 50 50 repetitions.

		$n=200$		$n=1000$		$n=2500$
		Avg FN	Avg FP	Avg FN	Avg FP	Avg FN	Avg FP
	0.05	0.06 (0.03)	84.76 (0.41)	0.1 (0.04)	34.62 (0.71)	0.14 (0.06)	8.2 (0.44)
COSCI	0.08	0.22 (0.06)	61.9 (0.57)	0.38 (0.08)	15.18 (0.54)	0.58 (0.11)	2.64 (0.25)
with	0.1	0.38 (0.08)	49.5 (0.64)	0.52 (0.09)	10.3 (0.40)	0.76 (0.12)	1.38 (0.18)
${\alpha }_0$	0.12	0.6 (0.09)	39.32 (0.62)	0.72 (0.12)	7.34 (0.30)	0.96 (0.11)	0.96 (0.13)
fixed	0.15	0.88 (0.1)	28.16 (0.62)	1.06 (0.13)	4.86 (0.26)	1.36 (0.12)	0.58 (0.11)
	0.2	1.4 (0.11)	15.9 (0.52)	1.74 (0.12)	2.46 (0.19)	1.86 (0.10)	0.34 (0.09)
	Data driven	1.66 (0.08)	2.04 (0.28)	0.86 (0.11)	6.98 (0.24)	0.18 (0.06)	7.04 (0.27)
	SpKM	4.36 (0.25)	29.04 (3.23)	5.02 (0.26)	18.42 (3.93)	5.5 (0.50)	14.00 (8.56)
Other	SpHC	3.52 (0.27)	38.6 (3.05)	1.54 (0.23)	61.9 (3.18)	—	—
methods	SAS	1.14 (0.15)	58.48 (1.91)	2.1 (0.28)	41.66 (2.95)	1.82 (0.27)	43.9 (2.95)
	Ex. Mass	2.36 (0.10)	0.00 (0.00)	2.00 (0.00)	0.00 (0.00)	2.00 (0.00)	0.00 (0.00)
	IF-PCA	1.96 (0.15)	18.46 (2.24)	0.3 (0.06)	38.84 (0.61)	0.14 (0.05)	42.36 (0.23)

Table 3. Table 3: False Negatives and False Positive rates for Simulation Experiment III. Here, ℐ S = { 1 , … , 7 } subscript ℐ 𝑆 1 … 7 \mathcal{I}_{S}=\{1,\ldots,7\} , p = 5000 𝑝 5000 p=5000 . The numbers in parenthesis are standard errors over 50 50 50 repetitions.

		$n=200$		$n=1000$		$n=2500$
		Avg FN	Avg FP	Avg FN	Avg FP	Avg FN	Avg FP
	0.05	0.12 (0.05)	4341.6 (3.13)	0.82 (0.08)	1,355.3 (3.84)	1.06 (0.07)	274.62 (2.56)
COSCI	0.08	0.48 (0.08)	2979.62 (5.06)	1.26 (0.09)	548.8 (3.11)	1.58 (0.11)	82.86 (1.30)
with	0.1	0.84 (0.11)	2279.44 (4.37)	1.44 (0.1)	349.68 (2.04)	1.76 (0.12)	47.5 (1.08)
${\alpha }_0$	0.12	1.22 (0.12)	1764.22 (5.09)	1.68 (0.12)	242.7 (1.75)	1.96 (0.11)	30.8 (0.83)
fixed	0.15	1.56 (0.13)	1215.96 (4.86)	2.04 (0.13)	151.74 (1.21)	2.36 (0.12)	17.64 (0.62)
	0.2	2.22 (0.12)	636.62 (3.68)	2.72 (0.13)	69.58 (0.96)	2.86 (0.10)	7.74 (0.44)
	Data driven	2.68 (0.12)	297.14 (6.50)	1.32 (0.10)	498.76 (2.86)	0.92 (0.07)	516.28 (4.19)
Other	Ex. Mass	4.86 (0.14)	0.00 (0.00)	4.00 (0.00)	0.00 (0.00)	4.00 (0.00)	0.00 (0.00)
methods	IF-PCA	5.68 (0.09)	1802.4 (4.29)	6.00 (0.00)	1943.6 (1.86)	6.00 (0.00)	1970.00 (0.88)

Table 4. Table 4: False Negatives and False Positive rates for Simulation Experiment IV. Here, ℐ S = { 1 , … , 9 } subscript ℐ 𝑆 1 … 9 \mathcal{I}_{S}=\{1,\ldots,9\} , p = 25,000 𝑝 25,000 p=\mbox{25,000} . The numbers in parenthesis are standard errors over 10 10 10 repetitions.

		$n=200$		$n=1000$		$n=2500$
		Avg FN	Avg FP	Avg FN	Avg FP	Avg FN	Avg FP
	0.05	0.12 (0.05)	18587.2 (8.03)	0.6 (0.22)	5310.00 (22.81)	1.00 (0.1)	1101.67 (11.5)
COSCI	0.08	0.48 (0.08)	12176.00 (10.6)	1.30 (0.3)	2170.00 (11.71)	1.47 (0.192)	339.2 (4.82)
with	0.1	0.84 (0.11)	9221.48 (10.01)	1.50 (0.27)	1397.5 (9.96)	1.6 (0.19)	196 (2.68)
${\alpha }_0$	0.12	1.22 (0.12)	7088.34 (9.18)	1.70 (0.26)	976.5 (8.16)	1.73 (0.18)	126 (2.99)
fixed	0.15	1.60 (0.13)	4862.02 (8.61)	2.10 (0.28)	616.40 (7.22)	2.00 (0.22)	72.87 (2.41)
	0.2	2.42 (0.12)	2537.76 (6.92)	2.60 (0.27)	279.5 (3.27)	2.8 (0.22)	29.07 (1.65)
	Data driven	2.84 (0.14)	1687.34 (17.77)	1.40 (0.30)	1928.3 (18.72)	1.00 (0.1)	1324.8 (19.27)
Other	Ex. Mass	7.16 (0.13)	0.00 (0.00)	6.00 (0.00)	0.00 (0.00)	5.2 (0.11)	0.00 (0.00)
methods	IF-PCA	9.00 (0.00)	4061.2 (17.71)	9.00 (0.00)	5286.5 (17.64)	9.00 (0.00)	5563.3 (12.97)

Table 5. Table 5: This table has \textcolor blackthree parts. In the first part of the table, we provide some preliminary information on the \textcolor blackthree datasets. The second part of the table has CER for the \textcolor blackthree competing methods. For IF-PCA, we also report the number of selected features after the ‘/’ symbol. In the last part of the table, CER for COSCI + “clustering method” is reported \textcolor blackalong with the number of selected features. The numbers in parenthesis are standard errors over 30 repetitions of k 𝑘 k -means and IF-PCA. The standard errors are only reported when they exceed 0.0005. The SAS algorithm faced \textcolor blackscalability issues (marked by dash) for RNASeq data

		Multi-tissue	Cardio	RNASeq
Prelim.	$n$	102	63	2730
info.	$p$	5565	20426	8716
	# clusters	4	2	19
	SpKM	0.698	0.492	0.157
Competing	KM	0.721 (0.002)	0.493	0.157
methods	SAS	0.554 / 1307	0.492 / 1071	—
	Ex. Mass	no selection	no selection	0.157 / 6814
	IF-PCA	no selection	0.486 / 523	0.213 (0.002) / 4358
	SpKM	0.396 / 4	0.505 / 33	0.152 / 2304
COSCI +	KM	0.397 / 4	0.505 / 33	0.153 / 2304
	IF-PCA	0.381 / 4	0.505 / 33	0.156 / 2304

Table 6. Table 6: False Negatives and False Positive rates for Simulation Experiment V. Here, ℐ S = { 1 , 2 , 3 , 4 } subscript ℐ 𝑆 1 2 3 4 \mathcal{I}_{S}=\{1,2,3,4\} , p = 25 𝑝 25 p=25 . The numbers in parenthesis are standard errors over 10 10 10 repetitions.

		$n=2000$ / $p=25$
		Avg FN		Avg FP
	0.05	0.00	(0.00)	21.00	(0.00)
COSCI	0.08	0.00	(0.00)	21.00	(0.00)
with	0.1	0.00	(0.00)	20.9	(0.10)
${\alpha }_0$	0.12	0.00	(0.00)	20.3	(0.26)
fixed	0.15	0.00	(0.00)	18.8	(0.39)
	0.2	0.00	(0.00)	13.7	(0.68)
	0.25	0.00	(0.00)	3.8	(0.78)
	Data driven	0.00 (0.00)	3.1 (0.23)

Table 7. Table 7: False Negatives and False Positive rates for COSCI in studying the impact of correlation. Here, ℐ S = { 1 , 2 , 3 , 4 } subscript ℐ 𝑆 1 2 3 4 \mathcal{I}_{S}=\{1,2,3,4\} , p = 25 𝑝 25 p=25 . The numbers in parenthesis are standard errors over 10 10 10 repetitions.

		$n$ = 2000 / $p$ = 25
		Avg FN	Avg FP
	0.05	0.00 (0.00)	21.0 (0.00)
COSCI	0.08	0.00 (0.00)	20.7 (0.15)
with	0.1	0.00 (0.00)	19.9 (0.18)
${\alpha }_0$	0.12	0.00 (0.00)	18.4 (0.40)
fixed	0.15	0.00 (0.00)	16.0 (0.60)
	0.2	0.00 (0.00)	10.5 (0.65)
	0.25	0.00 (0.00)	4.0 (0.79)
	Data driven	0.00 (0.00)	1.00 (0.79)

Table 8. Table 8: Clustering Errors for Simulation Experiment I (left) and II (right).

		$n=200$	$n=1000$	$n=2500$
	0.05	0.178	0.242	0.178
COSCI	0.08	0.238	0.229	0.165
with	0.1	0.264	0.229	0.138
${\alpha }_0$	0.12	0.305	0.222	0.122
fixed	0.15	0.338	0.211	0.114
	0.2	0.344	0.183	0.121
	Data driven	0.195	0.172	0.133
	SpKM	0.484	0.46	0.47
Other	SpHC	0.495	0.491	0.494
methods	SAS	0.495	0.494	0.478
	Ex. Mass	0.234	0.203	0.203
	IF-PCA	0.478	0.484	0.488

Table 9. Table 9: Minimum CER of the competing methods \textcolor blackalong with the method name. For COSCI, the smallest CER observed for COSCI + “clustering method” is shown where “clustering method” includes the \textcolor blackfive competing methods in Table 5 . Number of selected features are reported after the ‘/’ symbol.

					Min. CER		COSCI
Data Set	Source	$n$	$p$	$K$	Method	CER	Min. CER	Method
ProstateCancer	Singh et al. (2002)	102	6033	2	IF-PCA	0.477 / 1551	0.498 / 178	$k$-means
Lymphoma	Alizadeh et al. (2000)	62	4026	3	IF-PCA	0.103 / 42	0.285 / 22	IF-PCA
Brain	Pomeroy et al. (2002)	42	5597	5	IF-PCA	0.159 / 453	0.131 / 323	IF-PCA
Colon	Alon et al. (1999)	62	2000	2	IF-PCA	0.490 / 25	0.444 / 3	$k$-means
SRBCT	Kahn (2001)	63	2308	4	SAS	0.327 / 246	0.350 / 84	$k$-means
Leukemia	Golub et al. (1999)	72	3571	2	IF-PCA	0.131 / 213	0.430 / 112	IF-PCA
SuCancer	Su et al. (2001)	174	7909	2	SAS	0.501 / 526	0.500 / 211	$k$-means
LungCancer(1)	Gordon et al. (2002)	181	12533	2	IF-PCA	0.064 / 251	0.206 / 420	$k$-means
LungCancer(2)	Bhattacharjee et al. (2001)	203	12600	2	IF-PCA	0.341 / 418	0.502 / 260	$k$-means
BreastCancer	Wang et al. (2005)	276	22215	2	IF-PCA	0.484 / 721	0.492 / 1195	$k$-means
Michigan	Subramanium et al. (2005)	86	5217	2	Ex. Mass	0.427 / 4	0.479 /51	$k$-means

Table 10. Table 10: This table reports the % percent \% of cases where COSCI detects clusters across varying sample sizes n 𝑛 n and thresholds α 0 subscript 𝛼 0 \alpha_{0} .

				${\alpha }_0$
Unimodal density	$n$	0.01	0.02	0.05	0.1	0.15	0.2	0.25
	100	100	100	99	74	47	28	12
	500	100	100	67	33	17	11	7
	1000	100	98	49	22	13	6	2
$𝒩(0,1)$	2000	100	82	21	10	3	0	0
	5000	94	38	3	1	1	0	0
	10000	47	6	0	0	0	0	0
	100	100	100	59	19	11	3	1
	500	99	50	4	0	0	0	0
	1000	74	14	0	0	0	0	0
Student $t_{\left(1\right)}$	2000	27	0	0	0	0	0	0
	5000	0	0	0	0	0	0	0
	10000	0	0	0	0	0	0	0
	100	100	100	96	59	33	14	7
	500	100	99	47	10	2	2	0
	1000	100	84	7	2	0	0	0
$ℰ\left(1\right)$	2000	95	33	0	0	0	0	0
	5000	35	0	0	0	0	0	0
	10000	2	0	0	0	0	0	0
	100	100	99	52	18	10	2	1
	500	99	62	7	1	0	0	0
	1000	82	17	0	0	0	0	0
Cauchy	2000	36	3	0	0	0	0	0
	5000	0	0	0	0	0	0	0
	10000	0	0	0	0	0	0	0
	100	100	100	75	29	13	2	1
	500	100	81	11	1	0	0	0
	1000	98	40	1	0	0	0	0
$ℒ\left(1\right)$	2000	71	8	0	0	0	0	0
	5000	4	0	0	0	0	0	0
	10000	0	0	0	0	0	0	0
	100	100	100	80	28	16	7	2
GEV with	500	100	82	15	3	0	0	0
shape	1000	99	48	3	0	0	0	0
parameter	2000	68	12	0	0	0	0	0
= 0.8	5000	13	0	0	0	0	0	0
	10000	0	0	0	0	0	0	0
	100	100	100	99	78	50	32	19
$ℬ(1,3)$	500	100	100	74	28	16	11	2
	1000	100	98	37	5	3	1	0
	2000	99	80	14	2	0	0	0
	5000	88	24	2	0	0	0	0
	10000	43	1	0	0	0	0	0
	100	100	100	99	84	54	33	22
Triangle distbn.	500	100	100	79	34	15	7	0
$\in [0,1]$	1000	100	98	42	11	4	4	1
with	2000	100	88	28	5	1	0	0
mode at 0.8	5000	98	40	3	0	0	0	0
	10000	68	11	0	0	0	0	0

Table 11. Table 11: Finding the population splits for 2 2 2 –normal mixtures: p 1 ​ 𝒩 ​ ( μ 1 , 1 ) + p 2 ​ 𝒩 ​ ( μ 2 , 1 ) subscript 𝑝 1 𝒩 subscript 𝜇 1 1 subscript 𝑝 2 𝒩 subscript 𝜇 2 1 p_{1}\,\mathcal{N}(\mu_{1},1)+p_{2}\,\mathcal{N}(\mu_{2},1) .

CASE	$p_1$	$p_2$	${\mu }_1$	${\mu }_2$	$\underset{¯}{m}$	$s^{\star }$	$s_{\text{MC}}$	Excess MCE
$\left|{𝝁}_{\mathrm{𝟐}}-{𝝁}_{\mathrm{𝟏}}\right|=\mathrm{𝟗}$	0.50	0.50	–4.50	4.50	0.00	0.00	0.00	$0.000$
	0.45	0.55	–4.50	4.50	–0.02	–0.45	–0.02	$0.000$
	0.40	0.60	–4.50	4.50	–0.05	–0.90	–0.04	$0.000$
	0.35	0.65	–4.50	4.50	–0.07	–1.36	–0.07	$0.000$
	0.30	0.70	–4.50	4.50	–0.10	–1.82	–0.09	$0.001$
	0.25	0.75	–4.50	4.50	–0.13	–2.31	–0.12	$0.004$
	0.20	0.80	–4.50	4.50	–0.16	–2.90	–0.15	$0.011$
	0.15	0.85	–4.50	4.50	–0.20	–3.82	–0.19	$0.037$
	0.10	0.90	–4.50	4.50	–0.26	NO	–0.24	$0.100$
$\left|{𝝁}_{\mathrm{𝟐}}-{𝝁}_{\mathrm{𝟏}}\right|=\mathrm{𝟖}$	0.50	0.50	–4.00	4.00	0.00	0.00	0.00	$0.000$
	0.45	0.55	–4.00	4.00	–0.03	–0.40	–0.03	$0.000$
	0.40	0.60	–4.00	4.00	–0.05	–0.80	–0.05	$0.000$
	0.35	0.65	–4.00	4.00	–0.08	–1.22	–0.08	$0.001$
	0.30	0.70	–4.00	4.00	–0.11	–1.64	–0.11	$0.003$
	0.25	0.75	–4.00	4.00	–0.15	–2.12	–0.14	$0.008$
	0.20	0.80	–4.00	4.00	–0.18	–2.72	–0.17	$0.020$
	0.15	0.85	–4.00	4.00	–0.23	NO	–0.22	$0.150$
	0.10	0.90	–4.00	4.00	–0.29	NO	–0.28	$0.100$
$\left|{𝝁}_{\mathrm{𝟐}}-{𝝁}_{\mathrm{𝟏}}\right|=\mathrm{𝟕}$	0.50	0.50	–3.50	3.50	0.00	0.00	0.00	$0.000$
	0.45	0.55	–3.50	3.50	–0.03	–0.35	–0.03	$0.000$
	0.40	0.60	–3.50	3.50	–0.06	–0.71	–0.06	$0.001$
	0.35	0.65	–3.50	3.50	–0.10	–1.09	–0.09	$0.003$
	0.30	0.70	–3.50	3.50	–0.13	–1.49	–0.12	$0.007$
	0.25	0.75	–3.50	3.50	–0.17	–1.97	–0.16	$0.016$
	0.20	0.80	–3.50	3.50	–0.22	–2.66	–0.20	$0.040$
	0.15	0.85	–3.50	3.50	–0.27	NO	–0.25	$0.150$
	0.10	0.90	–3.50	3.50	–0.34	NO	–0.31	$0.100$
$\left|{𝝁}_{\mathrm{𝟐}}-{𝝁}_{\mathrm{𝟏}}\right|=\mathrm{𝟔}$	0.50	0.50	–3.00	3.00	0.00	0.00	0.00	0.000
	0.45	0.55	–3.00	3.00	–0.04	–0.32	–0.03	0.001
	0.40	0.60	–3.00	3.00	–0.08	–0.64	–0.07	0.004
	0.35	0.65	–3.00	3.00	–0.12	–0.99	–0.10	0.008
	0.30	0.70	–3.00	3.00	–0.16	–1.39	–0.14	0.016
	0.25	0.75	–3.00	3.00	–0.21	–1.91	–0.18	0.034
	0.20	0.80	–3.00	3.00	–0.26	NO	–0.23	0.200
	0.15	0.85	–3.00	3.00	–0.33	NO	–0.29	0.150
	0.10	0.90	–3.00	3.00	–0.41	NO	–0.37	0.100
$\left|{𝝁}_{\mathrm{𝟐}}-{𝝁}_{\mathrm{𝟏}}\right|=\mathrm{𝟓}$	0.50	0.50	–2.50	2.50	0.00	0.00	0.00	0.000
	0.45	0.55	–2.50	2.50	–0.05	–0.30	–0.04	0.005
	0.40	0.60	–2.50	2.50	–0.10	–0.61	–0.08	0.011
	0.35	0.65	–2.50	2.50	–0.15	–0.96	–0.12	0.021
	0.30	0.70	–2.50	2.50	–0.20	–1.41	–0.17	0.041
	0.25	0.75	–2.50	2.50	–0.26	NO	–0.22	0.250
	0.20	0.80	–2.50	2.50	–0.33	NO	–0.28	0.200
	0.15	0.85	–2.50	2.50	–0.41	NO	–0.35	0.149
	0.10	0.90	–2.50	2.50	–0.53	NO	–0.44	0.100
$\left|{𝝁}_{\mathrm{𝟐}}-{𝝁}_{\mathrm{𝟏}}\right|=\mathrm{𝟒}$	0.50	0.50	–2.00	2.00	0.00	0.00	0.00	0.000
	0.45	0.55	–2.00	2.00	–0.07	–0.32	–0.05	0.015
	0.40	0.60	–2.00	2.00	–0.14	–0.67	–0.10	0.034
	0.35	0.65	–2.00	2.00	–0.21	–1.12	–0.15	0.065
	0.30	0.70	–2.00	2.00	–0.28	NO	–0.21	0.298
	0.25	0.75	–2.00	2.00	–0.37	NO	–0.28	0.248
	0.20	0.80	–2.00	2.00	–0.47	NO	–0.35	0.198
	0.15	0.85	–2.00	2.00	–0.58	NO	–0.43	0.148
	0.10	0.90	–2.00	2.00	–0.74	NO	–0.55	0.097
$\left|{𝝁}_{\mathrm{𝟐}}-{𝝁}_{\mathrm{𝟏}}\right|=\mathrm{𝟑}$	0.50	0.50	–1.50	1.50	0.00	0.00	0.00	0.000
	0.45	0.55	–1.50	1.50	–0.12	–0.50	–0.07	0.057
	0.40	0.60	–1.50	1.50	–0.24	NO	–0.14	0.396
	0.35	0.65	–1.50	1.50	–0.38	NO	–0.21	0.344
	0.30	0.70	–1.50	1.50	–0.53	NO	–0.28	0.293
	0.25	0.75	–1.50	1.50	–0.71	NO	–0.37	0.241
	0.20	0.80	–1.50	1.50	–1.50	NO	–0.46	0.190
	0.15	0.85	–1.50	1.50	–1.50	NO	–0.58	0.139
	0.10	0.90	–1.50	1.50	–1.50	NO	–0.73	0.090

Equations

$$\underset{c_1,\dots ,c_n\in \mathbb{R}}{\min }\sum _{i=1}^n(x_i-c_i{)}^2+\lambda \sum _{1\le k<\ell \le n}|c_k-c_{\ell }|.$$

$$\alpha =n^{-1}\min (|{\mathcal{C}}_1|,|{\mathcal{C}}_2|)1\{n^{-1}(|{\mathcal{C}}_1|+|{\mathcal{C}}_2|)\ge 0.5\},$$

$$\underset{j\in {\mathcal{I}}_N}{\max }{S}_j(\tau )\le b\frac{\ln (p_N\vee n)}{n}$$

$$\mathrm{Pr}(\{j:S_j(\tau )>{\alpha }_0^{\mathrm{OR}}\}={\mathcal{I}}_S)\ge 1-c_1{p}^{-c_2},$$

$${\hat{\alpha }}_0=\underset{j\in {\hat{\mathcal{I}}}_S}{\min }{S}_j,$$

$$T_j={\hat{\pi }}_0{\hat{f}}_0({\psi }_j)/\hat{f}({\psi }_j)$$

$${\mu }_1=(0,0),{\mu }_2=(0,-4),{\mu }_3=(4,0),{\mu }_4={\mu }_3+{\mu }_2,$$

$${\Sigma }_1={\Sigma }_4=\left(\begin{array}{cc}1& -0.85\\ -0.85& 1\end{array}\right),{\Sigma }_2={\Sigma }_3=\left(\begin{array}{cc}1& 0.85\\ 0.85& 1\end{array}\right).$$

$${\hat{\mathcal{I}}}_S=\text{unique}(\{i:S_i\ge {\alpha }_0\}\cup \{(i,j):S_{(i,j)}\ge {\alpha }_0\}).$$

$${\hat{\mathcal{I}}}_S^{\mathrm{thr}}$$

$${\hat{\mathcal{I}}}_S^{\mathsf{mod}}$$

$${\mu }_1=(0.9,-0.9),{\mu }_2=-{\mu }_1,\Sigma =\left(\begin{array}{cc}1& 0.9\\ 0.9& 1\end{array}\right)$$

$$\underset{j\in {\mathcal{I}}_N}{\max }\Vert g_j{\Vert }_{\infty }+\Vert g_j^{\prime }{\Vert }_{\infty }$$

$$\underset{j\in {\mathcal{I}}_N}{\max }|q_j(\tau )|+|q_j(1-\tau )|$$

$$\underset{j\in {\mathcal{I}}_N}{\min }{g}_j\{q_j(\tau )\}\wedge g_j\{q_j(1-\tau )\}$$

$$\underset{j\in {\mathcal{I}}_N}{\max }\mathrm{E}|X_j|\le 1.$$

$$|\{{\hat{G}}_{L_j,R_j}^j(R_j)-{\hat{G}}_{L_j,R_j}^j(a_j)\}-\{G_{L_j,R_j}^j(R_j)-G_{L_j,R_j}^j(a_j)\}|\le {ϵ}_1(R_j-a_j)+c_3\frac{\ln (p\vee n)}{n},$$

$$\begin{gathered} |\{{\hat{G}}_{L_j,R_j}^j(L_j)-{\hat{G}}_{L_j,R_j}^j(a_j)\}-\{G_{L_j,R_j}^j(L_j)-G_{L_j,R_j}^j(a_j)\}| \\ \le {ϵ}_1(a_j-L_j)+c_3\frac{\ln (p\vee n)}{n}, \end{gathered}$$

$$|P_{nj}(R_j-a_j)-P_j(R_j-a_j)|\le {ϵ}_1(R_j-a_j)+c_3\frac{\ln (p\vee n)}{n},$$

$$|P_{nj}(a_j-L_j)-P_j(a_j-L_j)|\le {ϵ}_1(a_j-L_j)+c_3\frac{\ln (p\vee n)}{n},$$

$$\{G_{L_j,R_j}^j(L_j)-G_{L_j,R_j}^j(a_j)\}\vee \{G_{L_j,R_j}^j(R_j)-G_{L_j,R_j}^j(a_j)\}\ge c_0(R_j-a_j)\wedge (a_j-L_j).$$

$$0$$

$$(R_j-a_j)\wedge (a_j-L_j)\le 2c_0^{-1}{c}_3\frac{\ln (p\vee n)}{n}.$$

$$\Pr\Bigl{\{}\max_{j}\sup_{\ell<r}|P_{nj}(\ell,r)-P_{j}(\ell,r)|>\epsilon\Bigr{\}}\leq a_{1}\exp(-a_{2}n-a_{2}\ln n-a_{3}\ln p),$$

$$|{\hat{\mu }}_{a,R}^j-{\mu }_{a,R}^j|$$

$$|({\mu }_{\ell ,r}^j-{\mu }_{L,a}^j)-({\hat{\mu }}_{\ell ,r}^j-{\hat{\mu }}_{L,a}^j)|$$

$$\mathrm{Pr}(D_n>ϵ)\le a_1\exp (-a_{2}n-a_2\ln n-a_3\ln p),$$

$$|P_{nj}(\ell ,r)-P_j(\ell ,r)|\le {ϵ}_1(r-\ell )+\frac{\ln (p\vee n)}{n}{M}_n$$

$$\mathrm{Pr}(M_n>c_3)$$

$${\mu }_1=(0.9,-0.9),{\mu }_2=-{\mu }_1,\Sigma =\left(\begin{array}{cc}1& 0.9\\ 0.9& 1\end{array}\right).$$

$$T_j=\widehat{{\pi }_0}\widehat{f_0}({\psi }_j)/\hat{f}({\psi }_j)$$

Full text

1 Introduction

We consider the problem of feature screening in large scale cluster analysis. Clustering is one of the most popular unsupervised classification techniques; \textcolorblackit is widely used in a myriad of statistical applications for stratification and sub-population identification [21, 29, 45]. In recent years, due to massive advancements in the modern data collection and assimilation techniques, very big datasets, with both a large number of observations and a large number of features, have been generated with increasing frequency. Classical clustering methods (see, \textcolorblacke.g., Chap. 14 of Friedman et al. [22]) are either computationally challenging or ineffective for conducting segmentation analysis of such massive modern data. In many scientific applications, when the dimension of the data is very high, most of the coordinates (i.e., features) contain very little information regarding the grouping structure. Classical clustering methods, which do not reduce the dimension of the data, suffer, because the agglomerative effects of the large number of “noise” features conceal important clustering information available in a relatively smaller number of “signal” features.

Recently developed clustering algorithms, which exploit the underlying sparseness, are effective in dealing with high-dimensional data; see, e.g., [3, 10, 35, 56]. \textcolorblackWe propose herein a scalable computationally efficient approach, entitled COSCI (COnvex Screening for Cluster Information), that can efficiently weed out the features that are non-informative for clustering. As a \textcolorblacknonparametric approach, COSCI has competitive advantages over the popular Gaussian mixture based parametric techniques [3, 4]. Unlike the \textcolorblacknonparametric density estimation based screening techniques, COSCI is very scalable, and can successfully handle datasets with more than one million observations. Our proposed procedure discards non-informative features by first computing a clustering score for the clustering tree constructed for each feature, and then thresholding the resulting values. We provide the theoretical motivation for our approach by establishing uniform non-asymptotic bounds on the clustering scores of the noise features.

Our theoretical results are a significant extension of the univariate results in Radchenko and Mukherjee [43]. \textcolorblackWe rely herein on a more careful analysis of the empirical processes corresponding to the clustering trees constructed for each feature by the associated clustering procedure. We derive a stronger tail probability bound for the clustering score of each feature, which then allows us to establish a uniform bound for all the clustering scores of the non-informative features in the high-dimensional setting, where the number of such features can be extremely large. Using this uniform bound, we infer that under mild regularity conditions on the population distribution, the proposed COSCI algorithm will discard the non-informative features and perfectly select the informative features with very high probability.

\textcolor

blackHigh-dimensional datasets arising in modern biology, econometrics, engineering, text mining and signal processing [33] require a significant degree of feature screening for subsequent application of a clustering algorithm. We illustrate the applicability of COSCI for making scientific discoveries through the analysis of single cell biological data. Emerging technologies, such as single-cell mass cytometry [6], next-generation sequencing [39] and micro-fluidic methods [15, 55], have recently enabled us to collect gene and protein expression information for each cell [53]. The resulting datasets, which are not only very high-dimensional but also contain a large number of cellular observations, serve as invaluable resources for the characterization of the cellular hierarchy in multi-cellular organisms. Often, the biological question associated with these datasets is the identification of homogeneous cellular sub-populations based on the differential expression patterns of the genes. The composition of these sub-populations is subsequently analyzed to detect interesting structures.

Recently, several algorithms such as viSNE [1], Wanderlust [5], ECLAIR [26], SPADE [42], SLIDE [46] and Scaffold maps [49] have been developed for conducting such sub-population analysis for single cell data. Most of these methods conduct either dimension reduction through PCA related methods or handle large sample sizes via down-sampling. While these algorithms are widely used, they lack appropriate mathematical guarantees to show that the resulting sub-populations are not due to random fluctuations and will be reproducible across datasets generated from experiments conducted under similar conditions. Furthermore, the use of techniques such as PCA to reduce dimensionality in these settings can be called into question [11] because (a) the inferred sub-populations may not be sparse in the expression patterns of the genes, in which case PCA has been theoretically proven to produce inconsistent results [36]; (b) there is no guarantee that cluster information is aligned in the direction of maximum variance.

As a motivational example, consider the problem of cellular sub-population detection in a single-cell RNAseq data analyzed in Giecold et al. [26] (henceforth referred to as G16). This massive dataset holds the expression levels of $p=8716$ genes for $n=2730$ mice bone marrow cells [41]. \textcolorblackOne of the key scientific objectives is to infer the lineage pattern of the identified sub-populations based on $33$ lineage markers. If the sub-populations differ with respect to a relatively small subset of all the genes considered, then, from the statistical perspective, the problem reduces to screening out the non-informative features from the data consisting of vectors \textcolorblack$X_{1},\ldots,X_{p}\in\mathbb{R}^{n}$ and identifying a subset of features that retains the cluster information. As such, identifying the best subset of the genes with respect to the cluster information is an important statistical endeavor with critical biological implications; \textcolorblacksee [56] and the references therein.

Unlike the methods that use downsampling, COSCI can accommodate sample sizes on the order of $10^{6}$. COSCI first produces a score, $S_{j}\in(0,0.5]$, for each feature, which reflects its relative importance for clustering, and then screens out the features with lower scores. When applied to the aforementioned RNASeq data, COSCI orders the scores of the $8716$ genes and selects the top $2304$ genes, which include $12$ of the $33$ lineage markers (see Figure 1, bottom left; the selected genes are highlighted in solid black and the black dots are the $33$ lineage markers). The heatmap (Figure 1, bottom right) of the expressions for the $19$ sub-populations detected via $k$-means on the COSCI selected features shows that the inferred sub-populations differ significantly across the $12$ selected lineage markers. After applying $k$-means to detect the sub-populations on the $2304$ genes selected by COSCI we got a misclassification error (computed as CER = $1-$ Rand Index, [14, 44]) of approximately $0.15$. This error is significantly smaller than several other methods that are widely used for such clustering problems especially when the proportion of selected features is taken into account (see Figure 1, top right). We revisit this example with more details in Section 4.2.

1.1 Connections to related statistical literature

Within the statistical literature, a number of recently proposed clustering approaches execute feature screening as the first step and then rely on conventional clustering techniques, such as $k$-means, to cluster the remaining data. For example, Chan and Hall [10] proposed a \textcolorblacknonparametric feature screening method that is based on coordinate-wise \textcolorblackexcess mass tests [12]. They rank the features using the values of the corresponding test statistic. Feature selection then follows by identifying a kink in the plot of the \textcolorblackwithin-cluster sum of squares versus the number of identified clusters. Witten and Tibshirani [56] proposed the sparse $k$-means and sparse hierarchical clustering approaches, which employ $k$-means and hierarchical clustering, respectively, on a feature weighted dissimilarity matrix, where the weights are encouraged to be sparse. Their method is largely inspired by the popular COSA algorithm of [23] and is more adept at sparse clustering. Recently, Arias-Castro and Pu [2] proposed Sparse Alternate Sum (SAS) clustering, which uses a hill-climbing approach to solve the sparse $k$-means optimization problem.

On the parametric side, several \textcolorblackmodel-based clustering approaches have been introduced [40, 54, 57]. These techniques typically maximize a penalized likelihood under a Gaussian mixture model, where the penalization serves the purpose of implicit feature selection. Jin et al. [34] and Jin et al. [35] propose IF-PCA, which is a \textcolorblacktwo-step clustering method — the first step conducts coordinate wise feature selection, and the second step performs $k$-means clustering on the matrix of left singular vectors of the selected features. The feature selection step uses the \textcolorblackKolmogorov–Smirnov test for normality to rank the features, followed by the use of the Higher Criticism (HC) [16, 17] functional to finally select the features. Theoretical properties of clustering algorithms that combine feature selection with clustering have also been recently studied. For example, Azizyan et al. [4] provide information theoretic bounds on clustering accuracy of the \textcolorblackhigh-dimensional Gaussian mixtures, while Arias-Castro and Verzelen [3] establish minimax rates for the problems of mixture detection and feature selection under the sparsity assumption.

Our work is closer to the approaches of [10], [35] and [56], where the objective is to screen out the noise features. We analyze the problem of feature screening in large scale clustering and propose COSCI — a novel computationally efficient screening procedure with strong theoretical motivation. COSCI uses a \textcolorblacknonparametric approach to rank order the features by their clustering leverage. In this respect, it differs from the recently proposed screening techniques, such as IF-PCA, which rely on a parametric family as a point of reference to gauge feature strength for clustering.

1.2 Organization of the paper

In Section 2 we present and discuss our screening methodology. More specifically, Algorithm 1 provides the details of the implementation, while Section 2.1 contains the main theoretical motivation and associated results. Section 3 provides two approaches that aid the selection of the screening threshold. In Section 4 we conduct a detailed empirical analysis of our approach using both simulated data and real data from microarray experiments. Section 5 concludes the paper. Proofs and additional technical details are relegated to the Appendix.

2 Methodology and main results

We consider the problem of clustering $n$ observations based on $p$ features in the setting where most of the features contain no clustering information. Noting that these “noise” features have unimodal marginal distributions (which may differ across the features), we develop a univariate approach, which, based on the sample observations, evaluates whether the true underlying density is unimodal. The theoretical and empirical results provided in Sections 2.1 and 4 demonstrate that our approach is successful at screening out the noise and identifying the signal features even in very high-dimensional scenarios.

The key ingredients of our proposed COSCI methodology are a univariate merging algorithm, which constructs a clustering tree for each of the features, and a screening of the merge sizes, which identifies the signal features as those with at least one sufficiently large merge. We discuss each of these components in detail below.

Univariate merging algorithm. With the goal of checking unimodality for each of the feature coordinates, we consider the following univariate optimization problem:

[TABLE]

It is based on the observations $x_{i}$ and corresponds to minimizing the \textcolorblackwithin-cluster sum of squares under constraints on the \textcolorblack$L_{1}$-distance between the cluster centroids $c_{k}$. Here $\lambda$ is a non-negative penalty weight. The convexity of the objective function in (1) has been exploited to develop algorithms for efficiently producing the path of solutions as a function of the penalty weight [31, 32, 43]. Clustering algorithms based on fusion penalization of this type have become very popular in large scale clustering [13, 31, 43, 51, 58] and regression analysis [8, 37, 47, 48].

The entire path of solutions corresponding to the objective criterion (1) can be found by a simple merge algorithm in $O(n\ln n)$ operations. Starting with $n$ observations in $n$ clusters, we sequentially merge the nearest (in terms of the weighted distance as shown in Algorithm 1) adjacent centroids until we are left with just one cluster in the end.

Merge sizes. Given a merge of sub-clusters $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$, we calculate its size, which we define as

[TABLE]

where $|\cdot|$ denotes the cardinality of a set. The above thresholding of the mass after merge ensures that we only identify a merge as big if it results in a significantly large cluster. This protects us from the risk of discovering potential big merges on smaller fragments of the sample, where the nature of these merges can be very fragile due to sampling fluctuations.

The fundamental working principle of our proposed COSCI approach rests on the following property of the merge sizes: if $x_{1},\ldots,x_{n}$ are indeed generated by a non-informative density, then the sample merges $\alpha_{1},\ldots,\alpha_{n-1}$ will be uniformly small for sufficiently large $n$. A multivariate version of this property is formalized in Theorem 1 below. \textcolorblackIn contrast, if the underlying distribution contains a moderate amount of cluster information, then the merge sequence will have at least one merge that is big. The last fact is illustrated by Theorem 1 in [43]; see also the discussion above Corollary 1. Based on these properties of the merge sizes, we conduct the following screening procedure.

Screening the merges. Given a pre-defined threshold $\alpha_{0}$, we flag the feature as potential “signal”, if there exists a corresponding merge, say the \textcolorblack$i$th merge, such that $\alpha_{i}\geq\alpha_{0}$.

Our proposed methodology is formalized in Algorithm 1 below.

Algorithm 1: COSCI procedure for feature screening

INPUT: Data matrix ${X}_{n\times p}$ and tuning parameter $\alpha_{0}$.

FOR each $j\in\{1,\ldots,p\}$:

INITIALIZE:

Sort data in ascending order and store them as $x=\{x_{1},\ldots,x_{n}\}$.

Set $k$, the number of clusters, equal to n. For each $r\in\{1,\ldots,n\}$, set $c_{r}=\{x_{i}\}$.

REPEAT:

 Find consecutive adjacent centroid distances: $d(r,r+1)\leftarrow(c_{r+1}-c_{r})/(|c_{r}|+|c_{r+1}|)$.

 Find clusters with minimum merging distance: $r^{*}\leftarrow\arg\min_{r}d(r,r+1)$.

 Merge clusters $r^{*}$, $r^{*}+1$, re-label remaining clusters and set $k\leftarrow k-1$.

 Find the merge size, $\alpha_{n-k}^{j}$, using equation (2).

UNTIL $k=1$

STORE: Clustering score: $S_{j}=\max_{1\leq k\leq n-1}\alpha_{k}^{j}$

FEATURE SCREENING: $\hat{\mathcal{I}}_{S}=\{j:S_{j}\geq\alpha_{0}\}$

See A.3.1 for a detailed description of the computational steps involved in \textcolorblackAlgorithm 1. We note that a significant gain in computational time is achievable via a parallel implementation of the top for loop in Algorithm 1 that runs across the $p$ features.

2.1 Theoretical support: Perfect screening property

Let $\mathcal{I}_{S}$ and $\mathcal{I}_{N}$ be the index sets corresponding to the “signal” and the “noise” features, respectively. Define $p_{S}=|\mathcal{I}_{S}|$ and $p_{N}=|\mathcal{I}_{N}|$. Given feature $j$, we write $S_{j}(\tau)$ for the corresponding largest merge size, computed the same way as $S_{j}$ in Algorithm 1, but under the restriction that the midpoint between the two merged sub-clusters lies between the \textcolorblacksample quantiles of order $100\times\tau$ and $100\times(1-\tau)$, where $\tau$ is an arbitrarily small but positive number.

Theorem 1, stated below, establishes a uniform non-asymptotic bound on the merge sizes that are produced when our procedure is applied to the noise features. The regularity conditions, C1 and C2, are imposed on the family of marginal distributions of the noise features, and are fairly mild. In particular, they are satisfied for location-scale families of unimodal differentiable densities with a finite first moment. The proof of Theorem 1 is provided in A.1. \textcolorblackIn what follows, $a\vee b=\max(a,b)$ and $a\wedge b=\min(a,b)$, for any $a,b\in\mathbb{R}$.

Theorem 1**.**

Suppose that regularity conditions C1 and C2, stated in A.1, are satisfied. For each $\tau>0$ there exist positive constants $c_{1}$, $c_{2}$, $b$ and $\kappa$, whose choice does not depend on either $n$ or $p_{N}$, such that, as long as $p_{N}\leq\exp(\kappa n)$, inequalities

[TABLE]

hold with (high) probability bounded below by $1-c_{1}p_{N}^{-c_{2}}$.

The above theorem provides the theoretical justification for the screening step in our proposed procedure. The result is non-asymptotic, and the proof involves careful analysis of the empirical process associated with the merging algorithm for each feature. Under a mild restriction on the number of features relative to the sample size, the above theorem ensures that the clustering scores of all the noise features are uniformly very close to zero. Thus, if we use any arbitrarily small but prefixed value for the threshold $\alpha_{0}$, we have theoretical guarantees for perfectly screening out all the noise coordinates.

It is important to have a small value of $\alpha_{0}$ to avoid screening out the informative features, which have non-negligible clustering scores $S_{j}$. Theorem 1 suggests that $\alpha_{0}^{\rm OR}=b\ln(p\vee n)/n$ is a reasonable choice. Provided $n$ is sufficiently large, an approach using the above choice of $\alpha_{0}$ will not screen out the features identified as multi-cluster features by the population clustering procedure, defined in Section 2.2 of [43]. The next result, which is a consequence of Theorem 1 above and Theorem 1 in [43], formalizes this point. Note that [43] demonstrate, through simulations and theoretical analysis, that the population procedure generally classifies multi-modal distributions as multi-cluster, provided the corresponding sub-populations are of reasonable size and have a moderate amount of separation. To illustrate this fact, in A.3.4 we provide a detailed summary of how the population procedure performs on a wide variety of bimodal Gaussian distributions.

Corollary 1**.**

Suppose that regularity conditions C1 and C2, stated in A.1, are satisfied. Let the cardinality of the set $\mathcal{I}_{S}$ be bounded above by a universal constant. Suppose that the population clustering procedure identifies each feature in $\mathcal{I}_{S}$ as multi-cluster. Then, for all sufficiently small $\tau>0$ there exist positive constants $c_{1}$, $c_{2}$, $b$ and $\kappa$, whose choice does not depend on either $n$ or $p$, such that

[TABLE]

as long as $p\leq\exp(\kappa n)$.

We view $\alpha_{0}^{\rm OR}$ as the oracle choice of the threshold. However, it is difficult to evaluate it from the data, primarily because the constant $b$ depends on the marginal densities $g_{j}$ of the noise features. In the following section, we discuss several practical choices for the threshold parameter.

3 Estimation of hyperparameters

The COSCI procedure presented in Algorithm 1 requires only one tuning parameter, $\alpha_{0}$, as an input. The oracle threshold choice $\alpha_{0}^{\rm OR}$ is difficult to estimate from the data, as the marginal distributions of the features are typically unknown. In what follows, we present two approaches for estimating $\alpha_{0}$ that are adaptive to the sample size $n$, in the sense that a larger threshold is chosen for smaller sample sizes.

3.1 Simulation based

We generate data of varying sample sizes from several well-known unimodal distributions and use it to assess at which values of $\alpha_{0}$ COSCI will detect no clusters, screen out the corresponding non-informative feature. The eight unimodal distributions considered in Table 10 are meant to represent the spectrum of the noise coordinates that are commonly encountered in real data applications. They include symmetric densities with the support equal to the entire \textcolorblackreal line $\mathbb{R}$, densities with heavy tails and those with bounded support.

Table 10 in A.3.2 presents the results of this simulation exercise over $100$ repetitions. For example, when the noise coordinate is Gaussian, and the sample size is $n=500$, COSCI detects clusters in the majority of the $100$ repetitions when $\alpha_{0}\leq 0.05$. Thus, an appropriate threshold for this case should at least be greater than $0.05$. For larger sample sizes, COSCI detects no clusters with a relatively smaller $\alpha_{0}$. A general theme that emerges from this table is that when the underlying density is non-Gaussian with support over \textcolorblackall of $\mathbb{R}$ or $\mathbb{R}_{+}=(0,\infty)$, smaller thresholds seem to succeed at screening out the corresponding feature, when compared to the Gaussian case for the same sample size. Similarly to the Gaussian case, densities with bounded support, such as the Beta and the Triangular distribution, require a larger threshold to succeed. Our practical recommendation is to assume an underlying Gaussian noise distribution and use $\alpha_{0}$ as the smallest threshold that detects no clusters given the sample size $n$. Let $\hat{\alpha}_{0}$ be such a threshold. Then, the selected feature set is $\mathcal{\hat{I}}_{S}=\{j:S_{j}\geq\hat{\alpha}_{0}\}$.

3.2 Data driven

In this section, we use a data driven technique to estimate $\alpha_{0}$. We work under the large-scale multiple testing framework of Efron [19] and transform the problem of estimating $\alpha_{0}$ into a problem of feature selection using the merge sizes \textcolorblack$S_{1},\ldots,S_{p}$. One can then read off the optimal $\alpha_{0}$ from the selected features as

[TABLE]

where $\mathcal{\hat{I}}_{S}$ holds the indices of the selected features. Let $\psi_{j}=2S_{j}$ be the test statistic for testing the significance of cluster strength in feature $j$. Note that \textcolorblackthe $\psi_{j}$’s have a mixture density $f$ given by $f(\psi)=\pi_{0}f_{0}(\psi)+(1-\pi_{0})f_{1}(\psi)$, where $f_{0}$ is the theoretical null distribution and $\pi_{0}$ is the null prior probability. The exact distributional form of $f_{0}$ is, however, unknown, primarily because we do not know the underlying distribution of the noise coordinates that generate the $S_{j}$’s. Nonetheless, Table 10 ascertains that $f_{0}$ is \textcolorblackright-skewed on the support $[0,1]$, with the mass concentrated around zero for large $n$. We use the MLE method of [19] to estimate the empirical null distribution from the observed $\psi_{j}$’s as a \textcolorblackBeta distribution and obtain the estimated false discovery rate (fdr)

[TABLE]

under the prescribed assumption that $\pi_{0}\geq 0.9$, i.e., \textcolorblackat least $90\%$ of the $p$ tests are null and $f_{1}(\psi)=0$ on $\mathcal{A}$, where $\mathcal{A}=\{\psi_{(1)}\leq\cdots\leq\psi_{([0.9p])}\}$. This scheme works reasonably well in both the simulations and the real data examples that we considered. For additional information about fdr see, e.g., Chap. 5 of Efron [18]. For estimating the mixture density $f$, we adopt Lindsey’s method [20, 38] that models the histogram bin counts using Poisson regression, treating the bin centers as covariates. Finally, to select the features we adopt the two-stage approach to signal screening recently introduced in Cai and Sun [9], the details of which are relegated to A.3.3.

In the empirical analysis that follows, we estimate $\mathcal{\hat{I}}_{S}$ and $\hat{\alpha}_{0}$ using the method described above.

4 Empirical analysis

4.1 Simulations

We perform several simulation experiments to gauge the feature screening performance of COSCI under two scenarios: (i) $p<n$ and (ii) $p>n$. For scenario (i), we implement two experiments and use the following four competing approaches of feature screening in clustering to compare the performance of our proposed method:

Sparse $k$-means clustering (SpKM) and hierarchical clustering (SpHC) [56] — we use the R-package sparcl. 2. 2.

Sparse Alternate Sum (SAS) clustering [2] — we use the R codes available at the author’s site (see: https://github.com/victorpu/SAS_Hill_Climb). 3. 3.

Influential Features PCA (IF-PCA) [35] — we use the MATLAB codes available at the author’s site (see http://www.stat.cmu.edu/~jiashun/Research/software/HCClustering/). 4. 4.

Excess Mass method (Ex. Mass) [10] for which unfortunately a software implementation is not available in the public domain. Noting that the excess mass test and the dip test are equivalent in \textcolorblacka univariate setting [12], we implement a version of this method for coordinate-wise feature screening using Hartigan’s Dip test [27, 28]. For feature selection, we select those features for which the multiplicity adjusted [7] \textcolorblack$p$-values from this test are at most 0.05.

For each of the above methods, we are only interested in their feature selection capabilities and not on their clustering performance.

For simulation Experiment I, we consider features from a wide range of parametric distributions including correlated features and consider three different sample sizes from low to high. Simulation Experiment I represents scenario (i); we fix $p=50$ and consider a design matrix ${X}_{n\times 50}$ with $p_{S}=5$ and $p_{N}=45$. The $p_{N}$ noise coordinates are taken to be iid standard Gaussian, \textcolorblack$\mathcal{N}(0,1)$, while the $p_{S}$ signal coordinates are chosen as follows:

$X_{1}\sim 0.5\,\mathcal{B}(4,6)+0.5\,\mathcal{B}(7,3)$ 2. 2.

$X_{2}\sim 0.5\,\mathcal{LN}(0.2,0.35)+0.5\,\mathcal{N}(4,0.5)$ 3. 3.

$X_{3}\sim 0.5\,\mathcal{L}(3,1.5)+0.5\mathcal{L}(5,1.5)$, \textcolorblackwhere $\mathcal{L}$ refers to the Laplace or double exponential distribution, and, 4. 4.

$(X_{4},X_{5})\sim\sum_{i=1}^{4}\mathcal{N}({\mu}_{i},\Sigma_{i})/4$, where

[TABLE]

In this setting, the signal coordinates are all bi-modal and with the exception of $X_{1},X_{3}$, the separation between the adjacent medians is fairly large. For Experiment II, we fix $p=100$ and consider a design matrix ${X}_{n\times 100}$ with $p_{S}=6$ and $p_{N}=94$. We let half of the $p_{N}$ noise coordinates to be iid standard Gaussian and the other half to be iid \textcolorblackStudent $t$ random variables with 5 degrees of freedom. The $p_{S}$ signal coordinates are chosen as follows:

$(X_{1},\ldots,X_{5})\sim$ as in Experiment I 2. 2.

$X_{6}\sim\,0.3\,\mathcal{N}(-2.5,1)+0.3\,\mathcal{N}(0,1)+0.4\,\mathcal{N}(2.5,1)$ which is a non-symmetric, tri-modal density.

For each of the setups described above, columns of the data matrix ${X}$ are simulated for \textcolorblackeach $n\in\{200,1000,2500\}$. We analyze two variants of COSCI: (i) COSCI with $\alpha_{0}$ fixed over a grid of six values, \textcolorblackviz. $0.05,0.08,0.1,0.12,0.15,0.2$, and, (ii) COSCI with $\alpha_{0}$ estimated in a data driven fashion as discussed in Section 3.2. For each method we calculated two statistics averaged over $50$ repetitions: False Negative (FN) — the number of signal features incorrectly identified as noise; False Positive (FP) — the number of noise features incorrectly identified as signal.

From Table 1, it is evident that SpKM fails to detect \textcolorblackat least three out of the five signal features across the three sample sizes, while SpHC and SAS have relatively better FN performance. The SpHC algorithm faced \textcolorblackscalability issues (marked by dash in the table) for \textcolorblacksample size $n=2500$. In contrast, IF-PCA correctly identifies the five signal features. All of these four methods have high false positives. Variations in sample sizes do not seem to affect their performance in any significant manner. Ex. Mass consistently fails to identify multi-modality in $X_{1}$ and $X_{3}$ but has the best false positive performance across all the methods. This is not unexpected given that the noise features considered are uni-modal. For $\alpha_{0}$ small, COSCI identifies the five signal features (low FN rate) but also incorrectly includes many noise features as signals (high FP rate). As expected, this trend transitions into a high FN rate and low FP rate for larger $\alpha_{0}$’s. This is where the data driven approach to choose $\alpha_{0}$ is seen to be beneficial. For moderately large $n$, COSCI, coupled with the data driven approach, returns a FN rate comparable to IF-PCA but with a significantly smaller FP rate. The results in Table 2 reveal a similar picture. When both FN and FP rates are taken into consideration, COSCI delivers the best performance amongst all the competing methods for moderately large $n$. Even when $n=200$, COSCI returns a FN rate which is only slightly higher than IF-PCA and SAS but enjoys a far better FP rate.

The next two simulation experiments III and IV exemplify scenario (ii) where $p>n$. We consider noise features which include non-symmetric distributions like \textcolorblackthe Exponential with rate parameter 1, \textcolorblack$\mathcal{E}(1)$, and heavy-tailed distributions like \textcolorblackthe standard Cauchy. This presents an interesting setting especially for methods like IF-PCA that rely on statistical comparison with a fixed parametric distribution to determine feature importance for clustering. In experiment III, we have $p=5000$ and $p_{S}=7$ with

$(X_{1},\ldots,X_{6})\sim$ as in Experiment I 2. 2.

$X_{7}\sim 0.5\,\mathcal{N}(-1.1,1)+0.5\,\mathcal{N}(1.1,1)$.

Note that $X_{7}$ is a bi-modal Gaussian mixture but with a relatively small separation between the modes. The $p_{N}$ noise features consists of approximately $40\%$ iid \textcolorblack$\mathcal{E}(1)$ noise \textcolorblackalong with $30\%$ iid standard Gaussian and $30\%$ iid $t_{5}$ noise. For Experiment IV, we consider $p=\mbox{25,000}$ and add two more signal features so that $p_{S}=9$ with

$(X_{1},\ldots,X_{7})\sim$ as in Experiment III 2. 2.

$X_{8}\sim 0.3\mathcal{L}(-3,1)+0.35\mathcal{L}(0,1)+0.35\mathcal{L}(3,1)$. 3. 3.

$X_{9}\sim 0.3\,\mathcal{B}(8,2)+0.35\,\mathcal{B}(5,5)+0.35\,\mathcal{B}(2,8)$.

Here $X_{8},X_{9}$, once again, represent features with a relatively small separation between the modes and, thus, are particularly difficult examples for modality detection. The $p_{N}$ noise features here include approximately $28\%$ iid standard Cauchy noise \textcolorblackalong with $24\%$ iid standard Gaussian, $24\%$ iid t5 and $24\%$ iid $\mathcal{E}(1)$ noises. We keep all the other design parameters of Experiment IV identical to Experiment III. However, we do not include the performance of SpKM, SpHC and SAS in our analysis of Experiment IV as these algorithms are computationally very demanding and often exhibited convergence issues in this regime.

From Tables 3 and 4, a drop in the FN performance of both IF-PCA and Ex. Mass is conspicuous. For IF-PCA, the non-Gaussian noises are identified as signals which are ultimately selected in favor of the true signals under the HC functional. \textcolorblackIn contrast, Ex. Mass continues to conclude \textcolorblackthat some of the difficult multi-modal signals like $X_{1},X_{3},X_{7},X_{8}$ and $X_{9}$ \textcolorblackare unimodal. COSCI, with the data driven approach, once again returns the best performance across all the different sample size regimes considered in these two experiments.

4.2 Real data examples

We test the performance of COSCI on a number of real data examples. In all these datasets, the number of clusters / subpopulations is known \textcolorblacka priori. However, unlike the simulation study, we do not know the “true” feature set for these datasets and thus an estimate of screening performance based on False Negatives or False Positives is impossible. Instead, we use the following scheme to gauge the performance of COSCI. On each of the datasets considered in this section, we overlay COSCI with $k$-means (KM), Sparse $k$-means (SpKM) and IF-PCA. In other words, we allow COSCI to screen and select the best features and thereafter, we run classical $k$-means, Sparse $k$-means and IF-PCA on the selected features. While using IF-PCA on the COSCI screened features, we only use the clustering component of IF-PCA and not \textcolorblackits feature selection step. The classification error rates (CER) so obtained from the above \textcolorblackthree schemes are then compared to the CER of the competing methods, \textcolorblacki.e., $k$-means, Sparse $k$-means and IF-PCA, all without any COSCI screening. For $k$-means and IF-PCA, we report the average CER over 30 independent replications and the associated standard error whenever they are bigger than $0.0005$. Next, we describe in details the application of the aforementioned methods on three real \textcolorblackdatasets. In A.2.4 the classification results on eleven other datasets are also demonstrated.

Multi-tissue data: This is a microarray data on different mammalian tissue types. The data \textcolorblackwere produced by Su et al. [50] and \textcolorblackthey hold gene expression from human and mouse samples across a diverse array of tissues, organs and cell lines. There are $n=102$ samples and $p=5565$ genes in \textcolorblackthese data. The tissue types have four categories that are known and the goal is to identify the sub-populations that correspond to the four tissue types. \textcolorblackThese data can be publicly sourced from the R-package FABIA [30].

In \textcolorblackFigure 2, we plot the distribution of $\psi_{j}$ (left) \textcolorblackalong with the rank ordering of the 5565 features with respect to the scores $S_{j}$ (right). Using the data driven selection procedure, COSCI selects the top \textcolorblackfour features and returns an estimate of $\alpha_{0}$ almost close to $0.4$ which agrees well with the threshold choice prescribed in Table 10 for sample size 100 and Gaussian noise coordinates. CERs for the aforementioned classification methods are reported in Table 5.

Cardio data: This data has $n=63$ subjects of which 19 are cardiovascular patients and the rest are healthy controls. The genetic expression of each subject has been recorded for $p=\mbox{20,426}$ genes. The goal is to classify the subjects as healthy or cardiovascular based on expression levels of the $p$ genes. The dataset is publicly available on Brad Efron’s webpage. In \textcolorblackFigure 4 (A.2.3), we plot the distribution of $\psi_{j}$ (left) \textcolorblackalong with the rank ordering of the 20,426 features with respect to the scores $S_{j}$ (right). Using the data driven selection procedure, COSCI selects the top 33 features and results in an estimate of $\alpha_{0}$ equal to $0.428$. Due to the small sample size of this data, the data driven procedure returns an unrealistic estimate of null proportion $\hat{\pi}_{0}=0.99$ when the empirical null distribution is estimated on $\mathcal{A}=\{\psi_{(1)}\leq\cdots\leq\psi_{([0.9p])}\}$. We comment more on this observation after introducing the results in Table 5.

RNASeq data: The RNASeq data discussed in Section 1 is an example where both $n$ and $p$ are large. For such large values of $n$, one would expect an asymptotic regime to kick in and hope to see the distribution of $\psi_{j}$’s concentrated around a small value. In \textcolorblackFigure 5 (Appendix A.2.3) (left), we see that approximately $70\%$ of the coordinates have $\psi_{j}\leq 0.1$. Moreover, the distribution of $S_{j}$ (figure 5 right) has an explicit “elbow” at $0.1$, suggesting that we could take the cut-off on $\alpha_{0}$ to be 0.1 and select those features with $S_{j}\geq 0.1$. This is indeed one of the ways that we can approach the problem given this data and that would have left us with \textcolorblackapproximately $2800$ features. We, however, continue with the data driven approach but take $\mathcal{A}=\{\psi_{(1)}\leq\cdots\leq\psi_{([0.7p])}\}$. This leads to an estimated null proportion $\hat{\pi}_{0}=0.70$ with the top 2304 features as the selected ones. Estimated $\hat{\alpha}_{0}$ is 0.188 which (i) is not far from the “elbow” in the distribution of $S_{j}$ and (ii) once again agrees well with the theoretical choice prescribed under Gaussian noises in Table 10.

We make several comments on the CER’s reported in Table 5. The general theme of those results is compelling. These examples suggest that feature screening by COSCI can potentially lead to improvements in clustering error rates even when the underlying clustering algorithm is vanilla $k$-means. On the Multi-tissue data, IF-PCA and Ex. Mass fail to select any feature whereas the \textcolorblackfour COSCI screened features overlaid with IF-PCA clustering clearly demonstrates substantial improvement over the competing error rates. A similar observation follows when COSCI is overlaid with Sparse $k$-means and $k$-means. On the RNASeq data, $k$-means performs well but in this example too, COSCI screening leads to an overall improvement in the error rate with much fewer features. On the Cardio data, performance enhancement is not observed primarily due to the small sample size of this dataset. For such low values of $n$, the added dispersion in the empirical null often masks relatively weaker signals from being identified. In these cases, the empirical null distribution may be estimated on a slightly smaller set $\mathcal{A}=\{\psi_{(1)}\leq\cdots\leq\psi_{([0.85p])}\}$ to improve clustering accuracy at the expense of a few more features. For example, on the Cardio data with $\mathcal{A}$ as prescribed above, COSCI selects only 129 features but returns a comparable CER of 0.486 using $k$-means.

4.3 Beyond marginal screening: Two-way interactions

In this section we present an extension of COSCI that can successfully identify pairs of features that hold cluster information jointly but are un-informative marginally. We expect that these pairs of features will reveal their inherent cluster strength through a suitable linear combination of the form ${X}^{i,j}u$ where ${X}^{i,j}_{n\times 2}$ is the feature pair using the $(i,j)$ feature in ${X}$ and $u\in\mathbb{R}^{2}$ with $||u||_{2}=1$. To determine the optimal $u$ for each feature pair $(i,j)$, we use a grid search in $\mathbb{R}^{2}$ as follows.

Generate a uniform grid of $m$ points $u_{k}$ on the unit circle in $\mathbb{R}^{2}$. 2. 2.

For the feature pair $(i,j)$, define ${Y}_{n\times m}=\{{X}^{i,j}\,u_{k}:1\leq k\leq m\}$ and use Algorithm 1 to get $S_{(i,j)}(u_{1}),\ldots,S_{(i,j)}(u_{m})$ for the $m$ features in ${Y}$. 3. 3.

Obtain the feature score for the pair $(i,j)$ as $S_{(i,j)}=\max_{1\leq k\leq m}S_{(i,j)}(u_{k})$ and the optimal $u$ as $u^{*}_{(i,j)}=\arg\max_{1\leq k\leq m}$ $S_{(i,j)}(u_{k})$.

We repeat this procedure for all the \textcolorblack$p(p-1)/2$ feature pairs and choose

[TABLE]

Often a component of $u^{*}_{(i,j)}$ will be small indicating that one of the contributing features in the pair $(i,j)$ dominates the other in terms of cluster strength. In these scenarios, we may run into the problem of including a lot of redundant pairs in $\mathcal{\hat{I}}_{S}$ due to the strong effect of only one of the features. Such issues are easily resolved by using a naive thresholding rule on $u^{*}_{(i,j)}$ to select the features and modifying $\mathcal{\hat{I}}_{S}$ accordingly:

[TABLE]

To demonstrate the effectiveness of this approach, we consider a simple simulation setting (Experiment V) with $m=20$, $p_{S}=4$, $p_{N}=21$, and $n=2000$. We let the $p_{N}$ noise features be iid Standard Normal, \textcolorblack$\mathcal{N}(0,1)$. For the signal coordinates, we take

$(X_{1},X_{2})\sim\sum_{i=1}^{2}\mathcal{N}({\mu}_{i},\Sigma)/2$ 2. 2.

$X_{3}\sim 0.5\,\mathcal{B}(4,6)+0.5\,\mathcal{B}(7,3)$ 3. 3.

$X_{4}\sim 0.5\,\mathcal{LN}(0.2,0.35)+0.5\,\mathcal{N}(4,0.5)$, where

[TABLE]

In this setting features $X_{3}$ and $X_{4}$ are bi-modal but $(X_{1},X_{2})$ are only jointly bi-modal. Note that the effective dimensionality of the data in this example is \textcolorblack$p_{S}+p_{N}+m(p_{S}+p_{N})(p_{S}+p_{N}-1)/2=6025$. In \textcolorblackTable 6, we report the False Positive and False Negative proportions produced by applying the aforementioned extension of our proposed COSCI algorithm. It successfully identifies all the signal coordinates. As expected, higher values of the threshold parameter improve the FP rate, and the benefit of using the data driven approach to select the features is evident. Once again, the prescribed theoretical choice (Table 10), for the Gaussian noises (with $n=2000$) agrees with these results. Table 7 of the Appendix shows that the performance of the COSCI procedure is not affected even when the noises are significantly correlated with the features.

5 Discussion

We propose COSCI, a novel feature screening method for large scale cluster analysis problems that are characterized by both large sample sizes and high dimensionality of the observations. COSCI efficiently ranks the candidate features in a \textcolorblacknonparametric fashion and, under mild regularity conditions, is robust to the distributional form of the true noise coordinates. We establish theoretical results supporting ideal feature screening properties of our proposed procedure and provide a data driven approach for selecting the screening threshold parameter. Extensive simulation experiments and real data studies demonstrate encouraging performance of our proposed approach.

An interesting topic for future research is extending our marginal screening method by means of utilizing multivariate objective criteria, which are more potent in detecting multivariate cluster information among marginally unimodal features. Preliminary analysis of the corresponding $\ell_{2}$ fusion penalty based criterion, which, unlike the $\ell_{1}$ based approach used in this paper, is non-separable across dimensions, suggests that this criterion can provide a way to move beyond marginal screening.

Appendix A

A.1 Regularity conditions and proofs

We write $g_{j}$ for the marginal density of the standardized \textcolorblack$j$th feature, $(X_{j}-\mathrm{E}X_{j})/\mathrm{SD}(X_{j})$, and let $q_{j}$ denote the corresponding quantile function. Given feature $j$, an interval $(\ell,r)$ and a point $a\in(\ell,r)$, we define the corresponding population criterion function as $G^{j}_{\ell,r}(a)=\mu^{j}_{a,R}-\mu^{j}_{L,a}$, where $\mu^{j}_{\ell,r}$ is the conditional mean on $(\ell,r)$ under the population distribution of the standardized \textcolorblack$j$th feature. Let $P_{nj}$ denote the empirical measure associated with the observations in the \textcolorblack$j$th feature and let $P_{j}$ be the corresponding population distribution. We use the term population cluster to refer to all intervals that appear along the path of the corresponding univariate population splitting procedure, as defined in [43].

We now formally state the regularity conditions needed in Theorem 1. We suppose that there exist positive universal constants $C_{\tau}$ and $c_{\tau}$, which depend only on $\tau$, and a positive universal $C$, such that for all $j\in\mathcal{I}_{N}$ and $\tau\in(0,1)$:

1. C1:

The densities $g_{j}$ are differentiable and unimodal. Also,

[TABLE] 2. C2:

For each population cluster $(\ell,r)$ of the \textcolorblack$j$th feature, satisfying inequality $\int_{L}^{R}g_{j}(x)dx\geq 0.49$, and each $t\in\{L,R\}\cap\arg\max G_{\ell,r}^{j}\cap[q_{j}(\tau),q_{j}(1-\tau)]$, we have $|G^{j\,\prime}_{\ell,r}(t)|>c_{\tau}$.

Condition C1 ensures that we have uniform control over the noise densities. Condition C2 is an appropriate adaptation of a standard regularity condition in M-estimation. As we mentioned earlier, C1 and C2 are satisfied for location-scale families of unimodal differentiable densities with a finite first moment. The assumption that $g_{j}$ are differentiable can be slightly relaxed. However, we prefer to keep this assumption, as it simplifies the presentation of the results.

A.1.1 Proof of Theorem 1 and Corollary 1

We will use “$\lesssim$” to mean that inequality “$\leq$” holds when the \textcolorblackright-hand side is multiplied by a positive constant, which is chosen independently from the parameters $p$, $n$, $j$ and $c_{3}$ (note that the constant is allowed to depend on $\tau$ and $\epsilon_{1}$). To simplify the exposition, we will write $p$ for $p_{N}$, and have index $j$ always correspond to the (noise) coordinates in $\mathcal{I}_{N}$.

Note that the COSCI approach is invariant to linear transformations. We will assume throughout the proof, without loss of generality, that the underlying distribution of each feature has mean \textcolorblack[math] and variance \textcolorblack$1$. As a consequence of standardization, we have

[TABLE]

Let $L_{j}$ and $R_{j}$ denote the smallest and the largest value, respectively, of the (univariate) cluster formed in the \textcolorblack$j$th coordinate after the merge corresponding to $S_{j}(\tau)$. For the same merge, let $a_{j}$ be the midpoint between the closest representatives of the two sub-clusters. We define the empirical criterion functions as $\hat{G}^{j}_{\ell,r}(a)=\hat{\mu}^{j}_{a,R}-\hat{\mu}^{j}_{L,a}$, where $\hat{\mu}_{\ell,r}$ denotes the average of the observations on the \textcolorblack$j$th feature that fall in $[\ell,r]$. We need the following lemma, which is proved in A.1.2.

Lemma 1**.**

For every positive $\epsilon_{1}$ and $\tau$, there exist positive constants $c_{0}$, $c_{3}$, $c_{4}$, $c_{5}$ and $\kappa$, such that, as long as $p\leq\exp(\kappa n)$, the following inequalities simultaneously hold for all $j$ with probability bounded below by $1-c_{4}p^{-c_{5}}$:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

For the remainder of the proof we restrict our attention to the set on which inequalities (8)–(12) are valid. It follows directly from Proposition 2 in [43] that $a_{j}\in\arg\max\hat{G}^{j}_{L_{j},R_{j}}$. Taking $\epsilon_{1}=c_{0}/2$ in Lemma 1, we derive

[TABLE]

Hence, for all $j$, we have

[TABLE]

Because the densities are uniformly bounded, the left-hand side can be replaced by $P_{j}(R_{j}-a_{j})\wedge P_{j}(a_{j}-L_{j})$ at the cost of an additional universal multiplicative factor on the right-hand side. To complete the proof, it is only left to replace each $P_{j}$ with the corresponding empirical probability, $P_{nj}$. This replacement is justified (again, at the cost of an additional universal multiplicative factor) by applying inequalities (10) and (11).

A.1.2 Proof of Lemma 1

We will focus on inequality (8) and, for concreteness, suppose that $a_{j}\in[(L_{j}+R_{j})/2,R_{j}]$. The rest of the cases can be handled using analogous arguments.

Applying Theorem 2.14.9 in [52], together with the union bound, we note that, for each positive $\epsilon$, we have

[TABLE]

for some positive constants $a_{\ell}$ that only depend on $\epsilon$. Because of the assumption on the magnitude of $p$, the \textcolorblackright-hand side in the above display can be bounded \textcolorblackabove by $a_{1}p^{-a_{4}}$. Thus, taking into account condition (5), as well as the definition of $a_{j}$, we can conclude that, with the exception of the corresponding set of small probability, all $a_{j}$ are uniformly bounded. Moreover, applying Theorem 2.14.9 in [52] again and taking advantage of (7), (4) and (6), we can deduce that all $L_{j}$ and $R_{j}$ are restricted to a uniformly chosen bounded interval, on which all $\|g_{j}\|_{\infty}$ and $\|1/g_{j}\|_{\infty}$ are bounded.

From now on, we focus on all the triples $L\leq a\leq R$ in the aforementioned bounded interval, for which $\hat{G}_{\ell,r}^{j}(a)$ is well defined. Note that for all such triples, and all $j$, we have $R-a\lesssim P_{j}(a,R)$ and $P_{j}(a,R)\lesssim R-a$. Define $d_{\ell,r}(x)=(r-x)(r-\ell)^{-1}\mathbf{1}({\ell<x<r})$ and note that

[TABLE]

Let $h_{\ell,r}(x)=x\mathbf{1}(\ell<x<r)$. Observe that $|P_{nj}h_{a,R}-P_{j}h_{a,R}|\lesssim E_{2}$, and define $D_{n}=\max_{j}\sup_{\ell<r}(|P_{nj}h_{\ell,r}-P_{j}h_{\ell,r}|)+(|P_{nj}(\ell,r)-P_{j}(\ell,r)|)$. It follows that

[TABLE]

Applying Theorem 2.14.9 in [52], together with the union bound, we note that, for each positive $\epsilon$, we have

[TABLE]

for some positive constants $a_{l}$ that only depend on $\epsilon$. Because of the assumption on the magnitude of $p$, the \textcolorblackright-hand side in the above display can be bounded \textcolorblackabove by $a_{1}p^{-a_{4}}$. Thus, to establish the bound in inequality (8) we only need to verify that it holds for $E_{1}$ and $E_{2}$, uniformly over $L\leq a\leq R$ in the aforementioned bounded interval. We will focus on $E_{2}$, as $E_{1}$ can be handled with only minor modifications to the argument. We need to show that there exist positive constants $c_{2}$ and $c_{3}$ and a sequence of random variables $M_{n}$, such that $\mathrm{Pr}\,(M_{n}>c_{3})\lesssim p^{-c2}$ and inequalities

[TABLE]

hold for all $j$ and $(\ell,r)$ contained within the bounded interval. Let $M^{j}_{n}$ be the infimum of all those values for which equation (13) holds in the \textcolorblack$j$th coordinate, and define $M_{n}=\max_{j}M^{j}_{n}$. Recall that we have restricted our attention to a uniformly bounded interval, on which $c_{4}=\max_{j}\|1/g_{j}\|_{\infty}$ is positive. Let $\xi_{n}=\ln(p\vee n)/n$, and write $A_{k,n}^{j}$ for the set of intervals $(\ell,r)$ that lie inside the aforementioned uniformly bounded interval and satisfy inequalities $(2^{k-1}-1)\xi_{n}<P_{j}(\ell,r)\leq 2^{k}\xi_{n}$. In what follows, constants $c_{\ell}$ are positive and can be chosen independently from $c_{3}$, $p$, $n$, and $j$. Observe that

[TABLE]

We will bound each summand in the last expression by applying Theorem 2.14.25 in [52] with $\mu_{n}=2^{k/2}\xi_{n}\sqrt{n}$ and $\sigma^{2}_{\mathcal{F}}=2^{k}\xi_{n}$. Note that $\mu_{n}$ is an upper bound on $\sup_{A_{k,n}^{j}}\sqrt{n}\,|P_{nj}(\ell,r)-P_{j}(\ell,r)|$, up to some universal multiplicative factors, as demonstrated in the proof of Lemma 5 in [43]. It follows that

[TABLE]

where $c_{6}$ can be chosen to be positive as long as we take $c_{5}>1/c_{3}$.

We complete the proof by noting that the remaining bound, (12), is implied by conditions (4) and C2, the derivations in the proof of Theorem 1 in [43], and the aforementioned fact that all $L_{j}$ and $R_{j}$ lie in a bounded interval with high probability.

The result of Corollary 1 follows from our Theorem 1 and Theorem 1 in [43], due to the fact that the clustering scores, $S_{j}$, of the signal features are bounded away from zero. $\Box$

A.2 Further empirical results

A.2.1 Correlated data

Here, we assess the screening performance of COSCI when the data has significant correlation between the features. Consider the setup where pairs of features hold cluster information jointly but are un-informative marginally. The design of this experiment is similar to Experiment V presented in Section 4.3. We set $m=20$, $p_{S}=4$, $p_{N}=21$ and $n=2000$. We let the $p_{N}$ noise features be $\mathcal{N}(0,1)$. For the signal coordinates, we take

$(X_{1},X_{2})\sim\sum_{i=1}^{2}\mathcal{N}({\mu}_{i},\Sigma)/2$ 2. 2.

$X_{3}\sim 0.5\,\mathcal{\beta}(4,6)+0.5\,\mathcal{\beta}(7,3)$ 3. 3.

$X_{4}\sim 0.5\,\mathcal{LN}(0.2,0.35)+0.5\,\mathcal{N}(4,0.5)$, where

[TABLE]

In this setting features $X_{3}$ and $X_{4}$ are bi-modal but $(X_{1},X_{2})$ are only jointly bi-modal. Note that the effective dimensionality of the data in this example is \textcolorblack$p_{S}+p_{N}+m(p_{S}+p_{N})(p_{S}+p_{N}-1)/2=6025$. We introduce dependence between the signal features $(X_{1},X_{2})$ and the Gaussian noise features $(X_{5},\ldots,X_{14})$ through a Gaussian copula with a correlation matrix that has all off-diagonal elements equal to $0.9$. To introduce dependence between the signal features $(X_{3},X_{4})$ and the Gaussian noise features $(X_{15},\ldots,X_{25})$, we use a multivariate $t$ copula with $2$ degrees of freedom and with a correlation matrix with all off-diagonal elements equal to $0.8$. \textcolorblack(For additional information about copulas, see, e.g., [24, 25]). Figure 3 presents the distribution of observed average linear correlation between the $p_{S}$ signal features and the Gaussian noise features. The average is across $10$ repetitions of the simulated data for $n=2000$.

We report the False Positive and False Negative proportions for the aforementioned simulation study in \textcolorblackTable 7. COSCI successfully identifies all the signal coordinates. Comparing with \textcolorblackTable 6, we see that its performance is not affected due to correlation in the data. As expected, higher values of the threshold parameter improve the FP rate, and the benefit of using the data driven approach to select the features is evident. We find that correlation among features does not appear to impact the screening performance of COSCI.

A.2.2 Clustering errors

In the simulation experiments of Section 4.1 we evaluated the feature selection performance of COSCI and obtained encouraging results. Here, we check whether COSCI’s better feature selection performance also lead to a reduction of the clustering error rates in those experiments.

In \textcolorblackTable 8, we present the clustering errors for Simulation Experiments I and II. Recall that for Simulation Experiment I, the $p_{S}=5$ signal features each had two clusters and for Simulation Experiment II, the first five of the $p_{S}=6$ signal features had two clusters and the sixth one had thee clusters. Thus, the true number of clusters in those experiments is $32$ and $96$ respectively. As the true number of clusters are so large, traditional clustering dissimilarity measures (such as the classification error rates (CER) used in Witten and Tibshirani [56]) fail to distinguish between good and bad clustering if we compare with respect to all the product clusters. Here, for each signal coordinate we calculate the CER considering the true cluster labeling to be based on that signal only and report the average CER across all the signal coordinates as the clustering error. Each simulation experiment was repeated $50$ times and the average clustering error over the $50$ repetitions is reported in the tables. For Ex. Mass and COSCI, we use \textcolorblack$k$-means to perform clustering on the selected features.

In \textcolorblackTable 8, it was seen that COSCI coupled with the data driven approach continues to provide the best CER amongst all the competing methods. IF-PCA, SAS, SpKM and SpHC exhibit very high CER in these settings even though their screening performance, especially that of SAS and IF-PCA, was competitive (see Tables 1–2 in Section 4.1). COSCI and Ex. Mass return far better CER and as is expected, with $n$ small and increasing $\alpha_{0}$, COSCI screening begins to miss the signal features and therefore returns a higher CER. Similar characteristics in the clustering efficacies were also observed across the other simulation scenarios of Section 4.1.

A.2.3 COSCI Analysis on Cardio and RNA Seq Datasets

A.2.4 More real data examples

In Table 9, we present more real data examples to assess the performance of COSCI. We follow the theme presented in \textcolorblackSection 4.2 and overlay COSCI with $k$-means, Sparse $k$-means and IF-PCA for each of the eleven datasets considered below. The first ten data sets are sourced from Jiashun Jin’s webpage (see: http://www.stat.cmu.edu/~jiashun/Research/software/GenomicsData/ and Table 1 in [35] for more information on these data sets.). The last data set is available on Brad Efron’s webpage (see: http://statweb.stanford.edu/~ckirby/brad/LSI/datasets-and-programs/datasets.html).

IF-PCA returns a significantly smaller error rate on the Lymphoma, Leukemia and Lung Cancer datasets but barring these four datasets (Lymphoma, Leukemia, Lung Cancer (1) and Lung Cancer (2)), feature screening by COSCI returns comparable and often better error rates using fewer features when compared to the best competing method. In the datasets of Breast Cancer and Michigan COSCI selects more features than the best competing method but returns a comparable error rate. Detailed introspection on the four datasets where IF-PCA performs better than COSCI reveals that the overlap between these two methods in the screening step is minimal here. COSCI fails to pick up those signal features that pass through its merge size filter and are statistically far away from a Gaussian distribution. \textcolorblackIn contrast, IF-PCA precisely picks up these features which happen to be the best features to perform clustering in these four datasets. In simulation Experiments III–IV and in the real datasets considered in Section 4.2, the setting is, however, exactly opposite where the best signal features are not uni-modal and the noise features are not always Normally distributed. In these settings, COSCI prospers and by far returns the best screening performance than all the other competing methods considered.

A.3 Algorithm details

A.3.1 Detailed description of Algorithm 1

A.3.2 Estimation of Hyperparameters — simulation based

A.3.3 Two-stage approach to signal screening

We briefly describe the \textcolorblacktwo-stage approach to signal screening introduced in [9]. We start with Eq.  (3) where the estimated fdr has been obtained as

[TABLE]

Order the estimated fdr’s $T_{j}$ from smallest to largest so that $T_{(1)}\leq\cdots\leq T_{(p)}$. In the first stage, we estimate \textcolorblackStage 1 screening cutoff $k_{s}$ by

[TABLE]

with $\delta_{p}=\min\{p,(\ln p)^{-1}\}$ and select the features as

[TABLE]

In the second stage, the so-called Discovery stage, we take the \textcolorblackStage 1 screened fdr’s $l\{T_{j}:j\in\mathcal{S}_{s}\}$ and estimate the \textcolorblackStage 2 screening cutoff $k_{d}$ by

[TABLE]

and select the features as

[TABLE]

Proposition 2 in [9] guarantees that in \textcolorblackStage 1, $\mathcal{S}_{s}$ is the largest subset such that the missed discovery rate (MDR) is controlled at level $\delta_{p}$ while in \textcolorblackStage 2, $\mathcal{\hat{I}}_{S}$ is the smallest subset such that the false positive rate (FPR) is controlled at level $\delta_{p}$.

A.3.4 Analysis of the population procedure for Gaussian mixtures

In Table 11, which is an adapted version of Table 1 in the supplementary material for [43], we document the behavior of the population clustering procedure for a wide variety of mixtures of two Gaussian distributions on the real line. For \textcolorblackseven different levels of separation between the two means we consider \textcolorblacknine different mixing proportions, from the symmetric case of 50:50 mixing to the highly skewed 10:90 mixing. The behavior of the population splitting procedure in other cases can be interpolated using continuity arguments. In all of the scenarios, the population procedure identifies either two clusters or one. The latter happens only in the settings where the separation between the sub-populations, or the size of one sub-population, is very small. We also present the location of the split point, $s^{*}$ (“NO” denotes the cases where no splits are detected, and, thus, only one cluster is identified), the local minimum of the density, $\operatorname*{\underline{m}}$, and the split point minimizing the expected misclassification error, $s_{\text{MC}}$. Finally, we report, under Excess MCE, how much the misclassification error of the population clustering procedure exceeds that of the $s_{\text{MC}}$ based oracle rule.

Acknowledgments. We would like to thank the Editor, the Associate Editor and two anonymous referees for many helpful suggestions that improved the paper. We thank Sara Garcia and Gregory Giecold for sharing the RNASeq data and ECLAIR source codes with us. Mukherjee’s research was partially supported by the Zumberge individual award from the University of Southern California’s James H. Zumberge Faculty Research and Innovation Fund.

Supplementary material

The R code and the data sets used in this paper can be downloaded from the following link -

http://github.com/trambakbanerjee/COSCI.