Curso de Análisis de datos ecológicos en R - V Edición

Se abre ya el plazo de matriculación de la V Edición del Curso de Análisis de datos ecológicos en R, que desde hace años vengo impartiendo en el Centro Andaluz de Medio Ambiente de Granada. Este año será en octubre, del 14 al 18, con una duración de 32 horas.

El curso está dirigido a estudiantes, fundamentalmente de postgrado, profesores e investigadores con conocimientos básicos de estadística, que deseen aprender el manejo de R para el análisis de datos en ecología. El curso ofrece la oportunidad de aprender el manejo de una herramienta muy potente y versátil para el análisis de datos, por lo que contribuirá a la formación técnica de alumnos de postgrado con un perfil investigador, así como de investigadores postdoctorales y profesores.

Normalmente las plazas se cubren en el plazo de apenas unos días, así que recomiendo a quien esté interesado en asistir que contacte con la Fundación General Empresa - Universidad de Granada (cursos@fundacionugrempresa.es) para inscribirse.

Máster en caracterización y conservación de la diversidad biológica

Está abierto el plazo de preinscripción en el Máster Universitario en Técnicas de Caracterización y Conservación de la Diversidad Biológica, que se imparte en la Universidad Rey Juan Carlos.

Los objetivos generales de esta titulación son dos:

• Proporcionar una serie de metodologías y herramientas de uso frecuente en el desarrollo de tareas técnicas e investigadoras en el campo de la Ecología y Conservación de la Diversidad Biológica.

• Desarrollar las numerosas aplicaciones específicas que ofrecen las metodologías/ instrumentos, es decir, su utilidad práctica en el campo de la Ecología y Conservación (ej. las aplicaciones de modelos de distribución de especies en escenarios de calentamiento global, uso de modelos poblacionales y metapoblacionales para estudios de viabilidad de especies y desarrollo de planes de conservación).

El profesorado está constituido por personal de la URJC, además de profesores externos pertenecientes al GBIF y al CSIC.

Es un máster áltamente recomendable. Más información en: http://www.slideshare.net/txiriondo/presentacin-master-tccdv-2013-2014.

300

Hace tres años mi cuñado y yo nos compramos una cuba de acero inoxidable de 300 litros y empezamos a hacer vino. En estos tres años hemos aprendido mucho y ahora hemos aumentado la producción al comprar una cuba nueva de 1000 litros. La cosecha de este año ha sido bastante buena. Hemos producido dos vinos: un monovarietal de tempranillo y una mezcla de garnacha (80%) y tempranillo (20%). Hemos conseguido aromas a madera mediante la adición de chips de roble tostado, que añadimos al vino tras la fermentación alcohólica. Unos amigos que son diseñadores (Maria José y Amadeo, de La Clave Gráfica) nos han diseñado la etiqueta. El resultado lo podéis ver aquí. El vino se llama así: 300, en honor a nuestros comienzos en aquella cuba de 300 litros.

Y para terminar, una nota de cata elaborada por una amiga aficionada (más bien diría experta)... ¡gracias Carol!

Monovarietal de tempranillo: 

Vino joven, con tonos granates persistentes, lágrima en caída elegante y bien definida. Aromas primarios a madera tostada, secundarios a frutas rojas frescas (notas claras de frambuesa, grosella, fresa y algún toque de cereza). Perfil láctico proveniente de la fermentación bastante acentuado (toques de mantequilla y aromas a leche suaves). En boca redondo y suave al paladar. Terciarios no muy marcados y con 1-3 caudalías (relativamente normal en este tipo de vinos). Perfecto para acompañar ensaladas, quesos grasos (tipo camembert y similares), todo tipo de pescado, pollo, e incluso pato, conejo, codorniz/perdiz o similares.

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.

 

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.

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.

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.

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.

Si quieres apoyar esta iniciativa, firma tú también ¡Exige una casilla en tudeclaración de la renta para poder dedicar un 0,7% de tus impuestos a la Ciencia (I+D+i)!

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

Mayor riesgo de infección alimentaria por E. coli en ganado alimentado con grano

Aunque Escherichia coli es un organismo comensal que reside dentro del intestino del huésped, algunas cepas patógenas de E. coli pueden causar colitis hemorrágica en humanos. Una de las cepas más infecciosas de E. coli es la O157:H7. El ganado es un reservorio asintomático de esta cepa. Esto quiere decir que contiene la cepa, pero que ésta no produce síntomas de la enfermedad. Normalmente la cepa está presente en alrededor de un 30% del ganado y, en determinadas ocasiones, puede llegar a estar presente en hasta el 80% del del ganado.

 
Al ganado intensivo, ya sea en su aptitud productiva de carne o leche, se le alimenta frecuentemente con grano (cereales y leguminosas) con el objeto de aumentar la eficiencia alimentaria, ya que el grano contiene muchas más calorías que el pasto en forma de hidratos de carbono, fundamentalmente almidón. Cuando el ganado es alimentado con mucho grano, algo de almidón escapa a la degradación microbiana del rúmen y pasa al intestino, donde es fermentado. E. coli O157:H7 fermenta los azúcares que se liberan tras la ruptura del almidón en el colon. Esto implica que tiene alimento suficiente para reproducirse y aumentar sus poblaciones, por lo que aparece en mayores densidades en esta parte del aparato digestivo, lo que se traduce finalmente en una mayor presencia de esta cepa en las heces.

