La PRCF ¿qué es y cómo se calcula?

Recientemente me he puesto a analizar unos datos de series temporales de dinámica de poblaciones de procesionaria (estimadas a partir de mediciones del grado de infestación de rodales) en Mora de Rubielos, Teruel, en el marco de una colaboración a tres bandas entre el Laboratorio de Sanidad Vegetal de Mora de Rubielos (Rodólfo Sánchez, Gerárdo Sánchez), la Universidad de Granada (José Antonio Hódar, Regino Zamora) y la Universidad Rey Juan Carlos (en mi humilde persona).

Para poder discernir entre factores endógenos y exógenos que podrían estar causando las oscilaciones de las series, una de las primeras cosas que hay que hacer es tratar de determinar si hay dependencia de la densidad y con cuantos retardos. Es decir, si el hecho de que, por ejemplo, uno, dos o tres años atrás, la población estuviera en un máximo, va a influir en que ahora, en tiempo t, la población haya sufrido un descenso en el número de individuos como consecuencia del agotamiento de los recursos u otros factores que podrían estar operando a nivel intraespecífico. Pues bien, para determinar con cuantos retardos la propia variable tiene un efecto sobre la observación de esta misma variable en tiempo t, Alan Berryman (1999) sugiere el uso de lo que llama la PRCF (partial rate correlation function), que no traduzco porque no se si lo haría bien. Para aquellos que estéis un poco familiarizados con la jerga de series temporales, no os será desconocido la función de autocorrelación (ACF) y la función de autocorrelación parcial (PACF). Estas son herramientas exploratorias básicas para identificar la estructura de una serie temporal. Pues bien, en el caso de dinámica de poblaciones, Berryman sugiere que es más apropiado utilizar la PRCF... ¿y qué es esto de la PRCF? Se trata de correlacionar la tasa de cambio poblacional en tiempo t, con la densidad poblacional en tiempo t-1, t-2 y así sucesivamente. La tasa de cambio de la población es lo que Berryman denomina como variable Rt, y que en realidad es el log(Nt/Nt-1). Es decir, una medida de cómo la población cambia entre un periodo y el siguiente.

En mi periplo analítico, me encontré con que no tenía ni idea de cómo calcular la PRCF en R, que es ya el único software con el que trabajo. Busqué en los foros y nada. Lo más que encontré fue la página de un software (PAS) que diseñó Alan Berryman ya hace tiempo y que explicaba con fórmulas cómo calcular la PACF y, a partir de ahí, la PRCF. Pero honestamente, no entendí mucho. Así que di en mis indagaciones con un artículo de Mauricio Lima publicado en PLoS ONE donde decía explícitamente que ellos habían calculado la PRCF con un código escrito en R. Escribí a Mauricio y este rápidamente me dirigió a Sergio Estay, quien con mucha diligencia me pasó el código, que con permiso suyo y agradeciéndole de antemano su enorme ayuda, pongo aquí a disposición de todos aquellos interesados (he modificado un par de cosas para generalizarlo, pero en esencia, todo el código es suyo).

Código para calcular la función PRCF en R 

Al final, y viendo el código de Sergio (y con alguna aclaración suya) entendí que era esto de la PRCF... y es bastante simple. Para más de un retardo la PRCF es equivalente a la PACF. Y para un retardo es simplemente la correlación entre Rt (la tasa de cambio) y Nt-1. La (última) pregunta que yo me hice fue ¿por qué para dos o más retardos la PRCF es igual a la PACF? Y dándole vueltas mientras volvía a casa en coche, llegué a una conclusión, que expongo a continuación y que no se si será más o menos acertada (si alguien conoce una respuesta más acertada por favor que me corrija).

Básicamente lo que se está diciendo es que para dos retardos (y generalizable a otros retardos por encima de uno):

r(Rt ~ Nt-2) habiendo quitado el efecto de Nt-1 sobre Rt 
r(Nt ~ Nt-2) habiendo quitado el efecto de Nt-1 sobre Nt 

¿Tiene esto algún sentido? Lo tiene si pensamos que Rt se calcula como el log(Nt/Nt-1). Por tanto si quitamos el efecto de Nt-1 sobre Rt, lo que nos queda es básicamente Nt. Por tanto, con dos o más retardos, la PRCF será equivalente a la PACF.

Y ya para terminar, incluyo aquí los resultados de la PRCF para las 11 parcelas con las que estamos trabajando en Mora de Rubielos, y en donde se ve que la dinámica de la procesionaria está regulada por efectos de densidad con un único retardo.

 

Nueva versión del paquete TPL (v. 1.1)

No se ha hecho esperar la nueva versión 1.1. del paquete TPL (anteriormente v. 1.0) para realizar la estandarización taxonómica de nombres de plantas utilizando una conexión en línea al portal de The Plant List. Gracias a los comentarios de varios colegas, he corregido unos cuantos bugs, afinado los mecanismos de búsqueda de sinónimos, variedades y subespecies y, aprovechando la coyuntura, he reducido significativamente el código (en aproximádamente un 40%) sin perder funcionalidad. Esto último no es de interés para el usuario estándar, pero para aquellos que quieran utilizar el código fuente para modificarlo con posterioridad, sí será de gran utilidad, ya que el código ahora es algo más legible y está mejor estructurado.

He probado el paquete con varias bases de datos de plantas, una propia que he sacado del proyecto BIOTREE-NET, con más de 5000 nombres de árboles para toda Centroamérica, una lista de 3047 nombres de plantas (procedentes del Banco de Datos de Biodiversidad de la Generalitat Valenciana) que me ha pasado Jaume Tormo, una lista de 1122 nombres de briófitos que me ha pasado Íñigo de la Cerda-Granzow, una lista de 238 árboles de la Reserva del Triunfo, en México, que me ha pasado Neptalí Ramírez-Marcial, y un listado de 1188 nombres de árboles de la Amazonia procedentes de Kalle Ruokolainen (disponible en el paquete de R 'betaper'). El código ha funcionado bien para todas ellas.

