
\chapter{Smoothing and SiZer \label{chapN:Smoothing}}

Smoothing\index{Smoothing|(} methods are useful approaches to obtaining
information from low dimensional (often just 1-d) sets of data. They
are often called \emph{nonparametric\index{nonparametric|(}} because
they do not rely on typical parametric models, instead aiming to produce
a useful visual interface for data analysis. The two main smoothing
contexts are \emph{\index{density estimation|(}density estimation}
and \emph{\index{nonparametric regression}nonparametric regression}. 

Density estimation starts with a random sample\index{random sample}
(recall independent and identically distributed from Section \ref{secE.1:Fr=0000E9chet-mean})
of scalars $\widetilde{x}_{1},\cdots,\widetilde{x}_{n}$ (using notation
from Table \ref{tab7-1:TildeNotation}) from an unknown probability
density\index{probability density|(} curve $f(x)$. The goal is to
recover $f$ which as noted in Section \ref{subsecBA.2.1:ProbDistns}
often provides useful population structure insights. The first thing
most people would try is a histogram, but as noted in Section \ref{secN-A:NoHisto}
that approach has problems that are not all that widely understood. 

Regression starts with a random sample of pairs $\left(\widetilde{x}_{1},\widetilde{y}_{1}\right),\cdots,\left(\widetilde{x}_{n},\widetilde{y}_{n}\right)$
and the goal is to recover the conditional expected value (i.e. regression
function) $m(x)=E\left[\widetilde{y}|\widetilde{x}=x\right]$. Again
the smoothing approach is nonparametric in the sense that instead
of imposing an assumed parametric structure by something in the spirit
of a classical linear model, the data are allowed to speak for themselves.
Smoothing methods provide a much more flexible view, which can lead
to unexpected discoveries that can be totally obscured by an inappropriate
rigid parametric structure. An elegant example is the discovery of
pre-pubertal growth spurts in human growth data, in a series of papers
summarized in \citet{Mueller1988NonparaReg}. 

In both cases, for the target to be identifiable, some assumptions
must be made. Such assumptions are almost always some variant of the
theme of \emph{smoothness}. This can be defined in terms of continuity
and differentiability, or in terms of various function spaces such
as Sobolev\index{Sobolev space} and Besov space\index{Besov space}s.

Perhaps because of the intuitive appeal and elegant simplicity of
smoothing, there is a very large literature. Important monographs
on density estimation include \citet{TapiaThompson1978} who point
to a biblical reference for early density estimation, \citet{DevroyeGyorfi1985L1DensityEst}
who present compelling reasons why theoretical analysis should be
based on the $L^{1}$ norm, \citet{Silverman1986Density} who gave
a very accessible early account, \citet{Scott2015MultivarDensity}
whose treatment is very comprehensive and \citet{Klemela2009Smoothing}
who proposed some novel high dimensional visualization\index{visualization}s.
Monographs on regression are more diverse and tend to advocate specific
methods, such as \citet{Mueller1988NonparaReg} and \citet{Hardle1990AppNonparaReg}
for two different kernel smoothers, \citet{FanGijbels1996LocPolyBook}
who made a strong case for local polynomial methods, \citet{Wahba1990Splines},
\citet{Eubank1999NonparaReg}, \citet{GreenSilverman1994}, \citet{Wang2011Splines}
and \citet{Gu2013SplineANOVA} who advocated smoothing splines, \citet{StoneEtAl1997Splines}
who make a compelling case for regression splines, and \citet{Gyorfi2002NonParaReg}
who take a distribution free approach. Additional accounts of smoothing
generally can be found in \citet{PrakasaRao1983}, \citet{ThompsonTapia1990},
\citet{WandJones1995} who provide a very accessible introduction
to kernel smoothing methods, \citet{Simonoff1996SmoothingBook}, \citet{Hart1997Tests},
\citet{Efromovich1999CurveEst} and \citet{schimek2013smoothing}.

\section{Why not Histograms? - Hidalgo Stamps Data\label{secN-A:NoHisto}}

When faced with a density estimation task (i.e. the need to visualize
the probability distribution underlying a set of numbers), most people
construct a histogram. Many users appreciate the fact that features
discovered in a histogram are impacted by choice of the binwidth.
That point is illustrated in Figure \ref{figN.A1:StampsHistosBinW}
showing various histograms of the 1872 \emph{Hidalgo Stamps} data\index{Hidalgo Stamps data|(}
using different binwidths. That data set was brought to the statistical
literature by \citet{Izenman1988StampsData}. It consists of $n=485$
paper thicknesses of the Hidalgo postage stamp. The question of interest
to philatelists is how many paper thicknesses appear in the data,
which is unknown as the relevant records have been lost. This question
is sensibly rephrased as: how many modes are in the density? 

Figure \ref{figN.A1:StampsHistosBinW} makes it clear this is not
a simple issue, and in fact each panel reflects node numbers that
have been suggested in the literature. In particular, at least two
modes was the conclusion using bootstrap ideas in Section 16.5 of
\citet{EfronTib1993Bootstrap}. Three modes were found by \citet{Walther2002Modes},
four were flagged as statistically significant by \citet{Sommerfeld2017TDASiZer},
3-5 by \citet{ChaudhuriMarron1999SiZer} and five by \citet{Minnotte2010ModeTest}.
Seven modes in this data were concluded by \citet{Izenman1988StampsData},
\citet{Basford1997Modelling} and \citet{FisherMarron2001ExMass}.
\citet{MinnotteScott1993ModeTree} suggested there might be up to
10 modes. 

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterN/OODAbookChpNFigA1StampsHistBW}

\caption{Histograms of the Hidalgo Stamps data for various binwidths. Shows
that the number of modes is strongly driven by choice of binwidth.\label{figN.A1:StampsHistosBinW} }
\end{figure}
Note that several different choices for the vertical axis are commonly
used. One is simple bin counts, another is the relative frequency,
while a third is to make the total area of the bars equal to 1. For
investigating a single data set, the choice usually does not matter
much as sensible scaling gives the same impression as long as the
bins all have the same width. However for comparing data sets of different
sizes, the second and third options are generally preferred. For density
estimation, it is only the third (used in all histograms here) that
results in an actual probability density. The difference becomes visually
important in the case of unequally spaced bins (for example equally
spaced in terms of distribution quantile\index{quantile}s) where
it is most interpretable to choose \emph{areas} to correspond to amounts
of probability mass because areas are most naturally perceived by
the human visual system. 

Figure \ref{figN.A1:StampsHistosBinW} makes it clear that too large
a histogram binwidth can result in missing important discoveries in
data, and a small binwidth appears to be strongly driven by random
sampling variation. However, as pointed out by Mats Rudemo, the human
perceptual system is rather good at visually smoothing. So in the
hands of an experienced analyst who understands the issue of spurious
small scale spikes, a small binwidth can lead to good results in terms
of data discovery. It may be worth keeping in mind that this notion
does \emph{not} seem to be a factor in the default histogram choice
of most common software packages. For example, as noted in Section
3.3 of \citet{Scott2015MultivarDensity}, the \citet{Sturges1926Histo}
method tends to grossly oversmooth. Nor does Scott's recommendation
give an undersmoothed binwidth because it targets the uni-modal Gaussian
distribution.

While many users understand the impact of the binwidth on features
discovered using histograms, as illustrated in Figure \ref{figN.A1:StampsHistosBinW},
less well understood is the \emph{\index{bin edge effect}bin edge
effect}, which was elegantly pointed out in Section 4.3 of \citet{Scott2015MultivarDensity}.
This is demonstrated in Figure \ref{figN.A2:StampsHistosLoc}, which
shows histograms of the Hidalgo Stamps data that both have the same
binwidth of 0.005 (same as in the lower left panel of Figure \ref{figN.A1:StampsHistosBinW}).
The obvious difference between these histograms is entirely driven
by horizontal shift, which clearly has a major impact on the number
of modes (from 6 in the left panel, to 2 on the right).

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterN/OODAbookChpNFigA2StampsHistLoc}

\caption{Histograms of Hidalgo Stamps both with binwidth 0.005, but shifted
to the right by half a bin in the right panel. Shows number of modes
can also be very sensitive to bin location.\label{figN.A2:StampsHistosLoc} }
\end{figure}
One approach to handling this bin edge effect is to average the histograms
over many shifts, resulting in the Average Shifted Histogram, studied
in Chapter 5 of \citet{Scott2015MultivarDensity}. This can be represented
as a variant of the very popular kernel density estimator\index{kernel density estimator|(}.
Given a \emph{kernel} function $K$ with integral one (e.g. a probability
density\index{probability density|)}), and a \emph{bandwidth\index{bandwidth}}
$h>0$,the Kernel Density Estimator (KDE) is
\begin{equation}
\hat{f}(x)=n^{-1}\sum_{i=1}^{n}K_{h}\left(x-\widetilde{x}_{i}\right),\label{eqN.A:KDE}
\end{equation}
where $K_{h}(\cdot)$ denotes the area preserving rescaling $\frac{1}{h}K\left(\frac{\cdot}{h}\right)$.
Three KDEs, showing how the bandwidth $h$ controls the smoothness,
for the Hidalgo Stamps data appear in Figure \ref{figN.A3:StampsKDE}.
As done many times in earlier chapters, the density estimates are
overlayed with the raw data as a \emph{\index{jitter plot|(}jitter
plot} (as discussed in Section \ref{secB2:Data-Visualization}). Note
that for this data set, the jitter plot gives a reasonable impression
of regions of high density. The three bandwidths were manually chosen
to give an impression of oversmoothing (red, $h=0.005,$showing just
two modes), undersmoothing (blue, $h=0.0006,$with many spurious modes)
and a perhaps good amount of smoothing (green, $h=0.0016,$with 7
modes as found in several of the above references). The small kernel
functions at the bottom of Figure \ref{figN.A3:StampsKDE} also give
an indication of how the KDE is constructed as a sum of shifts of
such functions. These kernels have been vertically rescaled so they
are all visible, but essentially each KDE is a sum of such functions
shifted to be centered at each data point.

\begin{figure}[H]
\includegraphics[scale=0.7]{ChapterN/OODAbookChpNFigA3StampsKDE}

\caption{Illustration of kernel density estimation. Raw Hidalgo Stamps data
are shown using a jitter plot. Colored curves for three different
bandwidths showing over (red), under (blue) and about right (green)
levels of smoothing. Corresponding vertically rescaled component kernels
appear at the bottom. \label{figN.A3:StampsKDE}}
\end{figure}
Figure \ref{figN.A3:StampsKDE} shows that choice of bandwidth is
critical, and more about this choice can be found in Sections \ref{secN-C:SmooParSel}
and \ref{secN-D:InfSiZer}. Here bandwidth\index{bandwidth!selection|(}
selection was done using what Bradley Davis has termed the \emph{Goldilocks
principle}\index{Goldilocks bandwidth selection|(}\emph{\index{}}.
Following a theme of the popular children's story the idea is to find
bandwidths that are both clearly under and over smoothed, which then
provide a context in which the one in the middle seems to be about
right. 

Less critical is the choice of kernel, although there are many views
on that as well. Computational issues were once an important consideration,
but much less so now with the ready availability of fast computation.
There have been several approaches to optimal kernel choice, with
the most sensible based on decoupling the bandwidth and kernel choices
using the \emph{\index{canonical kernel}canonical \index{kernel!canonical}kernel}
idea (see \citet{MarronNolan1988CanonicalKernels}). For actual data
analysis the Gaussian kernel is recommended, both because it tends
to produce the fewest distracting small scale visual artifacts, and
because of the \index{variation diminishing property|(}variation
diminishing property (it is the only kernel where the number of modes
is a decreasing function of the bandwidth) as discussed in \citet{ChaudhuriMarron2001ScaleSpace}.

The KDE can also give clear insight into the perhaps surprising behavior
of histograms highlighted in Figure \ref{figN.A2:StampsHistosLoc}.
This is done in Figure \ref{figN.A4:StampsHistKDE}, which simply
adds a KDE overlay. The peaks in the KDE clearly show why there are
so many modes in the histogram in the left panel (the peaks in the
KDE fall nicely within histogram bins) and why there tend to be so
few in the right panel (the peaks are split between bins). Since the
KDE so clearly diagnoses the bin edge problem, it follows that the
KDE (not the standard histogram) should be the method of first choice
for visualizing univariate distributions of data. This convention,
together with the jitter plot overlay as in Figure \ref{figN.A3:StampsKDE},
has been followed at most points in this monograph. An exception is
where each bin has a specific meaning, such as for the heat-map\index{heat-map}
color distributions shown in Figures \ref{figF.F1:ComboViewInputData}
- \ref{figF.F4:ComboViewMode2}.

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterN/OODAbookChpNFigA4StampsHistKDE}

\caption{Understanding of the histogram bin edge effect shown for the Hidalgo
Stamps data in Figure \ref{figN.A2:StampsHistosLoc}, by overlaying
a KDE (black curve). Shows many histogram modes when the KDE peaks
fall mostly within single bins (left panel) and few histogram modes
when the peaks are split between bins (right). \label{figN.A4:StampsHistKDE}}
\end{figure}
KDEs have already been extensively used in earlier chapters, most
notably in the scatterplot matrix\index{scatterplot matrix}\index{scores!scatterplot}
visualization\index{visualization}s introduced in Figure \ref{fig:50d-ToyPCAScatMat}
(for the Twin Arches data\index{Twin Arches data}), in the \index{marginal distribution plot}marginal
distribution plots described in Section \ref{secD1:Preprocess-MargDists},
and even as early as Figure \ref{fig:Mortality-PC1} (Spanish Mortality
data\index{Spanish Mortality data}). An issue deserving a little
more discussion is that of subdensities as introduced in the \index{Lung Cancer data}Lung
Cancer PCA example of Figure \ref{fig:LungCancerPCABrushed}. Figure
\ref{figN.F:SubDensities} studies an alternative to those sub-densities,
which is to use full densities (recall curves with area one underneath)
for each of the sub-populations of the Pan Cancer data\index{Pan Cancer data}.
The left panel is essentially a replot of the top left panel of Figure
\ref{figB-R:PanCanDWD}. Recall that view studied the gene expression
of 6 different cancer types (highlighted using colors and symbols)
compared in the direction of the projection on the DWD direction trained
on HNSC (Head and Neck Squamous cell Carcinoma) vs. KIRC (KIdney Renal
clear cell Carcinoma). Note that the left panel gives a clear impression
of how the subpopulations combine to form the full population (whose
density is the black curve). In particular, the blue KIRC group forms
really all of the left hand bump (even overlaying the black curve).
While the green HNSC is most of the bump on the right, there are clearly
a few magenta BLCA (Bladder Urothelial Carcinoma) cases mixed into
that bump as well. It is also visually clear how the other classes
make up the main bump of the black density. Plotting each colored
curve as a density itself is much less useful. One problem is that,
because the curves are narrower, the constraint of having area one
makes them much taller than the black curve. A simple fix is to scale
them to all have the same height, which is shown in the right panel
of Figure \ref{figN.F:SubDensities}. This gives a generally much
less intuitive impression of how the subgroups relate to the full
population. Exceptions to this principle include when some subgroups
are much smaller than others (so the curves can be hard to see), and
when it is desirable to think of subgroups in a relative way (independent
of subgroup sizes)

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterN/OODAbookChpNFigFSubDensities}

\caption{Comparison of sub-density plots (left panel) and subset density plots
(right), for the same Pan Cancer data as in the top left panel of
Figure \ref{figB-R:PanCanDWD}. Indicates sub-density scale\index{scale!sub-density}
is usually more intuitive. \label{figN.F:SubDensities}}
\end{figure}
A perhaps non-obvious issue in overlaying density estimates is that
use of a common \index{bandwidth!common}bandwidth $h$ can be very
important. A critical case is sub-density displays as in the left
panel of Figure \ref{figN.F:SubDensities} which are most interpretable
when the subdensities sum to the overall density (this happens naturally
for a common bandwidth). Due to the high noise level inherent to bandwidth
selection discussed in Section \ref{secN-C:SmooParSel}, this is generally
far from the case using separate default bandwidth selectors for each.
\citet{MarronSchmitz1992SimKDE} proposed use of a sensible common
bandwidth $h$ which is the \index{geometric mean|(}geometric mean
of the individually selected bandwidths. A reasonable alternative
(especially when subgroup sizes are very unbalanced) is to use the
union of the data to choose the common bandwidth.

Some other key ideas about smoothing are reviewed in Section \ref{secN-B:SmoothBasics}.
Section \ref{secN-C:SmooParSel} contains an overview of data driven
bandwidth selection. An alternate approach, based on scale space\index{scale space|(}
ideas, together with a fundamental approach to statistical inference
appears in Section \ref{secN-D:InfSiZer}.\index{kernel density estimator|)}

\section{Smoothing Basics - Bralower Fossils Data\label{secN-B:SmoothBasics} }

As discussed above the main contexts of statistical smoothing are
density estimation and nonparametric regression\index{nonparametric regression|(}.
An example of the latter is the \index{Bralower Fossils data|(}Bralower
Fossils data in Figure \ref{figN.B:FossilsLocLin}. The data set comes
from \citet{Bralower1997}, and was collected for the study of global
climate over very long time periods. Each plus sign represents a microscopic
fossil shell found in a core drilled into the floor of the ocean.
The shells are dated by the surrounding material, and the ages range
over about 30 million years, around a hundred million years ago. Strontium
is an element that naturally occurs in two different isotopes with
about 70\% prevalence of one. These shells are climatic indicators
because during an ice age, much of the sea water is frozen in the
ice caps, so chemicals in the ocean are more concentrated. When the
temperature warms, melting results in dilution. The shells absorb
the strontium isotopes at differing rates under varying dilutions,
which results in the observed data. This may sound like a serious
stretch of the data (e.g. the numbers on the vertical axis only differ
in the 4th decimal place), but note that there are clear systematic
patterns in the data as opposed to being driven by pure random noise.
Furthermore, there is not any natural choice of parametric model here
so it is appropriate to approach the data with smoothing techniques. 

The Bralower Fossil data are smoothed in Figure \ref{figN.B:FossilsLocLin}
using a local linear smoother\index{local linear smoother}, which
is generally not far from a simple moving average (except that a line
is fit to the data within each window). Each fit is weighted by a
Gaussian window function (playing the same role as the kernel $K$
in (\ref{eqN.A:KDE})). As noted in \citet{FanGijbels1996LocPolyBook},
one difference between the local linear and local average can be seen
at the edges, where these smooths track the local slope quite effectively,
in contrast to local averages which tend to inappropriately level
off at the edges. As for the KDE, the bandwidth plays a critical role,
and the Goldilocks approach to manual selection has again been used
here. The red curve, using $h=4$, seems clearly oversmoothed. In
particular it does not follow the data well through the very large
valley around 115 million years ago, but is also somewhat low at the
peak around 105 million years ago. The blue curve, for $h=0.4$, has
many wiggles that seem to follow fine scale noise driven variation,
and thus appears undersmoothed. This curve also reveals another property
of the local linear smoother: it can leave the vertical range of the
data, e.g. at around 96 million years ago. This is a consequence of
sparsity\index{sparsity} in the data together with tail properties
of the Gaussian kernel. In large horizontal gaps between data points,
the local linear smoother tends to follow the lines determined by
the two closest data points on each side, and to smoothly connect
those regions in the middle of the gap. The green curve, based on
$h=1.3$ may be about the right amount of smoothing, in terms of following
the overall trends in the data, while not responding to apparent measurement
error. Note that this value of $h=1.3$ is between $h=4$ and $h=0.4$
in a \emph{multiplicative} (not additive) sense, i.e. it is close
to the geometric mean. Generally smoothing parameters are best thought
of in a multiplicative way as discussed around Figure \ref{figN.D:FossilsSiZer}.
An important question is which of these observed features are ``really
there'', in the sense of statistical significance? This critical
inferential issue is addressed in Section \ref{secN-D:InfSiZer}.

\begin{figure}[H]
\includegraphics[scale=0.7]{ChapterN/OODAbookChpNFigBFossilsLocLin}

\caption{Local linear analysis of the Bralower Fossils data. Plus symbols are
ratios of strontium isotopes as a function of age. Three curves are
estimates of the regression function, with bandwidths that are oversmoothed
(red), undersmoothed (blue), and about right (green) for seeing important
underlying structure in the data, in a situation where no parametric
model is available. \label{figN.B:FossilsLocLin}}
\end{figure}
As noted in the above referenced monographs, many papers have been
written that analyze smoothing methods using asymptotics as the sample
size $n\rightarrow\infty$. Appealing lessons from them include insightful
quantification of how excessive sampling variation (i.e. wiggliness
of the curves) is a consequence of a small bandwidth, while bias which
appears as under-estimation of peaks and valleys is exacerbated when
the bandwidth is too large. When these effects are jointly studied,
it is natural to combine them into a trade-off resulting in an optimality
theory, which then characterizes smoothing methods in terms of rates
of convergence. Early on it was noticed that for a given amount of
smoothness, various methods had the same rate of convergence. This
motivated the idea that these rates may be the best possible, which
has been elegantly established using \emph{lower bound} theory\index{lower bound theory}.
Landmark papers of that type include \citet{Farrell1972}, \citet{Stone1982OptimalRates}
and \citet{Pinsker1980Optimal}. Particularly insightful was the geometric
view of optimal rates discovered by \citet{DonohoLiu1991GeomtrizingRates}.

An interesting twist on these results was started in the statistical
community by \citet{DonohoJohnstone1994Biometrika,Donoho1995WaveletAsymptopia},
who essentially showed that conventional smoothing methods had poor
spatial adaptivity properties. They went on to show that wavelet\index{wavelet|(}
thresholding provided an elegantly straightforward way to overcome
this problem. This generated a large amount of research, including
the Bayesian wavelet approach of \citet{KohnMarron2007BayesWavelets}.

An important and sometimes controversial aspect of such mathematical
analyses has been choice of error criteria. For nonparametric regression
the most commonly used criteria tend to involve squared error such
as the $L^{2}$ norm, which is natural for a number of reasons including
the fact that the regression function itself minimizes the theoretical
squared error. Another advantage of squared error measures is that
their Hilbert space inner product structure allows for straightforward
analysis, including direct quantification of the smoothing trade-off
in terms of variance and bias. Less clear is the case of density estimation.
As noted above, \citet{DevroyeGyorfi1985L1DensityEst} present a number
of compelling reasons why the $L^{1}$ norm is a more sensible criterion.
However that comes at substantial cost in terms of technical tractability,
which has meant that most discoveries such as performance of bandwidth
selectors have first been discovered using $L^{2}$ and then with
much more effort adapted to $L^{1}$. A weakness of most error criteria
is that their optimal bandwidths can be quite far from what human
experts would select as discussed around Figure \ref{figE.G:Wasserstein}.

While asymptotic results have been very popular and insightful, they
do leave open the question of how large a sample size is needed before
they take effect (some of which can be dubious such as the so-called
\emph{higher order kernels}\index{higher order kernels}). A number
of such questions were answered using \emph{exact risk analysis\index{exact risk analysis}}
in the context of univariate density estimation by \citet{MarronWand1992ExactRisk},
and in the multivariate case by \citet{JoesWand1993ExactRisk}. See
\citet{MarronAdak1998ExactRiskWavelet} for exact risk analysis of
wavelet\index{wavelet|)}s.

With the advent of massive modern computer resources, computational
time is no longer the consideration that it once was. However for
large data sets it can still sometimes be an issue worth considering.
In the days of slower computers, substantial effort was invested in
this issue. A path-breaking approach was Silverman's \index{Fourier!Fast Transform}\index{Fast Fourier Transform}Fast
Fourier Transform computation of\index{smoothing computation!FFT}
the KDE, see \citet{Silverman1982FFT}. Another very interesting and
effective approach was the \emph{\index{smoothing computation!updating}updating}
idea of \citet{Seifert1994FastComp}. \citet{FanMarron1994FastImpl}
analyzed all of these, and discovered that surprisingly the speed
of Silverman's method was more due to the \emph{\index{smoothing computation!binning}binning}
pre-processing than to the Fast Fourier Transform computation.

\section{Smoothing Parameter Selection\label{secN-C:SmooParSel}}

Given the importance of the bandwidth to what is seen in a KDE as
demonstrated using the Hidalgo Stamps data in Figure \ref{figN.A3:StampsKDE}
(parallel issues exist for all other types of smoothing method), it
is natural to believe that the fundamental issue in smoothing is selection
of the tuning parameter (the bandwidth $h$ for the KDE). That idea
led to a large amount of research on data based bandwidth selection.
For surveys of that area see \citet{Marron1989BWSelSurvey}, \citet{JonesMarronSheather1996BriefSurvey,JonesMarronSheather1996Progress},
and many of the monographs referenced above, as well as \citet{CaoCuevasManteiga1994BWSurvey}
and \citet{ParkTurlach1992BWsurvey}.

As is clear from several of the contributions to the discussion edited
by \citet{HardleSchimek1996}, there was never a consensus among researchers
in the area as to how smoothing parameters should be chosen. Using
terminology from \citet{JonesMarronSheather1996BriefSurvey}, opinions
were divided mostly between \emph{first generation} methods which
included familiar ideas such as cross-validation\index{cross-validation}
and AIC and more sophisticated \emph{second generation} methods, such
as the recommended \citet{SheatherJones1991PlugIn} Plug In method,
and the bootstrap approach of \citet{JonesMarronPark1991RootN}. 

Strong motivation for the development of second generation methods
came from deep asymptotic analysis\index{asymptotic analysis} revealing
a surprisingly large amount of sampling variation, first discovered
by \citet{HallMarron1987a,HallMarron1987b}. For example, using the
notation $\hat{h}$ for a first generation data driven bandwidth,
and $h_{0}$ as some type of optimal bandwidth such as a minimizer
of the Integrated Squared Error, $\intop\left[\hat{f}(x)-f(x)\right]^{2}dx,$the
bandwidth variation is quantified as 
\[
\frac{\hat{h}-h_{0}}{h_{0}}\sim n^{-1/10},
\]
in the limit as $n\rightarrow\infty$. Many variants of this were
later developed. This very slow rate of convergence was very discouraging,
yet did confirm an empirically known fact that methods such as cross-validation
occasionally result in unacceptably poor bandwidths, which has also
been observed in many other contexts. This also motivated the development
of second generation methods. However the latter sometimes required
very large sample sizes for their benefits to be felt, and seemed
to be less useful in other contexts such as spline smoothing. This
apparently led to the lack of consensus discussed above.

A perhaps surprising property of cross-validation and other first
generation methods is that they have a tendency to actually be \emph{negatively
correlated} with various optimal bandwidths, such as the minimizer
of the integrated squared error. See \citet{ChiuMarron1990NegCorr}
for a detailed analysis of this phenomenon.

For practical bandwidth choice, the recommended bandwidth selection
method depends on the context. In situations where sufficient time
and effort are available, the best approach is subjective visual trial
and error by an experienced analyst. When there is analysis time but
less expertise available, the Goldilocks approach suggested in Section
\ref{secN-A:NoHisto} (and as applied in Figure \ref{figN.B:FossilsLocLin})
is recommended. In situations where insufficient human resources can
be employed (e.g. when many estimates need to be constructed, such
as for the marginal distribution plot\index{marginal distribution plot}s
introduced in Section \ref{secD1:Preprocess-MargDists}), automatic
choice is needed. Following the recommendation of \citet{JonesMarronSheather1996BriefSurvey},
the plug-in method of \citet{SheatherJones1991PlugIn} is used here
as a default here for KDEs. Similarly the \citet{RuppertSheatherWand1995}
plug-in method is used here for local linear regression, when a default
is required.

Smoothing parameter selection issues are similarly important for the
kernel approaches to classification\index{classification!kernel}
discussed in Section \ref{secJ-B:KernelMethods}, especially in the
case of the radial basis function\index{radial basis function}, an
example of which is a Gaussian kernel. See \citet{Ahn2010RBFbandwidth}
and \citet{Oliva2016KernelLearn} for interesting work in that direction.

\section{Statistical Inference and SiZer\label{secN-D:InfSiZer}}

As\index{SiZer|(} noted in Section \ref{secN-B:SmoothBasics}, and
visually clear from Figures \ref{figN.A3:StampsKDE} (Hidalgo Stamps
data\index{Hidalgo Stamps data|)}) and \ref{figN.B:FossilsLocLin}
(Bralower Fossils data), a central issue for any actual data analysis
using smoothing methods is which observed phenomena (usually peaks
and / or valleys) are ``really there''. More precisely, which represent
important and reproducible underlying structure (thus real scientific
discoveries), as opposed to which are spurious artifacts of the sampling
variation. 

A natural approach for statisticians is to attempt to generalize the
idea of confidence interval for studying the uncertainty in a parameter
estimate to some type of confidence band\index{confidence band|(}s.
In the context of smoothing methods, this has three fundamental flaws:
\begin{itemize}
\item Variation in a bundle of curves is almost always far richer than can
be understood through simply studying bands. For example think of
how the rich modes of variation\index{mode of variation} displayed
through PCAs at many points in this monograph, such as Figures \ref{fig:Mortality-Raw}
- \ref{fig:Mortality-Scatterplot} (Spanish Mortality data\index{Spanish Mortality data})
and \ref{fig:LungCancerLogCurves} - \ref{fig:LungCancerLogCurvesBrushed}
(Lung Cancer data\index{Lung Cancer data}), would be totally missed
by simply putting bands around the set of curves in the original data. 
\item As noted in \citet{HardleMarron1991SimErrorBars}, a major challenge
to correct coverage of confidence bands in smoothing is correctly
modeling the inherent bias. Doug Nychka pointed out that if one could
adequately model this bias, one could just subtract it from the estimator
and then be back in the same position. While a number of alternate
approaches have been proposed none has been completely satisfactory.
Despite this, confidence bands are still widely used. It is recommended
that their inherent crudeness as an approximation be kept in mind
when interpreting these.
\item All existing confidence bands are based on the choice of the smoothing
parameter, which is problematic as discussed in Section \ref{secN-C:SmooParSel}.
This issue could be mitigated by taking the Goldilocks\index{Goldilocks bandwidth selection}
approach to bandwidth selection and applying confidence bands to all
three choices of bandwidth, although this would result in a busy graphic
that may be hard to interpret. That challenge could be handled with
multiple panels, although the preference in the literature seems to
be to ignore this issue.
\end{itemize}
These three challenges motivated \citet{ChaudhuriMarron1999SiZer,ChaudhuriMarron2001ScaleSpace}
to seek a completely different approach to statistical inference in
smoothing via the SiZer (SIgnificance of ZERo crossings) method. SiZer
tackles the bandwidth choice issue by taking a \emph{\index{scale space}scale
space} approach. As elegantly detailed in \citet{Lindeberg2013ScaleSpace},
scale space was a model for human visual perception\index{visual perception}
that was explored in the early days of the field of computer vision.
The human perceptual ability to perceive both macroscopic big picture
aspects as well as fine scale details in the same visual field, was
modeled as a family of Gaussian kernel smooths. In particular, the
large bandwidths put the focus on large scale features while the small
bandwidths focus on small scale details. Note that from this perspective
\index{bandwidth!selection|)}bandwidth selection (as discussed in
Section \ref{secN-C:SmooParSel}) is not a sensible pursuit, because
it is the full collection of curves that contain the useful information,
not any particular individual member.

\citet{Lindeberg2013ScaleSpace} was also very notable for providing
a detailed exploration of the \emph{\index{variation diminishing property|)}variation
diminishing property}. Stated in kernel smoothing terms, this idea
is the property that the number of modes in a smooth should be a decreasing
function of the bandwidth. As noted in \citet{ChaudhuriMarron2001ScaleSpace},
\citet{Lindeberg2013ScaleSpace} showed this property always holds
if and only if the kernel is Gaussian, in two different ways. One
approach is using \index{total positivity}total positivity ideas,
and the other is based on a \index{heat equation}heat equation representation
of smoothing. Parallel applications of heat equation ideas are discussed
in Section \ref{subsecG-A-A:CLTon Manifolds} in the context of central
limit theory on manifolds.

\subsection{Case Study: British Family Incomes Data\label{subsecN-D.A:BritishFamilyIncomes}}

The top panel of Figure \ref{figN.C:IncomesSiZer} shows a scale space
KDE analysis of the 1985 \index{British Family Incomes data|(}British
Family Incomes data from \citet{MarronSchmitz1992SimKDE}. As is common
in income analysis, inflationary aspects are controlled by dividing
by the sample mean. The sample size $n=7145$ after 56 cases with
scaled income larger than 3 have been truncated (allowing more attention
to be paid to the main body of the distribution). The green points
are a jitter plot of the raw data, highlighting a weakness of the
jitter plot: visual impression is strongly affected by sample size,
e.g. in this case population structure is obscured by too many data
points. The family of blue curves are KDEs with bandwidths logarithmically
equally spaced from $h=0.01$ to $h=1$. Logarithmic spacing is natural
when working with bandwidths (note the evenly spaced visual impression
among the blue curves) because bandwidths work as \emph{scale factors}.
For example one more naturally considers \emph{doubling} the bandwidth
than \emph{adding two} to the bandwidth. 

As expected, the KDEs indicate that the distribution of incomes is
strongly right skewed: most people have small to moderate incomes,
while relatively few have far larger incomes. However the smaller
bandwidths suggest an unexpected feature: there appear to be two modes
to this distribution. This was a potentially ground-breaking discovery,
because the standard parametric model for incomes (there are several
of these, all of which were strongly rejected by goodness of fit tests)
are all entirely unimodal. The existence of this bimodal structure
was confirmed in the PhD dissertation of Heinz-Peter Schmitz, by stratifying
the population so that the first peak corresponded to families on
fixed government incomes and the second was driven by three other
income categories, each of which were reasonably fit by standard income
distributions. The resulting mixture distribution then revealed that
the bimodal structure is actually an inherent aspect of the population.
\citet{MarronSchmitz1992SimKDE} studied this behavior over time,
revealing an interesting time trend in these peaks. As noted in Section
\ref{secN-A:NoHisto} that paper also pointed out that when using
multiple data-driven bandwidth choices the noise inherent to those
choices (described in Section \ref{secN-C:SmooParSel}) can be usefully
mitigated by using the common geometric mean\index{geometric mean|)}
of the selected bandwidths in each smooth. 

\begin{figure}[H]
\includegraphics[scale=0.6]{ChapterN/OODAbookChpNFigCIncomeSiZer}

\caption{SiZer analysis of British Family Incomes data. Reveals that both apparent
large bumps represent important population structure worthy of deeper
investigation, while smaller bumps may be only spurious sampling artifacts.
\label{figN.C:IncomesSiZer}}
\end{figure}
While this example highlights the potential for scientific success
using smoothing methods, it also reveals a potential risk. Had the
discovered feature been a sampling artifact (not reproducible in other
data sets), much analytical time and effort could be wasted. This
underlines once again the importance of confirmatory analysis, i.e.
of doing proper \emph{statistical inference} when using smoothing
methods to make discoveries using real data. \index{British Family Incomes data|)}

The scale space way of thinking also provides a natural solution to
the second drawback of confidence bands: the bias inherent to smoothing.
This is done by regarding the family of data smooths as targeting
not the true underlying curve ($f$ for density estimation\index{density estimation|)}
or $m$ for regression), but instead targeting the \emph{curve at
the level of resolution\index{level of resolution|(}} of the bandwidth.
The natural multi-resolution scale space representation of the underlying
curve is its convolution with the kernel at each bandwidth, which
is also the expected value of the estimator, so bias free estimates
are available. This comes at the price of a revised interpretation
of the goal of the analysis. That interpretation has not yet been
embraced by theoretical statisticians, however it is quite natural
for scientists who actually analyze data and are thus quite accustomed
to understanding data simultaneously at various levels of resolution.

The first flaw of confidence band\index{confidence band|)}s (that
bands only poorly reflect variation in curves) is addressed by another
shift in the focus of the analysis. The key idea is that most important
scientific discoveries by smoothing result from finding unexpected
peaks and / or valleys in data, so that should be the direct focus
of the inference. There is a mode testing literature, started by \citet{GoodGaskins1980BumpHunting},
and including \citet{Silverman1981MultimodalTest}, \citet{MuellerSawtzki1991EME},
\citet{FisherMarron2001ExMass} and a number of others. 

But the SiZer approach proposed by \citet{ChaudhuriMarron1999SiZer}
directly connects the statistical inference with the scale space visualization\index{visualization}.
The main idea is that a peak is determined by the curve going up on
one side, and down on the other. Hence SiZer focuses on the derivative,
and showing when it is significantly different from 0. Zero crossings
between both types of significant regions indicate significant peaks
and valleys. This type of statistical inference is shown in the SiZer
map in the bottom panel of Figure \ref{figN.C:IncomesSiZer}. Each
row corresponds to one scale of resolution (bandwidth), and colors
across each row flag statistical inference (in terms of a hypothesis
test whose null is a zero derivative). In particular, red is used
for significantly decreasing, blue for significantly increasing, and
the intermediate color of magenta when the hypothesis test is indeterminate.
One more color used is gray, which indicates not enough (essentially
$<5$) points in the Gaussian window for reliable inference. The top
rows of the bottom panel show blue (increasing) followed by red, so
both the up and the down are statistically significant. At the finest
scales (wiggliest blue curves) many peaks and valleys are present,
yet the corresponding rows of the SiZer map are almost all magenta,
indicating that those bumps can not be distinguished from spurious
sampling variation. The black line across the SiZer map indicates
the bandwidth $h=0.08$, chosen to focus attention on the apparent
bimodal structure, which also gives the thick black KDE among the
family in the top panel. That row shows blue - red - blue - red indicating
that both peaks are statistically significant. While the magenta indicates
the small scale bumps are spurious, a common question is: ``does
this analysis show there are one or two bumps in the underlying distribution?''
The answer is ``yes'', because from the scale space perspective
modality depends entirely on the level of resolution. At the coarsest
scales there is one mode, at medium scales two significant modes,
and at the finest scales many insignificant modes.

Of course careful attention needs to be paid to multiple comparison
issues because SiZer is based on many simultaneous hypothesis tests.
The first version accomplished this using some relatively crude approximations,
that resulted in somewhat anti-conservative inference. Much more precise
inference was developed using deeper probability theory in \citet{HannigMarron2006SiZer},
which is the current default and used in all examples here. 

Insight as to the amount of smoothing that is done at each scale,
i.e. row of the SiZer map, is provided by the dashed white curves.
These indicate $\pm2$ standard deviations of a centered Gaussian
kernel at each scale. 

The choice of red for down and blue for up is arbitrary. That scheme
makes intuitive sense for say economic applications. However, for
climatological applications where red is associated with warm and
blue with cooler trends, it makes sense to reverse the colors, as
done by \citet{HolmstromErasto2002Climate,ErastoHolmstrom2005MultiScale}. 

\subsection{Case Study: Bralower Fossils Data\label{subsecN-D.B:SiZerFossils}}

SiZer has also been developed for \index{nonparametric|)}nonparametric
regression\index{nonparametric regression|)}, as demonstrated in
Figure \ref{figN.D:FossilsSiZer}, which is a SiZer analysis of the
Bralower Fossils data from Section \ref{secN-B:SmoothBasics}. This
time at the coarsest scales the blue curves are not far from linear,
but are slightly steeper on the left. The top rows of the SiZer map
reveal that the decrease is statistically significant on the left
side. The level of resolution\index{level of resolution|)} $h=1.3$
is highlighted using the horizontal line in the bottom panel, and
the thick black curve in the top. The SiZer colors show a significant
increase on the left, and also flag the importance of the large valley
that bottoms out around 115 million years ago. While the corresponding
smooth green curve in Figure \ref{figN.B:FossilsLocLin} seductively
suggests a second valley around 98 million years ago, the gray in
the SiZer map indicates that the data in this region are too sparse
to draw any such conclusion.

\begin{figure}[H]
\includegraphics[scale=0.6]{ChapterN/OODAbookChpNFigDFossilsSiZer}

\caption{SiZer analysis of the Bralower Fossils data\index{Bralower Fossils data|)}
from Figure \ref{figN.B:FossilsLocLin}. Indicates the main valley
at around 115 million years ago is statistically significant, but
the smaller valley around 98 might be an artifact of random variation.
\label{figN.D:FossilsSiZer}}
\end{figure}

\subsection{Case Study: Mass Flux Data\label{subsecN-D.C:SiZerMassFlux}}

The Mass Flux data\index{Mass Flux data|(} set was introduced in
Figure \ref{figK.J:MassFluxPCA}. Recall that the PC1 scores indicated
3 bumps in the distribution, that suggested 3 clusters. A SiZer analysis
of the scores is shown in Figure \ref{figN.D:MassFluxSiZer}. As in
the above analyses the input data are shown as green in a jitter plot.
But this time instead of dots, green circles are used to be more consistent
with Figure \ref{figK.J:MassFluxPCA}. The blue-red patterns in the
SiZer map show that all three modes are strong underlying features
of the data, which is consistent with the three cloud types as discussed
in Section \ref{secK-C:VisualClust}. Note also that different bumps
actually are significant at different bandwidths (scales). In particular
the center bump is only really there at bandwidths around $h=10^{0.4}\approx2.5$
to $10^{0.6}\approx4.0$, while the right bump needs $h\approx10^{0.75}\approx5.6$
(highlighted with the horizontal black line). Thus this example highlights
how the inference done by SiZer is truly multi-scale in nature.

\begin{figure}[H]
\includegraphics[scale=0.6]{ChapterN/OODAbookChpNFigEMassFluxSiZer}

\caption{SiZer analysis of the PC1 scores for the Mass Flux data. Shows three
bumps are important underlying features that are significant at different
scales. \label{figN.D:MassFluxSiZer}}
\end{figure}
Of course clustering\index{clustering!Mass Flux data} in the Mass
Flux data\index{Mass Flux data|)} could also be analyzed using the
formal clustering methods described in Chapter \ref{chapK:Clustering}.
Both 3-means clustering and also hierarchical clustering\index{hierarchical clustering}
(Euclidean distance\index{distance!Euclidean} with either average
or Ward's linkage\index{linkage!Ward's}\index{linkage!average})
gave very similar results. Consistent with the SiZer analysis in Figure
\ref{figN.D:MassFluxSiZer}, these are statistically significant using
the SigClust\index{SigClust|(} method described in Section \ref{secL2:SigClust}.

\subsection{Case Study: Kidney Cancer Data\label{subsecN-D.E:KidneySiZer}}

The Kidney Cancer part of the \index{Pan Cancer data}Pan Cancer data
introduced in Section \ref{subsecB.2.B:PanCancerExample} illustrates
a variation of SiZer. Unlike the analysis there (and also in Section
\ref{subsecL-B-A:HiDimSigClust}), instead of considering only a subset
of $n=50$, here the full $n=551$ data set of KIRC cases is analyzed.
The first mode of variation\index{mode of variation} scores from
the PCA of that data set are shown as green dots in the top panel
of Figure \ref{figN.G:KidneyCanSiZer}. The pattern in these dots,
as well as the overlaid kernel density estimates, suggest an apparent
cluster on the right which is a potential distinct subtype of Kidney
Cancer. The importance of cancer subtypes discovered by finding clusters
in gene expression was demonstrated by \citet{perou2000molecular},
who discovered clinically relevant subtypes that since motivated the
development of different treatment regimes.

\begin{figure}[H]
\includegraphics[scale=0.6]{ChapterN/OODAbookChpNFigGKidneyCanSiZer}

\caption{Slope and curvature \index{SiZer!curvature}SiZer analyses of PC1
scores of the full Kidney Cancer data. Jitter plot\index{Jitter plot|)}
and family of kernel density estimates in the top panel suggest bimodality.
Curvature analysis (bottom panel) confirms statistical significance
of bumps, suggesting subtypes that are not flagged as significant
using the slope version (middle panel). \label{figN.G:KidneyCanSiZer} }
\end{figure}
The middle panel of Figure \ref{figN.G:KidneyCanSiZer} investigates
this potential subtype using a Sizer analysis. The large peak in the
area of PC1 scores around -20 is clearly indicated as statistically
significant by the nearby blue and red regions in the SiZer map. However,
the smaller peak suggested by the cluster of PC1 scores near 100 is
not indicated as statistically significant. The reason is that the
density is generally sloping downwards, and the cluster is not large
enough for a statistically significant upwards slope to appear.

