martes, 13 de noviembre de 2012

R: Curvas de Respuesta Principal (PRC) 1. El lío del escalado

Este post va a ser un poco más técnico y menos divulgativo, aviso. Son un poco mis apuntes de R, paquete estadístico libre con muchas ventajas pero difícil de utilizar para los que venimos de programas con interfaces gráficas. Estos posts de R espero que además sirvan de ayuda para aquellos que andan dándose de cabezazos con él.

Estoy trabajando con unos datos de campo para hacer unas curvas de respuesta principal (PRC). Existe un paquete para R desarrollado principalmente por Jari Oksanen llamado vegan. Tenía las PRC hechas antes con el programa se suele usar para estos casos: CANOCO. Sin embargo, los resultados con vegan diferían bastante de los que había obtenido con CANOCO.

El paquete vegan contiene una colección de datos denominada pyrifos para hacer prácticas con él. Curiosamente los datos de pyrifos son de un estudio de van den Brink y ter Braak1 en el que se analizan precisamente con PRC. Además, ter Braak ha participado en el desarrollo de la función prc del paquete vegan junto con Oksanen. Lo más gracioso es que haciendo el análisis con el código propuesto por el archivo de ayuda de la función prc se obtienen también resultados distintos a los que se exponen en el artículo.

Lo primero que hay que recordar para entender esta disparidad de resultados es que la PRC es un análisis multivariante en el que las escalas pueden ser muy arbitrarias porque lo que interesa son las distancias relativas entre los elementos en estudio (normalmente parcelas y especies). Jari Oksanen ha generado varios documentos respecto a su paquete vegan para aclarar ciertas dudas. Uno de ellos es Design decisions and implementation details in vegan. En él explica cuáles son las constantes que ha utilizado para hacer los escalados de especies y parcelas en los análisis de redundancias (RDA, las PRC son un caso especial de estas) comparándolos con las usadas en CANOCO. Además, hay que observar que el escalado se puede hacer centrado en especies o en parcelas y que la función prc nos da la opción de especificarlo con el parámetro scaling.

Veamos qué es lo que ocurre siguiendo el código propuesto por el archivo de ayuda de la función prc:

##Cargamos el paquete
library (vegan)
##Cargamos en el espacio de trabajo los datos
data(pyrifos)
##Generamos un factor con las semanas de aplicación del tratamiento
week <- gl(11, 12, labels=c(-4, -1, 0.1, 1, 2, 4, 8, 12, 15, 19, 24))
##Generamos otro factor con los tratamientos
dose <- factor(rep(c(0.1, 0, 0, 0.9, 0, 44, 6, 0.1, 44, 0.9, 0, 6), 11))
##Hacemos el análisis prc
mod <- prc(pyrifos, dose, week)
##Metemos en un vector la suma de los datos para cada especie
logabu <- colSums(pyrifos)
##Hacemos la gráfica solo con las especies cuya suma de datos sea superior
##a 100
png ('prc.man.png')
plot(mod, select = logabu > 100)
dev.off()





Vemos que no solo no encontramos las parcelas en la parte negativa del eje, sino que además el rango de la escala sobrepasa las 2 unidades cuando en la figura 3 del trabajo de van den Brink y ter Braak no sobrepasa el valor 1.80. También en el eje de especies vemos cómo la más alejada, caenhora (Caenis horaria), se encuentra por debajo del valor -3, cuando en la figura del artículo se encuentra por debajo del 5. Esto aparte de lo feo de la figura y la saturación del eje de especies. Y esta gráfica obtenida con R es la que se repite en multitud de sitios. Hay que poner un poco más de cuidado.

Lo primero que hemos de saber es que el escalado puede estar centrado en las especies o en las parcelas, unidades de muestreo o como lo queramos llamar. Esto lo especificábamos en CANOCO y también lo podemos especificar en R con el parámetro scaling. Los valores para este parámetro están perfectamente explicados en el archivo de ayuda de la función plot.cca. En el caso de la PRC y siguiendo el método descrito en Šmilauer y Lepš2, centramos por parcelas, así que scaling=1.

png ('prc.scale1.png')
plot(mod, select = logabu > 100, scaling=1)
dev.off()



Aún así, vemos que los resultados siguen si parecerse mucho a los del artículo, sobre todo en las especies si prestamos atención a la más alejada. Y es que, como explica Oksanen en Design decisions and implementation details in vegan, las puntuaciones de especies y parcelas se multiplican por una constante distinta en CANOCO y en vegan. Para obtener al menos los mismos valores absolutos en CANOCO y en vegan, tenemos que dar un rodeo algo más largo, porque plot.cca no permite especificar dicha constante, pero scores, sí. Vamos a hacerlo con las especies para comprobarlo. En CANOCO esta constante es la raiz cuadrada del número de especies.

##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, const=sqrt(132), dis='sp')
##Extraemos solo los valores para las especies con más presencia,
##así como sus nombres
scores.sps.sel <- scores.sps [logabu>100]
names <- colnames (pyrifos)
names.sel <- names [logabu>100]
##Dibujamos el eje de las especies de la PRC
png ('spsCANOCO.png')
linestack (scores.sps.sel, axis=T, labels=names.sel)
dev.off()




Ahora, nuestra especie de referencia por ejemplo, caenhora, cae en la zona del eje que en valor absoluto más se asemeja a los resultados del artículo. Si además contáramos con las puntuaciones exactas de las especies en CANOCO, veríamos que coinciden con las obtenidas por scores con la constante establecida. Y lo mismo para las puntuaciones de cada combinación de semana y dosis, solo que aquí la constante es la raiz cuadrada del número de observaciones.

Quedaría por saber por qué narices en vegan obtenemos resultados con el signo cambiado, algo que no he podido entender y que espero que algún lector más avispado o entendido que yo u Oksanen mismo lo puedan explicar. Sin embargo, da un poco igual porque, como decíamos al principio, lo importante es la posición relativa de los objetos de estudio en el espacio, sea en la parte positiva o negativa del eje. Queda para un próximo post la solución que he encontrado para hacer el gráfico de las curvas un poco más elegante.

1Van den Brink, P. J. and Braak, C. J. F. T. (1999), Principal response curves: Analysis of time-dependent multivariate responses of biological community to stress. Environmental Toxicology and Chemistry, 18: 138–148.

2Jan Lepš, Petr Šmilauer. Multivariate Analysis of Ecological Data using CANOCO. Cambridge University Press 2003.