---
title: "Using Self-Organizing Maps with SOMbrero to cluster a numeric dataset"
---

<!--
%\VignetteEngine{knitr}
%\VignetteIndexEntry{SOM for numeric data}
%\VignettePackage{SOMbrero}
-->

## Basic package description

To be able to run the SOM algorithm, you have to load the package called 
```SOMbrero```. The function used to run it is called ```trainSOM()``` and is 
detailed below.

__This documentation only considers the case of numerical data.__

```{r loading, results='hide', echo=FALSE, warning=FALSE, message=FALSE}
library(SOMbrero)
```

### Arguments

The ```trainSOM``` function has several arguments, but only the first one is
required. This argument is ```x.data``` which is the dataset used to train the 
SOM. In this documentation, it is passed to the function as a matrix or a data
frame with numeric variables in columns and observations of these variables in
rows.

The other arguments are the same as the arguments passed to the ```initSOM```
function (they are parameters defining the algorithm, see ```help(initSOM)```
for further details).

### Outputs

The ```trainSOM``` function returns an object of class ```somRes``` (see 
```help(trainSOM)``` for further details on this class).

### Graphics

The following table indicates which graphics are available for a numeric SOM.

Type        |   Energy  |   Obs |   Prototypes  |   Add   |   Super Cluster
------------|:---------:|:-----:|:-------------:|:-------:|:----------------:
no type     |   x       |       |               |         |
hitmap      |           |   x   |               |         |         x
color       |           |   x   |       x       |   x     |         x2
lines       |           |   x   |       x       |   x     |         x2
barplot     |           |   x   |       x       |   x     |         x2
radar       |           |   x   |       x       |   x     |         x2
pie         |           |       |               |   x     |         x
boxplot     |           |   x   |               |   x     |         x2
3d          |           |       |       x       |         |
poly.dist   |           |       |       x       |         |         x
umatrix     |           |       |       x       |         |
smooth.dist |           |       |       x       |         |
words       |           |       |               |   x     |
names       |           |   x   |               |   x     |
graph       |           |       |               |   x     |         x
mds         |           |       |       x       |         |         x
grid.dist   |           |       |       x       |         |         
grid        |           |       |               |         |         x
dendrogram  |           |       |               |         |         x
dendro3d    |           |       |               |         |         x

The plots marked by "x2", in the "Super Cluster" column, are available for 
both data set variables and additional variables. 

## First case study: simulated data in $[0,1]^2$

The first case study show the clustering of points randomly distributed in the
square $[0,1]^2$. The data are generated by:
```{r dataGeneration}
set.seed(4031719)
the.data <- data.frame("x1"=runif(500), "x2"=runif(500))
plot(the.data, pch=19)
```

### Training the SOM
The numeric SOM algorithm is used to cluster the data:
```{r dataTrain}
set.seed(593)
# run the SOM algorithm with 10 intermediate backups and 2000 iterations
my.som <- trainSOM(x.data=the.data, dimension=c(5,5), nb.save=10, maxit=2000, 
                   scaling="none", radius.type="letremy")
```

The energy evolves as described in the following graphic:
```{r energy}
plot(my.som, what="energy")
```

It is stabilized during the last 500 iterations.

### Clustering

The resulting clustering distribution can be visualized by the hitmap:
```{r hitmapObs}
plot(my.som, what="obs", type="hitmap")
```

The observations are almost uniformly distributed on the map.

The clustering component allows us to plot the initial data according to the 
final clustering.

```{r clusteredData, echo=FALSE, cache=TRUE}
# prepare a vector of colors
my.colors <- rainbow(prod(my.som$parameters$the.grid$dim))[my.som$clustering]

# points depicted with the same color are in the same final cluster
plot(my.som$data[,1], my.som$data[,2], col=my.colors, pch=19, xlab="x1", 
     ylab="x2", main="Data according to final clustering")
```

### Clustering interpretation

The values of the prototypes can be represented with the plot function and help
interpret the clusters:
```{r colorProto, fig.width=5, fig.height=2.5}
par(mfrow=c(1,2))
plot(my.som, what="prototypes", type="color", var=1, main="prototypes - x1")
plot(my.som, what="prototypes", type="color", var=2, main="prototypes - x2")
```

Here, the interpretation is simple enough: high values of the first variables x1
are located at the right side of the map and small values at the left side of 
the map. Large values of x2 are located at the top of the map, whereas, small
values are located at the bottom of the map.

We obtain the same results with a similar plot on the observation mean values:
```{r colorObs, fig.width=5, fig.height=2.5}
par(mfrow=c(1,2))
plot(my.som, what="obs", type="color", var=1, main="obs mean values - x1")
plot(my.som, what="obs", type="color", var=2, main="obs mean values - x2")
```

