Clusteranalyse in R

Im ersten Teil des Blogs haben wir die theoretischen Grundlagen der Clusteranalyse näher beleuchtet. Im Folgenden geht es nun darum, die Theorie mithilfe der Statistikumgebung R in die Praxis umzusetzen. Für die Analyse in R werden wir die Variablen (Kraftstoffverbrauch in miles per gallon) und (Hubraum in Kubikinch) für 32 Automobile aus dem Datensatz verwenden. Eine ausführliche Beschreibung des Datensatzes finden Sie hier.

require(dplyr)
require(broom)
require(ggplot2)
require(ggdendro)
require(cluster)

Zunächst wird ein neuer Datensatz angelegt, der die genannten Variablen und für die 32 Autos enthält. Die Beobachtungen werden in einem Streudiagramm dargestellt.

car <- select(mtcars, mpg, disp)
ggplot(car, aes(mpg, disp)) + geom_point()

Streudiagramm

Nun führen wir ein hierarchisches Verfahren durch und schauen uns das Dendrogramm an. Auf der y-Achse des Dendrogramms ist die Distanz (euklidischer Abstand) abgetragen.

h.cluster <- car %>% dist(., method = "euclidean") %>% hclust(., method = "ward.D")
ggdendrogram(h.cluster)

Dendrogramm

Augenscheinlich sind 2 Cluster in diesem Falle sinnvoll. Als nächstes erstellen wir einen Screeplot, wobei für jedes ein k-means Verfahren durchgeführt wird. Der Screeplot zeigt auf der x-Achse die Anzahl der Cluster und auf der y-Achse die Varianz innerhalb der Cluster.

multi.clust <- data.frame(k = 1:6) %>% group_by(k) %>% do(clust = kmeans(car, .$k))
sumsq.clust <- multi.clust %>% group_by(k) %>% do(glance(.$clust[[1]]))

ggplot(sumsq.clust, aes(k, tot.withinss)) + geom_line() + geom_point()

Screeplot

Der Knick befindet sich bei , entsprechend legen wir uns auf 2 Cluster fest. Nun sehen wir uns die durch den k-means-Algorithmus erhaltene Gruppenstruktur für an, die Clusterzentren sind mit einem * gekennzeichnet.

p.cluster <- car %>% kmeans(., 2)
p.cluster$cluster <- as.factor(p.cluster$cluster)

ggplot(car, aes(mpg, disp, label = rownames(car))) + 
  scale_fill_discrete(name = "Cluster") + xlim(9,35) +
  geom_label(aes(fill = p.cluster$cluster), colour = "white", 
  fontface = "bold", size=2)

Cluster k = 2

An diesem Punkt wäre es interessant, vergleichsweise die Gruppenstruktur auch für , , etc. zu sehen.

multi.clust <- data.frame(k = 1:6) %>% group_by(k) %>% do(clust = kmeans(car, .$k))
multi.k <- multi.clust %>% group_by(k) %>% do(augment(.$clust[[1]], car))

ggplot(multi.k, aes(mpg, disp)) + geom_point(aes(color = .cluster)) + 
  facet_wrap(~k)

Wie im ersten Teil des Blogs zum Thema “Clusteranalyse” erwähnt, lässt sich die Gruppenzuordnung gut mit einem Silhouettenplot beurteilen.

sil <- car %>% dist(.) %>% silhouette(p.cluster$cluster, .)
sil.data <- data.frame(cluster = factor(sil[, 1]), sil_width = sil[, 3])

ggplot(sil.data, aes(x = row.names(sil.data), y = sil_width, fill = cluster)) +
  geom_bar(stat = "identity", width = 0.5) + coord_flip() + labs(x = "") +
  scale_x_discrete(limits = row.names(sil.data[order(sil.data$cluster, 
  sil.data$sil_width), ]))

Silhouettenplot

mean(sil.data$sil_width)

## [1] 0.687796

Keine der Silhouettenbreiten hat einen Wert . Darüber hinaus beträgt die durchschnittliche Silhouettenbreite 0.69, womit wir die Zuordnung der Beobachtungen zu den zwei Clustern als gut betiteln können.

Interpretation der Clusterlösung

Nachdem wir die Clusterlösung validiert haben, kommen wir nun zum eigentlich schwierigsten Teil der Analyse: der Interpretation. Für die Analyse haben wir den Kraftstoffverbrauch und den Hubraum herangezogen. Kurze Information für alle, die auch keine Ahnung von Autos haben: je größer der Hubraum eines Fahrzeugs, desto leistungsstärker ist der Motor.

Zur Hilfestellung nehmen wir die restlichen Informationen aus dem Datensatz . Dazu übergeben wir die Spalte mit der Clusterzuordnung für jedes Auto an den ursprünglichen Datensatz zurück.

clustcars <- cbind(mtcars, cluster = p.cluster$cluster)

Nun können wir die restlichen Variablen getrennt nach Cluster untersuchen.

round(clustcars %>% 
        group_by(cluster) %>% 
        summarise(avg_cyl = mean(cyl), 
                  avg_horsepower = mean(hp), 
                  avg_drat = mean(drat),
                  avg_weight = mean(wt),
                  avg_qsec = mean(qsec),
                  share_vmotor = sum(vs)/nrow(clustcars),
                  share_manualtransmission = sum(am)/nrow(clustcars),
                  avg_gears = mean(gear),
                  avg_carb = mean(carb)), 1)

  cluster avg_cyl avg.horsepower avg_drat avg_weight avg_qsec share_vmotor
    (dbl)   (dbl)          (dbl)    (dbl)      (dbl)    (dbl)        (dbl)
1       1     7.9          202.6      3.2        3.9     16.9          0.0
2       2     4.7           97.4      3.9        2.6     18.6          0.4

  cluster share_manual_transmission avg_gears avg_carb
    (dbl)                     (dbl)     (dbl)    (dbl)
1       1                       0.1       3.3      3.3
2       2                       0.3       4.1      2.4

Zunächst halten wir aus der oben gefundenen Clusterlösung fest, dass die Automobile, die dem Cluster 1 zugeordnet wurden, generell eine höhere Leistung und einen höheren Verbrauch als die Autos aus Cluster 2 haben.

Aus der Tabelle wird ersichtlich, dass Cluster 1-Autos durchschnittlich knapp 8 Zylinder haben, während im Cluster 2 sich typischerweise Autos mit 4 bis 6 Zylindern befinden. Mit den Zylindern steigt i.d.R. die Leistung (Variable - sie ist für die Autos im Cluster 1 im Durchschnitt doppelt so hoch wie für die im Cluster 2) , aber auch der Preis und der Kraftstoffverbrauch.

Mit den gesammelten Informationen interpretieren wir die zwei gefundenen Cluster als zwei Preisklassen. Das Cluster 1 enthält die luxuriöseren Autos, mit mehr Leistung, mehr Zylindern und höherem Verbrauch, man könnte die Gruppe als 'Fahrzeuge der Oberklasse' betiteln. Das Cluster 1 enthält demnach leistungsschwächere Autos, die günstiger im Preis sind und einen geringeren Verbrauch aufweisen. Um diese Interpretation abzusichern, haben wir die unverbindliche Preisempfehlung (Achtung: die Modelle stammen aus den 1970ern) für die Autos herausgesucht. In der unteren Tabelle sind die Durchschnittswerte für jedes der beiden Cluster angegeben, in Klammern steht zusätzlich der Median.