<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/hierarchical-clustering-nearest-neighbors-algorithm-in-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.aaronschlegel.com/hierarchical-clustering-nearest-neighbors-algorithm-r/"> R – Aaron Schlegel</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>Hierarchical clustering is a widely used and popular tool in statistics and data mining for grouping data into ‘clusters’ that exposes similarities or dissimilarities in the data. There are many approaches to hierarchical clustering as it is not possible to investigate all clustering possibilities. One set of approaches to hierarchical clustering is known as agglomerative, whereby in each step of the clustering process an observation or cluster is merged into another cluster. The first approach we will explore is known as the <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Single-linkage_clustering">single linkage</a> method, also known as nearest neighbors.</p>
<p>The data we will cluster is seven different measures of air pollution of 41 cities throughout the United States (qtd. Rencher 502). The data were obtained from the <a target="_blank" rel="nofollow">companion FTP site</a> of the book Methods of Multivariate Analysis and contains the variables as follows:</p>
<ul>
<li>SO2 content of air in mg/cm</li>
<li>Average annual temperature in F</li>
<li>Number of manufacturing plants with 20 or more workers</li>
<li>Population size as recorded by 1970 census in thousands</li>
<li>Average annual wind speed in mph</li>
<li>Average annual precipitation in inches</li>
<li>Average number of days with precipitation per year</li>
</ul>
<pre>
poll &lt;- read.table('T15_13_POLLUTION.dat', col.names = c('city', 'so2', 'tempf', 'manufacturing',
                                                         'pop', 'wind', 'precip', 'days.w.precip'))
