Clustering with Statistical Error Control

TL;DR

This paper introduces a clustering method that incorporates statistical error control, enabling rigorous control over overestimation and underestimation probabilities of true clusters, with theoretical guarantees and practical applications.

Contribution

It develops estimators for the number of clusters and the clusters themselves that incorporate a tunable significance level for error control, a novel approach in clustering.

Findings

The estimators effectively control overestimation probability.

The method provides asymptotic control over underestimation.

Applications demonstrate practical utility in temperature and gene expression data.

Abstract

This paper presents a clustering approach that allows for rigorous statistical error control similar to a statistical test. We develop estimators for both the unknown number of clusters and the clusters themselves. The estimators depend on a tuning parameter alpha which is similar to the significance level of a statistical hypothesis test. By choosing alpha, one can control the probability of overestimating the true number of clusters, while the probability of underestimation is asymptotically negligible. In addition, the probability that the estimated clusters differ from the true ones is controlled. In the theoretical part of the paper, formal versions of these statements on statistical error control are derived in a standard model setting with convex clusters. A simulation study and two applications to temperature and gene expression microarray data complement the theoretical analysis.

Tables

Table 1. Table 1: Estimates of the probabilities ℙ ​ ( K ^ 0 < K 0 ) ℙ subscript ^ 𝐾 0 subscript 𝐾 0 \mathbb{P}(\widehat{K}_{0}<K_{0}) , ℙ ​ ( K ^ 0 = K 0 ) ℙ subscript ^ 𝐾 0 subscript 𝐾 0 \mathbb{P}(\widehat{K}_{0}=K_{0}) and ℙ ​ ( K ^ 0 > K 0 ) ℙ subscript ^ 𝐾 0 subscript 𝐾 0 \mathbb{P}(\widehat{K}_{0}>K_{0}) in the three simulation scenarios with the noise-to-signal ratios NSR = 1 NSR 1 \textnormal{NSR}=1 , 1.5 1.5 1.5 and 2 2 2 .

	$\text{NSR}=1$	$\text{NSR}=1.5$	$\text{NSR}=2$
	0.000	0.274	0.631
$\mathbb{P}\left({\hat{K}}_0=K_0\right)$	0.941	0.549	0.242
$\mathbb{P}\left({\hat{K}}_0>K_0\right)$	0.059	0.177	0.127

Table 2. Table 2: Estimates of the probabilities ℙ ​ ( K ^ 0 < K 0 ) ℙ subscript ^ 𝐾 0 subscript 𝐾 0 \mathbb{P}(\widehat{K}_{0}<K_{0}) , ℙ ​ ( K ^ 0 = K 0 ) ℙ subscript ^ 𝐾 0 subscript 𝐾 0 \mathbb{P}(\widehat{K}_{0}=K_{0}) and ℙ ​ ( K ^ 0 > K 0 ) ℙ subscript ^ 𝐾 0 subscript 𝐾 0 \mathbb{P}(\widehat{K}_{0}>K_{0}) in the simulation scenario with NSR = 1.5 NSR 1.5 \textnormal{NSR}=1.5 and five different sample sizes ( n , p ) 𝑛 𝑝 (n,p) .

$(n,p)$	$(1000,30)$	$(1500,40)$	$(2000,50)$	$(2500,50)$	$(3000,60)$
	0.274	0.037	0.002	0.000	0.000
$\mathbb{P}\left({\hat{K}}_0=K_0\right)$	0.549	0.795	0.893	0.918	0.956
$\mathbb{P}\left({\hat{K}}_0>K_0\right)$	0.177	0.168	0.105	0.082	0.044

Table 3. Table 3: Results of the comparison study. The entries of the table give the numbers of simulations (out of a total of 100 100 100 ) for which a certain estimate of K 0 subscript 𝐾 0 K_{0} is obtained. The first line in part (a) of the table, for example, has to be read as follows: The CluStErr estimate K ^ 0 subscript ^ 𝐾 0 \widehat{K}_{0} is equal to the true K 0 = 10 subscript 𝐾 0 10 K_{0}=10 in 95 95 95 out of 100 100 100 simulations, and it is equal to K = 11 , 12 , 13 𝐾 11 12 13 K=11,12,13 in 2 , 1 , 2 2 1 2 2,1,2 simulations, respectively.

		Estimated number of clusters
NSR	Method	1	2	3	4	5	6	7	8	9	10	11	12	13	14	$\ge 15$
1	CluStErr	0	0	0	0	0	0	0	0	0	95	2	1	2	0	0
	Gap	3	5	0	0	1	3	11	25	36	16	0	0	0	0	0
	Silhouette	0	0	0	0	0	0	0	5	19	43	20	11	1	0	1
	Hartigan	0	0	0	0	0	0	0	1	4	40	31	10	5	6	3
1.5	CluStErr	0	0	0	0	0	0	0	1	22	59	10	3	3	1	1
	Gap	13	26	0	0	0	0	2	13	17	29	0	0	0	0	0
	Silhouette	0	0	0	0	0	0	0	2	12	54	23	8	1	0	0
	Hartigan	0	0	0	0	0	0	0	1	3	61	18	9	2	4	2
2	CluStErr	0	0	0	0	1	2	10	13	32	31	7	2	1	0	1
	Gap	22	26	2	0	0	0	1	6	9	34	0	0	0	0	0
	Silhouette	0	0	0	0	0	0	0	0	6	66	26	2	0	0	0
	Hartigan	0	0	0	0	0	0	0	0	10	68	18	2	2	0	0

Table 4. (a) balanced scenario

		Estimated number of clusters
NSR	Method	1	2	3	4	5	6	7	8	9	10	11	12	13	14	$\ge 15$
1	CluStErr	0	0	0	0	0	0	0	0	0	95	2	1	2	0	0
	Gap	3	5	0	0	1	3	11	25	36	16	0	0	0	0	0
	Silhouette	0	0	0	0	0	0	0	5	19	43	20	11	1	0	1
	Hartigan	0	0	0	0	0	0	0	1	4	40	31	10	5	6	3
1.5	CluStErr	0	0	0	0	0	0	0	1	22	59	10	3	3	1	1
	Gap	13	26	0	0	0	0	2	13	17	29	0	0	0	0	0
	Silhouette	0	0	0	0	0	0	0	2	12	54	23	8	1	0	0
	Hartigan	0	0	0	0	0	0	0	1	3	61	18	9	2	4	2
2	CluStErr	0	0	0	0	1	2	10	13	32	31	7	2	1	0	1
	Gap	22	26	2	0	0	0	1	6	9	34	0	0	0	0	0
	Silhouette	0	0	0	0	0	0	0	0	6	66	26	2	0	0	0
	Hartigan	0	0	0	0	0	0	0	0	10	68	18	2	2	0	0

Table 5. (b) unbalanced scenario

		Estimated number of clusters
NSR	Method	1	2	3	4	5	6	7	8	9	10	11	12	13	14	$\ge 15$
1	CluStErr	0	0	0	0	0	0	0	0	0	94	4	2	0	0	0
	Gap	0	1	4	9	17	11	21	20	16	1	0	0	0	0	0
	Silhouette	0	0	0	0	0	0	10	26	37	19	6	2	0	0	0
	Hartigan	0	0	1	2	7	8	21	21	12	11	7	2	2	1	5
1.5	CluStErr	0	0	0	0	0	0	0	3	15	55	20	5	2	0	0
	Gap	2	4	4	3	8	15	21	21	22	0	0	0	0	0	0
	Silhouette	0	0	0	0	0	1	7	20	47	24	1	0	0	0	0
	Hartigan	0	0	0	1	14	20	12	14	20	8	7	4	0	0	0
2	CluStErr	0	0	0	0	0	0	3	17	25	31	15	3	6	0	0
	Gap	14	4	5	1	6	6	19	26	19	0	0	0	0	0	0
	Silhouette	0	0	0	0	0	0	5	28	48	19	0	0	0	0	0
	Hartigan	0	0	1	4	15	15	15	18	17	7	5	0	0	2	1

Table 6. Table 4: Simulation results for the blocked maximum statistic.

	Estimated number of clusters
NSR	1	2	3	4	5	6	7	8	9	10	11	12	13	14
1	0	0	0	0	0	0	0	0	0	98	1	1	0	0
1.5	0	0	0	0	0	0	0	0	0	79	12	5	2	2
2	0	0	0	0	0	0	0	0	0	73	21	2	3	1

Equations

$$\displaystyle\mathbb{P}\big{(}\widehat{K}_{0}>K_{0}\big{)}$$

$$\displaystyle\mathbb{P}\big{(}\widehat{K}_{0}<K_{0}\big{)}$$

$${\mathbf{Y}}_i={\mathbf{\mu }}_i+{\mathbf{e}}_i,$$

$${\mathbf{\mu }}_i={\mathbf{m}}_k\text{ for all }i\in G_k,$$

$${\mathbf{Y}}_i={\mathbf{\mu }}_i+{\mathbf{\alpha }}_i+{\mathbf{\epsilon }}_i,$$

$$Y_{ij}={\mu }_i(t_j)+{\alpha }_i+{\epsilon }_{ij}(1\le j\le p),$$

$$Y_{ij}={\mu }_{ij}+{\alpha }_i+{\epsilon }_{ij}(1\le j\le p).$$

$$\mathbb{P}\Big{(}\big{\{}\widehat{G}_{k}^{[K_{0}]}:1\leq k\leq K_{0}\big{\}}=\big{\{}G_{k}:1\leq k\leq K_{0}\big{\}}\Big{)}\rightarrow 1.$$

$$\widehat{K}_{0}=\min\big{\{}K=1,2,\ldots\,\big{|}\,\widehat{\mathcal{H}}^{[K]}\leq q(\alpha)\big{\}}.$$

$$\widehat{K}_{0}=\min\big{\{}K=1,2,\ldots\,\big{|}\,\widehat{p}^{[K]}>\alpha\big{\}},$$

$$\mathbb{P}\Big{(}\widehat{\mathcal{H}}^{[K]}\leq q(\alpha)\Big{)}=\begin{cases}o(1)&\text{for }K<K_{0}\\ (1-\alpha)+o(1)&\text{for }K=K_{0}.\end{cases}$$

$${\pi }_{>}(\alpha )$$

$${\pi }_{<}(\alpha )$$

$$\mathbb{P}\big{(}\widehat{K}_{0}\neq K_{0}\big{)}=\alpha+o(1),$$

$$\mathbb{P}\Big{(}\big{\{}\widehat{G}_{k}:1\leq k\leq\widehat{K}_{0}\big{\}}\neq\big{\{}G_{k}:1\leq k\leq K_{0}\big{\}}\Big{)}=\alpha+o(1),$$

$${\hat{Y}}_{ij}$$

$$={\mu }_{ij}+{\epsilon }_{ij}+O_p(p^{-1/2})=Y_{ij}^{\ast }+O_p(p^{-1/2}).$$

$$\widehat{\Delta}_{i}^{[K]}=\frac{1}{\sqrt{p}}\sum\limits_{j=1}^{p}\Big{\{}\Big{(}\frac{\widehat{\varepsilon}_{ij}^{[K]}}{\widehat{\sigma}}\Big{)}^{2}-1\Big{\}}\Big{/}\widehat{\kappa}$$

$$\widehat{\Delta}_{i}^{[K]}=\frac{1}{\sqrt{p}}\sum\limits_{j=1}^{p}\Big{\{}\frac{(\varepsilon_{ij}+d_{ij})^{2}}{\sigma^{2}}-1\Big{\}}\Big{/}\kappa$$

$$\widehat{\Delta}_{i}^{[K_{0}]}=\frac{1}{\sqrt{p}}\sum\limits_{j=1}^{p}\Big{\{}\frac{\varepsilon_{ij}^{2}}{\sigma^{2}}-1\Big{\}}\Big{/}\kappa.$$

$${\hat{\Delta }}_i^{[K_0]}\sim \frac{{\chi }_p^2-p}{\sqrt{2p}}$$

$${\hat{\Delta }}_i^{[K]}$$

$$=\frac{1}{\sqrt{p}}\sum _{j=1}^p\frac{d_{ij}^2}{{\sigma }^2\kappa }+O_p(1)$$

$${\hat{\Delta }}_i^{[K]}\ge c\sqrt{p}\text{for some }c>0\text{ with prob. tending to }1,$$

$${\hat{\mathcal{H}}}^{[K]}=\underset{1\le i\le n}{\max }{\hat{\Delta }}_i^{[K]}$$

$$\mathbb{P}\Big{(}\widehat{\mathcal{H}}^{[K_{0}]}\leq q(\alpha)\Big{)}\approx(1-\alpha).$$

$$\mathbb{P}\Big{(}\widehat{\mathcal{H}}^{[K]}\leq q(\alpha)\Big{)}=o(1)$$

$$\widehat{\rho}(i,i^{\prime})=\frac{1}{p}\sum\limits_{j=1}^{p}\big{(}\widehat{Y}_{ij}-\widehat{Y}_{i^{\prime}j}\big{)}^{2}.$$

$${\hat{\sigma }}_{\text{Lip}}^2=\frac{1}{n(p-1)}\sum _{i=1}^n\sum _{j=2}^p\frac{(Y_{ij}-Y_{i,j-1}{)}^2}{2}.$$

$${\hat{\vartheta }}_{\text{Lip}}=\frac{1}{n(p-1)}\sum _{i=1}^n\sum _{j=2}^p\frac{(Y_{ij}-Y_{i,j-1}{)}^4}{2}-3({\hat{\sigma }}_{\text{Lip}}^2{)}^2,$$

$$\widehat{\kappa}_{\text{Lip}}=\Big{(}\frac{\widehat{\vartheta}_{\text{Lip}}}{(\widehat{\sigma}_{\text{Lip}}^{2})^{2}}-1\Big{)}^{1/2}.$$

$${\hat{\sigma }}_{\text{pc}}^2$$

Full text

Clustering with

Statistical Error Control

Michael Vogt111Corresponding author. Address: Department of Economics and Hausdorff Center for Mathematics, University of Bonn, 53113 Bonn, Germany. Email: michael.vogt@uni-bonn.de.

University of Bonn

Matthias Schmid222Address: Department of Medical Biometry, Informatics and Epidemiology, University of Bonn, 53105 Bonn, Germany. matthias.schmid@imbie.uni-bonn.de.

University of Bonn

