Feature screening for clustering analysis

Suppose that we have $n$ independent observations ${\bf x}_{i}=(x_{i1},\dots,x_{ip})^{\rm T}\in\mathbb{R}^{p}$ ($i=1,\dots,n$) from $G$ clusters and ${\bm{\alpha}}=(\alpha_{1},\dots,\alpha_{G})$ be the proportions of different clusters $(\sum_{g=1}^{G}\alpha_{g}=1,\alpha_{g}>0,g=1,\dots,G)$. We denote the unknown cluster labels as $g_{i}\in\{1,\dots,G\}$ ($i=1,\dots,n$). Assume that given the cluster label $g$, the conditional distribution $F_{j}(x|g)$ of $x_{ij}$ is from a known identifiable parametric distribution family $\mathcal{P}=\{f(x;\bm{\theta}):\bm{\theta}\in\Theta\subset\mathbb{R}^{d}\}$, where $f(x;\bm{\theta})$ is the density function with respect to a $\sigma$-finite measure $\mu$, on $\mathbb{R}$ parameterized by $\bm{\theta}$, and $\Theta\subset\mathbb{R}^{d}$ is a convex compact parameter space. Note that for count data, the measure $\mu$ can be taken as the counting measure of the nonnegative integers; for continuous data, $\mu$ is the Lebesgue measure on $\mathbb{R}$. Thus, our method and theory apply to both count and continuous data. Define $\Xi=\Theta^{G}$ as the product space of $\Theta$.

In high dimensional clustering problems, only a small portion of the $p$ features contain information about the cluster labels and the majority of them are irrelevant to the sample clusters. Our goal is to screen out the cluster-irrelevant features to facilitate downstream clustering analysis. Intuitively, if the $j$th random variable $x_{j}\in\mathbb{R}$ is unrelated with the cluster label $g$, the conditional distribution $F_{j}(x|g)$ of $x_{j}$ given the cluster label $g$ should be independent of the cluster label $g$, or $F_{j}(x|g=1)=\cdots=F_{j}(x|g=G)$. If, on the other hand, the $j$th random variable $x_{j}\in R$ is a cluster-relevant feature, there are at least two $g\neq g^{\prime}$ such that $F_{j}(x|g)\neq F_{j}(x|g^{\prime})$.

Let $f(x;\bm{\theta}_{jg})$ be the density function of the conditional distribution $F_{j}(x|g)$. The labels are unknown and the random variable $x_{j}$ should follow a mixture distribution $\varphi(x;\bm{\xi}_{j},{\bm{\alpha}})=\sum_{g=1}^{G}\alpha_{g}f(x;\bm{\theta}_{jg})$, where $\bm{\xi}_{j}=(\bm{\theta}_{j1}^{\rm T},\dots,\bm{\theta}_{jG}^{\rm T})^{\rm T}\in\Xi=\Theta^{G}$. Define the interior of the $G-1$ dimensional probability simplex as $\mathbb{S}^{G-1}=\{{\bm{\alpha}}\in\mathbb{R}^{G}:\sum_{g=1}^{G}\alpha_{g}=1,\alpha_{g}>0,\mbox{ for }g=1,\dots,G\}$ and the $G$-mixture distribution family as

$$\mathcal{P}^{G}=\left\{\sum_{g=1}^{G}\alpha_{g}f(x;\bm{\theta}_{g}):{\bm{\alpha}}\in\mathbb{S}^{G-1},\bm{\theta}_{g}\in\Theta,g=1,\dots,G\right\}.$$

We have $\varphi(x;\bm{\xi}_{j},{\bm{\alpha}})\in\mathcal{P}^{G}$. In this paper, we assume that $\mathcal{P}^{G}$ is an identifiable finite mixture, in other words, $\mathcal{P}$ is a linearly independent set over the field of real numbers (Yakowitz and Spragins, 1968). For a cluster-irrelevant feature $j$, $\bm{\theta}_{j1}=\cdots=\bm{\theta}_{jG}$ and thus $\varphi(x;\bm{\xi}_{j},{\bm{\alpha}})\in\mathcal{P}$. For a cluster-relevant feature $j$, there are at least two $g\neq g^{\prime}$ such that $\bm{\theta}_{jg}\neq\bm{\theta}_{jg^{\prime}}$ and $\varphi(x;\bm{\xi}_{j},{\bm{\alpha}})\in\mathcal{P}^{G}\backslash\mathcal{P}$. Therefore, we can consider the following hypothesis testing problems to screen for the cluster-relevant features.