The bottom panel of Figure \ref{figN.G:KidneyCanSiZer} addresses
this challenge by modifying the conventional SiZer analysis, from
inference about slopes (first derivatives) to a parallel inference
based on statistical significance of curvature (second derivatives),
in an analysis first done by Tom Keefe. In particular, using the method
proposed in \citet{Chaudhuri2002CurvatureSiZer}, regions in the scale
space map (with location on the horizontal axis and bandwidth\index{bandwidth|)}
scale, i.e. $\log_{10}h$, on the vertical axis) of significant concavity
($f''<0$) are colored cyan, regions of significant convexity ($f''>0$)
are orange shaded, and inconclusive regions are colored green. The
horizontal black line at the scale $h=20$ passes through two substantial
cyan regions, leading to the conclusion of two significant bumps (hence
subtypes) in the Kidney Cancer data, despite its on-going downward
tendency on the right side of the main peak. 

This result is consistent with the SigClust discovery of two clusters
in a random subsample of $n=50$ of these data points in Section \ref{subsecL-B-A:HiDimSigClust}.
An important point is that since SigClust\index{SigClust|)} directly
targets modes it is a more powerful test of modality. In particular
it found the subtypes in the Kidney Cancer data with only a small
subset ($n=50$) of the full data ($n=551)$. Subtypes in a larger
set of KIRC data from TCGA have been studied in detail by \citet{Ricketts2018TCGAKIRC}.\index{kernel density estimate|)}

