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.

Más sobre fumigaciones aéreas y procesionaria

Nuestro artículo publicado recientemente en Forest Ecology and Management sobre la inutilidad de las fumigaciones aéreas para el control de la procesionaria ha suscitado algunas críticas. La principal crítica (reiterada por el revisor de otro artículo que tenemos actualmente en revisión en Climatic Change utilizando esta misma base de datos) se refiere al hecho de que nuestro artículo sólo se refiere a rodales que sufren un grado de infestación de procesionaria alta (nivel 3 o superior sobre una escala de 5). En estos casos es lógico pensar que la plaga ya está al borde del colapso poblacional y, que por tanto, los rodales fumigados van a tener una respuesta similar a los no fumigados: esto es, un colapso da las poblaciones de procesionaria al siguiente año (se puede encontrar una explicación más detallada en esta otra entrada). Sin embargo, puede haber rodales que tengan un grado de infestación medio (por ejemplo 2) y que también sean fumigados. Se podría pensar que en estos rodales la respuesta no es la misma, y que mientras que en los fumigados se rompería el ciclo poblacional, los no fumigados tendrían más probabilidades de sufrir una superpoblación al siguiente año.

Explicaré a continuación qué parte de esta crítica tiene sustento y qué parte no la tiene. Cuando comenzamos este trabajo, tomamos la información de la Consejería de Medio Ambiente de la Junta de Andalucía. En las propias directrices de actuación de la Junta se establecía que sólo los rodales con un nivel de infestación de 3 o más son tratados, mientras que los que sufrían un nivel 2 sólo eran tratados si estaban próximos a un rodal con un grado de infestación de 3 o más. Esta información está disponible en el siguiente enlace (ver páginas 58 a 64). Como esta era información oficial, la dimos por válida sin cuestionarla. Más tarde, cuando surgieron las críticas, procedimos a verificar que, efectivamente, se fumigaban mayoritariamente rodales con un grado de infestación de 3 o más. La siguiente tabla muestra, para el período 2002-2005 (que es aquel para el que disponemos de la información más completa), el número de rodales fumigados (Trat.) y sin fumigar.

	>=3		2		1		0		Total
	Trat.	Sin	Trat.	Sin	Trat.	Sin	Trat.	Sin	Trat.	Sin
2002	47	271	23	414	44	985	29	1550	149	4240
2003	29	371	30	629	38	1074	36	1404	135	4254
2004	29	409	22	482	23	1010	39	1666	118	4271
2005	34	184	16	361	28	919	31	2132	111	4278

A partir de estos valores podemos calcular la odds ratio de la prevalencia de rodales tratados frente a no tratados para cada año entre las categorías de daño 3 o más y el resto. Para ello, haríamos el siguiente cálculo (como ejemplo tomaremos las categorías 3 o más y 2):

A = Rodales daño >= 3 Trat. / Rodales daño >= 3 Sin

B = Rodales daño = 2 Trat. / Rodales daño = 2 Sin

Odds ratio = A / B

Un valor positivo de la odds ratio indicaría que la presencia de tratamientos en rodales con un nivel daño de 3 o más ocurre con ese valor más frecuentemente que en rodales con un nivel de daño de 2. Por ejemplo, una odds ratio de 5 indicaría que los rodales con daño 3 o más se tratan 5 veces más frecuentemente que los rodales con daño 2. Hacemos estos cálculos comparando el nivel de daño de 3 o más con el resto. Los resultados se muestran a continuación.

	OR >=3 / 2	OR >=3 / 1	OR >=3 / 0
2002	3.12	3.88	9.26
2003	1.63	2.21	3.04
2004	1.55	3.11	3.03
2005	4.17	6.06	12.71