The prototypes coordinates are also registered for each intermediate backup so 
they can be displayed on different graphics to see the evolution in the
prototypes organization.
```{r protoEvoluation, fig.width=15, fig.height=6, echo=FALSE}
# find out the prototypes to be linked
tra <- NULL
for (i in c(1,6,11,16,21)){
  tra <- c(tra, i, i+1, i+1, i+2, i+2, i+3, i+3, i+4)
}
for (i in c(1:5)){
  tra <- c(tra, i, i+5, i+5, i+10, i+10, i+15, i+15, i+20)
}
tmp <- matrix(tra, ncol=2, byrow=TRUE)
# plot the prototypes
par(mfrow=c(2, 5),mar=c(3,2,2,1))
invisible(sapply(1:my.som$parameters$nb.save, function(ind){
  plot(my.som$backup$prototypes[[ind]][,1], my.som$backup$prototypes[[ind]][,2],
       xlab="", ylab="", main=c("iteration ", my.som$backup$steps[ind]))
  for (i in 1:nrow(tmp)){
    segments(x0=my.som$backup$prototypes[[ind]][tmp[i,1],1], 
             y0=my.som$backup$prototypes[[ind]][tmp[i,1],2],
             x1=my.som$backup$prototypes[[ind]][tmp[i,2],1], 
             y1=my.som$backup$prototypes[[ind]][tmp[i,2],2], 
             col="red", pch=19)
  }
}))
```

At the begining of the algorithm, the prototypes are randomly distributed in
[0,1]^2 and then, they organize as a regular rectangular grid in $[0,1]^2$.

## Second case study: the iris dataset

This second case study is performed on the famous (Fisher's or Anderson's) iris
data set that gives the measurements in centimeters of the variables sepal
length and width and petal length and width, respectively, for 50 flowers from
each of 3 species of iris (setosa, versicolor, and virginica).

### Training the SOM

The first four variables of the data set (that are the numeric variables) are 
used to map each flower on the SOM grid.
```{r irisTrain, cache=TRUE}
set.seed(255)
# run the SOM algorithm with verbose set to TRUE
iris.som <- trainSOM(x.data=iris[,1:4], verbose=TRUE, nb.save=5)
iris.som
```

As the energy is registered during the intermediate backups, we can have a look
at its evolution.
```{r energyIris}
plot(iris.som, what="energy")
```

Here the energy does not stabilize as in the case of ```dist.type="letremy"```
because the Gaussian annealing of the neighbourhood is continuous and not 
stepwise.


### Resulting clustering

The clustering component contains the final classification of the dataset. It is 
a vector with length equal to the number of rows of the input dataset.
```{r irisClusters}
iris.som$clustering
table(iris.som$clustering)
```
which can also be visualized by a hitmap plot:
```{r irisHitmap}
plot(iris.som, what="obs", type="hitmap")
```

To access the relevance of each explanatory variable in the definition of the 
clusters, the function ```summary``` includes an ANOVA with the factor being the
clustering, for each (numeric) input variable.
```{r irisSummary}
summary(iris.som)
```
Here, all variables have significantly different means among the different
clusters and can thus be considered to be relevant for the clustering
definition.

Another useful function is ```predict.somRes```. This function predicts the 
neuron to which a new observation would be assigned. The first argument must be 
a ```somRes``` object and the second one the new observation.
Let us have a try on the first observation of the iris data set:
```{r irisPred1}
# call predict.somRes
predict(iris.som, iris[1,1:4])
# check the result of the final clustering with the SOM algorithm
iris.som$clustering[1]
```

### Clustering interpretation

#### Graphics common to observations and prototypes

__NB: in all the following graphics, data are scaled in the preprocessing
stage.__

Some graphics are shared between observations and prototypes. They allow you 
to plot the level values, either of the prototypes or of the observation means. 
In this example these common graphics are plotted for the mean observation 
values.

```{r irisGraphOP}
par(mfrow=c(2,2))
plot(iris.som, what="obs", type="color", variable=1, print.title=TRUE, 
     main="Sepal length")
plot(iris.som, what="obs", type="color", variable=2, print.title=TRUE, 
     main="Sepal width")
plot(iris.som, what="obs", type="color", variable=3, print.title=TRUE, 
     main="Petal length")
plot(iris.som, what="obs", type="color", variable=4, print.title=TRUE, 
     main="Petal width")
plot(iris.som, what="prototypes", type="lines", print.title=TRUE)
plot(iris.som, what="obs", type="barplot", print.title=TRUE)
plot(iris.som, what="obs", type="radar", key.loc=c(-0.5,5), mar=c(0,10,2,0))
```

Empty squares catch the attention: none of the  observations are assigned to
clusters 4, 9, 14 and 19.

