2 Cluster Peces

Análisis multivariante en R: aplicación en ecología

Rosana Ferrero

5 de febrero de 2014

Índice

1. Análisis cluster con R: datos de peces (ambos métodos). 21.1. Análisis jerárquico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.1.1. PASO 1. Elección de variables (o casos) . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.1.2. PASO 2. Elección de la técnica cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.1.3. PASO 3: Determinar el número de conglomerados finales . . . . . . . . . . . . . . . . 51.1.4. PASO 4. Representación e interpretación de los resultados . . . . . . . . . . . . . . . . 8

1.2. Análisis no jerárquico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91.2.1. Partición por k-medias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91.2.2. k-medoides (PAM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

1.3. PLUS: Análisis extra para mejorar la interpretación y validación. . . . . . . . . . . . . . . . . 171.3.1. Comparación con datos ambientales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171.3.2. Ensables de especies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

1

1. Análisis cluster con R: datos de peces (ambos métodos).

Análisis cluster jerárquico y no jerárquico para datos de abundancia de peces a lo largo del río Doubs(francés-suizo) obtenidos por Verneaux (1973) y las características ambientales del sitio. Se quiere analizar:a) Si existen distintos ensambles de peces a lo largo del cauce del río. b) Cómo los cambios ambientalesregistrados en el río se relacionan con los cambios en la estructura de los ensambles de peces. Tenemos30 filas (sitios), 27 especies de peces (columnas del archivo spe), 11 variables ambientales (columnas delarchivo env) y las coordenadas espaciales de los sitios (columnas del archivo spa). Para obtener más detallesde los datos puedes preguntar en R sobre los datos doubs del paquete ade4.

1.1. Análisis jerárquico

# --------- Análisis jerárquico ---------

# cargamos los paquetes que vamos a necesitar (activar vegan luego de# ade4 para evitar conflictos)library(ade4)library(vegan)library(gclus)library(cluster)library(RColorBrewer)library(labdsv)library(mvpart)library(MVPARTwrap)

# cargamos algunas funciones que hemos de necesitar (los archivos# deben estar en nuestro directorio de trabajo)source("hcoplot.R")source("test.a.R")source("coldiss.R")

1.1.1. PASO 1. Elección de variables (o casos)

# PASO 1. Elección de variables (o casos) importamos los archivos de# datos (csv) (los archivos deben estar en nuestro directorio de# trabajo)spe <- read.csv("DoubsSpe.csv", row.names = 1)env <- read.csv("DoubsEnv.csv", row.names = 1)spa <- read.csv("DoubsSpa.csv", row.names = 1)# Remove empty site 8spe <- spe[-8, ]env <- env[-8, ]spa <- spa[-8, ]

# estandarizamos los datosspe.norm <- decostand(spe, "normalize")

1.1.2. PASO 2. Elección de la técnica cluster

2

# PASO 2. Elección de la técnica cluster I) Análisis de conglomerados# jerárquico de los datos de abundancias de especies ---

# ------ Determinar qué método de conglomerado es el mejor ----

# calculamos la matriz de distanciasspe.ch <- vegdist(spe.norm, "euc")

# calculamos el cluster por varios métodos para elegir luego el más# adecuadopar(mfrow = c(2, 3)) #para disponer juntos los gráficos

# Método ’single’spe.ch.single <- hclust(spe.ch, method = "single")plot(spe.ch.single) # graficamos el dendrograma

# Método ’complete’spe.ch.complete <- hclust(spe.ch, method = "complete")plot(spe.ch.complete)

# Métod ’UPGMA’spe.ch.UPGMA <- hclust(spe.ch, method = "average")plot(spe.ch.UPGMA)

# Método ’centroid’spe.ch.centroid <- hclust(spe.ch, method = "centroid")plot(spe.ch.centroid)

# Método de mínima varianza de Wardspe.ch.ward <- hclust(spe.ch, method = "ward")plot(spe.ch.ward)# transformación de la raíz cuadradaspe.ch.ward$height <- sqrt(spe.ch.ward$height)plot(spe.ch.ward)

dev.off() #cierro los gráficos

## null device## 1

# Correlaciones cofenéticas -------

# Singlespe.ch.single.coph <- cophenetic(spe.ch.single)cor(spe.ch, spe.ch.single.coph)

## [1] 0.5992

# Completespe.ch.comp.coph <- cophenetic(spe.ch.complete)cor(spe.ch, spe.ch.comp.coph)

3

## [1] 0.7656

# Averagespe.ch.UPGMA.coph <- cophenetic(spe.ch.UPGMA)cor(spe.ch, spe.ch.UPGMA.coph)

## [1] 0.8608

# Wardspe.ch.ward.coph <- cophenetic(spe.ch.ward)cor(spe.ch, spe.ch.ward.coph)

## [1] 0.7985

cor(spe.ch, spe.ch.ward.coph, method = "spearman")

## [1] 0.7661

# el mayor coeficiente de correlación cofenética corresponde al# método UPGMA

# Diagrama de Shepard ---- Para ilustrar las relaciones entre la# matriz de distancias original y la matriz de distancias cofenéticapar(mfrow = c(2, 2))plot(spe.ch, spe.ch.single.coph, xlab = "Chord distance", ylab = "Cophenetic distance",

asp = 1, xlim = c(0, sqrt(2)), ylim = c(0, sqrt(2)), main = c("Single linkage",paste("Cophenetic correlation =", round(cor(spe.ch, spe.ch.single.coph),

3))))abline(0, 1)lines(lowess(spe.ch, spe.ch.single.coph), col = "red")

plot(spe.ch, spe.ch.comp.coph, xlab = "Chord distance", ylab = "Cophenetic distance",asp = 1, xlim = c(0, sqrt(2)), ylim = c(0, sqrt(2)), main = c("Complete linkage",

paste("Cophenetic correlation =", round(cor(spe.ch, spe.ch.comp.coph),3))))

abline(0, 1)lines(lowess(spe.ch, spe.ch.comp.coph), col = "red")

plot(spe.ch, spe.ch.UPGMA.coph, xlab = "Chord distance", ylab = "Cophenetic distance",asp = 1, xlim = c(0, sqrt(2)), ylim = c(0, sqrt(2)), main = c("UPGMA",

paste("Cophenetic correlation =", round(cor(spe.ch, spe.ch.UPGMA.coph),3))))

abline(0, 1)lines(lowess(spe.ch, spe.ch.UPGMA.coph), col = "red")

plot(spe.ch, spe.ch.ward.coph, xlab = "Chord distance", ylab = "Cophenetic distance",asp = 1, xlim = c(0, sqrt(2)), ylim = c(0, max(spe.ch.ward$height)),main = c("Ward clustering", paste("Cophenetic correlation =", round(cor(spe.ch,

spe.ch.ward.coph), 3))))abline(0, 1)lines(lowess(spe.ch, spe.ch.ward.coph), col = "red")

# Distancia de Gower (1983) ------(gow.dist.single <- sum((spe.ch - spe.ch.single.coph)^2))

## [1] 95.41

4

(gow.dist.comp <- sum((spe.ch - spe.ch.comp.coph)^2))

## [1] 40.49

(gow.dist.UPGMA <- sum((spe.ch - spe.ch.UPGMA.coph)^2))

## [1] 11.67

(gow.dist.ward <- sum((spe.ch - spe.ch.ward.coph)^2))

## [1] 532

# según la distancia de Gower, el mejor método es UPGMA

# por tanto, seleccionaríamos los resultados del método UPGMA sin# embargo, para comparar el método con análisis posteriores, aquí# vamos a seleccionar el método Ward

95

119

3020

21 2229

26 27 2817 18

1615 10

1113 14

122 3 7

4 625

23 24

0.2

0.4

0.6

0.8

Cluster Dendrogram

hclust (*, "single")spe.ch

Hei

ght 1

122 3 7

