
\chapter{Robust Methods\label{chapO:Robust-Methods}}

As\index{robust methods|(} discussed at many points in previous chapters,
\emph{robustness} issues arise frequently in OODA. Good overview of
a large body of research done in this area can be found in the monographs
\citet{HampelEtAl2011Robust}, \citet{HuberRonchetti2009RobustStat}
and \citet{StaudteSheather1990Robust}, with a more recent survey
in \citet{Clarke2018RobustnessBook}. A general view of the area is
that it studies, and recommends remedies for, violation of various
classical statistical modeling assumptions. 

The great bulk of this work has been motivated by the fact that under
the very typical assumption of a Gaussian error distribution \emph{outliers\index{outliers}}
occur with vanishingly small probability. Yet they can be rather common
in real data, and have strong potential to seriously disrupt classical
statistical methods which tend to ignore them, as illustrated in Figure
\ref{figO.A:Outliers}. The left panel shows a toy example of $n=10$
simulated realizations from the $N_{2}\left(\left[\begin{array}{c}
2\\
-2
\end{array}\right],\left[\begin{array}{cc}
0.49 & 0.49\\
0.49 & 0.03
\end{array}\right]\right)$ distribution shown as + signs, with one point deliberately moved
far away. Note that the sample mean, shown as the thick + sign with
a circle around it, is a poor approximation of the center of the remaining
point cloud, because the outlier has a strong \index{influence|(}influence
on it (in particular with weight $\frac{1}{n}=0.1$). This pulls it
far outside the elliptical contour of the Gaussian density that contains
say 95\% of the probability mass, hence the sample mean is a very
poor estimate of $\mu=\left[\begin{array}{cc}
2 & -2\end{array}\right]^{t}$. 

A parallel impact of an outlier on PCA\index{Principal Component Analysis!and robustness|(}
is shown in the right panel of Figure \ref{figO.A:Outliers}. This
time $n=40$ points were drawn from the $N_{2}\left(\left[\begin{array}{c}
-2\\
-2
\end{array}\right],\left[\begin{array}{cc}
1 & 1\\
1 & 0.15^{2}
\end{array}\right]\right)$ distribution (also + signs), one of which was similarly moved away.
While the sample mean (at the intersection of the heavy dashed line
segments) is somewhat affected, it at least lies within the point
cloud of the bulk of the data. However the first PC direction (indicated
by the longer dashed line segment) clearly does not reflect the dominant
mode of variation\index{modes of variation|(} in the overwhelming
majority of the data set. That is because PC1 finds the direction
of maximal projected variance, and sample variance (which is driven
by squared distance\index{distance|(}s, recall how that had a strong
impact on the toy data set on $S^{2}$ in Figure \ref{figG.C:BadTangentPlane})
is notoriously impacted by outlying data points. Here this happened
to the extent of actually being larger than the sample variance along
the major axis of the remaining data cloud, thus pulling off the PC1
direction as seen. Scores plots (first seen in Figure \ref{fig:Mortality-PC1}
and used in many different ways above) provide a very effective diagnostic
for detecting undue influence from outliers, as seen in the right
hand column of Figure \ref{figO.C:10dToynOut}. The second PC in Figure
\ref{figO.A:Outliers} is shown as the shorter dashed line segment
to give an impression of these as axes of the PC coordinate system.
In robust statistics many approaches to dealing with such challenges
have been considered from a number of interesting perspectives. 

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterO/OODAbookChpOFigAOutliersEg}

\caption{Toy examples illustrating impact of outliers on the sample mean and
PCA. The single outlier in the left panel pulls the mean (circled
+ sign) actually outside of the convex hull of the remaining Gaussian
data points. On the right the single outlier pulls the PC1 direction
(longer dashed line) almost orthogonal to the major axis of the rest
of the again Gaussian point cloud. \label{figO.A:Outliers}}
\end{figure}
While the study of outliers has dominated research in the area, robustness
against other violations has been studied as well, such as the violation
of the typical \emph{independence} assumption, as studied in \citet{Beran1994LongMemory},
\citet{Beran2017Longmemory} and \citet{ClarkeHall2009RobustIndep}.
Another important violation of typical (i.e. Gaussian error model)
assumptions is treated in the area of robustness against \emph{data
heterogeneity}\index{data heterogeneity}, as discussed in Section
\ref{secJ-D:DWD}. The fundamental idea there is that modern large
data sets are often created by combining smaller data sets, that are
frequently not identically distributed. In such cases a Gaussian mixture
error model can be far more appropriate than the common single Gaussian.
See \citet{Marron2017BigData} for more discussion, and \citet{Buhlmann2015Inhomo,MeinshausenBuhlmann2015InhomoData}
for leadership research in that area.

An overview of some controversial robustness issues is in Section
\ref{secO.A:RobustControv}. Specific discussion of robust methods
that have proven to be useful in OODA, e.g. techniques well suited
for handling high dimensional data, appears in Section \ref{secO.B:RobustHiD}.
Two specific case studies, involving dealing with surprisingly different
types of outliers in high dimensions are presented in Sections \ref{secO.B.A:CorneaData}
and \ref{secO.B.B:GWAS}. A few other robustness topics are discussed
in Section \ref{secO.C:OtherRobust}. 

\section{Robustness Controversies\label{secO.A:RobustControv}}

Perhaps because of the many potential ways to approach it, the field
of robust statistics has seen its share of controversy and very robust
debate over the years. One topic is \emph{\index{outlier deletion}outlier
deletion} versus \emph{robust methodology}. From the former viewpoint,
it is the responsibility of the data analyst to find outliers, and
after careful consideration to delete them from the data set when
that is sensible. A number of approaches to this can be found in \citet{Barnett1994Outliers}.
The latter approach has been to shift away from explicit identification
of outliers towards the development of methods that are simply far
less sensitive to them. The large majority of researchers in robust
statistics have chosen this latter approach for a number of reasons,
including: 
\begin{itemize}
\item Outliers can be quite challenging to even define, making them difficult
to delete in an objective way.
\item Outlier can also be difficult to discover, especially in high dimensional
and complex data situations such as arise in OODA. 
\item Outlier deletion is a time consuming process in a world where all
too many analysts tend to be increasingly pressed for time, as is
apparent from the fact that so little data visualization\index{visualization}
is frequently done in many modern data analyses.
\item Outlier deletion may involve just too large a part of the data set,
including some data objects that still contain useful information,
as illustrated using the Cornea data in Section \ref{secO.B.A:CorneaData}. 
\end{itemize}
There has also been strong controversy even among those who advocate
development of robust statistical methods over outlier deletion. The
key to understanding the two sides in the debate is to consider the
potential genesis of the outliers. One natural source is recording
errors, which can easily happen while writing down numbers (e.g. missing
a decimal point), or else through unnoticed failure of a measuring
device. In that case, outliers are simply bad data (i.e. contain no
useful information), so the statistical methods should aim to completely
ignore them (i.e. give them weight 0) as done by outlier deletion.
But outliers also naturally appear as just atypically large values,
that are not a mistake but are an unusual yet genuine part of the
data being studied, e.g. as illustrated using the Drug Discovery data\index{Drug Discovery data}
in Figure \ref{figD-A:DrugPCARaw} and the Pan Can data in Figure
\ref{figD-L:TCGA-Log2vsRaw}. In such cases the outliers contain useful
information, but robust methods should keep them from unduly dominating
the analysis (as happened for both Gaussian toy examples in Figure
\ref{figO.A:Outliers}). The canonical example of such a method is
the sample \index{median|(}median in $\mathbb{R}^{1}$, which feels
a small but appropriate influence from every data point, even those
that are very far from the rest of the data. The former position of
implicit outlier removal was strongly advocated by the Hampel school
(e.g in early editions of \citet{HampelEtAl2011Robust}), while the
latter of accommodating yet downweighting outliers was equally strongly
pushed by Huber (e.g. in early editions of \citet{HuberRonchetti2009RobustStat}).
Of course from a broader perspective, each approach can be appropriate
depending on the context. The challenge for the data analyst is to
determine which to use in any particular case. 

\section{Robust Methods for OODA\label{secO.B:RobustHiD}}

As noted above the challenges of OODA have motivated the re-thinking
of some robust statistical ideas. An important example is the notion
of \emph{mean outliers\index{outliers!mean}} versus \emph{shape outliers\index{outliers!shape}}
in FDA\index{functional data analysis}, as introduced by \citet{DaiGenton2018FDAOutliers,DaiGenton2019DirnOutlyers,DaiGenton2020OutlierTaxonomy}
(who give several interesting examples clarifying how challenging
it can be to even define outliers). The latter concept is illustrated
in Figure \ref{figO.C:10dToynOut}. The toy data set there starts
with the $n=50$ 10-d Tilted Parabolas\index{Tilted Parabolas data|(}
in Figures \ref{fig:10d-ToyPCA} and \ref{fig:10d-ToyPCA-MinWhiteSpace}
shown in black. The additional red curve that has been added demonstrates
the concept of shape outlier. Note that the values of the red curve
lie entirely in the vertical range $[0,30]$ while the black curves
occupy a larger range of $[-10,40]$. However, instead of having that
generally very smooth parabolic shape it has a much higher frequency
(it is actually a coarsely discretized sinusoid) which puts it into
a far different region of the space of curves. Perhaps as expected
from the right panel of Figure \ref{figO.A:Outliers}, this red outlier
has a strong impact on PCA (especially compared to Figure \ref{fig:10d-ToyPCA-MinWhiteSpace}).
A careful look at the sample mean in the top center panel of Figure
\ref{figO.C:10dToynOut} shows that the outlier has some influence
there, in terms of some quite angular corners corresponding to the
valleys in the red curve compared to the mean of the Tilted Parabolas
data in Figure \ref{fig:10d-ToyPCA}. The first mode of variation\index{modes of variation}
is generally similar to the original analysis in terms of mostly being
an overall up and down mode, although the red outlier clearly has
a strong influence on the shape of the loadings vector. The second
mode of variation is completely dominated by the outlier. There is
hardly any variation among the black curves in the loadings plot in
the left panel of the third row. The outlier also stands out starkly
in the scores plot on the right of the third row, highlighting how
far the outlier is from the Gaussian bulk of the data in the 10-d
feature space\index{feature space}. As noted above, in general such
scores plots are a good diagnostic for indicating the presence of
the outlier. The third mode of variation (fourth row) is essentially
the random tilt mode which appears as the second mode in Figure \ref{fig:10d-ToyPCA-MinWhiteSpace},
although again there has been some corruption by the outlier. As in
the bottom row of the earlier analysis, the residuals from the first
three modes (bottom center panel) are relatively much smaller (recalling
the lesson of vertical axis choice conveyed by Figures \ref{fig:10d-ToyPCA}
and \ref{fig:10d-ToyPCA-MinWhiteSpace}).

\begin{figure}[H]
\includegraphics[scale=0.7]{ChapterO/OODAbookChpOFigC10dToynOut1}

\includegraphics[scale=0.7]{ChapterO/OODAbookChpOFigC10dToynOut2}

\caption{Parabolas and Outlier toy example in the same format as the PCA of
Figure \ref{fig:10d-ToyPCA-MinWhiteSpace}. Tilted Parabola data from
there are shown in black, with a red shape outlier added, illustrating
how outliers can have a major impact on PCA, and how scores plots
can provide useful outlier diagnostics. \label{figO.C:10dToynOut}}
\end{figure}
Among the many approaches to mitigation of outliers that have been
proposed in the literature, two approaches that are intuitively appealing
(and are seen in Section \ref{secO.B.A:CorneaData} to have good high
dimensional properties) are illustrated in Figure \ref{figO.B:OutliersnEsts}.
The toy data sets in the two panels are the same as in Figure \ref{figO.A:Outliers}.
The left panel shows a robust \index{median!robust}median that, as
noted in Section \ref{secE.1:Fr=0000E9chet-mean}, appears to have
been re-discovered and re-named several times. \citet{Haldane1948MultivarMedian}
named it the geometric median\index{median!geometric}\index{geometric median}.
It was shown to be unique in the case $d>1$ (for $d=1$ this is the
ordinary median which of course is not unique for $n$ even) by \citet{Milasevic1987UniqueMedian}.
That paper used the name spatial median\index{median!spatial}\index{spatial median},
which is common in the multivariate rank statistics\index{multivariate rank statistics}
literature, see \citet{Oja1995MultSignRank} and \citet{Oja2010MultivarRanks}.
In the robustness literature it is usually called Huber's \index{$L^{1}$ M-estimate}$L^{1}$
M-estimate, using terminology from \citet{HuberRonchetti2009RobustStat}.
From that perspective it is the $p=1$ special case of the $L^{p}$
\index{M-estimate|(}M-estimate defined for random vectors $\widetilde{\boldsymbol{x}}_{1},\cdots,\widetilde{\boldsymbol{x}}_{n}$
centered at $\boldsymbol{\theta}\in\mathbb{R}^{d}$ as 
\begin{equation}
\arg\underset{\boldsymbol{\theta}\in\mathbb{R}^{d}}{\min}\sum_{i=1}^{n}\left\Vert \widetilde{\boldsymbol{x}}_{i}-\boldsymbol{\theta}\right\Vert _{2}^{p}.\label{eqO.A:HubersLpMest}
\end{equation}
Note that for $p=2$ this is the Fréchet mean\index{Fréchet mean}
defined at (\ref{eq:ChpEEqBFrechetMean}) with respect to Euclidean
distance\index{distance|)}. Another important special case is $p=1$,
the Fréchet median\index{median!Fréchet}\index{Fréchet median} from
(\ref{eqE.H:FrechetMedian}).

A simple iterative algorithm due to \citet{Gower1966DS}, provides
insight into the workings of this estimator of centrality and why
it is robust. In the case $p=1$, it is based on setting the vector
of $d$ partial derivatives of the sum in (\ref{eqO.A:HubersLpMest})
to $\boldsymbol{0}$. Using the notation $\boldsymbol{\theta}=\left[\begin{array}{ccc}
\theta_{1} & \cdots & \theta_{d}\end{array}\right]^{t}$ and $\widetilde{\boldsymbol{x}}_{i}=\left[\begin{array}{ccc}
\widetilde{x}_{i1} & \cdots & \widetilde{x}_{id}\end{array}\right]^{t}$ this gives, for $j'=1,\cdots,d$ 
\[
0=\frac{\partial}{\partial\theta_{j'}}\sum_{i=1}^{n}\left(\sum_{j=1}^{d}\left(\widetilde{x}_{ij}-\theta_{j}\right)^{2}\right)^{1/2}=\sum_{i=1}^{n}\frac{1}{2}\left(\sum_{j=1}^{d}\left(\widetilde{x}_{ij}-\theta_{j}\right)^{2}\right)^{-1/2}2\left(\widetilde{x}_{ij'}-\theta_{j'}\right).
\]
Multiplying by $\frac{1}{n}$ and putting these back into vectors
gives 
\begin{equation}
\mathbf{0}=\frac{1}{n}\sum_{i=1}^{n}\frac{\widetilde{\boldsymbol{x}}_{i}-\boldsymbol{\boldsymbol{\theta}}}{\left\Vert \widetilde{\boldsymbol{x}}_{i}-\boldsymbol{\boldsymbol{\theta}}\right\Vert _{2}}.\label{eqO.B:SetDerivZero}
\end{equation}
Insight into this quantity comes from the example in the left panel
of Figure \ref{figO.B:OutliersnEsts}. The + signs are the points
$\widetilde{\boldsymbol{x}}_{1},\cdots,\widetilde{\boldsymbol{x}}_{n}$,
and a candidate value of $\boldsymbol{\theta}$ is the bold x sign
at $(1,2)$. Note that each $\widetilde{\boldsymbol{x}}_{i}-\boldsymbol{\theta}$
is the vector pointing from $\boldsymbol{\theta}$ to $\widetilde{\boldsymbol{x}}_{i}$,
and dividing by its norm (i.e. length) projects it onto the circle
of radius 1 centered at $\boldsymbol{\theta}$. For the case of $\boldsymbol{\theta}$
at that x sign, these $\boldsymbol{\theta}$ centered projections
are shown as small circles lying on the large dashed circle centered
at $\boldsymbol{\theta}$. The equation (\ref{eqO.B:SetDerivZero})
is solved by the choice of $\boldsymbol{\theta}$ that makes the average
of these projections, shown as a bold circle, equal to $\boldsymbol{\theta}$.
That is not achieved by $\boldsymbol{\theta}$ at that x sign. Gower's
algorithm is iterative movement of $\boldsymbol{\theta}$ to the bold
circle at each step, until convergence. The final result is shown
using the the bold circled x sign, and its corresponding large dashed
circle with projections is also shown. This gives a much more sensible
notion of center of this data set than the sample mean (shown as in
Figure \ref{figO.A:Outliers} using the bold circled + sign). 

The key to the robustness of this notion of center is that the outlier
becomes downweighted by projection onto the circle. Note that this
estimator is independent of the radius of the circle (equation (\ref{eqO.B:SetDerivZero})
can be arbitrarily rescaled), and in fact the circle shown in the
figure has a radius larger than $1$ to give a better visual impression
of the projected small circles. As noted in Section \ref{secE.1:Fr=0000E9chet-mean},
a version of this robust notion of center - the Fréchet median - that
is especially adapted to data objects on a manifold\index{manifold}
as discussed in Chapter \ref{chapG:Manifold-Data-Analysis}, can be
found in \citet{Fletcher2009GeometricMedian}.

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterO/OODAbookChpOFigBOutliersnEsts}

\caption{Illustration of Huber's $L^{1}$ M-estimate of the mean (left panel)
and Spherical PCA (right panel). Toy Gaussian data sets are the same
as Figure \ref{eqO.A:HubersLpMest}. \label{figO.B:OutliersnEsts} }
\end{figure}
The right panel of Figure \ref{figO.B:OutliersnEsts} shows how this
basic idea can be simply extended to give a robust version of PCA.
Recall that straightforward PCA applied to this data set gave a first
PC direction pointing towards the first quadrant in Figure \ref{figO.A:Outliers}
while the bulk of the data objects (the + signs) clearly suggest a
nearly orthogonal trend. That direction can be found by using the
same type of projection used to understand Huber's $L^{1}$ M-estimate
in the left panel, in particular onto a sphere centered there. Again
the small circles on the dashed circle are projections of all of the
data points. Note that the main sausage of data is projected to ice
caps on this one dimensional sphere. The outlier now has much reduced
influence, and is comparable to the projections of a few points near
the center. PCA on these projections finds a first direction driven
by the ice caps, thus following the main body of the data, as shown
by the long solid line. Again the PC2 direction is the shorter orthogonal
line indicating the PC coordinate system. This method was called \emph{spherical
PCA\index{spherical PCA|(}} by\index{Principal Component Analysis!spherical}
\citet{Locantore1999SphericalPCA}, and was the basis of a very effective
robust analysis of cornea image data objects as discussed in Section
\ref{secO.B.A:CorneaData}. Its perhaps surprisingly poor performance
in very high dimensions is explored using Genome Wide Association
data in Section \ref{secO.B.B:GWAS}. This same approach to robust
PCA was independently discovered and called \emph{\index{PCA!on the signs}PCA
on the signs} by \citet{Oja1995MultSignRank}.

Figure \ref{figO.D:10dToynOutSPCA} shows the results of a spherical
PCA of the Parabolas and Outlier toy data set in Figure \ref{figO.C:10dToynOut}.
The first mode of variation continues to be mostly curves moving up
and down together, but now the loadings vector is less influenced
by the outlier, hence closer to the first mode of variation of the
original Tilted Parabolas\index{Tilted Parabolas data} data in Figure
\ref{fig:10d-ToyPCA-MinWhiteSpace} (i.e. to that of the bulk of the
data). The strong robustness properties of spherical PCA are apparent
in the second mode of variation which is now driven by the random
tilt, again an important property of the bulk of the data as seen
in PC2 of Figure \ref{fig:10d-ToyPCA-MinWhiteSpace}. As the outlier
still represents interesting variation, it is appropriate that it
now appears in the third mode of variation, as it is clearly more
important than the small variation which remains. Some of the outlier
remains in the PC3 Residuals (bottom center panel) which is an artifact
of the nonlinear analysis done by spherical PCA. The robustness of
Huber's $L^{1}$ M-estimate\index{M-estimate|)} in the top center
panel appears as less angularity at the corners than in the sample
mean in Figure \ref{figO.C:10dToynOut}. 

\begin{figure}[H]
\includegraphics[scale=0.7]{ChapterO/OODAbookChpOFigD10dOutSPCA1}

\includegraphics[scale=0.7]{ChapterO/OODAbookChpOFigD10dOutSPCA2}

\caption{Spherical PCA of the Parabolas and Outlier data in (nearly) the format
of Figure \ref{figO.C:10dToynOut}. First and second modes of variation
are now much closer to those of the bulk of the data as shown in Figure
\ref{fig:10d-ToyPCA-MinWhiteSpace}. Lessened influence of the outlier
is reflected by it only appearing in the third mode. \label{figO.D:10dToynOutSPCA}}
\end{figure}
One difference between Figure \ref{figO.D:10dToynOutSPCA} and Figures
\ref{figO.C:10dToynOut} and \ref{fig:10d-ToyPCA-MinWhiteSpace} is
the scree plot\index{scree plot} in the upper right panel. Recall
from Figure \ref{figB-Y:Screeplots} that those showed a red curve
and circles indicating the fraction of the energy (measured using
sums of squares) explained by each mode of variation, i.e. PC component.
The blue dashed curve and + signs show the cumulative version of these
fractions. For spherical PCA as studied in Figure \ref{figO.D:10dToynOutSPCA}
there are two such energy decompositions available and of interest.
The first is the conventional PCA of the data projected to the sphere
(e.g. the small circles in the right panel of Figure \ref{figO.B:OutliersnEsts}).
This is shown as the red curve decomposition, which is monotone decreasing
because it reflects the eigen analysis\index{eigen analysis|(} underlying
this PCA (which determines the directions being used in the modes
of variation). But perhaps of even more interest may be the energy
in these spherical modes of variation in terms of the original data
objects, e.g. the + signs in the right panel of Figure \ref{figO.B:OutliersnEsts}.
These are also the sums of squares of the curves in the left columns
in Figure \ref{figO.B:OutliersnEsts}. This decomposition is shown
using red circles for the fractions, and blue + signs for the corresponding
cumulatives. These are quite different, and not even monotone. That
is because spherical PCA puts the outlier into the third mode of variation,
even though its projected variance is larger than the variance of
the random tilt mode of variation. Recall this was the reason that
the robust spherical PCA is useful for this data set.

Spherical PCA was analyzed using the HDLSS\index{High Dimension Low Sample Size|(}
asymptotic methods of Section \ref{secM-B:HDLSS} by \citet{ZhouMarron2015HDLSS}.
That paper first showed that under classical Gaussian assumptions,
spherical PCA has the same good \index{consistency}consistency /
strong inconsistency\index{inconsistency!strong} properties as conventional
PCA that were quantified in (\ref{eqM.F:ConsistStrongInconsist}).
Next they showed that for a sufficiently strong sequence of outliers,
PCA could be asymptotically pulled off in the fashion shown in the
right panel of Figure \ref{figO.A:Outliers}, while spherical PCA
gave the correct robust solution in that limit.

\subsection{Case Study: Cornea Curvature Data\label{secO.B.A:CorneaData}}

Spherical\index{data objects!images|(} PCA was actually developed
for the analysis of the Cornea Curvature data\index{Cornea Curvature data|(}
set described in \citet{Locantore1999SphericalPCA}. The cornea is
the outer surface of the eye, and its curvature is critical to vision
because most of the refraction of light entering the eye occurs there.
The data objects in that study were images as shown in Figure \ref{figO.E:RawCorneaData}.
Color codes radial curvature (i.e. along rays emanating from the center)
which is the curvature component with the most impact on visual acuity.
Stronger radial curvature is represented by warmer colors. Representation
of these functions on the disk was done using the \emph{Zernike orthogonal
basis}\index{orthogonal basis!Zernike}. That is a special case of
orthogonal basis data object representation as discussed in Section
\ref{secBA.2:OverviewObjFeat}. The Zernike basis for functions on
the disc is defined in polar coordinates as a tensor product of the
Fourier\index{Fourier!basis} basis in the angular direction, and
a system of \index{Jacobi polynomials}Jacobi\index{orthogonal basis!Jacobi}
polynomials (carefully chosen to avoid a singularity at the origin)
in the radial direction. See \citet{Schwiegerling1995Zernike} and
\citet{Born2013Optics} for further details. The $9$ images shown
in Figure \ref{figO.E:RawCorneaData}, which is Figure 3 of \citet{Locantore1999SphericalPCA},
were selected to give an impression of the full set of $n=43$. These
are reconstructions of least squares Zernike fits, based on the data
object choice of feature vectors having size $d=66$ (which gives
good noise reduction while maintaining image aspects of clinical interest).
The three cases on the bottom row show extreme edge effects. These
are artifacts of the data acquisition, which is done by collecting
reflected light from the cornea. This can be somewhat blocked by eyelids,
resulting in some missing data near the outer edge of the image. The
boundary of the missing regions is indicated by the thin white curves.
The wild red and blue values are the consequence of extrapolation
in the Zernike fits.

\begin{figure}[H]
\includegraphics{ChapterO/lmtzcf21}

\caption{Zernike basis representation of 9 data object images in the Cornea
Curvature data from \citet{Locantore1999SphericalPCA}. Warmer colors
are stronger radial curvature. Note strong edge effects are quite
common. Circle plus signs denote the image (hence basis) centerpoint,
and white curves indicate boundaries of missing data. \label{figO.E:RawCorneaData}}
\end{figure}
As at many places in this book, PCA provides useful population level
insights through decomposition into modes of variation. The object\index{object space}
- feature space\index{feature space} terminology of Section \ref{secBA.1:Foundations_Terminology}
is useful here, as the images are the data objects, and each is represented
by a feature vector of $d=66$ Zernike coefficients. Hence, PCA is
done in $\mathbb{R}^{66}$, but the resulting modes of variation are
interpreted in terms of images. This is done in Figure \ref{figO.F:CorneaPCA},
which is from Figures 6-8 of \citet{Locantore1999SphericalPCA}, using
an approach similar to that of Figure \ref{fig:BladProsRect-PGA},
where the modes of variation are back-projected into the space of
the original objects. The top row shows the first mode of variation
with the sample mean in the center, the reconstructed image at two
standard deviations below the mean is shown on the left, and on the
right is two standard deviations above the mean. Note that the overall
color is hotter on the left and cooler on the right, which is consistent
with the first optometric measurement of overall curvature. The light
orange figure 8 pattern in the mean is a consequence of astigmatism,
which is a ridge of high curvature (the 8 shape is the radial curvature
view of such ridges). The astigmatism is stronger when there is more
overall curvature, indicating a correlation that was not previously
known by the clinician co-authors of \citet{Locantore1999SphericalPCA}. 

The second and third of these modes of variation are visualized on
the second and third rows of Figure \ref{figO.F:CorneaPCA}, using
the same format (mean between minus and plus two standard deviations).
Both the second and third mode have strong artifacts that look similar
to those in the bottom row of Figure \ref{figO.E:RawCorneaData}.
In particular, the case on the lower right appears to be impacting
the PC2 direction in $\mathbb{R}^{66}$, in a manner similar to that
illustrated in the right panel of Figure \ref{figO.A:Outliers}. This
was confirmed by looking at the PC2 scores plot, which looked quite
similar to the PC2 scores of the Parabolas and Outlier data in Figure
\ref{figO.C:10dToynOut}, verifying that this direction was indeed
driven by that single outlying case. A natural approach was to delete
that case, which indeed eliminated the impact of that outlier, but
then PC2 was dominated by another outlier (not surprising from the
bottom row of Figure \ref{figO.F:CorneaPCA}. Even after sequentially
deleting $4$ such outliers, there were still strong outlier effects
on the PCA, which is a concern because that is almost 10\% of the
data. That motivated a robust approach, i.e. an analog of PCA that
downweighted the effect of the outliers, while still leaving them
in the data set. This is also sensible since the central parts of
all images contain useful information.

\begin{figure}[H]
\includegraphics{ChapterO/lmtzcf33}

\includegraphics{ChapterO/lmtzcf34} 

\includegraphics{ChapterO/lmtzcf35}

\caption{Conventional PCA of the Cornea Curvature data from \citet{Locantore1999SphericalPCA}.
Rows are first 3 modes of variation, center column shows the sample
mean (same for each row). Left (right, respectively) column is the
mean minus (plus) 2 standard deviations. First mode is clinically
interesting. Second and third appear to be dominated by outliers,
motivating a robust approach. \label{figO.F:CorneaPCA}}
\end{figure}
As discussed in Section \ref{secP-A:PCA-Views} there are a number
of approaches to PCA. Most of these give can motivate approaches to
robust PCA. One is via the characterization of PCA as eigen analysis\index{eigen analysis|)}
of the sample covariance matrix (\ref{eqB.F:SampleCovariance}), with
the variance and covariance estimates replaced by robust versions.
A major challenge to this approach is that the resulting estimated
covariance is typically no longer positive definite, resulting in
hard to interpret negative eigenvalues. The approach of \citet{LiChen1985ProjPursPCA}
followed the idea of PCA as directions of maximal projected variation,
doing an iterative search over directions employing more robust measures
of spread on the projections. Local optima are a challenge to that
approach, and there do not seem to be successful implementations for
dimensions as large as $d=66$. Another intuitively appealing approach
is the minimum volume ellipsoid approach of \citet{RosseeuwLeroy1987},
but that has the drawback of requiring $d<n$, which is sensible for
affine invariance reasons, but not workable for this Cornea Curvature
data set with $d=66$ and $n=43$. This dearth of widely-known robust
PCA methods for high dimensional situations motivated \citet{Locantore1999SphericalPCA}
to invent spherical PCA.

While spherical PCA gave a much improved analysis, the outliers still
had a substantial influence. This was due to some of the Zernike coefficients
being orders of magnitude different from others (a common issue with
orthogonal basis representations). \citet{Locantore1999SphericalPCA}
tackled that problem by extending the spherical PCA idea illustrated
in Figure \ref{figO.B:OutliersnEsts}, to \emph{elliptical PCA}\index{elliptical PCA}.
The main idea is to replace projection onto the sphere by projection
onto an ellipse whose axes are parallel to the coordinate axes, thus
effectively handling the wildly different scales of the Zernike coefficients.
The elliptical PCA analysis of the Cornea Curvature data\index{Cornea Curvature data|)}
is shown in Figure \ref{figO.G:CorneaEllPCA} (from Figures 19-21
of \citet{Locantore1999SphericalPCA}), using the same format as Figure
\ref{figO.F:CorneaPCA}. Note that the first mode of variation is
quite similar to that of the first PC above, with the same interpretation,
but without the small influence of outliers. The second mode now clearly
appears as the known (to the clinicians) mode of variation of steeper
on top versus steeper on the bottom. The third mode is also clinically
interpretable in terms of direction of astigmatism. While astigmatism
can appear in any direction, most of it is approximately vertical,
called ``with the rule'' in optometry. Most of the rest is horizontal,
and is called ``against the rule''. The third mode is with the rule
versus against the rule variation, although only one lobe of the figure
8 in the lower right plot actually appears above the yellow to orange
color threshold.

\begin{figure}[H]
\includegraphics{ChapterO/lmtzcf56}

\includegraphics{ChapterO/lmtzcf57}

\includegraphics{ChapterO/lmtzcf58}

\caption{Elliptical PCA of the same data (and in the same format, again from
\citet{Locantore1999SphericalPCA}) as Figure \ref{figO.F:CorneaPCA}.
Shows outlier effects have been removed. All three modes of variation
are now clinically interpretable. \label{figO.G:CorneaEllPCA}}
\end{figure}
As noted in Section 6.11.2 of \citet{MaronnaMartinYohai2019}, spherical
PCA has since become a standard approach to robust PCA.\index{data objects!images|)}

\subsection{Case Study: Genome Wide Association Data\label{secO.B.B:GWAS}}

Genome Wide Association Studies (GWAS) were discussed in Section \ref{secM-B:HDLSS}.
As introduced in \citet{Klein2005FirstGWAS}, the data objects are
vectors representing many (up to 5 million with technology available
at this writing) local variants in the DNA of humans at genomic locations
called Single Nucleotide Proteins (SNPs). The entries of such vectors
are typically binary, yet as seen in Figure \ref{figO.H:GWAS-PCA},
which is part of Figure 1 in \citet{ZhouMarron2016VL1PCA}, PCA still
shows a large amount of interesting population structure. This is
an extreme example of how PCA is useful even for non-Gaussian data.
In fact these data objects all lie at vertices of a very high dimensional
unit cube.

Colors are used in Figure \ref{figO.H:GWAS-PCA} to represent some
ethnic groups, which are clearly important aspects of the variation
in this population. The data come from the cystic fibrosis analysis
done by \citet{Wright2011GWAS}. That paper did careful preprocessing\index{preprocessing}
resulting in binary SNP vectors of length $d=21,000$ (chosen to give
a representative set of SNPs). The data set is quite large, so for
convenient display only a subset of $n=347$ people chosen as described
in \citet{ZhouMarron2016VL1PCA} is shown here. While the ethnic groups
stand out strongly, note that some of the PC directions (PCs 3-6)
are clearly dominated by just a few outliers. The above robust PCA
ideas suggest analyzing this data set using spherical PCA, which is
done in Figure \ref{figO.I:GWAS-SphericalPCA}. 

\begin{figure}[H]
\includegraphics[scale=0.6]{ChapterO/GWAS-PCA}

\caption{Conventional $L^{2}$ PCA of a subset of the Cystic Fibrosis GWAS
data set, from \citet{ZhouMarron2016VL1PCA}. Shows major ethnic groups
in the first two modes of variation. Others are driven by pairs of
outliers. \label{figO.H:GWAS-PCA}}
\end{figure}
Based on the ideas in Figure \ref{figO.B:OutliersnEsts}, it is natural
to expect that the influence of the outliers would be greatly reduced
by using spherical PCA. Yet that clearly has not happened in Figure
\ref{figO.I:GWAS-SphericalPCA} (the rest of Figure 1 in \citet{ZhouMarron2016VL1PCA}),
which looks very similar to Figure \ref{figO.H:GWAS-PCA}, particularly
in terms of the impact of the outliers. This perhaps surprising behavior
can be understood using the HDLSS geometric representation\index{geometric representation}
ideas of Section \ref{secM-B:HDLSS}. Those concepts suggest that
all the data should be equidistant around the population expected
value, i.e. essentially lying on the surface of a sphere of radius
$d^{1/2}$. Furthermore pairwise angles (from the population center)
are all approximately $90^{\circ}$. A strong exception to this overall
pairwise orthogonality is first degree relatives (e.g. siblings or
parents and children) whom share about half of their SNPs. This results
in a pairwise angle of approximately $\cos^{-1}\left(\frac{1}{2}\right)=60^{\circ}$,
i.e. much smaller than the typical $90^{\circ}$, which is enough
to make just a pair of points drive a PC direction. It is not easy
to comprehend that the clearly visible ethnic arms visible in Figure
\ref{figO.H:GWAS-PCA} are the consequences of pairwise angles only
slightly smaller than $90^{\circ}$, caused by groups of people having
a larger number of SNPs in common. As these data objects already lie
near the surface of the sphere centered at the population expected
value, it is actually not surprising that the spherical PCA device
of projecting to a sphere has no meaningful impact, i.e. Figure \ref{figO.H:GWAS-PCA}
is so similar to Figure \ref{figO.I:GWAS-SphericalPCA}.

\begin{figure}[H]
\includegraphics[scale=0.6]{ChapterO/GWAS-SPCA}

\caption{Spherical PCA\index{Spherical PCA|)} of the same data (also from
\citet{ZhouMarron2016VL1PCA}) as in Figure \ref{figO.H:GWAS-PCA}.
Shows essentially no impact of projecting data onto a sphere and then
computing PCA. The strong influence\index{influence|)} of the outlying
pairs, and the perhaps surprising lack of robustness, is understood
using HDLSS ideas from Section \ref{secM-B:HDLSS}. \label{figO.I:GWAS-SphericalPCA}}
\end{figure}
An effective solution to the first degree relative outlier problem
in GWAS data illustrated in Figures \ref{figO.H:GWAS-PCA} and \ref{figO.I:GWAS-SphericalPCA}
is the Visual $L^{1}$ PCA proposed by \citet{ZhouMarron2016VL1PCA}.
That analysis is shown in Figure \ref{figO.J:GWAS-VL1PCA}. The first
two PC modes of variation clearly capture the ethnic groups visible
in the first two components of the earlier analyses. But the remaining
modes appear to be much more useful, as they shift the focus from
outlying pairs to separate groups of people that appear to represent
various ethnic subgroups. 

\begin{figure}[H]
\includegraphics[scale=0.6]{ChapterO/GWAS-VL1PCA} 

\caption{Visual $L^{1}$ analysis (also from \citet{ZhouMarron2016VL1PCA})
of the same GWAS data in Figures \ref{figO.H:GWAS-PCA} and \ref{figO.I:GWAS-SphericalPCA}.
This view is much more robust against the pairs of outlying first
degree relatives, instead shifting the focus to additional groups
of people (different ethnic subgroups). \label{figO.J:GWAS-VL1PCA}}
\end{figure}
\index{modes of variation|)}Visual $L^{1}$ PCA is based on the robust
$L^{1}$ PCA of \citet{BrooksDulaBoone2013LIPCA}. As noted in Section
\ref{secG-C:BackwardsPCA}, that method uses a clever backwards\index{backwards methods}
implementation. While $L^{1}$ PCA\index{Principal Component Analysis!and robustness|)}
gives a robust and useful set of directions for understanding population
variation, the corresponding $L^{1}$ scores tend to be quite hard
to interpret as noted in \citet{ZhouMarron2016VL1PCA}, due to the
lack of rotation invariance\index{rotation invariance} of the $L^{1}$
norm. The Visual $L^{1}$ approach is based on the same robust direction
vectors of maximal $L^{1}$ variation, but computes scores using the
more interpretable $L^{2}$ projections. This gives the needed insensitivity
to outliers demonstrated in Figure \ref{figO.J:GWAS-VL1PCA}, which
is part of Figure 6 in \citet{ZhouMarron2016VL1PCA}.

\section{Other Robustness Areas\label{secO.C:OtherRobust}}

As discussed in Section \ref{secE.1:Fr=0000E9chet-mean}, there are
many notions of multivariate median\index{median!multivariate}. These
all have the property that in the case $d=1$ they are just the conventional
sample median. Good access to this area can be found in \citet{Small1990Medians}
and \citet{ChakrbortyChaudhuri1999Medians}. Especially appealing
are the \emph{\index{data depth}data depth} methods, see e.g. \citet{Liu1999DataDepth},
\citet{VardiZhang2000Median} and \citet{LopezRomo2006DepthFDA}.
A median type estimate of scale\index{scale!MAD} is the Median\index{median!absolute deviation}
Absolute Deviation\index{Median Absolute Deviation} (from the median),
defined in (\ref{eqL.A:MAD}).

A novel approach to finding outliers in RNAseq data\index{RNAseq data},
of the type illustrated in Sections \ref{subsecB.2.A:LungCancerData}
and \ref{subsecL-B-A:HiDimSigClust} can be found in \citet{Choi2018RNAseqOutliers}.
See \citet{Ahn2019Outliers} for interesting outlier detection methods
in HDLSS\index{High Dimension Low Sample Size|)} contexts.

Some discussion of robustness in landmark\index{landmark} shape analysis
has been given by \citet[Section 13.6]{DrydenMardia2016}. In particular
distinctions are made between three different types of robustness:
resistance to landmark outliers (some of the $d$ variables are unusual
within a particular object); resistance to object outliers (some of
the $n$ objects are unusual in the sample); and robustness to model
mis-specification. \citet{Drydwalk99} discuss methods for landmark
outlier resistance and adapt the S-estimator (\citet{Rousyoha84,RosseeuwLeroy1987})
and the Least Median of Squares estimator (\citet{Rousseeuw84}) for
landmark shape matching. It is important to balance having a high
breakdown with high statistical efficiency, where \emph{\index{breakdown}breakdown}
is the minimum percentage of points that can be moved arbitrarily
to achieve maximum discrepancy. It is also crucial that any approach
to robust shape analysis be equivariant, for example the same match
should be obtained regardless of the arbitrary rotations of the objects.
\citet{Drydwalk99} used such an S-estimator with 25\% breakdown,
which has a good trade-off between breakdown and efficiency. Least
Median\index{median|)} of Squares has almost 50\% breakdown, and
an example of its use for matching electrophoresis gels is given by
\citet[Section 13.6]{DrydenMardia2016}. These resistant shape analysis
techniques can be helpful for identifying outlier landmarks, which
are highlighted by very large residuals after resistant matching.\index{robust methods|)} 