Let us analyse the results for cluster 5 more precisely. On the ```"color"```
plots, cluster 5 has the following results: high values for
```Sepal.Width``` and low values for all the other variables.
These results are also noticeable in the other plots:

* on the ```"lines"``` plot, the second point shows a peak;
* on the ```"barplot"``` plot, this evolution is observed again;
* on the ```"radar"``` plot, the slice corresponding to ```Sepal.Width``` is 
very big whereas the two Petal slices are almost invisible.

Flowers with a long sepal and all the other variables being small are classified
in the top left corner of the map. Flowers with large sepals and petals (length
and width) are classified in the bottom right corner of the map whereas flowers
with small sepal length are classified in the top right corner of the map.

#### More graphics on observations

```{r irisObs}
plot(iris.som, what="obs", type="boxplot", print.title=TRUE)
rownames(iris)
plot(iris.som, what="obs", type="names", print.title=TRUE, scale=c(0.9,0.5))
```

* By default, the ```"boxplot"``` plot considers all numerical variables of the 
data set. As this example used the four numerical variables of the iris data 
set, four boxplots are plotted in each neuron.
* The ```"names"``` plot prints, for each neuron, the names of the observations 
which are assigned to its cluster. Here, the row names of iris data set are 
simply the row numbers.

#### More graphics on prototypes

Some more graphics handling prototypes have been implemented.

```{r irisProto}
par(mfrow=c(2,2))
plot(iris.som, what="prototypes", type="3d", variable=1, main="Sepal length")
plot(iris.som, what="prototypes", type="3d", variable=2, main="Sepal width")
plot(iris.som, what="prototypes", type="3d", variable=3, main="Petal length")
plot(iris.som, what="prototypes", type="3d", variable=4, main="Petal width")
```

```"3d"``` provides the same results as "color" but in a 3 dimensional graphic:
x is the x dimension of the grid, y is the y dimension of the grid and z is the
value of the prototype for the variable ```variable``` (by name or number in the
dataset) of the corresponding neuron.

Also, some graphics are provided to visualize the distance between prototypes on
the grid:
```{r irisDistProto}
plot(iris.som, what="prototypes", type="poly.dist")
plot(iris.som, what="prototypes", type="umatrix")
plot(iris.som, what="prototypes", type="smooth.dist")
plot(iris.som, what="prototypes", type="mds")
plot(iris.som, what="prototypes", type="grid.dist")
```

* ```"poly.dist"``` represents the distances between neighboring prototypes with
polygons plotted for each cell of the grid. The smaller the distance between 
a polygon's vertex and a cell border, the closer the pair of prototypes.
The colors indicates the number of observations in the neuron (white is used 
for empty neurons);

* ```"umatrix"``` fills the neurons of the grid using colors that represent
the average distance between the current prototype and its neighbors;

* ```"smooth.dist"``` plots the mean distance between the current prototype and 
its neighbors with a color gradation;

* ```mds``` plots the number of the neuron on a map according to a Multi
Dimensional Scaling (MDS) projection;

* ```grid.dist``` plots a point for each pair of prototypes, with x coordinates
representing the distance between the prototypes in the input space, and y
coordinates representing the distance between the corresponding neurons on the
grid.

These graphics show that there is a big gap (large distances) in the diagonal
that goes from the bottom left part of the map to its top right corner. Hence,
for instance, the flowers classified in cluster 15 are rather different from
those classified in cluster 20, even if these two clusters are adjacent on the
map. This is also emphasized on the MDS plot where clusters 15 and 20 are not
very close.

#### Graphics showing an additional variable

##### additional factor

Let us plot pies for all neurons according to the flower species of the 
observations.
```{r irisAdd1}
class(iris$Species)
levels(iris$Species)
plot(iris.som, what="add", type="pie", variable=iris$Species)
```

This figure shows that the clustering produced by the SOM is indeed relevant
to identify the three different species of iris: they are well separated on the
map and almost all clusters only contain one species of iris.

##### additional numerical vector

The ```"color"``` plot available for ```"add"``` is similar to the ```"obs"``` 
or ```"prototypes"``` cases. Here we choose as an additional variable the first 
variable of the iris data set, so we obtain the same graphic as above (see 
section __Graphics common to observations and prototypes__).
```{r irisAdd2}
plot(iris.som, what="add", type="color", variable=iris$Sepal.Length)
```

##### additional numerical matrix or data frame

The ```"lines"```, ```"barplot"```, ```"radar"``` and ```"boxplot"``` plots 
available for ```"add"``` are similar to the ```"obs"``` or ```"prototypes"``` 
cases.

