domingo, 8 de diciembre de 2013

Código R para diseccionar los efectos del ganado en una comunidad vegetal mediterránea

Como prometí, voy a ir intentando publicar en Muestra Aleatoria el código de R utilizado en los distintos análisis de nuestro último artículo publicado en PLoS ONE.

Sigo practicando con el paquete knitr y el lenguaje de marcado markdown, así que esta entrada la he hecho por entero en RStudio. Con ello quiero resaltar que algunos parámetros gráficos están adaptados para que salga bien en el blog y tendréis que cambiarlos seguramente para obtener otros resultados en los dispositivos que elijáis para dibujar en R.

Como los datos son míos y hago con ellos lo que quiera, tenéis las matrices colgadas en el servicio Google Drive. Luego, al ordenar a R leer los archivos, recordad cambiar las rutas de los mismos a las que tengan en vuestros equipos, claro:

  • Con_Sps_1 Este es el archivo con la matriz de coberturas de especies por tratamiento, parcela y año en formato CANOCO. Este análisis lo hice al principio con CANOCO y como el paquete vegan trae una función, read.cep (), para importar ficheros producidos por CANOCO, me decidí a hacerlo así para comparar resultados

  • names.csv Es el fichero con los nombres de las especies.

Y ahora vamos a ver cómo se hizo el análisis de Principal Response Curves:

library(vegan)  ##Cargamos el paquete vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.0-9
sps.table <- read.cep("~/Documentos/UAM/Experimentos/Analisis_Cerro San Pedro/Sampling/CANOCO/Sin_Xolantha/Con_Sps_1_Con_Lavandula/Sin_Control/Con_Sps_1", 
    force = T)  ##Metemos en un objeto los datos de la matriz
log.sps.table <- log1p(sps.table)  ##Transformamos los datos a su logaritmo
center.log.sps.table <- scale(log.sps.table, scale = F)  ##Centramos por especies

## Definimos dos objetos con los factores de tratamiento y año con sus
## respectivas etiquetas

Treatment <- factor(rep(c("C", "D", "T", "F", "DT", "DF", "FT"), 28))
Year <- gl(4, 49, labels = c("2005", "2006", "2007", "2008"))

## Hacemos el análisis PRC

curves1 <- prc(center.log.sps.table, Treatment, Year)

Dibujamos, por un lado, la gráfica sin las especies, puesto que, por los distintos escalados, luego es más fácil dibujar el eje de las especies aparte con la función linestack y pegarlo a la gráfica de las curvas con otro programa para que sea de fácil interpretación.

plot(curves1, scaling = 1, species = F, ylim = c(-0.6, 0.6), xlim = c(2004.5, 
    2008.5), pch = c(1:2, 15:18), col = "black", type = "b", ylab = c("PRC"), 
    cex = 2, legpos = NA)
legend(0.55, 0.43, legend = levels(Treatment)[-1], pch = c(1, 2, 15, 16, 17, 
    18), ncol = 2, cex = 1.5)

plot of chunk PRCsites

Ahora, dibujamos el eje de las especies sin Lavandula stoechas subsp. pedunculata. Esto lo tenemos que hacer así debido a la puntuación tan alta de esta especie, lo que nos haría perder la perspectiva en la escala del eje. Luego dibujaremos a Lavandula sola y, con algún programa de tratamiento gráfico, haremos un corte en el eje y la introduciremos.

Por otro lado, vamos a elegir solo las especies con una puntuación en el eje mayor de 0.2 en valor absoluto, puesto que si no, nos quedaría demasiado abigarrada y no podríamos leer los nombres de las especies. Consideramos que las especies con una menor puntuación en el eje no son relevantes para definir los cambios en las comunidades.

sps.scores <- scores(curves1, scaling = 1, choices = 1, display = "sp")
names.tot <- read.csv("names.csv", header = F)
sps.names.sel <- names.tot[abs(sps.scores) > 0.2]
sps.scores.sel2 <- sps.scores[abs(sps.scores) > 0.2]
linestack(sps.scores.sel2[-18], labels = sps.names.sel[-18], axis = T, air = 1.3, 
    hoff = 6, at = -1, font = 3)

plot of chunk Sps

Ahora, dibujamos a Lavandula.

linestack(sps.scores.sel2[18], labels = sps.names.sel[18], axis = T, air = 1.3, 
    hoff = 6, at = -1, font = 3)

plot of chunk Lav

Y finalmente hacemos los test de significatividad de los ejes mediante permutaciones de MonteCarlo. Las razones de hacer el análisis del segundo eje de forma distinta las tenéis explicadas en un post anterior.

## Significación del primer eje

