Hybridisation of K-Means Clustering Analysis and Principal Component Analysis

Perhaps the most important issue when performing K-Means Clustering Analysis, or any clustering analysis for that matter, would be how well clusters are formed against the given variables within the data set. In some cases, outliers may need to be eliminated before performing clustering analysis as they sometimes interfere with forming best clusters.

This hybridisation technique aims to manipulate the data so that more distinct clusters are formed without losing any data points. Principal components are used in multidimensional scaling when performing clustering analysis. In this technique, the variables are normalised before principal components are derived. This has an effect of congregating data points into clusters more densely, hence separating clusters more distinctly.

The technique used in this blog modified the technique described in the below paper.

A hybridized K-means clustering approach for high dimensional dataset; Rajashree Dash, Debahuti Mishra, Amiya Kumar Rath, Milu Acharya; International Journal of Engineering, Science and Technology Vol. 2, No. 2, 2010, pp. 59-66; http://ijest-ng.com/ijest-ng-vol2-no2-pp59-66.pdf

To demonstrate the effect of this technique, iris data set was used.

Dat <- iris[, 1:4] 

CTR <- aggregate(. ~ Species, iris, mean)

#clustering on the original data 

km2 <- kmeans(Dat, centers = CTR[, 2:5], iter.max = 100) 

#clustering on the principal components

PC <- princomp(Dat) 

PC2 <- as.data.frame(PC$score[, 1:2]) 

PC2$Species <- iris$Species  

CTR_PC2 <- aggregate(. ~ Species, PC2, mean) 

kmPC <- kmeans(PC2[, 1:2], centers = CTR_PC2[, 2:3], iter.max = 100)

The below graph compares clusters formed by the original data set to clusters formed against its principal components. It is clearly shown that clusters formed against principal components are better congregated and there are less overlaps between clusters in 2 dimensional display.

par(mfrow = c(3,1))

plot(Dat[, 1:2], col = km2$cluster, main = "Original Iris Data Set - 1st & 2nd variables") points(km2$centers[, 1:2], col = c(1:3), pch = 16, cex = 2) 

text(km2$centers[, 1:2], labels = c(1:3), col = c(1:3), pos = 3, cex = 2)

plot(Dat[, 3:4], col = km2$cluster, main = "Original Iris Data Set - 3rd & 4th variables") points(km2$centers[, 3:4], col = c(1:3), pch = 16, cex = 2) 

text(km2$centers[, 3:4], labels = c(1:3), col = c(1:3), pos = 3, cex = 2)

plot(PC2[, 1:2], col = kmPC$cluster, xlab = "PC1", ylab = "PC2", main = "1st & 2nd Principal Components") 

points(kmPC$centers, col = c(1:3), pch = 16, cex = 2) 

text(kmPC$centers, labels = c(1:3), col = c(1:3), pos = 3, cex = 2)

In this hybridisation technique, we are normalising the variables prior to deriving principal components as shown below. Then, only the principal components with respective eigen values greater than the average eigen values of all principal components are used in the clustering analysis.

V <- apply(Dat, 2, var) Input <- Dat[, which(V > 0)] Mean <- apply(Input, 2, mean) Sdev <- apply(Input, 2, sd) Adj.Data <- t(apply(Input, 1, function(x) (x - Mean))) Norm.Data <- t(apply(Adj.Data, 1, function(x) (x/Sdev))) Cov <- var(Norm.Data) Eig.Vec <- eigen(Cov) Featured <- Eig.Vec[[2]][, 1:2] temp <- t(Featured) %*% t(Norm.Data) test1 <- t(temp) DatN <- cbind(as.data.frame(test1), Species = iris$Species) CTR_PCN <- aggregate(. ~ Species, DatN, mean) kmPCN <- kmeans(test1, centers = CTR_PCN[, 2:3], iter.max = 100)

The below compares clusters formed by original variables, by principal components and by principal components of normalised variables. The main difference with normalisation is that the outliers become more distinct when clusters are formed in this example.

