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_1Este 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 resultadosnames.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)
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)
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)
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.
