<table width = 900 border = 0><tr><td>
# <b>NMST539 | Lab Session 10</b>

## <b>Cluster Analysis and Discriminant Analysis</b>
###  (application in R)

### <font color = "#880000"><b>LS 2017 | Tuesday 09/05/2017</b></font>
###### <a href = "cvicenie2017_10.Rmd">Rmd file (UTF8 encoding)</a>

**************
The R-software is available for download from the website: <a href = "https://www.r-project.org">https://www.r-project.org</a>

A user-friendly interface (one of many): <a href = "https://www.rstudio.com">RStudio</a>.<br><br>

<b>Manuals and introduction into <i>R</i>
(in Czech or English):</b><br>
    
<ul>
<li>Bína, V., Komárek, A. a Komárková, L.: <i>Jak na jazyk R.</i>  (<a href = "NMFM301/Rmanual2.pdf">PDF súbor</a>)</li>
<li>Komárek, A.: <i>Základy práce s R.</i>  (<a href = "http://msekce.karlin.mff.cuni.cz/~komarek/vyuka/dataRko/Rmanual1.pdf">PDF súbor</a>)</li>
<li>Kulich, M.: <i>Velmi stručný úvod do R.</i> (<a href = "NMFM301/uvodr.pdf" >PDF súbor</a>)</li>
<li>De Vries, A. a Meys, J.: <i>R for Dummies.</i> (ISBN-13: 978-1119055808) </li>  	
</ul>
**************


### <font color = "#880000"><b>1. Cluster Analysis in R</b></font>

The cluster analysis is a statistical method designed for distinguishing various objects, respectively grouping them into some disjoint sets with respect to their similarity/dissimilarity. The similarity/dissimilarty is usually measured with respect to some proximity or distance measure (see function `dist()` in the R environment - similarly, as it was used in case of the mutidimensional scaling technique). However, the cluster analysis is not just one specific method/approach but it is rather a whole set of various tools and algoritms, which can be used to solve the grouping problem. <br><br>

Considering the constructing proces we distinguish between two types of the clustering algorithms: partitioning algorithms where the assignments of objects into given 
groups can change during the algorithm and hierarchical algoritihms, where the assignments of objecs is kept fixed during the algorithm. The result of the clustering algorithm 
is mainly affected by the choice of the distance/proximity matrix and the clustering method used in the algorithm. 
Among others, we can include into the clustering approaches the following:
<ul>
<li>Partitioning clustering</li>
<li>Hierarchical clustering</li>
<li>Centroid-based clustering</li>
<li>Distribution-based clustering</li>
<li>Density-based clustering</li>
<li>and many others... </li>
</ul>

Let us consider the data on the consumption of automobiles in the United states in 80's (the dataset 'mtcars' available in R): 

```{r, warning = F, message = F}
rm(list = ls())
data <- mtcars
attach(data)
head(data)
```
<br><br>

The data consist of 32 observations (32 different cars on the US market available that time) and 11 covariates (some of them are continious, some of them are categorical). 
For an additional information on this dataset use the help session in R (type `?mtcars`). 

<br><br>

#### <font color = "#880000">Hierarchical Clustering in R</font>

Let us consider all covariates from the dataset. We will start by calculating the distance mantrix. In the following we 
are using the default option -  the Euclidian distance.<br><br>

In order to improve the visibility in the following dendograms (which we will construct later)
we will not consider the whole dataset ($n = 803$ students) but only a smaller subsample. A random subsample is created using the following piece of the R code (remmember to set the seed - command `set.seed()` in order to get the same data set each time you run the command). 

```{r}
D <- dist(mtcars) ### for the euclidian distance
```

Now, we can apply function `hclust()` (available under the standard R instalation) to run a hierarchical clustering approach based on the proximity matrix $D$. See the help session in R for more details (type `?hclust()`).

```{r, warning=F, message = F, fig.width=10, fig.height=8}
HC1 <- hclust(D)
plot(HC1, xlab = "Observations", ylab = "Proximity measure")
```

A little more fancy version available in the R software is the package called 'sparcl', which needs to be installed first (use `install.packages('sparcl')` for the instalation). After loading the library the corrresponding command is `ColorDendrogram()` (see the help session for more details). <br><br>

In order to specify groups (prespecified number) with respect to the proximity measure we can use 
the following R code: 

