11affiliationtext: Department of Statistics, Federal University of Rio Grande
do Sul**affiliationtext: Corresponding author: Marcio Valk,
<EMAIL_ADDRESS>
# Clustering inference in multiple groups

Debora Zava Bello  Marcio Valk  Gabriela Bettella Cybis

(Jun 2021)

###### Abstract

Inference in clustering is paramount to uncovering inherent group structure in
data. Clustering methods which assess statistical significance have recently
drawn attention owing to their importance for the identification of patterns
in high dimensional data with applications in many scientific fields. We
present here a U-statistics based approach, specially tailored for high-
dimensional data, that clusters the data into _three_ groups while assessing
the significance of such partitions. Because our approach stands on the
U-statistics based clustering framework of the methods in R package $uclust$,
it inherits its characteristics being a non-parametric method relying on very
few assumptions about the data, and thus can be applied to a wide range of
dataset. Furthermore our method aims to be a more powerful tool to find the
best partitions of the data into three groups when that particular structure
is present. In order to do so, we first propose an extension of the test
U-statistic and develop its asymptotic theory. Additionally we propose a
ternary non-nested significance clustering method. Our approach is tested
through multiple simulations and found to have more statistical power than
competing alternatives in all scenarios considered. Applications to peripheral
blood mononuclear cells and to image recognition shows the versatility of our
proposal, presenting a superior performance when compared with other
approaches.

##  1 Introduction

In clusters analysis the aim is to divide data into groups of similar items
and there are different ways to accomplish this task. A large number of
algorithms based on different measures have been proposed and each different
measure may lead to potentially different results (Euan et al. (2019)).
Clusters can be inherently present in the data like in phylogenetic analysis
(Rosenberg et al. (2002); Chen et al. (2015)) or they can be built when
clustering should take place regardless of whether innate cluster structure is
present as in customer segmentation (Motlagh et al. (2019); Hennig (2015)). In
order to evaluate clustering methods, it is necessary to consider the context,
the objectives of clustering and to have a suitable measure of dissimilarity
(Von Luxburg et al. (2012)). A critical issue is how to discover inherent
cluster structure in data, in other words, whether the clusters represent in
fact an important feature or are simply the result of sample variation. This
becomes even more challenging when considering the context of high dimensional
data. We present here a U-statistics based approach that clusters the data in
three groups while assessing the significance of such partitions. Our method
is specially tailored for high-dimensional data and adaptable to different
distance measures.

In a typical application of inference in clustering when the groups are
already defined and there is no need for an algorithm or method to find them,
the null hypothesis is that all groups are random samples from the same
population (overall sample homogeneity). In the multivariate analysis of
variance (M)ANOVA procedure, when presented in terms of a linear model, the
homogeneity of groups stands for equality of means between all groups.
Assumptions of independence and normality of the data, homoscedasticity of
variance and homogeneity in group are required for exact (finite sample)
inference. In addition, a large sample size, depending on the dimension of the
data is generally necessary. For the context where there is no information
about the existence of groups and the objective is to know if they exist and
what they are, some approaches have been proposed for addressing the problem
of assessing significance of partitions, or determining which clustering
layers represent actual population structure and which are simple consequence
of spurious random effects. To avoid resorting to heuristic criteria or the
researcher’s judgement to define which partition levels should be assigned
meaning these approaches proposes to assess statistical significance. However
the success of these methods depends on the underlying cluster structure
(Adolfsson et al. (2019)).

Several approaches have been proposed to assess statistical significance in
clustering, for example the procedure presented in McLachlan and Peel (2004)
which considers mixture models of distributions such as the Gaussian. A
maximum likelihood approach is used by Demidenko (2018) to test no-clusters
hypothesis. However, when the data are high dimensional and have small sample
sizes the problem becomes increasingly challenging, since it involves complete
parametric estimation, usually requiring costly matrix inversions. The works
of McShane et al. (2002); Helgeson et al. (2020) address this issue by using
reduction of dimensionality of the data matrix and sparse covariance
estimation. An approach inspired on the bootstrap strategy is proposed by
Shimodaira et al. (2004) which is implemented in the R package $pvclust$
(Suzuki and Shimodaira (2006)) and used in phylogenetics to assess confidence
in hierarchical clustering. Liu et al. (2008) proposes a statistical test to
assess the significance of clustering the data into $K$ groups, specifically
tailored to the high dimension low sample size (HDLSS) scenario, that has been
implemented in the R package $sigclust$. However, the implementation and
applications consider only two groups. Additionally, Kimes et al. (2017)
extend the method to assess significance in hierarchical clustering. However,
this approach requires that the data comes from a single multivariate normal
distribution, which can be an issue since rejection of the no cluster
hypothesis may be a simple consequence of non-normal data.

Our work focuses specifically on the HDLSS setting and extends the works of
Cybis et al. (2018); Valk and Cybis (2020) making it possible to
simultaneously test the homogeneity of _three_ groups, one of which may have
size _one_. The test statistic to compare _three_ groups, where one of them
may be an outlier, is a extension of the test statistic $B_{n}$ proposed by
Pinheiro et al. (2009). Here the hypotheses are similar to those of (M)ANOVA
where the null is that the elements in the _three_ groups come from the same
distribution (homogeneity, no-clusters) versus the alternative hypothesis that
the data distribution (not necessarily normal) of at least one of the groups
is different from the others. Asymptotic normality of the extended $B_{n}$ is
obtained using U-statistics theory. An estimator for the variance of the
extended $B_{n}$ is proposed. In addition, we have developed an algorithm
($uclust3$) that finds the best significant separation in _three_ groups.
Simulation studies show that our proposal presents coherent results, such as
control of Type I Error and the increased Power to identify clusters as they
become more separated. Furthermore, our comparative simulation study with
other methods shows that in the case where there are exactly _three_ groups,
the approach we are proposing has greater power, that is, greater ability to
correctly identify _three_ clusters. More accurate results of $uclust3$ are
found in an application to real image recognition data, corroborating the
better performance of our approach observed in the simulations. Although we
are using Euclidean distance and simulating data with normal distribution,
these aspects are not essential to the validity of the method properties.

The steps to developing our three groups clustering method are outlined as
follows. First, in Section 2.1 we review the U-statistics based theory of the
homogeneity test of Cybis et al. (2018) and present the U-statistics theory
for _three_ groups. In Section 2.2 we present the extension of the $B_{n}$
statistics proposed by Pinheiro et al. (2009) to contemplate _three_ groups in
which one may have size one, in order to devise a clustering algorithm that
can properly identify outlier elements. Additionally an investigation of
theoretical properties that show its compatibility with the previous framework
and asymptotic theory, is also presented. In Section 2.3 we explore the
variance aspects of the extended $B_{n}$ and propose an approach to estimate
this variance. In Section 3 we propose the $uclust3$ method which finds the
statistically significant data partition that better separates the sample into
three groups. The remainder of the paper focuses on evaluating the methodology
through simulation studies, in Section 4, and applications to real data in
Section 5. Finally, in Section 6 we discuss the overall results.

##  2 Methods

###  2.1 U-Statistics based test for three group separation

Let $\bf{X}=({\bf X}_{1},\dots,{\bf X}_{n})$ be a random sample of $n$
$L$-dimensional vectors divided in three groups $G_{1}$, $G_{2}$ and $G_{3}$
of sample sizes $n_{1}$, $n_{2}$ and $n_{3}$, respectively, where
$n=n_{1}+n_{2}+n_{3}$. In the $g$-th group, for $g\in\\{1,2,3\\}$,
observations ${\bf X}^{(g)}_{1},\dots,{\bf X}^{(g)}_{n_{g}}$ are assumed to be
independent and identically distributed with a $L$-variate distribution
$F_{g}$. Here, the distribution $F_{g}$ admits finite mean vector
$\bm{\mu}_{g}$ and positive definite dispersion matrix $\bm{\Sigma}_{g}$ (not
necessarily multi-normal). Following the approach of Sen (2006) and Pinheiro
et al. (2009), we define the functional distance
$\theta(F_{g},F_{g^{\prime}})$ as

$\theta(F_{g},F_{g^{\prime}})=\int\int\phi(x_{1},x_{2})dF_{g}(x_{1})dF_{g^{\prime}}(x_{2}),\quad
x_{1},x_{2}\in\mathbb{R}^{L},$ (1)

where $g,g^{\prime}\in\\{1,2,3\\}$ and $\phi(\cdot,\cdot)$ is a symmetric
kernel of order $2$. If we assume that $\theta(\cdot,\cdot)$ is a convex
linear function of its marginal components, then we have

$\theta(F_{g},F_{g^{\prime}})\geq\frac{1}{2}\,\\{\theta(F_{g},F_{g})+\theta(F_{g^{\prime}},F_{g^{\prime}})\\},$
(2)

for all distributions $F_{g}$ and $F_{g^{\prime}}$, with equality holding
whenever $\mu_{g}=\mu_{g^{\prime}}$.

Note that the functional $\theta(\cdot,\cdot)$ can be used to define both
distance within and between groups. It follows from U-statistics theory that
an unbiased estimator of this functional for within group distance
$\theta(F_{g},F_{g})$ is a generalized U-statistic Hoeffding (1948), with
kernel $\phi(\cdot,\cdot)$, defined as

$U_{n_{g}}^{(g)}=\dbinom{n_{g}}{2}^{-1}\sum_{1\leq i<j\leq n_{g}}\phi({\bf
X}^{(g)}_{i},{\bf X}^{(g)}_{j}),$ (3)

where $g\in\\{1,2,3\\}$. Analogously, the unbiased estimator for the between
group functional distance $\theta(F_{g},F_{g^{\prime}})$ is defined by

$U_{n_{g},n_{g^{\prime}}}^{(g,g^{\prime})}=\frac{1}{n_{g}n_{g^{\prime}}}\sum_{i=1}^{n_{g}}\sum_{j=1}^{n_{g^{\prime}}}\phi({\bf
X}^{(g)}_{i},{\bf X}^{(g^{\prime})}_{j}),$ (4)

