# Robustly fitting Gaussian graphical models:
the R-package robFitConGraph

Daniel Vogel, Stuart J. Watt, and Anna Wiedemann Corresponding author:
<EMAIL_ADDRESS>
(December 1, 2022)

###### Abstract

This paper gives a tutorial-style introduction to the R-package
robFitConGraph, which provides a robust goodness-of-fit test for Gaussian
graphical models. Its use is demonstrated at a data example on music
performance anxiety, which also illustrates _why_ one would want to fit a
Gaussian graphical model – and why one should do so robustly. The underlying
statistical theory is briefly explained.

The paper describes package version 0.4.1, available on CRAN from December 1,
2022. See https://CRAN.R-project.org/package=robFitConGraph.

Keywords: covariance selection model, deviance test, M-estimator, music
performance anxiety, partial correlation

##  1 Introduction

The first two Sections 1 _Introduction_ and 2 _A Case Study_ are intended for
a general audience, assuming neither deeper familiarity with graphical models
nor robustness. Section 3 _Background and Theory_ discusses some aspects in
detail.

###  1.1 Gaussian graphical modeling

Graphical models are an important tool for analyzing the dependence structure
of several variables. Gaussian graphical models are employed for continuously
distributed data, where a multivariate normal, or Gaussian, distribution is
adequate.

A graph encodes the dependence structure of a random vector
$\boldsymbol{X}=(X_{1},\ldots,X_{p})$ as follows: The nodes represent the
individual variables, and an edge between two nodes represents partial
correlation, or the absence of an edge zero partial correlation. Partial
correlation is a measure for conditional dependence between two variables,
conditional on all remaining variables, which, in a very loose sense, can be
understood as a measure for dependence _not_ explained by joint dependence on
other variables. As soon as more than two variables are considered,
conditional dependence is arguably more important than marginal dependence: It
suggests and helps to verify causal relationships. While ice cream sales and
criminal assault rates per day – at any given place in the temperate climate
zone – are certainly positively correlated, their conditional independence
given a suitable mediator variable, such as outside mean temperature, provides
evidence for why this may be so.

Graphical modeling refers to the statistical task of finding an appropriate
graph that describes the dependence structure of a given data set, i. e.,
identifying all zero partial correlations. The aim is to find a parsimonious
graph, i. e., one with few edges, which does not contradict the data. A full
graph, with all edges present, means no restriction and contains no structural
information. A completely empty graph, with no edges at all, means all
variables are independent.

Let $\Sigma$ denote the covariance matrix of $\boldsymbol{X}$. Throughout, we
assume that all second moments of $\boldsymbol{X}$ are finite, and that
furthermore $\Sigma$ is positive definite and thus may be inverted to yield
the _concentration matrix_ or _precision matrix_ $K=\Sigma^{-1}$. The
assumption of strict positive definiteness is a mild one and is equivalent to
the probability mass not being concentrated on a lower-dimensional affine
linear subspace of $\mathbb{R}^{p}$. We make the same assumption on the data
points $\mathbb{X}_{n}=(\boldsymbol{x}_{1},\ldots\boldsymbol{x}_{n})^{\top}$,
and hence the sample covariance matrix $\hat{\Sigma}_{n}$ is positive definite
and the sample concentration matrix $\hat{K}_{n}=\hat{\Sigma}_{n}^{-1}$
exists.

The basis for Gaussian graphical modeling is the following characterization:
The variables $X_{i}$ and $X_{j}$ are partially uncorrelated given all
remaining $p-2$ variables if and only if $K_{i,j}=0$. The _partial
correlation_ $p_{i,j}$ of $X_{i}$ and $X_{j}$ given the remaining components
of $\boldsymbol{X}$ is defined as the correlation between the residuals of
$X_{i}$ and $X_{j}$ when regressing both on the remaining components. Some
matrix calculus yields that

$p_{i,j}=-\frac{K_{i,j}}{\sqrt{K_{i,i}K_{j,j}}},$ (1)

see, e. g., Whittaker, (1990, Chapter 5). Thus the pairwise partial
correlations are obtained from the inverse covariance matrix $K$ in a very
similar fashion as the pairwise correlations are obtained from the covariance
matrix $\Sigma$ itself. The only difference is the minus sign. Consequently,
an absent edge in the graph means a zero entry in $K$, and finding the graph
for given data comes down to finding the zero-pattern in the inverse of the
true covariance matrix $\Sigma$. Three basic sub-tasks of graphical modeling
can be identified:

  1. (T1)

Finding an appropriate graph.

  2. (T2)

Determining if a given graph fits the data.

  3. (T3)

Estimating the (remaining) partial correlations under a given graph structure.

Whether (T1) or (T2) is considered more important may be debatable and mainly
depends on whether one pursues an explorative or an inferential analysis. Task
(T3) may appear of lesser importance, but it is intrinsically linked to (T2).
We briefly outline how these tasks are approached, in reverse order, starting
in with (T3).

#### Task (T3)

We require some mathematical notation. Define a graph $G=(V,E)$ as a set of
vertices $V=\\{1,\ldots,p\\}$ and a set of undirected edges
$E\subseteq\\{\,\\{i,j\\}:i,j=1,\ldots,p,\ i<j\\}$. Let $\mathscr{S}_{p}$
denote the set of all symmetric $p\times p$ matrices and $\mathscr{S}^{+}_{p}$
the set of all positive definite, symmetric $p\times p$ matrices. For any
graph $G=(V,E)$, let $\mathscr{S}^{+}_{p}(G)$ be the set of matrices
$A\in\mathscr{S}^{+}_{p}$ with zero entries at off-diagonal positions
specified by $G$, i. e., $A_{i,j}=0$ for all $i,j=1,\ldots,p$, $i\neq j$, with
$\\{i,j\\}\notin E$. We call any set of $p$-dimensional probability measures
with the common property that they possess a concentration matrix
$K\in\mathscr{S}^{+}_{p}(G)$ a _covariance selection model_ induced by $G$. We
call a covariance selection model consisting of all regular $p$-variate
Gaussian distributions a _Gaussian graphical model_ and denote it by
$N_{p}(G)$, i. e.,
$N_{p}(G)=\\{N_{p}(\boldsymbol{\mu},\Sigma):\boldsymbol{\mu}\in\mathbb{R}^{p},\Sigma^{-1}\in\mathscr{S}^{+}_{p}(G)\\}$.
The maximum-likelihood estimator $\hat{\Sigma}_{G}$ of $\Sigma$ within the
parametric family $N_{p}(G)$ is given as the solution to

$\hat{\Sigma}_{G}=\operatorname*{arg\,min}_{\Sigma^{-1}\in\mathscr{S}^{+}_{p}(G)}\left\\{\log\det\Sigma+\frac{1}{n}\sum_{i=1}^{n}\operatorname*{trace}\left(\hat{\Sigma}_{n}\,\Sigma^{-1}\right)\right\\}.$
(2)

This optimization problem leads to the estimation equations

$\begin{cases}\
\left[\hat{\Sigma}_{G}\right]_{i,j}=\left[\hat{\Sigma}_{n}\right]_{i,j}&\qquad\mbox{
for }\\{i,j\\}\in E\ \mbox{ or }\ i=j,\\\\[4.30554pt] \
\left[\hat{\Sigma}_{G}^{-1}\right]_{i,j}=0&\qquad\mbox{ for }\\{i,j\\}\notin
E\ \mbox{ and }\ i\neq j.\\\ \end{cases}$

The solution $\hat{\Sigma}_{G}$ depends on the data only through the sample
covariance matrix $\hat{\Sigma}_{n}$. This approach is due to Dempster,
(1972), and the optimization problem (2) has since been thoroughly studied (e.
g. Speed and Kiiveri,, 1986). Algorithms to compute $\hat{\Sigma}_{G}$ for
arbitrary graphs $G$ can be found, e. g., in Lauritzen, (1996, Chapter 5) or
Hastie et al., (2009, Chapter 17). For decomposable graphs $G$, there is also
an explicit solution, i. e., $\hat{\Sigma}_{G}$ can be computed in a finite
number of steps. For details see also Lauritzen, (1996, Chapter 5). In R (R
Core Team,, 2022), these algorithms are implemented in the function
`fitConGraph` in the package ggm (Marchetti et al.,, 2020). With a graph-
constrained111i. e., it obeys the zero-pattern in the inverse induced by $G$
estimate $\hat{\Sigma}_{G}$ for the covariance matrix available, the remaining
non-zero partial correlations are computed from $\hat{\Sigma}_{G}$ as
unconstrained partial correlations are computed from $\hat{\Sigma}_{n}$ by
virtue of (1).

#### Task (T2)

Within the parametric framework described above, the likelihood-ratio test for
testing $G$ against the full model is given by the test statistic

$D^{\Sigma}_{n}(G)=n\left(\log\det\hat{\Sigma}_{G}-\log\det\hat{\Sigma}_{n}\right),$
(3)

which, under the null hypothesis that $G$ is the true graph, converges to a
$\chi_{q}^{2}$ distribution as $n\to\infty$, where $q$ is the number of
missing edges in $G$. The quantity $D_{n}^{\Sigma}$ is also called _deviance_
and this likelihood-ratio test hence _deviance test_. It simultaneously tests
the absence of all edges not in $G$ avoiding any multiple-testing problems.
The deviance is also returned by `fitConGraph`.

#### Task (T1)

While statistical theory provides rather precise and unambiguous solutions to
tasks (T3) and (T2), this is not the case for (T1), which is already reflected
by the phrasing of _finding an appropriate graph_ rather than _finding the
best fitting graph_. Deciding on an appropriate graph may also be influenced
by interpretability aspects and relevant domain knowledge. A multitude of
approaches exist. A basic idea, also initiated by Dempster, (1972) is the
iterative application of the deviance test. For instance, one starts with the
full graph, then removes one or several edges (with small absolute partial
correlations), and keeps the new candidate graph if the deviance test accepts
it. This may be iterated until no further edge removal leads to an accepted
graph. The opposite search direction is also possible: One starts with the
empty graph and successively adds edges until a graph is obtained which is
accepted by the deviance test. For further reading see the textbooks by
Whittaker, (1990) or Edwards, (2000). Elaborate model search strategies have
been proposed (e. g. Edwards and Havránek,, 1985; Drton and Perlman,, 2008).

Many other model selection approaches use $L_{1}$-regularization and are aimed
at finding sparse graphs in high-dimensional settings. For instance,
Meinshausen and Bühlmann, (2006) propose a node-wise LASSO-regression. Yuan
and Lin, (2007) and Friedman et al., (2008) add an $L_{1}$-penalty for
$K=\Sigma^{-1}$ to the optimization problem (2). Various algorithms have since
been proposed for an efficient computation of such high-dimensional
optimization problems (e.g. Yuan,, 2010; Cai et al.,, 2011; Sun and Zhang,,
2013). Within this framework, the regularization parameter must be chosen,
usually be means of cross validation. Liu and Wang, (2017) particularly
address the latter issue.

###  1.2 Robustness

Robustness in general terms is the property of a statistical method to yield
sensible results if its assumptions are violated. In a more specific sense, it
means insensitivity to outliers. Starting with the pioneering work of Huber,
(1964) and Hampel, (1971, 1974), robust statistics has evolved into a large
research area, see, e. g., the textbooks by Huber and Ronchetti, (2009) or
Maronna et al., (2019).

Figure 1: The non-robust sample covariance matrix (solid line) and the robust
$t_{3}$ M-estimate (dashed line). The data on the left-hand and right-hand
panel differ only by one point.

For our purposes, the important fact to note is that the sample covariance
matrix $\hat{\Sigma}_{n}$ is not robust. In Figure 1, the left-hand panel
shows a small data set of 20 two-dimensional observations. The black ellipse
visualizes the sample covariance matrix, i. e., the 95% probability ellipse of
the thus fitted normal model. In the right-hand panel, one single observation
has been moved from the center to the upper right corner. The covariance
estimate has tremendously changed, suggesting even a positive rather than a
negative correlation. The dashed curve, in contrast, represents an
alternative, robust estimator of multivariate scatter (a $t_{3}$ M-estimator,
see below) and is little altered by the outlier.

With the sample covariance matrix $\hat{\Sigma}_{n}$ being the main ingredient
of essentially all graphical modeling tasks, they all inherit its lack of
robustness. The good news is: robust alternatives exist. The work on robust
multivariate location and scatter estimation as been originated by Maronna,
(1976), who developed Huber’s M-estimation approach for the multivariate
setting. Since then, many proposals have been made.222with numerous
contributions by David Tyler (e. g. Tyler,, 1987; Kent and Tyler,, 1996). Here
we consider only one rather simple and easy-to-compute robust scatter
estimator, the $t_{\nu}$ M-estimator, which is also already mentioned in
Maronna’s paper. This is an M-estimator, whose loss function stems from the
maximum-likelihood estimator within the elliptical $t_{\nu}$-model. The
parameter $\nu$ is referred to as the _degrees of freedom_ and may be any
positive real number. The smaller the $\nu$, the heavier-tailed the $t_{\nu}$
distribution and, consequently, the more outlier-resistant the corresponding
M-estimator. The parameter $\nu$ is usually not inferred from the data, but
selected by the data analyst. A common choice is $\nu=3$. This is not
extremely heavy-tailed, second moments are finite (and hence the covariance
matrix is properly defined), but it is sufficiently heavy-tailed to yield a
strongly outlier-resistant estimator.

