Example R programs and commands
			    28. Clustering 

# All lines preceded by the "#" character are comments or output.
# All other left-justified lines are my input


### Make sure that all graphics functions are available:
require(graphics)

### Generate a 2-cluster, 2-dimensional example,
#  with 50 bivariate normals per cluster in k=2 clusters.
xy <- rbind(matrix(rnorm(100, sd = 0.3), ncol = 2),
	  matrix(rnorm(100, mean = 1, sd = 0.3), ncol = 2))
colnames(xy) <- c("x", "y")
plot(xy)  # initial scatterplot, no clustering.
plot(xy, col=rep(1:2,each=50))   # Identify true clusters by color


########## Voronoi cell clustering around means, or barycenters
#
## Voronoi algorithm to update the means, assuming 2 clusters
# This chooses 2 initial centers at random from among the data:
k<-2;
cl <- kmeans(xy, centers=k)

plot(xy, col = cl$cluster)   # Identify clusters by color 1:k
points(cl$centers, col = 1:k, pch = 8) # Mark centers with *


################## Voronoi cell clustering around medoids 
#
install.packages("cluster")
library(cluster)    # for "pam()" and cluster graphing functions.

# Reuse the 2-d example from above

## Voronoi algorithm to update the medoids, assuming 2 clusters
# This chooses 2 initial centers at random from among the data:
k<-2;
pc <- pam(xy, k)

plot(xy, col = pc$cluster)   # Identify clusters by color 1:k
points(pc$medoids, col = 1:k, pch = 9, cex=2) # Mark centers with #

# Alternative plotting in principal component variables:
plot(pc)


## Compare means with medoids:
plot(xy)   # Scatterplot of clusters
points(cl$centers, pch = 8, cex=2) # Mark centers with big *
points(pc$medoids, pch = 9, cex=2) # Mark centers with big #


########### Agglomerative clustering
#
# Estimate the number of cluster from a dendrogram, or classification
# tree. 

library(cluster)  # for "agnes()" and "diana()" functions

### Agglomerative hierarchical clustering (agnes)
ah <- agnes(xy)  # notice that the number of clusters is not specified
plot(ah,which=2)  # plot only the dendrogram
#
# Commentary: the tree has an initial fork into two trunks with
# approximately equal numbers of branches, suggesting k=2 clusters.

### Divisive hierarchical clustering (diana)
dh <- diana(xy)  # notice that the number of clusters is not specified
plot(dh,which=2)  # plot only the dendrogram
#
# Commentary: again, the tree has an initial fork into two trunks with
# approximately equal numbers of branches, suggesting k=2 clusters.

# Add a third cluster and repeat:
xy3 <- rbind(xy, matrix(rnorm(100, mean=2,sd=0.3), ncol=2))
plot(xy3)

plot(agnes(xy3),which=2)   # 1+2 main trunks suggests k=3


plot(diana(xy3),which=2)   # 1+2 main trunks suggests k=3



############  Fisher Iris dataset -- 3 clusters in 4 dimensions

data(iris)
imeas <- iris[,1:4]   # just the 4 measurements

# Get a prior for the number of clusters
plot(agnes(imeas), which=2)   # unclear
plot(diana(imeas), which=2)   # more suggestive of k=3 clusters
#
# ... so use k=3 in the clustering functions
#
kmeans(imeas,3)    # 3 clusters of sizes 96, 33, 21 --- bad
summary(pam(imeas,3))  # 3 clusters of sizes 50, 62, 38 ---- better

# Misclassification rate:
truth <- rep(c(2,3,1), each=50);  # kmeans() names the species 2,3,1
guess <- kmeans(imeas, 3)$cluster
sum(truth!=guess)   ### OUTPUT: number of misclassifications, out of 150

truth <- rep(1:3, each=50);   # pam() names the species 1,2,3
guess <- pam(imeas, 3)$clustering
sum(truth!=guess)   ### OUTPUT: number of misclassifications, out of 150