1113 14

915 16

104 6

519

17 18 29 3020

21 2226

27 2825

23 240.2

0.6

1.0

1.4

Cluster Dendrogram

hclust (*, "complete")spe.ch

Hei

ght

20 3026 27 2829

21 2225

23 241 9

1011

13 144 6 12

2 3 75

15 16 1917 18

0.2

0.4

0.6

0.8

1.0

1.2

Cluster Dendrogram

hclust (*, "average")spe.ch

Hei

ght

19

520 30

26 29 21 22 27 2819

17 1825

23 2416

15 1011 13 14

6 412

23 7

0.1

0.3

0.5

0.7

Cluster Dendrogram

hclust (*, "centroid")spe.ch

Hei

ght

20 21 22 26 27 28 29 3025 23 24

5 915 16 19 17 18 2 3 7 10 4 6

113 14 11 12

02

46

Cluster Dendrogram

hclust (*, "ward")spe.ch

Hei

ght

2021 22 26 27 28 29 30

2523 24

5 915 16 19

17 18 2 3 710 4 6

113 14 11 12

0.0

0.5

1.0

1.5

2.0

2.5

Cluster Dendrogram

hclust (*, "ward")spe.ch

Hei

ght

1.1.3. PASO 3: Determinar el número de conglomerados finales

# PASO 3: Determinar el número de conglomerados finales Podemos# seleccionarlo de manera subjetiva mediante el examen visual del# dendrograma, o podemos elegir un criterio más objetivo para# predeterminar el número de cluster a retener.

# ------# Seleccionar el número de conglomerados final ------#

# Gráfico de nivel de fusión ------------ Warddev.off()

5

## null device## 1

plot(spe.ch.ward$height, nrow(spe):2, type = "S", main = "Fusion levels - Chord - Ward",ylab = "k (number of clusters)", xlab = "h (node height)", col = "grey")

text(spe.ch.ward$height, nrow(spe):2, nrow(spe):2, col = "red", cex = 0.8)

# Número óptimo de cluster según el ancho de silueta --------# (Rousseeuw quality index)

# graficamos el ancho medio de silueta para todas las particiones,# excepto para la trivial (k=1) creamos un vector nulo para ir# agregando los valoresasw <- numeric(nrow(spe))# creamos un loop para los cálculosfor (k in 2:(nrow(spe) - 1)) {

sil <- silhouette(cutree(spe.ch.ward, k = k), spe.ch)asw[k] <- summary(sil)$avg.width

}k.best <- which.max(asw)

# graficamos

plot(1:nrow(spe), asw, type = "h", main = "Silhouette-optimal number of clusters, Ward",xlab = "k (number of groups)", ylab = "Average silhouette width")

axis(1, k.best, paste("optimum", k.best, sep = "\n"), col = "red", font = 2,col.axis = "red")

points(k.best, max(asw), pch = 16, col = "red", cex = 1.5)cat("", "Silhouette-optimal number of clusters k =", k.best, "\n", "with an average silhouette width of",

max(asw), "\n")

## Silhouette-optimal number of clusters k = 2## with an average silhouette width of 0.3658

# según este criterio se deben seleccionar 2 clusters (k=2), sin# embargo, por razones ecológicas, resulta más interesante formar 4# particiones

# Estadístico de Mantel (Pearson) -----------------------

# Creamos una función para calcular la matriz de distancia binaria# con los gruposgrpdist <- function(X) {

require(cluster)gr <- as.data.frame(as.factor(X))distgr <- daisy(gr, "gower")distgr

}

# calculamos el estadísticos para el método de Wardkt <- data.frame(k = 1:nrow(spe), r = 0)for (i in 2:(nrow(spe) - 1)) {

gr <- cutree(spe.ch.ward, i)distgr <- grpdist(gr)mt <- cor(spe.ch, distgr, method = "pearson")

6

kt[i, 2] <- mt}kt

## k r## 1 1 0.0000## 2 2 0.6554## 3 3 0.7081## 4 4 0.7155## 5 5 0.6429## 6 6 0.6741## 7 7 0.6865## 8 8 0.6879## 9 9 0.6926## 10 10 0.6499## 11 11 0.6464## 12 12 0.6396## 13 13 0.6304## 14 14 0.5009## 15 15 0.4960## 16 16 0.4375## 17 17 0.4027## 18 18 0.3950## 19 19 0.3852## 20 20 0.3571## 21 21 0.3424## 22 22 0.3271## 23 23 0.3111## 24 24 0.2931## 25 25 0.2512## 26 26 0.1971## 27 27 0.1615## 28 28 0.1149## 29 29 0.0000

k.best <- which.max(kt$r)# graficamosplot(kt$k, kt$r, type = "h", main = "Mantel-optimal number of clusters - Ward",

xlab = "k (number of groups)", ylab = "Pearson’s correlation")axis(1, k.best, paste("optimum", k.best, sep = "\n"), col = "red", font = 2,

col.axis = "red")points(k.best, max(kt$r), pch = 16, col = "red", cex = 1.5)cat("", "Mantel-optimal number of clusters k =", k.best, "\n", "with a matrix linear correlation of",

max(kt$r), "\n")

## Mantel-optimal number of clusters k = 4## with a matrix linear correlation of 0.7155

# este criterio selecciona como óptimo el formar 4 conglomerados o# clusters (k=4).

Validación final

7

# Validación final Gráfico de silueta para la partición final# -------- Vemos si los miembros de los grupos están bien# clasificados

# para 4 clustersk <- 4# graficamos la siluetacutg <- cutree(spe.ch.ward, k = k)sil <- silhouette(cutg, spe.ch)rownames(sil) <- row.names(spe)

plot(sil, main = "Silhouette plot - Chord - Ward", cex.names = 0.8, col = 2:(k +1), nmax = 100)

# los grupos 1 y 3 son bastante coherentes, mientras que el grupo 2# contiene algunos miembros mal clasificados

25242320263028292722211519

518

9161710

61

144

1311

712

32

Silhouette width si

0.0 0.2 0.4 0.6 0.8 1.0

Silhouette plot − Chord − Ward

Average silhouette width : 0.35

n = 29 4 clusters Cj

j : nj | avei∈Cj si

1 : 11 | 0.37

2 : 7 | 0.04

3 : 8 | 0.60

4 : 3 | 0.34

1.1.4. PASO 4. Representación e interpretación de los resultados

# PASO 4. Representación e interpretación de los resultados Agregar# colores a los grupos hcoplot(spe.ch.ward, spe.ch, k=4)

# Gráfico con mapa del río Doubs (para datos espacialmente# explícitos) -----plot(spa, asp = 1, type = "n", main = "Four Ward clusters along the Doubs river",

xlab = "x coordinate (km)", ylab = "y coordinate (km)")lines(spa, col = "light blue")text(50, 10, "Upstream", cex = 1.2)

8

text(25, 115, "Downstream", cex = 1.2)

# Add the four groupsspebc.ward.g <- cutree(spe.ch.ward, k = 4)grw <- spebc.ward.gk <- length(levels(factor(grw)))for (i in 1:k) {

points(spa[grw == i, 1], spa[grw == i, 2], pch = i + 20, cex = 3, col = i +1, bg = i + 1)

}text(spa, row.names(spa), cex = 0.8, col = "white", font = 2)legend("bottomright", paste("Cluster", 1:k), pch = (1:k) + 20, col = 2:(k +

1), pt.bg = 2:(k + 1), pt.cex = 2, bty = "n")

−100 0 100 200 300

050

100

150

200

Four Ward clusters along the Doubs river

x coordinate (km)

y co

ordi

nate

(km

)

Upstream

Downstream

12 3

45

67

910

11

12

13

1415

161718

192021

22

2324

25

26

2728

29

30

Cluster 1Cluster 2Cluster 3Cluster 4

1.2. Análisis no jerárquico

