jueves, 15 de noviembre de 2012

R: Curvas de Respuesta Principal (PRC) 2. Una gráfica elegante

En el pasado post hablábamos de cómo conseguir que los resultados de la función prc del paquete vegan se parecieran lo más posible a los obtenidos con CANOCO. A pesar de ello la gráfica que se obtiene directamente con la función plot queda bastante fea. Aparte del tema de colores, tamaños de fuentes y ausencia de puntos, el principal problema es la diferencia de rangos de los valores de las puntuaciones de las parcelas (observaciones, puntos de muestreo,etc.) y de las especies. De hecho, Lepš y Šmilauer utilizan un método gráfico en CANOCO que separa los ejes de ambos, de forma que pueden así tener dos escalas distintas, sin importar más que luego coincidan los dos ceros en la gráfica general.

Esta es la misma solución que propongo para obtener una gráfica más elegante y más sencilla de interpretar visualmente. Lo podemos conseguir haciendo la gráfica con plot por un lado sin dibujar las especies con el parámetro species, y haciendo una gráfica independiente para estas con linestack. Ambas gráficas tendremos que fusionarlas luego con alguna aplicación externa. Vamos a ello.

Lo primero es dejar la parte de las parcelas visualmente aceptable. Para ello necesitaremos establecer varias parámetros gráficos previos con par, después dibujarla sin la leyenda con plot y, por último, añadirle la leyenda con legend en la posición más adecuada.

library (vegan)
##Extraemos los parámetros gráficos por defecto para no perderlos
opar <- par (no.readonly=T)
##Establecemos los parámetros gráficos para que nos haga todo en negro,
##aumente un poco el tamaño de las fuentes y las líneas,suprima los
##márgenes exteriores,nos ponga márgenes adecuados para los títulos
##de los ejes y los números de las escalas, nos deje los números de
##las escalas a 0.5 y en el sentido de lectura y establecemos el
##tamaño de las marcas de escala
png ('PRC.sites2.png')
par (bty='l', col='black', cex=2, cex.axis=1, cex.lab=1.3, lwd=1.8,
     xaxs='i',yaxs='i', mar=c(3,3,1,4), mgp=c(2,0.5,0), omi=c(0,0,0,0), 
     adj=0.5, las=1, tcl= -0.3)
##Dibujamos el gráfico sin especies y sin leyenda, con líneas y puntos y
##con límites de escala adecuados
plot (mod, scaling=1, species=F, col='black',type='b',ylab=c('PRC'),
      pch=c(15,16,17,18,19), ylim=c(-1,1.8),xlim=c(-5,25),legpos=NA)
##Dibujamos la leyenda sin el nivel de control, cuidando que los
##caracteres de los puntos coincidan con los graficados
legend (0.05,0.2,legend=levels(dose)[-1],pch=c(15,16,17,18,19),
        ncol=2,cex=1.5)
dev.off()


Ahora toca dibujar el eje de las especies. Si recordáis, en el post anterior extrajimos del análisis las puntuaciones de las especies a un vector con el escalado adecuado. Por si acaso, aquí está otra vez, solo que sin multiplicar por la constante correspondiente a CANOCO:

##Extraemos las puntuaciones de las especies para el primer eje
##RDA (choices=1), multiplicadas por la constante (const) que
##utiliza CANOCO
scores.sps <- scores (mod, choices=1, scaling=1,
                      dis='sp')

Vamos a hacer la gráfica con un poco menos de especies para que no salga tan abigarrada, así que estableceremos con colSums un límite inferior de 250, en vez de 100 como en el ejemplo de la función prc. Vamos a dibujar el eje de las especies con la función linestack, a la cual hay que especificarle que nos tiene que dibujar el eje con el parámetro axis. Con los parámetros air y hoff vamos a hacer que los nombres de las especies se separen un poco entre ellos y del eje. Por último, con el parámetro at colocaremos el eje en el centro de la gráfica. Podéis poner en el vector names los nombres completos de las especies, pero eso ya os lo dejo a vosotros.

##Dibujo del eje de especies
png ('PRCsps.png', 240, 480, res=100)
linestack (scores.sps[logabu>250], labels=names[logabu>250], axis=T,
           air=1.3, hoff=6, at=-1, font=3)
dev.off()

Y ya solo tenemos que unir los dos gráficos en uno con algún programa externo. La función linestack tiene un parámetro, add, con el que añadir el eje al gráfico actual, pero si hacemos esto, por la diferencia de escalas, no nos quedará bien colocado respecto al 0 de la figura de las observaciones. Por eso es mejor hacer los dos gráficos independientemente y luego juntarlos con un programa externo. Yo lo he hecho con Libre Office Draw.

Y, en principio, ya está. Sugerencias y comentarios serán bienvenidos.