1.1 Motivation: Clusteranalyse mit Hilfe Gaußscher Mischverteilungen

Als Datenbeispiel verwendent wir die palmerpenguins Daten (Horst, Hill, und Gorman 2020).

Diese Daten stammen aus Vermessungen von Pinguinpopulationen auf dem Palmer-Archipel (Antarktische Halbinsel). Pinguine sind oft schwer von einander zu unterscheiden (Abbildung 1.1). Wir werden versuchen, mit Hilfe einer Gaußschen Mischverteilung Gruppierungen in den Pinguindaten (Flossenlänge) zu finden. Um solche Mischverteilungen schätzen zu können, führen wir den EM Algorithmus ein.

Frecher Pinguin bei der Tat.

Abbildung 1.1: Frecher Pinguin bei der Tat.

Der folgende Code-Chunck bereitet die Daten auf.

Achtung: Wir haben zwar die Information zu den verschiedenen Pinguin-Arten (Penguine_Art) tun aber im Folgenden so, also ob wir diese Information nicht kennen würden. Wir wollen alleine auf Basis der Flossenlängen (Penguine_Flosse) die Gruppenzugehörigkeiten per Clusteranalyse bestimmen. (Im Nachhinein können wir dann mit Hilfe der Daten in Penguine_Art prüfen, wie gut unsere Clusteranalyse ist.)

library("palmerpenguins") # Pinguin-Daten
library("RColorBrewer")   # Hübsche Farben
library("scales")         # Für transparente Farben: alpha()

col_v <- RColorBrewer::brewer.pal(n = 3, name = "Set2")

## Vorbereitung der Daten:
Pinguine <- palmerpenguins::penguins %>% # Pinguin-Daten
  tidyr::as_tibble() %>%                 # Datenformat: 'tibble'-dataframe
  dplyr::filter(species!="Adelie") %>%   # Pinguin-Art 'Adelie' löschen 
  droplevels() %>%                       # Lösche das nicht mehr benötigte Adelie-Level
  tidyr::drop_na() %>%                   # NAs löschen
  dplyr::mutate(Art    = species,        # Variablen umbenennen
                Flosse = flipper_length_mm) %>% 
  dplyr::select(Art, Flosse)             # Variablen auswählen

##  
n      <- nrow(Pinguine)                 # Stichprobenumfang

## Variable 'Penguine_Art' aus Pinguine-Daten "herausziehen"
Penguine_Art    <- dplyr::pull(Pinguine, Art)

## Variable 'Penguine_Flosse' aus Pinguine-Daten "herausziehen"
Penguine_Flosse <- dplyr::pull(Pinguine, Flosse)

## Plot
## Histogramm:
hist(x = Penguine_Flosse, freq = FALSE, 
     xlab="Flossenlänge (mm)", main="Pinguine\n(Zwei Gruppen)",
     col=gray(.65,.5), border=gray(.35,.5), ylim=c(0.0003, 0.039))
## Stipchart hinzufügen:
stripchart(x = Penguine_Flosse, method = "jitter", 
           jitter = .0005, at = .001,
           pch = 21, col=alpha(col_v[3],.5), 
           bg=alpha(col_v[3],.5), cex=1.3, add = TRUE)

Das Clusterverfahren basierend auf Gaußschen Mischverteilungen:

  1. Gaußsche Mischverteilung (per EM Algorithmus) schätzen
  2. Die Datenpunkte \(x_i\) derjenigen Gruppe zuordnen, welche die „a-posteriori-Wahrscheinlichkeit“ maximiert (siehe Abbildung 1.2)


Clusteranalyse basierend auf einer Mischverteilung mit zwei gewichteten Normalverteilungen.

Abbildung 1.2: Clusteranalyse basierend auf einer Mischverteilung mit zwei gewichteten Normalverteilungen.

Abbildung 1.2 zeigt das Resultat einer Clusteranalyse basierend auf einer Mischverteilung zweier gewichteter Normalverteilungen. Cluster-Ergebnis: 95% der Pinguine konnten richtig zugeordnet werden - lediglich auf Basis ihrer Flossenlängen.




Mit Hilfe der folgenden R-Codes kann die obige Clusteranalyse und die Ergebnisgrafik repliziert werden: 

## mclust R-Paket:
## Clusteranalyse mit Hilfe von Gaußschen Mischmodellen
suppressMessages(library("mclust"))

## Anzahl der Gruppen
G <- 2 

## Schätzung des Gaußschen Mischmodells (per EM Algorithmus)
## und Clusteranalyse
mclust_obj <- mclust::Mclust(data = Penguine_Flosse, G=G, 
                              modelNames = "V", 
                              verbose = FALSE)

# summary(mclust_obj)
# str(mclust_obj)

## Geschätzte Gruppen-Zuordnungen (Cluster-Resultat)
class <- mclust_obj$classification

## Anteil der korrekten Zuordnungen:
# cbind(class, Penguine_Art)
round(sum(class == as.numeric(Penguine_Art))/n, 2)

## Geschätzte Mittelwerte 
mean_m <- t(mclust_obj$parameters$mean)

## Geschätzte Varianzen (und evtl. Kovarianzen) 
cov_l  <- list("Cov1" = mclust_obj$parameters$variance$sigmasq[1], 
               "Cov2" = mclust_obj$parameters$variance$sigmasq[2])

## Geschätzte Gewichte (a-priori-Wahrscheinlichkeiten) 
prop_v <- mclust_obj$parameters$pro

## Auswerten der Gaußsche Mischung-Dichtefunktion
np      <- 100 # Anzahl der Auswertungspunkte
xxd     <- seq(min(Penguine_Flosse)-3, max(Penguine_Flosse)+5, length.out = np)
## Mischungs-Dichte
yyd     <- dnorm(xxd, mean_m[1], sqrt(cov_l[[1]]))*prop_v[1] +
           dnorm(xxd, mean_m[2], sqrt(cov_l[[2]]))*prop_v[2]
## Einzel-Dichten
yyd1    <- dnorm(xxd, mean_m[1], sqrt(cov_l[[1]]))*prop_v[1]
yyd2    <- dnorm(xxd, mean_m[2], sqrt(cov_l[[2]]))*prop_v[2]

## Plot
hist(x = Penguine_Flosse, xlab="Flossenlänge (mm)", main="Pinguine\n(Zwei Gruppen)",
     col=gray(.65,.5), border=gray(.35,.5), freq = FALSE, ylim=c(0, 0.04))
lines(x = xxd, y=yyd, lwd=2, col=gray(.35,.75))
lines(x = xxd, y=yyd1, lwd=2, col=gray(.35,.75), lty=2)
lines(x = xxd, y=yyd2, lwd=2, col=gray(.35,.75), lty=2)
stripchart(Penguine_Flosse[class==1], method = "jitter", jitter = .0005, at = .001,
           pch = 21, col=alpha(col_v[1],.5), bg=alpha(col_v[1],.5), cex=1.3, add = TRUE)
stripchart(Penguine_Flosse[class==2], method = "jitter", jitter = .0005, at = .001,
           pch = 21, col=alpha(col_v[2],.5), bg=alpha(col_v[2],.5), cex=1.3, add = TRUE)

Literatur

Horst, Allison Marie, Alison Presmanes Hill, und Kristen B Gorman. 2020. palmerpenguins: Palmer Archipelago (Antarctica) Penguin Data. https://allisonhorst.github.io/palmerpenguins/.