This paper presents a clustering approach that allows for rigorous statistical error control similar to a statistical test. We develop estimators for both the unknown number of clusters and the clusters themselves. The estimators depend on a tuning parameter $\alpha$ which is similar to the significance level of a statistical hypothesis test. By choosing $\alpha$, one can control the probability of overestimating the true number of clusters, while the probability of underestimation is asymptotically negligible. In addition, the probability that the estimated clusters differ from the true ones is controlled. In the theoretical part of the paper, formal versions of these statements on statistical error control are derived in a standard model setting with convex clusters. A simulation study and two applications to temperature and gene expression microarray data complement the theoretical analysis.

Key words: Cluster analysis; number of clusters; multiple statistical testing; statistical error control; $k$-means clustering.

AMS 2010 subject classifications: 62H30; 62H15; 62E20.

1 Introduction

In a wide range of applications, the aim is to cluster a large number of subjects into a small number of groups. Prominent examples are the clustering of genes in microarray analysis (Jiang et al., 2004), the clustering of temperature curves using data recorded on a spatial grid (Fovell and Fovell, 1993; DeGaetano, 2001), and the clustering of consumer profiles on the basis of survey data (Wedel and Kamakura, 2000).

A major challenge in cluster analysis is to estimate the unknown number of groups $K_{0}$ from a sample of data. A common approach is to compute a criterion function which measures the quality of the clustering for different cluster numbers $K$. An estimator of $K_{0}$ is then obtained by optimizing the criterion function over $K$. Prominent examples of this approach are the Hartigan index (Hartigan, 1975), the silhouette statistic (Rousseeuw, 1987) and the gap statistic (Tibshirani et al., 2001).

Another common way to estimate $K_{0}$ is based on statistical test theory. Roughly speaking, one can distinguish between two types of test-based procedures: The first type relies on a statistical test which either checks whether some clusters can be merged or whether a cluster can be subdivided. Given a set of clusters, the test is repeatedly applied until no clusters can be merged or split any more. The number of remaining clusters serves as an estimator of $K_{0}$. Classical examples of methods that proceed in this way are discussed in Gordon (1999, Chapter 3.5) who terms them “local methods”. Obviously, these methods involve a multiple testing problem. However, the employed critical values do not properly control for the fact that multiple tests are performed. The significance level $\alpha$ used to carry out the tests thus cannot be interpreted strictly. Put differently, the procedures do not allow for rigorous statistical error control.

Test-based approaches of the second type proceed by sequentially testing a model with $K$ clusters against one with $K+1$ clusters. The smallest number $K$ for which the test does not reject serves as an estimator of $K_{0}$. Most work in this direction has been done in the framework of Gaussian mixture models; see McLachlan and Rathnayake (2014) for an overview. However, deriving a general theory for testing a mixture with $K$ components against one with $K^{\prime}>K$ components has turned out to be a very challenging problem; see Ghosh and Sen (1985) and Hartigan (1985) for a description of the main technical issues involved. Many results are therefore restricted to the special case of testing a homogeneous model against a mixture with $K=2$ clusters; see Liu and Shao (2004) and Li et al. (2009) among many others. More general test procedures often lack a complete theoretical foundation or are based on very restrictive conditions.

Only recently, there have been some advances in developing a general theory for testing $K$ against $K^{\prime}>K$ clusters under reasonably weak conditions. In a mixture model setup, Li and Chen (2010) and Chen et al. (2012) have constructed a new expectation-maximization (EM) procedure to approach this testing problem. Outside the mixture model context, Maitra et al. (2012) have developed a bootstrap procedure to test a model with $K$ groups against one with $K^{\prime}>K$ groups. These papers derive the theoretical properties of the proposed tests under the null hypothesis of $K$ clusters, where $K$ is a pre-specified fixed number. However, they do not formally investigate the properties of a procedure which estimates $K_{0}$ by sequentially applying the tests. In particular, they do not analyze whether such a sequential procedure may allow for a rigorous interpretation of the significance level $\alpha$ that is used to carry out the tests.

The main contribution of this paper is to construct an estimator $\widehat{K}_{0}$ of $K_{0}$ which allows for rigorous statistical error control in the following sense: For any pre-specified significance level $\alpha\in(0,1)$, the proposed estimator $\widehat{K}_{0}=\widehat{K}_{0}(\alpha)$ has the property that

[TABLE]

According to this, the probability of overestimating $K_{0}$ is controlled by the level $\alpha$, while the probability of underestimating $K_{0}$ is asymptotically negligible. By picking $\alpha$, we can thus control the probability of choosing too many clusters, while, on the other hand, we can ignore the probability of choosing too few clusters (at least asymptotically).

We show how to construct an estimator $\widehat{K}_{0}$ with the properties (1.1) and (1.2) in a standard model setting with convex clusters which is introduced in Section 2. Our estimation approach is developed in Section 3. As we will see, the proposed procedure does not only provide us with an estimator of $K_{0}$. It also yields estimators of the groups themselves which allow for statistical error control similarly to $\widehat{K}_{0}$. Our approach is based on the following general strategy:

1. (i)

Construct a statistical test which, for any given number $K$, checks the null hypothesis that there are $K$ clusters in the data. 2. (ii)

Starting with $K=1$, sequentially apply this test until it does not reject the null hypothesis of $K$ clusters any more. 3. (iii)

Define the estimator $\widehat{K}_{0}$ of $K_{0}$ as the smallest number $K$ for which the test does not reject the null.

This strategy is discussed in detail in Section 3.1. It is generic in the sense that it can be employed with different test statistics. For our theoretical analysis, we apply it with a specific statistic which is introduced in Section 3.2. For this specific choice, we derive the statements (1.1) and (1.2) on statistical error control under suitable regularity conditions. Some alternative choices of the test statistic are discussed in Section 6. In the following, we refer to our estimation procedure as CluStErr (“Clustering with Statistical Error Control”).

The theoretical properties of our estimators, in particular the statements (1.1) and (1.2), are derived in Section 4. As we will see there, our theory is valid under quite general conditions. First of all, as opposed to many other studies from the clustering literature including those from a Gaussian mixture context, we do not restrict the random variables in our model to be Gaussian. For our theory to work, we merely require them to satisfy a set of moment conditions. Secondly, our approach is essentially free of tuning parameters, the only choice parameter being the significance level $\alpha$. Thirdly, to apply our method, we of course need to compute critical values for the underlying test. However, as opposed to other test-based methods, we do not have to estimate or bootstrap the critical values by a complicated procedure. They can rather be easily computed analytically. This makes our method particularly simple to implement in practice.

We complement the theoretical analysis of the paper by a simulation study and two applications on temperature and microarray data in Section 5. The R code to reproduce the numerical examples is contained in the add-on package CluStErr (Lasota et al., 2017), which implements the CluStErr method and which is part of the supplemental materials of the paper.

2 Model

Suppose we measure $p$ features on $n$ different subjects. In particular, for each subject $i\in\{1,\ldots,n\}$, we observe the vector $\boldsymbol{Y}_{i}=(Y_{i1}\,\ldots,Y_{ip})^{\top}$, where $Y_{ij}$ denotes the measurement of the $j$-th feature for the $i$-th subject. Our data sample thus has the form $\{\boldsymbol{Y}_{i}:1\leq i\leq n\}$. Both the number of subjects $n$ and the number of features $p$ are assumed to tend to infinity, with $n$ diverging much faster than $p$. This reflects the fact that $n$ is much larger than $p$ in the applications we have in mind. When clustering the genes in a typical microarray data set, for instance, the number of genes $n$ is usually a few thousands, whereas the number of tissue samples $p$ is not more than a few tenths. The exact technical conditions on the sizes of $n$ and $p$ are laid out in Section 4.1.

The data vectors $\boldsymbol{Y}_{i}$ of the various subjects $i=1,\ldots,n$ are supposed to satisfy the model

[TABLE]

where $\boldsymbol{\mu}_{i}=(\mu_{i1},\ldots,\mu_{ip})^{\top}$ is a deterministic signal vector and $\boldsymbol{e}_{i}=(e_{i1},\ldots,e_{ip})^{\top}$ is the noise vector. The subjects in our sample are assumed to belong to $K_{0}$ different classes. More specifically, the set of subjects $\{1,\ldots,n\}$ can be partitioned into $K_{0}$ groups $G_{1},\ldots,G_{K_{0}}$ such that for each $k=1,\ldots,K_{0}$,

[TABLE]

where $\boldsymbol{m}_{k}\in\mathbb{R}^{p}$ are vectors with $\boldsymbol{m}_{k}\neq\boldsymbol{m}_{k^{\prime}}$ for $k\neq k^{\prime}$. Hence, the members of each group $G_{k}$ all have the same signal vector $\boldsymbol{m}_{k}$.

Equations (2.1) and (2.2) specify a model with convex spherical clusters which underlies the $k$-means and many other Euclidean distance-based clustering algorithms. This framework has been employed extensively in the literature and is useful in a wide range of applications, which is also illustrated by the examples in Section 5. It is thus a suitable baseline model for developing our ideas on clustering with statistical error control. We now discuss the two model equations (2.1) and (2.2) in detail.

Details on equation (2.1). The noise vector $\boldsymbol{e}_{i}=(e_{i1},\ldots,e_{ip})^{\top}$ is assumed to consist of entries $e_{ij}$ with the additive component structure $e_{ij}=\alpha_{i}+\varepsilon_{ij}$. Equation (2.1) for the $i$-th subject thus writes as

[TABLE]

where $\boldsymbol{\alpha}_{i}=(\alpha_{i},\ldots,\alpha_{i})^{\top}$ and $\boldsymbol{\varepsilon}_{i}=(\varepsilon_{i1},\ldots,\varepsilon_{ip})^{\top}$. Here, $\alpha_{i}$ is a subject-specific random intercept term. Moreover, the terms $\varepsilon_{ij}$ are standard idiosyncratic noise variables with $\mathbb{E}[\varepsilon_{ij}]=0$. We assume the error terms $\varepsilon_{ij}$ to be i.i.d. both across $i$ and $j$. The random intercepts $\alpha_{i}$, in contrast, are allowed to be dependent across subjects $i$ in an arbitrary way.

In general, the components of (2.3) may depend on the sample size $p$. The exact formulation of the model equation (2.3) for the $i$-th subject thus reads $\boldsymbol{Y}_{i,p}=\boldsymbol{\mu}_{i,p}+\boldsymbol{\alpha}_{i,p}+\boldsymbol{\varepsilon}_{i,p}$, where $\boldsymbol{Y}_{i,p}=(Y_{i1,p},\ldots,Y_{ip,p})^{\top}$, $\boldsymbol{\mu}_{i,p}=(\mu_{i1,p},\ldots,\mu_{ip,p})^{\top}$, $\boldsymbol{\alpha}_{i,p}=(\alpha_{i,p},\ldots,\alpha_{i,p})^{\top}$ and $\boldsymbol{\varepsilon}_{i,p}=(\varepsilon_{i1,p},\ldots,\varepsilon_{ip,p})^{\top}$. However, to keep the notation simple, we suppress this dependence on $p$ and write the model for the $i$-th subject as (2.3).

If we drop the random intercept $\boldsymbol{\alpha}_{i}$ from (2.3), the signal vector $\boldsymbol{\mu}_{i}$ is equal to the mean $\mathbb{E}[\boldsymbol{Y}_{i}]$. In the general equation (2.3) in contrast, $\boldsymbol{\mu}_{i}$ is only identified up to an additive constant. To identify $\boldsymbol{\mu}_{i}$ in (2.3), we impose the normalization constraint $p^{-1}\sum\nolimits_{j=1}^{p}\mu_{ij}=0$ for each $i$. We thus normalize the entries of $\boldsymbol{\mu}_{i}$ to be zero on average for each $i$. Under the technical conditions specified in Section 4.1, the constraint $p^{-1}\sum\nolimits_{j=1}^{p}\mu_{ij}=0$ implies that $\alpha_{i}=\lim_{p\rightarrow\infty}p^{-1}\sum\nolimits_{j=1}^{p}Y_{ij}$ almost surely, which in turn identifies the signal vector $\boldsymbol{\mu}_{i}$.

Details on equation (2.2). This equation specifies the group structure in our model. We assume the number of groups $K_{0}$ to be fixed, implying that the groups $G_{k}=G_{k,n}$ depend on the sample size $n$. Keeping the number of classes $K_{0}$ fixed while letting the size of the classes $G_{k,n}$ grow is a reasonable assumption: It reflects the fact that in most applications, we expect the number of groups $K_{0}$ to be very small as compared to the total number of subjects $n$. To keep the notation simple, we suppress the dependence of the classes $G_{k,n}$ on the sample size $n$ and denote them by $G_{k}$ throughout the paper.

In the remainder of this section, we discuss two special cases of model (2.1)–(2.2) which are relevant for our applications in Section 5.

A model for the clustering of time series data. Suppose we observe time series $\boldsymbol{Y}_{i}=(Y_{i1},\ldots,Y_{ip})^{\top}$ of length $p$ for $n$ different subjects $i$. The time series $\boldsymbol{Y}_{i}$ of the $i$-th subject is assumed to follow the time trend model

[TABLE]

where $\mu_{i}(\cdot)$ is an unknown nonparametric trend function and $t_{1}<\ldots<t_{p}$ are the observed time points. The deterministic design points $t_{j}$ are supposed to be the same across subjects $i$ and are normalized to lie in the unit interval. An important example is the equidistant design $t_{j}=j/p$. However, it is also possible to allow for non-equidistant designs. To identify the trend function $\mu_{i}(\cdot)$ in (2.4), we suppose that $\int_{0}^{1}\mu_{i}(w)dw=0$ for each $i$, which is a slight modification of the identification constraint stipulated in (2.3). Analogous to our general model, we impose a group structure on the observed time series: There are $K_{0}$ groups of time series $G_{1},\ldots,G_{K_{0}}$ such that $\mu_{i}(\cdot)=m_{k}(\cdot)$ for all $i\in G_{k}$. Hence, the members of each class $G_{k}$ all have the same time trend function $m_{k}(\cdot)$.

A model for the clustering of genes in microarray experiments. In a microarray experiment, the expression levels of $n$ different genes are often measured in $p$ different tissue samples (obtained, e.g., from $p$ different patients). For each gene $i$, we observe the vector $\boldsymbol{Y}_{i}=(Y_{i1},\ldots,Y_{ip})^{\top}$, where $Y_{ij}$ is the measured expression level of gene $i$ for tissue sample $j$. The vector $\boldsymbol{Y}_{i}$ of gene $i$ is supposed to satisfy the model equation (2.3), which componentwise reads as