</pre>
<h6>Variable Scaling Consideration</h6>
<p>Although recommended by many authors, the question of scaling in the context of hierarchical clustering (particularly using the Euclidean distance measure) is not so black and white (Rencher 2002, pp. 454). The variables that best separate clusters might not do so after scaling. We will scale the pollution data because the variables’ scale of measurement is quite different from each other, but it should be noted standardization should be thoughtfully applied (or not applied) rather than proceeding to scale automatically.</p>
<pre>
poll.scale &lt;- data.frame(scale(poll[,2:8]))
</pre>
<h6>Measures of Similarity (or Dissimilarity) for Cluster Analysis</h6>
<p>Before performing hierarchical clustering, we must find the similarity between each pair of observations, which is referred to as the distance. The distance measure is more of a measure of dissimilarity as it increases as the observations move farther away. As in agglomerative hierarchical clustering, there are many approaches to measuring the distance, the most common of which is the <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Euclidean_distance">Euclidean distance</a>. There is no ‘best’ distance measure, and one must consider the data at hand and the assumptions of each measure to select an appropriate method. Euclidean distance is the straight line between two pairs of observations and is defined as:</p>
<p><span data-display="true" class="wp-katex-eq"> d(x,y) = \sqrt{(x – y)^\prime (x – y)} = \sqrt{\sum^p_{j=1} (x_j – y_j)^2} </span></p>
<p>The following function implements the Euclidean distance calculations for each pair of observations in the dataset.</p>
<pre>
euclidean.distance &lt;- function(x) {

  n &lt;- nrow(x)
  dist.mat &lt;- matrix(0, n, n)
  xj &lt;- x[1,]
  for (i in 1:n) {
    for (j in 1:n) {
      yj &lt;- x[j,]
      d &lt;- sqrt(as.matrix(xj - yj) %*% as.matrix(t(xj - yj)))
      dist.mat[j,i] &lt;- d
    }
    xj &lt;- x[1+i,]
  }
  return(dist.mat)
}
</pre>
<pre>
euclid.dist &lt;- euclidean.distance(poll.scale)
</pre>
<p>The <code>dist()</code> function in R also calculates the distance.</p>
<pre>
poll.dist &lt;- dist(poll.scale, method = 'euclidean')
</pre>
<p>The output of the <code>dist()</code> function is an atomic vector. We convert it to a matrix here to compare our results with the function.</p>
<pre>
poll.dist.mat &lt;- as.matrix(poll.dist)
</pre>
<p>As the resulting matrices are <span data-display="false" class="wp-katex-eq">41 \times 41</span>, check the first and bottom three rows of each to verify.</p>
<pre>
cbind(head(poll.dist.mat[,1:3]), tail(poll.dist.mat[,1:3]))
##          1        2        3        1        2        3
## 1 0.000000 4.789018 3.171606 4.034076 3.264390 1.782405
## 2 4.789018 0.000000 3.009865 5.751432 2.007170 3.321471
## 3 3.171606 3.009865 0.000000 4.790368 1.171199 2.906303
## 4 3.871066 3.491389 1.262450 6.637764 3.208636 4.153093
## 5 6.230609 2.817075 3.800426 5.675892 2.526646 4.136598
## 6 5.305038 1.729939 2.964289 6.546541 4.008765 3.474188
</pre>
<pre>
cbind(head(euclid.dist[,1:3]), tail(euclid.dist[,1:3]))
##           [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
## [36,] 0.000000 4.789018 3.171606 4.034076 3.264390 1.782405
## [37,] 4.789018 0.000000 3.009865 5.751432 2.007170 3.321471
## [38,] 3.171606 3.009865 0.000000 4.790368 1.171199 2.906303
## [39,] 3.871066 3.491389 1.262450 6.637764 3.208636 4.153093
## [40,] 6.230609 2.817075 3.800426 5.675892 2.526646 4.136598
## [41,] 5.305038 1.729939 2.964289 6.546541 4.008765 3.474188
</pre>
<h6>Hierarchical Clustering with Single Linkage</h6>
<p><a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Johnson%27s_algorithm">Johnson’s algorithm</a> describes the general process of hierarchical clustering given <span data-display="false" class="wp-katex-eq">N</span> observations to be clustered and an <span data-display="false" class="wp-katex-eq">N \times N</span> distance matrix. The steps of Johnson’s algorithm as applied to hierarchical clustering is as follows:</p>
<ol>
<li>Begin with disjoint clustering with level <span data-display="false" class="wp-katex-eq">L(0) = 0</span> and <span data-display="false" class="wp-katex-eq">m = 0</span>.</li>
<li>In the case of single linkage, find the pair with the minimum distance, with pairs denoted as <span data-display="false" class="wp-katex-eq">r</span> and <span data-display="false" class="wp-katex-eq">s</span>, according to:</li>
</ol>
<ul>
<li><span data-display="false" class="wp-katex-eq">d[(r), (s)] = min (d[(i),(j)])</span></li>
</ul>
<ol start="3">
<li>Add one to <span data-display="false" class="wp-katex-eq">m</span>, <span data-display="false" class="wp-katex-eq">m = m + 1</span>. Merge clusters <span data-display="false" class="wp-katex-eq">r</span> and <span data-display="false" class="wp-katex-eq">s</span> into one cluster to form the next clustering at <span data-display="false" class="wp-katex-eq">m</span>. <span data-display="false" class="wp-katex-eq">L(m)</span> then becomes: <span data-display="false" class="wp-katex-eq">L(m) = d[(r), (s)]</span></li>
<li>The distance matrix is updated by removing the rows and columns corresponding to clusters <span data-display="false" class="wp-katex-eq">r</span> and <span data-display="false" class="wp-katex-eq">s</span> and inserting a row and column for the newly formed cluster. The distance between the newly formed cluster <span data-display="false" class="wp-katex-eq">(r,s)</span> and the old cluster <span data-display="false" class="wp-katex-eq">k</span> is calculated again with the minimum distance (in the case of single linkage):</li>
</ol>
<ul>
<li><span data-display="false" class="wp-katex-eq">d[(k), (r,s)] = min (d[(k),(r)], d[(k),(s)])</span></li>
</ul>
<ol start="5">
<li>Stop if all <span data-display="false" class="wp-katex-eq">N</span> observations are in one cluster, otherwise repeat starting at Step 2.</li>
</ol>
<h6>Implementing the Single Linkage Hierarchical Clustering Technique</h6>
<p>Although hierarchical clustering with a variety of different methods can be performed in R with the <code>hclust()</code> function, we can also replicate the routine to an extent to better understand how Johnson’s algorithm is applied to hierarchical clustering and how <code>hclust()</code> works. To plot our clustering results, we will rely on one output of the <code>hclust()</code> function, the <code>order</code> value, which is purely used for plotting. To plot and compare our clustering results, the best method (that I know of) is to create an <code>hclust</code> object, which according to <code>?hclust</code>, requires <code>merge</code>, <code>height</code> and the aforementioned <code>order</code> components.</p>
<p>Quoting from the <code>hclust</code> documentation of the required components:</p>
<blockquote><p>
merge is an n-1 by 2 matrix. Row i of merge describes the merging of clusters at step i of the clustering. If an element j in the row is negative, then observation -j was merged at this stage. If j is positive then the merge was with the cluster formed at the (earlier) stage j of the algorithm. Thus negative entries in merge indicate agglomerations of singletons, and positive entries indicate agglomerations of non-singletons.
</p>
</blockquote>
<blockquote><p>
height is a set of n-1 real values (non-decreasing for ultrametric trees). The clustering height: that is, the value of the criterion associated with the clustering method for the particular agglomeration.
</p>
</blockquote>
<blockquote><p>
Order is a vector giving the permutation of the original observations suitable for plotting, in the sense that a cluster plot using this ordering and matrix merge will not have crossings of the branches.
</p>
</blockquote>
<p>Therefore we need to build the <code>merge</code> and <code>height</code> objects. We will only use the <code>order</code> output of <code>hclust()</code> to plot.</p>
<p>To start, initialize the <code>merge</code> and <code>height</code> objects and set the diagonal of the distance matrix to <code>Inf</code> to avoid including 0s in our calculations.</p>
<pre>
Lm &lt;- data.frame(euclid.dist)
N &lt;- nrow(Lm)

