###################################################
### R Code for the book chapter:
###
### Friedrich Leisch. Visualizing cluster analysis and finite mixture
### models. In Chunhouh Chen, Wolfgang Härdle, and Antony Unwin,
### editors, Handbook of Data Visualization, Springer Handbooks of
### Computational Statistics. Springer Verlag, 2008.
###
### (c) 2008 Friedrich Leisch
###
### Figures may appear different than in the book due to changes in
### the corresponding R packages.
###################################################


###################################################
### chunk number 1: 
###################################################

options(digits=2)
library("flexclust")
library("vcd")
data("dentitio")

colnames(dentitio) <- sub("..\\.", "\\.", colnames(dentitio))
dentitio <- as.matrix(dentitio)

ahclust <- function(x) hclust(x, "average")
mdist <- function(x) dist(x, "manhattan")

h1 <- ahclust(mdist(dentitio))
d1 <- as.dendrogram(h1)

qdent <- quantile(as.vector(dentitio))


###################################################
### chunk number 2: Table 11.1
###################################################
head(dentitio, 10)


###################################################
### chunk number 3: 
###################################################
x2 <- bundestag(2005)
colnames(x2)[2] <- "CDU/CSU"
p2 <- prcomp(x2)

library("cluster")

set.seed(2409)
pam2.5 <- pam(x2, 5)
k2.5 <- as.kcca(pam2.5)


###################################################
### chunk number 4: Table 11.2
###################################################
round(head(x2, 10),2)

###################################################
### chunk number 5: Figure 11.1
###################################################
plot(p2$x, col=grey(.7))
projAxes(p2, radius=0.22)


###################################################
### chunk number 6: Figure 11.2
###################################################
plot(d1, horiz=TRUE)


###################################################
### chunk number 7: 
###################################################
hitze <- diverge_hcl(5)


###################################################
### chunk number 8: Figure 11.3 (Color)
###################################################
heatmap(dentitio, hclustfun=ahclust, distfun=mdist, 
        scale="none", Colv=1:8, breaks=c(-1,qdent), col=hitze)

###################################################
### chunk number 9: Figure 11.3 (Grey)
###################################################
hitze <- grey(seq(0,1,length=5))
heatmap(dentitio, hclustfun=ahclust, distfun=mdist, 
        scale="none", Colv=1:8, breaks=c(-1,qdent), col=hitze)

###################################################
### chunk number 10: Figure 11.4
###################################################
plot(k2.5, data=x2, project=p2, hull="convex")


###################################################
### chunk number 11: Figure 11.5
###################################################
set.seed(1)
image(k2.5, xlab="SPD", ylab="CDU/CSU", clwd=1, npoints=200, graph=FALSE, col=grey(c(.6,.8,1)))
points(x2)
plot(k2.5, add=T)


###################################################
### chunk number 12: Figure 11.6
###################################################
plot(silhouette(pam2.5), main="")


###################################################
### chunk number 13: Figure 11.7
###################################################
print(barchart(k2.5))


###################################################
### chunk number 14: Figure 11.8
###################################################
library(lattice)
cl <- factor(predict(k2.5))
parallel(~x2[,5:1]|cl, col="black", lty=1, as.table=TRUE)


###################################################
### chunk number 16:
###################################################
library(fpc)
p3 <- discrproj(x2, predict(k2.5)!=2, method="anc")
class(p3) <- "discrproj"
predict.discrproj <- function(object, data, ...)
{
    data %*% object$units
}


###################################################
### chunk number 17: Figure 11.9
###################################################
plot(k2.5, data=x2, project=p3, hull="convex")


###################################################
### chunk number 18: Table 11.3
###################################################
require("gplots")

state <- bundestag(2005, state=TRUE)
links <- tapply(x2[,"LINKE"], list(state), median)
ORD <- order(links)
ORD <- names(links)[ORD]
TAB <- table(state, pam2.5$clustering)[ORD,c(5,4,1,2,3)]

DTAB <- as.data.frame(as.table(TAB))
balloonplot(list(Cluster=DTAB[[2]]), list(State=DTAB[[1]]), DTAB[[3]],
            main="", rowmar=4, label.lines=FALSE,
            xaxs="i", yaxs="i", dotcolor="grey")


###################################################
### chunk number 19: Figure 11.11
###################################################
mosaic(TAB, shade=T, labeling_args =list(rot_labels=0,
                     just_labels="right"), 
       margin=c(left=12, top=2, bottom=2), 
       keep=FALSE)


###################################################
### chunk number 20: 
###################################################
somgridlines <- function(x, project=NULL, which=1:2)
{
    if(!is.null(project))
        x$codes <- predict(project, x$codes)
    
    X <- Y <- matrix(0, x$grid$xdim, x$grid$ydim)
    X[x$grid$pts] <- x$codes[,which[1]]
    Y[x$grid$pts] <- x$codes[,which[2]]
    X <- t(X)
    Y <- t(Y)

    M <- nrow(X)
    N <- ncol(X)
    
    for(m in 1:M){
        for(n in 1:N){
            if(m<M) segments(X[m,n],Y[m,n],X[m+1,n],Y[m+1,n], col=2)
            if(n<N) segments(X[m,n],Y[m,n],X[m,n+1],Y[m,n+1], col=2)
        }
    }

    points(x$codes[,which], col="red", pch=21, cex=3, bg="white")
    text(x$codes[,which], labels=1:(N*M))
}

somcodes <- function (x, ...) 
{
    symbols(x$grid$pts[, 1], x$grid$pts[, 2], 
            circles = rep(0.5, nrow(x$grid$pts)), 
            inches = FALSE, bg = "white",
            axes=FALSE, xlab="", ylab="",xpd=NA)

    codes <- x$codes
    rownames(codes) <- 1:nrow(codes)

    stars(codes, location = x$grid$pts, len = 0.4, 
          add = TRUE, col.segments = grey(seq(1,0,length=5)), 
          draw.segments = TRUE, xpd=NA)
    invisible(codes)
}




###################################################
### chunk number 21: 
###################################################
require("kohonen")
set.seed(1)
som2 = som(x2, somgrid(4,4))


###################################################
### chunk number 22:  Figure 11.12
###################################################
plot(p2$x, col=grey(.7))
somgridlines(som2, project=p2)


###################################################
### chunk number 23: Figure 11.13
###################################################
par(mar=rep(1,4))
somcodes(som2)


###################################################
### chunk number 24: 
###################################################
pie(rep(1,5), colnames(x2), col=grey(seq(1,0,length=5)))


###################################################
### chunk number 25: 
###################################################
library(flexmix)
set.seed(2509)
fl2.5 <- stepFlexmix(x2 ~ 1, k = 5, model = FLXmclust(diag = FALSE), nrep=10)


###################################################
### chunk number 26: 
###################################################
parameters(fl2.5, 3, simplify=FALSE)$center


###################################################
### chunk number 27: 
###################################################
parameters(fl2.5, 3, simplify=FALSE)$cov


###################################################
### chunk number 28: 
###################################################
cov2cor(parameters(fl2.5, 3, simplify=FALSE)$cov)


###################################################
### chunk number 29: Figure 11.14
###################################################
plotEll(fl2.5, data=x2, project=p2, eqscale=FALSE)


###################################################
### chunk number 30: Figure 11.15
###################################################
plot(fl2.5, main="", mark=3)


###################################################
### chunk number 31: Figure 11.16
###################################################
par(mfrow=c(2,1))
plotEll(fl2.5, data=x2, which=c(2,4), eqscale=F, points=F)
plotEll(fl2.5, data=x2, which=c(1,3), eqscale=F, points=F)