$$\mathbb{H}_{j0}:\varphi(x;\bm{\xi}_{j},{\bm{\alpha}})\in\mathcal{P}\mbox{ v.s. }\mathbb{H}_{j1}:\varphi(x;\bm{\xi}_{j},{\bm{\alpha}})\in\mathcal{P}^{G}\backslash\mathcal{P}.$$

(1)

We call the models under the null hypotheses $\mathbb{H}_{j0}$ homogeneous models, and those under the alternative hypotheses $\mathbb{H}_{j1}$ heterogeneous models. In real applications, the number of clusters $G$ is often unknown. However, we often can have a rough estimate of $G$ and can choose $G$ to be larger than the true number of clusters. In such cases, the null and alternative hypotheses still hold for the cluster-irrelevant and relevant features, respectively. Simulation shows that the choice of $G$ has little influence on the performance of EM-test, especially when $G$ is chosen to be larger than the true clusters (Supplementary Section D).

We use the EM-test statistic for feature screening of clustering analyses. Theoretical results of the EM-test statistic are developed for the hypothesis testing problem (1) under general settings with multiple parameters ($d\geq 1$), multiple components ($G\geq 2$) and both continuous and count data. Let ${\bf x}=(x_{1},\dots,x_{n})$ be a random sample of size $n$ from a $G$-mixture model

$$\varphi(x;\bm{\xi},{\bm{\alpha}})=\sum_{g=1}^{G}\alpha_{g}f(x;\bm{\theta}_{g}),$$

(2)

where $\bm{\theta}_{g}\in\Theta$, ($g=1,\ldots,G$), $\bm{\xi}=(\bm{\theta}_{1},\ldots,\bm{\theta}_{G})$ and ${\bm{\alpha}}=(\alpha_{1},\dots,\alpha_{G})$. Let $l_{n}(\bm{\xi},{\bm{\alpha}})=\sum_{i=1}^{n}\hbox{log}\ \varphi(x_{i};\bm{\xi},{\bm{\alpha}})$ be the log-likelihood function, and define the penalized log-likelihood function as

$$pl_{n}(\bm{\xi},{\bm{\alpha}})=\sum_{i=1}^{n}\hbox{log}\ \varphi(x_{i};\bm{\xi},{\bm{\alpha}})+p({\bm{\alpha}}),$$

(3)

where $p({\bm{\alpha}})=\lambda\left(\sum_{g=1}^{G}{\rm log}\alpha_{g}+G\hbox{log}(G)\right)$ is a penalty function, where $\lambda>0$ is a penalty parameter and is always set as 0.00001 in the simulation and real data analyses of this paper. Simulation shows that EM-test is robust to the choice of the penalty parameter $\lambda$ and $\lambda=0.00001$ gives very similar results to other choices of $\lambda$ (Supplementary Table S3). Largely speaking, the EM-test statistic is defined as the difference between the maximum penalized log-likelihoods of the heterogeneous and homogeneous models. The maximum penalized log-likelihood under the heterogeneous model is obtained using the EM algorithm. More specifically, we use the following procedure to calculate the EM-test statistic.

Suppose that $\hat{\bm{\xi}}_{0}=\left(\hat{\bm{\theta}}_{0},\ldots,\hat{\bm{\theta}}_{0}\right)$ is the estimator that maximizes the penalized log-likelihood function (3) under the homogeneous model. Under the heterogeneous model, given any initial value ${\bm{\alpha}}^{(0)}\in\mathbb{S}^{G-1}$ we first compute

$$\bm{\xi}^{(0)}=\mathop{\arg\max}_{\begin{subarray}{c}\bm{\xi}\in\Xi\end{subarray}}\sum_{i=1}^{n}\hbox{log}\ \varphi\left(x_{i};\bm{\xi},{\bm{\alpha}}^{(0)}\right)+p\left({\bm{\alpha}}^{(0)}\right).$$

(4)

Assume that ${\bm{\alpha}}^{(k)}$ and $\bm{\xi}^{(k)}$ are the estimators at the $k$-th iteration of the EM algorithm. The E-step updates the posterior probability of the $i$-th sample coming from the $g$-th component by

$$w_{gi}^{(k)}=\frac{\alpha^{(k)}_{g}f\left(x_{i};\bm{\theta}_{g}^{(k)}\right)}{\varphi\left(x_{i};\bm{\xi}^{(k)},{\bm{\alpha}}^{(k)}\right)}.$$

(5)

At the $k+1$-th iteration, the M-step updates ${\bm{\alpha}}$ and $\bm{\xi}$ such that

$${\bm{\alpha}}^{(k+1)}=\mathop{\arg\max}_{{{\bm{\alpha}}}\in\mathbb{S}^{G-1}}\sum_{g=1}^{G}\sum_{i=1}^{n}w_{gi}^{(k)}\hbox{log}(\alpha_{g})+p({\bm{\alpha}}),~{}\mbox{and}$$