1.2.1. Partición por k-medias

# Partición por k-medias El criterio de partición de k-medias es# similar al del método de Ward, ya que el método iterativamente# minimiza el total del error de la suma de cuadrados. Es un método# lineal, por loque no es apropiado para datos de abundancias de# especies con muchos ceros.

# con 4 gruposset.seed(1)spe.kmeans <- kmeans(spe.norm, centers = 4, nstart = 100)

9

# Comparación con la clasificación de Ward con 4 grupostable(spe.kmeans$cluster, spebc.ward.g)

## spebc.ward.g## 1 2 3 4## 1 0 0 0 3## 2 11 1 0 0## 3 0 6 0 0## 4 0 0 8 0

# vemos que solo se diferencian en la clasificación de 1 caso.

# Cuál es el mejor número de clusters a retener? ---- calculamos# particiones desde 2 a 10 clustersspe.KM.cascade <- cascadeKM(spe.norm, inf.gr = 2, sup.gr = 10, iter = 100,

criterion = "ssi")summary(spe.KM.cascade)

## Length Class Mode## partition 261 -none- numeric## results 18 -none- numeric## criterion 1 -none- character## size 90 -none- numeric

spe.KM.cascade$results

## 2 groups 3 groups 4 groups 5 groups 6 groups 7 groups 8 groups## SSE 8.2149 6.4768 5.0720 4.302 3.5856 2.9524 2.48405## ssi 0.1312 0.1676 0.1245 0.115 0.1282 0.1306 0.07276## 9 groups 10 groups## SSE 2.0522 1.75993## ssi 0.1267 0.07095

# el mínimo del valor SEE es el criterio para encontrar la agrupación# óptima de los objetos para un valor dado de k mientras que el# máximo del índice SSI encuentra el valor k óptimo.

# vemos que el criterio SSI se maximiza para 3 grupos.

# vemos la partición creada y la graficamoshead(spe.KM.cascade$partition)

## 2 groups 3 groups 4 groups 5 groups 6 groups 7 groups 8 groups## 1 2 3 4 4 4 5 4## 2 2 3 4 3 2 6 8## 3 2 3 4 3 2 6 8## 4 2 3 4 3 2 6 8## 5 1 1 3 2 5 2 3## 6 2 3 4 3 2 6 8## 9 groups 10 groups## 1 7 1## 2 2 10## 3 2 10## 4 2 10## 5 9 9## 6 2 4

10

plot(spe.KM.cascade, sortg = TRUE)

# Evaluamos los contenidos de los cluster, en este caso para 4# clusters. reordenamos según los resultados del clusterhead(spe[order(spe.kmeans$cluster), ])

## CHA TRU VAI LOC OMB BLA HOT TOX VAN CHE BAR SPI GOU BRO PER BOU## 23 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0## 24 0 0 0 0 0 0 1 0 0 2 0 0 1 0 0 0## 25 0 0 0 0 0 0 0 0 1 1 0 0 2 1 0 0## 1 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0## 2 0 5 4 3 0 0 0 0 0 0 0 0 0 0 0 0## 3 0 5 5 5 0 0 0 0 0 0 0 0 0 1 0 0## PSO ROT CAR TAN BCO PCH GRE GAR BBO ABL ANG## 23 0 0 0 0 0 0 0 1 0 2 0## 24 1 0 0 0 0 0 2 2 1 5 0## 25 0 1 0 0 0 0 1 1 0 3 0## 1 0 0 0 0 0 0 0 0 0 0 0## 2 0 0 0 0 0 0 0 0 0 0 0## 3 0 0 0 0 0 0 0 0 0 0 0

# reordenamos sitios y especies con la función vegemiteord.KM <- vegemite(spe, spe.kmeans$cluster)

#### 222 111111 111122222223## 34512346701234559678901267890## TRU ...3554351355542.321.......1.## CHA ..........12233..211.........## OMB ..........12342...11.......1.## VAI ....4554444455431343311....1.## LOC ....355554142452354352111111.## BLA ............234..5211......1.## CHE 121...121211.1325223132423443## VAN ..1....112....35.532222312233## ABL 253...............22355555555## GOU .12...11.1...122.212445534455## TOX .................4433222.1122## BRO ..1..121.....1.4.111123324354## GAR 121....1.......54122555545555## PER ......21.......4.123123411245## HOT .1................11222311121## TAN ......12......131111244435435## BAR .............12..233244524353## SPI .................143232111135## GRE .21................1123445555## ROT ..1............2....122211225## PSO .1...............111122323453## CAR .................111112312435## BOU ..................12233323455## BBO .1..................123424545## ANG ..................11122223445## BCO ....................113423345

11

## PCH ......................1212345## sites species## 29 27

head(spe[ord.KM$sites, ord.KM$species])

## TRU CHA OMB VAI LOC BLA CHE VAN ABL GOU TOX BRO GAR PER HOT TAN## 23 0 0 0 0 0 0 1 0 2 0 0 0 1 0 0 0## 24 0 0 0 0 0 0 2 0 5 1 0 0 2 0 1 0## 25 0 0 0 0 0 0 1 1 3 2 0 1 1 0 0 0## 1 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0## 2 5 0 0 4 3 0 0 0 0 0 0 0 0 0 0 0## 3 5 0 0 5 5 0 0 0 0 0 0 1 0 0 0 0## BAR SPI GRE ROT PSO CAR BOU BBO ANG BCO PCH## 23 0 0 0 0 0 0 0 0 0 0 0## 24 0 0 2 0 1 0 0 1 0 0 0## 25 0 0 1 1 0 0 0 0 0 0 0## 1 0 0 0 0 0 0 0 0 0 0 0## 2 0 0 0 0 0 0 0 0 0 0 0## 3 0 0 0 0 0 0 0 0 0 0 0

# los valores centrales (medios) de cada cluster, para cada variable,# son:head(spe.kmeans$centers[, 1:5])

## CHA TRU VAI LOC OMB## 1 0.00000 0.000000 0.00000 0.00000 0.000000## 2 0.10380 0.542301 0.50087 0.43326 0.114024## 3 0.06168 0.122088 0.26994 0.35943 0.032665## 4 0.00000 0.006691 0.02506 0.06987 0.006691

# si usamos clustIndex (para 2 culsters)cl <- cclust(as.matrix(spe.norm), 2, 100, verbose = TRUE, method = "kmeans")resultindexes <- clustIndex(cl, as.matrix(spe.norm), index = "all")resultindexes

12

5 10 15 20 25

K−means partitions comparison

Objects

Num

ber

of g

roup

s in

eac

h pa

rtiti

on

24

68

10

5 10 15 20 25 0.08 0.12 0.16

ssicriterion

Values

24

68

10

0.08 0.12 0.16

1.2.2. k-medoides (PAM)

# ------# Partición alrededor de medoides (PAM) ------- ------# El# objetivo es encontrar k objetos representativos que minimicen la# suma de las disimilaridades de las observaciones respecto al objeto# representativo más cercano. A diferencia del método de k-medias,# no minimiza la suma de los cuadrados de las distancias euclídeas# dentro de los grupos (método de cuadrados mínimos)

# Para determinar el número de clusters óptimo ------- Generamos un# loop para calcular los anchos de silueta media desde 2 a 28# clusters

asw <- numeric(nrow(spe))for (k in 2:(nrow(spe) - 1)) asw[k] <- pam(spe.ch, k, diss = TRUE)$silinfo$avg.widthk.best <- which.max(asw)cat("", "Silhouette-optimal number of clusters k =", k.best, "\n", "with an average silhouette width of",

max(asw), "\n")

## Silhouette-optimal number of clusters k = 2## with an average silhouette width of 0.3841

# graficamosplot(1:nrow(spe), asw, type = "h", main = "Choice of the number of clusters",