De todo esto se deduce lo siguiente. Primero, no todos los rodales que sufren un nivel de daño alto son tratados, como decíamos en nuestro trabajo. Sólo un porcentaje relativamente pequeño como vemos en la primera tabla. Esto es tranquilizador en parte, porque como hemos demostrado en nuestro trabajo, las fumigaciones producen exactamente los mismos resultados que el no hacer nada (la razón es que el propio insecto es controlado por la escasez de recursos alimenticios y el aumento de predadores y parasitoides). Segundo, es cierto que los rodales con un nivel de daño intenso son tratados con mayor frecuencia que los rodales que tienen un nivel de daño medio o incluso bajo, pero resulta curioso ver qué todavía hay muchos rodales con niveles de daño bajo que son tratados ¿qué criterios -más allá de la proximidad a rodales con un grado de infestación alto- utilizan los técnicos para decidir que estos rodales merecen ser tratados? Por último, alguien podría pensar que esto podría invalidar los resultados de nuestro trabajo, pero hemos repetido los análisis publicados en el artículo de Forest Ecology and Management, utilizando sólamente los rodales con un nivel de daño 2. Los resultados, que no muestro aquí por brevedad, muestran claramente que tampoco hay diferencias en la respuesta cuando comparamos rodales fumigados y no fumigados con un nivel de daño moderado. Por tanto, las conclusiones de nuestro estudio siguen siendo válidas.

Un nuevo método para identificar "puntos calientes" de diversidad

Uno de los principales intereses de los biólogos de la conservación es la de priorizar áreas para la conservación. De manera muy sencilla, ésto suele hacerse de la siguiente forma:

1. Se seleccionan las áreas que se quieren priorizar. Éstas pueden tener una forma irregular (ej. parcelas forestales con tamaño y forma distinta) o regular (ej. celdas de 1 x 1, 5 x 5 o 10 x 10 km en las que un determinado territorio es dividido);
2. Se caracterizan dichas áreas de acuerdo a una o varias variables. Generalmente se utiliza la riqueza de especies (de uno o varios grupos taxonómicos), aunque también se puede utilizar la rareza, la vulnerabilidad o una combinación de todos estos índices (Rey Benayas & de la Montaña 2004);
3. Se ordenan estas áreas en orden descendente, siendo las más importantes la que más riqueza, rareza o vulnerabilidad tienen.
4. Se define un criterio para seleccionar sólo un pequeño número de estas áreas. Este criterio suele indicar la selección de un porcentaje de las áreas más importantes para la conservación de acuerdo a la(s) variable(s) seleccionada(s) sobre el total de áreas (ej. 1%, 5% o 10%). Obviamente, este criterio es totalmente arbitrario y depende de cada caso de estudio y de las preferencias del o de los investigadores que lo realizan para resaltar un número mayor o menor de áreas importantes para la conservación.
   

Cuando la variable utilizada para hacer la priorización es la riqueza de especies, los científicos suelen hablar de identificar "puntos calientes" de diversidad, aunque en realidad lo que están haciendo es priorizar unas zonas sobre otras en función del número de especies que tengan. Uno de los artículos más paradigmáticos en este contexto es el de Myers et al. (2000) en el que se identifican los "puntos calientes" de diversidad a nivel global, aunque esto mismo se puede hacer a escala regional o de paisaje (ej. Cayuela et al. 2006, Altamirano et al. 2010).

Una propuesta interesante

Recientemente, Bartolino et al. (2010) en un artículo publicado en Population Ecology, proponen un nuevo método para identificar "puntos calientes" de diversidad de forma objetiva. Su aproximación es bastante sencilla y aparentemente robusta. Los pasos a seguir son:

1. Se cuantifica la riqueza de especies (generalmente de uno o unos pocos táxones) en todas las unidades muestrales (localidades, sitios, celdas de 10 x 10 km, etc).
2. Se calcula la frecuencia de unidades muestrales que caen en cada valor de riqueza de especies. Por ejemplo, hay 3 localidades con 15 especies, 6 con 16 especies, 9 con 17 especies, y así sucesivamente hasta el máximo, que pongamos que es 1 localidad que tiene 82 especies.
3. Se dibuja un gráfico de la frecuencia acumulada de especies (eje y) para cada valor de riqueza (eje x). Estos valores se van a relativizar. Cada valor de frecuencia se divide por el total de unidades muestrales y cada valor de riqueza se divide por la máxima riqueza de especies encontrada en toda el área de estudio (82 en nuestro ejemplo anterior). El gráfico resultante se muestra en la figura 1a) y d) para dos ejemplos propios con datos de riqueza de aves (a) y plantas herbáceas (d) en áreas naturales protegidas de Centroamérica y en celdas de 10 x 10 km en el Reino Unido, respectivamente.
4. Se calcula la pendiente en cada punto. Ésta se muestra en la gráfica 1c) y f) para los dos casos de estudio anteriores. Según Bartolino y colaboradores las localidades que tienen una pendiente menor a 45º representan zonas dónde la riqueza de especies se encuentra en densidades altas mientras que las que tienen una pendiente mayor a 45º representan zonas dónde la riqueza de especies se encuentra en densidades bajas.
5. Estos autores sugieren que los "puntos calientes" de riqueza pueden ser aquellos en los que haya una pendiente menor a 45º a partir del punto de riqueza mayor que cumpla que la pendiente es = 45º.
   

Figura 1. Curvas de frecuencia acumulada de la riqueza relativa de especies (a, d), frecuencia (no acumulada) de la riqueza de especies (b, e) y la pendiente de la tangente a la curva de frecuencia acumulada (c, f) para dos casos de estudio. El primero hace referencia a la riqueza de aves en áreas naturales protegidas de Centroamérica (desde Guatemala hasta Panamá) excluyendo las islas (a, b, c) (F. Suzart et al., en revisión). El segundo hace referencia a la riqueza de plantas herbáceas en celdas de 10 x 10 k en Reino Unido (d, e, f) (F. Suzart et al., en preparación).

Pero no es oro todo lo que reluce...

El enfoque de Bartolino y colaboradores parece realmente prometedor y una forma objetiva para realmente seleccionar "puntos calientes" de riqueza o diversidad. Sin embargo, cuando uno mira en detalle el estudio se da cuenta de que no es oro todo lo que reluce y de que hay cosas que no acaban de encajar.

Empecemos por la primera. Cuando se calcula la pendiente a la recta tangente a la curva que se genera (la que se muestra en la figura 1a) y d) respectivamente) nos damos cuenta de que esta condición se alcanza en muchos localidades, tanto con riquezas de especies altas como con riquezas bajas (figura 1c) y f)). Por tanto, la afirmación de Bartolino y colaboradores de que estas áreas son zonas en dónde la riqueza de especies se encuentra en altas densidades es falsa. En realidad, una pendiente menor a 45º sólo indica que la riqueza de especies aumenta en al menos 1 por cada nueva localidad añadida al total de localidades muestreadas. Y es por eso que Bartolino y colaboradores, de forma muy inteligente, sugieren que el umbral de corte sea aquel en el que se cumpla que la pendiente es = 45º pero para aquella localidad que tenga la máxima riqueza. Pero esta regla es arbitraria. Podríamos sugerir en su lugar que se tomase el primer punto que cumpla que la pendiente es = 45º una vez superada la media o la mediana, que es la que marca la "tipicidad" de un determinado patrón, en este caso riqueza de especies. O la que cumpliera esta condición pasada el cuartil 75, etc.

En el caso de las aves, el punto de corte que se obtiene utilizando este método es de 475 especies de aves. Si un área protegida tiene 475 especies de aves o maś, es un "punto caliente" de diversidad de aves. Y sino no. Pero encontramos localidades con una riqueza sólo ligeramente inferior (460, 465, 467) que además cumplen que su pendiente de la tangente a la curva de acumulación es menor a 45º y que, curiosamente, no serían clasificados como "puntos calientes" según el método aquí descrito... ¿por qué? La respuesta es: "por arbitrariedad".