The $t_{\nu}$ M-estimator of scatter $\hat{S}_{n}$, along with the
corresponding estimate of location $\hat{\boldsymbol{\mu}}_{n}$, is defined as
the solution to the optimization problem

$(\hat{\boldsymbol{\mu}}_{n},\hat{S}_{n})=\operatorname*{arg\,min}_{\boldsymbol{\mu}\in\mathbb{R}^{p},S\in\mathscr{S}^{+}_{p}}\left[\sum\nolimits_{i=1}^{n}\rho_{\nu,p}\left\\{(\boldsymbol{x}_{i}-\boldsymbol{\mu})^{\top}S^{-1}(\boldsymbol{x}_{i}-\boldsymbol{\mu})\right\\}+n\log\det
S\right]$ (4)

with $\rho_{\nu,p}(x)=(\nu+p)\log(1+x/\nu)$. This yields the estimation
equations

$\begin{cases}\displaystyle\
0=\sum\nolimits_{i=1}^{n}\psi_{\nu,p}(\hat{r}_{i})(\boldsymbol{x}_{i}-\hat{\boldsymbol{\mu}}_{n}),\\\\[8.0pt]
\displaystyle\
\hat{S}_{n}=n^{-1}\sum\nolimits_{i=1}^{n}\psi_{\nu,p}(\hat{r}_{i})(\boldsymbol{x}_{i}-\hat{\boldsymbol{\mu}}_{n})(\boldsymbol{x}_{i}-\hat{\boldsymbol{\mu}}_{n})^{\top},\\\
\end{cases}$

where $\psi_{\nu,p}(x)=\rho^{\prime}_{\nu,p}(x)=(\nu+p)/(\nu+x)$ and
$\hat{r}_{i}=(\boldsymbol{x}_{i}-\hat{\boldsymbol{\mu}}_{n})^{\top}\hat{S}_{n}^{-1}(\boldsymbol{x}_{i}-\hat{\boldsymbol{\mu}}_{n})$.
The $t_{\nu}$ M-estimator can be computed, e. g., by a fixed-point algorithm.
Implementations in R can be found, e. g., in the functions `cov.trob` from the
package MASS (Venables and Ripley,, 2002), `tM` from the package ICS
(Nordhausen et al.,, 2008), and `MVTMLE` from the package fastM (Dümbgen et
al.,, 2016, 2018). The latter uses a partial Newton-Raphson algorithm.

With robust alternatives $\hat{S}_{n}$ being available, an intuitive path to a
robust analysis is the plug-in approach: First solving (4) and then plugging
the thus obtained estimate $\hat{S}_{n}$ instead of $\hat{\Sigma}_{n}$ into
(2) to obtain a graph-constrained robust estimate $\hat{S}_{G}$ and the
corresponding robust (pseudo-)deviance

$D^{S}_{n}(G)=n\left(\log\det\hat{S}_{G}-\log\det\hat{S}_{n}\right).$ (5)

This is accomplished in the package robFitConGraph by the function of the same
name. Plug-in robustifications for the $L_{1}$-regularization methods have
equally been proposed (Öllerer and Croux,, 2015; Tarr et al.,, 2016).
Alternatively, Finegold and Drton, (2011) regularize the elliptical
$t_{\nu}$-log-density.

Before providing further details on the function `robFitConGraph` and the
underlying theory in Section 3, its use shall be demonstrated at a data
example on music performance anxiety.

##  2 A Case Study: Music Performance Anxiety

The fear about one’s ability to perform a specific task, such as giving a
presentation or sitting an exam, affects almost everyone. Pressure can be
particularly high in certain professions where performing in front of others
is an integral part of day-to-day life. While some levels of stress and
anxiety are normal and actually help us to achieve optimal performance, severe
levels of stress and anxiety are debilitating and can develop into a disorder.
Professional musicians are often exposed to extreme pressure where maintaining
top quality performances are not just essential to keeping their job, but to
progress in their careers. Music performance anxiety (MPA) can be understood
as a continuum ranging from low to high anxiety levels. The latter poses a
serious problem to the profession and is the subject of ongoing clinical
research, see, e. g. Fernholz et al., (2019) for a recent review.

MPA is often considered to be a form of social anxiety which, loosely
speaking, is the overwhelming fear of social situations (e. g. Kenny,, 2011;
Cox and Kenardy,, 1993; Nicholson et al.,, 2015; Dobos et al.,, 2019). This is
underlined by the fact that the Diagnostic and Statistical Manual of Mental
Disorders (DSM-5; American Psychiatric Association,, 2013) now acknowledges
evidence of individuals suffering exclusively from performance anxiety as a
distinct sub-type of social anxiety disorder. However, some researchers and
clinicians have questioned this description, as MPA is a complex phenomenon
caused by the interaction of many different factors, and the fear of social
judgment is not necessarily always the main problem.

Wiedemann et al., (2022) analyzed a data set consisting of eight numerical
variables measured at $n=82$ students at German music colleges: music
performance anxiety (MPA), agoraphobia (AG), generalized anxiety disorder
(GAD), panic disorder (PD), separation anxiety disorder (SEP), specific phobia
(SP), social anxiety disorder (SAD) as well as illness anxiety disorder (ILL).
Each variable is a summary score from a self-assessment inventory with Likert-
scale items. A higher value signifies a higher severity of the condition. MPA
was assessed using the German version of the Kenny Music Performance Anxiety
Inventory (K-MPAI; Kenny,, 2009) translated by Spahn et al., (2016). All other
anxieties were assessed using the German translation of the disorder-specific
anxiety measures (Beesdo-Baum et al.,, 2012; Lebeau et al.,, 2012) for the
dimensional anxiety scales of the DSM-5. The data set is also included in the
package robFitConGraph as `anxieties`.

Figure 2: Pairwise scatter plots of the anxieties data

Figure 2 shows the pairwise scatter plots of the data set. While most
participants score low on most anxiety scales – as we would expect and hope –,
there are a few very high values in all variables. The normality assumption
can be seen to be violated in a similar manner as in Figure 1 with the same
implications for any sample-covariance-based analysis. For instance, removing
the outlier in the variable AG reduces the Pearson correlation between AG and
SEP from 0.612 to 0.407. Computing the correlation coefficient from a $t_{3}$
M-estimator, we obtain the value 0.410, which reduces to 0.350 when removing
said outlier. A robust analysis is highly recommended for this data. All
results reported in the following are based on a $t_{3}$ M-estimator of
scatter. It yields the correlation coefficients given in Table 1. They can
also be obtained by the function `robFitConGraph` by supplying the full model
as adjacency matrix

    > library(robFitConGraph)
    > data(anxieties)
    > p <- ncol(anxieties)
    > Shat <- robFitConGraph(X = anxieties,
    +     amat = matrix(1, ncol = p, nrow = p),
    +     df = 3)$Shat
    > round(cov2cor(Shat), d = 2)

Table 1: Pairwise correlations computed from $t_{3}$-M-estimate |  MPA |  GAD |  SAD |  PD |  AG |  SP |  SEP |  ILL
---|---|---|---|---|---|---|---|---
MPA |  $\,\cdot\,$ |  0.62 |  0.37 |  0.43 |  0.17 |  0.26 |  0.38 |  0.32
GAD |  0.62 |  $\,\cdot\,$ |  0.66 |  0.65 |  0.40 |  0.35 |  0.60 |  0.47
SAD |  0.37 |  0.66 |  $\,\cdot\,$ |  0.50 |  0.51 |  0.36 |  0.64 |  0.42
PD |  0.43 |  0.65 |  0.50 |  $\,\cdot\,$ |  0.51 |  0.44 |  0.48 |  0.46
AG |  0.17 |  0.40 |  0.51 |  0.51 |  $\,\cdot\,$ |  0.49 |  0.41 |  0.36
SP |  0.26 |  0.35 |  0.36 |  0.44 |  0.49 |  $\,\cdot\,$ |  0.37 |  0.34
SEP |  0.38 |  0.60 |  0.64 |  0.48 |  0.41 |  0.37 |  $\,\cdot\,$ |  0.29
ILL |  0.32 |  0.47 |  0.42 |  0.46 |  0.36 |  0.34 |  0.29 |  $\,\cdot\,$

All variables are positively correlated. This is neither surprising nor
uncommon for multivariate data. The corresponding partial correlations, as
given in Table 2, shed further light on the dependence structure of the
variables.

    > Phat <- -cov2cor(solve(Shat))
    > diag(Phat) <- 1
    > round(Phat,d = 2)

Table 2: Pairwise _partial_ correlations computed from $t_{3}$-M-estimate

Many of the partial correlations are near zero, suggesting conditional
independences, i. e., their association may be fully mediated by other
variables in the data set.

###  2.1 Inferential analysis: MPA and social anxiety

One hypothesis examined by Wiedemann et al., (2022) is whether MPA is
primarily related to a social anxiety disorder (SAD). Based on the positive
correlation of $0.37$, a short-sighted analysis may deduce a strong connection
between MPA and SAD, which in light of other, even larger correlations, is
right away questionable. As we will see, the data in fact carries no evidence
for a particularly strong connection between MPA and SAD.

For that purpose one tests the hypothesis that MPA and SAD are conditionally
independent given all remaining six variables, i. e., testing a graphical
model with only one missing edge between MPA and SAD by means a robust pseudo-
deviance test. If this test rejects, there is evidence for a strong link
between MPA and SAD. Their dependence can then not be fully explained by their
associations with other anxiety types. However, this hypothesis is _not_
rejected (a p-value of 0.66), and we find no such evidence in the data. The
p-value is returned by the function `robFitConGraph` via the named list
element `pval`, see also the next code chunk below.

One can even further test the stronger hypothesis that MPA is, given GAD only,
conditionally independent of all other six specific anxiety scales. This is
done by testing the graph in Figure 3, which can equivalently be expressed by
its adjacency matrix on the right-hand side of Figure 3.

Almost unnoticed we make here use of a non-trivial result, which is the real
merit and the real beauty of graphical models: Each absent edge in the graph
denotes a conditional independence of two individual variables given the
respective remaining six variables. This is indeed equivalent to MPA and (SAD,
PD, AG, SP, SEP, ILL) being conditionally independent given GAD only, because,
in the graph, GAD separates MPA from (SAD, PD, AG, SP, SEP, ILL). This is
known as the equivalence between the local Markov property and the global
Markov property. For details see Lauritzen, (1996, Chapter 3).

Figure 3: The graph encoding the hypothesis MPA is conditionally independent
of all other variables given GAD (left) and the corresponding adjacency matrix
(right).

|  MPA |  GAD |  SAD |  PD |  AG |  SP |  SEP |  ILL
---|---|---|---|---|---|---|---|---
MPA |  1 |  1 |  0 |  0 |  0 |  0 |  0 |  0
GAD |  1 |  1 |  1 |  1 |  1 |  1 |  1 |  1
SAD |  0 |  1 |  1 |  1 |  1 |  1 |  1 |  1
PD |  0 |  1 |  1 |  1 |  1 |  1 |  1 |  1
AG |  0 |  1 |  1 |  1 |  1 |  1 |  1 |  1
SP |  0 |  1 |  1 |  1 |  1 |  1 |  1 |  1
SEP |  0 |  1 |  1 |  1 |  1 |  1 |  1 |  1
ILL |  0 |  1 |  1 |  1 |  1 |  1 |  1 |  1

The pseudo-deviance test based on the $t_{3}$ M-estimator for the graph in
Figure 3 is carried out as follows:

    > amat <- matrix(1, ncol = p, nrow = p)
    > rownames(amat) <- colnames(anxieties)
    > colnames(amat) <- colnames(anxieties)
    > amat["MPA", c("SAD", "PD", "AG", "SP", "SEP", "ILL")] <- 0
    > amat[c("SAD", "PD", "AG", "SP", "SEP", "ILL"), "MPA"] <- 0
    > robFitConGraph(X = anxieties, amat = amat, df = 3)$pval
    [1] 0.881159

With a p-value of 0.88, this hypothesis is also _not_ rejected. There is no
evidence against the null hypothesis of MPA being conditionally independent of
the specific anxiety types given generalized anxiety (GAD). Accepting the null
hypothesis, GAD is fully sufficient for predicting MPA. When already knowing a
person’s GAD score, i. e., how anxious the person generally is, additionally
knowing any of their specific anxieties scores provides no further information
about their MPA score. So this analysis does not support the hypothesis that
MPA is foremost related to social anxiety.

