\HeaderA{silhouette}{Compute or Extract Silhouette Information from Clustering}{silhouette}
\aliasA{plot.silhouette}{silhouette}{plot.silhouette}
\aliasA{print.summary.silhouette}{silhouette}{print.summary.silhouette}
\methaliasA{silhouette.clara}{silhouette}{silhouette.clara}
\methaliasA{silhouette.default}{silhouette}{silhouette.default}
\methaliasA{silhouette.partition}{silhouette}{silhouette.partition}
\aliasA{sortSilhouette}{silhouette}{sortSilhouette}
\aliasA{summary.silhouette}{silhouette}{summary.silhouette}
\keyword{cluster}{silhouette}
\begin{Description}\relax
Compute silhouette information according to a given clustering in
\eqn{k}{} clusters.
\end{Description}
\begin{Usage}
\begin{verbatim}
silhouette(x, ...)
## Default S3 method:
silhouette  (x, dist, dmatrix, ...)
## S3 method for class 'partition':
silhouette(x, ...)
## S3 method for class 'clara':
silhouette(x, full = FALSE, ...)

sortSilhouette(object, ...)
## S3 method for class 'silhouette':
summary(object, FUN = mean, ...)
## S3 method for class 'silhouette':
plot(x, nmax.lab = 40, max.strlen = 5,
     main = NULL, sub = NULL, xlab = expression("Silhouette width "* s[i]),
     col = "gray",  do.col.sort = length(col) > 1, border = 0,
     cex.names = par("cex.axis"), do.n.k = TRUE, do.clus.stat = TRUE, ...)
\end{verbatim}
\end{Usage}
\begin{Arguments}
\begin{ldescription}
\item[\code{x}] an object of appropriate class; for the \code{default}
method an integer vector with \eqn{k}{} different integer cluster
codes or a list with such an \code{x\$clustering}
component.  Note that silhouette statistics are only defined if
\eqn{2 \le k \le n-1}{2 <= k <= n-1}.
\item[\code{dist}] a dissimilarity object inheriting from class
\code{\LinkA{dist}{dist}} or coercible to one.  If not specified,
\code{dmatrix} must be.
\item[\code{dmatrix}] a symmetric dissimilarity matrix (\eqn{n \times n}{n * n}),
specified instead of \code{dist}, which can be more efficient.
\item[\code{full}] logical specifying if a \emph{full} silhouette should be
computed for \code{\LinkA{clara}{clara}} object.  Note that this requires
\eqn{O(n^2)}{} memory, since the full dissimilarity (see
\code{\LinkA{daisy}{daisy}}) is needed internally.
\item[\code{object}] an object of class \code{silhouette}.
\item[\code{...}] further arguments passed to and from methods.
\item[\code{FUN}] function used to summarize silhouette widths.
\item[\code{nmax.lab}] integer indicating the number of labels which is
considered too large for single-name labeling the silhouette plot.
\item[\code{max.strlen}] positive integer giving the length to which
strings are truncated in silhouette plot labeling.
\item[\code{main, sub, xlab}] arguments to \code{\LinkA{title}{title}}; have a
sensible non-NULL default here.
\item[\code{col, border, cex.names}] arguments passed
\code{\LinkA{barplot}{barplot}()}; note that the default used to be \code{col
      = heat.colors(n), border = par("fg")} instead.\\
\code{col} can also be a color vector of length \eqn{k}{} for
clusterwise coloring, see also \code{do.col.sort}:

\item[\code{do.col.sort}] logical indicating if the colors \code{col} should
be sorted \dQuote{along} the silhouette; this is useful for casewise or
clusterwise coloring.
\item[\code{do.n.k}] logical indicating if \eqn{n}{} and \eqn{k}{} \dQuote{title text}
should be written.
\item[\code{do.clus.stat}] logical indicating if cluster size and averages
should be written right to the silhouettes.
\end{ldescription}
\end{Arguments}
\begin{Details}\relax
For each observation i, the \emph{silhouette width} \eqn{s(i)}{} is
defined as follows: \\
Put a(i) = average dissimilarity between i and all other points of the
cluster to which i belongs (if i is the \emph{only} observation in
its cluster, \eqn{s(i) := 0}{} without further calculations).
For all \emph{other} clusters C, put \eqn{d(i,C)}{} = average
dissimilarity of i to all observations of C.  The smallest of these
\eqn{d(i,C)}{} is \eqn{b(i) := \min_C d(i,C)}{},
and can be seen as the dissimilarity between i and its \dQuote{neighbor}
cluster, i.e., the nearest one to which it does \emph{not} belong.
Finally, \deqn{s(i) := \frac{b(i) - a(i) }{max(a(i), b(i))}.}{s(i) := ( b(i) - a(i) ) / max( a(i), b(i) ).}

\code{silhouette.default()} is now based on C code donated by Romain
Francois (the R version being still available as
\code{cluster:::silhouette.default.R}).

