11  Actividad 4

11.1 Caso real: Contaminación en bentos por plataformas petroleras.

El siguiente caso de estudio fue publicado por Gray JS, Clarke KR, Warwick RM, Hobbs G (1990). Detection of initial effects of pollution on marine benthos: an example from the Ekofisk and Eldfisk oilfields, North Sea. Marine Ecology Progress Series 66, 285-299. Este estudio consiste en la evaluación del macrobentos y varios contaminantes del sedimento en 39 sitios dispuestos en un diseño radial alrededor de una plataforma de perforación petrolera en el mar del Norte, donde se espera que los contaminantes asociados a la actividad petrolera afecten la estructura del ecosistema. La disposición de los sitios es circular, alejándose cada ciertos kilómetros del centro de perforación. Para esta actividad ignoraremos la distancia de cada sitio respecto a la plataforma. Intentaremos clasificar las muestras según sus descriptores, es decir, dejaremos que los datos cuenten su propia historia en vez de forzarlos a un modelo específico.

11.1.1 Análisis de clasificación o Cluster analysis. Evaluación de los métodos de agrupación

  1. Importe la matriz de datos “sedimentos”

  2. ¿Cuáles son las dimensiones del objeto?

  3. ¿Cuáles son las variables que allí aparecen y cuáles deberían incluirse en el análisis?

  4. Seleccione sólo las variables que ambientalmente describan la muestra de sedimento, incluyendo contaminantes. Llame a esa selección “sed”.

  5. Evalue las escalas de magnitud de las variables, así como sus promedios. ¿Qué función de R podría ser útil para esto?

  6. Considerarndo la escala de magnitud de las variables, qué tipo de pretratamiento debe ser aplicado?

  7. Aplique el pretratamiento y llame a la matriz resultante “sed.stand”

  8. Aplique distancias Euclideanas a la matriz pretratada y llame a ese objeto “DE.sed”

  9. Genere un dendograma jerárquico con el método de agrupamiento simple (“single linkage”). Para esto será necesario aplicar la función hclust. Le sugerimoms leer la documentación en R sobre hclust del paquete “Stats”. Por ahora le adelantamos que los argumentos necesario son la matriz de distancias y el método de vinculación, que en este caso será el simple. Llame al cluster “single.c” y grafíquelo. ¿Cuántos clusters (grupos) identifica. ¿Es capáz de reconocer algún patrón lógico? ¿Cuáles son los sitios que más se parecen?¿es este gráfico evidencia de contaminación alrededor de la plataforma petrolera?

  10. Genere un nuevo dendograma jerárquico, pero ahora con el método completo (“complete linkage”). Llame al cluster “comp.c” y grafíquelo. Responda: ¿cuántos grupos identifica? ¿hay algún patrón más lógico? ¿cambió la clasificación respecto al método anterior?

  11. Genere un tercer dendograma pero ahora con el método promedio (“average linkage”). Llame al cluster “ave.c” y grafíquelo. Responda: ¿cuántos grupos identifica? ¿hay algún patrón más lógico? ¿cambió la clasificación respecto a los dos métodos anteriores? Para tener una mejor perspectiva, puede proyectar los tres dendogramas usando la función par con una fila y tres columnas, y llamando nuevamente a los tres dendogramas. Es recomendable asignar un título a cada plot para poder hacer comparaciones efectivas.

  12. Modifique la posición de las etiquetas agregando como atributo en cada plot “hang = -0.1”. ¿Prefiere esta salida o la anterior?

Código 
par(mfrow=c(1,3),mar=c(2,4,4,2))
plot(ave.c, main="average linkage", hang = -0.1)
plot(comp.c, main="Complete linkage", hang = -0.1)
plot(single.c, main="Single linkage", hang = -0.1)
  1. Hasta ahora, la identifiación de grupos es compeltamente arbitraria. Usaremos un método computacionalmente intensivo para “podar” el dendograma. Esta función es la simprof del paquete {clustsig}. Copie en su cónsola ??simprof para obtener ayuda de cómo ejecutar simprof. Asigne un nombre al resultado simprof e dentifique cuántos grupos realmente distinguibles hay explorando ese objeto.
Código 
sig.clust <-
  simprof(
    sed.stand,
    num.expected = 500,
    num.simulated = 99,
    method.cluster = "average",
    method.distance = "euclidean",
    method.transform = "identity",
    alpha = 0.05,
    sample.orientation = "row",
    const = 0,
    silent = TRUE,
    increment = 100
  )

Grafique su resultado con la función simprof.plot
par(mfrow = c(1, 1))
sig.clust.dend <- simprof.plot(sig.clust, leaflab = "perpendicular")