(6)

$$\bm{\xi}^{(k+1)}=\mathop{\arg\max}_{{\bm{\xi}\in\Xi}}\sum_{g=1}^{G}\sum_{i=1}^{n}w_{gi}^{(k)}\hbox{log}f(x_{i};\bm{\theta}_{g}).$$

(7)

Let $K>0$ be the maximum number of EM updates. We define ${M}_{n}^{(K)}({\bm{\alpha}}^{(0)})=2\{pl_{n}(\bm{\xi}^{(K)},{\bm{\alpha}}^{(K)})-pl_{n}(\hat{\bm{\xi}}_{0},{\bm{\alpha}}_{0})\},$ where ${\bm{\alpha}}_{0}=(1/G,\dots,1/G)^{\rm T}$. To improve the performance, we choose a set of initial values $\{{\bm{\alpha}}_{1},\dots,{\bm{\alpha}}_{T}\}$ and define the EM-test statistic as ${\rm EM}_{n}^{(K)}=\max\{M_{n}^{(K)}({\bm{\alpha}}_{t}),t=1,\dots,T\}.$ Intuitively, under the homogeneous model, $\bm{\xi}^{(K)}$ and $\hat{\bm{\xi}}_{0}$ are close to $\bm{\xi}_{0}$ and are close to each other, while under the heterogeneous model, $\bm{\xi}^{(K)}$ and $\hat{\bm{\xi}}_{0}$ are far away from each other. Hence, we reject the null hypothesis in (1) when ${\rm EM}_{n}^{(K)}$ is large. In this paper, we always assume that $K\geq 3$ is a fixed number. In simulation and real data analyses, we set $K=100$. Simulation shows that EM-test is robust to the choice of $K$. When $K$ is chosen too large, EM-test tends to have slightly more false positives and $K=100$ is a reasonable choice (Supplementary Section D).

With the EM-test statistic, we can use the following procedure to screen for the cluster-relevant features. Let ${\rm EM}_{nj}^{(K)}$ be the EM-test statistic corresponding to the $j$th hypothesis testing problem (1). Given a threshold $t_{n}>0$, if ${\rm EM}_{nj}^{(K)}<t_{n}$, we will screen out the $j$th feature; Otherwise, we will retain the $j$th feature as a cluster-relevant feature. Theoretical results in Section 3 show that if we choose $t_{n}=n^{\vartheta}$ ($0<\vartheta<1$), this feature screening procedure can have the sure independent screening property or even the consistency of selection property.

We use ${\rm diam}(\Xi)$ to represent the Euclidean diameter of $\Xi$. Denote by $|\cdot|$ the absolute value of a real number or cardinality of a set. For two sequences of random variable $\left\{a_{n}\right\}_{n=1}^{\infty}$ and $\left\{b_{n}\right\}_{n=1}^{\infty}$, we write $a_{n}=o_{p}\left(|b_{n}|\right)$ if $a_{n}/|b_{n}|\rightarrow 0$ in probability, and $a_{n}=O_{p}\left(|b_{n}|\right)$ if there exists a positive constant $C$ such that $a_{n}\leq C|b_{n}|$ in probability. For real numbers $a$ and $b$, let $a\wedge b=\min\{a,b\}$ and $a\vee b=\max\{a,b\}$. Define $\left\|{\bf a}\right\|_{2}=\sqrt{\sum_{i=1}^{n}a_{i}^{2}}$ as the $L_{2}$-norm of the vector ${\bf a}=(a_{1},\dots,a_{n})^{\rm T}\in\mathbb{R}^{n}$, and ${\rm vech}({\bf A})=\left(a_{11},a_{22},\ldots,a_{dd},a_{12},\ldots,a_{1d},\ldots,a_{d(d-1)}\right)^{\rm T}\in\mathbb{R}^{d(d+1)/2}$ as the vectorization of the symmetric $d$-dimensional matrix ${\bf A}=(a_{ij})$. We use ${\bf A}\succeq{\bf 0}$ to represent that the matrix ${\bf A}$ is positive semi-definite. For a random variable $X$, we define its sub-exponential norm as $\|X\|_{\psi_{1}}=\inf\{t>0:\mathbb{E}\left(\exp(|X|/t)\right)\leq 2\}$. If $X$ is a random variable from a homogeneous model $\mathbb{H}_{0}$ and $g(x)$ is a function, we define the $L_{m}$-norm of $g(X)$ as $\left\|g(X)\right\|_{L^{m}}=\left({\mathbb{E}}\left[g^{m}(X)\right]\right)^{1/m}$.