```"words"``` is only implemented for an additional variable. In this case, the 
additional variable must be a contingency matrix: the words used on the graph 
are the names of the colums and the presence or lack of the word is expressed 
by respectively 1 or 0. The size of the words on the grid depends on the rate 
of presence in the observations of the current neuron.
```{r irisMatCont, echo=FALSE}
my.cont.mat <- matrix(data=c(rep(c(rep(1,50), rep(0,150)), 2), rep(1,50)), 
                      nrow=150, ncol=3)
colnames(my.cont.mat) <- levels(iris$Species)
```
```{r irisAdd4}
# my.cont.mat is the contingency matrix corresponding to the variable 
# iris$Species - overview of the 5 first lines:
my.cont.mat[1:5,]
plot(iris.som, what="add", type="words", variable=my.cont.mat)
```

##### additional non-numerical vector

```"names"``` is similar to the ```"names"``` case implemented for ```"obs"```. 
Here we choose to give the argument ```variable``` the row names of the iris 
data set: so we obtain the same graphic as above (see 
__More graphics on observations__).
```{r irisAdd5}
plot(iris.som, what="add", type="names", variable=rownames(iris),
     scale=c(0.9,0.5))
```

Then, we can try to call this plot again but on the variable ```iris$Species```.
```{r irisAdd5bis}
plot(iris.som, what="add", type="names", variable=iris$Species)
```

We obtain exactly the same plot as we obtained above for ```"words"``` with the 
contingency matrix corresponding to the variable ```iris$Species```.

### Analyze the projection quality

```{r irisQuality}
quality(iris.som)
```
By default, the quality function calculates both quantization and topographic 
errors. It is also possible to specify which one you want using the 
argument ```quality.type```.

The topographic error value varies between 0 (good projection quality) and 1 
(poor projection quality). Here, the topographic quality of the mapping is 
equal to 0.05 which means that all observations have a second best unit which is
in the neighborhood of the best matching unit.

The quantization error is an unbounded positive number. The closer it is from 0,
the better the projection quality.


### Building super classes from the resulting SOM

In the SOM algorithm, the number of clusters is necessarily close to the number 
of neurons on the grid (not necessarily equal as some neurons may have no 
observations assigned to them). This - quite large - number may not suit the 
original data for a clustering purpose.

A usual way to address clustering with SOM is to perform a hierarchical
clustering on the prototypes. This clustering is directly available in the
package ```SOMbrero``` using the function ```superClass```. To do so, you can
first have a quick overview to decide on the number of super clusters which 
suits your data.
```{r irisSC}
plot(superClass(iris.som))
```

By default, the function plots both a dendrogram and the evolution of the
percentage of explained variance. Here, 3 super clusters seem to be the best
choice. The output of ```superClass``` is a ```somSC``` class object.
Basic functions have been defined for this class:
```{r irisSC3}
my.sc <- superClass(iris.som, k=3)
summary(my.sc)
plot(my.sc, plot.var=FALSE)
```

Like ```plot.somRes```, the function ```plot.somSC``` has an argument ```type```
which offers many different plots and can thus be combined with most of the
graphics produced by ```plot.somSC```:
```{r irisSCplot, fig.width=6, fig.height=4}
plot(my.sc, type="grid", plot.legend=TRUE)
```
```{r irisSCplot3d}
plot(my.sc, type="dendro3d")
```

Case ```"grid"``` fills the grid with colors according to the super clustering 
(and can provide a legend).
Case ```"dendro3d"``` plots a 3d dendrogram.

A couple of plots from ```plot.somRes``` are also available for the super 
clustering. Some identify the super clusters with colors:
```{r irisSCplot2, fig.width=6, fig.height=4}
plot(my.sc, type="hitmap", plot.legend=TRUE)
```

```{r irisSCplot2B}
plot(my.sc, type="lines", print.title=TRUE)
plot(my.sc, type="barplot", print.title=TRUE)
plot(my.sc, type="boxplot", print.title=TRUE)
```
```{r irisSCplot2C, fig.height=4, fig.width=6}
plot(my.sc, type="mds", plot.legend=TRUE, cex =2)
```

And some others identify the super clusters with titles:
```{r irisSCplot3}
plot(my.sc, type="color")
plot(my.sc, type="poly.dist")
plot(my.sc, type="pie", variable=iris$Species)
plot(my.sc, type="radar", key.loc=c(-0.5,5), mar=c(0,10,2,0))
```

It is also possible to consider an additional variable using the 
argument ```add.variable```:
```{r irisSCplot4}
plot(my.sc, type="color", add.type=TRUE, variable=iris$Sepal.Length)
```

All these graphics help to see that super-cluster number 2 is found in the top
right corner of the map and contains only flowers of the "setosa" species. 
Super cluster number 3, at the bottom right corner of the map, contains mostly
flowers of the "virginica" species. The last super-cluster is located at the
top right corner of the map and mostly contains flowers of the species 
"versicolor".