The analysis is complemented by testing if MPA and GAD are conditionally
independent given the remaining six variables. This hypothesis _is_ rejected
with a p-value below 0.01, which corresponds to the partial correlation of
0.44 in Table 2.333For a single missing edge, the deviance test can indeed be
expressed in terms of corresponding partial correlation only. So there is a
strong connection between MPA and GAD, this connection is not mediated by any
of the other variables, and GAD explains the connection between MPA and the
remaining variables. Altogether, GAD is the link between MPA and the other
anxieties.

###  2.2 Explorative analysis

Naturally, the question arises which other edges may be removed. For that
purpose we remove all edges which have absolute partial correlation below 0.15
(Table 2), corresponding to an individual p-value above 0.2. The resulting
graph is depicted in Figure LABEL:fig:4, which is not rejected by the pseudo-
deviance test with a p-value of 0.45.

Figure 4: A fitting graph with 9 edges based on a $t_{3}$ scatter estimator
with p-value 0.45.

The partial correlations along the edges in Figure LABEL:fig:4 are fitted
under the graph, i. e., they are estimated taking the graph structure into
account. They are different from the partial correlations given in Table 2.

    > amat <- abs(Phat) > 0.15
    > Shat_G <- robFitConGraph(X = anxieties,
    +     amat = amat, df = 3)$Shat
    > Phat_G <- -cov2cor(solve(Shat_G))
    > diag(Phat_G) <- 1
    > round(Phat_G, d = 2)ΨΨΨΨΨΨΨΨΨ

The information learned from such a fitted graph is, e. g., that GAD is
central within this set of variables: It has many edges to other vertices, it
has much explanatory power about the other variables. If one were to retain
only a single variable (as a very simple dimension reduction approach, say),
GAD would be a natural candidate based on this “vertex degree criterion”.

The fitted graph in Figure LABEL:fig:4 with 9 edges (out of 28 possible) is
one parsimonious graph that fits the data. Generally, there is no well-defined
_most parsimonious_ , _best-fitting_ graph, as obviously fit and parsimony are
contradicting goals. However, one may elaborate upon the initial search by
checking if other (equally or more) parsimonious graphs may fit as well. For
instance, increasing the initial partial correlation threshold to 0.17, say,
which results in a further removal of the edge ILL–PD, leads to a p-value of
0.01 and hence should be rejected. Alternatively, removing the “next smallest
edge” from the candidate graph in Figure LABEL:fig:4, which is GAD–SEP,
results in a p-value of 0.17, which is still acceptable.

Despite using deviance-test p-values as decision criterion whether to accept a
graph or not, this analysis is purely explorative. We have not _tested_ the
graph of Figure LABEL:fig:4. Hypotheses to be tested must be formed a-priori.
Testing for a graph which is the result of model selection procedure is as
prohibitive as testing if the observed sample mean is the population mean.

So far, we have used `robFitConGraph` with `df = 3`, which is also the default
setting in case `df` is not specified. Generally, the results are not very
sensitive to variations in $\nu$. A smaller value of $\nu$ downweights
outliers more strongly. In the present example, taking $\nu=1$, we find the
hypothesis of Figure 3 equally accepted with a p-value of 0.58. The
explorative graph in Figure LABEL:fig:4 is accepted with a p-value of 0.26.

###  2.3 The classical analysis

After having cautioned its use at the beginning, we close this section by
remarking that a classical sample-covariance-based analysis leads to
qualitatively the same findings. The classical deviance test gives a p-value
of 0.63 for the hypothesized graph of Figure 3. An explorative analysis
analogous to the one above leads to the graph in Figure LABEL:fig:5 with a
p-value of 0.51, which is also communicated by Wiedemann et al., (2022).

Figure 5: A fitting graph with 10 edges based on the sample covariance matrix
with p-value 0.51.

The graph is different, but it equally shows the centrality of GAD as well as
the marginality of MPA within the set of variables. It is generally advisable
to try varying parameter settings and different methods, robust and non-
robust, in any statistical analysis. Conclusions unanimously obtained by
several methods may be considered even more trustworthy.

##  3 Background and theory

The centerpiece of the package robFitConGraph is the function by the same
name. For a given data set $\mathbb{X}_{n}$ and a given graph $G$, it
simultaneously provides a robust graph-constrained matrix $\hat{S}_{G}$ and
the p-value of the corresponding pseudo-deviance test. But so far, the actual
worth of the function `robFitConGraph` has not become apparent: It appears,
its plug-in functionality is equally achieved by a combination of, say,
`cov.trob` and `fitConGraph`, as the latter takes the sample covariance matrix
as input. However, there are two good reasons for `robFitConGraph`, which both
require a slightly deeper foray into statistical theory.

  1. (1)

The limit distribution of the pseudo-deviance $D^{S}_{n}(G)$ does not follow
strictly a $\chi^{2}_{q}$-distribution under the graph $G$, but a
$\chi^{2}_{q}$-variate multiplied by a constant $\sigma_{1}$. This constant is
required to obtain p-values and depends on the data dimension $p$ and the
degrees of freedom $\nu$ of the estimator. It is computed by the function
`find_sigma1`.

  2. (2)

The plug-in estimator $\hat{S}_{G}$ is just one approach to a graph-
constrained robust scatter estimator. The function `robFitConGraph` also
provides an alternative approach, referred to as the _direct estimator_ and
denoted by $\tilde{S}_{G}$ below.

Another good reason is speed. Being implemented in C++, `robFitConGraph` is
considerably faster than the combination of `cov.trob` and `fitConGraph`. A
runtime comparison is given by Watt, (2019).

###  3.1 The Constant $\sigma_{1}$

Recall the class of all $p$-dimensional, continuous, elliptical distributions,
i. e., distributions possessing a $p$-dimensional Lebesgue density $f$ of the
form

$f(\boldsymbol{x})=\det(S)^{-\frac{1}{2}}g\big{\\{}(\boldsymbol{x}-\boldsymbol{\mu})^{\top}S^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\big{\\}}$
(6)

for some $\boldsymbol{\mu}\in\mathbb{R}$, $S\in\mathscr{S}^{+}_{p}$ and
$g:[0,\infty)\to[0,\infty)$, such that $f$ integrates to 1. Let
$E_{p}(\boldsymbol{\mu},S,g)$ denote the distribution described by (6). The
univariate function $g$ is called the _elliptical generator_ and $S$ the
_scatter_ or _shape matrix_ of $E_{p}(\boldsymbol{\mu},S,g)$. It is
proportional to the covariance matrix if $E_{p}(\boldsymbol{\mu},S,g)$ has
second moments. In the present paper, we only consider two examples of
elliptical distributions: the normal distribution with

$g_{N}(x)=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{x}{2}\right),\qquad 0\leq x,$

and the elliptical $t_{\nu}$-distribution with

$g_{\nu,p}(x)=c_{\nu,p}\left(1+\frac{x}{\nu}\right)^{-(\nu+p)/2},\qquad 0\leq
x.$

where the normalizing constant $c_{\nu,p}$ is given in the Technical Appendix
at the end of the chapter. Assume the independent observations
$\mathbb{X}_{n}=(\boldsymbol{x}_{1},\ldots\boldsymbol{x}_{n})^{\top}$ to stem
from an elliptical distribution $E_{p}(\boldsymbol{\mu},S,g)$ and let
$\hat{V}_{n}$ be an _arbitrary_ $p\times p$ scatter matrix estimator
fulfilling two fairly natural conditions:

  1. (1)

$\hat{V}_{n}$ is affine equivalent, i. e.,
$\hat{V}_{n}(\mathbb{X}_{n}A^{\top}+\boldsymbol{1}_{n}\boldsymbol{b}^{\top})=A\hat{S}_{n}(\mathbb{X}_{n})A^{\top}$
for any $\boldsymbol{b}\in\mathbb{R}^{p}$ and full rank
$A\in\mathbb{R}^{p\times p}$, where $\boldsymbol{1}_{n}$ is the $n$-vector
consisting of ones.

  2. (2)

$\hat{V}_{n}$ is $\sqrt{n}$-consistent and asymptotically normal at
$E_{p}(\boldsymbol{\mu},S,g)$, i. e., there exist matrices
$V\in\mathscr{S}_{p}$ and and $W\in\mathscr{S}_{p^{2}}$ such that

$\sqrt{n}\rm vec\\{\hat{V}_{n}(\mathbb{X}_{n})-V\\}\ \to\
N_{p^{2}}\\!\left(\boldsymbol{0},\ W\right)$

in distribution as $n\to\infty$.

Tyler, (1982) showed that under these two conditions, $V$ and $W$ take on
specific forms: $V=\eta S$ for some $\eta>0$, and

$W=2\eta^{2}\sigma_{1}\mathcal{M}_{p}(S\otimes S)\ +\ \eta^{2}\sigma_{2}\rm
vec(S)\rm vec(S)^{\top},$ (7)

where $\otimes$ is the Kronecker product, $\sigma_{1}\geq 0$ and
$\sigma_{2}\geq-2\sigma_{1}/p$ are scalar constants independent of
$\boldsymbol{\mu}$ and $S$, and $\mathcal{M}_{p}$ is a fixed $p^{2}\times
p^{2}$ matrix defined in the Technical Appendix. The latter formula greatly
simplifies the asymptotic efficiency comparison of any two affine equivariant
scatter estimators at elliptical distributions. In case of $\hat{V}_{n}$ being
an elliptical M-estimator, with a general loss function $\rho$ instead of
$\rho_{\nu,p}$ in (4), Tyler, (1983) gives the following expressions for the
scalars $\eta$, $\sigma_{1}$ and $\sigma_{2}$: Letting
$\psi(x)=\rho^{\prime}(x)$, $\phi(x)=x\psi(x)$, and
$R=(\boldsymbol{X}-\boldsymbol{\mu})^{\top}S^{-1}(\boldsymbol{X}-\boldsymbol{\mu})$
for $\boldsymbol{X}\sim E_{p}(\boldsymbol{\mu},S,g)$, the scalar $\eta$ is the
solution to $\mathrm{\mathbf{E}}\\{\phi(R/\eta)\\}=p$. Letting further

$\gamma_{1}=\frac{\mathrm{\mathbf{E}}\\{\phi^{2}(R/\eta)\\}}{p(p+2)},\qquad\gamma_{2}=\frac{1}{p}\mathrm{\mathbf{E}}\left\\{\frac{R}{\eta}\phi^{\prime}\left(\frac{R}{\eta}\right)\right\\},$

the scalars $\sigma_{1}$ and $\sigma_{2}$ are

$\sigma_{1}=\frac{(p+2)^{2}\gamma_{1}}{(2\gamma_{2}+p)^{2}},\qquad\sigma_{2}=\gamma_{2}^{-2}\left[\gamma_{1}-1-\frac{2\gamma_{1}(\gamma_{2}-1)\\{p+(p+4)\gamma_{2}\\}}{(2\gamma_{2}+p)^{2}}\right].$

Tyler, (1983) further showed that the asymptotic variance of any scale-
invariant, continuously differentiable function $h$ of $\hat{V}_{n}$ only
depends $\sigma_{1}$ and not on $\sigma_{2}$. A _scale-invariant_ function
$h:\mathscr{S}_{p}\to\mathbb{R}$ satisfies $h(\alpha V)=h(V)$ for any
$\alpha>0$. The pseudo-deviance test statistic (5) is such a scale-invariant
function of the scatter estimator $\hat{S}_{n}$. Dependence is an inherently
scale-free concept. So this equally applies to any aspect of multivariate
scatter that quantifies dependence in one way or another, may it be
correlations, partial correlations, canonical correlations, principal
components, etc.

In the package robFitConGraph, the functions `find_eta` and `find_sigma1`
compute the scalars $\eta$ and $\sigma_{1}$, respectively, for $t_{\nu_{1}}$
M-estimators in case the data stem from a normal or an elliptical
$t_{\nu_{2}}$ distribution. Note that the degrees of freedom $\nu_{1}$ of the
estimator loss function $\rho_{\nu_{1},p}$ and the degrees of freedom of the
population distribution $\nu_{2}$ can be generally different. If they
coincide, the $t_{\nu_{1}}$ M-estimator is a maximum-likelihood estimator.

The value of $\nu_{1}$ is specified by `df_est` and $\nu_{2}$ by `df_data`.
For both, `Inf` is allowed, which corresponds to the sample covariance matrix
and the normal distribution, respectively. For `find_sigma1`, the input
`df_est = 0` is also allowed. This corresponds to Tyler’s distribution-free
M-estimator of scatter (Tyler,, 1987). In this case $\sigma_{1}=1+2/p$
regardless of the elliptical population distribution. There is no $t_{0}$
distribution, hence `df_data = 0` is not allowed.