xlab = "k (number of clusters)", ylab = "Average silhouette width")axis(1, k.best, paste("optimum", k.best, sep = "\n"), col = "red", font = 2,

col.axis = "red")points(k.best, max(asw), pch = 16, col = "red", cex = 1.5)# la opción óptima para PAM es k=2 con asw=0.3841

# Si aún así elegimos k=4 clusters para comparar con los anteriores# métodos

13

spe.ch.pam <- pam(spe.ch, k = 4, diss = TRUE)summary(spe.ch.pam)

## Medoids:## ID## [1,] "11" "12"## [2,] "6" "6"## [3,] "16" "17"## [4,] "25" "26"## Clustering vector:## 1 2 3 4 5 6 7 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24## 1 1 1 2 2 2 1 2 2 1 1 1 1 1 3 3 3 3 4 4 4 4 4## 25 26 27 28 29 30## 4 4 4 4 4 4## Objective function:## build swap## 0.4426 0.4391#### Numerical information per cluster:## size max_diss av_diss diameter separation## [1,] 9 0.8770 0.4351 1.1169 0.3429## [2,] 5 0.8815 0.4962 1.0855 0.3429## [3,] 4 0.5989 0.3679 0.8037 0.4977## [4,] 11 0.8842 0.4422 1.1276 0.4977#### Isolated clusters:## L-clusters: character(0)## L*-clusters: character(0)#### Silhouette plot information:## cluster neighbor sil_width## 13 1 2 0.42989## 12 1 2 0.40576## 2 1 2 0.38341## 14 1 2 0.35676## 11 1 2 0.34189## 3 1 2 0.27044## 1 1 2 0.24590## 7 1 2 0.24136## 15 1 3 0.05062## 6 2 1 0.09344## 10 2 1 0.06179## 9 2 3 0.04917## 5 2 3 -0.03817## 4 2 1 -0.14709## 17 3 2 0.45358## 18 3 4 0.38401## 16 3 2 0.23661## 19 3 4 0.15365## 26 4 3 0.48031## 27 4 3 0.46381## 28 4 3 0.45850## 22 4 3 0.41331## 21 4 3 0.39376

14

## 29 4 3 0.36427## 30 4 3 0.36421## 24 4 3 0.31070## 25 4 3 0.26599## 20 4 3 0.23650## 23 4 3 0.21030## Average silhouette width per cluster:## [1] 0.302892 0.003827 0.306963 0.360150## Average silhouette width of total data set:## [1] 0.2736#### Available components:## [1] "medoids" "id.med" "clustering" "objective" "isolation"## [6] "clusinfo" "silinfo" "diss" "call"

spe.ch.pam.g <- spe.ch.pam$clusteringspe.ch.pam$silinfo$widths

## cluster neighbor sil_width## 13 1 2 0.42989## 12 1 2 0.40576## 2 1 2 0.38341## 14 1 2 0.35676## 11 1 2 0.34189## 3 1 2 0.27044## 1 1 2 0.24590## 7 1 2 0.24136## 15 1 3 0.05062## 6 2 1 0.09344## 10 2 1 0.06179## 9 2 3 0.04917## 5 2 3 -0.03817## 4 2 1 -0.14709## 17 3 2 0.45358## 18 3 4 0.38401## 16 3 2 0.23661## 19 3 4 0.15365## 26 4 3 0.48031## 27 4 3 0.46381## 28 4 3 0.45850## 22 4 3 0.41331## 21 4 3 0.39376## 29 4 3 0.36427## 30 4 3 0.36421## 24 4 3 0.31070## 25 4 3 0.26599## 20 4 3 0.23650## 23 4 3 0.21030

# Comparamos con las anteriores clasificacionestable(spe.ch.pam.g, spebc.ward.g)

## spebc.ward.g## spe.ch.pam.g 1 2 3 4

15

## 1 8 1 0 0## 2 3 2 0 0## 3 0 4 0 0## 4 0 0 8 3

table(spe.ch.pam.g, spe.kmeans$cluster)

#### spe.ch.pam.g 1 2 3 4## 1 0 9 0 0## 2 0 3 2 0## 3 0 0 4 0## 4 3 0 0 8

# los resultados de PAM difieren mucho de los métodos anteriores

# Qué solución (PAM o k-means) tiene mejor silhouette? -----par(mfrow = c(1, 2))k <- 4sil <- silhouette(spe.kmeans$cluster, spe.ch)rownames(sil) <- row.names(spe)plot(sil, main = "Silhouette plot - k-means", cex.names = 0.8, col = 2:(k +

1))plot(silhouette(spe.ch.pam), main = "Silhouette plot - PAM", cex.names = 0.8,

col = 2:(k + 1))

# vemos como dos métodos que tienen un mismo objetivo y clase# (no-jerárquico), producen resultados distintos el usuario debe# elegir cuál es la mejor clasificación, según alguna otra# información pertinente o según la interpretación de los datos# ambientales, por ejemplo.

# ----------------# bonus: graficar los cluster de 4 k-means en un# mapa con el río Doubs ---- ----------------#dev.off()

## null device## 1

plot(spa, asp = 1, type = "n", main = "Four k-means groups", xlab = "x coordinate (km)",ylab = "y coordinate (km)")

lines(spa, col = "light blue")text(50, 10, "Upstream", cex = 1.2)text(25, 115, "Downstream", cex = 1.2)grKM <- spe.kmeans$clusterk <- length(levels(factor(grKM)))for (i in 1:k) {

points(spa[grKM == i, 1], spa[grKM == i, 2], pch = i + 20, cex = 3,col = i + 1, bg = i + 1)

}text(spa, row.names(spa), cex = 0.8, col = "white", font = 2)legend("bottomright", paste("Cluster", 1:k), pch = (1:k) + 20, col = 2:(k +

1), pt.bg = 2:(k + 1), pt.cex = 2, bty = "n")

16

0 5 10 15 20 25 30

0.0

0.1

0.2

0.3

Choice of the number of clusters

k (number of clusters)

Ave

rage

silh

ouet

te w

idth

optimum2

2026302928272122195

916181715106

141413117

3122252423

Silhouette width si

0.0 0.2 0.4 0.6 0.8 1.0

Silhouette plot − k−means



j : nj | avei∈Cj si1 : 3 | 0.34

2 : 12 | 0.37

3 : 6 | 0.08

4 : 8 | 0.59232025243029212228272619161817

459106

15713

111421213

Silhouette width si

0.0 0.2 0.4 0.6 0.8 1.0

Silhouette plot − PAM



j : nj | avei∈Cj si

1 : 9 | 0.30

2 : 5 | 0.004

3 : 4 | 0.31

4 : 11 | 0.36

1.3. PLUS: Análisis extra para mejorar la interpretación y validación.

1.3.1. Comparación con datos ambientales.

# necesitaremos cargar los siguiente paquetesrequire(pgirmess)

# comparaciones Post hoc. para instalar el paquete:# https://rforge.net/NCStats/files/; http://www.rrforge.net/FSAdatas/# o intentar:# source(’http://www.rforge.net/NCStats/InstallNCStats.R’)# library(FSAdata)library(NCStats)

17

# -----------------# Comparación con datos ambientales---------# Relación entre los clusters de peces y 4 variables ambientales# basados en los resultados del cluster k-means -----------------#

# Validación externa:

# 1) Comparando la tipología con datos externos (aproximación ANOVA)# ----- Hemos visto que solamente con los datos de especies, tanto# con índices de silueta o calidad del cluster, no hemos podido# seleccionar la mejor partición de los sitios. Podemos confrontar# los resultados del cluster (respuesta) con datos externos (tal vez# generados con análisis discriminantes para crear variables# independientes). También podemos hacer el planteamiento opuesto,# considerar como factores los clusters obtenidos con los datos de# especies, para realizar un ANOVA.

attach(env)

