viernes, 24 de agosto de 2012

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.

 

lunes, 16 de julio de 2012

Material cursos de R para modelos mixtos y series temporales


Aunque con cierto retraso publico ahora el material del curso que impartí junto con Ana justel en mayo sobre series temporales en Granada. En concreto la parte aplicada de R, que es la que yo impartí. El material está disponible en la lista de enlaces que hay en el menú de la izquierda, bajo el epígrafe "Curso de análisis de datos ecológicos en R".

También hay disponible una actualización del curso sobre modelos mixtos, que impartí recientemente en la URJC, auspiciado por la AEET. Esta versión es bastante más completa e incluye una introducción a los GLMM (modelos lineales generalizados mixtos). También hay un vínculo nuevo a la wiki del curso (que se suma a la de dos cursos anteriores), en donde se pueden encontrar más ejemplos de GLM y GLMM resueltos por los alumnos.

miércoles, 15 de febrero de 2012

Curso de series temporales aplicadas a datos ambientales

En mayo, se impartirá la primera edición del curso titulado "Series temporales aplicadas a datos ambientales" en el Centro Andaluz de Medio Ambiente (CEAMA), en Granada. El curso lo impartiremos Ana Justel, profesora del Departamento de Matemáticas de la Universidad Autónoma de Madrid, y yo.



El curso ofrece la oportunidad de aprender los principales métodos y herramientas para el análisis estadístico de series temporales. Los objetivos principales de estos métodos son la predicción e interpretación de la evolución de fenómenos que se observan en intervalos regulares de tiempo. Mediante un enfoque práctico se pretende que a lo largo del curso los asistentes aprendan las principales técnicas para llevar a cabo un correcto análisis descriptivo de una serie temporal, trabajen con distintos métodos para el filtrado y extracción de tendencias y ciclos estacionales, y sean capaces de aplicar la metodología Box-Jenkins para la modelización ARIMA de series temporales.

Pronto publicaremos un manual con el material del curso, que se sumará a los ya publicados anteriormente sobre análisis de datos en ecología.

Esperamos que este curso sea de vuestro interés.

miércoles, 1 de febrero de 2012

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!

viernes, 27 de enero de 2012

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):


o *.zip (Windows):


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.

lunes, 16 de enero de 2012

Apoyo a la ciencia en la declaración de la Renta

El editorial de la revista científica Nature de Diciembre de 2011 llamaba a los nuevos gobiernos de España, Italia y Grecia a invertir más en Ciencia. Según los analistas, invertir en ciencia ahora hubiese traído beneficios desproporcionados. Impulsar la ciencia en el sur de Europa, decían, no sólo beneficiaría a estos países, si no que harían toda Europa más competitiva.

Dicho y hecho; continuando la tradición de los últimos años del gobierno Zapatero, el nuevo Gobierno de Rajoy en España ha anunciado un recorte de 600 millones de Euros en I+D+i. El presupuesto para 2011 (8600 millones) era ya un 8% menor que en 2010, a su vez 15% inferior que el de 2009. Queda claro entonces que este nuevo recorte deja a la Ciencia española en una situación de emergencia, en un contexto de crisis económica y política nacional. También hay que tener en cuenta al estudiar esas cifras que buena parte del I+D español financia investigación en tecnología militar, con lo que el dinero dedicado al tipo de investigación que querrían mayoritariamente los españoles es incluso menor. La investigación y la innovación son pilares fundamentales para el desarrollo de una sociedad moderna, y hemos demostrado muchas veces que los científicos españoles podemos ser tan buenos como los del resto de países si nos dan la oportunidad.  Por ello, el científico español Francisco J. Hernández ha impulsado una iniciativa para pedir a nuestros políticos que podamos elegir destinar el 0.7% de nuestros impuestos a la ciencia. Si esta iniciativa prospera, seríamos el primer país del mundo en el que los ciudadanos elegirían qué parte de sus impuestos se destinasen de forma explícita a I+D+i.

La revista Nature en su editorial de enero se hacía eco de esta iniciativa y clamaba que los científicos españoles habíamos conseguido en menos de una semana cerca de 32.000 firmas. Pues bien, el apoyo a esta iniciativa sigue creciendo y, tras poco más de dos semanas ya se cuenta con más de 138.000 firmas.

Buscar entradas