The value $\sigma_{1}$ can be given directly to `robFitConGraph` via the
optional argument `sigma1`. If none is provided, `robFitConGraph` calls
`find_sigma1` with `df_data = Inf` (i. e., assuming Gaussian data) and its
argument `df` being passed on as `df_est` to the function `find_sigma1`. The
argument `df` of `robFitConGraph` is optional with the default `df = 3`.

###  3.2 The Direct vs. the Plug-in Estimate

Instead of solving (4) and (2) sequentially, an alternative approach is to
directly solve the optimization problem

$(\tilde{\boldsymbol{\mu}}_{G},\tilde{S}_{G})=\operatorname*{arg\,min}_{\boldsymbol{\mu}\in\mathbb{R}^{p},S^{-1}\in\mathscr{S}^{+}_{p}(G)}\left[\sum\nolimits_{i=1}^{n}\rho_{\nu,p}\\!\left\\{(\boldsymbol{x}_{i}-\boldsymbol{\mu})^{\top}S^{-1}(\boldsymbol{x}_{i}-\boldsymbol{\mu})\right\\}+n\log\det
S\right],$ (8)

with $\rho_{\nu,p}$ as in (4), leading to the estimation equations

$\begin{cases}\displaystyle\
0=\sum\nolimits_{i=1}^{n}\psi_{\nu,p}(\tilde{r}_{i})(\boldsymbol{x}_{i}-\tilde{\boldsymbol{\mu}}_{G}),\\\\[8.0pt]
\displaystyle\
\left[\tilde{S}_{G}\right]_{j,k}=\left[n^{-1}\sum\nolimits_{i=1}^{n}\psi_{\nu,p}(\tilde{r}_{i})(\boldsymbol{x}_{i}-\tilde{\boldsymbol{\mu}}_{G})(\boldsymbol{x}_{i}-\tilde{\boldsymbol{\mu}}_{G})^{\top}\right]_{j,k},&\mbox{for
}\\{j,k\\}\in E\ \mbox{ or }\ j=k,\\\\[10.0pt] \displaystyle\
\left[\tilde{S}_{G}^{-1}\right]_{j,k}=0,&\mbox{for }\\{j,k\\}\notin E\ \mbox{
and }\ j\neq k,\end{cases}$ (9)

where
$\tilde{r}_{i}=(\boldsymbol{x}_{i}-\tilde{\boldsymbol{\mu}}_{G})^{\top}\tilde{S}_{n}^{-1}(\boldsymbol{x}_{i}-\tilde{\boldsymbol{\mu}}_{G})$
and, as before, $\psi_{\nu,p}(x)=\rho^{\prime}_{\nu,p}(x)$.

The estimator $\tilde{S}_{G}$ shall be called the _direct estimator_ , which
is short for _direct graph-constrained $t_{\nu}$ M-estimator_, and is an
alternative to the plug-in graph-constrained $t_{\nu}$ M-estimator. Using the
function `robFitConGraph`, the direct estimator is invoked by setting the
option `direct = TRUE` or `plug_in = FALSE`. In case of conflicting
specifications, `plug_in` has priority, and a message will be displayed.
Contrary to the plug-in estimator $\hat{S}_{G}$, the direct estimator
$\tilde{S}_{G}$ is not a function of the corresponding unconstrained estimate
$\hat{S}_{n}$ alone.

One main theoretical result of Vogel and Tyler, (2014) is the asymptotic
equivalence of $\hat{S}_{G}$ and $\tilde{S}_{G}$ under elliptical population
distributions. Hence the limiting distribution of the pseudo-deviance and the
constant $\sigma_{1}$ are the same in both cases. It may be argued that this
asymptotic equivalence result favors the plug-in estimator: Considering the
elliptical distribution (6) with the fixed generator $g_{\nu,p}$, one defines
the elliptical-$t_{\nu}$ graphical model $E_{p}(g_{\nu,p},G)$ analogously to
the Gaussian graphical model $N_{p}(G)$. The direct graphical $t_{\nu}$
M-estimator $(\tilde{\boldsymbol{\mu}}_{G},\tilde{S}_{G})$ with the
corresponding loss function $\rho_{\nu,p}$ is then the maximum-likelihood
estimator within this parametric family.444The parameters of interest are
$\boldsymbol{\mu}$ and $S$. The degrees of freedom $\nu$ is held fixed. So
direct graphical M-estimators generalize maximum-likelihood estimators, which
are known to be first-order efficient. Plug-in estimators are popular with
practitioners as they are easily applied, fast to compute, and now it turns
out that they also possess desirable asymptotic properties. Direct graphical
M-estimators are generally harder to compute. They may be solved by a double-
loop, iteratively reweighted least-squares algorithm, where the Gaussian model
fit is nested into the M-estimation loop.

The other half of the story is that direct graphical M-estimators can be
substantially more efficient in small samples as is demonstrated by
simulations in Vogel and Tyler, (2014).

###  3.3 Ellipticity vs. Normality

Multivariate data containing outliers may be modeled conceptually in two
different ways: Either by a corrupted Gaussian distribution, i. e., a few
observations are erroneous and stem from a different, outlier-generating
distribution, or by a heavy-tailed elliptical distribution, which generates
outliers itself. Using the $t_{\nu}$ M-estimator implicitly suggests the
latter viewpoint. However, two issues arise:

  1. (1)

Many data sets exhibit features such as the anxiety data set: It is clearly
not normal as it contains outliers, but it is clearly not elliptical either as
it is skewed. Hence we do adopt the viewpoint of corrupted normal data. In the
analysis, however, we apply outlier-resistant methods which have been derived
from considerations in the elliptical $t_{\nu}$ model. We work on the
plausible, but not formalized assumption that they provide outlier-resistance
regardless of the outliers being scattered symmetrically around the center or
not.

  2. (2)

An absent edge in a Gaussian graphical model, i. e., a zero partial
correlation, has the interpretation of conditional independence, a notion we
have repeatedly used in Section 2. Under ellipticity, an absent edge has the
slightly weaker interpretation of conditional uncorrelatedness. This is
occasionally mentioned as a limitation of elliptical graphical models.
However, this limitation is largely void. The conclusions of any statistical
analysis are always a combination of the information contained in the data and
the modeling assumptions. When performing a purely linear analysis based on
partial correlations, concluding conditional independences due the normality
assumption certainly falls into the latter category.

## Technical Appendix

The normalizing constant in the elliptical $t_{\nu}$-density is

$c_{\nu,p}=\frac{\Gamma\\{(\nu+p)/2\\}}{(\nu\pi)^{p/2}\Gamma(\nu/2)},$

where $\Gamma$ is the gamma function, see e. g. Bilodeau and Brenner, (1999,
Chapter 13).

The Kronecker product $A\otimes B$ of two matrices $A,B\in\mathbb{R}^{p\times
p}$ is defined as the $p^{2}\times p^{2}$ matrix with entry $a_{i,j}b_{k,l}$
at position $((i-1)p+k,(j-1)p+l)$. For a matrix
$A=(\boldsymbol{a}_{1},\ldots,\boldsymbol{a}_{p})\in\mathbb{R}^{p\times p}$,
the notation $\rm vec(A)$ means the $p^{2}$-vector obtained by stacking the
columns of $A$, i. e., $\rm
vec(A)=(\boldsymbol{a}_{1}^{\top},\ldots,\boldsymbol{a}_{p}^{\top})^{\top}$.
The matrix $\mathcal{M}_{p}$ in (7) is defined as

$\mathcal{M}_{p}=\frac{1}{2}\left(\mathcal{I}_{p^{2}}+\mathcal{K}_{p}\right),$

where $\mathcal{I}_{p^{2}}$ denotes the $p^{2}\times p^{2}$ identity matrix
and

$\mathcal{K}_{p}=\sum\nolimits_{i=1}^{p}\sum\nolimits_{j=1}^{p}\boldsymbol{e}_{i}\boldsymbol{e}_{j}^{\top}\otimes\boldsymbol{e}_{j}\boldsymbol{e}_{i}^{\top},$

where $\boldsymbol{e}_{1},\ldots,\boldsymbol{e}_{p}$ denote the Euclidean
basis vectors in $\mathbb{R}^{p}$. The matrix $\mathcal{K}_{p}$, commonly
referred to as the _commutation matrix_ , is orthogonal and corresponds to the
transpose operator $\mathcal{K}_{p}\rm vec(A)=\rm vec(A^{\top})$. The
idempotent matrix $\mathcal{M}_{p}$ is called the _symmetrization matrix_
since it maps $\rm vec(A)$ to $\rm vec(A+A^{\top})/2$.

## Acknowledgment

This article appears in the festschrift for David E. Tyler. We thank the
anonymous referees for many valuable suggestions and spotting mistakes in the
originally submitted manuscript.