En definitiva, el método de Bartolino et al. (2010), aunque elegante, sencillo y fácil de implementar, es igual de arbitrario que los métodos utilizados hasta el momento, con la diferencia de que nos exime de la responsabilidad y del ejercicio intelectual que supone tener que reflexionar sobre el contexto, las necesidades y las posibilidades de conservación de cada caso de estudio concreto. No puede haber métodos objetivos para la selección de "puntos calientes" de diversidad, pero esto no supone, en mi opinión, ningún impedimento para seguir investigando qué áreas tienen más o menos especies de cara a promover una conservación efectiva.

Este estudio forma parte del siguiente trabajo: Cayuela, L., de Albuquerque, F.S. & Gálvez, L. Are mathematical approaches to hotspot identification biologically meaningful? Reply to Bartolino et al. (2010). In prep.

Better species distribution modeling needed for the tropics

In order to conserve the world's biodiversity we need to know where species are found. We also need to predict where they might be found if the climate changes or human activity alters habitats. One way of gaining such knowledge is through field studies. Such work on the ground produces lists of species and adds to museum collections. However many tropical areas have not yet been visited by scientists. Even the most detailed studies from the best known areas of the tropics are far from exhaustive. This means that accurate distribution maps are not available for many tropical species. In order to address the problem increasingly sophisticated computer models have been designed that aim to predict where species might occur based on current knowledge. These models can often add a great deal of value to the limited information available. However, models are only as good as the data from which they are built.

A new study, published in the journal Tropical Conservation Science, looks closely at recent attempts to apply species distribution modeling to guide conservation in the tropics. A review of the literature suggested that models built for the most threatened species are still likely to be the least reliable. Cayuela and colleagues found that output from species distribution models is rarely being used when setting conservation priorities. The difficulty is due to a chronic shortfall in the quantity and quality of data used to build models. Although there have been many improvements in the algorithms used for modeling, these advances can not be expected to address underlying weakness of the data. The work points out the need to continue to work on improved frameworks for sharing scarce yet invaluable data on tropical biodiversity. The authors also suggest that a more systematic approach to future data generation is needed in order to fill key gaps in the knowledge base used for tropical conservation.

L. Cayuela, D.J. Golicher, A.C. Newton, M. Kolb, F.S. de Alburquerque, E.J.M.M. Arets, J.R.M. Alkemade, & A.M. Pérez (2009). Better species distribution modeling needed for the tropics Tropical Conservation Science, 2 (3), 319-352.

¿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.

Estima de los atributos funcionales de ecosistemas forestales a partir del NDVI

En 2006, Domingo Alcaraz publicaba un trabajo en Global Ecology and Biogeography sobre la caracterización de los tipos funcionales de ecosistemas de la Península Ibérica a partir del uso de índices de vegetación obtenidos de imágenes de satélite. Este trabajo consiste, básicamente, en definir unos atributos que permitan distinguir desde un punto funcional distintos tipos de ecosistemas. En su trabajo, Alcaraz et al. (2006) definen con estos fines tres atributos: la integral del índice de verdor NDVI (NDVI-I), la diferencia entre el máximo y el mínimo NDVI (RREL) y el mes de máximo NDVI (MMAX) (Figura 1). El primero se relaciona funcionalmente con la productividad de los ecosistemas, mientras que el segundo y el tercero están más vinculados a la estacionalidad. Utilizando estos tres atributos, calculados para una serie temporal larga obtenida a partir de imágenes NOAA/AVHRR a una resolución de 1 x 1 km, Alcaraz et al. (2006) clasificaron todos los tipos de ecosistemas en la Península Ibérica a lo largo de gradientes de productividad y estacionalidad.

Figura 1. Los tres atributos del NDVI empleados en la caracterización de los tipos funcionales de ecosistemas por Alcaraz et al. (2006): la integral del NDVI (NDVI.I), la diferencia entre el máximo y el mínimo NDVI (RREL) y el mes de máximo NDVI (MMAX). Ejemplo tomado de un pinar de Pinus nigra en la Sierra de los Filabres, Almería.