Observations with a large \eqn{s(i)}{} (almost 1) are very well
clustered, a small \eqn{s(i)}{} (around 0) means that the observation
lies between two clusters, and observations with a negative
\eqn{s(i)}{} are probably placed in the wrong cluster.
\end{Details}
\begin{Value}
\code{silhouette()} returns an object, \code{sil}, of class
\code{silhouette} which is an [n x 3] matrix with attributes.  For
each observation i, \code{sil[i,]} contains the cluster to which i
belongs as well as the neighbor cluster of i (the cluster, not
containing i, for which the average dissimilarity between its
observations and i is minimal), and the silhouette width \eqn{s(i)}{} of
the observation.  The \code{\LinkA{colnames}{colnames}} correspondingly are
\code{c("cluster", "neighbor", "sil\_width")}.

\code{summary(sil)} returns an object of class
\code{summary.silhouette}, a list with components
\begin{ldescription}
\item[\code{si.summary}] numerical \code{\LinkA{summary}{summary}} of the individual
silhouette widths \eqn{s(i)}{}.
\item[\code{clus.avg.widths}] numeric (rank 1) array of clusterwise
\emph{means} of silhouette widths where \code{mean = FUN} is used.
\item[\code{avg.width}] the total mean \code{FUN(s)} where \code{s} are the
individual silhouette widths.
\item[\code{clus.sizes}] \code{\LinkA{table}{table}} of the \eqn{k}{} cluster sizes.
\item[\code{call}] if available, the call creating \code{sil}.
\item[\code{Ordered}] logical identical to \code{attr(sil, "Ordered")}, see
below.
\end{ldescription}


\code{sortSilhouette(sil)} orders the rows of \code{sil} as in the
silhouette plot, by cluster (increasingly) and decreasing silhouette
width \eqn{s(i)}{}.
\\
\code{attr(sil, "Ordered")} is a logical indicating if \code{sil} \emph{is}
ordered as by \code{sortSilhouette()}. In that case,
\code{rownames(sil)} will contain case labels or numbers, and \\
\code{attr(sil, "iOrd")} the ordering index vector.
\end{Value}
\begin{Note}\relax
While \code{silhouette()} is \emph{intrinsic} to the
\code{\LinkA{partition}{partition}} clusterings, and hence has a (trivial) method
for these, it is straightforward to get silhouettes from hierarchical
clusterings from \code{silhouette.default()} with
\code{\LinkA{cutree}{cutree}()} and distance as input.

By default, for \code{\LinkA{clara}{clara}()} partitions, the silhouette is
just for the best random \emph{subset} used.  Use \code{full = TRUE}
to compute (and later possibly plot) the full silhouette.
\end{Note}
\begin{References}\relax
Rousseeuw, P.J. (1987)
Silhouettes: A graphical aid to the interpretation and validation of
cluster analysis. \emph{J. Comput. Appl. Math.}, \bold{20}, 53--65.

chapter 2 of Kaufman, L. and Rousseeuw, P.J. (1990), see
the references in \code{\LinkA{plot.agnes}{plot.agnes}}.
\end{References}
\begin{SeeAlso}\relax
\code{\LinkA{partition.object}{partition.object}}, \code{\LinkA{plot.partition}{plot.partition}}.
\end{SeeAlso}
\begin{Examples}
\begin{ExampleCode}
data(ruspini)
pr4 <- pam(ruspini, 4)
str(si <- silhouette(pr4))
(ssi <- summary(si))
plot(si) # silhouette plot

si2 <- silhouette(pr4$clustering, dist(ruspini, "canberra"))
summary(si2) # has small values: "canberra"'s fault
plot(si2, nmax= 80, cex.names=0.6)

op <- par(mfrow= c(3,2), oma= c(0,0, 3, 0),
          mgp= c(1.6,.8,0), mar= .1+c(4,2,2,2))
for(k in 2:6)
   plot(silhouette(pam(ruspini, k=k)), main = paste("k = ",k), do.n.k=FALSE)
mtext("PAM(Ruspini) as in Kaufman & Rousseeuw, p.101",
      outer = TRUE, font = par("font.main"), cex = par("cex.main"))
par(op)

## clara(): standard silhouette is just for the best random subset
data(xclara)
set.seed(7)
str(xc1k <- xclara[sample(nrow(xclara), size = 1000) ,])
cl3 <- clara(xc1k, 3)
plot(silhouette(cl3))# only of the "best" subset of 46
## The full silhouette: internally needs large (36 MB) dist object:
sf <- silhouette(cl3, full = TRUE) ## this is the same as
s.full <- silhouette(cl3$clustering, daisy(xc1k))
if(paste(R.version$major, R.version$minor, sep=".") >= "2.3.0")
   stopifnot(all.equal(sf, s.full, check.attributes = FALSE, tol = 0))
## color dependent on original "3 groups of each 1000":
plot(sf, col = 2+ as.integer(names(cl3$clustering) ) %/% 1000,
     main ="plot(silhouette(clara(.), full = TRUE))")

## Silhouette for a hierarchical clustering:
ar <- agnes(ruspini)
si3 <- silhouette(cutree(ar, k = 5), # k = 4 gave the same as pam() above
                   daisy(ruspini))
plot(si3, nmax = 80, cex.names = 0.5)
## 2 groups: Agnes() wasn't too good:
si4 <- silhouette(cutree(ar, k = 2), daisy(ruspini))
plot(si4, nmax = 80, cex.names = 0.5)
\end{ExampleCode}
\end{Examples}