colnames(Dat) <- paste("Col", c(1:4), sep = "") Dat$Species <- iris$Species Dat$Cluster <- km2$cluster Dat$Dist <- apply(Dat, 1, function(x) sqrt(sum((as.numeric(x[1:4]) - km2$center[as.numeric(x[6]), ])^2))) Dat$Ind <- "Orig" colnames(PC2)[1:2] <- paste("Col", c(1:2), sep = "") PC2$Cluster <- kmPC$cluster PC2$Dist <- apply(PC2, 1, function(x) sqrt(sum((as.numeric(x[1:2]) - kmPC$center[as.numeric(x[4]),])^2))) PC2$Col3 <- NA PC2$Col4 <- NA PC2$Ind <- "PCA-2D" colnames(DatN)[1:2] <- paste("Col", c(1:2), sep = "") DatN$Cluster <- kmPCN$cluster DatN$Dist <- apply(DatN, 1, function(x) sqrt(sum((as.numeric(x[1:2]) - kmPCN$center[as.numeric(x[4]),])^2))) DatN$Col3 <- NA DatN$Col4 <- NA DatN$Ind <- "PCA-Norm" DatX <- rbind(Dat, PC2, DatN) DatX1 <- split(DatX, DatX$Ind) DatX1 <- lapply(DatX1, function(x) { x$Dist <- x$Dist/max(x$Dist) x })

DatX1 <- ldply(DatX1, data.frame)

par(mfrow = c(3, 1)) plot(Dat[, 1:2], col = km2$cluster, main = "Original Iris Data Set - 1st & 2nd variables") points(km2$centers[, 1:2], col = c(1:3), pch = 16, cex = 2) text(km2$centers[, 1:2], labels = c(1:3), col = c(1:3), pos = 3, cex = 2)

plot(PC2[, 1:2], col = kmPC$cluster, xlab = "PC1", ylab = "PC2", main = "1st & 2nd Principal Components") points(kmPC$centers, col = c(1:3), pch = 16, cex = 2) text(kmPC$centers, labels = c(1:3), col = c(1:3), pos = 3, cex = 2)

plot(DatN[, 1:2], col = kmPCN$cluster, xlab = "PC1", ylab = "PC2", main = "Normalised Variables - 1st & 2nd Principal Components") points(kmPCN$centers, col = c(1:3), pch = 16, cex = 2) text(kmPCN$centers, labels = c(1:3), col = c(1:3), pos = 3, cex = 2)

The below examines the density of each cluster when different methods are applied. It can be seen that normalised principal components generally forms denser clusters.

direct.label(densityplot(~Dist | paste("Cluster", Cluster, sep = " "), DatX1, groups = DatX1$Ind,layout = c(3, 1)))

Drawing Great Circle in R - Flight Map

'Great Circle' is usually drawn on a map of world to display flight paths between cities. This is useful when there are multiple flights between same cities or in the nearby cities where the straight lines often overlap and do not display the fact that there are more than one flights represented.

Firstly, the packages I used include:

• 'maps'
• 'plyr'
• 'sp'
• 'geosphere'

The data I have used in the below examples came from http://openflights.org/data.html which provided information about airlines, airports and flight routes in separate files. The data is structured in a way that makes merging fairly easy.

airline <- read.csv("Data/Airlines.csv", sep = ",", header = TRUE) airport <- read.csv("Data/Airport.csv", sep = ",", header = TRUE) route <- read.csv("Data/Routes.csv", sep = ",", header = TRUE) A <- airport[, c("ID", "Airport", "City", "Country", "Lat", "Lon")] colnames(A) <- paste("Dep", colnames(A), sep = "_") Dat <- merge(route[, c("Airline_ID", "Dep_Airport_ID", "Arr_Airport_ID")], A, by.x = "Dep_Airport_ID", by.y = "Dep_ID") A <- airport[, c("ID", "Airport", "City", "Country", "Lat", "Lon")] colnames(A) <- paste("Arr", colnames(A), sep = "_") Dat <- merge(Dat, A, by.x = "Arr_Airport_ID", by.y = "Arr_ID") Dat <- merge(Dat, airline[, c("ID", "Airline", "Active")], by.x = "Airline_ID", by.y = "ID")