¡Ojo! Esto no quiere decir necesariamente que el código no cometa errores a la hora de buscar la información dentro de TPL (no es tan sencillo como parece,... sólo hay que mirar el código fuente para ver la gran variedad de situaciones que podemos encontrarnos). Aunque he hecho comprobaciones manuales (remitiéndome a TPL) con nombres elegidos al azar y aparentemente todo funciona bien, sería conveniente que los usuarios finales hagan verificaciones de este tipo, sobre todo en estas primeras etapas, y me comuniquen si encuentran alguna incoherencia entre los resulutados obtenidos en R y lo que realmente dice TPL.

Pero cuidado con lo que interpretamos como error. En la lista que me pasó amablemente Jaume Tormo, encontré algunos nombres que, aunque existían en TPL, la función los intepretaba como que no estaban allí (argumento 'Plant.Name.Index' = FALSE). Algunos ejemplos de ésto son Xanthium strumarium o Hyoseris scabra. ¿Qué ocurre entonces? La razón es, que al acceder a TPL desde R, lo que hacemos es leer un archivo *.csv. Y en algunos casos estos archivos están mal formateados (problema de entrada de los datos en TPL). En estos dos casos particulares, la información de 'Genus' y 'Species' estaban en las columnas de 'Genus.hybrid.marker' y 'Species.hybrid.marker' respectivamente. Con esto, la función acaba por no identificar el nombre y lo da como que no está en TPL. Lo malo es que tampoco es fácil identificar esto como un caso de tabla mal formateada de forma automática, por lo que la columna 'WFormat' arroja, de forma incorrecta, un valor FALSE. Se podría hacer algo al respecto pero la variabilidad en la forma en la que las tablas vienen mal formateadas es tan grande que incorporar todos estos posibles casos al código sería una pesadilla. En cualquier caso, son pocos nombres a los que les ocurre esto. Lo mejor es comprobar los casos que tienen el argumento 'Plant.Name.Index' = FALSE uno por uno y ver si alguno de ellos está realmente en TPL y no ha sido registrado por estos problemas.

Otra de las mejoras de esta nueva versión del paquete TPL es que el código corrige ahora con más exactitud los errores tipográficos en los nombres de las plantas. Los errores tipográficos son más comunes de lo que creemos. Es fácil poner una letra de más o de menos en un nombre complejo. A veces, es el subconsciente el que nos traiciona y acabamos llamando, por ejemplo, Marsilea bastardae a lo que viene siendo Marsilea batardae. De los listados revisados anteriormente, en todos había errores tipográficos, y la función TPL permite corregir con bastante exactitud muchos de ellos (siempre y cuando los errores estén en el epíteto específico). En el listado de más de 5000 plantas del proyecto BIOTREE-NET, identificamos, por ejemplo, 299 errores, que fueron corregidos automáticamente. Se pueden ver todos estos nombres y el consiguiente output de la función TPL aquí. Aparentemente todo está bien.

Bueno, creo que ya me he enrollado bastante. El paquete se puede descargar aquí en *.tar.gz (Linux): 

TPL_v.1.1.tar.gz 

o *.zip (Windows): 

TPL_v.1.1.zip 

Pues nada más... por favor, cualquier incidencia, sugerencia o comentario hacédmelo llegar. Y si todo funciona bien, también agradecería comentarios al respecto (e información sobre el número de nombres con el que se ha probado la función, etc).

¡Cruzad los dedos y... a correr el código!

Un método automatizado para la estandarización taxónomica de nombres de plantas con The Plant List (TPL)

Cuando se trabaja con grandes bases de datos de vegetación con procedencia muy diversa, la taxonomía puede jugarnos una mala pasada. En estas bases de datos es frecuente encontrar: (1) especies que se llaman de distinta forma pero que en realidad son la misma (sinónimos); (2) especies con el mismo nombre pero que en realidad son distintas (homónimos); (3) errores tipográficos, que hacen muchas veces que registremos como distintas, especies que ya han sido registradas en la misma base de datos. En consecuencia es necesario estandarizar la taxonomía. The PlantList (TPL) es una iniciativa que permite poner un poco de orden en todo este caos que suponen los sistemas nomenclaturales (aunque obviamente no está exenta de errores).

En TPL se puede ingresar el nombre de una planta y te dice si dicho nombre está aceptado, es sinónimo de otro o está todavía sin resolver. En caso de que no aparezca en TPL es muy probable que haya un error tipográfico, aunque a veces simplemente ocurre que el nombre en cuestión todavía no ha sido ingresado en la base de datos. La principal limitación con el uso de este portal web, es que la validación hay que hacerla nombre por nombre, lo que supone una carga de trabajo muy alta cuando tenemos listas enormes de nombres de plantas.

Pues bien, en el marco del proyecto BIOTREE-NET, una iniciativa que trata de compilar y estandarizar información de árboles en inventarios forestales para toda Centroamérica, he diseñado una función en R que te permite hacer todo este trabajo de forma automatizada. Esto ahorra mucho trabajo manual y puedes cotejar miles de nombres en poco tiempo (unas horas como mucho). La función (TPL) se presenta dentro de paquete de R con el mismo nombre, y utiliza a su vez otra función (TPLck) que hace el cotejo para nombres individuales.

El procedimiento está resumido en el siguiente esquema (que he preparado para una posible publicación):

El paquete se puede descargar aquí en *.tar.gz (Linux):