Partiendo de esta idea, pensamos en utilizar estos atributos, no para distinguir distintos tipos funcionales de ecosistemas, sino para intentar detectar el decaimiento en plantaciones de Pinus sylvestris y P. nigra en la Sierra de los Filabres. Para ello seguimos los siguientes pasos:

1. Se promediaron en cada fecha los valores de las series temporales de NDVI obtenidos a partir de las imágenes MODIS (250 x 250 m). Esto supone unos 8 a 9 datos por fecha (2000-2008). Cómo estamos trabajando con valores promedios no hace falta interpolar los datos faltantes como hicimos anteriormente.
2. Se calculó la integral total anual (NDVI-I) como la suma del área bajo la curva en cada intervalo.
3. Se calculó la diferencia entre el mínimo y el máximo NDVI (RREL) y el mes del máximo NDVI (MMAX). Estos cálculos son bastante intuitivos así que no hace falta explicarlos (ver Figura 1).
   

Para los datos de Filabres, usamos el siguiente código en R, que permite calcular estos atributos a partir de los datos de las series temporales de NDVI obtenidos previamente de las imágenes MODIS, y obtener gráficos de dichos atributos para cada una de las 76 parcelas de estudio. Como en ocasiones anteriores, el código permite acceder directamente a los datos y ejecutar las distintas funciones diseñadas para tal fin.

load(url("http://archivos-para-subir.googlegroups.com/web/R+step+1.Rob?gda=GQGm0T4AAACkOfwqHAVd4YqgfIB09GDRsFEbW00qADF89i9HmLDg8UMvt5QpYAf3GGSwB4Eu5X3jsKXVs-X7bdXZc5buSfmx"))
ndvi <- points_ndvi@data[,-c(1:4)]
julday <- as.numeric(substr(labels(ndvi)[[2]], 5, 7)) year <- as.numeric(substr(labels(ndvi)[[2]], 1, 4))