El ganado intensivo se hacina en pequeñas concentraciones parcelarias, dónde la acumulación deheces es alta. Como consecuencia de ello, se han producido intoxicaciones de E. coli O157:H7 en reses sacrificadas procedentes de dichas explotaciones. La contaminación se produce por el contacto exterior de los animales con las heces. Aunque en los mataderos hay sistemas de limpieza y desinfección, cuanto mayor sea la proporción de E. coli, mayor será el riesgo de fallo de estos sistemas y ésta acabará entrando en contacto con la carne procesada tras el sacrificio.

Otra forma de intoxicación se produce cuando se utilizan estas heces como abono, como ocurrió recientemente en la famosa crisis del pepino (aunque la cepa aquí era otra). En este caso la intoxicación se produce a través de alimentos de origen vegetal, aunque este tipo de intoxicaciones son más raras.

¿Qué se puede hacer? No es necesario hablar de sistemas alternativos a los sistemas de explotación ganaderos intensivos. Pero sin necesidad de movernos hacia posturas totalmente contrarias, una revisión no ya tan reciente (Callaway et al. 2003) demuestra que basta con que el ganado sea alimentado 5 días con pasto para que las poblaciones de E. coli en el colon disminuyan del orden de 1000 veces su densidad, y que la capacidad de las poblaciones de E. coli para sobrevivir a los ácidos gástricos del estómago humano descienda drásticamente.

Por supuesto que también es necesario aumentar los controles sanitarios en explotaciones intensivas para evitar este tipo de problemas.

Callaway TR, Elder RO, Keen JE, Anderson RC, & Nisbet DJ (2003). Forage feeding to reduce preharvest Escherichia coli populations in cattle, a review. Journal of dairy science, 86 (3), 852-60 PMID: 12703622

Niños delante del televisor: ¿qué hacemos los padres para potenciarlo?

Como padre en ejercicio me he sentido muy identificado con la lectura de un artículo de investigación recientemente publicado en la revista Public Health Nutrition titulado “Parental factors associated with screen time in pre-school children in primary-care practice: a TARGet Kids! Study”. Se trata de un estudio realizado en Toronto, Canadá, con una población muestral de niños de 3 años. Se entrevistaba a los padres y se recogían datos de la cantidad de tiempo que los niños pasaban delante del televisor diariamente, además de información relacionada con las condiciones socio-económicas, culturales y laborales de los padres. El estudio pone en evidencia una relación positiva entre tiempo de exposición a la televión y el hábito de comer viendo la tele, madres trabajadoras, y padres que pasan más tiempo delante de la tele. También detectó una relación negativa entre tiempo de exposición a la tele y la existencia de normas en la casa con respecto a cuándo se puede ver y no ver la tele.

En conclusión, como padres, además de tener que predicar con el ejemplo sobre el consumo moderado de programas televisivos, tenemos que poner algunas normas a nuestros hijos para que sepan que la tele no es un entretenimiento al que puedan recurrir a cada rato. En cualquiera de los casos, hay que saber que la televisión en niños menores de 3 años tiene siempre efectos negativos, a saber: retrasa el desarrollo del lenguaje, genera comportamientos agresivos, sobrepeso y obesidad (los niños delante de la tele no se dan cuenta de qué están comiendo y comen, por lo general, por encima de sus necesidades). La propia Sociedad Canadiense de Pediatría recomienda que los niños de pre-escolar (menores de 3 años) no vean la tele más de 1 hora diaría. No sé que dirá la Asociación Española de Pediatría. En un vistazo rápido no he encontrado nada en su página web, pero sí he encontrado algunos artículos más de investigaciones españolas que alertan del uso y abuso de la televisión.

Así que ¡tomemos nota y tengamos cuidado sobre lo que les damos a nuestros hijos!

Slow food: un nuevo paradigma en la forma de alimentarse

Ayer estuve en una charla del fundador del movimiento Slow Food, Carlo Petrini, que me dió mucho que pensar. Slow Food es un movimiento que nace en 1989 en París, en oposición al Fast Food o comida rápida. Reivindica un nuevo concepto de gastronomía, más multidisciplinar, que considera los alimentos no sólo como algo lúdico o social, sino desde la perspectiva de la física, la química, la biología, la agronomía (cómo se producen los alimentos), la genética, la historia, la antropología, la cultura, la economía, etc.

Carlo Petrini, fundador del movimiento Slow Food, considerado por la revista Time como una de las 50 personas más influyentes del mundo.

