---
lesson: Statistik6
knitr:
  opts_chunk: 
    collapse: false
---

# Statistik 8: Demo

[Demoscript herunterladen (.R)](../purl/Statistik8_Demo.R){.dld}

[Demoscript herunterladen (.qmd)](../purl/Statistik8_Demo.qmd){.dld}

-   Datensatz *Doubs_spe* von  @borcard2011



## Cluster-Analyse k-means

### Datenbeschreibung

Der Datensatz enthält Daten zum Vorkommen von Fischarten und den zugehörigen 
Umweltvariablen im Fluss Doubs (Jura). Es gibt 29 Probestellen (sites), an denen 
jeweils die Abundanzen von 27 Fischarten (auf einer Skalen von 0 bis 5) erhoben 
wurden. In dieser Demo verwenden wir die Artdaten.

```{r}
#| message: false
library("readr")
library("ggplot2")
library("tibble")
```


```{r}
#| message: false
# Daten importieren
spe <- read_delim("./datasets/stat/Doubs_species.csv", delim = ";") |>
  column_to_rownames(var = "Site")

str(spe)
summary(spe)
```

k-means ist eine lineare Methode und daher nicht für Artdaten geeignet .Darum müssen wir unsere Daten transformieren (für die meisten anderen Daten ist die Funktion „scale“, welche jede Variable so skaliert, dass sie einen Mittelwert von 0 und einen Standardabweichungswert von 1 hat), besser geeignet.


```{r}
library("vegan")
spe_norm <- decostand(spe, "normalize")
```

### k-means clustering mit Artdaten

```{r}
# k-means-clustering mit 4 Gruppen durchführen
set.seed(123)
kmeans_1 <- kmeans(spe_norm, centers = 4, nstart = 100)
kmeans_1$cluster

#  Visualisierung
library("factoextra")

fviz_cluster(kmeans_1, main = "", data = spe) +
  theme_classic()
```

Wie viele Cluster (Gruppen) sollen definiert werden? Oft haben wir eine Vorstellung über den 
Range der Anzahl Cluster. Value criterions wie der Simple Structure Index (ssi) können
eine zusätzliche Hilfe sein, um die geeignete Anzahl von Clustern zu finden.

```{r}
# k-means Partionierung, 2 bis 10 Gruppen
set.seed(123)
km_cascade <- cascadeKM(spe_norm,  inf.gr = 2, sup.gr = 10, iter = 100, criterion = "ssi")
km_cascade$results
km_cascade$partition

# Visualisierung citerion Simple Structure Index
plot(km_cascade, sortg = TRUE)

# k-means-Clustering mit 3 Gruppen durchführen
set.seed(123)
kmeans_2 <- kmeans(spe_norm, centers = 3, nstart = 100)


#  Clustering-Resultat in Ordinationsplots darstellen
fviz_cluster(kmeans_2, data = spe) +
  theme_classic() +
  theme(plot.title = element_blank())

# Resultat intepretieren
kmeans_2

# Zuordnung Sites zu den Clustern (separat)
kmeans_2$cluster

# Anzahl Sites pro Cluster
kmeans_2$size

# Mittlere Abundance für jede Variable (Fischart) in jedem Cluster (mit untransformierten Daten)
aggregate(spe, by = list(cluster = kmeans_2$cluster), mean)

# Mittlere Fisch-Artenzahl in jedem Cluster
aggregate( specnumber(spe), by = list(cluster = kmeans_2$cluster), mean)


# Unterschiede Mittlere Fisch-Artenzahl pro Cluster testen

# File für Anova erstellen
spe_2 <- data.frame(spe, 
                    "cluster" = as.factor(kmeans_2$cluster), 
                    "species_richness" = specnumber(spe))
str(spe_2)

aov_1 <- aov(species_richness~cluster, data = spe_2)
summary(aov_1)

library("agricolae")
( Tukey <- HSD.test(aov_1, "cluster") )

sig_letters <- Tukey$groups[order(row.names(Tukey$groups)), ]

ggplot(spe_2, aes(x = cluster,  y = species_richness)) + 
  geom_boxplot() +
  geom_text(data = sig_letters, 
            aes(label = groups, x = c(1:3), 
                y = max(spe_2$species_richness) * 1.2)) +
  theme_classic()
```
