
\chapter{High Dimensional Inference\label{chapL:Confirmatory-Analysis}}

\index{confirmatory analysis|(}Confirmatory analysis is a nearly
completely dominant part of classical statistics, as typically taught
in modern courses. Much of this is based on fitting parametric probability
distributions to data, and basing inference such as hypothesis tests
and confidence intervals, on such distributions. As noted in \citet{Marron2017BigData}
and \citet{CarmichaelMarron2017DataSciStat}, this approach lies at
the roots of the \emph{scientific method}\index{scientific method},
which has provided great benefits to science over the years.

So far this approach has not been well developed in OODA contexts,
often because in many OODA data settings such as the manifold data
objects of Chapter \ref{chapG:Manifold-Data-Analysis} and the \index{tree}tree
structured \index{data objects!tree structured}data objects in Chapter
\ref{chapI:Tree-Structured-Data}, suitable probability models are
not yet in common use. However as shown using the Overlapping Classes\index{Overlapping Classes data|(}
1000 dimensional Gaussian data set in Figures \ref{fig:HiDimToyPCA }
and \ref{fig:HiDimToyDWD}, statistical confirmation is a critically
important task, to ensure that non-spurious and reproducible structure
has been discovered. Two main approaches to this are discussed in
this chapter. The DiProPerm\index{DiProPerm|(} hypothesis test, for
testing the difference between two \emph{previously labeled} subgroups
is discussed in Section \ref{secL1:DiProPerm}. While it is tempting
to use DiProPerm to assess the significance of subgroups found by
clustering\index{clustering} methods, it is seen in Section \ref{secL2:SigClust}
that this can be quite inappropriate. A test for significance of clustering,
that is appropriate in the challenging high dimensional context, called
SigClust, is described in Section \ref{secL2:SigClust}.

\section{DiProPerm - Two Sample Testing\label{secL1:DiProPerm}}

A sensible way of understanding statistical significance of features
observed in a visualization\index{visualization|(} is to use a hypothesis
test directly related to that visualization. This goal is achieved
by the DiProPerm (for DIrection-PROjection-PERMutation) hypothesis
test proposed by \citet{Wei2016DiProPerm}, discussed in Section \ref{secB3:Confirmatory-Analysis}.
For example the PCA scatterplot\index{scores!scatterplot} analysis
of the $d=20,000$ dimensional Two Class Gaussian data\index{Two Class Gaussian data}
set in Figure \ref{figF.C2:ToyPCA} shows two distinct classes, which
are not at all visually apparent in the heat map view of the same
data shown in Figure \ref{figF.C1:ToyHeatMap}. Statistical significance
of this visual difference is investigated using DiProPerm in Figure
\ref{figL.A:ToyPCADiProPerm}. The upper left panel shows the first
step of DiProPerm: visualization based on the MD (recall the acronym
for Mean Difference\index{classifier!Mean Difference}\index{Mean Difference|(}
from Section \ref{secJ-A:ClassicalMethods}) direction and projection
of the data onto it, which clearly highlights this class difference.
The latter is quantified by subtracting the means of the projected
classes, which gives about 23.1, indicated using a vertical green
line in the top right panel. Statistical significance of this class
difference is assessed using a permutation test.

The permuted null distribution comes from randomly relabeling the
data, and for each of those recomputing both MD and the projections
followed by recomputation of the projected class mean differences.
Two realizations of this process are shown in the bottom panels of
Figure \ref{figL.A:ToyPCADiProPerm}. The colors and symbols indicate
how the permutation works. In the top left panel the magenta data
are all indicated with circles, with x signs for the green cases.
In the bottom panels the same symbols are used for each data point,
but about half of the magenta circles have had their label (i.e. class
color) switched to green and similarly for the x signs. Note that
the horizontal scales in the bottom panels are different from the
top left, and in fact the two permuted group means, 20.04 and 20.02
are substantially smaller than the 23.1 observed for the original
data. Those are two of the black dots shown in the \emph{permutation
plot\index{permutation plot}} in the upper right panel of Figure
\ref{figL.A:ToyPCADiProPerm}, which are projected class differences
from 1000 independent permutations. Note the black dots reflect permutation
differences that are all much smaller than the green line showing
this difference is strongly statistically significant. A quantification
of this is the \emph{empirical p-value}\index{empirical p-value|(},
which is the proportion of the black dots to the right of the green
line, which in this case is 0. In situations where it is desired to
compare the significance of several DiProPerm results, empirical p-values
are not useful when both are 0, e.g. as for the Pan Cancer data\index{Pan Cancer data|(}
shown in Figure \ref{figL.C:TCGADiProPerm}. This issue is usefully
addressed with the \emph{Z-score}\index{Z-score|(}, an alternate
measure of significance. The Z-score summarizes the distribution of
the black dots using its mean and standard deviation, and measures
significance of the DiProPerm test as the number of standard deviations
that the green line is above the mean, in this case about 30.5.

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterL/OODAbookChpLFigAToyPCADiProPerm}

\caption{DiProPerm analysis of the significance of the difference between the
means of the magenta and green classes (same symbols and colors) of
the Two Class Gaussian data shown in Figure \ref{figF.C2:ToyPCA}.
Upper left panel shows projections onto the MD direction, with class
contrast summarized by the difference of the class means (green line
in the permutation plot, top right panel). Bottom panels show two
realizations of permuted class labels, whose class mean differences
are summarized in the top right panel. Shows class difference is strongly
significant, despite not being visible in the heat map view of Figure
\ref{figF.C1:ToyHeatMap}. \label{figL.A:ToyPCADiProPerm}}
\end{figure}
The choice of MD direction made in Figure \ref{figL.A:ToyPCADiProPerm}
results in a hypothesis test which is equivalent to a standard permutation
mean test, i.e. projection is not needed except for the visualization.
However, as noted in Section \ref{secJ-D:DWD} DWD\index{Distance Weighted Discrimination|(}
provides a direction which can be better in a number of senses, including
data visualization, in which case the projection step of DiProPerm
is essential beyond just visualization. The potential benefit of using
the DWD direction in DiProPerm is illustrated in Figure \ref{figL.C:TCGADiProPerm}.

Figure \ref{figL.B:ToyDWDDiProPerm} shows a DWD DiProPerm analysis
of the Overlapping Classes data from Figures \ref{fig:HiDimToyPCA }
and \ref{fig:HiDimToyDWD}. Recall the classes shown there were both
sampled from the standard normal distribution, so there should be
no significant difference, yet the DWD projection suggested an apparent
strong difference in the upper left panel of Figure \ref{fig:HiDimToyDWD}
(repeated in the left panel of Figure \ref{figL.B:ToyDWDDiProPerm}).
The realizations of the permutations (analogs of the bottom panels
of Figure \ref{figL.A:ToyPCADiProPerm}) are quite similar to the
left panel and are thus not shown here (and in other DiProPerm examples
in this section) to save space. The right panel shows the permutation
plot with the mean differences of the projections onto the DWD direction
for $1000$ permutations. Note that this time the green line lies
in the middle of the null distribution of black dots, as quantified
by an empirical p-value of 0.82, showing this projected difference
is clearly not significant (again despite the contrary visual impression
in the left panel). As discussed in Section \ref{secB3:Confirmatory-Analysis},
this reveals how important it is for exploratory data analysis\index{exploratory analysis}
to be complemented by confirmatory analysis. 

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterL/OODAbookChpLFigBToyDWDDiProPerm}

\caption{DiProPerm analysis of the Overlapping Classes data\index{Overlapping Classes data|)}
shown in Figures \ref{fig:HiDimToyPCA } and \ref{fig:HiDimToyDWD},
using the DWD projection. Reveals that the apparent large gap between
the red and blue classes is not statistically significant. i.e. is
an artifact of natural sampling variation. \label{figL.B:ToyDWDDiProPerm}}
\end{figure}
Given the goal of connecting inference with visualization by DiProPerm,
there are several reasonable choices of direction vector. E.g. the
MD was used in Figure \ref{figL.A:ToyPCADiProPerm}, while Figure
\ref{figL.B:ToyDWDDiProPerm} was based on DWD. An advantage of the
MD is that it is faster to compute. However the sharper class distinctions
available from DWD may be expected to give stronger inference, when
there actually is a difference between the groups. An example of this
property is shown in Figure \ref{figL.C:TCGADiProPerm}, using the
same log transformed Pan Cancer data (same colors and symbols) as
in Figure \ref{figD-L:TCGA-Log2vsRaw}. The top left panel shows the
MD projection, which does separate the classes, but the bottom left
panel shows that projecting on the DWD direction gives a cleaner separation.
This can be seen from carefully studying the kernel density estimates
in the valley between the cancer types. This example also demonstrates
the value of the Z-score over the empirical p-value. In particular
both of the latter are 0, yet the much larger Z-score (about 28.8
for DWD versus 15.9 for MD) demonstrates that DWD gives much more
powerful inference in this case. An even more dramatic example of
DWD giving stronger DiProPerm results appears in the analysis of the
Lung Cancer data\index{Lung Cancer data|(} discussed around Figure
\ref{figL.L:ScatPlotsLungC}.

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterL/OODAbookChpLFigCDiProPermMD}

\bigskip{}

\includegraphics[scale=0.8]{ChapterL/OODAbookChpLFigCDiProPermDWD}

\caption{Comparison of MD and DWD directions in DiProPerm, using log Pan Cancer
data from Figure \ref{figD-L:TCGA-Log2vsRaw}. The larger Z-score
in the permutation plot (lower right panel) shows that sharper class
separation from DWD results in a more powerful hypothesis test. \label{figL.C:TCGADiProPerm}}
\end{figure}
Also important is the statistic used to summarize the difference between
the projected means, shown in green text in the projection plots and
as a green vertical line in the permutation plots. A perhaps natural
choice is the 2 sample t statistic. This is classical in conventional
hypothesis testing because it incorporates sample variability into
the statistic. However, such incorporation is not needed here, as
the permutation already deals with this variation. \citet{Wei2016DiProPerm}
showed that using the t statistic actually gives some instability
causing a potential loss of power, so the simple difference of the
means is recommended instead. There are other reasonable choices for
quantifying differences in 1-d distributions, such as the AUC from
the ROC curve, discussed in Section \ref{secD2:Preprocess-Trans}. 

Several DiProPerm analyses of the Drug Discovery data\index{Drug Discovery data}
studied in Chapter \ref{chapD:Preprocessing} are summarized in Table
\ref{tabL.A:DrugDiProPermZscores}. Recall from the full raw data
PCA plot in Figure \ref{figD-A:DrugPCARaw} a major goal for this
data set was to contrast actives (red) from in-actives (blue), so
these are contrasted using DiProPerm hypothesis tests. All empirical
p-values are 0, so again meaningful comparisons come from the Z-scores
shown in the right column of the table, with a relevant figure number
shown for quick reference in the middle. The second row shows the
initial result of removing all variables with any missings, which
were coded as -999 (recall the finding of these was highlighted in
Figure \ref{figD-D:DrugMargDistMeanSmall}). While the PCA for that
data set has not been shown here, as it is essentially the same as
Figure \ref{figD-A:DrugPCARaw}, this shows that removal of these
variables actually has a meaningful effect.

\begin{table}[H]
\begin{tabular}{|c|c|c|}
\hline 
Data & Figure Number & Z-score\tabularnewline
\hline 
\hline 
Raw & \ref{figD-A:DrugPCARaw} & 10.4\tabularnewline
\hline 
No Missings & \ref{figD-D:DrugMargDistMeanSmall} & 11.6\tabularnewline
\hline 
Binary Only & \ref{figD-J:DrugVDBinaryPCA} & 14.6\tabularnewline
\hline 
Standardized & \ref{figD-K:DrugVDnonBinaryPCA} & 17.3\tabularnewline
\hline 
Auto Log Trans'd & - & 17.8\tabularnewline
\hline 
\end{tabular}

\caption{Summaries of DiProPerm Z-scores for various versions of the Drug Discovery
data studied in Sections \ref{secD1:Preprocess-MargDists} - \ref{secD2:Preprocess-Trans}.
Shows improved separation of classes following improved preprocessing\index{preprocessing}.
\label{tabL.A:DrugDiProPermZscores}}
\end{table}
The third row of Table \ref{tabL.A:DrugDiProPermZscores} summarizes
the application of DiProPerm to only the binary data\index{binary data},
whose PCA appears in Figure \ref{figD-J:DrugVDBinaryPCA}. Perhaps
surprisingly this filtering of the data contains more active - inactive
information than the raw data, as evidenced by the larger Z-score.
The benefits of standardizing the non-binary part of the data are
clear from the larger Z-score in the fourth row (PCA in Figure \ref{figD-K:DrugVDnonBinaryPCA}).
Some (fairly marginal) improvement comes from application of the automatic
shifted log transformation\index{shifted log transformation}\index{transformation!log}
discussed in Section \ref{secD2:Preprocess-Trans}.

DiProPerm was motivated by the desire to quantify the reproducibility
of lessons learned from data visualizations. That goal is not standard
in the area of high dimensional hypothesis testing, where the focus
is usually on optimization of test power. Important references in
that area include \citet{BaiSarana1996}, \citet{SrivastavaDu2008}
and \citet{ChenQin2010}. A method quite similar to DiProPerm was
proposed by \citet{Ghosh2016Hyperplanes}. The statistical power of
DiProPerm has been substantially improved in strong signal cases by
using \emph{\index{balanced permutations}balanced permutations} in
\citet{Yang2021DiProPerm}.

\section{Statistical Significance in Clustering\label{secL2:SigClust}}

Another\index{clustering|(} important type of confirmatory analysis
is determining when clusters discovered using say a method from Chapter
\ref{chapK:Clustering} are ``really there'', i.e. represent important
and reproducible underlying population structure as opposed to being
spurious sampling artifacts. One might attempt to use a mean based
hypothesis test, such as the DiProPerm method of Section \ref{secL1:DiProPerm}
to the discovered classes for this purpose. However, as pointed out
by Andrew Nobel, this is inappropriate as illustrated in Figure \ref{figL.D:NotDiProPerm}.
That analysis is based on the $n=100$ two dimensional standard normal
data points displayed as a scatterplot in the lower right panel. The
colors and symbols used there show the result of an application of
2-means clustering, which essentially splits the data in half. These
can be considered to be two classes whose mean difference can be tested
with DiProPerm. Projection onto the MD direction is shown in the upper
left panel. One realization of the permutation test is shown in the
lower left panel, revealing a difference of projected class means
that is \emph{much} smaller than in the upper left panel. This is
born out in the 1000 replication summary (upper right panel) where
the Z-score of about 11.5 shows the red and blue class means are clearly
different. However, because the data are a single standard normal
cluster, it does not make sense to conclude that there are two distinct
clusters here. The lesson is that \index{DiProPerm|)}DiProPerm is
about the difference in means of \emph{given classes}, which is \emph{different}
from determination of clusters, so another approach is needed for
clustering applications. 

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterL/OODAbookChpLFigDNotDiProPerm}

\caption{Simple 2-means clustering example to illustrate how DiProPerm applied
in clustering contexts may inappropriately indicate clusters. Raw
data in the bottom right are clearly a single (in fact Gaussian) cluster,
yet the DiProPerm (applied to the two colored groups in the lower
right) result shown in the permutation plot (top right) is clearly
strongly significant. \label{figL.D:NotDiProPerm}}
\end{figure}
The first major hurdle to construction of a test for statistically
significant clusters is to define the notion of \emph{cluster}. \citet{Tibshirani2001GapStat}
suggest thinking about circular uniform distributions, as a conservative
approach, but there are not natural probability mechanisms which generate
that. A reasonable null hypothesis for naturally generated clusters
is the Gaussian distribution, proposed by \citet{Sarle1993MODECLUS}.
That idea underlies the \emph{SigClust\index{SigClust|(}} hypothesis
test proposed by \citet{Liu2008SigClust} and substantially improved
by \citet{Huang2015SigClust}. 

SigClust assesses how well a single Gaussian distribution fits the
data set, by measuring goodness of fit using the 2-means $CI$ (Cluster
Index, defined at (\ref{eqK-A:ClusterIndex})). A p-value (or Z-score)
can be estimated by comparing the $CI$ from the data with a null
population generated by simulating data sets of the same size $n$
from a fitted Gaussian distribution. It is worth noting there are
two useful SigClust tasks, which are handled in a similar way.
\begin{enumerate}
\item \emph{Confirmatory SigClust.\index{Sig-Clust!confimatory}} Testing
whether a \emph{given} clustering (i.e. set of labels) is significantly
different from a fitted Gaussian distribution. Such labels could come
from some extrinsic information, or just from a visual impression,
when there is interest in whether the clusters are ``really there''.
Here the focus is on whether or not the given $CI$ is smaller than
those of the null Gaussian distribution.
\item \emph{Exploratory\index{Exploratory analysis} SigClust}.\index{Sig-Clust!exploratory}
Investigating whether a data set should be considered to be more than
one cluster. In this case comparison is done between the 2-means $CI$
(recall the minimizer over all potential clusterings) for the data
and the null Gaussian $CI$ distribution.
\end{enumerate}
Such a p-value from this second approach, for the data in Figure \ref{figL.D:NotDiProPerm},
appears in the lower right panel. This is clearly insignificant as
expected from underlying Gaussian data, so SigClust gives an appropriate
answer of one cluster in this case.

The next major hurdle is estimation of an appropriate Gaussian null
distribution. In most situations it will be inappropriate to use a
spherically symmetric Gaussian distribution. Thinking about the \index{Four Cluster data|(}Four
Cluster data set in the left panel of Figure \ref{figK-A:ClusterToyData},
the blue circles and the cyan squares each seem to constitute separate
clusters. But what about the collection of green plus signs? Using
the definition in terms of coming from a single Gaussian this is taken
here to be a single cluster, although others might prefer a different
choice. Also thinking about the Gaussian based definition of cluster,
the two red points are regarded as their own small cluster, as single
sufficiently extreme outliers are also regarded as singleton clusters. 

Various subsets of the two dimensional Four Cluster data in the left
panel of Figure \ref{figK-A:ClusterToyData} are used to build intuition
for SigClust. The straightforward case of just the blue circles and
cyan squares is considered in Figure \ref{figL.E:SigClustTwoClust}.
Note that just those two clusters are in full color in the left panel,
while the other points are gray because they are not used in this
example. Also shown in the left panel is an elliptical contour of
the fit Gaussian density together with a sample of the same size drawn
from that density (black stars). Finally the $CI$ using the blue
and cyan labeling (hence using SigClust in confirmatory mode) is shown
near the top ($CI=0.085$), and also is the horizontal coordinate
of the green line in the right panel. As expected, that is much smaller
than the $CI$ of the black points shown near the bottom ($CI=0.462$).
The cluster labeling used for the black points comes from applying
two means clustering, which aims to minimize the $CI$ over all possible
class labelings, as discussed in Section \ref{secK-A:KMeansClust}.
The comparison of the data $CI$ with the null distribution is shown
in the right panel. The black dots are 1000 realizations of the $CI$
from data drawn from the Gaussian fit. Generally in SigClust, a real
data $CI$ (shown in green) that is \emph{smaller} than the simulated
null $CI$ population provides strong evidence of more than one cluster
(in contrast to DiProPerm where larger values are more significant).
In this example, the green line is much smaller than the simulated
null distribution, indicating strong significance. These SigClust
results are reported in a similar fashion as in Section \ref{secL1:DiProPerm},
where the empirical p-value is the proportion of black dots to the
left of the green line, and the Z-score represents the number of standard
deviations the green line is to the left of the black dot null population
(again being useful for comparing results with 0 p-values).

\begin{figure}[H]
\includegraphics[scale=0.6]{ChapterL/OODAbookChpLFigESigClustTwoClustLeft}\qquad{}\includegraphics[scale=0.5]{ChapterL/OODAbookChpLFigESigClustTwoClustRight}

\caption{SigClust analysis, using the given cluster labels of the blue and
cyan clusters in the Four Cluster data. Shows a strongly significant
result. \label{figL.E:SigClustTwoClust}}
\end{figure}
Potential splitting of the long green cluster is studied in Figure
\ref{figL.F:SigClustSplitLong}, using a similar graphic. This time
only the green points are not grayed out, and are used in the calculation.
Here SigClust is used in the exploratory mode described above, so
no input clustering is used, and instead SigClust only aims to assess
whether or not the long green cluster should be split. Because there
is no input clustering, no $CI$ is shown in the left panel. However,
everything else is similar to before, with an overlay of the fit density
contour, and a corresponding simulated Gaussian data set of the same
size. This time the distributional summary plot in the right panel
shows the data $CI$ (green line) is well in the middle of the null
population of black dots. Hence there is no evidence for more than
one cluster in this case, also as expected.

\begin{figure}[H]
\includegraphics[scale=0.6]{ChapterL/OODAbookChpLFigFSigClustSplitLongLeft}\qquad{}\includegraphics[scale=0.5]{ChapterL/OODAbookChpLFigFSigClustSplitLongRight}

\caption{SigClust investigation of whether to split the long green cluster
in the Four Cluster data. Large p-value indicates a statistically
significant split has not been found, i.e. the green plus signs are
a single cluster. \label{figL.F:SigClustSplitLong}}
\end{figure}
A much more challenging case is the union of the long green cluster
and the red outliers, shown in Figure \ref{figL.G:SigClustLongOutlier}.
The format is the same as the two figures above, but now the surprising
feature is that the data $CI$ is actually larger than the simulated
null population. This is caused by the severe imbalance between the
given cluster sizes. The 2-means clustering based null distribution
of $CI$s essentially starts with the black points in the left panel,
which are more spread than the green points because of the influence\index{influence}
of the red points on the sample covariance matrix estimate. These
black points are split in half (as shown in the right panel of Figure
\ref{figK-A:ClusterToyData}) and the corresponding $WCSS$ is compared
with the overall $TSS$ (recall this Sums of Squares notation from
(\ref{eqK-A:ClusterIndex})). On the other hand the given clustering
$CI$ (the vertical green line in the right panel) is relatively closer
to 1 because the very few red outliers have a relatively small impact
on either the $WCSS$ or the $TSS$, which are thus close to each
other. Note that the resulting hypothesis test in this case is still
valid, but it is important to realize that it can be quite conservative
in the case of strongly unbalanced cluster sizes. An interesting open
problem is the development of a modified SigClust that gives better
performance in the unbalanced case. 

\begin{figure}[H]
\includegraphics[scale=0.6]{ChapterL/OODAbookChpLFigGSigClustLongOutLeft}\qquad{}\includegraphics[scale=0.5]{ChapterL/OODAbookChpLFigGSigClustLongOutRight}

\caption{Challenging (subset of the \index{Four Cluster data|)}Four Cluster
data) case for SigClust, with the long green and small red clusters
input. Given $CI$ (based on green (plus) vs. red (x) labels) is actually
larger than the null distribution, because of the severely unbalanced
class sizes. \label{figL.G:SigClustLongOutlier}}
\end{figure}

\subsection{High Dimensional SigClust\label{subsecL-B-A:HiDimSigClust}}

The above two dimensional examples are relatively simple, since good
estimates of the fit Gaussian distribution, $N_{2}\left(\boldsymbol{\mu},\boldsymbol{\varSigma}\right)$
on $\mathbb{R}^{2}$, are available. But construction of decent estimates
is much more challenging in higher dimensional cases. Some serious
reduction of this problem can be made by exploiting invariance properties
of the Cluster Index, $CI$. The first one is \emph{shift invariance}\index{shift invariance}\index{invariance!shift},
which is the fact that if an entire data set is shifted in any direction
by a constant amount the $CI$ remains the same. For this reason,
there is no need to estimate the population expected value $\boldsymbol{\mu}$
of the $N_{d}\left(\boldsymbol{\mu},\boldsymbol{\varSigma}\right)$
distribution, instead the simulations can equally well be drawn from
the $N_{d}\left(\boldsymbol{0}_{d,1},\boldsymbol{\varSigma}\right)$
(using the zero matrix notation from (\ref{eqB.G:0sMatrix})) distribution.
The second key property of the $CI$ is \emph{rotation invariance\index{rotation invariance}}\index{invariance!rotation},
where a rigid rotation of the full data set similarly leaves the $CI$
unchanged (because all the line segments in the right panel of Figure
\ref{figK-A:ClusterToyData} will keep the same lengths in a rigid
rotation). The value of rotation invariance comes from the eigen (spectral)
analysis of the underlying population covariance matrix $\boldsymbol{\varSigma}=\boldsymbol{U}\boldsymbol{\varLambda}\boldsymbol{U}^{t}$
where $\boldsymbol{U}\in O(d)$ (the \index{orthogonal group}orthogonal
group of matrix transformation\index{transformation!matrix}s defined
in Section \ref{subsec:Covariance-matrix-distances}), and $\boldsymbol{\varLambda}$
is the diagonal matrix of eigenvalues. See Section \ref{subsecP-A-B:SVD}
for more discussion of eigenanalysis. Hence by a suitable rotation
of the data, it is enough to draw the null hypothesis data from the
$N_{d}\left(\boldsymbol{0}_{d,1},\boldsymbol{\varLambda}\right)$
distribution. Given a set of eigenvalues $\lambda_{1},\cdots,\lambda_{d}$
along the diagonal of $\boldsymbol{\varLambda}$, this simulation
is straightforward by drawing standard normal vectors and rescaling
the $j$th entry by $\lambda_{j}$ for $j=1,\cdots,d$. Note that
for high $d$ this is computationally much more efficient than using
packaged multivariate Gaussian data generation functions, because
those tend to entail working through a $d\times d$ covariance matrix,
as discussed in Section \ref{subsecP-A-D:PCA-Comput}. 

It remains to estimate $\boldsymbol{\varLambda}$, which is challenging
in the case of $d>n$. There are two major problems in this case.
The first is that only the first $n$ eigenvalues (or maybe $n-1$
when mean centering\index{centering} has been done) of the sample
covariance matrix are nonzero. The second is that the nonzero sample
eigenvalues tend to over estimate their corresponding underlying theoretical
eigenvalues, because the additional energy (quantified as sums of
squares about the mean as discussed in Section \ref{secBA.1:Foundations_Terminology})
generated by the theoretical eigenvalues that are not represented
gets added in. This seems to result in a conservative null distribution,
first pointed out in \citet{Liu2008SigClust}, because simulating
data in a lower dimensional subspace results in a smaller $CI$ despite
the larger eigenvalues. An open problem is theoretical verification
that use of the sample covariance matrix $\widehat{\boldsymbol{\Sigma}}$
(defined in (\ref{eqB.F:SampleCovariance})) always leads to \index{Sig-Clust!conservative}conservative
inference in SigClust.

Some assumptions need to be made to address this. A common assumption,
which makes sense in many (but it is important to note not all) real
data situations is that there is essentially a strong low rank signal
together with some relatively small scale noise that is \emph{isotropic}\index{noise!isotropic},
i.e. has the same variance in all directions. One model for isotropic
noise is a Gaussian distribution with all eigenvalues the same (i.e.
diagonal covariance matrix). However this holds for any multivariate
distribution generated by independent marginals with the same variance
(not even the same distribution is needed). To see why, note that
for independent random variables $X_{1},X_{2}$ each with variance
$\sigma^{2}$, the variance of the projection of $\left(X_{1},X_{2}\right)$
onto the arbitrary unit direction $\left(\cos\theta,\sin\theta\right)$
is 
\begin{equation}
var\left(X_{1}\cos\theta+X_{2}\sin\theta\right)=var\left(X_{1}\right)\cos^{2}\theta+var\left(X_{2}\right)\sin^{2}\theta=\sigma^{2}.\label{eqL.M:IsoVar}
\end{equation}
Similar calculations are straightforward for directions in higher
dimensions. In the isotropic errors case, the problem of estimating
$\boldsymbol{\varLambda}$ has been effectively solved by \citet{Huang2015SigClust}. 

An example where the assumption of a strong low rank signal plus additive
isotropic noise makes sense is the Pan Cancer gene expression data
studied in Section \ref{subsecB.2.B:PanCancerExample}. Here only
the $n=50$ Kidney Cancer cases (blue diamonds in Figure \ref{figB-Q:PanCanPCA})
are studied. The data objects are vectors where each of the $d=12478$
entries represents activity of a gene. These entries are essentially
counts measured with random variation. Those counts are large enough
that on the log scale\index{scale!log|(} the random variation is
roughly Gaussian. This (appropriate for high dimension eigenvalue
estimation by SigClust) behavior is verified using the diagnostic
plot in the left panel of Figure \ref{figL.H:SigClustPCkdeQQ}. That
graphic studies the distribution of all $50\times12478=623,900$ entries
of the data matrix (called \index{pixel}pixels in titles in Figures
\ref{figL.H:SigClustPCkdeQQ} and \ref{figL.K:SigClustLungC}). Because
the full jitter plot\index{jitter plot} of this data set would result
in massive overplotting, just a random sample\index{random sample}
of 5000 green dots is shown. Those points show a distribution which
is highly concentrated near 0 with some substantial outliers. The
latter is driven by the low rank signal structure (much of which is
appropriately captured by the first principal components), while the
smaller scale noise appears in the middle. This behavior is nicely
reflected by the blue curve which is a kernel density estimate (discussed
in more detail in Section \ref{secN-A:NoHisto}). That suggests the
random noise part of this distribution is approximately Gaussian (i.
e. suggests strong potential for an isotropic background noise distribution).
The variance of this background noise distribution is very useful
for essentially interpolating those zero sample eigenvalues. 

\begin{figure}[H]
\includegraphics[scale=0.5]{ChapterL/OODAbookChpLFigHSigClustPCkde}\qquad{}\qquad{}\includegraphics[scale=0.5]{ChapterL/OODAbookChpLFigHSigClustPCQQ}

\caption{SigClust diagnostics for the Kidney Cancer type from the Pan Cancer
data\index{Pan Cancer data|)}. Left panel studies the distribution
of all entries of the data matrix, suggesting an isotropic Gaussian
background noise distribution. Right panel is fit Gaussian Q-Q plot
indicating good fit of the suggested distribution within the quartiles.
\label{figL.H:SigClustPCkdeQQ}}
\end{figure}
The full sample standard deviation (about $1.15$ for this data set)
provides an inappropriate estimate of the background standard deviation,
because it strongly feels the relatively few, but quite large, signal
values. Scale measures that are much less affected by large values
have been developed in the field of \emph{robust statistics}, which
is overviewed in Chapter \ref{chapO:Robust-Methods}. For the estimation
of the background noise standard deviation, it is desirable to use
a method which nearly ignores values outside the quartiles such as
the (Gaussian rescaled) \index{Median Absolute Deviation}Median Absolute
Deviation (MAD) from the median, $\widehat{\sigma}_{MAD}$. For a
data set of size $m$ (in Figure \ref{figL.H:SigClustPCkdeQQ} $m=d\times n$),
$\widetilde{x}_{1},\cdots,\widetilde{x}_{m}$, this is defined as
\begin{equation}
\widehat{\sigma}_{MAD}=\underset{i=1,\cdots,m}{\mathrm{median}}\left(\left|\widetilde{x}_{i}-\underset{j=1,\cdots,m}{\mathrm{median}}\left(\widetilde{x}_{j}\right)\right|\right)/\varPhi^{-1}\left(0.75\right),\label{eqL.A:MAD}
\end{equation}
where $\varPhi$ denotes the standard Gaussian cumulative distribution
function. The purpose of that rescaling by $\varPhi^{-1}\left(0.75\right)$
(which is the theoretical $MAD$ of the standard Gaussian distribution)
is to put the $\widehat{\sigma}_{MAD}$ on a Gaussian scale\index{scale!Gaussian}
in the sense that when applied to a large Gaussian data set it appropriately
estimates the population standard deviation. The Gaussian density
with standard deviation taken to be $\widehat{\sigma}_{MAD}=0.781$
(and centered at the sample median $\widehat{\mu}_{MED}=-0.02$ which
similarly ignores the outlying values) is shown using the red dashed
curve in the left panel of Figure \ref{figL.H:SigClustPCkdeQQ}. Note
that it appears to fit the central ``background noise'' part of
the distribution quite well. A more precise view of this goodness
of fit appears in the QQ plots\index{QQ plots} shown in the right
panel of Figure \ref{figL.H:SigClustPCkdeQQ}. This version of the
QQ plot shows the sorted data $\widetilde{x}_{(1)},\cdots,\widetilde{x}_{(m)}$
on the vertical axis versus the corresponding estimated Gaussian quantile\index{quantile}s
$\widehat{\mu}_{MED}+\widehat{\sigma}_{MAD}\varPhi^{-1}\left(\frac{1}{m+1}\right),\cdots,\widehat{\mu}_{MED}+\widehat{\sigma}_{MAD}\varPhi^{-1}\left(\frac{m}{m+1}\right)$
on the horizontal axis, as the red curve. When the given distribution
fits the data, the QQ curve should roughly follow the green curve,
which is the 45 degree line. In this case the fit is good for the
part of the data between the quartiles, but quite bad for the outlying
data points (which again represent important signal in the data) suggesting
that the background noise standard deviation has been effectively
estimated. 

Use of this isotropic background noise estimate to provide a reasonable
full set of estimated eigenvalues is shown in Figure \ref{figL.I:SigClustPCplotEV}.
Each panel shows various eigenvalues on different scales. The top
two panels study the full set of $d=12478$ eigenvalues, while the
bottom panels are the same as the top, but are zoomed into just the
first 100 eigenvalues (as the structure of those important eigenvalues
is essentially squished into the vertical axis because $100\ll12478$).
Both panels on the left use the original variance scale on the vertical
axis, while the $\log_{10}$ scale is used in both right hand panels.
In all panels the sample eigenvalues (this is the scree plot\index{scree plot}
introduced in Figure \ref{figB-Y:Screeplots}) are shown using black
circles. Note that there a few extremely large eigenvalues and many
much smaller ones. Only $49=n-1$ of the log scale\index{scale!log|)}
sample eigenvalues appear (because the $\log_{10}0$ is undefined
and thus not plotted). 

The magenta horizontal line shows the background variance that is
estimating the variance of the noise part of the model. An important
goal of the eigenvalue estimation is to keep the total energy of the
estimated eigenvalues (i.e. their sum) equal to the total energy of
the data (i.e. the sum of squares about the mean, which is the sum
of the eigenvalues). Replacing the 0 sample eigenvalues with the estimated
background variance adds a lot to the total, which thus motivates
some shrinkage of the nonzero sample eigenvalues. Appropriate shrinkage
formulas have been developed in \citet{Huang2015SigClust}, and the
resulting estimated eigenvalues are shown as the red dashed curve
in each panel of Figure \ref{figL.I:SigClustPCplotEV}. Note that
on the right side (indeed most) of each panel, this dashed red curve
is the same as the magenta line. It rises for the first few eigenvalues,
but represents substantial shrinkage from the black circle sample
eigenvalues, again with the goal of making the total energy the same.
For some data sets, the total energy in the sample eigenvalues (their
sum) may be less than $d$ times the background variance. In this
case, the background noise is clearly not isotropic so this estimation
scheme cannot be used. Hence one must resort to simply using the sample
eigenvalues. As noted in \citet{Liu2008SigClust}, this results in
a conservative SigClust, so the result is valid but the test may be
less powerful in those situations.

\begin{figure}[H]
\includegraphics[scale=0.65]{ChapterL/OODAbookChpLFigISigClustPCplotEV}

\caption{Various estimated eigenvalues for the Kidney Cancer data. Sample eigenvalues
are shown as black circles. Improved estimates, that both incorporating
the estimated isotropic background noise variance (shown as the horizontal
magenta line) and also keep the total energy constant, are shown as
the red dashed curve. \label{figL.I:SigClustPCplotEV}}
\end{figure}
 

The results of using SigClust to investigate splitting the Kidney
Cancer data from Section \ref{subsecB.2.B:PanCancerExample} (thus
working in the second mode described above) are summarized in the
p-value plots of Figure \ref{figL.J:SigClustPCpvals}. The left panel
uses Gaussian data simulated from the carefully estimated eigenvalues
shown as the red dashed curve in Figure \ref{figL.I:SigClustPCplotEV}.
Note that none of the Gaussian simulated black dots lies to the left
of the vertical green line, giving an empirical p-value\index{empirical p-value|)}
of 0. This is strong evidence that there are at least two clusters
in the Kidney Cancer data indicating the existence of cancer subtypes,
which is consistent with the results of \citet{Ricketts2018TCGAKIRC}
and with the analysis in Section \ref{subsecN-D.E:KidneySiZer}. 

The right panel of Figure \ref{figL.J:SigClustPCpvals} shows the
result of replacing the careful eigenvalue estimates (the red dashed
curve in Figure \ref{figL.I:SigClustPCplotEV}) with the sample eigenvalues
(the black circles). With that estimate the p-value is only about
0.14, no longer statistically significant in the classical 0.05 sense.
So this is an example where the method of eigenvalue estimation clearly
matters and the sophisticated approach of \citet{Huang2015SigClust}
is worthwhile. It also is an example of the conservative nature (especially
in high dimensions) of the use of naive sample eigenvalues in SigClust
discussed above.

\begin{figure}[H]
\includegraphics[scale=0.5]{ChapterL/OODAbookChpLFigJSigClustPCpvals}\qquad{}\includegraphics[scale=0.5]{ChapterL/OODAbookChpLFigJSigClustPCpvalsSC} 

\caption{SigClust p-values for the Kidney Cancer data. Left panel shows that
careful eigenvalue estimation (red dashed curve in Figure \ref{figL.I:SigClustPCplotEV})
gives a statistically significant result, indicating this cancer type
has potentially important subtypes. Right panel reveals that naive
sample eigenvalue estimation (black circles) gives a too conservative
(not even significant) result in this same case. \label{figL.J:SigClustPCpvals}}
\end{figure}
For effective use of SigClust it is also important to realize there
are cases where the sophisticated shrinkage eigenvalue estimates demonstrated
in Figure \ref{figL.I:SigClustPCplotEV} are inappropriate. An example
of this is the Lung Cancer data studied in Section \ref{subsecB.2.A:LungCancerData}
where clusters were seen to be of keen interest. Figure \ref{figL.K:SigClustLungC}
investigates whether the red and blue clusters visually brush\index{brushing}ed
in Figure \ref{fig:LungCancerPCABrushed} are a statistically significant
separation using the confirmatory version of SigClust (i.e. based
on given, in this case manually selected, cluster labels). The left
panel is the first SigClust diagnostic plot showing the distribution
of all entries in the data matrix, which is far different from that
shown in Figure \ref{figL.H:SigClustPCkdeQQ}. While both distributions
are based on logs of counts (plus 1), here there are far more 0s and
many other very small discrete values (appearing at $\log_{10}$ of
small positive integers), so no Gaussian background distribution is
apparent. In this case it is recommended to use the conservative sample
covariance eigenvalues in the application of confirmatory SigClust.
One can plunge ahead with the formulas developed above for the isotropic
background variance case, but there are serious consequences for this
data set. In particular $d$ times the estimated background noise
is now larger than the total sum of squares about the mean (energy
in the data), so the estimated eigenvalues are all simply the background
variance, i.e. the dashed red curve in the analog of Figure \ref{figL.I:SigClustPCplotEV}
is completely flat. This puts the simulated Gaussian data at the heart
of the SigClust p-value computation into the HDLSS isotropic domain
that will be studied in Section \ref{secM-B:HDLSS}, where all simulated
data are essentially equidistant from each other, as characterized
in (\ref{eqM.C:PairNormDist}). This results in simulated $CI$s being
very similar to each other, so the Z-score (whose denominator is the
standard deviation of the black dots) is unrealistically negative
(in the thousands in this case, since the standard deviation of the
black dots is so small). Hence this violation of the isotropic background
noise distribution can result in a seriously anti-conservative SigClust
result. Thus the diagnostics are important in high dimensional situations. 

The result using the sample eigenvalues, recommended when the isotropic\index{noise!isotropic}
background diagnostics fail, is shown in the right panel of Figure
\ref{figL.K:SigClustLungC}. This shows that the red and blue clusters
in Figure \ref{fig:LungCancerPCABrushed} are clearly statistically
significant. This is consistent with the explanation of those clusters
from the discussion of Figure \ref{fig:LungCancerLogCurvesBrushed}.

\begin{figure}[H]
\includegraphics[scale=0.45]{ChapterL/OODAbookChpLFigKSigClustLungCkde}\qquad{}\includegraphics[scale=0.45]{ChapterL/OODAbookChpLFigKSigClustLungCpval}

\caption{SigClust for the red vs. blue clusters in the Lung Cancer data of
Figures \ref{fig:LungCancerLogCurves}-\ref{fig:LungCancerLogCurvesBrushed}.
Diagnostic plot in the left panel reveals no underlying Gaussian background
noise. SigClust using sample covariance eigenvalues shows these clusters
are statistically significant in the right panel. \label{figL.K:SigClustLungC}}
\end{figure}
Figure \ref{figL.L:ScatPlotsLungC} shows another analysis of the
Lung Cancer data studied in Figures \ref{fig:LungCancerLogCurves}-\ref{fig:LungCancerLogCurvesBrushed}.
An application of SigClust, to separate the blue cluster from the
union of the red and gold (shown as an intermediate orange color in
Figure \ref{figL.L:ScatPlotsLungC}) gave a p-value of 1. The p-value
plot looks much like the right panel of Figure \ref{figL.G:SigClustLongOutlier},
and hence is not shown here. The scatterplots of Figure \ref{figL.L:ScatPlotsLungC}
show that this was caused by the orange cluster really consisting
of two clusters so the WCSS is much bigger than the WCSS that results
from just splitting all of the data. This shows that SigClust\index{SigClust|)}
is best used together with data visualization\index{visualization|)},
not as a blind hypothesis test. This data set also provides an example
where the DWD\index{Distance Weighted Discrimination|)} version of
DiProPerm\index{DiProPerm} (which had a Z-score of 55.3) gave a much
stronger result than the MD version (Z-score of 9.2). The left panel
of Figure \ref{figL.L:ScatPlotsLungC} shows the projections of the
data onto the orthogonal MD direction. Note that \index{Mean Difference}MD
gives a rather poor separation of the two classes. On the other hand,
the projection onto the DWD projection in the right panel (again together
with the orthogonal PC1 direction) shows much better class separation
which is why the DWD Z-score\index{Z-score|)} is so much better.
Of course the use of DiProPerm is not really appropriate here as these
clusters were visually chosen. 

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterL/OODAbookChpLFigLScatPlotLungC}\caption{Scatterplots of projections of Lung Cancer data\index{Lung Cancer data|)}
from Figures \ref{fig:LungCancerLogCurves}-\ref{fig:LungCancerLogCurvesBrushed}.
\label{figL.L:ScatPlotsLungC}}

\end{figure}
Other approaches to testing significance of clustering can be found
in \citet{McShane2002ClusterSig}, \citet{Tibshirani2001GapStat,Tibshirani2005ClusterVal}
and \citet{Suzuki2006Pvclust}. A different approach is Bayesian inference
based on normal mixtures, as in \citet{RichardsonGreen1997MixCompNum}
and \citet{BouveyronRaftery2019MClust}. Yet another interesting way
to think of importance of clusters is via their stability using \emph{\index{consensus clustering}consensus
clustering\index{clustering|)}} proposed by \citet{Monti2003ConsensusClustering},
with a Bayesian variation by \citet{Lock2013ConsCluster}.\index{confirmatory analysis|)}