anova(curves1, strata = Year, first = T, perm.max = 499)
## Permutation test for rda under reduced model
## Permutations stratified within 'Year'
## 
## Model: prc(response = center.log.sps.table, treatment = Treatment, time = Year)
##           Df   Var    F N.Perm Pr(>F)   
## RDA1       1  2.31 13.4    199  0.005 **
## Residual 168 28.96                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Significación del segundo eje

SamE <- curves1$CCA$u
SamE <- SamE[, 1]
seconaxis <- rda(center.log.sps.table ~ Treatment * Year + Condition(Year) + 
    Condition(SamE))
anova(seconaxis, strata = Year, perm.max = 499, first = T)
## Permutation test for rda under reduced model
## Permutations stratified within 'Year'
## 
## Model: rda(formula = center.log.sps.table ~ Treatment * Year + Condition(Year) + Condition(SamE))
##           Df   Var    F N.Perm Pr(>F)
## RDA1       1  0.54 3.13     99   0.78
## Residual 168 28.96

Y ya está… bueno… no. Ahora faltaría tratar las gráficas obtenidas con algún software de manipulación de gráficos (sugiero GIMP que es Open Source) para obtener finalmente esto:

Espero que os sirva de ayuda. Si tenéis alguna duda, ya sabéis, utilizad los comentarios para que todos nos beneficiemos.

sábado, 30 de noviembre de 2013

Diseccionando los efectos del ganado sobre la vegetación


A pesar de no disponer de mucho tiempo y encontrarme en una situación algo complicada, vamos consiguiendo sacar adelante algunas cosas. Acabamos de publicar el que creo que ya será el último artículo con los resultados de mi tesis doctoral:


Durante cuatro años realizamos un experimento en el Cerro de San Pedro, en la Comunidad de Madrid, en el que simulamos tres actividades del ganado sobre un cantuesar sin pastoreo frecuente desde hace más de 50 años: defoliación, pisoteo y deposición de heces. El objetivo del estudio era comprobar los cambios producidos por cada una de las actividades en la composición florística y funcional en la vegetación e intentar dilucidar la importancia de los efectos indirectos (cambios en las características edáficas y lumínicas) y directos (destrucción de tejido vegetal) de las mismas.

Existía ya algún trabajo similar como los de Kohler et al. (2004 y 2006), pero ninguno sobre un sistema abandonado en el medio mediterráneo. Los resultados indicaron principalmente que no parecía haber ninguna diferencia en los cambios producidos en la vegetación entre los distintos tratamientos, excepto la deposición de heces, el cual no produjo ningún cambio notable ni en la composición específica ni en la funcional. Todos los tratamientos convergieron en comunidades similares al final del experimento, caracterizadas por la casi totalidad de desaparición del cantueso y una composición de rasgos funcionales habitual de los sitios pastoreados: anuales, con roseta basal, de poca altura y con semillas ligeras. La sorpresa surgió con el SLA (Specific Leaf Area), cuyos valores disminuyeron significativamente con los tratamientos.

Parece que todos estos cambios se debieron principalmente a los efectos directos de las actividades simuladas, es decir, a la destrucción de tejido vegetal y, en un segundo plano, a los cambios producidos en las condiciones lumínicas; las condiciones edáficas no sufrieron ningún cambio a lo largo del tiempo que duró el experimento.

Aparte de todo esto, quiero resaltar que analicé los datos mediante los modernos y complicados modelos mixtos, a los cuales llegué después de darle vueltas y vueltas a diferentes modelos de análisis de la varianza que conjugasen bloque, medidas repetidas y efectos fijos que no se ajustaban al algo complejo diseño experimental que teníamos. Lo digo por si alguien anda devanándose los sesos también con los modelos mixtos; creo que aquí tendrá un buen ejemplo de cómo y cuándo aplicarlos y, sobre todo, cómo presentar e interpretar los resultados. El campo, lo sabemos bien los ecólogos, nunca nos pone las cosas fáciles a la hora de luego trabajar con los datos y los modelos mixtos han venido para facilitarnos un poco la vida y tener resultados más precisos y fiables. No les tengáis miedo.


Por último, y también relacionado con los análisis estadísticos, utilicé el maravilloso y apasionante software R por primera vez en un artículo, lo que creo que también encontraréis interesante los que andéis empezando con él. No sé si tendré tiempo de publicar aquí los scripts de manera comprensible de los modelos mixtos, las PRC y las gráficas, aunque la intención la tengo. De todas formas, dentro del tiempo de que disponga, atenderé las dudas que alguien pueda tener sobre ello; utilizad los comentarios y sed pacientes, así nos beneficiaremos todos.