```{r, warning=F, message = F, fig.width=10, fig.height=8}
plot(HC1, xlab = "Observations", ylab = "Proximity measure")
groups <- cutree(HC1, k=3)
rect.hclust(HC1, k=3, border="red")
```

Different aglomerative approaches are possible with different settings of the 'method' parameter. 
The following options are possible: 

<ul>
<li>'method = 'single'</li>
<li>'method = 'complete'</li>
<li>'method = 'average'</li>
<li>'method = 'median'</li>
<li>'method = 'centroid'</li>
<li>'method = 'mcquitty'</li>
<li>'method = 'ward.D'</li>
</ul>

The default setting for the `hclust()` function is 'method = complete' which is also displayed on 
each dendogram figure. One can of course improve the overall impression by considering more fancy versions of dendograms. 

```{r, warning=F, message = F, fig.width=10, fig.height=8}
plot(HC1, xlab = "Observations", ylab = "Proximity measure", hang = -1)
groups <- cutree(HC1, k=3)
rect.hclust(HC1, k=3, border="red")
```

and even more sophisticated alternatives, e.g.: 

```{r, warning=F, message = F, fig.width=10, fig.height=8}
cutHC1 = cutree(HC1, 3)
labelColors = c("red", "green", "blue")

# function to get color labels
colLab <- function(n) {
    if (is.leaf(n)) {
        a <- attributes(n)
        labCol <- labelColors[cutHC1[which(names(cutHC1) == a$label)]]
        attr(n, "nodePar") <- c(a$nodePar, lab.col = labCol)
    }
    n
}
# using dendrapply
clusDendro = dendrapply(as.dendrogram(HC1), colLab)
plot(clusDendro, main = "Cool Dendrogram", type = "triangle")
```


The same idea of the code can be also used to define colors of branches with respect to some additional covariate (for instance, the covariate `gear`). 

```{r, warning=F, message = F, fig.width=10, fig.height=8}
local({
  colLab <<- function(n) {
    if(is.leaf(n)) {
      a <- attributes(n)
      i <<- i+1
      attr(n, "edgePar") <-
        c(a$nodePar, list(col = mycols[i], lab.font= i%%3))
    }
    n
  }
  mycols <- mtcars$gear
  i <- 0
})

x.clust.dend1 <- dendrapply(as.dendrogram(HC1), colLab)
plot(x.clust.dend1)
```



<br><br><br>