As it would make little sense to draw all flight paths around the world, I have filtered the data to selected countries to be used in the example.

DatX <- Dat[which(Dat$Dep_Country != Dat$Arr_Country), ] DatX <- DatX[which(DatX$Dep_Country %in% c("Canada", "Australia", "South Korea", "United States", "Japan") & DatX$Arr_Country %in% c("Canada", "Australia", "South Korea", "United States", "Japan")), ]

To distinguish between flight paths, I have used the colour degradation between red and white, as there are too many flight paths to colour code them with different colours.

col_gradient <- colorRampPalette(c("white", "red")) Colours <- col_gradient(length(unique(DatX$Airline_ID))) DatX$Colour <- Colours[as.numeric(factor(as.character(DatX$Airline_ID)))]

The below code calculates the projection of great circles to correspond to flight routes. In essence, I have assigned 100 points between the destination cities to connect in order to form a great circle, and repeated this for all flight destinations in the dataset. The tricky part is that depending on the centre of the world map (e.g. Pacific-centred or Atlantic-centred), you may need to recalculate the longitude as the map does not show negative longitudes and may result in disconnected or distorted great circles .

GC_Dat <- list()

for (j in 1:nrow(DatX)) {
    X <- as.data.frame(gcIntermediate(DatX[j, c("Arr_Lon", "Arr_Lat")], DatX[j, 
        c("Dep_Lon", "Dep_Lat")], n = 100, addStartEnd = TRUE))

    R <- nrow(X[which(X$lon < 0), ])
    if (R == 0) {
        RR <- X
    }
    if (R > 0) {
        a <- ifelse(X$lon >= 0, X$lon, X$lon + 360)

        DF <- diff(a)
        L <- length(DF) - length(DF[which(DF < 0)])

        if (L == 0 | L == length(DF)) {
            X$lon <- a
            RR <- X
        }
        if (L != 0 & L != length(DF)) {
            aa <- X[which(X$lon >= 0), ]
            b <- X[which(X$lon < 0), ]
            if (nrow(b) > 0) {
                B <- b[nrow(b), ]

                b$lon <- b$lon + 360
                B$lon <- 360

                b <- rbind(b, B)
                b$id <- b$id + 0.5
            }

            RR <- rbind(aa, b)
        }
    }

    RR$Airline_ID <- DatX[j, "Airline_ID"]

    GC_Dat[[j]] <- merge(DatX[j, ], RR, by = "Airline_ID")
}

Now that the data preparation is completed, we need to plot the world map as a below.

map("world2", col = "beige", bg = "black")



To make it look pretty, the following codes will show top 5% cities by population on the map.

data(world.cities) Cities <- world.cities Maj_Cities <- Cities[which(Cities$pop >= quantile(Cities$pop, prob = 0.95)), ] Maj_Cities$long <- ifelse(Maj_Cities$long >= 0, Maj_Cities$long, Maj_Cities$long + 360) map("world2", col = "beige", bg = "black") points(Maj_Cities$long, Maj_Cities$lat, pch = 16, cex = 0.3, col = "gold")



Finally, the flight routes as represented by the great circles are laid on top of the map.

map("world2", col = "beige", bg = "black") for (i in 1:length(GC_Dat)) { lines(lat ~ lon, GC_Dat[[i]], lty = 1, col = unique(GC_Dat[[i]]$Colour)) } points(Maj_Cities$long, Maj_Cities$lat, pch = 16, cex = 0.3, col = "gold") title(sub = "Flight Routes", col.sub = "white")