En su charla, Carlo recalcó varias veces que tenemos un sistema alimentario "criminal" (dicho con sus propias palabras), pero también nos recordó que nosotros somos responsables directos de ello. Planteó nuestro papel no como consumidores pasivos, sino como co-productores, en dónde nuestro consumo influye directamente en los sistemas de producción y condiciona el tipo de agricultura. También enfatizó el hecho de que nos encontramos en una crisis entrópica, como consecuencia principalmente de que nos cuesta mucho más producir los alimentos que la energía derivada de los propios alimentos. Por ejemplo, producir un 1 kg de carne cuesta aproximadamente 5,000 lts de agua, además de una gran inversión energética como consecuencia de la dependencia en combustibles fósiles que tienen los actuales sistemas de producción: el grano para alimentar el ganado viene normalmente de otras partes del mundo, los fertilizantes necesarios para producir este grano también vienen de otros países, y todo este movimiento se traduce en consumo de combustibles fósiles (ver también esta otra entrada).

Los sistemas de producción de alimentos no son sostenibles en la actualidad por diversos motivos:

1. Agotan la fertilidad del suelo.
2. Agotan los recursos hídricos. Recordemos que el 72% del agua dulce en el mundo se destina a agricultura. Estos usos también condicionan la calidad del agua residual, que en muchas partes del mundo está contaminada.
3. Producen pérdida de biodiversidad. La selección de variedades, razas y especies más productivas en agricultura y ganadería ha motivado la desaparición del 70% de éstas durante el último siglo.
4. Promueven la desaparición del campesinado, lo que implica en muchas ocasiones una pérdida cultural de gran importancia.
   

Los defensores de los sistemas de producción de alimentos que tenemos en la actualidad arguyen que no es posible mantener alimentada a la población mundial bajo otros esquemas productivos no intensivos. Sin embargo, los datos que arroja la FAO a este respecto son muy claros:

• A nivel mundial se producen alimentos para alimentar a 12 mil millones de personas.
• Hay actualmente cerca de 7 mil millones de personas en nuestro planeta.
• 1 mil millones de personas sufren de desnutrición.
• 1.7 mil millones de personas sufren obesidad.
• El 40% de la producción alimentaria acaba en la basura.

Todo ello implica que tenemos un sistema de producción de alimentos altamente ineficiente, que genera excedentes que acaban en la basura, que consume más energía que la que proporcionan los propios alimentos que se producen, y en dónde siguen existiendo grandes desigualdades en la distribución y acceso a los alimentos.

¿Podemos hacer algo? Claro que sí. Como decía Carlo, somos corresponsables de lo que ocurre, y está en nuestra mano cambiarlo. Slow Food nace con esa filosofía y hoy son más de 100,000 personas las que se han unido a este movimiento en todo el mundo. Para más información, vísita la página del proyecto.

Software libre y evolución

¿Qué tiene que ver el software libre con la evolución? Pues mucho. Mientras que el software privativo va creciendo en función de las decisiones tomadas por un grupo reducido de personas (dirigido, eso sí, en parte por las demandas de los usuarios), el software libre va siendo construido por los propios usuarios en función de sus necesidades. Aquellas herramientas que son más útiles para otros usuarios van a permanecer y mejorarse a merced de los cambios efectuados por otros usuarios, mientras que aquello que no se use va a acabar desapareciendo... ¡toda una muestra del poder de la evolución!

En esta línea argumental se mueve el artículo de Tufto & Cavallini titulado "Should wildlife biologists use free software?", publicado en la revista Wildlife Biology en 2005. Todo un alegato al uso de software libre. Y eso, añadiré, que muchas de las herramientas que se mencionan en este artículo, incluyendo el sistema operativo Linux o los software OpenOffice, R, GRASS y Quantum GIS, han mejorado enormemente en los últimos cinco años.

Pues eso ¡viva el software libre!

Tufto, J., & Cavallini, P. (2005). Should wildlife biologists use free software? Wildlife Biology, 11 (1), 67-76 DOI: 10.2981/0909-6396(2005)11[67:SWBUFS]2.0.CO;2

Manual de comunicación para investigadores

Aquí va un vínculo muy interesante sobre como divulgar la ciencia, algo que , por desgracia, pocos científicos saben hacer bien (entre los que me incluyo). Aquí va un extracto:

"... la divulgación científica no sólo se hace desde los medios de comunicación. Se realiza también desde los museos o a través de libros y conferencias, y suele incluir la participación más o menos directa del científico. Centrándonos únicamente en la divulgación hecha desde los medios, podemos encontrar productos eminentemente divulgativos (revistas, suplementos de la prensa escrita, programas de radio y televisión) o noticias científicas mezcladas con un punto de divulgación, que se considera necesaria para entender algo mejor. En ocasiones la divulgación profundiza, otras pone en contexto, algo que es vital para hacer entender muchos temas. Muchas veces, ya que en ciencia se trabaja casi siempre desde la extrema especialización, el periodismo científico desempeña precisamente ese papel de proporcionar una visión de conjunto que la audiencia pueda digerir. En general, la divulgación también tiene mayor afán didáctico y suele ser más atemporal. Y no olvidemos que, si queremos hablar con la sociedad, deberemos ser capaces de entretener al lector, conseguir que se quede con nosotros hasta el final."

Se puede leer más en http://comunicaciencia.unirioja.es/.

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.