#### <font color = "#880000"><b>To Do</b></font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
 <ul>
 <li>Consider different types of distances for function `dist()` and apply the aglomerative 
 clustering approch while selecting different methods for the cluster calculations (various settings 
 of the 'method = ' parameter when calling `hclust()` function. </li>
 
 <li>Use some dendogram pictures and compare results from different distances and different clustering approaches.</li>
 
  <li>How many clusters seem to be a reasonable choice and what would be the natural interpretation of them?</li>
 </ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>


Now we apply the same approach for much bigger dataset (for this purpose we consider the data from the geographical questionare in the Czech republic and Slovakia) distinguishing for three clusters and using colors to distinguish for different grades (grade1, grade2 and grade3). The idea is to
see whether clusters somehow corresponds with three available grades (as we had maybe an impression above) or rather not... 

```{r, warning=F, message = F, fig.width=10, fig.height=8}
dataGeo <- read.csv("http://msekce.karlin.mff.cuni.cz/~maciak/NMST539/geoData2017.csv", header = T, sep = " ", dec = ".")
n <- dim(dataGeo)[1]
gColor <- rep("lightblue", n)
gColor[dataGeo[,2] == "grade2"] <- "blue"
gColor[dataGeo[,2] == "grade3"] <- "darkred"

local({
  colLab <<- function(n) {
    if(is.leaf(n)) {
      a <- attributes(n)
      i <<- i+1
      attr(n, "edgePar") <-
        c(a$nodePar, list(col = mycols[i], lab.font= i%%3))
    }
    n
  }
  mycols <- gColor
  i <- 0
})

Dall <- dist(dataGeo[, 6:27])
HC2 <- hclust(Dall)

x.clust.dend2 <- dendrapply(as.dendrogram(HC2), colLab)
plot(x.clust.dend2)

groups2 <- cutree(HC2, k=3)
rect.hclust(HC2, k=3, border="red")
```

More sophisticated dendograms can be obtained with some additional libraries (e.g. library 'ape' which can be installed by `install.packages("ape")` or library 'sparcl' which can be installed by running `install.packages("sparcl")`).

```{r, fig.width=10, fig.height=8}
library(sparcl)
ColorDendrogram(HC1, y = cutHC1, labels = names(cutHC1), main = "Mtcars dataset", 
    branchlength = 80)
```

<br><br>

#### <font color = "#880000">Non-hierarchical Clustering in R</font>

The most common approache for non-hierarchical clustering is so-called <b>K means algorithm</b>. The idea behind the algorithm is  to partition the observations into some prespecified number of groups ($k$ groups - thus the name of the method) such that the sum of squares from points to the assigned cluster centres is minimized. 

In the R software the function which performs the non-hierarchical K-means clustering is `kmeans()` available under the standard R instalation (see the help session for more details).<br><br>

```{r, fig.width=10}
Kmeans <- kmeans(mtcars, centers = 3)
colVector <- as.numeric(Kmeans$cluster)

plot(mtcars$mpg ~ mtcars$hp, bg = colVector, xlab = "Horse Power", ylab = "Miles per Gallon", pch = 21, col = "black")
points(Kmeans$centers[,1] ~ Kmeans$centers[,4], col = 1:3, pch = 8, cex = 2)
```

#### <font color = "#880000">Comment</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
 <ul>
 <li>Some other methodologies and approaches are avialable in  different R packages and various
 R functions. </li>
 </ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

### <font color = "#880000"><b>2. Introduction to discriminant analysis in R</b></font>

Discriminant analysis is a (not purely a multivariate) statistical method which was designed for allocating various objects into similar (homogenous) groups. The idea is to use some available information (given covariates) in order to create some classification rule which can be used later to assign some new object into one of the existing groups. 

There are many different approaches how to define a descrimination rule. The most classical approaches include a linear discriminant function (in the R software implemented in `lda()` function available in the packages 'MASS') and quadratic discriminant function (command 'qda()' in package 'MASS'). 

Use the help sessions in R for both functions to find out more about there implementation in R. 

The following is a nice R Markdown available in the R studio prepared by G. Martos:

<a href = "http://rstudio-pubs-static.s3.amazonaws.com/35817_2552e05f1d4e4db8ba87b334101a43da.html">http://rstudio-pubs-static.s3.amazonaws.com/35817_2552e05f1d4e4db8ba87b334101a43da.html</a>


For instance, we can use some continuous covariates to try to classify the avialable cars by the number of cylinders.

```{r, fig.width = 10}
library(MASS)
lda1 <- lda(cyl ~ mpg + hp + qsec + drat + disp, data = mtcars)

pred <- predict(lda1)
ldahist(data = pred$x[,1], g=cyl)
```


<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>

#### <font color = "#880000">To Do</font>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
 <ul>
 <li>Use the R Markdown available on the given website and apply some linear and quadratic 
 distriminant functions to some dataset.</li>
 
 <li>Provide some graphical outputs and discuss obtained results.</li>
 
 <li>Use some additional datapoints (feel free to make up some) and try to use the discrimination 
 rules to assign the new points to the existing groups. </li>
 </ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>




<br><br>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<hr size="3" width="100%" noshade style="color:#440000" align="left" /><br><br> 




## <font color = "#880000"><b>Homework Assignment</b></font>
### (Deadline: sixth lab session / 11.04.2007)

Use the command below to instal the `SMSdata` packages with various multivariate datasets: 

```{r, eval = F}
install.packages(pkgs="http://www.karlin.mff.cuni.cz/~hlavka/sms2/SMSdata_1.0.tar.gz", repos=NULL, type="source")
```

Chose one dataset of your choise from the list of all available datasets in the package:

```{r, eval = F}
data(package = "SMSdata") 
```
 
 There are 21 different datasets and you can load each of them by typing and running 
`data(dataset_name)`, for instance:

```{r, eval = F}
data(plasma) 
```

<ul>
<li>Choose one dataset and apply the cluster analysis; </li>

<li>Discuss the optimal choice of the number of the clusters and interpret the results;</li>
</ul>




</td></tr></table>
<br><br><br>