# Boxplots para variables cualitativas ambientales: altitud,# pendiente, oxígeno y amoniopar(mfrow = c(2, 2))boxplot(sqrt(alt) ~ spe.kmeans$cluster, main = "Altitude", las = 1, ylab = "sqrt(alt)",

col = 2:5, varwidth = TRUE)boxplot(log(pen) ~ spe.kmeans$cluster, main = "Slope", las = 1, ylab = "log(pen)",

col = 2:5, varwidth = TRUE)boxplot(oxy ~ spe.kmeans$cluster, main = "Oxygen", las = 1, ylab = "oxy",

col = 2:5, varwidth = TRUE)boxplot(sqrt(amm) ~ spe.kmeans$cluster, main = "Ammonium", las = 1, ylab = "sqrt(amm)",

col = 2:5, varwidth = TRUE)

# Ponemos a prueba los supuestos del ANOVA Residuos normales . Ho:# la variable tiene distribución normalshapiro.test(resid(lm(sqrt(alt) ~ as.factor(spe.kmeans$cluster))))

## Shapiro-Wilk normality test with resid(lm(sqrt(alt) ~ as.factor(spe.kmeans$cluster)))## W = 0.9624, p-value = 0.3756

shapiro.test(resid(lm(log(pen) ~ as.factor(spe.kmeans$cluster))))

## Shapiro-Wilk normality test with resid(lm(log(pen) ~ as.factor(spe.kmeans$cluster)))## W = 0.9518, p-value = 0.204

shapiro.test(resid(lm(oxy ~ as.factor(spe.kmeans$cluster))))

## Shapiro-Wilk normality test with resid(lm(oxy ~ as.factor(spe.kmeans$cluster)))## W = 0.9394, p-value = 0.09674

shapiro.test(resid(lm(sqrt(amm) ~ as.factor(spe.kmeans$cluster))))

## Shapiro-Wilk normality test with resid(lm(sqrt(amm) ~ as.factor(spe.kmeans$cluster)))## W = 0.954, p-value = 0.2319

# no rechazamos normalidad

# Homogeneidad de varianza. Ho: las varianzas son iguales entre# gruposbartlett.test(sqrt(alt), as.factor(spe.kmeans$cluster))

18

## Bartlett test of homogeneity of variances with sqrt(alt) and as.factor(spe.kmeans$cluster)## Bartlett's K-squared = 16.95, df = 3, p-value = 0.0007227

bartlett.test(log(pen), as.factor(spe.kmeans$cluster))

## Bartlett test of homogeneity of variances with log(pen) and as.factor(spe.kmeans$cluster)## Bartlett's K-squared = 2.983, df = 3, p-value = 0.3942

bartlett.test(oxy, as.factor(spe.kmeans$cluster))

## Bartlett test of homogeneity of variances with oxy and as.factor(spe.kmeans$cluster)## Bartlett's K-squared = 1.826, df = 3, p-value = 0.6093

bartlett.test(sqrt(amm), as.factor(spe.kmeans$cluster))

## Bartlett test of homogeneity of variances with sqrt(amm) and as.factor(spe.kmeans$cluster)## Bartlett's K-squared = 5.72, df = 3, p-value = 0.126

# solo rechazamos homogeneidad de varianza para alt

# por tanto, alt lo analizamos con técnicas no-paramétricas, y las# demás con técnicas paramétricas

# Realizamos el ANOVA. Ho: mu_grupo1=mu_grupo2=.....summary(fm1 <- aov(log(pen) ~ as.factor(spe.kmeans$cluster)))

## Df Sum Sq Mean Sq F value Pr(>F)## as.factor(spe.kmeans$cluster) 3 18.1 6.02 7.28 0.0011 **## Residuals 25 20.7 0.83## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fm2 <- aov(oxy ~ as.factor(spe.kmeans$cluster)))

## Df Sum Sq Mean Sq F value Pr(>F)## as.factor(spe.kmeans$cluster) 3 103.5 34.5 26.3 6.8e-08 ***## Residuals 25 32.8 1.3## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(fm3 <- aov(sqrt(amm) ~ as.factor(spe.kmeans$cluster)))

## Df Sum Sq Mean Sq F value Pr(>F)## as.factor(spe.kmeans$cluster) 3 2.613 0.871 49.5 1.2e-10 ***## Residuals 25 0.439 0.018## ---## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# rechazamos igualdad de medias, para todos.

# Kruskal-Wallis test of variable alt. Ho: misma distribución en# todos los grupos.(fm4 <- kruskal.test(alt ~ as.factor(spe.kmeans$cluster)))

## Kruskal-Wallis rank sum test with alt by as.factor(spe.kmeans$cluster)## Kruskal-Wallis chi-squared = 21.53, df = 3, p-value =## 8.186e-05

19

# rechazamos igualdad de medias

# realizamos comparaciones múltiples post-hoc para determinar entre# qué cluster existen diferencias significativas. Ho:# mu_grupo1=mu_grupo2par(mfrow = c(2, 2))

TukeyHSD(fm1, ordered = TRUE)

## Tukey multiple comparisons of means## 95% family-wise confidence level## factor levels have been ordered#### Fit: aov(formula = log(pen) ~ as.factor(spe.kmeans$cluster))#### $às.factor(spe.kmeans$cluster)`## diff lwr upr p adj## 4-1 0.1064 -1.58638 1.799 0.9981## 3-1 0.5881 -1.17995 2.356 0.7971## 2-1 1.8129 0.19896 3.427 0.0235## 3-4 0.4817 -0.86865 1.832 0.7613## 2-4 1.7066 0.56530 2.848 0.0020## 2-3 1.2249 -0.02533 2.475 0.0563

plot(TukeyHSD(fm1))TukeyHSD(fm2, ordered = TRUE)

## Tukey multiple comparisons of means## 95% family-wise confidence level## factor levels have been ordered#### Fit: aov(formula = oxy ~ as.factor(spe.kmeans$cluster))#### $às.factor(spe.kmeans$cluster)`## diff lwr upr p adj## 4-1 3.188 1.0532 5.322 0.0020## 3-1 4.233 2.0042 6.463 0.0001## 2-1 6.083 4.0484 8.118 0.0000## 3-4 1.046 -0.6567 2.748 0.3500## 2-4 2.896 1.4569 4.335 0.0001## 2-3 1.850 0.2737 3.426 0.0171

plot(TukeyHSD(fm2))TukeyHSD(fm3, ordered = TRUE)

## Tukey multiple comparisons of means## 95% family-wise confidence level## factor levels have been ordered#### Fit: aov(formula = sqrt(amm) ~ as.factor(spe.kmeans$cluster))#### $às.factor(spe.kmeans$cluster)`## diff lwr upr p adj## 3-2 0.32984 0.1475 0.5122 0.0002## 4-2 0.36747 0.2010 0.5339 0.0000

20

## 1-2 1.00955 0.7742 1.2449 0.0000## 4-3 0.03763 -0.1593 0.2346 0.9521## 1-3 0.67971 0.4219 0.9376 0.0000## 1-4 0.64208 0.3952 0.8890 0.0000

plot(TukeyHSD(fm3))

kruskalmc(alt ~ as.factor(spe.kmeans$cluster))

## Multiple comparison test after Kruskal-Wallis## p.value: 0.05## Comparisons## obs.dif critical.dif difference## 1-2 15.333 14.50 TRUE## 1-3 9.833 15.88 FALSE## 1-4 1.375 15.21 FALSE## 2-3 5.500 11.23 FALSE## 2-4 16.708 10.25 TRUE## 3-4 11.208 12.13 FALSE

detach(env)

# El procedimiento opuesto. vamos a realizar cluster con las# variables ambientales (para obtener conjuntos de tipos de hábitats)# y probar si las especies responden al tipo de habitat. Utilizaremos# análisis de especies indicadoras, Otra alternativa sería probar la# relación especies-hábitat directamente, ver cap. 6 de Numerical# ecology with R Of course, the reverse procedure could be applied as# well.