merge &lt;- matrix(0, N-1, 2)
height &lt;- vector(length = N-1)
diag(Lm) &lt;- Inf
</pre>
<p>The row and column names of the distance matrix are set to <span data-display="false" class="wp-katex-eq">-1</span> increasing up to <span data-display="false" class="wp-katex-eq">-N</span>, which allows us to note if the merged cluster at a particular step was an individual cluster.</p>
<pre>
colnames(Lm) &lt;- -(1:N)
rownames(Lm) &lt;- -(1:N)
</pre>
<p>The following loop builds the <code>merge</code> and <code>height</code> objects needed for plotting the hierarchical clustering results. A <a target="_blank" rel="nofollow" href="https://home.deib.polimi.it/matteucc/Clustering/tutorial_html/hierarchical.html">visual example of how clustering is performed</a> shows how the rows and columns are merged with a simple dataset.</p>
<pre>
for (m in 1:(N-1)) { 
  cols &lt;- colnames(Lm)
  
  # Johnson's algorithm Step 2L Find the pair with the minimum distance
  
  # The which() function returns the row and column position of the pair
  d &lt;- which(Lm == min(Lm), arr.ind = TRUE)[1,,drop=FALSE]
  
  height[m] &lt;- min(Lm) # The height is the value of the pair with the minimum distance
  
  # The row and column position of the minimum pair is stored as sequence m in the merge object
  merge[m,] &lt;- as.numeric(cols[d])
  
  # Johnson's algorithm Step 3: The pair with the minimum distance is merged
  
  # The cluster object is used to find previous clusters that the pair belong to (if they exist)
  # Does this by finding any columns above 0 (since all column names are negative, a positive 
  # column value implies it has been clustered)
  cluster &lt;- c(d, which(cols %in% cols[d[1, cols[d] > 0]]))

  colnames(Lm)[cluster] &lt;- m # Rename the columns indicated by cluster to the sequence number, m

  # Merge the pairs according to Johnson's algorithm and the single linkage method
  sl &lt;- apply(Lm[d,], 2, min)
  
  # Johnson's algorithm Step 4: Remove column and row corresponding to old clusters and
  # insert a new column and row for newly formed cluster.
  
  # The insertion of the cluster is done by setting the first sequential row and column of the
  # minimum pair in the distance matrix (top to bottom, left to right) as the cluster resulting 
  # from the single linkage step
  Lm[min(d),] &lt;- sl
  Lm[,min(d)] &lt;- sl
  
  # Make sure the minimum distance pair is not used again by setting it to Inf
  Lm[min(d), min(d)] &lt;- Inf
  
  # The removal step is done by setting the second sequential row and column of the minimum pair
  # (farthest right, farthest down) to Inf
  Lm[max(d),] &lt;- Inf
  Lm[,max(d)] &lt;- Inf
}
</pre>
<h6>Plotting the Hierarchical Clustering as a Dendrogram</h6>
<p>We now have the <code>merge</code> and <code>height</code> components of an <code>hclust</code> object. The <code>order</code> component comes from the <code>hclust()</code> function. As noted, the <code>order</code> component is used just by <code>hclust</code> to plot and has no bearing on our understanding of cluster analysis.</p>
<pre>
poll.clust &lt;- hclust(poll.dist, method = 'single')
</pre>
<p>According to the documentation in <code>?hclust</code>, an hclust object is a list with class hclust of the above components. We construct the hclust class with the following:</p>
<pre>
hclust.obj &lt;- list() # Initialize an empty list
hclust.obj$merge &lt;- merge # Add the merge component obtained earlier
hclust.obj$order &lt;- poll.clust$order # Here the order component from hclust is added to the list
hclust.obj$height &lt;- height # The height component determines the lengths of the dendogram nodes
hclust.obj$labels &lt;- poll$city # Add the city names to the labels component
class(hclust.obj) &lt;- 'hclust' # The list is set to class hclust
</pre>
<pre>
plot(hclust.obj)
</pre>
<p><a target="_blank" rel="nofollow" href="http://www.aaronschlegel.com/?attachment_id=22498"><img sizes="(max-width: 1060px) 100vw, 1060px" alt="hierarchical cluster single linkage" class="wp-image-22498" src="https://i1.wp.com/www.aaronschlegel.com/wp-content/uploads/hier_clust_mine.png?resize=1060%2C757" data-recalc-dims="1"></a></p>
<p>Plot the hierarchical clustering as given by the <code>hclust()</code> function to verify our results.</p>
<pre>
plot(poll.clust, labels = poll$city)
</pre>
<p><a target="_blank" rel="nofollow" href="http://www.aaronschlegel.com/?attachment_id=22500"><img sizes="(max-width: 1060px) 100vw, 1060px" alt="hierarchical clustering single linkage from hclust" class="wp-image-22500" src="https://i0.wp.com/www.aaronschlegel.com/wp-content/uploads/hier_clust_hclust-2.png?resize=1060%2C757" data-recalc-dims="1"></a></p>
<h6>Interpretation of the Cluster Plot</h6>
<p>Moving left to the right along the x-axis, we see there Chicago, Phoenix, and Philadelphia are separated from other cities to a decreasing degree. Otherwise, there isn’t too much distinction between clusters. At the height of slightly more than two we have a three cluster solution (counting the number of lines crossed going left to right on the x-axis), but the third cluster would contain the bulk of the cities and thus isn’t particularly useful. Using a lower height increases the number of clusters significantly, for example, a height of about 1.9 gives a six-cluster solution. However, the final cluster still contains a large amount of the cities and therefore isn’t very descriptive of how the cities group themselves. This lack of separation in the clusters might indicate there isn’t too much variation in the cities’ pollution levels to begin with, the method of clustering used, or a combination of the two. The former point is further evidenced by small distances in the y-axis between many of the clusters.</p>
<p>We will see in later posts if using a different distance measure or clustering method improves the clustering of the cities. The hierarchical cluster could be ‘cut’ to the number of groups we are interested in, but since the clustering isn’t that great, we will save this step once for when a better clustering and distance method are used.</p>
<h6>Disadvantages of the Single Linkage Approach</h6>
<p>The single linkage (nearest neighbors) approach has several disadvantages compared to other clustering methods and is therefore not recommended by many authors (Rencher, 2002, pp. 475). Single linkage is rather prone to chaining (also known as space-contracting), which is the tendency for newly formed clusters to move closer to individual observations, so observations end up joining other clusters rather than another individual observation. We can see this markedly in the hierarchical cluster plot above. The single linkage method is also sensitive to errors in distances between observations.</p>
<p>It should be noted there is no clear ‘best’ clustering method and often a good approach is to try several different methods. If the resulting clusters are somewhat similar, that is evidence there may be natural clusters in the data. Many studies, however, conclude Ward’s method and the average linkage method are the overall best performers (Rencher, 2002, pp. 479).</p>
<h6>References</h6>
<p>Gan, Guojun, Chaoqun Ma, and Jianhong Wu, Data Clustering: Theory, Algorithms, and Applications, ASA-SIAM Series on Statistics and Applied Probability, SIAM, Philadelphia, ASA, Alexandria, VA, 2007.</p>
<p>Murtagh, F. (n.d.). Multivariate Data Analysis with Fortran C and Java. Queen’s University Belfast.</p>
<p>Rencher, A. (2002). Methods of Multivariate Analysis (2nd ed.). Brigham Young University: John Wiley &amp; Sons, Inc.</p>
<p><a target="_blank" rel="nofollow" href="http://www.cse.iitb.ac.in/infolab/Data/Courses/CS632/1999/clustering/node12.html#SECTION00033100000000000000">http://www.cse.iitb.ac.in/infolab/Data/Courses/CS632/1999/clustering/node12.html#SECTION00033100000000000000</a></p>
<p><a target="_blank" rel="nofollow" href="http://84.89.132.1/~michael/stanford/maeb7.pdf">http://84.89.132.1/~michael/stanford/maeb7.pdf</a></p>
<p><a target="_blank" rel="nofollow" href="http://www.stat.berkeley.edu/~spector/s133/Clus.html">http://www.stat.berkeley.edu/~spector/s133/Clus.html</a></p>
<p><a target="_blank" rel="nofollow" href="http://www.pbarrett.net/techpapers/euclid.pdf">http://www.pbarrett.net/techpapers/euclid.pdf</a></p>
<p>The post <a target="_blank" rel="nofollow" href="http://www.aaronschlegel.com/hierarchical-clustering-nearest-neighbors-algorithm-r/">Hierarchical Clustering Nearest Neighbors Algorithm in R</a> appeared first on <a target="_blank" rel="nofollow" href="http://www.aaronschlegel.com/">Aaron Schlegel</a>.</p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.aaronschlegel.com/hierarchical-clustering-nearest-neighbors-algorithm-r/"> R – Aaron Schlegel</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...
<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/wD0xJ90FyyI"></a></div>