## References

  * American Psychiatric Association, (2013) American Psychiatric Association (2013).  Diagnostic and statistical manual of mental disorders.  Arlington, VA: American Psychiatric Publishing, fifth edition.
  * Beesdo-Baum et al., (2012) Beesdo-Baum, K., Klotsche, J., Knappe, S., Craske, M. G., LeBeau, R. T., Hoyer, J., Strobel, A., Pieper, L., and Wittchen, H.-U. (2012).  Psychometric properties of the dimensional anxiety scales for DSM-V in an unselected sample of german treatment seeking patients.  Depression and Anxiety, 29(12):1014–1024.
  * Bilodeau and Brenner, (1999) Bilodeau, M. and Brenner, D. (1999).  Theory of Multivariate Statistics. Springer Texts in Statistics. New York: Springer.
  * Cai et al., (2011) Cai, T., Liu, W., and Luo, X. (2011).  A constrained $l_{1}$ minimization approach to sparse precision matrix estimation.  Journal of the American Statistical Association, 106(494):594–607.
  * Cox and Kenardy, (1993) Cox, W. J. and Kenardy, J. (1993).  Performance anxiety, social phobia, and setting effects in instrumental music students.  Journal of Anxiety Disorders, 7(1):49–60.
  * Dempster, (1972) Dempster, A. P. (1972).  Covariance Selection.  Biometrics, 28:157–175.
  * Dobos et al., (2019) Dobos, B., Piko, B. F., and Kenny, D. T. (2019).  Music performance anxiety and its relationship with social phobia and dimensions of perfectionism.  Research Studies in Music Education, 41(3):310–326.
  * Drton and Perlman, (2008) Drton, M. and Perlman, M. D. (2008).  A SINful approach to Gaussian graphical model selection.  Journal of Statistical Planning and Inference, 138(4):1179–1200.
  * Dümbgen et al., (2016) Dümbgen, L., Nordhausen, K., and Schuhmacher, H. (2016).  New algorithms for M-estimation of multivariate scatter and location.  Journal of Multivariate Analysis, 144:200–217.
  * Dümbgen et al., (2018) Dümbgen, L., Nordhausen, K., and Schuhmacher, H. (2018).  fastM: Fast Computation of Multivariate M-Estimators.  R package version 0.0-4.
  * Edwards, (2000) Edwards, D. (2000).  Introduction to graphical modelling. Springer Texts in Statistics. New York, NY: Springer.
  * Edwards and Havránek, (1985) Edwards, D. and Havránek, T. (1985).  A fast procedure for model search in multidimensional contingency tables.  Biometrika, 72:339–351.
  * Fernholz et al., (2019) Fernholz, I., Mumm, J. L., Plag, J., Noeres, K., Rotter, G., Willich, S. N., Ströhle, A., Berghöfer, A., and Schmidt, A. (2019).  Performance anxiety in professional musicians: a systematic review on prevalence, risk factors and clinical treatment effects.  Psychological Medicine, 49(14):2287–2306.
  * Finegold and Drton, (2011) Finegold, M. and Drton, M. (2011).  Robust graphical modeling of gene networks using classical and alternative $t$-distributions.  Annals of Applied Statistics, 5(2A):1057–1080.
  * Friedman et al., (2008) Friedman, J., Hastie, T., and Tibshirani, R. (2008).  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics, 9(3):432–441.
  * Hampel, (1971) Hampel, F. R. (1971).  A general qualitative definition of robustness.  Annals of Mathematical Statistics, 42:1887–1896.
  * Hampel, (1974) Hampel, F. R. (1974).  The influence curve and its role in robust estimation.  Journal of the American Statistical Association, 69:383–393.
  * Hastie et al., (2009) Hastie, T., Tibshirani, R., and Friedman, J. (2009).  The Elements of Statistical Learning: Data Mining, Inference, and Prediction.  Springer Series in Statistics. Springer, New York, second edition.
  * Huber, (1964) Huber, P. J. (1964).  Robust estimation of a location parameter.  Annals of Mathematical Statistics, 35(1):73–101.
  * Huber and Ronchetti, (2009) Huber, P. J. and Ronchetti, E. M. (2009).  Robust statistics. Wiley Series in Probability and Statistics. Hoboken, NJ: Wiley, sencond edition.
  * Kenny, (2009) Kenny, D. T. (2009).  The factor structure of the revised Kenny music performance anxiety inventory.  In International Symposium on performance science, pages 37–41. Association Européenne des Conservatoires Utrecht, The Netherlands.
  * Kenny, (2011) Kenny, D. T. (2011).  The Psychology of Music Performance Anxiety.  Oxford University Press.
  * Kent and Tyler, (1996) Kent, J. T. and Tyler, D. E. (1996).  Constrained M-estimation for multivariate location and scatter.  Annals of Statistics, 24(3):1346–1370.
  * Lauritzen, (1996) Lauritzen, S. L. (1996).  Graphical models. Oxford Statistical Science Series. Oxford University Press.
  * Lebeau et al., (2012) Lebeau, R. T., Glenn, D. E., Hanover, L. N., Beesdo-Baum, K., Wittchen, H.-U., and Craske, M. G. (2012).  A dimensional approach to measuring anxiety for DSM-5.  International Journal of Methods in Psychiatric Research, 21(4):258–272.
  * Liu and Wang, (2017) Liu, H. and Wang, L. (2017).  TIGER: A tuning-insensitive approach for optimally estimating gaussian graphical models.  Electronic Journal of Statistics, 11(1):241–294.
  * Marchetti et al., (2020) Marchetti, G. M., Drton, M., and Sadeghi, K. (2020).  ggm: Graphical Markov Models with Mixed Graphs.  R package version 2.5.
  * Maronna, (1976) Maronna, R. A. (1976).  Robust M-estimators of multivariate location and scatter.  Annals of Statistics, 4:51–67.
  * Maronna et al., (2019) Maronna, R. A., Martin, R. D., Yohai, V. J., and Salibián-Barrera, M. (2019).  Robust statistics: Theory and methods (with R).  John Wiley & Sons, second edition.
  * Meinshausen and Bühlmann, (2006) Meinshausen, N. and Bühlmann, P. (2006).  High-dimensional graphs and variable selection with the Lasso.  Annals of Statistics, 34(3):1436–1462.
  * Nicholson et al., (2015) Nicholson, D. R., Cody, M. W., and Beck, J. G. (2015).  Anxiety in musicians: On and off stage.  Psychology of Music, 43(3):438–449.
  * Nordhausen et al., (2008) Nordhausen, K., Oja, H., and Tyler, D. E. (2008).  Tools for exploring multivariate data: The package ICS.  Journal of Statistical Software, 28(6):1–31.
  * Öllerer and Croux, (2015) Öllerer, V. and Croux, C. (2015).  Robust high-dimensional precision matrix estimation.  In Nordhausen, K. and Taskinen, S., editors, Modern Nonparametric, Robust and Multivariate Methods., pages 325–350. Berlin, Heidelberg: Springer-Verlag.
  * R Core Team, (2022) R Core Team (2022).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Spahn et al., (2016) Spahn, C., Walther, J.-C., and Nusseck, M. (2016).  The effectiveness of a multimodal concept of audition training for music students in coping with music performance anxiety.  Psychology of Music, 44(4):893–909.
  * Speed and Kiiveri, (1986) Speed, T. P. and Kiiveri, H. T. (1986).  Gaussian Markov distributions over finite graphs.  Annals of Statistics, 14:138–150.
  * Sun and Zhang, (2013) Sun, T. and Zhang, C.-H. (2013).  Sparse matrix inversion with scaled lasso.  The Journal of Machine Learning Research, 14(1):3385–3418.
  * Tarr et al., (2016) Tarr, G., Müller, S., and Weber, N. C. (2016).  Robust estimation of precision matrices under cellwise contamination.  Computational Statistics & Data Analysis, 93:404–420.
  * Tyler, (1982) Tyler, D. E. (1982).  Radial estimates and the test for sphericity.  Biometrika, 69:429–436.
  * Tyler, (1983) Tyler, D. E. (1983).  Robustness and efficiency properties of scatter matrices.  Biometrika, 70:411–420.
  * Tyler, (1987) Tyler, D. E. (1987).  A Distribution-Free $M$-Estimator of Multivariate Scatter.  Annals of Statistics, 15(1):234 – 251.
  * Venables and Ripley, (2002) Venables, W. N. and Ripley, B. D. (2002).  Modern Applied Statistics with S.  Springer, New York, fourth edition.  ISBN 0-387-95457-0.
  * Vogel and Tyler, (2014) Vogel, D. and Tyler, D. E. (2014).  Robust estimators for non-decomposable elliptical graphical models.  Biometrika, 101(4):865–882.
  * Watt, (2019) Watt, S. J. (2019).  Algorithms for data analysis.  B.Sc. thesis, University of Aberdeen.
  * Whittaker, (1990) Whittaker, J. (1990).  Graphical Models in Applied Multivariate Statistics. Chichester, etc.: Wiley.
  * Wiedemann et al., (2022) Wiedemann, A., Vogel, D., Voss, C., and Hoyer, J. (2022).  How does music performance anxiety relate to other anxiety disorders?  Psychology of Music, 50:204–2017.
  * Yuan, (2010) Yuan, M. (2010).  High dimensional inverse covariance matrix estimation via linear programming.  The Journal of Machine Learning Research, 11:2261–2286.
  * Yuan and Lin, (2007) Yuan, M. and Lin, Y. (2007).  Model selection and estimation in the Gaussian graphical model.  Biometrika, 94(1):19–35.


# A self-adaptive and robust fission clustering algorithm via heat diffusion
and maximal turning angle

Yu Han1, Shizhan Lu2, Haiyan Xu1∗ 1College of Economics and Management,
Nanjing University of Aeronautics and Astronautics, Nanjing 211106, P.R. China
2School of Economics and Management, Nanjing University of Science and
Technology, Xiaolingwei 200, Nanjing 210094, P.R. China

###### Abstract

Cluster analysis, which focuses on the grouping and categorization of similar
elements, is widely used in various fields of research. A novel and fast
clustering algorithm, fission clustering algorithm, is proposed in recent
year. In this article, we propose a robust fission clustering (RFC) algorithm
and a self-adaptive noise identification method. The RFC and the self-adaptive
noise identification method are combine to propose a self-adaptive robust
fission clustering (SARFC) algorithm. Several frequently-used datasets were
applied to test the performance of the proposed clustering approach and to
compare the results with those of other algorithms. The comprehensive
comparisons indicate that the proposed method has advantages over other common
methods.

###### keywords:

Clustering, Density-based, Self-adaptive, Robust fission clustering algorithm.

##  1 Introduction

Clustering analysis, as a useful technique in data mining focusing on grouping
and categorizing similar elements in a data set, has been widely applied in
diverse research areas, ranging from climate predictions[1], gene
expression[2], bioinformatics[3], finance and economics[4, 5] to
neuroscience[6, 7], to list a few. Benefited from the modern automatic
monitoring technologies, people have gathered more and more data. However,
with poor information of this data we rarely exactly know in advance the
number of clusters before any analysis. Therefore, a clustering method that
does not need setting an input parameter or parameters is preferred, which is
what we aim for in this article.

There have been various different clustering methods in the literature. In
general, they can be classified into 6 categories, that is density-based
(c.f., DP[8], DP-HD[9], DBSCAN[10], NQ-DBSCAN[11], CSSub[12] and GDPC[13]),
grid-based (c.f., CLIQUE[14], Gridwave[15] and WaveCluster[16]), model-based
(c.f., Gaussian parsimonious[17], Gaussian mixture models[18] and Latent tree
models[19]), partitioning (c.f., K-means[20, 21, 22], K-partitioning[23] and
TLBO[24]), graph-based (SEGC[25], ProClust[26] and MCSSGC[27]), and
hierarchical (c.f., BIRCH[28] and CHAMELEON[29]) approaches.

Lu, Cheng and Mehmood[30] proposed a novel and interesting clustering
algorithm which can cluster data without any parameters, fission clustering
(FC) algorithm. FC has several excellent features to compare with some
existing algorithms: (F1), FC algorithm has stronger coping ability for the
extreme challenge data sets to compare with many of current clustering
algorithms; (F2), the time complexity of FC is superior to many of current
existing algorithms; (F3), FC can utilize computing group to reduce the
difficulty of data processing when it faces big data sets in the future. More
detailed descriptions are shown in section 2.

Figure 1: Distribution shapes of data sets Imbalance, Animals, Ring

However, FC only can apply in idealized data sets, i.e., the border points of
$i$th cluster are far away from the border points of $j$th cluster. When FC
deal with the general data sets, it needs two parameters (such as FC-KNN[30]).
The algorithms which need parameters are difficult applied to the data sets
that obtain by automatic monitoring devices. In the following, a self-adaptive
noise identification method and a robust fission clustering algorithm are
proposed, and then they are combined together to propose a self-adaptive and
robust fission clustering (SARFC) algorithm. SARFC retains the performance of
FC and it can deal with all the density-based data sets without any
parameters.

The remainder of this article is organized as follows: Section 2 presents the
relevant knowledge of the fission clustering algorithm, heat diffusion density
factor and self-adaptive noise identification method of linear regression.
Section 3 proposes a robust fission clustering algorithm, a self-adaptive
noise identification method of maximal turning angle, and a self-adaptive and
robust fission clustering (SARFC) algorithm. Section 4 demonstrates our
algorithms by some numerical experiments of both simulated and real data sets
and image segmentation data through comparison with some state-of-the-art
clustering algorithms. Section 5 makes a conclusion.

##  2 Literature Research

###  2.1 Fission clustering (FC)

FC that proposed by Lu, Cheng and Mehmood[30] has several excellent feature,
such as the following F1, F2 and F3.

F1. To compare with many of current clustering algorithms, FC algorithm has
stronger coping ability for the extreme challenge data sets. As shown in
Figure 1, there are three data sets that are constructed for increasing the
clustering difficulty. The Imbalance data set magnifies the differences of two
clusters’ densities to increase the clustering difficulty. Many of current
density-based algorithms cannot discover the cluster with small density and
identify Imbalance as one category data. Many of current grid-based algorithms
classify Imbalance as more than two clusters data. The points in the cluster
with small density are divided in different grid cells and cannot be combined,
because the huge differences of two clusters’ densities are out of the
identification ability for grid-based algorithms with the same size of the
grid. The Animals data set has several profiles of animals, it is constructed
for increasing the difficulty of pattern recognition. Many of current model-
based algorithms are difficult to learn the correct pattern in Animals data
set and cannot obtain the correct clustering result. As shown in Figure 1 (c),
there is no single point can be considered as the geometrical centroid of the
annulus in the Ring data set, the points of the ”S” cluster interleave through
other two cluster. Many of current partitioning algorithms cannot find out the
central point for every cluster in data set Ring, and obtain the wrong
clustering result. As shown in Figure 1, FC algorithm can correctly cluster
these three data sets with extreme challenge.

F2. As shown in Figure 2, the processes above clustering result in Figure 2
are the clustering processes of FC algorithm, they are similar to atomic
fissions. FC algorithm can obtain $k$ clusters after $k-1$ times subset
partitions, its time complexity is $\mathcal{O}(k)$. The processes under
clustering result in Figure 2 are similar to atomic fusions. The algorithms
used this kind of clustering processes need to combine data point by point,
their time complexities are at least $\mathcal{O}(n)$. FC algorithm has a big
advantage in time complexity.

Figure 2: The clustering processes of two different directions

F3. As shown in Figure 2, FC clusters data via splitting data sets. When FC
faces the sets with billions of samples, it can use computing group to reduce
difficulty of data processing. FC divides data set into subsets and then sends
the subsets to the members of computing group to deal with subsets
respectively, such as, computer 1 deals with the subset $C_{1}$, computer 2
deals with the subset $C_{2}$, and so on. The algorithms which cluster data
with the pattern of atomic fusion cannot use computing group to optimize the
procedures of data processing for big data sets.

Although FC has many advantages to compare with other algorithms, but it also
has two disadvantages. (i) If a category in the data set has more than two
points which keep away from other points in this category, FC will consider
these outliers as another category, then one category is considered as two
categories (as shown in Figure 4). (ii) If FC deals with the nonidealized data
sets, it needs two parameters $k,t$[30]. The algorithms needed parameters are
not befitting for the automatic monitoring data. The section 3 of this article
presents an improved FC algorithm to solve the issues (i) and (ii).

###  2.2 Density factor of heat diffusion

Nonparametric density estimation is an important tool for statistical analysis
of data. It is used to evaluate skewness, multi-modality, summarizing Bayesian
posteriors, discriminant analysis, and classification[31, 32]. Nonparametric
approaches are more flexible for modelling of datasets and are not affected by
specification bias[32], in contrast to the classical approach[31]. Kernel
density estimation (KDE) is the most commonly used nonparametric density
estimation method[31]. Mehmood et al[9] utilized KDE to proposed a method to
calculate the density of every point in data set. Mehmood’s method can
estimate the density based on the nature of the data set without any
parameter, it overcomes the aforementioned issues that method is self-adaptive
for various data sets obtained by automatic monitoring devices.