We denote ${\bm{\alpha}}^{*}=(\alpha_{1}^{*},\dots,\alpha_{G}^{*})$ as the true proportions of the $G$ clusters and $\bm{\xi}^{*}_{j}=(\bm{\theta}^{*}_{j1},\ldots,\bm{\theta}^{*}_{jG})$ as the true parameters of the mixture model corresponding to the $j$th feature. We assume that ${\bm{\alpha}}^{*}$ is fixed and $\min_{g}\alpha_{g}^{*}>0$. We use $\delta>0$ as a fixed very small constant. If the $j$th feature is from the homogeneous model $\mathbb{H}_{j0}$ (cluster-irrelevant), we write $\bm{\xi}^{*}_{j}=\bm{\xi}_{j0}=(\bm{\theta}_{j0},\ldots,\bm{\theta}_{j0})$ as its true parameters. When appropriate, we drop the subscript $j$ and use $\bm{\xi}^{*}$ as the true parameters of a general mixture model and $\bm{\xi}_{0}=(\bm{\theta}_{0},\ldots,\bm{\theta}_{0})$ as the true parameter of some general homogeneous model $\mathbb{H}_{0}$. We always assume that $\bm{\theta}_{0}$ is an interior point of $\Theta$ and use ${\bm{\alpha}}_{0}=(1/G,\dots,1/G)^{\rm T}$. Let

	$\displaystyle Y_{ih}=\frac{1}{f(x_{i},\bm{\theta}_{0})}\frac{\partial f(x_{i},\bm{\theta}_{0})}{\partial\theta_{h}},~{}Z_{ih}=\frac{1}{2f(x_{i},\bm{\theta}_{0})}\frac{\partial^{2}f(x_{i},\bm{\theta}_{0})}{\partial\theta_{h}^{2}},$
	$\displaystyle U_{ih\ell}=\frac{1}{f(x_{i},\bm{\theta}_{0})}\frac{\partial^{2}f(x_{i},\bm{\theta}_{0})}{\partial\theta_{h}\partial\theta_{\ell}}(h<\ell),~{}{\bf b}_{1i}=(Y_{i1},\ldots,Y_{id})^{\rm T},$		(8)
	$\displaystyle{\bf b}_{2i}=(Z_{i1},\ldots,Z_{id},U_{i12},\ldots,U_{i(d-1)d})^{\rm T},~{}\mbox{and}~{}{\bf b}_{i}=\left({\bf b}_{1i}^{\rm T},{\bf b}_{2i}^{\rm T}\right)^{\rm T}.$

We need the following regularity conditions before presenting the tail probability bounds of the EM-test statistic under $\mathbb{H}_{0}$.

• (C1)

  For every $\bm{\theta}_{0}\in\bm{\Theta}$ and sufficiently small ball $V\subset\Theta$ around $\bm{\theta}_{0}$, we assume that the function $\sup_{\bm{\theta}\in V}{f^{1/2}(x;\bm{\theta})}{f^{1/2}(x;\bm{\theta}_{0})}$ is measurable and $\int\sup_{\bm{\theta}\in V}{f^{1/2}(x;\bm{\theta})}{f^{1/2}(x;\bm{\theta}_{0})}\mu({\rm d}x)<\infty.$ In addition, for every sufficiently small ball $U\subset\Xi$ around $\bm{\xi}_{0}=(\bm{\theta}_{0},\cdots,\bm{\theta}_{0})$ and ${\bm{\alpha}}\in\mathbb{S}^{G-1}$, we assume that the function $\sup_{\bm{\xi}\in U}\hbox{log}\{\varphi(x;\bm{\xi},{\bm{\alpha}})\}$ is measurable and ${\mathbb{E}}(\sup_{\bm{\xi}\in U}\hbox{log}\{1+\varphi(x;\bm{\xi},{\bm{\alpha}})\})<+\infty.$

• (C2)

  The density function $f(x;\bm{\theta})$ has a common support for $\bm{\theta}\in\Theta$ and continuous 5th order partial derivatives with respect to $\bm{\theta}$.

• (C3)

  Let $m>0$ be an integer and $M>0$ be a constant. There are a function $g(x,\bm{\theta})>0$ with $\sup_{\bm{\theta}\in\Theta}\left\|g(x;{\bm{\theta}})\right\|_{L^{8m}}\leq M$, a function $r(x)>0$ with $\int r^{2}(x)\mu({\rm d}x)\leq M$ and a constant $\tau>0$, such that, for all $\bm{\theta}_{0}\in\Theta$, $h=1,\dots,5$ and ${j_{1}},\ldots,{j_{h}}\in\{1,\dots,d\}$,

  $$\sup_{\|\bm{\theta}-\bm{\theta}_{0}\|_{2}\leq\tau}\left|\frac{\partial^{h}f(x;\bm{\theta})}{\partial\theta_{j_{1}}\cdots\partial\theta_{j_{h}}}\bigg{/}f(x;\bm{\theta}_{0})\right|\leq g(x;{\bm{\theta}_{0}}),\mbox{ and }$$

  $$\sup_{\bm{\theta}\in\Theta}f(x;\bm{\theta})+\sup_{\bm{\theta}\in\Theta}\frac{1}{{f(x;\bm{\theta})}}\left\|\frac{\partial f(x;\bm{\theta})}{\partial\bm{\theta}}\right\|_{2}^{2}\leq r^{2}(x).$$