# 2) Tablas de continfencia de dos tipologías ------------ Si# queremos comparar la tipología generada con los datos de especies# con aquella que obtenemos con datos ambientales podemos generar una# tabla de contingencia para cruzar ambas tipologías y probar si# existe relación mediante la prueba de independencia Chi-cuadrado.

# Tipología basada en datos ambientalesenv2 <- env[, -1]env.de <- vegdist(scale(env2), "euc")set.seed(1)env.kmeans <- kmeans(env.de, centers = 4, nstart = 100)env.KM.4 <- env.kmeans$cluster

# Tabla de contingencia: columnas=env, filas=sp(tab <- table(as.factor(spe.kmeans$cluster), as.factor(env.kmeans$cluster)))

#### 1 2 3 4## 1 0 0 1 2## 2 7 4 0 1## 3 2 4 0 0## 4 0 3 5 0

# agregamos nombres

21

rownames(tab) <- c("sp_c1", "sp_c2", "sp_c3", "sp_c4")colnames(tab) <- c("env_c1", "env_c2", "env_c3", "env_c4")

# Prueba de independencia Chi-cuadrado. Ho: X e Y son variables# independientes.(M1 <- chisq.test(tab))

## Pearson's Chi-squared test with tab## X-squared = 30.23, df = 9, p-value = 0.0004014

# con test de permutación(M2 <- chisq.test(tab, simulate.p.value = TRUE))

## Pearson's Chi-squared test with simulated p-value## (based on 2000 replicates) with tab## X-squared = 30.23, df = NA, p-value = 0.0004998

# rechazamos la independencia, por lo que existe una relación entre# tipología ambiental y de especies.

# recordar: columnas=env, filas=spchisqPostHoc(chisq.test(t(tab))) # comparación entre columnas (env)

## Adjusted p-values used the fdr method.## comparison raw.p adj.p## 1 env_c1 vs. env_c2 NaN NaN## 2 env_c1 vs. env_c3 0.0018 0.0054## 3 env_c1 vs. env_c4 NaN NaN## 4 env_c2 vs. env_c3 0.0322 0.0322## 5 env_c2 vs. env_c4 0.0262 0.0322## 6 env_c3 vs. env_c4 NaN NaN

chisqPostHoc(M2) # comparación entre filas (sp)

## Adjusted p-values used the fdr method.## comparison raw.p adj.p## 1 sp_c1 vs. sp_c2 0.0127 0.0190## 2 sp_c1 vs. sp_c3 0.0293 0.0293## 3 sp_c1 vs. sp_c4 NaN NaN## 4 sp_c2 vs. sp_c3 NaN NaN## 5 sp_c2 vs. sp_c4 0.0050 0.0149## 6 sp_c3 vs. sp_c4 NaN NaN

# produce NaN cuando hay ceros en los datos Para columnas, vemos que# existen diferencias significativas entre 1-3, 1-4, 2-3. Para# filas, vemos que existen diferencias significativas entre 1-2, 1-3,# 2-4.

# para que sea más fácil la interpretación, graficamoslibrary(plotrix)barp(tab, main = "Clusters", ylab = "Value", names.arg = colnames(tab),

col = 1:4)# agregamos leyendaaddtable2plot(-0.5, 5.5, tab, bty = "o", display.rownames = TRUE, hlines = TRUE,

title = "The table")

22

1 2 3 4

15

20

25

30

Altitude

sqrt

(alt)

1 2 3 4

−1

0

1

2

3

4

Slope

log(

pen)

1 2 3 4

4

6

8

10

12

Oxygen

oxy

1 2 3 4

0.0

0.2

0.4

0.6

0.8

1.0

1.2

Ammonium

sqrt

(am

m)

−3 −2 −1 0 1 2 3

4−3

4−2

3−2

4−1

3−1

2−1

95% family−wise confidence level

Differences in mean levels of as.factor(spe.kmeans$cluster)

−4 −2 0 2 4 6 8

4−3

4−2

3−2

4−1

3−1

2−1



−1.0 −0.5 0.0 0.5

4−3

4−2

3−2

4−1

3−1

2−1



Clusters

Val

ue

env_c1 env_c2 env_c3 env_c4

02

46

env_c1 env_c2 env_c3 env_c4sp_c1 0 0 1 2sp_c2 7 4 0 1sp_c3 2 4 0 0sp_c4 0 3 5 0

The table

23

1.3.2. Ensables de especies

# -----------------# Ensambles de especies -----------------# El# objetivo es identificar asociaciones de especies en nuestros datos.

# 1) Abundancias medias en los clusters de sitios (k-means) -------# Calcularemos estadísticos simples (abundancias medias) a partir de# las tipologías obtenidas con el análisis cluster, y miraremos qué# especies están presentes en mayor abundancia o qué especies son# específicas de cada cluster de sitios

groups <- as.factor(spe.kmeans$cluster)spe.means <- matrix(0, ncol(spe), length(levels(groups)))row.names(spe.means) <- colnames(spe)for (i in 1:ncol(spe)) {

spe.means[i, ] <- tapply(spe[, i], spe.kmeans$cluster, mean)}# abundancias promedio de especies para los 4 gruposgroup1 <- round(sort(spe.means[, 1], decreasing = TRUE), 2)group2 <- round(sort(spe.means[, 2], decreasing = TRUE), 2)group3 <- round(sort(spe.means[, 3], decreasing = TRUE), 2)group4 <- round(sort(spe.means[, 4], decreasing = TRUE), 2)# especies con abundancias mayores que la abundancia media de# especies del grupogroup1.domin <- which(group1 > mean(group1))group1 #para el grupo 1

## ABL CHE GAR GOU GRE HOT VAN BRO PSO ROT BBO CHA TRU VAI## 3.33 1.33 1.33 1.00 1.00 0.33 0.33 0.33 0.33 0.33 0.33 0.00 0.00 0.00## LOC OMB BLA TOX BAR SPI PER BOU CAR TAN BCO PCH ANG## 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

group1.domin #especies con abundancias mayores que la media

## ABL CHE GAR GOU GRE## 1 2 3 4 5

# lo mismo podemos hacer con los demás gruposgroup2.domin <- which(group2 > mean(group2))group2 #para el grupo 2

## TRU VAI LOC OMB CHE CHA BLA VAN GOU BRO TAN BAR PER GAR## 4.00 4.00 3.58 1.00 1.00 0.92 0.75 0.58 0.50 0.42 0.33 0.25 0.25 0.08## HOT TOX SPI BOU PSO ROT CAR BCO PCH GRE BBO ABL ANG## 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00


## TRU VAI LOC OMB CHE CHA BLA## 1 2 3 4 5 6 7

group3.domin <- which(group3 > mean(group3))group3 #para el grupo 3

24

## LOC GAR VAI VAN CHE TOX GOU PER BAR SPI BLA TAN TRU BRO## 3.67 3.17 2.83 2.83 2.50 2.33 1.83 1.83 1.67 1.67 1.50 1.50 1.33 1.33## ABL BOU CHA HOT PSO CAR ROT ANG OMB GRE BCO BBO PCH## 1.17 0.83 0.67 0.67 0.67 0.67 0.50 0.50 0.33 0.33 0.17 0.17 0.00


## LOC GAR VAI VAN CHE TOX GOU PER BAR SPI BLA TAN## 1 2 3 4 5 6 7 8 9 10 11 12

group4.domin <- which(group4 > mean(group4))group4 #para el grupo 4

## ABL GAR GOU GRE TAN BAR BBO BOU BRO CHE BCO PSO ANG PER## 5.00 4.88 4.38 4.12 4.00 3.75 3.62 3.50 3.25 3.12 3.12 3.00 3.00 2.75## CAR VAN PCH SPI ROT HOT TOX LOC VAI TRU OMB BLA CHA## 2.62 2.25 2.25 2.12 2.12 1.62 1.50 1.00 0.38 0.12 0.12 0.12 0.00