[TABLE]

Here, $\mu_{ij}$ can be regarded as the true expression level of gene $i$ for tissue $j$, whereas $Y_{ij}$ is the measured expression level corrupted by the noise term $\alpha_{i}+\varepsilon_{ij}$.

Most microarray experiments involve different types of tissues, for example tumor “cases” versus healthy “controls”, or different tumor (sub)types. We therefore suppose that there are $T$ different types of tissues in our sample and order them according to their type (which is known by experimental design). More specifically, the tissues $j$ of type $t$ are labelled by $j_{t-1}\leq j<j_{t}$, where $1=j_{0}<j_{1}<\ldots<j_{T-1}<j_{T}=p+1$. If the patients from which tissues are obtained constitute samples of sufficiently homogeneous populations, it is natural to assume that the true expression level $\mu_{ij}$ of gene $i$ is the same for tissues $j$ of the same type, i.e., $\mu_{ij}=\mu_{ij^{\prime}}$ for $j_{t-1}\leq j,j^{\prime}<j_{t}$. The signal vector $\boldsymbol{\mu}_{i}$ thus has a piecewise constant structure for each $i$; see Figures 4 and 6 in Section 5 for an illustration.

As in our general model, we suppose that there are $K_{0}$ groups of genes $G_{1},\ldots,G_{K_{0}}$ such that $\boldsymbol{\mu}_{i}=\boldsymbol{m}_{k}$ for all $i\in G_{k}$ and some vector $\boldsymbol{m}_{k}$. The genes of each class $G_{k}$ thus have the same (co-)expression profile $\boldsymbol{m}_{k}$.

3 Estimation Method

We now present our approach to estimate the unknown groups $G_{1},\ldots,G_{K_{0}}$ and their unknown number $K_{0}$ in model (2.1)–(2.2). Section 3.1 gives an overview of the general method, while Sections 3.2–3.4 fill in the details.

3.1 The general method

To construct our method, we proceed in two steps: In the first step, we specify an algorithm that clusters the set of subjects $\{1,\ldots,n\}$ into $K$ groups for any given number $K$ (which may or may not coincide with the true number of classes $K_{0}$). Let $\{\widehat{G}_{k}^{[K]}:1\leq k\leq K\}$ be the $K$ clusters produced by the algorithm when the number of clusters is $K$. For $K=1$, we trivially set $\widehat{G}_{1}^{[1]}=\{1,\ldots,n\}$. For our theory to work, we require the clustering algorithm to consistently estimate the class structure $\{G_{k}:1\leq k\leq K_{0}\}$ when $K=K_{0}$. More specifically, we require the estimators $\{\widehat{G}_{k}^{[K_{0}]}:1\leq k\leq K_{0}\}$ to have the property that

[TABLE]

This is a quite weak restriction which is satisfied by a wide range of clustering algorithms under our regularity conditions. As shown in Section 3.3, it is for example satisfied by a $k$-means type algorithm. Moreover, it can be shown to hold for a number of hierarchical clustering algorithms, in particular for agglomerative algorithms with single, average and complete linkage. Our estimation method can be based on any clustering algorithm that has the consistency property (3.1).

In the second step, we construct a test for each $K$ which checks whether the data can be well described by the $K$ clusters $\widehat{G}_{1}^{[K]},\ldots,\widehat{G}_{K}^{[K]}$. We thereby test whether the number of clusters is equal to $K$. More formally, we use the $K$-cluster partition $\{\widehat{G}_{k}^{[K]}:1\leq k\leq K\}$ to construct a statistic $\widehat{\mathcal{H}}^{[K]}$ that allows us to test the hypothesis $H_{0}:K=K_{0}$ versus $H_{1}:K<K_{0}$. For any given number of clusters $K$, our test is defined as $T_{\alpha}^{[K]}=\boldsymbol{1}(\widehat{\mathcal{H}}^{[K]}>q(\alpha))$, where $q(\alpha)$ is the $(1-\alpha)$-quantile of a known distribution which will be specified later on. We reject $H_{0}$ at the level $\alpha$ if $T_{\alpha}^{[K]}=1$, i.e., if $\widehat{\mathcal{H}}^{[K]}>q(\alpha)$. A detailed construction of the statistic $\widehat{\mathcal{H}}^{[K]}$ along with a precise definition of the quantile $q(\alpha)$ is given in Section 3.2.

To estimate the classes $G_{1},\ldots,G_{K_{0}}$ and their number $K_{0}$, we proceed as follows: For each $K=1,2,\ldots$, we check whether $\widehat{\mathcal{H}}^{[K]}\leq q(\alpha)$ and stop as soon as this criterion is satisfied. Put differently, we carry out our test for each $K=1,2,\ldots$ until it does not reject $H_{0}$ any more. Our estimator of $K_{0}$ is defined as the smallest number $K$ for which $\widehat{\mathcal{H}}^{[K]}\leq q(\alpha)$, that is, for which the test does not reject $H_{0}$. Formally speaking, we define

[TABLE]

Moreover, we estimate the class structure $\{G_{k}:1\leq k\leq K_{0}\}$ by the partition $\{\widehat{G}_{k}:1\leq k\leq\widehat{K}_{0}\}$, where we set $\widehat{G}_{k}=\widehat{G}_{k}^{[\widehat{K}_{0}]}$. The definition (3.2) can equivalently be written as

[TABLE]

where $\widehat{p}^{[K]}$ is the $p$-value corresponding to the statistic $\widehat{\mathcal{H}}^{[K]}$. The heuristic idea behind (3.3) is as follows: Starting with $K=1$, we successively test whether the data can be well described by a model with $K$ clusters, in particular by the partition $\{\widehat{G}_{k}^{[K]}:1\leq k\leq K\}$. For each $K$, we compute the $p$-value $\widehat{p}^{[K]}$ which expresses our confidence in a model with $K$ clusters. We stop as soon as $\widehat{p}^{[K]}>\alpha$, that is, as soon as we have enough statistical confidence in a model with $K$ groups.

As shown in Section 4, under appropriate regularity conditions, our statistic $\widehat{\mathcal{H}}^{[K]}$ has the property that

[TABLE]

Put differently, $\mathbb{P}(T_{\alpha}^{[K]}=0)\rightarrow 1-\alpha$ for $K=K_{0}$ and $\mathbb{P}(T_{\alpha}^{[K]}=1)\rightarrow 1$ for $K<K_{0}$. Hence, our test is asymptotically of level $\alpha$. Moreover, it detects the alternative $H_{1}:K<K_{0}$ with probability tending to $1$, that is, its power against $H_{1}$ is asymptotically equal to $1$. From (3.4), it follows that

[TABLE]

Hence, the probability of overestimating $K_{0}$ is asymptotically bounded by $\alpha$, while the probability of underestimating $K_{0}$ is asymptotically negligible. By picking $\alpha$, we can thus control the probability of choosing too many clusters similarly to the type-I-error probability of a test. Moreover, we can asymptotically ignore the probability of choosing too few clusters similarly to the type-II-error probability of a test. In finite samples, there is of course a trade-off between the probabilities of under- and overestimating $K_{0}$: By decreasing the significance level $\alpha$, we can reduce the probability of overestimating $K_{0}$, since $\alpha^{\prime}+o(1)=\pi_{>}(\alpha^{\prime})\leq\pi_{>}(\alpha)=\alpha+o(1)$ for $\alpha^{\prime}<\alpha$. However, we pay for this by increasing the probability of underestimating $K_{0}$, since $\pi_{<}(\alpha^{\prime})\geq\pi_{<}(\alpha)$ for $\alpha^{\prime}<\alpha$. This can also be regarded as a trade-off between the size and the power of the test on which $\widehat{K}_{0}$ is based. Taken together, the two statements (3.5) and (3.6) yield that

[TABLE]

i.e., the probability that the estimated number of classes $\widehat{K}_{0}$ differs from the true number of classes $K_{0}$ is asymptotically equal to $\alpha$. With the help of (3.7) and the consistency property (3.1) of the estimated clusters, we can further show that

[TABLE]

i.e., the probability of making a classification error is asymptotically equal to $\alpha$ as well. The statements (3.5)–(3.8) give a mathematically precise description of the statistical error control that can be performed by our method.

3.2 Construction of the statistic $\boldsymbol{\widehat{\mathcal{H}}^{[K]}}$

To construct the statistic $\widehat{\mathcal{H}}^{[K]}$, we use the following notation:

1. (i)

Let $Y_{ij}^{*}=Y_{ij}-\alpha_{i}$ be the observations adjusted for the random intercepts $\alpha_{i}$ and set $\widehat{Y}_{ij}=Y_{ij}-\overline{Y}_{i}$ with $\overline{Y}_{i}=p^{-1}\sum\nolimits_{j=1}^{p}Y_{ij}$. The variables $\widehat{Y}_{ij}$ serve as approximations of $Y_{ij}^{*}$, since under standard regularity conditions

[TABLE] 2. (ii)