• (C4)

  The minimum eigenvalue $\lambda_{\rm min}({\bf B}(\bm{\theta}_{0}))$ of the covariance matrix ${\bf B}(\bm{\theta}_{0})=\mbox{cov}({\bf b}_{i})$ satisfies $\lambda_{\rm min}:=\inf_{\bm{\theta}_{0}\in\Theta}\lambda_{\rm min}({\bf B}(\bm{\theta}_{0}))>0.$

Condition (C1) is the Wald consistency condition, which can be founded in Van der Vaart (2000). It also ensures the continuity of the Hellinger distance which we will define below. Condition (C2) guarantees the smoothness of $f(x;\bm{\theta})$. Condition (C3) ia a technical condition on the partial derivatives. It guarantees that there is a dominating function of the remainder term in the Taylor expansion, and thus allows us to give polynomial tail probability bounds of the higher-order infinitesimal terms of the EM-test statistic under $\mathbb{H}_{0}$. Condition (C4) is the strong identifiability condition (Chen, 1995; Nguyen, 2013). Most of the commonly used one-parameter distributions, such as the Poisson distribution and the exponential distribution, satisfy Condition (C3) and (C4). Many multiple-parameter distributions including the negative binomial distribution and the gamma distribution also satisfy Condition (C3) and (C4).

Observe that when $n\rightarrow\infty$, $Cn^{-1/16}{\rm log}^{3/2}n$ approaches to zero. Therefore, roughly speaking, Theorem 2 shows that when $n$ is sufficiently large, under $\mathbb{H}_{0}$, the tail probability of the EM-test statistic greater than $t$ has a polynomial decay rate. To prove Theorem 2, we first derive the tail probability bound for the mixture parameter estimators $\bm{\xi}^{(k)}$ by analyzing the empirical processes indexed by $\bm{\xi}$ (Wong and Shen, 1995). Then, we analyze the Taylor expansion of ${\rm EM}_{n}^{(K)}$ and bound each term in the expansion using concentration inequalities (Wainwright, 2019). Details of the proof are given in the Supplementary material.

Our next goal is to show that the EM-test statistic diverges to infinity under $\mathbb{H}_{1}$ with a high probability. Recall that ${\bm{\alpha}}^{*}$ is the true proportion parameter and that under $\mathbb{H}_{1}$, $\bm{\xi}^{*}\in\Xi_{1}$, where $\Xi_{1}$ is defined in (85). We define ${\bm{\theta}}^{\dagger}_{0}=\mathop{\arg\max}_{\bm{\theta}\in{\Theta}}{\mathbb{E}}_{{\bm{\alpha}}^{*},\bm{\xi}^{*}}\left[\hbox{log}f(x;\bm{\theta})\right]$ as the parameter of the homogeneous model that is closest to the true heterogeneous model in terms of the Kullback-Leibler divergence. Denote $\bm{\xi}_{0}^{\dagger}=\left({\bm{\theta}}^{\dagger}_{0},\dots,{\bm{\theta}}^{\dagger}_{0}\right)$. Similarly, given an initial value ${\bm{\alpha}}^{(0)}$, we can find a heterogeneous model with a proportion parameter ${\bm{\alpha}}^{(0)}$ that is closest to the true heterogeneous model and denote its parameter as ${\bm{\xi}}^{\dagger}=\mathop{\arg\max}_{\bm{\xi}\in\Xi}{\mathbb{E}}_{{\bm{\alpha}}^{*},\bm{\xi}^{*}}\left[\hbox{log}\ \varphi\left(x;\bm{\xi},{\bm{\alpha}}^{(0)}\right)\right].$ Note that ${\bm{\xi}}^{\dagger}$ and ${\bm{\theta}}^{\dagger}_{0}$ depend on the true value ${\bm{\alpha}}^{*},\bm{\xi}^{*}$.

