In order to illustrate hierarchical clustering techniques and k-means, I did borrow François Husson‘s dataset, with monthly average temperature in several French cities.
> temp=read.table(
+ "http://freakonometrics.free.fr/FR_temp.txt",
+ header=TRUE,dec=",")
We have 15 cities, with monthly observations
> X=temp[,1:12]
> boxplot(X)
Since the variance seems to be rather stable, we will not ‘normalize’ the variables here,
> apply(X,2,sd)
Janv Fevr Mars Avri
2.007296 1.868409 1.529083 1.414820
Mai Juin juil Aout
1.504596 1.793507 2.128939 2.011988
Sept Octo Nove Dece
1.848114 1.829988 1.803753 1.958449
In order to get a hierarchical cluster analysis, use for instance
> h <- hclust(dist(X), method = "ward")
> plot(h, labels = rownames(X), sub = "")
An alternative is to use
> library(FactoMineR)
> h2=HCPC(X)
> plot(h2)
Here, we visualise observations with a principal components analysis. We have here also an automatic selection of the number of classes, here 3. We can get the description of the groups using
> h2$desc.ind
or directly
> cah=hclust(dist(X))
> groups.3 <- cutree(cah,3)
We can also visualise those classes by ourselves,
> acp=PCA(X,scale.unit=FALSE)
> plot(acp$ind$coord[,1:2],col="white")
> text(acp$ind$coord[,1],acp$ind$coord[,2],
+ rownames(acp$ind$coord),col=groups.3)
It is possible to plot the centroïds of those clusters
> PT=aggregate(acp$ind$coord,list(groups.3),mean)
> points(PT$Dim.1,PT$Dim.2,pch=19)
If we add Voroid sets around those centroïds, here we do not see them (actually, we see the point – in the middle – that is exactly at the intersection of the three regions),
> library(tripack)
> V <- voronoi.mosaic(PT$Dim.1,PT$Dim.2)
> plot(V,add=TRUE)
To visualize those regions, use
> p=function(x,y){
+ which.min((PT$Dim.1-x)^2+(PT$Dim.2-y)^2)
+ }
> vx=seq(-10,12,length=251)
> vy=seq(-6,8,length=251)
> z=outer(vx,vy,Vectorize(p))
> image(vx,vy,z,col=c(rgb(1,0,0,.2),
+ rgb(0,1,0,.2),rgb(0,0,1,.2)))
> CL=c("red","black","blue")
> text(acp$ind$coord[,1],acp$ind$coord[,2],
+ rownames(acp$ind$coord),col=CL[groups.3])
Actually, those three groups (and those three regions) are also the ones we obtain using a k-mean algorithm,
> km=kmeans(acp$ind$coord[,1:2],3)
> km
K-means clustering
with 3 clusters of sizes 3, 7, 5
(etc). But actually, since again we have some spatial data, it is possible to visualize them on a map
> library(maps)
> map("france")
> points(temp$Long,temp$Lati,col=groups.3,pch=19)
or, to visualize the regions, use e.g.
> library(car)
> for(i in 1:3)
+ dataEllipse(temp$Long[groups.3==i],
+ temp$Lati[groups.3==i], levels=.7,add=TRUE,
+ col=i+1,fill=TRUE)
Those three regions actually make sense, geographically speaking.