\subsection{Additional SiZer Applications and Variants\label{subsecN-D.D:AdditionalSiZer}}

There are a number of variants of Sizer, for example the inference
can be based on the 2nd derivative, instead of the 1st (i.e. curvature
not slope). The curvature version of SiZer found important population
structure (modality) that was not available from the slope version
in the analysis of tree\index{tree} structured data objects\index{data objects!tree structured}
of \citet{ShenTrees2014}. SiZer has been modified to handle the specific
case of jumps (i.e. change-points) by \citet{KimMarron2006JumpSiZer}.
\citet{Li2005LocLikSiZer} developed a local likelihood version of
SiZer. A smoothing spline version can be found in \citet{MarronZhang2005SiZerSplines}.
A wavelet\index{wavelet} version was invented by \citet{Park2007WaveletSiZer}.
An interesting application to establish the existence of \emph{code
decay} in software engineering can be found in \citet{Eick2001CodeDecay}.
Analysis of internet traffic data motivated variations of SiZer that
can handle dependent data, see \citet{Hannig2001Zooming}, \citet{Park2004DepSiZer}
and \citet{Rondonotti2007DepSiZer}.

SiZer type ideas have also been adapted to the two dimensional case,
which thus includes significance of features in images\index{images!significant features}
to give a methodology called \emph{\index{Statistical Significance in Scale-space}Statistical
Significance in Scale-space} by \citet{GodtliebsenMarron2002SSSdensity,GodtliebsenMarron2004SSSimage}.
The underlying inferential techniques are basically the same, but
the visual interface is much different, mostly because of the challenges
of extending the basic notions of increasing and decreasing to more
than one dimension.

The scale space\index{scale space|)} ideas that underpin SiZer are
not unrelated to the persistent homology filtration\index{filtration}s
that are the basis of Topological Data Analysis\index{Topological Data Analysis}
as discussed in Section \ref{subsecI-A-D:Persistent-Homology}. This
connection has been studied in \citet{Sommerfeld2017TDASiZer}, who
use it to propose an alternate approach to inference for bumps, based
on the bootstrap ideas of \citet{Fasy2014TDABootstrap}, that is seen
to be slightly more sensitive than SiZer\index{SiZer|)}.\index{smoothing}