Define $R(x;\bm{\xi}^{*})=\hbox{log}\ \varphi\left(x;{\bm{\xi}}^{\dagger},{\bm{\alpha}}^{(0)}\right)-\hbox{log}\ f\left(x;{\bm{\theta}}_{0}^{\dagger}\right)$ as the difference between two “working” log-likelihood $\hbox{log}\ \varphi\left(x;{\bm{\xi}}^{\dagger},{\bm{\alpha}}^{(0)}\right)$ and $\hbox{log}\ f\left(x;{\bm{\theta}}_{0}^{\dagger}\right)$. If the initial value ${\bm{\alpha}}^{(0)}$ is close to the true proportion ${\bm{\alpha}}^{*}$, the expectation of $R(x;\bm{\xi}^{*})$ would be bounded away from zero. So, the one-step EM-test statistic, and thus the EM-test statistic, would be large. Thus, we would correctly reject the null hypothesis with a high probability. Furthermore, denoting $D(\bm{\theta})={\mathbb{E}}_{{\bm{\alpha}}^{*},\bm{\xi}^{*}}[\hbox{log}f(x;\bm{\theta})]$, we define a mean-zero empirical process indexed by $\bm{\theta}\in\Theta$ as

$$Z_{\bm{\theta}}(\bm{\xi}^{*})=n^{-1/2}\sum_{i=1}^{n}\left\{\hbox{log}\ f(x_{i};\bm{\theta})-\hbox{log}\ f\left(x_{i};{\bm{\theta}}^{\dagger}_{0}\right)-\left[D(\bm{\theta})-D\left({\bm{\theta}}^{\dagger}_{0}\right)\right]\right\}.$$

We need the following conditions under $\mathbb{H}_{1}$.

• (C5)

  The initial value ${\bm{\alpha}}^{(0)}$ fulfills $\varrho=\inf_{\bm{\xi}^{*}\in\Xi_{1}}{\mathbb{E}}_{{\bm{\alpha}}^{*},\bm{\xi}^{*}}\left[R(x;\bm{\xi}^{*})\right]>0.$

• (C6)

  There exists a constant $M_{\psi_{1}}$ such that $\sup_{\bm{\xi}^{*}\in\Xi_{1}}\left\|R(x;\bm{\xi}^{*})-{\mathbb{E}}_{{\bm{\alpha}}^{*},\bm{\xi}^{*}}\left[R(x;\bm{\xi}^{*})\right]\right\|_{\psi_{1}}\leq M_{\psi_{1}}.$

• (C7)

  $\{Z_{\bm{\theta}}(\bm{\xi}^{*}):\bm{\theta}\in\Theta\}$ is a $\psi_{1}$- process such that for any $\bm{\theta},\bm{\theta}^{\prime}\in\Theta$, $\sup_{\bm{\xi}^{*}\in\Xi_{1}}\|Z_{\bm{\theta}}(\bm{\xi}^{*})-Z_{\bm{\theta}^{\prime}}(\bm{\xi}^{*})\|_{\psi_{1}}\leq C_{\rho}\|\bm{\theta}-\bm{\theta}^{\prime}\|_{2},$ where $C_{\rho}>0$ is a constant.

Condition (C5) is a key assumption. Because the EM algorithm cannot guarantee convergence to the global maximum, we need to choose an initial value ${\bm{\alpha}}^{(0)}$ such that the theoretical best heterogeneous model that we can achieve in one step EM update is uniformly closer to the true heterogeneous model than the best homogeneous model. Note that we always have $\varrho\geq 0$, but it is hard to give a necessary and sufficient condition for the choice of ${\bm{\alpha}}^{(0)}$ such that $\varrho>0$. However, we can show that if $\left\|{\bm{\alpha}}^{*}-{\bm{\alpha}}^{(0)}\right\|_{2}\leq\tau(\gamma)$ for some constant $\tau(\gamma)$, Condition (C5) holds (see Section C.1 in Supplementary material for more discussion). Condition (C6) and (C7) are two weaker conditions and hold for many commonly used distribution families.Under these conditions, we obtain the following tail probability bound of the EM-test statistic under $\mathbb{H}_{1}$.

Theorem 3 and Corollary 1 say that under $\mathbb{H}_{1}$, by selecting a suitable threshold, the cluster-relevant feature can be retained with high probability.

In many applications, giving a valid $p$-value of the retained feature is also crucial. In this section, we give the limiting distribution of the EM-test statistic under $\mathbb{H}_{0}$. To derive the limiting distribution, we only need the following weaker conditions in replacement of Condition (C3) and (C4).