TPL_v.1.0.tar.gz

o *.zip (Windows):

TPL_v.1.0.zip

El resultado de aplicar la función TPL a un listado de nombres de plantas es un arreglo de datos (data.frame) con información sobre el estatus del nombre según The Plant List, el nombre actualizado (en caso de sinónimos o errores tipográficos), la familia y la autoridad, entre otras cosas.

Las principales limitaciones hasta el momento son que: (1) no es posible resolver el problema de los homónimos; (2) las correcciones de errores tipográficos sólo se realizan cuando los errores existen en el epíteto específico. Si los errores se producen en el género o en el epíteto infraespecífico, no hay nada que hacer; (3) en el caso de que una especie sea sinónima de otra y mantenga el mismo nombre cambiando sólo la autoridad (ej. Bartramia pomiformis var. elongata Turner como sinónimo de Bartramia pomiformis var. elongata Hedw.), el procedimiento extráe finalmente la autoridad del sinónimo (Turner) y no del nombre aceptado (Hedw.); (4) en el caso de que haya muchos epítetos infraespecíficos y ninguno se corresponda con el nombre que estamos validando, la función busca por defecto el nombre de la especie que NO tenga epíteto infraespecífico. Si este nombre no existe en TPL el nombre se queda sin resolver (pero se da una advertencia al usuario para que pueda revisar este nombre a posteriori).

Si alguien prueba el paquete y detecta algún error o se le ocurre alguna mejoría posible, por favor que no dude en contactarme.

Cómo estimar parámetros en un modelo por máxima verosimilitud (versión fuerza bruta)

Una parte importante de los modelos estadísticos es la forma en la que se estiman los parámetros. En modelos lineales (regresión, ANOVA, ANCOVA), la estimación se hace frecuentemente por el método de mínimos cuadrados. En modelos lineales generalizados (GLM) la estimación se hace por el método de máxima verosimilitud. Cómo opera este método es algo normalmente desconocido al usuario estándar de programas estadísticos (incluido R). Aunque hay muchas formas de estimar los parámetros de un modelo utilizando máxima verosimilitud, la idea subyacente es la misma: encontrar los parámetros que maximizan la probabilidad de los datos observados. La forma más intuitiva de hacer esta estimación es por fuerza bruta, esto es, probando muchas posibles combinaciones de los valores de los parámetros y ver qué combinación maximiza la probabilidad de los datos. A continuación pongo un breve manual sóbre cómo hacer esto en R con objeto de poder entender un poco mejor este proceso. 

Cómo estimar los parámetros de un modelo por máxima verosimilitud en R

No obstante, esta forma iterativa es computacionalmente inpracticable, como se explica en el manual, y distintos métodos se han ideado para solucionar estos problemas y encontrar parámetros óptimos (locales y globales). No entraré a detallar los distintos métodos de estimación por máxima verosimilitud, pero si entendemos la idea de fondo, creo que será mucho más fácil entender el proceso por el que generamos nuestros modelos y algunos de los problemas (por ejemplo, convergencia) que a veces surgen por el camino.

Modelos a la carta: El modelo gaussiano

Los modelos lineales generalizados (GLM) tienen tres componentes: (1) la estructura de errores; (2) el estimador lineal y; (3) la función de vínculo. La función de vínculo linealiza la relación entre la respuesta y el estimador lineal. La función de vínculo, junto con el estimador lineal, constituyen lo que lo algunos autores denominan el "modelo científico". La máxima versatilidad en la modelización se alcanza cuando uno ajusta modelos a la carta, fuera de las restricciones que los GLM imponen. Ello implica que debemos escribir nuestras propias funciones, para después ajustarlas a los datos. El proceso se completa en el contexto de la selección de modelos por criterios de información cuando comparamos toda una batería de modelos plausibles ajustados por métodos de máxima verosimilitud. En R esto puede hacerse por medio del paquete 'likelihood' escrito por Lora Murphy. El paquete todavía no está disponible en CRAN, aunque es de esperar que lo esté pronto. Mientras tanto, este paquete se puede descargar de la página del curso de "Modelos y métodos de máxima verosimilitud" de Charlie Canham. Este curso fue por cierto impartido en abril en Granada y he de decir que ha sido uno de los mejores cursos a los que he asistido... ¡gracias Charlie!

Bueno, a lo que iba... en lo que respecta a los modelos científicos. Podemos usar muchos tipos de modelos dependiendo de los datos que tengamos. No tienen por qué ser lineales. Una familia de modelos muy versátiles que permiten ajustar respuestas que aumentan primero y disminuyen después son los modelos gaussianos. La función gaussiana se expresa así:

$y=a\cdot e^{(\frac{-{(X-b)}^2}{(2\cdot c^2)})}$

Dicha función consta de 3 parámetros, que llamaremos 'a', 'b' y 'c', pero que podríamos haber denominado de cualquier otra forma. El parámetro 'a' indica el valor máximo que alcanza la curva. El parámetro 'b' determina la forma de la curva: monotónica, creciente, decreciente. El parámetro 'c' va a condicionar la pendiente de la curva (lo que en inglés se denomina sharpness).

Veamos ahora cómo cambios en los valores de estos parámetros determinan cambios en la forma de las curvas tipo. Tomemos un ejemplo en dónde vamos a modelar la abundancia de arácnidos en sistemas agrícolas (y) en función del tiempo (X). La idea con el siguiente código en R es que nos familiaricemos con esta función y sus parámetros.

Si modificamos el parámetro 'a' podemos obtener algunas las siguientes curvas:

X <- 1:250
a <- 1000
b <- 125
c <- 0.01
Y1 <- a*exp(-((X-b)^2/2*c^2))
plot(Y1~X, type="l", ylim=c(0,1050), ylab="Abundancia de arácnidos", xlab="Tiempo", main="Función gausiana")
text(x=125, y=1025, "a=1000")
a <- 800
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=2)
text(x=125, y=825, "a=800")
a <- 600
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=3)
text(x=125, y=625, "a=600")
a <- 400
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=4)
text(x=125, y=425, "a=400")
a <- 200
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=125, y=225, "a=200")

Modificando el parámetro 'b' obtenemos cambios en la forma de la curva gaussiana como se observa en la siguiente gráfica:

X <- 1:250
a <- 1000
b <- 125
c <- 0.01
Y1 <- a*exp(-((X-b)^2/2*c^2))
plot(Y1~X, type="l", ylim=c(0,1050), ylab="Abundancia de arácnidos", xlab="Tiempo", main="Función gausiana")
text(x=125, y=1025, "b=125")
b <- 200
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=2)
text(x=200, y=1025, "b=200")
b <- 275
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=4)
text(x=240, y=825, "b=275")
b <- 50
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=50, y=1025, "b=50")
b <- -25
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=15, y=825, "b=-25")

Cambios en el parámetro 'c' van a provocar cambios en las pendientes de la curva, como se observa a continuación:

X <- 1:250
a <- 1000
b <- 125
c <- 0.01
Y1 <- a*exp(-((X-b)^2/2*c^2))
plot(Y1~X, type="l", ylim=c(0,1050), ylab="Abundancia de arácnidos", xlab="Tiempo", main="Función gausiana")
text(x=45, y=825, "c=0.01")
c <- 0.001
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=2)
text(x=125, y=1025, "c=0.001")
c <- 0.25
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=4)
text(x=125, y=150, "c=.025")
c <- 0.1
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=90, y=250, "c=0.1")
c <- 0.02
Y1 <- a*exp(-((X-b)^2/2*c^2))
lines(X, Y1, lty=5)
text(x=40, y=400, "c=0.02")

Como simular un bosque en 3D

Ahora que estoy trabajando con datos espaciales de infestación de pinos por procesionaria y muérdago, estaba buscando la mejor manera de representar los datos gráficamente. Tradicionalmente estos datos se representan en superficies de dos dimensiones, pero ¿se pueden representar de una forma eficiente en 3D? La respuesta es sí y la solución nos la da (posiblemente es sólo una de las muchas soluciones) el paquete scatterplot3d de R y la función con el mismo nombre.

Veamos unos datos que simulan los datos con los que yo estoy trabajando. Coordenadas x e y con la posición espacial de los datos. Una variable z con información de la altura de los datos. Y finalmente, una cuarta variable con información sobre, por ejemplo, el número de bolsones de procesionaria (pero podría ser cualquier otra cosa).

x <- rep(1:10, each=10) + rnorm(100, 0, 0.2)
y <- rep(1:10, 10) + rnorm(100, 0, 0.2)
z <- rnorm(100, 6, 1)
procesionaria <- rpois(100, lambda = 5)
procesionaria.col <- cut(procesionaria, breaks=c(0, 5, 10, 50, 100), labels=c("greenyellow", "green2", "forestgreen", "darkgreen"))
table(procesionaria.col)

library(scatterplot3d)
scatterplot3d(x = x, y = y, z = z, type="h", cex.symbols=5, pch=19, color=procesionaria.col, xlab="", ylab="", zlab="Altura", zlim=c(0,10))

Created by Pretty R at inside-R.org

Y el resultado sería este:

Modificando los argumentos de la función scatterplot3d() podemos hacer que el color de los árboles sea representativo del grado de infestación por procesionaria (más oscuro más infestado). También podríamos, si quisiéramos, hacer que el símbolo de copa representara la especie.

Por cierto, gracias Antonio por enseñarme Pretty R para incorporar código de R en mis entradas. ¡Con todo el tiempo que llevo en ésto y todavía no lo conocía!

Problemas de convergencia en el escalamiento multidimensional no métrico (NMDS) ¿Qué hacer?

El escalamiento multidimensional no métrico (NMS, MDS, NMDS o NMMDS) es una técnica multivariante de interdependencia que trata de representar en un espacio geométrico de pocas dimensiones las proximidades existentes entre un conjunto de objetos. El NMDS es un método de ordenación adecuado para datos que no son normales o que están en una escala discontinua o arbitraria. Una ventaja del NMDS frente a otras técnicas de ordenación es que, al estar basada en rangos de distancias, tiende a linealizar la relación entre las distancias ambientales y las distancias biológicas (esto es, calculadas a partir de una matriz de sitios x especies). Una de las desventajas de esta técnica es la di cultad para alcanzar una solución estable única. A pesar de ello, el NMDS es una técnica ampliamente utilizada en ecología para detectar gradientes en comunidades biológicas.

El NMDS se implementa de la siguiente forma:

1. Se calcula la matriz de disimilaridad X a partir de la matriz de datos de sitios x especies. Esta matriz nos indica cómo de iguales son cada par de sitios utilizando para ello la similaridad entre sus especies. Supongamos que tenemos tres especies (sp1, sp2, sp3) y tres sitios (A, B, C). El sitio A tiene sp1 = 3, sp2 = 0 y sp3 = 8. El sitio B tiene sp1 = 3, sp2 = 0 y sp3 = 6. El sitio C tiene sp1 = 0, sp2 = 5 y sp3 = 1. Por tanto, podemos calcular una matriz de disimilaridad que nos indique con números que los sitios A y B son muy iguales, mientras que los sitios A y C y B y C son muy distintos entre sí. Cuando se trata de datos biológicos la distancia más usada es la distancia de Sorensen (Bray-Curtis) en vez de la distancia Euclídea.
2. Se asignan los sitios (unidades muestrales) a una con guración inicial aleatoria en un espacio k-dimensional (dónde k es el número de especies), aunque en realidad, la ordenación se va a realizar principalmente sobre unas pocas dimensiones (2 o 3).
3. Se calculan las distancias sobre este nuevo espacio geométrico y se calcula una matriz de distancia Y .
4. Se comparan las matrices de distancia X e Y y se mide cómo son de parecidas entre ellas (stress).
5. A partir de la con guración inicial, se reasignan los sitios (unidades muestrales) para reducir las distancias con la matriz X.
6. Se repite este proceso de manera iterativa hasta que se consigue una solución óptima en dónde la matriz de distancias Y es muy parecida a la matriz de distancias X. Esto es, se minimiza el stress.