KDE demonstrates a prominent performance for mining the local dense families
of sample set. The KDE for an identical and independent data set
$X=\\{x_{1},x_{2},\cdots,x_{n}\\}$ drawn with an unknown probability density
function (PDF) is

$\hat{f}_{h}(x;h)=\frac{1}{n}\sum\limits_{j=1}^{n}K_{h}(x-x_{j}),$ (1)

where $K_{h}(x-x_{j})$ is the kernel function and the Gaussian kernel
$K(x,x_{j};h)=\frac{1}{\sqrt{2\pi h}}e^{-\frac{(x-x_{j})^{2}}{2h}}$ is
normally used to estimate the density. The interpretation of KDE via heat
diffusion derives from the concept of the Weiner process, interpreted as the
probability distribution function at time $t$ of this process.

$\hat{f}_{t}(x;t)=\frac{1}{2}\sum\limits_{j=1}^{n}\frac{1}{\sqrt{2\pi
t}}e^{-\frac{(x-x_{j})^{2}}{2t}}.$ (2)

Equation (2) is an iterative process that is similar to equation (1) with
bandwidth $h$, it has an approximate solution as below[9],

$\hat{f}(x;t)\approx\sum\limits_{k=0}^{n-1}a_{k}e^{-k^{2}\pi^{2}t/2}\cos(k\pi
x),$ (3)

where $n$ is a large positive integer, and
$a_{k}=\frac{1}{n}\sum_{i=1}^{n}\cos(k\pi x_{i}),\ k=0,1,2,\cdots,n-1$.

Botev et al[31] applied the improved Sheather-Jones algorithm to obtain a
unique solution of the nonlinear equation. This has achieved by finding
adaptively the optimal bandwidth $t$ for KDE,

$t=\xi\gamma^{[l]}(t).$ (4)

The detailed description of the improved Sheather-Jones algorithm is shown
in[31].

The role of Equation (3) is applied to obtain an indicator (denoted as
$\rho_{i}$) for every object in the clustering sample set. A sample $x_{i}$
which has a larger value of $\rho_{i}$ is surrounded by many samples in the
data set.

###  2.3 Self-adaptive method of noise identification

Lu and Cheng et al[33] created an interesting and effective method to
distinguish the border objects and dense area objects. They researched the
density variations of two kinds of objects, and then utilized smooth sequence
of density sequence to make linear regression for noise identification (border
objects identification), as shown in follows.

Figure 3: The distinctions of border objects and dense area objects showed in
$V$. ((a) shows the border objects (black points) and dense area objects
(colourful points) in data set. (b) shows two distinctly different
distribution trends of $\\{v_{i}\\}=V$ which can be used to distinguish the
border objects and dense area objects.)

They utilized grey correlation analysis to obtain the density sequence
$\rho=\\{\rho_{i}\\}_{n}$ for data set $X$, and then smoothed the density
sequence $\rho$ to sequence $V$. Let
$\rho^{\prime}=\\{\rho_{i}^{\prime}:1\leqslant i\leqslant n\\}$ be the
descending sequence of $\rho=\\{\rho_{j}:1\leqslant j\leqslant n\\}$, and
$V=\\{v_{i}:5\leqslant i\leqslant n\\}$ be the mean smoothing sequence of
$\rho^{\prime}$, where
$v_{i}=(\rho_{i}^{\prime}+\rho_{i-1}^{\prime}+\rho_{i-2}^{\prime}+\rho_{i-3}^{\prime}+\rho_{i-4}^{\prime})/5$.

As shown in Fig. 3 (b), the elements of $V$ have two distinctly different
distribution trends. The elements in $V$ are divided into two parts
$V_{p}^{-}=\\{v_{i}:5\leqslant i\leqslant p\\}$ and
$V_{p}^{+}=\\{v_{i}:p<i\leqslant n\\}$. Then $(I_{p}^{-},V_{p}^{-})$ and
$(I_{p}^{+},V_{p}^{+})$ are used to obtain two models of linear regressions,
where $I_{p}^{-}=\\{i:5\leqslant i\leqslant p\\}$ and
$I_{p}^{+}=\\{i:p<i\leqslant n\\}$. Let $f_{1}$ and $f_{2}$ be the regression
equations for $(I_{p}^{-},V_{p}^{-})$ and $(I_{p}^{+},V_{p}^{+})$,
respectively. Then define

