############################################################################### ## Inhalt: ## Kapitel 7.2 (Elemente Multivariater Statistik) des Kurses von ## Ruckdeschel & Kohl ############################################################################### ############################################################################### ## Nachtrag: (Dis-)Similarity-Matrix ############################################################################### library(hopach) data(golub) ## Genexpressionsdaten zu Leukämie golub.names <- ifelse(golub.cl, "AML", "ALL") # funktioniert nicht unter 2.7.0 Windows! #library(SLmisc) ## Similarity-Matrix basierend auf der Stichproben-Korrelation ## Datei corPlot.R einladen!!! corPlot(golub, new = TRUE, labels = golub.names, minCor = 0.45, title = "Similarity Matrix") ## oder auch farbig library(RColorBrewer) corPlot(golub, new = TRUE, labels = golub.names, minCor = 0.45, title = "Similarity Matrix", col = colorRampPalette(rev(brewer.pal(9,"Blues")))(50)) corPlot(golub, new = TRUE, labels = golub.names, minCor = 0.45, title = "Similarity Matrix", col = colorRampPalette(brewer.pal(9,"RdYlGn"))(50)) ############################################################################### ## Clusteranalyse (unüberwachtes statistisches Lernen) ## Lesen Sie Abschnitt 7.2.6 im Skript von Ruckdeschel & Kohl ## vgl. http://en.wikipedia.org/wiki/Cluster_analysis ############################################################################### ########################################################### ## Hierarchisches Clustering ## vgl. ## http://www.resample.com/xlminer/help/HClst/HClst_intro.htm ## http://www.jmu.edu/docs/sasdoc/sashtml/stat/chap23/sect12.htm ########################################################### ####################################### ## Agglomeratives (aufhäufendes) Hierarchisches Clustering ## Start: Jede Beobachtung ein eigener Cluster ## Sukzessives paarweises Zusammenfassen der Cluster mit kleinstem Abstand ## zueinander. ## Ende: Ein Cluster mit allen Beobachtungen ####################################### ## complete linkage (eukl. Abstand) ## default-Methode ## Abstand zweier Cluster zueinander ist der maximale paarweise Abstand der ## Mitglieder hc1 <- hclust(dist(t(golub)), method = "complete") hc1 str(hc1) plot(hc1, labels = golub.names, main = "complete linkage") ## Identifikation von 4 (?) Clustern cutree(hc1, k = 4) ## Scheiden des Baumes in der Höhe 52 (-> 4 Cluster) cutree(hc1, h = 52) ## auch möglich cutree(hc1, k = 2:5) ## oder cutree(hc1, h = c(40, 45, 50, 55)) ## single linkage (eukl. Abstand) ## Abstand zweier Cluster zueinander ist der minimale paarweise Abstand der ## Mitglieder hc2 <- hclust(dist(t(golub)), method = "single") plot(hc2, labels = golub.names, main = "single linkage") ## average linkage (eukl. Abstand) ## Abstand zweier Cluster zueinander ist der Mittelwert aller paarweisen Abstände ## der Mitglieder (liegt irgendwo zwischen "single" und "complete") hc3 <- hclust(dist(t(golub)), method = "average") plot(hc3, labels = golub.names, main = "average linkage") ## Ward (eukl. Abstand) ## Abstandskriterium ist die minimale Varianz ## (liegt irgendwo zwischen "single" und "complete") hc4 <- hclust(dist(t(golub)), method = "ward") plot(hc4, labels = golub.names, main = "Ward minimum variance method") ## McQuitty hc5 <- hclust(dist(t(golub)), method = "mcquitty") plot(hc5, labels = golub.names, main = "McQuitty") ## Median ## Kein monotones Abstandmaß hc6 <- hclust(dist(t(golub)), method = "median") plot(hc6, labels = golub.names, main = "Median") ## Centroid ## Kein monotones Abstandmaß hc7 <- hclust(dist(t(golub)), method = "centroid") plot(hc7, labels = golub.names, main = "centroid") ## verwende anderen Abstand, z.B. Pearson Stichprobenkorrelation ## Datei corDist.R einladen!!! hc8 <- hclust(corDist(t(golub)), method = "average") plot(hc8, labels = golub.names, main = "average linkage") ## eine leicht andere, flexiblere Darstellungsmöglichkeit erhält man durch ## Umwandlung in ein Dendrogramm D <- corDist(t(golub)) attr(D, "Labels") <- golub.names hc9 <- as.dendrogram(hclust(D, method = "average")) plot(hc9, main = "hclust mit average linkage", ylab = "Pearson Stichprobenkorrelationsabstand") ## Einfärben der Beschriftung des Dendrogramms group <- factor(golub.names) mypalette <- brewer.pal(4, "Dark2") colLab1 <- function(n){ if(is.leaf(n)){ a <- attributes(n) ind <- which(a$label == levels(group)) attr(n, "nodePar") <- c(a$nodePar, list(lab.col = mypalette[ind], lab.cex = 1, pch = "")) } n } dend91 <- dendrapply(hc9, colLab1) plot(dend91, ylab = "Pearson Stichprobenkorrelationsabstand", main = "hclust mit average linkage") ## Zusätzliches Einfärben des Dendrogramms colLab2 <- function(n){ if(is.leaf(n)){ a <- attributes(n) ind <- which(a$label == levels(group)) attr(n, "nodePar") <- c(a$nodePar, list(lab.col = mypalette[ind], lab.cex = 1, pch = 20, col = mypalette[ind])) attr(n, "edgePar") <- c(a$edgePar, list(col = mypalette[ind])) } n } ## durch zusätzliches Setzen von "hang" kann man eine Darstellung ## ähnlich zu plot.hclust erreichen hc92 <- as.dendrogram(hclust(D, method = "average"), hang = 0.1) dend92 <- dendrapply(hc92, colLab2) plot(dend92, ylab = "Pearson Stichprobenkorrelationsabstand", main = "hclust mit average linkage") ## Zusätzlich Berchnung von p-Werten für die Cluster mittels Bootstrap library(pvclust) hc10 <- pvclust(golub, method.hclust="average", method.dist="correlation", nboot = 20) # um Rechenzeit gering zu halten, nboot = 20 plot(hc10) pvrect(hc10) colnames(golub) <- golub.names hc11 <- pvclust(golub[,-21], method.hclust="average", method.dist="correlation", nboot = 20) # um Rechenzeit gering zu halten, nboot = 20 plot(hc11) pvrect(hc11) ####################################### ## Weitere agglomerative hierarchische Clusterverfahren sind z.B. in der ## Funktion "agnes" (AGglomerative NESting) im Paket "cluster", in der ## Funktion "hcluster" im Paket "amap" und in der Funktion "varclus" im ## Paket "Hmisc" implementiert. ## ## Weitere graphische Darstellungsmöglichkeiten finden sich im Paket ## "clusterfly". ####################################### ####################################### ## Divisives (Aufteilendes) Hierarchisches Clustering ## Start: Ein Cluster mit allen Beobachtungen ## Sukzessives Aufteilen des Clusters in Subcluster. ## Ende: Jede Beobachtung ein eigener Cluster ####################################### ## DIvisive ANAlysis Clustering library(cluster) hc12 <- diana(D) plot(hc12, which.plots = 2, main = "Diana") ## Sog. Bannerplot ## Dargestellt ist mittels eines Barplots der Durchmesser eines jeden Clusters, ## der aufgeteilt wurde. par(mfrow=c(1,2)) plot(hc12, which.plots = 1, main = "Bannerplot", nmax.lab = 38) plot(hc12, which.plots = 2, main = "Dendrogramm") ####################################### ## Hyprid-Algorithmen für Hierarchisches Clustering ## Mischung aus agglomerativem und divisivem Verfahren ####################################### ## Hierarchical Ordered Partitioning and Collapsing Hybrid (HOPACH) algorithm ## Verwendet den PAM (partitioning around medoids) Algorithmus; vlg. Funktion ## "pam" im Paket "cluster". library(hopach) hc13 <- hopach(t(golub), d = "cor") str(hc13) ## Anzahl cluster hc13$clustering$k ## Keine Darstellung mittels Dendrogramm möglich. Stattdessen durch nach dem ## Ergebnis des Clusterings geordnete Abstandsmatrix. dplot(distancematrix(t(golub), d = "cor"), hc13, labels = golub.names, main = "HOPACH", col = colorRampPalette(rev(brewer.pal(9,"RdYlGn")))(50)) ####################################### ## Modell-basiertes Hierarchisches Clustering ####################################### library(mclust) ## EII: spherical, equal volume hc14 <- hc(modelName = "EII", data = t(golub)) cbind(golub.names, hclass(hc14, 2:5)) ## EII: spherical, equal volume hc15 <- hc(modelName = "VII", data = t(golub)) cbind(golub.names, hclass(hc15, 2:5)) ## EEE: ellipsoidal, equal volume, shape, and orientation hc16 <- hc(modelName = "EEE", data = t(golub)) cbind(golub.names, hclass(hc16, 2:5)) ## VVV: ellipsoidal, varying volume, shape, and orientation hc17 <- hc(modelName = "EII", data = t(golub)) cbind(golub.names, hclass(hc17, 2:5)) ####################################### ## Eine spezielle Form der kombinierten Darstellung von Daten und Clustering- ## ergebnissen ist die sog. heatmap. ## Die Funktion "heatmap.2" aus dem Paket "gplots" bietet mehr Möglichkeiten ## als die Funktion "heatmap" aus dem Paket "stats". ####################################### library(gplots) ## Farbe für die Daten ## Datei heatmapCol.R einladen!!! farbe <- heatmapCol(data = golub, col = rev(colorRampPalette(brewer.pal(10, "RdYlBu"))(256)), lim = 1.5) ## Farbe für AML und ALL mypalette <- brewer.pal(3, "Dark2")[c(1,2)] farb.col <- as.integer(factor(golub.names)) dev.off() heatmap.2(golub, col = farbe, scale = "none", ColSideColors = mypalette[farb.col], dendrogram = "column", tracecol = "black", labRow = "", trace = "none", cexRow = 0.4, cexCol = 1, Rowv = FALSE, margins = c(4, 2), main = "Golub-Daten", hclustfun = function(d) hclust(d, method = "average"), dist = corDist) ## "Sortieren" der Zeilen (Gene) mittels hclust (jedoch kein dendrogram) ## Hohe Rechenzeit! heatmap.2(golub, col = farbe, scale = "none", ColSideColors = mypalette[farb.col], dendrogram = "column", tracecol = "black", labRow = "", trace = "none", cexRow = 0.4, cexCol = 1, Rowv = TRUE, margins = c(4, 2), main = "Golub-Daten", hclustfun = function(d) hclust(d, method = "average"), dist = corDist) ########################################################### ## k-means ## vgl. http://en.wikipedia.org/wiki/K-means_algorithm ########################################################### ## Animation zur Funktionsweise von kmeans library(animation) x <- matrix(runif(200), ncol = 2) kmeans.ani(x = x, centers = 2, interval = 0.5) ## Ergebnis hängt ab von den zufällig gewählten Startpunkten kmeans.ani(x = x, centers = 2, interval = 0.5) ## Mehr Cluster kmeans.ani(x = x, centers = 4, interval = 0.5) ## "Standardimplementation" ist Funktion "kmeans" aus dem Paket "stats" km1 <- kmeans(t(golub), centers = 3, iter.max = 20, nstart = 10) str(km1) km1$cluster ####################################### ## Wie viele Cluster sind in den Daten enthalten? ####################################### ## Calinski-Harabasz-Kriterium ## Funktioniert gut, wenn die Cluster alle etwa gleich groß sind library(vegan) casc1 <- cascadeKM(t(golub), inf.gr=2, sup.gr=10, iter = 10, criterion = "calinski") plot(casc1) ## simple structure index ## Funktioniert gut, wenn die Cluster alle etwa gleich groß sind casc2 <- cascadeKM(t(golub), inf.gr=2, sup.gr=12, iter = 10, criterion = "ssi") plot(casc2) ## Median (or Mean) Split Silhouette (MSS) ## Median library(hopach) msscheck(D, graph = TRUE) # 5 cluster ## Mean msscheck(D, within = "mean", between = "mean", graph = TRUE) # 3 cluster ## Silhouettes silcheck(D, graph = TRUE, diss = TRUE) # 4 cluster ## Für einige weitere Möglichkeiten siehe Paket "clusterSim" ####################################### ## Weiteres zu k-means: ## 1. k-means unter Verwendung verschiedener Distanzen: ## vgl. Funktion "Kmeans" im Paket "amap" ## ## 2. Trimmed k-means Clustering ## Versuch der Robustifizierung von k-means ## Implementiert in der Funktion "trimkmeans" im Paket "trimcluster". ## ## 3. k-means Clustering unter Verwendung von Kernen ## Problem: Das klassische k-means ist nicht in der Lage, Cluster zu erkennen, ## die nicht linear separabel im Eingaberaum sind. ## "kernel trick": Mittels spezieller Funktionen, sog. Kernfunktionen, werden ## die Daten transformiert bzw. in einen anderen Raum projeziert. In diesem ## neuen Raum sind dann mglweise nicht linear trennbare Cluster linear ## trennbar und können von k-means "erkannt" werden. ####################################### ## Kombination von hclust und kmeans im Rahmen von heatmap.2 km <- kmeans(golub, centers = 6, iter.max = 20, nstart = 10) golub.ord <- golub[order(km$cluster),] k <- length(table(km$cluster)) farb.row <- brewer.pal(k, "Set1")[sort(km$cluster)] heatmap.2(golub.ord, col = farbe, scale = "none", ColSideColors = mypalette[farb.col], RowSideColors = farb.row, dendrogram = "column", tracecol = "black", labRow = "", trace = "none", cexRow = 0.4, cexCol = 1, Rowv = FALSE, margins = c(4, 2), main = "Golub-Daten", labCol = golub.names, hclustfun = function(d) hclust(d, method = "average"), dist = corDist) ########################################################### ## k-medoids (PAM, clara) ## vgl. http://rosuda.org/lehre/SS07/MVstat-f/MSsClus2.pdf ########################################################### ## Betrachte Medoide, d.h., Repräsentanten für die k-verschiedenen Cluster ## Ordne die anderen Beobachtungen den Medoiden zu. med1 <- pam(D, k = 4) ## Für große Datensätze wird "clara" anstelle von "pam" empfohlen med2 <- clara(golub, k = 10) ## Anderer Algorithmus library(clue) med3 <- kmedoids(D, k = 4) str(med3) ########################################################### ## Self-Organizing Maps (SOM) ## vgl. ## http://de.wikipedia.org/wiki/Selbstorganisierende_Karte ## http://en.wikipedia.org/wiki/Self-organizing_map ########################################################### library(som) som1 <- som(t(golub), xdim = 1, ydim = 3) plot(som1) ## weitere Darstellungsmöglichkeit library(chemometrics) plotsom(som1, grp = as.factor(golub.names), type = "bar") som2 <- som(t(golub), xdim = 3, ydim = 2) plotsom(som2, grp = as.factor(golub.names), type = "bar") ########################################################### ## Es gibt eine Vielzahl weiterer Clusterverfahren. ## vgl. ## http://de.wikipedia.org/wiki/Clusteranalyse ## http://en.wikipedia.org/wiki/Cluster_analysis ########################################################### ############################################################################### ## Übungsaufgabe 3: ## Verwenden Sie den Datensatz "banknote" aus dem Paket "dr". ## (a) Lassen Sie die Variable Y aus dem Datensatz "banknote" weg und ## führen Sie für den Datensatz verschiedene hierarchische Clusteranalysen ## durch. ## (b) Stellen Sie eines der Ergebnisse mittels eines Dendrogramms graphisch ## dar, wobei Sie die Werte, die zu echten Banknoten (d.h., Y = 0) gehören ## grün einfärben und Werte, die zu gefälschten Banknoten gehören (d.h., Y = 1) ## rot einfärben. ## (c) Verwenden Sie eine Heatmap, um den Datensatz und das Clustering ## darzustellen. ## (d) Untersuchen Sie mit Hilfe geeigneter Kriterien (Indizes), wieviele ## Cluster in den Daten vorhanden sind. ## (e) Cluster Sie die Daten mittels k-means, k-medoids und self-organizing ## maps. ############################################################################### ############################################################################### ## Lösungen: ############################################################################### ########################################################### ## Übungsaufgabe 3 ########################################################### ####################################### ## Teil (a) ####################################### library(dr) data(banknote) D.dr <- dist(banknote[,-7]) attr(D.dr, "Labels") <- as.character(!as.logical(banknote[,7])) hc.dr <- hclust(D.dr, method = "average") plot(hc.dr, cex = 0.5) ## oder dend.dr <- as.dendrogram(hc.dr) plot(dend.dr, main = "hclust mit average linkage", ylab = "Euklidischer Abstand") ## oder library(hopach) library(RColorBrewer) hc.dr2 <- hopach(banknote[,-7], d = "euclid") dplot(distancematrix(banknote[,-7], d = "euclid"), hc.dr2, labels = as.character(!as.logical(banknote[,7])), main = "HOPACH", col = colorRampPalette(rev(brewer.pal(9,"RdYlGn")))(50), cex.axis = 0.5) ## oder hc.dr3 <- hc(modelName = "VII", data = banknote[,-7]) cbind(as.character(!as.logical(banknote[,7])), hclass(hc.dr3, 2:5)) ####################################### ## Teil (b) ####################################### group <- factor(as.character(!as.logical(banknote[,7]))) mypalette <- brewer.pal(4, "Dark2")[c(2,1)] colLab2 <- function(n){ if(is.leaf(n)){ a <- attributes(n) ind <- which(a$label == levels(group)) attr(n, "nodePar") <- c(a$nodePar, list(lab.col = mypalette[ind], lab.cex = 0.5, pch = "", col = mypalette[ind])) attr(n, "edgePar") <- c(a$edgePar, list(col = mypalette[ind])) } n } hc.dr4 <- as.dendrogram(hclust(D.dr, method = "average"), hang = 0.1) dend.dr2 <- dendrapply(hc.dr4, colLab2) plot(dend.dr2, ylab = "Euklidischer Abstand", main = "hclust mit average linkage") ####################################### ## Teil (c) ####################################### library(gplots) banknote.scal <- t(t(banknote[,-7])-colMeans(banknote[,-7])) farbe <- heatmapCol(data = banknote.scal, col = rev(colorRampPalette(brewer.pal(10, "RdYlBu"))(256)), lim = 1.5) heatmap.2(banknote.scal, col = farbe, scale = "none", dendrogram = "both", tracecol = "black", labRow = as.character(!as.logical(banknote[,7])), trace = "none", cexRow = 0.4, cexCol = 1, Rowv = TRUE, margins = c(4, 2), main = "Banknote-Daten", hclustfun = function(d) hclust(d, method = "average"), dist = dist) ####################################### ## Teil (d) ####################################### ## Calinski-Harabasz-Kriterium library(vegan) casc1 <- cascadeKM(banknote[,-7], inf.gr=2, sup.gr=10, iter = 10, criterion = "calinski") plot(casc1) # 2 cluster ## simple structure index ## Funktioniert gut, wenn die Cluster alle etwa gleich groß sind casc2 <- cascadeKM(banknote[,-7], inf.gr=2, sup.gr=10, iter = 10, criterion = "ssi") plot(casc2) # 9 cluster ## Median (or Mean) Split Silhouette (MSS) ## Median library(hopach) msscheck(D.dr, graph = TRUE) # 6 cluster ## Mean msscheck(D.dr, within = "mean", between = "mean", graph = TRUE) # 5 cluster ## Silhouettes silcheck(D.dr, graph = TRUE, diss = TRUE) # 7 cluster ####################################### ## Teil (e) ####################################### ## k-means km.dr <- kmeans(banknote[,-7], centers = 2, iter.max = 20, nstart = 10) ## k-medoids med.dr1 <- pam(D.dr, k = 2) ## SOM library(som) som.dr1 <- som(banknote[,-7], xdim = 1, ydim = 2) library(chemometrics) plotsom(som.dr1, grp = as.factor(banknote[,7]), type = "bar") som.dr2 <- som(banknote[,-7], xdim = 2, ydim = 2) plotsom(som.dr2, grp = as.factor(banknote[,7]), type = "bar")