NDVI.I <- rep(0,dim(ndvi)[1])
RREL <- rep(0, dim(ndvi)[1])
MMAX <- rep(0, dim(ndvi)[1])
for (i in 1:dim(ndvi)[1]) {
ndvi.mean <- tapply(as.numeric(ndvi[i, ]), as.factor(julday), mean)/10000
NDVI.i <- abs(ndvi.mean[1] - ndvi.mean[length(ndvi.mean)])/2 + min(c(ndvi.mean[1], ndvi.mean[length(ndvi.mean)]))
for (j in 1:(length(ndvi.mean)-1)) {
NDVI.i <- NDVI.i + abs(ndvi.mean[j+1] - ndvi.mean[j])/2 + min(c(ndvi.mean[j+1], ndvi.mean[j]))
}
NDVI.I[i] <- NDVI.i
RREL[i] <- max(ndvi.mean) - min(ndvi.mean)
MMAX[i] <- dimnames(ndvi.mean)[[1]][ndvi.mean == max(ndvi.mean)]
png(as.character(paste("NDVI.I", i, ".png", sep = ""), width = 550))
dimnames(ndvi.mean)[[1]] <- c("Jan", "", "Feb", "", "Mar", "", "Apr", "", "May", "", "Jun", "", "Jul", "", "Aug", "", "Sep", "", "Oct", "", "Nov", "", "Dec")
plot(c(0, ndvi.mean), type = "n", xlab = "", ylab = "NDVI", axes = F, ylim = c(0.1,0.9))
axis(1, at = c(2:24), labels = dimnames(ndvi.mean)[[1]], las = 2)
axis(2)
polygon(x = c(2:24, 24, 2), y = c(ndvi.mean, 0.1, 0.1), col = "grey90")
segments(x0 = 1, y0 = min(ndvi.mean), x1 = 1, y1 = max(ndvi.mean), col = "red", lwd = 3)
pos <- pmatch(ndvi.mean[ndvi.mean == max(ndvi.mean)], ndvi.mean)
segments(x0 = pos +1, y0 = 0.1, x1= pos+1, y1 = as.numeric(ndvi.mean[ndvi.mean == max(ndvi.mean)]), col = "blue", lwd = 3, lty = 3)
box()
title(main = paste(points_ndvi@data$Species[i], points_ndvi@data$Id[i]), sub =paste("Level of damage =", points_ndvi@data$Damage[i]))
text(x = 12, y = mean(c(max(ndvi.mean), 0.1)), labels = "NDVI.I", cex = 1.5)
text(x = 1.5, y =max(ndvi.mean) + 0.05, labels = "RREL", cex = 1.2)
text(x = ifelse(pos < y ="ifelse(pos" labels = "MMAX" cex =" 1.2)
dev.off()
}

El resultado de aplicar este código son tres vectores con los valores de estos atributos para cada una de las 76 parcelas (NDVI.I, RREL, MMAX) y 76 gráficas con la representación de dichos valores.

Estos valores se utilizarán junto a otras variables (p.e. insolación, altitud, etc) para intentar discriminar masas forestales con distintos grados de decaimiento. Pronto esperamos poder mostrar avances en esta línea, en la que llevamos trabajando ya varios meses.

Referencias:

Alcaraz, D., Paruelo, J.M. & Cabello, J. 2006. Identification of current ecosystem functional types in the Iberian Peninsula. Global Ecology and Biogeography 15: 200-212.

Presentación al 5º Congreso Forestal Español

Algunos de los resultados obtenidos hasta el momento en relación al uso de imágenes MODIS (250 x 250 m) para detectar el decaimiento en repoblaciones de pino silvestre y pino salgareño en la Sierra de los Filabres (Almería) van a ser presentados al 5º Congreso Forestal Español, que tendrá lugar en Ávila, del 21 al 25 de septiembre de 2009. A continuación se muestra el resumen y el texto completo en pdf.

Detección del decaimiento en masas de Pinus sylvestris L. y P. nigra Arnold. procedentes de repoblación en la Sierra de los Filabres (Almería) mediante el uso del sensor MODIS

Durante las últimas décadas se han detectado en Europa y Norte América numerosos casos de decaimiento en masas forestales asociados, de forma general, con la contaminación, la presencia de plagas forestales y el cambio climático. En España, se han descrito procesos de decaimiento en abetales, encinares y pinares. Se prevé que el cambio climático resulte en una intensificación de estos procesos, por lo que la detección temprana de este fenómeno es un paso esencial para la gestión forestal sostenible. El diseño de modelos con base biológica para llevar a cabo dicha detección en Andalucía es uno de los objetivos prioritarios del proyecto GESBOME.

En este estudio se muestran resultados preliminares sobre la aplicación de imágenes MODIS (resolución espacial ~ 250 m) como herramienta para el seguimiento del estado vegetativo de masas forestales en la Sierra de los Filabres (Almería). En esta localidad, el problema del decaimiento de pinares de Pinus nigra Arnold. y Pinus sylvestris L. se viene observando desde el año 2001.

Para ello, se seleccionaron 76 puntos de control en campo y se asignaron a tres niveles de daño: sin afectar, moderado e intenso. En cada parcela se extrajeron los valores del índice de vegetación NDVI (Normalized Difference Vegetation Index) de las imágenes MODIS cada 16 días desde comienzos del año 2000 hasta agosto de 2008. Las series temporales de los tres niveles de daños se compararon entre sí utilizando un índice de referencia basado en la mediana del total de las series temporales para todas las parcelas. Las señales espectrales permitieron discriminar correctamente el nivel de daño extremos de los daños intermedios y moderado, especialmente en pinares de Pinus sylvestris, por lo que pueden utilizarse como indicadores para la detección del decaimiento. El uso de series temporales basadas en índices espectrales puede contribuir además a explorar procesos causales del decaimiento forestal mediante la comparación con variables edafoclimáticas y/o fisiográficas.

Análisis de los factores de incitación del decaimiento en masas forestales

En la estructura inicial planteada para analizar el proceso de decaimiento en la Sierra de los Filabres, Almería, se propuso la exploración de los factores de incitación mediante el análisis correlacional de distintas variables (topográficas, estructurales, espaciales) con el porcentaje de decaimiento estimado de visu en campo.

Uno de los problemas técnicos que se planteó previo a la realización de dichos análisis fue el cálculo de las superficies de insolación, problema que se ha resuelto mediante la aplicación de modelos específicos de insolacióimplementados en GRASS. En una entrada reciente en este blog se proveen los detalles técnicos para la preparación de la bases de datos points. Esta base contiene la siguiente información:

1. Porcentaje de decaimiento observado in situ en campo.
2. Especie (Pinus nigra y P. sylvestris).
3. Área basimétrica del rodal (var. estructural).
4. Densidad de pies del rodal (var. estructural).
5. Coordenadas geográficas x e y (var. espaciales).
6. Elevación (var. topográfica).
7. Pendiente (var. topográfica).
8. Orientación (var. topográfica).
9. Insolación en el solsticio de invierno (var. topográfica).
10. Insolación en el solsticio de verano (var. topográfica).
11. Índice de humedad del suelo (var. topográfica).
    

La información referente a variables topográficas (6-10) ha sido recalculada a tres resoluciones distintas (10x10 m, 30x30 m, 50x50 m). El modelo digital de elevaciones, del que se deriva el resto de información, está a una resolución original de 10x10 m. Sin embargo, dado el grano de las imágenes MODIS utilizadas (250x250 m) y el posible error de localización de los puntos puede ser conveniente re-escalar la información a escalas más groseras como 30x30 m o 50x50 m.

¿Está esta información proyectada a distintas resoluciones muy correlacionada entre sí? Si es así (p.e. r = 0.95), no importará mucho a qué resolución midamos las variables. De lo contrario, habrá que investigar si el cambio de resolución tiene un efecto sobre los resultados de los análisis de correlación con la variable respuesta (i.e. decaimiento).

Las tres bases de datos (points10, points30, points50) con la información referente a variables topográficas calculada a tres resoluciones distintas (10x10 m, 30x30 m y 50x50 m respectivamente) se puede descargar en R con el siguiente código (sólo la base de datos de 30x30 contiene además la información referente a la variable índice de humedad del suelo).

load(url("http://archivos-para-subir.googlegroups.com/web/points_correlaciones.Rob?hl=es&gda=mttqMzwAAACkOfwqHAVd4YqgfIB09GDRafzhFXYXrkFb8NzylP1tVw7qabgqw0xDKbwB-h3MnSf9Wm-ajmzVoAFUlE7c_fAt"))
ls()

[1] "points10" "points30" "points50"

Representamos ahora la matriz de correlaciones.


Cómo podemos observar, el cambio de resolución en la elevación no tiene ningún efecto perceptible sobre los valores de esta variable. Sin embargo, en la pendiente y orientación sí que se observa un cambio importante en los valores de estas variables según se cambia la resolución, sobretodo cuando se cambia la resolución de 10x10 m o 30x30 m a 50x50 m.

Vamos a seleccionar una resolución de 30x30 m para las variables topográficas, ya que es más representativa de los procesos a escala de rodal que la de 10x10 m y, al mismo tiempo, no se pierde mucha información con respecto a una escala de más detalle, cosa que sí pasa con la resolución de 50x50 m. Representamos ahora gráficamente los datos.
Parece que las variables que influyen más como factores de incitación en el proceso de decaimiento son la elevación, sobretodo para el caso de Pinus sylvestris, y la insolación invernal. Relaciones con las variables topográficas reescaladas a otra resolución (10x10 m, 50x50 m) no cambian apenas los resultados obtenidos. Es necesario, no obstante, analizar en mayor profundidad estos datos.