11.1.2 Métodos de interpretación: Heat Maps para avariables biológicas

  1. Importe la matriz de datos “macrofauna.cvs” y renómbrela “matriz.macrofauna”

  2. Evalue el contenido de la matriz: unidades, escalas y dimensiones.

  3. Genere una nueva matriz pero filtrando las primeras dos columnas. Llámelas “macrofauna”

  4. Calcule descriptores univariados por cada sitio, use la riqueza de especies, la abundancia de individuos y el índice de diversidad de Simpson. Para ello genere un data.frame con las columnas como las variables, y filas como los sitios. Llame a la tabla generada “uni”. Trate de identificar patrones de variación en esto sestimadores según el sitio

  5. Para usar una aproximación multivariada, aplique el índice de similitud Bray-Curtis luego de transformar las abundancias a raíz cuadrada. Llame a la matriz de abundancias transformadas “r.c” y a la matriz Bray-Curtis resultante “bray1”. Trate de identifcar un patrón de asociación entre las muestras viendo a “bray1”.

  6. Genere un dendograma jerárquico con el método de agrupamiento “average linkage” o UPGMA. Para esto será necesario aplicar la función hclust. Llame al cluster “cluster.bray” y conviértalo en dendograma con nuevo nombre “dend.macro”. ¿Encuentra familar esa representacion? compárela con el dendograma ave.c de las variables ambientales en un solo plot. Identifique un patrón respecto a grupos y muestras que lo conforman. ¿Hay correspondencia en ambos análisis?

  7. Identifique la relación entre las matrices con la función Mantel. Esto nos dará una idea de cuán similares son las relaciones espaciales según sus descriptores. En teoría, si los contaminantes afectan el macrobentos, debe haber correlación entre las matrices. Ejecute la prueba Mantel usando el coeficiente de correlación Spearman y construya una hipótesis nula con 9999 permutaciones. mantel(bray1, DE.sed, method=“spearman”, permutations=9999)

  8. ¿Cuán importante cree usted que es la similitud entre ambas matrices? Sirve esta aproximación para asociar matrices de igual dimensión?

  9. Ahora apliquemos la prueba estadística “simprof” para podar el dendograma del macrobentos usando un método cuantitativo. ¿Cuántos grupos estadístcamente significativos identifica el análisis jerárquico usando el método promedio?

Código 
library(clustsig)
sig.clust.macro <-
  simprof(
    r.c,
    num.expected = 1000,
    num.simulated = 999,
    method.cluster = "average",
    method.distance = "braycurtis",
    method.transform = "identity",
    alpha = 0.05,
    sample.orientation = "row",
    const = 0,
    silent = TRUE,
    increment = 100
  )

sig.clust.dend.m <- simprof.plot(sig.clust.macro, leaflab = "perpendicular")

  1. Utilicemos ahora la función “heatmap” para apreciar la contribución de las especies a la estructuración de los grupos. Primero, como práctica, generemos un factor, con base en el grupo “distancia” en la matriz.macrofana.

  2. Identifique cuáles son los niveles del factor distancia generado llamando a ese objeto distancia.

  3. Lea en el menú de ayuda sobre la función “heatmap”. Esta función requiere definir varios argumentos, los más importantes son: la matríz con los valores de abundancia de especies y el dendograma final. Básicamente, la función reordena las filas en la matriz de abundancias, y asigna un color a las abundancias, que aumenta según la intensidad. Las columnas las asocia con base en un dendograma basado en correlacioens entre las especies. Luego genera un gráfico que facilita interpretar las especies que generan los grupos. Lo primero que debe hacer es convertir el dataframe “r.c” a una matriz de nombre “r.c.m”

Luego debe llamar a la función “heatmap”, incluyendo la matrix r.c.m y al dendograma “dend.macro”

  1. ¿Puede interpretar algo del HeatMap generado?

  2. Mejoresmos el heatmap incluyendo solo a las 10 especies que mejor contribuyen a diferenciar los grupos en las distancias 1 y 4. Para esto usaremos la herramienta SIMPER desarrollada por Clarke (1993), original del software PRIMER, pero incluída en la librearía Vegan. Esta rutina descompone las similitudes entres cada par de sitios y relativiza el peso de cada especie en contribuir a la disimilitud promedio entre dos grupos. Para más información:

Código 
??simper

Apliquemos “simper” a la matriz “r.c” usando el factor distancia. Para poder acceder al resultado, llamaremos a este “simper_1v4” y al resumen de los resultados sum1v4

  1. El objeto “sum1v4” es una lista, y como tal, podemos observar su elementos uno a uno. Identifiquemos las especies que mayor contribución a las diferencias entre los grupos 1 y 4 con el siguiente comando

  2. Generemos un vector de nombre “sp” con los nombres de las 20 especies con mayor contribución a las diferencias entre la distancia 1 y 4 y usémoslo como subsetting en la matriz “r.c”, que simultáneamente podemos convertir de data.frame a matriz con el nombre “r.c.sub”

  3. Genere un nuevo heatmap, pero usando la matriz de 20 especies. ¿Qué logra apreciar?

  4. Trate de mejorar el heatmap pero ahora permitiendo que las filas se agrupen según la distancia a la que pertenecen. ¿cómo haría esto?

  5. ¿Es posible reordenar las especies según su similitud? Ello se puede abordar usando el índice de Bray-Curtis luego de estandarizar las abundancias de las especies por sus totales y multiplicar por 100. Esto equivale a estimar el índice de asociación de Whittaker (1952), que es ideal para medir correlación entre especies cuyas abundancias están infladas por velores cero. Esto lo puede lograr con los siguientes códigos, los cuales debe analizar y describir.

Código 
total <- apply(macrofauna, MARGIN = 2, FUN = sum)

mac2 <- macrofauna

for (i in 1:nrow(macrofauna)){
  mac2[i,] <- 100*(macrofauna[i,]/total)  
}

d <- vegdist(t(mac2[,sp]), method = "bray")

cluster.sp <- hclust(d, method = "average")

dend.sp<-as.dendrogram(cluster.sp)

plot(dend.sp, main= "Similitud entre especies, Whittaker's index of association")

heatmap(r.c.sub, Rowv = distancia, labRow = distancia, Colv = dend.sp)