For any set $S\subseteq\{1,\ldots,n\}$, let $m_{j,S}=(\#S)^{-1}\sum\nolimits_{i\in S}\mu_{ij}$ be the average of the signals $\mu_{ij}$ with $i\in S$ and estimate it by $\widehat{m}_{j,S}=(\#S)^{-1}\sum\nolimits_{i\in S}\widehat{Y}_{ij}$. We use the notation $m_{j,k}^{[K]}=m_{j,\widehat{G}_{k}^{[K]}}$ and $\widehat{m}_{j,k}^{[K]}=\widehat{m}_{j,\widehat{G}_{k}^{[K]}}$ to denote the average of the signals in the cluster $\widehat{G}_{k}^{[K]}$ and its estimator, respectively. 3. (iii)

For any cluster $\widehat{G}_{k}^{[K]}$, we define cluster-specific residuals by setting $\widehat{\varepsilon}_{ij}^{[K]}=\widehat{Y}_{ij}-\widehat{m}_{j,k}^{[K]}$ for $i\in\widehat{G}_{k}^{[K]}$ and $1\leq j\leq p$. 4. (iv)

Let $\widehat{\sigma}^{2}$ be an estimator of the error variance $\sigma^{2}=\mathbb{E}[\varepsilon_{ij}^{2}]$. Moreover, let $\widehat{\kappa}$ be an estimator of the parameter $\kappa=(\mathbb{E}[\{(\varepsilon_{ij}/\sigma)^{2}-1\}^{2}])^{1/2}$, which serves as a normalization constant later on. See Section 3.4 for a detailed construction of the estimators $\widehat{\sigma}^{2}$ and $\widehat{\kappa}$.

With this notation at hand, we define the statistic

[TABLE]

for each subject $i$. This is essentially a scaled version of the residual sum of squares for the $i$-th subject when the number of clusters is $K$. Intuitively, $\widehat{\Delta}_{i}^{[K]}$ measures how well the data of the $i$-th subject are described when the sample of subjects is partitioned into the $K$ clusters $\widehat{G}_{1}^{[K]},\ldots,\widehat{G}_{K}^{[K]}$. The individual statistics $\widehat{\Delta}_{i}^{[K]}$ are the building blocks of the overall statistic $\widehat{\mathcal{H}}^{[K]}$.

Before we move on with the construction of $\widehat{\mathcal{H}}^{[K]}$, we have a closer look at the stochastic behaviour of the statistics $\widehat{\Delta}_{i}^{[K]}$. To do so, we consider the following stylized situation: We assume that the variables $\varepsilon_{ij}$ are i.i.d. normally distributed with mean [math] and variance $\sigma^{2}$. Moreover, we neglect the estimation error in the expressions $\widehat{Y}_{ij}$, $\widehat{m}_{j,k}^{[K]}$, $\widehat{\sigma}^{2}$ and $\widehat{\kappa}$. In this situation,

[TABLE]

for any $i\in S=\widehat{G}_{k}^{[K]}$, where $d_{ij}=\mu_{ij}-(\#S)^{-1}\sum\nolimits_{i^{\prime}\in S}\mu_{i^{\prime}j}$ is the difference between the signal $\mu_{ij}$ of the $i$-th subject and the average signal in the cluster $S$. We now give a heuristic discussion of the behaviour of $\widehat{\Delta}_{i}^{[K]}$ in the following two cases:

• $K=K_{0}$:

By condition (3.1), $\widehat{G}_{k}^{[K_{0}]}$ consistently estimates $G_{k}$. Neglecting the estimation error in $\widehat{G}_{k}^{[K_{0}]}$, we obtain that

[TABLE]

Since $\varepsilon_{ij}/\sigma$ is standard normal, $\kappa=\sqrt{2}$ and thus

[TABLE]

for each $i$. Hence, the individual statistics $\widehat{\Delta}_{i}^{[K_{0}]}$ all have a rescaled $\chi^{2}$-distribution.

• $K<K_{0}$:

If we pick $K$ smaller than the true number of classes $K_{0}$, the clusters $\{\widehat{G}_{k}^{[K]}:1\leq k\leq K\}$ cannot provide an appropriate approximation of the true class structure $\{G_{k}:1\leq k\leq K_{0}\}$. In particular, there is always a cluster $S=\widehat{G}_{k}^{[K]}$ which contains subjects from at least two different classes. For simplicity, let $S=G_{k_{1}}\cup G_{k_{2}}$. For any $i\in S$, it holds that

[TABLE]

under our regularity conditions from Section 4.1. Moreover, it is not difficult to see that for at least one $i\in S$, $p^{-1/2}\sum\nolimits_{j=1}^{p}d_{ij}^{2}\geq c\sqrt{p}$ for some small constant $c>0$. This implies that for some $i\in S$,

[TABLE]

i.e., the statistic $\widehat{\Delta}_{i}^{[K]}$ has an explosive behaviour.

According to these heuristic considerations, the statistics $\widehat{\Delta}_{i}^{[K]}$ exhibit a quite different behaviour depending on whether $K<K_{0}$ or $K=K_{0}$. When $K<K_{0}$, the statistic $\widehat{\Delta}_{i}^{[K]}$ has an explosive behaviour at least for some subjects $i$. This mirrors the fact that a partition with $K<K_{0}$ clusters cannot give a reasonable approximation to the true class structure. In particular, it cannot describe the data of all subjects $i$ in an accurate way, resulting in an explosive behaviour of the (rescaled) residual sum of squares $\widehat{\Delta}_{i}^{[K]}$ for some subjects $i$. When $K=K_{0}$ in contrast, $\{\widehat{\Delta}_{i}^{[K_{0}]}:1\leq i\leq n\}$ is a collection of (approximately) independent random variables that (approximately) have a rescaled $\chi^{2}$-distribution. Hence, all statistics $\widehat{\Delta}_{i}^{[K_{0}]}$ have a stable, non-explosive behaviour. This reflects the fact that the partition $\{\widehat{G}_{k}^{[K_{0}]}:1\leq k\leq K_{0}\}$ is an accurate estimate of the true class structure and thus yields a moderate residual sum of squares $\widehat{\Delta}_{i}^{[K_{0}]}$ for all subjects $i$.

Since the statistics $\widehat{\Delta}_{i}^{[K]}$ behave quite differently depending on whether $K=K_{0}$ or $K<K_{0}$, they can be used to test $H_{0}:K=K_{0}$ versus $H_{1}:K<K_{0}$. In particular, testing $H_{0}$ versus $H_{1}$ can be achieved by testing the hypothesis that $\widehat{\Delta}_{i}^{[K]}$ are i.i.d. variables with a rescaled $\chi^{2}$-distribution against the alternative that at least one $\widehat{\Delta}_{i}^{[K]}$ has an explosive behaviour. We now construct a statistic $\widehat{\mathcal{H}}^{[K]}$ for this testing problem. A natural approach is to take the maximum of the individual statistics $\widehat{\Delta}_{i}^{[K]}$: Define

[TABLE]

and let $q(\alpha)$ be the $(1-\alpha)$-quantile of $\mathcal{H}=\max_{1\leq i\leq n}Z_{i}$, where $Z_{i}$ are independent random variables with the distribution $(\chi_{p}^{2}-p)/\sqrt{2p}$.

Our heuristic discussion from above, in particular formula (3.10), suggests that for $K=K_{0}$,

[TABLE]

Moreover, for $K<K_{0}$, we can show with the help of (3.11) and some additional considerations that $\widehat{\mathcal{H}}^{[K]}\geq c\sqrt{p}$ for some $c>0$ with probability tending to $1$. The quantile $q(\alpha)$, in contrast, can be shown to grow at the rate $\sqrt{\log n}$. Since $\sqrt{\log n}=o(\sqrt{p})$ under our conditions from Section 4.1, $\widehat{\mathcal{H}}^{[K]}$ diverges faster than the quantile $q(\alpha)$, implying that

[TABLE]

for $K<K_{0}$. This suggests that $\widehat{\mathcal{H}}^{[K]}$ has the property (3.4) and thus is a reasonable statistic to test the hypothesis $H_{0}:K=K_{0}$ versus $H_{1}:K<K_{0}$.

In this paper, we restrict attention to the maximum statistic $\widehat{\mathcal{H}}^{[K]}$ defined in (3.12). In principle though, we may work with any statistic that satisfies the higher-order property (3.4). In Section 6, we discuss some alternative choices of $\widehat{\mathcal{H}}^{[K]}$.

3.3 A $\boldsymbol{k}$-means clustering algorithm

We now construct a $k$-means type clustering algorithm which has the consistency property (3.1). Since its introduction by Cox (1957) and Fisher (1958), the $k$-means algorithm has become one of the most popular tools in cluster analysis. Our version of the algorithm mainly differs from the standard one in the choice of the initial values. To ensure the consistency property (3.1), we pick initial clusters $\mathscr{C}_{1}^{[K]},\ldots,\mathscr{C}_{K}^{[K]}$ for each given $K$ as follows:

Choice of the starting values. Let $i_{1},\ldots,i_{K}$ be indices which (with probability tending to $1$) belong to $K$ different classes $G_{k_{1}},\ldots,G_{k_{K}}$ in the case that $K\leq K_{0}$ and to $K_{0}$ different classes in the case that $K>K_{0}$. We explain how to obtain such indices below. With these indices at hand, we compute the distance measures $\widehat{\rho}_{k}(i)=\widehat{\rho}(i_{k},i)$ for all $1\leq i\leq n$ and $1\leq k\leq K$, where

[TABLE]

The starting values $\mathscr{C}_{1}^{[K]},\ldots,\mathscr{C}_{K}^{[K]}$ are now defined by assigning the index $i$ to cluster $\mathscr{C}_{k}^{[K]}$ if $\widehat{\rho}_{k}(i)=\min_{1\leq k^{\prime}\leq K}\widehat{\rho}_{k^{\prime}}(i)$.

The indices $i_{1},\ldots,i_{K}$ in this construction are computed as follows: For $K=2$, pick any index $i_{1}\in\{1,\ldots,n\}$ and calculate $i_{2}=\text{arg}\max_{1\leq i\leq n}\widehat{\rho}(i_{1},i)$. Next suppose we have already constructed the indices $i_{1},\ldots,i_{K-1}$ for the case of $K-1$ clusters and compute the corresponding starting values $\mathscr{C}_{1}^{[K-1]},\ldots,\mathscr{C}_{K-1}^{[K-1]}$ as described above. Calculate the maximal within-cluster distance $\widehat{\rho}_{\max}(k)=\max_{i\in\mathscr{C}_{k}^{[K-1]}}\widehat{\rho}_{k}(i)$ for each $1\leq k\leq K-1$ and let $\mathscr{C}_{k^{*}}^{[K-1]}$ be a cluster with $\widehat{\rho}_{\max}(k^{*})\geq\widehat{\rho}_{\max}(k)$ for all $k$. Define $i_{K}=\text{arg}\max_{i\in\mathscr{C}_{k^{*}}^{[K-1]}}\widehat{\rho}_{k^{*}}(i)$.

The $\boldsymbol{k}$-means algorithm. Let the number of clusters $K$ be given and denote the starting values by $C_{k}^{(0)}:=\mathscr{C}_{k}^{[K]}$ for $1\leq k\leq K$. The $r$-th iteration of our $k$-means algorithm proceeds as follows:

• Step $r$:

Let $C_{1}^{(r-1)},\ldots,C_{K}^{(r-1)}$ be the clusters from the $(r-1)$-th iteration step. Compute cluster means $m_{j,k}^{(r)}=(\#C_{k}^{(r-1)})^{-1}\sum\nolimits_{i\in C_{k}^{(r-1)}}\widehat{Y}_{ij}$ and calculate the distance measures $\widehat{\rho}_{k}^{(r)}(i)=p^{-1}\sum\nolimits_{j=1}^{p}(\widehat{Y}_{ij}-m_{j,k}^{(r)})^{2}$ for all $1\leq i\leq n$ and $1\leq k\leq K$. Define updated groups $C_{1}^{(r)},\ldots,C_{K}^{(r)}$ by assigning the index $i$ to the cluster $C_{k}^{(r)}$ if $\widehat{\rho}_{k}^{(r)}(i)=\min_{1\leq k^{\prime}\leq K}\widehat{\rho}_{k^{\prime}}^{(r)}(i)$.

Repeat this algorithm until the estimated groups do not change any more. For a given sample of data, this is guaranteed to happen after finitely many steps. The resulting $k$-means estimators are denoted by $\{\widehat{G}_{k}^{[K]}:1\leq k\leq K\}$. In Section 4, we formally show that these estimators have the consistency property (3.1) under our regularity conditions.

3.4 Estimation of $\boldsymbol{\sigma^{2}}$ and $\boldsymbol{\kappa}$

In practice, the error variance $\sigma^{2}$ and the normalization constant $\kappa$ are unknown and need to be estimated from the data at hand. We distinguish between two different estimation approaches, namely a difference- and a residual-based approach.

Difference-based estimators. To start with, consider the time trend model from Section 2, where $\mu_{ij}=\mu_{i}(j/p)$ with some trend function $\mu_{i}(\cdot)$. Supposing that the functions $\mu_{i}(\cdot)$ are Lipschitz continuous, we get that $Y_{ij}-Y_{i,j-1}=\{\varepsilon_{ij}-\varepsilon_{i,j-1}\}+\{\mu_{i}(j/p)-\mu_{i}((j-1)/p)\}=\{\varepsilon_{ij}-\varepsilon_{i,j-1}\}+O(p^{-1})$. This motivates to estimate the error variance $\sigma^{2}=\mathbb{E}[\varepsilon_{ij}^{2}]$ by

[TABLE]

Similarly, the fourth moment $\vartheta=\mathbb{E}[\varepsilon_{ij}^{4}]$ can be estimated by

[TABLE]

which in turn allows us to estimate the parameter $\kappa$ by

[TABLE]

Difference-based estimators of this type have been considered in the context of nonparametric regression by Müller et al. (1988) and Hall et al. (1990) among others. Under the technical conditions (C1)–(C3) from Section 4.1, it is straightforward to show that $\widehat{\sigma}^{2}_{\text{Lip}}=\sigma^{2}+O_{p}((np)^{-1/2}+p^{-2})$ and $\widehat{\kappa}_{\text{Lip}}=\kappa+O_{p}((np)^{-1/2}+p^{-2})$. The estimators $\widehat{\sigma}^{2}_{\text{Lip}}$ and $\widehat{\kappa}_{\text{Lip}}$ are particularly suited for applications where the trend functions $\mu_{i}(\cdot)$ can be expected to be fairly smooth. This ensures that the unknown first differences $(\varepsilon_{ij}-\varepsilon_{i,j-1})$ can be sufficiently well approximated by the terms $(Y_{ij}-Y_{i,j-1})$.

A similar difference-based estimation strategy can be used in the model for gene expression microarray data from Section 2. In this setting, the signal vectors $\boldsymbol{\mu}_{i}$ have a piecewise constant structure. In particular, $\mu_{ij}=\mu_{ij^{\prime}}$ for $j_{t-1}\leq j,j^{\prime}<j_{t}$, where $j_{t}$ are known indices with $1=j_{0}<j_{1}<\ldots<j_{T-1}<j_{T}=p+1$. This implies that $Y_{ij}-Y_{i,j-1}=\varepsilon_{ij}-\varepsilon_{i,j-1}$ for $j_{t-1}<j<j_{t}$. Similarly as before, we may thus estimate $\sigma^{2}$, $\vartheta$ and $\kappa$ by

[TABLE]

and $\widehat{\kappa}_{\text{pc}}=(\widehat{\vartheta}_{\text{pc}}/(\widehat{\sigma}_{\text{pc}}^{2})^{2}-1)^{1/2}$, where $\boldsymbol{1}_{j}=\boldsymbol{1}(j\notin\{j_{0},j_{1},\ldots,j_{T}\})$. It is not difficult to see that under the conditions (C1)–(C3), $\widehat{\sigma}^{2}_{\text{pc}}=\sigma^{2}+O_{p}((np)^{-1/2})$ and $\widehat{\kappa}_{\text{pc}}=\kappa+O_{p}((np)^{-1/2})$.

Residual-based estimators. Let $\{\widehat{G}_{k}^{[K]}:1\leq k\leq K\}$ be the $k$-means estimators from Section 3.3 for a given $K$. Moreover, let $\widehat{\varepsilon}_{ij}^{[K]}$ be the cluster-specific residuals introduced at the beginning of Section 3.2 and denote the vector of residuals for the $i$-th subject by $\widehat{\boldsymbol{\varepsilon}}_{i}^{[K]}=(\widehat{\varepsilon}_{i1}^{[K]},\ldots,\widehat{\varepsilon}_{ip}^{[K]})^{\top}$. With this notation at hand, we define the residual sum of squares for $K$ clusters by

[TABLE]

where $\|\cdot\|$ denotes the usual Euclidean norm for vectors. $\text{RSS}(K)$ can be shown to be a consistent estimator of $\sigma^{2}$ for any fixed $K\geq K_{0}$. The reason is the following: For any $K\geq K_{0}$, the $k$-means estimators $\widehat{G}_{k}^{[K]}$ have the property that

[TABLE]

for $1\leq k\leq K$ under the technical conditions (C1)–(C3) from Section 4.1. Hence, with probability tending to $1$, the estimated clusters $\widehat{G}_{k}^{[K]}$ contain elements from only one class $G_{k^{\prime}}$. The residuals $\widehat{\varepsilon}_{ij}^{[K]}$ should thus give a reasonable approximation to the unknown error terms $\varepsilon_{ij}$. This in turn suggests that the residual sum of squares $\text{RSS}(K)$ should be a consistent estimator of $\sigma^{2}$ for $K\geq K_{0}$.

Now suppose we know that $K_{0}$ is not larger than some upper bound $K_{\max}$. In this situation, we may try to estimate $\sigma^{2}$ by $\widetilde{\sigma}^{2}_{\text{RSS}}=\text{RSS}(K_{\max})$. Even though consistent, this is a very poor estimator of $\sigma^{2}$. The issue is the following: The larger $K_{\max}$, the smaller the residual sum of squares $\text{RSS}(K_{\max})$ tends to be. This is a natural consequence of the way in which the $k$-means algorithm works. Hence, if $K_{\max}$ is much larger than $K_{0}$, then $\widetilde{\sigma}^{2}_{\text{RSS}}=\text{RSS}(K_{\max})$ tends to strongly underestimate $\sigma^{2}$. To avoid this issue, we replace the naive estimator $\widetilde{\sigma}^{2}_{\text{RSS}}$ by a refined version:

1. (i)

Split the data vector $\boldsymbol{Y}_{i}=(Y_{i1},\ldots,Y_{ip})^{\top}$ into the two parts $\boldsymbol{Y}_{i}^{A}=(Y_{i1},Y_{i3},\ldots)^{\top}$ and $\boldsymbol{Y}_{i}^{B}=(Y_{i2},Y_{i4},\ldots)^{\top}$. Moreover, let $\overline{Y}_{i}^{A}$ be the empirical mean of the entries in the vector $\boldsymbol{Y}_{i}^{A}$ and define $\widehat{\boldsymbol{Y}}_{i}^{A}=(Y_{i1}-\overline{Y}_{i}^{A},Y_{i3}-\overline{Y}_{i}^{A},\ldots)^{\top}$. Finally, set $\mathcal{Y}^{A}=\{\widehat{\boldsymbol{Y}}_{i}^{A}:1\leq i\leq n\}$ and analogously define $\mathcal{Y}^{B}=\{\widehat{\boldsymbol{Y}}_{i}^{B}:1\leq i\leq n\}$. Importantly, under our technical conditions, the random vectors in $\mathcal{Y}^{A}$ are independent from those in $\mathcal{Y}^{B}$. 2. (ii)

Apply the $k$-means algorithm with $K=K_{\max}$ to the sample $\mathcal{Y}^{A}$ and denote the resulting estimators by $\{\widehat{G}_{k}^{A}:1\leq k\leq K_{\max}\}$. These estimators can be shown to have the property (3.14), provided that we impose the following condition: Let $\boldsymbol{m}_{k}$ be the class-specific signal vector of the class $G_{k}$ and define the vectors $\boldsymbol{m}_{k}^{A}$ and $\boldsymbol{m}_{k}^{B}$ in the same way as above. Assume that

[TABLE]

According to this assumption, the signal vectors $\boldsymbol{m}_{k}$ and $\boldsymbol{m}_{k^{\prime}}$ of two different classes can be distinguished from each other only by looking at their odd entries $\boldsymbol{m}_{k}^{A}$ and $\boldsymbol{m}_{k^{\prime}}^{A}$. It goes without saying that this is not a very severe restriction. 3. (iii)

Compute cluster-specific residuals from the data sample $\mathcal{Y}^{B}$,

[TABLE]

and define

[TABLE]

In contrast to the naive estimator $\widetilde{\sigma}^{2}_{\text{RSS}}$, the refined version $\widehat{\sigma}^{2}_{\text{RSS}}$ does not tend to strongly underestimate $\sigma^{2}$. The main reason is that the residuals $\widehat{\boldsymbol{\varepsilon}}_{i}^{B}$ are computed from the random vectors $\widehat{\boldsymbol{Y}}_{i}^{B}$ which are independent of the estimated clusters $\widehat{G}_{k}^{A}$.

Writing $\widehat{\boldsymbol{\varepsilon}}_{i}^{B}=(\widehat{\varepsilon}_{i1}^{B},\ldots,\widehat{\varepsilon}_{i\lfloor p/2\rfloor}^{B})^{\top}$, we can analogously estimate the fourth error moment $\vartheta=\mathbb{E}[\varepsilon_{ij}^{4}]$ by

[TABLE]

and set $\widehat{\kappa}_{\text{RSS}}=(\widehat{\vartheta}_{\text{RSS}}/(\widehat{\sigma}^{2}_{\text{RSS}})^{2}-1)^{1/2}$. Under the conditions (C1)–(C3), it can be shown that

[TABLE]

A sketch of the proof is given in the Supplementary Material.

4 Asymptotics

In this section, we investigate the asymptotic properties of our estimators. We first list the assumptions needed for the analysis and then summarize the main results.

4.1 Assumptions

To formulate the technical conditions that we impose on model (2.1)–(2.2), we denote the size, i.e., the cardinality of the class $G_{k}$ by $n_{k}=\#G_{k}$. Moreover, we use the shorthand $a_{\nu}\ll b_{\nu}$ to express that $a_{\nu}/b_{\nu}\leq c\nu^{-\delta}$ for sufficiently large $\nu$ with some $c>0$ and a small $\delta>0$. Our assumptions read as follows:

1. (C1)

The errors $\varepsilon_{ij}$ are identically distributed and independent across both $i$ and $j$ with $\mathbb{E}[\varepsilon_{ij}]=0$ and $\mathbb{E}[|\varepsilon_{ij}|^{\theta}]\leq C<\infty$ for some $\theta>8$. 2. (C2)

The class-specific signal vectors $\boldsymbol{m}_{k}=(m_{1,k},\ldots,m_{p,k})^{\top}$ differ across groups in the following sense: There exists a constant $\delta_{0}>0$ such that

[TABLE]

for any pair of groups $G_{k}$ and $G_{k^{\prime}}$ with $k\neq k^{\prime}$. Moreover, $|m_{j,k}|\leq C$ for all $k$ and $j$, where $C>0$ is a sufficiently large constant. 3. (C3)

Both $n$ and $p$ tend to infinity. The group sizes $n_{k}=\#G_{k}$ are such that $p\ll n_{k}\ll p^{(\theta/4)-1}$ for all $1\leq k\leq K_{0}$, implying that $p\ll n\ll p^{(\theta/4)-1}$.

We briefly comment on the above conditions. By imposing (C1), we restrict the noise terms $\varepsilon_{ij}$ to be i.i.d. Yet the error terms $e_{ij}=\alpha_{i}+\varepsilon_{ij}$ of our model may be dependent across subjects $i$, as we do not impose any restrictions on the random intercepts $\alpha_{i}$. This is important, for instance, when clustering the genes in a microarray data set, where we may expect different genes $i$ to be correlated. (C2) is a fairly harmless condition, which requires the signal vectors to differ in an $L_{2}$-sense across groups. By (C3), the group sizes $n_{k}$ and thus the total number of subjects $n$ are supposed to grow faster than the number of features $p$. We thus focus attention on applications where $n$ is (much) larger than $p$. However, $n$ should not grow too quickly as compared to $p$. Specifically, it should not grow faster than $p^{(\theta/4)-1}$, where $\theta$ is the number of existing error moments $\mathbb{E}[|\varepsilon_{ij}|^{\theta}]<\infty$. As can be seen, the bound $p^{(\theta/4)-1}$ on the growth rate of $n$ gets larger with increasing $\theta$. In particular, if all moments of $\varepsilon_{ij}$ exist, $n$ may grow as quickly as any polynomial of $p$. Importantly, (C3) allows the group sizes $n_{k}$ to grow at different rates (between $p$ and $p^{(\theta/4)-1}$). Put differently, it allows for strongly heterogeneous group sizes. Our estimation methods are thus able to deal with situations where some groups are much smaller than others.

4.2 Main results

Our first result shows that the maximum statistic $\widehat{\mathcal{H}}^{[K]}=\max_{1\leq i\leq n}\widehat{\Delta}_{i}^{[K]}$ has the property (3.4) and thus is a reasonable statistic to test the hypothesis $H_{0}:K=K_{0}$ versus $H_{1}:K<K_{0}$.

Theorem 4.1**.**

Assume that the estimated clusters have the consistency property (3.1). Moreover, let $\widehat{\sigma}^{2}$ and $\widehat{\kappa}$ be any estimators with $\widehat{\sigma}^{2}=\sigma^{2}+O_{p}(p^{-(1/2+\delta)})$ and $\widehat{\kappa}=\kappa+O_{p}(p^{-\delta})$ for some $\delta>0$. Under (C1)–(C3), the statistic $\widehat{\mathcal{H}}^{[K]}$ has the property (3.4), that is,

[TABLE]

This theorem is the main stepping stone to derive the central result of the paper, which describes the asymptotic properties of the estimators $\widehat{K}_{0}$ and $\{\widehat{G}_{k}:1\leq k\leq\widehat{K}_{0}\}$.

Theorem 4.2**.**

Under the conditions of Theorem 4.1, it holds that

[TABLE]

implying that $\mathbb{P}(\widehat{K}_{0}\neq K_{0})=\alpha+o(1)$. Moreover,

[TABLE]

Theorem 4.2 holds true for any clustering algorithm with the consistency property (3.1). The next result shows that this property is fulfilled, for example, by the $k$-means algorithm from Section 3.3.

Theorem 4.3**.**

Under (C1)–(C3), the $k$-means estimators $\{\widehat{G}_{k}^{[K_{0}]}:1\leq k\leq K_{0}\}$ from Section 3.3 satisfy (3.1), that is,

[TABLE]

The proofs of Theorems 4.1–4.3 can be found in the Supplementary Material.

5 Applications and Simulation Study

5.1 Clustering of temperature curves

Our first application is concerned with the analysis of a data set on land surface temperatures that was collected by the investigators of the Berkeley Earth project (Rohde et al., 2013). The data, which are publicly available at http://berkeleyearth.org/data, contain measurements on a grid of worldwide locations that is defined on a one degree (longitude) by one degree (latitude) basis. For each grid point, the data set contains a monthly land surface temperature profile. This profile is a vector with twelve entries, the first entry specifying the average temperature of all Januaries from 1951 to 1980, the second entry specifying the average temperature of all Februaries from 1951 to 1980, and so on. The temperature profiles at various example locations on earth are shown in Figure 1. As grid points containing $100\%$ sea surface are not taken into account, the overall number of grid points is equal to $n=\text{24,311}$. A detailed description of the derivation of the data can be found in Rohde et al. (2013). Our analysis is based on the Berkeley Earth source file from April 19, 2016.

The aim of our analysis is to cluster the 24,311 grid points in order to obtain a set of climate regions characterized by distinct temperature profiles. For this purpose, we impose the time trend model (2.4) on the data and apply the CluStErr method to them, setting $n=\text{24,311}$, $p=12$ and $\alpha=0.05$. To estimate the error variance $\sigma^{2}$ and the normalization constant $\kappa$, we apply the difference-based estimators $\widehat{\sigma}^{2}_{\text{Lip}}$ and $\widehat{\kappa}_{\text{Lip}}$ from Section 3.4, thus making use of the smoothness of the temperature curves illustrated in Figure 1.

The estimation results are presented in Figures 2 and 3. Figure 2 depicts the $p$-values $\widehat{p}^{[K]}$ corresponding to the test statistics $\widehat{\mathcal{H}}^{[K]}$ for different numbers of clusters $K$. It shows that the $p$-value $\widehat{p}^{[K]}$ remains below the $\alpha=0.05$ threshold for any $K<26$ but jumps across this threshold for $K=26$. The CluStErr algorithm thus estimates the number of clusters to be equal to $\widehat{K}_{0}=26$, suggesting that there are $26$ distinct climate regions. The sizes of the estimated clusters range between $244$ and 2,110; the error variance is estimated to be $\widehat{\sigma}^{2}=16.25$. Figure 3 uses a spatial grid to visualize the 26 regions and demonstrates the plausibility of the obtained results. For example, mountain ranges such as the Himalayas and the South American Andes, but also tropical climates in Africa, South America and Indonesia are easily identified from the plot. Of note, the results presented in Figure 3 show a remarkable similarity to the most recent modification of the Köppen-Geiger classification, which is one of the most widely used classification systems in environmental research (Peel et al., 2007). In particular, the overall number of climate regions defined in Peel et al. (2007) is equal to 29, which is similar to the cluster number $\widehat{K}_{0}=26$ identified by the CluStErr algorithm. Thus, although our example is purely illustrative, and although expert classification systems account for additional characteristics such as precipitation and vegetation, Figure 3 confirms the usefulness of the CluStErr method.

5.2 Clustering of gene expression data

Our second application is concerned with the analysis of gene expression data, which has become a powerful tool for the understanding of disease processes in biomedical research (Jiang et al., 2004). A popular approach to measure gene expression is to carry out microarray experiments. These experiments simultaneously quantify the expression levels of $n$ genes across $p$ samples of patients with different clinical conditions, such as tumor stages or disease subtypes. In the analysis of microarray data, clustering of the $n$ genes is frequently used to detect genes with similar cellular function and to discover groups of “co-expressed” genes showing similar expression patterns across clinical conditions (Chipman et al., 2003; Jiang et al., 2004; D’haeseleer, 2005).

In what follows, we analyze a set of gene expression data that was collected during the first stage of the Microarray Innovations in Leukemia (MILE) study (Haferlach et al., 2010). The data set contains expression level measurements for 20,172 genes and is publicly available as part of the Bioconductor package leukemiasEset (Aibar et al., 2013). The gene expression levels were measured using Affymetrix HG-U133 Plus 2.0 microarrays. For statistical analysis, the raw expression data were normalized using the Robust Multichip Average (RMA) method, followed by an additional gene-wise standardization of the expression levels. For details on data collection and pre-processing, we refer to Haferlach et al. (2010) and Aibar et al. (2013).

The data of the MILE study were obtained from $p=60$ bone marrow samples of patients that were untreated at the time of diagnosis. Of these patients, $48$ were either diagnosed with acute lymphoblastic leukemia (ALL, 12 patients), acute myeloid leukemia (AML, 12 patients), chronic lymphocytic leukemia (CLL, 12 patients), or chronic myeloid leukemia (CML, 12 patients). The other $12$ samples were obtained from non-leukemia (NoL) patients. From a biomedical point of view, the main interest focuses on the set of “differentially expressed” genes, that is, on those genes that show a sufficient amount of variation in their expression levels across the five tissue types (ALL, AML, CLL, CML, NoL). To identify the set of these genes, we run a univariate ANOVA for each gene and discard those with Bonferroni-corrected $p$-values $\geq 0.01$ in the respective overall $F$-tests. Application of this procedure results in a sample of $n=\text{3,167}$ univariately significant genes.

The aim of our analysis is to cluster the $n=\text{3,167}$ genes into groups whose members have similar expression patterns across the five tissue types (ALL, AML, CLL, CML, NoL). To do so, we impose model (2.5) from Section 2 on the data. The measured expression profiles $\boldsymbol{Y}_{i}=(Y_{i1},\ldots,Y_{ip})^{\top}$ of the various genes $i=1,\ldots,n$ are thus assumed to follow the model equation $\boldsymbol{Y}_{i}=\boldsymbol{\mu}_{i}+\boldsymbol{\alpha}_{i}+\boldsymbol{\varepsilon}_{i}$. The signal vectors $\boldsymbol{\mu}_{i}$ are supposed to have a piecewise constant structure after the patients have been ordered according to their tissue type (ALL, AML, CLL, CML, NoL). For illustration, the expression profile $\boldsymbol{Y}_{i}$ of a randomly selected gene is plotted in Figure 4.

To cluster the genes, we apply the CluStErr algorithm with the significance level $\alpha=0.05$ and the difference-based estimators $\widehat{\sigma}^{2}_{\text{pc}}$ and $\widehat{\kappa}_{\text{pc}}$ from Section 3.4, thus exploiting the piecewise constant structure of the signal vectors. The estimation results are presented in Figures 5 and 6. The plot in Figure 5 depicts the $p$-values $\widehat{p}^{[K]}$ corresponding to the test statistics $\widehat{\mathcal{H}}^{[K]}$ as a function of the cluster number $K$. It shows that the estimated number of clusters is $\widehat{K}_{0}=14$. The estimated sizes of the $14$ clusters range between $58$ and $469$. Moreover, the estimated error variance is $\widehat{\sigma}^{2}=0.442$. In Figure 6, the cluster centres $\boldsymbol{\widehat{m}}_{k}=(\#\widehat{G}_{k})^{-1}\sum\nolimits_{i\in\widehat{G}_{k}}\boldsymbol{\widehat{Y}}_{i}$ are presented, which estimate the cluster-specific signal vectors $\boldsymbol{m}_{k}=(\#G_{k})^{-1}\sum\nolimits_{i\in G_{k}}\boldsymbol{\mu}_{i}$. All clusters show a distinct separation of at least one tissue type, supporting the assumption of piecewise constant signals $\boldsymbol{m}_{k}$ and indicating that the genes contained in the clusters are co-expressed differently across the five groups. For example, cluster $\#2$ separates CML and NoL samples from ALL, AML and CLL samples, whereas cluster $\#4$ separates CLL samples from the other tissue types. Thus, each of the $14$ clusters represents a specific pattern of co-expressed gene profiles.

5.3 Simulation study

To explore the properties of the CluStErr method more systematically, we carry out a simulation study which splits into two main parts. The first part investigates the finite sample behaviour of CluStErr, whereas the second part compares CluStErr with several competing methods. The simulation design is inspired by the analysis of the gene expression data in Section 5.2. It is based on model (2.5) from Section 2. The data vectors $\boldsymbol{Y}_{i}$ have the form $\boldsymbol{Y}_{i}=\boldsymbol{\mu}_{i}+\boldsymbol{\varepsilon}_{i}$ with piecewise constant signal profiles $\boldsymbol{\mu}_{i}$. We set the number of clusters to $K_{0}=10$ and define the cluster-specific signal vectors $\boldsymbol{m}_{k}$ by

[TABLE]

where $\boldsymbol{1}=(1,\ldots,1)$ and $\boldsymbol{0}=(0,\ldots,0)$ are vectors of length $p/5$. A graphical illustration of the signal vectors $\boldsymbol{m}_{k}$ is provided in Figure 7. The error terms $\varepsilon_{ij}$ are assumed to be i.i.d. normally distributed with mean [math] and variance $\sigma^{2}$. In the course of the simulation study, we consider different values of $n$, $p$ and $\sigma^{2}$ as well as different cluster sizes. To assess the noise level in the simulated data, we consider the ratios between the error variance $\sigma^{2}$ and the “variances” of the signals $\boldsymbol{m}_{k}$. In particular, we define the noise-to-signal ratios $\text{NSR}_{k}(\sigma^{2})=\sigma^{2}/\textnormal{Var}(\boldsymbol{m}_{k})$, where $\textnormal{Var}(\boldsymbol{m}_{k})$ denotes the empirical variance of the vector $\boldsymbol{m}_{k}$. Since $\textnormal{Var}(\boldsymbol{m}_{k})\approx 0.16$ is the same for all $k$ in our design, we obtain that $\text{NSR}_{k}(\sigma^{2})=\text{NSR}(\sigma^{2})\approx\sigma^{2}/0.16$ for all $k$.

Finite sample properties of CluStErr. In this part of the simulation study, we analyze a design with equally sized clusters and set the sample size to $(n,p)=(1000,30)$. Three different noise-to-signal ratios NSR are considered, in particular $\textnormal{NSR}=1$, $1.5$ and $2$. Since $\sigma^{2}\approx 0.16\,\text{NSR}$, the corresponding error variances amount to $\sigma^{2}\approx 0.16$, $0.25$ and $0.32$, respectively. The noise-to-signal ratio $\textnormal{NSR}=1$ mimics the noise level in the application on gene expression data from Section 5.2, where the estimated noise-to-signal ratios all lie between $0.6$ and $1$. The ratios $\textnormal{NSR}=1.5$ and $\textnormal{NSR}=2$ are used to investigate how the CluStErr method behaves when the noise level increases. We implement the CluStErr algorithm with $\alpha=0.05$ and the difference-based estimators $\widehat{\sigma}^{2}_{\text{pc}}$ and $\widehat{\kappa}_{\text{pc}}$ from Section 3.4. For each of the three noise-to-signal ratios under consideration, we simulate $B=1000$ samples and compute the estimate $\widehat{K}_{0}$ for each sample.

The simulation results are presented in Figure 8. Each panel shows a histogram of the estimates $\widehat{K}_{0}$ for a specific noise-to-signal ratio. For the ratio level $\textnormal{NSR}=1$, the CluStErr method produces very accurate results: About $95\%$ of the estimates are equal to the true value $K_{0}=10$ and most of the remaining estimates take the value $11$. For the ratio level $\textnormal{NSR}=1.5$, the estimation results are also quite precise: Most of the estimates take a value between $9$ and $11$ with around $55\%$ of them being equal to the true value $K_{0}=10$. Only for the highest noise-to-signal ratio $\textnormal{NSR}=2$, the estimation results are less accurate. In this case, the noise level in the data is too high for the method to produce precise results. As one can see, the estimates have a strong downward bias, which can be explained as follows: When there is too much noise in the data, the test procedure on which the estimator $\widehat{K}_{0}$ is based does not have enough power to detect the alternative $H_{1}:K<K_{0}$. As a result, our repeated test procedure stops too soon, thus underestimating the true number of clusters.

According to our theoretical results, the estimator $\widehat{K}_{0}$ allows for statistical error control in the following sense: It has the property that $\mathbb{P}(\widehat{K}_{0}>K_{0})=\alpha+o(1)$ and $\mathbb{P}(\widehat{K}_{0}<K_{0})=o(1)$, implying that $\mathbb{P}(\widehat{K}_{0}=K_{0})=(1-\alpha)+o(1)$. Setting $\alpha$ to $0.05$, we should thus observe that $\widehat{K}_{0}$ equals $K_{0}=10$ in approximately $95\%$ of the simulations and overestimates $K_{0}$ in about $5\%$ of them. Table 1 shows that this is indeed the case for the lowest noise-to-signal ratio $\textnormal{NSR}=1$. In this situation, the probability $\mathbb{P}(\widehat{K}_{0}>K_{0})$ of overestimating $K_{0}$ is around $5\%$, while the probability $\mathbb{P}(\widehat{K}_{0}<K_{0})$ of underestimating $K_{0}$ is $0\%$, implying that $\mathbb{P}(\widehat{K}_{0}=K_{0})$ is about $95\%$. For the two higher ratio levels $\textnormal{NSR}=1.5$ and $2$, in contrast, the estimated values of the probabilities $\mathbb{P}(\widehat{K}_{0}<K_{0})$, $\mathbb{P}(\widehat{K}_{0}=K_{0})$ and $\mathbb{P}(\widehat{K}_{0}>K_{0})$ do not accurately match the values predicted by the theory. This is due to the fact that the statistical error control of the CluStErr method is asymptotic in nature. Table 2 illustrates this fact by reporting the estimated values of the probabilities $\mathbb{P}(\widehat{K}_{0}<K_{0})$, $\mathbb{P}(\widehat{K}_{0}=K_{0})$ and $\mathbb{P}(\widehat{K}_{0}>K_{0})$ for the noise-to-signal ratio $\textnormal{NSR}=1.5$ and various sample sizes $(n,p)=(1000,30)$, $(1500,40)$, $(2000,50)$, $(2500,60)$, $(3000,70)$. As one can clearly see, the estimated probabilities gradually approach the values predicted by the theory as the sample size increases.

To summarize, our simulations on the finite sample behaviour of the CluStErr method indicate the following: (i) The method produces accurate estimates of $K_{0}$ as long as the noise level in the data is not too high. (ii) For sufficiently large sample sizes, it controls the probability of under- and overestimating the number of clusters $K_{0}$ as predicted by the theory. (iii) For smaller sample sizes, however, the error control is not fully accurate.

It is important to note that (iii) is not a big issue: Even in situations where the error control is not very precise, the CluStErr method may still produce accurate estimates of $K_{0}$. This is illustrated by our simulations. Inspecting the histogram of Figure 8 with $\textnormal{NSR}=1.5$, for example, the estimated probability $\mathbb{P}(\widehat{K}_{0}=K_{0})$ is seen to be only around $55\%$ rather than $95\%$. Nevertheless, most of the estimates take a value between $9$ and $11$. Hence, in most of the simulations, the CluStErr method yields a reasonable approximation to the true number of clusters. From a heuristic perspective, the CluStErr method can indeed be expected to produce satisfying estimation results even in smaller samples when the error control is not very precise. This becomes clear when regarding CluStErr as a thresholding procedure. For $K=1,2,\ldots$, it checks whether the statistic $\widehat{\mathcal{H}}^{[K]}$ is below a certain threshold level $q$ and stops as soon as this is the case. For this approach to work, it is crucial to pick the threshold level $q$ appropriately. Our theoretical results suggest that the choice $q=q(\alpha)$ for a common significance level such as $\alpha=0.05$ should be appropriate. Of course, this choice guarantees precise error control only for sufficiently large sample sizes. Nevertheless, in smaller samples, the threshold level $q=q(\alpha)$ can still be expected to be of the right order of magnitude, thus resulting in reasonable estimates of $K_{0}$.

Comparison of CluStErr with competing methods. We now compare the CluStErr method to other criteria for selecting the number of clusters $K_{0}$, in particular to (i) the gap statistic (Tibshirani et al., 2001), (ii) the silhouette statistic (Rousseeuw, 1987), and (iii) the Hartigan index (Hartigan, 1975). As before, we set the sample size to $(n,p)=(1000,30)$ and consider the three noise-to-signal ratios $\textnormal{NSR}=1$, $1.5$ and $2$. In addition to a “balanced” scenario with clusters of the same size $n/K_{0}$ each, we also consider an “unbalanced” scenario with clusters of sizes $1+18k$ for $k=1,\ldots,K_{0}$. For each design, we simulate $B=100$ samples and compare the estimated cluster numbers obtained from the CluStErr method with those produced by the three competing algorithms.

The CluStErr estimates are computed as described in the first part of the simulation study. The three competing methods are implemented with a $k$-means algorithm as the underlying clustering method. To compute the values of the gap statistic, we employ the clusGap function contained in the R package cluster (Maechler et al., 2016). The number of clusters is estimated by the function maxSE with the option Tibs2001SEmax. We thus determine the number of clusters as suggested in Tibshirani et al. (2001). To compute the silhouette and Hartigan statistics, we apply the R package NbClust (Charrad et al., 2015).

The results of the comparison study are presented in Table 3. Part (a) of the table provides the results for the balanced scenario with equal cluster sizes. As can be seen, the CluStErr method clearly outperforms its competitors in the setting with $\textnormal{NSR}=1$. In the scenario with $\textnormal{NSR}=1.5$, it also performs well in comparison to the other methods. Only for the highest noise-to-signal ratio $\textnormal{NSR}=2$, it produces estimates of $K_{0}$ with a strong downward bias and is outperformed by the silhouette and Hartigan statistics. Part (b) of Table 3 presents the results for the unbalanced scenario where the clusters strongly differ in size. In this scenario, all of the three competing methods substantially underestimate the number of clusters. The CluStErr method, in contrast, provides accurate estimates of $K_{0}$ in the two designs with $\textnormal{NSR}=1$ and $\textnormal{NSR}=1.5$. Only in the high-noise design with $\textnormal{NSR}=2$, it produces estimates with a substantial downward bias, which nevertheless is much less pronounced than that of its competitors.

To summarize, the main findings of our comparison study are as follows: (i) The CluStErr method performs well in comparison to its competitors as long as the noise-to-signal ratio is not too high. It is however outperformed by its competitors in a balanced setting when the noise level is high. In Section 6, we discuss some modifications of the CluStErr method to improve its behaviour in the case of high noise. (ii) The CluStErr method is able to deal with both balanced and unbalanced cluster sizes, whereas its competitors perform less adequately in unbalanced settings.

The findings (i) and (ii) can heuristically be explained as follows: The CluStErr method is based on the test statistic $\widehat{\mathcal{H}}^{[K]}=\max_{1\leq i\leq n}\widehat{\Delta}_{i}^{[K]}$, which is essentially the maximum over the residual sums of squares of the various individuals $i$. Its competitors, in contrast, are based on statistics which evaluate averages rather than maxima. Hartigan’s rule, for instance, relies on a statistic which is essentially a scaled version of the ratio $\text{RSS}(K)/\text{RSS}(K+1)$, where $\text{RSS}(K)$ is defined as in (3.13) and denotes the average residual sum of squares for a partition with $K$ clusters. Averaging the residual sums of squares reduces the noise in the data much more strongly than taking the maximum. This is the reason why Hartigan’s rule tends to perform better than the CluStErr method in a balanced setting with high noise. On the other hand, the average residual sum of squares hardly reacts to changes in the residual sums of squares of a few individuals that form a small cluster. Hence, small clusters are effectively ignored when taking the average of the residual sums of squares. This is the reason why Hartigan’s statistic is not able to deal adequately with unbalanced settings. Taking the maximum of the residual sums of squares instead allows us to handle even highly unbalanced cluster sizes.

6 Extensions

In this paper, we have developed an approach for estimating the number of clusters with statistical error control. We have derived a rigorous mathematical theory for a model with convex spherical clusters. This model is widely used in practice and is suitable for a large number of applications. Nevertheless, it of course has some limitations. In particular, it is not suitable for applications where the clusters have non-convex shapes. An interesting question is how to extend our ideas to the case of general, potentially non-convex clusters. Developing theory for this general case is a very challenging problem. We have made a first step into this direction by providing theory for the case of spherical clusters.

There are several ways to modify and extend our estimation methods in the model setting (2.1)–(2.2). So far, we have based our methods on the maximum statistic $\widehat{\mathcal{H}}^{[K]}=\max_{1\leq i\leq n}\widehat{\Delta}_{i}^{[K]}$. However, we are not bound to this choice. Our approach can be based on any test statistic $\widehat{\mathcal{H}}^{[K]}$ that fulfills the higher-order property (3.4). The maximum statistic serves as a baseline which may be modified and improved in several directions. The building blocks of the maximum statistic are the individual statistics $\widehat{\Delta}_{i}^{[K]}$. Their stochastic behaviour has been analyzed in detail in Section 3.2. Under the null hypothesis $H_{0}:K=K_{0}$, the statistics $\widehat{\Delta}_{i}^{[K]}$ are approximately independent and distributed as $(\chi_{p}^{2}-p)/\sqrt{2p}$ variables. Under the alternative $H_{1}:K<K_{0}$ in contrast, they have an explosive behaviour at least for some $i$. This difference in behaviour suggests to test $H_{0}$ by checking whether the hypothesis

[TABLE]

holds for all subjects $i=1,\ldots,n$. We are thus faced with a multiple testing problem. A maximum statistic is a classical tool to tackle this problem. However, as is well known from the field of multiple testing, maximum statistics tend to be fairly conservative. When the noise level in the data is high, a test based on the maximum statistic $\widehat{\mathcal{H}}^{[K]}=\max_{1\leq i\leq n}\widehat{\Delta}_{i}^{[K]}$ can thus be expected to have low power against the alternative $H_{1}:K<K_{0}$. As a consequence, the repeated test procedure on which the estimator $\widehat{K}_{0}$ is based tends to stop too soon, thus underestimating the true number of clusters. This is exactly what we have seen in the high-noise scenarios of the simulation study from Section 5.3. We now present two ways how to construct a statistic $\widehat{\mathcal{H}}^{[K]}$ with better power properties.

A blocked maximum statistic. Let $w_{0}=0$ and define $w_{k}=\sum\nolimits_{r=1}^{k}\#\widehat{G}_{r}^{[K]}$ for $1\leq k\leq K$. Moreover, write $\widehat{G}_{k}^{[K]}=\{i_{w_{k-1}+1},\ldots,i_{w_{k}}\}$ with $i_{w_{k-1}+1}<\ldots<i_{w_{k}}$ for any $k$. To start with, we order the indices $\{1,\ldots,n\}$ clusterwise. In particular, we write them as $\{i_{1},i_{2},\ldots,i_{n}\}$, which yields the ordering

[TABLE]

We next partition the ordered indices into blocks

[TABLE]

where $N$ is the block length and $L=\lceil n/N\rceil$ is the number of blocks. With this notation at hand, we construct blockwise averages

[TABLE]

of the individual statistics $\widehat{\Delta}_{i}^{[K]}$ and aggregate them by taking their maximum, thus defining

[TABLE]

In addition, we let $q_{B}(\alpha)$ be the $(1-\alpha)$-quantile of

[TABLE]

where $Z_{i}$ are i.i.d. variables with the distribution $(\chi_{p}^{2}-p)/\sqrt{2p}$. Note that this definition of $\widehat{\mathcal{H}}^{[K]}_{B}$ nests the maximum statistic $\widehat{\mathcal{H}}^{[K]}=\max_{1\leq i\leq n}\widehat{\Delta}_{i}^{[K]}$ as a special case with the block length $N=1$.

Under appropriate restrictions on the block length $N$, the estimators that result from applying the CluStErr method with the blocked statistic $\widehat{\mathcal{H}}^{[K]}_{B}$ can be shown to have the theoretical properties stated in Theorems 4.1–4.3. More specifically, Theorems 4.1–4.3 can be shown to hold true for the blocked statistic $\widehat{\mathcal{H}}^{[K]}_{B}$ if the following two restrictions are satisfied: (i) $N/p^{1-\eta}=O(1)$ for some small $\eta>0$, that is, the block length $N$ diverges more slowly than $p$. (ii) $\#G_{k}/n\rightarrow c_{k}>0$ for all $k$, that is, the cluster sizes $\#G_{k}$ all grow at the same rate. Condition (ii) essentially rules out strongly differing cluster sizes. It is not surprising that we require such a restriction: To construct the blocked statistic $\widehat{\mathcal{H}}^{[K]}_{B}$, we average over the individual statistics $\widehat{\Delta}_{i}^{[K]}$. As already discussed in the context of the simulation study of Section 5.3, averaging has the effect that small clusters are effectively ignored. Hence, in contrast to the maximum statistic $\widehat{\mathcal{H}}^{[K]}=\max_{1\leq i\leq n}\widehat{\Delta}_{i}^{[K]}$, the blocked statistic $\widehat{\mathcal{H}}^{[K]}_{B}$ with a large block size $N$ can be expected not to perform adequately when the cluster sizes are highly unbalanced.

In balanced settings, however, the blocked statistic $\widehat{\mathcal{H}}^{[K]}_{B}$ can be shown to have better power properties than the maximum statistic when the block size $N$ is chosen sufficiently large. To see this, we examine the behaviour of $\widehat{\mathcal{H}}^{[K]}_{B}$ for different block lengths $N$. Our heuristic discussion of the individual statistics $\widehat{\Delta}_{i}^{[K]}$ from Section 3.2 directly carries over to the blocked versions $\widehat{\Lambda}_{\ell}^{[K]}$: With the help of (3.10), it is easy to see that

[TABLE]

Moreover, (3.11) together with some additional arguments suggests that $\widehat{\mathcal{H}}_{B}^{[K]}$ has an explosive behaviour for $K<K_{0}$. Specifically, under our technical conditions from Section 4.1 and the two additional restrictions (i) and (ii) from above, we can show that

[TABLE]

As the quantile $q_{B}(\alpha)$ grows at the slower rate $\sqrt{\log L}\ (\leq\sqrt{\log n})$, we can conclude that

[TABLE]

for $K<K_{0}$. As a result, $\widehat{\mathcal{H}}_{B}^{[K]}$ should satisfy the higher-order property (3.4). Moreover, according to (6.1), the statistic $\widehat{\mathcal{H}}_{B}^{[K]}$ explodes at the rate $\sqrt{Np}$ for $K<K_{0}$. Hence, the faster the block size $N$ grows, the faster $\widehat{\mathcal{H}}_{B}^{[K]}$ diverges to infinity. Put differently, the larger $N$, the more power we have to detect that $K<K_{0}$. This suggests to select $N$ as large as possible. According to restriction (i) from above, we may choose any $N$ with $N/p^{1-\eta}=O(1)$ for some small $\eta>0$. Ideally, we would thus like to pick $N$ so large that it grows at the same rate as $p^{1-\eta}$. In practice, we neglect the small constant $\eta>0$ and set $N=p$ as a simple rule of thumb.

According to the heuristic arguments from above, the blocked maximum statistic $\widehat{\mathcal{H}}^{[K]}_{B}$ with block length $N=p$ should be more powerful than the maximum statistic $\widehat{\mathcal{H}}^{[K]}=\max_{1\leq i\leq n}\widehat{\Delta}_{i}^{[K]}$ in settings with balanced cluster sizes. We examine this claim with the help of some simulations. To do so, we return to the balanced scenario of the comparison study in Section 5.3. For each of the data samples that were simulated for this scenario, we compute the CluStErr estimate of $K_{0}$ based on the blocked statistic $\widehat{\mathcal{H}}^{[K]}_{B}$ with $N=p$. Table 4 presents the results. It shows that the blocked CluStErr method yields accurate estimates of $K_{0}$ for all three noise-to-signal ratios. Comparing the results to those in Table 3(a), the blocked method can be seen to perform very well in comparison to the other procedures even in the high-noise setting with $\textnormal{NSR}=2$. This clearly shows the gain in power induced by the block structure of the statistic.

An FDR-based statistic. There are several approaches in the literature how to construct multiple testing procedures that have better power properties than the classical maximum statistic. Prominent examples are methods that control the false discovery rate (FDR) or the higher criticism procedure by Donoho and Jin (2004). We may try to exploit ideas from these approaches to construct a more powerful statistic $\widehat{\mathcal{H}}^{[K]}$. As an example, we set up a test statistic which uses ideas from the FDR literature: Rather than only taking into account the maximum of the statistics $\widehat{\Delta}_{i}^{[K]}$, we may try to exploit the information in all of the ordered statistics $\widehat{\Delta}_{(1)}^{[K]}\geq\ldots\geq\widehat{\Delta}_{(n)}^{[K]}$. In particular, following Simes (1986) and Benjamini and Hochberg (1995), we may set up our test for a given number of clusters $K$ as follows: Reject $H_{0}:K=K_{0}$ if

[TABLE]

where $q_{\chi}(\beta)$ is the $(1-\beta)$-quantile of the distribution $(\chi_{p}^{2}-p)/\sqrt{2p}$. This procedure can be rephrased as follows: Define the statistic

[TABLE]

and reject $H_{0}$ if $\widehat{\mathcal{H}}^{[K]}_{\text{FDR}}>1$. Instead of taking the maximum over the original statistics $\widehat{\Delta}_{i}^{[K]}$, we thus take the maximum over rescaled versions of the ordered statistics $\widehat{\Delta}_{(i)}^{[K]}$. Developing theory for the FDR-type statistic $\widehat{\mathcal{H}}^{[K]}_{\text{FDR}}$ is a very interesting topic which is however far from trivial.

Auxiliary results

In the proofs of Theorems 4.1–4.3, we frequently make use of the following uniform convergence result.

Lemma S.1**.**

Let $\mathcal{Z}_{s}=\{Z_{st}:1\leq t\leq T\}$ be sequences of real-valued random variables for $1\leq s\leq S$ with the following properties: (i) for each $s$, the random variables in $\mathcal{Z}_{s}$ are independent of each other, and (ii) $\mathbb{E}[Z_{st}]=0$ and $\mathbb{E}[|Z_{st}|^{\phi}]\leq C<\infty$ for some $\phi>2$ and $C>0$ that depend neither on $s$ nor on $t$. Suppose that $S=T^{q}$ with $0\leq q<\phi/2-1$. Then

[TABLE]

where the constant $\eta>0$ can be chosen as small as desired.

Proof of Lemma S.1. Define $\tau_{S,T}=(ST)^{1/\{(2+\delta)(q+1)\}}$ with some sufficiently small $\delta>0$. In particular, let $\delta>0$ be so small that $(2+\delta)(q+1)<\phi$. Moreover, set

[TABLE]

and write

[TABLE]

In what follows, we show that

[TABLE]

for any fixed constant $C>0$. Combining (S.1) and (S.2) immediately yields the statement of Lemma S.1.

We start with the proof of (S.1): It holds that

[TABLE]

where

[TABLE]

and

[TABLE]

for $S$ and $T$ sufficiently large, since

[TABLE]

This yields (S.1).

We next turn to the proof of (S.2): We apply the crude bound

[TABLE]

and show that for any $1\leq s\leq S$,

[TABLE]

where $C_{0}$ is a fixed constant and $\rho>0$ can be chosen as large as desired by picking $\eta$ slightly larger than $1/2-1/(2+\delta)$. Since $S=O(T^{q})$, this immediately implies (S.2). To prove (S.3), we make use of the following facts:

1. (i)

For a random variable $Z$ and $\lambda>0$, Markov’s inequality says that

[TABLE] 2. (ii)

Since $|Z_{st}^{\leq}/\sqrt{T}|\leq 2\tau_{S,T}/\sqrt{T}$, it holds that $\lambda_{S,T}|Z_{st}^{\leq}/\sqrt{T}|\leq 1/2$, where we set $\lambda_{S,T}=\sqrt{T}/(4\tau_{S,T})$. As $\exp(x)\leq 1+x+x^{2}$ for $|x|\leq 1/2$, this implies that

[TABLE] 3. (iii)

By definition of $\lambda_{S,T}$, it holds that

[TABLE]

Using (i)–(iii) and writing $\mathbb{E}(Z_{st}^{\leq})^{2}\leq C_{Z}<\infty$, we obtain that

[TABLE]

where $\rho>0$ can be chosen arbitrarily large if we pick $\eta$ slightly larger than $1/2-1/(2+\delta)$. ∎

Proof of Theorem 4.1

We first prove that

[TABLE]

To do so, we derive a stochastic expansion of the individual statistics $\widehat{\Delta}_{i}^{[K_{0}]}$.

Lemma S.2**.**

It holds that

[TABLE]

where

[TABLE]

and the remainder $R_{i}^{[K_{0}]}$ has the property that

[TABLE]

for some $\xi>0$.

The proof of Lemma S.2 as well as those of the subsequent Lemmas S.3–S.5 are postponed until the proof of Theorem 4.1 is complete. With the help of Lemma S.2, we can bound the probability of interest

[TABLE]

as follows: Since

[TABLE]

it holds that

[TABLE]

where

[TABLE]

As the remainder $R_{i}^{[K_{0}]}$ has the property (S.5), we further obtain that

[TABLE]

where

[TABLE]

With the help of strong approximation theory, we can derive the following result on the asymptotic behaviour of the probabilities $P_{\alpha}^{\ll}$ and $P_{\alpha}^{\gg}$.

Lemma S.3**.**

It holds that

[TABLE]

Together with (S.6), this immediately yields that $P_{\alpha}=(1-\alpha)+o(1)$, thus completing the proof of (S.4).

We next show that for any $K<K_{0}$,

[TABLE]

Consider a fixed $K<K_{0}$ and let $S\in\{\widehat{G}_{k}^{[K]}:1\leq k\leq K\}$ be any cluster with the following property:

[TABLE]

It is not difficult to see that a cluster with the property (S.8) must always exist under our conditions. By $\mathscr{C}\subseteq\{\widehat{G}_{k}^{[K]}:1\leq k\leq K\}$, we denote the collection of clusters that have the property (S.8). With this notation at hand, we can derive the following stochastic expansion of the individual statistics $\widehat{\Delta}_{i}^{[K]}$.

Lemma S.4**.**

For any $i\in S$ and $S\in\mathscr{C}$, it holds that

[TABLE]

where $d_{ij}=\mu_{ij}-(\#S)^{-1}\sum\nolimits_{i^{\prime}\in S}\mu_{i^{\prime}j}$ and the remainder $R_{i}^{[K]}$ has the property that

[TABLE]

for some small $\xi>0$.

Using (S.9) and the fact that

[TABLE]

we obtain that

[TABLE]

The arguments from the proof of Lemma S.3, in particular (S.22), imply that $q(\alpha)\leq C\sqrt{\log n}$ for some fixed constant $C>0$ and sufficiently large $n$. Moreover, we can prove the following result.

Lemma S.5**.**

It holds that

[TABLE]

for some fixed constant $c>0$.

Since $q(\alpha)\leq C\sqrt{\log n}$ and $\sqrt{\log n}/\sqrt{p}=o(1)$ by (C3), Lemma S.5 allows us to infer that

[TABLE]

Together with (S.10), this yields that $\mathbb{P}(\widehat{\mathcal{H}}^{[K]}\leq q(\alpha))=o(1)$. ∎

Proof of Lemma S.2. Let $n_{k}=\#G_{k}$ and write $\overline{\varepsilon}_{i}=p^{-1}\sum\nolimits_{j=1}^{p}\varepsilon_{ij}$ along with $\overline{\mu}_{i}=p^{-1}\sum\nolimits_{j=1}^{p}\mu_{ij}$. Since

[TABLE]

by (3.1), we can ignore the estimation error in the clusters $\widehat{G}_{k}^{[K_{0}]}$ and replace them by the true classes $G_{k}$. For $i\in G_{k}$, we thus get

[TABLE]

where

[TABLE]

We now show that $\max_{i\in G_{k}}|R_{i,\ell}^{[K_{0}]}|=o_{p}(p^{-\xi})$ for any $k$ and $\ell=A,\ldots,D$. This implies that $\max_{1\leq i\leq n}|R_{i,\ell}^{[K_{0}]}|=\max_{1\leq k\leq K_{0}}\max_{i\in G_{k}}|R_{i,\ell}^{[K_{0}]}|=o_{p}(p^{-\xi})$ for $\ell=A,\ldots,D$, which in turn yields the statement of Lemma S.2. Throughout the proof, we use the symbol $\eta>0$ to denote a sufficiently small constant which results from applying Lemma S.1.

By assumption, $\widehat{\sigma}^{2}=\sigma^{2}+O_{p}(p^{-(1/2+\delta)})$ and $\widehat{\kappa}=\kappa+O_{p}(p^{-\delta})$ for some $\delta>0$. Applying Lemma S.1 and choosing $\xi>0$ such that $\xi<\delta-\eta$, we obtain that

[TABLE]

and

[TABLE]

We next show that

[TABLE]

To do so, we work with the decomposition $R_{i,C}^{[K_{0}]}=\{2\widehat{\kappa}^{-1}\widehat{\sigma}^{-2}\}\{R_{i,C,1}^{[K_{0}]}+R_{i,C,2}^{[K_{0}]}-R_{i,C,3}^{[K_{0}]}\}$, where

[TABLE]

With the help of Lemma S.1, we obtain that

[TABLE]

Moreover,

[TABLE]

since

[TABLE]

$p\ll n_{k}$ and

[TABLE]

(S.14) is an immediate consequence of Lemma S.1. (S.15) follows upon observing that for any constant $C_{0}>0$,

[TABLE]

the last inequality resulting from the fact that the mean $\mathbb{E}[\varepsilon_{i\,j_{1}}\ldots\varepsilon_{i\,j_{4}}\varepsilon_{i_{1}^{\prime}j_{1}}\ldots\varepsilon_{i_{4}^{\prime}j_{4}}]$ can only be non-zero if some of the index pairs $(i_{\ell}^{\prime},j_{\ell})$ for $\ell=1,\ldots,4$ are identical. Finally, with the help of Lemma S.1, we get that

[TABLE]

Combining (S.12), (S.13) and (S.16), we arrive at the statement (S.11) on the remainder $R_{i,C}^{[K_{0}]}$.

We finally show that

[TABLE]

For the proof, we write $R_{i,D}^{[K_{0}]}=\{\widehat{\kappa}^{-1}\widehat{\sigma}^{-2}\}\{R_{i,D,1}^{[K_{0}]}+R_{i,D,2}^{[K_{0}]}\}$, where

[TABLE]

With the help of Lemma S.1, we obtain that

[TABLE]

Moreover, straightforward calculations yield that

[TABLE]

(S.17) now follows upon combining (S.18) and (S.19). ∎

Proof of Lemma S.3. We make use of the following three results:

1. (R1)

Let $\{W_{i}:1\leq i\leq n\}$ be independent random variables with a standard normal distribution and define $a_{n}=1/\sqrt{2\log n}$ together with

[TABLE]

Then for any $w\in\mathbb{R}$,

[TABLE]

In particular, for $w(\alpha\pm\varepsilon)=-\log(-\log(1-\alpha\pm\varepsilon))$, we get

[TABLE]

The next result is known as Khintchine’s Theorem.

1. (R2)

Let $F_{n}$ be distribution functions and $G$ a non-degenerate distribution function. Moreover, let $\alpha_{n}>0$ and $\beta_{n}\in\mathbb{R}$ be such that

[TABLE]

for any continuity point $x$ of $G$. Then there are constants $\alpha_{n}^{\prime}>0$ and $\beta_{n}^{\prime}\in\mathbb{R}$ as well as a non-degenerate distribution function $G_{*}$ such that

[TABLE]

at any continuity point $x$ of $G_{*}$ if and only if

[TABLE]

The final result exploits strong approximation theory and is a direct consequence of the so-called KMT Theorems; see Komlós et al. (1975, 1976):

1. (R3)

Write

[TABLE]

and let $F$ denote the distribution function of $X_{ij}$. It is possible to construct i.i.d. random variables $\{\widetilde{X}_{ij}:1\leq i\leq n,\,1\leq j\leq p\}$ with the distribution function $F$ and independent standard normal random variables $\{Z_{ij}:1\leq i\leq n,\,1\leq j\leq p\}$ such that

[TABLE]

have the following property:

[TABLE]

for some arbitrarily small but fixed $\delta>0$ and some constant $C>0$ that does not depend on $i$, $p$ and $n$.

We now proceed as follows:

1. (i)

We show that for any $w\in\mathbb{R}$,

[TABLE]

This in particular implies that

[TABLE]

where $w_{n}(\alpha\pm\varepsilon)=a_{n}w(\alpha\pm\varepsilon)+b_{n}$ with $a_{n}$, $b_{n}$ and $w(\alpha\pm\varepsilon)$ as defined in (R1). The proof of (S.20) is postponed until the arguments for Lemma S.3 are complete. 2. (ii)

The statement (S.21) in particular holds in the special case that $\varepsilon_{ij}\sim N(0,\sigma^{2})$. In this case, $q(\alpha)$ is the $(1-\alpha)$-quantile of $\max_{1\leq i\leq n}\Delta_{i}^{[K_{0}]}$. Hence, we have

[TABLE]

which implies that

[TABLE]

for sufficiently large $n$. 3. (iii)

Since $p^{-\xi}/a_{n}=p^{-\xi}\sqrt{2\log n}=o(1)$ by (C3), we can use (S.20) together with (R2) to obtain that

[TABLE]

As $w_{n}(\alpha-\varepsilon)-p^{-\xi}\leq q(\alpha)-p^{-\xi}\leq q(\alpha)+p^{-\xi}\leq w_{n}(\alpha+\varepsilon)+p^{-\xi}$ for sufficiently large $n$, it holds that

[TABLE]

for large $n$. Moreover, since $P_{\alpha,\varepsilon}^{\ll}\rightarrow 1-\alpha-\varepsilon$ and $P_{\alpha,\varepsilon}^{\gg}\rightarrow 1-\alpha+\varepsilon$ for any fixed $\varepsilon>0$ by (S.23), we can conclude that $P_{\alpha}^{\ll}=(1-\alpha)+o(1)$ and $P_{\alpha}^{\gg}=(1-\alpha)+o(1)$, which is the statement of Lemma S.3.

It remains to prove (S.20): Using the notation from (R3) and the shorthand $w_{n}=a_{n}w+b_{n}$, we can write

[TABLE]

with

[TABLE]

The probabilities $\pi_{i}$ can be decomposed into two parts as follows:

[TABLE]

where

[TABLE]

With the help of (R3) and the assumption that $n\ll p^{(\theta/4)-1}$, we can show that

[TABLE]

where $R_{n}$ is a non-negative remainder term with

[TABLE]

Moreover, the probabilities $\pi_{i}^{\leq}$ can be bounded by

[TABLE]

the second line making use of (R3). From this, we obtain that

[TABLE]

where

[TABLE]

By combining (S.24)–(S.26), we arrive at the intermediate result that

[TABLE]

Since $p^{\frac{1}{2+\delta}-\frac{1}{2}}/a_{n}=p^{\frac{1}{2+\delta}-\frac{1}{2}}\sqrt{2\log n}=o(1)$, we can use (R1) together with (R2) to show that

[TABLE]

Plugging (S.28) into (S.27) immediately yields that

[TABLE]

which completes the proof. ∎

Proof of Lemma S.4. We use the notation $n_{S}=\#S$ along with $\overline{\varepsilon}_{i}=p^{-1}\sum\nolimits_{j=1}^{p}\varepsilon_{ij}$, $\overline{\mu}_{i}=p^{-1}\sum\nolimits_{j=1}^{p}\mu_{ij}$ and $\overline{d}_{i}=p^{-1}\sum\nolimits_{j=1}^{p}d_{ij}$. For any $i\in S$ and $S\in\mathscr{C}$, we can write

[TABLE]

where

[TABLE]

We now show that $\max_{S\in\mathscr{C}}\max_{i\in S}|R_{i,\ell}^{[K]}|=o_{p}(p^{1/2-\xi})$ for $\ell=A,\ldots,G$. This immediately yields the statement of Lemma S.4. Throughout the proof, $\eta>0$ denotes a sufficiently small constant that results from applying Lemma S.1.

With the help of Lemma S.1 and our assumptions on $\widehat{\sigma}^{2}$ and $\widehat{\kappa}$, it is straightforward to see that $\max_{S\in\mathscr{C}}\max_{i\in S}|R_{i,\ell}^{[K]}|\leq\max_{1\leq i\leq n}|R_{i,\ell}^{[K]}|=o_{p}(p^{1/2-\xi})$ for $\ell=A,B,C,D$ with some sufficiently small $\xi>0$. We next show that

[TABLE]

To do so, we write $R_{i,E}^{[K]}=\{2\widehat{\kappa}^{-1}\widehat{\sigma}^{-2}\}\{R_{i,E,1}^{[K]}+R_{i,E,2}^{[K]}-R_{i,E,3}^{[K]}\}$, where

[TABLE]

Lemma S.1 yields that $\max_{S\in\mathscr{C}}\max_{i\in S}|R_{i,E,1}^{[K]}|\leq\max_{1\leq i\leq n}|R_{i,E,1}^{[K]}|=O_{p}(p^{2\eta}/\sqrt{p})$. Moreover, it holds that

[TABLE]

since

[TABLE]

and

[TABLE]

which follows upon applying Lemma S.1. Finally,

[TABLE]

which can again be seen by applying Lemma S.1. Putting everything together, we arrive at (S.29). Similar arguments show that

[TABLE]

as well.

To analyze the term $R_{i,G}^{[K]}$, we denote the signal vector of the group $G_{k}$ by $\boldsymbol{m}_{k}=(m_{1,k},\ldots,m_{p,k})^{\top}$ and write

[TABLE]

with $\lambda_{S,k}=\#(S\cap G_{k})/n_{S}$. With this notation, we get

[TABLE]

where

[TABLE]

With the help of Lemma S.1, it can be shown that \max_{S\in\mathscr{C}}\max_{i\in S}\big{|}R_{i,G,\ell}^{[K]}\big{|}=O_{p}(p^{\eta}) for $\ell=1,\ldots,5$. For example, it holds that

[TABLE]

As a result, we obtain that

[TABLE]

This completes the proof. ∎

Proof of Lemma S.5. Let $S\in\mathscr{C}$. In particular, suppose that $S\cap G_{k_{1}}\neq\emptyset$ and $S\cap G_{k_{2}}\neq\emptyset$ for some $k_{1}\neq k_{2}$. We show the following claim: there exists some $i\in S$ such that

[TABLE]

where $c=(\sqrt{\delta_{0}}/2)^{2}$ with $\delta_{0}$ defined in assumption (C2). From this, the statement of Lemma S.5 immediately follows.

For the proof of (S.32), we denote the Euclidean distance between vectors $v=(v_{1},\ldots,v_{p})^{\top}$ and $w=(w_{1},\ldots,w_{p})^{\top}$ by $d(v,w)=(\sum\nolimits_{j=1}^{p}|v_{j}-w_{j}|^{2})^{1/2}$. Moreover, as in Lemma S.4, we use the notation

[TABLE]

where $n_{S}=\#S$, $\lambda_{S,k}=\#(S\cap G_{k})/n_{S}$ and $\boldsymbol{m}_{k}=(m_{1,k},\ldots,m_{p,k})^{\top}$ is the signal vector of the class $G_{k}$.

Take any $i\in S\cap G_{k_{1}}$. If

[TABLE]

the proof is finished, as (S.32) is satisfied for $i$. Next consider the case that

[TABLE]

By assumption (C2), it holds that $d(\boldsymbol{m}_{k},\boldsymbol{m}_{k^{\prime}})\geq\sqrt{\delta_{0}p}$ for $k\neq k^{\prime}$. Hence, by the triangle inequality,

[TABLE]

implying that

[TABLE]

This shows that the claim (S.32) is fulfilled for any $i^{\prime}\in S\cap G_{k_{2}}$. ∎

Proof of Theorem 4.2

By Theorem 4.1,

[TABLE]

and

[TABLE]

Moreover,

[TABLE]

since

[TABLE]

by the consistency property (3.1) and

[TABLE]

∎

Proof of Theorem 4.3

With the help of Lemma S.1, we can show that

[TABLE]

uniformly over $i$ and $i^{\prime}$. This together with (C2) allows us to prove the following claim:

[TABLE]

Now let $K=K_{0}$. With the help of (S.33) and (S.34), the starting values $\mathscr{C}_{1}^{[K_{0}]},\ldots$ $\ldots,\mathscr{C}_{K_{0}}^{[K_{0}]}$ can be shown to have the property that

[TABLE]

Together with Lemma S.1, (S.35) yields that

[TABLE]

uniformly over $i$ and $k$. Combined with (C2), this in turn implies that the $k$-means algorithm converges already after the first iteration step with probability tending to $1$ and $\widehat{G}_{k}^{[K_{0}]}$ are consistent estimators of the classes $G_{k}$ in the sense of (3.1). ∎

Proof of (3.16)

Suppose that (C1)–(C3) along with (3.15) are satisfied. As already noted in Section 3.4, the $k$-means estimators $\{\widehat{G}_{k}^{A}:1\leq k\leq K_{\max}\}$ can be shown to satisfy (3.14), that is,

[TABLE]

for any $k=1,\ldots,K_{\max}$. This can be proven by very similar arguments as the consistency property (3.1). We thus omit the details. Let $E^{A}$ be the event that

[TABLE]

holds for all clusters $\widehat{G}_{k}^{A}$ with $k=1,\ldots,K_{\max}$. $E^{A}$ can be regarded as the event that the partition $\{\widehat{G}_{k}^{A}:1\leq k\leq K_{\max}\}$ is a refinement of the class structure $\{G_{k}:1\leq k\leq K_{0}\}$. By (S.36), the event $E^{A}$ occurs with probability tending to $1$.

Now consider the estimator

[TABLE]

Since the random variables $\widehat{\boldsymbol{Y}}_{i}^{B}$ are independent of the estimators $\widehat{G}_{k}^{A}$, it is not difficult to verify the following: for any $\delta>0$, there exists a constant $C_{\delta}>0$ (that does not depend on $\{\widehat{G}_{k}^{A}:1\leq k\leq K_{\max}\}$) such that on the event $E^{A}$,

[TABLE]

From this, the first statement of (3.16) easily follows. The second statement can be obtained by similar arguments. ∎