En R tenemos una implementación de esta función (metaMDS) en el paquete vegan.

Problemas de convergencia

Uno de los problemas que pueden surgir a la hora de realizar un NMDS en R es el de la convergencia a la hora de encontrar una solución óptima. Esto ocurre porque existen dos sitios (o más) que tienen exactamente la misma composición de especies. Es decir, todos los sitios tienen que tener una composición de especies única (puede ser muy parecida pero no exactamente iguales). Todavía no se muy bien si esto es una premisa del NMDS como tal o simplemente un resultado de su implementación en R. Lo que está claro es que el NMDS está basado en distancias y puede darse el caso de que la distancia entre dos sitios sea 0 (independientemente del método que utilicemos para calcular dichas distancias), lo que implicaría que su composición de especies sea exactamente igual. Por lo tanto no es un problema de cálculo de distancias sino del algoritmo que realiza el NMDS.

Aunque desconozco todavía el motivo de este error, existe una forma de solucionar el problema. La idea es sencilla. Consiste en incluir un pequeño ruido aleatorio en los números que forman la matriz de datos. Este error debe de ser lo suficientemente pequeño como para no alterar las distancias multidimensionales entre pares de sitios, pero conseguirán que los valores de composición no sean exactamente iguales. Y aunque al tratarse de abundancias de especies los datos deberían de ser enteros, no supone un problema análitico el meter un ruido aleatorio que sea de tipo decimal. Veámoslo con un ejemplo.

Creamos una matriz de datos de abundancia con 10 sitios y 10 especies:

m <- matrix(rpois(100, 1), nrow=10, ncol=10)
colnames(m)<- paste("sp", c(1:10), sep="")
rownames(m) <- paste("sitio", c(1:10), sep="_")
m

sp1 sp2 sp3 sp4 sp5 sp6 sp7 sp8 sp9 sp10
sitio_1 1 1 0 4 4 2 1 1 1 0
sitio_2 2 1 2 1 0 3 1 0 1 0
sitio_3 0 0 2 1 0 1 0 0 0 1
sitio_4 2 1 2 0 2 1 1 0 0 0
sitio_5 0 1 0 1 0 0 0 1 0 0
sitio_6 1 1 0 0 0 1 2 1 1 0
sitio_7 2 1 1 1 0 2 0 2 1 0
sitio_8 1 0 3 1 0 2 2 1 1 0
sitio_9 0 1 1 2 1 3 0 1 0 0
sitio_10 0 2 6 1 1 2 0 0 1 1

La analizamos en R con la función metaMDS del paquete vegan.

library(vegan)
set.seed(0)
nmds1 <- metaMDS(m)
plot(nmds1)

Como vemos, todo va bien. Ahora vamos a crear una nueva fila (sitio) en la matriz m que tendrá la misma abundancia que el sitio_1.

m2 <- rbind(m, m[1,])
rownames(m2) <- paste("sitio", c(1:11), sep="_")
m2

sp1 sp2 sp3 sp4 sp5 sp6 sp7 sp8 sp9 sp10
sitio_1 1 1 0 4 4 2 1 1 1 0
sitio_2 2 1 2 1 0 3 1 0 1 0
sitio_3 0 0 2 1 0 1 0 0 0 1
sitio_4 2 1 2 0 2 1 1 0 0 0
sitio_5 0 1 0 1 0 0 0 1 0 0
sitio_6 1 1 0 0 0 1 2 1 1 0
sitio_7 2 1 1 1 0 2 0 2 1 0
sitio_8 1 0 3 1 0 2 2 1 1 0
sitio_9 0 1 1 2 1 3 0 1 0 0
sitio_10 0 2 6 1 1 2 0 0 1 1
sitio_11 1 1 0 4 4 2 1 1 1 0

Y analizamos esta nueva matriz de datos con la función metaMDS.

set.seed(0)
nmds2 <- metaMDS(m2)