• (WC3)

  For all $h=1,\dots,5$ and $\theta_{j_{1}},\ldots,\theta_{j_{h}}$, there exists a function $g(x;\bm{\theta}_{0})\geq 0$ and a constant $\tau>0$ such that

  $$\sup_{\|\bm{\theta}-\bm{\theta}_{0}\|\leq\tau}\left|\frac{\partial^{h}f(x,\bm{\theta})}{\partial\theta_{j_{1}}\cdots\partial\theta_{j_{h}}}\bigg{/}f(x,\bm{\theta}_{0})\right|\leq g(x;\bm{\theta}_{0})$$

  and $\left\|g(x;\bm{\theta}_{0})\right\|_{L^{3}}<\infty$.

• (WC4)

  The covariance matrix ${\bf B}(\bm{\theta}_{0})=\mbox{cov}({\bf b}_{i})$ of ${\bf b}_{i}$ is positive definite.

Let $r=\min(G-1,d)$ and

$$\mathcal{V}=\left\{{\rm vech}({\bf V}):{\bf V}\in\mathbb{R}^{d\times d}\mbox{ is symmetric},~{}{\rm rank}({\bf V})\leq r,{\bf V}\succeq 0\right\}.$$

(11)

For $j,k=1,2$, let ${\bf B}_{jk}={\mathbb{E}}_{\bm{\theta}_{0}}(\{{\bf b}_{ji}-{\mathbb{E}}({\bf b}_{ji})\}\{{\bf b}_{ki}-{\mathbb{E}}({\bf b}_{ki})\}^{\rm T})$ and $\tilde{{\bf b}}_{2i}={\bf b}_{2i}-{\bf B}_{21}{\bf B}_{11}^{-1}{\bf b}_{1i}$. The covariance matrix of $\tilde{{\bf b}}_{2i}$ is $\widetilde{{\bf B}}_{22}={{\bf B}}_{22}-{\bf B}_{21}{\bf B}_{11}^{-1}{\bf B}_{12}$.

If $d=1$, the limiting distribution in Theorem 4 is $0.5\chi^{2}_{1}+0.5\chi^{2}_{0}$, the same as the one in Li et al. (2009), while, if $G=2$, it is the distribution in Niu et al. (2011). When $G>d$, we have $r=d$ and the limiting distribution will be independent of the component number $G$. Generally, it is computationally difficult to calculate the limiting distribution in Theorem 4. When $G>d$, the feasible domain $\mathcal{V}=\left\{{\rm vech}({\bf V}):{\bf V}\in\mathbb{R}^{d\times d},{\bf V}\succeq 0\right\}$ is a positive semi-definite matrix cone. Computation of the limiting distribution in Theorem 4 becomes a classic cone quadratic program and can be solved using the algorithms reviewed in Vandenberghe (2010), but these algorithms are still computationally expensive. However, it can be easily shown that $\sup_{{\bf v}\in\mathcal{V}}2{\bf v}^{\rm T}{\bf{w}}-{\bf v}^{\rm T}\widetilde{{\bf B}}_{22}{\bf v}\leq{\bf w}^{\rm T}\widetilde{{\bf B}}_{22}^{-1}{\bf w}$, and thus ${\rm EM}_{n}^{(K)}$ is stochastically less than $\chi^{2}({d(d+1)/2})$. Therefore, we can always use $\chi^{2}({d(d+1)/2})$ as the limiting distribution. The test will be conservative, but our empirical studies show that the test still has a high power.

In the simulations, we set the number of clusters as $G=5$ and the proportions of the clusters as ${\bm{\alpha}}=(\alpha_{1},\dots,\alpha_{5})=(0.5,0.125,0.125,0.125,0.125)$. The sample size is set as $n=1000$. The dimension is set as $p=500$, $5000$ or $p=20,000$. Skmeans and COSCI are not evaluated for $p=20,000$ because they are computationally too expensive for this ultra-high dimensional setting. The number of cluster-relevant features is fixed at $s=20$. We always set the first $20$ features as the cluster-relevant and all other features as cluster-irrelevant.

More specifically, for the $i$th sample, we first randomly assign it to a cluster $g$ with the probability $\alpha_{g}$. Then, if the $j$th feature is cluster-relevant ($j=1,\dots,20$), we randomly sample $x_{ij}$ from ${\rm NB}(\mu_{gj},r_{j})$; If it is cluster-irrelevant ($j=21,\dots,p$), we randomly sample $x_{ij}$ from ${\rm NB}(\mu_{j},r_{j})$. The mean and over-dispersion parameters of the negative binomial distributions are randomly generated (see below).

We consider simulation setups of two noise levels (low or high) and three cluster signal strength levels (low, medium and high). The over-dispersion parameters $r_{j}$ represent the noise level of the data and the differences of the mean parameters $\mu_{gj}$ between clusters represent the cluster signal strength. The details of generating $r_{j}$ and $\mu_{gj}$ are given in the Supplementary material. Thus, in total, we have 18 different simulation setups (3 dimension setups $\times$ 2 noise levels $\times$ 3 signal levels). In each simulation setup, we generate 100 datasets.