$\left\\{\begin{aligned} e_{1}(i)=f_{1}(i)-v_{i},\ 5\leqslant i\leqslant p,\\\
e_{2}(i)=f_{2}(i)-v_{i},\ p<i\leqslant n,\end{aligned}\right.$ (5)

and $R_{p}=\sum\limits_{5\leqslant i\leqslant
p}|e_{1}(i)|+\sum\limits_{p<i\leqslant n}|e_{2}(i)|$, $R_{p}$ is the sum of
absolute residuals (SoAR), depending on $p$.

As illustrated in Fig. 3, to pick $p_{1}$, $p_{2}$ and $p_{3}$ in different
positions for linear regressions and obtain three SoARs $R_{p_{1}}$,
$R_{p_{2}}$ and $R_{p_{3}}$, where, as shown, $R_{p_{2}}<R_{p_{1}}$ and
$R_{p_{2}}<R_{p_{3}}$.

The different positions $p_{i}$ can be used to generate a sequence of SoARs,
$\mathfrak{R}=\\{R_{i}:5+5\leqslant i\leqslant n-5\\}$ (5 boundary points left
for boundary regression). If for some $p$, $R_{p}=min(\mathfrak{R})$, then the
object $x_{j}$ is considered as a member of the dense subset for the case
$\rho_{j}\geqslant\rho^{\prime}_{p}$. Then, an idealized dense subset is
obtained for FC algorithm.

##  3 Proposed methods

###  3.1 Robust fission clustering algorithm

Lu, Cheng and Mehmood[30] utilized the following Definition 1 to divide set
into two subsets, utilized Theorem 1 to stop dividing the subsets when $k$
subsets was obtained. However, FC algorithm is not enough robust. As shown in
Figure 4 (a), $x_{1}$ and $x_{2}$ are two points that tinily leave the left
cluster in (a), FC considers $x_{1}$ and $x_{2}$ as an individual cluster.
$x_{3}$, $x_{4}$ and $x_{5}$ are also considered as an individual cluster in
(b). These points that tinily leave cluster cannot distinctly consider as
individual cluster. We will propose a robust fission clustering (RFC)
algorithm in the following.

Definition 1[30]. $f:X\times X\rightarrow R$ is a distance (similarity)
function, where $X$ is a sample set, and $R$ is the real number set. For all
$x_{k}\in X$, if $f(x_{0},x_{k})\not\in(f(x_{0},x_{i}),f(x_{0},x_{j}))$ (or
$f(x_{0},x_{k})\not\in(f(x_{0},x_{j}),f(x_{0},x_{i}))$), we call
$|f(x_{0},x_{i})-f(x_{0},x_{j})|$ a crack of $(X,f)$, where
$x_{0},x_{i},x_{j}\in X$.

Theorem 1[30]. If the distance (similarity) function $f$ satisfies triangle
inequality and $C\subset X$ has a $d_{0}(C)$-road, then $MC(C)\leqslant
d_{0}(C)$, where $MC(C)$ is the $MC$ of $(C,f)$.

Figure 4: Some points that tinily leave their clusters

To develop the RFC, some symbols and Definition 2 are necessary.

$\hat{x_{i}}^{(r)}$ is denoted as the $r$th near neighbor of $x_{i}$. Let
$d_{0}^{(m)}(C)=max\\{f(x_{i},\hat{x_{i}}^{(m)}):x_{i}\in C\subset X\\}$, and
$d_{0}^{(m)}=max\\{f(x_{i},\hat{x_{i}}^{(m)}):x_{i}\in X\\}$. $S(X)$ is
denoted as the similarity matrix with respect to similarity function $f$.
$S_{1}(X)$ is obtained by sorting every row of $S(X)$ with ascending. The
$i$th column of $S_{1}(X)$ is subtracted from the $(i+1)$th column of
$S_{1}(X)$ to obtain the $i$th column of the matrix $S_{2}(X)$.
$MC^{(r)}(X)=max\\{S_{2}^{(r)}(X)(i,j):i=1,2,\cdots,n,j=1,2,\cdots,n-r\\}$.

Definition 2. If there is a connected component in $C$ such that there is a
path from $x$ to $y$ for all $x,y\in C$ and the distance of random $r$
adjacent points $x_{p+1},x_{p+2},\cdots,x_{p+r}$ in this path is less than or
equal $d_{0}^{(r)}(C)$ (i.e. $f(x_{p},x_{p+r})\leqslant d_{0}^{(r)}(C)$). This
connected component is called $d_{0}^{(r)}(C)$-connected component (briefly
$d_{0}^{(r)}(C)$-CC).

Theorem 2. If the similarity function $f$ satisfies triangle inequality and
$C\subset X$ has a $d_{0}^{(r)}(C)$-CC, then $MC^{(r)}(C)\leqslant
d_{0}^{(r)}(C)$.

Proof. If there are $i$ and $p$ such that
$S_{2}^{(r)}(C)(i,p)=|f(x_{i},x_{p+r})-f(x_{i},x_{p})|>d_{0}^{(r)}(C)$ is
hold. Since $f$ satisfies triangle inequality and $C$ has a
$d_{0}^{(r)}(C)$-CC, then $x_{p+r}$ is not the $r$th near neighbor of $x_{p}$.
There must exists a $x_{t}$ such that $x_{p+r}$ is the $r$th near neighbor of
$x_{t}$ in a path of $d_{0}^{(r)}(C)$-CC.

The $x_{t}$ can be described in two cases: (i)
$f(x_{i},x_{t})\in(f(x_{i},x_{p+j}),f(x_{i},x_{p+j+1}))$; (ii)
$f(x_{i},x_{t})=f(x_{i},x_{p+j})$, where $0\leqslant j\leqslant r-1$.

For the case (i). $MC^{(r)}(C)$ is constituted by $r$ adjacent $MC(C)$. By
Definition 1, there is no $x_{t}\in C$ such that
$f(x_{i},x_{t})\in(f(x_{i},x_{p+j}),f(x_{i},x_{p+j+1}))$. A contradiction.

For case (ii). If $f(x_{i},x_{t})=f(x_{i},x_{p+j})$, then
$|f(x_{i},x_{p+r})-f(x_{i},x_{p})|=S_{1}(C)(i,p+r+1)-S_{1}(C)(i,p)$. Although
the mappings $f(x_{i},x_{t})$ and $f(x_{i},x_{p+j})$ are coincident, but there
are $r+1$ mappings in $(f(x_{i},x_{p}),f(x_{i},x_{p+r}))$, it is a
contradiction with the construction of $S_{2}(C)$.

Figure 5: The turning angles in curve $V$

Hence, $S_{2}^{(r)}(C)(i,p)\leqslant d_{0}^{(r)}(C)$ is hold for all $i$ and
$p$, i.e., $MC^{(r)}(C)\leqslant d_{0}^{(r)}(C)$. $\blacksquare$

Note that, FC in[30] is only a specific case of our RFC. When $r=1$, RFC and
FC are same.

To look back the issue that describes in Figure 4. When $r=2$, RFC can
clusters the data set in Figure 4 (a) as two categories, not three categories;
but also clusters the data set in Figure 4 (b) as three categories. When
$r=3$, RFC can clusters the data sets in Figure 4 (a) and (b) as two
categories. The distinctions of RFC and FC are described as following.  FC RFC
Stop dividing subsets with the condition
$d_{0}=max\\{f(x_{i},\hat{x_{i}}):x_{i}\in X\\}$ Stop dividing subsets with
the condition $d_{0}^{(r)}=max\\{f(x_{i},\hat{x_{i}}^{(r)}):x_{i}\in X\\}$ The
construction of $MC(C)$, $S_{2}(C)(:,k)=\newline
S_{1}(C)(:,k+1)-S_{1}(C)(:,k)$ The construction of $MC^{(r)}(C)$,
$S_{2}^{(r)}(C)(:,k)=\newline S_{1}(C)(:,k+r)-S_{1}(C)(:,k)$ Allow $r$
($r\geqslant 2$) points as an autocephalous category Avoid $m$ ($m\leqslant
r$) points to be an autocephalous category

Algorithm 1: RFC algorithm.

Input: Similarity (distance) matrix $S(X)$.

Output: Clusters of $X$.

1\. If $|X|\leqslant 1000$, $r\leftarrow 1$; if $1000<|X|\leqslant 2000$,
$r\leftarrow 2$; if $|X|>2000$, $r\leftarrow 3$.

2\. $d_{0}^{(r)}\leftarrow max\\{f(x_{i},\hat{x_{i}}^{(r)}):x_{i}\in X\\}$.

3\. $C_{1}\leftarrow X$ (initial value).

4\. While There is a subset $C_{i}$ such that $MC^{(r)}(C_{i})>d_{0}^{(r)}$ do

5\. repeat

6\. Pick the subset $C_{i}$ if $MC^{(r)}(C_{i})>d_{0}^{(r)}$.

7\. Sort every row of $S(C_{i})$ to obtain $S_{1}(C_{i})$.

8\. $S_{2}^{(r)}(C_{i})(:,k)\leftarrow S_{1}(C_{i})(:,k+r)-S_{1}(C_{i})(:,k)$,
$k=1,2,\cdots,n-r$.

9\. $MC^{(r)}(C_{i})\leftarrow max\\{S_{2}^{(r)}(C_{i})(k,j):k=1,2,\cdots,n,\
j=1,2,\cdots,n-r\\}$.

10\. Find out $x_{i},x_{j}$ and $x_{k}$ which generate the $MC^{(r)}(C_{i})$
($|f(x_{i},x_{j})-f(x_{i},x_{k})|=MC^{(r)}(C_{i})$).

11\. If $f(x_{i},x_{t})\leqslant min\\{f(x_{i},x_{j}),f(x_{i},x_{k})\\}$ then
$x_{t}\in C_{count}$; otherwise, $x_{t}\in C_{count+1}$.

12\. until $max\\{MC^{(r)}(C_{i})\\}\leqslant d_{0}^{(r)}$.

13\. end while

The details of the RFC algorithm are as shown in follows, where
$S_{k}(C)(:,i)$ is the $i$th column of $S_{k}(C)$.

###  3.2 An improved method for self-adaptive noise identification

$V$ is described in subsection 2.3. As shown in Figure 5, $l_{1}$ is the
straight line which passes two points $(2,v_{2})$ and $(3,v_{3})$, $l_{2}$ is
the straight line which passes two points $(3,v_{3})$ and $(4,v_{4})$. The
intersection angle of $l_{1}$ and $l_{2}$ is $\angle\alpha$. $\angle\alpha$ is
a turning angles which describes the veer of $v_{2}$, $v_{3}$ and $v_{4}$. As
shown in Figure 5, the points $v_{i}$ of $V$ have two distinctly different
distribution trends. The position which has a maximal turning angle can be
used to distinguish the points of two different distribution trends. We can
applied the maximal turning angle position to distinguish the dense area
points and border points.

If the sequence $V=\\{v_{i}\\}$ is smooth enough, the turning angles will be
very small (as shown in Figure 5). It is ensure that the turning angles are
fall in interval $[0,\pi/2]$. When $0\leqslant\angle\alpha_{i}\leqslant\pi/2$,
$\angle\alpha_{i}$ is positively related with $\tan\angle\alpha_{i}$. The
self-adaptive algorithm that utilize maximal turning angle to recognize noise
is shown in follows.

Figure 6: The main processes of the SARFC

Algorithm 2: Self-adaptive noise identification algorithm.

Input: Similarity (distance) matrix $S(X)$.

Output: Dense subset $C$.

1\. Apply equation (3) to obtain $\rho_{i}$ for every object.

2\. Sort the sequence $\rho=\\{\rho_{i}\\}$ in descending order to obtain
sequence $\rho^{\prime}=\\{\rho^{\prime}_{i}\\}$.

3\. Smooth the sequence $\rho^{\prime}$ to obtain sequence $V=\\{v_{i}\\}$,
where
$v_{i}=(\rho_{i}^{\prime}+\rho_{i-1}^{\prime}+\rho_{i-2}^{\prime}+\rho_{i-3}^{\prime}+\rho_{i-4}^{\prime})/5$.

4\. Utilize the method in Section 2.3 to obtain a position $p_{r}$.

5\. Obtain sequence $H=\\{h_{i}\\}$, where $h_{i}=v_{i+1}-v_{i}$.

6\. Obtain sequence $K=\\{k_{i}\\}$, where $k_{i}=h_{i}/\Delta T$, $\Delta T$
is the length of horizontal axis between $\rho_{i+1}$ and $\rho_{i}$, $\Delta
T=1$ in this article.

7\. Calculate tangent sequence $\\{\tan\alpha_{i}\\}$ for turning angles,
where $\tan\alpha_{i}=|(k_{i}-k_{i-1})/(1+k_{i}k_{i-1})|$.

8\. $p_{max}\leftarrow 1$.

9\. While $p_{max}<p_{r}$ do

10\. repeat

11\. $\tan\alpha_{p_{max}}\leftarrow 0$, and update tangent sequence
$\\{\tan\alpha_{i}\\}$.

12\. Find out the maximum of $\\{\tan\alpha_{i}\\}$ and the position $p_{max}$
of the maximum.

13\. until $p_{max}>p_{r}$.

14\. end while

15\. If $\rho_{j}\geqslant\rho^{\prime}_{p_{max}}$, the object $x_{j}$ is
considered as a member of the dense subset, then the dense subset $C$ is
obtained.

###  3.3 A self-adaptive and robust fission clustering (SARFC) algorithm

The main steps of the SARFC algorithm are as shown in follows.

Algorithm 3: SARFC algorithm.

Input: Similarity (distance) matrix $S(X)$.

Output: Clusters of $X$.

1\. Apply Algorithm 2 to obtain a dense subset $C$.

2\. Cluster the subset $C$ by using Algorithm 1.

3\. While $X-C\neq\emptyset$, do

4\. repeat

5\. Assign $x^{\prime}_{j}$ into the cluster which contains the object
$x^{\prime}_{i}$, if
$f(x^{\prime}_{i},x^{\prime}_{j})=min\\{f(x_{i},x_{j}):x_{i}\in C,x_{j}\in
X-C\\}$.

6\. Update $C$ and $X-C$, $C\leftarrow C\cup\\{x^{\prime}_{j}\\}$,
$X-C\leftarrow X-C-\\{x^{\prime}_{j}\\}$.

7\. until $X-C=\emptyset$.

8\. end while.

Figure 7: Distribution shapes of data sets SuPole and SquCir

The main processes of the SARFC can be displayed as Figure 6. Figure 6 (a) to
(b) shows the processes that Algorithm 2 obtains a dense subset $C$ from the
data set $X$. Figure 6 (b) to (c) shows the processes that Algorithm 1
clusters the dense subset $C$. Figure 6 (c) to (d) shows the processes that
assign the object in $X-C$ to its nearest category.

Equation (3) can be solved using fast Fourier transform and takes
$\mathscr{O}(nlogn)$ operations. The methods of Section 2.3 and 3.2 obtain
$p_{r}$ and $p_{max}$ take $\mathscr{O}(n)$ operations, respectively. The time
complexity of Algorithm 2 is $\mathscr{O}(nlogn)$. Algorithm 1 splits the set
$X$ (or dense subset $C$) into subsets $C_{1},C_{2},\cdots,C_{k}$. Since
$|C_{i}|\ll|X|$, data processing will become faster and faster, accompanied by
the dividing courses of subsets. $k$ clusters are obtained after $k-1$ times
of dividing subsets, then, the time complexity of Algorithm 1 is
$\mathscr{O}(1)$. Moreover, if $|C|=m$, then $|X-C|=n-m$, and the time
complexity of assigning border points is hence $\mathscr{O}(n-m)$. The time
complexity of the whole SARFC algorithm is $\mathscr{O}(nlogn)$.

##  4 Experiments

Some classical or state-of-the-art density-based algorithms will be selected
to compare with the proposed density-based algorithm: SARFC algorithm. Such
as, DP[8], ADPC[34], RNN-DBSCAN[35], MC-FNGLVQ[36] and UNIC[37]. The classical
density-based algorithm DP[8], and DP’s state-of-the-art improved algorithm
ADPC[34]. RNN-DBSCAN[35] is the self-adaptive improved algorithm of the
classical density-based algorithm DBSCAN[10] (not self-adaptive). Moreover, an
up to date intelligence algorithm MC-FNGLVQ[36] and a base on statistical
algorithm UNIC[37] are selected to be the added comparison algorithm.

###  4.1 Experiments on synthetic and real data sets

####  4.1.1 Descriptions of data sets

In this part, some frequently-used data sets obtained from different
clustering references are applied to test the algorithms. Such as R15[38],
D31[38], Aggregation (Agg)[39], A1[40] and S1[41] etc. And then two datasets,
SuPole and SquCir (Figure 7), are constructed for the supplementary tests. All
the simulation data are points of two-dimensional Euclidean space.

Three real-world data sets
Irisaaahttp://archive.ics.uci.edu/ml/datasets.php[42, 43], Seeds1[44] and Wine
1[45] are applied to make data experiments for algorithms comparisons. Data
set Iris contains three kinds of irises: Setosa, Versicolour and Virginica,
each kind has 50 samples. Data set Seed is obtained from three places: Kama,
Rosa and Canadian, 70 seeds from each place. Data set Wine are the results of
a chemical analysis of wines grown in the same region in Italy but derived
from three different cultivars. The description of the data sets are shown in
Table 1.

Data |  Instances |  Features |  Clusters
---|---|---|---
R15 |  600 |  2 |  15
D31 |  3100 |  2 |  31
Agg |  788 |  2 |  7
A1 |  3000 |  2 |  20
S1 |  5000 |  2 |  15
SuPole |  513 |  2 |  2
SquCir |  50000 |  2 |  2
Iris |  150 |  4 |  3
Seeds |  210 |  7 |  3
Wine |  178 |  13 |  3

Table 1: The simple description of data sets

####  4.1.2 Presentation of results

Table 2 presents the number of clusters estimated by various methods. Table 3
shows the clustering results when compared with other various methods.

To evaluate and compare the performance of the clustering methods, we apply
the evaluation metrics: Accuracy (Acc), F1-Score (F1), Adjusted Rand Index
(ARI)[46] and Normalized Mutual Information (NMI)[47] in our experiments to do
a comprehensive evaluation. The higher the value, the better the clustering
performance for all these measures. Compared with the best results of other
algorithms, our method has relative advantages of 7.4%, 5.5%, 11.5% and 7.2%
(Table 3) with respect to Accuracy, F1-Score, ARI and NMI for the Iris
dataset, respectively.

In summary, our method achieves better results with respect to the estimation
of cluster number, Acc, F1, ARI and NMI, compared comprehensively with other
methods.

Data |  DP |  ADPC |  RNN-DBSCAN |  MC-FNGLVQ |  UNIC |  SARFC
---|---|---|---|---|---|---
R15 |  15 |  15 |  15 |  9 |  15 |  15
D31 |  31 |  31 |  31 |  23 |  31 |  31
Agg |  7 |  7 |  7 |  8 |  7 |  7
A1 |  20 |  20 |  18 |  18 |  17 |  20
S1 |  15 |  15 |  8 |  11 |  9 |  15
SuPole |  5 |  3 |  15 |  3 |  5 |  2
SquCir |  9 |  2 |  2 |  15 |  6 |  2
Iris |  2 |  2 |  4 |  6 |  4 |  3
Seeds |  3 |  3 |  4 |  5 |  4 |  3
Wine |  7 |  6 |  5 |  6 |  5 |  4

Table 2: Number of clusters estimated by various algorithms

####  4.1.3 Comparisons

We conduct the Friedman test with the post-hoc Nemenyi test ($\alpha$=
0.10)[12] to examine whether the difference between any two clustering
algorithms is significant in terms of their average ranks. The difference
between two algorithms is significant if the gap between their ranks is larger
than CD. There is a line between two algorithms if the rank gap between them
is smaller than CD. This test shows that SARFC is significantly better than
the DP, RNN-DBSCAN, MC-FNGLVQ and UNIC. Except the SARFC, other algorithms are
not significantly better than each other. The SARFC is the best performer of
these algorithms, followed by the ADPC (as seen in Figure 8).

Data |  |  DP |  ADPC |  RNN-D BSCAN |  MC-FN GLVQ |  UNIC |  SARFC
---|---|---|---|---|---|---|---
R15 |  Acc |  99.2 |  99.2 |  98.5 |  65.7 |  66.9 |  99.3
|  $F_{1}$ |  99.2 |  99.2 |  97.1 |  63.3 |  64.3 |  99.4
|  ARI |  98.2 |  98.2 |  98.4 |  65.2 |  65.3 |  98.6
|  NMI |  98.6 |  98.6 |  98.2 |  61.3 |  65.2 |  98.9
D31 |  Acc |  96.8 |  96.8 |  91.6 |  53.2 |  68.0 |  96.8
|  $F_{1}$ |  96.8 |  96.8 |  95.3 |  49.9 |  68.6 |  96.8
|  ARI |  93.5 |  93.5 |  89.6 |  50.9 |  66.1 |  93.5
|  NMI |  95.8 |  95.8 |  90.3 |  52.1 |  70.0 |  95.8
Agg |  Acc |  99.9 |  99.9 |  99.3 |  87.3 |  73.3 |  100
|  $F_{1}$ |  99.8 |  99.8 |  98.9 |  83.1 |  71.7 |  100
|  ARI |  99.8 |  99.8 |  97.9 |  90.1 |  75.3 |  100
|  NMI |  99.6 |  99.6 |  99.9 |  89.3 |  73.6 |  100
A1 |  Acc |  94.3 |  94.3 |  88.5 |  85.2 |  78.0 |  97.2
|  $F_{1}$ |  93.3 |  93.3 |  89.4 |  83.7 |  80.5 |  97.2
|  ARI |  92.1 |  92.1 |  90.3 |  87.2 |  77.9 |  94.6
|  NMI |  93.1 |  93.1 |  89.4 |  89.0 |  81.0 |  96.2
S1 |  Acc |  92.6 |  92.6 |  80.9 |  89.3 |  85.5 |  99.3
|  $F_{1}$ |  93.3 |  93.3 |  82.2 |  88.1 |  81.7 |  99.3
|  ARI |  89.2 |  89.2 |  81.0 |  86.9 |  82.6 |  98.6
|  NMI |  94.5 |  94.5 |  84.1 |  89.6 |  83.0 |  99.0
SuPole |  Acc |  47.6 |  77.3 |  14.3 |  71.7 |  37.1 |  100
|  $F_{1}$ |  36.2 |  77.8 |  14.1 |  65.2 |  30.2 |  100
|  ARI |  26.8 |  27.1 |  29.6 |  39.9 |  20.5 |  100
|  NMI |  27.5 |  31.4 |  15.5 |  28.8 |  20.9 |  100
SquCir |  Acc |  46.7 |  58.2 |  57.3 |  15.7 |  26.1 |  98.6
|  $F_{1}$ |  43.3 |  55.6 |  56.2 |  11.3 |  23.3 |  97.8
|  ARI |  39.8 |  51.1 |  52.8 |  9.8 |  20.8 |  98.1
|  NMI |  37.2 |  52.8 |  53.6 |  10.2 |  21.1 |  98.3
Iris |  Acc |  66.7 |  66.7 |  83.3 |  69.3 |  80.7 |  90.7
|  $F_{1}$ |  57.1 |  57.1 |  86.2 |  72.1 |  84.5 |  91.7
|  ARI |  56.8 |  56.8 |  64.4 |  33.6 |  60.9 |  75.9
|  NMI |  73.4 |  73.4 |  73.4 |  41.6 |  72.3 |  80.6
Seeds |  Acc |  81.6 |  81.6 |  85.7 |  73.3 |  76.2 |  88.6
|  $F_{1}$ |  83.1 |  83.1 |  86.9 |  76.5 |  78.5 |  89.0
|  ARI |  80.1 |  80.1 |  61.7 |  37.8 |  43.0 |  70.3
|  NMI |  79.0 |  79.0 |  62.6 |  35.9 |  27.3 |  69.8
Wine |  Acc |  39.9 |  44.4 |  53.9 |  55.1 |  56.2 |  61.2
|  $F_{1}$ |  53.9 |  58.2 |  61.2 |  64.9 |  63.8 |  58.7
|  ARI |  21.7 |  23.9 |  21.6 |  28.7 |  23.3 |  31.6
|  NMI |  19.8 |  23.4 |  14.0 |  16.5 |  13.2 |  40.3

Table 3: The results’ comparison for different methods (%)

Figure 8: Critical difference (CD) diagram of the post-hoc Nemenyi test (at a
significance level $\alpha$= 0.10)

##  5 Conclusion

Benefit from the automatic monitoring technique, we will obtain more and more
automatic monitoring data sets. Except the sample size and data
dimensionality, we do not know to much information for this kind of data. Many
current clustering algorithms need to input some parameters, but we do not
know how to set the parameters for the automatic monitoring data. Hence, self-
adaptive clustering algorithms are the future needs. In this article, a self-
adaptive robust fission clustering (SARFC) algorithm is proposed. And then
both the simulation and real datasets are applied to test the performance and
effectiveness of the proposed method. Our proposed algorithm is also compared
with several frequently-used clustering algorithms. The experiments indicate
that our method achieves better results, in terms of the evaluation metrics
(Table 3) and the estimated number of clusters (Table 2), than the other
methods under comparison.

REFERENCES

## References

  * [1] CP Kappe, M Böttinger, H Leitte. Analysis of Decadal Climate Predictions with User guided Hierarchical Ensemble Clustering. Computer Graphics Forum, 2019, 38(3): 505-515.
  * [2] T Wang, J Zhang, K Huang. Generalized gene co-expression analysis via subspace clustering using low-rank representation. BMC Bioinformatics, 2019, 20(7): 17-27.
  * [3] A Pessia, J Corander. Kpax3: Bayesian bi-clustering of large sequence datasets. Bioinformatics, 2018, 34(12): 2132-2133.
  * [4] JD Hamilton. A new approach to the economic analysis of nonstationary time series and the business cycle. Econometrica, 1989, 57(2): 357-384.
  * [5] G Leibon, S Pauls, D Rockmore, R Savell. Topological structures in the equities market network. Proceedings of the national academy of sciences of the united states of america, 2008, 105(52): 20589-20594.
  * [6] S Galbraith, JA Daniel, B. Vissel. A study of clustered data and approaches to its analysis. Journal of Neuroscience, 2010, 30(32): 10601-10608.
  * [7] HK Aljobouri, et al. Clustering fMRI data with a robust unsupervised learning algorithm for neuroscience data mining. Journal of Neuroscience Methods, 2018, 299: 45-54.
  * [8] A Rodriguez, A Laio. Clustering by fast search and find of density peaks. Science, 2014, 344(6191): 1492-1496.
  * [9] R Mehmood et al. Clustering by fast search and find of density peaks via heat diffusion. Neurocomputing, 2016, 208: 210-217.
  * [10] M Ester, HP Kriegel, J Sander, X Xu. A density-based algorithm for discovering clusters in large spatial databases with noise. Data Mining and Knowledge Discovery, 1996, 96(34): 226-231.
  * [11] Y Chen, S Tang, et al. A fast clustering algorithm based on pruning unnecessary distance computations in DBSCAN for high-dimensional data. Pattern Recognition, 2018, 83: 375-387.
  * [12] Y Zhu, KM Ting, MJ Carman. Grouping points by shared subspaces for effective subspace clustering. Pattern Recognition, 2018, 83: 230-244.
  * [13] J Jiang, D Hao, et al. GDPC: Gravitation-based Density Peaks Clustering algorithm. Physica A, 2018, 502: 345-355.
  * [14] R Agrawal, JE Gehrke, D Gunopulos, P Raghavan. Automatic subspace clustering of high dimensional data for data mining applications. in Proc. ACM SIGMOD Int. Conf. Manage. Data (SIGMOD), Seattle, WA, USA, 1998: 94-105.
  * [15] C Deng et al. Gridwave: a grid-based clustering algorithm for market transaction data based on spatial-temporal density-waves and synchronization. Multimedia Tools and Applications, 2018, 77(21): 1-15.
  * [16] G Sheikholeslami, S Chatterjee, A Zhang. WaveCluster: a wavelet-based clustering approach for spatial data in very large data bases. VLDB Journal, 2000, 8: 289-304.
  * [17] K Murphy, TB Murphy. Gaussian parsimonious clustering models with covariates and a noise component. Advances in Data Analysis and Classification, 2019: 1-33.
  * [18] A O’Hagan et al. Investigation of Parameter Uncertainty in Clustering Using a Gaussian Mixture Model Via Jackknife, Bootstrap and Weighted Likelihood Bootstrap. Computational Statistics, 2019, 34(4): 1779-1813.
  * [19] T Chen, NL Zhang, T Liu, KM Poon, Y Wang. Model-based multidimensional clustering of categorical data. Artificial Intelligence, 2012, 176(1): 2246-2269.
  * [20] J MacQueen. Some Methods for Classification and Analysis of Multivariate Observations. in Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, LM Le Cam, J Neyman, Eds. (Univ. California Press, Berkeley, CA), 1967, 1: 281-297.
  * [21] R Zhang, X Li, et al. Deep Fuzzy K-Means with Adaptive Loss and Entropy Regularization. IEEE Transactions on Fuzzy Systems, 2019, DOI:10.1109/TFUZZ.2019.2945232.
  * [22] R Zhang et al. Joint Learning of Fuzzy k-Means and Nonnegative Spectral Clustering With Side Information. IEEE Transactions on Image Processing, 2019, 28(5): 2152-2162.
  * [23] DW Choi, CW Chung. A K-partitioning algorithm for clustering large-scale spatio-textual data. Information Systems, 2017, 64: 1-11.
  * [24] K Lahari, MR Murty, SC Satapathy. Partition based clustering using genetic algorithm and teaching learning based optimization: performance analysis. Advances in Intelligent Systems and Computing, 2015, 338: 191-200.
  * [25] J Wang, W Zheng, Y Qian, J Liang. A Seed Expansion Graph Clustering Method for Protein Complexes Detection in Protein Interaction Networks. Molecules, 2017, 22(12): 1-19.
  * [26] P Pipenbacher, A Schliep, S Schneckener, et al. ProClust: improved clustering of protein sequences with an extended graph-based approach. Bioinformatics, 2002, 18(2): 182-191.
  * [27] VV Vu, HQ Du. Graph-based Clustering with Background Knowledge. SoICT 2017 Proceedings of the Eighth International Symposium on Information and Communication Technology, NY, USA, 2017: 167-172.
  * [28] T Zhang, R Ramakrishnan, M Livny. BIRCH: An efficient data clustering method for very large databases. in Proc. ACM SIGMOD Int. Conf. Manage. Data (SIGMOD), Montreal, QC, Canada, 1996, 25(2): 103-114.
  * [29] G Karypis, EH Han, V Kumar. CHAMELEON: A hierarchical clustering algorithm using dynamic modeling. IEEE Computer, 1999, 32(8): 68-75.
  * [30] SZ Lu, LS Cheng, R Mehmood. Clustering by Using the Way of Atomic Fission. IEEE Access, 2020, 8: 70997-71007.
  * [31] ZI Botev, JF Grotowski, DP Kroese. Kernel density estimation via diffusion. Annals of Statistics, 2010, 38(5): 2916-2957.
  * [32] LE Lehmann. Model specification: the views of Fisher and Neyman, and later developments, Statistical Science, 1990, 5(2): 160-168.
  * [33] SZ Lu, LS Cheng, ZD Lu, BA Khan. A Self-adaptive Grey DBSCAN Clustering. 2020. https://arxiv.org/abs/1912.11477.
  * [34] T Liu, H Li, X Zhao. Clustering by Search in Descending Order and Automatic Find of Density Peaks. IEEE Access, 2019, 7: 133772-133780.
  * [35] A Bryant, K Cios. RNN-DBSCAN: A Density-Based Clustering Algorithm Using Reverse Nearest Neighbor Density Estimates[J]. IEEE Transactions on Knowledge and Data Engineering, 2018, 30(6): 1109-1121.
  * [36] MA Másum. Intelligent Clustering and Dynamic Incremental Learning to Generate Multi-Codebook Fuzzy Neural Network for Multi-Modal Data Classification[J]. Symmetry, 2020, 12(4): 679.
  * [37] N Leopold, O Rose. UNIC: A fast nonparametric clustering[J]. Pattern Recognition, 2020, 100: 107-117.
  * [38] CJ Veenman, MJT Reinders, E Backer. A maximum variance cluster algorithm. Ieee Transactions On Pattern Analysis And Machine Intelligence, 2002, 24(9): 1273-1280.
  * [39] A Gionis, H Mannila, P Tsaparas. Clustering aggregation. ACM Transactions on Knowledge Discovery from Data, 2007, 1(1): 1-30.
  * [40] K Ismo, P Franti. Dynamic local search for clustering with unknown number of clusters. in: Proceedings of International Conference on Pattern Recogn., 2002, 2(16): 240-243.
  * [41] P Franti, O Virmajoki. Iterative shrinking method for clustering problems. Pattern Recognition, 2006, 39(5): 761-775.
  * [42] RA Fisher. The use of multiple measurements in taxonomic problems. Annual Eugenics, 1936, 7(2): 179-188.
  * [43] F Huang, X Li, S Zhang, J Zhang. Harmonious Genetic Clustering. IEEE Transactions on Cybernetics, 2018, 48(1): 199-214.
  * [44] M Charytanowicz, J Niewczas, et al. A Complete Gradient Clustering Algorithm for Features Analysis of X-ray Images. in: Information Technologies in Biomedicine, Ewa Pietka, Jacek Kawa (eds.), Springer-Verlag, Berlin-Heidelberg, 2010, 15-24.
  * [45] P Zhong, M Fukushima. Regularized nonsmooth Newton method for multi-class support vector machines. Optimization Methods and Software, 2007, 22(1): 225-236.
  * [46] L Du, Y Pan, X Luo. Robust spectral clustering via matrix aggregation. IEEE Access, 2018, 6: 53661-53670.
  * [47] S Abbasi, S Nejatian, et al. Clustering ensemble selection considering quality and diversity. Artificial Intelligence Review, 2019, 52: 1311-1340.