################################################ ### Klusteranalyse ### ################################################ ################################################## # A.- spielen wir ein bisschen mit Klusteranalyse # 1.- Daten laden: data(iris) colnames(iris) # 2. Matrix von Interdistanzen rechnen: dd = dist(iris[,-5]) # ACHTUNG: 5. Spalte ist ein Faktor, nicht mitnehmen # 3.- hyerarchische Klusteranalyse: ? hclust # Hilfeseite lesen hh = hclust(dd, method="ward.D") plot(hh) # 4.- Interpretation: grW = cutree(hh, k=3) # wieviele Gruppen soll man kriegen? table(grW, iris[,5]) plot(iris[,-5], col=grW, bg=iris[,5], pch=22) # Uebung: # Testen Sie andere "method"-Werte in hclust. # Speichern Sie die Gruppen jedes Mal in einem neuen Object, e.g. grW2, grCL, grSL usw. # Sind die Gruppen stabil? Welche Methode unterscheiden sind? Welche liefern ähnliche Ergebnisse? ################################################## # B.- jetzt mit simulierten Daten # 1.- Daten simulieren: require(mvtnorm) # wenn es nicht klappt, install.packages("mvtnorm") zuerst (S = matrix(c(1,-0.5,-0.5,1), ncol=2)) # ist S symmetrisch positivdefiniert? alpha = 0.25 set.seed(10054) # um sicher zu gehen, dass wir alle die gleich Ergebnisse bekommen X1 = rmvnorm(30, mean=c(2,2), sigma=alpha*S) X2 = rmvnorm(30, mean=c(0,0), sigma=alpha*S) X3 = rmvnorm(30, mean=c(-1,2), sigma=alpha*S) X = rbind(X1, X2, X3) gr0 = rep(1:3, each=30) plot(X, col=gr0) # Danach, spielen mit alpha, Sigma und mit mean von X3, z.B # X3 = rmvnorm(30, mean=c(1,2), sigma=alpha*abs(S)) # 2. Matrix von Interdistanzen rechnen: dd = dist(X) # 3.- massive hyerarchische Klusteranalyse: gr = list() par(mfrow=c(2,4)) for(mymethod in c("ward.D", "ward.D2", "single", "complete", "average", "mcquitty", "median","centroid")){ hh = hclust(dd, method=mymethod) plot(hh, main=mymethod) h0 = locator(1) gr[[mymethod]] = rect.hclust(hh, h=h0$y) } # 4.- Gruppierungen vergleichen: print(gr) ################################################## # C.- kmeans für die simulierten Daten # 1.- Daten simulieren: require(mvtnorm) # wenn es nicht klappt, install.packages("mvtnorm") zuerst (S = matrix(c(1,-0.5,-0.5,1), ncol=2)) # ist S symmetrisch positivdefiniert? alpha = 0.25 set.seed(10054) # um sicher zu gehen, dass wir alle die gleich Ergebnisse bekommen X1 = rmvnorm(30, mean=c(2,2), sigma=alpha*S) X2 = rmvnorm(30, mean=c(0,0), sigma=alpha*S) X3 = rmvnorm(30, mean=c(-1,2), sigma=alpha*S) X = rbind(X1, X2, X3) gr0 = rep(1:3, each=30) plot(X, col=gr0) points(2,2, col=1, pch=4, cex=2) points(0,0, col=2, pch=4, cex=2) points(-1,2, col=3, pch=4, cex=2) # 2.- kmeans algorithm km = kmeans(X, centers=3) ?km # Hilfeseite lesen!! # 3.- Interpretation table(gr0, km$cluster) points(km$centers, pch=19) plot(X, col=gr0, bg=km$cluster, pch=21, lwd=4, cex=1.5) # Schliessen sie die Kmeans-Analyse in eine Schleife um, um die verschiedene Algorithmen zu vergleichen par(mfrow=c(2,2)) for(myalgorithm in c(XXXX)){ } # TIPP: man braucht nur der letzte Plot (Zeile 105)