Error in metaMDSdist(comm, distance = distance, autotransform = autotransform, :
Zero dissimilarities are not allowed

Resolvemos el problema. Vamos a meter ruido únicamente en los valores > 0.

m3 <- m2
for(i in 1:dim(m2)[2]){
m3[, i] <- ifelse(m2[,i] > 0, jitter(m2[,i],
factor=1e-05), m2[,i])
}

Y ahora analizamos la matriz m3 con la función metaMDS().

set.seed(0)
nmds3 <- metaMDS(m3)
plot(nmds3)

El único problema a este enfoque es la replicabilidad exacta de los resultados. Al utilizar la función set.seed(0) hacemos que la configuración inicial del NMDS sea siempre la misma. Sin embargo, al meter un ruido aleatorio con la función jitter() cada vez que corremos el código vamos a crear matrices ligeramente distintas y el gráfico de ordenación que obtengamos será distinto cada vez, aunque la posición relativa de los puntos va a ser muy parecida siempre. Si queremos evitar este problema podríamos salvar la matriz m3 como un archivo de texto y leer esos mismos datos cada vez que vayamos a repetir el análisis. Esto solucionaría el problema.

Si alguien averigua por qué ocurre este problema de convergencia en R que por favor me lo diga. Se puede encontrar más información sobre análisis multivariante aquí.

Análisis de datos ecológicos en R

Ya están disponibles seis manuales para el análisis de datos ecológicos en R. Algunos de ellos habían sido publicados con anterioridad en entradas anteriores (Introducción a R, gráficos en R, modelos lineales, modelos lineales generalizados). Esta nueva edición ofrece una versión actualizada y ampliada de estos manuales, incluyendo como novedad el análisis de diversos diseños de experimentos (bloques aleatorizados, medidas repetidas, split-plot, diseños jerarquizados) por medio de los modelos lineales mixtos.

Los manuales pueden descargarse aquí:

1. Una introducción a R (2.0 Mb).
2. Modelos lineales: Regresión, ANOVA y ANCOVA (621 Kb).
3. Modelos lineales generalizados (GLM) (416 Kb).
4. Introducción al diseño de experimentos (13.7 Mb).
5. Modelos lineales mixtos en R (1.9 Mb).
6. Análisis multivariante (1.9 Mb).
   

Cursos de análisis de datos en ecología

Para aquellos que usen estos manuales con fines docentes, he aquí algunas recomendaciones. Se puede organizar un curso general de "Análisis de datos en ecología" que incluya una introducción a R (1), modelos lineales, (2) GLM (3) y multivariante (6). Este curso puede durar entre 25 y 30 horas, siempre y cuando los alumnos tengan un conocimiento medio-alto de estadística. No es recomendable impartir este curso a alumnos de pre-grado. El perfil más adecuado del alumnado serían alumnos de postgrado, fundamentalmente doctorandos, investigadores y profesores.

Otro curso puede ser el de "Diseño de experimentos y análisis de datos con modelos lineales mixtos". Es recomendable que este curso incluya una introducción a R (1), modelos lineales (2), diseño de experimentos (4) y modelos lineales mixtos (5). Este curso puede durar otras 25-30 horas y estaría enfocado, al igual que el curso anterior, a alumnos de postgrado, investigadores y profesores, fundamentalmente del área de biología.

Seis métodos para incluir la autocorrelación espacial en el análisis de datos espaciales

Los datos espacialmente explícitos (ej. datos que muestran la distribución de las especies) generalmente manifiestan autocorrelación espacial. La autocorrelación espacial ocurre cuando los valores de las variables muestreadas en puntos cercanos no son independientes entre sí o, dicho de otro modo, cuando muestras próximas entre sí exhiben valores más parecidos que con muestras más alejadas. La causa principal de la autocorrelación espacial es la relación existente entre la distancia y determinados procesos biológicos como la especiación, la extinción, la dispersión o las interacciones entre especies. Existen dos causas más que pueden producir dependencia espacial de los residuos del modelo (no autocorrelación espacial en sentido estricto), pero a efectos estadísticos, los efectos de dicha dependencia espacial suponen los mismos problemas que los de la autocorrelación espacial. Estas son: (1) el intento de modelar linealmente relaciones no lineales entre la variable respuesta y las variables ambientales; y (2) la ausencia en el modelo de variables ambientales que están espacialmente estructuradas y, que por tanto, causan una estructura espacial en la variable respuesta (ej. variables climáticas).

La autocorrelación espacial es a la vez una oportunidad para explicar determinados procesos (ej. procesos de contagio, dispersión geográfica, organización social, etc.) y un reto para el análisis de datos espaciales, ya que los residuos de los modelos no son totalmente independientes y esto conlleva un aumento del error de tipo I (esto es, rechazar la hipótesis nula siendo cierta). Por ello se han desarrollado en los últimos años una gran variedad de métodos para corregir los efectos de la autocorrelación espacial. En este artículo se presentan y explican seis métodos concretos: el mapeo espacial de vectores propios (spatial eigenvector mapping, SEVM), generalización de mínimos cuadrados (generalised least squares, GLS), modelos autorregresivos condicionales (conditional autoregressive models, CAR), modelos autorregresivos simultáneos (simultaneous autoregressive models, SAR), modelos lineales generalizados mixtos (generalised linear mixed models, GLMM) y ecuaciones de estimación generalizadas (generalised estimation equations, GEE). También se discute en qué condiciones el uso de uno u otro modelo es más adecuado y se provee el código en R para implementar estas funciones.

Todos estos modelos asumen la existencia de estacionariedad (spatial stationarity) e isotropía (isotropic spatial autocorrelation). La estacionariedad se refiere al hecho de que la autocorrelación espacial es constante en el espacio. Esto no siempre es necesariamente cierto. Por ejemplo, en el caso de la capacidad de dispersión de un organismo, ésta podría cambiar al pasar de la llanura a la montaña, en dónde el movimiento está más restringido. La isotropía se refiere a qué la autocorrelación espacial actúa de la misma forma en todas las direcciones. Algunos factores ambientales que podrían causar anisotropía son el viento (dando a un organismo que se dispersa con el viento una dirección de movimiento preferente), las corrientes de agua (ej. en el movimiento del plancton) o la direccionalidad en el transporte del suelo a favor de pendiente.

F. Dormann, C., M. McPherson, J., B. Araújo, M., Bivand, R., Bolliger, J., Carl, G., G. Davies, R., Hirzel, A., Jetz, W., Daniel Kissling, W., Kühn, I., Ohlemüller, R., R. Peres-Neto, P., Reineking, B., Schröder, B., M. Schurr, F., & Wilson, R. (2007). Methods to account for spatial autocorrelation in the analysis of species distributional data: a review Ecography, 30 (5), 609-628 DOI: 10.1111/j.2007.0906-7590.05171.x

Lyx y Sweave

Lyx es un procesador de texto que promueve la edición de textos basándose en la estructura del documento, y no sólamente en su apariencia. En la terminología anglosajona ésto es lo que se denomina WYSIWYM (What You See Is What You Mean), frente al enfoque adoptado por otros procesadores, como Word u OpenOffice, que se encuadran más dentro del concepto WYSIWYG (What You See Is What You Get).

¿Qué significa ésto? Básicamente, que no tenemos que preocuparnos por editar lo que escribimos. Lyx, de manera inteligente, asume el rol de editor del documento y nosotros sólo tenemos que definir qué parte del texto es qué. De esta manera, si a Lyx le decimos que una parte del texto es el resumen, no hace falta darle un formato determinado, Lyx lo hace por nosotros.

En Lyx tenemos distintas clases que vienen establecidas por defecto. Las clases definen los argumentos de ese documento. Así tenemos una clase que es 'Artículo', otra que es 'Libro', otra 'Carta', y así sucesivamente. Dentro de cada clase hay distintos ambientes (Environments). Un ambiente es una parte del texto que va a tener un determinado formato, por ejemplo: un título, un resumen, un autor, una sección, las referencias, etc.

Principales diferencias de Lyx con otros editores de texto

1. Una de las principales características de Lyx es que no se pueden introducir espacios de más con el tabulador o el Enter. Esto es porque Lyx asume que nosotros no debemos dejar espacios para separar una determinada sección de otra. De nuevo, esta es la misión del programa. Así se busca que el usuario se centre en los contenidos mientras que el programa se encarga del formato.
2. Las secciones y subsecciones, así como los listados, se enumeran automáticamente si introducimos o eliminamos cualquier elemento nuevo en el texto.
3. Las figuras se ajustan al texto de manera automática. No tenemos que preocuparnos por su ubicación exacta y ésta se actualiza si modificamos el texto asociado.

Cómo instalar Lyx

Lyx es un editor de texto basado en Latex. Por tanto, para instalar Lyx en Ubuntu es conveniente que instalemos todos los paquetes relacionados con Latex. Si queremos que funcionen los acentos cuando escribamos en Español es importante que vayamos a Documents-Settings y especifiquemos que Language = Spanish y Encoding = utf8.

Lyx y Sweave

Una de las ventajas principales de Lyx (y por la cúal he empezado a aprender a utilizarlo) es que funciona bien con Sweave, un programa que permite escribir y decodificar código en R dentro de Lyx. Esto es de gran utilidad a la hora de escribir manuales de R o informes que contengan código o el resultado de la implementación de un código en R. Podemos elegir si queremos que el código que escribimos en el documento se vea, que sólamente se vea el resultado o ambos. Al actualizar cualquier parte del código se actualizan automáticamente los resultados obtenidos. Sin embargo, si el código contiene algún error, entonces no es posible decodificar y, por tanto, exportar a pdf o a cualquier otro formato, el documento. en cuestión Esto actúa en cierta manera como control de calidad del código que escribamos dentro del documento. Al escribir manuales de R en OpenOffice me he encontrado a veces con errores que se hacen patentes cuando los alumnos tratan de implementar el código, simplemente porque no he verificado correctamente todo el código al escribir el manual. Esto no pasa con Lyx.

Instalar Sweave

Para instalar Sweave basta con seguir las instrucciones dadas en INSTALL en la siguiente dirección:

http://cran.r-project.org/contrib/extra/lyx/

En Linux no hay que preocuparse mucho por los pasos 4 y 6. Basta con mirar que, una vez que se ha reconfigurado Lyx (Tools-Reconfigure), tengamos en Documents-Settings-Document class la clase 'article (Sweave-noweb)'.

Para empezar a utilizar código en R dentro de los documentos conviene empezar con algo que ya esté escrito, para ver cómo hay que empezar y finalizar las líneas de código. Encontraremos algunos ejemplos en la página anterior (Sweave-test-1.lyx y test.lyx). En esencia, para escribir algo de código dentro de un documento de Lyx, la estructura que hay que seguir es la siguiente:

<<>>=
hist(rnorm(100))
@

Se puede encontrar más información y algunos ejemplos de manuales de R escritos en Lyx en la página de Duncan J. Golicher.

¿Cómo analizamos los datos cuándo no conocemos las identidades de todas las especies?

Un problema típico al que se enfrentan los ecólogos que trabajan en regiones tropicales es el no poder conocer las identidades de todas las especies con las que trabajan. Ello es debido, por un lado, a la gran diversidad que existe en estas regiones y, por otro, a la falta de estudios taxonómicos detallados. Muchas veces, biólogos y ecólogos tienen que trabajar sin claves taxónomicas o con claves incompletas. En el sur de México, en dónde realicé mi tesis doctoral, no existen por ejemplo claves taxonómicas específicas de plantas y hay que identificar las especies usando la Flora de Guatemala y consultando especímenes de herbarios. Todo ello hace que, en los listados de especies que se generan como consecuencia del trabajo de campo, haya muchas especies que estén identificadas sólo a nivel de género, a nivel de familia o de las que no se tenga ni la menor idea del grupo al que pertenecen.

A la hora de realizar análisis estadísticos específicos para comparar la composición de especies entre distintos sitios (p.e. test de Mantel, RDA, CCA, MANOVA semi-paramétrico), este problema se puede solventar fácilmente si las especies son identificadas a nivel de morfoespecies, es decir, que sabemos que la especie A es distinta del resto de las especies en función de atributos morfológicos (tipo de hoja, fuste del tronco, flor, fruto, etc.), aún sin saber qué especie es o, a veces, ni siquiera la familia o grupo a la que pertenece. Para ello es necesario cruzar muestras de todas las especies colectadas de todos los sitios muestreados (en inglés 'cross-checking'). Esto resulta tremendamente laborioso. Además, hay ocasiones en las que para identificar ciertas especies (p.e. las de la familia Lauraceae) hace falta tener información de atributos muy específicos, como la flor o el fruto, los cuales no están muchas veces disponibles a la hora de realizar el trabajo de campo. En este caso, puede surgir la incertidumbre de si la especie que hemos llamado Persea A en la muestra 1 no sea la misma que hemos llamado Persea liebmanii en la muestra 2. Esta situación puede ser mucho más crítica cuando trabajamos con muestras colectadas por distintos investigadores o técnicos. Este es el caso típico de trabajos a una escala más regional. En esta situación la incertidumbre taxonómica es muchísimo mayor ya que es seguro que no ha habido cruce de la información de las especies colectadas en las distintas muestras.

En estos casos, las dos aproximaciones más comunes al análisis de datos multivariantes han consistido, o bien en eliminar las morfoespecies y/o especies no identificadas, o bien en llevar a cabo el análisis a nivel de género, en dónde la incertidumbre taxonómica suele ser mucho menor. En un trabajo realizado recientemente, hemos propuesto una alternativa estadísticamente mucho más robusta al análisis de datos multivariantes cuando existe incertidumbre taxonómica. Nuestro enfoque supone permutar las identidades de las especies no identificadas dentro del nivel taxónomico en el que se encuentran e iteratuar este procedimiento n veces, generando así no una sóla matriz de muestras x especies, sino n matrices. Posteriormente, calcularíamos el parámetro específico del análisis deseado sobre cada una de estas n matrices, obteniendo un rango de parámetros estimados que nos indicaría los posibles valores que podría tomar dicho parámetro ante distintos escenarios plausibles de incertidumbre taxonómica. Por ejemplo, en el test de Mantel, dicho parámetro podría ser el coeficiente de correlación de Pearson, r. En el RDA/CCA podría ser la cantidad de variabilidad explicada por las variables ambientales, y así sucesivamente.

Para implementar dichas funciones, hemos creado un paquete en R, 'betaper', que permite, por un lado generar las n matrices permutando las especies no identificadas (función 'pertables') y, posteriormente, aplicar distintas funciones a cada una de estas matrices, calculando el rango de parámetros de interés deseado. Hasta el momento hemos implementado los siguientes métodos multivariantes disponibles, todos ellos, en el paquete 'vegan': análisis de la varianza multivariante semi-paramétrico (función 'adonis.pertables'), test de Mantel (función 'mantel.pertables'), CCA (función 'cca.pertables') y RDA (función 'rda.pertables'). Todas ellas tienen una sálida gráfica ('plot').

Tomemos como ejemplo un grupo de nueve inventarios muestreados por uno de los autores de este trabajo (Dr. Kalle Ruokolainen) en la Amazonia. Estos datos están disponibles en el paquete 'betaper'. Vamos a estimar el efecto de la incertidumbre taxonómica sobre la varianza explicada en la composición de especies por una serie de variables edáficas (cationes de Ca, K, Mg y Na).

install.packages("betaper")
library(betaper)

data(Amazonia)
data(soils)

# Definimos un nuevo índice que incluye los términos usados en la base de datos Amazonia para definir especies no identificadas a diferentes niveles taxonómicos

index.Amazon <- c(paste("sp.", rep(1:20), sep=""), "Indet.", "indet.")

# Generamos un objeto 'pertables' (i.e. una lista de matrices permutando las especies no identificadas o morfoespecies)

Amazonia100 <- pertables(Amazonia, index=index.Amazon, nsim=100)

# Y ahora comprobamos el efecto de la incertidumbre taxonómica sobre la varianza explicada en la composición de especies por las variables edáficas en un RDA

Amazonia.rda <- rda.pertables(Amazonia100 ~., data=soils) Amazonia.rda

Confidence intervals of R-squared and pseudo-F values for RDA under different taxonomic scenarios

Rsquared pseudoF
0% 0.4754156 0.9062709
0.5% 0.4768863 0.9116456
2.5% 0.4802825 0.9241262
50% 0.4936813 0.9750405
97.5% 0.5058685 1.0237546
99.5% 0.5088126 1.0358837
100% 0.5091075 1.0371060

plot(Amazonia.rda)

Vemos que, a pesar de la alta incertidumbre taxonómica de esta base de datos (casi el 50% de las especies no están identificadas a nivel de especie), los efectos que ésta tiene sobre el parámetro estimado en este caso (variabilidad explicada por las variables edáficas) no varía mucho (entre 47% y 51%). La gráfica muestra los valores que tomaría cada una de las muestras en los dos primeros ejes del RDA bajo cada uno de los 100 escenarios de reasignación de las identidades de las especies. Todos los puntos pertenecientes a la misma muestra están agrupados dentro de una elipse por lo que es posible analizar visualmente como afecta la incertidumbre taxonómica a los valores de cada muestra. Las cruces señalan los valores estimados bajo una de los enfoques tradicionales consistente en eliminar las especies no identificadas y las morfoespecies del análisis. Podemos observar que no siempre las cruces caen dentro de las elipses por lo que es fácil deducir que los resultados y conclusiones a las que se llegaría utilizando este enfoque no se corresponden con ninguno de los escenarios plausibles de identidad de las especies.

Los resultados de este trabajo se encuentran actualmente en segunda revisión en la revista Ecography.

Cayuela, L., de la Cruz, M. & Ruokolainen, K. 2009. A method to incorporate the effect of taxonomic uncertainty on multivariate analyses of ecological data. Ecography, in 2nd rev.