where $g,g^{\prime}\in\\{1,2,3\\}$ and $g\neq g^{\prime}$.

The combined sample U-statistic is usually decomposed as

$\displaystyle U_{n}$ $\displaystyle=$
$\displaystyle\sum_{g=1}^{3}\frac{n_{g}}{n}U_{n_{g}}^{(g)}+\sum_{1\leq
g<g^{\prime}\leq
3}\frac{n_{g}n_{g^{\prime}}}{n(n-1)}\left\\{2U_{n_{g},n_{g^{\prime}}}^{(g,g^{\prime})}-U_{n_{g}}^{(g)}-U_{n_{g^{\prime}}}^{(g^{\prime})}\right\\}$
(5) $\displaystyle=$ $\displaystyle W_{n}+B_{n}.$

Decomposition (5) leads to the statistic $B_{n}$, which provides the focal
point of our methodology,

$\displaystyle B_{n}$ $\displaystyle=$ $\displaystyle\sum_{1\leq
g<g^{\prime}\leq
3}\frac{n_{g}n_{g^{\prime}}}{n(n-1)}\left\\{2U_{n_{g},n_{g^{\prime}}}^{(g,g^{\prime})}-U_{n_{g}}^{(g)}-U_{n_{g^{\prime}}}^{(g^{\prime})}\right\\}.$
(6)

Here $U_{n_{g}}^{(g)}$ for $g\in\\{1,2,3\\}$ are U-statistics associated to
within group distances, as defined in (3), and
$U_{n_{g}n_{g^{\prime}}}^{(g,g^{\prime})}$, $g\neq g^{\prime}\in\\{1,2,3\\}$,
are the U-statistics associated to between group distances as defined in (4).
Note that the definition of $U_{n_{g}}^{(g)}$ require a minimum of 2 elements
in the group. This imposes minimum group sizes $n_{g}\geq 2$, for
$g\in\\{1,2,3\\}$ for proper definition of $B_{n}$.

The methodology proposed in Cybis et al. (2018) and Valk and Cybis (2020)
considers a group homogeneity test which verifies whether two groups in fact
constitute separated groups, or if they stem from the same distribution. In
this work, for data arranged in three groups $G_{1}$, $G_{2}$ and $G_{3}$, the
interest is in verifying whether the data are homogeneous or if there is at
least one group statistically separated. Thus, the null hypothesis $H_{0}$
states that $F_{1}=F_{2}=F_{3}$, while the alternative $H_{1}$ states that
there are $i\neq j$, $\in\\{1,2,3\\}$ where $F_{i}\neq F_{j}$. In cases where
groups $G_{1}$, $G_{2}$ and $G_{3}$ have more than two elements, the
asymptotic properties of $B_{n}$ are addressed in Pinheiro et al. (2009). The
statistics $B_{n}$ is in the class of degenerate U-statistics for which
asymptotic normality prevails and the convergence rates are $L$ and/or
$\sqrt{n}$. Additionally, under the null, we have $\mathbb{E}(B_{n})=0$ and
under the alternative, $\mathbb{E}(B_{n})>0$. The null hypothesis is rejected
for large values of standardized $B_{n}$, where the variance of $B_{n}$, under
$H_{0}$, is obtained by a resampling procedure Sen (2006).

###  2.2 The extension of test U-statistics for tree groups

The homogeneity test proposed in Cybis et al. (2018) presents an essential
concept for our clustering algorithm. However, the group size restriction
required by the definition of the U-statistic $B_{n}$ in (6) constrains this
method to cases where all subgroups have sizes $n_{i}\geq 2$, $i=1,2,3$, and
consequently clustering methods will fail in cases where the data has an
outlier. In order to build a clustering algorithm that admits groups of size 1
we propose an extension of $B_{n}$. We can assume, without loss of generality,
that only the group $G_{1}$ may have one element, and define