## ABL GAR GOU GRE TAN BAR BBO BOU BRO CHE BCO PSO ANG PER CAR## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

# 2) Coeficiente W de Kendall para evaluar la concordancia -------# *********************************************** El objetivo es# encontrar los ensambles más globales, es decir, el menor número de# grupos que contendan el mayor número de asociaciones positivas y# significativas de especies. Legendre (2005) propuso el coeficiente# de concordancia W de Kendall para, conjuntamente con test de# permutación, identificar ensambles de especies en datos de# abundancia (no para datos de presencia-ausencia). Primero se# utiliza un test de independencia de todas las especies si# rechazamos la hipótesis nua, miramos grupos de especies# correlacionadas y, dentro de cada grupo, evaluamos la constribución# de cada especie al estadístico global, usando un test de# permutación. aquí no usamos información de tipología de sitios.

# primero extraemos las especies más abundantes del conjunto de# datos, las clasificamos con el método de pertición de k-medias, y# corremos un test global (kendall.global) para conocer si todos los# grupos de especies (’judges’) están asociados globalmente de manera# significativa si este es el caso, corremos test a posteriori# (kendall.post) en las especies de cada grupo para verificar si# todas las especies en un grupo son concortantes

# Extraemos la especies más abundantessp.sum <- apply(spe, 2, sum)spe.sorted <- spe[, order(sp.sum, decreasing = TRUE)]spe.small <- spe.sorted[, 1:20]

# Transformamos los datos de especiesspe.small.hel <- decostand(spe.small, "hellinger")

25

spe.small.std <- decostand(spe.small.hel, "standardize")spe.small.t <- t(spe.small.std)

# Realizamos una partición de k-meansspe.t.kmeans.casc <- cascadeKM(spe.small.t, inf.gr = 2, sup.gr = 8, iter = 100,

criterion = "calinski")plot(spe.t.kmeans.casc, sortg = TRUE)# la mejor elección es k=2 nos quedamos con la partición de 2 gruposclusters <- spe.t.kmeans.casc$partition[, 1]clusters

## LOC VAI GAR TRU ABL CHE GOU TAN VAN BAR BRO GRE PER BOU BBO PSO SPI## 2 2 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1## ANG TOX BCO## 1 1 1

# Realizamos un análisis de concordancia(spe.kendall.global <- kendall.global(spe.small.hel, clusters))

## $Concordance_analysis## Group.1 Group.2## W 4.545e-01 8.364e-01## F 1.333e+01 1.022e+01## Prob.F 5.222e-43 3.487e-13## Corrected prob.F 1.044e-42 3.487e-13## Chi2 2.164e+02 7.026e+01## Prob.perm 1.000e-03 1.000e-03## Corrected prob.perm 2.000e-03 2.000e-03#### $Correction.type## [1] "holm"#### attr(,"class")## [1] "kendall.global"

# Miramos los p-valores corregidos para el test de permutación como# ambos grupos tienen p<0.05, decimos que todos los grupos son# significativos a nivel global i.e. de manera global, contienen# especies que son concordantes (al menos algunas especies son# concordantes) Si no fuera significativo para algunos grupos,# indicaría que esos grupos deberían ser subdivididos en menores# grupos.

# Realizamos test a posteriori para identificar las especies# significativamente concordantes en cada grupo:(spe.kendall.post <- kendall.post(spe.small.hel, clusters, nperm = 9999))

## $A_posteriori_tests_Group## $A_posteriori_tests_Group[[1]]## GAR ABL CHE GOU TAN VAN BAR## Spearman.mean 0.3777 0.4563 -0.09825 0.4707 0.4288 0.2341 0.4528## W.per.species 0.4143 0.4883 -0.03365 0.5019 0.4624 0.2792 0.4850## Prob 0.0018 0.0002 0.75900 0.0001 0.0001 0.0413 0.0004

26

## Corrected prob 0.0072 0.0020 0.75900 0.0020 0.0020 0.0826 0.0020## BRO GRE PER BOU BBO PSO SPI## Spearman.mean 0.3088 0.5156 0.4343 0.5823 0.5172 0.5637 0.5353## W.per.species 0.3494 0.5441 0.4675 0.6069 0.5456 0.5893 0.5626## Prob 0.0094 0.0001 0.0003 0.0001 0.0001 0.0001 0.0001## Corrected prob 0.0282 0.0020 0.0020 0.0020 0.0020 0.0020 0.0020## ANG TOX BCO## Spearman.mean 0.5848 0.4729 0.5419## W.per.species 0.6093 0.5040 0.5689## Prob 0.0001 0.0001 0.0001## Corrected prob 0.0020 0.0020 0.0020#### $A_posteriori_tests_Group[[2]]## LOC VAI TRU## Spearman.mean 0.7414 0.8232 0.6976## W.per.species 0.8276 0.8822 0.7984## Prob 0.0001 0.0001 0.0001## Corrected prob 0.0020 0.0020 0.0020###### $Correction.type## [1] "holm"#### attr(,"class")## [1] "kendall.post"

# miramos los coef. de cor. de Spearman y su p-valor corregido.# Dentro de cada grupo, las especies concordantes deben estar# positivamente asociadas con las otras Una especie con una cor.# Spearman media negativa con las demás especies de su cluster está# clasificada de manera incorrecta. lo que sugiere que debe quitarse# del grupo. Si muchas especies tienen cor neg, el grupo debe# dividirse. Especies con p-valores no sig. no están contribuyendo a# la concordancia global del grupo, y deben quitarse.

# --- 4) Indicadores de especies (Dufrene and Legendre)---- Dufrêne# and Legendre (1997) propusieron un indicador de valores de especies# dentro de los clusters de sitios. El índice IndVal combina la# abundancia media de especies con sus frecuencias de ocurrencias en# los grupos. Un valor alto indica gran abundancia media dentro del# grupo en comparación con los demás grupos (especificidad) y# presencia en la mayoría de sitios de ese grupo (fidelidad).

# Dave Roberts, resume el concepto: “The indval approach looks for# species that are both necessary and sufficient, i.e. if you find# that species you should be in that type, and if you are in that# type you should find that species.”

# Vamos a agrupar los sitios con datos independientes (variables# ambientales), así el indicador de especie indicará especies que# están muy relacionadas a las condiciones ecológicas de su grupo.# Con un test de permutación a posteriori evaluaremos la

27

# significación del indicador (la probabilidad de obtener por azar un# valor del indicador como el observado)

# Dividimos los sitios en 4 grupos según la distancia al orígen del# río (das). Nos preguntamos si la variable ’das’ puede subdividirse# en grupos donde detectar especies indicadoras.

das.D1 <- dist(data.frame(das = env[, 1], row.names = rownames(env)))set.seed(1)dasD1.kmeans <- kmeans(das.D1, centers = 4, nstart = 100)dasD1.kmeans$cluster

## 1 2 3 4 5 6 7 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24## 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3## 25 26 27 28 29 30## 3 4 4 4 4 4

# Indicador de especies para esta tipología de sitios(iva <- indval(spe, dasD1.kmeans$cluster))