To evaluate the performance of EM-test on continuous data, we also generate simulation data based on normal distributions. The simulation setups are detailed in Supplementary material. EM-test for normal models are used for these continuous data. To investigate the robustness of EM-test to model mis-specification, we further generate count data based on Poisson-truncated-normal and the binomial-Gamma distributions. EM-test for negative binomial model is used for these count data (Supplementary Section D). In the following, we only discuss the negative binomial simulations. For the continuous data simulation, we find that EM-test performs similar to other available methods under easier simulation setups and outperforms other methods under more difficult setups (Supplementary Section D). From the mis-specified count data simulations, we find that EM-test is robust to model mis-specification and outperforms other methods (Supplementary Section D).

We first evaluate the accuracy of feature screening. For different screening methods, we first rank the features by their corresponding test statistics / p-values or feature weights. Following previous researches about feature screening (Zhu et al., 2011; Li et al., 2012), let $\mathcal{S}$ be the minimum number of features needed to include all cluster-relevant features in a rank. Table 1 shows the mean and the standard deviation of $\mathcal{S}$ over the 100 replications. Table 1 does not include COSCI because it only reports a selected feature index but does not provide an order of all features. In the low dimensional cases ($p=500$), all count-data methods work well. EM-test only needs 20 or slightly more than 20 features to include all cluster-relevant features. In higher dimensions ($p=5000$ or $20,000$), the EM-test outperforms other methods, often by a large amount. For example, in the medium signal, high noise and $p=20,000$ case, the EM-test needs around 21 features to include all cluster-relevant features, while other methods need over a thousand features. SC-FS works well in lower dimensional cases, but its performance deteriorates in higher dimensional cases, especially in higher dimensional cases with lower signal to noise ratios. This is because SC-FS needs a pre-cluster label for feature screening. When $p$ is large, the pre-cluster results can be very inaccurate, leading to the inferior performance of SC-FS in these settings. The continuous methods (KS-test and Dip-test) do not perform well for these count data.

The minimum model size $\mathcal{S}$ measures the feature ranks given by different methods. However, in clustering analysis, simply having $\mathcal{S}$ close to $s$ is inadequate because we need a criterion to determine which features to retain. Therefore, we further compare the number of correctly retained cluster-relevant features (denoted as $\mathcal{R}$) and falsely retained cluster-irrelevant features (denoted as $\mathcal{F}$) by different methods. For SC-FS, COSCI and Skmeans, we use their default parameters to select the cluster-relevant features. For the EM-test, we select the features by the adjusted p-values (EM-adjust, adjusted p-value $<$ 0.01) and by the threshold $n^{0.35}$ (EM-0.35). The p-value is calculated using the $\chi^{2}(3)$-distribution, because compared with the limiting distribution in Theorem 4, the $\chi^{2}(3)$-distribution is computationally more efficient, achieves good sensitivity and false discovery rate (FDR) control (Supplementary Section D). The Benjamini-Hochberg (BH) procedure (Benjamini and Hochberg, 1995a) is used to adjust the p-values. We choose the threshold $n^{0.35}$ because the EM-test has a good balance between retaining cluster-relevant features and excluding cluster-irrelevant features at these cutoffs (Supplementary Section D). For the Chi-square goodness-of-fit test, KS-test and Dip-test, we use the BH-adjusted p-values ($<0.01$) to screen the cluster-relevant features.

Table 2 shows the numbers of correctly retained and falsely retained features by different methods. Similarly, we find that the EM-test methods (EM-adjust, EM-0.35) outperform other methods in most settings especially when $p$ is larger and different clusters are more similar to each other. In most cases, Skmeans is able to select all cluster-relevant features, but also falsely select many cluster-irrelevant features. SC-FS is conservative. In low dimensional settings ($p=500$), SC-FS could correctly select all cluster-relevant features with almost no false positives. However, in higher dimensions, SC-FS also has almost zero false positives, but its power is low. For example, in the $p=5000$, medium signal and high noise case, SC-FS only reports 2 cluster-relevant features. In the same case, EM-adjust reports all 20 features with almost zero false positives. The performance of two versions of EM-test are slightly different. EM-adjust is more conservative than EM-0.35. In the more difficult settings (with large $p$ and low signal to noise ratio), EM-adjust is still able to control the false positives, but detects less cluster-relevant features. In most cases, EM-0.35 can detect most cluster-relevant features, but also report some cluster-irrelevant features in the more difficult simulation settings. KS-test and Dip-test report many false positives and select almost all features as clustering-relevant features. COCSI is very conservative in this simulation and could not select any features.