$\displaystyle
B_{n}=\left\\{\begin{array}[]{ll}\frac{2n_{2}}{n(n-1)}\left(U_{1,n_{2}}^{(1,2)}-U_{n_{2}}^{(2)}\right)+\frac{2n_{3}}{n(n-1)}\left(U^{(1,3)}_{1,n_{3}}-U_{n_{3}}^{(3)}\right)\\\
\\\
+\frac{n_{2}n_{3}}{n(n-1)}\left(2U_{n_{2},n_{3}}^{(2,3)}-U_{n_{2}}^{(2)}-U_{n_{3}}^{(3)}\right),\hskip
14.22636pt\hbox{ if }n_{1}=1,\hbox{ and }n_{2},n_{3}>1\\\ \\\
\displaystyle\sum_{1\leq i<j\leq
3}\frac{n_{i}n_{j}}{n(n-1)}\left(2U_{n_{i},n_{j}}^{(i,j)}-U_{n_{i}}^{(i)}-U_{n_{j}}^{(j)}\right),\hskip
14.22636pt\hbox{if }n_{1},n_{2},n_{3}>1.\\\ \end{array}\right.$ (12)

where $U_{n_{g},n_{g^{\prime}}}^{(g,g^{\prime})}$ and $U_{n_{g}}^{(g)}$ are
defined, respectively, in (4) and (3).

This is a natural extension of $B_{n}$ considering data separation in three
groups, when allowing for clusters of size 1. This extension coincides with
that of expression (6) for group of sizes $n_{1},n_{2},n_{3}>1$, and thus all
properties mentioned above are still valid for the new definition in that
case. We ascertain the validity of these asymptotic properties or analogous
alternatives in the case of $n_{1}=1$.

Note that, when $G_{1}$ has size one, we can rewrite $B_{n}$ as

$\displaystyle B_{n}$ $\displaystyle=$
$\displaystyle\frac{2n_{2}}{n(n-1)}U^{(1,2)}_{1,n_{2}}+\frac{2n_{3}}{n(n-1)}U^{(1,3)}_{1,n_{3}}+\frac{2n_{2}n_{3}}{n(n-1)}U^{(2,3)}_{n_{2},n_{3}}$
$\displaystyle-\frac{n_{2}(2+n_{3})}{n(n-1)}U^{(2)}_{n_{2}}-\frac{n_{3}(2+n_{2})}{n(n-1)}U^{(3)}_{n_{3}}$

where $U_{1,g}^{(1,g)}$ and $U_{n_{g}}^{(g)}$, $g=2,3$ are as defined in (4)
and (3). If we consider the extension of $B_{n}$ in (12), then we can write
the combined sample U-statistics as

$\displaystyle U_{n}$ $\displaystyle=$ $\displaystyle B_{n}+W_{n}^{*}.$

where $W_{n}^{*}$ is an appropriate modification the term $W_{n}$. Thus,
$B_{n}$ still arises from the decomposition of the combined sample
U-statistics into $B_{n}$ and a modified term $W_{n}$. This extended
definition allows us to build a U-test when a group has size 1. We
conveniently labeled the data in order to arrange the groups as follows. Let
$G_{1}=\\{{\bf X}_{1}\\}$, $G_{2}=\\{{\bf X}_{2},\dots,{\bf X}_{n_{2}+1}\\}$
and $G_{3}=\\{{\bf X}_{n_{2}+2},\dots,{\bf X}_{n}\\}$, $n=1+n_{2}+n_{3}$. We
still have $\mathbb{E}[B_{n}]=0$, under the null hypothesis of overall group
homogeneity. Additionally, if we make the assumption that

$\displaystyle\theta_{gg^{\prime}}>\theta_{g},$ (13)

for $g\neq g^{\prime}\,\in\\{1,2,3\\}$ where
$\theta_{g}=\operatorname{E}\left[\phi(X_{g},X_{g})\right]$ and
$\theta_{gg^{\prime}}=\operatorname{E}\left[\phi(X_{g},X_{g^{\prime}})\right]$,
then under alternative we have that $\mathbb{E}[B_{n}]>0$. Note that this
assumption is usual and when (13) is valid then equation (2) is always
satisfied.

Asymptotic theory for the $B_{n}$ statistic for group sizes greater than 2 is
developed in the work of Pinheiro et al. (2009), where it is established that
$B_{n}$ is a degenerate U-statistic and asymptotic normality is provided. The
following theorems demonstrate that the extended $B_{n}$ is a non degenerated
U-statistics and establish the asymptotic distribution of the extended $B_{n}$
under $H_{0}$ for increasing dimension $L$ and sample size $n$, requiring
regularity conditions akin to those of the $n_{1},n_{2},n_{3}>1$ case. The
following Lemma is an important result required to demonstrate the asymptotic
convergence of the test statistic.

###### Lemma 2.1

Let $\frac{X}{\delta_{n}}\xrightarrow{D}\operatorname{N}(0,1)$,
$\delta_{n}=\operatorname{O}(1)$ and $\delta^{*}_{n}=\operatorname{O}(1)$.
Then, $\frac{X}{\delta^{*}_{n}}\xrightarrow{D}\operatorname{N}(0,M)$ where
$M=\lim_{n\rightarrow\infty}\left(\frac{\delta_{n}^{2}}{\delta_{n}^{*2}}\right)$.

Proof:  Note that

$\frac{X}{\delta^{*}_{n}}\frac{\delta_{n}}{\delta_{n}}=\frac{\delta_{n}}{\delta^{*}_{n}}\frac{X}{\delta_{n}}\xrightarrow{D}\operatorname{N}(0,\gamma),$

where

$\gamma=\operatorname{Var}\left(\frac{\delta_{n}}{\delta^{*}_{n}}\frac{X}{\delta_{n}}\right)\rightarrow\lim_{n\rightarrow\infty}\left(\frac{\delta_{n}}{\delta^{*}_{n}}\right)^{2}=M.$

###### Theorem 1

Let ${\bf X}_{1},{\bf X}_{2},\dots,{\bf X}_{n}$ be a sequence of i.i.d.
$L\times 1$ random vectors. Let $\phi(\cdot,\cdot)$ be a kernel of degree 2
satisfying $\mathbb{E}[\phi({\bf X}_{1},{\bf X}_{2})^{2}]<\infty$ and
$\operatorname{Var}[\mathbb{E}(\phi({\bf X}_{1},{\bf X}_{2})|{\bf
X}_{1})]=\sigma_{1}^{2}>0$. Consider definition (12) for $B_{n}$ when
$n_{1}=1$ and let $V_{n}=\operatorname{Var}(B_{n})$,
$\tau_{n}=(n/2)V_{n}^{1/2}$ and $W=J_{1}+J_{2}-J_{3}-J_{4}$, where
$\frac{\psi_{1}(X_{1})}{\tau_{n}}\xrightarrow{D}J_{1}$, and $J_{2},J_{3}\hbox{
and }J_{4}$ are random variables with normal distribution. Then

$\frac{(n/2)B_{n}}{\tau_{n}}\xrightarrow{D}W\hbox{ as n }\rightarrow\infty.$
(14)

Proof:

We are interested in the distribution of $B_{n}$ with fixed $L$ and
$n\rightarrow\infty$. Is is strightforward to show that
$\tau_{n}=\frac{n}{2}\sqrt{\operatorname{Var}(B_{n})}=O(1)$. From the
Hoeffding decomposition of $B_{n}$ we have:

$\frac{n}{2}B_{n}=W_{1}+W_{2}-W_{3}-W_{4}$ (15)

where

$\displaystyle W_{1}$ $\displaystyle=$
$\displaystyle\psi_{1}(X_{1})-\frac{1}{n-1}\sum_{i=1}^{n_{2}}\psi_{1}(X_{2i})-\frac{1}{n-1}\sum_{j=1}^{n_{3}}\psi_{!}(X_{3j})+$
(16)
$\displaystyle+\frac{1}{n-1}\sum_{i=1}^{n_{2}}\psi_{2}(X_{1},X_{2i})+\frac{1}{n-1}\sum_{j=1}^{n_{3}}\psi_{2}(X_{1},X_{3j})$
$\displaystyle W_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{n-1}\sum_{i=1}^{n_{2}}\sum_{j=1}^{n_{3}}\psi_{2}(X_{2i},X_{3j})$
(17) $\displaystyle W_{3}$ $\displaystyle=$
$\displaystyle\frac{2+n_{3}}{(n-1)(n_{2}-1)}\sum_{1\leq i<j\leq
n_{2}}\psi_{2}(X_{2i},X_{2j})$ (18) $\displaystyle W_{4}$ $\displaystyle=$
$\displaystyle\frac{2+n2}{(n-1)(n_{3}-1)}\sum_{1\leq i<j\leq
n_{3}}\psi_{2}(X_{3i},X_{3j}).$ (19)

Under the null hypothesis $\bm{X}_{1},\hbox{ }\bm{X}_{2}\hbox{ and
}\bm{X}_{3}$ are identically distributed, thus $W_{1}$ can be expressed as

$\displaystyle
W_{1}=\psi_{1}(X_{1})-\frac{1}{n-1}\displaystyle\sum_{i=2}^{n}\psi_{1}(X_{i})+\frac{1}{n-1}\displaystyle\sum_{j=2}^{n}\psi_{2}(X_{1},X_{j}).$
(20)

By the Law of Large Numbers (LLN) follows that

$\displaystyle\frac{1}{n-1}\displaystyle\sum_{i=2}^{n}\psi_{1}(X_{i})\xrightarrow{P}\operatorname{E}[\psi_{1}(X_{1})]=0$
(21)
$\displaystyle\frac{1}{n-1}\displaystyle\sum_{j=2}^{n}\psi_{2}(X_{1},X_{j})\xrightarrow{P}\operatorname{E}[\psi_{2}(X_{1},X_{2})]=0.$
(22)

Thereby,

$W_{1}\xrightarrow{P}\psi_{1}(X_{1}).$

As $\frac{\psi_{1}(X_{1})}{\tau_{n}}\xrightarrow{D}J_{1}$ and
$W_{1}\xrightarrow{P}\psi_{1}(X_{1})$, then, by Slutsky’s theorem,
$\frac{W_{1}}{\tau_{n}}\xrightarrow{D}J_{1}\hbox{ as n }\rightarrow\infty$.

From the Central Limit Theorem (TCL) we have

$\displaystyle\frac{W_{2}-E(W_{2})}{\sqrt{\operatorname{Var}(W_{2})}}=\frac{\frac{1}{n-1}\displaystyle\sum_{i=1}^{n_{2}}\sum_{j=1}^{n_{3}}\psi_{2}(X_{2i},X_{3j})}{\sqrt{\frac{n_{2}n_{3}}{{(n-1)^{2}}}\tau_{2}^{2}}}\xrightarrow{D}N(0,1)\hbox{
as n }\rightarrow\infty.$ (23)

Observe that $\sqrt{\frac{n_{2}n_{3}}{{(n-1)^{2}}}\tau_{2}^{2}}=O(1)$ and
$\tau_{n}=O(1)$. Then by Lemma 2.1 follows that

$\displaystyle\frac{W_{2}}{\tau_{n}}\xrightarrow{D}J_{2}\sim\operatorname{N}(0,M_{2})\hbox{,
where
}M_{2}=\lim_{n\rightarrow\infty}\left(\frac{\frac{n_{2}n_{3}}{{(n-1)^{2}}}\tau_{2}^{2}}{\tau^{2}_{n}}\right).$
(24)

Similarly,

$\displaystyle\frac{W_{3}-\operatorname{E}(W_{3})}{\sqrt{\operatorname{Var}(W_{3})}}=\frac{\frac{(2+n_{3})}{(n-1)(n_{2}-1)}\displaystyle\sum_{1\leq
i<j\leq
n_{2}}\psi_{2}(X_{2i},X_{2j})}{\sqrt{\frac{(2+n_{3})^{2}n_{2}\tau_{2}^{2}}{2(n-1)^{2}(n^{2}-1)}}}\xrightarrow{D}N(0,1)\hbox{
as n }\rightarrow\infty.$ (25)

Other properties are that
$\sqrt{\frac{(2+n_{3})^{2}n_{2}\tau_{2}^{2}}{2(n-1)^{2}(n^{2}-1)}}=O(1)$ and
$\tau_{n}=O(1)$, then by the Lemma 2.1

$\displaystyle\frac{W_{3}}{\tau_{n}}\xrightarrow{D}J_{3}\sim\operatorname{N}(0,M_{3})\hbox{,
where
}M_{3}=\lim_{n\rightarrow\infty}\left(\frac{\frac{(2+n_{3})^{2}n_{2}\tau_{2}^{2}}{2(n-1)^{2}(n^{2}-1)}}{\tau^{2}_{n}}\right).$
(26)

Analogously,

$\displaystyle\frac{W_{4}-\operatorname{E}(W_{4})}{\sqrt{\operatorname{Var}}(W_{4})}=\frac{\frac{(2+n_{2})}{(n-1)(n_{3}-1)}\displaystyle\sum_{1\leq
i<j\leq
n_{3}}\psi_{2}(X_{3i},X_{3j})}{\sqrt{\frac{(2+n_{2})^{2}n_{3}\tau_{2}^{2}}{2(n-1)^{2}(n_{3}-1)}}}\xrightarrow{D}N(0,1)\hbox{
as n }\rightarrow\infty.$ (27)

Once more,
$\sqrt{\frac{(2+n_{2})^{2}n_{3}\tau_{2}^{2}}{2(n-1)^{2}(n_{3}-1)}}=O(1)$ and
$\tau_{n}=O(1)$, then

$\displaystyle\frac{W_{4}}{\tau_{n}}\xrightarrow{D}J_{4}\sim\operatorname{N}(0,M_{4})\hbox{,
where
}M_{4}=\lim_{n\rightarrow\infty}\left(\frac{\frac{(2+n_{2})^{2}n_{3}\tau_{2}^{2}}{2(n-1)^{2}(n_{3}-1)}}{\tau^{2}_{n}}\right).$
(28)

Thus, applying Slutsky’s theorem we have

$\displaystyle\frac{(n/2)B_{n}}{\tau_{n}}$ $\displaystyle=$
$\displaystyle\frac{(n/2)B_{n}}{(n/2)V_{n}^{1/2}}=\frac{B_{n}}{\sqrt{\operatorname{Var}(B_{n})}}$
(29) $\displaystyle=$
$\displaystyle\frac{W_{1}+W_{2}-W_{3}-W_{4}}{\tau_{n}}\xrightarrow{D}J_{1}+J_{2}-J_{3}-J_{4}\hbox{
as n }\rightarrow\infty.$

This result shows that the test statistic asymptotically converges in $n$ to a
non-degenerate random variable whose limit distribution depends on the choice
of kernel $\phi(\cdot,\cdot)$.

###### Theorem 2

Let ${\bf X}_{1},{\bf X}_{2},\dots,{\bf X}_{n}$ be a sequence of i.i.d.
$L\times 1$ random vectors. Let $\phi(\cdot,\cdot)$ be a kernel of degree 2
such that

$\phi({\bf X}_{i},{\bf
X}_{j})=\frac{1}{L}\sum_{l=1}^{L}\phi^{*}(X_{il},X_{jl})$ (30)

for some kernel $\phi^{*}(\cdot,\cdot):\mathbb{R}^{2}\rightarrow\mathbb{R}$,
where $X_{il}$ is the $l$-th entry of ${\bf X}_{i}$. Define
$\phi_{1}^{*}(x_{il})=\mathbb{E}[\phi^{*}(X_{il},X_{jl})|X_{il}=x_{il}]$ and
suppose $\operatorname{Var}(\phi_{1}^{*}(X_{il}))>0$ and
$\operatorname{Var}(\phi^{*}(X_{il},X_{jl}))<\infty$. Let $B_{n}$ be defined
by (12) for the case where $n_{1}=1$, and assume that all conditions in
Theorem 1 hold. Suppose also that

$\sum_{1\leq l<m\leq
n}\mathbb{E}[\phi^{*}(X_{il},X_{jl})\phi^{*}(X_{im},X_{jm})]=O(L)$ (31)

and

$\sum_{1\leq l<m\leq
L}\mathbb{E}[\phi_{1}^{*}(X_{il})\phi_{1}^{*}(X_{jm})]=O(L).$ (32)

Then

$\frac{B_{n}}{\sqrt{\operatorname{Var}(B_{n})}}\xrightarrow{D}N(0,1)\quad\mbox{
as }\quad L\rightarrow\infty.$ (33)

Proof:  We start writing $\psi_{1}(X_{i})$ and $\psi_{2}(X_{i},X_{j})$ as a
function of $\phi^{*}_{1}(\cdot)$ and $\phi^{*}_{2}(\cdot,\cdot)$. Note that

$\displaystyle\psi_{1}\left(\mathbf{X}_{i}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{L}\sum_{l=1}^{L}\psi_{1}^{*}\left(X_{il}\right)$ (34)
$\displaystyle\psi_{2}\left(\mathbf{X}_{i},\mathbf{X}_{j}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{L}\sum_{l=1}^{L}\phi^{*}(X_{il},,X_{jl})-\frac{1}{L}\sum_{l=1}^{L}\psi^{*}_{1}(X_{il})$
(35) $\displaystyle-\frac{1}{L}\sum_{l=1}^{L}\psi^{*}_{1}(X_{jl})-\theta$

where

$\displaystyle\psi^{*}_{1}(X_{il})$ $\displaystyle=$
$\displaystyle\phi^{*}_{1}(X_{il})-\theta$ (36)
$\displaystyle\phi^{*}_{1}(x_{il})$ $\displaystyle=$
$\displaystyle\mathbb{E}[\phi^{*}(X_{il},X_{jl})\mid X_{il}=x_{il}]$ (37)
$\displaystyle\phi^{*}_{2}(x_{il},x_{jl})$ $\displaystyle=$
$\displaystyle\mathbb{E}[\phi^{*}(X_{il},X_{jl})\mid
X_{il}=x_{il},X_{jl}=x_{jl}].$ (38)

We can write $\psi_{1}(\cdot)$ as

$\displaystyle\psi_{1}(\boldmath{X}_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{L}\sum_{l=1}^{L}\left[\phi^{*}_{1}(X_{ij})-\theta\right],$
(39)

or

$\displaystyle\psi_{1}(\boldmath{X}_{i})$ $\displaystyle=$
$\displaystyle\frac{1}{L}\sum_{l=1}^{L}\psi_{1}^{*}(X_{ij}).$ (40)

Thus the variance of $\psi_{1}(\cdot)$ is given by

$\displaystyle\operatorname{Var}(\psi_{1}(\boldmath{X}_{i}))$ $\displaystyle=$
$\displaystyle\operatorname{Var}\left[\frac{1}{L}\sum_{l=1}^{L}\psi_{1}^{*}(X_{il})\right].$
(41)

By (31) we have that

$\displaystyle\operatorname{Var}\left(\psi_{1}\left(\mathbf{X}_{i}\right)\right)$
$\displaystyle=$
$\displaystyle\frac{1}{L^{2}}\left\\{\sum_{l=1}^{L}\operatorname{Var}\left[\psi_{1}^{*}\left(X_{il}\right)\right]\right.$
(42) $\displaystyle\left.+2\sum_{1\leq l<m\leq
L}\operatorname{Cov}\left(\psi_{1}^{*}\left(X_{il}\right),\psi_{1}^{*}\left(X_{im}\right)\right)\right\\}$
$\displaystyle=$ $\displaystyle O\left(L^{-1}\right)$

and by (32) the variance of $\psi_{2}(,)$ is

$\displaystyle\operatorname{Var}\left(\psi_{2}\left(\mathbf{X}_{i},\mathbf{X}_{j}\right)\right)=$
$\displaystyle\frac{1}{L^{2}}\left\\{\sum_{l=1}^{L}\operatorname{Var}\left(\phi^{*}\left(X_{il},X_{jl}\right)\right)+\right.$
(43) $\displaystyle+2\sum_{1\leq l<m\leq
L}\operatorname{Cov}\left(\phi^{*}\left(X_{il},X_{jl}\right),\phi^{*}\left(X_{im},X_{jm}\right)\right)$
$\displaystyle\left.+2\operatorname{Var}\left(\frac{1}{L}\sum_{l=1}^{L}\psi_{1}^{*}\left(X_{il}\right)\right)\right\\}$
$\displaystyle=$ $\displaystyle O\left(L^{-1}\right).$

Thus, for fixed $n$ and for $L\rightarrow\infty$ it follows that

$\displaystyle\frac{B_{n}}{\sqrt{\operatorname{Var}(B_{n})}}=V_{n}^{-1/2}B_{n}\xrightarrow{D}N(0,1).$
(44)

This result is fundamental to our inference procedure for clustering in the
HDLSS context.

###  2.3 Variance of $B_{n}$

In the $utest$ the estimation of $B_{n}$’s variance under $H_{0}$ plays an
essential role in hypothesis testing (see Cybis et al. (2018) ). As shown
below, even under $H_{0}$, the variance of $B_{n}$ depends on the particular
group configuration under consideration. For the homogeneity test of Section
3, we must evaluate this variance for the many group configurations visited in
an optimization algorithm. This variance estimation is performed through a
resampling procedure, however it becomes computationally expensive to perform
one resampling procedure for each individual group size configuration. To
circumvent this issue, Cybis et al. (2018) propose a reweighting scheme taking
advantage of analytic calculations for the variance for the case $K=2$ groups.
They are able to compute all variances from a single resampling procedure. In
this section we extend their argument to the case of $K=3$ groups.

In this Section we provide an estimator for the variance of $B_{n}$ under
$H_{0}$ based on U-statistics properties of $B_{n}$. For cases where all
groups have more than two elements, the Hoeffding decomposition of $B_{n}$ can
be found in Pinheiro et al. (2009) which is given by

$\displaystyle B_{n}$ $\displaystyle=$
$\displaystyle\left(\frac{2}{n(n-1)}\right)\sum_{1\leq i<j\leq
n}\eta_{nij}\psi_{2}(X_{i},X_{j}),$ (45)

where $\psi_{2}(\dot{,}\dot{)}$ is the second order term of the Hoeffding
decomposition of $B_{n}$ and

$\displaystyle\eta_{nij}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}1,&\hbox{ if }i\hbox{ and }j\hbox{
are from different groups}\\\ \\\ -\frac{(n-n_{g})}{n_{g}-1},&\hbox{ if
}i\hbox{ and }j\hbox{ are from the same group }g.\end{array}\right.$ (49)

Thereby,

$\displaystyle\operatorname{Var}(B_{n})=\left(\frac{2}{n(n-1)}\right)^{2}\tau_{2}^{2}\sum_{1\leq
i<j\leq n}\eta^{2}_{nij}.$ (50)

where $\tau_{2}^{2}=\operatorname{Var}(\psi_{2}(X_{1},X_{2}))$. From Pinheiro
et al. (2009) we also know that

$\displaystyle\sum_{1\leq i<j\leq n}\eta_{nij}^{2}={n\choose
2}(G-1)\left\\{1+\frac{1}{n}\sum_{g=1}^{G}\frac{n-n_{g}}{\left(n_{g}-1\right)(G-1)}\right\\}.$
(51)

For the case in which we have three groups, $G_{1}$, $G_{2}$ and $G_{3}$, with
sizes $n_{1}$, $n_{2}$ and $n_{3}$, respectively, where $n_{1}+n_{2}+n_{3}=n$,
it can be rewritten as

$\displaystyle C_{n}(n_{1},n_{2})=\sum_{1\leq i<j\leq
n}\eta_{nij}^{2}=2{n\choose
2}\left\\{1+\frac{1}{n}\sum_{g=1}^{3}\frac{n-n_{g}}{2\left(n_{g}-1\right)}\right\\},$
(52)

and therefore

$\displaystyle\operatorname{Var}(B_{n})$ $\displaystyle=$
$\displaystyle\left(\frac{2}{n(n-1)}\right)^{2}\tau_{2}^{2}C_{n}(n_{1},n_{2})=V_{n_{1},n_{2}}.$
(53)

Note that only $\tau_{2}^{2}$ depends on the probability distribution of the
data. Given three groups of sizes $n_{1}$, $n_{2}$ and $n_{3}$, the variance
of $B_{n}$ for this configuration is estimated through a resampling procedure.
For optimization purposes, it is not interesting to perform a resampling
procedure for each group configuration, so the idea is to use (the relation)
expression (53) to estimate $B_{n}$’s variance for any group configuration
from a single resampling procedure. Let $G_{1}^{*}$, $G_{2}^{*}$ and
$G_{3}^{*}$, with sizes $n_{1}^{*}$, $n_{2}^{*}$ and $n_{3}^{*}$,
respectively, where $n_{1}^{*}+n_{2}^{*}+n_{3}^{*}=n$, be an other group
configuration for the same data set. From (53) it follows that

$\displaystyle V_{n_{1}^{*},n_{2}^{*}}$ $\displaystyle=$
$\displaystyle\frac{C_{n}(n_{1}^{*},n_{2}^{*})}{C_{n}(n_{1},n_{2})}V_{n_{1},n_{2}}.$
(54)

Thus estimating $V_{n_{1},n_{2}}$ through a resampling procedure is sufficient
to estimate the variance of $B_{n}$ for any other group configuration.
Although the variance of $B_{n}$ is estimated under $H_{0}$, we note that the
choice of $n_{1}$ and $n_{2}$ may be important to reduce the bias of the
variance estimator. To understand the $C_{n}(\cdot,\cdot)$ function’s behavior
we plot (52) assuming that $n_{1},n_{2},n_{3}\geq 2$ and
$n=n_{1}+n_{2}+n_{3}$. As $\tau_{2}^{2}$ does not depend on group sizes, the
behavior of $C_{n}(\cdot,\cdot)$ governs the behavior of $B_{n}$’s variance
and Figure 1 shows that smaller values are obtained when groups have balanced
sizes, while larger values of $C_{n}(\cdot,\cdot)$ are obtained when group
sizes are unbalanced.

Figure 1:  $C_{n}(\cdot,\cdot)$ function behavior for $n_{1},n_{2},n_{3}\geq
2$ and $n=n_{1}+n_{2}+n_{3}$.

####  2.3.1 Variance of the extended $B_{n}$

We propose an extended statistic $B_{n}$ in (12) to accommodate cases in which
the data set is divided into three groups, one of which has size one. For
inference purposes it is essential establish a strategy to estimate the
variance of the extended $B_{n}$. Through the Hoeffding decomposition of (12)
(see Suplementary Material) we have that the variance of the extended $B_{n}$
is

$\displaystyle\operatorname{Var}(B_{n})=\zeta_{1}(n)\tau_{1}^{2}+\zeta_{2}(n,n_{2})\tau_{2}^{2},$
(55)

where $\tau_{1}^{2}=\operatorname{Var}(\psi_{1}(X_{1}))$ and
$\tau_{2}^{2}=\operatorname{Var}(\psi_{2}(X_{1},X_{2}))$ are, respectively,
the variance of the first and the second order terms of the Hoeffding
decomposition,

$\displaystyle\zeta_{1}(n)$ $\displaystyle=$ $\displaystyle\frac{4}{n(n-1)},$
$\displaystyle\zeta_{2}(n,n_{2})$ $\displaystyle=$
$\displaystyle\frac{4}{n^{2}(n-1)}+\frac{4n_{2}n_{3}}{n^{2}(n-1)^{2}}+\frac{2n_{2}(2+n_{3})^{2}}{n^{2}(n_{2}-1)(n-1)^{2}}$
$\displaystyle+\frac{2n_{3}(2+n_{2})^{2}}{n^{2}(n_{3}-1)(n-1)^{2}},$

$n_{1}=1$, and $n_{3}=n-n_{2}-1$. Note that in expression (55) the terms
$\tau_{1}^{2}$ and $\tau_{2}^{2}$ depend on the probability distribution of
the data, $\zeta_{1}(\cdot)$ depends only on $n$ and $\zeta_{2}(\cdot,\cdot)$
depends on $n$ and $n_{2}$ since $n_{3}=n-n_{2}-1$. Thus for another group
configuration keeping one of the groups with size one, the only change occurs
at $n_{2}$, say $n_{2}^{*}$. For this new group configuration, the extended
$B_{n}$ variance is given by

$\displaystyle\operatorname{Var}(B_{n})=\zeta_{1}(n)\tau_{1}^{2}+\zeta_{2}(n,n_{2}^{*})\tau_{2}^{2}.$
(57)

Again, the choice of $n_{2}$ may affect the variance of the estimator.
Denoting (55) by $V_{n_{2}}$ and (57) by $V_{n_{2}^{*}}$, we have from simple
algebra that

$\displaystyle V_{n_{2}^{*}}$ $\displaystyle=$ $\displaystyle
V_{n_{2}}+[\zeta_{2}(n,n_{2}^{*})-\zeta_{2}(n,n_{2})]\tau_{2}^{2}.$ (58)

For a given $n_{2}$ we can estimate $V_{n_{2}}$ from a resampling procedure.
Additionally, an estimate for $\tau_{2}^{2}$ can be obtained from the strategy
employed to estimate the variance of $B_{n}$ without outlier through
expression (53) as

$\displaystyle\widehat{\tau}_{2}^{2}$ $\displaystyle=$
$\displaystyle\frac{\widehat{V}_{n_{1},n_{2}}}{C(n_{1},n_{2})\left(\frac{2}{n(n-1)}\right)^{2}}.$
(59)

Thus we have a procedure to estimate the extended $B_{n}$´s variance for any
group configuration from only two independent resampling procedures, through
expression

$\displaystyle\widehat{V}_{n_{2}^{*}}$ $\displaystyle=$
$\displaystyle\widehat{V}_{n_{2}}+[\zeta_{2}(n,n_{2}^{*})-\zeta_{2}(n,n_{2})]\widehat{\tau}_{2}^{2},$
(60)

where $\widehat{\tau}_{2}^{2}$ is obtained from the resampling employed to
estimate the variance of $B_{n}$ without outlier and $\widehat{V}_{n_{2}}$ is
obtained from an additional resampling specific to $n_{1}=1$ case. Thus,
taking into account the resampling procedure performed to estimate the
variance of $B_{n}$ when the groups are larger than two and, with one more
resampling procedure for the size one group, we have an estimator for extended
$B_{n}$’s variance.

In Figure 2 we have the behavior of $\zeta_{2}(n,n_{2})$ as a function of
$n_{2}$.

Figure 2: Behavior of function $\zeta_{2}(n,n_{2})$ for a given $n$, with
$n_{1}=1$ and $n=1+n_{2}+n_{3}$.

These results are fundamental for the development of feasible algorithms that
find significant clusters which is computationally challenging problem.

##  3 Homogeneity test for three groups

Assessment of group homogeneity is a great challenge for standard statistics,
especially in the HDLSS context. The $uclust$ algorithm presented in Cybis et
al. (2018) and Valk and Cybis (2020) is effective to assess overall group
homogeneity by verifying whether there exists some significant partition of
the data in two groups. Here we are proposing an extension of the $uclust$
algorithm for data partitions in three groups $G_{1}$, $G_{2}$ and $G_{3}$. A
combinatorial procedure like the one proposed by Valk and Pinheiro (2012) in
which a $utest$ is applied for each possible partition of all group elements
into three subgroups has serious computational restrictions due to the
exponential increase in the number of tests that need to be performed.

###  3.1 Total of combinations

In order to develop the homogeneity test we require the number of different
group configurations that can be formed by separating $n$ elements,
$x_{1},x_{2},\dots,x_{n}$ into three groups, $G_{1},G_{2}$ e $G_{3}$. Follows
from Valk and Pinheiro (2012) that the number of combination of $n$ elements
into two groups is $p(n)=2^{n-1}-n-1$. Then if we divide $n$ elements into
three groups where one of them has size 1, it follows that the number of
combinations is

$\displaystyle\delta_{3}(n)$ $\displaystyle=$ $\displaystyle(2^{n-2}-n)n.$
(61)

Now we focus on the case where all groups have more than one element. We can
fix, without loss of generality, $x_{1}$ as an element that belongs to the
first group, $G_{1}$. Thus, we still have $n-1$ elements to be distributed
among the three groups. Since we cannot have a unitary group, we need at least
one more point for the first group. This group can have up to $n-4$
observations, since the remaining sets must necessarily have two elements
each. Thus, we then have the following number of possible first sets

$\displaystyle{n-1\choose 1}+{n-1\choose 2}+\cdots+{n-1\choose n-5}.$

For the remaining elements that need to be divided into two clusters, just
divide them into two groups with at least $2$ elements in each using the
function $p(\cdot)$. Combining these results, we have a number of different
configurations of non-unitary groups when we separate $n$ elements into $3$
groups given by

$\displaystyle S_{3}(n)$ $\displaystyle=$ $\displaystyle{n-1\choose
1}p(n-2)+{n-1\choose 2}p(n-3)+\cdots+{n-1\choose n-5}p(4)$ (62)
$\displaystyle=$ $\displaystyle\sum_{k=1}^{n-5}{n-1\choose k}p(n-k-1).$

We can still rewrite this equation on a recurring basis. Note that if we
already know how many configurations of groups we have with $n$ non-unitary
elements, and how many configurations with a unitary group, then it is
possible to calculate $S_{3}(n+1)$ as

$\displaystyle S_{3}(n+1)$ $\displaystyle=$ $\displaystyle
3S_{3}(n)+\delta_{3}(n).$ (63)

With such equations we can rewrite $S_{3}(n)$ as

$\displaystyle S_{3}(n)$ $\displaystyle=$
$\displaystyle\frac{233(3^{n-6})+1+n+n^{2}-(2+n)2^{n-1}}{2}.$ (64)

Thus, the number of different group configurations where at most one of them
has size one is given by

$\displaystyle\gamma_{3}(n)$ $\displaystyle=$
$\displaystyle\frac{233(3^{n-6})+1+n+n^{2}-(2+n)2^{n-1}}{2}+\delta_{3}(n)$
(65) $\displaystyle=$ $\displaystyle\frac{233(3^{n-6})+1+n-n^{2}-2^{n}}{2}.$

which becomes computationally onerous, especially for large sample size $n$.
To address this issue, we proceed similarly to Cybis et al. (2018) proposing
an optimization procedure to assess group homogeneity by finding the group
configuration $G_{1}$, $G_{2}$ and $G_{3}$ that maximizes the objective
function

$f(G_{1},G_{2},G_{3})=\frac{B_{n}}{\sqrt{\operatorname{Var}(B_{n})}}.$ (66)

By maximizing the standardized $B_{n}$ we must apply only one test. If this
three group partition is found significant, then there is at least one
subgroup that is significantly different from the others. However, if $H_{0}$
is not rejected for this partition, then all other three group partitions will
also be non-significant, and the whole data will be considered homogeneous.
While only the group configuration with maximum standardized $B_{n}$ is tested
we have to consider the distribution of $B_{n}$’s maximum under $H_{0}$.
Making the untrue, but useful, simplifying assumption that the $B_{n}$’s are
independent for different group configurations, the asymptotic cumulative
distribution function of the maximum standardized $B_{n}$ is given by

$F_{\mbox{max}}(x)=\mathbb{P}\left(\mbox{max}\left(\frac{B_{n}}{\sqrt{\operatorname{Var}(B_{n})}}\right)<x\right)=\Phi(x)^{n^{*}},$
(67)

where $n^{*}=\gamma_{3}(n)$, for $\gamma_{3}(n)$ defined in (65) and
$\Phi(\cdot)^{n^{*}}$ is the standard normal cumulative distribution function
at the power $n^{*}$. For $F_{\mbox{max}}(x)>1-\alpha$, we reject the null
hypothesis of overall group homogeneity with $\alpha$ significance level.

The number of tests increases rapidly, even for moderate sample size due to
the combinatorial nature of our approach. The maximum distribution in (67)
adequately accounts for multiple testing for reasonably small values of
$n^{*}$. However, this approach has some shortcomings since $n^{*}$ rapidly
increases. Proceeding similarly to Valk and Cybis (2020) and considering the
simplifying assumption that the $B_{n}$’s are independent, we use extreme
value theory and model it as Gumbel. However, the Gumbel approximation is only
valid for very large values of $n^{*}$. Thus, for small $n$ we employ the
standard max distribution of (67), and when $n^{*}\geq 2^{28}$ the Gumbel
distribution.

###  3.2 The clustering method $uclust3$

Our homogeneity test in the Section 3 is a method that finds the configuration
of _three_ subgroups that maximizes the standardized $B_{n}$. This is
appropriate for the context, since if the homogeneity test accepts the null
for this partition, then it would also be accepted for all other partitions.
However, the standardized $B_{n}$ might not be the best criteria to choose
between competing partitions when more than one significant group separation
exists. This issue is addressed in Cybis et al. (2018) and arises from the
fact that the variance of $B_{n}$ has different magnitudes depending on
subgroup sizes $n_{1}$ and $n_{2}$ (expression (52) dictates the relationship
between variances, which is shown in Figure 1). Consequently, this criteria
favours partitions with group sizes of smaller variance, namely
$n_{1},n_{2}\approx n/3$. We note that the magnitude of the variance is quite
different when we have a size one group, being much smaller in that case.
Again if we use the standardized $B_{n}$ statistic as a criterion, we will
have an effect of choosing groups of size one over the configurations of
groups that present greater variance according to the Figure 1.

Considering this issue, we proceed similarly to Valk and Cybis (2020) starting
by testing overall group homogeneity which is based on maximum of standardized
$B_{n}$. If the dataset is not homogeneous we adopt instead the maximum
$B_{n}$ as the criteria for finding the configuration that better divides the
sample into _three_ groups. Thus our significance clustering algorithm
$uclust3$ will find the partition with maximum $B_{n}$ among the universe of
all significant partitions in _three_ groups. This is sufficient to ensure
that the chosen configuration is statistically significant. However, it is not
efficient to find all arrangements of the data in _three_ groups that are
statistically significant. Furthermore, we cannot simply test the clusters
that maximizes $B_{n}$ since there are non-homogeneous samples for which this
maximal partition is not significant.

Based on these characteristics of the $B_{n}$ we propose a restricted search
algorithm, which is based on the behaviors of the $B_{n}$’s variances (see
Figure 1). It starts from the group configuration that maximizes $B_{n}$ and
if that partition is not significant, it searches for partitions whose
$B_{n}$’s variances are smaller than the previous one. This is suitable since
only for smaller variances, standardized $B_{n}$ can be significant. The
equation (54) is used to avoid a new resampling procedure to estimate the
$B_{n}$’s variance. As there is a difference in the magnitudes of the
$B_{n}$’s variances (see Figures 1 and 2) this algorithm treats separately the
cases when we have a group of size one and the cases with no outlier. The
detailed algorithm can be found in Section S3 of the supplementary materials.

##  4 Simulation Studies

In this section we present simulation studies in order to evaluate some
aspects of our proposed methodology. For that we simulate canonical data and
use the euclidean distance on our studies, but those are not mandatory for our
methods. As presented in Section 2.3, $B_{n}$’s variance has a behavior that
depends on the groups sizes. Moreover when we have a size one group, the order
of magnitude of the $B_{n}$’s variance is quite different when compared to
cases in which groups sizes are larger than one. For this reason, our
simulations studies typically have a configuration in which a group has size 1
and another configuration in which all groups have more than one element.
Figures 1 and 2 show that $B_{n}$’s variance is smaller at a central group
configuration, where the three groups have approximately the same number of
elements. Conversely, the variance is greater for extreme group
configurations, in which one of the groups has only two elements and the other
has $n/2$ elements (or $n-1-n_{2}$ elements for cases where we have a group of
size one). Naturally, the third group’s size is defined as
$n_{3}=n-n_{1}-n_{2}$. These scenarios are explored in our simulation studies.

In the Section 4.1 we evaluate the empirical size and power of the proposed
$utest$ for homogeneity of three groups. Section 4.2 present a simulation
study to evaluate the empirical properties of the homogeneity test $uclust3$.
The ability to find correct clusters of $uclust3$ and $kmeans$ clustering are
compared in Section 4.3.

###  4.1 Simulations for the $utest$

We present here a simulation study to evaluate the empirical performance of
the $utest$ for three groups. We simulate data from independent normally
distributed (i.i.d.) samples divided in three groups $G_{1}$, $G_{2}$ and
$G_{3}$. The elements of the $L$ dimensional vectors in $G_{1}$ are generated
from i.i.d. normal with mean $m_{1}=0$ and standard deviation equal to one.
The vectors in $G_{2}$ and $G_{3}$ have the same properties with mean $m_{2}$
and $m_{3}$, respectively. In order to allow a graphical representation of the
power of the test which is the proportion of rejection considering a
significance level $\alpha$ (the power curves), the groups were symmetrically
separated and on the x-axis the difference $m_{2}-m_{1}$ is reported. The
difference $m_{3}-m_{2}=m_{2}-m_{1}$. The sample size $n$ takes values in
$\\{10,20,50\\}$. Figure 3 presents power curves of the $utest$ for _three_
groups with separation degree $m_{2}-m_{1}$, where the vectors have dimension
$L=1000$ (gray) and $L=2000$ (black) and we have 100 replications of each
scenario. Furthermore group $G_{1}$ has size one and group $G_{2}$ was set to
have size $n_{2}=\lfloor n/3\rfloor$, where $\lfloor x\rfloor$ means the
integer part of $x$. Naturally the third group’s size is defined as
$n_{3}=n-1-n_{2}$. The significance level used to determine whether the test
rejects the null hypothesis that the elements in $G_{1}$, $G_{2}$ and $G_{3}$
have the same distribution was $\alpha=0.05$.

Figure 3: Power curves of $utest$ for two dimension $L=1000$ (gray) and
$L=2000$ (black) for 100 replications of each scenario of $n\in\\{10,20,50\\}$
with $\alpha=0.05$.

The empirical results obtained in this study reported in Figure 3 corroborate
the theoretical properties. As the $L$ increases, the rejection ratio also
increases and as the groups become more separated, the power increases. When
there is no separation, $m_{2}-m_{1}=0$, the rejection ratio is close to the
significance level $\alpha$ suggesting control of Type I error. Similar
results are found for cases where all groups have more than one element (see
Figure S1 in the Supplementary Material).

###  4.2 Simulations for the homogeneity test in $uclust3$

To evaluate the statistical properties of the homogeneity test $uclust3$
considering the max distribution (67) with the Gumbel correction when
appropriate, we simulate data with the same characteristics as the data in
Section 4.1. For each sample size $n$ in $\\{10,20,50\\}$, group $G_{1}$ has
size one and group $G_{2}$ was set to have size $n_{2}=2$ and $n_{2}=n/2$, and
consequently the third group’s size was defined as $n_{3}=n-1-n_{2}$. Table 1
shows the proportion of rejection of the null hypothesis for significance
level $\alpha=0.05$ considering two scenarios of $(m_{2},m_{3})$ and the
dimension $L$ taking values in $\\{1000,2000\\}$.

Table 1: Empirical power of the homogeneity test $uclust3$ with a group of size one | $(m_{2},\,\,m_{3})$ | ($n_{2}$) | Dimension $L$
---|---|---|---
$n$ | 1000 | 2000
10 | (0.25, 0.5) | 2 | 0.27 | 0.36
5 | 0.69 | 0.89
(0.5, 1) | 2 | 0.22 | 0.25
5 | 0.98 | 1
20 | (0.25, 0.5) | 2 | 0.93 | 1
10 | 1 | 1
(0.5, 1) | 2 | 0.9 | 0.89
10 | 0.92 | 1
50 | (0.25, 0.5) | 2 | 0.68 | 0.68
25 | 1 | 0.99
(0.5, 1) | 2 | 0.99 | 0.96
25 | 1 | 1

We can observe that even in an extreme group configuration, where the group
$G_{1}$ has size one and the group $G_{2}$ has size two, the method presents
consistent empirical power to reject the null hypothesis. The power increases
as $L$ and/or $n$ and/or the difference between $m_{2}$ and $m_{3}$ increases,
emphasizing the inherent properties of the method.

Supplementary Table S1 presents estimates of type I error rates for $uclust3$.
The significance level considered in this simulations was $\alpha=0.05$ and we
can observe that the method presents an adequate control of the Type I Error
for cases where $L>>n$ (typically HDLSS scenario). Supplementary Table S2
presents power of the $uclust3$ for group configurations of sizes greater than
1. For small sample size $n$ the test had more difficulty in finding the
correct clusters. However, for larger $n$ the method showed an excellent
performance.

###  4.3 Simulations for finding correct clusters

In order to evaluate the accuracy of our clustering method, we present
simulation studies comparing $uclust3$ with $kmeans$ clustering, one of the
most popular clustering algorithms. We refer the reader to the vastly cited
work of Jain (2010) for a general discussion about $kmeans$. The data were
simulated under the same distribution scheme of Section 4.2, with $Re=100$
replications and the methods were compared in terms of mean Adjusted Rand
Index (ARI) which measures the agreement of clustering results with simulation
scenarios, adjusting for randomness Hubert and Arabie (1985). An ARI of one
indicates perfect matching. No inference is used in this analysis. This is an
appropriate comparison as both methods are set to find exactly three groups.
Table 2 reports the results for three sample sizes $n\in\\{10,20,50\\}$, two
dimension $L\in\\{1000,2000\\}$ and three groups of sizes $n_{1}$, $n_{2}$ and
$n_{3}=n-n_{1}-n_{2}$. The data vectors in group $G_{1}$ have zero mean and
the data vectors in $G_{2}$ and $G_{3}$ have mean $m_{2}$ and $m_{3}$,
respectively. Note that the clustering method $uclust3$, based on the
maximization of $B_{n}$ is comparable to $kmeans$ to find the correct
clusters, considering this data configuration. However for larger sample
sizes, as the clusters become better defined, with greater separation between
the means, $uclust3$ outperforms $kmeans$. Table S3 shows that for the case
where $G_{1}$ has size one, $kmeans$ tends to perform slightly better for
smaller sample sizes.

Table 2: Comparison of mean ARI and standard deviation (Sd) of the accuracy in clustering of $kmeans$ and $uclust3$ methods. | $(m_{2},\,\,m_{3})$ | ($n_{1},\,\,n_{2}$) | Method | Dimension $L$
---|---|---|---|---
$n$ | 1000 | 2000
| Mean | Sd | Mean | Sd
10 | (0.25, 0.5) | (2, 5) | $kmeans$ | 0.59 | 0.05 | 0.73 | 0.06
$uclust3$ | 0.58 | 0.03 | 0.63 | 0.02
(3, 3) | $kmeans$ | 0.56 | 0.05 | 0.74 | 0.08
$uclust3$ | 0.52 | 0.05 | 0.6 | 0.05
(0.5, 1) | (2, 5) | $kmeans$ | 0.91 | 0.04 | 0.94 | 0.03
uclust3 | 0.74 | 0.01 | 0.74 | 0
(3, 3) | $kmeans$ | 0.9 | 0.05 | 0.87 | 0.07
$uclust3$ | 0.92 | 0.03 | 0.96 | 0.02
20 | (0.25, 0.5) | (2, 10) | $kmeans$ | 0.73 | 0.02 | 0.77 | 0.03
$uclust3$ | 0.7 | 0.02 | 0.74 | 0.02
(6, 6) | $kmeans$ | 0.74 | 0.05 | 0.94 | 0.03
$uclust3$ | 0.68 | 0.04 | 0.91 | 0.02
(0.5, 1) | (2, 10) | $kmeans$ | 0.96 | 0.01 | 0.94 | 0.02
$uclust3$ | 1 | 0 | 1 | 0
(6, 6) | $kmeans$ | 0.81 | 0.07 | 0.84 | 0.07
$uclust3$ | 1 | 0 | 1 | 0
50 | (0.25, 0.5) | (2, 25) | $kmeans$ | 0.76 | 0.01 | 0.79 | 0.01
$uclust3$ | 0.73 | 0 | 0.74 | 0.01
(16, 16) | $kmeans$ | 0.93 | 0.02 | 0.89 | 0.05
$uclust3$ | 0.94 | 0 | 1 | 0
(0.5 , 1) | (2 , 25) | $kmeans$ | 0.95 | 0.01 | 0.95 | 0.01
$uclust3$ | 1 | 0 | 1 | 0
(16, 16) | $kmeans$ | 0.8 | 0.07 | 0.81 | 0.07
$uclust3$ | 1 | 0 | 1 | 0

###  4.4 Finding correct clusters and comparing $uclust3$ and $uhclust$ in a
presence of an outlier

A simulation study similar to Section 4.1 was performed to compare our
$uclust3$ with the hierarchical methods $uhclust$ from Valk and Cybis (2020)
and $sigclust$ from Kimes et al. (2017); Kimes (2019) in terms of the ability
to correctly find statistically significant groups. The group $G_{1}$ has only
one element, the size of $G_{2}$ is $n_{2}=\lfloor n/3\rfloor$. For all three
methods the same level of significance $\alpha=0.05$ was considered. The
$sigclust$ method was not able to find the correct groups in any scenario,
with a proportion of correct answers equal to zero and for this reason it was
excluded from the analysis. Figures 4 and 5 report curves of proportion times
that the algorithms found significant separation and correct groups
considering different values of $m_{2}-m_{1}$ varying on the $x$ axis, with
sample size $n$ taking values in $\\{10,20,50\\}$ and dimension $L=1000$ and
$L=2000$ The results are based on 50 repetitions.

Figure 4: True cluster proportion curves of $uclust3$ (dark gray) and
$uhclust$ (light gray) for dimension $L=1000$ with 50 replications of each
scenario of $n$ with $\alpha=0.05$ and one outlier.  Figure 5: True cluster
proportion curves of $uclust3$ (dark gray) and $uhclust$ (light gray) for
dimension $L=2000$ with 50 replications of each scenario of $n$ with
$\alpha=0.05$ and one outlier.

The $uclust3$ method (dark grey) outperforms $uhclust$ method (light gray) in
all scenarios presenting greater ability to find the correct groups for less
separation. However, for $n=50$ these method are more competitive although the
method proposed here $uclust3$ still stands out for larger separations. The
conclusions do not change with the variation of dimension $L$. In Section S5
on the supplementary materials we present results of a simulation study for
the cases where there are no outlier. Supplementary Figures S2 and S3 shows
the true cluster proportion curves of $uclust3$ and $uhclust$ for dimension
$L=1000$ and $L=2000$. We note that the $uclust3$ method outperforms $uhclust$
in all scenarios.

##  5 Applications

###  5.1 Peripheral blood mononuclear cells

In order to illustrate of the applicability of the _utest_ we consider a one-
way MANOVA (multivariate analysis of variance) testing problem for high-
dimensional data. This issue was addressed in Zhang et al. (2017) by exploring
peripheral blood mononuclear cell (PBMC) data, consisting of 42 normal, 26
ulcerative colitis (UC) and 59 Crohn’s disease (CD) tissue samples ($n=127$),
each having $L=22,283$ gene expression level measurements. This dataset has
been studied by Burczynski et al. (2006) and is available at
_http://www.ncbi.nlm.nih.gov/gds_ with accession ID GDS1615. The classical
hypothesis test where the interest is to test whether the 3 mean vectors are
equal, can be described as follows: Let
$\mathbf{X}^{(g)}_{1},\ldots,\mathbf{X}^{(g)}_{n_{g}}$ be a sample of i.i.d.
vectors from the $L$-variate distribution ${F}_{g}$, with
$\mathbb{E}(\mathbf{X}^{(g)}_{1})=\bm{\mu}_{g}$ and
$\operatorname{cov}(\mathbf{X}^{(g)}_{1})=\bm{\Sigma},\quad$ for
$g=1,\ldots,3$ and $n=n_{1}+n_{2}+n_{3}$. Then, the null hypothesis is
$H_{0}:\bm{\mu}_{1}=\bm{\mu}_{2}=\bm{\mu}_{3}$. In our context, however, the
normality and variance homogeneity requirements are not necessary, and the
null hypothesis becomes the more general

$H_{0}:{F}_{1}={F}_{2}={F}_{3}.$

We apply the $utest$ for testing the equality of mean expression levels of the
normal, UC and CD groups of the PBMC data. The value of standardized $B_{n}$
statistic is 13.20997 (p-value$<<$0.001) with which we reject the null
hypothesis of equality of mean expression levels.

###  5.2 Image recognition

We consider a simple example of image recognition to illustrate the
applicability of our methodology. The data consists of images from three
public figures (Tony Blair, Colin Powell and George W. Bush) which were
selected from the Labeled Faces Wild (LFW) dataset (Huang et al. (2007)). The
data were run through OpenFace’s convolutional neural network (Amos et al.
(2016)), a procedure that outputs a 128-dimensional representation of the
faces which preserves Euclidean distances. In case the reader wants to know
more about how the OpenFace works, we recommend reading their website Amos et
al. (2016). In this illustrative application, we randomly select 10 images
from each public figure in the above cited dataset, run $uhclust$, $sigclust$
and $uclust3$ with significance level $\alpha=0.05$. Figure 6 presents the
hierarchical clustering dendrogram annotated with p-values for all tests
performed in the $uhclust$ method. We found 4 homogeneous groups, with a
significant division in the Bush image group and an ARI=0.8585. Figure 7
presents the dendrogram with corresponding $sigclust$ analysis of the same
data which produces six significant clusters, segregating Bush and Powell’s
images from the reminder and finding one outlier in Blair’s group. The ARI for
this case was 0.7788. Applying the $uhclust3$ method we found exactly 3
homogeneous groups, each corresponding to one of the public figures with
ARI=1.

In the Section S6 in the supplementary materials we consider the same dataset
and public figures to carry out an analysis with three groups in which one has
size one. Figures S4 and S5 in the supplementary materials present the
clustering dendrogram annotated with results of all tests performed in the
$uhclust$ and $sigclust$ methods. None of these methods were able to identify
the outlier and both methods achieved ARI of 0.8135593. However, when we
applied the $uclust3$ method we found the correct groups with ARI of 1,
supporting the best results $uclust3$ in the simulation study.

Figure 6: Annotated dendrogram of significance analysis for hierarchical
clustering $uhclust$ for 30 pictures of 3 public figures. P-values and
corrected significance levels $\alpha^{*}$ are shown for each test performed
at the corresponding node.  Figure 7: Annotated dendrogram of significance
analysis for hierarchical clustering $sigclust$ for 30 pictures of 3 public
figures. P-values and corrected significance levels $\alpha^{*}$ are shown for
each test performed at the corresponding node.

##  6 Discussion

We have developed a clustering method that separates a dataset specifically
into _three_ groups allowing the assessment of significance of this partition.
Our methodology is based on the U-statistics clustering framework proposed in
Pinheiro et al. (2009) and is an extension of the approach of Cybis et al.
(2018); Valk and Cybis (2020). Considering the $B_{n}$ statistic of Pinheiro
et al. (2009) that aims to test homogeneity of three predefined groups we
propose an extension of the $B_{n}$ statistic to allow for an outlier, namely
one of the groups has only one element ($n_{1}=1$). Additionally we verified
statistical properties that ensure the compatibility of this new definition
with the overall framework. We then considered group homogeneity testing with
this newly defined statistic, and explored empirical properties such as Type I
error control and power, showing adequate preformance. Afterwards, we extended
this framework to address the issue of partitioning a dataset into three
optimal statistically significant clusters, proposing a new clustering
criteria that defines the $uclust3$ method. This differs from previous methods
for instead of find and testing a two group separation, $uclust3$ finds the
best significant partitions in three clusters. This can pave the way for
inference in $K$ groups.

This U-statistics based methodology can be applied to a wide range of
problems, since they make very few assumptions about the distribution of the
data. Although in the simulation study and in the application we have used
Euclidean distance, this is not a necessary requirement for theory
development. Additionally, even if the data come from a non-normal
multivariate distribution, the required asymptotic normality is guaranteed as
long as the distances have finite variance and the sum of all distance
covariances do not grow too fast ($O(L)$ see Theorem 2). The clustering
procedures $uclust3$ proposed here require large $L$ since $B_{n}$ for
$n_{1}=1$ is only asymptotically normal in the dimension $L$. As verified in
previously work of Valk and Cybis (2020), for the settings in the simulation
studies, in practice our tests achieve good Type I error control having
difficulties only when $L$ is smaller than $10n$. This is, by excellence, the
HDLSS setting.

An important step for developing the homogeneity test is to establish the
number of possible configurations of $n$ elements separated in three groups. A
system of recursive equations was developed to solve this combinatorial
problem and the idea may be used to solve an equivalent problem involving
$K>3$ groups.

The significance clustering method $uclust3$ proposed here returns the
partition that better separates the data into three statistically significant
groups in terms of the $B_{n}$ statistic. Thus we can compare it with
$kmeans$, which is one of the most popular clustering method, regarding the
ability of correctly find three groups. A simulation study suggests that
$uclust3$ is competitive with $kmeans$ when we have a size one group and
outperforms $kmeans$ in the context in which groups having an underlying
cluster structure with more than 2 elements each and large sample sizes.

Since our methodology is a natural extension of the $uclust$ method proposed
by Valk and Cybis (2020) it inherits many helpful properties such as the
ability to avoid the hazards of directly estimating the covariance matrix, by
obtaining $\operatorname{Var}(B_{n})$ through resampling. However, they have
different purposes, while $uclust$ aims to find the best significant partition
in two groups, $uclust3$ aims to find the best significant separation in three
groups, so they are not directly comparable. To support the usefulness of the
$uclust3$, we carried out a simulation study to compare this method with the
hierarchical version of $uclust$ ($uhclsut$) and with another hierarchical
approach ($sigclust$), which both are able to find a significant partition
into three groups, when this partition exists. We simulated normal data with a
three group structure, separating these groups in terms of the means and use
the proportion of correct configurations found to compare the methods. In the
situations considered, $sigclust$ had serious difficulties in finding the
proper arrangement, while $uclust3$ performed better than $uhclust$ in all
scenarios. Additionally, in the applications we have shown the applicability
of this methodology, first with a one-way MANOVA testing problem without the
requirement of normality of data and variance homogeneity, and then with an
application to image recognition data where we select three public figures and
observe that the $uclust3$ method was the only one able to correctly find the
three groups of figures.

Finally the conclusion is that our $uclust3$ method is appropriate to separate
a high dimensional low sample size datasets into _three_ groups, being more
powerful than some other methods in the specific situation in which a
structure of three groups is present in the data.

## Supplementary material

Supplementary material:

Derivations, supplementary tables and figures (pdf)

Code:

R-functions containing all methods developed in this article (will be
available in the uclust package at CRAN).

Data:

Dataset used in the application and corresponding script (zip).

##  Acknowledgements

We would like to thanks Pedro Fusieger for the fruitful discussions about the
number of possible assignments of all $n$ elements in three subgroups.

##  Funding

Debora Zava Bello research was supported by Coordenação de Aperfeiçoamento de
Pessoal de Nível Superior (CAPES).

## References

  * Euan et al. [2019] Carolina Euan, Ying Sun, Hernando Ombao, et al.  Coherence-based time series clustering for statistical inference and visualization of brain connectivity.  _Annals of Applied Statistics_ , 13(2):990–1015, 2019.
  * Rosenberg et al. [2002] Noah A. Rosenberg, Jonathan K. Pritchard, James L. Weber, Howard M. Cann, Kenneth K. Kidd, Lev A. Zhivotovsky, and Marcus W. Feldman.  Genetic structure of human populations.  _Science_ , 298(5602):2381–2385, 2002.
  * Chen et al. [2015] Gary K. Chen, Eric C. Chi, John Michael O. Ranola, and Kenneth Lange.  Convex clustering: an attractive alternative to herarchical clustering.  _PLoS Computational Biology_ , 11(5):e1004228, 2015.
  * Motlagh et al. [2019] Omid Motlagh, Adam Berry, and Lachlan O’Neil.  Clustering of residential electricity customers using load time series.  _Applied energy_ , 237:11–24, 2019.
  * Hennig [2015] Christian Hennig.  What are the true clusters?  _Pattern Recognition Letters_ , 64:53–62, 2015.
  * Von Luxburg et al. [2012] Ulrike Von Luxburg, Robert C Williamson, and Isabelle Guyon.  Clustering: Science or art?  In _Proceedings of ICML workshop on unsupervised and transfer learning_ , pages 65–79. JMLR Workshop and Conference Proceedings, 2012.
  * Adolfsson et al. [2019] Andreas Adolfsson, Margareta Ackerman, and Naomi C Brownstein.  To cluster, or not to cluster: An analysis of clusterability methods.  _Pattern Recognition_ , 88:13–26, 2019.
  * McLachlan and Peel [2004] Geoffrey McLachlan and David Peel.  _Finite mixture models_.  John Wiley & Sons, 2004.
  * Demidenko [2018] Eugene Demidenko.  The next-generation k-means algorithm.  _Statistical Analysis and Data Mining: The ASA Data Science Journal_ , 11(4):153–166, 2018.
  * McShane et al. [2002] Lisa M McShane, Michael D Radmacher, Boris Freidlin, Ren Yu, Ming-Chung Li, and Richard Simon.  Methods for assessing reproducibility of clustering patterns observed in analyses of microarray data.  _Bioinformatics_ , 18(11):1462–1469, 2002.
  * Helgeson et al. [2020] Erika S Helgeson, David M Vock, and Eric Bair.  Nonparametric cluster significance testing with reference to a unimodal null distribution.  _Biometrics_ , 2020.
  * Shimodaira et al. [2004] Hidetoshi Shimodaira et al.  Approximately unbiased tests of regions using multistep-multiscale bootstrap resampling.  _The Annals of Statistics_ , 32(6):2616–2641, 2004.
  * Suzuki and Shimodaira [2006] Ryota Suzuki and Hidetoshi Shimodaira.  Pvclust: an R package for assessing the uncertainty in hierarchical clustering.  _Bioinformatics_ , 22(12):1540–1542, 2006.
  * Liu et al. [2008] Yufeng Liu, David Neil Hayes, Andrew Nobel, and JS Marron.  Statistical significance of clustering for high-dimension, low–sample size data.  _Journal of the American Statistical Association_ , 103(483):1281–1293, 2008.
  * Kimes et al. [2017] Patrick K Kimes, Yufeng Liu, David Neil Hayes, and James Stephen Marron.  Statistical significance for hierarchical clustering.  _Biometrics_ , 73(3):811–821, 2017.
  * Cybis et al. [2018] Gabriela B. Cybis, Marcio Valk, and Sílvia R. C. Lopes.  Clustering and classification problems in genetics through u-statistics.  _Journal of Statistical Computation and Simulation_ , 2018.
  * Valk and Cybis [2020] Marcio Valk and Gabriela Bettella Cybis.  U-statistical inference for hierarchical clustering.  _Journal of Computational and Graphical Statistics_ , 2020.
  * Pinheiro et al. [2009] Aluísio Pinheiro, Pranab Kumar Sen, and Hildete Prisco Pinheiro.  Decomposability of high-dimensional diversity measures: Quasi-u-statistics, martingales and nonstandard asymptotics.  _Journal of Multivariate Analysis_ , 2009.
  * Sen [2006] Pranab Kumar Sen.  Robust statistical inference for high-dimensional data models with application to genomics.  _Austrian journal of statistics_ , 35(2&3):197–214, 2006.
  * Hoeffding [1948] Wassily Hoeffding.  A class of statistics with asymptotically normal distribution.  _The Annals of Mathematical Statistics_ , pages 293–325, 1948.
  * Valk and Pinheiro [2012] Marcio Valk and Aluísio Pinheiro.  Time-series clustering via quasi u-statistics.  _Journal of Time Series Analysis_ , 2012.
  * Jain [2010] Anil K Jain.  Data clustering: 50 years beyond k-means.  _Pattern recognition letters_ , 31(8):651–666, 2010.
  * Hubert and Arabie [1985] Lawrence Hubert and Phipps Arabie.  Comparing partitions.  _Journal of classification_ , 2(1):193–218, 1985\.
  * Kimes [2019] Patrick Kimes.  pkimes/sigclust2 documentation, 2019.  URL https://rdrr.io/github/pkimes/sigclust2/man/.  (Accessed on 02/02/2021).
  * Zhang et al. [2017] Jin-Ting Zhang, Jia Guo, and Bu Zhou.  Linear hypothesis testing in high-dimensional one-way manova.  _Journal of Multivariate Analysis_ , 155:200–216, 2017\.
  * Burczynski et al. [2006] Michael E Burczynski, Ron L Peterson, Natalie C Twine, Krystyna A Zuberek, Brendan J Brodeur, Lori Casciotti, Vasu Maganti, Padma S Reddy, Andrew Strahs, Fred Immermann, et al.  Molecular classification of crohn’s disease and ulcerative colitis patients using transcriptional profiles in peripheral blood mononuclear cells.  _The journal of molecular diagnostics_ , 8(1):51–61, 2006.
  * Huang et al. [2007] Gary B. Huang, Manu Ramesh, Tamara Berg, and Erik Learned-Miller.  Labeled faces in the wild: A database for studying face recognition in unconstrained environments.  Technical Report 07-49, University of Massachusetts, Amherst, October 2007\.
  * Amos et al. [2016] Brandon Amos, Bartosz Ludwiczuk, and Mahadev Satyanarayanan.  OpenFace: A general-purpose face recognition library with mobile applications.  Technical report, CMU-CS-16-118, CMU School of Computer Science, 2016\.

See pages - of images/Appendix.pdf