### R code from vignette source 'chap11-clustering.Rnw'

###################################################
### code chunk number 1: chap11-clustering.Rnw:5-6
###################################################
options(show.signif.stars=FALSE,digits=4,width=90)


###################################################
### code chunk number 2: chap11-clustering.Rnw:34-38
###################################################
library(alr3)
head(BigMac2003)
library(psych)
describe(BigMac2003)[ , c(2:5, 9:10)]


###################################################
### code chunk number 3: chap11-clustering.Rnw:41-43
###################################################
BigMac2003 <- BigMac2003[,-9]
as.matrix(dist(scale(BigMac2003)))[1:7,1:7]


###################################################
### code chunk number 4: chap11-clustering.Rnw:46-47
###################################################
as.matrix(dist(BigMac2003))[1:7,1:7]


###################################################
### code chunk number 5: one
###################################################
hc <- hclust(dist(scale(BigMac2003)), "single")
plot(hc, cex=.55, xlab="2003 Big Mac Data, single-link clustering")


###################################################
### code chunk number 6: two
###################################################
hc2 <- hclust(dist(scale(BigMac2003)), "complete")
plot(hc2, cex=.55, xlab="2003 Big Mac Data, complete-link clustering")


###################################################
### code chunk number 7: twoa
###################################################
pc1 <- princomp(BigMac2003,cor=TRUE)
plot(pc1$scores[,1],pc1$scores[,2],type="n")
text(pc1$scores[,1],pc1$scores[,2],
      substr(rownames(BigMac2003),1,5),col=cutree(hc2,5))
legend("topleft",paste("Clust",1:5), col=1:5, pch="x", cex=0.9,inset=0.02)


###################################################
### code chunk number 8: twob
###################################################
summary(pc1 <- princomp(BigMac2003, cor=TRUE))
hc3 <- hclust(dist(pc1$scores[, 1:4]), "complete")
plot(hc3, cex=.6, xlab="2003 Big Mac Data, Mahalanobis Distance")


###################################################
### code chunk number 9: chap11-clustering.Rnw:88-89
###################################################
table(cutree(hc2, 5), cutree(hc3, 5))


###################################################
### code chunk number 10: three
###################################################
summary(BigPow <- powerTransform(BigMac2003, family="yjPower"))
BigTran <- yjPower(BigPow$y, coef(BigPow))
(pc2 <- princomp(BigTran, cor=TRUE))
hc4 <- hclust(dist(scale(pc2$scores[, 1:4])), "complete")
plot(hc4, cex=.6, xlab="Normalized Big Mac data")


###################################################
### code chunk number 11: threea
###################################################
plot(pc2$scores[,1], pc2$scores[,2], type="n")
text(pc2$scores[,1], pc2$scores[,2],
      substr(rownames(BigMac2003), 1, 5), col=cutree(hc4, 5))
legend("bottomleft", paste("Clust", 1:5), col=1:5, pch="x", cex=0.9)
pc2$loadings[, 1]


###################################################
### code chunk number 12: chap11-clustering.Rnw:115-116
###################################################
table(cutree(hc2, 5),cutree(hc4, 5))


###################################################
### code chunk number 13: chap11-clustering.Rnw:120-127
###################################################
data <- read.csv("Nowakskulls.csv", header=TRUE)
# remove duplicate rows 45:55, which are included as both
# Minnesota and NE Minnesota
data <- data[-(56:66), ]
data$Source <- factor(data$Source)
data$Source <- factor(data$Source,
  labels=c("1 Alg", "2 MN70s", "3 NEMNold", "4 WUS"))


###################################################
### code chunk number 14: chap11-clustering.Rnw:129-131
###################################################
describe(data)
xtabs( ~ Source, data)


###################################################
### code chunk number 15: wolf1
###################################################
p1 <- prcomp(data[, -1], scale=TRUE)
summary(p1)
plot(p1$x[, 1:2], col=as.numeric(data$Source),
   pch=as.character(as.numeric(data$Source)),
   main="Two PC solution")


###################################################
### code chunk number 16: wolf2
###################################################
plot(h2<-hclust(dist(scale(data[,-1], center=FALSE)), method="average"),
  labels=substr(data$Source, 1, 1),
  frame.plot=TRUE,
  main="Nowak Wolf Skulls, clustering of individuals",
  xlab="Average Link Clustering", sub="")
rect.hclust(h2, k=2, border="red")
# confusion matrix
cluster.number <- cutree(h2, k=2)
xtabs(~cluster.number + Source, data)


###################################################
### code chunk number 17: chap11-clustering.Rnw:168-172
###################################################
m1 <- lm(as.matrix(data[, -1]) ~ Source -1, data)
raw.means <- coef(m1)
rownames(raw.means) <- substr(rownames(raw.means),7,20)
t(raw.means)


###################################################
### code chunk number 18: wolf3
###################################################
plot(h1<-hclust(dist(scale(raw.means, center=FALSE)), method="average"),
  frame.plot=TRUE, main="",
  xlab="Average Link Clustering", sub="")


###################################################
### code chunk number 19: four
###################################################
dd <- as.dist(1 - cor(USJudgeRatings[,-1]))
round(1000 * dd) # (prints more nicely)
plot(hclust(dd),xlab="Complete linkage")


###################################################
### code chunk number 20: chap11-clustering.Rnw:228-231
###################################################
X <- pc2$scores[,1:4]
(initial <- tapply(X, list(rep(cutree(hc4, 5), ncol(X)), col(X)), mean))
km <- kmeans(X, initial)


###################################################
### code chunk number 21: chap11-clustering.Rnw:234-248
###################################################
pr <- function (x, ...)
{
    cat("K-means clustering with ", length(x$size), " clusters of sizes ",
        paste(x$size, collapse = ", "), "\n", sep = "")
    cat("\nCluster means:\n")
    print(x$centers, ...)
    cat("\nWithin cluster sum of squares by cluster:\n")
    print(x$withinss, ...)
    cat("\nAvailable components:\n")
    print(names(x))
    invisible(x)
}
pr(km)
table(km$cluster,cutree(hc4,5))


###################################################
### code chunk number 22: km
###################################################
pr(km1 <- kmeans(X, 5, nstart=25))
table(km1$cluster, km$cluster)
par(mfrow=c(1,2))
plot(X[, 1], X[, 2], type="n", main="Start=complete linkage")
text(X[, 1], X[, 2], cex=.8,
      substr(rownames(BigMac2003), 1, 5), col=km$cluster)
plot(X[, 1], X[, 2],type="n", main="25 random starts")
text(X[, 1], X[, 2], cex=.8,
      substr(rownames(BigMac2003), 1, 5), col=km1$cluster)