## $relfrq## 1 2 3 4## CHA 0.000 0.8 0.0000 0.0## TRU 0.875 0.9 0.0000 0.2## VAI 0.875 1.0 0.3333 0.2## LOC 0.875 1.0 0.5000 0.8## OMB 0.000 0.7 0.0000 0.2## BLA 0.000 0.7 0.0000 0.2## HOT 0.000 0.3 0.6667 1.0## TOX 0.000 0.4 0.5000 0.8## VAN 0.375 0.6 0.6667 1.0## CHE 0.625 0.9 1.0000 1.0## BAR 0.000 0.6 0.5000 1.0## SPI 0.000 0.4 0.5000 1.0## GOU 0.375 0.7 0.8333 1.0## BRO 0.500 0.5 0.6667 1.0## PER 0.375 0.4 0.5000 1.0## BOU 0.000 0.3 0.5000 1.0## PSO 0.000 0.4 0.6667 1.0## ROT 0.125 0.1 0.6667 1.0## CAR 0.000 0.4 0.5000 1.0## TAN 0.500 0.5 0.5000 1.0## BCO 0.000 0.1 0.5000 1.0## PCH 0.000 0.0 0.3333 1.0## GRE 0.000 0.2 0.8333 1.0## GAR 0.375 0.4 1.0000 1.0## BBO 0.000 0.1 0.6667 1.0## ABL 0.000 0.3 1.0000 1.0## ANG 0.000 0.3 0.5000 1.0#### $relabu## 1 2 3 4## CHA 0.00000 1.00000 0.00000 0.00000## TRU 0.52124 0.44788 0.00000 0.03089## VAI 0.42299 0.50759 0.04338 0.02603

28

## LOC 0.40385 0.42692 0.07692 0.09231## OMB 0.00000 0.87500 0.00000 0.12500## BLA 0.00000 0.90000 0.00000 0.10000## HOT 0.00000 0.13636 0.45455 0.40909## TOX 0.00000 0.38889 0.27778 0.33333## VAN 0.14325 0.27831 0.21828 0.36016## CHE 0.16484 0.19181 0.25974 0.38362## BAR 0.00000 0.18932 0.31553 0.49515## SPI 0.00000 0.23810 0.23810 0.52381## GOU 0.05660 0.14717 0.32075 0.47547## BRO 0.15152 0.07576 0.22727 0.54545## PER 0.15419 0.12335 0.26432 0.45815## BOU 0.00000 0.08621 0.25862 0.65517## PSO 0.00000 0.07792 0.25974 0.66234## ROT 0.06726 0.02691 0.31390 0.59193## CAR 0.00000 0.09091 0.22727 0.68182## TAN 0.11706 0.08027 0.26756 0.53512## BCO 0.00000 0.02069 0.27586 0.70345## PCH 0.00000 0.00000 0.14286 0.85714## GRE 0.00000 0.02857 0.28571 0.68571## GAR 0.12235 0.09788 0.30995 0.46982## BBO 0.00000 0.01734 0.28902 0.69364## ABL 0.00000 0.07095 0.42230 0.50676## ANG 0.00000 0.06122 0.20408 0.73469#### $indval## 1 2 3 4## CHA 0.000000 0.800000 0.00000 0.000000## TRU 0.456081 0.403089 0.00000 0.006178## VAI 0.370119 0.507592 0.01446 0.005206## LOC 0.353365 0.426923 0.03846 0.073846## OMB 0.000000 0.612500 0.00000 0.025000## BLA 0.000000 0.630000 0.00000 0.020000## HOT 0.000000 0.040909 0.30303 0.409091## TOX 0.000000 0.155556 0.13889 0.266667## VAN 0.053718 0.166985 0.14552 0.360164## CHE 0.103022 0.172627 0.25974 0.383616## BAR 0.000000 0.113592 0.15777 0.495146## SPI 0.000000 0.095238 0.11905 0.523810## GOU 0.021226 0.103019 0.26730 0.475472## BRO 0.075758 0.037879 0.15152 0.545455## PER 0.057819 0.049339 0.13216 0.458150## BOU 0.000000 0.025862 0.12931 0.655172## PSO 0.000000 0.031169 0.17316 0.662338## ROT 0.008408 0.002691 0.20927 0.591928## CAR 0.000000 0.036364 0.11364 0.681818## TAN 0.058528 0.040134 0.13378 0.535117## BCO 0.000000 0.002069 0.13793 0.703448## PCH 0.000000 0.000000 0.04762 0.857143## GRE 0.000000 0.005714 0.23810 0.685714## GAR 0.045881 0.039152 0.30995 0.469821## BBO 0.000000 0.001734 0.19268 0.693642## ABL 0.000000 0.021284 0.42230 0.506757## ANG 0.000000 0.018367 0.10204 0.734694

29

#### $maxcls## CHA TRU VAI LOC OMB BLA HOT TOX VAN CHE BAR SPI GOU BRO PER BOU PSO## 2 1 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4## ROT CAR TAN BCO PCH GRE GAR BBO ABL ANG## 4 4 4 4 4 4 4 4 4 4#### $indcls## CHA TRU VAI LOC OMB BLA HOT TOX VAN CHE## 0.8000 0.4561 0.5076 0.4269 0.6125 0.6300 0.4091 0.2667 0.3602 0.3836## BAR SPI GOU BRO PER BOU PSO ROT CAR TAN## 0.4951 0.5238 0.4755 0.5455 0.4581 0.6552 0.6623 0.5919 0.6818 0.5351## BCO PCH GRE GAR BBO ABL ANG## 0.7034 0.8571 0.6857 0.4698 0.6936 0.5068 0.7347#### $pval## CHA TRU VAI LOC OMB BLA HOT TOX VAN CHE BAR## 0.001 0.021 0.001 0.011 0.007 0.005 0.052 0.283 0.128 0.029 0.009## SPI GOU BRO PER BOU PSO ROT CAR TAN BCO PCH## 0.013 0.009 0.006 0.031 0.001 0.001 0.007 0.001 0.001 0.003 0.001## GRE GAR BBO ABL ANG## 0.002 0.015 0.001 0.009 0.001#### attr(,"class")## [1] "indval"

# relfrq=frecuencia relativa de especies en cada grupo# relabu=frecuencia relativa de abundancia de especies en los groupss# idval=valor indicador de cada especie

# Tabla de especies indicadoras significativasgr <- iva$maxcls[iva$pval <= 0.05]iv <- iva$indcls[iva$pval <= 0.05]pv <- iva$pval[iva$pval <= 0.05]fr <- apply(spe > 0, 2, sum)[iva$pval <= 0.05]fidg <- data.frame(group = gr, indval = iv, pvalue = pv, freq = fr)fidg <- fidg[order(fidg$group, -fidg$indval), ]fidg

## group indval pvalue freq## TRU 1 0.4561 0.021 17## CHA 2 0.8000 0.001 8## BLA 2 0.6300 0.005 8## OMB 2 0.6125 0.007 8## VAI 2 0.5076 0.001 20## LOC 2 0.4269 0.011 24## PCH 4 0.8571 0.001 7## ANG 4 0.7347 0.001 11## BCO 4 0.7034 0.003 9## BBO 4 0.6936 0.001 10## GRE 4 0.6857 0.002 12## CAR 4 0.6818 0.001 12## PSO 4 0.6623 0.001 13## BOU 4 0.6552 0.001 11

30

## ROT 4 0.5919 0.007 11## BRO 4 0.5455 0.006 18## TAN 4 0.5351 0.001 17## SPI 4 0.5238 0.013 12## ABL 4 0.5068 0.009 14## BAR 4 0.4951 0.009 14## GOU 4 0.4755 0.009 20## GAR 4 0.4698 0.015 18## PER 4 0.4581 0.031 15## CHE 4 0.3836 0.029 25

# para saber a qué especies refieren cada abreviación, podemos:# data(doubs) y preguntar sobre estos datos ?doubs Vemos que el# cluster 2 tiene como especie indicadora el TRU (Salmo trutta# fario), etc.

# Podemos exportar los resultados en un archiv csvwrite.csv(fidg, "IndVal-das.csv")

31

5 10 15 20

K−means partitions comparison

Objects

Num

ber

of g

roup

s in

eac

h pa

rtiti

on

23

45

67

8

5 10 15 20 8.0 9.0 10.0

calinskicriterion

Values

23

45

67

8

8.0 9.0 10.0

32

2 Cluster Peces

Documents

Transcript of 2 Cluster Peces